Add option clipping to get_native

CHANGES.rst

@@ -1,3 +1,7 @@
+vX.Y.Z
+======
+- Support for new option "clip_to_bounds" in ``pyodesys``.
+
 v0.6.4
 ======
 - Enhancements for ``.units.is_unitless`` & ``.units.get_physical_quantity``

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);
 """
 
@@ -133,6 +133,32 @@ 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 odesys.clip_to_bounds:
+        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'] |= {"<stdlib.h>"}
+        kw['namespace_override']['p_y_preprocessing'] = """
+        double * y = (double*)(this->user_data);
+        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
         kw['namespace_override']['p_nroots'] = ' return 1; '

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(

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))
-    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 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[:, 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
@@ -162,20 +166,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
@@ -227,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