Program Listing for File solver.h¶
↰ Return to documentation for file (include/amici/solver.h
)
#ifndef AMICI_SOLVER_H
#define AMICI_SOLVER_H
#include "amici/amici.h"
#include "amici/defines.h"
#include "amici/sundials_linsol_wrapper.h"
#include "amici/symbolic_functions.h"
#include "amici/vector.h"
#include <cmath>
#include <functional>
#include <memory>
namespace amici {
class ReturnData;
class ForwardProblem;
class BackwardProblem;
class Model;
class Solver;
class AmiciApplication;
extern AmiciApplication defaultContext;
} // namespace amici
// for serialization friend in Solver
namespace boost {
namespace serialization {
template <class Archive>
void serialize(Archive &ar, amici::Solver &s, unsigned int version);
}
} // namespace boost::serialization
namespace amici {
class Solver {
public:
Solver() = default;
Solver(AmiciApplication *app);
Solver(const Solver &other);
virtual ~Solver() = default;
virtual Solver *clone() const = 0;
int run(realtype tout) const;
int step(realtype tout) const;
void runB(realtype tout) const;
void setup(realtype t0, Model *model, const AmiVector &x0,
const AmiVector &dx0, const AmiVectorArray &sx0,
const AmiVectorArray &sdx0) const;
void setupB(int *which, realtype tf, Model *model, const AmiVector &xB0,
const AmiVector &dxB0, const AmiVector &xQB0) const;
void setupSteadystate(const realtype t0, Model *model, const AmiVector &x0,
const AmiVector &dx0, const AmiVector &xB0,
const AmiVector &dxB0, const AmiVector &xQ0) const;
void updateAndReinitStatesAndSensitivities(Model *model);
virtual void getRootInfo(int *rootsfound) const = 0;
virtual void calcIC(realtype tout1) const = 0;
virtual void calcICB(int which, realtype tout1) const = 0;
virtual void solveB(realtype tBout, int itaskB) const = 0;
virtual void turnOffRootFinding() const = 0;
SensitivityMethod getSensitivityMethod() const;
void setSensitivityMethod(SensitivityMethod sensi_meth);
SensitivityMethod getSensitivityMethodPreequilibration() const;
void setSensitivityMethodPreequilibration(SensitivityMethod sensi_meth_preeq);
void switchForwardSensisOff() const;
int getNewtonMaxSteps() const;
void setNewtonMaxSteps(int newton_maxsteps);
bool getPreequilibration() const;
void setPreequilibration(bool require_preequilibration);
int getNewtonMaxLinearSteps() const;
void setNewtonMaxLinearSteps(int newton_maxlinsteps);
NewtonDampingFactorMode getNewtonDampingFactorMode() const;
void setNewtonDampingFactorMode(NewtonDampingFactorMode dampingFactorMode);
double getNewtonDampingFactorLowerBound() const;
void setNewtonDampingFactorLowerBound(double dampingFactorLowerBound);
SensitivityOrder getSensitivityOrder() const;
void setSensitivityOrder(SensitivityOrder sensi);
double getRelativeTolerance() const;
void setRelativeTolerance(double rtol);
double getAbsoluteTolerance() const;
void setAbsoluteTolerance(double atol);
double getRelativeToleranceFSA() const;
void setRelativeToleranceFSA(double rtol);
double getAbsoluteToleranceFSA() const;
void setAbsoluteToleranceFSA(double atol);
double getRelativeToleranceB() const;
void setRelativeToleranceB(double rtol);
double getAbsoluteToleranceB() const;
void setAbsoluteToleranceB(double atol);
double getRelativeToleranceQuadratures() const;
void setRelativeToleranceQuadratures(double rtol);
double getAbsoluteToleranceQuadratures() const;
void setAbsoluteToleranceQuadratures(double atol);
double getRelativeToleranceSteadyState() const;
void setRelativeToleranceSteadyState(double rtol);
double getAbsoluteToleranceSteadyState() const;
void setAbsoluteToleranceSteadyState(double atol);
double getRelativeToleranceSteadyStateSensi() const;
void setRelativeToleranceSteadyStateSensi(double rtol);
double getAbsoluteToleranceSteadyStateSensi() const;
void setAbsoluteToleranceSteadyStateSensi(double atol);
long int getMaxSteps() const;
void setMaxSteps(long int maxsteps);
long int getMaxStepsBackwardProblem() const;
void setMaxStepsBackwardProblem(long int maxsteps);
LinearMultistepMethod getLinearMultistepMethod() const;
void setLinearMultistepMethod(LinearMultistepMethod lmm);
NonlinearSolverIteration getNonlinearSolverIteration() const;
void setNonlinearSolverIteration(NonlinearSolverIteration iter);
InterpolationType getInterpolationType() const;
void setInterpolationType(InterpolationType interpType);
int getStateOrdering() const;
void setStateOrdering(int ordering);
bool getStabilityLimitFlag() const;
void setStabilityLimitFlag(bool stldet);
LinearSolver getLinearSolver() const;
void setLinearSolver(LinearSolver linsol);
InternalSensitivityMethod getInternalSensitivityMethod() const;
void setInternalSensitivityMethod(InternalSensitivityMethod ism);
RDataReporting getReturnDataReportingMode() const;
void setReturnDataReportingMode(RDataReporting rdrm);
void writeSolution(realtype *t, AmiVector &x, AmiVector &dx,
AmiVectorArray &sx, AmiVector &xQ) const;
void writeSolutionB(realtype *t, AmiVector &xB, AmiVector &dxB,
AmiVector &xQB, int which) const;
const AmiVector &getState(realtype t) const;
const AmiVector &getDerivativeState(realtype t) const;
const AmiVectorArray &getStateSensitivity(realtype t) const;
const AmiVector &getAdjointState(int which, realtype t) const;
const AmiVector &getAdjointDerivativeState(int which, realtype t) const;
const AmiVector &getAdjointQuadrature(int which, realtype t) const;
const AmiVector &getQuadrature(realtype t) const;
virtual void reInit(realtype t0, const AmiVector &yy0,
const AmiVector &yp0) const = 0;
virtual void sensReInit(const AmiVectorArray &yyS0,
const AmiVectorArray &ypS0) const = 0;
virtual void sensToggleOff() const = 0;
virtual void reInitB(int which, realtype tB0, const AmiVector &yyB0,
const AmiVector &ypB0) const = 0;
virtual void quadReInitB(int which, const AmiVector &yQB0) const = 0;
realtype gett() const;
realtype getCpuTime() const;
realtype getCpuTimeB() const;
int nx() const;
int nplist() const;
int nquad() const;
bool computingFSA() const {
return getSensitivityOrder() >= SensitivityOrder::first &&
getSensitivityMethod() == SensitivityMethod::forward && nplist() > 0;
}
bool computingASA() const {
return getSensitivityOrder() >= SensitivityOrder::first &&
getSensitivityMethod() == SensitivityMethod::adjoint && nplist() > 0;
}
void resetDiagnosis() const;
void storeDiagnosis() const;
void storeDiagnosisB(int which) const;
std::vector<int> const& getNumSteps() const {
return ns_;
}
std::vector<int> const& getNumStepsB() const {
return nsB_;
}
std::vector<int> const& getNumRhsEvals() const {
return nrhs_;
}
std::vector<int> const& getNumRhsEvalsB() const {
return nrhsB_;
}
std::vector<int> const& getNumErrTestFails() const {
return netf_;
}
std::vector<int> const& getNumErrTestFailsB() const {
return netfB_;
}
std::vector<int> const& getNumNonlinSolvConvFails() const {
return nnlscf_;
}
std::vector<int> const& getNumNonlinSolvConvFailsB() const {
return nnlscfB_;
}
std::vector<int> const& getLastOrder() const {
return order_;
}
template <class Archive>
friend void boost::serialization::serialize(Archive &ar, Solver &s,
unsigned int version);
friend bool operator==(const Solver &a, const Solver &b);
AmiciApplication *app = &defaultContext;
protected:
virtual void setStopTime(realtype tstop) const = 0;
virtual int solve(realtype tout, int itask) const = 0;
virtual int solveF(realtype tout, int itask, int *ncheckPtr) const = 0;
virtual void reInitPostProcessF(realtype tnext) const = 0;
virtual void reInitPostProcessB(realtype tnext) const = 0;
virtual void getSens() const = 0;
virtual void getB(int which) const = 0;
virtual void getQuadB(int which) const = 0;
virtual void getQuad(realtype &t) const = 0;
virtual void init(realtype t0, const AmiVector &x0,
const AmiVector &dx0) const = 0;
virtual void initSteadystate(realtype t0, const AmiVector &x0,
const AmiVector &dx0) const = 0;
virtual void sensInit1(const AmiVectorArray &sx0,
const AmiVectorArray &sdx0) const = 0;
virtual void binit(int which, realtype tf, const AmiVector &xB0,
const AmiVector &dxB0) const = 0;
virtual void qbinit(int which, const AmiVector &xQB0) const = 0;
virtual void rootInit(int ne) const = 0;
void initializeNonLinearSolverSens(const Model *model) const;
virtual void setDenseJacFn() const = 0;
virtual void setSparseJacFn() const = 0;
virtual void setBandJacFn() const = 0;
virtual void setJacTimesVecFn() const = 0;
virtual void setDenseJacFnB(int which) const = 0;
virtual void setSparseJacFnB(int which) const = 0;
virtual void setBandJacFnB(int which) const = 0;
virtual void setJacTimesVecFnB(int which) const = 0;
virtual void setSparseJacFn_ss() const = 0;
virtual void allocateSolver() const = 0;
virtual void setSStolerances(double rtol, double atol) const = 0;
virtual void setSensSStolerances(double rtol, const double *atol) const = 0;
virtual void setSensErrCon(bool error_corr) const = 0;
virtual void setQuadErrConB(int which, bool flag) const = 0;
virtual void setQuadErrCon(bool flag) const = 0;
virtual void setErrHandlerFn() const = 0;
virtual void setUserData(Model *model) const = 0;
virtual void setUserDataB(int which, Model *model) const = 0;
virtual void setMaxNumSteps(long int mxsteps) const = 0;
virtual void setMaxNumStepsB(int which, long int mxstepsB) const = 0;
virtual void setStabLimDet(int stldet) const = 0;
virtual void setStabLimDetB(int which, int stldet) const = 0;
virtual void setId(const Model *model) const = 0;
virtual void setSuppressAlg(bool flag) const = 0;
virtual void setSensParams(const realtype *p, const realtype *pbar,
const int *plist) const = 0;
virtual void getDky(realtype t, int k) const = 0;
virtual void getDkyB(realtype t, int k, int which) const = 0;
virtual void getSensDky(realtype t, int k) const = 0;
virtual void getQuadDkyB(realtype t, int k, int which) const = 0;
virtual void getQuadDky(realtype t, int k) const = 0;
virtual void adjInit() const = 0;
virtual void quadInit(const AmiVector &xQ0) const = 0;
virtual void allocateSolverB(int *which) const = 0;
virtual void setSStolerancesB(int which, realtype relTolB,
realtype absTolB) const = 0;
virtual void quadSStolerancesB(int which, realtype reltolQB,
realtype abstolQB) const = 0;
virtual void quadSStolerances(realtype reltolQB,
realtype abstolQB) const = 0;
virtual void getNumSteps(const void *ami_mem, long int *numsteps) const = 0;
virtual void getNumRhsEvals(const void *ami_mem,
long int *numrhsevals) const = 0;
virtual void getNumErrTestFails(const void *ami_mem,
long int *numerrtestfails) const = 0;
virtual void
getNumNonlinSolvConvFails(const void *ami_mem,
long int *numnonlinsolvconvfails) const = 0;
virtual void getLastOrder(const void *ami_mem, int *order) const = 0;
void initializeLinearSolver(const Model *model) const;
void initializeNonLinearSolver() const;
virtual void setLinearSolver() const = 0;
virtual void setLinearSolverB(int which) const = 0;
virtual void setNonLinearSolver() const = 0;
virtual void setNonLinearSolverB(int which) const = 0;
virtual void setNonLinearSolverSens() const = 0;
void initializeLinearSolverB(const Model *model, int which) const;
void initializeNonLinearSolverB(int which) const;
virtual const Model *getModel() const = 0;
bool getInitDone() const;
bool getSensInitDone() const;
bool getAdjInitDone() const;
bool getInitDoneB(int which) const;
bool getQuadInitDoneB(int which) const;
bool getQuadInitDone() const;
virtual void diag() const = 0;
virtual void diagB(int which) const = 0;
void resetMutableMemory(int nx, int nplist, int nquad) const;
virtual void *getAdjBmem(void *ami_mem, int which) const = 0;
void applyTolerances() const;
void applyTolerancesFSA() const;
void applyTolerancesASA(int which) const;
void applyQuadTolerancesASA(int which) const;
void applyQuadTolerances() const;
void applySensitivityTolerances() const;
mutable std::unique_ptr<void, std::function<void(void *)>> solver_memory_;
mutable std::vector<std::unique_ptr<void, std::function<void(void *)>>>
solver_memory_B_;
InternalSensitivityMethod ism_ {InternalSensitivityMethod::simultaneous};
LinearMultistepMethod lmm_ {LinearMultistepMethod::BDF};
NonlinearSolverIteration iter_ {NonlinearSolverIteration::newton};
InterpolationType interp_type_ {InterpolationType::hermite};
long int maxsteps_ {10000};
mutable std::unique_ptr<SUNLinSolWrapper> linear_solver_;
mutable std::unique_ptr<SUNLinSolWrapper> linear_solver_B_;
mutable std::unique_ptr<SUNNonLinSolWrapper> non_linear_solver_;
mutable std::unique_ptr<SUNNonLinSolWrapper> non_linear_solver_B_;
mutable std::unique_ptr<SUNNonLinSolWrapper> non_linear_solver_sens_;
mutable bool solver_was_called_F_ {false};
mutable bool solver_was_called_B_ {false};
void setInitDone() const;
void setSensInitDone() const;
void setSensInitOff() const;
void setAdjInitDone() const;
void setInitDoneB(int which) const;
void setQuadInitDoneB(int which) const;
void setQuadInitDone() const;
void checkSensitivityMethod(const SensitivityMethod sensi_meth,
bool preequilibration) const;
mutable AmiVector x_ {0};
mutable AmiVector dky_ {0};
mutable AmiVector dx_ {0};
mutable AmiVectorArray sx_ {0, 0};
mutable AmiVectorArray sdx_ {0, 0};
mutable AmiVector xB_ {0};
mutable AmiVector dxB_ {0};
mutable AmiVector xQB_ {0};
mutable AmiVector xQ_ {0};
mutable realtype t_ {std::nan("")};
mutable bool force_reinit_postprocess_F_ {false};
mutable bool force_reinit_postprocess_B_ {false};
private:
void apply_max_num_steps() const;
void apply_max_num_steps_B() const;
SensitivityMethod sensi_meth_ {SensitivityMethod::forward};
SensitivityMethod sensi_meth_preeq_ {SensitivityMethod::forward};
booleantype stldet_ {true};
int ordering_ {static_cast<int>(SUNLinSolKLU::StateOrdering::AMD)};
long int newton_maxsteps_ {0L};
long int newton_maxlinsteps_ {0L};
NewtonDampingFactorMode newton_damping_factor_mode_
{NewtonDampingFactorMode::on};
realtype newton_damping_factor_lower_bound_ {1e-8};
bool requires_preequilibration_ {false};
LinearSolver linsol_ {LinearSolver::KLU};
realtype atol_ {1e-16};
realtype rtol_ {1e-8};
realtype atol_fsa_ {NAN};
realtype rtol_fsa_ {NAN};
realtype atolB_ {NAN};
realtype rtolB_ {NAN};
realtype quad_atol_ {1e-12};
realtype quad_rtol_ {1e-8};
realtype ss_atol_ {NAN};
realtype ss_rtol_ {NAN};
realtype ss_atol_sensi_ {NAN};
realtype ss_rtol_sensi_ {NAN};
RDataReporting rdata_mode_ {RDataReporting::full};
mutable realtype cpu_time_ {0.0};
mutable realtype cpu_timeB_ {0.0};
long int maxstepsB_ {0L};
SensitivityOrder sensi_ {SensitivityOrder::none};
mutable bool initialized_ {false};
mutable bool sens_initialized_ {false};
mutable bool adj_initialized_ {false};
mutable bool quad_initialized_ {false};
mutable std::vector<bool> initializedB_{false};
mutable std::vector<bool> initializedQB_{false};
mutable int ncheckPtr_ {0};
mutable std::vector<int> ns_;
mutable std::vector<int> nsB_;
mutable std::vector<int> nrhs_;
mutable std::vector<int> nrhsB_;
mutable std::vector<int> netf_;
mutable std::vector<int> netfB_;
mutable std::vector<int> nnlscf_;
mutable std::vector<int> nnlscfB_;
mutable std::vector<int> order_;
};
bool operator==(const Solver &a, const Solver &b);
void wrapErrHandlerFn(int error_code, const char *module, const char *function,
char *msg, void *eh_data);
} // namespace amici
#endif // AMICISOLVER_H