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
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
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
35
38
40 return arg[1].termFlag
41
43 return not arg[1].termFlag
44
46 return arg[1].activeFlag
47
49 return not arg[1].activeFlag
50
52 return arg[1].varlinked
53
55 return not arg[1].varlinked
56
58 return arg[1].preciseFlag
59
61 return not arg[1].preciseFlag
62
63
64
65
67 """A data structure to store and interface with multiple events."""
68
70 self.events = {}
71 self._keylist = ['term', 'nonterm', 'active', 'varlinked', 'precise',
72 'notprecise', 'notvarlinked', 'notactive', 'highlevel',
73 'lowlevel']
74 self._makeFilterDict()
75
76
77 self.resetEvtimes()
78
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
94
96
97
98 deepcopied_events = copy.deepcopy(self.events)
99 dcopy = copy.copy(self)
100 dcopy.events = deepcopied_events
101 return dcopy
102
104
105 for ev in self.events:
106 del ev
107
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
126 del(self.events[ename])
127
128 delete = __delitem__
129
131 return self.events[ename]
132
134 if eventlist is None:
135 return sortedDictKeys(self.events)
136 else:
137 eventlist.sort()
138 return eventlist
139
141 hlList = []
142 for epair in self.events.items():
143 if isinstance(epair[1], HighLevelEvent):
144 hlList.append(epair)
145 return hlList
146
148 llList = []
149 for epair in self.events.items():
150 if isinstance(epair[1], LowLevelEvent):
151 llList.append(epair)
152 return llList
153
155 return self.events.items()
156
158 teList = []
159 for epair in self.events.items():
160 if epair[1].termFlag:
161 teList.append(epair)
162 return teList
163
165 neList = []
166 for epair in self.events.items():
167 if not epair[1].termFlag:
168 neList.append(epair)
169 return neList
170
172 neList = []
173 for epair in self.events.items():
174 if epair[1].activeFlag:
175 neList.append(epair)
176 return neList
177
179 neList = []
180 for epair in self.events.items():
181 if not epair[1].activeFlag:
182 neList.append(epair)
183 return neList
184
186 neList = []
187 for epair in self.events.items():
188 if not epair[1].preciseFlag:
189 neList.append(epair)
190 return neList
191
193 neList = []
194 for epair in self.events.items():
195 if epair[1].preciseFlag:
196 neList.append(epair)
197 return neList
198
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
228
229 utils_info(self.__dict__, "EventStruct")
230
231
232 - def info(self, verboselevel=0):
233 if verboselevel > 0:
234
235 utils_info(self.__dict__, "EventStruct",
236 recurseDepthLimit=1+verboselevel)
237 else:
238 print self.__repr__()
239
240
242 return evname in self.events
243
244
245
246
247
248
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
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
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
281
282
283
284
285
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
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
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
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
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
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
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
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
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
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
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
414 if 'name' in kw:
415 self.name = kw['name']
416 else:
417 raise KeyError('Name must be supplied to event')
418
419 if 'eventtol' in kw:
420 self.eventtol = kw['eventtol']
421 else:
422 self.eventtol = 1e-9
423
424 if 'eventdelay' in kw:
425 self.eventdelay = kw['eventdelay']
426 else:
427 self.eventdelay = 1e-3
428
429 if 'eventinterval' in kw:
430 self.eventinterval = kw['eventinterval']
431 else:
432 self.eventinterval = 1e-3
433
434 if 'bisectlimit' in kw:
435 self.bisectlimit = kw['bisectlimit']
436 else:
437 self.bisectlimit = 100
438
439 if 'term' in kw:
440 self.termFlag = kw['term']
441 else:
442 self.termFlag = False
443
444 if 'active' in kw:
445 self.activeFlag = kw['active']
446 else:
447 self.activeFlag = True
448
449 if 'precise' in kw:
450 self.preciseFlag = kw['precise']
451 else:
452 self.preciseFlag = True
453
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
461
462 try:
463 self.dircode = kw['dircode']
464 assert self.dircode in [-1, 0, 1]
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
471 if 'starttime' in kw:
472 self.starttime = kw['starttime']
473 else:
474 self.starttime = 0
475
476
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
490
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
498 else:
499 self.varlinked = False
500 else:
501 raise KeyError('vars dictionary not present')
502
503
504
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
514
515 if not kw['noHighLevel']:
516 self.addMethods()
517 else:
518
519 self.addMethods()
520 if 'prevsign_IC' in kw:
521
522
523 self.prevsign = kw['prevsign_IC']
524 self.prevsign_IC = self.prevsign
525 else:
526 self.prevsign = None
527 self.prevsign_IC = None
528
529 self.prevprevsign = None
530 self.fval = None
531
532
533
534
535 self.initialconditions = {}
536
537
538 self.globalt0 = 0
539
540
541 self.queues = {}
542 self._sorted_queues = []
543
544
545 self.evpars = {}
546 if 'evpars' in kw:
547 self.evpars.update(kw['evpars'])
548
549 self.quadratic = None
550
551
552
553
554
555
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
571 self.queues[qname] = list(seq)
572 if sorted:
573 self._sorted_queues.append(qname)
574
576 return self.queues[qname].pop(0)
577
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
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):
609
610
613
614 __str__ = __repr__
615
616
618
619
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
635
636
637 for funcpair in self._fnspecs.values():
638
639
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
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
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
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
826 tlimits = None
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
838
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
871
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
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
897 try:
898 if self.varlinked:
899 try:
900 fvals = test_fn(tlist, parDict)
901 boollist = [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
911
912
913 self.starttime = trange[0]
914 boollist = [self_caller(t, parDict=parDict) for t in tlist]
915
916
917
918
919 boollist[0] = False
920 else:
921
922 vallist = [v(tlist) for v in sortedDictValues(varDict)]
923 varnames = sortedDictKeys(varDict)
924
925 self.starttime = trange[0]
926 if not eventdelay:
927
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
938
939
940
941 boollist[0] = False
942 if not eventdelay:
943
944 self.eventdelay = restore_val
945 except KeyError, e:
946
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
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
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:
983 break
984 return eventsfound
985
986
987 - def contains(self, interval, val, checklevel=2):
988
989 if checklevel == 0:
990
991
992
993
994 return True
995 elif checklevel in [1,2]:
996
997
998
999 if interval.contains(val) is not notcontained:
1000 return True
1001 else:
1002 return False
1003 else:
1004
1005 if val in interval:
1006 return True
1007 else:
1008 return False
1009
1011 d = copy.copy(self.__dict__)
1012 for fname in self._funcreg:
1013 del d[fname]
1014 return d
1015
1019
1021 pickledself = pickle.dumps(self)
1022 return pickle.loads(pickledself)
1023
1025 pickledself = pickle.dumps(self)
1026 return pickle.loads(pickledself)
1027
1028
1029
1031 """Event defined using python function code.
1032
1033 First argument to function must always be 'self'
1034 """
1036 assert 'LLfuncspec' not in kw, \
1037 "'LLfuncspec' key is invalid for purely high level events"
1038 Event.__init__(self, kw)
1039
1040
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."""
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
1064
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
1070
1071
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."""
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
1094
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
1100
1101
1102
1103
1104
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
1127
1128 try:
1129 exprname = argDict['name']
1130 except KeyError:
1131 raise PyDSTool_KeyError("name key must be present in argDict")
1132
1133 expQS = QuantSpec('__ev_expr__', str(expr))
1134
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
1142
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
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
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
1223 varnames = [str(v) for v in varnames]
1224 parnames = [str(p) for p in parnames]
1225
1226
1227
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
1235
1236
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
1243
1244 auxfns_temp = filteredDict(auxfns, ['Jacobian', 'Jacobian_pars', 'massMatrix'],
1245 neg=True)
1246
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
1263 parsedstr = parsedstr.replace("parsinps[","p[").replace("x[","v[")
1264 if 'parsinps' in parsedstr:
1265
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
1284
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
1317 expr = dummyfuncspec._processSpecialC(expr)
1318 elif targetlang == 'matlab':
1319 newargs.append('Matlabfuncspec')
1320
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
1333
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
1349
1350 LLreusestr, LLbody_processed_dict = processReusedC(['ev'],
1351 {'ev': LLfuncstr},
1352 copy.copy(reuseterms))
1353
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
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
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
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
1399
1400
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
1410
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
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
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
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
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
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
1540
1541
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
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
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
1604 rootival = (a,b)
1605
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
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
1646 rootival = (a,b)
1647
1648 return (None, rootival)
1649