Merge pull request #772 from drezap/feature/issue-756-rename-GP-functions

Feature/issue 756 rename cov_exp_quad -> gp_exp_quad_cov

stan/math/prim/mat/fun/cov_exp_quad.hpp

@@ -4,6 +4,7 @@
 #include <stan/math/prim/scal/meta/return_type.hpp>
 #include <stan/math/prim/scal/err/check_size_match.hpp>
 #include <stan/math/prim/mat/fun/Eigen.hpp>
+#include <stan/math/prim/mat/fun/gp_exp_quad_cov.hpp>
 #include <stan/math/prim/mat/fun/squared_distance.hpp>
 #include <stan/math/prim/scal/err/check_not_nan.hpp>
 #include <stan/math/prim/scal/err/check_positive.hpp>
@@ -16,19 +17,7 @@ namespace stan {
 namespace math {
 
 /**
- * Returns a squared exponential kernel.
- *
- * @tparam T_x type of std::vector of elements
- * @tparam T_sigma type of sigma
- * @tparam T_l type of length scale
- *
- * @param x std::vector of elements that can be used in square distance.
- *    This function assumes each element of x is the same size.
- * @param sigma standard deviation
- * @param length_scale length scale
- * @return squared distance
- * @throw std::domain_error if sigma <= 0, l <= 0, or
- *   x is nan or infinite
+ * @deprecated use <code>gp_exp_quad_cov</code>
  */
 template <typename T_x, typename T_sigma, typename T_l>
 inline
@@ -36,49 +25,11 @@ inline
                            Eigen::Dynamic, Eigen::Dynamic>
     cov_exp_quad(const std::vector<T_x>& x, const T_sigma& sigma,
                  const T_l& length_scale) {
-  using std::exp;
-  check_positive("cov_exp_quad", "marginal variance", sigma);
-  check_positive("cov_exp_quad", "length-scale", length_scale);
-  for (size_t n = 0; n < x.size(); ++n)
-    check_not_nan("cov_exp_quad", "x", x[n]);
-
-  Eigen::Matrix<typename stan::return_type<T_x, T_sigma, T_l>::type,
-                Eigen::Dynamic, Eigen::Dynamic>
-      cov(x.size(), x.size());
-
-  int x_size = x.size();
-  if (x_size == 0)
-    return cov;
-
-  T_sigma sigma_sq = square(sigma);
-  T_l neg_half_inv_l_sq = -0.5 / square(length_scale);
-
-  for (int j = 0; j < (x_size - 1); ++j) {
-    cov(j, j) = sigma_sq;
-    for (int i = j + 1; i < x_size; ++i) {
-      cov(i, j)
-          = sigma_sq * exp(squared_distance(x[i], x[j]) * neg_half_inv_l_sq);
-      cov(j, i) = cov(i, j);
-    }
-  }
-  cov(x_size - 1, x_size - 1) = sigma_sq;
-  return cov;
+  return gp_exp_quad_cov(x, sigma, length_scale);
 }
 
 /**
- * Returns a squared exponential kernel.
- *
- * @tparam T_x type of std::vector of elements
- * @tparam T_sigma type of sigma
- * @tparam T_l type of std::vector of length scale
- *
- * @param x std::vector of elements that can be used in square distance.
- *    This function assumes each element of x is the same size.
- * @param sigma standard deviation
- * @param length_scale std::vector length scale
- * @return squared distance
- * @throw std::domain_error if sigma <= 0, l <= 0, or
- *   x is nan or infinite
+ * @deprecated use <code>gp_exp_quad_cov</code>
  */
 template <typename T_x, typename T_sigma, typename T_l>
 inline
@@ -86,59 +37,11 @@ inline
                            Eigen::Dynamic, Eigen::Dynamic>
     cov_exp_quad(const std::vector<T_x>& x, const T_sigma& sigma,
                  const std::vector<T_l>& length_scale) {
-  using std::exp;
-
-  check_positive("cov_exp_quad", "marginal variance", sigma);
-  for (size_t n = 0; n < x.size(); ++n) {
-    check_not_nan("cov_exp_quad", "x", x[n]);
-  }
-  for (size_t n = 0; n < length_scale.size(); ++n) {
-    check_positive("cov_exp_quad", "length-scale", length_scale[n]);
-    check_not_nan("cov_exp_quad", "length-scale", length_scale[n]);
-  }
-
-  Eigen::Matrix<typename stan::return_type<T_x, T_sigma, T_l>::type,
-                Eigen::Dynamic, Eigen::Dynamic>
-      cov(x.size(), x.size());
-
-  int x_size = x.size();
-  int l_size = length_scale.size();
-  if (x_size == 0)
-    return cov;
-
-  T_sigma sigma_sq = square(sigma);
-  T_l temp_exp;
-
-  for (int j = 0; j < (x_size - 1); ++j) {
-    cov(j, j) = sigma_sq;
-    for (int i = j + 1; i < x_size; ++i) {
-      temp_exp = 0;
-      for (int k = 0; k < l_size; ++k) {
-        temp_exp += squared_distance(x[i], x[j]) / square(length_scale[k]);
-      }
-      cov(i, j) = sigma_sq * exp(-0.5 * temp_exp);
-      cov(j, i) = cov(i, j);
-    }
-  }
-  cov(x_size - 1, x_size - 1) = sigma_sq;
-  return cov;
+  return gp_exp_quad_cov(x, sigma, length_scale);
 }
 
 /**
- * Returns a squared exponential kernel.
- *
- * @tparam T_x1 type of first std::vector of elements
- * @tparam T_x2 type of second std::vector of elements
- * @tparam T_sigma type of sigma
- * @tparam T_l type of of length scale
- *
- * @param x1 std::vector of elements that can be used in square distance
- * @param x2 std::vector of elements that can be used in square distance
- * @param sigma standard deviation
- * @param length_scale length scale
- * @return squared distance
- * @throw std::domain_error if sigma <= 0, l <= 0, or
- *   x is nan or infinite
+ * @deprecated use <code>gp_exp_quad_cov</code>
  */
 template <typename T_x1, typename T_x2, typename T_sigma, typename T_l>
 inline typename Eigen::Matrix<
@@ -146,47 +49,11 @@ inline typename Eigen::Matrix<
     Eigen::Dynamic>
 cov_exp_quad(const std::vector<T_x1>& x1, const std::vector<T_x2>& x2,
              const T_sigma& sigma, const T_l& length_scale) {
-  using std::exp;
-  check_positive("cov_exp_quad", "marginal variance", sigma);
-  check_positive("cov_exp_quad", "length-scale", length_scale);
-  for (size_t n = 0; n < x1.size(); ++n)
-    check_not_nan("cov_exp_quad", "x1", x1[n]);
-  for (size_t n = 0; n < x2.size(); ++n)
-    check_not_nan("cov_exp_quad", "x2", x2[n]);
-
-  Eigen::Matrix<typename stan::return_type<T_x1, T_x2, T_sigma, T_l>::type,
-                Eigen::Dynamic, Eigen::Dynamic>
-      cov(x1.size(), x2.size());
-  if (x1.size() == 0 || x2.size() == 0)
-    return cov;
-
-  T_sigma sigma_sq = square(sigma);
-  T_l neg_half_inv_l_sq = -0.5 / square(length_scale);
-
-  for (size_t i = 0; i < x1.size(); ++i) {
-    for (size_t j = 0; j < x2.size(); ++j) {
-      cov(i, j)
-          = sigma_sq * exp(squared_distance(x1[i], x2[j]) * neg_half_inv_l_sq);
-    }
-  }
-  return cov;
+  return gp_exp_quad_cov(x1, x2, sigma, length_scale);
 }
 
 /**
- * Returns a squared exponential kernel.
- *
- * @tparam T_x1 type of first std::vector of elements
- * @tparam T_x2 type of second std::vector of elements
- * @tparam T_sigma type of sigma
- * @tparam T_l type of length scale
- *
- * @param x1 std::vector of elements that can be used in square distance
- * @param x2 std::vector of elements that can be used in square distance
- * @param sigma standard deviation
- * @param length_scale std::vector of length scale
- * @return squared distance
- * @throw std::domain_error if sigma <= 0, l <= 0, or
- *   x is nan or infinite
+ * @deprecated use <code>gp_exp_quad_cov</code>
  */
 template <typename T_x1, typename T_x2, typename T_sigma, typename T_l>
 inline typename Eigen::Matrix<
@@ -194,36 +61,7 @@ inline typename Eigen::Matrix<
     Eigen::Dynamic>
 cov_exp_quad(const std::vector<T_x1>& x1, const std::vector<T_x2>& x2,
              const T_sigma& sigma, const std::vector<T_l>& length_scale) {
-  using std::exp;
-  check_positive("cov_exp_quad", "marginal variance", sigma);
-  for (size_t n = 0; n < x1.size(); ++n)
-    check_not_nan("cov_exp_quad", "x1", x1[n]);
-  for (size_t n = 0; n < x2.size(); ++n)
-    check_not_nan("cov_exp_quad", "x2", x2[n]);
-  for (size_t n = 0; n < length_scale.size(); ++n)
-    check_positive("cov_exp_quad", "length-scale", length_scale[n]);
-  for (size_t n = 0; n < length_scale.size(); ++n)
-    check_not_nan("cov_exp_quad", "length-scale", length_scale[n]);
-
-  Eigen::Matrix<typename stan::return_type<T_x1, T_x2, T_sigma, T_l>::type,
-                Eigen::Dynamic, Eigen::Dynamic>
-      cov(x1.size(), x2.size());
-  if (x1.size() == 0 || x2.size() == 0)
-    return cov;
-
-  T_sigma sigma_sq = square(sigma);
-  T_l temp_exp;
-
-  for (size_t i = 0; i < x1.size(); ++i) {
-    for (size_t j = 0; j < x2.size(); ++j) {
-      temp_exp = 0;
-      for (size_t k = 0; k < length_scale.size(); ++k) {
-        temp_exp += squared_distance(x1[i], x2[j]) / square(length_scale[k]);
-      }
-      cov(i, j) = sigma_sq * exp(-0.5 * temp_exp);
-    }
-  }
-  return cov;
+  return gp_exp_quad_cov(x1, x2, sigma, length_scale);
 }
 }  // namespace math
 }  // namespace stan

stan/math/prim/mat/fun/dot_product.hpp

@@ -8,7 +8,6 @@
 
 namespace stan {
 namespace math {
-
 /**
  * Returns the dot product of the specified vectors.
  *
@@ -19,8 +18,8 @@ namespace math {
  * size or if they are both not vector dimensioned.
  */
 template <int R1, int C1, int R2, int C2>
-inline double dot_product(const Eigen::Matrix<double, R1, C1>& v1,
-                          const Eigen::Matrix<double, R2, C2>& v2) {
+inline double dot_product(const Eigen::Matrix<double, R1, C1> &v1,
+                          const Eigen::Matrix<double, R2, C2> &v2) {
   check_vector("dot_product", "v1", v1);
   check_vector("dot_product", "v2", v2);
   check_matching_sizes("dot_product", "v1", v1, "v2", v2);
@@ -32,7 +31,7 @@ inline double dot_product(const Eigen::Matrix<double, R1, C1>& v1,
  * @param v2 Second array.
  * @param length Length of both arrays.
  */
-inline double dot_product(const double* v1, const double* v2, size_t length) {
+inline double dot_product(const double *v1, const double *v2, size_t length) {
   double result = 0;
   for (size_t i = 0; i < length; i++)
     result += v1[i] * v2[i];
@@ -44,12 +43,11 @@ inline double dot_product(const double* v1, const double* v2, size_t length) {
  * @param v2 Second array.
  * @throw std::domain_error if the vectors are not the same size.
  */
-inline double dot_product(const std::vector<double>& v1,
-                          const std::vector<double>& v2) {
+inline double dot_product(const std::vector<double> &v1,
+                          const std::vector<double> &v2) {
   check_matching_sizes("dot_product", "v1", v1, "v2", v2);
   return dot_product(&v1[0], &v2[0], v1.size());
 }
-
 }  // namespace math
 }  // namespace stan
 #endif

stan/math/prim/mat/fun/gp_exp_quad_cov.hpp

+#ifndef STAN_MATH_PRIM_MAT_FUN_GP_EXP_QUAD_COV_HPP
+#define STAN_MATH_PRIM_MAT_FUN_GP_EXP_QUAD_COV_HPP
+
+#include <stan/math/prim/scal/meta/return_type.hpp>
+#include <stan/math/prim/scal/err/check_size_match.hpp>
+#include <stan/math/prim/mat/fun/Eigen.hpp>
+#include <stan/math/prim/mat/fun/squared_distance.hpp>
+#include <stan/math/prim/scal/err/check_not_nan.hpp>
+#include <stan/math/prim/scal/err/check_positive.hpp>
+#include <stan/math/prim/scal/fun/square.hpp>
+#include <stan/math/prim/scal/fun/exp.hpp>
+#include <vector>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns a squared exponential kernel.
+ *
+ * @tparam T_x type of std::vector of elements
+ * @tparam T_sigma type of sigma
+ * @tparam T_l type of length scale
+ *
+ * @param x std::vector of elements that can be used in square distance.
+ *    This function assumes each element of x is the same size.
+ * @param sigma standard deviation
+ * @param length_scale length scale
+ * @return squared distance
+ * @throw std::domain_error if sigma <= 0, l <= 0, or
+ *   x is nan or infinite
+ */
+template <typename T_x, typename T_sigma, typename T_l>
+inline
+    typename Eigen::Matrix<typename stan::return_type<T_x, T_sigma, T_l>::type,
+                           Eigen::Dynamic, Eigen::Dynamic>
+    gp_exp_quad_cov(const std::vector<T_x>& x, const T_sigma& sigma,
+                    const T_l& length_scale) {
+  using std::exp;
+  check_positive("gp_exp_quad_cov", "marginal variance", sigma);
+  check_positive("gp_exp_quad_cov", "length-scale", length_scale);
+  for (size_t n = 0; n < x.size(); ++n)
+    check_not_nan("gp_exp_quad_cov", "x", x[n]);
+
+  Eigen::Matrix<typename stan::return_type<T_x, T_sigma, T_l>::type,
+                Eigen::Dynamic, Eigen::Dynamic>
+      cov(x.size(), x.size());
+
+  int x_size = x.size();
+  if (x_size == 0)
+    return cov;
+
+  T_sigma sigma_sq = square(sigma);
+  T_l neg_half_inv_l_sq = -0.5 / square(length_scale);
+
+  for (int j = 0; j < (x_size - 1); ++j) {
+    cov(j, j) = sigma_sq;
+    for (int i = j + 1; i < x_size; ++i) {
+      cov(i, j)
+          = sigma_sq * exp(squared_distance(x[i], x[j]) * neg_half_inv_l_sq);
+      cov(j, i) = cov(i, j);
+    }
+  }
+  cov(x_size - 1, x_size - 1) = sigma_sq;
+  return cov;
+}
+
+/**
+ * Returns a squared exponential kernel.
+ *
+ * @tparam T_x type of std::vector of elements
+ * @tparam T_sigma type of sigma
+ * @tparam T_l type of std::vector of length scale
+ *
+ * @param x std::vector of elements that can be used in square distance.
+ *    This function assumes each element of x is the same size.
+ * @param sigma standard deviation
+ * @param length_scale std::vector length scale
+ * @return squared distance
+ * @throw std::domain_error if sigma <= 0, l <= 0, or
+ *   x is nan or infinite
+ */
+template <typename T_x, typename T_sigma, typename T_l>
+inline
+    typename Eigen::Matrix<typename stan::return_type<T_x, T_sigma, T_l>::type,
+                           Eigen::Dynamic, Eigen::Dynamic>
+    gp_exp_quad_cov(const std::vector<T_x>& x, const T_sigma& sigma,
+                    const std::vector<T_l>& length_scale) {
+  using std::exp;
+
+  check_positive("gp_exp_quad_cov", "marginal variance", sigma);
+  for (size_t n = 0; n < x.size(); ++n) {
+    check_not_nan("gp_exp_quad_cov", "x", x[n]);
+  }
+  for (size_t n = 0; n < length_scale.size(); ++n) {
+    check_positive("gp_exp_quad_cov", "length-scale", length_scale[n]);
+    check_not_nan("gp_exp_quad_cov", "length-scale", length_scale[n]);
+  }
+
+  Eigen::Matrix<typename stan::return_type<T_x, T_sigma, T_l>::type,
+                Eigen::Dynamic, Eigen::Dynamic>
+      cov(x.size(), x.size());
+
+  int x_size = x.size();
+  int l_size = length_scale.size();
+  if (x_size == 0)
+    return cov;
+
+  T_sigma sigma_sq = square(sigma);
+  T_l temp_exp;
+
+  for (int j = 0; j < (x_size - 1); ++j) {
+    cov(j, j) = sigma_sq;
+    for (int i = j + 1; i < x_size; ++i) {
+      temp_exp = 0;
+      for (int k = 0; k < l_size; ++k) {
+        temp_exp += squared_distance(x[i], x[j]) / square(length_scale[k]);
+      }
+      cov(i, j) = sigma_sq * exp(-0.5 * temp_exp);
+      cov(j, i) = cov(i, j);
+    }
+  }
+  cov(x_size - 1, x_size - 1) = sigma_sq;
+  return cov;
+}
+
+/**
+ * Returns a squared exponential kernel.
+ *
+ * @tparam T_x1 type of first std::vector of elements
+ * @tparam T_x2 type of second std::vector of elements
+ * @tparam T_sigma type of sigma
+ * @tparam T_l type of of length scale
+ *
+ * @param x1 std::vector of elements that can be used in square distance
+ * @param x2 std::vector of elements that can be used in square distance
+ * @param sigma standard deviation
+ * @param length_scale length scale
+ * @return squared distance
+ * @throw std::domain_error if sigma <= 0, l <= 0, or
+ *   x is nan or infinite
+ */
+template <typename T_x1, typename T_x2, typename T_sigma, typename T_l>
+inline typename Eigen::Matrix<
+    typename stan::return_type<T_x1, T_x2, T_sigma, T_l>::type, Eigen::Dynamic,
+    Eigen::Dynamic>
+gp_exp_quad_cov(const std::vector<T_x1>& x1, const std::vector<T_x2>& x2,
+                const T_sigma& sigma, const T_l& length_scale) {
+  using std::exp;
+  check_positive("gp_exp_quad_cov", "marginal variance", sigma);
+  check_positive("gp_exp_quad_cov", "length-scale", length_scale);
+  for (size_t n = 0; n < x1.size(); ++n)
+    check_not_nan("gp_exp_quad_cov", "x1", x1[n]);
+  for (size_t n = 0; n < x2.size(); ++n)
+    check_not_nan("gp_exp_quad_cov", "x2", x2[n]);
+
+  Eigen::Matrix<typename stan::return_type<T_x1, T_x2, T_sigma, T_l>::type,
+                Eigen::Dynamic, Eigen::Dynamic>
+      cov(x1.size(), x2.size());
+  if (x1.size() == 0 || x2.size() == 0)
+    return cov;
+
+  T_sigma sigma_sq = square(sigma);
+  T_l neg_half_inv_l_sq = -0.5 / square(length_scale);
+
+  for (size_t i = 0; i < x1.size(); ++i) {
+    for (size_t j = 0; j < x2.size(); ++j) {
+      cov(i, j)
+          = sigma_sq * exp(squared_distance(x1[i], x2[j]) * neg_half_inv_l_sq);
+    }
+  }
+  return cov;
+}
+
+/**
+ * Returns a squared exponential kernel.
+ *
+ * @tparam T_x1 type of first std::vector of elements
+ * @tparam T_x2 type of second std::vector of elements
+ * @tparam T_sigma type of sigma
+ * @tparam T_l type of length scale
+ *
+ * @param x1 std::vector of elements that can be used in square distance
+ * @param x2 std::vector of elements that can be used in square distance
+ * @param sigma standard deviation
+ * @param length_scale std::vector of length scale
+ * @return squared distance
+ * @throw std::domain_error if sigma <= 0, l <= 0, or
+ *   x is nan or infinite
+ */
+template <typename T_x1, typename T_x2, typename T_sigma, typename T_l>
+inline typename Eigen::Matrix<
+    typename stan::return_type<T_x1, T_x2, T_sigma, T_l>::type, Eigen::Dynamic,
+    Eigen::Dynamic>
+gp_exp_quad_cov(const std::vector<T_x1>& x1, const std::vector<T_x2>& x2,
+                const T_sigma& sigma, const std::vector<T_l>& length_scale) {
+  using std::exp;
+  check_positive("gp_exp_quad_cov", "marginal variance", sigma);
+  for (size_t n = 0; n < x1.size(); ++n)
+    check_not_nan("gp_exp_quad_cov", "x1", x1[n]);
+  for (size_t n = 0; n < x2.size(); ++n)
+    check_not_nan("gp_exp_quad_cov", "x2", x2[n]);
+  for (size_t n = 0; n < length_scale.size(); ++n)
+    check_positive("gp_exp_quad_cov", "length-scale", length_scale[n]);
+  for (size_t n = 0; n < length_scale.size(); ++n)
+    check_not_nan("gp_exp_quad_cov", "length-scale", length_scale[n]);
+
+  Eigen::Matrix<typename stan::return_type<T_x1, T_x2, T_sigma, T_l>::type,
+                Eigen::Dynamic, Eigen::Dynamic>
+      cov(x1.size(), x2.size());
+  if (x1.size() == 0 || x2.size() == 0)
+    return cov;
+
+  T_sigma sigma_sq = square(sigma);
+  T_l temp_exp;
+
+  for (size_t i = 0; i < x1.size(); ++i) {
+    for (size_t j = 0; j < x2.size(); ++j) {
+      temp_exp = 0;
+      for (size_t k = 0; k < length_scale.size(); ++k) {
+        temp_exp += squared_distance(x1[i], x2[j]) / square(length_scale[k]);
+      }
+      cov(i, j) = sigma_sq * exp(-0.5 * temp_exp);
+    }
+  }
+  return cov;
+}
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/rev/mat/fun/cov_exp_quad.hpp

@@ -20,17 +20,7 @@ namespace stan {
 namespace math {
 
 /**
- * This is a subclass of the vari class for precomputed
- * gradients of cov_exp_quad.
- *
- * The class stores the double values for the distance
- * matrix, pointers to the varis for the covariance
- * matrix, along with a pointer to the vari for sigma,
- * and the vari for l.
- *
- * @tparam T_x type of std::vector of elements
- * @tparam T_sigma type of sigma
- * @tparam T_l type of length scale
+ * @deprecated use <code>gp_exp_quad_cov_vari</code>
  */
 template <typename T_x, typename T_sigma, typename T_l>
 class cov_exp_quad_vari : public vari {
@@ -47,22 +37,7 @@ class cov_exp_quad_vari : public vari {
   vari** cov_diag_;
 
   /**
-   * Constructor for cov_exp_quad.
-   *
-   * All memory allocated in
-   * ChainableStack's stack_alloc arena.
-   *
-   * It is critical for the efficiency of this object
-   * that the constructor create new varis that aren't
-   * popped onto the var_stack_, but rather are
-   * popped onto the var_nochain_stack_. This is
-   * controlled to the second argument to
-   * vari's constructor.
-   *
-   * @param x std::vector input that can be used in square distance
-   *    Assumes each element of x is the same size
-   * @param sigma standard deviation
-   * @param l length scale
+   * @deprecated use <code>gp_exp_quad_cov_vari</code>
    */
   cov_exp_quad_vari(const std::vector<T_x>& x, const T_sigma& sigma,
                     const T_l& l)
@@ -115,16 +90,7 @@ class cov_exp_quad_vari : public vari {
 };
 
 /**
- * This is a subclass of the vari class for precomputed
- * gradients of cov_exp_quad.
- *
- * The class stores the double values for the distance
- * matrix, pointers to the varis for the covariance
- * matrix, along with a pointer to the vari for sigma,
- * and the vari for l.
- *
- * @tparam T_x type of std::vector of elements
- * @tparam T_l type of length scale
+ * @deprecated use <code>gp_exp_quad_cov_vari</code>
  */
 template <typename T_x, typename T_l>
 class cov_exp_quad_vari<T_x, double, T_l> : public vari {
@@ -140,22 +106,7 @@ class cov_exp_quad_vari<T_x, double, T_l> : public vari {
   vari** cov_diag_;
 
   /**
-   * Constructor for cov_exp_quad.
-   *
-   * All memory allocated in
-   * ChainableStack's stack_alloc arena.
-   *
-   * It is critical for the efficiency of this object
-   * that the constructor create new varis that aren't
-   * popped onto the var_stack_, but rather are
-   * popped onto the var_nochain_stack_. This is
-   * controlled to the second argument to
-   * vari's constructor.
-   *
-   * @param x std::vector input that can be used in square distance
-   *    Assumes each element of x is the same size
-   * @param sigma standard deviation
-   * @param l length scale
+   * @deprecated use <code>gp_exp_quad_cov_vari</code>
    */
   cov_exp_quad_vari(const std::vector<T_x>& x, double sigma, const T_l& l)
       : vari(0.0),
@@ -198,87 +149,25 @@ class cov_exp_quad_vari<T_x, double, T_l> : public vari {
 };
 
 /**
- * Returns a squared exponential kernel.
- *
- * @param x std::vector input that can be used in square distance
- *    Assumes each element of x is the same size
- * @param sigma standard deviation
- * @param l length scale
- * @return squared distance
- * @throw std::domain_error if sigma <= 0, l <= 0, or
- *   x is nan or infinite
+ * @deprecated use <code>gp_exp_quad_cov_vari</code>
  */
 template <typename T_x>
 inline typename boost::enable_if_c<
     boost::is_same<typename scalar_type<T_x>::type, double>::value,
     Eigen::Matrix<var, -1, -1> >::type
 cov_exp_quad(const std::vector<T_x>& x, const var& sigma, const var& l) {
-  check_positive("cov_exp_quad", "sigma", sigma);
-  check_positive("cov_exp_quad", "l", l);
-  size_t x_size = x.size();
-  for (size_t i = 0; i < x_size; ++i)
-    check_not_nan("cov_exp_quad", "x", x[i]);
-
-  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
-  if (x_size == 0)
-    return cov;
-
-  cov_exp_quad_vari<T_x, var, var>* baseVari
-      = new cov_exp_quad_vari<T_x, var, var>(x, sigma, l);
-
-  size_t pos = 0;
-  for (size_t j = 0; j < x_size - 1; ++j) {
-    for (size_t i = (j + 1); i < x_size; ++i) {
-      cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
-      cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
-      ++pos;
-    }
-    cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
-  }
-  cov.coeffRef(x_size - 1, x_size - 1).vi_ = baseVari->cov_diag_[x_size - 1];
-  return cov;
+  return gp_exp_quad_cov(x, sigma, l);
 }
 
 /**
- * Returns a squared exponential kernel.
- *
- * @param x std::vector input that can be used in square distance
- *    Assumes each element of x is the same size
- * @param sigma standard deviation
- * @param l length scale
- * @return squared distance
- * @throw std::domain_error if sigma <= 0, l <= 0, or
- *   x is nan or infinite
+ * @deprecated use <code>gp_exp_quad_cov_vari</code>
  */
 template <typename T_x>
 inline typename boost::enable_if_c<
     boost::is_same<typename scalar_type<T_x>::type, double>::value,
     Eigen::Matrix<var, -1, -1> >::type
 cov_exp_quad(const std::vector<T_x>& x, double sigma, const var& l) {
-  check_positive("cov_exp_quad", "marginal variance", sigma);
-  check_positive("cov_exp_quad", "length-scale", l);
-  size_t x_size = x.size();
-  for (size_t i = 0; i < x_size; ++i)
-    check_not_nan("cov_exp_quad", "x", x[i]);
-
-  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
-  if (x_size == 0)
-    return cov;
-
-  cov_exp_quad_vari<T_x, double, var>* baseVari
-      = new cov_exp_quad_vari<T_x, double, var>(x, sigma, l);
-
-  size_t pos = 0;
-  for (size_t j = 0; j < x_size - 1; ++j) {
-    for (size_t i = (j + 1); i < x_size; ++i) {
-      cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
-      cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
-      ++pos;
-    }
-    cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
-  }
-  cov.coeffRef(x_size - 1, x_size - 1).vi_ = baseVari->cov_diag_[x_size - 1];
-  return cov;
+  return gp_exp_quad_cov(x, sigma, l);
 }
 
 }  // namespace math

stan/math/rev/mat/fun/gp_exp_quad_cov.hpp

+#ifndef STAN_MATH_REV_MAT_FUN_GP_EXP_QUAD_COV_HPP
+#define STAN_MATH_REV_MAT_FUN_GP_EXP_QUAD_COV_HPP
+
+#include <boost/math/tools/promotion.hpp>
+#include <boost/type_traits.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <stan/math/prim/mat/fun/Eigen.hpp>
+#include <stan/math/prim/scal/err/check_not_nan.hpp>
+#include <stan/math/prim/scal/err/check_positive.hpp>
+#include <stan/math/prim/scal/fun/exp.hpp>
+#include <stan/math/prim/scal/fun/square.hpp>
+#include <stan/math/prim/scal/fun/squared_distance.hpp>
+#include <stan/math/prim/scal/meta/scalar_type.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/scal/fun/value_of.hpp>
+#include <vector>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * This is a subclass of the vari class for precomputed
+ * gradients of gp_exp_quad_cov.
+ *
+ * The class stores the double values for the distance
+ * matrix, pointers to the varis for the covariance
+ * matrix, along with a pointer to the vari for sigma,
+ * and the vari for length_scale.
+ *
+ * @tparam T_x type of std::vector of elements
+ * @tparam T_sigma type of sigma
+ * @tparam T_l type of length scale
+ */
+template <typename T_x, typename T_sigma, typename T_l>
+class gp_exp_quad_cov_vari : public vari {
+ public:
+  const size_t size_;
+  const size_t size_ltri_;
+  const double l_d_;
+  const double sigma_d_;
+  const double sigma_sq_d_;
+  double *dist_;
+  vari *l_vari_;
+  vari *sigma_vari_;
+  vari **cov_lower_;
+  vari **cov_diag_;
+
+  /**
+   * Constructor for gp_exp_quad_cov.
+   *
+   * All memory allocated in
+   * ChainableStack's stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled to the second argument to
+   * vari's constructor.
+   *
+   * @param x std::vector input that can be used in square distance
+   *    Assumes each element of x is the same size
+   * @param sigma standard deviation
+   * @param length_scale length scale
+   */
+  gp_exp_quad_cov_vari(const std::vector<T_x> &x, const T_sigma &sigma,
+                       const T_l &length_scale)
+      : vari(0.0),
+        size_(x.size()),
+        size_ltri_(size_ * (size_ - 1) / 2),
+        l_d_(value_of(length_scale)),
+        sigma_d_(value_of(sigma)),
+        sigma_sq_d_(sigma_d_ * sigma_d_),
+        dist_(ChainableStack::instance().memalloc_.alloc_array<double>(
+            size_ltri_)),
+        l_vari_(length_scale.vi_),
+        sigma_vari_(sigma.vi_),
+        cov_lower_(ChainableStack::instance().memalloc_.alloc_array<vari *>(
+            size_ltri_)),
+        cov_diag_(
+            ChainableStack::instance().memalloc_.alloc_array<vari *>(size_)) {
+    double inv_half_sq_l_d = 0.5 / (l_d_ * l_d_);
+    size_t pos = 0;
+    for (size_t j = 0; j < size_ - 1; ++j) {
+      for (size_t i = j + 1; i < size_; ++i) {
+        double dist_sq = squared_distance(x[i], x[j]);
+        dist_[pos] = dist_sq;
+        cov_lower_[pos] = new vari(
+            sigma_sq_d_ * std::exp(-dist_sq * inv_half_sq_l_d), false);
+        ++pos;
+      }
+    }
+    for (size_t i = 0; i < size_; ++i)
+      cov_diag_[i] = new vari(sigma_sq_d_, false);
+  }
+
+  virtual void chain() {
+    double adjl = 0;
+    double adjsigma = 0;
+
+    for (size_t i = 0; i < size_ltri_; ++i) {
+      vari *el_low = cov_lower_[i];
+      double prod_add = el_low->adj_ * el_low->val_;
+      adjl += prod_add * dist_[i];
+      adjsigma += prod_add;
+    }
+    for (size_t i = 0; i < size_; ++i) {
+      vari *el = cov_diag_[i];
+      adjsigma += el->adj_ * el->val_;
+    }
+    l_vari_->adj_ += adjl / (l_d_ * l_d_ * l_d_);
+    sigma_vari_->adj_ += adjsigma * 2 / sigma_d_;
+  }
+};
+
+/**
+ * This is a subclass of the vari class for precomputed
+ * gradients of gp_exp_quad_cov.
+ *
+ * The class stores the double values for the distance
+ * matrix, pointers to the varis for the covariance
+ * matrix, along with a pointer to the vari for sigma,
+ * and the vari for length_scale.
+ *
+ * @tparam T_x type of std::vector of elements
+ * @tparam T_l type of length scale
+ */
+template <typename T_x, typename T_l>
+class gp_exp_quad_cov_vari<T_x, double, T_l> : public vari {
+ public:
+  const size_t size_;
+  const size_t size_ltri_;
+  const double l_d_;
+  const double sigma_d_;
+  const double sigma_sq_d_;
+  double *dist_;
+  vari *l_vari_;
+  vari **cov_lower_;
+  vari **cov_diag_;
+
+  /**
+   * Constructor for gp_exp_quad_cov.
+   *
+   * All memory allocated in
+   * ChainableStack's stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled to the second argument to
+   * vari's constructor.
+   *
+   * @param x std::vector input that can be used in square distance
+   *    Assumes each element of x is the same size
+   * @param sigma standard deviation
+   * @param length_scale length scale
+   */
+  gp_exp_quad_cov_vari(const std::vector<T_x> &x, double sigma,
+                       const T_l &length_scale)
+      : vari(0.0),
+        size_(x.size()),
+        size_ltri_(size_ * (size_ - 1) / 2),
+        l_d_(value_of(length_scale)),
+        sigma_d_(value_of(sigma)),
+        sigma_sq_d_(sigma_d_ * sigma_d_),
+        dist_(ChainableStack::instance().memalloc_.alloc_array<double>(
+            size_ltri_)),
+        l_vari_(length_scale.vi_),
+        cov_lower_(ChainableStack::instance().memalloc_.alloc_array<vari *>(
+            size_ltri_)),
+        cov_diag_(
+            ChainableStack::instance().memalloc_.alloc_array<vari *>(size_)) {
+    double inv_half_sq_l_d = 0.5 / (l_d_ * l_d_);
+    size_t pos = 0;
+    for (size_t j = 0; j < size_ - 1; ++j) {
+      for (size_t i = j + 1; i < size_; ++i) {
+        double dist_sq = squared_distance(x[i], x[j]);
+        dist_[pos] = dist_sq;
+        cov_lower_[pos] = new vari(
+            sigma_sq_d_ * std::exp(-dist_sq * inv_half_sq_l_d), false);
+        ++pos;
+      }
+    }
+    for (size_t i = 0; i < size_; ++i)
+      cov_diag_[i] = new vari(sigma_sq_d_, false);
+  }
+
+  virtual void chain() {
+    double adjl = 0;
+
+    for (size_t i = 0; i < size_ltri_; ++i) {
+      vari *el_low = cov_lower_[i];
+      adjl += el_low->adj_ * el_low->val_ * dist_[i];
+    }
+    l_vari_->adj_ += adjl / (l_d_ * l_d_ * l_d_);
+  }
+};
+
+/**
+ * Returns a squared exponential kernel.
+ *
+ * @param x std::vector input that can be used in square distance
+ *    Assumes each element of x is the same size
+ * @param sigma standard deviation
+ * @param length_scale length scale
+ * @return squared distance
+ * @throw std::domain_error if sigma <= 0, l <= 0, or
+ *   x is nan or infinite
+ */
+template <typename T_x>
+inline typename boost::enable_if_c<
+    boost::is_same<typename scalar_type<T_x>::type, double>::value,
+    Eigen::Matrix<var, -1, -1>>::type
+gp_exp_quad_cov(const std::vector<T_x> &x, const var &sigma,
+                const var &length_scale) {
+  check_positive("gp_exp_quad_cov", "sigma", sigma);
+  check_positive("gp_exp_quad_cov", "length_scale", length_scale);
+  size_t x_size = x.size();
+  for (size_t i = 0; i < x_size; ++i)
+    check_not_nan("gp_exp_quad_cov", "x", x[i]);
+
+  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
+  if (x_size == 0)
+    return cov;
+
+  gp_exp_quad_cov_vari<T_x, var, var> *baseVari
+      = new gp_exp_quad_cov_vari<T_x, var, var>(x, sigma, length_scale);
+
+  size_t pos = 0;
+  for (size_t j = 0; j < x_size - 1; ++j) {
+    for (size_t i = (j + 1); i < x_size; ++i) {
+      cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
+      cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
+      ++pos;
+    }
+    cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
+  }
+  cov.coeffRef(x_size - 1, x_size - 1).vi_ = baseVari->cov_diag_[x_size - 1];
+  return cov;
+}
+
+/**
+ * Returns a squared exponential kernel.
+ *
+ * @param x std::vector input that can be used in square distance
+ *    Assumes each element of x is the same size
+ * @param sigma standard deviation
+ * @param length_scale length scale
+ * @return squared distance
+ * @throw std::domain_error if sigma <= 0, l <= 0, or
+ *   x is nan or infinite
+ */
+template <typename T_x>
+inline typename boost::enable_if_c<
+    boost::is_same<typename scalar_type<T_x>::type, double>::value,
+    Eigen::Matrix<var, -1, -1>>::type
+gp_exp_quad_cov(const std::vector<T_x> &x, double sigma,
+                const var &length_scale) {
+  check_positive("gp_exp_quad_cov", "marginal variance", sigma);
+  check_positive("gp_exp_quad_cov", "length-scale", length_scale);
+  size_t x_size = x.size();
+  for (size_t i = 0; i < x_size; ++i)
+    check_not_nan("gp_exp_quad_cov", "x", x[i]);
+
+  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
+  if (x_size == 0)
+    return cov;
+
+  gp_exp_quad_cov_vari<T_x, double, var> *baseVari
+      = new gp_exp_quad_cov_vari<T_x, double, var>(x, sigma, length_scale);
+
+  size_t pos = 0;
+  for (size_t j = 0; j < x_size - 1; ++j) {
+    for (size_t i = (j + 1); i < x_size; ++i) {
+      cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
+      cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
+      ++pos;
+    }
+    cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
+  }
+  cov.coeffRef(x_size - 1, x_size - 1).vi_ = baseVari->cov_diag_[x_size - 1];
+  return cov;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif