From 41ff9b918bdf2c9dc5e058c7d02f44218a479730 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 1 Mar 2022 23:21:49 +0100 Subject: [PATCH 01/24] Add support for observable-dependent sigmas Standard deviations may now depend on observables. Closes #609 --- include/amici/abstract_model.h | 6 ++- matlab/@amifun/getArgs.m | 4 +- models/model_calvetti/model_calvetti.h | 12 +++--- .../model_calvetti/model_calvetti_sigmay.cpp | 2 +- models/model_dirac/model_dirac.h | 12 +++--- models/model_dirac/model_dirac_sigmay.cpp | 2 +- models/model_events/model_events.h | 12 +++--- models/model_events/model_events_sigmay.cpp | 2 +- .../model_jakstat_adjoint.h | 16 +++---- .../model_jakstat_adjoint_dsigmaydp.cpp | 2 +- .../model_jakstat_adjoint_sigmay.cpp | 2 +- .../model_jakstat_adjoint_o2.h | 16 +++---- .../model_jakstat_adjoint_o2_dsigmaydp.cpp | 2 +- .../model_jakstat_adjoint_o2_sigmay.cpp | 2 +- .../model_nested_events/model_nested_events.h | 12 +++--- .../model_nested_events_sigmay.cpp | 2 +- models/model_neuron/model_neuron.h | 12 +++--- models/model_neuron/model_neuron_sigmay.cpp | 2 +- models/model_neuron_o2/model_neuron_o2.h | 12 +++--- .../model_neuron_o2_sigmay.cpp | 2 +- models/model_robertson/model_robertson.h | 12 +++--- .../model_robertson_sigmay.cpp | 2 +- models/model_steadystate/model_steadystate.h | 12 +++--- .../model_steadystate_sigmay.cpp | 2 +- python/amici/ode_export.py | 4 +- python/tests/test_sbml_import.py | 42 +++++++++++++++++++ src/abstract_model.cpp | 4 +- src/model.cpp | 8 +++- 28 files changed, 135 insertions(+), 85 deletions(-) diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h index 4ddacdb26b..99202ca0f1 100644 --- a/include/amici/abstract_model.h +++ b/include/amici/abstract_model.h @@ -561,9 +561,10 @@ class AbstractModel { * @param t current time * @param p parameter vector * @param k constant vector + * @param y model output at timepoint t */ virtual void fsigmay(realtype *sigmay, const realtype t, const realtype *p, - const realtype *k); + const realtype *k, const realtype *y); /** * @brief Model-specific implementation of fsigmay @@ -574,7 +575,8 @@ class AbstractModel { * @param ip sensitivity index */ virtual void fdsigmaydp(realtype *dsigmaydp, const realtype t, - const realtype *p, const realtype *k, int ip); + const realtype *p, const realtype *k, + const realtype *y, int ip); /** * @brief Model-specific implementation of fsigmaz diff --git a/matlab/@amifun/getArgs.m b/matlab/@amifun/getArgs.m index a8e3b69c9c..e59c10173c 100644 --- a/matlab/@amifun/getArgs.m +++ b/matlab/@amifun/getArgs.m @@ -80,9 +80,9 @@ case 'dxdotdp' this.argstr = ['(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip' dx ', const realtype *w, const realtype *dwdp)']; case 'sigma_y' - this.argstr = '(double *sigmay, const realtype t, const realtype *p, const realtype *k)'; + this.argstr = '(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y)'; case 'dsigma_ydp' - this.argstr = '(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip)'; + this.argstr = '(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip)'; case 'sigma_z' this.argstr = '(double *sigmaz, const realtype t, const realtype *p, const realtype *k)'; case 'dsigma_zdp' diff --git a/models/model_calvetti/model_calvetti.h b/models/model_calvetti/model_calvetti.h index 52675c57ac..41fe28a957 100644 --- a/models/model_calvetti/model_calvetti.h +++ b/models/model_calvetti/model_calvetti.h @@ -1,6 +1,6 @@ #ifndef _amici_model_calvetti_h #define _amici_model_calvetti_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -22,7 +22,7 @@ extern void dJydy_model_calvetti(double *dJydy, const int iy, const realtype *p, extern void dwdx_model_calvetti(realtype *dwdx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl); extern void dydx_model_calvetti(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx); extern void root_model_calvetti(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx); -extern void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void w_model_calvetti(realtype *w, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); extern void x0_model_calvetti(realtype *x0, const realtype t, const realtype *p, const realtype *k); extern void xdot_model_calvetti(realtype *xdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *w); @@ -67,7 +67,7 @@ class Model_model_calvetti : public amici::Model_DAE { amici::Model* clone() const override { return new Model_model_calvetti(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override { JSparse_model_calvetti(JSparse, t, x, p, k, h, cj, dx, w, dwdx); @@ -125,7 +125,7 @@ class Model_model_calvetti : public amici::Model_DAE { void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -161,8 +161,8 @@ class Model_model_calvetti : public amici::Model_DAE { void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_calvetti(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_calvetti(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_calvetti/model_calvetti_sigmay.cpp b/models/model_calvetti/model_calvetti_sigmay.cpp index a970a9b280..690d1a01a9 100644 --- a/models/model_calvetti/model_calvetti_sigmay.cpp +++ b/models/model_calvetti/model_calvetti_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_calvetti{ -void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = 1.0; sigmay[1] = 1.0; sigmay[2] = 1.0; diff --git a/models/model_dirac/model_dirac.h b/models/model_dirac/model_dirac.h index 03639ddef6..a379b3ea24 100644 --- a/models/model_dirac/model_dirac.h +++ b/models/model_dirac/model_dirac.h @@ -1,6 +1,6 @@ #ifndef _amici_model_dirac_h #define _amici_model_dirac_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -23,7 +23,7 @@ extern void deltax_model_dirac(double *deltax, const realtype t, const realtype extern void dxdotdp_model_dirac(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp); extern void dydx_model_dirac(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx); extern void root_model_dirac(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); -extern void sigmay_model_dirac(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_dirac(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void stau_model_dirac(double *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl, const realtype *sx, const int ip, const int ie); extern void xdot_model_dirac(realtype *xdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w); extern void y_model_dirac(double *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w); @@ -67,7 +67,7 @@ class Model_model_dirac : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_dirac(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_dirac(JSparse, t, x, p, k, h, w, dwdx); @@ -123,7 +123,7 @@ class Model_model_dirac : public amici::Model_ODE { void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -159,8 +159,8 @@ class Model_model_dirac : public amici::Model_ODE { void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_dirac(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_dirac(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_dirac/model_dirac_sigmay.cpp b/models/model_dirac/model_dirac_sigmay.cpp index 97538bff2a..9b5782f79f 100644 --- a/models/model_dirac/model_dirac_sigmay.cpp +++ b/models/model_dirac/model_dirac_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_dirac{ -void sigmay_model_dirac(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_dirac(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = 1.0; } diff --git a/models/model_events/model_events.h b/models/model_events/model_events.h index 2dc99f9399..c83b1dbe62 100644 --- a/models/model_events/model_events.h +++ b/models/model_events/model_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_events_h #define _amici_model_events_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -32,7 +32,7 @@ extern void dydx_model_events(double *dydx, const realtype t, const realtype *x, extern void dzdx_model_events(double *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h); extern void root_model_events(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); extern void rz_model_events(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h); -extern void sigmay_model_events(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_events(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void sigmaz_model_events(double *sigmaz, const realtype t, const realtype *p, const realtype *k); extern void srz_model_events(double *srz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *sx, const int ip); extern void stau_model_events(double *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl, const realtype *sx, const int ip, const int ie); @@ -81,7 +81,7 @@ class Model_model_events : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_events(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_events(JSparse, t, x, p, k, h, w, dwdx); @@ -143,7 +143,7 @@ class Model_model_events : public amici::Model_ODE { drzdx_model_events(drzdx, ie, t, x, p, k, h); } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -182,8 +182,8 @@ class Model_model_events : public amici::Model_ODE { rz_model_events(rz, ie, t, x, p, k, h); } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_events(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_events(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_events/model_events_sigmay.cpp b/models/model_events/model_events_sigmay.cpp index 42fac3dfbf..87bd4ffb30 100644 --- a/models/model_events/model_events_sigmay.cpp +++ b/models/model_events/model_events_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_events{ -void sigmay_model_events(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_events(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = 1.0; } diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint.h b/models/model_jakstat_adjoint/model_jakstat_adjoint.h index 69bd1eadc2..5ec0de560d 100644 --- a/models/model_jakstat_adjoint/model_jakstat_adjoint.h +++ b/models/model_jakstat_adjoint/model_jakstat_adjoint.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_h #define _amici_model_jakstat_adjoint_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -18,13 +18,13 @@ extern void JSparse_model_jakstat_adjoint(SUNMatrixContent_Sparse JSparse, const extern void Jy_model_jakstat_adjoint(double *nllh, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); extern void dJydsigma_model_jakstat_adjoint(double *dJydsigma, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); extern void dJydy_model_jakstat_adjoint(double *dJydy, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); -extern void dsigmaydp_model_jakstat_adjoint(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip); +extern void dsigmaydp_model_jakstat_adjoint(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip); extern void dwdp_model_jakstat_adjoint(realtype *dwdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl, const realtype *stcl); extern void dwdx_model_jakstat_adjoint(realtype *dwdx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl); extern void dxdotdp_model_jakstat_adjoint(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp); extern void dydp_model_jakstat_adjoint(double *dydp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp); extern void dydx_model_jakstat_adjoint(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx); -extern void sigmay_model_jakstat_adjoint(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_jakstat_adjoint(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void sx0_model_jakstat_adjoint(realtype *sx0, const realtype t,const realtype *x0, const realtype *p, const realtype *k, const int ip); extern void w_model_jakstat_adjoint(realtype *w, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); extern void x0_model_jakstat_adjoint(realtype *x0, const realtype t, const realtype *p, const realtype *k); @@ -70,7 +70,7 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_jakstat_adjoint(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_jakstat_adjoint(JSparse, t, x, p, k, h, w, dwdx); @@ -124,8 +124,8 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE { void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { - dsigmaydp_model_jakstat_adjoint(dsigmaydp, t, p, k, ip); + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { + dsigmaydp_model_jakstat_adjoint(dsigmaydp, t, p, k, y, ip); } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -163,8 +163,8 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE { void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_jakstat_adjoint(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_jakstat_adjoint(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint_dsigmaydp.cpp b/models/model_jakstat_adjoint/model_jakstat_adjoint_dsigmaydp.cpp index 38045122ba..7313fb731b 100644 --- a/models/model_jakstat_adjoint/model_jakstat_adjoint_dsigmaydp.cpp +++ b/models/model_jakstat_adjoint/model_jakstat_adjoint_dsigmaydp.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_jakstat_adjoint{ -void dsigmaydp_model_jakstat_adjoint(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) { +void dsigmaydp_model_jakstat_adjoint(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) { switch (ip) { case 14: { dsigmaydp[0] = 1.0; diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint_sigmay.cpp b/models/model_jakstat_adjoint/model_jakstat_adjoint_sigmay.cpp index 18172dff23..99dd28dff6 100644 --- a/models/model_jakstat_adjoint/model_jakstat_adjoint_sigmay.cpp +++ b/models/model_jakstat_adjoint/model_jakstat_adjoint_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_jakstat_adjoint{ -void sigmay_model_jakstat_adjoint(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_jakstat_adjoint(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = p[14]; sigmay[1] = p[15]; sigmay[2] = p[16]; diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h index 142c8ae8c6..d8a0167aad 100644 --- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h +++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_o2_h #define _amici_model_jakstat_adjoint_o2_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -18,13 +18,13 @@ extern void JSparse_model_jakstat_adjoint_o2(SUNMatrixContent_Sparse JSparse, co extern void Jy_model_jakstat_adjoint_o2(double *nllh, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); extern void dJydsigma_model_jakstat_adjoint_o2(double *dJydsigma, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); extern void dJydy_model_jakstat_adjoint_o2(double *dJydy, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); -extern void dsigmaydp_model_jakstat_adjoint_o2(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip); +extern void dsigmaydp_model_jakstat_adjoint_o2(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip); extern void dwdp_model_jakstat_adjoint_o2(realtype *dwdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl, const realtype *stcl); extern void dwdx_model_jakstat_adjoint_o2(realtype *dwdx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl); extern void dxdotdp_model_jakstat_adjoint_o2(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp); extern void dydp_model_jakstat_adjoint_o2(double *dydp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp); extern void dydx_model_jakstat_adjoint_o2(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx); -extern void sigmay_model_jakstat_adjoint_o2(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_jakstat_adjoint_o2(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void sx0_model_jakstat_adjoint_o2(realtype *sx0, const realtype t,const realtype *x0, const realtype *p, const realtype *k, const int ip); extern void w_model_jakstat_adjoint_o2(realtype *w, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); extern void x0_model_jakstat_adjoint_o2(realtype *x0, const realtype t, const realtype *p, const realtype *k); @@ -70,7 +70,7 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_jakstat_adjoint_o2(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_jakstat_adjoint_o2(JSparse, t, x, p, k, h, w, dwdx); @@ -124,8 +124,8 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE { void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { - dsigmaydp_model_jakstat_adjoint_o2(dsigmaydp, t, p, k, ip); + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { + dsigmaydp_model_jakstat_adjoint_o2(dsigmaydp, t, p, k, y, ip); } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -163,8 +163,8 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE { void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_jakstat_adjoint_o2(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_jakstat_adjoint_o2(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2_dsigmaydp.cpp b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2_dsigmaydp.cpp index 7ebea670ba..ac9cb8730d 100644 --- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2_dsigmaydp.cpp +++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2_dsigmaydp.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_jakstat_adjoint_o2{ -void dsigmaydp_model_jakstat_adjoint_o2(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) { +void dsigmaydp_model_jakstat_adjoint_o2(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) { switch (ip) { case 14: { dsigmaydp[0] = 1.0; diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2_sigmay.cpp b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2_sigmay.cpp index 26120455ae..6b162c6f03 100644 --- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2_sigmay.cpp +++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_jakstat_adjoint_o2{ -void sigmay_model_jakstat_adjoint_o2(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_jakstat_adjoint_o2(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = p[14]; sigmay[1] = p[15]; sigmay[2] = p[16]; diff --git a/models/model_nested_events/model_nested_events.h b/models/model_nested_events/model_nested_events.h index fcf703b914..bcc1969455 100644 --- a/models/model_nested_events/model_nested_events.h +++ b/models/model_nested_events/model_nested_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_nested_events_h #define _amici_model_nested_events_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -24,7 +24,7 @@ extern void deltax_model_nested_events(double *deltax, const realtype t, const r extern void dxdotdp_model_nested_events(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp); extern void dydx_model_nested_events(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx); extern void root_model_nested_events(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); -extern void sigmay_model_nested_events(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_nested_events(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void stau_model_nested_events(double *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl, const realtype *sx, const int ip, const int ie); extern void sx0_model_nested_events(realtype *sx0, const realtype t,const realtype *x0, const realtype *p, const realtype *k, const int ip); extern void x0_model_nested_events(realtype *x0, const realtype t, const realtype *p, const realtype *k); @@ -70,7 +70,7 @@ class Model_model_nested_events : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_nested_events(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_nested_events(JSparse, t, x, p, k, h, w, dwdx); @@ -127,7 +127,7 @@ class Model_model_nested_events : public amici::Model_ODE { void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -163,8 +163,8 @@ class Model_model_nested_events : public amici::Model_ODE { void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_nested_events(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_nested_events(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_nested_events/model_nested_events_sigmay.cpp b/models/model_nested_events/model_nested_events_sigmay.cpp index 4852d83e1f..0f928ca69a 100644 --- a/models/model_nested_events/model_nested_events_sigmay.cpp +++ b/models/model_nested_events/model_nested_events_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_nested_events{ -void sigmay_model_nested_events(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_nested_events(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = 1.0; } diff --git a/models/model_neuron/model_neuron.h b/models/model_neuron/model_neuron.h index aacdef0da6..d32f24c4fd 100644 --- a/models/model_neuron/model_neuron.h +++ b/models/model_neuron/model_neuron.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_h #define _amici_model_neuron_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -34,7 +34,7 @@ extern void dydx_model_neuron(double *dydx, const realtype t, const realtype *x, extern void dzdx_model_neuron(double *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h); extern void root_model_neuron(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); extern void rz_model_neuron(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h); -extern void sigmay_model_neuron(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_neuron(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void sigmaz_model_neuron(double *sigmaz, const realtype t, const realtype *p, const realtype *k); extern void srz_model_neuron(double *srz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *sx, const int ip); extern void stau_model_neuron(double *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl, const realtype *sx, const int ip, const int ie); @@ -84,7 +84,7 @@ class Model_model_neuron : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_neuron(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_neuron(JSparse, t, x, p, k, h, w, dwdx); @@ -149,7 +149,7 @@ class Model_model_neuron : public amici::Model_ODE { drzdx_model_neuron(drzdx, ie, t, x, p, k, h); } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -187,8 +187,8 @@ class Model_model_neuron : public amici::Model_ODE { rz_model_neuron(rz, ie, t, x, p, k, h); } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_neuron(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_neuron(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_neuron/model_neuron_sigmay.cpp b/models/model_neuron/model_neuron_sigmay.cpp index 092c626156..8bc0d292f8 100644 --- a/models/model_neuron/model_neuron_sigmay.cpp +++ b/models/model_neuron/model_neuron_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_neuron{ -void sigmay_model_neuron(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_neuron(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = 1.0; } diff --git a/models/model_neuron_o2/model_neuron_o2.h b/models/model_neuron_o2/model_neuron_o2.h index 8fb9ce900e..c3879526c3 100644 --- a/models/model_neuron_o2/model_neuron_o2.h +++ b/models/model_neuron_o2/model_neuron_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_o2_h #define _amici_model_neuron_o2_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -35,7 +35,7 @@ extern void dydx_model_neuron_o2(double *dydx, const realtype t, const realtype extern void dzdx_model_neuron_o2(double *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h); extern void root_model_neuron_o2(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); extern void rz_model_neuron_o2(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h); -extern void sigmay_model_neuron_o2(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_neuron_o2(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void sigmaz_model_neuron_o2(double *sigmaz, const realtype t, const realtype *p, const realtype *k); extern void srz_model_neuron_o2(double *srz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *sx, const int ip); extern void stau_model_neuron_o2(double *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl, const realtype *sx, const int ip, const int ie); @@ -86,7 +86,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_neuron_o2(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_neuron_o2(JSparse, t, x, p, k, h, w, dwdx); @@ -151,7 +151,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE { drzdx_model_neuron_o2(drzdx, ie, t, x, p, k, h); } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -190,8 +190,8 @@ class Model_model_neuron_o2 : public amici::Model_ODE { rz_model_neuron_o2(rz, ie, t, x, p, k, h); } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_neuron_o2(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_neuron_o2(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_neuron_o2/model_neuron_o2_sigmay.cpp b/models/model_neuron_o2/model_neuron_o2_sigmay.cpp index 5938f5e08f..3e5743cb70 100644 --- a/models/model_neuron_o2/model_neuron_o2_sigmay.cpp +++ b/models/model_neuron_o2/model_neuron_o2_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_neuron_o2{ -void sigmay_model_neuron_o2(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_neuron_o2(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = 1.0; } diff --git a/models/model_robertson/model_robertson.h b/models/model_robertson/model_robertson.h index 041e793d0a..7e0756ffef 100644 --- a/models/model_robertson/model_robertson.h +++ b/models/model_robertson/model_robertson.h @@ -1,6 +1,6 @@ #ifndef _amici_model_robertson_h #define _amici_model_robertson_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -23,7 +23,7 @@ extern void dwdp_model_robertson(realtype *dwdp, const realtype t, const realtyp extern void dwdx_model_robertson(realtype *dwdx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl); extern void dxdotdp_model_robertson(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *dx, const realtype *w, const realtype *dwdp); extern void dydx_model_robertson(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx); -extern void sigmay_model_robertson(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_robertson(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void w_model_robertson(realtype *w, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); extern void x0_model_robertson(realtype *x0, const realtype t, const realtype *p, const realtype *k); extern void xdot_model_robertson(realtype *xdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *w); @@ -68,7 +68,7 @@ class Model_model_robertson : public amici::Model_DAE { amici::Model* clone() const override { return new Model_model_robertson(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override { JSparse_model_robertson(JSparse, t, x, p, k, h, cj, dx, w, dwdx); @@ -126,7 +126,7 @@ class Model_model_robertson : public amici::Model_DAE { void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -163,8 +163,8 @@ class Model_model_robertson : public amici::Model_DAE { void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_robertson(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_robertson(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_robertson/model_robertson_sigmay.cpp b/models/model_robertson/model_robertson_sigmay.cpp index 6e57e3c50f..9e1e5e019c 100644 --- a/models/model_robertson/model_robertson_sigmay.cpp +++ b/models/model_robertson/model_robertson_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_robertson{ -void sigmay_model_robertson(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_robertson(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = 1.0; sigmay[1] = 1.0; sigmay[2] = 1.0; diff --git a/models/model_steadystate/model_steadystate.h b/models/model_steadystate/model_steadystate.h index f955d97c84..dfbddbd850 100644 --- a/models/model_steadystate/model_steadystate.h +++ b/models/model_steadystate/model_steadystate.h @@ -1,6 +1,6 @@ #ifndef _amici_model_steadystate_h #define _amici_model_steadystate_h -/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */ +/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */ #include #include #include "amici/defines.h" @@ -22,7 +22,7 @@ extern void dwdp_model_steadystate(realtype *dwdp, const realtype t, const realt extern void dwdx_model_steadystate(realtype *dwdx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl); extern void dxdotdp_model_steadystate(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp); extern void dydx_model_steadystate(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx); -extern void sigmay_model_steadystate(double *sigmay, const realtype t, const realtype *p, const realtype *k); +extern void sigmay_model_steadystate(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y); extern void w_model_steadystate(realtype *w, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl); extern void x0_model_steadystate(realtype *x0, const realtype t, const realtype *p, const realtype *k); extern void xdot_model_steadystate(realtype *xdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w); @@ -67,7 +67,7 @@ class Model_model_steadystate : public amici::Model_ODE { amici::Model* clone() const override { return new Model_model_steadystate(*this); }; - std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; }; + std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; }; void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { JSparse_model_steadystate(JSparse, t, x, p, k, h, w, dwdx); @@ -121,7 +121,7 @@ class Model_model_steadystate : public amici::Model_ODE { void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override { } void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { @@ -158,8 +158,8 @@ class Model_model_steadystate : public amici::Model_ODE { void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { } - void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { - sigmay_model_steadystate(sigmay, t, p, k); + void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override { + sigmay_model_steadystate(sigmay, t, p, k, y); } void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { diff --git a/models/model_steadystate/model_steadystate_sigmay.cpp b/models/model_steadystate/model_steadystate_sigmay.cpp index 8e9ebc4417..d895acd6ca 100644 --- a/models/model_steadystate/model_steadystate_sigmay.cpp +++ b/models/model_steadystate/model_steadystate_sigmay.cpp @@ -10,7 +10,7 @@ namespace amici { namespace model_model_steadystate{ -void sigmay_model_steadystate(double *sigmay, const realtype t, const realtype *p, const realtype *k) { +void sigmay_model_steadystate(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) { sigmay[0] = 1.0; sigmay[1] = 1.0; sigmay[2] = 1.0; diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index dd3e9d893a..238df8f629 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -172,12 +172,12 @@ class _FunctionInfo: 'dsigmaydp': _FunctionInfo( 'realtype *dsigmaydp, const realtype t, const realtype *p, ' - 'const realtype *k, const int ip', + 'const realtype *k, const realtype *y, const int ip', ), 'sigmay': _FunctionInfo( 'realtype *sigmay, const realtype t, const realtype *p, ' - 'const realtype *k', + 'const realtype *k, const realtype *y', ), 'sroot': _FunctionInfo( diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index 3385a2cdcb..e66f164127 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -97,6 +97,48 @@ def test_nosensi(simple_sbml_model): assert rdata.status == amici.AMICI_ERROR +def test_sbml2amici_observable_dependent_error(simple_sbml_model): + """Check gradients for model with observable-dependent error""" + sbml_doc, sbml_model = simple_sbml_model + # add parameter and rate rule + relative_sigma = sbml_model.createParameter() + relative_sigma.setId('relative_sigma') + relative_sigma.setValue(0.05) + rr = sbml_model.createRateRule() + rr.setVariable("S1") + rr.setMath(libsbml.parseL3Formula("1")) + sbml_model.getSpecies("S1").setInitialConcentration(1.0) + sbml_importer = SbmlImporter(sbml_source=sbml_model, + from_file=False) + + with TemporaryDirectory() as tmpdir: + sbml_importer.sbml2amici( + model_name="test", + output_dir=tmpdir, + observables={'observed_s1': {'formula': 'S1'}}, + sigmas={'observed_s1': '0.1 + relative_sigma * observed_s1'}, + ) + model_module = amici.import_model_module(module_name='test', + module_path=tmpdir) + model = model_module.getModel() + model.setTimepoints(np.linspace(0, 60, 61)) + solver = model.getSolver() + + # generate artificial data + rdata = amici.runAmiciSimulation(model, solver) + assert np.allclose(rdata.sigmay, 0.1 + 0.05 * rdata.y) + edata = amici.ExpData(rdata, 1.0, 0.0) + edata.setObservedDataStdDev(np.nan) + + solver.setSensitivityOrder(amici.SensitivityOrder.first) + for sensitivity_method in (amici.SensitivityMethod.forward, + amici.SensitivityMethod.adjoint): + solver.setSensitivityMethod(sensitivity_method) + check_derivatives(model, solver, edata) + rdata = amici.runAmiciSimulation(model, solver, edata) + assert rdata.sllh[1] != 0.0 + + @pytest.fixture def model_steadystate_module(): sbml_file = os.path.join(os.path.dirname(__file__), '..', diff --git a/src/abstract_model.cpp b/src/abstract_model.cpp index 486239e103..e207f78a54 100644 --- a/src/abstract_model.cpp +++ b/src/abstract_model.cpp @@ -339,7 +339,8 @@ void AbstractModel::fsigmay(realtype* /*sigmay*/, const realtype /*t*/, const realtype* /*p*/, - const realtype* /*k*/) + const realtype* /*k*/, + const realtype */*y*/) { throw AmiException("Requested functionality is not supported as %s is " "not implemented for this model!", @@ -351,6 +352,7 @@ AbstractModel::fdsigmaydp(realtype* /*dsigmaydp*/, const realtype /*t*/, const realtype* /*p*/, const realtype* /*k*/, + const realtype */*y*/, const int /*ip*/) { throw AmiException("Requested functionality is not supported as %s is " diff --git a/src/model.cpp b/src/model.cpp index b2dc5fc008..bfa20a86c2 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1447,8 +1447,11 @@ void Model::fsigmay(const int it, const ExpData *edata) { derived_state_.sigmay_.assign(ny, 0.0); - fsigmay(derived_state_.sigmay_.data(), getTimepoint(it), state_.unscaledParameters.data(), - state_.fixedParameters.data()); + fsigmay(derived_state_.sigmay_.data(), + getTimepoint(it), + state_.unscaledParameters.data(), + state_.fixedParameters.data(), + derived_state_.y_.data()); if (edata) { auto sigmay_edata = edata->getObservedDataStdDevPtr(it); @@ -1481,6 +1484,7 @@ void Model::fdsigmaydp(const int it, const ExpData *edata) { fdsigmaydp(&derived_state_.dsigmaydp_.at(ip * ny), getTimepoint(it), state_.unscaledParameters.data(), state_.fixedParameters.data(), + derived_state_.y_.data(), plist(ip)); // sigmas in edata override model-sigma -> for those sigmas, set dsigmaydp From 44690d6ad5c1e49c203acc1300ed16475233985f Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 1 Mar 2022 23:39:03 +0100 Subject: [PATCH 02/24] doc --- include/amici/abstract_model.h | 1 + 1 file changed, 1 insertion(+) diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h index 99202ca0f1..164ea5daa4 100644 --- a/include/amici/abstract_model.h +++ b/include/amici/abstract_model.h @@ -572,6 +572,7 @@ class AbstractModel { * @param t current time * @param p parameter vector * @param k constant vector + * @param y model output at timepoint t * @param ip sensitivity index */ virtual void fdsigmaydp(realtype *dsigmaydp, const realtype t, From b46e2bd3a2e7fc3376b608fa796c2044d8b8a80e Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 00:15:56 +0100 Subject: [PATCH 03/24] fix test case --- python/tests/test_sbml_import.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index e66f164127..9ba733c286 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -106,8 +106,10 @@ def test_sbml2amici_observable_dependent_error(simple_sbml_model): relative_sigma.setValue(0.05) rr = sbml_model.createRateRule() rr.setVariable("S1") - rr.setMath(libsbml.parseL3Formula("1")) + rr.setMath(libsbml.parseL3Formula("p1 * S1")) sbml_model.getSpecies("S1").setInitialConcentration(1.0) + sbml_model.getParameter("p1").setValue(0.2) + sbml_importer = SbmlImporter(sbml_source=sbml_model, from_file=False) @@ -130,12 +132,14 @@ def test_sbml2amici_observable_dependent_error(simple_sbml_model): edata = amici.ExpData(rdata, 1.0, 0.0) edata.setObservedDataStdDev(np.nan) + # check sensitivities solver.setSensitivityOrder(amici.SensitivityOrder.first) for sensitivity_method in (amici.SensitivityMethod.forward, amici.SensitivityMethod.adjoint): solver.setSensitivityMethod(sensitivity_method) check_derivatives(model, solver, edata) rdata = amici.runAmiciSimulation(model, solver, edata) + assert rdata.sllh[0] != 0.0 assert rdata.sllh[1] != 0.0 From 50d421191834f318f9e522ef2d47b386fd613350 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 00:16:57 +0100 Subject: [PATCH 04/24] fix test case --- python/tests/test_sbml_import.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index 9ba733c286..a785817d4b 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -101,14 +101,14 @@ def test_sbml2amici_observable_dependent_error(simple_sbml_model): """Check gradients for model with observable-dependent error""" sbml_doc, sbml_model = simple_sbml_model # add parameter and rate rule + sbml_model.getSpecies("S1").setInitialConcentration(1.0) + sbml_model.getParameter("p1").setValue(0.2) + rr = sbml_model.createRateRule() + rr.setVariable("S1") + rr.setMath(libsbml.parseL3Formula("p1")) relative_sigma = sbml_model.createParameter() relative_sigma.setId('relative_sigma') relative_sigma.setValue(0.05) - rr = sbml_model.createRateRule() - rr.setVariable("S1") - rr.setMath(libsbml.parseL3Formula("p1 * S1")) - sbml_model.getSpecies("S1").setInitialConcentration(1.0) - sbml_model.getParameter("p1").setValue(0.2) sbml_importer = SbmlImporter(sbml_source=sbml_model, from_file=False) From 69b8836b524c8e20f50b9c27d16c4b2e051d1938 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 10:53:18 +0100 Subject: [PATCH 05/24] observable symbols with assumptions --- python/amici/sbml_import.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index c4f713d11c..45d8444d0f 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -1214,6 +1214,10 @@ def _process_observables( # add user-provided observables or make all species, and compartments # with assignment rules, observable if observables: + # gather local symbols before parsing observable and sigma formulas + for obs in observables.keys(): + self.add_local_symbol(obs, symbol_with_assumptions(obs)) + self.symbols[SymbolId.OBSERVABLE] = { symbol_with_assumptions(obs): { 'name': definition.get('name', f'y{iobs}'), @@ -1230,13 +1234,9 @@ def _process_observables( for iobs, (obs, definition) in enumerate(observables.items()) } # check for nesting of observables (unsupported) - # performing string-based check, because lhs symbols are `real`, - # whereas rhs symbols are not. - observable_syms = { - str(obs) for obs in self.symbols[SymbolId.OBSERVABLE].keys() - } + observable_syms = set(self.symbols[SymbolId.OBSERVABLE].keys()) for obs in self.symbols[SymbolId.OBSERVABLE].values(): - if any(str(sym) in observable_syms + if any(sym in observable_syms for sym in obs['value'].free_symbols): raise ValueError( "Nested observables are not supported, " From 172bdab223f482cb4f16870f276e1c275e4d53ca Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 12:18:32 +0100 Subject: [PATCH 06/24] dsigmaysy --- include/amici/abstract_model.h | 16 +++++++++++++++- include/amici/model.h | 9 +++++++++ include/amici/model_state.h | 7 ++++++- python/amici/ode_export.py | 5 +++++ src/abstract_model.cpp | 13 +++++++++++++ src/model.cpp | 32 ++++++++++++++++++++++++++++++++ src/model_header.ODE_template.h | 3 +++ 7 files changed, 83 insertions(+), 2 deletions(-) diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h index 164ea5daa4..d333fab1d0 100644 --- a/include/amici/abstract_model.h +++ b/include/amici/abstract_model.h @@ -567,7 +567,7 @@ class AbstractModel { const realtype *k, const realtype *y); /** - * @brief Model-specific implementation of fsigmay + * @brief Model-specific implementation of fdsigmaydp * @param dsigmaydp partial derivative of standard deviation of measurements * @param t current time * @param p parameter vector @@ -578,6 +578,20 @@ class AbstractModel { virtual void fdsigmaydp(realtype *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, int ip); + /** + * @brief Model-specific implementation of fsigmay + * @param dsigmaydy partial derivative of standard deviation of measurements + * w.r.t. model outputs + * @param t current time + * @param p parameter vector + * @param k constant vector + * @param y model output at timepoint t + * @param ip sensitivity index + */ + virtual void fdsigmaydy(realtype *dsigmaydy, const realtype t, + const realtype *p, const realtype *k, + const realtype *y, int ip); + /** * @brief Model-specific implementation of fsigmaz diff --git a/include/amici/model.h b/include/amici/model.h index b9e8fb7b8a..1a06cb3f92 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -115,6 +115,7 @@ class Model : public AbstractModel, public ModelDimensions { using AbstractModel::fdrzdp; using AbstractModel::fdrzdx; using AbstractModel::fdsigmaydp; + using AbstractModel::fdsigmaydy; using AbstractModel::fdsigmazdp; using AbstractModel::fdwdp; using AbstractModel::fdwdp_colptrs; @@ -1396,6 +1397,14 @@ class Model : public AbstractModel, public ModelDimensions { */ void fdsigmaydp(int it, const ExpData *edata); + /** + * @brief Compute partial derivative of standard deviation of measurements + * w.r.t. model outputs. + * @param it Timepoint index + * @param edata pointer to `amici::ExpData` data instance holding sigma values + */ + void fdsigmaydy(int it, const ExpData *edata); + /** * @brief Compute negative log-likelihood of measurements \f$ y \f$. * diff --git a/include/amici/model_state.h b/include/amici/model_state.h index 578835283a..3a851fec3b 100644 --- a/include/amici/model_state.h +++ b/include/amici/model_state.h @@ -228,11 +228,16 @@ struct ModelStateDerived { /** data standard deviation for current timepoint (dimension: ny) */ std::vector sigmay_; - /** temporary storage for parameter derivative of data standard deviation, + /** temporary storage for parameter derivative of data standard deviation, * (dimension: ny x nplist, row-major) */ std::vector dsigmaydp_; + /** temporary storage for observable derivative of data standard deviation, + * (dimension: ny x nplist, row-major) + */ + std::vector dsigmaydy_; + /** temporary storage for event-resolved observable (dimension: nz) */ std::vector z_; diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index 238df8f629..863c9c9be3 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -169,6 +169,11 @@ class _FunctionInfo: 'const int ip, const realtype *w, const realtype *tcl, ' 'const realtype *dtcldp', ), + 'dsigmaydy': + _FunctionInfo( + 'realtype *dsigmaydy, const realtype t, const realtype *p, ' + 'const realtype *k, const realtype *y, const int ip', + ), 'dsigmaydp': _FunctionInfo( 'realtype *dsigmaydp, const realtype t, const realtype *p, ' diff --git a/src/abstract_model.cpp b/src/abstract_model.cpp index e207f78a54..859b9d8424 100644 --- a/src/abstract_model.cpp +++ b/src/abstract_model.cpp @@ -360,6 +360,19 @@ AbstractModel::fdsigmaydp(realtype* /*dsigmaydp*/, __func__); } +void +AbstractModel::fdsigmaydy(realtype */*dsigmaydy*/, + const realtype /*t*/, + const realtype */*p*/, + const realtype */*k*/, + const realtype */*y*/, + int /*ip*/) +{ + throw AmiException("Requested functionality is not supported as %s is " + "not implemented for this model!", + __func__); +} + void AbstractModel::fsigmaz(realtype* /*sigmaz*/, const realtype /*t*/, diff --git a/src/model.cpp b/src/model.cpp index bfa20a86c2..eba2236654 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1504,6 +1504,38 @@ void Model::fdsigmaydp(const int it, const ExpData *edata) { } } +void Model::fdsigmaydy(const int it, const ExpData *edata) { + if (!ny) + return; + + derived_state_.dsigmaydy_.assign(ny * nplist(), 0.0); + + for (int ip = 0; ip < nplist(); ip++) + // get dsigmaydy slice (ny) for current timepoint and parameter + fdsigmaydy(&derived_state_.dsigmaydy_.at(ip * ny), getTimepoint(it), + state_.unscaledParameters.data(), + state_.fixedParameters.data(), + derived_state_.y_.data(), + plist(ip)); + + // sigmas in edata override model-sigma -> for those sigmas, set dsigmaydy + // to zero + if (edata) { + for (int iy = 0; iy < nytrue; iy++) { + if (!edata->isSetObservedDataStdDev(it, iy)) + continue; + for (int ip = 0; ip < nplist(); ip++) { + derived_state_.dsigmaydy_.at(ip * ny + iy) = 0.0; + } + } + } + + if (always_check_finite_) { + app->checkFinite(derived_state_.dsigmaydy_, "dsigmaydy"); + } +} + + void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { if (!ny) return; diff --git a/src/model_header.ODE_template.h b/src/model_header.ODE_template.h index bcf3c8feff..c253d79efc 100644 --- a/src/model_header.ODE_template.h +++ b/src/model_header.ODE_template.h @@ -56,6 +56,7 @@ TPL_DYDX_DEF TPL_DYDP_DEF TPL_SIGMAY_DEF TPL_DSIGMAYDP_DEF +TPL_DSIGMAYDY_DEF TPL_W_DEF TPL_X0_DEF TPL_X0_FIXEDPARAMETERS_DEF @@ -293,6 +294,8 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { TPL_DSIGMAYDP_IMPL + TPL_DSIGMAYDY_IMPL + /** * @brief model specific implementation of fsigmaz * @param dsigmazdp partial derivative of standard deviation of event From 831fc50d1a4da4ab655a5f5f47fca1006c1de0b9 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 14:09:30 +0100 Subject: [PATCH 07/24] ssigmay --- include/amici/model.h | 8 +++++++- src/model.cpp | 29 ++++++++++++++++++++++++++++- src/rdata.cpp | 16 ++++++++++------ 3 files changed, 45 insertions(+), 8 deletions(-) diff --git a/include/amici/model.h b/include/amici/model.h index 1a06cb3f92..e94cbccc07 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -852,12 +852,18 @@ class Model : public AbstractModel, public ModelDimensions { * Total derivative (can be used with both adjoint and forward sensitivity). * * @param ssigmay Buffer (shape `ny` x `nplist`, row-major) + * @param sy Sensitivity of time-resolved observables for current timepoint * @param it Timepoint index * @param edata Pointer to experimental data instance (optional, pass * `nullptr` to ignore) + * @param x State variables + * @param sx State sensitivities */ void getObservableSigmaSensitivity(gsl::span ssigmay, - const int it, const ExpData *edata); + gsl::span sy, + const int it, const ExpData *edata, + const AmiVector &x, + const AmiVectorArray &sx); /** * @brief Add time-resolved measurement negative log-likelihood \f$ Jy \f$. diff --git a/src/model.cpp b/src/model.cpp index eba2236654..c4afe30fb3 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -849,9 +849,36 @@ void Model::getObservableSigma(gsl::span sigmay, const int it, } void Model::getObservableSigmaSensitivity(gsl::span ssigmay, - const int it, const ExpData *edata) { + gsl::span sy, + const int it, const ExpData *edata, + const AmiVector &x, + const AmiVectorArray &sx) { fdsigmaydp(it, edata); writeSlice(derived_state_.dsigmaydp_, ssigmay); + + if(pythonGenerated) { + // = dsigmaydy*(dydx_solver*sx+dydp)+dsigmaydp + // = dsigmaydy*sy+dsigmaydp + + fdydp(getTimepoint(it), x); + fdydx(getTimepoint(it), x); + fdsigmaydy(it, edata); + + derived_state_.sx_.resize(nx_solver * nplist()); + sx.flatten_to_vector(derived_state_.sx_); + + // compute 1.0 * dsigmaydp += 1.0*dsigmaydy*sy + // dsigmaydp C[ny,nplist] += dsigmaydy A[ny,ny] * sy B[ny,nplist] + // M N M K K N + // ldc lda ldb + amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, + BLASTranspose::noTrans, ny, nplist(), ny, 1.0, + derived_state_.dsigmaydy_.data(), ny, + sy.data(), nx_solver, 1.0, ssigmay.data(), ny); + } + + if (always_check_finite_) + checkFinite(ssigmay, "ssigmay"); } void Model::addObservableObjective(realtype &Jy, const int it, diff --git a/src/rdata.cpp b/src/rdata.cpp index 3eab8f5c13..f01393d209 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -303,16 +303,18 @@ void ReturnData::getDataOutput(int it, Model &model, ExpData const *edata) { if (sensi >= SensitivityOrder::first && nplist > 0) { - if (!ssigmay.empty()) - model.getObservableSigmaSensitivity(slice(ssigmay, it, nplist * ny), - it, edata); - if (sensi_meth == SensitivityMethod::forward) { getDataSensisFSA(it, model, edata); } else if (edata && !sllh.empty()) { model.addPartialObservableObjectiveSensitivity( sllh, s2llh, it, x_solver_, *edata); } + + if (!ssigmay.empty()) + model.getObservableSigmaSensitivity( + slice(ssigmay, it, nplist * ny), + sy.empty() ? gsl::span() : slice(sy, it, nplist * ny), + it, edata, x_solver_, sx_solver_); } } @@ -866,7 +868,8 @@ void ReturnData::fsres(const int it, Model &model, const ExpData &edata) { std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); std::vector ssigmay_it(ny * nplist, 0.0); - model.getObservableSigmaSensitivity(ssigmay_it, it, &edata); + model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata, x_solver_, + sx_solver_); auto observedData = edata.getObservedDataPtr(it); for (int iy = 0; iy < nytrue; ++iy) { @@ -904,7 +907,8 @@ void ReturnData::fFIM(int it, Model &model, const ExpData &edata) { std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); std::vector ssigmay_it(ny * nplist, 0.0); - model.getObservableSigmaSensitivity(ssigmay_it, it, &edata); + model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata, + x_solver_, sx_solver_); /* * https://www.wolframalpha.com/input/?i=d%2Fdu+d%2Fdv+log%28s%28u%2Cv%29%29+%2B+0.5+*+%28r%28u%2Cv%29%2Fs%28u%2Cv%29%29%5E2 From b8ddd8c2502c9eee5f9507863763208bd3eb0223 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 14:23:26 +0100 Subject: [PATCH 08/24] .. --- src/rdata.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/rdata.cpp b/src/rdata.cpp index f01393d209..1b182c5d2b 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -142,7 +142,6 @@ void ReturnData::initializeFullReporting(bool quadratic_llh) { srz.resize(nmaxevent * nz * nplist, 0.0); } - ssigmay.resize(nt * ny * nplist, 0.0); ssigmaz.resize(nmaxevent * nz * nplist, 0.0); if (sensi >= SensitivityOrder::second && sensi_meth == SensitivityMethod::forward) @@ -311,10 +310,10 @@ void ReturnData::getDataOutput(int it, Model &model, ExpData const *edata) { } if (!ssigmay.empty()) - model.getObservableSigmaSensitivity( - slice(ssigmay, it, nplist * ny), - sy.empty() ? gsl::span() : slice(sy, it, nplist * ny), - it, edata, x_solver_, sx_solver_); + model.getObservableSigmaSensitivity(slice(ssigmay, it, nplist * ny), + slice(sy, it, nplist * ny), + it, edata, x_solver_, + sx_solver_); } } From 5ea1ba340294d5a443fe46fe0116cdf2e2deba6b Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 17:41:17 +0100 Subject: [PATCH 09/24] doc --- include/amici/model_state.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/amici/model_state.h b/include/amici/model_state.h index 3a851fec3b..60b91f2ea9 100644 --- a/include/amici/model_state.h +++ b/include/amici/model_state.h @@ -234,7 +234,7 @@ struct ModelStateDerived { std::vector dsigmaydp_; /** temporary storage for observable derivative of data standard deviation, - * (dimension: ny x nplist, row-major) + * (dimension: ny x ny, row-major) */ std::vector dsigmaydy_; From b6e3286a0c80051b9c65d8799344acf24550fd8e Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 17:41:28 +0100 Subject: [PATCH 10/24] fdJydy --- src/model.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/model.cpp b/src/model.cpp index c4afe30fb3..bc197becc3 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1587,6 +1587,18 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { derived_state_.sigmay_.data(), edata.getObservedDataPtr(it)); + // add dJydsigma * dsigmaydy + fdJydsigma(it, x, edata); + fdsigmaydy(it, &edata); + + // dJydy_ += dJydsigma * dsigmaydy + // (nJ,ny) (nJ,ny) * (ny,ny) + amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, + BLASTranspose::noTrans, nJ, ny, ny, 1.0, + &derived_state_.dJydsigma_.at(iyt * nJ * ny), ny, + derived_state_.dsigmaydy_.data(), ny, 1.0, + derived_state_.dJydy_.at(iyt).data(), nJ); + if (always_check_finite_) { app->checkFinite( gsl::make_span(derived_state_.dJydy_.at(iyt).get()), From 7703cb63053a0f7f89f6dde8002ed54a877307cc Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 21:51:06 +0100 Subject: [PATCH 11/24] fix dgemm --- src/model.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/model.cpp b/src/model.cpp index bc197becc3..9a4786f5d6 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -874,7 +874,7 @@ void Model::getObservableSigmaSensitivity(gsl::span ssigmay, amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, ny, nplist(), ny, 1.0, derived_state_.dsigmaydy_.data(), ny, - sy.data(), nx_solver, 1.0, ssigmay.data(), ny); + sy.data(), ny, 1.0, ssigmay.data(), ny); } if (always_check_finite_) @@ -1595,7 +1595,7 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { // (nJ,ny) (nJ,ny) * (ny,ny) amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, nJ, ny, ny, 1.0, - &derived_state_.dJydsigma_.at(iyt * nJ * ny), ny, + &derived_state_.dJydsigma_.at(iyt * nJ * ny), nJ, derived_state_.dsigmaydy_.data(), ny, 1.0, derived_state_.dJydy_.at(iyt).data(), nJ); From 9ec612cd0ed3ab8f834c433244ecb3cd6f78ea20 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 22:13:26 +0100 Subject: [PATCH 12/24] matlab ssigmay --- src/returndata_matlab.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/returndata_matlab.cpp b/src/returndata_matlab.cpp index 0a1d79961a..739d23cfaf 100644 --- a/src/returndata_matlab.cpp +++ b/src/returndata_matlab.cpp @@ -82,7 +82,7 @@ mxArray *initMatlabReturnFields(ReturnData const *rdata) { } } - if (rdata->ny > 0) { + if (!rdata->ssigmay.empty()) { writeMatlabField3(matlabSolutionStruct, "ssigmay", rdata->ssigmay, rdata->nt, rdata->nplist, rdata->ny, perm2); } if ((rdata->nz > 0) & (rdata->ne > 0)) { From 6512245a936fc6859c0f7bbdda237c951e322e10 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 22:37:01 +0100 Subject: [PATCH 13/24] cleanup --- include/amici/model.h | 4 +--- src/model.cpp | 12 ++++-------- src/rdata.cpp | 9 ++++----- 3 files changed, 9 insertions(+), 16 deletions(-) diff --git a/include/amici/model.h b/include/amici/model.h index e94cbccc07..4c84c374e5 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -857,13 +857,11 @@ class Model : public AbstractModel, public ModelDimensions { * @param edata Pointer to experimental data instance (optional, pass * `nullptr` to ignore) * @param x State variables - * @param sx State sensitivities */ void getObservableSigmaSensitivity(gsl::span ssigmay, gsl::span sy, const int it, const ExpData *edata, - const AmiVector &x, - const AmiVectorArray &sx); + const AmiVector &x); /** * @brief Add time-resolved measurement negative log-likelihood \f$ Jy \f$. diff --git a/src/model.cpp b/src/model.cpp index 9a4786f5d6..2d28a67d7e 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -851,23 +851,19 @@ void Model::getObservableSigma(gsl::span sigmay, const int it, void Model::getObservableSigmaSensitivity(gsl::span ssigmay, gsl::span sy, const int it, const ExpData *edata, - const AmiVector &x, - const AmiVectorArray &sx) { + const AmiVector &x) { fdsigmaydp(it, edata); writeSlice(derived_state_.dsigmaydp_, ssigmay); if(pythonGenerated) { - // = dsigmaydy*(dydx_solver*sx+dydp)+dsigmaydp - // = dsigmaydy*sy+dsigmaydp + // ssigmay = dsigmaydy*(dydx_solver*sx+dydp)+dsigmaydp + // = dsigmaydy*sy+dsigmaydp fdydp(getTimepoint(it), x); fdydx(getTimepoint(it), x); fdsigmaydy(it, edata); - derived_state_.sx_.resize(nx_solver * nplist()); - sx.flatten_to_vector(derived_state_.sx_); - - // compute 1.0 * dsigmaydp += 1.0*dsigmaydy*sy + // compute ssigmay = 1.0 * dsigmaydp + 1.0 * dsigmaydy * sy // dsigmaydp C[ny,nplist] += dsigmaydy A[ny,ny] * sy B[ny,nplist] // M N M K K N // ldc lda ldb diff --git a/src/rdata.cpp b/src/rdata.cpp index 1b182c5d2b..a5f1ce689e 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -312,8 +312,7 @@ void ReturnData::getDataOutput(int it, Model &model, ExpData const *edata) { if (!ssigmay.empty()) model.getObservableSigmaSensitivity(slice(ssigmay, it, nplist * ny), slice(sy, it, nplist * ny), - it, edata, x_solver_, - sx_solver_); + it, edata, x_solver_); } } @@ -867,8 +866,8 @@ void ReturnData::fsres(const int it, Model &model, const ExpData &edata) { std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); std::vector ssigmay_it(ny * nplist, 0.0); - model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata, x_solver_, - sx_solver_); + model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata, + x_solver_); auto observedData = edata.getObservedDataPtr(it); for (int iy = 0; iy < nytrue; ++iy) { @@ -907,7 +906,7 @@ void ReturnData::fFIM(int it, Model &model, const ExpData &edata) { model.getObservableSigma(sigmay_it, it, &edata); std::vector ssigmay_it(ny * nplist, 0.0); model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata, - x_solver_, sx_solver_); + x_solver_); /* * https://www.wolframalpha.com/input/?i=d%2Fdu+d%2Fdv+log%28s%28u%2Cv%29%29+%2B+0.5+*+%28r%28u%2Cv%29%2Fs%28u%2Cv%29%29%5E2 From ed2a979d09e7688892784f980113a08017958729 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 2 Mar 2022 23:10:18 +0100 Subject: [PATCH 14/24] cleanup --- include/amici/model.h | 3 +-- src/model.cpp | 6 +----- src/rdata.cpp | 8 +++----- 3 files changed, 5 insertions(+), 12 deletions(-) diff --git a/include/amici/model.h b/include/amici/model.h index 4c84c374e5..cc053b8fee 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -860,8 +860,7 @@ class Model : public AbstractModel, public ModelDimensions { */ void getObservableSigmaSensitivity(gsl::span ssigmay, gsl::span sy, - const int it, const ExpData *edata, - const AmiVector &x); + const int it, const ExpData *edata); /** * @brief Add time-resolved measurement negative log-likelihood \f$ Jy \f$. diff --git a/src/model.cpp b/src/model.cpp index 2d28a67d7e..a12e10b2c6 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -850,8 +850,7 @@ void Model::getObservableSigma(gsl::span sigmay, const int it, void Model::getObservableSigmaSensitivity(gsl::span ssigmay, gsl::span sy, - const int it, const ExpData *edata, - const AmiVector &x) { + const int it, const ExpData *edata) { fdsigmaydp(it, edata); writeSlice(derived_state_.dsigmaydp_, ssigmay); @@ -859,8 +858,6 @@ void Model::getObservableSigmaSensitivity(gsl::span ssigmay, // ssigmay = dsigmaydy*(dydx_solver*sx+dydp)+dsigmaydp // = dsigmaydy*sy+dsigmaydp - fdydp(getTimepoint(it), x); - fdydx(getTimepoint(it), x); fdsigmaydy(it, edata); // compute ssigmay = 1.0 * dsigmaydp + 1.0 * dsigmaydy * sy @@ -1648,7 +1645,6 @@ void Model::fdJydp(const int it, const AmiVector &x, const ExpData &edata) { // dJydy nJ, nytrue x ny // dydp nplist * ny // dJydp nplist x nJ - // dJydsigma if (!ny) return; diff --git a/src/rdata.cpp b/src/rdata.cpp index a5f1ce689e..00b08e7353 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -312,7 +312,7 @@ void ReturnData::getDataOutput(int it, Model &model, ExpData const *edata) { if (!ssigmay.empty()) model.getObservableSigmaSensitivity(slice(ssigmay, it, nplist * ny), slice(sy, it, nplist * ny), - it, edata, x_solver_); + it, edata); } } @@ -866,8 +866,7 @@ void ReturnData::fsres(const int it, Model &model, const ExpData &edata) { std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); std::vector ssigmay_it(ny * nplist, 0.0); - model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata, - x_solver_); + model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata); auto observedData = edata.getObservedDataPtr(it); for (int iy = 0; iy < nytrue; ++iy) { @@ -905,8 +904,7 @@ void ReturnData::fFIM(int it, Model &model, const ExpData &edata) { std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); std::vector ssigmay_it(ny * nplist, 0.0); - model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata, - x_solver_); + model.getObservableSigmaSensitivity(ssigmay_it, sy_it, it, &edata); /* * https://www.wolframalpha.com/input/?i=d%2Fdu+d%2Fdv+log%28s%28u%2Cv%29%29+%2B+0.5+*+%28r%28u%2Cv%29%2Fs%28u%2Cv%29%29%5E2 From 30c0a98a3621b435456ed48081d925a022e615b6 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 3 Mar 2022 14:55:01 +0100 Subject: [PATCH 15/24] ny=2 --- python/tests/test_sbml_import.py | 55 +++++++++++++++++++------------- 1 file changed, 32 insertions(+), 23 deletions(-) diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index a785817d4b..4637e149ad 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -97,8 +97,8 @@ def test_nosensi(simple_sbml_model): assert rdata.status == amici.AMICI_ERROR -def test_sbml2amici_observable_dependent_error(simple_sbml_model): - """Check gradients for model with observable-dependent error""" +@pytest.fixture +def observable_dependent_error_model(simple_sbml_model): sbml_doc, sbml_model = simple_sbml_model # add parameter and rate rule sbml_model.getSpecies("S1").setInitialConcentration(1.0) @@ -117,30 +117,39 @@ def test_sbml2amici_observable_dependent_error(simple_sbml_model): sbml_importer.sbml2amici( model_name="test", output_dir=tmpdir, - observables={'observed_s1': {'formula': 'S1'}}, - sigmas={'observed_s1': '0.1 + relative_sigma * observed_s1'}, + observables={'observable_s1': {'formula': 'S1'}, + 'observable_s1_scaled': {'formula': '0.5 * S1'}}, + sigmas={'observable_s1': '0.1 + relative_sigma * observable_s1', + 'observable_s1_scaled': '0.02 * observable_s1_scaled'}, ) - model_module = amici.import_model_module(module_name='test', - module_path=tmpdir) - model = model_module.getModel() - model.setTimepoints(np.linspace(0, 60, 61)) - solver = model.getSolver() + yield amici.import_model_module(module_name='test', + module_path=tmpdir) - # generate artificial data - rdata = amici.runAmiciSimulation(model, solver) - assert np.allclose(rdata.sigmay, 0.1 + 0.05 * rdata.y) - edata = amici.ExpData(rdata, 1.0, 0.0) - edata.setObservedDataStdDev(np.nan) - # check sensitivities - solver.setSensitivityOrder(amici.SensitivityOrder.first) - for sensitivity_method in (amici.SensitivityMethod.forward, - amici.SensitivityMethod.adjoint): - solver.setSensitivityMethod(sensitivity_method) - check_derivatives(model, solver, edata) - rdata = amici.runAmiciSimulation(model, solver, edata) - assert rdata.sllh[0] != 0.0 - assert rdata.sllh[1] != 0.0 +def test_sbml2amici_observable_dependent_error(observable_dependent_error_model): + """Check gradients for model with observable-dependent error""" + model_module = observable_dependent_error_model + model = model_module.getModel() + model.setTimepoints(np.linspace(0, 60, 61)) + solver = model.getSolver() + + # generate artificial data + rdata = amici.runAmiciSimulation(model, solver) + assert np.allclose(rdata.sigmay[:, 0], 0.1 + 0.05 * rdata.y[:, 0]) + assert np.allclose(rdata.sigmay[:, 1], 0.02 * rdata.y[:, 1]) + edata = amici.ExpData(rdata, 1.0, 0.0) + edata.setObservedDataStdDev(np.nan) + + # check sensitivities + solver.setSensitivityOrder(amici.SensitivityOrder.first) + # FSA + solver.setSensitivityMethod(amici.SensitivityMethod.forward) + rdata = amici.runAmiciSimulation(model, solver, edata) + assert np.any(rdata.ssigmay != 0.0) + check_derivatives(model, solver, edata) + # ASA + solver.setSensitivityMethod(amici.SensitivityMethod.adjoint) + check_derivatives(model, solver, edata) @pytest.fixture From abd0ea33d1ebaee56f4dc4cd0a5994e26e56c7fa Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 3 Mar 2022 14:55:18 +0100 Subject: [PATCH 16/24] cleanup --- src/model.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/model.cpp b/src/model.cpp index a12e10b2c6..8a110c9fcb 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1564,6 +1564,9 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { fsigmay(it, &edata); if (pythonGenerated) { + fdJydsigma(it, x, edata); + fdsigmaydy(it, &edata); + for (int iyt = 0; iyt < nytrue; iyt++) { if (!derived_state_.dJydy_.at(iyt).capacity()) continue; @@ -1580,12 +1583,8 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { derived_state_.sigmay_.data(), edata.getObservedDataPtr(it)); - // add dJydsigma * dsigmaydy - fdJydsigma(it, x, edata); - fdsigmaydy(it, &edata); - // dJydy_ += dJydsigma * dsigmaydy - // (nJ,ny) (nJ,ny) * (ny,ny) + // C(nJ,ny) A(nJ,ny) * B(ny,ny) amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, nJ, ny, ny, 1.0, &derived_state_.dJydsigma_.at(iyt * nJ * ny), nJ, From 6339adcdf9076e2c2861fd399df53451be10b5c9 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 3 Mar 2022 15:08:10 +0100 Subject: [PATCH 17/24] doc --- include/amici/model.h | 1 - 1 file changed, 1 deletion(-) diff --git a/include/amici/model.h b/include/amici/model.h index cc053b8fee..18e5928be3 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -856,7 +856,6 @@ class Model : public AbstractModel, public ModelDimensions { * @param it Timepoint index * @param edata Pointer to experimental data instance (optional, pass * `nullptr` to ignore) - * @param x State variables */ void getObservableSigmaSensitivity(gsl::span ssigmay, gsl::span sy, From c12de8a65bf656ccbbc5a788bbac3284e3fec706 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 3 Mar 2022 21:55:25 +0100 Subject: [PATCH 18/24] Fix mul --- src/model.cpp | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/model.cpp b/src/model.cpp index 8a110c9fcb..d58b7f0165 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1566,6 +1566,7 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { if (pythonGenerated) { fdJydsigma(it, x, edata); fdsigmaydy(it, &edata); + SUNMatrixWrapper tmp_dense(nJ, ny); for (int iyt = 0; iyt < nytrue; iyt++) { if (!derived_state_.dJydy_.at(iyt).capacity()) @@ -1583,13 +1584,27 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { derived_state_.sigmay_.data(), edata.getObservedDataPtr(it)); - // dJydy_ += dJydsigma * dsigmaydy + // dJydy += dJydsigma * dsigmaydy // C(nJ,ny) A(nJ,ny) * B(ny,ny) + // sparse dense dense + tmp_dense.zero(); amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, nJ, ny, ny, 1.0, &derived_state_.dJydsigma_.at(iyt * nJ * ny), nJ, derived_state_.dsigmaydy_.data(), ny, 1.0, - derived_state_.dJydy_.at(iyt).data(), nJ); + tmp_dense.data(), nJ); + + auto tmp_sparse = SUNSparseFromDenseMatrix( + tmp_dense.get(), 0.0, CSC_MAT); + auto ret = SUNMatScaleAdd(1.0, derived_state_.dJydy_.at(iyt).get(), + tmp_sparse); + if(ret != SUNMAT_SUCCESS) { + SUNMatDestroy(tmp_sparse); + throw AmiException("SUNMatScaleAdd failed with status %d in %s", + ret, __func__); + } + SUNMatDestroy(tmp_sparse); + derived_state_.dJydy_.at(iyt).refresh(); if (always_check_finite_) { app->checkFinite( From 353e7f9106c6592921f191a0bc816280afdff3c0 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 3 Mar 2022 22:22:04 +0100 Subject: [PATCH 19/24] update tests/cpp/expectedResults.h5 --- tests/cpp/expectedResults.h5 | Bin 4198368 -> 4198368 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/tests/cpp/expectedResults.h5 b/tests/cpp/expectedResults.h5 index 3525c1b4b1f29ba0c7eef5de599bb5e865918c73..a24be851566e794a6dffd0f0564bed2660e06996 100644 GIT binary patch delta 377 zcmXZVIZnes6o6qRU>k!mge74O#y}Ew%)W(mT!4ZO$O$MCB#K0GVaaVGiYeo{0YXUQ zGPj_0gSY__9lsI9r;*-2PcP~D*#v%&39`XJF!-+BH8PsGZ)j3m$nCEj%SyJ?WY*XZ7E4?22v4v6C7{fSnn7|~aFpWG4n87UOFpmW+qKG9p yScZ!gtYQr%l#$+C#|AdBg>CF$7kluqj{_Xy2*)_VDb8??3tXa#E9J+vhudGk{-PxS delta 424 zcmXZWJ5R!J6vp9D?Y|T%)*E;OD=Hu=g5q6d^9wk*DSQG3FkyhCCMCEM6FWjT9W^Fp zbz^pM?PB~E#vjDt$r;Xj-W;!MN*+whXGzwWHU2K&RDHs}t;&wAe+vJY-jX}7(WjbR zcY8JY@{Q$BLO6EmH=s&c!Or&bLDlXgMdi_pt8OkOzOP6{8l*-QgVrdfXB5@`Xo!k! zgbHfA+!SuPsocfMesz`q?8A#PtTdCFYTx&iVWq!qI0}YSMF2t9aI{#M)1s^K%yU|m z2HNNuSF5+pzJ>{eFbN$~n8pldF^4cBn8yMZv4mwr5yJ`$tRjvztYZTSB$2`<(%8Z_ hc96j?vdCc%`^e(}hd9D9PH>7doTGpX*X$Qd*MH^Ht?d8+ From 24d59da25a0a42840283d2a3ebc5de002c2bf6f2 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 4 Mar 2022 10:34:38 +0100 Subject: [PATCH 20/24] address comments --- src/model.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/model.cpp b/src/model.cpp index d58b7f0165..f20497baa6 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1598,12 +1598,11 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { tmp_dense.get(), 0.0, CSC_MAT); auto ret = SUNMatScaleAdd(1.0, derived_state_.dJydy_.at(iyt).get(), tmp_sparse); + SUNMatDestroy(tmp_sparse); if(ret != SUNMAT_SUCCESS) { - SUNMatDestroy(tmp_sparse); throw AmiException("SUNMatScaleAdd failed with status %d in %s", ret, __func__); } - SUNMatDestroy(tmp_sparse); derived_state_.dJydy_.at(iyt).refresh(); if (always_check_finite_) { From 969d10c65c4a26472c7f1d8775812b554c586d46 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 4 Mar 2022 12:02:54 +0100 Subject: [PATCH 21/24] Duplicate ./python/examples/example_steadystate/model_steadystate_scaled.xml history in ./python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml history. --- ...steadystate_scaled_without_observables.xml | 199 ++++++++++++++++++ 1 file changed, 199 insertions(+) create mode 100644 python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml diff --git a/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml b/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml new file mode 100644 index 0000000000..dc5ae669f7 --- /dev/null +++ b/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml @@ -0,0 +1,199 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + x1 + + + + + x2 + + + + + x3 + + + + + + + scaling_x1 + x1 + + + + + + + + offset_x2 + x2 + + + + + + x1 + + + + + + + + + + + + + + + + p1 + + + x1 + 2 + + + + + + + + + + + + + + + + + + p2 + x1 + x2 + + + + + + + + + + + + + + + + p3 + x2 + + + + + + + + + + + + + + + + + p4 + x3 + + + + + + + + + + + + + k0 + x3 + + + + + + + + + + + p5 + + + + + + + From e179568eef86d7d85a5c285c10c374a140fcd244 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 4 Mar 2022 12:08:26 +0100 Subject: [PATCH 22/24] Fix redundant observables in examples --- documentation/INSTALL.md | 1 - ...steadystate_scaled_without_observables.xml | 1 + .../ExampleSteadystate.ipynb | 6 +-- ...steadystate_scaled_without_observables.xml | 47 ------------------- 4 files changed, 4 insertions(+), 51 deletions(-) delete mode 120000 documentation/INSTALL.md create mode 120000 documentation/model_steadystate_scaled_without_observables.xml diff --git a/documentation/INSTALL.md b/documentation/INSTALL.md deleted file mode 120000 index 71db8b4934..0000000000 --- a/documentation/INSTALL.md +++ /dev/null @@ -1 +0,0 @@ -../INSTALL.md \ No newline at end of file diff --git a/documentation/model_steadystate_scaled_without_observables.xml b/documentation/model_steadystate_scaled_without_observables.xml new file mode 120000 index 0000000000..be8991d93e --- /dev/null +++ b/documentation/model_steadystate_scaled_without_observables.xml @@ -0,0 +1 @@ +../python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml \ No newline at end of file diff --git a/python/examples/example_steadystate/ExampleSteadystate.ipynb b/python/examples/example_steadystate/ExampleSteadystate.ipynb index d1acfe687d..52f014dba8 100644 --- a/python/examples/example_steadystate/ExampleSteadystate.ipynb +++ b/python/examples/example_steadystate/ExampleSteadystate.ipynb @@ -16,7 +16,7 @@ "outputs": [], "source": [ "# SBML model we want to import\n", - "sbml_file = 'model_steadystate_scaled.xml'\n", + "sbml_file = 'model_steadystate_scaled_without_observables.xml'\n", "# Name of the model that will also be the name of the python module\n", "model_name = 'model_steadystate_scaled'\n", "# Directory to which the generated model code is written\n", @@ -133,7 +133,7 @@ "\n", "Specifying observables is beyond the scope of SBML. Here we define them manually.\n", "\n", - "If you are looking for a more scalable way for defining observables, then checkout [PEtab](https://github.com/PEtab-dev/PEtab). Another possibility is using SBML's [`AssignmentRule`s](http://sbml.org/Software/libSBML/5.13.0/docs//python-api/classlibsbml_1_1_rule.html) to specify model outputs within the SBML file." + "If you are looking for a more scalable way for defining observables, then checkout [PEtab](https://github.com/PEtab-dev/PEtab). Another possibility is using SBML's [`AssignmentRule`s](https://sbml.org/software/libsbml/5.18.0/docs/formatted/python-api/classlibsbml_1_1_assignment_rule.html) to specify model outputs within the SBML file." ] }, { @@ -1937,7 +1937,7 @@ ], "source": [ "# look at the States in rdata as DataFrame\n", - "amici.getSimulationStatesAsDataFrame(model, [edata], [rdata])" + "amici.getSimulationStatesAsDataFrame(model, [edata], [rdata])\n" ] } ], diff --git a/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml b/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml index dc5ae669f7..ac6c806770 100644 --- a/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml +++ b/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml @@ -43,56 +43,9 @@ - - - - - - - - - - - x1 - - - - - x2 - - - - - x3 - - - - - - - scaling_x1 - x1 - - - - - - - - offset_x2 - x2 - - - - - - x1 - - - From 6ddd41e28128bdf18e2ebc89c7385d51ab62c8ed Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 4 Mar 2022 12:51:49 +0100 Subject: [PATCH 23/24] RAII wrapper --- src/model.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/model.cpp b/src/model.cpp index f20497baa6..2bb2dd6fbb 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1594,11 +1594,9 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { derived_state_.dsigmaydy_.data(), ny, 1.0, tmp_dense.data(), nJ); - auto tmp_sparse = SUNSparseFromDenseMatrix( - tmp_dense.get(), 0.0, CSC_MAT); + auto tmp_sparse = SUNMatrixWrapper(tmp_dense, 0.0, CSC_MAT); auto ret = SUNMatScaleAdd(1.0, derived_state_.dJydy_.at(iyt).get(), - tmp_sparse); - SUNMatDestroy(tmp_sparse); + tmp_sparse.get()); if(ret != SUNMAT_SUCCESS) { throw AmiException("SUNMatScaleAdd failed with status %d in %s", ret, __func__); From 6be733269ba88dc80d76c50cd4d3bd2367879dde Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 4 Mar 2022 13:03:15 +0100 Subject: [PATCH 24/24] fixup --- .../model_steadystate_scaled_without_observables.xml | 1 + 1 file changed, 1 insertion(+) diff --git a/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml b/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml index ac6c806770..62c673e513 100644 --- a/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml +++ b/python/examples/example_steadystate/model_steadystate_scaled_without_observables.xml @@ -45,6 +45,7 @@ +