1 """Euler integrator for ODE systems, with no step refinement for events.
2 """
3 from __future__ import division
4
5 from allimports import *
6 from PyDSTool.Generator import ODEsystem as ODEsystem
7 from baseclasses import Generator, theGenSpecHelper, _pollInputs
8 from PyDSTool.utils import *
9 from PyDSTool.common import *
10
11
12 from numpy import Inf, NaN, isfinite, sometrue, alltrue, sign, all, any, \
13 array, zeros, less_equal, transpose, concatenate, asarray, linspace
14 try:
15 from numpy import unique
16 except ImportError:
17
18 from numpy import unique1d as unique
19 import math, random
20 import numpy as np
21 from copy import copy, deepcopy
22 import os, platform, shutil, sys, gc
23
24 try:
25
26 import psyco
27 HAVE_PSYCO = True
28 except ImportError:
29 HAVE_PSYCO = False
30
31
34 self.f = rhs
35 self.y = None
36 self.t = None
37
39 self.y = y0
40 self.t = t0
41
43 if extraparams is None:
44 self.f_params = []
45 else:
46 self.f_params = extraparams
47
49 """Jacobian is not used"""
50 pass
51
53 """Single step"""
54 self.t += dt
55 self.y = self.y + dt*self.f(self.t, self.y, self.f_params)
56 return 1
57
60
63
65 """Euler method. Fixed step.
66
67 Uses Python target language only for functional specifications."""
68
70 if 'user_func_beforestep' in kw:
71 self.ufunc_before = kw['user_func_beforestep']
72
73 del kw['user_func_beforestep']
74 else:
75 self.ufunc_before = _dummy_userfunc
76 if 'user_func_afterstep' in kw:
77 self.ufunc_after = kw['user_func_afterstep']
78
79 del kw['user_func_afterstep']
80 else:
81 self.ufunc_after = _dummy_userfunc
82 ODEsystem.__init__(self, kw)
83 self._paraminfo = {'init_step': 'Fixed step size for time mesh.'}
84 self.diagnostics._errorcodes = {1: 'Step OK'}
85 self.diagnostics.outputStatsInfo = {'errorStatus': 'Error status on completion.'}
86 algparams_def = {'poly_interp': False,
87 'init_step': 0.01,
88 'max_pts': 100000
89 }
90 for k, v in algparams_def.iteritems():
91 if k not in self.algparams:
92 self.algparams[k] = v
93
94
96
97 ODEsystem.addMethods(self, usePsyco=HAVE_PSYCO)
98
99 self._solver = euler_solver(getattr(self,self.funcspec.spec[1]))
100 self._funcreg['_solver'] = ('self', 'euler_solver(getattr(self,' \
101 + 'self.funcspec.spec[1]))')
102
103
105 ivals = [i(solver.t) for i in inputlist]
106 s = "\n***************\nNew t, x, inputs: " + " ".join([str(s) for s in (solver.t,solver.y,ivals)])
107 s += "\ndt="+str(dt)+" f_params="+str(solver.f_params)+" dx/dt="
108 s += str(solver.f(solver.t, solver.y, sortedDictValues(self.pars)+ivals))
109 return s
110
111 - def compute(self, trajname, dirn='f', ics=None):
112 continue_integ = ODEsystem.prepDirection(self, dirn)
113 if self._dircode == -1:
114 raise NotImplementedError('Backwards integration is not implemented')
115 if ics is not None:
116 self.set(ics=ics)
117 self.validateICs()
118 self.diagnostics.clearWarnings()
119 self.diagnostics.clearErrors()
120 pnames = sortedDictKeys(self.pars)
121 xnames = self._var_ixmap
122
123 self.checkInitialConditions()
124 haveJac = int(self.haveJacobian())
125 indepdom0, indepdom1 = self.indepvariable.depdomain.get()
126 if continue_integ:
127 if indepdom0 > self._solver.t:
128 print "Previous end time is %f"%self._solver.t
129 raise ValueError("Start time not correctly updated for "
130 "continuing orbit")
131 x0 = self._solver.y
132 indepdom0 = self._solver.t
133 self.indepvariable.depdomain.set((indepdom0,indepdom1))
134 else:
135 x0 = sortedDictValues(self.initialconditions,
136 self.funcspec.vars)
137 t0 = indepdom0
138 dt = self.algparams['init_step']
139
140 solver = self._solver
141 solver.set_initial_value(x0, t0)
142 solver.dt = dt
143
144 alltData = [t0]
145 allxDataDict = dict(zip(xnames, map(listid, x0)))
146 plist = sortedDictValues(self.pars)
147 extralist = copy(plist)
148 if self.inputs:
149
150 inames = sortedDictKeys(self.inputs)
151 listend = self.numpars + len(self.inputs)
152 inputVarList = sortedDictValues(self.inputs)
153 ilist = _pollInputs(inputVarList, alltData[0]+self.globalt0,
154 self.checklevel)
155 else:
156 ilist = []
157 inames = []
158 listend = self.numpars
159 inputVarList = []
160 extralist.extend(ilist)
161 solver.set_f_params(extralist)
162 if haveJac:
163 solver.set_jac_params(extralist)
164 do_poly = self.algparams['poly_interp']
165 if do_poly:
166 rhsfn = getattr(self, self.funcspec.spec[1])
167 dx0 = rhsfn(t0, x0, extralist)
168 alldxDataDict = dict(zip(xnames, map(listid, dx0)))
169 auxvarsfn = getattr(self,self.funcspec.auxspec[1])
170
171 if not all(isfinite(self.indepvariable.depdomain.get())):
172 print "Time domain was: ", self.indepvariable.depdomain.get()
173 raise ValueError("Ensure time domain is finite")
174 if dt == indepdom1 - indepdom0:
175
176
177
178 tmesh = [indepdom0, indepdom1]
179 else:
180 notDone = True
181 while notDone:
182 tmesh = self.indepvariable.depdomain.sample(dt,
183 strict=True,
184 avoidendpoints=True)
185 notDone = False
186 first_found_t = None
187 eventslist = self.eventstruct.query(['active', 'notvarlinked'])
188 termevents = self.eventstruct.query(['term'], eventslist)
189 tmesh.pop(0)
190 xnames = self.funcspec.vars
191
192 allaDataDict = {}
193 anames = self.funcspec.auxvars
194 avals = auxvarsfn(t0, x0, extralist)
195 for aix in range(len(anames)):
196 aname = anames[aix]
197 try:
198 allaDataDict[aname] = [avals[aix]]
199 except IndexError:
200 print "\nEuler generator: There was a problem evaluating " \
201 + "an auxiliary variable"
202 print "Debug info: avals (length", len(avals), ") was ", avals
203 print "Index out of range was ", aix
204 print self.funcspec.auxspec[1]
205 print hasattr(self, self.funcspec.auxspec[1])
206 print "Args were:", [t0, x0, extralist]
207 raise
208
209 self.setEventICs(self.initialconditions, self.globalt0)
210 if self.inputs:
211 parsinps = copy(self.pars)
212 parsinps.update(dict(zip(inames,ilist)))
213 else:
214 parsinps = self.pars
215 if eventslist != []:
216 dataDict = copy(self.initialconditions)
217 dataDict.update(dict(zip(anames, avals)))
218 dataDict['t'] = t0
219 evsflagged = self.eventstruct.pollHighLevelEvents(None,
220 dataDict,
221 parsinps,
222 eventslist)
223 if len(evsflagged) > 0:
224 raise RuntimeError("Some events flagged at initial condition")
225 if continue_integ:
226
227 self.eventstruct.resetHighLevelEvents(t0, eventslist, 'prev')
228 elif self._for_hybrid_DS:
229
230
231
232
233
234
235 for evname, ev in eventslist:
236 ev.starttime = t0
237 else:
238
239 self.eventstruct.resetHighLevelEvents(t0, eventslist, None)
240 self.eventstruct.validateEvents(self.funcspec.vars + \
241 self.funcspec.auxvars + \
242 self.funcspec.inputs + \
243 ['t'], eventslist)
244 evnames = [ev[0] for ev in eventslist]
245 lastevtime = {}.fromkeys(evnames, None)
246
247 Evtimes = {}
248 Evpoints = {}
249 if continue_integ:
250 for evname in evnames:
251 try:
252
253 lastevtime[evname] = self.eventstruct.Evtimes[evname][-1] \
254 - self.globalt0
255 except (IndexError, KeyError):
256
257
258 pass
259 for evname in evnames:
260 Evtimes[evname] = []
261 Evpoints[evname] = []
262
263 depdomains = dict(zip(range(self.dimension),
264 [self.variables[xn].depdomain for xn in xnames]))
265
266 num_points = 0
267 breakwhile = False
268 while not breakwhile:
269 try:
270 new_t = tmesh.pop(0)
271 except IndexError:
272
273 break
274
275 self.ufunc_before(self)
276 try:
277 errcode = solver.integrate(dt)
278 except:
279 print "Error calling right hand side function:"
280 self.showSpec()
281 print "Numerical traceback information (current state, " \
282 + "parameters, etc.)"
283 print "in generator dictionary 'traceback'"
284 self.traceback = {'vars': dict(zip(xnames,solver.y)),
285 'pars': dict(zip(pnames,plist)),
286 'inputs': dict(zip(inames,ilist)),
287 self.indepvariable.name: new_t}
288 raise
289 avals = auxvarsfn(new_t, solver.y, extralist)
290
291
292
293 if eventslist != []:
294 dataDict = dict(zip(xnames,solver.y))
295 dataDict.update(dict(zip(anames,avals)))
296 dataDict['t'] = new_t
297 if self.inputs:
298 parsinps = copy(self.pars)
299 parsinps.update(dict(zip(inames,
300 _pollInputs(inputVarList, new_t+self.globalt0,
301 self.checklevel))))
302 else:
303 parsinps = self.pars
304 evsflagged = self.eventstruct.pollHighLevelEvents(None,
305 dataDict,
306 parsinps,
307 eventslist)
308
309
310 termevsflagged = filter(lambda e: e in evsflagged, termevents)
311 nontermevsflagged = filter(lambda e: e not in termevsflagged,
312 evsflagged)
313
314
315
316 if len(nontermevsflagged) > 0:
317 evnames = [ev[0] for ev in nontermevsflagged]
318 precEvts = self.eventstruct.query(['precise'],
319 nontermevsflagged)
320
321
322 nonprecEvts = self.eventstruct.query(['notprecise'],
323 nontermevsflagged) + precEvts
324 nonprec_evnames = [e[0] for e in nonprecEvts] + [e[0] for e in precEvts]
325
326
327 if nonprec_evnames != []:
328 temp_names = []
329 for evname, ev in nonprecEvts:
330 prevevt_time = lastevtime[evname]
331 if prevevt_time is None:
332 ignore_ev = False
333 else:
334 if solver.t-prevevt_time < ev.eventinterval:
335 ignore_ev = True
336 else:
337 ignore_ev = False
338 if not ignore_ev:
339 temp_names.append(evname)
340 lastevtime[evname] = solver.t
341 self.diagnostics.warnings.append((W_NONTERMEVENT,
342 (solver.t, temp_names)))
343 for evname in temp_names:
344 Evtimes[evname].append(solver.t)
345 xv = solver.y
346 av = array(avals)
347 Evpoints[evname].append(concatenate((xv, av)))
348 do_termevs = []
349 if termevsflagged != []:
350
351
352 for e in termevsflagged:
353 prevevt_time = lastevtime[e[0]]
354
355
356
357
358 if prevevt_time is None:
359 ignore_ev = False
360 else:
361
362 if solver.t-dt-prevevt_time < e[1].eventinterval:
363 ignore_ev = True
364
365 else:
366 ignore_ev = False
367 if not ignore_ev:
368 do_termevs.append(e)
369 if len(do_termevs) > 0:
370
371 evnames = [ev[0] for ev in do_termevs]
372 self.diagnostics.warnings.append((W_TERMEVENT, \
373 (solver.t, evnames)))
374 first_found_t = solver.t
375 for evname in evnames:
376 Evtimes[evname].append(solver.t)
377 xv = solver.y
378 av = array(avals)
379 Evpoints[evname].append(concatenate((xv, av)))
380
381 breakwhile = True
382
383
384
385 if not breakwhile:
386
387 for xi in xrange(self.dimension):
388 if not self.contains(depdomains[xi],
389 solver.y[xi],
390 self.checklevel):
391 self.diagnostics.warnings.append((W_TERMSTATEBD,
392 (solver.t,
393 xnames[xi],solver.y[xi],
394 depdomains[xi].get())))
395 breakwhile = True
396 break
397 if breakwhile:
398 break
399 alltData.append(solver.t)
400 if do_poly:
401 dxvals = rhsfn(solver.t, solver.y, extralist)
402 for xi, xname in enumerate(xnames):
403 allxDataDict[xname].append(solver.y[xi])
404 alldxDataDict[xname].append(dxvals[xi])
405 else:
406 for xi, xname in enumerate(xnames):
407 allxDataDict[xname].append(solver.y[xi])
408 for aix, aname in enumerate(anames):
409 allaDataDict[aname].append(avals[aix])
410 num_points += 1
411 if not breakwhile:
412 try:
413 extralist[self.numpars:listend] = [f(solver.t+self.globalt0,
414 self.checklevel) \
415 for f in inputVarList]
416 except ValueError:
417 print 'External input call caused value out of range error:',\
418 't = ', solver.t
419 for f in inputVarList:
420 if f.diagnostics.hasWarnings():
421 print 'External input variable %s out of range:'%f.name
422 print ' t = ', repr(f.diagnostics.warnings[-1][0]), ', ', \
423 f.name, ' = ', repr(f.diagnostics.warnings[-1][1])
424 raise
425 except AssertionError:
426 print 'External input call caused t out of range error: t = ', \
427 solver.t
428 raise
429 solver.set_f_params(extralist)
430 breakwhile = not solver.successful()
431
432 self.ufunc_after(self)
433
434
435 if first_found_t is not None:
436
437 try:
438 if self.diagnostics.warnings[-1][0] not in [W_TERMEVENT,
439 W_TERMSTATEBD]:
440 print "t =", solver.t
441 print "state =", dict(zip(xnames,solver.y))
442 raise RuntimeError("Event finding code for terminal event "
443 "failed in Generator " + self.name + \
444 ": try decreasing eventdelay or "
445 "eventinterval below eventtol, or the "
446 "atol and rtol parameters")
447 except IndexError:
448 info(self.diagnostics.outputStats, "Output statistics")
449 print "t =", solver.t
450 print "x =", solver.y
451 raise RuntimeError("Event finding failed in Generator " + \
452 self.name + ": try decreasing eventdelay "
453 "or eventinterval below eventtol")
454
455
456 for f in inputVarList:
457 for winfo in f.diagnostics.warnings:
458 self.diagnostics.warnings.append((W_NONTERMSTATEBD,
459 (winfo[0], f.name, winfo[1],
460 f.depdomain.get())))
461
462 termcount = 0
463 for (w,i) in self.diagnostics.warnings:
464 if w == W_TERMEVENT or w == W_TERMSTATEBD:
465 termcount += 1
466 if termcount > 1:
467 self.diagnostics.errors.append((E_NONUNIQUETERM,
468 (alltData[-1], i[1])))
469
470
471
472
473
474
475
476
477 variables = copyVarDict(self.variables)
478 for x in xnames:
479 if len(alltData) > 1:
480 if do_poly:
481 xvals = array([allxDataDict[x], alldxDataDict[x]]).T
482 interp = PiecewisePolynomial(alltData, xvals, 2)
483 else:
484 xvals = allxDataDict[x]
485 interp = interp1d(alltData, xvals)
486 variables[x] = Variable(interp, 't', x, x)
487 else:
488 print "Error in Generator:", self.name
489 print "t = ", alltData
490 print "x = ", allxDataDict
491 raise PyDSTool_ValueError("Fewer than 2 data points computed")
492 for a in anames:
493 if len(alltData) > 1:
494 variables[a] = Variable(interp1d(alltData, allaDataDict[a]),
495 't', a, a)
496 else:
497 print "Error in Generator:", self.name
498 print "t = ", alltData
499 print "x = ", allxDataDict
500 raise PyDSTool_ValueError("Fewer than 2 data points computed")
501 self.diagnostics.outputStats = {'last_step': dt,
502 'num_fcns': num_points,
503 'num_steps': num_points,
504 'errorStatus': errcode
505 }
506 if solver.successful():
507
508 for evname, evtlist in Evtimes.iteritems():
509 try:
510 self.eventstruct.Evtimes[evname].extend([et+self.globalt0 \
511 for et in evtlist])
512 except KeyError:
513 self.eventstruct.Evtimes[evname] = [et+self.globalt0 \
514 for et in evtlist]
515
516 self.trajevents = {}
517 for (evname, ev) in eventslist:
518 evpt = Evpoints[evname]
519 if evpt == []:
520 self.trajevents[evname] = None
521 else:
522 evpt = transpose(array(evpt))
523 self.trajevents[evname] = Pointset({'coordnames': xnames+anames,
524 'indepvarname': 't',
525 'coordarray': evpt,
526 'indepvararray': Evtimes[evname],
527 'indepvartype': float})
528 self.defined = True
529 return Trajectory(trajname, variables.values(),
530 abseps=self._abseps, globalt0=self.globalt0,
531 checklevel=self.checklevel,
532 FScompatibleNames=self._FScompatibleNames,
533 FScompatibleNamesInv=self._FScompatibleNamesInv,
534 events=self.trajevents,
535 modelNames=self.name,
536 modelEventStructs=self.eventstruct)
537 else:
538 try:
539 self.diagnostics.errors.append((E_COMPUTFAIL, (solver.t,
540 self.diagnostics._errorcodes[errcode])))
541 except TypeError:
542
543 print "Error information: ", errcode
544 self.diagnostics.errors.append((E_COMPUTFAIL, (solver.t,
545 self.diagnostics._errorcodes[0])))
546 self.defined = False
547
548
549 - def Rhs(self, t, xdict, pdict=None, asarray=True):
550 """asarray is an unused, dummy argument for compatibility with Model.Rhs"""
551
552
553
554 x = [float(val) for val in sortedDictValues(filteredDict(self._FScompatibleNames(xdict),
555 self.funcspec.vars))]
556 if pdict is None:
557 pdict = self.pars
558
559 p = sortedDictValues(pdict)
560 else:
561 p = sortedDictValues(self._FScompatibleNames(pdict))
562 i = _pollInputs(sortedDictValues(self.inputs),
563 t, self.checklevel)
564 return apply(getattr(self, self.funcspec.spec[1]), [t, x, p+i])
565
566
567 - def Jacobian(self, t, xdict, pdict=None, asarray=True):
568 """asarray is an unused, dummy argument for compatibility with
569 Model.Jacobian"""
570 if self.haveJacobian():
571
572 x = [float(val) for val in sortedDictValues(filteredDict(self._FScompatibleNames(xdict),
573 self.funcspec.vars))]
574 if pdict is None:
575 pdict = self.pars
576
577 p = sortedDictValues(pdict)
578 else:
579 p = sortedDictValues(self._FScompatibleNames(pdict))
580 i = _pollInputs(sortedDictValues(self.inputs),
581 t, self.checklevel)
582 return apply(getattr(self, self.funcspec.auxfns["Jacobian"][1]), \
583 [t, x, p+i])
584 else:
585 raise PyDSTool_ExistError("Jacobian not defined")
586
587
588 - def JacobianP(self, t, xdict, pdict=None, asarray=True):
589 """asarray is an unused, dummy argument for compatibility with
590 Model.JacobianP"""
591 if self.haveJacobian_pars():
592
593 x = [float(val) for val in sortedDictValues(filteredDict(self._FScompatibleNames(xdict),
594 self.funcspec.vars))]
595 if pdict is None:
596 pdict = self.pars
597
598 p = sortedDictValues(pdict)
599 else:
600 p = sortedDictValues(self._FScompatibleNames(pdict))
601 i = _pollInputs(sortedDictValues(self.inputs),
602 t, self.checklevel)
603 return apply(getattr(self, self.funcspec.auxfns["Jacobian_pars"][1]), \
604 [t, x, p+i])
605 else:
606 raise PyDSTool_ExistError("Jacobian w.r.t. parameters not defined")
607
608
609 - def AuxVars(self, t, xdict, pdict=None, asarray=True):
610 """asarray is an unused, dummy argument for compatibility with
611 Model.AuxVars"""
612
613 x = [float(val) for val in sortedDictValues(filteredDict(self._FScompatibleNames(xdict),
614 self.funcspec.vars))]
615 if pdict is None:
616 pdict = self.pars
617
618 p = sortedDictValues(pdict)
619 else:
620 p = sortedDictValues(self._FScompatibleNames(pdict))
621 i = _pollInputs(sortedDictValues(self.inputs),
622 t, self.checklevel)
623 return apply(getattr(self, self.funcspec.auxspec[1]), [t, x, p+i])
624
625
628
629
630
631
632
633 symbolMapDict = {}
634
635
636 theGenSpecHelper.add(Euler_ODEsystem, symbolMapDict, 'python')
637