1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_ODE_SYSTEM_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_ODE_SYSTEM_HPP
23 const std::vector<double> theta_;
24 const std::vector<double>& x_;
25 const std::vector<int>& x_int_;
40 const std::vector<double>& x,
const std::vector<int>& x_int,
42 : f_(f), theta_(theta), x_(x), x_int_(x_int), msgs_(msgs) { }
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_);
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 {
72 using Eigen::RowVectorXd;
75 vector<double>
grad(y.size());
76 Map<RowVectorXd> grad_eig(&
grad[0], y.size());
79 vector<var> y_var(y.begin(), y.end());
80 vector<var> dy_dt_var = f_(t, y_var, theta_, x_, x_int_, msgs_);
81 for (
size_t i = 0; i < dy_dt_var.size(); ++i) {
82 dy_dt(i) = dy_dt_var[i].val();
84 dy_dt_var[i].grad(y_var,
grad);
87 }
catch (
const std::exception&
e) {
106 template <
typename Derived1,
typename Derived2>
107 void jacobian(
const double t,
const std::vector<double>& y,
108 Eigen::MatrixBase<Derived1>& dy_dt,
109 Eigen::MatrixBase<Derived2>& Jy,
110 Eigen::MatrixBase<Derived2>& Jtheta)
const {
111 using Eigen::Dynamic;
114 using Eigen::RowVectorXd;
117 vector<double>
grad(y.size() + theta_.size());
118 Map<RowVectorXd> grad_eig(&
grad[0], y.size() + theta_.size());
121 vector<var> y_var(y.begin(), y.end());
122 vector<var> theta_var(theta_.begin(), theta_.end());
124 z_var.reserve(y.size() + theta_.size());
125 z_var.insert(z_var.end(), y_var.begin(), y_var.end());
126 z_var.insert(z_var.end(), theta_var.begin(), theta_var.end());
127 vector<var> dy_dt_var = f_(t, y_var, theta_var, x_, x_int_, msgs_);
128 for (
size_t i = 0; i < dy_dt_var.size(); ++i) {
129 dy_dt(i) = dy_dt_var[i].val();
131 dy_dt_var[i].grad(z_var,
grad);
132 Jy.row(i) = grad_eig.leftCols(y.size());
133 Jtheta.row(i) = grad_eig.rightCols(theta_.size());
135 }
catch (
const std::exception&
e) {
void operator()(const double t, const std::vector< double > &y, std::vector< double > &dy_dt) const
Calculate the RHS of the ODE.
static void set_zero_all_adjoints_nested()
Reset all adjoint values in the top nested portion of the stack to zero.
Independent (input) and dependent (output) variables for gradients.
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.
static void grad(vari *vi)
Compute the gradient for all variables starting from the specified root variable implementation.
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...
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.
Internal representation of an ODE model object which provides convenient Jacobian functions to obtain...
double e()
Return the base of the natural logarithm.
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.