Merge remote-tracking branch 'lingium/feature/issue-3113-beta-neg-binomial-lccdf' into HEAD

Author: SteveBronder

makefile

@@ -87,6 +87,7 @@ doxygen:
 clean:
 	@echo '  removing generated test files'
 	@$(RM) $(wildcard test/prob/generate_tests$(EXE))
+	@$(RM) $(EXPRESSION_TESTS) $(call findfiles,test/expressions,*_test.cpp)
 	@$(RM) $(call findfiles,test/prob,*_generated_v_test.cpp)
 	@$(RM) $(call findfiles,test/prob,*_generated_vv_test.cpp)
 	@$(RM) $(call findfiles,test/prob,*_generated_fd_test.cpp)

runTests.py

@@ -383,7 +383,7 @@ def checkToolchainPathWindows():
             universal_newlines=True,
         )
         out, err = p1.communicate()
-        if re.search(" |\(|\)", out):
+        if re.search(r" |\(|\)", out):
             stopErr(
                 "The RTools toolchain is installed in a path with spaces or bracket. Please reinstall to a valid path.",
                 -1,

stan/math/prim/fun/grad_F32.hpp

@@ -49,9 +49,9 @@ namespace math {
  * @param[in] max_steps number of steps to take
  */
 template <bool grad_a1 = true, bool grad_a2 = true, bool grad_a3 = true,
-  bool grad_b1 = true, bool grad_b2 = true, bool grad_z = true,
-  typename T1, typename T2, typename T3, typename T4, typename T5,
-  typename T6, typename T7, typename T8 = double>
+          bool grad_b1 = true, bool grad_b2 = true, bool grad_z = true,
+          typename T1, typename T2, typename T3, typename T4, typename T5,
+          typename T6, typename T7, typename T8 = double>
 void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, const T5& b1,
               const T6& b2, const T7& z, const T8& precision = 1e-6,
               int max_steps = 1e5) {
@@ -87,47 +87,53 @@ void grad_F32(T1* g, const T2& a1, const T3& a2, const T4& a3, const T5& b1,
     log_t_new += log(fabs(p)) + log_z;
     log_t_new_sign = p >= 0.0 ? log_t_new_sign : -log_t_new_sign;
     if constexpr (grad_a1) {
-      term[0] = log_g_old_sign[0] * log_t_old_sign * exp(log_g_old[0] - log_t_old)
-              + inv(a1 + k);
+      term[0]
+          = log_g_old_sign[0] * log_t_old_sign * exp(log_g_old[0] - log_t_old)
+            + inv(a1 + k);
       log_g_old[0] = log_t_new + log(fabs(term[0]));
       log_g_old_sign[0] = term[0] >= 0.0 ? log_t_new_sign : -log_t_new_sign;
       g[0] += log_g_old_sign[0] * exp(log_g_old[0]);
     }
 
     if constexpr (grad_a2) {
-      term[1] = log_g_old_sign[1] * log_t_old_sign * exp(log_g_old[1] - log_t_old)
+      term[1]
+          = log_g_old_sign[1] * log_t_old_sign * exp(log_g_old[1] - log_t_old)
             + inv(a2 + k);
       log_g_old[1] = log_t_new + log(fabs(term[1]));
       log_g_old_sign[1] = term[1] >= 0.0 ? log_t_new_sign : -log_t_new_sign;
       g[1] += log_g_old_sign[1] * exp(log_g_old[1]);
     }
 
     if constexpr (grad_a3) {
-      term[2] = log_g_old_sign[2] * log_t_old_sign * exp(log_g_old[2] - log_t_old)
+      term[2]
+          = log_g_old_sign[2] * log_t_old_sign * exp(log_g_old[2] - log_t_old)
             + inv(a3 + k);
       log_g_old[2] = log_t_new + log(fabs(term[2]));
       log_g_old_sign[2] = term[2] >= 0.0 ? log_t_new_sign : -log_t_new_sign;
       g[2] += log_g_old_sign[2] * exp(log_g_old[2]);
     }
 
     if constexpr (grad_b1) {
-      term[3] = log_g_old_sign[3] * log_t_old_sign * exp(log_g_old[3] - log_t_old)
+      term[3]
+          = log_g_old_sign[3] * log_t_old_sign * exp(log_g_old[3] - log_t_old)
             - inv(b1 + k);
       log_g_old[3] = log_t_new + log(fabs(term[3]));
       log_g_old_sign[3] = term[3] >= 0.0 ? log_t_new_sign : -log_t_new_sign;
       g[3] += log_g_old_sign[3] * exp(log_g_old[3]);
     }
 
     if constexpr (grad_b2) {
-      term[4] = log_g_old_sign[4] * log_t_old_sign * exp(log_g_old[4] - log_t_old)
+      term[4]
+          = log_g_old_sign[4] * log_t_old_sign * exp(log_g_old[4] - log_t_old)
             - inv(b2 + k);
       log_g_old[4] = log_t_new + log(fabs(term[4]));
       log_g_old_sign[4] = term[4] >= 0.0 ? log_t_new_sign : -log_t_new_sign;
       g[4] += log_g_old_sign[4] * exp(log_g_old[4]);
     }
 
     if constexpr (grad_z) {
-      term[5] = log_g_old_sign[5] * log_t_old_sign * exp(log_g_old[5] - log_t_old)
+      term[5]
+          = log_g_old_sign[5] * log_t_old_sign * exp(log_g_old[5] - log_t_old)
             + inv(z);
       log_g_old[5] = log_t_new + log(fabs(term[5]));
       log_g_old_sign[5] = term[5] >= 0.0 ? log_t_new_sign : -log_t_new_sign;

stan/math/prim/fun/grad_pFq.hpp

@@ -89,10 +89,11 @@ template <bool calc_a = true, bool calc_b = true, bool calc_z = true,
           typename T_Rtn = return_type_t<Ta, Tb, Tz>,
           typename Ta_Rtn = promote_scalar_t<T_Rtn, plain_type_t<Ta>>,
           typename Tb_Rtn = promote_scalar_t<T_Rtn, plain_type_t<Tb>>>
-inline std::tuple<Ta_Rtn, Tb_Rtn, T_Rtn> grad_pFq(const TpFq& pfq_val, const Ta& a,
-                                           const Tb& b, const Tz& z,
-                                           double precision = 1e-14,
-                                           int max_steps = 1e6) {
+inline std::tuple<Ta_Rtn, Tb_Rtn, T_Rtn> grad_pFq(const TpFq& pfq_val,
+                                                  const Ta& a, const Tb& b,
+                                                  const Tz& z,
+                                                  double precision = 1e-14,
+                                                  int max_steps = 1e6) {
   using std::max;
   using Ta_Array = Eigen::Array<return_type_t<Ta>, -1, 1>;
   using Tb_Array = Eigen::Array<return_type_t<Tb>, -1, 1>;

stan/math/prim/fun/hypergeometric_3F2.hpp

@@ -19,12 +19,13 @@ namespace internal {
 template <typename Ta, typename Tb, typename Tz,
           require_all_vector_t<Ta, Tb>* = nullptr,
           require_stan_scalar_t<Tz>* = nullptr>
-inline return_type_t<Ta, Tb, Tz> hypergeometric_3F2_infsum(const Ta& a, const Tb& b, const Tz& z,
-                                   double precision = 1e-6,
-                                   int max_steps = 1e5) {
+inline return_type_t<Ta, Tb, Tz> hypergeometric_3F2_infsum(
+    const Ta& a, const Tb& b, const Tz& z, double precision = 1e-6,
+    int max_steps = 1e5) {
   using T_return = return_type_t<Ta, Tb, Tz>;
   Eigen::Array<scalar_type_t<Ta>, 3, 1> a_array = as_array_or_scalar(a);
-  Eigen::Array<scalar_type_t<Tb>, 3, 1> b_array = append_row(as_array_or_scalar(b), 1.0);
+  Eigen::Array<scalar_type_t<Tb>, 3, 1> b_array
+      = append_row(as_array_or_scalar(b), 1.0);
   check_3F2_converges("hypergeometric_3F2", a_array[0], a_array[1], a_array[2],
                       b_array[0], b_array[1], z);
 
@@ -141,7 +142,8 @@ inline auto hypergeometric_3F2(const Ta& a, const Tb& b, const Tz& z) {
 template <typename Ta, typename Tb, typename Tz,
           require_all_stan_scalar_t<Ta, Tb, Tz>* = nullptr>
 inline auto hypergeometric_3F2(const std::initializer_list<Ta>& a,
-                        const std::initializer_list<Tb>& b, const Tz& z) {
+                               const std::initializer_list<Tb>& b,
+                               const Tz& z) {
   return hypergeometric_3F2(std::vector<Ta>(a), std::vector<Tb>(b), z);
 }
 

stan/math/prim/prob/beta_neg_binomial_lccdf.hpp

@@ -61,7 +61,6 @@ inline return_type_t<T_r, T_alpha, T_beta> beta_neg_binomial_lccdf(
   check_positive_finite(function, "Prior success parameter", alpha_ref);
   check_positive_finite(function, "Prior failure parameter", beta_ref);
 
-
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_r_ref> r_vec(r_ref);
   scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
@@ -94,24 +93,23 @@ inline return_type_t<T_r, T_alpha, T_beta> beta_neg_binomial_lccdf(
     auto r_plus_n = r_dbl + n_dbl;
     auto a_plus_r = alpha_dbl + r_dbl;
     using a_t = return_type_t<decltype(b_plus_n), decltype(r_plus_n)>;
-    using b_t = return_type_t<decltype(n_dbl), decltype(a_plus_r), decltype(b_plus_n)>;
-    auto F
-        = hypergeometric_3F2(
-          std::initializer_list<a_t>{1.0, b_plus_n + 1.0, r_plus_n + 1.0},
-          std::initializer_list<b_t>{n_dbl + 2.0, a_plus_r + b_plus_n + 1.0}, 1.0);
+    using b_t = return_type_t<decltype(n_dbl), decltype(a_plus_r),
+                              decltype(b_plus_n)>;
+    auto F = hypergeometric_3F2(
+        std::initializer_list<a_t>{1.0, b_plus_n + 1.0, r_plus_n + 1.0},
+        std::initializer_list<b_t>{n_dbl + 2.0, a_plus_r + b_plus_n + 1.0},
+        1.0);
     auto C = lgamma(r_plus_n + 1.0) + lbeta(a_plus_r, b_plus_n + 1.0)
-                          - lgamma(r_dbl) - lbeta(alpha_dbl, beta_dbl)
-                          - lgamma(n_dbl + 2);
+             - lgamma(r_dbl) - lbeta(alpha_dbl, beta_dbl) - lgamma(n_dbl + 2);
     log_ccdf += C + stan::math::log(F);
 
     if constexpr (!is_constant_all<T_r, T_alpha, T_beta>::value) {
-      auto digamma_n_r_alpha_beta
-          = digamma(a_plus_r + b_plus_n + 1.0);
+      auto digamma_n_r_alpha_beta = digamma(a_plus_r + b_plus_n + 1.0);
       T_partials_return dF[6];
-      grad_F32<false, !is_constant<T_beta>::value,
-        !is_constant_all<T_r>::value, false, true, false>(dF, 1.0,
-          b_plus_n + 1.0, r_plus_n + 1.0, n_dbl + 2.0,
-          a_plus_r + b_plus_n + 1.0, 1.0, precision, max_steps);
+      grad_F32<false, !is_constant<T_beta>::value, !is_constant_all<T_r>::value,
+               false, true, false>(dF, 1.0, b_plus_n + 1.0, r_plus_n + 1.0,
+                                   n_dbl + 2.0, a_plus_r + b_plus_n + 1.0, 1.0,
+                                   precision, max_steps);
 
       if constexpr (!is_constant<T_r>::value || !is_constant<T_alpha>::value) {
         auto digamma_r_alpha = digamma(a_plus_r);

stan/math/prim/prob/wiener5_lpdf.hpp

@@ -679,12 +679,12 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0,
   if (!include_summand<propto, T_y, T_a, T_t0, T_w, T_v, T_sv>::value) {
     return ret_t(0.0);
   }
-  using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
-  using T_a_ref = ref_type_if_t<!is_constant<T_a>::value, T_a>;
-  using T_t0_ref = ref_type_if_t<!is_constant<T_t0>::value, T_t0>;
-  using T_w_ref = ref_type_if_t<!is_constant<T_w>::value, T_w>;
-  using T_v_ref = ref_type_if_t<!is_constant<T_v>::value, T_v>;
-  using T_sv_ref = ref_type_if_t<!is_constant<T_sv>::value, T_sv>;
+  using T_y_ref = ref_type_t<T_y>;
+  using T_a_ref = ref_type_t<T_a>;
+  using T_t0_ref = ref_type_t<T_t0>;
+  using T_w_ref = ref_type_t<T_w>;
+  using T_v_ref = ref_type_t<T_v>;
+  using T_sv_ref = ref_type_t<T_sv>;
 
   static constexpr const char* function_name = "wiener5_lpdf";
 

stan/math/prim/prob/wiener_full_lpdf.hpp

@@ -327,14 +327,14 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0,
     return ret_t(0);
   }
 
-  using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
-  using T_a_ref = ref_type_if_t<!is_constant<T_a>::value, T_a>;
-  using T_v_ref = ref_type_if_t<!is_constant<T_v>::value, T_v>;
-  using T_w_ref = ref_type_if_t<!is_constant<T_w>::value, T_w>;
-  using T_t0_ref = ref_type_if_t<!is_constant<T_t0>::value, T_t0>;
-  using T_sv_ref = ref_type_if_t<!is_constant<T_sv>::value, T_sv>;
-  using T_sw_ref = ref_type_if_t<!is_constant<T_sw>::value, T_sw>;
-  using T_st0_ref = ref_type_if_t<!is_constant<T_st0>::value, T_st0>;
+  using T_y_ref = ref_type_t<T_y>;
+  using T_a_ref = ref_type_t<T_a>;
+  using T_v_ref = ref_type_t<T_v>;
+  using T_w_ref = ref_type_t<T_w>;
+  using T_t0_ref = ref_type_t<T_t0>;
+  using T_sv_ref = ref_type_t<T_sv>;
+  using T_sw_ref = ref_type_t<T_sw>;
+  using T_st0_ref = ref_type_t<T_st0>;
 
   using T_partials_return
       = partials_return_t<T_y, T_a, T_t0, T_w, T_v, T_sv, T_sw, T_st0>;
@@ -449,6 +449,9 @@ inline auto wiener_lpdf(const T_y& y, const T_a& a, const T_t0& t0,
   // calculate density and partials
   for (size_t i = 0; i < N; i++) {
     if (sw_vec[i] == 0 && st0_vec[i] == 0) {
+      // note: because we're delegating to wiener5_lpdf,
+      // we need to make sure is_constant is consistent between
+      // our inputs and these
       result += wiener_lpdf<propto>(y_vec[i], a_vec[i], t0_vec[i], w_vec[i],
                                     v_vec[i], sv_vec[i], precision_derivatives);
       continue;

test/generate_expression_tests.py

@@ -29,6 +29,9 @@ def save_tests_in_files(N_files, tests):
     for i in range(N_files):
         start = i * len(tests) // N_files
         end = (i + 1) * len(tests) // N_files
+        if start >= end:
+            # don't try to compile an empty file
+            continue
         with open(src_folder + "tests%d_test.cpp" % i, "w") as out:
             out.write("#include <test/expressions/expression_test_helpers.hpp>\n\n")
             for test in tests[start:end]:
@@ -125,5 +128,5 @@ def main(functions=(), j=1):
                     code = cg.cpp(),
                 )
             )
-    
+
     save_tests_in_files(j, tests)

test/sig_utils.py

@@ -132,7 +132,7 @@ def get_cpp_type(stan_type):
     "uniform_lcdf": [None, 0.2, 0.9],
     "uniform_lpdf": [None, 0.2, 0.9],
     "uniform_rng": [0.2, 1.9, None],
-    "wiener_lpdf": [0.8, None, 0.4, None, None],
+    "wiener_lpdf": [0.8, None, 0.4, None, None, None, None, None],
 }
 
 # list of functions we do not test. These are mainly functions implemented in compiler