From 95e4e41560b93ba617b1ba583e3243a7c0b05c9d Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Sun, 27 Sep 2026 22:51:06 +0100 Subject: [PATCH] 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). --- scripts/make_oracle_gtcs.sh | 26 +- src/simurg/cli/main.py | 24 +- src/simurg/core/schemas.py | 49 +++ src/simurg/pipeline/normalize.py | 167 ++++++++++ src/simurg/signal/normalization/__init__.py | 1 + src/simurg/signal/normalization/apply.py | 69 +++++ src/simurg/signal/normalization/fit.py | 322 ++++++++++++++++++++ src/simurg/signal/normalization/numerics.py | 205 +++++++++++++ tests/conftest.py | 3 +- tests/oracle/test_p2_raw_xy.py | 2 +- tests/oracle/test_p3_normalization.py | 45 +++ tests/unit/test_normalization.py | 198 ++++++++++++ 12 files changed, 1097 insertions(+), 14 deletions(-) create mode 100644 src/simurg/pipeline/normalize.py create mode 100644 src/simurg/signal/normalization/__init__.py create mode 100644 src/simurg/signal/normalization/apply.py create mode 100644 src/simurg/signal/normalization/fit.py create mode 100644 src/simurg/signal/normalization/numerics.py create mode 100644 tests/oracle/test_p3_normalization.py create mode 100644 tests/unit/test_normalization.py diff --git a/scripts/make_oracle_gtcs.sh b/scripts/make_oracle_gtcs.sh index 5417165..ff814c9 100755 --- a/scripts/make_oracle_gtcs.sh +++ b/scripts/make_oracle_gtcs.sh @@ -1,5 +1,8 @@ #!/usr/bin/env bash -# Build bcftools + gtc2vcf plugins (MIT) under refs/tools and generate oracle GTC files from the P0 GSA demo IDATs. +# Build bcftools + gtc2vcf plugins (MIT) under refs/tools and generate oracle GTC files from the P0 GSA demo IDATs, +# once per normalization version (GenTrain 3 = normalization 1.2.0, GenTrain 2 = 1.1.2). +# The plugin is compiled with -ffp-contract=off: compilers may otherwise fuse a*b+c into FMA instructions (default +# on aarch64), which changes the last bits of the normalization transforms. Vendor .NET code and NumPy do not fuse. # Usage: bash scripts/make_oracle_gtcs.sh [chip-folder ...] (default: every chip in the demo set) set -euo pipefail cd "$(dirname "$0")/.." @@ -15,20 +18,21 @@ if [[ ! -x $B/bcftools || ! -f $B/plugins/idat2gtc.so ]]; then rm -f $B/plugins/{idat2gtc,gtc2vcf,affy2vcf}.{c,h} cp refs/gtc2vcf/{idat2gtc,gtc2vcf,affy2vcf}.c refs/gtc2vcf/*.h $B/plugins/ (cd $B/htslib-$V && ./configure --disable-libcurl --disable-bz2 --disable-lzma --disable-ref-cache >/dev/null) - (cd $B && ./configure --without-gsl >/dev/null && make -j"$(nproc)" >/dev/null) + (cd $B && ./configure --without-gsl >/dev/null && make -j"$(nproc)" CFLAGS="-g -O2 -ffp-contract=off" >/dev/null) fi export BCFTOOLS_PLUGINS=$B/plugins D=data/p0/gsa_v4_pgx IDATS=$D/idats_Demo_36/idats_Demo_36 -OUT=data/p0/oracle/gsa_gtc_gtc2vcf -mkdir -p "$OUT" chips=("$@") [[ ${#chips[@]} -eq 0 ]] && mapfile -t chips < <(find "$IDATS" -mindepth 1 -maxdepth 1 -type d | sort) -for chip in "${chips[@]}"; do - $B/bcftools +idat2gtc \ - --bpm $D/GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm \ - --egt $D/GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt \ - --idats "$chip" --output "$OUT" 2>&1 | grep -E "^(Call rate|Gender|Writing)" || true +for v in 3 2; do + OUT=data/p0/oracle/gsa_gt${v}_nofma + mkdir -p "$OUT" + for chip in "${chips[@]}"; do + $B/bcftools +idat2gtc --gentrain-version $v \ + --bpm $D/GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm \ + --egt $D/GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt \ + --idats "$chip" --output "$OUT" 2>&1 | grep -E "^(Call rate|Gender|Writing)" || true + done + echo "oracle GTCs in $OUT: $(ls "$OUT"/*.gtc | wc -l)" done -$B/bcftools +idat2gtc --version 2>&1 | head -2 || true -echo "oracle GTCs in $OUT: $(ls "$OUT"/*.gtc | wc -l)" diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 23fb9be..ef14a39 100644 --- a/src/simurg/cli/main.py +++ b/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() diff --git a/src/simurg/core/schemas.py b/src/simurg/core/schemas.py index 18caadb..72ff2af 100644 --- a/src/simurg/core/schemas.py +++ b/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: diff --git a/src/simurg/pipeline/normalize.py b/src/simurg/pipeline/normalize.py new file mode 100644 index 0000000..80cdcc5 --- /dev/null +++ b/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:: + + .transforms.parquet one row per normalization bin (GTC XForm fields) + .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 diff --git a/src/simurg/signal/normalization/__init__.py b/src/simurg/signal/normalization/__init__.py new file mode 100644 index 0000000..27a4176 --- /dev/null +++ b/src/simurg/signal/normalization/__init__.py @@ -0,0 +1 @@ +"""Per-bin affine intensity normalization (plan §7.2–§7.3).""" diff --git a/src/simurg/signal/normalization/apply.py b/src/simurg/signal/normalization/apply.py new file mode 100644 index 0000000..0725c02 --- /dev/null +++ b/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 diff --git a/src/simurg/signal/normalization/fit.py b/src/simurg/signal/normalization/fit.py new file mode 100644 index 0000000..3e1174f --- /dev/null +++ b/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 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 diff --git a/src/simurg/signal/normalization/numerics.py b/src/simurg/signal/normalization/numerics.py new file mode 100644 index 0000000..64e005b --- /dev/null +++ b/src/simurg/signal/normalization/numerics.py @@ -0,0 +1,205 @@ +"""Numerics with the exact rounding behaviour of the reference normalization (plan Appendix A). + +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) diff --git a/tests/conftest.py b/tests/conftest.py index 935bb64..63fbe0d 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -24,7 +24,8 @@ GDA_EGT = DATA / "gda_v1" / "infinium-global-diversity-array-8-v1-0-D1-cluster-f GDA_BPM = DATA / "gda_v1" / "infinium-global-diversity-array-8-v1-0-D1-manifest-file-bpm" / "GDA-8v1-0_D1.bpm" GDA_SHEET = DATA / "gda_v1" / "infinium-global-diversity-array-8-v1-0_D1_SampleSheet_Demo_36.csv" BIOC_IDATS = DATA / "bioc" / "IlluminaDataTestFiles_1.50.0" / "IlluminaDataTestFiles" / "inst" / "extdata" / "idat" -ORACLE_GTCS = DATA / "oracle" / "gsa_gtc_gtc2vcf" +ORACLE_GTCS = DATA / "oracle" / "gsa_gt3_nofma" # gtc2vcf, normalization 1.2.0, no FMA +ORACLE_GTCS_GT2 = DATA / "oracle" / "gsa_gt2_nofma" # normalization 1.1.2 def need(*paths: Path) -> None: diff --git a/tests/oracle/test_p2_raw_xy.py b/tests/oracle/test_p2_raw_xy.py index 2642bb4..40919c4 100644 --- a/tests/oracle/test_p2_raw_xy.py +++ b/tests/oracle/test_p2_raw_xy.py @@ -6,12 +6,12 @@ from __future__ import annotations import numpy as np import pytest +from tests.conftest import GSA_BPM, GSA_IDATS, ORACLE_GTCS, need from simurg.formats.bpm import read_bpm from simurg.formats.gtc import read_gtc from simurg.formats.idat import read_idat from simurg.signal.raw_xy import LocusMapCache, assemble_raw_xy, check_coverage -from tests.conftest import GSA_BPM, GSA_IDATS, ORACLE_GTCS, need pytestmark = pytest.mark.data diff --git a/tests/oracle/test_p3_normalization.py b/tests/oracle/test_p3_normalization.py new file mode 100644 index 0000000..4326d8c --- /dev/null +++ b/tests/oracle/test_p3_normalization.py @@ -0,0 +1,45 @@ +"""Phase P3 acceptance on P0 data: normalization transforms fitted from raw intensities are bit-identical to the +transforms stored in oracle GTCs (gtc2vcf built without FMA), for normalization 1.2.0 and 1.1.2. +The same check against vendor GTCs runs in scripts/validate_p1_large.py --p3. +""" + +from __future__ import annotations + +import numpy as np +import pytest +from tests.conftest import GSA_BPM, ORACLE_GTCS, ORACLE_GTCS_GT2, need + +from simurg.formats.bpm import read_bpm +from simurg.formats.gtc import read_gtc +from simurg.signal.normalization.fit import fit_transforms + +pytestmark = pytest.mark.data + +FIELDS = ("offset_x", "offset_y", "scale_x", "scale_y", "shear", "theta") + + +def as_bits(xf) -> np.ndarray: # type: ignore[no-untyped-def,type-arg] + v = [getattr(xf, f) for f in FIELDS] + [xf.nn12, xf.rr12, xf.taa, xf.tbb] + return np.array(v, dtype=np.float32).view(np.int32) + + +def ref_bits(rec) -> np.ndarray: # type: ignore[no-untyped-def,type-arg] + v = [rec[f] for f in FIELDS] + list(rec["reserved"][2:6]) + return np.array(v, dtype=np.float32).view(np.int32) + + +@pytest.mark.parametrize( + ("folder", "version"), [(ORACLE_GTCS, "1.2.0"), (ORACLE_GTCS_GT2, "1.1.2")], ids=["1.2.0", "1.1.2"] +) +def test_transforms_bit_identical(folder, version) -> None: # type: ignore[no-untyped-def] + need(GSA_BPM, folder) + m = read_bpm(GSA_BPM) + bins = m.norm_bin_ids() + gtcs = sorted(folder.glob("*.gtc")) + assert len(gtcs) >= 18 + for g in gtcs: + ref = read_gtc(g) + mine = fit_transforms(ref.raw_x, ref.raw_y, bins, version) + assert len(mine) == len(ref.normalization_transforms) + for k, (a, b) in enumerate(zip(mine, ref.normalization_transforms, strict=True)): + assert np.array_equal(as_bits(a), ref_bits(b)), f"{g.name} bin {k}" diff --git a/tests/unit/test_normalization.py b/tests/unit/test_normalization.py new file mode 100644 index 0000000..8286ada --- /dev/null +++ b/tests/unit/test_normalization.py @@ -0,0 +1,198 @@ +"""Phase P3: normalization numerics, fit and apply (plan §7.2–§7.4, Appendix A).""" + +from __future__ import annotations + +import numpy as np +import pytest +from typer.testing import CliRunner + +from simurg.signal.normalization import numerics as nm +from simurg.signal.normalization.apply import normalize, polar +from simurg.signal.normalization.fit import ( + ROBUST_THRESHOLD, + _closest_binned, + _closest_slow, + fit_transforms, + normalize_single_bin_single_channel, +) + +F32 = np.float32 + + +def test_percentile_vendor_rule() -> None: + v = np.arange(1, 11, dtype=np.float32) # 1..10 + # n=10, p=50: i1=5, f=0 -> i1=4; x1=45, x2=55; y = 5 + (6-5)/10*(50-45) = 5.5 + assert nm.percentile(v, 50) == F32(5.5) + assert nm.percentile(v, 0) == F32(1) and nm.percentile(v, 100) == F32(10) + assert np.isnan(nm.percentile(np.zeros(0, np.float32), 50)) + + +def test_trimmean_and_iqr() -> None: + # n=20: vendor bounds P10 = 2.5 and P90 = 18.5 -> keeps 3..18 (the 1000 outlier is trimmed) + v = np.array([*range(1, 20), 1000], dtype=np.float32) + assert nm.trimmean(v, 20) == pytest.approx(10.5) + # with only 5 values the vendor rule puts P90 at the maximum, so nothing is trimmed from the top + assert nm.trimmean(np.array([1, 2, 3, 4, 100], np.float32), 20) == pytest.approx(22.0) + assert nm.iqr(np.arange(0, 101, dtype=np.float32)) == pytest.approx(50.0, abs=1.0) + + +def test_seqsum_is_sequential() -> None: + a = np.array([0.1, 0.2, 0.3, 1e16, -1e16, 0.4]) + assert nm.seqsum(a) == 0.4 # left-to-right; a compensated sum would give 1.0 + + +@pytest.mark.parametrize("madsigma", [nm.madsigma_new, nm.madsigma_old]) +def test_robustfit_recovers_lines_despite_outliers(madsigma) -> None: # type: ignore[no-untyped-def] + rng = np.random.default_rng(3) + x = rng.uniform(0, 100, 300).astype(np.float32) + y = (2.5 * x + 7 + rng.normal(0, 0.5, 300)).astype(np.float32) + y[:20] += 500 # gross outliers + b, m = nm.robustfit1(x, y, madsigma) + assert float(m) == pytest.approx(2.5, abs=0.02) and float(b) == pytest.approx(7, abs=1.0) + m0 = nm.robustfit0(x, (3.0 * x).astype(np.float32), madsigma) + assert float(m0) == pytest.approx(3.0, rel=1e-6) + + +def test_linspace_and_unique() -> None: + s = nm.linspace_f(5, F32(0), F32(1)) + assert s.tolist() == [0.0, 0.25, 0.5, 0.75, 1.0] and s.dtype == np.float32 + assert nm.unique_indices(np.array([3, -1, 3, 0, -1])).tolist() == [0, 3] + + +def test_nearest_neighbour_searches() -> None: + rng = np.random.default_rng(0) + xr = rng.uniform(0, 1000, 5000).astype(np.float32) + yr = rng.uniform(0, 1000, 5000).astype(np.float32) + sweep = nm.linspace_f(400, F32(0), F32(xr.max())) + zeros = np.zeros(400, np.float32) + exact = _closest_slow(xr, yr, sweep, zeros) + binned = _closest_binned(xr, yr, sweep, zeros) + assert (binned >= -1).all() and (binned < 5000).all() + # the binned search is approximate but finds the exact neighbour for most sweep points + found = binned >= 0 + assert (binned[found] == exact[found]).mean() > 0.9 + + +def synth_bin( + n: int, rng: np.random.Generator, *, offset=(300.0, 500.0), scale=(8000.0, 12000.0), theta=0.05, shear=0.03 +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: # type: ignore[type-arg] + """Canonical AA/AB/BB clusters distorted by the inverse of the vendor transform. + + Like real data, each genotype cluster is a streak along its ray (intensity varies 0.3–1.6×), which is what + the vendor's axis-sweep line fits rely on. + """ + g = rng.choice(3, n, p=[0.4, 0.3, 0.3]) + inten = rng.uniform(0.3, 1.6, n) + cx = inten * np.select([g == 0, g == 1, g == 2], [1.0, 0.5, 0.02]) + rng.normal(0, 0.01, n) + cy = inten * np.select([g == 0, g == 1, g == 2], [0.02, 0.5, 1.0]) + rng.normal(0, 0.01, n) + tx3, ty2 = cx * scale[0], cy * scale[1] + tx2 = tx3 + shear * ty2 + c, s = np.cos(theta), np.sin(theta) + tx = c * tx2 - s * ty2 + ty = s * tx2 + c * ty2 + raw_x = np.clip(np.round(tx + offset[0]), 0, 65535).astype(np.uint16) + raw_y = np.clip(np.round(ty + offset[1]), 0, 65535).astype(np.uint16) + return raw_x, raw_y, g + + +@pytest.mark.parametrize("version", ["1.2.0", "1.1.2"]) +def test_fit_undoes_known_distortion(version) -> None: # type: ignore[no-untyped-def] + rng = np.random.default_rng(11) + rx, ry, g = synth_bin(3000, rng) + bins = np.zeros(3000, np.uint8) + xf = fit_transforms(rx, ry, bins, version) + assert len(xf) == 1 + # the AA streak itself is tilted by atan(0.02·12000/8000) ≈ 0.03 rad, which the fit also removes + assert 0.05 < float(xf[0].theta) < 0.1 + nx, ny = normalize(rx, ry, np.zeros(3000, np.int32), xf) + th, r = polar(nx, ny) + assert np.nanmedian(th[g == 0]) < 0.08 and np.nanmedian(th[g == 2]) > 0.92 + assert 0.35 < np.nanmedian(th[g == 1]) < 0.65 + assert 0.8 < np.nanmedian(r[g != 1]) < 1.3 + + +def test_bins_ordered_and_small_bins() -> None: + rng = np.random.default_rng(5) + rx1, ry1, _ = synth_bin(1000, rng) + rx2, ry2, _ = synth_bin(150, rng) # < 192 loci: offset/rotation/shear skipped + rx = np.concatenate([rx2, rx1]) + ry = np.concatenate([ry2, ry1]) + bins = np.concatenate([np.full(150, 7, np.uint8), np.full(1000, 3, np.uint8)]) + xf = fit_transforms(rx, ry, bins, "1.2.0") + assert len(xf) == 2 # ascending bin id: 3 then 7 + assert xf[0].theta != 0 and xf[1].theta == 0 and xf[1].offset_x == 0 and xf[1].shear == 0 + with pytest.raises(ValueError): + fit_transforms(rx[:5], ry[:5], np.zeros(5, np.uint8)) + with pytest.raises(ValueError): + fit_transforms(rx, ry, bins, "9.9") + + +def test_single_channel_model_for_small_infinium_i_bins() -> None: + rng = np.random.default_rng(8) + n = 120 + a = rng.uniform(2000, 9000, n) + b = rng.uniform(200, 600, n) + x = np.where(rng.random(n) < 0.5, a, b).astype(np.float32) + y = np.where(x > 1000, b, a).astype(np.float32) + xf = normalize_single_bin_single_channel(x, y) + assert xf.offset_x == xf.offset_y and xf.scale_x == xf.scale_y and xf.theta == 0 and xf.shear == 0 + assert xf.rr12 == 1 and np.isfinite(xf.scale_x) + bins = np.full(n, 105, np.uint8) # Infinium I bin (id >= 100), < 192 loci + rx, ry = x.astype(np.uint16), y.astype(np.uint16) + via_fit = fit_transforms(rx, ry, bins, "1.2.0")[0] + assert via_fit.scale_x == via_fit.scale_y # 1.2.0 → 2-parameter model + assert fit_transforms(rx, ry, bins, "1.1.2")[0].scale_x != via_fit.scale_x or n >= ROBUST_THRESHOLD + + +def test_polar_edge_cases() -> None: + nx = np.array([0, 1, 0, 1], np.float32) + ny = np.array([0, 0, 1, 1], np.float32) + th, r = polar(nx, ny) + assert np.isnan(th[0]) and np.isnan(r[0]) + assert th[1:].tolist() == [0.0, 1.0, 0.5] and r[1:].tolist() == [1.0, 1.0, 2.0] + + +def test_normalize_cli_on_ingest_output(tmp_path) -> None: # type: ignore[no-untyped-def] + import pyarrow.parquet as pq + + from simurg.cli.main import app + from simurg.synth.files import SynthIdat, SynthLocus, encode_bpm, encode_idat + + rng = np.random.default_rng(2) + n = 600 + rx, ry, _ = synth_bin(n, rng) + addr = np.arange(1, n + 1, dtype=np.int32) + for ch, mean in (("Red", rx), ("Grn", ry)): + s = SynthIdat( + address_ids=addr, mean=mean, n_beads=np.full(n, 12, np.uint8), barcode="200000000001", position="R01C01" + ) + (tmp_path / f"200000000001_R01C01_{ch}.idat").write_bytes(encode_idat(s)) + (tmp_path / "m.bpm").write_bytes(encode_bpm([SynthLocus(f"rs{i}", address_a=int(a)) for i, a in enumerate(addr)])) + out = tmp_path / "out" + runner = CliRunner() + assert ( + runner.invoke( + app, + [ + "ingest", + "--quiet", + "--bpm", + str(tmp_path / "m.bpm"), + "--idats", + str(tmp_path), + "--out", + str(out), + "--threads", + "1", + ], + ).exit_code + == 0 + ) + r = runner.invoke(app, ["normalize", "--quiet", "--in", str(out), "--threads", "1"]) + assert r.exit_code == 0, r.output + t = pq.read_table(next(out.rglob("*.normalized.parquet"))) + assert t.num_rows == n and set(t.column_names) == {"norm_x", "norm_y", "theta", "r"} + tr = pq.read_table(next(out.rglob("*.transforms.parquet"))) + assert tr.num_rows == 1 and tr.column("n_loci").to_pylist() == [n] + r = runner.invoke(app, ["normalize", "--quiet", "--in", str(out)]) + assert r.exit_code == 2 and "E_OUTPUT_EXISTS" in r.output