initial implementation, fixes #1760

include/amici/forwardproblem.h

@@ -290,9 +290,11 @@ class ForwardProblem {
      *
      * @param tlastroot pointer to the timepoint of the last event
      * @param seflag Secondary event flag
+     * @param initial_event initial event flag
      */
 
-    void handleEvent(realtype *tlastroot,bool seflag);
+    void handleEvent(realtype *tlastroot, bool seflag,
+                     bool initial_event);
 
     /**
      * @brief Extract output information for events

include/amici/model.h

@@ -88,6 +88,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param o2mode Second order sensitivity mode
      * @param idlist Indexes indicating algebraic components (DAE only)
      * @param z2event Mapping of event outputs to events
+     * @param r0 Root initial values
      * @param pythonGenerated Flag indicating matlab or python wrapping
      * @param ndxdotdp_explicit Number of nonzero elements in `dxdotdp_explicit`
      * @param ndxdotdx_explicit Number of nonzero elements in `dxdotdx_explicit`
@@ -97,7 +98,8 @@ class Model : public AbstractModel, public ModelDimensions {
           SimulationParameters simulation_parameters,
           amici::SecondOrderMode o2mode,
           std::vector<amici::realtype> idlist,
-          std::vector<int> z2event, bool pythonGenerated = false,
+          std::vector<int> z2event, std::vector<bool> r0,
+          bool pythonGenerated = false,
           int ndxdotdp_explicit = 0, int ndxdotdx_explicit = 0,
           int w_recursion_depth = 0);
 
@@ -202,9 +204,11 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param sdx Reference to time derivative of state sensitivities (DAE only)
      * @param computeSensitivities Flag indicating whether sensitivities are to
      * be computed
+     * @param roots_found boolean indicators indicating whether roots were found at t0 by this fun
      */
     void initialize(AmiVector &x, AmiVector &dx, AmiVectorArray &sx,
-                    AmiVectorArray &sdx, bool computeSensitivities);
+                    AmiVectorArray &sdx, bool computeSensitivities,
+                    std::vector<int> &roots_founds);
 
     /**
      * @brief Initialize model properties.
@@ -236,8 +240,10 @@ class Model : public AbstractModel, public ModelDimensions {
      *
      * @param x Reference to state variables
      * @param dx Reference to time derivative of states (DAE only)
+     * @param roots_found boolean indicators indicating whether roots were found at t0 by this fun
      */
-    void initHeaviside(const AmiVector &x, const AmiVector &dx);
+    void initEvents(const AmiVector &x, const AmiVector &dx,
+                    std::vector<int> &roots_founds);
 
     /**
      * @brief Get number of parameters wrt to which sensitivities are computed.
@@ -1864,6 +1870,11 @@ class Model : public AbstractModel, public ModelDimensions {
     /** vector of bools indicating whether state variables are to be assumed to
      * be positive */
     std::vector<bool> state_is_non_negative_;
+
+    /** Vector of booleans indicating the initial boolean value for every event trigger function. Events at t0
+     * can only trigger if the initial value is set to False/ Must be specified during model compilation by
+     * setting the `initialValue` attribute of an event trigger. */
+    std::vector<bool> root_initial_values_;
 
     /** boolean indicating whether any entry in stateIsNonNegative is `true` */
     bool any_state_non_negative_ {false};

python/amici/ode_export.py

@@ -2937,6 +2937,11 @@ def _write_model_header_cpp(self) -> None:
             'W_RECURSION_DEPTH': self.model._w_recursion_depth,
             'QUADRATIC_LLH': 'true'
                 if self.model._has_quadratic_nllh else 'false',
+            'ROOT_INITIAL_VALUES':
+                ','.join([
+                    'true' if event.get0() else 'false'
+                    for event in self.model._events
+                ])
         }
 
         for func_name, func_info in self.functions.items():

python/amici/ode_model.py

@@ -2,50 +2,27 @@
 
 
 import sympy as sp
-import numpy as np
-import re
-import shutil
-import subprocess
-import sys
-import os
-import copy
 import numbers
-import logging
-import itertools
-import contextlib
 
 try:
     import pysb
 except ImportError:
     pysb = None
 
 from typing import (
-    Callable, Optional, Union, List, Dict, Tuple, SupportsFloat, Sequence,
-    Set, Any
+    Optional, Union, Dict, SupportsFloat, Set
 )
-from dataclasses import dataclass
-from string import Template
-from sympy.matrices.immutable import ImmutableDenseMatrix
-from sympy.matrices.dense import MutableDenseMatrix
-from sympy.logic.boolalg import BooleanAtom
-from itertools import chain
-from .cxxcodeprinter import AmiciCxxCodePrinter, get_switch_statement
-
-from . import (
-    amiciSwigPath, amiciSrcPath, amiciModulePath, __version__, __commit__,
-    sbml_import
-)
-from .logging import get_logger, log_execution_time, set_log_level
-from .constants import SymbolId
-from .import_utils import smart_subs_dict, toposort_symbols, \
-    ObservableTransformation, generate_measurement_symbol, RESERVED_SYMBOLS
+
+from .import_utils import ObservableTransformation, \
+    generate_measurement_symbol, RESERVED_SYMBOLS
 from .import_utils import cast_to_sym
 
 __all__ = [
     'ConservationLaw', 'Constant', 'Event', 'Expression', 'LogLikelihood',
     'ModelQuantity', 'Observable', 'Parameter', 'SigmaY', 'State'
 ]
 
+
 class ModelQuantity:
     """
     Base class for model components
@@ -166,14 +143,6 @@ def __init__(self,
         self._ncoeff: sp.Expr = coefficients[state_id]
         super(ConservationLaw, self).__init__(identifier, name, value)
 
-    def get_state(self) -> sp.Symbol:
-        """
-        Get the identifier of the state that this conservation law replaces
-
-        :return: identifier of the state
-        """
-        return self._state_id
-
     def get_ncoeff(self, state_id) -> Union[sp.Expr, int, float]:
         """
         Computes the normalized coefficient a_i/a_j where i is the index of
@@ -211,10 +180,6 @@ class State(ModelQuantity):
         algebraic formula that defines the temporal derivative of this state
 
     """
-
-    _dt: Union[sp.Expr, None] = None
-    _conservation_law: Union[sp.Expr, None] = None
-
     def __init__(self,
                  identifier: sp.Symbol,
                  name: str,
@@ -276,7 +241,7 @@ def get_dt(self) -> sp.Expr:
         """
         return self._dt
 
-    def get_free_symbols(self) -> Set[sp.Symbol]:
+    def get_free_symbols(self) -> Set[sp.Basic]:
         """
         Gets the set of free symbols in time derivative and initial conditions
 
@@ -307,8 +272,9 @@ def get_x_rdata(self):
 
     def get_dx_rdata_dx_solver(self, state_id):
         """
-        Returns the expression that allows computation of ``dx_rdata_dx_solver`` for this
-        state, accounting for conservation laws.
+        Returns the expression that allows computation of
+        ``dx_rdata_dx_solver`` for this state, accounting for conservation
+        laws.
 
         :return: dx_rdata_dx_solver expression
         """
@@ -514,7 +480,8 @@ def __init__(self,
                  name: str,
                  value: sp.Expr,
                  state_update: Union[sp.Expr, None],
-                 event_observable: Union[sp.Expr, None]):
+                 event_observable: Union[sp.Expr, None],
+                 initial_value: Optional[bool] = True):
         """
         Create a new Event instance.
 
@@ -534,15 +501,29 @@ def __init__(self,
         :param event_observable:
             formula a potential observable linked to the event
             (None for Heaviside functions, empty events without observable)
+
+        :param initial_value:
+            initial boolean value of the trigger function at t0. If set to
+            `False`, events may trigger at t==0, otherwise not.
         """
         super(Event, self).__init__(identifier, name, value)
         # add the Event specific components
         self._state_update = state_update
         self._observable = event_observable
+        self._root0 = initial_value
+
+    def get0(self) -> bool:
+        """
+        Return the initial value for the root function.
+        :return:
+            initial value formula
+        """
+        return self._root0
 
     def __eq__(self, other):
         """
         Check equality of events at the level of trigger/root functions, as we
         need to collect unique root functions for ``roots.cpp``
         """
-        return self.get_val() == other.get_val()
+        return self.get_val() == other.get_val() and \
+            (self.get0() == other.get0())

python/amici/sbml_import.py

@@ -484,13 +484,6 @@ def check_event_support(self) -> None:
             if trigger_sbml.getMath() is None:
                 logger.warning(f'Event {event_id} trigger has no trigger '
                                'expression, so a dummy trigger will be set.')
-            if not trigger_sbml.getInitialValue():
-                # True: event not executed if triggered at time == 0
-                # (corresponding to AMICI default). Raise if set to False.
-                raise SBMLException(
-                    f'Event {event_id} has a trigger that has an initial '
-                    'value of False. This is currently not supported in AMICI.'
-                )
 
             if not trigger_sbml.getPersistent():
                 raise SBMLException(
@@ -1140,6 +1133,7 @@ def get_empty_bolus_value() -> sp.Float:
                 'value': trigger,
                 'state_update': sp.MutableDenseMatrix(bolus),
                 'event_observable': None,
+                'root0': trigger_sbml.getInitialValue(),
             }
 
     @log_execution_time('processing SBML observables', logger)

src/forwardproblem.cpp

@@ -52,15 +52,21 @@ void ForwardProblem::workForwardProblem() {
     if (!preequilibrated_)
         model->initialize(x_, dx_, sx_, sdx_,
                           solver->getSensitivityOrder() >=
-                          SensitivityOrder::first);
-    else if (model->ne)
-        model->initHeaviside(x_, dx_);
+                          SensitivityOrder::first,
+                          roots_found_);
+    else if (model->ne) {
+        model->initEvents(x_, dx_, roots_found_);
+    }
 
     /* compute initial time and setup solver for (pre-)simulation */
     auto t0 = model->t0();
     if (presimulate)
         t0 -= edata->t_presim;
     solver->setup(t0, model, x_, dx_, sx_, sdx_);
+
+    if (model->ne && std::any_of(roots_found_.begin(), roots_found_.end(),
+                                 [](int rf){return rf==1;}))
+        handleEvent(&t0, false, true);
 
     /* perform presimulation if necessary */
     if (presimulate) {
@@ -69,9 +75,14 @@ void ForwardProblem::workForwardProblem() {
                                " is currently not implemented.");
         handlePresimulation();
         t_ = model->t0();
-        if (model->ne)
-            model->initHeaviside(x_, dx_);
+        if (model->ne) {
+            model->initEvents(x_, dx_, roots_found_);
+            if (std::any_of(roots_found_.begin(), roots_found_.end(),
+                            [](int rf){return rf==1;}))
+                handleEvent(&t0, false, true);
+        }
     }
+
     /* when computing adjoint sensitivity analysis with presimulation,
      we need to store sx after the reinitialization after preequilibration
      but before reinitialization after presimulation. As presimulation with ASA
@@ -92,7 +103,6 @@ void ForwardProblem::workForwardProblem() {
     if (solver->computingFSA() || (solver->computingASA() && !presimulate ))
         sx_ = solver->getStateSensitivity(model->t0());
 
-
     /* store initial state and sensitivity*/
     initial_state_ = getSimulationState();
 
@@ -114,7 +124,7 @@ void ForwardProblem::workForwardProblem() {
                     /* clustering of roots => turn off rootfinding */
                     solver->turnOffRootFinding();
                 } else if (status == AMICI_ROOT_RETURN) {
-                    handleEvent(&tlastroot_, false);
+                    handleEvent(&tlastroot_, false, false);
                 }
             }
         }
@@ -138,22 +148,23 @@ void ForwardProblem::handlePresimulation()
 }
 
 
-void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
+void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag,
+                                 const bool initial_event) {
     /* store Heaviside information at event occurrence */
     model->froot(t_, x_, dx_, rootvals_);
 
     /* store timepoint at which the event occurred*/
     discs_.push_back(t_);
 
     /* extract and store which events occurred */
-    if (!seflag) {
+    if (!seflag && !initial_event) {
         solver->getRootInfo(roots_found_.data());
     }
     root_idx_.push_back(roots_found_);
 
     rval_tmp_ = rootvals_;
 
-    if (!seflag) {
+    if (!seflag && !initial_event) {
         /* only check this in the first event fired, otherwise this will always
          * be true */
         if (t_ == *tlastroot) {
@@ -165,39 +176,42 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
         *tlastroot = t_;
     }
 
-    if(model->nz>0)
+    if(model->nz > 0)
         storeEvent();
 
     /* if we need to do forward sensitivities later on we need to store the old
      * x and the old xdot */
     if (solver->getSensitivityOrder() >= SensitivityOrder::first) {
         /* store x and xdot to compute jump in sensitivities */
-        x_old_ = x_;
+        x_old_.copy(x_);
     }
     if (solver->computingFSA()) {
         model->fxdot(t_, x_, dx_, xdot_);
-        xdot_old_ = xdot_;
-        dx_old_ = dx_;
+        xdot_old_.copy(xdot_);
+        dx_old_.copy(dx_);
         /* compute event-time derivative only for primary events, we get
          * into trouble with multiple simultaneously firing events here (but
          * is this really well defined then?), in that case just use the
          * last ie and hope for the best. */
-        if (!seflag) {
+        if (!seflag && !initial_event) {
             for (int ie = 0; ie < model->ne; ie++) {
                 if (roots_found_.at(ie) == 1) {
                     /* only consider transitions false -> true */
                     model->getEventTimeSensitivity(stau_, t_, ie, x_, sx_);
                 }
             }
         }
+        if (initial_event) // t0 has no parameter dependency
+            std::fill(stau_.begin(), stau_.end(), 0.0);
     } else if (solver->computingASA()) {
         /* store x to compute jump in discontinuity */
         x_disc_.push_back(x_);
         xdot_disc_.push_back(xdot_);
         xdot_old_disc_.push_back(xdot_old_);
     }
 
-    model->updateHeaviside(roots_found_);
+    if (!initial_event)
+        model->updateHeaviside(roots_found_);
 
     applyEventBolus();
 
@@ -239,7 +253,7 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
                                  "Secondary event was triggered. Depending on "
                                  "the bolus of the secondary event, forward "
                                  "sensitivities can be incorrect.");
-        handleEvent(tlastroot, true);
+        handleEvent(tlastroot, true, false);
     }
 
     /* only reinitialise in the first event fired */