Stan Math Library  2.9.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 
8 
11 #include <stan/math/rev/core.hpp>
14 #include <ostream>
15 #include <vector>
16 
17 namespace stan {
18 
19  namespace math {
20 
21  // This code is in this directory because it includes var
22  // It is in namespace stan::math so that the partial template
23  // specializations are treated as such.
24 
25 
37  void add_initial_values(const std::vector<stan::math::var>& y0,
38  std::vector<std::vector<stan::math::var> >& y) {
39  for (size_t n = 0; n < y.size(); n++)
40  for (size_t m = 0; m < y0.size(); m++)
41  y[n][m] += y0[m];
42  }
43 
65  template <typename F>
66  struct coupled_ode_system <F, double, stan::math::var> {
67  const F& f_;
68  const std::vector<double>& y0_dbl_;
69  const std::vector<stan::math::var>& theta_;
70  std::vector<double> theta_dbl_;
71  const std::vector<double>& x_;
72  const std::vector<int>& x_int_;
73  const size_t N_;
74  const size_t M_;
75  const size_t size_;
76  std::ostream* msgs_;
77 
90  coupled_ode_system(const F& f,
91  const std::vector<double>& y0,
92  const std::vector<stan::math::var>& theta,
93  const std::vector<double>& x,
94  const std::vector<int>& x_int,
95  std::ostream* msgs)
96  : f_(f),
97  y0_dbl_(y0),
98  theta_(theta),
99  theta_dbl_(theta.size(), 0.0),
100  x_(x),
101  x_int_(x_int),
102  N_(y0.size()),
103  M_(theta.size()),
104  size_(N_ + N_ * M_),
105  msgs_(msgs) {
106  for (size_t m = 0; m < M_; m++)
107  theta_dbl_[m] = stan::math::value_of(theta[m]);
108  }
109 
127  void operator()(const std::vector<double>& z,
128  std::vector<double>& dz_dt,
129  double t) {
130  using std::vector;
131  using stan::math::var;
132 
133  vector<double> y(z.begin(), z.begin() + N_);
134  dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
135  stan::math::check_equal("coupled_ode_system",
136  "dz_dt", dz_dt.size(), N_);
137 
138  vector<double> coupled_sys(N_ * M_);
139  vector<double> grad(N_ + M_);
140 
141  try {
143 
144  vector<var> z_vars;
145  z_vars.reserve(N_ + M_);
146 
147  vector<var> y_vars(y.begin(), y.end());
148  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
149 
150  vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
151  z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
152 
153  vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
154 
155  for (size_t i = 0; i < N_; i++) {
157  dy_dt_vars[i].grad(z_vars, grad);
158 
159  for (size_t j = 0; j < M_; j++) {
160  // orders derivatives by equation (i.e. if there are 2 eqns
161  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
162  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
163  double temp_deriv = grad[N_ + j];
164  for (size_t k = 0; k < N_; k++)
165  temp_deriv += z[N_ + N_ * j + k] * grad[k];
166 
167  coupled_sys[i + j * N_] = temp_deriv;
168  }
169  }
170  } catch (const std::exception& e) {
172  throw;
173  }
175 
176  dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
177  }
178 
184  size_t size() const {
185  return size_;
186  }
187 
201  std::vector<double> initial_state() {
202  std::vector<double> state(size_, 0.0);
203  for (size_t n = 0; n < N_; n++)
204  state[n] = y0_dbl_[n];
205  return state;
206  }
207 
214  std::vector<std::vector<stan::math::var> >
215  decouple_states(const std::vector<std::vector<double> >& y) {
217  std::vector<stan::math::var> temp_vars;
218  std::vector<double> temp_gradients;
219  std::vector<std::vector<stan::math::var> > y_return(y.size());
220 
221  for (size_t i = 0; i < y.size(); i++) {
222  temp_vars.clear();
223 
224  // iterate over number of equations
225  for (size_t j = 0; j < N_; j++) {
226  temp_gradients.clear();
227 
228  // iterate over parameters for each equation
229  for (size_t k = 0; k < M_; k++)
230  temp_gradients.push_back(y[i][y0_dbl_.size()
231  + y0_dbl_.size() * k + j]);
232 
233  temp_vars.push_back(precomputed_gradients(y[i][j],
234  theta_,
235  temp_gradients));
236  }
237  y_return[i] = temp_vars;
238  }
239 
240  return y_return;
241  }
242  };
243 
270  template <typename F>
271  struct coupled_ode_system <F, stan::math::var, double> {
272  const F& f_;
273  const std::vector<stan::math::var>& y0_;
274  std::vector<double> y0_dbl_;
275  const std::vector<double>& theta_dbl_;
276  const std::vector<double>& x_;
277  const std::vector<int>& x_int_;
278  std::ostream* msgs_;
279  const size_t N_;
280  const size_t M_;
281  const size_t size_;
282 
296  coupled_ode_system(const F& f,
297  const std::vector<stan::math::var>& y0,
298  const std::vector<double>& theta,
299  const std::vector<double>& x,
300  const std::vector<int>& x_int,
301  std::ostream* msgs)
302  : f_(f),
303  y0_(y0),
304  y0_dbl_(y0.size(), 0.0),
305  theta_dbl_(theta),
306  x_(x),
307  x_int_(x_int),
308  msgs_(msgs),
309  N_(y0.size()),
310  M_(theta.size()),
311  size_(N_ + N_ * N_) {
312  for (size_t n = 0; n < N_; n++)
313  y0_dbl_[n] = stan::math::value_of(y0_[n]);
314  }
315 
332  void operator()(const std::vector<double>& z,
333  std::vector<double>& dz_dt,
334  double t) {
335  using std::vector;
336  using stan::math::var;
337 
338  std::vector<double> y(z.begin(), z.begin() + N_);
339  for (size_t n = 0; n < N_; n++)
340  y[n] += y0_dbl_[n];
341 
342  dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
343  stan::math::check_equal("coupled_ode_system",
344  "dz_dt", dz_dt.size(), N_);
345 
346  std::vector<double> coupled_sys(N_ * N_);
347  std::vector<double> grad(N_);
348 
349  try {
351 
352  vector<var> z_vars;
353  z_vars.reserve(N_);
354 
355  vector<var> y_vars(y.begin(), y.end());
356  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
357 
358  vector<var> dy_dt_vars = f_(t, y_vars, theta_dbl_, x_, x_int_, msgs_);
359 
360  for (size_t i = 0; i < N_; i++) {
362  dy_dt_vars[i].grad(z_vars, grad);
363 
364  for (size_t j = 0; j < N_; j++) {
365  // orders derivatives by equation (i.e. if there are 2 eqns
366  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
367  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
368  double temp_deriv = grad[j];
369  for (size_t k = 0; k < N_; k++)
370  temp_deriv += z[N_ + N_ * j + k] * grad[k];
371 
372  coupled_sys[i + j * N_] = temp_deriv;
373  }
374  }
375  } catch (const std::exception& e) {
377  throw;
378  }
380 
381  dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
382  }
383 
389  size_t size() const {
390  return size_;
391  }
392 
407  std::vector<double> initial_state() {
408  return std::vector<double>(size_, 0.0);
409  }
410 
418  std::vector<std::vector<stan::math::var> >
419  decouple_states(const std::vector<std::vector<double> >& y) {
421  using stan::math::var;
422  using std::vector;
423 
424  vector<var> temp_vars;
425  vector<double> temp_gradients;
426  vector<vector<var> > y_return(y.size());
427 
428  for (size_t i = 0; i < y.size(); i++) {
429  temp_vars.clear();
430 
431  // iterate over number of equations
432  for (size_t j = 0; j < N_; j++) {
433  temp_gradients.clear();
434 
435  // iterate over parameters for each equation
436  for (size_t k = 0; k < N_; k++)
437  temp_gradients.push_back(y[i][y0_.size() + y0_.size() * k + j]);
438 
439  temp_vars.push_back(precomputed_gradients(y[i][j],
440  y0_, temp_gradients));
441  }
442 
443  y_return[i] = temp_vars;
444  }
445 
446  add_initial_values(y0_, y_return);
447 
448  return y_return;
449  }
450  };
451 
487  template <typename F>
489  const F& f_;
490  const std::vector<stan::math::var>& y0_;
491  std::vector<double> y0_dbl_;
492  const std::vector<stan::math::var>& theta_;
493  std::vector<double> theta_dbl_;
494  const std::vector<double>& x_;
495  const std::vector<int>& x_int_;
496  const size_t N_;
497  const size_t M_;
498  const size_t size_;
499  std::ostream* msgs_;
500 
514  coupled_ode_system(const F& f,
515  const std::vector<stan::math::var>& y0,
516  const std::vector<stan::math::var>& theta,
517  const std::vector<double>& x,
518  const std::vector<int>& x_int,
519  std::ostream* msgs)
520  : f_(f),
521  y0_(y0),
522  y0_dbl_(y0.size(), 0.0),
523  theta_(theta),
524  theta_dbl_(theta.size(), 0.0),
525  x_(x),
526  x_int_(x_int),
527  N_(y0.size()),
528  M_(theta.size()),
529  size_(N_ + N_ * (N_ + M_)),
530  msgs_(msgs) {
531  for (size_t n = 0; n < N_; n++)
532  y0_dbl_[n] = stan::math::value_of(y0[n]);
533 
534  for (size_t m = 0; m < M_; m++)
535  theta_dbl_[m] = stan::math::value_of(theta[m]);
536  }
537 
554  void operator()(const std::vector<double>& z,
555  std::vector<double>& dz_dt,
556  double t) {
557  using std::vector;
558  using stan::math::var;
559 
560  vector<double> y(z.begin(), z.begin() + N_);
561  for (size_t n = 0; n < N_; n++)
562  y[n] += y0_dbl_[n];
563 
564  dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
565  stan::math::check_equal("coupled_ode_system",
566  "dz_dt", dz_dt.size(), N_);
567 
568  vector<double> coupled_sys(N_ * (N_ + M_));
569  vector<double> grad(N_ + M_);
570 
571  try {
573 
574  vector<var> z_vars;
575  z_vars.reserve(N_ + M_);
576 
577  vector<var> y_vars(y.begin(), y.end());
578  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
579 
580  vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
581  z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
582 
583  vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
584 
585  for (size_t i = 0; i < N_; i++) {
587  dy_dt_vars[i].grad(z_vars, grad);
588 
589  for (size_t j = 0; j < N_ + M_; j++) {
590  // orders derivatives by equation (i.e. if there are 2 eqns
591  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
592  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
593  double temp_deriv = grad[j];
594  for (size_t k = 0; k < N_; k++)
595  temp_deriv += z[N_ + N_ * j + k] * grad[k];
596 
597  coupled_sys[i + j * N_] = temp_deriv;
598  }
599  }
600  } catch (const std::exception& e) {
602  throw;
603  }
605 
606  dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
607  }
608 
614  size_t size() const {
615  return size_;
616  }
617 
629  std::vector<double> initial_state() {
630  return std::vector<double>(size_, 0.0);
631  }
632 
640  std::vector<std::vector<stan::math::var> >
641  decouple_states(const std::vector<std::vector<double> >& y) {
642  using std::vector;
643  using stan::math::var;
645 
646  vector<var> vars = y0_;
647  vars.insert(vars.end(), theta_.begin(), theta_.end());
648 
649  vector<var> temp_vars;
650  vector<double> temp_gradients;
651  vector<vector<var> > y_return(y.size());
652 
653  for (size_t i = 0; i < y.size(); i++) {
654  temp_vars.clear();
655 
656  // iterate over number of equations
657  for (size_t j = 0; j < N_; j++) {
658  temp_gradients.clear();
659 
660  // iterate over parameters for each equation
661  for (size_t k = 0; k < N_ + M_; k++)
662  temp_gradients.push_back(y[i][N_ + N_ * k + j]);
663 
664  temp_vars.push_back(precomputed_gradients(y[i][j],
665  vars, temp_gradients));
666  }
667  y_return[i] = temp_vars;
668  }
669  add_initial_values(y0_, y_return);
670  return y_return;
671  }
672  };
673 
674 
675  }
676 
677 }
678 
679 #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–2015, Stan Development Team.