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

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

    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;

    std::vector<int> preeq_numlinsteps;

    int preeq_numstepsB = 0;

    std::vector<int> posteq_numsteps;

    std::vector<int> posteq_numlinsteps;

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

  protected:

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

    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 */