Stan Math Library  2.15.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 
33  inline void add_initial_values(const std::vector<var>& y0,
34  std::vector<std::vector<var> >& y) {
35  for (size_t n = 0; n < y.size(); n++)
36  for (size_t m = 0; m < y0.size(); m++)
37  y[n][m] += y0[m];
38  }
39 
61  template <typename F>
62  struct coupled_ode_system <F, double, var> {
63  const F& f_;
64  const std::vector<double>& y0_dbl_;
65  const std::vector<var>& theta_;
66  const std::vector<double> theta_dbl_;
67  const std::vector<double>& x_;
68  const std::vector<int>& x_int_;
69  const size_t N_;
70  const size_t M_;
71  const size_t size_;
72  std::ostream* msgs_;
73 
86  coupled_ode_system(const F& f,
87  const std::vector<double>& y0,
88  const std::vector<var>& theta,
89  const std::vector<double>& x,
90  const std::vector<int>& x_int,
91  std::ostream* msgs)
92  : f_(f),
93  y0_dbl_(y0),
94  theta_(theta),
95  theta_dbl_(value_of(theta)),
96  x_(x),
97  x_int_(x_int),
98  N_(y0.size()),
99  M_(theta.size()),
100  size_(N_ + N_ * M_),
101  msgs_(msgs) {}
102 
120  void operator()(const std::vector<double>& z,
121  std::vector<double>& dz_dt,
122  double t) const {
123  using std::vector;
124 
125  vector<double> grad(N_ + M_);
126 
127  try {
128  start_nested();
129 
130  vector<var> z_vars;
131  z_vars.reserve(N_ + M_);
132 
133  vector<var> y_vars(z.begin(), z.begin() + N_);
134  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
135 
136  vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
137  z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
138 
139  vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
140 
141  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
142  "states", N_);
143 
144  for (size_t i = 0; i < N_; i++) {
145  dz_dt[i] = dy_dt_vars[i].val();
147  dy_dt_vars[i].grad(z_vars, grad);
148 
149  for (size_t j = 0; j < M_; j++) {
150  // orders derivatives by equation (i.e. if there are 2 eqns
151  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
152  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
153  double temp_deriv = grad[N_ + j];
154  for (size_t k = 0; k < N_; k++)
155  temp_deriv += z[N_ + N_ * j + k] * grad[k];
156 
157  dz_dt[N_ + i + j * N_] = temp_deriv;
158  }
159  }
160  } catch (const std::exception& e) {
162  throw;
163  }
165  }
166 
172  size_t size() const {
173  return size_;
174  }
175 
189  std::vector<double> initial_state() const {
190  std::vector<double> state(size_, 0.0);
191  for (size_t n = 0; n < N_; n++)
192  state[n] = y0_dbl_[n];
193  return state;
194  }
195 
202  std::vector<std::vector<var> >
203  decouple_states(const std::vector<std::vector<double> >& y) const {
204  std::vector<var> temp_vars(N_);
205  std::vector<double> temp_gradients(M_);
206  std::vector<std::vector<var> > y_return(y.size());
207 
208  for (size_t i = 0; i < y.size(); i++) {
209  // iterate over number of equations
210  for (size_t j = 0; j < N_; j++) {
211  // iterate over parameters for each equation
212  for (size_t k = 0; k < M_; k++)
213  temp_gradients[k] = y[i][y0_dbl_.size() + y0_dbl_.size() * k + j];
214 
215  temp_vars[j] = precomputed_gradients(y[i][j],
216  theta_,
217  temp_gradients);
218  }
219  y_return[i] = temp_vars;
220  }
221  return y_return;
222  }
223  };
224 
251  template <typename F>
252  struct coupled_ode_system <F, var, double> {
253  const F& f_;
254  const std::vector<var>& y0_;
255  const std::vector<double> y0_dbl_;
256  const std::vector<double>& theta_dbl_;
257  const std::vector<double>& x_;
258  const std::vector<int>& x_int_;
259  std::ostream* msgs_;
260  const size_t N_;
261  const size_t M_;
262  const size_t size_;
263 
277  coupled_ode_system(const F& f,
278  const std::vector<var>& y0,
279  const std::vector<double>& theta,
280  const std::vector<double>& x,
281  const std::vector<int>& x_int,
282  std::ostream* msgs)
283  : f_(f),
284  y0_(y0),
285  y0_dbl_(value_of(y0)),
286  theta_dbl_(theta),
287  x_(x),
288  x_int_(x_int),
289  msgs_(msgs),
290  N_(y0.size()),
291  M_(theta.size()),
292  size_(N_ + N_ * N_) {}
293 
310  void operator()(const std::vector<double>& z,
311  std::vector<double>& dz_dt,
312  double t) const {
313  using std::vector;
314 
315  std::vector<double> y(z.begin(), z.begin() + N_);
316  for (size_t n = 0; n < N_; n++)
317  y[n] += y0_dbl_[n];
318 
319  std::vector<double> grad(N_);
320 
321  try {
322  start_nested();
323 
324  vector<var> z_vars;
325  z_vars.reserve(N_);
326 
327  vector<var> y_vars(y.begin(), y.end());
328  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
329 
330  vector<var> dy_dt_vars = f_(t, y_vars, theta_dbl_, x_, x_int_, msgs_);
331 
332  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
333  "states", N_);
334 
335  for (size_t i = 0; i < N_; i++) {
336  dz_dt[i] = dy_dt_vars[i].val();
338  dy_dt_vars[i].grad(z_vars, grad);
339 
340  for (size_t j = 0; j < N_; j++) {
341  // orders derivatives by equation (i.e. if there are 2 eqns
342  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
343  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
344  double temp_deriv = grad[j];
345  for (size_t k = 0; k < N_; k++)
346  temp_deriv += z[N_ + N_ * j + k] * grad[k];
347 
348  dz_dt[N_ + i + j * N_] = temp_deriv;
349  }
350  }
351  } catch (const std::exception& e) {
353  throw;
354  }
356  }
357 
363  size_t size() const {
364  return size_;
365  }
366 
381  std::vector<double> initial_state() const {
382  return std::vector<double>(size_, 0.0);
383  }
384 
392  std::vector<std::vector<var> >
393  decouple_states(const std::vector<std::vector<double> >& y) const {
394  using std::vector;
395 
396  vector<var> temp_vars(N_);
397  vector<double> temp_gradients(N_);
398  vector<vector<var> > y_return(y.size());
399 
400  for (size_t i = 0; i < y.size(); i++) {
401  // iterate over number of equations
402  for (size_t j = 0; j < N_; j++) {
403  // iterate over parameters for each equation
404  for (size_t k = 0; k < N_; k++)
405  temp_gradients[k] = y[i][y0_.size() + y0_.size() * k + j];
406 
407  temp_vars[j] = precomputed_gradients(y[i][j],
408  y0_, temp_gradients);
409  }
410  y_return[i] = temp_vars;
411  }
412 
413  add_initial_values(y0_, y_return);
414 
415  return y_return;
416  }
417  };
418 
454  template <typename F>
456  const F& f_;
457  const std::vector<var>& y0_;
458  const std::vector<double> y0_dbl_;
459  const std::vector<var>& theta_;
460  const std::vector<double> theta_dbl_;
461  const std::vector<double>& x_;
462  const std::vector<int>& x_int_;
463  const size_t N_;
464  const size_t M_;
465  const size_t size_;
466  std::ostream* msgs_;
467 
481  coupled_ode_system(const F& f,
482  const std::vector<var>& y0,
483  const std::vector<var>& theta,
484  const std::vector<double>& x,
485  const std::vector<int>& x_int,
486  std::ostream* msgs)
487  : f_(f),
488  y0_(y0),
489  y0_dbl_(value_of(y0)),
490  theta_(theta),
491  theta_dbl_(value_of(theta)),
492  x_(x),
493  x_int_(x_int),
494  N_(y0.size()),
495  M_(theta.size()),
496  size_(N_ + N_ * (N_ + M_)),
497  msgs_(msgs) {}
498 
515  void operator()(const std::vector<double>& z,
516  std::vector<double>& dz_dt,
517  double t) const {
518  using std::vector;
519 
520  vector<double> y(z.begin(), z.begin() + N_);
521  for (size_t n = 0; n < N_; n++)
522  y[n] += y0_dbl_[n];
523 
524  vector<double> grad(N_ + M_);
525 
526  try {
527  start_nested();
528 
529  vector<var> z_vars;
530  z_vars.reserve(N_ + M_);
531 
532  vector<var> y_vars(y.begin(), y.end());
533  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
534 
535  vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
536  z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
537 
538  vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
539 
540  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
541  "states", N_);
542 
543  for (size_t i = 0; i < N_; i++) {
544  dz_dt[i] = dy_dt_vars[i].val();
546  dy_dt_vars[i].grad(z_vars, grad);
547 
548  for (size_t j = 0; j < N_ + M_; j++) {
549  // orders derivatives by equation (i.e. if there are 2 eqns
550  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
551  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
552  double temp_deriv = grad[j];
553  for (size_t k = 0; k < N_; k++)
554  temp_deriv += z[N_ + N_ * j + k] * grad[k];
555 
556  dz_dt[N_ + i + j * N_] = temp_deriv;
557  }
558  }
559  } catch (const std::exception& e) {
561  throw;
562  }
564  }
565 
571  size_t size() const {
572  return size_;
573  }
574 
586  std::vector<double> initial_state() const {
587  return std::vector<double>(size_, 0.0);
588  }
589 
597  std::vector<std::vector<var> >
598  decouple_states(const std::vector<std::vector<double> >& y) const {
599  using std::vector;
600 
601  vector<var> vars = y0_;
602  vars.insert(vars.end(), theta_.begin(), theta_.end());
603 
604  vector<var> temp_vars(N_);
605  vector<double> temp_gradients(N_ + M_);
606  vector<vector<var> > y_return(y.size());
607 
608  for (size_t i = 0; i < y.size(); i++) {
609  // iterate over number of equations
610  for (size_t j = 0; j < N_; j++) {
611  // iterate over parameters for each equation
612  for (size_t k = 0; k < N_ + M_; k++)
613  temp_gradients[k] = y[i][N_ + N_ * k + j];
614 
615  temp_vars[j] = precomputed_gradients(y[i][j],
616  vars, temp_gradients);
617  }
618  y_return[i] = temp_vars;
619  }
620  add_initial_values(y0_, y_return);
621  return y_return;
622  }
623  };
624 
625  }
626 }
627 #endif
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
Return the basic ODE solutions given the specified coupled system solutions, including the partials v...
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:16
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
static void set_zero_all_adjoints_nested()
Reset all adjoint values in the top nested portion of the stack to zero.
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:30
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
coupled_ode_system(const F &f, const std::vector< var > &y0, const std::vector< 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...
static void grad(vari *vi)
Compute the gradient for all variables starting from the specified root variable implementation.
Definition: grad.hpp:30
coupled_ode_system(const F &f, const std::vector< double > &y0, const std::vector< 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.
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
Return the solutions to the basic ODE system, including appropriate autodiff partial derivatives...
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
Populates the derivative vector with derivatives of the coupled ODE system state with respect to time...
var precomputed_gradients(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...
int M_
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
Assign the derivative vector with the system derivatives at the specified state and time...
size_t size_
Definition: dot_self.hpp:18
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
Calculates the derivative of the coupled ode system with respect to the state y at time t...
size_t size() const
Returns the size of the coupled system.
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:94
size_t size() const
Returns the size of the coupled system.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
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.
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
Returns the base ODE system state corresponding to the specified coupled system state.
coupled_ode_system(const F &f, const std::vector< 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 ...
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
void add_initial_values(const std::vector< var > &y0, std::vector< std::vector< var > > &y)
Increment the state derived from the coupled system in the with the original initial state...
int N_
static void start_nested()
Record the current position so that recover_memory_nested() can find it.
size_t size() const
Returns the size of the coupled system.

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