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..441c512 --- /dev/null +++ b/src/smokescreen/backends/ccl.py @@ -0,0 +1,146 @@ +# 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``. + +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 +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..ae1c470 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,156 @@ .. 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, 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 - """ - 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 + # 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 - # 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) + # 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) - # create concealed cosmology object: - self.__concealed_cosmo = self._create_concealed_cosmo(self.__shifts) + # shape guard: theory vector must align to the SACC rows + self.theory_vec_fid = self._checked_theory_vec(self.fiducial_params) - if 'debug' in kwargs and kwargs['debug']: - self._debug = True + if self._debug: 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') + print(f"[DEBUG] Concealed params: {self.__concealed_params}") - 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: + 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( - "Likelihood has covariance but user-provided SACC has None for covariance." + 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 _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. - - 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__ (shape guard); recompute + # nothing there --- reuse it and evaluate only the hidden point here. + self.theory_vec_conceal = self._checked_theory_vec(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 +181,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 +202,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 +256,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/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 diff --git a/src/smokescreen/param_shifts.py b/src/smokescreen/param_shifts.py index e6d91fd..fd62375 100644 --- a/src/smokescreen/param_shifts.py +++ b/src/smokescreen/param_shifts.py @@ -6,149 +6,118 @@ .. 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 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 ---------- - 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)`` 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 - 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") + + base_seed = _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") + 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, list)): + if len(envelope) != 2: + raise ValueError( + f"Shift envelope for {key!r} must have length 2 " + f"(lo, hi); got {envelope!r}" + ) + 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, list)): + raise ValueError( + f"Gaussian shift for {key} requires a float half-width, " + f"got {envelope!r}" + ) + 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) 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..72d26a0 100644 --- a/tests/test_datavector.py +++ b/tests/test_datavector.py @@ -1,1040 +1,240 @@ -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 theory_fn(): + return _synthetic_theory_fn(5) - def __getitem__(self, key): - return self._params[key] +# --- construction & length guard -------------------------------------------- -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} +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 can be instantiated with valid inputs - smokescreen = ConcealDataVector(cosmo, likelihood, - shifts_dict, sacc_data, systematics_dict) - assert isinstance(smokescreen, 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)} + ConcealDataVector(FIDUCIAL, SHIFTS, sacc_data, seed=SEED, theory_fn=wrong) - # 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") +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) - # Check that the concealing (blinding) factor is correct - assert concealing_factor == smokescreen.theory_vec_conceal - smokescreen.theory_vec_fid +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_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) +def test_seed_is_required(sacc_data, theory_fn): + # custody contract: no default seed — omitting it must be a TypeError 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, theory_fn=theory_fn) - # 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, 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 + + +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)) + # 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) + + 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..4a12ba3 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -2,163 +2,93 @@ 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_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 ) -@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..a664c7c 100644 --- a/tests/test_param_shifts.py +++ b/tests/test_param_shifts.py @@ -1,144 +1,138 @@ -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_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_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_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 -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_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"] -# 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 - shifts = draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed) +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) + } - assert shifts["Omega_c"] == 0.01 - assert shifts["sigma8"] == 0.02 +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 -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) +# --- delta semantics -------------------------------------------------------- - assert 0.01 <= shifts["Omega_c"] <= 0.02 - assert 0.02 <= shifts["sigma8"] <= 0.03 +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_invalid_key(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"invalid_param": 0.01} - seed = 1234 +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 - with pytest.raises(ValueError, match=r"Key invalid_param not in cosmology parameters"): - draw_flat_or_deterministic_param_shifts(cosmo, shifts_dict, seed) +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 -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 - 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) +# --- seed normalization ----------------------------------------------------- +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) -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) +# --- distributions ---------------------------------------------------------- - assert 0.01 <= shifts["Omega_c"] <= 0.02 - assert 0.02 <= shifts["sigma8"] <= 0.03 +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"] -# 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 +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) + ) - 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_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_gaussian_param_shifts_invalid_key(): - cosmo = ccl.CosmologyVanillaLCDM() - shifts_dict = {"invalid_param": (0.3, 0.01)} - seed = 1234 +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) - with pytest.raises(ValueError, match=r"Key invalid_param not in cosmology parameters"): - draw_gaussian_param_shifts(cosmo, shifts_dict, seed) +def test_unknown_distribution_raises(): + with pytest.raises(NotImplementedError): + draw_param_shifts({"Omega_c": 0.05}, 2112, shift_distr="lorentzian") -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) +# --- back-compat alias ------------------------------------------------------ - -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" + )