anomaly_scores.py

#!/usr/bin/env python3
"""Score rows of a numeric CSV with four unsupervised anomaly detectors and a rank-averaged ensemble.

Detectors (standard library only):
  robust_z     max over features of |x - median| / (1.4826 * MAD): one feature far from normal
  mahalanobis  sqrt((x - mu)' S^-1 (x - mu)) with a covariance fitted on the central 90% of rows by robust z
               (resists masking by the outliers themselves): unusual combinations of features
  knn          distance to the k-th nearest neighbour (standardised features): isolated points
  iforest      isolation forest (Liu et al. 2008): average path length over random trees on subsamples
               of 256 rows, score = 2^(-E[h]/c(256)); short paths = easy to isolate = anomalous
  ensemble     mean of the four rank-normalised scores (robust to the detectors' different scales)
For the top rows it shows which features deviate most (robust z per feature).

Fit and score can differ: --train defines "normal" (e.g. last month), --score is scored against it.

Usage:
    python anomaly_scores.py --train sensors.csv --drop reading_id is_anomaly --top 10
    python anomaly_scores.py --train ref.csv --score new.csv --drop id --out scores.csv
    python anomaly_scores.py --train data.csv --columns temperature pressure --k 10 --trees 200 --json

Exit codes: 0 scored, 2 usage or read error.
"""

import argparse
import csv
import json
import math
import random
import sys


def num(v):
    try:
        x = float(v)
        return x if math.isfinite(x) else None
    except (TypeError, ValueError):
        return None


def load(path, cols, drop):
    with open(path, newline="", encoding="utf-8-sig") as fh:
        rd = csv.DictReader(fh)
        header = rd.fieldnames or []
        rows = list(rd)
    cols = cols or [c for c in header if c not in drop and all(num(r[c]) is not None for r in rows if (r[c] or "").strip())]
    keep = [r for r in rows if all(num(r.get(c)) is not None for c in cols)]
    return cols, keep, [[num(r[c]) for c in cols] for r in keep]


def median(xs):
    s = sorted(xs)
    n = len(s)
    return (s[n // 2] + s[(n - 1) // 2]) / 2


def invert(M):
    n = len(M)
    A = [row[:] + [1.0 if i == j else 0.0 for j in range(n)] for i, row in enumerate(M)]
    for c in range(n):
        p = max(range(c, n), key=lambda r: abs(A[r][c]))
        if abs(A[p][c]) < 1e-12:
            A[p][c] += 1e-6                     # ridge for singular covariance
        A[c], A[p] = A[p], A[c]
        piv = A[c][c]
        A[c] = [v / piv for v in A[c]]
        for r in range(n):
            if r != c:
                f = A[r][c]
                A[r] = [v - f * w for v, w in zip(A[r], A[c])]
    return [row[n:] for row in A]


def c_factor(n):
    return 2 * (math.log(n - 1) + 0.5772156649) - 2 * (n - 1) / n if n > 2 else (1.0 if n == 2 else 0.0)


def build_tree(X, depth, limit, rng):
    if depth >= limit or len(X) <= 1:
        return ("leaf", len(X))
    d = len(X[0])
    f = rng.randrange(d)
    lo, hi = min(x[f] for x in X), max(x[f] for x in X)
    if lo == hi:
        return ("leaf", len(X))
    t = rng.uniform(lo, hi)
    left = [x for x in X if x[f] < t]
    right = [x for x in X if x[f] >= t]
    return ("node", f, t, build_tree(left, depth + 1, limit, rng), build_tree(right, depth + 1, limit, rng))


def path_len(tree, x, depth=0):
    if tree[0] == "leaf":
        return depth + c_factor(tree[1])
    _, f, t, l, r = tree
    return path_len(l if x[f] < t else r, x, depth + 1)


def ranks01(v):
    order = sorted(range(len(v)), key=v.__getitem__)
    r = [0.0] * len(v)
    for pos, i in enumerate(order):
        r[i] = pos / max(1, len(v) - 1)
    return r


def main():
    ap = argparse.ArgumentParser(description="Unsupervised anomaly scores (robust z, Mahalanobis, kNN, iForest).",
                                 epilog=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--train", required=True)
    ap.add_argument("--score")
    ap.add_argument("--columns", nargs="*")
    ap.add_argument("--drop", nargs="*", default=[])
    ap.add_argument("--k", type=int, default=5)
    ap.add_argument("--trees", type=int, default=100)
    ap.add_argument("--top", type=int, default=10)
    ap.add_argument("--out", help="write all scores to this CSV")
    ap.add_argument("--seed", type=int, default=0)
    ap.add_argument("--json", action="store_true")
    a = ap.parse_args()
    try:
        cols, _, Xtr = load(a.train, a.columns, set(a.drop))
        if a.score:
            _, srows, Xsc = load(a.score, cols, set())
        else:
            with open(a.train, newline="", encoding="utf-8-sig") as fh:
                srows = [r for r in csv.DictReader(fh) if all(num(r.get(c)) is not None for c in cols)]
            Xsc = Xtr
    except (OSError, KeyError) as e:
        print(f"error: {e}", file=sys.stderr)
        return 2
    if not cols or len(Xtr) < 20:
        print("error: need >= 20 complete numeric training rows", file=sys.stderr)
        return 2
    d = len(cols)
    rng = random.Random(a.seed)
    med = [median([x[j] for x in Xtr]) for j in range(d)]
    mad = [median([abs(x[j] - med[j]) for x in Xtr]) * 1.4826 or 1e-9 for j in range(d)]
    rz = lambda x: [(x[j] - med[j]) / mad[j] for j in range(d)]
    robust = [max(abs(v) for v in rz(x)) for x in Xsc]
    # covariance from the central 90% (by robust z) to reduce masking
    rz_tr = sorted(range(len(Xtr)), key=lambda i: max(abs(v) for v in rz(Xtr[i])))
    core = [Xtr[i] for i in rz_tr[:max(d + 2, int(0.9 * len(Xtr)))]]
    mu = [sum(x[j] for x in core) / len(core) for j in range(d)]
    S = [[sum((x[i] - mu[i]) * (x[j] - mu[j]) for x in core) / (len(core) - 1) for j in range(d)] for i in range(d)]
    Si = invert(S)
    maha = []
    for x in Xsc:
        dv = [x[j] - mu[j] for j in range(d)]
        maha.append(math.sqrt(max(0.0, sum(dv[i] * Si[i][j] * dv[j] for i in range(d) for j in range(d)))))
    sd = [math.sqrt(S[j][j]) or 1.0 for j in range(d)]
    Ztr = [[(x[j] - mu[j]) / sd[j] for j in range(d)] for x in Xtr]
    Zsc = [[(x[j] - mu[j]) / sd[j] for j in range(d)] for x in Xsc]
    same = a.score is None
    knn = []
    for i, z in enumerate(Zsc):
        ds = sorted(sum((p - q) ** 2 for p, q in zip(z, t)) for jj, t in enumerate(Ztr) if not (same and jj == i))
        knn.append(math.sqrt(ds[min(a.k, len(ds)) - 1]))
    psi = min(256, len(Xtr))
    limit = math.ceil(math.log2(psi))
    trees = [build_tree(rng.sample(Xtr, psi), 0, limit, rng) for _ in range(a.trees)]
    cn = c_factor(psi)
    iforest = [2 ** (-(sum(path_len(t, x) for t in trees) / len(trees)) / cn) for x in Xsc]
    ens = [sum(v) / 4 for v in zip(ranks01(robust), ranks01(maha), ranks01(knn), ranks01(iforest))]
    order = sorted(range(len(Xsc)), key=lambda i: -ens[i])
    top = []
    for i in order[:a.top]:
        dev = sorted(zip(cols, rz(Xsc[i])), key=lambda cv: -abs(cv[1]))[:3]
        top.append({"line": i + 2, "ensemble": round(ens[i], 4), "robust_z": round(robust[i], 2),
                    "mahalanobis": round(maha[i], 2), "knn": round(knn[i], 3), "iforest": round(iforest[i], 3),
                    "deviations": [{"feature": c, "robust_z": round(v, 2)} for c, v in dev]})
    if a.out:
        with open(a.out, "w", newline="", encoding="utf-8") as fh:
            keys = list(srows[0].keys()) if srows else []
            w = csv.writer(fh)
            w.writerow(keys + ["score_robust_z", "score_mahalanobis", "score_knn", "score_iforest", "score_ensemble"])
            for r, s1, s2, s3, s4, s5 in zip(srows, robust, maha, knn, iforest, ens):
                w.writerow([r[k] for k in keys] + [round(s1, 4), round(s2, 4), round(s3, 4), round(s4, 4), round(s5, 4)])
    out = {"train_rows": len(Xtr), "scored_rows": len(Xsc), "columns": cols, "top": top}
    if a.json:
        print(json.dumps(out, indent=2))
        return 0
    print(f"Fitted on {len(Xtr)} rows x {d} features; scored {len(Xsc)} rows (k={a.k}, {a.trees} trees, psi={psi})")
    print(f"{'line':>6} {'ensemble':>8} {'robust_z':>8} {'mahal':>7} {'knn':>7} {'iforest':>7}  most deviating features")
    for t in top:
        devs = ", ".join(f"{d_['feature']} {d_['robust_z']:+.1f}" for d_ in t["deviations"])
        print(f"{t['line']:>6} {t['ensemble']:>8.3f} {t['robust_z']:>8.2f} {t['mahalanobis']:>7.2f} {t['knn']:>7.3f} "
              f"{t['iforest']:>7.3f}  {devs}")
    if a.out:
        print(f"\nwrote all scores to {a.out} (evaluate with anomaly_eval.py)")
    print("High Mahalanobis with modest robust z = unusual combination; high robust z = one extreme feature.")
    return 0


if __name__ == "__main__":
    sys.exit(main())