Program Listing for File steadystateproblem.h

Return to documentation for file (include/amici/steadystateproblem.h)

#ifndef AMICI_STEADYSTATEPROBLEM_H
#define AMICI_STEADYSTATEPROBLEM_H

#include "amici/defines.h"
#include "amici/vector.h"
#include "amici/solver_cvodes.h"
#include "amici/forwardproblem.h"
#include <amici/newton_solver.h>

#include <nvector/nvector_serial.h>

#include <functional>
#include <memory>

namespace amici {

class ExpData;
class Solver;
class Model;

class SteadystateProblem {
  public:
    explicit SteadystateProblem(const Solver &solver,
                                const Model &model);

    void workSteadyStateProblem(Solver *solver, Model *model, int it);


    void workSteadyStateBackwardProblem(Solver *solver, Model *model,
                                        const BackwardProblem *bwd);

    void findSteadyState(Solver *solver,
                         NewtonSolver *newtonSolver,
                         Model *model, int it);

    void findSteadyStateByNewtonsMethod(NewtonSolver *newtonSolver,
                                        Model *model,
                                        bool newton_retry);

    void findSteadyStateBySimulation(const Solver *solver,
                                     Model *model,
                                     int it);

    void computeSteadyStateQuadrature(NewtonSolver *newtonSolver,
                                      const Solver *solver, Model *model);

    void getQuadratureByLinSolve(NewtonSolver *newtonSolver, Model *model);

    void getQuadratureBySimulation(const Solver *solver, Model *model);

    [[noreturn]] void handleSteadyStateFailure(const Solver *solver,
                                               Model *model);

    void writeErrorString(std::string *errorString, SteadyStateStatus
                          status) const;

    bool getSensitivityFlag(const Model *model, const Solver *solver, int it,
                            SteadyStateContext context);

    realtype getWrmsNorm(AmiVector const &x,
                         AmiVector const &xdot,
                         realtype atol,
                         realtype rtol,
                         AmiVector &ewt) const;

    bool checkConvergence(const Solver *solver, Model *model,
                          SensitivityMethod checkSensitivities);

    void applyNewtonsMethod(Model *model, NewtonSolver *newtonSolver,
                            bool newton_retry);

    void runSteadystateSimulation(const Solver *solver, Model *model, bool backward);

    std::unique_ptr<Solver> createSteadystateSimSolver(const Solver *solver,
                                                       Model *model,
                                                       bool forwardSensis,
                                                       bool backward) const;

    bool initializeBackwardProblem(Solver *solver, Model *model,
                                   const BackwardProblem *bwd);

    void computeQBfromQ(Model *model, const AmiVector &yQ, AmiVector &yQB) const;

    void storeSimulationState(Model *model, bool storesensi);

    const SimulationState &getFinalSimulationState() const {
        return state_;
    };

    const AmiVector &getEquilibrationQuadratures() const {
        return xQB_;
    }
    const AmiVector &getState() const {
        return x_;
    };


    const AmiVectorArray &getStateSensitivity() const {
        return sx_;
    };

    std::vector<realtype> const& getDJydx() const {
         return dJydx_;
     }

    double getCPUTime() const { return cpu_time_; }

    double getCPUTimeB() const { return cpu_timeB_; }

    std::vector<SteadyStateStatus> const& getSteadyStateStatus() const
    { return steady_state_status_; }

    realtype getSteadyStateTime() const { return t_; }

    realtype getResidualNorm() const { return wrms_; }

    const std::vector<int> &getNumSteps() const { return numsteps_; }

    int getNumStepsB() const { return numstepsB_; }

    const std::vector<int> &getNumLinSteps() const { return numlinsteps_; }

    void getAdjointUpdates(Model &model, const ExpData &edata);

    AmiVector const& getAdjointQuadrature() const { return xQB_; }

    bool hasQuadrature() const { return hasQuadrature_; }

    bool checkSteadyStateSuccess() const;

  private:
    realtype t_;
    AmiVector delta_;
    AmiVector ewt_;
    AmiVector ewtQB_;
    AmiVector rel_x_newton_;
    AmiVector x_newton_;
    AmiVector x_;
    AmiVector x_old_;
    AmiVector dx_;
    AmiVector xdot_;
    AmiVector xdot_old_;
    AmiVectorArray sx_;
    AmiVectorArray sdx_;
    AmiVector xB_;
    AmiVector xQ_;
    AmiVector xQB_;
    AmiVector xQBdot_;

    int max_steps_ {0};

    realtype wrms_ {NAN};

    std::vector<realtype> dJydx_;

    SimulationState state_;

    std::vector<int> numsteps_ {std::vector<int>(3, 0)};

    std::vector<int> numlinsteps_;

    int numstepsB_ {0};

    double cpu_time_ {0.0};

    double cpu_timeB_ {0.0};

    bool hasQuadrature_ {false};

    std::vector<SteadyStateStatus> steady_state_status_;
};

} // namespace amici
#endif // STEADYSTATEPROBLEM_H