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/vector.h"
#include "amici/model.h"
#include "amici/misc.h"
#include "amici/forwardproblem.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, const Model &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 = 0;
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;
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.get());
// 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, const ExpData &edata);
void fchi2(int it, const ExpData &edata);
void fsres(int it, Model &model, const ExpData &edata);
void fFIM(int it, Model &model, const ExpData &edata);
void invalidate(int it_start);
void invalidateLLH();
void invalidateSLLH();
void applyChainRuleFactorToSimulationResults(const Model &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, const std::vector<int> rootidx,
Model &model, ExpData const *edata);
void getEventSensisFSA(int ie, realtype t, Model &model,
ExpData const *edata);
void handleSx0Backward(const Model &model, SteadystateProblem const &preeq,
std::vector<realtype> &llhS0, AmiVector &xQB) const;
void handleSx0Forward(const Model &model,
std::vector<realtype> &llhS0,
AmiVector &xB) const;
};
class ModelContext : public ContextManager {
public:
explicit ModelContext(Model *model);
ModelContext &operator=(const ModelContext &other) = delete;
~ModelContext();
void restore();
private:
Model *model_ {nullptr};
ModelState original_state_;
};
} // namespace amici
#endif /* _MY_RDATA */