206 lines
7.0 KiB
Python
206 lines
7.0 KiB
Python
"""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)
|