Blinding: refuse tomography at the producers, check stamps on consumption, test the guarantees

src/sp_validation/blinding.py

@@ -77,7 +77,7 @@ def _shear_file(entry):
     shear = entry.get("shear") if isinstance(entry, dict) else None
     if not isinstance(shear, dict) or "path" not in shear:
         return None
-    return os.path.normpath(os.path.join(str(entry.get("subdir", "")), shear["path"]))
+    return os.path.realpath(os.path.join(str(entry.get("subdir", "")), shear["path"]))
 
 
 def blind_of(catalogues, entry):

src/sp_validation/cosmo_val/pseudo_cl.py

@@ -213,6 +213,8 @@ def calculate_pseudo_cl(self, compute_tomography=True, out_path=None):
         self.print_done("Done pseudo-Cl's")
 
     def calculate_pseudo_cl_map(self, ver, nside, out_path, tomo_bin_a, tomo_bin_b):
+        if (tomo_bin_a, tomo_bin_b) != ("all", "all"):
+            self._refuse_blinded_tomography(ver)
         assert (tomo_bin_a == "all" and tomo_bin_b == "all") or (
             isinstance(tomo_bin_a, (int, np.integer))
             and isinstance(tomo_bin_b, (int, np.integer))
@@ -312,6 +314,8 @@ def calculate_pseudo_cl_map(self, ver, nside, out_path, tomo_bin_a, tomo_bin_b):
         )
 
     def calculate_pseudo_cl_catalog(self, ver, out_path, tomo_bin_a, tomo_bin_b):
+        if (tomo_bin_a, tomo_bin_b) != ("all", "all"):
+            self._refuse_blinded_tomography(ver)
         assert (tomo_bin_a == "all" and tomo_bin_b == "all") or (
             isinstance(tomo_bin_a, (int, np.integer))
             and isinstance(tomo_bin_b, (int, np.integer))

src/sp_validation/cosmo_val/real_space.py

@@ -81,6 +81,7 @@ def calculate_2pcf_version(
         ggs = {f"tomo_bin_{b1}_tomo_bin_{b2}": None for b1, b2 in tomo_bin_pairs}
         held = self.xi_parts.get((ver, grid)) if reuse else None
         if ALL_PAIR in ggs and held is not None:
+            sacc_io.expect(held, blind)
             ggs[ALL_PAIR] = sacc_io.xi_correlation(held)
         to_compute = [
             (b1, b2)

src/sp_validation/sacc_io.py

@@ -983,6 +983,15 @@ def stamp(s):
     return s.metadata.get(STAMP_KEY, "none")
 
 
+def expect(part, blind):
+    """Refuse a part from a different blind before consuming its signal."""
+    name = blind if isinstance(blind, str) else blind.name
+    if stamp(part) != name:
+        raise ValueError(
+            f"part is stamped {stamp(part)}, but the catalogue is declared under {name}"
+        )
+
+
 def seal(s, blind, theory=None):
     """What :func:`save` writes for a birth, kept in memory: a copy of ``s``
     concealed under ``blind`` with ``theory`` and stamped with its name."""
@@ -1023,13 +1032,10 @@ def save(s, path, *, blind=None, theory=None, derived_from=None):
         if len(stamps) != 1:
             raise ValueError(f"input parts carry different blinds: {sorted(stamps)}")
         (name,) = stamps
-        if blind is not None and blind.name != name:
-            raise ValueError(
-                f"parts are stamped {name}, but the catalogue is declared under "
-                f"{blind.name}"
-            )
         out = s.copy()
         out.metadata[STAMP_KEY] = name
+        if blind is not None:
+            expect(out, blind)
     elif blind is not None:
         out = seal(s, blind, theory)
     else:
@@ -1038,9 +1044,12 @@ def save(s, path, *, blind=None, theory=None, derived_from=None):
     return out
 
 
-def load(path):
-    """Load the SACC at ``path``."""
-    return sacc.Sacc.load_fits(str(path))
+def load(path, blind=None):
+    """Load a SACC, checking its stamp against ``blind`` when supplied."""
+    part = sacc.Sacc.load_fits(str(path))
+    if blind is not None:
+        expect(part, blind)
+    return part
 
 
 # =============================================================================

src/sp_validation/tests/test_blind_consumers.py

@@ -0,0 +1,85 @@
+"""Check parts against the declared blind before derivations or reporting."""
+
+import runpy
+from pathlib import Path
+from types import SimpleNamespace
+
+import pytest
+
+from sp_validation import blinding, sacc_io
+
+SCRIPTS = Path(__file__).resolve().parents[3] / "workflow" / "scripts"
+
+
+@pytest.fixture
+def part_path(tmp_path):
+    part = sacc_io.new_sacc({0: ([0.1, 1.0], [1.0, 1.0])})
+    sacc_io.add_xi(part, (0, 0), [10.0], [1.0], [0.5], grid="integration")
+    path = tmp_path / "xi.sacc"
+    sacc_io.save(part, path, blind=blinding.NONE)
+    return path
+
+
+@pytest.mark.parametrize("expected", ["none", "toy"])
+def test_reporting_parts(part_path, monkeypatch, expected):
+    monkeypatch.syspath_prepend(str(SCRIPTS))
+    from cv_runner import take_reporting_parts
+
+    cv = SimpleNamespace(
+        versions=["v"], xi_parts={}, blind=lambda v: blinding.Blind(expected)
+    )
+    if expected == "toy":
+        with pytest.raises(ValueError, match="stamped none.*toy"):
+            take_reporting_parts(cv, [part_path])
+        assert not cv.xi_parts
+    else:
+        take_reporting_parts(cv, [part_path])
+        assert sacc_io.stamp(cv.xi_parts["v", "reporting"]) == "none"
+
+
+@pytest.mark.parametrize(
+    "script, input_key, reader",
+    [
+        ("cv_cosebis.py", "xi", "get_xi"),
+        ("cv_pure_eb.py", "xi_integration", "xi_correlation"),
+    ],
+)
+@pytest.mark.parametrize("expected", ["none", "toy"])
+def test_derived_script_checks_before_reading_signal(
+    part_path, monkeypatch, script, input_key, reader, expected
+):
+    monkeypatch.syspath_prepend(str(SCRIPTS))
+    import cv_runner
+
+    monkeypatch.setattr(cv_runner, "_unbuffer_streams", lambda: None)
+
+    class ReachedSignal(Exception):
+        pass
+
+    def read_signal(*args, **kwargs):
+        raise ReachedSignal
+
+    monkeypatch.setattr(sacc_io, reader, read_signal)
+    smk = SimpleNamespace(
+        params={"version": "v", "blind": expected, "fiducial_scale_cut": [1, 100]},
+        input={input_key: part_path},
+    )
+    error = ReachedSignal if expected == "none" else ValueError
+    with pytest.raises(error):
+        runpy.run_path(str(SCRIPTS / script), init_globals={"snakemake": smk})
+    assert not list(part_path.parent.glob("*.npz"))
+    assert not list(part_path.parent.glob("*.png"))
+
+
+@pytest.mark.parametrize("as_name", [False, True])
+@pytest.mark.parametrize("stamp", ["none", "toy"])
+def test_load_checks_expected_blind(part_path, as_name, stamp):
+    part = sacc_io.load(part_path)
+    part.metadata["blind"] = stamp
+    part.save_fits(str(part_path), overwrite=True)
+    expected = "toy" if as_name else blinding.Blind("toy")
+    if stamp == "none":
+        with pytest.raises(ValueError, match="stamped none.*toy"):
+            sacc_io.load(part_path, blind=expected)
+    else:
+        assert sacc_io.stamp(sacc_io.load(part_path, blind=expected)) == "toy"

src/sp_validation/tests/test_blinding.py

@@ -1,7 +1,7 @@
 """What a blind shifts, and what the catalogue config may declare.
 
-The shifts run the default theory, :func:`sp_validation.blinding.shear`, under
-a blind whose hidden point raises S8 above the fiducial.
+Small SACCs exercise the default theory and Smokescreen's parameter draw.
+The ``up`` fixture raises S8 above the fiducial for the E/B response tests.
 """
 
 import json
@@ -86,13 +86,143 @@ def test_pure_eb_under_a_blind_shifts_e_and_not_b(up):
     assert np.max(np.abs(response[2 * n : 4 * n])) < 0.1
 
 
+def _small_xi():
+    s = _sacc()
+    theta = np.array([10.0, 50.0])
+    sio.add_xi(s, (0, 0), theta, 0 * theta, 0 * theta, grid="reporting")
+    return s
+
+
+def test_distinct_blinds_and_determinism(tmp_path):
+    catalogues = {"paths": {"blinds": str(tmp_path / "blinds")}}
+    a, b = (bd.init(name, catalogues) for name in ("a", "b"))
+    s = _small_xi()
+    first = bd.conceal(s, a).mean
+    np.testing.assert_array_equal(first, bd.conceal(s, a).mean)
+    assert not np.array_equal(first, bd.conceal(s, b).mean)
+    np.testing.assert_array_equal(s.mean, np.zeros(len(s.mean)))
+
+
+def test_smokescreen_draw_stays_in_envelope(tmp_path):
+    from smokescreen.datavector import concealing_factor
+
+    blind = bd.init("envelope", {"paths": {"blinds": str(tmp_path)}})
+    record = blind.record()
+    keys = ("S8", "Omega_m")
+    fid = np.array([record["fiducial"][k] for k in keys])
+    envelope = np.array([record["envelope"][k] for k in keys])
+    for seed in [record["seed"], *range(1000)]:
+        # An identity theory observes the actual draw, without cosmology work.
+        delta = concealing_factor(
+            record["fiducial"],
+            record["envelope"],
+            seed=seed,
+            theory_fn=lambda p: np.array([p[k] for k in keys]),
+        )
+        hidden = fid + delta
+        assert np.all(hidden >= fid - envelope)
+        assert np.all(hidden <= fid + envelope)
+
+
+def test_multitracer_shift_matches_single_pair_theories(up, monkeypatch):
+    import pyccl as ccl
+    from smokescreen.datavector import draw_param_shifts
+
+    record = up.record()
+    shifts = draw_param_shifts(record["envelope"], record["seed"])
+    fid = record["fiducial"]
+    hidden = {k: value + shifts.get(k, 0.0) for k, value in fid.items()}
+    z = np.linspace(0.0, 3.0, 100)
+    s = sio.new_sacc(
+        {
+            i: (z, np.exp(-0.5 * ((z - centre) / 0.2) ** 2))
+            for i, centre in enumerate((0.5, 1.1))
+        }
+    )
+    pairs = ((0, 0), (0, 1), (1, 1))
+    for pair in pairs:
+        theta = np.array([10.0, 50.0])
+        sio.add_xi(s, pair, theta, 0 * theta, 0 * theta, grid="reporting")
+        ell = np.array([50.0, 100.0, 200.0])
+        sio.add_pseudo_cl(
+            s,
+            pair,
+            ell,
+            0 * ell,
+            0 * ell,
+            0 * ell,
+            window_ells=ell,
+            window_weights=np.eye(3),
+        )
+    # Observe CCL's actual inputs independently of blinding's pair selection.
+    make_tracer, angular_cl = ccl.WeakLensingTracer, ccl.angular_cl
+    nz_by_tracer, seen = {}, []
+
+    def tracer(*args, **kwargs):
+        result = make_tracer(*args, **kwargs)
+        nz_by_tracer[id(result)] = np.asarray(kwargs["dndz"][1]).copy()
+        return result
+
+    def spectrum(cosmo, a, b, ell, **kwargs):
+        seen.append((nz_by_tracer[id(a)], nz_by_tracer[id(b)]))
+        return angular_cl(cosmo, a, b, ell, **kwargs)
+
+    monkeypatch.setattr(ccl, "WeakLensingTracer", tracer)
+    monkeypatch.setattr(ccl, "angular_cl", spectrum)
+    shifted = bd.conceal(s, up)
+    for a, b in pairs:
+        nz_a, nz_b = s.tracers[f"source_{a}"].nz, s.tracers[f"source_{b}"].nz
+        # Each pair reaches angular_cl at fiducial and hidden points, for ξ and Cℓ.
+        assert (
+            sum(np.array_equal(x, nz_a) and np.array_equal(y, nz_b) for x, y in seen)
+            == 4
+        )
+    for a, b in pairs:
+        rows = s.indices(tracers=(f"source_{a}", f"source_{b}"))
+        single = s.copy()
+        single.keep_indices(rows)
+        expected = bd.shear(hidden, single) - bd.shear(fid, single)
+        np.testing.assert_allclose(
+            shifted.mean[rows] - s.mean[rows], expected, rtol=1e-12, atol=0
+        )
+        for kind in (sio.XI_PLUS, sio.XI_MINUS, sio.CL_EE):
+            assert np.any(expected[single.indices(kind)] != 0)
+
+
+def test_saved_blind_is_secret_free(tmp_path, capfd):
+    s = _small_xi()
+    path = tmp_path / "sealed.sacc"
+    blind = bd.init("secret-free", {"paths": {"blinds": str(tmp_path / "blinds")}})
+    params_seen = []
+
+    def theory(params, vector):
+        params_seen.append(dict(params))
+        return bd.shear(params, vector)
+
+    written = sio.save(s, path, blind=blind, theory=theory)
+    record = blind.record()
+    hidden = next(p for p in params_seen if p != record["fiducial"])
+    captured = capfd.readouterr()
+    secrets = [record["seed"]]
+    for key in ("S8", "Omega_m"):
+        value = hidden[key]
+        secrets.extend((repr(value), f"{value:.4f}", f"{value:.3f}"))
+    for secret in secrets:
+        assert secret.encode() not in path.read_bytes()
+        assert secret not in captured.out + captured.err
+    assert written.metadata == {**s.metadata, "blind": blind.name}
+    loaded = sio.load(path)
+    assert loaded.metadata == written.metadata
+    np.testing.assert_array_equal(loaded.mean, written.mean)
+
+
 def test_seal_refuses_a_sacc_already_stamped():
     sealed = sio.seal(_sacc(), bd.NONE)
     with pytest.raises(ValueError, match="already carries a blind stamp"):
         sio.seal(sealed, bd.NONE)
 
 
-def test_entries_reading_one_file_declare_one_blind():
+def test_entries_reading_one_file_declare_one_blind(tmp_path):
     catalogues = {
         "paths": {"blinds": "/blinds"},
         "Y3": {"blind": "y3", "subdir": "/d", "shear": {"path": "y3.fits"}},
@@ -102,3 +232,14 @@ def test_entries_reading_one_file_declare_one_blind():
     catalogues["TWIN"]["shear"]["path"] = "/d/sub/../y3.fits"
     with pytest.raises(ValueError, match="different blinds"):
         bd.blind_of(catalogues, "TWIN")
+
+    physical = tmp_path / "y3.fits"
+    physical.touch()
+    alias = tmp_path / "alias.fits"
+    alias.symlink_to(physical)
+    catalogues["Y3"]["shear"]["path"] = str(physical)
+    catalogues["TWIN"]["shear"]["path"] = str(alias)
+    with pytest.raises(ValueError, match="different blinds"):
+        bd.blind_of(catalogues, "TWIN")
+    catalogues["TWIN"]["blind"] = "y3"
+    assert bd.blind_of(catalogues, "TWIN") == "y3"

src/sp_validation/tests/test_cosmo_val.py

@@ -649,6 +649,39 @@ def test_a_blinded_catalogue_refuses_tomographic_signal(self, tmp_path):
             cv.calculate_2pcf_version(version, compute_tomography=True)
         with pytest.raises(blinding.BlindingError, match="tomographic"):
             cv.calculate_pseudo_cl(compute_tomography=True)
+        # Direct producers must refuse before opening even a missing catalogue.
+        cv.cc[version]["shear"]["path"] = str(tmp_path / "missing.fits")
+        out = tmp_path / "tomographic.fits"
+        with pytest.raises(blinding.BlindingError, match="tomographic"):
+            cv.calculate_pseudo_cl_map(version, 16, out, 1, 2)
+        with pytest.raises(blinding.BlindingError, match="tomographic"):
+            cv.calculate_pseudo_cl_catalog(version, out, 1, 2)
+        assert not out.exists()
+
+    @pytest.mark.parametrize("stamp", ["none", "toy"])
+    def test_held_xi_matches_declared_blind(self, tmp_path, stamp):
+        from sp_validation import blinding, sacc_io
+
+        params, version = self._write_synthetic_catalogs(tmp_path)
+        config = yaml.safe_load(open(params["catalog_config"]))
+        config["paths"]["blinds"] = str(tmp_path / "blinds")
+        config[version]["blind"] = "toy"
+        with open(params["catalog_config"], "w") as f:
+            yaml.dump(config, f)
+        blinding.init("toy", config)
+        cv = CosmologyValidation(versions=[version], npatch=1, **params)
+        part = sacc_io.new_sacc({0: ([0.1, 1.0], [1.0, 1.0])})
+        sacc_io.add_xi(part, (0, 0), [10.0], [1.0], [0.5], grid="reporting")
+        part.metadata["blind"] = stamp
+        cv.xi_parts[version, "reporting"] = part
+        if stamp == "none":
+            with pytest.raises(ValueError, match="stamped none.*toy"):
+                cv.calculate_2pcf_version(version)
+        else:
+            result = cv.calculate_2pcf_version(version)
+            np.testing.assert_array_equal(
+                result["tomo_bin_all_tomo_bin_all"].xip, [1.0]
+            )
 
     def test_treecorr_runs_on_the_cpus_the_process_holds(self, tmp_path):
         """By default TreeCorr takes the process's CPU affinity, not the node's count."""