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

Source Code for Module PyDSTool.Generator.Dopri_ODEsystem'

   1  # Dopri ODE system
 
   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  # Other imports
 
  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  #import distutils
 
  22  from numpy.distutils.core import setup, Extension 
  23  from distutils.sysconfig import get_python_inc 
  24  from time import clock, sleep 
  25  
 
  26  # path to the installation
 
  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  
 
35 -class dopri(integrator):
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 # check module's directory 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 # For a run, we want to ensure indices are set to 0 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 # For a continue, we do not set indices to 0 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 # interpret all greater values as 1 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
282 -class Dopri_ODEsystem(ODEsystem):
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
304 - def __init__(self, kw):
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 # delete because not covered in ODEsystem 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': [], # for events 346 'maxevtpts': 1000, # for events 347 'eventInt': [], # set using setEventInterval only 348 'eventDelay': [], # set using setEventDelay only 349 'eventTol': [], # set using setEventTol only 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 # verify that no additional keys are present in algparams, after 365 # defaults are added above 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
408 - def forceLibRefresh(self):
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 # (try to) free dopri module from namespace 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 ## del(self._integMod) 421 except NameError: 422 # modules weren't loaded, so nothing to do 423 delfiles = False 424 if delfiles: 425 gc.collect() 426 # still not able to delete these files!!!!! Argh! 427 ## if os.path.isfile(os.path.join(os.getcwd(), 428 ## "dop853"+self._vf_filename_ext+".py")): 429 ## os.remove(os.path.join(os.getcwd(), 430 ## "dop853"+self._vf_filename_ext+".py")) 431 ## if os.path.isfile(os.path.join(os.getcwd(), 432 ## "_dop853"+self._vf_filename_ext+self._dllext)): 433 ## os.remove(os.path.join(os.getcwd(), 434 ## "_dop853"+self._vf_filename_ext+self._dllext)) 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 ## self.makeLibSource() 439 440
441 - def _prepareEventSpecs(self):
442 eventActive = [] 443 eventTerm = [] 444 eventDir = [] 445 eventDelay = [] 446 eventTol = [] 447 maxbisect = [] 448 eventInt = [] 449 # convert event specs (term, active, etc.) into integparam specs 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 # if event 'precise' flags set to False then set their tolerances 456 # to be > max_step 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
488 - def makeLibSource(self, include=[]):
489 """makeLibSource generates the C source for the vector field specification. 490 It should be called only once per vector field.""" 491 492 # Make vector field (and event) file for compilation 493 assert isinstance(include, list), "includes must be in form of a list" 494 # codes for library types (default is USERLIB, since compiler will look in standard library d 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 # don't repeat libraries 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 # sorted version of var, par, and input names 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 # add to defines 558 pardefines += self.funcspec._defstr+" "+p+"\tp_["+str(i)+"]\n" 559 for i in xrange(self.dimension): 560 v = vnames[i] 561 # add to defines 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 # add to defines 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 # preprocess event code 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 # check fbody for calls to user-defined aux fns 588 # and add hidden p argument 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 # convert 'initcond(x)' to 'initcond("x")' for 595 # compatibility with C syntax 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 # add signature for auxiliary functions 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 # subs _p into auxfn-to-auxfn calls (but not to the signature) 633 fbody_parsed = addArgToCalls(fbody, 634 self.funcspec.auxfns.keys(), 635 "p_, wk_, xv_", notFirst=fname) 636 if 'initcond' in self.funcspec.auxfns: 637 # convert 'initcond(x)' to 'initcond("x")' for 638 # compatibility with C syntax, but don't affect the 639 # function signature! 640 fbody_parsed = wrapArgInCall(fbody_parsed, 641 'initcond', '"', notFirst=True) 642 allfilestr += "\n" + fbody_parsed + "\n\n" 643 # add auxiliary variables (shell of the function always present) 644 # add event functions 645 allfilestr += self.funcspec.auxspec[0] + allevs 646 # if jacobians or mass matrix not present, fill in dummy 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 """ #+ "\n/* Variable and parameter substitutions undefined:*/\n" + parundefines + varundefines + "\n" 665 # write out C file 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
676 - def compileLib(self, libsources=[], libdirs=[]):
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 # then DLL file already exists and we can't overwrite it at this 687 # time 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 # The following if statement attempts to avoid recompiling the SWIG wrapper 731 # if the files mentioned already exist, because in principle the SWIG interface 732 # only needs compiling once. But this step doesn't seem to work yet. 733 # Instead, it seems that SWIG always gets recompiled with everything else 734 # (at least on Win32). Maybe the list of files is incorrect... 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=.', #+os.getcwd(), # '-t/', 747 '-tdopri853_temp',#+self._compilation_tempdir, 748 '--build-base=dopri853_temp', '--build-purelib=dopri853_temp']#+self._compilation_sourcedir] 749 #script_args = ['-q', 'build', '--build-lib='+os.getcwd(), '-t/', 750 # '--build-base='+self._compilation_tempdir] 751 if self._compiler != '': 752 script_args.append('-c'+str(self._compiler)) 753 754 # include directories for libraries 755 narraydir = npy.get_numarray_include() 756 npydir = npy.get_include() 757 758 incdirs = [npydir, narraydir, os.getcwd(), self._compilation_sourcedir] #_compilation_tempdir] 759 incdirs.extend(libdirs) 760 # Use distutils to perform the compilation of the selected files 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 # library_dirs=['./'], 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 # Attempt to unload module through a shutdown() function being 776 # added to the SWIG module file. But it didn't work! 777 ## try: 778 ## modfilepy = open(os.path.join(self._compilation_tempdir, 779 ## "dop853"+self._vf_filename_ext+".py"), 'a') 780 ## extfilename = "_dop853"+self._vf_filename_ext 781 ## modfilepy.write("""# The following addition made by PyDSTool: 782 ##def shutdown(): 783 ## import sys 784 ## del sys.modules['""" + extfilename + """'] 785 ## """) 786 #### del """ + extfilename + """ 787 #### del new_doubleArray 788 #### del delete_doubleArray 789 #### del doubleArray_getitem 790 #### del doubleArray_setitem 791 #### del new_intArray 792 #### del delete_intArray 793 #### del intArray_getitem 794 #### del intArray_setitem 795 #### del Integrate 796 ## modfilepy.close() 797 ## except IOError: 798 ## print "dop853.py modifying error in dopri853 temp compilation " \ 799 ## + "directory" 800 ## raise 801 rout.start() # redirect stdout 802 try: 803 # move library files into the user's CWD 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 # temporary hack to fix numpy_distutils bug 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() # restore stdout
822 823 824 # def _ensureLoaded(self, modname): 825 ## if modname in sys.modules: 826 ## _integMod = reload(sys.modules[modname]) 827 ## else: 828 # try: 829 # _integMod = __import__(modname, globals()) 830 # except: 831 # print "Error in importing compiled vector field and integrator." 832 # print "Did you compile the RHS C code?" 833 # raise 834 # # Initialize integrator 835 # assert 'Integrate' in dir(_integMod), \ 836 # "dopri853 library does not contain Integrate()" 837 # return _integMod 838 839
840 - def compute(self, trajname, dirn='f', ics=None):
841 continue_integ = ODEsystem.prepDirection(self, dirn) 842 # _integMod = __import__("dop853"+self._vf_filename_ext, globals()) 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 # Check i.c.'s are well defined (finite) 862 self.checkInitialConditions() 863 self.setEventICs(self.initialconditions, self.globalt0) 864 # update event params in case changed since last run 865 self._prepareEventSpecs() 866 # Main integration 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 # _integMod = self._ensureLoaded("dop853"+self._vf_filename_ext) 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 # dopri does reverse time integration simply by switching t0 and t1 895 # and using negative steps 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 = [[],[]] # lower, then upper 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 # overwrite t0 from self.indepvariable.domain, but use its t1 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 # if self.inputs and self._extInputsChanged: 921 # self._extInputsChanged = False 922 # self._solver.setContParams(tend, plist, 923 # use_special, 924 # self.algparams['verbose'], 925 # True, deepcopy(self._inputVarList), 926 # deepcopy(self._inputTimeList)) 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 # self._extInputsChanged if global t0 changed so that can 940 # adjust times given to the integrator (it is blind to global t0 941 # when accesses input variable times) 942 self._ensure_inputs(self._extInputsChanged) 943 # hinit only set if not continue_integ 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 # record needed for bounds checking and truncation 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 # reset this now 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 # reverse the array object (no reverse method!) 1008 alltData = alltData[::-1] 1009 X = X[:,::-1] 1010 if anames != []: 1011 A = A[:,::-1] 1012 xnames = self._var_ixmap 1013 # Package up computed trajectory in Variable variables 1014 # Add external inputs warnings to self.diagnostics.warnings, if any 1015 # (not presently supported) 1016 ## for f in inputVarList: 1017 ## for winfo in f.diagnostics.warnings: 1018 ## self.diagnostics.warnings.append((W_NONTERMSTATEBD, 1019 ## (winfo[0], f.name, winfo[1], 1020 ## f.depdomain))) 1021 eventslist = self.eventstruct.query(['lowlevel', 'active']) 1022 termevents = self.eventstruct.query(['term'], eventslist) 1023 if self._eventNames != []: 1024 # build self.diagnostics.warnings because events happened -- 1025 # and keep a record of which times terminal events happened because 1026 # Model.py's event handling procedure assumes multiple events 1027 # happening at one time are listed in one warning 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 # For safety, we should assert that this event 1042 # also appears in termevents, but we don't 1043 if Evtimes[evix][0] in termevtimes.keys(): 1044 # append event name to this warning 1045 warning_ix = termevtimes[Evtimes[evix][0]] 1046 self.diagnostics.warnings[warning_ix][1][1].append(evname) 1047 else: 1048 # make new termevtime entry for the new warning 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 # append event name to this warning 1057 warning_ix = nontermevtimes[Evtimes[evix][ev]] 1058 self.diagnostics.warnings[warning_ix][1][1].append(evname) 1059 else: 1060 # make new nontermevtime entry for the new warning 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 # post-process check of variable bounds (if defined and algparams['checkBounds'] True) 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 # temp storage for repeatedly used object attributes (for lookup efficiency) 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 # no point in checking when the bounds are +/- infinity 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 # won't count as truncating unless the following checks 1100 # hold 1101 last_ix = next_last_ix 1102 offender_ix = xi 1103 if not isfinite(last_ix) and last_ix < 0: 1104 # only use +Inf hereon to flag no truncation needed 1105 last_ix = Inf 1106 elif last_ix >= 0 and last_ix < highest_ix: 1107 # truncate data 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 # Create variables (self.variables contains no actual data) 1115 variables = copyVarDict(self.variables) 1116 # build event pointset information (reset previous trajectory's) 1117 # don't include events after any truncation due to state bound violation 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 # A is empty 1134 ev_array = Evpoints[evix] 1135 if last_ix >= 0 and last_ix < highest_ix: 1136 # don't count last_ix = -1 which is the same as highest_ix 1137 last_ev_tix = npy.argmax(Evtimes[evix] >= alltData[last_ix]) 1138 if last_ev_tix == 0 and Evtimes[evix][0] >= last_t: 1139 # checks that there was actually a violation 1140 # - so no events to record 1141 self.trajevents[evname] = None 1142 else: 1143 # truncation needed 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 # no truncation needed 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 # truncate 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 # storage of all auxiliary variable data 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 # output OK 1184 if self.algparams['poly_interp']: 1185 rhsfn = self._solver.Rhs 1186 # when Dopri can output the Rhs values alongside variable 1187 # values then this won't be necessary 1188 dxvals = zeros((len(alltData),self.dimension),float) 1189 for tix, tval in enumerate(alltData): 1190 # solver's Rhs function already contains the inputs so no 1191 # need to recompute and provide here. 1192 #i = _pollInputs(sortedDictValues(self.inputs), tval, 1193 # self.checklevel) 1194 # X is the output variable array, but rhsfn demands a list 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 # final checks 1215 #self.validateSpec() 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 # errcode messed up from Dopri 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 # Did the solver run out of memory? 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 # avstep = (self.algparams['init_step']+self.diagnostics.outputStats['last_step'])/2. 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):
1256 """asarray is an unused, dummy argument for compatibility with Model.Rhs""" 1257 # must convert names to FS-compatible as '.' sorts before letters 1258 # while '_' sorts after! 1259 x = sortedDictValues(filteredDict(self._FScompatibleNames(xdict), 1260 self.funcspec.vars)) 1261 if pdict is None: 1262 pdict = self.pars 1263 # internal self.pars already is FS-compatible 1264 p = sortedDictValues(pdict) 1265 else: 1266 p = sortedDictValues(self._FScompatibleNames(pdict)) 1267 i = _pollInputs(sortedDictValues(self.inputs), 1268 t, self.checklevel) 1269 self._ensure_solver({'params': p, 't0': 0, 'tend': 1}) 1270 self._ensure_inputs() 1271 return self._solver.Rhs(t, x, p+i)[0]
1272 1273
1274 - def Jacobian(self, t, xdict, pdict=None, asarray=True):
1275 """asarray is an unused, dummy argument for compatibility with 1276 Model.Jacobian""" 1277 if self.haveJacobian(): 1278 x = sortedDictValues(filteredDict(self._FScompatibleNames(xdict), 1279 self.funcspec.vars)) 1280 if pdict is None: 1281 pdict = self.pars 1282 # internal self.pars already is FS-compatible 1283 p = sortedDictValues(pdict) 1284 else: 1285 p = sortedDictValues(self._FScompatibleNames(pdict)) 1286 i = _pollInputs(sortedDictValues(self.inputs), 1287 t, self.checklevel) 1288 self._ensure_solver({'params': p, 't0': 0, 'tend': 1}) 1289 self._ensure_inputs() 1290 return self._solver.Jacobian(t, x, p+i)[0] 1291 else: 1292 raise PyDSTool_ExistError("Jacobian not defined")
1293 1294
1295 - def JacobianP(self, t, xdict, pdict=None, asarray=True):
1296 """asarray is an unused, dummy argument for compatibility with 1297 Model.JacobianP""" 1298 if self.haveJacobian_pars(): 1299 x = sortedDictValues(filteredDict(self._FScompatibleNames(xdict), 1300 self.funcspec.vars)) 1301 if pdict is None: 1302 pdict = self.pars 1303 # internal self.pars already is FS-compatible 1304 p = sortedDictValues(pdict) 1305 else: 1306 p = sortedDictValues(self._FScompatibleNames(pdict)) 1307 i = _pollInputs(sortedDictValues(self.inputs), 1308 t, self.checklevel) 1309 self._ensure_solver({'params': p, 't0': 0, 'tend': 1}) 1310 self._ensure_inputs() 1311 return self._solver.JacobianP(t, x, p+i)[0] 1312 else: 1313 raise PyDSTool_ExistError("Jacobian w.r.t. parameters not defined")
1314 1315
1316 - def AuxVars(self, t, xdict, pdict=None, asarray=True):
1317 """asarray is an unused, dummy argument for compatibility with 1318 Model.AuxVars""" 1319 x = sortedDictValues(filteredDict(self._FScompatibleNames(xdict), 1320 self.funcspec.vars)) 1321 if pdict is None: 1322 pdict = self.pars 1323 # internal self.pars already is FS-compatible 1324 p = sortedDictValues(pdict) 1325 else: 1326 p = sortedDictValues(self._FScompatibleNames(pdict)) 1327 i = _pollInputs(sortedDictValues(self.inputs), 1328 t, self.checklevel) 1329 self._ensure_solver({'params': p, 't0': 0, 'tend': 1}) 1330 self._ensure_inputs() 1331 return self._solver.AuxFunc(t, x, p+i)[0]
1332 1333
1334 - def _ensure_solver(self, pars=None):
1335 if self._solver is None: 1336 sortedDictValues(filteredDict(self.initialconditions, self.funcspec.vars)) 1337 # _integMod = self._ensureLoaded("dopri853"+self._vf_filename_ext) 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 # tend value doesn't matter 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
1364 - def _ensure_inputs(self, force=False):
1365 if not self.inputs: 1366 return 1367 if force: 1368 listOK = False 1369 else: 1370 try: 1371 listOK = self._inputTimest0 == self.globalt0 1372 except AttributeError: 1373 # not yet defined, so proceed 1374 listOK = False 1375 if not listOK: 1376 self._inputVarList = [] 1377 self._inputTimeList = [] 1378 self._inputTimest0 = self.globalt0 1379 # inputVarList is a list of Variables or Pointsets 1380 for inp in sortedDictValues(self.inputs): 1381 if isinstance(inp, Variable): 1382 pts = inp.getDataPoints() 1383 if pts is None: 1384 raise TypeError("Can only pass external input Variable objects if based on" 1385 " an underlying mesh") 1386 else: 1387 tvals = copy(pts[inp.indepvarname]) 1388 tvals -= self.globalt0 1389 self._inputVarList.append(pts[inp.coordname].tolist()) 1390 self._inputTimeList.append(tvals.tolist()) 1391 elif isinstance(inp, Pointset): 1392 tvals = copy(inp.indepvararray) 1393 tvals -= self.globalt0 1394 self._inputVarList.append(inp[inp.coordname].tolist()) 1395 self._inputTimeList.append(tvals.tolist()) 1396 else: 1397 raise TypeError("Invalid type of input") 1398 if not self._solver.initExtInputs: 1399 self._solver.setExtInputs(True, deepcopy(self._inputVarList), 1400 deepcopy(self._inputTimeList)) 1401 elif not listOK: 1402 self._solver.clearExtInputs() 1403 self._solver.setExtInputs(True, deepcopy(self._inputVarList), 1404 deepcopy(self._inputTimeList)) 1405 self._solver.canContinue=True
1406 1407
1408 - def __del__(self):
1409 genDB.unregister(self) 1410 ODEsystem.__del__(self)
1411 1412 1413 1414 # Register this Generator with the database 1415 1416 symbolMapDict = {'abs': 'fabs', 'sign': 'signum', 'mod': 'fmod'} 1417 # in future, provide appropriate mappings for libraries math, 1418 # random, etc. (for now it's left to FuncSpec) 1419 theGenSpecHelper.add(Dopri_ODEsystem, symbolMapDict, 'c') 1420