Merge branch 'develop' into sum-to-zero-matrix

Author: WardBrian

lib/tbb_2020.3/STAN_CHANGES

@@ -10,6 +10,7 @@ This file documents changes done for the stan-math project
 - Support for Windows ARM64 with RTools:
   - build/Makefile.tbb
     - L94 Wrapped the use of `--version-script` export in conditional on non-WINARM64
-  - build/windows.gcc.ino
+  - build/windows.gcc.inc
     - L84 Wrapped the use of `-flifetime-dse` flag in conditional on non-WINARM64
+    - L101 Wrapped the use of `-msse` in conditional on non-WINARM64
 

lib/tbb_2020.3/build/windows.gcc.inc

@@ -98,7 +98,11 @@ CPLUS_FLAGS += -DUSE_WINTHREAD
 CPLUS_FLAGS += -D_WIN32_WINNT=$(_WIN32_WINNT)
 
 # MinGW specific
-CPLUS_FLAGS += -DMINGW_HAS_SECURE_API=1 -D__MSVCRT_VERSION__=0x0700 -msse -mthreads
+CPLUS_FLAGS += -DMINGW_HAS_SECURE_API=1 -D__MSVCRT_VERSION__=0x0700 -mthreads
+# Unused under ARM64 and error with LLVM19+
+ifeq (, $(WINARM64))
+    CPLUS_FLAGS +=  -msse 
+endif
 
 CONLY = gcc
 debugger = gdb

stan/math/fwd/fun/hypergeometric_1F0.hpp

@@ -31,10 +31,10 @@ namespace math {
 template <typename Ta, typename Tz, typename FvarT = return_type_t<Ta, Tz>,
           require_all_stan_scalar_t<Ta, Tz>* = nullptr,
           require_any_fvar_t<Ta, Tz>* = nullptr>
-FvarT hypergeometric_1f0(const Ta& a, const Tz& z) {
+FvarT hypergeometric_1F0(const Ta& a, const Tz& z) {
   partials_type_t<Ta> a_val = value_of(a);
   partials_type_t<Tz> z_val = value_of(z);
-  FvarT rtn = FvarT(hypergeometric_1f0(a_val, z_val), 0.0);
+  FvarT rtn = FvarT(hypergeometric_1F0(a_val, z_val), 0.0);
   if (!is_constant_all<Ta>::value) {
     rtn.d_ += forward_as<FvarT>(a).d() * -rtn.val() * log1m(z_val);
   }

stan/math/fwd/fun/hypergeometric_2F1.hpp

@@ -30,11 +30,11 @@ namespace math {
 template <typename Ta1, typename Ta2, typename Tb, typename Tz,
           require_all_stan_scalar_t<Ta1, Ta2, Tb, Tz>* = nullptr,
           require_any_fvar_t<Ta1, Ta2, Tb, Tz>* = nullptr>
-inline return_type_t<Ta1, Ta1, Tb, Tz> hypergeometric_2F1(const Ta1& a1,
+inline return_type_t<Ta1, Ta2, Tb, Tz> hypergeometric_2F1(const Ta1& a1,
                                                           const Ta2& a2,
                                                           const Tb& b,
                                                           const Tz& z) {
-  using fvar_t = return_type_t<Ta1, Ta1, Tb, Tz>;
+  using fvar_t = return_type_t<Ta1, Ta2, Tb, Tz>;
 
   auto a1_val = value_of(a1);
   auto a2_val = value_of(a2);

stan/math/fwd/fun/hypergeometric_pFq.hpp

@@ -7,6 +7,8 @@
 #include <stan/math/prim/fun/dot_product.hpp>
 #include <stan/math/prim/fun/grad_pFq.hpp>
 #include <stan/math/prim/fun/hypergeometric_pFq.hpp>
+#include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 
 namespace stan {
 namespace math {
@@ -30,33 +32,27 @@ template <typename Ta, typename Tb, typename Tz,
           bool grad_z = !is_constant<Tz>::value,
           require_all_vector_t<Ta, Tb>* = nullptr,
           require_fvar_t<FvarT>* = nullptr>
-inline FvarT hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) {
-  using PartialsT = partials_type_t<FvarT>;
-  using ARefT = ref_type_t<Ta>;
-  using BRefT = ref_type_t<Tb>;
-
-  ARefT a_ref = a;
-  BRefT b_ref = b;
+inline FvarT hypergeometric_pFq(Ta&& a, Tb&& b, Tz&& z) {
+  auto&& a_ref = to_ref(as_column_vector_or_scalar(a));
+  auto&& b_ref = to_ref(as_column_vector_or_scalar(b));
   auto&& a_val = value_of(a_ref);
   auto&& b_val = value_of(b_ref);
   auto&& z_val = value_of(z);
-  PartialsT pfq_val = hypergeometric_pFq(a_val, b_val, z_val);
+
+  partials_type_t<FvarT> pfq_val = hypergeometric_pFq(a_val, b_val, z_val);
   auto grad_tuple
       = grad_pFq<grad_a, grad_b, grad_z>(pfq_val, a_val, b_val, z_val);
 
   FvarT rtn = FvarT(pfq_val, 0.0);
 
-  if (grad_a) {
-    rtn.d_ += dot_product(forward_as<promote_scalar_t<FvarT, ARefT>>(a_ref).d(),
-                          std::get<0>(grad_tuple));
+  if constexpr (grad_a) {
+    rtn.d_ += dot_product(a_ref.d(), std::get<0>(grad_tuple));
   }
-  if (grad_b) {
-    rtn.d_ += dot_product(forward_as<promote_scalar_t<FvarT, BRefT>>(b_ref).d(),
-                          std::get<1>(grad_tuple));
+  if constexpr (grad_b) {
+    rtn.d_ += dot_product(b_ref.d(), std::get<1>(grad_tuple));
   }
-  if (grad_z) {
-    rtn.d_ += forward_as<promote_scalar_t<FvarT, Tz>>(z).d_
-              * std::get<2>(grad_tuple);
+  if constexpr (grad_z) {
+    rtn.d_ += z.d_ * std::get<2>(grad_tuple);
   }
 
   return rtn;

stan/math/prim/fun/hypergeometric_1F0.hpp

@@ -28,9 +28,8 @@ namespace math {
  * @return Hypergeometric 1F0 function
  */
 template <typename Ta, typename Tz, require_all_arithmetic_t<Ta, Tz>* = nullptr>
-return_type_t<Ta, Tz> hypergeometric_1f0(const Ta& a, const Tz& z) {
-  constexpr const char* function = "hypergeometric_1f0";
-  check_less("hypergeometric_1f0", "abs(z)", std::fabs(z), 1.0);
+return_type_t<Ta, Tz> hypergeometric_1F0(const Ta& a, const Tz& z) {
+  check_less("hypergeometric_1F0", "abs(z)", std::fabs(z), 1.0);
 
   return boost::math::hypergeometric_1F0(a, z, boost_policy_t<>());
 }

stan/math/prim/fun/hypergeometric_2F1.hpp

@@ -43,7 +43,7 @@ namespace internal {
  * @return Gauss hypergeometric function
  */
 template <typename Ta1, typename Ta2, typename Tb, typename Tz,
-          typename RtnT = boost::optional<return_type_t<Ta1, Ta1, Tb, Tz>>,
+          typename RtnT = boost::optional<return_type_t<Ta1, Ta2, Tb, Tz>>,
           require_all_arithmetic_t<Ta1, Ta2, Tb, Tz>* = nullptr>
 inline RtnT hyper_2F1_special_cases(const Ta1& a1, const Ta2& a2, const Tb& b,
                                     const Tz& z) {
@@ -148,10 +148,10 @@ inline RtnT hyper_2F1_special_cases(const Ta1& a1, const Ta2& a2, const Tb& b,
  * @return Gauss hypergeometric function
  */
 template <typename Ta1, typename Ta2, typename Tb, typename Tz,
-          typename ScalarT = return_type_t<Ta1, Ta1, Tb, Tz>,
+          typename ScalarT = return_type_t<Ta1, Ta2, Tb, Tz>,
           typename OptT = boost::optional<ScalarT>,
           require_all_arithmetic_t<Ta1, Ta2, Tb, Tz>* = nullptr>
-inline return_type_t<Ta1, Ta1, Tb, Tz> hypergeometric_2F1(const Ta1& a1,
+inline return_type_t<Ta1, Ta2, Tb, Tz> hypergeometric_2F1(const Ta1& a1,
                                                           const Ta2& a2,
                                                           const Tb& b,
                                                           const Tz& z) {

stan/math/prim/fun/hypergeometric_3F2.hpp

@@ -114,13 +114,20 @@ template <typename Ta, typename Tb, typename Tz,
           require_all_vector_t<Ta, Tb>* = nullptr,
           require_stan_scalar_t<Tz>* = nullptr>
 inline auto hypergeometric_3F2(const Ta& a, const Tb& b, const Tz& z) {
-  check_3F2_converges("hypergeometric_3F2", a[0], a[1], a[2], b[0], b[1], z);
+  check_size_match("hypergeometric_3F2", "a", a.size(), "3", 3);
+  check_size_match("hypergeometric_3F2", "b", b.size(), "2", 2);
+
+  auto a_ref = to_vector(a);
+  auto b_ref = to_vector(b);
+
+  check_3F2_converges("hypergeometric_3F2", a_ref[0], a_ref[1], a_ref[2],
+                      b_ref[0], b_ref[1], z);
   // Boost's pFq throws convergence errors in some cases, fallback to naive
   // infinite-sum approach (tests pass for these)
-  if (z == 1.0 && (sum(b) - sum(a)) < 0.0) {
-    return internal::hypergeometric_3F2_infsum(a, b, z);
+  if (z == 1.0 && (sum(b_ref) - sum(a_ref)) < 0.0) {
+    return internal::hypergeometric_3F2_infsum(a_ref, b_ref, z);
   }
-  return hypergeometric_pFq(to_vector(a), to_vector(b), z);
+  return hypergeometric_pFq(a_ref, b_ref, z);
 }
 
 /**

stan/math/prim/fun/hypergeometric_pFq.hpp

@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err/check_not_nan.hpp>
 #include <stan/math/prim/err/check_finite.hpp>
+#include <stan/math/prim/fun/to_row_vector.hpp>
 #include <boost/math/special_functions/hypergeometric_pFq.hpp>
 
 namespace stan {
@@ -14,10 +15,6 @@ namespace math {
  * input arguments:
  * \f$_pF_q(a_1,...,a_p;b_1,...,b_q;z)\f$
  *
- * This function is not intended to be exposed to end users, only
- * used for p & q values that are stable with the grad_pFq
- * implementation.
- *
  * See 'grad_pFq.hpp' for the derivatives wrt each parameter
  *
  * @param[in] a Vector of 'a' arguments to function
@@ -26,7 +23,7 @@ namespace math {
  * @return Generalized hypergeometric function
  */
 template <typename Ta, typename Tb, typename Tz,
-          require_all_eigen_st<std::is_arithmetic, Ta, Tb>* = nullptr,
+          require_all_vector_st<std::is_arithmetic, Ta, Tb>* = nullptr,
           require_arithmetic_t<Tz>* = nullptr>
 return_type_t<Ta, Tb, Tz> hypergeometric_pFq(const Ta& a, const Tb& b,
                                              const Tz& z) {
@@ -47,8 +44,9 @@ return_type_t<Ta, Tb, Tz> hypergeometric_pFq(const Ta& a, const Tb& b,
     std::stringstream msg;
     msg << "hypergeometric function pFq does not meet convergence "
         << "conditions with given arguments. "
-        << "a: " << a_ref << ", b: " << b_ref << ", "
-        << ", z: " << z;
+        << "a: " << to_row_vector(a_ref) << ", "
+        << "b: " << to_row_vector(b_ref) << ", "
+        << "z: " << z;
     throw std::domain_error(msg.str());
   }
 

stan/math/prim/fun/inv_logit.hpp

@@ -49,15 +49,14 @@ namespace math {
  * @return Inverse logit of argument.
  */
 inline double inv_logit(double a) {
-  using std::exp;
   if (a < 0) {
-    double exp_a = exp(a);
+    double exp_a = std::exp(a);
     if (a < LOG_EPSILON) {
       return exp_a;
     }
-    return exp_a / (1 + exp_a);
+    return exp_a / (1.0 + exp_a);
   }
-  return inv(1 + exp(-a));
+  return inv(1 + std::exp(-a));
 }
 
 /**
@@ -69,28 +68,46 @@ inline double inv_logit(double a) {
  */
 struct inv_logit_fun {
   template <typename T>
-  static inline auto fun(const T& x) {
-    return inv_logit(x);
+  static inline auto fun(T&& x) {
+    return inv_logit(std::forward<T>(x));
   }
 };
 
 /**
- * Vectorized version of inv_logit().
+ * Vectorized version of inv_logit() for containers containing ad types.
  *
- * @tparam T type of container
- * @param x container
+ * @tparam T type of std::vector
+ * @param x std::vector
  * @return Inverse logit applied to each value in x.
  */
-template <
-    typename T, require_not_var_matrix_t<T>* = nullptr,
-    require_all_not_nonscalar_prim_or_rev_kernel_expression_t<T>* = nullptr>
-inline auto inv_logit(const T& x) {
-  return apply_scalar_unary<inv_logit_fun, T>::apply(x);
+template <typename Container, require_ad_container_t<Container>* = nullptr,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              Container>* = nullptr,
+          require_not_rev_matrix_t<Container>* = nullptr>
+inline auto inv_logit(Container&& x) {
+  return apply_scalar_unary<inv_logit_fun, Container>::apply(
+      std::forward<Container>(x));
 }
 
-// TODO(Tadej): Eigen is introducing their implementation logistic() of this
-// in 3.4. Use that once we switch to Eigen 3.4
-
+/**
+ * Vectorized version of inv_logit() for containers with arithmetic scalar
+ * types.
+ *
+ * @tparam T A type of either `std::vector` or a type that directly inherits
+ * from `Eigen::DenseBase`. The inner scalar type must not have an auto diff
+ * scalar type.
+ * @param x Eigen expression
+ * @return Inverse logit applied to each value in x.
+ */
+template <typename Container,
+          require_container_bt<std::is_arithmetic, Container>* = nullptr,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              Container>* = nullptr>
+inline auto inv_logit(Container&& x) {
+  return apply_vector_unary<Container>::apply(
+      std::forward<Container>(x),
+      [](const auto& v) { return v.array().logistic(); });
+}
 }  // namespace math
 }  // namespace stan