monitor.anomaly

import { type Anomaly, type AnomalyMethod } from "./monitor_anomaly_types.ts";

function isqrt(n: bigint): bigint {
  if (n < 2n) return n;
  let x = n;
  let y = (x + 1n) / 2n;
  while (y < x) {
    x = y;
    y = (x + n / x) / 2n;
  }
  return x;
}

const abs = (n: bigint): bigint => (n < 0n ? -n : n);

// Half-up with halves away from zero, as math.round-div rounds.
function halfUp(n: bigint, d: bigint): bigint {
  const q = (2n * abs(n) + d) / (2n * d);
  return n < 0n ? -q : q;
}

function medianTwice(sorted: readonly bigint[]): bigint {
  const m = Math.floor(sorted.length / 2);
  return sorted.length % 2 === 1 ? 2n * sorted[m] : sorted[m - 1] + sorted[m];
}

const byValue = (a: bigint, b: bigint): number => (a < b ? -1 : a > b ? 1 : 0);

/**
 * Whether `value` is an outlier against `history`, by z-score (`stddev`) or by
 * Iglewicz and Hoaglin's modified z-score (`mad`). Exact integer arithmetic:
 * the threshold comparison never depends on floating-point rounding.
 */
export function detectAnomaly(history: readonly number[], value: number, method: AnomalyMethod, thresholdHundredths: number): Anomaly {
  if (history.length < 2) throw new RangeError(`history must have at least 2 values, received ${history.length}`);
  if (thresholdHundredths < 1) throw new RangeError(`thresholdHundredths must be at least 1, received ${thresholdHundredths}`);
  const xs = history.map((x) => BigInt(x));
  const v = BigInt(value);
  const t = BigInt(thresholdHundredths);
  const n = BigInt(xs.length);
  let anomaly: boolean;
  let signed: bigint;
  let center: bigint;
  let spread: bigint;
  let score: bigint | null;
  if (method === "stddev") {
    let s = 0n;
    let ss = 0n;
    for (const x of xs) {
      s += x;
      ss += x * x;
    }
    // n^2 * variance and n * (value - mean), both whole numbers.
    const q = n * ss - s * s;
    signed = n * v - s;
    center = halfUp(s * 1000n, n);
    spread = isqrt(1000000n * q) / n;
    if (q === 0n) {
      anomaly = signed !== 0n;
      score = null;
    } else {
      anomaly = 10000n * signed * signed > t * t * q;
      score = isqrt((10000n * signed * signed) / q);
    }
  } else if (method === "mad") {
    // Twice the median and four times the MAD are whole numbers.
    const med2 = medianTwice([...xs].sort(byValue));
    const mad4 = medianTwice(xs.map((x) => abs(2n * x - med2)).sort(byValue));
    signed = 2n * v - med2;
    center = med2 * 500n;
    spread = mad4 * 250n;
    if (mad4 === 0n) {
      anomaly = signed !== 0n;
      score = null;
    } else {
      // |M| = 0.6745 * |2v - med2| / 2 / (mad4 / 4) = 1349 * |signed| / (1000 * mad4)
      anomaly = 1349n * abs(signed) * 100n > t * 1000n * mad4;
      score = (1349n * abs(signed) * 100n) / (1000n * mad4);
    }
  } else {
    throw new RangeError(`unknown anomaly method: ${method as string}`);
  }
  return {
    anomaly,
    direction: !anomaly ? "normal" : signed > 0n ? "high" : "low",
    centerMilli: Number(center),
    spreadMilli: Number(spread),
    scoreHundredths: score === null ? null : Number(score),
  };
}