1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_IDAS_RESIDUAL_HPP 2 #define STAN_MATH_REV_MAT_FUNCTOR_IDAS_RESIDUAL_HPP 13 #include <idas/idas.h> 14 #include <nvector/nvector_serial.h> 18 #define CHECK_IDAS_CALL(call) idas_check(call, #call) 28 std::ostringstream ss;
29 ss << func <<
" failed with error flag " << flag;
30 throw std::runtime_error(ss.str());
42 const size_t& nv_size) {
44 size_t n = NV_LENGTH_S(nv[0]);
46 for (
size_t j = 0; j < m; ++j) {
47 auto nvp = N_VGetArrayPointer(nv[j]);
48 for (
size_t i = 0; i < n; ++i) {
63 const size_t& nv_size) {
65 size_t n = NV_LENGTH_S(nv[0]);
66 for (
size_t j = 0; j < m; ++j) {
67 auto nvp = N_VGetArrayPointer(nv[j]);
68 for (
size_t i = 0; i < n; ++i) {
89 template <
typename F,
typename Tyy,
typename Typ,
typename Tpar>
93 const std::vector<Tyy>&
yy_;
94 const std::vector<Typ>&
yp_;
98 const std::vector<double>&
x_r_;
99 const std::vector<int>&
x_i_;
134 const std::vector<Tyy>&
yy0,
const std::vector<Typ>&
yp0,
135 const std::vector<Tpar>&
theta,
const std::vector<double>& x_r,
136 const std::vector<int>& x_i, std::ostream* msgs)
147 ns_((is_var_yy0 ? N_ : 0) + (is_var_yp0 ? N_ : 0)
148 + (is_var_par ? M_ : 0)),
149 nv_yy_(N_VMake_Serial(N_, yy_val_.data())),
150 nv_yp_(N_VMake_Serial(N_, yp_val_.data())),
152 nv_rr_(N_VMake_Serial(N_, rr_val_.data())),
153 id_(N_VNew_Serial(N_)),
156 if (nv_yy_ == NULL || nv_yp_ == NULL)
157 throw std::runtime_error(
"N_VMake_Serial failed to allocate memory");
160 throw std::runtime_error(
"IDACreate failed to allocate memory");
162 static const char* caller =
"idas_system";
170 "derivative initial state", yp0);
172 "derivative-algebra id", eq_id);
176 for (
size_t i = 0; i <
N_; ++i)
177 NV_Ith_S(id_, i) = eq_id[i];
184 N_VDestroy_Serial(nv_yy_);
185 N_VDestroy_Serial(nv_yp_);
186 N_VDestroy_Serial(nv_rr_);
187 N_VDestroy_Serial(id_);
238 const std::vector<Tyy>&
yy0()
const {
return yy_; }
245 const std::vector<Typ>&
yp0()
const {
return yp_; }
262 std::vector<scalar_type>
vars()
const {
263 std::vector<scalar_type> res;
265 res.insert(res.end(),
yy0().begin(),
yy0().end());
268 res.insert(res.end(),
yp0().begin(),
yp0().end());
271 res.insert(res.end(),
theta().begin(),
theta().end());
280 const size_t n() {
return N_; }
290 const size_t n_sys() {
return N_ * (ns_ + 1); }
295 const size_t n_par() {
return theta_.size(); }
305 const F&
f() {
return f_; }
311 return [](
double t, N_Vector yy, N_Vector yp, N_Vector rr,
312 void* user_data) ->
int {
314 DAE* dae =
static_cast<DAE*
>(user_data);
316 size_t N = NV_LENGTH_S(yy);
317 auto yy_val = N_VGetArrayPointer(yy);
319 auto yp_val = N_VGetArrayPointer(yp);
321 auto res = dae->f_(t, yy_vec, yp_vec, dae->theta_, dae->x_r_, dae->x_i_,
323 for (
size_t i = 0; i < N; ++i)
331 const std::vector<double> theta_d(
value_of(theta_));
332 const std::vector<double> yy_d(
value_of(yy_));
333 const std::vector<double> yp_d(
value_of(yp_));
334 static const char* caller =
"idas_integrator";
335 std::vector<double> res(
f_(t0, yy_d, yp_d, theta_d, x_r_, x_i_, msgs_));
static constexpr bool is_var_yy0
const std::vector< Typ > & yp0() const
return reference to derivative initial condition
const std::vector< int > & x_i_
idas_system(const F &f, const std::vector< int > &eq_id, const std::vector< Tyy > &yy0, const std::vector< Typ > &yp0, const std::vector< Tpar > &theta, const std::vector< double > &x_r, const std::vector< int > &x_i, std::ostream *msgs)
Construct IDAS DAE system from initial condition and parameters.
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
static constexpr bool is_var_yp0
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_less_or_equal(const char *function, const char *name, const T_y &y, const T_high &high)
Check if y is less or equal to high.
fvar< T > sqrt(const fvar< T > &x)
Eigen::MatrixXd matrix_d_from_NVarray(const N_Vector *nv, const size_t &nv_size)
copy NV_Vector* array to Eigen::MatrixXd
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Defines a public enum named value which is defined to be false as the primitive scalar types cannot b...
IDAS DAE system that contains informtion on residual equation functor, sensitivity residual equation ...
std::vector< scalar_type > vars() const
return a vector of vars for that contains the initial condition and parameters in case they are vars...
std::vector< double > yp_val_
fvar< T > dot_self(const Eigen::Matrix< fvar< T >, R, C > &v)
const size_t n()
return number of unknown variables
void matrix_d_to_NVarray(const Eigen::MatrixXd &mat, N_Vector *nv, const size_t &nv_size)
copy Eigen::MatrixXd to NV_Vector* array.
void check_greater_or_equal(const char *function, const char *name, const T_y &y, const T_low &low)
Check if y is greater or equal than low.
~idas_system()
destructor to deallocate IDAS solution memory and workspace.
std::vector< double > yy_val_
const std::vector< Tyy > & yy0() const
return reference to initial condition
void idas_check(int flag, const char *func)
check IDAS return flag & throw runtime error
static constexpr bool is_var_par
void * mem()
return IDAS memory handle
boost::math::tools::promote_args< double, typename scalar_type< T >::type, typename return_type< Types_pack... >::type >::type type
const size_t n_sys()
return size of DAE system for primary and sensitivity unknowns
N_Vector & nv_yp()
return reference to current N_Vector of derivative variable
const size_t n_par()
return theta size
const std::vector< double > & x_r_
const size_t ns()
return number of sensitivity parameters
const std::vector< Tyy > & yy_
IDAResFn residual()
return a closure for IDAS residual callback
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > matrix_d
Type for matrix of double values.
N_Vector & id()
return reference to DAE variable IDs
static constexpr bool need_sens
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
const std::vector< Tpar > & theta_
std::vector< std::vector< scalar_type > > return_type
void check_ic_consistency(const double &t0, const double &tol)
N_Vector & nv_yy()
return reference to current N_Vector of unknown variable
const std::vector< Tpar > & theta() const
return reference to parameter
const std::vector< Typ > & yp_
const std::vector< double > & yp_val()
return reference to current solution derivative vector value
const std::vector< double > & yy_val()
return reference to current solution vector value
void check_consistent_sizes(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)
Check if the dimension of x1 is consistent with x2.
std::vector< double > rr_val_
N_Vector & nv_rr()
return reference to current N_Vector of residual workspace
typename stan::return_type< Tyy, Typ, Tpar >::type scalar_type
const F & f()
return reference to DAE functor