1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21 


22  using System;


23  using System.Collections.Generic;


24  using System.Diagnostics.Contracts;


25  using System.Linq;


26 


27  namespace HeuristicLab.Common {


28  public static class EnumerableStatisticExtensions {


29  /// <summary>


30  /// Calculates the median element of the enumeration.


31  /// </summary>


32  /// <param name="values"></param>


33  /// <returns></returns>


34  public static double Median(this IEnumerable<double> values) {


35  // See unit tests for comparison with naive implementation


36  return Quantile(values, 0.5);


37  }


38 


39  /// <summary>


40  /// Calculates the alphaquantile element of the enumeration.


41  /// </summary>


42  /// <param name="values"></param>


43  /// <returns></returns>


44  public static double Quantile(this IEnumerable<double> values, double alpha) {


45  // See unit tests for comparison with naive implementation


46  double[] valuesArr = values.ToArray();


47  int n = valuesArr.Length;


48  if (n == 0) throw new InvalidOperationException("Enumeration contains no elements.");


49 


50  // "When N is even, statistics books define the median as the arithmetic mean of the elements k = N/2


51  // and k = N/2 + 1 (that is, N/2 from the bottom and N/2 from the top).


52  // If you accept such pedantry, you must perform two separate selections to find these elements."


53 


54  // return the element at Math.Ceiling (if n*alpha is fractional) or the average of two elements if n*alpha is integer.


55  var pos = n * alpha;


56  Contract.Assert(pos >= 0);


57  Contract.Assert(pos < n);


58  bool isInteger = Math.Round(pos).IsAlmost(pos);


59  if (isInteger) {


60  return 0.5 * (Select((int)pos  1, valuesArr) + Select((int)pos, valuesArr));


61  } else {


62  return Select((int)Math.Ceiling(pos)  1, valuesArr);


63  }


64  }


65 


66  // Numerical Recipes in C++, §8.5 Selecting the Mth Largest, O(n)


67  // Giben k in [0..n1] returns an array value from array arr[0..n1] such that k array values are


68  // lee than or equal to the one returned. The input array will be rearranged to hav this value in


69  // location arr[k], with all smaller elements moved to arr[0..k1] (in arbitrary order) and all


70  // larger elements in arr[k+1..n1] (also in arbitrary order).


71  private static double Select(int k, double[] arr) {


72  Contract.Assert(arr.GetLowerBound(0) == 0);


73  Contract.Assert(k >= 0 && k < arr.Length);


74  int i, ir, j, l, mid, n = arr.Length;


75  double a;


76  l = 0;


77  ir = n  1;


78  for (; ; ) {


79  if (ir <= l + 1) {


80  // Active partition contains 1 or 2 elements.


81  if (ir == l + 1 && arr[ir] < arr[l]) {


82  // if (ir == l + 1 && arr[ir].CompareTo(arr[l]) < 0) {


83  // Case of 2 elements.


84  // SWAP(arr[l], arr[ir]);


85  double temp = arr[l];


86  arr[l] = arr[ir];


87  arr[ir] = temp;


88  }


89  return arr[k];


90  } else {


91  mid = (l + ir) >> 1; // Choose median of left, center, and right elements


92  {


93  // SWAP(arr[mid], arr[l + 1]); // as partitioning element a. Also


94  double temp = arr[mid];


95  arr[mid] = arr[l + 1];


96  arr[l + 1] = temp;


97  }


98 


99  if (arr[l] > arr[ir]) {


100  // if (arr[l].CompareTo(arr[ir]) > 0) { // rearrange so that arr[l] arr[ir] <= arr[l+1],


101  // SWAP(arr[l], arr[ir]); . arr[ir] >= arr[l+1]


102  double temp = arr[l];


103  arr[l] = arr[ir];


104  arr[ir] = temp;


105  }


106 


107  if (arr[l + 1] > arr[ir]) {


108  // if (arr[l + 1].CompareTo(arr[ir]) > 0) {


109  // SWAP(arr[l + 1], arr[ir]);


110  double temp = arr[l + 1];


111  arr[l + 1] = arr[ir];


112  arr[ir] = temp;


113  }


114  if (arr[l] > arr[l + 1]) {


115  //if (arr[l].CompareTo(arr[l + 1]) > 0) {


116  // SWAP(arr[l], arr[l + 1]);


117  double temp = arr[l];


118  arr[l] = arr[l + 1];


119  arr[l + 1] = temp;


120 


121  }


122  i = l + 1; // Initialize pointers for partitioning.


123  j = ir;


124  a = arr[l + 1]; // Partitioning element.


125  for (; ; ) { // Beginning of innermost loop.


126  do i++; while (arr[i] < a /* arr[i].CompareTo(a) < 0 */); // Scan up to find element > a.


127  do j; while (arr[j] > a /* arr[j].CompareTo(a) > 0 */); // Scan down to find element < a.


128  if (j < i) break; // Pointers crossed. Partitioning complete.


129  {


130  // SWAP(arr[i], arr[j]);


131  double temp = arr[i];


132  arr[i] = arr[j];


133  arr[j] = temp;


134  }


135  } // End of innermost loop.


136  arr[l + 1] = arr[j]; // Insert partitioning element.


137  arr[j] = a;


138  if (j >= k) ir = j  1; // Keep active the partition that contains the


139  if (j <= k) l = i; // kth element.


140  }


141  }


142  }


143 


144  /// <summary>


145  /// Calculates the range (max  min) of the enumeration.


146  /// </summary>


147  /// <param name="values"></param>


148  /// <returns></returns>


149  public static double Range(this IEnumerable<double> values) {


150  double min = double.PositiveInfinity;


151  double max = double.NegativeInfinity;


152  int i = 0;


153  foreach (var e in values) {


154  if (min > e) min = e;


155  if (max < e) max = e;


156  i++;


157  }


158  if (i < 1) throw new ArgumentException("The enumerable must contain at least two elements", "values");


159  return max  min;


160  }


161 


162 


163  /// <summary>


164  /// Calculates the standard deviation of values.


165  /// </summary>


166  /// <param name="values"></param>


167  /// <returns></returns>


168  public static double StandardDeviation(this IEnumerable<double> values) {


169  return Math.Sqrt(Variance(values));


170  }


171 


172  /// <summary>


173  /// Calculates the variance of values. (sum (x  x_mean)² / n)


174  /// </summary>


175  /// <param name="values"></param>


176  /// <returns></returns>


177  public static double Variance(this IEnumerable<double> values) {


178  int m_n = 0;


179  double m_oldM = 0.0;


180  double m_newM = 0.0;


181  double m_oldS = 0.0;


182  double m_newS = 0.0;


183  foreach (double x in values) {


184  m_n++;


185  if (m_n == 1) {


186  m_oldM = m_newM = x;


187  m_oldS = 0.0;


188  } else {


189  m_newM = m_oldM + (x  m_oldM) / m_n;


190  m_newS = m_oldS + (x  m_oldM) * (x  m_newM);


191 


192  // set up for next iteration


193  m_oldM = m_newM;


194  m_oldS = m_newS;


195  }


196  }


197  return ((m_n > 1) ? m_newS / (m_n  1) : 0.0);


198  }


199 


200  /// <summary>


201  /// Calculates the pth percentile of the values.


202  /// </summary>


203  public static double Percentile(this IEnumerable<double> values, double p) {


204  // iterate only once


205  double[] valuesArr = values.ToArray();


206  int n = valuesArr.Length;


207  if (n == 0) throw new InvalidOperationException("Enumeration contains no elements.");


208  if (n == 1) return values.ElementAt(0);


209 


210  if (p.IsAlmost(0.0)) return valuesArr[0];


211  if (p.IsAlmost(1.0)) return valuesArr[n  1];


212 


213  double t = p * (n  1);


214  int index = (int)Math.Floor(t);


215  double percentage = t  index;


216  return valuesArr[index] * (1  percentage) + valuesArr[index + 1] * percentage;


217  }


218 


219  public static IEnumerable<double> LimitToRange(this IEnumerable<double> values, double min, double max) {


220  if (min > max) throw new ArgumentException(string.Format("Minimum {0} is larger than maximum {1}.", min, max));


221  foreach (var x in values) {


222  if (double.IsNaN(x)) yield return (max + min) / 2.0;


223  else if (x < min) yield return min;


224  else if (x > max) yield return max;


225  else yield return x;


226  }


227  }


228  }


229  }

