1 /** 2 Numeric related utilities used by TSV Utilities. 3 4 Utilities in this file: 5 $(LIST 6 * [formatNumber] - An alternate print format for numbers, especially useful when 7 doubles are being used to represent integer and float values. 8 9 * [rangeMedian] - Finds the median value of a range. 10 11 * [quantile] - Generates quantile values for a data set. 12 ) 13 14 Copyright (c) 2016-2020, eBay Inc. 15 Initially written by Jon Degenhardt 16 17 License: Boost Licence 1.0 (http://boost.org/LICENSE_1_0.txt) 18 */ 19 20 module tsv_utils.common.numerics; 21 22 import std.traits : isFloatingPoint, isIntegral, Unqual; 23 24 /** 25 formatNumber is an alternate way to print numbers. It is especially useful when 26 representing both integral and floating point values with float point data types. 27 28 formatNumber was created for tsv-summarize, where all calculations are done as doubles, 29 but may be integers by nature. In addition, output may be either for human consumption 30 or for additional machine processing. Integers are printed normally. Floating point is 31 printed as follows: 32 $(LIST 33 * Values that are exact integral values are printed as integers, as long as they 34 are within the range of where all integers are represented exactly by the floating 35 point type. The practical effect is to avoid switching to exponential notion. 36 37 * If the specified floatPrecision is between 0 and readablePrecisionMax, then 38 floatPrecision is used to set the significant digits following the decimal point. 39 Otherwise, it is used to set total significant digits. This does not apply to 40 really large numbers, for doubles, those larger than 2^53. Trailing zeros are 41 chopped in all cases. 42 ) 43 */ 44 auto formatNumber(T, size_t readablePrecisionMax = 6)(T num, const size_t floatPrecision = 12) 45 if (isFloatingPoint!T || isIntegral!T) 46 { 47 alias UT = Unqual!T; 48 49 import std.conv : to; 50 import std.format : format; 51 52 static if (isIntegral!T) 53 { 54 return format("%d", num); // The easy case. 55 } 56 else 57 { 58 static assert(isFloatingPoint!T); 59 60 static if (!is(UT == float) && !is(UT == double)) 61 { 62 /* Not a double or float, but a floating point. Punt on refinements. */ 63 return format("%.*g", floatPrecision, num); 64 } 65 else 66 { 67 static assert(is(UT == float) || is(UT == double)); 68 69 if (floatPrecision <= readablePrecisionMax) 70 { 71 /* Print with a fixed precision beyond the decimal point (%.*f), but 72 * remove trailing zeros. Notes: 73 * - This handles integer values stored in floating point types. 74 * - Values like NaN and infinity also handled. 75 */ 76 immutable str = format("%.*f", floatPrecision, num); 77 size_t trimToLength = str.length; 78 79 if (floatPrecision != 0 && str.length > floatPrecision + 1) 80 { 81 import std.ascii : isDigit; 82 assert(str.length - floatPrecision - 1 > 0); 83 size_t decimalIndex = str.length - floatPrecision - 1; 84 85 if (str[decimalIndex] == '.' && str[decimalIndex - 1].isDigit) 86 { 87 size_t lastNonZeroDigit = str.length - 1; 88 assert(decimalIndex < lastNonZeroDigit); 89 while (str[lastNonZeroDigit] == '0') lastNonZeroDigit--; 90 trimToLength = (decimalIndex < lastNonZeroDigit) 91 ? lastNonZeroDigit + 1 92 : decimalIndex; 93 } 94 } 95 96 return str[0 .. trimToLength]; 97 } 98 else 99 { 100 /* Determine if the number is subject to special integer value printing. 101 * Goal is to avoid exponential notion for integer values that '%.*g' 102 * generates. Numbers within the significant digit range of floatPrecision 103 * will print as desired with '%.*g', whether there is a fractional part 104 * or not. The '%.*g' format, with exponential notation, is also used for 105 * really large numbers. "Really large" being numbers outside the range 106 * of integers exactly representable by the floating point type. 107 */ 108 109 enum UT maxConsecutiveUTInteger = 2.0^^UT.mant_dig; 110 enum bool maxUTIntFitsInLong = (maxConsecutiveUTInteger <= long.max); 111 112 import std.math : fabs; 113 immutable UT absNum = num.fabs; 114 115 if (!maxUTIntFitsInLong || 116 absNum < 10.0^^floatPrecision || 117 absNum > maxConsecutiveUTInteger) 118 { 119 /* Within signficant digits range or very large. */ 120 return format("%.*g", floatPrecision, num); 121 } 122 else 123 { 124 /* Check for integral values needing to be printed in decimal format. 125 * modf/modff are used to determine if the value has a non-zero 126 * fractional component. 127 */ 128 import core.stdc.math : modf, modff; 129 130 static if (is(UT == float)) alias modfUT = modff; 131 else static if (is(UT == double)) alias modfUT = modf; 132 else static assert(0); 133 134 UT integerPart; 135 136 if (modfUT(num, &integerPart) == 0.0) return format("%d", num.to!long); 137 else return format("%.*g", floatPrecision, num); 138 } 139 } 140 } 141 } 142 assert(0); 143 } 144 145 unittest // formatNumber unit tests 146 { 147 import std.conv : to; 148 import std.format : format; 149 150 /* Integers */ 151 assert(formatNumber(0) == "0"); 152 assert(formatNumber(1) == "1"); 153 assert(formatNumber(-1) == "-1"); 154 assert(formatNumber(999) == "999"); 155 assert(formatNumber(12345678912345) == "12345678912345"); 156 assert(formatNumber(-12345678912345) == "-12345678912345"); 157 158 size_t a1 = 10; assert(a1.formatNumber == "10"); 159 const int a2 = -33234; assert(a2.formatNumber == "-33234"); 160 immutable long a3 = -12345678912345; assert(a3.formatNumber == "-12345678912345"); 161 162 // Specifying precision should never matter for integer values. 163 assert(formatNumber(0, 0) == "0"); 164 assert(formatNumber(1, 0) == "1"); 165 assert(formatNumber(-1, 0) == "-1"); 166 assert(formatNumber(999, 0) == "999"); 167 assert(formatNumber(12345678912345, 0) == "12345678912345"); 168 assert(formatNumber(-12345678912345, 0) == "-12345678912345"); 169 170 assert(formatNumber(0, 3) == "0"); 171 assert(formatNumber(1, 3) == "1"); 172 assert(formatNumber(-1, 3 ) == "-1"); 173 assert(formatNumber(999, 3) == "999"); 174 assert(formatNumber(12345678912345, 3) == "12345678912345"); 175 assert(formatNumber(-12345678912345, 3) == "-12345678912345"); 176 177 assert(formatNumber(0, 9) == "0"); 178 assert(formatNumber(1, 9) == "1"); 179 assert(formatNumber(-1, 9 ) == "-1"); 180 assert(formatNumber(999, 9) == "999"); 181 assert(formatNumber(12345678912345, 9) == "12345678912345"); 182 assert(formatNumber(-12345678912345, 9) == "-12345678912345"); 183 184 /* Doubles */ 185 assert(formatNumber(0.0) == "0"); 186 assert(formatNumber(0.2) == "0.2"); 187 assert(formatNumber(0.123412, 0) == "0"); 188 assert(formatNumber(0.123412, 1) == "0.1"); 189 assert(formatNumber(0.123412, 2) == "0.12"); 190 assert(formatNumber(0.123412, 5) == "0.12341"); 191 assert(formatNumber(0.123412, 6) == "0.123412"); 192 assert(formatNumber(0.123412, 7) == "0.123412"); 193 assert(formatNumber(9.123412, 5) == "9.12341"); 194 assert(formatNumber(9.123412, 6) == "9.123412"); 195 assert(formatNumber(99.123412, 5) == "99.12341"); 196 assert(formatNumber(99.123412, 6) == "99.123412"); 197 assert(formatNumber(99.123412, 7) == "99.12341"); 198 assert(formatNumber(999.123412, 0) == "999"); 199 assert(formatNumber(999.123412, 1) == "999.1"); 200 assert(formatNumber(999.123412, 2) == "999.12"); 201 assert(formatNumber(999.123412, 3) == "999.123"); 202 assert(formatNumber(999.123412, 4) == "999.1234"); 203 assert(formatNumber(999.123412, 5) == "999.12341"); 204 assert(formatNumber(999.123412, 6) == "999.123412"); 205 assert(formatNumber(999.123412, 7) == "999.1234"); 206 assert(formatNumber!(double, 9)(999.12341234, 7) == "999.1234123"); 207 assert(formatNumber(9001.0) == "9001"); 208 assert(formatNumber(1234567891234.0) == "1234567891234"); 209 assert(formatNumber(1234567891234.0, 0) == "1234567891234"); 210 assert(formatNumber(1234567891234.0, 1) == "1234567891234"); 211 212 // Test round off cases 213 assert(formatNumber(0.6, 0) == "1"); 214 assert(formatNumber(0.6, 1) == "0.6"); 215 assert(formatNumber(0.06, 0) == "0"); 216 assert(formatNumber(0.06, 1) == "0.1"); 217 assert(formatNumber(0.06, 2) == "0.06"); 218 assert(formatNumber(0.06, 3) == "0.06"); 219 assert(formatNumber(9.49999, 0) == "9"); 220 assert(formatNumber(9.49999, 1) == "9.5"); 221 assert(formatNumber(9.6, 0) == "10"); 222 assert(formatNumber(9.6, 1) == "9.6"); 223 assert(formatNumber(99.99, 0) == "100"); 224 assert(formatNumber(99.99, 1) == "100"); 225 assert(formatNumber(99.99, 2) == "99.99"); 226 assert(formatNumber(9999.9996, 3) == "10000"); 227 assert(formatNumber(9999.9996, 4) == "9999.9996"); 228 assert(formatNumber(99999.99996, 4) == "100000"); 229 assert(formatNumber(99999.99996, 5) == "99999.99996"); 230 assert(formatNumber(999999.999996, 5) == "1000000"); 231 assert(formatNumber(999999.999996, 6) == "999999.999996"); 232 233 /* Turn off precision, the 'human readable' style. 234 * Note: Remains o if both are zero (first test). If it becomes desirable to support 235 * turning it off when for the precision equal zero case the simple extension is to 236 * allow the 'human readable' precision template parameter to be negative. 237 */ 238 assert(formatNumber!(double, 0)(999.123412, 0) == "999"); 239 assert(formatNumber!(double, 0)(999.123412, 1) == "1e+03"); 240 assert(formatNumber!(double, 0)(999.123412, 2) == "1e+03"); 241 assert(formatNumber!(double, 0)(999.123412, 3) == "999"); 242 assert(formatNumber!(double, 0)(999.123412, 4) == "999.1"); 243 244 // Default number printing 245 assert(formatNumber(1.2) == "1.2"); 246 assert(formatNumber(12.3) == "12.3"); 247 assert(formatNumber(12.34) == "12.34"); 248 assert(formatNumber(123.45) == "123.45"); 249 assert(formatNumber(123.456) == "123.456"); 250 assert(formatNumber(1234.567) == "1234.567"); 251 assert(formatNumber(1234.5678) == "1234.5678"); 252 assert(formatNumber(12345.6789) == "12345.6789"); 253 assert(formatNumber(12345.67891) == "12345.67891"); 254 assert(formatNumber(123456.78912) == "123456.78912"); 255 assert(formatNumber(123456.789123) == "123456.789123"); 256 assert(formatNumber(1234567.891234) == "1234567.89123"); 257 assert(formatNumber(12345678.912345) == "12345678.9123"); 258 assert(formatNumber(123456789.12345) == "123456789.123"); 259 assert(formatNumber(1234567891.2345) == "1234567891.23"); 260 assert(formatNumber(12345678912.345) == "12345678912.3"); 261 assert(formatNumber(123456789123.45) == "123456789123"); 262 assert(formatNumber(1234567891234.5) == "1.23456789123e+12"); 263 assert(formatNumber(12345678912345.6) == "1.23456789123e+13"); 264 assert(formatNumber(123456789123456.0) == "123456789123456"); 265 assert(formatNumber(0.3) == "0.3"); 266 assert(formatNumber(0.03) == "0.03"); 267 assert(formatNumber(0.003) == "0.003"); 268 assert(formatNumber(0.0003) == "0.0003"); 269 assert(formatNumber(0.00003) == "3e-05" || formatNumber(0.00003) == "3e-5"); 270 assert(formatNumber(0.000003) == "3e-06" || formatNumber(0.000003) == "3e-6"); 271 assert(formatNumber(0.0000003) == "3e-07" || formatNumber(0.0000003) == "3e-7"); 272 273 // Large number inside and outside the contiguous integer representation range 274 double dlarge = 2.0^^(double.mant_dig - 2) - 10.0; 275 double dhuge = 2.0^^(double.mant_dig + 1) + 1000.0; 276 277 assert(dlarge.formatNumber == format("%d", dlarge.to!long)); 278 assert(dhuge.formatNumber!(double) == format("%.12g", dhuge)); 279 280 // Negative values - Repeat most of above tests. 281 assert(formatNumber(-0.0) == "-0" || formatNumber(-0.0) == "0"); 282 assert(formatNumber(-0.2) == "-0.2"); 283 assert(formatNumber(-0.123412, 0) == "-0"); 284 assert(formatNumber(-0.123412, 1) == "-0.1"); 285 assert(formatNumber(-0.123412, 2) == "-0.12"); 286 assert(formatNumber(-0.123412, 5) == "-0.12341"); 287 assert(formatNumber(-0.123412, 6) == "-0.123412"); 288 assert(formatNumber(-0.123412, 7) == "-0.123412"); 289 assert(formatNumber(-9.123412, 5) == "-9.12341"); 290 assert(formatNumber(-9.123412, 6) == "-9.123412"); 291 assert(formatNumber(-99.123412, 5) == "-99.12341"); 292 assert(formatNumber(-99.123412, 6) == "-99.123412"); 293 assert(formatNumber(-99.123412, 7) == "-99.12341"); 294 assert(formatNumber(-999.123412, 0) == "-999"); 295 assert(formatNumber(-999.123412, 1) == "-999.1"); 296 assert(formatNumber(-999.123412, 2) == "-999.12"); 297 assert(formatNumber(-999.123412, 3) == "-999.123"); 298 assert(formatNumber(-999.123412, 4) == "-999.1234"); 299 assert(formatNumber(-999.123412, 5) == "-999.12341"); 300 assert(formatNumber(-999.123412, 6) == "-999.123412"); 301 assert(formatNumber(-999.123412, 7) == "-999.1234"); 302 assert(formatNumber!(double, 9)(-999.12341234, 7) == "-999.1234123"); 303 assert(formatNumber(-9001.0) == "-9001"); 304 assert(formatNumber(-1234567891234.0) == "-1234567891234"); 305 assert(formatNumber(-1234567891234.0, 0) == "-1234567891234"); 306 assert(formatNumber(-1234567891234.0, 1) == "-1234567891234"); 307 308 // Test round off cases 309 assert(formatNumber(-0.6, 0) == "-1"); 310 assert(formatNumber(-0.6, 1) == "-0.6"); 311 assert(formatNumber(-0.06, 0) == "-0"); 312 assert(formatNumber(-0.06, 1) == "-0.1"); 313 assert(formatNumber(-0.06, 2) == "-0.06"); 314 assert(formatNumber(-0.06, 3) == "-0.06"); 315 assert(formatNumber(-9.49999, 0) == "-9"); 316 assert(formatNumber(-9.49999, 1) == "-9.5"); 317 assert(formatNumber(-9.6, 0) == "-10"); 318 assert(formatNumber(-9.6, 1) == "-9.6"); 319 assert(formatNumber(-99.99, 0) == "-100"); 320 assert(formatNumber(-99.99, 1) == "-100"); 321 assert(formatNumber(-99.99, 2) == "-99.99"); 322 assert(formatNumber(-9999.9996, 3) == "-10000"); 323 assert(formatNumber(-9999.9996, 4) == "-9999.9996"); 324 assert(formatNumber(-99999.99996, 4) == "-100000"); 325 assert(formatNumber(-99999.99996, 5) == "-99999.99996"); 326 assert(formatNumber(-999999.999996, 5) == "-1000000"); 327 assert(formatNumber(-999999.999996, 6) == "-999999.999996"); 328 329 assert(formatNumber!(double, 0)(-999.123412, 0) == "-999"); 330 assert(formatNumber!(double, 0)(-999.123412, 1) == "-1e+03"); 331 assert(formatNumber!(double, 0)(-999.123412, 2) == "-1e+03"); 332 assert(formatNumber!(double, 0)(-999.123412, 3) == "-999"); 333 assert(formatNumber!(double, 0)(-999.123412, 4) == "-999.1"); 334 335 // Default number printing 336 assert(formatNumber(-1.2) == "-1.2"); 337 assert(formatNumber(-12.3) == "-12.3"); 338 assert(formatNumber(-12.34) == "-12.34"); 339 assert(formatNumber(-123.45) == "-123.45"); 340 assert(formatNumber(-123.456) == "-123.456"); 341 assert(formatNumber(-1234.567) == "-1234.567"); 342 assert(formatNumber(-1234.5678) == "-1234.5678"); 343 assert(formatNumber(-12345.6789) == "-12345.6789"); 344 assert(formatNumber(-12345.67891) == "-12345.67891"); 345 assert(formatNumber(-123456.78912) == "-123456.78912"); 346 assert(formatNumber(-123456.789123) == "-123456.789123"); 347 assert(formatNumber(-1234567.891234) == "-1234567.89123"); 348 349 assert(formatNumber(-12345678.912345) == "-12345678.9123"); 350 assert(formatNumber(-123456789.12345) == "-123456789.123"); 351 assert(formatNumber(-1234567891.2345) == "-1234567891.23"); 352 assert(formatNumber(-12345678912.345) == "-12345678912.3"); 353 assert(formatNumber(-123456789123.45) == "-123456789123"); 354 assert(formatNumber(-1234567891234.5) == "-1.23456789123e+12"); 355 assert(formatNumber(-12345678912345.6) == "-1.23456789123e+13"); 356 assert(formatNumber(-123456789123456.0) == "-123456789123456"); 357 358 assert(formatNumber(-0.3) == "-0.3"); 359 assert(formatNumber(-0.03) == "-0.03"); 360 assert(formatNumber(-0.003) == "-0.003"); 361 assert(formatNumber(-0.0003) == "-0.0003"); 362 assert(formatNumber(-0.00003) == "-3e-05" || formatNumber(-0.00003) == "-3e-5"); 363 assert(formatNumber(-0.000003) == "-3e-06" || formatNumber(-0.000003) == "-3e-6"); 364 assert(formatNumber(-0.0000003) == "-3e-07" || formatNumber(-0.0000003) == "-3e-7"); 365 366 const double dlargeNeg = -2.0^^(double.mant_dig - 2) + 10.0; 367 immutable double dhugeNeg = -2.0^^(double.mant_dig + 1) - 1000.0; 368 369 assert(dlargeNeg.formatNumber == format("%d", dlargeNeg.to!long)); 370 assert(dhugeNeg.formatNumber!(double) == format("%.12g", dhugeNeg)); 371 372 // Type qualifiers 373 const double b1 = 0.0; assert(formatNumber(b1) == "0"); 374 const double b2 = 0.2; assert(formatNumber(b2) == "0.2"); 375 const double b3 = 0.123412; assert(formatNumber(b3, 2) == "0.12"); 376 immutable double b4 = 99.123412; assert(formatNumber(b4, 5) == "99.12341"); 377 immutable double b5 = 99.123412; assert(formatNumber(b5, 7) == "99.12341"); 378 379 // Special values 380 assert(formatNumber(double.nan) == "nan"); 381 assert(formatNumber(double.nan, 0) == "nan"); 382 assert(formatNumber(double.nan, 1) == "nan"); 383 assert(formatNumber(double.nan, 9) == "nan"); 384 assert(formatNumber(double.infinity) == "inf"); 385 assert(formatNumber(double.infinity, 0) == "inf"); 386 assert(formatNumber(double.infinity, 1) == "inf"); 387 assert(formatNumber(double.infinity, 9) == "inf"); 388 389 // Float. Mix negative and type qualifiers in. 390 assert(formatNumber(0.0f) == "0"); 391 assert(formatNumber(0.5f) == "0.5"); 392 assert(formatNumber(0.123412f, 0) == "0"); 393 assert(formatNumber(0.123412f, 1) == "0.1"); 394 assert(formatNumber(-0.123412f, 2) == "-0.12"); 395 assert(formatNumber(9.123412f, 5) == "9.12341"); 396 assert(formatNumber(9.123412f, 6) == "9.123412"); 397 assert(formatNumber(-99.123412f, 5) == "-99.12341"); 398 assert(formatNumber(99.123412f, 7) == "99.12341"); 399 assert(formatNumber(-999.123412f, 5) == "-999.12341"); 400 401 float c1 = 999.123412f; assert(formatNumber(c1, 7) == "999.1234"); 402 float c2 = 999.1234f; assert(formatNumber!(float, 9)(c2, 3) == "999.123"); 403 const float c3 = 9001.0f; assert(formatNumber(c3) == "9001"); 404 const float c4 = -12345678.0f; assert(formatNumber(c4) == "-12345678"); 405 immutable float c5 = 12345678.0f; assert(formatNumber(c5, 0) == "12345678"); 406 immutable float c6 = 12345678.0f; assert(formatNumber(c6, 1) == "12345678"); 407 408 float flarge = 2.0^^(float.mant_dig - 2) - 10.0; 409 float fhuge = 2.0^^(float.mant_dig + 1) + 1000.0; 410 411 assert(flarge.formatNumber == format("%d", flarge.to!long)); 412 assert(fhuge.formatNumber!(float, 12) == format("%.12g", fhuge)); 413 414 // Reals - No special formatting 415 real d1 = 2.0^^(double.mant_dig) - 1000.0; assert(formatNumber(d1) == format("%.12g", d1)); 416 real d2 = 123456789.12341234L; assert(formatNumber(d2, 12) == format("%.12g", d2)); 417 } 418 419 /** 420 rangeMedian. Finds the median. Modifies the range via topN or sort in the process. 421 422 Note: topN is the preferred algorithm, but the version prior to Phobos 2.073 423 is pathologically slow on certain data sets. Use topN in 2.073 and later, 424 sort in earlier versions. 425 426 See: https://issues.dlang.org/show_bug.cgi?id=16517 427 https://github.com/dlang/phobos/pull/4815 428 http://forum.dlang.org/post/ujuugklmbibuheptdwcn@forum.dlang.org 429 */ 430 static if (__VERSION__ >= 2073) 431 { 432 version = rangeMedianViaTopN; 433 } 434 else 435 { 436 version = rangeMedianViaSort; 437 } 438 439 auto rangeMedian (Range) (Range r) 440 if (isRandomAccessRange!Range && hasLength!Range && hasSlicing!Range) 441 { 442 version(rangeMedianViaSort) 443 { 444 version(rangeMedianViaTopN) 445 { 446 assert(0, "Both rangeMedianViaSort and rangeMedianViaTopN assigned as versions. Assign only one."); 447 } 448 } 449 else version(rangeMedianViaTopN) 450 { 451 } 452 else 453 { 454 static assert(0, "A version of rangeMedianViaSort or rangeMedianViaTopN must be assigned."); 455 } 456 457 import std.traits : isFloatingPoint; 458 459 ElementType!Range median; 460 461 if (r.length > 0) 462 { 463 size_t medianIndex = r.length / 2; 464 465 version(rangeMedianViaSort) 466 { 467 import std.algorithm : sort; 468 sort(r); 469 median = r[medianIndex]; 470 471 static if (isFloatingPoint!(ElementType!Range)) 472 { 473 if (r.length % 2 == 0) 474 { 475 /* Even number of values. Split the difference. */ 476 median = (median + r[medianIndex - 1]) / 2.0; 477 } 478 } 479 } 480 else version(rangeMedianViaTopN) 481 { 482 import std.algorithm : maxElement, topN; 483 topN(r, medianIndex); 484 median = r[medianIndex]; 485 486 static if (isFloatingPoint!(ElementType!Range)) 487 { 488 if (r.length % 2 == 0) 489 { 490 /* Even number of values. Split the difference. */ 491 if (r[medianIndex - 1] < median) 492 { 493 median = (median + r[0..medianIndex].maxElement) / 2.0; 494 } 495 } 496 } 497 } 498 else 499 { 500 static assert(0, "A version of rangeMedianViaSort or rangeMedianViaTopN must be assigned."); 501 } 502 } 503 504 return median; 505 } 506 507 /* rangeMedian unit tests. */ 508 @safe unittest 509 { 510 import std.math : isNaN; 511 import std.algorithm : all, permutations; 512 513 // Median of empty range is (type).init. Zero for int, nan for floats/doubles 514 assert(rangeMedian(new int[0]) == int.init); 515 assert(rangeMedian(new double[0]).isNaN && double.init.isNaN); 516 assert(rangeMedian(new string[0]) == ""); 517 518 assert(rangeMedian([3]) == 3); 519 assert(rangeMedian([3.0]) == 3.0); 520 assert(rangeMedian([3.5]) == 3.5); 521 assert(rangeMedian(["aaa"]) == "aaa"); 522 523 /* Even number of elements: Split the difference for floating point, but not other types. */ 524 assert(rangeMedian([3, 4]) == 4); 525 assert(rangeMedian([3.0, 4.0]) == 3.5); 526 527 assert(rangeMedian([3, 6, 12]) == 6); 528 assert(rangeMedian([3.0, 6.5, 12.5]) == 6.5); 529 530 // Do the rest with permutations 531 assert([4, 7].permutations.all!(x => (x.rangeMedian == 7))); 532 assert([4.0, 7.0].permutations.all!(x => (x.rangeMedian == 5.5))); 533 assert(["aaa", "bbb"].permutations.all!(x => (x.rangeMedian == "bbb"))); 534 535 assert([4, 7, 19].permutations.all!(x => (x.rangeMedian == 7))); 536 assert([4.5, 7.5, 19.5].permutations.all!(x => (x.rangeMedian == 7.5))); 537 assert(["aaa", "bbb", "ccc"].permutations.all!(x => (x.rangeMedian == "bbb"))); 538 539 assert([4.5, 7.5, 19.5, 21.0].permutations.all!(x => (x.rangeMedian == 13.5))); 540 assert([4.5, 7.5, 19.5, 20.5, 36.0].permutations.all!(x => (x.rangeMedian == 19.5))); 541 assert([4.5, 7.5, 19.5, 24.0, 24.5, 25.0].permutations.all!(x => (x.rangeMedian == 21.75))); 542 assert([1.5, 3.25, 3.55, 4.5, 24.5, 25.0, 25.6].permutations.all!(x => (x.rangeMedian == 4.5))); 543 } 544 545 /// Quantiles 546 547 /** The different quantile interpolation methods. 548 * See: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/quantile.html 549 */ 550 enum QuantileInterpolation 551 { 552 R1 = 1, /// R quantile type 1 553 R2 = 2, /// R quantile type 2 554 R3 = 3, /// R quantile type 3 555 R4 = 4, /// R quantile type 4 556 R5 = 5, /// R quantile type 5 557 R6 = 6, /// R quantile type 6 558 R7 = 7, /// R quantile type 7 559 R8 = 8, /// R quantile type 8 560 R9 = 9, /// R quantile type 9 561 } 562 563 564 import std.traits : isFloatingPoint, isNumeric, Unqual; 565 import std.range; 566 567 /** 568 Returns the quantile in a data vector for a cumulative probability. 569 570 Takes a data vector and a probability and returns the quantile cut point for the 571 probability. The vector must be sorted and the probability in the range [0.0, 1.0]. 572 The interpolation methods available are the same as in R and available in a number 573 of statistical packages. See the R documentation or wikipedia for details 574 (https://en.wikipedia.org/wiki/Quantile). 575 576 Examples: 577 ---- 578 double data = [22, 57, 73, 97, 113]; 579 double median = quantile(0.5, data); // 73 580 auto q1 = [0.25, 0.5, 0.75].map!(p => p.quantile(data)); // 57, 73, 97 581 auto q2 = [0.25, 0.5, 0.75].map!(p => p.quantile(data), QuantileInterpolation.R8); //45.3333, 73, 102.333 582 ---- 583 */ 584 double quantile(ProbType, Range) 585 (const ProbType prob, Range data, QuantileInterpolation method = QuantileInterpolation.R7) 586 if (isRandomAccessRange!Range && hasLength!Range && isNumeric!(ElementType!Range) && 587 isFloatingPoint!ProbType) 588 in 589 { 590 import std.algorithm : isSorted; 591 assert(0.0 <= prob && prob <= 1.0); 592 assert(method >= QuantileInterpolation.min && method <= QuantileInterpolation.max); 593 assert(data.isSorted); 594 } 595 do 596 { 597 import core.stdc.math : modf; 598 import std.algorithm : max, min; 599 import std.conv : to; 600 import std.math : ceil, lrint; 601 602 /* Note: In the implementation below, 'h1' is the 1-based index into the data vector. 603 * This follows the wikipedia notation for the interpolation methods. One will be 604 * subtracted before the vector is accessed. 605 */ 606 607 double q = double.nan; // The return value. 608 609 if (data.length == 1) q = data[0].to!double; 610 else if (data.length > 1) 611 { 612 if (method == QuantileInterpolation.R1) 613 { 614 q = data[((data.length * prob).ceil - 1.0).to!long.max(0)].to!double; 615 } 616 else if (method == QuantileInterpolation.R2) 617 { 618 immutable double h1 = data.length * prob + 0.5; 619 immutable size_t lo = ((h1 - 0.5).ceil.to!long - 1).max(0); 620 immutable size_t hi = ((h1 + 0.5).to!size_t - 1).min(data.length - 1); 621 q = (data[lo].to!double + data[hi].to!double) / 2.0; 622 } 623 else if (method == QuantileInterpolation.R3) 624 { 625 /* Implementation notes: 626 * - R3 uses 'banker's rounding', where 0.5 is rounded to the nearest even 627 * value. The 'lrint' routine does this. 628 * - DMD will sometimes choose the incorrect 0.5 rounding if the calculation 629 * is done as a single step. The separate calculation of 'h1' avoids this. 630 */ 631 immutable double h1 = data.length * prob; 632 q = data[h1.lrint.max(1) - 1].to!double; 633 } 634 else if ((method == QuantileInterpolation.R4) || 635 (method == QuantileInterpolation.R5) || 636 (method == QuantileInterpolation.R6) || 637 (method == QuantileInterpolation.R7) || 638 (method == QuantileInterpolation.R8) || 639 (method == QuantileInterpolation.R9)) 640 { 641 /* Methods 4-9 have different formulas for generating the real-valued index, 642 * but work the same after that, choosing the final value by linear interpolation. 643 */ 644 double h1; 645 switch (method) 646 { 647 case QuantileInterpolation.R4: h1 = data.length * prob; break; 648 case QuantileInterpolation.R5: h1 = data.length * prob + 0.5; break; 649 case QuantileInterpolation.R6: h1 = (data.length + 1) * prob; break; 650 case QuantileInterpolation.R7: h1 = (data.length - 1) * prob + 1.0; break; 651 case QuantileInterpolation.R8: h1 = (data.length.to!double + 1.0/3.0) * prob + 1.0/3.0; break; 652 case QuantileInterpolation.R9: h1 = (data.length + 0.25) * prob + 3.0/8.0; break; 653 default: assert(0); 654 } 655 656 double h1IntegerPart; 657 immutable double h1FractionPart = modf(h1, &h1IntegerPart); 658 immutable size_t lo = (h1IntegerPart - 1.0).to!long.max(0).min(data.length - 1); 659 q = data[lo]; 660 if (h1FractionPart > 0.0) 661 { 662 immutable size_t hi = h1IntegerPart.to!long.min(data.length - 1); 663 q += h1FractionPart * (data[hi].to!double - data[lo].to!double); 664 } 665 } 666 else assert(0); 667 } 668 return q; 669 } 670 671 unittest 672 { 673 import std.algorithm : equal, map; 674 import std.array : array; 675 import std.traits : EnumMembers; 676 677 /* A couple simple tests. */ 678 assert(quantile(0.5, [22, 57, 73, 97, 113]) == 73); 679 assert(quantile(0.5, [22.5, 57.5, 73.5, 97.5, 113.5]) == 73.5); 680 assert([0.25, 0.5, 0.75].map!(p => p.quantile([22, 57, 73, 97, 113])).array == [57.0, 73.0, 97.0]); 681 assert([0.25, 0.5, 0.75].map!(p => p.quantile([22, 57, 73, 97, 113], QuantileInterpolation.R1)).array == [57.0, 73.0, 97.0]); 682 683 /* Data arrays. */ 684 double[] d1 = []; 685 double[] d2 = [5.5]; 686 double[] d3 = [0.0, 1.0]; 687 double[] d4 = [-1.0, 1.0]; 688 double[] d5 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]; 689 double[] d6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0]; 690 double[] d7 = [ 31.79, 64.19, 81.77]; 691 double[] d8 = [-94.43, -74.55, -50.81, 27.45, 78.79]; 692 double[] d9 = [-89.17, 20.93, 38.51, 48.03, 76.43, 77.02]; 693 double[] d10 = [-99.53, -76.87, -76.69, -67.81, -40.26, -11.29, 21.02]; 694 double[] d11 = [-78.32, -52.22, -50.86, 13.45, 15.96, 17.25, 46.35, 85.00]; 695 double[] d12 = [-81.36, -70.87, -53.56, -42.14, -9.18, 7.23, 49.52, 80.43, 98.50]; 696 double[] d13 = [ 38.37, 44.36, 45.70, 50.69, 51.36, 55.66, 56.91, 58.95, 62.01, 65.25]; 697 698 /* Spot check a few other data types. Same expected outputs.*/ 699 int[] d3Int = [0, 1]; 700 int[] d4Int = [-1, 1]; 701 int[] d5Int = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; 702 size_t[] d6Size_t = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]; 703 float[] d7Float = [ 31.79f, 64.19f, 81.77f]; 704 float[] d8Float = [-94.43f, -74.55f, -50.81f, 27.45f, 78.79f]; 705 float[] d9Float = [-89.17f, 20.93f, 38.51f, 48.03f, 76.43f, 77.02f]; 706 float[] d10Float = [-99.53f, -76.87f, -76.69f, -67.81f, -40.26f, -11.29f, 21.02f]; 707 708 /* Probability values. */ 709 double[] probs = [0.0, 0.05, 0.1, 0.25, 0.4, 0.49, 0.5, 0.51, 0.75, 0.9, 0.95, 0.98, 1.0]; 710 711 /* Expected values for each data array, for 'probs'. One expected result for each of the nine methods. 712 * The expected values were generated by R and Octave. 713 */ 714 double[13][9] d1_expected; // All values double.nan, the default. 715 double[13][9] d2_expected = [ 716 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 717 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 718 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 719 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 720 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 721 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 722 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 723 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 724 [5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5], 725 ]; 726 double[13][9] d3_expected = [ 727 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 728 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 729 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0], 730 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.02, 0.5, 0.8, 0.9, 0.96, 1.0], 731 [0.0, 0.0, 0.0, 0.0, 0.3, 0.48, 0.5, 0.52, 1.0, 1.0, 1.0, 1.0, 1.0], 732 [0.0, 0.0, 0.0, 0.0, 0.2, 0.47, 0.5, 0.53, 1.0, 1.0, 1.0, 1.0, 1.0], 733 [0.0, 0.05, 0.1, 0.25, 0.4, 0.49, 0.5, 0.51, 0.75, 0.9, 0.95, 0.98, 1.0], 734 [0.0, 0.0, 0.0, 0.0, 0.2666667, 0.4766667, 0.5, 0.5233333, 1.0, 1.0, 1.0, 1.0, 1.0], 735 [0.0, 0.0, 0.0, 0.0, 0.275, 0.4775, 0.5, 0.5225, 1.0, 1.0, 1.0, 1.0, 1.0], 736 ]; 737 double[13][9] d4_expected = [ 738 [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 739 [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 740 [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 741 [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -0.96, 0.0, 0.6, 0.8, 0.92, 1.0], 742 [-1.0, -1.0, -1.0, -1.0, -0.4, -0.04, 0.0, 0.04, 1.0, 1.0, 1.0, 1.0, 1.0], 743 [-1.0, -1.0, -1.0, -1.0, -0.6, -0.06, 0.0, 0.06, 1.0, 1.0, 1.0, 1.0, 1.0], 744 [-1.0, -0.9, -0.8, -0.5, -0.2, -0.02, 0.0, 0.02, 0.5, 0.8, 0.9, 0.96, 1.0], 745 [-1.0, -1.0, -1.0, -1.0, -0.4666667, -0.04666667, -4.440892e-16, 0.04666667, 1.0, 1.0, 1.0, 1.0, 1.0], 746 [-1.0, -1.0, -1.0, -1.0, -0.45, -0.045, 0.0, 0.045, 1.0, 1.0, 1.0, 1.0, 1.0], 747 ]; 748 double[13][9] d5_expected = [ 749 [0.0, 0.0, 1.0, 2.0, 4.0, 5.0, 5.0, 5.0, 8.0, 9.0, 10.0, 10.0, 10.0], 750 [0.0, 0.0, 1.0, 2.0, 4.0, 5.0, 5.0, 5.0, 8.0, 9.0, 10.0, 10.0, 10.0], 751 [0.0, 0.0, 0.0, 2.0, 3.0, 4.0, 5.0, 5.0, 7.0, 9.0, 9.0, 10.0, 10.0], 752 [0.0, 0.0, 0.1, 1.75, 3.4, 4.39, 4.5, 4.61, 7.25, 8.9, 9.45, 9.78, 10.0], 753 [0.0, 0.05, 0.6, 2.25, 3.9, 4.89, 5.0, 5.11, 7.75, 9.4, 9.95, 10.0, 10.0], 754 [0.0, 0.0, 0.2, 2.0, 3.8, 4.88, 5.0, 5.12, 8.0, 9.8, 10.0, 10.0, 10.0], 755 [0.0, 0.5, 1.0, 2.5, 4.0, 4.9, 5.0, 5.1, 7.5, 9.0, 9.5, 9.8, 10.0], 756 [0.0, 0.0, 0.4666667, 2.166667, 3.866667, 4.886667, 5.0, 5.113333, 7.833333, 9.533333, 10.0, 10.0, 10.0], 757 [0.0, 0.0, 0.5, 2.1875, 3.875, 4.8875, 5.0, 5.1125, 7.8125, 9.5, 10.0, 10.0, 10.0], 758 ]; 759 double[13][9] d6_expected = [ 760 [0.0, 0.0, 1.0, 2.0, 4.0, 5.0, 5.0, 6.0, 8.0, 10.0, 11.0, 11.0, 11.0], 761 [0.0, 0.0, 1.0, 2.5, 4.0, 5.0, 5.5, 6.0, 8.5, 10.0, 11.0, 11.0, 11.0], 762 [0.0, 0.0, 0.0, 2.0, 4.0, 5.0, 5.0, 5.0, 8.0, 10.0, 10.0, 11.0, 11.0], 763 [0.0, 0.0, 0.2, 2.0, 3.8, 4.88, 5.0, 5.12, 8.0, 9.8, 10.4, 10.76, 11.0], 764 [0.0, 0.1, 0.7, 2.5, 4.3, 5.38, 5.5, 5.62, 8.5, 10.3, 10.9, 11.0, 11.0], 765 [0.0, 0.0, 0.3, 2.25, 4.2, 5.37, 5.5, 5.63, 8.75, 10.7, 11.0, 11.0, 11.0], 766 [0.0, 0.55, 1.1, 2.75, 4.4, 5.39, 5.5, 5.61, 8.25, 9.9, 10.45, 10.78, 11.0], 767 [0.0, 0.0, 0.5666667, 2.416667, 4.266667, 5.376667, 5.5, 5.623333, 8.583333, 10.43333, 11.0, 11.0, 11.0], 768 [0.0, 0.0, 0.6, 2.4375, 4.275, 5.3775, 5.5, 5.6225, 8.5625, 10.4, 11.0, 11.0, 11.0], 769 ]; 770 double[13][9] d7_expected = [ 771 [31.79, 31.79, 31.79, 31.79, 64.19, 64.19, 64.19, 64.19, 81.77, 81.77, 81.77, 81.77, 81.77], 772 [31.79, 31.79, 31.79, 31.79, 64.19, 64.19, 64.19, 64.19, 81.77, 81.77, 81.77, 81.77, 81.77], 773 [31.79, 31.79, 31.79, 31.79, 31.79, 31.79, 64.19, 64.19, 64.19, 81.77, 81.77, 81.77, 81.77], 774 [31.79, 31.79, 31.79, 31.79, 38.27, 47.018, 47.99, 48.962, 68.585, 76.496, 79.133, 80.7152, 81.77], 775 [31.79, 31.79, 31.79, 39.89, 54.47, 63.218, 64.19, 64.7174, 77.375, 81.77, 81.77, 81.77, 81.77], 776 [31.79, 31.79, 31.79, 31.79, 51.23, 62.894, 64.19, 64.8932, 81.77, 81.77, 81.77, 81.77, 81.77], 777 [31.79, 35.03, 38.27, 47.99, 57.71, 63.542, 64.19, 64.5416, 72.98, 78.254, 80.012, 81.0668, 81.77], 778 [31.79, 31.79, 31.79, 37.19, 53.39, 63.11, 64.19, 64.776, 78.84, 81.77, 81.77, 81.77, 81.77], 779 [31.79, 31.79, 31.79, 37.865, 53.66, 63.137, 64.19, 64.76135, 78.47375, 81.77, 81.77, 81.77, 81.77], 780 ]; 781 double[13][9] d8_expected = [ 782 [-94.43, -94.43, -94.43, -74.55, -74.55, -50.81, -50.81, -50.81, 27.45, 78.79, 78.79, 78.79, 78.79], 783 [-94.43, -94.43, -94.43, -74.55, -62.68, -50.81, -50.81, -50.81, 27.45, 78.79, 78.79, 78.79, 78.79], 784 [-94.43, -94.43, -94.43, -94.43, -74.55, -74.55, -74.55, -50.81, 27.45, 27.45, 78.79, 78.79, 78.79], 785 [-94.43, -94.43, -94.43, -89.46, -74.55, -63.867, -62.68, -61.493, 7.885, 53.12, 65.955, 73.656, 78.79], 786 [-94.43, -94.43, -94.43, -79.52, -62.68, -51.997, -50.81, -46.897, 40.285, 78.79, 78.79, 78.79, 78.79], 787 [-94.43, -94.43, -94.43, -84.49, -65.054, -52.2344, -50.81, -46.1144, 53.12, 78.79, 78.79, 78.79, 78.79], 788 [-94.43, -90.454, -86.478, -74.55, -60.306, -51.7596, -50.81, -47.6796, 27.45, 58.254, 68.522, 74.6828, 78.79], 789 [-94.43, -94.43, -94.43, -81.17667, -63.47133, -52.07613, -50.81, -46.63613, 44.56333, 78.79, 78.79, 78.79, 78.79], 790 [-94.43, -94.43, -94.43, -80.7625, -63.2735, -52.05635, -50.81, -46.70135, 43.49375, 78.79, 78.79, 78.79, 78.79], 791 ]; 792 double[13][9] d9_expected = [ 793 [-89.17, -89.17, -89.17, 20.93, 38.51, 38.51, 38.51, 48.03, 76.43, 77.02, 77.02, 77.02, 77.02], 794 [-89.17, -89.17, -89.17, 20.93, 38.51, 38.51, 43.27, 48.03, 76.43, 77.02, 77.02, 77.02, 77.02], 795 [-89.17, -89.17, -89.17, 20.93, 20.93, 38.51, 38.51, 38.51, 48.03, 76.43, 77.02, 77.02, 77.02], 796 [-89.17, -89.17, -89.17, -34.12, 27.962, 37.4552, 38.51, 39.0812, 62.23, 76.666, 76.843, 76.9492, 77.02], 797 [-89.17, -89.17, -78.16, 20.93, 36.752, 42.6988, 43.27, 43.8412, 76.43, 76.961, 77.02, 77.02, 77.02], 798 [-89.17, -89.17, -89.17, -6.595, 34.994, 42.6036, 43.27, 43.9364, 76.5775, 77.02, 77.02, 77.02, 77.02], 799 [-89.17, -61.645, -34.12, 25.325, 38.51, 42.794, 43.27, 43.746, 69.33, 76.725, 76.8725, 76.961, 77.02], 800 [-89.17, -89.17, -89.17, 11.755, 36.166, 42.66707, 43.27, 43.87293, 76.47917, 77.02, 77.02, 77.02, 77.02], 801 [-89.17, -89.17, -89.17, 14.04875, 36.3125, 42.675, 43.27, 43.865, 76.46688, 77.02, 77.02, 77.02, 77.02], 802 ]; 803 double[13][9] d10_expected = [ 804 [-99.53, -99.53, -99.53, -76.87, -76.69, -67.81, -67.81, -67.81, -11.29, 21.02, 21.02, 21.02, 21.02], 805 [-99.53, -99.53, -99.53, -76.87, -76.69, -67.81, -67.81, -67.81, -11.29, 21.02, 21.02, 21.02, 21.02], 806 [-99.53, -99.53, -99.53, -76.87, -76.69, -76.69, -67.81, -67.81, -40.26, -11.29, 21.02, 21.02, 21.02], 807 [-99.53, -99.53, -99.53, -82.535, -76.726, -72.8716, -72.25, -71.6284, -33.0175, -1.597, 9.7115, 16.4966, 21.02], 808 [-99.53, -99.53, -94.998, -76.825, -74.026, -68.4316, -67.81, -65.8815, -18.5325, 14.558, 21.02, 21.02, 21.02], 809 [-99.53, -99.53, -99.53, -76.87, -74.914, -68.5204, -67.81, -65.606, -11.29, 21.02, 21.02, 21.02, 21.02], 810 [-99.53, -92.732, -85.934, -76.78, -73.138, -68.3428, -67.81, -66.157, -25.775, 1.634, 11.327, 17.1428, 21.02], 811 [-99.53, -99.53, -98.01933, -76.84, -74.322, -68.4612, -67.81, -65.78967, -16.11833, 18.866, 21.02, 21.02, 21.02], 812 [-99.53, -99.53, -97.264, -76.83625, -74.248, -68.4538, -67.81, -65.81263, -16.72187, 17.789, 21.02, 21.02, 21.02], 813 ]; 814 double[13][9] d11_expected = [ 815 [-78.32, -78.32, -78.32, -52.22, 13.45, 13.45, 13.45, 15.96, 17.25, 85.0, 85.0, 85.0, 85.0], 816 [-78.32, -78.32, -78.32, -51.54, 13.45, 13.45, 14.705, 15.96, 31.8, 85.0, 85.0, 85.0, 85.0], 817 [-78.32, -78.32, -78.32, -52.22, -50.86, 13.45, 13.45, 13.45, 17.25, 46.35, 85.0, 85.0, 85.0], 818 [-78.32, -78.32, -78.32, -52.22, -37.998, 8.3052, 13.45, 13.6508, 17.25, 54.08, 69.54, 78.816, 85.0], 819 [-78.32, -78.32, -70.49, -51.54, -5.843, 14.5042, 14.705, 14.9058, 31.8, 73.405, 85.0, 85.0, 85.0], 820 [-78.32, -78.32, -78.32, -51.88, -12.274, 14.4791, 14.705, 14.9309, 39.075, 85.0, 85.0, 85.0, 85.0], 821 [-78.32, -69.185, -60.05, -51.2, 0.588, 14.5293, 14.705, 14.8807, 24.525, 57.945, 71.4725, 79.589, 85.0], 822 [-78.32, -78.32, -73.97, -51.65333, -7.986667, 14.49583, 14.705, 14.91417, 34.225, 78.55833, 85.0, 85.0, 85.0], 823 [-78.32, -78.32, -73.1, -51.625, -7.45075, 14.49792, 14.705, 14.91208, 33.61875, 77.27, 85.0, 85.0, 85.0], 824 ]; 825 double[13][9] d12_expected = [ 826 [-81.36, -81.36, -81.36, -53.56, -42.14, -9.18, -9.18, -9.18, 49.52, 98.5, 98.5, 98.5, 98.5], 827 [-81.36, -81.36, -81.36, -53.56, -42.14, -9.18, -9.18, -9.18, 49.52, 98.5, 98.5, 98.5, 98.5], 828 [-81.36, -81.36, -81.36, -70.87, -42.14, -42.14, -42.14, -9.18, 49.52, 80.43, 98.5, 98.5, 98.5], 829 [-81.36, -81.36, -81.36, -66.5425, -46.708, -28.6264, -25.66, -22.6936, 38.9475, 82.237, 90.3685, 95.2474, 98.5], 830 [-81.36, -81.36, -77.164, -57.8875, -38.844, -12.1464, -9.18, -7.7031, 57.2475, 91.272, 98.5, 98.5, 98.5], 831 [-81.36, -81.36, -81.36, -62.215, -42.14, -12.476, -9.18, -7.539, 64.975, 98.5, 98.5, 98.5, 98.5], 832 [-81.36, -77.164, -72.968, -53.56, -35.548, -11.8168, -9.18, -7.8672, 49.52, 84.044, 91.272, 95.6088, 98.5], 833 [-81.36, -81.36, -78.56267, -59.33, -39.94267, -12.25627, -9.18, -7.6484, 59.82333, 93.68133, 98.5, 98.5, 98.5], 834 [-81.36, -81.36, -78.213, -58.96938, -39.668, -12.2288, -9.18, -7.662075, 59.17938, 93.079, 98.5, 98.5, 98.5], 835 ]; 836 double[13][9] d13_expected = [ 837 [38.37, 38.37, 38.37, 45.7, 50.69, 51.36, 51.36, 55.66, 58.95, 62.01, 65.25, 65.25, 65.25], 838 [38.37, 38.37, 41.365, 45.7, 51.025, 51.36, 53.51, 55.66, 58.95, 63.63, 65.25, 65.25, 65.25], 839 [38.37, 38.37, 38.37, 44.36, 50.69, 51.36, 51.36, 51.36, 58.95, 62.01, 65.25, 65.25, 65.25], 840 [38.37, 38.37, 38.37, 45.03, 50.69, 51.293, 51.36, 51.79, 57.93, 62.01, 63.63, 64.602, 65.25], 841 [38.37, 38.37, 41.365, 45.7, 51.025, 53.08, 53.51, 53.94, 58.95, 63.63, 65.25, 65.25, 65.25], 842 [38.37, 38.37, 38.969, 45.365, 50.958, 53.037, 53.51, 53.983, 59.715, 64.926, 65.25, 65.25, 65.25], 843 [38.37, 41.0655, 43.761, 46.9475, 51.092, 53.123, 53.51, 53.897, 58.44, 62.334, 63.792, 64.6668, 65.25], 844 [38.37, 38.37, 40.56633, 45.58833, 51.00267, 53.06567, 53.51, 53.95433, 59.205, 64.062, 65.25, 65.25, 65.25], 845 [38.37, 38.37, 40.766, 45.61625, 51.00825, 53.06925, 53.51, 53.95075, 59.14125, 63.954, 65.25, 65.25, 65.25], 846 ]; 847 848 void compareResults(const double[] actual, const double[] expected, string dataset, QuantileInterpolation method) 849 { 850 import std.conv : to; 851 import std.format : format; 852 import std.math : approxEqual, isNaN; 853 import std.range : lockstep; 854 855 foreach (i, actualValue, expectedValue; lockstep(actual, expected)) 856 { 857 assert(actualValue.approxEqual(expectedValue) || (actualValue.isNaN && expectedValue.isNaN), 858 format("Quantile unit test failure, dataset %s, method: %s, index: %d, expected: %g, actual: %g", 859 dataset, method.to!string, i, expectedValue, actualValue)); 860 } 861 } 862 863 foreach(methodIndex, method; EnumMembers!QuantileInterpolation) 864 { 865 compareResults(probs.map!(p => p.quantile(d1, method)).array, d1_expected[methodIndex], "d1", method); 866 compareResults(probs.map!(p => p.quantile(d2, method)).array, d2_expected[methodIndex], "d2", method); 867 compareResults(probs.map!(p => p.quantile(d3, method)).array, d3_expected[methodIndex], "d3", method); 868 compareResults(probs.map!(p => p.quantile(d3Int, method)).array, d3_expected[methodIndex], "d3Int", method); 869 compareResults(probs.map!(p => p.quantile(d4, method)).array, d4_expected[methodIndex], "d4", method); 870 compareResults(probs.map!(p => p.quantile(d4Int, method)).array, d4_expected[methodIndex], "d4Int", method); 871 compareResults(probs.map!(p => p.quantile(d5, method)).array, d5_expected[methodIndex], "d5", method); 872 compareResults(probs.map!(p => p.quantile(d5Int, method)).array, d5_expected[methodIndex], "d5Int", method); 873 compareResults(probs.map!(p => p.quantile(d6, method)).array, d6_expected[methodIndex], "d6", method); 874 compareResults(probs.map!(p => p.quantile(d6Size_t, method)).array, d6_expected[methodIndex], "d6Size_t", method); 875 compareResults(probs.map!(p => p.quantile(d7, method)).array, d7_expected[methodIndex], "d7", method); 876 compareResults(probs.map!(p => p.quantile(d7Float, method)).array, d7_expected[methodIndex], "d7Float", method); 877 compareResults(probs.map!(p => p.quantile(d8, method)).array, d8_expected[methodIndex], "d8", method); 878 compareResults(probs.map!(p => p.quantile(d8Float, method)).array, d8_expected[methodIndex], "d8Float", method); 879 compareResults(probs.map!(p => p.quantile(d9, method)).array, d9_expected[methodIndex], "d9", method); 880 compareResults(probs.map!(p => p.quantile(d9Float, method)).array, d9_expected[methodIndex], "d9Float", method); 881 compareResults(probs.map!(p => p.quantile(d10, method)).array, d10_expected[methodIndex], "d10", method); 882 compareResults(probs.map!(p => p.quantile(d10Float, method)).array, d10_expected[methodIndex], "d10Float", method); 883 compareResults(probs.map!(p => p.quantile(d11, method)).array, d11_expected[methodIndex], "d11", method); 884 compareResults(probs.map!(p => p.quantile(d12, method)).array, d12_expected[methodIndex], "d12", method); 885 compareResults(probs.map!(p => p.quantile(d13, method)).array, d13_expected[methodIndex], "d13", method); 886 } 887 }