account for heaviside functions in expressions

documentation/implementation_discontinuities.rst

@@ -72,11 +72,20 @@ respective root function as argument. These will be automatically updated
 during events and take either 0 or 1 values as appropriate pre/post event
 limits.
 
-In order to fully support SBML events, AMICI uses the SUNDIALS functionality to
-only track zero crossings from negative to positive. Accordingly, two root
-functions are necessary to keep track of Heaviside functions and two
-Heaviside function helper variables will be created, where one corresponds
-to the value of `Heaviside(...)` and one to the value of `1-Heaviside(...)`.
+In order to fully support SBML events and Piecewise functions, AMICI uses 
+the SUNDIALS functionality to only track zero crossings from negative to 
+positive. Accordingly, two root functions are necessary to keep track of 
+Heaviside functions and two Heaviside function helper variables will be 
+created, where one corresponds to the value of `Heaviside(...)` and one 
+to the value of `1-Heaviside(...)`. To ensure that Heaviside functions are 
+correctly evaluated at the beginning of the simulation, Heaviside functions 
+are implement as unit steps that evaluate to `1` at `0`. The arguments of 
+Heaviside functions are normalized such that respective properties of 
+Piecewise functions are conserved for the first Heaviside function variable.
+Accordingly, the value of of the second helper variable is incorrect when
+simulation starts when the respective Heaviside function evaluates to zero
+at initialization and should generally not be used.
+
 
 
 Right-Hand-Side Removable Discontinuities

python/amici/ode_export.py

@@ -392,6 +392,15 @@ def get_val(self) -> sp.Expr:
         """
         return self._value
 
+    def set_val(self, val: sp.Expr):
+        """
+        Set ModelQuantity value
+
+        :return:
+            value of the ModelQuantity
+        """
+        self._value = cast_to_sym(val, 'value')
+
 
 class State(ModelQuantity):
     """
@@ -927,7 +936,7 @@ def __init__(self, verbose: Optional[Union[bool, int]] = False,
             'xB': self.num_states_solver,
             'sigmay': self.num_obs,
         }
-        self._eqs: Dict[str, sp.Matrix] = dict()
+        self._eqs: Dict[str, Union[sp.Matrix, List[sp.Matrix]]] = dict()
         self._sparseeqs: Dict[str, Union[sp.Matrix, List[sp.Matrix]]] = dict()
         self._vals: Dict[str, List[float]] = dict()
         self._names: Dict[str, List[str]] = dict()
@@ -1592,6 +1601,9 @@ def parse_events(self) -> None:
         for state in self._states:
             state.set_dt(_process_heavisides(state.get_dt(), roots))
 
+        for expr in self._expressions:
+            expr.set_val(_process_heavisides(expr.get_val(), roots))
+
         # Now add the found roots to the model components
         for root in roots:
             self.add_component(root)
@@ -1792,7 +1804,8 @@ def _compute_equation(self, name: str) -> None:
             # backsubstitution of optimized right hand side terms into RHS
             # calling subs() is costly. Due to looping over events though, the
             # following lines are only evaluated if a model has events
-            tmp_xdot = self._eqs['xdot'].subs(self._syms['w'], self._eqs['w'])
+            tmp_xdot = self._eqs['xdot'].subs(zip(self._syms['w'],
+                                                  self._eqs['w']))
             self._eqs[name] = smart_multiply(self.eq('drootdx'), tmp_xdot) + \
                               self.eq('drootdt')
 
@@ -3685,7 +3698,7 @@ def _process_heavisides(dxdt: sp.Expr, roots: List[Event]) -> sp.Expr:
         # we want unique identifiers for the roots
         tmp_new = _get_unique_root(tmp_old, roots)
         # For Heavisides, we need to add the negative function as well
-        _get_unique_root(sp.sympify(-1 * tmp_old), roots)
+        _get_unique_root(sp.sympify(- tmp_old), roots)
         heavisides.append((sp.Heaviside(tmp_old), tmp_new))
 
     if heavisides:

python/amici/sbml_import.py

@@ -1602,11 +1602,7 @@ def _parse_piecewise_to_heaviside(args: Iterable[sp.Expr]) -> sp.Expr:
         if trigger == sp.false:
             continue
 
-        # we now need to convert the relational >, >=, ... expression into
-        # a trigger function which will be a Heaviside argument
-        root = _parse_trigger(trigger)
-
-        tmp = sp.Heaviside(root)
+        tmp = _parse_trigger(trigger)
         formula += coeff * sp.simplify(not_condition * tmp)
         not_condition *= (1-tmp)
 
@@ -1623,17 +1619,32 @@ def _parse_trigger(trigger: sp.Expr) -> sp.Expr:
     """
     if trigger.is_Relational:
         root = trigger.args[0] - trigger.args[1]
-        if isinstance(trigger, (sp.core.relational.StrictLessThan,
-                                sp.core.relational.LessThan)):
-            root *= -1
-        return root
 
+        # normalize such that we always implement <,
+        # this ensures that we can correctly evaluate the condition if
+        # simulation starts at H(0). This is achieved by translating
+        # conditionals into Heaviside functions H that is implemented as unit
+        # step with H(0) = 1
+        if isinstance(trigger, sp.core.relational.StrictLessThan):
+            # x < y => x - y < 0 => r < 0
+            return 1 - sp.Heaviside(root)
+        if isinstance(trigger, sp.core.relational.LessThan):
+            # x <= y => not(y < x) => not(y - x < 0) => not -r < 0
+            return sp.Heaviside(-root)
+        if isinstance(trigger, sp.core.relational.StrictGreaterThan):
+            # y > x => y - x < 0 => -r < 0
+            return 1 - sp.Heaviside(-root)
+        if isinstance(trigger, sp.core.relational.GreaterThan):
+            # y >= x => not(x < y) => not(x - y < 0) => not r < 0
+            return sp.Heaviside(root)
+
+    # or(x,y) = not(and(not(x),not(y))
     if isinstance(trigger, sp.Or):
-        return sp.Max(*[_parse_trigger(arg)
+        return 1-sp.Mul(*[1-_parse_trigger(arg)
                         for arg in trigger.args])
 
     if isinstance(trigger, sp.And):
-        return sp.Min(*[_parse_trigger(arg)
+        return sp.Mul(*[_parse_trigger(arg)
                         for arg in trigger.args])
 
     raise SBMLException('AMICI can not parse piecewise functions '

src/backwardproblem.cpp

@@ -9,6 +9,7 @@
 #include "amici/misc.h"
 
 #include <cstring>
+#include <cassert>
 
 namespace amici {
 
@@ -89,7 +90,8 @@ void BackwardProblem::workBackwardProblem() {
             }
 
             /* handle discontinuity */
-            if (tnext > model_->getTimepoint(it)) {
+            if (!discs_.empty() && tnext == discs_.back()) {
+                discs_.pop_back();
                 handleEventB();
             }
 
@@ -174,9 +176,18 @@ void BackwardProblem::handleDataPointB(const int it) {
 }
 
 realtype BackwardProblem::getTnext(const int it) {
-    if (discs_.size() > 0 && discs_.back() > model_->getTimepoint(it)) {
+    if (it < 0 && discs_.empty()) {
+        throw AmiException(
+            "No more timepoints (it=%d, ie=%d) available at %f. This should "
+            "not happen, please report a bug including this stacktrace at "
+            "https://github.com/AMICI-dev/AMICI/issues/new/choose",
+            it, discs_.size(), this->t_
+        );
+    }
+
+    if (!discs_.empty() &&
+        (it < 0 || discs_.back() > model_->getTimepoint(it))) {
         double tdisc = discs_.back();
-        discs_.pop_back();
         return tdisc;
     }
 

src/forwardproblem.cpp

@@ -53,6 +53,8 @@ void ForwardProblem::workForwardProblem() {
         model->initialize(x_, dx_, sx_, sdx_,
                           solver->getSensitivityOrder() >=
                           SensitivityOrder::first);
+    else if (model->ne)
+        model->initHeaviside(x_, dx_);
 
     /* compute initial time and setup solver for (pre-)simulation */
     auto t0 = model->t0();
@@ -67,6 +69,8 @@ void ForwardProblem::workForwardProblem() {
                                " is currently not implemented.");
         handlePresimulation();
         t_ = model->t0();
+        if (model->ne)
+            model->initHeaviside(x_, dx_);
     }
     /* when computing adjoint sensitivity analysis with presimulation,
      we need to store sx after the reinitialization after preequilibration

src/model.cpp

@@ -264,12 +264,6 @@ void Model::initHeaviside(AmiVector const &x, AmiVector const &dx) {
     for (int ie = 0; ie < ne; ie++) {
         if (rootvals.at(ie) < 0) {
             state_.h.at(ie) = 0.0;
-        } else if (rootvals.at(ie) == 0) {
-            throw AmiException(
-                "Simulation started in an event. This could lead to "
-                "unexpected results, aborting simulation! Please "
-                "specify an earlier simulation start via "
-                "options.t0");
         } else {
             state_.h.at(ie) = 1.0;
         }

src/symbolic_functions.cpp

@@ -88,7 +88,7 @@ double dirac(double x) {
  *
  */
 double heaviside(double x) {
-    if (x <= 0.0) {
+    if (x < 0.0) {
         return 0.0;
     }
     return 1.0;

tests/benchmark-models/test_benchmark_collection.sh

@@ -16,12 +16,9 @@ Sneyd_PNAS2002
 Zheng_PNAS2012
 Weber_BMC2015"
 
-# Model needs fixing:
-# Chen_MSB2009
 #
 # Not matching reference for unclear reasons
 # Lucarelli_CellSystems2018
-# Weber_BMC2015
 #
 # PEtab needs fixing: Bachmann_MSB2011
 #
@@ -59,6 +56,9 @@ Weber_BMC2015"
 #
 # state-dependent sigmas:
 # Raia_CancerResearch2011
+#
+# Evaluation is known to be inconsistent:
+# Chen_MSB2009
 
 set -e
 

tests/cpputest/unittests/tests1.cpp

@@ -198,7 +198,7 @@ TEST(symbolicFunctions, testSign)
 TEST(symbolicFunctions, testHeaviside)
 {
     CHECK_EQUAL(0, heaviside(-1));
-    CHECK_EQUAL(0, heaviside(0));
+    CHECK_EQUAL(1, heaviside(0));
     CHECK_EQUAL(1, heaviside(1));
 }