...

Author: dweindl

include/amici/model.h

@@ -1362,6 +1362,14 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void updateHeaviside(std::vector<int> const& rootsfound);
 
+    /**
+     * @brief Disable the event with index `ie` because it just triggered.
+     *
+     * Not to be called by user code.
+     *
+     * @param ie Event index.
+     */
+    void disable_event(int const ie) { state_.root_enabled.at(ie) = false; }
     /**
      * @brief Check if the given array has only finite elements.
      *

include/amici/model_state.h

@@ -29,6 +29,7 @@ struct ModelState {
         stotal_cl.resize((dim.nx_rdata - dim.nx_solver) * dim.np, 0.0);
         unscaledParameters.resize(dim.np);
         fixedParameters.resize(dim.nk);
+        root_enabled.resize(dim.ne, true);
     }
 
     /**
@@ -56,14 +57,20 @@ struct ModelState {
      * (dimension: nplist)
      */
     std::vector<int> plist;
+
+    /**
+     * Flags indicating whether a root function element is enabled
+     * (dimension: `ne`)
+     */
+    std::vector<bool> root_enabled;
 };
 
 inline bool operator==(ModelState const& a, ModelState const& b) {
     return is_equal(a.h, b.h) && is_equal(a.total_cl, b.total_cl)
            && is_equal(a.stotal_cl, b.stotal_cl)
            && is_equal(a.unscaledParameters, b.unscaledParameters)
            && is_equal(a.fixedParameters, b.fixedParameters)
-           && a.plist == b.plist;
+           && a.plist == b.plist && a.root_enabled == b.root_enabled;
 }
 
 /**

include/amici/serialization.h

@@ -145,6 +145,7 @@ void serialize(Archive& ar, amici::Model& m, unsigned int const /*version*/) {
     ar & m.state_.unscaledParameters;
     ar & m.state_.fixedParameters;
     ar & m.state_.plist;
+    ar & m.state_.root_enabled;
     ar & m.x0data_;
     ar & m.sx0data_;
     ar & m.nmaxevent_;

python/tests/test_events.py

@@ -1072,3 +1072,47 @@ def test_event_uses_values_from_trigger_time(tempdir):
     )
 
     # TODO: test ASA after https://github.com/AMICI-dev/AMICI/pull/1539
+
+
+def test_simultaneous_events(tempdir):
+    """Test simultaneously firing events with different trigger functions."""
+    from amici.antimony_import import antimony2amici
+
+    model_name = "test_simultaneous_events"
+    antimony2amici(
+        r"""
+        target1_0 = 1
+        target1 = target1_0
+        one = 1
+        target1' = one
+        two = 2
+        target2_0 = two
+        target2 = target2_0
+        target2' = 1
+        some_time = time
+        some_time' = 1
+        trigger_time = 1000
+
+        E1: at some_time >= trigger_time, priority=10, fromTrigger=false:
+            target1 = target1 + 10;
+        E2: at time >= trigger_time, priority=20, fromTrigger=false:
+            target2 = target2 + 10;
+        """,
+        model_name=model_name,
+        output_dir=tempdir,
+    )
+
+    model_module = import_model_module(model_name, tempdir)
+
+    model = model_module.get_model()
+    model.setTimepoints([0, 2])
+    solver = model.getSolver()
+    solver.setRelativeTolerance(1e-6)
+    solver.setAbsoluteTolerance(1e-6)
+    solver.setSensitivityOrder(SensitivityOrder.first)
+    solver.setSensitivityMethod(SensitivityMethod.forward)
+
+    rdata = amici.runAmiciSimulation(model, solver)
+    assert rdata.status == amici.AMICI_SUCCESS
+    assert_allclose(rdata.by_id("target1"), [1.0, 13.0])
+    assert_allclose(rdata.by_id("target2"), [2.0, 14.0])

src/forwardproblem.cpp

@@ -359,6 +359,7 @@ void EventHandlingSimulator::handle_event(
                         ? std::optional<SimulationState>(get_simulation_state())
                         : std::nullopt)}
             );
+            model_->disable_event(ie);
         }
     }
 

src/model.cpp

@@ -437,6 +437,7 @@ void Model::initEvents(
                 roots_found.at(ie) = 1;
             }
         }
+        state_.root_enabled[ie] = rootvals.at(ie) < 0;
     }
 }
 

src/model_ode.cpp

@@ -97,6 +97,19 @@ void Model_ODE::froot(realtype t, const_N_Vector x, gsl::span<realtype> root) {
         state_.unscaledParameters.data(), state_.fixedParameters.data(),
         state_.h.data(), state_.total_cl.data()
     );
+
+    for (int ie = 0; ie < ne; ++ie) {
+        if (!state_.root_enabled[ie]) {
+            if (root[ie] < 0.0) {
+                // If the disabled root function becomes negative,
+                // re-enable it.
+                state_.root_enabled[ie] = true;
+            } else {
+                // If the root function is disabled, mask it
+                root[ie] = 1.0;
+            }
+        }
+    }
 }
 
 void Model_ODE::fxdot(