1#ifndef STAN_MATH_PRIM_FUNCTOR_ODE_CKRK_HPP
2#define STAN_MATH_PRIM_FUNCTOR_ODE_CKRK_HPP
9#include <boost/numeric/odeint.hpp>
54template <
typename F,
typename T_y0,
typename T_t0,
typename T_ts,
55 typename... Args, require_eigen_vector_t<T_y0>* =
nullptr>
59 T_t0 t0,
const std::vector<T_ts>& ts,
60 double relative_tolerance,
double absolute_tolerance,
61 long int max_num_steps,
62 std::ostream* msgs,
const Args&... args) {
63 using boost::numeric::odeint::integrate_times;
64 using boost::numeric::odeint::make_dense_output;
65 using boost::numeric::odeint::max_step_checker;
66 using boost::numeric::odeint::no_progress_error;
67 using boost::numeric::odeint::runge_kutta_cash_karp54;
68 using boost::numeric::odeint::vector_space_algebra;
72 Eigen::Matrix<T_y0_t0, Eigen::Dynamic, 1> y0
73 = y0_arg.template cast<T_y0_t0>();
79 std::tuple<ref_type_t<Args>...> args_ref_tuple(args...);
82 [&](
const auto&... args_ref) {
84 std::vector<int> unused_temp{
86 (
check_finite(function_name,
"ode parameters and data", args_ref),
94 check_less(function_name,
"initial time", t0, ts[0]);
105 [&](
const auto&... args_ref) {
112 std::vector<double> ts_vec(ts.size() + 1);
114 for (
size_t i = 0; i < ts.size(); ++i)
117 std::vector<Eigen::Matrix<return_t, Eigen::Dynamic, 1>> y;
118 y.reserve(ts.size());
119 bool observer_initial_recorded =
false;
120 size_t time_index = 0;
124 auto filtered_observer
125 = [&](
const std::vector<double>& coupled_state,
double t) ->
void {
126 if (!observer_initial_recorded) {
127 observer_initial_recorded =
true;
131 [&](
const auto&... args_ref) {
133 f, coupled_state, y0, t0, ts[time_index], msgs, args_ref...));
140 std::vector<double> initial_coupled_state = coupled_system.initial_state();
142 const double step_size = 0.1;
145 make_controlled(absolute_tolerance, relative_tolerance,
146 runge_kutta_cash_karp54<std::vector<double>,
double,
147 std::vector<double>,
double>()),
148 std::ref(coupled_system), initial_coupled_state, std::begin(ts_vec),
149 std::end(ts_vec), step_size, filtered_observer,
150 max_step_checker(max_num_steps));
151 }
catch (
const no_progress_error&
e) {
153 "Failed to integrate to next output time (",
154 ") in less than max_num_steps steps");
195template <
typename F,
typename T_y0,
typename T_t0,
typename T_ts,
200 const std::vector<T_ts>& ts,
double relative_tolerance,
201 double absolute_tolerance,
202 long int max_num_steps,
203 std::ostream* msgs,
const Args&... args) {
205 relative_tolerance, absolute_tolerance,
206 max_num_steps, msgs, args...);
241template <
typename F,
typename T_y0,
typename T_t0,
typename T_ts,
245ode_ckrk(
const F& f,
const T_y0& y0, T_t0 t0,
const std::vector<T_ts>& ts,
246 std::ostream* msgs,
const Args&... args) {
247 double relative_tolerance = 1
e-6;
248 double absolute_tolerance = 1
e-6;
249 long int max_num_steps = 1e6;
252 absolute_tolerance, max_num_steps, msgs, args...);
require_t< is_eigen_vector< std::decay_t< T > > > require_eigen_vector_t
Require type satisfies is_eigen_vector.
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
std::vector< Eigen::Matrix< stan::return_type_t< T_y0, T_t0, T_ts, Args... >, Eigen::Dynamic, 1 > > ode_ckrk_tol(const F &f, const T_y0 &y0_arg, T_t0 t0, const std::vector< T_ts > &ts, double relative_tolerance, double absolute_tolerance, long int max_num_steps, std::ostream *msgs, const Args &... args)
Solve the ODE initial value problem y' = f(t, y), y(t0) = y0 at a set of times, { t1,...
std::vector< Eigen::Matrix< stan::return_type_t< T_y0, T_t0, T_ts, Args... >, Eigen::Dynamic, 1 > > ode_ckrk(const F &f, const T_y0 &y0, T_t0 t0, const std::vector< T_ts > &ts, std::ostream *msgs, const Args &... args)
Solve the ODE initial value problem y' = f(t, y), y(t0) = y0 at a set of times, { t1,...
static constexpr double e()
Return the base of the natural logarithm.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
void throw_domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
Eigen::VectorXd ode_store_sensitivities(const F &f, const std::vector< double > &coupled_state, const Eigen::Matrix< T_y0_t0, Eigen::Dynamic, 1 > &y0, T_t0 t0, T_t t, std::ostream *msgs, const Args &... args)
When all arguments are arithmetic, there are no sensitivities to store, so the function just returns ...
void check_nonzero_size(const char *function, const char *name, const T_y &y)
Check if the specified matrix/vector is of non-zero size.
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
void check_sorted(const char *function, const char *name, const EigVec &y)
Check if the specified vector is sorted into increasing order (repeated values are okay).
std::vector< Eigen::Matrix< stan::return_type_t< T_y0, T_t0, T_ts, Args... >, Eigen::Dynamic, 1 > > ode_ckrk_tol_impl(const char *function_name, const F &f, const T_y0 &y0_arg, T_t0 t0, const std::vector< T_ts > &ts, double relative_tolerance, double absolute_tolerance, long int max_num_steps, std::ostream *msgs, const Args &... args)
Solve the ODE initial value problem y' = f(t, y), y(t0) = y0 at a set of times, { t1,...
void check_less(const char *function, const char *name, const T_y &y, const T_high &high, Idxs... idxs)
Throw an exception if y is not strictly less than high.
constexpr decltype(auto) apply(F &&f, Tuple &&t, PreArgs &&... pre_args)
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...