Program Listing for File steadystateproblem.h

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

#ifndef AMICI_STEADYSTATE_PROBLEM_H
#define AMICI_STEADYSTATE_PROBLEM_H

#include <amici/defines.h>
#include <amici/model_state.h>
#include <amici/newton_solver.h>
#include <amici/vector.h>

namespace amici {

class ExpData;
class Solver;
class Model;
class BackwardProblem;
struct FwdSimWorkspace;
struct BwdSimWorkspace;

class WRMSComputer {
  public:
    WRMSComputer(
        int const n, SUNContext const sunctx, realtype const atol,
        realtype const rtol, AmiVector mask
    )
        : ewt_(n, sunctx)
        , rtol_(rtol)
        , atol_(atol)
        , mask_(mask) {}

    realtype wrms(AmiVector const& x, AmiVector const& x_ref);

  private:
    AmiVector ewt_;
    realtype rtol_;
    realtype atol_;
    AmiVector mask_;
};

class NewtonsMethod {
  public:
    NewtonsMethod(
        gsl::not_null<Model*> model, SUNContext sunctx,
        gsl::not_null<NewtonSolver*> solver,
        NewtonDampingFactorMode damping_factor_mode,
        realtype damping_factor_lower_bound, int max_steps, bool check_delta
    );

    void
    run(AmiVector& xdot, DEStateView const& state, WRMSComputer& wrms_computer);

    void compute_step(AmiVector const& xdot, DEStateView const& state);

    [[nodiscard]] AmiVector const& get_delta() const { return delta_; }

    [[nodiscard]] int get_num_steps() const { return i_step_; }

    [[nodiscard]] realtype get_wrms() const { return wrms_; }

  private:

    bool update_damping_factor(bool step_successful, double& gamma);

    realtype compute_wrms(
        AmiVector const& xdot, DEStateView const& state,
        WRMSComputer& wrms_computer
    );

    bool has_converged(
        AmiVector& xdot, DEStateView const& state, WRMSComputer& wrms_computer
    );

    static constexpr realtype conv_thresh = 1.0;

    gsl::not_null<Model*> model_;

    int max_steps_{0};

    NewtonDampingFactorMode damping_factor_mode_{NewtonDampingFactorMode::on};

    realtype damping_factor_lower_bound_{1e-8};

    bool check_delta_;

    gsl::not_null<NewtonSolver*> solver_;

    AmiVector delta_;

    AmiVector delta_old_;

    AmiVector x_old_;

    realtype wrms_ = INFINITY;

    int i_step_ = 0;
};

} // namespace amici
#endif // AMICI_STEADYSTATE_PROBLEM_H