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).
This commit is contained in:
kemal committed
1 parent
85229fba84
commit
d163b34eca
12 files changed
+1097
-14
No files matched your search
@@ -0,0 +1 @@
|
||||
"""Per-bin affine intensity normalization (plan §7.2–§7.3)."""
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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)
|
||||
Reference in new issue
Block a user