Merge remote-tracking branch 'origin/develop' into feature/threadsafe…

…-matrixcl

stan/math/prim/fun/promote_scalar.hpp

@@ -59,9 +59,11 @@ inline auto promote_scalar(UnPromotedType&& x) {
 
 // Forward decl for iterating over tuples used in std::vector<tuple>
 template <typename PromotionScalars, typename UnPromotedTypes,
-          require_all_tuple_t<PromotionScalars, UnPromotedTypes>* = nullptr>
+          require_all_tuple_t<PromotionScalars, UnPromotedTypes>* = nullptr,
+          require_not_same_t<PromotionScalars, UnPromotedTypes>* = nullptr>
 inline constexpr promote_scalar_t<PromotionScalars, UnPromotedTypes>
 promote_scalar(UnPromotedTypes&& x);
+
 /**
  * Promote the scalar type of an standard vector to the requested type.
  *
@@ -93,7 +95,8 @@ inline auto promote_scalar(UnPromotedType&& x) {
  * @param x input
  */
 template <typename PromotionScalars, typename UnPromotedTypes,
-          require_all_tuple_t<PromotionScalars, UnPromotedTypes>*>
+          require_all_tuple_t<PromotionScalars, UnPromotedTypes>*,
+          require_not_same_t<PromotionScalars, UnPromotedTypes>*>
 inline constexpr promote_scalar_t<PromotionScalars, UnPromotedTypes>
 promote_scalar(UnPromotedTypes&& x) {
   return index_apply<std::tuple_size<std::decay_t<UnPromotedTypes>>::value>(

stan/math/prim/fun/stan_print.hpp

@@ -52,6 +52,10 @@ void stan_print(std::ostream* o, const EigMat& x) {
   *o << ']';
 }
 
+// forward decl to allow the next two overloads to call each other
+template <typename T, require_tuple_t<T>* = nullptr>
+void stan_print(std::ostream* o, const T& x);
+
 template <typename T, require_std_vector_t<T>* = nullptr>
 void stan_print(std::ostream* o, const T& x) {
   *o << '[';
@@ -64,7 +68,7 @@ void stan_print(std::ostream* o, const T& x) {
   *o << ']';
 }
 
-template <typename T, require_tuple_t<T>* = nullptr>
+template <typename T, require_tuple_t<T>*>
 void stan_print(std::ostream* o, const T& x) {
   *o << '(';
   constexpr auto tuple_size = std::tuple_size<std::decay_t<T>>::value;

stan/math/prim/prob/bernoulli_logit_glm_lpmf.hpp

@@ -8,10 +8,11 @@
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/fun/isfinite.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 <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <cmath>
 
@@ -54,11 +55,16 @@ return_type_t<T_x, T_alpha, T_beta> bernoulli_logit_glm_lpmf(
   using Eigen::log1p;
   using Eigen::Matrix;
   using std::exp;
+  using std::isfinite;
   constexpr int T_x_rows = T_x::RowsAtCompileTime;
+  using T_xbeta_partials = partials_return_t<T_x, T_beta>;
   using T_partials_return = partials_return_t<T_y, T_x, T_alpha, T_beta>;
   using T_ytheta_tmp =
       typename std::conditional_t<T_x_rows == 1, T_partials_return,
                                   Array<T_partials_return, Dynamic, 1>>;
+  using T_xbeta_tmp =
+      typename std::conditional_t<T_x_rows == 1, T_xbeta_partials,
+                                  Array<T_xbeta_partials, Dynamic, 1>>;
   using T_x_ref = ref_type_if_t<!is_constant<T_x>::value, T_x>;
   using T_alpha_ref = ref_type_if_t<!is_constant<T_alpha>::value, T_alpha>;
   using T_beta_ref = ref_type_if_t<!is_constant<T_beta>::value, T_beta>;
@@ -86,11 +92,10 @@ return_type_t<T_x, T_alpha, T_beta> bernoulli_logit_glm_lpmf(
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
 
-  const auto& y_val = value_of_rec(y_ref);
-  const auto& x_val
-      = to_ref_if<!is_constant<T_beta>::value>(value_of_rec(x_ref));
-  const auto& alpha_val = value_of_rec(alpha_ref);
-  const auto& beta_val = value_of_rec(beta_ref);
+  const auto& y_val = value_of(y_ref);
+  const auto& x_val = to_ref_if<!is_constant<T_beta>::value>(value_of(x_ref));
+  const auto& alpha_val = value_of(alpha_ref);
+  const auto& beta_val = value_of(beta_ref);
 
   const auto& y_val_vec = as_column_vector_or_scalar(y_val);
   const auto& alpha_val_vec = as_column_vector_or_scalar(alpha_val);
@@ -103,7 +108,7 @@ return_type_t<T_x, T_alpha, T_beta> bernoulli_logit_glm_lpmf(
   Array<T_partials_return, Dynamic, 1> ytheta(N_instances);
   if (T_x_rows == 1) {
     T_ytheta_tmp ytheta_tmp
-        = forward_as<T_ytheta_tmp>((x_val * beta_val_vec)(0, 0));
+        = forward_as<T_xbeta_tmp>((x_val * beta_val_vec)(0, 0));
     ytheta = signs * (ytheta_tmp + as_array_or_scalar(alpha_val_vec));
   } else {
     ytheta = (x_val * beta_val_vec).array();
@@ -120,7 +125,7 @@ return_type_t<T_x, T_alpha, T_beta> bernoulli_logit_glm_lpmf(
           .select(-exp_m_ytheta,
                   (ytheta < -cutoff).select(ytheta, -log1p(exp_m_ytheta))));
 
-  if (!std::isfinite(logp)) {
+  if (!isfinite(logp)) {
     check_finite(function, "Weight vector", beta);
     check_finite(function, "Intercept", alpha);
     check_finite(function, "Matrix of independent variables", ytheta);
@@ -133,7 +138,7 @@ return_type_t<T_x, T_alpha, T_beta> bernoulli_logit_glm_lpmf(
         = (ytheta > cutoff)
               .select(-exp_m_ytheta,
                       (ytheta < -cutoff)
-                          .select(signs,
+                          .select(signs * T_partials_return(1.0),
                                   signs * exp_m_ytheta / (exp_m_ytheta + 1)));
     if (!is_constant_all<T_beta>::value) {
       if (T_x_rows == 1) {

stan/math/prim/prob/categorical_logit_glm_lpmf.hpp

@@ -6,12 +6,13 @@
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/fun/isfinite.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.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 <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <cmath>
@@ -50,11 +51,13 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::exp;
+  using std::isfinite;
   using std::log;
   using T_y_ref = ref_type_t<T_y>;
   using T_x_ref = ref_type_if_t<!is_constant<T_x>::value, T_x>;
   using T_alpha_ref = ref_type_if_t<!is_constant<T_alpha>::value, T_alpha>;
   using T_beta_ref = ref_type_if_t<!is_constant<T_beta>::value, T_beta>;
+  using T_beta_partials = partials_type_t<scalar_type_t<T_beta>>;
   constexpr int T_x_rows = T_x::RowsAtCompileTime;
 
   const size_t N_instances = T_x_rows == 1 ? stan::math::size(y) : x.rows();
@@ -82,11 +85,10 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
 
-  const auto& x_val
-      = to_ref_if<!is_constant<T_beta>::value>(value_of_rec(x_ref));
-  const auto& alpha_val = value_of_rec(alpha_ref);
+  const auto& x_val = to_ref_if<!is_constant<T_beta>::value>(value_of(x_ref));
+  const auto& alpha_val = value_of(alpha_ref);
   const auto& beta_val
-      = to_ref_if<!is_constant<T_x>::value>(value_of_rec(beta_ref));
+      = to_ref_if<!is_constant<T_x>::value>(value_of(beta_ref));
 
   const auto& alpha_val_vec = as_column_vector_or_scalar(alpha_val).transpose();
 
@@ -117,7 +119,7 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
   // when we have newer Eigen  T_partials_return logp =
   // lin(Eigen::all,y-1).sum() + log(inv_sum_exp_lin).sum() - lin_max.sum();
 
-  if (!std::isfinite(logp)) {
+  if (!isfinite(logp)) {
     check_finite(function, "Weight vector", beta_ref);
     check_finite(function, "Intercept", alpha_ref);
     check_finite(function, "Matrix of independent variables", x_ref);
@@ -128,7 +130,7 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
 
   if (!is_constant_all<T_x>::value) {
     if (T_x_rows == 1) {
-      Array<double, 1, Dynamic> beta_y = beta_val.col(y_seq[0] - 1);
+      Array<T_beta_partials, 1, Dynamic> beta_y = beta_val.col(y_seq[0] - 1);
       for (int i = 1; i < N_instances; i++) {
         beta_y += beta_val.col(y_seq[i] - 1).array();
       }
@@ -137,7 +139,8 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
             - (exp_lin.matrix() * beta_val.transpose()).array().colwise()
                   * inv_sum_exp_lin * N_instances;
     } else {
-      Array<double, Dynamic, Dynamic> beta_y(N_instances, N_attributes);
+      Array<T_beta_partials, Dynamic, Dynamic> beta_y(N_instances,
+                                                      N_attributes);
       for (int i = 0; i < N_instances; i++) {
         beta_y.row(i) = beta_val.col(y_seq[i] - 1);
       }
@@ -166,12 +169,11 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
       }
     }
     if (!is_constant_all<T_beta>::value) {
-      Matrix<T_partials_return, Dynamic, Dynamic> beta_derivative;
+      Matrix<T_partials_return, Dynamic, Dynamic> beta_derivative
+          = x_val.transpose().template cast<T_partials_return>()
+            * neg_softmax_lin.matrix();
       if (T_x_rows == 1) {
-        beta_derivative
-            = x_val.transpose() * neg_softmax_lin.matrix() * N_instances;
-      } else {
-        beta_derivative = x_val.transpose() * neg_softmax_lin.matrix();
+        beta_derivative *= N_instances;
       }
 
       for (int i = 0; i < N_instances; i++) {

stan/math/prim/prob/neg_binomial_2_log_glm_lpmf.hpp

@@ -15,7 +15,7 @@
 #include <stan/math/prim/fun/size.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/fun/value_of.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <vector>
 #include <cmath>
@@ -72,6 +72,7 @@ return_type_t<T_x, T_alpha, T_beta, T_precision> neg_binomial_2_log_glm_lpmf(
   using Eigen::log1p;
   using Eigen::Matrix;
   constexpr int T_x_rows = T_x::RowsAtCompileTime;
+  using T_xbeta_partials = partials_return_t<T_x, T_beta>;
   using T_partials_return
       = partials_return_t<T_y, T_x, T_alpha, T_beta, T_precision>;
   using T_precision_val = typename std::conditional_t<
@@ -85,6 +86,9 @@ return_type_t<T_x, T_alpha, T_beta, T_precision> neg_binomial_2_log_glm_lpmf(
   using T_theta_tmp =
       typename std::conditional_t<T_x_rows == 1, T_partials_return,
                                   Array<T_partials_return, Dynamic, 1>>;
+  using T_xbeta_tmp =
+      typename std::conditional_t<T_x_rows == 1, T_xbeta_partials,
+                                  Array<T_xbeta_partials, Dynamic, 1>>;
   using T_x_ref = ref_type_if_t<!is_constant<T_x>::value, T_x>;
   using T_alpha_ref = ref_type_if_t<!is_constant<T_alpha>::value, T_alpha>;
   using T_beta_ref = ref_type_if_t<!is_constant<T_beta>::value, T_beta>;
@@ -103,8 +107,8 @@ return_type_t<T_x, T_alpha, T_beta, T_precision> neg_binomial_2_log_glm_lpmf(
   check_consistent_size(function, "Vector of intercepts", alpha, N_instances);
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
-  const auto& beta_val = value_of_rec(beta_ref);
-  const auto& alpha_val = value_of_rec(alpha_ref);
+  const auto& beta_val = value_of(beta_ref);
+  const auto& alpha_val = value_of(alpha_ref);
   const auto& beta_val_vec = to_ref(as_column_vector_or_scalar(beta_val));
   const auto& alpha_val_vec = to_ref(as_column_vector_or_scalar(alpha_val));
   check_finite(function, "Weight vector", beta_val_vec);
@@ -117,8 +121,8 @@ return_type_t<T_x, T_alpha, T_beta, T_precision> neg_binomial_2_log_glm_lpmf(
   const auto& y_ref = to_ref(y);
   T_phi_ref phi_ref = phi;
 
-  const auto& y_val = value_of_rec(y_ref);
-  const auto& phi_val = value_of_rec(phi_ref);
+  const auto& y_val = value_of(y_ref);
+  const auto& phi_val = value_of(phi_ref);
 
   const auto& y_val_vec = to_ref(as_column_vector_or_scalar(y_val));
   const auto& phi_val_vec = to_ref(as_column_vector_or_scalar(phi_val));
@@ -131,16 +135,15 @@ return_type_t<T_x, T_alpha, T_beta, T_precision> neg_binomial_2_log_glm_lpmf(
 
   T_x_ref x_ref = x;
 
-  const auto& x_val
-      = to_ref_if<!is_constant<T_beta>::value>(value_of_rec(x_ref));
+  const auto& x_val = to_ref_if<!is_constant<T_beta>::value>(value_of(x_ref));
 
   const auto& y_arr = as_array_or_scalar(y_val_vec);
   const auto& phi_arr = as_array_or_scalar(phi_val_vec);
 
   Array<T_partials_return, Dynamic, 1> theta(N_instances);
   if (T_x_rows == 1) {
     T_theta_tmp theta_tmp
-        = forward_as<T_theta_tmp>((x_val * beta_val_vec)(0, 0));
+        = forward_as<T_xbeta_tmp>((x_val * beta_val_vec)(0, 0));
     theta = theta_tmp + as_array_or_scalar(alpha_val_vec);
   } else {
     theta = (x_val * beta_val_vec).array();
@@ -171,10 +174,11 @@ return_type_t<T_x, T_alpha, T_beta, T_precision> neg_binomial_2_log_glm_lpmf(
         logp += multiply_log(phi_vec[n], phi_vec[n]) - lgamma(phi_vec[n]);
       }
     } else {
+      using T_phi_scalar = scalar_type_t<decltype(phi_val_vec)>;
       logp += N_instances
-              * (multiply_log(forward_as<double>(phi_val),
-                              forward_as<double>(phi_val))
-                 - lgamma(forward_as<double>(phi_val)));
+              * (multiply_log(forward_as<T_phi_scalar>(phi_val),
+                              forward_as<T_phi_scalar>(phi_val))
+                 - lgamma(forward_as<T_phi_scalar>(phi_val)));
     }
   }
   logp -= sum(y_plus_phi * logsumexp_theta_logphi);

stan/math/prim/prob/normal_id_glm_lpdf.hpp

@@ -6,12 +6,13 @@
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/isfinite.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #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/fun/value_of.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <cmath>
 
@@ -59,6 +60,7 @@ return_type_t<T_y, T_x, T_alpha, T_beta, T_scale> normal_id_glm_lpdf(
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using Eigen::VectorXd;
+  using std::isfinite;
   constexpr int T_x_rows = T_x::RowsAtCompileTime;
   using T_partials_return
       = partials_return_t<T_y, T_x, T_alpha, T_beta, T_scale>;
@@ -86,7 +88,7 @@ return_type_t<T_y, T_x, T_alpha, T_beta, T_scale> normal_id_glm_lpdf(
                         N_instances);
   check_consistent_size(function, "Vector of intercepts", alpha, N_instances);
   T_sigma_ref sigma_ref = sigma;
-  const auto& sigma_val = value_of_rec(sigma_ref);
+  const auto& sigma_val = value_of(sigma_ref);
   const auto& sigma_val_vec = to_ref(as_column_vector_or_scalar(sigma_val));
   check_positive_finite(function, "Scale vector", sigma_val_vec);
 
@@ -102,11 +104,10 @@ return_type_t<T_y, T_x, T_alpha, T_beta, T_scale> normal_id_glm_lpdf(
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
 
-  const auto& y_val = value_of_rec(y_ref);
-  const auto& x_val
-      = to_ref_if<!is_constant<T_beta>::value>(value_of_rec(x_ref));
-  const auto& alpha_val = value_of_rec(alpha_ref);
-  const auto& beta_val = value_of_rec(beta_ref);
+  const auto& y_val = value_of(y_ref);
+  const auto& x_val = to_ref_if<!is_constant<T_beta>::value>(value_of(x_ref));
+  const auto& alpha_val = value_of(alpha_ref);
+  const auto& beta_val = value_of(beta_ref);
 
   const auto& y_val_vec = as_column_vector_or_scalar(y_val);
   const auto& alpha_val_vec = as_column_vector_or_scalar(alpha_val);
@@ -116,7 +117,7 @@ return_type_t<T_y, T_x, T_alpha, T_beta, T_scale> normal_id_glm_lpdf(
   T_scale_val inv_sigma = 1.0 / as_array_or_scalar(sigma_val_vec);
 
   // the most efficient way to calculate this depends on template parameters
-  double y_scaled_sq_sum;
+  T_partials_return y_scaled_sq_sum;
 
   Array<T_partials_return, Dynamic, 1> y_scaled(N_instances);
   if (T_x_rows == 1) {
@@ -178,7 +179,8 @@ return_type_t<T_y, T_x, T_alpha, T_beta, T_scale> normal_id_glm_lpdf(
       } else {
         y_scaled_sq_sum = sum(y_scaled * y_scaled);
         partials<4>(ops_partials)[0]
-            = (y_scaled_sq_sum - N_instances) * forward_as<double>(inv_sigma);
+            = (y_scaled_sq_sum - N_instances)
+              * forward_as<partials_return_t<T_sigma_ref>>(inv_sigma);
       }
     } else {
       y_scaled_sq_sum = sum(y_scaled * y_scaled);
@@ -187,7 +189,7 @@ return_type_t<T_y, T_x, T_alpha, T_beta, T_scale> normal_id_glm_lpdf(
     y_scaled_sq_sum = sum(y_scaled * y_scaled);
   }
 
-  if (!std::isfinite(y_scaled_sq_sum)) {
+  if (!isfinite(y_scaled_sq_sum)) {
     check_finite(function, "Vector of dependent variables", y_val_vec);
     check_finite(function, "Weight vector", beta_val_vec);
     check_finite(function, "Intercept", alpha_val_vec);
@@ -204,7 +206,8 @@ return_type_t<T_y, T_x, T_alpha, T_beta, T_scale> normal_id_glm_lpdf(
     if (is_vector<T_scale>::value) {
       logp -= sum(log(sigma_val_vec));
     } else {
-      logp -= N_instances * log(forward_as<double>(sigma_val_vec));
+      logp -= N_instances
+              * log(forward_as<partials_return_t<T_sigma_ref>>(sigma_val_vec));
     }
   }
   logp -= 0.5 * y_scaled_sq_sum;