Stan Math Library  2.6.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros
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 
5 #include <stan/math/rev/core.hpp>
8 #include <ostream>
9 #include <vector>
10 
11 namespace stan {
12 
13  namespace math {
14 
15  // This code is in this directory because it includes var
16  // It is in namespace stan::math so that the partial template
17  // specializations are treated as such.
18 
19 
31  void add_initial_values(const std::vector<stan::math::var>& y0,
32  std::vector<std::vector<stan::math::var> >& y) {
33  for (size_t n = 0; n < y.size(); n++)
34  for (size_t m = 0; m < y0.size(); m++)
35  y[n][m] += y0[m];
36  }
37 
59  template <typename F>
60  struct coupled_ode_system <F, double, stan::math::var> {
61  const F& f_;
62  const std::vector<double>& y0_dbl_;
63  const std::vector<stan::math::var>& theta_;
64  std::vector<double> theta_dbl_;
65  const std::vector<double>& x_;
66  const std::vector<int>& x_int_;
67  const size_t N_;
68  const size_t M_;
69  const size_t size_;
70  std::ostream* msgs_;
71 
84  coupled_ode_system(const F& f,
85  const std::vector<double>& y0,
86  const std::vector<stan::math::var>& theta,
87  const std::vector<double>& x,
88  const std::vector<int>& x_int,
89  std::ostream* msgs)
90  : f_(f),
91  y0_dbl_(y0),
92  theta_(theta),
93  theta_dbl_(theta.size(), 0.0),
94  x_(x),
95  x_int_(x_int),
96  N_(y0.size()),
97  M_(theta.size()),
98  size_(N_ + N_ * M_),
99  msgs_(msgs) {
100  for (size_t m = 0; m < M_; m++)
101  theta_dbl_[m] = stan::math::value_of(theta[m]);
102  }
103 
118  void operator()(const std::vector<double>& y,
119  std::vector<double>& dy_dt,
120  double t) {
121  using std::vector;
122  using stan::math::var;
123 
124  vector<double> y_base(y.begin(), y.begin()+N_);
125  dy_dt = f_(t, y_base, theta_dbl_, x_, x_int_, msgs_);
126  stan::math::check_equal("coupled_ode_system",
127  "dy_dt", dy_dt.size(), N_);
128 
129  vector<double> coupled_sys(N_ * M_);
130  vector<var> theta_temp;
131  vector<var> y_temp;
132  vector<var> dy_dt_temp;
133  vector<double> grad;
134  vector<var> vars;
135 
136  for (size_t i = 0; i < N_; i++) {
137  theta_temp.clear();
138  y_temp.clear();
139  dy_dt_temp.clear();
140  grad.clear();
141  vars.clear();
142  try {
144  for (size_t j = 0; j < N_; j++) {
145  y_temp.push_back(y[j]);
146  vars.push_back(y_temp[j]);
147  }
148 
149  for (size_t j = 0; j < M_; j++) {
150  theta_temp.push_back(theta_dbl_[j]);
151  vars.push_back(theta_temp[j]);
152  }
153  dy_dt_temp = f_(t, y_temp, theta_temp, x_, x_int_, msgs_);
154  dy_dt_temp[i].grad(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[y_temp.size() + j];
161  for (size_t k = 0; k < N_; k++)
162  temp_deriv += y[N_ + N_ * j + k] * grad[k];
163 
164  coupled_sys[i + j * N_] = temp_deriv;
165  }
166  } catch (const std::exception& e) {
168  throw;
169  }
171  }
172 
173  dy_dt.insert(dy_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 
205 
212  std::vector<std::vector<stan::math::var> >
213  decouple_states(const std::vector<std::vector<double> >& y) {
215  std::vector<stan::math::var> temp_vars;
216  std::vector<double> temp_gradients;
217  std::vector<std::vector<stan::math::var> > y_return(y.size());
218 
219  for (size_t i = 0; i < y.size(); i++) {
220  temp_vars.clear();
221 
222  // iterate over number of equations
223  for (size_t j = 0; j < N_; j++) {
224  temp_gradients.clear();
225 
226  // iterate over parameters for each equation
227  for (size_t k = 0; k < M_; k++)
228  temp_gradients.push_back(y[i][y0_dbl_.size()
229  + y0_dbl_.size() * k + j]);
230 
231  temp_vars.push_back(precomputed_gradients(y[i][j],
232  theta_,
233  temp_gradients));
234  }
235  y_return[i] = temp_vars;
236  }
237 
238  return y_return;
239  }
240  };
241 
242 
243 
244 
245 
246 
273  template <typename F>
274  struct coupled_ode_system <F, stan::math::var, double> {
275  const F& f_;
276  const std::vector<stan::math::var>& y0_;
277  std::vector<double> y0_dbl_;
278  const std::vector<double>& theta_dbl_;
279  const std::vector<double>& x_;
280  const std::vector<int>& x_int_;
281  std::ostream* msgs_;
282  const size_t N_;
283  const size_t M_;
284  const size_t size_;
285 
299  coupled_ode_system(const F& f,
300  const std::vector<stan::math::var>& y0,
301  const std::vector<double>& theta,
302  const std::vector<double>& x,
303  const std::vector<int>& x_int,
304  std::ostream* msgs)
305  : f_(f),
306  y0_(y0),
307  y0_dbl_(y0.size(), 0.0),
308  theta_dbl_(theta),
309  x_(x),
310  x_int_(x_int),
311  msgs_(msgs),
312  N_(y0.size()),
313  M_(theta.size()),
314  size_(N_ + N_ * N_) {
315  for (size_t n = 0; n < N_; n++)
316  y0_dbl_[n] = stan::math::value_of(y0_[n]);
317  }
318 
332  void operator()(const std::vector<double>& y,
333  std::vector<double>& dy_dt,
334  double t) {
335  std::vector<double> y_base(y.begin(), y.begin() + N_);
336  for (size_t n = 0; n < N_; n++)
337  y_base[n] += y0_dbl_[n];
338 
339  dy_dt = f_(t, y_base, theta_dbl_, x_, x_int_, msgs_);
340  stan::math::check_equal("coupled_ode_system",
341  "dy_dt", dy_dt.size(), N_);
342 
343  std::vector<double> coupled_sys(N_ * N_);
344 
345  std::vector<stan::math::var> y_temp;
346  std::vector<stan::math::var> dy_dt_temp;
347  std::vector<double> grad;
348  std::vector<stan::math::var> vars;
349 
350  for (size_t i = 0; i < N_; i++) {
351  y_temp.clear();
352  dy_dt_temp.clear();
353  grad.clear();
354  vars.clear();
355  try {
357  for (size_t j = 0; j < N_; j++) {
358  y_temp.push_back(y[j] + y0_dbl_[j]);
359  vars.push_back(y_temp[j]);
360  }
361 
362  dy_dt_temp = f_(t, y_temp, theta_dbl_, x_, x_int_, msgs_);
363  dy_dt_temp[i].grad(vars, grad);
364 
365  for (size_t j = 0; j < N_; j++) {
366  // orders derivatives by equation (i.e. if there are 2 eqns
367  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
368  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
369  double temp_deriv = grad[j];
370  for (size_t k = 0; k < N_; k++)
371  temp_deriv += y[N_ + N_ * j + k] * grad[k];
372 
373  coupled_sys[i+j*N_] = temp_deriv;
374  }
375  } catch (const std::exception& e) {
377  throw;
378  }
380  }
381 
382  dy_dt.insert(dy_dt.end(), coupled_sys.begin(), coupled_sys.end());
383  }
384 
390  size_t size() const {
391  return size_;
392  }
393 
408  std::vector<double> initial_state() {
409  return std::vector<double>(size_, 0.0);
410  }
411 
419  std::vector<std::vector<stan::math::var> >
420  decouple_states(const std::vector<std::vector<double> >& y) {
422  using stan::math::var;
423  using std::vector;
424 
425  vector<var> temp_vars;
426  vector<double> temp_gradients;
427  vector<vector<var> > y_return(y.size());
428 
429  for (size_t i = 0; i < y.size(); i++) {
430  temp_vars.clear();
431 
432  // iterate over number of equations
433  for (size_t j = 0; j < N_; j++) {
434  temp_gradients.clear();
435 
436  // iterate over parameters for each equation
437  for (size_t k = 0; k < N_; k++)
438  temp_gradients.push_back(y[i][y0_.size() + y0_.size() * k + j]);
439 
440  temp_vars.push_back(precomputed_gradients(y[i][j],
441  y0_, temp_gradients));
442  }
443 
444  y_return[i] = temp_vars;
445  }
446 
447  add_initial_values(y0_, y_return);
448 
449  return y_return;
450  }
451  };
452 
453 
454 
455 
456 
457 
458 
494  template <typename F>
496  const F& f_;
497  const std::vector<stan::math::var>& y0_;
498  std::vector<double> y0_dbl_;
499  const std::vector<stan::math::var>& theta_;
500  std::vector<double> theta_dbl_;
501  const std::vector<double>& x_;
502  const std::vector<int>& x_int_;
503  const size_t N_;
504  const size_t M_;
505  const size_t size_;
506  std::ostream* msgs_;
507 
521  coupled_ode_system(const F& f,
522  const std::vector<stan::math::var>& y0,
523  const std::vector<stan::math::var>& theta,
524  const std::vector<double>& x,
525  const std::vector<int>& x_int,
526  std::ostream* msgs)
527  : f_(f),
528  y0_(y0),
529  y0_dbl_(y0.size(), 0.0),
530  theta_(theta),
531  theta_dbl_(theta.size(), 0.0),
532  x_(x),
533  x_int_(x_int),
534  N_(y0.size()),
535  M_(theta.size()),
536  size_(N_ + N_ * (N_ + M_)),
537  msgs_(msgs) {
538  for (size_t n = 0; n < N_; n++)
539  y0_dbl_[n] = stan::math::value_of(y0[n]);
540 
541  for (size_t m = 0; m < M_; m++)
542  theta_dbl_[m] = stan::math::value_of(theta[m]);
543  }
544 
558  void operator()(const std::vector<double>& y,
559  std::vector<double>& dy_dt,
560  double t) {
561  using std::vector;
562  using stan::math::var;
563 
564  vector<double> y_base(y.begin(), y.begin()+N_);
565  for (size_t n = 0; n < N_; n++)
566  y_base[n] += y0_dbl_[n];
567 
568  dy_dt = f_(t, y_base, theta_dbl_, x_, x_int_, msgs_);
569  stan::math::check_equal("coupled_ode_system",
570  "dy_dt", dy_dt.size(), N_);
571 
572  vector<double> coupled_sys(N_ * (N_ + M_));
573  vector<var> theta_temp;
574  vector<var> y_temp;
575  vector<var> dy_dt_temp;
576  vector<double> grad;
577  vector<var> vars;
578 
579  for (size_t i = 0; i < N_; i++) {
580  theta_temp.clear();
581  y_temp.clear();
582  dy_dt_temp.clear();
583  grad.clear();
584  vars.clear();
585  try {
587 
588  for (size_t j = 0; j < N_; j++) {
589  y_temp.push_back(y[j] + y0_dbl_[j]);
590  vars.push_back(y_temp[j]);
591  }
592 
593  for (size_t j = 0; j < M_; j++) {
594  theta_temp.push_back(theta_dbl_[j]);
595  vars.push_back(theta_temp[j]);
596  }
597 
598  dy_dt_temp = f_(t, y_temp, theta_temp, x_, x_int_, msgs_);
599  dy_dt_temp[i].grad(vars, grad);
600 
601  for (size_t j = 0; j < N_+M_; j++) {
602  // orders derivatives by equation (i.e. if there are 2 eqns
603  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
604  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
605  double temp_deriv = grad[j];
606  for (size_t k = 0; k < N_; k++)
607  temp_deriv += y[N_ + N_ * j + k] * grad[k];
608 
609  coupled_sys[i + j * N_] = temp_deriv;
610  }
611  } catch (const std::exception& e) {
613  throw;
614  }
616  }
617  dy_dt.insert(dy_dt.end(), coupled_sys.begin(), coupled_sys.end());
618  }
619 
625  size_t size() const {
626  return size_;
627  }
628 
640  std::vector<double> initial_state() {
641  return std::vector<double>(size_, 0.0);
642  }
643 
651  std::vector<std::vector<stan::math::var> >
652  decouple_states(const std::vector<std::vector<double> >& y) {
653  using std::vector;
654  using stan::math::var;
656 
657  vector<var> vars = y0_;
658  vars.insert(vars.end(), theta_.begin(), theta_.end());
659 
660  vector<var> temp_vars;
661  vector<double> temp_gradients;
662  vector<vector<var> > y_return(y.size());
663 
664  for (size_t i = 0; i < y.size(); i++) {
665  temp_vars.clear();
666 
667  // iterate over number of equations
668  for (size_t j = 0; j < N_; j++) {
669  temp_gradients.clear();
670 
671  // iterate over parameters for each equation
672  for (size_t k = 0; k < N_ + M_; k++)
673  temp_gradients.push_back(y[i][N_ + N_ * k + j]);
674 
675  temp_vars.push_back(precomputed_gradients(y[i][j],
676  vars, temp_gradients));
677  }
678  y_return[i] = temp_vars;
679  }
680  add_initial_values(y0_, y_return);
681  return y_return;
682  }
683  };
684 
685 
686  }
687 
688 }
689 
690 #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 grad(chainable *vi)
Compute the gradient for all variables starting from the specified root variable implementation.
Definition: grad.hpp:30
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:32
size_t size() const
Returns the size of the coupled system.
int M_
void operator()(const std::vector< double > &y, std::vector< double > &dy_dt, double t)
Assign the derivative vector with the system derivatives at the specified state and time...
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 ...
void operator()(const std::vector< double > &y, std::vector< double > &dy_dt, double t)
Populates the derivative vector with derivatives of the coupled ODE system state with respect to time...
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...
void operator()(const std::vector< double > &y, std::vector< double > &dy_dt, double t)
Calculates the derivative of the coupled ode system with respect to the state y at time t...
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)
Definition: size.hpp:11
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...
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.