Automatic Differentiation
 
Loading...
Searching...
No Matches
ode_ckrk.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUNCTOR_ODE_CKRK_HPP
2#define STAN_MATH_PRIM_FUNCTOR_ODE_CKRK_HPP
3
9#include <boost/numeric/odeint.hpp>
10#include <ostream>
11#include <tuple>
12#include <vector>
13
14namespace stan {
15namespace math {
16
54template <typename F, typename T_y0, typename T_t0, typename T_ts,
55 typename... Args, require_eigen_vector_t<T_y0>* = nullptr>
56std::vector<Eigen::Matrix<stan::return_type_t<T_y0, T_t0, T_ts, Args...>,
57 Eigen::Dynamic, 1>>
58ode_ckrk_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
59 T_t0 t0, const std::vector<T_ts>& ts,
60 double relative_tolerance, double absolute_tolerance,
61 long int max_num_steps, // NOLINT(runtime/int)
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;
69
70 using T_y0_t0 = return_type_t<T_y0, T_t0>;
71
72 Eigen::Matrix<T_y0_t0, Eigen::Dynamic, 1> y0
73 = y0_arg.template cast<T_y0_t0>();
74
75 check_finite(function_name, "initial state", y0);
76 check_finite(function_name, "initial time", t0);
77 check_finite(function_name, "times", ts);
78
79 std::tuple<ref_type_t<Args>...> args_ref_tuple(args...);
80
82 [&](const auto&... args_ref) {
83 // Code from https://stackoverflow.com/a/17340003
84 std::vector<int> unused_temp{
85 0,
86 (check_finite(function_name, "ode parameters and data", args_ref),
87 0)...};
88 },
89 args_ref_tuple);
90
91 check_nonzero_size(function_name, "initial state", y0);
92 check_nonzero_size(function_name, "times", ts);
93 check_sorted(function_name, "times", ts);
94 check_less(function_name, "initial time", t0, ts[0]);
95
96 check_positive_finite(function_name, "relative_tolerance",
97 relative_tolerance);
98 check_positive_finite(function_name, "absolute_tolerance",
99 absolute_tolerance);
100 check_positive(function_name, "max_num_steps", max_num_steps);
101
102 using return_t = return_type_t<T_y0, T_t0, T_ts, Args...>;
103 // creates basic or coupled system by template specializations
104 auto&& coupled_system = math::apply(
105 [&](const auto&... args_ref) {
107 args_ref...);
108 },
109 args_ref_tuple);
110
111 // first time in the vector must be time of initial state
112 std::vector<double> ts_vec(ts.size() + 1);
113 ts_vec[0] = value_of(t0);
114 for (size_t i = 0; i < ts.size(); ++i)
115 ts_vec[i + 1] = value_of(ts[i]);
116
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;
121
122 // avoid recording of the initial state which is included by the
123 // conventions of odeint in the output
124 auto filtered_observer
125 = [&](const std::vector<double>& coupled_state, double t) -> void {
126 if (!observer_initial_recorded) {
127 observer_initial_recorded = true;
128 return;
129 }
131 [&](const auto&... args_ref) {
132 y.emplace_back(ode_store_sensitivities(
133 f, coupled_state, y0, t0, ts[time_index], msgs, args_ref...));
134 },
135 args_ref_tuple);
136 time_index++;
137 };
138
139 // the coupled system creates the coupled initial state
140 std::vector<double> initial_coupled_state = coupled_system.initial_state();
141
142 const double step_size = 0.1;
143 try {
144 integrate_times(
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) {
152 throw_domain_error(function_name, "", ts_vec[time_index + 1],
153 "Failed to integrate to next output time (",
154 ") in less than max_num_steps steps");
155 }
156
157 return y;
158}
159
195template <typename F, typename T_y0, typename T_t0, typename T_ts,
196 typename... Args, require_eigen_vector_t<T_y0>* = nullptr>
197std::vector<Eigen::Matrix<stan::return_type_t<T_y0, T_t0, T_ts, Args...>,
198 Eigen::Dynamic, 1>>
199ode_ckrk_tol(const F& f, const T_y0& y0_arg, T_t0 t0,
200 const std::vector<T_ts>& ts, double relative_tolerance,
201 double absolute_tolerance,
202 long int max_num_steps, // NOLINT(runtime/int)
203 std::ostream* msgs, const Args&... args) {
204 return ode_ckrk_tol_impl("ode_ckrk_tol", f, y0_arg, t0, ts,
205 relative_tolerance, absolute_tolerance,
206 max_num_steps, msgs, args...);
207}
208
241template <typename F, typename T_y0, typename T_t0, typename T_ts,
242 typename... Args, require_eigen_vector_t<T_y0>* = nullptr>
243std::vector<Eigen::Matrix<stan::return_type_t<T_y0, T_t0, T_ts, Args...>,
244 Eigen::Dynamic, 1>>
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 = 1e-6;
248 double absolute_tolerance = 1e-6;
249 long int max_num_steps = 1e6; // NOLINT(runtime/int)
250
251 return ode_ckrk_tol_impl("ode_ckrk", f, y0, t0, ts, relative_tolerance,
252 absolute_tolerance, max_num_steps, msgs, args...);
253}
254
255} // namespace math
256} // namespace stan
257#endif
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,...
Definition ode_ckrk.hpp:199
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,...
Definition ode_ckrk.hpp:245
static constexpr double e()
Return the base of the natural logarithm.
Definition constants.hpp:20
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
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,...
Definition ode_ckrk.hpp:58
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)
Definition apply.hpp:52
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 ...