Catch time overflow during steady-state simulation

* Remove the old and faulty check for finiteness of `x` in `Model::fxdot`. Closes #3076.
* The root cause of those NaNs in `x` was not CVODES.
  The problem is uncaught integer overflow during steady-state simulations.
  Catch that.
* Add a dedicated error code
* Add debug-level expectations to all integration functions of `Solver`

Author: dweindl

include/amici/defines.h

@@ -93,6 +93,7 @@ constexpr int AMICI_SINGULAR_JACOBIAN=     -9987;
 constexpr int AMICI_NOT_IMPLEMENTED=        -999;
 constexpr int AMICI_MAX_TIME_EXCEEDED  =   -1000;
 constexpr int AMICI_NOT_RUN=               -1001;
+constexpr int AMICI_T_OVERFLOW=            -1002;
 constexpr int AMICI_SUCCESS=                   0;
 constexpr int AMICI_DATA_RETURN=               1;
 constexpr int AMICI_ROOT_RETURN=               2;

python/tests/test_preequilibration.py

@@ -1,17 +1,18 @@
 """Tests for pre- and post-equilibration"""
 
+import inspect
 import itertools
 
 import amici
 import numpy as np
 import pytest
 from amici.debugging import get_model_for_preeq
+from amici.importers.antimony import antimony2amici
 from amici.sim.sundials import (
     AMICI_ERROR,
     AMICI_SUCCESS,
     ExpData,
     LogSeverity,
-    LogSeverity_debug,
     ParameterScaling,
     SensitivityMethod,
     SensitivityOrder,
@@ -592,7 +593,7 @@ def test_simulation_errors(preeq_fixture):
     for e in [edata, edata_preeq]:
         rdata = run_simulation(model, solver, e)
         assert rdata["status"] != AMICI_SUCCESS
-        assert rdata._swigptr.messages[0].severity == LogSeverity_debug
+        assert rdata._swigptr.messages[0].severity == LogSeverity.debug
         assert rdata._swigptr.messages[0].identifier == "EQUILIBRATION_FAILURE"
         assert (
             "exceeded maximum number of integration steps"
@@ -601,26 +602,51 @@ def test_simulation_errors(preeq_fixture):
         assert rdata._swigptr.messages[1].severity == LogSeverity.error
         assert rdata._swigptr.messages[1].identifier == "OTHER"
 
-    # too long simulations
-    solver.set_max_steps(int(1e4))
-    solver.set_relative_tolerance_steady_state(0.0)
-    solver.set_absolute_tolerance_steady_state(0.0)
-    # preeq & posteq
-    for e in [edata_preeq, edata]:
-        rdata = run_simulation(model, solver, e)
-        assert rdata["status"] != AMICI_SUCCESS
-        messages = []
-        # remove repeated RHSFUNC_FAIL messages
-        for message in rdata._swigptr.messages:
-            if not messages or message.message != messages[-1].message:
-                messages.append(message)
-        assert messages[0].severity == LogSeverity.debug
-        assert messages[0].identifier.endswith(":RHSFUNC_FAIL")
-        assert messages[1].severity == LogSeverity.debug
-        assert messages[1].identifier == "EQUILIBRATION_FAILURE"
-        assert "exceedingly long simulation time" in messages[1].message
-        assert messages[2].severity == LogSeverity.error
-        assert messages[2].identifier == "OTHER"
+
+def test_t_overflow():
+    """Test that equilibration fails with an informative error message
+    upon time overflow."""
+    module_name = inspect.stack()[0].function
+    ant_str = f"""
+    model {module_name}
+        # Constant increase so the solver will take large steps;
+        #  small enough to let `t` to overflow before `x`.
+        dxx_dt = 1e-16
+        xx' = dxx_dt
+        xx = 0
+    end
+    """
+    with TemporaryDirectory(prefix=module_name) as outdir:
+        model = antimony2amici(
+            ant_str,
+            model_name=module_name,
+            output_dir=outdir,
+            fixed_parameters=["dxx_dt"],
+        )
+        model.set_steady_state_computation_mode(
+            SteadyStateComputationMode.integrationOnly
+        )
+
+        # require simulation until forever
+        solver = model.create_solver()
+        solver.set_relative_tolerance_steady_state(0)
+        solver.set_absolute_tolerance_steady_state(0)
+
+        edata_preeq = ExpData(model)
+        edata_preeq.set_timepoints([np.inf])
+        edata_preeq.fixed_parameters_pre_equilibration = [1e-16]
+        edata_posteq = ExpData(model)
+        edata_posteq.set_timepoints([float("inf")])
+
+        for edata in (edata_preeq, edata_posteq):
+            rdata = run_simulation(model, solver, edata=edata)
+            assert rdata.status == AMICI_ERROR
+            for msg in rdata.messages:
+                if "exceedingly long simulation time" in msg.message:
+                    assert msg.identifier == "EQUILIBRATION_FAILURE"
+                    break
+            else:
+                assert False, list(rdata.messages)
 
 
 def test_get_model_for_preeq(preeq_fixture):

src/amici.cpp

@@ -59,6 +59,7 @@ std::map<int, std::string> simulation_status_to_str_map = {
     {AMICI_UNREC_SRHSFUNC_ERR, "AMICI_UNREC_SRHSFUNC_ERR"},
     {AMICI_RTFUNC_FAIL, "AMICI_RTFUNC_FAIL"},
     {AMICI_LINESEARCH_FAIL, "AMICI_LINESEARCH_FAIL"},
+    {AMICI_T_OVERFLOW, "AMICI_T_OVERFLOW"},
 };
 
 std::unique_ptr<ReturnData> run_simulation(

src/forwardproblem.cpp

@@ -197,6 +197,12 @@ void EventHandlingSimulator::run_steady_state(
         auto tout = std::isfinite(next_t_event)
                         ? next_t_event
                         : std::max(ws_->sol.t, 1.0) * 10;
+
+        if(!std::isfinite(tout)) {
+            // tout overflowed
+            throw IntegrationFailure(AMICI_T_OVERFLOW, tout);
+        }
+
         auto status = solver_->step(tout);
         solver_->write_solution(ws_->sol);
 
@@ -914,7 +920,7 @@ SteadyStateStatus SteadyStateProblem::find_steady_state_by_simulation(
                     ex.time
                 );
             return SteadyStateStatus::failed_convergence;
-        case AMICI_RHSFUNC_FAIL:
+        case amici::AMICI_T_OVERFLOW:
             if (model_->get_logger())
                 model_->get_logger()->log(
                     LogSeverity::debug, "EQUILIBRATION_FAILURE",

src/solver.cpp

@@ -5,6 +5,7 @@
 #include "amici/model.h"
 
 #include <sundials/sundials_context.h>
+#include <gsl/gsl-lite.hpp>
 
 #include <cstdio>
 #include <filesystem>
@@ -138,6 +139,8 @@ void Solver::apply_max_num_steps_B() const {
 }
 
 int Solver::run(realtype const tout) const {
+    gsl_ExpectsDebug(std::isfinite(tout));
+
     set_stop_time(tout);
     CpuTimer const cpu_timer;
     int status = AMICI_SUCCESS;
@@ -157,6 +160,8 @@ int Solver::run(realtype const tout) const {
 }
 
 int Solver::step(realtype const tout) const {
+    gsl_ExpectsDebug(std::isfinite(tout));
+
     int status = AMICI_SUCCESS;
 
     apply_max_num_steps();
@@ -173,6 +178,8 @@ int Solver::step(realtype const tout) const {
 }
 
 void Solver::run_b(realtype const tout) const {
+    gsl_ExpectsDebug(std::isfinite(tout));
+
     CpuTimer cpu_timer;
 
     apply_max_num_steps_B();

src/solver_cvodes.cpp

@@ -1178,14 +1178,6 @@ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void* user_data) {
         return AMICI_MAX_TIME_EXCEEDED;
     }
 
-    if (t > 1e200
-        && !model->check_finite(gsl::make_span(x), ModelQuantity::xdot, t)) {
-        /* when t is large (typically ~1e300), CVODES may pass all NaN x
-           to fxdot from which we typically cannot recover. To save time
-           on normal execution, we do not always want to check finiteness
-           of x, but only do so when t is large and we expect problems. */
-        return AMICI_UNRECOVERABLE_ERROR;
-    }
 
     model->fxdot(t, x, xdot);
     return model->check_finite(gsl::make_span(xdot), ModelQuantity::xdot, t);