From a884dd255944124180f8168f21502c90f793d757 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 18 Jul 2025 16:42:05 +0200 Subject: [PATCH] Handle state-independent, parameter-dependent event triggers without root-finding Since #2227, events that trigger at an explicit timepoint are handled without root-finding. Events with state-independent but parameter-dependent triggers were still handled with root-finding. With these changes, those event triggers will also be handled without root-finding. Closes #2185. --- include/amici/abstract_model.h | 11 +++++ include/amici/model.h | 28 +++++++++-- include/amici/model_dae.h | 7 +-- include/amici/model_ode.h | 7 +-- include/amici/serialization.h | 2 +- include/amici/sundials_matrix_wrapper.h | 1 - models/model_calvetti_py/CMakeLists.txt | 1 + models/model_calvetti_py/explicit_roots.cpp | 24 +++++++++ models/model_calvetti_py/model_calvetti_py.h | 12 +++-- models/model_dirac_py/CMakeLists.txt | 1 + models/model_dirac_py/explicit_roots.cpp | 21 ++++++++ models/model_dirac_py/model_dirac_py.h | 14 ++++-- models/model_events_py/CMakeLists.txt | 1 + models/model_events_py/explicit_roots.cpp | 25 ++++++++++ models/model_events_py/model_events_py.h | 14 ++++-- .../model_jakstat_adjoint_py.h | 8 +-- models/model_nested_events_py/CMakeLists.txt | 1 + models/model_nested_events_py/deltaqB.cpp | 2 +- models/model_nested_events_py/deltasx.cpp | 10 ++-- models/model_nested_events_py/deltax.cpp | 2 +- models/model_nested_events_py/deltaxB.cpp | 4 +- .../model_nested_events_py/explicit_roots.cpp | 21 ++++++++ models/model_nested_events_py/h.h | 6 +-- .../model_nested_events_py.h | 18 ++++--- models/model_nested_events_py/root.cpp | 6 +-- models/model_nested_events_py/stau.cpp | 14 +++--- models/model_neuron_py/model_neuron_py.h | 8 +-- .../model_robertson_py/model_robertson_py.h | 8 +-- .../model_steadystate_py.h | 8 +-- python/sdist/amici/_codegen/cxx_functions.py | 21 +++++++- python/sdist/amici/_codegen/model_class.py | 25 ++-------- python/sdist/amici/de_export.py | 49 ++++++++++++++----- python/sdist/amici/de_model.py | 13 +++-- python/sdist/amici/de_model_components.py | 17 ++++++- python/tests/test_events.py | 37 +++++++++++--- src/forwardproblem.cpp | 13 ++--- src/model.cpp | 40 ++++++++++++--- src/model_header.template.h | 7 +-- tests/sbml/testSBMLSuite.py | 11 +++-- 39 files changed, 376 insertions(+), 142 deletions(-) create mode 100644 models/model_calvetti_py/explicit_roots.cpp create mode 100644 models/model_dirac_py/explicit_roots.cpp create mode 100644 models/model_events_py/explicit_roots.cpp create mode 100644 models/model_nested_events_py/explicit_roots.cpp diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h index 2ebed26d73..b24b9a2e43 100644 --- a/include/amici/abstract_model.h +++ b/include/amici/abstract_model.h @@ -1170,6 +1170,17 @@ class AbstractModel { virtual void fdspline_slopesdp( realtype* dspline_slopesdp, realtype const* p, realtype const* k, int ip ); + + /** + * @brief Compute explicit roots of the model. + * @param p parameter vector + * @param k constant vector + * @return A vector of length ne_solver, each containing a vector of + * explicit roots for the corresponding event. + */ + virtual std::vector> fexplicit_roots( + [[maybe_unused]] realtype const* p, [[maybe_unused]] realtype const* k + ) = 0; }; } // namespace amici diff --git a/include/amici/model.h b/include/amici/model.h index 91bb288a6b..d60222442a 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -112,15 +112,12 @@ class Model : public AbstractModel, public ModelDimensions { * @param idlist Indexes indicating algebraic components (DAE only) * @param z2event Mapping of event outputs to events * @param events Vector of events - * @param state_independent_events Map of events with state-independent - * triggers functions, mapping trigger timepoints to event indices. */ Model( ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode o2mode, std::vector idlist, std::vector z2event, - std::vector events = {}, - std::map> state_independent_events = {} + std::vector events = {} ); /** Destructor. */ @@ -301,6 +298,13 @@ class Model : public AbstractModel, public ModelDimensions { std::vector& roots_found ); + /** + * @brief Re-compute the explicit roots. + * + * Re-compute the explicit roots based on the current model parameters. + */ + void reinit_explicit_roots(); + /** * @brief Get number of parameters wrt to which sensitivities are computed. * @return Length of sensitivity index vector @@ -1514,6 +1518,8 @@ class Model : public AbstractModel, public ModelDimensions { /** * @brief Get trigger times for events that don't require root-finding. * + * To be called only after Model::initialize. + * * @return List of unique trigger points for events that don't require * root-finding (i.e. that trigger at predetermined timepoints), * in ascending order. @@ -1562,6 +1568,18 @@ class Model : public AbstractModel, public ModelDimensions { return any_state_non_negative_; } + [[nodiscard]] std::vector> fexplicit_roots( + [[maybe_unused]] realtype const* p, [[maybe_unused]] realtype const* k + ) override { + if (ne != ne_solver) { + throw AmiException( + "ne!=ne_solver, but 'fexplicit_roots' is not implemented for " + "this model." + ); + } + return {}; + } + /** * Flag indicating whether for * `amici::Solver::sensi_` == `amici::SensitivityOrder::second` @@ -1579,7 +1597,7 @@ class Model : public AbstractModel, public ModelDimensions { * @brief Map of trigger timepoints to event indices for events that don't * require root-finding. */ - std::map> state_independent_events_ = {}; + std::map> explicit_roots_ = {}; protected: /** diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h index 0c83f85e6b..4d0604ebe4 100644 --- a/include/amici/model_dae.h +++ b/include/amici/model_dae.h @@ -34,19 +34,16 @@ class Model_DAE : public Model { * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events * @param events Vector of events - * @param state_independent_events Map of events with state-independent - * triggers functions, mapping trigger timepoints to event indices. */ Model_DAE( ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode const o2mode, std::vector const& idlist, - std::vector const& z2event, std::vector events = {}, - std::map> state_independent_events = {} + std::vector const& z2event, std::vector events = {} ) : Model( model_dimensions, simulation_parameters, o2mode, idlist, z2event, - events, state_independent_events + events ) { SUNContext sunctx = derived_state_.sunctx_; derived_state_.M_ = SUNMatrixWrapper(nx_solver, nx_solver, sunctx); diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h index 190e1e59eb..201d9ce5b0 100644 --- a/include/amici/model_ode.h +++ b/include/amici/model_ode.h @@ -32,19 +32,16 @@ class Model_ODE : public Model { * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events * @param events Vector of events - * @param state_independent_events Map of events with state-independent - * triggers functions, mapping trigger timepoints to event indices. */ Model_ODE( ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode const o2mode, std::vector const& idlist, - std::vector const& z2event, std::vector events = {}, - std::map> state_independent_events = {} + std::vector const& z2event, std::vector events = {} ) : Model( model_dimensions, simulation_parameters, o2mode, idlist, z2event, - events, state_independent_events + events ) {} void diff --git a/include/amici/serialization.h b/include/amici/serialization.h index 4aa00c975d..5d4466650c 100644 --- a/include/amici/serialization.h +++ b/include/amici/serialization.h @@ -154,7 +154,7 @@ void serialize(Archive& ar, amici::Model& m, unsigned int const /*version*/) { ar & m.sigma_res_; ar & m.steadystate_computation_mode_; ar & m.steadystate_sensitivity_mode_; - ar & m.state_independent_events_; + ar & m.explicit_roots_; ar & m.steadystate_mask_; } diff --git a/include/amici/sundials_matrix_wrapper.h b/include/amici/sundials_matrix_wrapper.h index 26dbb9adee..0686b12eb9 100644 --- a/include/amici/sundials_matrix_wrapper.h +++ b/include/amici/sundials_matrix_wrapper.h @@ -591,7 +591,6 @@ class SUNMatrixWrapper { bool ownmat = true; }; - /** * @brief Output formatter for SUNMatrixWrapper. * @param os output stream diff --git a/models/model_calvetti_py/CMakeLists.txt b/models/model_calvetti_py/CMakeLists.txt index f7a3f39f4f..0f97d655a2 100644 --- a/models/model_calvetti_py/CMakeLists.txt +++ b/models/model_calvetti_py/CMakeLists.txt @@ -72,6 +72,7 @@ dxdotdw.h dxdotdx_explicit.cpp dxdotdx_explicit.h dydx.cpp +explicit_roots.cpp h.h k.h model_calvetti_py.cpp diff --git a/models/model_calvetti_py/explicit_roots.cpp b/models/model_calvetti_py/explicit_roots.cpp new file mode 100644 index 0000000000..fb85bacadc --- /dev/null +++ b/models/model_calvetti_py/explicit_roots.cpp @@ -0,0 +1,24 @@ +#include "amici/symbolic_functions.h" +#include "amici/defines.h" +#include "sundials/sundials_types.h" + +#include +#include + +#include +#include "k.h" + +namespace amici { +namespace model_model_calvetti_py { + +std::vector> explicit_roots_model_calvetti_py(const realtype *p, const realtype *k){ + return { + {10}, + {10}, + {12}, + {12} + }; +} + +} // namespace model_model_calvetti_py +} // namespace amici diff --git a/models/model_calvetti_py/model_calvetti_py.h b/models/model_calvetti_py/model_calvetti_py.h index c1b4d67a23..2c9876f0ff 100644 --- a/models/model_calvetti_py/model_calvetti_py.h +++ b/models/model_calvetti_py/model_calvetti_py.h @@ -99,7 +99,7 @@ extern void x_solver_model_calvetti_py(realtype *x_solver, const realtype *x_rda extern std::vector create_splines_model_calvetti_py(const realtype *p, const realtype *k); - +extern std::vector> explicit_roots_model_calvetti_py(const realtype *p, const realtype *k); /** * @brief AMICI-generated model subclass. */ @@ -155,8 +155,7 @@ class Model_model_calvetti_py : public amici::Model_DAE { Event("Heaviside_1", true, true, NAN), Event("Heaviside_2", true, true, NAN), Event("Heaviside_3", true, true, NAN) - }, // events - {{10.0, {0, 1}}, {12.0, {2, 3}}} // state-independent events + } // events ) { } @@ -413,6 +412,11 @@ class Model_model_calvetti_py : public amici::Model_DAE { void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {} + std::vector> fexplicit_roots(const realtype *p, const realtype *k) override { + return explicit_roots_model_calvetti_py(p, k); + } + + std::string getName() const override { return "model_calvetti_py"; } @@ -554,7 +558,7 @@ class Model_model_calvetti_py : public amici::Model_DAE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "bcedb951ddf674996b269489d74f9b86112038ff"; + return "d587a622b8295dff051b8cb45d009d9b037f7012"; } bool hasQuadraticLLH() const override { diff --git a/models/model_dirac_py/CMakeLists.txt b/models/model_dirac_py/CMakeLists.txt index 5ae41874a9..fcb4addf05 100644 --- a/models/model_dirac_py/CMakeLists.txt +++ b/models/model_dirac_py/CMakeLists.txt @@ -70,6 +70,7 @@ dxdotdp_explicit.h dxdotdx_explicit.cpp dxdotdx_explicit.h dydx.cpp +explicit_roots.cpp h.h model_dirac_py.cpp model_dirac_py.h diff --git a/models/model_dirac_py/explicit_roots.cpp b/models/model_dirac_py/explicit_roots.cpp new file mode 100644 index 0000000000..b0ff22a3e3 --- /dev/null +++ b/models/model_dirac_py/explicit_roots.cpp @@ -0,0 +1,21 @@ +#include "amici/symbolic_functions.h" +#include "amici/defines.h" +#include "sundials/sundials_types.h" + +#include +#include + +#include +#include "p.h" + +namespace amici { +namespace model_model_dirac_py { + +std::vector> explicit_roots_model_dirac_py(const realtype *p, const realtype *k){ + return { + {p2} + }; +} + +} // namespace model_model_dirac_py +} // namespace amici diff --git a/models/model_dirac_py/model_dirac_py.h b/models/model_dirac_py/model_dirac_py.h index 5162fc0637..3d9ef2d9f8 100644 --- a/models/model_dirac_py/model_dirac_py.h +++ b/models/model_dirac_py/model_dirac_py.h @@ -99,7 +99,7 @@ extern void x_solver_model_dirac_py(realtype *x_solver, const realtype *x_rdata) extern std::vector create_splines_model_dirac_py(const realtype *p, const realtype *k); - +extern std::vector> explicit_roots_model_dirac_py(const realtype *p, const realtype *k); /** * @brief AMICI-generated model subclass. */ @@ -123,7 +123,7 @@ class Model_model_dirac_py : public amici::Model_ODE { 0, // nz 0, // nztrue 1, // nevent - 1, // nevent_solver + 0, // nevent_solver 0, // nspl 1, // nobjective 1, // nw @@ -152,8 +152,7 @@ class Model_model_dirac_py : public amici::Model_ODE { std::vector{}, // z2events std::vector{ Event("_E0", true, true, NAN) - }, // events - {} // state-independent events + } // events ) { } @@ -400,6 +399,11 @@ class Model_model_dirac_py : public amici::Model_ODE { void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {} + std::vector> fexplicit_roots(const realtype *p, const realtype *k) override { + return explicit_roots_model_dirac_py(p, k); + } + + std::string getName() const override { return "model_dirac_py"; } @@ -541,7 +545,7 @@ class Model_model_dirac_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "bcedb951ddf674996b269489d74f9b86112038ff"; + return "d587a622b8295dff051b8cb45d009d9b037f7012"; } bool hasQuadraticLLH() const override { diff --git a/models/model_events_py/CMakeLists.txt b/models/model_events_py/CMakeLists.txt index a4bc612bfe..c3f53e961d 100644 --- a/models/model_events_py/CMakeLists.txt +++ b/models/model_events_py/CMakeLists.txt @@ -79,6 +79,7 @@ dxdotdx_explicit.h dydp.cpp dydx.cpp dzdx.cpp +explicit_roots.cpp h.h k.h model_events_py.cpp diff --git a/models/model_events_py/explicit_roots.cpp b/models/model_events_py/explicit_roots.cpp new file mode 100644 index 0000000000..1a3f81a508 --- /dev/null +++ b/models/model_events_py/explicit_roots.cpp @@ -0,0 +1,25 @@ +#include "amici/symbolic_functions.h" +#include "amici/defines.h" +#include "sundials/sundials_types.h" + +#include +#include + +#include +#include "p.h" +#include "k.h" + +namespace amici { +namespace model_model_events_py { + +std::vector> explicit_roots_model_events_py(const realtype *p, const realtype *k){ + return { + {p4}, + {p4}, + {4}, + {4} + }; +} + +} // namespace model_model_events_py +} // namespace amici diff --git a/models/model_events_py/model_events_py.h b/models/model_events_py/model_events_py.h index 7fa56b11e3..00d6a2272c 100644 --- a/models/model_events_py/model_events_py.h +++ b/models/model_events_py/model_events_py.h @@ -99,7 +99,7 @@ extern void x_solver_model_events_py(realtype *x_solver, const realtype *x_rdata extern std::vector create_splines_model_events_py(const realtype *p, const realtype *k); - +extern std::vector> explicit_roots_model_events_py(const realtype *p, const realtype *k); /** * @brief AMICI-generated model subclass. */ @@ -123,7 +123,7 @@ class Model_model_events_py : public amici::Model_ODE { 2, // nz 2, // nztrue 6, // nevent - 4, // nevent_solver + 2, // nevent_solver 0, // nspl 1, // nobjective 1, // nw @@ -157,8 +157,7 @@ class Model_model_events_py : public amici::Model_ODE { Event("Heaviside_3", true, true, NAN), Event("Heaviside_4", true, true, NAN), Event("Heaviside_5", true, true, NAN) - }, // events - {{4.0, {4, 5}}} // state-independent events + } // events ) { } @@ -435,6 +434,11 @@ class Model_model_events_py : public amici::Model_ODE { void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {} + std::vector> fexplicit_roots(const realtype *p, const realtype *k) override { + return explicit_roots_model_events_py(p, k); + } + + std::string getName() const override { return "model_events_py"; } @@ -576,7 +580,7 @@ class Model_model_events_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "bcedb951ddf674996b269489d74f9b86112038ff"; + return "d587a622b8295dff051b8cb45d009d9b037f7012"; } bool hasQuadraticLLH() const override { diff --git a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h index c023f4e367..35681de3e1 100644 --- a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h +++ b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h @@ -150,8 +150,7 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE { amici::SecondOrderMode::none, // o2mode std::vector{1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}, // idlist std::vector{}, // z2events - std::vector{}, // events - {} // state-independent events + std::vector{} // events ) { } @@ -410,6 +409,9 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE { void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {} + std::vector> fexplicit_roots(const realtype *p, const realtype *k) override { return {}; } + + std::string getName() const override { return "model_jakstat_adjoint_py"; } @@ -551,7 +553,7 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "bcedb951ddf674996b269489d74f9b86112038ff"; + return "d587a622b8295dff051b8cb45d009d9b037f7012"; } bool hasQuadraticLLH() const override { diff --git a/models/model_nested_events_py/CMakeLists.txt b/models/model_nested_events_py/CMakeLists.txt index d162a9e239..45629d6d44 100644 --- a/models/model_nested_events_py/CMakeLists.txt +++ b/models/model_nested_events_py/CMakeLists.txt @@ -71,6 +71,7 @@ dxdotdp_explicit.h dxdotdx_explicit.cpp dxdotdx_explicit.h dydx.cpp +explicit_roots.cpp h.h model_nested_events_py.cpp model_nested_events_py.h diff --git a/models/model_nested_events_py/deltaqB.cpp b/models/model_nested_events_py/deltaqB.cpp index f8833fb8ea..e391deef58 100644 --- a/models/model_nested_events_py/deltaqB.cpp +++ b/models/model_nested_events_py/deltaqB.cpp @@ -18,7 +18,7 @@ namespace model_model_nested_events_py { void deltaqB_model_nested_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){ switch(ie) { - case 0: + case 2: switch(ip) { case 1: deltaqB[0] = xB0; diff --git a/models/model_nested_events_py/deltasx.cpp b/models/model_nested_events_py/deltasx.cpp index 7cc2e5c2bf..edd73525fd 100644 --- a/models/model_nested_events_py/deltasx.cpp +++ b/models/model_nested_events_py/deltasx.cpp @@ -20,28 +20,28 @@ namespace model_model_nested_events_py { void deltasx_model_nested_events_py(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old){ switch(ie) { case 0: + case 1: switch(ip) { case 0: + case 1: case 2: case 3: case 4: deltasx[0] = stau0*(dVirusdt - xdot_old0); break; - case 1: - deltasx[0] = stau0*(dVirusdt - xdot_old0) + 1; - break; } break; - case 1: case 2: switch(ip) { case 0: - case 1: case 2: case 3: case 4: deltasx[0] = stau0*(dVirusdt - xdot_old0); break; + case 1: + deltasx[0] = stau0*(dVirusdt - xdot_old0) + 1; + break; } break; } diff --git a/models/model_nested_events_py/deltax.cpp b/models/model_nested_events_py/deltax.cpp index 91b845e21d..148ea2c536 100644 --- a/models/model_nested_events_py/deltax.cpp +++ b/models/model_nested_events_py/deltax.cpp @@ -16,7 +16,7 @@ namespace model_model_nested_events_py { void deltax_model_nested_events_py(double *deltax, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old){ switch(ie) { - case 0: + case 2: deltax[0] = V_0_inject - Virus + x_old0; break; } diff --git a/models/model_nested_events_py/deltaxB.cpp b/models/model_nested_events_py/deltaxB.cpp index 606aa33b08..103e21bad6 100644 --- a/models/model_nested_events_py/deltaxB.cpp +++ b/models/model_nested_events_py/deltaxB.cpp @@ -18,10 +18,10 @@ namespace model_model_nested_events_py { void deltaxB_model_nested_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl){ switch(ie) { - case 1: + case 0: deltaxB[0] = xB0*(dVirusdt - xdot_old0)/(Heaviside_1*Virus*rho_V - Virus*delta_V); break; - case 2: + case 1: deltaxB[0] = -xB0*(dVirusdt - xdot_old0)/(-Heaviside_1*Virus*rho_V + Virus*delta_V); break; } diff --git a/models/model_nested_events_py/explicit_roots.cpp b/models/model_nested_events_py/explicit_roots.cpp new file mode 100644 index 0000000000..b217936b66 --- /dev/null +++ b/models/model_nested_events_py/explicit_roots.cpp @@ -0,0 +1,21 @@ +#include "amici/symbolic_functions.h" +#include "amici/defines.h" +#include "sundials/sundials_types.h" + +#include +#include + +#include +#include "p.h" + +namespace amici { +namespace model_model_nested_events_py { + +std::vector> explicit_roots_model_nested_events_py(const realtype *p, const realtype *k){ + return { + {t_0} + }; +} + +} // namespace model_model_nested_events_py +} // namespace amici diff --git a/models/model_nested_events_py/h.h b/models/model_nested_events_py/h.h index 87eb526963..2145df2f11 100644 --- a/models/model_nested_events_py/h.h +++ b/models/model_nested_events_py/h.h @@ -1,3 +1,3 @@ -#define injection h[0] -#define Heaviside_1 h[1] -#define Heaviside_2 h[2] +#define Heaviside_1 h[0] +#define Heaviside_2 h[1] +#define injection h[2] diff --git a/models/model_nested_events_py/model_nested_events_py.h b/models/model_nested_events_py/model_nested_events_py.h index ed084befb3..8670f38e53 100644 --- a/models/model_nested_events_py/model_nested_events_py.h +++ b/models/model_nested_events_py/model_nested_events_py.h @@ -99,7 +99,7 @@ extern void x_solver_model_nested_events_py(realtype *x_solver, const realtype * extern std::vector create_splines_model_nested_events_py(const realtype *p, const realtype *k); - +extern std::vector> explicit_roots_model_nested_events_py(const realtype *p, const realtype *k); /** * @brief AMICI-generated model subclass. */ @@ -123,7 +123,7 @@ class Model_model_nested_events_py : public amici::Model_ODE { 0, // nz 0, // nztrue 3, // nevent - 3, // nevent_solver + 2, // nevent_solver 0, // nspl 1, // nobjective 1, // nw @@ -151,11 +151,10 @@ class Model_model_nested_events_py : public amici::Model_ODE { std::vector{1.0}, // idlist std::vector{}, // z2events std::vector{ - Event("injection", true, true, NAN), Event("Heaviside_1", true, true, NAN), - Event("Heaviside_2", true, true, NAN) - }, // events - {} // state-independent events + Event("Heaviside_2", true, true, NAN), + Event("injection", true, true, NAN) + } // events ) { } @@ -408,6 +407,11 @@ class Model_model_nested_events_py : public amici::Model_ODE { void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {} + std::vector> fexplicit_roots(const realtype *p, const realtype *k) override { + return explicit_roots_model_nested_events_py(p, k); + } + + std::string getName() const override { return "model_nested_events_py"; } @@ -549,7 +553,7 @@ class Model_model_nested_events_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "bcedb951ddf674996b269489d74f9b86112038ff"; + return "d587a622b8295dff051b8cb45d009d9b037f7012"; } bool hasQuadraticLLH() const override { diff --git a/models/model_nested_events_py/root.cpp b/models/model_nested_events_py/root.cpp index f76cc5c298..68004445bd 100644 --- a/models/model_nested_events_py/root.cpp +++ b/models/model_nested_events_py/root.cpp @@ -13,9 +13,9 @@ namespace amici { namespace model_model_nested_events_py { void root_model_nested_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){ - root[0] = t - t_0; - root[1] = Virus - 1; - root[2] = 1 - Virus; + root[0] = Virus - 1; + root[1] = 1 - Virus; + root[2] = t - t_0; } } // namespace model_model_nested_events_py diff --git a/models/model_nested_events_py/stau.cpp b/models/model_nested_events_py/stau.cpp index a5df0dc8b1..2c81b50be7 100644 --- a/models/model_nested_events_py/stau.cpp +++ b/models/model_nested_events_py/stau.cpp @@ -17,8 +17,12 @@ void stau_model_nested_events_py(realtype *stau, const realtype t, const realtyp switch(ie) { case 0: switch(ip) { + case 0: + case 1: case 2: - stau[0] = -1; + case 3: + case 4: + stau[0] = sx0/(Heaviside_1*Virus*rho_V - Virus*delta_V); break; } break; @@ -29,18 +33,14 @@ void stau_model_nested_events_py(realtype *stau, const realtype t, const realtyp case 2: case 3: case 4: - stau[0] = sx0/(Heaviside_1*Virus*rho_V - Virus*delta_V); + stau[0] = -sx0/(-Heaviside_1*Virus*rho_V + Virus*delta_V); break; } break; case 2: switch(ip) { - case 0: - case 1: case 2: - case 3: - case 4: - stau[0] = -sx0/(-Heaviside_1*Virus*rho_V + Virus*delta_V); + stau[0] = -1; break; } break; diff --git a/models/model_neuron_py/model_neuron_py.h b/models/model_neuron_py/model_neuron_py.h index 143f5de248..385bf1a150 100644 --- a/models/model_neuron_py/model_neuron_py.h +++ b/models/model_neuron_py/model_neuron_py.h @@ -152,8 +152,7 @@ class Model_model_neuron_py : public amici::Model_ODE { std::vector{1}, // z2events std::vector{ Event("event_1", false, true, NAN) - }, // events - {} // state-independent events + } // events ) { } @@ -432,6 +431,9 @@ class Model_model_neuron_py : public amici::Model_ODE { void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {} + std::vector> fexplicit_roots(const realtype *p, const realtype *k) override { return {}; } + + std::string getName() const override { return "model_neuron_py"; } @@ -573,7 +575,7 @@ class Model_model_neuron_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "bcedb951ddf674996b269489d74f9b86112038ff"; + return "d587a622b8295dff051b8cb45d009d9b037f7012"; } bool hasQuadraticLLH() const override { diff --git a/models/model_robertson_py/model_robertson_py.h b/models/model_robertson_py/model_robertson_py.h index 7f188e3849..7c7d021a74 100644 --- a/models/model_robertson_py/model_robertson_py.h +++ b/models/model_robertson_py/model_robertson_py.h @@ -150,8 +150,7 @@ class Model_model_robertson_py : public amici::Model_DAE { amici::SecondOrderMode::none, // o2mode std::vector{1.0, 1.0, 0.0}, // idlist std::vector{}, // z2events - std::vector{}, // events - {} // state-independent events + std::vector{} // events ) { } @@ -394,6 +393,9 @@ class Model_model_robertson_py : public amici::Model_DAE { void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {} + std::vector> fexplicit_roots(const realtype *p, const realtype *k) override { return {}; } + + std::string getName() const override { return "model_robertson_py"; } @@ -535,7 +537,7 @@ class Model_model_robertson_py : public amici::Model_DAE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "bcedb951ddf674996b269489d74f9b86112038ff"; + return "d587a622b8295dff051b8cb45d009d9b037f7012"; } bool hasQuadraticLLH() const override { diff --git a/models/model_steadystate_py/model_steadystate_py.h b/models/model_steadystate_py/model_steadystate_py.h index 289857eb02..0b3f8d89c3 100644 --- a/models/model_steadystate_py/model_steadystate_py.h +++ b/models/model_steadystate_py/model_steadystate_py.h @@ -150,8 +150,7 @@ class Model_model_steadystate_py : public amici::Model_ODE { amici::SecondOrderMode::none, // o2mode std::vector{1.0, 1.0, 1.0}, // idlist std::vector{}, // z2events - std::vector{}, // events - {} // state-independent events + std::vector{} // events ) { } @@ -394,6 +393,9 @@ class Model_model_steadystate_py : public amici::Model_ODE { void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {} + std::vector> fexplicit_roots(const realtype *p, const realtype *k) override { return {}; } + + std::string getName() const override { return "model_steadystate_py"; } @@ -535,7 +537,7 @@ class Model_model_steadystate_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "bcedb951ddf674996b269489d74f9b86112038ff"; + return "d587a622b8295dff051b8cb45d009d9b037f7012"; } bool hasQuadraticLLH() const override { diff --git a/python/sdist/amici/_codegen/cxx_functions.py b/python/sdist/amici/_codegen/cxx_functions.py index b7079c0c55..a90829090d 100644 --- a/python/sdist/amici/_codegen/cxx_functions.py +++ b/python/sdist/amici/_codegen/cxx_functions.py @@ -2,7 +2,7 @@ from __future__ import annotations -from dataclasses import dataclass +from dataclasses import dataclass, field import re @@ -28,6 +28,11 @@ class _FunctionInfo: indicates whether a model-specific implementation is to be generated :ivar body: the actual function body. will be filled later + :ivar default_return_value: + The value to return if the function does not generate a body. + :ivar header: + List of header lines to include in the generated C++ file for this + function such as `#include` statements. """ ode_arguments: str = "" @@ -37,6 +42,8 @@ class _FunctionInfo: sparse: bool = False generate_body: bool = True body: str = "" + default_return_value: str = "" + header: list[str] = field(default_factory=list) def arguments(self, ode: bool = True) -> str: """Get the arguments for the ODE or DAE function""" @@ -143,6 +150,10 @@ def var_in_signature(self, varname: str, ode: bool = True) -> bool: "create_splines": _FunctionInfo( "const realtype *p, const realtype *k", return_type="std::vector", + header=[ + '#include "amici/splinefunctions.h"', + "#include ", + ], ), "spl": _FunctionInfo(generate_body=False), "sspl": _FunctionInfo(generate_body=False), @@ -379,6 +390,14 @@ def var_in_signature(self, varname: str, ode: bool = True) -> bool: "realtype *rz, const int ie, const realtype t, const realtype *x, " "const realtype *p, const realtype *k, const realtype *h" ), + "explicit_roots": _FunctionInfo( + "const realtype *p, const realtype *k", + return_type="std::vector>", + default_return_value="{}", + header=[ + "#include ", + ], + ), } #: list of sparse functions diff --git a/python/sdist/amici/_codegen/model_class.py b/python/sdist/amici/_codegen/model_class.py index 87f070f772..6440fe358c 100644 --- a/python/sdist/amici/_codegen/model_class.py +++ b/python/sdist/amici/_codegen/model_class.py @@ -3,9 +3,6 @@ from __future__ import annotations from .cxx_functions import functions, multiobs_functions -from ..cxxcodeprinter import get_initializer_list - -from ..de_model_components import Event def get_function_extern_declaration(fun: str, name: str, ode: bool) -> str: @@ -72,7 +69,9 @@ def get_model_override_implementation( """ func_info = functions[fun] body = ( - "" + f" return {func_info.default_return_value}; " + if nobody and func_info.default_return_value + else "" if nobody else "\n{ind8}{maybe_return}{fun}_{name}({eval_signature});\n{ind4}".format( ind4=" " * 4, @@ -164,21 +163,3 @@ def remove_argument_types(signature: str) -> str: signature = signature.replace(type_str, "") return signature - - -def get_state_independent_event_intializer(events: list[Event]) -> str: - """Get initializer list for state independent events in amici::Model.""" - map_time_to_event_idx = {} - for event_idx, event in enumerate(events): - if not event.triggers_at_fixed_timepoint(): - continue - trigger_time = float(event.get_trigger_time()) - try: - map_time_to_event_idx[trigger_time].append(event_idx) - except KeyError: - map_time_to_event_idx[trigger_time] = [event_idx] - - return ", ".join( - f"{{{trigger_time}, {get_initializer_list(event_idxs)}}}" - for trigger_time, event_idxs in map_time_to_event_idx.items() - ) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index bf94a8e983..7c6de98d9f 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -47,7 +47,6 @@ get_sunindex_extern_declaration, get_model_override_implementation, get_sunindex_override_implementation, - get_state_independent_event_intializer, ) from ._codegen.template import apply_template from .cxxcodeprinter import ( @@ -433,7 +432,7 @@ def _write_function_file(self, function: str) -> None: ): # Not required. Will create empty function body. equations = sp.Matrix() - elif function == "create_splines": + elif function in ("create_splines", "explicit_roots"): # nothing to do pass else: @@ -442,6 +441,8 @@ def _write_function_file(self, function: str) -> None: # function body if function == "create_splines": body = self._get_create_splines_body() + elif function == "explicit_roots": + body = self._get_explicit_roots_body() else: body = self._get_function_body(function, equations) if not body: @@ -455,6 +456,8 @@ def _write_function_file(self, function: str) -> None: else: lines = [] + func_info = self.functions[function] + # function header lines.extend( [ @@ -465,15 +468,9 @@ def _write_function_file(self, function: str) -> None: "#include ", "#include ", "", + *func_info.header, ] ) - if function == "create_splines": - lines += [ - '#include "amici/splinefunctions.h"', - "#include ", - ] - - func_info = self.functions[function] # extract symbols that need definitions from signature # don't add includes for files that won't be generated. @@ -961,6 +958,37 @@ def _get_create_splines_body(self): body.append("};") return [" " + line for line in body] + def _get_explicit_roots_body(self) -> list[str]: + events = self.model.events() + lines = [] + constant_syms = set(self.model.sym("k")) | set(self.model.sym("p")) + + for event_idx, event in enumerate(events): + if not ( + tigger_times := { + tt + for tt in event.get_trigger_times() + if tt.free_symbols.issubset(constant_syms) + } + ): + continue + + line = ( + " {" + + ", ".join( + self._code_printer.doprint(tt) for tt in tigger_times + ) + + "}," + ) + if event_idx == len(events) - 1: + line = line[:-1] # remove trailing comma for last event + lines.append(line) + + if not lines: + return [] + + return [" return {", *(" " * 4 + line for line in lines), " };"] + def _write_wrapfunctions_cpp(self) -> None: """ Write model-specific 'wrapper' file (``wrapfunctions.cpp``). @@ -1111,9 +1139,6 @@ def event_initializer_list() -> str: ), "EVENT_LIST_INITIALIZER": event_initializer_list(), "Z2EVENT": ", ".join(map(str, self.model._z2event)), - "STATE_INDEPENDENT_EVENTS": get_state_independent_event_intializer( - self.model.events() - ), "ID": ", ".join( str(float(isinstance(s, DifferentialState))) for s in self.model.states() diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index e39bcff0b3..829b99fe56 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -799,8 +799,10 @@ def num_events_solver(self) -> int: :return: number of event symbols (length of the root vector in AMICI) """ + constant_syms = set(self.sym("k")) | set(self.sym("p")) return sum( - not event.triggers_at_fixed_timepoint() for event in self.events() + not event.has_explicit_trigger_times(constant_syms) + for event in self.events() ) def sym(self, name: str) -> sp.Matrix: @@ -1305,12 +1307,17 @@ def parse_events(self) -> None: self.add_component(root) # re-order events - first those that require root tracking, then the others + constant_syms = set(self.sym("k")) | set(self.sym("p")) self._events = list( chain( itertools.filterfalse( - Event.triggers_at_fixed_timepoint, self._events + lambda e: e.has_explicit_trigger_times(constant_syms), + self._events, + ), + filter( + lambda e: e.has_explicit_trigger_times(constant_syms), + self._events, ), - filter(Event.triggers_at_fixed_timepoint, self._events), ) ) diff --git a/python/sdist/amici/de_model_components.py b/python/sdist/amici/de_model_components.py index bfea153009..4f12ec70be 100644 --- a/python/sdist/amici/de_model_components.py +++ b/python/sdist/amici/de_model_components.py @@ -836,13 +836,26 @@ def get_trigger_time(self) -> sp.Float: ) return self._t_root[0] - def has_explicit_trigger_times(self) -> bool: + def has_explicit_trigger_times( + self, allowed_symbols: set[sp.Symbol] | None = None + ) -> bool: """Check whether the event has explicit trigger times. Explicit trigger times do not require root finding to determine the time points at which the event triggers. + + :param allowed_symbols: + The set of symbols that are allowed in the trigger time + expressions. If `None`, any symbols are allowed. + If empty, only numeric values are allowed. """ - return len(self._t_root) > 0 + if allowed_symbols is None: + return len(self._t_root) > 0 + + return len(self._t_root) > 0 and all( + t.is_Number or t.free_symbols.issubset(allowed_symbols) + for t in self._t_root + ) def get_trigger_times(self) -> set[sp.Expr]: """Get the time points at which the event triggers. diff --git a/python/tests/test_events.py b/python/tests/test_events.py index 87b99d3070..132c38df42 100644 --- a/python/tests/test_events.py +++ b/python/tests/test_events.py @@ -735,16 +735,25 @@ def test_handling_of_fixed_time_point_event_triggers(tempdir): """Test handling of events without solver-tracked root functions.""" ant_model = """ model test_events_time_based + one = 1 + two = 2 + three = 3 + four = 4 + five = 5 event_target = 0 bolus = 1 - at (time > 1): event_target = 1 - at (time > 2): event_target = event_target + bolus - at (time > 3): event_target = 3 + at (time > one): event_target = one + at (time > two): event_target = event_target + bolus + at (time > three): event_target = three + at (time > four): event_target = four + at (time > five): event_target = five end """ module_name = "test_events_time_based" antimony2amici( ant_model, + # test with constant parameters and non-constant parameters! + constant_parameters=["four"], model_name=module_name, output_dir=tempdir, ) @@ -752,18 +761,30 @@ def test_handling_of_fixed_time_point_event_triggers(tempdir): module_name=module_name, module_path=tempdir ) amici_model = model_module.getModel() - assert amici_model.ne == 3 + assert amici_model.ne == 5 assert amici_model.ne_solver == 0 - amici_model.setTimepoints(np.linspace(0, 4, 200)) + + amici_model.setTimepoints(np.linspace(0, 10, 20)) amici_solver = amici_model.getSolver() rdata = amici.runAmiciSimulation(amici_model, amici_solver) assert rdata.status == amici.AMICI_SUCCESS assert (rdata.x[rdata.ts < 1] == 0).all() assert (rdata.x[(rdata.ts >= 1) & (rdata.ts < 2)] == 1).all() assert (rdata.x[(rdata.ts >= 2) & (rdata.ts < 3)] == 2).all() - assert (rdata.x[(rdata.ts >= 3)] == 3).all() - - check_derivatives(amici_model, amici_solver, edata=None) + assert (rdata.x[(rdata.ts >= 3) & (rdata.ts < 4)] == 3).all() + assert (rdata.x[(rdata.ts >= 4) & (rdata.ts < 5)] == 4).all() + assert (rdata.x[(rdata.ts >= 5)] == 5).all() + assert rdata.x[-1, :] == 5 + + edata = amici.ExpData(rdata, 1, 0, 0) + + for sens_meth in ( + SensitivityMethod.forward, + SensitivityMethod.adjoint, + ): + amici_solver.setSensitivityMethod(sens_meth) + amici_solver.setSensitivityOrder(SensitivityOrder.first) + check_derivatives(amici_model, amici_solver, edata=edata) @skip_on_valgrind diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index da8639f5d7..e964916100 100644 --- a/src/forwardproblem.cpp +++ b/src/forwardproblem.cpp @@ -73,7 +73,8 @@ void EventHandlingSimulator::run( auto trigger_timepoints = std::ranges::views::filter( trigger_timepoints_tmp, [this, timepoints](auto t) { - return t > ws_->sol.t && t <= timepoints.at(timepoints.size() - 1); + return t > ws_->sol.t && !timepoints.empty() + && t <= timepoints.at(timepoints.size() - 1); } ); auto it_trigger_timepoints = trigger_timepoints.begin(); @@ -122,7 +123,7 @@ void EventHandlingSimulator::run( // if so, set the root-found flag if (ws_->sol.t == next_t_event) { for (auto const ie : - model_->state_independent_events_[ws_->sol.t]) { + model_->explicit_roots_[ws_->sol.t]) { // determine the direction of root crossing from // root function value at the previous event ws_->roots_found[ie] @@ -207,8 +208,7 @@ void EventHandlingSimulator::run_steady_state( // check if we are at a trigger timepoint. // if so, set the root-found flag if (ws_->sol.t == next_t_event) { - for (auto const ie : - model_->state_independent_events_[ws_->sol.t]) { + for (auto const ie : model_->explicit_roots_[ws_->sol.t]) { // determine the direction of root crossing from // root function value at the previous event ws_->roots_found[ie] = std::copysign(1, -ws_->rootvals[ie]); @@ -678,10 +678,7 @@ ForwardProblem::getAdjointUpdates(Model& model, ExpData const& edata) { } SimulationState EventHandlingSimulator::get_simulation_state() { - return { - .sol = ws_->sol, - .mod = model_->getModelState() - }; + return {.sol = ws_->sol, .mod = model_->getModelState()}; } std::vector compute_nroots( diff --git a/src/model.cpp b/src/model.cpp index 91b5a3b3c9..f193931ccb 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -177,13 +177,11 @@ Model::Model( ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode const o2mode, std::vector idlist, std::vector z2event, - std::vector events, - std::map> state_independent_events + std::vector events ) : ModelDimensions(model_dimensions) , o2mode(o2mode) , idlist(std::move(idlist)) - , state_independent_events_(std::move(state_independent_events)) , state_(*this) , derived_state_(*this) , z2event_(std::move(z2event)) @@ -438,6 +436,36 @@ void Model::initEvents( } } } + + // re-compute parameter-dependent but state-independent roots + reinit_explicit_roots(); +} + +void Model::reinit_explicit_roots() { + explicit_roots_.clear(); + + // evaluate timepoints + auto const exp_roots = fexplicit_roots( + state_.unscaledParameters.data(), state_.fixedParameters.data() + ); + Expects(exp_roots.size() == gsl::narrow(ne - ne_solver)); + + // group events by timepoints + for (decltype(exp_roots)::size_type iee = 0; iee < exp_roots.size(); + ++iee) { + // index within all events / root functions (not just explicit ones) + int const ie = ne_solver + gsl::narrow(iee); + auto const& cur_roots = exp_roots[iee]; + Expects(!cur_roots.empty()); + for (auto const& root : cur_roots) { + auto it = explicit_roots_.find(root); + if (it != explicit_roots_.end()) { + it->second.push_back(ie); + } else { + explicit_roots_[root] = {ie}; + } + } + } } int Model::nplist() const { return gsl::narrow(state_.plist.size()); } @@ -3147,13 +3175,11 @@ void Model::fstotal_cl( } std::vector Model::get_trigger_timepoints() const { - std::vector trigger_timepoints( - state_independent_events_.size(), 0.0 - ); + std::vector trigger_timepoints(explicit_roots_.size(), 0.0); // collect keys from state_independent_events_ which are the trigger // timepoints auto it = trigger_timepoints.begin(); - for (auto const& kv : state_independent_events_) { + for (auto const& kv : explicit_roots_) { *(it++) = kv.first; } std::ranges::sort(trigger_timepoints); diff --git a/src/model_header.template.h b/src/model_header.template.h index 154e2a1d4a..118b1be8bb 100644 --- a/src/model_header.template.h +++ b/src/model_header.template.h @@ -99,7 +99,7 @@ TPL_DTOTAL_CLDX_RDATA_ROWVALS_DEF TPL_CREATE_SPLINES_DEF TPL_DSPLINE_VALUESDP_DEF TPL_DSPLINE_SLOPESDP_DEF - +TPL_EXPLICIT_ROOTS_DEF /** * @brief AMICI-generated model subclass. */ @@ -150,8 +150,7 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER { TPL_O2MODE, // o2mode std::vector{TPL_ID}, // idlist std::vector{TPL_Z2EVENT}, // z2events - std::vector{TPL_EVENT_LIST_INITIALIZER}, // events - {TPL_STATE_INDEPENDENT_EVENTS} // state-independent events + std::vector{TPL_EVENT_LIST_INITIALIZER} // events ) { } @@ -283,6 +282,8 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER { TPL_DTOTAL_CLDX_RDATA_COLPTRS_IMPL TPL_DTOTAL_CLDX_RDATA_ROWVALS_IMPL + TPL_EXPLICIT_ROOTS_IMPL + std::string getName() const override { return "TPL_MODELNAME"; } diff --git a/tests/sbml/testSBMLSuite.py b/tests/sbml/testSBMLSuite.py index b9aa163360..163696fee2 100755 --- a/tests/sbml/testSBMLSuite.py +++ b/tests/sbml/testSBMLSuite.py @@ -77,10 +77,13 @@ def test_sbml_testsuite_case(test_id, result_path, sbml_semantic_cases_dir): atol, rtol = apply_settings(settings, solver, model, test_id) - if test_id != "00885": - # 00885: root-after-reinitialization with FSA - solver.setSensitivityOrder(amici.SensitivityOrder.first) - solver.setSensitivityMethod(amici.SensitivityMethod.forward) + solver.setSensitivityOrder(amici.SensitivityOrder.first) + solver.setSensitivityMethod(amici.SensitivityMethod.forward) + + if test_id == "00885": + # 00885: root-after-reinitialization with FSA with default settings + solver.setAbsoluteTolerance(1e-16) + solver.setRelativeTolerance(1e-15) # simulate model rdata = amici.runAmiciSimulation(model, solver)