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

Source Code for Module PyDSTool.Model

   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  ## PyDSTool imports
 
  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  ## Other imports
 
  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  
 
68 -class boundary_containment(MProject.qt_feature_leaf):
69 # not implemented using metrics because the metrics are trivial 70 # and cause a lot of overhead for this often-evaluated feature
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 # tolerance for small rounding errors 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 # test all coords at once 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
103 -class boundary_containment_by_event(boundary_containment):
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 # assume that the supplied model will correspond to the source of 111 # trajectories in evaluate method 112 try: 113 self.pars.model 114 except AttributeError: 115 raise ValueError("Missing model associated with event")
116 # have not verified that event is present in model 117
118 - def evaluate(self, traj):
119 # verify whether event exists and was flagged in associated model 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 # no events of this kind were found, so passed feature eval test 131 # dereferencing None (unsubscriptable object) 132 if self.pars.abseps > 0: 133 # would like to re-evaluate event at its threshold+abseps, but 134 # leave for now 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 # event found 145 if self.pars.abseps > 0: 146 # would like to re-evaluate event at its threshold+abseps, but 147 # leave for now 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] # only use first event (in case not Terminal) 151 self.results.uncertain = False 152 else: 153 self.results.output = evpt[0] # only use first event (in case not Terminal) 154 self.results.uncertain = False 155 return False
156
157 - def _find_idx(self):
158 """Helper function for finding index in trajectory meshpoints 159 at which containment first failed.""" 160 if self.results.satisfied: 161 # Trajectory satisfied constraint! 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 # especially for discrete domains 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 # if p values are *outside* thresh by up to abseps amount 185 # then flag this as 'uncertain' for use by domain_test class's 186 # transversality testing 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 # Trajectory satisfied constraint! 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 # starts with 0 already 207 adjusted_res = list(res != 0) 208 # find first index at which value is non-zero 209 # should never raise ValueError because this method is 210 # only run if there was a sign change found 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 # tolerance for small rounding errors 223 self.pars.abseps 224 except AttributeError: 225 if isinstance(self.pars.interval, Interval): 226 # Interval type passed, inherit its abseps 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 # assume a singleton numeric type passed 232 self.pars.interval = [self.pars.interval, self.pars.interval] 233 elif len(self.pars.interval)==1: 234 # singleton passed, so copy the value for both 235 # "endpoints" so that xlo_bc etc. below will work 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 # multiply interior directions by the integer value of not self.isdiscrete 249 # in order to set them to zero when the "interval" is actually a singleton 250 # value for a discrete domain. That fixes the boundary containment evaluation 251 # code which compares the sign of coord differences with that interior 252 # direction value. 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 # Arg is a traj! 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 # accept uncertain case for discrete domain 295 xlo_test = True 296 else: 297 # check transversality at critical (boundary) value of domain 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 # accept uncertain case for discrete domain 306 xhi_test = True 307 else: 308 # check transversality at critical (boundary) value of domain 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 # Trajectory satified domain conditions! 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 # dxdt transversality test was not run so ignore 328 continue 329 if sf.results.satisfied: 330 continue 331 if self.pars.verbose_level > 0: 332 print res 333 # Find first index at which value is non-zero. 334 # Will not raise ValueError because test satisfaction 335 # was already checked, so must have a zero crossing 336 idx = res.index(False) 337 if idx < lowest_idx: 338 lowest_idx = idx 339 return lowest_idx
340 341 342 # ------------------- 343 344
345 -class Model(object):
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 # query keys for 'query' method 354 _querykeys = ['pars', 'parameters', 'events', 'submodels', 355 'ics', 'initialconditions', 'vars', 'variables', 356 'auxvariables', 'auxvars', 'vardomains', 'pardomains', 357 'abseps'] 358 # valid keys for 'set' method 359 _setkeys = ['pars', 'algparams', 'checklevel', 'abseps', 360 'ics', 'inputs', 'tdata', 'restrictDSlist', 361 'globalt0', 'verboselevel', 'inputs_t0'] 362
363 - def __init__(self, legit, *a, **kw):
364 # legit is a way to ensure that instances of this abstract class 365 # are not created directly 366 if not legit==True: 367 # use explicit comparison to True otherwise kw argument will 368 # eval to True, which is not what we want 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 # modelInfo is a dict mapping model names --> a dict of: 382 # 'dsi': GeneratorInterface or ModelInterface object 383 # 'swRules': dict of switching rules (transitions between 384 # trajectory segments) 385 # 'globalConRules': list of global consistency DS names 386 # 'domainTests': dictionary of variable name -> domain 387 self.modelInfo = kw['modelInfo'] 388 except KeyError: 389 raise KeyError('Necessary keys missing from argument') 390 _foundKeys = len(self._needKeys) 391 392 # by default, 'observables' are all variables that are common 393 # to the Model/Generator objects in self.modelInfo, and the 394 # 'internals' are those remaining generate obsvars, intvars... 395 # aux vars are those aux vars common to *all* sub-models. 396 # create self.obsvars, self.intvars, self.auxvars 397 self.defaultVars() 398 399 # trajectories is a dict of trajectory segments (in sequence) 400 # ... can only add one at a time! 401 self.trajectories = {} 402 self.trajectory_defining_args = {} 403 # registry of generators or models (depending on sub-class) 404 # Using registry provides a shortcut for accessing a sub-model regardless 405 # of whether it's a Generator or a Model class 406 self.registry = {} 407 for name, infodict in self.modelInfo.iteritems(): 408 # set super model tag of ds object (which is either a 409 # ModelInterface or Generator) 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 # set initial conditions if specified already (if not, they must 418 # be specified before or during when compute() is called) 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 # impose this absolute epsilon (small scale) on all components 437 self._abseps = kw['abseps'] 438 _foundKeys += 1 439 else: 440 # All components use their defaults 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 # If not already created, a True result for 463 # self.haveJacobian() means that the Model will have a 464 # Jacobian method made during compute(), which 465 # references the appropriate _auxfn_Jac function in the Generator 466 # objects. 467 # 468 # Use this dict to record if any external input t0 time shift 469 # values are linked to a parameter value, for automatic 470 # updating before trajectory computation. (Internal use) 471 self._inputt0_par_links = {}
472 473
474 - def __len__(self):
475 """Return number of sub-models""" 476 return len(self.registry)
477
478 - def sub_models(self):
479 """Return a list of all sub-model instances (model interfaces or generators)""" 480 return self.registry.values()
481
482 - def _makeDefaultVarNames(self):
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 # first ds to have auxvars, so just add them all 494 auxvars.extend(auxvarnames) 495 else: 496 auxvars = intersect(auxvars, auxvarnames) 497 if obsvars == []: 498 # first ds, so add them all 499 obsvars = varnames 500 else: 501 obsvars = intersect(obsvars, varnames) 502 intvars = remain(all_known_varnames, obsvars) 503 return (obsvars, intvars, auxvars)
504
505 - def _generateParamInfo(self):
506 """Record parameter info locally, for future queries. 507 Internal use only. 508 """ 509 # use query method in case model in registry is a wrapped Generator 510 # that uses _ versions of hierarchical names that are used natively 511 # here. 512 self.pars = {} 513 for model in self.registry.values(): 514 try: 515 self.pars.update(model.query('pars')) 516 except AttributeError: 517 # no pars present 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
550 - def showSpec(self):
551 for ds in self.registry.values(): 552 ds.showSpec()
553
554 - def showAuxSpec(self):
555 for ds in self.registry.values(): 556 ds.showAuxSpec()
557
558 - def showAuxFnSpec(self):
559 for ds in self.registry.values(): 560 ds.showAuxFnSpec()
561
562 - def showEventSpec(self):
563 for ds in self.registry.values(): 564 ds.showEventSpec()
565
566 - def current_defining_args(self):
567 return args(pars=self.pars, ics=self.icdict, 568 tdata=self.tdata)
569
570 - def has_exact_traj(self, trajname, info):
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 # not even a trajectory of this name 579 return False
580
581 - def _prepareCompute(self, trajname, **kw):
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 # print "tdata in kw of %s (%s):"%(self.name, type(self)), tdata 598 else: 599 tdata = self.tdata 600 # print "tdata from self of %s (%s):"%(self.name, type(self)), tdata 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 # ds is a ModelInterface, not a Generator 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 # accumulate domains from each sub-model for regular variables 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 # singleton 669 vdom_lo = vdom.get() 670 vdom_hi = vdom_lo 671 else: 672 # range 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 # have to keep abseps the tightest 681 # of any of the instances for safety 682 result[vname]._abseps = vdom._abseps 683 result[vname] = vdom 684 # remaining vars are promoted aux vars 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 # accumulate domains from each sub-model for regular variables 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 # singleton 698 pdom_lo = pdom.get() 699 pdom_hi = pdom_lo 700 else: 701 # range 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 # have to keep abseps the tightest 710 # of any of the instances for safety 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
719 - def getEventMappings(self, dsName):
720 try: 721 return self.modelInfo[dsName]['swRules'] 722 except KeyError: 723 raise NameError("Sub-model %s not found in model"%dsName)
724
725 - def setPars(self, p, val):
726 # process multirefs first, then hierarchical names 727 # calls itself recursively to resolve all names 728 if isMultiRef(p): 729 # e.g to change a group of numerically indexed parameter names 730 # such as p0 - p9 or comp1.p0.g - comp1.p9.g 731 # extract numeric range of pars 732 # [ and ] are guaranteed to be present, from isMultiRef() 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 # no more of p after this multireference 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 # call setPars for each resolved name (these may include further 751 # multi references or hierarchical names 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 # find out: is root of p a valid 'type' in model spec? 759 # find all occurrences of last p 760 allFoundNames = [] 761 for mspecinfo in self._mspecdict.values(): 762 foundNames = searchModelSpec(mspecinfo['modelspec'], p) 763 # don't add duplicates 764 allFoundNames.extend(remain(foundNames,allFoundNames)) 765 # if allFoundNames == [] then either the hierarchical name was 766 # a specific reference, or the type matching is invalid. 767 # either way, we can just call set(p) and let that resolve the issue 768 # (note that all multi-refs will have been dealt with by this point) 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
777 - def setICs(self, p, val):
778 # process multirefs first, then hierarchical names 779 # calls itself recursively to resolve all names 780 if isMultiRef(p): 781 # e.g to change a group of numerically indexed parameter names 782 # such as p0 - p9 or comp1.p0.g - comp1.p9.g 783 # extract numeric range of pars 784 # [ and ] are guaranteed to be present, from isMultiRef() 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 # no more of p after this multireference 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 # call setICs for each resolved name (these may include further 803 # multi references or hierarchical names 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 # find out: is root of p a valid 'type' in model spec? 811 # find all occurrences of last p 812 allFoundNames = [] 813 for mspecinfo in self._mspecdict.values(): 814 foundNames = searchModelSpec(mspecinfo['modelspec'], p) 815 # don't add duplicates 816 allFoundNames.extend(remain(foundNames,allFoundNames)) 817 # if allFoundNames == [] then either the hierarchical name was 818 # a specific reference, or the type matching is invalid. 819 # either way, we can just call set(p) and let that resolve the issue 820 # (note that all multi-refs will have been dealt with by this point) 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 # Handle initial conditions here, because compute will pass 848 # the values on to the appropriate sub-models when they are called. 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 # Propagate values to sub-models. 858 # If any values are strings they must refer to a parameter, so 859 # here we evaluate them. We keep a record of these for later 860 # automated updating. 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 # no longer parameter-linked 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 # generator doesn't involve inputs 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 # propagate to sub-models 889 try: 890 model.set(verboselevel=self.verboselevel) 891 except KeyError: 892 # generator doesn't support verboselevel 893 pass 894 if restrictDSlist == []: 895 restrictDSlist = self.registry.keys() 896 # For the remaining keys, must propagate parameter changes to all 897 # sub-models throughout modelInfo structure. 898 # 899 # Changed by WES 10FEB06 to handle problem of 'new' pars being 900 # added if the parameter names do not exist in any generators 901 dsis = self.modelInfo.values() 902 numDSs = len(dsis) 903 # loop over keywords 904 for key, value in filteredDict(kw, ['ics', 'tdata', 905 'inputs_t0', 'restrictDSlist', 'globalt0'], 906 neg=True).iteritems(): 907 # keep track of the number of errors on this keyword 908 if isinstance(value, dict): 909 # keep track of entry errors for this key 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 # try setting pars in each sub-model 919 for infodict in dsis: 920 callparsDict = {} 921 # select out the ones relevant to this sub-model 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 # only apply these keys to the restricted list 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 # key not valid for this type of ds 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 # key not valid for this type of ds 952 pass 953 # Check that none of the entries in the dictionary caused errors 954 # in each sub-model 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 # can't think of other ways for this error to crop up 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
978 - def __getitem__(self, trajname):
979 try: 980 return self.trajectories[trajname] 981 except KeyError: 982 raise ValueError('No such trajectory.')
983
984 - def __delitem__(self, trajname):
985 self._delTraj(trajname)
986
987 - def _delTraj(self, trajname):
988 """Delete a named trajectory from the database.""" 989 try: 990 traj = self.trajectories[trajname] 991 except KeyError: 992 # a trajectory piece may have been created without 993 # the top-level trajectory ever being completed 994 # (e.g. after an unxpected error or ^C interruption) 995 ##raise ValueError('No such trajectory.') 996 l = len(trajname) 997 for m in self.registry.values(): 998 # delete all matching pieces (of form trajname + '_' + <digits>) 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 # propagate deletions down through registry 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
1025 - def getEndPoint(self, trajname, end=1):
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 # auxiliary var didn't need calling 1041 pass 1042 except AttributeError: 1043 # non-point based output attributes of a Variable need 1044 # to be called ... 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
1057 - def getEndTime(self, trajname, end=1):
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
1069 - def _validateVarNames(self, names):
1070 """Check types and uniqueness of variable names.""" 1071 namelist = [] # records those seen so far 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
1079 - def forceObsVars(self, varnames):
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 # then there are names given that are not known as 1085 # obs, int, or aux 1086 raise ValueError("Unknown variable names: "+str(r)) 1087 for v in remain(varnames, self.obsvars): 1088 # only include names that are not already observables 1089 self.obsvars.append(v) 1090 # remove any vars that are now observables 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
1101 - def forceIntVars(self, varnames):
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 # then there are names given that are not known as 1107 # obs, int, or aux 1108 raise ValueError("Unknown variable names: "+str(r)) 1109 for v in remain(varnames, self.intvars): 1110 # only include names that are not already internals 1111 self.intvars.append(v) 1112 # remove any vars that are now internals 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
1123 - def defaultVars(self):
1124 """(Re)set to default observable and internal variable names.""" 1125 obsvars, intvars, auxvars = self._makeDefaultVarNames() 1126 self._validateVarNames(obsvars + intvars + auxvars) 1127 # OK to store these permanently after varname validation 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):
1140 print self._infostr(verboselevel)
1141 1142
1143 - def __repr__(self):
1144 return self._infostr(verbose=0)
1145 1146 __str__ = __repr__ 1147 1148
1149 - def __copy__(self):
1150 pickledself = pickle.dumps(self) 1151 return pickle.loads(pickledself)
1152 1153
1154 - def __deepcopy__(self, memo=None, _nil=[]):
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
1175 - def getTrajModelName(self, trajname, t=None):
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 # return list of Generators for associated hybrid trajectory 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
1194 - def getTrajEventTimes(self, trajname, events=None):
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 #raise ValueError('No such trajectory name') 1207
1208 - def getTrajEvents(self, trajname, events=None):
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
1221 - def getTrajEventStruct(self, trajname):
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
1230 - def getTrajTimeInterval(self, trajname):
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
1238 - def getTrajTimePartitions(self, trajname):
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
1246 - def getDSAlgPars(self, target, par, idx=None):
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
1273 - def setDSAlgPars(self, target, par, val):
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 # If target is a list, make sure that the input list is of the 1296 # appropriate length, otherwise warn and skip 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 # Set every member of the list to that value 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 # return more information for a single sub-model 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
1363 - def showDSEventInfo(self, target, verbosity=1, ics=None, t=0):
1364 # call to eventstruct prints info to stdout 1365 if ics is None: 1366 ics = self.icdict 1367 estruct = self.modelInfo[target]['dsi'].get('eventstruct', ics, t) 1368 estruct(verbosity)
1369
1370 - def getDSEventTerm(self, target, flagVal=True, ics=None, t=0):
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
1408 - def getDSEventActive(self, target, flagVal=True, ics=None, t=0):
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
1426 - def setDSEventActive(self, target, eventTarget, flagVal, ics=None, t=0):
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
1468 - def setDSEventPrecise(self, target, eventTarget, flagVal, ics=None, t=0):
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
1487 - def setDSEventDelay(self, target, eventTarget, val, ics=None, t=0):
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
1506 - def setDSEventInterval(self, target, eventTarget, val, ics=None, t=0):
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
1563 - def setDSEventStartTime(self, target, eventTarget, val, ics=None, t=0):
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
1582 - def setDSEventBisect(self, target, eventTarget, val, ics=None, t=0):
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
1601 - def resetEventTimes(self):
1602 """Internal method""" 1603 for ds in self.registry.values(): 1604 ds.resetEventTimes()
1605
1606 - def _set_for_hybrid_DS(self, state):
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
1614 - def searchForNames(self, template):
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 # HACK to force compatibility of mspecinfo being a dict/args vs. a GDescriptor 1623 if isinstance(mspecinfo, dict): 1624 foundNames = searchModelSpec(mspecinfo['modelspec'], template) 1625 else: 1626 # args or GDescriptor 1627 foundNames = searchModelSpec(mspecinfo.modelspec, template) 1628 result[modelName] = foundNames 1629 return result
1630 1631
1632 - def searchForVars(self, template):
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 # all sub-models must define the same variables, so need only 1639 # to look at one 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 # for non-hybrid models 1647 fspec = self.registry[a_ds_name].registry.values()[0].get('funcspec') 1648 # funcspec won't have its internal names converted to 1649 # hierarchical form, so do it here on the attributes we 1650 # need 1651 return intersect(self.registry[a_ds_name]._FScompatibleNamesInv(fspec.vars + fspec.auxvars), 1652 foundNames)
1653 1654 1655 1656
1657 -class NonHybridModel(Model):
1658 - def __init__(self, *a, **kw):
1659 Model.__init__(self, True, *a, **kw) 1660 # collect parameter info from all modelInfo objects 1661 self._generateParamInfo() 1662 self._validateRegistry(self.obsvars, self.intvars)
1663
1664 - def _findTrajInitiator(self, end_reasons, partition_num, t0, xdict, 1665 gi=None, swRules=None):
1666 # return GeneratorInterface object and include return 1667 # elements for compatibility with HybridModel._findTrajInitiator: 1668 # swRules, globalConRules, nextModelName, reused (True b/c always same), 1669 # epochStateMaps, notDone 1670 infodict = self.modelInfo.values()[0] 1671 return infodict['dsi'], infodict['swRules'], \ 1672 infodict['globalConRules'], infodict['dsi'].model.name, \ 1673 True, None, True
1674
1675 - def cleanupMemory(self):
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
1681 - def haveJacobian(self):
1682 """Returns True iff all objects in modelInfo have 1683 defined Jacobians.""" 1684 return self.registry.values()[0].haveJacobian()
1685
1686 - def haveJacobian_pars(self):
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 # get Generator as 'ds' 1703 ds = self.registry.values()[0] 1704 fscm = ds._FScompatibleNames 1705 fscmInv = ds._FScompatibleNamesInv 1706 vars = ds.get('funcspec').vars # FS compatible 1707 x_fs = filteredDict(fscm(xdict), vars) 1708 # in case ds i.c.'s are different or not set yet (NaN's) 1709 # for purposes of auxiliary function 'initcond' calls during Rhs 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 # some may not be in the same order, so must inefficiently ensure 1717 # re-ordered according to alphabetical order of non-FS compatible names 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 # FS compatible 1745 x_fs = filteredDict(fscm(xdict), vars) 1746 # in case ds i.c.'s are different or not set yet (NaN's) 1747 # for purposes of auxiliary function 'initcond' calls during Rhs 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 # in case ds i.c.'s are different or not set yet (NaN's) 1782 # for purposes of auxiliary function 'initcond' calls during Rhs 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 # FS compatible 1816 x_fs = filteredDict(fscm(xdict), vars) 1817 # in case ds i.c.'s are different or not set yet (NaN's) 1818 # for purposes of auxiliary function 'initcond' calls during Rhs 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 # in case ds i.c.'s are different or not set yet (NaN's) 1851 # for purposes of auxiliary function 'initcond' calls during Rhs 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 # just get i.c. from the single generator, 1886 # making sure it's non-empty 1887 self.icdict = self.registry.values()[0].get('initialconditions') 1888 if self.icdict == {}: 1889 # gen's i.c.s were empty too 1890 raise PyDSTool_ExistError("No initial conditions specified") 1891 xdict = {} 1892 for xname, value in self.icdict.iteritems(): 1893 #if xname not in self.obsvars+self.intvars: 1894 # raise ValueError("Invalid variable name in initial " 1895 # "conditions: " + xname) 1896 xdict[xname] = ensurefloat(value) 1897 # clean up self.icdict 1898 self.icdict = xdict.copy() 1899 1900 ## compute trajectory segment until an event or t1 is reached 1901 gen = self.registry.values()[0] # sole entry, a Generator 1902 # print "N-h model before fiddling: t0_global =", t0_global, "t1_global =", t1_global 1903 t0_global += gen.globalt0 1904 t1_global += gen.globalt0 1905 t0 = t0_global # initial value 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 # For sub-models that support setPars() and 1917 # compute() methods, set the sub-model's global time 1918 # reference (in case it needs it), and verify that t=0 is valid 1919 # in the sub-model's time domain 1920 # print "\nNonhybrid model %s 'compute'"%self.name 1921 # print " generator global t0 =", gen.globalt0 1922 # print " model t0_global=", t0_global, "t1_global=",t1_global 1923 # print " original model tdata = [%f,%f], gen tdata = [0,%f]"%(tdata[0],tdata[1],t1-t0), " and globalt0=",t0 1924 1925 # add remaining pars for system 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 # system does not support this integration parameter 1940 pass 1941 # if external inputs used, make sure any parameter-bound 1942 # global t0 offsets are updated in case epoch state map changed them 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 #### Compute a trajectory segment as a NonHybrid Trajectory 1953 #print "Entering generator %s at t0=%f, t1=%f"%(gen.name, t0, t1) 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 # update original copy of the generator with new event times during 1992 # this trajectory, for high level event detection to be able to check 1993 # eventinterval for terminal events between hybrid steps 1994 # ????? 1995 ##self.modelInfo[gen.name]['dsi'].eventstruct = copy.copy(gen.eventstruct) 1996 1997 #### Post-process trajectory segment: 1998 # look at warnings etc. for any terminating events that occurred 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 # NB. end reason will always be 'time' for ExplicitFnGen 2006 # Default time interval after Generator completed, 2007 # in case it was truncated by a terminal event 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 # w[1] always has t as first entry, and 2022 # either a state variable name 2023 # or a list of terminal event names as second. 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 #print "End reason: ", end_reasons[0] 2034 self.diagnostics.update(gen.diagnostics) 2035 2036 # Apply global consistency conditions, and truncate 2037 # trajectory to the last OK index in the independent variable 2038 # if any fail, for those conditions that do not apply to the 2039 # whole trajectory in an accept/reject fashion. If any fail 2040 # without specifying an index then the whole trajectory is being 2041 # rejected. 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 # if global consistency fails then the features *must* provide 2053 # the last OK position in the trajectory to truncate to, 2054 # otherwise the whole trajectory FAILs. 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 # this idx may be Inf if globalDS passed all points OK 2072 # - but it may have still changed (refined) the end reason 2073 if isfinite(smallest_ok_idx): 2074 # truncate in place 2075 traj.truncate_to_idx(smallest_ok_idx) 2076 # overwrite current end reason 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 # DEBUG 2081 # print "Nonhybrid traj info:", traj.indepdomain.get(), traj.globalt0, type(traj) 2082 # print "traj vars:", traj.variables.keys(), "\n" 2083 2084 #### Clean up 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 # this Generator has no eventstruct 2093 pass 2094 2095 ## New code in Trajectory.__init__ takes care of this now 2096 # Take a copy of the events as they were at the time that the 2097 # trajectory was computed, for future reference 2098 #traj.events = epochEvents 2099 #traj._createEventTimes() 2100 # TEMP 2101 #print "Nonhybridmodel - temp removal of eventstruct non-copy" 2102 #traj.modelEventStructs = gen.eventstruct #copy.deepcopy(gen.eventstruct) 2103 self.trajectories[trajname] = traj
2104 2105
2106 - def _validateRegistry(self, obsvars, intvars):
2107 """Validate Model's modelInfo attribute.""" 2108 # ensure that modelInfo is a single Generator object only 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
2116 - def _infostr(self, verbose=1):
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
2131 -class HybridModel(Model):
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 """
2142 - def __init__(self, *a, **kw):
2143 Model.__init__(self, True, *a, **kw) 2144 # collect parameter info from all modelInfo objects 2145 self._generateParamInfo() 2146 # Ensure all Generators provided to build trajectory share the same 2147 # observables, and that modelInfo switch rules mappings are OK 2148 self._validateRegistry(self.obsvars, self.intvars) 2149 # Check that all pars of the same name are equal 2150 self._validateParameters()
2151
2152 - def _prepareICs(self, xdict, traj, ds, t0, ti_1):
2153 # ds_old no longer needed 2154 # get new initial conditions for this epoch (partition) 2155 # DEBUG 2156 # print "\n\n_prepareICs for HybridModel %s and old traj %s"%(self.name,traj.name) 2157 # print "t0=%f, ti_1=%f"%(t0,ti_1) 2158 # print "old traj globalt0 =", traj.globalt0 2159 # print "t domain =", traj.indepdomain.get() 2160 # print dict(traj(traj.indepdomain[1], asGlobalTime=False)), "\n" 2161 for xname in self.allvars+self.auxvars: 2162 # update final condition of epoch to be used to generate 2163 # next initial condition when we re-enter the loop 2164 xname_compat = traj._FScompatibleNames(xname) 2165 try: 2166 xdict[xname] = traj.variables[xname_compat].output.datapoints[1][-1] 2167 except KeyError: 2168 # e.g. auxiliary var didn't need calling 2169 pass 2170 except AttributeError: 2171 # non-point based output attributes of a Variable need 2172 # to be 'called' ... 2173 try: 2174 xdict[xname] = traj.variables[xname_compat](ti_1) 2175 except RuntimeError: 2176 # traj is a hybrid traj 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 # the mapping might be the ID fn 2199 num_maps = len(epochStateMaps) 2200 # apply all maps to the initial condition 2201 # but check that the ordering doesn't matter (they are 2202 # commutative), otherwise tell the user we can't safely 2203 # proceed [FOR NOW, JUST DON'T ALLOW SIMULTANEOUS EVENTS] 2204 if num_maps > 1: 2205 # ensure that all the maps are the same for multiple 2206 # simultaneous events that point to the same generator 2207 # (the latter was already verified above) 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 # Now work out new inputs state 2214 dsinps = model_interface.get('inputs', xdict, t0) 2215 if dsinps != {}: 2216 # build dictionary of initial values for the 2217 # external inputs 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 # epochStateMap might update xdict_temp and pars_dict 2235 epochStateMap.evmapping(xdict_temp, pars_temp, 2236 extinputs_ic, eventstruct, t0) 2237 # use update method to change dicts *in place* 2238 pars.update(traj._FScompatibleNamesInv(pars_temp)) 2239 # update next MI's parameters if mapped by the event mapping 2240 # (the one corresponding to the original xdict) 2241 try: 2242 next_model_interface.set('pars', pars, xdict, t0) 2243 except (PyDSTool_AttributeError, PyDSTool_ValueError): 2244 # next model may not have these params, so just ignore 2245 pass 2246 xdict.update(traj._FScompatibleNamesInv(xdict_temp)) 2247 # No longer need the aux vars 2248 for xname in self.auxvars: 2249 try: 2250 del xdict[xname] 2251 except KeyError: 2252 # aux var value never put here 2253 pass
2254
2255 - def cleanupMemory(self):
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 # Initial values for these 2267 epochStateMaps = None 2268 notDone = True 2269 reused = False 2270 if end_reasons is None: 2271 # first time in the while loop, so set up xnames 2272 assert partition_num == 0, ("end_reasons was None on a " 2273 "non-initial epoch") 2274 xnames = self.icdict.keys() 2275 xnames.sort() 2276 # only for the initial Generator of a trajectory 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 # find new Generator using switching rules 2304 num_reasons = len(end_reasons) 2305 # num_reasons > 1 will occur most often 2306 # for discrete time systems 2307 if num_reasons > 1: 2308 # then all next generators must be the same! 2309 candidateNextModel = "" 2310 mapinfos = [] 2311 for ri in range(num_reasons): 2312 try: 2313 mapinfos.append(swRules[end_reasons[ri]]) 2314 except KeyError: 2315 # no rule for this termination reason 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 # No next generators found 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 # will not be needed 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 # terminates compute method 2346 elif nextModelName == mi.model.name: 2347 # Reuse the same model copy (take advantage of 2348 # not having to re-copy the model from modelInfo) 2349 # Currently we can't use the continue integration 2350 # option of ODE solvers because of the way that each 2351 # new hybrid partition is solved in "local time", i.e. 2352 # from t=0. So this would mess up the solvers` 2353 # 'continue' option. 2354 reused = True # not used in compute 2355 # swRules is untouched 2356 # return same globalConRules as had before 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 #genEvStructs = dict(zip(modelNames, 2377 # [copy.copy(self.modelInfo[gn]['dsi'].eventstruct) \ 2378 # for gn in modelNames])) 2379 # TEMP 2380 #print "Model: avoid copying eventstructs?" 2381 genEvStructs = dict(zip(modelNames, 2382 [self.modelInfo[gn]['dsi'].eventstruct \ 2383 for gn in modelNames])) 2384 # time_range is whole time range over sequence 2385 time_range = [None, None] 2386 indepvarname = '' 2387 genEventTimes = {} # keyed by event name 2388 genEvents = {} # keyed by event name 2389 # collate all FScompatible name maps from trajectories 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 ## else: 2398 ## if indepvarname != traj.indepvarname: 2399 ## print indepvarname, traj.indepvarname 2400 ## raise ValueError("Independent " 2401 ## "variable name mismatch in trajectory") 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 # in extreme cases with highest precision events, 2413 # may have an "equal" time being added 2414 if genEventTimes[evname][-1] - val['t'][0] > 0: 2415 # some other ValueError 2416 raise 2417 else: 2418 # don't add this point, it's identical 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 # first partition of sequence 2436 DS_t_numtype = traj.indepdomain.type 2437 if isinputcts(traj): 2438 DS_t_typestr = 'continuous' 2439 # continuous time must be float 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 # discrete time can be floats or ints 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 # Pick one observable variable ('varname') to test with -- they 2450 # must all have the same type. This variable must be present in 2451 # all trajectory segments. 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 ## assert DS_x_numtype == self.obsvars[varname], \ 2458 ## ('Mismatch between declared type of variable and ' 2459 ## 'that found in trajectory segment') 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 # DEBUG 2468 # temp 2469 # print "\n---- %s _addTraj case 1: "%self.name 2470 # print traj.indepdomain.get() 2471 # print time_range 2472 # print traj.globalt0 2473 time_partitions = [(traj.indepdomain, \ 2474 traj.globalt0, \ 2475 traj.checklevel)] # initial value 2476 else: 2477 # remove the following line to support trajectories that 2478 # only have partially-defined variable information, esp. 2479 # if part is a map (vars defined only at discrete times) 2480 if not compareNumTypes(DS_t_numtype, traj.indepdomain.type): 2481 raise TypeError('Mismatched time types for hybrid ' 2482 'trajectory sequence') 2483 # print "\n---- %s _addTraj case 2: "%self.name 2484 # print traj.indepdomain.get() 2485 # temp = copy.copy(time_range) 2486 # temp[1] += traj.indepdomain[1] 2487 # print temp 2488 # print traj.globalt0 2489 if not traj.indepdomain.atEndPoint(time_range[1] \ 2490 - traj.globalt0, 'lo'): 2491 # temp 2492 # print "\n***", time_range 2493 # print traj.indepdomain.get() 2494 # print traj.globalt0 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 # Full time interval of the trajectory 2504 time_interval = Interval(indepvarname, DS_t_numtype, time_range, 2505 abseps=self._abseps) 2506 # Add to trajectory dictionary, using hierarchical event names for 2507 # genEvents and genEventTimes (which are called directly by user methods 2508 # of Model) 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
2523 - def haveJacobian(self):
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
2531 - def haveJacobian_pars(self):
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 # get returns FS compatible names 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 # initially expect to compute over a global time interval, 2699 # [t0_global, t1_global], which is truncated if it extends beyond 2700 # DS's independent variable domain to give [t0, t1]. 2701 # DS tdata is then set to compute over relative time interval, 2702 # [0, t1-t0]. 2703 tdata, t0_global, t1_global, force_overwrite = \ 2704 Model._prepareCompute(self, trajname, **kw) 2705 # Set initial reason for an epoch to end to None to initiate 2706 # search for eligible Generator or Model for first epoch 2707 end_reasons = None 2708 # check initial condition specification in icdict 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 # clean up self.icdict 2719 self.icdict = xdict.copy() 2720 2721 # initial values 2722 notDone = True # finished computing trajectory segment? 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 # flag for re-use of a model from one hybrid segment to the next 2734 reused = False 2735 # reset persistent storage of event times 2736 self.resetEventTimes() 2737 2738 # From t0 (if non-autonomous system), icdict, and switching rules 2739 # (self-consistency validity conditions), determine which model 2740 # applies for the initial portion of the trajectory 2741 while notDone: 2742 # find appropriate model to compute trajectory segment 2743 # (reused flag indicates whether MI is the same as on previous loop 2744 # but is not used at the moment) 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 # convenient shorthand 2753 model = MI.model 2754 2755 if partition_num > 0: 2756 ## map previous ds state 2757 # to new ds state using epoch mapping 2758 # and current ds's pars and external inputs 2759 # (previous) traj is well defined from previous run 2760 pars_copy = copy.copy(self.pars) 2761 self._applyStateMap(epochStateMaps, MI_prev, MI, traj, xdict, t0) 2762 ## if changed last MI's event delay, then change it back. 2763 for event, delay in event_delay_record.items(): 2764 event.eventdelay = delay 2765 2766 ## temporarily set terminal events from previous in this MI to have 2767 # event_delay = event_interval because event finding code 2768 # does not know that an event just occurred in previous MI and 2769 # partition. 2770 event_delay_record = {} 2771 dummy1, dummy2, targetMI = MI._get_initiator_cache(xdict, t0) 2772 evstruct = targetMI.get('eventstruct', xdict, t0) 2773 # record change for later resetting 2774 for reason in end_reasons: 2775 # should only be one! 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 # way to force partition 0 to definitely update any 2785 # parameter-linked input t0's 2786 pars_copy = {} 2787 event_delay_record = {} 2788 2789 # if external inputs used, make sure any parameter-bound 2790 # global t0 offsets are updated in case epoch state map changed them 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 ## compute trajectory segment until an event or t1 is reached 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 # For Models that support setPars() and 2811 # compute() methods, set the global time 2812 # reference (in case it needs it), and verify that t=0 is valid 2813 # in the time domain 2814 ## next line redundant (and changed) by setup_pars['tdata'] below 2815 # print "\nHybrid model %s: set tdata=[0, %f] in %s (type %s)"%(self.name, t1-t0, MI.model.name, str((type(MI), type(MI.model)))) 2816 # print " H model t0, t1 =", t0, t1 2817 # print " and globalt0 = ", t0 2818 MI.set('tdata', [0, t1-t0], xdict, t0) 2819 MI.set('globalt0', t0, xdict, t0) 2820 2821 # add remaining pars for system 2822 setup_pars = {'ics': xdict, 'algparams': {}} 2823 #if self._abseps is not None: 2824 # setup_pars['abseps'] = self._abseps 2825 # setup_pars['algparams'].update({'abseps': self._abseps}) 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 # system does not support this integration parameter 2835 pass 2836 else: 2837 # we know there's compatibility to set verbose level 2838 setup_pars['algparams']['verbose'] = self.verboselevel 2839 model.set(**setup_pars) 2840 # ensure that if reusing same model as previous segment, that 2841 # any high level events are not reset: especially for VODE and map system 2842 # (!!! unsure why not needed for Dopri & Radau) 2843 try: 2844 if modelNames[-1] == model.name: 2845 # reusing same model again 2846 model._set_for_hybrid_DS(True) 2847 except IndexError: 2848 # first hybrid partition, so doesn't apply 2849 pass 2850 model.diagnostics.clearWarnings() 2851 model.diagnostics.clearErrors() 2852 # assert model._abseps == self._abseps 2853 2854 #### Compute a trajectory segment (as a Hybrid Trajectory) 2855 # print "DEBUG: Entering model %s at t0=%f, t1=%f"%(model.name, t0, t1) 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 # ensure this is reset to False 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 # rename from ModelInterface default name 2887 model.renameTraj(MProject.ModelInterface._trajname, 2888 trajname+'_'+str(partition_num), 2889 force=force_overwrite) 2890 # if self._abseps is not None: 2891 # assert traj.indepdomain._abseps == self._abseps 2892 # for v in traj.depdomain.values(): 2893 # assert v._abseps == self._abseps 2894 if not isparameterized(traj): 2895 raise ValueError("Model " + model.name + " produced a " \ 2896 + " non-parameterized Trajectory") 2897 time_interval = traj.indepdomain 2898 # ti_1 means time interval index 1 (i.e., the endpoint) 2899 # this is a local time 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 # update original copy of the model with new event times during 2908 # this trajectory, for high level event detection to be able to check 2909 # eventinterval for terminal events between hybrid steps 2910 # ????? 2911 ##MI.eventstruct = copy.copy(model.eventstruct) 2912 2913 # append record of the sequence of used generators 2914 modelNames.append(model.name) 2915 2916 #### Post-process trajectory segment: 2917 # look at warnings etc. for any terminating events that occurred 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 # default end reason is time (may be overwritten) 2924 end_reasons = ['time'] 2925 # Default time interval after traj computation completed, 2926 # in case it was truncated by a terminal event 2927 # DEBUG 2928 # print "Hybrid Traj %s t0, ti_1 ="%traj.name, traj.indepdomain[0], ti_1 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 # w[1] always has t as first entry, and 2941 # either a state variable name 2942 # or a list of terminal event names as second. 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 #print "End reason: ", end_reasons[0] 2953 self.diagnostics.update(model.diagnostics) 2954 2955 # Apply global consistency conditions, and truncate 2956 # trajectory to the last OK index in the independent variable 2957 # if any fail, for those conditions that do not apply to the 2958 # whole trajectory in an accept/reject fashion. If any fail 2959 # without specifying an index then the whole trajectory is 2960 # rejected. 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 # if global consistency fails then the features *must* provide 2971 # the last OK position in the trajectory to truncate to, 2972 # otherwise the whole trajectory FAILs. 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 # this idx may be Inf if globalDS passed all points OK 2991 # - but it may have still changed (refined) the end reason 2992 if isfinite(smallest_ok_idx): 2993 # truncate in place 2994 traj.truncate_to_idx(smallest_ok_idx) 2995 # overwrite current end reason only if it's non-empty 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 # DEBUG 3003 # print "\nTraj segment %i info:"%partition_num, traj.indepdomain.get(), traj.globalt0, type(traj) 3004 # print "\n" 3005 3006 #### Prepare for next partition & clean up 3007 partition_num += 1 3008 if time_interval.atEndPoint(t1_global-t0, 'hi'): 3009 # we reached end of desired computation interval 3010 notDone = False 3011 if ti_1 > t1_global-t0: 3012 # We computed too far for the time interval 3013 print "Warning: Segment time interval exceeds prescribed limits:" 3014 print " ... ", ti_1, " > ", t1_global-t0 3015 notDone = False 3016 # last Model's time == t1 should hold after successful traj computation 3017 # otherwise t1 needs to be updated according to last terminating 3018 # event time point 3019 if ti_1 < t1-t0: 3020 if end_reasons[0] not in swRules: 3021 # there's prescribed time left to compute but no 3022 # eligible DS to transfer control to 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 # update xdict with final point in traj 3044 self._prepareICs(xdict, traj, MI, t0, ti_1) 3045 # store trajectory segment in the sequence 3046 trajseq.append(traj) 3047 # before moving on, update t0 for next epoch / regime 3048 t0 = t1 3049 try: 3050 epochEvents.append(traj.getEvents()) 3051 except AttributeError: 3052 # no eventstruct present 3053 epochEvents.append({}) 3054 MI_prev = MI 3055 # end while loop 3056 3057 # Take a copy of the events as they were at the time that the 3058 # trajectory was computed, for future reference 3059 self._addTraj(trajname, trajseq, 3060 epochEvents, modelNames, force_overwrite) 3061 # Record defining data for this trajectory (especially useful for interfaces) 3062 self.trajectory_defining_args[trajname] = args(ics=copy.copy(self.icdict), 3063 pars=copy.copy(self.pars), 3064 tdata=copy.copy(tdata))
3065
3066 - def _validateRegistry(self, obsvars, intvars):
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 ## # ensure that modelInfo consists of more than one ModelInterface object 3074 ## # unless it's a single model punctuated by discrete events 3075 ## if len(self.modelInfo) == 1: 3076 ## infodict = self.modelInfo.values()[0] 3077 ## swmapping = infodict['swRules'] 3078 ## if all([outcome[0] == 'terminate' \ 3079 ## for (reason, outcome) in swmapping.iteritems()]): 3080 ## # then all point to 'terminate' and user needed a non-hybrid 3081 ## # model class, otherwise has discrete event mappings 3082 ## # that make the model technically "hybrid" 3083 ## raise AssertionError("Use a non-hybrid Model class") 3084 ## for infodict in self.modelInfo.values(): 3085 ## if not isinstance(infodict['dsi'], MProject.ModelInterface): 3086 ## raise TypeError("Must provide ModelInterface objects") 3087 ## allDSnames = self.modelInfo.keys() 3088 ### print "\n_validateRegistry: ", obsvars, intvars 3089 ## for modelName, infodict in self.modelInfo.iteritems(): 3090 ## # dynamical system for trajectory 3091 ## miface = infodict['dsi'] 3092 ## # dict of reason -> (ds_name, epmapping) pairs 3093 ## swmapping = infodict['swRules'] 3094 ## allvarnames = miface.model.query('variables') 3095 ### print "Model %s had vars "%modelName, allvarnames 3096 ## special_reasons = ['time'] + allvarnames 3097 ## assert miface.model.name == modelName, ('DS`s name does not ' 3098 ## 'correspond to the key in the modelInfo dict') 3099 ## assert miface.model.name in allDSnames, ('DS`s name not in list ' 3100 ## 'of all available DS names!') 3101 ## assert len(intersect(obsvars, allvarnames)) ==\ 3102 ## len(obsvars), 'DS `'+miface.model.name+'` did not contain ' \ 3103 ## + ' required observable variables' 3104 ## assert len(intersect(intvars, allvarnames)) ==\ 3105 ## len(intvars), 'DS '+miface.model.name+' did not contain ' \ 3106 ## + ' required internal variables' 3107 ## try: 3108 ## allEndReasonNames = miface.model.query('events').keys() \ 3109 ## + special_reasons 3110 ## except AttributeError: 3111 ## allEndReasonNames = special_reasons 3112 ## seenReasons = [] 3113 ## # can't assert this when non-event based reasons are used 3114 ## # (e.g. for DSSRT post-processing which creates the reason based on 3115 ## # generic events) 3116 ### assert len(swmapping) == len(allEndReasonNames), ('Incorrect number' 3117 ### ' of map pairs given in argument') 3118 ## for (reason, outcome) in swmapping.iteritems(): 3119 ## targetName = outcome[0] 3120 ## epmapping = outcome[1] 3121 ## # no checks are made on epmapping here 3122 ## assert reason not in seenReasons, ('reasons cannot appear more' 3123 ## ' than once in map domain') 3124 ## seenReasons.append(reason) 3125 ## assert reason in allEndReasonNames, ('name `'+reason+'` in map ' 3126 ## 'domain is invalid') 3127 ## if targetName != 'terminate': 3128 ## assert targetName in allDSnames, ('name `'+targetName+\ 3129 ## '` in map range is invalid') 3130 3131 3132
3133 - def _validateParameters(self):
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 # gen doesn't have any pars 3141 continue # for loop 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
3156 - def _infostr(self, verbose=1):
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 #### Private functions 3172
3173 -def getAuxVars(dsi, t, icdict, pardict):
3174 """Return auxiliary variable values evaluated at t, icdict, 3175 pardict. 3176 """ 3177 # prepare auxiliary variables eval'd at i.c. 3178 icdict_local = copy.copy(icdict) 3179 # to restore these to their former values after domain tests are done, 3180 # (in case events refer to auxiliary function initcond() 3181 old_icdict = copy.copy(dsi.model.icdict) 3182 # update icdict with auxiliary variable vals eval'd at t, for 3183 # use in event function calls for domain test. 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 # no auxiliary variables for this model 3196 # e.g. ExplicitFnGen has its auxvars defined through the 3197 # main variables 3198 auxdict = {} 3199 # restore model's original initial conditions 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 # if external inputs then include initial value of those in pardict 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 # don't overwrite any aux vars that are treated as regular variables 3240 # at this level, but fetch other aux vars and external inputs 3241 icdict.update(filteredDict(getAuxVars(MI, t, icdict, pardict), 3242 icdict.keys(), neg=True)) 3243 # override icdict with any finite-valued generator 3244 # initial condition (deliberately preset in generator definition) 3245 # for non-internal variables, and any auxiliary variables that aren't 3246 # already defined in icdict from above call to getAuxVars (must do this 3247 # to prevent NaNs being passed to numeric_to_traj. 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 # ignore 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 # ic (for internal variable) not defined for this 3260 # generator, so ignore 3261 pass 3262 # set initial conditions of MI in case auxiliary function initcond 3263 # is called during Rhs call 3264 model.set(ics=icdict) 3265 x_test = True # initial value 3266 try: 3267 # derivatives of variables are just Rhs values for ODE 3268 dxdt = dict(MI.Rhs(t, icdict, pardict)) 3269 except AttributeError: 3270 # no Rhs for the underlying generator 3271 dxdt = dxdt_zeros 3272 else: 3273 for vname, val in dxdt.items(): 3274 # hack to fix NaNs showing up for VODE Rhs function 3275 # but very large finite values come back from Dopri/Radau 3276 # and Nan's mess up call to numeric_to_traj below 3277 if not isfinite(val): 3278 dxdt[vname] = 1e308 3279 for aname in fs.auxvars: 3280 # set derivatives for remaining auxvars to dummy value 3281 dxdt[aname] = 0 3282 # print "Temp in findTrajInitiator for %s: x, dxdt ="%MI.model.name 3283 # print icdict 3284 # print dxdt 3285 for xname, x in icdict.iteritems(): 3286 # don't count time 't' or any internal variables 3287 if xname == 't' or xname in intvars or xname not in domtests: 3288 # 't' is a special value inserted into icdict, for use 3289 # in auxiliary variable evaluation at initial time, but not 3290 # a regular variable so don't do domain check here. 3291 continue 3292 xtraj = numeric_to_traj([[x], [dxdt[xname]]], 'test', [xname, 'D_'+xname], t) 3293 # print "Dom test for %s in %s"%(xname, MI.model.name) 3294 # print " discrete? ", domtests[xname].isdiscrete 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 # initial value 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 # ensure that events are cleared in case global con rules in a MI 3306 # check for events (for when used after trajectories computed) 3307 # don't use MI.get('diagnostics', xdict, t) as that will return a copy 3308 # so that clearWarnings will not propagate to underlying sub-model 3309 xdict, t, I = MI._get_initiator_cache(xdict, t) 3310 while True: 3311 # descend into sub-models until leaf found 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 # !!! Currently assumes globalConRules is a list, not a Context 3319 # (should be upgraded to that) 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 # this t, xdict is a valid initial state for this Generator 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 # only remaining possibility is a single eligible model 3352 return eligibleMI[0]
3353