Let prim functions return expressions

stan/math/fwd.hpp

@@ -5,9 +5,9 @@
 
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/meta.hpp>
-#include <stan/math/prim.hpp>
-
 #include <stan/math/fwd/fun.hpp>
 #include <stan/math/fwd/functor.hpp>
 
+#include <stan/math/prim.hpp>
+
 #endif

stan/math/fwd/fun/sum.hpp

@@ -1,9 +1,10 @@
 #ifndef STAN_MATH_FWD_FUN_SUM_HPP
 #define STAN_MATH_FWD_FUN_SUM_HPP
 
-#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/sum.hpp>
+#include <stan/math/fwd/core.hpp>
 #include <vector>
 
 namespace stan {

stan/math/mix.hpp

@@ -2,11 +2,19 @@
 #define STAN_MATH_MIX_HPP
 
 #include <stan/math/mix/meta.hpp>
-
-#include <stan/math/prim.hpp>
-#include <stan/math/fwd.hpp>
-#include <stan/math/rev.hpp>
 #include <stan/math/mix/fun.hpp>
 #include <stan/math/mix/functor.hpp>
 
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun.hpp>
+#include <stan/math/rev/functor.hpp>
+
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/fun.hpp>
+#include <stan/math/fwd/functor.hpp>
+
+#include <stan/math/prim.hpp>
+
 #endif

stan/math/prim/err/elementwise_check.hpp

@@ -142,8 +142,9 @@ inline void elementwise_check(const F& is_good, const char* function,
  */
 template <typename F, typename T, typename... Indexings,
           require_eigen_t<T>* = nullptr,
-          std::enable_if_t<static_cast<bool>(Eigen::internal::traits<T>::Flags&(
-              Eigen::LinearAccessBit | Eigen::DirectAccessBit))>* = nullptr>
+          std::enable_if_t<(Eigen::internal::traits<T>::Flags
+                            & Eigen::LinearAccessBit)
+                           || T::IsVectorAtCompileTime>* = nullptr>
 inline void elementwise_check(const F& is_good, const char* function,
                               const char* name, const T& x, const char* must_be,
                               const Indexings&... indexings) {
@@ -189,13 +190,13 @@ inline void elementwise_check(const F& is_good, const char* function,
  * @throws `std::domain_error` if `is_good` returns `false` for the value
  * of any element in `x`
  */
-template <
-    typename F, typename T, typename... Indexings,
-    require_eigen_t<T>* = nullptr,
-    std::enable_if_t<!(Eigen::internal::traits<T>::Flags
-                       & (Eigen::LinearAccessBit | Eigen::DirectAccessBit))
-                     && !(Eigen::internal::traits<T>::Flags
-                          & Eigen::RowMajorBit)>* = nullptr>
+template <typename F, typename T, typename... Indexings,
+          require_eigen_t<T>* = nullptr,
+          std::enable_if_t<!(Eigen::internal::traits<T>::Flags
+                             & Eigen::LinearAccessBit)
+                           && !T::IsVectorAtCompileTime
+                           && !(Eigen::internal::traits<T>::Flags
+                                & Eigen::RowMajorBit)>* = nullptr>
 inline void elementwise_check(const F& is_good, const char* function,
                               const char* name, const T& x, const char* must_be,
                               const Indexings&... indexings) {
@@ -230,13 +231,13 @@ inline void elementwise_check(const F& is_good, const char* function,
  * @throws `std::domain_error` if `is_good` returns `false` for the value
  * of any element in `x`
  */
-template <
-    typename F, typename T, typename... Indexings,
-    require_eigen_t<T>* = nullptr,
-    std::enable_if_t<!(Eigen::internal::traits<T>::Flags
-                       & (Eigen::LinearAccessBit | Eigen::DirectAccessBit))
-                     && static_cast<bool>(Eigen::internal::traits<T>::Flags
-                                          & Eigen::RowMajorBit)>* = nullptr>
+template <typename F, typename T, typename... Indexings,
+          require_eigen_t<T>* = nullptr,
+          std::enable_if_t<
+              !(Eigen::internal::traits<T>::Flags & Eigen::LinearAccessBit)
+              && !T::IsVectorAtCompileTime
+              && static_cast<bool>(Eigen::internal::traits<T>::Flags
+                                   & Eigen::RowMajorBit)>* = nullptr>
 inline void elementwise_check(const F& is_good, const char* function,
                               const char* name, const T& x, const char* must_be,
                               const Indexings&... indexings) {

stan/math/prim/fun/add.hpp

@@ -42,7 +42,7 @@ template <typename Mat1, typename Mat2,
           require_all_not_st_var<Mat1, Mat2>* = nullptr>
 inline auto add(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("add", "m1", m1, "m2", m2);
-  return (m1 + m2).eval();
+  return m1 + m2;
 }
 
 /**
@@ -58,7 +58,7 @@ 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 add(const Mat& m, const Scal c) {
-  return (m.array() + c).matrix().eval();
+  return (m.array() + c).matrix();
 }
 
 /**
@@ -74,7 +74,7 @@ template <typename Scal, typename Mat, require_stan_scalar_t<Scal>* = nullptr,
           require_eigen_t<Mat>* = nullptr,
           require_all_not_st_var<Scal, Mat>* = nullptr>
 inline auto add(const Scal c, const Mat& m) {
-  return (c + m.array()).matrix().eval();
+  return (c + m.array()).matrix();
 }
 
 }  // namespace math

stan/math/prim/fun/block.hpp

@@ -19,8 +19,7 @@ namespace math {
  * @throw std::out_of_range if either index is out of range.
  */
 template <typename T, require_eigen_t<T>* = nullptr>
-inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> block(
-    const T& m, size_t i, size_t j, size_t nrows, size_t ncols) {
+inline auto block(const T& m, size_t i, size_t j, size_t nrows, size_t ncols) {
   check_row_index("block", "i", m, i);
   check_row_index("block", "i+nrows-1", m, i + nrows - 1);
   check_column_index("block", "j", m, j);

stan/math/prim/fun/col.hpp

@@ -23,7 +23,7 @@ namespace math {
 template <typename T, typename = require_eigen_t<T>>
 inline auto col(const T& m, size_t j) {
   check_column_index("col", "j", m, j);
-  return m.col(j - 1).eval();
+  return m.col(j - 1);
 }
 
 }  // namespace math

stan/math/prim/fun/corr_matrix_free.hpp

@@ -40,16 +40,16 @@ Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, 1> corr_matrix_free(const T& y) {
 
   Eigen::Index k = y.rows();
   Eigen::Index k_choose_2 = (k * (k - 1)) / 2;
-  Array<value_type_t<T>, Dynamic, 1> x(k_choose_2);
+  Eigen::Matrix<value_type_t<T>, Dynamic, 1> x(k_choose_2);
   Array<value_type_t<T>, Dynamic, 1> sds(k);
-  bool successful = factor_cov_matrix(y, x, sds);
+  bool successful = factor_cov_matrix(y, x.array(), sds);
   if (!successful) {
     throw_domain_error("corr_matrix_free", "factor_cov_matrix failed on y", y,
                        "");
   }
   check_bounded("corr_matrix_free", "log(sd)", sds, -CONSTRAINT_TOLERANCE,
                 CONSTRAINT_TOLERANCE);
-  return x.matrix();
+  return x;
 }
 
 }  // namespace math

stan/math/prim/fun/diag_post_multiply.hpp

@@ -12,7 +12,7 @@ template <typename T1, typename T2, require_eigen_t<T1>* = nullptr,
 auto diag_post_multiply(const T1& m1, const T2& m2) {
   check_size_match("diag_post_multiply", "m2.size()", m2.size(), "m1.cols()",
                    m1.cols());
-  return (m1 * m2.asDiagonal()).eval();
+  return m1 * m2.asDiagonal();
 }
 
 }  // namespace math

stan/math/prim/fun/diag_pre_multiply.hpp

@@ -12,7 +12,7 @@ template <typename T1, typename T2, require_eigen_vector_t<T1>* = nullptr,
 auto diag_pre_multiply(const T1& m1, const T2& m2) {
   check_size_match("diag_pre_multiply", "m1.size()", m1.size(), "m2.rows()",
                    m2.rows());
-  return (m1.asDiagonal() * m2).eval();
+  return m1.asDiagonal() * m2;
 }
 
 }  // namespace math

stan/math/prim/fun/diagonal.hpp

@@ -17,7 +17,7 @@ namespace math {
  */
 template <typename T, typename = require_eigen_t<T>>
 inline auto diagonal(const T& m) {
-  return m.diagonal().eval();
+  return m.diagonal();
 }
 
 }  // namespace math

stan/math/prim/fun/divide.hpp

@@ -19,7 +19,7 @@ namespace math {
  * @return Scalar divided by the scalar.
  */
 template <typename Scal1, typename Scal2,
-          typename = require_all_stan_scalar_t<Scal1, Scal2>>
+          require_all_stan_scalar_t<Scal1, Scal2>* = nullptr>
 inline return_type_t<Scal1, Scal2> divide(const Scal1& x, const Scal2& y) {
   return x / y;
 }
@@ -41,10 +41,10 @@ inline int divide(int x, int y) {
  * @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>>
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_all_not_var_t<scalar_type_t<Mat>, Scal>* = nullptr>
 inline auto divide(const Mat& m, Scal c) {
-  return (m / c).eval();
+  return m / c;
 }
 
 }  // namespace math

stan/math/prim/fun/elt_divide.hpp

@@ -24,7 +24,7 @@ template <typename Mat1, typename Mat2,
           require_all_not_st_var<Mat1, Mat2>* = nullptr>
 auto elt_divide(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("elt_divide", "m1", m1, "m2", m2);
-  return (m1.array() / m2.array()).matrix().eval();
+  return (m1.array() / m2.array()).matrix();
 }
 
 /**
@@ -38,8 +38,8 @@ auto elt_divide(const Mat1& m1, const Mat2& m2) {
  * @param s scalar
  * @return Elementwise division of a scalar by matrix.
  */
-template <typename Mat, typename Scal, typename = require_matrix_t<Mat>,
-          typename = require_stan_scalar_t<Scal>>
+template <typename Mat, typename Scal, require_matrix_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr>
 auto elt_divide(const Mat& m, Scal s) {
   return divide(m, s);
 }
@@ -55,10 +55,10 @@ auto elt_divide(const Mat& m, Scal s) {
  * @param m matrix or expression
  * @return Elementwise division of a scalar by matrix.
  */
-template <typename Scal, typename Mat, typename = require_stan_scalar_t<Scal>,
-          typename = require_eigen_t<Mat>>
+template <typename Scal, typename Mat, require_stan_scalar_t<Scal>* = nullptr,
+          require_eigen_t<Mat>* = nullptr>
 auto elt_divide(Scal s, const Mat& m) {
-  return (s / m.array()).matrix().eval();
+  return (s / m.array()).matrix();
 }
 
 template <typename Scal1, typename Scal2,

stan/math/prim/fun/elt_multiply.hpp

@@ -24,7 +24,7 @@ template <typename Mat1, typename Mat2,
           require_all_not_st_var<Mat1, Mat2>* = nullptr>
 auto elt_multiply(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("elt_multiply", "m1", m1, "m2", m2);
-  return m1.cwiseProduct(m2).eval();
+  return m1.cwiseProduct(m2);
 }
 
 /**

stan/math/prim/fun/head.hpp

@@ -23,7 +23,7 @@ inline auto head(const T& v, size_t n) {
   if (n != 0) {
     check_vector_index("head", "n", v, n);
   }
-  return v.head(n).eval();
+  return v.head(n);
 }
 
 /**

stan/math/prim/fun/identity_matrix.hpp

@@ -14,7 +14,7 @@ namespace math {
  * @return An identity matrix of size K.
  * @throw std::domain_error if K is negative.
  */
-inline Eigen::MatrixXd identity_matrix(int K) {
+inline auto identity_matrix(int K) {
   check_nonnegative("identity_matrix", "size", K);
   return Eigen::MatrixXd::Identity(K, K);
 }

stan/math/prim/fun/inv_square.hpp

@@ -4,8 +4,6 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/square.hpp>
-#include <stan/math/prim/fun/inv_square.hpp>
-#include <stan/math/prim/functor/apply_scalar_unary.hpp>
 #include <stan/math/prim/functor/apply_vector_unary.hpp>
 
 namespace stan {

stan/math/prim/fun/linspaced_row_vector.hpp

@@ -22,7 +22,7 @@ namespace math {
  * @throw std::domain_error if K is negative, if low is nan or infinite,
  * if high is nan or infinite, or if high is less than low.
  */
-inline Eigen::RowVectorXd linspaced_row_vector(int K, double low, double high) {
+inline auto linspaced_row_vector(int K, double low, double high) {
   static const char* function = "linspaced_row_vector";
   check_nonnegative(function, "size", K);
   check_finite(function, "low", low);

stan/math/prim/fun/linspaced_vector.hpp

@@ -22,7 +22,7 @@ namespace math {
  * @throw std::domain_error if K is negative, if low is nan or infinite,
  * if high is nan or infinite, or if high is less than low.
  */
-inline Eigen::VectorXd linspaced_vector(int K, double low, double high) {
+inline auto linspaced_vector(int K, double low, double high) {
   static const char* function = "linspaced_vector";
   check_nonnegative(function, "size", K);
   check_finite(function, "low", low);

stan/math/prim/fun/log1p.hpp

@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/is_nan.hpp>
-#include <stan/math/prim/fun/log1p.hpp>
 #include <stan/math/prim/functor/apply_scalar_unary.hpp>
 #include <cmath>
 

stan/math/prim/fun/log_softmax.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/functor/apply_vector_unary.hpp>
 
 namespace stan {
@@ -37,18 +38,16 @@ namespace math {
  * @tparam Container type of input vector to transform
  * @param[in] x vector to transform
  * @return log unit simplex result of the softmax transform of the vector.
- *
- * Note: The return must be evaluated otherwise the Ref object falls out
- * of scope
  */
-template <typename Container,
-          require_arithmetic_t<scalar_type_t<Container>>* = nullptr>
+template <typename Container, require_st_arithmetic<Container>* = nullptr>
 inline auto log_softmax(const Container& x) {
-  return apply_vector_unary<Container>::apply(x, [](const auto& v) {
-    const Eigen::Ref<const plain_type_t<decltype(v)>>& v_ref = v;
-    check_nonzero_size("log_softmax", "v", v_ref);
-    return (v_ref.array() - log_sum_exp(v_ref)).eval();
-  });
+  check_nonzero_size("log_softmax", "v", x);
+  return make_holder(
+      [](const auto& a) {
+        return apply_vector_unary<ref_type_t<Container>>::apply(
+            a, [](const auto& v) { return v.array() - log_sum_exp(v); });
+      },
+      to_ref(x));
 }
 
 }  // namespace math

stan/math/prim/fun/logit.hpp

@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/log.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/functor/apply_scalar_unary.hpp>
 #include <stan/math/prim/functor/apply_vector_unary.hpp>
 #include <cmath>
@@ -102,10 +103,13 @@ inline auto logit(const Container& x) {
 template <typename Container,
           require_container_st<std::is_arithmetic, Container>* = nullptr>
 inline auto logit(const Container& x) {
-  return apply_vector_unary<Container>::apply(x, [](const auto& v) {
-    const Eigen::Ref<const plain_type_t<decltype(v)>>& v_ref = v;
-    return (v_ref.array() / (1 - v_ref.array())).log().eval();
-  });
+  return make_holder(
+      [](const auto& v_ref) {
+        return apply_vector_unary<ref_type_t<Container>>::apply(
+            v_ref,
+            [](const auto& v) { return (v.array() / (1 - v.array())).log(); });
+      },
+      to_ref(x));
 }
 
 }  // namespace math

stan/math/prim/fun/minus.hpp

@@ -14,7 +14,7 @@ namespace math {
  * @return Negation of subtrahend.
  */
 template <typename T>
-inline plain_type_t<T> minus(const T& x) {
+inline auto minus(const T& x) {
   return -x;
 }