Native odesys ignores provided parameters if it is based on SymbolicSys created by chempy about pyodesys HOT 4 OPEN

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.

from pyodesys.

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))

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]

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.

from pyodesys.

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.

from pyodesys.

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).

from pyodesys.