From 777deb0a1a3e50d39b6f6e1e172cf754514d62f3 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 5 Mar 2024 18:49:36 +0100 Subject: [PATCH] Always include timepoints in NaN/Inf warnings (#2347) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Always include timepoints in NaN/Inf warnings Closes #2328 --------- Co-authored-by: Fabian Fröhlich --- include/amici/model.h | 7 +++-- src/amici.cpp | 4 +-- src/model.cpp | 69 +++++++++++++++++++++++++------------------ src/model_dae.cpp | 3 +- src/model_ode.cpp | 3 +- src/solver_cvodes.cpp | 34 +++++++++++---------- src/solver_idas.cpp | 30 +++++++++++-------- 7 files changed, 89 insertions(+), 61 deletions(-) diff --git a/include/amici/model.h b/include/amici/model.h index 13de007e00..29b98aa913 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -1334,10 +1334,12 @@ class Model : public AbstractModel, public ModelDimensions { * * @param array * @param model_quantity The model quantity `array` corresponds to + * @param t Current timepoint * @return */ int checkFinite( - gsl::span array, ModelQuantity model_quantity + gsl::span array, ModelQuantity model_quantity, + realtype t ) const; /** * @brief Check if the given array has only finite elements. @@ -1347,11 +1349,12 @@ class Model : public AbstractModel, public ModelDimensions { * @param array Flattened matrix * @param model_quantity The model quantity `array` corresponds to * @param num_cols Number of columns of the non-flattened matrix + * @param t Current timepoint * @return */ int checkFinite( gsl::span array, ModelQuantity model_quantity, - size_t num_cols + size_t num_cols, realtype t ) const; /** diff --git a/src/amici.cpp b/src/amici.cpp index ef6d640e13..cc8c2f00ab 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -223,8 +223,8 @@ std::unique_ptr runAmiciSimulation( try { rdata->processSimulationObjects( - preeq.get(), fwd.get(), bwd_success ? bwd.get() : nullptr, posteq.get(), - model, solver, edata + preeq.get(), fwd.get(), bwd_success ? bwd.get() : nullptr, + posteq.get(), model, solver, edata ); } catch (std::exception const& ex) { rdata->status = AMICI_ERROR; diff --git a/src/model.cpp b/src/model.cpp index b671de3c73..cefdf1ac97 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -372,7 +372,7 @@ void Model::initializeStates(AmiVector& x) { std::copy(x0_solver.cbegin(), x0_solver.cend(), x.data()); } - checkFinite(x.getVector(), ModelQuantity::x0); + checkFinite(x.getVector(), ModelQuantity::x0, t0()); } void Model::initializeSplines() { @@ -1482,7 +1482,7 @@ void Model::addStateEventUpdate( ); if (always_check_finite_) { - checkFinite(derived_state_.deltax_, ModelQuantity::deltax); + checkFinite(derived_state_.deltax_, ModelQuantity::deltax, t); } // update @@ -1536,7 +1536,7 @@ void Model::addAdjointStateEventUpdate( ); if (always_check_finite_) { - checkFinite(derived_state_.deltaxB_, ModelQuantity::deltaxB); + checkFinite(derived_state_.deltaxB_, ModelQuantity::deltaxB, t); } // apply update @@ -1582,7 +1582,7 @@ void Model::updateHeavisideB(int const* rootsfound) { } int Model::checkFinite( - gsl::span array, ModelQuantity model_quantity + gsl::span array, ModelQuantity model_quantity, realtype t ) const { auto it = std::find_if(array.begin(), array.end(), [](realtype x) { return !std::isfinite(x); @@ -1654,23 +1654,27 @@ int Model::checkFinite( gsl_ExpectsDebug(false); model_quantity_str = std::to_string(static_cast(model_quantity)); } - if (logger) + if (logger) { + auto t_msg = std::isfinite(t) + ? std::string(" at t=" + std::to_string(t) + " ") + : std::string(); + logger->log( LogSeverity::warning, msg_id, - "AMICI encountered a %s value for %s[%i] (%s)", + "AMICI encountered a %s value for %s[%i] (%s)%s", non_finite_type.c_str(), model_quantity_str.c_str(), - gsl::narrow(flat_index), element_id.c_str() + gsl::narrow(flat_index), element_id.c_str(), t_msg.c_str() ); - + } // check upstream, without infinite recursion if (model_quantity != ModelQuantity::k && model_quantity != ModelQuantity::p && model_quantity != ModelQuantity::ts) { - checkFinite(state_.fixedParameters, ModelQuantity::k); - checkFinite(state_.unscaledParameters, ModelQuantity::p); - checkFinite(simulation_parameters_.ts_, ModelQuantity::ts); + checkFinite(state_.fixedParameters, ModelQuantity::k, t); + checkFinite(state_.unscaledParameters, ModelQuantity::p, t); + checkFinite(simulation_parameters_.ts_, ModelQuantity::ts, t); if (!always_check_finite_ && model_quantity != ModelQuantity::w) { // don't check twice if always_check_finite_ is true - checkFinite(derived_state_.w_, ModelQuantity::w); + checkFinite(derived_state_.w_, ModelQuantity::w, t); } } return AMICI_RECOVERABLE_ERROR; @@ -1678,7 +1682,7 @@ int Model::checkFinite( int Model::checkFinite( gsl::span array, ModelQuantity model_quantity, - size_t num_cols + size_t num_cols, realtype t ) const { auto it = std::find_if(array.begin(), array.end(), [](realtype x) { return !std::isfinite(x); @@ -1768,19 +1772,25 @@ int Model::checkFinite( model_quantity_str = std::to_string(static_cast(model_quantity)); } - if (logger) + if (logger) { + auto t_msg = std::isfinite(t) + ? std::string(" at t=" + std::to_string(t) + " ") + : std::string(); + logger->log( LogSeverity::warning, msg_id, - "AMICI encountered a %s value for %s[%i] (%s, %s)", + "AMICI encountered a %s value for %s[%i] (%s, %s)%s", non_finite_type.c_str(), model_quantity_str.c_str(), - gsl::narrow(flat_index), row_id.c_str(), col_id.c_str() + gsl::narrow(flat_index), row_id.c_str(), col_id.c_str(), + t_msg.c_str() ); + } // check upstream - checkFinite(state_.fixedParameters, ModelQuantity::k); - checkFinite(state_.unscaledParameters, ModelQuantity::p); - checkFinite(simulation_parameters_.ts_, ModelQuantity::ts); - checkFinite(derived_state_.w_, ModelQuantity::w); + checkFinite(state_.fixedParameters, ModelQuantity::k, t); + checkFinite(state_.unscaledParameters, ModelQuantity::p, t); + checkFinite(simulation_parameters_.ts_, ModelQuantity::ts, t); + checkFinite(derived_state_.w_, ModelQuantity::w, t); return AMICI_RECOVERABLE_ERROR; } @@ -1868,10 +1878,10 @@ int Model::checkFinite(SUNMatrix m, ModelQuantity model_quantity, realtype t) ); // check upstream - checkFinite(state_.fixedParameters, ModelQuantity::k); - checkFinite(state_.unscaledParameters, ModelQuantity::p); - checkFinite(simulation_parameters_.ts_, ModelQuantity::ts); - checkFinite(derived_state_.w_, ModelQuantity::w); + checkFinite(state_.fixedParameters, ModelQuantity::k, t); + checkFinite(state_.unscaledParameters, ModelQuantity::p, t); + checkFinite(simulation_parameters_.ts_, ModelQuantity::ts, t); + checkFinite(derived_state_.w_, ModelQuantity::w, t); return AMICI_RECOVERABLE_ERROR; } @@ -1895,7 +1905,7 @@ void Model::fx0(AmiVector& x) { state_.unscaledParameters.data(), state_.fixedParameters.data() ); - checkFinite(derived_state_.x_rdata_, ModelQuantity::x0_rdata); + checkFinite(derived_state_.x_rdata_, ModelQuantity::x0_rdata, t0()); } void Model::fx0_fixedParameters(AmiVector& x) { @@ -1982,7 +1992,10 @@ void Model::fx_rdata(AmiVector& x_rdata, AmiVector const& x) { state_.unscaledParameters.data(), state_.fixedParameters.data() ); if (always_check_finite_) - checkFinite(x_rdata.getVector(), ModelQuantity::x_rdata); + checkFinite( + x_rdata.getVector(), ModelQuantity::x_rdata, + std::numeric_limits::quiet_NaN() + ); } void Model::fsx_rdata( @@ -2072,7 +2085,7 @@ void Model::fy(realtype const t, AmiVector const& x) { if (always_check_finite_) { checkFinite( - gsl::make_span(derived_state_.y_.data(), ny), ModelQuantity::y + gsl::make_span(derived_state_.y_.data(), ny), ModelQuantity::y, t ); } } @@ -2875,7 +2888,7 @@ void Model::fw(realtype const t, realtype const* x, bool include_static) { state_.spl_.data(), include_static); if (always_check_finite_) { - checkFinite(derived_state_.w_, ModelQuantity::w); + checkFinite(derived_state_.w_, ModelQuantity::w, t); } } diff --git a/src/model_dae.cpp b/src/model_dae.cpp index ad397bc82f..22025e1181 100644 --- a/src/model_dae.cpp +++ b/src/model_dae.cpp @@ -138,7 +138,8 @@ void Model_DAE::fJDiag( fJSparse(t, 0.0, x.getNVector(), dx.getNVector(), derived_state_.J_); derived_state_.J_.refresh(); derived_state_.J_.to_diag(JDiag.getNVector()); - if (checkFinite(JDiag.getVector(), ModelQuantity::JDiag) != AMICI_SUCCESS) + if (checkFinite(JDiag.getVector(), ModelQuantity::JDiag, t) + != AMICI_SUCCESS) throw AmiException("Evaluation of fJDiag failed!"); } diff --git a/src/model_ode.cpp b/src/model_ode.cpp index e3e401e1bf..787df210ef 100644 --- a/src/model_ode.cpp +++ b/src/model_ode.cpp @@ -122,7 +122,8 @@ void Model_ODE::fJDiag( AmiVector const& x, AmiVector const& /*dx*/ ) { fJDiag(t, JDiag.getNVector(), x.getNVector()); - if (checkFinite(JDiag.getVector(), ModelQuantity::JDiag) != AMICI_SUCCESS) + if (checkFinite(JDiag.getVector(), ModelQuantity::JDiag, t) + != AMICI_SUCCESS) throw AmiException("Evaluation of fJDiag failed!"); } diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp index 4fe97720d8..0952724457 100644 --- a/src/solver_cvodes.cpp +++ b/src/solver_cvodes.cpp @@ -927,7 +927,7 @@ fJB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, SUNMatrix JB, Expects(model); model->fJB(t, x, xB, xBdot, JB); - return model->checkFinite(gsl::make_span(JB), ModelQuantity::JB); + return model->checkFinite(gsl::make_span(JB), ModelQuantity::JB, t); } /** @@ -978,7 +978,7 @@ static int fJSparseB( Expects(model); model->fJSparseB(t, x, xB, xBdot, JB); - return model->checkFinite(gsl::make_span(JB), ModelQuantity::JB); + return model->checkFinite(gsl::make_span(JB), ModelQuantity::JB, t); } /** @@ -1041,7 +1041,7 @@ fJv(N_Vector v, N_Vector Jv, realtype t, N_Vector x, N_Vector /*xdot*/, Expects(model); model->fJv(v, Jv, t, x); - return model->checkFinite(gsl::make_span(Jv), ModelQuantity::Jv); + return model->checkFinite(gsl::make_span(Jv), ModelQuantity::Jv, t); } /** @@ -1067,7 +1067,7 @@ static int fJvB( Expects(model); model->fJvB(vB, JvB, t, x, xB); - return model->checkFinite(gsl::make_span(JvB), ModelQuantity::JvB); + return model->checkFinite(gsl::make_span(JvB), ModelQuantity::JvB, t); } /** @@ -1094,7 +1094,7 @@ static int froot(realtype t, N_Vector x, realtype* root, void* user_data) { model->froot(t, x, gsl::make_span(root, model->ne_solver)); } return model->checkFinite( - gsl::make_span(root, model->ne_solver), ModelQuantity::root + gsl::make_span(root, model->ne_solver), ModelQuantity::root, t ); } @@ -1119,7 +1119,7 @@ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void* user_data) { } if (t > 1e200 - && !model->checkFinite(gsl::make_span(x), ModelQuantity::xdot)) { + && !model->checkFinite(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 @@ -1128,7 +1128,7 @@ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void* user_data) { } model->fxdot(t, x, xdot); - return model->checkFinite(gsl::make_span(xdot), ModelQuantity::xdot); + return model->checkFinite(gsl::make_span(xdot), ModelQuantity::xdot, t); } /** @@ -1154,7 +1154,7 @@ fxBdot(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot, void* user_data) { } model->fxBdot(t, x, xB, xBdot); - return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot); + return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot, t); } /** @@ -1174,7 +1174,7 @@ fqBdot(realtype t, N_Vector x, N_Vector xB, N_Vector qBdot, void* user_data) { Expects(model); model->fqBdot(t, x, xB, qBdot); - return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot); + return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot, t); } /** @@ -1193,7 +1193,9 @@ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector xBdot, void* user_data) { Expects(model); model->fxBdot_ss(t, xB, xBdot); - return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot_ss); + return model->checkFinite( + gsl::make_span(xBdot), ModelQuantity::xBdot_ss, t + ); } /** @@ -1212,7 +1214,9 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector qBdot, void* user_data) { Expects(model); model->fqBdot_ss(t, xB, qBdot); - return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot_ss); + return model->checkFinite( + gsl::make_span(qBdot), ModelQuantity::qBdot_ss, t + ); } /** @@ -1228,8 +1232,8 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector qBdot, void* user_data) { * @return status flag indicating successful execution */ static int fJSparseB_ss( - realtype /*t*/, N_Vector /*x*/, N_Vector xBdot, SUNMatrix JB, - void* user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/ + realtype t, N_Vector /*x*/, N_Vector xBdot, SUNMatrix JB, void* user_data, + N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/ ) { auto typed_udata = static_cast(user_data); Expects(typed_udata); @@ -1238,7 +1242,7 @@ static int fJSparseB_ss( model->fJSparseB_ss(JB); return model->checkFinite( - gsl::make_span(xBdot), ModelQuantity::JSparseB_ss + gsl::make_span(xBdot), ModelQuantity::JSparseB_ss, t ); } @@ -1267,7 +1271,7 @@ static int fsxdot( Expects(model); model->fsxdot(t, x, ip, sx, sxdot); - return model->checkFinite(gsl::make_span(sxdot), ModelQuantity::sxdot); + return model->checkFinite(gsl::make_span(sxdot), ModelQuantity::sxdot, t); } bool operator==(CVodeSolver const& a, CVodeSolver const& b) { diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp index 1c69a041db..fa08bc29f9 100644 --- a/src/solver_idas.cpp +++ b/src/solver_idas.cpp @@ -1065,7 +1065,7 @@ int fJv( Expects(model); model->fJv(t, x, dx, v, Jv, cj); - return model->checkFinite(gsl::make_span(Jv), ModelQuantity::Jv); + return model->checkFinite(gsl::make_span(Jv), ModelQuantity::Jv, t); } /** @@ -1096,7 +1096,7 @@ int fJvB( Expects(model); model->fJvB(t, x, dx, xB, dxB, vB, JvB, cj); - return model->checkFinite(gsl::make_span(JvB), ModelQuantity::JvB); + return model->checkFinite(gsl::make_span(JvB), ModelQuantity::JvB, t); } /** @@ -1118,7 +1118,7 @@ int froot( model->froot(t, x, dx, gsl::make_span(root, model->ne)); return model->checkFinite( - gsl::make_span(root, model->ne), ModelQuantity::root + gsl::make_span(root, model->ne), ModelQuantity::root, t ); } @@ -1144,7 +1144,7 @@ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot, void* user_data) { } if (t > 1e200 - && !model->checkFinite(gsl::make_span(x), ModelQuantity::xdot)) { + && !model->checkFinite(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 @@ -1153,7 +1153,7 @@ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot, void* user_data) { } model->fxdot(t, x, dx, xdot); - return model->checkFinite(gsl::make_span(xdot), ModelQuantity::xdot); + return model->checkFinite(gsl::make_span(xdot), ModelQuantity::xdot, t); } /** @@ -1183,7 +1183,7 @@ int fxBdot( } model->fxBdot(t, x, dx, xB, dxB, xBdot); - return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot); + return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot, t); } /** @@ -1208,7 +1208,7 @@ int fqBdot( Expects(model); model->fqBdot(t, x, dx, xB, dxB, qBdot); - return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot); + return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot, t); } /** @@ -1230,7 +1230,9 @@ static int fxBdot_ss( Expects(model); model->fxBdot_ss(t, xB, dxB, xBdot); - return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot_ss); + return model->checkFinite( + gsl::make_span(xBdot), ModelQuantity::xBdot_ss, t + ); } /** @@ -1252,7 +1254,9 @@ static int fqBdot_ss( Expects(model); model->fqBdot_ss(t, xB, dxB, qBdot); - return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot_ss); + return model->checkFinite( + gsl::make_span(qBdot), ModelQuantity::qBdot_ss, t + ); } /** @@ -1270,7 +1274,7 @@ static int fqBdot_ss( * @return status flag indicating successful execution */ static int fJSparseB_ss( - realtype /*t*/, realtype /*cj*/, N_Vector /*x*/, N_Vector /*dx*/, + realtype t, realtype /*cj*/, N_Vector /*x*/, N_Vector /*dx*/, N_Vector xBdot, SUNMatrix JB, void* user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/ ) { @@ -1281,7 +1285,7 @@ static int fJSparseB_ss( model->fJSparseB_ss(JB); return model->checkFinite( - gsl::make_span(xBdot), ModelQuantity::JSparseB_ss + gsl::make_span(xBdot), ModelQuantity::JSparseB_ss, t ); } @@ -1314,7 +1318,9 @@ int fsxdot( for (int ip = 0; ip < model->nplist(); ip++) { model->fsxdot(t, x, dx, ip, sx[ip], sdx[ip], sxdot[ip]); - if (model->checkFinite(gsl::make_span(sxdot[ip]), ModelQuantity::sxdot) + if (model->checkFinite( + gsl::make_span(sxdot[ip]), ModelQuantity::sxdot, t + ) != AMICI_SUCCESS) return AMICI_RECOVERABLE_ERROR; }