Stan Math Library  2.20.0
reverse mode automatic differentiation
gp_periodic_cov.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_GP_PERIODIC_COV_HPP
2 #define STAN_MATH_REV_MAT_FUN_GP_PERIODIC_COV_HPP
3 
4 #include <stan/math/rev/meta.hpp>
11 #include <stan/math/rev/core.hpp>
13 #include <cmath>
14 #include <vector>
15 #include <type_traits>
16 
17 namespace stan {
18 namespace math {
19 
54 template <typename T_x, typename T_sigma, typename T_l, typename T_p>
55 class gp_periodic_cov_vari : public vari {
56  public:
57  const size_t size_;
58  const size_t size_ltri_;
59  const double l_d_;
60  const double sigma_d_;
61  const double p_d_;
62  const double sigma_sq_d_;
63  double *dist_;
64  double *sin_2_dist_;
65  double *sin_dist_sq_;
71 
91  gp_periodic_cov_vari(const std::vector<T_x> &x, const T_sigma &sigma,
92  const T_l &l, const T_p &p)
93  : vari(0.0),
94  size_(x.size()),
95  size_ltri_(size_ * (size_ - 1) / 2),
96  l_d_(value_of(l)),
97  sigma_d_(value_of(sigma)),
98  p_d_(value_of(p)),
99  sigma_sq_d_(sigma_d_ * sigma_d_),
100  dist_(ChainableStack::instance_->memalloc_.alloc_array<double>(
101  size_ltri_)),
102  sin_2_dist_(ChainableStack::instance_->memalloc_.alloc_array<double>(
103  size_ltri_)),
104  sin_dist_sq_(ChainableStack::instance_->memalloc_.alloc_array<double>(
105  size_ltri_)),
106  l_vari_(l.vi_),
107  sigma_vari_(sigma.vi_),
108  p_vari_(p.vi_),
109  cov_lower_(ChainableStack::instance_->memalloc_.alloc_array<vari *>(
110  size_ltri_)),
111  cov_diag_(
112  ChainableStack::instance_->memalloc_.alloc_array<vari *>(size_)) {
113  double neg_two_inv_l_sq = -2.0 / (l_d_ * l_d_);
114  double pi_div_p = pi() / p_d_;
115 
116  size_t pos = 0;
117  for (size_t j = 0; j < size_; ++j) {
118  for (size_t i = j + 1; i < size_; ++i) {
119  double dist = distance(x[i], x[j]);
120  double sin_dist = sin(pi_div_p * dist);
121  double sin_dist_sq = square(sin_dist);
122  dist_[pos] = dist;
123  sin_2_dist_[pos] = sin(2.0 * pi_div_p * dist);
124  sin_dist_sq_[pos] = sin_dist_sq;
125  cov_lower_[pos] = new vari(
126  sigma_sq_d_ * std::exp(sin_dist_sq * neg_two_inv_l_sq), false);
127  ++pos;
128  }
129  cov_diag_[j] = new vari(sigma_sq_d_, false);
130  }
131  }
132 
133  virtual void chain() {
134  double adjl = 0;
135  double adjsigma = 0;
136  double adjp = 0;
137 
138  for (size_t i = 0; i < size_ltri_; ++i) {
139  vari *el_low = cov_lower_[i];
140  double prod_add = el_low->adj_ * el_low->val_;
141  adjl += prod_add * sin_dist_sq_[i];
142  adjsigma += prod_add;
143  adjp += prod_add * sin_2_dist_[i] * dist_[i];
144  }
145  for (size_t i = 0; i < size_; ++i) {
146  vari *el = cov_diag_[i];
147  adjsigma += el->adj_ * el->val_;
148  }
149  double l_d_sq = l_d_ * l_d_;
150  l_vari_->adj_ += adjl * 4 / (l_d_sq * l_d_);
151  sigma_vari_->adj_ += adjsigma * 2 / sigma_d_;
152  p_vari_->adj_ += adjp * 2 * pi() / l_d_sq / (p_d_ * p_d_);
153  }
154 };
155 
188 template <typename T_x, typename T_l, typename T_p>
189 class gp_periodic_cov_vari<T_x, double, T_l, T_p> : public vari {
190  public:
191  const size_t size_;
192  const size_t size_ltri_;
193  const double l_d_;
194  const double sigma_d_;
195  const double p_d_;
196  const double sigma_sq_d_;
197  double *dist_;
198  double *sin_2_dist_;
199  double *sin_dist_sq_;
204 
224  gp_periodic_cov_vari(const std::vector<T_x> &x, double sigma, const T_l &l,
225  const T_p &p)
226  : vari(0.0),
227  size_(x.size()),
228  size_ltri_(size_ * (size_ - 1) / 2),
229  l_d_(value_of(l)),
230  sigma_d_(sigma),
231  p_d_(value_of(p)),
232  sigma_sq_d_(sigma_d_ * sigma_d_),
233  dist_(ChainableStack::instance_->memalloc_.alloc_array<double>(
234  size_ltri_)),
235  sin_2_dist_(ChainableStack::instance_->memalloc_.alloc_array<double>(
236  size_ltri_)),
237  sin_dist_sq_(ChainableStack::instance_->memalloc_.alloc_array<double>(
238  size_ltri_)),
239  l_vari_(l.vi_),
240  p_vari_(p.vi_),
241  cov_lower_(ChainableStack::instance_->memalloc_.alloc_array<vari *>(
242  size_ltri_)),
243  cov_diag_(
244  ChainableStack::instance_->memalloc_.alloc_array<vari *>(size_)) {
245  double neg_two_inv_l_sq = -2.0 / (l_d_ * l_d_);
246  double pi_div_p = pi() / p_d_;
247 
248  size_t pos = 0;
249  for (size_t j = 0; j < size_; ++j) {
250  for (size_t i = j + 1; i < size_; ++i) {
251  double dist = distance(x[i], x[j]);
252  double sin_dist = sin(pi_div_p * dist);
253  double sin_dist_sq = square(sin_dist);
254  dist_[pos] = dist;
255  sin_2_dist_[pos] = sin(2.0 * pi_div_p * dist);
256  sin_dist_sq_[pos] = sin_dist_sq;
257  cov_lower_[pos] = new vari(
258  sigma_sq_d_ * std::exp(sin_dist_sq * neg_two_inv_l_sq), false);
259  ++pos;
260  }
261  cov_diag_[j] = new vari(sigma_sq_d_, false);
262  }
263  }
264 
265  virtual void chain() {
266  double adjl = 0;
267  double adjp = 0;
268 
269  for (size_t i = 0; i < size_ltri_; ++i) {
270  vari *el_low = cov_lower_[i];
271  double prod_add = el_low->adj_ * el_low->val_;
272  adjl += prod_add * sin_dist_sq_[i];
273  adjp += prod_add * sin_2_dist_[i] * dist_[i];
274  }
275  double l_d_sq = l_d_ * l_d_;
276  l_vari_->adj_ += adjl * 4 / (l_d_sq * l_d_);
277  p_vari_->adj_ += adjp * 2 * pi() / l_d_sq / (p_d_ * p_d_);
278  }
279 };
280 
300 template <typename T_x>
301 inline typename std::enable_if<
302  std::is_same<typename scalar_type<T_x>::type, double>::value,
303  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>>::type
304 gp_periodic_cov(const std::vector<T_x> &x, const var &sigma, const var &l,
305  const var &p) {
306  const char *fun = "gp_periodic_cov";
307  check_positive(fun, "signal standard deviation", sigma);
308  check_positive(fun, "length-scale", l);
309  check_positive(fun, "period", p);
310  size_t x_size = x.size();
311  for (size_t i = 0; i < x_size; ++i)
312  check_not_nan(fun, "element of x", x[i]);
313 
314  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
315  if (x_size == 0)
316  return cov;
317 
319  = new gp_periodic_cov_vari<T_x, var, var, var>(x, sigma, l, p);
320 
321  size_t pos = 0;
322  for (size_t j = 0; j < x_size; ++j) {
323  for (size_t i = (j + 1); i < x_size; ++i) {
324  cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
325  cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
326  ++pos;
327  }
328  cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
329  }
330  return cov;
331 }
332 
352 template <typename T_x>
353 inline typename std::enable_if<
354  std::is_same<typename scalar_type<T_x>::type, double>::value,
355  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>>::type
356 gp_periodic_cov(const std::vector<T_x> &x, double sigma, const var &l,
357  const var &p) {
358  const char *fun = "gp_periodic_cov";
359  check_positive(fun, "signal standard deviation", sigma);
360  check_positive(fun, "length-scale", l);
361  check_positive(fun, "period", p);
362  size_t x_size = x.size();
363  for (size_t i = 0; i < x_size; ++i)
364  check_not_nan(fun, "element of x", x[i]);
365 
366  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
367  if (x_size == 0)
368  return cov;
369 
371  = new gp_periodic_cov_vari<T_x, double, var, var>(x, sigma, l, p);
372 
373  size_t pos = 0;
374  for (size_t j = 0; j < x_size - 1; ++j) {
375  for (size_t i = (j + 1); i < x_size; ++i) {
376  cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
377  cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
378  ++pos;
379  }
380  cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
381  }
382  cov.coeffRef(x_size - 1, x_size - 1).vi_ = baseVari->cov_diag_[x_size - 1];
383  return cov;
384 }
385 } // namespace math
386 } // namespace stan
387 #endif
This is a subclass of the vari class for precomputed gradients of gp_periodic_cov.
gp_periodic_cov_vari(const std::vector< T_x > &x, double sigma, const T_l &l, const T_p &p)
Constructor for gp_periodic_cov.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
The variable implementation base class.
Definition: vari.hpp:30
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:33
friend class var
Definition: vari.hpp:32
const double val_
The value of this variable.
Definition: vari.hpp:38
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
fvar< T > sin(const fvar< T > &x)
Definition: sin.hpp:11
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:11
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
gp_periodic_cov_vari(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &l, const T_p &p)
Constructor for gp_periodic_cov.
vari(double x)
Construct a variable implementation from a value.
Definition: vari.hpp:58
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
boost::math::tools::promote_args< T1, T2 >::type distance(const Eigen::Matrix< T1, R1, C1 > &v1, const Eigen::Matrix< T2, R2, C2 > &v2)
Returns the distance between the specified vectors.
Definition: distance.hpp:23
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
double pi()
Return the value of pi.
Definition: constants.hpp:80
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
This struct always provides access to the autodiff stack using the singleton pattern.
Eigen::Matrix< typename stan::return_type< T_x, T_sigma, T_l, T_p >::type, Eigen::Dynamic, Eigen::Dynamic > gp_periodic_cov(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &l, const T_p &p)
Returns a periodic covariance matrix using the input .

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