diff --git a/src/sp_validation/blinding.py b/src/sp_validation/blinding.py index 09c49718..619b060a 100644 --- a/src/sp_validation/blinding.py +++ b/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): diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index a728d3eb..361369ac 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/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)) diff --git a/src/sp_validation/cosmo_val/real_space.py b/src/sp_validation/cosmo_val/real_space.py index 8678b828..a486e8b9 100644 --- a/src/sp_validation/cosmo_val/real_space.py +++ b/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) diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index 2c2e7919..c4bdc5c7 100644 --- a/src/sp_validation/sacc_io.py +++ b/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 # ============================================================================= diff --git a/src/sp_validation/tests/test_blind_consumers.py b/src/sp_validation/tests/test_blind_consumers.py new file mode 100644 index 00000000..3d3b7a8d --- /dev/null +++ b/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" diff --git a/src/sp_validation/tests/test_blinding.py b/src/sp_validation/tests/test_blinding.py index d148e88b..79f5647c 100644 --- a/src/sp_validation/tests/test_blinding.py +++ b/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" diff --git a/src/sp_validation/tests/test_cosmo_val.py b/src/sp_validation/tests/test_cosmo_val.py index 25cf0865..9d4ca719 100644 --- a/src/sp_validation/tests/test_cosmo_val.py +++ b/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.""" diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk index 81d305d1..08526643 100644 --- a/workflow/rules/cosmo_val.smk +++ b/workflow/rules/cosmo_val.smk @@ -361,6 +361,7 @@ rule cv_plot_pseudo_cl: **cv_pseudo_cl_figures(), params: versions=CV_VERSIONS, + blinds=[blind_of(v) for v in CV_VERSIONS], # Style is per catalogue, so the derived variants take their parent's. markers=[CATALOG_CONFIG[base_version(v)]["marker"] for v in CV_VERSIONS], colours=[CATALOG_CONFIG[base_version(v)]["colour"] for v in CV_VERSIONS], @@ -390,6 +391,7 @@ rule cv_pure_eb: **cv_pure_eb_figures("{version}"), params: version="{version}", + blind=lambda w: blind_of(w.version), min_sep=CV["theta_min"], max_sep=CV["theta_max"], nbins=CV["nbins"], @@ -417,6 +419,7 @@ rule cv_cosebis: **cv_cosebis_figures("{version}"), params: version="{version}", + blind=lambda w: blind_of(w.version), min_sep=XI_GRIDS["integration"]["min_sep"], max_sep=XI_GRIDS["integration"]["max_sep"], nbins=XI_GRIDS["integration"]["nbins"], @@ -448,6 +451,7 @@ rule cv_summarize_bmodes: summary_json=str(COSMO_VAL / "bmode_summary.json"), params: versions=CV_VERSIONS, + blinds=[blind_of(v) for v in CV_VERSIONS], fiducial_scale_cut=CV["fiducial_scale_cut"], include_pseudo_cl=CV.get("include_pseudo_cl", False), resources: diff --git a/workflow/scripts/assemble_sacc.py b/workflow/scripts/assemble_sacc.py index d756bce2..034720aa 100644 --- a/workflow/scripts/assemble_sacc.py +++ b/workflow/scripts/assemble_sacc.py @@ -136,7 +136,8 @@ def assemble_sacc( path = part_paths.get(name) if path is None: continue - parts.append(_attach_cov(sacc_io.load(path), name, xi_cov, pseudo_cl_cov)) + part = sacc_io.load(path, blind=blind) + parts.append(_attach_cov(part, name, xi_cov, pseudo_cl_cov)) if not parts: raise ValueError(f"no parts found for {version}: {part_paths}") s = sacc_io.save( diff --git a/workflow/scripts/cv_cosebis.py b/workflow/scripts/cv_cosebis.py index ee9d8aef..c7176f82 100644 --- a/workflow/scripts/cv_cosebis.py +++ b/workflow/scripts/cv_cosebis.py @@ -27,7 +27,7 @@ version = p["version"] fiducial_scale_cut = tuple(p["fiducial_scale_cut"]) -part = sacc_io.load(snakemake.input["xi"]) +part = sacc_io.load(snakemake.input["xi"], blind=p["blind"]) theta, xip, xim = sacc_io.get_xi(part, (0, 0), grid="integration") edges = log_bin_edges(p["min_sep"], p["max_sep"], p["nbins"]) diff --git a/workflow/scripts/cv_plot_pseudo_cl.py b/workflow/scripts/cv_plot_pseudo_cl.py index f9aeeda9..dc1b32a1 100644 --- a/workflow/scripts/cv_plot_pseudo_cl.py +++ b/workflow/scripts/cv_plot_pseudo_cl.py @@ -18,7 +18,7 @@ spectra = {} for i, version in enumerate(p["versions"]): - part = sacc_io.load(snakemake.input["pseudo_cl"][i]) + part = sacc_io.load(snakemake.input["pseudo_cl"][i], blind=p["blinds"][i]) ell, ee, bb, eb, _window = sacc_io.get_pseudo_cl(part, (0, 0)) with fits.open(snakemake.input["pseudo_cl_cov"][i]) as hdul: covs = { diff --git a/workflow/scripts/cv_pure_eb.py b/workflow/scripts/cv_pure_eb.py index 3b582699..25ee3f72 100644 --- a/workflow/scripts/cv_pure_eb.py +++ b/workflow/scripts/cv_pure_eb.py @@ -29,7 +29,7 @@ version = p["version"] fiducial_scale_cut = tuple(p["fiducial_scale_cut"]) -part = sacc_io.load(snakemake.input["xi_integration"]) +part = sacc_io.load(snakemake.input["xi_integration"], blind=p["blind"]) gg_int = sacc_io.xi_correlation(part, grid="integration") # A part stores bin centres; the edges come from the grid it was measured on. grid = p["integration"] diff --git a/workflow/scripts/cv_runner.py b/workflow/scripts/cv_runner.py index bb67cdfb..b4ffcea2 100644 --- a/workflow/scripts/cv_runner.py +++ b/workflow/scripts/cv_runner.py @@ -41,7 +41,8 @@ def take_reporting_parts(cv, paths): from sp_validation import sacc_io for version, path in zip(cv.versions, paths, strict=True): - cv.xi_parts[version, "reporting"] = sacc_io.load(path) + part = sacc_io.load(path, blind=cv.blind(version)) + cv.xi_parts[version, "reporting"] = part def verify_outputs(snakemake): diff --git a/workflow/scripts/cv_summarize_bmodes.py b/workflow/scripts/cv_summarize_bmodes.py index 3f9d9e82..802de42c 100644 --- a/workflow/scripts/cv_summarize_bmodes.py +++ b/workflow/scripts/cv_summarize_bmodes.py @@ -48,7 +48,7 @@ if p["include_pseudo_cl"]: from astropy.io import fits - part = sacc_io.load(snakemake.input["pseudo_cl"][i]) + part = sacc_io.load(snakemake.input["pseudo_cl"][i], blind=p["blinds"][i]) _ell, _ee, bb, _eb, _window = sacc_io.get_pseudo_cl(part, (0, 0)) with fits.open(snakemake.input["pseudo_cl_cov"][i]) as hdul: cov_bb = np.asarray(hdul["COVAR_BB_BB"].data, float)