diff --git a/src/sp_validation/__init__.py b/src/sp_validation/__init__.py index aca7a9da..3d1ad5f1 100644 --- a/src/sp_validation/__init__.py +++ b/src/sp_validation/__init__.py @@ -12,6 +12,7 @@ "glass_mock", "masks", "plots", + "pseudo_cl", "rho_tau", "statistics", "survey", diff --git a/src/sp_validation/cosmo_val/catalog_characterization.py b/src/sp_validation/cosmo_val/catalog_characterization.py index c44caacc..bb5f503d 100644 --- a/src/sp_validation/cosmo_val/catalog_characterization.py +++ b/src/sp_validation/cosmo_val/catalog_characterization.py @@ -17,6 +17,14 @@ from astropy.io import fits from cs_util import plots as cs_plots +from ..survey import ( + additive_bias, + area_from_coords, + effective_survey_stats, + n_eff_density, +) +from ..survey import ellipticity_dispersion as ellipticity_dispersion_stat + class CatalogCharacterizationMixin: def compute_survey_stats( @@ -79,10 +87,6 @@ def compute_survey_stats( w = np.asarray(data[weight_column], dtype=float) - sum_w = float(np.sum(w)) - sum_w2 = float(np.sum(w**2)) - sum_w2_e2 = float(np.sum((w**2) * (e1**2 + e2**2))) - if mask_path is not None: if not os.path.exists(mask_path): raise FileNotFoundError(f"Mask path not found: {mask_path}") @@ -106,23 +110,22 @@ def compute_survey_stats( f"Unable to determine survey area for {ver}. Provide mask_path or nside." ) - area_arcmin2 = area_deg2 * 3600.0 - - n_eff = (sum_w**2) / (area_arcmin2 * sum_w2) if sum_w2 > 0 else 0.0 - sigma_e = np.sqrt(sum_w2_e2 / sum_w2) if sum_w2 > 0 else 0.0 + stats = effective_survey_stats(e1, e2, w, area_deg2) results = { "area_deg2": area_deg2, - "n_eff": n_eff, - "sigma_e": sigma_e, - "sum_w": sum_w, - "sum_w2": sum_w2, + "n_eff": stats["n_eff"], + "sigma_e": stats["sigma_e"], + "sum_w": stats["sum_w"], + "sum_w2": stats["sum_w2"], "catalog_size": n_rows, } if overwrite_config: self.cc[ver].setdefault("cov_th", {}).update( - A=float(area_deg2), n_e=float(n_eff), sigma_e=float(sigma_e) + A=float(area_deg2), + n_e=float(stats["n_eff"]), + sigma_e=float(stats["sigma_e"]), ) self._write_catalog_config() @@ -132,11 +135,7 @@ def _area_from_catalog(self, catalog_path, nside): data = fits.getdata(catalog_path, memmap=True) ra = np.asarray(data["RA"], dtype=float) dec = np.asarray(data["Dec"], dtype=float) - theta = np.radians(90.0 - dec) - phi = np.radians(ra) - pix = hp.ang2pix(nside, theta, phi, lonlat=False) - unique_pix = np.unique(pix) - return float(unique_pix.size * hp.nside2pixarea(nside, degrees=True)) + return area_from_coords(ra, dec, nside) def _area_from_mask(self, mask_map_path): mask = hp.read_map(mask_map_path, dtype=np.float64) @@ -210,9 +209,7 @@ def calculate_n_eff_gal(self): self.print_magenta(ver) with self.results[ver].temporarily_read_data(): w = self._read_shear_cols(ver, "w_col") - n_eff_gal[ver] = ( - 1 / (self.area[ver] * 60 * 60) * np.sum(w) ** 2 / np.sum(w**2) - ) + n_eff_gal[ver] = n_eff_density(w, self.area[ver]) print(f"n_eff_gal = {n_eff_gal[ver]:.2f} gal./arcmin^-2") self._n_eff_gal = n_eff_gal @@ -225,13 +222,7 @@ def calculate_ellipticity_dispersion(self): self.print_magenta(ver) with self.results[ver].temporarily_read_data(): e1, e2, w = self._read_shear_cols(ver, "e1_col", "e2_col", "w_col") - ellipticity_dispersion[ver] = np.sqrt( - 0.5 - * ( - np.average(e1**2, weights=w**2) - + np.average(e2**2, weights=w**2) - ) - ) + ellipticity_dispersion[ver] = ellipticity_dispersion_stat(e1, e2, w) print(f"Ellipticity dispersion = {ellipticity_dispersion[ver]:.4f}") self._ellipticity_dispersion = ellipticity_dispersion @@ -363,18 +354,9 @@ def calculate_additive_bias(self): for ver in self.versions: self.print_magenta(ver) R = self.cc[ver]["shear"]["R"] - e1_col, e2_col, w_col = [ - self.cc[ver]["shear"][k] for k in ["e1_col", "e2_col", "w_col"] - ] with self.results[ver].temporarily_read_data(): - self._c1[ver] = np.average( - self.results[ver].dat_shear[e1_col] / R, - weights=self.results[ver].dat_shear[w_col], - ) - self._c2[ver] = np.average( - self.results[ver].dat_shear[e2_col] / R, - weights=self.results[ver].dat_shear[w_col], - ) + e1, e2, w = self._read_shear_cols(ver, "e1_col", "e2_col", "w_col") + self._c1[ver], self._c2[ver] = additive_bias(e1, e2, w, R) self.print_done("Finished additive bias calculation.") @property diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 91042853..ee65f658 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -91,6 +91,9 @@ class CosmologyValidation( Apply polarization correction factor in pseudo-C_ell calculations. nrandom_cell : int, default 10 Number of random realizations for C_ell error estimation. + cell_seed : int, default 8192 + Seed for the random-rotation noise realizations in the pseudo-C_ell + noise debiasing, making those realizations reproducible run-to-run. cosmo_params : dict, optional Cosmological parameters to pass to get_cosmo(). If None, uses Planck 2018. @@ -215,6 +218,7 @@ def __init__( noise_bias_method="analytic", fiducial_input_inka="coupled", nrandom_cell=10, + cell_seed=8192, path_onecovariance=None, cosmo_params=None, blind=None, @@ -240,6 +244,7 @@ def __init__( self.ell_step = ell_step self.pol_factor = pol_factor self.nrandom_cell = nrandom_cell + self.cell_seed = cell_seed self.cell_method = cell_method self.noise_bias_method = noise_bias_method self.fiducial_input_inka = fiducial_input_inka diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index 9b51b208..e4733eea 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -17,6 +17,13 @@ from astropy.io import fits from ..cosmology import get_theo_c_ell +from ..pseudo_cl import ( + apply_random_rotation, + get_n_gal_map, + get_pseudo_cls_catalog, + get_pseudo_cls_map, + make_namaster_bin, +) from ..rho_tau import get_params_rho_tau from ..statistics import chi2_and_pte, cov_from_one_covariance @@ -36,51 +43,16 @@ def pseudo_cls_onecov(self): return self._pseudo_cls_onecov def get_namaster_bin(self, lmin, lmax, b_lmax): - """Build NaMaster binning object. - - Parameters - ---------- - lmin, lmax : int - Multipole range. - b_lmax : int - Maximum multipole for the NmtBin object. - - Returns - ------- - nmt.NmtBin - """ - ells = np.arange(lmin, lmax + 1) - - if self.binning == "linear": - bpws = (ells - lmin) // self.ell_step - bpws = np.minimum(bpws, bpws[-1]) - b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax) - elif self.binning == "logspace": - # Start geomspace at ell_min_log (>= lmin) to avoid - # sub-multipole bins at low ell that destabilize the MCM. - # All ell below ell_min_log go into bin 0 as padding. - ell_min_log = max(lmin, 50) - bins_ell = np.geomspace(ell_min_log, lmax, self.n_ell_bins + 1) - bpws = np.digitize(ells.astype(float), bins_ell) - 1 - bpws = np.clip(bpws, 0, self.n_ell_bins - 1) - b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax) - elif self.binning == "powspace": - start = np.power(lmin, self.power) - end = np.power(lmax, self.power) - bins_ell = np.power( - np.linspace(start, end, self.n_ell_bins + 1), 1 / self.power - ) - bpws = np.digitize(ells.astype(float), bins_ell) - 1 - bpws[0] = 0 - bpws[-1] = self.n_ell_bins - 1 - b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax) - else: - raise ValueError( - f"Unknown binning '{self.binning}'. " - "Choose from 'linear', 'logspace', 'powspace'." - ) - - return b + """Build NaMaster binning object (thin wrapper, state -> primitive).""" + return make_namaster_bin( + lmin, + lmax, + b_lmax, + self.binning, + ell_step=self.ell_step, + n_ell_bins=self.n_ell_bins, + power=self.power, + ) def get_variance_map(self, nside, e1, e2, w, unique_pix, idx_rep): """ @@ -197,6 +169,7 @@ def calculate_pseudo_cl_eb_cov(self): n_gal, unique_pix, idx_rep, + np.random.default_rng(self.cell_seed), ) noise_bias_cl = np.mean(cl_noise, axis=0) @@ -542,12 +515,13 @@ def calculate_pseudo_cl_map(self, ver, nside, out_path): ell_eff, cl_shear, wsp = self.get_pseudo_cls_map(shear_map, n_gal_map) cl_noise = np.zeros_like(cl_shear) + rng = np.random.default_rng(self.cell_seed) for i in range(self.nrandom_cell): noise_map_e1 = np.zeros(hp.nside2npix(nside)) noise_map_e2 = np.zeros(hp.nside2npix(nside)) - e1_rot, e2_rot = self.apply_random_rotation(e1, e2) + e1_rot, e2_rot = self.apply_random_rotation(e1, e2, rng) noise_map_e1[unique_pix] += np.bincount(idx_rep, weights=e1_rot * w) noise_map_e2[unique_pix] += np.bincount(idx_rep, weights=e2_rot * w) @@ -569,8 +543,7 @@ def calculate_pseudo_cl_map(self, ver, nside, out_path): pass del n_gal_map - # This is a problem because the measurement depends on the seed. To be fixed. - # cl_shear = cl_shear - np.mean(cl_noise, axis=1, keepdims=True) + # Noise realizations are now reproducible (seeded rng from self.cell_seed). cl_shear = cl_shear - cl_noise self.print_cyan("Saving pseudo-Cl's...") @@ -596,29 +569,19 @@ def calculate_pseudo_cl_catalog(self, ver, out_path): self._pseudo_cls[ver]["pseudo_cl"] = cl_shear def get_n_gal_map(self, params, nside, cat_gal): - """ - Compute the galaxy number density map. - """ - ra = cat_gal[params["ra_col"]] - dec = cat_gal[params["dec_col"]] - w = cat_gal[params["w_col"]] - - theta = (90.0 - dec) * np.pi / 180.0 - phi = ra * np.pi / 180.0 - pix = hp.ang2pix(nside, theta, phi) - - unique_pix, idx, idx_rep = np.unique( - pix, return_index=True, return_inverse=True + """Weighted galaxy number-density map (thin wrapper -> primitive).""" + return get_n_gal_map( + nside, + cat_gal[params["ra_col"]], + cat_gal[params["dec_col"]], + weights=cat_gal[params["w_col"]], ) - n_gal = np.zeros(hp.nside2npix(nside)) - n_gal[unique_pix] = np.bincount(idx_rep, weights=w) - return n_gal, unique_pix, idx, idx_rep def get_gaussian_real( - self, params, nside, lmax, cat_gal, n_gal, mask, unique_pix, idx_rep + self, params, nside, lmax, cat_gal, n_gal, mask, unique_pix, idx_rep, rng=None ): e1_rot, e2_rot = self.apply_random_rotation( - cat_gal[params["e1_col"]], cat_gal[params["e2_col"]] + cat_gal[params["e1_col"]], cat_gal[params["e2_col"]], rng ) noise_map_e1 = np.zeros(hp.nside2npix(nside)) noise_map_e2 = np.zeros(hp.nside2npix(nside)) @@ -632,10 +595,20 @@ def get_gaussian_real( return noise_map_e1 + 1j * noise_map_e2 def get_sample( - self, params, nside, lmax, b, cat_gal, n_gal, mask, unique_pix, idx_rep + self, + params, + nside, + lmax, + b, + cat_gal, + n_gal, + mask, + unique_pix, + idx_rep, + rng=None, ): noise_map = self.get_gaussian_real( - params, nside, lmax, cat_gal, n_gal, mask, unique_pix, idx_rep + params, nside, lmax, cat_gal, n_gal, mask, unique_pix, idx_rep, rng ) f = nmt.NmtField(mask=mask, maps=[noise_map.real, noise_map.imag], lmax=lmax) @@ -648,85 +621,39 @@ def get_sample( return cl_noise, f, wsp def get_pseudo_cls_map(self, map, mask, wsp=None): - """ - Compute the pseudo-cl for a given map. - """ - - lmin = 8 - lmax = 2 * self.nside - b_lmax = lmax - 1 - - b = self.get_namaster_bin(lmin, lmax, b_lmax) - - ell_eff = b.get_effective_ells() - - factor = -1 if self.pol_factor else 1 - - f_all = nmt.NmtField(mask=mask, maps=[map.real, factor * map.imag], lmax=b_lmax) - if wsp is None: - wsp = nmt.NmtWorkspace.from_fields(f_all, f_all, b) - - cl_coupled = nmt.compute_coupled_cell(f_all, f_all) - cl_all = wsp.decouple_cell(cl_coupled) - - return ell_eff, cl_all, wsp + """Map-based pseudo-cl (thin wrapper, state -> primitive).""" + return get_pseudo_cls_map( + map, + mask, + self.nside, + self.binning, + pol_factor=self.pol_factor, + wsp=wsp, + ell_step=self.ell_step, + n_ell_bins=self.n_ell_bins, + power=self.power, + ) def get_pseudo_cls_catalog(self, catalog, params, wsp=None): - """ - Compute the pseudo-cl for a given catalog. - """ - - lmin = 8 - lmax = 2 * self.nside - b_lmax = lmax - 1 - - b = self.get_namaster_bin(lmin, lmax, b_lmax) - - ell_eff = b.get_effective_ells() - - factor = -1 if self.pol_factor else 1 - - f_all = nmt.NmtFieldCatalog( - positions=[catalog[params["ra_col"]], catalog[params["dec_col"]]], - weights=catalog[params["w_col"]], - field=[catalog[params["e1_col"]], factor * catalog[params["e2_col"]]], - lmax=b_lmax, - lmax_mask=b_lmax, - spin=2, - lonlat=True, + """Catalog-based pseudo-cl (thin wrapper, state -> primitive).""" + return get_pseudo_cls_catalog( + catalog, + params, + self.nside, + self.binning, + pol_factor=self.pol_factor, + wsp=wsp, + ell_step=self.ell_step, + n_ell_bins=self.n_ell_bins, + power=self.power, ) - if wsp is None: - wsp = nmt.NmtWorkspace.from_fields(f_all, f_all, b) - - cl_coupled = nmt.compute_coupled_cell(f_all, f_all) - cl_all = wsp.decouple_cell(cl_coupled) - - return ell_eff, cl_all, wsp + def apply_random_rotation(self, e1, e2, rng=None): + """Random ellipticity rotation (thin wrapper -> primitive). - def apply_random_rotation(self, e1, e2): - """ - Apply a random rotation to the ellipticity components e1 and e2. - - Parameters - ---------- - e1 : np.array - First component of the ellipticity. - e2 : np.array - Second component of the ellipticity. - - Returns - ------- - np.array - First component of the rotated ellipticity. - np.array - Second component of the rotated ellipticity. + Pass a seeded ``rng`` for reproducible noise realizations. """ - np.random.seed() - rot_angle = np.random.rand(len(e1)) * 2 * np.pi - e1_out = e1 * np.cos(rot_angle) + e2 * np.sin(rot_angle) - e2_out = -e1 * np.sin(rot_angle) + e2 * np.cos(rot_angle) - return e1_out, e2_out + return apply_random_rotation(e1, e2, rng) def save_pseudo_cl(self, ell_eff, pseudo_cl, out_path): """ diff --git a/src/sp_validation/glass_mock.py b/src/sp_validation/glass_mock.py index e7a184b8..7ba0057b 100644 --- a/src/sp_validation/glass_mock.py +++ b/src/sp_validation/glass_mock.py @@ -34,7 +34,6 @@ "growth_factor", "ia_convergence", "create_mask_from_catalogue", - "powspace_bins", "compute_two_point_xi", "compute_two_point_cl", "compute_two_point_cl_map", @@ -336,28 +335,21 @@ def create_mask_from_catalogue(nside, path, output, ra_col="RA", dec_col="DEC"): } -def powspace_bins(nside=1024, lmin=8, n_bins=32): - """NaMaster powspace bandpower binning used across the mock harmonic stats. +def _mock_powspace_bin(nside, lmin, n_bins): + """Powspace NaMaster binning for the mock harmonic stats, via the primitive. - Returns ``(bin, ell_eff, lmax, b_lmax)``: the NaMaster ``NmtBin`` built on - square-root-spaced bandpowers between ``lmin`` and ``lmax = 2 * nside``, - plus the effective ells and the two lmax values the callers reuse. + Returns ``(bin, ell_eff, lmax, b_lmax)`` — the 4-tuple the mock harmonic + helpers consume. The binning itself is the shared + ``pseudo_cl.make_namaster_bin`` powspace scheme (``power=0.5`` reproduces the + old square-root spacing exactly); ``lmax = 2 * nside`` / ``b_lmax = lmax - 1`` + come from ``pseudo_cl.pseudo_cl_geometry``. Note the mock's ``lmin`` floor is + a call argument here, whereas the geometry helper fixes ``lmin = 8``; the two + coincide for the production default, and this wrapper keeps the override. """ - import pymaster as nmt - - lmax = 2 * nside - b_lmax = lmax - 1 - - ells = np.arange(lmin, lmax + 1) - start = np.power(lmin, 1 / 2) - end = np.power(lmax, 1 / 2) - bins_ell = np.power(np.linspace(start, end, n_bins + 1), 2) + from sp_validation.pseudo_cl import make_namaster_bin, pseudo_cl_geometry - bpws = np.digitize(ells.astype(float), bins_ell) - 1 - bpws[0] = 0 - bpws[-1] = n_bins - 1 - - b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax) + _lmin, lmax, b_lmax = pseudo_cl_geometry(nside) + b = make_namaster_bin(lmin, lmax, b_lmax, "powspace", n_ell_bins=n_bins, power=0.5) return b, b.get_effective_ells(), lmax, b_lmax @@ -380,23 +372,29 @@ def compute_two_point_xi(cat, config=None): def get_n_gal_map(nside, ra, dec): - """Galaxy-count HEALPix map plus the unique-pixel bookkeeping arrays.""" - import healpy as hp + """Galaxy-count HEALPix map plus the unique-pixel bookkeeping arrays. - theta = (90.0 - dec) * np.pi / 180.0 - phi = ra * np.pi / 180.0 - pix = hp.ang2pix(nside, theta, phi) + The unweighted twin of the pseudo-Cl ``get_n_gal_map`` primitive: delegates + to it with ``weights=None`` (counts). Imported lazily so this module keeps + resolving in CAMB-only environments without the harmonic stack. + """ + from sp_validation.pseudo_cl import get_n_gal_map as _get_n_gal_map - unique_pix, idx, idx_rep = np.unique(pix, return_index=True, return_inverse=True) - n_gal = np.zeros(hp.nside2npix(nside)) - n_gal[unique_pix] = np.bincount(idx_rep) - return n_gal, unique_pix, idx, idx_rep + return _get_n_gal_map(nside, ra, dec) def compute_two_point_cl(cat, nside=1024, lmin=8, n_bins=32): """Catalogue-based NaMaster pseudo-Cl for a mock shear field. Returns ``(cl_coupled, cl_all)``, each with the ell axis prepended. + + Deliberately keeps its own ``NmtFieldCatalog`` construction rather than + delegating to the ``pseudo_cl.get_pseudo_cls_catalog`` primitive: the mock + field is UNWEIGHTED (``weights = ones``, no per-object shear weight) and this + function returns the COUPLED pseudo-Cl alongside the decoupled one (both with + an ell axis prepended), which the primitive's ``(ell_eff, cl_all, wsp)`` + contract does not expose. The shared piece — the powspace binning — is taken + from ``_mock_powspace_bin`` so the bandpower drift is gone. """ import pymaster as nmt @@ -408,7 +406,7 @@ def compute_two_point_cl(cat, nside=1024, lmin=8, n_bins=32): ra[ra < 0] += 360 - b, ell_eff, lmax, b_lmax = powspace_bins(nside=nside, lmin=lmin, n_bins=n_bins) + b, ell_eff, lmax, b_lmax = _mock_powspace_bin(nside, lmin, n_bins) factor = -1 f_all = nmt.NmtFieldCatalog( @@ -436,6 +434,15 @@ def compute_two_point_cl_map(cat, nside=1024, lmin=8, n_bins=32): Bins the galaxies onto a HEALPix shear map (count-weighted mean per pixel) and runs the MCM on the resulting masked field. Returns ``(cl_coupled, cl_all)`` with the ell axis prepended. + + Deliberately keeps its own ``NmtField`` construction rather than delegating + to the ``pseudo_cl.get_pseudo_cls_map`` primitive: the mock shear map is + built from UNWEIGHTED galaxy counts (the mask is the count map and the + per-pixel mean divides by counts, not summed shear weights), and this + function returns the COUPLED pseudo-Cl alongside the decoupled one (both with + an ell axis prepended), which the primitive's ``(ell_eff, cl_all, wsp)`` + contract does not expose. The shared pieces — the powspace binning and the + galaxy-count map — come from ``_mock_powspace_bin`` and ``get_n_gal_map``. """ import healpy as hp import pymaster as nmt @@ -448,7 +455,7 @@ def compute_two_point_cl_map(cat, nside=1024, lmin=8, n_bins=32): ra[ra < 0] += 360 - b, ell_eff, lmax, b_lmax = powspace_bins(nside=nside, lmin=lmin, n_bins=n_bins) + b, ell_eff, lmax, b_lmax = _mock_powspace_bin(nside, lmin, n_bins) factor = -1 n_gal_map, unique_pix, _idx, idx_rep = get_n_gal_map(nside, ra, dec) @@ -498,7 +505,7 @@ def compute_leakage_harmony(cat, cat_star, nside=1024, lmin=8, n_bins=32): ra[ra < 0] += 360 ra_star[ra_star < 0] += 360 - b, ell_eff, _lmax, b_lmax = powspace_bins(nside=nside, lmin=lmin, n_bins=n_bins) + b, ell_eff, _lmax, b_lmax = _mock_powspace_bin(nside, lmin, n_bins) factor = -1 f_psf = nmt.NmtFieldCatalog( diff --git a/src/sp_validation/pseudo_cl.py b/src/sp_validation/pseudo_cl.py new file mode 100644 index 00000000..c9355ec9 --- /dev/null +++ b/src/sp_validation/pseudo_cl.py @@ -0,0 +1,282 @@ +"""Pseudo-Cl / harmonic-space estimator primitives for cosmology validation. + +Stateless shared primitives for the pseudo-Cl (harmonic-space) estimator, +mirroring ``b_modes.py`` and ``rho_tau.py``: the orchestrator mixin in +``sp_validation.cosmo_val.pseudo_cl`` (and analysis scripts directly) call these +free functions. Everything here is pure computation -- NaMaster binning, +weighted galaxy number-density maps, random-rotation noise debiasing, and the +map-/catalog-based pseudo-Cl estimators. Depends on pymaster (NaMaster) and +healpy. + +The harmonic geometry the estimators use is fixed by ``nside``: ``lmin = 8``, +``lmax = 2 * nside``, ``b_lmax = lmax - 1``. ``pseudo_cl_geometry`` returns that +triple so the binning and the fields stay in lockstep. +""" + +import healpy as hp +import numpy as np +import pymaster as nmt + +# Lowest multipole retained by the pseudo-Cl estimators. +LMIN = 8 + + +def pseudo_cl_geometry(nside): + """Return ``(lmin, lmax, b_lmax)`` for the pseudo-Cl estimator at ``nside``. + + ``lmax = 2 * nside`` is the NaMaster band-power ceiling; ``b_lmax = lmax - 1`` + is the field/binning ``lmax``. ``lmin`` is the fixed low-ell floor. + """ + lmax = 2 * nside + return LMIN, lmax, lmax - 1 + + +def make_namaster_bin( + lmin, lmax, b_lmax, binning, *, ell_step=10, n_ell_bins=32, power=0.5 +): + """Build a NaMaster binning object for one of the supported schemes. + + Parameters + ---------- + lmin, lmax : int + Multipole range. + b_lmax : int + Maximum multipole for the ``NmtBin`` object. + binning : {'linear', 'logspace', 'powspace'} + Binning scheme. + ell_step : int, optional + Bin width in ell for ``'linear'`` binning. + n_ell_bins : int, optional + Number of ell bins for ``'logspace'`` / ``'powspace'`` binning. + power : float, optional + Exponent for ``'powspace'`` binning. + + Returns + ------- + nmt.NmtBin + """ + ells = np.arange(lmin, lmax + 1) + + if binning == "linear": + bpws = (ells - lmin) // ell_step + bpws = np.minimum(bpws, bpws[-1]) + b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax) + elif binning == "logspace": + # Start geomspace at ell_min_log (>= lmin) to avoid + # sub-multipole bins at low ell that destabilize the MCM. + # All ell below ell_min_log go into bin 0 as padding. + ell_min_log = max(lmin, 50) + bins_ell = np.geomspace(ell_min_log, lmax, n_ell_bins + 1) + bpws = np.digitize(ells.astype(float), bins_ell) - 1 + bpws = np.clip(bpws, 0, n_ell_bins - 1) + b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax) + elif binning == "powspace": + start = np.power(lmin, power) + end = np.power(lmax, power) + bins_ell = np.power(np.linspace(start, end, n_ell_bins + 1), 1 / power) + bpws = np.digitize(ells.astype(float), bins_ell) - 1 + bpws[0] = 0 + bpws[-1] = n_ell_bins - 1 + b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax) + else: + raise ValueError( + f"Unknown binning '{binning}'. " + "Choose from 'linear', 'logspace', 'powspace'." + ) + + return b + + +def get_n_gal_map(nside, ra, dec, weights=None): + """Weighted galaxy number-density HEALPix map plus pixel bookkeeping. + + Bins ``(ra, dec)`` (degrees) onto an ``nside`` HEALPix grid. With + ``weights=None`` each object counts as 1 (galaxy counts); pass per-object + ``weights`` for a weight-summed occupancy map. + + Returns + ------- + n_gal : np.ndarray + Map of summed weights (or counts) per pixel, shape ``(npix,)``. + unique_pix : np.ndarray + Sorted unique occupied pixel indices. + idx : np.ndarray + First-occurrence indices into the input from ``np.unique``. + idx_rep : np.ndarray + Inverse map: pixel-group index for each input object. + """ + theta = (90.0 - dec) * np.pi / 180.0 + phi = ra * np.pi / 180.0 + pix = hp.ang2pix(nside, theta, phi) + + unique_pix, idx, idx_rep = np.unique(pix, return_index=True, return_inverse=True) + n_gal = np.zeros(hp.nside2npix(nside)) + n_gal[unique_pix] = np.bincount(idx_rep, weights=weights) + return n_gal, unique_pix, idx, idx_rep + + +def apply_random_rotation(e1, e2, rng=None): + """Apply a uniform random rotation to ellipticity components. + + Parameters + ---------- + e1, e2 : np.ndarray + Ellipticity components. + rng : np.random.Generator, optional + Random generator for the rotation angles. Pass a seeded generator + (e.g. ``np.random.default_rng(seed)``) for reproducible draws; when + ``None`` a fresh entropy-seeded generator is used (non-reproducible). + + Returns + ------- + e1_out, e2_out : np.ndarray + Rotated ellipticity components. + """ + if rng is None: + rng = np.random.default_rng() + rot_angle = rng.random(len(e1)) * 2 * np.pi + e1_out = e1 * np.cos(rot_angle) + e2 * np.sin(rot_angle) + e2_out = -e1 * np.sin(rot_angle) + e2 * np.cos(rot_angle) + return e1_out, e2_out + + +def get_pseudo_cls_map( + shear_map, + mask, + nside, + binning, + *, + pol_factor=True, + wsp=None, + ell_step=10, + n_ell_bins=32, + power=0.5, +): + """Map-based pseudo-Cl for a complex shear map. + + Parameters + ---------- + shear_map : np.ndarray + Complex shear map (``e1 + 1j * e2``). + mask : np.ndarray + Field mask (the galaxy number-density map). + nside : int + HEALPix resolution; fixes the harmonic geometry. + binning : str + Binning scheme passed to :func:`make_namaster_bin`. + pol_factor : bool, optional + If ``True`` flip the sign of the imaginary (e2) component. + wsp : nmt.NmtWorkspace, optional + Reuse a coupling workspace; built from the field if ``None``. + ell_step, n_ell_bins, power : optional + Binning-scheme parameters forwarded to :func:`make_namaster_bin`. + + Returns + ------- + ell_eff : np.ndarray + Effective multipoles of the bandpowers. + cl_all : np.ndarray + Decoupled EE/EB/BE/BB spectra, shape ``(4, n_bands)``. + wsp : nmt.NmtWorkspace + The coupling workspace (newly built or the one passed in). + """ + lmin, lmax, b_lmax = pseudo_cl_geometry(nside) + + b = make_namaster_bin( + lmin, + lmax, + b_lmax, + binning, + ell_step=ell_step, + n_ell_bins=n_ell_bins, + power=power, + ) + ell_eff = b.get_effective_ells() + + factor = -1 if pol_factor else 1 + + f_all = nmt.NmtField( + mask=mask, maps=[shear_map.real, factor * shear_map.imag], lmax=b_lmax + ) + if wsp is None: + wsp = nmt.NmtWorkspace.from_fields(f_all, f_all, b) + + cl_coupled = nmt.compute_coupled_cell(f_all, f_all) + cl_all = wsp.decouple_cell(cl_coupled) + + return ell_eff, cl_all, wsp + + +def get_pseudo_cls_catalog( + catalog, + params, + nside, + binning, + *, + pol_factor=True, + wsp=None, + ell_step=10, + n_ell_bins=32, + power=0.5, +): + """Catalog-based pseudo-Cl via NaMaster's ``NmtFieldCatalog``. + + Parameters + ---------- + catalog : np.ndarray + Structured catalog array with the columns named in ``params``. + params : dict + Column-name mapping (``ra_col``, ``dec_col``, ``w_col``, ``e1_col``, + ``e2_col``). + nside : int + HEALPix resolution; fixes the harmonic geometry. + binning : str + Binning scheme passed to :func:`make_namaster_bin`. + pol_factor : bool, optional + If ``True`` flip the sign of the e2 component. + wsp : nmt.NmtWorkspace, optional + Reuse a coupling workspace; built from the field if ``None``. + ell_step, n_ell_bins, power : optional + Binning-scheme parameters forwarded to :func:`make_namaster_bin`. + + Returns + ------- + ell_eff : np.ndarray + Effective multipoles of the bandpowers. + cl_all : np.ndarray + Decoupled EE/EB/BE/BB spectra, shape ``(4, n_bands)``. + wsp : nmt.NmtWorkspace + The coupling workspace (newly built or the one passed in). + """ + lmin, lmax, b_lmax = pseudo_cl_geometry(nside) + + b = make_namaster_bin( + lmin, + lmax, + b_lmax, + binning, + ell_step=ell_step, + n_ell_bins=n_ell_bins, + power=power, + ) + ell_eff = b.get_effective_ells() + + factor = -1 if pol_factor else 1 + + f_all = nmt.NmtFieldCatalog( + positions=[catalog[params["ra_col"]], catalog[params["dec_col"]]], + weights=catalog[params["w_col"]], + field=[catalog[params["e1_col"]], factor * catalog[params["e2_col"]]], + lmax=b_lmax, + lmax_mask=b_lmax, + spin=2, + lonlat=True, + ) + + if wsp is None: + wsp = nmt.NmtWorkspace.from_fields(f_all, f_all, b) + + cl_coupled = nmt.compute_coupled_cell(f_all, f_all) + cl_all = wsp.decouple_cell(cl_coupled) + + return ell_eff, cl_all, wsp diff --git a/src/sp_validation/survey.py b/src/sp_validation/survey.py index 9b2d10f2..fe35cde9 100644 --- a/src/sp_validation/survey.py +++ b/src/sp_validation/survey.py @@ -10,6 +10,9 @@ import os from collections import Counter +import healpy as hp +import numpy as np + def get_area(dd, area_tile, verbose=False): """Get area. @@ -249,3 +252,155 @@ def get_footprint(patch, ra, dec): else: return dec > dec_min + + +def area_from_coords(ra, dec, nside): + """Survey area from galaxy coordinates via HEALPix pixel counting. + + Bins the galaxy positions onto a HEALPix grid at resolution ``nside`` and + sums the area of every occupied pixel. This ignores partial-pixel coverage + and so overestimates the true observed area, but needs no external mask. + + Parameters + ---------- + ra : array of float + Right-ascension coordinates in degrees. + dec : array of float + Declination coordinates in degrees. + nside : int + HEALPix resolution parameter. + + Returns + ------- + float + Area in square degrees. + """ + theta = np.radians(90.0 - np.asarray(dec, dtype=float)) + phi = np.radians(np.asarray(ra, dtype=float)) + pix = hp.ang2pix(nside, theta, phi, lonlat=False) + unique_pix = np.unique(pix) + return float(unique_pix.size * hp.nside2pixarea(nside, degrees=True)) + + +def n_eff_density(w, area_deg2): + """Effective galaxy number density per square arcminute. + + Uses the weighted definition ``n_eff = (Σw)² / (A · Σw²)`` where ``A`` is the + survey area in square arcminutes (``area_deg2 · 3600``). Returns ``0.0`` when + ``Σw² == 0`` to avoid division by zero. + + Parameters + ---------- + w : array of float + Per-galaxy weights. + area_deg2 : float + Survey area in square degrees. + + Returns + ------- + float + Effective number density in galaxies per square arcminute. + """ + w = np.asarray(w, dtype=float) + sum_w = float(np.sum(w)) + sum_w2 = float(np.sum(w**2)) + area_arcmin2 = area_deg2 * 3600.0 + return (sum_w**2) / (area_arcmin2 * sum_w2) if sum_w2 > 0 else 0.0 + + +def ellipticity_dispersion(e1, e2, w): + """Per-component weighted ellipticity dispersion. + + Computes ``sqrt(0.5 · (⟨e1²⟩ + ⟨e2²⟩))`` where each component average is + weighted by ``w²``. This is the *per-component* RMS convention; it differs by + a factor ``√2`` from the summed-component ``sigma_e`` returned by + :func:`effective_survey_stats`. Keep the two conventions distinct. + + Parameters + ---------- + e1, e2 : array of float + Ellipticity components. + w : array of float + Per-galaxy weights. + + Returns + ------- + float + Per-component ellipticity dispersion. + """ + e1 = np.asarray(e1, dtype=float) + e2 = np.asarray(e2, dtype=float) + w = np.asarray(w, dtype=float) + return float( + np.sqrt( + 0.5 * (np.average(e1**2, weights=w**2) + np.average(e2**2, weights=w**2)) + ) + ) + + +def additive_bias(e1, e2, w, R): + """Weighted additive-bias estimates ``(c1, c2)``. + + Returns the weighted mean of the response-corrected ellipticities + ``⟨e1 / R⟩`` and ``⟨e2 / R⟩``, weighted by ``w``. + + Parameters + ---------- + e1, e2 : array of float + Ellipticity components (before response correction). + w : array of float + Per-galaxy weights. + R : float + Multiplicative shear response. + + Returns + ------- + tuple of float + ``(c1, c2)`` additive-bias estimates. + """ + e1 = np.asarray(e1, dtype=float) + e2 = np.asarray(e2, dtype=float) + w = np.asarray(w, dtype=float) + c1 = float(np.average(e1 / R, weights=w)) + c2 = float(np.average(e2 / R, weights=w)) + return c1, c2 + + +def effective_survey_stats(e1, e2, w, area_deg2): + """Effective survey statistics for a shear catalogue. + + Bundles the area-dependent number density and the summed-component shape + noise. Here ``sigma_e = sqrt(Σ[w²(e1² + e2²)] / Σw²)`` sums over *both* + components (no ``0.5`` factor); this is distinct from the per-component + :func:`ellipticity_dispersion`, differing by a factor ``√2``. + + Parameters + ---------- + e1, e2 : array of float + Ellipticity components. + w : array of float + Per-galaxy weights. + area_deg2 : float + Survey area in square degrees. + + Returns + ------- + dict + Keys ``n_eff``, ``sigma_e``, ``sum_w``, ``sum_w2``. + """ + e1 = np.asarray(e1, dtype=float) + e2 = np.asarray(e2, dtype=float) + w = np.asarray(w, dtype=float) + + sum_w = float(np.sum(w)) + sum_w2 = float(np.sum(w**2)) + sum_w2_e2 = float(np.sum((w**2) * (e1**2 + e2**2))) + + sigma_e = np.sqrt(sum_w2_e2 / sum_w2) if sum_w2 > 0 else 0.0 + + return { + "n_eff": n_eff_density(w, area_deg2), + "sigma_e": sigma_e, + "sum_w": sum_w, + "sum_w2": sum_w2, + } diff --git a/src/sp_validation/tests/test_pseudo_cl.py b/src/sp_validation/tests/test_pseudo_cl.py new file mode 100644 index 00000000..45a2366d --- /dev/null +++ b/src/sp_validation/tests/test_pseudo_cl.py @@ -0,0 +1,592 @@ +"""GOLDEN-VALUE CHARACTERIZATION TESTS FOR THE PSEUDO-Cl ESTIMATOR. + +This module pins the numeric behavior of the pseudo-Cl (harmonic-space) +primitives in ``sp_validation.cosmo_val.pseudo_cl.PseudoClMixin`` against a +fixed, deterministic, in-memory synthetic catalog (seeded ``np.random``, a +small ``nside=64`` HEALPix grid -- no cluster data, no catalogue files). Every +pinned literal was produced by an actual run of the estimator inside the +sp_validation container. + +WHY THIS EXISTS +--------------- +The pseudo-Cl estimator is a Paper-II deliverable that, before this file, had +ZERO numeric pins. Track B2 will extract these mixin methods into free +functions in a new ``sp_validation/pseudo_cl.py`` module. These golden values +are the safety net: an extraction that silently changes the numbers fails here +rather than staying green on a smoke test. + +The pinned primitives (the exact functions Track B2 moves): + +* ``get_namaster_bin`` -- NaMaster binning, all three schemes +* ``get_n_gal_map`` -- weighted galaxy number-density map +* ``get_pseudo_cls_map`` -- map-based pseudo-Cl (decoupled EE/EB/BE/BB) +* ``get_pseudo_cls_catalog`` -- catalog-based pseudo-Cl (NmtFieldCatalog) +* ``calculate_pseudo_cl_catalog`` -- the deterministic end-to-end catalog path +* ``apply_random_rotation`` -- pinned by INVARIANT + seed reproducibility (see below) + +REPRODUCIBILITY / TOLERANCE (measured over 3 independent container processes) +---------------------------------------------------------------------------- +* ``get_namaster_bin`` / ``get_n_gal_map`` / ``get_pseudo_cls_map`` are + BITWISE-identical across processes (pure binning math; map-based NaMaster + reduction is deterministic). Pinned at ``rtol=1e-9`` -- far inside the + float-noise floor of 0, yet a hair of margin so an exact-equality assertion + doesn't make the suite brittle to a future NaMaster build. +* ``get_pseudo_cls_catalog`` (the ``NmtFieldCatalog`` path) and the end-to-end + ``calculate_pseudo_cl_catalog`` are NOT bitwise-reproducible: their reduction + order drifts at the ~2e-12 relative level (worst case measured 2.1e-12, + absolute 2.7e-19). Pinned at ``rtol=1e-6 / atol=1e-12`` -- ~6 orders of + magnitude above that float-noise floor (no flakiness margin consumed) yet + tight enough that a sub-percent change in any band bites. This matches the + precedent in ``test_b_modes.py`` and the pure-E/B pins in ``test_cosmo_val``. +* ``apply_random_rotation`` takes an optional ``rng``: with a seeded generator + the draw is reproducible, with ``rng=None`` it is non-deterministic. It is + pinned by the rotation INVARIANT (magnitude conservation + ``e1_rot^2 + e2_rot^2 == e1^2 + e2^2`` and shape, both exact) plus a + reproducibility test (same seed -> identical, different seed -> differs). + The noise debiasing now threads a seeded ``rng`` (``CosmologyValidation``'s + ``cell_seed``) so the pseudo-Cl realizations are reproducible run-to-run. + +:Author: claude-opus on behalf of Cail + +""" + +import os + +import numpy as np +import numpy.testing as npt +import pytest +import yaml + +from sp_validation.cosmo_val import CosmologyValidation +from sp_validation.rho_tau import get_params_rho_tau + +# These tests need the full harmonic-space stack (pymaster/NaMaster + healpy), +# i.e. they run inside the sp_validation container. Skip cleanly otherwise so +# the rest of the suite still collects. +pymaster = pytest.importorskip("pymaster") +healpy = pytest.importorskip("healpy") + +from astropy.io import fits # noqa: E402 (after importorskip) + +# --- Fixed experiment configuration ----------------------------------------- +NSIDE = 64 +SEED = 1234 +N_ELL_BINS = 8 + + +# --------------------------------------------------------------------------- +# Synthetic-catalog fixture (deterministic; no cluster data) +# --------------------------------------------------------------------------- +def _write_synthetic_config(tmp_path): + """Write a tiny seeded shear catalog + dndz, return (config_path, version). + + A 5000-object random-shear catalog over a 20x20 deg patch, plus the minimal + config blocks the pseudo-Cl seam reads (``shear`` for column names, ``psf`` + so ``get_params_rho_tau`` is satisfied). Fully determined by ``SEED``. + """ + from astropy.table import Table + + rng = np.random.default_rng(SEED) + version = "TestCatalog" + + cat_dir = tmp_path / "catalog" + nz_dir = tmp_path / "nz" + output_dir = tmp_path / "output" + for d in (cat_dir, nz_dir, output_dir): + d.mkdir() + + n_gal = 5000 + ra = rng.uniform(10.0, 30.0, n_gal) + dec = rng.uniform(10.0, 30.0, n_gal) + e1 = rng.normal(0, 0.25, n_gal) + e2 = rng.normal(0, 0.25, n_gal) + w = rng.uniform(0.5, 1.0, n_gal) + Table({"RA": ra, "Dec": dec, "e1": e1, "e2": e2, "w": w}).write( + cat_dir / "shear.fits", overwrite=True + ) + + z_edges = np.linspace(0.05, 3.0, 31) + dndz = np.exp(-(((z_edges - 0.7) / 0.3) ** 2)) + lines = ["# z dn_dz"] + [f"{zz} {nn}" for zz, nn in zip(z_edges, dndz)] + (nz_dir / "dndz_SP_A.txt").write_text("\n".join(lines) + "\n") + + shear_cfg = { + "path": "shear.fits", + "w_col": "w", + "e1_col": "e1", + "e2_col": "e2", + "R": 1.0, + "e1_col_corrected": "e1", + "e2_col_corrected": "e2", + "ra_col": "RA", + "dec_col": "Dec", + } + # Minimal psf block so get_params_rho_tau() succeeds; the pseudo-Cl + # primitives only read the shear-side keys, but the production call path + # routes through get_params_rho_tau, so we mirror it. + psf_cfg = { + "path": "shear.fits", + "ra_col": "RA", + "dec_col": "Dec", + "e1_PSF_col": "e1", + "e2_PSF_col": "e2", + "e1_star_col": "e1", + "e2_star_col": "e2", + "PSF_size": "w", + "star_size": "w", + "PSF_flag": "w", + "star_flag": "w", + } + config_data = { + "nz": {"subdir": str(nz_dir), "dndz": {"blind": "A", "path": "dndz"}}, + "paths": {"output": str(output_dir)}, + version: { + "subdir": str(cat_dir), + "pipeline": "SP", + "shear": shear_cfg, + "star": {**psf_cfg}, + "psf": psf_cfg, + }, + } + config_path = tmp_path / "config.yaml" + config_path.write_text(yaml.dump(config_data, sort_keys=False)) + return str(config_path), version + + +@pytest.fixture +def cv(tmp_path): + """CosmologyValidation on the synthetic catalog, pinned pseudo-Cl config.""" + config_path, version = _write_synthetic_config(tmp_path) + cv = CosmologyValidation( + versions=[version], + catalog_config=config_path, + output_dir=str(tmp_path / "output"), + nside=NSIDE, + binning="powspace", + power=0.5, + n_ell_bins=N_ELL_BINS, + pol_factor=True, + ) + cv._test_version = version + return cv + + +@pytest.fixture +def cat_and_params(cv): + """(catalog ndarray, params dict) for the synthetic version.""" + ver = cv._test_version + params = get_params_rho_tau(cv.cc[ver], survey=ver) + cat_gal = fits.getdata(cv.cc[ver]["shear"]["path"]) + return cat_gal, params + + +# Tolerances ---------------------------------------------------------------- +# Bitwise-stable primitives (binning math, n_gal map, map-based pseudo-Cl). +RTOL_DET = 1e-9 +ATOL_DET = 1e-30 +# Catalog (NmtFieldCatalog) path: drifts ~2e-12 cross-process; pin well above. +RTOL_CAT = 1e-6 +ATOL_CAT = 1e-12 + +LMIN = 8 +LMAX = 2 * NSIDE +B_LMAX = LMAX - 1 + + +# =========================================================================== +# get_namaster_bin -- all three binning schemes +# =========================================================================== +class TestGetNamasterBin: + """Pin the NaMaster binning object for each supported scheme.""" + + def test_powspace(self, cv): + cv.binning, cv.power, cv.n_ell_bins = "powspace", 0.5, 8 + b = cv.get_namaster_bin(LMIN, LMAX, B_LMAX) + assert b.get_n_bands() == 8 + npt.assert_allclose( + b.get_effective_ells(), + np.array([11.5, 20.0, 30.5, 43.5, 58.5, 75.5, 95.0, 116.5]), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_array_equal( + np.asarray(b.get_ell_list(0)), + np.array([8, 9, 10, 11, 12, 13, 14, 15]), + ) + + def test_logspace(self, cv): + cv.binning, cv.n_ell_bins = "logspace", 8 + b = cv.get_namaster_bin(LMIN, LMAX, B_LMAX) + assert b.get_n_bands() == 8 + npt.assert_allclose( + b.get_effective_ells(), + np.array([32.0, 60.0, 67.5, 75.5, 84.5, 95.5, 107.5, 120.5]), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + # Low-ell padding: all multipoles below ell_min_log(=50) fall in bin 0. + npt.assert_array_equal( + np.asarray(b.get_ell_list(0)), np.arange(8, 57, dtype=np.int32) + ) + + def test_linear(self, cv): + cv.binning, cv.ell_step = "linear", 10 + b = cv.get_namaster_bin(LMIN, LMAX, B_LMAX) + assert b.get_n_bands() == 12 + npt.assert_allclose( + b.get_effective_ells(), + np.array( + [ + 12.5, + 22.5, + 32.5, + 42.5, + 52.5, + 62.5, + 72.5, + 82.5, + 92.5, + 102.5, + 112.5, + 122.5, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_array_equal( + np.asarray(b.get_ell_list(0)), + np.array([8, 9, 10, 11, 12, 13, 14, 15, 16, 17]), + ) + + def test_unknown_binning_raises(self, cv): + cv.binning = "nonsense" + with pytest.raises(ValueError, match="Unknown binning"): + cv.get_namaster_bin(LMIN, LMAX, B_LMAX) + + +# =========================================================================== +# get_n_gal_map -- weighted galaxy number-density map +# =========================================================================== +def test_get_n_gal_map(cv, cat_and_params): + cat_gal, params = cat_and_params + n_gal, unique_pix, idx, idx_rep = cv.get_n_gal_map(params, NSIDE, cat_gal) + + # Structural invariants of the HEALPix occupancy map. + assert n_gal.shape == (healpy.nside2npix(NSIDE),) + assert int(np.count_nonzero(n_gal)) == 485 + assert unique_pix.size == 485 + assert idx_rep.size == 5000 # one entry per galaxy + # The map is supported exactly on the occupied pixels. + npt.assert_array_equal(np.nonzero(n_gal)[0], np.sort(unique_pix)) + + # Pinned scalar summaries: total weight is conserved (sum of weights), + # peak occupancy, and the pixel index bookkeeping. + npt.assert_allclose(n_gal.sum(), 3760.657282591494, rtol=RTOL_DET) + npt.assert_allclose(n_gal.max(), 16.063410433711447, rtol=RTOL_DET) + npt.assert_array_equal(unique_pix[:5], [12167, 12168, 12169, 12170, 12171]) + assert int(unique_pix.sum()) == 7849989 + npt.assert_allclose( + n_gal[unique_pix][:5], + np.array( + [ + 3.733980174372703, + 3.4379134578473773, + 4.573348360592046, + 7.926927899839571, + 4.5074633637396575, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + + +# =========================================================================== +# get_pseudo_cls_map -- map-based pseudo-Cl (bitwise-stable) +# =========================================================================== +def _build_shear_map(cv, cat_gal, params): + """Replicate calculate_pseudo_cl_map's weighted shear-map construction.""" + n_gal, unique_pix, _idx, idx_rep = cv.get_n_gal_map(params, NSIDE, cat_gal) + w = cat_gal[params["w_col"]] + e1 = cat_gal[params["e1_col"]] + e2 = cat_gal[params["e2_col"]] + mask = n_gal != 0 + m1 = np.zeros(n_gal.size) + m2 = np.zeros(n_gal.size) + m1[unique_pix] += np.bincount(idx_rep, weights=e1 * w) + m2[unique_pix] += np.bincount(idx_rep, weights=e2 * w) + m1[mask] /= n_gal[mask] + m2[mask] /= n_gal[mask] + return m1 + 1j * m2, n_gal + + +def test_get_pseudo_cls_map(cv, cat_and_params): + cat_gal, params = cat_and_params + shear_map, n_gal = _build_shear_map(cv, cat_gal, params) + ell_eff, cl_all, wsp = cv.get_pseudo_cls_map(shear_map, n_gal) + + assert cl_all.shape == (4, N_ELL_BINS) + npt.assert_allclose( + ell_eff, + np.array([11.5, 20.0, 30.5, 43.5, 58.5, 75.5, 95.0, 116.5]), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_allclose( + cl_all[0], # EE + np.array( + [ + -3.5657848980278575e-06, + 2.834003490842137e-06, + 5.141388691413763e-07, + 1.4864803917859117e-06, + 1.529668502385154e-06, + 9.126381881185905e-07, + 8.60858982887623e-07, + 2.324190517614256e-06, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_allclose( + cl_all[1], # EB + np.array( + [ + 1.1955866959532166e-05, + -2.422720314939584e-06, + -9.9636714751786e-07, + -6.033940105986599e-07, + 7.694521848269152e-07, + 2.2127016030540707e-07, + -1.716732270067554e-07, + 1.290575887877422e-07, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_allclose( + cl_all[3], # BB + np.array( + [ + 1.0742046525039353e-05, + -8.20401788619865e-07, + 1.1220545996996594e-06, + 9.428235070512487e-07, + 1.8882197714480114e-06, + 1.5461915098836886e-06, + 1.8834735730847215e-06, + 1.8974610324105966e-06, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + # BE (index 2) is the transpose-symmetric partner of EB for an auto-spectrum. + npt.assert_allclose(cl_all[2], cl_all[1], rtol=RTOL_DET, atol=1e-18) + + +# =========================================================================== +# get_pseudo_cls_catalog -- NmtFieldCatalog path (drifts ~2e-12) +# =========================================================================== +def test_get_pseudo_cls_catalog(cv, cat_and_params): + cat_gal, params = cat_and_params + ell_eff, cl_all, wsp = cv.get_pseudo_cls_catalog(catalog=cat_gal, params=params) + + assert cl_all.shape == (4, N_ELL_BINS) + # Effective ells share the binning math with the map path: bitwise-stable. + npt.assert_allclose( + ell_eff, + np.array([11.5, 20.0, 30.5, 43.5, 58.5, 75.5, 95.0, 116.5]), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_allclose( + cl_all[0], # EE + np.array( + [ + -8.480351557857007e-06, + 3.1313332030157945e-06, + -1.7823125201406074e-06, + 3.9616637779304446e-07, + -2.620506952433997e-07, + -3.110249362821697e-07, + -6.168246601737946e-07, + 4.523379773412076e-08, + ] + ), + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + npt.assert_allclose( + cl_all[1], # EB + np.array( + [ + 1.0851111241279543e-05, + -2.3080564471133627e-06, + -8.888571491169305e-07, + -3.261182548134984e-07, + 5.997911542118897e-07, + 1.0323861486573016e-07, + -2.120110521972082e-07, + 5.976985156531123e-09, + ] + ), + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + npt.assert_allclose( + cl_all[3], # BB + np.array( + [ + 1.1206499651601968e-05, + -3.114109101436142e-06, + 1.4336907374239745e-07, + -4.6001797302624655e-07, + 7.817610027040966e-07, + -4.18006575112908e-07, + -1.5043730205251898e-08, + -1.251687104580433e-07, + ] + ), + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + # BE (index 2) equals EB to within the catalog-path float noise. + npt.assert_allclose(cl_all[2], cl_all[1], rtol=RTOL_CAT, atol=ATOL_CAT) + + +# =========================================================================== +# apply_random_rotation -- invariant + reproducibility +# =========================================================================== +def test_apply_random_rotation_preserves_magnitude(cv, cat_and_params): + """The rotation conserves shape magnitude and array shape. + + The rotation is orthogonal, so per-object magnitude (e1^2 + e2^2) and array + shape are exactly preserved regardless of the random angle. + """ + cat_gal, params = cat_and_params + e1 = np.asarray(cat_gal[params["e1_col"]], dtype=np.float64) + e2 = np.asarray(cat_gal[params["e2_col"]], dtype=np.float64) + + e1_rot, e2_rot = cv.apply_random_rotation(e1, e2) + + assert e1_rot.shape == e1.shape + assert e2_rot.shape == e2.shape + # Rotation is orthogonal -> per-object magnitude is exactly preserved. + npt.assert_allclose(e1_rot**2 + e2_rot**2, e1**2 + e2**2, rtol=1e-12, atol=1e-15) + + +def test_apply_random_rotation_reproducible_with_seed(cv, cat_and_params): + """A seeded generator reproduces the rotation; rng=None stays random. + + The noise-debiasing used to call np.random.seed() with no argument, + re-seeding from entropy every call. apply_random_rotation now takes an + optional rng: the same seeded generator reproduces the draw, a different + seed differs, and rng=None remains non-deterministic. + """ + cat_gal, params = cat_and_params + e1 = np.asarray(cat_gal[params["e1_col"]], dtype=np.float64) + e2 = np.asarray(cat_gal[params["e2_col"]], dtype=np.float64) + + a1, a2 = cv.apply_random_rotation(e1, e2, np.random.default_rng(42)) + b1, b2 = cv.apply_random_rotation(e1, e2, np.random.default_rng(42)) + npt.assert_array_equal(a1, b1) + npt.assert_array_equal(a2, b2) + + c1, _ = cv.apply_random_rotation(e1, e2, np.random.default_rng(7)) + assert not np.allclose(a1, c1) + + d1, _ = cv.apply_random_rotation(e1, e2) + f1, _ = cv.apply_random_rotation(e1, e2) + assert not np.allclose(d1, f1) + + +# =========================================================================== +# calculate_pseudo_cl_catalog -- deterministic end-to-end catalog path +# =========================================================================== +def test_calculate_pseudo_cl_catalog_end_to_end(cv, tmp_path): + """End-to-end catalog path: FITS round-trip of ell + EE/EB/BB. + + The catalog method has no random noise debiasing, so it is reproducible to + the same ~2e-12 catalog-path float noise. save_pseudo_cl stores ELL/EE/EB/BB + (it drops the BE row); we pin the round-tripped table. + """ + ver = cv._test_version + cv._pseudo_cls = {ver: {}} + out_path = cv._output_path(f"pseudo_cl_cat_{ver}.fits") + cv.calculate_pseudo_cl_catalog(ver, out_path) + + assert os.path.exists(out_path) + d = fits.getdata(out_path) + # FITS gives big-endian f8; normalize for value comparison. + ell = np.asarray(d["ELL"], dtype=np.float64) + ee = np.asarray(d["EE"], dtype=np.float64) + eb = np.asarray(d["EB"], dtype=np.float64) + bb = np.asarray(d["BB"], dtype=np.float64) + + npt.assert_allclose( + ell, + np.array([11.5, 20.0, 30.5, 43.5, 58.5, 75.5, 95.0, 116.5]), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_allclose( + ee, + np.array( + [ + -8.480351557857063e-06, + 3.131333203015815e-06, + -1.7823125201406123e-06, + 3.961663777930465e-07, + -2.6205069524340103e-07, + -3.1102493628216637e-07, + -6.168246601737958e-07, + 4.523379773413286e-08, + ] + ), + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + npt.assert_allclose( + eb, + np.array( + [ + 1.0851111241279614e-05, + -2.3080564471133826e-06, + -8.888571491169258e-07, + -3.2611825481350053e-07, + 5.997911542118896e-07, + 1.0323861486573065e-07, + -2.1201105219720712e-07, + 5.976985156542028e-09, + ] + ), + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + npt.assert_allclose( + bb, + np.array( + [ + 1.1206499651602029e-05, + -3.114109101436161e-06, + 1.43369073742402e-07, + -4.600179730262474e-07, + 7.817610027040972e-07, + -4.180065751129079e-07, + -1.5043730205253605e-08, + -1.2516871045803156e-07, + ] + ), + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + # The end-to-end catalog EE matches the primitive get_pseudo_cls_catalog EE + # (same computation, FITS round-trip) -- consistency, not an independent pin. + cat_gal = fits.getdata(cv.cc[ver]["shear"]["path"]) + params = get_params_rho_tau(cv.cc[ver], survey=ver) + _, cl_prim, _ = cv.get_pseudo_cls_catalog(catalog=cat_gal, params=params) + npt.assert_allclose(ee, cl_prim[0], rtol=RTOL_CAT, atol=ATOL_CAT)