1 #ifndef STAN_MATH_REV_ARR_FUNCTOR_COUPLED_ODE_SYSTEM_HPP 2 #define STAN_MATH_REV_ARR_FUNCTOR_COUPLED_ODE_SYSTEM_HPP 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++)
67 const std::vector<double>&
x_;
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,
121 std::vector<double>& dz_dt,
125 vector<double>
grad(N_ + M_);
131 z_vars.reserve(N_ + M_);
133 vector<var> y_vars(z.begin(), z.begin() +
N_);
134 z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
136 vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
137 z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
139 vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
144 for (
size_t i = 0; i <
N_; i++) {
145 dz_dt[i] = dy_dt_vars[i].val();
147 dy_dt_vars[i].grad(z_vars, grad);
149 for (
size_t j = 0; j <
M_; j++) {
153 double temp_deriv = grad[N_ + j];
154 for (
size_t k = 0; k <
N_; k++)
155 temp_deriv += z[N_ + N_ * j + k] * grad[k];
157 dz_dt[N_ + i + j *
N_] = temp_deriv;
160 }
catch (
const std::exception&
e) {
190 std::vector<double> state(size_, 0.0);
191 for (
size_t n = 0; n <
N_; n++)
192 state[n] = y0_dbl_[n];
202 std::vector<std::vector<var> >
204 std::vector<var> temp_vars(N_);
205 std::vector<double> temp_gradients(M_);
206 std::vector<std::vector<var> > y_return(y.size());
208 for (
size_t i = 0; i < y.size(); i++) {
210 for (
size_t j = 0; j <
N_; j++) {
212 for (
size_t k = 0; k <
M_; k++)
213 temp_gradients[k] = y[i][y0_dbl_.size() + y0_dbl_.size() * k + j];
219 y_return[i] = temp_vars;
251 template <
typename F>
254 const std::vector<var>&
y0_;
257 const std::vector<double>&
x_;
278 const std::vector<var>& y0,
279 const std::vector<double>& theta,
280 const std::vector<double>& x,
281 const std::vector<int>& x_int,
292 size_(N_ + N_ * N_) {}
311 std::vector<double>& dz_dt,
315 std::vector<double> y(z.begin(), z.begin() +
N_);
316 for (
size_t n = 0; n <
N_; n++)
319 std::vector<double>
grad(N_);
327 vector<var> y_vars(y.begin(), y.end());
328 z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
330 vector<var> dy_dt_vars = f_(t, y_vars, theta_dbl_, x_, x_int_, msgs_);
335 for (
size_t i = 0; i <
N_; i++) {
336 dz_dt[i] = dy_dt_vars[i].val();
338 dy_dt_vars[i].grad(z_vars, grad);
340 for (
size_t j = 0; j <
N_; j++) {
344 double temp_deriv = grad[j];
345 for (
size_t k = 0; k <
N_; k++)
346 temp_deriv += z[N_ + N_ * j + k] * grad[k];
348 dz_dt[N_ + i + j *
N_] = temp_deriv;
351 }
catch (
const std::exception&
e) {
382 return std::vector<double>(
size_, 0.0);
392 std::vector<std::vector<var> >
396 vector<var> temp_vars(N_);
397 vector<double> temp_gradients(N_);
398 vector<vector<var> > y_return(y.size());
400 for (
size_t i = 0; i < y.size(); i++) {
402 for (
size_t j = 0; j <
N_; j++) {
404 for (
size_t k = 0; k <
N_; k++)
405 temp_gradients[k] = y[i][y0_.size() + y0_.size() * k + j];
408 y0_, temp_gradients);
410 y_return[i] = temp_vars;
454 template <
typename F>
457 const std::vector<var>&
y0_;
461 const std::vector<double>&
x_;
482 const std::vector<var>& y0,
483 const std::vector<var>& theta,
484 const std::vector<double>& x,
485 const std::vector<int>& x_int,
496 size_(N_ + N_ * (N_ + M_)),
516 std::vector<double>& dz_dt,
520 vector<double> y(z.begin(), z.begin() +
N_);
521 for (
size_t n = 0; n <
N_; n++)
524 vector<double>
grad(N_ + M_);
530 z_vars.reserve(N_ + M_);
532 vector<var> y_vars(y.begin(), y.end());
533 z_vars.insert(z_vars.end(), y_vars.begin(), y_vars.end());
535 vector<var> theta_vars(theta_dbl_.begin(), theta_dbl_.end());
536 z_vars.insert(z_vars.end(), theta_vars.begin(), theta_vars.end());
538 vector<var> dy_dt_vars = f_(t, y_vars, theta_vars, x_, x_int_, msgs_);
543 for (
size_t i = 0; i <
N_; i++) {
544 dz_dt[i] = dy_dt_vars[i].val();
546 dy_dt_vars[i].grad(z_vars, grad);
548 for (
size_t j = 0; j < N_ +
M_; j++) {
552 double temp_deriv = grad[j];
553 for (
size_t k = 0; k <
N_; k++)
554 temp_deriv += z[N_ + N_ * j + k] * grad[k];
556 dz_dt[N_ + i + j *
N_] = temp_deriv;
559 }
catch (
const std::exception&
e) {
587 return std::vector<double>(
size_, 0.0);
597 std::vector<std::vector<var> >
601 vector<var> vars = y0_;
602 vars.insert(vars.end(), theta_.begin(), theta_.end());
604 vector<var> temp_vars(N_);
605 vector<double> temp_gradients(N_ + M_);
606 vector<vector<var> > y_return(y.size());
608 for (
size_t i = 0; i < y.size(); i++) {
610 for (
size_t j = 0; j <
N_; j++) {
612 for (
size_t k = 0; k < N_ +
M_; k++)
613 temp_gradients[k] = y[i][N_ + N_ * k + j];
616 vars, temp_gradients);
618 y_return[i] = temp_vars;
const std::vector< double > y0_dbl_
const std::vector< double > theta_dbl_
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
Return the basic ODE solutions given the specified coupled system solutions, including the partials v...
const std::vector< var > & y0_
const std::vector< double > & x_
T value_of(const fvar< T > &v)
Return the value of the specified variable.
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
const std::vector< int > & x_int_
static void set_zero_all_adjoints_nested()
Reset all adjoint values in the top nested portion of the stack to zero.
const std::vector< var > & y0_
Independent (input) and dependent (output) variables for gradients.
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
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...
static void grad(vari *vi)
Compute the gradient for all variables starting from the specified root variable implementation.
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.
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
Return the solutions to the basic ODE system, including appropriate autodiff partial derivatives...
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
Populates the derivative vector with derivatives of the coupled ODE system state with respect to time...
var precomputed_gradients(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...
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
Assign the derivative vector with the system derivatives at the specified state and time...
void operator()(const std::vector< double > &z, std::vector< double > &dz_dt, double t) const
Calculates the derivative of the coupled ode system with respect to the state y at time t...
const std::vector< int > & x_int_
const std::vector< double > & y0_dbl_
const std::vector< var > & theta_
const std::vector< double > & theta_dbl_
size_t size() const
Returns the size of the coupled system.
const std::vector< double > y0_dbl_
double e()
Return the base of the natural logarithm.
const std::vector< int > & x_int_
size_t size() const
Returns the size of the coupled system.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
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.
std::vector< std::vector< var > > decouple_states(const std::vector< std::vector< double > > &y) const
Returns the base ODE system state corresponding to the specified coupled system state.
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 ...
const std::vector< double > & x_
std::vector< double > initial_state() const
Returns the initial state of the coupled system.
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...
const std::vector< double > & x_
const std::vector< double > theta_dbl_
static void start_nested()
Record the current position so that recover_memory_nested() can find it.
const std::vector< var > & theta_
size_t size() const
Returns the size of the coupled system.