feat: use more performant implementations for calculating mean, mad,aad and variance (#413)

Author: Uzlopak

src/utils.ts

@@ -253,21 +253,6 @@ export const isFnAsyncResource = (fn: Fn | null | undefined): boolean => {
   }
 }
 
-/**
- * Computes the average of a sample.
- * @param samples - the sample
- * @returns the average of the sample
- */
-const average = (samples: Samples) => {
-  let result = 0
-
-  for (const sample of samples) {
-    result += sample
-  }
-
-  return result / samples.length
-}
-
 /**
  * A type representing a samples-array with at least one number.
  */
@@ -286,14 +271,6 @@ export const isValidSamples = (
   return Array.isArray(value) && value.length !== 0
 }
 
-/**
- * Sorts samples and returns a new sorted array.
- * @param samples - samples to sort
- * @returns new sorted samples
- */
-export const toSortedSamples = (samples: Samples): SortedSamples =>
-  [...samples].sort(sortFn) as SortedSamples
-
 /**
  * Sorts samples in place.
  * @param samples - samples to sort
@@ -304,21 +281,29 @@ export function sortSamples (
   samples.sort(sortFn)
 }
 
-/**
- * Computes the variance of a sample with Bessel's correction.
- * @param samples - the sample
- * @param avg - the average of the sample
- * @returns the variance of the sample
- */
-const variance = (samples: Samples, avg = average(samples)) => {
-  if (samples.length === 1) {
-    return 0
+export const meanAndVariance = (samples: Samples): { mean: number; vr: number } => {
+  const len = samples.length
+  if (len === 1) {
+    return { mean: samples[0], vr: 0 }
+  }
+
+  let mean = 0
+  let m = 0
+  let x = 0
+  let d = 0
+  let i = 0
+
+  while (i < len) {
+    x = samples[i++]! // eslint-disable-line @typescript-eslint/no-non-null-assertion
+    d = x - mean
+    mean += d / i
+    m += d * (x - mean)
   }
-  let sumSq = 0
-  for (const sample of samples) {
-    sumSq += (sample - avg) ** 2
+
+  return {
+    mean,
+    vr: m / (len - 1)
   }
-  return sumSq / (samples.length - 1)
 }
 
 type Quantile = 0.5 | 0.75 | 0.99 | 0.995 | 0.999
@@ -342,47 +327,78 @@ const quantileSorted = (samples: SortedSamples, q: Quantile): number => {
     : samples[baseIndex]!
 }
 
-/**
- * Computes the median of a sorted sample.
- * @param samples - the sorted sample
- * @returns the median of the sample
- */
-const medianSorted = (samples: SortedSamples) => quantileSorted(samples, 0.5)
-
 /**
  * A sort function to be passed to Array.prototype.sort for numbers.
  * @param a - first number
  * @param b - second number
  * @returns a number indicating the sort order
  */
-const sortFn = (a: number, b: number) => a - b
+export const sortFn = (a: number, b: number) => a - b
 
 /**
- * Computes the median of an unsorted sample.
+ * Computes the average absolute deviation from the mean.
  * @param samples - the sample
- * @returns the median of the sample
+ * @param mean - the mean of the sample
+ * @returns the average absolute deviation
  */
-const median = (samples: Samples) => medianSorted(toSortedSamples(samples))
+export const absoluteDeviationMean = (samples: Samples, mean: number): number => {
+  let result = 0
+  const len = samples.length
+
+  let i = 0
+
+  while (i < len) {
+    result += (Math.abs(samples[i++]! - mean) - result) / i // eslint-disable-line @typescript-eslint/no-non-null-assertion
+  }
+
+  return result
+}
 
 /**
- * Computes the absolute deviation of a sample given an aggregation.
- * @param samples - the sample
- * @param aggFn - the aggregation function to use
- * @param aggValue - the aggregated value to use
- * @returns the absolute deviation of the sample given the aggregation
+ * Computes the median absolute deviation from the median.
+ * @param samples - the sorted sample
+ * @param median - the median of the sample
+ * @returns the median absolute deviation
  */
-const absoluteDeviation = <S extends Samples = Samples>(
-  samples: S,
-  aggFn: (arr: S) => number,
-  aggValue = aggFn(samples)
-) => {
-  const absoluteDeviations: S = [] as unknown as S
-
-  for (const sample of samples) {
-    absoluteDeviations.push(Math.abs(sample - aggValue))
-  }
+export function absoluteDeviationMedian (samples: SortedSamples, median: number): number {
+  const len = samples.length
+  if (len === 1) return 0
+
+  const mid = len >> 1
+  const halfLen = (len + 1) >> 1
+
+  let low = 0
+  let high = mid
+  let c1, c2, l1, l2, r1, r2
 
-  return aggFn(absoluteDeviations)
+  while (low <= high) {
+    c1 = (low + high) >> 1
+    c2 = halfLen - c1
+
+    // eslint-disable-next-line @typescript-eslint/no-non-null-assertion
+    l1 = c1 === 0 ? -Infinity : median - samples[mid - c1]!
+
+    // eslint-disable-next-line @typescript-eslint/no-non-null-assertion
+    r1 = c1 === mid ? Infinity : median - samples[mid - c1 - 1]!
+
+    // eslint-disable-next-line @typescript-eslint/no-non-null-assertion
+    l2 = c2 === 0 ? -Infinity : samples[mid + c2 - 1]! - median
+    // eslint-disable-next-line @typescript-eslint/no-non-null-assertion
+    r2 = c2 === len - mid ? Infinity : samples[mid + c2]! - median
+
+    if (l1 <= r2 && l2 <= r1) {
+      return len & 1 // check for odd length
+        ? Math.max(l1, l2)
+        : (Math.max(l1, l2) + Math.min(r1, r2)) / 2
+    }
+
+    if (l1 > r2) {
+      high = c1 - 1
+    } else {
+      low = c1 + 1
+    }
+  }
+  return 0 // should never reach here
 }
 
 /**
@@ -392,22 +408,21 @@ const absoluteDeviation = <S extends Samples = Samples>(
  * @returns the statistics of the sample
  */
 export const getStatisticsSorted = (samples: SortedSamples): Statistics => {
-  const mean = average(samples)
-  const vr = variance(samples, mean)
+  const { mean, vr } = meanAndVariance(samples)
   const sd = Math.sqrt(vr)
   const sem = sd / Math.sqrt(samples.length)
   const df = samples.length - 1
   const critical = tTable[df || 1] ?? tTable[0]
   const moe = sem * critical
   const absMean = Math.abs(mean)
   const rme = absMean === 0 ? Infinity : (moe / absMean) * 100
-  const p50 = medianSorted(samples)
+  const p50 = quantileSorted(samples, 0.5)
 
   return {
-    aad: absoluteDeviation(samples, average, mean),
+    aad: absoluteDeviationMean(samples, mean),
     critical,
     df,
-    mad: absoluteDeviation(samples, median, p50),
+    mad: absoluteDeviationMedian(samples, p50),
     // eslint-disable-next-line @typescript-eslint/no-non-null-assertion
     max: samples[df]!,
     mean,

test/utils-absolute-deviation-median.test.ts

@@ -0,0 +1,119 @@
+import { describe, expect, it } from 'vitest'
+
+import { absoluteDeviationMedian, Samples, type SortedSamples } from '../src/utils'
+import { toSortedSamples } from './utils'
+
+// Helper: calculate median of a sorted array
+const medianFn = (samples: SortedSamples): number => {
+  const len = samples.length
+  const mid = len >> 1
+  return len & 1
+    ? samples[mid]! // eslint-disable-line @typescript-eslint/no-non-null-assertion
+    : (samples[mid - 1]! + samples[mid]!) / 2 // eslint-disable-line @typescript-eslint/no-non-null-assertion
+}
+
+// Reference implementation: median of absolute deviations
+const absoluteDeviationMedianTrivial = (samples: SortedSamples): number => {
+  const median = medianFn(samples)
+  const deviations = samples.map(v => Math.abs(v - median)) as Samples
+  return medianFn(toSortedSamples(deviations))
+}
+
+describe('absoluteDeviationMedian()', () => {
+  it('Simple odd length', () => {
+    const samples = toSortedSamples([1, 2, 3, 4, 5, 6, 7, 8, 9])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Simple even length', () => {
+    const samples = toSortedSamples([1, 2, 3, 4, 5, 6, 7, 8])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('With outliers', () => {
+    const samples = toSortedSamples([1, 2, 3, 100, 200])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('All same', () => {
+    const samples = toSortedSamples([5, 5, 5, 5, 5])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Two elements', () => {
+    const samples = toSortedSamples([1, 9])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Single element', () => {
+    const samples = toSortedSamples([42])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Symmetric', () => {
+    const samples = toSortedSamples([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Large spread', () => {
+    const samples = toSortedSamples([1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Duplicates at start', () => {
+    const samples = toSortedSamples([1, 1, 1, 5, 10, 15, 20])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Duplicates at end', () => {
+    const samples = toSortedSamples([1, 5, 10, 15, 20, 20, 20])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Duplicates around median', () => {
+    const samples = toSortedSamples([1, 2, 5, 5, 5, 8, 9])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Many duplicates', () => {
+    const samples = toSortedSamples([1, 2, 2, 3, 3, 3, 4, 4, 5])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Alternating duplicates', () => {
+    const samples = toSortedSamples([1, 1, 2, 2, 3, 3, 4, 4, 5, 5])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Almost all same with outlier', () => {
+    const samples = toSortedSamples([5, 5, 5, 5, 5, 5, 5, 100])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Two values repeated', () => {
+    const samples = toSortedSamples([1, 1, 1, 1, 9, 9, 9, 9])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('Complex duplicates', () => {
+    const samples = toSortedSamples([1, 1, 2, 3, 3, 3, 4, 5, 5, 6, 6, 6, 6, 7])
+    expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(absoluteDeviationMedianTrivial(samples))
+  })
+
+  it('fuzzing test', () => {
+    const rounds = 1000
+    const len = 10
+
+    for (let j = 0; j < rounds; ++j) {
+      const samplesArray: Samples = new Array(len) as unknown as Samples
+      for (let i = 0; i < len; i++) {
+        samplesArray[i] = (Math.random() * 10)
+      }
+
+      const samples = toSortedSamples(samplesArray)
+      expect(absoluteDeviationMedian(samples, medianFn(samples))).toBe(
+        absoluteDeviationMedianTrivial(samples)
+      )
+    }
+  })
+})

test/utils-default-convert-task-result-for-console-table.test.ts

@@ -5,8 +5,8 @@ import type { Task } from '../src'
 import {
   defaultConvertTaskResultForConsoleTable,
   getStatisticsSorted,
-  toSortedSamples
 } from '../src/utils'
+import { toSortedSamples } from './utils'
 
 test('defaultConvertTaskResultForConsoleTable - not-started', () => {
   expect(defaultConvertTaskResultForConsoleTable({