Stan Math Library  2.12.0
reverse mode automatic differentiation
ode_system.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_ODE_SYSTEM_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_ODE_SYSTEM_HPP
3 
4 #include <stan/math/rev/core.hpp>
6 #include <iostream>
7 #include <vector>
8 
9 namespace stan {
10  namespace math {
11 
20  template <typename F>
21  class ode_system {
22  const F& f_;
23  const std::vector<double> theta_;
24  const std::vector<double>& x_;
25  const std::vector<int>& x_int_;
26  std::ostream* msgs_;
27 
28  public:
39  ode_system(const F& f, const std::vector<double> theta,
40  const std::vector<double>& x, const std::vector<int>& x_int,
41  std::ostream* msgs)
42  : f_(f), theta_(theta), x_(x), x_int_(x_int), msgs_(msgs) { }
43 
51  inline void operator()(const double t, const std::vector<double>& y,
52  std::vector<double>& dy_dt) const {
53  dy_dt = f_(t, y, theta_, x_, x_int_, msgs_);
54  }
55 
66  template <typename Derived1, typename Derived2>
67  void jacobian(const double t, const std::vector<double>& y,
68  Eigen::MatrixBase<Derived1>& dy_dt,
69  Eigen::MatrixBase<Derived2>& Jy) const {
70  using Eigen::Matrix;
71  using Eigen::Map;
72  using Eigen::RowVectorXd;
73  using std::vector;
74  vector<double> grad(y.size());
75  Map<RowVectorXd> grad_eig(&grad[0], y.size());
76  try {
77  start_nested();
78  vector<var> y_var(y.begin(), y.end());
79  vector<var> dy_dt_var = f_(t, y_var, theta_, x_, x_int_, msgs_);
80  for (size_t i = 0; i < dy_dt_var.size(); ++i) {
81  dy_dt(i) = dy_dt_var[i].val();
83  dy_dt_var[i].grad(y_var, grad);
84  Jy.row(i) = grad_eig;
85  }
86  } catch (const std::exception& e) {
88  throw;
89  }
91  }
92 
105  template <typename Derived1, typename Derived2>
106  void jacobian(const double t, const std::vector<double>& y,
107  Eigen::MatrixBase<Derived1>& dy_dt,
108  Eigen::MatrixBase<Derived2>& Jy,
109  Eigen::MatrixBase<Derived2>& Jtheta) const {
110  using Eigen::Dynamic;
111  using Eigen::Map;
112  using Eigen::Matrix;
113  using Eigen::RowVectorXd;
114  using std::vector;
115  vector<double> grad(y.size() + theta_.size());
116  Map<RowVectorXd> grad_eig(&grad[0], y.size() + theta_.size());
117  try {
118  start_nested();
119  vector<var> y_var(y.begin(), y.end());
120  vector<var> theta_var(theta_.begin(), theta_.end());
121  vector<var> z_var;
122  z_var.reserve(y.size() + theta_.size());
123  z_var.insert(z_var.end(), y_var.begin(), y_var.end());
124  z_var.insert(z_var.end(), theta_var.begin(), theta_var.end());
125  vector<var> dy_dt_var = f_(t, y_var, theta_var, x_, x_int_, msgs_);
126  for (size_t i = 0; i < dy_dt_var.size(); ++i) {
127  dy_dt(i) = dy_dt_var[i].val();
129  dy_dt_var[i].grad(z_var, grad);
130  Jy.row(i) = grad_eig.leftCols(y.size());
131  Jtheta.row(i) = grad_eig.rightCols(theta_.size());
132  }
133  } catch (const std::exception& e) {
135  throw;
136  }
138  }
139  };
140 
141  }
142 }
143 #endif
void operator()(const double t, const std::vector< double > &y, std::vector< double > &dy_dt) const
Calculate the RHS of the ODE.
Definition: ode_system.hpp:51
static void set_zero_all_adjoints_nested()
Reset all adjoint values in the top nested portion of the stack to zero.
void jacobian(const double t, const std::vector< double > &y, Eigen::MatrixBase< Derived1 > &dy_dt, Eigen::MatrixBase< Derived2 > &Jy) const
Calculate the Jacobian of the ODE RHS wrt to states y.
Definition: ode_system.hpp:67
static void grad(vari *vi)
Compute the gradient for all variables starting from the specified root variable implementation.
Definition: grad.hpp:30
ode_system(const F &f, const std::vector< double > theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs)
Construct an ODE model with the specified base ODE system, parameters, data, and a message stream...
Definition: ode_system.hpp:39
void jacobian(const double t, const std::vector< double > &y, Eigen::MatrixBase< Derived1 > &dy_dt, Eigen::MatrixBase< Derived2 > &Jy, Eigen::MatrixBase< Derived2 > &Jtheta) const
Calculate the Jacobian of the ODE RHS wrt to states y and parameters theta.
Definition: ode_system.hpp:106
Internal representation of an ODE model object which provides convenient Jacobian functions to obtain...
Definition: ode_system.hpp:21
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:94
static void recover_memory_nested()
Recover only the memory used for the top nested call.
static void start_nested()
Record the current position so that recover_memory_nested() can find it.

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