cn build-model: compute the per-probe LRR spread in column chunks (peak memory hit the 5 GB cap on 120 GSA arrays)

src/simurg/cn/model.py

@@ -94,10 +94,14 @@ def reference_lrr_spread(lrr: npt.NDArray[np.float32], ref_lrr: npt.NDArray[np.f
     fewer than 5 finite values."""
     import warnings
 
-    x = np.where(np.isfinite(lrr), lrr, np.nan).astype(np.float32)
+    out = np.full(lrr.shape[1], np.nan, dtype=np.float32)
     with warnings.catch_warnings():
         warnings.simplefilter("ignore", RuntimeWarning)
-        x = x - np.nanmedian(x, axis=1, keepdims=True) - ref_lrr[None, :]
+        smed = np.nanmedian(np.where(np.isfinite(lrr), lrr, np.nan), axis=1, keepdims=True).astype(np.float32)
-        mad = np.nanmedian(np.abs(x - np.nanmedian(x, axis=0)), axis=0)
+        for j in range(0, lrr.shape[1], 20_000):  # column chunks keep the temporaries small
-    n = np.isfinite(x).sum(axis=0)
+            x = np.where(np.isfinite(lrr[:, j : j + 20_000]), lrr[:, j : j + 20_000], np.nan) - smed
-    return np.where(n >= 5, 1.4826 * mad, np.nan).astype(np.float32)
+            x -= ref_lrr[None, j : j + 20_000]
+            mad = np.nanmedian(np.abs(x - np.nanmedian(x, axis=0)), axis=0)
+            n = np.isfinite(x).sum(axis=0)
+            out[j : j + 20_000] = np.where(n >= 5, 1.4826 * mad, np.nan)
+    return out