| Stan Math Library 2.12.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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

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