1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_IDAS_FORWARD_SYSTEM_HPP 2 #define STAN_MATH_REV_MAT_FUNCTOR_IDAS_FORWARD_SYSTEM_HPP 11 #include <idas/idas.h> 12 #include <nvector/nvector_serial.h> 27 template <
typename F,
typename Tyy,
typename Typ,
typename Tpar>
48 const std::vector<Tyy>&
yy0,
const std::vector<Typ>&
yp0,
49 const std::vector<Tpar>&
theta,
50 const std::vector<double>& x_r,
51 const std::vector<int>& x_i, std::ostream* msgs)
52 :
idas_system<F, Tyy, Typ, Tpar>(f, eq_id, yy0, yp0, theta, x_r, x_i,
55 nv_yys_ = N_VCloneVectorArray(this->
ns_, this->
nv_yy_);
56 nv_yps_ = N_VCloneVectorArray(this->
ns_, this->
nv_yp_);
57 for (
size_t is = 0; is < this->
ns_; ++is) {
58 N_VConst(RCONST(0.0), nv_yys_[is]);
59 N_VConst(RCONST(0.0), nv_yps_[is]);
69 N_VDestroyVectorArray_Serial(this->nv_yys_, this->
ns_);
70 N_VDestroyVectorArray_Serial(this->nv_yps_, this->
ns_);
77 N_Vector*
nv_yys() {
return nv_yys_; }
82 N_Vector*
nv_yps() {
return nv_yps_; }
88 return static_cast<void*
>(
this);
95 return [](
int ns,
double t, N_Vector yy, N_Vector yp, N_Vector res,
96 N_Vector* yys, N_Vector* yps, N_Vector* ress,
void* user_data,
97 N_Vector temp1, N_Vector temp2, N_Vector temp3) {
99 using Eigen::MatrixXd;
100 using Eigen::VectorXd;
101 using Eigen::Dynamic;
104 DAE* dae =
static_cast<DAE*
>(user_data);
106 static const char* caller =
"sensitivity_residual";
109 const size_t& N = dae->N_;
110 const size_t& M = dae->M_;
112 Eigen::Map<VectorXd> vec_yy(N_VGetArrayPointer(yy), N);
113 Eigen::Map<VectorXd> vec_yp(N_VGetArrayPointer(yp), N);
114 std::vector<double> vyy(vec_yy.data(), vec_yy.data() + N);
115 std::vector<double> vyp(vec_yp.data(), vec_yp.data() + N);
116 std::vector<double> vtheta =
value_of(dae->theta());
118 std::vector<double> vpar =
value_of(dae->theta_);
119 Eigen::Map<VectorXd> vec_par(vpar.data(), vpar.size());
132 std::vector<var> yy(x.data(), x.data() + N);
134 = dae->f_(t, yy, vyp, vtheta, dae->x_r_, dae->x_i_, dae->msgs_);
135 Eigen::Map<vector_v> res(eval.data(), N);
143 std::vector<var> yp(x.data(), x.data() + N);
145 = dae->f_(t, vyy, yp, vtheta, dae->x_r_, dae->x_i_, dae->msgs_);
146 Eigen::Map<vector_v> res(eval.data(), N);
152 if (dae->is_var_par) {
155 std::vector<var> par(x.data(), x.data() + M);
157 = dae->f_(t, vyy, vyp, par, dae->x_r_, dae->x_i_, dae->msgs_);
158 Eigen::Map<vector_v> res(eval.data(), N);
162 r.block(0, (dae->is_var_yy0 ? N : 0) + (dae->is_var_yp0 ? N : 0), N,
168 }
catch (
const std::exception&
e) {
const std::vector< Typ > & yp0() const
return reference to derivative initial condition
N_Vector * nv_yps()
return N_Vector pointer array of sensitivity time derivative
Eigen::Matrix< var, Eigen::Dynamic, 1 > vector_v
The type of a (column) vector holding var values.
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.
IDAS DAE system that contains informtion on residual equation functor, sensitivity residual equation ...
IDAS DAE system with forward sensitivity calculation.
idas_forward_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.
Eigen::Matrix< var, Eigen::Dynamic, Eigen::Dynamic > matrix_v
The type of a matrix holding var values.
void jacobian(const F &f, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &fx, Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &J)
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 * to_user_data()
convert to void pointer for IDAS callbacks
const std::vector< Tyy > & yy0() const
return reference to initial condition
const size_t ns()
return number of sensitivity parameters
~idas_forward_system()
destructor to deallocate IDAS solution memory and workspace.
static constexpr bool need_sens
double e()
Return the base of the natural logarithm.
static void recover_memory_nested()
Recover only the memory used for the top nested call.
void check_greater(const char *function, const char *name, const T_y &y, const T_low &low)
Check if y is strictly greater than low.
const std::vector< Tpar > & theta() const
return reference to parameter
IDASensResFn sensitivity_residual() const
return a lambda for sensitivity residual callback.
static void start_nested()
Record the current position so that recover_memory_nested() can find it.
const F & f()
return reference to DAE functor
N_Vector * nv_yys()
return N_Vector pointer array of sensitivity