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/model_state.h>
#include <amici/newton_solver.h>
#include <amici/vector.h>

#include <nvector/nvector_serial.h>

#include <memory>

namespace amici {

class ExpData;
class Solver;
class Model;

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

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

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

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

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

    AmiVectorArray const& 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_; }

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

    int getNumStepsB() const { return numstepsB_; }

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

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

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

    bool hasQuadrature() const { return hasQuadrature_; }

    bool checkSteadyStateSuccess() const;

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

    void findSteadyStateByNewtonsMethod(Model& model, bool newton_retry);

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

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

    void getQuadratureByLinSolve(Model& model);

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

    [[noreturn]] void handleSteadyStateFailure(
        bool tried_newton_1, bool tried_simulation, bool tried_newton_2
    );

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

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

    realtype getWrmsNorm(
        AmiVector const& x, AmiVector const& xdot, AmiVector const& mask,
        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(Solver const& solver, Model& model, bool backward);

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

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

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

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

    bool makePositiveAndCheckConvergence(Model& model);

    bool updateDampingFactor(bool step_successful);

    void flagUpdatedState();

    void updateSensiSimulation(Solver const& 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