Native odesys ignores provided parameters if it is based on SymbolicSys created by chempy

[niklastoe]

I tried your example on native odesys. Providing different parameters leads to different concentrations as it should.

Now I noticed that if I construct the odesys using chempy, the results stay the same even if I change the parameters:

from chempy import ReactionSystem  # The rate constants below are arbitrary
from chempy.kinetics.ode import get_odesys
from collections import defaultdict
import numpy as np

from pyodesys.native import native_sys

rsys = ReactionSystem.from_string("""2 Fe+2 + H2O2 -> 2 Fe+3 + 2 OH-; 42
    2 Fe+3 + H2O2 -> 2 Fe+2 + O2 + 2 H+; 17
    H+ + OH- -> H2O; 1e10
    H2O -> H+ + OH-; 1e-4""")  # "[H2O]" = 1.0 (actually 55.4 at RT)
chempy_odesys, extra = get_odesys(rsys)

tout = sorted(np.concatenate((np.linspace(0, 23), np.logspace(-8, 1))))
c0 = defaultdict(float, {'Fe+2': 0.05, 'H2O2': 0.1, 'H2O': 1.0, 'H+': 1e-2, 'OH-': 1e-12})

print('Results from symbolic odesys')
result = chempy_odesys.integrate(tout, c0, atol=1e-12, rtol=1e-14)
np.array(result.at(tout))[-1,0]

chempy_native = native_sys['cvode'].from_other(chempy_odesys)

print('Results from native odesys without providing parameters')
result_chempy_native = chempy_native.integrate(tout, c0)
np.array(result_chempy_native.at(tout))[-1,0]

chempy_org_params = rsys.params()

print('Results from native odesys providing old parameters: ' + str(chempy_org_params))
result_chempy_native_org_params = chempy_native.integrate(tout, c0, chempy_org_params)
np.array(result_chempy_native_org_params.at(tout))[-1,0]

chempy_new_params = [4.2, 1.7, 1e5, 1e-2]
print('Results from native odesys providing new parameters: ' + str(chempy_new_params))
result_chempy_native_new_params = chempy_native.integrate(tout, c0, chempy_new_params)
np.array(result_chempy_native_new_params.at(tout))[-1,0]

print(chempy_org_params)
print(chempy_new_params)
print(np.array(chempy_org_params) / np.array(chempy_new_params))

The output:

Results from symbolic odesys

array([1.94487014e-02, 3.05512986e-02, 7.38097945e-12, 1.05023516e+00,
       4.44891879e-02, 2.01175814e-02, 2.05512986e-02])

Results from native odesys without providing parameters

INFO:pyodesys.native._base:In "/tmp/tmpGbkXsw_pycodeexport_pyodesys_NativeCvodeCode", executing:
"/usr/bin/g++ -c -std=c++11 -Wall -Wextra -fPIC -O2 -ffast-math -funroll-loops -fopenmp -o ./odesys_anyode.o -I/home/niklas/anaconda2/lib/python2.7/site-packages/numpy/core/include -I/home/niklas/anaconda2/lib/python2.7/site-packages/pyodesys/native/sources -I/home/niklas/anaconda2/lib/python2.7/site-packages/pycvodes/include odesys_anyode.cpp"
INFO:pyodesys.native._base:In "/tmp/tmpGbkXsw_pycodeexport_pyodesys_NativeCvodeCode", executing:
"/usr/bin/g++ -pthread -shared -std=c++11 -Wall -Wextra -fPIC -O2 -ffast-math -funroll-loops -fopenmp -o /tmp/tmpGbkXsw_pycodeexport_pyodesys_NativeCvodeCode/_cvode_wrapper.so -I/home/niklas/anaconda2/lib/python2.7/site-packages/numpy/core/include -I/home/niklas/anaconda2/lib/python2.7/site-packages/pyodesys/native/sources -I/home/niklas/anaconda2/lib/python2.7/site-packages/pycvodes/include odesys_anyode.o _cvode_wrapper.o -lsundials_cvodes -lsundials_nvecserial -lsundials_sunlinsollapackdense -lsundials_sunlinsollapackband -lopenblas -lpython2.7"

array([ 1.94486898e-02,  3.05513102e-02, -2.05513102e-02,  1.07078646e+00,
        4.44891934e-02,  2.01175757e-02, -7.42631650e-12])

Results from native odesys providing old parameters: [42, 17, 10000000000.0, 0.0001]

array([ 1.94486898e-02,  3.05513102e-02, -2.05513102e-02,  1.07078646e+00,
        4.44891934e-02,  2.01175757e-02, -7.42631650e-12])

Results from native odesys providing new parameters: [4.2, 1.7, 100000.0, 0.01]

array([ 1.94486898e-02,  3.05513102e-02, -2.05513102e-02,  1.07078646e+00,
        4.44891934e-02,  2.01175757e-02, -7.42631650e-12])

[42, 17, 10000000000.0, 0.0001]
[4.2, 1.7, 100000.0, 0.01]
[1.e+01 1.e+01 1.e+05 1.e-02]

There are minor numerical differences between the results from the symbolic odesys and the native one. In the latter, the ordering of species is slightly different and there are some negative concentrations, though. Can their absolute values be safely used or this there something wrong with them?
The provided parameters change nothing, though. Is there a way to make this work? Thanks in advance!

chempy 0.6.15
pyodesys 0.11.17
python 2.7.15

same behaviour with
chempy 0.7.6
pyodesys 0.12.4

[bjodah]

Thank you for reporting this. Indeed something's not quite right here. I'll try to find some time this weekend to take a closer look.

[bjodah]

OK I've found some time and looked closer at this. The API is not fully finished here (hence the conventional leading underscore in _native indicating its unofficial status). To ensure that you can change the parameters at integration you'll need to pass: include_params=False to get_odesys. And when I tried this quickly I also needed to name the parameters (that shouldn't be needed but perhaps it's good from a readbility standpoint anyway):

from chempy import ReactionSystem  # The rate constants below are arbitrary
from chempy.kinetics.ode import get_odesys
from collections import defaultdict
import numpy as np

from pyodesys.native import native_sys

PARAMETER_KEYS = 'k_Fe2p_H2O2 k_Fe3p_H2O2 k_Hp_OHm k_H2O'.split()
ORIGINAL_PARAMETERS = [42, 17, 1e10, 1e-4]
rsys = ReactionSystem.from_string("""2 Fe+2 + H2O2 -> 2 Fe+3 + 2 OH-; 'k_Fe2p_H2O2'
    2 Fe+3 + H2O2 -> 2 Fe+2 + O2 + 2 H+; 'k_Fe3p_H2O2'
    H+ + OH- -> H2O; 'k_Hp_OHm'
    H2O -> H+ + OH-; 'k_H2O'""")  # "[H2O]" = 1.0 (actually 55.4 at RT)
chempy_odesys, extra = get_odesys(rsys, include_params=False)

tout = sorted(np.concatenate((np.linspace(0, 23), np.logspace(-8, 1))))
c0 = defaultdict(float, {'Fe+2': 0.05, 'H2O2': 0.1, 'H2O': 1.0, 'H+': 1e-2, 'OH-': 1e-12})

print('Results from symbolic odesys')
result = chempy_odesys.integrate(tout, c0, dict(zip(PARAMETER_KEYS, ORIGINAL_PARAMETERS)), atol=1e-12, rtol=1e-14)
np.array(result.at(tout))[-1,0]

chempy_native = native_sys['cvode'].from_other(chempy_odesys)

print('Results from native odesys without providing parameters')
result_chempy_native = chempy_native.integrate(tout, c0, dict(zip(PARAMETER_KEYS, ORIGINAL_PARAMETERS)))
np.array(result_chempy_native.at(tout))[-1,0]

#chempy_org_params = rsys.params()

print('Results from native odesys providing old parameters: ' + str(chempy_org_params))
result_chempy_native_org_params = chempy_native.integrate(tout, c0, dict(zip(PARAMETER_KEYS, ORIGINAL_PARAMETERS)))
print(np.array(result_chempy_native_org_params.at(tout))[-1,0])

chempy_new_params = [4.2, 1.7, 1e5, 1e-2]
print('Results from native odesys providing new parameters: ' + str(chempy_new_params))
result_chempy_native_new_params = chempy_native.integrate(tout, c0, dict(zip(PARAMETER_KEYS, chempy_new_params)))
print(np.array(result_chempy_native_new_params.at(tout))[-1,0])

print(chempy_org_params)
print(chempy_new_params)
print(np.array(chempy_org_params) / np.array(chempy_new_params))

output

Results from symbolic odesys
INFO:pyodesys.native._base:In "/tmp/tmpvhonaiww_pycodeexport_pyodesys_NativeCvodeCode", executing:

"/usr/bin/g++ -c -std=c++11 -Wall -Wextra -fPIC -O2 -ffast-math -funroll-loops -fopenmp -o ./odesys_anyode.o -DPYCVODES_NO_KLU=0 -DPYCVODES_NO_LAPACK=0 -DANYODE_NO_LAPACK=0 -I/opt/cpython-3.7.3/lib/python3.7/site-packages/numpy/core/include -I/home/bjorn/vc/pyodesys/pyodesys/native/sources -I/home/bjorn/vc/pycvodes/pycvodes/include odesys_anyode.cpp"
Results from native odesys without providing parameters
INFO:pyodesys.native._base:In "/tmp/tmpvhonaiww_pycodeexport_pyodesys_NativeCvodeCode", executing:

"/usr/bin/g++ -pthread -shared -std=c++11 -Wall -Wextra -fPIC -O2 -ffast-math -funroll-loops -fopenmp -o /tmp/tmpvhonaiww_pycodeexport_pyodesys_NativeCvodeCode/_cvode_wrapper.cpython-37m-x86_64-linux-gnu.so -DPYCVODES_NO_KLU=0 -DPYCVODES_NO_LAPACK=0 -DANYODE_NO_LAPACK=0 -I/opt/cpython-3.7.3/lib/python3.7/site-packages/numpy/core/include -I/home/bjorn/vc/pyodesys/pyodesys/native/sources -I/home/bjorn/vc/pycvodes/pycvodes/include odesys_anyode.o _cvode_wrapper.o -lsundials_nvecserial -lsundials_cvodes -lsundials_sunlinsollapackdense -lsundials_sunlinsollapackband -lsundials_sunlinsolklu -lblas -llapack /opt/cpython-3.7.3/lib/libpython3.7m.so.1.0"
Results from native odesys providing old parameters: [42, 17, 10000000000.0, 0.0001]

[ 1.94486784e-02  3.05513216e-02 -2.05513216e-02  1.07078647e+00

4.44891923e-02  2.01175735e-02 -7.44019011e-12]

Results from native odesys providing new parameters: [4.2, 1.7, 100000.0, 0.01]

[2.76590861e-02 2.23409139e-02 8.30926852e-06 1.01166717e+00

8.71540680e-02 8.37737524e-04 1.23492232e-02]

[42, 17, 10000000000.0, 0.0001]

[4.2, 1.7, 100000.0, 0.01]

[1.e+01 1.e+01 1.e+05 1.e-02]

I tried using uppercase where I made the changes to highlight them.

If you are interested in polishing the API of chempy here (and/or pyodesys) I'd be more than happy to assist/review any proposal. Otherwise it might be a while since I don't currently have too much time so I've prioritized bug-fixes to adding new features.

[niklastoe]

Your fix solves my problem. Thank you for the quick reply, I really appreciate it. Btw, thanks for both chempy and pyodesys, they are very convenient to build kinetic models!

Personally, I think naming the parameters is very advantageous for readability and avoidance of disordering the parameters.

To be honest, I wouldn't know where to start in order to make this work with include_params=True. Raising a Warning during the creation of the native odesys if the parameters have been included by get_odesys should point other users in the right direction. Since include_params=False fixes the issue, it seems to me that it's more a matter of convenience. It is unclear to me though how instances of SymbolicSys (created with include_params=True and include_params=False resepectively) differ/what could trigger this warning. I compared their attributes and couldn't find anything.

[bjodah]

I'm happy to hear that you have found them useful! They aren't perfect but I've done my best to build a large test suite and (hopefully) decent API/documentation. Modelling kinetics and exploring different mechanisms is what I mainly use ChemPy for myself -- and it turns out that it has catched quite a few errors in some published models as well (unbalanced reactions, incorrect units etc.) which I miss even though I'm actively looking for them while transcribing from some article..

Yes, don't worry about adding warnings, but if you think of one which would be helpful, just let me know and we can try to figure out how to incorporate it. (One way that might work would be to store e.g. nparams in the SymbolicSys instance, and have the solve method check that we're passing the right number of parameters and emit a warning otherwise).