diff --git a/include/amici/defines.h b/include/amici/defines.h index a437abd7b3..eb7b01dd57 100644 --- a/include/amici/defines.h +++ b/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; diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index 71b1ba35d3..28bfe89489 100644 --- a/python/tests/test_preequilibration.py +++ b/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): diff --git a/src/amici.cpp b/src/amici.cpp index 76e49ca1e6..f41c4c44fe 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -59,6 +59,7 @@ std::map 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 run_simulation( diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index fd4239adb2..711ebd9e1a 100644 --- a/src/forwardproblem.cpp +++ b/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", diff --git a/src/solver.cpp b/src/solver.cpp index a53b3d4ed2..f5b6f7c37b 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -5,6 +5,7 @@ #include "amici/model.h" #include +#include #include #include @@ -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(); diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp index 291af6b786..39697533af 100644 --- a/src/solver_cvodes.cpp +++ 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);