CN: per-probe heterozygous BAF spread from the batch reference; --mosaic-prior/--switch-distance; Pi memory guard

scripts/pi_run.sh runs heavy jobs in their own process group and stops them when
MemAvailable drops below a floor (Raspberry Pi OS boots with cgroup_disable=memory, so a
systemd MemoryMax cap is not available); a 360-sample genotype run with cohort QC had
exhausted memory and the OOM killer took down the desktop session.

src/simurg/cn/hmm.py

@@ -54,7 +54,7 @@ def _d_split(kind: str, f: float) -> float:
     return f / 2
 
 
-def diploid_states() -> list[State]:
+def diploid_states(mosaic_prior: float = MOSAIC_PRIOR) -> list[State]:
     s = [
         State("CN0", 0, False, "", 1.0, "DEL", LRR_MEAN[0], (), 0.01),
         State("CN1", 1, False, "", 1.0, "DEL", LRR_MEAN[1], (("h0", 0.0), ("h1", 1.0)), 0.06),
@@ -71,11 +71,11 @@ def diploid_states() -> list[State]:
             d = _d_split(kind, f)
             cn = {"loss": 1, "gain": 3, "cnloh": 2}[kind]
             bands = (("aa", 0.0), ("ab-", 0.5 - d), ("ab+", 0.5 + d), ("bb", 1.0))
-            s.append(State(f"M{kind}{f:.1f}", cn, kind == "cnloh", kind, f, call, mu_end * f, bands, MOSAIC_PRIOR / 12))
+            s.append(State(f"M{kind}{f:.1f}", cn, kind == "cnloh", kind, f, call, mu_end * f, bands, mosaic_prior / 12))
     return s
 
 
-def haploid_states() -> list[State]:
+def haploid_states(mosaic_prior: float = MOSAIC_PRIOR) -> list[State]:
     s = [
         State("H0", 0, False, "", 1.0, "DEL", LRR_MEAN[0], (), 0.02),
         State("H1", 1, False, "", 1.0, "", 0.0, (("h0", 0.0), ("h1", 1.0)), 0.80),
@@ -83,9 +83,9 @@ def haploid_states() -> list[State]:
     ]
     for f in FRACTIONS:
         s.append(State(f"Hgain{f:.1f}", 2, False, "gain", f, "MDUP", 0.45 * f, (("h0", 0.0), ("h1", 1.0)),
-                       MOSAIC_PRIOR / 8))  # fmt: skip
+                       mosaic_prior / 8))  # fmt: skip
         s.append(State(f"Hloss{f:.1f}", 0, False, "loss", f, "MDEL", -1.0 * f, (("h0", 0.0), ("h1", 1.0)),
-                       MOSAIC_PRIOR / 8))  # fmt: skip
+                       mosaic_prior / 8))  # fmt: skip
     return s
 
 
@@ -100,7 +100,7 @@ def _geno_weights(key: str, p: F64) -> F64:
     return table[key]
 
 
-def _norm(x: F64, mu: float, sd: float) -> F64:
+def _norm(x: F64, mu: float, sd: float | F64) -> F64:
     return np.exp(-0.5 * ((x - mu) / sd) ** 2) / (sd * math.sqrt(2 * math.pi))
 
 
@@ -114,8 +114,9 @@ class Params:
     extra: dict[str, float] = field(default_factory=dict)
 
 
-def emissions(states: list[State], lrr: F64, baf: F64, pfb: F64, par: Params) -> F64:
-    """Per-probe, per-state emission likelihood (n × S); missing LRR/BAF contribute 1."""
+def emissions(states: list[State], lrr: F64, baf: F64, pfb: F64, par: Params, het_sd: F64 | None = None) -> F64:
+    """Per-probe, per-state emission likelihood (n × S); missing LRR/BAF contribute 1. ``het_sd``: per-probe spread
+    of heterozygous BAF in the reference; heterozygous components use max(σ_b, het_sd)."""
     n = len(lrr)
     out = np.ones((n, len(states)), dtype=np.float64)
     lok = np.isfinite(lrr)
@@ -123,13 +124,17 @@ def emissions(states: list[State], lrr: F64, baf: F64, pfb: F64, par: Params) ->
     x = np.where(lok, lrr, 0.0)
     b = np.clip(np.where(bok, baf, 0.5), 0.0, 1.0)
     p = np.clip(np.where(np.isfinite(pfb), pfb, 0.5), 0.01, 0.99)
+    sd_het: float | F64 = (
+        np.maximum(par.baf_sd, np.nan_to_num(het_sd, nan=par.baf_sd)) if het_sd is not None else par.baf_sd
+    )
     for j, st in enumerate(states):
         sd = 2.0 if st.cn == 0 and not st.mosaic else par.lrr_sd
         le = (1 - par.eps) * _norm(x, st.lrr_mu + par.lrr_offset, sd) + par.eps / 8.0
         if st.baf:
             be = np.zeros(n)
             for key, mean in st.baf:
-                be += _geno_weights(key, p) * _norm(b, mean, par.baf_sd)
+                bsd = par.baf_sd if mean in (0.0, 1.0) or het_sd is None else sd_het
+                be += _geno_weights(key, p) * _norm(b, mean, bsd)
             be = (1 - par.eps) * be + par.eps
         else:
             be = np.ones(n)

src/simurg/cn/model.py

@@ -71,3 +71,18 @@ def reference_stats(
         np.nan_to_num(het_c, nan=0.5).astype(np.float32),
         n_het,
     )
+
+
+def reference_het_spread(baf: npt.NDArray[np.float32]) -> npt.NDArray[np.float32]:
+    """Per-probe robust SD (1.4826·MAD) of heterozygous BAF (0.2 < BAF < 0.8) across reference samples; NaN when
+    fewer than 5 heterozygous samples."""
+    het = np.where((baf > 0.2) & (baf < 0.8), baf, np.nan)
+    with np.errstate(all="ignore"):
+        import warnings
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", RuntimeWarning)
+            med = np.nanmedian(het, axis=0)
+            mad = np.nanmedian(np.abs(het - med), axis=0)
+    n = np.isfinite(het).sum(axis=0)
+    return np.where(n >= 5, 1.4826 * mad, np.nan).astype(np.float32)