1 """General purpose (hybrid) model class, and associated hybrid trajectory
2 and variable classes.
3
4 Robert Clewley, March 2005.
5
6 A Model object's hybrid trajectory can be treated as a curve, or as
7 a mapping. Call the model object with the trajectory name, time(s), and
8 set the 'asmap' argument to be True to use an integer time to select the
9 trajectory segment. These are numbered from zero.
10
11 A trajectory value in a Model object's 'trajectories' dictionary
12 attribute is a HybridTrajectory object, having the following
13 attributes (among others):
14
15 timeInterval is the entire time interval for the trajectory.
16
17 timePartitions is a sequence of time_intervals (for each trajectory
18 segment in trajseq), and
19
20 trajSeq is a list of epoch or regime trajectory segments [traj_0, traj_1,
21 ..., traj_(R-1)],
22
23 where traj_i is a callable Trajectory or HybridTrajectory object.
24
25 eventStruct is the event structure used to determine that trajectory.
26
27 events is a dictionary of event names -> list of times at which that
28 event took place.
29
30 modelNames is a list of the generators used for the trajectory (one per
31 partition).
32
33 variables is a dictionary that mimics the variables of the trajectory.
34 """
35
36
37
38
39 import Generator, Events, MProject
40 from utils import *
41 from common import *
42 from errors import *
43 from Interval import *
44 from Trajectory import *
45 from Variable import *
46 from Points import *
47 from ModelSpec import *
48 from Symbolic import isMultiRef
49 from parseUtils import isHierarchicalName, NAMESEP, mapNames, symbolMapClass
50
51
52 import math, sys
53 from numpy import Inf, NaN, isfinite, sign, abs, array, arange, \
54 zeros, concatenate, transpose, shape
55 from numpy import sometrue, alltrue, any, all
56 import copy
57 from time import clock
58 import pprint
59
60 __all__ = ['Model', 'HybridModel', 'NonHybridModel',
61 'boundary_containment', 'boundary_containment_by_postproc',
62 'boundary_containment_by_event',
63 'domain_test']
64
65
66
67
69
70
71 - def __init__(self, name, description='', pars=None):
72 MProject.qt_feature_leaf.__init__(self, name, description, pars)
73 try:
74 pars.thresh
75 except AttributeError:
76 raise ValueError("Missing threshold specification")
77 try:
78
79 pars.abseps
80 except AttributeError:
81 self.pars.abseps = 0
82 try:
83 assert pars.interior_dirn in [-1, 0, 1]
84 except AttributeError:
85 raise ValueError("Missing interior direction specification")
86 except AssertionError:
87 raise ValueError("Invalid interior direction specification value "
88 " use -1 for 'below', 1 for 'above', or 0 for 'discrete domain'")
89 try:
90 self.pars.coordname
91 except AttributeError:
92
93 self.pars.coordname = None
94 else:
95 assert isinstance(self.pars.coordname, str), \
96 "Coordinate name must be a string"
97
98 def evaluate(self, target):
99 raise NotImplementedError("Only call this method on a concrete "
100 "sub-class")
101
102
104 - def __init__(self, name, description='', pars=None):
105 boundary_containment.__init__(self, name, description, pars)
106 try:
107 self.pars.bd_eventname
108 except AttributeError:
109 raise ValueError("Missing boundary event name")
110
111
112 try:
113 self.pars.model
114 except AttributeError:
115 raise ValueError("Missing model associated with event")
116
117
119
120 try:
121 evpts = traj.getEvents(self.pars.bd_eventname)
122 except ValueError, errinfo:
123 print errinfo
124 raise RuntimeError("Could not find flagged events for this trajectory")
125 try:
126 evpt = evpts[self.pars.coordname]
127 except KeyError:
128 raise ValueError("No such coordinate %s in the defined event"%self.pars.coordname)
129 except TypeError:
130
131
132 if self.pars.abseps > 0:
133
134
135 print "Warning -- Boundary containment feature %s:"%self.name
136 print " Check for uncertain case using events not implemented"
137 self.results.output = None
138 self.results.uncertain = False
139 else:
140 self.results.output = None
141 self.results.uncertain = False
142 return True
143 else:
144
145 if self.pars.abseps > 0:
146
147
148 print "Warning -- Boundary containment feature %s:"%self.name
149 print " Check for uncertain case using events not implemented"
150 self.results.output = evpt[0]
151 self.results.uncertain = False
152 else:
153 self.results.output = evpt[0]
154 self.results.uncertain = False
155 return False
156
158 """Helper function for finding index in trajectory meshpoints
159 at which containment first failed."""
160 if self.results.satisfied:
161
162 return None
163 return len(self.results.output)
164
165
166 -class boundary_containment_by_postproc(boundary_containment):
167 - def evaluate(self, traj):
168 diffs = [p - self.pars.thresh for p in \
169 traj.sample(coords=self.pars.coordname)]
170 if self.pars.verbose_level > 1:
171 print "%s diffs in coord %s ="%(self.name,self.pars.coordname), diffs
172 res_strict = array([sign(d) \
173 == self.pars.interior_dirn for d in diffs])
174 satisfied_strict = alltrue(res_strict)
175 if self.pars.abseps > 0:
176 if self.pars.interior_dirn == 0:
177
178 res_loose = array([abs(d) < self.pars.abseps for d in diffs])
179 else:
180 res_loose = array([sign(d + self.pars.interior_dirn*self.pars.abseps) \
181 == self.pars.interior_dirn for d in diffs])
182 satisfied_loose = alltrue(res_loose)
183 self.results.output = res_loose
184
185
186
187 self.results.uncertain = satisfied_loose and not satisfied_strict
188 return satisfied_loose
189 else:
190 self.results.output = res_strict
191 self.results.uncertain = sometrue(array(diffs)==0)
192 return alltrue(res_strict)
193
194 - def _find_idx(self):
195 """Helper function for finding index in trajectory meshpoints
196 at which containment first failed"""
197 if self.results.satisfied:
198
199 return None
200 res = self.results.output
201 if res[0] == -1:
202 adjusted_res = list((res + 1) != 0)
203 elif res[0] == 1:
204 adjusted_res = list((res - 1) != 0)
205 else:
206
207 adjusted_res = list(res != 0)
208
209
210
211 return adjusted_res.index(True)
212
213
214 -class domain_test(MProject.qt_feature_node):
215 - def __init__(self, name, description='', pars=None):
216 MProject.qt_feature_node.__init__(self, name, description, pars)
217 try:
218 self.pars.interval
219 except AttributeError:
220 raise ValueError("Missing domain interval specification")
221 try:
222
223 self.pars.abseps
224 except AttributeError:
225 if isinstance(self.pars.interval, Interval):
226
227 self.pars.abseps = self.pars.interval._abseps
228 else:
229 self.pars.abseps = 0
230 if not isinstance(self.pars.interval, (tuple, list)):
231
232 self.pars.interval = [self.pars.interval, self.pars.interval]
233 elif len(self.pars.interval)==1:
234
235
236 self.pars.interval = [self.pars.interval[0],
237 self.pars.interval[0]]
238 self.isdiscrete = self.pars.interval[0] == self.pars.interval[1]
239 try:
240 self.pars.coordname
241 except AttributeError:
242 raise ValueError("Missing coordinate name")
243 try:
244 self.pars.derivname
245 except AttributeError:
246 raise ValueError("Missing coordinate derivative name")
247
248
249
250
251
252
253 xlo_bc = boundary_containment_by_postproc('x_test_lo',
254 description='Test x lower bound',
255 pars=args(thresh=self.pars.interval[0],
256 interior_dirn=1*int(not self.isdiscrete),
257 abseps=self.pars.abseps,
258 coordname=self.pars.coordname))
259 xhi_bc = boundary_containment_by_postproc('x_test_hi',
260 description='Test x upper bound',
261 pars=args(thresh=self.pars.interval[1],
262 interior_dirn=-1*int(not self.isdiscrete),
263 abseps=self.pars.abseps,
264 coordname=self.pars.coordname))
265 dxlo_bc = boundary_containment_by_postproc('dx_test_lo',
266 description='Test dx at lower x bound',
267 pars=args(thresh=0,
268 interior_dirn=1*int(not self.isdiscrete),
269 abseps=0,
270 coordname=self.pars.derivname))
271 dxhi_bc = boundary_containment_by_postproc('dx_test_hi',
272 description='Test dx at upper x bound',
273 pars=args(thresh=0,
274 interior_dirn=-1*int(not self.isdiscrete),
275 abseps=0,
276 coordname=self.pars.derivname))
277 self.subfeatures = {'x_test_lo': xlo_bc,
278 'dx_test_lo': dxlo_bc,
279 'x_test_hi': xhi_bc,
280 'dx_test_hi': dxhi_bc}
281
282 - def evaluate(self, traj):
283
284 for sf in self.subfeatures.values():
285 self.propagate_verbosity(sf)
286 sf.super_pars = self.pars
287 sf.super_results = self.results
288 xlo_bc = self.subfeatures['x_test_lo']
289 xlo_test = xlo_bc(traj)
290 if xlo_bc.results.uncertain:
291 if self.pars.verbose_level > 0:
292 print "Lo bd uncertain"
293 if self.isdiscrete:
294
295 xlo_test = True
296 else:
297
298 xlo_test = self.subfeatures['dx_test_lo'](traj)
299 xhi_bc = self.subfeatures['x_test_hi']
300 xhi_test = xhi_bc(traj)
301 if xhi_bc.results.uncertain:
302 if self.pars.verbose_level > 0:
303 print "Hi bd uncertain"
304 if self.isdiscrete:
305
306 xhi_test = True
307 else:
308
309 xhi_test = self.subfeatures['dx_test_hi'](traj)
310 for sf in self.subfeatures.values():
311 self.results[sf.name] = sf.results
312 return xlo_test and xhi_test
313
314 - def _find_idx(self):
315 """Helper function for finding lowest index in trajectory meshpoints
316 at which domain test first failed"""
317 if self.results.satisfied:
318
319 return None
320 lowest_idx = Inf
321 for sfname, sf in self.subfeatures.items():
322 if self.pars.verbose_level > 0:
323 print "\n", sfname, sf.results
324 try:
325 res = list(self.results[sfname].output)
326 except AttributeError:
327
328 continue
329 if sf.results.satisfied:
330 continue
331 if self.pars.verbose_level > 0:
332 print res
333
334
335
336 idx = res.index(False)
337 if idx < lowest_idx:
338 lowest_idx = idx
339 return lowest_idx
340
341
342
343
344
346 """
347 General-purpose Hybrid and Non-Hybrid Model abstract class.
348
349 """
350 _needKeys = ['name', 'modelInfo']
351 _optionalKeys = ['ics', 'mspecdict', 'verboselevel',
352 'norm', 'tdata', 'eventPars', 'abseps']
353
354 _querykeys = ['pars', 'parameters', 'events', 'submodels',
355 'ics', 'initialconditions', 'vars', 'variables',
356 'auxvariables', 'auxvars', 'vardomains', 'pardomains',
357 'abseps']
358
359 _setkeys = ['pars', 'algparams', 'checklevel', 'abseps',
360 'ics', 'inputs', 'tdata', 'restrictDSlist',
361 'globalt0', 'verboselevel', 'inputs_t0']
362
364
365
366 if not legit==True:
367
368
369 raise RuntimeError("Only use HybridModel or NonHybridModel classes")
370 if len(a) > 0:
371 if len(a) == 1 and isinstance(a[0], dict):
372 if intersect(a[0].keys(),kw.keys()) != []:
373 raise ValueError("Cannot have initialization keys "
374 "common to both dictionary and keyword arguments")
375 kw.update(a[0])
376 else:
377 raise ValueError("Non-keyword arguments must be a single "
378 "dictionary")
379 try:
380 self.name = kw['name']
381
382
383
384
385
386
387 self.modelInfo = kw['modelInfo']
388 except KeyError:
389 raise KeyError('Necessary keys missing from argument')
390 _foundKeys = len(self._needKeys)
391
392
393
394
395
396
397 self.defaultVars()
398
399
400
401 self.trajectories = {}
402 self.trajectory_defining_args = {}
403
404
405
406 self.registry = {}
407 for name, infodict in self.modelInfo.iteritems():
408
409
410 try:
411 infodict['dsi']._supermodel = self.name
412 except KeyError:
413 raise TypeError("Invalid modelInfo entry found with name %s"%name)
414 self.registry[name] = infodict['dsi'].model
415
416 self.diagnostics = Diagnostics()
417
418
419 self.icdict = {}
420 if 'ics' in kw:
421 self.icdict = dict(kw['ics'])
422 _foundKeys += 1
423 else:
424 self.icdict = {}.fromkeys(self.allvars, NaN)
425 if 'tdata' in kw:
426 self.tdata = kw['tdata']
427 _foundKeys += 1
428 else:
429 self.tdata = None
430 if 'norm' in kw:
431 self._normord = kw['norm']
432 _foundKeys += 1
433 else:
434 self._normord = 2
435 if 'abseps' in kw:
436
437 self._abseps = kw['abseps']
438 _foundKeys += 1
439 else:
440
441 self._abseps = self.query('abseps')
442 if 'verboselevel' in kw:
443 if kw['verboselevel'] in [0,1,2]:
444 self.verboselevel = kw['verboselevel']
445 else:
446 raise ValueError("Verbosity level value must be 0, 1, or 2")
447 _foundKeys += 1
448 else:
449 self.verboselevel = 0
450 if 'mspecdict' in kw:
451 self._mspecdict = kw['mspecdict']
452 _foundKeys += 1
453 else:
454 self._mspecdict = None
455 if 'eventPars' in kw:
456 self._eventPars = kw['eventPars']
457 _foundKeys += 1
458 else:
459 self._eventPars = {}
460 if _foundKeys < len(kw):
461 raise KeyError('Invalid keys found in arguments')
462
463
464
465
466
467
468
469
470
471 self._inputt0_par_links = {}
472
473
475 """Return number of sub-models"""
476 return len(self.registry)
477
479 """Return a list of all sub-model instances (model interfaces or generators)"""
480 return self.registry.values()
481
483 """Return default observable, internal, and auxiliary variable names
484 from modelInfo."""
485 obsvars = []
486 auxvars = []
487 all_known_varnames = []
488 for infodict in self.modelInfo.values():
489 varnames = infodict['dsi'].query('variables')
490 all_known_varnames.extend(varnames)
491 auxvarnames = infodict['dsi'].query('auxvariables')
492 if auxvars == []:
493
494 auxvars.extend(auxvarnames)
495 else:
496 auxvars = intersect(auxvars, auxvarnames)
497 if obsvars == []:
498
499 obsvars = varnames
500 else:
501 obsvars = intersect(obsvars, varnames)
502 intvars = remain(all_known_varnames, obsvars)
503 return (obsvars, intvars, auxvars)
504
506 """Record parameter info locally, for future queries.
507 Internal use only.
508 """
509
510
511
512 self.pars = {}
513 for model in self.registry.values():
514 try:
515 self.pars.update(model.query('pars'))
516 except AttributeError:
517
518 pass
519
520 - def showDef(self, target=None, type=''):
521 """type = 'spec', 'auxspec', 'auxfnspec', 'events', or 'modelspec'
522 (leave blank for the first *four* together).
523
524 'spec', 'auxspec' and 'auxfnspec', 'events' refer to the compiled
525 target language code for the specifications. 'modelspec' refers to
526 the pre-compiled abstract specifications of the model."""
527 if target is None:
528 print "Use showInfo() to find names of defined sub-models"
529 return
530 else:
531 showAll = type==''
532 try:
533 if type=='spec' or showAll:
534 self.registry[target].showSpec()
535 if type=='auxspec' or showAll:
536 self.registry[target].showAuxSpec()
537 if type=='auxfnspec' or showAll:
538 self.registry[target].showAuxFnSpec()
539 if type=='events' or showAll:
540 self.registry[target].showEventSpec()
541 if type=='modelspec':
542 if self._mspecdict is None:
543 raise PyDSTool_ExistError("Cannot use this function "
544 "for models not defined through ModelSpec")
545 info(self._mspecdict[target]['modelspec'].flattenSpec(\
546 [self.modelInfo[target]['dsi'].get('indepvariable').name]))
547 except KeyError:
548 raise ValueError("Model named %s is not known"%target)
549
553
557
561
565
567 return args(pars=self.pars, ics=self.icdict,
568 tdata=self.tdata)
569
571 """Compare self.pars, self.icdict and self.tdata
572 against what's stored for a previously computed trajectory,
573 so that re-computation can be avoided.
574 """
575 try:
576 return info == self.trajectory_defining_args[trajname]
577 except KeyError:
578
579 return False
580
582 foundKeys = 0
583 if 'verboselevel' in kw:
584 self.set(verboselevel=kw['verboselevel'])
585 foundKeys += 1
586 else:
587 self.set(verboselevel=0)
588 if 'ics' in kw:
589 self.icdict = dict(kw['ics'])
590 foundKeys += 1
591 if 'pars' in kw:
592 self.set(pars=kw['pars'])
593 foundKeys += 1
594 if 'tdata' in kw:
595 tdata = kw['tdata']
596 foundKeys += 1
597
598 else:
599 tdata = self.tdata
600
601 if 'force' in kw:
602 force_overwrite = kw['force']
603 foundKeys += 1
604 else:
605 force_overwrite = False
606 if len(kw) != foundKeys:
607 raise PyDSTool_KeyError('Invalid argument keys passed to compute()')
608 if tdata is None:
609 raise ValueError("tdata must be specified")
610 if len(tdata) == 1:
611 assert isinstance(tdata, float) or isinstance(tdata, int), \
612 'tdata must be either a single number or a pair'
613 t0_global = tdata[0]
614 t1_global = Inf
615 elif len(tdata) == 2:
616 t0_global = tdata[0]
617 t1_global = tdata[1]
618 else:
619 raise ValueError('tdata argument key may be either a single '
620 'float or a pair of floats')
621 if not force_overwrite:
622 assert trajname not in self.trajectories, \
623 'Trajectory name already exists'
624 assert self.modelInfo != {}, \
625 'No Generator or Model objects defined for this model'
626 return tdata, t0_global, t1_global, force_overwrite
627
628 - def query(self, querykey=''):
629 """Return info about Model set-up.
630 Valid query key: 'pars', 'parameters', 'events', 'submodels',
631 'ics', 'initialconditions', 'vars', 'variables',
632 'auxvars', 'auxvariables', 'vardomains', 'pardomains', 'abseps'
633 """
634 assert isinstance(querykey, str), \
635 ("Query argument must be a single string")
636 if querykey not in self._querykeys:
637 print 'Valid query keys are:', self._querykeys
638 print "('events' key only queries model-level events, not those"
639 print " inside sub-models)"
640 if querykey != '':
641 raise TypeError('Query key '+querykey+' is not valid')
642 if querykey in ['pars', 'parameters']:
643 result = copy.copy(self.pars)
644 elif querykey in ['ics', 'initialconditions']:
645 result = copy.copy(self.icdict)
646 elif querykey == 'events':
647 result = {}
648 for dsName, model in self.registry.iteritems():
649 try:
650 result.update(model.eventstruct.events)
651 except AttributeError:
652
653 result.update(model.query('events'))
654 elif querykey == 'submodels':
655 result = self.registry
656 elif querykey in ['vars', 'variables']:
657 result = copy.copy(self.allvars)
658 elif querykey in ['vardomains', 'xdomains']:
659 result = {}
660
661 for model in self.registry.values():
662 vardoms = model.query('vardomains')
663 if len(result)==0:
664 result.update(vardoms)
665 else:
666 for vname, vdom in result.iteritems():
667 if vdom.issingleton:
668
669 vdom_lo = vdom.get()
670 vdom_hi = vdom_lo
671 else:
672
673 vdom_lo = vdom[0]
674 vdom_hi = vdom[1]
675 if vardoms[vname][0] < vdom[0]:
676 vdom[0] = vardoms[vname][0]
677 if vardoms[vname][1] > vdom[1]:
678 vdom[1] = vardoms[vname][1]
679 if vdom._abseps < result[vname]._abseps:
680
681
682 result[vname]._abseps = vdom._abseps
683 result[vname] = vdom
684
685 for vname in remain(self.allvars, result.keys()):
686 result[vname] = Interval(vname, float, [-Inf, Inf])
687 elif querykey in ['pardomains', 'pdomains']:
688 result = {}
689
690 for model in self.registry.values():
691 pardoms = model.query('pardomains')
692 if len(result)==0:
693 result.update(pardoms)
694 else:
695 for pname, pdom in result.iteritems():
696 if pdom.issingleton:
697
698 pdom_lo = pdom.get()
699 pdom_hi = pdom_lo
700 else:
701
702 pdom_lo = pdom[0]
703 pdom_hi = pdom[1]
704 if pardoms[pname][0] < pdom[0]:
705 pdom[0] = pardoms[pname][0]
706 if pardoms[pname][1] > pdom[1]:
707 pdom[1] = pardoms[pname][1]
708 if pdom._abseps < result[pname]._abseps:
709
710
711 result[pname]._abseps = pdom._abseps
712 result[pname] = pdom
713 elif querykey in ['auxvars', 'auxvariables']:
714 result = copy.copy(self.auxvars)
715 elif querykey == 'abseps':
716 result = min([ds.query('abseps') for ds in self.registry.values()])
717 return result
718
720 try:
721 return self.modelInfo[dsName]['swRules']
722 except KeyError:
723 raise NameError("Sub-model %s not found in model"%dsName)
724
726
727
728 if isMultiRef(p):
729
730
731
732
733 lbrace = p.find('[')
734 rbrace = p.find(']')
735 if rbrace < lbrace:
736 raise ValueError("Invalid multiple reference to pars")
737 rootname = p[:lbrace]
738 rangespec = p[lbrace+1:rbrace].split(',')
739 try:
740 remainder = p[rbrace+1:]
741 except KeyError:
742
743 remainder = ''
744 if len(rangespec) != 2:
745 raise ValueError("Invalid multiple reference to pars")
746 loix = int(rangespec[0])
747 hiix = int(rangespec[1])
748 if loix >= hiix:
749 raise ValueError("Invalid multiple reference to pars")
750
751
752 for ix in range(loix, hiix+1):
753 self.setPars(rootname+str(ix), val)
754 elif isHierarchicalName(p):
755 if self._mspecdict is None:
756 raise PyDSTool_ExistError("Cannot use this functionality for"
757 " models not defined through ModelSpec")
758
759
760 allFoundNames = []
761 for mspecinfo in self._mspecdict.values():
762 foundNames = searchModelSpec(mspecinfo['modelspec'], p)
763
764 allFoundNames.extend(remain(foundNames,allFoundNames))
765
766
767
768
769 if allFoundNames == []:
770 self.set(pars={p: val})
771 else:
772 self.set(pars={}.fromkeys(allFoundNames, val))
773 else:
774 self.set(pars={p: val})
775
776
778
779
780 if isMultiRef(p):
781
782
783
784
785 lbrace = p.find('[')
786 rbrace = p.find(']')
787 if rbrace < lbrace:
788 raise ValueError("Invalid multiple reference to initial conditions")
789 rootname = p[:lbrace]
790 rangespec = p[lbrace+1:rbrace].split(',')
791 try:
792 remainder = p[rbrace+1:]
793 except KeyError:
794
795 remainder = ''
796 if len(rangespec) != 2:
797 raise ValueError("Invalid multiple reference to initial conditions")
798 loix = int(rangespec[0])
799 hiix = int(rangespec[1])
800 if loix >= hiix:
801 raise ValueError("Invalid multiple reference to initial conditions")
802
803
804 for ix in range(loix, hiix+1):
805 self.setICs(rootname+str(ix), val)
806 elif isHierarchicalName(p):
807 if self._mspecdict is None:
808 raise PyDSTool_ExistError("Cannot use this functionality for"
809 " models not defined through ModelSpec")
810
811
812 allFoundNames = []
813 for mspecinfo in self._mspecdict.values():
814 foundNames = searchModelSpec(mspecinfo['modelspec'], p)
815
816 allFoundNames.extend(remain(foundNames,allFoundNames))
817
818
819
820
821 if allFoundNames == []:
822 self.set(ics={p: val})
823 else:
824 self.set(ics={}.fromkeys(allFoundNames, val))
825 else:
826 self.set(ics={p: val})
827
828
829 - def set(self, **kw):
830 """Set specific parameters of Model. These will get passed on to
831 all Generators/sub-models that support these keys unless the
832 restrictDSList argument is set (only applies to keys algparams,
833 checklevel, and abseps).
834
835 Permitted keys: 'pars', 'algparams', 'checklevel', 'abseps',
836 'ics', 'inputs', 'tdata', 'restrictDSlist',
837 'globalt0', 'verboselevel', 'inputs_t0'
838 """
839 for key in kw:
840 if key not in self._setkeys:
841 raise KeyError("Not a permitted parameter argument: %s"%key + \
842 ". Allowed keys: "+str(self._setkeys))
843 if 'restrictDSlist' in kw:
844 restrictDSlist = kw['restrictDSlist']
845 else:
846 restrictDSlist = []
847
848
849 if 'ics' in kw:
850 self.icdict.update(filteredDict(dict(kw['ics']), self.icdict.keys()))
851 if 'tdata' in kw:
852 self.tdata = kw['tdata']
853 if 'abseps' in kw:
854 self._abseps = kw['abseps']
855 if 'inputs_t0' in kw:
856 for model in self.registry.values():
857
858
859
860
861 t0val_dict = kw['inputs_t0']
862 for inp, val in t0val_dict.items():
863 if isinstance(val, str):
864 try:
865 new_val = self.pars[val]
866 except KeyError:
867 raise ValueError("Input t0 to parameter link "
868 "invalid: no such parameter "+val)
869 t0val_dict[inp] = new_val
870 self._inputt0_par_links.update({inp: val})
871 elif isinstance(val, _num_types):
872 if inp in self._inputt0_par_links:
873
874 del self._inputt0_par_links[inp]
875 else:
876 raise TypeError("Invalid type for input t0 value")
877 try:
878 model.set(inputs_t0=t0val_dict)
879 except AssertionError:
880
881 pass
882 if 'verboselevel' in kw:
883 if kw['verboselevel'] in [0,1,2]:
884 self.verboselevel = kw['verboselevel']
885 else:
886 raise ValueError("Verbosity level value must be 0, 1, or 2")
887 for model in self.registry.values():
888
889 try:
890 model.set(verboselevel=self.verboselevel)
891 except KeyError:
892
893 pass
894 if restrictDSlist == []:
895 restrictDSlist = self.registry.keys()
896
897
898
899
900
901 dsis = self.modelInfo.values()
902 numDSs = len(dsis)
903
904 for key, value in filteredDict(kw, ['ics', 'tdata',
905 'inputs_t0', 'restrictDSlist', 'globalt0'],
906 neg=True).iteritems():
907
908 if isinstance(value, dict):
909
910 entry_err_attr = {}
911 entry_err_val = {}
912 for entrykey, entryval in value.iteritems():
913 entry_err_attr[entrykey] = 0
914 entry_err_val[entrykey] = 0
915 else:
916 entry_err_attr = 0
917 entry_err_val = 0
918
919 for infodict in dsis:
920 callparsDict = {}
921
922 ds = infodict['dsi'].model
923 if hasattr(ds, '_validKeys'):
924 options = ds._validKeys
925 elif hasattr(ds, '_setkeys'):
926 options = ds._setkeys
927 if key in options:
928 if key in ['algparams', 'checklevel', 'abseps'] and \
929 ds.name not in restrictDSlist:
930
931 continue
932 if isinstance(value, dict):
933 for entrykey, entryval in value.iteritems():
934 try:
935 ds.set(**{key:{entrykey:entryval}})
936 except PyDSTool_AttributeError:
937 entry_err_attr[entrykey] += 1
938 except PyDSTool_ValueError:
939 entry_err_val[entrykey] += 1
940 except AssertionError:
941
942 pass
943 else:
944 try:
945 ds.set(**{key: value})
946 except PyDSTool_AttributeError:
947 entry_err_attr += 1
948 except PyDSTool_ValueError:
949 entry_err_val += 1
950 except AssertionError:
951
952 pass
953
954
955 if isinstance(value, dict):
956 for entrykey, entryval in value.iteritems():
957 if entry_err_attr[entrykey] == numDSs:
958 raise PyDSTool_AttributeError('Parameter does not' +\
959 ' exist in any sub-model: %s = %f'%(entrykey,
960 entryval))
961 if entry_err_val[entrykey] == numDSs:
962 raise PyDSTool_ValueError('Parameter value error in' +\
963 ' every sub-model: %s = %f'%(entrykey, entryval))
964 else:
965
966 pass
967 else:
968 if entry_err_attr == numDSs:
969 raise PyDSTool_AttributeError('Parameter does not exist' +\
970 ' in any sub-model: %s'%key)
971 if entry_err_val == numDSs:
972 raise PyDSTool_ValueError('Parameter value error in' +\
973 ' every sub-model: %s'%key)
974 del(entry_err_attr)
975 del(entry_err_val)
976 self._generateParamInfo()
977
979 try:
980 return self.trajectories[trajname]
981 except KeyError:
982 raise ValueError('No such trajectory.')
983
986
988 """Delete a named trajectory from the database."""
989 try:
990 traj = self.trajectories[trajname]
991 except KeyError:
992
993
994
995
996 l = len(trajname)
997 for m in self.registry.values():
998
999 for n in m.trajectories.keys():
1000 if n[:l] == trajname and n[l] == '_' and n[l+1:].isdigit():
1001 m._delTraj(trajname)
1002 else:
1003
1004 if not isinstance(traj.modelNames, str):
1005 for i, model_name_i in enumerate(traj.modelNames):
1006 del(self.registry[model_name_i][trajname+'_%i'%i])
1007 del(self.trajectories[trajname])
1008
1009
1010 - def __call__(self, trajname, t, coords=None, asGlobalTime=True,
1011 asmap=False):
1012 """Evaluate position of hybrid trajectory at time t.
1013 if optional argument asmap == True then t must be an integer in
1014 [0, #segments].
1015 """
1016 try:
1017 traj = self.trajectories[trajname]
1018 except KeyError:
1019 raise ValueError("trajectory '"+trajname+"' unknown")
1020 else:
1021 return traj(t, coords=coords, asGlobalTime=asGlobalTime,
1022 asmap=asmap)
1023
1024
1026 """Returns endpoint of specified trajectory as Point.
1027 trajname: name of selected trajectory
1028 end: (default=1) index of trajectory endpoint.
1029 0 => first, 1 => last
1030 """
1031 xdict = {}
1032 if end not in [0,1]:
1033 raise ValueError("end must be 0, 1")
1034 endtraj = self.trajectories[trajname].trajSeq[-end]
1035 for xname in endtraj.coordnames:
1036 try:
1037 xdict[endtraj._FScompatibleNamesInv(xname)] = \
1038 endtraj.variables[xname].output.datapoints[1][-end]
1039 except KeyError:
1040
1041 pass
1042 except AttributeError:
1043
1044
1045 tend = endtraj.indepdomain[end]
1046 xdict[endtraj._FScompatibleNamesInv(xname)] = \
1047 endtraj.variables[xname](tend)
1048 except PyDSTool_BoundsError:
1049 print "Value out of bounds in variable call:"
1050 print " variable '%s' was called at time %f"%(xname, tend)
1051 raise
1052 return Point({'coorddict': xdict,
1053 'coordnames': endtraj._FScompatibleNamesInv(endtraj.coordnames),
1054 'coordtype': float,
1055 'norm': self._normord})
1056
1058 """Returns end time of specified trajectory.
1059 trajname: name of selected trajectory
1060 end: (default=1) index of trajectory endpoint.
1061 0 => first, 1 => last
1062 """
1063 if end not in [0,1]:
1064 raise ValueError("end must be 0, 1")
1065 endtraj = self.trajectories[trajname].trajSeq[-end]
1066 tend = endtraj.indepdomain[end]
1067 return tend
1068
1070 """Check types and uniqueness of variable names."""
1071 namelist = []
1072 for vname in names:
1073 assert isinstance(vname, str), 'variable name must be a string'
1074 assert vname not in namelist, ('variable names must be unique for'
1075 ' model')
1076 namelist.append(vname)
1077
1078
1080 """Force variables to be the observables in the Model.
1081 May also promote auxiliary variables."""
1082 r = remain(varnames, self.obsvars+self.intvars+self.auxvars)
1083 if len(r) > 0:
1084
1085
1086 raise ValueError("Unknown variable names: "+str(r))
1087 for v in remain(varnames, self.obsvars):
1088
1089 self.obsvars.append(v)
1090
1091 self.intvars = remain(self.intvars, varnames)
1092 self.auxvars = remain(self.auxvars, varnames)
1093 self.obsvars.sort()
1094 self.intvars.sort()
1095 self.auxvars.sort()
1096 self.allvars = self.obsvars + self.intvars
1097 self.allvars.sort()
1098 self.dimension = len(self.allvars)
1099
1100
1102 """Force variables to become internal variables in the Model.
1103 May also promote auxiliary variables."""
1104 r = remain(varnames, self.obsvars+self.intvars+self.auxvars)
1105 if len(r) > 0:
1106
1107
1108 raise ValueError("Unknown variable names: "+str(r))
1109 for v in remain(varnames, self.intvars):
1110
1111 self.intvars.append(v)
1112
1113 self.obsvars = remain(self.obsvars, varnames)
1114 self.auxvars = remain(self.auxvars, varnames)
1115 self.obsvars.sort()
1116 self.intvars.sort()
1117 self.auxvars.sort()
1118 self.allvars = self.obsvars + self.intvars
1119 self.allvars.sort()
1120 self.dimension = len(self.allvars)
1121
1122
1124 """(Re)set to default observable and internal variable names."""
1125 obsvars, intvars, auxvars = self._makeDefaultVarNames()
1126 self._validateVarNames(obsvars + intvars + auxvars)
1127
1128 self.obsvars = obsvars
1129 self.intvars = intvars
1130 self.auxvars = auxvars
1131 self.obsvars.sort()
1132 self.intvars.sort()
1133 self.auxvars.sort()
1134 self.allvars = self.obsvars + self.intvars
1135 self.allvars.sort()
1136 self.dimension = len(self.allvars)
1137
1138
1139 - def info(self, verboselevel=1):
1141
1142
1145
1146 __str__ = __repr__
1147
1148
1150 pickledself = pickle.dumps(self)
1151 return pickle.loads(pickledself)
1152
1153
1155 pickledself = pickle.dumps(self)
1156 return pickle.loads(pickledself)
1157
1158
1159 - def renameTraj(self, trajname, newname, force=False):
1160 """Rename stored trajectory. Force option (default False)
1161 will overwrite any existing trajectory with the new name.
1162 """
1163 try:
1164 traj = self.trajectories[trajname]
1165 except KeyError:
1166 raise ValueError('No such trajectory name %s'%trajname)
1167 if trajname != newname:
1168 if newname not in self.trajectories or force:
1169 self.trajectories[newname] = traj
1170 del self.trajectories[trajname]
1171 traj.name = newname
1172 else:
1173 raise ValueError("Name %s already exists"%newname)
1174
1176 """Return the named trajectory's associated sub-model(s) used
1177 to create it, specific to time t if given."""
1178 try:
1179 modelNames = self.trajectories[trajname].modelNames
1180 except KeyError:
1181 raise ValueError('No such trajectory name %s'%trajname)
1182 if t is None:
1183
1184 return modelNames
1185 else:
1186 parts = self.getTrajTimePartitions(trajname)
1187 pix = 0
1188 for pinterval in parts:
1189 if pinterval.contains(t) is not notcontained:
1190 return modelNames[pix]
1191 else:
1192 pix += 1
1193
1195 """Return the named trajectory's Generator-flagged event times.
1196
1197 events argument can be singleton string name of an event,
1198 returning the event data, or events can be a list of event names,
1199 returning a dictionary of event name -> event data.
1200 Event names should use hierarchical naming convention, if
1201 applicable."""
1202 try:
1203 return self.trajectories[trajname].getEventTimes(events)
1204 except KeyError:
1205 raise
1206
1207
1209 """Return the named trajectory's Generator-flagged events.
1210
1211 events argument can be singleton string name of an event,
1212 returning the event data, or events can be a list of event names,
1213 returning a dictionary of event name -> event data.
1214 Event names should use hierarchical naming convention, if
1215 applicable."""
1216 try:
1217 return self.trajectories[trajname].getEvents(events)
1218 except KeyError:
1219 raise ValueError('No such trajectory name')
1220
1222 """Return the named trajectory's model event structure (representing
1223 external constraints), as present when it was computed.
1224 """
1225 try:
1226 return self.trajectories[trajname].modelEventStructs
1227 except KeyError:
1228 raise ValueError('No such trajectory name')
1229
1231 """Return the named trajectory's time domain,
1232 over which it is defined, in a single interval."""
1233 try:
1234 return self.trajectories[trajname].indepdomain
1235 except KeyError:
1236 raise ValueError('No such trajectory name')
1237
1239 """Return the named trajectory's time domain,
1240 over which it is defined, as a list of time interval partitions."""
1241 try:
1242 return self.trajectories[trajname].timePartitions
1243 except KeyError:
1244 raise ValueError('No such trajectory name')
1245
1247 """
1248 Returns value of given algorithmic parameter for selected sub-model.
1249 target -- name of sub-model in model (cannot be list).
1250 par -- name of algorithmic parameter.
1251 idx -- (optional) index into value if algorithmic parameter val is a
1252 list of values.
1253 """
1254 if target in self.registry.keys():
1255 algpars = self.registry[target].get('algparams')
1256 if par in algpars.keys():
1257 if isinstance(algpars[par], list):
1258 if idx is not None:
1259 if isinstance(idx, list):
1260 val = [algpars[par][x] for x in idx]
1261 else:
1262 val = algpars[par][idx]
1263 else:
1264 val = algpars[par]
1265 else:
1266 val = algpars[par]
1267 else:
1268 val = None
1269 else:
1270 raise ValueError("Target sub-model name not found")
1271 return val
1272
1274 """
1275 Set value of algorithmic parameter in a specific generator.
1276 target -- name or list of generators in model.
1277 par -- name of algorithmic parameter is to be set.
1278 val -- value to which the algorithmic parameter is to be set.
1279
1280 if target is a list, then algorithmic pararameter 'par' is
1281 set to 'val' for every generator in the list, if par exists for that
1282 generator.
1283
1284 WARNING: THIS FUNCTION IS NOT 'SAFE' -- IT DOES NOT CHECK THAT VALS
1285 ARE APPROPRIATE TO PARAMETERS!!!
1286 """
1287
1288 if isinstance(target, list):
1289 subModelList = target
1290 else:
1291 subModelList = [target]
1292 for dsName in subModelList:
1293 algpars = self.registry[dsName].get('algparams')
1294 if par in algpars.keys():
1295
1296
1297 if isinstance(algpars[par], list):
1298 if isinstance(val, list):
1299 if len(algpars[par]) != len(val):
1300 print "Warning: par %s list len (%d) in generator %s doesn't match val list len (%d). Skipping."%(par, len(algpars[par]), dsName, len(val))
1301 continue
1302 else:
1303 algpars[par] = val
1304 else:
1305
1306 for x in range(len(algpars[par])):
1307 algpars[par][x] = val
1308 else:
1309 if isinstance(val, list):
1310 print "Warning: par %s type (%s) in generator %s doesn't match val type (%s). Skipping."%(par, type(algpars[par]), dsName, type(val))
1311 else:
1312 algpars[par] = val
1313 else:
1314 pass
1315
1316
1317 - def sample(self, trajname, coords=None, dt=None,
1318 tlo=None, thi=None, doEvents=True, precise=False):
1319 """Uniformly sample the named trajectory over range indicated,
1320 including any event points. (e.g. use this call for plotting purposes.)
1321 Outputs a Pointset from the trajectory over a given range.
1322
1323 Arguments:
1324
1325 trajname Name of stored trajectory to sample
1326 coords (optional) list of variable names to include in output
1327 dt (optional) step size to use in sampling of independent variable
1328 If not given, the underlying time mesh is used, if available.
1329 tlo (optional) Start value for independent variable, default 0
1330 thi (optional) End value for independent variable, default last value
1331 doEvents (optional) include any event points in output, default True
1332 precise (optional) The default value, False, causes an attempt to use
1333 the underlying mesh of the trajectory to return a Pointset more
1334 quickly. Currently, this can only be used for trajectories that
1335 have a single segment (non-hybrid).
1336 """
1337 try:
1338 return self.trajectories[trajname].sample(coords, dt, tlo, thi,
1339 doEvents, precise)
1340 except KeyError:
1341 raise
1342
1343
1344 - def _infostr(self, dsName=None, verbosity=0):
1345 """Return string information about named sub-model (if given by dsName)
1346 at given verbosity level (default 0)"""
1347 pp = pprint.PrettyPrinter(indent=3)
1348 if dsName is None:
1349 result = {}
1350 res_str = "Sub-models defined in model %s:\n" % self.name
1351 for name, infodict in self.modelInfo.iteritems():
1352 result[name] = infodict['dsi'].model._infostr(verbosity)
1353 return res_str + pp.pprint(result)
1354 else:
1355
1356 try:
1357 result = {dsName: self.modelInfo[dsName]['dsi'].model._infostr(verbosity)+"\n"+\
1358 "Event mapping info:\n" + str(self.modelInfo[dsName]['swRules'])}
1359 except KeyError:
1360 raise NameError("Sub-model %s not found in model"%dsName)
1361 return "Sub-model %s:\n" % (dsName) + pp.pprint(result)
1362
1364
1365 if ics is None:
1366 ics = self.icdict
1367 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1368 estruct(verbosity)
1369
1371 """
1372 List of events in target which are terminal/non-terminal according to
1373 value of flagVal.
1374 target -- name of generators in model (cannot be list)
1375 flagVal -- True (terminal) or False (non-terminal)
1376 """
1377 if ics is None:
1378 ics = self.icdict
1379 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1380 if isinstance(target, str):
1381 if flagVal:
1382 return [x[0] for x in estruct.getTermEvents()]
1383 else:
1384 return [x[0] for x in estruct.getNonTermEvents()]
1385 else:
1386 return []
1387
1388 - def setDSEventTerm(self, target, eventTarget, flagVal, ics=None, t=0):
1389 """
1390 Set event in a specific generator to be (non)terminal.
1391 target -- name or list of names of generators in model.
1392 eventTarget -- name or list of names of events in specified generator(s)
1393 flagVal -- True (event terminal) or False (event non-terminal)
1394
1395 flagVal is applied to every listed event in every listed generator, if
1396 events and generators exist.
1397 """
1398 if ics is None:
1399 ics = self.icdict
1400 if isinstance(target, list):
1401 for targName in target:
1402 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1403 estruct.setTermFlag(eventTarget, flagVal)
1404 else:
1405 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1406 estruct.setTermFlag(eventTarget, flagVal)
1407
1409 """
1410 List of events in target which are active/inactive according to
1411 value of flagVal.
1412 target -- name of generators in model (cannot be list)
1413 flagVal -- True (active) or False (inactive)
1414 """
1415 if ics is None:
1416 ics = self.icdict
1417 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1418 if isinstance(target, str):
1419 if flagVal:
1420 return [x[0] for x in estruct.getActiveEvents()]
1421 else:
1422 return [x[0] for x in estruct.getNonActiveEvents()]
1423 else:
1424 return []
1425
1427 """
1428 Set event in a specific generator to be (in)active.
1429 target -- name or list of names of generators in model.
1430 eventTarget -- name or list of names of events in specified generator(s)
1431 flagVal -- True (event active) or False (event inactive)
1432
1433 flagVal is applied to every listed event in every listed generator, if
1434 events and generators exist.
1435 """
1436 if ics is None:
1437 ics = self.icdict
1438 if isinstance(target, list):
1439 for targName in target:
1440 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1441 estruct.setActiveFlag(eventTarget, flagVal)
1442 else:
1443 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1444 estruct.setActiveFlag(eventTarget, flagVal)
1445
1446 - def setDSEventICs(self, target, eventTarget, val, ics=None, t=0):
1447 """
1448 Set event in a specific generator to have specified initial conditions.
1449 target -- name or list of names of generators in model.
1450 eventTarget -- name or list of names of events in specified generator(s)
1451 val -- dictionary of varnames and initial condition values
1452
1453 val is applied to every listed event in every listed generator, if
1454 events and generators exist.
1455 If a varname in the val dict does not exist in a specified
1456 generator/event, then it is skipped.
1457 """
1458 if ics is None:
1459 ics = self.icdict
1460 if isinstance(target, list):
1461 for targName in target:
1462 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1463 estruct.setICs(eventTarget, val)
1464 else:
1465 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1466 estruct.setICs(eventTarget, val)
1467
1469 """
1470 Set event in a specific generator to be (im)precise.
1471 target -- name or list of names of generators in model.
1472 eventTarget -- name or list of names of events in specified generator(s)
1473 flagVal -- True (event precise) or False (event imprecise)
1474
1475 flagVal is applied to every listed event in every listed generator, if events and generators exist.
1476 """
1477 if ics is None:
1478 ics = self.icdict
1479 if isinstance(target, list):
1480 for targName in target:
1481 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1482 estruct.setPreciseFlag(eventTarget, flagVal)
1483 else:
1484 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1485 estruct.setPreciseFlag(eventTarget, flagVal)
1486
1488 """
1489 Set event in a specific generator to have specified event delay.
1490 target -- name or list of names of generators in model.
1491 eventTarget -- name or list of names of events in specified generator(s)
1492 val -- event delay (int, float >= 0)
1493
1494 val is applied to every listed event in every listed generator, if events and generators exist.
1495 """
1496 if ics is None:
1497 ics = self.icdict
1498 if isinstance(target, list):
1499 for targName in target:
1500 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1501 estruct.setEventDelay(eventTarget, val)
1502 else:
1503 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1504 estruct.setEventDelay(eventTarget, val)
1505
1507 """
1508 Set event in a specific generator to have specified event interval.
1509 target -- name or list of names of generators in model.
1510 eventTarget -- name or list of names of events in specified generator(s)
1511 val -- event interval (int, float >= 0)
1512
1513 val is applied to every listed event in every listed generator, if events and generators exist.
1514 """
1515 if ics is None:
1516 ics = self.icdict
1517 if isinstance(target, list):
1518 for targName in target:
1519 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1520 estruct.setEventInterval(eventTarget, val)
1521 else:
1522 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1523 estruct.setEventInterval(eventTarget, val)
1524
1525 - def setDSEventTol(self, target, eventTarget, val, ics=None, t=0):
1526 """
1527 Set event in a specific generator to have specified event tolerance.
1528 target -- name or list of names of generators in model.
1529 eventTarget -- name or list of names of events in specified generator(s)
1530 val -- event tolerance (int, float >= 0)
1531
1532 val is applied to every listed event in every listed generator, if events and generators exist.
1533 """
1534 if ics is None:
1535 ics = self.icdict
1536 if isinstance(target, list):
1537 for targName in target:
1538 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1539 estruct.setEventTol(eventTarget, val)
1540 else:
1541 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1542 estruct.setEventTol(eventTarget, val)
1543
1544 - def setDSEventDir(self, target, eventTarget, val, ics=None, t=0):
1545 """
1546 Set event in a specific generator to have specified direction code.
1547 target -- name or list of names of generators in model.
1548 eventTarget -- name or list of names of events in specified generator(s)
1549 val -- direction code (-1: decreasing, 1: increasing, 0: either direction)
1550
1551 val is applied to every listed event in every listed generator, if events and generators exist.
1552 """
1553 if ics is None:
1554 ics = self.icdict
1555 if isinstance(target, list):
1556 for targName in target:
1557 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1558 estruct.setEventDir(eventTarget, val)
1559 else:
1560 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1561 estruct.setEventDir(eventTarget, val)
1562
1564 """
1565 Set event in a specific generator to have specified start time.
1566 target -- name or list of names of generators in model.
1567 eventTarget -- name or list of names of events in specified generator(s)
1568 val -- start time (float, int)
1569
1570 val is applied to every listed event in every listed generator, if events and generators exist.
1571 """
1572 if ics is None:
1573 ics = self.icdict
1574 if isinstance(target, list):
1575 for targName in target:
1576 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1577 estruct.setStartTime(eventTarget, val)
1578 else:
1579 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1580 estruct.setStartTime(eventTarget, val)
1581
1583 """
1584 Set event in a specific generator to have specified bisect limit.
1585 target -- name or list of names of generators in model.
1586 eventTarget -- name or list of names of events in specified generator(s)
1587 val -- bisect limit (int > 0)
1588
1589 val is applied to every listed event in every listed generator, if events and generators exist.
1590 """
1591 if ics is None:
1592 ics = self.icdict
1593 if isinstance(target, list):
1594 for targName in target:
1595 estruct = self.modelInfo[targName]['dsi'].get('eventstruct', ics, t)
1596 estruct.setBisect(eventTarget, val)
1597 else:
1598 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t)
1599 estruct.setBisect(eventTarget, val)
1600
1605
1607 """Internal method to set all sub-models flag for being part of a hybrid
1608 trajectory computation.
1609 Useful to indicate that high level events should not be reset when a
1610 Generator is reused between hybrid trajectory segments."""
1611 for ds in self.registry.values():
1612 ds._set_for_hybrid_DS(state)
1613
1615 """Find parameter / variables names matching template in
1616 the model's generators or sub-models: the returned result is a list."""
1617 if self._mspecdict is None:
1618 raise PyDSTool_ExistError("Cannot use this function for models"
1619 " not defined through ModelSpec")
1620 result = {}
1621 for modelName, mspecinfo in self._mspecdict.iteritems():
1622
1623 if isinstance(mspecinfo, dict):
1624 foundNames = searchModelSpec(mspecinfo['modelspec'], template)
1625 else:
1626
1627 foundNames = searchModelSpec(mspecinfo.modelspec, template)
1628 result[modelName] = foundNames
1629 return result
1630
1631
1633 """Find variable and auxiliary variable names that have to be in every
1634 generator or sub-model: the returned result is a list."""
1635 if self._mspecdict is None:
1636 raise PyDSTool_ExistError("Cannot use this function for models"
1637 " not defined through ModelSpec")
1638
1639
1640 a_ds_name = self._mspecdict.keys()[0]
1641 a_ds_mspec = self._mspecdict[a_ds_name]['modelspec']
1642 foundNames = searchModelSpec(a_ds_mspec, template)
1643 try:
1644 fspec = self.registry[a_ds_name].get('funcspec')
1645 except AttributeError:
1646
1647 fspec = self.registry[a_ds_name].registry.values()[0].get('funcspec')
1648
1649
1650
1651 return intersect(self.registry[a_ds_name]._FScompatibleNamesInv(fspec.vars + fspec.auxvars),
1652 foundNames)
1653
1654
1655
1656
1663
1664 - def _findTrajInitiator(self, end_reasons, partition_num, t0, xdict,
1665 gi=None, swRules=None):
1666
1667
1668
1669
1670 infodict = self.modelInfo.values()[0]
1671 return infodict['dsi'], infodict['swRules'], \
1672 infodict['globalConRules'], infodict['dsi'].model.name, \
1673 True, None, True
1674
1676 """Clean up memory usage from past runs of a solver that is interfaced through
1677 a dynamic link library. This will prevent the 'continue' integration option from
1678 being accessible and will delete other data about the last integration run."""
1679 self.registry.values()[0].cleanupMemory()
1680
1682 """Returns True iff all objects in modelInfo have
1683 defined Jacobians."""
1684 return self.registry.values()[0].haveJacobian()
1685
1687 """Returns True iff all objects in modelInfo have
1688 defined Jacobians."""
1689 return self.registry.values()[0].haveJacobian_pars()
1690
1691 - def Rhs(self, t, xdict, pdict=None, asarray=False):
1692 """Direct access to a generator's Rhs function.
1693 Parameters:
1694
1695 t time (can use 0 for an autonomous system)
1696 xdict state dictionary or Point.
1697 pdict parameter dictionary or Point
1698 (optional, default current parameters)
1699 asarray (Bool, optional, default False) If true, will return an array
1700 in state name alphabetical order, else a Point
1701 """
1702
1703 ds = self.registry.values()[0]
1704 fscm = ds._FScompatibleNames
1705 fscmInv = ds._FScompatibleNamesInv
1706 vars = ds.get('funcspec').vars
1707 x_fs = filteredDict(fscm(xdict), vars)
1708
1709
1710 old_ics = ds.initialconditions.copy()
1711 ds.initialconditions.update(x_fs)
1712 x_fs.update(ds.initialconditions)
1713 if pdict is None:
1714 pdict = self.pars
1715 if asarray:
1716
1717
1718 rhsval = array(ds.Rhs(t, x_fs, pdict))[fscm.reorder()]
1719 else:
1720 rhsval = Point({'coorddict': dict(zip(fscmInv(vars),
1721 ds.Rhs(t, x_fs, pdict))),
1722 'coordtype': float,
1723 'norm': self._normord})
1724 ds.initialconditions.update(old_ics)
1725 return rhsval
1726
1727 - def Jacobian(self, t, xdict, pdict=None, asarray=False):
1728 """Direct access to a generator's Jacobian function (if defined).
1729
1730 Arguments:
1731
1732 t time (can use 0 for an autonomous system)
1733 xdict state dictionary or Point.
1734 pdict parameter dictionary or Point
1735 (optional, default current parameters)
1736 asarray (Bool, optional, default False) If true, will return an array
1737 in state name alphabetical order, else a Point
1738 """
1739 ds = self.registry.values()[0]
1740 if not ds.haveJacobian():
1741 raise PyDSTool_ExistError("Jacobian not defined")
1742 fscm = ds._FScompatibleNames
1743 fscmInv = ds._FScompatibleNamesInv
1744 vars = ds.get('funcspec').vars
1745 x_fs = filteredDict(fscm(xdict), vars)
1746
1747
1748 old_ics = ds.initialconditions.copy()
1749 ds.initialconditions.update(x_fs)
1750 x_fs.update(ds.initialconditions)
1751 if pdict is None:
1752 pdict = self.pars
1753 if asarray:
1754 J = array(ds.Jacobian(t, x_fs, pdict))[fscm.reorder()]
1755 else:
1756 J = Pointset({'coorddict': dict(zip(fscmInv(vars),
1757 ds.Jacobian(t, x_fs, pdict))),
1758 'coordtype': float,
1759 'norm': self._normord})
1760 ds.initialconditions.update(old_ics)
1761 return J
1762
1763 - def JacobianP(self, t, xdict, pdict=None, asarray=False):
1764 """Direct access to a generator's JacobianP function (if defined).
1765
1766 Arguments:
1767
1768 t time (can use 0 for an autonomous system)
1769 xdict state dictionary or Point.
1770 pdict parameter dictionary or Point
1771 (optional, default current parameters)
1772 asarray (Bool, optional, default False) If true, will return an array
1773 in state name alphabetical order, else a Point
1774 """
1775 ds = self.registry.values()[0]
1776 if not ds.haveJacobian_pars():
1777 raise PyDSTool_ExistError("Jacobian w.r.t. pars not defined")
1778 fscm = ds._FScompatibleNames
1779 fscmInv = ds._FScompatibleNamesInv
1780 x_fs = filteredDict(fscm(xdict), ds.get('funcspec').vars)
1781
1782
1783 old_ics = ds.initialconditions.copy()
1784 ds.initialconditions.update(x_fs)
1785 x_fs.update(ds.initialconditions)
1786 if pdict is None:
1787 pdict = self.pars
1788 if asarray:
1789 Jp = array(ds.JacobianP(t, x_fs, pdict))[fscm.reorder()]
1790 else:
1791 Jp = Pointset({'coorddict': dict(zip(fscmInv(ds.get('funcspec').pars),
1792 ds.JacobianP(t, x_fs, pdict))),
1793 'coordtype': float,
1794 'norm': self._normord})
1795 ds.initialconditions.update(old_ics)
1796 return Jp
1797
1798 - def MassMatrix(self, t, xdict, pdict=None, asarray=False):
1799 """Direct access to a generator's MassMatrix function (if defined).
1800
1801 Arguments:
1802
1803 t time (can use 0 for an autonomous system)
1804 xdict state dictionary or Point.
1805 pdict parameter dictionary or Point
1806 (optional, default current parameters)
1807 asarray (Bool, optional, default False) If true, will return an array
1808 in state name alphabetical order, else a Point
1809 """
1810 ds = self.registry.values()[0]
1811 if not ds.haveMass():
1812 raise PyDSTool_ExistError("Mass matrix not defined")
1813 fscm = ds._FScompatibleNames
1814 fscmInv = ds._FScompatibleNamesInv
1815 vars = ds.get('funcspec').vars
1816 x_fs = filteredDict(fscm(xdict), vars)
1817
1818
1819 old_ics = ds.initialconditions.copy()
1820 ds.initialconditions.update(x_fs)
1821 x_fs.update(ds.initialconditions)
1822 if pdict is None:
1823 pdict = self.pars
1824 if asarray:
1825 M = array(ds.MassMatrix(t, x_fs, pdict))[fscm.reorder()]
1826 else:
1827 M = Point({'coorddict': dict(zip(fscmInv(vars),
1828 ds.MassMatrix(t, x_fs, pdict))),
1829 'coordtype': float,
1830 'norm': self._normord})
1831 ds.initialconditions.update(old_ics)
1832 return M
1833
1834 - def AuxVars(self, t, xdict, pdict=None, asarray=False):
1835 """Direct access to a generator's auxiliary variables
1836 definition (if defined).
1837
1838 Arguments:
1839
1840 t time (can use 0 for an autonomous system)
1841 xdict state dictionary or Point.
1842 pdict parameter dictionary or Point
1843 (optional, default current parameters)
1844 asarray (Bool, optional, default False) If true, will return an array
1845 in state name alphabetical order, else a Point
1846 """
1847 fscm = ds._FScompatibleNames
1848 fscmInv = ds._FScompatibleNamesInv
1849 x_fs = filteredDict(fscm(xdict), ds.get('funcspec').vars)
1850
1851
1852 old_ics = ds.initialconditions.copy()
1853 ds.initialconditions.update(x_fs)
1854 x_fs.update(ds.initialconditions)
1855 if pdict is None:
1856 pdict = self.pars
1857 if asarray:
1858 A = array(ds.AuxVars(t, x_fs, pdict))[fscm.reorder()]
1859 else:
1860 A = Point({'coorddict': dict(zip(fscmInv(ds.get('funcspec').auxvars),
1861 ds.AuxVars(t, x_fs, pdict))),
1862 'coordtype': float,
1863 'norm': self._normord})
1864 ds.initialconditions.update(old_ics)
1865 return A
1866
1867
1868 - def compute(self, trajname, **kw):
1869 """Compute a non-hybrid trajectory. Returns a Trajectory object.
1870
1871 Arguments (non-keyword):
1872 trajname Name of trajectory to create (string)
1873
1874 Arguments (keyword only, all optional):
1875 force (Bool, default False) - force overwrite of any trajectory
1876 stored in this object with the same name
1877 verboselevel (int, default 0)
1878 ics initial conditions dict or Point
1879 pars parameters dict or Point
1880 tdata time data (interval as sequence of 2 numeric values)
1881 """
1882 tdata, t0_global, t1_global, force_overwrite = \
1883 Model._prepareCompute(self, trajname, **kw)
1884 if self.icdict == {}:
1885
1886
1887 self.icdict = self.registry.values()[0].get('initialconditions')
1888 if self.icdict == {}:
1889
1890 raise PyDSTool_ExistError("No initial conditions specified")
1891 xdict = {}
1892 for xname, value in self.icdict.iteritems():
1893
1894
1895
1896 xdict[xname] = ensurefloat(value)
1897
1898 self.icdict = xdict.copy()
1899
1900
1901 gen = self.registry.values()[0]
1902
1903 t0_global += gen.globalt0
1904 t1_global += gen.globalt0
1905 t0 = t0_global
1906 if t1_global > gen.indepvariable.indepdomain[1]+t0:
1907 if isfinite(gen.indepvariable.indepdomain[1]):
1908 if self.verboselevel > 0:
1909 print "Warning: end time was truncated to max size of specified independent variable domain"
1910 t1 = gen.indepvariable.indepdomain[1]+t0
1911 else:
1912 t1 = t1_global
1913 else:
1914 t1 = t1_global
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926 setup_pars = {'ics': xdict, 'tdata': [0, t1-t0]}
1927 if self._abseps is not None:
1928 setup_pars['abseps'] = self._abseps
1929 if hasattr(gen, 'algparams'):
1930 setup_pars['algparams'] = {'verbose': self.verboselevel}
1931 try:
1932 if gen.algparams['init_step'] > t1-t0:
1933 if self.verboselevel > 0:
1934 print "Warning: time step too large for remaining time"\
1935 + " interval. Temporarily reducing time step to " \
1936 + "1/10th of its previous value"
1937 setup_pars['algparams'].update({'init_step': (t1-t0)/10})
1938 except (AttributeError, KeyError):
1939
1940 pass
1941
1942
1943 if self._inputt0_par_links != {}:
1944 t0val_dict = {}
1945 for inp, val in self._inputt0_par_links.items():
1946 t0val_dict[inp] = self.pars[val]
1947 gen.set(inputs_t0=t0val_dict)
1948 gen.set(**setup_pars)
1949 gen.diagnostics.clearWarnings()
1950 gen.diagnostics.clearErrors()
1951
1952
1953
1954 try:
1955 traj = gen.compute(trajname+'_0')
1956 except PyDSTool_ValueError, e:
1957 print "\nError in Generator:", gen.name
1958 gen.diagnostics.showWarnings()
1959 gen.diagnostics.showErrors()
1960 print "Are the constituent generator absolute epsilon "
1961 print " tolerances too small? -- this abseps =", gen._abseps
1962 raise
1963 except KeyboardInterrupt:
1964 raise
1965 except:
1966 print "\nError in Generator:", gen.name
1967 gen.diagnostics.showWarnings()
1968 gen.diagnostics.showErrors()
1969 self.diagnostics.traceback = {}
1970 for k,v in gen.diagnostics.traceback.iteritems():
1971 if isinstance(v, dict):
1972 dentry = {}
1973 for vk, vv in v.iteritems():
1974 dentry[gen._FScompatibleNamesInv(vk)] = vv
1975 else:
1976 dentry = v
1977 self.diagnostics.traceback[k] = dentry
1978 if self.diagnostics.traceback != {}:
1979 print "Traceback dictionary copied to model"
1980 raise
1981 if traj is None:
1982 raise ValueError("Generator %s failed to create a trajectory"%gen.name)
1983 else:
1984 if self._abseps is not None:
1985 assert traj.indepdomain._abseps == self._abseps
1986 for v in traj.depdomain.values():
1987 assert v._abseps == self._abseps
1988 if not isparameterized(traj):
1989 raise ValueError("Generator " + gen.name + " produced a " \
1990 + " non-parameterized Trajectory")
1991
1992
1993
1994
1995
1996
1997
1998
1999 if gen.diagnostics.hasErrors() and self.verboselevel > 0:
2000 for e in gen.diagnostics.errors:
2001 print 'Generator ' + gen.name + \
2002 ' in trajectory segment had errors'
2003 gen.diagnostics.showErrors()
2004 end_reasons = ['time']
2005
2006
2007
2008 time_interval = traj.indepdomain
2009 ti_1 = time_interval[1]
2010 if gen.diagnostics.hasWarnings():
2011 if self.verboselevel >= 1:
2012 gen.diagnostics.showWarnings()
2013 for w in gen.diagnostics.warnings:
2014 if w[0] == Generator.W_TERMEVENT or \
2015 w[0] == Generator.W_TERMSTATEBD:
2016 if self.verboselevel > 1:
2017 print "Time:", ti_1+t0
2018 print 'Generator ' + gen.name + \
2019 ' had a terminal event. Details (in local ' + \
2020 'system time) ...\n ' + str(w[1])
2021
2022
2023
2024 if isinstance(w[1][1], list):
2025 if len(w[1][1]) > 1:
2026 if self.verboselevel > 0:
2027 print 'Warning: More than one terminal event found.'
2028 print ' Consider altering event pars.'
2029 end_reasons = [w[1][1][ri] \
2030 for ri in range(len(w[1][1]))]
2031 else:
2032 end_reasons = [w[1][1]]
2033
2034 self.diagnostics.update(gen.diagnostics)
2035
2036
2037
2038
2039
2040
2041
2042 global_end_reasons = {}
2043 global_end_ixs = []
2044 dsi = self.modelInfo.values()[0]['dsi']
2045 for globalDS in self.modelInfo.values()[0]['globalConRules']:
2046 if globalDS(dsi):
2047 global_end_reasons[Inf] = \
2048 globalDS.conditions.collate_results('reasons',
2049 merge_lists=True)
2050 global_end_ixs.append(Inf)
2051 else:
2052
2053
2054
2055 try:
2056 final_ok_idx = globalDS.conditions._find_idx()
2057 except (AttributeError, RuntimeError):
2058 print "Trajectory creation failed..."
2059 globalDS.conditions.results.info()
2060 raise PyDSTool_ExistError("Global consistency checks failed for"
2061 " model interface %s for trajectory %s"%(str(globalDS),
2062 traj.name))
2063 else:
2064 if final_ok_idx is not None:
2065 global_end_ixs.append(final_ok_idx)
2066 global_end_reasons[final_ok_idx] = \
2067 globalDS.conditions.collate_results('reasons',
2068 merge_lists=True)
2069 if len(global_end_reasons) > 0:
2070 smallest_ok_idx = min(global_end_ixs)
2071
2072
2073 if isfinite(smallest_ok_idx):
2074
2075 traj.truncate_to_idx(smallest_ok_idx)
2076
2077 end_reasons = global_end_reasons[smallest_ok_idx]
2078 if self.verboselevel > 1:
2079 print "End reason in %s was %s"%(self.name,str(end_reasons))
2080
2081
2082
2083
2084
2085 if ti_1 > t1_global-t0:
2086 print "Warning: Generator time interval exceeds prescribed limits:"
2087 print " ... ", ti_1, " > ", t1_global-t0
2088 epochEvents = {}
2089 try:
2090 epochEvents.update(gen.getEvents(asGlobalTime=True))
2091 except (AttributeError, TypeError):
2092
2093 pass
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103 self.trajectories[trajname] = traj
2104
2105
2107 """Validate Model's modelInfo attribute."""
2108
2109 assert len(self.modelInfo) == 1, \
2110 "Non-hybrid model must contain a single Generator"
2111 infodict = self.modelInfo.values()[0]
2112 if not isinstance(infodict['dsi'], MProject.GeneratorInterface):
2113 raise TypeError("Must provide a single Generator object"
2114 " wrapped in a GeneratorInterface")
2115
2117 if verbose > 0:
2118 outputStr = 'Non-Hybrid Model '+self.name+" containing components:"
2119 print "Observable variables:", self.obsvars
2120 print "Internal variables:", self.intvars
2121 print "Auxiliary variables:", self.auxvars
2122 name, infodict = self.modelInfo.items()[0]
2123 outputStr += "\n--- Generator: "+name
2124 outputStr += "\n "+infodict['dsi'].model._infostr(verbose-1)
2125 else:
2126 outputStr = 'Non-Hybrid Model '+self.name
2127 return outputStr
2128
2129
2130
2132 """
2133 obsvars specifies the observable
2134 variables for this model, which must be present in all
2135 trajectory segments (regardless of which other variables are
2136 specified in those segments). intvars specifies
2137 non-observable (internal) variables for the model, which
2138 may or may not be present in trajectory segments, depending
2139 on reductions, etc., that are present in the sub-model (e.g.
2140 generator) that determines the trajectory segment.
2141 """
2151
2153
2154
2155
2156
2157
2158
2159
2160
2161 for xname in self.allvars+self.auxvars:
2162
2163
2164 xname_compat = traj._FScompatibleNames(xname)
2165 try:
2166 xdict[xname] = traj.variables[xname_compat].output.datapoints[1][-1]
2167 except KeyError:
2168
2169 pass
2170 except AttributeError:
2171
2172
2173 try:
2174 xdict[xname] = traj.variables[xname_compat](ti_1)
2175 except RuntimeError:
2176
2177 print "**********"
2178 print xname, " dep domain:", traj.variables[xname_compat].depdomain.get()
2179 print traj.depdomain[xname_compat].get()
2180 print "indep domain:", traj.variables[xname_compat].indepdomain.get()
2181 print "Var val at t-eps %f = "%(ti_1-1e-4), \
2182 traj.variables[xname_compat](ti_1-1e-4)
2183 raise
2184 except:
2185 print "**********"
2186 print xname, " dep domain:", traj.variables[xname_compat].depdomain.get()
2187 print traj.depdomain[xname_compat].get()
2188 print traj.trajSeq[-1].variables[xname_compat].depdomain.get()
2189 raise
2190 except PyDSTool_BoundsError:
2191 print "Value out of bounds in variable call:"
2192 print " variable '%s' was called at time %f"%(xname, ti_1)
2193 raise
2194
2195
2196 - def _applyStateMap(self, epochStateMaps, model_interface,
2197 next_model_interface, traj, xdict, t0):
2198
2199 num_maps = len(epochStateMaps)
2200
2201
2202
2203
2204 if num_maps > 1:
2205
2206
2207
2208 for mapix in range(1,num_maps):
2209 if epochStateMaps[0] != epochStateMaps[mapix]:
2210 raise RuntimeError("PyDSTool does not yet allow "
2211 "truly simultaneous events that do not point "
2212 "to the same model with the same mapping")
2213
2214 dsinps = model_interface.get('inputs', xdict, t0)
2215 if dsinps != {}:
2216
2217
2218 extinputs_ic = {}.fromkeys(dsinps)
2219 try:
2220 for k in extinputs_ic:
2221 extinputs_ic[k] = dsinps[k](t0)
2222 except PyDSTool_BoundsError:
2223 print "Cannot proceed - sub-model '" + model_interface.model.name \
2224 + "' had an external input undefined " \
2225 + "at t =", t0
2226 raise
2227 else:
2228 extinputs_ic = {}
2229 eventstruct = model_interface.get('eventstruct', xdict, t0)
2230 pars = model_interface.get('pars', xdict, t0)
2231 pars_temp = traj._FScompatibleNames(pars)
2232 xdict_temp = traj._FScompatibleNames(xdict)
2233 for epochStateMap in epochStateMaps:
2234
2235 epochStateMap.evmapping(xdict_temp, pars_temp,
2236 extinputs_ic, eventstruct, t0)
2237
2238 pars.update(traj._FScompatibleNamesInv(pars_temp))
2239
2240
2241 try:
2242 next_model_interface.set('pars', pars, xdict, t0)
2243 except (PyDSTool_AttributeError, PyDSTool_ValueError):
2244
2245 pass
2246 xdict.update(traj._FScompatibleNamesInv(xdict_temp))
2247
2248 for xname in self.auxvars:
2249 try:
2250 del xdict[xname]
2251 except KeyError:
2252
2253 pass
2254
2256 """Clean up memory usage from past runs of a solver that is interfaced
2257 through a dynamic link library. This will prevent the 'continue'
2258 integration option from being accessible and will delete other data
2259 about the last integration run.
2260 """
2261 for MI in self.registry.values():
2262 MI.model.cleanupMemory()
2263
2264 - def _findTrajInitiator(self, end_reasons, partition_num, t0, xdict,
2265 mi=None, swRules=None):
2266
2267 epochStateMaps = None
2268 notDone = True
2269 reused = False
2270 if end_reasons is None:
2271
2272 assert partition_num == 0, ("end_reasons was None on a "
2273 "non-initial epoch")
2274 xnames = self.icdict.keys()
2275 xnames.sort()
2276
2277 try:
2278 infodict = findTrajInitiator(self.modelInfo, t0, xdict,
2279 self.pars, self.intvars,
2280 self.verboselevel)
2281 except PyDSTool_ValueError, errinfo:
2282 print errinfo
2283 raise PyDSTool_ExistError('No unique eligible Model found:'
2284 ' cannot continue (check active terminal event definitions'
2285 ' or error message above)')
2286 mi = infodict['dsi']
2287 swRules = infodict['swRules']
2288 globalConRules = infodict['globalConRules']
2289 missing = remain(mi.get('variables', ics=xdict, t0=t0).keys(),
2290 xnames+self.auxvars+self.intvars)
2291 if missing != []:
2292 raise AssertionError('Missing initial condition specifications'
2293 ' for %s' % str(missing))
2294 for entry, val in xdict.iteritems():
2295 if not isfinite(val):
2296 print "Warning: %s initial condition for "%mi.model.name \
2297 + str(entry) + " = " + str(val)
2298 if self.verboselevel > 1:
2299 print "\nStarting partition 0"
2300 print "Chose initiator '%s'"%mi.model.name
2301 nextModelName = None
2302 else:
2303
2304 num_reasons = len(end_reasons)
2305
2306
2307 if num_reasons > 1:
2308
2309 candidateNextModel = ""
2310 mapinfos = []
2311 for ri in range(num_reasons):
2312 try:
2313 mapinfos.append(swRules[end_reasons[ri]])
2314 except KeyError:
2315
2316 continue
2317 nextModelName = mapinfos[-1][0]
2318 if candidateNextModel == "":
2319 candidateNextModel = nextModelName
2320 else:
2321 if candidateNextModel != nextModelName:
2322 raise RuntimeError("The multiple reasons "
2323 "for terminating last trajectory segment"
2324 " did not point to the same generator to"
2325 " continue the calculation.")
2326 if len(mapinfos) == 0:
2327
2328 nextModelName = 'terminate'
2329 else:
2330 epochStateMaps = [m[1] for m in mapinfos]
2331 else:
2332 if end_reasons[0] in swRules:
2333 mapinfo = swRules[end_reasons[0]]
2334 nextModelName = mapinfo[0]
2335 epochStateMaps = [mapinfo[1]]
2336 else:
2337 nextModelName = 'terminate'
2338 epochStateMaps = None
2339 if nextModelName == 'terminate':
2340 if self.verboselevel > 0:
2341 print 'Trajectory calculation for Model `'+mi.model.name\
2342 +'` terminated without a DS specified in switch'\
2343 +' rules with which to continue.'
2344 globalConRules = None
2345 notDone = False
2346 elif nextModelName == mi.model.name:
2347
2348
2349
2350
2351
2352
2353
2354 reused = True
2355
2356
2357 globalConRules = self.modelInfo[nextModelName]['globalConRules']
2358 else:
2359 next = self.modelInfo[nextModelName]
2360 mi = next['dsi']
2361 swRules = next['swRules']
2362 globalConRules = next['globalConRules']
2363 if self.verboselevel > 1:
2364 print "\nStarting partition #%i"%partition_num
2365 print "Chose model '%s'"%mi.model.name
2366 return mi, swRules, globalConRules, nextModelName, reused, \
2367 epochStateMaps, notDone
2368
2369 - def _addTraj(self, trajname, trajseq, epochEvents,
2370 modelNames, force_overwrite=False):
2371 """Add a computed trajectory to database."""
2372
2373 if not force_overwrite:
2374 assert trajname not in self.trajectories, \
2375 'Trajectory name already used'
2376
2377
2378
2379
2380
2381 genEvStructs = dict(zip(modelNames,
2382 [self.modelInfo[gn]['dsi'].eventstruct \
2383 for gn in modelNames]))
2384
2385 time_range = [None, None]
2386 indepvarname = ''
2387 genEventTimes = {}
2388 genEvents = {}
2389
2390 FScompatibleNames = symbolMapClass()
2391 FScompatibleNamesInv = symbolMapClass()
2392 for traj_ix, traj in enumerate(trajseq):
2393 assert isinstance(traj, Trajectory), \
2394 ('traj must contain trajectories')
2395 if indepvarname == '':
2396 indepvarname = traj.indepvarname
2397
2398
2399
2400
2401
2402 FScompatibleNames.update(traj._FScompatibleNames)
2403 FScompatibleNamesInv.update(traj._FScompatibleNamesInv)
2404 for evname in epochEvents[traj_ix]:
2405 if evname in genEvents and genEvents[evname] is not None:
2406 val = mapNames(traj._FScompatibleNamesInv,
2407 epochEvents[traj_ix][evname])
2408 if val is not None:
2409 try:
2410 genEvents[evname].append(val)
2411 except ValueError:
2412
2413
2414 if genEventTimes[evname][-1] - val['t'][0] > 0:
2415
2416 raise
2417 else:
2418
2419 pass
2420 else:
2421 try:
2422 genEventTimes[evname].extend( \
2423 val.indepvararray.tolist())
2424 except AttributeError:
2425 pass
2426 else:
2427 genEvents[evname] = mapNames(traj._FScompatibleNamesInv,
2428 epochEvents[traj_ix][evname])
2429 try:
2430 genEventTimes[evname] = \
2431 epochEvents[traj_ix][evname].indepvararray.tolist()
2432 except AttributeError:
2433 pass
2434 if time_range[0] is None:
2435
2436 DS_t_numtype = traj.indepdomain.type
2437 if isinputcts(traj):
2438 DS_t_typestr = 'continuous'
2439
2440 if not compareNumTypes(DS_t_numtype, _all_float):
2441 raise TypeError('continuous time inconsistent with '
2442 'non-float type')
2443 else:
2444 DS_t_typestr = 'discrete'
2445
2446 traj_vars = traj._FScompatibleNamesInv(traj.variables)
2447 assert intersect(self.allvars, traj_vars.keys()), \
2448 "variable name for traj not present in sub-model's variables"
2449
2450
2451
2452 varname = self.obsvars[0]
2453 if varname in traj_vars:
2454 DS_x_numtype = traj_vars[varname].depdomain.type
2455 else:
2456 raise ValueError("varname not known")
2457
2458
2459
2460 time_range = traj.indepdomain.get()
2461 try:
2462 time_range[0] += traj.globalt0
2463 time_range[1] += traj.globalt0
2464 except IndexError:
2465 raise ValueError('time interval of segment in hybrid '
2466 'trajectory cannot be singleton')
2467
2468
2469
2470
2471
2472
2473 time_partitions = [(traj.indepdomain, \
2474 traj.globalt0, \
2475 traj.checklevel)]
2476 else:
2477
2478
2479
2480 if not compareNumTypes(DS_t_numtype, traj.indepdomain.type):
2481 raise TypeError('Mismatched time types for hybrid '
2482 'trajectory sequence')
2483
2484
2485
2486
2487
2488
2489 if not traj.indepdomain.atEndPoint(time_range[1] \
2490 - traj.globalt0, 'lo'):
2491
2492
2493
2494
2495 raise ValueError('Hybrid trajectory sequence time intervals'
2496 ' must be contiguous: ' + \
2497 str(traj.indepdomain[0]) + ' vs. ' + \
2498 str(time_range[1]-traj.globalt0))
2499 time_range[1] += traj.indepdomain[1]
2500 time_partitions.append((traj.indepdomain, \
2501 traj.globalt0, \
2502 traj.checklevel))
2503
2504 time_interval = Interval(indepvarname, DS_t_numtype, time_range,
2505 abseps=self._abseps)
2506
2507
2508
2509 self.trajectories.update({trajname: \
2510 HybridTrajectory(trajname, trajseq,
2511 timePartitions=time_partitions,
2512 timeInterval=time_interval,
2513 eventTimes=genEventTimes,
2514 events=genEvents,
2515 modelEventStructs=genEvStructs,
2516 modelNames=modelNames,
2517 FScompatibleNames=FScompatibleNames,
2518 FScompatibleNamesInv=FScompatibleNamesInv,
2519 abseps=self._abseps,
2520 globalt0=0, norm=self._normord)
2521 })
2522
2524 """Returns True iff all objects in modelInfo have
2525 defined Jacobians."""
2526 result = True
2527 for model in self.registry.values():
2528 result = result and model.haveJacobian()
2529 return result
2530
2532 """Returns True iff all objects in modelInfo have
2533 defined Jacobians."""
2534 result = True
2535 for model in self.registry.values():
2536 result = result and model.haveJacobian_pars()
2537 return result
2538
2539 - def Rhs(self, dsName, t, xdict, pdict=None, asarray=False):
2540 """Direct access to a sub-model generator's Rhs function.
2541 Arguments:
2542
2543 dsName Name of a sub-model
2544 t time (can use 0 for an autonomous system)
2545 xdict state dictionary or Point.
2546 pdict parameter dictionary or Point
2547 (optional, default current parameters)
2548 asarray (Bool, optional, default False) If true, will return an array
2549 in state name alphabetical order, else a Point
2550 """
2551 try:
2552 dsi = self.modelInfo[dsName]['dsi']
2553 except KeyError:
2554 raise ValueError("No DS named %s was found"%dsName)
2555 if pdict is None:
2556 pdict = self.pars
2557 if asarray:
2558 return dsi.Rhs(t, xdict, pdict, asarray=True)
2559 else:
2560
2561 varnames = dsi.get('funcspec', xdict, t).vars
2562 return Point({'coorddict': dict(zip(varnames,
2563 dsi.Rhs(t, xdict, pdict))),
2564 'coordtype': float,
2565 'norm': self._normord})
2566
2567 - def Jacobian(self, dsName, t, xdict, pdict=None, asarray=False):
2568 """Direct access to a sub-model generator's Jacobian function (if defined).
2569 Arguments:
2570
2571 dsName Name of a sub-model
2572 t time (can use 0 for an autonomous system)
2573 xdict state dictionary or Point.
2574 pdict parameter dictionary or Point
2575 (optional, default current parameters)
2576 asarray (Bool, optional, default False) If true, will return an array
2577 in state name alphabetical order, else a Point
2578 """
2579 try:
2580 dsi = self.modelInfo[dsName]['dsi']
2581 except KeyError:
2582 raise ValueError("No DS named %s was found"%dsName)
2583 if pdict is None:
2584 pdict = self.pars
2585 if dsi.haveJacobian():
2586 if asarray:
2587 return dsi.Jacobian(t, xdict, pdict, asarray=True)
2588 else:
2589 varnames = dsi.get('funcspec', xdict, t).vars
2590 return Pointset({'coorddict': dict(zip(varnames,
2591 dsi.Jacobian(t, xdict, pdict))),
2592 'coordtype': float,
2593 'norm': self._normord})
2594 else:
2595 raise PyDSTool_ExistError("Jacobian not defined")
2596
2597 - def JacobianP(self, dsName, t, xdict, pdict=None, asarray=False):
2598 """Direct access to a generator's JacobianP function (if defined).
2599 Arguments:
2600
2601 dsName Name of a sub-model
2602 t time (can use 0 for an autonomous system)
2603 xdict state dictionary or Point.
2604 pdict parameter dictionary or Point
2605 (optional, default current parameters)
2606 asarray (Bool, optional, default False) If true, will return an array
2607 in state name alphabetical order, else a Point
2608 """
2609 try:
2610 dsi = self.modelInfo[dsName]['dsi']
2611 except KeyError:
2612 raise ValueError("No DS named %s was found"%dsName)
2613 if dsi.haveJacobian_pars():
2614 if pdict is None:
2615 pdict = self.pars
2616 if asarray:
2617 return dsi.JacobianP(t, xdict, pdict, asarray=True)
2618 else:
2619 parnames = dsi.get('funcspec', xdict, t).pars
2620 return Pointset({'coorddict': dict(zip(parnames,
2621 dsi.JacobianP(t, xdict, pdict))),
2622 'coordtype': float,
2623 'norm': self._normord})
2624 else:
2625 raise PyDSTool_ExistError("Jacobian w.r.t. pars not defined")
2626
2627 - def MassMatrix(self, dsName, t, xdict, pdict=None, asarray=False):
2628 """Direct access to a generator's MassMatrix function (if defined).
2629 Arguments:
2630
2631 dsName Name of a sub-model
2632 t time (can use 0 for an autonomous system)
2633 xdict state dictionary or Point.
2634 pdict parameter dictionary or Point
2635 (optional, default current parameters)
2636 asarray (Bool, optional, default False) If true, will return an array
2637 in state name alphabetical order, else a Point
2638 """
2639 try:
2640 dsi = self.modelInfo[dsName]['dsi']
2641 except KeyError:
2642 raise ValueError("No DS named %s was found"%dsName)
2643 if pdict is None:
2644 pdict = self.pars
2645 if asarray:
2646 dsi.MassMatrix(t, xdict, pdict, asarray=True)
2647 else:
2648 varnames = dsi.get('funcspec', xdict, t).vars
2649 return Point({'coorddict': dict(zip(varnames,
2650 dsi.MassMatrix(t, xdict, pdict))),
2651 'coordtype': float,
2652 'norm': self._normord})
2653
2654 - def AuxVars(self, dsName, t, xdict, pdict=None, asarray=False):
2655 """Direct access to a generator's auxiliary variables
2656 definition (if defined).
2657
2658 Arguments:
2659
2660 dsName Name of a sub-model
2661 t time (can use 0 for an autonomous system)
2662 xdict state dictionary or Point.
2663 pdict parameter dictionary or Point
2664 (optional, default current parameters)
2665 asarray (Bool, optional, default False) If true, will return an array
2666 in state name alphabetical order, else a Point
2667 """
2668 try:
2669 dsi = self.modelInfo[dsName]['dsi']
2670 except KeyError:
2671 raise ValueError("No DS named %s was found"%dsName)
2672 if pdict is None:
2673 pdict = self.pars
2674 if asarray:
2675 return dsi.AuxVars(t, xdict, pdict, asarray=True)
2676 else:
2677 auxvarnames = dsi.get('funcspec', xdict, t).auxvars
2678 return Point({'coorddict': dict(zip(auxvarnames,
2679 dsi.AuxVars(t, xdict, pdict))),
2680 'coordtype': float,
2681 'norm': self._normord})
2682
2683 - def compute(self, trajname, **kw):
2684 """Compute a hybrid trajectory and store it internally in the 'trajectories'
2685 attribute.
2686
2687 Arguments (non-keyword):
2688 trajname Name of trajectory to create (string)
2689
2690 Arguments (keyword only, all optional):
2691 force (Bool, default False) - force overwrite of any trajectory
2692 stored in this object with the same name
2693 verboselevel (int, default 0)
2694 ics initial conditions dict or Point
2695 pars parameters dict or Point
2696 tdata time data (interval as sequence of 2 numeric values)
2697 """
2698
2699
2700
2701
2702
2703 tdata, t0_global, t1_global, force_overwrite = \
2704 Model._prepareCompute(self, trajname, **kw)
2705
2706
2707 end_reasons = None
2708
2709 if self.icdict == {}:
2710 raise PyDSTool_ExistError("No initial conditions specified")
2711
2712 xdict = {}
2713 for xname, value in self.icdict.iteritems():
2714 if xname not in self.allvars:
2715 raise ValueError("Invalid variable name in initial "
2716 "conditions: " + xname)
2717 xdict[xname] = ensurefloat(value)
2718
2719 self.icdict = xdict.copy()
2720
2721
2722 notDone = True
2723 t0 = t0_global
2724 partition_num = 0
2725 trajseq = []
2726 traj = None
2727 ti_1 = None
2728 modelNames = []
2729 epochEvents = []
2730 MI_prev = None
2731 MI = None
2732 swRules = None
2733
2734 reused = False
2735
2736 self.resetEventTimes()
2737
2738
2739
2740
2741 while notDone:
2742
2743
2744
2745 MI, swRules, globalConRules, nextModelName, reused, \
2746 epochStateMaps, notDone = self._findTrajInitiator(end_reasons,
2747 partition_num, t0, xdict,
2748 MI, swRules)
2749 if not notDone:
2750 continue
2751
2752
2753 model = MI.model
2754
2755 if partition_num > 0:
2756
2757
2758
2759
2760 pars_copy = copy.copy(self.pars)
2761 self._applyStateMap(epochStateMaps, MI_prev, MI, traj, xdict, t0)
2762
2763 for event, delay in event_delay_record.items():
2764 event.eventdelay = delay
2765
2766
2767
2768
2769
2770 event_delay_record = {}
2771 dummy1, dummy2, targetMI = MI._get_initiator_cache(xdict, t0)
2772 evstruct = targetMI.get('eventstruct', xdict, t0)
2773
2774 for reason in end_reasons:
2775
2776 try:
2777 event = evstruct.events[reason]
2778 except KeyError:
2779 pass
2780 else:
2781 event_delay_record[event] = event.eventdelay
2782 event.eventdelay = event.eventinterval
2783 else:
2784
2785
2786 pars_copy = {}
2787 event_delay_record = {}
2788
2789
2790
2791 if self._inputt0_par_links != {} and pars_copy != self.pars:
2792 t0val_dict = {}
2793 for inp, val in self._inputt0_par_links.items():
2794 t0val_dict[inp] = self.pars[val]
2795 model.set(inputs_t0=t0val_dict)
2796
2797
2798 indepvar = MI.get('indepvariable', xdict, t0)
2799 if t1_global > indepvar.indepdomain[1]+t0:
2800 if isfinite(indepvar.indepdomain[1]):
2801 if self.verboselevel > 0:
2802 print "Warning: end time was truncated to max size " + \
2803 "of specified independent variable domain"
2804 t1 = indepvar.indepdomain[1]+t0
2805 else:
2806 t1 = t1_global
2807 else:
2808 t1 = t1_global
2809
2810
2811
2812
2813
2814
2815
2816
2817
2818 MI.set('tdata', [0, t1-t0], xdict, t0)
2819 MI.set('globalt0', t0, xdict, t0)
2820
2821
2822 setup_pars = {'ics': xdict, 'algparams': {}}
2823
2824
2825
2826 try:
2827 if MI.get('algparams', xdict, t0)['init_step'] > t1-t0:
2828 if self.verboselevel > 0:
2829 print "Warning: time step too large for remaining time"\
2830 + " interval. Temporarily reducing time step to " \
2831 + "1/10th of its previous value"
2832 setup_pars['algparams'].update({'init_step': (t1-t0)/10})
2833 except (AttributeError, KeyError, ValueError):
2834
2835 pass
2836 else:
2837
2838 setup_pars['algparams']['verbose'] = self.verboselevel
2839 model.set(**setup_pars)
2840
2841
2842
2843 try:
2844 if modelNames[-1] == model.name:
2845
2846 model._set_for_hybrid_DS(True)
2847 except IndexError:
2848
2849 pass
2850 model.diagnostics.clearWarnings()
2851 model.diagnostics.clearErrors()
2852
2853
2854
2855
2856 try:
2857 traj = MI.get_test_traj(force=True)
2858 except PyDSTool_ValueError, e:
2859 print "\nError in Model:", model.name
2860 model.diagnostics.showWarnings()
2861 model.diagnostics.showErrors()
2862 print "Are the constituent generator absolute epsilon "
2863 print " tolerances too small? -- this abseps =", model._abseps
2864 raise
2865 except KeyboardInterrupt:
2866 raise
2867 except:
2868 print "\nError in Model:", model.name
2869 self.diagnostics.traceback = {}
2870 for k,v in model.diagnostics.traceback.iteritems():
2871 if isinstance(v, dict):
2872 dentry = {}
2873 for vk, vv in v.iteritems():
2874 dentry[traj._FScompatibleNamesInv(vk)] = vv
2875 else:
2876 dentry = v
2877 self.diagnostics.traceback[k] = dentry
2878 if self.diagnostics.traceback != {}:
2879 print "Traceback dictionary copied to model"
2880 raise
2881
2882 model._set_for_hybrid_DS(False)
2883 if traj is None:
2884 raise ValueError("Model %s failed to create a trajectory"%model.name)
2885 else:
2886
2887 model.renameTraj(MProject.ModelInterface._trajname,
2888 trajname+'_'+str(partition_num),
2889 force=force_overwrite)
2890
2891
2892
2893
2894 if not isparameterized(traj):
2895 raise ValueError("Model " + model.name + " produced a " \
2896 + " non-parameterized Trajectory")
2897 time_interval = traj.indepdomain
2898
2899
2900 ti_1 = time_interval[1]
2901 if self.verboselevel > 1:
2902 print "\nt0=", t0, "traj.globalt0=", traj.globalt0
2903 print "traj.indepdomain = ", traj.indepdomain.get()
2904 for vn, v in traj.variables.items():
2905 print vn + ":", v.trajirange, v.indepdomain.get()
2906 print "Last point of traj %s was: "%traj.name, traj(ti_1)
2907
2908
2909
2910
2911
2912
2913
2914 modelNames.append(model.name)
2915
2916
2917
2918 if model.diagnostics.hasErrors() and self.verboselevel > 0:
2919 for e in model.diagnostics.errors:
2920 print 'Model ' + model.name + \
2921 ' in trajectory segment had errors'
2922 model.diagnostics.showErrors()
2923
2924 end_reasons = ['time']
2925
2926
2927
2928
2929 if model.diagnostics.hasWarnings():
2930 if self.verboselevel >= 1:
2931 model.diagnostics.showWarnings()
2932 for w in model.diagnostics.warnings:
2933 if w[0] == Generator.W_TERMEVENT or \
2934 w[0] == Generator.W_TERMSTATEBD:
2935 if self.verboselevel > 1:
2936 print "Time:", ti_1+t0
2937 print 'Model ' + model.name + \
2938 ' had a terminal event. Details (in local ' + \
2939 'system time) ...\n ' + str(w[1])
2940
2941
2942
2943 if isinstance(w[1][1], list):
2944 if len(w[1][1]) > 1:
2945 if self.verboselevel > 0:
2946 print 'Warning: More than one terminal event found.'
2947 print ' Consider altering event pars.'
2948 end_reasons = [w[1][1][ri] \
2949 for ri in range(len(w[1][1]))]
2950 else:
2951 end_reasons = [w[1][1]]
2952
2953 self.diagnostics.update(model.diagnostics)
2954
2955
2956
2957
2958
2959
2960
2961 global_end_reasons = {}
2962 global_end_ixs = []
2963 for globalDS in globalConRules:
2964 if globalDS(MI):
2965 global_end_reasons[Inf] = \
2966 globalDS.conditions.collate_results('reasons',
2967 merge_lists=True)
2968 global_end_ixs.append(Inf)
2969 else:
2970
2971
2972
2973 try:
2974 final_ok_idx = globalDS.conditions._find_idx()
2975 except (AttributeError, RuntimeError):
2976 print "Trajectory creation failed..."
2977 globalDS.conditions.results.info()
2978 raise PyDSTool_ExistError("Global consistency checks failed for"
2979 " model interface %s for trajectory %s"%(str(globalDS),
2980 traj.name))
2981 else:
2982 if final_ok_idx is None:
2983 final_ok_idx = Inf
2984 global_end_ixs.append(final_ok_idx)
2985 global_end_reasons[final_ok_idx] = \
2986 globalDS.conditions.collate_results('reasons',
2987 merge_lists=True)
2988 if len(global_end_reasons) > 0:
2989 smallest_ok_idx = min(global_end_ixs)
2990
2991
2992 if isfinite(smallest_ok_idx):
2993
2994 traj.truncate_to_idx(smallest_ok_idx)
2995
2996 new_reason = global_end_reasons[smallest_ok_idx]
2997 if new_reason != []:
2998 end_reasons = new_reason
2999 if self.verboselevel > 1:
3000 print "End reason in %s was %s"%(self.name,str(end_reasons))
3001
3002
3003
3004
3005
3006
3007 partition_num += 1
3008 if time_interval.atEndPoint(t1_global-t0, 'hi'):
3009
3010 notDone = False
3011 if ti_1 > t1_global-t0:
3012
3013 print "Warning: Segment time interval exceeds prescribed limits:"
3014 print " ... ", ti_1, " > ", t1_global-t0
3015 notDone = False
3016
3017
3018
3019 if ti_1 < t1-t0:
3020 if end_reasons[0] not in swRules:
3021
3022
3023 print 'Trajectory calculation for Model `'+model.name+'` ' \
3024 +'terminated without a DS specified in the ' \
3025 +'switching rules from which to continue.'
3026 print 'Perhaps a variable went out of bounds:'
3027 model.diagnostics.showWarnings()
3028 print 'Last reasons for stopping were: ', end_reasons
3029 print "Trajectory calculated: "
3030 dt_sample = (ti_1-time_interval[0])/10.
3031 print traj.sample(dt=dt_sample)
3032 notDone = False
3033 raise PyDSTool_ValueError("Premature termination of "
3034 "trajectory calculation")
3035 elif swRules[end_reasons[0]][0] == 'terminate':
3036 if self.verboselevel > 0:
3037 print 'Trajectory calculation for Model `'+model.name+'` ' \
3038 +'terminated.'
3039 notDone = False
3040 t1 = ti_1+t0
3041 else:
3042 t1 = ti_1+t0
3043
3044 self._prepareICs(xdict, traj, MI, t0, ti_1)
3045
3046 trajseq.append(traj)
3047
3048 t0 = t1
3049 try:
3050 epochEvents.append(traj.getEvents())
3051 except AttributeError:
3052
3053 epochEvents.append({})
3054 MI_prev = MI
3055
3056
3057
3058
3059 self._addTraj(trajname, trajseq,
3060 epochEvents, modelNames, force_overwrite)
3061
3062 self.trajectory_defining_args[trajname] = args(ics=copy.copy(self.icdict),
3063 pars=copy.copy(self.pars),
3064 tdata=copy.copy(tdata))
3065
3067 """Validate Model's modelInfo attribute.
3068 Now that a HybridModel is invariably built using ModelConstructor,
3069 which already has validation built in, and also has more complex ways
3070 of construction, this method is defunct."""
3071 pass
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110
3111
3112
3113
3114
3115
3116
3117
3118
3119
3120
3121
3122
3123
3124
3125
3126
3127
3128
3129
3130
3131
3132
3134 pars_found = {}
3135 for infodict in self.modelInfo.values():
3136 try:
3137 pars_recognized = intersect(pars_found.keys(),
3138 infodict['dsi'].pars)
3139 except AttributeError:
3140
3141 continue
3142 pars_unknown = remain(infodict['dsi'].pars, pars_found.keys())
3143 for parname in pars_recognized:
3144 if infodict['dsi'].pars[parname] != pars_found[parname]:
3145 print "Inconsistent parameter values between model objects inside model"
3146 print "probably means that they were changed during trajectory computation"
3147 print "at discrete events: you should reset the pars when calling compute()"
3148 raise ValueError("Inconsistency between parameter values"
3149 " of same name ('%s') in different constiutent Generator"%parname +\
3150 " objects: reset during call to compute")
3151 new_vals = dict(zip(pars_unknown, [infodict['dsi'].pars[p] \
3152 for p in pars_unknown]))
3153 pars_found.update(new_vals)
3154
3155
3157 if verbose > 0:
3158 outputStr = 'Hybrid Model '+self.name+" containing components:"
3159 outputStr += "Observable variables: " + ",".join(self.obsvars)
3160 outputStr += "Internal variables: " + ",".join(self.intvars)
3161 outputStr += "Auxiliary variables: " + ",".join(self.auxvars)
3162 for name, infodict in self.modelInfo.iteritems():
3163 outputStr += "\n--- Sub-model: "+name
3164 outputStr += "\n "+infodict['dsi'].model._infostr(verbose-1)
3165 else:
3166 outputStr = 'Hybrid Model '+self.name
3167 return outputStr
3168
3169
3170
3171
3172
3174 """Return auxiliary variable values evaluated at t, icdict,
3175 pardict.
3176 """
3177
3178 icdict_local = copy.copy(icdict)
3179
3180
3181 old_icdict = copy.copy(dsi.model.icdict)
3182
3183
3184 dsi.model.icdict.update(icdict_local)
3185 icdict_local['t'] = t
3186 fspec = dsi.get('funcspec', icdict, t)
3187 if fspec.auxvars == []:
3188 auxdict = {}
3189 else:
3190 try:
3191 auxdict = dict(zip(fspec.auxvars,
3192 apply(dsi.get('AuxVars', icdict, t),
3193 (t, icdict_local, pardict)) ))
3194 except (ValueError, TypeError):
3195
3196
3197
3198 auxdict = {}
3199
3200 dsi.model.icdict = old_icdict
3201 return auxdict
3202
3203
3204 -def findTrajInitiator(modelInfo, t, vardict, pardict, intvars,
3205 verboselevel=0):
3206 """Find eligible Model to begin computation of trajectory.
3207 Cannot depend on any internal variables (e.g. those not common to
3208 all sub-models).
3209 """
3210 eligibleMI = []
3211 if len(modelInfo) == 1:
3212 return modelInfo.values()[0]
3213 outcome = {}
3214 for infodict in modelInfo.values():
3215 MI = infodict['dsi']
3216 MI_vars = MI.query('vars')
3217 xdict = filteredDict(vardict, MI_vars)
3218 dxdt_zeros = dict(zip(xdict.keys(),[0.]*len(xdict)))
3219 domtests = infodict['domainTests']
3220 globalConRules = infodict['globalConRules']
3221 model = MI.model
3222 outcome[model.name] = {}
3223 t_test = MI.get('indepvariable', xdict, t).indepdomain.contains(t) \
3224 is not notcontained
3225 fs = MI.get('funcspec', xdict, t)
3226 allvars = fs.vars + fs.auxvars
3227
3228 dsinps = MI.get('inputs', xdict, t)
3229 if dsinps != {}:
3230 try:
3231 input_ics = [dsinps[in_name](t) for in_name in dsinps.keys()]
3232 except KeyboardInterrupt:
3233 raise
3234 except:
3235 raise RuntimeError("Error in evaluating inputs for generator at"
3236 " initial condition")
3237 pardict.update(dict(zip(dsinps.keys(), input_ics)))
3238 icdict = copy.copy(xdict)
3239
3240
3241 icdict.update(filteredDict(getAuxVars(MI, t, icdict, pardict),
3242 icdict.keys(), neg=True))
3243
3244
3245
3246
3247
3248 for xname in allvars:
3249 if xname not in xdict or xname in intersect(intvars, xdict):
3250 if xname in xdict and xname in intvars:
3251
3252 continue
3253 if xname in fs.auxvars and xname in icdict and isfinite(icdict[xname]):
3254 continue
3255 try:
3256 if isfinite(model.icdict[xname]):
3257 icdict[xname] = model.icdict[xname]
3258 except KeyError:
3259
3260
3261 pass
3262
3263
3264 model.set(ics=icdict)
3265 x_test = True
3266 try:
3267
3268 dxdt = dict(MI.Rhs(t, icdict, pardict))
3269 except AttributeError:
3270
3271 dxdt = dxdt_zeros
3272 else:
3273 for vname, val in dxdt.items():
3274
3275
3276
3277 if not isfinite(val):
3278 dxdt[vname] = 1e308
3279 for aname in fs.auxvars:
3280
3281 dxdt[aname] = 0
3282
3283
3284
3285 for xname, x in icdict.iteritems():
3286
3287 if xname == 't' or xname in intvars or xname not in domtests:
3288
3289
3290
3291 continue
3292 xtraj = numeric_to_traj([[x], [dxdt[xname]]], 'test', [xname, 'D_'+xname], t)
3293
3294
3295 newtest = domtests[xname](xtraj)
3296 if verboselevel >=2:
3297 print "\nstate dom test for '%s' @ value %f:"%(xname, icdict[xname])
3298 print "depdomain is ", MI.get('variables', xdict, t)[xname].depdomain.get()
3299 print "-> test result is ", newtest
3300 x_test = x_test and newtest
3301 g_test = True
3302 xdict = filteredDict(icdict, intvars, neg=True)
3303 MI.test_traj = numeric_to_traj(array([xdict.values()]).T, 'ic_trajpt',
3304 xdict.keys(), t)
3305
3306
3307
3308
3309 xdict, t, I = MI._get_initiator_cache(xdict, t)
3310 while True:
3311
3312 if isinstance(I, MProject.GeneratorInterface):
3313 I.model.diagnostics.clearWarnings()
3314 break
3315 else:
3316 xdict, t, I = I._get_initiator_cache(xdict, t)
3317
3318
3319
3320 for gc in globalConRules:
3321 try:
3322 globtest = gc(MI)
3323 except KeyboardInterrupt:
3324 raise
3325 except:
3326 globtest = False
3327 else:
3328 if verboselevel >= 2:
3329 print "\nglobal con test for '%s' @ value %f:"%(str(gc), icdict[xname])
3330 print "-> test result is ", globtest
3331 g_test = g_test and globtest
3332 if verboselevel >= 1:
3333 print "\nModel '%s' tests..."%model.name
3334 print " ...for initial time: " + str(t_test)
3335 print " ...for initial state: " + str(x_test)
3336 print " ...for global conditions: " + str(g_test)
3337 outcome[model.name]['time dom. test'] = t_test
3338 outcome[model.name]['state dom. test'] = x_test
3339 outcome[model.name]['global con. test'] = g_test
3340 if t_test and x_test and g_test:
3341
3342 eligibleMI.append(infodict)
3343 if eligibleMI == []:
3344 info(outcome, "\nOutcomes for eligibility tests, by model")
3345 raise PyDSTool_ValueError('No eligible models from'
3346 ' which to begin trajectory computation')
3347 if len(eligibleMI) > 1:
3348 info(outcome, "\nOutcomes for eligibility tests, by model")
3349 raise PyDSTool_ValueError('Too many eligible models'
3350 ' to start trajectory computation')
3351
3352 return eligibleMI[0]
3353