diff --git a/applications/dynacell/src/dynacell/evaluation/README.md b/applications/dynacell/src/dynacell/evaluation/README.md index a8a2ba757..db704f777 100644 --- a/applications/dynacell/src/dynacell/evaluation/README.md +++ b/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/.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 `; `--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/.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 `; `--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. diff --git a/applications/dynacell/src/dynacell/evaluation/cp_reference.py b/applications/dynacell/src/dynacell/evaluation/cp_reference.py index 2b973e417..7dfc3e68f 100644 --- a/applications/dynacell/src/dynacell/evaluation/cp_reference.py +++ b/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,8 +552,8 @@ 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(), } @@ -530,6 +561,70 @@ def _provenance(fit: DatasetFit) -> dict[str, Any]: 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. diff --git a/applications/dynacell/src/dynacell/evaluation/cp_reference_test.py b/applications/dynacell/src/dynacell/evaluation/cp_reference_test.py index e4449df3a..995663267 100644 --- a/applications/dynacell/src/dynacell/evaluation/cp_reference_test.py +++ b/applications/dynacell/src/dynacell/evaluation/cp_reference_test.py @@ -12,6 +12,7 @@ from dynacell.evaluation import cp_reference from dynacell.evaluation.cache import StaleCacheError from dynacell.evaluation.cp_reference import ( + CP_Z_CLIP, STD_FLOOR_FRACTION, DatasetFit, fit_cp_reference, @@ -551,3 +552,128 @@ def test_old_sidecar_survives_a_rebuild_that_leaves_its_dataset_unchanged(tmp_pa _reference(tmp_path, [_fit("set-a", a * 1.1), _fit("set-b", b)]) # set-a's scaler moves with pytest.raises(ValueError, match="CP space of set-a .* changed since the eval"): sidecar_cp_space(eval_dir) + + +def test_build_records_the_z_clip_and_the_gt_survey(two_sets) -> None: + """The clip and the GT |z| survey that justifies it are hashed criteria; every dataset records its own.""" + ref, payload, cells = two_sets + criteria = payload["criteria"] + assert criteria["z_clip"] == CP_Z_CLIP + survey = criteria["gt_abs_z_survey"] + z = np.abs(ref.for_dataset("set-b").transform(cells["set-b"])) + assert payload["fit"]["datasets"]["set-b"]["gt_abs_z"]["max"] == pytest.approx(z.max()) + assert survey["max"] == max( + r["gt_abs_z"]["max"] for r in [*payload["fit"]["datasets"].values(), *payload["fit"]["lite"].values()] + ) + assert survey["max"] < CP_Z_CLIP + edited = json.loads(json.dumps(payload)) + edited["criteria"]["z_clip"] = 1e9 + assert payload_sha256(edited) != payload["sha256"] + + +def test_build_refuses_a_dataset_whose_gt_clip_fraction_exceeds_the_bound(monkeypatch) -> None: + """If more than CP_GT_CLIP_FRAC_MAX of one dataset's GT cells would be clipped, the reference is not built.""" + monkeypatch.setattr(cp_reference, "CP_Z_CLIP", 2.0) # many GT cells of a normal sample exceed |z| = 2 + with pytest.raises(ValueError, match=r"GT clip fraction above the bound 0\.001 at z_clip 2\.0 for \{'set-a'"): + fit_cp_reference([_fit("set-a", _gt())], target_name="er", feature_names=_NAMES, cp_identity={}) + + +def test_build_accepts_a_few_gt_cells_beyond_the_clip_under_the_bound() -> None: + """One GT cell beyond +-z_clip in 5000 (2e-4 <= 1e-3) is allowed, and recorded per dataset.""" + gt = _gt(n=5000) + gt[0, 0] += 1e3 # one extreme cell on f0 (kept): |z| far above the clip + payload = fit_cp_reference([_fit("set-a", gt)], target_name="er", feature_names=_NAMES, cp_identity={}) + rec = payload["fit"]["datasets"]["set-a"]["gt_abs_z"] + assert rec["gt_clip_frac"] == pytest.approx(1 / 5000) + assert rec["max"] > CP_Z_CLIP and rec["max_feature"] == "f0" + criteria = payload["criteria"] + assert criteria["gt_clip_frac_max"] == 1e-3 and "c >= 20" in criteria["rank_stability"] + assert criteria["gt_abs_z_survey"]["max_gt_clip_frac"] == pytest.approx(1 / 5000) + + +def test_lite_clip_fraction_is_recorded_but_not_enforced(monkeypatch) -> None: + """A lite set over the bound builds (recorded, ``enforced: false``); its parent stays enforced. + + The exemption is not "lite is unchecked". Lite GT cells are a subset of the + parent's cells, scored in the parent's scaler, so a lite set's clipped-cell COUNT + is at most its parent's. The enforced parent bound therefore covers the lite set; + only its FRACTION can exceed the bound, as small-sample noise on fewer cells. Here + one parent outlier is 1/5000 of the parent (under the bound) and 1/500 of the lite + set (over it): the build must pass, and must refuse once the parent itself is over. + """ + parent = _gt(n=5000) + parent[0, 0] += 1e3 # 1 of 5000 parent cells beyond the clip: 2e-4, under the bound + lite_cells = parent[:500] # the same outlier in 500 lite cells: 2e-3, over the bound + fits = [_fit("set-a", parent), _fit("set-a-lite", lite_cells, in_mask_fit=False, parent="set-a")] + payload = fit_cp_reference(fits, target_name="er", feature_names=_NAMES, cp_identity={}) + lite = payload["fit"]["lite"]["set-a-lite"]["gt_abs_z"] + assert lite["gt_clip_frac"] == pytest.approx(1 / 500) + assert lite["gt_clip_frac"] > cp_reference.CP_GT_CLIP_FRAC_MAX + assert lite["enforced"] is False and payload["fit"]["datasets"]["set-a"]["gt_abs_z"]["enforced"] is True + survey = payload["criteria"]["gt_abs_z_survey"] + assert survey["max_gt_clip_frac"] == pytest.approx(1 / 5000) + # Full test sets only: the lite set's copy of the parent's outlier is not counted twice. + assert (survey["n_cells"], survey["n_cells_beyond_clip"]) == (5000, 1) + + monkeypatch.setattr(cp_reference, "CP_GT_CLIP_FRAC_MAX", 1e-4) # now the parent itself is over + with pytest.raises(ValueError, match=r"for \{'set-a': 0\.0002\}"): + fit_cp_reference(fits, target_name="er", feature_names=_NAMES, cp_identity={}) + + +def _space_with_gt_beyond_the_clip(tmp_path: Path): + """A reference whose own GT holds one cell beyond +-z_clip (1 of 2000, under the bound); returns (space, gt).""" + gt = _gt(n=2000) + gt[0, 0] += 1e3 # one GT cell far out on f0 (kept) + ref, _ = _reference(tmp_path, [_fit("set-a", gt)]) + space = ref.for_dataset("set-a") + assert np.abs(space.transform(gt)).max() > space.z_clip + return space, gt + + +def test_staged_arrays_clip_both_pred_and_gt(tmp_path: Path) -> None: + """Dataset-level metric arrays are ``transform_clipped`` of BOTH sides, with values beyond the clip on both.""" + space, gt = _space_with_gt_beyond_the_clip(tmp_path) + pred = gt.copy() + pred[5, 0] += 1e4 # a pred-only outlier too + staged = _stage_cp_dataset_inputs(_lists(pred, gt), space, _blocks(gt), tmp_path) + np.testing.assert_array_equal(staged[1], space.transform_clipped(pred)) + np.testing.assert_array_equal(staged[2], space.transform_clipped(gt)) + assert not np.array_equal(staged[2], space.transform(gt)) # the GT clip is really active + assert np.abs(staged[1]).max() == np.abs(staged[2]).max() == space.z_clip + + +def test_row_features_clip_both_pred_and_gt(tmp_path: Path) -> None: + """Per-row CP inputs are ``transform_clipped`` of BOTH sides, with values beyond the clip on both.""" + space, gt = _space_with_gt_beyond_the_clip(tmp_path) + pred = gt.copy() + pred[5, 0] += 1e4 + p, g = _cp_row_features(pred[:10], gt[:10], space) + np.testing.assert_array_equal(p, space.transform_clipped(pred[:10])) + np.testing.assert_array_equal(g, space.transform_clipped(gt[:10])) + assert np.abs(g).max() == space.z_clip and np.abs(p).max() == space.z_clip + + +def test_clip_and_its_bound_are_part_of_every_binding(tmp_path: Path) -> None: + """Changing z_clip or the GT clip-fraction bound moves ``binding_sha256`` (the eval stamp) of every dataset.""" + fits = [ + _fit("set-a", _gt()), + _fit("set-b", _gt(seed=1)), + _fit("set-a-lite", _gt()[:100], in_mask_fit=False, parent="set-a"), + ] + ref, payload = _reference(tmp_path, fits) + names = ("set-a", "set-b", "set-a-lite") + before = {d: ref.for_dataset(d).binding_sha256 for d in names} + for key, value in (("z_clip", 25.0), ("gt_clip_frac_max", 2e-3)): + moved = json.loads(json.dumps(payload)) + moved["criteria"][key] = value + moved["sha256"] = payload_sha256(moved) + (tmp_path / "er.json").write_text(json.dumps(moved)) + again = load_cp_reference(tmp_path / "er.json", target_name="er") + assert all(again.for_dataset(d).binding_sha256 != before[d] for d in names), key + + +def test_empty_feature_mask_is_refused() -> None: + """A GT pool whose every CP column is constant keeps no feature: the build refuses with a clear error.""" + constant = np.ones((64, len(_NAMES))) + with pytest.raises(ValueError, match="kept no CP features"): + fit_cp_reference([_fit("set-a", constant)], target_name="er", feature_names=_NAMES, cp_identity={}) diff --git a/applications/dynacell/src/dynacell/evaluation/pipeline.py b/applications/dynacell/src/dynacell/evaluation/pipeline.py index a50dca4de..35727d2f2 100644 --- a/applications/dynacell/src/dynacell/evaluation/pipeline.py +++ b/applications/dynacell/src/dynacell/evaluation/pipeline.py @@ -104,14 +104,15 @@ def _build_focus_slabs_map(config: DictConfig, gt_positions) -> dict[str, list[s def _cp_row_features(pred_cp: np.ndarray, gt_cp: np.ndarray, cp_space: DatasetCPSpace) -> tuple[np.ndarray, np.ndarray]: - """Per-(FOV, timepoint) CP inputs: the reference mask + this dataset's GT scaler on both sides. + """Per-(FOV, timepoint) CP inputs: the reference mask + this dataset's GT scaler, clipped, on both sides. - The same transform as the dataset-level stage, so a per-row CP value is + The same transform as the dataset-level stage + (:meth:`DatasetCPSpace.transform_clipped`), so a per-row CP value is comparable across models and FOVs. Empty inputs (no cells) pass through. """ if not (pred_cp.size and gt_cp.size): return pred_cp, gt_cp - return cp_space.transform(pred_cp), cp_space.transform(gt_cp) + return cp_space.transform_clipped(pred_cp), cp_space.transform_clipped(gt_cp) def _real_vs_pred_probe( @@ -334,7 +335,10 @@ def _extend_backbone( def _gate_and_record_cp_space( - cp_space: DatasetCPSpace, gt_cp_blocks: dict[tuple[str, int], np.ndarray], save_dir: Path + cp_space: DatasetCPSpace, + gt_cp_blocks: dict[tuple[str, int], np.ndarray], + save_dir: Path, + pred_clip: dict[str, Any], ) -> None: """Run the GT content gate and write the CP sidecar for this eval dir. @@ -350,9 +354,25 @@ def _gate_and_record_cp_space( ``{(position, t): GT-finite CP rows}`` the run scored. save_dir : pathlib.Path Eval dir receiving the sidecar. + pred_clip : dict + :meth:`DatasetCPSpace.clip_fraction` of the run's pred CP cells, recorded in + the sidecar so what the z-clip discarded stays inspectable per feature. """ cp_space.check_gt_cells(gt_cp_blocks) - (save_dir / CP_SIDECAR_FILENAME).write_text(json.dumps(cp_sidecar_payload(cp_space), indent=2)) + + def _json_frac(v: float) -> float | None: + # No pred CP rows -> NaN fraction; strict JSON has no NaN, so record null. + return None if np.isnan(v) else v + + payload = { + **cp_sidecar_payload(cp_space), + "clip": { + "z_clip": cp_space.z_clip, + "pred_clip_frac": _json_frac(pred_clip["any"]), + "pred_clip_frac_per_feature": {k: _json_frac(v) for k, v in pred_clip["per_feature"].items()}, + }, + } + (save_dir / CP_SIDECAR_FILENAME).write_text(json.dumps(payload, indent=2, allow_nan=False)) def _stage_cp_dataset_inputs( @@ -360,11 +380,12 @@ def _stage_cp_dataset_inputs( cp_space: DatasetCPSpace, gt_cp_blocks: dict[tuple[str, int], np.ndarray], save_dir: Path, + pred_clip: dict[str, Any] | None = None, ) -> tuple[str, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """Stage the dataset-level CP inputs in the reference feature space. KID/FID/cosine get the reference-masked features standardized by this - dataset's GT scaler (identical for pred and GT); the linear probe gets the + dataset's GT scaler and clipped to +-z_clip (identical for pred and GT); the linear probe gets the reference-masked raw features, since it fits its own per-fold scaler. Writes ``cp_selected_feature_mask.json`` recording which reference and scaler were applied. @@ -381,19 +402,24 @@ def _stage_cp_dataset_inputs( (:meth:`DatasetCPSpace.check_gt_cells`) before anything is written. save_dir : pathlib.Path Eval dir receiving the sidecar. + pred_clip : dict, optional + :meth:`DatasetCPSpace.clip_fraction` of the pred cells, if the caller already + computed it; computed here otherwise. Returns ------- tuple ``("CP", pred_metric, target_metric, pred_probe, target_probe, pred_fovs, target_fovs)``. """ - _gate_and_record_cp_space(cp_space, gt_cp_blocks, save_dir) pred_cp_raw = np.concatenate(cp.pred_feats, axis=0) + _gate_and_record_cp_space( + cp_space, gt_cp_blocks, save_dir, pred_clip if pred_clip is not None else cp_space.clip_fraction(pred_cp_raw) + ) target_cp_raw = np.concatenate(cp.gt_feats, axis=0) return ( "CP", - cp_space.transform(pred_cp_raw), - cp_space.transform(target_cp_raw), + cp_space.transform_clipped(pred_cp_raw), + cp_space.transform_clipped(target_cp_raw), cp_space.select(pred_cp_raw), cp_space.select(target_cp_raw), np.concatenate(cp.pred_fovs, axis=0), @@ -978,6 +1004,7 @@ def _nuclei_seeds(segment_fn, **infer) -> np.ndarray | None: gt_cp_blocks[t] = gt_cp_per_t[t][np.isfinite(gt_cp_per_t[t]).all(axis=1)] pred_cp, gt_cp_t = drop_paired_nonfinite_rows(pred_cp, gt_cp_per_t[t]) pred_cp_z, gt_cp_z = _cp_row_features(pred_cp, gt_cp_t, cp_space) + cp_clip_frac = cp_space.clip_fraction(pred_cp)["any"] if pred_cp.size else float("nan") # Prefixes come from _COLUMN_PREFIX, not literals: these per-timepoint # columns must match the dataset-level ones the parent derives from the # same map, or the two halves of a _* family drift apart. @@ -985,6 +1012,7 @@ def _nuclei_seeds(segment_fn, **infer) -> np.ndarray | None: **compute_feature_similarity_pairwise( pred_cp_z, gt_cp_z, _COLUMN_PREFIX["cp"], compute_fid=compute_fid ), + f"{_COLUMN_PREFIX['cp']}_clip_frac": cp_clip_frac, **compute_feature_similarity_pairwise( pred_dinov3, gt_dinov3_per_t[t], _COLUMN_PREFIX["dinov3"], compute_fid=compute_fid ), @@ -1784,10 +1812,21 @@ def _aggregate(result: FovResult) -> None: # MADScaler can normalize per-fold. prefix_inputs: list[tuple[str, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]] = [] + # What the shared-space z-clip discards: the fraction of pred CP cells with any + # kept feature beyond +-z_clip (NaN without pred cells), next to Dataset_CP_KID. + pred_cp_all = ( + np.concatenate(parent_lists["cp"].pred_feats, axis=0) + if parent_lists["cp"].pred_feats + else np.empty((0, len(cp_space.feature_names))) + ) + pred_clip = cp_space.clip_fraction(pred_cp_all) + dataset_row[f"Dataset_{_COLUMN_PREFIX['cp']}_clip_frac"] = pred_clip["any"] if parent_lists["cp"].pred_feats: - prefix_inputs.append(_stage_cp_dataset_inputs(parent_lists["cp"], cp_space, gt_cp_blocks, save_dir)) + prefix_inputs.append( + _stage_cp_dataset_inputs(parent_lists["cp"], cp_space, gt_cp_blocks, save_dir, pred_clip) + ) else: # no finite pred CP rows: CP metrics are NaN, but the GT gate and sidecar still apply - _gate_and_record_cp_space(cp_space, gt_cp_blocks, save_dir) + _gate_and_record_cp_space(cp_space, gt_cp_blocks, save_dir, pred_clip) for key in deep_kinds: # cp is handled above (reference space) bb = parent_lists[key] @@ -2019,6 +2058,11 @@ def _final_metrics_cache_valid(config: DictConfig) -> bool: ] if any(f"Dataset_{prefix}_{suffix}" not in columns for prefix in prefixes for suffix in suffixes): return False + # CP rows scored before the shared-space z-clip carry no clip-fraction columns: + # their KID was computed unclipped, so they must be recomputed, not reused. + cp = _COLUMN_PREFIX["cp"] + if cp in prefixes and not {f"Dataset_{cp}_clip_frac", f"{cp}_clip_frac"} <= set(columns): + return False return pixel_ok and mask_ok and feature_ok diff --git a/applications/dynacell/tests/test_evaluation_cp_end_to_end.py b/applications/dynacell/tests/test_evaluation_cp_end_to_end.py index ba18ba140..1f040db1b 100644 --- a/applications/dynacell/tests/test_evaluation_cp_end_to_end.py +++ b/applications/dynacell/tests/test_evaluation_cp_end_to_end.py @@ -14,6 +14,7 @@ from __future__ import annotations +import json import tempfile import zlib from pathlib import Path @@ -312,7 +313,13 @@ def nan_pred_side(ctx, pos_name, image, cell_segmentation): monkeypatch.setattr(harness.pipeline, "fov_cp_features", nan_pred_side) row, _ = harness.run("nanpred", [g.copy() for g in harness.gt]) assert np.isnan(row["Dataset_CP_KID"]) - assert (harness.root / "nanpred" / "cp_selected_feature_mask.json").exists() + sidecar_text = (harness.root / "nanpred" / "cp_selected_feature_mask.json").read_text() + + def _reject(const: str) -> None: + raise ValueError(f"non-strict JSON constant {const}") + + sidecar = json.loads(sidecar_text, parse_constant=_reject) # strict JSON: no NaN/Infinity + assert sidecar["clip"]["pred_clip_frac"] is None # A GT recompute from a changed GT store bypasses the up-front cache check, so only the # post-staging gate can refuse it -- and it must, even with no finite pred CP rows. @@ -322,3 +329,57 @@ def nan_pred_side(ctx, pos_name, image, cell_segmentation): with pytest.raises(Exception, match="differ from the CP reference fit") as err: harness.run("nanpred_moved", [g.copy() for g in harness.gt], gt_cp=True) assert type(err.value).__name__ == "StaleCacheError" + + +def test_clip_fraction_is_reported_per_dataset_row_and_feature(harness: Harness) -> None: + """Dataset_CP_clip_frac, per-row CP_clip_frac and the sidecar's per-feature fractions match a direct count. + + The prediction replaces 4 of the 16 cells with a voxel checkerboard, whose texture + features lie far outside the GT's, so only those cells can leave the +-z_clip band. + """ + import json + + blown = np.isin(harness.seg, [1, 2, 3, 4]) + checker = (np.indices(harness.seg.shape).sum(axis=0) % 2).astype(np.float32) # max-texture voxel pattern + pred = [np.where(blown, checker, g).astype(np.float32) for g in harness.gt] + row, config = harness.run("clipped", pred) + space = harness.pipeline.eval_cp_space(config) + save_dir = Path(config.save.save_dir) + emb = np.load(save_dir / "embeddings" / "pred_cp_single_cell_embeddings.npz") + over = np.abs(space.transform(emb["embeddings"])) > space.z_clip + assert 0 < row["Dataset_CP_clip_frac"] == pytest.approx(over.any(axis=1).mean()) + sidecar = json.loads((save_dir / "cp_selected_feature_mask.json").read_text())["clip"] + assert sidecar["z_clip"] == space.z_clip and sidecar["pred_clip_frac"] == row["Dataset_CP_clip_frac"] + kept = [n for n, k in zip(space.feature_names, space.keep_mask, strict=True) if k] + assert sidecar["pred_clip_frac_per_feature"] == pytest.approx(dict(zip(kept, over.mean(axis=0), strict=True))) + _, _, rows = harness.pipeline.evaluate_predictions( + config, models=_stub_models(), cp_space=space + ) # per-row fractions: recomputed per (FOV, t) from that row's pred cells + per_row = {(r["FOV"], r["Timepoint"]): r["CP_clip_frac"] for r in rows} + for (fov, t), frac in per_row.items(): + sel = (emb["fov"] == fov) & (emb["timepoint"] == t) + assert frac == pytest.approx(over[sel].any(axis=1).mean()) + + # KID is computed on the CLIPPED z, dataset level and per row: it equals KID on + # transform_clipped, and the unclipped value (1e4x larger here) would dominate. + from dynacell.evaluation.feature_metrics import compute_feature_similarity, compute_feature_similarity_pairwise + + gt_emb = np.load(save_dir / "embeddings" / "gt_cp_single_cell_embeddings.npz")["embeddings"] + flags = {"compute_fid": False, "compute_prc": False, "compute_mind": False} + clipped = compute_feature_similarity( + space.transform_clipped(emb["embeddings"]), space.transform_clipped(gt_emb), "CP", **flags + )["CP_KID"] + unclipped = compute_feature_similarity(space.transform(emb["embeddings"]), space.transform(gt_emb), "CP", **flags)[ + "CP_KID" + ] + assert row["Dataset_CP_KID"] == pytest.approx(clipped, rel=1e-9) + assert abs(unclipped) > 1e3 * abs(clipped) + for r in rows: + sel = (emb["fov"] == r["FOV"]) & (emb["timepoint"] == r["Timepoint"]) + expected = compute_feature_similarity_pairwise( + space.transform_clipped(emb["embeddings"][sel]), + space.transform_clipped(gt_emb[sel]), + "CP", + compute_fid=False, + )["CP_KID"] + assert r["CP_KID"] == pytest.approx(expected, rel=1e-9) diff --git a/applications/dynacell/tests/test_evaluation_timepoint_guard.py b/applications/dynacell/tests/test_evaluation_timepoint_guard.py index 42310b816..13e6e17b0 100644 --- a/applications/dynacell/tests/test_evaluation_timepoint_guard.py +++ b/applications/dynacell/tests/test_evaluation_timepoint_guard.py @@ -124,8 +124,8 @@ def _write_final_caches(pipeline, save_dir: Path, feature_row: dict, config) -> def _dataset_row(prefixes: tuple[str, ...], families: tuple[str, ...]) -> dict: - """``Dataset__`` columns plus the always-present KID/cosine.""" - row: dict[str, float] = {"CP_KID": 0.0} + """``Dataset__`` columns plus the always-present KID/cosine (and CP's clip fractions).""" + row: dict[str, float] = {"CP_KID": 0.0, "CP_clip_frac": 0.0, "Dataset_CP_clip_frac": 0.0} for prefix in prefixes: for family in ("KID", "KID_std", "Median_Cosine_Similarity", *families): row[f"Dataset_{prefix}_{family}"] = 0.0 @@ -165,10 +165,25 @@ def test_real_full_benchmark_columns_stay_valid(tmp_path: Path) -> None: header = f.readline().rstrip("\n").split(",") pipeline = live_pipeline_module() config, _ = _feature_cache_config(tmp_path) - _write_final_caches(pipeline, tmp_path, dict.fromkeys(header, 0.0), config) + # A real pre-clip dir lacks the clip-fraction columns (see the next test); add them to + # check that the rest of the real column set still passes. + _write_final_caches( + pipeline, tmp_path, dict.fromkeys([*header, "CP_clip_frac", "Dataset_CP_clip_frac"], 0.0), config + ) assert pipeline._final_metrics_cache_valid(config) +def test_cache_scored_before_the_z_clip_is_invalid(tmp_path: Path) -> None: + """A feature cache without the CP clip-fraction columns was scored unclipped and must be recomputed.""" + pipeline = live_pipeline_module() + config, _ = _feature_cache_config(tmp_path) + for missing in ("Dataset_CP_clip_frac", "CP_clip_frac"): + row = _dataset_row(("CP", "DINOv3", "DynaCLR"), _ALL_FAMILIES) + del row[missing] + _write_final_caches(pipeline, tmp_path, row, config) + assert not pipeline._final_metrics_cache_valid(config), missing + + def test_cache_scored_in_another_cp_reference_is_invalid(tmp_path: Path) -> None: """Rebuilding the CP reference invalidates every final-metrics cache stamped with the old one.""" pipeline = live_pipeline_module() diff --git a/applications/dynacell/tools/lite_checks.py b/applications/dynacell/tools/lite_checks.py index a03cad7f4..e7c7d36e8 100644 --- a/applications/dynacell/tools/lite_checks.py +++ b/applications/dynacell/tools/lite_checks.py @@ -205,9 +205,9 @@ def check_e( Y = np.asarray(np.load(g, allow_pickle=True)["embeddings"], dtype=np.float64) keep = np.isfinite(X).all(1) & np.isfinite(Y).all(1) X, Y = X[keep], Y[keep] - if tok == "cp": # the pipeline's CP space: reference mask + this dataset's GT scaler + if tok == "cp": # the pipeline's CP space: reference mask + dataset GT scaler, clipped space = sidecar_cp_space(path) - X, Y = space.transform(X), space.transform(Y) + X, Y = space.transform_clipped(X), space.transform_clipped(Y) n_all = X.shape[0] if tok != "cp": # deep extractors: precompute cell-level kernels once, slice per draw K_XX, K_YY, K_XY = poly_kernel(X, X), poly_kernel(Y, Y), poly_kernel(X, Y) diff --git a/applications/dynacell/tools/lite_checks_test.py b/applications/dynacell/tools/lite_checks_test.py index 34b97f6dc..1efb2ad2f 100644 --- a/applications/dynacell/tools/lite_checks_test.py +++ b/applications/dynacell/tools/lite_checks_test.py @@ -121,6 +121,7 @@ def test_check_e_scores_cp_in_the_reference_space(tmp_path: Path) -> None: (path / "embeddings").mkdir(parents=True) gt = rng.normal(size=(n, d)) pred = 0.5 * gt + 0.2 * rng.normal(size=gt.shape) + pred[:3, 0] += 1e3 # a few pred cells far beyond the z-clip, so an unclipped check E diverges fov = np.array([f"0/0/fov{i % 6:04d}" for i in range(n)]) for side, arr in (("gt", gt), ("pred", pred)): np.savez(path / "embeddings" / f"{side}_cp_single_cell_embeddings.npz", embeddings=arr, fov=fov) @@ -132,11 +133,12 @@ def test_check_e_scores_cp_in_the_reference_space(tmp_path: Path) -> None: mask = np.array(payload["keep_mask"]) mean, std = np.array(payload["scalers"]["ds"]["mean"]), np.array(payload["scalers"]["ds"]["std"]) - expected = mmd2_from_features( - poly3_features((pred[:, mask] - mean) / std), poly3_features((gt[:, mask] - mean) / std) - ) + c = payload["criteria"]["z_clip"] + zp, zg = (pred[:, mask] - mean) / std, (gt[:, mask] - mean) / std + expected = mmd2_from_features(poly3_features(np.clip(zp, -c, c)), poly3_features(np.clip(zg, -c, c))) full = result.loc[result.extractor == "CP", "full"].iloc[0] assert full == pytest.approx(expected, rel=1e-10) + assert abs(mmd2_from_features(poly3_features(zp), poly3_features(zg))) > 10 * abs(full) # unclipped diverges def zscore(x: np.ndarray) -> np.ndarray: return (x - x.mean(0)) / (x.std(0) + 1e-8) diff --git a/applications/dynacell/tools/lite_subset_sim.py b/applications/dynacell/tools/lite_subset_sim.py index cdfdc5ff8..d7e5ceac2 100644 --- a/applications/dynacell/tools/lite_subset_sim.py +++ b/applications/dynacell/tools/lite_subset_sim.py @@ -16,7 +16,7 @@ equals the pipeline's ``100 x min(1000, n)`` subset mean exactly when n <= 1000 and in expectation above that. GLCM+ (``cp``) is scored in the pipeline's CP space: the target's CP reference mask plus the eval dataset's GT scaler (the parent's for a -lite dataset), both read via the run's ``cp_selected_feature_mask.json`` sidecar. +lite dataset), clipped to the reference's z_clip, both read via the run's ``cp_selected_feature_mask.json`` sidecar. That transform is fixed, so it does not depend on the cell subset; CP is still scored directly through an exact degree-3 feature map (22-dim: cheap). @@ -193,9 +193,9 @@ def __init__(self, path: Path, blocks: list[tuple[str, int]]): ok = cell_block >= 0 X, Y, cell_block = X[ok], Y[ok], cell_block[ok] if tok == "cp": - # pipeline: the CP reference mask + this dataset's GT scaler, on both sides. + # pipeline: the CP reference mask + this dataset's GT scaler, clipped, on both sides. space = sidecar_cp_space(path) - self.cp = dict(X=space.transform(X), Y=space.transform(Y), cell_block=cell_block) + self.cp = dict(X=space.transform_clipped(X), Y=space.transform_clipped(Y), cell_block=cell_block) continue # one-hot block membership (n x nb) -> block sums via M^T K M M = np.zeros((len(cell_block), nb)) diff --git a/applications/dynacell/tools/lite_subset_sim_test.py b/applications/dynacell/tools/lite_subset_sim_test.py index 14ab2ac05..8b84400c8 100644 --- a/applications/dynacell/tools/lite_subset_sim_test.py +++ b/applications/dynacell/tools/lite_subset_sim_test.py @@ -60,10 +60,12 @@ def _write_reference(path: Path) -> dict: return payload -def _reference_transform(payload: dict, x: np.ndarray) -> np.ndarray: - """Apply the pipeline's CP transform spelled out from the payload: mask, then this dataset's scaler.""" +def _reference_transform(payload: dict, x: np.ndarray, clip: bool = True) -> np.ndarray: + """Apply the pipeline's CP transform spelled out from the payload: mask, this dataset's scaler, then the z-clip.""" scaler = payload["scalers"][_DATASET] - return (x[:, np.array(payload["keep_mask"])] - np.array(scaler["mean"])) / np.array(scaler["std"]) + z = (x[:, np.array(payload["keep_mask"])] - np.array(scaler["mean"])) / np.array(scaler["std"]) + c = payload["criteria"]["z_clip"] + return np.clip(z, -c, c) if clip else z def _pipeline_kid(pred: np.ndarray, gt: np.ndarray) -> float: @@ -91,6 +93,8 @@ def _write_eval_dir(path: Path, seed: int, reference: Path) -> dict[str, dict[st for tok, d in (("dinov3", _D_DEEP), ("cp", _D_CP)): gt = rng.normal(size=(len(fov), d)) pred = 1.5 * gt + 0.8 + 0.3 * rng.normal(size=gt.shape) + if tok == "cp": + pred[:4, 0] += 1e3 # a few pred cells far beyond the z-clip, so an unclipped tool diverges for side, arr in (("gt", gt), ("pred", pred)): np.savez( path / "embeddings" / f"{side}_{tok}_single_cell_embeddings.npz", embeddings=arr, fov=fov, timepoint=tp @@ -155,6 +159,11 @@ def test_block_sum_kid_matches_pipeline_kid(system, sel: np.ndarray) -> None: mask = np.array(payload["keep_mask"]) old = _pipeline_kid(_per_side_zscore(raw_pred[:, mask]), _per_side_zscore(raw_gt[:, mask])) assert out["CP_KID"] != pytest.approx(old, rel=0.1) + if sel.tolist() == list(range(len(_BLOCKS))): # the outlier cells are in block 0: only the full set holds them + unclipped = _pipeline_kid( + _reference_transform(payload, raw_pred, clip=False), _reference_transform(payload, raw_gt, clip=False) + ) + assert abs(unclipped) > 10 * abs(out["CP_KID"]) def test_multiset_weights_equal_brute_force_duplication(system) -> None: