1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_CVODES_ODE_DATA_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_CVODES_ODE_DATA_HPP
6 #include <cvodes/cvodes.h>
7 #include <cvodes/cvodes_band.h>
8 #include <cvodes/cvodes_dense.h>
9 #include <nvector/nvector_serial.h>
26 template<
typename F,
typename T_initial,
typename T_param>
28 const std::vector<T_initial>& y0_;
29 const std::vector<T_param>& theta_;
32 const size_t param_var_ind_;
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,
63 param_var_ind_(initial_var::value ? N_ : 0),
64 ode_system_(f,
stan::math::
value_of(theta), x, x_int, msgs) { }
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);
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_,
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);
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);
102 const std::vector<double> y_vec(y, y + N_);
103 Eigen::VectorXd fy(N_);
105 Eigen::Map<Eigen::MatrixXd> Jy_map(J->data, N_, N_);
106 ode_system_.
jacobian(t, y_vec, fy, Jy_map);
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;
119 inline void rhs_sens_param(
const Eigen::MatrixXd& Jy,
120 const Eigen::MatrixXd& Jtheta,
121 N_Vector *yS, N_Vector *ySdot)
const {
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);
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);
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);
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);
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 {
T value_of(const fvar< T > &v)
Return the value of the specified variable.
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.
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...
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.
static int ode_rhs(double t, N_Vector y, N_Vector ydot, void *user_data)