libcxx/test/std/numerics/rand/rand.dist/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp

   1 //===----------------------------------------------------------------------===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8 //
   9 // REQUIRES: long_tests
  10
  11 // <random>
  12
  13 // template<class IntType = int>
  14 // class poisson_distribution
  15
  16 // template<class _URNG> result_type operator()(_URNG& g);
  17
  18 #include <cassert>
  19 #include <cmath>
  20 #include <cstdint>
  21 #include <limits>
  22 #include <numeric>
  23 #include <random>
  24 #include <vector>
  25
  26 #include "test_macros.h"
  27
  28 template <class T>
  29 T sqr(T x) {
  30   return x * x;
  31 }
  32
  33 void test_bad_ranges() {
  34   // Test cases where the mean is around the largest representable integer for
  35   // `result_type`. These cases don't generate valid poisson distributions, but
  36   // at least they don't blow up.
  37   std::mt19937 eng;
  38
  39   {
  40     std::poisson_distribution<std::int16_t> distribution(32710.9);
  41     for (int i=0; i < 1000; ++i) {
  42       volatile std::int16_t res = distribution(eng);
  43       ((void)res);
  44     }
  45   }
  46   {
  47     std::poisson_distribution<std::int16_t> distribution(std::numeric_limits<std::int16_t>::max());
  48     for (int i=0; i < 1000; ++i) {
  49       volatile std::int16_t res = distribution(eng);
  50       ((void)res);
  51     }
  52   }
  53   {
  54     std::poisson_distribution<std::int16_t> distribution(
  55     static_cast<double>(std::numeric_limits<std::int16_t>::max()) + 10);
  56     for (int i=0; i < 1000; ++i) {
  57       volatile std::int16_t res = distribution(eng);
  58       ((void)res);
  59     }
  60   }
  61   {
  62     std::poisson_distribution<std::int16_t> distribution(
  63       static_cast<double>(std::numeric_limits<std::int16_t>::max()) * 2);
  64       for (int i=0; i < 1000; ++i) {
  65         volatile std::int16_t res = distribution(eng);
  66         ((void)res);
  67       }
  68   }
  69   {
  70     // We convert `INF` to `DBL_MAX` otherwise the distribution will hang.
  71     std::poisson_distribution<std::int16_t> distribution(std::numeric_limits<double>::infinity());
  72     for (int i=0; i < 1000; ++i) {
  73       volatile std::int16_t res = distribution(eng);
  74       ((void)res);
  75     }
  76   }
  77   {
  78     std::poisson_distribution<std::int16_t> distribution(0);
  79     for (int i=0; i < 1000; ++i) {
  80       volatile std::int16_t res = distribution(eng);
  81       ((void)res);
  82     }
  83   }
  84   {
  85     // We convert `INF` to `DBL_MAX` otherwise the distribution will hang.
  86     std::poisson_distribution<std::int16_t> distribution(-100);
  87     for (int i=0; i < 1000; ++i) {
  88       volatile std::int16_t res = distribution(eng);
  89       ((void)res);
  90     }
  91   }
  92 }
  93
  94 template <class T>
  95 void tests() {
  96   {
  97     typedef std::poisson_distribution<T> D;
  98     typedef std::minstd_rand G;
  99     G g;
 100     D d(2);
 101     const int N = 100000;
 102     std::vector<double> u;
 103     for (int i = 0; i < N; ++i)
 104     {
 105       typename D::result_type v = d(g);
 106       assert(d.min() <= v && v <= d.max());
 107       u.push_back(v);
 108     }
 109     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
 110     double var = 0;
 111     double skew = 0;
 112     double kurtosis = 0;
 113     for (unsigned i = 0; i < u.size(); ++i)
 114     {
 115       double dbl = (u[i] - mean);
 116       double d2 = sqr(dbl);
 117       var += d2;
 118       skew += dbl * d2;
 119       kurtosis += d2 * d2;
 120     }
 121     var /= u.size();
 122     double dev = std::sqrt(var);
 123     skew /= u.size() * dev * var;
 124     kurtosis /= u.size() * var * var;
 125     kurtosis -= 3;
 126     double x_mean = d.mean();
 127     double x_var = d.mean();
 128     double x_skew = 1 / std::sqrt(x_var);
 129     double x_kurtosis = 1 / x_var;
 130     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 131     assert(std::abs((var - x_var) / x_var) < 0.01);
 132     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 133     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
 134   }
 135   {
 136     typedef std::poisson_distribution<T> D;
 137     typedef std::minstd_rand G;
 138     G g;
 139     D d(0.75);
 140     const int N = 100000;
 141     std::vector<double> u;
 142     for (int i = 0; i < N; ++i)
 143     {
 144       typename D::result_type v = d(g);
 145       assert(d.min() <= v && v <= d.max());
 146       u.push_back(v);
 147     }
 148     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
 149     double var = 0;
 150     double skew = 0;
 151     double kurtosis = 0;
 152     for (unsigned i = 0; i < u.size(); ++i)
 153     {
 154       double dbl = (u[i] - mean);
 155       double d2 = sqr(dbl);
 156       var += d2;
 157       skew += dbl * d2;
 158       kurtosis += d2 * d2;
 159     }
 160     var /= u.size();
 161     double dev = std::sqrt(var);
 162     skew /= u.size() * dev * var;
 163     kurtosis /= u.size() * var * var;
 164     kurtosis -= 3;
 165     double x_mean = d.mean();
 166     double x_var = d.mean();
 167     double x_skew = 1 / std::sqrt(x_var);
 168     double x_kurtosis = 1 / x_var;
 169     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 170     assert(std::abs((var - x_var) / x_var) < 0.01);
 171     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 172     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04);
 173   }
 174   {
 175     typedef std::poisson_distribution<T> D;
 176     typedef std::mt19937 G;
 177     G g;
 178     D d(20);
 179     const int N = 1000000;
 180     std::vector<double> u;
 181     for (int i = 0; i < N; ++i)
 182     {
 183       typename D::result_type v = d(g);
 184       assert(d.min() <= v && v <= d.max());
 185       u.push_back(v);
 186     }
 187     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
 188     double var = 0;
 189     double skew = 0;
 190     double kurtosis = 0;
 191     for (unsigned i = 0; i < u.size(); ++i)
 192     {
 193       double dbl = (u[i] - mean);
 194       double d2 = sqr(dbl);
 195       var += d2;
 196       skew += dbl * d2;
 197       kurtosis += d2 * d2;
 198     }
 199     var /= u.size();
 200     double dev = std::sqrt(var);
 201     skew /= u.size() * dev * var;
 202     kurtosis /= u.size() * var * var;
 203     kurtosis -= 3;
 204     double x_mean = d.mean();
 205     double x_var = d.mean();
 206     double x_skew = 1 / std::sqrt(x_var);
 207     double x_kurtosis = 1 / x_var;
 208     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 209     assert(std::abs((var - x_var) / x_var) < 0.01);
 210     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 211     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
 212   }
 213 }
 214
 215 int main(int, char**) {
 216   test_bad_ranges();
 217
 218   tests<short>();
 219   tests<int>();
 220   tests<long>();
 221   tests<long long>();
 222
 223   tests<unsigned short>();
 224   tests<unsigned int>();
 225   tests<unsigned long>();
 226   tests<unsigned long long>();
 227
 228 #if defined(_LIBCPP_VERSION) // extension
 229   // TODO: std::poisson_distribution currently doesn't work reliably with small types.
 230   // tests<int8_t>();
 231   // tests<uint8_t>();
 232 #if !defined(TEST_HAS_NO_INT128)
 233   tests<__int128_t>();
 234   tests<__uint128_t>();
 235 #endif
 236 #endif
 237
 238   return 0;
 239 }