ceph/src/boost/boost/math/statistics/bivariate_statistics.hpp

   1 //  (C) Copyright Nick Thompson 2018.
   2 //  (C) Copyright Matt Borland 2021.
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_STATISTICS_BIVARIATE_STATISTICS_HPP
   8 #define BOOST_MATH_STATISTICS_BIVARIATE_STATISTICS_HPP
   9
  10 #include <iterator>
  11 #include <tuple>
  12 #include <type_traits>
  13 #include <stdexcept>
  14 #include <vector>
  15 #include <algorithm>
  16 #include <cmath>
  17 #include <cstddef>
  18 #include <boost/math/tools/assert.hpp>
  19 #include <boost/math/tools/config.hpp>
  20
  21 // Support compilers with P0024R2 implemented without linking TBB
  22 // https://en.cppreference.com/w/cpp/compiler_support
  23 #if !defined(BOOST_NO_CXX17_HDR_EXECUTION) && defined(BOOST_HAS_THREADS)
  24 #include <execution>
  25 #include <future>
  26 #include <thread>
  27 #define EXEC_COMPATIBLE
  28 #endif
  29
  30 namespace boost{ namespace math{ namespace statistics { namespace detail {
  31
  32 // See Equation III.9 of "Numerically Stable, Single-Pass, Parallel Statistics Algorithms", Bennet et al.
  33 template<typename ReturnType, typename ForwardIterator>
  34 ReturnType means_and_covariance_seq_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)
  35 {
  36     using Real = typename std::tuple_element<0, ReturnType>::type;
  37
  38     Real cov = 0;
  39     ForwardIterator u_it = u_begin;
  40     ForwardIterator v_it = v_begin;
  41     Real mu_u = *u_it++;
  42     Real mu_v = *v_it++;
  43     std::size_t i = 1;
  44
  45     while(u_it != u_end && v_it != v_end)
  46     {
  47         Real u_temp = (*u_it++ - mu_u)/(i+1);
  48         Real v_temp = *v_it++ - mu_v;
  49         cov += i*u_temp*v_temp;
  50         mu_u = mu_u + u_temp;
  51         mu_v = mu_v + v_temp/(i+1);
  52         i = i + 1;
  53     }
  54
  55     if(u_it != u_end || v_it != v_end)
  56     {
  57         throw std::domain_error("The size of each sample set must be the same to compute covariance");
  58     }
  59
  60     return std::make_tuple(mu_u, mu_v, cov/i, Real(i));
  61 }
  62
  63 #ifdef EXEC_COMPATIBLE
  64
  65 // Numerically stable parallel computation of (co-)variance
  66 // https://dl.acm.org/doi/10.1145/3221269.3223036
  67 template<typename ReturnType, typename ForwardIterator>
  68 ReturnType means_and_covariance_parallel_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)
  69 {
  70     using Real = typename std::tuple_element<0, ReturnType>::type;
  71
  72     const auto u_elements = std::distance(u_begin, u_end);
  73     const auto v_elements = std::distance(v_begin, v_end);
  74
  75     if(u_elements != v_elements)
  76     {
  77         throw std::domain_error("The size of each sample set must be the same to compute covariance");
  78     }
  79
  80     const unsigned max_concurrency = std::thread::hardware_concurrency() == 0 ? 2u : std::thread::hardware_concurrency();
  81     unsigned num_threads = 2u;
  82
  83     // 5.16 comes from benchmarking. See boost/math/reporting/performance/bivariate_statistics_performance.cpp
  84     // Threading is faster for: 10 + 5.16e-3 N/j <= 5.16e-3N => N >= 10^4j/5.16(j-1).
  85     const auto parallel_lower_bound = 10e4*max_concurrency/(5.16*(max_concurrency-1));
  86     const auto parallel_upper_bound = 10e4*2/5.16; // j = 2
  87
  88     // https://lemire.me/blog/2020/01/30/cost-of-a-thread-in-c-under-linux/
  89     if(u_elements < parallel_lower_bound)
  90     {
  91         return means_and_covariance_seq_impl<ReturnType>(u_begin, u_end, v_begin, v_end);
  92     }
  93     else if(u_elements >= parallel_upper_bound)
  94     {
  95         num_threads = max_concurrency;
  96     }
  97     else
  98     {
  99         for(unsigned i = 3; i < max_concurrency; ++i)
 100         {
 101             if(parallel_lower_bound < 10e4*i/(5.16*(i-1)))
 102             {
 103                 num_threads = i;
 104                 break;
 105             }
 106         }
 107     }
 108
 109     std::vector<std::future<ReturnType>> future_manager;
 110     const auto elements_per_thread = std::ceil(static_cast<double>(u_elements)/num_threads);
 111
 112     ForwardIterator u_it = u_begin;
 113     ForwardIterator v_it = v_begin;
 114
 115     for(std::size_t i = 0; i < num_threads - 1; ++i)
 116     {
 117         future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [u_it, v_it, elements_per_thread]() -> ReturnType
 118         {
 119             return means_and_covariance_seq_impl<ReturnType>(u_it, std::next(u_it, elements_per_thread), v_it, std::next(v_it, elements_per_thread));
 120         }));
 121         u_it = std::next(u_it, elements_per_thread);
 122         v_it = std::next(v_it, elements_per_thread);
 123     }
 124
 125     future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [u_it, u_end, v_it, v_end]() -> ReturnType
 126     {
 127         return means_and_covariance_seq_impl<ReturnType>(u_it, u_end, v_it, v_end);
 128     }));
 129
 130     ReturnType temp = future_manager[0].get();
 131     Real mu_u_a = std::get<0>(temp);
 132     Real mu_v_a = std::get<1>(temp);
 133     Real cov_a = std::get<2>(temp);
 134     Real n_a = std::get<3>(temp);
 135
 136     for(std::size_t i = 1; i < future_manager.size(); ++i)
 137     {
 138         temp = future_manager[i].get();
 139         Real mu_u_b = std::get<0>(temp);
 140         Real mu_v_b = std::get<1>(temp);
 141         Real cov_b = std::get<2>(temp);
 142         Real n_b = std::get<3>(temp);
 143
 144         const Real n_ab = n_a + n_b;
 145         const Real delta_u = mu_u_b - mu_u_a;
 146         const Real delta_v = mu_v_b - mu_v_a;
 147
 148         cov_a = cov_a + cov_b + (-delta_u)*(-delta_v)*((n_a*n_b)/n_ab);
 149         mu_u_a = mu_u_a + delta_u*(n_b/n_ab);
 150         mu_v_a = mu_v_a + delta_v*(n_b/n_ab);
 151         n_a = n_ab;
 152     }
 153
 154     return std::make_tuple(mu_u_a, mu_v_a, cov_a, n_a);
 155 }
 156
 157 #endif // EXEC_COMPATIBLE
 158
 159 template<typename ReturnType, typename ForwardIterator>
 160 ReturnType correlation_coefficient_seq_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)
 161 {
 162     using Real = typename std::tuple_element<0, ReturnType>::type;
 163     using std::sqrt;
 164
 165     Real cov = 0;
 166     ForwardIterator u_it = u_begin;
 167     ForwardIterator v_it = v_begin;
 168     Real mu_u = *u_it++;
 169     Real mu_v = *v_it++;
 170     Real Qu = 0;
 171     Real Qv = 0;
 172     std::size_t i = 1;
 173
 174     while(u_it != u_end && v_it != v_end)
 175     {
 176         Real u_tmp = *u_it++ - mu_u;
 177         Real v_tmp = *v_it++ - mu_v;
 178         Qu = Qu + (i*u_tmp*u_tmp)/(i+1);
 179         Qv = Qv + (i*v_tmp*v_tmp)/(i+1);
 180         cov += i*u_tmp*v_tmp/(i+1);
 181         mu_u = mu_u + u_tmp/(i+1);
 182         mu_v = mu_v + v_tmp/(i+1);
 183         ++i;
 184     }
 185
 186
 187     // If one dataset is constant, then the correlation coefficient is undefined.
 188     // See https://stats.stackexchange.com/questions/23676/normalized-correlation-with-a-constant-vector
 189     // Thanks to zbjornson for pointing this out.
 190     if (Qu == 0 || Qv == 0)
 191     {
 192         return std::make_tuple(mu_u, Qu, mu_v, Qv, cov, std::numeric_limits<Real>::quiet_NaN(), Real(i));
 193     }
 194
 195     // Make sure rho in [-1, 1], even in the presence of numerical noise.
 196     Real rho = cov/sqrt(Qu*Qv);
 197     if (rho > 1) {
 198         rho = 1;
 199     }
 200     if (rho < -1) {
 201         rho = -1;
 202     }
 203
 204     return std::make_tuple(mu_u, Qu, mu_v, Qv, cov, rho, Real(i));
 205 }
 206
 207 #ifdef EXEC_COMPATIBLE
 208
 209 // Numerically stable parallel computation of (co-)variance:
 210 // https://dl.acm.org/doi/10.1145/3221269.3223036
 211 //
 212 // Parallel computation of variance:
 213 // http://i.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf
 214 template<typename ReturnType, typename ForwardIterator>
 215 ReturnType correlation_coefficient_parallel_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)
 216 {
 217     using Real = typename std::tuple_element<0, ReturnType>::type;
 218
 219     const auto u_elements = std::distance(u_begin, u_end);
 220     const auto v_elements = std::distance(v_begin, v_end);
 221
 222     if(u_elements != v_elements)
 223     {
 224         throw std::domain_error("The size of each sample set must be the same to compute covariance");
 225     }
 226
 227     const unsigned max_concurrency = std::thread::hardware_concurrency() == 0 ? 2u : std::thread::hardware_concurrency();
 228     unsigned num_threads = 2u;
 229
 230     // 3.25 comes from benchmarking. See boost/math/reporting/performance/bivariate_statistics_performance.cpp
 231     // Threading is faster for: 10 + 3.25e-3 N/j <= 3.25e-3N => N >= 10^4j/3.25(j-1).
 232     const auto parallel_lower_bound = 10e4*max_concurrency/(3.25*(max_concurrency-1));
 233     const auto parallel_upper_bound = 10e4*2/3.25; // j = 2
 234
 235     // https://lemire.me/blog/2020/01/30/cost-of-a-thread-in-c-under-linux/
 236     if(u_elements < parallel_lower_bound)
 237     {
 238         return correlation_coefficient_seq_impl<ReturnType>(u_begin, u_end, v_begin, v_end);
 239     }
 240     else if(u_elements >= parallel_upper_bound)
 241     {
 242         num_threads = max_concurrency;
 243     }
 244     else
 245     {
 246         for(unsigned i = 3; i < max_concurrency; ++i)
 247         {
 248             if(parallel_lower_bound < 10e4*i/(3.25*(i-1)))
 249             {
 250                 num_threads = i;
 251                 break;
 252             }
 253         }
 254     }
 255
 256     std::vector<std::future<ReturnType>> future_manager;
 257     const auto elements_per_thread = std::ceil(static_cast<double>(u_elements)/num_threads);
 258
 259     ForwardIterator u_it = u_begin;
 260     ForwardIterator v_it = v_begin;
 261
 262     for(std::size_t i = 0; i < num_threads - 1; ++i)
 263     {
 264         future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [u_it, v_it, elements_per_thread]() -> ReturnType
 265         {
 266             return correlation_coefficient_seq_impl<ReturnType>(u_it, std::next(u_it, elements_per_thread), v_it, std::next(v_it, elements_per_thread));
 267         }));
 268         u_it = std::next(u_it, elements_per_thread);
 269         v_it = std::next(v_it, elements_per_thread);
 270     }
 271
 272     future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [u_it, u_end, v_it, v_end]() -> ReturnType
 273     {
 274         return correlation_coefficient_seq_impl<ReturnType>(u_it, u_end, v_it, v_end);
 275     }));
 276
 277     ReturnType temp = future_manager[0].get();
 278     Real mu_u_a = std::get<0>(temp);
 279     Real Qu_a = std::get<1>(temp);
 280     Real mu_v_a = std::get<2>(temp);
 281     Real Qv_a = std::get<3>(temp);
 282     Real cov_a = std::get<4>(temp);
 283     Real n_a = std::get<6>(temp);
 284
 285     for(std::size_t i = 1; i < future_manager.size(); ++i)
 286     {
 287         temp = future_manager[i].get();
 288         Real mu_u_b = std::get<0>(temp);
 289         Real Qu_b = std::get<1>(temp);
 290         Real mu_v_b = std::get<2>(temp);
 291         Real Qv_b = std::get<3>(temp);
 292         Real cov_b = std::get<4>(temp);
 293         Real n_b = std::get<6>(temp);
 294
 295         const Real n_ab = n_a + n_b;
 296         const Real delta_u = mu_u_b - mu_u_a;
 297         const Real delta_v = mu_v_b - mu_v_a;
 298
 299         cov_a = cov_a + cov_b + (-delta_u)*(-delta_v)*((n_a*n_b)/n_ab);
 300         mu_u_a = mu_u_a + delta_u*(n_b/n_ab);
 301         mu_v_a = mu_v_a + delta_v*(n_b/n_ab);
 302         Qu_a = Qu_a + Qu_b + delta_u*delta_u*((n_a*n_b)/n_ab);
 303         Qv_b = Qv_a + Qv_b + delta_v*delta_v*((n_a*n_b)/n_ab);
 304         n_a = n_ab;
 305     }
 306
 307     // If one dataset is constant, then the correlation coefficient is undefined.
 308     // See https://stats.stackexchange.com/questions/23676/normalized-correlation-with-a-constant-vector
 309     // Thanks to zbjornson for pointing this out.
 310     if (Qu_a == 0 || Qv_a == 0)
 311     {
 312         return std::make_tuple(mu_u_a, Qu_a, mu_v_a, Qv_a, cov_a, std::numeric_limits<Real>::quiet_NaN(), n_a);
 313     }
 314
 315     // Make sure rho in [-1, 1], even in the presence of numerical noise.
 316     Real rho = cov_a/sqrt(Qu_a*Qv_a);
 317     if (rho > 1) {
 318         rho = 1;
 319     }
 320     if (rho < -1) {
 321         rho = -1;
 322     }
 323
 324     return std::make_tuple(mu_u_a, Qu_a, mu_v_a, Qv_a, cov_a, rho, n_a);
 325 }
 326
 327 #endif // EXEC_COMPATIBLE
 328
 329 } // namespace detail
 330
 331 #ifdef EXEC_COMPATIBLE
 332
 333 template<typename ExecutionPolicy, typename Container, typename Real = typename Container::value_type>
 334 inline auto means_and_covariance(ExecutionPolicy&& exec, Container const & u, Container const & v)
 335 {
 336     if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)
 337     {
 338         if constexpr (std::is_integral_v<Real>)
 339         {
 340             using ReturnType = std::tuple<double, double, double, double>;
 341             ReturnType temp = detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 342             return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));
 343         }
 344         else
 345         {
 346             using ReturnType = std::tuple<Real, Real, Real, Real>;
 347             ReturnType temp = detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 348             return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));
 349         }
 350     }
 351     else
 352     {
 353         if constexpr (std::is_integral_v<Real>)
 354         {
 355             using ReturnType = std::tuple<double, double, double, double>;
 356             ReturnType temp = detail::means_and_covariance_parallel_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 357             return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));
 358         }
 359         else
 360         {
 361             using ReturnType = std::tuple<Real, Real, Real, Real>;
 362             ReturnType temp = detail::means_and_covariance_parallel_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 363             return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));
 364         }
 365     }
 366 }
 367
 368 template<typename Container>
 369 inline auto means_and_covariance(Container const & u, Container const & v)
 370 {
 371     return means_and_covariance(std::execution::seq, u, v);
 372 }
 373
 374 template<typename ExecutionPolicy, typename Container>
 375 inline auto covariance(ExecutionPolicy&& exec, Container const & u, Container const & v)
 376 {
 377     return std::get<2>(means_and_covariance(exec, u, v));
 378 }
 379
 380 template<typename Container>
 381 inline auto covariance(Container const & u, Container const & v)
 382 {
 383     return covariance(std::execution::seq, u, v);
 384 }
 385
 386 template<typename ExecutionPolicy, typename Container, typename Real = typename Container::value_type>
 387 inline auto correlation_coefficient(ExecutionPolicy&& exec, Container const & u, Container const & v)
 388 {
 389     if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)
 390     {
 391         if constexpr (std::is_integral_v<Real>)
 392         {
 393             using ReturnType = std::tuple<double, double, double, double, double, double, double>;
 394             return std::get<5>(detail::correlation_coefficient_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));
 395         }
 396         else
 397         {
 398             using ReturnType = std::tuple<Real, Real, Real, Real, Real, Real, Real>;
 399             return std::get<5>(detail::correlation_coefficient_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));
 400         }
 401     }
 402     else
 403     {
 404         if constexpr (std::is_integral_v<Real>)
 405         {
 406             using ReturnType = std::tuple<double, double, double, double, double, double, double>;
 407             return std::get<5>(detail::correlation_coefficient_parallel_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));
 408         }
 409         else
 410         {
 411             using ReturnType = std::tuple<Real, Real, Real, Real, Real, Real, Real>;
 412             return std::get<5>(detail::correlation_coefficient_parallel_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));
 413         }
 414     }
 415 }
 416
 417 template<typename Container, typename Real = typename Container::value_type>
 418 inline auto correlation_coefficient(Container const & u, Container const & v)
 419 {
 420     return correlation_coefficient(std::execution::seq, u, v);
 421 }
 422
 423 #else // C++11 and single threaded bindings
 424
 425 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 426 inline auto means_and_covariance(Container const & u, Container const & v) -> std::tuple<double, double, double>
 427 {
 428     using ReturnType = std::tuple<double, double, double, double>;
 429     ReturnType temp = detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 430     return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));
 431 }
 432
 433 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 434 inline auto means_and_covariance(Container const & u, Container const & v) -> std::tuple<Real, Real, Real>
 435 {
 436     using ReturnType = std::tuple<Real, Real, Real, Real>;
 437     ReturnType temp = detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 438     return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));
 439 }
 440
 441 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 442 inline double covariance(Container const & u, Container const & v)
 443 {
 444     using ReturnType = std::tuple<double, double, double, double>;
 445     return std::get<2>(detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));
 446 }
 447
 448 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 449 inline Real covariance(Container const & u, Container const & v)
 450 {
 451     using ReturnType = std::tuple<Real, Real, Real, Real>;
 452     return std::get<2>(detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));
 453 }
 454
 455 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 456 inline double correlation_coefficient(Container const & u, Container const & v)
 457 {
 458     using ReturnType = std::tuple<double, double, double, double, double, double, double>;
 459     return std::get<5>(detail::correlation_coefficient_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));
 460 }
 461
 462 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 463 inline Real correlation_coefficient(Container const & u, Container const & v)
 464 {
 465     using ReturnType = std::tuple<Real, Real, Real, Real, Real, Real, Real>;
 466     return std::get<5>(detail::correlation_coefficient_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));
 467 }
 468
 469 #endif
 470
 471 }}} // namespace boost::math::statistics
 472
 473 #endif