1
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
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
18 import psyco
19 HAVE_PSYCO = True
20 except ImportError:
21 HAVE_PSYCO = False
22
23
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
39 ctsGen.__init__(self, kw)
40 dispatch_list = ['varspecs', 'tdomain', 'tdata', 'xtype', 'xdomain',
41 'ics', 'allvars', 'reuseterms', 'pars', 'pdomain',
42 'fnspecs', 'target']
43
44
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
72
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
80 """Add Python-specific functions to this object's methods,
81 accelerating them with psyco, if it is available."""
82
83
84 for auxfnname in self.funcspec._pyauxfns:
85 fninfo = self.funcspec._pyauxfns[auxfnname]
86 if not hasattr(self, fninfo[1]):
87
88
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
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
111 pass
112
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
119 pass
120
121
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
136
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
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
160
161
162
163
164
165
166
167
168
169
170
171
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
178
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
187
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
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
207
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
225
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
237
238
239
240
241
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
252
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
278
279 if Evtimes[evname][0] in termevtimes.keys():
280
281 warning_ix = termevtimes[Evtimes[evname][0]]
282 self.diagnostics.warnings[warning_ix][1][1].append(evname)
283 else:
284
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
294 warning_ix = nontermevtimes[Evtimes[evname][ev]]
295 self.diagnostics.warnings[warning_ix][1][1].append(evname)
296 else:
297
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
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
322
323
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
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
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
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
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
383 ctsGen.set(self, globalt0=kw['globalt0'])
384 if 'checklevel' in kw:
385
386 ctsGen.set(self, checklevel=kw['checklevel'])
387 if 'abseps' in kw:
388
389 ctsGen.set(self, abseps=kw['abseps'])
390
391
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
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
526
527
528
529
530
531 symbolMapDict = {}
532
533
534 theGenSpecHelper.add(ExplicitFnGen, symbolMapDict, 'python', 'ExpFuncSpec')
535