Release 0.30.1

CHANGELOG.md

@@ -4,6 +4,27 @@ See also our [versioning policy](https://amici.readthedocs.io/en/latest/versioni
 
 ## v0.X Series
 
+### v0.30.1 (2025-02-18)
+
+Bugfix-only release.
+
+* Removed `eqx.debug.nan`, fixes #2629
+  by @FFroehlich in https://github.com/AMICI-dev/AMICI/pull/2630
+* Fixes an SBML import issue that led to incorrect results for models with
+  species-dependent initial assignments (fixes #2642)
+  by @FFroehlich in https://github.com/AMICI-dev/AMICI/pull/2643
+* Fixed `CVodeGetSensDky` error message
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2644
+* Disabled `CvodeF` checkpointing to prevent certain rare crashes when forward
+  integration takes exactly `maxsteps` integration steps, plus some additional
+  yet unclear condition.
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2645
+* Fixed rare crashes due to uncaught exceptions in `~FinalStateStorer`
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2647
+
+**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.30.0...v0.30.1
+
+
 ### v0.30.0 (2024-12-10)
 
 *Please note that the amici JAX model generation introduced in v0.29.0 is

include/amici/forwardproblem.h

@@ -441,19 +441,37 @@
     /**
      * @brief destructor, stores simulation state
      */
-    ~FinalStateStorer() {
+    ~FinalStateStorer() noexcept(false) {
         if (fwd_) {
-            fwd_->final_state_ = fwd_->getSimulationState();
-            // if there is an associated output timepoint, also store it in
-            // timepoint_states if it's not present there.
-            // this may happen if there is an error just at
-            // (or indistinguishably before) an output timepoint
-            auto final_time = fwd_->getFinalTime();
-            auto const timepoints = fwd_->model->getTimepoints();
-            if (!fwd_->timepoint_states_.count(final_time)
-                && std::find(timepoints.cbegin(), timepoints.cend(), final_time)
-                       != timepoints.cend()) {
-                fwd_->timepoint_states_[final_time] = fwd_->final_state_;
+            try {
+                // This may throw in `CVodeSolver::getSens`
+                // due to https://github.com/LLNL/sundials/issues/82.
+                // Therefore, this dtor must be `noexcept(false)` to avoid
+                // programm termination.
+                fwd_->final_state_ = fwd_->getSimulationState();
+                // if there is an associated output timepoint, also store it in
+                // timepoint_states if it's not present there.
+                // this may happen if there is an error just at
+                // (or indistinguishably before) an output timepoint
+                auto final_time = fwd_->getFinalTime();
+                auto const timepoints = fwd_->model->getTimepoints();
+                if (!fwd_->timepoint_states_.count(final_time)
+                    && std::find(timepoints.cbegin(), timepoints.cend(), final_time)
+                           != timepoints.cend()) {
+                    fwd_->timepoint_states_[final_time] = fwd_->final_state_;
+                }
+            } catch (std::exception const&) {
+                // We must not throw in case we are already in the stack
+                // unwinding phase due to some other active exception, otherwise
+                // this will also lead to termination.
+                //
+                // In case there is another active exception,
+                // `fwd_->{final_state_,timepoint_states_}` won't be set,
+                // and we assume that they are either not accessed anymore, or
+                // that there is appropriate error handling in place.
+                if(!std::uncaught_exceptions()) {
+                    throw;
+                }
             }
         }
     }

python/sdist/amici/jax/petab.py

@@ -500,7 +500,7 @@ def run_simulation(
             simulation_condition[0], p
         )
         return self.model.simulate_condition(
-            p=eqx.debug.backward_nan(p),
+            p=p,
             ts_init=jax.lax.stop_gradient(jnp.array(ts_preeq)),
             ts_dyn=jax.lax.stop_gradient(jnp.array(ts_dyn)),
             ts_posteq=jax.lax.stop_gradient(jnp.array(ts_posteq)),

python/sdist/amici/sbml_import.py

@@ -2292,9 +2292,29 @@
 
         sym_math, rateof_to_dummy = _rateof_to_dummy(sym_math)
 
-        for species_id, species in self.symbols[SymbolId.SPECIES].items():
-            if "init" in species:
-                sym_math = smart_subs(sym_math, species_id, species["init"])
+        # we can't rely on anything else being properly initialized at this point, so we need to
+        # compute all initial values from scratch, recursively
+        for var in sym_math.free_symbols:
+            element_id = str(var)
+            # already recursive since _get_element_initial_assignment calls _make_initial
+            if (
+                ia := self._get_element_initial_assignment(element_id)
+            ) is not None:
+                sym_math = sym_math.subs(var, ia)
+            elif (species := self.sbml.getSpecies(element_id)) is not None:
+                # recursive!
+                init = self._make_initial(get_species_initial(species))
+                sym_math = sym_math.subs(var, init)
+            elif var in self.symbols[SymbolId.SPECIES]:
+                sym_math = sym_math.subs(
+                    var, self.symbols[SymbolId.SPECIES][var]["init"]
+                )
+            elif (
+                element := self.sbml.getElementBySId(element_id)
+            ) and self.is_rate_rule_target(element):
+                # no need to recurse here, as value is numeric
+                init = sp.Float(element.getValue())
+                sym_math = sym_math.subs(var, init)
 
         sym_math = smart_subs(sym_math, sbml_time_symbol, sp.Float(0))
 

python/tests/sbml_models/regression_2642.xml

@@ -0,0 +1,130 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
+  <model metaid="debug" id="debug" substanceUnits="substance" timeUnits="time_unit" volumeUnits="volume" areaUnits="area" lengthUnits="length">
+    <listOfUnitDefinitions>
+      <unitDefinition id="length">
+        <listOfUnits>
+          <unit kind="metre" exponent="1" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="area">
+        <listOfUnits>
+          <unit kind="metre" exponent="2" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="volume">
+        <listOfUnits>
+          <unit kind="litre" exponent="1" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="time_unit" name="time">
+        <listOfUnits>
+          <unit kind="second" exponent="1" scale="1" multiplier="6"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="substance">
+        <listOfUnits>
+          <unit kind="mole" exponent="1" scale="-9" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+    </listOfUnitDefinitions>
+    <listOfCompartments>
+      <compartment id="cytosol" spatialDimensions="3" size="0.0005" constant="true"/>
+    </listOfCompartments>
+    <listOfSpecies>
+      <species id="MR_i" compartment="cytosol" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="MR_m" compartment="cytosol" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+    </listOfSpecies>
+    <listOfParameters>
+      <parameter id="synthesis_MR" constant="true"/>
+      <parameter id="recycling" value="0.0000045" constant="true"/>
+      <parameter id="binding" constant="true"/>
+      <parameter id="MR_i_t0" value="1400000" constant="true"/>
+      <parameter id="MR_m_t0" value="600000" constant="true"/>
+    </listOfParameters>
+    <listOfInitialAssignments>
+      <initialAssignment symbol="MR_i">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <ci> MR_i_t0 </ci>
+        </math>
+      </initialAssignment>
+      <initialAssignment symbol="MR_m">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <ci> MR_m_t0 </ci>
+        </math>
+      </initialAssignment>
+      <initialAssignment symbol="synthesis_MR">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <divide/>
+            <cn> 0.00585177319092733 </cn>
+            <cn> 0.00606 </cn>
+          </apply>
+        </math>
+      </initialAssignment>
+      <initialAssignment symbol="binding">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <divide/>
+            <apply>
+              <minus/>
+              <apply>
+                <times/>
+                <ci> recycling </ci>
+                <ci> MR_i </ci>
+              </apply>
+              <ci> synthesis_MR </ci>
+            </apply>
+            <ci> MR_m </ci>
+          </apply>
+        </math>
+      </initialAssignment>
+    </listOfInitialAssignments>
+    <listOfReactions>
+      <reaction id="_J0" reversible="true">
+        <listOfProducts>
+          <speciesReference species="MR_i" stoichiometry="1" constant="true"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <ci> synthesis_MR </ci>
+          </math>
+        </kineticLaw>
+      </reaction>
+      <reaction id="_J1" reversible="true">
+        <listOfReactants>
+          <speciesReference species="MR_i" stoichiometry="1" constant="true"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference species="MR_m" stoichiometry="1" constant="true"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> recycling </ci>
+              <ci> MR_i </ci>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+      <reaction id="_J2" reversible="true">
+        <listOfReactants>
+          <speciesReference species="MR_m" stoichiometry="1" constant="true"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference species="MR_i" stoichiometry="1" constant="true"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> binding </ci>
+              <ci> MR_m </ci>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+    </listOfReactions>
+  </model>
+</sbml>

python/tests/test_sbml_import.py

@@ -852,3 +852,26 @@ def test_import_same_model_name():
 
             assert model_module_1c.get_model().getParameters()[0] == 1.0
             assert model_module_1c.get_model().module is model_module_1c
+
+
+@skip_on_valgrind
+def test_regression_2642():
+    sbml_file = Path(__file__).parent / "sbml_models" / "regression_2642.xml"
+    sbml_importer = amici.SbmlImporter(sbml_file)
+    model_name = "regression_2642"
+    with TemporaryDirectory(prefix="regression_2642") as outdir:
+        sbml_importer.sbml2amici(
+            model_name=model_name,
+            output_dir=outdir,
+        )
+        module = amici.import_model_module(
+            module_name=model_name, module_path=outdir
+        )
+        model = module.getModel()
+        solver = model.getSolver()
+        model.setTimepoints(np.linspace(0, 1, 3))
+        r = amici.runAmiciSimulation(model, solver)
+        assert (
+            len(np.unique(r.w[:, model.getExpressionIds().index("binding")]))
+            == 1
+        )

python/tests/valgrind-python.supp

@@ -1006,3 +1006,33 @@
    fun:_ZN4rpds*
    ...
 }
+
+{
+   Python
+   Memcheck:Leak
+   match-leak-kinds: definite
+   fun:realloc
+   fun:_PyUnicodeWriter_Finish
+   obj:/usr/bin/python3.*
+   ...
+}
+
+{
+   Something matplotlib. Cannot reproduce.
+   Memcheck:Leak
+   match-leak-kinds: definite
+   fun:_Znwm
+   fun:_ZL14PyFT2Font_initN8pybind116objectElSt8optionalISt6vectorIP9PyFT2FontSaIS4_EEEi.lto_priv.0
+   fun:_ZZN8pybind1112cpp_function10initializeIZNOS_6detail8initimpl7factoryIPFP9PyFT2FontNS_6objectElSt8optionalISt6vectorIS6_SaIS6_EEEiEPFNS2_9void_typeEvESD_SG_E7executeINS_6class_IS5_JEEEJNS_3argENS_5arg_vENS_7kw_onlyESN_SN_PKcEEEvRT_DpRKT0_EUlRNS2_16value_and_holderES7_lSC_iE_vJSY_S7_lSC_iEJNS_4nameENS_9is_methodENS_7siblingENS2_24is_new_style_constructorESM_SN_SO_SN_SN_SQ_EEEvOSR_PFT0_DpT1_EDpRKT2_ENUlRNS2_13function_callEE1_4_FUNES1F_
+   fun:_ZN8pybind1112cpp_function10dispatcherEP7_objectS2_S2_
+   fun:cfunction_call
+   fun:_PyObject_MakeTpCall
+   fun:_PyObject_VectorcallTstate
+   fun:_PyObject_VectorcallTstate
+   fun:method_vectorcall
+   fun:slot_tp_init
+   fun:type_call
+   fun:pybind11_meta_call
+   fun:_PyObject_MakeTpCall
+   fun:_PyEval_EvalFrameDefault
+}

scripts/run-valgrind-py.sh

@@ -23,7 +23,7 @@ if [ $# -eq 0 ]
     # for <class 'numpy.longdouble'> does not match any known type: falling back to type probe function.
 else
     # Run whatever was passed as arguments
-    command=($@)
+    command=("$@")
 fi
 
 

src/solver_cvodes.cpp

@@ -667,7 +667,7 @@
     int status
         = CVodeGetSensDky(solver_memory_.get(), t, k, sx_.getNVectorArray());
     if (status != CV_SUCCESS)
-        throw CvodeException(status, "CVodeGetSens");
+        throw CvodeException(status, "CVodeGetSensDky");
 }
 
 void CVodeSolver::getDkyB(realtype const t, int const k, int const which)
@@ -716,7 +716,7 @@
         status = CVodeAdjReInit(solver_memory_.get());
     } else {
         status = CVodeAdjInit(
-            solver_memory_.get(), static_cast<int>(maxsteps_),
+            solver_memory_.get(), static_cast<int>(maxsteps_ + 1),
             static_cast<int>(interp_type_)
         );
         setAdjInitDone();

version.txt

@@ -1 +1 @@
-0.30.0
+0.30.1