Stan Math Library  2.20.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 
4 #include <stan/math/rev/meta.hpp>
9 #include <stan/math/rev/core.hpp>
10 #include <ostream>
11 #include <stdexcept>
12 #include <vector>
13 
14 namespace stan {
15 namespace math {
16 
58 template <typename F>
59 struct coupled_ode_system<F, double, var> {
60  const F& f_;
61  const std::vector<double>& y0_dbl_;
62  const std::vector<var>& theta_;
63  std::vector<var> theta_nochain_;
64  const std::vector<double>& x_;
65  const std::vector<int>& x_int_;
66  const size_t N_;
67  const size_t M_;
68  const size_t size_;
69  std::ostream* msgs_;
70 
83  coupled_ode_system(const F& f, const std::vector<double>& y0,
84  const std::vector<var>& theta,
85  const std::vector<double>& x,
86  const std::vector<int>& x_int, std::ostream* msgs)
87  : f_(f),
88  y0_dbl_(y0),
89  theta_(theta),
90  x_(x),
91  x_int_(x_int),
92  N_(y0.size()),
93  M_(theta.size()),
94  size_(N_ + N_ * M_),
95  msgs_(msgs) {
96  for (const var& p : theta)
97  theta_nochain_.emplace_back(var(new vari(p.val(), false)));
98  }
99 
114  void operator()(const std::vector<double>& z, std::vector<double>& dz_dt,
115  double t) const {
116  using std::vector;
117 
118  try {
119  start_nested();
120 
121  const vector<var> y_vars(z.begin(), z.begin() + N_);
122 
123  vector<var> dy_dt_vars = f_(t, y_vars, theta_nochain_, x_, x_int_, msgs_);
124 
125  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
126  "states", N_);
127 
128  for (size_t i = 0; i < N_; i++) {
129  dz_dt[i] = dy_dt_vars[i].val();
130  dy_dt_vars[i].grad();
131 
132  for (size_t j = 0; j < M_; j++) {
133  // orders derivatives by equation (i.e. if there are 2 eqns
134  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
135  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
136  double temp_deriv = theta_nochain_[j].adj();
137  const size_t offset = N_ + N_ * j;
138  for (size_t k = 0; k < N_; k++)
139  temp_deriv += z[offset + k] * y_vars[k].adj();
140 
141  dz_dt[offset + i] = temp_deriv;
142  }
143 
145  // Parameters stored on the outer (non-nested) nochain stack are not
146  // reset to zero by the last call. This is done as a separate step here.
147  // See efficiency note above on template specalization for more details
148  // on this.
149  for (size_t j = 0; j < M_; ++j)
150  theta_nochain_[j].vi_->set_zero_adjoint();
151  }
152  } catch (const std::exception& e) {
154  throw;
155  }
157  }
158 
164  size_t size() const { return size_; }
165 
181  std::vector<double> initial_state() const {
182  std::vector<double> state(size_, 0.0);
183  for (size_t n = 0; n < N_; n++)
184  state[n] = y0_dbl_[n];
185  return state;
186  }
187 };
188 
216 template <typename F>
217 struct coupled_ode_system<F, var, double> {
218  const F& f_;
219  const std::vector<var>& y0_;
220  const std::vector<double>& theta_dbl_;
221  const std::vector<double>& x_;
222  const std::vector<int>& x_int_;
223  std::ostream* msgs_;
224  const size_t N_;
225  const size_t M_;
226  const size_t size_;
227 
240  coupled_ode_system(const F& f, const std::vector<var>& y0,
241  const std::vector<double>& theta,
242  const std::vector<double>& x,
243  const std::vector<int>& x_int, std::ostream* msgs)
244  : f_(f),
245  y0_(y0),
246  theta_dbl_(theta),
247  x_(x),
248  x_int_(x_int),
249  msgs_(msgs),
250  N_(y0.size()),
251  M_(theta.size()),
252  size_(N_ + N_ * N_) {}
253 
268  void operator()(const std::vector<double>& z, std::vector<double>& dz_dt,
269  double t) const {
270  using std::vector;
271 
272  try {
273  start_nested();
274 
275  const vector<var> y_vars(z.begin(), z.begin() + N_);
276 
277  vector<var> dy_dt_vars = f_(t, y_vars, theta_dbl_, x_, x_int_, msgs_);
278 
279  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
280  "states", N_);
281 
282  for (size_t i = 0; i < N_; i++) {
283  dz_dt[i] = dy_dt_vars[i].val();
284  dy_dt_vars[i].grad();
285 
286  for (size_t j = 0; j < N_; j++) {
287  // orders derivatives by equation (i.e. if there are 2 eqns
288  // (y1, y2) and 2 initial conditions (y0_a, y0_b), dy_dt will be
289  // ordered as: dy1_dt, dy2_dt, dy1_d{y0_a}, dy2_d{y0_a}, dy1_d{y0_b},
290  // dy2_d{y0_b}
291  double temp_deriv = 0;
292  const size_t offset = N_ + N_ * j;
293  for (size_t k = 0; k < N_; k++)
294  temp_deriv += z[offset + k] * y_vars[k].adj();
295 
296  dz_dt[offset + i] = temp_deriv;
297  }
298 
300  }
301  } catch (const std::exception& e) {
303  throw;
304  }
306  }
307 
313  size_t size() const { return size_; }
314 
328  std::vector<double> initial_state() const {
329  std::vector<double> initial(size_, 0.0);
330  for (size_t i = 0; i < N_; i++)
331  initial[i] = value_of(y0_[i]);
332  for (size_t i = 0; i < N_; i++)
333  initial[N_ + i * N_ + i] = 1.0;
334  return initial;
335  }
336 };
337 
386 template <typename F>
388  const F& f_;
389  const std::vector<var>& y0_;
390  const std::vector<var>& theta_;
391  std::vector<var> theta_nochain_;
392  const std::vector<double>& x_;
393  const std::vector<int>& x_int_;
394  const size_t N_;
395  const size_t M_;
396  const size_t size_;
397  std::ostream* msgs_;
398 
411  coupled_ode_system(const F& f, const std::vector<var>& y0,
412  const std::vector<var>& theta,
413  const std::vector<double>& x,
414  const std::vector<int>& x_int, std::ostream* msgs)
415  : f_(f),
416  y0_(y0),
417  theta_(theta),
418  x_(x),
419  x_int_(x_int),
420  N_(y0.size()),
421  M_(theta.size()),
422  size_(N_ + N_ * (N_ + M_)),
423  msgs_(msgs) {
424  for (const var& p : theta)
425  theta_nochain_.emplace_back(var(new vari(p.val(), false)));
426  }
427 
442  void operator()(const std::vector<double>& z, std::vector<double>& dz_dt,
443  double t) const {
444  using std::vector;
445 
446  try {
447  start_nested();
448 
449  const vector<var> y_vars(z.begin(), z.begin() + N_);
450 
451  vector<var> dy_dt_vars = f_(t, y_vars, theta_nochain_, x_, x_int_, msgs_);
452 
453  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
454  "states", N_);
455 
456  for (size_t i = 0; i < N_; i++) {
457  dz_dt[i] = dy_dt_vars[i].val();
458  dy_dt_vars[i].grad();
459 
460  for (size_t j = 0; j < N_; j++) {
461  // orders derivatives by equation (i.e. if there are 2 eqns
462  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
463  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
464  double temp_deriv = 0;
465  const size_t offset = N_ + N_ * j;
466  for (size_t k = 0; k < N_; k++)
467  temp_deriv += z[offset + k] * y_vars[k].adj();
468 
469  dz_dt[offset + i] = temp_deriv;
470  }
471 
472  for (size_t j = 0; j < M_; j++) {
473  double temp_deriv = theta_nochain_[j].adj();
474  const size_t offset = N_ + N_ * N_ + N_ * j;
475  for (size_t k = 0; k < N_; k++)
476  temp_deriv += z[offset + k] * y_vars[k].adj();
477 
478  dz_dt[offset + i] = temp_deriv;
479  }
480 
482  // Parameters stored on the outer (non-nested) nochain stack are not
483  // reset to zero by the last call. This is done as a separate step here.
484  // See efficiency note above on template specalization for more details
485  // on this.
486  for (size_t j = 0; j < M_; ++j)
487  theta_nochain_[j].vi_->set_zero_adjoint();
488  }
489  } catch (const std::exception& e) {
491  throw;
492  }
494  }
495 
501  size_t size() const { return size_; }
502 
521  std::vector<double> initial_state() const {
522  std::vector<double> initial(size_, 0.0);
523  for (size_t i = 0; i < N_; i++)
524  initial[i] = value_of(y0_[i]);
525  for (size_t i = 0; i < N_; i++)
526  initial[N_ + i * N_ + i] = 1.0;
527  return initial;
528  }
529 };
530 
531 } // namespace math
532 } // namespace stan
533 #endif
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
static void set_zero_all_adjoints_nested()
Reset all adjoint values in the top nested portion of the stack to zero.
The variable implementation base class.
Definition: vari.hpp:30
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:33
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< 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 from the base system function, initial state of the base system...
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
Calculates the derivative of the coupled ode system with respect to time.
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 time.
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 time.
size_t size() const
Returns the size of the coupled system.
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:87
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
The coupled_ode_system represents the coupled ode system, which is the base ode and the sensitivities...
static void recover_memory_nested()
Recover only the memory used for the top nested call.
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 from the base system function, initial state of the base system...
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
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 from the base system function, initial state of the base system...
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–2018, Stan Development Team.