Simplifying existing autodiff to prepare for adding var_value<Matrix<T, R, C>> support to rev/fun

stan/math/prim/fun/divide.hpp

@@ -40,9 +40,9 @@ inline int divide(int x, int y) {
  * @param[in] c specified scalar
  * @return matrix divided by the scalar
  */
-template <typename Mat, typename Scal, typename = require_eigen_t<Mat>,
-          typename = require_stan_scalar_t<Scal>,
-          typename = require_all_not_var_t<scalar_type_t<Mat>, Scal>>
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_all_not_st_var<Mat, Scal>* = nullptr>
 inline auto divide(const Mat& m, Scal c) {
   return (m / c).eval();
 }

stan/math/prim/fun/eigenvalues_sym.hpp

@@ -19,7 +19,8 @@ namespace math {
  * @param m Specified matrix.
  * @return Eigenvalues of matrix.
  */
-template <typename EigMat, require_eigen_t<EigMat>* = nullptr>
+template <typename EigMat, require_eigen_t<EigMat>* = nullptr,
+          require_not_vt_var<EigMat>* = nullptr>
 Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, 1> eigenvalues_sym(
     const EigMat& m) {
   using PlainMat = plain_type_t<EigMat>;

stan/math/prim/fun/eigenvectors_sym.hpp

@@ -8,7 +8,8 @@
 namespace stan {
 namespace math {
 
-template <typename EigMat, require_eigen_t<EigMat>* = nullptr>
+template <typename EigMat, require_eigen_t<EigMat>* = nullptr,
+          require_not_vt_var<EigMat>* = nullptr>
 Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic>
 eigenvectors_sym(const EigMat& m) {
   using PlainMat = plain_type_t<EigMat>;

stan/math/prim/fun/gp_exp_quad_cov.hpp

@@ -169,12 +169,13 @@ gp_exp_quad_cov(const std::vector<Eigen::Matrix<T_x1, -1, 1>> &x1,
  * @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>
+template <typename T_x, typename T_sigma, typename T_l,
+          require_all_not_st_var<T_x, T_sigma, T_l> * = nullptr>
 inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
                               Eigen::Dynamic>
 gp_exp_quad_cov(const std::vector<T_x> &x, const T_sigma &sigma,
                 const T_l &length_scale) {
-  check_positive("gp_exp_quad_cov", "magnitude", sigma);
+  check_positive("gp_exp_quad_cov", "marginal standard deviation", sigma);
   check_positive("gp_exp_quad_cov", "length scale", length_scale);
 
   const size_t x_size = x.size();
@@ -214,7 +215,8 @@ inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
                               Eigen::Dynamic>
 gp_exp_quad_cov(const std::vector<Eigen::Matrix<T_x, -1, 1>> &x,
                 const T_sigma &sigma, const std::vector<T_l> &length_scale) {
-  check_positive_finite("gp_exp_quad_cov", "magnitude", sigma);
+  check_positive_finite("gp_exp_quad_cov", "marginal standard deviation",
+                        sigma);
   check_positive_finite("gp_exp_quad_cov", "length scale", length_scale);
 
   size_t x_size = x.size();
@@ -258,7 +260,7 @@ inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l>,
 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) {
   const char *function_name = "gp_exp_quad_cov";
-  check_positive(function_name, "magnitude", sigma);
+  check_positive(function_name, "marginal standard deviation", sigma);
   check_positive(function_name, "length scale", length_scale);
 
   const size_t x1_size = x1.size();
@@ -326,7 +328,7 @@ gp_exp_quad_cov(const std::vector<Eigen::Matrix<T_x1, -1, 1>> &x1,
   for (size_t i = 0; i < x2_size; ++i) {
     check_not_nan(function_name, "x2", x2[i]);
   }
-  check_positive_finite(function_name, "magnitude", sigma);
+  check_positive_finite(function_name, "marginal standard deviation", sigma);
   check_positive_finite(function_name, "length scale", length_scale);
   check_size_match(function_name, "x dimension", x1[0].size(),
                    "number of length scales", l_size);
@@ -354,7 +356,7 @@ inline Eigen::MatrixXd gp_exp_quad_cov(const std::vector<double> &x,
                                        const double &sigma,
                                        const double &length_scale) {
   const char *function_name = "gp_exp_quad_cov";
-  check_positive(function_name, "magnitude", sigma);
+  check_positive(function_name, "marginal standard deviation", sigma);
   check_positive(function_name, "length scale", length_scale);
 
   const auto x_size = x.size();
@@ -397,7 +399,7 @@ inline Eigen::MatrixXd gp_exp_quad_cov(const std::vector<Eigen::VectorXd> &x,
                                        const double &sigma,
                                        const double &length_scale) {
   const char *function_name = "gp_exp_quad_cov";
-  check_positive(function_name, "magnitude", sigma);
+  check_positive(function_name, "marginal standard deviation", sigma);
   check_positive(function_name, "length scale", length_scale);
 
   const size_t x_size = x.size();
@@ -441,7 +443,7 @@ inline Eigen::MatrixXd gp_exp_quad_cov(
     const std::vector<Eigen::VectorXd> &x, const double &sigma,
     const std::vector<double> &length_scale) {
   const char *function_name = "gp_exp_quad_cov";
-  check_positive_finite(function_name, "magnitude", sigma);
+  check_positive_finite(function_name, "marginal standard deviation", sigma);
   check_positive_finite(function_name, "length scale", length_scale);
 
   const size_t x_size = x.size();
@@ -489,7 +491,7 @@ inline typename Eigen::MatrixXd gp_exp_quad_cov(const std::vector<double> &x1,
                                                 const double &sigma,
                                                 const double &length_scale) {
   const char *function_name = "gp_exp_quad_cov";
-  check_positive_finite(function_name, "magnitude", sigma);
+  check_positive_finite(function_name, "marginal standard deviation", sigma);
   check_positive_finite(function_name, "length scale", length_scale);
 
   Eigen::MatrixXd cov(x1.size(), x2.size());
@@ -542,7 +544,7 @@ inline typename Eigen::MatrixXd gp_exp_quad_cov(
   const char *function_name = "gp_exp_quad_cov";
   const int x1_size = x1.size();
   const int x2_size = x2.size();
-  check_positive_finite(function_name, "magnitude", sigma);
+  check_positive_finite(function_name, "marginal standard deviation", sigma);
   check_positive_finite(function_name, "length scale", length_scale);
 
   Eigen::MatrixXd cov(x1.size(), x2.size());
@@ -607,7 +609,7 @@ inline typename Eigen::MatrixXd gp_exp_quad_cov(
   const int x1_inner_size = x1[0].size();
   const int x2_inner_size = x1[0].size();
   const char *function_name = "gp_exp_quad_cov";
-  check_positive_finite(function_name, "magnitude", sigma);
+  check_positive_finite(function_name, "marginal standard deviation", sigma);
   check_positive_finite(function_name, "length scale", length_scale);
   check_size_match(function_name, "x dimension", x1[0].size(),
                    "number of length scales", l_size);

stan/math/prim/fun/gp_periodic_cov.hpp

@@ -41,7 +41,8 @@ namespace math {
  * @throw std::domain_error if sigma <= 0, l <= 0, p <= 0 or
  *   x is nan or infinite
  */
-template <typename T_x, typename T_sigma, typename T_l, typename T_p>
+template <typename T_x, typename T_sigma, typename T_l, typename T_p,
+          require_all_not_st_var<T_x, T_sigma, T_l, T_p> * = nullptr>
 inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l, T_p>,
                               Eigen::Dynamic, Eigen::Dynamic>
 gp_periodic_cov(const std::vector<T_x> &x, const T_sigma &sigma, const T_l &l,

stan/math/prim/fun/log_determinant_ldlt.hpp

@@ -17,7 +17,7 @@ namespace math {
  * @return the log(abs(det(A))
  */
 template <int R, int C, typename T>
-inline T log_determinant_ldlt(LDLT_factor<T, R, C> &A) {
+inline T log_determinant_ldlt(const LDLT_factor<T, R, C> &A) {
   if (A.rows() == 0) {
     return 0;
   }

stan/math/prim/fun/log_sum_exp.hpp

@@ -48,7 +48,7 @@ namespace math {
  * @param b the second variable
  */
 template <typename T1, typename T2, require_all_not_st_var<T1, T2>* = nullptr>
-inline return_type_t<T1, T2> log_sum_exp(const T2& a, const T1& b) {
+inline return_type_t<T1, T2> log_sum_exp(const T1& a, const T2& b) {
   if (a == NEGATIVE_INFTY) {
     return b;
   }

stan/math/prim/fun/matrix_power.hpp

@@ -29,11 +29,9 @@ matrix_power(const EigMat& M, const int n) {
   constexpr int C = EigMat::ColsAtCompileTime;
 
   check_square("matrix_power", "M", M);
-  if (n < 0)
-    invalid_argument("matrix_power", "n", n, "is ", ", but must be >= 0!");
+  check_nonnegative("matrix_power", "n", n);
   if (M.rows() == 0)
-    invalid_argument("matrix_power", "M.rows()", M.rows(), "is ",
-                     ", but must be > 0!");
+    return {};
   Eigen::Matrix<T, R, C> MM = M;
   check_finite("matrix_power", "M", MM);
   if (n == 0)

stan/math/prim/fun/ordered_constrain.hpp

@@ -19,7 +19,8 @@ namespace math {
  * @return Positive, increasing ordered vector.
  * @tparam T Type of scalar.
  */
-template <typename EigVec, require_eigen_col_vector_t<EigVec>* = nullptr>
+template <typename EigVec, require_eigen_col_vector_t<EigVec>* = nullptr,
+          require_not_st_var<EigVec>* = nullptr>
 plain_type_t<EigVec> ordered_constrain(const EigVec& x) {
   using Eigen::Dynamic;
   using Eigen::Matrix;

stan/math/prim/fun/positive_ordered_constrain.hpp

@@ -18,7 +18,8 @@ namespace math {
  * @param x Free vector of scalars.
  * @return Positive, increasing ordered vector.
  */
-template <typename EigVec, require_eigen_col_vector_t<EigVec>* = nullptr>
+template <typename EigVec, require_eigen_col_vector_t<EigVec>* = nullptr,
+          require_not_st_var<EigVec>* = nullptr>
 Eigen::Matrix<value_type_t<EigVec>, Eigen::Dynamic, 1>
 positive_ordered_constrain(const EigVec& x) {
   using Eigen::Dynamic;

stan/math/prim/fun/simplex_constrain.hpp

@@ -24,7 +24,8 @@ namespace math {
  * @param y Free vector input of dimensionality K - 1.
  * @return Simplex of dimensionality K.
  */
-template <typename ColVec, require_eigen_col_vector_t<ColVec>* = nullptr>
+template <typename ColVec, require_eigen_col_vector_t<ColVec>* = nullptr,
+          require_not_st_var<ColVec>* = nullptr>
 auto simplex_constrain(const ColVec& y) {
   // cut & paste simplex_constrain(Eigen::Matrix, T) w/o Jacobian
   using Eigen::Dynamic;

stan/math/prim/fun/trace_inv_quad_form_ldlt.hpp

@@ -25,7 +25,7 @@ namespace math {
  *
  */
 template <typename T, typename EigMat2, int R, int C,
-          typename = require_all_not_st_var<T, EigMat2>>
+          require_all_not_st_var<T, EigMat2>* = nullptr>
 inline return_type_t<T, EigMat2> trace_inv_quad_form_ldlt(
     const LDLT_factor<T, R, C>& A, const EigMat2& B) {
   check_multiplicable("trace_inv_quad_form_ldlt", "A", A, "B", B);

stan/math/prim/fun/unit_vector_constrain.hpp

@@ -23,12 +23,13 @@ namespace math {
  * @param y vector of K unrestricted variables
  * @return Unit length vector of dimension K
  */
-template <typename T, int R, int C>
-Eigen::Matrix<T, R, C> unit_vector_constrain(const Eigen::Matrix<T, R, C>& y) {
+template <typename T, require_eigen_t<T>* = nullptr,
+          require_not_vt_var<T>* = nullptr>
+inline plain_type_t<T> unit_vector_constrain(const T& y) {
   using std::sqrt;
   check_vector("unit_vector_constrain", "y", y);
   check_nonzero_size("unit_vector_constrain", "y", y);
-  T SN = dot_self(y);
+  scalar_type_t<T> SN = dot_self(y);
   check_positive_finite("unit_vector_constrain", "norm", SN);
   return y / sqrt(SN);
 }
@@ -45,13 +46,13 @@ Eigen::Matrix<T, R, C> unit_vector_constrain(const Eigen::Matrix<T, R, C>& y) {
  * @return Unit length vector of dimension K
  * @param lp Log probability reference to increment.
  */
-template <typename T, int R, int C>
-Eigen::Matrix<T, R, C> unit_vector_constrain(const Eigen::Matrix<T, R, C>& y,
-                                             T& lp) {
+template <typename T1, typename T2, require_eigen_t<T1>* = nullptr,
+          require_all_not_vt_var<T1, T2>* = nullptr>
+inline plain_type_t<T1> unit_vector_constrain(const T1& y, T2& lp) {
   using std::sqrt;
   check_vector("unit_vector_constrain", "y", y);
   check_nonzero_size("unit_vector_constrain", "y", y);
-  T SN = dot_self(y);
+  scalar_type_t<T1> SN = dot_self(y);
   check_positive_finite("unit_vector_constrain", "norm", SN);
   lp -= 0.5 * SN;
   return y / sqrt(SN);

stan/math/prim/fun/variance.hpp

@@ -4,6 +4,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/mean.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <vector>
 
 namespace stan {
@@ -45,20 +46,24 @@ inline return_type_t<T> variance(const std::vector<T>& v) {
  * @return sample variance of coefficients
  * @throw <code>std::invalid_argument</code> if the matrix has size zero
  */
-template <typename T, int R, int C>
-inline return_type_t<T> variance(const Eigen::Matrix<T, R, C>& m) {
+template <typename T, require_eigen_t<T>* = nullptr,
+          require_not_st_var<T>* = nullptr>
+inline return_type_t<T> variance(const T& m) {
   check_nonzero_size("variance", "m", m);
 
   if (m.size() == 1) {
     return 0.0;
   }
-  return_type_t<T> mn(mean(m));
+
+  const auto& m_ref = to_ref(m);
+
+  return_type_t<T> mn(mean(m_ref));
   return_type_t<T> sum_sq_diff(0);
-  for (int i = 0; i < m.size(); ++i) {
-    return_type_t<T> diff = m(i) - mn;
+  for (int i = 0; i < m_ref.size(); ++i) {
+    return_type_t<T> diff = m_ref(i) - mn;
     sum_sq_diff += diff * diff;
   }
-  return sum_sq_diff / (m.size() - 1);
+  return sum_sq_diff / (m_ref.size() - 1);
 }
 
 }  // namespace math

stan/math/prim/prob/inv_wishart_lpdf.hpp

@@ -53,8 +53,8 @@ return_type_t<T_y, T_dof, T_scale> inv_wishart_lpdf(const T_y& W,
   static const char* function = "inv_wishart_lpdf";
   check_size_match(function, "Rows of random variable", W.rows(),
                    "columns of scale parameter", S.rows());
-  check_square(function, "random variable", W);
-  check_square(function, "scale parameter", S);
+  check_symmetric(function, "random variable", W);
+  check_symmetric(function, "scale parameter", S);
   Eigen::Index k = S.rows();
   T_nu_ref nu_ref = nu;
   T_S_ref S_ref = S;
@@ -87,8 +87,7 @@ return_type_t<T_y, T_dof, T_scale> inv_wishart_lpdf(const T_y& W,
     //      &S(0), S.size(), 1);
     //    lp -= 0.5 * dot_product(S_vec, W_inv_vec); // trace(S * W^-1)
     Eigen::Matrix<return_type_t<T_y, T_scale>, Eigen::Dynamic, Eigen::Dynamic>
-        Winv_S(mdivide_left_ldlt(
-            ldlt_W, S_ref.template selfadjointView<Eigen::Lower>()));
+        Winv_S(mdivide_left_ldlt(ldlt_W, S_ref));
     lp -= 0.5 * trace(Winv_S);
   }
   if (include_summand<propto, T_dof, T_scale>::value) {

stan/math/prim/prob/inv_wishart_rng.hpp

@@ -16,7 +16,7 @@ inline Eigen::MatrixXd inv_wishart_rng(double nu, const Eigen::MatrixXd& S,
   static const char* function = "inv_wishart_rng";
   index_type_t<MatrixXd> k = S.rows();
   check_greater(function, "degrees of freedom > dims - 1", nu, k - 1);
-  check_square(function, "scale parameter", S);
+  check_symmetric(function, "scale parameter", S);
 
   MatrixXd S_inv = MatrixXd::Identity(k, k);
   S_inv = S.ldlt().solve(S_inv);

stan/math/rev/fun.hpp

@@ -32,7 +32,6 @@
 #include <stan/math/rev/fun/conj.hpp>
 #include <stan/math/rev/fun/cos.hpp>
 #include <stan/math/rev/fun/cosh.hpp>
-#include <stan/math/rev/fun/cov_exp_quad.hpp>
 #include <stan/math/rev/fun/determinant.hpp>
 #include <stan/math/rev/fun/digamma.hpp>
 #include <stan/math/rev/fun/dims.hpp>
@@ -58,6 +57,7 @@
 #include <stan/math/rev/fun/from_var_value.hpp>
 #include <stan/math/rev/fun/gamma_p.hpp>
 #include <stan/math/rev/fun/gamma_q.hpp>
+#include <stan/math/rev/fun/gp_exp_quad_cov.hpp>
 #include <stan/math/rev/fun/gp_periodic_cov.hpp>
 #include <stan/math/rev/fun/grad.hpp>
 #include <stan/math/rev/fun/grad_inc_beta.hpp>
@@ -150,6 +150,7 @@
 #include <stan/math/rev/fun/trigamma.hpp>
 #include <stan/math/rev/fun/trunc.hpp>
 #include <stan/math/rev/fun/typedefs.hpp>
+#include <stan/math/rev/fun/unit_vector_constrain.hpp>
 #include <stan/math/rev/fun/value_of.hpp>
 #include <stan/math/rev/fun/value_of_rec.hpp>
 #include <stan/math/rev/fun/variance.hpp>

stan/math/rev/fun/LDLT_alloc.hpp

@@ -36,7 +36,7 @@ class LDLT_alloc : public chainable_alloc {
    **/
   inline void compute(const Eigen::Matrix<var, R, C> &A) {
     N_ = A.rows();
-    variA_ = A.vi();
+    arena_A_ = A;
     ldlt_.compute(A.val());
   }
 
@@ -46,8 +46,8 @@ class LDLT_alloc : public chainable_alloc {
   }
 
   size_t N_;
-  Eigen::LDLT<Eigen::Matrix<double, R, C> > ldlt_;
-  Eigen::Matrix<vari *, R, C> variA_;
+  Eigen::LDLT<Eigen::Matrix<double, R, C>> ldlt_;
+  arena_matrix<Eigen::Matrix<var, R, C>> arena_A_;
 };
 
 }  // namespace math