From 9d02f52ea7508317a314b8ee1a4a11a9bbfc357e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 23:46:20 +0200 Subject: [PATCH 1/7] Refactor blinding onto a theory_fn protocol with a default CCL backend Replace the firecrown-likelihood theory path with a single theory_fn(cosmo_params) -> np.ndarray protocol. ConcealDataVector now takes (fiducial_params, shifts_dict, sacc_data, *, seed, theory_fn=None, shift_distr): both theory vectors come from theory_fn (the shipped default CCL cosmic-shear backend when None, else the caller's), the concealed point is the fiducial mapping with drawn deltas overlaid, and the save path writes the blinded mean back through SACC's own API with the covariance carried over unchanged. - smokescreen/backends/ccl.py: default cosmic-shear theory_fn built from a SACC's n(z) tracers and cosmic-shear rows (Cl_ee, xi_plus/minus). pyccl is imported here, so `import smokescreen` imports no theory backend. - param_shifts.draw_param_shifts: pure, order-independent (sorted keys, one scalar draw per key) local-RNG draw returning deltas; seed int|str with SHA-256 string normalization; no global np.random.seed. - Drop _verify_sacc_consistency, _create_concealed_cosmo, systm_dict, save_to_sacc/get_sacc_indices, and the firecrown/ParamsMap imports. The sole construction guard is a theory-vs-mean length check. - utils: pyccl/sacc imports made lazy; remove firecrown-only helper. - __main__: CLI reworked onto the new signature + default backend. The firecrown cosmic-shear example likelihood is retained in examples/ as inherited-upstream, unsupported code; the package imports it nowhere. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/smokescreen/__main__.py | 95 ++--- src/smokescreen/backends/__init__.py | 10 + src/smokescreen/backends/ccl.py | 132 ++++++ src/smokescreen/datavector.py | 574 ++++++--------------------- src/smokescreen/param_shifts.py | 188 ++++----- src/smokescreen/utils.py | 35 +- 6 files changed, 372 insertions(+), 662 deletions(-) create mode 100644 src/smokescreen/backends/__init__.py create mode 100644 src/smokescreen/backends/ccl.py diff --git a/src/smokescreen/__main__.py b/src/smokescreen/__main__.py index 5810275..39d4102 100644 --- a/src/smokescreen/__main__.py +++ b/src/smokescreen/__main__.py @@ -3,103 +3,84 @@ from typing import Union, Dict, Tuple from jsonargparse import CLI from jsonargparse.typing import Path_drw, Path_fr -from pyccl import Cosmology as CosmologyType -import pyccl as ccl # warnings related to sacc files import warnings from smokescreen import ConcealDataVector from smokescreen.encryption import encrypt_file, decrypt_file -from smokescreen.utils import load_cosmology_from_partial_dict, load_sacc_file +from smokescreen.utils import load_sacc_file from . import __version__ warnings.filterwarnings("ignore") # banner to be printed in the terminal banner = rf""" - ( - )\ ) ) -(()/( ) ( /( ( ( ( ( - /(_)) ( ( )\()) ))\ ( ( )( ))\ ))\ ( -(_)) )\ ' )\ ((_)\ /((_))\ )\(()\ /((_)/((_) )\ ) -/ __| _((_)) ((_)| |(_)(_)) ((_) ((_)((_)(_)) (_)) _(_/( -\__ \| ' \()/ _ \| / / / -_)(_- None: - r"""Main function to conceal a SACC file using a firecrown likelihood. + r"""Conceal a cosmic-shear SACC file with the default CCL theory backend. Args: - path_to_sacc (str): Path to the sacc file to blind. - likelihood_path (str): Path to the firecrown likelihood module file. - shifts_dict (dict): Dictionary with fixed values for the firecrown shifts parameters. - Example: {"Omega_c": (0.20, 0.39), "sigma8": (0.70, 0.90)} - shift_type (str): Type of shift to apply to the data vector. - Options are 'add' and 'mult'. Defaults to 'add'. - systematics (dict): Dictionary with fixed values for the firecrown systematics parameters. - shift_distribution (str): Distribution type for the parameter shifts. - Options are 'flat' and 'gaussian'. Defaults to 'flat'. - seed (int, str): Seed for the blinding process. Defaults to 2112. - reference_cosmology (Union[CosmologyType, dict]): - Cosmology object or dictionary with cosmological - parameters you want different than the VanillaLCDM as reference cosmology. - Defaults to ccl.CosmologyVanillaLCDM(). - path_to_output (str): Path to save the blinded sacc file. Defaults to None. - keep_original_sacc (bool): If True, keeps the original sacc file. - Defaults to False [keeps only the encrypted file]. + path_to_sacc (str): Path to the SACC file to blind. It must contain + exactly the cosmic-shear rows the default CCL backend models + (galaxy_shear_cl_ee and/or galaxy_shear_xi_plus/minus), with + weak-lensing tracers carrying n(z). + fiducial_params (dict): Fiducial cosmological parameters (CCL-native + names, e.g. {"sigma8": 0.8, "Omega_c": 0.25, "Omega_b": 0.05, + "h": 0.67, "n_s": 0.96}). + shifts_dict (dict): Shift envelopes, interpreted as deltas about zero. + Example: {"Omega_c": (-0.05, 0.05), "sigma8": 0.05} + seed (int, str): Seed for the blinding process (no default; must be a + deliberate, secret choice). + shift_type (str): Concealing factor type, 'add' or 'mult'. Default 'add'. + shift_distribution (str): 'flat' or 'gaussian'. Default 'flat'. + path_to_output (str): Directory to save the blinded SACC. Default None + (uses the input file's directory). + keep_original_sacc (bool): If True, keeps the original SACC file. + Default False (keeps only the encrypted file). output_suffix (str): Custom suffix for the output file name. - Defaults to None (uses 'concealed_data_vector'). """ print(banner) - if isinstance(reference_cosmology, dict): - cosmo = load_cosmology_from_partial_dict(reference_cosmology) - else: - cosmo = reference_cosmology - # tests if the sacc file exists assert os.path.exists(path_to_sacc), f"File {path_to_sacc} does not exist." - assert os.path.exists(likelihood_path), f"File {likelihood_path} does not exist." # reads the sacc file (returns sacc object and detected format) sacc_data, input_format = load_sacc_file(path_to_sacc) - # creates the smokescreen object - smoke = ConcealDataVector(cosmo, likelihood_path, shifts_dict, sacc_data, systematics, seed, - shift_distr=shift_distribution, input_format=input_format) + # creates the smokescreen object with the default CCL backend + smoke = ConcealDataVector(fiducial_params, shifts_dict, sacc_data, + seed=seed, shift_distr=shift_distribution, + input_format=input_format) # blinds the sacc file smoke.calculate_concealing_factor(factor_type=shift_type) - # applies the blinding factor to the sacc file smoke.apply_concealing_to_likelihood_datavec() print(f">> User {getpass.getuser()}", f"used Smokescreen on {path_to_sacc} ... it is super effective!") # get root name of the input file root_name = os.path.splitext(os.path.basename(path_to_sacc))[0] # saves the blinded sacc file - if path_to_output is not None: - smoke.save_concealed_datavector(path_to_output, root_name, - output_format=input_format, - suffix=output_suffix) - else: - # get the input file directory + if path_to_output is None: path_to_output = os.path.dirname(path_to_sacc) - smoke.save_concealed_datavector(path_to_output, root_name, - output_format=input_format, - suffix=output_suffix) - # Determine extension based on format + smoke.save_concealed_datavector(path_to_output, root_name, + output_format=input_format, + suffix=output_suffix) ext = '.hdf5' if input_format == 'hdf5' else '.fits' _suffix = output_suffix if output_suffix is not None else "concealed_data_vector" outprintfile = f"{path_to_output}/{root_name}_{_suffix}{ext}" diff --git a/src/smokescreen/backends/__init__.py b/src/smokescreen/backends/__init__.py new file mode 100644 index 0000000..67d6bee --- /dev/null +++ b/src/smokescreen/backends/__init__.py @@ -0,0 +1,10 @@ +# author: Arthur Loureiro +# license: BSD 3-Clause +''' +Theory backends (:mod:`smokescreen.backends`) +============================================= + +Built-in ``theory_fn`` implementations. The default CCL cosmic-shear backend +lives in :mod:`smokescreen.backends.ccl`; it is imported lazily (importing this +package does not import ``pyccl``). +''' diff --git a/src/smokescreen/backends/ccl.py b/src/smokescreen/backends/ccl.py new file mode 100644 index 0000000..079a068 --- /dev/null +++ b/src/smokescreen/backends/ccl.py @@ -0,0 +1,132 @@ +# author: Arthur Loureiro +# license: BSD 3-Clause +''' +Default CCL cosmic-shear backend (:mod:`smokescreen.backends.ccl`) +================================================================= + +.. currentmodule:: smokescreen.backends.ccl + +Builds a default ``theory_fn`` that computes a cosmic-shear theory vector from +CCL for a standard weak-lensing SACC file. ``pyccl`` is imported here, so it +enters ``sys.modules`` only once this backend is constructed --- importing +``smokescreen`` does not import it. + +The backend accepts CCL-native cosmological-parameter names +(``sigma8``, ``Omega_c``, ``Omega_b``, ``h``, ``n_s``, ...). It reads the +weak-lensing tracers' n(z) and the SACC row layout, and returns, for a given +parameter mapping, the theory vector aligned element-for-element to +``sacc_data.mean``. + +.. autofunction:: build_ccl_theory_fn +''' +import numpy as np +import pyccl as ccl + + +# SACC data types the default backend understands, and the CCL correlation +# 'type' string (Fourier -> real space) each ξ± measurement maps to. +_XI_TYPES = { + "galaxy_shear_xi_plus": "GG+", + "galaxy_shear_xi_minus": "GG-", +} +_CL_TYPES = {"galaxy_shear_cl_ee"} + + +def _tracer_nz(sacc_data, name): + """Return (z, nz) for a SACC tracer carrying a redshift distribution.""" + tracer = sacc_data.tracers[name] + if not (hasattr(tracer, "z") and hasattr(tracer, "nz")): + raise ValueError( + f"Default CCL backend requires n(z) tracers; tracer '{name}' " + f"({tracer.tracer_type}) carries no redshift distribution." + ) + return np.asarray(tracer.z), np.asarray(tracer.nz) + + +def _plan_blocks(sacc_data): + """ + Decompose the SACC into cosmic-shear blocks aligned to ``sacc_data.mean``. + + Returns a list of ``(data_type, tracer_pair, angle_values, row_indices)``, + one per (data_type, tracer-pair) group the default backend supports. Raises + if the SACC carries a data type this backend does not model. + """ + supported = set(_XI_TYPES) | _CL_TYPES + present = set(sacc_data.get_data_types()) + unsupported = present - supported + if unsupported: + raise ValueError( + f"Default CCL backend only models cosmic-shear data types " + f"{sorted(supported)}; SACC carries unsupported {sorted(unsupported)}. " + f"Supply an explicit theory_fn for this SACC." + ) + + blocks = [] + for data_type in present: + angle_key = "theta" if data_type in _XI_TYPES else "ell" + for pair in sacc_data.get_tracer_combinations(data_type): + idx = sacc_data.indices(data_type, pair) + angles = np.array( + [sacc_data.data[i].get_tag(angle_key) for i in idx] + ) + blocks.append((data_type, pair, angles, np.asarray(idx))) + return blocks + + +def build_ccl_theory_fn(sacc_data): + """ + Build a default cosmic-shear ``theory_fn`` from a SACC file. + + Parameters + ---------- + sacc_data : sacc.sacc.Sacc + SACC carrying weak-lensing tracers with n(z) and cosmic-shear rows + (``galaxy_shear_cl_ee`` and/or ``galaxy_shear_xi_plus`` / + ``galaxy_shear_xi_minus``). The row layout it exposes is the layout the + returned callable aligns to. + + Returns + ------- + Callable[[Mapping[str, float]], np.ndarray] + A ``theory_fn``: given a mapping of CCL-native cosmological-parameter + names to values, it constructs a :class:`pyccl.Cosmology`, builds a + :class:`pyccl.WeakLensingTracer` per bin, and returns the cosmic-shear + theory vector aligned to ``sacc_data.mean``. + """ + blocks = _plan_blocks(sacc_data) + n_rows = len(sacc_data.mean) + nz_cache = { + name: _tracer_nz(sacc_data, name) + for block in blocks + for name in block[1] + } + # ℓ grid for the Cℓ used by the real-space correlation transform. + ell_grid = np.unique(np.geomspace(2, 30000, 512).astype(int)) + + def theory_fn(cosmo_params): + # Default to pyccl-native transfer/power (no CAMB/CLASS dependency), so + # the shipped backend runs against a bare pyccl install. A caller-supplied + # theory_fn is free to route through a Boltzmann code. + params = {"transfer_function": "eisenstein_hu", + "matter_power_spectrum": "halofit"} + params.update(cosmo_params) + cosmo = ccl.Cosmology(**params) + wl = { + name: ccl.WeakLensingTracer(cosmo, dndz=nz_cache[name]) + for name in nz_cache + } + vec = np.empty(n_rows) + for data_type, (t1, t2), angles, idx in blocks: + cl = ccl.angular_cl(cosmo, wl[t1], wl[t2], ell_grid) + if data_type in _CL_TYPES: + values = np.interp(angles, ell_grid, cl) + else: + theta_deg = angles / 60.0 # arcmin -> degrees + values = ccl.correlation( + cosmo, ell=ell_grid, C_ell=cl, theta=theta_deg, + type=_XI_TYPES[data_type], + ) + vec[idx] = values + return vec + + return theory_fn diff --git a/src/smokescreen/datavector.py b/src/smokescreen/datavector.py index ec8e97a..ae71d5c 100644 --- a/src/smokescreen/datavector.py +++ b/src/smokescreen/datavector.py @@ -7,7 +7,14 @@ .. currentmodule:: smokescreen.datavector The :mod:`smokescreen.datavector` module provides functionalities -to conceal data vectors in the context of cosmological analysis. +to conceal (blind) data vectors in the context of cosmological analysis. + +Blinding adds a theory difference to the data vector, +``d -> d + t(hidden) - t(fiducial)``. The theory source is a single +``theory_fn(cosmo_params) -> np.ndarray`` protocol: any callable returning a +vector aligned to the SACC rows being concealed is a valid backend. The fork +ships a default CCL cosmic-shear backend (see +:mod:`smokescreen.backends.ccl`); power users supply their own callable. Conceal Data Vector Class @@ -16,491 +23,142 @@ .. autoclass:: ConcealDataVector :members: :undoc-members: - :inherited-members: ''' -import os -import types -import inspect import datetime import getpass from copy import deepcopy -from packaging.version import Version + import numpy as np -import pyccl as ccl import sacc -import firecrown - -# Handle different Firecrown versions -if Version(firecrown.__version__) >= Version("1.15.0a0"): - # New structure (1.15.0a0+) - from firecrown.likelihood import ( - load_likelihood, - load_likelihood_from_module_type, - NamedParameters - ) -else: - # Old structure (< 1.15.0a0) - from firecrown.likelihood.likelihood import ( - load_likelihood, - load_likelihood_from_module_type, - NamedParameters - ) - -from firecrown.parameters import ParamsMap -from firecrown.updatable import get_default_params_map -from firecrown.utils import save_to_sacc -from firecrown.ccl_factory import PoweSpecAmplitudeParameter - - -from smokescreen.param_shifts import draw_flat_or_deterministic_param_shifts -from smokescreen.param_shifts import draw_gaussian_param_shifts -from smokescreen.utils import load_module_from_path, modify_default_params + +from smokescreen.param_shifts import draw_param_shifts class ConcealDataVector(): """ - Class for calling a smokescreen on the measured data-vector. + Conceal (blind) a measured data vector by adding a theory difference. - FIXME: Only cosmological parameters are supported for now for the shifts + Both theory vectors (fiducial and hidden) come from a single ``theory_fn``: + the shipped default CCL backend when none is supplied, else the caller's. + There is no likelihood, no ``pyccl.Cosmology`` object on this path, and no + systematics dictionary transiting the class --- systematics, if a backend + needs them, are closed over inside that backend's ``theory_fn``. Parameters ---------- - cosmo : pyccl.Cosmology - Cosmology object from CCL with a fiducial cosmology. - likelihood : str or module - path to the likelihood or a module containing the likelihood - must contain both `build_likelihood` and `compute_theory_vector` methods - shifts_dict : dict - Dictionary of parameter names and corresponding shift widths. If the - shifts are single values, the dictionary values should be the shift - widths. If the shifts are tuples of values, the dictionary values - should be the (lower, upper) bounds of the shift widths. + fiducial_params : Mapping[str, float] + The fiducial point, a plain mapping of cosmological-parameter name to + value over a free-form parameter space (CCL-native names, e.g. + ``sigma8``, ``Omega_c``). Interpreting the names is the ``theory_fn``'s + job; this is not a ``pyccl.Cosmology``. + shifts_dict : Mapping[str, float | tuple[float, float]] + Maps a parameter name (a key of ``fiducial_params``) to its shift + envelope, interpreted as a *delta* about zero. A ``float`` ``h`` is the + symmetric delta envelope ``(-h, +h)``; a ``(lo, hi)`` tuple is a delta + box. The drawn delta is added to ``fiducial_params[k]``. sacc_data : sacc.sacc.Sacc - Data-vector to be concealed (blinded). - If None, the data-vector will be loaded from the likelihood. - systm_dict : dict - Dictionary of systematics names and corresponding fiducial values. - Default is None, which means that the systematics will be loaded - from firecrown defaults by the likelihood. + Data vector to be concealed. It must contain exactly the rows the + ``theory_fn`` returns, in ``sacc_data.mean`` order (extract-then-blind: + the caller scopes to the to-be-blinded block before construction). seed : int or str - Random seed. + Random seed. No default: blinding custody depends on the seed being a + deliberate, secret choice. A ``str`` is normalized to an ``int``. + theory_fn : Callable[[Mapping[str, float]], np.ndarray], optional + Theory backend. When ``None``, the default CCL cosmic-shear backend is + built from ``sacc_data``. + shift_distr : str + ``"flat"`` (uniform over the delta envelope) or ``"gaussian"`` + (zero-mean Gaussian with sigma = the ``float`` half-width). Default + ``"flat"``. Keyword Arguments ----------------- - shift_type : str - Type of shift to be applied. Default is "flat". + input_format : str + Original SACC file format for output preservation ('fits' or 'hdf5'). debug : bool If True, prints debug information. Default is False. - - """ - def __init__(self, cosmo, likelihood, shifts_dict, sacc_data, systm_dict=None, - seed="2112", **kwargs): - """ - unit - """ - # save the cosmology - self.cosmo = cosmo - # save the systematics dictionary - self.systematics_dict = systm_dict - # save the data-vector + def __init__(self, fiducial_params, shifts_dict, sacc_data, *, + seed, theory_fn=None, shift_distr="flat", **kwargs): + assert isinstance(sacc_data, sacc.sacc.Sacc), "sacc_data must be a sacc object" + + self.fiducial_params = dict(fiducial_params) + self.shifts_dict = shifts_dict self.sacc_data = sacc_data - # seed for the random number generator self.seed = seed - # checks if the sacc_data is in the correct format: - assert isinstance(self.sacc_data, sacc.sacc.Sacc), "sacc_data must be a sacc object" - # save the shifts - self.shifts_dict = shifts_dict + self.shift_distr = shift_distr + self._debug = bool(kwargs.get('debug', False)) # detect original file format for output preservation - # input_format can be passed via kwargs or detected from path_to_sacc - if 'input_format' in kwargs: - self._input_format = kwargs['input_format'] + self._input_format = kwargs.get( + 'input_format', + getattr(sacc_data, '_smokescreen_input_format', 'fits'), + ) + + # theory backend: caller-supplied, or the default CCL cosmic-shear one + # (imported lazily so `import smokescreen` never imports pyccl). + if theory_fn is None: + from smokescreen.backends.ccl import build_ccl_theory_fn + self.theory_fn = build_ccl_theory_fn(sacc_data) else: - # For backwards compatibility - try to get the sacc_data's source info - self._input_format = getattr(sacc_data, '_smokescreen_input_format', 'fits') + self.theory_fn = theory_fn - # load the shifts - # Check for 'shift_type' keyword argument - if 'shift_distr' in kwargs: - self.__shifts = self._load_shifts(seed, shift_distr=kwargs['shift_distr']) - else: - self.__shifts = self._load_shifts(seed) - - # create concealed cosmology object: - self.__concealed_cosmo = self._create_concealed_cosmo(self.__shifts) + # draw the hidden shift and overlay it on the fiducial point + self.__shifts = draw_param_shifts(self.shifts_dict, seed, + shift_distr=shift_distr) + self.__concealed_params = self._overlay_shifts(self.__shifts) - if 'debug' in kwargs and kwargs['debug']: - self._debug = True - print(f"[DEBUG] Shifts: {self.__shifts}") - print(f"[DEBUG] Concealed Cosmology: {self.__concealed_cosmo}") - else: - self._debug = False - - # load the likelihood - self.likelihood, self.tools = self._load_likelihood(likelihood, - self.sacc_data) - # # create the smokescreen data-vector - # self.smokescreen_data = self.create_smokescreen_data() - - # load the systematics - if self.systematics_dict is None: - self.systematics = self._load_default_systematics(self.likelihood) - else: - self.systematics = self._load_systematics(self.systematics_dict, self.likelihood) - - def _load_likelihood(self, likelihood, sacc_data): - """ - Loads the likelihood either from a python module or from a file. - - Parameters - ---------- - likelihood : str or module - path to the likelihood or a module containing the likelihood - must contain both `build_likelihood` and `compute_theory_vector` methods - """ - - build_parameters = NamedParameters({'sacc_data': sacc_data}) - - if type(likelihood) is str: - # check if the file can be found - if not os.path.isfile(likelihood): - raise FileNotFoundError(f'Could not find file {likelihood}') - - # test the likelihood - self._test_likelihood(likelihood, 'str') - # load the likelihood from the file - likelihood, tools = load_likelihood(likelihood, build_parameters) - # because now firecrown needs to know the amplitude parameter - # before we build the likelihood, need to check if we are - # concealing the correct parameter - self._check_amplitude_parameter(tools) - - # check if the likelihood has a compute_vector method - if not hasattr(likelihood, 'compute_theory_vector'): # pragma: no cover - raise AttributeError('Likelihood does not have a compute_vector method') - - # Verify SACC consistency after loading likelihood - self._verify_sacc_consistency(likelihood) - - return likelihood, tools - - elif isinstance(likelihood, types.ModuleType): - # test the likelihood - self._test_likelihood(likelihood, 'module') - - # tries to load the likelihood from the module - likelihood, tools = load_likelihood_from_module_type(likelihood, - build_parameters) - # because now firecrown needs to know the amplitude parameter - # before we build the likelihood, need to check if we are - # concealing the correct parameter - self._check_amplitude_parameter(tools) - # check if the likelihood has a compute_vector method - if not hasattr(likelihood, 'compute_theory_vector'): # pragma: no cover - raise AttributeError('Likelihood does not have a compute_vector method') - - # Verify SACC consistency after loading likelihood - self._verify_sacc_consistency(likelihood) - return likelihood, tools - else: - raise TypeError('Likelihood must be a string path to a likelihood module or a module') - - def _verify_sacc_consistency(self, likelihood): - """ - Verifies that the user-provided SACC data vector and covariance match - what the likelihood internally uses. - - After loading the likelihood, this method compares: - - self.sacc_data.mean (user's data vector) vs self.likelihood.get_data_vector() - - self.sacc_data.covariance (user's covariance) vs self.likelihood.get_cov() - - Raises - ------ - ValueError - If the data vector or covariance matrix don't match between - the user-provided SACC file and the likelihood's internal values. - """ - # Get the internal data vector and covariance from the likelihood - internal_data_vector = likelihood.get_data_vector() - internal_covariance = likelihood.get_cov() - - # Get the user-provided SACC data - user_data_vector = self.sacc_data.mean - - # Handle covariance - it could be None or a dense matrix - if self.sacc_data.covariance is not None: - user_covariance = self.sacc_data.covariance.dense - else: - user_covariance = None - - # Check data vector consistency - if not np.allclose(user_data_vector, internal_data_vector, rtol=1e-10, atol=1e-10): - # Calculate sum of absolute differences for reporting - data_diff_sum = np.sum(np.abs(user_data_vector - internal_data_vector)) - - raise ValueError( - f"Data vector mismatch between user-provided SACC and likelihood. " - f"Expected shape {internal_data_vector.shape}, got {user_data_vector.shape}. " - f"Sum of absolute differences: {data_diff_sum:.6e}" - ) - - # Check covariance consistency - if user_covariance is not None and internal_covariance is not None: - if not np.allclose(user_covariance, internal_covariance, rtol=1e-10, atol=1e-10): - # Calculate norm of difference for reporting - cov_diff_norm = np.linalg.norm(user_covariance - internal_covariance) - - raise ValueError( - f"Covariance matrix mismatch between user-provided SACC and likelihood. " - f"Expected shape {internal_covariance.shape}, got {user_covariance.shape}. " - f"Norm of difference: {cov_diff_norm:.6e}" - ) - elif user_covariance is not None and internal_covariance is None: - raise ValueError( - "User-provided SACC has covariance but likelihood returns None for covariance." - ) - elif user_covariance is None and internal_covariance is not None: + # length guard: theory vector must align to the SACC rows + self.theory_vec_fid = np.asarray(self.theory_fn(self.fiducial_params)) + if len(self.theory_vec_fid) != len(self.sacc_data.mean): raise ValueError( - "Likelihood has covariance but user-provided SACC has None for covariance." + f"theory_fn returned {len(self.theory_vec_fid)} elements but " + f"sacc_data.mean has {len(self.sacc_data.mean)}; the SACC must " + f"contain exactly the rows theory_fn returns." ) - def _test_likelihood(self, likelihood, like_type): - """ - Tests if the likelihood has the required methods. - - Parameters - ---------- - likelihood : str or module - path to the likelihood or a module containing the likelihood - must contain both `build_likelihood` and `compute_theory_vector` methods - like_type : str - Type of likelihood. Can be either 'str' or 'module'. - """ - if like_type == "str": - likelihood = load_module_from_path(likelihood) - else: - likelihood = likelihood - - # check if the module has a build_likelihood method - if not hasattr(likelihood, 'build_likelihood'): - raise AttributeError('Likelihood does not have a build_likelihood method') - - if self.sacc_data is not None: - sig = inspect.signature(likelihood.build_likelihood) - likefunc_params = sig.parameters - assert len(likefunc_params) >= 1, ("A sacc was provided, ", - "the likelihood must require a", - "build_parameters NamedParameters object!") - - def _check_amplitude_parameter(self, tools): - """ - Checks if the amplitude parameter is set in the tools is the same - as the one in the cosmology and in the concealing dictionary. - If not, raises an error. - - Parameters - ---------- - tools : firecrown.ccl_factory.CCLFactory - CCLFactory object with the cosmology and the amplitude parameter. - - Raises - ------ - ValueError - If the amplitude parameter is not supported or if the required parameter - is not in the cosmology or in the shifts dictionary. - """ - _amplitude_param = tools.ccl_factory.amplitude_parameter - - if _amplitude_param is PoweSpecAmplitudeParameter.SIGMA8: - _required_param = 'sigma8' - elif _amplitude_param is PoweSpecAmplitudeParameter.AS: - _required_param = 'A_s' - else: - raise ValueError(f"Amplitude parameter {_amplitude_param} not supported") - - _error_msg = "\n You probably need to set the amplitude parameter [A_s/sigma8] " - _error_msg += "that you want to conceal when calling ModelingTools in your likelihood " - _error_msg += "module. \n The amplitude parameter is currently set to" - _error_msg += f" {_amplitude_param} and Firecrown won't let Smokescreen change that." - - # check if the required parameter is in the cosmology - if _required_param not in self.cosmo.to_dict().keys(): - error_msg = f"Cosmology does not have the required parameter {_required_param}" - error_msg += _error_msg - raise ValueError(error_msg) - - # check if the required parameter is in the shifts dictionary - if any(param in self.shifts_dict for param in ['A_s', 'sigma8']): - if _required_param not in self.shifts_dict.keys(): - error_msg = "Shifts dictionary does not have the required parameter " - error_msg += f"{_required_param}" - error_msg += _error_msg - raise ValueError(error_msg) - - def _load_default_systematics(self, likelihood): - """ - Loads the default systematics from the likelihood. - - Parameters - ---------- - likelihood : firecrown.likelihood.Likelihood - Likelihood object with required systematics. - - Returns - ------- - ParamsMap : firecrown.parameters.ParamsMap - A ParamsMap object with the default values of the required systematics. - """ - # get the required systematics from the likelihood - req_systematics = likelihood.required_parameters() - default_systematics = req_systematics.get_default_values() - return ParamsMap(default_systematics) - - def _load_systematics(self, systematics_dict, likelihood): - """ - Loads the systematics from the systematics dictionary. - - Parameters - ---------- - systematics_dict : dict - Dictionary of systematics names and corresponding fiducial values. - - likelihood : firecrown.likelihood.Likelihood - Likelihood object with required systematics. - - Returns - ------- - ParamsMap : firecrown.parameters.ParamsMap - A ParamsMap object with the systematics values from the dictionary. - """ - likelihood_req_systematics = list(likelihood.required_parameters().get_params_names()) if self._debug: - print(f"[DEBUG] Likelihood requires systematics: {likelihood_req_systematics}") - # test if all keys in the systematics_dict are in the likelihood systematics: - for key in likelihood_req_systematics: - if key not in systematics_dict.keys(): - raise ValueError(f"Systematic {key} not in likelihood systematics") - return ParamsMap(systematics_dict) - - def _load_shifts(self, seed, shift_distr="flat"): - """ - Loads the shifts from the shifts dictionary. - - Parameters - ---------- - seed : int or str - Seed for the random number generator. If a string is provided, - it is converted to an integer using a hash function. - - shifts_dict : dict - Dictionary of parameter names and corresponding shift widths. If the - shifts are single values, it does a deterministic shift: PARAM = FIDUCIAL + SHIFT - If the shifts are tuples of values, the dictionary values - should be the (lower, upper) bounds of the shift widths: PARAM = U(a, b) - If the first valuee is negative, it is assumed that the parameter - is to be shifted from the fiducial value: PARAM = FIDUCIAL + U(-a, b) - - Returns - ------- - dict - Dictionary of parameter names and corresponding shifts. - """ - if shift_distr == "flat": - shifts_internal = draw_flat_or_deterministic_param_shifts(self.cosmo, self.shifts_dict, seed) - return shifts_internal - elif shift_distr == "gaussian": - return draw_gaussian_param_shifts(self.cosmo, self.shifts_dict, seed) - else: - raise NotImplementedError('Only flat and gaussian shifts are implemented') - - def _create_concealed_cosmo(self, shifts): - """ - Creates a blinded cosmology object with the shifts applied. - - Parameters - ---------- - shifts : dict - Dictionary of parameter names and corresponding shifts. + print(f"[DEBUG] Shifts: {self.__shifts}") + print(f"[DEBUG] Concealed params: {self.__concealed_params}") - Returns - ------- - pyccl.Cosmology - A Cosmology object with the shifts applied. - """ - concealed_cosmo_dict = deepcopy(self.cosmo.to_dict()) - # sometimes we have this extra paramters that can cause problems: - try: - del concealed_cosmo_dict['extra_parameters'] - except KeyError: # pragma: no cover - pass - for k in shifts.keys(): - concealed_cosmo_dict[k] = shifts[k] - concealed_cosmo = ccl.Cosmology(**concealed_cosmo_dict) - return concealed_cosmo + def _overlay_shifts(self, shifts): + """Overlay drawn deltas on the fiducial point: fiducial[k] + shift[k].""" + concealed = deepcopy(self.fiducial_params) + for k, delta in shifts.items(): + concealed[k] = self.fiducial_params[k] + delta + return concealed def calculate_concealing_factor(self, factor_type="add"): r""" - Calculates the concealing (blinding) factor for the data-vector, - according to Muir et al. 2019: + Calculate the concealing (blinding) factor, per Muir et al. 2019. Parameters ---------- factor_type : str - Type of concealing (blinding) factor to be calculated. Default is ``add``. + ``"add"`` (default) or ``"mult"``. Returns ------- - .__concealing_factor : np.ndarray - Concealing factor. + np.ndarray + Concealing factor (returned only in debug mode). Notes ----- type="add": - .. math:: - f^{\rm add} = d(\theta_{\rm blind}) - d(\theta_{\rm fid}) + .. math:: f^{\rm add} = t(\theta_{\rm hidden}) - t(\theta_{\rm fid}) type="mult": - .. math:: - f^{\rm mult} = d(\theta_{\rm blind}) / d(\theta_{\rm fid}) + .. math:: f^{\rm mult} = t(\theta_{\rm hidden}) / t(\theta_{\rm fid}) """ self.factor_type = factor_type - # need to get the defaults from firecrown: - _firecrown_defaults = get_default_params_map(self.tools, self.likelihood) - - _params_reference = modify_default_params(_firecrown_defaults, - self.cosmo.to_dict(), - self.systematics_dict) - # update the tools: - self.tools.update(_params_reference) - # prepare the original cosmology tools: - self.tools.prepare() - # update the likelihood with the systematics parameters: - self.likelihood.update(_params_reference) - # fiducial theory vector: - self.theory_vec_fid = self.likelihood.compute_theory_vector(self.tools) - # resets the likelihood and tools - self.likelihood.reset() - self.tools.reset() - - # now calculates the shifted theory vector: - # updating the default params with the concealed cosmology: - __params_concealed = modify_default_params(_firecrown_defaults, - self.__concealed_cosmo.to_dict(), - self.systematics_dict) - # update the tools: - self.tools.update(__params_concealed) - # prepare the original cosmology tools: - self.tools.prepare() - # update the likelihood with the systematics parameters: - self.likelihood.update(__params_concealed) - # concealed theory vector: - self.theory_vec_conceal = self.likelihood.compute_theory_vector(self.tools) + # fiducial vector computed once in __init__ (length guard); recompute + # nothing there --- reuse it and evaluate only the hidden point here. + self.theory_vec_conceal = np.asarray(self.theory_fn(self.__concealed_params)) - if self.factor_type == "add": + if factor_type == "add": self.__concealing_factor = self.theory_vec_conceal - self.theory_vec_fid - elif self.factor_type == "mult": + elif factor_type == "mult": self.__concealing_factor = self.theory_vec_conceal / self.theory_vec_fid else: raise NotImplementedError('Only "add" and "mult" concealing factor is implemented') @@ -509,21 +167,15 @@ def calculate_concealing_factor(self, factor_type="add"): def apply_concealing_to_likelihood_datavec(self): r""" - Applies the concealing (blinding) factor to the data-vector. + Apply the concealing (blinding) factor to the SACC data vector. Returns ------- - self.concealed_data_vector : np.ndarray - Concealed (blinded) data-vector. - - Notes - ----- - The data-vector is concealed by adding or multiplying the - concealing factor to the data-vector, depending on the - `factor_type` specified during the calculation of the - concealing factor. + np.ndarray + Concealed (blinded) data vector: ``sacc_data.mean + factor`` + (``"add"``) or ``sacc_data.mean * factor`` (``"mult"``). """ - self.data_vector = self.likelihood.get_data_vector() + self.data_vector = self.sacc_data.mean if self.factor_type == "add": self.concealed_data_vector = self.data_vector + self.__concealing_factor elif self.factor_type == "mult": @@ -536,50 +188,49 @@ def save_concealed_datavector(self, path_to_save, file_root, return_sacc=False, output_format=None, suffix=None): """ - Saves the concealed (blinded) data-vector to a file. + Save the concealed (blinded) data vector to a SACC file. - Saves the blinded data-vector to a file with the appropriate extension - based on the input format: ``.fits`` for FITS files or ``.hdf5`` for HDF5 files. + The blinded vector overwrites the mean of a deep copy of ``sacc_data``; + the covariance is carried over unchanged (blinding shifts the mean and + never touches the covariance). Metadata is stamped with the concealed + flag, creator, timestamp, and seed. Writing goes through SACC's own + ``save_fits`` / ``save_hdf5``. Parameters ---------- path_to_save : str - Path to save the blinded data-vector. + Directory to save the blinded data vector. file_root : str Root of the file name. return_sacc : bool - If True, returns the sacc object with the blinded data-vector. + If True, returns the sacc object with the blinded data vector. output_format : str, optional - Output format to use. If None, uses the detected input format. + Output format ('fits' or 'hdf5'). Defaults to the detected input + format. suffix : str, optional - Suffix for the output file name. Defaults to 'concealed_data_vector'. + Suffix for the output file name. Defaults to + 'concealed_data_vector'. Returns ------- sacc.sacc.Sacc or None - If `return_sacc` is True, returns the sacc object with - the blinded data-vector. Otherwise, returns None. """ if suffix is None: suffix = "concealed_data_vector" - # Determine output format: use specified format or fall back to input format if output_format is None: output_format = getattr(self, '_input_format', 'fits') - idx = self.likelihood.get_sacc_indices() - concealed_sacc = save_to_sacc(self.sacc_data, - self.concealed_data_vector, - idx) - # copies the metadata from the original sacc file: - concealed_sacc.metadata = self.sacc_data.metadata - # adds metadata to the sacc file: + concealed_sacc = deepcopy(self.sacc_data) + concealed_sacc.mean = self.concealed_data_vector + + # copy metadata from the original sacc file, then stamp custody info: + concealed_sacc.metadata = dict(self.sacc_data.metadata) concealed_sacc.metadata['concealed'] = True concealed_sacc.metadata['creator'] = getpass.getuser() concealed_sacc.metadata['creation'] = datetime.datetime.now().isoformat() concealed_sacc.metadata['info'] = 'Concealed (blinded) data-vector, created by Smokescreen.' concealed_sacc.metadata['seed_smokescreen'] = self.seed - # Determine file extension based on format if output_format == 'hdf5': ext = '.hdf5' save_method = concealed_sacc.save_hdf5 @@ -591,5 +242,4 @@ def save_concealed_datavector(self, path_to_save, file_root, save_method(output_path, overwrite=True) if return_sacc: return concealed_sacc - else: - return None + return None diff --git a/src/smokescreen/param_shifts.py b/src/smokescreen/param_shifts.py index e6d91fd..baf648f 100644 --- a/src/smokescreen/param_shifts.py +++ b/src/smokescreen/param_shifts.py @@ -6,149 +6,111 @@ .. currentmodule:: smokescreen.param_shifts -The :mod:`smokescreen.param_shifts` module provides modules -to perform shifts in the cosmological parameters. +The :mod:`smokescreen.param_shifts` module draws the hidden parameter +shifts that drive blinding. The shift for a given parameter is a *delta* +to be added to the fiducial value, drawn from a local, seeded RNG in a +way that does not depend on the iteration order of the shifts mapping. Parameter Shifts ----------------- -.. autofunction:: draw_flat_param_shifts -.. autofunction:: draw_flat_or_deterministic_param_shifts -.. autofunction:: draw_gaussian_param_shifts +.. autofunction:: draw_param_shifts ''' +import hashlib +from collections.abc import Mapping + import numpy as np -from .utils import string_to_seed -def draw_flat_param_shifts(shift_dict, seed): +def _normalize_seed(seed): """ - Draw flat parameter shifts from a dictionary of parameter names and - corresponding shift widths. - - Parameters - ---------- - shift_dict : dict - Dictionary of parameter names and corresponding shift widths. If the - shifts are single values, the dictionary values should be the shift - widths. If the shifts are tuples of values, the dictionary values - should be the (lower, upper) bounds of the shift widths. - seed : int or str - Random seed. + Reduce a seed to an integer suitable for ``numpy.random.default_rng``. - Returns - ------- - dict - Dictionary of parameter names and corresponding flat parameter shifts. - """ - if type(seed) is str: - seed = string_to_seed(seed) - np.random.seed(seed) - # check the if the shifts are single value or a tuple of values - if type(list(shift_dict.values())[0]) is tuple: - return {par: np.random.uniform(shift_dict[par][0], shift_dict[par][1]) - for par in shift_dict} - else: - return {par: np.random.uniform(-shift_dict[par], shift_dict[par]) - for par in shift_dict} - - -def draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed): - """ - Draw flat or deterministic parameter shifts from a dictionary of parameter - names and corresponding shift widths. + An ``int`` is returned as-is. A ``str`` is reduced deterministically to a + 64-bit integer via SHA-256 of its UTF-8 bytes; the mapping from a given + string to its integer seed is fixed and reproducible across runs and + machines. Parameters ---------- - cosmo : pyccl.Cosmology - Cosmology object. - shift_dict : dict - Dictionary of parameter names and corresponding shift. If the - shifts are single values, it does a deterministic shift: PARAM = FIDUCIAL + SHIFT - If the shifts are tuples of values, the dictionary values - should be the (lower, upper) bounds of the shift widths: PARAM = U(a, b) seed : int or str - Random seed. + Seed value. Returns ------- - dict - Dictionary of parameter names and corresponding flat or deterministic - parameter shifts. + int + Normalized integer seed. """ - if type(seed) is str: - seed = string_to_seed(seed) - np.random.seed(seed) - - # check if the keys in the shifts_dict are in the cosmology parameters - for key in shifts_dict.keys(): - try: - cosmo._params[key] - except (AttributeError, KeyError) as error: - # remove the key from the shifts_dict - # print(f"Key {key} not in cosmology parameters") - # failed_keys.append(key) - raise ValueError(f"[{error}]Key {key} not in cosmology parameters") - shifts = {} - # we loop over the ccl dict keys to ensure params are drawn in the same order! - for key in cosmo.to_dict().keys(): - if key in shifts_dict.keys(): - if isinstance(shifts_dict[key], tuple): - # check if the tuple is of length 2 - if len(shifts_dict[key]) == 2: - shifts[key] = np.random.uniform(shifts_dict[key][0], shifts_dict[key][1]) - else: - raise ValueError(f"Tuple {shifts_dict[key]} has to be of length 2") - else: - shifts[key] = shifts_dict[key] - else: - pass - return shifts + if isinstance(seed, str): + digest = hashlib.sha256(seed.encode("utf-8")).digest() + return int.from_bytes(digest[:8], "big") + if isinstance(seed, (int, np.integer)): + return int(seed) + raise TypeError(f"seed must be int or str, got {type(seed).__name__}") -def draw_gaussian_param_shifts(cosmo, shifts_dict, seed): +def draw_param_shifts(shifts_dict, seed, shift_distr="flat"): """ - Draw Gaussian parameter shifts from a dictionary of parameter names and - corresponding shift widths. + Draw a delta for each named parameter from a local RNG seeded by ``seed``. + + Each returned value is a shift to be *added* to the fiducial value, not an + absolute parameter value. The draw is order-independent: the delta for a + given key depends only on ``(key, seed, shift_distr)``, not on the + iteration order of ``shifts_dict``. This is achieved by iterating the keys + in sorted order and consuming exactly one scalar RNG call per key. Parameters ---------- - cosmo : pyccl.Cosmology - Cosmology object. - shift_dict : dict - Dictionary of parameter names and corresponding shift widths. The - dictionary values should be the (mean, std) of the Gaussian distribution. + shifts_dict : Mapping[str, float | tuple[float, float]] + Mapping of parameter name to its shift envelope, interpreted as a delta + about zero. A ``float`` ``h`` is the symmetric delta envelope + ``(-h, +h)``; a ``(lo, hi)`` tuple is a delta box, the drawn delta + lies in ``[lo, hi]`` (typically straddling zero). seed : int or str - Random seed. + Seed for the random number generator. A ``str`` is normalized to an + ``int`` (see :func:`_normalize_seed`). + shift_distr : str + ``"flat"`` draws uniformly over the delta envelope. ``"gaussian"`` + draws from a zero-mean Gaussian with ``sigma`` equal to the ``float`` + half-width and accepts only the symmetric ``float`` form. Returns ------- - shifts : dict - Dictionary of parameter names and corresponding Gaussian parameter shifts. + dict[str, float] + Mapping ``{param_name: drawn_delta}``, keys a subset of + ``shifts_dict``'s keys, values plain floats. Overlay these on the + fiducial parameters: ``concealed[k] = fiducial[k] + delta[k]``. """ - if type(seed) is str: - seed = string_to_seed(seed) - np.random.seed(seed) - for key in shifts_dict.keys(): - try: - cosmo._params[key] - except (AttributeError, KeyError) as error: - # remove the key from the shifts_dict - # print(f"Key {key} not in cosmology parameters") - # failed_keys.append(key) - raise ValueError(f"[{error}]Key {key} not in cosmology parameters") + if shift_distr not in ("flat", "gaussian"): + raise NotImplementedError("Only flat and gaussian shifts are implemented") + + rng = np.random.default_rng(_normalize_seed(seed)) + shifts = {} - # we loop over the ccl dict keys to ensure params are drawn in the same order! - for key in cosmo.to_dict().keys(): - if key in shifts_dict.keys(): - if isinstance(shifts_dict[key], tuple): - # check if the tuple is of length 2 - if len(shifts_dict[key]) == 2: - shifts[key] = np.random.normal(shifts_dict[key][0], shifts_dict[key][1]) - else: - raise ValueError(f"Tuple {shifts_dict[key]} has to be of length 2") + # Sorted iteration makes the per-key draw independent of dict order. + for key in sorted(shifts_dict): + envelope = shifts_dict[key] + if shift_distr == "flat": + if isinstance(envelope, tuple): + if len(envelope) != 2: + raise ValueError(f"Tuple {envelope} has to be of length 2") + lo, hi = envelope else: - raise ValueError(f"Value {shifts_dict[key]} has to be a tuple of length 2") - else: - pass + lo, hi = -envelope, envelope + shifts[key] = float(rng.uniform(lo, hi)) + else: # gaussian + if isinstance(envelope, tuple): + raise ValueError( + f"Gaussian shift for {key} requires a float half-width, " + f"got tuple {envelope}" + ) + shifts[key] = float(rng.normal(0.0, envelope)) return shifts + + +# Backwards-compatibility alias: the draw returns a plain mapping regardless of +# whether callers reach for it by the historical or the current name. +def draw_flat_param_shifts(shift_dict: Mapping, seed): + """Deprecated alias for :func:`draw_param_shifts` with ``shift_distr="flat"``.""" + return draw_param_shifts(shift_dict, seed, shift_distr="flat") diff --git a/src/smokescreen/utils.py b/src/smokescreen/utils.py index c932f5f..d8f2574 100644 --- a/src/smokescreen/utils.py +++ b/src/smokescreen/utils.py @@ -20,8 +20,6 @@ ''' import hashlib import importlib.util -import pyccl as ccl -import sacc def load_cosmology_from_partial_dict(cosmo_dict): @@ -39,6 +37,8 @@ def load_cosmology_from_partial_dict(cosmo_dict): Cosmology Cosmology object with the specified parameters. """ + import pyccl as ccl + # sets the default values for the cosmological parameters cosmo_dict_default = ccl.CosmologyVanillaLCDM().to_dict() @@ -116,34 +116,7 @@ def string_to_seed(seedstring): return int(int(hashlib.md5(seedstring.encode('utf-8')).hexdigest(), 16) % 1.e8) -def modify_default_params(default_params, ccl_cosmology, systematics=None): - """ - Modify the default parameters with the values from the CCL cosmology and - systematics if provided. - - Parameters - ---------- - default_params : dict - Dictionary with the default parameters. - ccl_cosmology : dict - Dictionary with the CCL cosmology parameters. - systematics : dict, optional - Dictionary with the systematics parameters to override the defaults. - - Returns - ------- - dict - Dictionary with the modified default parameters. - """ - for key, value in default_params.items(): - if key in ccl_cosmology: - default_params[key] = ccl_cosmology[key] - elif systematics is not None and key in systematics: - default_params[key] = systematics[key] - return default_params - - -def load_sacc_file(path_to_sacc: str) -> tuple[sacc.Sacc, str]: +def load_sacc_file(path_to_sacc: str): """ Load a SACC file, automatically detecting if it's FITS or HDF5 format. @@ -166,6 +139,8 @@ def load_sacc_file(path_to_sacc: str) -> tuple[sacc.Sacc, str]: ValueError If the file cannot be loaded as either format """ + import sacc + # Try HDF5 first (more specific format check) try: sacc_obj = sacc.Sacc.load_hdf5(path_to_sacc) From 76ded3f267a6a259d790af1b930d015eb8032b14 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 23:51:19 +0200 Subject: [PATCH 2/7] Retain the inherited firecrown concealment path as a lazy, separate module MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Preserve upstream's firecrown-likelihood blinding class (pre-refactor ConcealDataVector) as FirecrownConcealDataVector in its own module, so the fork does not amputate the upstream firecrown integration a later DESC PR might restore. All firecrown and pyccl imports are function-local (via a _import_firecrown() helper and local imports), so importing this module — or importing smokescreen — imports no firecrown. It is not wired into the default flow (no main module imports it) and no fork test exercises it; the module docstring records it as inherited-upstream and unsupported here. The supported path remains smokescreen.datavector.ConcealDataVector. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/smokescreen/firecrown_datavector.py | 435 ++++++++++++++++++++++++ 1 file changed, 435 insertions(+) create mode 100644 src/smokescreen/firecrown_datavector.py diff --git a/src/smokescreen/firecrown_datavector.py b/src/smokescreen/firecrown_datavector.py new file mode 100644 index 0000000..6864433 --- /dev/null +++ b/src/smokescreen/firecrown_datavector.py @@ -0,0 +1,435 @@ +# author: Arthur Loureiro +# license: BSD 3-Clause +''' +Firecrown concealment path (:mod:`smokescreen.firecrown_datavector`) +=================================================================== + +.. currentmodule:: smokescreen.firecrown_datavector + +Inherited-from-upstream firecrown-likelihood concealment. This module carries +the original DESC Smokescreen data-vector-blinding class, which drives theory +from a firecrown likelihood + CCL cosmology. + +**This path is inherited from upstream and is not supported in this fork.** +It is retained so the fork does not amputate upstream's firecrown integration. +It is *not* wired into the fork's default flow: nothing in +:mod:`smokescreen`'s main modules imports it, no test in this fork exercises +it, and firecrown is neither declared as a dependency nor installed. All +firecrown (and ``pyccl``) imports are function-local, so importing this module +does not import firecrown --- and importing :mod:`smokescreen` does not import +this module. The supported concealment path is +:class:`smokescreen.datavector.ConcealDataVector`, which drives theory from a +``theory_fn`` protocol and ships a default CCL backend. + + +Firecrown Conceal Data Vector Class +----------------------------------- + +.. autoclass:: FirecrownConcealDataVector + :members: + :undoc-members: +''' +import os +import types +import inspect +import datetime +import getpass +from copy import deepcopy + +import numpy as np +import sacc + + +def _import_firecrown(): + """ + Lazily import the firecrown symbols the concealment path needs. + + Kept function-local so that importing this module (and, transitively, + :mod:`smokescreen`) never imports firecrown. Firecrown is inherited from + upstream and unsupported in this fork; it is not installed here. + """ + from packaging.version import Version + import firecrown + + # Handle different Firecrown versions + if Version(firecrown.__version__) >= Version("1.15.0a0"): + from firecrown.likelihood import ( + load_likelihood, + load_likelihood_from_module_type, + NamedParameters, + ) + else: + from firecrown.likelihood.likelihood import ( + load_likelihood, + load_likelihood_from_module_type, + NamedParameters, + ) + from firecrown.parameters import ParamsMap + from firecrown.updatable import get_default_params_map + from firecrown.utils import save_to_sacc + from firecrown.ccl_factory import PoweSpecAmplitudeParameter + + return { + "load_likelihood": load_likelihood, + "load_likelihood_from_module_type": load_likelihood_from_module_type, + "NamedParameters": NamedParameters, + "ParamsMap": ParamsMap, + "get_default_params_map": get_default_params_map, + "save_to_sacc": save_to_sacc, + "PoweSpecAmplitudeParameter": PoweSpecAmplitudeParameter, + } + + +def _load_module_from_path(path): + """Load a module from a given filesystem path.""" + import importlib.util + + spec = importlib.util.spec_from_file_location("module.name", path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def _modify_default_params(default_params, ccl_cosmology, systematics=None): + """Override firecrown defaults with the CCL cosmology and systematics.""" + for key in default_params: + if key in ccl_cosmology: + default_params[key] = ccl_cosmology[key] + elif systematics is not None and key in systematics: + default_params[key] = systematics[key] + return default_params + + +def _draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed): + """ + Inherited draw: single values are deterministic absolute shifts, tuples are + absolute ``(lo, hi)`` uniform bounds. Validated against CCL parameter names. + Uses the process-global ``np.random`` for reproducibility with upstream. + """ + from smokescreen.utils import string_to_seed + + if type(seed) is str: + seed = string_to_seed(seed) + np.random.seed(seed) + + for key in shifts_dict.keys(): + try: + cosmo._params[key] + except (AttributeError, KeyError) as error: + raise ValueError(f"[{error}]Key {key} not in cosmology parameters") + shifts = {} + # loop over the ccl dict keys to ensure params are drawn in the same order + for key in cosmo.to_dict().keys(): + if key in shifts_dict.keys(): + if isinstance(shifts_dict[key], tuple): + if len(shifts_dict[key]) == 2: + shifts[key] = np.random.uniform(shifts_dict[key][0], shifts_dict[key][1]) + else: + raise ValueError(f"Tuple {shifts_dict[key]} has to be of length 2") + else: + shifts[key] = shifts_dict[key] + return shifts + + +def _draw_gaussian_param_shifts(cosmo, shifts_dict, seed): + """ + Inherited Gaussian draw: ``(mean, std)`` tuples, validated against CCL + parameter names. Uses the process-global ``np.random``. + """ + from smokescreen.utils import string_to_seed + + if type(seed) is str: + seed = string_to_seed(seed) + np.random.seed(seed) + + for key in shifts_dict.keys(): + try: + cosmo._params[key] + except (AttributeError, KeyError) as error: + raise ValueError(f"[{error}]Key {key} not in cosmology parameters") + shifts = {} + for key in cosmo.to_dict().keys(): + if key in shifts_dict.keys(): + if isinstance(shifts_dict[key], tuple): + if len(shifts_dict[key]) == 2: + shifts[key] = np.random.normal(shifts_dict[key][0], shifts_dict[key][1]) + else: + raise ValueError(f"Tuple {shifts_dict[key]} has to be of length 2") + else: + raise ValueError(f"Value {shifts_dict[key]} has to be a tuple of length 2") + return shifts + + +class FirecrownConcealDataVector(): + """ + Inherited-from-upstream concealment driven by a firecrown likelihood. + + **Unsupported in this fork** --- retained from upstream, not tested, not + installed, not wired into the default flow. Use + :class:`smokescreen.datavector.ConcealDataVector` instead. + + FIXME: Only cosmological parameters are supported for now for the shifts + + Parameters + ---------- + cosmo : pyccl.Cosmology + Cosmology object from CCL with a fiducial cosmology. + likelihood : str or module + path to the likelihood or a module containing the likelihood + must contain both `build_likelihood` and `compute_theory_vector` methods + shifts_dict : dict + Dictionary of parameter names and corresponding shift widths. + sacc_data : sacc.sacc.Sacc + Data-vector to be concealed (blinded). + systm_dict : dict + Dictionary of systematics names and corresponding fiducial values. + seed : int or str + Random seed. + + Keyword Arguments + ----------------- + shift_distr : str + Type of shift to be applied. Default is "flat". + debug : bool + If True, prints debug information. Default is False. + """ + def __init__(self, cosmo, likelihood, shifts_dict, sacc_data, systm_dict=None, + seed="2112", **kwargs): + self.cosmo = cosmo + self.systematics_dict = systm_dict + self.sacc_data = sacc_data + self.seed = seed + assert isinstance(self.sacc_data, sacc.sacc.Sacc), "sacc_data must be a sacc object" + self.shifts_dict = shifts_dict + + # detect original file format for output preservation + if 'input_format' in kwargs: + self._input_format = kwargs['input_format'] + else: + self._input_format = getattr(sacc_data, '_smokescreen_input_format', 'fits') + + # load the shifts + if 'shift_distr' in kwargs: + self.__shifts = self._load_shifts(seed, shift_distr=kwargs['shift_distr']) + else: + self.__shifts = self._load_shifts(seed) + + # create concealed cosmology object: + self.__concealed_cosmo = self._create_concealed_cosmo(self.__shifts) + + if 'debug' in kwargs and kwargs['debug']: + self._debug = True + print(f"[DEBUG] Shifts: {self.__shifts}") + print(f"[DEBUG] Concealed Cosmology: {self.__concealed_cosmo}") + else: + self._debug = False + + # load the likelihood + self.likelihood, self.tools = self._load_likelihood(likelihood, + self.sacc_data) + + # load the systematics + if self.systematics_dict is None: + self.systematics = self._load_default_systematics(self.likelihood) + else: + self.systematics = self._load_systematics(self.systematics_dict, self.likelihood) + + def _load_likelihood(self, likelihood, sacc_data): + """ + Loads the likelihood either from a python module or from a file. + """ + fc = _import_firecrown() + build_parameters = fc["NamedParameters"]({'sacc_data': sacc_data}) + + if type(likelihood) is str: + if not os.path.isfile(likelihood): + raise FileNotFoundError(f'Could not find file {likelihood}') + self._test_likelihood(likelihood, 'str') + likelihood, tools = fc["load_likelihood"](likelihood, build_parameters) + self._check_amplitude_parameter(tools) + if not hasattr(likelihood, 'compute_theory_vector'): # pragma: no cover + raise AttributeError('Likelihood does not have a compute_vector method') + return likelihood, tools + + elif isinstance(likelihood, types.ModuleType): + self._test_likelihood(likelihood, 'module') + likelihood, tools = fc["load_likelihood_from_module_type"](likelihood, + build_parameters) + self._check_amplitude_parameter(tools) + if not hasattr(likelihood, 'compute_theory_vector'): # pragma: no cover + raise AttributeError('Likelihood does not have a compute_vector method') + return likelihood, tools + else: + raise TypeError('Likelihood must be a string path to a likelihood module or a module') + + def _test_likelihood(self, likelihood, like_type): + """Tests if the likelihood has the required methods.""" + if like_type == "str": + likelihood = _load_module_from_path(likelihood) + + if not hasattr(likelihood, 'build_likelihood'): + raise AttributeError('Likelihood does not have a build_likelihood method') + + if self.sacc_data is not None: + sig = inspect.signature(likelihood.build_likelihood) + likefunc_params = sig.parameters + assert len(likefunc_params) >= 1, ("A sacc was provided, ", + "the likelihood must require a", + "build_parameters NamedParameters object!") + + def _check_amplitude_parameter(self, tools): + """Checks the amplitude parameter set in the tools matches the concealed one.""" + fc = _import_firecrown() + PoweSpecAmplitudeParameter = fc["PoweSpecAmplitudeParameter"] + _amplitude_param = tools.ccl_factory.amplitude_parameter + + if _amplitude_param is PoweSpecAmplitudeParameter.SIGMA8: + _required_param = 'sigma8' + elif _amplitude_param is PoweSpecAmplitudeParameter.AS: + _required_param = 'A_s' + else: + raise ValueError(f"Amplitude parameter {_amplitude_param} not supported") + + _error_msg = "\n You probably need to set the amplitude parameter [A_s/sigma8] " + _error_msg += "that you want to conceal when calling ModelingTools in your likelihood " + _error_msg += "module. \n The amplitude parameter is currently set to" + _error_msg += f" {_amplitude_param} and Firecrown won't let Smokescreen change that." + + if _required_param not in self.cosmo.to_dict().keys(): + error_msg = f"Cosmology does not have the required parameter {_required_param}" + error_msg += _error_msg + raise ValueError(error_msg) + + if any(param in self.shifts_dict for param in ['A_s', 'sigma8']): + if _required_param not in self.shifts_dict.keys(): + error_msg = "Shifts dictionary does not have the required parameter " + error_msg += f"{_required_param}" + error_msg += _error_msg + raise ValueError(error_msg) + + def _load_default_systematics(self, likelihood): + """Loads the default systematics from the likelihood.""" + fc = _import_firecrown() + req_systematics = likelihood.required_parameters() + default_systematics = req_systematics.get_default_values() + return fc["ParamsMap"](default_systematics) + + def _load_systematics(self, systematics_dict, likelihood): + """Loads the systematics from the systematics dictionary.""" + fc = _import_firecrown() + likelihood_req_systematics = list(likelihood.required_parameters().get_params_names()) + if self._debug: + print(f"[DEBUG] Likelihood requires systematics: {likelihood_req_systematics}") + for key in likelihood_req_systematics: + if key not in systematics_dict.keys(): + raise ValueError(f"Systematic {key} not in likelihood systematics") + return fc["ParamsMap"](systematics_dict) + + def _load_shifts(self, seed, shift_distr="flat"): + """Loads the shifts from the shifts dictionary.""" + if shift_distr == "flat": + return _draw_flat_or_deterministic_param_shifts(self.cosmo, self.shifts_dict, seed) + elif shift_distr == "gaussian": + return _draw_gaussian_param_shifts(self.cosmo, self.shifts_dict, seed) + else: + raise NotImplementedError('Only flat and gaussian shifts are implemented') + + def _create_concealed_cosmo(self, shifts): + """Creates a blinded cosmology object with the shifts applied.""" + import pyccl as ccl + + concealed_cosmo_dict = deepcopy(self.cosmo.to_dict()) + try: + del concealed_cosmo_dict['extra_parameters'] + except KeyError: # pragma: no cover + pass + for k in shifts.keys(): + concealed_cosmo_dict[k] = shifts[k] + return ccl.Cosmology(**concealed_cosmo_dict) + + def calculate_concealing_factor(self, factor_type="add"): + r""" + Calculates the concealing (blinding) factor, per Muir et al. 2019. + + Notes + ----- + type="add": :math:`f^{\rm add} = d(\theta_{\rm blind}) - d(\theta_{\rm fid})` + + type="mult": :math:`f^{\rm mult} = d(\theta_{\rm blind}) / d(\theta_{\rm fid})` + """ + fc = _import_firecrown() + self.factor_type = factor_type + + _firecrown_defaults = fc["get_default_params_map"](self.tools, self.likelihood) + + _params_reference = _modify_default_params(_firecrown_defaults, + self.cosmo.to_dict(), + self.systematics_dict) + self.tools.update(_params_reference) + self.tools.prepare() + self.likelihood.update(_params_reference) + self.theory_vec_fid = self.likelihood.compute_theory_vector(self.tools) + self.likelihood.reset() + self.tools.reset() + + __params_concealed = _modify_default_params(_firecrown_defaults, + self.__concealed_cosmo.to_dict(), + self.systematics_dict) + self.tools.update(__params_concealed) + self.tools.prepare() + self.likelihood.update(__params_concealed) + self.theory_vec_conceal = self.likelihood.compute_theory_vector(self.tools) + + if self.factor_type == "add": + self.__concealing_factor = self.theory_vec_conceal - self.theory_vec_fid + elif self.factor_type == "mult": + self.__concealing_factor = self.theory_vec_conceal / self.theory_vec_fid + else: + raise NotImplementedError('Only "add" and "mult" concealing factor is implemented') + if self._debug: + return self.__concealing_factor + + def apply_concealing_to_likelihood_datavec(self): + r"""Applies the concealing (blinding) factor to the data-vector.""" + self.data_vector = self.likelihood.get_data_vector() + if self.factor_type == "add": + self.concealed_data_vector = self.data_vector + self.__concealing_factor + elif self.factor_type == "mult": + self.concealed_data_vector = self.data_vector * self.__concealing_factor + else: + raise NotImplementedError('Only "add" and "mult" blinding factor is implemented') + return self.concealed_data_vector + + def save_concealed_datavector(self, path_to_save, file_root, + return_sacc=False, output_format=None, + suffix=None): + """Saves the concealed (blinded) data-vector to a file.""" + fc = _import_firecrown() + if suffix is None: + suffix = "concealed_data_vector" + if output_format is None: + output_format = getattr(self, '_input_format', 'fits') + + idx = self.likelihood.get_sacc_indices() + concealed_sacc = fc["save_to_sacc"](self.sacc_data, + self.concealed_data_vector, + idx) + concealed_sacc.metadata = self.sacc_data.metadata + concealed_sacc.metadata['concealed'] = True + concealed_sacc.metadata['creator'] = getpass.getuser() + concealed_sacc.metadata['creation'] = datetime.datetime.now().isoformat() + concealed_sacc.metadata['info'] = 'Concealed (blinded) data-vector, created by Smokescreen.' + concealed_sacc.metadata['seed_smokescreen'] = self.seed + + if output_format == 'hdf5': + ext = '.hdf5' + save_method = concealed_sacc.save_hdf5 + else: # default to FITS + ext = '.fits' + save_method = concealed_sacc.save_fits + + output_path = f"{path_to_save}/{file_root}_{suffix}{ext}" + save_method(output_path, overwrite=True) + if return_sacc: + return concealed_sacc + return None From d00c7a7598b9131693586939da083024bc7d7fe4 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 23:52:20 +0200 Subject: [PATCH 3/7] Migrate the test suite to the theory_fn protocol; drop firecrown tests Mechanics tests (construction, draw, concealing factor, apply, save) now run against a synthetic theory_fn closure, importing neither pyccl nor firecrown. Adds the order-independence / delta / local-RNG / string-seed draw tests, subprocess import-discipline tests (pyccl and firecrown absent from sys.modules after a fresh `import smokescreen`; pyccl present only once the CCL backend is built), a small default-CCL-backend end-to-end test, and the length-guard test. CLI tests updated to the new datavector_main signature. Every firecrown/likelihood/_verify_sacc_consistency test and the mock-likelihood fixture are removed; no test depends on firecrown. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- tests/test_data/mock_likelihood.py | 18 - tests/test_datavector.py | 1188 +++++----------------------- tests/test_import_discipline.py | 38 + tests/test_main.py | 163 ++-- tests/test_param_shifts.py | 185 ++--- 5 files changed, 333 insertions(+), 1259 deletions(-) delete mode 100644 tests/test_data/mock_likelihood.py create mode 100644 tests/test_import_discipline.py diff --git a/tests/test_data/mock_likelihood.py b/tests/test_data/mock_likelihood.py deleted file mode 100644 index f3eb9d8..0000000 --- a/tests/test_data/mock_likelihood.py +++ /dev/null @@ -1,18 +0,0 @@ -import types -import pyccl as ccl -from firecrown.modeling_tools import ModelingTools -from .test_datavector import EmptyLikelihood - -ccl.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION = False - -COSMO = ccl.CosmologyVanillaLCDM() - - -class MockLikelihoodModule(types.ModuleType): - def build_likelihood(self, *args, **kwargs): - self.mocktools = ModelingTools() - self.mocklike = EmptyLikelihood() - return self.mocklike, self.mocktools - - def compute_theory_vector(self, ModellingTools): - return self.mocklike.compute_theory_vector(ModellingTools) diff --git a/tests/test_datavector.py b/tests/test_datavector.py index bbd81ae..335fc54 100644 --- a/tests/test_datavector.py +++ b/tests/test_datavector.py @@ -1,1040 +1,212 @@ -import pytest # noqa: F401 -import types -import os -import datetime -from unittest.mock import patch, MagicMock # noqa: F401 -from packaging.version import Version +"""Blinding-mechanics tests for ConcealDataVector. + +All mechanics tests inject a SYNTHETIC ``theory_fn`` (a pure, deterministic +closure over the cosmo params) so they never import pyccl or firecrown. The one +default-CCL-backend test imports pyccl locally. +""" import numpy as np +import pytest import sacc -import pyccl as ccl -import firecrown -import shutil - -# Handle different Firecrown versions -if Version(firecrown.__version__) >= Version("1.15.0a0"): - from firecrown.likelihood import Likelihood -else: - from firecrown.likelihood.likelihood import Likelihood -from firecrown.modeling_tools import ModelingTools from smokescreen.datavector import ConcealDataVector -from smokescreen.utils import load_sacc_file +from smokescreen.param_shifts import draw_param_shifts -ccl.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION = False -COSMO = ccl.CosmologyVanillaLCDM() +FIDUCIAL = {"sigma8": 0.8, "Omega_c": 0.25} +SHIFTS = {"sigma8": 0.1, "Omega_c": (-0.05, 0.05)} +SEED = 2112 -@pytest.fixture -def fits_sacc_file(tmp_path): - """Create a FITS-format SACC file for testing.""" - sacc_data = sacc.Sacc() - # Add tracers matching cosmic shear example (src0, src1) - sacc_data.add_tracer('sp', 'src0') - sacc_data.add_tracer('sp', 'src1') - # Add data points with realistic values - n_ell = 5 - ell = np.logspace(np.log10(10), np.log10(1000), n_ell) - for e in ell: - sacc_data.add_data_point('galaxy_shear_cl_ee', ('src0', 'src0'), 1e-7, ell=e) - sacc_data.add_data_point('galaxy_shear_cl_ee', ('src0', 'src1'), 5e-8, ell=e) - sacc_data.add_data_point('galaxy_shear_cl_ee', ('src1', 'src1'), 2e-7, ell=e) +# --- fixtures --------------------------------------------------------------- - # Set mean (3 tracers * 5 ell values = 15 data points) - n_data = len(ell) * 3 - sacc_data.mean = np.ones(n_data) * 1e-7 - sacc_data.add_covariance(np.eye(n_data) * 1e-14) +def _make_sacc(n=5): + """A minimal SACC with n cl_ee rows and a dense covariance.""" + s = sacc.Sacc() + s.add_tracer("misc", "src") + for i in range(n): + s.add_data_point("galaxy_shear_cl_ee", ("src", "src"), + float(i + 1), ell=10 * (i + 1)) + cov = np.diag(np.arange(1, n + 1, dtype=float)) + 0.01 + s.add_covariance(cov) + return s - # Save to FITS file - fits_file = tmp_path / "test_sacc.fits" - sacc_data.save_fits(str(fits_file)) - return str(fits_file), sacc_data +def _synthetic_theory_fn(n): + """Pure length-n theory vector, deterministic in the cosmo params.""" + base = np.arange(1, n + 1, dtype=float) + return lambda p: base * p["sigma8"] + p["Omega_c"] @pytest.fixture -def hdf5_sacc_file(tmp_path): - """Create an HDF5-format SACC file for testing.""" - sacc_data = sacc.Sacc() - # Add tracers matching cosmic shear example (src0, src1) - sacc_data.add_tracer('sp', 'src0') - sacc_data.add_tracer('sp', 'src1') - # Add data points with realistic values - n_ell = 5 - ell = np.logspace(np.log10(10), np.log10(1000), n_ell) - for e in ell: - sacc_data.add_data_point('galaxy_shear_cl_ee', ('src0', 'src0'), 1e-7, ell=e) - sacc_data.add_data_point('galaxy_shear_cl_ee', ('src0', 'src1'), 5e-8, ell=e) - sacc_data.add_data_point('galaxy_shear_cl_ee', ('src1', 'src1'), 2e-7, ell=e) - - # Set mean (3 tracers * 5 ell values = 15 data points) - n_data = len(ell) * 3 - sacc_data.mean = np.ones(n_data) * 1e-7 - sacc_data.add_covariance(np.eye(n_data) * 1e-14) - - # Save to HDF5 file - hdf5_file = tmp_path / "test_sacc.hdf5" - sacc_data.save_hdf5(str(hdf5_file)) - - return str(hdf5_file), sacc_data +def sacc_data(): + return _make_sacc(5) @pytest.fixture -def cosmic_shear_resources(tmp_path): - """Copy cosmic shear example files to a temp directory for testing.""" - # Copy the likelihood file and SACC files - source_dir = os.path.join(os.path.dirname(__file__), '..', 'examples', 'cosmic_shear') - - # Copy likelihood file - likelihood_src = os.path.join(source_dir, 'cosmicshear_likelihood.py') - likelihood_dst = str(tmp_path / 'cosmicshear_likelihood.py') - shutil.copy2(likelihood_src, likelihood_dst) - - # Copy FITS SACC file - fits_src = os.path.join(source_dir, 'cosmicshear_sacc.fits') - fits_dst = str(tmp_path / 'cosmicshear_sacc.fits') - shutil.copy2(fits_src, fits_dst) - - # Copy HDF5 SACC file - hdf5_src = os.path.join(source_dir, 'cosmicshear_sacc.hdf5') - hdf5_dst = str(tmp_path / 'cosmicshear_sacc.hdf5') - shutil.copy2(hdf5_src, hdf5_dst) - - return { - 'likelihood': likelihood_dst, - 'fits_sacc': fits_dst, - 'hdf5_sacc': hdf5_dst - } - - -class EmptyLikelihood(Likelihood): - """ - empty mock likelihood based on: - https://github.com/LSSTDESC/firecrown/blob/master/tests/likelihood/lkdir/lkmodule.py - """ - def __init__(self): - self.nothing = 1.0 - self._data_vector = np.array([1.0, 2.0, 3.0]) - self._covariance = np.eye(3) * 0.1 - super().__init__() - - def read(self, sacc_data: sacc.Sacc): - pass - - def compute_loglike(self, ModellingTools): - return -self.nothing*2.0 - - def compute_theory_vector(self, ModellingTools): - return self.nothing - - def get_data_vector(self): - return self._data_vector - - def get_cov(self): - return self._covariance - - -class MockLikelihoodModule(types.ModuleType): - def build_likelihood(self, *args, **kwargs): - self.mocktools = ModelingTools() - self.mocklike = EmptyLikelihood() - return self.mocklike, self.mocktools - - def compute_theory_vector(self, ModellingTools): - return self.mocklike.compute_theory_vector(ModellingTools) - - -class MockCosmo: - def __init__(self, params): - self._params = params - - def __getitem__(self, key): - return self._params[key] +def theory_fn(): + return _synthetic_theory_fn(5) -def test_smokescreen_init(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(n_data) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} +# --- construction & length guard -------------------------------------------- - # Check that Smokescreen can be instantiated with valid inputs - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - assert isinstance(smokescreen, ConcealDataVector) +def test_construction_succeeds(sacc_data, theory_fn): + cdv = ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, + theory_fn=theory_fn) + assert isinstance(cdv, ConcealDataVector) - # Check that Smokescreen raises an error when given an invalid likelihood - with pytest.raises(AttributeError): - invalid_likelihood = types.ModuleType("invalid_likelihood") - ConcealDataVector(cosmo, invalid_likelihood, - shifts_dict, sacc_data, systematics_dict) - # check if breaks if given a shift with a key not in the cosmology parameters +def test_length_guard_raises(sacc_data): + wrong = lambda p: np.ones(3) # noqa: E731 -- wrong length with pytest.raises(ValueError): - invalid_shifts_dict = {"Omega_c": 1, "invalid_key": 1} - ConcealDataVector(cosmo, likelihood, - invalid_shifts_dict, sacc_data, systematics_dict) - - -def test_load_shifts(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1, "Omega_b": (-1, 2), "sigma8": (2, 3)} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - - # Call load_shifts and get the result - shifts = smokescreen._load_shifts(seed="2112") - - # Check that the shifts are correct - assert shifts["Omega_c"] == 1 - assert shifts["Omega_c"] >= 0 and shifts["Omega_b"] <= 3 - assert shifts["sigma8"] >= 2 and shifts["sigma8"] <= 3 - - # Check that an error is raised for an invalid tuple - smokescreen.shifts_dict["Omega_c"] = (1,) - with pytest.raises(ValueError): - smokescreen._load_shifts(seed="2112") - - # check if breaks if given a shift with a key not in the cosmology parameters - smokescreen.shifts_dict["invalid_key"] = 1 - with pytest.raises(ValueError): - smokescreen._load_shifts(seed="2112") - - # check if break if a a shift type is not flat - with pytest.raises(NotImplementedError): - smokescreen._load_shifts(seed="2112", shift_distr="invalid") - with pytest.raises(NotImplementedError): - smokescreen = ConcealDataVector(cosmo, likelihood, shifts_dict, sacc_data, - systematics_dict, - **{'shift_distr': 'invalid'}) - - -def test_load_shifts_gaussian(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": (0.3, 0.02), "Omega_b": (0.05, 0.002), "sigma8": (0.82, 0.02)} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict, - **{'shift_distr': 'gaussian'}) - - # Call load_shifts and get the result - shifts = smokescreen._load_shifts(seed="2112", shift_distr="gaussian") - - # Check that the shifts are correct - assert shifts["Omega_c"] >= 0.1 and shifts["Omega_c"] <= 0.4 - assert shifts["Omega_b"] >= 0.01 and shifts["Omega_b"] <= 0.05 - assert shifts["sigma8"] >= 0.5 and shifts["sigma8"] <= 1.2 - - -def test_verify_sacc_consistency_matching(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - - # Create a mock likelihood with matching data - mock_likelihood = MagicMock() - mock_likelihood.get_data_vector.return_value = np.array([1.0, 2.0, 3.0]) - mock_likelihood.get_cov.return_value = np.eye(3) * 0.1 - - # Test that no error is raised for matching data - smokescreen._verify_sacc_consistency(mock_likelihood) - - -def test_verify_sacc_consistency_mismatch_data_vector(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - - # Create a mock likelihood with mismatched data vector - mock_likelihood = MagicMock() - mock_likelihood.get_data_vector.return_value = np.array([2.0, 3.0, 4.0]) # Different values - mock_likelihood.get_cov.return_value = np.eye(3) * 0.1 - - # Test that ValueError is raised for mismatched data vector - with pytest.raises(ValueError) as exc_info: - smokescreen._verify_sacc_consistency(mock_likelihood) - - assert "Data vector mismatch" in str(exc_info.value) - - -def test_verify_sacc_consistency_mismatch_covariance(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - - # Create a mock likelihood with mismatched covariance - mock_likelihood = MagicMock() - mock_likelihood.get_data_vector.return_value = np.array([1.0, 2.0, 3.0]) - mock_likelihood.get_cov.return_value = np.eye(3) * 0.5 # Different variance - - # Test that ValueError is raised for mismatched covariance - with pytest.raises(ValueError) as exc_info: - smokescreen._verify_sacc_consistency(mock_likelihood) - - assert "Covariance matrix mismatch" in str(exc_info.value) - - -class EmptyLikelihoodNoCov(Likelihood): - """Empty mock likelihood that returns None for covariance.""" - def __init__(self): - self.nothing = 1.0 - self._data_vector = np.array([1.0, 2.0, 3.0]) - super().__init__() - - def read(self, sacc_data: sacc.Sacc): - pass - - def compute_loglike(self, ModellingTools): - return -self.nothing*2.0 - - def compute_theory_vector(self, ModellingTools): - return self.nothing - - def get_data_vector(self): - return self._data_vector - - def get_cov(self): - return None - - -def test_verify_sacc_consistency_none_covariance(): - # Create mock inputs where user has None covariance but likelihood has one - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - # Explicitly set covariance to None - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.covariance = None - - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Create a mock module that has build_likelihood returning EmptyLikelihoodNoCov - # (which returns None for covariance, so __init__ will succeed) - # Must accept build_parameters argument as required by _test_likelihood - mock_module = types.ModuleType("mock_likelihood") - mock_module.build_likelihood = lambda bp: (EmptyLikelihoodNoCov(), ModelingTools()) - - smokescreen = ConcealDataVector(cosmo, mock_module, - shifts_dict, sacc_data, systematics_dict) - - # Test that ValueError is raised when user has None but likelihood has covariance - # Use a mock likelihood with covariance (different from the one used in __init__) - mock_likelihood = MagicMock() - mock_likelihood.get_data_vector.return_value = np.array([1.0, 2.0, 3.0]) - mock_likelihood.get_cov.return_value = np.eye(3) * 0.1 - - with pytest.raises(ValueError) as exc_info: - smokescreen._verify_sacc_consistency(mock_likelihood) - - assert "Likelihood has covariance but user-provided SACC" in str(exc_info.value) - - -def test_verify_sacc_consistency_none_covariance_reverse(): - # Create mock inputs where user has covariance but likelihood returns None - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Add a proper covariance to the SACC object - n_data = 3 - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(n_data) * 0.1) - - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - - # Create a mock likelihood that returns None for covariance - mock_likelihood = MagicMock() - mock_likelihood.get_data_vector.return_value = np.array([1.0, 2.0, 3.0]) - mock_likelihood.get_cov.return_value = None - - # Test that ValueError is raised when likelihood has None but user has covariance - with pytest.raises(ValueError) as exc_info: - smokescreen._verify_sacc_consistency(mock_likelihood) - - assert "User-provided SACC has covariance but likelihood returns None for covariance." in str(exc_info.value) - - -def test_debug_mode(capfd): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Check that Smokescreen can be instantiated with valid inputs - _ = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict, - **{'debug': True}) - # Capture the output - out, err = capfd.readouterr() - - # Check that the debug output is correct - assert "[DEBUG] Shifts: " in out - assert "[DEBUG] Concealed Cosmology: " in out - assert f"{shifts_dict}" in out - - -def test_calculate_concealing_factor_add(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict, - **{'debug': True}) - - # Call calculate_concealing_factor with type="add" - concealing_factor = smokescreen.calculate_concealing_factor(factor_type="add") - - # Check that the concealing (blinding) factor is correct - assert concealing_factor == smokescreen.theory_vec_conceal - smokescreen.theory_vec_fid - - -def test_calculate_concealing_factor_add_gaussian(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": (0.3, 0.02), "Omega_b": (0.05, 0.002), "sigma8": (0.82, 0.02)} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict, - **{'debug': True, 'shift_distr': 'gaussian'}) - - # Call calculate_concealing_factor with type="add" - concealing_factor = smokescreen.calculate_concealing_factor(factor_type="add") - - # Check that the concealing (blinding) factor is correct - assert concealing_factor == smokescreen.theory_vec_conceal - smokescreen.theory_vec_fid - - -def test_calculate_concealing_factor_mult(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict, - **{'debug': True}) - - # Call calculate_concealing_factor with type="add" - concealing_factor = smokescreen.calculate_concealing_factor(factor_type="mult") - - # Check that the concealing (blinding) factor is correct - assert concealing_factor == smokescreen.theory_vec_conceal / smokescreen.theory_vec_fid - - -def test_calculate_concealing_factor_invalid_type(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - - # Call calculate_concealing_factor with an invalid type - with pytest.raises(NotImplementedError): - smokescreen.calculate_concealing_factor(factor_type="invalid") - - -def test_apply_concealing_to_likelihood_datavec_add(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict, - **{'debug': True}) - - # Set the concealing (blinding) factor and type - # Call calculate_concealing_factor with type="add" - concealing_factor = smokescreen.calculate_concealing_factor(factor_type="add") - - # Call apply_blinding_to_likelihood_datavec - concealed_data_vector = smokescreen.apply_concealing_to_likelihood_datavec() - expected_concealed = smokescreen.likelihood.get_data_vector() + concealing_factor - - # Check that the blinded data vector is correct - np.testing.assert_array_equal(concealed_data_vector, expected_concealed) - - -def test_apply_concealing_to_likelihood_datavec_mult(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict, - **{'debug': True}) - - # Set the concealing (blinding) factor and type - # Call calculate_concealing_factor with type="add" - concealing_factor = smokescreen.calculate_concealing_factor(factor_type="mult") - - # Call apply_concealing_to_likelihood_datavec - concealed_data_vector = smokescreen.apply_concealing_to_likelihood_datavec() - expected_concealing = smokescreen.likelihood.get_data_vector() * concealing_factor - - # Check that the concealing (blinding) data vector is correct - np.testing.assert_array_equal(concealed_data_vector, expected_concealing) - - -def test_apply_concealing_to_likelihood_datavec_invalid_type(): - # Create mock inputs - cosmo = COSMO - sacc_data = sacc.Sacc() - # Add a misc tracer and data points to set up the SACC object properly - sacc_data.add_tracer('misc', 'test') - n_data = 3 - for i in range(n_data): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - # Set mean to match EmptyLikelihood's get_data_vector() return value - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Instantiate Smokescreen - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict, - **{'debug': True}) - - # Set the expected_concealing factor and type - # Call calculate_concealing_factor with type="add" - _ = smokescreen.calculate_concealing_factor(factor_type="mult") - # Set an invalid type - smokescreen.factor_type = "invalid" - with pytest.raises(NotImplementedError): - smokescreen.apply_concealing_to_likelihood_datavec() - - -def test_load_likelihood(): - # Create mock inputs using a 3-element SACC file to match EmptyLikelihood - cosmo = COSMO - sacc_data = sacc.Sacc() - sacc_data.add_tracer('misc', 'test') - for i in range(3): - sacc_data.add_data_point('galaxy_shear_cl_ee', ('test', 'test'), 1.0, ell=10) - sacc_data.mean = np.array([1.0, 2.0, 3.0]) - sacc_data.add_covariance(np.eye(3) * 0.1) - - likelihood = MockLikelihoodModule("mock_likelihood") - systematics_dict = {"systematic1": 0.1} - shifts_dict = {"Omega_c": 1} - - # Create Smokescreen instance (this works because data vectors match) - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - - # Test with an invalid likelihood (neither a module nor a file path) - with pytest.raises(TypeError): - smokescreen._load_likelihood(123, sacc_data) - - # Test with a non-existent likelihood file path - with pytest.raises(FileNotFoundError): - smokescreen._load_likelihood("/path/to/nonexistent/file.py", sacc_data) - - # Test with a module that doesn't have a 'build_likelihood' method - invalid_likelihood = types.ModuleType("invalid_likelihood") - with pytest.raises(AttributeError): - smokescreen._load_likelihood(invalid_likelihood, sacc_data) - - # Test with a module that doesn't have a 'build_parameters' input - invalid_likelihood = types.ModuleType("invalid_likelihood") - invalid_likelihood.build_likelihood = lambda: None - with pytest.raises(AssertionError): - smokescreen._load_likelihood(invalid_likelihood, sacc_data) - - -@patch('src.smokescreen.datavector.getpass.getuser', return_value='test_user') -def test_save_concealed_datavector(mock_getuser, cosmic_shear_resources, tmp_path): - # Create mock inputs using fixture resources - cosmo = COSMO - likelihood = cosmic_shear_resources['likelihood'] - syst_dict = { - "trc1_delta_z": 0.1, - "trc0_delta_z": 0.1, - } - shift_dict = {"Omega_c": 0.34, "sigma8": 0.85} - sacc_data = sacc.Sacc.load_fits(cosmic_shear_resources['fits_sacc']) - sck = ConcealDataVector(cosmo, likelihood, - shift_dict, sacc_data, syst_dict, seed=1234) - - # Calculate the concealing factor and apply it to the likelihood data vector - sck.calculate_concealing_factor() - blinded_dv = sck.apply_concealing_to_likelihood_datavec() - - # Save the blinded data vector to a temporary file in tmp_path - temp_file_path = str(tmp_path) - temp_file_root = "temp_sacc" - temp_file_name = f"{temp_file_path}/{temp_file_root}_concealed_data_vector.fits" - returned_sacc = sck.save_concealed_datavector(temp_file_path, - temp_file_root, - return_sacc=True) - # checks if the return is a sacc object - assert isinstance(returned_sacc, sacc.Sacc) - - returned_sacc = sck.save_concealed_datavector(temp_file_path, - temp_file_root, - return_sacc=False) - - # Check that the return is None - assert returned_sacc is None - - # Check that the file was created - assert os.path.exists(temp_file_name) - - # Load the file and check that the data vector matches the blinded data vector - loaded_sacc = sacc.Sacc.load_fits(temp_file_name) - np.testing.assert_array_equal(loaded_sacc.mean, blinded_dv) - - info_str = 'Concealed (blinded) data-vector, created by Smokescreen.' - assert loaded_sacc.metadata['concealed'] is True - assert loaded_sacc.metadata['creator'] == mock_getuser.return_value - assert loaded_sacc.metadata['creation'][:10] == datetime.date.today().isoformat() - assert loaded_sacc.metadata['info'] == info_str - assert loaded_sacc.metadata['seed_smokescreen'] == 1234 - # File cleanup handled by tmp_path context manager - - -@patch('src.smokescreen.datavector.getpass.getuser', return_value='test_user') -def test_save_concealed_datavector_hdf5(mock_getuser): - # Create mock inputs using an HDF5 SACC file to ensure input_format is 'hdf5' - cosmo = COSMO - likelihood = "./examples/cosmic_shear/cosmicshear_likelihood.py" - syst_dict = { - "trc1_delta_z": 0.1, - "trc0_delta_z": 0.1, - } - shift_dict = {"Omega_c": 0.34, "sigma8": 0.85} - - # Load from HDF5 source - sacc_data = sacc.Sacc.load_hdf5("./examples/cosmic_shear/cosmicshear_sacc.hdf5") - - # Create ConcealDataVector with input_format='hdf5' - sck = ConcealDataVector(cosmo, likelihood, - shift_dict, sacc_data, syst_dict, seed=1234, - **{'input_format': 'hdf5'}) - - # Calculate the concealing factor and apply it to the likelihood data vector - sck.calculate_concealing_factor() - blinded_dv = sck.apply_concealing_to_likelihood_datavec() - - # Save the blinded data vector to a temporary HDF5 file - temp_file_path = "./tests/" - temp_file_root = "temp_sacc_hdf5" - temp_file_name = f"{temp_file_path}{temp_file_root}_concealed_data_vector.hdf5" - - # Save with output_format='hdf5' to test HDF5 output - returned_sacc = sck.save_concealed_datavector(temp_file_path, - temp_file_root, - return_sacc=True, - output_format='hdf5') - - # Check that the return is a sacc object - assert isinstance(returned_sacc, sacc.Sacc) - - # Check that the file was created with .hdf5 extension - assert os.path.exists(temp_file_name) - # Verify it's an HDF5 file by checking extension in path - assert temp_file_name.endswith('.hdf5') - - # Load the HDF5 file and check that the data vector matches - loaded_sacc = sacc.Sacc.load_hdf5(temp_file_name) - np.testing.assert_array_equal(loaded_sacc.mean, blinded_dv) - - # Check metadata - info_str = 'Concealed (blinded) data-vector, created by Smokescreen.' - assert loaded_sacc.metadata['concealed'] is True - assert loaded_sacc.metadata['creator'] == mock_getuser.return_value - assert loaded_sacc.metadata['info'] == info_str - - # Clean up the temporary file - os.remove(temp_file_name) - - -@patch('src.smokescreen.datavector.getpass.getuser', return_value='test_user') -def test_save_concealed_datavector_hdf5_from_fits_input(mock_getuser): - # Test that when input format is FITS but output format is explicitly set to HDF5, - # the file is saved with .hdf5 extension - cosmo = COSMO - likelihood = "./examples/cosmic_shear/cosmicshear_likelihood.py" - syst_dict = { - "trc1_delta_z": 0.1, - "trc0_delta_z": 0.1, - } - shift_dict = {"Omega_c": 0.34, "sigma8": 0.85} - - # Load from FITS source - sacc_data = sacc.Sacc.load_fits("./examples/cosmic_shear/cosmicshear_sacc.fits") - - # Create ConcealDataVector with input_format='fits' - sck = ConcealDataVector(cosmo, likelihood, - shift_dict, sacc_data, syst_dict, seed=1234, - **{'input_format': 'fits'}) - - # Calculate the concealing factor and apply it - sck.calculate_concealing_factor() - blinded_dv = sck.apply_concealing_to_likelihood_datavec() - - # Save with explicit HDF5 output format - temp_file_path = "./tests/" - temp_file_root = "temp_sacc_hdf5_from_fits" - temp_file_name = f"{temp_file_path}{temp_file_root}_concealed_data_vector.hdf5" - - returned_sacc = sck.save_concealed_datavector(temp_file_path, - temp_file_root, - return_sacc=True, - output_format='hdf5') - - assert isinstance(returned_sacc, sacc.Sacc) - assert os.path.exists(temp_file_name) - - # Verify it can be loaded as HDF5 - loaded_sacc = sacc.Sacc.load_hdf5(temp_file_name) - np.testing.assert_array_equal(loaded_sacc.mean, blinded_dv) - - os.remove(temp_file_name) - - -@patch('src.smokescreen.datavector.getpass.getuser', return_value='test_user') -def test_save_concealed_datavector_default_format_uses_input_format(mock_getuser): - # Test that when output_format is not specified, it defaults to input format - cosmo = COSMO - likelihood = "./examples/cosmic_shear/cosmicshear_likelihood.py" - syst_dict = { - "trc1_delta_z": 0.1, - "trc0_delta_z": 0.1, - } - shift_dict = {"Omega_c": 0.34, "sigma8": 0.85} - - # Test with HDF5 input - output should also be HDF5 (.hdf5) - sacc_data_hdf5, _ = load_sacc_file("./examples/cosmic_shear/cosmicshear_sacc.hdf5") - sck = ConcealDataVector(cosmo, likelihood, - shift_dict, sacc_data_hdf5, syst_dict, seed=1234) - - sck.calculate_concealing_factor() - blinded_dv = sck.apply_concealing_to_likelihood_datavec() - - temp_file_path = "./tests/" - temp_file_root = "temp_sacc_default_hdf5" - - # Don't specify output_format - should use input format (hdf5) - returned_sacc = sck.save_concealed_datavector(temp_file_path, - temp_file_root, - return_sacc=True) - - assert isinstance(returned_sacc, sacc.Sacc) - - # Check that the file has .hdf5 extension (from HDF5 input format) - expected_hdf5_name = f"{temp_file_path}{temp_file_root}_concealed_data_vector.hdf5" - assert os.path.exists(expected_hdf5_name) - - # Verify it can be loaded as HDF5 - loaded_sacc = sacc.Sacc.load_hdf5(expected_hdf5_name) - np.testing.assert_array_equal(loaded_sacc.mean, blinded_dv) - - os.remove(expected_hdf5_name) - - -@patch('src.smokescreen.datavector.getpass.getuser', return_value='test_user') -def test_save_concealed_datavector_custom_suffix(mock_getuser): - cosmo = COSMO - likelihood = "./examples/cosmic_shear/cosmicshear_likelihood.py" - syst_dict = { - "trc1_delta_z": 0.1, - "trc0_delta_z": 0.1, - } - shift_dict = {"Omega_c": 0.34, "sigma8": 0.85} - sacc_data, _ = load_sacc_file("./examples/cosmic_shear/cosmicshear_sacc.fits") - sck = ConcealDataVector(cosmo, likelihood, shift_dict, sacc_data, syst_dict, seed=1234) - sck.calculate_concealing_factor() - sck.apply_concealing_to_likelihood_datavec() - - temp_file_path = "./tests/" - temp_file_root = "temp_sacc_custom_suffix" - - sck.save_concealed_datavector(temp_file_path, temp_file_root, suffix="my_blind") - - expected_path = f"{temp_file_path}{temp_file_root}_my_blind.fits" - default_path = f"{temp_file_path}{temp_file_root}_concealed_data_vector.fits" - - assert os.path.exists(expected_path) - assert not os.path.exists(default_path) - - os.remove(expected_path) - - -@patch('src.smokescreen.datavector.getpass.getuser', return_value='test_user') -def test_save_concealed_datavector_default_suffix(mock_getuser): - cosmo = COSMO - likelihood = "./examples/cosmic_shear/cosmicshear_likelihood.py" - syst_dict = { - "trc1_delta_z": 0.1, - "trc0_delta_z": 0.1, - } - shift_dict = {"Omega_c": 0.34, "sigma8": 0.85} - sacc_data, _ = load_sacc_file("./examples/cosmic_shear/cosmicshear_sacc.fits") - sck = ConcealDataVector(cosmo, likelihood, shift_dict, sacc_data, syst_dict, seed=1234) - sck.calculate_concealing_factor() - sck.apply_concealing_to_likelihood_datavec() - - temp_file_path = "./tests/" - temp_file_root = "temp_sacc_default_suffix" - - sck.save_concealed_datavector(temp_file_path, temp_file_root) - - expected_path = f"{temp_file_path}{temp_file_root}_concealed_data_vector.fits" - assert os.path.exists(expected_path) - - os.remove(expected_path) - - -def test_flat_distribution_and_deterministic_blinding_equivalence(cosmic_shear_resources): - """ - Test that blinding from a flat distribution with a given seed produces the same - data vector as blinding deterministically with the registered sampled shifts. - """ - cosmo = COSMO - likelihood = cosmic_shear_resources['likelihood'] - syst_dict = { - "trc1_delta_z": 0.1, - "trc0_delta_z": 0.1, - } - seed = 1234 - # Distribution-based shifts: sample Omega_c and sigma8 from uniform ranges - distribution_shifts_dict = {"Omega_c": (0.20, 0.39), "sigma8": (0.6, 0.9)} - - # Step 1: Create smokescreen with distribution-based shifts and get blinded data vector - sacc_data = sacc.Sacc.load_fits(cosmic_shear_resources['fits_sacc']) - sck_distribution = ConcealDataVector(cosmo, likelihood, - distribution_shifts_dict, sacc_data, - syst_dict, seed=seed) - sck_distribution.calculate_concealing_factor() - blinded_dv_distribution = sck_distribution.apply_concealing_to_likelihood_datavec() + ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, theory_fn=wrong) - # Step 2: Register the sampled shift values by re-running the draw with the same seed - sampled_shifts = sck_distribution._load_shifts(seed=seed) - # Step 3: Create smokescreen with deterministic shifts equal to the registered sample values - sacc_data2 = sacc.Sacc.load_fits(cosmic_shear_resources['fits_sacc']) - sck_deterministic = ConcealDataVector(cosmo, likelihood, - sampled_shifts, sacc_data2, - syst_dict, seed=seed) - sck_deterministic.calculate_concealing_factor() - blinded_dv_deterministic = sck_deterministic.apply_concealing_to_likelihood_datavec() +# --- concealing factor (add) bit-for-bit ------------------------------------ - # Step 4: Both strategies must produce the same blinded data vector - np.testing.assert_array_almost_equal(blinded_dv_distribution, blinded_dv_deterministic) +def test_concealing_factor_add_bit_for_bit(sacc_data, theory_fn): + # debug mode returns the factor; compare bit-for-bit to an independent recompute + cdv_dbg = ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, + theory_fn=theory_fn, debug=True) + factor = cdv_dbg.calculate_concealing_factor(factor_type="add") + deltas = draw_param_shifts(SHIFTS, SEED, shift_distr="flat") + concealed = {k: FIDUCIAL[k] + deltas.get(k, 0.0) for k in FIDUCIAL} + expected = theory_fn(concealed) - theory_fn(FIDUCIAL) + np.testing.assert_array_equal(factor, expected) -def test_gaussian_distribution_and_deterministic_blinding_equivalence(cosmic_shear_resources): - """ - Test that blinding from a Gaussian distribution with a given seed produces the same - data vector as blinding deterministically with the registered sampled shifts. - """ - cosmo = COSMO - likelihood = cosmic_shear_resources['likelihood'] - syst_dict = { - "trc1_delta_z": 0.1, - "trc0_delta_z": 0.1, - } - seed = 1234 - # Gaussian distribution shifts: (mean, std) tuples - distribution_shifts_dict = {"Omega_c": (0.3, 0.02), "sigma8": (0.82, 0.02)} - # Step 1: Create smokescreen with Gaussian distribution shifts and get blinded data vector - sacc_data = sacc.Sacc.load_fits(cosmic_shear_resources['fits_sacc']) - sck_distribution = ConcealDataVector(cosmo, likelihood, - distribution_shifts_dict, sacc_data, - syst_dict, seed=seed, - **{'shift_distr': 'gaussian'}) - sck_distribution.calculate_concealing_factor() - blinded_dv_distribution = sck_distribution.apply_concealing_to_likelihood_datavec() +def test_apply_add_returns_mean_plus_factor(sacc_data, theory_fn): + cdv = ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, + theory_fn=theory_fn, debug=True) + factor = cdv.calculate_concealing_factor(factor_type="add") + blinded = cdv.apply_concealing_to_likelihood_datavec() + np.testing.assert_array_equal(blinded, sacc_data.mean + factor) - # Step 2: Register the sampled shift values by re-running the Gaussian draw with the same seed - sampled_shifts = sck_distribution._load_shifts(seed=seed, shift_distr='gaussian') - # Step 3: Create smokescreen with deterministic shifts equal to the registered sample values - sacc_data2 = sacc.Sacc.load_fits(cosmic_shear_resources['fits_sacc']) - sck_deterministic = ConcealDataVector(cosmo, likelihood, - sampled_shifts, sacc_data2, - syst_dict, seed=seed) - sck_deterministic.calculate_concealing_factor() - blinded_dv_deterministic = sck_deterministic.apply_concealing_to_likelihood_datavec() +# --- concealing factor (mult) ----------------------------------------------- - # Step 4: Both strategies must produce the same blinded data vector - np.testing.assert_array_almost_equal(blinded_dv_distribution, blinded_dv_deterministic) +def test_concealing_factor_mult(sacc_data, theory_fn): + cdv = ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, + theory_fn=theory_fn, debug=True) + factor = cdv.calculate_concealing_factor(factor_type="mult") + + deltas = draw_param_shifts(SHIFTS, SEED, shift_distr="flat") + concealed = {k: FIDUCIAL[k] + deltas.get(k, 0.0) for k in FIDUCIAL} + expected = theory_fn(concealed) / theory_fn(FIDUCIAL) + np.testing.assert_array_equal(factor, expected) + + blinded = cdv.apply_concealing_to_likelihood_datavec() + np.testing.assert_array_equal(blinded, sacc_data.mean * factor) + + +# --- save path -------------------------------------------------------------- + +def _blinded_cdv(sacc_data, theory_fn): + cdv = ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, + theory_fn=theory_fn) + cdv.calculate_concealing_factor(factor_type="add") + cdv.apply_concealing_to_likelihood_datavec() + return cdv + + +def test_save_fits_and_cov_unchanged(tmp_path, sacc_data, theory_fn, monkeypatch): + monkeypatch.setattr("smokescreen.datavector.getpass.getuser", + lambda: "test_user") + cov_before = sacc_data.covariance.dense.copy() + cdv = _blinded_cdv(sacc_data, theory_fn) + + out = cdv.save_concealed_datavector(str(tmp_path), "root", return_sacc=True) + assert isinstance(out, sacc.Sacc) + np.testing.assert_array_equal(out.mean, cdv.concealed_data_vector) + np.testing.assert_allclose(out.covariance.dense, cov_before) + # input SACC covariance untouched too + np.testing.assert_allclose(sacc_data.covariance.dense, cov_before) + + md = out.metadata + assert md["concealed"] is True + assert md["creator"] == "test_user" + assert "info" in md + assert md["seed_smokescreen"] == SEED + + assert (tmp_path / "root_concealed_data_vector.fits").exists() + + +def test_save_return_sacc_false(tmp_path, sacc_data, theory_fn, monkeypatch): + monkeypatch.setattr("smokescreen.datavector.getpass.getuser", + lambda: "test_user") + cdv = _blinded_cdv(sacc_data, theory_fn) + out = cdv.save_concealed_datavector(str(tmp_path), "root", return_sacc=False) + assert out is None + + +def test_save_custom_suffix(tmp_path, sacc_data, theory_fn, monkeypatch): + monkeypatch.setattr("smokescreen.datavector.getpass.getuser", + lambda: "test_user") + cdv = _blinded_cdv(sacc_data, theory_fn) + cdv.save_concealed_datavector(str(tmp_path), "root", suffix="mysuffix") + assert (tmp_path / "root_mysuffix.fits").exists() + + +def test_save_hdf5_format(tmp_path, sacc_data, theory_fn, monkeypatch): + pytest.importorskip("h5py") # sacc.save_hdf5 needs h5py + monkeypatch.setattr("smokescreen.datavector.getpass.getuser", + lambda: "test_user") + cdv = _blinded_cdv(sacc_data, theory_fn) + cdv.save_concealed_datavector(str(tmp_path), "root", output_format="hdf5") + assert (tmp_path / "root_concealed_data_vector.hdf5").exists() + + +# --- end-to-end synthetic (default backend never constructed) --------------- + +def test_end_to_end_synthetic(tmp_path, sacc_data, theory_fn, monkeypatch): + monkeypatch.setattr("smokescreen.datavector.getpass.getuser", + lambda: "test_user") + cdv = ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, + theory_fn=theory_fn) + cdv.calculate_concealing_factor(factor_type="add") + blinded = cdv.apply_concealing_to_likelihood_datavec() + out = cdv.save_concealed_datavector(str(tmp_path), "root", return_sacc=True) + assert not np.array_equal(blinded, sacc_data.mean) # actually blinded + # Import discipline (pyccl never entering sys.modules on the synthetic path) + # is covered by the subprocess suite in test_import_discipline.py. + assert isinstance(out, sacc.Sacc) + + +# --- default CCL backend ---------------------------------------------------- + +def _make_cosmic_shear_sacc(): + """Small cosmic-shear SACC with 2 NZ tracers and xi_plus/xi_minus rows.""" + s = sacc.Sacc() + z = np.linspace(0.05, 2.0, 50) + for name, zmean in [("src0", 0.5), ("src1", 1.0)]: + nz = np.exp(-0.5 * ((z - zmean) / 0.2) ** 2) + s.add_tracer("NZ", name, z, nz) + thetas = np.array([5.0, 20.0, 60.0, 120.0]) # arcmin + pairs = [("src0", "src0"), ("src0", "src1"), ("src1", "src1")] + for dt in ("galaxy_shear_xi_plus", "galaxy_shear_xi_minus"): + for p in pairs: + for th in thetas: + s.add_data_point(dt, p, 1e-6, theta=th) + n = len(s.mean) + s.add_covariance(np.eye(n) * 1e-12) + return s + + +def test_default_ccl_backend(tmp_path, monkeypatch): + pytest.importorskip("pyccl") + monkeypatch.setattr("smokescreen.datavector.getpass.getuser", + lambda: "test_user") + s = _make_cosmic_shear_sacc() + cov_before = s.covariance.dense.copy() + fiducial = {"sigma8": 0.81, "Omega_c": 0.26, "Omega_b": 0.045, + "h": 0.67, "n_s": 0.96} + shifts = {"sigma8": 0.05, "Omega_c": (-0.03, 0.03)} + + cdv = ConcealDataVector(fiducial, shifts, s, seed=SEED, theory_fn=None) + assert np.all(np.isfinite(cdv.theory_vec_fid)) + cdv.calculate_concealing_factor(factor_type="add") + blinded = cdv.apply_concealing_to_likelihood_datavec() + assert not np.array_equal(blinded, s.mean) + + out = cdv.save_concealed_datavector(str(tmp_path), "cs", return_sacc=True) + np.testing.assert_allclose(out.covariance.dense, cov_before) diff --git a/tests/test_import_discipline.py b/tests/test_import_discipline.py new file mode 100644 index 0000000..480e562 --- /dev/null +++ b/tests/test_import_discipline.py @@ -0,0 +1,38 @@ +"""Import-discipline guarantees, verified in fresh subprocesses. + +`import smokescreen` and `import smokescreen.datavector` must not pull in pyccl +or firecrown; pyccl enters sys.modules only once the default CCL backend is +constructed. +""" +import subprocess +import sys + + +def _run(code): + res = subprocess.run([sys.executable, "-c", code], + capture_output=True, text=True) + assert res.returncode == 0, f"subprocess failed:\n{res.stdout}\n{res.stderr}" + + +def test_import_smokescreen_no_pyccl_no_firecrown(): + _run( + "import sys, smokescreen\n" + "assert 'pyccl' not in sys.modules, 'pyccl imported by smokescreen'\n" + "assert 'firecrown' not in sys.modules, 'firecrown imported by smokescreen'\n" + ) + + +def test_import_datavector_no_pyccl_no_firecrown(): + _run( + "import sys, smokescreen.datavector\n" + "assert 'pyccl' not in sys.modules, 'pyccl imported by datavector'\n" + "assert 'firecrown' not in sys.modules, 'firecrown imported by datavector'\n" + ) + + +def test_building_ccl_backend_imports_pyccl(): + _run( + "import sys\n" + "from smokescreen.backends.ccl import build_ccl_theory_fn\n" + "assert 'pyccl' in sys.modules, 'pyccl not imported by ccl backend'\n" + ) diff --git a/tests/test_main.py b/tests/test_main.py index 952600a..6b64835 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -2,163 +2,92 @@ import os from unittest.mock import patch, MagicMock from cryptography.fernet import Fernet -from pyccl import CosmologyVanillaLCDM -from smokescreen.utils import load_cosmology_from_partial_dict from smokescreen import __main__ from smokescreen.__main__ import encrypt_main, decrypt_main -@patch('builtins.print') -@patch('smokescreen.__main__.ConcealDataVector') -@patch('smokescreen.__main__.load_sacc_file') -def test_main(mock_load_sacc, mock_smokescreen, mock_print): - # Arrange - path_to_sacc = "./examples/cosmic_shear/cosmicshear_sacc.fits" - likelihood_path = "./tests/test_data/mock_likelihood.py" - systematics = {} - shifts_dict = {"Omega_c": [-0.1, 0.2], "sigma8": [-0.1, 0.1]} - shift_type = 'add' - shift_distribution = 'flat' - seed = 2112 - reference_cosmology = CosmologyVanillaLCDM() - path_to_output = "./tests/test_data/" - keep_original_sacc = True - - mock_smokescreen_instance = MagicMock() - mock_smokescreen.return_value = mock_smokescreen_instance +FIDUCIAL = {"sigma8": 0.8, "Omega_c": 0.25, "Omega_b": 0.05, "h": 0.67, "n_s": 0.96} +SHIFTS = {"Omega_c": (-0.1, 0.2), "sigma8": (-0.1, 0.1)} - sacc_file = MagicMock() # Create a mock sacc_file - mock_load_sacc.return_value = (sacc_file, 'fits') # Make load_sacc return tuple of (sacc, format) - # Act - __main__.datavector_main(path_to_sacc, likelihood_path, shifts_dict, systematics, - shift_type, shift_distribution, seed, reference_cosmology, - path_to_output, keep_original_sacc) - - # Assert - mock_load_sacc.assert_called_once_with(path_to_sacc) - mock_smokescreen.assert_called_once_with(reference_cosmology, likelihood_path, - shifts_dict, sacc_file, systematics, seed, - shift_distr=shift_distribution, input_format='fits') - mock_smokescreen_instance.calculate_concealing_factor.assert_called_once() - mock_smokescreen_instance.apply_concealing_to_likelihood_datavec.assert_called_once() - mock_smokescreen_instance.save_concealed_datavector.assert_called_once_with( - path_to_output, 'cosmicshear_sacc', output_format='fits', suffix=None - ) - - -@patch('builtins.print') +@patch('smokescreen.__main__.encrypt_file') @patch('smokescreen.__main__.ConcealDataVector') @patch('smokescreen.__main__.load_sacc_file') -def test_main_loads_cosmology_from_dict(mock_load_sacc, mock_smokescreen, mock_print): - # Arrange +def test_datavector_main(mock_load_sacc, mock_smokescreen, mock_encrypt, mock_print=None): path_to_sacc = "./examples/cosmic_shear/cosmicshear_sacc.fits" - likelihood_path = "./tests/test_data/mock_likelihood.py" - systematics = {} - shifts_dict = {"Omega_c": [-0.1, 0.2], "sigma8": [-0.1, 0.1]} + seed = 2112 shift_type = 'add' shift_distribution = 'flat' - seed = 2112 - reference_cosmology = {'sigma8': 0.888} path_to_output = "./tests/test_data/" - keep_original_sacc = True - - mock_smokescreen_instance = MagicMock() - mock_smokescreen.return_value = mock_smokescreen_instance - sacc_file = MagicMock() # Create a mock sacc_file - mock_load_sacc.return_value = (sacc_file, 'fits') # Make load_sacc return tuple of (sacc, format) + mock_smoke_inst = MagicMock() + mock_smokescreen.return_value = mock_smoke_inst + sacc_file = MagicMock() + mock_load_sacc.return_value = (sacc_file, 'fits') + mock_encrypt.return_value = (b'enc', b'key') - # Act - __main__.datavector_main(path_to_sacc, likelihood_path, shifts_dict, systematics, shift_type, - shift_distribution, seed, reference_cosmology, - path_to_output, keep_original_sacc) + with patch('os.path.exists', return_value=True), patch('builtins.print'): + __main__.datavector_main(path_to_sacc, FIDUCIAL, SHIFTS, seed, + shift_type, shift_distribution, + path_to_output, keep_original_sacc=True) - # Assert - mod_ref_cosmo = load_cosmology_from_partial_dict(reference_cosmology) mock_load_sacc.assert_called_once_with(path_to_sacc) - mock_smokescreen.assert_called_once_with(mod_ref_cosmo, likelihood_path, shifts_dict, - sacc_file, systematics, seed, - shift_distr=shift_distribution, input_format='fits') - mock_smokescreen_instance.calculate_concealing_factor.assert_called_once() - mock_smokescreen_instance.apply_concealing_to_likelihood_datavec.assert_called_once() - mock_smokescreen_instance.save_concealed_datavector.assert_called_once_with( + mock_smokescreen.assert_called_once_with( + FIDUCIAL, SHIFTS, sacc_file, seed=seed, + shift_distr=shift_distribution, input_format='fits') + mock_smoke_inst.calculate_concealing_factor.assert_called_once() + mock_smoke_inst.apply_concealing_to_likelihood_datavec.assert_called_once() + mock_smoke_inst.save_concealed_datavector.assert_called_once_with( path_to_output, 'cosmicshear_sacc', output_format='fits', suffix=None ) -@patch('builtins.print') +@patch('smokescreen.__main__.encrypt_file') @patch('smokescreen.__main__.ConcealDataVector') @patch('smokescreen.__main__.load_sacc_file') -def test_main_gaussian_shift(mock_load_sacc, mock_smokescreen, mock_print): - # Arrange +def test_datavector_main_gaussian_shift(mock_load_sacc, mock_smokescreen, mock_encrypt): path_to_sacc = "./examples/cosmic_shear/cosmicshear_sacc.fits" - likelihood_path = "./tests/test_data/mock_likelihood.py" - systematics = {} - shifts_dict = {"Omega_c": [-0.1, 0.2], "sigma8": [-0.1, 0.1]} - shift_type = 'add' - shift_distribution = 'gaussian' seed = 2112 - reference_cosmology = CosmologyVanillaLCDM() path_to_output = "./tests/test_data/" - keep_original_sacc = True - mock_smokescreen_instance = MagicMock() - mock_smokescreen.return_value = mock_smokescreen_instance - - sacc_file = MagicMock() # Create a mock sacc_file - mock_load_sacc.return_value = (sacc_file, 'fits') # Make load_sacc return tuple of (sacc, format) + mock_smoke_inst = MagicMock() + mock_smokescreen.return_value = mock_smoke_inst + sacc_file = MagicMock() + mock_load_sacc.return_value = (sacc_file, 'fits') + mock_encrypt.return_value = (b'enc', b'key') - # Act - __main__.datavector_main(path_to_sacc, likelihood_path, shifts_dict, - systematics, shift_type, shift_distribution, seed, - reference_cosmology, - path_to_output, keep_original_sacc) + with patch('os.path.exists', return_value=True), patch('builtins.print'): + __main__.datavector_main(path_to_sacc, FIDUCIAL, SHIFTS, seed, + shift_type='add', shift_distribution='gaussian', + path_to_output=path_to_output, keep_original_sacc=True) - # Assert - mock_load_sacc.assert_called_once_with(path_to_sacc) - mock_smokescreen.assert_called_once_with(reference_cosmology, - likelihood_path, shifts_dict, sacc_file, - systematics, seed, - shift_distr=shift_distribution, input_format='fits') - mock_smokescreen_instance.calculate_concealing_factor.assert_called_once() - mock_smokescreen_instance.apply_concealing_to_likelihood_datavec.assert_called_once() - mock_smokescreen_instance.save_concealed_datavector.assert_called_once_with( - path_to_output, 'cosmicshear_sacc', output_format='fits', suffix=None - ) + mock_smokescreen.assert_called_once_with( + FIDUCIAL, SHIFTS, sacc_file, seed=seed, + shift_distr='gaussian', input_format='fits') -@patch('builtins.print') +@patch('smokescreen.__main__.encrypt_file') @patch('smokescreen.__main__.ConcealDataVector') @patch('smokescreen.__main__.load_sacc_file') -def test_datavector_main_custom_suffix(mock_load_sacc, mock_smokescreen, mock_print): - # Arrange +def test_datavector_main_custom_suffix(mock_load_sacc, mock_smokescreen, mock_encrypt): path_to_sacc = "./examples/cosmic_shear/cosmicshear_sacc.fits" - likelihood_path = "./tests/test_data/mock_likelihood.py" - systematics = {} - shifts_dict = {"Omega_c": [-0.1, 0.2], "sigma8": [-0.1, 0.1]} - shift_type = 'add' - shift_distribution = 'flat' seed = 2112 - reference_cosmology = CosmologyVanillaLCDM() path_to_output = "./tests/test_data/" - keep_original_sacc = True - output_suffix = "my_suffix" - - mock_smokescreen_instance = MagicMock() - mock_smokescreen.return_value = mock_smokescreen_instance + mock_smoke_inst = MagicMock() + mock_smokescreen.return_value = mock_smoke_inst sacc_file = MagicMock() mock_load_sacc.return_value = (sacc_file, 'fits') + mock_encrypt.return_value = (b'enc', b'key') - # Act - __main__.datavector_main(path_to_sacc, likelihood_path, shifts_dict, systematics, - shift_type, shift_distribution, seed, reference_cosmology, - path_to_output, keep_original_sacc, - output_suffix=output_suffix) + with patch('os.path.exists', return_value=True), patch('builtins.print'): + __main__.datavector_main(path_to_sacc, FIDUCIAL, SHIFTS, seed, + shift_type='add', shift_distribution='flat', + path_to_output=path_to_output, + keep_original_sacc=True, + output_suffix="my_suffix") - # Assert - mock_smokescreen_instance.save_concealed_datavector.assert_called_once_with( + mock_smoke_inst.save_concealed_datavector.assert_called_once_with( path_to_output, 'cosmicshear_sacc', output_format='fits', suffix="my_suffix" ) diff --git a/tests/test_param_shifts.py b/tests/test_param_shifts.py index 126a74c..10fb96c 100644 --- a/tests/test_param_shifts.py +++ b/tests/test_param_shifts.py @@ -1,144 +1,97 @@ -import pytest # noqa F401 -import pyccl as ccl -from smokescreen.param_shifts import draw_flat_param_shifts -from smokescreen.param_shifts import draw_flat_or_deterministic_param_shifts -from smokescreen.param_shifts import draw_gaussian_param_shifts +import numpy as np +import pytest +from smokescreen.param_shifts import ( + draw_param_shifts, + draw_flat_param_shifts, + _normalize_seed, +) -# tests for draw_flat_param_shifts -def test_single_value_shifts(): - shift_dict = {'param1': 1, 'param2': 2} - seed = 123 - result = draw_flat_param_shifts(shift_dict, seed) - assert isinstance(result, dict) - assert set(result.keys()) == set(shift_dict.keys()) - for key, value in result.items(): - assert -shift_dict[key] <= value <= shift_dict[key] +# --- reproducibility & structure -------------------------------------------- -def test_tuple_value_shifts(): - shift_dict = {'param1': (1, 2), 'param2': (2, 3)} - seed = 123 - result = draw_flat_param_shifts(shift_dict, seed) - assert isinstance(result, dict) - assert set(result.keys()) == set(shift_dict.keys()) - for key, value in result.items(): - assert -shift_dict[key][0] <= value <= shift_dict[key][1] +def test_fixed_seed_reproduces_fixed_shift(): + shifts = {"Omega_c": 0.05, "sigma8": (-0.1, 0.1)} + a = draw_param_shifts(shifts, 2112) + b = draw_param_shifts(shifts, 2112) + assert a == b + assert set(a) == set(shifts) -def test_string_seed(): - shift_dict = {'param1': 1, 'param2': 2} - seed = 'random_seed' - result = draw_flat_param_shifts(shift_dict, seed) - assert isinstance(result, dict) - assert set(result.keys()) == set(shift_dict.keys()) - for key, value in result.items(): - assert -shift_dict[key] <= value <= shift_dict[key] +def test_key_order_independence(): + forward = {"Omega_c": 0.05, "sigma8": 0.1, "h": (-0.02, 0.03)} + reversed_dict = {"h": (-0.02, 0.03), "sigma8": 0.1, "Omega_c": 0.05} + a = draw_param_shifts(forward, 2112) + b = draw_param_shifts(reversed_dict, 2112) + assert a == b # per-key deltas identical regardless of insertion order -# tests for draw_flat_or_deterministic_param_shifts -def test_draw_flat_or_deterministic_param_shifts_deterministic(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"Omega_c": 0.01, "sigma8": 0.02} - seed = 1234 +def test_draw_does_not_perturb_global_np_random(): + shifts = {"Omega_c": 0.05, "sigma8": 0.1} + np.random.seed(0) + a = np.random.random() + np.random.seed(0) + draw_param_shifts(shifts, 2112) + b = np.random.random() + assert a == b # local RNG never touched the global state - shifts = draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed) - assert shifts["Omega_c"] == 0.01 - assert shifts["sigma8"] == 0.02 +# --- delta semantics -------------------------------------------------------- +def test_flat_float_is_delta_within_symmetric_envelope(): + h = 0.05 + shifts = {"Omega_c": h} + for seed in range(20): + delta = draw_param_shifts(shifts, seed)["Omega_c"] + assert -h <= delta <= h -def test_draw_flat_or_deterministic_param_shifts_flat(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"Omega_c": (0.01, 0.02), "sigma8": (0.02, 0.03)} - seed = 1234 - shifts = draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed) +def test_flat_tuple_is_delta_within_box(): + lo, hi = -0.03, 0.07 + shifts = {"Omega_c": (lo, hi)} + for seed in range(20): + delta = draw_param_shifts(shifts, seed)["Omega_c"] + assert lo <= delta <= hi - assert 0.01 <= shifts["Omega_c"] <= 0.02 - assert 0.02 <= shifts["sigma8"] <= 0.03 +# --- seed normalization ----------------------------------------------------- -def test_draw_flat_or_deterministic_param_shifts_invalid_key(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"invalid_param": 0.01} - seed = 1234 +def test_string_seed_matches_normalized_int(): + shifts = {"Omega_c": 0.05, "sigma8": (-0.1, 0.1)} + str_seed = "my_secret_seed" + int_seed = _normalize_seed(str_seed) + assert draw_param_shifts(shifts, str_seed) == draw_param_shifts(shifts, int_seed) - with pytest.raises(ValueError, match=r"Key invalid_param not in cosmology parameters"): - draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed) +# --- distributions ---------------------------------------------------------- -def test_draw_flat_or_deterministic_param_shifts_invalid_tuple_length(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"Omega_c": (0.01, 0.02, 0.03)} - seed = 1234 +def test_gaussian_float(): + shifts = {"Omega_c": 0.05} + result = draw_param_shifts(shifts, 2112, shift_distr="gaussian") + assert isinstance(result["Omega_c"], float) - with pytest.raises(ValueError, match=r"Tuple \(0.01, 0.02, 0.03\) has to be of length 2"): - draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed) +def test_gaussian_tuple_raises(): + shifts = {"Omega_c": (-0.05, 0.05)} + with pytest.raises(ValueError): + draw_param_shifts(shifts, 2112, shift_distr="gaussian") -def test_draw_flat_or_deterministic_param_shifts_string_seed(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"Omega_c": (0.01, 0.02), "sigma8": (0.02, 0.03)} - seed = "1234" - shifts = draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed) +def test_flat_tuple_wrong_length_raises(): + shifts = {"Omega_c": (0.0, 0.05, 0.1)} + with pytest.raises(ValueError): + draw_param_shifts(shifts, 2112) - assert 0.01 <= shifts["Omega_c"] <= 0.02 - assert 0.02 <= shifts["sigma8"] <= 0.03 +def test_unknown_distribution_raises(): + with pytest.raises(NotImplementedError): + draw_param_shifts({"Omega_c": 0.05}, 2112, shift_distr="lorentzian") -# tests for draw_gaussian_param_shifts -def test_draw_gaussian_param_shifts_single_value(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"Omega_c": (0.3, 0.01), "sigma8": (0.8, 0.02)} - seed = 1234 - shifts = draw_gaussian_param_shifts(cosmo, shifts_dict, seed) +# --- back-compat alias ------------------------------------------------------ - assert isinstance(shifts, dict) - assert set(shifts.keys()) == set(shifts_dict.keys()) - for key, value in shifts.items(): - mean, std = shifts_dict[key] - assert mean - 3 * std <= value <= mean + 3 * std - - -def test_draw_gaussian_param_shifts_invalid_key(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"invalid_param": (0.3, 0.01)} - seed = 1234 - - with pytest.raises(ValueError, match=r"Key invalid_param not in cosmology parameters"): - draw_gaussian_param_shifts(cosmo, shifts_dict, seed) - - -def test_draw_gaussian_param_shifts_invalid_tuple_length(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"Omega_c": (0.3, 0.01, 0.02)} - seed = 1234 - - with pytest.raises(ValueError, match=r"Tuple \(0.3, 0.01, 0.02\) has to be of length 2"): - draw_gaussian_param_shifts(cosmo, shifts_dict, seed) - - -def test_draw_gaussian_param_shifts_non_tuple_value(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"Omega_c": 0.3} - seed = 1234 - - with pytest.raises(ValueError, match=r"Value 0.3 has to be a tuple of length 2"): - draw_gaussian_param_shifts(cosmo, shifts_dict, seed) - - -def test_draw_gaussian_param_shifts_string_seed(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"Omega_c": (0.3, 0.01), "sigma8": (0.8, 0.02)} - seed = "1234" - - shifts = draw_gaussian_param_shifts(cosmo, shifts_dict, seed) - - assert isinstance(shifts, dict) - assert set(shifts.keys()) == set(shifts_dict.keys()) - for key, value in shifts.items(): - mean, std = shifts_dict[key] - assert mean - 3 * std <= value <= mean + 3 * std +def test_draw_flat_param_shifts_alias(): + shifts = {"Omega_c": 0.05, "sigma8": (-0.1, 0.1)} + assert draw_flat_param_shifts(shifts, 2112) == draw_param_shifts( + shifts, 2112, shift_distr="flat" + ) From 8b80bf2dd88d6b45aad7a651b026fd5e200a32cb Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 10 Jul 2026 23:52:42 +0200 Subject: [PATCH 4/7] Document the default backend's transfer-function and theta-arcmin conventions Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/smokescreen/backends/ccl.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/smokescreen/backends/ccl.py b/src/smokescreen/backends/ccl.py index 079a068..441c512 100644 --- a/src/smokescreen/backends/ccl.py +++ b/src/smokescreen/backends/ccl.py @@ -17,6 +17,20 @@ parameter mapping, the theory vector aligned element-for-element to ``sacc_data.mean``. +Two conventions of this convenience backend, part of its row-to-theory +contract (Smokescreen does not verify them): + +- **Transfer function / power spectrum.** It defaults to CCL's native + ``transfer_function="eisenstein_hu"`` and ``matter_power_spectrum="halofit"`` + so the shipped backend runs against a bare ``pyccl`` install with no + Boltzmann code (CAMB/CLASS). Both are overridable through ``cosmo_params``; + a caller who wants a Boltzmann-code power spectrum can supply their own + ``theory_fn`` instead. +- **ξ± angle units.** The SACC ``theta`` tag of ``galaxy_shear_xi_plus`` / + ``galaxy_shear_xi_minus`` rows is interpreted as **arcminutes** and converted + to degrees for :func:`pyccl.correlation`. ``galaxy_shear_cl_ee`` rows use + their ``ell`` tag directly. + .. autofunction:: build_ccl_theory_fn ''' import numpy as np From 59609ad7da236212ed7727513c251e25196ffc5e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 00:05:25 +0200 Subject: [PATCH 5/7] Make param-shift draws per-key independent; pin golden values; accept list envelopes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Each key now gets its own RNG derived from (seed, key), so the docstring contract — a delta depends only on (key, seed, shift_distr) — is actually true: adding or removing parameters from the shifts envelope no longer changes other parameters' hidden shifts. Golden-value tests pin the seed-normalization and both distributions cross-run/cross-machine; a subset-invariance test guards the per-key contract; gaussian tests now discriminate the distribution (differs from flat, escapes the flat bound). [lo, hi] lists (the natural YAML/JSON form) are accepted alongside tuples. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/smokescreen/param_shifts.py | 35 +++++++++++++++++----------- tests/test_param_shifts.py | 41 +++++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 14 deletions(-) diff --git a/src/smokescreen/param_shifts.py b/src/smokescreen/param_shifts.py index baf648f..fd62375 100644 --- a/src/smokescreen/param_shifts.py +++ b/src/smokescreen/param_shifts.py @@ -55,18 +55,21 @@ def draw_param_shifts(shifts_dict, seed, shift_distr="flat"): Draw a delta for each named parameter from a local RNG seeded by ``seed``. Each returned value is a shift to be *added* to the fiducial value, not an - absolute parameter value. The draw is order-independent: the delta for a - given key depends only on ``(key, seed, shift_distr)``, not on the - iteration order of ``shifts_dict``. This is achieved by iterating the keys - in sorted order and consuming exactly one scalar RNG call per key. + absolute parameter value. The draw is per-key independent: the delta for a + given key depends only on ``(key, seed, shift_distr)`` --- not on the + iteration order of ``shifts_dict``, and not on which *other* keys are + present. Each key gets its own RNG derived from ``(seed, key)``, so adding + or removing a parameter from the shifts envelope never changes any other + parameter's drawn delta. Parameters ---------- shifts_dict : Mapping[str, float | tuple[float, float]] Mapping of parameter name to its shift envelope, interpreted as a delta about zero. A ``float`` ``h`` is the symmetric delta envelope - ``(-h, +h)``; a ``(lo, hi)`` tuple is a delta box, the drawn delta - lies in ``[lo, hi]`` (typically straddling zero). + ``(-h, +h)``; a ``(lo, hi)`` pair (tuple or list, e.g. as loaded from + YAML/JSON) is a delta box, the drawn delta lies in ``[lo, hi]`` + (typically straddling zero). seed : int or str Seed for the random number generator. A ``str`` is normalized to an ``int`` (see :func:`_normalize_seed`). @@ -85,25 +88,29 @@ def draw_param_shifts(shifts_dict, seed, shift_distr="flat"): if shift_distr not in ("flat", "gaussian"): raise NotImplementedError("Only flat and gaussian shifts are implemented") - rng = np.random.default_rng(_normalize_seed(seed)) + base_seed = _normalize_seed(seed) shifts = {} - # Sorted iteration makes the per-key draw independent of dict order. - for key in sorted(shifts_dict): - envelope = shifts_dict[key] + for key, envelope in shifts_dict.items(): + # Per-key derived RNG: the draw for a key depends only on (seed, key), + # never on the other keys in the mapping or on iteration order. + rng = np.random.default_rng([base_seed, _normalize_seed(key)]) if shift_distr == "flat": - if isinstance(envelope, tuple): + if isinstance(envelope, (tuple, list)): if len(envelope) != 2: - raise ValueError(f"Tuple {envelope} has to be of length 2") + raise ValueError( + f"Shift envelope for {key!r} must have length 2 " + f"(lo, hi); got {envelope!r}" + ) lo, hi = envelope else: lo, hi = -envelope, envelope shifts[key] = float(rng.uniform(lo, hi)) else: # gaussian - if isinstance(envelope, tuple): + if isinstance(envelope, (tuple, list)): raise ValueError( f"Gaussian shift for {key} requires a float half-width, " - f"got tuple {envelope}" + f"got {envelope!r}" ) shifts[key] = float(rng.normal(0.0, envelope)) return shifts diff --git a/tests/test_param_shifts.py b/tests/test_param_shifts.py index 10fb96c..a664c7c 100644 --- a/tests/test_param_shifts.py +++ b/tests/test_param_shifts.py @@ -26,6 +26,27 @@ def test_key_order_independence(): assert a == b # per-key deltas identical regardless of insertion order +def test_key_subset_independence(): + # a key's delta depends only on (key, seed, shift_distr) — adding or + # removing OTHER keys must not change it (staged-blinding reproducibility) + full = draw_param_shifts({"sigma8": 0.05, "Omega_c": 0.03, "h": 0.01}, 42) + alone = draw_param_shifts({"sigma8": 0.05}, 42) + assert full["sigma8"] == alone["sigma8"] + + +def test_golden_values(): + # Pin the cross-run/cross-machine reproducibility contract: these exact + # values are how historical blinding shifts are re-derived at unblinding. + # Any RNG/normalization/derivation change MUST fail this test. + assert _normalize_seed("my_secret_seed") == 3734715429001406524 + assert draw_param_shifts({"Omega_c": 0.05}, 2112) == { + "Omega_c": pytest.approx(-0.034665567512199846, abs=0, rel=1e-15) + } + assert draw_param_shifts({"Omega_c": 0.05}, 2112, shift_distr="gaussian") == { + "Omega_c": pytest.approx(0.008622042117174187, abs=0, rel=1e-15) + } + + def test_draw_does_not_perturb_global_np_random(): shifts = {"Omega_c": 0.05, "sigma8": 0.1} np.random.seed(0) @@ -54,6 +75,13 @@ def test_flat_tuple_is_delta_within_box(): assert lo <= delta <= hi +def test_flat_list_envelope_equivalent_to_tuple(): + # [lo, hi] is what YAML/JSON configs naturally produce + as_list = draw_param_shifts({"Omega_c": [-0.03, 0.07]}, 2112) + as_tuple = draw_param_shifts({"Omega_c": (-0.03, 0.07)}, 2112) + assert as_list == as_tuple + + # --- seed normalization ----------------------------------------------------- def test_string_seed_matches_normalized_int(): @@ -69,6 +97,19 @@ def test_gaussian_float(): shifts = {"Omega_c": 0.05} result = draw_param_shifts(shifts, 2112, shift_distr="gaussian") assert isinstance(result["Omega_c"], float) + # gaussian must actually be a different distribution than flat + flat = draw_param_shifts(shifts, 2112, shift_distr="flat") + assert result["Omega_c"] != flat["Omega_c"] + + +def test_gaussian_escapes_flat_envelope(): + # gaussian is unbounded; flat is confined to [-h, h] — a cheap + # discriminator that the gaussian branch is not silently flat + h = 0.05 + assert any( + abs(draw_param_shifts({"Omega_c": h}, s, shift_distr="gaussian")["Omega_c"]) > h + for s in range(200) + ) def test_gaussian_tuple_raises(): From ea864eb0ba0045716d7702497a30554559b2d122 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 00:05:39 +0200 Subject: [PATCH 6/7] Validate shift keys and theory_fn shape at construction; test seed custody and the Cl branch Unknown shifts_dict keys now raise a clear ValueError instead of a bare KeyError; the length guard becomes a full shape check (via a shared _checked_theory_vec, also applied to the concealed-point evaluation), so an (n,1) theory_fn can no longer broadcast the blinded mean into an (n,n) matrix and silently corrupt the saved SACC. New tests: column-vector rejection, unknown-key rejection, seed-is-required TypeError, and cl_ee rows in the default-CCL-backend fixture so the np.interp Cl branch is exercised. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- src/smokescreen/datavector.py | 34 ++++++++++++++++++++++++---------- tests/test_datavector.py | 30 +++++++++++++++++++++++++++++- 2 files changed, 53 insertions(+), 11 deletions(-) diff --git a/src/smokescreen/datavector.py b/src/smokescreen/datavector.py index ae71d5c..ae1c470 100644 --- a/src/smokescreen/datavector.py +++ b/src/smokescreen/datavector.py @@ -84,6 +84,15 @@ def __init__(self, fiducial_params, shifts_dict, sacc_data, *, self.fiducial_params = dict(fiducial_params) self.shifts_dict = shifts_dict + + # every shifted parameter must exist at the fiducial point + unknown = set(shifts_dict) - set(self.fiducial_params) + if unknown: + raise ValueError( + f"shifts_dict keys {sorted(unknown)} are not keys of " + f"fiducial_params {sorted(self.fiducial_params)}; every " + f"shifted parameter must have a fiducial value." + ) self.sacc_data = sacc_data self.seed = seed self.shift_distr = shift_distr @@ -108,19 +117,24 @@ def __init__(self, fiducial_params, shifts_dict, sacc_data, *, shift_distr=shift_distr) self.__concealed_params = self._overlay_shifts(self.__shifts) - # length guard: theory vector must align to the SACC rows - self.theory_vec_fid = np.asarray(self.theory_fn(self.fiducial_params)) - if len(self.theory_vec_fid) != len(self.sacc_data.mean): - raise ValueError( - f"theory_fn returned {len(self.theory_vec_fid)} elements but " - f"sacc_data.mean has {len(self.sacc_data.mean)}; the SACC must " - f"contain exactly the rows theory_fn returns." - ) + # shape guard: theory vector must align to the SACC rows + self.theory_vec_fid = self._checked_theory_vec(self.fiducial_params) if self._debug: print(f"[DEBUG] Shifts: {self.__shifts}") print(f"[DEBUG] Concealed params: {self.__concealed_params}") + def _checked_theory_vec(self, params): + """Evaluate theory_fn and enforce a 1-D vector aligned to sacc_data.mean.""" + vec = np.asarray(self.theory_fn(params)) + if vec.shape != np.shape(self.sacc_data.mean): + raise ValueError( + f"theory_fn returned shape {vec.shape} but sacc_data.mean has " + f"shape {np.shape(self.sacc_data.mean)}; theory_fn must return " + f"a 1-D vector aligned to the SACC rows." + ) + return vec + def _overlay_shifts(self, shifts): """Overlay drawn deltas on the fiducial point: fiducial[k] + shift[k].""" concealed = deepcopy(self.fiducial_params) @@ -152,9 +166,9 @@ def calculate_concealing_factor(self, factor_type="add"): """ self.factor_type = factor_type - # fiducial vector computed once in __init__ (length guard); recompute + # fiducial vector computed once in __init__ (shape guard); recompute # nothing there --- reuse it and evaluate only the hidden point here. - self.theory_vec_conceal = np.asarray(self.theory_fn(self.__concealed_params)) + self.theory_vec_conceal = self._checked_theory_vec(self.__concealed_params) if factor_type == "add": self.__concealing_factor = self.theory_vec_conceal - self.theory_vec_fid diff --git a/tests/test_datavector.py b/tests/test_datavector.py index 335fc54..72d26a0 100644 --- a/tests/test_datavector.py +++ b/tests/test_datavector.py @@ -61,6 +61,27 @@ def test_length_guard_raises(sacc_data): ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, theory_fn=wrong) +def test_shape_guard_rejects_column_vector(sacc_data): + # (n, 1) has len() == n but would broadcast mean(n,) + factor(n,1) into an + # (n, n) matrix — must be rejected, not silently corrupt the blinded SACC + wrong = lambda p: np.ones((5, 1)) # noqa: E731 + with pytest.raises(ValueError, match="shape"): + ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, theory_fn=wrong) + + +def test_unknown_shift_key_raises(sacc_data, theory_fn): + # a typo'd shift key must fail loudly at construction + with pytest.raises(ValueError, match="omega_c"): + ConcealDataVector(FIDUCIAL, {"omega_c": 0.05}, sacc_data, + seed=SEED, theory_fn=theory_fn) + + +def test_seed_is_required(sacc_data, theory_fn): + # custody contract: no default seed — omitting it must be a TypeError + with pytest.raises(TypeError): + ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, theory_fn=theory_fn) + + # --- concealing factor (add) bit-for-bit ------------------------------------ def test_concealing_factor_add_bit_for_bit(sacc_data, theory_fn): @@ -175,18 +196,22 @@ def test_end_to_end_synthetic(tmp_path, sacc_data, theory_fn, monkeypatch): # --- default CCL backend ---------------------------------------------------- def _make_cosmic_shear_sacc(): - """Small cosmic-shear SACC with 2 NZ tracers and xi_plus/xi_minus rows.""" + """Small cosmic-shear SACC with 2 NZ tracers, xi± and cl_ee rows.""" s = sacc.Sacc() z = np.linspace(0.05, 2.0, 50) for name, zmean in [("src0", 0.5), ("src1", 1.0)]: nz = np.exp(-0.5 * ((z - zmean) / 0.2) ** 2) s.add_tracer("NZ", name, z, nz) thetas = np.array([5.0, 20.0, 60.0, 120.0]) # arcmin + ells = [50, 200, 800] pairs = [("src0", "src0"), ("src0", "src1"), ("src1", "src1")] for dt in ("galaxy_shear_xi_plus", "galaxy_shear_xi_minus"): for p in pairs: for th in thetas: s.add_data_point(dt, p, 1e-6, theta=th) + for p in pairs: + for ell in ells: + s.add_data_point("galaxy_shear_cl_ee", p, 1e-9, ell=ell) n = len(s.mean) s.add_covariance(np.eye(n) * 1e-12) return s @@ -204,6 +229,9 @@ def test_default_ccl_backend(tmp_path, monkeypatch): cdv = ConcealDataVector(fiducial, shifts, s, seed=SEED, theory_fn=None) assert np.all(np.isfinite(cdv.theory_vec_fid)) + # the Cl branch (np.interp over the ell grid) must produce real power + cl_idx = s.indices("galaxy_shear_cl_ee") + assert np.all(cdv.theory_vec_fid[cl_idx] > 0) cdv.calculate_concealing_factor(factor_type="add") blinded = cdv.apply_concealing_to_likelihood_datavec() assert not np.array_equal(blinded, s.mean) From fea3bafcb40de88ac2b85b58b13a991fc6468b7f Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 11 Jul 2026 00:05:39 +0200 Subject: [PATCH 7/7] Verify factor_type wiring in the datavector_main CLI test assert_called_once_with(factor_type='add') instead of a bare call check, matching the shift_distr verification the gaussian variant already does. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01KpaRHk3QwN13myduQ3hJyf --- tests/test_main.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_main.py b/tests/test_main.py index 6b64835..4a12ba3 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -35,7 +35,8 @@ def test_datavector_main(mock_load_sacc, mock_smokescreen, mock_encrypt, mock_pr mock_smokescreen.assert_called_once_with( FIDUCIAL, SHIFTS, sacc_file, seed=seed, shift_distr=shift_distribution, input_format='fits') - mock_smoke_inst.calculate_concealing_factor.assert_called_once() + mock_smoke_inst.calculate_concealing_factor.assert_called_once_with( + factor_type='add') mock_smoke_inst.apply_concealing_to_likelihood_datavec.assert_called_once() mock_smoke_inst.save_concealed_datavector.assert_called_once_with( path_to_output, 'cosmicshear_sacc', output_format='fits', suffix=None