Improve debugging info in ReturnData

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,

include/amici/rdata.h

@@ -104,12 +104,12 @@ class ReturnData : public ModelDimensions {
      */
     std::vector<realtype> ts;
 
-    /** time derivative (shape `nx`) */
+    /** time derivative (shape `nx`) evaluated at `t_last`. */
     std::vector<realtype> xdot;
 
     /**
      * Jacobian of differential equation right hand side (shape `nx` x `nx`,
-     * row-major)
+     * row-major) evaluated at `t_last`.
      */
     std::vector<realtype> J;
 
@@ -456,6 +456,9 @@ class ReturnData : public ModelDimensions {
     /** log messages */
     std::vector<LogItem> messages;
 
+    /** The final internal time of the solver. */
+    realtype t_last{std::numeric_limits<realtype>::quiet_NaN()};
+
   protected:
     /** offset for sigma_residuals */
     realtype sigma_offset;

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]):

src/amici.cpp

@@ -236,6 +236,8 @@ std::unique_ptr<ReturnData> runAmiciSimulation(
         );
     }
 
+    rdata->t_last = solver.gett();
+
     rdata->cpu_time_total = cpu_timer.elapsed_milliseconds();
 
     // verify that reported CPU times are plausible

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_;

src/hdf5.cpp

@@ -532,6 +532,10 @@
         &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