1 /**
2   This module is part of bencher package, a simple rewrite
3   in the D programming language of the Rust standard benchmark
4   module.
5 
6   Original Rust source : https://github.com/rust-lang/rust/blob/master/src/libtest/stats.rs
7 
8   This module contains data structures and algorithms for doing 
9   statistics on benchmarks execution. The focus is put more on 
10   precision than execution speed. The algorithms are written to
11   fit for a particular use case and not intended to be used
12   elsewhere.
13 
14   Authors: Nouredine Hussain
15   Copyright: Copyright Nouredine Hussain 2017.
16   License: $(LINK3 http://www.boost.org/LICENSE_1_0.txt, Boost Software License - Version 1.0).
17 */
18 
19 module stats;
20 
21 import std.typecons : Tuple, tuple;
22 import std.range : empty;
23 import std.json;
24 
25 /** 
26   A structure that provides simple descriptive statistics
27   on a univariate set of numeric samples.
28 */
29 public struct Summary {
30 
31   /// Original rust comment : 
32   /// Note: this method sacrifices performance at the altar of accuracy
33   /// Depends on IEEE-754 arithmetic guarantees. See proof of correctness at:
34   /// ["Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates"]
35   /// (http://www.cs.cmu.edu/~quake-papers/robust-arithmetic.ps)
36   public double sum,
37     /// Original rust comment : 
38     /// Minimum value of the samples.
39     min,
40     /// Original rust comment : 
41     /// Maximum value of the samples.
42     max,
43     /// Original rust comment : 
44     /// Arithmetic mean (average) of the samples: sum divided by sample-count.
45     ///
46     /// See: https://en.wikipedia.org/wiki/Arithmetic_mean
47     mean,
48     /// Original rust comment : 
49     /// Median of the samples: value separating the lower half of the samples from the higher half.
50     /// Equal to `percentile(50.0)`.
51     ///
52     /// See: https://en.wikipedia.org/wiki/Median
53     median,
54     /// Original rust comment : 
55     /// Variance of the samples: bias-corrected mean of the squares of the differences of each
56     /// sample from the sample mean. Note that this calculates the _sample variance_ rather than the
57     /// population variance, which is assumed to be unknown. It therefore corrects the `(n-1)/n`
58     /// bias that would appear if we calculated a population variance, by dividing by `(n-1)` rather
59     /// than `n`.
60     ///
61     /// See: https://en.wikipedia.org/wiki/Variance
62     var,
63     /// Standard deviation: the square root of the sample variance.
64     ///
65     /// Note: this is not a robust statistic for non-normal distributions. Prefer the
66     /// `median_abs_dev` for unknown distributions.
67     ///
68     /// See: https://en.wikipedia.org/wiki/Standard_deviation
69     std_dev,
70     /// Original rust comment : 
71     /// Standard deviation as a percent of the mean value. See `std_dev` and `mean`.
72     ///
73     /// Note: this is not a robust statistic for non-normal distributions. Prefer the
74     /// `median_abs_dev_pct` for unknown distributions.
75     std_dev_pct,
76     /// Original rust comment : 
77     /// Scaled median of the absolute deviations of each sample from the sample median. This is a
78     /// robust (distribution-agnostic) estimator of sample variability. Use this in preference to
79     /// `std_dev` if you cannot assume your sample is normally distributed. Note that this is scaled
80     /// by the constant `1.4826` to allow its use as a consistent estimator for the standard
81     /// deviation.
82     ///
83     /// See: http://en.wikipedia.org/wiki/Median_absolute_deviation
84     median_abs_dev,
85     /// Original rust comment : 
86     /// Median absolute deviation as a percent of the median. See `median_abs_dev` and `median`.
87     median_abs_dev_pct,
88     /// Original rust comment : 
89     /// Inter-quartile range: the difference between the 25th percentile (1st quartile) and the 75th
90     /// percentile (3rd quartile). See `quartiles`.
91     ///
92     /// See also: https://en.wikipedia.org/wiki/Interquartile_range
93     iqr;
94 
95   /// Original rust comment : 
96   /// Quartiles of the sample: three values that divide the sample into four equal groups, each
97   /// with 1/4 of the data. The middle value is the median. See `median` and `percentile`. This
98   /// function may calculate the 3 quartiles more efficiently than 3 calls to `percentile`, but
99   /// is otherwise equivalent.
100   ///
101   /// See also: https://en.wikipedia.org/wiki/Quartile
102   public Tuple!(double, double, double) quartiles;
103 
104   /**
105     Construct a new summary of a sample set.
106   */
107   this(double[] samples){
108     import std.algorithm.iteration: reduce;
109     import std.algorithm.sorting: sort;
110     import std.math : sqrt, isNaN;
111     import std.array;
112 
113     double _min(double a, double b) {
114       if(a.isNaN)
115         return b;
116       if(b.isNaN)
117         return a;
118       return a < b ? a:b;
119     }
120 
121     double _max(double a, double b) {
122       if(a.isNaN)
123         return b;
124       if(b.isNaN)
125         return a;
126       return a > b ? a:b;
127     }
128 
129     bool _less(double a, double b) {
130       import std.math : isNaN;
131 
132       if(a.isNaN)
133         return false;
134       if(b.isNaN)
135         return true;
136       return a < b;
137 
138     }
139 
140     double[] sorted_samples = samples.sort!(_less).array;
141     sum = samples.sum();
142     min = samples.reduce!(_min);
143     max = samples.reduce!(_max);
144     mean = sum / samples.length;
145     median = sorted_samples.median_of_sorted();
146     var = samples.var(mean);
147     std_dev = var.sqrt();
148     std_dev_pct = (std_dev/mean) * 100.0;
149     median_abs_dev = sorted_samples.median_abs_dev_of_sorted(median);
150     median_abs_dev_pct = median_abs_dev.median_abs_dev_pct(median);
151     quartiles = sorted_samples.quartiles_of_sorted();
152     iqr = quartiles[2] - quartiles[0];
153   }
154 
155   /** 
156     Serialize the summary in JSON format.
157   */
158 
159   public JSONValue toJSON(){
160 
161     JSONValue jj = ["sum": sum];
162     jj.object["min"] = JSONValue(min);
163     jj.object["max"] = JSONValue(max);
164     jj.object["mean"] = JSONValue(mean);
165     jj.object["median"] = JSONValue(median);
166     jj.object["var"] = JSONValue(var);
167     jj.object["std_dev"] = JSONValue(std_dev);
168     jj.object["std_dev_pct"] = JSONValue(std_dev_pct);
169     jj.object["median_abs_dev"] = JSONValue(median_abs_dev);
170     jj.object["median_abs_dev_pct"] = JSONValue(median_abs_dev_pct);
171     jj.object["iqr"] = JSONValue(iqr);
172     jj.object["quartiles"] = JSONValue([quartiles[0], quartiles[1], quartiles[2]]);
173     return jj;
174   }
175 
176 }
177 
178 private alias sumKahan sum;
179 
180 // Kahan algo http://en.wikipedia.org/wiki/Kahan_summation_algorithm
181 private double sumKahan(double[] arr) {
182     import std.math : abs;
183     import std.algorithm.mutation: swap;
184     import std.algorithm.iteration: fold;
185     double[] partials;
186 
187     foreach(_, x; arr) {
188         ulong j = 0;
189         for(int i=0; i < partials.length; ++i) {
190             auto y = partials[i];
191             if(x.abs() < y.abs()) {
192                 swap(x, y);
193             }
194             immutable hi = x + y;
195             immutable lo = y - (hi - x);
196             if(lo != 0.0) {
197               partials[j] = lo;
198               j++;
199             }
200             x = hi;
201         }
202         if(j >= partials.length){
203             partials ~= x;
204         } else {
205             partials[j] = x;
206             partials.length = j+1;
207         }
208     }
209     return partials.fold!((a, b) => a+b);
210 }
211 
212 
213 /// Original rust comment : 
214 /// Percentile: the value below which `pct` percent of the values in the samples fall. For example,
215 /// percentile(95.0) will return the value `v` such that 95% of the samples `s` in the samples
216 /// satisfy `s <= v`.
217 ///
218 /// Calculated by linear interpolation between closest ranks.
219 ///
220 /// See: http://en.wikipedia.org/wiki/Percentile
221 private double percentile_of_sorted(double[] sorted_arr, double pct){
222   import std.math: floor;
223 
224   assert(!sorted_arr.empty);
225 
226   if(sorted_arr.length==1)
227     return sorted_arr[0];
228   if(pct == 100.0)
229     return sorted_arr[$-1];
230   assert(0.0 <= pct && pct <= 100.0);
231   immutable length = sorted_arr.length -1;
232   immutable rank = (pct / 100.0) * length;
233   immutable lrank = rank.floor();
234   immutable d = rank - lrank;
235   immutable n = cast(size_t)lrank;
236   immutable lo = sorted_arr[n];
237   immutable hi = sorted_arr[n+1];
238   return lo + (hi -lo) *d;
239 
240 
241 }
242 
243 private double median_of_sorted(double[] sorted_arr) {
244   return percentile_of_sorted(sorted_arr, 50.0);
245 }
246 
247 private double var(double[] arr, immutable double mean) {
248   if(arr.length < 2)
249     return 0.0;
250   double v = 0;
251   foreach(_, ref e; arr){
252     auto x = e - mean;
253     v = v + x*x;
254   }
255   auto denom = arr.length - 1;
256 
257   return v / denom;
258 }
259 
260 
261 private double median_abs_dev_of_sorted(double[] sorted, immutable double median) {
262   import std.algorithm.iteration: map;
263   import std.algorithm.sorting: sort;
264   import std.math: abs;
265   import std.array;
266 
267   static bool _less(double a, double b) {
268     import std.math : isNaN;
269 
270     if(a.isNaN)
271       return false;
272     if(b.isNaN)
273       return true;
274     return a < b;
275 
276   }
277   auto abs_devs = sorted.map!(a => abs(median-a))
278                         .array
279                         .sort!(_less)
280                         .array;
281 
282   // from rust dev comment:
283   // This constant is derived by smarter statistics brains than me, but it is
284   // consistent with how R and other packages treat the MAD.
285   immutable number = 1.4826;
286   return abs_devs.median_of_sorted() * number;
287 }
288 
289 private double median_abs_dev_pct(double median_abs_dev, double median) {
290   return (median_abs_dev/median)*100.0;
291 }
292 
293 private double iqr(double[] samples) {
294   return 0;
295 }
296 
297 private auto quartiles_of_sorted(double[] sorted) {
298   auto a = sorted.percentile_of_sorted(25.0);
299   auto b = sorted.percentile_of_sorted(50.0);
300   auto c = sorted.percentile_of_sorted(75.0);
301 
302   return tuple(a, b, c);
303 }
304 
305 /// Original rust comment
306 /// Winsorize a set of samples, replacing values above the `100-pct` percentile
307 /// and below the `pct` percentile with those percentiles themselves. This is a
308 /// way of minimizing the effect of outliers, at the cost of biasing the sample.
309 /// It differs from trimming in that it does not change the number of samples,
310 /// just changes the values of those that are outliers.
311 ///
312 /// See: http://en.wikipedia.org/wiki/Winsorising
313 public void winsorize(double[] arr, double pct){
314   import std.algorithm.sorting: sort;
315   import std.array;
316 
317   bool _less(double a, double b) {
318     import std.math : isNaN;
319 
320     if(a.isNaN)
321       return false;
322     if(b.isNaN)
323       return true;
324     return a < b;
325 
326   }
327 
328   auto sorted = arr.array.sort!(_less).array();
329   immutable lo = percentile_of_sorted(sorted , pct);
330   immutable hi = percentile_of_sorted(sorted, 100.0-pct);
331   foreach(_, ref e; arr){
332     if(e > hi){
333       e = hi;
334     } else if(e < lo) {
335       e = lo;
336     }
337   }
338 }
339 
340 
341 unittest 
342 {
343   double[5] samples = [1.0, 2.0, double.nan, 4.0, 3.0];
344 
345   auto summ = Summary(samples);
346 
347   assert(summ.min == 1.0);
348   assert(summ.max == 4.0);
349 }
350 
351 private void check(double[] samples, ref Summary summ) {
352     import std.math: approxEqual;
353     auto summ2 = Summary(samples);
354 
355     assert(summ.sum == summ2.sum);
356     assert(summ.min == summ2.min);
357     assert(summ.max == summ2.max);
358     assert(summ.mean == summ2.mean);
359     assert(summ.median == summ2.median);
360     assert(summ.quartiles == summ2.quartiles);
361     assert(summ.iqr == summ2.iqr);
362 
363     assert(approxEqual(summ.var, summ2.var));
364     assert(approxEqual(summ.std_dev, summ2.std_dev));
365     assert(approxEqual(summ.std_dev_pct, summ2.std_dev_pct));
366     assert(approxEqual(summ.median_abs_dev, summ2.median_abs_dev));
367     assert(approxEqual(summ.median_abs_dev_pct, summ2.median_abs_dev_pct));
368 }
369 
370 unittest
371 {
372   // norm2
373   auto samples = [958.0000000000, 
374                   924.0000000000];
375 
376   auto summ = Summary();
377   summ.sum= 1882.0000000000;
378   summ.min= 924.0000000000;
379   summ.max= 958.0000000000;
380   summ.mean= 941.0000000000;
381   summ.median= 941.0000000000;
382   summ.var= 578.0000000000;
383   summ.std_dev= 24.0416305603;
384   summ.std_dev_pct= 2.5549022912;
385   summ.median_abs_dev= 25.2042000000;
386   summ.median_abs_dev_pct= 2.6784484591;
387   summ.quartiles= tuple(932.5000000000, 941.0000000000, 949.5000000000);
388   summ.iqr= 17.0000000000;
389 
390   check(samples, summ);
391 }
392 
393 unittest
394 {
395   // norm10narrow
396   auto samples = [966.0000000000,
397                   985.0000000000,
398                   1110.0000000000,
399                   848.0000000000,
400                   821.0000000000,
401                   975.0000000000,
402                   962.0000000000,
403                   1157.0000000000,
404                   1217.0000000000,
405                   955.0000000000];
406 
407   auto summ = Summary();
408   summ.sum= 9996.0000000000;
409   summ.min= 821.0000000000;
410   summ.max= 1217.0000000000;
411   summ.mean= 999.6000000000;
412   summ.median= 970.5000000000;
413   summ.var= 16050.7111111111;
414   summ.std_dev= 126.6914010938;
415   summ.std_dev_pct= 12.6742097933;
416   summ.median_abs_dev= 102.2994000000;
417   summ.median_abs_dev_pct= 10.5408964451;
418   summ.quartiles= tuple(956.7500000000, 970.5000000000, 1078.7500000000);
419   summ.iqr= 122.0000000000;
420 
421   check(samples, summ);
422 }
423 
424 
425 unittest
426 {
427   // norm10medium
428   auto samples = [954.0000000000,
429                   1064.0000000000,
430                   855.0000000000,
431                   1000.0000000000,
432                   743.0000000000,
433                   1084.0000000000,
434                   704.0000000000,
435                   1023.0000000000,
436                   357.0000000000,
437                   869.0000000000];
438 
439   auto summ = Summary();
440   summ.sum= 8653.0000000000;
441   summ.min= 357.0000000000;
442   summ.max= 1084.0000000000;
443   summ.mean= 865.3000000000;
444   summ.median= 911.5000000000;
445   summ.var= 48628.4555555556;
446   summ.std_dev= 220.5186059170;
447   summ.std_dev_pct= 25.4846418487;
448   summ.median_abs_dev= 195.7032000000;
449   summ.median_abs_dev_pct= 21.4704552935;
450   summ.quartiles= tuple(771.0000000000, 911.5000000000, 1017.2500000000);
451   summ.iqr= 246.2500000000;
452 
453   check(samples, summ);
454 }
455 
456 unittest
457 {
458   // norm10wide
459   auto samples = [505.0000000000,
460                   497.0000000000,
461                   1591.0000000000,
462                   887.0000000000,
463                   1026.0000000000,
464                   136.0000000000,
465                   1580.0000000000,
466                   940.0000000000,
467                   754.0000000000,
468                   1433.0000000000];
469 
470   auto summ = Summary();
471   summ.sum= 9349.0000000000;
472   summ.min= 136.0000000000;
473   summ.max= 1591.0000000000;
474   summ.mean= 934.9000000000;
475   summ.median= 913.5000000000;
476   summ.var= 239208.9888888889;
477   summ.std_dev= 489.0899599142;
478   summ.std_dev_pct= 52.3146817750;
479   summ.median_abs_dev= 611.5725000000;
480   summ.median_abs_dev_pct= 66.9482758621;
481   summ.quartiles= tuple(567.2500000000, 913.5000000000, 1331.2500000000);
482   summ.iqr= 764.0000000000;
483 
484   check(samples, summ);
485 }
486 
487 unittest
488 {
489   // norm25verynarrow
490   auto samples = [991.0000000000,
491                   1018.0000000000,
492                   998.0000000000,
493                   1013.0000000000,
494                   974.0000000000,
495                   1007.0000000000,
496                   1014.0000000000,
497                   999.0000000000,
498                   1011.0000000000,
499                   978.0000000000,
500                   985.0000000000,
501                   999.0000000000,
502                   983.0000000000,
503                   982.0000000000,
504                   1015.0000000000,
505                   1002.0000000000,
506                   977.0000000000,
507                   948.0000000000,
508                   1040.0000000000,
509                   974.0000000000,
510                   996.0000000000,
511                   989.0000000000,
512                   1015.0000000000,
513                   994.0000000000,
514                   1024.0000000000];
515 
516   auto summ = Summary();
517   summ.sum= 24926.0000000000;
518   summ.min= 948.0000000000;
519   summ.max= 1040.0000000000;
520   summ.mean= 997.0400000000;
521   summ.median= 998.0000000000;
522   summ.var= 393.2066666667;
523   summ.std_dev= 19.8294393937;
524   summ.std_dev_pct= 1.9888308788;
525   summ.median_abs_dev= 22.2390000000;
526   summ.median_abs_dev_pct= 2.2283567134;
527   summ.quartiles= tuple(983.0000000000, 998.0000000000, 1013.0000000000);
528   summ.iqr= 30.0000000000;
529 
530   check(samples, summ);
531 }
532 
533 unittest
534 {
535   // exp10a
536   auto samples = [23.0000000000,
537                   11.0000000000,
538                   2.0000000000,
539                   57.0000000000,
540                   4.0000000000,
541                   12.0000000000,
542                   5.0000000000,
543                   29.0000000000,
544                   3.0000000000,
545                   21.0000000000];
546 
547   auto summ = Summary();
548   summ.sum= 167.0000000000;
549   summ.min= 2.0000000000;
550   summ.max= 57.0000000000;
551   summ.mean= 16.7000000000;
552   summ.median= 11.5000000000;
553   summ.var= 287.7888888889;
554   summ.std_dev= 16.9643416875;
555   summ.std_dev_pct= 101.5828843560;
556   summ.median_abs_dev= 13.3434000000;
557   summ.median_abs_dev_pct= 116.0295652174;
558   summ.quartiles= tuple(4.2500000000, 11.5000000000, 22.5000000000);
559   summ.iqr= 18.2500000000;
560         
561   check(samples, summ);
562 }
563 
564 unittest
565 {
566   // exp10b
567   auto samples = [24.0000000000,
568                   17.0000000000,
569                   6.0000000000,
570                   38.0000000000,
571                   25.0000000000,
572                   7.0000000000,
573                   51.0000000000,
574                   2.0000000000,
575                   61.0000000000,
576                   32.0000000000];
577 
578   auto summ = Summary();
579   summ.sum= 263.0000000000;
580   summ.min= 2.0000000000;
581   summ.max= 61.0000000000;
582   summ.mean= 26.3000000000;
583   summ.median= 24.5000000000;
584   summ.var= 383.5666666667;
585   summ.std_dev= 19.5848580967;
586   summ.std_dev_pct= 74.4671410520;
587   summ.median_abs_dev= 22.9803000000;
588   summ.median_abs_dev_pct= 93.7971428571;
589   summ.quartiles= tuple(9.5000000000, 24.5000000000, 36.5000000000);
590   summ.iqr= 27.0000000000;
591         
592   check(samples, summ);
593 }
594 
595 unittest
596 {
597   // exp10c
598   auto samples = [71.0000000000,
599                   2.0000000000,
600                   32.0000000000,
601                   1.0000000000,
602                   6.0000000000,
603                   28.0000000000,
604                   13.0000000000,
605                   37.0000000000,
606                   16.0000000000,
607                   36.0000000000];
608 
609   auto summ = Summary();
610   summ.sum= 242.0000000000;
611   summ.min= 1.0000000000;
612   summ.max= 71.0000000000;
613   summ.mean= 24.2000000000;
614   summ.median= 22.0000000000;
615   summ.var= 458.1777777778;
616   summ.std_dev= 21.4050876611;
617   summ.std_dev_pct= 88.4507754589;
618   summ.median_abs_dev= 21.4977000000;
619   summ.median_abs_dev_pct= 97.7168181818;
620   summ.quartiles= tuple(7.7500000000, 22.0000000000, 35.0000000000);
621   summ.iqr= 27.2500000000;
622         
623   check(samples, summ);
624 }
625 
626 unittest
627 {
628   // exp25
629   auto samples = [3.0000000000,
630                   24.0000000000,
631                   1.0000000000,
632                   19.0000000000,
633                   7.0000000000,
634                   5.0000000000,
635                   30.0000000000,
636                   39.0000000000,
637                   31.0000000000,
638                   13.0000000000,
639                   25.0000000000,
640                   48.0000000000,
641                   1.0000000000,
642                   6.0000000000,
643                   42.0000000000,
644                   63.0000000000,
645                   2.0000000000,
646                   12.0000000000,
647                   108.0000000000,
648                   26.0000000000,
649                   1.0000000000,
650                   7.0000000000,
651                   44.0000000000,
652                   25.0000000000,
653                   11.0000000000];
654 
655   auto summ = Summary();
656   summ.sum= 593.0000000000;
657   summ.min= 1.0000000000;
658   summ.max= 108.0000000000;
659   summ.mean= 23.7200000000;
660   summ.median= 19.0000000000;
661   summ.var= 601.0433333333;
662   summ.std_dev= 24.5161851301;
663   summ.std_dev_pct= 103.3565983562;
664   summ.median_abs_dev= 19.2738000000;
665   summ.median_abs_dev_pct= 101.4410526316;
666   summ.quartiles= tuple(6.0000000000, 19.0000000000, 31.0000000000);
667   summ.iqr= 25.0000000000;
668         
669   check(samples, summ);
670 }
671 
672 unittest
673 {
674   // binom25
675   auto samples = [18.0000000000,
676                   17.0000000000,
677                   27.0000000000,
678                   15.0000000000,
679                   21.0000000000,
680                   25.0000000000,
681                   17.0000000000,
682                   24.0000000000,
683                   25.0000000000,
684                   24.0000000000,
685                   26.0000000000,
686                   26.0000000000,
687                   23.0000000000,
688                   15.0000000000,
689                   23.0000000000,
690                   17.0000000000,
691                   18.0000000000,
692                   18.0000000000,
693                   21.0000000000,
694                   16.0000000000,
695                   15.0000000000,
696                   31.0000000000,
697                   20.0000000000,
698                   17.0000000000,
699                   15.0000000000];
700 
701   auto summ = Summary();
702   summ.sum= 514.0000000000;
703   summ.min= 15.0000000000;
704   summ.max= 31.0000000000;
705   summ.mean= 20.5600000000;
706   summ.median= 20.0000000000;
707   summ.var= 20.8400000000;
708   summ.std_dev= 4.5650848842;
709   summ.std_dev_pct= 22.2037202539;
710   summ.median_abs_dev= 5.9304000000;
711   summ.median_abs_dev_pct= 29.6520000000;
712   summ.quartiles= tuple(17.0000000000, 20.0000000000, 24.0000000000);
713   summ.iqr= 7.0000000000;
714 
715   check(samples, summ);
716 }
717 
718 unittest
719 {
720   // pois25lambda30
721   auto samples = [27.0000000000,
722                   33.0000000000,
723                   34.0000000000,
724                   34.0000000000,
725                   24.0000000000,
726                   39.0000000000,
727                   28.0000000000,
728                   27.0000000000,
729                   31.0000000000,
730                   28.0000000000,
731                   38.0000000000,
732                   21.0000000000,
733                   33.0000000000,
734                   36.0000000000,
735                   29.0000000000,
736                   37.0000000000,
737                   32.0000000000,
738                   34.0000000000,
739                   31.0000000000,
740                   39.0000000000,
741                   25.0000000000,
742                   31.0000000000,
743                   32.0000000000,
744                   40.0000000000,
745                   24.0000000000];
746 
747   auto summ = Summary();
748   summ.sum= 787.0000000000;
749   summ.min= 21.0000000000;
750   summ.max= 40.0000000000;
751   summ.mean= 31.4800000000;
752   summ.median= 32.0000000000;
753   summ.var= 26.5933333333;
754   summ.std_dev= 5.1568724372;
755   summ.std_dev_pct= 16.3814245145;
756   summ.median_abs_dev= 5.9304000000;
757   summ.median_abs_dev_pct= 18.5325000000;
758   summ.quartiles= tuple(28.0000000000, 32.0000000000, 34.0000000000);
759   summ.iqr= 6.0000000000;
760 
761   check(samples, summ);
762 }
763 
764 unittest
765 {
766   // pois25lambda40
767   auto samples = [42.0000000000,
768                   50.0000000000,
769                   42.0000000000,
770                   46.0000000000,
771                   34.0000000000,
772                   45.0000000000,
773                   34.0000000000,
774                   49.0000000000,
775                   39.0000000000,
776                   28.0000000000,
777                   40.0000000000,
778                   35.0000000000,
779                   37.0000000000,
780                   39.0000000000,
781                   46.0000000000,
782                   44.0000000000,
783                   32.0000000000,
784                   45.0000000000,
785                   42.0000000000,
786                   37.0000000000,
787                   48.0000000000,
788                   42.0000000000,
789                   33.0000000000,
790                   42.0000000000,
791                   48.0000000000];
792 
793   auto summ = Summary();
794   summ.sum= 1019.0000000000;
795   summ.min= 28.0000000000;
796   summ.max= 50.0000000000;
797   summ.mean= 40.7600000000;
798   summ.median= 42.0000000000;
799   summ.var= 34.4400000000;
800   summ.std_dev= 5.8685603004;
801   summ.std_dev_pct= 14.3978417577;
802   summ.median_abs_dev= 5.9304000000;
803   summ.median_abs_dev_pct= 14.1200000000;
804   summ.quartiles= tuple(37.0000000000, 42.0000000000, 45.0000000000);
805   summ.iqr= 8.0000000000;
806 
807   check(samples, summ);
808 }
809 
810 
811 unittest
812 {
813   // pois25lambda50
814   auto samples = [45.0000000000,
815                   43.0000000000,
816                   44.0000000000,
817                   61.0000000000,
818                   51.0000000000,
819                   53.0000000000,
820                   59.0000000000,
821                   52.0000000000,
822                   49.0000000000,
823                   51.0000000000,
824                   51.0000000000,
825                   50.0000000000,
826                   49.0000000000,
827                   56.0000000000,
828                   42.0000000000,
829                   52.0000000000,
830                   51.0000000000,
831                   43.0000000000,
832                   48.0000000000,
833                   48.0000000000,
834                   50.0000000000,
835                   42.0000000000,
836                   43.0000000000,
837                   42.0000000000,
838                   60.0000000000];
839 
840   auto summ = Summary();
841   summ.sum= 1235.0000000000;
842   summ.min= 42.0000000000;
843   summ.max= 61.0000000000;
844   summ.mean= 49.4000000000;
845   summ.median= 50.0000000000;
846   summ.var= 31.6666666667;
847   summ.std_dev= 5.6273143387;
848   summ.std_dev_pct= 11.3913245723;
849   summ.median_abs_dev= 4.4478000000;
850   summ.median_abs_dev_pct= 8.8956000000;
851   summ.quartiles= tuple(44.0000000000, 50.0000000000, 52.0000000000);
852   summ.iqr= 8.0000000000;
853 
854   check(samples, summ);
855 }
856 
857 
858 unittest
859 {
860   // unif25
861   auto samples = [99.0000000000,
862                   55.0000000000,
863                   92.0000000000,
864                   79.0000000000,
865                   14.0000000000,
866                   2.0000000000,
867                   33.0000000000,
868                   49.0000000000,
869                   3.0000000000,
870                   32.0000000000,
871                   84.0000000000,
872                   59.0000000000,
873                   22.0000000000,
874                   86.0000000000,
875                   76.0000000000,
876                   31.0000000000,
877                   29.0000000000,
878                   11.0000000000,
879                   41.0000000000,
880                   53.0000000000,
881                   45.0000000000,
882                   44.0000000000,
883                   98.0000000000,
884                   98.0000000000,
885                   7.0000000000];
886 
887   auto summ = Summary();
888   summ.sum= 1242.0000000000;
889   summ.min= 2.0000000000;
890   summ.max= 99.0000000000;
891   summ.mean= 49.6800000000;
892   summ.median= 45.0000000000;
893   summ.var= 1015.6433333333;
894   summ.std_dev= 31.8691595957;
895   summ.std_dev_pct= 64.1488719719;
896   summ.median_abs_dev= 45.9606000000;
897   summ.median_abs_dev_pct= 102.1346666667;
898   summ.quartiles= tuple(29.0000000000, 45.0000000000, 79.0000000000);
899   summ.iqr= 50.0000000000;
900 
901   check(samples, summ);
902 }
903 
904 unittest
905 {
906   assert([0.5, 3.2321, 1.5678].sum() == 5.2999);
907 }
908 
909 unittest
910 {
911   import std.stdio : writeln;
912 
913   // writeln(1e30 + 1.2);
914   // writeln(1.2 + 1e30);
915 
916   // writeln(1e30 - 1.2);
917   // writeln(1.2 - 1e30);
918 
919 
920   assert([1e30, 1.2, -1e30].sum() == 1.2);
921 }