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