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).

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv

src/simurg/cli/main.py

@@ -1,6 +1,7 @@
 """``simurg`` command-line interface (spec §2, plan §14).
 
-Available: ``version``, ``inspect`` (P1), ``ingest`` (P2). Later phases add the other commands of CLI-001.
+Available: ``version``, ``inspect`` (P1), ``ingest`` (P2), ``normalize`` (P3).
+Later phases add the other commands of CLI-001.
 """
 
 from __future__ import annotations
@@ -252,5 +253,26 @@ def ingest_cmd(
     run_guarded(rep, lambda: run_ingest(opts, rep))
 
 
+# --------------------------------------------------------------------------------------------------------
+@app.command("normalize")
+def normalize_cmd(
+    indir: Annotated[Path, typer.Option("--in", help="Output folder of `simurg ingest`.")],
+    norm_version: Annotated[str, typer.Option("--norm-version", help="1.2.0 (current) or 1.1.2 (older)")] = "1.2.0",
+    threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0,
+    force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False,
+    log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False,
+    debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False,
+    quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False,
+) -> None:
+    """Fit per-bin normalization transforms and write normalized X/Y, theta and R (phase P3)."""
+    from simurg.pipeline.normalize import NormalizeOptions, run_normalize
+
+    rep = Reporter(
+        command="normalize", level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet
+    )
+    opts = NormalizeOptions(indir=indir, version=norm_version, threads=threads, force=force, argv=sys.argv[1:])
+    run_guarded(rep, lambda: run_normalize(opts, rep))
+
+
 if __name__ == "__main__":  # pragma: no cover
     app()

src/simurg/core/schemas.py

@@ -20,6 +20,8 @@ from simurg.formats.bpm import Manifest
 LOCI_SCHEMA_ID = "simurg.loci/1"
 INTENSITY_SCHEMA_ID = "simurg.sample_intensity/1"
 CONTROLS_SCHEMA_ID = "simurg.sample_controls/1"
+NORMALIZED_SCHEMA_ID = "simurg.sample_normalized/1"
+TRANSFORMS_SCHEMA_ID = "simurg.sample_transforms/1"
 
 LOCI_SCHEMA = pa.schema(
     [
@@ -56,6 +58,33 @@ CONTROLS_SCHEMA = pa.schema(
     metadata={"simurg.schema": CONTROLS_SCHEMA_ID},
 )
 
+NORMALIZED_SCHEMA = pa.schema(
+    [("norm_x", pa.float32()), ("norm_y", pa.float32()), ("theta", pa.float32()), ("r", pa.float32())],
+    metadata={"simurg.schema": NORMALIZED_SCHEMA_ID},
+)
+
+TRANSFORMS_SCHEMA = pa.schema(
+    [
+        ("bin_index", pa.int32()),
+        ("norm_bin_id", pa.uint8()),
+        ("n_loci", pa.int32()),
+        ("version", pa.int32()),
+        ("offset_x", pa.float32()),
+        ("offset_y", pa.float32()),
+        ("scale_x", pa.float32()),
+        ("scale_y", pa.float32()),
+        ("shear", pa.float32()),
+        ("theta", pa.float32()),
+        ("cvx", pa.float32()),
+        ("cvy", pa.float32()),
+        ("nn12", pa.float32()),
+        ("rr12", pa.float32()),
+        ("taa", pa.float32()),
+        ("tbb", pa.float32()),
+    ],
+    metadata={"simurg.schema": TRANSFORMS_SCHEMA_ID},
+)
+
 _WRITE_OPTS: dict[str, Any] = {"compression": "zstd", "use_dictionary": True, "write_statistics": True}
 
 
@@ -122,6 +151,26 @@ def controls_table(addresses: np.ndarray, red: np.ndarray, green: np.ndarray, me
     return with_metadata(t, meta)
 
 
+def normalized_table(
+    nx: np.ndarray, ny: np.ndarray, theta: np.ndarray, r: np.ndarray, meta: dict[str, Any]
+) -> pa.Table:
+    t = pa.table(
+        {
+            "norm_x": pa.array(nx, pa.float32()),
+            "norm_y": pa.array(ny, pa.float32()),
+            "theta": pa.array(theta, pa.float32()),
+            "r": pa.array(r, pa.float32()),
+        },
+        schema=NORMALIZED_SCHEMA,
+    )
+    return with_metadata(t, meta)
+
+
+def transforms_table(rows: list[dict[str, Any]], meta: dict[str, Any]) -> pa.Table:
+    t = pa.Table.from_pylist(rows, schema=TRANSFORMS_SCHEMA)
+    return with_metadata(t, meta)
+
+
 def write_parquet(table: pa.Table, path: Path) -> None:
     """Atomic, deterministic Parquet write."""
     with atomic_open(Path(path), "wb") as fh:

src/simurg/pipeline/normalize.py

@@ -0,0 +1,167 @@
+"""``simurg normalize``: ingest output → per-sample normalization transforms and normalized X/Y, θ, R (phase P3).
+
+Reads ``loci.parquet`` and every ``*.intensity.parquet`` under ``--in`` and writes, next to each input::
+
+    <stem>.transforms.parquet   one row per normalization bin (GTC XForm fields)
+    <stem>.normalized.parquet   norm_x, norm_y, theta, r per locus
+
+plus ``normalize_summary.tsv`` and ``normalize_run_manifest.json`` in ``--in``.
+"""
+
+from __future__ import annotations
+
+import multiprocessing as mp
+import os
+from concurrent.futures import ProcessPoolExecutor
+from dataclasses import dataclass, field
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+from simurg import __version__
+from simurg.core.atomic import atomic_write_text, check_outputs_absent
+from simurg.core.errors import ExitCode, SimurgError
+from simurg.core.log import Reporter
+from simurg.core.provenance import RunManifest
+from simurg.core.schemas import (
+    INTENSITY_SCHEMA_ID,
+    LOCI_SCHEMA_ID,
+    normalized_table,
+    read_metadata,
+    read_table,
+    transforms_table,
+    write_parquet,
+)
+from simurg.formats.discovery import walk_files
+from simurg.signal.normalization.apply import normalize, polar
+from simurg.signal.normalization.fit import VERSIONS, fit_transforms
+
+SUFFIX = ".intensity.parquet"
+
+
+@dataclass
+class NormalizeOptions:
+    indir: Path
+    version: str = "1.2.0"
+    threads: int = 0
+    force: bool = False
+    argv: list[str] = field(default_factory=list)
+
+
+_STATE: dict[str, Any] = {}
+
+
+def _init(bins: np.ndarray, lookup: np.ndarray, version: str) -> None:
+    _STATE.update(bins=bins, lookup=lookup, version=version)
+
+
+def _outputs(p: Path) -> tuple[Path, Path]:
+    stem = p.name[: -len(SUFFIX)]
+    return p.with_name(f"{stem}.transforms.parquet"), p.with_name(f"{stem}.normalized.parquet")
+
+
+def _process(path: str) -> dict[str, Any]:
+    p = Path(path)
+    row: dict[str, Any] = {"file": path, "status": "OK", "fail_reason": ""}
+    try:
+        t = read_table(p, INTENSITY_SCHEMA_ID)
+        md = read_metadata(p)
+        row["sample_id"] = md.get("sample_id", "")
+        raw_x = t.column("raw_x").to_numpy().astype(np.uint16)
+        raw_y = t.column("raw_y").to_numpy().astype(np.uint16)
+        bins, lookup, version = _STATE["bins"], _STATE["lookup"], _STATE["version"]
+        if len(raw_x) != len(bins):
+            raise SimurgError("E_BPM_FORMAT", f"{p}: {len(raw_x)} loci but the manifest has {len(bins)}")
+        xforms = fit_transforms(raw_x, raw_y, bins, version)
+        nx, ny = normalize(raw_x, raw_y, lookup, xforms)
+        theta, r = polar(nx, ny)
+        uniq, counts = np.unique(bins, return_counts=True)
+        rows = []
+        for k, (bid, cnt, xf) in enumerate(zip(uniq, counts, xforms, strict=True)):
+            rows.append(
+                {
+                    "bin_index": k,
+                    "norm_bin_id": int(bid),
+                    "n_loci": int(cnt),
+                    "version": xf.version,
+                    **{
+                        f: float(getattr(xf, f))
+                        for f in (
+                            "offset_x",
+                            "offset_y",
+                            "scale_x",
+                            "scale_y",
+                            "shear",
+                            "theta",
+                            "cvx",
+                            "cvy",
+                            "nn12",
+                            "rr12",
+                            "taa",
+                            "tbb",
+                        )
+                    },
+                }
+            )
+        meta = {
+            "sample_id": md.get("sample_id", ""),
+            "sample_name": md.get("sample_name", ""),
+            "normalization_version": version,
+            "simurg_version": __version__,
+            "source_sha256_grn": md.get("grn_sha256", ""),
+            "source_sha256_red": md.get("red_sha256", ""),
+        }
+        tp, np_ = _outputs(p)
+        write_parquet(transforms_table(rows, meta), tp)
+        write_parquet(normalized_table(nx, ny, theta, r, meta), np_)
+        row["n_bins"] = len(rows)
+    except SimurgError as e:
+        row.update(status="FAILED", fail_reason=f"{e.code}: {e.message}", error_code=e.code)
+    except Exception as e:
+        row.update(status="FAILED", fail_reason=f"E_INTERNAL: {type(e).__name__}: {e}", error_code="E_INTERNAL")
+    return row
+
+
+def run_normalize(opts: NormalizeOptions, rep: Reporter) -> int:
+    if opts.version not in VERSIONS:
+        raise SimurgError("E_CONFIG_KEY", f"unknown normalization version {opts.version!r}; use {VERSIONS}")
+    loci = opts.indir / "loci.parquet"
+    if not loci.exists():
+        raise SimurgError("E_FILE_NOT_FOUND", f"{loci} not found — run `simurg ingest` first", path=str(loci))
+    rm = RunManifest(command="normalize", argv=opts.argv)
+    rep.run_id = rm.run_id
+    rm.set_option("normalization_version", opts.version, "cli")
+    lt = read_table(loci, LOCI_SCHEMA_ID)
+    bins = lt.column("norm_bin_id").to_numpy().astype(np.uint8)
+    lookup = lt.column("norm_lookup").to_numpy().astype(np.int32)
+    files = [str(p) for p in walk_files(opts.indir, (SUFFIX,))]
+    if not files:
+        raise SimurgError("E_NO_INPUT", f"no *{SUFFIX} files under {opts.indir}")
+    check_outputs_absent([q for f in files for q in _outputs(Path(f))], force=opts.force)
+
+    workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(files)))
+    rep.info(f"normalizing {len(files)} sample(s) with version {opts.version}, {workers} worker(s)")
+    if workers == 1:
+        _init(bins, lookup, opts.version)
+        results = [_process(f) for f in files]
+    else:
+        with ProcessPoolExecutor(
+            workers, mp_context=mp.get_context("spawn"), initializer=_init, initargs=(bins, lookup, opts.version)
+        ) as ex:
+            results = list(ex.map(_process, files))
+    for r in results:
+        if r["status"] != "OK":
+            rep.error(r.get("error_code", "E_INTERNAL"), f"{r['file']}: {r['fail_reason']}")
+            rm.failed_samples.append({"sample_id": r.get("sample_id", ""), "code": r.get("error_code", "")})
+    cols = ["sample_id", "status", "fail_reason", "n_bins", "file"]
+    lines = ["\t".join(cols)] + ["\t".join(str(r.get(c, "")) for c in cols) for r in results]
+    atomic_write_text(opts.indir / "normalize_summary.tsv", "\n".join(lines) + "\n")
+    ok = sum(r["status"] == "OK" for r in results)
+    rm.summary = {"samples": len(results), "succeeded": ok, "warnings": rep.warning_counts()}
+    rm.status = "completed"
+    rm.write(opts.indir / "normalize_run_manifest.json")
+    rep.info(f"done: {ok}/{len(results)} samples normalized")
+    if ok == len(results):
+        return ExitCode.OK
+    return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED

src/simurg/signal/normalization/__init__.py

@@ -0,0 +1 @@
+"""Per-bin affine intensity normalization (plan §7.2–§7.3)."""

src/simurg/signal/normalization/apply.py

@@ -0,0 +1,69 @@
+"""Apply normalization transforms and convert to polar coordinates — plan §7.3–§7.4.
+
+Exact mode follows the reference implementation used to write GTC files (gtc2vcf ``raw_x_y2norm_x_y`` /
+``norm_x_y2ilmn_theta_r``, MIT): cos θ and sin θ are computed in double and rounded to float once per bin, the
+affine transform runs in float32, negative results are clamped to 0 **before** scaling, and R = |x| + |y|.
+BeadArrayFiles (BSD-2) documents the same transform; it differs only in using double precision internally.
+"""
+
+from __future__ import annotations
+
+import math
+from collections.abc import Sequence
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg.signal.normalization.fit import M_2_PI, XForm
+from simurg.signal.normalization.numerics import TINY
+
+FloatArr = npt.NDArray[np.float32]
+XFormRecords = npt.NDArray[np.void]  # GTC normalization_transforms record array
+
+
+def transforms_to_arrays(xforms: Sequence[XForm] | XFormRecords) -> dict[str, FloatArr]:
+    """Per-bin parameter arrays from XForm objects or a GTC ``normalization_transforms`` record array."""
+    if isinstance(xforms, np.ndarray):
+        get = {
+            k: xforms[k].astype(np.float32) for k in ("offset_x", "offset_y", "scale_x", "scale_y", "shear", "theta")
+        }
+    else:
+        get = {
+            k: np.array([getattr(x, k) for x in xforms], dtype=np.float32)
+            for k in ("offset_x", "offset_y", "scale_x", "scale_y", "shear", "theta")
+        }
+    get["cos"] = np.array([math.cos(float(t)) for t in get["theta"]], dtype=np.float64).astype(np.float32)
+    get["sin"] = np.array([math.sin(float(t)) for t in get["theta"]], dtype=np.float64).astype(np.float32)
+    return get
+
+
+def normalize(
+    raw_x: npt.NDArray[np.uint16],
+    raw_y: npt.NDArray[np.uint16],
+    lookup: npt.NDArray[np.int32],
+    xforms: Sequence[XForm] | XFormRecords,
+) -> tuple[FloatArr, FloatArr]:
+    """Normalized X/Y per locus; ``lookup[i]`` is locus i's transform index (Manifest.normalization_lookups)."""
+    p = transforms_to_arrays(xforms)
+    ox, oy = p["offset_x"][lookup], p["offset_y"][lookup]
+    c, s = p["cos"][lookup], p["sin"][lookup]
+    sh, sx, sy = p["shear"][lookup], p["scale_x"][lookup], p["scale_y"][lookup]
+    tx = raw_x.astype(np.float32) - ox
+    ty = raw_y.astype(np.float32) - oy
+    tx2 = c * tx + s * ty
+    ty2 = (-s) * tx + c * ty
+    tx3 = tx2 - sh * ty2
+    with np.errstate(divide="ignore", invalid="ignore"):
+        nx = np.where(tx3 < 0, np.float32(0), tx3 / sx).astype(np.float32)
+        ny = np.where(ty2 < 0, np.float32(0), ty2 / sy).astype(np.float32)
+    return nx, ny
+
+
+def polar(norm_x: FloatArr, norm_y: FloatArr) -> tuple[FloatArr, FloatArr]:
+    """(θ, R): θ = (2/π)·atan2(y, x), R = |x| + |y|; both NaN where x = y = 0 (plan §7.4)."""
+    theta = (np.arctan2(norm_y.astype(np.float64), norm_x.astype(np.float64)) * M_2_PI).astype(np.float32)
+    r = np.where((norm_x < 0) & (norm_y < 0), TINY, np.abs(norm_x) + np.abs(norm_y)).astype(np.float32)
+    zero = (norm_x == 0) & (norm_y == 0)
+    theta[zero] = np.nan
+    r[zero] = np.nan
+    return theta, r

src/simurg/signal/normalization/fit.py

@@ -0,0 +1,322 @@
+"""Per-sample, per-bin normalization fit — plan §7.2 (versions 1.1.2 and 1.2.0).
+
+Ported from gtc2vcf ``idat2gtc.c`` (MIT, © Giulio Genovese), itself an emulation of the vendor's
+``NormalizationInfinium``/``Normalization111(_Dragonfish)``/``NormalizationDragonfish`` classes; concepts from
+Kermani, US 7,035,740 (expired) and Peiffer et al. 2006. Deliberately reproduced vendor quirks are marked
+``QUIRK``. All intermediate values keep the reference's float32/float64 types (see ``numerics``).
+"""
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg.signal.normalization.numerics import (
+    F32,
+    TINY,
+    FloatArr,
+    MadSigma,
+    iqr,
+    linspace_f,
+    madsigma_new,
+    madsigma_old,
+    nanmean_f,
+    percentile,
+    robustfit0,
+    robustfit1,
+    trimmean,
+    trunc_int,
+    unique_indices,
+)
+
+SAMPLING = 400
+ROBUST_THRESHOLD = 192
+NN_SLOW_LIMIT = 2000  # below this many points the exact O(n²) search is used
+M_1_PI = 0.318309886183790671538  # C <math.h> literals (not 1/np.pi, which may differ in the last bit)
+M_2_PI = 0.636619772367581343076
+VERSIONS = ("1.2.0", "1.1.2")
+
+IdxArr = npt.NDArray[np.int64]
+ZERO = F32(0)
+
+
+@dataclass
+class XForm:
+    """One normalization transform (GTC 52-byte record: version + 12 floats)."""
+
+    version: int = 1
+    offset_x: np.float32 = ZERO
+    offset_y: np.float32 = ZERO
+    scale_x: np.float32 = ZERO
+    scale_y: np.float32 = ZERO
+    shear: np.float32 = ZERO
+    theta: np.float32 = ZERO
+    cvx: np.float32 = ZERO
+    cvy: np.float32 = ZERO
+    nn12: np.float32 = ZERO
+    rr12: np.float32 = ZERO
+    taa: np.float32 = ZERO
+    tbb: np.float32 = ZERO
+
+    def as_record(self) -> tuple[int, float, float, float, float, float, float, list[float]]:
+        return (
+            self.version,
+            float(self.offset_x),
+            float(self.offset_y),
+            float(self.scale_x),
+            float(self.scale_y),
+            float(self.shear),
+            float(self.theta),
+            [float(self.cvx), float(self.cvy), float(self.nn12), float(self.rr12), float(self.taa), float(self.tbb)],
+        )
+
+
+# ---------------------------------------------------------------------------------------------------------
+# nearest neighbours of sweep points
+# ---------------------------------------------------------------------------------------------------------
+def closest_points(xref: FloatArr, yref: FloatArr, x: FloatArr, y: FloatArr) -> IdxArr:
+    """For each sweep point (x[i], y[i]) the index of the closest data point (xref, yref); -1 if none."""
+    if len(xref) < NN_SLOW_LIMIT:
+        return _closest_slow(xref, yref, x, y)
+    return _closest_binned(xref, yref, x, y)
+
+
+def _closest_slow(xref: FloatArr, yref: FloatArr, x: FloatArr, y: FloatArr) -> IdxArr:
+    out = np.empty(len(x), dtype=np.int64)
+    for i in range(len(x)):
+        dx = x[i] - xref
+        dy = y[i] - yref
+        d = dx * dx + dy * dy  # float32, as in the reference
+        out[i] = int(np.argmin(d))  # first minimum wins (strict <)
+    return out
+
+
+def _closest_binned(raw_x: FloatArr, raw_y: FloatArr, axis_x: FloatArr, axis_y: FloatArr) -> IdxArr:
+    """QUIRK: the vendor's binned approximate search (``findClosestSitesToPointsAlongAxis``)."""
+    if (axis_x > F32(0.0001)).any():
+        if (axis_y > F32(0.0001)).any():
+            return np.full(len(axis_x), -1, dtype=np.int64)
+        raw_a, raw_b, axis_a = raw_x, raw_y, axis_x
+    else:
+        raw_a, raw_b, axis_a = raw_y, raw_x, axis_y
+    bin_width = F32(axis_a[-1] / F32(12.0))
+    eligible = np.flatnonzero(raw_b.astype(np.float64) <= float(bin_width))
+    qa = (raw_a[eligible] / bin_width).astype(np.float32)
+    b1 = np.clip(trunc_int(qa), 0, 11)
+    b2 = np.clip(trunc_int((qa - F32(0.5)).astype(np.float32)), 0, 10)
+    bins = [eligible[b1 == k] for k in range(12)]
+    shifted = [eligible[b2 == k] for k in range(11)]
+    out = np.full(len(axis_a), -1, dtype=np.int64)
+    for i in range(len(axis_a)):
+        q = F32(axis_a[i] / bin_width)
+        k = int(np.trunc(q))
+        rem = F32(q - F32(k))
+        k = min(max(k, 0), 11)
+        if F32(0.25) <= rem <= F32(0.75):
+            cand = bins[k]
+        elif rem < F32(0.25):
+            cand = bins[k] if k == 0 else shifted[k - 1]
+        elif k == 11:
+            cand = bins[k]
+        else:
+            cand = shifted[k]
+        if len(cand) == 0:
+            continue
+        dx = raw_x[cand] - axis_x[i]
+        dy = raw_y[cand] - axis_y[i]
+        d = dx * dx + dy * dy
+        j = int(np.argmin(d))
+        if d[j] < 1e20:
+            out[i] = int(cand[j])
+    return out
+
+
+# ---------------------------------------------------------------------------------------------------------
+# steps of the six-parameter model
+# ---------------------------------------------------------------------------------------------------------
+def remove_outliers(x: FloatArr, y: FloatArr) -> tuple[FloatArr, FloatArr]:
+    n = len(x)
+    if n < SAMPLING:
+        return x, y
+    ts = (y / ((TINY + x) + y)).astype(np.float32)
+    xs, ys, tss = np.sort(x), np.sort(y), np.sort(ts)
+    m = 5
+    tcut1 = min(tss[m - 1], percentile(tss, 1))
+    tcut2 = max(tss[n - m], percentile(tss, 99))
+    xcut1 = min(xs[m - 1], percentile(xs, 1))
+    xcut2 = max(xs[n - m], percentile(xs, 99))
+    ycut1 = min(ys[m - 1], percentile(ys, 1))
+    ycut2 = max(ys[n - m], percentile(ys, 99))
+    keep = ~((y <= ycut1) | (x <= xcut1) | (y >= ycut2) | (x >= xcut2))
+    with np.errstate(invalid="ignore", divide="ignore"):
+        t = y.astype(np.float64) / (y + x).astype(np.float32).astype(np.float64)
+    keep &= ~((t <= float(tcut1)) | (t >= float(tcut2)))
+    return x[keep], y[keep]
+
+
+def remove_offset(x: FloatArr, y: FloatArr, madsigma: MadSigma) -> tuple[np.float32, np.float32, IdxArr, IdxArr]:
+    mx = F32(np.nanmin(x))
+    my = F32(np.nanmin(y))
+    xt = (x - mx).astype(np.float32)
+    yt = (y - my).astype(np.float32)
+    xsweep = linspace_f(SAMPLING, F32(0), F32(np.nanmax(xt)))
+    ysweep = linspace_f(SAMPLING, F32(0), F32(np.nanmax(yt)))
+    zeros = np.zeros(SAMPLING, dtype=np.float32)
+    iaa = unique_indices(closest_points(xt, yt, xsweep, zeros))
+    ibb = unique_indices(closest_points(xt, yt, zeros, ysweep))
+    baa, maa = robustfit1(xt[iaa], yt[iaa], madsigma)
+    bbb, mbb = robustfit1(yt[ibb], xt[ibb], madsigma)
+    ox = F32((bbb + mbb * baa) / (F32(1.0) - mbb * maa))
+    oy = F32(baa + maa * ox)
+    return F32(ox + mx), F32(oy + my), iaa, ibb
+
+
+def nn12_rr12(x: FloatArr, y: FloatArr) -> tuple[np.float32, np.float32]:
+    n = len(x)
+    xthr = F32(F32(0.2) * percentile(np.sort(x), 95))
+    ythr = F32(F32(0.2) * percentile(np.sort(y), 95))
+    low = (x < xthr) & (y < ythr)
+    count = int(low.sum())
+    nn12 = F32(F32(count) / F32(n))
+    if not count:
+        return nn12, F32(1.0)
+    xy = (x + y).astype(np.float32)
+    return nn12, F32(nanmean_f(xy[low]) / nanmean_f(xy))
+
+
+def handle_rotation(
+    x: FloatArr, y: FloatArr, iaa: IdxArr, ibb: IdxArr, madsigma: MadSigma
+) -> tuple[FloatArr, FloatArr, np.float32, IdxArr, IdxArr]:
+    xsweep = linspace_f(SAMPLING, F32(np.nanmin(x)), F32(np.nanmax(x)))
+    ysweep = linspace_f(SAMPLING, F32(np.nanmin(y)), F32(np.nanmax(y)))
+    zeros = np.zeros(SAMPLING, dtype=np.float32)
+    iaa = unique_indices(np.concatenate([iaa, closest_points(x, y, xsweep, zeros)]))
+    ibb = unique_indices(np.concatenate([ibb, closest_points(x, y, zeros, ysweep)]))
+    m = robustfit0(x[iaa], y[iaa], madsigma)
+    taa = math.atan(float(m))  # libm, as in the reference (NumPy's SIMD kernels may differ in the last bit)
+    ct, st = math.cos(taa), math.sin(taa)
+    xd, yd = x.astype(np.float64), y.astype(np.float64)
+    nx = (ct * xd + st * yd).astype(np.float32)
+    ny = ((0.0 - st) * xd + ct * yd).astype(np.float32)
+    return nx, ny, F32(taa), iaa, ibb
+
+
+def handle_shear(x: FloatArr, y: FloatArr, ibb: IdxArr, madsigma: MadSigma) -> tuple[FloatArr, np.float32]:
+    ysweep = linspace_f(SAMPLING, F32(0), F32(np.nanmax(y)))
+    zeros = np.zeros(SAMPLING, dtype=np.float32)
+    ibb = unique_indices(np.concatenate([ibb, closest_points(x, y, zeros, ysweep)]))
+    m = robustfit0(y[ibb], x[ibb], madsigma)
+    shy = math.tan(math.atan(float(m)))
+    nx = (x.astype(np.float64) - shy * y.astype(np.float64)).astype(np.float32)
+    return nx, F32(shy)
+
+
+def _atan2_f32(y: FloatArr, x: FloatArr, factor: float) -> FloatArr:
+    """``(float)(factor * atan2((double)y, (double)x))`` element-wise with libm (only used on small bins)."""
+    return np.array(
+        [factor * math.atan2(float(b), float(a)) for a, b in zip(x, y, strict=True)], dtype=np.float64
+    ).astype(np.float32)
+
+
+def base_handle_scale(x: FloatArr, y: FloatArr, version: str) -> tuple[np.float32, np.float32]:
+    """Scale for bins with < 192 loci (and the 1.1.x base path)."""
+    x = np.maximum(x, F32(0)).astype(np.float32)
+    y = np.maximum(y, F32(0)).astype(np.float32)
+    n = len(x)
+    if n < ROBUST_THRESHOLD:
+        ang = _atan2_f32(y, x, 180.0 * M_1_PI)
+        valid = (x > 0) | (y > 0) | (version == "1.1.2")  # QUIRK: 1.1.2 counts failed loci as AA
+        t = np.where(valid, ang, np.float32(np.nan))
+        with np.errstate(invalid="ignore"):
+            xaa, ybb = x[t < F32(10.0)], y[t > F32(80.0)]
+        return trimmean(xaa, 20), trimmean(ybb, 20)
+    xsweep = linspace_f(SAMPLING, F32(0), F32(np.nanmax(x)))
+    ysweep = linspace_f(SAMPLING, F32(0), F32(np.nanmax(y)))
+    zeros = np.zeros(SAMPLING, dtype=np.float32)
+    iaa = unique_indices(closest_points(x, y, xsweep, zeros))
+    ibb = unique_indices(closest_points(x, y, zeros, ysweep))
+    xaa, ybb = x[iaa], y[ibb]
+    sx = F32(F32(0.5) * trimmean(xaa, 50) + F32(0.7413) * iqr(xaa))
+    sy = F32(F32(0.5) * trimmean(ybb, 50) + F32(0.7413) * iqr(ybb))
+    return sx, sy
+
+
+def handle_scale(x: FloatArr, y: FloatArr, version: str) -> tuple[np.float32, np.float32]:
+    if len(x) < ROBUST_THRESHOLD:
+        return base_handle_scale(x, y, version)
+    # QUIRK: the vendor calls its percentile routine on the *unsorted* intensities here
+    xthr = F32(F32(0.1) * percentile(x, 99))
+    ythr = F32(F32(0.1) * percentile(y, 99))
+    aa = (x > F32(5.0) * y) & (x > xthr)
+    bb = (y > F32(5.0) * x) & (y > ythr)
+    return trimmean(x[aa], 50), trimmean(y[bb], 50)
+
+
+def normalize_single_bin(x: FloatArr, y: FloatArr, version: str) -> XForm:
+    madsigma = madsigma_new if version == "1.2.0" else madsigma_old
+    xf = XForm()
+    x, y = remove_outliers(x, y)
+    if len(x) >= ROBUST_THRESHOLD:
+        xf.offset_x, xf.offset_y, iaa, ibb = remove_offset(x, y, madsigma)
+        x = (x - xf.offset_x).astype(np.float32)
+        y = (y - xf.offset_y).astype(np.float32)
+        xf.nn12, xf.rr12 = nn12_rr12(x, y)
+        x, y, xf.theta, iaa, ibb = handle_rotation(x, y, iaa, ibb, madsigma)
+        x, xf.shear = handle_shear(x, y, ibb, madsigma)
+    else:
+        xf.nn12, xf.rr12 = nn12_rr12(x, y)
+    xf.scale_x, xf.scale_y = handle_scale(x, y, version)
+    xf.taa = F32(float(F32(xf.theta * F32(180.0))) * M_1_PI)
+    xf.tbb = F32((math.atan(float(xf.shear)) - float(xf.theta)) * 180.0 * M_1_PI)
+    return xf
+
+
+def normalize_single_bin_single_channel(x: FloatArr, y: FloatArr) -> XForm:
+    """Two-parameter model for small Infinium I bins (normalization 1.2.0 / GenTrain 3.0 tech note)."""
+    swap = y > x
+    x, y = np.where(swap, y, x).astype(np.float32), np.where(swap, x, y).astype(np.float32)
+    ythrsh = percentile(np.sort(y), 2)
+    x = (x - ythrsh).astype(np.float32)
+    y = (y - ythrsh).astype(np.float32)
+    both0 = (x == 0) & (y == 0)
+    r = np.where((x < 0) & (y < 0), TINY, np.abs(x) + np.abs(y)).astype(np.float32)
+    t = _atan2_f32(y, x, M_2_PI)
+    r[both0] = np.nan
+    t[both0] = np.nan
+    with np.errstate(invalid="ignore"):
+        aa = (t < F32(0.1)) & ~np.isnan(r) & (r != TINY)
+    ymean = F32(trimmean(y[aa], 20) - F32(50.0))
+    # QUIRK: NaNs are compacted only within the first j entries and the tail is zero-padded
+    j = int((~np.isnan(r)).sum())
+    head = r[:j]
+    kept = head[~np.isnan(head)]
+    padded = np.concatenate([kept, np.zeros(j - len(kept), dtype=np.float32)])
+    rmean = F32(trimmean(padded, 20) - F32(2.0) * ymean)
+    off = F32(ythrsh + ymean)
+    return XForm(version=1, offset_x=off, offset_y=off, scale_x=rmean, scale_y=rmean, rr12=F32(1.0))
+
+
+def fit_transforms(
+    raw_x: npt.NDArray[np.uint16],
+    raw_y: npt.NDArray[np.uint16],
+    norm_bin_ids: npt.NDArray[np.uint8],
+    version: str = "1.2.0",
+) -> list[XForm]:
+    """Fit one XForm per normalization bin, in ascending bin-id order (= GTC order; plan §6.2)."""
+    if version not in VERSIONS:
+        raise ValueError(f"unknown normalization version {version!r}")
+    out: list[XForm] = []
+    for bid in np.unique(norm_bin_ids):
+        idx = np.flatnonzero(norm_bin_ids == bid)  # loci keep manifest order within a bin
+        if len(idx) < 10:
+            raise ValueError(f"normalization bin {int(bid)} has only {len(idx)} loci (need >= 10)")
+        x = raw_x[idx].astype(np.float32)
+        y = raw_y[idx].astype(np.float32)
+        if version == "1.2.0" and int(bid) >= 100 and len(idx) < ROBUST_THRESHOLD:
+            out.append(normalize_single_bin_single_channel(x, y))
+        else:
+            out.append(normalize_single_bin(x, y, version))
+    return out

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)