P3: normalization fit (1.2.0 and 1.1.2), apply, polar, simurg normalize

- signal/normalization/numerics.py: vendor percentile, trimmean, MAD sigma
  (old/new), bisquare robust fits, float32/float64 typing and sequential
  sums reproduced exactly (ported from gtc2vcf, MIT)
- signal/normalization/fit.py: outlier removal, offset, rotation, shear,
  scale, nn12/rr12, single-channel model for small Infinium I bins,
  exact and binned nearest-neighbour sweeps; vendor quirks marked QUIRK
- signal/normalization/apply.py: affine transform and theta/R
- pipeline/normalize.py + CLI `simurg normalize`
- oracle GTCs rebuilt with -ffp-contract=off (aarch64 gcc fuses a*b+c by
  default, shifting the last bits); the old FMA oracle differed by <=4 ULP

Result: transforms bit-identical to the oracle for all 36 GSA samples x
48 bins, for both normalization versions (3456/3456 bins, 0 ULP).

src/simurg/signal/normalization/numerics.py

@@ -0,0 +1,205 @@
+"""Numerics with the exact rounding behaviour of the reference normalization (plan Appendix C).
+
+Ported from gtc2vcf ``idat2gtc.c`` (MIT, © Giulio Genovese), which reproduces the vendor's MATLAB-derived routines.
+Rules that matter for bit-exactness:
+
+* values that are ``float`` in the reference are ``np.float32`` here; ``double`` accumulators are ``np.float64``;
+* every running sum is **sequential** (``np.cumsum``), never pairwise (``np.sum``) or compensated (``sum``);
+* the reference is compiled without fused multiply-add, so every product is rounded before the next addition.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Callable
+
+import numpy as np
+import numpy.typing as npt
+
+F32 = np.float32
+F64 = np.float64
+FloatArr = npt.NDArray[np.float32]
+DblArr = npt.NDArray[np.float64]
+TINY = F32(np.finfo(np.float32).tiny) * F32(np.finfo(np.float32).eps)  # FLT_MIN * FLT_EPSILON (denormal)
+ADJ0 = 100.0 + 24832 * float(np.finfo(np.float64).eps)  # robustfit0 leverage factor as used by the vendor
+
+
+def seqsum(a: npt.NDArray[np.float64]) -> float:
+    """Left-to-right double-precision sum (C loop ``s += a[i]``)."""
+    return float(np.cumsum(a, dtype=np.float64)[-1]) if len(a) else 0.0
+
+
+def percentile(sorted_vals: FloatArr, p: int) -> np.float32:
+    """Vendor percentile on an array assumed sorted (the reference also calls it on unsorted data — see callers)."""
+    n = len(sorted_vals)
+    if n == 0:
+        return F32(np.nan)
+    i1 = (n * p) // 100
+    f = F32(n * p) / F32(100.0) - F32(i1)
+    if f < F32(0.5):
+        i1 -= 1
+    if i1 < 0:
+        return F32(sorted_vals[0])
+    if i1 >= n - 1:
+        return F32(sorted_vals[n - 1])
+    x1 = F32(100.0) * (F32(i1) + F32(0.5)) / F32(n)
+    x2 = F32(100.0) * (F32(i1 + 1) + F32(0.5)) / F32(n)
+    y1 = F32(sorted_vals[i1])
+    y2 = F32(sorted_vals[i1 + 1])
+    m = (y2 - y1) / (x2 - x1)
+    return F32(y1 + m * (F32(p) - x1))
+
+
+def trimmean(vals: FloatArr, percent: int) -> np.float32:
+    """MATLAB ``trimmean`` with the vendor percentile bounds; mean in double, final division in float."""
+    s = np.sort(np.asarray(vals, dtype=np.float32))
+    high = percentile(s, 100 - percent // 2)
+    low = percentile(s, percent // 2)
+    keep = s[(s >= low) & (s <= high)]
+    total = seqsum(keep.astype(np.float64))
+    with np.errstate(invalid="ignore", divide="ignore"):
+        return F32(F32(total) / F32(len(keep)))
+
+
+def iqr(vals: FloatArr) -> np.float32:
+    s = np.sort(np.asarray(vals, dtype=np.float32))
+    return F32(percentile(s, 75) - percentile(s, 25))
+
+
+def median_sorted(s: FloatArr) -> np.float32:
+    n = len(s)
+    if n == 0:
+        return F32(0.0)
+    if n % 2 == 1:
+        return F32(s[n // 2])
+    return F32((s[n // 2 - 1] + s[n // 2]) * F32(0.5))
+
+
+def mean_f(vals: FloatArr) -> np.float32:
+    if len(vals) == 0:
+        return F32(np.nan)
+    return F32(seqsum(np.asarray(vals, dtype=np.float64)) / float(len(vals)))
+
+
+def nanmean_f(vals: FloatArr) -> np.float32:
+    v = np.asarray(vals, dtype=np.float32)
+    v = v[~np.isnan(v)]
+    if len(vals) == 0:
+        return F32(np.nan)
+    with np.errstate(invalid="ignore", divide="ignore"):
+        return F32(seqsum(v.astype(np.float64)) / float(len(v)))
+
+
+def madsigma_new(r: DblArr, p: int) -> float:
+    """MATLAB ``madsigma`` (MAD from 0) — normalization 1.2.0."""
+    rs = np.sort(np.abs(r).astype(np.float32))
+    s = float(median_sorted(rs[p - 1 :])) / 0.6745
+    if s == 0.0:
+        s = 0.5 * float(mean_f(rs))
+    return s
+
+
+def madsigma_old(r: DblArr, p: int) -> float:
+    """MATLAB ``madsigma`` before 2002 (MAD from the median) — normalization 1.1.2."""
+    rs = np.sort(r.astype(np.float32))
+    m = median_sorted(rs)
+    rs = np.sort(np.abs(rs - m).astype(np.float32))
+    if abs(m) > rs[-1]:
+        rs = np.sort(np.abs(r.astype(np.float32)))
+    s = float(median_sorted(rs[p - 1 :])) / 0.6745
+    if s == 0.0:
+        s = 0.5 * float(mean_f(rs))
+    return s
+
+
+MadSigma = Callable[[DblArr, int], float]
+
+
+def _converged(new: float, old: float) -> bool:
+    return abs(new - old) <= 1e-6 * float(max(F32(abs(new)), F32(abs(old))))
+
+
+def robustfit0(x: FloatArr, y: FloatArr, madsigma: MadSigma) -> np.float32:
+    """``robustfit(x, y, 'bisquare', 4.685, 'off')`` — slope through the origin, vendor leverage quirk."""
+    xd = x.astype(np.float64)
+    yd = y.astype(np.float64)
+    sumx2 = seqsum(xd * xd)
+    if sumx2 == 0:
+        raise ValueError("linsolve0: singular")
+    m = seqsum(xd * yd) / sumx2
+    it = 0
+    while True:
+        r = (yd - m * xd) * ADJ0
+        s = madsigma(r, 1)
+        if s == 0.0:
+            s = 1.0
+        r = r * (1.0 / (s * 4.685))
+        t = 1.0 - r * r
+        w = np.where(np.abs(r) < 1, t * t, 0.0)
+        m0 = m
+        wsumx2 = seqsum(w * (xd * xd))
+        if wsumx2 == 0:
+            raise ValueError("wfit0: singular")
+        m = seqsum(w * xd * yd) / wsumx2
+        it += 1
+        if it >= 50 or _converged(m, m0):
+            break
+    return F32(m)
+
+
+def robustfit1(x: FloatArr, y: FloatArr, madsigma: MadSigma) -> tuple[np.float32, np.float32]:
+    """``robustfit(x, y, 'bisquare', 4.685, 'on')`` → (intercept, slope)."""
+    n = len(x)
+    xd = x.astype(np.float64)
+    yd = y.astype(np.float64)
+    sumx2 = seqsum(xd * xd)
+    sumxy = seqsum(yd * xd)
+    sumx = seqsum(xd)
+    sumy = seqsum(yd)
+    denom = float(n) * sumx2 - sumx * sumx
+    if denom == 0:
+        raise ValueError("linsolve1: singular (wrong manifest?)")
+    m = (n * sumxy - sumx * sumy) / denom
+    b = (sumy * sumx2 - sumx * sumxy) / denom
+    denom2 = float(n) * sumx2 - sumx * sumx
+    h = np.minimum(0.9999, (float(n) * (xd * xd) - 2.0 * sumx * xd + sumx2) / denom2)
+    adj = 1.0 / np.sqrt(1.0 - h)
+    it = 0
+    while True:
+        r = (yd - b - m * xd) * adj
+        s = madsigma(r, 2)
+        if s == 0.0:
+            s = 1.0
+        r = r * (1.0 / (s * 4.685))
+        t = 1.0 - r * r
+        w = np.where(np.abs(r) < 1, t * t, 0.0)
+        b0, m0 = b, m
+        wsumx2 = seqsum(w * (xd * xd))
+        wsumxy = seqsum(w * xd * yd)
+        wsumx = seqsum(w * xd)
+        wsumy = seqsum(w * yd)
+        wsum = seqsum(w)
+        d = wsum * wsumx2 - wsumx * wsumx
+        if d == 0:
+            raise ValueError("wfit1: singular")
+        m = (wsum * wsumxy - wsumx * wsumy) / d
+        b = (wsumy * wsumx2 - wsumx * wsumxy) / d
+        it += 1
+        if it >= 50 or (_converged(b, b0) and _converged(m, m0)):
+            break
+    return F32(b), F32(m)
+
+
+def linspace_f(n: int, lo: np.float32, hi: np.float32) -> FloatArr:
+    """``lo + (hi - lo) * i / (n - 1)`` evaluated in float32, left to right."""
+    i = np.arange(n, dtype=np.float32)
+    return (F32(lo) + ((F32(hi) - F32(lo)) * i) / F32(n - 1)).astype(np.float32)
+
+
+def unique_indices(idx: npt.NDArray[np.int64]) -> npt.NDArray[np.int64]:
+    """``matlab_unique`` on nearest-neighbour results: sorted, deduplicated, with -1 (no neighbour) dropped."""
+    return np.unique(idx[idx >= 0])
+
+
+def trunc_int(v: FloatArr) -> npt.NDArray[np.int64]:
+    """C ``(int)`` cast of float32 values (truncation toward zero); inputs are finite here."""
+    return np.trunc(v).astype(np.int64)