1# -*- coding: utf-8 -*- 2""" 3Core functionality for ODESys. 4 5Note that it is possible to use new custom ODE integrators with pyodesys by 6providing a module with two functions named ``integrate_adaptive`` and 7``integrate_predefined``. See the ``pyodesys.integrators`` module for examples. 8""" 9 10from __future__ import absolute_import, division, print_function 11 12import copy 13import os 14import warnings 15from collections import defaultdict 16 17import numpy as np 18 19from .plotting import plot_result, plot_phase_plane 20from .results import Result 21from .util import _ensure_4args, _default 22 23 24class RecoverableError(Exception): 25 pass 26 27 28class ODESys(object): 29 """ Object representing an ODE system. 30 31 ``ODESys`` provides unified interface to: 32 33 - scipy.integarte.ode 34 - pygslodeiv2 35 - pyodeint 36 - pycvodes 37 38 The numerical integration can be performed either in an :meth:`adaptive` 39 or :meth:`predefined` mode. Where locations to report the solution is 40 chosen by the stepper or the user respectively. For convenience in user 41 code one may use :meth:`integrate` which automatically chooses between 42 the two based on the length of ``xout`` provided by the user. 43 44 Parameters 45 ---------- 46 f : callback 47 first derivatives of dependent variables (y) with respect to 48 dependent variable (x). Signature is any of: 49 - ``rhs(x, y[:]) -> f[:]`` 50 - ``rhs(x, y[:], p[:]) -> f[:]`` 51 - ``rhs(x, y[:], p[:], backend=math) -> f[:]``. 52 jac : callback 53 Jacobian matrix (dfdy). Required for implicit methods. 54 Signature should be either of: 55 - ``jac(x, y[:]) -> J`` 56 - ``jac(x, y[:], p[:]) -J``. 57 If ``nnz < 0``, ``J`` should be a 2D array-like object if ``nnz < 0`` 58 corresponding to a dense or banded jacobian (see also ``band``). 59 If ``nnz >= 0``, ``J`` should be an instance of ``scipy.sparse.csc_matrix``. 60 dfdx : callback 61 Signature ``dfdx(x, y[:], p[:]) -> out[:]`` (used by e.g. GSL) 62 jtimes : callback 63 Jacobian-vector product (Jv). Signature is ```jtimes(x, y[:], v[:]) -> Jv[:]``` 64 This is supported only by ``cvode``. 65 first_step_cb : callback 66 Signature ``step1st(x, y[:], p[:]) -> dx0`` (pass first_step==0 to use). 67 This is available for ``cvode``, ``odeint`` & ``gsl``, but not for ``scipy``. 68 roots_cb : callback 69 Signature ``roots_cb(x, y[:], p[:]=(), backend=math) -> discr[:]``. 70 nroots : int 71 Length of return vector from ``roots_cb``. 72 band : tuple of 2 integers or None (default: None) 73 If jacobian is banded: number of sub- and super-diagonals 74 names : iterable of strings (default : None) 75 Names of variables, used for referencing dependent variables by name 76 and for labels in plots. 77 param_names : iterable of strings (default: None) 78 Names of the parameters, used for referencing parameters by name. 79 indep_name : str 80 Name of the independent variable 81 dep_by_name : bool 82 When ``True`` :meth:`integrate` expects a dictionary as input for y0. 83 par_by_name : bool 84 When ``True`` :meth:`integrate` expects a dictionary as input for params. 85 latex_names : iterable of strings (default : None) 86 Names of variables in LaTeX format (e.g. for labels in plots). 87 latex_param_names : iterable of strings (default : None) 88 Names of parameters in LaTeX format (e.g. for labels in plots). 89 latex_indep_name : str 90 LaTeX formatted name of independent variable. 91 taken_names : iterable of str 92 Names of dependent variables which are calculated in pre_processors 93 pre_processors : iterable of callables (optional) 94 signature: f(x1[:], y1[:], params1[:]) -> x2[:], y2[:], params2[:]. 95 When modifying: insert at beginning. 96 post_processors : iterable of callables (optional) 97 signature: f(x2[:], y2[:, :], params2[:]) -> x1[:], y1[:, :], 98 params1[:] 99 When modifying: insert at end. 100 append_iv : bool 101 See :attr:`append_iv`. 102 autonomous_interface : bool (optional) 103 If given, sets the :attr:`autonomous_interface` to indicate whether 104 the system appears autonomous or not upon call to :meth:`integrate`. 105 autonomous_exprs : bool 106 Describes whether the independent variable appears in the rhs expressions. 107 If set to ``True`` the underlying solver is allowed to shift the 108 independent variable during integration. 109 nnz : int (default : -1) 110 Maximum number of non-zero entries in sparse jacobian. When 111 non-negative, jacobian is assumed to be dense or banded. 112 113 Attributes 114 ---------- 115 f_cb : callback 116 For evaluating the vector of derivatives. 117 j_cb : callback or None 118 For evaluating the Jacobian matrix of f. 119 dfdx_cb : callback or None 120 For evaluating the second order derivatives. 121 jtimes_cb : callback or None 122 For evaluating the jacobian-vector product. 123 first_step_cb : callback or None 124 For calculating the first step based on x0, y0 & p. 125 roots_cb : callback 126 nroots : int 127 nnz : int 128 names : tuple of strings 129 param_names : tuple of strings 130 description : str 131 dep_by_name : bool 132 par_by_name : bool 133 latex_names : tuple of str 134 latex_param_names : tuple of str 135 pre_processors : iterable of callbacks 136 post_processors : iterable of callbacks 137 append_iv : bool 138 If ``True`` params[:] passed to :attr:`f_cb`, :attr:`jac_cb` will contain 139 initial values of y. Note that this happens after pre processors have been 140 applied. 141 autonomous_interface : bool or None 142 Indicates whether the system appears autonomous upon call to 143 :meth:`integrate`. ``None`` indicates that it is unknown. 144 145 Examples 146 -------- 147 >>> odesys = ODESys(lambda x, y, p: p[0]*x + p[1]*y[0]*y[0]) 148 >>> yout, info = odesys.predefined([1], [0, .2, .5], [2, 1]) 149 >>> print(info['success']) 150 True 151 152 153 Notes 154 ----- 155 Banded jacobians are supported by "scipy" and "cvode" integrators. 156 157 """ 158 159 def __init__(self, f, jac=None, dfdx=None, jtimes=None, first_step_cb=None, roots_cb=None, nroots=None, 160 band=None, names=(), param_names=(), indep_name=None, description=None, dep_by_name=False, 161 par_by_name=False, latex_names=(), latex_param_names=(), latex_indep_name=None, 162 taken_names=None, pre_processors=None, post_processors=None, append_iv=False, 163 autonomous_interface=None, to_arrays_callbacks=None, autonomous_exprs=None, 164 _indep_autonomous_key=None, numpy=None, nnz=-1, **kwargs): 165 self.f_cb = _ensure_4args(f) 166 self.j_cb = _ensure_4args(jac) if jac is not None else None 167 self.jtimes_cb = _ensure_4args(jtimes) if jtimes is not None else None 168 self.dfdx_cb = dfdx 169 self.first_step_cb = first_step_cb 170 self.roots_cb = roots_cb 171 self.nroots = nroots or 0 172 if band is not None: 173 if not band[0] >= 0 or not band[1] >= 0: 174 raise ValueError("bands needs to be > 0 if provided") 175 self.band = band 176 self.nnz = nnz 177 self.names = tuple(names or ()) 178 self.param_names = tuple(param_names or ()) 179 self.indep_name = indep_name 180 self.description = description 181 self.dep_by_name = dep_by_name 182 self.par_by_name = par_by_name 183 self.latex_names = tuple(latex_names or ()) 184 self.latex_param_names = tuple(latex_param_names or ()) 185 self.latex_indep_name = latex_indep_name 186 self.taken_names = tuple(taken_names or ()) 187 self.pre_processors = pre_processors or [] 188 self.post_processors = post_processors or [] 189 self.append_iv = append_iv 190 self.autonomous_exprs = autonomous_exprs 191 if hasattr(self, 'autonomous_interface'): 192 if autonomous_interface is not None and autonomous_interface != self.autonomous_interface: 193 raise ValueError("Got conflicting autonomous_interface infomation.") 194 else: 195 if (autonomous_interface is None and self.autonomous_exprs and 196 len(self.post_processors) == 0 and len(self.pre_processors) == 0): 197 self.autonomous_interface = True 198 else: 199 self.autonomous_interface = autonomous_interface 200 201 if self.autonomous_interface not in (True, False, None): 202 raise ValueError("autonomous_interface needs to be a boolean value or None.") 203 self._indep_autonomous_key = _indep_autonomous_key 204 self.to_arrays_callbacks = to_arrays_callbacks 205 self.numpy = numpy or np 206 if len(kwargs) > 0: 207 raise ValueError("Unknown kwargs: %s" % str(kwargs)) 208 209 @staticmethod 210 def _array_from_dict(d, keys, numpy=np): 211 vals = [d[k] for k in keys] 212 lens = [len(v) for v in vals if hasattr(v, '__len__') and getattr(v, 'ndim', 1) > 0] 213 if len(lens) == 0: 214 return vals, True 215 else: 216 if not all(l == lens[0] for l in lens): 217 raise ValueError("Mixed lenghts in dictionary.") 218 out = numpy.empty((lens[0], len(vals)), dtype=object) 219 for idx, v in enumerate(vals): 220 if getattr(v, 'ndim', -1) == 0: 221 for j in range(lens[0]): 222 out[j, idx] = v 223 else: 224 try: 225 for j in range(lens[0]): 226 out[j, idx] = v[j] 227 except TypeError: 228 out[:, idx] = v 229 return out, False 230 231 def _conditional_from_dict(self, cont, by_name, names): 232 if isinstance(cont, dict): 233 if not by_name: 234 raise ValueError("not by name, yet a dictionary was passed.") 235 cont, tp = self._array_from_dict(cont, names, numpy=self.numpy) 236 else: 237 tp = False 238 return cont, tp 239 240 def to_arrays(self, x, y, p, callbacks=None, reshape=True): 241 try: 242 nx = len(x) 243 except TypeError: 244 _x = 0*x, x 245 else: 246 _x = (0*x[0], x[0]) if nx == 0 else x 247 248 _names = [n for n in self.names if n not in self.taken_names] 249 _y, tp_y = self._conditional_from_dict(y, self.dep_by_name, _names) 250 _p, tp_p = self._conditional_from_dict(p, self.par_by_name, self.param_names) 251 del _names 252 253 callbacks = callbacks or self.to_arrays_callbacks 254 if callbacks is not None: # e.g. dedimensionalisation 255 if len(callbacks) != 3: 256 raise ValueError("Need 3 callbacks/None values.") 257 _x, _y, _p = [e if cb is None else cb(e) for cb, e in zip(callbacks, [_x, _y, _p])] 258 _y = self.numpy.atleast_1d(_y) 259 if self._indep_autonomous_key: 260 if _y.shape[-1] == self.ny: 261 pass 262 elif _y.shape[-1] == self.ny - 1: 263 _y = self.numpy.concatenate((_y, _x[0]*self.numpy.ones(_y.shape[:-1] + (1,))), axis=-1) 264 else: 265 raise ValueError("y of incorrect shape") 266 267 arrs = [arr.T if tp else arr for tp, arr in 268 zip([False, tp_y, tp_p], map(self.numpy.atleast_1d, (_x, _y, _p)))] 269 if reshape: 270 extra_shape = None 271 for a in arrs: 272 if a.ndim == 1: 273 continue 274 elif a.ndim == 2: 275 if extra_shape is None: 276 extra_shape = a.shape[0] 277 else: 278 if extra_shape != a.shape[0]: 279 raise ValueError("Size mismatch!") 280 else: 281 raise NotImplementedError("Only 2 dimensions currently supported.") 282 if extra_shape is not None: 283 arrs = [a if a.ndim == 2 else self.numpy.tile(a, (extra_shape, 1)) for a in arrs] 284 return arrs 285 286 def pre_process(self, xout, y0, params=()): 287 """ Transforms input to internal values, used internally. """ 288 for pre_processor in self.pre_processors: 289 xout, y0, params = pre_processor(xout, y0, params) 290 return [self.numpy.atleast_1d(arr) for arr in (xout, y0, params)] 291 292 def post_process(self, xout, yout, params): 293 """ Transforms internal values to output, used internally. """ 294 for post_processor in self.post_processors: 295 xout, yout, params = post_processor(xout, yout, params) 296 return xout, yout, params 297 298 def adaptive(self, y0, x0, xend, params=(), **kwargs): 299 """ Integrate with integrator chosen output. 300 301 Parameters 302 ---------- 303 integrator : str 304 See :meth:`integrate`. 305 y0 : array_like 306 See :meth:`integrate`. 307 x0 : float 308 Initial value of the independent variable. 309 xend : float 310 Final value of the independent variable. 311 params : array_like 312 See :meth:`integrate`. 313 \\*\\*kwargs : 314 See :meth:`integrate`. 315 316 Returns 317 ------- 318 Same as :meth:`integrate` 319 """ 320 return self.integrate((x0, xend), y0, 321 params=params, **kwargs) 322 323 def predefined(self, y0, xout, params=(), **kwargs): 324 """ Integrate with user chosen output. 325 326 Parameters 327 ---------- 328 integrator : str 329 See :meth:`integrate`. 330 y0 : array_like 331 See :meth:`integrate`. 332 xout : array_like 333 params : array_like 334 See :meth:`integrate`. 335 \\*\\*kwargs: 336 See :meth:`integrate` 337 338 Returns 339 ------- 340 Length 2 tuple : (yout, info) 341 See :meth:`integrate`. 342 """ 343 xout, yout, info = self.integrate(xout, y0, params=params, 344 force_predefined=True, **kwargs) 345 return yout, info 346 347 def integrate(self, x, y0, params=(), atol=1e-8, rtol=1e-8, **kwargs): 348 """ Integrate the system of ordinary differential equations. 349 350 Solves the initial value problem (IVP). 351 352 Parameters 353 ---------- 354 x : array_like or pair (start and final time) or float 355 if float: 356 make it a pair: (0, x) 357 if pair or length-2 array: 358 initial and final value of the independent variable 359 if array_like: 360 values of independent variable report at 361 y0 : array_like 362 Initial values at x[0] for the dependent variables. 363 params : array_like (default: tuple()) 364 Value of parameters passed to user-supplied callbacks. 365 integrator : str or None 366 Name of integrator, one of: 367 - 'scipy': :meth:`_integrate_scipy` 368 - 'gsl': :meth:`_integrate_gsl` 369 - 'odeint': :meth:`_integrate_odeint` 370 - 'cvode': :meth:`_integrate_cvode` 371 372 See respective method for more information. 373 If ``None``: ``os.environ.get('PYODESYS_INTEGRATOR', 'scipy')`` 374 atol : float 375 Absolute tolerance 376 rtol : float 377 Relative tolerance 378 with_jacobian : bool or None (default) 379 Whether to use the jacobian. When ``None`` the choice is 380 done automatically (only used when required). This matters 381 when jacobian is derived at runtime (high computational cost). 382 with_jtimes : bool (default: False) 383 Whether to use the jacobian-vector product. This is only supported 384 by ``cvode`` and only when ``linear_solver`` is one of: gmres', 385 'gmres_classic', 'bicgstab', 'tfqmr'. See the documentation 386 for ``pycvodes`` for more information. 387 force_predefined : bool (default: False) 388 override behaviour of ``len(x) == 2`` => :meth:`adaptive` 389 \\*\\*kwargs : 390 Additional keyword arguments for ``_integrate_$(integrator)``. 391 392 Returns 393 ------- 394 Length 3 tuple: (x, yout, info) 395 x : array of values of the independent variable 396 yout : array of the dependent variable(s) for the different 397 values of x. 398 info : dict ('nfev' is guaranteed to be a key) 399 """ 400 arrs = self.to_arrays(x, y0, params) 401 _x, _y, _p = _arrs = self.pre_process(*arrs) 402 ndims = [a.ndim for a in _arrs] 403 if ndims == [1, 1, 1]: 404 twodim = False 405 elif ndims == [2, 2, 2]: 406 twodim = True 407 else: 408 raise ValueError("Pre-processor made ndims inconsistent?") 409 410 if self.append_iv: 411 _p = self.numpy.concatenate((_p, _y), axis=-1) 412 413 if hasattr(self, 'ny'): 414 if _y.shape[-1] != self.ny: 415 raise ValueError("Incorrect shape of intern_y0") 416 if isinstance(atol, dict): 417 kwargs['atol'] = [atol[k] for k in self.names] 418 else: 419 kwargs['atol'] = atol 420 kwargs['rtol'] = rtol 421 422 integrator = kwargs.pop('integrator', None) 423 if integrator is None: 424 integrator = os.environ.get('PYODESYS_INTEGRATOR', 'scipy') 425 426 args = tuple(map(self.numpy.atleast_2d, (_x, _y, _p))) 427 428 self._current_integration_kwargs = kwargs 429 if isinstance(integrator, str): 430 nfo = getattr(self, '_integrate_' + integrator)(*args, **kwargs) 431 else: 432 kwargs['with_jacobian'] = getattr(integrator, 'with_jacobian', None) 433 nfo = self._integrate(integrator.integrate_adaptive, 434 integrator.integrate_predefined, 435 *args, **kwargs) 436 if twodim: 437 _xout = [d['internal_xout'] for d in nfo] 438 _yout = [d['internal_yout'] for d in nfo] 439 _params = [d['internal_params'] for d in nfo] 440 res = [Result(*(self.post_process(_xout[i], _yout[i], _params[i]) + (nfo[i], self))) 441 for i in range(len(nfo))] 442 else: 443 _xout = nfo[0]['internal_xout'] 444 _yout = nfo[0]['internal_yout'] 445 446 self._internal = _xout.copy(), _yout.copy(), _p.copy() 447 nfo = nfo[0] 448 res = Result(*(self.post_process(_xout, _yout, _p) + (nfo, self))) 449 return res 450 451 def chained_parameter_variation(self, *args, **kwargs): 452 """ See :func:`chained_parameter_variation`. """ 453 return chained_parameter_variation(self, *args, **kwargs) 454 455 def _integrate_scipy(self, intern_xout, intern_y0, intern_p, 456 atol=1e-8, rtol=1e-8, first_step=None, with_jacobian=None, 457 force_predefined=False, name=None, **kwargs): 458 """ Do not use directly (use ``integrate('scipy', ...)``). 459 460 Uses `scipy.integrate.ode <http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html>`_ 461 462 Parameters 463 ---------- 464 \\*args : 465 See :meth:`integrate`. 466 name : str (default: 'lsoda'/'dopri5' when jacobian is available/not) 467 What integrator wrapped in scipy.integrate.ode to use. 468 \\*\\*kwargs : 469 Keyword arguments passed onto `set_integrator(...) < 470 http://docs.scipy.org/doc/scipy/reference/generated/ 471 scipy.integrate.ode.set_integrator.html#scipy.integrate.ode.set_integrator>`_ 472 473 Returns 474 ------- 475 See :meth:`integrate`. 476 """ 477 from scipy.integrate import ode 478 ny = intern_y0.shape[-1] 479 nx = intern_xout.shape[-1] 480 results = [] 481 for _xout, _y0, _p in zip(intern_xout, intern_y0, intern_p): 482 if name is None: 483 if self.j_cb is None: 484 name = 'dopri5' 485 else: 486 name = 'lsoda' 487 if with_jacobian is None: 488 if name == 'lsoda': # lsoda might call jacobian 489 with_jacobian = True 490 elif name in ('dop853', 'dopri5'): 491 with_jacobian = False # explicit steppers 492 elif name == 'vode': 493 with_jacobian = kwargs.get('method', 'adams') == 'bdf' 494 495 def rhs(t, y, p=()): 496 rhs.ncall += 1 497 return self.f_cb(t, y, p) 498 rhs.ncall = 0 499 500 if self.j_cb is not None: 501 def jac(t, y, p=()): 502 jac.ncall += 1 503 return self.j_cb(t, y, p) 504 jac.ncall = 0 505 506 r = ode(rhs, jac=jac if with_jacobian else None) 507 if 'lband' in kwargs or 'uband' in kwargs or 'band' in kwargs: 508 raise ValueError("lband and uband set locally (set `band` at initialization instead)") 509 if self.band is not None: 510 kwargs['lband'], kwargs['uband'] = self.band 511 r.set_integrator(name, atol=atol, rtol=rtol, **kwargs) 512 if len(_p) > 0: 513 r.set_f_params(_p) 514 r.set_jac_params(_p) 515 r.set_initial_value(_y0, _xout[0]) 516 if nx == 2 and not force_predefined: 517 mode = 'adaptive' 518 if name in ('vode', 'lsoda'): 519 warnings.warn("'adaptive' mode with SciPy's integrator (vode/lsoda) may overshoot (itask=2)") 520 warnings.warn("'adaptive' mode with SciPy's integrator is unreliable, consider using e.g. cvode") 521 # vode itask 2 (may overshoot) 522 ysteps = [_y0] 523 xsteps = [_xout[0]] 524 while r.t < _xout[1]: 525 r.integrate(_xout[1], step=True) 526 if not r.successful(): 527 raise RuntimeError("failed") 528 xsteps.append(r.t) 529 ysteps.append(r.y) 530 else: 531 xsteps, ysteps = [], [] 532 533 def solout(x, y): 534 xsteps.append(x) 535 ysteps.append(y) 536 r.set_solout(solout) 537 r.integrate(_xout[1]) 538 if not r.successful(): 539 raise RuntimeError("failed") 540 _yout = np.array(ysteps) 541 _xout = np.array(xsteps) 542 543 else: # predefined 544 mode = 'predefined' 545 _yout = np.empty((nx, ny)) 546 _yout[0, :] = _y0 547 for idx in range(1, nx): 548 r.integrate(_xout[idx]) 549 if not r.successful(): 550 raise RuntimeError("failed") 551 _yout[idx, :] = r.y 552 info = { 553 'internal_xout': _xout, 554 'internal_yout': _yout, 555 'internal_params': _p, 556 'success': r.successful(), 557 'nfev': rhs.ncall, 558 'n_steps': -1, # don't know how to obtain this number 559 'name': name, 560 'mode': mode, 561 'atol': atol, 562 'rtol': rtol 563 } 564 if self.j_cb is not None: 565 info['njev'] = jac.ncall 566 results.append(info) 567 return results 568 569 def _integrate(self, adaptive, predefined, intern_xout, intern_y0, intern_p, 570 atol=1e-8, rtol=1e-8, first_step=0.0, with_jacobian=None, 571 force_predefined=False, **kwargs): 572 nx = intern_xout.shape[-1] 573 results = [] 574 for _xout, _y0, _p in zip(intern_xout, intern_y0, intern_p): 575 new_kwargs = dict(dx0=first_step, atol=atol, 576 rtol=rtol, check_indexing=False) 577 new_kwargs.update(kwargs) 578 579 def _f(x, y, fout): 580 try: 581 if len(_p) > 0: 582 fout[:] = np.asarray(self.f_cb(x, y, _p)) 583 else: 584 fout[:] = np.asarray(self.f_cb(x, y)) 585 except RecoverableError: 586 return 1 # recoverable error 587 588 if with_jacobian is None: 589 raise Exception("Must pass with_jacobian") 590 elif with_jacobian is True: 591 if self.nnz >= 0: 592 def _j(x, y, data, colptrs, rowvals): 593 if len(_p) > 0: 594 J = self.j_cb(x, y, _p) 595 else: 596 J = self.j_cb(x, y) 597 J = J.asformat("csc") 598 data[:] = J.data 599 colptrs[:] = J.indptr 600 rowvals[:] = J.indices 601 602 new_kwargs['nnz'] = self.nnz 603 else: 604 def _j(x, y, jout, dfdx_out=None, fy=None): 605 if len(_p) > 0: 606 jout[:, :] = np.asarray(self.j_cb(x, y, _p)) 607 else: 608 jout[:, :] = np.asarray(self.j_cb(x, y)) 609 if dfdx_out is not None: 610 if len(_p) > 0: 611 dfdx_out[:] = np.asarray(self.dfdx_cb(x, y, _p)) 612 else: 613 dfdx_out[:] = np.asarray(self.dfdx_cb(x, y)) 614 else: 615 _j = None 616 617 with_jtimes = kwargs.pop('with_jtimes', False) 618 if with_jtimes is True: 619 def _jtimes(v, Jv, x, y, fy=None): 620 yv = np.concatenate((y, v)) 621 if len(_p) > 0: 622 Jv[:] = np.asarray(self.jtimes_cb(x, yv, _p)) 623 else: 624 Jv[:] = np.asarray(self.jtimes_cb(x, yv)) 625 new_kwargs['jtimes'] = _jtimes 626 627 if self.first_step_cb is not None: 628 def _first_step(x, y): 629 if len(_p) > 0: 630 return self.first_step_cb(x, y, _p) 631 else: 632 return self.first_step_cb(x, y) 633 if 'dx0cb' in new_kwargs: 634 raise ValueError("cannot override dx0cb") 635 else: 636 new_kwargs['dx0cb'] = _first_step 637 638 if self.roots_cb is not None: 639 def _roots(x, y, out): 640 if len(_p) > 0: 641 out[:] = np.asarray(self.roots_cb(x, y, _p)) 642 else: 643 out[:] = np.asarray(self.roots_cb(x, y)) 644 if 'roots' in new_kwargs: 645 raise ValueError("cannot override roots") 646 else: 647 new_kwargs['roots'] = _roots 648 if 'nroots' in new_kwargs: 649 raise ValueError("cannot override nroots") 650 new_kwargs['nroots'] = self.nroots 651 652 if nx == 2 and not force_predefined: 653 _xout, yout, info = adaptive(_f, _j, _y0, *_xout, **new_kwargs) 654 info['mode'] = 'adaptive' 655 else: 656 yout, info = predefined(_f, _j, _y0, _xout, **new_kwargs) 657 info['mode'] = 'predefined' 658 659 info['internal_xout'] = _xout 660 info['internal_yout'] = yout 661 info['internal_params'] = _p 662 results.append(info) 663 return results 664 665 def _integrate_gsl(self, *args, **kwargs): 666 """ Do not use directly (use ``integrate(..., integrator='gsl')``). 667 668 Uses `GNU Scientific Library <http://www.gnu.org/software/gsl/>`_ 669 (via `pygslodeiv2 <https://pypi.python.org/pypi/pygslodeiv2>`_) 670 to integrate the ODE system. 671 672 Parameters 673 ---------- 674 \\*args : 675 see :meth:`integrate` 676 method : str (default: 'bsimp') 677 what stepper to use, see :py:attr:`gslodeiv2.steppers` 678 \\*\\*kwargs : 679 keyword arguments passed onto 680 :py:func:`gslodeiv2.integrate_adaptive`/:py:func:`gslodeiv2.integrate_predefined` 681 682 Returns 683 ------- 684 See :meth:`integrate` 685 """ 686 import pygslodeiv2 # Python interface GSL's "odeiv2" integrators 687 kwargs['with_jacobian'] = kwargs.get( 688 'method', 'bsimp') in pygslodeiv2.requires_jac 689 return self._integrate(pygslodeiv2.integrate_adaptive, 690 pygslodeiv2.integrate_predefined, 691 *args, **kwargs) 692 693 def _integrate_odeint(self, *args, **kwargs): 694 """ Do not use directly (use ``integrate(..., integrator='odeint')``). 695 696 Uses `Boost.Numeric.Odeint <http://www.odeint.com>`_ 697 (via `pyodeint <https://pypi.python.org/pypi/pyodeint>`_) to integrate 698 the ODE system. 699 """ 700 import pyodeint # Python interface to boost's odeint integrators 701 kwargs['with_jacobian'] = kwargs.get( 702 'method', 'rosenbrock4') in pyodeint.requires_jac 703 return self._integrate(pyodeint.integrate_adaptive, 704 pyodeint.integrate_predefined, 705 *args, **kwargs) 706 707 def _integrate_cvode(self, *args, **kwargs): 708 """ Do not use directly (use ``integrate(..., integrator='cvode')``). 709 710 Uses CVode from CVodes in 711 `SUNDIALS <https://computation.llnl.gov/casc/sundials/>`_ 712 (via `pycvodes <https://pypi.python.org/pypi/pycvodes>`_) 713 to integrate the ODE system. """ 714 import pycvodes # Python interface to SUNDIALS's cvodes integrators 715 kwargs['with_jacobian'] = kwargs.get('method', 'bdf') in pycvodes.requires_jac 716 if 'lband' in kwargs or 'uband' in kwargs or 'band' in kwargs: 717 raise ValueError("lband and uband set locally (set at" 718 " initialization instead)") 719 if self.band is not None: 720 kwargs['lband'], kwargs['uband'] = self.band 721 kwargs['autonomous_exprs'] = self.autonomous_exprs 722 723 return self._integrate(pycvodes.integrate_adaptive, 724 pycvodes.integrate_predefined, 725 *args, **kwargs) 726 727 def _plot(self, cb, internal_xout=None, internal_yout=None, 728 internal_params=None, **kwargs): 729 kwargs = kwargs.copy() 730 if 'x' in kwargs or 'y' in kwargs or 'params' in kwargs: 731 raise ValueError("x and y from internal_xout and internal_yout") 732 733 _internal = getattr(self, '_internal', [None]*3) 734 x, y, p = (_default(internal_xout, _internal[0]), 735 _default(internal_yout, _internal[1]), 736 _default(internal_params, _internal[2])) 737 for post_processor in self.post_processors: 738 x, y, p = post_processor(x, y, p) 739 740 if 'names' not in kwargs: 741 kwargs['names'] = getattr(self, 'names', None) 742 else: 743 if 'indices' not in kwargs and getattr(self, 'names', None) is not None: 744 kwargs['indices'] = [self.names.index(n) for n in kwargs['names']] 745 kwargs['names'] = self.names 746 return cb(x, y, **kwargs) 747 748 def plot_result(self, **kwargs): 749 """ Plots the integrated dependent variables from last integration. 750 751 This method will be deprecated. Please use :meth:`Result.plot`. 752 See :func:`pyodesys.plotting.plot_result` 753 """ 754 return self._plot(plot_result, **kwargs) 755 756 def plot_phase_plane(self, indices=None, **kwargs): 757 """ Plots a phase portrait from last integration. 758 759 This method will be deprecated. Please use :meth:`Result.plot_phase_plane`. 760 See :func:`pyodesys.plotting.plot_phase_plane` 761 """ 762 return self._plot(plot_phase_plane, indices=indices, **kwargs) 763 764 def _jac_eigenvals_svd(self, xval, yvals, intern_p): 765 from scipy.linalg import svd 766 J = self.j_cb(xval, yvals, intern_p) 767 return svd(J, compute_uv=False) 768 769 def stiffness(self, xyp=None, eigenvals_cb=None): 770 """ [DEPRECATED] Use :meth:`Result.stiffness`, stiffness ration 771 772 Running stiffness ratio from last integration. 773 Calculate sittness ratio, i.e. the ratio between the largest and 774 smallest absolute eigenvalue of the jacobian matrix. The user may 775 supply their own routine for calculating the eigenvalues, or they 776 will be calculated from the SVD (singular value decomposition). 777 Note that calculating the SVD for any but the smallest Jacobians may 778 prove to be prohibitively expensive. 779 780 Parameters 781 ---------- 782 xyp : length 3 tuple (default: None) 783 internal_xout, internal_yout, internal_params, taken 784 from last integration if not specified. 785 eigenvals_cb : callback (optional) 786 Signature (x, y, p) (internal variables), when not provided an 787 internal routine will use ``self.j_cb`` and ``scipy.linalg.svd``. 788 789 """ 790 if eigenvals_cb is None: 791 if self.band is not None: 792 raise NotImplementedError 793 eigenvals_cb = self._jac_eigenvals_svd 794 795 if xyp is None: 796 x, y, intern_p = self._internal 797 else: 798 x, y, intern_p = self.pre_process(*xyp) 799 800 singular_values = [] 801 for xval, yvals in zip(x, y): 802 singular_values.append(eigenvals_cb(xval, yvals, intern_p)) 803 804 return (np.abs(singular_values).max(axis=-1) / 805 np.abs(singular_values).min(axis=-1)) 806 807 808class OdeSys(ODESys): 809 """ DEPRECATED, use ODESys instead. """ 810 pass 811 812 813def _new_x(xout, x, guaranteed_autonomous): 814 if guaranteed_autonomous: 815 return 0, abs(x[-1] - xout[-1]) # rounding 816 else: 817 return xout[-1], x[-1] 818 819 820def integrate_auto_switch(odes, kw, x, y0, params=(), **kwargs): 821 """ Auto-switching between formulations of ODE system. 822 823 In case one has a formulation of a system of ODEs which is preferential in 824 the beginning of the integration, this function allows the user to run the 825 integration with this system where it takes a user-specified maximum number 826 of steps before switching to another formulation (unless final value of the 827 independent variables has been reached). Number of systems used i returned 828 as ``nsys`` in info dict. 829 830 Parameters 831 ---------- 832 odes : iterable of :class:`OdeSy` instances 833 kw : dict mapping kwarg to iterables of same legnth as ``odes`` 834 x : array_like 835 y0 : array_like 836 params : array_like 837 \\*\\*kwargs: 838 See :meth:`ODESys.integrate` 839 840 Notes 841 ----- 842 Plays particularly well with :class:`symbolic.TransformedSys`. 843 844 """ 845 x_arr = np.asarray(x) 846 if x_arr.shape[-1] > 2: 847 raise NotImplementedError("Only adaptive support return_on_error for now") 848 multimode = False if x_arr.ndim < 2 else x_arr.shape[0] 849 nfo_keys = ('nfev', 'njev', 'time_cpu', 'time_wall') 850 851 next_autonomous = getattr(odes[0], 'autonomous_interface', False) == True # noqa (np.True_) 852 if multimode: 853 tot_x = [np.array([0] if next_autonomous else [x[_][0]]) for _ in range(multimode)] 854 tot_y = [np.asarray([y0[_]]) for _ in range(multimode)] 855 tot_nfo = [defaultdict(int) for _ in range(multimode)] 856 glob_x = [_[0] for _ in x] if next_autonomous else [0.0]*multimode 857 else: 858 tot_x, tot_y, tot_nfo = np.array([0 if next_autonomous else x[0]]), np.asarray([y0]), defaultdict(int) 859 glob_x = x[0] if next_autonomous else 0.0 860 861 for oi in range(len(odes)): 862 if oi < len(odes) - 1: 863 next_autonomous = getattr(odes[oi+1], 'autonomous_interface', False) == True # noqa (np.True_) 864 _int_kw = kwargs.copy() 865 for k, v in kw.items(): 866 _int_kw[k] = v[oi] 867 res = odes[oi].integrate(x, y0, params, **_int_kw) 868 869 if multimode: 870 for idx in range(multimode): 871 tot_x[idx] = np.concatenate((tot_x[idx], res[idx].xout[1:] + glob_x[idx])) 872 tot_y[idx] = np.concatenate((tot_y[idx], res[idx].yout[1:, :])) 873 for k in nfo_keys: 874 if k in res[idx].info: 875 tot_nfo[idx][k] += res[idx].info[k] 876 tot_nfo[idx]['success'] = res[idx].info['success'] 877 else: 878 tot_x = np.concatenate((tot_x, res.xout[1:] + glob_x)) 879 tot_y = np.concatenate((tot_y, res.yout[1:, :])) 880 for k in nfo_keys: 881 if k in res.info: 882 tot_nfo[k] += res.info[k] 883 tot_nfo['success'] = res.info['success'] 884 885 if multimode: 886 if all([r.info['success'] for r in res]): 887 break 888 else: 889 if res.info['success']: 890 break 891 if oi < len(odes) - 1: 892 if multimode: 893 _x, y0 = [], [] 894 for idx in range(multimode): 895 _x.append(_new_x(res[idx].xout, x[idx], next_autonomous)) 896 y0.append(res[idx].yout[-1, :]) 897 if next_autonomous: 898 glob_x[idx] += res[idx].xout[-1] 899 x = _x 900 else: 901 x = _new_x(res.xout, x, next_autonomous) 902 y0 = res.yout[-1, :] 903 if next_autonomous: 904 glob_x += res.xout[-1] 905 if multimode: # don't return defaultdict 906 tot_nfo = [dict(nsys=oi+1, **_nfo) for _nfo in tot_nfo] 907 return [Result(tot_x[idx], tot_y[idx], res[idx].params, tot_nfo[idx], odes[0]) 908 for idx in range(len(res))] 909 else: 910 tot_nfo = dict(nsys=oi+1, **tot_nfo) 911 return Result(tot_x, tot_y, res.params, tot_nfo, odes[0]) 912 913 914integrate_chained = integrate_auto_switch # deprecated name 915 916 917def chained_parameter_variation(subject, durations, y0, varied_params, default_params=None, 918 integrate_kwargs=None, x0=None, npoints=1, numpy=None): 919 """ Integrate an ODE-system for a serie of durations with some parameters changed in-between 920 921 Parameters 922 ---------- 923 subject : function or ODESys instance 924 If a function: should have the signature of :meth:`pyodesys.ODESys.integrate` 925 (and resturn a :class:`pyodesys.results.Result` object). 926 If a ODESys instance: the ``integrate`` method will be used. 927 durations : iterable of floats 928 Spans of the independent variable. 929 y0 : dict or array_like 930 varied_params : dict mapping parameter name (or index) to array_like 931 Each array_like need to be of same length as durations. 932 default_params : dict or array_like 933 Default values for the parameters of the ODE system. 934 integrate_kwargs : dict 935 Keyword arguments passed on to ``integrate``. 936 x0 : float-like 937 First value of independent variable. default: 0. 938 npoints : int 939 Number of points per sub-interval. 940 941 Examples 942 -------- 943 >>> odesys = ODESys(lambda t, y, p: [-p[0]*y[0]]) 944 >>> int_kw = dict(integrator='cvode', method='adams', atol=1e-12, rtol=1e-12) 945 >>> kwargs = dict(default_params=[0], integrate_kwargs=int_kw) 946 >>> res = chained_parameter_variation(odesys, [2, 3], [42], {0: [.7, .1]}, **kwargs) 947 >>> mask1 = res.xout <= 2 948 >>> import numpy as np 949 >>> np.allclose(res.yout[mask1, 0], 42*np.exp(-.7*res.xout[mask1])) 950 True 951 >>> mask2 = 2 <= res.xout 952 >>> np.allclose(res.yout[mask2, 0], res.yout[mask2, 0][0]*np.exp(-.1*(res.xout[mask2] - res.xout[mask2][0]))) 953 True 954 955 """ 956 assert len(durations) > 0, 'need at least 1 duration (preferably many)' 957 assert npoints > 0, 'need at least 1 point per duration' 958 for k, v in varied_params.items(): 959 if len(v) != len(durations): 960 raise ValueError("Mismathced lengths of durations and varied_params") 961 962 if isinstance(subject, ODESys): 963 integrate = subject.integrate 964 numpy = numpy or subject.numpy 965 else: 966 integrate = subject 967 numpy = numpy or np 968 969 default_params = default_params or {} 970 integrate_kwargs = integrate_kwargs or {} 971 972 def _get_idx(cont, idx): 973 if isinstance(cont, dict): 974 return {k: (v[idx] if hasattr(v, '__len__') and getattr(v, 'ndim', 1) > 0 else v) 975 for k, v in cont.items()} 976 else: 977 return cont[idx] 978 979 durations = numpy.cumsum(durations) 980 for idx_dur in range(len(durations)): 981 params = copy.copy(default_params) 982 for k, v in varied_params.items(): 983 params[k] = v[idx_dur] 984 if idx_dur == 0: 985 if x0 is None: 986 x0 = durations[0]*0 987 out = integrate(numpy.linspace(x0, durations[0], npoints + 1), y0, params, **integrate_kwargs) 988 else: 989 if isinstance(out, Result): 990 out.extend_by_integration(durations[idx_dur], params, npoints=npoints, **integrate_kwargs) 991 else: 992 for idx_res, r in enumerate(out): 993 r.extend_by_integration(durations[idx_dur], _get_idx(params, idx_res), 994 npoints=npoints, **integrate_kwargs) 995 996 return out 997