Automatic Differentiation
 
Loading...
Searching...
No Matches
dae_system.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_REV_FUNCTOR_DAE_RESIDUAL_HPP
2#define STAN_MATH_REV_FUNCTOR_DAE_RESIDUAL_HPP
3
14#include <idas/idas.h>
15#include <nvector/nvector_serial.h>
16#include <ostream>
17#include <vector>
18#include <cmath>
19
20namespace stan {
21namespace math {
22
35template <typename F, typename Tyy, typename Typ, typename... T_par>
37 using dae_type = dae_system<F, Tyy, Typ, T_par...>;
38
39 protected:
40 const F& f_;
41 std::tuple<decltype(deep_copy_vars(std::declval<const T_par&>()))...>
43 std::tuple<plain_type_t<decltype(value_of(std::declval<const T_par&>()))>...>
45 std::tuple<const T_par&...> args_tuple_;
46 std::ostream* msgs_;
47
48 public:
49 Tyy const& yy;
50 Typ const& yp;
51 Eigen::VectorXd dbl_yy;
52 Eigen::VectorXd dbl_yp;
53 const size_t N;
54 const size_t M;
55 const size_t ns;
57 std::vector<stan::math::var> all_vars;
58 Eigen::VectorXd dbl_rr;
59
60 static constexpr bool is_var_yy0 = stan::is_var<return_type_t<Tyy>>::value;
61 static constexpr bool is_var_yp0 = stan::is_var<return_type_t<Typ>>::value;
62 static constexpr bool is_var_par
63 = stan::is_var<return_type_t<T_par...>>::value;
64 static constexpr bool use_fwd_sens = is_var_yy0 || is_var_yp0 || is_var_par;
65
66 using scalar_t = stan::return_type_t<Tyy, Typ, T_par...>;
67 using return_t = std::vector<Eigen::Matrix<scalar_t, -1, 1>>;
68
81 dae_system(const F& f, const Tyy& yy0, const Typ& yp0, std::ostream* msgs,
82 const T_par&... args)
83 : f_(f),
85 dbl_args_tuple_(value_of(args)...),
86 args_tuple_(std::forward_as_tuple(args...)),
87 msgs_(msgs),
88 yy(yy0),
89 yp(yp0),
90 dbl_yy(stan::math::value_of(yy0)),
91 dbl_yp(stan::math::value_of(yp0)),
92 N(yy0.size()),
93 M(count_vars(args...)),
94 ns((is_var_yy0 ? N : 0) + (is_var_yp0 ? N : 0) + M),
95 varis(ChainableStack::instance_->memalloc_.alloc_array<vari*>(ns)),
96 all_vars(ns),
97 dbl_rr(N) {
98 save_varis(varis, yy0, yp0, args...);
99 for (int i = 0; i < ns; ++i) {
100 all_vars[i].vi_ = *(varis + i);
101 }
102 }
103
104 void eval_residual(double t) {
106 [&](auto&&... args) { return f_(t, dbl_yy, dbl_yp, msgs_, args...); },
108 }
109
113 static int idas_res(double t, N_Vector yy, N_Vector yp, N_Vector rr,
114 void* user_data) {
115 dae_type* dae = static_cast<dae_type*>(user_data);
116 dae->dbl_yy = Eigen::Map<Eigen::VectorXd>(NV_DATA_S(yy), dae->N);
117 dae->dbl_yp = Eigen::Map<Eigen::VectorXd>(NV_DATA_S(yp), dae->N);
118 dae->eval_residual(t);
119 Eigen::Map<Eigen::VectorXd>(NV_DATA_S(rr), dae->N) = dae->dbl_rr;
120
121 return 0;
122 }
123
127 static int idas_sens_res(int ns, double t, N_Vector yy, N_Vector yp,
128 N_Vector res, N_Vector* yys, N_Vector* yps,
129 N_Vector* ress, void* user_data, N_Vector temp1,
130 N_Vector temp2, N_Vector temp3) {
131 dae_type* dae = static_cast<dae_type*>(user_data);
132
133 const int n = dae->N;
134 const int m = dae->M;
135
136 for (int i = 0; i < ns; ++i) {
137 N_VConst(0.0, ress[i]);
138 }
139
140 // Run nested autodiff in this scope
142
143 Eigen::Matrix<var, -1, 1> yy_var(n), yp_var(n);
144 for (int i = 0; i < n; ++i) {
145 yy_var.coeffRef(i) = NV_Ith_S(yy, i);
146 yp_var.coeffRef(i) = NV_Ith_S(yp, i);
147 }
148
149 Eigen::VectorXd g(m);
150 Eigen::Matrix<var, -1, 1> fy = math::apply(
151 [&](auto&&... args) {
152 return dae->f_(t, yy_var, yp_var, dae->msgs_, args...);
153 },
154 dae->local_args_tuple_);
155
156 for (int i = 0; i < n; ++i) {
157 if (i > 0) {
158 nested.set_zero_all_adjoints();
159 }
160 fy[i].grad();
161
162 for (int j = 0; j < ns; ++j) {
163 auto yysp = N_VGetArrayPointer(yys[j]);
164 auto ypsp = N_VGetArrayPointer(yps[j]);
165 auto ressp = N_VGetArrayPointer(ress[j]);
166 for (int k = 0; k < n; ++k) {
167 ressp[i] += yy_var[k].adj() * yysp[k] + yp_var[k].adj() * ypsp[k];
168 }
169 }
170
171 if (is_var_par) {
172 g.fill(0);
174 [&](auto&&... args) {
175 stan::math::accumulate_adjoints(g.data(), args...);
176 },
177 dae->local_args_tuple_);
178 for (int j = 0; j < m; ++j) {
179 auto ressp = N_VGetArrayPointer(ress[ns - m + j]);
180 ressp[i] += g[j];
181 }
182
184 dae->local_args_tuple_);
185 }
186 }
187 return 0;
188 }
189};
190
191} // namespace math
192} // namespace stan
193
194#endif
Eigen::VectorXd dbl_rr
std::tuple< const T_par &... > args_tuple_
stan::return_type_t< Tyy, Typ, T_par... > scalar_t
static constexpr bool is_var_par
static int idas_sens_res(int ns, double t, N_Vector yy, N_Vector yp, N_Vector res, N_Vector *yys, N_Vector *yps, N_Vector *ress, void *user_data, N_Vector temp1, N_Vector temp2, N_Vector temp3)
Evaluate DAE sensitivity residual according to IDAS signature.
static int idas_res(double t, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
Evaluate DAE residual according to IDAS signature.
void eval_residual(double t)
std::vector< stan::math::var > all_vars
Tyy const & yy
state variable y
Eigen::VectorXd dbl_yy
const F & f_
residual functor
std::tuple< plain_type_t< decltype(value_of(std::declval< const T_par & >()))> dbl_args_tuple_
dae_system(const F &f, const Tyy &yy0, const Typ &yp0, std::ostream *msgs, const T_par &... args)
Construct IDAS DAE system from initial condition and parameters.
Eigen::VectorXd dbl_yp
Typ const & yp
time derivatives of y
std::vector< Eigen::Matrix< scalar_t, -1, 1 > > return_t
static constexpr bool use_fwd_sens
std::tuple< decltype(deep_copy_vars(std::declval< const T_par & >()))... > local_args_tuple_
static constexpr bool is_var_yp0
static constexpr bool is_var_yy0
std::ostream * msgs_
IDAS DAE system that contains information on residual equation functor, sensitivity residual equation...
void set_zero_all_adjoints()
Reset all adjoint values in this nested stack to zero.
A class following the RAII idiom to start and recover nested autodiff scopes.
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
int64_t size(const T &m)
Returns the size (number of the elements) of a matrix_cl or var_value<matrix_cl<T>>.
Definition size.hpp:19
constexpr auto for_each(F &&f, T &&t)
Apply a function to each element of a tuple.
Definition for_each.hpp:66
fvar< T > arg(const std::complex< fvar< T > > &z)
Return the phase angle of the complex argument.
Definition arg.hpp:19
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
Arith deep_copy_vars(Arith &&arg)
Forward arguments that do not contain vars.
size_t count_vars(Pargs &&... args)
Count the number of vars in the input argument list.
vari ** save_varis(vari **dest, const var &x, Pargs &&... args)
Save the vari pointer in x into the memory pointed to by dest, increment the dest storage pointer,...
std::vector< Eigen::Matrix< stan::return_type_t< T_yy, T_yp, T_Args... >, -1, 1 > > dae(const F &f, const T_yy &yy0, const T_yp &yp0, double t0, const std::vector< double > &ts, std::ostream *msgs, const T_Args &... args)
Solve the DAE initial value problem f(t, y, y')=0, y(t0) = yy0, y'(t0)=yp0 at a set of times,...
Definition dae.hpp:167
void zero_adjoints() noexcept
End of recursion for set_zero_adjoints.
var_value< double > var
Definition var.hpp:1187
double * accumulate_adjoints(double *dest, const var &x, Pargs &&... args)
Accumulate adjoints from x into storage pointed to by dest, increment the adjoint storage pointer,...
constexpr decltype(auto) apply(F &&f, Tuple &&t, PreArgs &&... pre_args)
Definition apply.hpp:52
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
STL namespace.
Defines a static member named value which is defined to be false as the primitive scalar types cannot...
Definition is_var.hpp:14
This struct always provides access to the autodiff stack using the singleton pattern.