Stan Math Library  2.11.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 
6 #include <cvodes/cvodes.h>
7 #include <cvodes/cvodes_band.h>
8 #include <cvodes/cvodes_dense.h>
9 #include <nvector/nvector_serial.h>
10 #include <algorithm>
11 #include <vector>
12 
13 namespace stan {
14 
15  namespace math {
16 
26  template<typename F, typename T_initial, typename T_param>
28  const std::vector<T_initial>& y0_;
29  const std::vector<T_param>& theta_;
30  const size_t N_;
31  const size_t M_;
32  const size_t param_var_ind_;
33  const ode_system<F> ode_system_;
34 
37 
38  public:
53  cvodes_ode_data(const F& f,
54  const std::vector<T_initial>& y0,
55  const std::vector<T_param>& theta,
56  const std::vector<double>& x,
57  const std::vector<int>& x_int,
58  std::ostream* msgs)
59  : y0_(y0),
60  theta_(theta),
61  N_(y0.size()),
62  M_(theta.size()),
63  param_var_ind_(initial_var::value ? N_ : 0),
64  ode_system_(f, stan::math::value_of(theta), x, x_int, msgs) { }
65 
66  static int ode_rhs(double t, N_Vector y, N_Vector ydot, void* user_data) {
67  const ode_data* explicit_ode
68  = static_cast<const ode_data*>(user_data);
69  explicit_ode->rhs(NV_DATA_S(y), NV_DATA_S(ydot), t);
70  return 0;
71  }
72 
73  static int ode_rhs_sens(int Ns, realtype t,
74  N_Vector y, N_Vector ydot,
75  N_Vector *yS, N_Vector *ySdot, void *user_data,
76  N_Vector tmp1, N_Vector tmp2) {
77  const ode_data* explicit_ode = static_cast<const ode_data*>(user_data);
78  const std::vector<double> y_vec(NV_DATA_S(y),
79  NV_DATA_S(y) + explicit_ode->N_);
80  explicit_ode->rhs_sens(explicit_ode->y0_, explicit_ode->theta_,
81  t, y_vec, yS, ySdot);
82  return 0;
83  }
84 
85  static int dense_jacobian(long int N, // NOLINT(runtime/int)
86  realtype t, N_Vector y, N_Vector fy,
87  DlsMat J, void *user_data,
88  N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) {
89  const ode_data* explicit_ode = static_cast<const ode_data*>(user_data);
90  return explicit_ode->dense_jacobian(NV_DATA_S(y), J, t);
91  }
92 
93  private:
94  void rhs(const double y[], double dy_dt[], double t) const {
95  const std::vector<double> y_vec(y, y + N_);
96  std::vector<double> dy_dt_vec(N_);
97  ode_system_(t, y_vec, dy_dt_vec);
98  std::copy(dy_dt_vec.begin(), dy_dt_vec.end(), dy_dt);
99  }
100 
101  int dense_jacobian(const double* y, DlsMat J, double t) const {
102  const std::vector<double> y_vec(y, y + N_);
103  Eigen::VectorXd fy(N_);
104  // Eigen and CVODES use column major addressing
105  Eigen::Map<Eigen::MatrixXd> Jy_map(J->data, N_, N_);
106  ode_system_.jacobian(t, y_vec, fy, Jy_map);
107  return 0;
108  }
109 
110  inline void rhs_sens_initial(const Eigen::MatrixXd& Jy,
111  N_Vector *yS, N_Vector *ySdot) const {
112  for (size_t m = 0; m < N_; ++m) {
113  Eigen::Map<Eigen::VectorXd> yS_eig(NV_DATA_S(yS[m]), N_);
114  Eigen::Map<Eigen::VectorXd> ySdot_eig(NV_DATA_S(ySdot[m]), N_);
115  ySdot_eig = Jy * yS_eig;
116  }
117  }
118 
119  inline void rhs_sens_param(const Eigen::MatrixXd& Jy,
120  const Eigen::MatrixXd& Jtheta,
121  N_Vector *yS, N_Vector *ySdot) const {
122  using Eigen::Map;
123  using Eigen::VectorXd;
124  for (size_t m = 0; m < M_; ++m) {
125  Map<VectorXd> yS_eig(NV_DATA_S(yS[param_var_ind_ + m]), N_);
126  Map<VectorXd> ySdot_eig(NV_DATA_S(ySdot[param_var_ind_ + m]), N_);
127  ySdot_eig = Jy * yS_eig + Jtheta.col(m);
128  }
129  }
130 
142  void rhs_sens(const std::vector<stan::math::var>& initial,
143  const std::vector<stan::math::var>& param,
144  const double t, const std::vector<double>& y,
145  N_Vector *yS, N_Vector *ySdot) const {
146  Eigen::VectorXd dy_dt(N_);
147  Eigen::MatrixXd Jy(N_, N_);
148  Eigen::MatrixXd Jtheta(N_, M_);
149  ode_system_.jacobian(t, y, dy_dt, Jy, Jtheta);
150  rhs_sens_initial(Jy, yS, ySdot);
151  rhs_sens_param(Jy, Jtheta, yS, ySdot);
152  }
153 
165  void rhs_sens(const std::vector<double>& initial,
166  const std::vector<stan::math::var>& param,
167  const double t, const std::vector<double>& y,
168  N_Vector *yS, N_Vector *ySdot) const {
169  Eigen::VectorXd dy_dt(N_);
170  Eigen::MatrixXd Jy(N_, N_);
171  Eigen::MatrixXd Jtheta(N_, M_);
172  ode_system_.jacobian(t, y, dy_dt, Jy, Jtheta);
173  rhs_sens_param(Jy, Jtheta, yS, ySdot);
174  }
175 
187  void rhs_sens(const std::vector<stan::math::var>& initial,
188  const std::vector<double>& param,
189  const double t, const std::vector<double>& y,
190  N_Vector *yS, N_Vector *ySdot) const {
191  Eigen::VectorXd dy_dt(N_);
192  Eigen::MatrixXd Jy(N_, N_);
193  ode_system_.jacobian(t, y, dy_dt, Jy);
194  rhs_sens_initial(Jy, yS, ySdot);
195  }
196 
208  void rhs_sens(const std::vector<double>& initial,
209  const std::vector<double>& param,
210  const double t, const std::vector<double>& y,
211  N_Vector *yS, N_Vector *ySdot) const {
212  }
213  };
214 
215  }
216 }
217 #endif
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:16
static int ode_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)
void jacobian(const double t, const std::vector< double > &y, Eigen::MatrixBase< Derived1 > &dy_dt, Eigen::MatrixBase< Derived2 > &Jy) const
Calculate the Jacobian of the ODE RHS wrt to states y.
Definition: ode_system.hpp:67
CVODES ode data holder object which is used during CVODES integration for CVODES callbacks.
static int dense_jacobian(long int N, realtype t, N_Vector y, N_Vector fy, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Internal representation of an ODE model object which provides convenient Jacobian functions to obtain...
Definition: ode_system.hpp:21
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 int ode_rhs(double t, N_Vector y, N_Vector ydot, void *user_data)

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