1
2 from __future__ import division
3
4 from allimports import *
5 from PyDSTool.Generator import ODEsystem as ODEsystem
6 from baseclasses import Generator, theGenSpecHelper, genDB, _pollInputs
7 from PyDSTool.utils import *
8 from PyDSTool.common import *
9 from PyDSTool.ModelSpec import QuantSpec
10 from PyDSTool.integrator import integrator
11 from PyDSTool.parseUtils import addArgToCalls, wrapArgInCall
12 import PyDSTool.Redirector as redirc
13 import numpy as npy
14
15
16 from numpy import Inf, NaN, isfinite, int, int32, float, float64, \
17 sometrue, alltrue, any, all, concatenate, transpose, array, zeros
18 import math, random
19 import operator
20 from copy import copy, deepcopy
21 import os, platform, shutil, sys, gc
22
23
24 from numpy.distutils.core import setup, Extension
25 from distutils.sysconfig import get_python_inc
26
27
28 import PyDSTool
29 _pydstool_path = PyDSTool.__path__[0]
30
31
32 rout = redirc.Redirector(redirc.STDOUT)
33 rerr = redirc.Redirector(redirc.STDERR)
34
35
36
38 """Radau 5 specialization of the basic integrator class."""
39
40 - def __init__(self, modname, rhs='default_name', phaseDim=0, paramDim=0,
41 nAux=0, nEvents=0, nExtInputs=0, hasJac=0, hasJacP=0,
42 hasMass=0, extraSpace=0, defaultBound=1e8):
43 integrator.__init__(self, rhs=rhs, phaseDim=phaseDim, paramDim=paramDim,
44 nAux=nAux, nEvents=nEvents, nExtInputs=nExtInputs, hasJac=hasJac,
45 hasJacP=hasJacP, hasMass=hasMass, extraSpace=extraSpace,
46 defaultBound=defaultBound)
47 self.modname = modname
48 try:
49 self._integMod = __import__(modname, globals())
50 except:
51 print "Error in importing compiled vector field and integrator."
52 print "Did you compile the RHS C code?"
53 raise
54
55 assert 'Integrate' in dir(self._integMod), \
56 "radau library does not contain Integrate()"
57
58 self.safety = []
59 self.jacRecompute = []
60 self.newtonStop = []
61 self.stepChangeLB = []
62 self.stepChangeUB = []
63 self.stepSizeLB = []
64 self.stepSizeUB = []
65 self.hessenberg = []
66 self.maxNewton = []
67 self.newtonStart = []
68 self.index1dim = []
69 self.index2dim = []
70 self.index3dim = []
71 self.stepSizeStrategy = []
72 self.DAEstructureM1 = []
73 self.DAEstructureM2 = []
74
75 retval = self._integMod.InitBasic(self.phaseDim, self.paramDim, self.nAux,
76 self.nEvents, self.nExtInputs, self.hasJac,
77 self.hasJacP, self.hasMass, self.extraSpace)
78
79 if retval[0] != 1:
80 raise PyDSTool_InitError('Call to InitBasic failed! (radau)')
81
82 self.initBasic = True
83
84
85 - def Run(self, hinit=0, hmax=1.0, checkAux=0, calcSpecTimes=0, verbose=0,
86 safety=0.9, jacRecompute=0.001, newtonStop=-1, stepChangeLB=1,
87 stepChangeUB=1.2, stepSizeLB=0.2, stepSizeUB=8.0, hessenberg=0,
88 maxNewton=7, newtonStart=0, index1dim=-1, index2dim=0, index3dim=0,
89 stepSizeStrategy=1, DAEstructureM1=0, DAEstructureM2=0, useJac=0, useMass=0):
90 if not self.initBasic:
91 raise PyDSTool_InitError('initBasic is False (radau)')
92 if not self.initEvents:
93 raise PyDSTool_InitError('initEvents is False (radau)')
94 if not self.initIntegrate:
95 raise PyDSTool_InitError('initInteg is False (radau)')
96 if not self.setParams:
97 raise PyDSTool_InitError('setParams is False (radau)')
98 if self.nExtInputs > 0 and not self.initExtInputs:
99 raise PyDSTool_InitError('initExtInputs is False (radau)')
100
101 self.setRadauParams(hinit=hinit, hmax=hmax, checkAux=checkAux,
102 calcSpecTimes=calcSpecTimes,
103 verbose=verbose, safety=safety,
104 jacRecompute=jacRecompute, newtonStop=newtonStop,
105 stepChangeLB=stepChangeLB, stepChangeUB=stepChangeUB,
106 stepSizeLB=stepSizeLB, stepSizeUB=stepSizeUB,
107 hessenberg=hessenberg,maxNewton=maxNewton,
108 newtonStart=newtonStart, index1dim=index1dim,
109 index2dim=index2dim, index3dim=index3dim,
110 stepSizeStrategy=stepSizeStrategy,
111 DAEstructureM1=DAEstructureM1,
112 DAEstructureM2=DAEstructureM2,
113 useJac=useJac,useMass=useMass)
114
115
116 self.Reset()
117 T, P, A, Stats, H, Err, EvtT, EvtP = self._integMod.Integrate(self.ic,
118 self.t0,
119 self.hinit,
120 self.hmax,
121 self.safety,
122 self.jacRecompute,
123 self.newtonStop,
124 self.stepChangeLB,
125 self.stepChangeUB,
126 self.stepSizeLB,
127 self.stepSizeUB,
128 self.hessenberg,
129 self.maxNewton,
130 self.newtonStart,
131 self.index1dim,
132 self.index2dim,
133 self.index3dim,
134 self.stepSizeStrategy,
135 self.DAEstructureM1,
136 self.DAEstructureM2,
137 self.useJac,
138 self.useMass,
139 self.verbose,
140 self.checkAux,
141 self.calcSpecTimes)
142 self.points = P
143 self.times = T
144 self.auxPoints = A
145 self.eventTimes = EvtT
146 self.eventPoints = EvtP
147 self.errors = Err
148 self.stats = Stats
149 self.step = H
150
151 try:
152 self.lastTime = self.times[-1]
153 self.lastPoint = [self.points[i][-1] for i in range(self.phaseDim)]
154 self.lastStep = self.step
155 except IndexError:
156 self.lastTime = self.t0
157 self.lastPoint = self.ic
158 self.lastStep = self.hinit
159 self.numRuns += 1
160 self.canContinue = True
161
162 return T, P, A, Stats, H, Err, EvtT, EvtP
163
164
165 - def Continue(self, tend, params=[], calcSpecTimes=0, verbose=0,
166 extInputChanged=False, extInputVals=[], extInputTimes=[],
167 bounds=[]):
168 if not self.initBasic:
169 raise PyDSTool_InitError('initBasic is False (radau)')
170 if not self.initEvents:
171 raise PyDSTool_InitError('initEvents is False (radau)')
172 if not self.initIntegrate:
173 raise PyDSTool_InitError('initInteg is False (radau)')
174 if not self.setParams:
175 raise PyDSTool_InitError('setParams is False (radau)')
176 if self.nExtInputs > 0 and not self.initExtInputs:
177 raise PyDSTool_InitError('initExtInputs is False (radau)')
178
179 if not self.canContinue:
180 raise PyDSTool_ContError('Unable to continue trajectory -- '
181 'have you run the integrator and reset events, etc?')
182
183 self.setContParams(tend=tend, params=copy(params),
184 calcSpecTimes=calcSpecTimes, verbose=verbose, extInputChanged=extInputChanged,
185 extInputVals=copy(extInputVals), extInputTimes=copy(extInputTimes),
186 bounds=copy(bounds))
187
188
189 T, P, A, Stats, H, Err, EvtT, EvtP = \
190 self._integMod.Integrate(self.lastPoint,
191 self.lastTime,
192 self.lastStep, self.hmax,
193 self.safety,
194 self.jacRecompute,
195 self.newtonStop,
196 self.stepChangeLB,
197 self.stepChangeUB,
198 self.stepSizeLB,
199 self.stepSizeUB,
200 self.hessenberg,
201 self.maxNewton,
202 self.newtonStart,
203 self.index1dim,
204 self.index2dim,
205 self.index3dim,
206 self.stepSizeStrategy,
207 self.DAEstructureM1,
208 self.DAEstructureM2,
209 self.useJac,
210 self.useMass,
211 self.verbose,
212 self.checkAux,
213 self.calcSpecTimes)
214
215 self.points = P
216 self.times = T
217 self.auxPoints = A
218 self.eventTimes = EvtT
219 self.eventPoints = EvtP
220 self.errors = Err
221 self.stats = Stats
222 self.step = H
223
224 try:
225 self.lastTime = self.times[-1]
226 self.lastPoint = [self.points[i][-1] for i in range(self.phaseDim)]
227 self.lastStep = self.step
228 except IndexError:
229 self.lastTime = self.t0
230 self.lastPoint = self.ic
231 self.lastStep = self.hinit
232 self.numRuns += 1
233 self.numContinues += 1
234 self.canContinue = True
235
236 return T, P, A, Stats, H, Err, EvtT, EvtP
237
238
239 - def setRadauParams(self, hinit, hmax, checkAux, calcSpecTimes,
240 verbose, safety, jacRecompute, newtonStop,
241 stepChangeLB, stepChangeUB, stepSizeLB, stepSizeUB,
242 hessenberg, maxNewton, newtonStart, index1dim,
243 index2dim, index3dim, stepSizeStrategy,
244 DAEstructureM1, DAEstructureM2, useJac, useMass):
245 useJac = int(useJac)
246 useMass = int(useMass)
247 checkAux = int(checkAux)
248 calcSpecTimes = int(calcSpecTimes)
249 hessenberg = int(hessenberg)
250
251 if not isinstance(hinit, _num_types):
252 raise TypeError("hinit must be int, float")
253
254 if not isinstance(hmax, _num_types):
255 raise TypeError("hmax must be int, float")
256
257 if abs(hinit) > abs(hmax):
258 raise ValueError("Abs value of hinit (%g) must be less than hmax (%g)"%(hinit,hmax))
259
260 if not isinstance(checkAux, _int_types):
261 raise TypeError("checkAux must be int")
262 if checkAux not in (0,1):
263 raise TypeError("checkAux must be 0 or 1")
264 if checkAux == 1 and self.nAux <= 0:
265 raise ValueError("checkAux cannot be 1 if nAux is 0")
266
267 if not isinstance(verbose, _int_types):
268 raise TypeError("verbose must be int")
269 if verbose not in (0,1):
270 if verbose >= 2:
271
272 verbose = 1
273 else:
274 raise TypeError("verbose must be 0 or 1")
275
276 if not isinstance(calcSpecTimes, _int_types):
277 raise TypeError("calcSpecTimes must be int")
278 if calcSpecTimes not in (0,1):
279 raise TypeError("calcSpecTimes must be 0 or 1")
280 if calcSpecTimes == 1 and len(self.specTimes) <= 0:
281 raise ValueError("calcSpecTimes cannot be 1 if specTimes is empty")
282
283 if safety < 0:
284 raise ValueError("safety must be non-negative")
285 if jacRecompute <= 0.0:
286 raise ValueError("jacRecompute must be positive")
287 if newtonStop < 0:
288 newtonStop = 0
289 if stepChangeLB <= 0:
290 raise ValueError("stepChangeLB must be positive")
291 if stepChangeUB <= 0:
292 raise ValueError("stepChangeUB must be positive")
293 if stepSizeLB <= 0:
294 raise ValueError("stepSizeLB must be positive")
295 if stepSizeUB <= 0:
296 raise ValueError("stepSizeUB must be positive")
297
298 if stepChangeLB > stepChangeUB:
299 raise ValueError("stepChangeLB must be less than stepChangeUB")
300 if stepSizeLB >= stepSizeUB:
301 raise ValueError("stepSizeLB must be less than stepSizeUB")
302
303 if hessenberg not in (0,1):
304 raise ValueError("hessenberg must be 0 or 1")
305 if hessenberg == 1 and useMass != 0:
306 raise ValueError("hessenberg form cannot be used for implicit systems (mass matrix)")
307 if not isinstance(maxNewton, _int_types):
308 raise TypeError("maxNewton must be int")
309 if maxNewton <= 0:
310 raise ValueError("maxNewton must be positive")
311
312 if newtonStart not in (0,1):
313 raise ValueError("newtonStart must be 0 or 1")
314
315 if index1dim <= 0:
316 index1dim = self.phaseDim
317 if index2dim != 0:
318 raise ValueError("Currently index2dim must be 0")
319 if index3dim != 0:
320 raise ValueError("Currently index3dim must be 0")
321
322 if stepSizeStrategy not in (1,2):
323 raise ValueError("stepSizeStrategy must be 1 or 2")
324
325 if DAEstructureM1 != 0:
326 raise ValueError("Currently DAEstructureM1 must be 0")
327 if DAEstructureM2 != 0:
328 raise ValueError("Currently DAEstructureM2 must be 0")
329
330 if useJac not in (0,1):
331 raise ValueError("useJac must be 0 or 1")
332 if useMass not in (0,1):
333 raise ValueError("useMass must be 0 or 1")
334
335 if useJac == 1 and self.hasJac != 1:
336 raise ValueError("useJac must be 0 if hasJac is not 1")
337 if useMass == 1 and self.hasMass != 1:
338 raise ValueError("useMass must be 0 if hasMass is not 1")
339
340 self.hinit = hinit
341 self.hmax = hmax
342 self.safety = safety
343 self.jacRecompute = jacRecompute
344 self.newtonStop = newtonStop
345 self.stepChangeLB = stepChangeLB
346 self.stepChangeUB = stepChangeUB
347 self.stepSizeLB = stepSizeLB
348 self.stepSizeUB = stepSizeUB
349 self.hessenberg = hessenberg
350 self.maxNewton = maxNewton
351 self.newtonStart = newtonStart
352 self.index1dim = index1dim
353 self.index2dim = index2dim
354 self.index3dim = index3dim
355 self.stepSizeStrategy = stepSizeStrategy
356 self.DAEstructureM1 = DAEstructureM1
357 self.DAEstructureM2 = DAEstructureM2
358 self.useJac = useJac
359 self.useMass = useMass
360 self.verbose = verbose
361 self.checkAux = checkAux
362 self.calcSpecTimes = calcSpecTimes
363
364
365
367 """Wrapper for Radau integrator (with support for differential-algebraic equations).
368
369 Uses C target language only for functional specifications"""
370 _paraminfo = {'rtol': 'Relative error tolerance.',
371 'atol': 'Absolute error tolerance.',
372 'safety': 'Safety factor in the step size prediction, default 0.9.',
373 'max_step': 'Maximal step size, default tend-tstart.',
374 'init_step': 'Initial step size, default is a guess computed by the function init_step.',
375 'fac1': 'Parameter for step size selection; the new step size is chosen subject to the restriction fac1 <= new_step/old_step <= fac2. Default value is 1.0.',
376 'fac2': 'Parameter for step size selection; the new step size is chosen subject to the restriction fac1 <= new_step/old_step <= fac2. Default value is 1.2.',
377 'stepLB': '',
378 'stepUB': '',
379 'refine': 'Refine output by adding points interpolated using the RK4 polynomial (0, 1 or 2).',
380 'step_strategy': """Switch for step size strategy;
381 If step_strategy=1 mod. predictive controller (Gustafsson).
382 If step_strategy=2 classical step size control.
383 The default value (for step_strategy=0) is step_strategy=1.
384 the choice step_strategy=1 seems to produce safer results;
385 for simple problems, the choice step_strategy=2 produces
386 often slightly faster runs.""",
387 'jac_recompute': """Decides whether the Jacobian should be recomputed;
388 increase jac_recompute to 0.1 say, when Jacobian evaluations
389 are costly. for small systems jac_recompute should be smaller
390 (0.001, say). negative jac_recompute forces the code to
391 compute the Jacobian after every accepted step.
392 Default 0.001.""",
393 'newton_start': "",
394 'newton_stop': "",
395 'max_newton': "Maximum number of Newton iterations to take in solving the implicit system at each step (default 7)",
396 'DAEstructureM1': "",
397 'DAEstructureM2': "",
398 'hessenberg': "",
399 'index1dim': "",
400 'index2dim': "",
401 'index3dim': "",
402 'use_special': "Switch for using special times",
403 'specialtimes': "List of special times to use during integration",
404 'check_aux': "Switch",
405 'extraspace': ""
406 }
407
409 """Use the nobuild key to postpone building of the library, e.g. in
410 order to provide additional build options to makeLibSource and
411 compileLib methods or to make changes to the C code by hand.
412 No build options can be specified otherwise."""
413
414 if 'nobuild' in kw:
415 nobuild = kw['nobuild']
416
417 del kw['nobuild']
418 else:
419 nobuild = False
420 ODEsystem.__init__(self, kw)
421 self._solver = None
422 self.diagnostics._errorcodes = {
423 0: 'Unrecognized error code returned (see stderr output)',
424 -1 : 'input is not consistent',
425 -2 : 'larger nmax is needed',
426 2 : 'larger nmax or maxevtpts is probably needed (error raised by solout)',
427 -3 : 'step size becomes too small',
428 -4 : 'the matrix is repeatedly singular (interrupted)',
429 -8 : 'The solution exceeded a magbound (poor choice of initial step)'}
430 self.diagnostics.outputStatsInfo = {
431 'last_step': 'Predicted step size of the last accepted step (useful for a subsequent call to radau).',
432 'num_steps': 'Number of used steps.',
433 'num_accept': 'Number of accepted steps.',
434 'num_reject': 'Number of rejected steps.',
435 'num_fcns': 'Number of function evaluations.',
436 'num_jacs': 'Number of Jacobian evaluations.',
437 'num_dec': 'Number of LU-decompositions',
438 'num_subs': 'Number of forward-backward substitutions',
439 'errorStatus': 'Error status on completion.'
440 }
441
442
443 algparams_def = {'poly_interp': False,
444 'init_step': 0,
445 'max_step': 0,
446 'rtol': [1e-9 for i in range(self.dimension)],
447 'atol': [1e-12 for i in range(self.dimension)],
448 'fac1': 1.0,
449 'fac2': 1.2,
450 'stepLB': 0.2,
451 'stepUB': 8.0,
452 'safety': 0.9,
453 'max_pts': 10000,
454 'refine': 0,
455 'maxbisect': [],
456 'maxevtpts': 1000,
457 'eventInt': [],
458 'eventDelay': [],
459 'eventTol': [],
460 'use_special': 0,
461 'specialtimes': [],
462 'check_aux': 1,
463 'extraspace': 100,
464 'verbose': 0,
465 'jac_recompute': 0.001,
466 'step_strategy': 1,
467 'index1dim': -1,
468 'index2dim': 0,
469 'index3dim': 0,
470 'DAEstructureM1': 0,
471 'DAEstructureM2': 0,
472 'hessenberg': 0,
473 'newton_start': 0,
474 'newton_stop': -1,
475 'max_newton': 7,
476 'hasJac': 0,
477 'hasJacP': 0,
478 'checkBounds': self.checklevel
479 }
480 for k, v in algparams_def.iteritems():
481 if k not in self.algparams:
482 self.algparams[k] = v
483
484
485 if len(self.algparams) != len(algparams_def):
486 raise ValueError("Invalid keys present in algparams argument: " \
487 + str(remain(self.algparams.keys(),algparams_def.keys())))
488
489 if self.haveMass():
490 mspec = self.funcspec.auxfns['massMatrix']
491 lensig = len(mspec[1])
492 body_str = mspec[0][lensig:].replace('\n','')
493 qbody = QuantSpec('__body__', body_str, treatMultiRefs=False,
494 ignoreSpecial=['[',']','{','}'])
495 self._const_massmat = intersect(['Y_','t'], qbody.usedSymbols) == []
496 else:
497 self._const_massmat = True
498
499 thisplatform = platform.system()
500 if thisplatform == 'Windows':
501 self._dllext = ".pyd"
502 elif thisplatform in ['Linux', 'IRIX', 'Solaris', 'SunOS', 'MacOS', 'Darwin']:
503 self._dllext = '.so'
504 else:
505 print "Shared library extension not tested on this platform."
506 print "If this process fails please report the errors to the"
507 print "developers."
508 self._dllext = '.so'
509 self._compilation_tempdir = os.path.join(os.getcwd(),
510 "radau5_temp")
511 if not os.path.isdir(self._compilation_tempdir):
512 try:
513 assert not os.path.isfile(self._compilation_tempdir), \
514 "A file already exists with the same name"
515 os.mkdir(self._compilation_tempdir)
516 except:
517 print "Could not create compilation temp directory " + \
518 self._compilation_tempdir
519 raise
520 self._compilation_sourcedir = os.path.join(_pydstool_path,"integrator")
521 self._vf_file = self.name+"_vf.c"
522 self._vf_filename_ext = "_"+self._vf_file[:-2]
523 self._prepareEventSpecs()
524 if not (os.path.isfile(os.path.join(os.getcwd(),
525 "radau5"+self._vf_filename_ext+".py")) and \
526 os.path.isfile(os.path.join(os.getcwd(),
527 "_radau5"+self._vf_filename_ext+self._dllext))):
528 if not nobuild:
529 self.makeLibSource()
530 self.compileLib()
531 else:
532 print "Build the library using the makeLib method, or in "
533 print "stages using the makeLibSource and compileLib methods."
534 self._inputVarList = []
535 self._inputTimeList = []
536
537
539 """forceLibRefresh should be called after event contents are changed,
540 or alterations are made to the right-hand side of the ODEs.
541
542 Currently this function does NOT work!"""
543
544
545 delfiles = True
546 self._solver = None
547 try:
548 del(sys.modules["_radau5"+self._vf_filename_ext])
549 del(sys.modules["radau5"+self._vf_filename_ext])
550
551 except NameError:
552
553 delfiles = False
554 if delfiles:
555 gc.collect()
556
557
558
559
560
561
562
563
564
565 print "Cannot rebuild library without restarting session. Sorry."
566 print "Try asking the Python developers to make a working module"
567 print "unimport function!"
568
569
570
572 eventActive = []
573 eventTerm = []
574 eventDir = []
575 eventDelay = []
576 eventTol = []
577 maxbisect = []
578 eventInt = []
579
580 self._eventNames = self.eventstruct.sortedEventNames()
581 for evname in self._eventNames:
582 ev = self.eventstruct.events[evname]
583 assert isinstance(ev, LowLevelEvent), ("Radau can only "
584 "accept low level events")
585
586
587 maxstep = self.algparams['max_step']
588 for evname in self._eventNames:
589 ev = self.eventstruct.events[evname]
590 eventActive.append(int(ev.activeFlag))
591 eventTerm.append(int(ev.termFlag))
592 eventDir.append(ev.dircode)
593 eventInt.append(ev.eventinterval)
594 eventDelay.append(ev.eventdelay)
595 if ev.preciseFlag:
596 eventTol.append(ev.eventtol)
597 maxbisect.append(ev.bisectlimit)
598 else:
599 eventTol.append(maxstep*1.5)
600 maxbisect.append(1)
601 self.algparams['eventTol'] = eventTol
602 self.algparams['eventDelay'] = eventDelay
603 self.algparams['eventInt'] = eventInt
604 self.algparams['maxbisect'] = maxbisect
605 self.algparams['eventActive'] = eventActive
606 self.algparams['eventTerm'] = eventTerm
607 self.algparams['eventDir'] = eventDir
608
609
610 - def makeLib(self, libsources=[], libdirs=[], include=[]):
611 """makeLib calls makeLibSource and then the compileLib method.
612 To postpone compilation of the source to a DLL, call makelibsource()
613 separately."""
614 self.makeLibSource(include)
615 self.compileLib(libsources, libdirs)
616
617
619 """makeLibSource generates the C source for the vector field specification.
620 It should be called only once per vector field."""
621
622
623 assert isinstance(include, list), "includes must be in form of a list"
624
625 STDLIB = 0
626 USERLIB = 1
627 libinclude = dict([('Python.h', STDLIB), ('math.h', STDLIB), ('stdio.h', STDLIB), ('stdlib.h', STDLIB),
628 ('string.h', STDLIB), ('vfield.h', USERLIB), ('events.h', USERLIB),
629 ('signum.h', USERLIB), ('maxmin.h', USERLIB)])
630 include_str = ''
631 for libstr, libtype in libinclude.iteritems():
632 if libtype == STDLIB:
633 quoteleft = '<'
634 quoteright = '>'
635 else:
636 quoteleft = '"'
637 quoteright = '"'
638 include_str += "#include " + quoteleft + libstr + quoteright + "\n"
639 if include != []:
640 assert isUniqueSeq(include), "list of library includes must not contain repeats"
641 for libstr in include:
642 if libstr in libinclude:
643
644 print "Warning: library '" + libstr + "' already appears in list"\
645 + " of imported libraries"
646 else:
647 include_str += "#include " + '"' + libstr + '"\n'
648 allfilestr = "/* Vector field function and events for Radau integrator.\n" \
649 + " This code was automatically generated by PyDSTool, but may be modified " \
650 + "by hand. */\n\n" + include_str + """
651 extern double *gICs;
652 extern double **gBds;
653 extern double globalt0;
654
655 static double pi = 3.1415926535897931;
656
657 double signum(double x)
658 {
659 if (x<0) {
660 return -1;
661 }
662 else if (x==0) {
663 return 0;
664 }
665 else if (x>0) {
666 return 1;
667 }
668 else {
669 /* must be that x is Not-a-Number */
670 return x;
671 }
672 }
673
674 """
675 pardefines = ""
676 vardefines = ""
677 auxvardefines = ""
678 inpdefines = ""
679
680 vnames = self._var_ixmap
681 pnames = self.funcspec.pars
682 inames = self.funcspec.inputs
683 pnames.sort()
684 inames.sort()
685 for i in xrange(self.numpars):
686 p = pnames[i]
687
688 pardefines += self.funcspec._defstr+" "+p+"\tp_["+str(i)+"]\n"
689 for i in xrange(self.dimension):
690 v = vnames[i]
691
692 vardefines += self.funcspec._defstr+" "+v+"\tY_["+str(i)+"]\n"
693 for i, v in enumerate(self.funcspec.auxvars):
694 auxvardefines += self.funcspec._defstr+" "+v+"\t("+self.funcspec._auxdefs_parsed[v]+")\n"
695 for i in xrange(len(self.funcspec.inputs)):
696 inp = inames[i]
697
698 inpdefines += self.funcspec._defstr+" "+inp+"\txv_["+str(i)+"]\n"
699 allfilestr += "\n/* Variable, aux variable, parameter, and input definitions: */ \n" \
700 + pardefines + vardefines + auxvardefines + inpdefines + "\n"
701
702 allevs = ""
703 if self._eventNames == []:
704 numevs = 0
705 else:
706 numevs = len(self._eventNames)
707 for evname in self._eventNames:
708 ev = self.eventstruct.events[evname]
709 evfullfn = ""
710 assert isinstance(ev, LowLevelEvent), ("Radau can only "
711 "accept low level events")
712 evsig = ev._LLreturnstr + " " + ev.name + ev._LLargstr
713 assert ev._LLfuncstr.index(';') > 1, ("Event function code "
714 "error: Have you included a ';' character at the end of"
715 "your 'return' statement?")
716 fbody = ev._LLfuncstr
717
718
719 if self.funcspec.auxfns:
720 fbody_parsed = addArgToCalls(fbody,
721 self.funcspec.auxfns.keys(),
722 "p_, wk_, xv_")
723 if 'initcond' in self.funcspec.auxfns:
724
725
726 fbody_parsed = wrapArgInCall(fbody_parsed,
727 'initcond', '"')
728 else:
729 fbody_parsed = fbody
730 evbody = " {\n" + fbody_parsed + "\n}\n\n"
731 allevs += evsig + evbody
732 allfilestr += evsig + ";\n"
733
734 if self.funcspec.auxfns:
735 allfilestr += "\n"
736 for finfo in self.funcspec.auxfns.values():
737 allfilestr += finfo[1] + ";\n"
738 assignEvBody = ""
739 for evix in range(numevs):
740 assignEvBody += "events[%d] = &%s;\n"%(evix,self._eventNames[evix])
741 allfilestr += "\nint N_EVENTS = " + str(numevs) + ";\nvoid assignEvents(" \
742 + "EvFunType *events){\n " + assignEvBody \
743 + "\n}\n\nvoid auxvars(unsigned, unsigned, double, double*, double*, " \
744 + "double*, unsigned, double*, unsigned, double*);\n" \
745 + """void jacobian(unsigned, unsigned, double, double*, double*, double**, unsigned, double*, unsigned, double*);
746 void jacobianParam(unsigned, unsigned, double, double*, double*, double**, unsigned, double*, unsigned, double*);
747 """
748 if self.funcspec.auxvars == []:
749 allfilestr += "int N_AUXVARS = 0;\n\n\n"
750 else:
751 allfilestr += "int N_AUXVARS = " + str(len(self.funcspec.auxvars)) \
752 + ";\n\n\n"
753 if self.funcspec.inputs == []:
754 allfilestr += "int N_EXTINPUTS = 0;\n\n\n"
755 else:
756 allfilestr += "int N_EXTINPUTS = " + str(len(self.funcspec.inputs)) \
757 + ";\n\n\n"
758 allfilestr += self.funcspec.spec[0] + "\n\n"
759 if self.funcspec.auxfns:
760 for fname, finfo in self.funcspec.auxfns.iteritems():
761 fbody = finfo[0]
762
763 fbody_parsed = addArgToCalls(fbody,
764 self.funcspec.auxfns.keys(),
765 "p_, wk_, xv_", notFirst=fname)
766 if 'initcond' in self.funcspec.auxfns:
767
768
769
770 fbody_parsed = wrapArgInCall(fbody_parsed,
771 'initcond', '"', notFirst=True)
772 allfilestr += "\n" + fbody_parsed + "\n\n"
773
774
775 allfilestr += self.funcspec.auxspec[0] + allevs
776
777 if not self.haveMass():
778 allfilestr += """
779 void massMatrix(unsigned n_, unsigned np_, double t, double *Y_, double *p_, double **f_, unsigned wkn_, double *wk_, unsigned xvn_, double *xv_) {
780 }
781 """
782 if not self.haveJacobian():
783 allfilestr += """
784 void jacobian(unsigned n_, unsigned np_, double t, double *Y_, double *p_, double **f_, unsigned wkn_, double *wk_, unsigned xvn_, double *xv_) {
785 }
786 """
787 if not self.haveJacobian_pars():
788 allfilestr += """
789 void jacobianParam(unsigned n_, unsigned np_, double t, double *Y_, double *p_, double **f_, unsigned wkn_, double *wk_, unsigned xvn_, double *xv_) {
790 }
791 """
792
793 vffile = os.path.join(self._compilation_tempdir, self._vf_file)
794 try:
795 file = open(vffile, 'w')
796 file.write(allfilestr)
797 file.close()
798 except IOError, e:
799 print "Error opening file "+self._vf_file+" for writing"
800 raise IOError, e
801
802
804 """compileLib generates a python extension DLL with integrator and vector
805 field compiled and linked.
806
807 libsources list allows additional library sources to be linked.
808 libdirs list allows additional directories to be searched for
809 precompiled libraries."""
810
811 if os.path.isfile(os.path.join(os.getcwd(),
812 "_radau5"+self._vf_filename_ext+self._dllext)):
813
814
815 proceed = False
816 print "\n"
817 print "-----------------------------------------------------------"
818 print "Present limitation of Python: Cannot rebuild library"
819 print "without exiting Python and deleting the shared library"
820 print " " + str(os.path.join(os.getcwd(),
821 "_radau5"+self._vf_filename_ext+self._dllext))
822 print "by hand! If you made any changes to the system you should"
823 print "not proceed with running the integrator until you quit"
824 print "and rebuild."
825 print "-----------------------------------------------------------"
826 print "\n"
827 else:
828 proceed = True
829 if not proceed:
830 print "Did not compile shared library."
831 return
832 if self._solver is not None:
833 self.forceLibRefresh()
834 vffile = os.path.join(self._compilation_tempdir, self._vf_file)
835 try:
836 ifacefile_orig = open(os.path.join(self._compilation_sourcedir,
837 "radau5.i"), 'r')
838 ifacefile_copy = open(os.path.join(self._compilation_tempdir,
839 "radau5_"+self._vf_file[:-2]+".i"), 'w')
840 firstline = ifacefile_orig.readline()
841 ifacefile_copy.write('%module radau5_'+self._vf_file[:-2]+'\n')
842 iffilestr = ifacefile_orig.read()
843 ifacefile_copy.write(iffilestr)
844 ifacefile_orig.close()
845 ifacefile_copy.close()
846 except IOError:
847 print "radau5.i copying error in radau compilation directory"
848 raise
849 swigfile = os.path.join(self._compilation_tempdir,
850 "radau5"+self._vf_filename_ext+".i")
851 radwrapfile = os.path.join(self._compilation_sourcedir, "radau5mod.c")
852 if self._const_massmat:
853 radfile = os.path.join(self._compilation_sourcedir, "radau5.f")
854 else:
855
856 radfile = os.path.join(self._compilation_sourcedir, "radau5v.f")
857 dc_lapackfile = os.path.join(self._compilation_sourcedir, "dc_lapack.f")
858 lapackfile = os.path.join(self._compilation_sourcedir, "lapack.f")
859 lapackcfile = os.path.join(self._compilation_sourcedir, "lapackc.f")
860 integfile = os.path.join(self._compilation_sourcedir, "integration.c")
861 interfacefile = os.path.join(self._compilation_sourcedir, "interface.c")
862 eventfile = os.path.join(self._compilation_sourcedir, "eventFinding.c")
863 memfile = os.path.join(self._compilation_sourcedir, "memory.c")
864
865
866
867
868
869 if not (all([os.path.isfile(os.path.join(self._compilation_tempdir,
870 sf)) for sf in ['radau5'+self._vf_filename_ext+'_wrap.o',
871 'lib_radau5'+self._vf_filename_ext+'.a',
872 'radau5'+self._vf_filename_ext+'.py',
873 '_radau5'+self._vf_filename_ext+'.def']])):
874 modfilelist = [swigfile]
875 else:
876 modfilelist = []
877
878 modfilelist.extend([radwrapfile, vffile, integfile, eventfile,
879 interfacefile, memfile])
880 modfilelist.extend(libsources)
881 fortfilelist=[radfile, lapackcfile, dc_lapackfile, lapackfile]
882
883
884 script_args = ['-q', 'build', '--build-lib=.',
885 '-tradau5_temp',
886 '--build-base=radau5_temp',
887 '--build-platlib=radau5_temp', '--build-purelib=radau5_temp']
888 if self._compiler != '':
889 script_args.append('-c'+str(self._compiler))
890
891
892 narraydir = npy.get_numarray_include()
893 npydir = npy.get_include()
894
895 incdirs = [npydir, narraydir, os.getcwd(), self._compilation_tempdir,
896 self._compilation_sourcedir]
897 incdirs.extend(libdirs)
898 radlibdirs = [self._compilation_sourcedir, self._compilation_tempdir,
899 os.getcwd()]
900
901 rout.start()
902 try:
903 distobject = setup(name = "Radau5 integrator",
904 author = "PyDSTool (automatically generated)",
905 script_args = script_args,
906 ext_modules = [Extension("_radau5"+self._vf_filename_ext,
907 sources=modfilelist,
908 include_dirs=incdirs,
909 library_dirs=radlibdirs,
910 libraries=['radau5'],
911 extra_compile_args=['-w', '-D__RADAU__', '-m32'],
912 extra_link_args=['-w', '-m32'])],
913
914 libraries=[('radau5',{'sources': fortfilelist,
915 'library_dirs': radlibdirs+['./']})])
916 except:
917 rout.stop()
918 print "\nError occurred in generating Radau system..."
919 print sys.exc_info()[0], sys.exc_info()[1]
920 raise RuntimeError
921 rout.stop()
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948 try:
949
950 distdestdir = distutil_destination()
951 if swigfile in modfilelist or not \
952 os.path.isfile(os.path.join(self._compilation_tempdir,
953 "radau5"+self._vf_filename_ext+".py")):
954
955 shutil.move(os.path.join(os.getcwd(),
956 self._compilation_tempdir, distdestdir,
957 "radau5_temp",
958 "radau5"+self._vf_filename_ext+".py"),
959 os.path.join(os.getcwd(),
960 "radau5"+self._vf_filename_ext+".py"))
961 except:
962 print "\nError occurred in generating Radau system"
963 print "(while moving library extension modules to CWD)"
964
965 raise
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984 - def compute(self, trajname, dirn='f', ics=None):
985 continue_integ = ODEsystem.prepDirection(self, dirn)
986
987 if ics is not None:
988 self.set(ics=ics)
989 self.validateICs()
990 self.diagnostics.clearWarnings()
991 self.diagnostics.clearErrors()
992 if isinstance(self.algparams['rtol'], list):
993 if len(self.algparams['rtol']) != self.dimension:
994 raise ValueError('rtol list must have same length as phase dimension')
995 else:
996 rtol = self.algparams['rtol']
997 self.algparams['rtol'] = [rtol for i in range(self.dimension)]
998 if isinstance(self.algparams['atol'], list):
999 if len(self.algparams['atol']) != self.dimension:
1000 raise ValueError('atol list must have same length as phase dimension')
1001 else:
1002 atol = self.algparams['atol']
1003 self.algparams['atol'] = [atol for i in range(self.dimension)]
1004 anames = self.funcspec.auxvars
1005
1006 self.checkInitialConditions()
1007 self.setEventICs(self.initialconditions, self.globalt0)
1008
1009 self._prepareEventSpecs()
1010
1011 t0 = self.indepvariable.depdomain[0]
1012 t1 = self.indepvariable.depdomain[1]
1013 plist = sortedDictValues(self.pars)
1014 self.algparams['hasJac'] = self.haveJacobian()
1015 self.algparams['hasJacP'] = self.haveJacobian_pars()
1016 self._ensure_solver()
1017 if self._dircode == 1:
1018 tbegin = t0
1019 tend = t1
1020 elif self._dircode == -1:
1021
1022 tbegin = t1
1023 tend = t0
1024 if len(self.algparams['specialtimes'])>0:
1025 use_special = self.algparams['use_special']
1026 else:
1027 use_special = 0
1028 bounds = [[],[]]
1029 for v in self.funcspec.vars:
1030 bds = self.xdomain[v]
1031 try:
1032 bounds[0].append(bds[0])
1033 bounds[1].append(bds[1])
1034 except TypeError:
1035 print v, type(bds), bds
1036 print self.xdomain
1037 raise
1038 for p in self.funcspec.pars:
1039 bds = self.pdomain[p]
1040 try:
1041 bounds[0].append(bds[0])
1042 bounds[1].append(bds[1])
1043 except TypeError:
1044 print type(bds), bds
1045 raise
1046 if continue_integ:
1047 x0 = self._solver.lastPoint
1048
1049 tbegin = self._solver.lastTime
1050 if abs(self._solver.lastStep) < abs(self.algparams['init_step']):
1051 self.algparams['init_step'] = self._solver.lastStep
1052 if abs(t1-tbegin) < abs(self.algparams['init_step']):
1053 raise ValueError("Integration end point too close to initial "
1054 "point")
1055
1056
1057
1058
1059
1060
1061
1062 else:
1063 if self._solver.numRuns > 0:
1064 self._solver.clearAll()
1065 x0 = sortedDictValues(self.initialconditions, self.funcspec.vars)
1066 self._solver.setInteg(maxpts=self.algparams['max_pts'],
1067 rtol=self.algparams['rtol'], atol=self.algparams['atol'])
1068 self._solver.setRunParams(ic=x0, params=plist,
1069 t0=tbegin, tend=tend, gt0=self.globalt0,
1070 refine=self.algparams['refine'],
1071 specTimes=self.algparams['specialtimes'],
1072 bounds=bounds)
1073 if self.inputs:
1074
1075
1076
1077 self._ensure_inputs(self._extInputsChanged)
1078
1079 if len(anames)>0:
1080 check_aux = self.algparams['check_aux']
1081 else:
1082 check_aux = 0
1083 if self.algparams['max_step'] == 0:
1084 max_step = abs(tend-tbegin)
1085 else:
1086 max_step = self.algparams['max_step']
1087 init_step = self.algparams['init_step']
1088 if self._dircode == 1:
1089 if init_step < 0:
1090 init_step = -init_step
1091 if max_step < 0:
1092 max_step = -max_step
1093 else:
1094 if init_step > 0:
1095 init_step = -init_step
1096 if max_step > 0:
1097 max_step = -max_step
1098 if continue_integ:
1099
1100 old_highest_ix = self._solver.points.shape[1]
1101 alltData, X, A, Stats, H, Err, Evtimes, \
1102 Evpoints = self._solver.Continue(tend, plist,
1103 use_special, self.algparams['verbose'],
1104 self._extInputsChanged,
1105 deepcopy(self._inputVarList),
1106 deepcopy(self._inputTimeList), bounds)
1107 else:
1108 old_highest_ix = 0
1109 self._solver.setEvents(eventActive=self.algparams['eventActive'],
1110 eventTerm=self.algparams['eventTerm'],
1111 eventDir=self.algparams['eventDir'],
1112 eventDelay=self.algparams['eventDelay'],
1113 eventInt=self.algparams['eventInt'],
1114 eventTol=self.algparams['eventTol'],
1115 maxevtpts=self.algparams['maxevtpts'],
1116 maxbisect=self.algparams['maxbisect'])
1117 alltData, X, A, Stats, H, Err, Evtimes, \
1118 Evpoints = self._solver.Run(init_step,
1119 max_step,
1120 check_aux,
1121 use_special,
1122 self.algparams['verbose'],
1123 self.algparams['safety'],
1124 self.algparams['jac_recompute'],
1125 self.algparams['newton_stop'],
1126 self.algparams['fac1'],
1127 self.algparams['fac2'],
1128 self.algparams['stepLB'],
1129 self.algparams['stepUB'],
1130 self.algparams['hessenberg'],
1131 self.algparams['max_newton'],
1132 self.algparams['newton_start'],
1133 self.algparams['index1dim'],
1134 self.algparams['index2dim'],
1135 self.algparams['index3dim'],
1136 self.algparams['step_strategy'],
1137 self.algparams['DAEstructureM1'],
1138 self.algparams['DAEstructureM2'],
1139 self.haveJacobian(),
1140 self.haveMass())
1141 self._extInputsChanged = False
1142 self.diagnostics.outputStats = {'last_step': H,
1143 'last_time': self._solver.lastTime,
1144 'last_point': self._solver.lastPoint,
1145 'num_fcns': Stats[0],
1146 'num_jacs': Stats[1],
1147 'num_steps': Stats[2],
1148 'num_accept': Stats[3],
1149 'num_reject': Stats[4],
1150 'num_dec': Stats[5],
1151 'num_subs': Stats[6],
1152 'errorStatus': Err
1153 }
1154 if self._dircode == -1:
1155
1156 alltData = alltData[::-1]
1157 X = X[:,::-1]
1158 if anames != []:
1159 A = A[:,::-1]
1160 xnames = self._var_ixmap
1161
1162
1163
1164
1165
1166
1167
1168 eventslist = self.eventstruct.query(['lowlevel', 'active'])
1169 termevents = self.eventstruct.query(['term'], eventslist)
1170 if self._eventNames != []:
1171
1172
1173
1174
1175 termevtimes = {}
1176 nontermevtimes = {}
1177 try:
1178 for evix in range(len(self._eventNames)):
1179 if Evpoints[evix] is None:
1180 continue
1181 evname = self._eventNames[evix]
1182 numevs = len(Evtimes[evix])
1183 if self.algparams['eventTerm'][evix]:
1184 if numevs > 1:
1185 print "Event info:", Evpoints, Evtimes
1186 assert numevs <= 1, ("Internal error: more than one "
1187 "terminal event of same type found")
1188
1189
1190 if Evtimes[evix][0] in termevtimes.keys():
1191
1192 warning_ix = termevtimes[Evtimes[evix][0]]
1193 self.diagnostics.warnings[warning_ix][1][1].append(evname)
1194 else:
1195
1196 termevtimes[Evtimes[evix][0]] = \
1197 len(self.diagnostics.warnings)
1198 self.diagnostics.warnings.append((W_TERMEVENT,
1199 (Evtimes[evix][0],
1200 [evname])))
1201 else:
1202 for ev in range(numevs):
1203 if Evtimes[evix][ev] in nontermevtimes.keys():
1204
1205 warning_ix = nontermevtimes[Evtimes[evix][ev]]
1206 self.diagnostics.warnings[warning_ix][1][1].append(evname)
1207 else:
1208
1209 nontermevtimes[Evtimes[evix][ev]] = \
1210 len(self.diagnostics.warnings)
1211 self.diagnostics.warnings.append((W_NONTERMEVENT,
1212 (Evtimes[evix][ev],
1213 [evname])))
1214 except IndexError:
1215 print "Events returned from integrator are the wrong size."
1216 print " Did you change the system and not refresh the C " \
1217 + "library using the forcelibrefresh() method?"
1218 raise
1219 termcount = 0
1220 for (w,i) in self.diagnostics.warnings:
1221 if w == W_TERMEVENT or w == W_TERMSTATEBD:
1222 if termcount > 0:
1223 raise ValueError("Internal error: more than one terminal "
1224 "event found")
1225 termcount += 1
1226
1227 if self._dircode > 0:
1228 compare = operator.lt
1229 last_ix = Inf
1230 else:
1231 compare = operator.gt
1232 last_ix = -Inf
1233 highest_ix = X.shape[1]-1
1234 last_t = Inf
1235 if self.algparams['checkBounds'] > 0:
1236
1237 depdomains = dict(zip(range(self.dimension),
1238 [self.variables[xn].depdomain for xn in xnames]))
1239 offender_ix = None
1240 for xi in xrange(self.dimension):
1241 if not any(depdomains[xi].isfinite()):
1242
1243 continue
1244 next_last_ix = array_bounds_check(X[xi][old_highest_ix:],
1245 depdomains[xi], self._dircode) + old_highest_ix
1246 if compare(next_last_ix, last_ix):
1247
1248
1249 last_ix = next_last_ix
1250 offender_ix = xi
1251 if not isfinite(last_ix) and last_ix < 0:
1252
1253 last_ix = Inf
1254 elif last_ix >= 0 and last_ix < highest_ix:
1255
1256 last_t = alltData[last_ix]
1257 print "Warning; domain bound reached (because algparams['checkBounds'] > 0)"
1258 self.diagnostics.warnings.append((W_TERMSTATEBD,
1259 (last_t, xnames[offender_ix],
1260 X[offender_ix, last_ix],
1261 depdomains[offender_ix].get())))
1262
1263 variables = copyVarDict(self.variables)
1264
1265
1266 self.trajevents = {}
1267 for evix in range(len(self._eventNames)):
1268 evname = self._eventNames[evix]
1269 if Evpoints[evix] is None:
1270 self.trajevents[evname] = None
1271 else:
1272 try:
1273 ev_a_list = []
1274 for t in Evtimes[evix]:
1275 tix = find(alltData, t)
1276 ev_a_list.append(A[:,tix])
1277 ev_array = concatenate((Evpoints[evix],
1278 transpose(array(ev_a_list, 'd'))))
1279 del ev_a_list, tix
1280 except TypeError:
1281
1282 ev_array = Evpoints[evix]
1283 if last_ix >= 0 and last_ix < highest_ix:
1284
1285 last_ev_tix = npy.argmax(Evtimes[evix] >= alltData[last_ix])
1286 if last_ev_tix == 0 and Evtimes[evix][0] >= last_t:
1287
1288
1289 self.trajevents[evname] = None
1290 else:
1291
1292 ev_array = ev_array[:, :last_ev_tix+1]
1293 ev_times = Evtimes[evix][:last_ev_tix+1]
1294 self.trajevents[evname] = Pointset({'coordnames': xnames+anames,
1295 'indepvarname': 't',
1296 'coordarray': ev_array,
1297 'indepvararray': ev_times})
1298 else:
1299
1300 self.trajevents[evname] = Pointset({'coordnames': xnames+anames,
1301 'indepvarname': 't',
1302 'coordarray': ev_array,
1303 'indepvararray': Evtimes[evix]})
1304 if last_ix >= 0 and last_ix < highest_ix:
1305
1306 X = X[:, :last_ix]
1307 alltData = alltData[:last_ix]
1308 try:
1309 allxDataDict = dict(zip(xnames,X))
1310 except IndexError:
1311 print "Integration returned variable values of unexpected dimensions."
1312 print " Did you change the system and not refresh the C library" \
1313 + " using the forcelibrefresh() method?"
1314 raise
1315
1316 anames = self.funcspec.auxvars
1317 try:
1318 if anames != []:
1319 if last_ix < highest_ix:
1320 A = A[:, :last_ix]
1321 try:
1322 allaDataDict = dict(zip(anames,A))
1323 except TypeError:
1324 print "Internal error! Type of A: ", type(A)
1325 raise
1326 except IndexError:
1327 print "Integration returned auxiliary values of unexpected dimensions."
1328 print " Did you change the system and not refresh the C library" \
1329 + " using the forcelibrefresh() method?"
1330 raise
1331 if int(Err) == 1 or (int(Err) == 2 and termcount == 1):
1332
1333 if self.algparams['poly_interp']:
1334 rhsfn = self._solver.Rhs
1335
1336
1337 dxvals = zeros((len(alltData),self.dimension),float)
1338 for tix, tval in enumerate(alltData):
1339
1340
1341
1342
1343
1344 dxvals[tix] = rhsfn(tval, list(X[:,tix]), plist)[0]
1345 for xi, x in enumerate(xnames):
1346 if len(alltData) > 1:
1347 if self.algparams['poly_interp']:
1348 interp = PiecewisePolynomial(alltData,
1349 array([allxDataDict[x], dxvals[:,xi]]).T, 2)
1350 else:
1351 interp = interp1d(alltData, allxDataDict[x])
1352 variables[x] = Variable(interp, 't', x, x)
1353 else:
1354 raise PyDSTool_ValueError("Fewer than 2 data points computed")
1355 for a in anames:
1356 if len(alltData) > 1:
1357 variables[a] = Variable(interp1d(alltData,allaDataDict[a]),
1358 't', a, a)
1359 else:
1360 raise PyDSTool_ValueError("Fewer than 2 data points computed")
1361
1362
1363 self.defined = True
1364 return Trajectory(trajname, variables.values(),
1365 abseps=self._abseps, globalt0=self.globalt0,
1366 checklevel=self.checklevel,
1367 FScompatibleNames=self._FScompatibleNames,
1368 FScompatibleNamesInv=self._FScompatibleNamesInv,
1369 events=self.trajevents,
1370 modelNames=self.name,
1371 modelEventStructs=self.eventstruct)
1372 else:
1373 try:
1374 diagnost_info = self.diagnostics._errorcodes[int(Err)]
1375 except TypeError:
1376
1377 print "Error code: ", Err
1378 diagnost_info = self.diagnostics._errorcodes[0]
1379 if self._solver.verbose:
1380 info(self.diagnostics.outputStats, "Output statistics")
1381 self.defined = False
1382
1383 if (len(alltData) == self.algparams['max_pts'] or \
1384 self.diagnostics.outputStats['num_steps'] >= self.algparams['max_pts']) \
1385 and alltData[-1] < tend:
1386 print "max_pts algorithmic parameter too small: current " + \
1387 "value is %i"%self.algparams['max_pts']
1388
1389 if self.diagnostics.outputStats['last_time']-tbegin > 0:
1390 ms = str(int(round(self.algparams['max_pts'] / \
1391 (self.diagnostics.outputStats['last_time'] - \
1392 tbegin)*(tend-tbegin))))
1393 else:
1394 ms = 'Inf'
1395 print "(recommended value for this trajectory segment is " + \
1396 "estimated to be %s (saved in diagnostics.errors attribute))"%str(ms)
1397 diagnost_info += " -- recommended value is " + ms
1398 self.diagnostics.errors.append((E_COMPUTFAIL,
1399 (self._solver.lastTime, diagnost_info)))
1400 raise PyDSTool_ExistError("No trajectory created")
1401
1402
1403 - def Rhs(self, t, xdict, pdict=None, asarray=True):
1420
1421
1422 - def Jacobian(self, t, xdict, pdict=None, asarray=True):
1441
1442
1443 - def JacobianP(self, t, xdict, pdict=None, asarray=True):
1462
1463
1464 - def AuxVars(self, t, xdict, pdict=None, asarray=True):
1480
1481
1482 - def MassMatrix(self, t, xdict, pdict=None, asarray=True):
1501
1502
1504 if self._solver is None:
1505 x0 = sortedDictValues(filteredDict(self.initialconditions, self.funcspec.vars))
1506
1507 self._solver = radau("radau5"+self._vf_filename_ext,
1508 rhs=self.name, phaseDim=self.dimension,
1509 paramDim=self.numpars,
1510 nAux=len(self.funcspec.auxvars),
1511 nEvents=len(self._eventNames),
1512 nExtInputs=len(self.inputs),
1513 hasJac=self.haveJacobian(),
1514 hasJacP=self.haveJacobian_pars(),
1515 hasMass=self.haveMass(),
1516 extraSpace=self.algparams['extraspace'])
1517 try:
1518 genDB.register(self)
1519 except PyDSTool_KeyError:
1520 errstr = "Generator " + self.name + ": this vector field's " +\
1521 "DLL is already in use"
1522 raise RuntimeError(errstr)
1523 if pars is not None:
1524
1525 self._solver.setRunParams(
1526 ic=sortedDictValues(filteredDict(self.initialconditions,
1527 self.funcspec.vars)),
1528 params=pars['params'],
1529 t0=pars['t0'], tend=pars['tend'],
1530 gt0=self.globalt0,
1531 refine=0, specTimes=[])
1532
1575
1576
1580
1581
1582
1583
1584
1585 symbolMapDict = {'abs': 'fabs', 'sign': 'signum', 'mod': 'fmod'}
1586
1587
1588 theGenSpecHelper.add(Radau_ODEsystem, symbolMapDict, 'c')
1589