Use perfect forwarding for functions that use apply_*_* functions

stan/math/fwd/fun/Phi.hpp

@@ -12,10 +12,10 @@
 namespace stan {
 namespace math {
 
-template <typename T>
-inline fvar<T> Phi(const fvar<T>& x) {
-  T xv = x.val_;
-  return fvar<T>(Phi(xv), x.d_ * exp(xv * xv / -2.0) * INV_SQRT_TWO_PI);
+template <typename T, require_fvar_t<T>* = nullptr>
+inline auto Phi(T&& x) {
+  return std::decay_t<T>(Phi(x.val_),
+                         x.d_ * exp(x.val_ * x.val_ / -2.0) * INV_SQRT_TWO_PI);
 }
 
 }  // namespace math

stan/math/fwd/fun/cbrt.hpp

@@ -16,9 +16,9 @@ namespace math {
  * @param x Argument.
  * @return Cube root of argument.
  */
-template <typename T>
-inline fvar<T> cbrt(const fvar<T>& x) {
-  return fvar<T>(cbrt(x.val_), x.d_ / (3 * square(cbrt(x.val_))));
+template <typename T, require_fvar_t<T>* = nullptr>
+inline auto cbrt(T&& x) {
+  return std::decay_t<T>(cbrt(x.val_), x.d_ / (3 * square(cbrt(x.val_))));
 }
 
 }  // namespace math

stan/math/fwd/fun/digamma.hpp

@@ -4,7 +4,6 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/fun/trigamma.hpp>
-
 #include <stan/math/prim/fun/digamma.hpp>
 
 namespace stan {
@@ -19,9 +18,9 @@ namespace math {
  * @return derivative of the log gamma function at the specified
  * argument
  */
-template <typename T>
-inline fvar<T> digamma(const fvar<T>& x) {
-  return fvar<T>(digamma(x.val_), x.d_ * trigamma(x.val_));
+template <typename T, require_fvar_t<T>* = nullptr>
+inline auto digamma(T&& x) {
+  return std::decay_t<T>(digamma(x.val_), x.d_ * trigamma(x.val_));
 }
 
 }  // namespace math

stan/math/fwd/fun/log_softmax.hpp

@@ -6,6 +6,7 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun/softmax.hpp>
 #include <stan/math/prim/fun/log_softmax.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 
 namespace stan {
 namespace math {
@@ -19,24 +20,24 @@ namespace math {
  * @throw std::domain_error If the input vector is size 0.
  */
 template <typename T, require_vector_st<is_fvar, T>* = nullptr>
-inline auto log_softmax(const T& x) {
-  return apply_vector_unary<T>::apply(x, [&](const auto& alpha) {
+inline auto log_softmax(T&& x) {
+  return apply_vector_unary<T>::apply(std::forward<T>(x), [](auto&& alpha) {
     using T_alpha = decltype(alpha);
     using T_fvar = value_type_t<T_alpha>;
     using T_fvar_inner = typename T_fvar::Scalar;
 
-    const Eigen::Ref<const plain_type_t<T_alpha>>& alpha_ref = alpha;
-    Eigen::Matrix<T_fvar_inner, -1, 1> alpha_t = alpha_ref.val();
-    Eigen::Matrix<T_fvar_inner, -1, 1> softmax_alpha_t = softmax(alpha_t);
+    auto&& alpha_ref = to_ref(std::forward<decltype(alpha)>(alpha));
+    Eigen::Matrix<T_fvar_inner, -1, 1> alpha_v = alpha_ref.val();
+    Eigen::Matrix<T_fvar_inner, -1, 1> softmax_alpha_v = softmax(alpha_v);
 
-    Eigen::Matrix<T_fvar, -1, 1> log_softmax_alpha(alpha.size());
-    log_softmax_alpha.val() = log_softmax(alpha_t);
+    Eigen::Matrix<T_fvar, -1, 1> log_softmax_alpha(alpha_ref.size());
+    log_softmax_alpha.val() = log_softmax(alpha_v);
     log_softmax_alpha.d().setZero();
 
-    for (int m = 0; m < alpha.size(); ++m) {
+    for (int m = 0; m < alpha_ref.size(); ++m) {
       T_fvar_inner negative_alpha_m_d_times_softmax_alpha_t_m
-          = -alpha_ref.coeff(m).d_ * softmax_alpha_t(m);
-      for (int k = 0; k < alpha.size(); ++k) {
+          = -alpha_ref.coeff(m).d_ * softmax_alpha_v(m);
+      for (int k = 0; k < alpha_ref.size(); ++k) {
         if (m == k) {
           log_softmax_alpha(k).d_
               += alpha_ref.coeff(m).d_

stan/math/fwd/fun/log_sum_exp.hpp

@@ -50,9 +50,9 @@ inline fvar<T> log_sum_exp(const fvar<T>& x1, double x2) {
  * @return The log of the sum of the exponentiated vector values.
  */
 template <typename T, require_container_st<is_fvar, T>* = nullptr>
-inline auto log_sum_exp(const T& x) {
+inline auto log_sum_exp(T&& x) {
   return apply_vector_unary<ref_type_t<T>>::reduce(
-      to_ref(x), [&](const auto& v) {
+      to_ref(std::forward<T>(x)), [&](const auto& v) {
         using T_fvar_inner = typename value_type_t<decltype(v)>::Scalar;
         using mat_type = Eigen::Matrix<T_fvar_inner, -1, -1>;
         mat_type vals = v.val();

stan/math/fwd/fun/norm1.hpp

@@ -20,9 +20,9 @@ namespace math {
  * @return L1 norm of x.
  */
 template <typename Container, require_eigen_vt<is_fvar, Container>* = nullptr>
-inline auto norm1(const Container& x) {
+inline auto norm1(Container&& x) {
   return apply_vector_unary<ref_type_t<Container>>::reduce(
-      to_ref(x), [&](const auto& v) {
+      to_ref(std::forward<Container>(x)), [&](const auto& v) {
         using T_fvar_inner = typename value_type_t<decltype(v)>::Scalar;
         return fvar<T_fvar_inner>(norm1(v.val()),
                                   v.d().cwiseProduct(sign(v.val())).sum());

stan/math/fwd/fun/norm2.hpp

@@ -19,9 +19,9 @@ namespace math {
  * @return L2 norm of x.
  */
 template <typename Container, require_eigen_vt<is_fvar, Container>* = nullptr>
-inline auto norm2(const Container& x) {
+inline auto norm2(Container&& x) {
   return apply_vector_unary<ref_type_t<Container>>::reduce(
-      to_ref(x), [&](const auto& v) {
+      to_ref(std::forward<Container>(x)), [&](const auto& v) {
         using T_fvar_inner = typename value_type_t<decltype(v)>::Scalar;
         T_fvar_inner res = norm2(v.val());
         return fvar<T_fvar_inner>(res,

stan/math/fwd/fun/pow.hpp

@@ -73,9 +73,13 @@ inline auto pow(const T1& x1, const T2& x2) {
 template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr,
           require_all_not_matrix_st<is_var, T1, T2>* = nullptr,
           require_any_fvar_t<base_type_t<T1>, base_type_t<T2>>* = nullptr>
-inline auto pow(const T1& a, const T2& b) {
+inline auto pow(T1&& a, T2&& b) {
   return apply_scalar_binary(
-      [](const auto& c, const auto& d) { return stan::math::pow(c, d); }, a, b);
+      [](auto&& c, auto&& d) {
+        return stan::math::pow(std::forward<decltype(c)>(c),
+                               std::forward<decltype(d)>(d));
+      },
+      std::forward<T1>(a), std::forward<T2>(b));
 }
 
 }  // namespace math

stan/math/fwd/functor/apply_scalar_unary.hpp

@@ -18,20 +18,23 @@ namespace math {
  * autodiff variable.
  */
 template <typename F, typename T>
-struct apply_scalar_unary<F, fvar<T> > {
+struct apply_scalar_unary<F, T, require_fvar_t<T>> {
   /**
    * Function return type, which is same as the argument type for
    * the function, <code>fvar&lt;T&gt;</code>.
    */
-  using return_t = fvar<T>;
+  using return_t = std::decay_t<T>;
 
   /**
    * Apply the function specified by F to the specified argument.
    *
    * @param x Argument variable.
    * @return Function applied to the variable.
    */
-  static inline return_t apply(const fvar<T>& x) { return F::fun(x); }
+  template <typename T2>
+  static inline auto apply(T2&& x) {
+    return F::fun(std::forward<T2>(x));
+  }
 };
 
 }  // namespace math

stan/math/prim/constraint/cholesky_corr_constrain.hpp

@@ -14,15 +14,15 @@ namespace math {
 
 template <typename EigVec, require_eigen_col_vector_t<EigVec>* = nullptr>
 inline Eigen::Matrix<value_type_t<EigVec>, Eigen::Dynamic, Eigen::Dynamic>
-cholesky_corr_constrain(const EigVec& y, int K) {
+cholesky_corr_constrain(EigVec&& y, int K) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::sqrt;
   using T_scalar = value_type_t<EigVec>;
   int k_choose_2 = (K * (K - 1)) / 2;
   check_size_match("cholesky_corr_constrain", "constrain size", y.size(),
                    "k_choose_2", k_choose_2);
-  Matrix<T_scalar, Dynamic, 1> z = corr_constrain(y);
+  Matrix<T_scalar, Dynamic, 1> z = corr_constrain(std::forward<EigVec>(y));
   Matrix<T_scalar, Dynamic, Dynamic> x(K, K);
   if (K == 0) {
     return x;
@@ -47,15 +47,15 @@ template <typename EigVec, typename Lp,
           require_eigen_vector_t<EigVec>* = nullptr,
           require_convertible_t<return_type_t<EigVec>, Lp>* = nullptr>
 inline Eigen::Matrix<value_type_t<EigVec>, Eigen::Dynamic, Eigen::Dynamic>
-cholesky_corr_constrain(const EigVec& y, int K, Lp& lp) {
+cholesky_corr_constrain(EigVec&& y, int K, Lp& lp) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::sqrt;
   using T_scalar = value_type_t<EigVec>;
   int k_choose_2 = (K * (K - 1)) / 2;
   check_size_match("cholesky_corr_constrain", "y.size()", y.size(),
                    "k_choose_2", k_choose_2);
-  Matrix<T_scalar, Dynamic, 1> z = corr_constrain(y, lp);
+  Matrix<T_scalar, Dynamic, 1> z = corr_constrain(std::forward<EigVec>(y), lp);
   Matrix<T_scalar, Dynamic, Dynamic> x(K, K);
   if (K == 0) {
     return x;
@@ -87,9 +87,10 @@ cholesky_corr_constrain(const EigVec& y, int K, Lp& lp) {
  * @param K The size of the matrix to return
  */
 template <typename T, require_std_vector_t<T>* = nullptr>
-inline auto cholesky_corr_constrain(const T& y, int K) {
-  return apply_vector_unary<T>::apply(
-      y, [K](auto&& v) { return cholesky_corr_constrain(v, K); });
+inline auto cholesky_corr_constrain(T&& y, int K) {
+  return apply_vector_unary<T>::apply(std::forward<T>(y), [K](auto&& v) {
+    return cholesky_corr_constrain(std::forward<decltype(v)>(v), K);
+  });
 }
 
 /**
@@ -107,9 +108,10 @@ inline auto cholesky_corr_constrain(const T& y, int K) {
  */
 template <typename T, typename Lp, require_std_vector_t<T>* = nullptr,
           require_convertible_t<return_type_t<T>, Lp>* = nullptr>
-inline auto cholesky_corr_constrain(const T& y, int K, Lp& lp) {
-  return apply_vector_unary<T>::apply(
-      y, [&lp, K](auto&& v) { return cholesky_corr_constrain(v, K, lp); });
+inline auto cholesky_corr_constrain(T&& y, int K, Lp& lp) {
+  return apply_vector_unary<T>::apply(std::forward<T>(y), [&lp, K](auto&& v) {
+    return cholesky_corr_constrain(std::forward<decltype(v)>(v), K, lp);
+  });
 }
 
 /**
@@ -132,11 +134,11 @@ inline auto cholesky_corr_constrain(const T& y, int K, Lp& lp) {
  */
 template <bool Jacobian, typename T, typename Lp,
           require_convertible_t<return_type_t<T>, Lp>* = nullptr>
-inline auto cholesky_corr_constrain(const T& y, int K, Lp& lp) {
+inline auto cholesky_corr_constrain(T&& y, int K, Lp& lp) {
   if constexpr (Jacobian) {
-    return cholesky_corr_constrain(y, K, lp);
+    return cholesky_corr_constrain(std::forward<T>(y), K, lp);
   } else {
-    return cholesky_corr_constrain(y, K);
+    return cholesky_corr_constrain(std::forward<T>(y), K);
   }
 }
 

stan/math/prim/constraint/cholesky_corr_free.hpp

@@ -11,19 +11,19 @@ namespace stan {
 namespace math {
 
 template <typename T, require_eigen_t<T>* = nullptr>
-inline auto cholesky_corr_free(const T& x) {
+inline auto cholesky_corr_free(T&& x) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::sqrt;
 
   check_square("cholesky_corr_free", "x", x);
   // should validate lower-triangular, unit lengths
 
-  const auto& x_ref = to_ref(x);
-  int K = (x.rows() * (x.rows() - 1)) / 2;
+  auto&& x_ref = to_ref(std::forward<T>(x));
+  int K = (x_ref.rows() * (x_ref.rows() - 1)) / 2;
   Matrix<value_type_t<T>, Dynamic, 1> z(K);
   int k = 0;
-  for (int i = 1; i < x.rows(); ++i) {
+  for (int i = 1; i < x_ref.rows(); ++i) {
     z.coeffRef(k++) = corr_free(x_ref.coeff(i, 0));
     auto sum_sqs = square(x_ref.coeff(i, 0));
     for (int j = 1; j < i; ++j) {
@@ -42,9 +42,10 @@ inline auto cholesky_corr_free(const T& x) {
  * @param x The standard vector to untransform.
  */
 template <typename T, require_std_vector_t<T>* = nullptr>
-inline auto cholesky_corr_free(const T& x) {
-  return apply_vector_unary<T>::apply(
-      x, [](auto&& v) { return cholesky_corr_free(v); });
+inline auto cholesky_corr_free(T&& x) {
+  return apply_vector_unary<T>::apply(std::forward<T>(x), [](auto&& v) {
+    return cholesky_corr_free(std::forward<decltype(v)>(v));
+  });
 }
 
 }  // namespace math

stan/math/prim/constraint/cholesky_factor_constrain.hpp

@@ -27,7 +27,7 @@ namespace math {
  */
 template <typename T, require_eigen_col_vector_t<T>* = nullptr>
 inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
-cholesky_factor_constrain(const T& x, int M, int N) {
+cholesky_factor_constrain(T&& x, int M, int N) {
   using std::exp;
   using T_scalar = value_type_t<T>;
   check_greater_or_equal("cholesky_factor_constrain",
@@ -39,7 +39,7 @@ cholesky_factor_constrain(const T& x, int M, int N) {
   Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic> y(M, N);
   int pos = 0;
 
-  const auto& x_ref = to_ref(x);
+  auto&& x_ref = to_ref(std::forward<T>(x));
   for (int m = 0; m < N; ++m) {
     y.row(m).head(m) = x_ref.segment(pos, m);
     pos += m;
@@ -73,17 +73,17 @@ cholesky_factor_constrain(const T& x, int M, int N) {
 template <typename T, typename Lp, require_eigen_vector_t<T>* = nullptr,
           require_convertible_t<return_type_t<T>, Lp>* = nullptr>
 inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
-cholesky_factor_constrain(const T& x, int M, int N, Lp& lp) {
+cholesky_factor_constrain(T&& x, int M, int N, Lp& lp) {
   check_size_match("cholesky_factor_constrain", "x.size()", x.size(),
                    "((N * (N + 1)) / 2 + (M - N) * N)",
                    ((N * (N + 1)) / 2 + (M - N) * N));
   int pos = 0;
-  const auto& x_ref = to_ref(x);
+  auto&& x_ref = to_ref(std::forward<T>(x));
   for (int n = 0; n < N; ++n) {
     pos += n;
     lp += x_ref.coeff(pos++);
   }
-  return cholesky_factor_constrain(x_ref, M, N);
+  return cholesky_factor_constrain(std::forward<decltype(x_ref)>(x_ref), M, N);
 }
 
 /**
@@ -101,9 +101,10 @@ cholesky_factor_constrain(const T& x, int M, int N, Lp& lp) {
  * @return Cholesky factor
  */
 template <typename T, require_std_vector_t<T>* = nullptr>
-inline auto cholesky_factor_constrain(const T& x, int M, int N) {
-  return apply_vector_unary<T>::apply(
-      x, [M, N](auto&& v) { return cholesky_factor_constrain(v, M, N); });
+inline auto cholesky_factor_constrain(T&& x, int M, int N) {
+  return apply_vector_unary<T>::apply(std::forward<T>(x), [M, N](auto&& v) {
+    return cholesky_factor_constrain(std::forward<decltype(v)>(v), M, N);
+  });
 }
 
 /**
@@ -125,9 +126,10 @@ inline auto cholesky_factor_constrain(const T& x, int M, int N) {
  */
 template <typename T, typename Lp, require_std_vector_t<T>* = nullptr,
           require_convertible_t<return_type_t<T>, Lp>* = nullptr>
-inline auto cholesky_factor_constrain(const T& x, int M, int N, Lp& lp) {
-  return apply_vector_unary<T>::apply(x, [&lp, M, N](auto&& v) {
-    return cholesky_factor_constrain(v, M, N, lp);
+inline auto cholesky_factor_constrain(T&& x, int M, int N, Lp& lp) {
+  return apply_vector_unary<T>::apply(std::forward<T>(x), [&lp, M,
+                                                           N](auto&& v) {
+    return cholesky_factor_constrain(std::forward<decltype(v)>(v), M, N, lp);
   });
 }
 
@@ -155,11 +157,11 @@ inline auto cholesky_factor_constrain(const T& x, int M, int N, Lp& lp) {
  */
 template <bool Jacobian, typename T, typename Lp,
           require_convertible_t<return_type_t<T>, Lp>* = nullptr>
-inline auto cholesky_factor_constrain(const T& x, int M, int N, Lp& lp) {
+inline auto cholesky_factor_constrain(T&& x, int M, int N, Lp& lp) {
   if constexpr (Jacobian) {
-    return cholesky_factor_constrain(x, M, N, lp);
+    return cholesky_factor_constrain(std::forward<T>(x), M, N, lp);
   } else {
-    return cholesky_factor_constrain(x, M, N);
+    return cholesky_factor_constrain(std::forward<T>(x), M, N);
   }
 }