diff --git a/.drone.star b/.drone.star index fdb0db327d..089e13a959 100644 --- a/.drone.star +++ b/.drone.star @@ -16,11 +16,12 @@ windowsglobalimage="cppalliance/dronevs2019" def main(ctx): things_to_test = [ "special_fun", "distribution_tests", "mp", "misc", "interpolators", "quadrature", "autodiff", "long-running-tests", "float128_tests" ] - sanitizer_test = [ "special_fun", "distribution_tests", "misc", "interpolators", "quadrature", "autodiff", "float128_tests" ] - gnu_5_stds = [ "gnu++11" ] - gnu_6_stds = [ "gnu++11", "gnu++14" ] - gnu_8_stds = [ "gnu++11", "gnu++14", "gnu++17" ] - gnu_10_stds = [ "gnu++11", "gnu++14", "gnu++17", "gnu++20" ] + sanitizer_test = [ "special_fun", "distribution_tests", "misc", "interpolators", "quadrature", "float128_tests" ] + gnu_5_stds = [ "gnu++11", "c++11", "gnu++14", "c++14" ] + gnu_6_stds = [ "gnu++11", "c++11", "gnu++14", "c++14", "gnu++17", "c++17" ] + clang_6_stds = [ "c++11", "c++14", "c++17" ] + gnu_9_stds = [ "gnu++14", "c++14", "gnu++17", "c++17", "gnu++2a", "c++2a" ] + clang_10_stds = [ "c++14", "c++17", "c++2a" ] result = [] @@ -38,12 +39,18 @@ def main(ctx): result.append(linux_cxx("Ubunti g++-5 " + cxx + " " + suite, "g++-5", packages="g++-5", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-5', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_6_stds: result.append(linux_cxx("Ubunti g++-6 " + cxx + " " + suite, "g++-6", packages="g++-6", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-6', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - for cxx in gnu_8_stds: + result.append(linux_cxx("Ubunti g++-7 " + cxx + " " + suite, "g++-7", packages="g++-7", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-7', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubunti g++-8 " + cxx + " " + suite, "g++-8", packages="g++-8", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubunti g++-9 " + cxx + " " + suite, "g++-9", packages="g++-9", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + for cxx in clang_6_stds: + result.append(linux_cxx("Ubunti clang++-6 " + cxx + " " + suite, "clang++-6.0", packages="clang-6.0", llvm_os="xenial", llvm_ver="6.0", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-6.0', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubunti clang++-7 " + cxx + " " + suite, "clang++-7", packages="clang-7", llvm_os="xenial", llvm_ver="7", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-7', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubunti clang++-8 " + cxx + " " + suite, "clang++-8", packages="clang-8", llvm_os="xenial", llvm_ver="8", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubunti clang++-9 " + cxx + " " + suite, "clang++-9", packages="clang-9", llvm_os="xenial", llvm_ver="9", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - for cxx in gnu_10_stds: + for cxx in gnu_9_stds: result.append(linux_cxx("Ubunti g++-10 " + cxx + " " + suite, "g++-10", packages="g++-10", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubunti g++-11 " + cxx + " " + suite, "g++-11", packages="g++-11", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-11', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + for cxx in clang_10_stds: result.append(linux_cxx("Ubunti clang++-10 " + cxx + " " + suite, "clang++-10", packages="clang-10", llvm_os="xenial", llvm_ver="10", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) return result diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2c8b047f14..612fcebbd3 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -13,14 +13,14 @@ on: release: types: [published, created, edited] jobs: - ubuntu-focal: - runs-on: ubuntu-20.04 + ubuntu-jammy: + runs-on: ubuntu-22.04 strategy: fail-fast: false matrix: - compiler: [ g++-9, g++-11, clang++-10 ] - standard: [ c++11, c++14, c++17, c++2a ] - suite: [ float128_tests, special_fun, distribution_tests, misc, quadrature, mp, interpolators, autodiff, ../example//examples, ../tools ] + compiler: [ g++-12, clang++-14 ] + standard: [ c++14, c++17, c++20 ] + suite: [ github_ci_block_1, github_ci_block_2 ] steps: - uses: actions/checkout@v2 with: @@ -47,7 +47,7 @@ jobs: if: steps.retry1.outcome=='failure' run: sudo apt-add-repository -y "ppa:ubuntu-toolchain-r/test" - name: Install packages - run: sudo apt install g++-9 g++-11 clang-9 clang-10 libgmp-dev libmpfr-dev libfftw3-dev + run: sudo apt install g++-12 clang-14 libgmp-dev libmpfr-dev libfftw3-dev - name: Checkout main boost run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root - name: Update tools/boostdep @@ -75,16 +75,15 @@ jobs: run: ./config_info_travis working-directory: ../boost-root/libs/config/test - name: Test - run: ../../../b2 toolset=$TOOLSET ${{ matrix.suite }} define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER + run: ../../../b2 -j2 toolset=$TOOLSET ${{ matrix.suite }} define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER working-directory: ../boost-root/libs/math/test - ubuntu-bionic: - runs-on: ubuntu-18.04 + ubuntu-focal-no-eh: + runs-on: ubuntu-20.04 strategy: fail-fast: false matrix: - compiler: [ g++-6, clang++-6.0, g++-7, g++-8, clang++-7, clang++-8 ] - standard: [ c++11, c++14, c++17 ] - suite: [ float128_tests, special_fun, distribution_tests, misc, quadrature, mp, interpolators, autodiff, ../example//examples, ../tools ] + compiler: [ g++-9, g++-11, clang++-10 ] + standard: [ c++11, c++14, c++17, c++2a ] steps: - uses: actions/checkout@v2 with: @@ -111,7 +110,7 @@ jobs: if: steps.retry1.outcome=='failure' run: sudo apt-add-repository -y "ppa:ubuntu-toolchain-r/test" - name: Install packages - run: sudo apt install g++-6 g++-7 g++-8 clang-6.0 clang-7 clang-8 libgmp-dev libmpfr-dev libfftw3-dev + run: sudo apt install g++-9 g++-11 clang-9 clang-10 libgmp-dev libmpfr-dev libfftw3-dev - name: Checkout main boost run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root - name: Update tools/boostdep @@ -139,7 +138,7 @@ jobs: run: ./config_info_travis working-directory: ../boost-root/libs/config/test - name: Test - run: ../../../b2 toolset=$TOOLSET ${{ matrix.suite }} define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER + run: ../../../b2 toolset=$TOOLSET no_eh_tests exception-handling=off rtti=off define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER working-directory: ../boost-root/libs/math/test macos: runs-on: macos-latest @@ -147,8 +146,8 @@ jobs: fail-fast: false matrix: toolset: [ clang ] - standard: [ 11, 14, 17, 2a ] - suite: [ float128_tests, special_fun, distribution_tests, misc, quadrature, mp, interpolators, autodiff, ../example//examples, ../tools ] + standard: [ 11, 14, 17, 20 ] + suite: [ github_ci_block_1, github_ci_block_2 ] steps: - uses: actions/checkout@v2 with: @@ -194,8 +193,56 @@ jobs: strategy: fail-fast: false matrix: - toolset: [ gcc, msvc-14.0, msvc-14.2 ] - standard: [ 11, 14, 17 ] + toolset: [ msvc-14.0, msvc-14.2 ] + standard: [ 14, 17 ] + suite: [ github_ci_block_1, github_ci_block_2 ] + steps: + - uses: actions/checkout@v2 + with: + fetch-depth: '0' + - uses: mstachniuk/ci-skip@v1 + with: + commit-filter: '[skip ci];[ci skip];[CI SKIP];[SKIP CI];***CI SKIP***;***SKIP CI***;[apple];[Apple];[APPLE];[linux];[Linux];[LINUX];[standalone];[STANDALONE];[cygwin];[CYGWIN]' + commit-filter-separator: ';' + fail-fast: true + - name: Checkout main boost + run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root + - name: Update tools/boostdep + run: git submodule update --init tools/boostdep + working-directory: ../boost-root + - name: Copy files + run: xcopy /s /e /q %GITHUB_WORKSPACE% libs\math + working-directory: ../boost-root + - name: Install deps + run: python tools/boostdep/depinst/depinst.py math + working-directory: ../boost-root + - name: Bootstrap + run: bootstrap + working-directory: ../boost-root + - name: Generate headers + run: b2 headers + working-directory: ../boost-root + - name: Config info install + run: ..\..\..\b2 config_info_travis_install %ARGS% + working-directory: ../boost-root/libs/config/test + - name: Config info + run: config_info_travis + working-directory: ../boost-root/libs/config/test + - name: Test + run: ..\..\..\b2 --hash %ARGS% define=CI_SUPPRESS_KNOWN_ISSUES ${{ matrix.suite }} + working-directory: ../boost-root/libs/math/test + windows_gcc: + runs-on: windows-2019 + defaults: + run: + shell: cmd + env: + ARGS: toolset=${{ matrix.toolset }} address-model=64 cxxstd=${{ matrix.standard }} + strategy: + fail-fast: false + matrix: + toolset: [ gcc ] + standard: [ 14, 17 ] suite: [ float128_tests, special_fun, distribution_tests, misc, quadrature, mp, interpolators, autodiff, ../example//examples, ../tools ] steps: - uses: actions/checkout@v2 @@ -243,7 +290,7 @@ jobs: fail-fast: false matrix: standard: [ 14, 17, 20 ] - suite: [ float128_tests, special_fun, distribution_tests, misc, quadrature, mp, interpolators, autodiff, ../example//examples, ../tools ] + suite: [ github_ci_block_1, github_ci_block_2 ] steps: - uses: actions/checkout@v2 with: @@ -415,8 +462,8 @@ jobs: fail-fast: false matrix: compiler: [ g++-10 ] - standard: [ c++11, c++14, c++17, c++2a ] - suite: [ float128_tests, special_fun, distribution_tests, misc, quadrature, interpolators, autodiff, ../example//examples, ../tools ] + standard: [ c++14, c++17, c++20 ] + suite: [ github_ci_block_1, github_ci_block_2 ] steps: - uses: actions/checkout@v2 with: diff --git a/.gitignore b/.gitignore index 901283a3f7..5ff21efb64 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,6 @@ Makefile **CTestTestfile.cmake DartConfiguration.tcl cmake-build-debug/* +.cmake/* +build.ninja +.ninja* diff --git a/README.md b/README.md index d6f6b920f6..b9dba9a9a9 100644 --- a/README.md +++ b/README.md @@ -2,8 +2,8 @@ Boost Math Library [![Build Status](https://drone.cpp.al/api/badges/boostorg/math/status.svg)](https://drone.cpp.al/boostorg/math)[![Build Status](https://github.com/boostorg/math/workflows/CI/badge.svg?branch=develop)](https://github.com/boostorg/math/actions) ================== ->ANNOUNCEMENT: Support for C++03 is now deprecated in this library and will be supported in existing features ->only until March 2021. New features will require *at least* C++11, as will existing features from next year. +>ANNOUNCEMENT: Support for C++11 will be deprecated in this library starting in July 2023 (Boost 1.82). +>New features will require *at least* C++14, as will existing features starting with the deprecation release. This library is divided into several interconnected parts: diff --git a/doc/sf/ccmath.qbk b/doc/sf/ccmath.qbk index 0a88ce6310..4a2b0b2fbf 100644 --- a/doc/sf/ccmath.qbk +++ b/doc/sf/ccmath.qbk @@ -182,6 +182,13 @@ All of the following functions require C++17 or greater. template inline constexpr bool isunordered(T x, T y) noexcept + template + inline constexpr Real fma(Real x, Real y, Real z) noexcept + Requires compiling with fma flag + + template + inline constepxr Promoted fma(Arithmetic1 x, Arithmetic2 y, Arithmetic3 z) noexcept + } // Namespaces [endsect] [/section:ccmath Constexpr CMath] diff --git a/doc/statistics/chatterjee_correlation.qbk b/doc/statistics/chatterjee_correlation.qbk new file mode 100644 index 0000000000..520ab0cc6b --- /dev/null +++ b/doc/statistics/chatterjee_correlation.qbk @@ -0,0 +1,74 @@ +[/ + Copyright 2022 Matt Borland + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] + +[section:chatterjee_correlation Chatterjee Correlation] + +[heading Synopsis] + +`` +#include + +namespace boost::math::statistics { + + C++17: + template + auto chatterjee_correlation(ExecutionPolicy&& exec, const Container& u, const Container& v); + + C++11: + template + auto chatterjee_correlation(const Container& u, const Container& v); +} +`` + +[heading Description] + +The classical correlation coefficients like the Pearson's correlation are useful primarily for distinguishing when one dataset depends linearly on another. +However, Pearson's correlation coefficient has a known weakness in that when the dependent variable has an obvious functional relationship with the independent variable, the value of the correlation coefficient can take on any value. +As Chatterjee says: + +> Ideally, one would like a coefficient that approaches +its maximum value if and only if one variable looks more and more like a +noiseless function of the other, just as Pearson correlation is close to its maximum value if and only if one variable is close to being a noiseless linear function of the other. + +This is the problem Chatterjee's coefficient solves. +Let X and Y be random variables, where Y is not constant, and let (X_i, Y_i) be samples from this distribution. +Rearrange these samples so that X_(0) < X_{(1)} < ... X_{(n-1)} and create (R(X_{(i)}), R(Y_{(i)})). +The Chatterjee correlation is then given by + +[$../equations/chatterjee_correlation.svg] + +In the limit of an infinite amount of i.i.d data, the statistic lies in [0, 1]. +However, if the data is not infinite, the statistic may be negative. +If X and Y are independent, the value is zero, and if Y is a measurable function of X, then the statistic is unity. +The complexity is O(n log n). + +An example is given below: + + std::vector X{1,2,3,4,5}; + std::vector Y{1,2,3,4,5}; + using boost::math::statistics::chatterjee_correlation; + double coeff = chatterjee_correlation(X, Y); + +The implementation follows [@https://arxiv.org/pdf/1909.10140.pdf Chatterjee's paper]. + +/Nota bene:/ If the input is an integer type the output will be a double precision type. + +[heading Invariants] + +The function expects at least two samples, a non-constant vector Y, and the same number of X's as Y's. +If Y is constant, the result is a quiet NaN. +The data set must be sorted by X values. +If there are ties in the values of X, then the statistic is random due to the random breaking of ties. +Of course, random numbers are not used internally, but the result is not guaranteed to be identical on different systems. + +[heading References] + +* Chatterjee, Sourav. "A new coefficient of correlation." Journal of the American Statistical Association 116.536 (2021): 2009-2022. + +[endsect] +[/section:chatterjee_correlation Chatterjee Correlation] diff --git a/example/dot_net_example/distribution_explorer/readme.txt b/example/dot_net_example/distribution_explorer/readme.txt index d483173d73..04fec50f21 100644 --- a/example/dot_net_example/distribution_explorer/readme.txt +++ b/example/dot_net_example/distribution_explorer/readme.txt @@ -3,7 +3,7 @@ Statistical Distribution Explorer Paul A. Bristow John Maddock -Copyright © 2008 , 2009, 2010, 2012 Paul A. Bristow, John Maddock +Copyright (C) 2008, 2009, 2010, 2012 Paul A. Bristow, John Maddock Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) diff --git a/example/policy_eg_1.cpp b/example/policy_eg_1.cpp index 7f8d9fbb44..01f2f2797c 100644 --- a/example/policy_eg_1.cpp +++ b/example/policy_eg_1.cpp @@ -21,7 +21,7 @@ using boost::math::policies::policy; using boost::math::policies::evaluation_error; using boost::math::policies::domain_error; using boost::math::policies::overflow_error; -using boost::math::policies::domain_error; +using boost::math::policies::underflow_error; using boost::math::policies::pole_error; // Actions on error (in enum error_policy_type): using boost::math::policies::errno_on_error; diff --git a/include/boost/math/ccmath/ccmath.hpp b/include/boost/math/ccmath/ccmath.hpp index 72c49922f1..2749ec7b28 100644 --- a/include/boost/math/ccmath/ccmath.hpp +++ b/include/boost/math/ccmath/ccmath.hpp @@ -38,5 +38,6 @@ #include #include #include +#include #endif // BOOST_MATH_CCMATH_HPP diff --git a/include/boost/math/ccmath/fma.hpp b/include/boost/math/ccmath/fma.hpp new file mode 100644 index 0000000000..3056f76d47 --- /dev/null +++ b/include/boost/math/ccmath/fma.hpp @@ -0,0 +1,128 @@ +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_CCMATH_FMA_HPP +#define BOOST_MATH_CCMATH_FMA_HPP + +#include +#include +#include +#include +#include +#include + +namespace boost::math::ccmath { + +namespace detail { + +template +constexpr T fma_imp(const T x, const T y, const T z) noexcept +{ + #if defined(__GNUC__) && !defined(__clang__) && !defined(__INTEL_COMPILER) && !defined(__INTEL_LLVM_COMPILER) + if constexpr (std::is_same_v) + { + return __builtin_fmaf(x, y, z); + } + else if constexpr (std::is_same_v) + { + return __builtin_fma(x, y, z); + } + else if constexpr (std::is_same_v) + { + return __builtin_fmal(x, y, z); + } + #endif + + // If we can't use compiler intrinsics hope that -fma flag optimizes this call to fma instruction + return (x * y) + z; +} + +} // Namespace detail + +template , bool> = true> +constexpr Real fma(Real x, Real y, Real z) noexcept +{ + if (BOOST_MATH_IS_CONSTANT_EVALUATED(x)) + { + if (x == 0 && boost::math::ccmath::isinf(y)) + { + return std::numeric_limits::quiet_NaN(); + } + else if (y == 0 && boost::math::ccmath::isinf(x)) + { + return std::numeric_limits::quiet_NaN(); + } + else if (boost::math::ccmath::isnan(x)) + { + return std::numeric_limits::quiet_NaN(); + } + else if (boost::math::ccmath::isnan(y)) + { + return std::numeric_limits::quiet_NaN(); + } + else if (boost::math::ccmath::isnan(z)) + { + return std::numeric_limits::quiet_NaN(); + } + + return boost::math::ccmath::detail::fma_imp(x, y, z); + } + else + { + using std::fma; + return fma(x, y, z); + } +} + +template +constexpr auto fma(T1 x, T2 y, T3 z) noexcept +{ + if (BOOST_MATH_IS_CONSTANT_EVALUATED(x)) + { + // If the type is an integer (e.g. epsilon == 0) then set the epsilon value to 1 so that type is at a minimum + // cast to double + constexpr auto T1p = std::numeric_limits::epsilon() > 0 ? std::numeric_limits::epsilon() : 1; + constexpr auto T2p = std::numeric_limits::epsilon() > 0 ? std::numeric_limits::epsilon() : 1; + constexpr auto T3p = std::numeric_limits::epsilon() > 0 ? std::numeric_limits::epsilon() : 1; + + using promoted_type = + #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + std::conditional_t>>>>>; + #else + >>>; + #endif + + return boost::math::ccmath::fma(promoted_type(x), promoted_type(y), promoted_type(z)); + } + else + { + using std::fma; + return fma(x, y, z); + } +} + +constexpr float fmaf(float x, float y, float z) noexcept +{ + return boost::math::ccmath::fma(x, y, z); +} + +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS +constexpr long double fmal(long double x, long double y, long double z) noexcept +{ + return boost::math::ccmath::fma(x, y, z); +} +#endif + +} // Namespace boost::math::ccmath + +#endif // BOOST_MATH_CCMATH_FMA_HPP diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 5cb2c05f9b..a8fcbbc05f 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -29,6 +29,7 @@ #ifndef BOOST_MATH_DIST_ARCSINE_HPP #define BOOST_MATH_DIST_ARCSINE_HPP +#include #include #include // complements. #include // error checks. diff --git a/include/boost/math/distributions/chi_squared.hpp b/include/boost/math/distributions/chi_squared.hpp index e97bee7ce7..7b65a0da68 100644 --- a/include/boost/math/distributions/chi_squared.hpp +++ b/include/boost/math/distributions/chi_squared.hpp @@ -23,10 +23,10 @@ template > class chi_squared_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - chi_squared_distribution(RealType i) : m_df(i) + explicit chi_squared_distribution(RealType i) : m_df(i) { RealType result; detail::check_df( @@ -53,7 +53,7 @@ class chi_squared_distribution RealType m_df; // degrees of freedom is a positive real number. }; // class chi_squared_distribution -typedef chi_squared_distribution chi_squared; +using chi_squared = chi_squared_distribution; #ifdef __cpp_deduction_guides template @@ -66,7 +66,7 @@ chi_squared_distribution(RealType)->chi_squared_distribution -inline const std::pair range(const chi_squared_distribution& /*dist*/) +inline std::pair range(const chi_squared_distribution& /*dist*/) { // Range of permissible values for random variable x. if (std::numeric_limits::has_infinity) { @@ -84,7 +84,7 @@ inline const std::pair range(const chi_squared_distribution< #endif template -inline const std::pair support(const chi_squared_distribution& /*dist*/) +inline std::pair support(const chi_squared_distribution& /*dist*/) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. return std::pair(static_cast(0), tools::max_value()); // 0 to + infinity. @@ -224,17 +224,6 @@ inline RealType mode(const chi_squared_distribution& dist) { RealType df = dist.degrees_of_freedom(); static const char* function = "boost::math::mode(const chi_squared_distribution<%1%>&)"; - // Most sources only define mode for df >= 2, - // but for 0 <= df <= 2, the pdf maximum actually occurs at random variate = 0; - // So one could extend the definition of mode thus: - //if(df < 0) - //{ - // return policies::raise_domain_error( - // function, - // "Chi-Squared distribution only has a mode for degrees of freedom >= 0, but got degrees of freedom = %1%.", - // df, Policy()); - //} - //return (df <= 2) ? 0 : df - 2; if(df < 2) return policies::raise_domain_error( @@ -244,25 +233,12 @@ inline RealType mode(const chi_squared_distribution& dist) return df - 2; } -//template -//inline RealType median(const chi_squared_distribution& dist) -//{ // Median is given by Quantile[dist, 1/2] -// RealType df = dist.degrees_of_freedom(); -// if(df <= 1) -// return tools::domain_error( -// BOOST_CURRENT_FUNCTION, -// "The Chi-Squared distribution only has a mode for degrees of freedom >= 2, but got degrees of freedom = %1%.", -// df); -// return df - RealType(2)/3; -//} -// Now implemented via quantile(half) in derived accessors. - template inline RealType skewness(const chi_squared_distribution& dist) { BOOST_MATH_STD_USING // For ADL RealType df = dist.degrees_of_freedom(); - return sqrt (8 / df); // == 2 * sqrt(2 / df); + return sqrt (8 / df); } template diff --git a/include/boost/math/distributions/detail/derived_accessors.hpp b/include/boost/math/distributions/detail/derived_accessors.hpp index 9101c8ae1f..e2eca511bf 100644 --- a/include/boost/math/distributions/detail/derived_accessors.hpp +++ b/include/boost/math/distributions/detail/derived_accessors.hpp @@ -110,6 +110,13 @@ inline typename Distribution::value_type pdf(const Distribution& dist, const Rea return pdf(dist, static_cast(x)); } template +inline typename Distribution::value_type logpdf(const Distribution& dist, const RealType& x) +{ + using std::log; + typedef typename Distribution::value_type value_type; + return log(pdf(dist, static_cast(x))); +} +template inline typename Distribution::value_type cdf(const Distribution& dist, const RealType& x) { typedef typename Distribution::value_type value_type; diff --git a/include/boost/math/distributions/exponential.hpp b/include/boost/math/distributions/exponential.hpp index ba7eae927f..5214575a64 100644 --- a/include/boost/math/distributions/exponential.hpp +++ b/include/boost/math/distributions/exponential.hpp @@ -60,10 +60,10 @@ template > class exponential_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - exponential_distribution(RealType l_lambda = 1) + explicit exponential_distribution(RealType l_lambda = 1) : m_lambda(l_lambda) { RealType err; @@ -76,7 +76,7 @@ class exponential_distribution RealType m_lambda; }; -typedef exponential_distribution exponential; +using exponential = exponential_distribution; #ifdef __cpp_deduction_guides template @@ -84,7 +84,7 @@ exponential_distribution(RealType)->exponential_distribution -inline const std::pair range(const exponential_distribution& /*dist*/) +inline std::pair range(const exponential_distribution& /*dist*/) { // Range of permissible values for random variable x. if (std::numeric_limits::has_infinity) { @@ -98,7 +98,7 @@ inline const std::pair range(const exponential_distribution< } template -inline const std::pair support(const exponential_distribution& /*dist*/) +inline std::pair support(const exponential_distribution& /*dist*/) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -127,6 +127,24 @@ inline RealType pdf(const exponential_distribution& dist, cons return result; } // pdf +template +inline RealType logpdf(const exponential_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::logpdf(const exponential_distribution<%1%>&, %1%)"; + + RealType lambda = dist.lambda(); + RealType result = -std::numeric_limits::infinity(); + if(0 == detail::verify_lambda(function, lambda, &result, Policy())) + return result; + if(0 == detail::verify_exp_x(function, x, &result, Policy())) + return result; + + result = log(lambda) - lambda * x; + return result; +} // logpdf + template inline RealType cdf(const exponential_distribution& dist, const RealType& x) { diff --git a/include/boost/math/distributions/extreme_value.hpp b/include/boost/math/distributions/extreme_value.hpp index d503b31bc9..b2e8bea966 100644 --- a/include/boost/math/distributions/extreme_value.hpp +++ b/include/boost/math/distributions/extreme_value.hpp @@ -53,10 +53,10 @@ template > class extreme_value_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - extreme_value_distribution(RealType a = 0, RealType b = 1) + explicit extreme_value_distribution(RealType a = 0, RealType b = 1) : m_a(a), m_b(b) { RealType err; @@ -68,10 +68,11 @@ class extreme_value_distribution RealType scale()const { return m_b; } private: - RealType m_a, m_b; + RealType m_a; + RealType m_b; }; -typedef extreme_value_distribution extreme_value; +using extreme_value = extreme_value_distribution; #ifdef __cpp_deduction_guides template @@ -81,7 +82,7 @@ extreme_value_distribution(RealType,RealType)->extreme_value_distribution -inline const std::pair range(const extreme_value_distribution& /*dist*/) +inline std::pair range(const extreme_value_distribution& /*dist*/) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair( @@ -90,7 +91,7 @@ inline const std::pair range(const extreme_value_distributio } template -inline const std::pair support(const extreme_value_distribution& /*dist*/) +inline std::pair support(const extreme_value_distribution& /*dist*/) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -122,6 +123,31 @@ inline RealType pdf(const extreme_value_distribution& dist, co return result; } // pdf +template +inline RealType logpdf(const extreme_value_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::logpdf(const extreme_value_distribution<%1%>&, %1%)"; + + RealType a = dist.location(); + RealType b = dist.scale(); + RealType result = -std::numeric_limits::infinity(); + if(0 == detail::verify_scale_b(function, b, &result, Policy())) + return result; + if(0 == detail::check_finite(function, a, &result, Policy())) + return result; + if((boost::math::isinf)(x)) + return 0.0f; + if(0 == detail::check_x(function, x, &result, Policy())) + return result; + RealType e = (a - x) / b; + if(e < tools::log_max_value()) + result = log(1/b) + e - exp(e); + // else.... result *must* be zero since exp(e) is infinite... + return result; +} // logpdf + template inline RealType cdf(const extreme_value_distribution& dist, const RealType& x) { diff --git a/include/boost/math/distributions/gamma.hpp b/include/boost/math/distributions/gamma.hpp index 8a3414d312..28b7c55b0b 100644 --- a/include/boost/math/distributions/gamma.hpp +++ b/include/boost/math/distributions/gamma.hpp @@ -17,6 +17,7 @@ #include #include +#include namespace boost{ namespace math { @@ -71,10 +72,10 @@ template > class gamma_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - gamma_distribution(RealType l_shape, RealType l_scale = 1) + explicit gamma_distribution(RealType l_shape, RealType l_scale = 1) : m_shape(l_shape), m_scale(l_scale) { RealType result; @@ -108,14 +109,14 @@ gamma_distribution(RealType,RealType)->gamma_distribution -inline const std::pair range(const gamma_distribution& /* dist */) +inline std::pair range(const gamma_distribution& /* dist */) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(0), max_value()); } template -inline const std::pair support(const gamma_distribution& /* dist */) +inline std::pair support(const gamma_distribution& /* dist */) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -147,6 +148,33 @@ inline RealType pdf(const gamma_distribution& dist, const Real return result; } // pdf +template +inline RealType logpdf(const gamma_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + using boost::math::lgamma; + + static const char* function = "boost::math::logpdf(const gamma_distribution<%1%>&, %1%)"; + + RealType k = dist.shape(); + RealType theta = dist.scale(); + + RealType result = -std::numeric_limits::infinity(); + if(false == detail::check_gamma(function, theta, k, &result, Policy())) + return result; + if(false == detail::check_gamma_x(function, x, &result, Policy())) + return result; + + if(x == 0) + { + return std::numeric_limits::quiet_NaN(); + } + + result = -k*log(theta) + (k-1)*log(x) - lgamma(k) - (x/theta); + + return result; +} // logpdf + template inline RealType cdf(const gamma_distribution& dist, const RealType& x) { diff --git a/include/boost/math/distributions/inverse_gamma.hpp b/include/boost/math/distributions/inverse_gamma.hpp index 9266bc22f6..8c9e4763d5 100644 --- a/include/boost/math/distributions/inverse_gamma.hpp +++ b/include/boost/math/distributions/inverse_gamma.hpp @@ -28,6 +28,7 @@ #include #include +#include namespace boost{ namespace math { @@ -88,10 +89,10 @@ template > class inverse_gamma_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1) + explicit inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1) : m_shape(l_shape), m_scale(l_scale) { RealType result; @@ -117,10 +118,9 @@ class inverse_gamma_distribution RealType m_scale; // distribution scale }; -typedef inverse_gamma_distribution inverse_gamma; +using inverse_gamma = inverse_gamma_distribution; // typedef - but potential clash with name of inverse gamma *function*. -// but there is a typedef for gamma -// typedef boost::math::gamma_distribution gamma; +// but there is a typedef for the gamma distribution (gamma) #ifdef __cpp_deduction_guides template @@ -132,14 +132,14 @@ inverse_gamma_distribution(RealType,RealType)->inverse_gamma_distribution -inline const std::pair range(const inverse_gamma_distribution& /* dist */) +inline std::pair range(const inverse_gamma_distribution& /* dist */) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(0), max_value()); } template -inline const std::pair support(const inverse_gamma_distribution& /* dist */) +inline std::pair support(const inverse_gamma_distribution& /* dist */) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -190,11 +190,47 @@ inline RealType pdf(const inverse_gamma_distribution& dist, co } result /= (x * x); } - // better than naive - // result = (pow(scale, shape) * pow(x, (-shape -1)) * exp(-scale/x) ) / tgamma(shape); + return result; } // pdf +template +inline RealType logpdf(const inverse_gamma_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + using boost::math::lgamma; + + static const char* function = "boost::math::logpdf(const inverse_gamma_distribution<%1%>&, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = -std::numeric_limits::infinity(); + if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy())) + { // distribution parameters bad. + return result; + } + if(x == 0) + { // Treat random variate zero as a special case. + return result; + } + else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy())) + { // x bad. + return result; + } + result = scale / x; + if(result < tools::min_value()) + return result; // random variable is infinite or so close as to make no difference. + + // x * x may under or overflow, likewise our result + if (!(boost::math::isfinite)(x*x)) + { + return policies::raise_overflow_error(function, "PDF is infinite.", Policy()); + } + + return shape * log(scale) + (-shape-1)*log(x) - lgamma(shape) - (scale/x); +} // pdf + template inline RealType cdf(const inverse_gamma_distribution& dist, const RealType& x) { @@ -269,7 +305,6 @@ inline RealType cdf(const complemented2_type #include #include // for gamma function -// using boost::math::gamma_p; #include -//using std::tr1::tuple; -//using std::tr1::make_tuple; #include -//using boost::math::tools::newton_raphson_iterate; #include @@ -71,10 +67,10 @@ template > class inverse_gaussian_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - inverse_gaussian_distribution(RealType l_mean = 1, RealType l_scale = 1) + explicit inverse_gaussian_distribution(RealType l_mean = 1, RealType l_scale = 1) : m_mean(l_mean), m_scale(l_scale) { // Default is a 1,1 inverse_gaussian distribution. static const char* function = "boost::math::inverse_gaussian_distribution<%1%>::inverse_gaussian_distribution"; @@ -113,7 +109,7 @@ class inverse_gaussian_distribution RealType m_scale; // distribution standard deviation or scale, aka lambda. }; // class normal_distribution -typedef inverse_gaussian_distribution inverse_gaussian; +using inverse_gaussian = inverse_gaussian_distribution; #ifdef __cpp_deduction_guides template @@ -123,14 +119,14 @@ inverse_gaussian_distribution(RealType,RealType)->inverse_gaussian_distribution< #endif template -inline const std::pair range(const inverse_gaussian_distribution& /*dist*/) +inline std::pair range(const inverse_gaussian_distribution& /*dist*/) { // Range of permissible values for random variable x, zero to max. using boost::math::tools::max_value; return std::pair(static_cast(0.), max_value()); // - to + max value. } template -inline const std::pair support(const inverse_gaussian_distribution& /*dist*/) +inline std::pair support(const inverse_gaussian_distribution& /*dist*/) { // Range of supported values for random variable x, zero to max. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -174,6 +170,43 @@ inline RealType pdf(const inverse_gaussian_distribution& dist, return result; } // pdf +template +inline RealType logpdf(const inverse_gaussian_distribution& dist, const RealType& x) +{ // Probability Density Function + BOOST_MATH_STD_USING // for ADL of std functions + + RealType scale = dist.scale(); + RealType mean = dist.mean(); + RealType result = -std::numeric_limits::infinity(); + static const char* function = "boost::math::logpdf(const inverse_gaussian_distribution<%1%>&, %1%)"; + if(false == detail::check_scale(function, scale, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, mean, &result, Policy())) + { + return result; + } + if(false == detail::check_x_gt0(function, mean, &result, Policy())) + { + return result; + } + if(false == detail::check_positive_x(function, x, &result, Policy())) + { + return result; + } + + if (x == 0) + { + return std::numeric_limits::quiet_NaN(); // Convenient, even if not defined mathematically. log(0) + } + + const RealType two_pi = boost::math::constants::two_pi(); + + result = (-scale*pow(mean - x, RealType(2))/(mean*mean*x) + log(scale) - 3*log(x) - log(two_pi)) / 2; + return result; +} // pdf + template inline RealType cdf(const inverse_gaussian_distribution& dist, const RealType& x) { // Cumulative Density Function. @@ -203,10 +236,7 @@ inline RealType cdf(const inverse_gaussian_distribution& dist, { return 0; // Convenient, even if not defined mathematically. } - // Problem with this formula for large scale > 1000 or small x, - //result = 0.5 * (erf(sqrt(scale / x) * ((x / mean) - 1) / constants::root_two(), Policy()) + 1) - // + exp(2 * scale / mean) / 2 - // * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two(), Policy())); + // Problem with this formula for large scale > 1000 or small x // so use normal distribution version: // Wikipedia CDF equation http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution. @@ -270,22 +300,20 @@ namespace detail template inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1) { // guess at random variate value x for inverse gaussian quantile. - BOOST_MATH_STD_USING - using boost::math::policies::policy; - // Error type. - using boost::math::policies::overflow_error; - // Action. - using boost::math::policies::ignore_error; + BOOST_MATH_STD_USING + using boost::math::policies::policy; + // Error type. + using boost::math::policies::overflow_error; + // Action. + using boost::math::policies::ignore_error; - typedef policy< - overflow_error // Ignore overflow (return infinity) - > no_overthrow_policy; + using no_overthrow_policy = policy>; RealType x; // result is guess at random variate value x. RealType phi = lambda / mu; if (phi > 2.) { // Big phi, so starting to look like normal Gaussian distribution. - // x=(qnorm(p,0,1,true,false) - 0.5 * sqrt(mu/lambda)) / sqrt(lambda/mu); + // // Whitmore, G.A. and Yalovsky, M. // A normalising logarithmic transformation for inverse Gaussian random variables, // Technometrics 20-2, 207-208 (1978), but using expression from @@ -300,14 +328,12 @@ namespace detail using boost::math::gamma_distribution; // Define the distribution, using gamma_nooverflow: - typedef gamma_distribution gamma_nooverflow; + using gamma_nooverflow = gamma_distribution; gamma_nooverflow g(static_cast(0.5), static_cast(1.)); - // gamma_nooverflow g(static_cast(0.5), static_cast(1.)); - // R qgamma(0.2, 0.5, 1) 0.0320923 + // R qgamma(0.2, 0.5, 1) = 0.0320923 RealType qg = quantile(complement(g, p)); - //RealType qg1 = qgamma(1.- p, 0.5, 1.0, true, false); x = lambda / (qg * 2); // if (x > mu/2) // x > mu /2? @@ -350,7 +376,7 @@ inline RealType quantile(const inverse_gaussian_distribution& { // overflow result = policies::raise_overflow_error(function, "probability parameter is 1, but must be < 1!", Policy()); - return result; // std::numeric_limits::infinity(); + return result; // infinity; } RealType guess = detail::guess_ig(p, dist.mean(), dist.scale()); @@ -378,21 +404,7 @@ inline RealType cdf(const complemented2_type::has_infinity && x == std::numeric_limits::infinity()) - //{ // cdf complement +infinity is zero. - // return 0; - //} - //if(std::numeric_limits::has_infinity && x == -std::numeric_limits::infinity()) - //{ // cdf complement -infinity is unity. - // return 1; - //} + RealType result = 0; if(false == detail::check_scale(function, scale, &result, Policy())) return result; diff --git a/include/boost/math/distributions/laplace.hpp b/include/boost/math/distributions/laplace.hpp index 9b268f3f33..cc922a879a 100644 --- a/include/boost/math/distributions/laplace.hpp +++ b/include/boost/math/distributions/laplace.hpp @@ -36,13 +36,13 @@ class laplace_distribution // ---------------------------------- // public Types // ---------------------------------- - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; // ---------------------------------- // Constructor(s) // ---------------------------------- - laplace_distribution(RealType l_location = 0, RealType l_scale = 1) + explicit laplace_distribution(RealType l_location = 0, RealType l_scale = 1) : m_location(l_location), m_scale(l_scale) { RealType result; @@ -78,7 +78,7 @@ class laplace_distribution // // Convenient type synonym for double. -typedef laplace_distribution laplace; +using laplace = laplace_distribution; #ifdef __cpp_deduction_guides template @@ -90,9 +90,9 @@ laplace_distribution(RealType,RealType)->laplace_distribution -inline const std::pair range(const laplace_distribution&) +inline std::pair range(const laplace_distribution&) { - if (std::numeric_limits::has_infinity) + if (std::numeric_limits::has_infinity) { // Can use infinity. return std::pair(-std::numeric_limits::infinity(), std::numeric_limits::infinity()); // - to + infinity. } @@ -105,7 +105,7 @@ inline const std::pair range(const laplace_distribution -inline const std::pair support(const laplace_distribution&) +inline std::pair support(const laplace_distribution&) { if (std::numeric_limits::has_infinity) { // Can Use infinity. @@ -150,6 +150,48 @@ inline RealType pdf(const laplace_distribution& dist, const Re return result; } // pdf +template +inline RealType logpdf(const laplace_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + // Checking function argument + RealType result = -std::numeric_limits::infinity(); + const char* function = "boost::math::logpdf(const laplace_distribution<%1%>&, %1%))"; + + // Check scale and location. + if (false == dist.check_parameters(function, &result)) + { + return result; + } + // Special pdf values. + if((boost::math::isinf)(x)) + { + return result; // pdf + and - infinity is zero so logpdf is -INF + } + if (false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + const RealType mu = dist.scale(); + const RealType b = dist.location(); + + // if b is 0 avoid divde by 0 error + if(abs(b) < std::numeric_limits::epsilon()) + { + result = log(pdf(dist, x)); + } + else + { + // General case + const RealType log2 = boost::math::constants::ln_two(); + result = -abs(x-mu)/b - log(b) - log2; + } + + return result; +} // logpdf + template inline RealType cdf(const laplace_distribution& dist, const RealType& x) { @@ -201,14 +243,14 @@ inline RealType quantile(const laplace_distribution& dist, con { result = policies::raise_overflow_error(function, "probability parameter is 0, but must be > 0!", Policy()); - return -result; // -std::numeric_limits::infinity(); + return -result; // -inf } if(p == 1) { result = policies::raise_overflow_error(function, "probability parameter is 1, but must be < 1!", Policy()); - return result; // std::numeric_limits::infinity(); + return result; // inf } // Calculate Quantile RealType scale( dist.scale() ); @@ -238,8 +280,6 @@ inline RealType cdf(const complemented2_type #include +#include namespace boost{ namespace math{ @@ -28,10 +29,10 @@ template > class normal_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - normal_distribution(RealType l_mean = 0, RealType sd = 1) + explicit normal_distribution(RealType l_mean = 0, RealType sd = 1) : m_mean(l_mean), m_sd(sd) { // Default is a 'standard' normal distribution N01. static const char* function = "boost::math::normal_distribution<%1%>::normal_distribution"; @@ -69,7 +70,7 @@ class normal_distribution RealType m_sd; // distribution standard deviation or scale. }; // class normal_distribution -typedef normal_distribution normal; +using normal = normal_distribution; // // Deduction guides, note we don't check the @@ -91,7 +92,7 @@ normal_distribution(RealType)->normal_distribution -inline const std::pair range(const normal_distribution& /*dist*/) +inline std::pair range(const normal_distribution& /*dist*/) { // Range of permissible values for random variable x. if (std::numeric_limits::has_infinity) { @@ -105,7 +106,7 @@ inline const std::pair range(const normal_distribution -inline const std::pair support(const normal_distribution& /*dist*/) +inline std::pair support(const normal_distribution& /*dist*/) { // This is range values for random variable x where cdf rises from 0 to 1, and outside it, the pdf is zero. if (std::numeric_limits::has_infinity) { @@ -145,11 +146,6 @@ inline RealType pdf(const normal_distribution& dist, const Rea { return 0; // pdf + and - infinity is zero. } - // Below produces MSVC 4127 warnings, so the above used instead. - //if(std::numeric_limits::has_infinity && abs(x) == std::numeric_limits::infinity()) - //{ // pdf + and - infinity is zero. - // return 0; - //} if(false == detail::check_x(function, x, &result, Policy())) { return result; @@ -165,6 +161,42 @@ inline RealType pdf(const normal_distribution& dist, const Rea return result; } // pdf +template +inline RealType logpdf(const normal_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + const RealType sd = dist.standard_deviation(); + const RealType mean = dist.mean(); + + static const char* function = "boost::math::logpdf(const normal_distribution<%1%>&, %1%)"; + + RealType result = -std::numeric_limits::infinity(); + if(false == detail::check_scale(function, sd, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, mean, &result, Policy())) + { + return result; + } + if((boost::math::isinf)(x)) + { + return result; // pdf + and - infinity is zero so logpdf is -inf + } + if(false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + const RealType pi = boost::math::constants::pi(); + const RealType half = boost::math::constants::half(); + + result = -log(sd) - half*log(2*pi) - (x-mean)*(x-mean)/(2*sd*sd); + + return result; +} + template inline RealType cdf(const normal_distribution& dist, const RealType& x) { @@ -187,15 +219,6 @@ inline RealType cdf(const normal_distribution& dist, const Rea if(x < 0) return 0; // -infinity return 1; // + infinity } - // These produce MSVC 4127 warnings, so the above used instead. - //if(std::numeric_limits::has_infinity && x == std::numeric_limits::infinity()) - //{ // cdf +infinity is unity. - // return 1; - //} - //if(std::numeric_limits::has_infinity && x == -std::numeric_limits::infinity()) - //{ // cdf -infinity is zero. - // return 0; - //} if(false == detail::check_x(function, x, &result, Policy())) { return result; @@ -249,15 +272,6 @@ inline RealType cdf(const complemented2_type::has_infinity && x == std::numeric_limits::infinity()) - //{ // cdf complement +infinity is zero. - // return 0; - //} - //if(std::numeric_limits::has_infinity && x == -std::numeric_limits::infinity()) - //{ // cdf complement -infinity is unity. - // return 1; - //} if(false == detail::check_x(function, x, &result, Policy())) return result; diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp index 533c31a80d..5507360e82 100644 --- a/include/boost/math/distributions/poisson.hpp +++ b/include/boost/math/distributions/poisson.hpp @@ -47,6 +47,7 @@ #include #include +#include namespace boost { @@ -144,10 +145,10 @@ namespace boost class poisson_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - poisson_distribution(RealType l_mean = 1) : m_l(l_mean) // mean (lambda). + explicit poisson_distribution(RealType l_mean = 1) : m_l(l_mean) // mean (lambda). { // Expected mean number of events that occur during the given interval. RealType r; poisson_detail::check_dist( @@ -165,7 +166,7 @@ namespace boost RealType m_l; // mean number of occurrences. }; // template class poisson_distribution - typedef poisson_distribution poisson; // Reserved name of type double. + using poisson = poisson_distribution; // Reserved name of type double. #ifdef __cpp_deduction_guides template @@ -175,14 +176,14 @@ namespace boost // Non-member functions to give properties of the distribution. template - inline const std::pair range(const poisson_distribution& /* dist */) + inline std::pair range(const poisson_distribution& /* dist */) { // Range of permissible values for random variable k. using boost::math::tools::max_value; return std::pair(static_cast(0), max_value()); // Max integer? } template - inline const std::pair support(const poisson_distribution& /* dist */) + inline std::pair support(const poisson_distribution& /* dist */) { // Range of supported values for random variable k. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -202,15 +203,7 @@ namespace boost return floor(dist.mean()); } - //template - //inline RealType median(const poisson_distribution& dist) - //{ // median = approximately lambda + 1/3 - 0.2/lambda - // RealType l = dist.mean(); - // return dist.mean() + static_cast(0.3333333333333333333333333333333333333333333333) - // - static_cast(0.2) / l; - //} // BUT this formula appears to be out-by-one compared to quantile(half) - // Query posted on Wikipedia. - // Now implemented via quantile(half) in derived accessors. + // Median now implemented via quantile(half) in derived accessors. template inline RealType variance(const poisson_distribution& dist) @@ -218,7 +211,6 @@ namespace boost return dist.mean(); } - // RealType standard_deviation(const poisson_distribution& dist) // standard_deviation provided by derived accessors. template @@ -281,6 +273,42 @@ namespace boost return boost::math::gamma_p_derivative(k+1, mean, Policy()); } // pdf + template + RealType logpdf(const poisson_distribution& dist, const RealType& k) + { + BOOST_FPU_EXCEPTION_GUARD + + BOOST_MATH_STD_USING // for ADL of std functions. + using boost::math::lgamma; + + RealType mean = dist.mean(); + // Error check: + RealType result = -std::numeric_limits::infinity(); + if(false == poisson_detail::check_dist_and_k( + "boost::math::pdf(const poisson_distribution<%1%>&, %1%)", + mean, + k, + &result, Policy())) + { + return result; + } + + // Special case of mean zero, regardless of the number of events k. + if (mean == 0) + { // Probability for any k is zero. + return std::numeric_limits::quiet_NaN(); + } + + // Special case where k and lambda are both positive + if(k > 0 && mean > 0) + { + return -lgamma(k+1) + k*log(mean) - mean; + } + + result = log(pdf(dist, k)); + return result; + } + template RealType cdf(const poisson_distribution& dist, const RealType& k) { // Cumulative Distribution Function Poisson. @@ -320,8 +348,8 @@ namespace boost return 0; } if (k == 0) - { // return pdf(dist, static_cast(0)); - // but mean (and k) have already been checked, + { + // mean (and k) have already been checked, // so this avoids unnecessary repeated checks. return exp(-mean); } @@ -414,9 +442,10 @@ namespace boost { return policies::raise_overflow_error(function, 0, Policy()); } - typedef typename Policy::discrete_quantile_type discrete_type; + using discrete_type = typename Policy::discrete_quantile_type; std::uintmax_t max_iter = policies::get_max_root_iterations(); - RealType guess, factor = 8; + RealType guess; + RealType factor = 8; RealType z = dist.mean(); if(z < 1) guess = z; @@ -484,9 +513,10 @@ namespace boost { return 0; // Exact result regardless of discrete-quantile Policy } - typedef typename Policy::discrete_quantile_type discrete_type; + using discrete_type = typename Policy::discrete_quantile_type; std::uintmax_t max_iter = policies::get_max_root_iterations(); - RealType guess, factor = 8; + RealType guess; + RealType factor = 8; RealType z = dist.mean(); if(z < 1) guess = z; diff --git a/include/boost/math/distributions/rayleigh.hpp b/include/boost/math/distributions/rayleigh.hpp index cbe934471f..cbbf39876d 100644 --- a/include/boost/math/distributions/rayleigh.hpp +++ b/include/boost/math/distributions/rayleigh.hpp @@ -56,10 +56,10 @@ template > class rayleigh_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - rayleigh_distribution(RealType l_sigma = 1) + explicit rayleigh_distribution(RealType l_sigma = 1) : m_sigma(l_sigma) { RealType err; @@ -75,7 +75,7 @@ class rayleigh_distribution RealType m_sigma; }; // class rayleigh_distribution -typedef rayleigh_distribution rayleigh; +using rayleigh = rayleigh_distribution; #ifdef __cpp_deduction_guides template @@ -83,14 +83,14 @@ rayleigh_distribution(RealType)->rayleigh_distribution -inline const std::pair range(const rayleigh_distribution& /*dist*/) +inline std::pair range(const rayleigh_distribution& /*dist*/) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(0), std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : max_value()); } template -inline const std::pair support(const rayleigh_distribution& /*dist*/) +inline std::pair support(const rayleigh_distribution& /*dist*/) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -122,6 +122,32 @@ inline RealType pdf(const rayleigh_distribution& dist, const R return result; } // pdf +template +inline RealType logpdf(const rayleigh_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std function exp. + + const RealType sigma = dist.sigma(); + RealType result = -std::numeric_limits::infinity(); + static const char* function = "boost::math::logpdf(const rayleigh_distribution<%1%>&, %1%)"; + + if(false == detail::verify_sigma(function, sigma, &result, Policy())) + { + return result; + } + if(false == detail::verify_rayleigh_x(function, x, &result, Policy())) + { + return result; + } + if((boost::math::isinf)(x)) + { + return result; + } + + result = -(x*x)/(2*sigma*sigma) - 2*log(sigma) + log(x); + return result; +} // logpdf + template inline RealType cdf(const rayleigh_distribution& dist, const RealType& x) { @@ -265,31 +291,26 @@ inline RealType median(const rayleigh_distribution& dist) template inline RealType skewness(const rayleigh_distribution& /*dist*/) { - // using namespace boost::math::constants; return static_cast(0.63111065781893713819189935154422777984404221106391L); // Computed using NTL at 150 bit, about 50 decimal digits. - // return 2 * root_pi() * pi_minus_three() / pow23_four_minus_pi(); + // 2 * sqrt(pi) * (pi-3) / pow(4, 2/3) - pi } template inline RealType kurtosis(const rayleigh_distribution& /*dist*/) { - // using namespace boost::math::constants; return static_cast(3.2450893006876380628486604106197544154170667057995L); // Computed using NTL at 150 bit, about 50 decimal digits. - // return 3 - (6 * pi() * pi() - 24 * pi() + 16) / - // (four_minus_pi() * four_minus_pi()); + // 3 - (6*pi*pi - 24*pi + 16) / pow(4-pi, 2) } template inline RealType kurtosis_excess(const rayleigh_distribution& /*dist*/) { - //using namespace boost::math::constants; - // Computed using NTL at 150 bit, about 50 decimal digits. return static_cast(0.2450893006876380628486604106197544154170667057995L); - // return -(6 * pi() * pi() - 24 * pi() + 16) / - // (four_minus_pi() * four_minus_pi()); -} // kurtosis + // Computed using NTL at 150 bit, about 50 decimal digits. + // -(6*pi*pi - 24*pi + 16) / pow(4-pi,2) +} // kurtosis_excess template inline RealType entropy(const rayleigh_distribution& dist) diff --git a/include/boost/math/distributions/weibull.hpp b/include/boost/math/distributions/weibull.hpp index 76246e4ff2..724cce6ed8 100644 --- a/include/boost/math/distributions/weibull.hpp +++ b/include/boost/math/distributions/weibull.hpp @@ -70,10 +70,10 @@ template > class weibull_distribution { public: - typedef RealType value_type; - typedef Policy policy_type; + using value_type = RealType; + using policy_type = Policy; - weibull_distribution(RealType l_shape, RealType l_scale = 1) + explicit weibull_distribution(RealType l_shape, RealType l_scale = 1) : m_shape(l_shape), m_scale(l_scale) { RealType result; @@ -97,7 +97,7 @@ class weibull_distribution RealType m_scale; // distribution scale }; -typedef weibull_distribution weibull; +using weibull = weibull_distribution; #ifdef __cpp_deduction_guides template @@ -107,14 +107,14 @@ weibull_distribution(RealType,RealType)->weibull_distribution -inline const std::pair range(const weibull_distribution& /*dist*/) +inline std::pair range(const weibull_distribution& /*dist*/) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(0), max_value()); } template -inline const std::pair support(const weibull_distribution& /*dist*/) +inline std::pair support(const weibull_distribution& /*dist*/) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; @@ -157,6 +157,40 @@ inline RealType pdf(const weibull_distribution& dist, const Re return result; } +template +inline RealType logpdf(const weibull_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + static const char* function = "boost::math::logpdf(const weibull_distribution<%1%>, %1%)"; + + RealType shape = dist.shape(); + RealType scale = dist.scale(); + + RealType result = 0; + if(false == detail::check_weibull(function, scale, shape, &result, Policy())) + return result; + if(false == detail::check_weibull_x(function, x, &result, Policy())) + return result; + + if(x == 0) + { + if(shape == 1) + { + return log(1 / scale); + } + if(shape > 1) + { + return 0; + } + return policies::raise_overflow_error(function, 0, Policy()); + } + + result = log(shape) - shape * log(scale) + log(x) * (shape - 1) - pow(x / scale, shape); + + return result; +} + template inline RealType cdf(const weibull_distribution& dist, const RealType& x) { @@ -342,12 +376,11 @@ inline RealType skewness(const weibull_distribution& dist) { return result; } - RealType g1, g2, g3, d; - g1 = boost::math::tgamma(1 + 1 / shape, Policy()); - g2 = boost::math::tgamma(1 + 2 / shape, Policy()); - g3 = boost::math::tgamma(1 + 3 / shape, Policy()); - d = pow(g2 - g1 * g1, RealType(1.5)); + RealType g1 = boost::math::tgamma(1 + 1 / shape, Policy()); + RealType g2 = boost::math::tgamma(1 + 2 / shape, Policy()); + RealType g3 = boost::math::tgamma(1 + 3 / shape, Policy()); + RealType d = pow(g2 - g1 * g1, RealType(1.5)); result = (2 * g1 * g1 * g1 - 3 * g1 * g2 + g3) / d; return result; @@ -367,15 +400,13 @@ inline RealType kurtosis_excess(const weibull_distribution& di if(false == detail::check_weibull(function, scale, shape, &result, Policy())) return result; - RealType g1, g2, g3, g4, d, g1_2, g1_4; - - g1 = boost::math::tgamma(1 + 1 / shape, Policy()); - g2 = boost::math::tgamma(1 + 2 / shape, Policy()); - g3 = boost::math::tgamma(1 + 3 / shape, Policy()); - g4 = boost::math::tgamma(1 + 4 / shape, Policy()); - g1_2 = g1 * g1; - g1_4 = g1_2 * g1_2; - d = g2 - g1_2; + RealType g1 = boost::math::tgamma(1 + 1 / shape, Policy()); + RealType g2 = boost::math::tgamma(1 + 2 / shape, Policy()); + RealType g3 = boost::math::tgamma(1 + 3 / shape, Policy()); + RealType g4 = boost::math::tgamma(1 + 4 / shape, Policy()); + RealType g1_2 = g1 * g1; + RealType g1_4 = g1_2 * g1_2; + RealType d = g2 - g1_2; d *= d; result = -6 * g1_4 + 12 * g1_2 * g2 - 3 * g2 * g2 - 4 * g1 * g3 + g4; diff --git a/include/boost/math/policies/error_handling.hpp b/include/boost/math/policies/error_handling.hpp index 4ee6877e27..3af1f976b5 100644 --- a/include/boost/math/policies/error_handling.hpp +++ b/include/boost/math/policies/error_handling.hpp @@ -8,19 +8,23 @@ #ifndef BOOST_MATH_POLICY_ERROR_HANDLING_HPP #define BOOST_MATH_POLICY_ERROR_HANDLING_HPP +#include #include #include #include +#ifndef BOOST_NO_RTTI #include +#endif #include #include #include #include -#include -#include #include #include +#ifndef BOOST_NO_EXCEPTIONS +#include #include +#endif #ifdef _MSC_VER # pragma warning(push) // Quiet warnings in boost/format.hpp @@ -36,6 +40,8 @@ namespace boost{ namespace math{ +#ifndef BOOST_NO_EXCEPTIONS + class evaluation_error : public std::runtime_error { public: @@ -48,6 +54,8 @@ class rounding_error : public std::runtime_error rounding_error(const std::string& s) : std::runtime_error(s){} }; +#endif + namespace policies{ // // Forward declarations of user error handlers, @@ -120,6 +128,7 @@ inline const char* name_of() } #endif +#ifndef BOOST_NO_EXCEPTIONS template void raise_error(const char* pfunction, const char* message) { @@ -169,6 +178,7 @@ void raise_error(const char* pfunction, const char* pmessage, const T& val) E e(msg); BOOST_MATH_THROW_EXCEPTION(e) } +#endif template inline T raise_domain_error( @@ -177,9 +187,13 @@ inline T raise_domain_error( const T& val, const ::boost::math::policies::domain_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else raise_error(function, message, val); // we never get here: return std::numeric_limits::quiet_NaN(); +#endif } template @@ -224,7 +238,11 @@ inline T raise_pole_error( const T& val, const ::boost::math::policies::pole_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else return boost::math::policies::detail::raise_domain_error(function, message, val, ::boost::math::policies::domain_error< ::boost::math::policies::throw_on_error>()); +#endif } template @@ -257,16 +275,19 @@ inline T raise_pole_error( return user_pole_error(function, message, val); } - template inline T raise_overflow_error( const char* function, const char* message, const ::boost::math::policies::overflow_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else raise_error(function, message ? message : "numeric overflow"); // We should never get here: return std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(); +#endif } template @@ -276,9 +297,13 @@ inline T raise_overflow_error( const T& val, const ::boost::math::policies::overflow_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else raise_error(function, message ? message : "numeric overflow", val); // We should never get here: return std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(); +#endif } template @@ -358,9 +383,13 @@ inline T raise_underflow_error( const char* message, const ::boost::math::policies::underflow_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else raise_error(function, message ? message : "numeric underflow"); // We should never get here: return 0; +#endif } template @@ -402,9 +431,13 @@ inline T raise_denorm_error( const T& /* val */, const ::boost::math::policies::denorm_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else raise_error(function, message ? message : "denormalised result"); // we never get here: return T(0); +#endif } template @@ -449,9 +482,13 @@ inline T raise_evaluation_error( const T& val, const ::boost::math::policies::evaluation_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else raise_error(function, message, val); // we never get here: return T(0); +#endif } template @@ -497,9 +534,13 @@ inline TargetType raise_rounding_error( const TargetType&, const ::boost::math::policies::rounding_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else raise_error(function, message, val); // we never get here: return TargetType(0); +#endif } template @@ -564,9 +605,13 @@ inline T raise_indeterminate_result_error( const R& , const ::boost::math::policies::indeterminate_result_error< ::boost::math::policies::throw_on_error>&) { +#ifdef BOOST_NO_EXCEPTIONS + static_assert(sizeof(T) == 0, "Error handler called with throw_on_error and BOOST_NO_EXCEPTIONS set."); +#else raise_error(function, message, val); // we never get here: return std::numeric_limits::quiet_NaN(); +#endif } template diff --git a/include/boost/math/special_functions/chebyshev.hpp b/include/boost/math/special_functions/chebyshev.hpp index 66356ee012..7aff765e6f 100644 --- a/include/boost/math/special_functions/chebyshev.hpp +++ b/include/boost/math/special_functions/chebyshev.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #if (__cplusplus > 201103) || (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610)) # define BOOST_MATH_CHEB_USE_STD_ACOSH @@ -271,7 +272,7 @@ inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, co { if (x < a || x > b) { - throw std::domain_error("x in [a, b] is required."); + BOOST_MATH_THROW_EXCEPTION(std::domain_error("x in [a, b] is required.")); } if (length < 2) { diff --git a/include/boost/math/special_functions/detail/bernoulli_details.hpp b/include/boost/math/special_functions/detail/bernoulli_details.hpp index 91f1747c01..58015fec14 100644 --- a/include/boost/math/special_functions/detail/bernoulli_details.hpp +++ b/include/boost/math/special_functions/detail/bernoulli_details.hpp @@ -184,17 +184,22 @@ struct fixed_vector : private std::allocator const T& operator[](unsigned n)const { BOOST_MATH_ASSERT(n < m_used); return m_data[n]; } unsigned size()const { return m_used; } unsigned size() { return m_used; } - void resize(unsigned n, const T& val) + bool resize(unsigned n, const T& val) { if(n > m_capacity) { +#ifndef BOOST_NO_EXCEPTIONS BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers.")); +#else + return false; +#endif } for(unsigned i = m_used; i < n; ++i) new (m_data + i) T(val); m_used = n; + return true; } - void resize(unsigned n) { resize(n, T()); } + bool resize(unsigned n) { return resize(n, T()); } T* begin() { return m_data; } T* end() { return m_data + m_used; } T* begin()const { return m_data; } @@ -217,10 +222,14 @@ class bernoulli_numbers_cache typedef fixed_vector container_type; - void tangent(std::size_t m) + bool tangent(std::size_t m) { static const std::size_t min_overflow_index = b2n_overflow_limit() - 1; - tn.resize(static_cast(m), T(0U)); + + if (!tn.resize(static_cast(m), T(0U))) + { + return false; + } BOOST_MATH_INSTRUMENT_VARIABLE(min_overflow_index); @@ -268,17 +277,20 @@ class bernoulli_numbers_cache BOOST_MATH_INSTRUMENT_VARIABLE(i); BOOST_MATH_INSTRUMENT_VARIABLE(tn[static_cast(i)]); } + return true; } - void tangent_numbers_series(const std::size_t m) + bool tangent_numbers_series(const std::size_t m) { BOOST_MATH_STD_USING static const std::size_t min_overflow_index = b2n_overflow_limit() - 1; typename container_type::size_type old_size = bn.size(); - tangent(m); - bn.resize(static_cast(m)); + if (!tangent(m)) + return false; + if (!bn.resize(static_cast(m))) + return false; if(!old_size) { @@ -321,6 +333,7 @@ class bernoulli_numbers_cache bn[static_cast(i)] = ((!b_neg) ? b : T(-b)); } + return true; } template @@ -379,7 +392,10 @@ class bernoulli_numbers_cache if(start + n >= bn.size()) { std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity())); - tangent_numbers_series(new_size); + if (!tangent_numbers_series(new_size)) + { + return std::fill_n(out, n, policies::raise_evaluation_error("boost::math::bernoulli_b2n<%1%>(std::size_t)", "Unable to allocate Bernoulli numbers cache for %1% values", T(start + n), pol)); + } } for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n::value + 1), start); i < start + n; ++i) @@ -413,7 +429,8 @@ class bernoulli_numbers_cache if(start + n >= bn.size()) { std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity())); - tangent_numbers_series(new_size); + if (!tangent_numbers_series(new_size)) + return std::fill_n(out, n, policies::raise_evaluation_error("boost::math::bernoulli_b2n<%1%>(std::size_t)", "Unable to allocate Bernoulli numbers cache for %1% values", T(new_size), pol)); } m_counter.store(static_cast(bn.size()), std::memory_order_release); } @@ -483,7 +500,8 @@ class bernoulli_numbers_cache if(start + n >= bn.size()) { std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity())); - tangent_numbers_series(new_size); + if (!tangent_numbers_series(new_size)) + return std::fill_n(out, n, policies::raise_evaluation_error("boost::math::bernoulli_b2n<%1%>(std::size_t)", "Unable to allocate Bernoulli numbers cache for %1% values", T(start + n), pol)); } for(std::size_t i = start; i < start + n; ++i) @@ -527,7 +545,8 @@ class bernoulli_numbers_cache if(start + n >= bn.size()) { std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity())); - tangent_numbers_series(new_size); + if (!tangent_numbers_series(new_size)) + return std::fill_n(out, n, policies::raise_evaluation_error("boost::math::bernoulli_b2n<%1%>(std::size_t)", "Unable to allocate Bernoulli numbers cache for %1% values", T(start + n), pol)); } m_counter.store(static_cast(bn.size()), std::memory_order_release); } diff --git a/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp b/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp index 6f1461ffc9..07cb236db1 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp @@ -291,6 +291,14 @@ ak += 2; integer_part -= 2; } + if (ak - 1 == b) + { + // When ak - 1 == b are recursion coefficients dissappear to zero and + // we end up with a NaN result. Reduce the recursion steps by 1 to + // avoid this. We rely on |b| small and therefore no infinite recursion. + ak -= 1; + integer_part += 1; + } if (-integer_part > static_cast(policies::get_max_series_iterations())) return policies::raise_evaluation_error(function, "1F1 arguments sit in a range with a so negative that we have no evaluation method, got a = %1%", std::numeric_limits::quiet_NaN(), pol); diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp index 3745f7e24d..8ed8c8086f 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp @@ -132,6 +132,7 @@ Real scaling_factor = exp(Real(log_scaling_factor)); Real term_m1; long long local_scaling = 0; + bool have_no_correct_bits = false; if ((aj.size() == 1) && (bj.size() == 0)) { @@ -238,10 +239,20 @@ break; if (abs_result * tol > abs(result)) { - // We have no correct bits in the result... just give up! - result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the reuslt are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol); - return std::make_pair(result, result); + // Check if result is so small compared to abs_resuslt that there are no longer any + // correct bits... we require two consecutive passes here before aborting to + // avoid false positives when result transiently drops to near zero then rebounds. + if (have_no_correct_bits) + { + // We have no correct bits in the result... just give up! + result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the reuslt are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol); + return std::make_pair(result, result); + } + else + have_no_correct_bits = true; } + else + have_no_correct_bits = false; term0 = term; } //std::cout << "result = " << result << std::endl; diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp index e3dc0ba4fa..195f6abe18 100644 --- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp +++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp @@ -688,13 +688,17 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py) T bet = 0; T xg; bool overflow = false; +#ifndef BOOST_NO_EXCEPTIONS try { +#endif bet = boost::math::beta(a, b, pol); +#ifndef BOOST_NO_EXCEPTIONS } catch (const std::runtime_error&) { overflow = true; } +#endif if (overflow || !(boost::math::isfinite)(bet)) { xg = exp((boost::math::lgamma(a + 1, pol) + boost::math::lgamma(b, pol) - boost::math::lgamma(a + b, pol) + log(p)) / a); diff --git a/include/boost/math/special_functions/detail/unchecked_factorial.hpp b/include/boost/math/special_functions/detail/unchecked_factorial.hpp index 3bdc45fe5f..5fa9f0b6b4 100644 --- a/include/boost/math/special_functions/detail/unchecked_factorial.hpp +++ b/include/boost/math/special_functions/detail/unchecked_factorial.hpp @@ -1009,7 +1009,7 @@ struct max_factorial #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION template -const unsigned max_factorial::value; +constexpr unsigned max_factorial::value; #endif } // namespace math diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index b94dcbbc32..13cd25e087 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1465,14 +1465,18 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, result = pow(x, a) / (a); else { +#ifndef BOOST_NO_EXCEPTIONS try { +#endif result = pow(x, a) / boost::math::tgamma(a + 1, pol); +#ifndef BOOST_NO_EXCEPTIONS } catch (const std::overflow_error&) { result = 0; } +#endif } result *= 1 - a * x / (a + 1); if (p_derivative) diff --git a/include/boost/math/special_functions/owens_t.hpp b/include/boost/math/special_functions/owens_t.hpp index f96d6648e5..152b5cbdda 100644 --- a/include/boost/math/special_functions/owens_t.hpp +++ b/include/boost/math/special_functions/owens_t.hpp @@ -829,7 +829,7 @@ namespace boost val = owens_t_T6(h,a, pol); break; default: - BOOST_MATH_THROW_EXCEPTION(std::logic_error("selection routine in Owen's T function failed")); + val = policies::raise_evaluation_error("boost::math::owens_t", "selection routine in Owen's T function failed with h = %1%", h, pol); } return val; } diff --git a/include/boost/math/statistics/bivariate_statistics.hpp b/include/boost/math/statistics/bivariate_statistics.hpp index 1854898139..d6f91faa93 100644 --- a/include/boost/math/statistics/bivariate_statistics.hpp +++ b/include/boost/math/statistics/bivariate_statistics.hpp @@ -18,13 +18,10 @@ #include #include -// Support compilers with P0024R2 implemented without linking TBB -// https://en.cppreference.com/w/cpp/compiler_support -#if !defined(BOOST_NO_CXX17_HDR_EXECUTION) && defined(BOOST_HAS_THREADS) +#ifdef BOOST_MATH_EXEC_COMPATIBLE #include #include #include -#define EXEC_COMPATIBLE #endif namespace boost{ namespace math{ namespace statistics { namespace detail { @@ -60,7 +57,7 @@ ReturnType means_and_covariance_seq_impl(ForwardIterator u_begin, ForwardIterato return std::make_tuple(mu_u, mu_v, cov/i, Real(i)); } -#ifdef EXEC_COMPATIBLE +#ifdef BOOST_MATH_EXEC_COMPATIBLE // Numerically stable parallel computation of (co-)variance // https://dl.acm.org/doi/10.1145/3221269.3223036 @@ -154,7 +151,7 @@ ReturnType means_and_covariance_parallel_impl(ForwardIterator u_begin, ForwardIt return std::make_tuple(mu_u_a, mu_v_a, cov_a, n_a); } -#endif // EXEC_COMPATIBLE +#endif // BOOST_MATH_EXEC_COMPATIBLE template ReturnType correlation_coefficient_seq_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end) @@ -204,7 +201,7 @@ ReturnType correlation_coefficient_seq_impl(ForwardIterator u_begin, ForwardIter return std::make_tuple(mu_u, Qu, mu_v, Qv, cov, rho, Real(i)); } -#ifdef EXEC_COMPATIBLE +#ifdef BOOST_MATH_EXEC_COMPATIBLE // Numerically stable parallel computation of (co-)variance: // https://dl.acm.org/doi/10.1145/3221269.3223036 @@ -324,11 +321,11 @@ ReturnType correlation_coefficient_parallel_impl(ForwardIterator u_begin, Forwar return std::make_tuple(mu_u_a, Qu_a, mu_v_a, Qv_a, cov_a, rho, n_a); } -#endif // EXEC_COMPATIBLE +#endif // BOOST_MATH_EXEC_COMPATIBLE } // namespace detail -#ifdef EXEC_COMPATIBLE +#ifdef BOOST_MATH_EXEC_COMPATIBLE template inline auto means_and_covariance(ExecutionPolicy&& exec, Container const & u, Container const & v) diff --git a/include/boost/math/statistics/chatterjee_correlation.hpp b/include/boost/math/statistics/chatterjee_correlation.hpp new file mode 100644 index 0000000000..ad0b33a429 --- /dev/null +++ b/include/boost/math/statistics/chatterjee_correlation.hpp @@ -0,0 +1,159 @@ +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_STATISTICS_CHATTERJEE_CORRELATION_HPP +#define BOOST_MATH_STATISTICS_CHATTERJEE_CORRELATION_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_MATH_EXEC_COMPATIBLE +#include +#include +#include +#endif + +namespace boost { namespace math { namespace statistics { + +namespace detail { + +template +std::size_t chatterjee_transform(BDIter begin, BDIter end) +{ + std::size_t sum = 0; + + while(++begin != end) + { + if(*begin > *std::prev(begin)) + { + sum += *begin - *std::prev(begin); + } + else + { + sum += *std::prev(begin) - *begin; + } + } + + return sum; +} + +template +ReturnType chatterjee_correlation_seq_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end) +{ + using std::abs; + + BOOST_MATH_ASSERT_MSG(std::is_sorted(u_begin, u_end), "The x values must be sorted in order to use this functionality"); + + const std::vector rank_vector = rank(v_begin, v_end); + + std::size_t sum = chatterjee_transform(rank_vector.begin(), rank_vector.end()); + + ReturnType result = static_cast(1) - (static_cast(3 * sum) / static_cast(rank_vector.size() * rank_vector.size() - 1)); + + // If the result is 1 then Y is constant and all the elements must be ties + if (abs(result - static_cast(1)) < std::numeric_limits::epsilon()) + { + return std::numeric_limits::quiet_NaN(); + } + + return result; +} + +} // Namespace detail + +template ::value, double, Real>::type> +inline ReturnType chatterjee_correlation(const Container& u, const Container& v) +{ + return detail::chatterjee_correlation_seq_impl(std::begin(u), std::end(u), std::begin(v), std::end(v)); +} + +}}} // Namespace boost::math::statistics + +#ifdef BOOST_MATH_EXEC_COMPATIBLE + +namespace boost::math::statistics { + +namespace detail { + +template +ReturnType chatterjee_correlation_par_impl(ExecutionPolicy&& exec, ForwardIterator u_begin, ForwardIterator u_end, + ForwardIterator v_begin, ForwardIterator v_end) +{ + using std::abs; + BOOST_MATH_ASSERT_MSG(std::is_sorted(std::forward(exec), u_begin, u_end), "The x values must be sorted in order to use this functionality"); + + auto rank_vector = rank(std::forward(exec), v_begin, v_end); + + const auto num_threads = std::thread::hardware_concurrency() == 0 ? 2u : std::thread::hardware_concurrency(); + std::vector> future_manager {}; + const auto elements_per_thread = std::ceil(static_cast(rank_vector.size()) / num_threads); + + auto it = rank_vector.begin(); + auto end = rank_vector.end(); + for(std::size_t i {}; i < num_threads - 1; ++i) + { + future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [it, elements_per_thread]() -> std::size_t + { + return chatterjee_transform(it, std::next(it, elements_per_thread)); + })); + it = std::next(it, elements_per_thread - 1); + } + + future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [it, end]() -> std::size_t + { + return chatterjee_transform(it, end); + })); + + std::size_t sum {}; + for(std::size_t i {}; i < future_manager.size(); ++i) + { + sum += future_manager[i].get(); + } + + ReturnType result = static_cast(1) - (static_cast(3 * sum) / static_cast(rank_vector.size() * rank_vector.size() - 1)); + + // If the result is 1 then Y is constant and all the elements must be ties + if (abs(result - static_cast(1)) < std::numeric_limits::epsilon()) + { + return std::numeric_limits::quiet_NaN(); + } + + return result; +} + +} // Namespace detail + +template , double, Real>> +inline ReturnType chatterjee_correlation(ExecutionPolicy&& exec, const Container& u, const Container& v) +{ + if constexpr (std::is_same_v, decltype(std::execution::seq)>) + { + return detail::chatterjee_correlation_seq_impl(std::cbegin(u), std::cend(u), + std::cbegin(v), std::cend(v)); + } + else + { + return detail::chatterjee_correlation_par_impl(std::forward(exec), + std::cbegin(u), std::cend(u), + std::cbegin(v), std::cend(v)); + } +} + +} // Namespace boost::math::statistics + +#endif + +#endif // BOOST_MATH_STATISTICS_CHATTERJEE_CORRELATION_HPP diff --git a/include/boost/math/statistics/detail/rank.hpp b/include/boost/math/statistics/detail/rank.hpp new file mode 100644 index 0000000000..4e5211607e --- /dev/null +++ b/include/boost/math/statistics/detail/rank.hpp @@ -0,0 +1,140 @@ +// (C) Copyright Matt Borland 2022 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_STATISTICS_DETAIL_RANK_HPP +#define BOOST_MATH_STATISTICS_DETAIL_RANK_HPP + +#include +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_MATH_EXEC_COMPATIBLE +#include +#endif + +namespace boost { namespace math { namespace statistics { namespace detail { + +struct pair_equal +{ + template + bool operator()(const std::pair& a, const std::pair& b) const + { + return a.first == b.first; + } +}; + +}}}} // Namespaces + +#ifndef BOOST_MATH_EXEC_COMPATIBLE + +namespace boost { namespace math { namespace statistics { namespace detail { + +template ::value_type> +auto rank(ForwardIterator first, ForwardIterator last) -> std::vector +{ + std::size_t elements = std::distance(first, last); + + std::vector> rank_vector(elements); + std::size_t i = 0; + while (first != last) + { + rank_vector[i] = std::make_pair(*first, i); + ++i; + ++first; + } + + std::sort(rank_vector.begin(), rank_vector.end()); + + // Remove duplicates + rank_vector.erase(std::unique(rank_vector.begin(), rank_vector.end(), pair_equal()), rank_vector.end()); + elements = rank_vector.size(); + + std::pair rank; + std::vector result(elements); + for (i = 0; i < elements; ++i) + { + if (rank_vector[i].first != rank.first) + { + rank = std::make_pair(rank_vector[i].first, i); + } + result[rank_vector[i].second] = rank.second; + } + + return result; +} + +template +inline auto rank(const Container& c) -> std::vector +{ + return rank(std::begin(c), std::end(c)); +} + +}}}} // Namespaces + +#else + +namespace boost::math::statistics::detail { + +template ::value_type> +auto rank(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) +{ + std::size_t elements = std::distance(first, last); + + std::vector> rank_vector(elements); + std::size_t i = 0; + while (first != last) + { + rank_vector[i] = std::make_pair(*first, i); + ++i; + ++first; + } + + std::sort(exec, rank_vector.begin(), rank_vector.end()); + + // Remove duplicates + rank_vector.erase(std::unique(exec, rank_vector.begin(), rank_vector.end(), pair_equal()), rank_vector.end()); + elements = rank_vector.size(); + + std::pair rank; + std::vector result(elements); + for (i = 0; i < elements; ++i) + { + if (rank_vector[i].first != rank.first) + { + rank = std::make_pair(rank_vector[i].first, i); + } + result[rank_vector[i].second] = rank.second; + } + + return result; +} + +template +inline auto rank(ExecutionPolicy&& exec, const Container& c) +{ + return rank(exec, std::cbegin(c), std::cend(c)); +} + +template ::value_type> +inline auto rank(ForwardIterator first, ForwardIterator last) +{ + return rank(std::execution::seq, first, last); +} + +template +inline auto rank(const Container& c) +{ + return rank(std::execution::seq, std::cbegin(c), std::cend(c)); +} + +} // Namespaces + +#endif // BOOST_MATH_EXEC_COMPATIBLE + +#endif // BOOST_MATH_STATISTICS_DETAIL_RANK_HPP diff --git a/include/boost/math/statistics/univariate_statistics.hpp b/include/boost/math/statistics/univariate_statistics.hpp index 8a366384d7..7da6e5ba4f 100644 --- a/include/boost/math/statistics/univariate_statistics.hpp +++ b/include/boost/math/statistics/univariate_statistics.hpp @@ -20,9 +20,7 @@ #include #include -// Support compilers with P0024R2 implemented without linking TBB -// https://en.cppreference.com/w/cpp/compiler_support -#if !defined(BOOST_NO_CXX17_HDR_EXECUTION) && defined(BOOST_HAS_THREADS) +#ifdef BOOST_MATH_EXEC_COMPATIBLE #include namespace boost::math::statistics { @@ -699,7 +697,7 @@ inline auto mode(Container & v) } // Namespace boost::math::statistics -#else // Backwards compatible bindings for C++11 +#else // Backwards compatible bindings for C++11 or execution is not implemented namespace boost { namespace math { namespace statistics { diff --git a/include/boost/math/tools/color_maps.hpp b/include/boost/math/tools/color_maps.hpp index 073f65d1ef..e3cc8b9b09 100644 --- a/include/boost/math/tools/color_maps.hpp +++ b/include/boost/math/tools/color_maps.hpp @@ -10,6 +10,7 @@ #include // for table data #include // for std::floor #include // fixed width integer types +#include #if __has_include("lodepng.h") @@ -1863,12 +1864,21 @@ static constexpr std::array, 256> viridis_data_ = { template inline std::array color_map_(Real scalar, std::array, 256> const &table) { - static_assert(std::is_floating_point::value, + static_assert(std::is_floating_point_v, "Color tables are only implemented in floating point " "arithmetic. If you require bytes please submit an issue or " "pull request"); - scalar = std::clamp(scalar, static_cast(0), static_cast(1)); + using boost::math::isnan; + + if ((isnan)(scalar)) + { + scalar = static_cast(0); + } + else + { + scalar = std::clamp(scalar, static_cast(0), static_cast(1)); + } if (scalar == static_cast(1)) { return table.back(); @@ -1877,7 +1887,7 @@ color_map_(Real scalar, std::array, 256> const &table) { Real s = (table.size() - 1) * scalar; Real ii = std::floor(s); Real t = s - ii; - std::size_t i = static_cast(ii); + auto i = static_cast(ii); auto const &rgb0 = table[i]; auto const &rgb1 = table[i + 1]; return {(1 - t) * rgb0[0] + t * rgb1[0], (1 - t) * rgb0[1] + t * rgb1[1], @@ -1917,7 +1927,7 @@ std::array extended_kindlmann(Real x) { template std::array to_8bit_rgba(const std::array &v) { using std::sqrt; - std::array pixel; + std::array pixel {}; for (auto i = 0; i < 3; ++i) { // Apply gamma correction here: Real u = sqrt(v[i]); diff --git a/include/boost/math/tools/condition_numbers.hpp b/include/boost/math/tools/condition_numbers.hpp index a5801996a5..8754d12828 100644 --- a/include/boost/math/tools/condition_numbers.hpp +++ b/include/boost/math/tools/condition_numbers.hpp @@ -94,15 +94,18 @@ Real evaluation_condition_number(F const & f, Real const & x) } bool caught_exception = false; Real fp; +#ifndef BOOST_NO_EXCEPTIONS try { +#endif fp = finite_difference_derivative(f, x); +#ifndef BOOST_NO_EXCEPTIONS } catch(...) { caught_exception = true; } - +#endif if (isnan(fp) || caught_exception) { // Check if the right derivative exists: diff --git a/include/boost/math/tools/config.hpp b/include/boost/math/tools/config.hpp index 2a054fc261..09db298a88 100644 --- a/include/boost/math/tools/config.hpp +++ b/include/boost/math/tools/config.hpp @@ -13,6 +13,17 @@ #include +// Minimum language standard transition +#ifdef _MSVC_LANG +# if _MSVC_LANG < 201402L +# pragma warning("The minimum language standard to use Boost.Math will be C++14 starting in July 2023 (Boost 1.82 release)"); +# endif +#else +# if __cplusplus < 201402L +# warning "The minimum language standard to use Boost.Math will be C++14 starting in July 2023 (Boost 1.82 release)" +# endif +#endif + #ifndef BOOST_MATH_STANDALONE #include @@ -77,8 +88,33 @@ # define BOOST_NO_CXX11_THREAD_LOCAL #endif // BOOST_DISABLE_THREADS +#ifdef __GNUC__ +# if !defined(__EXCEPTIONS) && !defined(BOOST_NO_EXCEPTIONS) +# define BOOST_NO_EXCEPTIONS +# endif + // + // Make sure we have some std lib headers included so we can detect __GXX_RTTI: + // +# include // for min and max +# include +# ifndef __GXX_RTTI +# ifndef BOOST_NO_TYPEID +# define BOOST_NO_TYPEID +# endif +# ifndef BOOST_NO_RTTI +# define BOOST_NO_RTTI +# endif +# endif +#endif + #endif // BOOST_MATH_STANDALONE +// Support compilers with P0024R2 implemented without linking TBB +// https://en.cppreference.com/w/cpp/compiler_support +#if !defined(BOOST_NO_CXX17_HDR_EXECUTION) && defined(BOOST_HAS_THREADS) +# define BOOST_MATH_EXEC_COMPATIBLE +#endif + #include // for min and max #include #include diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index c2528b1da7..17e4362de2 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -360,7 +360,7 @@ class polynomial size_type degree() const { if (size() == 0) - throw std::logic_error("degree() is undefined for the zero polynomial."); + BOOST_MATH_THROW_EXCEPTION(std::logic_error("degree() is undefined for the zero polynomial.")); return m_data.size() - 1; } value_type& operator[](size_type i) diff --git a/include/boost/math/tools/random_vector.hpp b/include/boost/math/tools/random_vector.hpp index 8e3a3289aa..d1954ee3b5 100644 --- a/include/boost/math/tools/random_vector.hpp +++ b/include/boost/math/tools/random_vector.hpp @@ -12,8 +12,8 @@ namespace boost { namespace math { // To stress test, set global_seed = 0, global_size = huge. -static const std::size_t global_seed = 0; -static const std::size_t global_size = 128; +static constexpr std::size_t global_seed = 0; +static constexpr std::size_t global_size = 128; template::value, bool>::type = true> std::vector generate_random_vector(std::size_t size, std::size_t seed) @@ -35,6 +35,28 @@ std::vector generate_random_vector(std::size_t size, std::size_t seed) return v; } +template::value, bool>::type = true> +std::vector generate_random_uniform_vector(std::size_t size, std::size_t seed, T lower_bound = T(0), T upper_bound = T(1)) +{ + if (seed == 0) + { + std::random_device rd; + seed = rd(); + } + std::vector v(size); + + std::mt19937 gen(seed); + + std::uniform_real_distribution dis(lower_bound, upper_bound); + + for (auto& i : v) + { + i = dis(gen); + } + + return v; +} + template::value, bool>::type = true> std::vector generate_random_vector(std::size_t size, std::size_t seed, T mean, T stddev) { diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 5b06e5bc97..0f2129d1f7 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -372,6 +372,8 @@ namespace detail { T bracket_root_towards_max(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { using std::fabs; + if(count < 2) + return guess - (max + min) / 2; // Not enough counts left to do anything!! // // Move guess towards max until we bracket the root, updating min and max as we go: // @@ -427,6 +429,8 @@ namespace detail { T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { using std::fabs; + if (count < 2) + return guess - (max + min) / 2; // Not enough counts left to do anything!! // // Move guess towards min until we bracket the root, updating min and max as we go: // diff --git a/reporting/performance/chatterjee_correlation_performance.cpp b/reporting/performance/chatterjee_correlation_performance.cpp new file mode 100644 index 0000000000..866edb1048 --- /dev/null +++ b/reporting/performance/chatterjee_correlation_performance.cpp @@ -0,0 +1,34 @@ +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include + +using boost::math::generate_random_vector; + +template +void chatterjee_correlation(benchmark::State& state) +{ + constexpr std::size_t seed {}; + const std::size_t size = state.range(0); + std::vector u = generate_random_vector(size, seed); + std::vector v = generate_random_vector(size, seed); + + std::sort(u.begin(), u.end()); + + for (auto _ : state) + { + benchmark::DoNotOptimize(boost::math::statistics::chatterjee_correlation(u, v)); + } + state.SetComplexityN(state.range(0)); +} + +BENCHMARK_TEMPLATE(chatterjee_correlation, float)->RangeMultiplier(2)->Range(1 << 6, 1 << 20)->Complexity()->UseRealTime(); +BENCHMARK_TEMPLATE(chatterjee_correlation, double)->RangeMultiplier(2)->Range(1 << 6, 1 << 20)->Complexity()->UseRealTime(); + +BENCHMARK_MAIN(); diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 5deced8cb0..2708b28096 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -92,8 +92,13 @@ cpp-pch pch_light : pch_light.hpp : ../../test/build//boost_unit_test_frame lib compile_test_main : compile_test/main.cpp ; test-suite special_fun : - [ run test_1F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] ] # hypergeometric_pFq_checked_series.hpp uses auto, the rest are from quadrature tests. - [ run test_2F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] ] + [ run test_1F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] TEST=1 : test_1F0_1 ] # hypergeometric_pFq_checked_series.hpp uses auto, the rest are from quadrature tests. + [ run test_1F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] TEST=2 : test_1F0_2 ] # hypergeometric_pFq_checked_series.hpp uses auto, the rest are from quadrature tests. + [ run test_1F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] TEST=3 : test_1F0_3 ] # hypergeometric_pFq_checked_series.hpp uses auto, the rest are from quadrature tests. + [ run test_2F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] TEST=1 : test_2F0_1 ] + [ run test_2F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] TEST=2 : test_2F0_2 ] + [ run test_2F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] TEST=3 : test_2F0_3 ] + [ run test_2F0.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] TEST=4 : test_2F0_4 ] [ run test_0F1.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] TEST=1 : test_0F1_1 ] [ run test_0F1.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] TEST=2 : test_0F1_2 ] @@ -154,6 +159,7 @@ test-suite special_fun : [ run ccmath_isless_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] [ run ccmath_islessequal_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] [ run ccmath_isunordered_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] + [ run ccmath_fma_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] [ run log1p_expm1_test.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run powm1_sqrtp1m1_test.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run git_issue_705.cpp ../../test/build//boost_unit_test_framework ] @@ -1012,6 +1018,8 @@ test-suite misc : [ run bivariate_statistics_test.cpp : : : [ requires cxx11_hdr_forward_list cxx11_hdr_atomic cxx11_hdr_thread cxx11_hdr_tuple cxx11_hdr_future cxx11_sfinae_expr ] [ check-target-builds ../config//is_cygwin_run "Cygwin CI run" : no ] ] [ run linear_regression_test.cpp : : : [ requires cxx11_hdr_forward_list cxx11_hdr_atomic cxx11_hdr_thread cxx11_hdr_tuple cxx11_hdr_future cxx11_sfinae_expr ] ] [ run test_runs_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] + [ run test_chatterjee_correlation.cpp ../../test/build//boost_unit_test_framework ] + [ run test_rank.cpp ../../test/build//boost_unit_test_framework ] [ run lanczos_smoothing_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run condition_number_test.cpp ../../test/build//boost_unit_test_framework : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : "-Bstatic -lquadmath -Bdynamic" ] [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run test_real_concept.cpp ../../test/build//boost_unit_test_framework ] @@ -1504,3 +1512,132 @@ rule get_float128_tests } test-suite float128_tests : [ get_float128_tests ] ; + +# +# Things that we can test with exceptions and RTTI turned off: +# +alias no_eh_tests : + compl_abs_incl_test + compl_acos_incl_test + compl_acosh_incl_test + compl_asin_incl_test + compl_asinh_incl_test + compl_atan_incl_test + compl_atanh_incl_test + sf_acosh_incl_test + sf_asinh_incl_test + sf_atanh_incl_test + sf_beta_incl_test + sf_bernoulli_incl_test + sf_bessel_incl_test + sf_bessel_deriv_incl_test + sf_binomial_incl_test + sf_cbrt_incl_test + sf_cos_pi_incl_test + sf_digamma_incl_test + sf_polygamma_incl_test + sf_ellint_1_incl_test + sf_ellint_2_incl_test + sf_ellint_3_incl_test + sf_ellint_d_incl_test + sf_jacobi_theta_incl_test + sf_jacobi_zeta_incl_test + sf_heuman_lambda_incl_test + sf_ellint_rc_incl_test + sf_ellint_rd_incl_test + sf_ellint_rf_incl_test + sf_ellint_rj_incl_test + sf_ellint_rg_incl_test + sf_erf_incl_test + sf_expint_incl_test + sf_expm1_incl_test + sf_factorials_incl_test + sf_fpclassify_incl_test + sf_gamma_incl_test + sf_hermite_incl_test + sf_hypot_incl_test + sf_laguerre_incl_test + sf_lanczos_incl_test + sf_legendre_incl_test + #sf_legendre_stieltjes_incl_test + sf_log1p_incl_test + sf_math_fwd_incl_test + sf_modf_incl_test + sf_next_incl_test + sf_powm1_incl_test + sf_prime_incl_test + sf_relative_distance_incl_test + sf_round_incl_test + sf_sign_incl_test + sf_sin_pi_incl_test + sf_sinc_incl_test + sf_sinhc_incl_test + sf_sph_harm_incl_test + sf_sqrt1pm1_incl_test + sf_trunc_incl_test + sf_ulp_incl_test + sf_zeta_incl_test + sf_chebyshev_incl_test + #sf_chebyshev_transform_incl_test + sf_fibonacci_incl_test + #sf_gegenbauer_incl_test + sf_lambert_w_incl_test + sf_nonfinite_num_facets_incl_test + sf_airy_incl_test + sf_hankel_incl_test + sf_jacobi_incl_test + sf_owens_t_incl_test + dist_skew_norm_incl_test + constants_incl_test + quad_trapezoidal_incl_test + test_traits + tools_config_inc_test + tools_fraction_inc_test + tools_minima_inc_test + tools_polynomial_inc_test + tools_precision_inc_test + tools_rational_inc_test + tools_real_cast_inc_test + tools_remez_inc_test + tools_roots_inc_test + tools_series_inc_test + tools_solve_inc_test + tools_stats_inc_test + tools_test_data_inc_test + #tools_test_inc_test + tools_toms748_inc_test + tools_agm_incl_test + tools_assert_incl_test + tools_atomic_incl_test + tools_big_constant_incl_test + #tools_centered_continued_fraction_incl_test + tools_cohen_acceleration_incl_test + tools_complex_incl_test + tools_condition_numbers_incl_test + tools_convert_from_string_incl_test + tools_cxx03_warn_incl_test + #tools_engel_expansion_incl_test + tools_header_deprecated_incl_test + tools_is_detected_incl_test + #tools_luroth_expansion_incl_test + tools_mp_incl_test + tools_norms_incl_test + tools_polynomial_gcd_incl_test + tools_promotion_incl_test + tools_random_vector_incl_test + #tools_simple_continued_fraction_incl_test + tools_test_value_incl_test + tools_throw_exception_incl_test + tools_traits_incl_test + #tools_ulps_plot_incl_test + tools_workaround_incl_test +; + +explicit no_eh_tests ; + +# Some aliases which group blocks of tests for CI testing: + +alias github_ci_block_1 : special_fun float128_tests distribution_tests mp misc ; +alias github_ci_block_2 : quadrature interpolators autodiff ../example//examples ../tools ; +explicit github_ci_block_1 ; +explicit github_ci_block_2 ; diff --git a/test/ccmath_fma_test.cpp b/test/ccmath_fma_test.cpp new file mode 100644 index 0000000000..a5aa74914d --- /dev/null +++ b/test/ccmath_fma_test.cpp @@ -0,0 +1,73 @@ +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_HAS_FLOAT128 +#include +#endif + +#if !defined(BOOST_MATH_NO_CONSTEXPR_DETECTION) && !defined(BOOST_MATH_USING_BUILTIN_CONSTANT_P) +template +constexpr void test() +{ + // Error handling + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::fma(std::numeric_limits::infinity(), T(0), T(1)))); + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::fma(T(0), std::numeric_limits::infinity(), T(1)))); + + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::fma(std::numeric_limits::infinity(), T(0), std::numeric_limits::quiet_NaN()))); + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::fma(T(0), std::numeric_limits::infinity(), std::numeric_limits::quiet_NaN()))); + + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::fma(std::numeric_limits::quiet_NaN(), T(1), T(1)))); + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::fma(T(1), std::numeric_limits::quiet_NaN(), T(1)))); + + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::fma(T(1), T(1), std::numeric_limits::quiet_NaN()))); + + // Functionality + static_assert(boost::math::ccmath::fma(T(1), T(2), T(3)) == T(5)); + static_assert(boost::math::ccmath::fma(T(2), T(3), T(1)) == T(7)); + + // Correct promoted types + if constexpr (!std::is_same_v) + { + constexpr auto test_type = boost::math::ccmath::fma(T(1), 1.0, 1.0f); + static_assert(std::is_same_v>); + } + else + { + constexpr auto test_type = boost::math::ccmath::fma(1.0f, 1, 1.0); + static_assert(std::is_same_v>); + } +} + +int main() +{ + test(); + test(); + + #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test(); + #endif + + #ifdef BOOST_HAS_FLOAT128 + test(); + #endif + + return 0; +} +#else +int main() +{ + return 0; +} +#endif diff --git a/test/compile_test/ccmath_fma_incl_test.cpp b/test/compile_test/ccmath_fma_incl_test.cpp new file mode 100644 index 0000000000..ef035bd8b6 --- /dev/null +++ b/test/compile_test/ccmath_fma_incl_test.cpp @@ -0,0 +1,16 @@ +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include "test_compile_result.hpp" + +void compile_and_link_test() +{ + check_result(boost::math::ccmath::fma(1.0f, 1.0f, 1.0f)); + check_result(boost::math::ccmath::fma(1.0, 1.0, 1.0)); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + check_result(boost::math::ccmath::fma(1.0l, 1.0l, 1.0l)); +#endif +} diff --git a/test/compile_test/stats_chaterjee_incl_test.cpp b/test/compile_test/stats_chaterjee_incl_test.cpp new file mode 100644 index 0000000000..84c235c244 --- /dev/null +++ b/test/compile_test/stats_chaterjee_incl_test.cpp @@ -0,0 +1,9 @@ +// Copyright Matt Borland 2021. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +// Basic sanity check that header +// #includes all the files that it needs to. +// +#include diff --git a/test/compile_test/test_compile_result.hpp b/test/compile_test/test_compile_result.hpp index a86d390fac..0afdced913 100644 --- a/test/compile_test/test_compile_result.hpp +++ b/test/compile_test/test_compile_result.hpp @@ -9,14 +9,12 @@ // detect missing includes). // -static const float f = 0; -static const double d = 0; -static const long double l = 0; -static const unsigned u = 0; -static const int i = 0; - -//template -//inline void check_result_imp(T, T){} +static constexpr float f = 0; +static constexpr double d = 0; +static constexpr long double l = 0; +static constexpr unsigned u = 0; +static constexpr int i = 0; +static constexpr unsigned long li = 1; inline void check_result_imp(float, float){} inline void check_result_imp(double, double){} @@ -38,12 +36,12 @@ inline void check_result_imp(bool, bool){} template struct local_is_same { - enum{ value = false }; + static constexpr bool value = false; }; template struct local_is_same { - enum{ value = true }; + static constexpr bool value = true; }; template @@ -54,7 +52,7 @@ inline void check_result_imp(T1, T2) #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-local-typedefs" #endif - typedef BOOST_MATH_ASSERT_UNUSED_ATTRIBUTE int static_assertion[local_is_same::value ? 1 : 0]; + using static_assertion = int[local_is_same::value ? 1 : 0]; #if defined(__GNUC__) && ((__GNUC__ > 4) || ((__GNUC__ == 4) && (__GNUC_MINOR__ >= 8))) #pragma GCC diagnostic pop #endif @@ -68,23 +66,12 @@ inline void check_result(T2) return check_result_imp(a, b); } -union max_align_type -{ - char c; - short s; - int i; - long l; - double d; - long double ld; - long long ll; -}; - template struct DistributionConcept { static void constraints() { - typedef typename Distribution::value_type value_type; + using value_type = typename Distribution::value_type; const Distribution& dist = DistributionConcept::get_object(); @@ -93,6 +80,7 @@ struct DistributionConcept check_result(cdf(dist, x)); check_result(cdf(complement(dist, x))); check_result(pdf(dist, x)); + check_result(logpdf(dist, x)); check_result(quantile(dist, x)); check_result(quantile(complement(dist, x))); check_result(mean(dist)); @@ -116,6 +104,7 @@ struct DistributionConcept check_result(cdf(dist, f)); check_result(cdf(complement(dist, f))); check_result(pdf(dist, f)); + check_result(logpdf(dist, f)); check_result(quantile(dist, f)); check_result(quantile(complement(dist, f))); check_result(hazard(dist, f)); @@ -123,6 +112,7 @@ struct DistributionConcept check_result(cdf(dist, d)); check_result(cdf(complement(dist, d))); check_result(pdf(dist, d)); + check_result(logpdf(dist, d)); check_result(quantile(dist, d)); check_result(quantile(complement(dist, d))); check_result(hazard(dist, d)); @@ -130,6 +120,7 @@ struct DistributionConcept check_result(cdf(dist, l)); check_result(cdf(complement(dist, l))); check_result(pdf(dist, l)); + check_result(logpdf(dist, l)); check_result(quantile(dist, l)); check_result(quantile(complement(dist, l))); check_result(hazard(dist, l)); @@ -137,20 +128,32 @@ struct DistributionConcept check_result(cdf(dist, i)); check_result(cdf(complement(dist, i))); check_result(pdf(dist, i)); + check_result(logpdf(dist, i)); check_result(quantile(dist, i)); check_result(quantile(complement(dist, i))); check_result(hazard(dist, i)); check_result(chf(dist, i)); - unsigned long li = 1; check_result(cdf(dist, li)); check_result(cdf(complement(dist, li))); check_result(pdf(dist, li)); + check_result(logpdf(dist, li)); check_result(quantile(dist, li)); check_result(quantile(complement(dist, li))); check_result(hazard(dist, li)); check_result(chf(dist, li)); } private: + union max_align_type + { + char c; + short s; + int i; + long l; + double d; + long double ld; + long long ll; + }; + static void* storage() { static max_align_type storage[sizeof(Distribution)]; diff --git a/test/math_unit_test.hpp b/test/math_unit_test.hpp index 6a6af9fc15..0a526a15a3 100644 --- a/test/math_unit_test.hpp +++ b/test/math_unit_test.hpp @@ -384,6 +384,8 @@ int report_errors() #define CHECK_ULP_CLOSE(X, Y, Z) boost::math::test::check_ulp_close((X), (Y), (Z), __FILE__, __func__, __LINE__) +#define CHECK_GE(X, Y) boost::math::test::check_le((Y), (X), __FILE__, __func__, __LINE__) + #define CHECK_LE(X, Y) boost::math::test::check_le((X), (Y), __FILE__, __func__, __LINE__) #define CHECK_NAN(X) boost::math::test::check_nan((X), __FILE__, __func__, __LINE__) diff --git a/test/multiprc_concept_check_2.cpp b/test/multiprc_concept_check_2.cpp index 500e1fc4f3..30311fda24 100644 --- a/test/multiprc_concept_check_2.cpp +++ b/test/multiprc_concept_check_2.cpp @@ -35,7 +35,8 @@ void foo() int main() { - BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept)); + boost::math::concepts::RealTypeConcept checker; + checker.constraints(); } diff --git a/test/test_1F0.cpp b/test/test_1F0.cpp index f99881da70..c32f253928 100644 --- a/test/test_1F0.cpp +++ b/test/test_1F0.cpp @@ -11,6 +11,7 @@ BOOST_AUTO_TEST_CASE( test_main ) { +#if !defined(TEST) || (TEST == 1) #ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS test_spots(0.0F); #endif @@ -21,8 +22,13 @@ BOOST_AUTO_TEST_CASE( test_main ) test_spots(boost::math::concepts::real_concept(0.1)); #endif #endif +#endif +#if !defined(TEST) || (TEST == 2) test_spots(boost::multiprecision::cpp_bin_float_quad()); +#endif +#if !defined(TEST) || (TEST == 3) test_spots(boost::multiprecision::cpp_dec_float_50()); +#endif } diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp index 2e7805d442..5edf8cdc34 100644 --- a/test/test_1F1.hpp +++ b/test/test_1F1.hpp @@ -162,7 +162,7 @@ void test_spots5(T, const char* type_name) template void test_spots6(T, const char* type_name) { - static const std::array, 91> hypergeometric_1F1_bugs = { { + static const std::array, 93> hypergeometric_1F1_bugs = { { { { static_cast(17955.561660766602), static_cast(9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }}, { { static_cast(17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }}, { { static_cast(-17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(82.406154185533524), SC_(-42897094853118832762870100.8669248353530950866) }} , @@ -287,6 +287,11 @@ void test_spots6(T, const char* type_name) { { (T)std::ldexp((double)-9751199809536000, -45), (T)std::ldexp((double)-17654191685632000, -47), (T)std::ldexp((double)10587451850752000, -47), SC_(-1.9601415510439595625538337964298353914980331018955e+68) }}, { { (T)std::ldexp((double)-15233620754432000, -45), (T)std::ldexp((double)-12708283072512000, -46), (T)std::ldexp((double)10255461007360000, -46), SC_(-5.4344106361679075861859567858016187271235441673635e+125) }}, { { (T)std::ldexp((double)-11241354149888000, -45), (T)std::ldexp((double)-9580579905536000, -45), (T)std::ldexp((double)12224976846848000, -47), SC_(12046856548470067405870726490464935201150430438.035) }}, + // + // Bugs found while testing color maps: + // + { { SC_(0.078125000000000000), SC_(-0.039062500000000000), SC_(0.5), SC_(-0.3371910410915676603577770246237158427221) }}, + { { SC_(-19.492187500000000), SC_(0.50781250000000000), SC_(0.5), SC_(1.2551298228307647570646714060395253358015) }}, } }; static const std::array, 2> hypergeometric_1F1_big_bugs = { { #if DBL_MAX_EXP == LDBL_MAX_EXP diff --git a/test/test_2F0.cpp b/test/test_2F0.cpp index 66b2fd9fb1..3726a61840 100644 --- a/test/test_2F0.cpp +++ b/test/test_2F0.cpp @@ -82,21 +82,33 @@ BOOST_AUTO_TEST_CASE( test_main ) { expected_results(); BOOST_MATH_CONTROL_FP; - +#if !defined(TEST) || (TEST == 1) #ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS test_spots(0.0F, "float"); #endif test_spots(0.0, "double"); +#endif +#if !defined(TEST) || (TEST == 2) #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_spots(0.0L, "long double"); #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS test_spots(boost::math::concepts::real_concept(0.1), "real_concept"); #endif #endif +#endif + +#if defined(__GNUC__) && (__GNUC__ == 12) + // gcc-12 runs the machine out of memory: +#define BOOST_MATH_NO_MP_TESTS +#endif #ifndef BOOST_MATH_NO_MP_TESTS using dec_40 = boost::multiprecision::number>; +#if !defined(TEST) || (TEST == 3) test_spots(boost::multiprecision::cpp_bin_float_quad(), "cpp_bin_float_quad"); +#endif +#if !defined(TEST) || (TEST == 4) test_spots(dec_40(), "dec_40"); #endif +#endif } diff --git a/test/test_arcsine.cpp b/test/test_arcsine.cpp index 68036686e5..0c2d847a96 100644 --- a/test/test_arcsine.cpp +++ b/test/test_arcsine.cpp @@ -105,6 +105,14 @@ void test_ignore_policy(RealType) BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(0, 1), static_cast (+2)))); // x > x_max BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(-1, 1), static_cast (+2)))); // x > x_max + // Logpdf + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), std::numeric_limits::infinity()))); // x == infinity + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), std::numeric_limits::infinity()))); // x == infinity + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), static_cast (-2)))); // x < xmin + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), static_cast (-2)))); // x < xmin + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), static_cast (+2)))); // x > x_max + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), static_cast (+2)))); // x > x_max + // Mean BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(-nan, 0)))); BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(+nan, 0)))); @@ -239,6 +247,7 @@ void test_spots(RealType) using boost::math::arcsine_distribution; using ::boost::math::cdf; using ::boost::math::pdf; + using ::boost::math::logpdf; using ::boost::math::complement; using ::boost::math::quantile; @@ -292,6 +301,16 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION(pdf(arcsine_01, static_cast(1) - tolerance), 1 /(sqrt(tolerance) * boost::math::constants::pi()), 2 * tolerance); // + // Log PDF + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.000001), static_cast(5.7630258931329868780772138043668005779060097243996L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.000005), static_cast(4.9583089369219367114435788047327747268154560240604L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.05), static_cast(0.37878289812137058928728250884555529541061717942415L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.5), static_cast(-0.45158270528945486472619522989488214357179467855506L), tolerance); + // Note loss of significance when x is near x_max. + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.95), static_cast(0.37878289812137058928728250884555529541061717942415L), 8 * tolerance); // Less accurate. + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.999995), static_cast(4.9583089369219367114435788047327747268154560240604L), 50000 * tolerance); // Much less accurate. + BOOST_CHECK_CLOSE_FRACTION(logpdf(arcsine_01, 0.999999), static_cast(5.7630258931329868780772138043668005779060097243996L), 100000 * tolerance);// Even less accurate. + // CDF BOOST_CHECK_CLOSE_FRACTION(cdf(arcsine_01, 0.000001), static_cast(0.00063661987847092448418377367957384866092127786060574L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(arcsine_01, 0.000005), static_cast(0.0014235262731079289297302426454125318201831474507326L), tolerance); @@ -353,6 +372,10 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.5), static_cast(0.36755259694786136634088433220864629426492432024443L), tolerance); BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.95), static_cast(1.0194074882503562519812229448639426942621591013381L), 2 * tolerance); // Less accurate. + BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.05), static_cast(-1.1434783207403409089630164813372974217316704642782L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.5), static_cast(-1.0008888496235097104238178483561449958955399574664L), tolerance); + BOOST_CHECK_CLOSE_FRACTION(logpdf(as_m11, 0.95), static_cast(0.019221564639767605567429885545559909302927558782238L), 100 * tolerance); // Less accurate. + BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.05), static_cast(0.51592213323666034437274347433261364289389772737836L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.5), static_cast(0.66666666666666666666666666666666666666666666666667L), 2 * tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m11, 0.95), static_cast(0.89891737589574013042121018491729701360300248368629L), tolerance); // Not less accurate. @@ -445,6 +468,31 @@ void test_spots(RealType) arcsine_distribution(static_cast(0), static_cast(1)), // bad x > 1. static_cast(999)), std::domain_error); + BOOST_MATH_CHECK_THROW( // For various bad arguments. + logpdf( + arcsine_distribution(static_cast(+1), static_cast(-1)), // min_x > max_x + static_cast(1)), std::domain_error); + + BOOST_MATH_CHECK_THROW( + logpdf( + arcsine_distribution(static_cast(1), static_cast(0)), // bad constructor parameters. + static_cast(1)), std::domain_error); + + BOOST_MATH_CHECK_THROW( + logpdf( + arcsine_distribution(static_cast(1), static_cast(-1)), // bad constructor parameters. + static_cast(1)), std::domain_error); + + BOOST_MATH_CHECK_THROW( + logpdf( + arcsine_distribution(static_cast(1), static_cast(1)), // equal constructor parameters. + static_cast(-1)), std::domain_error); + + BOOST_MATH_CHECK_THROW( + logpdf( + arcsine_distribution(static_cast(0), static_cast(1)), // bad x > 1. + static_cast(999)), std::domain_error); + // Checks on things that are errors. // Construction with 'bad' parameters. @@ -453,6 +501,7 @@ void test_spots(RealType) arcsine_distribution<> dist; BOOST_MATH_CHECK_THROW(pdf(dist, -1), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(dist, -1), std::domain_error); BOOST_MATH_CHECK_THROW(cdf(dist, -1), std::domain_error); BOOST_MATH_CHECK_THROW(cdf(complement(dist, -1)), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(dist, -1), std::domain_error); @@ -463,6 +512,8 @@ void test_spots(RealType) // Various combinations of bad constructor and member function parameters. BOOST_MATH_CHECK_THROW(pdf(boost::math::arcsine_distribution(0, 1), -1), std::domain_error); BOOST_MATH_CHECK_THROW(pdf(boost::math::arcsine_distribution(-1, 1), +2), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(boost::math::arcsine_distribution(0, 1), -1), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(boost::math::arcsine_distribution(-1, 1), +2), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(boost::math::arcsine_distribution(1, 1), -1), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(boost::math::arcsine_distribution(1, 1), 2), std::domain_error); @@ -484,6 +535,7 @@ void test_spots(RealType) arcsine_distribution w(RealType(-1), RealType(+1)); // NaN parameters to member functions should throw. BOOST_MATH_CHECK_THROW(pdf(w, +nan), std::domain_error); // x = NaN + BOOST_MATH_CHECK_THROW(logpdf(w, +nan), std::domain_error); // x = NaN BOOST_MATH_CHECK_THROW(cdf(w, +nan), std::domain_error); // x = NaN BOOST_MATH_CHECK_THROW(cdf(complement(w, +nan)), std::domain_error); // x = + nan BOOST_MATH_CHECK_THROW(quantile(w, +nan), std::domain_error); // p = + nan @@ -511,6 +563,7 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(arcsine_distribution(1, inf), std::domain_error); #endif BOOST_MATH_CHECK_THROW(pdf(w, +inf), std::domain_error); // x = inf + BOOST_MATH_CHECK_THROW(logpdf(w, +inf), std::domain_error); // x = inf BOOST_MATH_CHECK_THROW(cdf(w, +inf), std::domain_error); // x = inf BOOST_MATH_CHECK_THROW(cdf(complement(w, +inf)), std::domain_error); // x = + inf BOOST_MATH_CHECK_THROW(quantile(w, +inf), std::domain_error); // p = + inf diff --git a/test/test_chatterjee_correlation.cpp b/test/test_chatterjee_correlation.cpp new file mode 100644 index 0000000000..d39a1a33a2 --- /dev/null +++ b/test/test_chatterjee_correlation.cpp @@ -0,0 +1,159 @@ +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "math_unit_test.hpp" + +// The Chatterjee correlation is invariant under: +// - Shuffles. (X_i, Y_i) -> (X_sigma(i), Y_sigma(i)), where sigma is a permutation. +// - Strictly monotone transformations: (X_i, Y_i) -> (f(X_i), g(Y_i)) where f' > 0 and g' > 0. +// + +using boost::math::statistics::chatterjee_correlation; + +template +void properties() +{ + std::size_t vector_size = 256; + std::mt19937_64 mt(123521); + std::uniform_real_distribution unif(-1, 1); + std::vector X(vector_size); + std::vector Y(vector_size); + + for (std::size_t i = 0; i < vector_size; ++i) + { + X[i] = unif(mt); + Y[i] = unif(mt); + } + std::sort(X.begin(), X.end()); + Real coeff1 = chatterjee_correlation(X, Y); + // The minimum possible value of En(X, Y) is -1/2 + O(1/n) + CHECK_GE(coeff1, Real(-0.5)); + CHECK_LE(coeff1, Real(1)); + + // Now apply a monotone function to the data + for (std::size_t i = 0; i < vector_size; ++i) + { + X[i] = Real(2.3)*X[i] - Real(7.3); + Y[i] = Real(7.6)*Y[i] - Real(8.6); + } + auto coeff3 = chatterjee_correlation(X, Y); + CHECK_EQUAL(coeff1, coeff3); + + // If there are no ties among the Yis, the maximum possible value of Xi(X, Y) is (n - 2)/(n + 1), which is attained if Yi = Xi for all i + auto coeff = chatterjee_correlation(X, X); + // These floating point numbers are computed by two different methods, so we can expect some floating point error: + const auto n = X.size(); + CHECK_ULP_CLOSE(coeff, Real(n-2)/Real(n+1), 1); + std::sort(Y.begin(), Y.end()); + coeff = chatterjee_correlation(Y, Y); + CHECK_ULP_CLOSE(coeff, Real(n-2)/Real(n+1), 1); +} + +template +void test_spots() +{ + // Rank Order: Result will be 1 - 3*3 / (4^2 - 1) = 1 - 9/15 = 0.6 + std::vector x = {1, 2, 3, 4}; + std::vector y = {1, 2, 3, 4}; + CHECK_ULP_CLOSE(chatterjee_correlation(x, y), 1 - Real(9)/15, 1); + + // Reverse rank order should be the same as above + y = {4, 3, 2, 1}; + CHECK_ULP_CLOSE(chatterjee_correlation(x, y), 1 - Real(9)/15, 1); + + // Alternating order: 1 - 3*5 / (4^2 - 1) = 1 - 15/15 = 0 + y = {1, 3, 2, 4}; + CHECK_ULP_CLOSE(chatterjee_correlation(x, y), Real(0), 1); + + // All ties will yield quiet NaN + y = {1, 1, 1, 1}; + CHECK_NAN(chatterjee_correlation(x, y)); +} + +#ifdef BOOST_MATH_EXEC_COMPATIBLE + +template +void test_threaded(ExecutionPolicy&& exec) +{ + std::vector x = boost::math::generate_random_vector(1024, 0); + std::vector y = boost::math::generate_random_vector(1024, 1); + + std::sort(std::forward(exec), x.begin(), x.end()); + + auto seq_ans = chatterjee_correlation(x, y); + auto par_ans = chatterjee_correlation(exec, x, y); + + CHECK_ULP_CLOSE(seq_ans, par_ans, 1); +}; + +#endif // BOOST_MATH_EXEC_COMPATIBLE + +template +void test_paper() +{ + constexpr Real two_pi = boost::math::constants::two_pi(); + + // Page 9 figure (a) y = x + std::vector x = boost::math::generate_random_uniform_vector(100, 0, -two_pi, two_pi); + std::sort(x.begin(), x.end()); + auto result = chatterjee_correlation(x, x); + CHECK_MOLLIFIED_CLOSE(result, Real(0.970), 0.005); + + // Page 9 figure (d) y = x^2 + std::vector y = x; + for (auto& i : y) + { + i *= i; + } + + result = chatterjee_correlation(x, y); + CHECK_MOLLIFIED_CLOSE(result, Real(0.941), 0.005); + + // Page 9 figure (g) y = sin(x) + for (std::size_t i {}; i < x.size(); ++i) + { + y[i] = std::sin(x[i]); + } + + result = chatterjee_correlation(x, y); + CHECK_MOLLIFIED_CLOSE(result, Real(0.885), 0.012); +} + +int main(void) +{ + properties(); + properties(); + properties(); + + test_spots(); + test_spots(); + test_spots(); + + #ifdef BOOST_MATH_EXEC_COMPATIBLE + + test_threaded(std::execution::par); + test_threaded(std::execution::par); + test_threaded(std::execution::par); + test_threaded(std::execution::par_unseq); + test_threaded(std::execution::par_unseq); + test_threaded(std::execution::par_unseq); + + #endif // BOOST_MATH_EXEC_COMPATIBLE + + test_paper(); + test_paper(); + test_paper(); + + return boost::math::test::report_errors(); +} diff --git a/test/test_chi_squared.cpp b/test/test_chi_squared.cpp index 41e4dd5e50..cc7747a6c0 100644 --- a/test/test_chi_squared.cpp +++ b/test/test_chi_squared.cpp @@ -36,6 +36,8 @@ using std::cout; using std::endl; #include using std::numeric_limits; +#include +using std::log; template RealType naive_pdf(RealType df, RealType x) @@ -59,6 +61,8 @@ void test_spot( cdf(dist, cs), P, tol); BOOST_CHECK_CLOSE( pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), cs), tol); + BOOST_CHECK_CLOSE( + logpdf(dist, cs), log(pdf(dist, cs)), tol); if((P < 0.99) && (Q < 0.99)) { // diff --git a/test/test_exponential_dist.cpp b/test/test_exponential_dist.cpp index c0423bf0de..12d9fcad19 100644 --- a/test/test_exponential_dist.cpp +++ b/test/test_exponential_dist.cpp @@ -189,6 +189,31 @@ void test_spots(RealType T) static_cast(9.0799859524969703071183031121101e-5L), // probability. tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + exponential_distribution(0.5), + static_cast(0.125)), // x + log(static_cast(0.46970653140673789305985541231115L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + exponential_distribution(0.5), + static_cast(5)), // x + log(static_cast(0.04104249931194939758476433723358L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + exponential_distribution(2), + static_cast(0.125)), // x + log(static_cast(1.5576015661428097364903405339566L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + exponential_distribution(2), + static_cast(5)), // x + log(static_cast(9.0799859524969703071183031121101e-5L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( ::boost::math::mean( exponential_distribution(2)), diff --git a/test/test_extreme_value.cpp b/test/test_extreme_value.cpp index f61731647e..2493f466df 100644 --- a/test/test_extreme_value.cpp +++ b/test/test_extreme_value.cpp @@ -120,6 +120,25 @@ void test_spots(RealType) static_cast(0.11522236828583456431277265757312L), // probability. tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + extreme_value_distribution(0.5, 2), + static_cast(0.125)), // x + log(static_cast(0.18052654830890205978204427757846L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + extreme_value_distribution(1, 3), + static_cast(5)), // x + log(static_cast(0.0675057324099851209129017326286L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + extreme_value_distribution(1, 3), + static_cast(0)), // x + log(static_cast(0.11522236828583456431277265757312L)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( ::boost::math::mean( extreme_value_distribution(2, 3)), diff --git a/test/test_gamma_dist.cpp b/test/test_gamma_dist.cpp index 682535c283..b7776c79cb 100644 --- a/test/test_gamma_dist.cpp +++ b/test/test_gamma_dist.cpp @@ -88,6 +88,14 @@ void check_gamma(RealType shape, RealType scale, RealType x, RealType p, RealTyp x), // random variable. NaivePDF(shape, scale, x), // PDF tol); // %tolerance. + + // LOGPDF: + BOOST_CHECK_CLOSE( + boost::math::logpdf( + gamma_distribution(shape, scale), // distribution. + x), // random variable. + log(boost::math::pdf(gamma_distribution(shape, scale), x)), // PDF + tol); // %tolerance. } template diff --git a/test/test_hyperexponential_dist.cpp b/test/test_hyperexponential_dist.cpp index 7a30e77de2..c7d3c9843b 100644 --- a/test/test_hyperexponential_dist.cpp +++ b/test/test_hyperexponential_dist.cpp @@ -374,6 +374,15 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(special_cases, RealT, test_types) BOOST_CHECK_CLOSE(boost::math::mode(hexp2), boost::math::mode(exp2), tol); } +// Test C++20 ranges (Currently only GCC10 has full support to P0896R4) +#if (__cplusplus > 202000L || _MSVC_LANG > 202000L) && __has_include() && __GNUC__ >= 10 +// Support for ranges is broken using gcc 11.1 +#if __GNUC__ != 11 +#include +#include +#endif +#endif + BOOST_AUTO_TEST_CASE_TEMPLATE(error_cases, RealT, test_types) { typedef boost::math::hyperexponential_distribution dist_t; @@ -393,8 +402,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(error_cases, RealT, test_types) #if (__cplusplus > 202000L || _MSVC_LANG > 202000L) && __has_include() && __GNUC__ >= 10 // Support for ranges is broken using gcc 11.1 #if __GNUC__ != 11 - #include - #include std::array probs_array {1,2}; std::array rates_array {1,2,3}; diff --git a/test/test_inverse_gamma_distribution.cpp b/test/test_inverse_gamma_distribution.cpp index bad1bb158b..68b238fbc8 100644 --- a/test/test_inverse_gamma_distribution.cpp +++ b/test/test_inverse_gamma_distribution.cpp @@ -72,6 +72,9 @@ void test_spot( BOOST_CHECK_CLOSE_FRACTION( // Compare to naive formula (might be less accurate). pdf(dist, x), naive_pdf(dist.shape(), dist.scale(), x), tol); + BOOST_CHECK_CLOSE_FRACTION( // Compare direct logpdf to naive log(pdf()) + logpdf(dist, x), log(pdf(dist,x)), tol); + BOOST_CHECK_CLOSE_FRACTION( // Compare to expected CDF. cdf(dist, x), P, tol); diff --git a/test/test_inverse_gaussian.cpp b/test/test_inverse_gaussian.cpp index 2bb5303122..68012d48a2 100644 --- a/test/test_inverse_gaussian.cpp +++ b/test/test_inverse_gaussian.cpp @@ -36,6 +36,8 @@ using std::endl; using std::setprecision; #include using std::numeric_limits; +#include +using std::log; template void check_inverse_gaussian(RealType mean, RealType scale, RealType x, RealType p, RealType q, RealType tol) @@ -207,21 +209,29 @@ BOOST_AUTO_TEST_CASE( test_main ) // formatC(SuppDists::dinverse_gaussian(1, 1, 1), digits=17) ... BOOST_CHECK_CLOSE_FRACTION( // x = 1 pdf(w11, 1.), static_cast(0.3989422804014327), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( // x = 1 + logpdf(w11, 1.), static_cast(log(0.3989422804014327)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 1.), static_cast(0.66810200122317065), 10 * tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 0.1), static_cast(0.21979480031862672), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 0.1), static_cast(log(0.21979480031862672)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 0.1), static_cast(0.0040761113207110162), 10 * tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( // small x pdf(w11, 0.01), static_cast(2.0811768202028392e-19), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( // small x + logpdf(w11, 0.01), static_cast(log(2.0811768202028392e-19)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 0.01), static_cast(4.122313403318778e-23), 10 * tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( // smaller x pdf(w11, 0.001), static_cast(2.4420044378793562e-213), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( // smaller x + logpdf(w11, 0.001), static_cast(log(2.4420044378793562e-213)), tolfeweps); // pdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 0.001), static_cast(4.8791443010851493e-219), 1000 * tolfeweps); // cdf // 4.8791443010859224e-219 versus 4.8791443010851493e-219 so still 14 decimal digits. @@ -240,25 +250,35 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 0.5), static_cast(0.87878257893544476), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 0.5), static_cast(log(0.87878257893544476)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 0.5), static_cast(0.3649755481729598), tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 2), static_cast(0.10984782236693059), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 2), static_cast(log(0.10984782236693059)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 2), static_cast(.88547542598600637), tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 10), static_cast(0.00021979480031862676), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 10), static_cast(log(0.00021979480031862676)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 10), static_cast(0.99964958546279115), tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 100), static_cast(2.0811768202028246e-25), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 100), static_cast(log(2.0811768202028246e-25)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 100), static_cast(1), tolfeweps); // cdf BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 1000), static_cast(2.4420044378793564e-222), 10 * tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w11, 1000), static_cast(log(2.4420044378793564e-222)), 10 * tolfeweps); // pdf BOOST_CHECK_CLOSE_FRACTION( cdf(w11, 1000), static_cast(1.), tolfeweps); // cdf @@ -295,26 +315,37 @@ BOOST_AUTO_TEST_CASE( test_main ) // =================== BOOST_CHECK_CLOSE_FRACTION( // formatC(SuppDists::dinvGauss(1, 2, 3), digits=17) "0.47490884963330904" pdf(w23, 1.), static_cast(0.47490884963330904), tolfeweps ); // pdf - + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w23, 1.), static_cast(log(0.47490884963330904)), tolfeweps ); // logpdf BOOST_CHECK_CLOSE_FRACTION( pdf(w23, 0.1), static_cast(2.8854207087665401e-05), tolfeweps * 2); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w23, 0.1), static_cast(log(2.8854207087665401e-05)), tolfeweps * 2); // logpdf //2.8854207087665452e-005 2.8854207087665401e-005 BOOST_CHECK_CLOSE_FRACTION( pdf(w23, 10.), static_cast(0.0019822751498574636), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w23, 10.), static_cast(log(0.0019822751498574636)), tolfeweps); // logpdf BOOST_CHECK_CLOSE_FRACTION( pdf(w23, 10.), static_cast(0.0019822751498574636), tolfeweps); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w23, 10.), static_cast(log(0.0019822751498574636)), tolfeweps); // logpdf // Bigger changes in mean and scale. inverse_gaussian w012(0.1, 2); BOOST_CHECK_CLOSE_FRACTION( pdf(w012, 1.), static_cast(3.7460367141230404e-36), tolfeweps ); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w012, 1.), static_cast(log(3.7460367141230404e-36)), tolfeweps ); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w012, 1.), static_cast(1), tolfeweps ); // pdf inverse_gaussian w0110(0.1, 10); BOOST_CHECK_CLOSE_FRACTION( pdf(w0110, 1.), static_cast(1.6279643678071011e-176), 100 * tolfeweps ); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w0110, 1.), static_cast(log(1.6279643678071011e-176)), 100 * tolfeweps ); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w0110, 1.), static_cast(1), tolfeweps ); // cdf BOOST_CHECK_CLOSE_FRACTION( @@ -323,6 +354,8 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_CLOSE_FRACTION( pdf(w0110, 0.1), static_cast(39.894228040143268), tolfeweps ); // pdf + BOOST_CHECK_CLOSE_FRACTION( + logpdf(w0110, 0.1), static_cast(log(39.894228040143268)), tolfeweps ); // logpdf BOOST_CHECK_CLOSE_FRACTION( cdf(w0110, 0.1), static_cast(0.51989761564832704), 10 * tolfeweps ); // cdf diff --git a/test/test_laplace.cpp b/test/test_laplace.cpp index 9167547f91..716c738736 100644 --- a/test/test_laplace.cpp +++ b/test/test_laplace.cpp @@ -67,6 +67,8 @@ Test 8: test_extreme_function_arguments() #include #include "test_out_of_range.hpp" using boost::math::laplace_distribution; +#include +using std::log; /* #include @@ -211,6 +213,12 @@ void test_pdf_cdf_ocatave() static_cast(0.067667641618306345946999747486242201703815773119812L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(-2.L)), + // static_cast(0.06766764161831L), + log(static_cast(0.067667641618306345946999747486242201703815773119812L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(-2.L)), //static_cast(0.06766764161831L), @@ -223,6 +231,12 @@ void test_pdf_cdf_ocatave() static_cast(0.18393972058572116079776188508073043372290556554506L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(-1.L)), + //static_cast(0.18393972058572L), + log(static_cast(0.18393972058572116079776188508073043372290556554506L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(-1.L)), static_cast(0.18393972058572L), @@ -245,6 +259,11 @@ void test_pdf_cdf_ocatave() static_cast(0.5L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(0.0L)), + log(static_cast(0.5L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(0.0L)), static_cast(0.5L), @@ -256,6 +275,11 @@ void test_pdf_cdf_ocatave() static_cast(0.30326532985631671180189976749559022672095906778368L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(0.5L)), + log(static_cast(0.30326532985631671180189976749559022672095906778368L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(0.5L)), // static_cast(0.69673467014368L), @@ -268,6 +292,12 @@ void test_pdf_cdf_ocatave() static_cast(0.18393972058572116079776188508073043372290556554506L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(1.0L)), + // static_cast(0.18393972058572L), + log(static_cast(0.18393972058572116079776188508073043372290556554506L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(1.00000000000000L)), // static_cast(0.81606027941428L), @@ -280,6 +310,12 @@ void test_pdf_cdf_ocatave() static_cast(0.067667641618306345946999747486242201703815772944649L), tolerance); + BOOST_CHECK_CLOSE( + logpdf(laplace_distribution(), static_cast(2.0L)), + // static_cast(0.06766764161831L), + log(static_cast(0.067667641618306345946999747486242201703815772944649L)), + tolerance); + BOOST_CHECK_CLOSE( cdf(laplace_distribution(), static_cast(2.0L)), // static_cast(0.93233235838169L), diff --git a/test/test_normal.cpp b/test/test_normal.cpp index b8d0024872..ef984d5e63 100644 --- a/test/test_normal.cpp +++ b/test/test_normal.cpp @@ -41,6 +41,8 @@ using std::setprecision; #include using std::numeric_limits; +#include + using std::log; template RealType NaivePDF(RealType mean, RealType sd, RealType x) @@ -126,6 +128,7 @@ void test_spots(RealType) { // No longer allow x to be NaN, then these tests should throw. BOOST_MATH_CHECK_THROW(pdf(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // x = NaN + BOOST_MATH_CHECK_THROW(logpdf(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // x = NaN BOOST_MATH_CHECK_THROW(cdf(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // x = NaN BOOST_MATH_CHECK_THROW(cdf(complement(N01, +std::numeric_limits::quiet_NaN())), std::domain_error); // x = + infinity BOOST_MATH_CHECK_THROW(quantile(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // p = + infinity @@ -215,6 +218,31 @@ void test_spots(RealType) static_cast(0.3989422804014326779399460599343818684759L / 5), tolerance); + // + // Tests for logpdf + // + RealType temp_tol = tolerance; + + BOOST_IF_CONSTEXPR (std::is_same::value) + { + tolerance *= 100; + } + + BOOST_CHECK_CLOSE( + logpdf(normal_distribution(), static_cast(0)), + log(static_cast(0.3989422804014326779399460599343818684759L)), // 1/sqrt(2*pi) + tolerance); + BOOST_CHECK_CLOSE( + logpdf(normal_distribution(3), static_cast(3)), + log(static_cast(0.3989422804014326779399460599343818684759L)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(normal_distribution(3, 5), static_cast(3)), + log(static_cast(0.3989422804014326779399460599343818684759L / 5)), + tolerance); + + tolerance = temp_tol; + // // Spot checks for mean = -5, sd = 6: // @@ -307,6 +335,8 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(pdf(normal_distribution(0, 0), 0), std::domain_error); BOOST_MATH_CHECK_THROW(pdf(normal_distribution(0, -1), 0), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(normal_distribution(0, 0), 0), std::domain_error); + BOOST_MATH_CHECK_THROW(logpdf(normal_distribution(0, -1), 0), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(normal_distribution(0, 1), -1), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(normal_distribution(0, 1), 2), std::domain_error); } // template void test_spots(RealType) diff --git a/test/test_poisson.cpp b/test/test_poisson.cpp index 934395fe8e..9b75ce162f 100644 --- a/test/test_poisson.cpp +++ b/test/test_poisson.cpp @@ -208,6 +208,31 @@ void test_spots(RealType) static_cast(20)), // K>> mean static_cast(8.277463646553730E-009), // probability. tolerance); + + // LOGPDF + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(4)), // mean 4. + static_cast(0)), + log(static_cast(1.831563888873410E-002)), // probability. + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(4)), // mean 4. + static_cast(2)), + log(static_cast(1.465251111098740E-001)), // probability. + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(20)), // mean big. + static_cast(1)), // k small + log(static_cast(4.122307244877130E-008)), // probability. + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(4)), // mean 4. + static_cast(20)), // K>> mean + log(static_cast(8.277463646553730E-009)), // probability. + tolerance); // CDF BOOST_CHECK_CLOSE( diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index e13ceb6196..41a3a511de 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -7,13 +7,17 @@ #define BOOST_TEST_MAIN #include #include +#ifndef BOOST_MATH_STANDALONE #include +#endif #include #include #include +#ifndef BOOST_MATH_STANDALONE #include #include #include +#endif #include #include #include @@ -271,7 +275,7 @@ typedef boost::mpl::list large_integral_test_types; typedef boost::mpl::list<> mp_integral_test_types; #elif defined(TEST2) typedef boost::mpl::list< -#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) +#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) && !defined(BOOST_MATH_STANDALONE) boost::multiprecision::cpp_int #endif > integral_test_types; @@ -287,13 +291,13 @@ typedef large_integral_test_types mp_integral_test_types; typedef boost::mpl::list non_integral_test_types; #elif defined(TEST2) typedef boost::mpl::list< -#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) +#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) && !defined(BOOST_MATH_STANDALONE) boost::multiprecision::cpp_rational #endif > non_integral_test_types; #elif defined(TEST3) typedef boost::mpl::list< -#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) +#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) && !defined(BOOST_MATH_STANDALONE) boost::multiprecision::cpp_bin_float_single, boost::multiprecision::cpp_dec_float_50 #endif > non_integral_test_types; diff --git a/test/test_rank.cpp b/test/test_rank.cpp new file mode 100644 index 0000000000..73c022121d --- /dev/null +++ b/test/test_rank.cpp @@ -0,0 +1,102 @@ +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include "math_unit_test.hpp" + +template +void test() +{ + std::vector test_vals {T(1.0), T(3.2), T(2.4), T(5.6), T(4.1)}; + auto rank_vector = boost::math::statistics::detail::rank(test_vals.begin(), test_vals.end()); + + CHECK_EQUAL(static_cast(0), rank_vector[0]); + CHECK_EQUAL(static_cast(2), rank_vector[1]); + CHECK_EQUAL(static_cast(1), rank_vector[2]); + CHECK_EQUAL(static_cast(4), rank_vector[3]); + CHECK_EQUAL(static_cast(3), rank_vector[4]); + + // Remove duplicates + test_vals.push_back(T(4.1)); + test_vals.push_back(T(2.4)); + rank_vector = boost::math::statistics::detail::rank(test_vals.begin(), test_vals.end()); + + // Check the size is correct and the ordering is not disrupted + CHECK_EQUAL(static_cast(5), rank_vector.size()); + CHECK_EQUAL(static_cast(0), rank_vector[0]); + CHECK_EQUAL(static_cast(2), rank_vector[1]); + CHECK_EQUAL(static_cast(1), rank_vector[2]); + CHECK_EQUAL(static_cast(4), rank_vector[3]); + CHECK_EQUAL(static_cast(3), rank_vector[4]); +} + +template +void container_test() +{ + std::vector test_vals {T(1.0), T(3.2), T(2.4), T(5.6), T(4.1)}; + auto rank_vector = boost::math::statistics::detail::rank(test_vals); + + CHECK_EQUAL(static_cast(0), rank_vector[0]); + CHECK_EQUAL(static_cast(2), rank_vector[1]); + CHECK_EQUAL(static_cast(1), rank_vector[2]); + CHECK_EQUAL(static_cast(4), rank_vector[3]); + CHECK_EQUAL(static_cast(3), rank_vector[4]); +} + +#ifdef BOOST_MATH_EXEC_COMPATIBLE + +#include + +template +void execution_test(ExecutionPolicy&& exec) +{ + std::vector test_vals {T(1.0), T(3.2), T(2.4), T(5.6), T(4.1)}; + auto rank_vector = boost::math::statistics::detail::rank(exec, test_vals.begin(), test_vals.end()); + + CHECK_EQUAL(static_cast(0), rank_vector[0]); + CHECK_EQUAL(static_cast(2), rank_vector[1]); + CHECK_EQUAL(static_cast(1), rank_vector[2]); + CHECK_EQUAL(static_cast(4), rank_vector[3]); + CHECK_EQUAL(static_cast(3), rank_vector[4]); + + // Remove duplicates + test_vals.push_back(T(4.1)); + test_vals.push_back(T(2.4)); + rank_vector = boost::math::statistics::detail::rank(exec, test_vals.begin(), test_vals.end()); + + // Check the size is correct and the ordering is not disrupted + CHECK_EQUAL(static_cast(5), rank_vector.size()); + CHECK_EQUAL(static_cast(0), rank_vector[0]); + CHECK_EQUAL(static_cast(2), rank_vector[1]); + CHECK_EQUAL(static_cast(1), rank_vector[2]); + CHECK_EQUAL(static_cast(4), rank_vector[3]); + CHECK_EQUAL(static_cast(3), rank_vector[4]); +} + +#endif // BOOST_MATH_EXEC_COMPATIBLE + +int main(void) +{ + test(); + test(); + test(); + + container_test(); + container_test(); + container_test(); + + #ifdef BOOST_MATH_EXEC_COMPATIBLE + + execution_test(std::execution::par); + execution_test(std::execution::par); + execution_test(std::execution::par); + + #endif // BOOST_MATH_EXEC_COMPATIBLE + + return boost::math::test::report_errors(); +} diff --git a/test/test_rayleigh.cpp b/test/test_rayleigh.cpp index 388adc68fc..6aa2d2fb83 100644 --- a/test/test_rayleigh.cpp +++ b/test/test_rayleigh.cpp @@ -26,6 +26,8 @@ using std::cout; using std::endl; using std::setprecision; +#include + using std::log; template void test_spot(RealType s, RealType x, RealType p, RealType q, RealType tolerance) @@ -159,6 +161,25 @@ void test_spots(RealType T) static_cast(exp_minus_half() /2), // probability. tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + rayleigh_distribution(1.L), + static_cast(1.L)), // x + log(static_cast(exp_minus_half())), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + rayleigh_distribution(0.5L), + static_cast(0.5L)), // x + log(static_cast(2 * exp_minus_half())), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( + ::boost::math::logpdf( + rayleigh_distribution(2.L), + static_cast(2.L)), // x + log(static_cast(exp_minus_half() /2)), // probability. + tolerance); // % + BOOST_CHECK_CLOSE( ::boost::math::mean( rayleigh_distribution(1.L)), @@ -255,6 +276,27 @@ BOOST_AUTO_TEST_CASE( test_main ) static_cast(exp_minus_half() /2 ), // p 1e-15); // % + BOOST_CHECK_CLOSE_FRACTION( + ::boost::math::logpdf( + rayleigh_distribution(1.), + static_cast(1)), // x + log(static_cast(exp_minus_half())), // p + 1e-15); // % + + BOOST_CHECK_CLOSE_FRACTION( + ::boost::math::logpdf( + rayleigh_distribution(0.5), + static_cast(0.5)), // x + log(static_cast(2 * exp_minus_half())), // p + 1e-15); // % + + BOOST_CHECK_CLOSE_FRACTION( + ::boost::math::logpdf( + rayleigh_distribution(2.), + static_cast(2)), // x + log(static_cast(exp_minus_half() /2 )), // p + 1e-15); // % + BOOST_CHECK_CLOSE_FRACTION( ::boost::math::cdf( rayleigh_distribution(1.), diff --git a/test/test_weibull.cpp b/test/test_weibull.cpp index 17e2bffb04..6cacf5402a 100644 --- a/test/test_weibull.cpp +++ b/test/test_weibull.cpp @@ -29,6 +29,8 @@ using std::setprecision; #include using std::numeric_limits; +#include + using std::log; template void check_weibull(RealType shape, RealType scale, RealType x, RealType p, RealType q, RealType tol) @@ -244,6 +246,50 @@ void test_spots(RealType) static_cast(0.551819), tolerance); + // + // Tests for logpdf + // + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.25, 0.5), static_cast(0.1)), + log(static_cast(0.856579)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.25, 0.5), static_cast(0.5)), + log(static_cast(0.183940)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.25, 0.5), static_cast(5)), + log(static_cast(0.015020)), + tolerance * 10); // fewer digits in test value + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.5, 2), static_cast(0.1)), + log(static_cast(0.894013)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.5, 2), static_cast(0.5)), + log(static_cast(0.303265)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(0.5, 2), static_cast(1)), + log(static_cast(0.174326)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(2, 0.25), static_cast(0.1)), + log(static_cast(2.726860)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(2, 0.25), static_cast(0.5)), + log(static_cast(0.293050)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(3, 2), static_cast(1)), + log(static_cast(0.330936)), + tolerance); + BOOST_CHECK_CLOSE( + logpdf(weibull_distribution(3, 2), static_cast(2)), + log(static_cast(0.551819)), + tolerance); + // // These test values were obtained using the formulas at // http://en.wikipedia.org/wiki/Weibull_distribution