Program Listing for File rdata.h
↰ Return to documentation for file (include/amici/rdata.h
)
#ifndef AMICI_RDATA_H
#define AMICI_RDATA_H
#include "amici/defines.h"
#include "amici/logging.h"
#include "amici/misc.h"
#include "amici/model.h"
#include "amici/vector.h"
#include <vector>
namespace amici {
class ReturnData;
class Solver;
class ExpData;
class ForwardProblem;
class BackwardProblem;
class SteadystateProblem;
} // namespace amici
namespace boost {
namespace serialization {
template <class Archive>
void serialize(Archive& ar, amici::ReturnData& r, unsigned int version);
}
} // namespace boost
namespace amici {
class ReturnData : public ModelDimensions {
public:
ReturnData() = default;
ReturnData(
std::vector<realtype> ts, ModelDimensions const& model_dimensions,
int nplist, int nmaxevent, int nt, int newton_maxsteps,
std::vector<ParameterScaling> pscale, SecondOrderMode o2mode,
SensitivityOrder sensi, SensitivityMethod sensi_meth,
RDataReporting rdrm, bool quadratic_llh, bool sigma_res,
realtype sigma_offset
);
ReturnData(Solver const& solver, Model const& model);
~ReturnData() = default;
void processSimulationObjects(
SteadystateProblem const* preeq, ForwardProblem const* fwd,
BackwardProblem const* bwd, SteadystateProblem const* posteq,
Model& model, Solver const& solver, ExpData const* edata
);
std::string id;
std::vector<realtype> ts;
std::vector<realtype> xdot;
std::vector<realtype> J;
std::vector<realtype> w;
std::vector<realtype> z;
std::vector<realtype> sigmaz;
std::vector<realtype> sz;
std::vector<realtype> ssigmaz;
std::vector<realtype> rz;
std::vector<realtype> srz;
std::vector<realtype> s2rz;
std::vector<realtype> x;
std::vector<realtype> sx;
std::vector<realtype> y;
std::vector<realtype> sigmay;
std::vector<realtype> sy;
std::vector<realtype> ssigmay;
std::vector<realtype> res;
std::vector<realtype> sres;
std::vector<realtype> FIM;
std::vector<int> numsteps;
std::vector<int> numstepsB;
std::vector<int> numrhsevals;
std::vector<int> numrhsevalsB;
std::vector<int> numerrtestfails;
std::vector<int> numerrtestfailsB;
std::vector<int> numnonlinsolvconvfails;
std::vector<int> numnonlinsolvconvfailsB;
std::vector<int> order;
double cpu_time = 0.0;
double cpu_timeB = 0.0;
double cpu_time_total = 0.0;
std::vector<SteadyStateStatus> preeq_status;
double preeq_cpu_time = 0.0;
double preeq_cpu_timeB = 0.0;
std::vector<SteadyStateStatus> posteq_status;
double posteq_cpu_time = 0.0;
double posteq_cpu_timeB = 0.0;
std::vector<int> preeq_numsteps;
int preeq_numstepsB = 0;
std::vector<int> posteq_numsteps;
int posteq_numstepsB = 0;
realtype preeq_t = NAN;
realtype preeq_wrms = NAN;
realtype posteq_t = NAN;
realtype posteq_wrms = NAN;
std::vector<realtype> x0;
std::vector<realtype> x_ss;
std::vector<realtype> sx0;
std::vector<realtype> sx_ss;
realtype llh = 0.0;
realtype chi2 = 0.0;
std::vector<realtype> sllh;
std::vector<realtype> s2llh;
int status = AMICI_NOT_RUN;
int nx{0};
int nxtrue{0};
int nplist{0};
int nmaxevent{0};
int nt{0};
int newton_maxsteps{0};
std::vector<ParameterScaling> pscale;
SecondOrderMode o2mode{SecondOrderMode::none};
SensitivityOrder sensi{SensitivityOrder::none};
SensitivityMethod sensi_meth{SensitivityMethod::none};
RDataReporting rdata_reporting{RDataReporting::full};
template <class Archive>
friend void boost::serialization::serialize(
Archive& ar, ReturnData& r, unsigned int version
);
bool sigma_res;
std::vector<LogItem> messages;
realtype t_last{std::numeric_limits<realtype>::quiet_NaN()};
protected:
realtype sigma_offset;
realtype t_;
AmiVector x_solver_;
AmiVector dx_solver_;
AmiVectorArray sx_solver_;
AmiVector x_rdata_;
AmiVectorArray sx_rdata_;
std::vector<int> nroots_;
void initializeLikelihoodReporting(bool quadratic_llh);
void initializeResidualReporting(bool enable_res);
void initializeFullReporting(bool enable_fim);
void initializeObjectiveFunction(bool enable_chi2);
void processPreEquilibration(SteadystateProblem const& preeq, Model& model);
void processPostEquilibration(
SteadystateProblem const& posteq, Model& model, ExpData const* edata
);
void processForwardProblem(
ForwardProblem const& fwd, Model& model, ExpData const* edata
);
void processBackwardProblem(
ForwardProblem const& fwd, BackwardProblem const& bwd,
SteadystateProblem const* preeq, Model& model
);
void processSolver(Solver const& solver);
template <class T>
void
storeJacobianAndDerivativeInReturnData(T const& problem, Model& model) {
readSimulationState(problem.getFinalSimulationState(), model);
AmiVector xdot(nx_solver);
if (!this->xdot.empty() || !this->J.empty())
model.fxdot(t_, x_solver_, dx_solver_, 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);
// 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++)
this->J.at(ix * model.nx_solver + jx)
= J.data()[ix + model.nx_solver * jx];
}
}
void readSimulationState(SimulationState const& state, Model& model);
void fres(int it, Model& model, ExpData const& edata);
void fchi2(int it, ExpData const& edata);
void fsres(int it, Model& model, ExpData const& edata);
void fFIM(int it, Model& model, ExpData const& edata);
void invalidate(int it_start);
void invalidateLLH();
void invalidateSLLH();
void applyChainRuleFactorToSimulationResults(Model const& model);
bool computingFSA() const {
return (
sensi_meth == SensitivityMethod::forward
&& sensi >= SensitivityOrder::first
);
}
void getDataOutput(int it, Model& model, ExpData const* edata);
void getDataSensisFSA(int it, Model& model, ExpData const* edata);
void getEventOutput(
realtype t, std::vector<int> const rootidx, Model& model,
ExpData const* edata
);
void
getEventSensisFSA(int ie, realtype t, Model& model, ExpData const* edata);
void handleSx0Backward(
Model const& model, SteadystateProblem const& preeq,
std::vector<realtype>& llhS0, AmiVector& xQB
) const;
void handleSx0Forward(
Model const& model, std::vector<realtype>& llhS0, AmiVector& xB
) const;
};
class ModelContext : public ContextManager {
public:
explicit ModelContext(Model* model);
ModelContext& operator=(ModelContext const& other) = delete;
~ModelContext();
void restore();
private:
Model* model_{nullptr};
ModelState original_state_;
};
} // namespace amici
#endif /* _MY_RDATA */