Stan Math Library  2.20.0
reverse mode automatic differentiation
idas_system.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_IDAS_RESIDUAL_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_IDAS_RESIDUAL_HPP
3 
4 #include <stan/math/rev/meta.hpp>
13 #include <idas/idas.h>
14 #include <nvector/nvector_serial.h>
15 #include <ostream>
16 #include <vector>
17 
18 #define CHECK_IDAS_CALL(call) idas_check(call, #call)
19 
26 inline void idas_check(int flag, const char* func) {
27  if (flag < 0) {
28  std::ostringstream ss;
29  ss << func << " failed with error flag " << flag;
30  throw std::runtime_error(ss.str());
31  }
32 }
33 
41 inline Eigen::MatrixXd matrix_d_from_NVarray(const N_Vector* nv,
42  const size_t& nv_size) {
43  size_t m = nv_size;
44  size_t n = NV_LENGTH_S(nv[0]);
45  stan::math::matrix_d res(n, m);
46  for (size_t j = 0; j < m; ++j) {
47  auto nvp = N_VGetArrayPointer(nv[j]);
48  for (size_t i = 0; i < n; ++i) {
49  res(i, j) = nvp[i];
50  }
51  }
52  return res;
53 }
54 
62 inline void matrix_d_to_NVarray(const Eigen::MatrixXd& mat, N_Vector* nv,
63  const size_t& nv_size) {
64  size_t m = nv_size;
65  size_t n = NV_LENGTH_S(nv[0]);
66  for (size_t j = 0; j < m; ++j) {
67  auto nvp = N_VGetArrayPointer(nv[j]);
68  for (size_t i = 0; i < n; ++i) {
69  nvp[i] = mat(i, j);
70  }
71  }
72 }
73 
74 namespace stan {
75 namespace math {
76 
89 template <typename F, typename Tyy, typename Typ, typename Tpar>
90 class idas_system {
91  protected:
92  const F& f_;
93  const std::vector<Tyy>& yy_;
94  const std::vector<Typ>& yp_;
95  std::vector<double> yy_val_; // workspace
96  std::vector<double> yp_val_; // workspace
97  const std::vector<Tpar>& theta_;
98  const std::vector<double>& x_r_;
99  const std::vector<int>& x_i_;
100  const size_t N_;
101  const size_t M_;
102  const size_t ns_; // nb. of sensi params
103  N_Vector nv_yy_;
104  N_Vector nv_yp_;
105  std::vector<double> rr_val_; // workspace
106  N_Vector nv_rr_;
107  N_Vector id_;
108  void* mem_;
109  std::ostream* msgs_;
110 
111  public:
112  static constexpr bool is_var_yy0 = stan::is_var<Tyy>::value;
113  static constexpr bool is_var_yp0 = stan::is_var<Typ>::value;
114  static constexpr bool is_var_par = stan::is_var<Tpar>::value;
115  static constexpr bool need_sens = is_var_yy0 || is_var_yp0 || is_var_par;
116 
118  using return_type = std::vector<std::vector<scalar_type> >;
119 
133  idas_system(const F& f, const std::vector<int>& eq_id,
134  const std::vector<Tyy>& yy0, const std::vector<Typ>& yp0,
135  const std::vector<Tpar>& theta, const std::vector<double>& x_r,
136  const std::vector<int>& x_i, std::ostream* msgs)
137  : f_(f),
138  yy_(yy0),
139  yp_(yp0),
140  yy_val_(value_of(yy0)),
141  yp_val_(value_of(yp0)),
142  theta_(theta),
143  x_r_(x_r),
144  x_i_(x_i),
145  N_(yy0.size()),
146  M_(theta.size()),
147  ns_((is_var_yy0 ? N_ : 0) + (is_var_yp0 ? N_ : 0)
148  + (is_var_par ? M_ : 0)),
149  nv_yy_(N_VMake_Serial(N_, yy_val_.data())),
150  nv_yp_(N_VMake_Serial(N_, yp_val_.data())),
151  rr_val_(N_, 0.0),
152  nv_rr_(N_VMake_Serial(N_, rr_val_.data())),
153  id_(N_VNew_Serial(N_)),
154  mem_(IDACreate()),
155  msgs_(msgs) {
156  if (nv_yy_ == NULL || nv_yp_ == NULL)
157  throw std::runtime_error("N_VMake_Serial failed to allocate memory");
158 
159  if (mem_ == NULL)
160  throw std::runtime_error("IDACreate failed to allocate memory");
161 
162  static const char* caller = "idas_system";
163  check_finite(caller, "initial state", yy0);
164  check_finite(caller, "derivative initial state", yp0);
165  check_finite(caller, "parameter vector", theta);
166  check_finite(caller, "continuous data", x_r);
167  check_nonzero_size(caller, "initial state", yy0);
168  check_nonzero_size(caller, "derivative initial state", yp0);
169  check_consistent_sizes(caller, "initial state", yy0,
170  "derivative initial state", yp0);
171  check_consistent_sizes(caller, "initial state", yy0,
172  "derivative-algebra id", eq_id);
173  check_greater_or_equal(caller, "derivative-algebra id", eq_id, 0);
174  check_less_or_equal(caller, "derivative-algebra id", eq_id, 1);
175 
176  for (size_t i = 0; i < N_; ++i)
177  NV_Ith_S(id_, i) = eq_id[i];
178  }
179 
184  N_VDestroy_Serial(nv_yy_);
185  N_VDestroy_Serial(nv_yp_);
186  N_VDestroy_Serial(nv_rr_);
187  N_VDestroy_Serial(id_);
188  IDAFree(&mem_);
189  }
190 
196  N_Vector& nv_yy() { return nv_yy_; }
197 
203  N_Vector& nv_yp() { return nv_yp_; }
204 
210  N_Vector& nv_rr() { return nv_rr_; }
211 
217  N_Vector& id() { return id_; }
218 
224  const std::vector<double>& yy_val() { return yy_val_; }
225 
231  const std::vector<double>& yp_val() { return yp_val_; }
232 
238  const std::vector<Tyy>& yy0() const { return yy_; }
239 
245  const std::vector<Typ>& yp0() const { return yp_; }
246 
252  const std::vector<Tpar>& theta() const { return theta_; }
253 
262  std::vector<scalar_type> vars() const {
263  std::vector<scalar_type> res;
264  if (is_var_yy0) {
265  res.insert(res.end(), yy0().begin(), yy0().end());
266  }
267  if (is_var_yp0) {
268  res.insert(res.end(), yp0().begin(), yp0().end());
269  }
270  if (is_var_par) {
271  res.insert(res.end(), theta().begin(), theta().end());
272  }
273 
274  return res;
275  }
276 
280  const size_t n() { return N_; }
281 
285  const size_t ns() { return ns_; }
286 
290  const size_t n_sys() { return N_ * (ns_ + 1); }
291 
295  const size_t n_par() { return theta_.size(); }
296 
300  void* mem() { return mem_; }
301 
305  const F& f() { return f_; }
306 
310  IDAResFn residual() { // a non-capture lambda
311  return [](double t, N_Vector yy, N_Vector yp, N_Vector rr,
312  void* user_data) -> int {
313  using DAE = idas_system<F, Tyy, Typ, Tpar>;
314  DAE* dae = static_cast<DAE*>(user_data);
315 
316  size_t N = NV_LENGTH_S(yy);
317  auto yy_val = N_VGetArrayPointer(yy);
318  std::vector<double> yy_vec(yy_val, yy_val + N);
319  auto yp_val = N_VGetArrayPointer(yp);
320  std::vector<double> yp_vec(yp_val, yp_val + N);
321  auto res = dae->f_(t, yy_vec, yp_vec, dae->theta_, dae->x_r_, dae->x_i_,
322  dae->msgs_);
323  for (size_t i = 0; i < N; ++i)
324  NV_Ith_S(rr, i) = value_of(res[i]);
325 
326  return 0;
327  };
328  }
329 
330  void check_ic_consistency(const double& t0, const double& tol) {
331  const std::vector<double> theta_d(value_of(theta_));
332  const std::vector<double> yy_d(value_of(yy_));
333  const std::vector<double> yp_d(value_of(yp_));
334  static const char* caller = "idas_integrator";
335  std::vector<double> res(f_(t0, yy_d, yp_d, theta_d, x_r_, x_i_, msgs_));
336  double res0 = std::sqrt(dot_self(res));
337  check_less_or_equal(caller, "DAE residual at t0", res0, tol);
338  }
339 };
340 
341 // TODO(yizhang): adjoint system construction
342 
343 } // namespace math
344 } // namespace stan
345 
346 #endif
static constexpr bool is_var_yy0
const std::vector< Typ > & yp0() const
return reference to derivative initial condition
const std::vector< int > & x_i_
Definition: idas_system.hpp:99
idas_system(const F &f, const std::vector< int > &eq_id, const std::vector< Tyy > &yy0, const std::vector< Typ > &yp0, const std::vector< Tpar > &theta, const std::vector< double > &x_r, const std::vector< int > &x_i, std::ostream *msgs)
Construct IDAS DAE system from initial condition and parameters.
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
static constexpr bool is_var_yp0
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.
void check_less_or_equal(const char *function, const char *name, const T_y &y, const T_high &high)
Check if y is less or equal to high.
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:13
Eigen::MatrixXd matrix_d_from_NVarray(const N_Vector *nv, const size_t &nv_size)
copy NV_Vector* array to Eigen::MatrixXd
Definition: idas_system.hpp:41
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
Defines a public enum named value which is defined to be false as the primitive scalar types cannot b...
Definition: is_var.hpp:10
IDAS DAE system that contains informtion on residual equation functor, sensitivity residual equation ...
Definition: idas_system.hpp:90
std::vector< scalar_type > vars() const
return a vector of vars for that contains the initial condition and parameters in case they are vars...
std::vector< double > yp_val_
Definition: idas_system.hpp:96
fvar< T > dot_self(const Eigen::Matrix< fvar< T >, R, C > &v)
Definition: dot_self.hpp:13
const size_t n()
return number of unknown variables
void matrix_d_to_NVarray(const Eigen::MatrixXd &mat, N_Vector *nv, const size_t &nv_size)
copy Eigen::MatrixXd to NV_Vector* array.
Definition: idas_system.hpp:62
void check_greater_or_equal(const char *function, const char *name, const T_y &y, const T_low &low)
Check if y is greater or equal than low.
~idas_system()
destructor to deallocate IDAS solution memory and workspace.
std::vector< double > yy_val_
Definition: idas_system.hpp:95
const std::vector< Tyy > & yy0() const
return reference to initial condition
void idas_check(int flag, const char *func)
check IDAS return flag & throw runtime error
Definition: idas_system.hpp:26
static constexpr bool is_var_par
void * mem()
return IDAS memory handle
boost::math::tools::promote_args< double, typename scalar_type< T >::type, typename return_type< Types_pack... >::type >::type type
Definition: return_type.hpp:36
const size_t n_sys()
return size of DAE system for primary and sensitivity unknowns
N_Vector & nv_yp()
return reference to current N_Vector of derivative variable
const size_t n_par()
return theta size
const std::vector< double > & x_r_
Definition: idas_system.hpp:98
const size_t ns()
return number of sensitivity parameters
const std::vector< Tyy > & yy_
Definition: idas_system.hpp:93
IDAResFn residual()
return a closure for IDAS residual callback
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > matrix_d
Type for matrix of double values.
Definition: typedefs.hpp:19
N_Vector & id()
return reference to DAE variable IDs
static constexpr bool need_sens
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
const std::vector< Tpar > & theta_
Definition: idas_system.hpp:97
std::vector< std::vector< scalar_type > > return_type
void check_ic_consistency(const double &t0, const double &tol)
N_Vector & nv_yy()
return reference to current N_Vector of unknown variable
const std::vector< Tpar > & theta() const
return reference to parameter
const std::vector< Typ > & yp_
Definition: idas_system.hpp:94
const std::vector< double > & yp_val()
return reference to current solution derivative vector value
const std::vector< double > & yy_val()
return reference to current solution vector value
void check_consistent_sizes(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)
Check if the dimension of x1 is consistent with x2.
std::vector< double > rr_val_
N_Vector & nv_rr()
return reference to current N_Vector of residual workspace
typename stan::return_type< Tyy, Typ, Tpar >::type scalar_type
const F & f()
return reference to DAE functor

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