Merge remote-tracking branch 'origin' into update/eigen_5.0.1

Author: SteveBronder

Jenkinsfile

@@ -487,6 +487,12 @@ pipeline {
             steps {
                 script {
                     retry(3) { checkout scm }
+
+                    if (params.withRowVector || isBranch('develop') || isBranch('master')) {
+                        sh "echo CXXFLAGS+=-DSTAN_TEST_ROW_VECTORS >> make/local"
+                        sh "echo CXXFLAGS+=-DSTAN_PROB_TEST_ALL >> make/local"
+                    }
+
                     if (params.runAllDistributions || isBranch('develop') || isBranch('master')) {
                         changedDistributionTests = sh(script:"python3 test/prob/getDependencies.py --pretend-all", returnStdout:true).trim().readLines()
                     } else {
@@ -511,12 +517,17 @@ pipeline {
             steps {
                 script {
                     def tests = [:]
-                    for (f in changedDistributionTests.collate(24)) {
+
+                    def executors = 10
+                    def tests_per_executor = Math.ceil((changedDistributionTests.size() / executors).doubleValue()).toInteger()
+                    def idx = 0
+                    for (f in changedDistributionTests.collate(tests_per_executor)) {
+                        idx = idx + 1
                         def names = f.join(" ")
-                        tests["Distribution Tests: ${names}"] = { node ("linux && docker && 8core") {
+                        tests["Distribution Tests: ${idx} / ${executors}"] = { node ("linux && docker && 8core") {
                             deleteDir()
                             docker.image('stanorg/ci:gpu-cpp17').inside {
-                                catchError {
+                                catchError(buildResult: "FAILURE", stageResult: "FAILURE") {
                                     unstash 'MathSetup'
                                     sh """
                                         echo CXX=${CLANG_CXX} > make/local

stan/math/fwd/meta/is_fvar.hpp

@@ -21,8 +21,5 @@ struct is_fvar<T,
                std::enable_if_t<internal::is_fvar_impl<std::decay_t<T>>::value>>
     : std::true_type {};
 
-template <typename T>
-inline constexpr bool is_fvar_v = is_fvar<T>::value;
-
 }  // namespace stan
 #endif

stan/math/prim/fun/log_gamma_q_dgamma.hpp

@@ -0,0 +1,134 @@
+#ifndef STAN_MATH_PRIM_FUN_LOG_GAMMA_Q_DGAMMA_HPP
+#define STAN_MATH_PRIM_FUN_LOG_GAMMA_Q_DGAMMA_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/digamma.hpp>
+#include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/fun/fabs.hpp>
+#include <stan/math/prim/fun/gamma_p.hpp>
+#include <stan/math/prim/fun/gamma_q.hpp>
+#include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
+#include <stan/math/prim/fun/lgamma.hpp>
+#include <stan/math/prim/fun/log.hpp>
+#include <stan/math/prim/fun/log1m.hpp>
+#include <stan/math/prim/fun/tgamma.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+
+constexpr double LOG_Q_GAMMA_CF_PRECISION = 1.49012e-12;
+
+/**
+ * Compute log(Q(a,z)) using continued fraction expansion for upper incomplete
+ * gamma function.
+ *
+ * @tparam T_a Type of shape parameter a (double or fvar types)
+ * @tparam T_z Type of value parameter z (double or fvar types)
+ * @param a Shape parameter
+ * @param z Value at which to evaluate
+ * @param precision Convergence threshold, default of sqrt(machine_epsilon)
+ * @param max_steps Maximum number of continued fraction iterations
+ * @return log(Q(a,z)) with the return type of T_a and T_z
+ */
+template <typename T_a, typename T_z>
+inline return_type_t<T_a, T_z> log_q_gamma_cf(const T_a& a, const T_z& z,
+                                              double precision
+                                              = LOG_Q_GAMMA_CF_PRECISION,
+                                              int max_steps = 250) {
+  using T_return = return_type_t<T_a, T_z>;
+  const T_return log_prefactor = a * log(z) - z - lgamma(a);
+
+  T_return b_init = z + 1.0 - a;
+  T_return C = (fabs(value_of_rec(b_init)) >= EPSILON)
+                   ? b_init
+                   : std::decay_t<decltype(b_init)>(EPSILON);
+  T_return D = 0.0;
+  T_return f = C;
+  for (int i = 1; i <= max_steps; ++i) {
+    T_a an = -i * (i - a);
+    const T_return b = b_init + 2.0 * i;
+    D = b + an * D;
+    D = (fabs(value_of_rec(D)) >= EPSILON) ? D
+                                           : std::decay_t<decltype(D)>(EPSILON);
+    C = b + an / C;
+    C = (fabs(value_of_rec(C)) >= EPSILON) ? C
+                                           : std::decay_t<decltype(C)>(EPSILON);
+    D = inv(D);
+    const T_return delta = C * D;
+    f *= delta;
+    const double delta_m1 = fabs(value_of_rec(delta) - 1.0);
+    if (delta_m1 < precision) {
+      break;
+    }
+  }
+  return log_prefactor - log(f);
+}
+
+}  // namespace internal
+
+/**
+ * Compute log(Q(a,z)) and its gradient with respect to a using continued
+ * fraction expansion, where Q(a,z) = Gamma(a,z) / Gamma(a) is the regularized
+ * upper incomplete gamma function.
+ *
+ * This uses a continued fraction representation for numerical stability when
+ * computing the upper incomplete gamma function in log space, along with
+ * analytical gradient computation.
+ *
+ * @tparam T_a type of the shape parameter
+ * @tparam T_z type of the value parameter
+ * @param a shape parameter (must be positive)
+ * @param z value parameter (must be non-negative)
+ * @param precision convergence threshold, default of sqrt(machine_epsilon)
+ * @param max_steps maximum iterations for continued fraction
+ * @return structure containing log(Q(a,z)) and d/da log(Q(a,z))
+ */
+template <typename T_a, typename T_z>
+inline std::pair<return_type_t<T_a, T_z>, return_type_t<T_a, T_z>>
+log_gamma_q_dgamma(const T_a& a, const T_z& z,
+                   double precision = internal::LOG_Q_GAMMA_CF_PRECISION,
+                   int max_steps = 250) {
+  using T_return = return_type_t<T_a, T_z>;
+  const double a_val = value_of(a);
+  const double z_val = value_of(z);
+  // For z > a + 1, use continued fraction for better numerical stability
+  if (z_val > a_val + 1.0) {
+    std::pair<T_return, T_return> result{
+        internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0};
+    // For gradient, use: d/da log(Q) = (1/Q) * dQ/da
+    // grad_reg_inc_gamma computes dQ/da
+    const T_return Q_val = exp(result.first);
+    const double dQ_da
+        = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val));
+    result.second = dQ_da / Q_val;
+    return result;
+  } else {
+    // For z <= a + 1, use log1m(P(a,z)) for better numerical accuracy
+    const double P_val = gamma_p(a_val, z_val);
+    std::pair<T_return, T_return> result{log1m(P_val), 0.0};
+    // Gradient: d/da log(Q) = (1/Q) * dQ/da
+    // grad_reg_inc_gamma computes dQ/da
+    const T_return Q_val = exp(result.first);
+    if (Q_val > 0) {
+      const double dQ_da
+          = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val));
+      result.second = dQ_da / Q_val;
+    } else {
+      // Fallback if Q rounds to zero - use asymptotic approximation
+      result.second = log(z_val) - digamma(a_val);
+    }
+    return result;
+  }
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/prim/meta/is_fvar.hpp

@@ -14,6 +14,9 @@ namespace stan {
 template <typename T, typename = void>
 struct is_fvar : std::false_type {};
 
+template <typename T>
+inline constexpr bool is_fvar_v = is_fvar<T>::value;
+
 /** \ingroup type_trait
  * Specialization for pointers returns the underlying value the pointer is
  * pointing to.