diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index d471c49a..acde638e 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -21,7 +21,7 @@ requirements: run: - python - numpy - - sym >=0.3.3 + - sym >=0.3.4 - sympy >=1.1.1,!=1.2 - scipy - matplotlib diff --git a/examples/_sparse_matrix.ipynb b/examples/_sparse_matrix.ipynb new file mode 100644 index 00000000..c1ce74a0 --- /dev/null +++ b/examples/_sparse_matrix.ipynb @@ -0,0 +1,131 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import sympy\n", + "from pyodesys.symbolic import SymbolicSys\n", + "from pyodesys.tests._robertson import get_ode_exprs\n", + "sympy.init_printing()\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "logc, logt, reduced, powsimp = False, False, False, False\n", + "odesys1 = SymbolicSys.from_callback(get_ode_exprs(logc, logt, reduced)[0],\n", + " 2 if reduced else 3, 6 if reduced else 3)\n", + "odesys1.exprs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys1._jac" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Converting to sparse representation of jacobian:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sparse = SymbolicSys.from_other(odesys1, sparse=True)\n", + "sparse.nnz, sparse._rowvals, sparse._colptrs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sparse._jac" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Converting back to dense representation of jacobian:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys2 = SymbolicSys.from_other(sparse, sparse=False)\n", + "assert odesys2.nnz == -1\n", + "odesys2._jac" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Solving the initial value problem numerically:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def solve_and_plot(odesys, integrate_kw, tout=[0, 1e11], y0=[1,0,0],\n", + " params=[0.04, 1e4, 3e7], ax=None):\n", + " result = odesys.integrate(tout, y0, params, nsteps=5000,\n", + " integrator='cvode', **integrate_kw)\n", + " result.plot(ax=ax, title_info=2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import pycvodes\n", + "has_klu = pycvodes.config['NO_KLU'] == '0'\n", + "systems = [odesys1, sparse, odesys2]\n", + "kws = [{}, dict(linear_solver='klu'), {}]\n", + "fig, axes = plt.subplots(1, len(systems), figsize=(5*len(systems), 4),\n", + " subplot_kw=dict(xscale='log', yscale='log'))\n", + "for odesys, ax, kw in zip(systems, axes, kws):\n", + " if kw.get('linear_solver', '') == 'klu' and not has_klu:\n", + " continue\n", + " solve_and_plot(odesys, kw, ax=ax)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/pyodesys/core.py b/pyodesys/core.py index 6b54fa57..c4d71c77 100644 --- a/pyodesys/core.py +++ b/pyodesys/core.py @@ -48,9 +48,15 @@ class ODESys(object): dependent variable (x). Signature is any of: - ``rhs(x, y[:]) -> f[:]`` - ``rhs(x, y[:], p[:]) -> f[:]`` - - ``rhs(x, y[:], p[:], backend=math) -> f[:]`` + - ``rhs(x, y[:], p[:], backend=math) -> f[:]``. jac : callback Jacobian matrix (dfdy). Required for implicit methods. + Signature should be either of: + - ``jac(x, y[:]) -> J`` + - ``jac(x, y[:], p[:]) -J``. + If ``nnz < 0``, ``J`` should be a 2D array-like object if ``nnz < 0`` + corresponding to a dense or banded jacobian (see also ``band``). + If ``nnz >= 0``, ``J`` should be an instance of ``scipy.sparse.csc_matrix``. dfdx : callback Signature ``dfdx(x, y[:], p[:]) -> out[:]`` (used by e.g. GSL) jtimes : callback @@ -100,6 +106,9 @@ class ODESys(object): Describes whether the independent variable appears in the rhs expressions. If set to ``True`` the underlying solver is allowed to shift the independent variable during integration. + nnz : int (default : -1) + Maximum number of non-zero entries in sparse jacobian. When + non-negative, jacobian is assumed to be dense or banded. Attributes ---------- @@ -115,6 +124,7 @@ class ODESys(object): For calculating the first step based on x0, y0 & p. roots_cb : callback nroots : int + nnz : int names : tuple of strings param_names : tuple of strings description : str @@ -151,7 +161,7 @@ def __init__(self, f, jac=None, dfdx=None, jtimes=None, first_step_cb=None, root par_by_name=False, latex_names=(), latex_param_names=(), latex_indep_name=None, taken_names=None, pre_processors=None, post_processors=None, append_iv=False, autonomous_interface=None, to_arrays_callbacks=None, autonomous_exprs=None, - _indep_autonomous_key=None, numpy=None, **kwargs): + _indep_autonomous_key=None, numpy=None, nnz=-1, **kwargs): self.f_cb = _ensure_4args(f) self.j_cb = _ensure_4args(jac) if jac is not None else None self.jtimes_cb = _ensure_4args(jtimes) if jtimes is not None else None @@ -163,6 +173,7 @@ def __init__(self, f, jac=None, dfdx=None, jtimes=None, first_step_cb=None, root if not band[0] >= 0 or not band[1] >= 0: raise ValueError("bands needs to be > 0 if provided") self.band = band + self.nnz = nnz self.names = tuple(names or ()) self.param_names = tuple(param_names or ()) self.indep_name = indep_name @@ -577,16 +588,29 @@ def _f(x, y, fout): if with_jacobian is None: raise Exception("Must pass with_jacobian") elif with_jacobian is True: - def _j(x, y, jout, dfdx_out=None, fy=None): - if len(_p) > 0: - jout[:, :] = np.asarray(self.j_cb(x, y, _p)) - else: - jout[:, :] = np.asarray(self.j_cb(x, y)) - if dfdx_out is not None: + if self.nnz >= 0: + def _j(x, y, data, colptrs, rowvals): + if len(_p) > 0: + J = self.j_cb(x, y, _p) + else: + J = self.j_cb(x, y) + J = J.asformat("csc") + data[:] = J.data + colptrs[:] = J.indptr + rowvals[:] = J.indices + + new_kwargs['nnz'] = self.nnz + else: + def _j(x, y, jout, dfdx_out=None, fy=None): if len(_p) > 0: - dfdx_out[:] = np.asarray(self.dfdx_cb(x, y, _p)) + jout[:, :] = np.asarray(self.j_cb(x, y, _p)) else: - dfdx_out[:] = np.asarray(self.dfdx_cb(x, y)) + jout[:, :] = np.asarray(self.j_cb(x, y)) + if dfdx_out is not None: + if len(_p) > 0: + dfdx_out[:] = np.asarray(self.dfdx_cb(x, y, _p)) + else: + dfdx_out[:] = np.asarray(self.dfdx_cb(x, y)) else: _j = None diff --git a/pyodesys/native/_base.py b/pyodesys/native/_base.py index a443d7fa..f8c8d947 100644 --- a/pyodesys/native/_base.py +++ b/pyodesys/native/_base.py @@ -125,11 +125,16 @@ def variables(self): all_invar = tuple(self.odesys.all_invariants()) ninvar = len(all_invar) jac = self.odesys.get_jac() + nnz = self.odesys.nnz all_exprs = self.odesys.exprs + all_invar - if jac is not False: + if jac is not False and nnz < 0: jac_dfdx = list(reduce(add, jac.tolist() + self.odesys.get_dfdx().tolist())) all_exprs += tuple(jac_dfdx) nj = len(jac_dfdx) + elif jac is not False and nnz >= 0: + jac_dfdx = list(reduce(add, jac.tolist())) + all_exprs += tuple(jac_dfdx) + nj = len(jac_dfdx) else: nj = 0 @@ -233,12 +238,18 @@ def _ccode(expr): 'cses': [(symb.name, _ccode(expr)) for symb, expr in jtimes_cses], 'exprs': list(map(_ccode, jtimes_exprs)) }, - p_jac=None if jac is False else { + p_jac_dense=None if jac is False or nnz >= 0 else { 'cses': [(symb.name, _ccode(expr)) for symb, expr in jac_cses], 'exprs': {(idx//ny, idx % ny): _ccode(expr) for idx, expr in enumerate(jac_exprs[:ny*ny])}, 'dfdt_exprs': list(map(_ccode, jac_exprs[ny*ny:])) }, + p_jac_sparse=None if jac is False or nnz < 0 else { + 'cses': [(symb.name, _ccode(expr)) for symb, expr in jac_cses], + 'exprs': list(map(_ccode, jac_exprs[:nj])), + 'colptrs': self.odesys._colptrs, + 'rowvals': self.odesys._rowvals + }, p_first_step=None if first_step is None else { 'cses': first_step_cses, 'expr': _ccode(first_step_exprs[0]), diff --git a/pyodesys/native/sources/odesys_anyode_iterative.hpp b/pyodesys/native/sources/odesys_anyode_iterative.hpp index aa2d462d..2eae06b4 100644 --- a/pyodesys/native/sources/odesys_anyode_iterative.hpp +++ b/pyodesys/native/sources/odesys_anyode_iterative.hpp @@ -19,6 +19,7 @@ namespace odesys_anyode { bool, Real_t, std::vector); int nrev=0; // number of calls to roots Index_t get_ny() const override; + Index_t get_nnz() const override; int get_nquads() const override; int get_nroots() const override; Real_t get_dx0(Real_t, const Real_t * const) override; @@ -39,6 +40,12 @@ namespace odesys_anyode { Real_t * const __restrict__ jac, long int ldim, Real_t * const __restrict__ dfdt=nullptr) override; + AnyODE::Status sparse_jac_csc(Real_t t, + const Real_t * const __restrict__ y, + const Real_t * const __restrict__ fy, + Real_t * const __restrict__ data, + Index_t * const __restrict__ colptrs, + Index_t * const __restrict__ rowvals) override; AnyODE::Status jtimes(const Real_t * const __restrict__ vec, Real_t * const __restrict__ out, Real_t t, diff --git a/pyodesys/native/sources/odesys_anyode_template.cpp b/pyodesys/native/sources/odesys_anyode_template.cpp index d8ff18e8..ba16c5f8 100644 --- a/pyodesys/native/sources/odesys_anyode_template.cpp +++ b/pyodesys/native/sources/odesys_anyode_template.cpp @@ -50,6 +50,7 @@ namespace odesys_anyode { std::vector); int nrev=0; // number of calls to roots indextype get_ny() const override; + indextype get_nnz() const override; int get_nquads() const override; int get_nroots() const override; realtype get_dx0(realtype, const realtype * const) override; @@ -70,6 +71,12 @@ namespace odesys_anyode { realtype * const __restrict__ jac, long int ldim, realtype * const __restrict__ dfdt=nullptr) override; + AnyODE::Status sparse_jac_csc(realtype t, + const realtype * const __restrict__ y, + const realtype * const __restrict__ fy, + realtype * const __restrict__ data, + indextype * const __restrict__ colptrs, + indextype * const __restrict__ rowvals) override; AnyODE::Status jtimes(const realtype * const __restrict__ vec, realtype * const __restrict__ out, realtype t, @@ -117,6 +124,10 @@ namespace odesys_anyode { return ${p_odesys.ny}; } + indextype OdeSys::get_nnz() const { + return ${p_odesys.nnz}; + } + int OdeSys::get_nquads() const { return 0; // Not implemeneted yet (cvodes from Sundials supports this) } @@ -227,9 +238,9 @@ namespace odesys_anyode { realtype * const __restrict__ jac, long int ldim, realtype * const __restrict__ dfdt) { - %if p_jac is not None: - %if order in p_jac: - ${p_jac[order]} + %if p_jac_dense is not None: + %if order in p_jac_dense: + ${p_jac_dense[order]} %else: // The AnyODE::ignore(...) calls below are used to generate code free from false compiler warnings. AnyODE::ignore(fy); // Currently we are not using fy (could be done through extensive pattern matching) @@ -237,14 +248,14 @@ namespace odesys_anyode { ${'AnyODE::ignore(y);' if (not any([yi in p_odesys.get_jac().free_symbols for yi in p_odesys.dep]) and not any([yi in p_odesys.get_dfdx().free_symbols for yi in p_odesys.dep])) else ''} - %for cse_token, cse_expr in p_jac['cses']: + %for cse_token, cse_expr in p_jac_dense['cses']: const auto ${cse_token} = ${cse_expr}; %endfor %for i_major in range(p_odesys.ny): %for i_minor in range(p_odesys.ny): <% - curr_expr = p_jac['exprs'][i_minor, i_major] if order == 'cmaj' else p_jac['exprs'][i_major, i_minor] + curr_expr = p_jac_dense['exprs'][i_minor, i_major] if order == 'cmaj' else p_jac_dense['exprs'][i_major, i_minor] if curr_expr == '0' and p_jacobian_set_to_zero_by_solver: continue %> jac[ldim*${i_major} + ${i_minor}] = ${curr_expr}; @@ -252,7 +263,7 @@ namespace odesys_anyode { %endfor if (dfdt){ - %for idx, expr in enumerate(p_jac['dfdt_exprs']): + %for idx, expr in enumerate(p_jac_dense['dfdt_exprs']): dfdt[${idx}] = ${expr}; %endfor } @@ -285,6 +296,41 @@ namespace odesys_anyode { %endif } + AnyODE::Status OdeSys::sparse_jac_csc(realtype x, + const realtype * const __restrict__ y, + const realtype * const __restrict__ fy, + realtype * const __restrict__ data, + indextype * const __restrict__ colptrs, + indextype * const __restrict__ rowvals) { + %if p_jac_sparse is not None: + AnyODE::ignore(fy); // Currently we are not using fy (could be done through extensive pattern matching) + ${'AnyODE::ignore(x);' if p_odesys.autonomous_exprs else ''} + ${'AnyODE::ignore(y);' if (not any([yi in p_odesys.get_jac().free_symbols for yi in p_odesys.dep]) and + not any([yi in p_odesys.get_dfdx().free_symbols for yi in p_odesys.dep])) else ''} + %for cse_token, cse_expr in p_jac_sparse['cses']: + const auto ${cse_token} = ${cse_expr}; + %endfor + + %for i in range(p_odesys.nnz): + data[${i}] = ${p_jac_sparse['exprs'][i]}; + %endfor + + %for i in range(p_odesys.nnz): + rowvals[${i}] = ${p_jac_sparse['rowvals'][i]}; + %endfor + + %for i in range(p_odesys.ny + 1): + colptrs[${i}] = ${p_jac_sparse['colptrs'][i]}; + %endfor + this->njev++; + return AnyODE::Status::success; + %else: + AnyODE::ignore(x); AnyODE::ignore(y); AnyODE::ignore(fy); + AnyODE::ignore(data); AnyODE::ignore(colptrs); AnyODE::ignore(rowvals); + return AnyODE::Status::unrecoverable_error; + %endif + } + realtype OdeSys::get_dx_max(realtype x, const realtype * const y) { %if p_get_dx_max is False: AnyODE::ignore(x); AnyODE::ignore(y); // avoid compiler warning about unused parameter. diff --git a/pyodesys/native/tests/test_cvode.py b/pyodesys/native/tests/test_cvode.py index 0fa74932..df86cee0 100644 --- a/pyodesys/native/tests/test_cvode.py +++ b/pyodesys/native/tests/test_cvode.py @@ -4,7 +4,7 @@ import numpy as np import pytest -from pyodesys.util import requires, pycvodes_double +from pyodesys.util import requires, pycvodes_double, pycvodes_klu from pyodesys.symbolic import SymbolicSys, PartiallySolvedSystem from ._tests import ( @@ -239,6 +239,28 @@ def f(t, y, p): assert np.allclose(yout[:, 0], ref) +@pytest.mark.veryslow +@requires('sym', 'pycvodes') +@pycvodes_klu +def test_sparse_jac_native_cvode(nu=0.01, k=1.0, m=1.0, x0=1.0, atol=1.0e-12, rtol=1.0e-12): + # Damped harmonic oscillator + w0 = (k/m)**0.5 + + def f(t, y, p): + return [y[1], -w0**2 * y[0] - nu * y[1]] + + odesys = NativeSys.from_callback(f, 2, 0, jac=True, sparse=True) + tout, yout, info = odesys.integrate(100, [x0, 0], integrator='cvode', linear_solver='klu', + with_jacobian=True, atol=atol, rtol=rtol, nsteps=20000) + + w = (w0**2 - nu**2/4.0)**0.5 + a = (x0**2 + (nu * x0/2)**2/w**2)**0.5 + phi = np.arctan(nu * x0 / (2 * x0 * w)) + ref = a * np.exp(-nu * tout/2) * np.cos(w * tout - phi) + assert info['njev'] > 0 + assert np.allclose(yout[:, 0], ref) + + @pytest.mark.slow @requires('pycvodes', 'sympy') def test_render_native_cse_regression(): diff --git a/pyodesys/symbolic.py b/pyodesys/symbolic.py index 7773a73b..bfb590ca 100644 --- a/pyodesys/symbolic.py +++ b/pyodesys/symbolic.py @@ -177,6 +177,10 @@ class SymbolicSys(ODESys): When ``True`` construct using ``be.Symbol``. See also :attr:`init_indep`. init_dep : tuple of Symbols, ``True`` or ``None`` When ``True`` construct using ``be.Symbol``. See also :attr:`init_dep`. + sparse : bool (default ``False``) + When ``True``, jacobian will be derived and stored (if ``jac = True``) + in compressed sparse column (CSC) format and the jacobian callback + will return an instance of ``scipy.sparse.csc_matrix``. \*\*kwargs: See :py:class:`ODESys` @@ -228,7 +232,7 @@ def __init__(self, dep_exprs, indep=None, params=None, jac=True, dfdx=True, jtimes=False, first_step_expr=None, roots=None, backend=None, lower_bounds=None, upper_bounds=None, linear_invariants=None, nonlinear_invariants=None, linear_invariant_names=None, nonlinear_invariant_names=None, steady_state_root=False, - init_indep=None, init_dep=None, **kwargs): + init_indep=None, init_dep=None, sparse=False, **kwargs): self.dep, self.exprs = zip(*dep_exprs.items()) if isinstance(dep_exprs, dict) else zip(*dep_exprs) self.indep = indep if params is True or params is None: @@ -278,6 +282,7 @@ def __init__(self, dep_exprs, indep=None, params=None, jac=True, dfdx=True, if _param_names is True: kwargs['param_names'] = [p.name for p in self.params] + self.sparse = sparse # needed by get_j_ty_callback self.band = kwargs.get('band', None) # needed by get_j_ty_callback # bounds needed by get_f_ty_callback: self.lower_bounds = None if lower_bounds is None else np.array(lower_bounds)*np.ones(self.ny) @@ -294,6 +299,8 @@ def __init__(self, dep_exprs, indep=None, params=None, jac=True, dfdx=True, autonomous_exprs=_is_autonomous(self.indep, self.exprs), **kwargs) + if self.sparse: + self.nnz = len(self._rowvals) self.linear_invariants = linear_invariants self.nonlinear_invariants = nonlinear_invariants self.linear_invariant_names = linear_invariant_names or None @@ -630,11 +637,13 @@ def autonomous_post_processor(x, y, p): def get_jac(self): """ Derives the jacobian from ``self.exprs`` and ``self.dep``. """ if self._jac is True: - if self.band is None: + if self.sparse is True: + self._jac, self._colptrs, self._rowvals = self.be.sparse_jacobian_csc(self.exprs, self.dep) + elif self.band is not None: # Banded + self._jac = self.be.banded_jacobian(self.exprs, self.dep, *self.band) + else: f = self.be.Matrix(1, self.ny, self.exprs) self._jac = f.jacobian(self.be.Matrix(1, self.ny, self.dep)) - else: # Banded - self._jac = self.be.banded_jacobian(self.exprs, self.dep, *self.band) elif self._jac is False: return False @@ -704,7 +713,17 @@ def get_j_ty_callback(self): j_exprs = self.get_jac() if j_exprs is False: return None - return self._callback_factory(j_exprs) + cb = self._callback_factory(j_exprs) + if self.sparse: + from scipy.sparse import csc_matrix + + def sparse_cb(x, y, p=()): + data = cb(x, y, p).flatten() + return csc_matrix((data, self._rowvals, self._colptrs)) + + return sparse_cb + else: + return cb def get_dfdx_callback(self): """ Generate a callback for evaluating derivative of ``self.exprs`` """ diff --git a/pyodesys/tests/test_core.py b/pyodesys/tests/test_core.py index 3536ee80..8ede8d7b 100644 --- a/pyodesys/tests/test_core.py +++ b/pyodesys/tests/test_core.py @@ -8,7 +8,7 @@ import numpy as np from .. import ODESys, OdeSys, chained_parameter_variation # OdeSys deprecated from ..core import integrate_chained -from ..util import requires +from ..util import requires, pycvodes_klu @requires('pycvodes') @@ -231,6 +231,17 @@ def sine_jac(t, y, p): return Jmat +def sine_jac_sparse(t, y, p): + from scipy.sparse import csc_matrix + k = p[0] + Jmat = np.zeros((2, 2)) + Jmat[0, 0] = 0 + Jmat[0, 1] = 1 + Jmat[1, 0] = -k**2 + Jmat[1, 1] = 0 + return csc_matrix(Jmat) + + def sine_dfdt(t, y, p): return [0, 0] @@ -308,6 +319,22 @@ def test_integrate_multiple_adaptive__pygslodeiv2(): integrator='gsl', method='bsimp') +@requires('pycvodes', 'scipy') +@pycvodes_klu +def test_sparse_jac(): + odesys = ODESys(sine, sine_jac_sparse, nnz=2) + A, k = 2, 3 # np.array(3) does not support len() + xout, yout, info = odesys.integrate(np.linspace(0, 1), [0, A*k], [k], + integrator='cvode', linear_solver='klu') + assert info['success'] + ref = [ + A*np.sin(k*(xout - xout[0])), + A*np.cos(k*(xout - xout[0]))*k + ] + assert np.allclose(yout[:, 0], ref[0], atol=1e-5, rtol=1e-5) + assert np.allclose(yout[:, 1], ref[1], atol=1e-5, rtol=1e-5) + + def decay_factory(ny): def f(t, y, p): diff --git a/pyodesys/tests/test_symbolic.py b/pyodesys/tests/test_symbolic.py index 562129a0..5af1200f 100644 --- a/pyodesys/tests/test_symbolic.py +++ b/pyodesys/tests/test_symbolic.py @@ -21,7 +21,7 @@ from .. import ODESys from ..core import integrate_auto_switch, chained_parameter_variation from ..symbolic import SymbolicSys, ScaledSys, symmetricsys, PartiallySolvedSystem, get_logexp, _group_invariants -from ..util import requires, pycvodes_double +from ..util import requires, pycvodes_double, pycvodes_klu from .bateman import bateman_full # analytic, never mind the details from .test_core import vdp_f from . import _cetsa @@ -136,6 +136,18 @@ def test_SymbolicSys__jacobian_singular(): assert odesys.jacobian_singular() +@requires('sym', 'pycvodes') +@pycvodes_klu +def test_SymbolicSys_jacobian_sparse(): + k = (4, 3) + y0 = (5, 4, 2) + odesys = SymbolicSys.from_callback(decay_dydt_factory(k), len(k) + 1, sparse=True) + xout, yout, info = odesys.integrate([0, 2], y0, integrator='cvode', linear_solver='klu', + atol=1e-12, rtol=1e-12, nsteps=1000) + ref = np.array(bateman_full(y0, k+(0,), xout - xout[0], exp=np.exp)).T + assert np.allclose(yout, ref, rtol=1e-10, atol=1e-10) + + @requires('sym', 'pycvodes') def test_TransformedSys_liny_linx(): _test_TransformedSys(idty2, idty2, 1e-11, 1e-11, 0, 15) diff --git a/pyodesys/util.py b/pyodesys/util.py index 3d1ef70f..86e9a2db 100644 --- a/pyodesys/util.py +++ b/pyodesys/util.py @@ -9,6 +9,7 @@ from pkg_resources import parse_requirements, parse_version import numpy as np +import pytest def stack_1d_on_left(x, y): @@ -225,7 +226,6 @@ def __init__(self, *reqs): self.incomp.append(str(req)) def __call__(self, cb): - import pytest r = 'Unfulfilled requirements.' if self.missing: r += " Missing modules: %s." % ', '.join(self.missing) @@ -235,7 +235,6 @@ def __call__(self, cb): def pycvodes_double(cb): - import pytest try: from pycvodes._config import env prec = env.get('SUNDIALS_PRECISION', 'double') @@ -245,6 +244,16 @@ def pycvodes_double(cb): return pytest.mark.skipif(prec != "double", reason=r)(cb) +def pycvodes_klu(cb): + try: + from pycvodes._config import env + no_klu = env['NO_KLU'] == '1' + except (ModuleNotFoundError, ImportError): + no_klu = True + return pytest.mark.skipif(no_klu, + reason="Sparse jacobian tests require pycvodes and sundials with KLU enabled.")(cb) + + class MissingImport(object): def __init__(self, modname, exc): diff --git a/setup.py b/setup.py index 38adbb56..21d31ad4 100755 --- a/setup.py +++ b/setup.py @@ -101,7 +101,7 @@ def _path_under_setup(*args): license=license, packages=[pkg_name] + submodules + tests, include_package_data=True, - install_requires=['numpy>=1.8.0', 'pytest>=2.9.2', 'scipy>=0.19.1', 'sym>=0.3.3', + install_requires=['numpy>=1.8.0', 'pytest>=2.9.2', 'scipy>=0.19.1', 'sym>=0.3.4', 'sympy>=1.1.1,!=1.2', 'matplotlib>=2.0.2', 'jupyter'], extras_require=extras_req )