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