From d121409b53bf037d760acf4b8e6b05a2c7f5b256 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 26 Sep 2024 13:29:40 +0200 Subject: [PATCH] ReturnData: remove protected data members (#2503) * ReturnData: remove protected data members * .. --- include/amici/model.h | 7 +- include/amici/model_state.h | 14 ++- include/amici/rdata.h | 117 +++++++++++++------------ src/model.cpp | 8 +- src/rdata.cpp | 170 ++++++++++++++++++++---------------- 5 files changed, 173 insertions(+), 143 deletions(-) diff --git a/include/amici/model.h b/include/amici/model.h index 1c020400ac..cfaca9526d 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -1427,19 +1427,20 @@ class Model : public AbstractModel, public ModelDimensions { * @param x_solver State variables with conservation laws applied * (solver returns this) */ - void fx_rdata(AmiVector& x_rdata, AmiVector const& x_solver); + void fx_rdata(gsl::span x_rdata, AmiVector const& x_solver); /** * @brief Expand conservation law for state sensitivities. * @param sx_rdata Output buffer for state variables sensitivities with - * conservation laws expanded (stored in `amici::ReturnData`). + * conservation laws expanded + * (stored in `amici::ReturnData` shape `nplist` x `nx`, row-major). * @param sx_solver State variables sensitivities with conservation laws * applied (solver returns this) * @param x_solver State variables with conservation laws * applied (solver returns this) */ void fsx_rdata( - AmiVectorArray& sx_rdata, AmiVectorArray const& sx_solver, + gsl::span sx_rdata, AmiVectorArray const& sx_solver, AmiVector const& x_solver ); diff --git a/include/amici/model_state.h b/include/amici/model_state.h index c6c46df3de..7c8d4f6b1e 100644 --- a/include/amici/model_state.h +++ b/include/amici/model_state.h @@ -332,11 +332,19 @@ struct ModelStateDerived { struct SimulationState { /** timepoint */ realtype t; - /** state variables */ + /** + * partial state vector, excluding states eliminated from conservation laws + */ AmiVector x; - /** state variables */ + /** + * partial time derivative of state vector, excluding states eliminated + * from conservation laws + */ AmiVector dx; - /** state variable sensitivity */ + /** + * partial sensitivity state vector array, excluding states eliminated from + * conservation laws + */ AmiVectorArray sx; /** state of the model that was used for simulation */ ModelState state; diff --git a/include/amici/rdata.h b/include/amici/rdata.h index b5b4c269b4..44f45ceef6 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`) evaluated at `t_last`. */ + /** time derivative (shape `nx_solver`) evaluated at `t_last`. */ std::vector xdot; /** - * Jacobian of differential equation right hand side (shape `nx` x `nx`, - * row-major) evaluated at `t_last`. + * Jacobian of differential equation right hand side (shape `nx_solver` x + * `nx_solver`, row-major) evaluated at `t_last`. */ std::vector J; @@ -156,11 +156,12 @@ class ReturnData : public ModelDimensions { */ std::vector s2rz; - /** state (shape `nt` x `nx`, row-major) */ + /** state (shape `nt` x `nx_rdata`, row-major) */ std::vector x; /** - * parameter derivative of state (shape `nt` x `nplist` x `nx`, row-major) + * parameter derivative of state (shape `nt` x `nplist` x `nx_rdata`, + * row-major) */ std::vector sx; @@ -358,18 +359,18 @@ class ReturnData : public ModelDimensions { */ realtype posteq_wrms = NAN; - /** initial state (shape `nx`) */ + /** initial state (shape `nx_rdata`) */ std::vector x0; - /** preequilibration steady state (shape `nx`) */ + /** preequilibration steady state (shape `nx_rdata`) */ std::vector x_ss; - /** initial sensitivities (shape `nplist` x `nx`, row-major) */ + /** initial sensitivities (shape `nplist` x `nx_rdata`, row-major) */ std::vector sx0; /** * preequilibration sensitivities - * (shape `nplist` x `nx`, row-major) + * (shape `nplist` x `nx_rdata`, row-major) */ std::vector sx_ss; @@ -463,28 +464,6 @@ class ReturnData : public ModelDimensions { /** offset for sigma_residuals */ realtype sigma_offset; - /** timepoint for model evaluation*/ - realtype t_; - - /** partial state vector, excluding states eliminated from conservation laws - */ - AmiVector x_solver_; - - /** partial time derivative of state vector, excluding states eliminated - * from conservation laws */ - AmiVector dx_solver_; - - /** partial sensitivity state vector array, excluding states eliminated from - * conservation laws */ - AmiVectorArray sx_solver_; - - /** full state vector, including states eliminated from conservation laws */ - AmiVector x_rdata_; - - /** full sensitivity state vector array, including states eliminated from - * conservation laws */ - AmiVectorArray sx_rdata_; - /** array of number of found roots for a certain event type * (shape `ne`) */ std::vector nroots_; @@ -568,18 +547,25 @@ class ReturnData : public ModelDimensions { template void storeJacobianAndDerivativeInReturnData(T const& problem, Model& model) { - readSimulationState(problem.getFinalSimulationState(), model); + auto simulation_state = problem.getFinalSimulationState(); + model.setModelState(simulation_state.state); AmiVector xdot(nx_solver); if (!this->xdot.empty() || !this->J.empty()) - model.fxdot(t_, x_solver_, dx_solver_, xdot); + model.fxdot( + simulation_state.t, simulation_state.x, simulation_state.dx, + xdot + ); if (!this->xdot.empty()) writeSlice(xdot, this->xdot); if (!this->J.empty()) { SUNMatrixWrapper J(nx_solver, nx_solver); - model.fJ(t_, 0.0, x_solver_, dx_solver_, xdot, J); + model.fJ( + simulation_state.t, 0.0, simulation_state.x, + simulation_state.dx, xdot, J + ); // CVODES uses colmajor, so we need to transform to rowmajor for (int ix = 0; ix < model.nx_solver; ix++) for (int jx = 0; jx < model.nx_solver; jx++) @@ -587,21 +573,18 @@ class ReturnData : public ModelDimensions { = J.data()[ix + model.nx_solver * jx]; } } - /** - * @brief sets member variables and model state according to provided - * simulation state - * @param state simulation state provided by Problem - * @param model model that was used for forward/backward simulation - */ - void readSimulationState(SimulationState const& state, Model& model); /** * @brief Residual function * @param it time index * @param model model that was used for forward/backward simulation + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance containing observable data */ - void fres(int it, Model& model, ExpData const& edata); + void fres( + int it, Model& model, SimulationState const& simulation_state, + ExpData const& edata + ); /** * @brief Chi-squared function @@ -614,17 +597,25 @@ class ReturnData : public ModelDimensions { * @brief Residual sensitivity function * @param it time index * @param model model that was used for forward/backward simulation + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance containing observable data */ - void fsres(int it, Model& model, ExpData const& edata); + void fsres( + int it, Model& model, SimulationState const& simulation_state, + ExpData const& edata + ); /** * @brief Fisher information matrix function * @param it time index * @param model model that was used for forward/backward simulation + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance containing observable data */ - void fFIM(int it, Model& model, ExpData const& edata); + void fFIM( + int it, Model& model, SimulationState const& simulation_state, + ExpData const& edata + ); /** * @brief Set likelihood, state variables, outputs and respective @@ -665,46 +656,58 @@ class ReturnData : public ModelDimensions { /** * @brief Extracts output information for data-points, expects that - * x_solver_ and sx_solver_ were set appropriately + * the model state was set appropriately * @param it timepoint index * @param model model that was used in forward solve + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance carrying experimental data */ - void getDataOutput(int it, Model& model, ExpData const* edata); + void getDataOutput( + int it, Model& model, SimulationState const& simulation_state, + ExpData const* edata + ); /** * @brief Extracts data information for forward sensitivity analysis, - * expects that x_solver_ and sx_solver_ were set appropriately + * expects that the model state was set appropriately * @param it index of current timepoint * @param model model that was used in forward solve + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance carrying experimental data */ - void getDataSensisFSA(int it, Model& model, ExpData const* edata); + void getDataSensisFSA( + int it, Model& model, SimulationState const& simulation_state, + ExpData const* edata + ); /** - * @brief Extracts output information for events, expects that x_solver_ - * and sx_solver_ were set appropriately + * @brief Extracts output information for events, expects that the model + * state was set appropriately * @param t event timepoint * @param rootidx information about which roots fired * (1 indicating fired, 0/-1 for not) * @param model model that was used in forward solve + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance carrying experimental data */ void getEventOutput( - realtype t, std::vector const rootidx, Model& model, - ExpData const* edata + realtype t, std::vector const& rootidx, Model& model, + SimulationState const& simulation_state, ExpData const* edata ); /** * @brief Extracts event information for forward sensitivity analysis, - * expects that x_solver_ and sx_solver_ were set appropriately + * expects the model state was set appropriately * @param ie index of event type * @param t event timepoint * @param model model that was used in forward solve + * @param simulation_state simulation state the timepoint `it` * @param edata ExpData instance carrying experimental data */ - void - getEventSensisFSA(int ie, realtype t, Model& model, ExpData const* edata); + void getEventSensisFSA( + int ie, realtype t, Model& model, + SimulationState const& simulation_state, ExpData const* edata + ); /** * @brief Updates contribution to likelihood from quadratures (xQB), @@ -725,12 +728,14 @@ class ReturnData : public ModelDimensions { * (llhS0), if no preequilibration was run or if forward sensitivities were * used * @param model model that was used for forward/backward simulation + * @param simulation_state simulation state the timepoint `it` * @param llhS0 contribution to likelihood for initial state sensitivities * @param xB vector with final adjoint state * (excluding conservation laws) */ void handleSx0Forward( - Model const& model, std::vector& llhS0, AmiVector& xB + Model const& model, SimulationState const& simulation_state, + std::vector& llhS0, AmiVector& xB ) const; }; diff --git a/src/model.cpp b/src/model.cpp index c27a5310d9..6719bedf9e 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1985,20 +1985,20 @@ void Model::fsx0_fixedParameters(AmiVectorArray& sx, AmiVector const& x) { void Model::fsdx0() {} -void Model::fx_rdata(AmiVector& x_rdata, AmiVector const& x) { +void Model::fx_rdata(gsl::span x_rdata, AmiVector const& x) { fx_rdata( x_rdata.data(), computeX_pos(x), state_.total_cl.data(), state_.unscaledParameters.data(), state_.fixedParameters.data() ); if (always_check_finite_) checkFinite( - x_rdata.getVector(), ModelQuantity::x_rdata, + x_rdata, ModelQuantity::x_rdata, std::numeric_limits::quiet_NaN() ); } void Model::fsx_rdata( - AmiVectorArray& sx_rdata, AmiVectorArray const& sx, + gsl::span sx_rdata, AmiVectorArray const& sx, AmiVector const& x_solver ) { realtype* stcl = nullptr; @@ -2006,7 +2006,7 @@ void Model::fsx_rdata( if (ncl() > 0) stcl = &state_.stotal_cl.at(plist(ip) * ncl()); fsx_rdata( - sx_rdata.data(ip), sx.data(ip), stcl, + &sx_rdata[ip * nx_rdata], sx.data(ip), stcl, state_.unscaledParameters.data(), state_.fixedParameters.data(), x_solver.data(), state_.total_cl.data(), plist(ip) ); diff --git a/src/rdata.cpp b/src/rdata.cpp index 649b420ccb..6d016298bd 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -46,10 +46,6 @@ ReturnData::ReturnData( , rdata_reporting(rdrm) , sigma_res(sigma_res) , sigma_offset(sigma_offset) - , x_solver_(nx_solver) - , sx_solver_(nx_solver, nplist) - , x_rdata_(nx) - , sx_rdata_(nx, nplist) , nroots_(ne) { switch (rdata_reporting) { case RDataReporting::full: @@ -194,16 +190,14 @@ void ReturnData::processSimulationObjects( void ReturnData::processPreEquilibration( SteadystateProblem const& preeq, Model& model ) { - readSimulationState(preeq.getFinalSimulationState(), model); + auto simulation_state = preeq.getFinalSimulationState(); + model.setModelState(simulation_state.state); if (!x_ss.empty()) { - model.fx_rdata(x_rdata_, x_solver_); - writeSlice(x_rdata_, x_ss); + model.fx_rdata(x_ss, simulation_state.x); } if (!sx_ss.empty() && sensi >= SensitivityOrder::first) { - model.fsx_rdata(sx_rdata_, sx_solver_, x_solver_); - for (int ip = 0; ip < nplist; ip++) - writeSlice(sx_rdata_[ip], slice(sx_ss, ip, nx)); + model.fsx_rdata(sx_ss, simulation_state.sx, simulation_state.x); } /* Get cpu time for pre-equilibration in milliseconds */ preeq_cpu_time = preeq.getCPUTime(); @@ -223,8 +217,9 @@ void ReturnData::processPostEquilibration( for (int it = 0; it < nt; it++) { auto t = model.getTimepoint(it); if (std::isinf(t)) { - readSimulationState(posteq.getFinalSimulationState(), model); - getDataOutput(it, model, edata); + auto simulation_state = posteq.getFinalSimulationState(); + model.setModelState(simulation_state.state); + getDataOutput(it, model, simulation_state, edata); } } /* Get cpu time for Newton solve in milliseconds */ @@ -248,25 +243,24 @@ void ReturnData::processForwardProblem( auto initialState = fwd.getInitialSimulationState(); if (initialState.x.getLength() == 0 && model.nx_solver > 0) return; // if x wasn't set forward problem failed during initialization - readSimulationState(initialState, model); + + model.setModelState(initialState.state); if (!x0.empty()) { - model.fx_rdata(x_rdata_, x_solver_); - writeSlice(x_rdata_, x0); + model.fx_rdata(x0, initialState.x); } if (!sx0.empty()) { - model.fsx_rdata(sx_rdata_, sx_solver_, x_solver_); - for (int ip = 0; ip < nplist; ip++) - writeSlice(sx_rdata_[ip], slice(sx0, ip, nx)); + model.fsx_rdata(sx0, initialState.sx, initialState.x); } // process timepoint data realtype tf = fwd.getFinalTime(); for (int it = 0; it < model.nt(); it++) { if (model.getTimepoint(it) <= tf) { - readSimulationState(fwd.getSimulationStateTimepoint(it), model); - getDataOutput(it, model, edata); + auto simulation_state = fwd.getSimulationStateTimepoint(it); + model.setModelState(simulation_state.state); + getDataOutput(it, model, simulation_state, edata); } else { // check for integration failure but consider postequilibration if (!std::isinf(model.getTimepoint(it))) @@ -278,38 +272,44 @@ void ReturnData::processForwardProblem( if (nz > 0) { auto rootidx = fwd.getRootIndexes(); for (int iroot = 0; iroot <= fwd.getEventCounter(); iroot++) { - readSimulationState(fwd.getSimulationStateEvent(iroot), model); - getEventOutput(t_, rootidx.at(iroot), model, edata); + auto simulation_state = fwd.getSimulationStateEvent(iroot); + model.setModelState(simulation_state.state); + getEventOutput( + simulation_state.t, rootidx.at(iroot), model, simulation_state, + edata + ); } } } -void ReturnData::getDataOutput(int it, Model& model, ExpData const* edata) { +void ReturnData::getDataOutput( + int it, Model& model, SimulationState const& simulation_state, + ExpData const* edata +) { if (!x.empty()) { - model.fx_rdata(x_rdata_, x_solver_); - writeSlice(x_rdata_, slice(x, it, nx)); + model.fx_rdata(slice(x, it, nx), simulation_state.x); } if (!w.empty()) - model.getExpression(slice(w, it, nw), ts[it], x_solver_); + model.getExpression(slice(w, it, nw), ts[it], simulation_state.x); if (!y.empty()) - model.getObservable(slice(y, it, ny), ts[it], x_solver_); + model.getObservable(slice(y, it, ny), ts[it], simulation_state.x); if (!sigmay.empty()) model.getObservableSigma(slice(sigmay, it, ny), it, edata); if (edata) { if (!isNaN(llh)) - model.addObservableObjective(llh, it, x_solver_, *edata); - fres(it, model, *edata); + model.addObservableObjective(llh, it, simulation_state.x, *edata); + fres(it, model, simulation_state, *edata); fchi2(it, *edata); } if (sensi >= SensitivityOrder::first && nplist > 0) { if (sensi_meth == SensitivityMethod::forward) { - getDataSensisFSA(it, model, edata); + getDataSensisFSA(it, model, simulation_state, edata); } else if (edata && !sllh.empty()) { model.addPartialObservableObjectiveSensitivity( - sllh, s2llh, it, x_solver_, *edata + sllh, s2llh, it, simulation_state.x, *edata ); } @@ -321,32 +321,37 @@ void ReturnData::getDataOutput(int it, Model& model, ExpData const* edata) { } } -void ReturnData::getDataSensisFSA(int it, Model& model, ExpData const* edata) { +void ReturnData::getDataSensisFSA( + int it, Model& model, SimulationState const& simulation_state, + ExpData const* edata +) { if (!sx.empty()) { - model.fsx_rdata(sx_rdata_, sx_solver_, x_solver_); - for (int ip = 0; ip < nplist; ip++) { - writeSlice(sx_rdata_[ip], slice(sx, it * nplist + ip, nx)); - } + model.fsx_rdata( + gsl::span(&sx.at(it * nplist * nx), nplist * nx), + simulation_state.sx, simulation_state.x + ); } if (!sy.empty()) { model.getObservableSensitivity( - slice(sy, it, nplist * ny), ts[it], x_solver_, sx_solver_ + slice(sy, it, nplist * ny), ts[it], simulation_state.x, + simulation_state.sx ); } if (edata) { if (!sllh.empty()) model.addObservableObjectiveSensitivity( - sllh, s2llh, it, x_solver_, sx_solver_, *edata + sllh, s2llh, it, simulation_state.x, simulation_state.sx, *edata ); - fsres(it, model, *edata); - fFIM(it, model, *edata); + fsres(it, model, simulation_state, *edata); + fFIM(it, model, simulation_state, *edata); } } void ReturnData::getEventOutput( - realtype t, std::vector rootidx, Model& model, ExpData const* edata + realtype t, std::vector const& rootidx, Model& model, + SimulationState const& simulation_state, ExpData const* edata ) { for (int ie = 0; ie < ne; ie++) { @@ -355,13 +360,15 @@ void ReturnData::getEventOutput( /* get event output */ if (!z.empty()) - model.getEvent(slice(z, nroots_.at(ie), nz), ie, t, x_solver_); + model.getEvent( + slice(z, nroots_.at(ie), nz), ie, t, simulation_state.x + ); /* if called from fillEvent at last timepoint, then also get the root function value */ if (t == model.getTimepoint(nt - 1)) if (!rz.empty()) model.getEventRegularization( - slice(rz, nroots_.at(ie), nz), ie, t, x_solver_ + slice(rz, nroots_.at(ie), nz), ie, t, simulation_state.x ); if (edata) { @@ -372,24 +379,25 @@ void ReturnData::getEventOutput( ); if (!isNaN(llh)) model.addEventObjective( - llh, ie, nroots_.at(ie), t, x_solver_, *edata + llh, ie, nroots_.at(ie), t, simulation_state.x, *edata ); /* if called from fillEvent at last timepoint, add regularization based on rz */ if (t == model.getTimepoint(nt - 1) && !isNaN(llh)) { model.addEventObjectiveRegularization( - llh, ie, nroots_.at(ie), t, x_solver_, *edata + llh, ie, nroots_.at(ie), t, simulation_state.x, *edata ); } } if (sensi >= SensitivityOrder::first) { if (sensi_meth == SensitivityMethod::forward) { - getEventSensisFSA(ie, t, model, edata); + getEventSensisFSA(ie, t, model, simulation_state, edata); } else if (edata && !sllh.empty()) { model.addPartialEventObjectiveSensitivity( - sllh, s2llh, ie, nroots_.at(ie), t, x_solver_, *edata + sllh, s2llh, ie, nroots_.at(ie), t, simulation_state.x, + *edata ); } } @@ -398,7 +406,8 @@ void ReturnData::getEventOutput( } void ReturnData::getEventSensisFSA( - int ie, realtype t, Model& model, ExpData const* edata + int ie, realtype t, Model& model, SimulationState const& simulation_state, + ExpData const* edata ) { if (t == model.getTimepoint(nt - 1)) { // call from fillEvent at last timepoint @@ -408,18 +417,20 @@ void ReturnData::getEventSensisFSA( ); if (!srz.empty()) model.getEventRegularizationSensitivity( - slice(srz, nroots_.at(ie), nz * nplist), ie, t, x_solver_, - sx_solver_ + slice(srz, nroots_.at(ie), nz * nplist), ie, t, + simulation_state.x, simulation_state.sx ); } else if (!sz.empty()) { model.getEventSensitivity( - slice(sz, nroots_.at(ie), nz * nplist), ie, t, x_solver_, sx_solver_ + slice(sz, nroots_.at(ie), nz * nplist), ie, t, simulation_state.x, + simulation_state.sx ); } if (edata && !sllh.empty()) { model.addEventObjectiveSensitivity( - sllh, s2llh, ie, nroots_.at(ie), t, x_solver_, sx_solver_, *edata + sllh, s2llh, ie, nroots_.at(ie), t, simulation_state.x, + simulation_state.sx, *edata ); } } @@ -430,7 +441,9 @@ void ReturnData::processBackwardProblem( ) { if (sllh.empty()) return; - readSimulationState(fwd.getInitialSimulationState(), model); + + auto simulation_state = fwd.getInitialSimulationState(); + model.setModelState(simulation_state.state); std::vector llhS0(model.nJ * model.nplist(), 0.0); auto xB = bwd.getAdjointState(); @@ -439,7 +452,7 @@ void ReturnData::processBackwardProblem( if (preeq && preeq->hasQuadrature()) { handleSx0Backward(model, *preeq, llhS0, xQB); } else { - handleSx0Forward(model, llhS0, xB); + handleSx0Forward(model, simulation_state, llhS0, xB); } for (int iJ = 0; iJ < model.nJ; iJ++) { @@ -486,7 +499,8 @@ void ReturnData::handleSx0Backward( } void ReturnData::handleSx0Forward( - Model const& model, std::vector& llhS0, AmiVector& xB + Model const& model, SimulationState const& simulation_state, + std::vector& llhS0, AmiVector& xB ) const { /* If preequilibration is run in forward mode or is not needed, then adjoint sensitivity analysis still needs the state sensitivities at t=0 (sx0), @@ -497,7 +511,7 @@ void ReturnData::handleSx0Forward( for (int ip = 0; ip < model.nplist(); ++ip) { llhS0[ip] = 0.0; for (int ix = 0; ix < model.nxtrue_solver; ++ix) { - llhS0[ip] += xB[ix] * sx_solver_.at(ix, ip); + llhS0[ip] += xB[ix] * simulation_state.sx.at(ix, ip); } } } else { @@ -506,9 +520,9 @@ void ReturnData::handleSx0Forward( for (int ix = 0; ix < model.nxtrue_solver; ++ix) { llhS0[ip + iJ * model.nplist()] += xB[ix + iJ * model.nxtrue_solver] - * sx_solver_.at(ix, ip) + * simulation_state.sx.at(ix, ip) + xB[ix] - * sx_solver_.at( + * simulation_state.sx.at( ix + iJ * model.nxtrue_solver, ip ); } @@ -575,17 +589,6 @@ void ReturnData::processSolver(Solver const& solver) { } } -void ReturnData::readSimulationState( - SimulationState const& state, Model& model -) { - x_solver_ = state.x; - dx_solver_ = state.dx; - if (computingFSA() || state.t == model.t0()) - sx_solver_ = state.sx; - t_ = state.t; - model.setModelState(state.state); -} - void ReturnData::invalidate(int const it_start) { if (it_start >= nt) return; @@ -817,12 +820,15 @@ static realtype fres_error(realtype sigma_y, realtype sigma_offset) { return sqrt(2 * log(sigma_y) + sigma_offset); } -void ReturnData::fres(int const it, Model& model, ExpData const& edata) { +void ReturnData::fres( + int const it, Model& model, SimulationState const& simulation_state, + ExpData const& edata +) { if (res.empty()) return; std::vector y_it(ny, 0.0); - model.getObservable(y_it, ts[it], x_solver_); + model.getObservable(y_it, ts[it], simulation_state.x); std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); @@ -877,14 +883,19 @@ fsres_error(realtype sigma_y, realtype ssigma_y, realtype sigma_offset) { return ssigma_y / (fres_error(sigma_y, sigma_offset) * sigma_y); } -void ReturnData::fsres(int const it, Model& model, ExpData const& edata) { +void ReturnData::fsres( + int const it, Model& model, SimulationState const& simulation_state, + ExpData const& edata +) { if (sres.empty()) return; std::vector y_it(ny, 0.0); - model.getObservable(y_it, ts[it], x_solver_); + model.getObservable(y_it, ts[it], simulation_state.x); std::vector sy_it(ny * nplist, 0.0); - model.getObservableSensitivity(sy_it, ts[it], x_solver_, sx_solver_); + model.getObservableSensitivity( + sy_it, ts[it], simulation_state.x, simulation_state.sx + ); std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata); @@ -917,14 +928,19 @@ void ReturnData::fsres(int const it, Model& model, ExpData const& edata) { } } -void ReturnData::fFIM(int it, Model& model, ExpData const& edata) { +void ReturnData::fFIM( + int it, Model& model, SimulationState const& simulation_state, + ExpData const& edata +) { if (FIM.empty()) return; std::vector y_it(ny, 0.0); - model.getObservable(y_it, ts[it], x_solver_); + model.getObservable(y_it, ts[it], simulation_state.x); std::vector sy_it(ny * nplist, 0.0); - model.getObservableSensitivity(sy_it, ts[it], x_solver_, sx_solver_); + model.getObservableSensitivity( + sy_it, ts[it], simulation_state.x, simulation_state.sx + ); std::vector sigmay_it(ny, 0.0); model.getObservableSigma(sigmay_it, it, &edata);