Let value_of and value_of_rec return expressions (#1872)

* let value_of and value_of_rec return expressions

* [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final)

* added missing include

* added missing include and generalized hmm_marginal_lpdf_val

* added more missing includes

* addressed review comments

* [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final)

Co-authored-by: Stan Jenkins <mc.stanislaw@gmail.com>

Author: t4c1

stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp

@@ -72,8 +72,6 @@ return_type_t<T_alpha_scalar, T_beta_scalar> categorical_logit_glm_lpmf(
   const auto& beta_val = value_of_rec(beta);
   const auto& alpha_val = value_of_rec(alpha);
 
-  const auto& alpha_val_vec = as_column_vector_or_scalar(alpha_val).transpose();
-
   const int local_size
       = opencl_kernels::categorical_logit_glm.get_option("LOCAL_SIZE_");
   const int wgs = (N_instances + local_size - 1) / local_size;

stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp

@@ -100,8 +100,6 @@ return_type_t<T_alpha, T_beta, T_precision> neg_binomial_2_log_glm_lpmf(
   const auto& alpha_val_vec = as_column_vector_or_scalar(alpha_val);
   const auto& phi_val_vec = as_column_vector_or_scalar(phi_val);
 
-  const auto& phi_arr = as_array_or_scalar(phi_val_vec);
-
   const int local_size
       = opencl_kernels::neg_binomial_2_log_glm.get_option("LOCAL_SIZE_");
   const int wgs = (N + local_size - 1) / local_size;

stan/math/opencl/prim/normal_id_glm_lpdf.hpp

@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/sum.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of_rec.hpp>
 #include <stan/math/prim/prob/normal_id_glm_lpdf.hpp>
 #include <stan/math/opencl/copy.hpp>
@@ -96,7 +97,7 @@ return_type_t<T_alpha, T_beta, T_scale> normal_id_glm_lpdf(
 
   const auto &beta_val_vec = as_column_vector_or_scalar(beta_val);
   const auto &alpha_val_vec = as_column_vector_or_scalar(alpha_val);
-  const auto &sigma_val_vec = as_column_vector_or_scalar(sigma_val);
+  const auto &sigma_val_vec = to_ref(as_column_vector_or_scalar(sigma_val));
 
   T_scale_val inv_sigma = 1 / as_array_or_scalar(sigma_val_vec);
   Matrix<T_partials_return, Dynamic, 1> y_minus_mu_over_sigma_mat(N);

stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp

@@ -83,9 +83,6 @@ return_type_t<T_beta_scalar, T_cuts_scalar> ordered_logistic_glm_lpmf(
   const auto& beta_val = value_of_rec(beta);
   const auto& cuts_val = value_of_rec(cuts);
 
-  const auto& beta_val_vec = as_column_vector_or_scalar(beta_val);
-  const auto& cuts_val_vec = as_column_vector_or_scalar(cuts_val);
-
   operands_and_partials<Eigen::Matrix<T_beta_scalar, Eigen::Dynamic, 1>,
                         Eigen::Matrix<T_cuts_scalar, Eigen::Dynamic, 1>>
       ops_partials(beta, cuts);

stan/math/prim/fun.hpp

@@ -306,6 +306,7 @@
 #include <stan/math/prim/fun/to_array_1d.hpp>
 #include <stan/math/prim/fun/to_array_2d.hpp>
 #include <stan/math/prim/fun/to_matrix.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/to_row_vector.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
 #include <stan/math/prim/fun/trace.hpp>

stan/math/prim/fun/log_mix.hpp

@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
 #include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <vector>
 #include <cmath>
@@ -86,8 +87,8 @@ return_type_t<T_theta, T_lam> log_mix(const T_theta& theta,
   check_finite(function, "theta", theta);
   check_consistent_sizes(function, "theta", theta, "lambda", lambda);
 
-  const auto& theta_dbl = value_of(as_column_vector_or_scalar(theta));
-  const auto& lam_dbl = value_of(as_column_vector_or_scalar(lambda));
+  const auto& theta_dbl = to_ref(value_of(as_column_vector_or_scalar(theta)));
+  const auto& lam_dbl = to_ref(value_of(as_column_vector_or_scalar(lambda)));
 
   T_partials_return logp = log_sum_exp(log(theta_dbl) + lam_dbl);
 
@@ -158,7 +159,7 @@ return_type_t<T_theta, std::vector<T_lam>> log_mix(
     check_consistent_sizes(function, "theta", theta, "lambda", lambda[n]);
   }
 
-  const auto& theta_dbl = value_of(as_column_vector_or_scalar(theta));
+  const auto& theta_dbl = to_ref(value_of(as_column_vector_or_scalar(theta)));
 
   T_partials_mat lam_dbl(M, N);
   for (int n = 0; n < N; ++n) {

stan/math/prim/fun/to_ref.hpp

@@ -0,0 +1,34 @@
+#ifndef STAN_MATH_PRIM_FUN_TO_REF_HPP
+#define STAN_MATH_PRIM_FUN_TO_REF_HPP
+
+#include <stan/math/prim/meta.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * No-op that should be optimized away.
+ * @tparam T non-Eigen argument type
+ * @param a argument
+ * @return argument
+ */
+template <typename T, require_not_eigen_t<T>* = nullptr>
+inline T to_ref(T&& a) {
+  return std::forward<T>(a);
+}
+
+/**
+ * Converts Eigen argument into `Eigen::Ref`. This evaluate expensive
+ * expressions.
+ * @tparam T argument type (Eigen expression)
+ * @param a argument
+ * @return argument converted to `Eigen::Ref`
+ */
+template <typename T, require_eigen_t<T>* = nullptr>
+inline Eigen::Ref<const plain_type_t<T>> to_ref(T&& a) {
+  return std::forward<T>(a);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/prim/fun/value_of.hpp

@@ -99,22 +99,15 @@ inline Vec value_of(Vec&& x) {
  * T must implement value_of. See
  * test/math/fwd/fun/value_of.cpp for fvar and var usage.
  *
- * @tparam T type of elements in the matrix
- * @tparam R number of rows in the matrix, can be Eigen::Dynamic
- * @tparam C number of columns in the matrix, can be Eigen::Dynamic
+ * @tparam EigMat type of the matrix
  *
  * @param[in] M Matrix to be converted
  * @return Matrix of values
  **/
 template <typename EigMat, require_eigen_t<EigMat>* = nullptr,
           require_not_vt_double_or_int<EigMat>* = nullptr>
-inline Eigen::Matrix<typename child_type<value_type_t<EigMat>>::type,
-                     EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime>
-value_of(const EigMat& M) {
-  return M.array()
-      .unaryExpr([](const auto& scal) { return value_of(scal); })
-      .matrix()
-      .eval();
+inline auto value_of(const EigMat& M) {
+  return M.unaryExpr([](const auto& scal) { return value_of(scal); });
 }
 
 /**
@@ -125,8 +118,7 @@ value_of(const EigMat& M) {
  *
  * <p>This inline pass-through no-op should be compiled away.
  *
- * @tparam R number of rows in the matrix, can be Eigen::Dynamic
- * @tparam C number of columns in the matrix, can be Eigen::Dynamic
+ * @tparam EigMat type of the matrix
  *
  * @param x Specified matrix.
  * @return Specified matrix.

stan/math/prim/fun/value_of_rec.hpp

@@ -65,7 +65,7 @@ inline std::complex<double> value_of_rec(const std::complex<T>& x) {
  * @param[in] x std::vector to be converted
  * @return std::vector of values
  **/
-template <typename T>
+template <typename T, require_not_same_t<double, T>* = nullptr>
 inline std::vector<double> value_of_rec(const std::vector<T>& x) {
   size_t x_size = x.size();
   std::vector<double> result(x_size);
@@ -86,8 +86,10 @@ inline std::vector<double> value_of_rec(const std::vector<T>& x) {
  * @param x Specified std::vector.
  * @return Specified std::vector.
  */
-inline const std::vector<double>& value_of_rec(const std::vector<double>& x) {
-  return x;
+template <typename T, require_std_vector_t<T>* = nullptr,
+          require_vt_same<double, T>* = nullptr>
+inline T value_of_rec(T&& x) {
+  return std::forward<T>(x);
 }
 
 /**
@@ -120,8 +122,8 @@ inline auto value_of_rec(const T& M) {
  */
 template <typename T, typename = require_st_same<T, double>,
           typename = require_eigen_t<T>>
-inline const T& value_of_rec(const T& x) {
-  return x;
+inline T value_of_rec(T&& x) {
+  return std::forward<T>(x);
 }
 }  // namespace math
 }  // namespace stan

stan/math/prim/prob/bernoulli_logit_glm_lpmf.hpp

@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of_rec.hpp>
 #include <cmath>
 
@@ -82,13 +83,13 @@ return_type_t<T_x_scalar, T_alpha, T_beta> bernoulli_logit_glm_lpmf(
   }
 
   T_partials_return logp(0);
-  const auto &x_val = value_of_rec(x);
+  const auto &x_val = to_ref(value_of_rec(x));
   const auto &y_val = value_of_rec(y);
   const auto &beta_val = value_of_rec(beta);
   const auto &alpha_val = value_of_rec(alpha);
 
   const auto &y_val_vec = as_column_vector_or_scalar(y_val);
-  const auto &beta_val_vec = as_column_vector_or_scalar(beta_val);
+  const auto &beta_val_vec = to_ref(as_column_vector_or_scalar(beta_val));
   const auto &alpha_val_vec = as_column_vector_or_scalar(alpha_val);
 
   T_y_val signs = 2 * as_array_or_scalar(y_val_vec) - 1;