Stan Math Library  2.20.0
reverse mode automatic differentiation
idas_forward_system.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_IDAS_FORWARD_SYSTEM_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_IDAS_FORWARD_SYSTEM_HPP
3 
4 #include <stan/math/rev/meta.hpp>
11 #include <idas/idas.h>
12 #include <nvector/nvector_serial.h>
13 #include <ostream>
14 #include <vector>
15 
16 namespace stan {
17 namespace math {
18 
27 template <typename F, typename Tyy, typename Typ, typename Tpar>
28 class idas_forward_system : public idas_system<F, Tyy, Typ, Tpar> {
29  N_Vector* nv_yys_;
30  N_Vector* nv_yps_;
31 
32  public:
47  idas_forward_system(const F& f, const std::vector<int>& eq_id,
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,
53  msgs) {
54  if (this->need_sens) {
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]);
60  }
61  }
62  }
63 
68  if (this->need_sens) {
69  N_VDestroyVectorArray_Serial(this->nv_yys_, this->ns_);
70  N_VDestroyVectorArray_Serial(this->nv_yps_, this->ns_);
71  }
72  }
73 
77  N_Vector* nv_yys() { return nv_yys_; }
78 
82  N_Vector* nv_yps() { return nv_yps_; }
83 
87  void* to_user_data() { // prepare to inject DAE info
88  return static_cast<void*>(this);
89  }
90 
94  IDASensResFn sensitivity_residual() const {
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) {
98  using Eigen::Matrix;
99  using Eigen::MatrixXd;
100  using Eigen::VectorXd;
101  using Eigen::Dynamic;
102 
104  DAE* dae = static_cast<DAE*>(user_data);
105 
106  static const char* caller = "sensitivity_residual";
107  check_greater(caller, "number of parameters", ns, 0);
108 
109  const size_t& N = dae->N_;
110  const size_t& M = dae->M_;
111 
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());
117 
118  std::vector<double> vpar = value_of(dae->theta_);
119  Eigen::Map<VectorXd> vec_par(vpar.data(), vpar.size());
120 
121  auto yys_mat = matrix_d_from_NVarray(yys, ns);
122  auto yps_mat = matrix_d_from_NVarray(yps, ns);
123 
124  try {
126 
127  MatrixXd J, r;
128  VectorXd f_val;
129 
130  auto fyy
131  = [&t, &vyp, &vtheta, &N, &dae](const matrix_v& x) -> vector_v {
132  std::vector<var> yy(x.data(), x.data() + N);
133  auto eval
134  = dae->f_(t, yy, vyp, vtheta, dae->x_r_, dae->x_i_, dae->msgs_);
135  Eigen::Map<vector_v> res(eval.data(), N);
136  return res;
137  };
138  stan::math::jacobian(fyy, vec_yy, f_val, J);
139  r = J * yys_mat;
140 
141  auto fyp
142  = [&t, &vyy, &vtheta, &N, &dae](const matrix_v& x) -> vector_v {
143  std::vector<var> yp(x.data(), x.data() + N);
144  auto eval
145  = dae->f_(t, vyy, yp, vtheta, dae->x_r_, dae->x_i_, dae->msgs_);
146  Eigen::Map<vector_v> res(eval.data(), N);
147  return res;
148  };
149  stan::math::jacobian(fyp, vec_yp, f_val, J);
150  r += J * yps_mat;
151 
152  if (dae->is_var_par) {
153  auto fpar
154  = [&t, &vyy, &vyp, &N, &M, &dae](const matrix_v& x) -> vector_v {
155  std::vector<var> par(x.data(), x.data() + M);
156  auto eval
157  = dae->f_(t, vyy, vyp, par, dae->x_r_, dae->x_i_, dae->msgs_);
158  Eigen::Map<vector_v> res(eval.data(), N);
159  return res;
160  };
161  stan::math::jacobian(fpar, vec_par, f_val, J);
162  r.block(0, (dae->is_var_yy0 ? N : 0) + (dae->is_var_yp0 ? N : 0), N,
163  M)
164  += J; // only for theta
165  }
166 
167  matrix_d_to_NVarray(r, ress, ns);
168  } catch (const std::exception& e) {
170  throw;
171  }
172 
174 
175  return 0;
176  };
177  }
178 };
179 
180 } // namespace math
181 } // namespace stan
182 
183 #endif
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.
Definition: typedefs.hpp:23
Eigen::MatrixXd matrix_d_from_NVarray(const N_Vector *nv, const size_t &nv_size)
copy NV_Vector* array to Eigen::MatrixXd
Definition: idas_system.hpp:41
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
IDAS DAE system that contains informtion on residual equation functor, sensitivity residual equation ...
Definition: idas_system.hpp:90
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.
Definition: typedefs.hpp:17
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)
Definition: jacobian.hpp:11
void matrix_d_to_NVarray(const Eigen::MatrixXd &mat, N_Vector *nv, const size_t &nv_size)
copy Eigen::MatrixXd to NV_Vector* array.
Definition: idas_system.hpp:62
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.
Definition: constants.hpp:87
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

     [ Stan Home Page ] © 2011–2018, Stan Development Team.