fix initial events, fixes #1760 (#1789)

* initial implementation, fixes #1760

* fixup

* fixup

* fixup

* fix doc

* fix doc

* fixup swig

* fixup sbml import

* add sensitivity check, fixup time-independent triggers

* fixup

* fixup

* fixup

* fixup parameter event assignments

* update test stats

* Apply suggestions from code review

Co-authored-by: Daniel Weindl <[email protected]>

Co-authored-by: Daniel Weindl <[email protected]>

Author: Fabian Fröhlich

documentation/python_interface.rst

@@ -26,7 +26,7 @@ AMICI can import :term:`SBML` models via the
 Status of SBML support in Python-AMICI
 ++++++++++++++++++++++++++++++++++++++
 
-Python-AMICI currently **passes 1014 out of the 1821 (~56%) test cases** from
+Python-AMICI currently **passes 1030 out of the 1821 (~57%) test cases** from
 the semantic
 `SBML Test Suite <https://github.com/sbmlteam/sbml-test-suite/>`_
 (`current status <https://github.com/AMICI-dev/AMICI/actions>`_).

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

@@ -97,7 +97,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,
+          bool pythonGenerated = false,
           int ndxdotdp_explicit = 0, int ndxdotdx_explicit = 0,
           int w_recursion_depth = 0);
 
@@ -202,9 +203,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_found);
 
     /**
      * @brief Initialize model properties.
@@ -236,8 +239,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_found);
 
     /**
      * @brief Get number of parameters wrt to which sensitivities are computed.
@@ -1864,6 +1869,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

@@ -822,7 +822,7 @@ def transform_dxdt_to_concentration(species_id, dxdt):
             else:
                 args += ['value']
             if symbol_name == SymbolId.EVENT:
-                args += ['state_update', 'event_observable']
+                args += ['state_update', 'event_observable', 'initial_value']
             if symbol_name == SymbolId.OBSERVABLE:
                 args += ['transformation']
 
@@ -1589,34 +1589,43 @@ def _compute_equation(self, name: str) -> None:
         elif name == 'stau':
             self._eqs[name] = [
                 -self.eq('sroot')[ie, :] / self.eq('drootdt_total')[ie]
+                if not self.eq('drootdt_total')[ie].is_zero else
+                sp.zeros(*self.eq('sroot')[ie, :].shape)
                 for ie in range(self.num_events())
             ]
 
         elif name == 'deltasx':
             event_eqs = []
             for ie, event in enumerate(self._events):
-                if event._state_update is not None:
-                    # ====== chain rule for the state variables ===============
-                    # get xdot with expressions back-substituted
-                    tmp_eq = smart_multiply(
+
+                tmp_eq = sp.zeros(self.num_states_solver(), self.num_par())
+
+                # only add stau part if trigger is time-dependent
+                if not self.eq('drootdt_total')[ie].is_zero:
+                    tmp_eq += smart_multiply(
                         (self.sym('xdot_old') - self.sym('xdot')),
                         self.eq('stau')[ie])
-                    # construct an enhanced state sensitivity, which accounts
-                    # for the time point sensitivity as well
+
+                # only add deltax part if there is state update
+                if event._state_update is not None:
+                    # partial derivative for the parameters
+                    tmp_eq += self.eq('ddeltaxdp')[ie]
+
+                    # initial part of chain rule state variables
                     tmp_dxdp = self.sym('sx') * sp.ones(1, self.num_par())
-                    tmp_dxdp += smart_multiply(self.sym('xdot'),
-                                               self.eq('stau')[ie])
+
+                    # only add stau part if trigger is time-dependent
+                    if not self.eq('drootdt_total')[ie].is_zero:
+                        # chain rule for the time point
+                        tmp_eq += smart_multiply(self.eq('ddeltaxdt')[ie],
+                                                 self.eq('stau')[ie])
+
+                        # additional part of chain rule state variables
+                        tmp_dxdp += smart_multiply(self.sym('xdot'),
+                                                   self.eq('stau')[ie])
+                    # finish chain rule for the state variables
                     tmp_eq += smart_multiply(self.eq('ddeltaxdx')[ie],
                                              tmp_dxdp)
-                    # ====== chain rule for the time point ====================
-                    tmp_eq += smart_multiply(self.eq('ddeltaxdt')[ie],
-                                             self.eq('stau')[ie])
-                    # ====== partial derivative for the parameters ============
-                    tmp_eq += self.eq('ddeltaxdp')[ie]
-                else:
-                    tmp_eq = smart_multiply(
-                        (self.eq('xdot_old') - self.eq('xdot')),
-                        self.eq('stau')[ie])
 
                 event_eqs.append(tmp_eq)
 
@@ -2937,6 +2946,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.get_initial_value() 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,30 @@ 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==t0``, otherwise not.
         """
         super(Event, self).__init__(identifier, name, value)
         # add the Event specific components
         self._state_update = state_update
         self._observable = event_observable
+        self._initial_value = initial_value
+
+    def get_initial_value(self) -> bool:
+        """
+        Return the initial value for the root function.
+
+        :return:
+            initial value formula
+        """
+        return self._initial_value
 
     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.get_initial_value() == other.get_initial_value())

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(
@@ -1001,10 +994,16 @@ def _convert_event_assignment_parameter_targets_to_species(self):
                     parameter_def = \
                         self.symbols[symbol_id].pop(parameter_target)
             if parameter_def is None:
-                raise AssertionError(
-                    'Unexpected error. The parameter target of an event '
-                    'assignment could not be found.'
+                # this happens for parameters that have initial assignments
+                # or are assignment rule targets
+                par = self.sbml.getElementBySId(str(parameter_target))
+                ia_init = self._get_element_initial_assignment(
+                    par.getId()
                 )
+                parameter_def = {
+                    'name': par.getName() if par.isSetName() else par.getId(),
+                    'value': par.getValue() if ia_init is None else ia_init
+                }
             # Fixed parameters are added as species such that they can be
             # targets of events.
             self.symbols[SymbolId.SPECIES][parameter_target] = {
@@ -1140,6 +1139,9 @@ def get_empty_bolus_value() -> sp.Float:
                 'value': trigger,
                 'state_update': sp.MutableDenseMatrix(bolus),
                 'event_observable': None,
+                'initial_value':
+                    trigger_sbml.getInitialValue() if trigger_sbml is not None
+                    else True,
             }
 
     @log_execution_time('processing SBML observables', logger)