From 33baa4579ca6ead51d045a8b71f5fa986f189e2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bj=C3=B6rn=20Dahlgren?= Date: Thu, 1 Feb 2018 16:03:00 +0100 Subject: [PATCH 1/2] Add option clipping to get_native --- CHANGES.rst | 4 +++ chempy/kinetics/_native.py | 22 +++++++++++++- chempy/kinetics/tests/test__native.py | 43 ++++++++++++++------------- 3 files changed, 47 insertions(+), 22 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index c36dbdb6..889e128d 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,7 @@ +vX.Y.Z +====== +- New option "clippin" in ``chempy.kinetics.ode._native.get_native`` + v0.6.4 ====== - Enhancements for ``.units.is_unitless`` & ``.units.get_physical_quantity`` diff --git a/chempy/kinetics/_native.py b/chempy/kinetics/_native.py index 9a5f8720..6f6d79be 100644 --- a/chempy/kinetics/_native.py +++ b/chempy/kinetics/_native.py @@ -107,7 +107,7 @@ def _get_subst_comp(rsys, odesys, comp_keys, skip_keys): return subst_comp -def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False, conc_roots=None): +def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False, conc_roots=None, clipping=None): comp_keys = Substance.composition_keys(rsys.substances.values(), skip_keys=skip_keys) if isinstance(odesys, PartiallySolvedSystem): init_conc = '&m_p[%d]' % (len(odesys.params) - len(odesys.original_dep)) @@ -133,6 +133,26 @@ def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False raise ValueError("integrator '%s' does not support roots." % integrator) if odesys.roots is not None: raise ValueError("roots already set") + if clipping: + if 'p_constructor' not in ns_extend: + ns_extend['p_constructor'] = [] + ns_extend['p_constructor'] += [ + 'this->user_data = malloc(sizeof(double)*%(ny)d);' % dict(ny=odesys.ny) + ] + if 'p_destructor' not in ns_extend: + ns_extend['p_destructor'] = [] + ns_extend['p_destructor'] += [ + 'free(this->user_data);' + ] + if 'p_includes' not in ns_extend: + ns_extend['p_includes'] = set() + ns_extend['p_includes'] |= {""} + kw['namespace_override']['p_y_preprocessing'] = """ + double * y = (double*)(this->user_data); + for (int i=0; i<%(ny)d; ++i){ + y[i] = (y_[i] < 0) ? 0 : y_[i]; + } + """ % dict(ny=odesys.ny) if steady_state_root: assert conc_roots is None kw['namespace_override']['p_nroots'] = ' return 1; ' diff --git a/chempy/kinetics/tests/test__native.py b/chempy/kinetics/tests/test__native.py index 305270de..2607a6ea 100644 --- a/chempy/kinetics/tests/test__native.py +++ b/chempy/kinetics/tests/test__native.py @@ -96,13 +96,13 @@ def test_get_native__a_substance_no_composition(solve): if len(solve) > 0: from pyodesys.symbolic import PartiallySolvedSystem odesys = PartiallySolvedSystem(odesys, extra['linear_dependencies'](solve)) - odesys = get_native(rsys, odesys, 'gsl') - xout, yout, info = odesys.integrate(1, c0, atol=1e-15, rtol=1e-15, integrator='gsl') - c_reac = c0['H2O+'], c0['e-(aq)'] - H2O_ref = binary_rev(xout, 1e10, 1e-4, c0['H2O'], max(c_reac), min(c_reac)) - assert np.allclose(yout[:, odesys.names.index('H2O')], H2O_ref) - assert np.allclose(yout[:, odesys.names.index('H2O+')], c0['H2O+'] + c0['H2O'] - H2O_ref) - assert np.allclose(yout[:, odesys.names.index('e-(aq)')], c0['e-(aq)'] + c0['H2O'] - H2O_ref) + for odesys in [get_native(rsys, odesys, 'gsl', clipping=clipping) for clipping in [True, False]]: + xout, yout, info = odesys.integrate(1, c0, atol=1e-15, rtol=1e-15, integrator='gsl') + c_reac = c0['H2O+'], c0['e-(aq)'] + H2O_ref = binary_rev(xout, 1e10, 1e-4, c0['H2O'], max(c_reac), min(c_reac)) + assert np.allclose(yout[:, odesys.names.index('H2O')], H2O_ref) + assert np.allclose(yout[:, odesys.names.index('H2O+')], c0['H2O+'] + c0['H2O'] - H2O_ref) + assert np.allclose(yout[:, odesys.names.index('e-(aq)')], c0['e-(aq)'] + c0['H2O'] - H2O_ref) @pytest.mark.veryslow @@ -162,20 +162,21 @@ def test_get_native__conc_roots(): odesys, extra = get_odesys(rsys, include_params=False, unit_registry=u_reg) c0 = {'O3': 4.2e-3*M, 'O2': 0*M} cr = ['O2'] - native = get_native(rsys, odesys, 'cvode', conc_roots=cr) - tend = 1e5*u.s - params = {'k2': logspace_from_lin(1e-3/M/s, 1e3/M/s, 14)} - tgt_O2 = 1e-3*M - results = native.integrate(tend, c0, params, integrator='native', return_on_root=True, - special_settings=[unitless_in_registry(tgt_O2, u_reg)]) - assert len(results) == params['k2'].size - # dydt = -p*y**2 - # 1/y0 - 1/y = -2*pt - # t = 1/2/p*(1/y - 1/y0) - tgt_O3 = c0['O3'] - 2/3 * tgt_O2 - for r in results: - ref = rescale(1/2/r.named_param('k2')*(1/tgt_O3 - 1/c0['O3']), u.s) - assert allclose(r.xout[-1], ref, rtol=1e-6) + + for native in [get_native(rsys, odesys, 'cvode', conc_roots=cr, clipping=clipping) for clipping in [True, False]]: + tend = 1e5*u.s + params = {'k2': logspace_from_lin(1e-3/M/s, 1e3/M/s, 14)} + tgt_O2 = 1e-3*M + results = native.integrate(tend, c0, params, integrator='native', return_on_root=True, + special_settings=[unitless_in_registry(tgt_O2, u_reg)]) + assert len(results) == params['k2'].size + # dydt = -p*y**2 + # 1/y0 - 1/y = -2*pt + # t = 1/2/p*(1/y - 1/y0) + tgt_O3 = c0['O3'] - 2/3 * tgt_O2 + for r in results: + ref = rescale(1/2/r.named_param('k2')*(1/tgt_O3 - 1/c0['O3']), u.s) + assert allclose(r.xout[-1], ref, rtol=1e-6) @pytest.mark.veryslow From 571431b43746ebaa76bc0f894bc2ce13aa38becc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bj=C3=B6rn=20Dahlgren?= Date: Sat, 3 Feb 2018 15:09:57 +0100 Subject: [PATCH 2/2] [WIP] clip_to_bounds for ODE system --- CHANGES.rst | 2 +- chempy/kinetics/_native.py | 18 +++++--- chempy/kinetics/ode.py | 9 +++- chempy/kinetics/tests/test__native.py | 60 ++++++++++++++++++++------- 4 files changed, 65 insertions(+), 24 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 889e128d..93ede6a0 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,6 +1,6 @@ vX.Y.Z ====== -- New option "clippin" in ``chempy.kinetics.ode._native.get_native`` +- Support for new option "clip_to_bounds" in ``pyodesys``. v0.6.4 ====== diff --git a/chempy/kinetics/_native.py b/chempy/kinetics/_native.py index 6f6d79be..b95119bb 100644 --- a/chempy/kinetics/_native.py +++ b/chempy/kinetics/_native.py @@ -54,7 +54,7 @@ _first_step = """ m_upper_bounds = upper_conc_bounds(${init_conc}); - m_lower_bounds.resize(${odesys.ny}); + m_lower_bounds.resize(${odesys.ny}, 0); return m_rtol*std::min(get_dx_max(x, y), 1.0); """ @@ -107,7 +107,7 @@ def _get_subst_comp(rsys, odesys, comp_keys, skip_keys): return subst_comp -def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False, conc_roots=None, clipping=None): +def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False, conc_roots=None): comp_keys = Substance.composition_keys(rsys.substances.values(), skip_keys=skip_keys) if isinstance(odesys, PartiallySolvedSystem): init_conc = '&m_p[%d]' % (len(odesys.params) - len(odesys.original_dep)) @@ -133,7 +133,7 @@ def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False raise ValueError("integrator '%s' does not support roots." % integrator) if odesys.roots is not None: raise ValueError("roots already set") - if clipping: + if odesys.clip_to_bounds: if 'p_constructor' not in ns_extend: ns_extend['p_constructor'] = [] ns_extend['p_constructor'] += [ @@ -149,9 +149,15 @@ def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False ns_extend['p_includes'] |= {""} kw['namespace_override']['p_y_preprocessing'] = """ double * y = (double*)(this->user_data); - for (int i=0; i<%(ny)d; ++i){ - y[i] = (y_[i] < 0) ? 0 : y_[i]; - } + for (int i=0; i<%(ny)d; ++i) y[i] = y_[i]; + if (m_lower_bounds.size() > 0) + for (int i=0; i<%(ny)d; ++i) + if (y[i] < m_lower_bounds[i]) + y[i] = m_lower_bounds[i]; + if (m_upper_bounds.size() > 0) + for (int i=0; i<%(ny)d; ++i) + if (y[i] > m_upper_bounds[i]) + y[i] = m_upper_bounds[i]; """ % dict(ny=odesys.ny) if steady_state_root: assert conc_roots is None diff --git a/chempy/kinetics/ode.py b/chempy/kinetics/ode.py index 755e2ce9..07515d12 100644 --- a/chempy/kinetics/ode.py +++ b/chempy/kinetics/ode.py @@ -105,7 +105,8 @@ def dCdt_list(rsys, rates): def get_odesys(rsys, include_params=True, substitutions=None, SymbolicSys=None, unit_registry=None, - output_conc_unit=None, output_time_unit=None, cstr=False, constants=None, **kwargs): + output_conc_unit=None, output_time_unit=None, cstr=False, constants=None, + lower_bounds=False, upper_bounds=False, **kwargs): """ Creates a :class:`pyneqsys.SymbolicSys` from a :class:`ReactionSystem` The parameters passed to RateExpr will contain the key ``'time'`` corresponding to the @@ -130,6 +131,10 @@ def get_odesys(rsys, include_params=True, substitutions=None, SymbolicSys=None, See :func:`chempy.units.get_derived_units`. output_conc_unit : unit (Optional) output_time_unit : unit (Optional) + lower_bounds : bool + Generate expressions for lower bounds based on composition. + upper_bounds : bool + Generate expressions for upper bounds based on composition. cstr : bool Generate expressions for continuously stirred tank reactor. constants : module @@ -324,6 +329,8 @@ def reaction_rates(t, y, p, backend=math): latex_names=latex_names, param_names=param_names_for_odesys, linear_invariants=None if len(compo_vecs) == 0 else compo_vecs, linear_invariant_names=None if len(compo_names) == 0 else list(map(str, compo_names)), + lower_bounds_cb=lambda t, y, p, be: [0]*odesys.ny if lower_bounds else None, + upper_bounds_cb=lambda t, y, p, be: rsys.upper_conc_bounds(y, min_=be.Min, dtype=object), **kwargs) symbolic_ratexs = reaction_rates( diff --git a/chempy/kinetics/tests/test__native.py b/chempy/kinetics/tests/test__native.py index 2607a6ea..ddd82da4 100644 --- a/chempy/kinetics/tests/test__native.py +++ b/chempy/kinetics/tests/test__native.py @@ -3,6 +3,7 @@ from collections import defaultdict from functools import reduce +from itertools import product from operator import mul try: @@ -46,13 +47,7 @@ def decay_get_Cref(k, y0, tout): min(3, len(k)+1))]) -@pytest.mark.veryslow -@requires('pycvodes', 'pyodesys') -@pytest.mark.parametrize('solve', [(), ('CNO',)]) -def test_get_native__first_step(solve): - integrator = 'cvode' - forgive = 20 - +def _mk_decay_rsys(): def k(num): return "MassAction(unique_keys=('k%d',))" % num @@ -61,7 +56,15 @@ def k(num): "ONC -> NCO; %s" % k(2), "NCO -> CON; %s" % k(3) ] - rsys = ReactionSystem.from_string('\n'.join(lines), 'CNO ONC NCO CON') + return ReactionSystem.from_string('\n'.join(lines), 'CNO ONC NCO CON') + +@pytest.mark.veryslow +@requires('pycvodes', 'pyodesys') +@pytest.mark.parametrize('solve', [(), ('CNO',)]) +def test_get_native__first_step(solve): + integrator = 'cvode' + forgive = 20 + rsys = _mk_decay_rsys() odesys, extra = get_odesys(rsys, include_params=False) if len(solve) > 0: from pyodesys.symbolic import PartiallySolvedSystem @@ -88,21 +91,22 @@ def k(num): @pytest.mark.veryslow @requires('pygslodeiv2', 'pyodesys') -@pytest.mark.parametrize('solve', [(), ('H2O',)]) -def test_get_native__a_substance_no_composition(solve): +@pytest.mark.parametrize('solve_clip', list(product([(), ('H2O',)], [True, False]))) +def test_get_native__a_substance_no_composition(solve_clip): + solve, clip = solve_clip rsys = ReactionSystem.from_string('\n'.join(['H2O -> H2O+ + e-(aq); 1e-8', 'e-(aq) + H2O+ -> H2O; 1e10'])) - odesys, extra = get_odesys(rsys) + odesys, extra = get_odesys(rsys, clip_to_bounds=clip) c0 = {'H2O': 0, 'H2O+': 2e-9, 'e-(aq)': 3e-9} if len(solve) > 0: from pyodesys.symbolic import PartiallySolvedSystem odesys = PartiallySolvedSystem(odesys, extra['linear_dependencies'](solve)) - for odesys in [get_native(rsys, odesys, 'gsl', clipping=clipping) for clipping in [True, False]]: - xout, yout, info = odesys.integrate(1, c0, atol=1e-15, rtol=1e-15, integrator='gsl') + for odes in [odesys, get_native(rsys, odesys, 'gsl')]: + xout, yout, info = odes.integrate(1, c0, atol=1e-15, rtol=1e-15, integrator='gsl') c_reac = c0['H2O+'], c0['e-(aq)'] H2O_ref = binary_rev(xout, 1e10, 1e-4, c0['H2O'], max(c_reac), min(c_reac)) - assert np.allclose(yout[:, odesys.names.index('H2O')], H2O_ref) - assert np.allclose(yout[:, odesys.names.index('H2O+')], c0['H2O+'] + c0['H2O'] - H2O_ref) - assert np.allclose(yout[:, odesys.names.index('e-(aq)')], c0['e-(aq)'] + c0['H2O'] - H2O_ref) + assert np.allclose(yout[:, odes.names.index('H2O')], H2O_ref) + assert np.allclose(yout[:, odes.names.index('H2O+')], c0['H2O+'] + c0['H2O'] - H2O_ref) + assert np.allclose(yout[:, odes.names.index('e-(aq)')], c0['e-(aq)'] + c0['H2O'] - H2O_ref) @pytest.mark.veryslow @@ -228,3 +232,27 @@ def analytic_H(t, p, k, H0): ref_H2_uM = to_unitless(c0['H2'], u.micromolar) + to_unitless(c0['H'], u.micromolar)/2 + t_ul*p_ul/2 - ref_H_uM/2 assert np.allclose(to_unitless(result.named_dep('H'), u.micromolar), ref_H_uM) assert np.allclose(to_unitless(result.named_dep('H2'), u.micromolar), ref_H2_uM) + + +@pytest.mark.veryslow +@requires('pyodeint', 'pygslodeiv2', 'pycvodes', 'pyodesys') +@pytest.mark.parametrize('integrator', 'odeint gsl cvode'.split()) +def test_get_native_clipping(integrator): + forgive = 1.0 + rsys = _mk_decay_rsys() + odesys, extra = get_odesys(rsys, include_params=False) + native = get_native(rsys, odesys, integrator, clipping=True) + c0 = dict(zip(rsys.substances, [-5, 42.0, 17.0, 0.0])) + rate_coeffs = 13, 17, 21 + tend = 7 + args = (tend, c0, dict(zip(native.param_names, rate_coeffs))) + kwargs = dict(integrator=integrator, atol=1e-9, rtol=1e-9) + + for ode in [odesys, native]: + result = ode.integrate(*args, **kwargs) + ref = decay_get_Cref((0,) + rate_coeffs[1:], [c0[key] for key in ode.names], result.xout) + assert np.allclose(result.yout[:, :3], ref, + atol=kwargs['atol']*forgive, + rtol=kwargs['rtol']*forgive) + assert result.info['success'] + assert result.info['nfev'] > 10 and result.info['njev'] > 1