Package PyDSTool :: Package Generator :: Module ExplicitFnGen'
[hide private]
[frames] | no frames]

Source Code for Module PyDSTool.Generator.ExplicitFnGen'

  1  # Explicit function generator
 
  2  from __future__ import division 
  3  
 
  4  from allimports import * 
  5  from baseclasses import ctsGen, theGenSpecHelper, \
 
  6       auxfn_container, _pollInputs 
  7  from PyDSTool.utils import * 
  8  from PyDSTool.common import * 
  9  from PyDSTool.Interval import uncertain 
 10  
 
 11  # Other imports
 
 12  from numpy import Inf, NaN, isfinite, sometrue, alltrue, array, arange, \
 
 13       transpose, shape 
 14  import math, random, types 
 15  from copy import copy, deepcopy 
 16  try: 
 17      # use pscyo JIT byte-compiler optimization, if available
 
 18      import psyco 
 19      HAVE_PSYCO = True 
 20  except ImportError: 
 21      HAVE_PSYCO = False 
 22  
 
 23  
 
24 -class ExplicitFnGen(ctsGen):
25 """Explicit functional form specifying a trajectory. 26 27 E.g. for an external input. This class allows parametric 28 forms of the function, but with no dependence on x or its 29 own external inputs.""" 30 _validKeys = ['globalt0', 'xdomain', 'tdata', 'tdomain', 31 'ics', 'pars', 'checklevel', 'pdomain', 'abseps'] 32 _needKeys = ctsGen._needKeys + ['varspecs'] 33 _optionalKeys = ctsGen._optionalKeys + ['tdomain', 'pars', 'pdomain', 'xdomain', 34 'xtype', 'ics', 'auxvars', 'vars', 'events', 35 'fnspecs', 'tdata', 'enforcebounds', 36 'activatedbounds', 'reuseterms'] 37
38 - def __init__(self, kw):
39 ctsGen.__init__(self, kw) 40 dispatch_list = ['varspecs', 'tdomain', 'tdata', 'xtype', 'xdomain', 41 'ics', 'allvars', 'reuseterms', 'pars', 'pdomain', 42 'fnspecs', 'target'] 43 # allow inputs only if it's empty, for compatibility with ModelConstructor 44 # which might put an empty dictionary in for this key 45 if 'inputs' in kw: 46 if kw['inputs'] != {}: 47 raise PyDSTool_KeyError('inputs option invalid for ExplicitFnGen ' 48 'class') 49 self.funcspec = ExpFuncSpec(self._kw_process_dispatch(dispatch_list, 50 kw)) 51 self.indepvartype = float 52 for s in self.funcspec.spec[0]: 53 if s.find('x[') > -1: 54 raise ValueError('Variable values cannot depend on ' 55 'other variables in explicit function specs -- ' 56 'in function:\n'+s) 57 self._kw_process_events(kw) 58 self.checkArgs(kw) 59 self.indepvariable = Variable(listid, Interval('t_domain', 60 self.indepvartype, 61 self.tdomain, self._abseps), 62 Interval('t', self.indepvartype, self.tdata, 63 self._abseps), 't') 64 self._register(self.indepvariable) 65 for x in self.funcspec.vars + self.funcspec.auxvars: 66 try: 67 xinterval=Interval(x, self.xtype[x], self.xdomain[x], self._abseps) 68 except KeyError, e: 69 raise PyDSTool_KeyError('Mismatch between declared variables ' 70 'and xspecs: ' + str(e)) 71 # placeholder variable so that this class can be 72 # copied before it is defined (listid function is a dummy) 73 self.variables[x] = Variable(None, self.indepvariable.depdomain, 74 xinterval, x) 75 self._generate_ixmaps() 76 self.auxfns = auxfn_container(self) 77 self.addMethods(usePsyco=HAVE_PSYCO)
78
79 - def addMethods(self, usePsyco=False):
80 """Add Python-specific functions to this object's methods, 81 accelerating them with psyco, if it is available.""" 82 83 # Add the auxiliary function specs to this Generator's namespace 84 for auxfnname in self.funcspec._pyauxfns: 85 fninfo = self.funcspec._pyauxfns[auxfnname] 86 if not hasattr(self, fninfo[1]): 87 # user-defined auxiliary functions 88 # (built-ins are provided explicitly) 89 try: 90 exec fninfo[0] 91 except: 92 print 'Error in supplied auxiliary function code' 93 self._funcreg[fninfo[1]] = ('self', fninfo[0]) 94 setattr(self, fninfo[1], types.MethodType(locals()[fninfo[1]], 95 self, 96 self.__class__)) 97 # user auxiliary function interface wrapper 98 try: 99 uafi_code = self.funcspec._user_auxfn_interface[auxfnname] 100 try: 101 exec uafi_code 102 except: 103 print 'Error in auxiliary function wrapper' 104 raise 105 setattr(self.auxfns, auxfnname, 106 types.MethodType(locals()[auxfnname], self.auxfns, 107 auxfn_container)) 108 self._funcreg[auxfnname] = ('', uafi_code) 109 except KeyError: 110 # not a user-defined aux fn 111 pass 112 # bind all the auxfns here 113 if HAVE_PSYCO and usePsyco: 114 psyco.bind(getattr(self, fninfo[1])) 115 try: 116 psyco.bind(self.auxfns[auxfnname]) 117 except KeyError: 118 # not a user-defined aux fn 119 pass 120 # Add the spec function to this Generator's namespace if 121 # target language is python (otherwise integrator exposes it anyway) 122 if self.funcspec.targetlang == 'python': 123 fninfo = self.funcspec.spec 124 try: 125 exec fninfo[0] 126 except: 127 print 'Error in supplied functional specification code' 128 raise 129 self._funcreg[fninfo[1]] = ('self', fninfo[0]) 130 setattr(self, fninfo[1], types.MethodType(locals()[fninfo[1]], 131 self, 132 self.__class__)) 133 if HAVE_PSYCO and usePsyco: 134 psyco.bind(getattr(self, fninfo[1])) 135 # Add the auxiliary spec function (if present) to this 136 # Generator's namespace 137 if self.funcspec.auxspec != '': 138 fninfo = self.funcspec.auxspec 139 try: 140 exec fninfo[0] 141 except: 142 print 'Error in supplied auxiliary variable code' 143 raise 144 self._funcreg[fninfo[1]] = ('self', fninfo[0]) 145 setattr(self, fninfo[1], types.MethodType(locals()[fninfo[1]], 146 self, 147 self.__class__)) 148 if HAVE_PSYCO and usePsyco: 149 psyco.bind(getattr(self, fninfo[1]))
150 151
152 - def compute(self, trajname, ics=None):
153 """Attach specification functions to callable interface.""" 154 # repeat this check (made in __init__) in case events were added since 155 assert self.eventstruct.getLowLevelEvents() == [], \ 156 "Can only pass high level events to ExplicitFnGen objects" 157 assert self.eventstruct.query(['highlevel', 'varlinked']) == [], \ 158 "Only non-variable linked events are valid for this class" 159 ## icdict_local = copy(self.initialconditions) 160 ## t0 = self.indepvariable.depdomain[0] 161 ## icdict_local['t'] = t0 162 ## have_aux = len(self.funcspec.auxvars)>0 163 ## for a in self.funcspec.auxvars: 164 ## # these functions are intended to be methods in their target 165 ## # Variable object, so expect a first argument 'self' 166 ## exec self.funcspec.auxspec[0] 167 ## if have_aux: 168 ## self.initialconditions.update(dict(zip(self.funcspec.auxvars, 169 ## apply(locals()[self.funcspec.auxspec[1]], 170 ## (None, t0, sortedDictValues(icdict_local), 171 ## sortedDictValues(self.pars))) ))) 172 if ics is not None: 173 self.set(ics=ics) 174 self.setEventICs(self.initialconditions, self.globalt0) 175 tempfs = deepcopy(self.funcspec) 176 tempvars = copyVarDict(self.variables) 177 # make unique fn for this trajectory: function definition gets executed 178 # finally in Variable.addMethods() method 179 tempspec = makeUniqueFn(copy(tempfs.spec[0]), 7, self.name) 180 tempfs.spec = tempspec 181 for x in self.funcspec.vars: 182 x_ix = self.funcspec.vars.index(x) 183 funcname = "_mapspecfn_" + x + "_" + timestamp(7) 184 funcstr = "def " + funcname + "(self, t):\n\treturn " 185 if len(self.funcspec.vars) == 1: 186 # this clause is unnecessary if [0] is ever dropped 187 # i.e. if spec would return plain scalar in 1D case 188 funcstr += tempfs.spec[1] + "(self, t, [0], " \ 189 + repr(sortedDictValues(self.pars)) + ")[0]\n" 190 else: 191 funcstr += tempfs.spec[1] + "(self, t, [0], " \ 192 + repr(sortedDictValues(self.pars)) + ")[" \ 193 + str(x_ix) + "]\n" 194 tempvars[x].setOutput((funcname, funcstr), tempfs, 195 self.globalt0, self._var_namemap, 196 copy(self.initialconditions)) 197 if self.funcspec.auxvars != []: 198 # make unique fn for this trajectory 199 tempauxspec = makeUniqueFn(copy(tempfs.auxspec[0]), 7, self.name) 200 tempfs.auxspec = tempauxspec 201 for a in self.funcspec.auxvars: 202 a_ix = self.funcspec.auxvars.index(a) 203 funcname = "_mapspecfn_" + a + "_" + timestamp(7) 204 funcstr = "def " + funcname + "(self, t):\n\treturn " 205 if len(self.funcspec.auxvars) == 1: 206 # this clause is unnecessary if [0] is ever dropped 207 # i.e. if auxspec would return plain scalar in 1D case 208 funcstr += tempfs.auxspec[1] + "(self, t, [v(t) " \ 209 + "for v in self._refvars], " \ 210 + repr(sortedDictValues(self.pars)) \ 211 + ")[0]\n" 212 else: 213 funcstr += tempfs.auxspec[1] + "(self, t, [v(t) " \ 214 + "for v in self._refvars], " \ 215 + repr(sortedDictValues(self.pars)) \ 216 + ")[" + str(a_ix) + "]\n" 217 tempvars[a].setOutput((funcname, funcstr), tempfs, 218 self.globalt0, self.funcspec.auxvars, 219 copy(self.initialconditions), 220 sortedDictValues(tempvars, 221 self.funcspec.vars)) 222 self.diagnostics.clearWarnings() 223 self.diagnostics.clearErrors() 224 # Find any events in tdomain, and adjust tdomain in case they 225 # are terminal 226 eventslist = self.eventstruct.query(['highlevel', 'active', 227 'notvarlinked']) 228 termevents = self.eventstruct.query(['term'], eventslist) 229 Evtimes = {} 230 Evpoints = {} 231 for evname, ev in eventslist: 232 Evtimes[evname] = [] 233 Evpoints[evname] = [] 234 if eventslist != []: 235 if self._for_hybrid_DS: 236 # self._for_hybrid_DS is set internally by HybridModel class 237 # to ensure not to reset events, because they may be about to 238 # flag on first step if previous hybrid state was the same 239 # generator and, for example, two variables are synchronizing 240 # so that their events get very close together. 241 # Just reset the starttimes of these events 242 for evname, ev in eventslist: 243 ev.starttime = t0 244 else: 245 self.eventstruct.resetHighLevelEvents(self.indepvariable.depdomain[0], 246 eventslist) 247 self.eventstruct.validateEvents(self.funcspec.vars + \ 248 self.funcspec.auxvars + \ 249 ['t'], eventslist) 250 for evname, ev in eventslist: 251 # select only continuous-valued variables for event detection 252 # (in case of indicator variables used in hybrid systems) 253 evsfound = ev.searchForEvents(self.indepvariable.depdomain.get(), 254 parDict=self.pars, 255 vars=copyVarDict(tempvars, only_cts=True), 256 checklevel=self.checklevel) 257 tvals = sortedDictValues(tempvars) 258 for evinfo in evsfound: 259 Evtimes[evname].append(evinfo[0]) 260 Evpoints[evname].append(array([v(evinfo[0]) for v in tvals])) 261 self.eventstruct.resetHighLevelEvents(self.indepvariable.depdomain[0], 262 eventslist) 263 self.eventstruct.validateEvents(self.funcspec.vars + \ 264 self.funcspec.auxvars + \ 265 ['t'], eventslist) 266 termevtimes = {} 267 nontermevtimes = {} 268 for evname, ev in eventslist: 269 numevs = shape(Evtimes[evname])[-1] 270 if numevs == 0: 271 continue 272 if ev.activeFlag: 273 if numevs > 1: 274 print "Event info:", Evtimes[evname] 275 assert numevs <= 1, ("Internal error: more than one " 276 "terminal event of same type found") 277 # For safety, we should assert that this event 278 # also appears in termevents, but we don't 279 if Evtimes[evname][0] in termevtimes.keys(): 280 # append event name to this warning 281 warning_ix = termevtimes[Evtimes[evname][0]] 282 self.diagnostics.warnings[warning_ix][1][1].append(evname) 283 else: 284 # make new termevtime entry for the new warning 285 termevtimes[Evtimes[evname][0]] = \ 286 len(self.diagnostics.warnings) 287 self.diagnostics.warnings.append((W_TERMEVENT, 288 (Evtimes[evname][0], 289 [evname]))) 290 else: 291 for ev in range(numevs): 292 if Evtimes[evname][ev] in nontermevtimes.keys(): 293 # append event name to this warning 294 warning_ix = nontermevtimes[Evtimes[evname][ev]] 295 self.diagnostics.warnings[warning_ix][1][1].append(evname) 296 else: 297 # make new nontermevtime entry for the new warning 298 nontermevtimes[Evtimes[evname][ev]] = \ 299 len(self.diagnostics.warnings) 300 self.diagnostics.warnings.append((W_NONTERMEVENT, 301 (Evtimes[evname][ev], 302 [evname]))) 303 termcount = 0 304 earliest_termtime = self.indepvariable.depdomain[1] 305 for (w,i) in self.diagnostics.warnings: 306 if w == W_TERMEVENT or w == W_TERMSTATEBD: 307 termcount += 1 308 if i[0] < earliest_termtime: 309 earliest_termtime = i[0] 310 # now delete any events found after the earliest terminal event, if any 311 if termcount > 0: 312 warn_temp = [] 313 for (w,i) in self.diagnostics.warnings: 314 if i[0] <= earliest_termtime: 315 warn_temp.append((w,i)) 316 self.diagnostics.warnings = warn_temp 317 self.indepvariable.depdomain.set([self.indepvariable.depdomain[0], 318 earliest_termtime]) 319 for v in tempvars.values(): 320 v.indepdomain.set(self.indepvariable.depdomain.get()) 321 ## print 'Time interval adjusted according to %s: %s' % \ 322 ## (self._warnmessages[w], str(i[0])+", "+ str(i[1])) 323 # build event pointset information (reset previous trajectory's) 324 self.trajevents = {} 325 for (evname, ev) in eventslist: 326 evpt = Evpoints[evname] 327 if evpt == []: 328 self.trajevents[evname] = None 329 else: 330 evpt = transpose(array(evpt)) 331 self.trajevents[evname] = Pointset({ 332 'coordnames': sortedDictKeys(tempvars), 333 'indepvarname': 't', 334 'coordarray': evpt, 335 'indepvararray': Evtimes[evname], 336 'indepvartype': self.indepvartype}) 337 if not self.defined: 338 self._register(self.variables) 339 self.validateSpec() 340 self.defined = True 341 return Trajectory(trajname, tempvars.values(), 342 abseps=self._abseps, globalt0=self.globalt0, 343 checklevel=self.checklevel, 344 FScompatibleNames=self._FScompatibleNames, 345 FScompatibleNamesInv=self._FScompatibleNamesInv, 346 events=self.trajevents, 347 modelNames=self.name, 348 modelEventStructs=self.eventstruct)
349 350
351 - def AuxVars(self, t, xdict, pdict=None, asarray=True):
352 """asarray is an unused, dummy argument for compatibility with 353 Model.AuxVars""" 354 # also, ensure xdict doesn't contain elements like array([4.1]) instead of 4 355 x = [float(val) for val in sortedDictValues(filteredDict(self._FScompatibleNames(xdict), 356 self.funcspec.vars))] 357 if pdict is None: 358 pdict = self.pars 359 # internal self.pars already is FS-compatible 360 p = sortedDictValues(pdict) 361 else: 362 p = sortedDictValues(self._FScompatibleNames(pdict)) 363 i = _pollInputs(sortedDictValues(self.inputs), t, self.checklevel) 364 return apply(getattr(self, self.funcspec.auxspec[1]), [t, x, p+i])
365
366 - def haveJacobian_pars(self):
367 """Report whether generator has an explicit user-specified Jacobian 368 with respect to pars associated with it.""" 369 return 'Jacobian_pars' in self.funcspec.auxfns
370
371 - def haveJacobian(self):
372 """Report whether generator has an explicit user-specified Jacobian 373 associated with it.""" 374 return 'Jacobian' in self.funcspec.auxfns
375 376
377 - def set(self, **kw):
378 """Set ExplicitFnGen parameters""" 379 if remain(kw.keys(), self._validKeys) != []: 380 raise KeyError("Invalid keys in argument") 381 if 'globalt0' in kw: 382 # pass up to generic treatment for this 383 ctsGen.set(self, globalt0=kw['globalt0']) 384 if 'checklevel' in kw: 385 # pass up to generic treatment for this 386 ctsGen.set(self, checklevel=kw['checklevel']) 387 if 'abseps' in kw: 388 # pass up to generic treatment for this 389 ctsGen.set(self, abseps=kw['abseps']) 390 # optional keys for this call are 391 # ['pars', 'tdomain', 'xdomain', 'pdomain'] 392 if 'xdomain' in kw: 393 for k_temp, v in kw['xdomain'].iteritems(): 394 k = self._FScompatibleNames(k_temp) 395 if k in self.funcspec.vars+self.funcspec.auxvars: 396 if isinstance(v, _seq_types): 397 assert len(v) == 2, \ 398 "Invalid size of domain specification for "+k 399 if v[0] >= v[1]: 400 raise PyDSTool_ValueError('xdomain values must be' 401 'in order of increasing ' 402 'size') 403 elif isinstance(v, _num_types): 404 pass 405 else: 406 raise PyDSTool_TypeError('Invalid type for xdomain spec' 407 ' '+k) 408 self.xdomain[k] = v 409 else: 410 raise ValueError('Illegal variable name') 411 try: 412 self.variables[k].depdomain.set(v) 413 except TypeError: 414 raise TypeError('xdomain must be a dictionary of variable' 415 ' names -> valid interval 2-tuples or ' 416 'singletons') 417 for ev in self.eventstruct.events.values(): 418 ev.xdomain[k] = v 419 if 'pdomain' in kw: 420 for k_temp, v in kw['pdomain'].iteritems(): 421 k = self._FScompatibleNames(k_temp) 422 if k in self.funcspec.pars: 423 if isinstance(v, _seq_types): 424 assert len(v) == 2, \ 425 "Invalid size of domain specification for "+k 426 if v[0] >= v[1]: 427 raise PyDSTool_ValueError('pdomain values must be' 428 'in order of increasing ' 429 'size') 430 else: 431 self.pdomain[k] = copy(v) 432 elif isinstance(v, _num_types): 433 self.pdomain[k] = [v, v] 434 else: 435 raise PyDSTool_TypeError('Invalid type for pdomain spec' 436 ' '+k) 437 else: 438 raise ValueError('Illegal parameter name') 439 try: 440 self.parameterDomains[k].depdomain.set(v) 441 except TypeError: 442 raise TypeError('pdomain must be a dictionary of parameter' 443 ' names -> valid interval 2-tuples or ' 444 'singletons') 445 for ev in self.eventstruct.events.values(): 446 ev.pdomain[k] = v 447 if 'tdata' in kw: 448 self.tdata = kw['tdata'] 449 if 'tdomain' in kw: 450 self.tdomain = kw['tdomain'] 451 self.indepvariable.indepdomain.set(self.tdomain) 452 if self.tdomain[0] > self.tdata[0]: 453 if self.indepvariable.indepdomain.contains(self.tdata[0]) == uncertain: 454 self.diagnostics.warnings.append((W_UNCERTVAL, 455 (self.tdata[0],self.tdomain))) 456 else: 457 print 'tdata cannot be specified below smallest '\ 458 'value in tdomain\n (possibly due to uncertain bounding).'\ 459 ' It has been automatically adjusted from\n ', self.tdata[0], \ 460 'to', self.tdomain[0], '(difference of', \ 461 self.tdomain[0]-self.tdata[0], ')' 462 self.tdata[0] = self.tdomain[0] 463 if self.tdomain[1] < self.tdata[1]: 464 if self.indepvariable.indepdomain.contains(self.tdata[1]) == uncertain: 465 self.diagnostics.warnings.append((W_UNCERTVAL, 466 (self.tdata[1],self.tdomain))) 467 else: 468 print 'tdata cannot be specified above largest '\ 469 'value in tdomain\n (possibly due to uncertain bounding).'\ 470 ' It has been automatically adjusted from\n ', \ 471 self.tdomain[1], 'to', \ 472 self.tdomain[1], '(difference of', \ 473 self.tdata[1]-self.tdomain[1], ')' 474 self.tdata[1] = self.tdomain[1] 475 self.indepvariable.depdomain.set(self.tdata) 476 if 'ics' in kw: 477 for k_temp, v in kw['ics'].iteritems(): 478 k = self._FScompatibleNames(k_temp) 479 if k in self.funcspec.vars+self.funcspec.auxvars: 480 self._xdatadict[k] = ensurefloat(v) 481 else: 482 raise ValueError('Illegal variable name') 483 self.initialconditions.update(self._xdatadict) 484 if 'pars' in kw: 485 if not self.pars: 486 raise ValueError('No pars were declared for this object' 487 ' at initialization.') 488 for k_temp, v in kw['pars'].iteritems(): 489 k = self._FScompatibleNames(k_temp) 490 if k in self.pars: 491 cval = self.parameterDomains[k].contains(v) 492 if self.checklevel < 3: 493 if cval is not notcontained: 494 self.pars[k] = ensurefloat(v) 495 if cval is uncertain and self.checklevel == 2: 496 print 'Warning: Parameter value at bound' 497 else: 498 raise PyDSTool_ValueError('Parameter value out of ' 499 'bounds') 500 else: 501 if cval is contained: 502 self.pars[k] = ensurefloat(v) 503 elif cval is uncertain: 504 raise PyDSTool_UncertainValueError('Parameter value' 505 ' at bound') 506 else: 507 raise PyDSTool_ValueError('Parameter value out of' 508 ' bounds') 509 else: 510 raise PyDSTool_AttributeError('Illegal parameter name')
511 512
513 - def validateSpec(self):
514 ctsGen.validateSpec(self) 515 try: 516 for v in self.variables.values(): 517 assert isinstance(v, Variable) 518 assert not self.inputs 519 except AssertionError: 520 print 'Invalid system specification' 521 raise
522 523
524 - def __del__(self):
525 ctsGen.__del__(self)
526 527 528 529 # Register this Generator with the database 530 531 symbolMapDict = {} 532 # in future, provide appropriate mappings for libraries math, 533 # random, etc. (for now it's left to FuncSpec) 534 theGenSpecHelper.add(ExplicitFnGen, symbolMapDict, 'python', 'ExpFuncSpec') 535