fix(dynacell/eval): clip the CP shared-space z at ±20 before KID/FID/cosine

applications/dynacell/src/dynacell/evaluation/README.md

@@ -24,7 +24,7 @@ End-to-end evaluation pipeline for virtual staining predictions against fluoresc
 - `io.gt_path` — fluorescence ground truth (channel: `io.gt_channel_name`)
 - `io.cell_segmentation_path` — *optional* precomputed cell segmentation HCS OME-Zarr. Required when `compute_feature_metrics=true` or when building CP/DINOv3/DynaCLR/CELL-DINO cache entries. Position layout must match GT/pred 1:1.
 - `io.gt_cache_dir`, `io.pred_cache_dir` — *optional* artifact cache directories; must be distinct. See [Caches](#caches).
-- `feature_metrics.cp.reference_path` — the target's CP (GLCM+) reference, required when `compute_feature_metrics=true`. `null` (default) resolves to `DATA_ROOT/cp_reference/<target_name>.json` for full and lite evals alike. The reference holds one feature mask chosen on GT cells only, pooled over iPSC + A549 mock/denv/zikv, and one GT mean/std per test set (keyed by `benchmark.dataset_ref.dataset`; HEK gets its own, lite reuses its parent's). Every eval applies the mask and its dataset's scaler to pred AND GT. The eval refuses a partial position walk, and refuses to score unless its staged GT CP cells match the reference's fit for that dataset: the same GT-finite cell count exactly, and every CP feature's raw GT mean/std within a relative tolerance of 1e-6 (`GT_MOMENT_RTOL`). A re-cache that reproduces the cells passes, including a GPU recompute (reproducible only to ~1e-15) and a CPU-vs-GPU recompute (up to ~4e-8, measured). It is a moments gate, not a content hash: a change that moves any feature's GT mean/std by more than the tolerance fails, but permuting values within a column, moving cells between blocks, or a single-cell change below about n × 1e-6 std passes. When the GT CP cache is already complete, the check also runs before the FOV loop. Build the reference from existing GT caches with `tools/build_cp_reference.py --target <target>`; `--verify` re-applies the same check to the current caches. `metrics_provenance.json` stamps the per-dataset `cp_space_sha256`: the mask, feature names, recipe, and that dataset's scaler, recorded GT position set, GT cell count and GT moments. A rebuild therefore invalidates only the eval dirs whose own dataset's CP space changed; adding or refitting another dataset leaves them reusable. The `cp_selected_feature_mask.json` sidecar records the same per-dataset hash.
+- `feature_metrics.cp.reference_path` — the target's CP (GLCM+) reference, required when `compute_feature_metrics=true`. `null` (default) resolves to `DATA_ROOT/cp_reference/<target_name>.json` for full and lite evals alike. The reference holds one feature mask chosen on GT cells only, pooled over iPSC + A549 mock/denv/zikv, and one GT mean/std per test set (keyed by `benchmark.dataset_ref.dataset`; HEK gets its own, lite reuses its parent's). Every eval applies the mask and its dataset's scaler to pred AND GT. The eval refuses a partial position walk, and refuses to score unless its staged GT CP cells match the reference's fit for that dataset: the same GT-finite cell count exactly, and every CP feature's raw GT mean/std within a relative tolerance of 1e-6 (`GT_MOMENT_RTOL`). A re-cache that reproduces the cells passes, including a GPU recompute (reproducible only to ~1e-15) and a CPU-vs-GPU recompute (up to ~4e-8, measured). It is a moments gate, not a content hash: a change that moves any feature's GT mean/std by more than the tolerance fails, but permuting values within a column, moving cells between blocks, or a single-cell change below about n × 1e-6 std passes. When the GT CP cache is already complete, the check also runs before the FOV loop. Build the reference from existing GT caches with `tools/build_cp_reference.py --target <target>`; `--verify` re-applies the same check to the current caches. `metrics_provenance.json` stamps the per-dataset `cp_space_sha256`: the mask, feature names, recipe, and that dataset's scaler, recorded GT position set, GT cell count and GT moments. A rebuild therefore invalidates only the eval dirs whose own dataset's CP space changed; adding or refitting another dataset leaves them reusable. The `cp_selected_feature_mask.json` sidecar records the same per-dataset hash. Before KID/FID/cosine, the shared-space z is clipped to ±`z_clip` (`CP_Z_CLIP` = 20) on both sides, at dataset level and per row, so one heavy-tailed predicted feature cannot set the cubic KID kernel's magnitude. The clip may touch at most `CP_GT_CLIP_FRAC_MAX` = 1e-3 of each FULL test set's GT cells (any |z| > 20 in its own scaler); lite test sets record the same fraction but are not held to the bound (see below). Measured on the registry's full test sets, the GT max |z| is 26.55 (er / a549-mantis-sec61b-denv / kurtosis), 2 of 23,970 distinct GT cells exceed 20, and the highest enforced per-dataset fraction is 5.3e-4 (er a549-mantis-sec61b-denv and -mock, one cell each). The unenforced lite record reaches 9.5e-4 (a549-mantis-sec61b-denv-lite, the same denv cell in a smaller set). The builder re-checks the bound on every build and refuses to build a reference that violates it. Lite datasets are exempt (their fraction is recorded with `enforced: false`): their GT cells are a subset of the parent's, scored in the parent's scaler, so the lite fraction is small-sample noise on the parent's, which is enforced. c = 20 comes from an offline sweep over c ∈ {5, 8, 10, 15, 20, 30, 50, none} on the nucleus-lite pilot: the in-domain model order was the same at every c except one near-tie (zikv at c = 8), and the cross-domain model order was identical for all c ≥ 20. Caveat: the cubic KID kernel keeps growing with c, so cross-domain CP KID magnitudes are not comparable to in-domain ones; only orderings are meaningful. What the clip discards is reported: `Dataset_CP_clip_frac` and the per-row `CP_clip_frac` give the fraction of pred cells with any kept feature beyond the clip, and the sidecar's `clip` block gives the same per feature. The linear probe still uses the masked, unscaled features.
   - GLCM+ (CP) space in two lines: one feature mask selected on GT cells only, pooled over the target's iPSC + A549 mock/denv/zikv test sets; then each test set's own GT mean/std standardizes pred and GT alike, so offsets and scale errors relative to that test set's GT stay visible.
   - GLCM+ median cosine is invariant to contraction toward the GT mean: for `pred = mean + 0.5·(GT − mean)`, every pred/GT cell pair is parallel after the shared scaler and the median cosine is exactly 1.0000, while KID catches it (reviewer-measured on real cells: KID 0.194 vs a pred==GT floor of −0.041; reproduced on synthetic cells in `cp_reference_test.py`). Cosine measures direction only; over-smoothing toward the mean shows up in KID, not cosine.
 

applications/dynacell/src/dynacell/evaluation/cp_reference.py

@@ -75,7 +75,7 @@
 from dynacell.evaluation.pipeline_cache import cp_recipe_identity
 
 #: Schema of the reference JSON. Bump on any change to its keys or their meaning.
-CP_REFERENCE_SCHEMA = 5
+CP_REFERENCE_SCHEMA = 6
 
 #: Sidecar each eval writes beside its metrics, naming the CP reference it scored in.
 CP_SIDECAR_FILENAME = "cp_selected_feature_mask.json"
@@ -90,6 +90,27 @@
 #: must both be <= GT_MOMENT_ZERO_STD_ATOL * (1 + |mean_fit|).
 GT_MOMENT_ZERO_STD_ATOL = 1e-12
 
+#: Shared-space z is clipped to +-CP_Z_CLIP on both sides before KID/FID/cosine
+#: (:meth:`DatasetCPSpace.transform_clipped`). The cubic KID kernel otherwise lets one
+#: heavy-tailed predicted feature set the metric's magnitude (a pilot A549->iPSC nucleus
+#: eval had pred glcm_ASM at z up to 239, KID 4e9). c = 20 was chosen by the paper owner
+#: from an offline sweep over c in {5, 8, 10, 15, 20, 30, 50, none} on the nucleus-lite
+#: pilot (see :data:`CP_Z_CLIP_RANK_STABILITY`). Surveying every registry GT cell of the
+#: full test sets in its own test set's scaler, 2 of 23,970 distinct cells exceed it (max
+#: |z| 26.55: er, a549-mantis-sec61b-denv, kurtosis); enforced max clip fraction 5.3e-4.
+CP_Z_CLIP = 20.0
+
+#: Per-dataset bound on the GT cells the clip may touch: the fraction of a dataset's GT
+#: cells with any kept-feature |z| > CP_Z_CLIP in its own scaler. :func:`fit_cp_reference`
+#: refuses to build a reference that violates it for any NON-lite dataset, so it is a
+#: checked invariant. Lite datasets record their fraction with ``enforced: false``.
+CP_GT_CLIP_FRAC_MAX = 1e-3
+
+#: Why c = 20: recorded in the hashed criteria next to the clip.
+CP_Z_CLIP_RANK_STABILITY = (
+    "cross-domain model order on the nucleus-lite pilot identical for all c >= 20 (swept 5,8,10,15,20,30,50,none)"
+)
+
 #: A per-dataset std below this fraction of the feature's pooled GT std is floored to it.
 STD_FLOOR_FRACTION = 1e-3
 
@@ -471,6 +492,11 @@ def fit_cp_reference(
         pooled, freq_cut=DEFAULT_FREQ_CUT, unique_cut=DEFAULT_UNIQUE_CUT, corr_threshold=DEFAULT_CORR_THRESHOLD
     )
     kept_names = [n for n, k in zip(feature_names, keep_mask, strict=True) if k]
+    if not kept_names:
+        raise ValueError(
+            f"GT-only feature selection kept no CP features out of {len(feature_names)} "
+            f"({pooled.shape[0]} pooled GT cells); a reference with an empty mask cannot score CP."
+        )
     pooled_std = pooled[:, keep_mask].std(axis=0)
     if not (pooled_std > 0).all():
         raise ValueError(f"pooled fit kept zero-variance features: {np.array(kept_names)[pooled_std == 0]}")
@@ -496,6 +522,7 @@ def _provenance(fit: DatasetFit) -> dict[str, Any]:
             "floored_features": [n for n, f in zip(kept_names, floored, strict=True) if f],
         }
     lite = [f for f in fits if f.parent is not None]
+    survey = _gt_abs_z_survey(fits, scalers, keep_mask, kept_names)
 
     payload: dict[str, Any] = {
         "schema": CP_REFERENCE_SCHEMA,
@@ -510,6 +537,10 @@ def _provenance(fit: DatasetFit) -> dict[str, Any]:
             "mask_fit_on": "gt_only_pooled",
             "scaler_fit_on": "gt_only_per_dataset",
             "std_floor_fraction_of_pooled": STD_FLOOR_FRACTION,
+            "z_clip": CP_Z_CLIP,
+            "gt_clip_frac_max": CP_GT_CLIP_FRAC_MAX,
+            "rank_stability": CP_Z_CLIP_RANK_STABILITY,
+            "gt_abs_z_survey": {k: v for k, v in survey.items() if k != "datasets"},
         },
         "cp_identity": cp_identity,
         "scalers": scalers,
@@ -521,15 +552,79 @@ def _provenance(fit: DatasetFit) -> dict[str, Any]:
                 "gt_matrix_sha256": gt_matrix_sha256(pooled),
                 "pooled_std": [float(v) for v in pooled_std],
             },
-            "datasets": {f.dataset: _provenance(f) for f in scaled},
-            "lite": {f.dataset: _provenance(f) for f in lite},
+            "datasets": {f.dataset: {**_provenance(f), "gt_abs_z": survey["datasets"][f.dataset]} for f in scaled},
+            "lite": {f.dataset: {**_provenance(f), "gt_abs_z": survey["datasets"][f.dataset]} for f in lite},
         },
         "created_at": datetime.now(UTC).isoformat(),
     }
     payload["sha256"] = payload_sha256(payload)
     return payload
 
 
+def _gt_abs_z_survey(
+    fits: list[DatasetFit], scalers: dict[str, dict[str, Any]], keep_mask: np.ndarray, kept_names: list[str]
+) -> dict[str, Any]:
+    """Survey |z| of every fit dataset's GT cells in its OWN test set's scaler (lite: the parent's).
+
+    Returns, over the FULL (non-lite) test sets -- lite cells are subsets of their
+    parent's, so pooling them would count cells twice -- the max |z|, the p99.99
+    (pooled over every GT cell and kept feature), where the max sits, the largest
+    enforced per-dataset GT clip fraction, the number of GT cells and of those with any
+    |z| > :data:`CP_Z_CLIP`; plus, for every dataset (lite included, unenforced),
+    ``{gt_clip_frac, enforced, max, p9999, max_feature}``.
+
+    Raises
+    ------
+    ValueError
+        If any non-lite dataset's ``gt_clip_frac`` (fraction of its GT cells with any
+        |z| > :data:`CP_Z_CLIP`) exceeds :data:`CP_GT_CLIP_FRAC_MAX`: the clip would
+        then discard more GT signal than the bound allows. Lite datasets are recorded
+        with ``enforced: false``: their cells are a subset of the parent's, in the
+        parent's scaler, so one extra cell on ~1000 must not brick every build.
+    """
+    per: dict[str, dict[str, Any]] = {}
+    pooled: list[np.ndarray] = []  # full (non-lite) datasets only: lite cells duplicate their parent's
+    n_full = n_full_beyond = 0
+    for fit in fits:
+        scaler = scalers[fit.parent or fit.dataset]
+        z = np.abs(
+            (np.asarray(fit.cells, dtype=np.float64)[:, keep_mask] - np.asarray(scaler["mean"]))
+            / np.asarray(scaler["std"])
+        )
+        col = int(np.unravel_index(np.argmax(z), z.shape)[1])
+        per[fit.dataset] = {
+            "gt_clip_frac": float((z > CP_Z_CLIP).any(axis=1).mean()),
+            # Lite GT is a subset of the parent's cells in the parent's scaler, so its fraction
+            # is small-sample noise on the parent's: recorded, never enforced.
+            "enforced": fit.parent is None,
+            "max": float(z.max()),
+            "p9999": float(np.quantile(z, 0.9999)),
+            "max_feature": kept_names[col],
+        }
+        if fit.parent is None:
+            pooled.append(z.ravel())
+            n_full += z.shape[0]
+            n_full_beyond += int((z > CP_Z_CLIP).any(axis=1).sum())
+    over = {d: r["gt_clip_frac"] for d, r in per.items() if r["enforced"] and r["gt_clip_frac"] > CP_GT_CLIP_FRAC_MAX}
+    if over:
+        raise ValueError(
+            f"GT clip fraction above the bound {CP_GT_CLIP_FRAC_MAX} at z_clip {CP_Z_CLIP} for {over}: the clip "
+            "would discard too much GT signal. Revisit CP_Z_CLIP (and re-survey) or investigate the GT cache."
+        )
+    worst = max((d for d in per if per[d]["enforced"]), key=lambda d: per[d]["max"])
+    return {
+        # Over FULL test sets only, so no GT cell is counted twice through a lite subset.
+        "max": per[worst]["max"],
+        "p9999": float(np.quantile(np.concatenate(pooled), 0.9999)),
+        "max_dataset": worst,
+        "max_feature": per[worst]["max_feature"],
+        "max_gt_clip_frac": max(r["gt_clip_frac"] for r in per.values() if r["enforced"]),
+        "n_cells": n_full,
+        "n_cells_beyond_clip": n_full_beyond,
+        "datasets": per,
+    }
+
+
 def write_cp_reference(payload: dict[str, Any], path: Path, *, force: bool = False) -> bool:
     """Write a reference payload atomically (tmp file in the same dir, then ``os.replace``).
 
@@ -621,7 +716,9 @@ def binding_sha256(self) -> str:
         """
         body = {
             "cp_identity": self.cp_identity,
-            "criteria": self.criteria,
+            # The GT |z| survey summarizes EVERY dataset of the target; it is hashed in the
+            # reference but left out here, so adding a dataset keeps this binding.
+            "criteria": {k: v for k, v in self.criteria.items() if k != "gt_abs_z_survey"},
             "feature_names": list(self.feature_names),
             "keep_mask": [bool(b) for b in self.keep_mask],
             "dataset": self.dataset,
@@ -658,6 +755,45 @@ def transform(self, x: np.ndarray) -> np.ndarray:
         """Return the kept columns of a raw CP matrix, standardized by this dataset's GT scaler."""
         return (self.select(x) - self.mean) / self.std
 
+    @property
+    def z_clip(self) -> float:
+        """The reference's shared-space clip, :data:`CP_Z_CLIP` at build time."""
+        return float(self.criteria["z_clip"])
+
+    def transform_clipped(self, x: np.ndarray) -> np.ndarray:
+        """Return :meth:`transform` clipped to ``+-z_clip``: the space KID/FID/cosine score in.
+
+        Applied identically to pred and GT, at dataset level and per row. The build
+        survey guarantees at most ``CP_GT_CLIP_FRAC_MAX`` of any dataset's GT cells are
+        clipped, so this mainly bounds how far one heavy-tailed predicted feature can
+        push the cubic KID kernel.
+        """
+        return np.clip(self.transform(x), -self.z_clip, self.z_clip)
+
+    def clip_fraction(self, x: np.ndarray) -> dict[str, Any]:
+        """Return what :meth:`transform_clipped` discards on ``x``: the fraction of cells clipped.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Raw ``(n_cells, n_features)`` CP matrix (the prediction's, in the eval).
+
+        Returns
+        -------
+        dict
+            ``{"any": fraction of cells with any kept-feature |z| > z_clip,
+            "per_feature": {kept feature: fraction of cells clipped on it}}``; NaN
+            fractions for zero cells.
+        """
+        kept_names = [n for n, k in zip(self.feature_names, self.keep_mask, strict=True) if k]
+        if x.shape[0] == 0:
+            return {"any": float("nan"), "per_feature": dict.fromkeys(kept_names, float("nan"))}
+        over = np.abs(self.transform(x)) > self.z_clip
+        return {
+            "any": float(over.any(axis=1).mean()),
+            "per_feature": {n: float(v) for n, v in zip(kept_names, over.mean(axis=0), strict=True)},
+        }
+
     def check_positions(self, positions: list[str]) -> None:
         """Refuse an eval whose GT positions are not exactly the ones the fit recorded.