Stan Math Library  2.20.0
reverse mode automatic differentiation
idas_integrator.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_IDAS_INTEGRATOR_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_IDAS_INTEGRATOR_HPP
3 
4 #include <stan/math/rev/meta.hpp>
10 #include <idas/idas.h>
11 #include <sunmatrix/sunmatrix_dense.h>
12 #include <sunlinsol/sunlinsol_dense.h>
13 #include <nvector/nvector_serial.h>
14 #include <ostream>
15 #include <vector>
16 #include <algorithm>
17 
19 
20 namespace stan {
21 namespace math {
22 
27  const double rtol_;
28  const double atol_;
29  const int64_t max_num_steps_;
33  template <typename Dae>
34  void init_sensitivity(Dae& dae);
35 
42  template <typename F>
43  void init_sensitivity(idas_forward_system<F, double, double, double>& dae) {}
44 
45  // /**
46  // * idas adjoint sens calculation requires different initialization
47  // *
48  // * @tparam F type of DAE RHS functor
49  // * @tparam Tyy type of DAE primary unknowns
50  // * @tparam Typ type of DAE derivative unknowns
51  // * @tparam Tpar type of DAE parameters.
52  // * @param[out] dae DAE system
53  // * @param[in] t0 initial time.
54  // * @param[in] ts times of the desired solutions
55  // * @param[out] res_yy DAE solutions
56  // */
57  // template <typename F, typename Tyy, typename Typ, typename Tpar>
58  // void init_sensitivity(idas_adjoint_system<F, Tyy, Typ, Tpar>& dae) {
59  // // TODO(yizhang): adjoint sensitivity initialization
60  // }
61 
62  template <typename F>
64  const double& t0, const std::vector<double>& ts,
65  std::vector<std::vector<double> >& res_yy);
66 
67  template <typename Dae>
68  void solve(Dae& dae, const double& t0, const std::vector<double>& ts,
69  typename Dae::return_type& res_yy);
70 
71  // TODO(yizhang): adjoint sensitivity solver
72 
73  public:
74  static constexpr int IDAS_MAX_STEPS = 500;
81  idas_integrator(const double rtol, const double atol,
82  const int64_t max_num_steps = IDAS_MAX_STEPS)
83  : rtol_(rtol), atol_(atol), max_num_steps_(max_num_steps) {
84  if (rtol_ <= 0)
85  invalid_argument("idas_integrator", "relative tolerance,", rtol_, "",
86  ", must be greater than 0");
87  if (rtol_ > 1.0E-3)
88  invalid_argument("idas_integrator", "relative tolerance,", rtol_, "",
89  ", must be less than 1.0E-3");
90  if (atol_ <= 0)
91  invalid_argument("idas_integrator", "absolute tolerance,", atol_, "",
92  ", must be greater than 0");
93  if (max_num_steps_ <= 0)
94  invalid_argument("idas_integrator", "max_num_steps,", max_num_steps_, "",
95  ", must be greater than 0");
96  }
97 
113  template <typename Dae>
114  typename Dae::return_type integrate(Dae& dae, double t0,
115  const std::vector<double>& ts) {
116  using Eigen::Dynamic;
117  using Eigen::Matrix;
118  using Eigen::MatrixXd;
119  using Eigen::VectorXd;
120 
121  static const char* caller = "idas_integrator";
122  check_finite(caller, "initial time", t0);
123  check_finite(caller, "times", ts);
124  check_ordered(caller, "times", ts);
125  check_nonzero_size(caller, "times", ts);
126  check_less(caller, "initial time", t0, ts.front());
127 
128  auto mem = dae.mem();
129  auto yy = dae.nv_yy();
130  auto yp = dae.nv_yp();
131  const size_t n = dae.n();
132 
133  typename Dae::return_type res_yy(
134  ts.size(), std::vector<typename Dae::scalar_type>(n, 0));
135 
136  auto A = SUNDenseMatrix(n, n);
137  auto LS = SUNDenseLinearSolver(yy, A);
138 
139  try {
140  CHECK_IDAS_CALL(IDASetUserData(mem, dae.to_user_data()));
141 
142  CHECK_IDAS_CALL(IDAInit(mem, dae.residual(), t0, yy, yp));
143  CHECK_IDAS_CALL(IDASetLinearSolver(mem, LS, A));
144  CHECK_IDAS_CALL(IDASStolerances(mem, rtol_, atol_));
145  CHECK_IDAS_CALL(IDASetMaxNumSteps(mem, max_num_steps_));
146 
147  init_sensitivity(dae);
148 
149  solve(dae, t0, ts, res_yy);
150  } catch (const std::exception& e) {
151  SUNLinSolFree(LS);
152  SUNMatDestroy(A);
153  throw;
154  }
155 
156  SUNLinSolFree(LS);
157  SUNMatDestroy(A);
158 
159  return res_yy;
160  }
161 }; // idas integrator
162 
171 template <typename Dae>
172 void idas_integrator::init_sensitivity(Dae& dae) {
173  if (Dae::need_sens) {
174  auto mem = dae.mem();
175  auto yys = dae.nv_yys();
176  auto yps = dae.nv_yps();
177  auto n = dae.n();
178 
179  if (Dae::is_var_yy0) {
180  for (size_t i = 0; i < n; ++i)
181  NV_Ith_S(yys[i], i) = 1.0;
182  }
183  if (Dae::is_var_yp0) {
184  for (size_t i = 0; i < n; ++i)
185  NV_Ith_S(yps[i + n], i) = 1.0;
186  }
187  CHECK_IDAS_CALL(IDASensInit(mem, dae.ns(), IDA_SIMULTANEOUS,
188  dae.sensitivity_residual(), yys, yps));
189  CHECK_IDAS_CALL(IDASensEEtolerances(mem));
190  CHECK_IDAS_CALL(IDAGetSensConsistentIC(mem, yys, yps));
191  }
192 }
193 
203 template <typename F>
204 void idas_integrator::solve(idas_forward_system<F, double, double, double>& dae,
205  const double& t0, const std::vector<double>& ts,
206  std::vector<std::vector<double> >& res_yy) {
207  double t1 = t0;
208  size_t i = 0;
209  auto mem = dae.mem();
210  auto yy = dae.nv_yy();
211  auto yp = dae.nv_yp();
212 
213  std::for_each(ts.begin(), ts.end(), [&](double t2) {
214  CHECK_IDAS_CALL(IDASolve(mem, t2, &t1, yy, yp, IDA_NORMAL));
215  std::copy(dae.yy_val().begin(), dae.yy_val().end(), res_yy[i].begin());
216  ++i;
217  });
218 }
219 
230 template <typename Dae>
231 void idas_integrator::solve(Dae& dae, const double& t0,
232  const std::vector<double>& ts,
233  typename Dae::return_type& res_yy) {
234  double t1 = t0;
235  size_t i = 0;
236  auto mem = dae.mem();
237  auto yy = dae.nv_yy();
238  auto yp = dae.nv_yp();
239  auto yys = dae.nv_yys();
240  const auto n = dae.n();
241  const auto ns = dae.ns();
242  auto vars = dae.vars();
243 
244  std::vector<stan::math::var> sol_t(n);
245  std::vector<double> sol_grad(ns);
246 
247  std::for_each(ts.begin(), ts.end(), [&](double t2) {
248  CHECK_IDAS_CALL(IDASolve(mem, t2, &t1, yy, yp, IDA_NORMAL));
249  CHECK_IDAS_CALL(IDAGetSens(mem, &t1, yys));
250  for (size_t k = 0; k < n; ++k) {
251  for (size_t j = 0; j < ns; ++j) {
252  sol_grad[j] = NV_Ith_S(yys[j], k);
253  }
254  sol_t[k]
255  = stan::math::precomputed_gradients(NV_Ith_S(yy, k), vars, sol_grad);
256  }
257  res_yy[i] = sol_t;
258  ++i;
259  });
260 }
261 } // namespace math
262 } // namespace stan
263 
264 #endif
IDAS DAE integrator.
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
void check_nonzero_size(const char *function, const char *name, const T_y &y)
Check if the specified matrix/vector is of non-zero size.
Dae::return_type integrate(Dae &dae, double t0, const std::vector< double > &ts)
Return the solutions for the specified DAE given the specified initial state, initial times...
void check_ordered(const char *function, const char *name, const std::vector< T_y > &y)
Check if the specified vector is sorted into strictly increasing order.
idas_integrator(const double rtol, const double atol, const int64_t max_num_steps=IDAS_MAX_STEPS)
constructor
IDAS DAE system with forward sensitivity calculation.
IDAS_SENSITIVITY
static constexpr int IDAS_MAX_STEPS
var precomputed_gradients(double value, const std::vector< var > &operands, const std::vector< double > &gradients)
This function returns a var for an expression that has the specified value, vector of operands...
void * mem()
return IDAS memory handle
N_Vector & nv_yp()
return reference to current N_Vector of derivative variable
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw an invalid_argument exception with a consistently formatted message.
#define CHECK_IDAS_CALL(call)
Definition: idas_system.hpp:18
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:87
const double E
The base of the natural logarithm, .
Definition: constants.hpp:19
N_Vector & nv_yy()
return reference to current N_Vector of unknown variable
const std::vector< double > & yy_val()
return reference to current solution vector value
void check_less(const char *function, const char *name, const T_y &y, const T_high &high)
Check if y is strictly less than high.
Definition: check_less.hpp:63
const kernel_cl< in_buffer, out_buffer, int, int > copy("copy", {indexing_helpers, copy_kernel_code})
See the docs for copy() .

     [ Stan Home Page ] © 2011–2018, Stan Development Team.