From 8a8ea67262593135af803c57111c628e7c496a0f Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Tue, 14 Jan 2020 11:36:33 +0100 Subject: [PATCH 01/26] Failing test --- .../fun/binomial_coefficient_log_test.cpp | 65 ++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp index 23cea698637..b4b9fb9b41b 100644 --- a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -7,7 +8,8 @@ template void test_binom_coefficient(const T_N& N, const T_n& n) { using stan::math::binomial_coefficient_log; EXPECT_FLOAT_EQ(lgamma(N + 1) - lgamma(n + 1) - lgamma(N - n + 1), - binomial_coefficient_log(N, n)); + binomial_coefficient_log(N, n)) + << "N = " << N << ", n = " << n; } TEST(MathFunctions, binomial_coefficient_log) { @@ -19,6 +21,15 @@ TEST(MathFunctions, binomial_coefficient_log) { EXPECT_FLOAT_EQ(29979.16, binomial_coefficient_log(100000, 91116)); + EXPECT_EQ(binomial_coefficient_log(-1, 0), 0); // Needed for neg_binomial_2 + EXPECT_EQ(binomial_coefficient_log(50, 0), 0); + EXPECT_EQ(binomial_coefficient_log(10000, 0), 0); + + EXPECT_EQ(binomial_coefficient_log(10, 11), + -std::numeric_limits::infinity()); + EXPECT_EQ(binomial_coefficient_log(10, -1), + -std::numeric_limits::infinity()); + for (int n = 0; n < 1010; ++n) { test_binom_coefficient(1010, n); test_binom_coefficient(1010.0, n); @@ -38,3 +49,55 @@ TEST(MathFunctions, binomial_coefficient_log_nan) { EXPECT_TRUE(std::isnan(stan::math::binomial_coefficient_log(nan, 2.0))); EXPECT_TRUE(std::isnan(stan::math::binomial_coefficient_log(nan, nan))); } + +TEST(MathFunctions, binomial_coefficient_log_errors_edge_cases) { + using stan::math::binomial_coefficient_log; + double inf = std::numeric_limits::infinity(); + EXPECT_NO_THROW(binomial_coefficient_log(10, 11)); + EXPECT_THROW(binomial_coefficient_log(10, 11.01), std::domain_error); + EXPECT_THROW(binomial_coefficient_log(10, -1.1), std::domain_error); + EXPECT_THROW(binomial_coefficient_log(-1, 0.3), std::domain_error); + EXPECT_NO_THROW(binomial_coefficient_log(-0.5, 0.49)); + EXPECT_NO_THROW(binomial_coefficient_log(10, -0.9)); + + EXPECT_FLOAT_EQ(binomial_coefficient_log(0, -1), -inf); + EXPECT_FLOAT_EQ(binomial_coefficient_log(-1, 0), 0); + EXPECT_FLOAT_EQ(binomial_coefficient_log(-1, -0.3), inf); + EXPECT_FLOAT_EQ(binomial_coefficient_log(0.3, -1), -inf); + +} + +namespace binomial_coefficient_test_internal { +struct TestValue { + double n; + double k; + double val; +}; + +std::vector testValues = { {-0.1,-0.9,-2.1152525390850903}, {-0.1,-1.0000000000000001e-11,1.7771062399418724e-12}, {-0.1,-1.0000000000000002e-6,1.7770950132278696e-7}, {-0.1,-0.001,0.00017592768030333383}, {-0.1,-0.020000000000000004,0.0028416866729211523}, {-0.1,-0.05,0.0044386492587971844}, {-0.1,-0.097,0.0005170837568402953}, {-0.1,-0.09999999400000001,1.0662676433991578e-9}, {-0.1,-0.09999999970000001,5.331331690149795e-11}, {-0.1,0.8,-2.1152525390850907}, {0.00003,-0.9,-2.21375637737528044964112}, {0.00003,3.e-15,1.48040820535563428588846e-19}, {0.00003,3.e-10,1.48039340142161568666212e-14}, {0.00003,3.e-7,1.46560412344434242311062e-11}, {0.00003,6.e-6,2.36865312869367134774661e-10}, {0.00003,0.000015,3.70102051348524044535716e-10}, {0.00003,0.0000291,4.30798787797857754693695e-11}, {0.00003,0.0000299999982,8.88244870007509649977929e-17}, {0.00003,0.00002999999991,4.44122460318735146597007e-18}, {0.00003,0.90003,-2.21375637737528044964112}, {0.002,-0.9,-2.21559326412971099686943}, {0.002,2.e-13,6.57013709556564677684856e-16}, {0.002,2.e-8,6.57007139476544583161173e-11}, {0.002,0.00002,6.50443564072189575550994e-8}, {0.002,0.0004,1.05122171458287350859763e-6}, {0.002,0.001,1.64253373496215313253469e-6}, {0.002,0.00194,1.91190982195918976356429e-7}, {0.002,0.00199999988,3.94208202120835082737684e-13}, {0.002,0.001999999994,1.97104112295366725515452e-14}, {0.002,0.902,-2.21559326412971099686943}, {1,-0.9,-2.85558226198351740582195}, {1,1.e-10,9.9999999988550659331851e-11}, {1,0.00001,9.99988550692664559909352e-6}, {1,0.01,0.00988583703486131052627978}, {1,0.2,0.156457962917688016707705}, {1,0.5,0.241564475270490444691037}, {1,0.97,0.028978328236256312960776}, {1,0.99999994,5.99999958782374313463811e-8}, {1,0.999999997,2.99999998969559340736596e-9}, {1,1.9,-2.85558226198351740582195}, {8,-0.9,-4.22528965320883461943031}, {8,8.e-10,2.17428571372173153982474e-9}, {8,0.00008,0.000217422931805073420417006}, {8,0.08,0.211982267378255838975509}, {8,1.6,2.90678606291134283426263}, {8,4,4.24849524204935898912334}, {8,7.76,0.606274586245453651115361}, {8,7.99999952,1.30457122553768403613331e-6}, {8,7.999999976,6.52285709209869625945566e-8}, {8,8.9,-4.22528965320883461943031}, {1325,-0.9,-8.72360867216657209762532}, {1325,1.325e-7,1.02909578020477960435539e-6}, {1325,0.01325,0.102766042691370430370992}, {1325,13.25,71.9898321274090629975055}, {1325,265,659.435649329029419323398}, {1325,662.5,914.599450340845275100724}, {1325,1285.25,175.786260651191862665015}, {1325,1324.9999205,0.000617452276410452190170437}, {1325,1324.999996025,0.0000308728608380968862741097}, {1325,1325.9,-8.72360867216657209762532}, {845000,-0.9,-14.5350963792733464918229}, {845000,0.0000845,0.00120194816738712136581358}, {845000,8.45,103.738303827736743600251}, {845000,8450,47315.8616457576200611209}, {845000,169000,422833.221695496506553128}, {845000,422500,585702.318235552114086514}, {845000,819650,113851.158132678562120685}, {845000,844999.9493,0.719108776819481762797449}, {845000,844999.997465,0.036053342347290003917417}, {845000,845000.9,-14.5350963792733464918229}, {3000000000000000.0,3,105.120406581508328854183}, {3000000000000000.0,100,3199.99949280231435502243}, {3000000000000000.0,12895,350387.5243605883687667}, {100000000000000000000.0,3,136.363346110414686040237}, {100000000000000000000.0,100,4241.4308104325278778424}, {100000000000000000000.0,12895,484680.092769031900296878},}; +} // namespace binomial_coefficient_test_internal + +TEST(MathFunctions, binomial_coefficient_log_precomputed) { + using binomial_coefficient_test_internal::TestValue; + using binomial_coefficient_test_internal::testValues; + using stan::test::expect_near_rel; + + for (TestValue t : testValues) { + std::ostringstream msg; + msg << std::setprecision(22) << "n = " << t.n << ", k = " << t.k; + + double val = stan::math::binomial_coefficient_log(t.n, t.k); + expect_near_rel(msg.str(), val, t.val); + + if (t.k > t.n / 10) { + // Testing the mirrored binomial coefficient is not performed for + // small k, as we may loose so much precision computing k2 + // that the test becomes invalid + std::ostringstream msg2; + double k2 = t.n - t.k; + msg2 << std::setprecision(22) << "n = " << t.n << ", k = " << k2; + double val2 = stan::math::binomial_coefficient_log(t.n, k2); + expect_near_rel(msg2.str(), val2, t.val); + } + } +} From 7991f40b478dc0766325153ef93bdec4152a6798 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Tue, 14 Jan 2020 11:36:47 +0100 Subject: [PATCH 02/26] Fixes #1592 --- .../prim/fun/binomial_coefficient_log.hpp | 51 ++++++++++++++----- 1 file changed, 38 insertions(+), 13 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 89600a35c55..f825a229f37 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -1,13 +1,17 @@ -#ifndef STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP -#define STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP +#ifndef STAN_MATH_PRIM_SCAL_FUN_BINOMIAL_COEFFICIENT_LOG_HPP +#define STAN_MATH_PRIM_SCAL_FUN_BINOMIAL_COEFFICIENT_LOG_HPP +#include #include -#include #include -#include +#include +#include +#include +#include namespace stan { namespace math { + /** * Return the log of the binomial coefficient for the specified * arguments. @@ -23,11 +27,15 @@ namespace math { * \f$ \log {N \choose n} * = \log \ \Gamma(N+1) - \log \Gamma(n+1) - \log \Gamma(N-n+1)\f$. * + * + * TODO[martinmodrak] figure out the cases for x < 0 and for partials \f[ \mbox{binomial\_coefficient\_log}(x, y) = \begin{cases} - \textrm{error} & \mbox{if } y > x \textrm{ or } y < 0\\ - \ln\Gamma(x+1) & \mbox{if } 0\leq y \leq x \\ + \textrm{error} & \mbox{if } y > x + 1 \textrm{ or } y < -1 \textrm{ or } x + < -1\\ + \textrm{-\infty} & \mbox{if } y = x + 1 \textrm{ or } y = -1\\ + \ln\Gamma(x+1) & \mbox{if } -1 < y < x + 1 \\ \quad -\ln\Gamma(y+1)& \\ \quad -\ln\Gamma(x-y+1)& \\[6pt] \textrm{NaN} & \mbox{if } x = \textrm{NaN or } y = \textrm{NaN} @@ -54,22 +62,39 @@ namespace math { \end{cases} \f] * + * This function is numerically more stable than naive evaluation via lgamma + * * @param N total number of objects. * @param n number of objects chosen. * @return log (N choose n). */ + template inline return_type_t binomial_coefficient_log(const T_N N, const T_n n) { - const double CUTOFF = 1000; - if (N - n < CUTOFF) { - const T_N N_plus_1 = N + 1; + if (is_nan(value_of_rec(N)) || is_nan(value_of_rec(n))) { + return std::numeric_limits::quiet_NaN(); + } + + // For some uses it is important this works even when N < 0 and therefore + // it is before checks + if (n == 0) { + return 0; + } + const T_N N_plus_1 = N + 1; + + static const char* function = "binomial_coefficient_log"; + check_greater_or_equal(function, "first argument", N, -1); + check_greater_or_equal(function, "second argument", n, -1); + check_greater_or_equal(function, "(first argument - second argument + 1)", + N - n + 1, 0.0); + + if (N / 2 < n) { + return binomial_coefficient_log(N, N - n); + } else if (N_plus_1 < lgamma_stirling_diff_useful) { return lgamma(N_plus_1) - lgamma(n + 1) - lgamma(N_plus_1 - n); } else { - return_type_t N_minus_n = N - n; - const double one_twelfth = inv(12); - return multiply_log(n, N_minus_n) + multiply_log((N + 0.5), N / N_minus_n) - + one_twelfth / N - n - one_twelfth / N_minus_n - lgamma(n + 1); + return -lbeta(N - n + 1, n + 1) - log(N_plus_1); } } From 77d1e85b85bf5672af3c0fa5e54f9c4aeb096b63 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Tue, 14 Jan 2020 13:27:17 +0100 Subject: [PATCH 03/26] Test for derivatives, docs --- .../prim/fun/binomial_coefficient_log.hpp | 12 +-- .../fun/binomial_coefficient_log_test.cpp | 39 +--------- .../rev/fun/binomial_coefficient_log_test.cpp | 73 +++++++++++++++++++ 3 files changed, 81 insertions(+), 43 deletions(-) create mode 100644 test/unit/math/rev/fun/binomial_coefficient_log_test.cpp diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index f825a229f37..260d2ce25bb 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -1,13 +1,13 @@ #ifndef STAN_MATH_PRIM_SCAL_FUN_BINOMIAL_COEFFICIENT_LOG_HPP #define STAN_MATH_PRIM_SCAL_FUN_BINOMIAL_COEFFICIENT_LOG_HPP -#include #include #include #include #include #include #include +#include namespace stan { namespace math { @@ -28,13 +28,11 @@ namespace math { * = \log \ \Gamma(N+1) - \log \Gamma(n+1) - \log \Gamma(N-n+1)\f$. * * - * TODO[martinmodrak] figure out the cases for x < 0 and for partials \f[ \mbox{binomial\_coefficient\_log}(x, y) = \begin{cases} \textrm{error} & \mbox{if } y > x + 1 \textrm{ or } y < -1 \textrm{ or } x < -1\\ - \textrm{-\infty} & \mbox{if } y = x + 1 \textrm{ or } y = -1\\ \ln\Gamma(x+1) & \mbox{if } -1 < y < x + 1 \\ \quad -\ln\Gamma(y+1)& \\ \quad -\ln\Gamma(x-y+1)& \\[6pt] @@ -45,7 +43,8 @@ namespace math { \f[ \frac{\partial\, \mbox{binomial\_coefficient\_log}(x, y)}{\partial x} = \begin{cases} - \textrm{error} & \mbox{if } y > x \textrm{ or } y < 0\\ + \textrm{error} & \mbox{if } y > x + 1 \textrm{ or } y < -1 \textrm{ or } x + < -1\\ \Psi(x+1) & \mbox{if } 0\leq y \leq x \\ \quad -\Psi(x-y+1)& \\[6pt] \textrm{NaN} & \mbox{if } x = \textrm{NaN or } y = \textrm{NaN} @@ -55,14 +54,15 @@ namespace math { \f[ \frac{\partial\, \mbox{binomial\_coefficient\_log}(x, y)}{\partial y} = \begin{cases} - \textrm{error} & \mbox{if } y > x \textrm{ or } y < 0\\ + \textrm{error} & \mbox{if } y > x + 1 \textrm{ or } y < -1 \textrm{ or } x + < -1\\ -\Psi(y+1) & \mbox{if } 0\leq y \leq x \\ \quad +\Psi(x-y+1)& \\[6pt] \textrm{NaN} & \mbox{if } x = \textrm{NaN or } y = \textrm{NaN} \end{cases} \f] * - * This function is numerically more stable than naive evaluation via lgamma + * This function is numerically more stable than naive evaluation via lgamma. * * @param N total number of objects. * @param n number of objects chosen. diff --git a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp index b4b9fb9b41b..1daec08e457 100644 --- a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp @@ -64,40 +64,5 @@ TEST(MathFunctions, binomial_coefficient_log_errors_edge_cases) { EXPECT_FLOAT_EQ(binomial_coefficient_log(-1, 0), 0); EXPECT_FLOAT_EQ(binomial_coefficient_log(-1, -0.3), inf); EXPECT_FLOAT_EQ(binomial_coefficient_log(0.3, -1), -inf); - -} - -namespace binomial_coefficient_test_internal { -struct TestValue { - double n; - double k; - double val; -}; - -std::vector testValues = { {-0.1,-0.9,-2.1152525390850903}, {-0.1,-1.0000000000000001e-11,1.7771062399418724e-12}, {-0.1,-1.0000000000000002e-6,1.7770950132278696e-7}, {-0.1,-0.001,0.00017592768030333383}, {-0.1,-0.020000000000000004,0.0028416866729211523}, {-0.1,-0.05,0.0044386492587971844}, {-0.1,-0.097,0.0005170837568402953}, {-0.1,-0.09999999400000001,1.0662676433991578e-9}, {-0.1,-0.09999999970000001,5.331331690149795e-11}, {-0.1,0.8,-2.1152525390850907}, {0.00003,-0.9,-2.21375637737528044964112}, {0.00003,3.e-15,1.48040820535563428588846e-19}, {0.00003,3.e-10,1.48039340142161568666212e-14}, {0.00003,3.e-7,1.46560412344434242311062e-11}, {0.00003,6.e-6,2.36865312869367134774661e-10}, {0.00003,0.000015,3.70102051348524044535716e-10}, {0.00003,0.0000291,4.30798787797857754693695e-11}, {0.00003,0.0000299999982,8.88244870007509649977929e-17}, {0.00003,0.00002999999991,4.44122460318735146597007e-18}, {0.00003,0.90003,-2.21375637737528044964112}, {0.002,-0.9,-2.21559326412971099686943}, {0.002,2.e-13,6.57013709556564677684856e-16}, {0.002,2.e-8,6.57007139476544583161173e-11}, {0.002,0.00002,6.50443564072189575550994e-8}, {0.002,0.0004,1.05122171458287350859763e-6}, {0.002,0.001,1.64253373496215313253469e-6}, {0.002,0.00194,1.91190982195918976356429e-7}, {0.002,0.00199999988,3.94208202120835082737684e-13}, {0.002,0.001999999994,1.97104112295366725515452e-14}, {0.002,0.902,-2.21559326412971099686943}, {1,-0.9,-2.85558226198351740582195}, {1,1.e-10,9.9999999988550659331851e-11}, {1,0.00001,9.99988550692664559909352e-6}, {1,0.01,0.00988583703486131052627978}, {1,0.2,0.156457962917688016707705}, {1,0.5,0.241564475270490444691037}, {1,0.97,0.028978328236256312960776}, {1,0.99999994,5.99999958782374313463811e-8}, {1,0.999999997,2.99999998969559340736596e-9}, {1,1.9,-2.85558226198351740582195}, {8,-0.9,-4.22528965320883461943031}, {8,8.e-10,2.17428571372173153982474e-9}, {8,0.00008,0.000217422931805073420417006}, {8,0.08,0.211982267378255838975509}, {8,1.6,2.90678606291134283426263}, {8,4,4.24849524204935898912334}, {8,7.76,0.606274586245453651115361}, {8,7.99999952,1.30457122553768403613331e-6}, {8,7.999999976,6.52285709209869625945566e-8}, {8,8.9,-4.22528965320883461943031}, {1325,-0.9,-8.72360867216657209762532}, {1325,1.325e-7,1.02909578020477960435539e-6}, {1325,0.01325,0.102766042691370430370992}, {1325,13.25,71.9898321274090629975055}, {1325,265,659.435649329029419323398}, {1325,662.5,914.599450340845275100724}, {1325,1285.25,175.786260651191862665015}, {1325,1324.9999205,0.000617452276410452190170437}, {1325,1324.999996025,0.0000308728608380968862741097}, {1325,1325.9,-8.72360867216657209762532}, {845000,-0.9,-14.5350963792733464918229}, {845000,0.0000845,0.00120194816738712136581358}, {845000,8.45,103.738303827736743600251}, {845000,8450,47315.8616457576200611209}, {845000,169000,422833.221695496506553128}, {845000,422500,585702.318235552114086514}, {845000,819650,113851.158132678562120685}, {845000,844999.9493,0.719108776819481762797449}, {845000,844999.997465,0.036053342347290003917417}, {845000,845000.9,-14.5350963792733464918229}, {3000000000000000.0,3,105.120406581508328854183}, {3000000000000000.0,100,3199.99949280231435502243}, {3000000000000000.0,12895,350387.5243605883687667}, {100000000000000000000.0,3,136.363346110414686040237}, {100000000000000000000.0,100,4241.4308104325278778424}, {100000000000000000000.0,12895,484680.092769031900296878},}; -} // namespace binomial_coefficient_test_internal - -TEST(MathFunctions, binomial_coefficient_log_precomputed) { - using binomial_coefficient_test_internal::TestValue; - using binomial_coefficient_test_internal::testValues; - using stan::test::expect_near_rel; - - for (TestValue t : testValues) { - std::ostringstream msg; - msg << std::setprecision(22) << "n = " << t.n << ", k = " << t.k; - - double val = stan::math::binomial_coefficient_log(t.n, t.k); - expect_near_rel(msg.str(), val, t.val); - - if (t.k > t.n / 10) { - // Testing the mirrored binomial coefficient is not performed for - // small k, as we may loose so much precision computing k2 - // that the test becomes invalid - std::ostringstream msg2; - double k2 = t.n - t.k; - msg2 << std::setprecision(22) << "n = " << t.n << ", k = " << k2; - double val2 = stan::math::binomial_coefficient_log(t.n, k2); - expect_near_rel(msg2.str(), val2, t.val); - } - } -} + EXPECT_FLOAT_EQ(binomial_coefficient_log(5.0, 6.0), -inf); +} \ No newline at end of file diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp new file mode 100644 index 00000000000..1aee5b9dae3 --- /dev/null +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -0,0 +1,73 @@ +#include +#include +#include +#include +#include +#include +#include + +namespace binomial_coefficient_log_test_internal { +struct TestValue { + double n; + double k; + double val; + double dn; + double dk; +}; + +constexpr double NaN = std::numeric_limits::quiet_NaN(); +// Test values generated in Mathematica, reproducible notebook at +// https://www.wolframcloud.com/obj/martin.modrak/Published/binomial_coefficient_log.nb +// Mathematica Code reproduced below for convenience: +// TODO + + +std::vector testValues = { {-0.1,-0.9,-2.1152525390850903,-1.039918383240913,10.708746373704937}, {-0.1,-1.0000000000000001e-11,1.7771062399418724e-12,-1.922551007282891e-11,-0.17771128500984368}, {-0.1,-1.0000000000000002e-6,1.7770950132278696e-7,-1.9225383585119715e-6,-0.1777077175718912}, {-0.1,-0.001,0.00017592768030333383,-0.0019209406823855746,-0.17414420715602041}, {-0.1,-0.020000000000000004,0.0028416866729211523,-0.037823152614322564,-0.10649980051731356}, {-0.1,-0.05,0.0044386492587971844,-0.09231733374172624,0.}, {-0.1,-0.097,0.0005170837568402953,-0.17276563502587905,0.16701237954910164}, {-0.1,-0.09999999400000001,1.0662676433991578e-9,-0.17771127517591423,0.17771126364067447}, {-0.1,-0.09999999970000001,5.331331690149795e-11,-0.17771128455203855,0.17771128397527658}, {-0.1,0.8,-2.1152525390850907,9.668827990464031,-10.708746373704944}, {0.00003,-0.9,-2.21375637737528044964112,-0.933371118080918307851728,10.7799597405306456982503}, {0.00003,3.e-15,1.48040820535563428588846e-19,4.9345858390686036409092e-15,0.0000493469401735864488431421}, {0.00003,3.e-10,1.48039340142161568666212e-14,4.93458584015035637326682e-10,0.0000493459532446518755643083}, {0.00003,3.e-7,1.46560412344434242311062e-11,4.93458692083243471796274e-7,0.0000483600013795032175416619}, {0.00003,6.e-6,2.36865312869367134774661e-10,9.86921494891292997054538e-6,0.0000296081641072682815641386}, {0.00003,0.000015,3.70102051348524044535716e-10,0.0000246731996398828628243033,0}, {0.00003,0.0000291,4.30798787797857754693695e-11,0.0000478665004969577356362445,-0.0000463861237716491741348152}, {0.00003,0.0000299999982,8.88244870007509649977929e-17,0.0000493469372225745165598424,-0.0000493469342618230131147924}, {0.00003,0.00002999999991,4.44122460318735146597007e-18,0.0000493469400354117708664247,-0.0000493469398873741956943572}, {0.00003,0.90003,-2.21375637737528044964112,9.84658862244972739039859,-10.7799597405306456982503}, {0.002,-0.9,-2.21559326412971099686943,-0.931489785799666354340725,10.7813141298573216195769}, {0.002,2.e-13,6.57013709556564677684856e-16,3.2802775880282756733109e-13,0.00328506854745431610233807}, {0.002,2.e-8,6.57007139476544583161173e-11,3.28027763585130905255152e-8,0.00328500284665411601964384}, {0.002,0.00002,6.50443564072189575550994e-8,0.0000328032541239782363672409,0.00321936709346449075328701}, {0.002,0.0004,1.05122171458287350859763e-6,0.000656246880415250668008776,0.00197104030081914294655284}, {0.002,0.001,1.64253373496215313253469e-6,0.00164133545687894154061287,0}, {0.002,0.00194,1.91190982195918976356429e-7,0.00318637683127152008217264,-0.00308796419928270612987841}, {0.002,0.00199999988,3.94208202120835082737684e-13,0.00328506835071924267003294,-0.00328506815390258737114863}, {0.002,0.001999999994,1.97104112295366725515452e-14,0.00328506853824172627346441,-0.00328506852840089350933798}, {0.002,0.902,-2.21559326412971099686943,9.84982434405765526523621,-10.7813141298573216195769}, {1,-0.9,-2.85558226198351740582195,-0.459715615539276790357555,11.3062548910488207249193}, {1,1.e-10,9.9999999988550659331851e-11,6.44934066868432126789198e-11,0.99999999977101318664035}, {1,0.00001,9.99988550692664559909352e-6,6.44936087425490392714312e-6,0.999977101418661870834803}, {1,0.01,0.00988583703486131052627978,0.00646962905305218281638255,0.977200163914089454058067}, {1,0.2,0.156457962917688016707705,0.137792901804605598784786,0.57403132988604983615591}, {1,0.5,0.241564475270490444691037,0.386294361119890618834464,0}, {1,0.97,0.028978328236256312960776,0.951705422383897641079507,-0.932173296099201807475084}, {1,0.99999994,5.99999958782374313463811e-8,0.999999901303960316511031,-0.999999862607915578212576}, {1,0.999999997,2.99999998969559340736596e-9,0.999999995065197810273833,-0.999999993130395607910641}, {1,1.9,-2.85558226198351740582195,10.8465392755095439345617,-11.3062548910488207249193}, {8,-0.9,-4.22528965320883461943031,-0.100538838650771402252215,12.6649352570174581939568}, {8,8.e-10,2.17428571372173153982474e-9,9.4009611759639002250559e-11,2.71785714144718599267395}, {8,0.00008,0.000217422931805073420417006,9.4010053144199442029858e-6,2.71771615481909065227871}, {8,0.08,0.211982267378255838975509,0.00944537775278953163286055,2.58399543820453334937797}, {8,1.6,2.90678606291134283426263,0.208248082071609202561846,1.1813459431105245420685}, {8,4,4.24849524204935898912334,0.634523809523809523809524,0}, {8,7.76,0.606274586245453651115361,2.38013515708154679097544,-2.35152741320850411987141}, {8,7.99999952,1.30457122553768403613331e-6,2.71785635328906772378496,-2.71785629688329908165944}, {8,7.999999976,6.52285709209869625945566e-8,2.71785710337872594517018,-2.71785710055843758854095}, {8,8.9,-4.22528965320883461943031,12.5643964183666867917046,-12.6649352570174581939568}, {1325,-0.9,-8.72360867216657209762532,-0.000678758619778166713080174,17.6139787484752694858826}, {1325,1.325e-7,1.02909578020477960435539e-6,9.99622736492339027462876e-11,7.76676049629224235717761}, {1325,0.01325,0.102766042691370430370992,9.99627732703511703537678e-6,7.74516389164454804001553}, {1325,13.25,71.9898321274090629975055,0.0100465251153381679080642,4.55823951796927536337622}, {1325,265,659.435649329029419323398,0.223049238391457549660772,1.38488037927413480754602}, {1325,662.5,914.599450340845275100724,0.692769964468769118810416,0}, {1325,1285.25,175.786260651191862665015,3.49440932918478008636475,-3.46396178959974953213922}, {1325,1324.9999205,0.000617452276410452190170437,7.76662994968438931022129,-7.76662988970702332503659}, {1325,1324.999996025,0.0000308728608380968862741097,7.76675417575202484528552,-7.76675417275315663146178}, {1325,1325.9,-8.72360867216657209762532,17.6132999898554913191696,-17.6139787484752694858826}, {845000,-0.9,-14.5350963792733464918229,-1.06508755996070778708647e-6,24.0708485035538072767331}, {845000,0.0000845,0.00120194816738712136581358,9.99999408334257028001921e-11,14.2241691745103886002024}, {845000,8.45,103.738303827736743600251,0.0000100000440831167345507131,11.454909927882486225397}, {845000,8450,47315.8616457576200611209,0.0100503298765747570022521,4.59506127739683575088585}, {845000,169000,422833.221695496506553128,0.223143403385281321887929,1.3862921421877147496422}, {845000,422500,585702.318235552114086514,0.693146588844319105852931,0}, {845000,819650,113851.158132678562120685,3.50653876529964054840506,-3.476079576115418785061}, {845000,844999.9493,0.719108776819481762797449,14.1438653554064391526676,-14.1438652954064728556142}, {845000,844999.997465,0.036053342347290003917417,14.2201459621965153509889,-14.2201459591965171216362}, {845000,845000.9,-14.5350963792733464918229,24.0708474384662473160253,-24.0708485035538072767331}, {3.0000000000000005e15,3,105.120406581508329354183,NaN,NaN}, {3.0000000000000005e15,100,3199.9994928023143716891,NaN,NaN}, {3.0000000000000005e15,12895,350387.524360588370915867,NaN,NaN}, {1.e20,3,136.363346110414686040252,NaN,NaN}, {1.e20,100,4241.4308104325278778429,NaN,NaN}, {1.e20,12895,484680.092769031900296942,NaN,NaN},}; + +} // namespace binomial_coefficient_log_test_internal + +TEST(MathFunctions, binomial_coefficient_log_precomputed) { + using binomial_coefficient_log_test_internal::TestValue; + using binomial_coefficient_log_test_internal::testValues; + using stan::math::is_nan; + using stan::math::value_of; + using stan::math::var; + using stan::test::expect_near_rel; + + for (TestValue t : testValues) { + std::ostringstream msg; + msg << std::setprecision(22) << "n = " << t.n << ", k = " << t.k; + + var n(t.n); + var k(t.k); + var val = stan::math::binomial_coefficient_log(n, k); + + std::vector vars; + vars.push_back(n); + vars.push_back(k); + + std::vector gradients; + val.grad(vars, gradients); + + for (int i = 0; i < 2; ++i) { + EXPECT_FALSE(is_nan(gradients[i])); + } + + expect_near_rel(msg.str(), value_of(val), t.val); + + double tol_grad; + if (n < 1 || k < 1) { + tol_grad = 1e-7; + } else { + tol_grad = 1e-8; + } + if (!is_nan(t.dn)) { + expect_near_rel(std::string("dn: ") + msg.str(), gradients[0], t.dn, + tol_grad); + } + if (!is_nan(t.dk)) { + expect_near_rel(std::string("dk: ") + msg.str(), gradients[1], t.dk, + tol_grad); + } + } +} From 118adb2b17f5cb8b62bcf76af8675e9cd93b5f82 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 14 Jan 2020 12:28:12 +0000 Subject: [PATCH 04/26] [Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (tags/RELEASE_500/final) --- .../rev/fun/binomial_coefficient_log_test.cpp | 143 +++++++++++++++++- 1 file changed, 141 insertions(+), 2 deletions(-) diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 1aee5b9dae3..bae569e01b4 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -21,8 +21,147 @@ constexpr double NaN = std::numeric_limits::quiet_NaN(); // Mathematica Code reproduced below for convenience: // TODO - -std::vector testValues = { {-0.1,-0.9,-2.1152525390850903,-1.039918383240913,10.708746373704937}, {-0.1,-1.0000000000000001e-11,1.7771062399418724e-12,-1.922551007282891e-11,-0.17771128500984368}, {-0.1,-1.0000000000000002e-6,1.7770950132278696e-7,-1.9225383585119715e-6,-0.1777077175718912}, {-0.1,-0.001,0.00017592768030333383,-0.0019209406823855746,-0.17414420715602041}, {-0.1,-0.020000000000000004,0.0028416866729211523,-0.037823152614322564,-0.10649980051731356}, {-0.1,-0.05,0.0044386492587971844,-0.09231733374172624,0.}, {-0.1,-0.097,0.0005170837568402953,-0.17276563502587905,0.16701237954910164}, {-0.1,-0.09999999400000001,1.0662676433991578e-9,-0.17771127517591423,0.17771126364067447}, {-0.1,-0.09999999970000001,5.331331690149795e-11,-0.17771128455203855,0.17771128397527658}, {-0.1,0.8,-2.1152525390850907,9.668827990464031,-10.708746373704944}, {0.00003,-0.9,-2.21375637737528044964112,-0.933371118080918307851728,10.7799597405306456982503}, {0.00003,3.e-15,1.48040820535563428588846e-19,4.9345858390686036409092e-15,0.0000493469401735864488431421}, {0.00003,3.e-10,1.48039340142161568666212e-14,4.93458584015035637326682e-10,0.0000493459532446518755643083}, {0.00003,3.e-7,1.46560412344434242311062e-11,4.93458692083243471796274e-7,0.0000483600013795032175416619}, {0.00003,6.e-6,2.36865312869367134774661e-10,9.86921494891292997054538e-6,0.0000296081641072682815641386}, {0.00003,0.000015,3.70102051348524044535716e-10,0.0000246731996398828628243033,0}, {0.00003,0.0000291,4.30798787797857754693695e-11,0.0000478665004969577356362445,-0.0000463861237716491741348152}, {0.00003,0.0000299999982,8.88244870007509649977929e-17,0.0000493469372225745165598424,-0.0000493469342618230131147924}, {0.00003,0.00002999999991,4.44122460318735146597007e-18,0.0000493469400354117708664247,-0.0000493469398873741956943572}, {0.00003,0.90003,-2.21375637737528044964112,9.84658862244972739039859,-10.7799597405306456982503}, {0.002,-0.9,-2.21559326412971099686943,-0.931489785799666354340725,10.7813141298573216195769}, {0.002,2.e-13,6.57013709556564677684856e-16,3.2802775880282756733109e-13,0.00328506854745431610233807}, {0.002,2.e-8,6.57007139476544583161173e-11,3.28027763585130905255152e-8,0.00328500284665411601964384}, {0.002,0.00002,6.50443564072189575550994e-8,0.0000328032541239782363672409,0.00321936709346449075328701}, {0.002,0.0004,1.05122171458287350859763e-6,0.000656246880415250668008776,0.00197104030081914294655284}, {0.002,0.001,1.64253373496215313253469e-6,0.00164133545687894154061287,0}, {0.002,0.00194,1.91190982195918976356429e-7,0.00318637683127152008217264,-0.00308796419928270612987841}, {0.002,0.00199999988,3.94208202120835082737684e-13,0.00328506835071924267003294,-0.00328506815390258737114863}, {0.002,0.001999999994,1.97104112295366725515452e-14,0.00328506853824172627346441,-0.00328506852840089350933798}, {0.002,0.902,-2.21559326412971099686943,9.84982434405765526523621,-10.7813141298573216195769}, {1,-0.9,-2.85558226198351740582195,-0.459715615539276790357555,11.3062548910488207249193}, {1,1.e-10,9.9999999988550659331851e-11,6.44934066868432126789198e-11,0.99999999977101318664035}, {1,0.00001,9.99988550692664559909352e-6,6.44936087425490392714312e-6,0.999977101418661870834803}, {1,0.01,0.00988583703486131052627978,0.00646962905305218281638255,0.977200163914089454058067}, {1,0.2,0.156457962917688016707705,0.137792901804605598784786,0.57403132988604983615591}, {1,0.5,0.241564475270490444691037,0.386294361119890618834464,0}, {1,0.97,0.028978328236256312960776,0.951705422383897641079507,-0.932173296099201807475084}, {1,0.99999994,5.99999958782374313463811e-8,0.999999901303960316511031,-0.999999862607915578212576}, {1,0.999999997,2.99999998969559340736596e-9,0.999999995065197810273833,-0.999999993130395607910641}, {1,1.9,-2.85558226198351740582195,10.8465392755095439345617,-11.3062548910488207249193}, {8,-0.9,-4.22528965320883461943031,-0.100538838650771402252215,12.6649352570174581939568}, {8,8.e-10,2.17428571372173153982474e-9,9.4009611759639002250559e-11,2.71785714144718599267395}, {8,0.00008,0.000217422931805073420417006,9.4010053144199442029858e-6,2.71771615481909065227871}, {8,0.08,0.211982267378255838975509,0.00944537775278953163286055,2.58399543820453334937797}, {8,1.6,2.90678606291134283426263,0.208248082071609202561846,1.1813459431105245420685}, {8,4,4.24849524204935898912334,0.634523809523809523809524,0}, {8,7.76,0.606274586245453651115361,2.38013515708154679097544,-2.35152741320850411987141}, {8,7.99999952,1.30457122553768403613331e-6,2.71785635328906772378496,-2.71785629688329908165944}, {8,7.999999976,6.52285709209869625945566e-8,2.71785710337872594517018,-2.71785710055843758854095}, {8,8.9,-4.22528965320883461943031,12.5643964183666867917046,-12.6649352570174581939568}, {1325,-0.9,-8.72360867216657209762532,-0.000678758619778166713080174,17.6139787484752694858826}, {1325,1.325e-7,1.02909578020477960435539e-6,9.99622736492339027462876e-11,7.76676049629224235717761}, {1325,0.01325,0.102766042691370430370992,9.99627732703511703537678e-6,7.74516389164454804001553}, {1325,13.25,71.9898321274090629975055,0.0100465251153381679080642,4.55823951796927536337622}, {1325,265,659.435649329029419323398,0.223049238391457549660772,1.38488037927413480754602}, {1325,662.5,914.599450340845275100724,0.692769964468769118810416,0}, {1325,1285.25,175.786260651191862665015,3.49440932918478008636475,-3.46396178959974953213922}, {1325,1324.9999205,0.000617452276410452190170437,7.76662994968438931022129,-7.76662988970702332503659}, {1325,1324.999996025,0.0000308728608380968862741097,7.76675417575202484528552,-7.76675417275315663146178}, {1325,1325.9,-8.72360867216657209762532,17.6132999898554913191696,-17.6139787484752694858826}, {845000,-0.9,-14.5350963792733464918229,-1.06508755996070778708647e-6,24.0708485035538072767331}, {845000,0.0000845,0.00120194816738712136581358,9.99999408334257028001921e-11,14.2241691745103886002024}, {845000,8.45,103.738303827736743600251,0.0000100000440831167345507131,11.454909927882486225397}, {845000,8450,47315.8616457576200611209,0.0100503298765747570022521,4.59506127739683575088585}, {845000,169000,422833.221695496506553128,0.223143403385281321887929,1.3862921421877147496422}, {845000,422500,585702.318235552114086514,0.693146588844319105852931,0}, {845000,819650,113851.158132678562120685,3.50653876529964054840506,-3.476079576115418785061}, {845000,844999.9493,0.719108776819481762797449,14.1438653554064391526676,-14.1438652954064728556142}, {845000,844999.997465,0.036053342347290003917417,14.2201459621965153509889,-14.2201459591965171216362}, {845000,845000.9,-14.5350963792733464918229,24.0708474384662473160253,-24.0708485035538072767331}, {3.0000000000000005e15,3,105.120406581508329354183,NaN,NaN}, {3.0000000000000005e15,100,3199.9994928023143716891,NaN,NaN}, {3.0000000000000005e15,12895,350387.524360588370915867,NaN,NaN}, {1.e20,3,136.363346110414686040252,NaN,NaN}, {1.e20,100,4241.4308104325278778429,NaN,NaN}, {1.e20,12895,484680.092769031900296942,NaN,NaN},}; +std::vector testValues = { + {-0.1, -0.9, -2.1152525390850903, -1.039918383240913, 10.708746373704937}, + {-0.1, -1.0000000000000001e-11, 1.7771062399418724e-12, + -1.922551007282891e-11, -0.17771128500984368}, + {-0.1, -1.0000000000000002e-6, 1.7770950132278696e-7, + -1.9225383585119715e-6, -0.1777077175718912}, + {-0.1, -0.001, 0.00017592768030333383, -0.0019209406823855746, + -0.17414420715602041}, + {-0.1, -0.020000000000000004, 0.0028416866729211523, -0.037823152614322564, + -0.10649980051731356}, + {-0.1, -0.05, 0.0044386492587971844, -0.09231733374172624, 0.}, + {-0.1, -0.097, 0.0005170837568402953, -0.17276563502587905, + 0.16701237954910164}, + {-0.1, -0.09999999400000001, 1.0662676433991578e-9, -0.17771127517591423, + 0.17771126364067447}, + {-0.1, -0.09999999970000001, 5.331331690149795e-11, -0.17771128455203855, + 0.17771128397527658}, + {-0.1, 0.8, -2.1152525390850907, 9.668827990464031, -10.708746373704944}, + {0.00003, -0.9, -2.21375637737528044964112, -0.933371118080918307851728, + 10.7799597405306456982503}, + {0.00003, 3.e-15, 1.48040820535563428588846e-19, + 4.9345858390686036409092e-15, 0.0000493469401735864488431421}, + {0.00003, 3.e-10, 1.48039340142161568666212e-14, + 4.93458584015035637326682e-10, 0.0000493459532446518755643083}, + {0.00003, 3.e-7, 1.46560412344434242311062e-11, + 4.93458692083243471796274e-7, 0.0000483600013795032175416619}, + {0.00003, 6.e-6, 2.36865312869367134774661e-10, + 9.86921494891292997054538e-6, 0.0000296081641072682815641386}, + {0.00003, 0.000015, 3.70102051348524044535716e-10, + 0.0000246731996398828628243033, 0}, + {0.00003, 0.0000291, 4.30798787797857754693695e-11, + 0.0000478665004969577356362445, -0.0000463861237716491741348152}, + {0.00003, 0.0000299999982, 8.88244870007509649977929e-17, + 0.0000493469372225745165598424, -0.0000493469342618230131147924}, + {0.00003, 0.00002999999991, 4.44122460318735146597007e-18, + 0.0000493469400354117708664247, -0.0000493469398873741956943572}, + {0.00003, 0.90003, -2.21375637737528044964112, 9.84658862244972739039859, + -10.7799597405306456982503}, + {0.002, -0.9, -2.21559326412971099686943, -0.931489785799666354340725, + 10.7813141298573216195769}, + {0.002, 2.e-13, 6.57013709556564677684856e-16, 3.2802775880282756733109e-13, + 0.00328506854745431610233807}, + {0.002, 2.e-8, 6.57007139476544583161173e-11, 3.28027763585130905255152e-8, + 0.00328500284665411601964384}, + {0.002, 0.00002, 6.50443564072189575550994e-8, + 0.0000328032541239782363672409, 0.00321936709346449075328701}, + {0.002, 0.0004, 1.05122171458287350859763e-6, 0.000656246880415250668008776, + 0.00197104030081914294655284}, + {0.002, 0.001, 1.64253373496215313253469e-6, 0.00164133545687894154061287, + 0}, + {0.002, 0.00194, 1.91190982195918976356429e-7, 0.00318637683127152008217264, + -0.00308796419928270612987841}, + {0.002, 0.00199999988, 3.94208202120835082737684e-13, + 0.00328506835071924267003294, -0.00328506815390258737114863}, + {0.002, 0.001999999994, 1.97104112295366725515452e-14, + 0.00328506853824172627346441, -0.00328506852840089350933798}, + {0.002, 0.902, -2.21559326412971099686943, 9.84982434405765526523621, + -10.7813141298573216195769}, + {1, -0.9, -2.85558226198351740582195, -0.459715615539276790357555, + 11.3062548910488207249193}, + {1, 1.e-10, 9.9999999988550659331851e-11, 6.44934066868432126789198e-11, + 0.99999999977101318664035}, + {1, 0.00001, 9.99988550692664559909352e-6, 6.44936087425490392714312e-6, + 0.999977101418661870834803}, + {1, 0.01, 0.00988583703486131052627978, 0.00646962905305218281638255, + 0.977200163914089454058067}, + {1, 0.2, 0.156457962917688016707705, 0.137792901804605598784786, + 0.57403132988604983615591}, + {1, 0.5, 0.241564475270490444691037, 0.386294361119890618834464, 0}, + {1, 0.97, 0.028978328236256312960776, 0.951705422383897641079507, + -0.932173296099201807475084}, + {1, 0.99999994, 5.99999958782374313463811e-8, 0.999999901303960316511031, + -0.999999862607915578212576}, + {1, 0.999999997, 2.99999998969559340736596e-9, 0.999999995065197810273833, + -0.999999993130395607910641}, + {1, 1.9, -2.85558226198351740582195, 10.8465392755095439345617, + -11.3062548910488207249193}, + {8, -0.9, -4.22528965320883461943031, -0.100538838650771402252215, + 12.6649352570174581939568}, + {8, 8.e-10, 2.17428571372173153982474e-9, 9.4009611759639002250559e-11, + 2.71785714144718599267395}, + {8, 0.00008, 0.000217422931805073420417006, 9.4010053144199442029858e-6, + 2.71771615481909065227871}, + {8, 0.08, 0.211982267378255838975509, 0.00944537775278953163286055, + 2.58399543820453334937797}, + {8, 1.6, 2.90678606291134283426263, 0.208248082071609202561846, + 1.1813459431105245420685}, + {8, 4, 4.24849524204935898912334, 0.634523809523809523809524, 0}, + {8, 7.76, 0.606274586245453651115361, 2.38013515708154679097544, + -2.35152741320850411987141}, + {8, 7.99999952, 1.30457122553768403613331e-6, 2.71785635328906772378496, + -2.71785629688329908165944}, + {8, 7.999999976, 6.52285709209869625945566e-8, 2.71785710337872594517018, + -2.71785710055843758854095}, + {8, 8.9, -4.22528965320883461943031, 12.5643964183666867917046, + -12.6649352570174581939568}, + {1325, -0.9, -8.72360867216657209762532, -0.000678758619778166713080174, + 17.6139787484752694858826}, + {1325, 1.325e-7, 1.02909578020477960435539e-6, + 9.99622736492339027462876e-11, 7.76676049629224235717761}, + {1325, 0.01325, 0.102766042691370430370992, 9.99627732703511703537678e-6, + 7.74516389164454804001553}, + {1325, 13.25, 71.9898321274090629975055, 0.0100465251153381679080642, + 4.55823951796927536337622}, + {1325, 265, 659.435649329029419323398, 0.223049238391457549660772, + 1.38488037927413480754602}, + {1325, 662.5, 914.599450340845275100724, 0.692769964468769118810416, 0}, + {1325, 1285.25, 175.786260651191862665015, 3.49440932918478008636475, + -3.46396178959974953213922}, + {1325, 1324.9999205, 0.000617452276410452190170437, + 7.76662994968438931022129, -7.76662988970702332503659}, + {1325, 1324.999996025, 0.0000308728608380968862741097, + 7.76675417575202484528552, -7.76675417275315663146178}, + {1325, 1325.9, -8.72360867216657209762532, 17.6132999898554913191696, + -17.6139787484752694858826}, + {845000, -0.9, -14.5350963792733464918229, -1.06508755996070778708647e-6, + 24.0708485035538072767331}, + {845000, 0.0000845, 0.00120194816738712136581358, + 9.99999408334257028001921e-11, 14.2241691745103886002024}, + {845000, 8.45, 103.738303827736743600251, 0.0000100000440831167345507131, + 11.454909927882486225397}, + {845000, 8450, 47315.8616457576200611209, 0.0100503298765747570022521, + 4.59506127739683575088585}, + {845000, 169000, 422833.221695496506553128, 0.223143403385281321887929, + 1.3862921421877147496422}, + {845000, 422500, 585702.318235552114086514, 0.693146588844319105852931, 0}, + {845000, 819650, 113851.158132678562120685, 3.50653876529964054840506, + -3.476079576115418785061}, + {845000, 844999.9493, 0.719108776819481762797449, 14.1438653554064391526676, + -14.1438652954064728556142}, + {845000, 844999.997465, 0.036053342347290003917417, + 14.2201459621965153509889, -14.2201459591965171216362}, + {845000, 845000.9, -14.5350963792733464918229, 24.0708474384662473160253, + -24.0708485035538072767331}, + {3.0000000000000005e15, 3, 105.120406581508329354183, NaN, NaN}, + {3.0000000000000005e15, 100, 3199.9994928023143716891, NaN, NaN}, + {3.0000000000000005e15, 12895, 350387.524360588370915867, NaN, NaN}, + {1.e20, 3, 136.363346110414686040252, NaN, NaN}, + {1.e20, 100, 4241.4308104325278778429, NaN, NaN}, + {1.e20, 12895, 484680.092769031900296942, NaN, NaN}, +}; } // namespace binomial_coefficient_log_test_internal From 3578ff611f0b82afd07d3d5ce83bcc490e2d8a22 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Tue, 14 Jan 2020 13:37:06 +0100 Subject: [PATCH 05/26] Lint, test sources --- .../fun/binomial_coefficient_log_test.cpp | 2 +- .../rev/fun/binomial_coefficient_log_test.cpp | 39 ++++++++++++++++++- 2 files changed, 39 insertions(+), 2 deletions(-) diff --git a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp index 1daec08e457..4b973091b5e 100644 --- a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp @@ -65,4 +65,4 @@ TEST(MathFunctions, binomial_coefficient_log_errors_edge_cases) { EXPECT_FLOAT_EQ(binomial_coefficient_log(-1, -0.3), inf); EXPECT_FLOAT_EQ(binomial_coefficient_log(0.3, -1), -inf); EXPECT_FLOAT_EQ(binomial_coefficient_log(5.0, 6.0), -inf); -} \ No newline at end of file +} diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index bae569e01b4..9b02298ed74 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -19,8 +19,45 @@ constexpr double NaN = std::numeric_limits::quiet_NaN(); // Test values generated in Mathematica, reproducible notebook at // https://www.wolframcloud.com/obj/martin.modrak/Published/binomial_coefficient_log.nb // Mathematica Code reproduced below for convenience: -// TODO +// bclog[n_,k_]:= LogGamma[n + 1] - LogGamma[k + 1] - LogGamma[n + 1 -k ] +// dbclogdn[n_,k_]= D[bclog[n,k],n]; +// dbclogdk[n_,k_]= D[bclog[n,k],k]; + +// singleTest[x_,y_]:= Module[{val, cdn,cdk},{ +// val = N[bclog[x,y],24]; +// cdn = If[x > 10^6 || y > 10^6,"NaN", CForm[N[dbclogdn[x,y],24]]]; +// cdk = If[x > 10^6 || y > 10^6,"NaN", CForm[N[dbclogdk[x,y],24]]]; +// WriteString[out," {",CForm[x],",",CForm[y],",", +// CForm[val],",",cdn,",",cdk,"},"] +// }]; + +// out = OpenWrite["bc_test.txt"]; +// ns= {-0.1,3*10^-5,2*10^-3,1,8, 1325,845*10^3}; +// ratios = {-1,10^- 10,10^-5,10^-2,1/5,1/2,1-3*10^-2,1-6*10^-8, 1 -3*10^-9,2}; +// WriteString[out, "std::vector testValues = {"]; +// For[i = 1, i <= Length[ns], i++, { +// For[j = 1, j <= Length[ratios], j++, { +// cn = ns[[i]]; +// ck = If[ratios[[j]] < 0,-9/10, +// If[ratios[[j]] > 1,cn + 9/10,cn * ratios[[j]] ]]; +// singleTest[cn,ck]; +// }] +// }] + +// extremeNs = {3*10^15+1/2,10^20 + 1/2}; +// lowKs = {3, 100, 12895}; +// For[i = 1, i <= Length[extremeNs], i++, { +// For[j = 1, j <= Length[lowKs], j++, { +// cn = extremeNs[[i]]; +// ck = lowKs[[j]]; +// singleTest[cn,ck]; +// }] +// }] + +// WriteString[out,"};"]; +// Close[out]; +// FilePrint[%] std::vector testValues = { {-0.1, -0.9, -2.1152525390850903, -1.039918383240913, 10.708746373704937}, {-0.1, -1.0000000000000001e-11, 1.7771062399418724e-12, From e03205d853bfce17ecbf40e28fecde1d98faeed3 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Tue, 14 Jan 2020 15:09:05 +0100 Subject: [PATCH 06/26] Mismatched define --- stan/math/prim/fun/binomial_coefficient_log.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 260d2ce25bb..1a112518577 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -1,5 +1,5 @@ -#ifndef STAN_MATH_PRIM_SCAL_FUN_BINOMIAL_COEFFICIENT_LOG_HPP -#define STAN_MATH_PRIM_SCAL_FUN_BINOMIAL_COEFFICIENT_LOG_HPP +#ifndef STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP +#define STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP #include #include From 035212df693a40b8b49be45c68d2712024a9fe60 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Tue, 14 Jan 2020 16:09:59 +0100 Subject: [PATCH 07/26] Improved code style --- .../prim/fun/binomial_coefficient_log.hpp | 33 +++++++++---------- .../fun/binomial_coefficient_log_test.cpp | 14 ++++---- .../rev/fun/binomial_coefficient_log_test.cpp | 7 ++-- 3 files changed, 26 insertions(+), 28 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 1a112518577..cd39e80b0fe 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -2,12 +2,11 @@ #define STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP #include -#include -#include +#include #include -#include -#include -#include +#include +#include +#include namespace stan { namespace math { @@ -64,6 +63,8 @@ namespace math { * * This function is numerically more stable than naive evaluation via lgamma. * + * @tparam T_N type of N. + * @tparam T_n type of n. * @param N total number of objects. * @param n number of objects chosen. * @return log (N choose n). @@ -72,30 +73,28 @@ namespace math { template inline return_type_t binomial_coefficient_log(const T_N N, const T_n n) { - if (is_nan(value_of_rec(N)) || is_nan(value_of_rec(n))) { - return std::numeric_limits::quiet_NaN(); + if (is_any_nan(N, n)) { + return stan::math::NOT_A_NUMBER; } - // For some uses it is important this works even when N < 0 and therefore - // it is before checks - if (n == 0) { - return 0; - } const T_N N_plus_1 = N + 1; static const char* function = "binomial_coefficient_log"; check_greater_or_equal(function, "first argument", N, -1); check_greater_or_equal(function, "second argument", n, -1); check_greater_or_equal(function, "(first argument - second argument + 1)", - N - n + 1, 0.0); + N_plus_1 - n, 0.0); + if (n == 0) { + return 0; + } if (N / 2 < n) { return binomial_coefficient_log(N, N - n); - } else if (N_plus_1 < lgamma_stirling_diff_useful) { + } + if (N_plus_1 < lgamma_stirling_diff_useful) { return lgamma(N_plus_1) - lgamma(n + 1) - lgamma(N_plus_1 - n); - } else { - return -lbeta(N - n + 1, n + 1) - log(N_plus_1); - } + } + return -lbeta(N - n + 1, n + 1) - log(N_plus_1); } } // namespace math diff --git a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp index 4b973091b5e..e372b6815d9 100644 --- a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp @@ -1,8 +1,7 @@ -#include #include +#include #include #include -#include template void test_binom_coefficient(const T_N& N, const T_n& n) { @@ -52,7 +51,8 @@ TEST(MathFunctions, binomial_coefficient_log_nan) { TEST(MathFunctions, binomial_coefficient_log_errors_edge_cases) { using stan::math::binomial_coefficient_log; - double inf = std::numeric_limits::infinity(); + using stan::math::INFTY; + EXPECT_NO_THROW(binomial_coefficient_log(10, 11)); EXPECT_THROW(binomial_coefficient_log(10, 11.01), std::domain_error); EXPECT_THROW(binomial_coefficient_log(10, -1.1), std::domain_error); @@ -60,9 +60,9 @@ TEST(MathFunctions, binomial_coefficient_log_errors_edge_cases) { EXPECT_NO_THROW(binomial_coefficient_log(-0.5, 0.49)); EXPECT_NO_THROW(binomial_coefficient_log(10, -0.9)); - EXPECT_FLOAT_EQ(binomial_coefficient_log(0, -1), -inf); + EXPECT_FLOAT_EQ(binomial_coefficient_log(0, -1), -INFTY); EXPECT_FLOAT_EQ(binomial_coefficient_log(-1, 0), 0); - EXPECT_FLOAT_EQ(binomial_coefficient_log(-1, -0.3), inf); - EXPECT_FLOAT_EQ(binomial_coefficient_log(0.3, -1), -inf); - EXPECT_FLOAT_EQ(binomial_coefficient_log(5.0, 6.0), -inf); + EXPECT_FLOAT_EQ(binomial_coefficient_log(-1, -0.3), INFTY); + EXPECT_FLOAT_EQ(binomial_coefficient_log(0.3, -1), -INFTY); + EXPECT_FLOAT_EQ(binomial_coefficient_log(5.0, 6.0), -INFTY); } diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 9b02298ed74..5b906159cab 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -1,7 +1,6 @@ -#include #include +#include #include -#include #include #include #include @@ -15,7 +14,7 @@ struct TestValue { double dk; }; -constexpr double NaN = std::numeric_limits::quiet_NaN(); +constexpr double NaN = stan::math::NOT_A_NUMBER; // Test values generated in Mathematica, reproducible notebook at // https://www.wolframcloud.com/obj/martin.modrak/Published/binomial_coefficient_log.nb // Mathematica Code reproduced below for convenience: @@ -211,7 +210,7 @@ TEST(MathFunctions, binomial_coefficient_log_precomputed) { using stan::test::expect_near_rel; for (TestValue t : testValues) { - std::ostringstream msg; + std::stringstream msg; msg << std::setprecision(22) << "n = " << t.n << ", k = " << t.k; var n(t.n); From 7922dd5045db851e500b114c084b9dd66a8df31c Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 14 Jan 2020 15:10:57 +0000 Subject: [PATCH 08/26] [Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (tags/RELEASE_500/final) --- stan/math/prim/fun/binomial_coefficient_log.hpp | 4 ++-- test/unit/math/prim/fun/binomial_coefficient_log_test.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index cd39e80b0fe..7365f9aab55 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -90,10 +90,10 @@ inline return_type_t binomial_coefficient_log(const T_N N, } if (N / 2 < n) { return binomial_coefficient_log(N, N - n); - } + } if (N_plus_1 < lgamma_stirling_diff_useful) { return lgamma(N_plus_1) - lgamma(n + 1) - lgamma(N_plus_1 - n); - } + } return -lbeta(N - n + 1, n + 1) - log(N_plus_1); } diff --git a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp index e372b6815d9..ef379939828 100644 --- a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp @@ -50,8 +50,8 @@ TEST(MathFunctions, binomial_coefficient_log_nan) { } TEST(MathFunctions, binomial_coefficient_log_errors_edge_cases) { - using stan::math::binomial_coefficient_log; using stan::math::INFTY; + using stan::math::binomial_coefficient_log; EXPECT_NO_THROW(binomial_coefficient_log(10, 11)); EXPECT_THROW(binomial_coefficient_log(10, 11.01), std::domain_error); From 6d6ac35f8f764948be1da2f82c790f7709802923 Mon Sep 17 00:00:00 2001 From: martinmodrak Date: Wed, 15 Jan 2020 14:57:08 +0100 Subject: [PATCH 09/26] Fixed lint error --- test/unit/math/prim/fun/binomial_coefficient_log_test.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp index ef379939828..29c2366badf 100644 --- a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp @@ -25,9 +25,9 @@ TEST(MathFunctions, binomial_coefficient_log) { EXPECT_EQ(binomial_coefficient_log(10000, 0), 0); EXPECT_EQ(binomial_coefficient_log(10, 11), - -std::numeric_limits::infinity()); + stan::math::NEGATIVE_INFTY); EXPECT_EQ(binomial_coefficient_log(10, -1), - -std::numeric_limits::infinity()); + stan::math::NEGATIVE_INFTY); for (int n = 0; n < 1010; ++n) { test_binom_coefficient(1010, n); @@ -42,7 +42,7 @@ TEST(MathFunctions, binomial_coefficient_log) { } TEST(MathFunctions, binomial_coefficient_log_nan) { - double nan = std::numeric_limits::quiet_NaN(); + double nan = stan::math::NOT_A_NUMBER; EXPECT_TRUE(std::isnan(stan::math::binomial_coefficient_log(2.0, nan))); EXPECT_TRUE(std::isnan(stan::math::binomial_coefficient_log(nan, 2.0))); From 3a320df6619b3b107991214214bdaf3f43ff3646 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 15 Jan 2020 14:04:43 +0000 Subject: [PATCH 10/26] [Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (tags/RELEASE_500/final) --- test/unit/math/prim/fun/binomial_coefficient_log_test.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp index 29c2366badf..e5df58947b9 100644 --- a/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/prim/fun/binomial_coefficient_log_test.cpp @@ -24,10 +24,8 @@ TEST(MathFunctions, binomial_coefficient_log) { EXPECT_EQ(binomial_coefficient_log(50, 0), 0); EXPECT_EQ(binomial_coefficient_log(10000, 0), 0); - EXPECT_EQ(binomial_coefficient_log(10, 11), - stan::math::NEGATIVE_INFTY); - EXPECT_EQ(binomial_coefficient_log(10, -1), - stan::math::NEGATIVE_INFTY); + EXPECT_EQ(binomial_coefficient_log(10, 11), stan::math::NEGATIVE_INFTY); + EXPECT_EQ(binomial_coefficient_log(10, -1), stan::math::NEGATIVE_INFTY); for (int n = 0; n < 1010; ++n) { test_binom_coefficient(1010, n); From 6ebef3ed00361153d7a9200be5002cf19834ac2c Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Thu, 16 Jan 2020 13:40:20 +0100 Subject: [PATCH 11/26] Removed problematic constexpr --- test/unit/math/rev/fun/binomial_coefficient_log_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 5b906159cab..dc7fb0a3091 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -14,7 +14,7 @@ struct TestValue { double dk; }; -constexpr double NaN = stan::math::NOT_A_NUMBER; +const double NaN = stan::math::NOT_A_NUMBER; // Test values generated in Mathematica, reproducible notebook at // https://www.wolframcloud.com/obj/martin.modrak/Published/binomial_coefficient_log.nb // Mathematica Code reproduced below for convenience: From 74d28243baa38f89f4fc7e361759415947577b04 Mon Sep 17 00:00:00 2001 From: martinmodrak Date: Wed, 29 Jan 2020 09:22:30 +0100 Subject: [PATCH 12/26] Updated and cleaned up test generating code --- .../rev/fun/binomial_coefficient_log_test.cpp | 245 +++++++----------- 1 file changed, 87 insertions(+), 158 deletions(-) diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index dc7fb0a3091..c7ca71df168 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -19,184 +19,113 @@ const double NaN = stan::math::NOT_A_NUMBER; // https://www.wolframcloud.com/obj/martin.modrak/Published/binomial_coefficient_log.nb // Mathematica Code reproduced below for convenience: -// bclog[n_,k_]:= LogGamma[n + 1] - LogGamma[k + 1] - LogGamma[n + 1 -k ] -// dbclogdn[n_,k_]= D[bclog[n,k],n]; -// dbclogdk[n_,k_]= D[bclog[n,k],k]; - +// toCString[x_] := ToString[CForm[N[x, 24]]]; // singleTest[x_,y_]:= Module[{val, cdn,cdk},{ -// val = N[bclog[x,y],24]; -// cdn = If[x > 10^6 || y > 10^6,"NaN", CForm[N[dbclogdn[x,y],24]]]; -// cdk = If[x > 10^6 || y > 10^6,"NaN", CForm[N[dbclogdk[x,y],24]]]; -// WriteString[out," {",CForm[x],",",CForm[y],",", -// CForm[val],",",cdn,",",cdk,"},"] +// val = toCString[bclog[x,y]]; +// cdn = If[x > 10^6 || y > 10^6,"NaN", toCString[dbclogdn[x,y]]]; +// cdk = If[x > 10^6 || y > 10^6,"NaN", toCString[dbclogdk[x,y]]]; +// StringJoin[" {",toCString[x],",",toCString[y],",", +// val,",",cdn,",",cdk,"},\n"] // }]; - -// out = OpenWrite["bc_test.txt"]; // ns= {-0.1,3*10^-5,2*10^-3,1,8, 1325,845*10^3}; // ratios = {-1,10^- 10,10^-5,10^-2,1/5,1/2,1-3*10^-2,1-6*10^-8, 1 -3*10^-9,2}; -// WriteString[out, "std::vector testValues = {"]; +// out = "std::vector testValues = {\n"; // For[i = 1, i <= Length[ns], i++, { // For[j = 1, j <= Length[ratios], j++, { // cn = ns[[i]]; // ck = If[ratios[[j]] < 0,-9/10, // If[ratios[[j]] > 1,cn + 9/10,cn * ratios[[j]] ]]; -// singleTest[cn,ck]; +// out = StringJoin[out, singleTest[cn,ck]]; // }] // }] - // extremeNs = {3*10^15+1/2,10^20 + 1/2}; // lowKs = {3, 100, 12895}; // For[i = 1, i <= Length[extremeNs], i++, { // For[j = 1, j <= Length[lowKs], j++, { // cn = extremeNs[[i]]; // ck = lowKs[[j]]; -// singleTest[cn,ck]; +// out = StringJoin[out,singleTest[cn,ck]]; // }] // }] - -// WriteString[out,"};"]; -// Close[out]; -// FilePrint[%] +// out = StringJoin[out,"};\n"]; +// out std::vector testValues = { - {-0.1, -0.9, -2.1152525390850903, -1.039918383240913, 10.708746373704937}, - {-0.1, -1.0000000000000001e-11, 1.7771062399418724e-12, - -1.922551007282891e-11, -0.17771128500984368}, - {-0.1, -1.0000000000000002e-6, 1.7770950132278696e-7, - -1.9225383585119715e-6, -0.1777077175718912}, - {-0.1, -0.001, 0.00017592768030333383, -0.0019209406823855746, - -0.17414420715602041}, - {-0.1, -0.020000000000000004, 0.0028416866729211523, -0.037823152614322564, - -0.10649980051731356}, - {-0.1, -0.05, 0.0044386492587971844, -0.09231733374172624, 0.}, - {-0.1, -0.097, 0.0005170837568402953, -0.17276563502587905, - 0.16701237954910164}, - {-0.1, -0.09999999400000001, 1.0662676433991578e-9, -0.17771127517591423, - 0.17771126364067447}, - {-0.1, -0.09999999970000001, 5.331331690149795e-11, -0.17771128455203855, - 0.17771128397527658}, - {-0.1, 0.8, -2.1152525390850907, 9.668827990464031, -10.708746373704944}, - {0.00003, -0.9, -2.21375637737528044964112, -0.933371118080918307851728, - 10.7799597405306456982503}, - {0.00003, 3.e-15, 1.48040820535563428588846e-19, - 4.9345858390686036409092e-15, 0.0000493469401735864488431421}, - {0.00003, 3.e-10, 1.48039340142161568666212e-14, - 4.93458584015035637326682e-10, 0.0000493459532446518755643083}, - {0.00003, 3.e-7, 1.46560412344434242311062e-11, - 4.93458692083243471796274e-7, 0.0000483600013795032175416619}, - {0.00003, 6.e-6, 2.36865312869367134774661e-10, - 9.86921494891292997054538e-6, 0.0000296081641072682815641386}, - {0.00003, 0.000015, 3.70102051348524044535716e-10, - 0.0000246731996398828628243033, 0}, - {0.00003, 0.0000291, 4.30798787797857754693695e-11, - 0.0000478665004969577356362445, -0.0000463861237716491741348152}, - {0.00003, 0.0000299999982, 8.88244870007509649977929e-17, - 0.0000493469372225745165598424, -0.0000493469342618230131147924}, - {0.00003, 0.00002999999991, 4.44122460318735146597007e-18, - 0.0000493469400354117708664247, -0.0000493469398873741956943572}, - {0.00003, 0.90003, -2.21375637737528044964112, 9.84658862244972739039859, - -10.7799597405306456982503}, - {0.002, -0.9, -2.21559326412971099686943, -0.931489785799666354340725, - 10.7813141298573216195769}, - {0.002, 2.e-13, 6.57013709556564677684856e-16, 3.2802775880282756733109e-13, - 0.00328506854745431610233807}, - {0.002, 2.e-8, 6.57007139476544583161173e-11, 3.28027763585130905255152e-8, - 0.00328500284665411601964384}, - {0.002, 0.00002, 6.50443564072189575550994e-8, - 0.0000328032541239782363672409, 0.00321936709346449075328701}, - {0.002, 0.0004, 1.05122171458287350859763e-6, 0.000656246880415250668008776, - 0.00197104030081914294655284}, - {0.002, 0.001, 1.64253373496215313253469e-6, 0.00164133545687894154061287, - 0}, - {0.002, 0.00194, 1.91190982195918976356429e-7, 0.00318637683127152008217264, - -0.00308796419928270612987841}, - {0.002, 0.00199999988, 3.94208202120835082737684e-13, - 0.00328506835071924267003294, -0.00328506815390258737114863}, - {0.002, 0.001999999994, 1.97104112295366725515452e-14, - 0.00328506853824172627346441, -0.00328506852840089350933798}, - {0.002, 0.902, -2.21559326412971099686943, 9.84982434405765526523621, - -10.7813141298573216195769}, - {1, -0.9, -2.85558226198351740582195, -0.459715615539276790357555, - 11.3062548910488207249193}, - {1, 1.e-10, 9.9999999988550659331851e-11, 6.44934066868432126789198e-11, - 0.99999999977101318664035}, - {1, 0.00001, 9.99988550692664559909352e-6, 6.44936087425490392714312e-6, - 0.999977101418661870834803}, - {1, 0.01, 0.00988583703486131052627978, 0.00646962905305218281638255, - 0.977200163914089454058067}, - {1, 0.2, 0.156457962917688016707705, 0.137792901804605598784786, - 0.57403132988604983615591}, - {1, 0.5, 0.241564475270490444691037, 0.386294361119890618834464, 0}, - {1, 0.97, 0.028978328236256312960776, 0.951705422383897641079507, - -0.932173296099201807475084}, - {1, 0.99999994, 5.99999958782374313463811e-8, 0.999999901303960316511031, - -0.999999862607915578212576}, - {1, 0.999999997, 2.99999998969559340736596e-9, 0.999999995065197810273833, - -0.999999993130395607910641}, - {1, 1.9, -2.85558226198351740582195, 10.8465392755095439345617, - -11.3062548910488207249193}, - {8, -0.9, -4.22528965320883461943031, -0.100538838650771402252215, - 12.6649352570174581939568}, - {8, 8.e-10, 2.17428571372173153982474e-9, 9.4009611759639002250559e-11, - 2.71785714144718599267395}, - {8, 0.00008, 0.000217422931805073420417006, 9.4010053144199442029858e-6, - 2.71771615481909065227871}, - {8, 0.08, 0.211982267378255838975509, 0.00944537775278953163286055, - 2.58399543820453334937797}, - {8, 1.6, 2.90678606291134283426263, 0.208248082071609202561846, - 1.1813459431105245420685}, - {8, 4, 4.24849524204935898912334, 0.634523809523809523809524, 0}, - {8, 7.76, 0.606274586245453651115361, 2.38013515708154679097544, - -2.35152741320850411987141}, - {8, 7.99999952, 1.30457122553768403613331e-6, 2.71785635328906772378496, - -2.71785629688329908165944}, - {8, 7.999999976, 6.52285709209869625945566e-8, 2.71785710337872594517018, - -2.71785710055843758854095}, - {8, 8.9, -4.22528965320883461943031, 12.5643964183666867917046, - -12.6649352570174581939568}, - {1325, -0.9, -8.72360867216657209762532, -0.000678758619778166713080174, - 17.6139787484752694858826}, - {1325, 1.325e-7, 1.02909578020477960435539e-6, - 9.99622736492339027462876e-11, 7.76676049629224235717761}, - {1325, 0.01325, 0.102766042691370430370992, 9.99627732703511703537678e-6, - 7.74516389164454804001553}, - {1325, 13.25, 71.9898321274090629975055, 0.0100465251153381679080642, - 4.55823951796927536337622}, - {1325, 265, 659.435649329029419323398, 0.223049238391457549660772, - 1.38488037927413480754602}, - {1325, 662.5, 914.599450340845275100724, 0.692769964468769118810416, 0}, - {1325, 1285.25, 175.786260651191862665015, 3.49440932918478008636475, - -3.46396178959974953213922}, - {1325, 1324.9999205, 0.000617452276410452190170437, - 7.76662994968438931022129, -7.76662988970702332503659}, - {1325, 1324.999996025, 0.0000308728608380968862741097, - 7.76675417575202484528552, -7.76675417275315663146178}, - {1325, 1325.9, -8.72360867216657209762532, 17.6132999898554913191696, - -17.6139787484752694858826}, - {845000, -0.9, -14.5350963792733464918229, -1.06508755996070778708647e-6, - 24.0708485035538072767331}, - {845000, 0.0000845, 0.00120194816738712136581358, - 9.99999408334257028001921e-11, 14.2241691745103886002024}, - {845000, 8.45, 103.738303827736743600251, 0.0000100000440831167345507131, - 11.454909927882486225397}, - {845000, 8450, 47315.8616457576200611209, 0.0100503298765747570022521, - 4.59506127739683575088585}, - {845000, 169000, 422833.221695496506553128, 0.223143403385281321887929, - 1.3862921421877147496422}, - {845000, 422500, 585702.318235552114086514, 0.693146588844319105852931, 0}, - {845000, 819650, 113851.158132678562120685, 3.50653876529964054840506, - -3.476079576115418785061}, - {845000, 844999.9493, 0.719108776819481762797449, 14.1438653554064391526676, - -14.1438652954064728556142}, - {845000, 844999.997465, 0.036053342347290003917417, - 14.2201459621965153509889, -14.2201459591965171216362}, - {845000, 845000.9, -14.5350963792733464918229, 24.0708474384662473160253, - -24.0708485035538072767331}, - {3.0000000000000005e15, 3, 105.120406581508329354183, NaN, NaN}, - {3.0000000000000005e15, 100, 3199.9994928023143716891, NaN, NaN}, - {3.0000000000000005e15, 12895, 350387.524360588370915867, NaN, NaN}, - {1.e20, 3, 136.363346110414686040252, NaN, NaN}, - {1.e20, 100, 4241.4308104325278778429, NaN, NaN}, - {1.e20, 12895, 484680.092769031900296942, NaN, NaN}, + {-0.1,-0.9,-2.1152525390850903,-1.039918383240913,10.708746373704937}, + {-0.1,-1.0000000000000001e-11,1.7771062399418724e-12,-1.922551007282891e-11,-0.17771128500984368}, + {-0.1,-1.0000000000000002e-6,1.7770950132278696e-7,-1.9225383585119715e-6,-0.1777077175718912}, + {-0.1,-0.001,0.00017592768030333383,-0.0019209406823855746,-0.17414420715602041}, + {-0.1,-0.020000000000000004,0.0028416866729211523,-0.037823152614322564,-0.10649980051731356}, + {-0.1,-0.05,0.0044386492587971844,-0.09231733374172624,0.}, + {-0.1,-0.097,0.0005170837568402953,-0.17276563502587905,0.16701237954910164}, + {-0.1,-0.09999999400000001,1.0662676433991578e-9,-0.17771127517591423,0.17771126364067447}, + {-0.1,-0.09999999970000001,5.331331690149795e-11,-0.17771128455203855,0.17771128397527658}, + {-0.1,0.8,-2.1152525390850907,9.668827990464031,-10.708746373704944}, + {0.00003,-0.9,-2.21375637737528044964112,-0.933371118080918307851728,10.7799597405306456982503}, + {0.00003,3.e-15,1.48040820535563428588846e-19,4.9345858390686036409092e-15,0.0000493469401735864488431421}, + {0.00003,3.e-10,1.48039340142161568666212e-14,4.93458584015035637326682e-10,0.0000493459532446518755643083}, + {0.00003,3.e-7,1.46560412344434242311062e-11,4.93458692083243471796274e-7,0.0000483600013795032175416619}, + {0.00003,6.e-6,2.36865312869367134774661e-10,9.86921494891292997054538e-6,0.0000296081641072682815641386}, + {0.00003,0.000015,3.70102051348524044535716e-10,0.0000246731996398828628243033,0}, + {0.00003,0.0000291,4.30798787797857754693695e-11,0.0000478665004969577356362445,-0.0000463861237716491741348152}, + {0.00003,0.0000299999982,8.88244870007509649977929e-17,0.0000493469372225745165598424,-0.0000493469342618230131147924}, + {0.00003,0.00002999999991,4.44122460318735146597007e-18,0.0000493469400354117708664247,-0.0000493469398873741956943572}, + {0.00003,0.90003,-2.21375637737528044964112,9.84658862244972739039859,-10.7799597405306456982503}, + {0.002,-0.9,-2.21559326412971099686943,-0.931489785799666354340725,10.7813141298573216195769}, + {0.002,2.e-13,6.57013709556564677684856e-16,3.2802775880282756733109e-13,0.00328506854745431610233807}, + {0.002,2.e-8,6.57007139476544583161173e-11,3.28027763585130905255152e-8,0.00328500284665411601964384}, + {0.002,0.00002,6.50443564072189575550994e-8,0.0000328032541239782363672409,0.00321936709346449075328701}, + {0.002,0.0004,1.05122171458287350859763e-6,0.000656246880415250668008776,0.00197104030081914294655284}, + {0.002,0.001,1.64253373496215313253469e-6,0.00164133545687894154061287,0}, + {0.002,0.00194,1.91190982195918976356429e-7,0.00318637683127152008217264,-0.00308796419928270612987841}, + {0.002,0.00199999988,3.94208202120835082737684e-13,0.00328506835071924267003294,-0.00328506815390258737114863}, + {0.002,0.001999999994,1.97104112295366725515452e-14,0.00328506853824172627346441,-0.00328506852840089350933798}, + {0.002,0.902,-2.21559326412971099686943,9.84982434405765526523621,-10.7813141298573216195769}, + {1.,-0.9,-2.85558226198351740582195,-0.459715615539276790357555,11.3062548910488207249193}, + {1.,1.e-10,9.9999999988550659331851e-11,6.44934066868432126789198e-11,0.99999999977101318664035}, + {1.,0.00001,9.99988550692664559909352e-6,6.44936087425490392714312e-6,0.999977101418661870834803}, + {1.,0.01,0.00988583703486131052627978,0.00646962905305218281638255,0.977200163914089454058067}, + {1.,0.2,0.156457962917688016707705,0.137792901804605598784786,0.57403132988604983615591}, + {1.,0.5,0.241564475270490444691037,0.386294361119890618834464,0}, + {1.,0.97,0.028978328236256312960776,0.951705422383897641079507,-0.932173296099201807475084}, + {1.,0.99999994,5.99999958782374313463811e-8,0.999999901303960316511031,-0.999999862607915578212576}, + {1.,0.999999997,2.99999998969559340736596e-9,0.999999995065197810273833,-0.999999993130395607910641}, + {1.,1.9,-2.85558226198351740582195,10.8465392755095439345617,-11.3062548910488207249193}, + {8.,-0.9,-4.22528965320883461943031,-0.100538838650771402252215,12.6649352570174581939568}, + {8.,8.e-10,2.17428571372173153982474e-9,9.4009611759639002250559e-11,2.71785714144718599267395}, + {8.,0.00008,0.000217422931805073420417006,9.4010053144199442029858e-6,2.71771615481909065227871}, + {8.,0.08,0.211982267378255838975509,0.00944537775278953163286055,2.58399543820453334937797}, + {8.,1.6,2.90678606291134283426263,0.208248082071609202561846,1.1813459431105245420685}, + {8.,4.,4.24849524204935898912334,0.634523809523809523809524,0}, + {8.,7.76,0.606274586245453651115361,2.38013515708154679097544,-2.35152741320850411987141}, + {8.,7.99999952,1.30457122553768403613331e-6,2.71785635328906772378496,-2.71785629688329908165944}, + {8.,7.999999976,6.52285709209869625945566e-8,2.71785710337872594517018,-2.71785710055843758854095}, + {8.,8.9,-4.22528965320883461943031,12.5643964183666867917046,-12.6649352570174581939568}, + {1325.,-0.9,-8.72360867216657209762532,-0.000678758619778166713080174,17.6139787484752694858826}, + {1325.,1.325e-7,1.02909578020477960435539e-6,9.99622736492339027462876e-11,7.76676049629224235717761}, + {1325.,0.01325,0.102766042691370430370992,9.99627732703511703537678e-6,7.74516389164454804001553}, + {1325.,13.25,71.9898321274090629975055,0.0100465251153381679080642,4.55823951796927536337622}, + {1325.,265.,659.435649329029419323398,0.223049238391457549660772,1.38488037927413480754602}, + {1325.,662.5,914.599450340845275100724,0.692769964468769118810416,0}, + {1325.,1285.25,175.786260651191862665015,3.49440932918478008636475,-3.46396178959974953213922}, + {1325.,1324.9999205,0.000617452276410452190170437,7.76662994968438931022129,-7.76662988970702332503659}, + {1325.,1324.999996025,0.0000308728608380968862741097,7.76675417575202484528552,-7.76675417275315663146178}, + {1325.,1325.9,-8.72360867216657209762532,17.6132999898554913191696,-17.6139787484752694858826}, + {845000.,-0.9,-14.5350963792733464918229,-1.06508755996070778708647e-6,24.0708485035538072767331}, + {845000.,0.0000845,0.00120194816738712136581358,9.99999408334257028001921e-11,14.2241691745103886002024}, + {845000.,8.45,103.738303827736743600251,0.0000100000440831167345507131,11.454909927882486225397}, + {845000.,8450.,47315.8616457576200611209,0.0100503298765747570022521,4.59506127739683575088585}, + {845000.,169000.,422833.221695496506553128,0.223143403385281321887929,1.3862921421877147496422}, + {845000.,422500.,585702.318235552114086514,0.693146588844319105852931,0}, + {845000.,819650.,113851.158132678562120685,3.50653876529964054840506,-3.476079576115418785061}, + {845000.,844999.9493,0.719108776819481762797449,14.1438653554064391526676,-14.1438652954064728556142}, + {845000.,844999.997465,0.036053342347290003917417,14.2201459621965153509889,-14.2201459591965171216362}, + {845000.,845000.9,-14.5350963792733464918229,24.0708474384662473160253,-24.0708485035538072767331}, + {3.0000000000000005e15,3.,105.120406581508329354183,NaN,NaN}, + {3.0000000000000005e15,100.,3199.9994928023143716891,NaN,NaN}, + {3.0000000000000005e15,12895.,350387.524360588370915867,NaN,NaN}, + {1.000000000000000000005e20,3.,136.363346110414686040252,NaN,NaN}, + {1.000000000000000000005e20,100.,4241.4308104325278778429,NaN,NaN}, + {1.000000000000000000005e20,12895.,484680.092769031900296942,NaN,NaN}, }; } // namespace binomial_coefficient_log_test_internal From dd59b06891c47cd4e7087d2bacb3d724f5933ee1 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 29 Jan 2020 08:25:42 +0000 Subject: [PATCH 13/26] [Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (tags/RELEASE_500/final) --- .../rev/fun/binomial_coefficient_log_test.cpp | 216 ++++++++++++------ 1 file changed, 140 insertions(+), 76 deletions(-) diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index c7ca71df168..34fb5344c91 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -50,82 +50,146 @@ const double NaN = stan::math::NOT_A_NUMBER; // out = StringJoin[out,"};\n"]; // out std::vector testValues = { - {-0.1,-0.9,-2.1152525390850903,-1.039918383240913,10.708746373704937}, - {-0.1,-1.0000000000000001e-11,1.7771062399418724e-12,-1.922551007282891e-11,-0.17771128500984368}, - {-0.1,-1.0000000000000002e-6,1.7770950132278696e-7,-1.9225383585119715e-6,-0.1777077175718912}, - {-0.1,-0.001,0.00017592768030333383,-0.0019209406823855746,-0.17414420715602041}, - {-0.1,-0.020000000000000004,0.0028416866729211523,-0.037823152614322564,-0.10649980051731356}, - {-0.1,-0.05,0.0044386492587971844,-0.09231733374172624,0.}, - {-0.1,-0.097,0.0005170837568402953,-0.17276563502587905,0.16701237954910164}, - {-0.1,-0.09999999400000001,1.0662676433991578e-9,-0.17771127517591423,0.17771126364067447}, - {-0.1,-0.09999999970000001,5.331331690149795e-11,-0.17771128455203855,0.17771128397527658}, - {-0.1,0.8,-2.1152525390850907,9.668827990464031,-10.708746373704944}, - {0.00003,-0.9,-2.21375637737528044964112,-0.933371118080918307851728,10.7799597405306456982503}, - {0.00003,3.e-15,1.48040820535563428588846e-19,4.9345858390686036409092e-15,0.0000493469401735864488431421}, - {0.00003,3.e-10,1.48039340142161568666212e-14,4.93458584015035637326682e-10,0.0000493459532446518755643083}, - {0.00003,3.e-7,1.46560412344434242311062e-11,4.93458692083243471796274e-7,0.0000483600013795032175416619}, - {0.00003,6.e-6,2.36865312869367134774661e-10,9.86921494891292997054538e-6,0.0000296081641072682815641386}, - {0.00003,0.000015,3.70102051348524044535716e-10,0.0000246731996398828628243033,0}, - {0.00003,0.0000291,4.30798787797857754693695e-11,0.0000478665004969577356362445,-0.0000463861237716491741348152}, - {0.00003,0.0000299999982,8.88244870007509649977929e-17,0.0000493469372225745165598424,-0.0000493469342618230131147924}, - {0.00003,0.00002999999991,4.44122460318735146597007e-18,0.0000493469400354117708664247,-0.0000493469398873741956943572}, - {0.00003,0.90003,-2.21375637737528044964112,9.84658862244972739039859,-10.7799597405306456982503}, - {0.002,-0.9,-2.21559326412971099686943,-0.931489785799666354340725,10.7813141298573216195769}, - {0.002,2.e-13,6.57013709556564677684856e-16,3.2802775880282756733109e-13,0.00328506854745431610233807}, - {0.002,2.e-8,6.57007139476544583161173e-11,3.28027763585130905255152e-8,0.00328500284665411601964384}, - {0.002,0.00002,6.50443564072189575550994e-8,0.0000328032541239782363672409,0.00321936709346449075328701}, - {0.002,0.0004,1.05122171458287350859763e-6,0.000656246880415250668008776,0.00197104030081914294655284}, - {0.002,0.001,1.64253373496215313253469e-6,0.00164133545687894154061287,0}, - {0.002,0.00194,1.91190982195918976356429e-7,0.00318637683127152008217264,-0.00308796419928270612987841}, - {0.002,0.00199999988,3.94208202120835082737684e-13,0.00328506835071924267003294,-0.00328506815390258737114863}, - {0.002,0.001999999994,1.97104112295366725515452e-14,0.00328506853824172627346441,-0.00328506852840089350933798}, - {0.002,0.902,-2.21559326412971099686943,9.84982434405765526523621,-10.7813141298573216195769}, - {1.,-0.9,-2.85558226198351740582195,-0.459715615539276790357555,11.3062548910488207249193}, - {1.,1.e-10,9.9999999988550659331851e-11,6.44934066868432126789198e-11,0.99999999977101318664035}, - {1.,0.00001,9.99988550692664559909352e-6,6.44936087425490392714312e-6,0.999977101418661870834803}, - {1.,0.01,0.00988583703486131052627978,0.00646962905305218281638255,0.977200163914089454058067}, - {1.,0.2,0.156457962917688016707705,0.137792901804605598784786,0.57403132988604983615591}, - {1.,0.5,0.241564475270490444691037,0.386294361119890618834464,0}, - {1.,0.97,0.028978328236256312960776,0.951705422383897641079507,-0.932173296099201807475084}, - {1.,0.99999994,5.99999958782374313463811e-8,0.999999901303960316511031,-0.999999862607915578212576}, - {1.,0.999999997,2.99999998969559340736596e-9,0.999999995065197810273833,-0.999999993130395607910641}, - {1.,1.9,-2.85558226198351740582195,10.8465392755095439345617,-11.3062548910488207249193}, - {8.,-0.9,-4.22528965320883461943031,-0.100538838650771402252215,12.6649352570174581939568}, - {8.,8.e-10,2.17428571372173153982474e-9,9.4009611759639002250559e-11,2.71785714144718599267395}, - {8.,0.00008,0.000217422931805073420417006,9.4010053144199442029858e-6,2.71771615481909065227871}, - {8.,0.08,0.211982267378255838975509,0.00944537775278953163286055,2.58399543820453334937797}, - {8.,1.6,2.90678606291134283426263,0.208248082071609202561846,1.1813459431105245420685}, - {8.,4.,4.24849524204935898912334,0.634523809523809523809524,0}, - {8.,7.76,0.606274586245453651115361,2.38013515708154679097544,-2.35152741320850411987141}, - {8.,7.99999952,1.30457122553768403613331e-6,2.71785635328906772378496,-2.71785629688329908165944}, - {8.,7.999999976,6.52285709209869625945566e-8,2.71785710337872594517018,-2.71785710055843758854095}, - {8.,8.9,-4.22528965320883461943031,12.5643964183666867917046,-12.6649352570174581939568}, - {1325.,-0.9,-8.72360867216657209762532,-0.000678758619778166713080174,17.6139787484752694858826}, - {1325.,1.325e-7,1.02909578020477960435539e-6,9.99622736492339027462876e-11,7.76676049629224235717761}, - {1325.,0.01325,0.102766042691370430370992,9.99627732703511703537678e-6,7.74516389164454804001553}, - {1325.,13.25,71.9898321274090629975055,0.0100465251153381679080642,4.55823951796927536337622}, - {1325.,265.,659.435649329029419323398,0.223049238391457549660772,1.38488037927413480754602}, - {1325.,662.5,914.599450340845275100724,0.692769964468769118810416,0}, - {1325.,1285.25,175.786260651191862665015,3.49440932918478008636475,-3.46396178959974953213922}, - {1325.,1324.9999205,0.000617452276410452190170437,7.76662994968438931022129,-7.76662988970702332503659}, - {1325.,1324.999996025,0.0000308728608380968862741097,7.76675417575202484528552,-7.76675417275315663146178}, - {1325.,1325.9,-8.72360867216657209762532,17.6132999898554913191696,-17.6139787484752694858826}, - {845000.,-0.9,-14.5350963792733464918229,-1.06508755996070778708647e-6,24.0708485035538072767331}, - {845000.,0.0000845,0.00120194816738712136581358,9.99999408334257028001921e-11,14.2241691745103886002024}, - {845000.,8.45,103.738303827736743600251,0.0000100000440831167345507131,11.454909927882486225397}, - {845000.,8450.,47315.8616457576200611209,0.0100503298765747570022521,4.59506127739683575088585}, - {845000.,169000.,422833.221695496506553128,0.223143403385281321887929,1.3862921421877147496422}, - {845000.,422500.,585702.318235552114086514,0.693146588844319105852931,0}, - {845000.,819650.,113851.158132678562120685,3.50653876529964054840506,-3.476079576115418785061}, - {845000.,844999.9493,0.719108776819481762797449,14.1438653554064391526676,-14.1438652954064728556142}, - {845000.,844999.997465,0.036053342347290003917417,14.2201459621965153509889,-14.2201459591965171216362}, - {845000.,845000.9,-14.5350963792733464918229,24.0708474384662473160253,-24.0708485035538072767331}, - {3.0000000000000005e15,3.,105.120406581508329354183,NaN,NaN}, - {3.0000000000000005e15,100.,3199.9994928023143716891,NaN,NaN}, - {3.0000000000000005e15,12895.,350387.524360588370915867,NaN,NaN}, - {1.000000000000000000005e20,3.,136.363346110414686040252,NaN,NaN}, - {1.000000000000000000005e20,100.,4241.4308104325278778429,NaN,NaN}, - {1.000000000000000000005e20,12895.,484680.092769031900296942,NaN,NaN}, + {-0.1, -0.9, -2.1152525390850903, -1.039918383240913, 10.708746373704937}, + {-0.1, -1.0000000000000001e-11, 1.7771062399418724e-12, + -1.922551007282891e-11, -0.17771128500984368}, + {-0.1, -1.0000000000000002e-6, 1.7770950132278696e-7, + -1.9225383585119715e-6, -0.1777077175718912}, + {-0.1, -0.001, 0.00017592768030333383, -0.0019209406823855746, + -0.17414420715602041}, + {-0.1, -0.020000000000000004, 0.0028416866729211523, -0.037823152614322564, + -0.10649980051731356}, + {-0.1, -0.05, 0.0044386492587971844, -0.09231733374172624, 0.}, + {-0.1, -0.097, 0.0005170837568402953, -0.17276563502587905, + 0.16701237954910164}, + {-0.1, -0.09999999400000001, 1.0662676433991578e-9, -0.17771127517591423, + 0.17771126364067447}, + {-0.1, -0.09999999970000001, 5.331331690149795e-11, -0.17771128455203855, + 0.17771128397527658}, + {-0.1, 0.8, -2.1152525390850907, 9.668827990464031, -10.708746373704944}, + {0.00003, -0.9, -2.21375637737528044964112, -0.933371118080918307851728, + 10.7799597405306456982503}, + {0.00003, 3.e-15, 1.48040820535563428588846e-19, + 4.9345858390686036409092e-15, 0.0000493469401735864488431421}, + {0.00003, 3.e-10, 1.48039340142161568666212e-14, + 4.93458584015035637326682e-10, 0.0000493459532446518755643083}, + {0.00003, 3.e-7, 1.46560412344434242311062e-11, + 4.93458692083243471796274e-7, 0.0000483600013795032175416619}, + {0.00003, 6.e-6, 2.36865312869367134774661e-10, + 9.86921494891292997054538e-6, 0.0000296081641072682815641386}, + {0.00003, 0.000015, 3.70102051348524044535716e-10, + 0.0000246731996398828628243033, 0}, + {0.00003, 0.0000291, 4.30798787797857754693695e-11, + 0.0000478665004969577356362445, -0.0000463861237716491741348152}, + {0.00003, 0.0000299999982, 8.88244870007509649977929e-17, + 0.0000493469372225745165598424, -0.0000493469342618230131147924}, + {0.00003, 0.00002999999991, 4.44122460318735146597007e-18, + 0.0000493469400354117708664247, -0.0000493469398873741956943572}, + {0.00003, 0.90003, -2.21375637737528044964112, 9.84658862244972739039859, + -10.7799597405306456982503}, + {0.002, -0.9, -2.21559326412971099686943, -0.931489785799666354340725, + 10.7813141298573216195769}, + {0.002, 2.e-13, 6.57013709556564677684856e-16, 3.2802775880282756733109e-13, + 0.00328506854745431610233807}, + {0.002, 2.e-8, 6.57007139476544583161173e-11, 3.28027763585130905255152e-8, + 0.00328500284665411601964384}, + {0.002, 0.00002, 6.50443564072189575550994e-8, + 0.0000328032541239782363672409, 0.00321936709346449075328701}, + {0.002, 0.0004, 1.05122171458287350859763e-6, 0.000656246880415250668008776, + 0.00197104030081914294655284}, + {0.002, 0.001, 1.64253373496215313253469e-6, 0.00164133545687894154061287, + 0}, + {0.002, 0.00194, 1.91190982195918976356429e-7, 0.00318637683127152008217264, + -0.00308796419928270612987841}, + {0.002, 0.00199999988, 3.94208202120835082737684e-13, + 0.00328506835071924267003294, -0.00328506815390258737114863}, + {0.002, 0.001999999994, 1.97104112295366725515452e-14, + 0.00328506853824172627346441, -0.00328506852840089350933798}, + {0.002, 0.902, -2.21559326412971099686943, 9.84982434405765526523621, + -10.7813141298573216195769}, + {1., -0.9, -2.85558226198351740582195, -0.459715615539276790357555, + 11.3062548910488207249193}, + {1., 1.e-10, 9.9999999988550659331851e-11, 6.44934066868432126789198e-11, + 0.99999999977101318664035}, + {1., 0.00001, 9.99988550692664559909352e-6, 6.44936087425490392714312e-6, + 0.999977101418661870834803}, + {1., 0.01, 0.00988583703486131052627978, 0.00646962905305218281638255, + 0.977200163914089454058067}, + {1., 0.2, 0.156457962917688016707705, 0.137792901804605598784786, + 0.57403132988604983615591}, + {1., 0.5, 0.241564475270490444691037, 0.386294361119890618834464, 0}, + {1., 0.97, 0.028978328236256312960776, 0.951705422383897641079507, + -0.932173296099201807475084}, + {1., 0.99999994, 5.99999958782374313463811e-8, 0.999999901303960316511031, + -0.999999862607915578212576}, + {1., 0.999999997, 2.99999998969559340736596e-9, 0.999999995065197810273833, + -0.999999993130395607910641}, + {1., 1.9, -2.85558226198351740582195, 10.8465392755095439345617, + -11.3062548910488207249193}, + {8., -0.9, -4.22528965320883461943031, -0.100538838650771402252215, + 12.6649352570174581939568}, + {8., 8.e-10, 2.17428571372173153982474e-9, 9.4009611759639002250559e-11, + 2.71785714144718599267395}, + {8., 0.00008, 0.000217422931805073420417006, 9.4010053144199442029858e-6, + 2.71771615481909065227871}, + {8., 0.08, 0.211982267378255838975509, 0.00944537775278953163286055, + 2.58399543820453334937797}, + {8., 1.6, 2.90678606291134283426263, 0.208248082071609202561846, + 1.1813459431105245420685}, + {8., 4., 4.24849524204935898912334, 0.634523809523809523809524, 0}, + {8., 7.76, 0.606274586245453651115361, 2.38013515708154679097544, + -2.35152741320850411987141}, + {8., 7.99999952, 1.30457122553768403613331e-6, 2.71785635328906772378496, + -2.71785629688329908165944}, + {8., 7.999999976, 6.52285709209869625945566e-8, 2.71785710337872594517018, + -2.71785710055843758854095}, + {8., 8.9, -4.22528965320883461943031, 12.5643964183666867917046, + -12.6649352570174581939568}, + {1325., -0.9, -8.72360867216657209762532, -0.000678758619778166713080174, + 17.6139787484752694858826}, + {1325., 1.325e-7, 1.02909578020477960435539e-6, + 9.99622736492339027462876e-11, 7.76676049629224235717761}, + {1325., 0.01325, 0.102766042691370430370992, 9.99627732703511703537678e-6, + 7.74516389164454804001553}, + {1325., 13.25, 71.9898321274090629975055, 0.0100465251153381679080642, + 4.55823951796927536337622}, + {1325., 265., 659.435649329029419323398, 0.223049238391457549660772, + 1.38488037927413480754602}, + {1325., 662.5, 914.599450340845275100724, 0.692769964468769118810416, 0}, + {1325., 1285.25, 175.786260651191862665015, 3.49440932918478008636475, + -3.46396178959974953213922}, + {1325., 1324.9999205, 0.000617452276410452190170437, + 7.76662994968438931022129, -7.76662988970702332503659}, + {1325., 1324.999996025, 0.0000308728608380968862741097, + 7.76675417575202484528552, -7.76675417275315663146178}, + {1325., 1325.9, -8.72360867216657209762532, 17.6132999898554913191696, + -17.6139787484752694858826}, + {845000., -0.9, -14.5350963792733464918229, -1.06508755996070778708647e-6, + 24.0708485035538072767331}, + {845000., 0.0000845, 0.00120194816738712136581358, + 9.99999408334257028001921e-11, 14.2241691745103886002024}, + {845000., 8.45, 103.738303827736743600251, 0.0000100000440831167345507131, + 11.454909927882486225397}, + {845000., 8450., 47315.8616457576200611209, 0.0100503298765747570022521, + 4.59506127739683575088585}, + {845000., 169000., 422833.221695496506553128, 0.223143403385281321887929, + 1.3862921421877147496422}, + {845000., 422500., 585702.318235552114086514, 0.693146588844319105852931, + 0}, + {845000., 819650., 113851.158132678562120685, 3.50653876529964054840506, + -3.476079576115418785061}, + {845000., 844999.9493, 0.719108776819481762797449, + 14.1438653554064391526676, -14.1438652954064728556142}, + {845000., 844999.997465, 0.036053342347290003917417, + 14.2201459621965153509889, -14.2201459591965171216362}, + {845000., 845000.9, -14.5350963792733464918229, 24.0708474384662473160253, + -24.0708485035538072767331}, + {3.0000000000000005e15, 3., 105.120406581508329354183, NaN, NaN}, + {3.0000000000000005e15, 100., 3199.9994928023143716891, NaN, NaN}, + {3.0000000000000005e15, 12895., 350387.524360588370915867, NaN, NaN}, + {1.000000000000000000005e20, 3., 136.363346110414686040252, NaN, NaN}, + {1.000000000000000000005e20, 100., 4241.4308104325278778429, NaN, NaN}, + {1.000000000000000000005e20, 12895., 484680.092769031900296942, NaN, NaN}, }; } // namespace binomial_coefficient_log_test_internal From 5bda4f53101be4649c98c871f70aa305f8bb085b Mon Sep 17 00:00:00 2001 From: martinmodrak Date: Wed, 29 Jan 2020 11:11:24 +0100 Subject: [PATCH 14/26] Improved precomputed tests, first pass at formula test --- .../prim/fun/binomial_coefficient_log.hpp | 3 +- .../rev/fun/binomial_coefficient_log_test.cpp | 301 +++++++++--------- 2 files changed, 155 insertions(+), 149 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 7365f9aab55..240a473bbfd 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -94,7 +95,7 @@ inline return_type_t binomial_coefficient_log(const T_N N, if (N_plus_1 < lgamma_stirling_diff_useful) { return lgamma(N_plus_1) - lgamma(n + 1) - lgamma(N_plus_1 - n); } - return -lbeta(N - n + 1, n + 1) - log(N_plus_1); + return -lbeta(N - n + 1, n + 1) - log1p(N); } } // namespace math diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 34fb5344c91..bd1564f9eb2 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -1,10 +1,81 @@ #include -#include +#include #include +#include +#include #include #include #include +TEST(MathFunctions, binomial_coefficient_log_identities) { + using stan::math::is_nan; + using stan::math::value_of; + using stan::math::var; + using stan::math::log; + using stan::math::log_sum_exp; + using stan::math::binomial_coefficient_log; + + std::vector n_to_test +// = {-0.1, 0, 1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 1e3, 1e30, 1e100}; + = {15, 1e3, 1e30, 1e100}; + + std::vector k_ratios_to_test +// = { -0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10 }; + = {1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5}; + + //Recurrence relation + for(double n_dbl : n_to_test) { + for(double k_ratio : k_ratios_to_test) { + double k_dbl = n_dbl * k_ratio; + if(n_dbl <= 0 || k_dbl <= 0) { + continue; + } + // The redundant -1 +1 is necessary as this copies the loss of precision + // for very small n_dbl + if((n_dbl - 1) + 1 - k_dbl <= 0) { + continue; + } + + var n(n_dbl); + var k(k_dbl); + var val; + + // if(n_dbl > 0 && k_dbl > 0) { + // val = log(binomial_coefficient_log(n, k)) - + // log(binomial_coefficient_log(n - 1, k) + + // binomial_coefficient_log(n - 1, k - 1)); + // } else { + // val = binomial_coefficient_log(n, k) - + // binomial_coefficient_log(n - 1, k) - + // binomial_coefficient_log(n - 1, k - 1); + // } + val = binomial_coefficient_log(n, k) / + log_sum_exp(binomial_coefficient_log(n - 1, k), + binomial_coefficient_log(n - 1, k - 1)); + + std::vector vars; + vars.push_back(n); + vars.push_back(k); + + std::vector gradients; + val.grad(vars, gradients); + + for (int i = 0; i < 2; ++i) { + EXPECT_FALSE(is_nan(gradients[i])); + } + + std::stringstream msg; + msg << std::setprecision(22) << " recurrence: n = " << n << ", k = " << k + << std::endl << "val = " << binomial_coefficient_log(n_dbl, k_dbl); + + EXPECT_NEAR(value_of(val), 1, 1e-8) << "val" << msg.str(); + EXPECT_NEAR(gradients[0], 0, 1e-8) << "dn" << msg.str(); + EXPECT_NEAR(gradients[1], 0, 1e-8) << "dx" << msg.str(); + } + } + +} + namespace binomial_coefficient_log_test_internal { struct TestValue { double n; @@ -50,146 +121,82 @@ const double NaN = stan::math::NOT_A_NUMBER; // out = StringJoin[out,"};\n"]; // out std::vector testValues = { - {-0.1, -0.9, -2.1152525390850903, -1.039918383240913, 10.708746373704937}, - {-0.1, -1.0000000000000001e-11, 1.7771062399418724e-12, - -1.922551007282891e-11, -0.17771128500984368}, - {-0.1, -1.0000000000000002e-6, 1.7770950132278696e-7, - -1.9225383585119715e-6, -0.1777077175718912}, - {-0.1, -0.001, 0.00017592768030333383, -0.0019209406823855746, - -0.17414420715602041}, - {-0.1, -0.020000000000000004, 0.0028416866729211523, -0.037823152614322564, - -0.10649980051731356}, - {-0.1, -0.05, 0.0044386492587971844, -0.09231733374172624, 0.}, - {-0.1, -0.097, 0.0005170837568402953, -0.17276563502587905, - 0.16701237954910164}, - {-0.1, -0.09999999400000001, 1.0662676433991578e-9, -0.17771127517591423, - 0.17771126364067447}, - {-0.1, -0.09999999970000001, 5.331331690149795e-11, -0.17771128455203855, - 0.17771128397527658}, - {-0.1, 0.8, -2.1152525390850907, 9.668827990464031, -10.708746373704944}, - {0.00003, -0.9, -2.21375637737528044964112, -0.933371118080918307851728, - 10.7799597405306456982503}, - {0.00003, 3.e-15, 1.48040820535563428588846e-19, - 4.9345858390686036409092e-15, 0.0000493469401735864488431421}, - {0.00003, 3.e-10, 1.48039340142161568666212e-14, - 4.93458584015035637326682e-10, 0.0000493459532446518755643083}, - {0.00003, 3.e-7, 1.46560412344434242311062e-11, - 4.93458692083243471796274e-7, 0.0000483600013795032175416619}, - {0.00003, 6.e-6, 2.36865312869367134774661e-10, - 9.86921494891292997054538e-6, 0.0000296081641072682815641386}, - {0.00003, 0.000015, 3.70102051348524044535716e-10, - 0.0000246731996398828628243033, 0}, - {0.00003, 0.0000291, 4.30798787797857754693695e-11, - 0.0000478665004969577356362445, -0.0000463861237716491741348152}, - {0.00003, 0.0000299999982, 8.88244870007509649977929e-17, - 0.0000493469372225745165598424, -0.0000493469342618230131147924}, - {0.00003, 0.00002999999991, 4.44122460318735146597007e-18, - 0.0000493469400354117708664247, -0.0000493469398873741956943572}, - {0.00003, 0.90003, -2.21375637737528044964112, 9.84658862244972739039859, - -10.7799597405306456982503}, - {0.002, -0.9, -2.21559326412971099686943, -0.931489785799666354340725, - 10.7813141298573216195769}, - {0.002, 2.e-13, 6.57013709556564677684856e-16, 3.2802775880282756733109e-13, - 0.00328506854745431610233807}, - {0.002, 2.e-8, 6.57007139476544583161173e-11, 3.28027763585130905255152e-8, - 0.00328500284665411601964384}, - {0.002, 0.00002, 6.50443564072189575550994e-8, - 0.0000328032541239782363672409, 0.00321936709346449075328701}, - {0.002, 0.0004, 1.05122171458287350859763e-6, 0.000656246880415250668008776, - 0.00197104030081914294655284}, - {0.002, 0.001, 1.64253373496215313253469e-6, 0.00164133545687894154061287, - 0}, - {0.002, 0.00194, 1.91190982195918976356429e-7, 0.00318637683127152008217264, - -0.00308796419928270612987841}, - {0.002, 0.00199999988, 3.94208202120835082737684e-13, - 0.00328506835071924267003294, -0.00328506815390258737114863}, - {0.002, 0.001999999994, 1.97104112295366725515452e-14, - 0.00328506853824172627346441, -0.00328506852840089350933798}, - {0.002, 0.902, -2.21559326412971099686943, 9.84982434405765526523621, - -10.7813141298573216195769}, - {1., -0.9, -2.85558226198351740582195, -0.459715615539276790357555, - 11.3062548910488207249193}, - {1., 1.e-10, 9.9999999988550659331851e-11, 6.44934066868432126789198e-11, - 0.99999999977101318664035}, - {1., 0.00001, 9.99988550692664559909352e-6, 6.44936087425490392714312e-6, - 0.999977101418661870834803}, - {1., 0.01, 0.00988583703486131052627978, 0.00646962905305218281638255, - 0.977200163914089454058067}, - {1., 0.2, 0.156457962917688016707705, 0.137792901804605598784786, - 0.57403132988604983615591}, + {-0.1, -0.9, -2.11525253908509081592028, -1.0399183832409129390763, 10.7087463737049383316859}, + {-0.1, -1.0000000000000001e-11, 1.77711285027681189779234e-12, -1.92253995946119474331752e-11, -0.177711285009843801688898}, + {-0.1, -2.0000000000000003e-6, 3.55415435144048059050025e-7, -3.8450735153733054435704e-6, -0.177704150099061235922007}, + {-0.1, -0.003, 0.000517083756840281579810978, -0.00575325547677734075296938, -0.167012379549101675417438}, + {-0.1, -0.020000000000000004, 0.00284168667292114343243947, -0.0378231526143224327191161, -0.106499800517313678664015}, + {-0.1, -0.05, 0.0044386492587971776182016, -0.0923173337417260614732854, 0}, + {-0.1, -0.097, 0.000517083756840282014391421, -0.172765635025879011871306, 0.167012379549101666140604}, + {-0.1, -0.09999999400000001, 1.06626764549734262955223e-9, -0.177711275175914102783629, 0.177711263640674409624395}, + {-0.1, -0.0999999998, 3.55422574122942521222662e-11, -0.17771128471653172414194, 0.177711284332023727176803}, + {-0.1, 0.8, -2.11525253908509135092963, 9.66882799046403046022088, -10.7087463737049434361164}, + {0.00003, -0.9, -2.21375637737528044964183, -0.933371118080918307851001, 10.7799597405306456982508}, + {0.00003, 3.0000000000000002e-15, 1.4804082053556344384283e-19, 4.93458583906860402434769e-15, 0.0000493469401735864500932795}, + {0.00003, 6.000000000000001e-10, 2.96075719467911309956804e-14, 9.86917168246424148279385e-10, 0.000049344966305847915513204}, + {0.00003, 9.e-7, 4.30798787797857747053402e-11, 1.48037672530856143443777e-6, 0.0000463861237716491755189368}, + {0.00003, 6.e-6, 2.36865312869367146776112e-10, 9.86921494891293022056408e-6, 0.0000296081641072682823142211}, + {0.00003, 0.000015, 3.70102051348524063287983e-10, 0.0000246731996398828634493583, 0}, + {0.00003, 0.0000291, 4.30798787797858483751685e-11, 0.0000478665004969577343409155, -0.0000463861237716491702941263}, + {0.00003, 0.000029999998200000002, 8.88244869467749046930373e-17, 0.0000493469372225745196092216, -0.0000493469342618230179633344}, + {0.00003, 0.00002999999994, 2.96081652032227041059628e-18, 0.0000493469400847597902807008, -0.0000493469399860680696581904}, + {0.00003, 0.90003, -2.2137563773752804140164, 9.84658862244972705518472, -10.7799597405306453607622}, + {0.002, -0.9, -2.21559326412971099690821, -0.931489785799666354301045, 10.7813141298573216196055}, + {0.002, 2.e-13, 6.57013709556564711297722e-16, 3.28027758802827577274611e-13, 0.00328506854745431617062256}, + {0.002, 4.e-8, 1.31400113866164680767665e-10, 6.56055536734964490768367e-8, 0.00328493714519690660826286}, + {0.002, 0.00006, 1.91190982195918985147542e-7, 0.0000984126319888139547815645, 0.0030879641992827061931754}, + {0.002, 0.0004, 1.05122171458287357370166e-6, 0.000656246880415250699426597, 0.0019710403008191429519067}, + {0.002, 0.001, 1.64253373496215320086901e-6, 0.00164133545687894157470528, 0}, + {0.002, 0.0019399999999999999, 1.91190982195919466419307e-7, 0.00318637683127151989160979, -0.00308796419928270568118357}, + {0.002, 0.00199999988, 3.94208201970328159413461e-13, 0.00328506835071924281368086, -0.00328506815390258758994027}, + {0.002, 0.001999999996, 1.31402741136587958866517e-14, 0.00328506854153159450171264, -0.0032850685349710393518526}, + {0.002, 0.902, -2.21559326412971125500408, 9.84982434405765769353491, -10.7813141298573240642835}, + {1., -0.9, -2.85558226198351740582195, -0.459715615539276790357555, 11.3062548910488207249193}, + {1., 1.e-10, 9.9999999988550662975071e-11, 6.44934066868432150285563e-11, 0.99999999977101318664035}, + {1., 0.00002, 0.0000199995420290398824268595, 0.0000128987621603843854005126, 0.999954203037316753927052}, + {1., 0.03, 0.0289783282362563119258558, 0.0195321262846958328746912, 0.932173296099201809954111}, + {1., 0.2, 0.156457962917688023080733, 0.137792901804605606966842, 0.57403132988604981390314}, {1., 0.5, 0.241564475270490444691037, 0.386294361119890618834464, 0}, - {1., 0.97, 0.028978328236256312960776, 0.951705422383897641079507, - -0.932173296099201807475084}, - {1., 0.99999994, 5.99999958782374313463811e-8, 0.999999901303960316511031, - -0.999999862607915578212576}, - {1., 0.999999997, 2.99999998969559340736596e-9, 0.999999995065197810273833, - -0.999999993130395607910641}, - {1., 1.9, -2.85558226198351740582195, 10.8465392755095439345617, - -11.3062548910488207249193}, - {8., -0.9, -4.22528965320883461943031, -0.100538838650771402252215, - 12.6649352570174581939568}, - {8., 8.e-10, 2.17428571372173153982474e-9, 9.4009611759639002250559e-11, - 2.71785714144718599267395}, - {8., 0.00008, 0.000217422931805073420417006, 9.4010053144199442029858e-6, - 2.71771615481909065227871}, - {8., 0.08, 0.211982267378255838975509, 0.00944537775278953163286055, - 2.58399543820453334937797}, - {8., 1.6, 2.90678606291134283426263, 0.208248082071609202561846, - 1.1813459431105245420685}, + {1., 0.97, 0.0289783282362563377988622, 0.951705422383897599096425, -0.932173296099201747978444}, + {1., 0.99999994, 5.99999958466560848176267e-8, 0.999999901303960368460267, -0.999999862607915650529701}, + {1., 0.999999998, 2.00000004987870302153356e-9, 0.999999996710131781531233, -0.999999995420263611904449}, + {1., 1.9, -2.85558226198351640162479, 10.846539275509534925475, -11.3062548910488116793316}, + {8., -0.9, -4.22528965320883461943031, -0.100538838650771402252215, 12.6649352570174581939568}, + {8., 8.e-10, 2.17428571372173161903874e-9, 9.40096117596390056755358e-11, 2.71785714144718599267395}, + {8., 0.00016, 0.000434834585178913876603697, 0.0000188020989077387807855132, 2.71757518207576360648536}, + {8., 0.24, 0.606274586245453630229602, 0.0286077438730426700300604, 2.35152741320850413169898}, + {8., 1.6, 2.90678606291134293918723, 0.208248082071609215411629, 1.18134594311052448766911}, {8., 4., 4.24849524204935898912334, 0.634523809523809523809524, 0}, - {8., 7.76, 0.606274586245453651115361, 2.38013515708154679097544, - -2.35152741320850411987141}, - {8., 7.99999952, 1.30457122553768403613331e-6, 2.71785635328906772378496, - -2.71785629688329908165944}, - {8., 7.999999976, 6.52285709209869625945566e-8, 2.71785710337872594517018, - -2.71785710055843758854095}, - {8., 8.9, -4.22528965320883461943031, 12.5643964183666867917046, - -12.6649352570174581939568}, - {1325., -0.9, -8.72360867216657209762532, -0.000678758619778166713080174, - 17.6139787484752694858826}, - {1325., 1.325e-7, 1.02909578020477960435539e-6, - 9.99622736492339027462876e-11, 7.76676049629224235717761}, - {1325., 0.01325, 0.102766042691370430370992, 9.99627732703511703537678e-6, - 7.74516389164454804001553}, - {1325., 13.25, 71.9898321274090629975055, 0.0100465251153381679080642, - 4.55823951796927536337622}, - {1325., 265., 659.435649329029419323398, 0.223049238391457549660772, - 1.38488037927413480754602}, - {1325., 662.5, 914.599450340845275100724, 0.692769964468769118810416, 0}, - {1325., 1285.25, 175.786260651191862665015, 3.49440932918478008636475, - -3.46396178959974953213922}, - {1325., 1324.9999205, 0.000617452276410452190170437, - 7.76662994968438931022129, -7.76662988970702332503659}, - {1325., 1324.999996025, 0.0000308728608380968862741097, - 7.76675417575202484528552, -7.76675417275315663146178}, - {1325., 1325.9, -8.72360867216657209762532, 17.6132999898554913191696, - -17.6139787484752694858826}, - {845000., -0.9, -14.5350963792733464918229, -1.06508755996070778708647e-6, - 24.0708485035538072767331}, - {845000., 0.0000845, 0.00120194816738712136581358, - 9.99999408334257028001921e-11, 14.2241691745103886002024}, - {845000., 8.45, 103.738303827736743600251, 0.0000100000440831167345507131, - 11.454909927882486225397}, - {845000., 8450., 47315.8616457576200611209, 0.0100503298765747570022521, - 4.59506127739683575088585}, - {845000., 169000., 422833.221695496506553128, 0.223143403385281321887929, - 1.3862921421877147496422}, - {845000., 422500., 585702.318235552114086514, 0.693146588844319105852931, - 0}, - {845000., 819650., 113851.158132678562120685, 3.50653876529964054840506, - -3.476079576115418785061}, - {845000., 844999.9493, 0.719108776819481762797449, - 14.1438653554064391526676, -14.1438652954064728556142}, - {845000., 844999.997465, 0.036053342347290003917417, - 14.2201459621965153509889, -14.2201459591965171216362}, - {845000., 845000.9, -14.5350963792733464918229, 24.0708474384662473160253, - -24.0708485035538072767331}, - {3.0000000000000005e15, 3., 105.120406581508329354183, NaN, NaN}, - {3.0000000000000005e15, 100., 3199.9994928023143716891, NaN, NaN}, - {3.0000000000000005e15, 12895., 350387.524360588370915867, NaN, NaN}, - {1.000000000000000000005e20, 3., 136.363346110414686040252, NaN, NaN}, - {1.000000000000000000005e20, 100., 4241.4308104325278778429, NaN, NaN}, - {1.000000000000000000005e20, 12895., 484680.092769031900296942, NaN, NaN}, + {8., 7.76, 0.606274586245454152373578, 2.38013515708154653288918, -2.35152741320850383600988}, + {8., 7.99999952, 1.30457122485101544957266e-6, 2.71785635328906813937859, -2.71785629688329952694258}, + {8., 7.999999984, 4.34857152442032476701103e-8, 2.71785711653819737865347, -2.71785711465800509058726}, + {8., 8.9, -4.22528965320883911891918, 12.5643964183667228280515, -12.664935257017494268063}, + {1325.45, -0.9, -8.72391406172695433576549, -0.000678528341300029218848988, 17.6143179550958346212948}, + {1325.45, 1.32545e-7, 1.02949027509089702892375e-6, 9.99622864543865370572321e-11, 7.76709993311726359860762}, + {1325.45, 0.026509000000000005, 0.205327156306863379978836, 0.0000199926571417065267265905, 7.72429965627466779820619}, + {1325.45, 39.7635, 175.846725556752161664879, 0.0304475435453562149043387, 3.46396589228722918392492}, + {1325.45, 265.09000000000003, 659.660660749231586923507, 0.223049270402325103613771, 1.38488085895377544787304}, + {1325.45, 662.725, 914.911196853663854616454, 0.692770092488070560492316, 0}, + {1325.45, 1285.6865, 175.846725556752235503753, 3.49441343583258486943696, -3.46396589228722863795938}, + {1325.45, 1325.449920473, 0.000617688970012744110414011, 7.76696934217533185815885, -7.76696928219795817077374}, + {1325.45, 1325.4499973491, 0.0000205898008981706567184859, 7.76709579069753161291118, -7.76709578869828579554802}, + {1325.45, 1326.3500000000001, -8.72391406172855634865104, 17.613639426763759896892, -17.6143179551050599946563}, + {845000.3, -0.9, -14.5350966987995578090701, -1.06508718182367185039981e-6, 24.0708488585827418941125}, + {845000.3, 0.00008450003000000001, 0.00120194862411218431712068, 9.99999408334467162987851e-11, 14.2241695294903594263306}, + {845000.3, 16.900006, 197.416639012626785645761, 0.0000200001881681193593499904, 10.790464758593674854536}, + {845000.3, 25350.009000000002, 113851.198555151510249396, 0.0304591891842282610860051, 3.47607957612220465494725}, + {845000.3, 169000.06000000003, 422833.371816046100941544, 0.223143403385333866823948, 1.38629214218850240580417}, + {845000.3, 422500.15, 585702.52617952879969091, 0.693146588844529182207756, 0}, + {845000.3, 819650.291, 113851.198555151775813377, 3.50653876530642990238351, -3.47607957612220154809006}, + {845000.3, 845000.2492999821, 0.71910904869572166682554, 14.1438656829570673122744, -14.1438656229571010748933}, + {845000.3, 845000.2983099994, 0.0240367434486935944287041, 14.2215320063338616170949, -14.2215320043338627454076}, + {845000.3, 845001.2000000001, -14.5350966993600009324015, 24.0708477958572381039047, -24.0708488609444199551305}, + {3.0000000000000005e15, 3.1, 108.557127724329303822723, NaN, NaN}, + {3.0000000000000005e15, 100.2, 3206.2047392970248044977, NaN, NaN}, + {3.0000000000000005e15, 12895.6, 350403.227999624864153782, NaN, NaN}, + {1.e20, 3.1, 140.841498570865873374959, NaN, NaN}, + {1.e20, 100.2, 4249.7189195624987706026, NaN, NaN}, + {1.e20, 12895.6, 484702.044995974181173534, NaN, NaN}, }; } // namespace binomial_coefficient_log_test_internal @@ -200,7 +207,6 @@ TEST(MathFunctions, binomial_coefficient_log_precomputed) { using stan::math::is_nan; using stan::math::value_of; using stan::math::var; - using stan::test::expect_near_rel; for (TestValue t : testValues) { std::stringstream msg; @@ -221,21 +227,20 @@ TEST(MathFunctions, binomial_coefficient_log_precomputed) { EXPECT_FALSE(is_nan(gradients[i])); } - expect_near_rel(msg.str(), value_of(val), t.val); + double tol_val = std::max(1e-14 * fabs(t.val), 1e-14); + EXPECT_NEAR(value_of(val), t.val, tol_val) << msg.str(); - double tol_grad; + std::function tol_grad; if (n < 1 || k < 1) { - tol_grad = 1e-7; + tol_grad = [](double x) { return std::max(fabs(x) * 1e-8, 1e-7); }; } else { - tol_grad = 1e-8; + tol_grad = [](double x) { return std::max(fabs(x) * 1e-10, 1e-8); }; } if (!is_nan(t.dn)) { - expect_near_rel(std::string("dn: ") + msg.str(), gradients[0], t.dn, - tol_grad); + EXPECT_NEAR(gradients[0], t.dn, tol_grad(t.dn)) << "dn: " << msg.str(); } if (!is_nan(t.dk)) { - expect_near_rel(std::string("dk: ") + msg.str(), gradients[1], t.dk, - tol_grad); + EXPECT_NEAR(gradients[1], t.dk, tol_grad(t.dk)) << "dk: " << msg.str(); } } } From cd3381a3abb76fdcc30a4d5278d5e8b71ca419d8 Mon Sep 17 00:00:00 2001 From: martinmodrak Date: Wed, 29 Jan 2020 15:59:45 +0100 Subject: [PATCH 15/26] Work towards formula tests (failing) --- .../rev/fun/binomial_coefficient_log_test.cpp | 15 +-- test/unit/math/rev/fun/lbeta_test.cpp | 94 +++++++++++++++++++ 2 files changed, 97 insertions(+), 12 deletions(-) diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index bd1564f9eb2..2ce6b471176 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -40,18 +40,8 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { var k(k_dbl); var val; - // if(n_dbl > 0 && k_dbl > 0) { - // val = log(binomial_coefficient_log(n, k)) - - // log(binomial_coefficient_log(n - 1, k) + - // binomial_coefficient_log(n - 1, k - 1)); - // } else { - // val = binomial_coefficient_log(n, k) - - // binomial_coefficient_log(n - 1, k) - - // binomial_coefficient_log(n - 1, k - 1); - // } val = binomial_coefficient_log(n, k) / - log_sum_exp(binomial_coefficient_log(n - 1, k), - binomial_coefficient_log(n - 1, k - 1)); + (binomial_coefficient_log(n - 1, k - 1) + log(n) - log(k)) ; std::vector vars; vars.push_back(n); @@ -65,7 +55,8 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { } std::stringstream msg; - msg << std::setprecision(22) << " recurrence: n = " << n << ", k = " << k + msg << std::setprecision(22) << " (n - 1) choose (k - 1): n = " << n + << ", k = " << k << std::endl << "val = " << binomial_coefficient_log(n_dbl, k_dbl); EXPECT_NEAR(value_of(val), 1, 1e-8) << "val" << msg.str(); diff --git a/test/unit/math/rev/fun/lbeta_test.cpp b/test/unit/math/rev/fun/lbeta_test.cpp index 780cd8186f2..05d08c1d479 100644 --- a/test/unit/math/rev/fun/lbeta_test.cpp +++ b/test/unit/math/rev/fun/lbeta_test.cpp @@ -8,6 +8,100 @@ #include #include +struct relative_tolerance { + double tol; + double tol_min; + + double operator()(const double x) const { + return std::max(tol * fabs(x), tol_min); + } + + double operator()(const double x, const double y) const { + return std::max(tol * 0.5 * (fabs(x) + fabs(y)), tol_min); + } +}; + +struct identity_tolerances { + relative_tolerance value; + relative_tolerance gradient; +}; + +template void expect_identity( + const std::string& msg, const identity_tolerances& tolerances, + const F1 lh, const F2 rh, double x_dbl, double y_dbl) { + using stan::math::var; + + stan::math::start_nested(); + + var x(x_dbl); + var y(y_dbl); + + std::vector vars = {x, y}; + + var left = lh(x, y); + double left_dbl = value_of(left); + std::vector gradients_left; + left.grad(vars, gradients_left); + + stan::math::set_zero_all_adjoints_nested(); + + var right = rh(x, y); + double right_dbl = value_of(right); + std::vector gradients_right; + right.grad(vars, gradients_right); + + std::stringstream args; + args << std::setprecision(22) << "args = [" << x << "," << y << "]"; + double tol_val = tolerances.value(left_dbl, right_dbl); + EXPECT_NEAR(left_dbl, right_dbl, tol_val) << "value, " << args.str() << ": " << msg; + + for(size_t i = 0; i < gradients_left.size(); ++i) { + double tol_grad = tolerances.gradient(gradients_left[i], gradients_right[i]); + EXPECT_NEAR(gradients_left[i], gradients_right[i], tol_grad) << "grad_" << i << ", " << args.str() << ": " << msg; + } + + stan::math::recover_memory_nested(); +} + +TEST(MathFunctions, lbeta_identities_gradient) { + using stan::math::lbeta; + using stan::math::pi; + using stan::math::var; + + std::vector to_test + = {1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 1e3, 1e30, 1e100}; + + identity_tolerances tol { {1e-15, 1e-15}, {1e-10,1e-10}}; + + for (double x : to_test) { + for (double y : to_test) { + auto rh = [](const var& a, const var& b) { + return stan::math::log_sum_exp(lbeta(a + 1, b), lbeta(a, b + 1)); + }; + expect_identity("succesors", tol, static_cast(lbeta), rh, x, y); + } + } + + for (double x : to_test) { + if (x < 1) { + std::stringstream msg; + msg << std::setprecision(22) << "sin: x = " << x; + double lh = lbeta(x, 1.0 - x); + double rh = log(pi()) - log(sin(pi() * x)); + EXPECT_NEAR(lh, rh, tol.value(lh, rh)) << msg.str(); + } + } + + for (double x : to_test) { + std::stringstream msg; + msg << std::setprecision(22) << "inv: x = " << x; + double lh = lbeta(x, 1.0); + double rh = -log(x); + EXPECT_NEAR(lh, rh, tol.value(lh, rh)) << msg.str(); + } +} + + namespace lbeta_test_internal { struct TestValue { double x; From d014464dc72e9a4cdb98b1b28ca83d47e7424ba7 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 29 Jan 2020 10:00:52 -0500 Subject: [PATCH 16/26] [Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/stable/2017-11-14) --- .../rev/fun/binomial_coefficient_log_test.cpp | 235 +++++++++++------- test/unit/math/rev/fun/lbeta_test.cpp | 27 +- 2 files changed, 167 insertions(+), 95 deletions(-) diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 2ce6b471176..00773da66d8 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -8,40 +8,41 @@ #include TEST(MathFunctions, binomial_coefficient_log_identities) { + using stan::math::binomial_coefficient_log; using stan::math::is_nan; - using stan::math::value_of; - using stan::math::var; using stan::math::log; using stan::math::log_sum_exp; - using stan::math::binomial_coefficient_log; + using stan::math::value_of; + using stan::math::var; std::vector n_to_test -// = {-0.1, 0, 1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 1e3, 1e30, 1e100}; + // = {-0.1, 0, 1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 1e3, 1e30, 1e100}; = {15, 1e3, 1e30, 1e100}; std::vector k_ratios_to_test -// = { -0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10 }; + // = { -0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10 + // }; = {1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5}; - //Recurrence relation - for(double n_dbl : n_to_test) { - for(double k_ratio : k_ratios_to_test) { + // Recurrence relation + for (double n_dbl : n_to_test) { + for (double k_ratio : k_ratios_to_test) { double k_dbl = n_dbl * k_ratio; - if(n_dbl <= 0 || k_dbl <= 0) { + if (n_dbl <= 0 || k_dbl <= 0) { continue; } // The redundant -1 +1 is necessary as this copies the loss of precision // for very small n_dbl - if((n_dbl - 1) + 1 - k_dbl <= 0) { + if ((n_dbl - 1) + 1 - k_dbl <= 0) { continue; } var n(n_dbl); var k(k_dbl); var val; - - val = binomial_coefficient_log(n, k) / - (binomial_coefficient_log(n - 1, k - 1) + log(n) - log(k)) ; + + val = binomial_coefficient_log(n, k) + / (binomial_coefficient_log(n - 1, k - 1) + log(n) - log(k)); std::vector vars; vars.push_back(n); @@ -55,16 +56,15 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { } std::stringstream msg; - msg << std::setprecision(22) << " (n - 1) choose (k - 1): n = " << n - << ", k = " << k - << std::endl << "val = " << binomial_coefficient_log(n_dbl, k_dbl); + msg << std::setprecision(22) << " (n - 1) choose (k - 1): n = " << n + << ", k = " << k << std::endl + << "val = " << binomial_coefficient_log(n_dbl, k_dbl); EXPECT_NEAR(value_of(val), 1, 1e-8) << "val" << msg.str(); EXPECT_NEAR(gradients[0], 0, 1e-8) << "dn" << msg.str(); EXPECT_NEAR(gradients[1], 0, 1e-8) << "dx" << msg.str(); } } - } namespace binomial_coefficient_log_test_internal { @@ -112,76 +112,143 @@ const double NaN = stan::math::NOT_A_NUMBER; // out = StringJoin[out,"};\n"]; // out std::vector testValues = { - {-0.1, -0.9, -2.11525253908509081592028, -1.0399183832409129390763, 10.7087463737049383316859}, - {-0.1, -1.0000000000000001e-11, 1.77711285027681189779234e-12, -1.92253995946119474331752e-11, -0.177711285009843801688898}, - {-0.1, -2.0000000000000003e-6, 3.55415435144048059050025e-7, -3.8450735153733054435704e-6, -0.177704150099061235922007}, - {-0.1, -0.003, 0.000517083756840281579810978, -0.00575325547677734075296938, -0.167012379549101675417438}, - {-0.1, -0.020000000000000004, 0.00284168667292114343243947, -0.0378231526143224327191161, -0.106499800517313678664015}, + {-0.1, -0.9, -2.11525253908509081592028, -1.0399183832409129390763, + 10.7087463737049383316859}, + {-0.1, -1.0000000000000001e-11, 1.77711285027681189779234e-12, + -1.92253995946119474331752e-11, -0.177711285009843801688898}, + {-0.1, -2.0000000000000003e-6, 3.55415435144048059050025e-7, + -3.8450735153733054435704e-6, -0.177704150099061235922007}, + {-0.1, -0.003, 0.000517083756840281579810978, -0.00575325547677734075296938, + -0.167012379549101675417438}, + {-0.1, -0.020000000000000004, 0.00284168667292114343243947, + -0.0378231526143224327191161, -0.106499800517313678664015}, {-0.1, -0.05, 0.0044386492587971776182016, -0.0923173337417260614732854, 0}, - {-0.1, -0.097, 0.000517083756840282014391421, -0.172765635025879011871306, 0.167012379549101666140604}, - {-0.1, -0.09999999400000001, 1.06626764549734262955223e-9, -0.177711275175914102783629, 0.177711263640674409624395}, - {-0.1, -0.0999999998, 3.55422574122942521222662e-11, -0.17771128471653172414194, 0.177711284332023727176803}, - {-0.1, 0.8, -2.11525253908509135092963, 9.66882799046403046022088, -10.7087463737049434361164}, - {0.00003, -0.9, -2.21375637737528044964183, -0.933371118080918307851001, 10.7799597405306456982508}, - {0.00003, 3.0000000000000002e-15, 1.4804082053556344384283e-19, 4.93458583906860402434769e-15, 0.0000493469401735864500932795}, - {0.00003, 6.000000000000001e-10, 2.96075719467911309956804e-14, 9.86917168246424148279385e-10, 0.000049344966305847915513204}, - {0.00003, 9.e-7, 4.30798787797857747053402e-11, 1.48037672530856143443777e-6, 0.0000463861237716491755189368}, - {0.00003, 6.e-6, 2.36865312869367146776112e-10, 9.86921494891293022056408e-6, 0.0000296081641072682823142211}, - {0.00003, 0.000015, 3.70102051348524063287983e-10, 0.0000246731996398828634493583, 0}, - {0.00003, 0.0000291, 4.30798787797858483751685e-11, 0.0000478665004969577343409155, -0.0000463861237716491702941263}, - {0.00003, 0.000029999998200000002, 8.88244869467749046930373e-17, 0.0000493469372225745196092216, -0.0000493469342618230179633344}, - {0.00003, 0.00002999999994, 2.96081652032227041059628e-18, 0.0000493469400847597902807008, -0.0000493469399860680696581904}, - {0.00003, 0.90003, -2.2137563773752804140164, 9.84658862244972705518472, -10.7799597405306453607622}, - {0.002, -0.9, -2.21559326412971099690821, -0.931489785799666354301045, 10.7813141298573216196055}, - {0.002, 2.e-13, 6.57013709556564711297722e-16, 3.28027758802827577274611e-13, 0.00328506854745431617062256}, - {0.002, 4.e-8, 1.31400113866164680767665e-10, 6.56055536734964490768367e-8, 0.00328493714519690660826286}, - {0.002, 0.00006, 1.91190982195918985147542e-7, 0.0000984126319888139547815645, 0.0030879641992827061931754}, - {0.002, 0.0004, 1.05122171458287357370166e-6, 0.000656246880415250699426597, 0.0019710403008191429519067}, - {0.002, 0.001, 1.64253373496215320086901e-6, 0.00164133545687894157470528, 0}, - {0.002, 0.0019399999999999999, 1.91190982195919466419307e-7, 0.00318637683127151989160979, -0.00308796419928270568118357}, - {0.002, 0.00199999988, 3.94208201970328159413461e-13, 0.00328506835071924281368086, -0.00328506815390258758994027}, - {0.002, 0.001999999996, 1.31402741136587958866517e-14, 0.00328506854153159450171264, -0.0032850685349710393518526}, - {0.002, 0.902, -2.21559326412971125500408, 9.84982434405765769353491, -10.7813141298573240642835}, - {1., -0.9, -2.85558226198351740582195, -0.459715615539276790357555, 11.3062548910488207249193}, - {1., 1.e-10, 9.9999999988550662975071e-11, 6.44934066868432150285563e-11, 0.99999999977101318664035}, - {1., 0.00002, 0.0000199995420290398824268595, 0.0000128987621603843854005126, 0.999954203037316753927052}, - {1., 0.03, 0.0289783282362563119258558, 0.0195321262846958328746912, 0.932173296099201809954111}, - {1., 0.2, 0.156457962917688023080733, 0.137792901804605606966842, 0.57403132988604981390314}, + {-0.1, -0.097, 0.000517083756840282014391421, -0.172765635025879011871306, + 0.167012379549101666140604}, + {-0.1, -0.09999999400000001, 1.06626764549734262955223e-9, + -0.177711275175914102783629, 0.177711263640674409624395}, + {-0.1, -0.0999999998, 3.55422574122942521222662e-11, + -0.17771128471653172414194, 0.177711284332023727176803}, + {-0.1, 0.8, -2.11525253908509135092963, 9.66882799046403046022088, + -10.7087463737049434361164}, + {0.00003, -0.9, -2.21375637737528044964183, -0.933371118080918307851001, + 10.7799597405306456982508}, + {0.00003, 3.0000000000000002e-15, 1.4804082053556344384283e-19, + 4.93458583906860402434769e-15, 0.0000493469401735864500932795}, + {0.00003, 6.000000000000001e-10, 2.96075719467911309956804e-14, + 9.86917168246424148279385e-10, 0.000049344966305847915513204}, + {0.00003, 9.e-7, 4.30798787797857747053402e-11, + 1.48037672530856143443777e-6, 0.0000463861237716491755189368}, + {0.00003, 6.e-6, 2.36865312869367146776112e-10, + 9.86921494891293022056408e-6, 0.0000296081641072682823142211}, + {0.00003, 0.000015, 3.70102051348524063287983e-10, + 0.0000246731996398828634493583, 0}, + {0.00003, 0.0000291, 4.30798787797858483751685e-11, + 0.0000478665004969577343409155, -0.0000463861237716491702941263}, + {0.00003, 0.000029999998200000002, 8.88244869467749046930373e-17, + 0.0000493469372225745196092216, -0.0000493469342618230179633344}, + {0.00003, 0.00002999999994, 2.96081652032227041059628e-18, + 0.0000493469400847597902807008, -0.0000493469399860680696581904}, + {0.00003, 0.90003, -2.2137563773752804140164, 9.84658862244972705518472, + -10.7799597405306453607622}, + {0.002, -0.9, -2.21559326412971099690821, -0.931489785799666354301045, + 10.7813141298573216196055}, + {0.002, 2.e-13, 6.57013709556564711297722e-16, + 3.28027758802827577274611e-13, 0.00328506854745431617062256}, + {0.002, 4.e-8, 1.31400113866164680767665e-10, 6.56055536734964490768367e-8, + 0.00328493714519690660826286}, + {0.002, 0.00006, 1.91190982195918985147542e-7, + 0.0000984126319888139547815645, 0.0030879641992827061931754}, + {0.002, 0.0004, 1.05122171458287357370166e-6, 0.000656246880415250699426597, + 0.0019710403008191429519067}, + {0.002, 0.001, 1.64253373496215320086901e-6, 0.00164133545687894157470528, + 0}, + {0.002, 0.0019399999999999999, 1.91190982195919466419307e-7, + 0.00318637683127151989160979, -0.00308796419928270568118357}, + {0.002, 0.00199999988, 3.94208201970328159413461e-13, + 0.00328506835071924281368086, -0.00328506815390258758994027}, + {0.002, 0.001999999996, 1.31402741136587958866517e-14, + 0.00328506854153159450171264, -0.0032850685349710393518526}, + {0.002, 0.902, -2.21559326412971125500408, 9.84982434405765769353491, + -10.7813141298573240642835}, + {1., -0.9, -2.85558226198351740582195, -0.459715615539276790357555, + 11.3062548910488207249193}, + {1., 1.e-10, 9.9999999988550662975071e-11, 6.44934066868432150285563e-11, + 0.99999999977101318664035}, + {1., 0.00002, 0.0000199995420290398824268595, + 0.0000128987621603843854005126, 0.999954203037316753927052}, + {1., 0.03, 0.0289783282362563119258558, 0.0195321262846958328746912, + 0.932173296099201809954111}, + {1., 0.2, 0.156457962917688023080733, 0.137792901804605606966842, + 0.57403132988604981390314}, {1., 0.5, 0.241564475270490444691037, 0.386294361119890618834464, 0}, - {1., 0.97, 0.0289783282362563377988622, 0.951705422383897599096425, -0.932173296099201747978444}, - {1., 0.99999994, 5.99999958466560848176267e-8, 0.999999901303960368460267, -0.999999862607915650529701}, - {1., 0.999999998, 2.00000004987870302153356e-9, 0.999999996710131781531233, -0.999999995420263611904449}, - {1., 1.9, -2.85558226198351640162479, 10.846539275509534925475, -11.3062548910488116793316}, - {8., -0.9, -4.22528965320883461943031, -0.100538838650771402252215, 12.6649352570174581939568}, - {8., 8.e-10, 2.17428571372173161903874e-9, 9.40096117596390056755358e-11, 2.71785714144718599267395}, - {8., 0.00016, 0.000434834585178913876603697, 0.0000188020989077387807855132, 2.71757518207576360648536}, - {8., 0.24, 0.606274586245453630229602, 0.0286077438730426700300604, 2.35152741320850413169898}, - {8., 1.6, 2.90678606291134293918723, 0.208248082071609215411629, 1.18134594311052448766911}, + {1., 0.97, 0.0289783282362563377988622, 0.951705422383897599096425, + -0.932173296099201747978444}, + {1., 0.99999994, 5.99999958466560848176267e-8, 0.999999901303960368460267, + -0.999999862607915650529701}, + {1., 0.999999998, 2.00000004987870302153356e-9, 0.999999996710131781531233, + -0.999999995420263611904449}, + {1., 1.9, -2.85558226198351640162479, 10.846539275509534925475, + -11.3062548910488116793316}, + {8., -0.9, -4.22528965320883461943031, -0.100538838650771402252215, + 12.6649352570174581939568}, + {8., 8.e-10, 2.17428571372173161903874e-9, 9.40096117596390056755358e-11, + 2.71785714144718599267395}, + {8., 0.00016, 0.000434834585178913876603697, 0.0000188020989077387807855132, + 2.71757518207576360648536}, + {8., 0.24, 0.606274586245453630229602, 0.0286077438730426700300604, + 2.35152741320850413169898}, + {8., 1.6, 2.90678606291134293918723, 0.208248082071609215411629, + 1.18134594311052448766911}, {8., 4., 4.24849524204935898912334, 0.634523809523809523809524, 0}, - {8., 7.76, 0.606274586245454152373578, 2.38013515708154653288918, -2.35152741320850383600988}, - {8., 7.99999952, 1.30457122485101544957266e-6, 2.71785635328906813937859, -2.71785629688329952694258}, - {8., 7.999999984, 4.34857152442032476701103e-8, 2.71785711653819737865347, -2.71785711465800509058726}, - {8., 8.9, -4.22528965320883911891918, 12.5643964183667228280515, -12.664935257017494268063}, - {1325.45, -0.9, -8.72391406172695433576549, -0.000678528341300029218848988, 17.6143179550958346212948}, - {1325.45, 1.32545e-7, 1.02949027509089702892375e-6, 9.99622864543865370572321e-11, 7.76709993311726359860762}, - {1325.45, 0.026509000000000005, 0.205327156306863379978836, 0.0000199926571417065267265905, 7.72429965627466779820619}, - {1325.45, 39.7635, 175.846725556752161664879, 0.0304475435453562149043387, 3.46396589228722918392492}, - {1325.45, 265.09000000000003, 659.660660749231586923507, 0.223049270402325103613771, 1.38488085895377544787304}, - {1325.45, 662.725, 914.911196853663854616454, 0.692770092488070560492316, 0}, - {1325.45, 1285.6865, 175.846725556752235503753, 3.49441343583258486943696, -3.46396589228722863795938}, - {1325.45, 1325.449920473, 0.000617688970012744110414011, 7.76696934217533185815885, -7.76696928219795817077374}, - {1325.45, 1325.4499973491, 0.0000205898008981706567184859, 7.76709579069753161291118, -7.76709578869828579554802}, - {1325.45, 1326.3500000000001, -8.72391406172855634865104, 17.613639426763759896892, -17.6143179551050599946563}, - {845000.3, -0.9, -14.5350966987995578090701, -1.06508718182367185039981e-6, 24.0708488585827418941125}, - {845000.3, 0.00008450003000000001, 0.00120194862411218431712068, 9.99999408334467162987851e-11, 14.2241695294903594263306}, - {845000.3, 16.900006, 197.416639012626785645761, 0.0000200001881681193593499904, 10.790464758593674854536}, - {845000.3, 25350.009000000002, 113851.198555151510249396, 0.0304591891842282610860051, 3.47607957612220465494725}, - {845000.3, 169000.06000000003, 422833.371816046100941544, 0.223143403385333866823948, 1.38629214218850240580417}, - {845000.3, 422500.15, 585702.52617952879969091, 0.693146588844529182207756, 0}, - {845000.3, 819650.291, 113851.198555151775813377, 3.50653876530642990238351, -3.47607957612220154809006}, - {845000.3, 845000.2492999821, 0.71910904869572166682554, 14.1438656829570673122744, -14.1438656229571010748933}, - {845000.3, 845000.2983099994, 0.0240367434486935944287041, 14.2215320063338616170949, -14.2215320043338627454076}, - {845000.3, 845001.2000000001, -14.5350966993600009324015, 24.0708477958572381039047, -24.0708488609444199551305}, + {8., 7.76, 0.606274586245454152373578, 2.38013515708154653288918, + -2.35152741320850383600988}, + {8., 7.99999952, 1.30457122485101544957266e-6, 2.71785635328906813937859, + -2.71785629688329952694258}, + {8., 7.999999984, 4.34857152442032476701103e-8, 2.71785711653819737865347, + -2.71785711465800509058726}, + {8., 8.9, -4.22528965320883911891918, 12.5643964183667228280515, + -12.664935257017494268063}, + {1325.45, -0.9, -8.72391406172695433576549, -0.000678528341300029218848988, + 17.6143179550958346212948}, + {1325.45, 1.32545e-7, 1.02949027509089702892375e-6, + 9.99622864543865370572321e-11, 7.76709993311726359860762}, + {1325.45, 0.026509000000000005, 0.205327156306863379978836, + 0.0000199926571417065267265905, 7.72429965627466779820619}, + {1325.45, 39.7635, 175.846725556752161664879, 0.0304475435453562149043387, + 3.46396589228722918392492}, + {1325.45, 265.09000000000003, 659.660660749231586923507, + 0.223049270402325103613771, 1.38488085895377544787304}, + {1325.45, 662.725, 914.911196853663854616454, 0.692770092488070560492316, + 0}, + {1325.45, 1285.6865, 175.846725556752235503753, 3.49441343583258486943696, + -3.46396589228722863795938}, + {1325.45, 1325.449920473, 0.000617688970012744110414011, + 7.76696934217533185815885, -7.76696928219795817077374}, + {1325.45, 1325.4499973491, 0.0000205898008981706567184859, + 7.76709579069753161291118, -7.76709578869828579554802}, + {1325.45, 1326.3500000000001, -8.72391406172855634865104, + 17.613639426763759896892, -17.6143179551050599946563}, + {845000.3, -0.9, -14.5350966987995578090701, -1.06508718182367185039981e-6, + 24.0708488585827418941125}, + {845000.3, 0.00008450003000000001, 0.00120194862411218431712068, + 9.99999408334467162987851e-11, 14.2241695294903594263306}, + {845000.3, 16.900006, 197.416639012626785645761, + 0.0000200001881681193593499904, 10.790464758593674854536}, + {845000.3, 25350.009000000002, 113851.198555151510249396, + 0.0304591891842282610860051, 3.47607957612220465494725}, + {845000.3, 169000.06000000003, 422833.371816046100941544, + 0.223143403385333866823948, 1.38629214218850240580417}, + {845000.3, 422500.15, 585702.52617952879969091, 0.693146588844529182207756, + 0}, + {845000.3, 819650.291, 113851.198555151775813377, 3.50653876530642990238351, + -3.47607957612220154809006}, + {845000.3, 845000.2492999821, 0.71910904869572166682554, + 14.1438656829570673122744, -14.1438656229571010748933}, + {845000.3, 845000.2983099994, 0.0240367434486935944287041, + 14.2215320063338616170949, -14.2215320043338627454076}, + {845000.3, 845001.2000000001, -14.5350966993600009324015, + 24.0708477958572381039047, -24.0708488609444199551305}, {3.0000000000000005e15, 3.1, 108.557127724329303822723, NaN, NaN}, {3.0000000000000005e15, 100.2, 3206.2047392970248044977, NaN, NaN}, {3.0000000000000005e15, 12895.6, 350403.227999624864153782, NaN, NaN}, diff --git a/test/unit/math/rev/fun/lbeta_test.cpp b/test/unit/math/rev/fun/lbeta_test.cpp index 05d08c1d479..e02d77a038a 100644 --- a/test/unit/math/rev/fun/lbeta_test.cpp +++ b/test/unit/math/rev/fun/lbeta_test.cpp @@ -26,11 +26,12 @@ struct identity_tolerances { relative_tolerance gradient; }; -template void expect_identity( - const std::string& msg, const identity_tolerances& tolerances, - const F1 lh, const F2 rh, double x_dbl, double y_dbl) { +template +void expect_identity(const std::string& msg, + const identity_tolerances& tolerances, const F1 lh, + const F2 rh, double x_dbl, double y_dbl) { using stan::math::var; - + stan::math::start_nested(); var x(x_dbl); @@ -53,11 +54,14 @@ template void expect_identity( std::stringstream args; args << std::setprecision(22) << "args = [" << x << "," << y << "]"; double tol_val = tolerances.value(left_dbl, right_dbl); - EXPECT_NEAR(left_dbl, right_dbl, tol_val) << "value, " << args.str() << ": " << msg; + EXPECT_NEAR(left_dbl, right_dbl, tol_val) + << "value, " << args.str() << ": " << msg; - for(size_t i = 0; i < gradients_left.size(); ++i) { - double tol_grad = tolerances.gradient(gradients_left[i], gradients_right[i]); - EXPECT_NEAR(gradients_left[i], gradients_right[i], tol_grad) << "grad_" << i << ", " << args.str() << ": " << msg; + for (size_t i = 0; i < gradients_left.size(); ++i) { + double tol_grad + = tolerances.gradient(gradients_left[i], gradients_right[i]); + EXPECT_NEAR(gradients_left[i], gradients_right[i], tol_grad) + << "grad_" << i << ", " << args.str() << ": " << msg; } stan::math::recover_memory_nested(); @@ -71,14 +75,16 @@ TEST(MathFunctions, lbeta_identities_gradient) { std::vector to_test = {1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 1e3, 1e30, 1e100}; - identity_tolerances tol { {1e-15, 1e-15}, {1e-10,1e-10}}; + identity_tolerances tol{{1e-15, 1e-15}, {1e-10, 1e-10}}; for (double x : to_test) { for (double y : to_test) { auto rh = [](const var& a, const var& b) { return stan::math::log_sum_exp(lbeta(a + 1, b), lbeta(a, b + 1)); }; - expect_identity("succesors", tol, static_cast(lbeta), rh, x, y); + expect_identity("succesors", tol, + static_cast(lbeta), rh, + x, y); } } @@ -101,7 +107,6 @@ TEST(MathFunctions, lbeta_identities_gradient) { } } - namespace lbeta_test_internal { struct TestValue { double x; From 51ef4911d768831f739d88ecc23d1e26bc2720f2 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Sat, 7 Mar 2020 20:57:16 +0100 Subject: [PATCH 17/26] Passing identity tests, ad test --- .../prim/fun/binomial_coefficient_log.hpp | 7 +- stan/math/rev/fun.hpp | 1 + .../math/rev/fun/binomial_coefficient_log.hpp | 85 +++++++++++++++++ .../rev/fun/binomial_coefficient_log_test.cpp | 81 +++++++++++----- test/unit/math/rev/fun/lbeta_test.cpp | 95 +++++++++---------- 5 files changed, 192 insertions(+), 77 deletions(-) create mode 100644 stan/math/rev/fun/binomial_coefficient_log.hpp diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 54904bb3ac7..938ee6dc42a 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -62,12 +62,11 @@ namespace math { \end{cases} \f] * - * @tparam T_N type of the first argument - * @tparam T_n type of the second argument * This function is numerically more stable than naive evaluation via lgamma. * - * @tparam T_N type of N. - * @tparam T_n type of n. + * @tparam T_N type of the first argument + * @tparam T_n type of the second argument + * * @param N total number of objects. * @param n number of objects chosen. * @return log (N choose n). diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 1deb92f9738..d692cdf0513 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/binomial_coefficient_log.hpp b/stan/math/rev/fun/binomial_coefficient_log.hpp new file mode 100644 index 00000000000..a7cadfbe58c --- /dev/null +++ b/stan/math/rev/fun/binomial_coefficient_log.hpp @@ -0,0 +1,85 @@ +#ifndef STAN_MATH_REV_FUN_BINOMIAL_COEFFICIENT_LOG_HPP +#define STAN_MATH_REV_FUN_BINOMIAL_COEFFICIENT_LOG_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +namespace internal { +class binomial_coefficient_log_vv_vari : public op_vv_vari { + public: + binomial_coefficient_log_vv_vari(vari* avi, vari* bvi) + : op_vv_vari(binomial_coefficient_log(avi->val_, bvi->val_), avi, bvi) {} + void chain() { + double digamma_ambp1 = digamma(avi_->val_ - bvi_->val_ + 1); + + avi_->adj_ += adj_ * (digamma(avi_->val_ + 1) - digamma_ambp1); + bvi_->adj_ += adj_ * (digamma_ambp1 - digamma(bvi_->val_ + 1)); + } +}; + +class binomial_coefficient_log_vd_vari : public op_vd_vari { + public: + binomial_coefficient_log_vd_vari(vari* avi, double b) + : op_vd_vari(binomial_coefficient_log(avi->val_, b), avi, b) {} + void chain() { + avi_->adj_ += adj_ * (digamma(avi_->val_ + 1) - digamma(avi_->val_ - bd_ + 1)); + } +}; + +class binomial_coefficient_log_dv_vari : public op_dv_vari { + public: + binomial_coefficient_log_dv_vari(double a, vari* bvi) + : op_dv_vari(binomial_coefficient_log(a, bvi->val_), a, bvi) {} + void chain() { + bvi_->adj_ += adj_ * (digamma(ad_ - bvi_->val_ + 1) - digamma(bvi_->val_ + 1)); + } +}; +} // namespace internal + +/** + * Return the log of the binomial coefficient for the specified + * arguments and its gradients. + * + * See the docs for the prim version for all relevant formulae. + * @param a var Argument + * @param b var Argument + * @return Result of log (a choose b) + */ +inline var binomial_coefficient_log(const var& a, const var& b) { + return var(new internal::binomial_coefficient_log_vv_vari(a.vi_, b.vi_)); +} + +/** + * Return the log of the binomial coefficient for the specified + * arguments and its gradients. + * + * See the docs for the prim version for all relevant formulae. + * @param a var Argument + * @param b double Argument + * @return Result of log (a choose b) + */ +inline var binomial_coefficient_log(const var& a, double b) { + return var(new internal::binomial_coefficient_log_vd_vari(a.vi_, b)); +} + +/** + * Return the log of the binomial coefficient for the specified + * arguments and its gradients. + * + * See the docs for the prim version for all relevant formulae. + * @param a double Argument + * @param b var Argument + * @return Result of log (a choose b) + */ +inline var binomial_coefficient_log(double a, const var& b) { + return var(new internal::binomial_coefficient_log_dv_vari(a, b.vi_)); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 00773da66d8..7a7179ce4cd 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include #include #include #include @@ -14,15 +16,14 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { using stan::math::log_sum_exp; using stan::math::value_of; using stan::math::var; + using stan::test::expect_near_rel; std::vector n_to_test - // = {-0.1, 0, 1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 1e3, 1e30, 1e100}; - = {15, 1e3, 1e30, 1e100}; + = {-0.1, 0, 1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 15, 10, 1e3, 1e30, 1e100}; std::vector k_ratios_to_test - // = { -0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10 - // }; - = {1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5}; + = { -0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10 + }; // Recurrence relation for (double n_dbl : n_to_test) { @@ -37,36 +38,69 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { continue; } + stan::math::nested_rev_autodiff nested; var n(n_dbl); var k(k_dbl); - var val; - val = binomial_coefficient_log(n, k) - / (binomial_coefficient_log(n - 1, k - 1) + log(n) - log(k)); + // TODO(martinmodrak) Use the framework for testing identities, once it is ready + var val_left = binomial_coefficient_log(n, k); + var val_right_partial; + var val_right; + // Choose the more stable identity + if(n_dbl > 1 && k_dbl > 1 && (n_dbl - 1) + 1 - k_dbl > 0 ) { + val_right_partial = binomial_coefficient_log(n - 1, k - 1); + val_right = val_right_partial + log(n) - log(k); + } else { + val_right_partial = binomial_coefficient_log(n + 1, k + 1); + val_right = val_right_partial - log(n + 1) + log(k + 1); + } std::vector vars; vars.push_back(n); vars.push_back(k); - std::vector gradients; - val.grad(vars, gradients); + std::vector gradients_left; + val_left.grad(vars, gradients_left); + + nested.set_zero_all_adjoints(); + + std::vector gradients_right; + val_right.grad(vars, gradients_right); for (int i = 0; i < 2; ++i) { - EXPECT_FALSE(is_nan(gradients[i])); + EXPECT_FALSE(is_nan(gradients_left[i])); + EXPECT_FALSE(is_nan(gradients_right[i])); } std::stringstream msg; - msg << std::setprecision(22) << " (n - 1) choose (k - 1): n = " << n + msg << std::setprecision(22) << " successor: n = " << n << ", k = " << k << std::endl - << "val = " << binomial_coefficient_log(n_dbl, k_dbl); + << "val = " << val_left + << ", val2 = " << val_right_partial << std::endl + << ", logn = " << log(n) + << ", logk = " << log(k); + - EXPECT_NEAR(value_of(val), 1, 1e-8) << "val" << msg.str(); - EXPECT_NEAR(gradients[0], 0, 1e-8) << "dn" << msg.str(); - EXPECT_NEAR(gradients[1], 0, 1e-8) << "dx" << msg.str(); + expect_near_rel(std::string("val") + msg.str(), value_of(val_left), value_of(val_right)); + expect_near_rel(std::string("dn") + msg.str(), gradients_left[0], gradients_right[0]); + expect_near_rel(std::string("dk") + msg.str(), gradients_left[1], gradients_right[1]); } } } +TEST(MathFunctions, binomial_coefficient_log_ad) { + using stan::test::expect_ad; + + auto f = [](const auto& n, const auto& k) { + return stan::math::binomial_coefficient_log(n, k); + }; + + expect_ad(f, 5, 3); + expect_ad(f, 1, 0); + expect_ad(f, 0, 1); + expect_ad(f, -0.3, 0.5); +} + namespace binomial_coefficient_log_test_internal { struct TestValue { double n; @@ -265,6 +299,8 @@ TEST(MathFunctions, binomial_coefficient_log_precomputed) { using stan::math::is_nan; using stan::math::value_of; using stan::math::var; + using stan::test::expect_near_rel; + using stan::test::relative_tolerance; for (TestValue t : testValues) { std::stringstream msg; @@ -285,20 +321,19 @@ TEST(MathFunctions, binomial_coefficient_log_precomputed) { EXPECT_FALSE(is_nan(gradients[i])); } - double tol_val = std::max(1e-14 * fabs(t.val), 1e-14); - EXPECT_NEAR(value_of(val), t.val, tol_val) << msg.str(); + expect_near_rel(msg.str(), value_of(val), t.val, relative_tolerance(1e-14, 1e-14)); - std::function tol_grad; + relative_tolerance tol_grad; if (n < 1 || k < 1) { - tol_grad = [](double x) { return std::max(fabs(x) * 1e-8, 1e-7); }; + tol_grad = relative_tolerance(1e-8, 1e-7); } else { - tol_grad = [](double x) { return std::max(fabs(x) * 1e-10, 1e-8); }; + tol_grad = relative_tolerance(1e-10, 1e-8); } if (!is_nan(t.dn)) { - EXPECT_NEAR(gradients[0], t.dn, tol_grad(t.dn)) << "dn: " << msg.str(); + expect_near_rel(std::string("dn: ") + msg.str(), gradients[0], t.dn, tol_grad); } if (!is_nan(t.dk)) { - EXPECT_NEAR(gradients[1], t.dk, tol_grad(t.dk)) << "dk: " << msg.str(); + expect_near_rel(std::string("dk: ") + msg.str(), gradients[1], t.dk, tol_grad); } } } diff --git a/test/unit/math/rev/fun/lbeta_test.cpp b/test/unit/math/rev/fun/lbeta_test.cpp index e02d77a038a..3772fbe0f44 100644 --- a/test/unit/math/rev/fun/lbeta_test.cpp +++ b/test/unit/math/rev/fun/lbeta_test.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include #include #include #include @@ -8,81 +10,74 @@ #include #include -struct relative_tolerance { - double tol; - double tol_min; +namespace lbeta_test_internal { + // TODO(martinmodrak) the function here should be replaced by helpers for testing identities once those are available - double operator()(const double x) const { - return std::max(tol * fabs(x), tol_min); - } + struct identity_tolerances { + stan::test::relative_tolerance value; + stan::test::relative_tolerance gradient; + }; - double operator()(const double x, const double y) const { - return std::max(tol * 0.5 * (fabs(x) + fabs(y)), tol_min); - } -}; + template + void expect_identity(const std::string& msg, + const identity_tolerances& tolerances, const F1 lh, + const F2 rh, double x_dbl, double y_dbl) { + using stan::math::var; + using stan::test::expect_near_rel; -struct identity_tolerances { - relative_tolerance value; - relative_tolerance gradient; -}; + stan::math::nested_rev_autodiff nested; -template -void expect_identity(const std::string& msg, - const identity_tolerances& tolerances, const F1 lh, - const F2 rh, double x_dbl, double y_dbl) { - using stan::math::var; - - stan::math::start_nested(); + var x(x_dbl); + var y(y_dbl); - var x(x_dbl); - var y(y_dbl); + std::vector vars = {x, y}; - std::vector vars = {x, y}; + var left = lh(x, y); + double left_dbl = value_of(left); + std::vector gradients_left; + left.grad(vars, gradients_left); - var left = lh(x, y); - double left_dbl = value_of(left); - std::vector gradients_left; - left.grad(vars, gradients_left); + nested.set_zero_all_adjoints(); - stan::math::set_zero_all_adjoints_nested(); + var right = rh(x, y); + double right_dbl = value_of(right); + std::vector gradients_right; + right.grad(vars, gradients_right); - var right = rh(x, y); - double right_dbl = value_of(right); - std::vector gradients_right; - right.grad(vars, gradients_right); + std::stringstream args; + args << std::setprecision(22) << "args = [" << x << "," << y << "]"; + expect_near_rel(std::string() + args.str() + std::string(": ") + msg, left_dbl, right_dbl, tolerances.value); - std::stringstream args; - args << std::setprecision(22) << "args = [" << x << "," << y << "]"; - double tol_val = tolerances.value(left_dbl, right_dbl); - EXPECT_NEAR(left_dbl, right_dbl, tol_val) - << "value, " << args.str() << ": " << msg; + for (size_t i = 0; i < gradients_left.size(); ++i) { + std::stringstream grad_msg; + grad_msg << "grad_" << i << ", " << args.str() << ": " << msg; + expect_near_rel(grad_msg.str(), gradients_left[i], gradients_right[i], tolerances.gradient); + } - for (size_t i = 0; i < gradients_left.size(); ++i) { - double tol_grad - = tolerances.gradient(gradients_left[i], gradients_right[i]); - EXPECT_NEAR(gradients_left[i], gradients_right[i], tol_grad) - << "grad_" << i << ", " << args.str() << ": " << msg; } - - stan::math::recover_memory_nested(); } TEST(MathFunctions, lbeta_identities_gradient) { using stan::math::lbeta; using stan::math::pi; using stan::math::var; + using stan::test::expect_near_rel; std::vector to_test - = {1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 1e3, 1e30, 1e100}; + = {1e-100, 1e-8, 1e-1, 1, 2, 1 + 1e-6, 1e3, 1e30, 1e100}; - identity_tolerances tol{{1e-15, 1e-15}, {1e-10, 1e-10}}; + lbeta_test_internal::identity_tolerances tol{{1e-15, 1e-15}, {1e-10, 1e-10}}; for (double x : to_test) { for (double y : to_test) { + // TODO(martinmodrak) this restriction on testing should be lifted once the log_sum_exp bug (#1679) is resolved + if(x > 1e10 || y > 1e10) { + continue; + } auto rh = [](const var& a, const var& b) { return stan::math::log_sum_exp(lbeta(a + 1, b), lbeta(a, b + 1)); }; - expect_identity("succesors", tol, + lbeta_test_internal::expect_identity("succesors", tol, static_cast(lbeta), rh, x, y); } @@ -94,7 +89,7 @@ TEST(MathFunctions, lbeta_identities_gradient) { msg << std::setprecision(22) << "sin: x = " << x; double lh = lbeta(x, 1.0 - x); double rh = log(pi()) - log(sin(pi() * x)); - EXPECT_NEAR(lh, rh, tol.value(lh, rh)) << msg.str(); + expect_near_rel(msg.str(), lh, rh, tol.value); } } @@ -103,7 +98,7 @@ TEST(MathFunctions, lbeta_identities_gradient) { msg << std::setprecision(22) << "inv: x = " << x; double lh = lbeta(x, 1.0); double rh = -log(x); - EXPECT_NEAR(lh, rh, tol.value(lh, rh)) << msg.str(); + expect_near_rel(msg.str(), lh, rh, tol.value); } } From 215cad07e23b01e8806d2dbae78db4db9106a6c7 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 7 Mar 2020 14:58:28 -0500 Subject: [PATCH 18/26] [Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/stable/2017-11-14) --- .../prim/fun/binomial_coefficient_log.hpp | 2 +- .../math/rev/fun/binomial_coefficient_log.hpp | 12 +-- .../rev/fun/binomial_coefficient_log_test.cpp | 38 +++++---- test/unit/math/rev/fun/lbeta_test.cpp | 79 ++++++++++--------- 4 files changed, 70 insertions(+), 61 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 938ee6dc42a..80a86bdac82 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -66,7 +66,7 @@ namespace math { * * @tparam T_N type of the first argument * @tparam T_n type of the second argument - * + * * @param N total number of objects. * @param n number of objects chosen. * @return log (N choose n). diff --git a/stan/math/rev/fun/binomial_coefficient_log.hpp b/stan/math/rev/fun/binomial_coefficient_log.hpp index a7cadfbe58c..a79735704a3 100644 --- a/stan/math/rev/fun/binomial_coefficient_log.hpp +++ b/stan/math/rev/fun/binomial_coefficient_log.hpp @@ -27,7 +27,8 @@ class binomial_coefficient_log_vd_vari : public op_vd_vari { binomial_coefficient_log_vd_vari(vari* avi, double b) : op_vd_vari(binomial_coefficient_log(avi->val_, b), avi, b) {} void chain() { - avi_->adj_ += adj_ * (digamma(avi_->val_ + 1) - digamma(avi_->val_ - bd_ + 1)); + avi_->adj_ + += adj_ * (digamma(avi_->val_ + 1) - digamma(avi_->val_ - bd_ + 1)); } }; @@ -36,7 +37,8 @@ class binomial_coefficient_log_dv_vari : public op_dv_vari { binomial_coefficient_log_dv_vari(double a, vari* bvi) : op_dv_vari(binomial_coefficient_log(a, bvi->val_), a, bvi) {} void chain() { - bvi_->adj_ += adj_ * (digamma(ad_ - bvi_->val_ + 1) - digamma(bvi_->val_ + 1)); + bvi_->adj_ + += adj_ * (digamma(ad_ - bvi_->val_ + 1) - digamma(bvi_->val_ + 1)); } }; } // namespace internal @@ -44,7 +46,7 @@ class binomial_coefficient_log_dv_vari : public op_dv_vari { /** * Return the log of the binomial coefficient for the specified * arguments and its gradients. - * + * * See the docs for the prim version for all relevant formulae. * @param a var Argument * @param b var Argument @@ -57,7 +59,7 @@ inline var binomial_coefficient_log(const var& a, const var& b) { /** * Return the log of the binomial coefficient for the specified * arguments and its gradients. - * + * * See the docs for the prim version for all relevant formulae. * @param a var Argument * @param b double Argument @@ -70,7 +72,7 @@ inline var binomial_coefficient_log(const var& a, double b) { /** * Return the log of the binomial coefficient for the specified * arguments and its gradients. - * + * * See the docs for the prim version for all relevant formulae. * @param a double Argument * @param b var Argument diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 7a7179ce4cd..499cf3a6e08 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -22,8 +22,7 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { = {-0.1, 0, 1e-100, 1e-8, 1e-1, 1, 1 + 1e-6, 15, 10, 1e3, 1e30, 1e100}; std::vector k_ratios_to_test - = { -0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10 - }; + = {-0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10}; // Recurrence relation for (double n_dbl : n_to_test) { @@ -42,12 +41,13 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { var n(n_dbl); var k(k_dbl); - // TODO(martinmodrak) Use the framework for testing identities, once it is ready + // TODO(martinmodrak) Use the framework for testing identities, once it is + // ready var val_left = binomial_coefficient_log(n, k); var val_right_partial; var val_right; // Choose the more stable identity - if(n_dbl > 1 && k_dbl > 1 && (n_dbl - 1) + 1 - k_dbl > 0 ) { + if (n_dbl > 1 && k_dbl > 1 && (n_dbl - 1) + 1 - k_dbl > 0) { val_right_partial = binomial_coefficient_log(n - 1, k - 1); val_right = val_right_partial + log(n) - log(k); } else { @@ -73,17 +73,18 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { } std::stringstream msg; - msg << std::setprecision(22) << " successor: n = " << n - << ", k = " << k << std::endl - << "val = " << val_left - << ", val2 = " << val_right_partial << std::endl - << ", logn = " << log(n) - << ", logk = " << log(k); + msg << std::setprecision(22) << " successor: n = " << n << ", k = " << k + << std::endl + << "val = " << val_left << ", val2 = " << val_right_partial + << std::endl + << ", logn = " << log(n) << ", logk = " << log(k); - - expect_near_rel(std::string("val") + msg.str(), value_of(val_left), value_of(val_right)); - expect_near_rel(std::string("dn") + msg.str(), gradients_left[0], gradients_right[0]); - expect_near_rel(std::string("dk") + msg.str(), gradients_left[1], gradients_right[1]); + expect_near_rel(std::string("val") + msg.str(), value_of(val_left), + value_of(val_right)); + expect_near_rel(std::string("dn") + msg.str(), gradients_left[0], + gradients_right[0]); + expect_near_rel(std::string("dk") + msg.str(), gradients_left[1], + gradients_right[1]); } } } @@ -321,7 +322,8 @@ TEST(MathFunctions, binomial_coefficient_log_precomputed) { EXPECT_FALSE(is_nan(gradients[i])); } - expect_near_rel(msg.str(), value_of(val), t.val, relative_tolerance(1e-14, 1e-14)); + expect_near_rel(msg.str(), value_of(val), t.val, + relative_tolerance(1e-14, 1e-14)); relative_tolerance tol_grad; if (n < 1 || k < 1) { @@ -330,10 +332,12 @@ TEST(MathFunctions, binomial_coefficient_log_precomputed) { tol_grad = relative_tolerance(1e-10, 1e-8); } if (!is_nan(t.dn)) { - expect_near_rel(std::string("dn: ") + msg.str(), gradients[0], t.dn, tol_grad); + expect_near_rel(std::string("dn: ") + msg.str(), gradients[0], t.dn, + tol_grad); } if (!is_nan(t.dk)) { - expect_near_rel(std::string("dk: ") + msg.str(), gradients[1], t.dk, tol_grad); + expect_near_rel(std::string("dk: ") + msg.str(), gradients[1], t.dk, + tol_grad); } } } diff --git a/test/unit/math/rev/fun/lbeta_test.cpp b/test/unit/math/rev/fun/lbeta_test.cpp index 3772fbe0f44..fb24c8879e7 100644 --- a/test/unit/math/rev/fun/lbeta_test.cpp +++ b/test/unit/math/rev/fun/lbeta_test.cpp @@ -11,51 +11,53 @@ #include namespace lbeta_test_internal { - // TODO(martinmodrak) the function here should be replaced by helpers for testing identities once those are available +// TODO(martinmodrak) the function here should be replaced by helpers for +// testing identities once those are available - struct identity_tolerances { - stan::test::relative_tolerance value; - stan::test::relative_tolerance gradient; - }; - - template - void expect_identity(const std::string& msg, - const identity_tolerances& tolerances, const F1 lh, - const F2 rh, double x_dbl, double y_dbl) { - using stan::math::var; - using stan::test::expect_near_rel; +struct identity_tolerances { + stan::test::relative_tolerance value; + stan::test::relative_tolerance gradient; +}; - stan::math::nested_rev_autodiff nested; +template +void expect_identity(const std::string& msg, + const identity_tolerances& tolerances, const F1 lh, + const F2 rh, double x_dbl, double y_dbl) { + using stan::math::var; + using stan::test::expect_near_rel; - var x(x_dbl); - var y(y_dbl); + stan::math::nested_rev_autodiff nested; - std::vector vars = {x, y}; + var x(x_dbl); + var y(y_dbl); - var left = lh(x, y); - double left_dbl = value_of(left); - std::vector gradients_left; - left.grad(vars, gradients_left); + std::vector vars = {x, y}; - nested.set_zero_all_adjoints(); + var left = lh(x, y); + double left_dbl = value_of(left); + std::vector gradients_left; + left.grad(vars, gradients_left); - var right = rh(x, y); - double right_dbl = value_of(right); - std::vector gradients_right; - right.grad(vars, gradients_right); + nested.set_zero_all_adjoints(); - std::stringstream args; - args << std::setprecision(22) << "args = [" << x << "," << y << "]"; - expect_near_rel(std::string() + args.str() + std::string(": ") + msg, left_dbl, right_dbl, tolerances.value); + var right = rh(x, y); + double right_dbl = value_of(right); + std::vector gradients_right; + right.grad(vars, gradients_right); - for (size_t i = 0; i < gradients_left.size(); ++i) { - std::stringstream grad_msg; - grad_msg << "grad_" << i << ", " << args.str() << ": " << msg; - expect_near_rel(grad_msg.str(), gradients_left[i], gradients_right[i], tolerances.gradient); - } + std::stringstream args; + args << std::setprecision(22) << "args = [" << x << "," << y << "]"; + expect_near_rel(std::string() + args.str() + std::string(": ") + msg, + left_dbl, right_dbl, tolerances.value); + for (size_t i = 0; i < gradients_left.size(); ++i) { + std::stringstream grad_msg; + grad_msg << "grad_" << i << ", " << args.str() << ": " << msg; + expect_near_rel(grad_msg.str(), gradients_left[i], gradients_right[i], + tolerances.gradient); } } +} // namespace lbeta_test_internal TEST(MathFunctions, lbeta_identities_gradient) { using stan::math::lbeta; @@ -70,16 +72,17 @@ TEST(MathFunctions, lbeta_identities_gradient) { for (double x : to_test) { for (double y : to_test) { - // TODO(martinmodrak) this restriction on testing should be lifted once the log_sum_exp bug (#1679) is resolved - if(x > 1e10 || y > 1e10) { + // TODO(martinmodrak) this restriction on testing should be lifted once + // the log_sum_exp bug (#1679) is resolved + if (x > 1e10 || y > 1e10) { continue; } auto rh = [](const var& a, const var& b) { return stan::math::log_sum_exp(lbeta(a + 1, b), lbeta(a, b + 1)); }; - lbeta_test_internal::expect_identity("succesors", tol, - static_cast(lbeta), rh, - x, y); + lbeta_test_internal::expect_identity( + "succesors", tol, static_cast(lbeta), + rh, x, y); } } From 71cbd40817499f76e4428007c4623e1956d016b0 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Sat, 7 Mar 2020 23:06:50 +0100 Subject: [PATCH 19/26] Resolved int and fwd issues --- stan/math/prim/fun/binomial_coefficient_log.hpp | 7 ++----- .../math/mix/fun/binomial_coefficient_log_test.cpp | 4 ++++ .../math/rev/fun/binomial_coefficient_log_test.cpp | 14 -------------- 3 files changed, 6 insertions(+), 19 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 938ee6dc42a..81a531726ca 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -87,13 +87,10 @@ inline return_type_t binomial_coefficient_log(const T_N N, check_greater_or_equal(function, "(first argument - second argument + 1)", N_plus_1 - n, 0.0); - if (n == 0) { - return 0; - } - if (N / 2 < n) { + if (N > 0 && N / 2 + 1 < n) { return binomial_coefficient_log(N, N - n); } - if (N_plus_1 < lgamma_stirling_diff_useful) { + if (N_plus_1 < lgamma_stirling_diff_useful || n == 0) { return lgamma(N_plus_1) - lgamma(n + 1) - lgamma(N_plus_1 - n); } return -lbeta(N - n + 1, n + 1) - log1p(N); diff --git a/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp b/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp index 84defe83a3e..a47c850ece6 100644 --- a/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp @@ -6,5 +6,9 @@ TEST(mathMixScalFun, binomialCoefficientLog) { }; stan::test::expect_ad(f, 3, 2); stan::test::expect_ad(f, 24.0, 12.0); + stan::test::expect_ad(f, 1.0, 0.0); + stan::test::expect_ad(f, 0.0, 1.0); + stan::test::expect_ad(f, -0.3, 0.5); + stan::test::expect_common_nonzero_binary(f); } diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 7a7179ce4cd..75664e6de4a 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -1,7 +1,6 @@ #include #include #include -#include #include #include #include @@ -88,19 +87,6 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { } } -TEST(MathFunctions, binomial_coefficient_log_ad) { - using stan::test::expect_ad; - - auto f = [](const auto& n, const auto& k) { - return stan::math::binomial_coefficient_log(n, k); - }; - - expect_ad(f, 5, 3); - expect_ad(f, 1, 0); - expect_ad(f, 0, 1); - expect_ad(f, -0.3, 0.5); -} - namespace binomial_coefficient_log_test_internal { struct TestValue { double n; From 4fa9c52fcb3497f8232e31f614146d95d9965ae9 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Mon, 9 Mar 2020 09:56:01 +0100 Subject: [PATCH 20/26] Attempt to handle edge cases for derivatives --- .../prim/fun/binomial_coefficient_log.hpp | 5 ++++- .../math/rev/fun/binomial_coefficient_log.hpp | 19 ++++++++++++++++--- .../mix/fun/binomial_coefficient_log_test.cpp | 3 +++ .../rev/fun/binomial_coefficient_log_test.cpp | 13 ++++++++++++- 4 files changed, 35 insertions(+), 5 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index ce88cebdfcb..9cd7626bae0 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -87,10 +87,13 @@ inline return_type_t binomial_coefficient_log(const T_N N, check_greater_or_equal(function, "(first argument - second argument + 1)", N_plus_1 - n, 0.0); + if (n == 0) { + return 0; + } if (N > 0 && N / 2 + 1 < n) { return binomial_coefficient_log(N, N - n); } - if (N_plus_1 < lgamma_stirling_diff_useful || n == 0) { + if (N_plus_1 < lgamma_stirling_diff_useful) { return lgamma(N_plus_1) - lgamma(n + 1) - lgamma(N_plus_1 - n); } return -lbeta(N - n + 1, n + 1) - log1p(N); diff --git a/stan/math/rev/fun/binomial_coefficient_log.hpp b/stan/math/rev/fun/binomial_coefficient_log.hpp index a79735704a3..94f055089c3 100644 --- a/stan/math/rev/fun/binomial_coefficient_log.hpp +++ b/stan/math/rev/fun/binomial_coefficient_log.hpp @@ -15,10 +15,23 @@ class binomial_coefficient_log_vv_vari : public op_vv_vari { binomial_coefficient_log_vv_vari(vari* avi, vari* bvi) : op_vv_vari(binomial_coefficient_log(avi->val_, bvi->val_), avi, bvi) {} void chain() { - double digamma_ambp1 = digamma(avi_->val_ - bvi_->val_ + 1); + if(bvi_->val_ == 0) { + // Ignoring avi, the derivative is zero + if(avi_->val_ == -1) { + bvi_->adj_ += adj_ * stan::math::NEGATIVE_INFTY; + } else { + bvi_->adj_ += adj_ * (boost::math::constants::euler() + digamma(avi_->val_ + 1)); + } + } else { + double digamma_ambp1 = digamma(avi_->val_ - bvi_->val_ + 1); - avi_->adj_ += adj_ * (digamma(avi_->val_ + 1) - digamma_ambp1); - bvi_->adj_ += adj_ * (digamma_ambp1 - digamma(bvi_->val_ + 1)); + avi_->adj_ += adj_ * (digamma(avi_->val_ + 1) - digamma_ambp1); + if(bvi_->val_ == -1) { + bvi_->adj_ += adj_ * stan::math::INFTY; + } else { + bvi_->adj_ += adj_ * (digamma_ambp1 - digamma(bvi_->val_ + 1)); + } + } } }; diff --git a/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp b/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp index a47c850ece6..5c6665c80da 100644 --- a/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp @@ -10,5 +10,8 @@ TEST(mathMixScalFun, binomialCoefficientLog) { stan::test::expect_ad(f, 0.0, 1.0); stan::test::expect_ad(f, -0.3, 0.5); + stan::test::expect_ad(f, -1.0, 0.0); + stan::test::expect_ad(f, 1.3, -1.0); + stan::test::expect_common_nonzero_binary(f); } diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 2e2474c2765..2fdc350004f 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -97,6 +97,14 @@ struct TestValue { double dk; }; +// Hand-checked edge cases +std::vector testValuesEdge = { + {-1, 0, 0, 0, stan::math::NEGATIVE_INFTY}, + {-1, 0, 0, 0, stan::math::NEGATIVE_INFTY}, + {1, 0 , 0 , 0, 1}, + {3, - 1, stan::math::NEGATIVE_INFTY, -0.25, stan::math::INFTY} +}; + const double NaN = stan::math::NOT_A_NUMBER; // Test values generated in Mathematica, reproducible notebook at // https://www.wolframcloud.com/obj/martin.modrak/Published/binomial_coefficient_log.nb @@ -283,13 +291,16 @@ std::vector testValues = { TEST(MathFunctions, binomial_coefficient_log_precomputed) { using binomial_coefficient_log_test_internal::TestValue; using binomial_coefficient_log_test_internal::testValues; + using binomial_coefficient_log_test_internal::testValuesEdge; using stan::math::is_nan; using stan::math::value_of; using stan::math::var; using stan::test::expect_near_rel; using stan::test::relative_tolerance; - for (TestValue t : testValues) { + std::vector allTestValues = testValues; + allTestValues.insert(allTestValues.end(), testValuesEdge.begin(), testValuesEdge.end()); + for (TestValue t : allTestValues) { std::stringstream msg; msg << std::setprecision(22) << "n = " << t.n << ", k = " << t.k; From eab19c350aa10f63c40b017124b00e655af3cbc6 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Mon, 9 Mar 2020 13:18:40 +0100 Subject: [PATCH 21/26] Moved to operands_and_partials, fixed edge cases --- .../prim/fun/binomial_coefficient_log.hpp | 97 ++++++++++++----- stan/math/rev/fun.hpp | 1 - .../math/rev/fun/binomial_coefficient_log.hpp | 100 ------------------ .../mix/fun/binomial_coefficient_log_test.cpp | 6 +- .../rev/fun/binomial_coefficient_log_test.cpp | 9 +- 5 files changed, 78 insertions(+), 135 deletions(-) delete mode 100644 stan/math/rev/fun/binomial_coefficient_log.hpp diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 9cd7626bae0..7f652a322c9 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -1,13 +1,16 @@ #ifndef STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP #define STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP +#include #include #include #include +#include #include #include #include #include +#include namespace stan { namespace math { @@ -16,16 +19,16 @@ namespace math { * Return the log of the binomial coefficient for the specified * arguments. * - * The binomial coefficient, \f${N \choose n}\f$, read "N choose n", is - * defined for \f$0 \leq n \leq N\f$ by + * The binomial coefficient, \f${n \choose k}\f$, read "n choose k", is + * defined for \f$0 \leq k \leq n\f$ by * - * \f${N \choose n} = \frac{N!}{n! (N-n)!}\f$. + * \f${n \choose k} = \frac{n!}{k! (n-k)!}\f$. * * This function uses Gamma functions to define the log - * and generalize the arguments to continuous N and n. + * and generalize the arguments to continuous n and k. * - * \f$ \log {N \choose n} - * = \log \ \Gamma(N+1) - \log \Gamma(n+1) - \log \Gamma(N-n+1)\f$. + * \f$ \log {n \choose k} + * = \log \ \Gamma(n+1) - \log \Gamma(k+1) - \log \Gamma(n-k+1)\f$. * * \f[ @@ -64,39 +67,79 @@ namespace math { * * This function is numerically more stable than naive evaluation via lgamma. * - * @tparam T_N type of the first argument - * @tparam T_n type of the second argument + * @tparam T_n type of the first argument + * @tparam T_k type of the second argument * - * @param N total number of objects. - * @param n number of objects chosen. - * @return log (N choose n). + * @param n total number of objects. + * @param k number of objects chosen. + * @return log (n choose k). */ -template -inline return_type_t binomial_coefficient_log(const T_N N, - const T_n n) { - if (is_any_nan(N, n)) { +template +inline return_type_t binomial_coefficient_log(const T_n n, + const T_k k) { + if (is_any_nan(n, k)) { return stan::math::NOT_A_NUMBER; } - const T_N N_plus_1 = N + 1; + // Choosing the more stable of the symmetric branches + if( n > 0 && k > value_of_rec(n) / 2.0 + 1e-8) { + return binomial_coefficient_log(n, n - k); + } + + using T_partials_return = partials_return_t; + + const T_partials_return n_ = value_of(n); + const T_partials_return k_ = value_of(k); + const T_partials_return n_plus_1 = n_ + 1; + const T_partials_return n_plus_1_mk = n_plus_1 - k_; static const char* function = "binomial_coefficient_log"; - check_greater_or_equal(function, "first argument", N, -1); - check_greater_or_equal(function, "second argument", n, -1); + check_greater_or_equal(function, "first argument", n, -1); + check_greater_or_equal(function, "second argument", k, -1); check_greater_or_equal(function, "(first argument - second argument + 1)", - N_plus_1 - n, 0.0); + n_plus_1_mk, 0.0); - if (n == 0) { - return 0; - } - if (N > 0 && N / 2 + 1 < n) { - return binomial_coefficient_log(N, N - n); + operands_and_partials ops_partials(n, k); + + T_partials_return value; + if (k_ == 0) { + value = 0; + } else if (n_plus_1 < lgamma_stirling_diff_useful) { + value = lgamma(n_plus_1) - lgamma(k_ + 1) - lgamma(n_plus_1_mk); + } else { + value = -lbeta(n_plus_1_mk, k_ + 1) - log1p(n_); } - if (N_plus_1 < lgamma_stirling_diff_useful) { - return lgamma(N_plus_1) - lgamma(n + 1) - lgamma(N_plus_1 - n); + + if (!is_constant_all::value) { + // Branching on all the edge cases. + // In direct computation many of those would be NaN + // But one-sided limits from within the domain exist. + T_partials_return digamma_n_plus_1_mk = digamma(n_plus_1_mk); + + if (!is_constant_all::value) { + if(n_ == -1.0) { + if(k_ == 0) { + ops_partials.edge1_.partials_[0] = 0; + } else { + ops_partials.edge1_.partials_[0] = stan::math::NEGATIVE_INFTY;; + } + } else { + ops_partials.edge1_.partials_[0] = (digamma(n_plus_1) - digamma_n_plus_1_mk); + } + } + if (!is_constant_all::value) { + if(k_ == 0 && n_ == -1.0) { + ops_partials.edge2_.partials_[0] = stan::math::NEGATIVE_INFTY; + } else if(k_ == -1) { + ops_partials.edge2_.partials_[0] = stan::math::INFTY; + } else { + ops_partials.edge2_.partials_[0] = (digamma_n_plus_1_mk - digamma(k_ + 1)); + } + } } - return -lbeta(N - n + 1, n + 1) - log1p(N); + + return ops_partials.build(value); } } // namespace math diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index d692cdf0513..1deb92f9738 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -22,7 +22,6 @@ #include #include #include -#include #include #include #include diff --git a/stan/math/rev/fun/binomial_coefficient_log.hpp b/stan/math/rev/fun/binomial_coefficient_log.hpp deleted file mode 100644 index 94f055089c3..00000000000 --- a/stan/math/rev/fun/binomial_coefficient_log.hpp +++ /dev/null @@ -1,100 +0,0 @@ -#ifndef STAN_MATH_REV_FUN_BINOMIAL_COEFFICIENT_LOG_HPP -#define STAN_MATH_REV_FUN_BINOMIAL_COEFFICIENT_LOG_HPP - -#include -#include -#include -#include - -namespace stan { -namespace math { - -namespace internal { -class binomial_coefficient_log_vv_vari : public op_vv_vari { - public: - binomial_coefficient_log_vv_vari(vari* avi, vari* bvi) - : op_vv_vari(binomial_coefficient_log(avi->val_, bvi->val_), avi, bvi) {} - void chain() { - if(bvi_->val_ == 0) { - // Ignoring avi, the derivative is zero - if(avi_->val_ == -1) { - bvi_->adj_ += adj_ * stan::math::NEGATIVE_INFTY; - } else { - bvi_->adj_ += adj_ * (boost::math::constants::euler() + digamma(avi_->val_ + 1)); - } - } else { - double digamma_ambp1 = digamma(avi_->val_ - bvi_->val_ + 1); - - avi_->adj_ += adj_ * (digamma(avi_->val_ + 1) - digamma_ambp1); - if(bvi_->val_ == -1) { - bvi_->adj_ += adj_ * stan::math::INFTY; - } else { - bvi_->adj_ += adj_ * (digamma_ambp1 - digamma(bvi_->val_ + 1)); - } - } - } -}; - -class binomial_coefficient_log_vd_vari : public op_vd_vari { - public: - binomial_coefficient_log_vd_vari(vari* avi, double b) - : op_vd_vari(binomial_coefficient_log(avi->val_, b), avi, b) {} - void chain() { - avi_->adj_ - += adj_ * (digamma(avi_->val_ + 1) - digamma(avi_->val_ - bd_ + 1)); - } -}; - -class binomial_coefficient_log_dv_vari : public op_dv_vari { - public: - binomial_coefficient_log_dv_vari(double a, vari* bvi) - : op_dv_vari(binomial_coefficient_log(a, bvi->val_), a, bvi) {} - void chain() { - bvi_->adj_ - += adj_ * (digamma(ad_ - bvi_->val_ + 1) - digamma(bvi_->val_ + 1)); - } -}; -} // namespace internal - -/** - * Return the log of the binomial coefficient for the specified - * arguments and its gradients. - * - * See the docs for the prim version for all relevant formulae. - * @param a var Argument - * @param b var Argument - * @return Result of log (a choose b) - */ -inline var binomial_coefficient_log(const var& a, const var& b) { - return var(new internal::binomial_coefficient_log_vv_vari(a.vi_, b.vi_)); -} - -/** - * Return the log of the binomial coefficient for the specified - * arguments and its gradients. - * - * See the docs for the prim version for all relevant formulae. - * @param a var Argument - * @param b double Argument - * @return Result of log (a choose b) - */ -inline var binomial_coefficient_log(const var& a, double b) { - return var(new internal::binomial_coefficient_log_vd_vari(a.vi_, b)); -} - -/** - * Return the log of the binomial coefficient for the specified - * arguments and its gradients. - * - * See the docs for the prim version for all relevant formulae. - * @param a double Argument - * @param b var Argument - * @return Result of log (a choose b) - */ -inline var binomial_coefficient_log(double a, const var& b) { - return var(new internal::binomial_coefficient_log_dv_vari(a, b.vi_)); -} - -} // namespace math -} // namespace stan -#endif diff --git a/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp b/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp index 5c6665c80da..09f6e96fef9 100644 --- a/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp @@ -1,5 +1,6 @@ #include + TEST(mathMixScalFun, binomialCoefficientLog) { auto f = [](const auto& x1, const auto& x2) { return stan::math::binomial_coefficient_log(x1, x2); @@ -9,9 +10,6 @@ TEST(mathMixScalFun, binomialCoefficientLog) { stan::test::expect_ad(f, 1.0, 0.0); stan::test::expect_ad(f, 0.0, 1.0); stan::test::expect_ad(f, -0.3, 0.5); - - stan::test::expect_ad(f, -1.0, 0.0); - stan::test::expect_ad(f, 1.3, -1.0); - + stan::test::expect_common_nonzero_binary(f); } diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 2fdc350004f..6cbcd2300a4 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -97,12 +97,15 @@ struct TestValue { double dk; }; -// Hand-checked edge cases +// Hand-checked edge cases. Using one-sided limits from +// within the function domain where the value doesn't exist std::vector testValuesEdge = { {-1, 0, 0, 0, stan::math::NEGATIVE_INFTY}, - {-1, 0, 0, 0, stan::math::NEGATIVE_INFTY}, + {0, -1, stan::math::NEGATIVE_INFTY, -1.0, stan::math::INFTY}, + {3, - 1, stan::math::NEGATIVE_INFTY, -0.25, stan::math::INFTY}, + {-1, -0.2, stan::math::INFTY, stan::math::NEGATIVE_INFTY, -4.324031329886049836 }, + {4.0, 5.0, stan::math::NEGATIVE_INFTY, stan::math::INFTY, stan::math::NEGATIVE_INFTY}, {1, 0 , 0 , 0, 1}, - {3, - 1, stan::math::NEGATIVE_INFTY, -0.25, stan::math::INFTY} }; const double NaN = stan::math::NOT_A_NUMBER; From 75a6f735104fe8aebede3cbdeacfaeff4f51adc5 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 9 Mar 2020 08:19:50 -0400 Subject: [PATCH 22/26] [Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/stable/2017-11-14) --- .../prim/fun/binomial_coefficient_log.hpp | 21 +++++++++++-------- .../mix/fun/binomial_coefficient_log_test.cpp | 3 +-- .../rev/fun/binomial_coefficient_log_test.cpp | 17 ++++++++------- 3 files changed, 23 insertions(+), 18 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 7f652a322c9..62807e695cc 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -83,9 +83,9 @@ inline return_type_t binomial_coefficient_log(const T_n n, } // Choosing the more stable of the symmetric branches - if( n > 0 && k > value_of_rec(n) / 2.0 + 1e-8) { + if (n > 0 && k > value_of_rec(n) / 2.0 + 1e-8) { return binomial_coefficient_log(n, n - k); - } + } using T_partials_return = partials_return_t; @@ -118,23 +118,26 @@ inline return_type_t binomial_coefficient_log(const T_n n, T_partials_return digamma_n_plus_1_mk = digamma(n_plus_1_mk); if (!is_constant_all::value) { - if(n_ == -1.0) { - if(k_ == 0) { + if (n_ == -1.0) { + if (k_ == 0) { ops_partials.edge1_.partials_[0] = 0; } else { - ops_partials.edge1_.partials_[0] = stan::math::NEGATIVE_INFTY;; + ops_partials.edge1_.partials_[0] = stan::math::NEGATIVE_INFTY; + ; } } else { - ops_partials.edge1_.partials_[0] = (digamma(n_plus_1) - digamma_n_plus_1_mk); + ops_partials.edge1_.partials_[0] + = (digamma(n_plus_1) - digamma_n_plus_1_mk); } } if (!is_constant_all::value) { - if(k_ == 0 && n_ == -1.0) { + if (k_ == 0 && n_ == -1.0) { ops_partials.edge2_.partials_[0] = stan::math::NEGATIVE_INFTY; - } else if(k_ == -1) { + } else if (k_ == -1) { ops_partials.edge2_.partials_[0] = stan::math::INFTY; } else { - ops_partials.edge2_.partials_[0] = (digamma_n_plus_1_mk - digamma(k_ + 1)); + ops_partials.edge2_.partials_[0] + = (digamma_n_plus_1_mk - digamma(k_ + 1)); } } } diff --git a/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp b/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp index 09f6e96fef9..a47c850ece6 100644 --- a/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/mix/fun/binomial_coefficient_log_test.cpp @@ -1,6 +1,5 @@ #include - TEST(mathMixScalFun, binomialCoefficientLog) { auto f = [](const auto& x1, const auto& x2) { return stan::math::binomial_coefficient_log(x1, x2); @@ -10,6 +9,6 @@ TEST(mathMixScalFun, binomialCoefficientLog) { stan::test::expect_ad(f, 1.0, 0.0); stan::test::expect_ad(f, 0.0, 1.0); stan::test::expect_ad(f, -0.3, 0.5); - + stan::test::expect_common_nonzero_binary(f); } diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 6cbcd2300a4..167cc4f6327 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -100,12 +100,14 @@ struct TestValue { // Hand-checked edge cases. Using one-sided limits from // within the function domain where the value doesn't exist std::vector testValuesEdge = { - {-1, 0, 0, 0, stan::math::NEGATIVE_INFTY}, - {0, -1, stan::math::NEGATIVE_INFTY, -1.0, stan::math::INFTY}, - {3, - 1, stan::math::NEGATIVE_INFTY, -0.25, stan::math::INFTY}, - {-1, -0.2, stan::math::INFTY, stan::math::NEGATIVE_INFTY, -4.324031329886049836 }, - {4.0, 5.0, stan::math::NEGATIVE_INFTY, stan::math::INFTY, stan::math::NEGATIVE_INFTY}, - {1, 0 , 0 , 0, 1}, + {-1, 0, 0, 0, stan::math::NEGATIVE_INFTY}, + {0, -1, stan::math::NEGATIVE_INFTY, -1.0, stan::math::INFTY}, + {3, -1, stan::math::NEGATIVE_INFTY, -0.25, stan::math::INFTY}, + {-1, -0.2, stan::math::INFTY, stan::math::NEGATIVE_INFTY, + -4.324031329886049836}, + {4.0, 5.0, stan::math::NEGATIVE_INFTY, stan::math::INFTY, + stan::math::NEGATIVE_INFTY}, + {1, 0, 0, 0, 1}, }; const double NaN = stan::math::NOT_A_NUMBER; @@ -302,7 +304,8 @@ TEST(MathFunctions, binomial_coefficient_log_precomputed) { using stan::test::relative_tolerance; std::vector allTestValues = testValues; - allTestValues.insert(allTestValues.end(), testValuesEdge.begin(), testValuesEdge.end()); + allTestValues.insert(allTestValues.end(), testValuesEdge.begin(), + testValuesEdge.end()); for (TestValue t : allTestValues) { std::stringstream msg; msg << std::setprecision(22) << "n = " << t.n << ", k = " << t.k; From 47c79624767fa26390c9549a2a39725c0532d975 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Mon, 9 Mar 2020 16:27:45 +0100 Subject: [PATCH 23/26] Fixed lint --- stan/math/prim/fun/binomial_coefficient_log.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 62807e695cc..db93d947018 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -123,7 +123,6 @@ inline return_type_t binomial_coefficient_log(const T_n n, ops_partials.edge1_.partials_[0] = 0; } else { ops_partials.edge1_.partials_[0] = stan::math::NEGATIVE_INFTY; - ; } } else { ops_partials.edge1_.partials_[0] From cf05d408a6bfb8bc72ccf8e7967e6086fb486fdc Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Thu, 12 Mar 2020 10:56:25 +0100 Subject: [PATCH 24/26] Handling review comments --- .../prim/fun/binomial_coefficient_log.hpp | 42 ++++++++++--------- .../rev/fun/binomial_coefficient_log_test.cpp | 3 +- test/unit/math/rev/fun/lbeta_test.cpp | 4 ++ 3 files changed, 29 insertions(+), 20 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index db93d947018..9e012bcb83a 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -1,7 +1,6 @@ #ifndef STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP #define STAN_MATH_PRIM_FUN_BINOMIAL_COEFFICIENT_LOG_HPP -#include #include #include #include @@ -78,8 +77,10 @@ namespace math { template inline return_type_t binomial_coefficient_log(const T_n n, const T_k k) { + using T_partials_return = partials_return_t; + if (is_any_nan(n, k)) { - return stan::math::NOT_A_NUMBER; + return NOT_A_NUMBER; } // Choosing the more stable of the symmetric branches @@ -87,12 +88,11 @@ inline return_type_t binomial_coefficient_log(const T_n n, return binomial_coefficient_log(n, n - k); } - using T_partials_return = partials_return_t; - const T_partials_return n_ = value_of(n); - const T_partials_return k_ = value_of(k); - const T_partials_return n_plus_1 = n_ + 1; - const T_partials_return n_plus_1_mk = n_plus_1 - k_; + const T_partials_return n_dbl = value_of(n); + const T_partials_return k_dbl = value_of(k); + const T_partials_return n_plus_1 = n_dbl + 1; + const T_partials_return n_plus_1_mk = n_plus_1 - k_dbl; static const char* function = "binomial_coefficient_log"; check_greater_or_equal(function, "first argument", n, -1); @@ -103,26 +103,30 @@ inline return_type_t binomial_coefficient_log(const T_n n, operands_and_partials ops_partials(n, k); T_partials_return value; - if (k_ == 0) { + if (k_dbl == 0) { value = 0; } else if (n_plus_1 < lgamma_stirling_diff_useful) { - value = lgamma(n_plus_1) - lgamma(k_ + 1) - lgamma(n_plus_1_mk); + value = lgamma(n_plus_1) - lgamma(k_dbl + 1) - lgamma(n_plus_1_mk); } else { - value = -lbeta(n_plus_1_mk, k_ + 1) - log1p(n_); + value = -lbeta(n_plus_1_mk, k_dbl + 1) - log1p(n_dbl); } if (!is_constant_all::value) { // Branching on all the edge cases. // In direct computation many of those would be NaN - // But one-sided limits from within the domain exist. + // But one-sided limits from within the domain exist, all of the below + // follows from lim x->0 from above digamma(x) == -Inf + // + // Note that we have k < n / 2 (see the first branch in this function) + // se we can ignore the n == k - 1 edge case. T_partials_return digamma_n_plus_1_mk = digamma(n_plus_1_mk); if (!is_constant_all::value) { - if (n_ == -1.0) { - if (k_ == 0) { + if (n_dbl == -1.0) { + if (k_dbl == 0) { ops_partials.edge1_.partials_[0] = 0; } else { - ops_partials.edge1_.partials_[0] = stan::math::NEGATIVE_INFTY; + ops_partials.edge1_.partials_[0] = NEGATIVE_INFTY; } } else { ops_partials.edge1_.partials_[0] @@ -130,13 +134,13 @@ inline return_type_t binomial_coefficient_log(const T_n n, } } if (!is_constant_all::value) { - if (k_ == 0 && n_ == -1.0) { - ops_partials.edge2_.partials_[0] = stan::math::NEGATIVE_INFTY; - } else if (k_ == -1) { - ops_partials.edge2_.partials_[0] = stan::math::INFTY; + if (k_dbl == 0 && n_dbl == -1.0) { + ops_partials.edge2_.partials_[0] = NEGATIVE_INFTY; + } else if (k_dbl == -1) { + ops_partials.edge2_.partials_[0] = INFTY; } else { ops_partials.edge2_.partials_[0] - = (digamma_n_plus_1_mk - digamma(k_ + 1)); + = (digamma_n_plus_1_mk - digamma(k_dbl + 1)); } } } diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 167cc4f6327..0a05e8832cd 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -23,7 +23,8 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { std::vector k_ratios_to_test = {-0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10}; - // Recurrence relation + // Recurrence relation: binomial_coefficient_log(n, k) == + // binomial_coefficient_log(n - 1, k - 1) + log(n) - log(k) for (double n_dbl : n_to_test) { for (double k_ratio : k_ratios_to_test) { double k_dbl = n_dbl * k_ratio; diff --git a/test/unit/math/rev/fun/lbeta_test.cpp b/test/unit/math/rev/fun/lbeta_test.cpp index fb24c8879e7..f4ea4ccf6bb 100644 --- a/test/unit/math/rev/fun/lbeta_test.cpp +++ b/test/unit/math/rev/fun/lbeta_test.cpp @@ -70,6 +70,8 @@ TEST(MathFunctions, lbeta_identities_gradient) { lbeta_test_internal::identity_tolerances tol{{1e-15, 1e-15}, {1e-10, 1e-10}}; + // All identities from https://en.wikipedia.org/wiki/Beta_function#Properties + // Successors: beta(a,b) = beta(a + 1, b) + beta(a, b + 1) for (double x : to_test) { for (double y : to_test) { // TODO(martinmodrak) this restriction on testing should be lifted once @@ -86,6 +88,7 @@ TEST(MathFunctions, lbeta_identities_gradient) { } } + // Sin: beta(x, 1 - x) == pi / sin(pi * x) for (double x : to_test) { if (x < 1) { std::stringstream msg; @@ -96,6 +99,7 @@ TEST(MathFunctions, lbeta_identities_gradient) { } } + // Inv: beta(1, x) == 1 / x for (double x : to_test) { std::stringstream msg; msg << std::setprecision(22) << "inv: x = " << x; From c84a259cf877ffcaa35eb30ae536f1ced674b050 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 12 Mar 2020 06:00:03 -0400 Subject: [PATCH 25/26] [Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/stable/2017-11-14) --- stan/math/prim/fun/binomial_coefficient_log.hpp | 1 - test/unit/math/rev/fun/binomial_coefficient_log_test.cpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index 9e012bcb83a..b9e3b87964e 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -88,7 +88,6 @@ inline return_type_t binomial_coefficient_log(const T_n n, return binomial_coefficient_log(n, n - k); } - const T_partials_return n_dbl = value_of(n); const T_partials_return k_dbl = value_of(k); const T_partials_return n_plus_1 = n_dbl + 1; diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 0a05e8832cd..907c49244a1 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -23,7 +23,7 @@ TEST(MathFunctions, binomial_coefficient_log_identities) { std::vector k_ratios_to_test = {-0.1, 1e-10, 1e-5, 1e-3, 1e-1, 0.5, 0.9, 1 - 1e-5, 1 - 1e-10}; - // Recurrence relation: binomial_coefficient_log(n, k) == + // Recurrence relation: binomial_coefficient_log(n, k) == // binomial_coefficient_log(n - 1, k - 1) + log(n) - log(k) for (double n_dbl : n_to_test) { for (double k_ratio : k_ratios_to_test) { From 6aed316f5f8e562e86250a1406a98c474492e9f2 Mon Sep 17 00:00:00 2001 From: martin_cerny Date: Thu, 12 Mar 2020 12:25:02 +0100 Subject: [PATCH 26/26] One more edge case-test+fix --- stan/math/prim/fun/binomial_coefficient_log.hpp | 2 +- test/unit/math/rev/fun/binomial_coefficient_log_test.cpp | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/stan/math/prim/fun/binomial_coefficient_log.hpp b/stan/math/prim/fun/binomial_coefficient_log.hpp index b9e3b87964e..c7d7bf4f16f 100644 --- a/stan/math/prim/fun/binomial_coefficient_log.hpp +++ b/stan/math/prim/fun/binomial_coefficient_log.hpp @@ -84,7 +84,7 @@ inline return_type_t binomial_coefficient_log(const T_n n, } // Choosing the more stable of the symmetric branches - if (n > 0 && k > value_of_rec(n) / 2.0 + 1e-8) { + if (n > -1 && k > value_of_rec(n) / 2.0 + 1e-8) { return binomial_coefficient_log(n, n - k); } diff --git a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp index 907c49244a1..c70e37a2c72 100644 --- a/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp +++ b/test/unit/math/rev/fun/binomial_coefficient_log_test.cpp @@ -106,6 +106,8 @@ std::vector testValuesEdge = { {3, -1, stan::math::NEGATIVE_INFTY, -0.25, stan::math::INFTY}, {-1, -0.2, stan::math::INFTY, stan::math::NEGATIVE_INFTY, -4.324031329886049836}, + {-0.5, 0.5, stan::math::NEGATIVE_INFTY, stan::math::INFTY, + stan::math::NEGATIVE_INFTY}, {4.0, 5.0, stan::math::NEGATIVE_INFTY, stan::math::INFTY, stan::math::NEGATIVE_INFTY}, {1, 0, 0, 0, 1},