Merge branch 'develop' into feature/scalar_binary_complex

Author: Andrew Johnson

stan/math/fwd/core/operator_division.hpp

@@ -31,7 +31,7 @@ inline fvar<T> operator/(const fvar<T>& x1, const fvar<T>& x2) {
  * @param x2 second argument
  * @return first argument divided by second argument
  */
-template <typename T, typename U, require_arithmetic_t<U>...>
+template <typename T, typename U, require_arithmetic_t<U>* = nullptr>
 inline fvar<T> operator/(const fvar<T>& x1, U x2) {
   return fvar<T>(x1.val_ / x2, x1.d_ / x2);
 }
@@ -44,7 +44,7 @@ inline fvar<T> operator/(const fvar<T>& x1, U x2) {
  * @param x2 second argument
  * @return first argument divided by second argument
  */
-template <typename T, typename U, require_arithmetic_t<U>...>
+template <typename T, typename U, require_arithmetic_t<U>* = nullptr>
 inline fvar<T> operator/(U x1, const fvar<T>& x2) {
   return fvar<T>(x1 / x2.val_, -x1 * x2.d_ / (x2.val_ * x2.val_));
 }
@@ -54,7 +54,7 @@ inline std::complex<fvar<T>> operator/(const std::complex<fvar<T>>& x1,
                                        const std::complex<fvar<T>>& x2) {
   return internal::complex_divide(x1, x2);
 }
-template <typename T, typename U, require_arithmetic_t<U>...>
+template <typename T, typename U, require_arithmetic_t<U>* = nullptr>
 inline std::complex<fvar<T>> operator/(const std::complex<fvar<T>>& x1,
                                        const std::complex<U>& x2) {
   return internal::complex_divide(x1, x2);
@@ -64,17 +64,17 @@ inline std::complex<fvar<T>> operator/(const std::complex<fvar<T>>& x1,
                                        const fvar<T>& x2) {
   return internal::complex_divide(x1, x2);
 }
-template <typename T, typename U, require_arithmetic_t<U>...>
+template <typename T, typename U, require_arithmetic_t<U>* = nullptr>
 inline std::complex<fvar<T>> operator/(const std::complex<fvar<T>>& x1, U x2) {
   return internal::complex_divide(x1, x2);
 }
 
-template <typename T, typename U, require_arithmetic_t<U>...>
+template <typename T, typename U, require_arithmetic_t<U>* = nullptr>
 inline std::complex<fvar<T>> operator/(const std::complex<U>& x1,
                                        const std::complex<fvar<T>>& x2) {
   return internal::complex_divide(x1, x2);
 }
-template <typename T, typename U, require_arithmetic_t<U>...>
+template <typename T, typename U, require_arithmetic_t<U>* = nullptr>
 inline std::complex<fvar<T>> operator/(const std::complex<U>& x1,
                                        const fvar<T>& x2) {
   return internal::complex_divide(x1, x2);
@@ -92,7 +92,7 @@ inline std::complex<fvar<T>> operator/(const fvar<T>& x1,
   return internal::complex_divide(x1, x2);
 }
 
-template <typename T, typename U, require_arithmetic_t<U>...>
+template <typename T, typename U, require_arithmetic_t<U>* = nullptr>
 inline std::complex<fvar<T>> operator/(U x1, const std::complex<fvar<T>>& x2) {
   return internal::complex_divide(x1, x2);
 }

stan/math/fwd/fun/is_nan.hpp

@@ -18,7 +18,7 @@ namespace math {
  * @param x Value to test.
  * @return <code>1</code> if the value is NaN and <code>0</code> otherwise.
  */
-template <typename T, require_fvar_t<T>...>
+template <typename T, require_fvar_t<T>* = nullptr>
 inline bool is_nan(T&& x) {
   return is_nan(std::forward<decltype(x.val())>(x.val()));
 }

stan/math/fwd/fun/log_mix.hpp

@@ -114,7 +114,7 @@ inline fvar<T> log_mix(const fvar<T>& theta, const fvar<T>& lambda1,
   }
 }
 
-template <typename T, typename P, require_all_arithmetic_t<P>...>
+template <typename T, typename P, require_all_arithmetic_t<P>* = nullptr>
 inline fvar<T> log_mix(const fvar<T>& theta, const fvar<T>& lambda1,
                        P lambda2) {
   if (lambda1.val_ > lambda2) {
@@ -132,7 +132,7 @@ inline fvar<T> log_mix(const fvar<T>& theta, const fvar<T>& lambda1,
   }
 }
 
-template <typename T, typename P, require_all_arithmetic_t<P>...>
+template <typename T, typename P, require_all_arithmetic_t<P>* = nullptr>
 inline fvar<T> log_mix(const fvar<T>& theta, P lambda1,
                        const fvar<T>& lambda2) {
   if (lambda1 > lambda2.val_) {
@@ -150,7 +150,7 @@ inline fvar<T> log_mix(const fvar<T>& theta, P lambda1,
   }
 }
 
-template <typename T, typename P, require_all_arithmetic_t<P>...>
+template <typename T, typename P, require_all_arithmetic_t<P>* = nullptr>
 inline fvar<T> log_mix(P theta, const fvar<T>& lambda1,
                        const fvar<T>& lambda2) {
   if (lambda1.val_ > lambda2.val_) {
@@ -169,7 +169,7 @@ inline fvar<T> log_mix(P theta, const fvar<T>& lambda1,
 }
 
 template <typename T, typename P1, typename P2,
-          require_all_arithmetic_t<P1, P2>...>
+          require_all_arithmetic_t<P1, P2>* = nullptr>
 inline fvar<T> log_mix(const fvar<T>& theta, P1 lambda1, P2 lambda2) {
   if (lambda1 > lambda2) {
     fvar<T> partial_deriv_array[1];
@@ -185,7 +185,7 @@ inline fvar<T> log_mix(const fvar<T>& theta, P1 lambda1, P2 lambda2) {
 }
 
 template <typename T, typename P1, typename P2,
-          require_all_arithmetic_t<P1, P2>...>
+          require_all_arithmetic_t<P1, P2>* = nullptr>
 inline fvar<T> log_mix(P1 theta, const fvar<T>& lambda1, P2 lambda2) {
   if (lambda1.val_ > lambda2) {
     fvar<T> partial_deriv_array[1];
@@ -201,7 +201,7 @@ inline fvar<T> log_mix(P1 theta, const fvar<T>& lambda1, P2 lambda2) {
 }
 
 template <typename T, typename P1, typename P2,
-          require_all_arithmetic_t<P1, P2>...>
+          require_all_arithmetic_t<P1, P2>* = nullptr>
 inline fvar<T> log_mix(P1 theta, P2 lambda1, const fvar<T>& lambda2) {
   if (lambda1 > lambda2.val_) {
     fvar<T> partial_deriv_array[1];

stan/math/fwd/fun/log_sum_exp.hpp

@@ -50,7 +50,7 @@ inline fvar<T> log_sum_exp(const fvar<T>& x1, double x2) {
  * @param[in] x Matrix of specified values.
  * @return The log of the sum of the exponentiated vector values.
  */
-template <typename T, require_container_st<is_fvar, T>...>
+template <typename T, require_container_st<is_fvar, T>* = nullptr>
 inline auto log_sum_exp(const T& x) {
   return apply_vector_unary<T>::reduce(x, [&](const auto& v) {
     using T_fvar_inner = typename value_type_t<decltype(v)>::Scalar;

stan/math/opencl/copy.hpp

@@ -41,7 +41,7 @@ namespace math {
  * @return matrix_cl with a copy of the data in the source matrix
  */
 template <typename Mat, typename Mat_scalar = scalar_type_t<Mat>,
-          require_eigen_vt<std::is_arithmetic, Mat>...>
+          require_eigen_vt<std::is_arithmetic, Mat>* = nullptr>
 inline matrix_cl<Mat_scalar> to_matrix_cl(Mat&& src) {
   return matrix_cl<Mat_scalar>(std::forward<Mat>(src));
 }
@@ -60,7 +60,7 @@ inline matrix_cl<Mat_scalar> to_matrix_cl(Mat&& src) {
  * @return matrix_cl with a copy of the data in the source matrix
  */
 template <typename Vec, typename Vec_scalar = scalar_type_t<Vec>,
-          require_std_vector_vt<std::is_arithmetic, Vec>...>
+          require_std_vector_vt<std::is_arithmetic, Vec>* = nullptr>
 inline matrix_cl<Vec_scalar> to_matrix_cl(Vec&& src) {
   return matrix_cl<Vec_scalar>(std::forward<Vec>(src));
 }
@@ -158,7 +158,7 @@ inline std::vector<T> packed_copy(const matrix_cl<T>& src) {
  */
 template <matrix_cl_view matrix_view, typename Vec,
           typename Vec_scalar = scalar_type_t<Vec>,
-          require_vector_vt<std::is_arithmetic, Vec>...>
+          require_vector_vt<std::is_arithmetic, Vec>* = nullptr>
 inline matrix_cl<Vec_scalar> packed_copy(Vec&& src, int rows) {
   const int packed_size = rows * (rows + 1) / 2;
   check_size_match("copy (packed std::vector -> OpenCL)", "src.size()",

stan/math/opencl/kernel_cl.hpp

@@ -112,7 +112,7 @@ inline void assign_event(const cl::Event& e,
   helper.set(e, m);
 }
 
-template <typename T, require_same_t<T, cl::Event>...>
+template <typename T, require_same_t<T, cl::Event>* = nullptr>
 inline void assign_events(const T&) {}
 
 /** \ingroup kernel_executor_opencl

stan/math/opencl/kernel_generator/holder_cl.hpp

@@ -74,35 +74,6 @@ auto holder_cl(T&& a, Ptrs*... ptrs) {
 }
 
 namespace internal {
-/**
- * Handles single element (moving rvalues to heap) for construction of
- * `holder_cl` from a functor. For lvalues just sets the `res` pointer.
- * @tparam T type of the element
- * @param a element to handle
- * @param res resulting pointer to element
- * @return tuple of pointer allocated on heap (empty).
- */
-template <typename T>
-auto holder_cl_handle_element(const T& a, const T*& res) {
-  res = &a;
-  return std::make_tuple();
-}
-
-/**
- * Handles single element (moving rvalues to heap) for construction of
- * `holder_cl` from a functor. Rvalue is moved to heap and the pointer to heap
- * memory is assigned to res and returned in a tuple.
- * @tparam T type of the element
- * @param a element to handle
- * @param res resulting pointer to element
- * @return tuple of pointer allocated on heap (containing single pointer).
- */
-template <typename T>
-auto holder_cl_handle_element(std::remove_reference_t<T>&& a, const T*& res) {
-  res = new T(std::move(a));
-  return std::make_tuple(res);
-}
-
 /**
  * Second step in implementation of construction `holder_cl` from a functor.
  * @tparam T type of the result expression
@@ -114,8 +85,8 @@ auto holder_cl_handle_element(std::remove_reference_t<T>&& a, const T*& res) {
  * @return `holder_cl` referencing given expression
  */
 template <typename T, std::size_t... Is, typename... Args>
-auto make_holder_cl_impl2(T&& expr, std::index_sequence<Is...>,
-                          const std::tuple<Args*...>& ptrs) {
+auto make_holder_cl_impl_step2(T&& expr, std::index_sequence<Is...>,
+                               const std::tuple<Args*...>& ptrs) {
   return holder_cl(std::forward<T>(expr), std::get<Is>(ptrs)...);
 }
 
@@ -129,12 +100,12 @@ auto make_holder_cl_impl2(T&& expr, std::index_sequence<Is...>,
  * @return `holder_cl` referencing given expression
  */
 template <typename T, std::size_t... Is, typename... Args>
-auto make_holder_cl_impl(const T& func, std::index_sequence<Is...>,
-                         Args&&... args) {
-  std::tuple<const std::remove_reference_t<Args>*...> res;
+auto make_holder_cl_impl_step1(const T& func, std::index_sequence<Is...>,
+                               Args&&... args) {
+  std::tuple<std::remove_reference_t<Args>*...> res;
   auto ptrs = std::tuple_cat(
-      holder_cl_handle_element(std::forward<Args>(args), std::get<Is>(res))...);
-  return make_holder_cl_impl2(
+      holder_handle_element(std::forward<Args>(args), std::get<Is>(res))...);
+  return make_holder_cl_impl_step2(
       func(*std::get<Is>(res)...),
       std::make_index_sequence<std::tuple_size<decltype(ptrs)>::value>(), ptrs);
 }
@@ -156,7 +127,7 @@ template <typename T, typename... Args,
               decltype(std::declval<T>()(std::declval<Args&>()...)),
               Args...>* = nullptr>
 auto make_holder_cl(const T& func, Args&&... args) {
-  return internal::make_holder_cl_impl(
+  return internal::make_holder_cl_impl_step1(
       func, std::make_index_sequence<sizeof...(Args)>(),
       std::forward<Args>(args)...);
 }

stan/math/opencl/matrix_cl.hpp

@@ -233,8 +233,8 @@ class matrix_cl<T, require_arithmetic_t<T>> {
    * @throw <code>std::system_error</code> if the memory on the device could not
    * be allocated
    */
-  template <typename Vec, require_std_vector_vt<is_eigen, Vec>...,
-            require_st_same<Vec, T>...>
+  template <typename Vec, require_std_vector_vt<is_eigen, Vec>* = nullptr,
+            require_st_same<Vec, T>* = nullptr>
   explicit matrix_cl(Vec&& A) try : rows_(A.empty() ? 0 : A[0].size()),
                                     cols_(A.size()) {
     if (this->size() == 0) {
@@ -305,7 +305,8 @@ class matrix_cl<T, require_arithmetic_t<T>> {
    * @throw <code>std::system_error</code> if the memory on the device could not
    * be allocated
    */
-  template <typename Mat, require_eigen_t<Mat>..., require_vt_same<Mat, T>...>
+  template <typename Mat, require_eigen_t<Mat>* = nullptr,
+            require_vt_same<Mat, T>* = nullptr>
   explicit matrix_cl(Mat&& A,
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(A.rows()), cols_(A.cols()), view_(partial_view) {
@@ -360,8 +361,8 @@ class matrix_cl<T, require_arithmetic_t<T>> {
    * @throw <code>std::system_error</code> if the memory on the device could not
    * be allocated
    */
-  template <typename Vec, require_std_vector_t<Vec>...,
-            require_vt_same<Vec, T>...>
+  template <typename Vec, require_std_vector_t<Vec>* = nullptr,
+            require_vt_same<Vec, T>* = nullptr>
   explicit matrix_cl(Vec&& A,
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : matrix_cl(std::forward<Vec>(A), A.size(), 1) {}
@@ -384,8 +385,8 @@ class matrix_cl<T, require_arithmetic_t<T>> {
    * @throw <code>std::system_error</code> if the memory on the device could not
    * be allocated
    */
-  template <typename Vec, require_std_vector_t<Vec>...,
-            require_vt_same<Vec, T>...>
+  template <typename Vec, require_std_vector_t<Vec>* = nullptr,
+            require_vt_same<Vec, T>* = nullptr>
   explicit matrix_cl(Vec&& A, const int& R, const int& C,
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(R), cols_(C), view_(partial_view) {
@@ -409,7 +410,7 @@ class matrix_cl<T, require_arithmetic_t<T>> {
    * @throw <code>std::system_error</code> if the memory on the device could not
    * be allocated
    */
-  template <typename U, require_same_t<T, U>...>
+  template <typename U, require_same_t<T, U>* = nullptr>
   explicit matrix_cl(const U* A, const int& R, const int& C,
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(R), cols_(C), view_(partial_view) {

stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp

@@ -49,7 +49,8 @@ return_type_t<T_alpha, T_beta> bernoulli_logit_glm_lpmf(
     const T_alpha& alpha, const T_beta& beta) {
   static const char* function = "bernoulli_logit_glm_lpmf(OpenCL)";
   using T_partials_return = partials_return_t<T_alpha, T_beta>;
-
+  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 Eigen::Dynamic;
   using Eigen::Matrix;
 
@@ -71,15 +72,12 @@ return_type_t<T_alpha, T_beta> bernoulli_logit_glm_lpmf(
   if (N == 0) {
     return 0;
   }
-
   if (!include_summand<propto, T_alpha, T_beta>::value) {
     return 0;
   }
 
-  T_partials_return logp(0);
-
-  const auto& beta_ref = to_ref_if<!is_constant<T_beta>::value>(beta);
-  const auto& alpha_ref = to_ref_if<!is_constant<T_alpha>::value>(alpha);
+  T_beta_ref beta_ref = beta;
+  T_alpha_ref alpha_ref = alpha;
 
   const auto& beta_val = value_of_rec(beta_ref);
   const auto& alpha_val = value_of_rec(alpha_ref);
@@ -112,7 +110,6 @@ return_type_t<T_alpha, T_beta> bernoulli_logit_glm_lpmf(
                                  (exp_m_ytheta_expr + 1))));
 
   const int wgs = logp_expr.rows();
-
   matrix_cl<double> logp_cl(wgs, 1);
   constexpr bool need_theta_derivative
       = !is_constant_all<T_beta, T_alpha>::value;
@@ -126,8 +123,7 @@ return_type_t<T_alpha, T_beta> bernoulli_logit_glm_lpmf(
       logp_expr, calc_if<need_theta_derivative>(theta_derivative_expr),
       calc_if<need_theta_derivative_sum>(colwise_sum(theta_derivative_expr)));
 
-  logp += sum(from_matrix_cl<Eigen::Dynamic, 1>(logp_cl));
-
+  T_partials_return logp = sum(from_matrix_cl<Eigen::Dynamic, 1>(logp_cl));
   if (!std::isfinite(logp)) {
     check_bounded(function, "Vector of dependent variables",
                   from_matrix_cl(y_cl), 0, 1);
@@ -137,8 +133,8 @@ return_type_t<T_alpha, T_beta> bernoulli_logit_glm_lpmf(
                  from_matrix_cl(x_cl));
   }
 
-  operands_and_partials<decltype(alpha_ref), decltype(beta_ref)> ops_partials(
-      alpha_ref, beta_ref);
+  operands_and_partials<T_alpha_ref, T_beta_ref> ops_partials(alpha_ref,
+                                                              beta_ref);
   // Compute the necessary derivatives.
   if (!is_constant_all<T_alpha>::value) {
     if (is_alpha_vector) {

stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp

@@ -45,8 +45,8 @@ return_type_t<T_alpha, T_beta> categorical_logit_glm_lpmf(
     const matrix_cl<int>& y_cl, const matrix_cl<double>& x_cl,
     const T_alpha& alpha, const T_beta& beta) {
   using T_partials_return = partials_return_t<T_alpha, T_beta>;
-  static const char* function = "categorical_logit_glm_lpmf";
-
+  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 Eigen::Array;
   using Eigen::Dynamic;
   using Eigen::Matrix;
@@ -55,6 +55,7 @@ return_type_t<T_alpha, T_beta> categorical_logit_glm_lpmf(
   const size_t N_attributes = x_cl.cols();
   const size_t N_classes = beta.cols();
 
+  static const char* function = "categorical_logit_glm_lpmf";
   if (y_cl.size() != 1) {
     check_size_match(function, "x.rows()", N_instances, "y.size()",
                      y_cl.size());
@@ -66,13 +67,12 @@ return_type_t<T_alpha, T_beta> categorical_logit_glm_lpmf(
   if (N_instances == 0 || N_classes <= 1) {
     return 0;
   }
-
   if (!include_summand<propto, T_alpha, T_beta>::value) {
     return 0;
   }
 
-  const auto& alpha_ref = to_ref_if<!is_constant<T_alpha>::value>(alpha);
-  const auto& beta_ref = to_ref_if<!is_constant<T_beta>::value>(beta);
+  T_alpha_ref alpha_ref = alpha;
+  T_beta_ref beta_ref = beta;
 
   const auto& alpha_val = value_of_rec(alpha_ref);
   const auto& beta_val = value_of_rec(beta_ref);
@@ -118,8 +118,8 @@ return_type_t<T_alpha, T_beta> categorical_logit_glm_lpmf(
                  from_matrix_cl(x_cl));
   }
 
-  operands_and_partials<decltype(alpha_ref), decltype(beta_ref)> ops_partials(
-      alpha_ref, beta_ref);
+  operands_and_partials<T_alpha_ref, T_beta_ref> ops_partials(alpha_ref,
+                                                              beta_ref);
   if (!is_constant_all<T_alpha>::value) {
     ops_partials.edge1_.partials_
         = from_matrix_cl(alpha_derivative_cl).colwise().sum();