1 /* 2 * Created by Martin on 15/06/2019. 3 * Adapted from donated nonius code. 4 * 5 * Distributed under the Boost Software License, Version 1.0. (See accompanying 6 * file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 7 */ 8 9 // Statistical analysis tools 10 11 #if defined(CATCH_CONFIG_ENABLE_BENCHMARKING) 12 13 #include "catch_stats.hpp" 14 15 #include "../../catch_compiler_capabilities.h" 16 17 #include <cassert> 18 #include <random> 19 20 21 #if defined(CATCH_CONFIG_USE_ASYNC) 22 #include <future> 23 #endif 24 25 namespace { erf_inv(double x)26 double erf_inv(double x) { 27 // Code accompanying the article "Approximating the erfinv function" in GPU Computing Gems, Volume 2 28 double w, p; 29 30 w = -log((1.0 - x) * (1.0 + x)); 31 32 if (w < 6.250000) { 33 w = w - 3.125000; 34 p = -3.6444120640178196996e-21; 35 p = -1.685059138182016589e-19 + p * w; 36 p = 1.2858480715256400167e-18 + p * w; 37 p = 1.115787767802518096e-17 + p * w; 38 p = -1.333171662854620906e-16 + p * w; 39 p = 2.0972767875968561637e-17 + p * w; 40 p = 6.6376381343583238325e-15 + p * w; 41 p = -4.0545662729752068639e-14 + p * w; 42 p = -8.1519341976054721522e-14 + p * w; 43 p = 2.6335093153082322977e-12 + p * w; 44 p = -1.2975133253453532498e-11 + p * w; 45 p = -5.4154120542946279317e-11 + p * w; 46 p = 1.051212273321532285e-09 + p * w; 47 p = -4.1126339803469836976e-09 + p * w; 48 p = -2.9070369957882005086e-08 + p * w; 49 p = 4.2347877827932403518e-07 + p * w; 50 p = -1.3654692000834678645e-06 + p * w; 51 p = -1.3882523362786468719e-05 + p * w; 52 p = 0.0001867342080340571352 + p * w; 53 p = -0.00074070253416626697512 + p * w; 54 p = -0.0060336708714301490533 + p * w; 55 p = 0.24015818242558961693 + p * w; 56 p = 1.6536545626831027356 + p * w; 57 } else if (w < 16.000000) { 58 w = sqrt(w) - 3.250000; 59 p = 2.2137376921775787049e-09; 60 p = 9.0756561938885390979e-08 + p * w; 61 p = -2.7517406297064545428e-07 + p * w; 62 p = 1.8239629214389227755e-08 + p * w; 63 p = 1.5027403968909827627e-06 + p * w; 64 p = -4.013867526981545969e-06 + p * w; 65 p = 2.9234449089955446044e-06 + p * w; 66 p = 1.2475304481671778723e-05 + p * w; 67 p = -4.7318229009055733981e-05 + p * w; 68 p = 6.8284851459573175448e-05 + p * w; 69 p = 2.4031110387097893999e-05 + p * w; 70 p = -0.0003550375203628474796 + p * w; 71 p = 0.00095328937973738049703 + p * w; 72 p = -0.0016882755560235047313 + p * w; 73 p = 0.0024914420961078508066 + p * w; 74 p = -0.0037512085075692412107 + p * w; 75 p = 0.005370914553590063617 + p * w; 76 p = 1.0052589676941592334 + p * w; 77 p = 3.0838856104922207635 + p * w; 78 } else { 79 w = sqrt(w) - 5.000000; 80 p = -2.7109920616438573243e-11; 81 p = -2.5556418169965252055e-10 + p * w; 82 p = 1.5076572693500548083e-09 + p * w; 83 p = -3.7894654401267369937e-09 + p * w; 84 p = 7.6157012080783393804e-09 + p * w; 85 p = -1.4960026627149240478e-08 + p * w; 86 p = 2.9147953450901080826e-08 + p * w; 87 p = -6.7711997758452339498e-08 + p * w; 88 p = 2.2900482228026654717e-07 + p * w; 89 p = -9.9298272942317002539e-07 + p * w; 90 p = 4.5260625972231537039e-06 + p * w; 91 p = -1.9681778105531670567e-05 + p * w; 92 p = 7.5995277030017761139e-05 + p * w; 93 p = -0.00021503011930044477347 + p * w; 94 p = -0.00013871931833623122026 + p * w; 95 p = 1.0103004648645343977 + p * w; 96 p = 4.8499064014085844221 + p * w; 97 } 98 return p * x; 99 } 100 standard_deviation(std::vector<double>::iterator first,std::vector<double>::iterator last)101 double standard_deviation(std::vector<double>::iterator first, std::vector<double>::iterator last) { 102 auto m = Catch::Benchmark::Detail::mean(first, last); 103 double variance = std::accumulate(first, last, 0., [m](double a, double b) { 104 double diff = b - m; 105 return a + diff * diff; 106 }) / (last - first); 107 return std::sqrt(variance); 108 } 109 110 } 111 112 namespace Catch { 113 namespace Benchmark { 114 namespace Detail { 115 weighted_average_quantile(int k,int q,std::vector<double>::iterator first,std::vector<double>::iterator last)116 double weighted_average_quantile(int k, int q, std::vector<double>::iterator first, std::vector<double>::iterator last) { 117 auto count = last - first; 118 double idx = (count - 1) * k / static_cast<double>(q); 119 int j = static_cast<int>(idx); 120 double g = idx - j; 121 std::nth_element(first, first + j, last); 122 auto xj = first[j]; 123 if (g == 0) return xj; 124 125 auto xj1 = *std::min_element(first + (j + 1), last); 126 return xj + g * (xj1 - xj); 127 } 128 129 erfc_inv(double x)130 double erfc_inv(double x) { 131 return erf_inv(1.0 - x); 132 } 133 normal_quantile(double p)134 double normal_quantile(double p) { 135 static const double ROOT_TWO = std::sqrt(2.0); 136 137 double result = 0.0; 138 assert(p >= 0 && p <= 1); 139 if (p < 0 || p > 1) { 140 return result; 141 } 142 143 result = -erfc_inv(2.0 * p); 144 // result *= normal distribution standard deviation (1.0) * sqrt(2) 145 result *= /*sd * */ ROOT_TWO; 146 // result += normal disttribution mean (0) 147 return result; 148 } 149 150 outlier_variance(Estimate<double> mean,Estimate<double> stddev,int n)151 double outlier_variance(Estimate<double> mean, Estimate<double> stddev, int n) { 152 double sb = stddev.point; 153 double mn = mean.point / n; 154 double mg_min = mn / 2.; 155 double sg = std::min(mg_min / 4., sb / std::sqrt(n)); 156 double sg2 = sg * sg; 157 double sb2 = sb * sb; 158 159 auto c_max = [n, mn, sb2, sg2](double x) -> double { 160 double k = mn - x; 161 double d = k * k; 162 double nd = n * d; 163 double k0 = -n * nd; 164 double k1 = sb2 - n * sg2 + nd; 165 double det = k1 * k1 - 4 * sg2 * k0; 166 return (int)(-2. * k0 / (k1 + std::sqrt(det))); 167 }; 168 169 auto var_out = [n, sb2, sg2](double c) { 170 double nc = n - c; 171 return (nc / n) * (sb2 - nc * sg2); 172 }; 173 174 return std::min(var_out(1), var_out(std::min(c_max(0.), c_max(mg_min)))) / sb2; 175 } 176 177 analyse_samples(double confidence_level,int n_resamples,std::vector<double>::iterator first,std::vector<double>::iterator last)178 bootstrap_analysis analyse_samples(double confidence_level, int n_resamples, std::vector<double>::iterator first, std::vector<double>::iterator last) { 179 CATCH_INTERNAL_START_WARNINGS_SUPPRESSION 180 CATCH_INTERNAL_SUPPRESS_GLOBALS_WARNINGS 181 static std::random_device entropy; 182 CATCH_INTERNAL_STOP_WARNINGS_SUPPRESSION 183 184 auto n = static_cast<int>(last - first); // seriously, one can't use integral types without hell in C++ 185 186 auto mean = &Detail::mean<std::vector<double>::iterator>; 187 auto stddev = &standard_deviation; 188 189 #if defined(CATCH_CONFIG_USE_ASYNC) 190 auto Estimate = [=](double(*f)(std::vector<double>::iterator, std::vector<double>::iterator)) { 191 auto seed = entropy(); 192 return std::async(std::launch::async, [=] { 193 std::mt19937 rng(seed); 194 auto resampled = resample(rng, n_resamples, first, last, f); 195 return bootstrap(confidence_level, first, last, resampled, f); 196 }); 197 }; 198 199 auto mean_future = Estimate(mean); 200 auto stddev_future = Estimate(stddev); 201 202 auto mean_estimate = mean_future.get(); 203 auto stddev_estimate = stddev_future.get(); 204 #else 205 auto Estimate = [=](double(*f)(std::vector<double>::iterator, std::vector<double>::iterator)) { 206 auto seed = entropy(); 207 std::mt19937 rng(seed); 208 auto resampled = resample(rng, n_resamples, first, last, f); 209 return bootstrap(confidence_level, first, last, resampled, f); 210 }; 211 212 auto mean_estimate = Estimate(mean); 213 auto stddev_estimate = Estimate(stddev); 214 #endif // CATCH_USE_ASYNC 215 216 double outlier_variance = Detail::outlier_variance(mean_estimate, stddev_estimate, n); 217 218 return { mean_estimate, stddev_estimate, outlier_variance }; 219 } 220 } // namespace Detail 221 } // namespace Benchmark 222 } // namespace Catch 223 224 #endif // CATCH_CONFIG_ENABLE_BENCHMARKING 225