Automatic Differentiation
 
Loading...
Searching...
No Matches
idas_integrator.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_REV_FUNCTOR_IDAS_INTEGRATOR_HPP
2#define STAN_MATH_REV_FUNCTOR_IDAS_INTEGRATOR_HPP
3
8#include <idas/idas.h>
9#include <sunmatrix/sunmatrix_dense.h>
10#include <sunlinsol/sunlinsol_dense.h>
11#include <nvector/nvector_serial.h>
12#include <ostream>
13#include <vector>
14#include <algorithm>
15
16namespace stan {
17namespace math {
18
23 sundials::Context sundials_context_;
24
25 const double rtol_;
26 const double atol_;
27 const int64_t max_num_steps_;
28
29 public:
36 idas_integrator(const double rtol, const double atol,
37 const int64_t max_num_steps)
38 : rtol_(rtol), atol_(atol), max_num_steps_(max_num_steps) {}
39
56 template <typename dae_type>
57 typename dae_type::return_t operator()(const char* func, dae_type& dae,
58 double t0,
59 const std::vector<double>& ts) {
61
62 void* mem = serv.mem;
63 N_Vector& yy = serv.nv_yy;
64 N_Vector& yp = serv.nv_yp;
65 N_Vector* yys = serv.nv_yys;
66 N_Vector* yps = serv.nv_yps;
67 const size_t n = dae.N;
68
69 CHECK_IDAS_CALL(IDASStolerances(mem, rtol_, atol_));
70 CHECK_IDAS_CALL(IDASetMaxNumSteps(mem, max_num_steps_));
71
72 if (dae_type::use_fwd_sens) {
73 CHECK_IDAS_CALL(IDASensEEtolerances(mem));
74 CHECK_IDAS_CALL(IDAGetSensConsistentIC(mem, yys, yps));
75 }
76
77 size_t nt = ts.size();
78 typename dae_type::return_t res_yy(
79 nt, Eigen::Matrix<typename dae_type::scalar_t, -1, 1>::Zero(n));
80 double t1 = t0;
81 for (size_t i = 0; i < nt; ++i) {
82 CHECK_IDAS_CALL(IDASolve(mem, ts[i], &t1, yy, yp, IDA_NORMAL));
83 if (dae_type::use_fwd_sens) {
84 CHECK_IDAS_CALL(IDAGetSens(mem, &t1, yys));
85 }
86 collect(yy, yys, dae, res_yy[i]);
87 }
88
89 return res_yy;
90 }
91
101 template <typename dae_type>
102 void collect(N_Vector const& yy, N_Vector const* yys, dae_type& dae,
103 Eigen::Matrix<double, -1, 1>& y) {
104 for (int i = 0; i < dae.N; ++i) {
105 y.coeffRef(i) = NV_Ith_S(yy, i);
106 }
107 }
108
109 template <typename dae_type>
110 void collect(N_Vector const& yy, N_Vector const* yys, dae_type& dae,
111 Eigen::Matrix<stan::math::var, -1, 1>& y) {
112 std::vector<double> g(dae.ns);
113 for (size_t i = 0; i < dae.N; ++i) {
114 for (size_t j = 0; j < dae.ns; ++j) {
115 g[j] = NV_Ith_S(yys[j], i);
116 }
117 y[i] = precomputed_gradients(NV_Ith_S(yy, i), dae.all_vars, g);
118 }
119 }
120};
121
122} // namespace math
123} // namespace stan
124
125#endif
#define CHECK_IDAS_CALL(call)
idas_integrator(const double rtol, const double atol, const int64_t max_num_steps)
constructor
sundials::Context sundials_context_
dae_type::return_t operator()(const char *func, dae_type &dae, double t0, const std::vector< double > &ts)
Return the solutions for the specified DAE given the specified initial state, initial times,...
void collect(N_Vector const &yy, N_Vector const *yys, dae_type &dae, Eigen::Matrix< stan::math::var, -1, 1 > &y)
void collect(N_Vector const &yy, N_Vector const *yys, dae_type &dae, Eigen::Matrix< double, -1, 1 > &y)
Solve DAE system, no sensitivity.
std::vector< Eigen::Matrix< stan::return_type_t< T_yy, T_yp, T_Args... >, -1, 1 > > dae(const F &f, const T_yy &yy0, const T_yp &yp0, double t0, const std::vector< double > &ts, std::ostream *msgs, const T_Args &... args)
Solve the DAE initial value problem f(t, y, y')=0, y(t0) = yy0, y'(t0)=yp0 at a set of times,...
Definition dae.hpp:167
var precomputed_gradients(Arith value, const VecVar &operands, const VecArith &gradients, const std::tuple< ContainerOperands... > &container_operands=std::tuple<>(), const std::tuple< ContainerGradients... > &container_gradients=std::tuple<>())
This function returns a var for an expression that has the specified value, vector of operands,...
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
For each type of Ode(with different rhs functor F and senstivity parameters), we allocate mem and wor...