Stan Math Library  2.11.0
reverse mode automatic differentiation
coupled_ode_system.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_ARR_FUNCTOR_COUPLED_ODE_SYSTEM_HPP
2 #define STAN_MATH_REV_ARR_FUNCTOR_COUPLED_ODE_SYSTEM_HPP
3 
10 #include <stan/math/rev/core.hpp>
11 #include <ostream>
12 #include <stdexcept>
13 #include <vector>
14 
15 namespace stan {
16  namespace math {
17 
18  // This code is in this directory because it includes var
19  // It is in namespace stan::math so that the partial template
20  // specializations are treated as such.
21 
22 
34  inline void add_initial_values(const std::vector<stan::math::var>& y0,
35  std::vector<std::vector<stan::math::var> >& y) {
36  for (size_t n = 0; n < y.size(); n++)
37  for (size_t m = 0; m < y0.size(); m++)
38  y[n][m] += y0[m];
39  }
40 
62  template <typename F>
63  struct coupled_ode_system <F, double, stan::math::var> {
64  const F& f_;
65  const std::vector<double>& y0_dbl_;
66  const std::vector<stan::math::var>& theta_;
67  std::vector<double> theta_dbl_;
68  const std::vector<double>& x_;
69  const std::vector<int>& x_int_;
70  const size_t N_;
71  const size_t M_;
72  const size_t size_;
73  std::ostream* msgs_;
74 
87  coupled_ode_system(const F& f,
88  const std::vector<double>& y0,
89  const std::vector<stan::math::var>& theta,
90  const std::vector<double>& x,
91  const std::vector<int>& x_int,
92  std::ostream* msgs)
93  : f_(f),
94  y0_dbl_(y0),
95  theta_(theta),
96  theta_dbl_(theta.size(), 0.0),
97  x_(x),
98  x_int_(x_int),
99  N_(y0.size()),
100  M_(theta.size()),
101  size_(N_ + N_ * M_),
102  msgs_(msgs) {
103  for (size_t m = 0; m < M_; m++)
104  theta_dbl_[m] = stan::math::value_of(theta[m]);
105  }
106 
124  void operator()(const std::vector<double>& z,
125  std::vector<double>& dz_dt,
126  double t) {
127  using std::vector;
128  using stan::math::var;
129 
130  vector<double> y(z.begin(), z.begin() + N_);
131  dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
132  stan::math::check_equal("coupled_ode_system",
133  "dz_dt", dz_dt.size(), N_);
134 
135  vector<double> coupled_sys(N_ * M_);
136  vector<double> grad(N_ + M_);
137 
138  try {
140 
141  vector<var> z_vars;
142  z_vars.reserve(N_ + M_);
143 
144  vector<var> y_vars(y.begin(), y.end());
145  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
146 
147  vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
148  z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
149 
150  vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
151 
152  for (size_t i = 0; i < N_; i++) {
154  dy_dt_vars[i].grad(z_vars, grad);
155 
156  for (size_t j = 0; j < M_; j++) {
157  // orders derivatives by equation (i.e. if there are 2 eqns
158  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
159  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
160  double temp_deriv = grad[N_ + j];
161  for (size_t k = 0; k < N_; k++)
162  temp_deriv += z[N_ + N_ * j + k] * grad[k];
163 
164  coupled_sys[i + j * N_] = temp_deriv;
165  }
166  }
167  } catch (const std::exception& e) {
169  throw;
170  }
172 
173  dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
174  }
175 
181  size_t size() const {
182  return size_;
183  }
184 
198  std::vector<double> initial_state() {
199  std::vector<double> state(size_, 0.0);
200  for (size_t n = 0; n < N_; n++)
201  state[n] = y0_dbl_[n];
202  return state;
203  }
204 
211  std::vector<std::vector<stan::math::var> >
212  decouple_states(const std::vector<std::vector<double> >& y) {
214  std::vector<stan::math::var> temp_vars(N_);
215  std::vector<double> temp_gradients(M_);
216  std::vector<std::vector<stan::math::var> > y_return(y.size());
217 
218  for (size_t i = 0; i < y.size(); i++) {
219  // iterate over number of equations
220  for (size_t j = 0; j < N_; j++) {
221  // iterate over parameters for each equation
222  for (size_t k = 0; k < M_; k++)
223  temp_gradients[k] = y[i][y0_dbl_.size() + y0_dbl_.size() * k + j];
224 
225  temp_vars[j] = precomputed_gradients(y[i][j],
226  theta_,
227  temp_gradients);
228  }
229  y_return[i] = temp_vars;
230  }
231  return y_return;
232  }
233  };
234 
261  template <typename F>
262  struct coupled_ode_system <F, stan::math::var, double> {
263  const F& f_;
264  const std::vector<stan::math::var>& y0_;
265  std::vector<double> y0_dbl_;
266  const std::vector<double>& theta_dbl_;
267  const std::vector<double>& x_;
268  const std::vector<int>& x_int_;
269  std::ostream* msgs_;
270  const size_t N_;
271  const size_t M_;
272  const size_t size_;
273 
287  coupled_ode_system(const F& f,
288  const std::vector<stan::math::var>& y0,
289  const std::vector<double>& theta,
290  const std::vector<double>& x,
291  const std::vector<int>& x_int,
292  std::ostream* msgs)
293  : f_(f),
294  y0_(y0),
295  y0_dbl_(y0.size(), 0.0),
296  theta_dbl_(theta),
297  x_(x),
298  x_int_(x_int),
299  msgs_(msgs),
300  N_(y0.size()),
301  M_(theta.size()),
302  size_(N_ + N_ * N_) {
303  for (size_t n = 0; n < N_; n++)
304  y0_dbl_[n] = stan::math::value_of(y0_[n]);
305  }
306 
323  void operator()(const std::vector<double>& z,
324  std::vector<double>& dz_dt,
325  double t) {
326  using std::vector;
327  using stan::math::var;
328 
329  std::vector<double> y(z.begin(), z.begin() + N_);
330  for (size_t n = 0; n < N_; n++)
331  y[n] += y0_dbl_[n];
332 
333  dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
334  stan::math::check_equal("coupled_ode_system",
335  "dz_dt", dz_dt.size(), N_);
336 
337  std::vector<double> coupled_sys(N_ * N_);
338  std::vector<double> grad(N_);
339 
340  try {
342 
343  vector<var> z_vars;
344  z_vars.reserve(N_);
345 
346  vector<var> y_vars(y.begin(), y.end());
347  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
348 
349  vector<var> dy_dt_vars = f_(t, y_vars, theta_dbl_, x_, x_int_, msgs_);
350 
351  for (size_t i = 0; i < N_; i++) {
353  dy_dt_vars[i].grad(z_vars, grad);
354 
355  for (size_t j = 0; j < N_; j++) {
356  // orders derivatives by equation (i.e. if there are 2 eqns
357  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
358  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
359  double temp_deriv = grad[j];
360  for (size_t k = 0; k < N_; k++)
361  temp_deriv += z[N_ + N_ * j + k] * grad[k];
362 
363  coupled_sys[i + j * N_] = temp_deriv;
364  }
365  }
366  } catch (const std::exception& e) {
368  throw;
369  }
371 
372  dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
373  }
374 
380  size_t size() const {
381  return size_;
382  }
383 
398  std::vector<double> initial_state() {
399  return std::vector<double>(size_, 0.0);
400  }
401 
409  std::vector<std::vector<stan::math::var> >
410  decouple_states(const std::vector<std::vector<double> >& y) {
412  using stan::math::var;
413  using std::vector;
414 
415  vector<var> temp_vars(N_);
416  vector<double> temp_gradients(N_);
417  vector<vector<var> > y_return(y.size());
418 
419  for (size_t i = 0; i < y.size(); i++) {
420  // iterate over number of equations
421  for (size_t j = 0; j < N_; j++) {
422  // iterate over parameters for each equation
423  for (size_t k = 0; k < N_; k++)
424  temp_gradients[k] = y[i][y0_.size() + y0_.size() * k + j];
425 
426  temp_vars[j] = precomputed_gradients(y[i][j],
427  y0_, temp_gradients);
428  }
429  y_return[i] = temp_vars;
430  }
431 
432  add_initial_values(y0_, y_return);
433 
434  return y_return;
435  }
436  };
437 
473  template <typename F>
475  const F& f_;
476  const std::vector<stan::math::var>& y0_;
477  std::vector<double> y0_dbl_;
478  const std::vector<stan::math::var>& theta_;
479  std::vector<double> theta_dbl_;
480  const std::vector<double>& x_;
481  const std::vector<int>& x_int_;
482  const size_t N_;
483  const size_t M_;
484  const size_t size_;
485  std::ostream* msgs_;
486 
500  coupled_ode_system(const F& f,
501  const std::vector<stan::math::var>& y0,
502  const std::vector<stan::math::var>& theta,
503  const std::vector<double>& x,
504  const std::vector<int>& x_int,
505  std::ostream* msgs)
506  : f_(f),
507  y0_(y0),
508  y0_dbl_(y0.size(), 0.0),
509  theta_(theta),
510  theta_dbl_(theta.size(), 0.0),
511  x_(x),
512  x_int_(x_int),
513  N_(y0.size()),
514  M_(theta.size()),
515  size_(N_ + N_ * (N_ + M_)),
516  msgs_(msgs) {
517  for (size_t n = 0; n < N_; n++)
518  y0_dbl_[n] = stan::math::value_of(y0[n]);
519 
520  for (size_t m = 0; m < M_; m++)
521  theta_dbl_[m] = stan::math::value_of(theta[m]);
522  }
523 
540  void operator()(const std::vector<double>& z,
541  std::vector<double>& dz_dt,
542  double t) {
543  using std::vector;
544  using stan::math::var;
545 
546  vector<double> y(z.begin(), z.begin() + N_);
547  for (size_t n = 0; n < N_; n++)
548  y[n] += y0_dbl_[n];
549 
550  dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
551  stan::math::check_equal("coupled_ode_system",
552  "dz_dt", dz_dt.size(), N_);
553 
554  vector<double> coupled_sys(N_ * (N_ + M_));
555  vector<double> grad(N_ + M_);
556 
557  try {
559 
560  vector<var> z_vars;
561  z_vars.reserve(N_ + M_);
562 
563  vector<var> y_vars(y.begin(), y.end());
564  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
565 
566  vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
567  z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
568 
569  vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
570 
571  for (size_t i = 0; i < N_; i++) {
573  dy_dt_vars[i].grad(z_vars, grad);
574 
575  for (size_t j = 0; j < N_ + M_; j++) {
576  // orders derivatives by equation (i.e. if there are 2 eqns
577  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
578  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
579  double temp_deriv = grad[j];
580  for (size_t k = 0; k < N_; k++)
581  temp_deriv += z[N_ + N_ * j + k] * grad[k];
582 
583  coupled_sys[i + j * N_] = temp_deriv;
584  }
585  }
586  } catch (const std::exception& e) {
588  throw;
589  }
591 
592  dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
593  }
594 
600  size_t size() const {
601  return size_;
602  }
603 
615  std::vector<double> initial_state() {
616  return std::vector<double>(size_, 0.0);
617  }
618 
626  std::vector<std::vector<stan::math::var> >
627  decouple_states(const std::vector<std::vector<double> >& y) {
628  using std::vector;
629  using stan::math::var;
631 
632  vector<var> vars = y0_;
633  vars.insert(vars.end(), theta_.begin(), theta_.end());
634 
635  vector<var> temp_vars(N_);
636  vector<double> temp_gradients(N_ + M_);
637  vector<vector<var> > y_return(y.size());
638 
639  for (size_t i = 0; i < y.size(); i++) {
640  // iterate over number of equations
641  for (size_t j = 0; j < N_; j++) {
642  // iterate over parameters for each equation
643  for (size_t k = 0; k < N_ + M_; k++)
644  temp_gradients[k] = y[i][N_ + N_ * k + j];
645 
646  temp_vars[j] = precomputed_gradients(y[i][j],
647  vars, temp_gradients);
648  }
649  y_return[i] = temp_vars;
650  }
651  add_initial_values(y0_, y_return);
652  return y_return;
653  }
654  };
655  } // math
656 } // stan
657 
658 #endif
var precomputed_gradients(const double value, const std::vector< var > &operands, const std::vector< double > &gradients)
This function returns a var for an expression that has the specified value, vector of operands...
std::vector< std::vector< stan::math::var > > decouple_states(const std::vector< std::vector< double > > &y)
Returns the base ODE system state corresponding to the specified coupled system state.
std::vector< std::vector< stan::math::var > > decouple_states(const std::vector< std::vector< double > > &y)
Return the solutions to the basic ODE system, including appropriate autodiff partial derivatives...
coupled_ode_system(const F &f, const std::vector< double > &y0, const std::vector< stan::math::var > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs)
Construct a coupled ODE system with the specified base ODE system, base initial state, parameters, data, and a message stream.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:16
static void set_zero_all_adjoints_nested()
Reset all adjoint values in the top nested portion of the stack to zero.
std::vector< double > initial_state()
Returns the initial state of the coupled system.
std::vector< std::vector< stan::math::var > > decouple_states(const std::vector< std::vector< double > > &y)
Return the basic ODE solutions given the specified coupled system solutions, including the partials v...
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:31
static void grad(vari *vi)
Compute the gradient for all variables starting from the specified root variable implementation.
Definition: grad.hpp:30
size_t size() const
Returns the size of the coupled system.
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t)
Calculates the derivative of the coupled ode system with respect to the state y at time t...
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t)
Assign the derivative vector with the system derivatives at the specified state and time...
int M_
size_t size_
Definition: dot_self.hpp:18
std::vector< double > initial_state()
Returns the initial state of the coupled system.
coupled_ode_system(const F &f, const std::vector< stan::math::var > &y0, const std::vector< double > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs)
Construct a coupled ODE system for an unknown initial state and known parameters givne the specified ...
bool check_equal(const char *function, const char *name, const T_y &y, const T_eq &eq)
Return true if y is equal to eq.
Definition: check_equal.hpp:90
void add_initial_values(const std::vector< stan::math::var > &y0, std::vector< std::vector< stan::math::var > > &y)
Increment the state derived from the coupled system in the with the original initial state...
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:95
std::vector< double > initial_state()
Returns the initial state of the coupled system.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
coupled_ode_system(const F &f, const std::vector< stan::math::var > &y0, const std::vector< stan::math::var > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs)
Construct a coupled ODE system with unknown initial value and known parameters, given the base ODE sy...
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t)
Populates the derivative vector with derivatives of the coupled ODE system state with respect to time...
Base template class for a coupled ordinary differential equation system, which adds sensitivities to ...
static void recover_memory_nested()
Recover only the memory used for the top nested call.
size_t size() const
Returns the size of the coupled system.
int N_
size_t size() const
Returns the size of the coupled system.
static void start_nested()
Record the current position so that recover_memory_nested() can find it.

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