Handle state-independent, parameter-dependent event triggers without root-finding

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<std::vector<realtype>> fexplicit_roots(
+        [[maybe_unused]] realtype const* p, [[maybe_unused]] realtype const* k
+    ) = 0;
 };
 
 } // namespace amici

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<realtype> idlist, std::vector<int> z2event,
-        std::vector<Event> events = {},
-        std::map<realtype, std::vector<int>> state_independent_events = {}
+        std::vector<Event> events = {}
     );
 
     /** Destructor. */
@@ -301,6 +298,13 @@ class Model : public AbstractModel, public ModelDimensions {
         std::vector<int>& 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<std::vector<realtype>> 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<realtype, std::vector<int>> state_independent_events_ = {};
+    std::map<realtype, std::vector<int>> explicit_roots_ = {};
 
   protected:
     /**

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<realtype> const& idlist,
-        std::vector<int> const& z2event, std::vector<Event> events = {},
-        std::map<realtype, std::vector<int>> state_independent_events = {}
+        std::vector<int> const& z2event, std::vector<Event> 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);

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<realtype> const& idlist,
-        std::vector<int> const& z2event, std::vector<Event> events = {},
-        std::map<realtype, std::vector<int>> state_independent_events = {}
+        std::vector<int> const& z2event, std::vector<Event> events = {}
     )
         : Model(
               model_dimensions, simulation_parameters, o2mode, idlist, z2event,
-              events, state_independent_events
+              events
           ) {}
 
     void

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_;
 }
 

include/amici/sundials_matrix_wrapper.h

@@ -591,7 +591,6 @@ class SUNMatrixWrapper {
     bool ownmat = true;
 };
 
-
 /**
  * @brief Output formatter for SUNMatrixWrapper.
  * @param os output stream

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 

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 <gsl/gsl-lite.hpp>
+#include <algorithm>
+
+#include <vector>
+#include "k.h"
+
+namespace amici {
+namespace model_model_calvetti_py {
+
+std::vector<std::vector<realtype>> explicit_roots_model_calvetti_py(const realtype *p, const realtype *k){
+    return {
+        {10},
+        {10},
+        {12},
+        {12}
+    };
+}
+
+} // namespace model_model_calvetti_py
+} // namespace amici

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<HermiteSpline> create_splines_model_calvetti_py(const realtype *p, const realtype *k);
 
 
-
+extern std::vector<std::vector<realtype>> 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<std::vector<realtype>> 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 {

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 

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 <gsl/gsl-lite.hpp>
+#include <algorithm>
+
+#include <vector>
+#include "p.h"
+
+namespace amici {
+namespace model_model_dirac_py {
+
+std::vector<std::vector<realtype>> explicit_roots_model_dirac_py(const realtype *p, const realtype *k){
+    return {
+        {p2}
+    };
+}
+
+} // namespace model_model_dirac_py
+} // namespace amici

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<HermiteSpline> create_splines_model_dirac_py(const realtype *p, const realtype *k);
 
 
-
+extern std::vector<std::vector<realtype>> 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<int>{},               // z2events
               std::vector<Event>{
                   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<std::vector<realtype>> 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 {

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 

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 <gsl/gsl-lite.hpp>
+#include <algorithm>
+
+#include <vector>
+#include "p.h"
+#include "k.h"
+
+namespace amici {
+namespace model_model_events_py {
+
+std::vector<std::vector<realtype>> explicit_roots_model_events_py(const realtype *p, const realtype *k){
+    return {
+        {p4},
+        {p4},
+        {4},
+        {4}
+    };
+}
+
+} // namespace model_model_events_py
+} // namespace amici

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<HermiteSpline> create_splines_model_events_py(const realtype *p, const realtype *k);
 
 
-
+extern std::vector<std::vector<realtype>> 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<std::vector<realtype>> 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 {