Spot an anomaly in a metric

"""
Flag anomalies in a metric with a robust threshold on a sliding window.

Rung N0. Median and median absolute deviation, standard library only, and a
verdict that can be read at three in the morning.

The mean and the standard deviation are the wrong tools here. One incident
drags both, so a large enough spike widens the very band that was supposed to
catch it. The median and the median absolute deviation ignore up to half the
window, which is exactly the property an alert needs.

Every verdict carries the numbers it was made of. "Anomaly at 03:12" leaves
whoever was woken up to reconstruct the reasoning before they can act;
"measured 4800, usual 1200, allowed up to 2100" is already half the
diagnosis, and it is the same three numbers the threshold itself used.
"""

from dataclasses import dataclass
from statistics import median

# Scaling that puts the median absolute deviation on the same footing as a
# standard deviation for normally distributed data, so that a threshold of
# 3.5 keeps the meaning it has everywhere else.
NORMAL_SCALE = 1.4826


@dataclass
class Verdict:
    """One judged point, and the whole reasoning behind the judgement."""

    index: int
    value: float
    usual: float  # median of the window that came before
    deviation: float  # how far the point sits from that median
    limit: float  # how far it was allowed to sit
    is_anomaly: bool


def scan(series: list[float], window: int = 24, threshold: float = 3.5) -> list[Verdict]:
    """
    Judge every point against the `window` points that precede it.

    The first `window` points get no verdict at all. A point cannot be
    compared with a history that does not exist yet, and saying nothing is
    more honest than comparing it with a shorter, noisier window.
    """
    verdicts = []
    for index in range(window, len(series)):
        reference = series[index - window : index]
        usual = median(reference)
        spread = NORMAL_SCALE * median([abs(value - usual) for value in reference])
        deviation = abs(series[index] - usual)
        limit = threshold * spread
        verdicts.append(Verdict(index, series[index], usual, deviation, limit, deviation > limit))
    return verdicts


def anomalies(series: list[float], window: int = 24, threshold: float = 3.5) -> list[Verdict]:
    """The subset a pager should see, each one still carrying its numbers."""
    return [verdict for verdict in scan(series, window, threshold) if verdict.is_anomaly]

/**
 * Flag anomalies in a metric with a robust threshold on a sliding window.
 *
 * Rung N0. Median and median absolute deviation, no dependency, and a verdict
 * that can be read at three in the morning.
 *
 * The mean and the standard deviation are the wrong tools here. One incident
 * drags both, so a large enough spike widens the very band that was supposed
 * to catch it. The median and the median absolute deviation ignore up to half
 * the window, which is exactly the property an alert needs.
 *
 * Every verdict carries the numbers it was made of. "Anomaly at 03:12" leaves
 * whoever was woken up to reconstruct the reasoning before they can act;
 * "measured 4800, usual 1200, allowed up to 2100" is already half the
 * diagnosis, and it is the same three numbers the threshold itself used.
 */

// Scaling that puts the median absolute deviation on the same footing as a
// standard deviation for normally distributed data, so that a threshold of
// 3.5 keeps the meaning it has everywhere else.
const NORMAL_SCALE = 1.4826;

/** Middle value, or the average of the two middle ones on an even count. */
function median(values) {
  const sorted = [...values].sort((a, b) => a - b);
  const middle = Math.floor(sorted.length / 2);
  return sorted.length % 2 ? sorted[middle] : (sorted[middle - 1] + sorted[middle]) / 2;
}

/**
 * Judge every point against the `window` points that precede it.
 *
 * The first `window` points get no verdict at all. A point cannot be compared
 * with a history that does not exist yet, and saying nothing is more honest
 * than comparing it with a shorter, noisier window.
 *
 * Each verdict is the whole reasoning: the value, what the window called
 * usual, how far the value sat from it, and how far it was allowed to sit.
 */
export function scan(series, { window = 24, threshold = 3.5 } = {}) {
  const verdicts = [];
  for (let index = window; index < series.length; index += 1) {
    const reference = series.slice(index - window, index);
    const usual = median(reference);
    const spread = NORMAL_SCALE * median(reference.map((value) => Math.abs(value - usual)));
    const deviation = Math.abs(series[index] - usual);
    const limit = threshold * spread;
    verdicts.push({
      index,
      value: series[index],
      usual,
      deviation,
      limit,
      isAnomaly: deviation > limit,
    });
  }
  return verdicts;
}

/** The subset a pager should see, each one still carrying its numbers. */
export function anomalies(series, options = {}) {
  return scan(series, options).filter((verdict) => verdict.isAnomaly);
}

"""
Spot anomalies in several metrics at once with an isolation forest.

Rung N1. The robust threshold of N0 watches one metric at a time, and some
incidents are invisible that way: every metric stays inside its usual range,
and only the combination is impossible. Errors at their normal ceiling while
traffic sits at its normal middle is one such minute, and no single-series
threshold will ever ring for it.

An isolation forest cuts the space at random and measures how few cuts it
takes to leave a point on its own. A point in the middle of the crowd needs
many; a point out on its own needs three or four. There is nothing to label,
and the only real knob is the size of the forest.

What this rung costs is the thing N0 was good at. The output is a rank
between zero and one, not a measured gap against a stated limit. Whoever is
woken up is told that the minute was unusual, not in what way.
"""

from sklearn.ensemble import IsolationForest

# Trees are grown on random cuts, so the seed is part of the contract: an
# alert that changes between two runs on the same data is not an alert.
SEED = 0


def train(rows: list[list[float]], trees: int = 100, seed: int = SEED) -> IsolationForest:
    """
    `rows` is one observation per minute, the metrics always in the same order.

    Nothing is labelled and nothing is scaled: the cuts are drawn between the
    smallest and the largest value of each metric, so a metric counted in
    milliseconds and one counted in requests weigh the same.
    """
    model = IsolationForest(n_estimators=trees, random_state=seed)
    return model.fit(rows)


def score(model: IsolationForest, row: list[float]) -> float:
    """
    Between zero and one. Above one half, the point took fewer cuts to isolate
    than the crowd did, which is the whole definition of an anomaly here.
    """
    return float(-model.score_samples([row])[0])


def anomalies(model: IsolationForest, rows: list[list[float]], threshold: float = 0.6) -> list[int]:
    """
    The indices worth looking at, and the threshold is yours to set.

    Move it towards one to be woken up less often and miss more; towards zero
    for the opposite trade. There is no value of it that turns the score into
    an explanation.
    """
    return [index for index, row in enumerate(rows) if score(model, row) > threshold]

/**
 * Spot anomalies in several metrics at once with an isolation forest.
 *
 * Rung N1. The robust threshold of N0 watches one metric at a time, and some
 * incidents are invisible that way: every metric stays inside its usual
 * range, and only the combination is impossible. Night-time traffic with
 * daytime errors is one such minute, and no single-series threshold will ever
 * ring for it.
 *
 * The forest cuts the space at random and measures how few cuts it takes to
 * leave a point on its own. A point in the middle of the crowd needs many; a
 * point out on its own needs three or four. There is nothing to label, and
 * the only real knob is the size of the forest.
 *
 * Written out in full rather than pulled from a library, because that is the
 * whole argument of this rung: the classical tool is small enough to read.
 * What it costs is the thing N0 was good at — the output is a rank between
 * zero and one, not a measured gap against a stated limit.
 */

/** Deterministic pseudo-random numbers: an alert that changes between two runs is not an alert. */
function generator(seed) {
  let state = seed >>> 0;
  return () => {
    state = (state + 0x6d2b79f5) >>> 0;
    let t = Math.imul(state ^ (state >>> 15), 1 | state);
    t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}

/** Average depth of an unsuccessful search in a binary tree of n points. */
function averageDepth(n) {
  if (n <= 1) return 0;
  const EULER = 0.5772156649015329;
  return 2 * (Math.log(n - 1) + EULER) - (2 * (n - 1)) / n;
}

/** One tree: cut on a random metric at a random value until the points are alone. */
function grow(rows, random, depth, maxDepth) {
  if (depth >= maxDepth || rows.length <= 1) return { size: rows.length };
  const metric = Math.floor(random() * rows[0].length);
  const values = rows.map((row) => row[metric]);
  const low = Math.min(...values);
  const high = Math.max(...values);
  if (low === high) return { size: rows.length };
  const cut = low + random() * (high - low);
  return {
    metric,
    cut,
    below: grow(rows.filter((row) => row[metric] < cut), random, depth + 1, maxDepth),
    above: grow(rows.filter((row) => row[metric] >= cut), random, depth + 1, maxDepth),
  };
}

/** How deep this row falls, plus what an unfinished leaf would still have cost. */
function depthOf(node, row, depth) {
  if (node.below === undefined) return depth + averageDepth(node.size);
  return depthOf(row[node.metric] < node.cut ? node.below : node.above, row, depth + 1);
}

/** A subsample without replacement, so every tree sees a different crowd. */
function subsample(rows, size, random) {
  const shuffled = [...rows];
  for (let i = 0; i < size; i += 1) {
    const j = i + Math.floor(random() * (shuffled.length - i));
    [shuffled[i], shuffled[j]] = [shuffled[j], shuffled[i]];
  }
  return shuffled.slice(0, size);
}

/**
 * `rows` is one observation per minute, the metrics always in the same order.
 *
 * Nothing is labelled and nothing is scaled: the cuts are drawn between the
 * smallest and the largest value of each metric, so a metric counted in
 * milliseconds and one counted in requests weigh the same.
 */
export function train(rows, { trees = 100, seed = 0 } = {}) {
  const random = generator(seed);
  const size = Math.min(256, rows.length);
  const maxDepth = Math.ceil(Math.log2(size));
  const forest = [];
  for (let i = 0; i < trees; i += 1) {
    forest.push(grow(subsample(rows, size, random), random, 0, maxDepth));
  }
  return { forest, normaliser: averageDepth(size) };
}

/**
 * Between zero and one. Above one half, the point took fewer cuts to isolate
 * than the crowd did, which is the whole definition of an anomaly here.
 */
export function score(model, row) {
  const total = model.forest.reduce((sum, tree) => sum + depthOf(tree, row, 0), 0);
  return 2 ** (-total / model.forest.length / model.normaliser);
}

/**
 * The indices worth looking at, and the threshold is yours to set.
 *
 * Move it towards one to be woken up less often and miss more; towards zero
 * for the opposite trade. There is no value of it that turns the score into
 * an explanation.
 */
export function anomalies(model, rows, threshold = 0.6) {
  return rows.map((_, index) => index).filter((index) => score(model, rows[index]) > threshold);
}