math/statistics/chatterjee_correlation.hpp

0001 //  (C) Copyright Matt Borland 2022.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_STATISTICS_CHATTERJEE_CORRELATION_HPP
0007 #define BOOST_MATH_STATISTICS_CHATTERJEE_CORRELATION_HPP
0008
0009 #include <cstdint>
0010 #include <cmath>
0011 #include <algorithm>
0012 #include <iterator>
0013 #include <vector>
0014 #include <limits>
0015 #include <utility>
0016 #include <type_traits>
0017 #include <boost/math/tools/assert.hpp>
0018 #include <boost/math/tools/config.hpp>
0019 #include <boost/math/statistics/detail/rank.hpp>
0020
0021 #ifdef BOOST_MATH_EXEC_COMPATIBLE
0022 #include <execution>
0023 #include <future>
0024 #include <thread>
0025 #endif
0026
0027 namespace boost { namespace math { namespace statistics {
0028
0029 namespace detail {
0030
0031 template <typename BDIter>
0032 std::size_t chatterjee_transform(BDIter begin, BDIter end)
0033 {
0034     std::size_t sum = 0;
0035
0036     while(++begin != end)
0037     {
0038         if(*begin > *std::prev(begin))
0039         {
0040             sum += *begin - *std::prev(begin);
0041         }
0042         else
0043         {
0044             sum += *std::prev(begin) - *begin;
0045         }
0046     }
0047
0048     return sum;
0049 }
0050
0051 template <typename ReturnType, typename ForwardIterator>
0052 ReturnType chatterjee_correlation_seq_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)
0053 {
0054     using std::abs;
0055
0056     BOOST_MATH_ASSERT_MSG(std::is_sorted(u_begin, u_end), "The x values must be sorted in order to use this functionality");
0057
0058     const std::vector<std::size_t> rank_vector = rank(v_begin, v_end);
0059
0060     std::size_t sum = chatterjee_transform(rank_vector.begin(), rank_vector.end());
0061
0062     ReturnType result = static_cast<ReturnType>(1) - (static_cast<ReturnType>(3 * sum) / static_cast<ReturnType>(rank_vector.size() * rank_vector.size() - 1));
0063
0064     // If the result is 1 then Y is constant and all the elements must be ties
0065     if (abs(result - static_cast<ReturnType>(1)) < std::numeric_limits<ReturnType>::epsilon())
0066     {
0067         return std::numeric_limits<ReturnType>::quiet_NaN();
0068     }
0069
0070     return result;
0071 }
0072
0073 } // Namespace detail
0074
0075 template <typename Container, typename Real = typename Container::value_type,
0076           typename ReturnType = typename std::conditional<std::is_integral<Real>::value, double, Real>::type>
0077 inline ReturnType chatterjee_correlation(const Container& u, const Container& v)
0078 {
0079     return detail::chatterjee_correlation_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));
0080 }
0081
0082 }}} // Namespace boost::math::statistics
0083
0084 #ifdef BOOST_MATH_EXEC_COMPATIBLE
0085
0086 namespace boost::math::statistics {
0087
0088 namespace detail {
0089
0090 template <typename ReturnType, typename ExecutionPolicy, typename ForwardIterator>
0091 ReturnType chatterjee_correlation_par_impl(ExecutionPolicy&& exec, ForwardIterator u_begin, ForwardIterator u_end,
0092                                                                    ForwardIterator v_begin, ForwardIterator v_end)
0093 {
0094     using std::abs;
0095     BOOST_MATH_ASSERT_MSG(std::is_sorted(std::forward<ExecutionPolicy>(exec), u_begin, u_end), "The x values must be sorted in order to use this functionality");
0096
0097     auto rank_vector = rank(std::forward<ExecutionPolicy>(exec), v_begin, v_end);
0098
0099     const auto num_threads = std::thread::hardware_concurrency() == 0 ? 2u : std::thread::hardware_concurrency();
0100     std::vector<std::future<std::size_t>> future_manager {};
0101     const auto elements_per_thread = std::ceil(static_cast<double>(rank_vector.size()) / num_threads);
0102
0103     auto it = rank_vector.begin();
0104     auto end = rank_vector.end();
0105     for(std::size_t i {}; i < num_threads - 1; ++i)
0106     {
0107         future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [it, elements_per_thread]() -> std::size_t
0108         {
0109             return chatterjee_transform(it, std::next(it, elements_per_thread));
0110         }));
0111         it = std::next(it, elements_per_thread - 1);
0112     }
0113
0114     future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [it, end]() -> std::size_t
0115     {
0116         return chatterjee_transform(it, end);
0117     }));
0118
0119     std::size_t sum {};
0120     for(std::size_t i {}; i < future_manager.size(); ++i)
0121     {
0122         sum += future_manager[i].get();
0123     }
0124
0125     ReturnType result = static_cast<ReturnType>(1) - (static_cast<ReturnType>(3 * sum) / static_cast<ReturnType>(rank_vector.size() * rank_vector.size() - 1));
0126
0127     // If the result is 1 then Y is constant and all the elements must be ties
0128     if (abs(result - static_cast<ReturnType>(1)) < std::numeric_limits<ReturnType>::epsilon())
0129     {
0130         return std::numeric_limits<ReturnType>::quiet_NaN();
0131     }
0132
0133     return result;
0134 }
0135
0136 } // Namespace detail
0137
0138 template <typename ExecutionPolicy, typename Container, typename Real = typename Container::value_type,
0139           typename ReturnType = std::conditional_t<std::is_integral_v<Real>, double, Real>>
0140 inline ReturnType chatterjee_correlation(ExecutionPolicy&& exec, const Container& u, const Container& v)
0141 {
0142     if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)
0143     {
0144         return detail::chatterjee_correlation_seq_impl<ReturnType>(std::cbegin(u), std::cend(u),
0145                                                                    std::cbegin(v), std::cend(v));
0146     }
0147     else
0148     {
0149         return detail::chatterjee_correlation_par_impl<ReturnType>(std::forward<ExecutionPolicy>(exec),
0150                                                                    std::cbegin(u), std::cend(u),
0151                                                                    std::cbegin(v), std::cend(v));
0152     }
0153 }
0154
0155 } // Namespace boost::math::statistics
0156
0157 #endif
0158
0159 #endif // BOOST_MATH_STATISTICS_CHATTERJEE_CORRELATION_HPP