Stan Math Library  2.11.0
reverse mode automatic differentiation
integrate_ode_bdf.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_INTEGRATE_ODE_BDF_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_INTEGRATE_ODE_BDF_HPP
3 
14 #include <cvodes/cvodes.h>
15 #include <cvodes/cvodes_band.h>
16 #include <cvodes/cvodes_dense.h>
17 #include <nvector/nvector_serial.h>
18 #include <algorithm>
19 #include <ostream>
20 #include <vector>
21 
22 namespace stan {
23  namespace math {
24 
34  inline void free_cvodes_memory(N_Vector& cvodes_state,
35  N_Vector* cvodes_state_sens,
36  void* cvodes_mem, size_t S) {
37  N_VDestroy_Serial(cvodes_state);
38  if (cvodes_state_sens != NULL)
39  N_VDestroyVectorArray_Serial(cvodes_state_sens, S);
40  CVodeFree(&cvodes_mem);
41  }
42 
80  template <typename F, typename T_initial, typename T_param>
81  std::vector<std::vector<typename stan::return_type<T_initial,
82  T_param>::type> >
83  integrate_ode_bdf(const F& f,
84  const std::vector<T_initial>& y0,
85  const double t0,
86  const std::vector<double>& ts,
87  const std::vector<T_param>& theta,
88  const std::vector<double>& x,
89  const std::vector<int>& x_int,
90  std::ostream* msgs = 0,
91  double relative_tolerance = 1e-10,
92  double absolute_tolerance = 1e-10,
93  long int max_num_steps = 1e8) { // NOLINT(runtime/int)
94  typedef stan::is_var<T_initial> initial_var;
95  typedef stan::is_var<T_param> param_var;
96 
97  stan::math::check_finite("integrate_ode_bdf", "initial state", y0);
98  stan::math::check_finite("integrate_ode_bdf", "initial time", t0);
99  stan::math::check_finite("integrate_ode_bdf", "times", ts);
100  stan::math::check_finite("integrate_ode_bdf", "parameter vector", theta);
101  stan::math::check_finite("integrate_ode_bdf", "continuous data", x);
102  stan::math::check_nonzero_size("integrate_ode_bdf", "times", ts);
103  stan::math::check_nonzero_size("integrate_ode_bdf", "initial state", y0);
104  stan::math::check_ordered("integrate_ode_bdf", "times", ts);
105  stan::math::check_less("integrate_ode_bdf", "initial time", t0, ts[0]);
106  if (relative_tolerance <= 0)
107  invalid_argument("integrate_ode_bdf",
108  "relative_tolerance,", relative_tolerance,
109  "", ", must be greater than 0");
110  if (absolute_tolerance <= 0)
111  invalid_argument("integrate_ode_bdf",
112  "absolute_tolerance,", absolute_tolerance,
113  "", ", must be greater than 0");
114  if (max_num_steps <= 0)
115  invalid_argument("integrate_ode_bdf",
116  "max_num_steps,", max_num_steps,
117  "", ", must be greater than 0");
118 
119  const size_t N = y0.size();
120  const size_t M = theta.size();
121  // total number of sensitivities for initial values and params
122  const size_t S = (initial_var::value ? N : 0)
123  + (param_var::value ? M : 0);
124  const size_t size = N * (S + 1); // size of the coupled system
125  std::vector<double> state(value_of(y0));
126  N_Vector cvodes_state(N_VMake_Serial(N, &state[0]));
127  N_Vector* cvodes_state_sens = NULL;
128 
130  ode_data cvodes_data(f, y0, theta, x, x_int, msgs);
131 
132  void* cvodes_mem = CVodeCreate(CV_BDF, CV_NEWTON);
133  if (cvodes_mem == NULL)
134  throw std::runtime_error("CVodeCreate failed to allocate memory");
135 
136  std::vector<std::vector<double> >
137  y_coupled(ts.size(), std::vector<double>(size, 0));
138 
139  try {
140  cvodes_check_flag(CVodeInit(cvodes_mem, &ode_data::ode_rhs,
141  t0, cvodes_state),
142  "CVodeInit");
143 
144  // Assign pointer to this as user data
145  cvodes_check_flag(CVodeSetUserData(cvodes_mem,
146  reinterpret_cast<void*>(&cvodes_data)),
147  "CVodeSetUserData");
148 
149  cvodes_set_options(cvodes_mem,
150  relative_tolerance, absolute_tolerance,
151  max_num_steps);
152 
153  // for the stiff solvers we need to reserve additional
154  // memory and provide a Jacobian function call
155  cvodes_check_flag(CVDense(cvodes_mem, N), "CVDense");
156  cvodes_check_flag(CVDlsSetDenseJacFn(cvodes_mem,
157  &ode_data::dense_jacobian),
158  "CVDlsSetDenseJacFn");
159 
160  // initialize forward sensitivity system of CVODES as needed
161  if (S > 0) {
162  cvodes_state_sens = N_VCloneVectorArray_Serial(S, cvodes_state);
163  for (size_t s = 0; s < S; s++)
164  N_VConst(RCONST(0.0), cvodes_state_sens[s]);
165 
166  // for varying initials, first N sensitivity systems
167  // are for initials which have as initial the identity matrix
168  if (initial_var::value) {
169  for (size_t n = 0; n < N; n++)
170  NV_Ith_S(cvodes_state_sens[n], n) = 1.0;
171  }
172  cvodes_check_flag(CVodeSensInit(cvodes_mem, static_cast<int>(S),
173  CV_STAGGERED,
174  &ode_data::ode_rhs_sens,
175  cvodes_state_sens),
176  "CVodeSensInit");
177 
178  cvodes_check_flag(CVodeSensEEtolerances(cvodes_mem),
179  "CVodeSensEEtolerances");
180  }
181 
182  double t_init = t0;
183  for (size_t n = 0; n < ts.size(); ++n) {
184  double t_final = ts[n];
185  if (t_final != t_init)
186  cvodes_check_flag(CVode(cvodes_mem, t_final, cvodes_state,
187  &t_init, CV_NORMAL),
188  "CVode");
189  std::copy(state.begin(), state.end(), y_coupled[n].begin());
190  if (S > 0) {
191  cvodes_check_flag(CVodeGetSens(cvodes_mem, &t_init,
192  cvodes_state_sens),
193  "CVodeGetSens");
194  for (size_t s = 0; s < S; s++)
195  std::copy(NV_DATA_S(cvodes_state_sens[s]),
196  NV_DATA_S(cvodes_state_sens[s]) + N,
197  y_coupled[n].begin() + N + s * N);
198  }
199  t_init = t_final;
200  }
201  } catch (const std::exception& e) {
202  free_cvodes_memory(cvodes_state, cvodes_state_sens, cvodes_mem, S);
203  throw;
204  }
205 
206  free_cvodes_memory(cvodes_state, cvodes_state_sens, cvodes_mem, S);
207 
208  return decouple_ode_states(y_coupled, y0, theta);
209  }
210 
211  }
212 }
213 #endif
bool check_less(const char *function, const char *name, const T_y &y, const T_high &high)
Return true if y is strictly less than high.
Definition: check_less.hpp:81
void cvodes_set_options(void *cvodes_mem, double rel_tol, double abs_tol, long int max_num_steps)
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:16
Metaprogram to calculate the base scalar return type resulting from promoting all the scalar types of...
Definition: return_type.hpp:19
std::vector< std::vector< typename stan::return_type< T_initial, T_param >::type > > integrate_ode_bdf(const F &f, const std::vector< T_initial > &y0, const double t0, const std::vector< double > &ts, const std::vector< T_param > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs=0, double relative_tolerance=1e-10, double absolute_tolerance=1e-10, long int max_num_steps=1e8)
Return the solutions for the specified system of ordinary differential equations given the specified ...
void free_cvodes_memory(N_Vector &cvodes_state, N_Vector *cvodes_state_sens, void *cvodes_mem, size_t S)
Free memory allocated for CVODES state, sensitivity, and general memory.
bool check_nonzero_size(const char *function, const char *name, const T_y &y)
Return true if the specified matrix/vector is of non-zero size.
bool check_ordered(const char *function, const char *name, const std::vector< T_y > &y)
Return true if the specified vector is sorted into strictly increasing order.
CVODES ode data holder object which is used during CVODES integration for CVODES callbacks.
void cvodes_check_flag(int flag, const std::string &func_name)
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw an invalid_argument exception with a consistently formatted message.
std::vector< std::vector< typename stan::return_type< T_initial, T_param >::type > > decouple_ode_states(const std::vector< std::vector< double > > &y, const std::vector< T_initial > &y0, const std::vector< T_param > &theta)
Takes sensitivity output from integrators and returns results in precomputed_gradients format...
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:95
bool check_finite(const char *function, const char *name, const T_y &y)
Return true if y is finite.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17

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