From 81c67a0e33b68936d9f3eb1551ab1cb3d54dda27 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Tue, 13 Jan 2026 15:55:27 -0500 Subject: [PATCH 1/3] refactor gamma_lccdf --- stan/math/prim/prob/gamma_lccdf.hpp | 176 +++++++++++++++------------- 1 file changed, 97 insertions(+), 79 deletions(-) diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index ad685ab2137..aab4590dd96 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -18,21 +18,95 @@ #include #include #include -#include +#include #include #include #include namespace stan { namespace math { +namespace internal { + template + struct Q_eval { + T log_Q{0.0}; + T dlogQ_dalpha{0.0}; + bool ok{false}; + }; + + /** + * Computes log q and d(log q) / d(alpha) using continued fraction. + */ + template + static inline Q_eval eval_q_cf(const T& alpha, + const T& beta_y) { + Q_eval out; + if constexpr (!any_fvar && is_autodiff_v) { + auto log_q_result = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); + out.log_Q = log_q_result.log_q; + out.dlogQ_dalpha = log_q_result.dlog_q_da; + } else { + out.log_Q = internal::log_q_gamma_cf(alpha, beta_y); + if constexpr (is_autodiff_v) { + if constexpr (!partials_fvar) { + out.dlogQ_dalpha + = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), + digamma(alpha)) / exp(out.log_Q); + } else { + T alpha_unit = alpha; + alpha_unit.d_ = 1; + T beta_y_unit = beta_y; + beta_y_unit.d_ = 0; + T log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); + out.dlogQ_dalpha = log_Q_fvar.d_; + } + } + } + + out.ok = std::isfinite(value_of_rec(out.log_Q)); + return out; + } + + /** + * Computes log q and d(log q) / d(alpha) using log1m. + */ + template + static inline Q_eval eval_q_log1m(const T& alpha, + const T& beta_y) { + Q_eval out; + out.log_Q = log1m(gamma_p(alpha, beta_y)); + + if (!std::isfinite(value_of_rec(out.log_Q))) { + out.ok = false; + return out; + } + + if constexpr (is_autodiff_v) { + if constexpr (partials_fvar) { + T alpha_unit = alpha; + alpha_unit.d_ = 1; + T beta_unit = beta_y; + beta_unit.d_ = 0; + T log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); + out.dlogQ_dalpha = log_Q_fvar.d_; + } else { + out.dlogQ_dalpha = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.log_Q); + } + } + + out.ok = true; + return out; + } +} template -return_type_t gamma_lccdf(const T_y& y, - const T_shape& alpha, - const T_inv_scale& beta) { - using T_partials_return = partials_return_t; +inline return_type_t gamma_lccdf(const T_y& y, + const T_shape& alpha, + const T_inv_scale& beta) { using std::exp; using std::log; + using T_partials_return = partials_return_t; using T_y_ref = ref_type_t; using T_alpha_ref = ref_type_t; using T_beta_ref = ref_type_t; @@ -58,7 +132,6 @@ return_type_t gamma_lccdf(const T_y& y, scalar_seq_view beta_vec(beta_ref); const size_t N = max_size(y, alpha, beta); - constexpr bool need_y_beta_deriv = !is_constant_all::value; constexpr bool any_fvar = is_fvar>::value || is_fvar>::value || is_fvar>::value; @@ -83,90 +156,35 @@ return_type_t gamma_lccdf(const T_y& y, return ops_partials.build(negative_infinity()); } - bool use_cf = beta_y > alpha_dbl + 1.0; - T_partials_return log_Qn; - [[maybe_unused]] T_partials_return dlogQ_dalpha = 0.0; - - // Branch by autodiff type first, then handle use_cf logic inside each path - if constexpr (!any_fvar && is_autodiff_v) { - // var-only path: use log_gamma_q_dgamma which computes both log_q - // and its gradient analytically with double inputs - const double beta_y_dbl = value_of(value_of(beta_y)); - const double alpha_dbl_val = value_of(value_of(alpha_dbl)); - - if (use_cf) { - auto log_q_result = log_gamma_q_dgamma(alpha_dbl_val, beta_y_dbl); - log_Qn = log_q_result.log_q; - dlogQ_dalpha = log_q_result.dlog_q_da; - } else { - const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); - log_Qn = log1m(Pn); - const T_partials_return Qn = exp(log_Qn); - - // Check if we need to fallback to continued fraction - bool need_cf_fallback - = !std::isfinite(value_of(value_of(log_Qn))) || Qn <= 0.0; - if (need_cf_fallback && beta_y > 0.0) { - auto log_q_result = log_gamma_q_dgamma(alpha_dbl_val, beta_y_dbl); - log_Qn = log_q_result.log_q; - dlogQ_dalpha = log_q_result.dlog_q_da; - } else { - dlogQ_dalpha = -grad_reg_lower_inc_gamma(alpha_dbl, beta_y) / Qn; - } - } - } else if constexpr (partials_fvar && is_autodiff_v) { - // fvar path: use unit derivative trick to compute gradients - auto alpha_unit = alpha_dbl; - alpha_unit.d_ = 1; - auto beta_unit = beta_y; - beta_unit.d_ = 0; - - if (use_cf) { - log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); - auto log_Qn_fvar = internal::log_q_gamma_cf(alpha_unit, beta_unit); - dlogQ_dalpha = log_Qn_fvar.d_; - } else { - const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); - log_Qn = log1m(Pn); - - if (!std::isfinite(value_of(value_of(log_Qn))) && beta_y > 0.0) { - // Fallback to continued fraction - log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); - auto log_Qn_fvar = internal::log_q_gamma_cf(alpha_unit, beta_unit); - dlogQ_dalpha = log_Qn_fvar.d_; - } else { - auto log_Qn_fvar = log1m(gamma_p(alpha_unit, beta_unit)); - dlogQ_dalpha = log_Qn_fvar.d_; - } - } + const bool use_continued_fraction = beta_y > alpha_dbl + 1.0; + internal::Q_eval result; + if (use_continued_fraction) { + result = internal::eval_q_cf(alpha_dbl, beta_y); } else { - // No alpha derivative needed (alpha is constant or double-only) - if (use_cf) { - log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); - } else { - const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); - log_Qn = log1m(Pn); + result = internal::eval_q_log1m(alpha_dbl, beta_y); - if (!std::isfinite(value_of(value_of(log_Qn))) && beta_y > 0.0) { - log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); - } + if (!result.ok && beta_y > 0.0) { + // Fallback to continued fraction if log1m fails + result = internal::eval_q_cf(alpha_dbl, beta_y); } } - if (!std::isfinite(value_of(value_of(log_Qn)))) { + if (!result.ok) { return ops_partials.build(negative_infinity()); } - P += log_Qn; - if constexpr (need_y_beta_deriv) { + P += result.log_Q; + + if constexpr (is_autodiff_v || is_autodiff_v) { const T_partials_return log_y = log(y_dbl); - const T_partials_return log_beta = log(beta_dbl); - const T_partials_return lgamma_alpha = lgamma(alpha_dbl); const T_partials_return alpha_minus_one = fma(alpha_dbl, log_y, -log_y); const T_partials_return log_pdf - = alpha_dbl * log_beta - lgamma_alpha + alpha_minus_one - beta_y; + = alpha_dbl * log(beta_dbl) - lgamma(alpha_dbl) + alpha_minus_one - beta_y; - const T_partials_return hazard = exp(log_pdf - log_Qn); // f/Q + const T_partials_return hazard = exp(log_pdf - result.log_Q); // f/Q if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] -= hazard; @@ -176,7 +194,7 @@ return_type_t gamma_lccdf(const T_y& y, } } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] += dlogQ_dalpha; + partials<1>(ops_partials)[n] += result.dlogQ_dalpha; } } return ops_partials.build(P); From df0101ed85999d609fe182bf0549edc559dac64c Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 13 Jan 2026 15:58:02 -0500 Subject: [PATCH 2/3] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/prob/gamma_lccdf.hpp | 147 ++++++++++++++-------------- 1 file changed, 73 insertions(+), 74 deletions(-) diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index aab4590dd96..6c7a0ebd9d2 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -26,84 +26,81 @@ namespace stan { namespace math { namespace internal { - template - struct Q_eval { - T log_Q{0.0}; - T dlogQ_dalpha{0.0}; - bool ok{false}; - }; - - /** - * Computes log q and d(log q) / d(alpha) using continued fraction. - */ - template - static inline Q_eval eval_q_cf(const T& alpha, - const T& beta_y) { - Q_eval out; - if constexpr (!any_fvar && is_autodiff_v) { - auto log_q_result = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); - out.log_Q = log_q_result.log_q; - out.dlogQ_dalpha = log_q_result.dlog_q_da; - } else { - out.log_Q = internal::log_q_gamma_cf(alpha, beta_y); - if constexpr (is_autodiff_v) { - if constexpr (!partials_fvar) { - out.dlogQ_dalpha - = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), - digamma(alpha)) / exp(out.log_Q); - } else { - T alpha_unit = alpha; - alpha_unit.d_ = 1; - T beta_y_unit = beta_y; - beta_y_unit.d_ = 0; - T log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); - out.dlogQ_dalpha = log_Q_fvar.d_; - } - } - } - - out.ok = std::isfinite(value_of_rec(out.log_Q)); - return out; - } - - /** - * Computes log q and d(log q) / d(alpha) using log1m. - */ - template - static inline Q_eval eval_q_log1m(const T& alpha, - const T& beta_y) { - Q_eval out; - out.log_Q = log1m(gamma_p(alpha, beta_y)); - - if (!std::isfinite(value_of_rec(out.log_Q))) { - out.ok = false; - return out; - } - +template +struct Q_eval { + T log_Q{0.0}; + T dlogQ_dalpha{0.0}; + bool ok{false}; +}; + +/** + * Computes log q and d(log q) / d(alpha) using continued fraction. + */ +template +static inline Q_eval eval_q_cf(const T& alpha, const T& beta_y) { + Q_eval out; + if constexpr (!any_fvar && is_autodiff_v) { + auto log_q_result + = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); + out.log_Q = log_q_result.log_q; + out.dlogQ_dalpha = log_q_result.dlog_q_da; + } else { + out.log_Q = internal::log_q_gamma_cf(alpha, beta_y); if constexpr (is_autodiff_v) { - if constexpr (partials_fvar) { + if constexpr (!partials_fvar) { + out.dlogQ_dalpha + = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), digamma(alpha)) + / exp(out.log_Q); + } else { T alpha_unit = alpha; alpha_unit.d_ = 1; - T beta_unit = beta_y; - beta_unit.d_ = 0; - T log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); + T beta_y_unit = beta_y; + beta_y_unit.d_ = 0; + T log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); out.dlogQ_dalpha = log_Q_fvar.d_; - } else { - out.dlogQ_dalpha = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.log_Q); } } + } + + out.ok = std::isfinite(value_of_rec(out.log_Q)); + return out; +} - out.ok = true; +/** + * Computes log q and d(log q) / d(alpha) using log1m. + */ +template +static inline Q_eval eval_q_log1m(const T& alpha, const T& beta_y) { + Q_eval out; + out.log_Q = log1m(gamma_p(alpha, beta_y)); + + if (!std::isfinite(value_of_rec(out.log_Q))) { + out.ok = false; return out; } + + if constexpr (is_autodiff_v) { + if constexpr (partials_fvar) { + T alpha_unit = alpha; + alpha_unit.d_ = 1; + T beta_unit = beta_y; + beta_unit.d_ = 0; + T log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); + out.dlogQ_dalpha = log_Q_fvar.d_; + } else { + out.dlogQ_dalpha + = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.log_Q); + } + } + + out.ok = true; + return out; } +} // namespace internal template -inline return_type_t gamma_lccdf(const T_y& y, - const T_shape& alpha, - const T_inv_scale& beta) { +inline return_type_t gamma_lccdf( + const T_y& y, const T_shape& alpha, const T_inv_scale& beta) { using std::exp; using std::log; using T_partials_return = partials_return_t; @@ -159,16 +156,17 @@ inline return_type_t gamma_lccdf(const T_y& y, const bool use_continued_fraction = beta_y > alpha_dbl + 1.0; internal::Q_eval result; if (use_continued_fraction) { - result = internal::eval_q_cf(alpha_dbl, beta_y); - } else { - result = internal::eval_q_log1m(alpha_dbl, beta_y); + } else { + result + = internal::eval_q_log1m( + alpha_dbl, beta_y); if (!result.ok && beta_y > 0.0) { // Fallback to continued fraction if log1m fails - result = internal::eval_q_cf(alpha_dbl, beta_y); + result = internal::eval_q_cf(alpha_dbl, beta_y); } } if (!result.ok) { @@ -181,8 +179,9 @@ inline return_type_t gamma_lccdf(const T_y& y, const T_partials_return log_y = log(y_dbl); const T_partials_return alpha_minus_one = fma(alpha_dbl, log_y, -log_y); - const T_partials_return log_pdf - = alpha_dbl * log(beta_dbl) - lgamma(alpha_dbl) + alpha_minus_one - beta_y; + const T_partials_return log_pdf = alpha_dbl * log(beta_dbl) + - lgamma(alpha_dbl) + alpha_minus_one + - beta_y; const T_partials_return hazard = exp(log_pdf - result.log_Q); // f/Q From ac6d518c4cdafca712ef0a1cc154ac627267bc3a Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Thu, 5 Feb 2026 10:20:43 -0500 Subject: [PATCH 3/3] update --- stan/math/prim/prob/gamma_lccdf.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index 67c0e209c66..c03204dbb6f 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -14,12 +14,8 @@ #include #include #include -<<<<<<< HEAD #include #include -======= -#include ->>>>>>> develop #include #include