Stan Math Library  2.20.0
reverse mode automatic differentiation
cvodes_ode_data.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_CVODES_ODE_DATA_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_CVODES_ODE_DATA_HPP
3 
4 #include <stan/math/rev/meta.hpp>
7 #include <cvodes/cvodes.h>
8 #include <sunmatrix/sunmatrix_dense.h>
9 #include <sunlinsol/sunlinsol_dense.h>
10 #include <nvector/nvector_serial.h>
11 #include <algorithm>
12 #include <vector>
13 
14 namespace stan {
15 namespace math {
16 
26 template <typename F, typename T_initial, typename T_param>
28  const F& f_;
29  const std::vector<T_initial>& y0_;
30  const std::vector<T_param>& theta_;
31  const std::vector<double> theta_dbl_;
32  const size_t N_;
33  const size_t M_;
34  const std::vector<double>& x_;
35  const std::vector<int>& x_int_;
36  std::ostream* msgs_;
37  const size_t S_;
38 
42 
43  public:
45  std::vector<double> coupled_state_;
46  N_Vector nv_state_;
47  N_Vector* nv_state_sens_;
48  SUNMatrix A_;
49  SUNLinearSolver LS_;
50 
73  cvodes_ode_data(const F& f, const std::vector<T_initial>& y0,
74  const std::vector<T_param>& theta,
75  const std::vector<double>& x, const std::vector<int>& x_int,
76  std::ostream* msgs)
77  : f_(f),
78  y0_(y0),
79  theta_(theta),
80  theta_dbl_(value_of(theta)),
81  N_(y0.size()),
82  M_(theta.size()),
83  x_(x),
84  x_int_(x_int),
85  msgs_(msgs),
86  S_((initial_var::value ? N_ : 0) + (param_var::value ? M_ : 0)),
87  coupled_ode_(f, y0, theta, x, x_int, msgs),
88  coupled_state_(coupled_ode_.initial_state()),
89  nv_state_(N_VMake_Serial(N_, &coupled_state_[0])),
90  nv_state_sens_(nullptr),
91  A_(SUNDenseMatrix(N_, N_)),
92  LS_(SUNDenseLinearSolver(nv_state_, A_)) {
93  if (S_ > 0) {
94  nv_state_sens_ = N_VCloneVectorArrayEmpty_Serial(S_, nv_state_);
95  for (std::size_t i = 0; i < S_; i++) {
96  NV_DATA_S(nv_state_sens_[i]) = &coupled_state_[N_] + i * N_;
97  }
98  }
99  }
100 
102  SUNLinSolFree(LS_);
103  SUNMatDestroy(A_);
104  N_VDestroy_Serial(nv_state_);
105  if (S_ > 0)
106  N_VDestroyVectorArray_Serial(nv_state_sens_, S_);
107  }
108 
113  static int cv_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data) {
114  const ode_data* explicit_ode = static_cast<const ode_data*>(user_data);
115  explicit_ode->rhs(t, NV_DATA_S(y), NV_DATA_S(ydot));
116  return 0;
117  }
118 
123  static int cv_rhs_sens(int Ns, realtype t, N_Vector y, N_Vector ydot,
124  N_Vector* yS, N_Vector* ySdot, void* user_data,
125  N_Vector tmp1, N_Vector tmp2) {
126  const ode_data* explicit_ode = static_cast<const ode_data*>(user_data);
127  explicit_ode->rhs_sens(t, NV_DATA_S(y), yS, ySdot);
128  return 0;
129  }
130 
137  static int cv_jacobian_states(realtype t, N_Vector y, N_Vector fy,
138  SUNMatrix J, void* user_data, N_Vector tmp1,
139  N_Vector tmp2, N_Vector tmp3) {
140  const ode_data* explicit_ode = static_cast<const ode_data*>(user_data);
141  return explicit_ode->jacobian_states(t, NV_DATA_S(y), J);
142  }
143 
144  private:
149  inline void rhs(double t, const double y[], double dy_dt[]) const {
150  const std::vector<double> y_vec(y, y + N_);
151  const std::vector<double>& dy_dt_vec
152  = f_(t, y_vec, theta_dbl_, x_, x_int_, msgs_);
153  check_size_match("cvodes_ode_data", "dz_dt", dy_dt_vec.size(), "states",
154  N_);
155  std::move(dy_dt_vec.begin(), dy_dt_vec.end(), dy_dt);
156  }
157 
165  inline int jacobian_states(double t, const double y[], SUNMatrix J) const {
166  start_nested();
167  const std::vector<var> y_vec_var(y, y + N_);
168  coupled_ode_system<F, var, double> ode_jacobian(f_, y_vec_var, theta_dbl_,
169  x_, x_int_, msgs_);
170  std::vector<double>&& jacobian_y = std::vector<double>(ode_jacobian.size());
171  ode_jacobian(ode_jacobian.initial_state(), jacobian_y, t);
172  std::move(jacobian_y.begin() + N_, jacobian_y.end(), SM_DATA_D(J));
174  return 0;
175  }
176 
183  inline void rhs_sens(double t, const double y[], N_Vector* yS,
184  N_Vector* ySdot) const {
185  std::vector<double> z(coupled_state_.size());
186  std::vector<double>&& dz_dt = std::vector<double>(coupled_state_.size());
187  std::copy(y, y + N_, z.begin());
188  for (std::size_t s = 0; s < S_; s++)
189  std::copy(NV_DATA_S(yS[s]), NV_DATA_S(yS[s]) + N_,
190  z.begin() + (s + 1) * N_);
191  coupled_ode_(z, dz_dt, t);
192  for (std::size_t s = 0; s < S_; s++)
193  std::move(dz_dt.begin() + (s + 1) * N_, dz_dt.begin() + (s + 2) * N_,
194  NV_DATA_S(ySdot[s]));
195  }
196 };
197 
198 } // namespace math
199 } // namespace stan
200 #endif
The coupled_ode_system template specialization for unknown initial values and known parameters...
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
static int cv_rhs_sens(int Ns, realtype t, N_Vector y, N_Vector ydot, N_Vector *yS, N_Vector *ySdot, void *user_data, N_Vector tmp1, N_Vector tmp2)
Implements the function of type CVSensRhsFn which is the RHS of the sensitivity ODE system...
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
std::vector< double > coupled_state_
static int cv_rhs(realtype t, N_Vector y, N_Vector ydot, void *user_data)
Implements the function of type CVRhsFn which is the user-defined ODE RHS passed to CVODES...
CVODES ode data holder object which is used during CVODES integration for CVODES callbacks.
const coupled_ode_system< F, T_initial, T_param > coupled_ode_
cvodes_ode_data(const F &f, const std::vector< T_initial > &y0, const std::vector< T_param > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs)
Construct CVODES ode data object to enable callbacks from CVODES during ODE integration.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
static void recover_memory_nested()
Recover only the memory used for the top nested call.
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
static int cv_jacobian_states(realtype t, N_Vector y, N_Vector fy, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Implements the function of type CVDlsJacFn which is the user-defined callback for CVODES to calculate...
static void start_nested()
Record the current position so that recover_memory_nested() can find it.
size_t size() const
Returns the size of the coupled system.
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.