| Stan Math Library 2.12.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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/prim/arr/meta/get.hpp>
5 #include <stan/math/prim/arr/meta/length.hpp>
6 #include <stan/math/prim/scal/err/check_equal.hpp>
7 #include <stan/math/prim/arr/functor/coupled_ode_system.hpp>
8 #include <stan/math/rev/scal/fun/value_of_rec.hpp>
9 #include <stan/math/rev/scal/fun/value_of.hpp>
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 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_(theta.size(), 0.0),
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 for (size_t m = 0; m < M_; m++)
103 theta_dbl_[m] = value_of(theta[m]);
104 }
105
123 void operator()(const std::vector<double>& z,
124 std::vector<double>& dz_dt,
125 double t) {
126 using std::vector;
127
128 vector<double> y(z.begin(), z.begin() + N_);
129 dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
130 check_equal("coupled_ode_system",
131 "dz_dt", dz_dt.size(), N_);
132
133 vector<double> coupled_sys(N_ * M_);
134 vector<double> grad(N_ + M_);
135
136 try {
137 start_nested();
138
139 vector<var> z_vars;
140 z_vars.reserve(N_ + M_);
141
142 vector<var> y_vars(y.begin(), y.end());
143 z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
144
145 vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
146 z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
147
148 vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
149
150 for (size_t i = 0; i < N_; i++) {
151 set_zero_all_adjoints_nested();
152 dy_dt_vars[i].grad(z_vars, grad);
153
154 for (size_t j = 0; j < M_; j++) {
155 // orders derivatives by equation (i.e. if there are 2 eqns
156 // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
157 // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
158 double temp_deriv = grad[N_ + j];
159 for (size_t k = 0; k < N_; k++)
160 temp_deriv += z[N_ + N_ * j + k] * grad[k];
161
162 coupled_sys[i + j * N_] = temp_deriv;
163 }
164 }
165 } catch (const std::exception& e) {
166 recover_memory_nested();
167 throw;
168 }
169 recover_memory_nested();
170
171 dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
172 }
173
179 size_t size() const {
180 return size_;
181 }
182
196 std::vector<double> initial_state() {
197 std::vector<double> state(size_, 0.0);
198 for (size_t n = 0; n < N_; n++)
199 state[n] = y0_dbl_[n];
200 return state;
201 }
202
209 std::vector<std::vector<var> >
210 decouple_states(const std::vector<std::vector<double> >& y) {
211 std::vector<var> temp_vars(N_);
212 std::vector<double> temp_gradients(M_);
213 std::vector<std::vector<var> > y_return(y.size());
214
215 for (size_t i = 0; i < y.size(); i++) {
216 // iterate over number of equations
217 for (size_t j = 0; j < N_; j++) {
218 // iterate over parameters for each equation
219 for (size_t k = 0; k < M_; k++)
220 temp_gradients[k] = y[i][y0_dbl_.size() + y0_dbl_.size() * k + j];
221
222 temp_vars[j] = precomputed_gradients(y[i][j],
223 theta_,
224 temp_gradients);
225 }
226 y_return[i] = temp_vars;
227 }
228 return y_return;
229 }
230 };
231
258 template <typename F>
259 struct coupled_ode_system <F, var, double> {
260 const F& f_;
261 const std::vector<var>& y0_;
262 std::vector<double> y0_dbl_;
263 const std::vector<double>& theta_dbl_;
264 const std::vector<double>& x_;
265 const std::vector<int>& x_int_;
266 std::ostream* msgs_;
267 const size_t N_;
268 const size_t M_;
269 const size_t size_;
270
284 coupled_ode_system(const F& f,
285 const std::vector<var>& y0,
286 const std::vector<double>& theta,
287 const std::vector<double>& x,
288 const std::vector<int>& x_int,
289 std::ostream* msgs)
290 : f_(f),
291 y0_(y0),
292 y0_dbl_(y0.size(), 0.0),
293 theta_dbl_(theta),
294 x_(x),
295 x_int_(x_int),
296 msgs_(msgs),
297 N_(y0.size()),
298 M_(theta.size()),
299 size_(N_ + N_ * N_) {
300 for (size_t n = 0; n < N_; n++)
301 y0_dbl_[n] = value_of(y0_[n]);
302 }
303
320 void operator()(const std::vector<double>& z,
321 std::vector<double>& dz_dt,
322 double t) {
323 using std::vector;
324
325 std::vector<double> y(z.begin(), z.begin() + N_);
326 for (size_t n = 0; n < N_; n++)
327 y[n] += y0_dbl_[n];
328
329 dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
330 check_equal("coupled_ode_system",
331 "dz_dt", dz_dt.size(), N_);
332
333 std::vector<double> coupled_sys(N_ * N_);
334 std::vector<double> grad(N_);
335
336 try {
337 start_nested();
338
339 vector<var> z_vars;
340 z_vars.reserve(N_);
341
342 vector<var> y_vars(y.begin(), y.end());
343 z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
344
345 vector<var> dy_dt_vars = f_(t, y_vars, theta_dbl_, x_, x_int_, msgs_);
346
347 for (size_t i = 0; i < N_; i++) {
348 set_zero_all_adjoints_nested();
349 dy_dt_vars[i].grad(z_vars, grad);
350
351 for (size_t j = 0; j < N_; j++) {
352 // orders derivatives by equation (i.e. if there are 2 eqns
353 // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
354 // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
355 double temp_deriv = grad[j];
356 for (size_t k = 0; k < N_; k++)
357 temp_deriv += z[N_ + N_ * j + k] * grad[k];
358
359 coupled_sys[i + j * N_] = temp_deriv;
360 }
361 }
362 } catch (const std::exception& e) {
363 recover_memory_nested();
364 throw;
365 }
366 recover_memory_nested();
367
368 dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
369 }
370
376 size_t size() const {
377 return size_;
378 }
379
394 std::vector<double> initial_state() {
395 return std::vector<double>(size_, 0.0);
396 }
397
405 std::vector<std::vector<var> >
406 decouple_states(const std::vector<std::vector<double> >& y) {
407 using std::vector;
408
409 vector<var> temp_vars(N_);
410 vector<double> temp_gradients(N_);
411 vector<vector<var> > y_return(y.size());
412
413 for (size_t i = 0; i < y.size(); i++) {
414 // iterate over number of equations
415 for (size_t j = 0; j < N_; j++) {
416 // iterate over parameters for each equation
417 for (size_t k = 0; k < N_; k++)
418 temp_gradients[k] = y[i][y0_.size() + y0_.size() * k + j];
419
420 temp_vars[j] = precomputed_gradients(y[i][j],
421 y0_, temp_gradients);
422 }
423 y_return[i] = temp_vars;
424 }
425
426 add_initial_values(y0_, y_return);
427
428 return y_return;
429 }
430 };
431
467 template <typename F>
468 struct coupled_ode_system <F, var, stan::math::var> {
469 const F& f_;
470 const std::vector<var>& y0_;
471 std::vector<double> y0_dbl_;
472 const std::vector<var>& theta_;
473 std::vector<double> theta_dbl_;
474 const std::vector<double>& x_;
475 const std::vector<int>& x_int_;
476 const size_t N_;
477 const size_t M_;
478 const size_t size_;
479 std::ostream* msgs_;
480
494 coupled_ode_system(const F& f,
495 const std::vector<var>& y0,
496 const std::vector<var>& theta,
497 const std::vector<double>& x,
498 const std::vector<int>& x_int,
499 std::ostream* msgs)
500 : f_(f),
501 y0_(y0),
502 y0_dbl_(y0.size(), 0.0),
503 theta_(theta),
504 theta_dbl_(theta.size(), 0.0),
505 x_(x),
506 x_int_(x_int),
507 N_(y0.size()),
508 M_(theta.size()),
509 size_(N_ + N_ * (N_ + M_)),
510 msgs_(msgs) {
511 for (size_t n = 0; n < N_; n++)
512 y0_dbl_[n] = value_of(y0[n]);
513
514 for (size_t m = 0; m < M_; m++)
515 theta_dbl_[m] = value_of(theta[m]);
516 }
517
534 void operator()(const std::vector<double>& z,
535 std::vector<double>& dz_dt,
536 double t) {
537 using std::vector;
538
539 vector<double> y(z.begin(), z.begin() + N_);
540 for (size_t n = 0; n < N_; n++)
541 y[n] += y0_dbl_[n];
542
543 dz_dt = f_(t, y, theta_dbl_, x_, x_int_, msgs_);
544 check_equal("coupled_ode_system",
545 "dz_dt", dz_dt.size(), N_);
546
547 vector<double> coupled_sys(N_ * (N_ + M_));
548 vector<double> grad(N_ + M_);
549
550 try {
551 start_nested();
552
553 vector<var> z_vars;
554 z_vars.reserve(N_ + M_);
555
556 vector<var> y_vars(y.begin(), y.end());
557 z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
558
559 vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
560 z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
561
562 vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
563
564 for (size_t i = 0; i < N_; i++) {
565 set_zero_all_adjoints_nested();
566 dy_dt_vars[i].grad(z_vars, grad);
567
568 for (size_t j = 0; j < N_ + M_; j++) {
569 // orders derivatives by equation (i.e. if there are 2 eqns
570 // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
571 // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
572 double temp_deriv = grad[j];
573 for (size_t k = 0; k < N_; k++)
574 temp_deriv += z[N_ + N_ * j + k] * grad[k];
575
576 coupled_sys[i + j * N_] = temp_deriv;
577 }
578 }
579 } catch (const std::exception& e) {
580 recover_memory_nested();
581 throw;
582 }
583 recover_memory_nested();
584
585 dz_dt.insert(dz_dt.end(), coupled_sys.begin(), coupled_sys.end());
586 }
587
593 size_t size() const {
594 return size_;
595 }
596
608 std::vector<double> initial_state() {
609 return std::vector<double>(size_, 0.0);
610 }
611
619 std::vector<std::vector<var> >
620 decouple_states(const std::vector<std::vector<double> >& y) {
621 using std::vector;
622
623 vector<var> vars = y0_;
624 vars.insert(vars.end(), theta_.begin(), theta_.end());
625
626 vector<var> temp_vars(N_);
627 vector<double> temp_gradients(N_ + M_);
628 vector<vector<var> > y_return(y.size());
629
630 for (size_t i = 0; i < y.size(); i++) {
631 // iterate over number of equations
632 for (size_t j = 0; j < N_; j++) {
633 // iterate over parameters for each equation
634 for (size_t k = 0; k < N_ + M_; k++)
635 temp_gradients[k] = y[i][N_ + N_ * k + j];
636
637 temp_vars[j] = precomputed_gradients(y[i][j],
638 vars, temp_gradients);
639 }
640 y_return[i] = temp_vars;
641 }
642 add_initial_values(y0_, y_return);
643 return y_return;
644 }
645 };
646
647 }
648 }
649 #endif
stan::math::precomputed_gradients
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...
Definition: precomputed_gradients.hpp:95
stan::math::coupled_ode_system< F, double, var >::size_
const size_t size_
Definition: coupled_ode_system.hpp:71
stan::math::coupled_ode_system< F, var, double >::M_
const size_t M_
Definition: coupled_ode_system.hpp:268
stan::math::coupled_ode_system< F, double, var >::msgs_
std::ostream * msgs_
Definition: coupled_ode_system.hpp:72
stan::math::coupled_ode_system< F, var, stan::math::var >::operator()
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...
Definition: coupled_ode_system.hpp:534
stan::math::coupled_ode_system< F, var, stan::math::var >::initial_state
std::vector< double > initial_state()
Returns the initial state of the coupled system.
Definition: coupled_ode_system.hpp:608
stan::math::coupled_ode_system< F, var, double >::size
size_t size() const
Returns the size of the coupled system.
Definition: coupled_ode_system.hpp:376
value_of_rec.hpp
stan::math::coupled_ode_system< F, var, stan::math::var >::y0_dbl_
std::vector< double > y0_dbl_
Definition: coupled_ode_system.hpp:471
stan::math::coupled_ode_system< F, var, stan::math::var >::y0_
const std::vector< var > & y0_
Definition: coupled_ode_system.hpp:470
core.hpp
stan::math::coupled_ode_system< F, var, double >::f_
const F & f_
Definition: coupled_ode_system.hpp:260
stan::math::coupled_ode_system< F, var, stan::math::var >::msgs_
std::ostream * msgs_
Definition: coupled_ode_system.hpp:479
stan::math::coupled_ode_system< F, double, var >::x_
const std::vector< double > & x_
Definition: coupled_ode_system.hpp:67
stan::math::coupled_ode_system< F, double, var >::initial_state
std::vector< double > initial_state()
Returns the initial state of the coupled system.
Definition: coupled_ode_system.hpp:196
stan::math::coupled_ode_system< F, var, double >::y0_dbl_
std::vector< double > y0_dbl_
Definition: coupled_ode_system.hpp:262
stan
Definition: log_sum_exp.hpp:8
stan::math::value_of
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:16
stan::math::coupled_ode_system< F, var, double >::msgs_
std::ostream * msgs_
Definition: coupled_ode_system.hpp:266
stan::math::coupled_ode_system< F, var, stan::math::var >::theta_dbl_
std::vector< double > theta_dbl_
Definition: coupled_ode_system.hpp:473
stan::math::coupled_ode_system< F, var, double >::x_int_
const std::vector< int > & x_int_
Definition: coupled_ode_system.hpp:265
stan::math::set_zero_all_adjoints_nested
static void set_zero_all_adjoints_nested()
Reset all adjoint values in the top nested portion of the stack to zero.
Definition: set_zero_all_adjoints_nested.hpp:17
stan::math::coupled_ode_system< F, double, var >::size
size_t size() const
Returns the size of the coupled system.
Definition: coupled_ode_system.hpp:179
stan::math::coupled_ode_system< F, var, stan::math::var >::size
size_t size() const
Returns the size of the coupled system.
Definition: coupled_ode_system.hpp:593
stan::math::coupled_ode_system< F, var, stan::math::var >::f_
const F & f_
Definition: coupled_ode_system.hpp:469
stan::math::coupled_ode_system< F, var, stan::math::var >::size_
const size_t size_
Definition: coupled_ode_system.hpp:478
stan::math::coupled_ode_system< F, var, double >::y0_
const std::vector< var > & y0_
Definition: coupled_ode_system.hpp:261
stan::math::var
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:30
stan::math::coupled_ode_system< F, var, stan::math::var >::coupled_ode_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 with unknown initial value and known parameters, given the base ODE sy...
Definition: coupled_ode_system.hpp:494
stan::math::grad
static void grad(vari *vi)
Compute the gradient for all variables starting from the specified root variable implementation.
Definition: grad.hpp:30
stan::math::coupled_ode_system< F, var, stan::math::var >::M_
const size_t M_
Definition: coupled_ode_system.hpp:477
stan::math::coupled_ode_system< F, double, var >::coupled_ode_system
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.
Definition: coupled_ode_system.hpp:86
stan::math::coupled_ode_system< F, var, double >::decouple_states
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y)
Return the solutions to the basic ODE system, including appropriate autodiff partial derivatives...
Definition: coupled_ode_system.hpp:406
stan::math::coupled_ode_system< F, double, var >::theta_dbl_
std::vector< double > theta_dbl_
Definition: coupled_ode_system.hpp:66
stan::math::coupled_ode_system< F, var, stan::math::var >::N_
const size_t N_
Definition: coupled_ode_system.hpp:476
M_
int M_
Definition: mdivide_left.hpp:19
size_
size_t size_
Definition: dot_self.hpp:18
stan::math::coupled_ode_system< F, double, var >::decouple_states
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y)
Returns the base ODE system state corresponding to the specified coupled system state.
Definition: coupled_ode_system.hpp:210
check_equal.hpp
stan::math::coupled_ode_system< F, var, stan::math::var >::x_int_
const std::vector< int > & x_int_
Definition: coupled_ode_system.hpp:475
stan::math::coupled_ode_system< F, var, double >::size_
const size_t size_
Definition: coupled_ode_system.hpp:269
stan::math::coupled_ode_system< F, double, var >::y0_dbl_
const std::vector< double > & y0_dbl_
Definition: coupled_ode_system.hpp:64
stan::math::coupled_ode_system< F, var, double >::N_
const size_t N_
Definition: coupled_ode_system.hpp:267
stan::math::coupled_ode_system< F, double, var >::theta_
const std::vector< var > & theta_
Definition: coupled_ode_system.hpp:65
stan::math::check_equal
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:89
stan::math::coupled_ode_system< F, var, double >::theta_dbl_
const std::vector< double > & theta_dbl_
Definition: coupled_ode_system.hpp:263
stan::math::coupled_ode_system< F, double, var >::f_
const F & f_
Definition: coupled_ode_system.hpp:63
stan::math::coupled_ode_system< F, double, var >::N_
const size_t N_
Definition: coupled_ode_system.hpp:69
stan::math::e
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:94
stan::math::coupled_ode_system< F, var, double >::initial_state
std::vector< double > initial_state()
Returns the initial state of the coupled system.
Definition: coupled_ode_system.hpp:394
stan::math::coupled_ode_system< F, double, var >::x_int_
const std::vector< int > & x_int_
Definition: coupled_ode_system.hpp:68
stan::math::size
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
stan::math::coupled_ode_system
Base template class for a coupled ordinary differential equation system, which adds sensitivities to ...
Definition: coupled_ode_system.hpp:25
stan::math::recover_memory_nested
static void recover_memory_nested()
Recover only the memory used for the top nested call.
Definition: recover_memory_nested.hpp:20
coupled_ode_system.hpp
stan::math::coupled_ode_system< F, var, double >::coupled_ode_system
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 ...
Definition: coupled_ode_system.hpp:284
get.hpp
stan::math::coupled_ode_system< F, double, var >::M_
const size_t M_
Definition: coupled_ode_system.hpp:70
stan::math::coupled_ode_system< F, var, double >::operator()
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...
Definition: coupled_ode_system.hpp:320
stan::math::coupled_ode_system< F, var, double >::x_
const std::vector< double > & x_
Definition: coupled_ode_system.hpp:264
stan::math::add_initial_values
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...
Definition: coupled_ode_system.hpp:33
stan::math::coupled_ode_system< F, var, stan::math::var >::decouple_states
std::vector< std::vector< 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...
Definition: coupled_ode_system.hpp:620
N_
int N_
Definition: mdivide_left.hpp:20
stan::math::coupled_ode_system< F, var, stan::math::var >::x_
const std::vector< double > & x_
Definition: coupled_ode_system.hpp:474
value_of.hpp
stan::math::start_nested
static void start_nested()
Record the current position so that recover_memory_nested() can find it.
Definition: start_nested.hpp:13
stan::math::coupled_ode_system< F, double, var >::operator()
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...
Definition: coupled_ode_system.hpp:123
stan::math::coupled_ode_system< F, var, stan::math::var >::theta_
const std::vector< var > & theta_
Definition: coupled_ode_system.hpp:472
length.hpp

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