From b9965c478d082dcb05841c37c610bbe83ee3aa23 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 5 Mar 2024 08:39:31 +0100 Subject: [PATCH 1/2] Improve debugging info in ReturnData * Add `ReturnData::t_last` that holds the last solver timepoint (Closes #2246) * `ReturnData::J` is now the evaluted at `t_last` instead of the last successfully reached output timepoint (Closes #2334) * `ReturnData::xdot` is now the evaluted at `t_last` instead of the last successfully reached output timepoint --- include/amici/forwardproblem.h | 2 +- include/amici/rdata.h | 7 +++++-- python/sdist/amici/numpy.py | 1 + src/amici.cpp | 2 ++ src/forwardproblem.cpp | 5 ++++- src/hdf5.cpp | 4 ++++ 6 files changed, 17 insertions(+), 4 deletions(-) diff --git a/include/amici/forwardproblem.h b/include/amici/forwardproblem.h index 8c500bf73c..67658ebee4 100644 --- a/include/amici/forwardproblem.h +++ b/include/amici/forwardproblem.h @@ -316,7 +316,7 @@ class ForwardProblem { * @brief Creates a carbon copy of the current simulation state variables * @return state */ - SimulationState getSimulationState() const; + SimulationState getSimulationState(); /** array of index vectors (dimension ne) indicating whether the respective * root has been detected for all so far encountered discontinuities, diff --git a/include/amici/rdata.h b/include/amici/rdata.h index 6357d24748..b5b4c269b4 100644 --- a/include/amici/rdata.h +++ b/include/amici/rdata.h @@ -104,12 +104,12 @@ class ReturnData : public ModelDimensions { */ std::vector ts; - /** time derivative (shape `nx`) */ + /** time derivative (shape `nx`) evaluated at `t_last`. */ std::vector xdot; /** * Jacobian of differential equation right hand side (shape `nx` x `nx`, - * row-major) + * row-major) evaluated at `t_last`. */ std::vector J; @@ -456,6 +456,9 @@ class ReturnData : public ModelDimensions { /** log messages */ std::vector messages; + /** The final internal time of the solver. */ + realtype t_last{std::numeric_limits::quiet_NaN()}; + protected: /** offset for sigma_residuals */ realtype sigma_offset; diff --git a/python/sdist/amici/numpy.py b/python/sdist/amici/numpy.py index f75d927b7b..c1aef949c6 100644 --- a/python/sdist/amici/numpy.py +++ b/python/sdist/amici/numpy.py @@ -239,6 +239,7 @@ class ReturnDataView(SwigPtrView): "cpu_timeB", "cpu_time_total", "messages", + "t_last", ] def __init__(self, rdata: Union[ReturnDataPtr, ReturnData]): diff --git a/src/amici.cpp b/src/amici.cpp index cc8c2f00ab..6166ab9352 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -236,6 +236,8 @@ std::unique_ptr runAmiciSimulation( ); } + rdata->t_last = solver.gett(); + rdata->cpu_time_total = cpu_timer.elapsed_milliseconds(); // verify that reported CPU times are plausible diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index b4bcd9740f..7fb66d7340 100644 --- a/src/forwardproblem.cpp +++ b/src/forwardproblem.cpp @@ -416,7 +416,10 @@ void ForwardProblem::getAdjointUpdates(Model& model, ExpData const& edata) { } } -SimulationState ForwardProblem::getSimulationState() const { +SimulationState ForwardProblem::getSimulationState() { + if(std::isfinite(solver->gett())) { + solver->writeSolution(&t_, x_, dx_, sx_, dx_); + } auto state = SimulationState(); state.t = t_; state.x = x_; diff --git a/src/hdf5.cpp b/src/hdf5.cpp index 4efed799a8..0dffb533bd 100644 --- a/src/hdf5.cpp +++ b/src/hdf5.cpp @@ -532,6 +532,10 @@ void writeReturnDataDiagnosis( &rdata.cpu_time_total, 1 ); + H5LTset_attribute_double( + file.getId(), hdf5Location.c_str(), "t_last", &rdata.t_last, 1 + ); + if (!rdata.J.empty()) createAndWriteDouble2DDataset( file, hdf5Location + "/J", rdata.J, rdata.nx, rdata.nx From 68d18e12fc9c4470fb9da8c1875132e06931d9c1 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 6 Mar 2024 23:09:57 +0100 Subject: [PATCH 2/2] skip xdot/J checks in case of failures --- src/forwardproblem.cpp | 2 +- tests/cpp/testfunctions.cpp | 23 +++++++++++++---------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index 7fb66d7340..ef4585f9e9 100644 --- a/src/forwardproblem.cpp +++ b/src/forwardproblem.cpp @@ -417,7 +417,7 @@ void ForwardProblem::getAdjointUpdates(Model& model, ExpData const& edata) { } SimulationState ForwardProblem::getSimulationState() { - if(std::isfinite(solver->gett())) { + if (std::isfinite(solver->gett())) { solver->writeSolution(&t_, x_, dx_, sx_, dx_); } auto state = SimulationState(); diff --git a/tests/cpp/testfunctions.cpp b/tests/cpp/testfunctions.cpp index 3e03eedfc4..8a8054f147 100644 --- a/tests/cpp/testfunctions.cpp +++ b/tests/cpp/testfunctions.cpp @@ -192,13 +192,6 @@ void verifyReturnData(std::string const& hdffile, std::string const& resultPath, // CHECK_EQUAL(AMICI_O2MODE_FULL, udata->o2mode); - if(hdf5::locationExists(file, resultPath + "/diagnosis/J")) { - expected = hdf5::getDoubleDataset2D(file, resultPath + "/diagnosis/J", m, n); - checkEqualArray(expected, rdata->J, atol, rtol, "J"); - } else { - ASSERT_TRUE(rdata->J.empty()); - } - if(hdf5::locationExists(file, resultPath + "/y")) { expected = hdf5::getDoubleDataset2D(file, resultPath + "/y", m, n); checkEqualArray(expected, rdata->y, atol, rtol, "y"); @@ -234,12 +227,22 @@ void verifyReturnData(std::string const& hdffile, std::string const& resultPath, ASSERT_TRUE(rdata->sigmaz.empty()); } - expected = hdf5::getDoubleDataset1D(file, resultPath + "/diagnosis/xdot"); - checkEqualArray(expected, rdata->xdot, atol, rtol, "xdot"); - expected = hdf5::getDoubleDataset1D(file, resultPath + "/x0"); checkEqualArray(expected, rdata->x0, atol, rtol, "x0"); + if(rdata->status == AMICI_SUCCESS) { + // for the failed cases, the stored results don't match + // since https://github.com/AMICI-dev/AMICI/pull/2349 + expected = hdf5::getDoubleDataset1D(file, resultPath + "/diagnosis/xdot"); + checkEqualArray(expected, rdata->xdot, atol, rtol, "xdot"); + + if(hdf5::locationExists(file, resultPath + "/diagnosis/J")) { + expected = hdf5::getDoubleDataset2D(file, resultPath + "/diagnosis/J", m, n); + checkEqualArray(expected, rdata->J, atol, rtol, "J"); + } else { + ASSERT_TRUE(rdata->J.empty()); + } + } if(rdata->sensi >= SensitivityOrder::first) { verifyReturnDataSensitivities(file, resultPath, rdata, model, atol, rtol); } else {