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/forwardproblem.h"
#include "amici/solver_cvodes.h"
#include "amici/vector.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(const Solver &solver, Model &model, int it);

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

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

    const AmiVector &getEquilibrationQuadratures() const { return xQB_; }
    const AmiVector &getState() const { return state_.x; };

    const AmiVectorArray &getStateSensitivity() const { return state_.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 state_.t; }

    realtype getResidualNorm() const { return wrms_; }

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

    int getNumStepsB() const { return numstepsB_; }

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

    AmiVector const &getAdjointState() const { return xB_; }

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

    bool hasQuadrature() const { return hasQuadrature_; }

    bool checkSteadyStateSuccess() const;

  private:
    void findSteadyState(const Solver &solver, Model &model, int it);

    void findSteadyStateByNewtonsMethod(Model &model, bool newton_retry);

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

    void computeSteadyStateQuadrature(const Solver &solver, Model &model);

    void getQuadratureByLinSolve(Model &model);

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

    [[noreturn]] void handleSteadyStateFailure();

    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;

    realtype getWrms(Model &model, SensitivityMethod sensi_method);

    realtype getWrmsFSA(Model &model);

    void applyNewtonsMethod(Model &model, 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;

    void initializeForwardProblem(int it, const Solver &solver, Model &model);

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

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

    bool makePositiveAndCheckConvergence(Model &model);

    bool updateDampingFactor(bool step_successful);

    void flagUpdatedState();

    void updateSensiSimulation(const Solver &solver);

    void updateRightHandSide(Model &model);

    void getNewtonStep(Model &model);

    AmiVector delta_;
    AmiVector delta_old_;
    AmiVector ewt_;
    AmiVector ewtQB_;
    AmiVector x_old_;
    AmiVector xdot_;
    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)};

    int numstepsB_{0};

    double cpu_time_{0.0};

    double cpu_timeB_{0.0};

    bool hasQuadrature_{false};

    double gamma_{1.0};

    std::vector<SteadyStateStatus> steady_state_status_;

    realtype atol_{NAN};
    realtype rtol_{NAN};
    realtype atol_sensi_{NAN};
    realtype rtol_sensi_{NAN};
    realtype atol_quad_{NAN};
    realtype rtol_quad_{NAN};

    std::unique_ptr<NewtonSolver> newton_solver_{nullptr};

    NewtonDampingFactorMode damping_factor_mode_{NewtonDampingFactorMode::on};
    realtype damping_factor_lower_bound_{1e-8};
    bool newton_step_conv_ {false};
    bool check_sensi_conv_ {true};

    bool xdot_updated_ {false};
    bool delta_updated_ {false};
    bool sensis_updated_ {false};

};

} // namespace amici
#endif // STEADYSTATEPROBLEM_H