1# -*- coding: utf-8 -*- 2 3from __future__ import print_function, absolute_import, division 4 5from collections import defaultdict, OrderedDict 6from itertools import product 7import math 8 9try: 10 from time import clock as perf_counter 11except ImportError: 12 from time import perf_counter 13 14import numpy as np 15import pytest 16 17try: 18 import sym 19except ImportError: 20 sym = None 21 sym_backends = [] 22else: 23 sym_backends = sym.Backend.backends.keys() 24 25from .. import ODESys 26from ..core import integrate_auto_switch, chained_parameter_variation 27from ..symbolic import SymbolicSys, ScaledSys, symmetricsys, PartiallySolvedSystem, get_logexp, _group_invariants 28from ..util import requires, pycvodes_double, pycvodes_klu 29from .bateman import bateman_full # analytic, never mind the details 30from .test_core import vdp_f 31from . import _cetsa 32 33 34def identity(x): 35 return x 36 37idty2 = (identity, identity) 38 39 40def _decay3(x, y, p): 41 return [ 42 -p[0]*y[0], 43 p[0]*y[0] - p[1]*y[1], 44 p[1]*y[1] - p[2]*y[2] 45 ] 46 47 48def _get_decay3(**kwargs): 49 return SymbolicSys.from_callback(_decay3, 3, 3, **kwargs) 50 51 52def _get_decay3_names(yn, pn, **kwargs): 53 def f(x, y, p): 54 y = [y[n] for n in yn] 55 p = [p[n] for n in pn] 56 return dict(zip(yn, [ 57 -p[0]*y[0], 58 p[0]*y[0] - p[1]*y[1], 59 p[1]*y[1] - p[2]*y[2] 60 ])) 61 return SymbolicSys.from_callback(f, names=yn, param_names=pn, dep_by_name=True, 62 par_by_name=True, indep_name='t', **kwargs) 63 64 65@requires('sym', 'scipy') 66def test_SymbolicSys(): 67 from pyodesys.integrators import RK4_example_integrator 68 odesys = SymbolicSys.from_callback(lambda x, y, p, be: [-y[0], y[0]], 2, 69 names=['foo', 'bar']) 70 assert odesys.autonomous_interface is True 71 assert isinstance(odesys.exprs, tuple) 72 with pytest.raises(ValueError): 73 odesys.integrate(1, [0]) 74 75 odesys2 = SymbolicSys.from_callback(lambda x, y, p, be: {'foo': -y['foo'], 76 'bar': y['foo']}, 2, 77 names=['foo', 'bar'], dep_by_name=True) 78 for system, y0 in zip([odesys, odesys2], [[2, 3], {'foo': 2, 'bar': 3}]): 79 xout, yout, info = system.integrate(1, y0, integrator=RK4_example_integrator, first_step=1e-3) 80 assert np.allclose(yout[:, 0], 2*np.exp(-xout)) 81 assert np.allclose(yout[:, 1], 3 + 2*(1 - np.exp(-xout))) 82 83 with pytest.raises(ValueError): 84 SymbolicSys.from_callback(lambda x, y, p, be: None, 2, names=['foo', 'bar']) 85 86 with pytest.raises(ValueError): 87 SymbolicSys.from_callback(lambda x, y, p, be: [], 2, names=['foo', 'bar']) 88 89 90@requires('sym', 'scipy') 91def test_SymbolicSys__indep_name(): 92 odesys = SymbolicSys.from_callback( 93 lambda t, y, p: { 94 'x': -p['a']*y['x'], 95 'y': -p['b']*y['y'] + p['a']*y['x'], 96 'z': p['b']*y['y'] 97 }, names='xyz', param_names='ab', dep_by_name=True, par_by_name=True) 98 pars = {'a': [11, 17, 19], 'b': 13} 99 results = odesys.integrate([42, 43, 44], {'x': 7, 'y': 5, 'z': 3}, pars) 100 for r, a in zip(results, pars['a']): 101 assert np.allclose(r.named_dep('x'), 7*np.exp(-a*(r.xout - r.xout[0]))) 102 103 104@requires('sym') 105def test_SymbolicSys__init_indep__init_dep(): 106 odesys = SymbolicSys.from_callback(lambda x, y, p, be: [-y[0], y[0]], 2, names=['foo', 'bar'], indep_name='t', 107 init_indep=True, init_dep=True) 108 assert odesys.init_indep.name == 'i_t' 109 assert [dep.name for dep in odesys.init_dep] == ['i_foo', 'i_bar'] 110 111 112def decay_rhs(t, y, k): 113 ny = len(y) 114 dydt = [0]*ny 115 for idx in range(ny): 116 if idx < ny-1: 117 dydt[idx] -= y[idx]*k[idx] 118 if idx > 0: 119 dydt[idx] += y[idx-1]*k[idx-1] 120 return dydt 121 122 123def _test_TransformedSys(dep_tr, indep_tr, rtol, atol, first_step, forgive=1, y_zero=1e-20, t_zero=1e-12, **kwargs): 124 k = [7., 3, 2] 125 ts = symmetricsys(dep_tr, indep_tr).from_callback( 126 decay_rhs, len(k)+1, len(k)) 127 y0 = [y_zero]*(len(k)+1) 128 y0[0] = 1 129 xout, yout, info = ts.integrate( 130 [t_zero, 1], y0, k, integrator='cvode', atol=atol, rtol=rtol, 131 first_step=first_step, **kwargs) 132 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 133 assert np.allclose(yout, ref, rtol=rtol*forgive, atol=atol*forgive) 134 135 136@requires('sym') 137def test_SymbolicSys__jacobian_singular(): 138 k = (4, 3) 139 odesys = SymbolicSys.from_callback(decay_dydt_factory(k), len(k)+1) 140 assert odesys.jacobian_singular() 141 142 143@requires('sym', 'pycvodes') 144@pycvodes_klu 145def test_SymbolicSys_jacobian_sparse(): 146 k = (4, 3) 147 y0 = (5, 4, 2) 148 odesys = SymbolicSys.from_callback(decay_dydt_factory(k), len(k) + 1, sparse=True) 149 xout, yout, info = odesys.integrate([0, 2], y0, integrator='cvode', linear_solver='klu', 150 atol=1e-12, rtol=1e-12, nsteps=1000) 151 ref = np.array(bateman_full(y0, k+(0,), xout - xout[0], exp=np.exp)).T 152 assert np.allclose(yout, ref, rtol=1e-10, atol=1e-10) 153 154 155@requires('sym', 'pycvodes') 156def test_TransformedSys_liny_linx(): 157 _test_TransformedSys(idty2, idty2, 1e-11, 1e-11, 0, 15) 158 159 160@requires('sym', 'pycvodes') 161def test_TransformedSys_logy_logx(): 162 _test_TransformedSys(get_logexp(), get_logexp(), 1e-7, 1e-7, 1e-4, 150, nsteps=800) 163 164 165@requires('sym', 'pycvodes', 'sympy') 166def test_TransformedSys_logy_logx_scaled_shifted(): 167 import sympy as sp 168 em16 = (sp.S.One*10)**-16 169 _test_TransformedSys(get_logexp(42, em16), get_logexp(42, em16), 1e-7, 1e-7, 1e-4, 170 150, y_zero=0, t_zero=0, nsteps=800) 171 172 173@requires('sym', 'pycvodes') 174def test_TransformedSys_logy_linx(): 175 _test_TransformedSys(get_logexp(), idty2, 1e-8, 1e-8, 0, 150, nsteps=1700) 176 177 178@requires('sym', 'pycvodes') 179def test_TransformedSys_liny_logx(): 180 _test_TransformedSys(idty2, get_logexp(), 1e-9, 1e-9, 0, 150) 181 182 183@requires('sym', 'pycvodes') 184def test_ScaledSys(): 185 import sympy as sp 186 k = k0, k1, k2 = [7., 3, 2] 187 y0, y1, y2, y3 = sp.symbols('y0 y1 y2 y3', real=True, positive=True) 188 # this is actually a silly example since it is linear 189 l = [ 190 (y0, -7*y0), 191 (y1, 7*y0 - 3*y1), 192 (y2, 3*y1 - 2*y2), 193 (y3, 2*y2) 194 ] 195 ss = ScaledSys(l, dep_scaling=1e8) 196 y0 = [0]*(len(k)+1) 197 y0[0] = 1 198 xout, yout, info = ss.integrate([1e-12, 1], y0, integrator='cvode', 199 atol=1e-12, rtol=1e-12, nsteps=1000) 200 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 201 assert np.allclose(yout, ref, rtol=2e-11, atol=2e-11) 202 203 204@requires('sym', 'scipy', 'pycvodes') 205@pytest.mark.parametrize('nbody', [2, 3, 4, 5]) 206def test_ScaledSysByName(nbody): 207 sfact = nbody * 7 208 kwargs = dict(names=['foo', 'bar'], dep_scaling=sfact) 209 210 def nmerization(x, y, p): 211 return [-nbody*p[0]*y[0]**nbody, nbody*p[0]*y[0]**nbody] 212 213 odesys = ScaledSys.from_callback(nmerization, 2, 1, **kwargs) 214 assert odesys.autonomous_interface is True 215 with pytest.raises(TypeError): 216 odesys.integrate(1, [0]) 217 218 def nmerization_name(x, y, p): 219 return {'foo': -nbody*p[0]*y['foo']**nbody, 'bar': nbody*p[0]*y['foo']**nbody} 220 221 odesys2 = ScaledSys.from_callback(nmerization_name, dep_by_name=True, nparams=1, **kwargs) 222 assert odesys2.autonomous_interface is True 223 k = 5 224 foo0 = 2 225 for system, y0 in zip([odesys, odesys2], [[foo0, 3], {'foo': foo0, 'bar': 3}]): 226 xout, yout, info = system.integrate(1, y0, [k], integrator='cvode', nsteps=707*1.01, 227 first_step=1e-3, atol=1e-10, rtol=1e-10) 228 _r = (1/(foo0**(1-nbody) + nbody*k*xout*(nbody-1)))**(1/(nbody-1)) 229 assert np.allclose(yout[:, 0], _r, atol=1e-9, rtol=1e-9) 230 assert np.allclose(yout[:, 1], 3 + 2 - _r, atol=1e-9, rtol=1e-9) 231 232 233@requires('sym', 'scipy') 234def test_ScaledSys_from_callback(): 235 # this is actually a silly example since it is linear 236 def f(t, x, k): 237 return [-k[0]*x[0], 238 k[0]*x[0] - k[1]*x[1], 239 k[1]*x[1] - k[2]*x[2], 240 k[2]*x[2]] 241 odesys = ScaledSys.from_callback(f, 4, 3, 3.14e8) 242 k = [7, 3, 2] 243 y0 = [0]*(len(k)+1) 244 y0[0] = 1 245 xout, yout, info = odesys.integrate([1e-12, 1], y0, k, integrator='scipy') 246 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 247 assert np.allclose(yout, ref, rtol=3e-11, atol=3e-11) 248 249 with pytest.raises(TypeError): 250 odesys.integrate([1e-12, 1], [0]*len(k), k, integrator='scipy') 251 252 253@requires('sym', 'scipy') 254def test_ScaledSys_from_callback__exprs(): 255 def f(t, x, k): 256 return [-k[0]*x[0]*x[0]*t] 257 x, y, nfo = SymbolicSys.from_callback(f, 1, 1).integrate( 258 [0, 1], [3.14], [2.78]) 259 xs, ys, nfos = ScaledSys.from_callback(f, 1, 1, 100).integrate( 260 [0, 1], [3.14], [2.78]) 261 from scipy.interpolate import interp1d 262 cb = interp1d(x, y[:, 0]) 263 cbs = interp1d(xs, ys[:, 0]) 264 t = np.linspace(0, 1) 265 assert np.allclose(cb(t), cbs(t)) 266 267 268def timeit(callback, *args, **kwargs): 269 t0 = perf_counter() 270 result = callback(*args, **kwargs) 271 return perf_counter() - t0, result 272 273 274@requires('sym', 'pyodeint') 275@pytest.mark.parametrize('method', ['bs', 'rosenbrock4']) 276def test_exp(method): 277 import sympy as sp 278 x = sp.Symbol('x') 279 symsys = SymbolicSys([(x, sp.exp(x))]) 280 tout = [0, 1e-9, 1e-7, 1e-5, 1e-3, 0.1] 281 xout, yout, info = symsys.integrate( 282 tout, [1], method=method, integrator='odeint', atol=1e-12, rtol=1e-12) 283 e = math.e 284 ref = -math.log(1/e - 0.1) 285 assert abs(yout[-1, 0] - ref) < 4e-8 286 287 288# @pytest.mark.xfail 289def _test_mpmath(): # too slow 290 import sympy as sp 291 x = sp.Symbol('x') 292 symsys = SymbolicSys([(x, sp.exp(x))]) 293 tout = [0, 1e-9, 1e-7, 1e-5, 1e-3, 0.1] 294 xout, yout, info = symsys.integrate(tout, [1], integrator='mpmath') 295 e = math.e 296 ref = -math.log(1/e - 0.1) 297 assert abs(yout[-1, 0] - ref) < 4e-8 298 299 300def decay_dydt_factory(k, names=None): 301 # Generates a callback for evaluating a dydt-callback for 302 # a chain of len(k) + 1 species with len(k) decays 303 # with corresponding decay constants k 304 ny = len(k) + 1 305 306 def dydt(t, y): 307 exprs = [] 308 for idx in range(ny): 309 expr = 0 310 curr_key = idx 311 prev_key = idx - 1 312 if names is not None: 313 curr_key = names[curr_key] 314 prev_key = names[prev_key] 315 if idx < ny-1: 316 expr -= y[curr_key]*k[curr_key] 317 if idx > 0: 318 expr += y[prev_key]*k[prev_key] 319 exprs.append(expr) 320 if names is None: 321 return exprs 322 else: 323 return dict(zip(names, exprs)) 324 325 return dydt 326 327 328# Short decay chains, using Bateman's equation 329# -------------------------------------------- 330@requires('sym', 'scipy') 331@pytest.mark.parametrize('band', [(1, 0), None]) 332def test_SymbolicSys__from_callback_bateman(band): 333 # Decay chain of 3 species (2 decays) 334 # A --[k0=4]--> B --[k1=3]--> C 335 tend, k, y0 = 2, [4, 3], (5, 4, 2) 336 atol, rtol = 1e-11, 1e-11 337 odesys = SymbolicSys.from_callback(decay_dydt_factory(k), len(k)+1, 338 band=band) 339 xout, yout, info = odesys.integrate( 340 tend, y0, atol=atol, integrator='scipy', rtol=rtol) 341 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 342 assert np.allclose(yout, ref, rtol=rtol, atol=atol) 343 344 345def _test_bateman(SymbSys, **kwargs): 346 import sympy as sp 347 tend, k, y0 = 2, [4, 3], (5, 4, 2) 348 y = sp.symarray('y', len(k)+1) 349 dydt = decay_dydt_factory(k) 350 f = dydt(0, y) 351 odesys = SymbSys(zip(y, f), **kwargs) 352 xout, yout, info = odesys.integrate(tend, y0, integrator='scipy') 353 ref = np.array(bateman_full(y0, k+[0], xout-xout[0], exp=np.exp)).T 354 assert np.allclose(yout, ref) 355 356 357@requires('sym', 'scipy') 358@pytest.mark.parametrize('band', [(1, 0), None]) 359def test_SymbolicSys_bateman(band): 360 _test_bateman(SymbolicSys, band=band) 361 362 363@requires('sym', 'scipy') 364@pytest.mark.parametrize('band', [(1, 0), None]) 365def test_ScaledSys_bateman(band): 366 _test_bateman(ScaledSys, band=band, dep_scaling=1e3) 367 368 369# Longer chains with careful choice of parameters 370# ----------------------------------------------- 371 372 373def analytic1(i, p, a): 374 assert i > 0 375 assert p >= 0 376 assert a > 0 377 from scipy.special import binom 378 return binom(p+i-1, p) * a**(i-1) * (a+1)**(-i-p) 379 380 381def check(vals, n, p, a, atol, rtol, forgiveness=1): 382 # Check solution vs analytic reference: 383 for i in range(n-1): 384 val = vals[i] 385 ref = analytic1(i+1, p, a) 386 diff = val - ref 387 acceptance = (atol + abs(ref)*rtol)*forgiveness 388 assert abs(diff) < acceptance 389 390 391def get_special_chain(n, p, a, **kwargs): 392 assert n > 1 393 assert p >= 0 394 assert a > 0 395 y0 = np.zeros(n) 396 y0[0] = 1 397 k = [(i+p+1)*math.log(a+1) for i in range(n-1)] 398 dydt = decay_dydt_factory(k) 399 return y0, k, SymbolicSys.from_callback(dydt, n, **kwargs) 400 401 402@requires('sym', 'scipy') 403@pytest.mark.parametrize('p', [0, 1, 2, 3]) 404def test_check(p): 405 n, a = 7, 5 406 y0, k, _odesys = get_special_chain(n, p, a) 407 vals = bateman_full(y0, k+[0], 1, exp=np.exp) 408 check(vals, n, p, a, atol=1e-12, rtol=1e-12) 409 410 411# adaptive stepsize with vode is performing ridiculously 412# poorly for this problem 413@requires('sym', 'scipy') 414@pytest.mark.parametrize('name,forgive', zip( 415 'dopri5 dop853 vode'.split(), (1, 1, (3, 3e6)))) 416def test_scipy(name, forgive): 417 n, p, a = 13, 1, 13 418 atol, rtol = 1e-10, 1e-10 419 y0, k, odesys_dens = get_special_chain(n, p, a) 420 if name == 'vode': 421 tout = [0]+[10**i for i in range(-10, 1)] 422 xout, yout, info = odesys_dens.integrate( 423 tout, y0, integrator='scipy', name=name, atol=atol, rtol=rtol) 424 check(yout[-1, :], n, p, a, atol, rtol, forgive[0]) 425 426 xout, yout, info = odesys_dens.integrate( 427 1, y0, integrator='scipy', name=name, atol=atol, rtol=rtol) 428 check(yout[-1, :], n, p, a, atol, rtol, forgive[1]) 429 430 else: 431 xout, yout, info = odesys_dens.integrate( 432 1, y0, integrator='scipy', name=name, atol=atol, rtol=rtol) 433 check(yout[-1, :], n, p, a, atol, rtol, forgive) 434 assert yout.shape[0] > 2 435 436 437# (dopri5, .2), (bs, .03) <-- works in boost 1.59 438@requires('sym', 'pyodeint') 439@pytest.mark.parametrize('method,forgive', zip( 440 'rosenbrock4 dopri5 bs'.split(), (.2, .2, .04))) 441def test_odeint(method, forgive): 442 n, p, a = 13, 1, 13 443 atol, rtol = 1e-10, 1e-10 444 y0, k, odesys_dens = get_special_chain(n, p, a) 445 dydt = decay_dydt_factory(k) 446 odesys_dens = SymbolicSys.from_callback(dydt, len(k)+1) 447 # adaptive stepper fails to produce the accuracy asked for. 448 xout, yout, info = odesys_dens.integrate( 449 [10**i for i in range(-15, 1)], y0, integrator='odeint', 450 method=method, atol=atol, rtol=rtol) 451 check(yout[-1, :], n, p, a, atol, rtol, forgive) 452 453 454def _gsl(tout, method, forgive): 455 n, p, a = 13, 1, 13 456 atol, rtol = 1e-10, 1e-10 457 y0, k, odesys_dens = get_special_chain(n, p, a) 458 dydt = decay_dydt_factory(k) 459 odesys_dens = SymbolicSys.from_callback(dydt, len(k)+1) 460 # adaptive stepper fails to produce the accuracy asked for. 461 xout, yout, info = odesys_dens.integrate( 462 tout, y0, method=method, integrator='gsl', atol=atol, rtol=rtol) 463 check(yout[-1, :], n, p, a, atol, rtol, forgive) 464 465 466@requires('sym', 'pygslodeiv2') 467@pytest.mark.parametrize('method,forgive', zip( 468 'msadams msbdf rkck bsimp'.split(), (5, 14, 0.2, 0.02))) 469def test_gsl_predefined(method, forgive): 470 _gsl([10**i for i in range(-14, 1)], method, forgive) 471 472 473@requires('sym', 'pygslodeiv2') 474@pytest.mark.parametrize('method,forgive', zip( 475 'bsimp msadams msbdf rkck'.split(), (0.01, 4, 14, 0.21))) 476def test_gsl_adaptive(method, forgive): 477 _gsl(1, method, forgive) 478 479 480def _cvode(tout, method, forgive): 481 n, p, a = 13, 1, 13 482 atol, rtol = 1e-10, 1e-10 483 y0, k, odesys_dens = get_special_chain(n, p, a) 484 dydt = decay_dydt_factory(k) 485 odesys_dens = SymbolicSys.from_callback(dydt, len(k)+1) 486 # adaptive stepper fails to produce the accuracy asked for. 487 xout, yout, info = odesys_dens.integrate( 488 tout, y0, method=method, integrator='cvode', atol=atol, rtol=rtol) 489 check(yout[-1, :], n, p, a, atol, rtol, forgive) 490 491 492@requires('sym', 'pycvodes') 493@pytest.mark.parametrize('method,forgive', zip( 494 'adams bdf'.split(), (2.4, 5.0))) 495def test_cvode_predefined(method, forgive): 496 _cvode([10**i for i in range(-15, 1)], method, forgive) 497 498 499# cvode performs significantly better than vode: 500@requires('sym', 'pycvodes') 501@pytest.mark.parametrize('method,forgive', zip( 502 'adams bdf'.split(), (2.4, 5))) 503def test_cvode_adaptive(method, forgive): 504 _cvode(1, method, forgive) 505 506 507@requires('sym', 'scipy') 508@pytest.mark.parametrize('n,forgive', [(4, 1), (17, 1), (42, 7)]) 509def test_long_chain_dense(n, forgive): 510 p, a = 0, n 511 y0, k, odesys_dens = get_special_chain(n, p, a) 512 atol, rtol = 1e-12, 1e-12 513 tout = 1 514 xout, yout, info = odesys_dens.integrate( 515 tout, y0, integrator='scipy', atol=atol, rtol=rtol) 516 check(yout[-1, :], n, p, a, atol, rtol, forgive) 517 518 519@requires('sym', 'scipy') 520@pytest.mark.parametrize('n', [29]) # something maxes out at 31 521def test_long_chain_banded_scipy(n): 522 p, a = 0, n 523 y0, k, odesys_dens = get_special_chain(n, p, a) 524 y0, k, odesys_band = get_special_chain(n, p, a, band=(1, 0)) 525 atol, rtol = 1e-7, 1e-7 526 tout = np.logspace(-10, 0, 10) 527 528 def mk_callback(odesys): 529 def callback(*args, **kwargs): 530 return odesys.integrate(*args, integrator='scipy', **kwargs) 531 return callback 532 min_time_dens, min_time_band = float('inf'), float('inf') 533 for _ in range(3): # warmup 534 time_dens, (xout_dens, yout_dens, info) = timeit( 535 mk_callback(odesys_dens), tout, y0, atol=atol, rtol=rtol, 536 name='vode', method='bdf', first_step=1e-10) 537 assert info['njev'] > 0 538 min_time_dens = min(min_time_dens, time_dens) 539 for _ in range(3): # warmup 540 time_band, (xout_band, yout_band, info) = timeit( 541 mk_callback(odesys_band), tout, y0, atol=atol, rtol=rtol, 542 name='vode', method='bdf', first_step=1e-10) 543 assert info['njev'] > 0 544 min_time_band = min(min_time_band, time_band) 545 check(yout_dens[-1, :], n, p, a, atol, rtol, 1.5) 546 check(yout_band[-1, :], n, p, a, atol, rtol, 1.5) 547 assert min_time_dens*2 > min_time_band # (2x: fails sometimes due to load) 548 549 550@requires('sym', 'pycvodes') 551@pytest.mark.parametrize('n', [29, 79]) 552def test_long_chain_banded_cvode(n): 553 p, a = 0, n 554 y0, k, odesys_dens = get_special_chain(n, p, a) 555 y0, k, odesys_band = get_special_chain(n, p, a, band=(1, 0)) 556 atol, rtol = 1e-9, 1e-9 557 558 def mk_callback(odesys): 559 def callback(*args, **kwargs): 560 return odesys.integrate(*args, integrator='cvode', **kwargs) 561 return callback 562 for _ in range(2): # warmup 563 time_band, (xout_band, yout_band, info) = timeit( 564 mk_callback(odesys_band), 1, y0, atol=atol, rtol=rtol) 565 assert info['njev'] > 0 566 for _ in range(2): # warmup 567 time_dens, (xout_dens, yout_dens, info) = timeit( 568 mk_callback(odesys_dens), 1, y0, atol=atol, rtol=rtol) 569 assert info['njev'] > 0 570 check(yout_dens[-1, :], n, p, a, atol, rtol, 7) 571 check(yout_band[-1, :], n, p, a, atol, rtol, 25) # suspicious 572 assert info['njev'] > 0 573 try: 574 assert time_dens > time_band 575 except AssertionError: 576 pass # will fail sometimes due to load 577 578 579@requires('sym', 'pycvodes', 'pygslodeiv2') 580@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 581def test_no_diff_adaptive_auto_switch_single(integrator): 582 odesys = _get_decay3() 583 tout, y0, k = [3, 5], [3, 2, 1], [3.5, 2.5, 1.5] 584 xout1, yout1, info1 = odesys.integrate(tout, y0, k, integrator=integrator) 585 ref = np.array(bateman_full(y0, k, xout1 - xout1[0], exp=np.exp)).T 586 assert info1['success'] 587 assert xout1.size > 10 588 assert xout1.size == yout1.shape[0] 589 assert np.allclose(yout1, ref) 590 591 xout2, yout2, info2 = integrate_auto_switch([odesys], {}, tout, y0, k, integrator=integrator) 592 assert info1['success'] 593 assert xout2.size == xout1.size 594 assert np.allclose(yout2, ref) 595 596 597@requires('sym', 'pycvodes', 'pygslodeiv2') 598@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 599def test_no_diff_adaptive_auto_switch_single__multimode(integrator): 600 odesys = _get_decay3() 601 tout = [[3, 5], [4, 6], [6, 8], [9, 11]] 602 _y0 = [3, 2, 1] 603 y0 = [_y0]*4 604 _k = [3.5, 2.5, 1.5] 605 k = [_k]*4 606 res1 = odesys.integrate(tout, y0, k, integrator=integrator, first_step=1e-14) 607 for res in res1: 608 xout1, yout1, info1 = res.xout, res.yout, res.info 609 ref = np.array(bateman_full(_y0, _k, xout1 - xout1[0], exp=np.exp)).T 610 assert info1['success'] 611 assert xout1.size > 10 612 assert xout1.size == yout1.shape[0] 613 assert np.allclose(yout1, ref) 614 615 res2 = integrate_auto_switch([odesys], {}, tout, y0, k, integrator=integrator, first_step=1e-14) 616 for res in res2: 617 xout2, yout2, info2 = res.xout, res.yout, res.info 618 assert info2['success'] 619 assert xout2.size == xout1.size 620 assert np.allclose(yout2, ref) 621 622 623@requires('sym', 'pycvodes', 'pygslodeiv2') 624@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 625def test_PartiallySolvedSystem(integrator): 626 odesys = _get_decay3(lower_bounds=[0, 0, 0]) 627 partsys = PartiallySolvedSystem(odesys, lambda x0, y0, p0, be: { 628 odesys.dep[0]: y0[0]*be.exp(-p0[0]*(odesys.indep-x0)) 629 }) 630 y0 = [3, 2, 1] 631 k = [3.5, 2.5, 1.5] 632 xout, yout, info = partsys.integrate( 633 [0, 1], y0, k, integrator=integrator) 634 ref = np.array(bateman_full(y0, k, xout - xout[0], exp=np.exp)).T 635 assert info['success'] 636 assert np.allclose(yout, ref) 637 638 639@requires('sym', 'pycvodes') 640def test_PartiallySolvedSystem_ScaledSys(): 641 odesys = _get_decay3(lower_bounds=[0, 0, 0]) 642 partsys = PartiallySolvedSystem(odesys, lambda x0, y0, p0, be: { 643 odesys.dep[0]: y0[0]*be.exp(-p0[0]*(odesys.indep-x0)) 644 }) 645 y0 = [3, 2, 1] 646 k = [3.5, 2.5, 1.5] 647 648 def _check(res): 649 ref = np.array(bateman_full(y0, k, res.xout - res.xout[0], exp=np.exp)).T 650 assert res.info['success'] 651 assert np.allclose(res.yout, ref) 652 args = [0, 1], y0, k 653 kwargs = dict(integrator='cvode') 654 _check(odesys.integrate(*args, **kwargs)) 655 _check(partsys.integrate(*args, **kwargs)) 656 scaledsys = ScaledSys.from_other(partsys, dep_scaling=42, indep_scaling=17) 657 _check(scaledsys.integrate(*args, **kwargs)) 658 659 660@requires('sym', 'pycvodes', 'pygslodeiv2') 661@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 662def test_PartiallySolvedSystem_multi(integrator): 663 odesys = _get_decay3() 664 665 def _get_analytic(x0, y0, p0, be): 666 a0 = y0[0]*be.exp(-p0[0]*(odesys.indep - x0)) 667 a1 = y0[0] + y0[1] + y0[2] - a0 - odesys.dep[2] 668 return [a0, a1] 669 670 def subst(x0, y0, p0, be): 671 a0, a1 = _get_analytic(x0, y0, p0, be) 672 return {odesys.dep[0]: a0, odesys.dep[1]: a1} 673 674 partsys = PartiallySolvedSystem(odesys, subst) 675 a0, a1 = _get_analytic(partsys.init_indep, 676 partsys.init_dep, 677 odesys.params, 678 odesys.be) 679 assert partsys.ny == 1 680 assert partsys.exprs[0].subs(odesys.params[2], 0) - odesys.params[1]*a1 == 0 681 682 683@requires('sym', 'pycvodes', 'pygslodeiv2') 684@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 685def test_PartiallySolvedSystem__using_y(integrator): 686 odesys = _get_decay3() 687 partsys = PartiallySolvedSystem(odesys, lambda x0, y0, p0: { 688 odesys.dep[2]: y0[0] + y0[1] + y0[2] - odesys.dep[0] - odesys.dep[1] 689 }) 690 y0 = [3, 2, 1] 691 k = [3.5, 2.5, 0] 692 xout, yout, info = partsys.integrate([0, 1], y0, k, integrator=integrator) 693 ref = np.array(bateman_full(y0, k, xout - xout[0], exp=np.exp)).T 694 assert info['success'] 695 assert np.allclose(yout, ref) 696 assert np.allclose(np.sum(yout, axis=1), sum(y0)) 697 698 699@requires('sym', 'pycvodes', 'pygslodeiv2') 700@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 701def test_PartiallySolvedSystem_multiple_subs(integrator): 702 odesys = _get_decay3(lower_bounds=[0, 0, 0]) 703 704 def substitutions(x0, y0, p0, be): 705 analytic0 = y0[0]*be.exp(-p0[0]*(odesys.indep-x0)) 706 analytic2 = y0[0] + y0[1] + y0[2] - analytic0 - odesys.dep[1] 707 return {odesys.dep[0]: analytic0, odesys.dep[2]: analytic2} 708 709 partsys = PartiallySolvedSystem(odesys, substitutions) 710 y0 = [3, 2, 1] 711 k = [3.5, 2.5, 0] 712 xout, yout, info = partsys.integrate([0, 1], y0, k, integrator=integrator) 713 ref = np.array(bateman_full(y0, k, xout - xout[0], exp=np.exp)).T 714 assert info['success'] 715 assert np.allclose(yout, ref) 716 717 718@requires('sym', 'pycvodes', 'pygslodeiv2') 719@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 720def test_PartiallySolvedSystem_multiple_subs__transformed(integrator): 721 odesys = _get_decay3(lower_bounds=[0, 0, 0]) 722 723 def substitutions(x0, y0, p0, be): 724 analytic0 = y0[0]*be.exp(-p0[0]*(odesys.indep-x0)) 725 analytic2 = y0[0] + y0[1] + y0[2] - analytic0 - odesys.dep[1] 726 return {odesys.dep[0]: analytic0, odesys.dep[2]: analytic2} 727 728 partsys = PartiallySolvedSystem(odesys, substitutions) 729 LogLogSys = symmetricsys(get_logexp(), get_logexp()) 730 loglogpartsys = LogLogSys.from_other(partsys) 731 y0 = [3, 2, 1] 732 k = [3.5, 2.5, 0] 733 tend = 1 734 for system, ny_internal in [(odesys, 3), (partsys, 1), (loglogpartsys, 1)]: 735 xout, yout, info = system.integrate([1e-12, tend], y0, k, integrator=integrator, 736 first_step=1e-14, nsteps=1000) 737 ref = np.array(bateman_full(y0, k, xout - xout[0], exp=np.exp)).T 738 assert info['success'] 739 assert info['internal_yout'].shape[-1] == ny_internal 740 if system == loglogpartsys: 741 assert info['internal_yout'][-1, 0] < 0 # ln(y[1]) 742 assert np.allclose(yout, ref) 743 744 745def _get_transf_part_system(): 746 import sympy as sp 747 odesys = _get_decay3() 748 partsys = PartiallySolvedSystem(odesys, lambda x0, y0, p0: { 749 odesys.dep[0]: y0[0]*sp.exp(-p0[0]*(odesys.indep-x0)) 750 }) 751 LogLogSys = symmetricsys(get_logexp(), get_logexp()) 752 return LogLogSys.from_other(partsys) 753 754 755@requires('sym', 'pycvodes', 'pygslodeiv2') 756@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 757def test_PartiallySolvedSystem__symmetricsys(integrator): 758 trnsfsys = _get_transf_part_system() 759 y0 = [3., 2., 1.] 760 k = [3.5, 2.5, 0] 761 xout, yout, info = trnsfsys.integrate([1e-10, 1], y0, k, integrator=integrator, first_step=1e-14) 762 ref = np.array(bateman_full(y0, k, xout - xout[0], exp=np.exp)).T 763 assert info['success'] 764 assert np.allclose(yout, ref) 765 assert np.allclose(np.sum(yout, axis=1), sum(y0)) 766 767 768@requires('sym', 'pycvodes', 'pygslodeiv2') 769@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 770def test_PartiallySolvedSystem__symmetricsys__multi(integrator): 771 trnsfsys = _get_transf_part_system() 772 y0s = [[3., 2., 1.], [3.1, 2.1, 1.1], [3.2, 2.3, 1.2], [3.6, 2.4, 1.3]] 773 ks = [[3.5, 2.5, 0], [3.3, 2.4, 0], [3.2, 2.1, 0], [3.3, 2.4, 0]] 774 results = trnsfsys.integrate([(1e-10, 1)]*len(ks), y0s, ks, integrator=integrator, first_step=1e-14) 775 for i, (y0, k) in enumerate(zip(y0s, ks)): 776 xout, yout, info = results[i] 777 ref = np.array(bateman_full(y0, k, xout - xout[0], exp=np.exp)).T 778 assert info['success'] and info['nfev'] > 10 779 assert info['nfev'] > 1 and info['time_cpu'] < 100 780 assert np.allclose(yout, ref) and np.allclose(np.sum(yout, axis=1), sum(y0)) 781 782 783def _get_nonlin(**kwargs): 784 return SymbolicSys.from_callback( 785 lambda x, y, p: [ 786 -p[0]*y[0]*y[1] + p[1]*y[2], 787 -p[0]*y[0]*y[1] + p[1]*y[2], 788 p[0]*y[0]*y[1] - p[1]*y[2] 789 ], 3, 2, **kwargs) 790 791 792def _get_nonlin_part_system(): 793 odesys = _get_nonlin() 794 return PartiallySolvedSystem(odesys, lambda x0, y0, p0: { 795 odesys.dep[0]: y0[0] + y0[2] - odesys.dep[2] 796 }) 797 798 799def _ref_nonlin(y0, k, t): 800 X, Y, Z = y0[2], max(y0[:2]), min(y0[:2]) 801 kf, kb = k 802 x0 = Y*kf 803 x1 = Z*kf 804 x2 = 2*X*kf 805 x3 = -kb - x0 - x1 806 x4 = -x2 + x3 807 x5 = np.sqrt(-4*kf*(X**2*kf + X*x0 + X*x1 + Z*x0) + x4**2) 808 x6 = kb + x0 + x1 + x5 809 x7 = (x3 + x5)*np.exp(-t*x5) 810 x8 = x3 - x5 811 return (x4*x8 + x5*x8 + x7*(x2 + x6))/(2*kf*(x6 + x7)) 812 813 814@requires('sym', 'pycvodes', 'pygslodeiv2') 815@pytest.mark.parametrize('integrator', ['cvode', 'gsl']) 816def test_PartiallySolvedSystem__symmetricsys__nonlinear(integrator): 817 partsys = _get_nonlin_part_system() 818 logexp = get_logexp(7, partsys.indep**0/10**7) 819 trnsfsys = symmetricsys(logexp, logexp).from_other(partsys) 820 y0 = [3., 2., 1.] 821 k = [9.351, 2.532] 822 tend = 1.7 823 atol, rtol = 1e-12, 1e-13 824 for odesys, forgive in [(partsys, 21), (trnsfsys, 298)]: 825 xout, yout, info = odesys.integrate(tend, y0, k, integrator=integrator, 826 first_step=1e-14, atol=atol, rtol=rtol, 827 nsteps=1000) 828 assert info['success'] 829 yref = np.empty_like(yout) 830 yref[:, 2] = _ref_nonlin(y0, k, xout - xout[0]) 831 yref[:, 0] = y0[0] + y0[2] - yref[:, 2] 832 yref[:, 1] = y0[1] + y0[2] - yref[:, 2] 833 assert np.allclose(yout, yref, atol=forgive*atol, rtol=forgive*rtol) 834 835 836@requires('sym', 'scipy') 837def test_SymbolicSys_from_other(): 838 scaled = ScaledSys.from_callback(lambda x, y: [y[0]*y[0]], 1, 839 dep_scaling=101) 840 LogLogSys = symmetricsys(get_logexp(), get_logexp()) 841 transformed_scaled = LogLogSys.from_other(scaled) 842 tout = np.array([0, .2, .5]) 843 y0 = [1.] 844 ref, nfo1 = ODESys(lambda x, y: y[0]*y[0]).predefined( 845 y0, tout, first_step=1e-14) 846 analytic = 1/(1-tout.reshape(ref.shape)) 847 assert np.allclose(ref, analytic) 848 yout, nfo0 = transformed_scaled.predefined(y0, tout+1) 849 assert np.allclose(yout, analytic) 850 851 852@requires('sym', 'scipy') 853def test_backend(): 854 855 def f(x, y, p, backend=math): 856 return [backend.exp(p[0]*y[0])] 857 858 def analytic(x, p, y0): 859 # dydt = exp(p*y(t)) 860 # y(t) = - log(p*(c1-t))/p 861 # 862 # y(0) = - log(p*c1)/p 863 # p*y(0) = -log(p) -log(c1) 864 # c1 = exp(-log(p)-p*y(0)) 865 # c1 = 866 # 867 # y(t) = -log(p*(exp(-p*y(0))/p - t))/p 868 return -np.log(p*(np.exp(-p*y0)/p - x))/p 869 870 y0, tout, p = .07, [0, .1, .2], .3 871 ref = analytic(tout, p, y0) 872 873 def _test_odesys(odesys): 874 yout, info = odesys.predefined([y0], tout, [p]) 875 assert np.allclose(yout.flatten(), ref) 876 877 _test_odesys(ODESys(f)) 878 _test_odesys(SymbolicSys.from_callback(f, 1, 1)) 879 880 881@requires('sym', 'scipy') 882def _test_SymbolicSys_from_callback__backend(backend): 883 ss = SymbolicSys.from_callback(vdp_f, 2, 1, backend=backend) 884 xout, yout, info = ss.integrate([0, 1, 2], [1, 0], params=[2.0]) 885 # blessed values: 886 ref = [[1, 0], [0.44449086, -1.32847148], [-1.89021896, -0.71633577]] 887 assert np.allclose(yout, ref) 888 assert info['nfev'] > 0 889 890 891@requires('sym', 'sympy') 892def test_SymbolicSys_from_callback__sympy(): 893 _test_SymbolicSys_from_callback__backend('sympy') 894 895 896@requires('sym', 'symengine') 897def test_SymbolicSys_from_callback__symengine(): 898 _test_SymbolicSys_from_callback__backend('symengine') 899 900 901@requires('sym', 'symcxx') 902def test_SymbolicSys_from_callback__symcxx(): 903 _test_SymbolicSys_from_callback__backend('symcxx') 904 905 906@requires('sym', 'pycvodes', 'pygslodeiv2') 907@pytest.mark.parametrize('integrator,method', [('cvode', 'adams'), ('gsl', 'msadams')]) 908def test_integrate_auto_switch(integrator, method): 909 for p in (0, 1, 2, 3): 910 n, a = 7, 5 911 atol, rtol = 1e-10, 1e-10 912 y0, k, linsys = get_special_chain(n, p, a) 913 y0 += 1e-10 914 LogLogSys = symmetricsys(get_logexp(), get_logexp()) 915 logsys = LogLogSys.from_other(linsys) 916 tout = [10**-12, 1] 917 kw = dict( 918 integrator=integrator, method=method, atol=atol, rtol=rtol, 919 ) 920 forgive = (5+p)*1.2 921 922 xout, yout, info = integrate_auto_switch([logsys, linsys], {'nsteps': [1, 1]}, tout, y0, 923 return_on_error=True, **kw) 924 assert info['success'] == False # noqa 925 ntot = 400 926 nlinear = 60*(p+3) 927 928 xout, yout, info = integrate_auto_switch([logsys, linsys], { 929 'nsteps': [ntot - nlinear, nlinear], 930 'first_step': [30.0, 1e-5], 931 'return_on_error': [True, False] 932 }, tout, y0, **kw) 933 assert info['success'] == True # noqa 934 check(yout[-1, :], n, p, a, atol, rtol, forgive) 935 936 937def _test_cetsa(y0, params, extra=False, stepx=1, **kwargs): 938 # real-world based test-case 939 from ._cetsa import _get_cetsa_odesys 940 molar_unitless = 1e9 941 t0, tend = 1e-16, 180 942 odesys = _get_cetsa_odesys(molar_unitless, False) 943 tsys = _get_cetsa_odesys(molar_unitless, True) 944 if y0.ndim == 1: 945 tout = [t0, tend] 946 elif y0.ndim == 2: 947 tout = np.asarray([(t0, tend)]*y0.shape[0]) 948 949 comb_res = integrate_auto_switch([tsys, odesys], {'nsteps': [500*stepx, 20*stepx]}, tout, y0/molar_unitless, 950 params, return_on_error=True, autorestart=2, **kwargs) 951 if isinstance(comb_res, list): 952 for r in comb_res: 953 assert r.info['success'] 954 assert r.info['nfev'] > 10 955 else: 956 assert comb_res.info['success'] 957 assert comb_res.info['nfev'] > 10 958 959 if extra: 960 with pytest.raises(RuntimeError): # (failure) 961 odesys.integrate(np.linspace(t0, tend, 20), y0/molar_unitless, params, atol=1e-7, rtol=1e-7, 962 nsteps=500, first_step=1e-14, **kwargs) 963 964 res = odesys.integrate(np.linspace(t0, tend, 20), y0/molar_unitless, params, nsteps=int(38*1.1), 965 first_step=1e-14, **kwargs) 966 assert np.min(res.yout[-1, :]) < -1e-6 # crazy! (failure of the linear formulation) 967 tres = tsys.integrate([t0, tend], y0/molar_unitless, params, nsteps=int(1345*1.1), **kwargs) 968 assert tres.info['success'] is True 969 assert tres.info['nfev'] > 100 970 971 972@requires('sym', 'pycvodes', 'pygslodeiv2', 'pyodeint') 973@pytest.mark.parametrize('integrator', ['cvode', 'gsl', 'odeint']) 974def test_cetsa(integrator): 975 _test_cetsa(_cetsa.ys[1], _cetsa.ps[1], integrator=integrator, first_step=1e-14, 976 stepx=2 if integrator == 'odeint' else 1) 977 if integrator == 'cvode': 978 _test_cetsa(_cetsa.ys[0], _cetsa.ps[0], extra=True, integrator=integrator) 979 980 981@requires('sym', 'pycvodes', 'pygslodeiv2', 'pyodeint') 982@pytest.mark.parametrize('integrator', ['cvode', 'gsl', 'odeint']) 983def test_cetsa_multi(integrator): 984 _test_cetsa(np.asarray(_cetsa.ys), np.asarray(_cetsa.ps), integrator=integrator, first_step=1e-14, 985 stepx=2 if integrator == 'odeint' else 1) 986 987 988@requires('sym', 'pycvodes') 989def test_dep_by_name(): 990 def _sin(t, y, p): 991 return {'prim': y['bis'], 'bis': -p[0]**2 * y['prim']} 992 odesys = SymbolicSys.from_callback(_sin, names=['prim', 'bis'], nparams=1, dep_by_name=True) 993 A, k = 2, 3 994 for y0 in ({'prim': 0, 'bis': A*k}, [0, A*k]): 995 xout, yout, info = odesys.integrate(np.linspace(0, 1), y0, [k], 996 integrator='cvode', method='adams') 997 assert info['success'] 998 assert xout.size > 7 999 ref = [ 1000 A*np.sin(k*(xout - xout[0])), 1001 A*np.cos(k*(xout - xout[0]))*k 1002 ] 1003 assert np.allclose(yout[:, 0], ref[0], atol=1e-5, rtol=1e-5) 1004 assert np.allclose(yout[:, 1], ref[1], atol=1e-5, rtol=1e-5) 1005 1006 1007def _get_cetsa_isothermal(): 1008 # tests rhs returning dict, integrate with y0 being a dict & multiple solved variables 1009 names = ('NL', 'N', 'L', 'U', 'A') 1010 k_names = ('dis', 'as', 'un', 'fo', 'ag') 1011 1012 def i(n): 1013 return names.index(n) 1014 1015 def k(n): 1016 return k_names.index(n) 1017 1018 def rhs(x, y, p): 1019 r = { 1020 'diss': p['dis']*y['NL'], 1021 'asso': p['as']*y['N']*y['L'], 1022 'unfo': p['un']*y['N'], 1023 'fold': p['fo']*y['U'], 1024 'aggr': p['ag']*y['U'] 1025 } 1026 return { 1027 'NL': r['asso'] - r['diss'], 1028 'N': r['diss'] - r['asso'] + r['fold'] - r['unfo'], 1029 'L': r['diss'] - r['asso'], 1030 'U': r['unfo'] - r['fold'] - r['aggr'], 1031 'A': r['aggr'] 1032 } 1033 1034 return SymbolicSys.from_callback(rhs, dep_by_name=True, par_by_name=True, names=names, param_names=k_names) 1035 1036 1037@requires('sym', 'scipy') 1038def test_cetsa_isothermal(): 1039 odesys = _get_cetsa_isothermal() 1040 tout = (0, 300) 1041 par = {'dis': 10.0, 'as': 1e9, 'un': 0.1, 'fo': 2.0, 'ag': 0.05} 1042 conc0 = defaultdict(float, {'NL': 1, 'L': 5}) 1043 xout, yout, nfo = odesys.integrate(tout, conc0, par) 1044 assert nfo['success'] 1045 1046 1047@requires('sym', 'sympy', 'pycvodes') 1048def test_SymbolicSys__first_step_expr(): 1049 import sympy 1050 tend, k, y0 = 5, [1e23, 3], (.7, .0, .0) 1051 kwargs = dict(integrator='cvode', atol=1e-8, rtol=1e-8) 1052 factory = decay_dydt_factory(k) 1053 dep = sympy.symbols('y0 y1 y2', real=True) 1054 exprs = factory(k, dep) 1055 odesys = SymbolicSys(zip(dep, exprs), jac=True, first_step_expr=dep[0]*1e-30) 1056 xout, yout, info = odesys.integrate(tend, y0, **kwargs) 1057 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 1058 assert np.allclose(yout, ref, atol=10*kwargs['atol'], rtol=10*kwargs['rtol']) 1059 1060 1061@requires('sym', 'pygslodeiv2') 1062def test_SymbolicSys__from_callback__first_step_expr(): 1063 tend, k, y0 = 5, [1e23, 3], (.7, .0, .0) 1064 kwargs = dict(integrator='gsl', atol=1e-8, rtol=1e-8) 1065 factory = decay_dydt_factory(k) 1066 odesys = SymbolicSys.from_callback(factory, 3, first_step_factory=lambda x, y, p: y[0]*1e-30) 1067 xout, yout, info = odesys.integrate(tend, y0, **kwargs) 1068 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 1069 assert np.allclose(yout, ref, atol=10*kwargs['atol'], rtol=10*kwargs['rtol']) 1070 1071 1072@requires('sym', 'pycvodes') 1073def test_SymbolicSys__from_callback__first_step_expr__by_name(): 1074 kwargs = dict(integrator='cvode', atol=1e-8, rtol=1e-8) 1075 names = ['foo', 'bar', 'baz'] 1076 par_names = 'first second third'.split() 1077 odesys = SymbolicSys.from_callback( 1078 lambda x, y, p, be: { 1079 'foo': -p['first']*y['foo'], 1080 'bar': p['first']*y['foo'] - p['second']*y['bar'], 1081 'baz': p['second']*y['bar'] - p['third']*y['baz'] 1082 }, names=names, param_names=par_names, 1083 dep_by_name=True, par_by_name=True, 1084 first_step_factory=lambda x0, ic: 1e-30*ic['foo']) 1085 y0 = {'foo': .7, 'bar': 0, 'baz': 0} 1086 p = {'first': 1e23, 'second': 2, 'third': 3} 1087 result = odesys.integrate(5, y0, p, **kwargs) 1088 assert result.info['success'] 1089 ref = bateman_full([y0[k] for k in names], [p[k] for k in par_names], result.xout - result.xout[0], exp=np.exp) 1090 for i, k in enumerate(odesys.names): 1091 assert np.allclose(result.named_dep(k), ref[i], atol=10*kwargs['atol'], rtol=10*kwargs['rtol']) 1092 for k, v in p.items(): 1093 assert result.named_param(k) == v 1094 1095 1096@requires('sym', 'pyodeint') 1097def test_PartiallySolvedSystem__by_name(): 1098 k = [math.log(2)/(138.4*24*3600)] 1099 names = 'Po-210 Pb-206'.split() 1100 with pytest.raises(ValueError): 1101 odesys = SymbolicSys.from_callback(decay_dydt_factory({'Po-210': k[0]}, names=names), 1102 dep_by_name=True, par_by_name=True, names=names, param_names=names) 1103 odesys = SymbolicSys.from_callback(decay_dydt_factory({'Po-210': k[0]}, names=names), 1104 dep_by_name=True, names=names) 1105 1106 assert odesys.ny == 2 1107 partsys = PartiallySolvedSystem(odesys, lambda x0, y0, p0, be=None: { 1108 odesys['Pb-206']: y0[odesys['Pb-206']] + y0[odesys['Po-210']] - odesys['Po-210'] 1109 }) 1110 assert partsys.free_names == ['Po-210'] 1111 assert partsys.ny == 1 1112 assert (partsys['Pb-206'] - partsys.init_dep[partsys.names.index('Pb-206')] - 1113 partsys.init_dep[partsys.names.index('Po-210')] + odesys['Po-210']) == 0 1114 duration = 7*k[0] 1115 atol, rtol, forgive = 1e-9, 1e-9, 10 1116 y0 = [1e-20]*(len(k)+1) 1117 y0[0] = 1 1118 for system in (odesys, partsys): 1119 xout, yout, info = system.integrate(duration, y0, integrator='odeint', rtol=rtol, atol=atol) 1120 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 1121 assert np.allclose(yout, ref, rtol=rtol*forgive, atol=atol*forgive) 1122 assert yout.shape[1] == 2 1123 assert xout.shape[0] == yout.shape[0] 1124 assert yout.ndim == 2 and xout.ndim == 1 1125 1126 1127@requires('sym', 'pycvodes') 1128def test_PartiallySolvedSystem__by_name_2(): 1129 yn, pn = 'x y z'.split(), 'p q r'.split() 1130 odesys = _get_decay3_names(yn, pn) 1131 partsys = PartiallySolvedSystem(odesys, lambda x0, y0, p0, be: { 1132 odesys['x']: y0[odesys['x']]*be.exp(-p0['p']*(odesys.indep-x0)) 1133 }) 1134 y0 = [3, 2, 1] 1135 k = [3.5, 2.5, 1.5] 1136 1137 def _check(res): 1138 ref = np.array(bateman_full(y0, k, res.xout - res.xout[0], exp=np.exp)).T 1139 assert res.info['success'] 1140 assert np.allclose(res.yout, ref) 1141 args = [0, 1], dict(zip(yn, y0)), dict(zip(pn, k)) 1142 kwargs = dict(integrator='cvode') 1143 _check(odesys.integrate(*args, **kwargs)) 1144 _check(partsys.integrate(*args, **kwargs)) 1145 scaledsys = ScaledSys.from_other(partsys, dep_scaling=42, indep_scaling=17) 1146 _check(scaledsys.integrate(*args, **kwargs)) 1147 1148 1149@requires('sym') 1150def test_symmetricsys__invariants(): 1151 yn, pn = 'x y z'.split(), 'a b'.split() 1152 odesys = SymbolicSys.from_callback( 1153 lambda t, y, p: { 1154 'x': -p['a']*y['x'], 1155 'y': -p['b']*y['y'] + p['a']*y['x'], 1156 'z': p['b']*y['y'] 1157 }, names=yn, param_names=pn, dep_by_name=True, par_by_name=True, 1158 linear_invariants=[[1, 1, 1]], linear_invariant_names=['mass-conservation'], 1159 indep_name='t') 1160 assert odesys.linear_invariants.tolist() == [[1, 1, 1]] 1161 assert odesys.linear_invariant_names == ['mass-conservation'] 1162 assert odesys.nonlinear_invariants is None 1163 assert odesys.nonlinear_invariant_names is None 1164 1165 logexp = get_logexp() 1166 LogLogSys = symmetricsys(logexp, logexp) 1167 tsys = LogLogSys.from_other(odesys) 1168 assert tsys.linear_invariants is None 1169 assert tsys.linear_invariant_names is None 1170 assert len(tsys.nonlinear_invariants) == 1 1171 E = odesys.be.exp 1172 assert tsys.nonlinear_invariants[0] - sum(E(odesys[k]) for k in yn) == 0 1173 assert tsys.nonlinear_invariant_names == ['mass-conservation'] 1174 1175 1176@requires('sym', 'pycvodes') 1177@pycvodes_double 1178def test_SymbolicSys__roots(): 1179 def f(t, y): 1180 return [y[0]] 1181 1182 def roots(t, y, p, backend): 1183 return [y[0] - backend.exp(1)] 1184 odesys = SymbolicSys.from_callback(f, 1, roots_cb=roots) 1185 kwargs = dict(first_step=1e-14, atol=1e-14, rtol=1e-14, method='adams', integrator='cvode') 1186 xout, yout, info = odesys.integrate(2, [1], **kwargs) 1187 assert len(info['root_indices']) == 1 1188 assert np.min(np.abs(xout - 1)) < 1e-11 1189 1190 1191@requires('sym', 'pyodeint') 1192@pytest.mark.parametrize('method', ['bs', 'rosenbrock4']) 1193def test_SymbolicSys__reference_parameters_using_symbols(method): 1194 be = sym.Backend('sympy') 1195 x, p = map(be.Symbol, 'x p'.split()) 1196 symsys = SymbolicSys([(x, -p*x)], params=True) 1197 tout = [0, 1e-9, 1e-7, 1e-5, 1e-3, 0.1] 1198 for y_symb in [False, True]: 1199 for p_symb in [False, True]: 1200 xout, yout, info = symsys.integrate( 1201 tout, {x: 2} if y_symb else [2], {p: 3} if p_symb else [3], 1202 method=method, integrator='odeint', atol=1e-12, rtol=1e-12) 1203 assert np.allclose(yout[:, 0], 2*np.exp(-3*xout)) 1204 1205 1206@requires('sym', 'pygslodeiv2') 1207@pytest.mark.parametrize('method', ['rkck', 'rk4imp']) 1208def test_SymbolicSys__reference_parameters_using_symbols_from_callback(method): 1209 be = sym.Backend('sympy') 1210 k = be.Symbol('p') 1211 1212 def dydt(t, y): # external symbolic parameter 'k', should be allowed 1213 return [-k*y[0]] # even though reminiscent of global variables. 1214 1215 odesys1 = SymbolicSys.from_callback(dydt, 1, backend=be, params=True) 1216 odesys2 = SymbolicSys.from_callback(dydt, 1, backend=be, par_by_name=True, param_names=[], params=True) 1217 tout = [0, 1e-9, 1e-7, 1e-5, 1e-3, 0.1] 1218 for symsys in (odesys1, odesys2): 1219 for y_symb in [False, True]: 1220 for p_symb in [False, True]: 1221 xout, yout, info = symsys.integrate( 1222 tout, {symsys.dep[0]: 2} if y_symb else [2], {k: 3} if p_symb else [3], 1223 method=method, integrator='gsl', atol=1e-12, rtol=1e-12) 1224 assert xout.size > 4 1225 assert np.allclose(yout[:, 0], 2*np.exp(-3*xout)) 1226 1227 1228@requires('sym', 'pycvodes') 1229@pytest.mark.parametrize('scaled', [False, True]) 1230def test_PartiallySolvedSystem__from_linear_invariants(scaled): 1231 atol, rtol, forgive = 1e-11, 1e-11, 20 1232 k = [7., 3, 2] 1233 _ss = SymbolicSys.from_callback(decay_rhs, len(k)+1, len(k), 1234 linear_invariants=[[1]*(len(k)+1)], 1235 linear_invariant_names=['tot_amount']) 1236 if scaled: 1237 ss = ScaledSys.from_other(_ss, dep_scaling=1e3) 1238 else: 1239 ss = _ss 1240 1241 y0 = [0]*(len(k)+1) 1242 y0[0] = 1 1243 1244 def check_formulation(odesys): 1245 xout, yout, info = odesys.integrate( 1246 [0, 1], y0, k, integrator='cvode', atol=atol, rtol=rtol, nsteps=800) 1247 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 1248 assert np.allclose(yout, ref, rtol=rtol*forgive, atol=atol*forgive) 1249 1250 check_formulation(ss) 1251 1252 ps = PartiallySolvedSystem.from_linear_invariants(ss) 1253 assert ps.ny == ss.ny - 1 1254 check_formulation(ps) 1255 1256 ps2 = PartiallySolvedSystem(ss, lambda x0, y0, p0, be: { 1257 ss.dep[0]: y0[0]*be.exp(-p0[0]*(ss.indep-x0))}) 1258 assert ps2.ny == ss.ny - 1 1259 check_formulation(ps2) 1260 1261 1262@requires('sym', 'pyodeint') 1263def test_PartiallySolvedSystem__by_name__from_linear_invariants(): 1264 k = [math.log(2)/(138.4*24*3600)] 1265 names = 'Po-210 Pb-206'.split() 1266 odesys = SymbolicSys.from_callback( 1267 decay_dydt_factory({'Po-210': k[0]}, names=names), 1268 dep_by_name=True, names=names, linear_invariants=[[1, 1]]) 1269 assert odesys.ny == 2 1270 partsys1 = PartiallySolvedSystem.from_linear_invariants(odesys) 1271 partsys2 = PartiallySolvedSystem.from_linear_invariants(odesys, ['Pb-206']) 1272 partsys3 = PartiallySolvedSystem.from_linear_invariants(odesys, ['Po-210']) 1273 1274 assert partsys1.free_names in (['Po-210'], ['Pb-206']) 1275 assert partsys2.free_names == ['Po-210'] 1276 assert partsys3.free_names == ['Pb-206'] 1277 assert partsys1.ny == partsys2.ny == partsys3.ny == 1 1278 1279 assert (partsys2['Pb-206'] - partsys2.init_dep[partsys2.names.index('Pb-206')] - 1280 partsys2.init_dep[partsys2.names.index('Po-210')] + odesys['Po-210']) == 0 1281 duration = 7*k[0] 1282 atol, rtol, forgive = 1e-9, 1e-9, 10 1283 y0 = [1e-20]*(len(k)+1) 1284 y0[0] = 1 1285 for system in (odesys, partsys1, partsys2, partsys3): 1286 xout, yout, info = system.integrate(duration, y0, integrator='odeint', rtol=rtol, atol=atol) 1287 ref = np.array(bateman_full(y0, k+[0], xout - xout[0], exp=np.exp)).T 1288 assert np.allclose(yout, ref, rtol=rtol*forgive, atol=atol*forgive) 1289 assert yout.shape[1] == 2 1290 assert xout.shape[0] == yout.shape[0] 1291 assert yout.ndim == 2 and xout.ndim == 1 1292 1293 1294@requires('sym') 1295def test_SymbolicSys__indep_in_exprs(): 1296 def dydt(t, y, p): 1297 return [t*p[0]*y[0]] 1298 be = sym.Backend('sympy') 1299 t, y, p = map(be.Symbol, 't y p'.split()) 1300 odesys = SymbolicSys([(y, dydt(t, [y], [p])[0])], t, params=True) 1301 fout = odesys.f_cb(2, [3], [4]) 1302 assert len(fout) == 1 1303 assert abs(fout[0] - 2*3*4) < 1e-14 1304 1305 1306@requires('sym', 'pycvodes') 1307@pytest.mark.parametrize('idx', [0, 1, 2]) 1308def test_PartiallySolvedSystem__roots(idx): 1309 import sympy as sp 1310 t, x, y, z, p, q = sp.symbols('t x y z, p, q') 1311 odesys = SymbolicSys({x: -p*x, y: p*x - q*y, z: q*y}, t, params=(p, q), roots=([x - y], [x - z], [y - z])[idx]) 1312 _p, _q, tend = 7, 3, 0.7 1313 dep0 = {x: 1, y: 0, z: 0} 1314 ref = [0.11299628093544488, 0.20674119231833346, 0.3541828705348678] # determined in notebook: 1315 # test_symbolic__test_PartiallySolvedSystem__roots.ipynb 1316 1317 def check(odesys): 1318 res = odesys.integrate(tend, [dep0[k] for k in getattr(odesys, 'original_dep', odesys.dep)], (_p, _q), 1319 integrator='cvode', return_on_root=True) 1320 assert abs(res.xout[-1] - ref[idx]) < 1e-7 1321 1322 check(odesys) 1323 psys = PartiallySolvedSystem(odesys, lambda t0, xyz, par0: {x: xyz[odesys.dep.index(x)]*sp.exp(-p*(t-t0))}) 1324 check(psys) 1325 1326 1327@requires('sym', 'pycvodes') 1328@pytest.mark.parametrize('idx1,idx2,scaled,b2', product([0, 1, 2], [0, 1, 2], [True, False], [None, 0])) 1329def test_TransformedSys__roots(idx1, idx2, scaled, b2): 1330 def f(x, y, p): 1331 return [-p[0]*y[0], p[0]*y[0] - p[1]*y[1], p[1]*y[1]] 1332 1333 def roots(x, y): 1334 return ([y[0] - 3*y[1]], [y[0] - 3*y[2]], [3*y[1] - y[2]])[idx1] 1335 1336 if scaled: 1337 orisys = SymbolicSys.from_callback(f, 3, 2, roots_cb=roots) 1338 else: 1339 orisys = ScaledSys.from_callback(f, 3, 2, roots_cb=roots, dep_scaling=42) 1340 _p, _q, tend = 7, 3, 0.7 1341 dep0 = (1, .1, 0) 1342 ref = [0.02969588399749174, 0.1241509730780618, 0.6110670818418275] # determined in notebook: 1343 # test_symbolic__test_PartiallySolvedSystem__roots.ipynb 1344 1345 def check(odesys): 1346 res = odesys.integrate(tend, dep0, (_p, _q), 1347 integrator='cvode', return_on_root=True) 1348 assert abs(res.xout[-1] - ref[idx1]) < 6e-7 1349 1350 logexp = get_logexp(1, 1e-20, b2=None) 1351 LogLogSys = symmetricsys(logexp, logexp, check_transforms=False) 1352 1353 if idx2 == 0: # no need to rerun this code more than once 1354 check(orisys) 1355 loglog = LogLogSys.from_other(orisys) 1356 check(loglog) 1357 1358 psys1 = PartiallySolvedSystem(orisys, lambda t0, xyz, par0, be: { 1359 orisys.dep[0]: xyz[0]*be.exp(-par0[0]*(orisys.indep-t0))}) 1360 check(psys1) 1361 ploglog1 = LogLogSys.from_other(psys1) 1362 check(ploglog1) 1363 1364 psys2 = PartiallySolvedSystem(orisys, lambda t0, iv, p0: { 1365 orisys.dep[idx2]: iv[0] + iv[1] + iv[2] - sum(orisys.dep[j] for j in range(3) if j != idx2) 1366 }) 1367 ploglog2 = LogLogSys.from_other(psys2) 1368 check(ploglog2) 1369 1370 1371@requires('sym', 'sympy') 1372def test__group_invariants(): 1373 be = sym.Backend('sympy') 1374 x, y, z = symbs = be.symbols('x y z') 1375 1376 coeff1 = [3, 2, -1] 1377 expr1 = 3*x + 2*y - z 1378 lin, nonlin = _group_invariants([expr1], symbs, be) 1379 assert lin == [coeff1] 1380 assert nonlin == [] 1381 1382 expr2 = 3*x*x + 2*y - z 1383 lin, nonlin = _group_invariants([expr2], symbs, be) 1384 assert lin == [] 1385 assert nonlin == [expr2] 1386 1387 lin, nonlin = _group_invariants([expr1, expr2], symbs, be) 1388 assert lin == [coeff1] 1389 assert nonlin == [expr2] 1390 1391 lin, nonlin = _group_invariants([x + be.exp(y)], symbs, be) 1392 assert lin == [] 1393 assert nonlin == [x+be.exp(y)] 1394 1395 1396@requires('sym', 'pycvodes') 1397def test_SymbolicSys_as_autonomous(): 1398 import sympy 1399 1400 def rhs(t, y, p, backend=math): 1401 return [y[1], backend.sin(t)-p[0]*y[0]] 1402 odesys = SymbolicSys.from_callback(rhs, 2, 1) 1403 1404 def analytic(tout, init_y, p): 1405 t, (k,) = odesys.indep, odesys.params 1406 c1, c2 = sympy.symbols('c1 c2') 1407 sqk = sympy.sqrt(k) 1408 f = c1*sympy.cos(sqk*t) + c2*sympy.sin(sqk*t) + sympy.sin(t)/(k-1) 1409 dfdt = f.diff(t) 1410 t0 = tout[0] 1411 sol, = sympy.solve([f.subs(t, t0) - init_y[0], 1412 dfdt.subs(t, t0) - init_y[1]], 1413 [c1, c2], dict=True) 1414 sol[k] = p[0] 1415 exprs = [f.subs(sol), dfdt.subs(sol)] 1416 cb = sympy.lambdify([t], exprs) 1417 return np.array(cb(tout)).T 1418 1419 def integrate_and_check(system): 1420 init_y = [0, 0] 1421 p = [2] 1422 result = system.integrate([0, 80], init_y, p, integrator='cvode', nsteps=5000) 1423 yref = analytic(result.xout, init_y, p) 1424 assert np.all(result.yout - yref < 1.6e-5) 1425 1426 integrate_and_check(odesys) 1427 assert len(odesys.dep) == 2 1428 assert not odesys.autonomous_interface 1429 assert not odesys.autonomous_exprs 1430 asys = odesys.as_autonomous() 1431 integrate_and_check(asys) 1432 assert len(asys.dep) == 3 1433 assert not asys.autonomous_interface 1434 assert asys.autonomous_exprs 1435 1436 1437@requires('sym', 'pycvodes', 'quantities') 1438@pytest.mark.parametrize('auto', [False, True]) 1439def test_SymbolicSys_as_autonomous__to_arrays(auto): 1440 import quantities as pq 1441 1442 def rhs(t, y, p): 1443 k = t**p[0] 1444 return [-k*y[0], k*y[0]] 1445 1446 def analytic(tout, init_y, params): 1447 y0ref = init_y[0]*np.exp(-tout**(params[0]+1)/(params[0]+1)) 1448 return np.array([y0ref, init_y[0] - y0ref + init_y[1]]).T 1449 1450 odes = SymbolicSys.from_callback(rhs, 2, 1, to_arrays_callbacks=( 1451 lambda t: [_t.rescale(pq.s).magnitude for _t in t], 1452 lambda y: [_y.rescale(pq.molar).magnitude for _y in y], 1453 None 1454 )) 1455 odesys = odes.as_autonomous() if auto else odes 1456 for from_other in [False, True]: 1457 if from_other: 1458 odesys = SymbolicSys.from_other(odesys) 1459 result = odesys.integrate(4*pq.s, [5*pq.molar, 2*pq.molar], [3], integrator='cvode') 1460 ref = analytic(result.xout, result.yout[0, :], result.params) 1461 assert np.allclose(result.yout, ref, atol=1e-6) 1462 1463 1464@requires('sym', 'pycvodes') 1465def test_SymbolicSys_as_autonomous__linear_invariants(): 1466 def rhs(t, y, p): 1467 k = t**p[0] 1468 return [-k*y[0], k*y[0]] 1469 1470 def analytic(tout, init_y, params): 1471 y0ref = init_y[0]*np.exp(-tout**(params[0]+1)/(params[0]+1)) 1472 return np.array([y0ref, init_y[0] - y0ref + init_y[1]]).T 1473 1474 odes = SymbolicSys.from_callback(rhs, 2, 1, linear_invariants=[[1, 1]]) 1475 for odesys in [odes, odes.as_autonomous()]: 1476 result = odesys.integrate(4, [5, 2], [3], integrator='cvode') 1477 ref = analytic(result.xout, result.yout[0, :], result.params) 1478 assert np.allclose(result.yout, ref, atol=1e-6) 1479 1480 invar_viol = result.calc_invariant_violations() 1481 assert np.allclose(invar_viol, 0) 1482 1483 1484@requires('sym', 'pycvodes') 1485def test_SymbolicSys__by_name__as_autonomous(): 1486 def f(t, y, p): 1487 k = t**p['e'] 1488 return { 1489 'a': -k*y['a'], 1490 'b': +k*y['a'] 1491 } 1492 1493 def analytic(tout, init_y, p): 1494 y0ref = init_y[0]*np.exp(-tout**(p[0]+1)/(p[0]+1)) 1495 return np.array([y0ref, init_y[0] - y0ref + init_y[1]]).T 1496 1497 odes = SymbolicSys.from_callback( 1498 f, names='ab', param_names='e', dep_by_name=True, par_by_name=True, 1499 linear_invariants=[{'a': 1, 'b': 1}] 1500 ) 1501 1502 for odesys in [odes, odes.as_autonomous()]: 1503 result = odesys.integrate(3, {'a': 2, 'b': 1}, {'e': 2}, integrator='cvode') 1504 ref = analytic(result.xout, result.yout[0, :], result.params) 1505 assert np.allclose(result.yout, ref, atol=1e-6) 1506 1507 invar_viol = result.calc_invariant_violations() 1508 assert np.allclose(invar_viol, 0) 1509 1510 1511@requires('sym', 'pycvodes') 1512def test_SymbolicSys_as_autonomous__scaling(): 1513 1514 # 2 HNO2 -> H2O + NO + NO2; MassAction(EyringHS.fk('dH1', 'dS1')) 1515 # 2 NO2 -> N2O4; MassAction(EyringHS.fk('dH2', 'dS2')) 1516 # 1517 # HNO2 H2O NO NO2 N2O4 1518 def get_odesys(scaling=1): 1519 def rhs(t, y, p, backend=math): 1520 HNO2, H2O, NO, NO2, N2O4 = y 1521 dH1, dS1, dH2, dS2 = p 1522 R = 8.314 1523 T = 300 + 10*backend.sin(0.2*math.pi*t - math.pi/2) 1524 kB_h = 2.08366e10 1525 k1 = kB_h*T*backend.exp(dS1/R - dH1/(R*T))/scaling 1526 k2 = kB_h*T*backend.exp(dS2/R - dH2/(R*T))/scaling 1527 r1 = k1*HNO2**2 1528 r2 = k2*NO2**2 1529 return [-2*r1, r1, r1, r1 - 2*r2, r2] 1530 1531 return SymbolicSys.from_callback( 1532 rhs, 5, 4, names='HNO2 H2O NO NO2 N2O4'.split(), 1533 param_names='dH1 dS1 dH2 dS2'.split(), 1534 ) 1535 1536 def check(system, scaling=1): 1537 init_y = [1*scaling, 55*scaling, 0, 0, 0] 1538 p = [85e3, 10, 70e3, 20] 1539 return system.integrate(np.linspace(0, 60, 200), init_y, p, integrator='cvode', nsteps=5000) 1540 1541 def compare_autonomous(scaling): 1542 odesys = get_odesys(scaling) 1543 autsys = odesys.as_autonomous() 1544 copsys = SymbolicSys.from_other(autsys) 1545 res1 = check(odesys, scaling=scaling) 1546 res2 = check(autsys, scaling=scaling) 1547 res3 = check(copsys, scaling=scaling) 1548 assert np.allclose(res1.yout, res2.yout, atol=1e-6) 1549 assert np.allclose(res1.yout, res3.yout, atol=1e-6) 1550 1551 compare_autonomous(1) 1552 compare_autonomous(1000) 1553 1554 1555@requires('sym', 'pycvodes') 1556def test_SymbolicSys_as_autonomous__scaling__by_name(): 1557 1558 # 2 HNO2 -> H2O + NO + NO2; MassAction(EyringHS.fk('dH1', 'dS1')) 1559 # 2 NO2 -> N2O4; MassAction(EyringHS.fk('dH2', 'dS2')) 1560 # 1561 # HNO2 H2O NO NO2 N2O4 1562 def get_odesys(scaling=1): 1563 def rhs(t, y, p, backend=math): 1564 R = 8.314 1565 T = 300 + 10*backend.sin(0.2*math.pi*t - math.pi/2) 1566 kB_h = 2.08366e10 1567 k1 = kB_h*T*backend.exp(p['dS1']/R - p['dH1']/(R*T))/scaling # bimolecular => scaling**-1 1568 k2 = kB_h*T*backend.exp(p['dS2']/R - p['dH2']/(R*T))/scaling # bimolecular => scaling**-1 1569 r1 = k1*y['HNO2']**2 1570 r2 = k2*y['NO2']**2 1571 return {'HNO2': -2*r1, 'H2O': r1, 'NO': r1, 'NO2': r1 - 2*r2, 'N2O4': r2} 1572 1573 return SymbolicSys.from_callback( 1574 rhs, 5, 4, names='HNO2 H2O NO NO2 N2O4'.split(), 1575 param_names='dH1 dS1 dH2 dS2'.split(), dep_by_name=True, par_by_name=True, 1576 to_arrays_callbacks=( 1577 None, 1578 lambda y: [_y*scaling for _y in y], 1579 None 1580 )) 1581 1582 def check(system): 1583 init_y = {'HNO2': 1, 'H2O': 55, 'NO': 0, 'NO2': 0, 'N2O4': 0} 1584 p = {'dH1': 85e3, 'dS1': 10, 'dH2': 70e3, 'dS2': 20} 1585 return system.integrate(np.linspace(0, 60, 200), init_y, p, integrator='cvode', nsteps=5000) 1586 1587 def compare_autonomous(scaling): 1588 odesys = get_odesys(scaling) 1589 autsys = odesys.as_autonomous() 1590 copsys = SymbolicSys.from_other(autsys) 1591 res1 = check(odesys) 1592 res2 = check(autsys) 1593 res3 = check(copsys) 1594 assert np.allclose(res1.yout, res2.yout, atol=1e-6) 1595 assert np.allclose(res1.yout, res3.yout, atol=1e-6) 1596 1597 compare_autonomous(1) 1598 compare_autonomous(1000) 1599 1600 1601@requires('sym', 'pycvodes') 1602def _test_chained_parameter_variation(odesys_proc): 1603 def f(t, y, p): 1604 k = t**p['e'] 1605 return { 1606 'a': -k*y['a'], 1607 'b': +k*y['a'] 1608 } 1609 1610 def analytic(tout, init_y, p): 1611 y0ref = init_y['a']*np.exp(-tout**(p['e']+1)/(p['e']+1)) 1612 return np.array([y0ref, init_y['a'] - y0ref + init_y['b']]).T 1613 1614 odes = SymbolicSys.from_callback( 1615 f, names='ab', param_names='e', dep_by_name=True, par_by_name=True, 1616 linear_invariants=[{'a': 1, 'b': 1}] 1617 ) 1618 1619 durs = dur1, dur2 = 1.0, 3.0 1620 ndurs = 2 1621 e2 = 0 1622 1623 def _check_result(result): 1624 e1 = result.params[0] # only the initial parameter is stored in result.params 1625 mask1 = result.xout <= dur1 1626 mask2 = result.xout >= dur1 1627 x1 = result.xout[mask1] 1628 x2 = result.xout[mask2] - dur1 1629 y1 = result.yout[mask1, :] 1630 y2 = result.yout[mask2, :] 1631 ref1 = analytic(x1, dict(zip('ab', y1[0, :])), {'e': e1}) 1632 assert np.allclose(y1, ref1, atol=1e-6) 1633 1634 ref2 = analytic(x2, dict(zip('ab', y2[0, :])), {'e': e2}) 1635 assert np.allclose(y2, ref2, atol=1e-6) 1636 1637 invar_viol = result.calc_invariant_violations() 1638 assert np.allclose(invar_viol, 0) 1639 1640 for odesys in map(odesys_proc, [odes, odes.as_autonomous()]): 1641 for e1 in 2, 3: 1642 pars = {'e': [e1, e2]} 1643 result = chained_parameter_variation(odesys, durs, {'a': 2, 'b': 1}, pars, 1644 integrate_kwargs=dict(integrator='cvode')) 1645 _check_result(result) 1646 1647 for force_p in [False, True]: 1648 results = chained_parameter_variation( 1649 odesys, durs, {'a': [2, 3, 4], 'b': 1}, pars, 1650 integrate_kwargs=dict(integrator='cvode', force_predefined=force_p)) 1651 for res in results: 1652 assert res.info['mode'] == ['predefined']*ndurs if force_p else ['adaptive']*ndurs 1653 _check_result(res) 1654 1655 1656@requires('sym', 'pycvodes') 1657def test_chained_parameter_variation(): 1658 _test_chained_parameter_variation(lambda x: x) 1659 1660 1661@requires('sym', 'pycvodes') 1662def test_SymbolicSys_from_other_new_params(): 1663 odesys = _get_decay3() 1664 assert len(odesys.params) == 3 1665 p0, p1, p2 = odesys.params 1666 p3, = odesys.be.real_symarray('p', len(odesys.params)+1)[-1:] 1667 newode, extra = SymbolicSys.from_other_new_params(odesys, OrderedDict([ 1668 (p0, p3 - 1), 1669 (p1, p3 + 1) 1670 ]), (p3,)) 1671 assert len(newode.params) == 2 1672 tout = np.array([0.3, 0.4, 0.7, 0.9, 1.3, 1.7, 1.8, 2.1]) 1673 y0 = [7, 5, 2] 1674 p_vals = [5, 3] 1675 res1 = newode.integrate(tout, y0, p_vals) 1676 k = [p_vals[1] - 1, p_vals[1] + 1, p_vals[0]] 1677 ref1 = np.array(bateman_full(y0, k, res1.xout - res1.xout[0], exp=np.exp)).T 1678 assert np.allclose(res1.yout, ref1) 1679 orip = extra['recalc_params'](res1.xout, res1.yout, res1.params) 1680 assert np.allclose(orip, np.atleast_2d([3-1, 3+1, 5])) 1681 1682 1683@requires('sym', 'pycvodes') 1684def test_SymbolicSys_from_other_new_params_by_name(): 1685 yn, pn = 'x y z'.split(), 'p q r'.split() 1686 odesys = _get_decay3_names(yn, pn) 1687 assert len(odesys.params) == 3 1688 newode, extra = SymbolicSys.from_other_new_params_by_name( 1689 odesys, OrderedDict([ 1690 ('p', lambda x, y, p, backend: p['s'] - 1), 1691 ('q', lambda x, y, p, backend: p['s'] + 1) 1692 ]), ('s',) 1693 ) 1694 assert len(newode.params) == 2 1695 tout = np.array([0.3, 0.4, 0.7, 0.9, 1.3, 1.7, 1.8, 2.1]) 1696 y0 = {'x': 7, 'y': 5, 'z': 2} 1697 p_vals = {'r': 5, 's': 3} 1698 res1 = newode.integrate(tout, y0, p_vals) 1699 k = [p_vals['s'] - 1, p_vals['s'] + 1, p_vals['r']] 1700 ref1 = np.array(bateman_full( 1701 [y0[n] for n in newode.names], 1702 [p_vals['s'] - 1, p_vals['s'] + 1, p_vals['r']], 1703 res1.xout - res1.xout[0], exp=np.exp 1704 )).T 1705 assert np.allclose(res1.yout, ref1) 1706 orip = extra['recalc_params'](res1.xout, res1.yout, res1.params) 1707 assert np.allclose(orip, np.atleast_2d(k)) 1708