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

stan/math/fwd/fun/log_softmax.hpp

@@ -19,36 +19,38 @@ 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) {
-    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);
-
-    Eigen::Matrix<T_fvar, -1, 1> log_softmax_alpha(alpha.size());
-    log_softmax_alpha.val() = log_softmax(alpha_t);
-    log_softmax_alpha.d().setZero();
-
-    for (int m = 0; m < alpha.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) {
-        if (m == k) {
-          log_softmax_alpha(k).d_
-              += alpha_ref.coeff(m).d_
-                 + negative_alpha_m_d_times_softmax_alpha_t_m;
-        } else {
-          log_softmax_alpha(k).d_ += negative_alpha_m_d_times_softmax_alpha_t_m;
+inline auto log_softmax(T&& x) {
+  return apply_vector_unary<T>::apply(
+      std::forward<T>(x), [&](const 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);
+
+        Eigen::Matrix<T_fvar, -1, 1> log_softmax_alpha(alpha.size());
+        log_softmax_alpha.val() = log_softmax(alpha_t);
+        log_softmax_alpha.d().setZero();
+
+        for (int m = 0; m < alpha.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) {
+            if (m == k) {
+              log_softmax_alpha(k).d_
+                  += alpha_ref.coeff(m).d_
+                     + negative_alpha_m_d_times_softmax_alpha_t_m;
+            } else {
+              log_softmax_alpha(k).d_
+                  += negative_alpha_m_d_times_softmax_alpha_t_m;
+            }
+          }
         }
-      }
-    }
 
-    return log_softmax_alpha;
-  });
+        return log_softmax_alpha;
+      });
 }
 
 }  // namespace math

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(const T2& x) {
+    return F::fun(x);
+  }
 };
 
 }  // namespace math

stan/math/prim/constraint/cholesky_corr_constrain.hpp

@@ -87,9 +87,11 @@ 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<std::decay_t<T>>::apply(
+      std::forward<T>(y), [K](auto&& v) {
+        return cholesky_corr_constrain(std::forward<decltype(v)>(v), K);
+      });
 }
 
 /**
@@ -107,9 +109,11 @@ 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<std::decay_t<T>>::apply(
+      std::forward<T>(y), [&lp, K](auto&& v) {
+        return cholesky_corr_constrain(std::forward<decltype(v)>(v), K, lp);
+      });
 }
 
 /**
@@ -132,11 +136,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

@@ -42,9 +42,11 @@ 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<std::decay_t<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

@@ -101,9 +101,11 @@ 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<std::decay_t<T>>::apply(
+      std::forward<T>(x), [M, N](auto&& v) {
+        return cholesky_factor_constrain(std::forward<decltype(v)>(v), M, N);
+      });
 }
 
 /**
@@ -125,10 +127,12 @@ 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<std::decay_t<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 +159,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);
   }
 }
 

stan/math/prim/constraint/cholesky_factor_free.hpp

@@ -56,9 +56,11 @@ Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, 1> cholesky_factor_free(
  * @param x The standard vector to untransform.
  */
 template <typename T, require_std_vector_t<T>* = nullptr>
-inline auto cholesky_factor_free(const T& x) {
-  return apply_vector_unary<T>::apply(
-      x, [](auto&& v) { return cholesky_factor_free(v); });
+inline auto cholesky_factor_free(T&& x) {
+  return apply_vector_unary<std::decay_t<T>>::apply(
+      std::forward<T>(x), [](auto&& v) {
+        return cholesky_factor_free(std::forward<decltype(v)>(v));
+      });
 }
 
 }  // namespace math

stan/math/prim/constraint/corr_constrain.hpp

@@ -24,8 +24,8 @@ namespace math {
  * @return tanh transform
  */
 template <typename T>
-inline plain_type_t<T> corr_constrain(const T& x) {
-  return tanh(x);
+inline plain_type_t<T> corr_constrain(T&& x) {
+  return tanh(std::forward<T>(x));
 }
 
 /**
@@ -43,7 +43,7 @@ inline plain_type_t<T> corr_constrain(const T& x) {
  * @param[in,out] lp log density accumulator
  */
 template <typename T_x, typename T_lp>
-inline auto corr_constrain(const T_x& x, T_lp& lp) {
+inline auto corr_constrain(T_x&& x, T_lp& lp) {
   plain_type_t<T_x> tanh_x = tanh(x);
   lp += sum(log1m(square(tanh_x)));
   return tanh_x;
@@ -65,11 +65,11 @@ inline auto corr_constrain(const T_x& x, T_lp& lp) {
  * @param[in,out] lp log density accumulator
  */
 template <bool Jacobian, typename T_x, typename T_lp>
-inline auto corr_constrain(const T_x& x, T_lp& lp) {
+inline auto corr_constrain(T_x&& x, T_lp& lp) {
   if constexpr (Jacobian) {
-    return corr_constrain(x, lp);
+    return corr_constrain(std::forward<T_x>(x), lp);
   } else {
-    return corr_constrain(x);
+    return corr_constrain(std::forward<T_x>(x));
   }
 }
 

stan/math/prim/constraint/corr_free.hpp

@@ -25,9 +25,9 @@ namespace math {
  * @return free scalar that transforms to the specified input
  */
 template <typename T>
-inline plain_type_t<T> corr_free(const T& y) {
+inline plain_type_t<T> corr_free(T&& y) {
   check_bounded("lub_free", "Correlation variable", y, -1.0, 1.0);
-  return atanh(y);
+  return atanh(std::forward<T>(y));
 }
 
 }  // namespace math

stan/math/prim/constraint/corr_matrix_constrain.hpp

@@ -38,11 +38,11 @@ namespace math {
  */
 template <typename T, require_eigen_col_vector_t<T>* = nullptr>
 inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
-corr_matrix_constrain(const T& x, Eigen::Index k) {
-  Eigen::Index k_choose_2 = (k * (k - 1)) / 2;
+corr_matrix_constrain(T&& x, Eigen::Index k) {
+  const Eigen::Index k_choose_2 = (k * (k - 1)) / 2;
   check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2",
                    k_choose_2);
-  return read_corr_matrix(corr_constrain(x), k);
+  return read_corr_matrix(corr_constrain(std::forward<T>(x)), k);
 }
 
 /**
@@ -70,11 +70,11 @@ corr_matrix_constrain(const T& x, Eigen::Index k) {
 template <typename T, typename Lp, require_eigen_col_vector_t<T>* = nullptr,
           require_convertible_t<return_type_t<T>, Lp>* = nullptr>
 inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
-corr_matrix_constrain(const T& x, Eigen::Index k, Lp& lp) {
-  Eigen::Index k_choose_2 = (k * (k - 1)) / 2;
+corr_matrix_constrain(T&& x, Eigen::Index k, Lp& lp) {
+  const Eigen::Index k_choose_2 = (k * (k - 1)) / 2;
   check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2",
                    k_choose_2);
-  return read_corr_matrix(corr_constrain(x, lp), k, lp);
+  return read_corr_matrix(corr_constrain(std::forward<T>(x), lp), k, lp);
 }
 
 /**
@@ -91,9 +91,11 @@ corr_matrix_constrain(const T& x, Eigen::Index k, Lp& lp) {
  * @param K Dimensionality of returned correlation matrix
  */
 template <typename T, require_std_vector_t<T>* = nullptr>
-inline auto corr_matrix_constrain(const T& y, int K) {
-  return apply_vector_unary<T>::apply(
-      y, [K](auto&& v) { return corr_matrix_constrain(v, K); });
+inline auto corr_matrix_constrain(T&& y, int K) {
+  return apply_vector_unary<std::decay_t<T>>::apply(
+      std::forward<T>(y), [K](auto&& v) {
+        return corr_matrix_constrain(std::forward<decltype(v)>(v), K);
+      });
 }
 
 /**
@@ -114,9 +116,11 @@ inline auto corr_matrix_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 corr_matrix_constrain(const T& y, int K, Lp& lp) {
-  return apply_vector_unary<T>::apply(
-      y, [&lp, K](auto&& v) { return corr_matrix_constrain(v, K, lp); });
+inline auto corr_matrix_constrain(T&& y, int K, Lp& lp) {
+  return apply_vector_unary<std::decay_t<T>>::apply(
+      std::forward<T>(y), [&lp, K](auto&& v) {
+        return corr_matrix_constrain(std::forward<decltype(v)>(v), K, lp);
+      });
 }
 
 /**
@@ -142,11 +146,11 @@ inline auto corr_matrix_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 corr_matrix_constrain(const T& x, Eigen::Index k, Lp& lp) {
+inline auto corr_matrix_constrain(T&& x, Eigen::Index k, Lp& lp) {
   if constexpr (Jacobian) {
-    return corr_matrix_constrain(x, k, lp);
+    return corr_matrix_constrain(std::forward<T>(x), k, lp);
   } else {
-    return corr_matrix_constrain(x, k);
+    return corr_matrix_constrain(std::forward<T>(x), k);
   }
 }