Package PyDSTool :: Module Events
[hide private]
[frames] | no frames]

Source Code for Module PyDSTool.Events

   1  """Event handling for python-based computations, and specification for
 
   2  both python and externally compiled code. (Externally compiled code
 
   3  may include its own event determination implementation.)
 
   4  
 
   5  "High-level" events are built in native Python function format.
 
   6  "Low-level" events are built for external platforms, e.g. C or Matlab code.
 
   7  
 
   8      Robert Clewley, October 2005.
 
   9  """ 
  10  
 
  11  # PyDSTool imports
 
  12  from Variable import * 
  13  from utils import * 
  14  from common import * 
  15  from parseUtils import * 
  16  from errors import * 
  17  from Interval import * 
  18  from utils import info as utils_info 
  19  from Symbolic import QuantSpec, Var, ensureStrArgDict 
  20  import FuncSpec 
  21  
 
  22  # Other imports
 
  23  import scipy, numpy, scipy, scipy.special 
  24  import math, random 
  25  import copy, types 
  26  
 
  27  __all__ = ['EventStruct', 'Event',
 
  28              'HighLevelEvent', 'LowLevelEvent', 'MatlabEvent',
 
  29              'makePythonStateZeroCrossEvent', 'makeZeroCrossEvent'] 
  30  
 
  31  
 
  32  # helper functions
 
33 -def _highlevel(arg):
34 return isinstance(arg[1], HighLevelEvent)
35
36 -def _lowlevel(arg):
37 return isinstance(arg[1], LowLevelEvent)
38
39 -def _term(arg):
40 return arg[1].termFlag
41
42 -def _nonterm(arg):
43 return not arg[1].termFlag
44
45 -def _active(arg):
46 return arg[1].activeFlag
47
48 -def _notactive(arg):
49 return not arg[1].activeFlag
50
51 -def _varlinked(arg):
52 return arg[1].varlinked
53
54 -def _notvarlinked(arg):
55 return not arg[1].varlinked
56
57 -def _precise(arg):
58 return arg[1].preciseFlag
59
60 -def _notprecise(arg):
61 return not arg[1].preciseFlag
62 63 64 # --------------------------------------------------------------------------- 65
66 -class EventStruct(object):
67 """A data structure to store and interface with multiple events.""" 68
69 - def __init__(self):
70 self.events = {} 71 self._keylist = ['term', 'nonterm', 'active', 'varlinked', 'precise', 72 'notprecise', 'notvarlinked', 'notactive', 'highlevel', 73 'lowlevel'] 74 self._makeFilterDict() 75 # record of recent events (mainly for checking event interval 76 # in high-level events) 77 self.resetEvtimes()
78
79 - def _makeFilterDict(self):
80 self._filterDict = {} 81 self._filterDict['highlevel'] = _highlevel 82 self._filterDict['lowlevel'] = _lowlevel 83 self._filterDict['term'] = _term 84 self._filterDict['nonterm'] = _nonterm 85 self._filterDict['active'] = _active 86 self._filterDict['notactive'] = _notactive 87 self._filterDict['varlinked'] = _varlinked 88 self._filterDict['notvarlinked'] = _notvarlinked 89 self._filterDict['precise'] = _precise 90 self._filterDict['notprecise'] = _notprecise
91
92 - def resetEvtimes(self):
93 self.Evtimes = {}
94
95 - def __deepcopy__(self, dummy):
96 # bug in deepcopy concerning _filterDict attribute, so work around 97 # - self.events is the only thing to really deep copy anyway 98 deepcopied_events = copy.deepcopy(self.events) 99 dcopy = copy.copy(self) 100 dcopy.events = deepcopied_events 101 return dcopy
102
103 - def __del__(self):
104 # not sure if this is necessary in order to evoke Event.__del__ 105 for ev in self.events: 106 del ev
107
108 - def __setitem__(self, ev):
109 if compareBaseClass(ev, Event): 110 if ev.name not in self.events: 111 self.events[ev.name] = ev 112 else: 113 print self 114 raise ValueError("Event name '"+ev.name+"' already present " 115 "in database") 116 elif isinstance(ev,list): 117 for ev_item in ev: 118 self.__setitem__(ev_item) 119 else: 120 raise TypeError('Argument must be an Event: received type ' 121 '%s'%str(type(ev)))
122 123 add = __setitem__ 124
125 - def __delitem__(self, ename):
126 del(self.events[ename])
127 128 delete = __delitem__ 129
130 - def __getitem__(self, ename):
131 return self.events[ename]
132
133 - def sortedEventNames(self, eventlist=None):
134 if eventlist is None: 135 return sortedDictKeys(self.events) 136 else: 137 eventlist.sort() 138 return eventlist
139
140 - def getHighLevelEvents(self):
141 hlList = [] 142 for epair in self.events.items(): 143 if isinstance(epair[1], HighLevelEvent): 144 hlList.append(epair) 145 return hlList
146
147 - def getLowLevelEvents(self):
148 llList = [] 149 for epair in self.events.items(): 150 if isinstance(epair[1], LowLevelEvent): 151 llList.append(epair) 152 return llList
153
154 - def getAllEvents(self):
155 return self.events.items()
156
157 - def getTermEvents(self):
158 teList = [] 159 for epair in self.events.items(): 160 if epair[1].termFlag: 161 teList.append(epair) 162 return teList
163
164 - def getNonTermEvents(self):
165 neList = [] 166 for epair in self.events.items(): 167 if not epair[1].termFlag: 168 neList.append(epair) 169 return neList
170
171 - def getActiveEvents(self):
172 neList = [] 173 for epair in self.events.items(): 174 if epair[1].activeFlag: 175 neList.append(epair) 176 return neList
177
178 - def getNonActiveEvents(self):
179 neList = [] 180 for epair in self.events.items(): 181 if not epair[1].activeFlag: 182 neList.append(epair) 183 return neList
184
185 - def getNonPreciseEvents(self):
186 neList = [] 187 for epair in self.events.items(): 188 if not epair[1].preciseFlag: 189 neList.append(epair) 190 return neList
191
192 - def getPreciseEvents(self):
193 neList = [] 194 for epair in self.events.items(): 195 if epair[1].preciseFlag: 196 neList.append(epair) 197 return neList
198
199 - def setglobalt0(self, t0):
200 for epair in self.events.items(): 201 epair[1].globalt0 = t0
202
203 - def query(self, keylist, eventlist=None):
204 """Return eventlist with results of queries corresponding 205 to self._keylist keys. 206 207 Multiple keys permitted in the query. 208 """ 209 if not isinstance(keylist, list): 210 raise TypeError('Query argument must be a list of keys') 211 for key in keylist: 212 if key not in self._keylist: 213 raise TypeError('Query keys must be in _keylist attribute') 214 filterFuncs = [self._filterDict[key] for key in keylist] 215 if eventlist is None: 216 if self.events == {}: 217 return [] 218 else: 219 eventlist = self.events.items() 220 if filterFuncs == []: 221 return [] 222 for f in filterFuncs: 223 eventlist = filter(f, eventlist) 224 return eventlist
225 226
227 - def __call__(self):
228 # info is defined in utils.py 229 utils_info(self.__dict__, "EventStruct")
230 231
232 - def info(self, verboselevel=0):
233 if verboselevel > 0: 234 # info is defined in utils.py 235 utils_info(self.__dict__, "EventStruct", 236 recurseDepthLimit=1+verboselevel) 237 else: 238 print self.__repr__()
239 240
241 - def __contains__(self, evname):
242 return evname in self.events
243 244 # If low level events are passed to this method, there will be an 245 # exception because they do not have __call__ method defined! 246 # This function is used for non-varlinked events, e.g. for use in 247 # integration loops. 248 # varDict must be supplied (not None) for this polling to work!
249 - def pollHighLevelEvents(self, tval=None, varDict=None, parDict=None, 250 eventlist=None):
251 if eventlist is None: 252 eventlist = self.query(['highlevel', 'active', 'notvarlinked']) 253 return filter(lambda (n,e): e(t=tval, varDict=varDict, 254 parDict=parDict), eventlist)
255 256
257 - def resetHighLevelEvents(self, t0, eventlist=None, state=None):
258 if eventlist is None: 259 eventlist = self.query(['highlevel']) 260 if type(state)==list: 261 for i, (ev, s) in enumerate(zip(eventlist, state)): 262 ev[1].reset(s) 263 ev[1].starttime = t0 264 else: 265 for ev in eventlist: 266 ev[1].reset(state) 267 ev[1].starttime = t0
268 269
270 - def validateEvents(self, database, eventlist):
271 """validateEvents is only used for high level events.""" 272 assert eventlist != [], 'Empty event list passed to validateEvents' 273 for ev in eventlist: 274 if isinstance(ev[1], HighLevelEvent): 275 if not reduce(bool.__and__, [key in database for key in \ 276 ev[1].vars.keys()]): 277 ek = ev[1].vars.keys() 278 print "Missing keys: ", remain(ek, database) 279 raise RuntimeError("Invalid keys in event '%s'" % ev[0])
280 #else: 281 # print "Warning: Low level events should not be passed to " \ 282 # + "validateEvents()" 283 # print " (event '", ev[0], "')" 284 285
286 - def setTermFlag(self, eventTarget, flagval):
287 if isinstance(flagval, bool): 288 if isinstance(eventTarget, list): 289 for evTarg in intersect(eventTarget, self.events.keys()): 290 self.events[evTarg].termFlag = flagval 291 else: 292 if eventTarget in self.events.keys(): 293 self.events[eventTarget].termFlag = flagval 294 else: 295 raise TypeError("Invalid flag type")
296 297
298 - def setActiveFlag(self, eventTarget, flagval):
299 if isinstance(flagval, bool): 300 if isinstance(eventTarget, list): 301 for evTarg in intersect(eventTarget, self.events.keys()): 302 self.events[evTarg].activeFlag = flagval 303 else: 304 if eventTarget in self.events.keys(): 305 self.events[eventTarget].activeFlag = flagval 306 else: 307 raise TypeError("Invalid flag type")
308
309 - def setPreciseFlag(self, eventTarget, flagval):
310 if isinstance(flagval, bool): 311 if isinstance(eventTarget, list): 312 for evTarg in intersect(eventTarget, self.events.keys()): 313 self.events[evTarg].preciseFlag = flagval 314 else: 315 if eventTarget in self.events.keys(): 316 self.events[eventTarget].preciseFlag = flagval 317 else: 318 raise TypeError("Invalid flag type")
319
320 - def setEventICs(self, eventTarget, val):
321 if isinstance(val, dict): 322 if isinstance(eventTarget, list): 323 for evTarg in intersect(eventTarget, self.events.keys()): 324 for name in val.keys(): 325 if name in self.events[evTarg].initialconditions.keys(): 326 self.events[evTarg].initialconditions[name] = val[name] 327 else: 328 if eventTarget in self.events.keys(): 329 for name in val.keys(): 330 if name in self.events[eventTarget].initialconditions.keys(): 331 self.events[eventTarget].initialconditions[name] = val[name] 332 else: 333 raise TypeError("Invalid ICs type -- must be dict of varname, value pairs")
334 335
336 - def setEventDelay(self, eventTarget, val):
337 if isinstance(val, int) or isinstance(val, float): 338 if isinstance(eventTarget, list): 339 for evTarg in intersect(eventTarget, self.events.keys()): 340 self.events[evTarg].eventdelay = val 341 else: 342 if eventTarget in self.events.keys(): 343 self.events[eventTarget].eventdelay = val 344 else: 345 raise TypeError("Invalid eventdelay type/value")
346
347 - def setEventInterval(self, eventTarget, val):
348 if isinstance(val, int) or isinstance(val, float): 349 if isinstance(eventTarget, list): 350 for evTarg in intersect(eventTarget, self.events.keys()): 351 self.events[evTarg].eventinterval = val 352 else: 353 if eventTarget in self.events.keys(): 354 self.events[eventTarget].eventinterval = val 355 else: 356 raise TypeError("Invalid eventinterval type/value")
357
358 - def setEventTol(self, eventTarget, val):
359 if val > 0: 360 if isinstance(eventTarget, list): 361 for evTarg in intersect(eventTarget, self.events.keys()): 362 self.events[evTarg].eventtol = val 363 else: 364 if eventTarget in self.events.keys(): 365 self.events[eventTarget].eventtol = val 366 else: 367 raise TypeError("Invalid eventtol type/value")
368
369 - def setEventDir(self, eventTarget, val):
370 if val in [-1,0,1] and isinstance(val, int): 371 if isinstance(eventTarget, list): 372 for evTarg in intersect(eventTarget, self.events.keys()): 373 self.events[evTarg].dircode = val 374 else: 375 if eventTarget in self.events.keys(): 376 self.events[eventTarget].dircode = val 377 else: 378 raise TypeError("Invalid eventdir type/value")
379
380 - def setStartTime(self, eventTarget, val):
381 if isinstance(val, int) or isinstance(val, float): 382 if isinstance(eventTarget, list): 383 for evTarg in intersect(eventTarget, self.events.keys()): 384 self.events[evTarg].starttime = val 385 else: 386 if eventTarget in self.events.keys(): 387 self.events[eventTarget].starttime = val 388 else: 389 raise TypeError("Invalid starttime type")
390 391
392 - def setBisect(self, eventTarget, val):
393 if isinstance(val, int) and val > 0: 394 if isinstance(eventTarget, list): 395 for evTarg in intersect(eventTarget, self.events.keys()): 396 self.events[evTarg].bisectlimit = val 397 else: 398 if eventTarget in self.events.keys(): 399 self.events[eventTarget].bisectlimit = val 400 else: 401 raise TypeError("Invalid bisectlimit type/value")
402 403 404
405 -class Event(object):
406 """Generic Event. 407 408 Possible keys in argument dictionary at initialization: 409 name, eventtol, eventdelay, starttime, bisectlimit, term, active, 410 precise, vars, expr. 411 """ 412
413 - def __init__(self, kw):
414 if 'name' in kw: 415 self.name = kw['name'] 416 else: 417 raise KeyError('Name must be supplied to event') 418 # absolute tolerance for event location (in dependent variable/expression) 419 if 'eventtol' in kw: 420 self.eventtol = kw['eventtol'] 421 else: 422 self.eventtol = 1e-9 423 # time interval before event detection begins on each run 424 if 'eventdelay' in kw: 425 self.eventdelay = kw['eventdelay'] 426 else: 427 self.eventdelay = 1e-3 428 # time interval between event detections restart 429 if 'eventinterval' in kw: 430 self.eventinterval = kw['eventinterval'] 431 else: 432 self.eventinterval = 1e-3 433 # number of bisection steps to take when finding events 434 if 'bisectlimit' in kw: 435 self.bisectlimit = kw['bisectlimit'] 436 else: 437 self.bisectlimit = 100 438 # terminating event flag 439 if 'term' in kw: 440 self.termFlag = kw['term'] 441 else: 442 self.termFlag = False 443 # active event flag 444 if 'active' in kw: 445 self.activeFlag = kw['active'] 446 else: 447 self.activeFlag = True 448 # determines whether event must be computed precisely 449 if 'precise' in kw: 450 self.preciseFlag = kw['precise'] 451 else: 452 self.preciseFlag = True 453 # store 'plain text' definition of event as string, if provided 454 if 'expr' in kw: 455 assert isinstance(kw['expr'],str), \ 456 "Invalid type for event definition string" 457 self._expr = kw['expr'] 458 else: 459 self._expr = None 460 # list of error structures 461 ## self.errors = [] 462 try: 463 self.dircode = kw['dircode'] 464 assert self.dircode in [-1, 0, 1] # direction codes 465 except AssertionError: 466 print 'invalid value for direction code -- must be -1, 0, or 1' 467 raise 468 except KeyError: 469 self.dircode = 0 470 # effective time zero for searches (used for eventdelay) 471 if 'starttime' in kw: 472 self.starttime = kw['starttime'] 473 else: 474 self.starttime = 0 475 # optional variable and parameter bounds information, in case an event 476 # wishes to refer to them 477 if 'xdomain' in kw: 478 assert type(kw['xdomain'])==dict, \ 479 "Invalid type for variable bounds information" 480 self.xdomain = kw['xdomain'] 481 else: 482 self.xdomain = {} 483 if 'pdomain' in kw: 484 assert type(kw['pdomain'])==dict, \ 485 "Invalid type for parameter bounds information" 486 self.pdomain = kw['pdomain'] 487 else: 488 self.pdomain = {} 489 # var dictionary (can be just dict of keys for individual 490 # value calls only) -- only for purely high level events 491 if 'vars' in kw: 492 self.vars = kw['vars'] 493 assert len(self.vars) > 0, 'vars dictionary must be non-empty' 494 if reduce(bool.__and__, [isinstance(var, Variable) for \ 495 var in self.vars.itervalues()]): 496 self.varlinked = True 497 # doesn't check that only argument is spec'd for _fn method 498 else: 499 self.varlinked = False 500 else: 501 raise KeyError('vars dictionary not present') 502 # _funcreg is a register of dynamically created Event method names 503 # in case of object copying (requiring destruction and re-creation 504 # of dynamically created methods) 505 self._funcreg = [] 506 self._funcstr = kw['funcspec'][0] 507 self._funcname = kw['funcspec'][1] 508 if 'auxfnspec' in kw: 509 self._fnspecs = ensureStrArgDict(kw['auxfnspec']) 510 else: 511 self._fnspecs = {} 512 if 'noHighLevel' in kw: 513 # Boolean to indicate whether the non-Python event type can make 514 # a high-level image of the event function 515 if not kw['noHighLevel']: 516 self.addMethods() 517 else: 518 # assume is high level 519 self.addMethods() 520 if 'prevsign_IC' in kw: 521 # set initial value -- useful for one-off tests or as part of 522 # map-based hybrid models with only one integer timestep 523 self.prevsign = kw['prevsign_IC'] 524 self.prevsign_IC = self.prevsign 525 else: 526 self.prevsign = None 527 self.prevsign_IC = None 528 # additional history for previous sign, in case of resetting 529 self.prevprevsign = None 530 self.fval = None 531 # placeholder for variables' initial conditions. If events use 532 # auxiliary functions that access 'initcond' auxiliary function 533 # then caller had better assign the current initial conditions 534 # to this dictionary first. 535 self.initialconditions = {} 536 # placeholder for global independent variable value, for use in 537 # hybrid systems 538 self.globalt0 = 0 539 # event queue information (used for discrete delays, recording 540 # features, etc.) 541 self.queues = {} 542 self._sorted_queues = [] 543 # event internal parameter information (for temp usage by 544 # event mappings) 545 self.evpars = {} 546 if 'evpars' in kw: 547 self.evpars.update(kw['evpars']) 548 # used for quadratic interpolation, if requested in searchForEvents 549 self.quadratic = None #fit_quadratic()
550 551 552 # Queues are intended to be an ordered sequence of numeric types, e.g. 553 # for times at which next terminal event should occur. 554 # Queues are first in, first out, and may be sorted or not. 555 # !! This feature is in development.
556 - def addToQ(self, qname, item):
557 try: 558 self.queues[qname].append(item) 559 except KeyError: 560 raise PyDSTool_ExistError("Queue %s was not declared"%qname) 561 else: 562 if qname in self._sorted_queues: 563 self.queues[qname].sort()
564
565 - def createQ(self, qname, sorted=True, seq=None):
566 """Also use to reset a queue.""" 567 if seq is None: 568 self.queues[qname] = [] 569 else: 570 # ensure list argument (for pop to work) 571 self.queues[qname] = list(seq) 572 if sorted: 573 self._sorted_queues.append(qname)
574
575 - def popFromQ(self, qname):
576 return self.queues[qname].pop(0)
577
578 - def deleteQ(self, qname):
579 if qname in self._sorted_queues: 580 i = self._sorted_queues.index(qname) 581 del self._sorted_queues[i] 582 try: 583 del self.queues[qname] 584 except KeyError: 585 pass
586
587 - def _infostr(self, verbose=1):
588 dirstr = ["decreasing", "either", "increasing"] 589 if verbose <= 0: 590 outputStr = "Event "+self.name 591 elif verbose > 0: 592 outputStr = "Event "+self.name+"\n" + \ 593 " active: " + str(self.activeFlag) + "\n" + \ 594 " terminal: " + str(self.termFlag) + "\n" + \ 595 " precise: " + str(self.preciseFlag) + "\n" + \ 596 " direction: " + dirstr[self.dircode+1] + "\n" + \ 597 " event tol: " + str(self.eventtol) + "\n" + \ 598 " definition: " + self._expr 599 if verbose >= 2: 600 outputStr += "\n bisect limit: " + str(self.bisectlimit) \ 601 + "\n" + \ 602 " event delay: " + str(self.eventdelay) + "\n" +\ 603 " event interval: " + str(self.eventinterval) 604 return outputStr
605 606
607 - def info(self, verboselevel=1):
608 print self._infostr(verboselevel)
609 610
611 - def __repr__(self):
612 return self._infostr(verbose=0)
613 614 __str__ = __repr__ 615 616
617 - def addMethods(self):
618 # clean up FuncSpec usage of parsinps and x for pars/inputs and 619 # variables 620 try: 621 exec self._funcstr 622 except: 623 print 'Invalid event function definition:' 624 print self._funcstr 625 raise 626 try: 627 setattr(self, '_fn', types.MethodType(locals()[self._funcname], 628 self, self.__class__)) 629 except KeyError: 630 print 'Must pass objective function for event at initialization' 631 raise 632 if '_fn' not in self._funcreg: 633 self._funcreg.append('_fn') 634 # _fnspecs keys are convenient, user-readable, short-hand names 635 # for the functions, whose actual function names are the second entry 636 # of the pair funcstr 637 for funcpair in self._fnspecs.values(): 638 # clean up FuncSpec usage of parsinps and x for pars/inputs and 639 # variables 640 try: 641 exec funcpair[0] 642 except: 643 print 'Invalid auxiliary function definition:' 644 print funcpair[0] 645 raise 646 try: 647 setattr(self, funcpair[1], types.MethodType(locals()[funcpair[1]], 648 self, 649 self.__class__)) 650 except KeyError: 651 print 'Must pass objective function for event at initialization' 652 raise 653 if funcpair[1] not in self._funcreg: 654 self._funcreg.append(funcpair[1])
655 656
657 - def reset(self, state=None):
658 """Reset event`s prevsign attribute to a certain state (defaults to None) 659 """ 660 self.fval = None 661 if state is None or self.dircode == 0: 662 if self.prevsign_IC is None and state is not None: 663 self.prevsign = self.prevprevsign 664 else: 665 self.prevsign = self.prevsign_IC 666 elif state == 'prev': 667 self.prevsign = self.prevprevsign 668 elif state == 'on': 669 self.prevsign = self.dircode 670 elif state == 'off': 671 self.prevsign = -self.dircode 672 else: 673 raise ValueError("Invalid state passed to event reset method")
674 675
676 - def __call__(self, t=None, varDict=None, parDict=None):
677 """Report on correct sign change. 678 For external inputs, add input names and vales at time t to parDict 679 """ 680 assert self.activeFlag, "Event cannot be called when inactivated" 681 if varDict is None: 682 if t is not None: 683 assert self.varlinked, ('wrong type of call for non var-' 684 ' linked event') 685 if t < self.starttime + self.eventdelay: 686 return False 687 else: 688 raise ValueError('t must be specified for this type of call') 689 if self.prevsign is None: 690 try: 691 self.fval = self._fn(t, parDict) 692 except: 693 print "Error in event", self.name 694 info(parDict, "Parameters") 695 raise 696 self.prevsign = scipy.sign(self.fval) 697 self.prevprevsign = None 698 return False 699 else: 700 try: 701 self.fval = self._fn(t, parDict) 702 except: 703 print "Error in event", self.name 704 info(parDict, "Parameters") 705 raise 706 sval = scipy.sign(self.fval) 707 if self.dircode == 0: 708 result = self.prevsign != sval 709 else: 710 result = self.prevsign != sval and \ 711 self.prevsign * self.dircode < 0 712 self.prevprevsign = self.prevsign 713 self.prevsign = sval 714 return result 715 else: 716 # then calling as individual value, not var object 717 varDict_temp = dict(varDict) 718 if t is None: 719 t = varDict['t'] 720 else: 721 assert 't' not in varDict, ('`t` key already present in ' 722 'varDict argument. Cannot redefine.') 723 varDict_temp['t'] = t 724 if t < self.starttime + self.eventdelay: 725 return False 726 assert not self.varlinked, ('wrong type of call for var-linked' 727 ' event') 728 if self.prevsign is None: 729 try: 730 self.fval = self._fn(varDict_temp, parDict) 731 except: 732 print "Error in event", self.name 733 info(parDict, "Parameters") 734 print "\n" 735 info(varDict_temp, "Variables") 736 raise 737 self.prevsign = scipy.sign(self.fval) 738 self.prevprevsign = None 739 return False 740 else: 741 try: 742 self.fval = self._fn(varDict_temp, parDict) 743 except: 744 print "Error in event", self.name 745 info(parDict, "Parameters") 746 print "\n" 747 info(varDict_temp, "Variables") 748 raise 749 sval = scipy.sign(self.fval) 750 if self.dircode == 0: 751 result = self.prevsign != sval 752 else: 753 result = self.prevsign != sval and \ 754 self.prevsign * self.dircode < 0 755 self.prevprevsign = self.prevsign 756 self.prevsign = sval 757 return result
758 759
760 - def searchForEvents(self, trange=None, dt=None, checklevel=2, 761 parDict=None, vars=None, inputs=None, 762 abseps=1e-13, eventdelay=True, globalt0=0):
763 """Search a variable-linked event, or an event with supplied vars 764 dictionary and relevant parameters, for zero crossings. 765 766 (Variable-linked search not applicable to low level events.) 767 768 trange=None, dt=None, checklevel=2, parDict=None, vars=None, inputs=None, 769 abseps=1e-13, eventdelay=True -> (ev_t, (ev_tlo, ev_thi)) 770 where the lo-hi tuple is the smallest bound around ev_t (in 771 case it is None because event was not found accurately). 772 773 dt will default to 1e-3 * the time interval of the variables. 774 'eventinterval' inherited from the event will be used to separate 775 detected events. 776 777 Only pass vars dictionary when event.varlinked is False. 778 """ 779 780 # flags indicating whether continuous and discrete variables are present 781 discretevars = False 782 continuousvars = False 783 if self.varlinked: 784 assert vars is None, ("event is variable-linked already -- " 785 "do not pass vars argument") 786 varDict = self.vars 787 else: 788 assert vars is not None, ("event is not variable-linked " 789 "already -- require a vars argument") 790 varDict = vars 791 precise = self.preciseFlag 792 try: 793 for var in varDict.itervalues(): 794 if isdiscrete(var): 795 discretevars = True 796 elif iscontinuous(var): 797 continuousvars = True 798 else: 799 raise TypeError('varDict must consist of Variable objects') 800 except AttributeError: 801 print 'event must contain a dictionary of vars' 802 raise 803 if discretevars and continuousvars: 804 raise TypeError('Cannot mix discrete and continuous Variable types') 805 if discretevars: 806 if precise: 807 print 'argument precise cannot be used for discrete Variable objects' 808 precise = False 809 varnames = varDict.keys() 810 if dt is not None: 811 print 'argument dt is unused for discrete Variable objects' 812 if trange is not None: 813 # trange had better be contained in all var.indepdomain ranges 814 assert len(varDict) > 0, 'varDict was empty!' 815 if not reduce(bool.__and__, [trange[0] in var.indepdomain \ 816 and trange[1] in var.indepdomain for var \ 817 in varDict.itervalues()]): 818 raise ValueError('trange not contained in all var ranges') 819 else: 820 raise ValueError('trange must be defined for discrete' 821 ' Variable objects') 822 tlist = varDict[varnames[0]].indepdomain 823 tlist = tlist[tlist.index(trange[0]):tlist.index(trange[1])] 824 elif continuousvars: 825 # trange had better be contained in all var.indepdomain ranges 826 tlimits = None # initial value 827 for var in varDict.itervalues(): 828 if tlimits is None: 829 tlimits = var.indepdomain.get() 830 else: 831 temp = var.indepdomain.get() 832 if temp[0] < tlimits[0]: 833 tlimits[0] = temp[0] 834 if temp[1] > tlimits[1]: 835 tlimits[1] = temp[1] 836 if trange is not None: 837 # compare trange and tlimits using Interval inclusion, 838 # so that can make use of uncertain values at boundaries 839 trange_int = Interval('trange', float, trange, abseps) 840 tlimits_int = Interval('tlimits', float, tlimits, abseps) 841 if not self.contains(tlimits_int, trange_int, checklevel): 842 raise ValueError('trange not contained in all var ranges') 843 else: 844 trange = tlimits 845 if dt is None: 846 dt = max(self.eventtol, 1e-3*(trange[1]-trange[0])) 847 if dt > trange[1]-trange[0]: 848 raise ValueError('dt (eventtol if not specified) is too large' 849 ' for trange in event %s'%self.name) 850 ttemp = copy.copy(var.indepdomain) 851 ttemp.set(trange) 852 tlist = ttemp.sample(dt, avoidendpoints=True) 853 else: 854 raise TypeError('var must be a Variable object') 855 if inputs is None: 856 test_fn = self._fn 857 self_caller = self.__call__ 858 else: 859 def test_fn(ts, ps): 860 if ps is None: 861 pidict = {} 862 else: 863 pidict = copy.copy(ps) 864 try: 865 for t in ts: 866 idict = dict([(n,i(t+globalt0)) for n,i in inputs.iteritems()]) 867 pidict.update(idict) 868 fvals.append(self._fn(t, pidict)) 869 except TypeError: 870 # Iteration over non sequence, so assume ts is just a single 871 # numeric value 872 try: 873 tvals = array(ts)+globalt0 874 except: 875 print "\n Found type", type(ts) 876 raise TypeError("t values invalid") 877 idict = dict([(n,i(tvals)) for n,i in inputs.iteritems()]) 878 pidict.update(idict) 879 fvals = self._fn(ts, pidict) 880 return fvals
881 def self_caller(t=None, varDict=None, parDict=None): 882 if parDict is None: 883 pidict = dict([(n,i(t+globalt0)) for n,i in inputs.iteritems()]) 884 else: 885 pidict = copy.copy(parDict) 886 pidict.update(dict([(n,i(t+globalt0)) for n,i in inputs.iteritems()])) 887 return self(t, varDict=varDict, parDict=pidict)
888 # before search, check that start after 'eventdelay', if option set 889 if eventdelay: 890 while tlist[0] <= self.eventdelay: 891 try: 892 tlist.pop(0) 893 except IndexError: 894 raise ValueError('eventdelay too large -- ' 895 'no search interval') 896 # now do the search 897 try: 898 if self.varlinked: 899 try: 900 fvals = test_fn(tlist, parDict) 901 boollist = [False] # initial point is always False 902 if self.dircode == 0: 903 boollist.extend([fvals[i] != fvals[i+1] for i in \ 904 range(len(fvals)-1)]) 905 else: 906 boollist.extend([fvals[i] * fvals[i+1] < 0 and \ 907 fvals[i] * self.dircode < 0 for i in \ 908 range(len(fvals)-1)]) 909 except: 910 # event does not support vectorized t calls 911 ## print "Warning: event did not support vectorized t calls" 912 # set start time for self-calling purposes 913 self.starttime = trange[0] 914 boollist = [self_caller(t, parDict=parDict) for t in tlist] 915 # first bool may be True because event hasn't been called 916 # before and has prevsign unset (or mis-set)... 917 # we can safely overwrite it (bit of a hack) because 918 # it's always really False 919 boollist[0] = False 920 else: 921 # vallist is a list of 1D lists 922 vallist = [v(tlist) for v in sortedDictValues(varDict)] 923 varnames = sortedDictKeys(varDict) 924 # set start time for self-calling purposes 925 self.starttime = trange[0] 926 if not eventdelay: 927 # switch off event delay temporarily 928 restore_val = self.eventdelay 929 self.eventdelay = 0. 930 boollist = [] 931 for tix in xrange(len(tlist)): 932 boollist.append(self_caller(t=tlist[tix], 933 varDict=dict(zip(varnames, 934 [vallist[i][tix] for\ 935 i in range(len(vallist))])), 936 parDict=parDict)) 937 # first bool may be True because event hasn't been called 938 # before and has prevsign unset (or mis-set)... 939 # we can safely overwrite it (bit of a hack) because 940 # it's always really False 941 boollist[0] = False 942 if not eventdelay: 943 # restore value for future use 944 self.eventdelay = restore_val 945 except KeyError, e: 946 # for discretevars if not all var.indepdomains equal over trange 947 if discretevars: 948 print 'Note: All discrete Variables must have identical ' \ 949 'independent variable domains over trange' 950 print "Check that all variable references in events are legitimate" 951 raise 952 # loop through boollist and find all events unless terminating at tpos! 953 tpos = -1 954 eventsfound = [] 955 t_last = -numpy.inf 956 t_interval = self.eventinterval 957 while True: 958 try: 959 tpos += boollist[tpos+1:].index(True)+1 960 except (ValueError, IndexError): 961 # no event found in this slice, at this dt 962 break 963 if precise and continuousvars: 964 if self.varlinked: 965 result = findpreciseroot(self, tlist[tpos-1], tlist[tpos], 966 parDict, inputs=inputs, 967 globalt0=globalt0, 968 quadratic_interp=self.quadratic) 969 else: 970 result = findpreciseroot(self, tlist[tpos-1], tlist[tpos], 971 parDict, varDict, inputs, 972 globalt0=globalt0, 973 quadratic_interp=self.quadratic) 974 if result is not None and result[0] > t_last + t_interval: 975 eventsfound.append(result) 976 t_last = result[0] 977 else: 978 result = tlist[tpos] 979 if result > t_last + t_interval: 980 eventsfound.append((result, (tlist[tpos-1], result))) 981 t_last = result 982 if self.termFlag: # quit searching now! 983 break 984 return eventsfound 985 986
987 - def contains(self, interval, val, checklevel=2):
988 # NB. val may be another interval 989 if checklevel == 0: 990 # level 0 -- no bounds checking at all 991 # code should avoid calling this function with checklevel = 0 992 # if possible, but this case is left here for completeness and 993 # consistency 994 return True 995 elif checklevel in [1,2]: 996 # level 1 -- ignore uncertain cases (treat as contained) 997 # level 2 -- warn on uncertain and continue (but warnings not 998 # used in Events, so treat as #1) 999 if interval.contains(val) is not notcontained: 1000 return True 1001 else: 1002 return False 1003 else: 1004 # level 3 -- exception will be raised for uncertain case 1005 if val in interval: 1006 return True 1007 else: 1008 return False
1009
1010 - def __getstate__(self):
1011 d = copy.copy(self.__dict__) 1012 for fname in self._funcreg: 1013 del d[fname] 1014 return d
1015
1016 - def __setstate__(self, state):
1017 self.__dict__.update(state) 1018 self.addMethods()
1019
1020 - def __copy__(self):
1021 pickledself = pickle.dumps(self) 1022 return pickle.loads(pickledself)
1023
1024 - def __deepcopy__(self, memo=None, _nil=[]):
1025 pickledself = pickle.dumps(self) 1026 return pickle.loads(pickledself)
1027 1028 1029
1030 -class HighLevelEvent(Event):
1031 """Event defined using python function code. 1032 1033 First argument to function must always be 'self' 1034 """
1035 - def __init__(self, kw):
1036 assert 'LLfuncspec' not in kw, \ 1037 "'LLfuncspec' key is invalid for purely high level events" 1038 Event.__init__(self, kw)
1039 1040
1041 -class LowLevelEvent(Event):
1042 """Event defined using externally-compiled and linked function code (i.e. C) 1043 1044 Specification of the function body as a string must include 1045 necessary temporary variable declarations and a return statement."""
1046 - def __init__(self, kw):
1047 if 'vars' in kw: 1048 for v in kw['vars'].itervalues(): 1049 if v is not None: 1050 raise TypeError("Low level events cannot be linked to a" 1051 " variable") 1052 Event.__init__(self, kw) 1053 assert isinstance(kw['LLfuncspec'], str), ("For low level events, must " 1054 "pass string for 'LLfuncspec' in initialization") 1055 LLfuncstr = kw['LLfuncspec'] 1056 dummyQ = QuantSpec('dummy', LLfuncstr, preserveSpace=True) 1057 dummyQ.mapNames({'abs': 'fabs', 'sign': 'signum', 'mod': 'fmod'}) 1058 if dummyQ[0] != 'return': 1059 print "Found: ", dummyQ[0] 1060 print "in event specification: ", dummyQ() 1061 raise ValueError("'return' must be first token in low level " 1062 "event specification") 1063 # take out 'return' from tokenized because whitespace will be lost 1064 # (parser not designed for non-math expressions) 1065 self._LLfuncstr = "return " + "".join(dummyQ[1:]) 1066 self._LLfuncname = self.name 1067 self._LLreturnstr = "double" 1068 self._LLargstr = "(unsigned n_, double t, double *Y_, double *p_, unsigned wkn_, double *wk_, unsigned xvn_, double *xv_)" 1069 self.varlinked = False # for compatibility with query calls
1070 1071
1072 -class MatlabEvent(LowLevelEvent):
1073 """Event defined using MATLAB syntax for use with ADMC++ 1074 1075 Specification of the function body as a string must include 1076 necessary temporary variable declarations."""
1077 - def __init__(self, kw):
1078 if 'vars' in kw: 1079 for v in kw['vars'].itervalues(): 1080 if v is not None: 1081 raise TypeError("Low level events cannot be linked to a variable") 1082 kw['noHighLevel'] = True 1083 Event.__init__(self, kw) 1084 assert isinstance(kw['Matlabfuncspec'], str), ("For low level events, must " 1085 "pass string for 'LLfuncspec' in initialization") 1086 LLfuncstr = kw['Matlabfuncspec'] 1087 dummyQ = QuantSpec('dummy', LLfuncstr, preserveSpace=True) 1088 if dummyQ[0] != 'return': 1089 print "Found: ", dummyQ[0] 1090 print "in event specification: ", dummyQ() 1091 raise ValueError("'return' must be first token in low level " 1092 "event specification") 1093 # take out 'return' from tokenized because whitespace will be lost 1094 # (parser not designed for non-math expressions) 1095 self._LLfuncstr = "Y_ = " + "".join(dummyQ[1:]) 1096 self._LLfuncname = self.name 1097 self._LLreturnstr = "Y_ = " 1098 self._LLargstr = "(t_, x_, p_)" 1099 self.varlinked = False # for compatibility with query calls
1100 1101 1102 # ------------------------------------------------------------- 1103 1104 ## Public exported functions
1105 -def makeZeroCrossEvent(expr, dircode, argDict, varnames=[], parnames=[], 1106 inputnames=[], fnspecs={}, targetlang='python', 1107 reuseterms={}, flatspec=None):
1108 """Target language-independent user-defined event involving coordinates, 1109 parameters, and time. Returns a non variable-linked event only. 1110 1111 List of used variable, parameter, and input names defaults to 1112 the empty list and can be omitted if these are not referenced. If 1113 variable names are omitted then the expression must only depend on time 1114 't' and any declared parameters. 1115 1116 Auxiliary function dictionary is required if the event accesses them. 1117 1118 'targetlang' argument defaults to 'python'. 1119 The expression 'expr' may not use intermediate temporary variables (for 1120 which you should specify the body of the event function by hand. 1121 1122 Optional argument reuseterms is a dictionary of terms used in expr 1123 that map to their definitions in terms of the state variables, time, 1124 parameters, and inputs. 1125 """ 1126 # preprocess expression for varnames -> v['varname'] and similarly 1127 # for parameter names 1128 try: 1129 exprname = argDict['name'] 1130 except KeyError: 1131 raise PyDSTool_KeyError("name key must be present in argDict") 1132 # support ModelSpec definitions, so take str() of expr 1133 expQS = QuantSpec('__ev_expr__', str(expr)) 1134 # make copies of arguments to prevent defaults getting messed up 1135 varnames = copy.copy(varnames) 1136 parnames = copy.copy(parnames) 1137 inputnames = copy.copy(inputnames) 1138 auxfns = copy.copy(fnspecs) 1139 auxVarDefMap = {} 1140 if flatspec is not None: 1141 # transfer used definitions in flatpsec into varnames, parnames, auxfns 1142 # etc. 1143 assert varnames==parnames==[], 'flatspec argument requires that varnames and parnames lists are not used' 1144 assert fnspecs=={}, 'flatspec argument requires that fnspecs is not used' 1145 try: 1146 allvars = flatspec['vars'] 1147 except KeyError: 1148 allvars = {} 1149 try: 1150 allinputs = flatspec['inputs'] 1151 except KeyError: 1152 allinputs = {} 1153 try: 1154 allpars = flatspec['pars'] 1155 except KeyError: 1156 allpars = {} 1157 1158 try: 1159 allfuns = flatspec['auxfns'] 1160 except KeyError: 1161 allfuns = {} 1162 collectReused(expQS, allvars, allpars, allinputs, allfuns, 1163 varnames, parnames, inputnames, auxfns, 1164 auxVarDefMap, flatspec) 1165 ## for symb in expQS.freeSymbols: 1166 ## if symb in allvars: 1167 ## spectype = flatspec['spectypes'][symb] 1168 ## if spectype == 'RHSfuncSpec': 1169 ## if symb not in varnames: 1170 ## varnames.append(symb) 1171 ## elif spectype == 'ExpFuncSpec': 1172 ## auxVarDef = allvars[symb] 1173 ## auxVarDefMap[symb] = auxVarDef 1174 ## # only need to look one level deeper because auxvars 1175 ## # cannot depend on each other, in a well-defined system. 1176 ## symbQS = QuantSpec('__symb_expr__', auxVarDef) 1177 ## for auxsymb in symbQS.freeSymbols: 1178 ## if auxsymb in allvars: 1179 ## if auxsymb not in varnames: 1180 ## varnames.append(auxsymb) 1181 ## spectype = flatspec['spectypes'][auxsymb] 1182 ## if spectype == 'ExpFuncSpec': 1183 ## auxVarDefMap[auxsymb] = allvars[auxsymb] 1184 ## elif auxsymb in allpars: 1185 ## if auxsymb not in parnames: 1186 ## parnames.append(auxsymb) 1187 ## elif auxsymb in allinputs: 1188 ## if auxsymb not in inputnames: 1189 ## inputnames.append(auxsymb) 1190 ## elif auxsymb in allfuns: 1191 ## if auxsymb not in auxfns: 1192 ## auxfns[auxsymb] = allfuns[auxsymb] 1193 ## elif symb in allpars: 1194 ## if symb not in parnames: 1195 ## parnames.append(symb) 1196 ## elif symb in allinputs: 1197 ## if symb not in inputnames: 1198 ## inputnames.append(symb) 1199 ## elif symb in allfuns: 1200 ## if symb not in auxfns: 1201 ## auxfns[symb] = allfuns[symb] 1202 ## # Add any dependencies of the function 1203 ## symbQS = QuantSpec('__symb_expr__', allfuns[symb][1]) 1204 ## for fnsymb in symbQS.freeSymbols: 1205 ## if fnsymb in allpars: 1206 ## if fnsymb not in parnames: 1207 ## parnames.append(fnsymb) 1208 ## if fnsymb in allinputs: 1209 ## if fnsymb not in inputnames: 1210 ## inputnames.append(fnsymb) 1211 ## elif fnsymb in allfuns: 1212 ## raise "Ask Rob to add support for function-to-function calls which he didn't anticipate here!" 1213 expr = '' 1214 done = False 1215 while not done: 1216 expQS = expQS.eval(auxVarDefMap) 1217 expr_new = str(expQS) 1218 if expr_new == expr: 1219 done = True 1220 else: 1221 expr = expr_new 1222 # support ModelSpec definitions, so take str() of list elements 1223 varnames = [str(v) for v in varnames] 1224 parnames = [str(p) for p in parnames] 1225 ## if inputnames != [] and targetlang != 'c': 1226 ## raise NotImplementedError("Inputs to non C-based events are not " 1227 ## "yet supported.") 1228 if exprname in varnames: 1229 raise ValueError("Expression name %s coincides with a variable name"%exprname) 1230 if exprname in parnames: 1231 raise ValueError("Expression name %s coincides with a parameter name"%exprname) 1232 if varnames == [] and inputnames == [] and 't' not in expQS: 1233 raise ValueError("In absence of variable or input names, time must appear in expression: %s"%expr) 1234 # zeros for t and varnames are unused placeholders. the names have 1235 # to be declared to funcspec with some value in order for the 1236 # object to be properly created. this is a kludge! 1237 specdict = {exprname: expr, 't': '0'} 1238 specdict.update({}.fromkeys(varnames, '0')) 1239 if 't' not in varnames: 1240 varnames.append('t') 1241 varnames.append(exprname) 1242 # Thus, if Jacobian etc. listed in funcspec then need to make a copy of the auxfns 1243 # and filter these out 1244 auxfns_temp = filteredDict(auxfns, ['Jacobian', 'Jacobian_pars', 'massMatrix'], 1245 neg=True) 1246 # Hack: targetLang must be python if (actually) c or python; matlab otherwise 1247 if targetlang == 'matlab': 1248 dummytlang = 'matlab' 1249 else: 1250 dummytlang = 'python' 1251 dummyfs = {'varspecs': specdict, 1252 'name': exprname, 'vars': varnames, 1253 'fnspecs': auxfns_temp, 1254 'pars': parnames, 'inputs': inputnames, 1255 'targetlang': dummytlang 1256 } 1257 dummyfuncspec = FuncSpec.RHSfuncSpec(dummyfs) 1258 parsedstr = dummyfuncspec._specStrParse([exprname], 1259 dummyfuncspec.varspecs, 1260 noreturndefs=True, 1261 forexternal=True) 1262 # alter par and var name dictionaries used by Events.py 1263 parsedstr = parsedstr.replace("parsinps[","p[").replace("x[","v[") 1264 if 'parsinps' in parsedstr: 1265 # then there are aux fns that use it 1266 funcstr = "parsinps=sortedDictValues(p,%s)+sortedDictValues(p,%s)"%(str(parnames),str(inputnames)) +\ 1267 "\n\treturn "+parsedstr+"\n" 1268 else: 1269 funcstr = "return "+parsedstr+"\n" 1270 if reuseterms != {}: 1271 illegalterms = ['globalindepvar', 'initcond', 'getindex', 'Jacobian', 1272 'Jacobian_pars'] 1273 reusestr, body_processed_dict = processReusedPy([exprname], 1274 {exprname: funcstr}, 1275 copy.copy(reuseterms), 1276 dummyfuncspec, 1277 illegal=illegalterms) 1278 funcstr_processed = (len(reusestr)>0)*"# local definitions\n" \ 1279 + reusestr + (len(reusestr)>0)*"\n" \ 1280 + body_processed_dict[exprname] 1281 else: 1282 funcstr_processed = funcstr 1283 # find used variable names (assumes matching braces in supplied spec!) 1284 # so use original parsedstr (before being processed for reused terms) 1285 varnames_found = [] 1286 currpos = 0 1287 done = False 1288 while not done: 1289 relfindpos = parsedstr[currpos:].find("v[") 1290 if relfindpos >= 0: 1291 findpos = relfindpos + currpos 1292 lbpos = findpos+1 1293 rbpos = findEndBrace(parsedstr[lbpos:], '[', ']')+lbpos 1294 varnames_found.append(parsedstr[lbpos+2:rbpos-1]) 1295 currpos = rbpos 1296 else: 1297 done = True 1298 inputnames_found = [] 1299 currpos = 0 1300 done = False 1301 while not done and inputnames != []: 1302 relfindpos = parsedstr[currpos:].find("p[") 1303 if relfindpos >= 0: 1304 findpos = relfindpos + currpos 1305 lbpos = findpos+1 1306 rbpos = findEndBrace(parsedstr[lbpos:], '[', ']')+lbpos 1307 in_name = parsedstr[lbpos+2:rbpos-1] 1308 if in_name in inputnames: 1309 inputnames_found.append(in_name) 1310 currpos = rbpos 1311 else: 1312 done = True 1313 newargs = ['dircode', 'funcspec', 'auxfnspec'] 1314 if targetlang in ['c', 'C']: 1315 newargs.append('LLfuncspec') 1316 # convert any special C-specific functions 1317 expr = dummyfuncspec._processSpecialC(expr) 1318 elif targetlang == 'matlab': 1319 newargs.append('Matlabfuncspec') 1320 # convert any special C-specific functions 1321 expr = dummyfuncspec._processSpecialC(expr) 1322 for arg in newargs: 1323 if arg in argDict: 1324 print 'Warning: `' + arg + '` already appears in argDict!' 1325 print ' This value will be overwritten.' 1326 argDict_out = copy.copy(argDict) 1327 funcname = "_f_"+exprname+"_ud" 1328 if parnames == []: 1329 pdefstr = "=None" 1330 else: 1331 pdefstr = "" 1332 # 'ds' plays role of 'self' (named for compatibility with FuncSpec's 1333 # automatic name resolution of auxiliary functions during parsing. 1334 if targetlang == 'matlab': 1335 funcstr_full = funcstr_processed 1336 else: 1337 funcstr_full = makeUniqueFn("def "+funcname+"(ds, v, p%s):\n\t"%pdefstr + \ 1338 funcstr_processed, idstr="event") 1339 if varnames_found+inputnames_found != []: 1340 argDict_out['vars'] = {}.fromkeys(varnames_found+inputnames_found, None) 1341 argDict_out.update({'funcspec': funcstr_full, 1342 'auxfnspec': dummyfuncspec.auxfns, 1343 'expr': expr, 1344 'dircode': dircode}) 1345 if targetlang in ['c', 'C']: 1346 LLfuncstr = "return "+expr+";\n" 1347 if reuseterms != {}: 1348 ## illegalterms = ['globalindepvar', 'initcond', 'getindex', 'Jacobian', 1349 ## 'Jacobian_pars'] 1350 LLreusestr, LLbody_processed_dict = processReusedC(['ev'], 1351 {'ev': LLfuncstr}, 1352 copy.copy(reuseterms)) #, 1353 ## illegal=illegalterms) 1354 LLfuncstr_processed = (len(LLreusestr)>0)*"/* local definitions */\n" \ 1355 + LLreusestr + (len(LLreusestr)>0)*"\n" + LLbody_processed_dict['ev'] 1356 else: 1357 LLfuncstr_processed = LLfuncstr 1358 argDict_out['LLfuncspec'] = LLfuncstr_processed 1359 return LowLevelEvent(argDict_out) 1360 elif targetlang == 'matlab': 1361 LLfuncstr = "return "+expr+";\n" 1362 if reuseterms != {}: 1363 LLreusestr, LLbody_processed_dict = processReusedMatlab(['ev'], 1364 {'ev': LLfuncstr}, 1365 copy.copy(reuseterms)) 1366 1367 LLfuncstr_processed = (len(LLreusestr)>0)*"% local definitions \n" \ 1368 + LLreusestr + (len(LLreusestr)>0)*"\n" + LLbody_processed_dict['ev'] 1369 else: 1370 LLfuncstr_processed = LLfuncstr 1371 argDict_out['Matlabfuncspec'] = LLfuncstr_processed 1372 return MatlabEvent(argDict_out) 1373 else: 1374 return HighLevelEvent(argDict_out)
1375 1376
1377 -def makePythonStateZeroCrossEvent(varname, targetvalue, dircode, argDict, 1378 var=None):
1379 """Python function-specified zero-crossing event in coordinate, or in 1380 time. Use 'var' argument to create a variable-linked event. 1381 1382 varname may be a Quantity object. dircode is -1, 0, or 1. 1383 varname may be the reserved word 't', for the independent variable. 1384 """ 1385 # set `var` only for non-var linked events 1386 newargs = ['dircode', 'funcspec', 'vars'] 1387 for arg in newargs: 1388 if arg in argDict: 1389 print 'Warning: `' + arg + '` already appears in argDict!' 1390 print ' This value will be overwritten.' 1391 if isinstance(varname, Var): 1392 # supporting a Variable object 1393 varname = varname.name 1394 elif not isinstance(varname, str): 1395 raise TypeError("Invalid type for event variable") 1396 funcname = "_f_"+varname+"_zc" 1397 if isinstance(targetvalue, str): 1398 # this is assumed to be a parameter name if not a numeric value, 1399 # but user preferred to add the p[' '] wrapping around a 1400 # parameter name 1401 if not (targetvalue[0:3] == "p['" and targetvalue[-2:] == "']"): 1402 targstr = "p['"+targetvalue+"']" 1403 else: 1404 targstr = targetvalue 1405 pdefstr = "=None" 1406 else: 1407 targstr = str(targetvalue) 1408 pdefstr = "" 1409 # 'ds' plays role of 'self' (named for compatibility with FuncSpec's 1410 # automatic name resolution of auxiliary functions during parsing. 1411 if isinstance(var, Variable): 1412 if not isinstance(targetvalue, str): 1413 assert targetvalue in var.depdomain, 'targetvalue not in var.depdomain' 1414 funcstr = "def "+funcname+"(ds, t, p%s):\n\treturn "%pdefstr + \ 1415 "ds.vars['"+varname+"'](t) - "+targstr+"\n" 1416 elif var is None: 1417 funcstr = "def "+funcname+"(ds, v, p%s):\n\treturn v['"%pdefstr + \ 1418 varname+"'] - " + targstr + "\n" 1419 else: 1420 raise ValueError('var must be a Variable object or None') 1421 argDict_out = copy.copy(argDict) 1422 argDict_out.update({'vars': {varname: var}, 1423 'funcspec': makeUniqueFn(funcstr, idstr="event"), 1424 'expr': varname + ' - ' + targstr, 1425 'dircode': dircode}) 1426 return HighLevelEvent(argDict_out)
1427 1428
1429 -def collectReused(quant, allvars, allpars, allinputs, allfuns, 1430 varnames, parnames, inputnames, auxfns, auxVarDefMap, flatspec):
1431 for symb in quant.freeSymbols: 1432 if symb in allvars: 1433 spectype = flatspec['spectypes'][symb] 1434 if spectype == 'RHSfuncSpec': 1435 if symb not in varnames: 1436 varnames.append(symb) 1437 elif spectype == 'ExpFuncSpec': 1438 auxVarDef = allvars[symb] 1439 auxVarDefMap[symb] = auxVarDef 1440 # need to recurse in case auxvar depends on another 1441 symbQS = QuantSpec('__symb_expr__', auxVarDef) 1442 collectReused(symbQS, allvars, allpars, allinputs, allfuns, 1443 varnames, parnames, inputnames, auxfns, 1444 auxVarDefMap, flatspec) 1445 elif symb in allpars: 1446 if symb not in parnames: 1447 parnames.append(symb) 1448 elif symb in allinputs: 1449 if symb not in inputnames: 1450 inputnames.append(symb) 1451 elif symb in allfuns: 1452 if symb not in auxfns: 1453 auxfns[symb] = allfuns[symb] 1454 # Add any dependencies of the function 1455 symbQS = QuantSpec('__symb_expr__', allfuns[symb][1]) 1456 collectReused(symbQS, allvars, allpars, allinputs, allfuns, 1457 varnames, parnames, inputnames, auxfns, 1458 auxVarDefMap, flatspec)
1459 1460
1461 -def processReusedPy(specnames, specdict, reuseterms, fspec, specials=[], 1462 dovars=True, dopars=True, doinps=True, illegal=[]):
1463 """Process reused subexpression terms for Python code. 1464 (Similar to function of similar name in FuncSpec.py) 1465 """ 1466 reused, specupdated, new_protected, order = FuncSpec._processReused(specnames, 1467 specdict, 1468 reuseterms, 1469 _indentstr) 1470 fspec._protected_reusenames = new_protected 1471 fspec.varspecs.update(specupdated) 1472 # symbols to parse are at indices 2 and 4 of 'reused' dictionary 1473 reusedParsed = fspec._parseReusedTermsPy(reused, [2,4], 1474 specials=specials, dovars=dovars, 1475 dopars=dopars, doinps=doinps, 1476 illegal=illegal) 1477 reusedefs = {}.fromkeys(new_protected) 1478 for vname, deflist in reusedParsed.iteritems(): 1479 for d in deflist: 1480 reusedefs[d[2]] = d 1481 return (concatStrDict(reusedefs, intersect(order,reusedefs.keys())), 1482 specupdated)
1483 1484
1485 -def processReusedC(specnames, specdict, reuseterms):
1486 """Process reused subexpression terms for C code. 1487 (Similar to function processReusedC in FuncSpec.py) 1488 """ 1489 reused, specupdated, new_protected, order = FuncSpec._processReused(specnames, 1490 specdict, 1491 reuseterms, 1492 '', 'double', ';') 1493 reusedefs = {}.fromkeys(new_protected) 1494 for vname, deflist in reused.iteritems(): 1495 for d in deflist: 1496 reusedefs[d[2]] = d 1497 return (concatStrDict(reusedefs, intersect(order, reusedefs.keys())), 1498 specupdated)
1499 1500
1501 -def processReusedMatlab(specnames, specdict, reuseterms):
1502 """Process reused subexpression terms for matlab code. 1503 (Similar to function of similar name in FuncSpec.py) 1504 """ 1505 reused, specupdated, new_protected, order = FuncSpec._processReused(specnames, 1506 specdict, 1507 reuseterms, 1508 '', '', ';') 1509 reusedefs = {}.fromkeys(new_protected) 1510 for vname, deflist in reused.iteritems(): 1511 for d in deflist: 1512 reusedefs[d[2]] = d 1513 return (concatStrDict(reusedefs, intersect(order, reusedefs.keys())), 1514 specupdated)
1515 1516
1517 -def findpreciseroot(ev, tlo, thi, parDict=None, vars=None, inputs=None, 1518 globalt0=0, quadratic_interp=None):
1519 """Find root more accurately from a Variable object using bisection. 1520 1521 (Adapted from scipy.optimize.minpack.bisection code to make use of 1522 quadratic interpolation, which assumes that tlo and thi are already known 1523 to be close enough together that the variable's curve is purely concave 1524 up or down in the neighbourhood, and so can be fitted accurately with a 1525 single quadratic). 1526 1527 To use quadratic interpolation, pass a fit_quadratic instance as the 1528 quadratic_interp argument. Interpolation will also be done on any inputs 1529 provided (**not yet implemented**). 1530 """ 1531 if tlo >= thi: 1532 raise ValueError('time limits are not finitely separated' 1533 ' or are given in wrong order') 1534 if inputs is None: 1535 if quadratic_interp is not None: 1536 raise NotImplementedError 1537 # 1538 q = quadratic_interp 1539 dt = thi-tlo #??? # may not be variable's mesh dt 1540 # and so far Event has not assumed that variable is either defined 1541 # using a mesh or not !!! 1542 ts = linspace(tlo, thi, 8) 1543 res = smooth_pts(ts, [ev._fn(t) for t in ts]) 1544 test_fn = res.results.f 1545 else: 1546 test_fn = ev._fn 1547 else: 1548 # only going to be calling this with individual t values 1549 if ev.varlinked: 1550 if quadratic_interp is not None: 1551 raise NotImplementedError 1552 else: 1553 def test_fn(x, p): 1554 if p is None: 1555 pidict = dict([(n,i(t+globalt0)) for n,i in inputs.iteritems()]) 1556 else: 1557 pidict = copy.copy(p) 1558 pidict.update(dict([(n,i(t+globalt0)) for n,i in inputs.iteritems()])) 1559 return ev._fn(t, pidict)
1560 else: 1561 if quadratic_interp is not None: 1562 raise NotImplementedError 1563 else: 1564 def test_fn(x, p): 1565 if p is None: 1566 pidict = dict([(n,i(x['t']+globalt0)) for n,i in inputs.iteritems()]) 1567 else: 1568 pidict = copy.copy(p) 1569 pidict.update(dict([(n,i(x['t']+globalt0)) for n,i in inputs.iteritems()])) 1570 return ev._fn(x, pidict) 1571 if ev.varlinked: 1572 assert vars is None, ("Only pass vars argument when event is not " 1573 "linked to a variable") 1574 assert reduce(bool.__and__, [iscontinuous(var) \ 1575 for var in ev.vars.itervalues()]), \ 1576 'Only pass continously-defined Variables in event' 1577 elo = test_fn(tlo, parDict) 1578 ehi = test_fn(thi, parDict) 1579 if elo == 0: 1580 return (tlo, (tlo, thi)) 1581 elif ehi == 0: 1582 return (thi, (tlo, thi)) 1583 elif elo * ehi > 0: 1584 # event cannot be present 1585 return None 1586 a = tlo 1587 b = thi 1588 i = 1 1589 eva = test_fn(a, parDict) 1590 rootival = (a,b) 1591 while i <= ev.bisectlimit: 1592 d = (b-a)/2.0 1593 p = a + d 1594 evp = test_fn(p, parDict) 1595 if abs(evp-eva) < ev.eventtol or evp == 0: 1596 return (p, rootival) 1597 i += 1 1598 if evp*eva > 0: 1599 a = p 1600 eva = evp 1601 else: 1602 b = p 1603 # do this at end of while loop in case i > bisectlimit 1604 rootival = (a,b) 1605 # search failed after bisectlimit 1606 return (None, rootival) 1607 else: 1608 assert vars is not None, ("vars argument is required when event is not" 1609 " linked to a variable") 1610 assert reduce(bool.__and__, [iscontinuous(var) \ 1611 for var in vars.itervalues()]), \ 1612 'Only pass continously-defined Variables' 1613 varnames = sortedDictKeys(vars)+['t'] 1614 dlo = dict(zip(varnames, [v(tlo) for v in sortedDictValues(vars)]+[tlo])) 1615 dhi = dict(zip(varnames, [v(thi) for v in sortedDictValues(vars)]+[thi])) 1616 elo = test_fn(dlo, parDict) 1617 ehi = test_fn(dhi, parDict) 1618 if elo == 0: 1619 return (tlo, (tlo, thi)) 1620 elif ehi == 0: 1621 return (thi, (tlo, thi)) 1622 elif elo * ehi > 0: 1623 # event cannot be present 1624 return None 1625 a = tlo 1626 b = thi 1627 i = 1 1628 da = dict(zip(varnames, [v(a) for v in sortedDictValues(vars)]+[a])) 1629 db = dict(zip(varnames, [v(b) for v in sortedDictValues(vars)]+[b])) 1630 eva = test_fn(da, parDict) 1631 rootival = (a,b) 1632 while i <= ev.bisectlimit: 1633 d = (b-a)/2.0 1634 p = a + d 1635 dp = dict(zip(varnames, [v(p) for v in sortedDictValues(vars)]+[p])) 1636 evp = test_fn(dp, parDict) 1637 if abs(evp-eva) < ev.eventtol or evp == 0: 1638 return (p, rootival) 1639 i += 1 1640 if evp*eva > 0: 1641 a = p 1642 eva = evp 1643 else: 1644 b = p 1645 # do this at end of while loop in case i > bisectlimit 1646 rootival = (a,b) 1647 # search failed after bisectlimit 1648 return (None, rootival) 1649