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 
22 /**
23  * The coupled ODE system for known initial values and unknown
24  * parameters.
25  *
26  * <p>If the base ODE state is size N and there are M parameters,
27  * the coupled system has N + N * M states.
28  * <p>The first N states correspond to the base system's N states:
29  * \f$ \frac{d x_n}{dt} \f$
30  *
31  * <p>The next M states correspond to the sensitivities of the
32  * parameters with respect to the first base system equation:
33  * \f[
34  * \frac{d x_{N+m}}{dt}
35  * = \frac{d}{dt} \frac{\partial x_1}{\partial \theta_m}
36  * \f]
37  *
38  * <p>The final M states correspond to the sensitivities with respect
39  * to the second base system equation, etc.
40  *
41  * @tparam F type of functor for the base ode system.
42  */
43 template <typename F>
44 struct coupled_ode_system<F, double, var> {
45  const F& f_;
46  const std::vector<double>& y0_dbl_;
47  const std::vector<var>& theta_;
48  const std::vector<double> theta_dbl_;
49  const std::vector<double>& x_;
50  const std::vector<int>& x_int_;
51  const size_t N_;
52  const size_t M_;
53  const size_t size_;
54  std::ostream* msgs_;
55 
56  /**
57  * Construct a coupled ODE system with the specified base
58  * ODE system, base initial state, parameters, data, and a
59  * message stream.
60  *
61  * @param[in] f the base ODE system functor.
62  * @param[in] y0 the initial state of the base ode.
63  * @param[in] theta parameters of the base ode.
64  * @param[in] x real data.
65  * @param[in] x_int integer data.
66  * @param[in, out] msgs stream to which messages are printed.
67  */
68  coupled_ode_system(const F& f, const std::vector<double>& y0,
69  const std::vector<var>& theta,
70  const std::vector<double>& x,
71  const std::vector<int>& x_int, std::ostream* msgs)
72  : f_(f),
73  y0_dbl_(y0),
74  theta_(theta),
75  theta_dbl_(value_of(theta)),
76  x_(x),
77  x_int_(x_int),
78  N_(y0.size()),
79  M_(theta.size()),
80  size_(N_ + N_ * M_),
81  msgs_(msgs) {}
82 
83  /**
84  * Assign the derivative vector with the system derivatives at
85  * the specified state and time.
86  *
87  * <p>The input state must be of size <code>size()</code>, and
88  * the output produced will be of the same size.
89  *
90  * @param[in] z state of the coupled ode system.
91  * @param[out] dz_dt populated with the derivatives of
92  * the coupled system at the specified state and time.
93  * @param[in] t time.
94  * @throw exception if the system function does not return the
95  * same number of derivatives as the state vector size.
96  *
97  * y is the base ODE system state
98  *
99  */
100  void operator()(const std::vector<double>& z, std::vector<double>& dz_dt,
101  double t) const {
102  using std::vector;
103 
104  vector<double> grad(N_ + M_);
105 
106  try {
107  start_nested();
108 
109  vector<var> z_vars;
110  z_vars.reserve(N_ + M_);
111 
112  vector<var> y_vars(z.begin(), z.begin() + N_);
113  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
114 
115  vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
116  z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
117 
118  vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
119 
120  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
121  "states", N_);
122 
123  for (size_t i = 0; i < N_; i++) {
124  dz_dt[i] = dy_dt_vars[i].val();
126  dy_dt_vars[i].grad(z_vars, grad);
127 
128  for (size_t j = 0; j < M_; j++) {
129  // orders derivatives by equation (i.e. if there are 2 eqns
130  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
131  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
132  double temp_deriv = grad[N_ + j];
133  for (size_t k = 0; k < N_; k++)
134  temp_deriv += z[N_ + N_ * j + k] * grad[k];
135 
136  dz_dt[N_ + i + j * N_] = temp_deriv;
137  }
138  }
139  } catch (const std::exception& e) {
141  throw;
142  }
144  }
145 
146  /**
147  * Returns the size of the coupled system.
148  *
149  * @return size of the coupled system.
150  */
151  size_t size() const { return size_; }
152 
153  /**
154  * Returns the initial state of the coupled system. Because the
155  * initial values are known, the initial state of the coupled
156  * system is the same as the initial state of the base ODE
157  * system.
158  *
159  * <p>This initial state returned is of size <code>size()</code>
160  * where the first N (base ODE system size) parameters are the
161  * initial conditions of the base ode system and the rest of the
162  * initial condition elements are 0.
163  *
164  * @return the initial condition of the coupled system.
165  */
166  std::vector<double> initial_state() const {
167  std::vector<double> state(size_, 0.0);
168  for (size_t n = 0; n < N_; n++)
169  state[n] = y0_dbl_[n];
170  return state;
171  }
172 
173  /**
174  * Returns the base ODE system state corresponding to the
175  * specified coupled system state.
176  *
177  * @param y coupled states after solving the ode
178  */
179  std::vector<std::vector<var> > decouple_states(
180  const std::vector<std::vector<double> >& y) const {
181  std::vector<var> temp_vars(N_);
182  std::vector<double> temp_gradients(M_);
183  std::vector<std::vector<var> > y_return(y.size());
184 
185  for (size_t i = 0; i < y.size(); i++) {
186  // iterate over number of equations
187  for (size_t j = 0; j < N_; j++) {
188  // iterate over parameters for each equation
189  for (size_t k = 0; k < M_; k++)
190  temp_gradients[k] = y[i][y0_dbl_.size() + y0_dbl_.size() * k + j];
191 
192  temp_vars[j] = precomputed_gradients(y[i][j], theta_, temp_gradients);
193  }
194  y_return[i] = temp_vars;
195  }
196  return y_return;
197  }
198 };
199 
200 /**
201  * The coupled ODE system for unknown initial values and known
202  * parameters.
203  *
204  * <p>If the original ODE has states of size N, the
205  * coupled system has N + N * N states. (derivatives of each
206  * state with respect to each initial value)
207  *
208  * <p>The coupled system has N + N * N states, where N is the size of
209  * the state vector in the base system.
210  *
211  * <p>The first N states correspond to the base system's N states:
212  * \f$ \frac{d x_n}{dt} \f$
213  *
214  * <p>The next N states correspond to the sensitivities of the initial
215  * conditions with respect to the to the first base system equation:
216  * \f[
217  * \frac{d x_{N+n}}{dt}
218  * = \frac{d}{dt} \frac{\partial x_1}{\partial y0_n}
219  * \f]
220  *
221  * <p>The next N states correspond to the sensitivities with respect
222  * to the second base system equation, etc.
223  *
224  * @tparam F type of base ODE system functor
225  */
226 template <typename F>
227 struct coupled_ode_system<F, var, double> {
228  const F& f_;
229  const std::vector<var>& y0_;
230  const std::vector<double> y0_dbl_;
231  const std::vector<double>& theta_dbl_;
232  const std::vector<double>& x_;
233  const std::vector<int>& x_int_;
234  std::ostream* msgs_;
235  const size_t N_;
236  const size_t M_;
237  const size_t size_;
238 
239  /**
240  * Construct a coupled ODE system for an unknown initial state
241  * and known parameters givne the specified base system functor,
242  * base initial state, parameters, data, and an output stream
243  * for messages.
244  *
245  * @param[in] f base ODE system functor.
246  * @param[in] y0 initial state of the base ODE.
247  * @param[in] theta system parameters.
248  * @param[in] x real data.
249  * @param[in] x_int integer data.
250  * @param[in, out] msgs output stream for messages.
251  */
252  coupled_ode_system(const F& f, const std::vector<var>& y0,
253  const std::vector<double>& theta,
254  const std::vector<double>& x,
255  const std::vector<int>& x_int, std::ostream* msgs)
256  : f_(f),
257  y0_(y0),
258  y0_dbl_(value_of(y0)),
259  theta_dbl_(theta),
260  x_(x),
261  x_int_(x_int),
262  msgs_(msgs),
263  N_(y0.size()),
264  M_(theta.size()),
265  size_(N_ + N_ * N_) {}
266 
267  /**
268  * Calculates the derivative of the coupled ode system
269  * with respect to the state y at time t.
270  *
271  * @param[in] z the current state of the coupled, shifted ode
272  * system. This is a a vector of double of length size().
273  * @param[out] dz_dt a vector of length size() with the
274  * derivatives of the coupled system evaluated with state y and
275  * time t.
276  * @param[in] t time.
277  * @throw exception if the system functor does not return a
278  * derivative vector of the same size as the state vector.
279  *
280  * y is the base ODE system state
281  *
282  */
283  void operator()(const std::vector<double>& z, std::vector<double>& dz_dt,
284  double t) const {
285  using std::vector;
286 
287  std::vector<double> grad(N_);
288 
289  try {
290  start_nested();
291 
292  vector<var> y_vars(z.begin(), z.begin() + N_);
293 
294  vector<var> dy_dt_vars = f_(t, y_vars, theta_dbl_, x_, x_int_, msgs_);
295 
296  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
297  "states", N_);
298 
299  for (size_t i = 0; i < N_; i++) {
300  dz_dt[i] = dy_dt_vars[i].val();
302  dy_dt_vars[i].grad(y_vars, grad);
303 
304  for (size_t j = 0; j < N_; j++) {
305  // orders derivatives by equation (i.e. if there are 2 eqns
306  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
307  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
308  double temp_deriv = 0;
309  for (size_t k = 0; k < N_; k++)
310  temp_deriv += z[N_ + N_ * j + k] * grad[k];
311 
312  dz_dt[N_ + i + j * N_] = temp_deriv;
313  }
314  }
315  } catch (const std::exception& e) {
317  throw;
318  }
320  }
321 
322  /**
323  * Returns the size of the coupled system.
324  *
325  * @return size of the coupled system.
326  */
327  size_t size() const { return size_; }
328 
329  /**
330  * Returns the initial state of the coupled system.
331  *
332  * <p>Because the starting state is unknown, the coupled system
333  * incorporates the initial conditions as parameters. The
334  * initial conditions for the coupled part of the system are set
335  * to zero along with the rest of the initial state, because the
336  * value of the initial state has been moved into the
337  * parameters.
338  *
339  * @return the initial condition of the coupled system.
340  * This is a vector of length size() where all elements
341  * are 0.
342  */
343  std::vector<double> initial_state() const {
344  std::vector<double> initial(size_, 0.0);
345  for (size_t i = 0; i < N_; i++)
346  initial[i] = y0_dbl_[i];
347  for (size_t i = 0; i < N_; i++)
348  initial[N_ + i * N_ + i] = 1.0;
349  return initial;
350  }
351 
352  /**
353  * Return the solutions to the basic ODE system, including
354  * appropriate autodiff partial derivatives, given the specified
355  * coupled system solution.
356  *
357  * @param y the vector of the coupled states after solving the ode
358  */
359  std::vector<std::vector<var> > decouple_states(
360  const std::vector<std::vector<double> >& y) const {
361  using std::vector;
362 
363  vector<var> temp_vars(N_);
364  vector<double> temp_gradients(N_);
365  vector<vector<var> > y_return(y.size());
366 
367  for (size_t i = 0; i < y.size(); i++) {
368  // iterate over number of equations
369  for (size_t j = 0; j < N_; j++) {
370  // iterate over parameters for each equation
371  for (size_t k = 0; k < N_; k++)
372  temp_gradients[k] = y[i][y0_.size() + y0_.size() * k + j];
373 
374  temp_vars[j] = precomputed_gradients(y[i][j], y0_, temp_gradients);
375  }
376  y_return[i] = temp_vars;
377  }
378 
379  return y_return;
380  }
381 };
382 
383 /**
384  * The coupled ode system for unknown intial values and unknown
385  * parameters.
386  *
387  * <p>The coupled system has N + N * (N + M) states, where N is
388  * size of the base ODE state vector and M is the number of
389  * parameters.
390  *
391  * <p>The first N states correspond to the base system's N states:
392  * \f$ \frac{d x_n}{dt} \f$
393  *
394  * <p>The next N+M states correspond to the sensitivities of the
395  * initial conditions, then to the parameters with respect to the
396  * to the first base system equation:
397  *
398  * \f[
399  * \frac{d x_{N + n}}{dt}
400  * = \frac{d}{dt} \frac{\partial x_1}{\partial y0_n}
401  * \f]
402  *
403  * \f[
404  * \frac{d x_{N+N+m}}{dt}
405  * = \frac{d}{dt} \frac{\partial x_1}{\partial \theta_m}
406  * \f]
407  *
408  * <p>The next N+M states correspond to the sensitivities with
409  * respect to the second base system equation, etc.
410  *
411  * <p>If the original ode has a state vector of size N states and
412  * a parameter vector of size M, the coupled system has N + N * (N
413  * + M) states. (derivatives of each state with respect to each
414  * initial value and each theta)
415  *
416  * @tparam F the functor for the base ode system
417  */
418 template <typename F>
420  const F& f_;
421  const std::vector<var>& y0_;
422  const std::vector<double> y0_dbl_;
423  const std::vector<var>& theta_;
424  const std::vector<double> theta_dbl_;
425  const std::vector<double>& x_;
426  const std::vector<int>& x_int_;
427  const size_t N_;
428  const size_t M_;
429  const size_t size_;
430  std::ostream* msgs_;
431 
432  /**
433  * Construct a coupled ODE system with unknown initial value and
434  * known parameters, given the base ODE system functor, the
435  * initial state of the base ODE, the parameters, data, and an
436  * output stream to which to write messages.
437  *
438  * @param[in] f the base ode system functor.
439  * @param[in] y0 the initial state of the base ode.
440  * @param[in] theta parameters of the base ode.
441  * @param[in] x real data.
442  * @param[in] x_int integer data.
443  * @param[in, out] msgs output stream to which to print messages.
444  */
445  coupled_ode_system(const F& f, const std::vector<var>& y0,
446  const std::vector<var>& theta,
447  const std::vector<double>& x,
448  const std::vector<int>& x_int, std::ostream* msgs)
449  : f_(f),
450  y0_(y0),
451  y0_dbl_(value_of(y0)),
452  theta_(theta),
453  theta_dbl_(value_of(theta)),
454  x_(x),
455  x_int_(x_int),
456  N_(y0.size()),
457  M_(theta.size()),
458  size_(N_ + N_ * (N_ + M_)),
459  msgs_(msgs) {}
460 
461  /**
462  * Populates the derivative vector with derivatives of the
463  * coupled ODE system state with respect to time evaluated at the
464  * specified state and specified time.
465  *
466  * @param[in] z the current state of the coupled, shifted ode system,
467  * of size <code>size()</code>.
468  * @param[in, out] dz_dt populate with the derivatives of the
469  * coupled system evaluated at the specified state and time.
470  * @param[in] t time.
471  * @throw exception if the base system does not return a
472  * derivative vector of the same size as the state vector.
473  *
474  * y is the base ODE system state
475  *
476  */
477  void operator()(const std::vector<double>& z, std::vector<double>& dz_dt,
478  double t) const {
479  using std::vector;
480 
481  vector<double> grad(N_ + M_);
482 
483  try {
484  start_nested();
485 
486  vector<var> z_vars;
487  z_vars.reserve(N_ + M_);
488 
489  vector<var> y_vars(z.begin(), z.begin() + N_);
490  z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
491 
492  vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
493  z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
494 
495  vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
496 
497  check_size_match("coupled_ode_system", "dz_dt", dy_dt_vars.size(),
498  "states", N_);
499 
500  for (size_t i = 0; i < N_; i++) {
501  dz_dt[i] = dy_dt_vars[i].val();
503  dy_dt_vars[i].grad(z_vars, grad);
504 
505  for (size_t j = 0; j < N_ + M_; j++) {
506  // orders derivatives by equation (i.e. if there are 2 eqns
507  // (y1, y2) and 2 parameters (a, b), dy_dt will be ordered as:
508  // dy1_dt, dy2_dt, dy1_da, dy2_da, dy1_db, dy2_db
509  double temp_deriv = j < N_ ? 0 : grad[j];
510  for (size_t k = 0; k < N_; k++)
511  temp_deriv += z[N_ + N_ * j + k] * grad[k];
512 
513  dz_dt[N_ + i + j * N_] = temp_deriv;
514  }
515  }
516  } catch (const std::exception& e) {
518  throw;
519  }
521  }
522 
523  /**
524  * Returns the size of the coupled system.
525  *
526  * @return size of the coupled system.
527  */
528  size_t size() const { return size_; }
529 
530  /**
531  * Returns the initial state of the coupled system.
532  *
533  * Because the initial state is unknown, the coupled system
534  * incorporates the initial condition offset from zero as
535  * a parameter, and hence the return of this function is a
536  * vector of zeros.
537  *
538  * @return the initial condition of the coupled system. This is
539  * a vector of length size() where all elements are 0.
540  */
541  std::vector<double> initial_state() const {
542  std::vector<double> initial(size_, 0.0);
543  for (size_t i = 0; i < N_; i++)
544  initial[i] = y0_dbl_[i];
545  for (size_t i = 0; i < N_; i++)
546  initial[N_ + i * N_ + i] = 1.0;
547  return initial;
548  }
549 
550  /**
551  * Return the basic ODE solutions given the specified coupled
552  * system solutions, including the partials versus the
553  * parameters encoded in the autodiff results.
554  *
555  * @param y the vector of the coupled states after solving the ode
556  */
557  std::vector<std::vector<var> > decouple_states(
558  const std::vector<std::vector<double> >& y) const {
559  using std::vector;
560 
561  vector<var> vars = y0_;
562  vars.insert(vars.end(), theta_.begin(), theta_.end());
563 
564  vector<var> temp_vars(N_);
565  vector<double> temp_gradients(N_ + M_);
566  vector<vector<var> > y_return(y.size());
567 
568  for (size_t i = 0; i < y.size(); i++) {
569  // iterate over number of equations
570  for (size_t j = 0; j < N_; j++) {
571  // iterate over parameters for each equation
572  for (size_t k = 0; k < N_ + M_; k++)
573  temp_gradients[k] = y[i][N_ + N_ * k + j];
574 
575  temp_vars[j] = precomputed_gradients(y[i][j], vars, temp_gradients);
576  }
577  y_return[i] = temp_vars;
578  }
579 
580  return y_return;
581  }
582 };
583 
584 } // namespace math
585 } // namespace stan
586 #endif
#define F(x, y, z)
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)
T value_of(const fvar< T > &v)
Definition: value_of.hpp:16
static void set_zero_all_adjoints_nested()
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
static void grad(vari *vi)
Definition: grad.hpp:30
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
var precomputed_gradients(double value, const std::vector< var > &operands, const std::vector< double > &gradients)
int M_
size_t size_
Definition: dot_self.hpp:18
const std::map< std::pair< std::string, std::string >, Variable > vars
const double j
Definition: BetheBloch.cxx:29
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
z
Definition: test.py:28
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
double e()
Definition: constants.hpp:89
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
int size(const std::vector< T > &x)
Definition: size.hpp:17
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)
static void recover_memory_nested()
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)
int N_
static void start_nested()