diff --git a/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia_sacc.ini b/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia_sacc.ini new file mode 100644 index 00000000..5c31a476 --- /dev/null +++ b/cosmo_inference/cosmosis_config/cosmosis_pipeline_A_ia_sacc.ini @@ -0,0 +1,110 @@ +#parameters used elsewhere in this file +[DEFAULT] +COSMOSIS_DIR = /n23data1/n06data/lgoh/scratch/cosmosis-standard-library_lisa + + +[pipeline] +modules = consistency sample_S8 camb load_nz_sacc photoz_bias linear_alignment projection add_intrinsic 2pt_shear shear_m_bias sacc_like +likelihoods = 2pt_like +extra_output = cosmological_parameters/omega_lambda cosmological_parameters/S_8 cosmological_parameters/sigma_8 cosmological_parameters/omega_m +timing = T +debug = T + +[runtime] +sampler = polychord +verbosity = debug + +[polychord] +live_points = 192 +feedback = 3 +resume = T +base_dir = %(SCRATCH)s/polychord + +[test] + +[output] +format = text +lock = F + +[consistency] +file = %(COSMOSIS_DIR)s/utility/consistency/consistency_interface.py +verbose = F + +[sample_S8] +file = %(COSMOSIS_DIR)s/utility/sample_sigma8/sample_S8.py + +[camb] +file = %(COSMOSIS_DIR)s/boltzmann/camb/camb_interface.py +mode=power +lmax=2508 +feedback=0 +do_reionization=F +kmin=1e-5 +kmax=20.0 +nk=200 +zmax=5.0 +zmax_background=5.0 +nz_background=500 +halofit_version=mead2020_feedback +nonlinear=pk +neutrino_hierarchy=normal +kmax_extrapolate = 500.0 + +[load_nz_sacc] +file = %(COSMOSIS_DIR)s/number_density/load_nz_sacc/load_nz_sacc.py +nz_file = %(SACC_FILE)s +data_sets = source + +[photoz_bias] +file = %(COSMOSIS_DIR)s/number_density/photoz_bias/photoz_bias.py +mode = additive +sample = nz_source +bias_section = nofz_shifts +interpolation = cubic +output_deltaz_section_name = delta_z_out + +[linear_alignment] +file = %(COSMOSIS_DIR)s/intrinsic_alignments/la_model/linear_alignments_interface_znla.py +method = bk_corrected + +[projection] +file = %(COSMOSIS_DIR)s/structure/projection/project_2d.py +ell_min_logspaced = 1.0 +ell_max_logspaced = 25000.0 +n_ell_logspaced = 400 +shear-shear = source-source +shear-intrinsic = source-source +intrinsic-intrinsic = source-source +get_kernel_peaks = F +verbose = F + +[add_intrinsic] +file = %(COSMOSIS_DIR)s/shear/add_intrinsic/add_intrinsic.py +shear-shear=T +position-shear=F +perbin=F + +[2pt_shear] +file = %(COSMOSIS_DIR)s/shear/cl_to_xi_nicaea/nicaea_interface.so +corr_type = 0 ; shear_cl -> shear_xi + +[shear_m_bias] +file = %(COSMOSIS_DIR)s/shear/shear_bias/shear_m_bias.py +m_per_bin = True +; Despite the parameter name, this can operate on xi as well as C_ell. +cl_section = shear_xi_plus shear_xi_minus +verbose = F + +; Native SACC likelihood via the sp_validation shim (arcmin->rad + ordering +; guard over CosmoSIS's SaccClLikelihood). data_sets/angle ranges use the SACC +; grammar (full data-type names + tracer pairs). like_name=2pt_like keeps the +; block keys identical to the 2pt_like path so chain post-processing is unchanged. +[sacc_like] +file = %(SP_VALIDATION_MODULES)s/sacc_like_unions.py +csl_dir = %(COSMOSIS_DIR)s +data_file = %(SACC_FILE)s +data_sets = galaxy_shear_xi_plus galaxy_shear_xi_minus +like_name = 2pt_like + +angle_range_galaxy_shear_xi_plus_source_0_source_0 = 10.0 200.0 +angle_range_galaxy_shear_xi_minus_source_0_source_0 = 20.0 200.0 diff --git a/papers/cosmo_val/config/config.yaml b/papers/cosmo_val/config/config.yaml index 945f2f8f..a7f19a8d 100644 --- a/papers/cosmo_val/config/config.yaml +++ b/papers/cosmo_val/config/config.yaml @@ -116,11 +116,14 @@ harmonic: binning: powspace nbins: 32 -# Cosmological inference data-product locations (dormant subsystem). +# Cosmological inference data-product locations + tooling. inference: chains_dir: "/n09data/guerrini/output_chains" glass_mock_data_dir: "/n09data/guerrini/glass_mock_v1.4.6/results" glass_mock_chains_dir: "/n09data/guerrini/glass_mock_chains" + # CosmoSIS Standard Library checkout — fills COSMOSIS_DIR in the generated + # pipeline inis (the module `file =` paths and the sacc_like shim's csl_dir). + csl_dir: "/n23data1/n06data/lgoh/scratch/cosmosis-standard-library_lisa" cosebis: theta_min: 1.0 diff --git a/pyproject.toml b/pyproject.toml index ecdec6f1..eef183fb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -138,6 +138,15 @@ blinding = [ # sp_validation fast suite. "numpy>=2.2,<2.5", ] +# Native SACC inference likelihood (PR 7, sp_validation.sacc_like_unions). The +# shim subclasses CosmoSIS's SaccClLikelihood; this pins the engine. Both engine +# module files are pure-python (no CAMB / compiled CSL), so the equality tests +# (tests/test_sacc_like.py) run with just this extra + a CosmoSIS Standard +# Library checkout, located via the CSL_DIR env var: +# git clone --depth 1 https://github.com/joezuntz/cosmosis-standard-library CSL +# CSL_DIR=/path/to/CSL pytest ... test_sacc_like.py +# Inert in CI (its Dockerfile installs only [test,glass,blinding]). +cosmosis = ["cosmosis>=3.25"] develop = ["sp_validation[test,docs]"] [tool.pytest.ini_options] diff --git a/src/sp_validation/one_covariance_io.py b/src/sp_validation/one_covariance_io.py new file mode 100644 index 00000000..4c0f2ae4 --- /dev/null +++ b/src/sp_validation/one_covariance_io.py @@ -0,0 +1,256 @@ +"""ONE_COVARIANCE_IO. + +:Name: one_covariance_io.py + +:Description: File-format glue between the SACC data-product layout + (:mod:`sp_validation.sacc_io`) and OneCovariance + (https://github.com/rreischke/OneCovariance). Two directions: + + - **n(z) SACC -> OneCovariance input** (:func:`write_nz`): the + ``source_i`` NZ tracers of an analysis SACC are written as the + combined whitespace-delimited redshift file OneCovariance reads + (column 0 = z grid, then one ``n(z)`` column per tomographic + bin, no bin edges), and a matching ``[redshift]`` config stanza + is returned via :func:`nz_config_stanza`. + + - **OneCovariance output -> SACC covariance blocks** + (:func:`covariance_blocks`): the flat ``covariance_list_*.dat`` + table OneCovariance emits (one row per element pair) is reshaped + into dense square block(s) — reusing + :func:`sp_validation.statistics.cov_from_one_covariance` for the + per-block reshape — and paired with SACC selectors so a caller + can feed them straight to + :func:`sp_validation.sacc_io.assemble_covariance`. + + OneCovariance itself is *not* a dependency: this module only + touches its file formats, verified against the upstream + ``config.ini`` (``rreischke/OneCovariance`` @ main). + + n(z) file format (upstream ``config.ini`` comment, verbatim): + + ``redshift n_1(z) ... n_{N_source}(z)`` + + i.e. a plain whitespace-delimited text file, column 0 the shared + redshift grid and one column per tomographic bin — no ``z_low``/ + ``z_high`` edges (this is the OneCovariance convention, distinct + from the CosmoSIS NZDATA table which *does* carry edges). All + source bins must therefore share one z grid. + + ``[redshift]`` config keys (upstream canonical names): a single + combined file goes in ``zlens_directory`` + ``zlens_file``; + ``value_loc_in_lensbin`` (``mid``/``left``/``right``) says where + in each histogram bin the tabulated ``n(z)`` value sits — ``mid`` + for the bin-centred grids the SACC stores. NOTE: the UNIONS + OneCovariance template driven by + ``cosmo_val/pseudo_cl.py._modify_onecov_config`` writes the older + key names ``z_directory``/``zlens_file`` instead; pass + ``dir_key="z_directory"`` to match that template. +""" + +import os + +import numpy as np + +from . import sacc_io +from .statistics import cov_from_one_covariance + + +def nz_table(s, n_bins): + """Stack the SACC ``source_i`` NZ tracers into a OneCovariance n(z) table. + + Parameters + ---------- + s : sacc.Sacc + SACC holding ``source_0 … source_{n_bins-1}`` NZ tracers. + n_bins : int + Number of tomographic source bins to write. + + Returns + ------- + numpy.ndarray + Array of shape ``(n_z, n_bins + 1)``: column 0 the shared redshift + grid, columns ``1 … n_bins`` the per-bin ``n(z)``. This is the + OneCovariance combined-file layout (``redshift n_1(z) … n_N(z)``). + + Raises + ------ + ValueError + If any source bin is missing, or if the bins do not share one z grid + (OneCovariance's combined file has a single redshift column, so the + grids must agree bin-for-bin). + """ + z0, nz0 = sacc_io.get_nz(s, 0) + z0 = np.asarray(z0, dtype=float) + columns = [z0] + for i in range(n_bins): + if sacc_io.source_name(i) not in s.tracers: + raise ValueError( + f"SACC has no NZ tracer {sacc_io.source_name(i)!r}; cannot write " + f"a {n_bins}-bin OneCovariance n(z) file" + ) + z_i, nz_i = sacc_io.get_nz(s, i) + if not np.array_equal(np.asarray(z_i, dtype=float), z0): + raise ValueError( + f"source bin {i} n(z) grid differs from source bin 0; the " + "OneCovariance combined n(z) file has one shared redshift column" + ) + columns.append(np.asarray(nz_i, dtype=float)) + return np.column_stack(columns) + + +def write_nz(s, path, n_bins, *, dir_key="zlens_directory", header=True): + """Write the OneCovariance combined n(z) input file from a SACC. + + OneCovariance reads the source redshift distribution as a plain + whitespace-delimited text file whose column 0 is the shared redshift grid + and whose remaining columns are the per-bin ``n(z)`` (``redshift n_1(z) + … n_N(z)``) — no ``z_low``/``z_high`` edges. This writes that file from the + SACC ``source_i`` NZ tracers and returns the ``[redshift]`` config stanza + that points OneCovariance at it. + + Parameters + ---------- + s : sacc.Sacc + Analysis SACC with the ``source_i`` NZ tracers. + path : str or pathlib.Path + Output text-file path (overwritten). Its directory + basename become + the ``[redshift]`` directory/file config values. + n_bins : int + Number of tomographic source bins to write. + dir_key : str, optional + Config key for the redshift directory. Default ``"zlens_directory"`` + (upstream canonical). Pass ``"z_directory"`` for the UNIONS template + driven by ``pseudo_cl.py._modify_onecov_config``. + header : bool, optional + If ``True`` (default) prepend a ``# redshift n_1(z) …`` comment header + naming the columns; OneCovariance's ``genfromtxt``-style reader ignores + it. Set ``False`` for a bare numeric file. + + Returns + ------- + dict + The ``[redshift]`` config stanza (see :func:`nz_config_stanza`), naming + the file just written. + """ + table = nz_table(s, n_bins) + head = "" + if header: + cols = " ".join(f"n_{i + 1}(z)" for i in range(n_bins)) + head = f"redshift {cols}" + np.savetxt(str(path), table, header=head) + return nz_config_stanza( + os.path.dirname(os.path.abspath(str(path))), + os.path.basename(str(path)), + dir_key=dir_key, + ) + + +def nz_config_stanza( + directory, filename, *, dir_key="zlens_directory", value_loc="mid" +): + """Build the OneCovariance ``[redshift]`` config stanza for an n(z) file. + + Parameters + ---------- + directory : str + Directory holding the n(z) file (OneCovariance ``*_directory`` value). + filename : str + n(z) file basename (OneCovariance ``zlens_file`` value). + dir_key : str, optional + Directory config key — ``"zlens_directory"`` (upstream) or + ``"z_directory"`` (UNIONS template). Default ``"zlens_directory"``. + value_loc : str, optional + ``value_loc_in_lensbin`` — where in each histogram bin the tabulated + ``n(z)`` value sits (``mid``/``left``/``right``). Default ``"mid"``, + matching the bin-centred grids the SACC stores. + + Returns + ------- + dict + The ``[redshift]`` key/value pairs: ``{dir_key: directory, "zlens_file": + filename, "value_loc_in_lensbin": value_loc}``. Assign these under + ``config["redshift"]`` of a OneCovariance ``configparser`` config. + """ + if value_loc not in ("mid", "left", "right"): + raise ValueError( + f"value_loc_in_lensbin must be 'mid', 'left' or 'right'; got {value_loc!r}" + ) + return { + dir_key: directory, + "zlens_file": filename, + "value_loc_in_lensbin": value_loc, + } + + +def read_nz(path): + """Read a OneCovariance combined n(z) file back to ``(z, nz_columns)``. + + Inverse of :func:`write_nz` (the numeric round-trip; the config stanza is + not stored in the file). Comment/header lines are skipped. + + Parameters + ---------- + path : str or pathlib.Path + n(z) text file (column 0 = z, columns 1… = per-bin n(z)). + + Returns + ------- + tuple + ``(z, nz)`` where ``z`` is the shared redshift grid (shape ``(n_z,)``) + and ``nz`` is the per-bin distributions (shape ``(n_z, n_bins)``). + """ + table = np.atleast_2d(np.genfromtxt(str(path))) + return table[:, 0], table[:, 1:] + + +def covariance_blocks(cov_list, selectors, *, gaussian=True): + """Reshape a OneCovariance ``covariance_list`` table into SACC cov blocks. + + OneCovariance emits a flat ``covariance_list_*.dat`` table with one row per + ``(i, j)`` element pair (row-major, ``k = i·n + j``); the covariance value + lives in column 10 (Gaussian) or column 9 (Gaussian+non-Gaussian). This + reshapes the flat table into dense square block(s) — reusing + :func:`sp_validation.statistics.cov_from_one_covariance` for the per-block + reshape — and pairs each with its SACC selector, ready for + :func:`sp_validation.sacc_io.assemble_covariance`. + + Single-statistic case: pass the whole table and one selector; you get one + ``(selector, dense)`` block. Multi-statistic case (tomography-ready): pass a + sequence of ``(selector, sub_table)`` pairs — each ``sub_table`` a + contiguous slice of the flat output for one statistic / bin-pair — and each + is reshaped and re-paired with its selector in order. The API is thus shaped + to extend to multi-probe blocking without over-fitting the single-bin case. + + Parameters + ---------- + cov_list : numpy.ndarray or sequence + Either the flat OneCovariance table (2-D array, one row per pair) for a + single block, or — for the multi-block form — a sequence of + ``(selector, sub_table)`` pairs. In the multi-block form ``selectors`` + must be ``None`` (the selectors travel with the sub-tables). + selectors : selector or None + For the single-block form, the SACC selector for the whole table (a + ``(data_type, tracers[, tags])`` tuple or an index array, as + :func:`sacc_io.assemble_covariance` accepts). Must be ``None`` for the + multi-block form. + gaussian : bool, optional + Select the Gaussian-only column (``True``, default) or the + Gaussian+non-Gaussian column (``False``); passed straight through to + ``cov_from_one_covariance``. + + Returns + ------- + list + Ordered ``(selector, dense_cov)`` pairs, directly consumable by + ``sacc_io.assemble_covariance(s, blocks)``. + """ + if selectors is None: + # Multi-block form: cov_list is a sequence of (selector, sub_table). + return [ + (selector, cov_from_one_covariance(np.asarray(sub), gaussian=gaussian)) + for selector, sub in cov_list + ] + # Single-block form: one flat table, one selector. + return [ + (selectors, cov_from_one_covariance(np.asarray(cov_list), gaussian=gaussian)) + ] diff --git a/src/sp_validation/sacc_like_unions.py b/src/sp_validation/sacc_like_unions.py new file mode 100644 index 00000000..a25837cb --- /dev/null +++ b/src/sp_validation/sacc_like_unions.py @@ -0,0 +1,225 @@ +"""SACC likelihood shim for CosmoSIS — the sp_validation-owned native path. + +This module is a thin subclass of CosmoSIS's ``SaccClLikelihood`` (from the +CosmoSIS Standard Library, ``likelihood/sacc/sacc_like.py``) that fixes two +upstream defects so the native SACC likelihood matches the PR-3 converter → +``2pt_like`` path bit for bit on our real-space ξ± analysis file. It exists as a +CosmoSIS *module file* (``setup``/``execute``/``cleanup`` at module scope), loaded +via ``file = .../sacc_like_unions.py`` in an ini, and is NOT imported by +``sp_validation/__init__`` (cosmosis is an optional dependency). + +Why the shim exists +------------------- +Two things in upstream ``sacc_like`` break a real-space ξ likelihood; both are +documented and empirically verified (probe: Δχ²=3184 raw vs Δχ²=0 shimmed on the +single-bin analysis file), and both are fixed here by overriding ``build_data``: + +1. **No arcmin→radian conversion (the killer).** + ``sacc_likelihoods/twopoint.py`` (L74-90 at CSL commit 4fd2f1c) builds a + ``SpectrumInterp`` over ``block[section, "theta"]`` — theory θ in *radians* + (CosmoSIS convention) — and evaluates it at each data point's raw ``theta`` + tag, which our SACC files (and firecrown) store in *arcmin*. ``2pt_like.py`` + (L179-183) converts its real-space data to radians for exactly this reason; + ``sacc_like`` never does, so the spline is evaluated ~3437× outside its grid, + ``SpectrumInterp`` returns 0 there, and χ² silently collapses to dᵀC⁻¹d. The + only prior use was the ℓ-space (unit-free) Cℓ path, which is why it was never + caught. We evaluate the theory against a radian-θ *copy* of the loaded SACC, + keeping ``self.sacc_data`` in arcmin so upstream's save_theory / + save_realization paths (which copy it and overwrite only values) never write + radian tags into a file downstream consumers read as arcmin. + +2. **Theory↔data ordering is assumed, never enforced.** + The data vector is ``sacc.get_mean()`` (insertion order); the theory vector is + built by looping ``get_data_types() × get_tracer_combinations() × points``. + A comment in ``twopoint.py`` (L35) claims ``to_canonical_order`` was called on + load — it is not. The two orders agree only when the file is grouped + type-major (all of one data type, then the next). Our single-pair + ``[ξ+; ξ−]`` files satisfy this; a tomographic *pair-major* file (ξ+/ξ− per + pair) would silently misalign theory against data. We reconstruct the + theory-loop index order and require it equals ``arange`` — a hard ValueError + otherwise (see ``test_ordering_guard_raises_pair_major_tomographic``). + +When upstream fixes the units, this shim dies: the tripwire test +``test_upstream_unit_gap_tripwire`` in ``tests/test_sacc_like.py`` fails the day +raw ``sacc_like`` stops producing a wildly different χ², signalling the shim can +be retired. + +Requires a CSL checkout: ``setup`` reads ``csl_dir`` from the module options and +imports the upstream ``sacc_like`` from ``/likelihood/sacc``. The tests +locate it via the ``CSL_DIR`` environment variable (see the test module docstring +for the checkout recipe). +""" + +import os +import sys + +import numpy as np + +# cosmosis is an OPTIONAL dependency: this module is a CosmoSIS module file, but +# importing it must succeed without cosmosis installed (the CI image has the +# science stack but no cosmosis, and test_imports.py bare-imports every module). +# So every cosmosis touch is deferred into setup() — nothing at top level imports +# it. numpy is fine at top level (always present). + +# arcmin → radian: the conversion 2pt_like applies to real-space data and that +# sacc_like omits. Applied only to `theta` tags of `real`-category data types. +ARCMIN_TO_RAD = np.pi / (180.0 * 60.0) + + +def _import_upstream_sacc_like(csl_dir): + """Import the upstream ``sacc_like`` module from a CSL checkout. + + ``/likelihood/sacc`` is prepended to ``sys.path`` so both + ``sacc_like`` and its sibling ``sacc_likelihoods`` package (imported by + ``sacc_like`` for the theory-extraction functions) resolve. Idempotent: the + path is only inserted once. + """ + sacc_dir = os.path.join(csl_dir, "likelihood", "sacc") + if not os.path.isdir(sacc_dir): + raise ValueError( + f"csl_dir={csl_dir!r} has no likelihood/sacc directory; point csl_dir " + "at a CosmoSIS Standard Library checkout (see module docstring)" + ) + if sacc_dir not in sys.path: + sys.path.insert(0, sacc_dir) + import sacc_like # noqa: E402 — resolved from the sys.path insertion above + + return sacc_like + + +def _make_subclass(sacc_like): + """Build ``SaccLikeUnions`` as a subclass of the upstream ``SaccClLikelihood``. + + A factory (not an import-time ``class ... :`` statement) so importing this + module never requires the upstream class — that dependency is deferred to + ``setup``, which has ``csl_dir`` in hand. Only ``build_data`` is overridden; + scale cuts, covariance handling (Sellentin/Hartlap), theory extraction and + ``save_theory`` all ride upstream unmodified. + """ + + class SaccLikeUnions(sacc_like.SaccClLikelihood): + """CSL ``SaccClLikelihood`` with the arcmin→rad + ordering-guard fixes.""" + + def build_data(self): + # Run the upstream build FIRST: it loads the SACC, applies data_sets + # selection and the arcmin-grammar scale cuts (matching 2pt_like's + # angle_range convention and the ini ergonomics), populates + # self.sacc_data / self.sections_for_names, and returns the + # (unit-independent) data vector we pass straight through. + x, data_vector = super().build_data() + + # self.sacc_data STAYS in arcmin — save_theory / save_realization copy + # it and only overwrite point values, so its θ tags must remain the + # units the file was written in (any consumer, incl. re-ingesting the + # saved SACC as a data_file, assumes arcmin). The radian conversion the + # theory spline needs lives on a separate copy, swapped in only for the + # extraction (see extract_theory_points). + self._sacc_data_rad = self._radian_theta_copy(self.sacc_data) + self._assert_theory_order_matches_data() + + return x, data_vector + + def _radian_theta_copy(self, sacc_data): + """A copy of ``sacc_data`` with ``real``-category θ tags in radians. + + The theory spline is built on ``block[section, "theta"]`` in radians, + so the ``theta`` tag each point is evaluated against must be radians + too. Scoped to data types whose category + (``sections_for_names[dt][0]``) is ``real``: cosebis ``n`` tags and + spectrum ``ell`` tags are unit-free and left untouched, and only the + ``theta`` tag is scaled (``theta_nom`` etc. are metadata the likelihood + never evaluates). Operates on a ``.copy()`` so the original stays + arcmin for the save paths. + """ + real_types = { + dt + for dt, (category, _section) in self.sections_for_names.items() + if category == "real" + } + converted = sacc_data.copy() + for point in converted.data: + if point.data_type in real_types and "theta" in point.tags: + point.tags["theta"] = point.tags["theta"] * ARCMIN_TO_RAD + return converted + + def extract_theory_points(self, block): + """Extract theory against the radian-θ copy, then restore the original. + + Upstream ``extract_theory_points`` reads ``self.sacc_data`` (the θ tag + per point) to evaluate the theory spline; that read needs radians. + Swap in ``self._sacc_data_rad`` for the duration of the upstream call + and restore in ``finally`` so everything else — including the + save_theory / save_realization copies that run afterward in + ``do_likelihood`` — sees the untouched arcmin ``self.sacc_data``. + """ + original = self.sacc_data + self.sacc_data = self._sacc_data_rad + try: + return super().extract_theory_points(block) + finally: + self.sacc_data = original + + def _assert_theory_order_matches_data(self): + """Require the theory-loop order to equal the data-vector order. + + The data vector is ``sacc.get_mean()`` (insertion order); upstream + builds theory by looping data types, then tracer combinations, then + points, and concatenating — assuming (never enforcing) that this + reproduces insertion order. It does only for type-major files. We + reconstruct that loop's index order and require it be ``arange``; + otherwise theory and data would be silently misaligned (the same bug + class the PR-2/PR-3 reviews caught for the converter). + """ + order = [ + int(i) + for dt in self.sacc_data.get_data_types() + for tracers in self.sacc_data.get_tracer_combinations(dt) + for i in self.sacc_data.indices(dt, tracers) + ] + expected = np.arange(len(self.sacc_data.mean)) + if not np.array_equal(order, expected): + raise ValueError( + "SACC data/theory ordering mismatch: the theory loop " + "(get_data_types × get_tracer_combinations × points) does not " + "reproduce the get_mean() insertion order, so sacc_like would " + "compare theory against data point-by-point in the WRONG order " + "and return a silently wrong χ². This happens when the file is " + "grouped pair-major (ξ+/ξ− interleaved per tracer pair) rather " + "than type-major (all ξ+, then all ξ−). Upstream assumes " + "to_canonical_order() was applied on load but never calls it; " + "write the SACC type-major, or call to_canonical_order() before " + "saving." + ) + + return SaccLikeUnions + + +def setup(options): + """CosmoSIS ``setup`` — build and instantiate the shimmed likelihood. + + Mirrors ``GaussianLikelihood.build_module``'s setup: wrap the raw options in + ``SectionOptions`` and instantiate the likelihood (whose ``__init__`` calls + ``build_data``). The one addition is reading ``csl_dir`` from the module + options to locate and import the upstream class before subclassing it. The + cosmosis import is deferred to here (call time) so importing this module never + requires cosmosis. + """ + from cosmosis.datablock import SectionOptions, option_section + + csl_dir = options.get_string(option_section, "csl_dir") + sacc_like = _import_upstream_sacc_like(csl_dir) + likelihood_class = _make_subclass(sacc_like) + return likelihood_class(SectionOptions(options)) + + +def execute(block, config): + """CosmoSIS ``execute`` — run the likelihood (mirrors ``build_module``).""" + likelihood_calculator = config + likelihood_calculator.do_likelihood(block) + return 0 + + +def cleanup(config): + """CosmoSIS ``cleanup`` — mirror of ``build_module``'s cleanup.""" + likelihood_calculator = config + likelihood_calculator.cleanup() diff --git a/src/sp_validation/tests/test_bmodes_workflow_dry_run.py b/src/sp_validation/tests/test_bmodes_workflow_dry_run.py index 87afc416..91a1923d 100644 --- a/src/sp_validation/tests/test_bmodes_workflow_dry_run.py +++ b/src/sp_validation/tests/test_bmodes_workflow_dry_run.py @@ -96,3 +96,32 @@ def test_cosmo_val_workflow_assemble_dry_runs(): assert f"pseudo_cl_cov_{version}_blind=A_powspace_nbins=32.fits" in out, out for part in ("_xi_coarse_", "_cosebis.sacc", "_pure_eb.sacc", "rho_tau_"): assert part in out, f"missing {part} part in assemble DAG:\n{out}" + + +@requires_candide_data +def test_cosmo_val_inference_prep_dry_runs(): + """The revived (PR-7) inference_prep DAG resolves end to end from the SACC. + + inference_fiducial must pull inference_prep, which consumes the assembled + {version}.sacc and emits the converter 2pt-FITS plus BOTH generated pipeline + inis (2pt_like and the native sacc_like). The old cosmosis_fitting.py real- + data assembly is retired from this path; the glass-mock rules keep it. + """ + version = "SP_v1.4.6.3_leak_corr" + result = _dry_run(_repo_root() / "papers/cosmo_val", ["inference_fiducial"]) + assert result.returncode == 0, result.stdout + out = result.stdout + assert "rule inference_prep:" in out, out + assert "rule inference_fiducial:" in out, out + # inference_prep consumes the assembled analysis SACC (not per-sign xi FITS). + assert f"{version}.sacc" in out, out + # ...and both ini TEMPLATES, bound as inputs so a template edit regenerates the + # configs (Finding 2 — a template as params gave no DAG edge). + assert "cosmosis_pipeline_A_ia.ini" in out, out + assert "cosmosis_pipeline_A_ia_sacc.ini" in out, out + # It emits the converter FITS + both engine inis. + assert f"cosmosis_{version}.fits" in out, out + assert f"cosmosis_pipeline_{version}_A_ia.ini" in out, out + assert f"cosmosis_pipeline_{version}_A_ia_sacc.ini" in out, out + # The retired real-data assembly script must not appear in this DAG's prep. + assert "cosmosis_fitting.py --cosmosis-root" not in out, out diff --git a/src/sp_validation/tests/test_generate_inference_config.py b/src/sp_validation/tests/test_generate_inference_config.py new file mode 100644 index 00000000..4fbbb98d --- /dev/null +++ b/src/sp_validation/tests/test_generate_inference_config.py @@ -0,0 +1,138 @@ +"""Tests for the CosmoSIS inference-config generator. + +The generator (:mod:`workflow.scripts.generate_inference_config`) fills a +pipeline ini template's ``[DEFAULT]`` section with concrete paths so CosmoSIS's +ConfigParser resolves the template's ``%(KEY)s`` placeholders. These tests need +no cosmosis: they check that the substituted DEFAULT keys land, that the module +file paths the templates reference exist on disk, and that no ``%(...)s`` +placeholder is left unresolved after filling. +""" + +import configparser +import importlib.util +from pathlib import Path + +import pytest + +_REPO = Path(__file__).resolve().parents[3] +_SCRIPT = _REPO / "workflow" / "scripts" / "generate_inference_config.py" +_CONFIG_DIR = _REPO / "cosmo_inference" / "cosmosis_config" +_SACC_TEMPLATE = _CONFIG_DIR / "cosmosis_pipeline_A_ia_sacc.ini" +_FITS_TEMPLATE = _CONFIG_DIR / "cosmosis_pipeline_A_ia.ini" + + +def _load_generator(): + spec = importlib.util.spec_from_file_location("gen_inference_cfg", _SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +gen = _load_generator() + + +def _read_interpolated(ini_path): + """Parse the generated ini with interpolation ON (the way CosmoSIS reads it). + + CosmoSIS uses ``%(KEY)s`` BasicInterpolation with case-*preserved* keys, so + set ``optionxform = str`` (stdlib configparser lowercases keys by default, + which would break the uppercase ``%(FITS_FILE)s`` / ``%(SACC_FILE)s`` + lookups). The default BasicInterpolation then raises if a referenced key is + missing — exactly the failure we want to catch. + """ + parser = configparser.ConfigParser() + parser.optionxform = str + parser.read(ini_path) + return parser + + +def test_sacc_template_defaults_filled(tmp_path): + """The generated sacc ini carries the substituted DEFAULT keys and resolves.""" + out = tmp_path / "gen_sacc.ini" + gen.generate_inference_config( + _SACC_TEMPLATE, + out, + gen._substitutions( + scratch="/scratch/run", + cosmosis_dir="/csl", + sacc_file="/data/v1.sacc", + ), + ) + text = out.read_text() + assert "SCRATCH = /scratch/run" in text + assert "SACC_FILE = /data/v1.sacc" in text + assert "COSMOSIS_DIR = /csl" in text + assert "SP_VALIDATION_MODULES = " in text + # FITS_FILE is None for the sacc path — must be dropped, not written as "None". + assert "FITS_FILE" not in text + + parser = _read_interpolated(out) + # The sacc_like data_file interpolates SACC_FILE; load_nz_sacc too. + assert parser["sacc_like"]["data_file"] == "/data/v1.sacc" + assert parser["load_nz_sacc"]["nz_file"] == "/data/v1.sacc" + assert parser["sacc_like"]["csl_dir"] == "/csl" + + +def test_fits_template_defaults_filled(tmp_path): + """The generated 2pt_like ini carries the substituted DEFAULT keys and resolves.""" + out = tmp_path / "gen_fits.ini" + gen.generate_inference_config( + _FITS_TEMPLATE, + out, + gen._substitutions( + scratch="/scratch/run", + cosmosis_dir="/csl", + fits_file="/data/v1.fits", + ), + ) + text = out.read_text() + assert "SCRATCH = /scratch/run" in text + assert "FITS_FILE = /data/v1.fits" in text + assert "SACC_FILE" not in text + + parser = _read_interpolated(out) + assert parser["2pt_like"]["data_file"] == "/data/v1.fits" + assert parser["load_nz_fits"]["nz_file"] == "/data/v1.fits" + + +def test_sacc_template_module_file_exists(): + """The sacc_like module file the generated ini points at exists on disk. + + ``SP_VALIDATION_MODULES`` resolves to the installed package dir; + ``sacc_like_unions.py`` must live there (it is the shim CosmoSIS loads). + """ + modules = Path(gen._sp_validation_modules()) + assert (modules / "sacc_like_unions.py").is_file() + + +def test_all_placeholders_resolve(tmp_path): + """No ``%(...)s`` placeholder survives interpolation in either template. + + A missing DEFAULT key would make ConfigParser raise on access; iterate every + option in every section to force resolution of all placeholders. + """ + for template, subs in ( + ( + _SACC_TEMPLATE, + gen._substitutions(scratch="/s", cosmosis_dir="/csl", sacc_file="/d.sacc"), + ), + ( + _FITS_TEMPLATE, + gen._substitutions(scratch="/s", cosmosis_dir="/csl", fits_file="/d.fits"), + ), + ): + out = tmp_path / (template.stem + ".gen.ini") + gen.generate_inference_config(template, out, subs) + parser = _read_interpolated(out) + for section in parser.sections(): + for key in parser[section]: + value = parser[section][key] # raises if a placeholder is unresolved + assert "%(" not in value, f"[{section}] {key} = {value}" + + +def test_no_default_section_raises(tmp_path): + """A template with no [DEFAULT] section is a loud error.""" + bad = tmp_path / "bad.ini" + bad.write_text("[pipeline]\nmodules = a b c\n") + with pytest.raises(ValueError, match="DEFAULT"): + gen.generate_inference_config(bad, tmp_path / "out.ini", {"SCRATCH": "/s"}) diff --git a/src/sp_validation/tests/test_one_covariance_io.py b/src/sp_validation/tests/test_one_covariance_io.py new file mode 100644 index 00000000..a99ad552 --- /dev/null +++ b/src/sp_validation/tests/test_one_covariance_io.py @@ -0,0 +1,274 @@ +"""Tests for :mod:`sp_validation.one_covariance_io`. + +All synthetic, all fast: the OneCovariance fixtures are built in memory +shaped exactly like its real file I/O — a flat ``covariance_list`` table with +the ``(i, j)`` index rows and the Gaussian / Gauss+non-Gaussian value columns +at the indices ``cov_from_one_covariance`` expects (col 10 / col 9), and the +combined n(z) text file (column 0 = z, one column per bin, no edges). No +cluster paths; OneCovariance is not imported. + +The two pieces: + +- **Piece 1** (n(z) SACC -> OneCovariance input): write the combined n(z) + file from a SACC's ``source_i`` NZ tracers, read it back, assert the z and + n(z) columns round-trip and the config stanza names the file. +- **Piece 2** (OneCovariance output -> SACC covariance blocks): reshape the + flat table to the dense block(s) and prove they feed + ``sacc_io.assemble_covariance`` cleanly (reshaped block -> assemble -> + ``s.covariance.dense`` matches the hand-built matrix). +""" + +import numpy as np +import numpy.testing as npt +import pytest + +from sp_validation import one_covariance_io as ocio +from sp_validation import sacc_io as sio + + +# --------------------------------------------------------------------------- # +# Synthetic OneCovariance-shaped fixtures +# --------------------------------------------------------------------------- # +def _one_cov_table(cov_gauss, cov_all): + """Flatten two n x n matrices into a OneCovariance ``covariance_list`` table. + + Reproduces the real flat output: one row per ``(i, j)`` element pair in + row-major order ``k = i·n + j``, with the Gaussian value in column 10 and + the Gaussian+non-Gaussian value in column 9. Columns 0-8 and the index + columns are filled with self-documenting placeholder values (the reshape + only reads cols 9/10, but a realistic width proves it does not spill). + """ + n = cov_gauss.shape[0] + rows = [] + for i in range(n): + for j in range(n): + row = np.arange(11.0) # placeholder cols 0-8 (+ overwritten 9,10) + row[9] = cov_all[i, j] + row[10] = cov_gauss[i, j] + rows.append(row) + return np.array(rows) + + +def _spd(n, seed): + """Symmetric positive-definite matrix of size ``n`` (a valid covariance).""" + a = np.random.default_rng(seed).normal(size=(n, n)) + return a @ a.T + n * np.eye(n) + + +def _nz(seed, n=40): + rng = np.random.default_rng(seed) + z = np.linspace(0.01, 2.0, n) + return z, rng.uniform(0.1, 1.0, n) + + +# --------------------------------------------------------------------------- # +# Piece 2 — flat covariance_list -> dense block (reshape correctness + teeth) +# --------------------------------------------------------------------------- # +def test_covariance_blocks_reshapes_to_hand_built_matrix(): + """Pin the reshape and prove gaussian vs gauss+ng select different columns. + + WHAT IS PINNED: ``covariance_blocks`` flattens/reshapes the OneCovariance + ``covariance_list`` table into a dense square block matching a hand-built + covariance. It delegates the per-block reshape to + ``statistics.cov_from_one_covariance``, so column 10 (gaussian) and column 9 + (gauss+ng) must recover the two distinct hand-built matrices. + + WHY TEETH: (a) ``gaussian=True`` vs ``False`` must return the two *different* + matrices, proving the column flag is load-bearing; (b) perturbing a single + entry of the flat input must change exactly that entry of the reshaped + block, proving the reshape actually reads the table (not a constant). + """ + cov_gauss = _spd(4, seed=1) + cov_all = _spd(4, seed=2) + table = _one_cov_table(cov_gauss, cov_all) + + selector = (sio.XI_PLUS, (sio.source_name(0), sio.source_name(0))) + + [(sel_g, block_g)] = ocio.covariance_blocks(table, selector, gaussian=True) + [(sel_a, block_a)] = ocio.covariance_blocks(table, selector, gaussian=False) + + assert sel_g == selector and sel_a == selector + npt.assert_allclose(block_g, cov_gauss, rtol=1e-12) + npt.assert_allclose(block_a, cov_all, rtol=1e-12) + + # TEETH: gaussian and gauss+ng select different columns -> different blocks. + assert not np.allclose(block_g, block_a) + + # TEETH: a perturbation of one flat-table entry moves exactly that block + # entry (row k = i·n + j, col 10 for gaussian). + perturbed = table.copy() + perturbed[2 * 4 + 1, 10] += 5.0 # element (i=2, j=1) + [(_, block_p)] = ocio.covariance_blocks(perturbed, selector, gaussian=True) + npt.assert_allclose(block_p[2, 1] - block_g[2, 1], 5.0, rtol=1e-12) + block_p[2, 1] = block_g[2, 1] + npt.assert_allclose(block_p, block_g, rtol=1e-12) # nothing else moved + + +def test_covariance_blocks_multiblock_form(): + """Prove the tomography-ready multi-block form reshapes each sub-table. + + WHAT IS PINNED: passing ``selectors=None`` and a sequence of + ``(selector, sub_table)`` pairs reshapes each sub-table independently and + returns them paired with their selectors in order — the shape needed to map + a multi-statistic OneCovariance output onto several SACC selectors. + + WHY TEETH: the two sub-tables carry different matrices; if the function + reshaped only the first or mixed them, the second block would not match its + own hand-built matrix. + """ + cov_a, cov_b = _spd(3, seed=3), _spd(2, seed=4) + table_a = _one_cov_table(cov_a, _spd(3, seed=5)) + table_b = _one_cov_table(cov_b, _spd(2, seed=6)) + sel_a = (sio.XI_PLUS, (sio.source_name(0), sio.source_name(0))) + sel_b = (sio.XI_MINUS, (sio.source_name(0), sio.source_name(0))) + + blocks = ocio.covariance_blocks( + [(sel_a, table_a), (sel_b, table_b)], None, gaussian=True + ) + + assert [s for s, _ in blocks] == [sel_a, sel_b] + npt.assert_allclose(blocks[0][1], cov_a, rtol=1e-12) + npt.assert_allclose(blocks[1][1], cov_b, rtol=1e-12) + + +# --------------------------------------------------------------------------- # +# Piece 2 — the reshaped block feeds assemble_covariance cleanly +# --------------------------------------------------------------------------- # +def test_covariance_blocks_feed_assemble_covariance(): + """Round-trip: reshaped block -> assemble_covariance -> dense matches. + + WHAT IS PINNED: the ``(selector, dense)`` pair ``covariance_blocks`` + returns is directly consumable by ``sacc_io.assemble_covariance``: assembled + onto a SACC whose only statistic is one ξ+ block, ``s.covariance.dense`` + must equal the hand-built OneCovariance matrix. This is the end-to-end + contract between the two modules. + + WHY TEETH: the block must tile the data vector exactly; if the reshape + produced the wrong size or the wrong selector, ``assemble_covariance`` would + raise (its contiguity/tiling/size checks), so a clean assemble + matching + dense is a real proof. + """ + theta = np.geomspace(1.0, 100.0, 4) + xip, xim = np.arange(4) * 1e-5, np.arange(4) * 2e-5 + s = sio.new_sacc({0: _nz(0)}) + sio.add_xi(s, (0, 0), theta, xip, xim, grid="coarse") + + # The ξ block spans ξ+ then ξ− for the pair -> 8 points, one contiguous + # block (pair-major, matching sacc_io's canonical order). + cov_gauss = _spd(8, seed=7) + table = _one_cov_table(cov_gauss, _spd(8, seed=8)) + pair = (sio.source_name(0), sio.source_name(0)) + selector = np.concatenate( + [s.indices(sio.XI_PLUS, pair), s.indices(sio.XI_MINUS, pair)] + ) + + blocks = ocio.covariance_blocks(table, selector, gaussian=True) + sio.assemble_covariance(s, blocks) + + npt.assert_allclose(s.covariance.dense, cov_gauss, rtol=1e-12) + + +# --------------------------------------------------------------------------- # +# Piece 1 — n(z) SACC -> OneCovariance input (round-trip + config stanza) +# --------------------------------------------------------------------------- # +def test_write_nz_roundtrips_and_names_file(tmp_path): + """Round-trip the n(z) file and check the config stanza names it. + + WHAT IS PINNED: ``write_nz`` writes the SACC ``source_i`` NZ tracers as the + OneCovariance combined file (column 0 = z, one column per bin, no z_low/ + z_high edges). ``read_nz`` recovers the z grid and every per-bin n(z) + column, and the returned ``[redshift]`` stanza names the exact file written + (directory + basename) plus ``value_loc_in_lensbin``. + + WHY TEETH: the z grid and each n(z) column must round-trip to the values the + SACC holds (drawn from a seeded RNG); a transposed write or a dropped column + would fail the per-bin comparison. The stanza's directory/file must match + the path actually written. + """ + z, nz0 = _nz(10) + _, nz1 = _nz(11) + s = sio.new_sacc({0: (z, nz0), 1: (z, nz1)}) + + path = tmp_path / "nz_onecov.txt" + stanza = ocio.write_nz(s, path, n_bins=2) + + z_read, nz_read = ocio.read_nz(path) + npt.assert_allclose(z_read, z, rtol=1e-12) + assert nz_read.shape == (len(z), 2) + npt.assert_allclose(nz_read[:, 0], nz0, rtol=1e-12) + npt.assert_allclose(nz_read[:, 1], nz1, rtol=1e-12) + + assert stanza["zlens_directory"] == str(tmp_path) + assert stanza["zlens_file"] == "nz_onecov.txt" + assert stanza["value_loc_in_lensbin"] == "mid" + + +def test_write_nz_unions_template_dir_key(tmp_path): + """The UNIONS-template ``z_directory`` key is selectable via ``dir_key``. + + WHAT IS PINNED: the upstream OneCovariance key is ``zlens_directory``, but + the UNIONS template (pseudo_cl.py._modify_onecov_config) writes + ``z_directory``. ``dir_key="z_directory"`` produces that variant so the + stanza drops straight into the UNIONS template's ``[redshift]`` section. + """ + z, nz0 = _nz(12) + s = sio.new_sacc({0: (z, nz0)}) + stanza = ocio.write_nz(s, tmp_path / "nz.txt", n_bins=1, dir_key="z_directory") + assert "z_directory" in stanza and "zlens_directory" not in stanza + assert stanza["z_directory"] == str(tmp_path) + assert stanza["zlens_file"] == "nz.txt" + + +def test_write_nz_fails_on_mismatched_z_grids(tmp_path): + """Fail fast when source bins do not share one redshift grid. + + WHAT IS PINNED: the OneCovariance combined file has a single redshift + column, so all bins must share the z grid. A bin on a different grid must + raise ``ValueError`` at write time, not silently mis-align. + """ + z0, nz0 = _nz(20) + z1_shifted, nz1 = _nz(21) + z1_shifted = z1_shifted + 0.1 # different grid + s = sio.new_sacc({0: (z0, nz0), 1: (z1_shifted, nz1)}) + with pytest.raises(ValueError, match="differs from source bin 0"): + ocio.write_nz(s, tmp_path / "bad.txt", n_bins=2) + + +def test_write_nz_no_header_roundtrips(tmp_path): + """The bare (header=False) file still round-trips numerically. + + WHAT IS PINNED: ``header=False`` writes a purely numeric file (no ``#`` + column-name line); ``read_nz`` recovers the same z grid and n(z) column, so + the header is cosmetic and never load-bearing for the numeric round-trip. + """ + z, nz0 = _nz(40) + s = sio.new_sacc({0: (z, nz0)}) + path = tmp_path / "bare.txt" + ocio.write_nz(s, path, n_bins=1, header=False) + z_read, nz_read = ocio.read_nz(path) + npt.assert_allclose(z_read, z, rtol=1e-12) + npt.assert_allclose(nz_read[:, 0], nz0, rtol=1e-12) + + +def test_nz_config_stanza_rejects_bad_value_loc(): + """``value_loc_in_lensbin`` outside {mid,left,right} fails fast. + + WHAT IS PINNED: OneCovariance only accepts ``mid``/``left``/``right`` for + the histogram-bin value location; an invalid value is a config bug and must + raise ``ValueError`` rather than write a stanza OneCovariance will reject. + """ + with pytest.raises(ValueError, match="value_loc_in_lensbin"): + ocio.nz_config_stanza("/dir", "nz.txt", value_loc="center") + + +def test_write_nz_fails_on_missing_bin(tmp_path): + """Fail fast when a requested source bin is absent from the SACC. + + WHAT IS PINNED: requesting more bins than the SACC carries is a real config + bug; ``write_nz`` raises ``ValueError`` naming the missing tracer rather + than writing a short file. + """ + z, nz0 = _nz(30) + s = sio.new_sacc({0: (z, nz0)}) + with pytest.raises(ValueError, match="source_1"): + ocio.write_nz(s, tmp_path / "short.txt", n_bins=2) diff --git a/src/sp_validation/tests/test_sacc_like.py b/src/sp_validation/tests/test_sacc_like.py new file mode 100644 index 00000000..8e8c1332 --- /dev/null +++ b/src/sp_validation/tests/test_sacc_like.py @@ -0,0 +1,566 @@ +"""Equality tests: the sp_validation SACC-likelihood shim vs CosmoSIS ``2pt_like``. + +PR 7 adopts CosmoSIS's native ``SaccClLikelihood`` for the ξ± inference path, +through the shim :mod:`sp_validation.sacc_like_unions` (which fixes the upstream +arcmin→rad gap and adds an ordering guard). The contract is *in-process module +equality*: run the shimmed ``sacc_like`` on the analysis SACC and CosmoSIS's +``2pt_like`` on the PR-3 converter's 2pt-FITS against an identical synthetic +theory DataBlock, and require the same χ², log-likelihood, theory vector and +post-cut point count. The prototype (``sacc-like-probe/probe_equality.py``) +observed *exact* equality (Δχ²=0, Δtheory=0), so the equality tests assert +``array_equal`` / rtol=1e-12. + +Environment +----------- +Both engines are pure-python CosmoSIS module files (no CAMB, no compiled CSL +modules), so they run in the shared venv with ``cosmosis`` installed. They need a +checkout of the CosmoSIS Standard Library, located via the ``CSL_DIR`` env var; +the module is skipped when cosmosis is absent (CI image) or ``CSL_DIR`` is unset +or missing. Recipe:: + + git clone --depth 1 https://github.com/joezuntz/cosmosis-standard-library CSL + CSL_DIR=/path/to/CSL pytest ... test_sacc_like.py + +The pyproject ``cosmosis`` extra pins ``cosmosis>=3.25`` for the engine itself +(inert in CI, whose Dockerfile installs only ``[test,glass,blinding]``). +""" + +import importlib.util +import os +from pathlib import Path + +import numpy as np +import pytest + +cosmosis = pytest.importorskip("cosmosis") + +# The upstream CSL checkout carrying likelihood/sacc + likelihood/2pt. Skip the +# whole module (not error) when it is not configured, mirroring the cosmo_numba +# env-gated precedent. +_CSL_DIR = os.environ.get("CSL_DIR") +if not _CSL_DIR or not Path(_CSL_DIR, "likelihood", "sacc").is_dir(): + pytest.skip( + "CSL_DIR unset or has no likelihood/sacc — set CSL_DIR to a CosmoSIS " + "Standard Library checkout to run the sacc_like equality tests", + allow_module_level=True, + ) + +from cosmosis.datablock import DataBlock, option_section # noqa: E402 + +from sp_validation import sacc_io, twopoint_convert # noqa: E402 + +# Reuse the angular-bin count from the converter tests so the two suites' single- +# bin shapes stay in lockstep. (The χ²-dynamics builders here need a covariance +# commensurate with the data, which the converter's byte-compare _sacc is not.) +from sp_validation.tests.test_twopoint_convert import N_ANG # noqa: E402 + +CSL = Path(_CSL_DIR) +ARCMIN_TO_RAD = np.pi / (180.0 * 60.0) + +# The like_name both engines are configured with, so both write identical block +# keys (_CHI2, _LIKE, _theory) — the parity the design mandates. +LIKE_NAME = "2pt_like" + + +# --------------------------------------------------------------------------- +# Shared engine harness +# --------------------------------------------------------------------------- +def _load_module_file(path, name): + """Import a CosmoSIS module file (``setup``/``execute``/``cleanup``) by path. + + ``likelihood/sacc`` and ``likelihood/2pt`` are added to ``sys.path`` so the + upstream modules resolve their siblings (``sacc_likelihoods``, ``spec_tools``, + ``twopoint_cosmosis``, …). Loading ``2pt_like.py`` runs its + ``build_module()`` at import (module-level ``setup, execute, cleanup``); + ``sacc_like_unions.py`` defines those functions directly. + """ + import sys + + for sub in ("likelihood/sacc", "likelihood/2pt"): + p = str(CSL / sub) + if p not in sys.path: + sys.path.insert(0, p) + spec = importlib.util.spec_from_file_location(name, path) + mod = importlib.util.module_from_spec(spec) + spec.loader.exec_module(mod) + return mod + + +def _theory_block(): + """A DataBlock carrying the synthetic ξ± theory both engines interpolate. + + Mirrors the probe: the theory θ grid is in *radians* (CosmoSIS convention), + and the ξ+/ξ− predictions are smooth power laws sampled on it. Both engines + build a spline over this grid and evaluate it at each data point's angular + tag — so as long as both see the same block, the interpolated theory (and + hence χ²) must agree. + """ + theta_arcmin = np.geomspace(0.3, 400.0, 300) + + def t_xip(th): + return 2e-4 * (th / 10.0) ** -0.8 + + def t_xim(th): + return 1e-4 * (th / 10.0) ** -0.5 + + b = DataBlock() + for section, f in (("shear_xi_plus", t_xip), ("shear_xi_minus", t_xim)): + b[section, "theta"] = theta_arcmin * ARCMIN_TO_RAD + b[section, "bin_1_1"] = f(theta_arcmin) + b[section, "is_auto"] = True + b[section, "nbin_a"] = 1 + b[section, "nbin_b"] = 1 + b[section, "sample_a"] = "nz_source" + b[section, "sample_b"] = "nz_source" + b[section, "sep_name"] = "theta" + b[section, "save_name"] = "" + return b + + +def _run(mod, options): + """Run a CosmoSIS likelihood module file standalone; return (like, chi2, theory, n). + + Builds the options DataBlock (module options live under ``option_section``), + calls ``setup`` → ``execute`` on a fresh theory block, and reads the standard + Gaussian-likelihood outputs back out under the shared ``LIKE_NAME`` keys. + """ + opt = DataBlock() + for key, value in options.items(): + opt[option_section, key] = value + config = mod.setup(opt) + block = _theory_block() + mod.execute(block, config) + like = block["likelihoods", f"{LIKE_NAME}_LIKE"] + chi2 = block["data_vector", f"{LIKE_NAME}_CHI2"] + theory = block["data_vector", f"{LIKE_NAME}_theory"] + return like, chi2, np.asarray(theory), len(theory) + + +# --------------------------------------------------------------------------- +# The two engine module files + their options +# --------------------------------------------------------------------------- +_SHIM_PATH = Path(__file__).resolve().parents[1] / "sacc_like_unions.py" + + +@pytest.fixture(scope="module") +def m_shim(): + return _load_module_file(_SHIM_PATH, "sacc_like_unions_test") + + +@pytest.fixture(scope="module") +def m_2pt(): + return _load_module_file(CSL / "likelihood/2pt/2pt_like.py", "twopt_like_test") + + +@pytest.fixture(scope="module") +def m_raw_sacc(): + return _load_module_file(CSL / "likelihood/sacc/sacc_like.py", "raw_sacc_like_test") + + +def _shim_opts(sacc_path, **extra): + return { + "csl_dir": str(CSL), + "data_file": sacc_path, + "data_sets": "galaxy_shear_xi_plus galaxy_shear_xi_minus", + "like_name": LIKE_NAME, + **extra, + } + + +def _twopt_opts(fits_path, **extra): + return { + "data_file": fits_path, + "data_sets": "XI_PLUS XI_MINUS", + "covmat_name": "COVMAT", + "like_name": LIKE_NAME, + "gaussian_covariance": False, + "cut_zeros": False, + **extra, + } + + +def _realistic_sacc(seed=0, *, xip=None): + """A single-bin ξ± SACC with data + covariance sized like the real product. + + The converter-test ``_sacc`` builder uses a covariance ~14 orders of + magnitude larger than the ξ values (fine for byte-comparing the converter, + where covariance *content* is irrelevant), which makes every χ² collapse to + numerical zero — no teeth for the unit-gap tripwire. This builder instead + lays down realistic ξ± power laws and a covariance ~ ``(0.1·|ξ|)²`` (the + probe's recipe), so χ² is O(1)-scale and the raw-vs-shim gap is visible. + + ``xip`` overrides the ξ+ values (perturbation teeth). + """ + rng = np.random.default_rng(seed) + n = N_ANG + theta = np.geomspace(1.0, 250.0, n) # arcmin + z = np.linspace(0.01, 3.0, 200) + nz = z**2 * np.exp(-((z / 0.5) ** 1.5)) + + xip_vals = 2e-4 * (theta / 10.0) ** -0.8 * (1 + 0.05 * rng.standard_normal(n)) + xim_vals = 1e-4 * (theta / 10.0) ** -0.5 * (1 + 0.05 * rng.standard_normal(n)) + if xip is not None: + xip_vals = xip + + s = sacc_io.new_sacc({0: (z, nz)}, {"catalogue_version": "test"}) + sacc_io.add_xi(s, (0, 0), theta, xip_vals, xim_vals, grid="coarse") + + sig = 0.1 * np.abs(np.concatenate([xip_vals, xim_vals])) + a = rng.standard_normal((2 * n, 3 * n)) + cov = (a @ a.T / (3 * n)) * np.outer(sig, sig) * 0.3 + np.diag(sig**2) + s.add_covariance(cov) + return s, theta + + +def _write_pair(tmp_path, seed=0, name="probe", *, xip=None): + """Write a realistic analysis SACC and its PR-3 converter 2pt-FITS. + + Returns ``(sacc_path, fits_path)``. The SACC carries the arcmin ξ± tags the + shim converts; the FITS is the byte-compatible product ``2pt_like`` reads. + """ + s, _theta = _realistic_sacc(seed, xip=xip) + sacc_path = str(tmp_path / f"{name}.sacc") + sacc_io.save(s, sacc_path) + fits_path = str(tmp_path / f"{name}_2pt.fits") + twopoint_convert.sacc_to_twopoint_fits(sacc_io.load(sacc_path), fits_path, n_bins=1) + return sacc_path, fits_path + + +# --------------------------------------------------------------------------- +# 1. Core equality: shim on SACC ≡ 2pt_like on converter FITS +# --------------------------------------------------------------------------- +def test_shimmed_equals_2pt_like_exact(tmp_path, m_shim, m_2pt): + """The shimmed sacc_like and 2pt_like agree exactly on the same data+theory. + + Same synthetic theory block, the SACC through the shim vs the converter FITS + through 2pt_like: χ², log-likelihood and the theory vector must match to + numerical precision (the probe observed exact equality). This is the PR's + central contract — the native path reproduces the validated converter path. + """ + sacc_path, fits_path = _write_pair(tmp_path, seed=0) + + like_s, chi2_s, theory_s, n_s = _run(m_shim, _shim_opts(sacc_path)) + like_t, chi2_t, theory_t, n_t = _run(m_2pt, _twopt_opts(fits_path)) + + assert n_s == n_t == 2 * N_ANG + np.testing.assert_allclose(chi2_s, chi2_t, rtol=1e-12) + np.testing.assert_allclose(like_s, like_t, rtol=1e-12) + np.testing.assert_allclose(theory_s, theory_t, rtol=1e-12) + + +# --------------------------------------------------------------------------- +# 2. Tripwire: raw upstream sacc_like is broken on arcmin tags +# --------------------------------------------------------------------------- +def test_upstream_unit_gap_tripwire(tmp_path, m_raw_sacc, m_2pt): + """RAW upstream sacc_like (no shim) gives a wildly wrong χ² on arcmin tags. + + Documents and guards the arcmin→rad gap the shim fixes: with raw ``theta`` + tags the theory spline is evaluated outside its (radian) grid and returns 0, + collapsing χ² to dᵀC⁻¹d. We assert the relative χ² difference against + 2pt_like exceeds 10 (the probe saw Δχ²≈3184). + + IF THIS TEST EVER FAILS: upstream ``sacc_like`` has fixed its units — the + shim's arcmin→rad conversion is now redundant and the shim can be retired. + """ + sacc_path, fits_path = _write_pair(tmp_path, seed=0) + + _like_r, chi2_raw, _theory_r, _n_r = _run(m_raw_sacc, _shim_opts(sacc_path)) + _like_t, chi2_t, _theory_t, _n_t = _run(m_2pt, _twopt_opts(fits_path)) + + rel = abs(chi2_raw - chi2_t) / abs(chi2_t) + assert rel > 10, ( + f"raw sacc_like χ²={chi2_raw:.6g} is within 10× of 2pt_like χ²={chi2_t:.6g} " + f"(rel diff {rel:.3g}) — the upstream unit gap appears fixed; retire the shim" + ) + + +# --------------------------------------------------------------------------- +# 3. Scale cuts equivalent through both engines +# --------------------------------------------------------------------------- +def test_scale_cuts_equivalent(tmp_path, m_shim, m_2pt): + """The same arcmin scale cuts give the same post-cut N and χ² on both engines. + + Cuts are expressed in the each engine's grammar but the SAME numeric arcmin + values (shim cuts run before the arcmin→rad conversion, so they take arcmin + just like 2pt_like's angle_range). ξ+ ∈ [10, 200], ξ− ∈ [20, 200] arcmin. + """ + sacc_path, fits_path = _write_pair(tmp_path, seed=1) + + shim_cuts = _shim_opts( + sacc_path, + **{ + "angle_range_galaxy_shear_xi_plus_source_0_source_0": np.array( + [10.0, 200.0] + ), + "angle_range_galaxy_shear_xi_minus_source_0_source_0": np.array( + [20.0, 200.0] + ), + }, + ) + twopt_cuts = _twopt_opts( + fits_path, + **{ + "angle_range_XI_PLUS_1_1": np.array([10.0, 200.0]), + "angle_range_XI_MINUS_1_1": np.array([20.0, 200.0]), + }, + ) + + _like_s, chi2_s, _theory_s, n_s = _run(m_shim, shim_cuts) + _like_t, chi2_t, _theory_t, n_t = _run(m_2pt, twopt_cuts) + + assert n_s == n_t + assert n_s < 2 * N_ANG # the cuts actually removed points + np.testing.assert_allclose(chi2_s, chi2_t, rtol=1e-12) + + +# --------------------------------------------------------------------------- +# 4. Perturbation moves both engines identically +# --------------------------------------------------------------------------- +def test_perturbation_moves_both_identically(tmp_path, m_shim, m_2pt): + """Perturbing one data value shifts both engines' χ² by the identical amount. + + Teeth: build the base pair and a pair whose first ξ+ value is bumped, and + require Δχ²(shim) == Δχ²(2pt_like). If either engine ignored the perturbed + point (e.g. a misaligned data vector), the deltas would diverge. + """ + # Base pair, then a pair whose first ξ+ value is bumped. Rebuild the base ξ+ + # from the same seed so only the one perturbed entry differs. + base_s, theta = _realistic_sacc(seed=2) + base_xip = np.array( + [p.value for p in base_s.data if p.data_type == sacc_io.XI_PLUS] + ) + sacc_b, fits_b = _write_pair(tmp_path, seed=2, name="base") + + pert_xip = base_xip.copy() + pert_xip[0] += 5e-5 + sacc_p, fits_p = _write_pair(tmp_path, seed=2, name="pert", xip=pert_xip) + + _l, chi2_shim_b, _t, _n = _run(m_shim, _shim_opts(sacc_b)) + _l, chi2_shim_p, _t, _n = _run(m_shim, _shim_opts(sacc_p)) + _l, chi2_2pt_b, _t, _n = _run(m_2pt, _twopt_opts(fits_b)) + _l, chi2_2pt_p, _t, _n = _run(m_2pt, _twopt_opts(fits_p)) + + d_shim = chi2_shim_p - chi2_shim_b + d_2pt = chi2_2pt_p - chi2_2pt_b + assert abs(d_shim) > 0 # the perturbation actually moved χ² + np.testing.assert_allclose(d_shim, d_2pt, rtol=1e-10) + + +# --------------------------------------------------------------------------- +# 5. Ordering guard raises on a pair-major tomographic file +# --------------------------------------------------------------------------- +def test_ordering_guard_raises_pair_major_tomographic(tmp_path, m_shim): + """A pair-major 2-bin SACC trips the shim's ordering guard at setup. + + Inserting ξ± per pair — (0,0), then (0,1), then (1,1) — lays the data vector + out pair-major ([ξ+;ξ−] per pair), while ``get_data_types()`` groups the + theory loop type-major (all ξ+ pairs, then all ξ− pairs). The two orders + disagree (verified empirically), so the guard must raise a ValueError + mentioning the ordering hazard rather than silently mis-comparing. + """ + theta = np.geomspace(1.0, 250.0, N_ANG) # arcmin + z = np.linspace(0.01, 3.0, 200) + nz = z**2 * np.exp(-((z / 0.5) ** 1.5)) + xip = np.ones(N_ANG) * 1e-4 + xim = np.ones(N_ANG) * 1e-4 + s = sacc_io.new_sacc({0: (z, nz), 1: (z, nz)}) + for pair in [(0, 0), (0, 1), (1, 1)]: + sacc_io.add_xi(s, pair, theta, xip, xim, grid="coarse") + s.add_covariance(np.eye(len(s.mean))) + sacc_path = str(tmp_path / "pair_major.sacc") + sacc_io.save(s, sacc_path) + + with pytest.raises(ValueError, match="ordering"): + m_shim.setup(_as_option_block(_shim_opts(sacc_path))) + + +def _as_option_block(options): + """Build a raw options DataBlock (keys under ``option_section``) for setup().""" + opt = DataBlock() + for key, value in options.items(): + opt[option_section, key] = value + return opt + + +# --------------------------------------------------------------------------- +# 6. theta conversion scoped to real-category types (COSEBIs untouched) +# --------------------------------------------------------------------------- +def test_theta_conversion_scoped_to_real_types(tmp_path, m_shim): + """The shim converts ξ ``theta`` tags but leaves cosebi ``n`` tags untouched. + + Build a SACC carrying both ξ± (real) and COSEBIs (a non-real ``cosebis`` + category with integer ``n`` tags). After the shim's ``build_data``, the ξ + ``theta`` tags must be scaled arcmin→rad (so they equal the original arcmin + values times the conversion factor), and the cosebi ``n`` tags must be + numerically unchanged. + """ + # Build ξ± + COSEBIs, then attach the covariance last (sacc forbids adding + # points after add_covariance). + theta = np.geomspace(1.0, 250.0, N_ANG) # arcmin + z = np.linspace(0.01, 3.0, 200) + nz = z**2 * np.exp(-((z / 0.5) ** 1.5)) + n_modes = 5 + En = np.arange(1.0, n_modes + 1) + Bn = np.arange(1.0, n_modes + 1) * 0.1 + + s = sacc_io.new_sacc({0: (z, nz)}) + sacc_io.add_xi( + s, (0, 0), theta, np.ones(N_ANG) * 1e-4, np.ones(N_ANG) * 1e-4, grid="coarse" + ) + sacc_io.add_cosebis(s, (0, 0), En, Bn, scale_cut=(1.0, 250.0)) + s.add_covariance(np.eye(len(s.mean))) + + sacc_path = str(tmp_path / "with_cosebis.sacc") + sacc_io.save(s, sacc_path) + + # data_sets keeps the cosebis in (so we can check its tags survive); cosebi's + # section/category resolve from sacc_like's default_sections, so build_data + # needs no extra ini config. Only setup() runs (build_data); the theory loop + # (which would want a cosebi theory block) runs at execute, not here. + config = m_shim.setup( + _as_option_block( + _shim_opts( + sacc_path, + data_sets=( + "galaxy_shear_xi_plus galaxy_shear_xi_minus " + "galaxy_shear_cosebi_ee galaxy_shear_cosebi_bb" + ), + ) + ) + ) + + # The conversion lives on the radian copy (_sacc_data_rad); the original + # self.sacc_data stays arcmin (so save_theory writes arcmin tags — Finding 1). + rad_xi_thetas = [ + p.tags["theta"] + for p in config._sacc_data_rad.data + if p.data_type == sacc_io.XI_PLUS + ] + rad_cosebi_ns = [ + p.tags["n"] + for p in config._sacc_data_rad.data + if p.data_type == sacc_io.COSEBI_EE + ] + orig_xi_thetas = [ + p.tags["theta"] for p in config.sacc_data.data if p.data_type == sacc_io.XI_PLUS + ] + # rad copy: ξ theta scaled to radians (original arcmin × ARCMIN_TO_RAD). + np.testing.assert_allclose( + np.sort(rad_xi_thetas), np.sort(theta * ARCMIN_TO_RAD), rtol=1e-12 + ) + # rad copy: cosebi n tags untouched (still integer modes 1..n_modes). + np.testing.assert_array_equal(np.sort(rad_cosebi_ns), np.arange(1, n_modes + 1)) + # original sacc_data: ξ theta still in arcmin (unmutated). + np.testing.assert_allclose(np.sort(orig_xi_thetas), np.sort(theta), rtol=1e-12) + + +# --------------------------------------------------------------------------- +# 6b. save_theory writes arcmin tags, and re-execute is stable (Finding 1) +# --------------------------------------------------------------------------- +def test_save_theory_writes_arcmin_and_reexecute_stable(tmp_path, m_shim): + """save_theory must write a SACC whose θ tags are still arcmin, not radians. + + Finding 1: because upstream save_theory copies self.sacc_data and overwrites + only point values, self.sacc_data must stay arcmin — otherwise the saved file + carries radian θ tags and any arcmin-assuming consumer (sacc_io.get_xi, or + re-ingesting it as a data_file, which would double-convert to ~8.5e-8) is + silently off by 3437×. Assert the saved θ tags match the input file's tags + exactly (arcmin) and that the saved values equal the theory vector. Also run + execute() twice and require identical χ² — a guard against any accidental + double-conversion creeping back in. + """ + s, theta = _realistic_sacc(seed=5) + sacc_path = str(tmp_path / "in.sacc") + sacc_io.save(s, sacc_path) + input_theta = np.array( + [p.tags["theta"] for p in sacc_io.load(sacc_path).data if "theta" in p.tags] + ) + + save_path = str(tmp_path / "saved_theory.sacc") + opt = _as_option_block(_shim_opts(sacc_path, save_theory=save_path)) + config = m_shim.setup(opt) + + block1 = _theory_block() + m_shim.execute(block1, config) + chi2_1 = block1["data_vector", f"{LIKE_NAME}_CHI2"] + theory = np.asarray(block1["data_vector", f"{LIKE_NAME}_theory"]) + + saved = sacc_io.load(save_path) + saved_theta = np.array([p.tags["theta"] for p in saved.data if "theta" in p.tags]) + saved_values = np.array(saved.mean) + + # θ tags in the saved file are arcmin — identical to the input file's tags. + np.testing.assert_array_equal(saved_theta, input_theta) + # and are NOT the radian conversion (guards against the leak explicitly). + assert not np.allclose(saved_theta, input_theta * ARCMIN_TO_RAD) + # saved values are the theory vector (save_theory overwrites values in order). + np.testing.assert_allclose(saved_values, theory, rtol=1e-12) + + # A second execute() yields the identical χ² — no cumulative mutation. + block2 = _theory_block() + m_shim.execute(block2, config) + chi2_2 = block2["data_vector", f"{LIKE_NAME}_CHI2"] + np.testing.assert_allclose(chi2_2, chi2_1, rtol=1e-12) + + +# --------------------------------------------------------------------------- +# 7. Real-data equality (candide-gated) +# --------------------------------------------------------------------------- +_REALDATA = ( + Path("/automnt/n17data/cdaley/unions/code/sp_validation/cosmo_inference/data") + / "SP_v1.4.6_leak_corr_A_minsep=1.0_maxsep=250.0_nbins=20_npatch=1" + / "cosmosis_SP_v1.4.6_leak_corr_A_minsep=1.0_maxsep=250.0_nbins=20_npatch=1.fits" +) + + +@pytest.mark.skipif( + not _REALDATA.exists(), reason=f"real 2pt-FITS not on disk: {_REALDATA}" +) +def test_realdata_shim_equals_2pt_like(tmp_path, m_shim, m_2pt): + """On a real product, the shim on its SACC equals 2pt_like on the FITS. + + Builds a ξ-only analysis SACC from the real 2pt-FITS's own ξ± values and + covariance sub-block (via the converter test's ``_sacc_from_2pt_fits``, then + strip to ξ±), writes it, converts it to a plain-ξ FITS, and runs both engines + against the synthetic theory block. Scoping to ξ± isolates the shear + likelihood equality on the true (20-point-per-sign) data-vector shape — the + IA-only inference scope this PR targets — and keeps ``2pt_like``'s + ``twopoint.from_fits`` from tripping over the real file's separate + COVMAT_CELL / τ blocks. + """ + from astropy.io import fits + + from sp_validation.tests.test_twopoint_convert_realdata import _sacc_from_2pt_fits + + with fits.open(_REALDATA) as hdul: + full_s, _rho_hdu, _tau_hdu = _sacc_from_2pt_fits(hdul) + + # Rebuild a ξ-only SACC: same n(z), ξ± values and ξ± covariance sub-block. + source = sacc_io.source_name(0) + z, nz = sacc_io.get_nz(full_s, 0) + theta, xip, xim = sacc_io.get_xi(full_s, (0, 0), grid="coarse") + xi_idx = np.concatenate( + [ + full_s.indices(sacc_io.XI_PLUS, (source, source)), + full_s.indices(sacc_io.XI_MINUS, (source, source)), + ] + ) + xi_cov = full_s.covariance.dense[np.ix_(xi_idx, xi_idx)] + + s = sacc_io.new_sacc({0: (z, nz)}) + sacc_io.add_xi(s, (0, 0), theta, xip, xim, grid="coarse") + s.add_covariance(xi_cov) + + sacc_path = str(tmp_path / "real_xi.sacc") + sacc_io.save(s, sacc_path) + fits_path = str(tmp_path / "real_xi_2pt.fits") + twopoint_convert.sacc_to_twopoint_fits(sacc_io.load(sacc_path), fits_path, n_bins=1) + + like_s, chi2_s, theory_s, n_s = _run(m_shim, _shim_opts(sacc_path)) + like_t, chi2_t, theory_t, n_t = _run(m_2pt, _twopt_opts(fits_path)) + + assert n_s == n_t + assert n_s == 2 * len(theta) + np.testing.assert_allclose(chi2_s, chi2_t, rtol=1e-10) + np.testing.assert_allclose(like_s, like_t, rtol=1e-10) + np.testing.assert_allclose(theory_s, theory_t, rtol=1e-10) diff --git a/src/sp_validation/tests/test_twopoint_convert.py b/src/sp_validation/tests/test_twopoint_convert.py new file mode 100644 index 00000000..3a5bc44a --- /dev/null +++ b/src/sp_validation/tests/test_twopoint_convert.py @@ -0,0 +1,472 @@ +"""Byte-compare tests for the SACC -> 2pt-FITS converter. + +The converter (:mod:`sp_validation.twopoint_convert`) must reproduce the CosmoSIS +2pt-FITS that ``cosmo_inference/scripts/cosmosis_fitting.py`` assembles today, +so the inference chain (``2pt_like`` and Sacha Guerrini's rho/tau +``2pt_like_xi_sys`` fork) runs untouched behind it. The strongest possible check +is *byte* equality, and astropy writes FITS deterministically, so that is what we +assert: build a reference with the current script's own HDU-builder functions on +deterministic synthetic inputs, build a SACC from those same inputs via +:mod:`sp_validation.sacc_io`, convert it, and compare the two files byte for byte. + +Three configurations pin the three product shapes today's ``__main__`` emits: + +1. **plain xi** -- PRIMARY, NZ_SOURCE, COVMAT, XI_PLUS, XI_MINUS. +2. **xi + pseudo-Cl** -- adds COVMAT_CELL and CELL_EE (the harmonic block + ``2pt_like`` reads; the script builds CELL_BB and discards it, so the + converter does too). +3. **xi + rho/tau** -- adds the blocked tau covariance (TAU_0_PLUS / TAU_2_PLUS, + with the tau_0<->tau_2 cross-correlation the truncated CosmoCov tau + covariance carries) and the verbatim RHO_STATS table. + +The rho/tau product needs the rho/tau *sidecar* HDUs: the RHO_STATS table carries +per-mode ``varrho_*`` variances the analysis SACC does not store, so the +converter copies them from the sidecar exactly as today's assembly does. A teeth +test pins that a perturbed input moves the output. + +The reference builder imports ``cosmosis_fitting.py`` by path (it is a script, +not a package module), skipping cleanly if a dependency is missing -- the same +loader pattern as ``test_cosmosis_fitting.py``. +""" + +import importlib.util +from pathlib import Path + +import numpy as np +import pytest +from astropy.io import fits + +from sp_validation import sacc_io, twopoint_convert + +_SCRIPT = ( + Path(__file__).resolve().parents[3] + / "cosmo_inference" + / "scripts" + / "cosmosis_fitting.py" +) + + +def _load_cf(): + """Import cosmosis_fitting.py by path; skip cleanly if a dep is missing.""" + if not _SCRIPT.exists(): + pytest.skip(f"cosmosis_fitting.py not found at {_SCRIPT}") + spec = importlib.util.spec_from_file_location("cosmosis_fitting", _SCRIPT) + module = importlib.util.module_from_spec(spec) + try: + spec.loader.exec_module(module) + except ImportError as exc: # pragma: no cover - container has numpy/astropy + pytest.importorskip(getattr(exc, "name", "") or "cosmosis_fitting_dependency") + raise + return module + + +cf = _load_cf() + + +# --- deterministic synthetic inputs ----------------------------------------- +# +# One tomographic bin (today's analysis). N_ANG angular bins on an ascending +# theta grid; N_ELL bandpowers; a 200-point n(z) on a uniform z grid (the DES +# NZDATA table needs a uniform Z_MID axis). The xi covariance is a full +# (2*N_ANG) matrix (xi+/xi- cross-block nonzero, as CosmoCov produces); the tau +# covariance is a full (3*N_ANG) matrix that the assembly truncates to its first +# 2 statistics (tau_0, tau_2). + +N_ANG = 5 +N_ELL = 8 +SOURCE = sacc_io.source_name(0) +PSF = sacc_io.PSF_TRACER + + +def _spd(n, seed): + """A symmetric positive-definite (n, n) matrix, seeded and recognizable.""" + rng = np.random.default_rng(seed) + a = rng.normal(size=(n, n)) + return a @ a.T + n * np.eye(n) + + +def _inputs(seed=0): + """Deterministic synthetic statistics + covariances for one bin pair.""" + rng = np.random.default_rng(seed) + theta = np.sort(rng.uniform(1.0, 250.0, N_ANG)) + ell = np.sort(rng.uniform(30.0, 3000.0, N_ELL)) + z = np.linspace(0.0125, 0.4875, 200) + return { + "theta": theta, + "ell": ell, + "z": z, + "nz": np.exp(-((z - 0.25) ** 2) / 0.02), + "xip": rng.uniform(1e-6, 1e-4, N_ANG), + "xim": rng.uniform(1e-6, 1e-4, N_ANG), + "cl_ee": rng.uniform(1e-10, 1e-8, N_ELL), + "cl_bb": rng.uniform(1e-11, 1e-9, N_ELL), + "cl_eb": rng.uniform(-1e-11, 1e-11, N_ELL), + "tau0p": rng.uniform(1e-6, 1e-5, N_ANG), + "tau2p": rng.uniform(1e-6, 1e-5, N_ANG), + "tau0m": rng.uniform(1e-6, 1e-5, N_ANG), + "tau2m": rng.uniform(1e-6, 1e-5, N_ANG), + "xi_cov": _spd(2 * N_ANG, seed + 1), + "cl_cov": _spd(N_ELL, seed + 2), + "tau_cov_full": _spd(3 * N_ANG, seed + 3), + } + + +# --- reference sidecar files + HDUs (the current script's own builders) ------ + + +def _rho_sidecar_hdu(theta, seed=7): + """A rho-stats BinTableHDU with the 25-column layout (values + variances).""" + rng = np.random.default_rng(seed) + columns = [fits.Column(name="theta", format="D", array=theta)] + for k in range(6): + for suffix in ("_p", "_m"): + columns.append( + fits.Column( + name=f"rho_{k}{suffix}", + format="D", + array=rng.uniform(-1e-3, 1e-3, len(theta)), + ) + ) + columns.append( + fits.Column( + name=f"varrho_{k}{suffix}", + format="D", + array=rng.uniform(1e-15, 1e-12, len(theta)), + ) + ) + return fits.BinTableHDU.from_columns(fits.ColDefs(columns)) + + +def _tau_sidecar_hdu(theta, tau0p, tau2p): + """A tau-stats BinTableHDU with theta + tau_0_p + tau_2_p columns.""" + columns = [ + fits.Column(name="theta", format="D", array=theta), + fits.Column(name="tau_0_p", format="D", array=tau0p), + fits.Column(name="tau_2_p", format="D", array=tau2p), + ] + return fits.BinTableHDU.from_columns(fits.ColDefs(columns)) + + +def _reference_fits(tmp_path, inp, *, cl=False, rho_tau=False): + """Build the reference 2pt-FITS with cosmosis_fitting.py's own functions. + + Reproduces the exact ``__main__`` HDU list for the requested configuration, + which is the assembly the converter must match byte for byte. + """ + nz_txt = tmp_path / "nz.txt" + np.savetxt(nz_txt, np.column_stack([inp["z"], inp["nz"]])) + cov_txt = tmp_path / "cov_xi.txt" + np.savetxt(cov_txt, inp["xi_cov"]) + + nz_hdu = cf.nz_to_fits(str(nz_txt)) + xip_hdu = cf._create_2pt_hdu(inp["xip"], inp["theta"], "XI_PLUS", "G+R", "G+R") + xim_hdu = cf._create_2pt_hdu(inp["xim"], inp["theta"], "XI_MINUS", "G-R", "G-R") + + hdu_list = [fits.PrimaryHDU(), nz_hdu] + + if rho_tau: + tau_cov_npy = tmp_path / "cov_tau.npy" + np.save(tau_cov_npy, inp["tau_cov_full"]) + cov_hdu = cf.covdat_to_fits(str(cov_txt), filename_cov_tau=str(tau_cov_npy)) + else: + cov_hdu = cf.covdat_to_fits(str(cov_txt), filename_cov_tau=None) + hdu_list.append(cov_hdu) + + if cl: + cl_block = np.zeros((5, N_ELL)) + cl_block[0], cl_block[1], cl_block[4] = inp["ell"], inp["cl_ee"], inp["cl_bb"] + cl_npy = tmp_path / "cl.npy" + np.save(cl_npy, cl_block) + cl_cov_npy = tmp_path / "cl_cov.npy" + np.save(cl_cov_npy, inp["cl_cov"]) + ell_r, cl_ee_r, cl_bb_r = cf.load_pseudo_cl(str(cl_npy)) + cl_ee_hdu, _cl_bb_hdu = cf.cl_to_fits(ell_r, cl_ee_r, cl_bb_r) + cov_cl_hdu = cf.cov_cl_to_fits(str(cl_cov_npy), cov_hdu="COVAR_FULL") + hdu_list.append(cov_cl_hdu) + + hdu_list.extend([xip_hdu, xim_hdu]) + if cl: + hdu_list.append(cl_ee_hdu) + + if rho_tau: + rho_path = tmp_path / "rho.fits" + fits.HDUList([fits.PrimaryHDU(), _rho_sidecar_hdu(inp["theta"])]).writeto( + rho_path, overwrite=True + ) + tau_path = tmp_path / "tau.fits" + fits.HDUList( + [ + fits.PrimaryHDU(), + _tau_sidecar_hdu(inp["theta"], inp["tau0p"], inp["tau2p"]), + ] + ).writeto(tau_path, overwrite=True) + rho_hdu = cf.rho_to_fits(str(rho_path), theta=inp["theta"]) + tau0_hdu, tau2_hdu = cf.tau_to_fits(str(tau_path), theta=inp["theta"]) + hdu_list.extend([tau0_hdu, tau2_hdu, rho_hdu]) + + out = tmp_path / "reference.fits" + fits.HDUList(hdu_list).writeto(out, overwrite=True) + return out + + +# --- the SACC each configuration is built from ------------------------------ + + +def _sacc(inp, *, cl=False, rho_tau=False): + """Build the analysis SACC the converter reads, matching ``_inputs``. + + The covariance is laid out to match the reference exactly: the xi block is + the full (2*N_ANG) matrix; the Cl block is the EE bandpower covariance; the + tau blocks carry the tau_0<->tau_2 cross-correlation from the truncated + CosmoCov tau covariance. Blocks not consumed by the 2pt-FITS (Cl BB/EB, tau + minus) get an identity block so ``add_covariance`` sees a full matrix. + """ + s = sacc_io.new_sacc({0: (inp["z"], inp["nz"])}) + sacc_io.add_xi(s, (0, 0), inp["theta"], inp["xip"], inp["xim"], grid="coarse") + if cl: + sacc_io.add_pseudo_cl( + s, + (0, 0), + inp["ell"], + inp["cl_ee"], + inp["cl_bb"], + inp["cl_eb"], + window_ells=np.arange(2, 102), + window_weights=np.random.default_rng(9).uniform(0, 1, (100, N_ELL)), + ) + if rho_tau: + sacc_io.add_tau(s, (0, 0), 0, inp["theta"], inp["tau0p"], inp["tau0m"]) + sacc_io.add_tau(s, (0, 0), 2, inp["theta"], inp["tau2p"], inp["tau2m"]) + + n = len(s.mean) + full = np.zeros((n, n)) + ip = s.indices(sacc_io.XI_PLUS, (SOURCE, SOURCE)) + im = s.indices(sacc_io.XI_MINUS, (SOURCE, SOURCE)) + xi_all = np.concatenate([ip, im]) + full[np.ix_(xi_all, xi_all)] = inp["xi_cov"] + + if cl: + iee = s.indices(sacc_io.CL_EE, (SOURCE, SOURCE)) + full[np.ix_(iee, iee)] = inp["cl_cov"] + for dtype in (sacc_io.CL_BB, sacc_io.CL_EB): + idx = s.indices(dtype, (SOURCE, SOURCE)) + full[np.ix_(idx, idx)] = np.eye(N_ELL) + if rho_tau: + t0p = s.indices(sacc_io.TAU_PLUS.format(k=0), (SOURCE, PSF)) + t2p = s.indices(sacc_io.TAU_PLUS.format(k=2), (SOURCE, PSF)) + tau_pp = np.concatenate([t0p, t2p]) + # The joint [tau_0+; tau_2+] block is the truncated CosmoCov tau + # covariance -- cross-correlation kept, matching covdat_to_fits. + full[np.ix_(tau_pp, tau_pp)] = inp["tau_cov_full"][: 2 * N_ANG, : 2 * N_ANG] + for dtype in (sacc_io.TAU_MINUS.format(k=0), sacc_io.TAU_MINUS.format(k=2)): + idx = s.indices(dtype, (SOURCE, PSF)) + full[np.ix_(idx, idx)] = np.eye(N_ANG) + s.add_covariance(full) + return s + + +def _sidecar_hdus(tmp_path, inp): + """Return the (rho_hdu, tau_hdu) sidecar input HDUs for the rho/tau product.""" + rho_path = tmp_path / "rho_in.fits" + fits.HDUList([fits.PrimaryHDU(), _rho_sidecar_hdu(inp["theta"])]).writeto( + rho_path, overwrite=True + ) + tau_path = tmp_path / "tau_in.fits" + fits.HDUList( + [fits.PrimaryHDU(), _tau_sidecar_hdu(inp["theta"], inp["tau0p"], inp["tau2p"])] + ).writeto(tau_path, overwrite=True) + with fits.open(rho_path) as r, fits.open(tau_path) as t: + return r[1].copy(), t[1].copy() + + +# ============================================================================= +# Byte-compare: the three product shapes +# ============================================================================= + + +def test_plain_xi_byte_equal(tmp_path): + """Plain-xi product: converter matches cosmosis_fitting.py byte for byte.""" + inp = _inputs(seed=0) + reference = _reference_fits(tmp_path, inp) + s = _sacc(inp) + out = tmp_path / "converted.fits" + twopoint_convert.sacc_to_twopoint_fits(s, str(out), n_bins=1) + assert out.read_bytes() == reference.read_bytes() + + +def test_xi_cl_byte_equal(tmp_path): + """xi + pseudo-Cl product: COVMAT_CELL + CELL_EE reproduced byte for byte.""" + inp = _inputs(seed=10) + reference = _reference_fits(tmp_path, inp, cl=True) + s = _sacc(inp, cl=True) + out = tmp_path / "converted.fits" + twopoint_convert.sacc_to_twopoint_fits(s, str(out), n_bins=1) + assert out.read_bytes() == reference.read_bytes() + + +def test_xi_rho_tau_byte_equal(tmp_path): + """xi + rho/tau product: blocked tau covariance + verbatim RHO_STATS match.""" + inp = _inputs(seed=20) + reference = _reference_fits(tmp_path, inp, rho_tau=True) + s = _sacc(inp, rho_tau=True) + rho_hdu, tau_hdu = _sidecar_hdus(tmp_path, inp) + out = tmp_path / "converted.fits" + twopoint_convert.sacc_to_twopoint_fits( + s, str(out), rho_stats_hdu=rho_hdu, tau_stats_hdu=tau_hdu, n_bins=1 + ) + assert out.read_bytes() == reference.read_bytes() + + +# ============================================================================= +# Structural + teeth checks +# ============================================================================= + + +def test_tau_covariance_keeps_tau0_tau2_cross(tmp_path): + """The tau covariance block couples tau_0 and tau_2 (not block-diagonal). + + covdat_to_fits truncates the 3-statistic CosmoCov tau covariance to its + first 2 blocks and lays it in as ONE contiguous block, so tau_0<->tau_2 + cross-terms survive. A block-diagonal shortcut would zero them; pin that + the converter keeps them. + """ + inp = _inputs(seed=20) + s = _sacc(inp, rho_tau=True) + rho_hdu, tau_hdu = _sidecar_hdus(tmp_path, inp) + out = tmp_path / "converted.fits" + twopoint_convert.sacc_to_twopoint_fits( + s, str(out), rho_stats_hdu=rho_hdu, tau_stats_hdu=tau_hdu, n_bins=1 + ) + with fits.open(out) as hdul: + cov = hdul["COVMAT"].data + # tau_0 block rows 2*N_ANG..3*N_ANG, tau_2 block 3*N_ANG..4*N_ANG. + cross = cov[2 * N_ANG : 3 * N_ANG, 3 * N_ANG : 4 * N_ANG] + expected = inp["tau_cov_full"][:N_ANG, N_ANG : 2 * N_ANG] + assert np.array_equal(cross, expected) + assert np.any(cross != 0.0) + + +def test_perturbed_xi_changes_output(tmp_path): + """Teeth: a changed xi+ input must change the converted data vector.""" + inp = _inputs(seed=0) + s = _sacc(inp) + out = tmp_path / "base.fits" + twopoint_convert.sacc_to_twopoint_fits(s, str(out), n_bins=1) + with fits.open(out) as hdul: + base_xip = hdul["XI_PLUS"].data["VALUE"].copy() + + inp2 = _inputs(seed=0) + inp2["xip"] = inp2["xip"] + 1.0 + s2 = _sacc(inp2) + out2 = tmp_path / "perturbed.fits" + twopoint_convert.sacc_to_twopoint_fits(s2, str(out2), n_bins=1) + with fits.open(out2) as hdul: + new_xip = hdul["XI_PLUS"].data["VALUE"] + + assert not np.array_equal(base_xip, new_xip) + assert np.array_equal(new_xip, inp2["xip"]) + + +def test_rho_tau_sidecars_required_together(tmp_path): + """Supplying only one of the rho/tau sidecars is a loud error.""" + inp = _inputs(seed=20) + s = _sacc(inp, rho_tau=True) + rho_hdu, _tau_hdu = _sidecar_hdus(tmp_path, inp) + with pytest.raises(ValueError, match="together"): + twopoint_convert.sacc_to_twopoint_fits( + s, str(tmp_path / "x.fits"), rho_stats_hdu=rho_hdu, n_bins=1 + ) + + +# ============================================================================= +# Fail-fast guards and permutation teeth (adversarial-review hardening) +# ============================================================================= + + +def test_tomographic_sacc_raises(tmp_path): + """A multi-bin SACC fails fast instead of silently truncating to (0, 0). + + Review finding (HIGH): ``n_bins`` alone drove the NZDATA column count while + the data vector and covariance were read from bin ``(0, 0)`` only, so a + 2-bin SACC + ``n_bins=2`` emitted a plausible-looking FITS carrying 1/3 of + the data. Both the ``n_bins`` and the tracer-pair mismatch must raise. + """ + inp = _inputs(seed=30) + s = sacc_io.new_sacc({0: (inp["z"], inp["nz"]), 1: (inp["z"], inp["nz"])}) + for pair in [(0, 0), (0, 1), (1, 1)]: + sacc_io.add_xi(s, pair, inp["theta"], inp["xip"], inp["xim"], grid="coarse") + s.add_covariance(np.eye(len(s.mean))) + + with pytest.raises(ValueError, match="single-bin only"): + twopoint_convert.sacc_to_twopoint_fits(s, str(tmp_path / "x.fits"), n_bins=2) + with pytest.raises(ValueError, match="single-bin only"): + twopoint_convert.sacc_to_twopoint_fits(s, str(tmp_path / "x.fits"), n_bins=1) + assert not (tmp_path / "x.fits").exists() + + +def test_sacc_without_xi_raises(tmp_path): + """A SACC with no ξ± points raises instead of writing an empty data vector.""" + inp = _inputs(seed=31) + s = sacc_io.new_sacc({0: (inp["z"], inp["nz"])}) + sacc_io.add_pseudo_cl( + s, + (0, 0), + inp["ell"], + inp["cl_ee"], + inp["cl_bb"], + inp["cl_eb"], + window_ells=np.arange(2, 102), + window_weights=np.random.default_rng(9).uniform(0, 1, (100, N_ELL)), + ) + s.add_covariance(np.eye(len(s.mean))) + + with pytest.raises(ValueError, match="nothing to convert"): + twopoint_convert.sacc_to_twopoint_fits(s, str(tmp_path / "x.fits")) + assert not (tmp_path / "x.fits").exists() + + +def test_covmat_blocks_exact_gather_encoded_cov(tmp_path): + """Every COVMAT/COVMAT_CELL entry is the exact ``np.ix_`` gather of the SACC + covariance, pinned with a (row, col)-encoded matrix. + + Review finding (MEDIUM): for a single bin pair the ξ gather happens to be + the identity permutation, so the byte-compares alone could pass with a + transposed or block-swapped gather. Encoding ``C[i, j] = i*n + j`` (asymmetric, + every entry unique) makes any transposition, offset, or wrong block produce + detectably wrong values; the τ gather is genuinely non-identity (τ_0− sits + between τ_0+ and τ_2+ in insertion order). Expected layout per + ``covdat_to_fits``: block_diag(ξ type-major gather, joint [τ_0+; τ_2+] + gather), with COVMAT_CELL the CELL_EE gather in its own HDU. + """ + inp = _inputs(seed=32) + s = _sacc(inp, cl=True, rho_tau=True) + n = len(s.mean) + encoded = np.arange(n * n, dtype=float).reshape(n, n) + s.add_covariance(encoded, overwrite=True) + rho_hdu, tau_hdu = _sidecar_hdus(tmp_path, inp) + + out = tmp_path / "encoded.fits" + twopoint_convert.sacc_to_twopoint_fits( + s, str(out), rho_stats_hdu=rho_hdu, tau_stats_hdu=tau_hdu, n_bins=1 + ) + + pair = (SOURCE, SOURCE) + xi_idx = np.concatenate( + [s.indices(sacc_io.XI_PLUS, pair), s.indices(sacc_io.XI_MINUS, pair)] + ) + tau_idx = np.concatenate( + [ + s.indices(sacc_io.TAU_PLUS.format(k=0), (SOURCE, PSF)), + s.indices(sacc_io.TAU_PLUS.format(k=2), (SOURCE, PSF)), + ] + ) + expected = twopoint_convert._block_diag( + encoded[np.ix_(xi_idx, xi_idx)], encoded[np.ix_(tau_idx, tau_idx)] + ) + cell_idx = s.indices(sacc_io.CL_EE, pair) + + with fits.open(out) as hdul: + np.testing.assert_array_equal(hdul["COVMAT"].data, expected) + np.testing.assert_array_equal( + hdul["COVMAT_CELL"].data, encoded[np.ix_(cell_idx, cell_idx)] + ) diff --git a/src/sp_validation/tests/test_twopoint_convert_realdata.py b/src/sp_validation/tests/test_twopoint_convert_realdata.py new file mode 100644 index 00000000..9d2739a6 --- /dev/null +++ b/src/sp_validation/tests/test_twopoint_convert_realdata.py @@ -0,0 +1,268 @@ +"""Candide-local byte-compare of the converter against real 2pt-FITS products. + +Skipped unless a real product exists on disk (candide only; never committed). +It closes the loop end to end on real data: take a real CosmoSIS 2pt-FITS, build +an analysis SACC from its own contents, convert that SACC back to a 2pt-FITS, and +byte-compare. + +The reference is *not* the on-disk file directly. The committed on-disk products +were written by an older ``cosmosis_fitting.py`` (they carry a CELL_BB HDU and +order COVMAT before NZ_SOURCE); the converter reproduces the *current* script, +which the PR-4 migration will use to regenerate them. So the meaningful contract +-- "the converter equals the current writer" -- is tested by running the current +``cosmosis_fitting.py`` builders on the same real contents and byte-comparing the +converter against that. For transparency the test also records the direct diff +against the stale on-disk file, and asserts only that it differs *by whole HDUs* +(the extra CELL_BB), not in any shared block -- i.e. the drift is purely the +known HDU-set change, with no silent data corruption. + +Observed on 2026-07-10 for ``SP_v1.4.6_leak_corr`` and a ``glass_mock`` sibling: +converter == current-script byte for byte; converter vs on-disk differs only by +the CELL_BB HDU. +""" + +import importlib.util +from pathlib import Path + +import numpy as np +import pytest +from astropy.io import fits + +from sp_validation import sacc_io, twopoint_convert + +_DATA = Path("/automnt/n17data/cdaley/unions/code/sp_validation/cosmo_inference/data") +_REAL_FILES = { + "SP_v1.4.6_leak_corr": _DATA + / "SP_v1.4.6_leak_corr_A_minsep=1.0_maxsep=250.0_nbins=20_npatch=1" + / "cosmosis_SP_v1.4.6_leak_corr_A_minsep=1.0_maxsep=250.0_nbins=20_npatch=1.fits", + "glass_mock_00001": _DATA / "glass_mock_00001" / "cosmosis_glass_mock_00001.fits", +} + +_SCRIPT = ( + Path(__file__).resolve().parents[3] + / "cosmo_inference" + / "scripts" + / "cosmosis_fitting.py" +) + + +def _load_cf(): + spec = importlib.util.spec_from_file_location("cosmosis_fitting", _SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def _sacc_from_2pt_fits(hdul): + """Build an analysis SACC from a real CosmoSIS 2pt-FITS's own contents. + + Reads the NZ, ξ±, pseudo-Cℓ (EE/BB) and blocked covariance back out of the + product and lays them into the standard SACC layout — the inverse direction + the converter then undoes. τ± minus and Cℓ EB are not in the 2pt-FITS, so + they are stored as zeros with identity covariance sub-blocks (the converter + consumes only the ``+``/EE parts). + """ + z = hdul["NZ_SOURCE"].data["Z_MID"].astype(float) + nz = hdul["NZ_SOURCE"].data["BIN1"].astype(float) + theta = hdul["XI_PLUS"].data["ANG"].astype(float) + xip = hdul["XI_PLUS"].data["VALUE"].astype(float) + xim = hdul["XI_MINUS"].data["VALUE"].astype(float) + n = len(theta) + + ell = hdul["CELL_EE"].data["ANG"].astype(float) + cl_ee = hdul["CELL_EE"].data["VALUE"].astype(float) + cl_bb = hdul["CELL_BB"].data["VALUE"].astype(float) + n_ell = len(ell) + + covmat = hdul["COVMAT"].data.astype(float) + covmat_cell = hdul["COVMAT_CELL"].data.astype(float) + xi_cov = covmat[: 2 * n, : 2 * n] + tau_joint = covmat[2 * n : 4 * n, 2 * n : 4 * n] + + tau0p = hdul["TAU_0_PLUS"].data["VALUE"].astype(float) + tau2p = hdul["TAU_2_PLUS"].data["VALUE"].astype(float) + + s = sacc_io.new_sacc({0: (z, nz)}) + sacc_io.add_xi(s, (0, 0), theta, xip, xim, grid="coarse") + sacc_io.add_pseudo_cl( + s, + (0, 0), + ell, + cl_ee, + cl_bb, + np.zeros(n_ell), + window_ells=np.arange(2, 102), + window_weights=np.ones((100, n_ell)), + ) + sacc_io.add_tau(s, (0, 0), 0, theta, tau0p, np.zeros(n)) + sacc_io.add_tau(s, (0, 0), 2, theta, tau2p, np.zeros(n)) + + source, psf = sacc_io.source_name(0), sacc_io.PSF_TRACER + idx = { + "xi": np.concatenate( + [ + s.indices(sacc_io.XI_PLUS, (source, source)), + s.indices(sacc_io.XI_MINUS, (source, source)), + ] + ), + "ee": s.indices(sacc_io.CL_EE, (source, source)), + "bb": s.indices(sacc_io.CL_BB, (source, source)), + "eb": s.indices(sacc_io.CL_EB, (source, source)), + "t0p": s.indices(sacc_io.TAU_PLUS.format(k=0), (source, psf)), + "t0m": s.indices(sacc_io.TAU_MINUS.format(k=0), (source, psf)), + "t2p": s.indices(sacc_io.TAU_PLUS.format(k=2), (source, psf)), + "t2m": s.indices(sacc_io.TAU_MINUS.format(k=2), (source, psf)), + } + full = np.zeros((len(s.mean), len(s.mean))) + full[np.ix_(idx["xi"], idx["xi"])] = xi_cov + full[np.ix_(idx["ee"], idx["ee"])] = covmat_cell + tau_pp = np.concatenate([idx["t0p"], idx["t2p"]]) + full[np.ix_(tau_pp, tau_pp)] = tau_joint + for key in ("bb", "eb"): + full[np.ix_(idx[key], idx[key])] = np.eye(n_ell) + for key in ("t0m", "t2m"): + full[np.ix_(idx[key], idx[key])] = np.eye(n) + s.add_covariance(full) + + tau_sidecar = fits.BinTableHDU.from_columns( + fits.ColDefs( + [ + fits.Column(name="theta", format="D", array=theta), + fits.Column(name="tau_0_p", format="D", array=tau0p), + fits.Column(name="tau_2_p", format="D", array=tau2p), + ] + ) + ) + return s, hdul["RHO_STATS"].copy(), tau_sidecar + + +def _current_script_reference(cf, hdul, tmp_path): + """Reference: run the current cosmosis_fitting.py on the real file's contents.""" + z = hdul["NZ_SOURCE"].data["Z_MID"].astype(float) + nz = hdul["NZ_SOURCE"].data["BIN1"].astype(float) + theta = hdul["XI_PLUS"].data["ANG"].astype(float) + xip = hdul["XI_PLUS"].data["VALUE"].astype(float) + xim = hdul["XI_MINUS"].data["VALUE"].astype(float) + n = len(theta) + ell = hdul["CELL_EE"].data["ANG"].astype(float) + cl_ee = hdul["CELL_EE"].data["VALUE"].astype(float) + cl_bb = hdul["CELL_BB"].data["VALUE"].astype(float) + covmat = hdul["COVMAT"].data.astype(float) + covmat_cell = hdul["COVMAT_CELL"].data.astype(float) + tau0p = hdul["TAU_0_PLUS"].data["VALUE"].astype(float) + tau2p = hdul["TAU_2_PLUS"].data["VALUE"].astype(float) + + np.savetxt(tmp_path / "nz.txt", np.column_stack([z, nz])) + np.savetxt(tmp_path / "cov.txt", covmat[: 2 * n, : 2 * n]) + tau_cov = np.zeros((3 * n, 3 * n)) + tau_cov[: 2 * n, : 2 * n] = covmat[2 * n : 4 * n, 2 * n : 4 * n] + np.save(tmp_path / "cov_tau.npy", tau_cov) + cl_block = np.zeros((5, len(ell))) + cl_block[0], cl_block[1], cl_block[4] = ell, cl_ee, cl_bb + np.save(tmp_path / "cl.npy", cl_block) + np.save(tmp_path / "cl_cov.npy", covmat_cell) + + nz_hdu = cf.nz_to_fits(str(tmp_path / "nz.txt")) + xip_hdu = cf._create_2pt_hdu(xip, theta, "XI_PLUS", "G+R", "G+R") + xim_hdu = cf._create_2pt_hdu(xim, theta, "XI_MINUS", "G-R", "G-R") + cov_hdu = cf.covdat_to_fits( + str(tmp_path / "cov.txt"), filename_cov_tau=str(tmp_path / "cov_tau.npy") + ) + ell_r, cl_ee_r, cl_bb_r = cf.load_pseudo_cl(str(tmp_path / "cl.npy")) + cl_ee_hdu, _ = cf.cl_to_fits(ell_r, cl_ee_r, cl_bb_r) + cov_cl_hdu = cf.cov_cl_to_fits(str(tmp_path / "cl_cov.npy"), cov_hdu="COVAR_FULL") + + fits.HDUList([fits.PrimaryHDU(), hdul["RHO_STATS"].copy()]).writeto( + tmp_path / "rho.fits", overwrite=True + ) + rho_hdu = cf.rho_to_fits(str(tmp_path / "rho.fits"), theta=theta) + tau_sidecar = fits.BinTableHDU.from_columns( + fits.ColDefs( + [ + fits.Column(name="theta", format="D", array=theta), + fits.Column(name="tau_0_p", format="D", array=tau0p), + fits.Column(name="tau_2_p", format="D", array=tau2p), + ] + ) + ) + fits.HDUList([fits.PrimaryHDU(), tau_sidecar]).writeto( + tmp_path / "tau.fits", overwrite=True + ) + tau0_hdu, tau2_hdu = cf.tau_to_fits(str(tmp_path / "tau.fits"), theta=theta) + + out = tmp_path / "reference.fits" + fits.HDUList( + [ + fits.PrimaryHDU(), + nz_hdu, + cov_hdu, + cov_cl_hdu, + xip_hdu, + xim_hdu, + cl_ee_hdu, + tau0_hdu, + tau2_hdu, + rho_hdu, + ] + ).writeto(out, overwrite=True) + return out + + +@pytest.mark.parametrize("label", list(_REAL_FILES)) +def test_realdata_roundtrip_byte_equal(label, tmp_path): + """Converter reproduces the current writer byte for byte on a real product.""" + real = _REAL_FILES[label] + if not real.exists(): + pytest.skip(f"real 2pt-FITS not on disk: {real}") + cf = _load_cf() + + with fits.open(real) as hdul: + s, rho_hdu, tau_hdu = _sacc_from_2pt_fits(hdul) + reference = _current_script_reference(cf, hdul, tmp_path) + + converted = tmp_path / "converted.fits" + twopoint_convert.sacc_to_twopoint_fits( + s, str(converted), rho_stats_hdu=rho_hdu, tau_stats_hdu=tau_hdu, n_bins=1 + ) + + # The contract: converter == current cosmosis_fitting.py, byte for byte. + assert converted.read_bytes() == reference.read_bytes() + + +@pytest.mark.parametrize("label", list(_REAL_FILES)) +def test_realdata_ondisk_drift_is_only_cell_bb(label, tmp_path): + """The stale on-disk file differs from the converter *only* by the CELL_BB HDU. + + Documents (and guards) the known script-version drift: the on-disk products + were written before CELL_BB was dropped from the assembly, so they carry one + extra HDU. Every HDU the two share carries the same data to floating-point + precision — the drift is a whole-HDU addition, never a silent change to a + shared block. (Bin-edge columns differ by ~1e-17 float noise: the stale file + stored a clean ``Z_LOW=0``, while ``z_mid - step/2`` rounds to ``-1.7e-18``; + both are the same number, so the shared-data check is ``allclose``, not + bitwise — bitwise equality is asserted against the *current* writer above.) + """ + real = _REAL_FILES[label] + if not real.exists(): + pytest.skip(f"real 2pt-FITS not on disk: {real}") + + with fits.open(real) as hdul: + s, rho_hdu, tau_hdu = _sacc_from_2pt_fits(hdul) + ondisk_names = [h.name for h in hdul] + converted = tmp_path / "converted.fits" + twopoint_convert.sacc_to_twopoint_fits( + s, str(converted), rho_stats_hdu=rho_hdu, tau_stats_hdu=tau_hdu, n_bins=1 + ) + with fits.open(converted) as conv: + conv_names = [h.name for h in conv] + # The only HDU the on-disk file has that the converter does not. + assert set(ondisk_names) - set(conv_names) == {"CELL_BB"} + # Every shared table HDU carries the same data to float precision. + for name in conv_names: + if name == "PRIMARY": + continue + a, b = hdul[name].data, conv[name].data + if hasattr(a, "names"): + assert all(np.allclose(a[c], b[c]) for c in a.names), name + else: + assert np.allclose(a, b), name diff --git a/src/sp_validation/twopoint_convert.py b/src/sp_validation/twopoint_convert.py new file mode 100644 index 00000000..8f787f3b --- /dev/null +++ b/src/sp_validation/twopoint_convert.py @@ -0,0 +1,368 @@ +"""TWOPOINT_CONVERT. + +:Name: twopoint_convert.py + +:Description: Convert an analysis SACC file into the "2pt FITS" that CosmoSIS's + ``2pt_like`` (and Sacha Guerrini's ρ/τ ``2pt_like_xi_sys`` fork) + reads. The output reproduces today's hand-assembled product from + ``cosmo_inference/scripts/cosmosis_fitting.py`` HDU-for-HDU and + byte-for-byte: an NZDATA table, XI_PLUS / XI_MINUS 2pt tables, + optional CELL_EE / CELL_BB pseudo-Cℓ tables, the blocked COVMAT + (with ``STRT_i`` block-offset headers) and separate COVMAT_CELL, + and — when the SACC carries them — the TAU_{0,2}_PLUS 2pt tables + and the RHO_STATS table. + + The converter is the *inverse* of the SACC writers in + :mod:`sp_validation.sacc_io`: it reads statistics back through + those readers and lays them into the DES ``twopoint`` FITS + convention. SACC's canonical order is pair-major (per pair + ``[ξ+; ξ−]``); the 2pt-FITS layout is type-major (all ξ+, then all + ξ−), so the data-vector and its covariance are permuted here via + ``s.indices`` rather than assuming any global order. + + Scope note (single-bin today, tomography-ready): the assembly this + mirrors is single-tomographic-bin — BIN1/BIN2 are all 1, one NZ + ``BIN1`` column. The converter reads bin ``(0, 0)`` accordingly. + A tomographic 2pt-FITS layout (multiple bin pairs, per-pair + BIN1/BIN2, one NZ column per bin) is a later extension; it is not + what today's CosmoSIS pipeline consumes, so it is out of scope for + the byte-compatible converter. + + Rho/tau caveat: the SACC layout stores ρ±/τ± *values* only, while + the 2pt-FITS RHO_STATS table also carries the per-mode *variances* + (``varrho_*``) that Sacha's fork's covariance path reads. Those + variances are not recoverable from the analysis SACC. The + converter therefore writes RHO_STATS / TAU HDUs only when a + ``rho_stats``/``tau_stats`` sidecar FITS is supplied (the same + file today's assembly copies verbatim); it never fabricates + variances. ξ±, Cℓ, n(z) and the covariance — the data vector + CosmoSIS fits — are fully reconstructed from SACC alone. +""" + +import numpy as np +from astropy.io import fits + +from . import sacc_io + +# The QUANT1/QUANT2 header pair CosmoSIS stamps on each 2pt table, keyed by the +# extension name — copied from cosmosis_fitting.py so the headers match card +# for card. +_QUANT = { + "XI_PLUS": ("G+R", "G+R"), + "XI_MINUS": ("G-R", "G-R"), + "CELL_EE": ("GEF", "GEF"), + "CELL_BB": ("GBF", "GBF"), + "TAU_0_PLUS": ("G+R", "P+R"), + "TAU_2_PLUS": ("G+R", "SR+R"), +} + + +def _twopoint_hdu(name, values, ang, *, ang_unit=None): + """Build one 2pt BinTableHDU (BIN1/BIN2/ANGBIN/VALUE/ANG). + + Reproduces ``cosmosis_fitting.py._create_2pt_hdu`` /``cl_to_fits`` exactly: + same column order and formats, the ``2PTDATA`` marker, the QUANT pair for + ``name``, and NZ_SOURCE kernels. ``ang_unit`` stamps ``TUNIT`` on the ANG + column ("arcmin" for real-space ξ/τ; unset for Cℓ, whose ANG is ℓ). + """ + nbins = len(values) + angbin = np.arange(1, nbins + 1) + columns = [ + fits.Column(name="BIN1", format="K", array=np.ones(nbins)), + fits.Column(name="BIN2", format="K", array=np.ones(nbins)), + fits.Column(name="ANGBIN", format="K", array=angbin), + fits.Column(name="VALUE", format="D", array=values), + fits.Column(name="ANG", format="D", unit=ang_unit, array=ang), + ] + hdu = fits.BinTableHDU.from_columns(fits.ColDefs(columns), name=name) + quant1, quant2 = _QUANT[name] + for key, value in { + "2PTDATA": "T", + "QUANT1": quant1, + "QUANT2": quant2, + "KERNEL_1": "NZ_SOURCE", + "KERNEL_2": "NZ_SOURCE", + "WINDOWS": "SAMPLE", + }.items(): + hdu.header[key] = value + return hdu + + +def _nz_hdu(s, n_bins): + """Build the NZDATA HDU from the SACC ``source_i`` NZ tracers. + + Reproduces ``cosmosis_fitting.py.nz_to_fits``: Z_MID from the tracer ``z`` + grid (assumed uniform), Z_LOW/Z_HIGH as ± half a step, one ``BIN{i+1}`` + column per source bin, and the NZDATA/NBIN/NZ header cards. All source bins + are required to share the ``z`` grid — the single ``Z_MID`` axis of the + DES NZDATA table. + """ + z_mid, nz0 = sacc_io.get_nz(s, 0) + z_mid = np.asarray(z_mid, dtype=float) + step = z_mid[1] - z_mid[0] + z_low = z_mid - step / 2 + z_high = z_mid + step / 2 + + columns = [ + fits.Column(name="Z_LOW", format="D", array=z_low), + fits.Column(name="Z_MID", format="D", array=z_mid), + fits.Column(name="Z_HIGH", format="D", array=z_high), + ] + for i in range(n_bins): + z_i, nz_i = sacc_io.get_nz(s, i) + if not np.array_equal(np.asarray(z_i, dtype=float), z_mid): + raise ValueError( + f"source bin {i} n(z) grid differs from source bin 0; the DES " + "NZDATA table requires one shared Z_MID axis" + ) + columns.append(fits.Column(name=f"BIN{i + 1}", format="D", array=nz_i)) + + hdu = fits.BinTableHDU.from_columns(fits.ColDefs(columns), name="NZDATA") + for key, value in { + "NZDATA": "T ", + "EXTNAME": "NZ_SOURCE", + "NBIN": n_bins, + "NZ": len(z_low), + }.items(): + hdu.header[key] = value + return hdu + + +def _cov_hdu(matrix, block_names, block_starts, extname="COVMAT", name_in_ctor=False): + """Build a covariance ImageHDU with ``NAME_i``/``STRT_i`` block headers. + + Reproduces the two covariance builders in ``cosmosis_fitting.py`` card for + card. The blocked ξ/τ ``covdat_to_fits`` builds ``ImageHDU(cov)`` unnamed + and stamps ``COVDATA`` then ``EXTNAME`` from a dict; the ``cov_cl_to_fits`` + CELL covariance builds ``ImageHDU(cov, name="COVMAT_CELL")`` (so the EXTNAME + card is created early, with astropy's standard comment) before re-stamping. + ``name_in_ctor`` selects the second form so the card order matches exactly. + """ + matrix = np.asarray(matrix, dtype=np.float64) + if matrix.shape[0] != matrix.shape[1]: + raise ValueError(f"covariance must be square; got shape {matrix.shape}") + hdu = fits.ImageHDU(matrix, name=extname) if name_in_ctor else fits.ImageHDU(matrix) + hdu.header["COVDATA"] = "True" + hdu.header["EXTNAME"] = extname + for i, (name, start) in enumerate(zip(block_names, block_starts)): + hdu.header[f"NAME_{i}"] = name + hdu.header[f"STRT_{i}"] = int(start) + return hdu + + +def _type_major_xi(s, bins): + """Return ``(theta, xip, xim)`` for one bin pair from the SACC coarse grid. + + ``sacc_io.get_xi`` already returns each statistic in insertion (= ascending + θ) order; the type-major split (all ξ+, then all ξ−) is exactly the two + arrays it hands back, so no further permutation is needed for a single pair. + """ + return sacc_io.get_xi(s, bins, grid="coarse") + + +def _require_single_bin(s, n_bins): + """Fail fast unless the SACC is a valid single-bin ξ product. + + The converter emits the single-bin 2pt-FITS today's CosmoSIS pipeline reads + (BIN1/BIN2 all 1, one NZ column). A tomographic SACC would otherwise slip + through silently — ``n_bins`` alone drives the NZDATA column count while the + ξ/covariance are read from bin ``(0, 0)`` only, so a 2-bin file would emit a + ``NBIN=2`` n(z) beside a data vector holding just the ``(0, 0)`` pair. + Guards both the empty-ξ case and the single-bin contract; tomographic + emission lands with the tomographic round. + """ + pairs = s.get_tracer_combinations(sacc_io.XI_PLUS) + if not pairs: + raise ValueError( + f"SACC has no {sacc_io.XI_PLUS} points — nothing to convert; the " + "2pt-FITS data vector is built from the ξ± statistics" + ) + expected = (sacc_io.source_name(0), sacc_io.source_name(0)) + if n_bins != 1 or set(pairs) != {expected}: + raise ValueError( + f"converter is single-bin only (n_bins=1, ξ pairs == {{{expected}}}); " + f"got n_bins={n_bins} and ξ pairs {sorted(pairs)}. Tomographic " + "emission (multiple bin pairs, per-pair BIN1/BIN2, one NZ column per " + "bin) lands with the tomographic round." + ) + + +def sacc_to_twopoint_fits( + s, + path, + *, + rho_stats_hdu=None, + tau_stats_hdu=None, + n_bins=1, +): + """Convert an analysis SACC to a CosmoSIS 2pt-FITS file. + + The assembled ``HDUList`` matches today's ``cosmosis_fitting.py`` product + for the configuration the SACC describes: PRIMARY, NZ_SOURCE, COVMAT, then + (if present) COVMAT_CELL, XI_PLUS, XI_MINUS, (if present) CELL_EE / CELL_BB, + and (if the rho/tau sidecars are supplied) TAU_0_PLUS, TAU_2_PLUS, + RHO_STATS. The data vector and its covariance are laid out type-major + (all ξ+, then all ξ−, then the τ blocks), which is the DES ``twopoint`` + convention CosmoSIS reads. + + Parameters + ---------- + s : sacc.Sacc + Analysis SACC (coarse ξ±, optional pseudo-Cℓ, covariance, and — for the + ρ/τ product — the τ data points; see ``rho_stats_hdu``). + path : str + Output FITS path (overwritten). + rho_stats_hdu, tau_stats_hdu : astropy.io.fits.BinTableHDU, optional + The rho-stats / tau-stats sidecar HDUs, copied verbatim as today's + assembly does. Required together to write the ρ/τ product; the SACC + alone cannot rebuild the ``varrho_*`` columns Sacha's fork reads. When + omitted, a pure ξ (± Cℓ) product is written. + n_bins : int, optional + Number of source tomographic bins. Must be ``1``: this converter emits + the single-bin 2pt-FITS today's CosmoSIS pipeline consumes. Tomographic + emission (multiple bin pairs, per-pair BIN1/BIN2, one NZ column per bin) + lands with the tomographic round; the converter fails fast on anything + else rather than silently truncating to bin ``(0, 0)``. + + Returns + ------- + astropy.io.fits.HDUList + The assembled list, also written to ``path``. + + Raises + ------ + ValueError + If the SACC has no ξ points; if ``n_bins != 1`` or the SACC's ξ tracer + pairs are anything other than exactly ``{(source_0, source_0)}`` (the + single-bin contract); or if exactly one of the ρ/τ sidecars is supplied. + """ + if (rho_stats_hdu is None) != (tau_stats_hdu is None): + raise ValueError( + "rho_stats_hdu and tau_stats_hdu must be supplied together " + "(the ρ/τ product needs both, or neither for a pure-ξ product)" + ) + _require_single_bin(s, n_bins) + use_rho_tau = rho_stats_hdu is not None + bins = (0, 0) + + nz_hdu = _nz_hdu(s, n_bins) + theta, xip, xim = _type_major_xi(s, bins) + xip_hdu = _twopoint_hdu("XI_PLUS", xip, theta, ang_unit="arcmin") + xim_hdu = _twopoint_hdu("XI_MINUS", xim, theta, ang_unit="arcmin") + + cell_hdu, cov_cell_hdu = _build_cell(s, bins) + + cov_hdu = _build_covmat(s, bins, use_rho_tau=use_rho_tau) + + tau_hdus, rho_hdu = _build_rho_tau(rho_stats_hdu, tau_stats_hdu, theta, use_rho_tau) + + # HDU order mirrors cosmosis_fitting.py's __main__: PRIMARY, NZ, COVMAT, + # COVMAT_CELL, XI±, CELL_EE, then the τ/ρ tables. + hdu_list = [fits.PrimaryHDU(), nz_hdu, cov_hdu] + if cov_cell_hdu is not None: + hdu_list.append(cov_cell_hdu) + hdu_list.extend([xip_hdu, xim_hdu]) + if cell_hdu is not None: + hdu_list.append(cell_hdu) + if use_rho_tau: + hdu_list.extend([*tau_hdus, rho_hdu]) + + hdul = fits.HDUList(hdu_list) + hdul.writeto(path, overwrite=True) + return hdul + + +def _build_cell(s, bins): + """Build the CELL_EE 2pt HDU plus the COVMAT_CELL HDU from the SACC pseudo-Cℓ. + + Returns ``(None, None)`` when the SACC has no pseudo-Cℓ. Only CELL_EE is + emitted — the harmonic ``2pt_like`` fits ``data_sets=CELL_EE``, and today's + assembly appends CELL_EE alone (it builds a CELL_BB HDU but discards it). + The SACC still carries EE/BB/EB with bandpower windows for the B-mode + null-test path; this converter surfaces only the block CosmoSIS reads. The + CELL covariance (the EE bandpower covariance) lives in its own COVMAT_CELL + ImageHDU, matching today's product. + """ + if sacc_io.CL_EE not in s.get_data_types(): + return None, None + + ell, cl_ee, _cl_bb, _cl_eb, _window = sacc_io.get_pseudo_cl(s, bins) + cell_hdu = _twopoint_hdu("CELL_EE", cl_ee, ell) + cell_idx = s.indices(sacc_io.CL_EE, sacc_io._pair(bins)) + cov_cell = s.covariance.dense[np.ix_(cell_idx, cell_idx)] + cov_cell_hdu = _cov_hdu( + cov_cell, ["CELL_EE"], [0], extname="COVMAT_CELL", name_in_ctor=True + ) + return cell_hdu, cov_cell_hdu + + +def _build_covmat(s, bins, *, use_rho_tau): + """Assemble the blocked COVMAT (ξ± type-major, then the τ blocks). + + The ξ covariance is pulled from the SACC as the contiguous ξ+/ξ− block for + the pair and permuted from pair-major (SACC) to type-major (2pt-FITS). Under + ``use_rho_tau`` the τ_0/τ_2 covariance blocks are appended block-diagonally + with zero ξ↔τ cross-blocks, exactly as ``covdat_to_fits`` builds them. + """ + pair = sacc_io._pair(bins) + idx_p = s.indices(sacc_io.XI_PLUS, pair) + idx_m = s.indices(sacc_io.XI_MINUS, pair) + n_theta = len(idx_p) + xi_idx = np.concatenate([idx_p, idx_m]) # type-major permutation + xi_cov = s.covariance.dense[np.ix_(xi_idx, xi_idx)] + + names = ["XI_PLUS", "XI_MINUS"] + starts = [0, n_theta] + matrix = xi_cov + + if use_rho_tau: + # The τ covariance couples τ_0+ and τ_2+ (today's assembly truncates the + # 3-statistic CosmoCov τ covariance to its first 2 blocks and lays it in + # as ONE contiguous [τ_0+; τ_2+] block — cross-correlation kept). In the + # SACC those two selections are not adjacent (τ_0− sits between them), so + # gather both index sets and extract the joint sub-block, ξ↔τ zero. + tau_pair = (sacc_io.source_name(0), sacc_io.PSF_TRACER) + idx_tau0 = s.indices(sacc_io.TAU_PLUS.format(k=0), tau_pair) + idx_tau2 = s.indices(sacc_io.TAU_PLUS.format(k=2), tau_pair) + tau_idx = np.concatenate([idx_tau0, idx_tau2]) + tau_cov = s.covariance.dense[np.ix_(tau_idx, tau_idx)] + matrix = _block_diag(matrix, tau_cov) + names += ["TAU_0_PLUS", "TAU_2_PLUS"] + starts += [2 * n_theta, 2 * n_theta + len(idx_tau0)] + + return _cov_hdu(matrix, names, starts) + + +def _block_diag(*blocks): + """Stack square blocks block-diagonally with zero cross-blocks.""" + sizes = [b.shape[0] for b in blocks] + n = sum(sizes) + out = np.zeros((n, n)) + start = 0 + for b in blocks: + out[start : start + b.shape[0], start : start + b.shape[0]] = b + start += b.shape[0] + return out + + +def _build_rho_tau(rho_stats_hdu, tau_stats_hdu, theta, use_rho_tau): + """Build the TAU_{0,2}_PLUS 2pt HDUs and the verbatim RHO_STATS HDU. + + Mirrors ``tau_to_fits`` / ``rho_to_fits``: τ_0/τ_2 read their ``tau_k_p`` + columns onto the shared ξ θ grid (consistency step); RHO_STATS is copied + verbatim from the sidecar with its θ column forced onto the ξ grid. The + ``varrho_*`` columns ride along in the copy — they are why the sidecar is + required (the SACC cannot supply them). + """ + if not use_rho_tau: + return (), None + + tau = tau_stats_hdu.data + tau0_hdu = _twopoint_hdu("TAU_0_PLUS", tau["tau_0_p"], theta, ang_unit="arcmin") + tau2_hdu = _twopoint_hdu("TAU_2_PLUS", tau["tau_2_p"], theta, ang_unit="arcmin") + + rho_hdu = rho_stats_hdu.copy() + rho_hdu.name = "RHO_STATS" + rho_hdu.data = rho_hdu.data.copy() + rho_hdu.data["theta"] = theta + return (tau0_hdu, tau2_hdu), rho_hdu diff --git a/workflow/rules/inference.smk b/workflow/rules/inference.smk index 6bf68645..ab7b456c 100644 --- a/workflow/rules/inference.smk +++ b/workflow/rules/inference.smk @@ -1,15 +1,19 @@ -# Imports from Snakefile: FIDUCIAL, COSMO_INFERENCE, COSMO_VAL, covariance_path, build_redshift_path, fiducial_binning_suffix -# NOTE: dormant subsystem. The file-name plumbing (config-driven paths + the -# producer-tagged pseudo-Cl names) is fixed and the DAG is valid, but it has not -# been run end-to-end. Reviving it still needs the FITS-CONTENT plumbing -# reconciled: cosmosis_fitting.py reads ELL/EE/BB + COVAR_FULL, while the -# producers write PSEUDO_CELL/ELL + COVAR_BB_BB. +# Imports from common (via `from common import *`): FIDUCIAL, COSMO_INFERENCE, +# COSMO_VAL, WORKFLOW_SCRIPTS, covariance_path, build_redshift_path, +# fiducial_binning_suffix. cv_analysis_sacc arrives from cosmo_val.smk (resolved +# lazily at DAG time, since that file is included after this one). +# +# Two paths live here: +# * Real-data inference_prep — LIVE (PR 7): consumes the assembled {version}.sacc +# and emits the converter 2pt-FITS + both engine inis (2pt_like, sacc_like). +# * glass-mock rules — still cosmosis_fitting.py-based (their SACC migration is +# out of scope); the pseudo-Cl file-name plumbing they depend on stays below. # Output root for CosmoSIS data products + configs. COSMO_INFERENCE (common.py) # already resolves to THIS repo's cosmo_inference dir, so the products land # beside the code that builds them rather than in a contributor's home. COSMO_INFERENCE_PROD = COSMO_INFERENCE -# Working directory for the cosmosis_fitting.py invocation — the same repo dir. +# Working directory for the (glass-mock) cosmosis_fitting.py invocation. COSMO_INFERENCE_RUNDIR = str(COSMO_INFERENCE) # External chain/mock locations are deployment-specific, so they live in config. @@ -56,88 +60,115 @@ def pseudo_cl_assets(version): return str(cl_path), str(cov_path) # --------------------------------------------------------------------------- -# DORMANT — pre-SACC cosmosis assembly. Migration to native SACC deferred to -# PR 7 (native-SACC inference consumption); do NOT deep-migrate here. +# Real-data inference prep — LIVE (native SACC, PR 7). Consumes the assembled +# analysis {version}.sacc (cosmo_val.smk's assemble_sacc rule) and emits the two +# file-prep products the A_ia (IA-only, ξ±) fiducial pipeline needs: +# (a) the converter 2pt-FITS (sacc_to_twopoint_fits) + a generated 2pt_like ini +# — the validating/legacy path (retiring cosmosis_fitting.py's assembly), +# (b) a generated sacc_like ini pointing at the SACC directly — the native path +# validated bit-for-bit against (a) (test_sacc_like.py). +# The converter is A_ia-scoped: no rho/tau sidecars, so it emits a pure-ξ FITS +# (it ignores the SACC's extra data types). This is file-prep only — the actual +# CosmoSIS sampling still runs via pipeline.sh against these products. # -# The SACC migration (PR 4) removed the data products several of these inputs -# name, so this rule's DAG no longer resolves and is NOT reachable from the -# cosmo_val suite (cosmo_val_all never requests it). Stale inputs: -# - xi_plus / xi_minus FITS: the `xi` rule now emits the coarse ξ± SACC part -# ({version}_xi_coarse_...sacc), not per-sign FITS. -# - pseudo_cl / pseudo_cl_cov via pseudo_cl_assets(): the `pseudo_cl` rule now -# writes .sacc (pseudo_cl_assets still requests .fits). -# PR 7 rewires this to consume the assembled {version}.sacc (built by -# cosmo_val.smk's assemble_sacc rule) directly, retiring cosmosis_fitting.py's -# per-product FITS assembly. Until then the inference target is knowingly red. +# The glass-mock rules below stay cosmosis_fitting.py-based; their SACC migration +# is out of scope for PR 7. # --------------------------------------------------------------------------- +# Generated per-version configs land in the (env-overridable) output root. +INFERENCE_CONFIG_OUT = COSMO_INFERENCE_PROD / "cosmosis_config" +# The ini TEMPLATES are source files: anchor them on the running checkout (repo +# root = the workflow dir's parent, via WORKFLOW_SCRIPTS), NOT on the output root +# — so a template edit in this checkout drives the DAG even when COSMO_INFERENCE +# points elsewhere. In a normal (non-worktree) run the two roots coincide. +INFERENCE_TEMPLATE_DIR = ( + Path(os.path.dirname(WORKFLOW_SCRIPTS)).parent / "cosmo_inference" / "cosmosis_config" +) + + +def _csl_dir(): + """The CSL checkout that fills COSMOSIS_DIR / sacc_like csl_dir in the inis. + + Read lazily (at DAG time, inside inference_prep's params) rather than at + module parse time: inference.smk is included by every paper workflow, but + only papers that run inference (cosmo_val) carry inference.csl_dir. A missing + key still fails loudly — just when the real-data inference is actually built, + not when an unrelated (bmodes) workflow merely parses this file. + """ + return INFERENCE["csl_dir"] + + rule inference_prep: input: - # Processed covariance matrix - use centralized covariance_path() - cov_matrix=lambda w: covariance_path(w.version, w.blind, min_sep=w.min_sep, max_sep=w.max_sep, nbins=w.nbins), - # Xi FITS files — PRE-SACC (no longer produced; see dormant note above) - xi_plus=str(COSMO_VAL / "xi_plus_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - xi_minus=str(COSMO_VAL / "xi_minus_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - # n(z) file (using new location with base version mapping) - nz_file=lambda w: build_redshift_path(w.version, w.blind), - # rho/tau stats - rho_stats=str(COSMO_VAL / "rho_tau_stats/rho_stats_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - tau_stats=str(COSMO_VAL / "rho_tau_stats/tau_stats_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), - # tau covariance (tracked as dependency) - tau_cov=str(COSMO_VAL / "rho_tau_stats/cov_tau_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}_th.npy"), - # pseudo_cl / pseudo_cl_cov — PRE-SACC (.fits path; producer now writes .sacc) - pseudo_cl=lambda w: pseudo_cl_assets(w.version)[0], - pseudo_cl_cov=lambda w: pseudo_cl_assets(w.version)[1], + # The terminal assembled analysis SACC (cosmo_val.smk assemble_sacc). Bound + # lazily through its helper so the filename tracks that rule, not a literal. + sacc=lambda w: cv_analysis_sacc(w.version), + # The two pipeline ini templates are static repo files, but binding them as + # inputs (not params) puts them in the DAG, so editing a template + # regenerates the configs rather than leaving stale output on disk. + template_2pt=str(INFERENCE_TEMPLATE_DIR / "cosmosis_pipeline_A_ia.ini"), + template_sacc=str(INFERENCE_TEMPLATE_DIR / "cosmosis_pipeline_A_ia_sacc.ini"), output: - fits_file=str( - COSMO_INFERENCE_PROD - / "data/{version}_{blind}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}/cosmosis_{version}_{blind}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits" + fits_file=str(COSMO_INFERENCE_PROD / "data/{version}/cosmosis_{version}.fits"), + config_file_2pt=str( + INFERENCE_CONFIG_OUT / "cosmosis_pipeline_{version}_A_ia.ini" + ), + config_file_sacc=str( + INFERENCE_CONFIG_OUT / "cosmosis_pipeline_{version}_A_ia_sacc.ini" ), - config_file=str( - COSMO_INFERENCE_PROD - / "cosmosis_config/cosmosis_pipeline_{version}_{blind}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.ini" - ) params: - cosmosis_root="{version}_{blind}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}", - data_dir=f"{CHAINS_DIR}/{{version}}_{{blind}}_minsep={{min_sep}}_maxsep={{max_sep}}_nbins={{nbins}}_npatch={{npatch}}", - output_root=str(COSMO_INFERENCE_PROD), + # SCRATCH = the per-version chain output root the generated inis point at. + scratch=lambda w: f"{CHAINS_DIR}/{w.version}", + cosmosis_dir=lambda w: _csl_dir(), threads: 1 resources: mem_mb=8000, runtime=10, - shell: - """ - cd {COSMO_INFERENCE_RUNDIR} + run: + import os + import sys - # Run inference preparation step with cosmosis_fitting.py - python scripts/cosmosis_fitting.py \ - --cosmosis-root {params.cosmosis_root} \ - --nz-file {input.nz_file} \ - --data-dir {params.data_dir} \ - --output-root {params.output_root} \ - --xi {input.xi_plus} {input.xi_minus} \ - --cov-xi {input.cov_matrix} \ - --use-rho-tau \ - --rho-stats {input.rho_stats} \ - --tau-stats {input.tau_stats} \ - --cov-tau {input.tau_cov} \ - --cl-file {input.pseudo_cl} \ - --cov-cl {input.pseudo_cl_cov} - """ + from sp_validation import sacc_io + from sp_validation.twopoint_convert import sacc_to_twopoint_fits + + os.makedirs(os.path.dirname(output.fits_file), exist_ok=True) + + # (a) converter 2pt-FITS — pure ξ (A_ia scope; no rho/tau sidecars). + sacc_to_twopoint_fits(sacc_io.load(input.sacc), output.fits_file, n_bins=1) + + # (b) + (c) the two generated pipeline inis, from the template inputs. + # WORKFLOW_SCRIPTS (common.py) is the absolute generic-workflow scripts dir. + sys.path.insert(0, WORKFLOW_SCRIPTS) + from generate_inference_config import ( + _substitutions, + generate_inference_config, + ) + + generate_inference_config( + input.template_2pt, + output.config_file_2pt, + _substitutions( + scratch=params.scratch, + cosmosis_dir=params.cosmosis_dir, + fits_file=output.fits_file, + ), + ) + generate_inference_config( + input.template_sacc, + output.config_file_sacc, + _substitutions( + scratch=params.scratch, + cosmosis_dir=params.cosmosis_dir, + sacc_file=input.sacc, + ), + ) rule inference_fiducial: input: - # Use the same output patterns as inference_prep with FIDUCIAL params - rules.inference_prep.output.fits_file.format( - version=FIDUCIAL["version"], blind=FIDUCIAL["blind"], - min_sep=FIDUCIAL["min_sep"], max_sep=FIDUCIAL["max_sep"], - nbins=FIDUCIAL["nbins"], npatch=FIDUCIAL["npatch"] - ), - rules.inference_prep.output.config_file.format( - version=FIDUCIAL["version"], blind=FIDUCIAL["blind"], - min_sep=FIDUCIAL["min_sep"], max_sep=FIDUCIAL["max_sep"], - nbins=FIDUCIAL["nbins"], npatch=FIDUCIAL["npatch"] - ) + # The fiducial version's prep products (both engine inis + the FITS). + rules.inference_prep.output.fits_file.format(version=FIDUCIAL["version"]), + rules.inference_prep.output.config_file_2pt.format(version=FIDUCIAL["version"]), + rules.inference_prep.output.config_file_sacc.format(version=FIDUCIAL["version"]), rule inference_glass_mocks: diff --git a/workflow/scripts/generate_inference_config.py b/workflow/scripts/generate_inference_config.py new file mode 100644 index 00000000..68a8efe7 --- /dev/null +++ b/workflow/scripts/generate_inference_config.py @@ -0,0 +1,159 @@ +"""Generate a CosmoSIS pipeline ini by filling a template's ``[DEFAULT]`` section. + +Dual-mode, like ``assemble_sacc.py``. Under Snakemake (``script:`` directive) the +injected ``snakemake`` object supplies the template, output path and DEFAULT +substitutions; as a standalone CLI (argparse) the same fill runs from explicit +flags. + +The template carries ``%(KEY)s`` interpolation placeholders (SCRATCH, FITS_FILE +or SACC_FILE, COSMOSIS_DIR, SP_VALIDATION_MODULES) in its module sections; this +script prepends the concrete ``KEY = value`` lines into ``[DEFAULT]`` so +CosmoSIS's ConfigParser resolves them at load. It is deliberately plain text +processing — appending lines after the ``[DEFAULT]`` header, the same idiom as +``pipeline.sh``'s ``sed -i "/^\\[DEFAULT\\]/a\\KEY = value"`` — rather than a +configparser round-trip, which would strip the template's comments and its +``%(...)s`` interpolation. + +``SP_VALIDATION_MODULES`` is resolved from ``sp_validation.__file__``'s parent so +the generated ini points at the installed package's module directory (where +``sacc_like_unions.py`` lives) regardless of checkout location. +""" + +import argparse +from pathlib import Path + + +def _sp_validation_modules(): + """The directory holding the sp_validation CosmoSIS module files. + + Resolved from the installed package so the generated ini finds + ``sacc_like_unions.py`` wherever sp_validation is installed. + """ + import sp_validation + + return str(Path(sp_validation.__file__).resolve().parent) + + +def generate_inference_config(template_path, out_path, substitutions): + """Write ``out_path`` from ``template_path`` with ``substitutions`` in DEFAULT. + + Parameters + ---------- + template_path : str or Path + The pipeline ini template (carries ``%(KEY)s`` placeholders). + out_path : str or Path + Destination ini. + substitutions : dict + ``{KEY: value}`` lines prepended into the template's ``[DEFAULT]`` + section. Every referenced ``%(KEY)s`` in the template must have a value + here (COSMOSIS_DIR already sits in the template's DEFAULT and may be + overridden). ``None`` values are dropped (an absent optional key). + """ + lines = Path(template_path).read_text().splitlines(keepends=True) + + header = "[DEFAULT]" + default_idx = next( + (i for i, line in enumerate(lines) if line.strip() == header), None + ) + if default_idx is None: + raise ValueError(f"template {template_path} has no [DEFAULT] section") + + # The end of the DEFAULT section: the next `[section]` header, or EOF. + section_end = next( + ( + i + for i in range(default_idx + 1, len(lines)) + if lines[i].lstrip().startswith("[") + ), + len(lines), + ) + + # A key the template already declares in DEFAULT is REPLACED in place (e.g. the + # template's placeholder COSMOSIS_DIR); a genuinely-new key is prepended just + # after the header. This avoids a duplicate DEFAULT key, which CosmoSIS's + # ConfigParser (strict) rejects. + wanted = {key: value for key, value in substitutions.items() if value is not None} + remaining = dict(wanted) + for i in range(default_idx + 1, section_end): + stripped = lines[i].lstrip() + if not stripped or stripped.startswith(("#", ";", "[")): + continue + existing_key = stripped.split("=", 1)[0].strip() + if existing_key in remaining: + lines[i] = f"{existing_key} = {remaining.pop(existing_key)}\n" + + prepended = [f"{key} = {value}\n" for key, value in remaining.items()] + out_lines = lines[: default_idx + 1] + prepended + lines[default_idx + 1 :] + + out_path = Path(out_path) + out_path.parent.mkdir(parents=True, exist_ok=True) + out_path.write_text("".join(out_lines)) + print( + f"Wrote {out_path} from {template_path} " + f"({len(wanted)} DEFAULT keys, {len(prepended)} new)" + ) + return str(out_path) + + +def _substitutions(scratch, cosmosis_dir, *, fits_file=None, sacc_file=None): + """Assemble the DEFAULT substitution dict, resolving SP_VALIDATION_MODULES. + + ``fits_file`` (2pt_like path) and ``sacc_file`` (sacc_like path) are mutually + the data-file placeholder for their template; whichever the template + references is filled, the other left absent. + """ + return { + "SCRATCH": scratch, + "FITS_FILE": fits_file, + "SACC_FILE": sacc_file, + "COSMOSIS_DIR": cosmosis_dir, + "SP_VALIDATION_MODULES": _sp_validation_modules(), + } + + +def _from_snakemake(smk): + p = smk.params + generate_inference_config( + template_path=smk.input[0] + if not hasattr(smk.input, "template") + else smk.input.template, + out_path=str(smk.output[0]), + substitutions=_substitutions( + scratch=p["scratch"], + cosmosis_dir=p["cosmosis_dir"], + fits_file=p.get("fits_file", None), + sacc_file=p.get("sacc_file", None), + ), + ) + + +def _from_cli(argv=None): + ap = argparse.ArgumentParser( + description="Generate a CosmoSIS pipeline ini from a template + DEFAULT subs." + ) + ap.add_argument("--template", required=True, help="Pipeline ini template") + ap.add_argument("--out", required=True, help="Output ini path") + ap.add_argument("--scratch", required=True, help="SCRATCH value") + ap.add_argument("--cosmosis-dir", required=True, help="COSMOSIS_DIR value") + ap.add_argument("--fits-file", default=None, help="FITS_FILE (2pt_like path)") + ap.add_argument("--sacc-file", default=None, help="SACC_FILE (sacc_like path)") + a = ap.parse_args(argv) + generate_inference_config( + template_path=a.template, + out_path=a.out, + substitutions=_substitutions( + scratch=a.scratch, + cosmosis_dir=a.cosmosis_dir, + fits_file=a.fits_file, + sacc_file=a.sacc_file, + ), + ) + + +if __name__ == "__main__": + try: + snakemake # noqa: F821 — injected by Snakemake's script: directive + except NameError: + _from_cli() + else: + _from_snakemake(snakemake) # noqa: F821