From ab93109124af173d00c5900d4b34659138bb977b Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Thu, 9 Jul 2026 10:18:33 +0200 Subject: [PATCH 1/8] Update the import in core comos_val of cs_util get_cosmo function. --- src/sp_validation/cosmo_val/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 5b88aacf..780fade1 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -9,13 +9,13 @@ import numpy as np import yaml from astropy.io import fits +from cs_util.cosmo import get_cosmo from shear_psf_leakage import run_object, run_scale from ..b_modes import ( _get_pte_from_scale_cut, find_conservative_scale_cut_key, ) -from ..cosmology import get_cosmo from ..statistics import chi2_and_pte from .catalog_characterization import CatalogCharacterizationMixin from .cosebis import CosebisMixin From 6eabbc6b466c944ffe33dc83de5a1f68a9a9a96e Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Thu, 16 Jul 2026 15:49:31 +0200 Subject: [PATCH 2/8] Add unit tests for the tomographic pseudo cls --- .../tests/data/test_cl_catalog.npz | Bin 0 -> 6903 bytes src/sp_validation/tests/test_pseudo_cl.py | 195 ++++++++++++++++-- 2 files changed, 181 insertions(+), 14 deletions(-) create mode 100644 src/sp_validation/tests/data/test_cl_catalog.npz diff --git a/src/sp_validation/tests/data/test_cl_catalog.npz b/src/sp_validation/tests/data/test_cl_catalog.npz new file mode 100644 index 0000000000000000000000000000000000000000..fdba323b63704311ff63bd9048ff9fad12c63c33 GIT binary patch literal 6903 zcmd5>3v?9K8J^8MP!NkCc!-9FB)n3BP&~pJmXJivhA|}|SRTVZW^;Gh%q%;z^duoB z6ha70wILJ0pgbBA{(ijo|HZ1*KSmYB5isig6VNlwY9_^cu?`5XeC&a`@D0nN+6>J`v+@|f|7 z#fkshr%!2Ab;nQxk~&!Q=X(oMc)R2iQu6&C#f6|-3pDjdyic%8qOAD*c15ew)M(vO zQ)5MkVN#`8m@;ja$)=TSbF^fOX$*)|qa8|tSHSiWPd%zttTS#!s-~41$8lLWMEb!p zL(yy%i?w1cSskQE`4Vq)iF_*Gb&77O@m^W*J0#xj(LjWHzZId85L%a%;&dspen2QS zH7@hXxw#lHS{Q=@u>_JO4qfy`xMJ{bd=7L!fC=2gC99Dc+1UtdMg~5lp+g!zOH0!f z3=Iqiebor(Z5v^GY=l7@u`p{jFrud|760myPxNR~s8<-;LaCa%@| z7L+Q>037}_vx#BMNd^_z!)OJ9Ak zbTsVv@r}cctzQPl6JOW0l^?4D!s178d;b1RA;I&DEF8IRoWQiel#yy^92nIS7Uip`17rb1~S!tykVDV%%DnY zO4q1MRK~mlva^`f81e|dt;R@R#U&wvF{l`NsW(oIM`e%N;W|d^bI)R-nx+RSYihjR zDoUcuZuRIW{5Cb3_lqt?yo$~K<@ekC-0S9B7U;_+zPA5kC*C0>|CN9UG?J?SH6G*06!$7 ziD;v?Nsmib z1X1wmD_CW_%cWKLZ9p{g36>)X@ByV;zCeM+1Pbv4SP6>3vm7V^rC?Rh=!3F5g#wi0 z$B_354y@LvE5sFoPZ1m%deMMj z^-4~sEGVdz(ywG!p@39~S8T>v94|Uug2y3i=>5#IdQpMk{rIg@=fEH&G(+x`vA*EJ@6~FZzvHb4);{4z>(h`daii59k|l1mE`(%>8?BCzEbr22 z!K+mz1GU{l*xEH4E12guz{90a*Ol(S#a{c#=tVcROH}jObBT`^b^K!=Jd~AN!mM7- zwzc-3zVq~3@WT^*I_H^B(w+de9@%MMTH<8wS@+wd!Ikh7J2$myLN45QvvB>pdp6RZ zfqf(Cv#cGqU2N{Dfk*%S`W5yyX2j6B+poe+)t9enoAYQ-0c#gdnmjwPojoukJ<~t+ z4E)`aH%k})Wg`2V&hy&9{iU?$U{|h=S5M4;9ahhLbm^|VpRq3<*_c{b^)x%bdHb@Q z`|@Z{0$+J^)W4j+5ZJoo**8<8%3*y{Lta8-4lKX&@%7wsOK8sm3z^ApzkFjMESs4X zF*a@>4~+zsES;exMUaD46If^T4O zd>6t6-_YRrI>H6tL*4KR-~)G*acy)0;7{NLI0@RoDexgU4L<6r;gh<)X!$>CxH}(X z=<$t7&~HJIjJw59jggYOrBIELn7h598Y4Y-FAvrDZj%1Vh_B%3iodcQ^X`4CxuKbi zE$oRuULAWn^(>Wi*m}IA?-Xz^eD7NMGtPIT$QZ*O9Fp)&Ovyb|(%BDM*VcXgQUh$- zY4?944<+MZSo%@I{x6HR1SGxsO69Hkiv_R|bk-+$&yw*B+_%}Jv~b%4l3vHO4@-Et znLTja|M->XR5EU0k3E#}x6c;zr;-kL700=ow#BgZ_333?{5U$oO`WyNsMh+EQ{XriclVxBpg!e~QBYe@C literal 0 HcmV?d00001 diff --git a/src/sp_validation/tests/test_pseudo_cl.py b/src/sp_validation/tests/test_pseudo_cl.py index cdaef9d4..61677827 100644 --- a/src/sp_validation/tests/test_pseudo_cl.py +++ b/src/sp_validation/tests/test_pseudo_cl.py @@ -51,6 +51,7 @@ """ import os +from pathlib import Path import numpy as np import numpy.testing as npt @@ -74,6 +75,8 @@ SEED = 1234 N_ELL_BINS = 8 +REFERENCE_TOMO = Path(__file__).parent / "data" / "test_cl_catalog.npz" + # --------------------------------------------------------------------------- # Synthetic-catalog fixture (deterministic; no cluster data) @@ -102,9 +105,10 @@ def _write_synthetic_config(tmp_path): 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 - ) + tomo_bin_id = rng.integers(1, 3, n_gal) # Create a two bin catalogue + Table( + {"RA": ra, "Dec": dec, "e1": e1, "e2": e2, "w": w, "tomo_bin_id": tomo_bin_id} + ).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)) @@ -121,6 +125,7 @@ def _write_synthetic_config(tmp_path): "e2_col_corrected": "e2", "ra_col": "RA", "dec_col": "Dec", + "tomo_bin_col": "tomo_bin_id", } # Minimal psf block so get_params_rho_tau() succeeds; the pseudo-Cl # primitives only read the shear-side keys, but the production call path @@ -181,6 +186,11 @@ def cat_and_params(cv): return cat_gal, params +def test_reference_exists(): + """Guard the guard: a missing reference must fail loudly, not skip.""" + assert REFERENCE_TOMO.exists(), f"committed reference missing: {REFERENCE_TOMO}" + + # Tolerances ---------------------------------------------------------------- # Bitwise-stable primitives (binning math, n_gal map, map-based pseudo-Cl). RTOL_DET = 1e-9 @@ -328,18 +338,11 @@ def _build_shear_map(cv, cat_gal, params): """Replicate calculate_pseudo_cl_map's weighted shear-map construction.""" unique_pix, _idx, idx_rep = cv.get_pixels(params, NSIDE, cat_gal) n_gal = cv.get_n_gal_map( - params, NSIDE, cat_gal, unique_pix=unique_pix, idx=None, idx_rep=idx_rep + params, NSIDE, cat_gal, unique_pix=unique_pix, idx=_idx, idx_rep=idx_rep + ) + m1, m2 = cv.get_shear_map( + params, NSIDE, cat_gal, unique_pix=unique_pix, idx=_idx, idx_rep=idx_rep ) - 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 @@ -347,6 +350,9 @@ 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) + ell_eff_2, cl_all_2, wsp_2 = cv.get_pseudo_cls_map( + shear_map, n_gal, shear_map_b=shear_map, mask_b=n_gal + ) assert cl_all.shape == (4, N_ELL_BINS) npt.assert_allclose( @@ -409,6 +415,80 @@ def test_get_pseudo_cls_map(cv, cat_and_params): # 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) + # Assert that running the same map against itself and against itself as a second map gives the same result. + npt.assert_allclose(cl_all, cl_all_2, rtol=RTOL_DET, atol=ATOL_DET) + npt.assert_allclose(ell_eff, ell_eff_2, rtol=RTOL_DET, atol=ATOL_DET) + + +def test_get_pseudo_cls_map_with_tomo(cv, cat_and_params): + cat_gal, params = cat_and_params + cat_gal_tomo1 = cat_gal[cat_gal[params["tomo_bin_col"]] == 1] + cat_gal_tomo2 = cat_gal[cat_gal[params["tomo_bin_col"]] == 2] + shear_map_1, n_gal_1 = _build_shear_map(cv, cat_gal_tomo1, params) + shear_map_2, n_gal_2 = _build_shear_map(cv, cat_gal_tomo2, params) + ell_eff, cl_all, wsp = cv.get_pseudo_cls_map( + shear_map_1, n_gal_1, shear_map_b=shear_map_2, mask_b=n_gal_2 + ) + + 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( + [ + -8.648947250571888e-06, + 4.386882223005456e-06, + -2.3193526251107696e-06, + 1.0651397314115168e-06, + -2.722591256637468e-07, + -3.862270431501105e-07, + -1.40420405782524e-06, + 2.1531429566476774e-07, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_allclose( + cl_all[1], # EB + np.array( + [ + -2.8032134416262087e-07, + -1.271298363785575e-06, + 1.7112863041724342e-08, + -2.899403854283714e-07, + 1.360217254538336e-06, + 7.682614612656511e-08, + 6.620365648170535e-07, + -3.164479179586416e-07, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_allclose( + cl_all[3], # BB + np.array( + [ + 9.387700364892905e-06, + -2.9045253374872504e-06, + -1.7006849005134948e-06, + -3.6869001010423035e-07, + 4.4986878855942184e-07, + 6.021221048891767e-07, + 8.413198154662088e-08, + -5.782957119801217e-07, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) + # =========================================================================== # get_pseudo_cls_catalog -- NmtFieldCatalog path (drifts ~2e-12) @@ -616,3 +696,90 @@ def test_calculate_pseudo_cl_catalog_end_to_end(cv, tmp_path): catalog=cat_gal, params=params, tomo_bin_a="all", tomo_bin_b="all" ) npt.assert_allclose(ee, cl_prim[0], rtol=RTOL_CAT, atol=ATOL_CAT) + + +def test_calculate_pseudo_cl_catalog_end_to_end_tomo(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: { + "tomo_bin_1_tomo_bin_1": {}, + "tomo_bin_1_tomo_bin_2": {}, + "tomo_bin_2_tomo_bin_2": {}, + } + } + out_path = cv._output_path(f"pseudo_cl_cat_{ver}.fits") + tomo_bin_ids, tomo_bin_pairs = cv._get_tomo_bins(ver) + + result_to_compare = np.load(REFERENCE_TOMO, allow_pickle=True)["arr_0"].item() + for tomo_bin_a, tomo_bin_b in tomo_bin_pairs: + out_path = cv._output_path( + f"pseudo_cl_cat_{ver}_{tomo_bin_a}_{tomo_bin_b}.fits" + ) + cv.calculate_pseudo_cl_catalog( + ver, out_path, tomo_bin_a=tomo_bin_a, tomo_bin_b=tomo_bin_b + ) + + 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) + be = np.asarray(d["BE"], dtype=np.float64) + bb = np.asarray(d["BB"], dtype=np.float64) + + npt.assert_allclose( + ell, + result_to_compare[f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ + "pseudo_cl" + ]["ELL"], + rtol=RTOL_DET, + atol=ATOL_DET, + ) + npt.assert_allclose( + ee, + result_to_compare[f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ + "pseudo_cl" + ]["EE"], + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + npt.assert_allclose( + eb, + result_to_compare[f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ + "pseudo_cl" + ]["EB"], + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + npt.assert_allclose( + be, + result_to_compare[f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ + "pseudo_cl" + ]["BE"], + rtol=RTOL_CAT, + atol=ATOL_CAT, + ) + npt.assert_allclose( + bb, + result_to_compare[f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ + "pseudo_cl" + ]["BB"], + 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, tomo_bin_a=tomo_bin_a, tomo_bin_b=tomo_bin_b + ) + npt.assert_allclose(ee, cl_prim[0], rtol=RTOL_CAT, atol=ATOL_CAT) From dee0141259278d70c370811326de8c6f23c60e60 Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Fri, 17 Jul 2026 08:46:43 +0200 Subject: [PATCH 3/8] Fix non-tomo branch bug in the _merge_iNKA_covariance_function --- src/sp_validation/cosmo_val/pseudo_cl.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index d8c1caf8..11c6572b 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -1079,7 +1079,7 @@ def _merge_iNKA_covariance(self, ver, tomography): else: block_path = self._output_path_iNKA_block_cov( - ver, tomo_bin_quad=(bin_key_a1, bin_key_a2, bin_key_b1, bin_key_b2) + ver, tomo_bin_quad=("all", "all", "all", "all") ) covar = fits.open(block_path) From 21863bc15c113a8caec2b758384b1dfbc62ea146 Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Fri, 17 Jul 2026 08:50:36 +0200 Subject: [PATCH 4/8] Fix small bugs identified by Fable --- src/sp_validation/cosmo_val/pseudo_cl.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index 11c6572b..13f747ad 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -223,7 +223,7 @@ def calculate_pseudo_cl_inka_cov( ) # Save in the dictionnaries - if not f"W{bin_key1}" not in n_gal_map_dict: + if f"W{bin_key1}" not in n_gal_map_dict: n_gal_map_dict[f"W{bin_key1}"] = n_gal_map_a if f"W{bin_key2}" not in n_gal_map_dict: n_gal_map_dict[f"W{bin_key2}"] = n_gal_map_b @@ -311,7 +311,7 @@ def calculate_pseudo_cl_inka_cov( ) input_cl_a1_b2 = ( fiducial_cl[f"W{bin_key_a1}xW{bin_key_b2}"] - if bin_key_a1 <= bin_key_a2 + if bin_key_a1 <= bin_key_b2 else fiducial_cl[f"W{bin_key_b2}xW{bin_key_a1}"] ) input_cl_a2_b1 = ( @@ -1024,15 +1024,13 @@ def _merge_iNKA_covariance(self, ver, tomography): ver, method="iNKA", tomography=tomography ) - print(self._pseudo_cls[ver]["tomo_bin_all_tomo_bin_all"]) - if tomography: # Merge the tomographic covariance matrices tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) # Get the number of bins from the pseudo_cls attribute # The non-tomographic pseudo-cl are computed from the call - # to this attribute. + # to this attribute if not already computed. n_ell = self._pseudo_cls[ver]["tomo_bin_all_tomo_bin_all"]["pseudo_cl"][ "ELL" ].shape[0] From 8996879d9d60dcf4548898fa125feba7b2683315 Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Fri, 17 Jul 2026 09:08:20 +0200 Subject: [PATCH 5/8] Update the cl tests accounting for Fable comments --- .../generate_test_cl_catalog_reference.py | 137 ++++++++++++++++++ src/sp_validation/tests/test_pseudo_cl.py | 45 +++++- 2 files changed, 181 insertions(+), 1 deletion(-) create mode 100644 src/sp_validation/tests/data/generate_test_cl_catalog_reference.py diff --git a/src/sp_validation/tests/data/generate_test_cl_catalog_reference.py b/src/sp_validation/tests/data/generate_test_cl_catalog_reference.py new file mode 100644 index 00000000..7c41d564 --- /dev/null +++ b/src/sp_validation/tests/data/generate_test_cl_catalog_reference.py @@ -0,0 +1,137 @@ +# %% +from pathlib import Path + +import IPython +import numpy as np +import yaml +from astropy.io import fits +from astropy.table import Table + +from sp_validation.cosmo_val import CosmologyValidation +from sp_validation.rho_tau import get_params_rho_tau + +ipython = IPython.get_ipython() + +# %% +NSIDE = 64 +SEED = 1234 +N_ELL_BINS = 8 + +rng = np.random.default_rng(SEED) +version = "TestCatalog" + +cat_dir = Path("./catalog") +nz_dir = Path("./nz") +output_dir = Path("./output_test") + +for d in (cat_dir, nz_dir, output_dir): + d.mkdir(exist_ok=True) + +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) +tomo_bin_id = rng.integers(1, 3, n_gal) # Create a two bin catalogue +Table( + {"RA": ra, "Dec": dec, "e1": e1, "e2": e2, "w": w, "tomo_bin_id": tomo_bin_id} +).write(cat_dir / "shear.fits", overwrite=True) + +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", + "tomo_bin_col": "tomo_bin_id", +} + +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 = Path("./config.yaml") +config_path.write_text(yaml.dump(config_data, sort_keys=False)) + +# %% +cv = CosmologyValidation( + versions=[version], + catalog_config=config_path, + output_dir=str(output_dir), + nside=NSIDE, + binning="powspace", + power=0.5, + n_ell_bins=N_ELL_BINS, + pol_factor=-1, +) +cv._test_version = version + +ver = cv._test_version +params = get_params_rho_tau(cv.cc[ver], survey=ver) +cat_gal = fits.getdata(cv.cc[ver]["shear"]["path"]) + +cat_gal_tomo_bin_1 = cat_gal[cat_gal[params["tomo_bin_col"]] == 1] +cat_gal_tomo_bin_2 = cat_gal[cat_gal[params["tomo_bin_col"]] == 2] +n_gal_map_a = cv.get_n_gal_map(params, NSIDE, cat_gal_tomo_bin_1) +shear_map_a_e1, shear_map_a_e2 = cv.get_shear_map(params, NSIDE, cat_gal_tomo_bin_1) +n_gal_map_b = cv.get_n_gal_map(params, NSIDE, cat_gal_tomo_bin_2) +shear_map_b_e1, shear_map_b_e2 = cv.get_shear_map(params, NSIDE, cat_gal_tomo_bin_2) + +shear_map_a = shear_map_a_e1 + 1j * shear_map_a_e2 +shear_map_b = shear_map_b_e1 + 1j * shear_map_b_e2 + + +# %% +ell_eff, cl_all, wsp = cv.get_pseudo_cls_map( + shear_map_a, n_gal_map_a, shear_map_b=shear_map_b, mask_b=n_gal_map_b +) + +# %% +# Get the catalogue output +ver = cv._test_version +cv._pseudo_cls = { + ver: { + "tomo_bin_1_tomo_bin_1": {}, + "tomo_bin_1_tomo_bin_2": {}, + "tomo_bin_2_tomo_bin_2": {}, + } +} +out_path = cv._output_path(f"pseudo_cl_cat_{ver}.fits") +tomo_bin_ids, tomo_bin_pairs = cv._get_tomo_bins(ver) +for tomo_bin_a, tomo_bin_b in tomo_bin_pairs: + out_path = cv._output_path(f"pseudo_cl_cat_{ver}_{tomo_bin_a}_{tomo_bin_b}.fits") + cv.calculate_pseudo_cl_catalog( + ver, out_path, tomo_bin_a=tomo_bin_a, tomo_bin_b=tomo_bin_b + ) + + +# %% +np.savez("./test_cl_catalog", cv._pseudo_cls[ver]) +# %% +test_result = np.load("./test_cl_catalog.npz", allow_pickle=True) diff --git a/src/sp_validation/tests/test_pseudo_cl.py b/src/sp_validation/tests/test_pseudo_cl.py index 61677827..c07afd00 100644 --- a/src/sp_validation/tests/test_pseudo_cl.py +++ b/src/sp_validation/tests/test_pseudo_cl.py @@ -340,10 +340,36 @@ def _build_shear_map(cv, cat_gal, params): n_gal = cv.get_n_gal_map( params, NSIDE, cat_gal, unique_pix=unique_pix, idx=_idx, idx_rep=idx_rep ) + 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_build_shear_map(cv, cat_and_params): + """Pin the weighted shear map construction from the synthetic catalog.""" + cat_gal, params = cat_and_params + unique_pix, _idx, idx_rep = cv.get_pixels(params, NSIDE, cat_gal) + n_gal = cv.get_n_gal_map( + params, NSIDE, cat_gal, unique_pix=unique_pix, idx=_idx, idx_rep=idx_rep + ) m1, m2 = cv.get_shear_map( params, NSIDE, cat_gal, unique_pix=unique_pix, idx=_idx, idx_rep=idx_rep ) - return m1 + 1j * m2, n_gal + + m = m1 + 1j * m2 + + m_replicate, n_gal_replicate = _build_shear_map(cv, cat_gal, params) + + npt.assert_allclose(m, m_replicate, rtol=RTOL_DET, atol=ATOL_DET) + npt.assert_allclose(n_gal, n_gal_replicate, rtol=RTOL_DET, atol=ATOL_DET) def test_get_pseudo_cls_map(cv, cat_and_params): @@ -471,6 +497,23 @@ def test_get_pseudo_cls_map_with_tomo(cv, cat_and_params): rtol=RTOL_DET, atol=ATOL_DET, ) + npt.assert_allclose( + cl_all[2], # BE + np.array( + [ + -3.83934537153214e-06, + 6.218772676611559e-06, + -5.542920986484232e-06, + 8.810792449470091e-07, + -3.2973293732012423e-07, + -1.3155110757359602e-08, + -1.420039280183438e-07, + 3.09464700816748e-07, + ] + ), + rtol=RTOL_DET, + atol=ATOL_DET, + ) npt.assert_allclose( cl_all[3], # BB np.array( From 9c948cd15ce92d5828a18152ebb09558d2e34bb9 Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Fri, 17 Jul 2026 09:27:46 +0200 Subject: [PATCH 6/8] Reapply "Tomographic pseudo-cl" This reverts commit f499ac8ddba01c8c5694d38805c75c104233fb0e. --- .gitignore | 4 +- cosmo_val/cat_config.yaml | 46 ++ src/sp_validation/cosmo_val/core.py | 51 +- src/sp_validation/cosmo_val/pseudo_cl.py | 4 - src/sp_validation/pseudo_cl.py | 858 ++++++++++++++++++++-- src/sp_validation/rho_tau.py | 4 + src/sp_validation/tests/test_pseudo_cl.py | 55 +- 7 files changed, 953 insertions(+), 69 deletions(-) diff --git a/.gitignore b/.gitignore index f8eec899..a568c1ca 100644 --- a/.gitignore +++ b/.gitignore @@ -197,4 +197,6 @@ papers/catalog/plots/*.pdf papers/cosmo_val/logs/ # Ignore scratch notebooks -scratch/*/*.ipynb \ No newline at end of file +scratch/*/*.ipynb +scratch/guerrini/work_notebooks +scratch/guerrini/launch_scripts \ No newline at end of file diff --git a/cosmo_val/cat_config.yaml b/cosmo_val/cat_config.yaml index cc1326df..ea36a81d 100644 --- a/cosmo_val/cat_config.yaml +++ b/cosmo_val/cat_config.yaml @@ -1259,3 +1259,49 @@ SP_v1.6.6: e1_col: e1 e2_col: e2 path: unions_shapepipe_star_2024_v1.6.a.fits + +GLASS_mock_validation: + subdir: /n09data/guerrini/glass_mock_test/results/ + pipeline: SP + colour: violet + getdist_colour: 0.0, 0.5, 1.0 + ls: dashed + marker: '*' + cov_th: + A: 2405.3892055695346 + n_e: 6.128201234871523 + n_psf: 0.752316232272063 + sigma_e: 0.379587601488189 + mask: /home/guerrini/sp_validation/cosmo_inference/data/mask/mask_map_footprint_nside_4096.fits + psf: + PSF_flag: FLAG_PSF_HSM + PSF_size: SIGMA_PSF_HSM + square_size: true + star_flag: FLAG_STAR_HSM + star_size: SIGMA_STAR_HSM + hdu: 1 + path: unions_shapepipe_psf_2024_v1.6.a.fits + ra_col: RA + dec_col: Dec + e1_PSF_col: E1_PSF_HSM + e1_star_col: E1_STAR_HSM + e2_PSF_col: E2_PSF_HSM + e2_star_col: E2_STAR_HSM + shear: + R: 1.0 + path: tomo_test_2_glass_sim_00001_1024.fits + redshift_path: /home/guerrini/sp_validation_cosmostat/config/glass_mock/test_data/redshift_distribution_tomo.txt + w_col: w + e1_col: e1 + e1_PSF_col: e1_PSF + e2_col: e2 + e2_PSF_col: e2_PSF + cols: RA,Dec + tomo_bin_col: tom_bin_id + star: + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + path: unions_shapepipe_star_2024_v1.6.a.fits + diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 48251bbb..066bf4a1 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -1,5 +1,6 @@ # %% import copy +import itertools import os import re from pathlib import Path @@ -7,7 +8,7 @@ import colorama import numpy as np import yaml -from cs_util.cosmo import get_cosmo +from astropy.io import fits from shear_psf_leakage import run_object, run_scale from ..b_modes import ( @@ -88,7 +89,7 @@ class CosmologyValidation( Number of ell bins for pseudo-C_ell analysis (used with binning='powspace'). ell_step : int, default 10 Bin width in ell for linear binning (used with binning='linear'). - pol_factor : bool, default True + pol_factor : int, default -1 Apply polarization correction factor in pseudo-C_ell calculations. nrandom_cell : int, default 10 Number of random realizations for C_ell error estimation. @@ -97,6 +98,10 @@ class CosmologyValidation( 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. + compute_tomography : bool, default False + Whether to compute tomographic correlation functions and pseudo-C_ell. + force_run : bool, default False + If True, forces re-computation of results even if cached outputs exist. Attributes ---------- @@ -214,7 +219,7 @@ def __init__( power=1 / 2, n_ell_bins=32, ell_step=10, - pol_factor=True, + pol_factor=-1, cell_method="map", noise_bias_method="analytic", fiducial_input_inka="coupled", @@ -223,6 +228,8 @@ def __init__( path_onecovariance=None, cosmo_params=None, blind=None, + compute_tomography=False, + force_run=False, ): self.rho_tau_method = rho_tau_method self.cov_estimate_method = cov_estimate_method @@ -243,7 +250,10 @@ def __init__( self.power = power self.n_ell_bins = n_ell_bins self.ell_step = ell_step + + assert pol_factor in (-1, 1), "The polarisatio factor must be -1 or 1." self.pol_factor = pol_factor + self.nrandom_cell = nrandom_cell self.cell_seed = cell_seed self.cell_method = cell_method @@ -252,6 +262,8 @@ def __init__( self.nside_mask = nside_mask self.path_onecovariance = path_onecovariance self.blind = blind + self.compute_tomography = compute_tomography + self.force_run = force_run assert self.cell_method in ["map", "catalog"], ( "cell_method must be 'map' or 'catalog'" @@ -636,3 +648,36 @@ def summarize_bmodes(self, fiducial_scale_cut=(12, 83), versions=None): print() return summary + + def _get_tomo_bins(self, version): + """ + Return the tomo_bin_ids for a given version. If the version does not have tomography, return None. + + Returns + ------- + tomo_bin_ids : list or None + List of unique tomographic bin IDs for the version, or None if no tomography is available + tomo_bin_pairs : list of tuples or None + List of unique pairs of tomographic bin IDs (including self-pairs) for the version, or None if no tomography is available + """ + if "tomo_bin_col" in self.cc[version]["shear"]: + self.print_cyan( + f"Extracting tomography information from version {version}." + ) + cat_gal = fits.getdata(self.cc[version]["shear"]["path"]) + tomo_bin = cat_gal[self.cc[version]["shear"]["tomo_bin_col"]] + tomo_bin_ids = np.unique(tomo_bin) + tomo_bin_ids = tomo_bin_ids[ + tomo_bin_ids > 0 + ] # Exclude zero or negative bins + self.print_cyan( + f"Found {len(tomo_bin_ids)} tomographic bins for version {version}: {tomo_bin_ids}." + ) + + tomo_bin_pairs = list( + itertools.combinations_with_replacement(tomo_bin_ids, 2) + ) + return tomo_bin_ids, tomo_bin_pairs + else: + self.print_cyan(f"Version {version} does not have tomography information.") + return None, None diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index 13f747ad..e2537630 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -1085,11 +1085,7 @@ def _merge_iNKA_covariance(self, ver, tomography): return covar - # ---------------- Plotting functions for pseudo-Cl's ---------------- # def plot_pseudo_cl( - self, - pol_list, - versions=None, ell_factor="ell", cov_type="iNKA", offset=0.15, diff --git a/src/sp_validation/pseudo_cl.py b/src/sp_validation/pseudo_cl.py index c9355ec9..3508e804 100644 --- a/src/sp_validation/pseudo_cl.py +++ b/src/sp_validation/pseudo_cl.py @@ -16,11 +16,13 @@ import healpy as hp import numpy as np import pymaster as nmt +from cs_util.cosmo import get_theo_c_ell # Lowest multipole retained by the pseudo-Cl estimators. LMIN = 8 +# ---------------------- Binning utility functions ---------------------- def pseudo_cl_geometry(nside): """Return ``(lmin, lmax, b_lmax)`` for the pseudo-Cl estimator at ``nside``. @@ -87,7 +89,38 @@ def make_namaster_bin( return b -def get_n_gal_map(nside, ra, dec, weights=None): +# ---------------------- Map computation utility functions ---------------------- +def get_pixels(ra, dec, nside): + """ + Get the HEALPix pixel indices for given RA and Dec. + + Parameters + ---------- + ra : np.ndarray + Right ascension in degrees. + dec : np.ndarray + Declination in degrees. + nside : int + HEALPix nside parameter. + + Returns + ------- + unique_pix : np.ndarray + Sorted unique 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. + """ + pixels = hp.ang2pix(nside, theta=np.radians(90 - dec), phi=np.radians(ra)) + + unique_pix, idx, idx_rep = np.unique(pixels, return_index=True, return_inverse=True) + return unique_pix, idx, idx_rep + + +def get_n_gal_map( + nside, ra, dec, weights=None, unique_pix=None, idx=None, idx_rep=None +): """Weighted galaxy number-density HEALPix map plus pixel bookkeeping. Bins ``(ra, dec)`` (degrees) onto an ``nside`` HEALPix grid. With @@ -98,21 +131,106 @@ def get_n_gal_map(nside, ra, dec, weights=None): ------- 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) + if unique_pix is None or idx is None or idx_rep is None: + unique_pix, idx, idx_rep = get_pixels(ra, dec, nside) - 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 + return n_gal + + +def get_shear_map( + ra, dec, e1, e2, w, nside, unique_pix=None, idx=None, idx_rep=None, n_gal_map=None +): + """Weighted shear HEALPix maps plus pixel bookkeeping. + + Bins ``(ra, dec)`` (degrees) onto an ``nside`` HEALPix grid. The shear + components ``(e1, e2)`` are weighted by ``w`` and summed per pixel. If + ``unique_pix``, ``idx``, and ``idx_rep`` are provided, they are used to + avoid recomputing the pixel indices. + If ``n_gal_map`` is provided, it is used to normalize the shear maps by the galaxy density. + + Returns + ------- + e1_map : np.ndarray + Weighted sum of E1 per pixel, shape ``(npix,)``. + e2_map : np.ndarray + Weighted sum of E2 per pixel, shape ``(npix,)``. + """ + if unique_pix is None or idx is None or idx_rep is None: + unique_pix, idx, idx_rep = get_pixels(ra, dec, nside) + + if n_gal_map is None: + n_gal_map = get_n_gal_map( + nside, ra, dec, weights=w, unique_pix=unique_pix, idx=idx, idx_rep=idx_rep + ) + + npix = hp.nside2npix(nside) + e1_map = np.zeros(npix) + e2_map = np.zeros(npix) + + e1_map[unique_pix] = np.bincount(idx_rep, weights=e1 * w) + e2_map[unique_pix] = np.bincount(idx_rep, weights=e2 * w) + + non_zero = n_gal_map > 0 + e1_map[non_zero] /= n_gal_map[non_zero] + e2_map[non_zero] /= n_gal_map[non_zero] + + return e1_map, e2_map + + +def get_variance_map( + nside, ra, dec, e1, e2, w, unique_pix=None, idx=None, idx_rep=None +): + """Compute the variance map of the shear components. + + The variance is computed as the weighted variance of the shear components in each pixel. + + Returns + ------- + variance_map : np.ndarray + Variance map of the shear components, shape ``(npix,)``. + """ + if unique_pix is None or idx is None or idx_rep is None: + unique_pix, idx, idx_rep = get_pixels(ra, dec, nside) + + npix = hp.nside2npix(nside) + variance_map = np.zeros(npix) + + variance_map[unique_pix] = np.bincount( + idx_rep, weights=0.5 * (e1**2 + e2**2) * w**2 + ) + + return variance_map + + +def get_noise_bias_analytical( + ra, dec, e1, e2, w, lmax, nside=1024, unique_pix=None, idx=None, idx_rep=None +): + """ + Compute the analytical noise bias for shear power spectrum. + """ + variance_map = get_variance_map( + nside=nside, + ra=ra, + dec=dec, + e1=e1, + e2=e2, + w=w, + unique_pix=unique_pix, + idx=idx, + idx_rep=idx_rep, + ) + + noise_bias = hp.nside2pixarea(nside) * np.mean(variance_map) + + noise_bias_cl = np.zeros((4, lmax)) + + noise_bias_cl[0, :] = noise_bias # EE + noise_bias_cl[3, :] = noise_bias # BB + + return noise_bias_cl def apply_random_rotation(e1, e2, rng=None): @@ -140,13 +258,493 @@ def apply_random_rotation(e1, e2, rng=None): return e1_out, e2_out +def get_noise_realisation( + ra, + dec, + e1, + e2, + w, + nside, + unique_pix=None, + idx=None, + idx_rep=None, + n_gal_map=None, + rng=None, +): + """ + Generate a random noise realisation of the shear maps by applying a random rotation to the ellipticity components. + + Parameters + ---------- + ra, dec : np.ndarray + Right ascension and declination of the sources. + e1, e2 : np.ndarray + Ellipticity components. + w : np.ndarray + Weights of the sources. + nside : int + HEALPix resolution. + unique_pix, idx, idx_rep : np.ndarray, optional + Pixel indices and bookkeeping arrays. If not provided, they will be computed. + n_gal_map : np.ndarray, optional + Galaxy number density map. If not provided, it will be computed. + + Returns + ------- + noise_map_e1, noise_map_e2 : np.ndarray + Noise map for ellipticity components. + """ + # Apply random rotation to the ellipticity components + e1_rot, e2_rot = apply_random_rotation(e1, e2, rng=rng) + + # Compute the noise maps using the rotated ellipticity components + noise_map_e1, noise_map_e2 = get_shear_map( + ra=ra, + dec=dec, + e1=e1_rot, + e2=e2_rot, + w=w, + nside=nside, + unique_pix=unique_pix, + idx=idx, + idx_rep=idx_rep, + n_gal_map=n_gal_map, + ) + + return noise_map_e1, noise_map_e2 + + +def get_noise_bias_from_gaussian_real( + ra, + dec, + e1, + e2, + w, + nside, + nrandom_cell, + binning, + ell_step=10, + n_ell_bins=32, + power=0.5, + unique_pix=None, + idx=None, + idx_rep=None, + n_gal_map=None, + wsp=None, + seed=42, +): + """ + Compute the power spectrum of the noise bias from random realisations. + + Parameters + ---------- + ra, dec : np.ndarray + Right ascension and declination of the sources. + e1, e2 : np.ndarray + Ellipticity components. + w : np.ndarray + Weights of the sources. + nside : int + HEALPix resolution. + nrandom_cell : int + Number of random cells to use for the noise estimation. + binning, ell_step, n_ell_bins, power : str, int, int, float + Binning scheme and parameters. + unique_pix, idx, idx_rep : np.ndarray, optional + Pixel indices and bookkeeping arrays. If not provided, they will be computed. + n_gal_map : np.ndarray, optional + Galaxy number density map. If not provided, it will be computed. + + Returns + ------- + noise_bias_cl : np.ndarray + Power spectrum of the noise bias. + """ + 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() + noise_bias_cl = np.zeros((4, ell_eff.size)) + + rng = np.random.default_rng(seed) + + if unique_pix is None or idx is None or idx_rep is None: + unique_pix, idx, idx_rep = get_pixels(ra, dec, nside) + + if n_gal_map is None: + n_gal_map = get_n_gal_map( + nside, ra, dec, weights=w, unique_pix=unique_pix, idx=idx, idx_rep=idx_rep + ) + + if wsp is None: + _, _, wsp = get_field_and_workspace_from_map(b, mask_a=n_gal_map) + + for _ in range(nrandom_cell): + noise_map_e1, noise_map_e2 = get_noise_realisation( + ra, + dec, + e1, + e2, + w, + nside, + unique_pix=unique_pix, + idx=idx, + idx_rep=idx_rep, + n_gal_map=n_gal_map, + rng=rng, + ) + + noise_map = noise_map_e1 + 1j * noise_map_e2 + del noise_map_e1, noise_map_e2 + + _, cl_noise_, _ = get_pseudo_cls_map( + noise_map, + n_gal_map, + nside, + binning, + ell_step=ell_step, + n_ell_bins=n_ell_bins, + power=power, + wsp=wsp, + ) + + noise_bias_cl += cl_noise_ + + noise_bias_cl /= nrandom_cell + + return noise_bias_cl + + +def get_noise_bias( + ra, + dec, + e1, + e2, + w, + nside, + noise_bias_method, + binning, + *, + ell_step=10, + n_ell_bins=32, + power=0.5, + nrandom_cell=100, + seed=42, +): + """ + Compute the noise bias from object positions and ellipticities. + + Parameters + ---------- + ra, dec : np.ndarray + Right ascension and declination of the sources. + e1, e2 : np.ndarray + Ellipticity components. + w : np.ndarray + Weights of the sources. + nside : int + HEALPix resolution. + noise_bias_method : {'randoms', 'analytic'} + Method used to estimate the noise bias. + 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. + nrandom_cell : int, optional + Number of random cells to use for the noise estimation. (only for the `randoms` method) + seed : int, optional + Random seed for reproducibility. (only for the `randoms` method) + + Returns + ------- + noise_bias_cl : np.ndarray + Power spectrum of the noise bias. + """ + if noise_bias_method not in ["randoms", "analytic"]: + raise ValueError("noise_bias_method must be 'randoms' or 'analytic'") + + 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, + ) + + unique_pix, idx, idx_rep = get_pixels(ra, dec, nside) + + if noise_bias_method == "analytic": + noise_bias_cl = get_noise_bias_analytical( + ra, + dec, + e1, + e2, + w, + lmax, + nside, + unique_pix=unique_pix, + idx=idx, + idx_rep=idx_rep, + ) + + elif noise_bias_method == "randoms": + noise_bias_cl = get_noise_bias_from_gaussian_real( + ra, + dec, + e1, + e2, + w, + nside, + nrandom_cell, + binning, + ell_step=ell_step, + n_ell_bins=n_ell_bins, + power=power, + unique_pix=unique_pix, + idx=idx, + idx_rep=idx_rep, + seed=seed, + ) + + noise_bias_cl = b.unbin_cell(noise_bias_cl) + else: + raise ValueError( + f"Invalid noise bias method `{noise_bias_method}`. Must be 'analytic' or 'randoms'." + ) + + return noise_bias_cl + + +# ---------------------- Cl computation functions ---------------------- +def get_field_and_workspace_from_map( + b, + mask_a, + e1_map_a=None, + e2_map_a=None, + mask_b=None, + e1_map_b=None, + e2_map_b=None, + pol_factor=-1, + return_wsp=True, +): + """Compute a NaMaster field and workspace object from the input maps. + + If the shear maps are None, returns field objects but only the workspace objects is relevant and contains the mixing matrix. + If the second mask and shear maps (indexed b) are provided, the mixing matrix is computed between the two fields. + + Parameters + ---------- + b : nmt.NmtBin + NaMaster binning object. + mask_a : np.ndarray + Field mask for the first map. + e1_map_a : np.ndarray, optional + E1 map for the first field. + e2_map_a : np.ndarray, optional + E2 map for the first field. + mask_b : np.ndarray, optional + Field mask for the second map. + e1_map_b : np.ndarray, optional + E1 map for the second field. + e2_map_b : np.ndarray, optional + E2 map for the second field. + pol_factor : float, optional + Polarization factor to apply to the E2 map. + return_wsp : bool, optional + If True, return the NaMaster workspace object containing the mixing matrix. + + Returns + ------- + field_a : nmt.NmtField + NaMaster field object for the first map. + field_b : nmt.NmtField + NaMaster field object for the second map (if provided, same than the first map otherwise). + wsp : nmt.NmtWorkspace + NaMaster workspace object containing the mixing matrix. + + """ + nside = hp.npix2nside(len(mask_a)) + lmax = b.lmax + if e1_map_a is None or e2_map_a is None: + e1_map_a = np.zeros(hp.nside2npix(nside)) + e2_map_a = np.zeros(hp.nside2npix(nside)) + + # Create NaMaster field + field_a = nmt.NmtField( + mask=mask_a, maps=[e1_map_a, pol_factor * e2_map_a], lmax=lmax + ) + + if mask_b is not None: + if e1_map_b is None or e2_map_b is None: + e1_map_b = np.zeros(hp.nside2npix(nside)) + e2_map_b = np.zeros(hp.nside2npix(nside)) + + field_b = nmt.NmtField( + mask=mask_b, maps=[e1_map_b, pol_factor * e2_map_b], lmax=lmax + ) + else: + field_b = field_a + + if return_wsp: + # Create NaMaster workspace + wsp = nmt.NmtWorkspace.from_fields(field_a, field_b, b) + + return field_a, field_b, wsp + else: + return field_a, field_b, None + + +def get_field_and_workspace_from_catalog( + b, + ra_a, + dec_a, + e1_a, + e2_a, + w_a, + ra_b=None, + dec_b=None, + e1_b=None, + e2_b=None, + w_b=None, + pol_factor=-1, + return_wsp=True, + same_bin=False, +): + """Create a NaMaster field and workspace from the input catalog. + + If the second catalog is provided, the mixing matrix is computed between the two fields. + + Parameters + ---------- + b : nmt.NmtBin + NaMaster binning object. + ra_a : np.ndarray + Right ascension of sources in the first catalog. + dec_a : np.ndarray + Declination of sources in the first catalog. + e1_a : np.ndarray + E1 shear component of sources in the first catalog. + e2_a : np.ndarray + E2 shear component of sources in the first catalog. + w_a : np.ndarray + Weights of sources in the first catalog. + ra_b : np.ndarray, optional + Right ascension of sources in the second catalog. + dec_b : np.ndarray, optional + Declination of sources in the second catalog. + e1_b : np.ndarray, optional + E1 shear component of sources in the second catalog. + e2_b : np.ndarray, optional + E2 shear component of sources in the second catalog. + w_b : np.ndarray, optional + Weights of sources in the second catalog. + pol_factor : float, optional + Polarization factor to apply to the E2 component. + return_wsp : bool, optional + If True, return the NaMaster workspace object containing the mixing matrix. + + Returns + ------- + field_a : nmt.NmtFieldCatalog + NaMaster field object for the first catalog. + field_b : nmt.NmtFieldCatalog + NaMaster field object for the second catalog (if provided, same as the first catalog otherwise). + wsp : nmt.NmtWorkspace + NaMaster workspace object containing the mixing matrix. + + """ + lmax = b.lmax + # Get field for input catalog a + field_a = nmt.NmtFieldCatalog( + positions=[ra_a, dec_a], + weights=w_a, + field=[e1_a, pol_factor * e2_a], + lmax=lmax, + lmax_mask=lmax, + spin=2, + lonlat=True, + ) + + if ( + ra_b is not None + and dec_b is not None + and e1_b is not None + and e2_b is not None + and w_b is not None + and not same_bin + ): + field_b = nmt.NmtFieldCatalog( + positions=[ra_b, dec_b], + weights=w_b, + field=[e1_b, pol_factor * e2_b], + lmax=lmax, + lmax_mask=lmax, + spin=2, + lonlat=True, + ) + else: + field_b = field_a + + if return_wsp: + wsp = nmt.NmtWorkspace.from_fields(field_a, field_b, b) + return field_a, field_b, wsp + else: + return field_a, field_b, None + + +def compute_cl_from_field_and_workspace(field_a, field_b, wsp, b): + """Compute the angular power spectrum from the input NaMaster field and workspace + + Parameters + ---------- + field_a : nmt.NmtField + NaMaster field object for the first catalog. + field_b : nmt.NmtField + NaMaster field object for the second catalog. + wsp : nmt.NmtWorkspace + NaMaster workspace object containing the mixing matrix. + b : nmt.NmtBin + NaMaster binning object. + + Returns + ------- + cl_coupled : np.ndarray + Coupled angular power spectrum. + cl_decoupled : np.ndarray + Decoupled angular power spectrum. + """ + cl_coupled = nmt.compute_coupled_cell(field_a, field_b) + cl_decoupled = wsp.decouple_cell(cl_coupled) + + return cl_coupled, cl_decoupled + + def get_pseudo_cls_map( - shear_map, - mask, + shear_map_a, + mask_a, nside, binning, *, - pol_factor=True, + shear_map_b=None, + mask_b=None, + pol_factor=-1, wsp=None, ell_step=10, n_ell_bins=32, @@ -156,16 +754,20 @@ def get_pseudo_cls_map( Parameters ---------- - shear_map : np.ndarray + shear_map_a : np.ndarray Complex shear map (``e1 + 1j * e2``). - mask : np.ndarray + mask_a : 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. + shear_map_b : np.ndarray, optional + Complex shear map for the second field (``e1 + 1j * e2``). + mask_b : np.ndarray, optional + Field mask for the second field (the galaxy number-density map). + pol_factor : float, optional + Polarization factor to apply to the E2 component. wsp : nmt.NmtWorkspace, optional Reuse a coupling workspace; built from the field if ``None``. ell_step, n_ell_bins, power : optional @@ -180,6 +782,27 @@ def get_pseudo_cls_map( wsp : nmt.NmtWorkspace The coupling workspace (newly built or the one passed in). """ + # First do some assertion checks + if shear_map_b is not None: + assert mask_b is not None, "mask_b must be provided if shear_map_b is provided" + assert shear_map_a.shape == shear_map_b.shape, ( + "shear_map_a and shear_map_b must have the same shape" + ) + assert mask_a.shape == mask_b.shape, ( + "mask_a and mask_b must have the same shape" + ) + + if mask_b is not None: + assert shear_map_b is not None, ( + "shear_map_b must be provided if mask_b is provided" + ) + assert shear_map_a.shape == shear_map_b.shape, ( + "shear_map_a and shear_map_b must have the same shape" + ) + assert mask_a.shape == mask_b.shape, ( + "mask_a and mask_b must have the same shape" + ) + lmin, lmax, b_lmax = pseudo_cl_geometry(nside) b = make_namaster_bin( @@ -193,18 +816,36 @@ def get_pseudo_cls_map( ) 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) + field_a, field_b, wsp = get_field_and_workspace_from_map( + b, + mask_a, + e1_map_a=shear_map_a.real, + e2_map_a=shear_map_a.imag, + mask_b=mask_b, + e1_map_b=shear_map_b.real if shear_map_b is not None else None, + e2_map_b=shear_map_b.imag if shear_map_b is not None else None, + pol_factor=pol_factor, + return_wsp=True, + ) + else: + field_a, field_b, _ = get_field_and_workspace_from_map( + b, + mask_a, + e1_map_a=shear_map_a.real, + e2_map_a=shear_map_a.imag, + mask_b=mask_b, + e1_map_b=shear_map_b.real if shear_map_b is not None else None, + e2_map_b=shear_map_b.imag if shear_map_b is not None else None, + pol_factor=pol_factor, + return_wsp=False, + ) - cl_coupled = nmt.compute_coupled_cell(f_all, f_all) - cl_all = wsp.decouple_cell(cl_coupled) + cl_coupled, cl_decoupled = compute_cl_from_field_and_workspace( + field_a, field_b, wsp, b + ) - return ell_eff, cl_all, wsp + return ell_eff, cl_decoupled, wsp def get_pseudo_cls_catalog( @@ -213,7 +854,9 @@ def get_pseudo_cls_catalog( nside, binning, *, - pol_factor=True, + tomo_bin_a="all", + tomo_bin_b="all", + pol_factor=-1, wsp=None, ell_step=10, n_ell_bins=32, @@ -232,8 +875,10 @@ def get_pseudo_cls_catalog( 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. + tomo_bin_a, tomo_bin_b : str or int or None, optional + Tomographic bin IDs for the two fields. + pol_factor : int, optional + Polarization factor to apply to the E2 component. wsp : nmt.NmtWorkspace, optional Reuse a coupling workspace; built from the field if ``None``. ell_step, n_ell_bins, power : optional @@ -248,6 +893,11 @@ def get_pseudo_cls_catalog( wsp : nmt.NmtWorkspace The coupling workspace (newly built or the one passed in). """ + # First make some assertion checks reagarding the run mode + assert (tomo_bin_a == "all" and tomo_bin_b == "all") or ( + tomo_bin_a != "all" and tomo_bin_b != "all" + ), "Both tomo_bin_a and tomo_bin_b must be provided or both must be 'all'" + lmin, lmax, b_lmax = pseudo_cl_geometry(nside) b = make_namaster_bin( @@ -261,22 +911,140 @@ def get_pseudo_cls_catalog( ) ell_eff = b.get_effective_ells() - factor = -1 if pol_factor else 1 + is_tomography = tomo_bin_a != "all" and tomo_bin_b != "all" + if is_tomography: + assert params["tomo_bin_col"] is not None, ( + "The column of tomographic bin ids is not specified." + ) + mask_tomo_a = catalog[params["tomo_bin_col"]] == tomo_bin_a + mask_tomo_b = catalog[params["tomo_bin_col"]] == tomo_bin_b + catalog_a = catalog[mask_tomo_a] + catalog_b = catalog[mask_tomo_b] + same_bin = tomo_bin_a == tomo_bin_b + else: + catalog_a = catalog + catalog_b = catalog + same_bin = True - 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: + field_a, field_b, wsp = get_field_and_workspace_from_catalog( + b, + ra_a=catalog_a[params["ra_col"]], + dec_a=catalog_a[params["dec_col"]], + e1_a=catalog_a[params["e1_col"]], + e2_a=catalog_a[params["e2_col"]], + w_a=catalog_a[params["w_col"]], + ra_b=catalog_b[params["ra_col"]], + dec_b=catalog_b[params["dec_col"]], + e1_b=catalog_b[params["e1_col"]], + e2_b=catalog_b[params["e2_col"]], + w_b=catalog_b[params["w_col"]], + pol_factor=pol_factor, + return_wsp=True, + same_bin=same_bin, + ) + else: + field_a, field_b, _ = get_field_and_workspace_from_catalog( + b, + ra_a=catalog_a[params["ra_col"]], + dec_a=catalog_a[params["dec_col"]], + e1_a=catalog_a[params["e1_col"]], + e2_a=catalog_a[params["e2_col"]], + w_a=catalog_a[params["w_col"]], + ra_b=catalog_b[params["ra_col"]], + dec_b=catalog_b[params["dec_col"]], + e1_b=catalog_b[params["e1_col"]], + e2_b=catalog_b[params["e2_col"]], + w_b=catalog_b[params["w_col"]], + pol_factor=pol_factor, + return_wsp=False, + same_bin=same_bin, + ) + + cl_coupled, cl_decoupled = compute_cl_from_field_and_workspace( + field_a, field_b, wsp, b ) - if wsp is None: - wsp = nmt.NmtWorkspace.from_fields(f_all, f_all, b) + return ell_eff, cl_decoupled, wsp + + +# ---------------------- Covariance computation functions ---------------------- +def get_fiducial_cl(z, dndz, lmax, cosmo, backend="camb"): + """ + Get the fiducial Cl's using the redshift distribution. + Cosmology is determined by the input cosmo object. + """ + ell = np.arange(1, lmax + 1) + + fiducial_cl = get_theo_c_ell(ell=ell, z=z, nz=dndz, backend=backend, cosmo=cosmo) + + return fiducial_cl + + +def get_pseudo_cl_iNKA_covariance( + input_cl_a1_b1, + input_cl_a1_b2, + input_cl_a2_b1, + input_cl_a2_b2, + field_a1, + field_a2, + field_b1, + field_b2, + wsp_a, + wsp_b, + b, +): + """Compute the iNKA covariance for pseudo-Cl. + + Parameters + ---------- + input_cl_a1_b1 : np.ndarray + Input Cl for field a1 and b1. + input_cl_a1_b2 : np.ndarray + Input Cl for field a1 and b2. + input_cl_a2_b1 : np.ndarray + Input Cl for field a2 and b1. + input_cl_a2_b2 : np.ndarray + Input Cl for field a2 and b2. + field_a1 : nmt.NmtField + NaMaster field object for the first catalog (a1). + field_a2 : nmt.NmtField + NaMaster field object for the second catalog (a2). + field_b1 : nmt.NmtField + NaMaster field object for the first catalog (b1). + field_b2 : nmt.NmtField + NaMaster field object for the second catalog (b2). + wsp_a : nmt.NmtWorkspace + NaMaster workspace object containing the mixing matrix for fields a. + wsp_b : nmt.NmtWorkspace + NaMaster workspace object containing the mixing matrix for fields b. + b : nmt.NmtBin + NaMaster binning object. + + Returns + ------- + cov_matrix : np.ndarray + Covariance matrix of the pseudo-Cl, shape ``(n_bins, n_bins)``. + """ + # Compute the coupling coefficients for the covariance + cw = nmt.NmtCovarianceWorkspace.from_fields(field_a1, field_a2, field_b1, field_b2) + + # Get actual number of ell bins from binning scheme + n_ell_actual = b.get_n_bands() - cl_coupled = nmt.compute_coupled_cell(f_all, f_all) - cl_all = wsp.decouple_cell(cl_coupled) + # Compute the covariance using NaMaster's built-in function + cov_matrix = nmt.gaussian_covariance( + cw, + 2, + 2, + 2, + 2, + input_cl_a1_b1, + input_cl_a1_b2, + input_cl_a2_b1, + input_cl_a2_b2, + wsp_a, + wb=wsp_b, + ).reshape([n_ell_actual, 4, n_ell_actual, 4]) - return ell_eff, cl_all, wsp + return cov_matrix diff --git a/src/sp_validation/rho_tau.py b/src/sp_validation/rho_tau.py index 6874a08c..95a720f9 100644 --- a/src/sp_validation/rho_tau.py +++ b/src/sp_validation/rho_tau.py @@ -51,6 +51,10 @@ def get_params_rho_tau(cat, survey="other"): params["w_col"] = cat["shear"]["w_col"] params["e1_col"] = cat["shear"]["e1_col"] params["e2_col"] = cat["shear"]["e2_col"] + try: + params["tomo_bin_col"] = cat["shear"]["tomo_bin_col"] + except KeyError: + params["tomo_bin_col"] = None params["R11"] = cat["shear"].get("R11") params["R22"] = cat["shear"].get("R22") diff --git a/src/sp_validation/tests/test_pseudo_cl.py b/src/sp_validation/tests/test_pseudo_cl.py index 0fed63d1..c07afd00 100644 --- a/src/sp_validation/tests/test_pseudo_cl.py +++ b/src/sp_validation/tests/test_pseudo_cl.py @@ -59,6 +59,7 @@ import yaml from sp_validation.cosmo_val import CosmologyValidation +from sp_validation.pseudo_cl import apply_random_rotation from sp_validation.rho_tau import get_params_rho_tau # These tests need the full harmonic-space stack (pymaster/NaMaster + healpy), @@ -170,7 +171,7 @@ def cv(tmp_path): binning="powspace", power=0.5, n_ell_bins=N_ELL_BINS, - pol_factor=True, + pol_factor=-1, ) cv._test_version = version return cv @@ -275,18 +276,38 @@ def test_unknown_binning_raises(self, cv): cv.get_namaster_bin(LMIN, LMAX, B_LMAX) +# ========================================================================== +# get_pixels -- fetch the pixels from ra, dec, and nside (HEALPix) +# ========================================================================== +def test_get_pixels(cv, cat_and_params): + """Pin the HEALPix pixelization of the synthetic catalog.""" + cat_gal, params = cat_and_params + unique_pix, idx, idx_rep = cv.get_pixels(params, NSIDE, cat_gal) + + # Structural invariants of the HEALPix occupancy map. + assert unique_pix.size == 485 + assert idx_rep.size == 5000 # one entry per galaxy + + # Pinned scalar summaries: total weight is conserved (sum of weights), + # peak occupancy, and the pixel index bookkeeping. + npt.assert_array_equal(unique_pix[:5], [12167, 12168, 12169, 12170, 12171]) + assert int(unique_pix.sum()) == 7849989 + + # =========================================================================== # 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) + + unique_pix, idx, idx_rep = cv.get_pixels(params, NSIDE, cat_gal) + n_gal = cv.get_n_gal_map( + params, NSIDE, cat_gal, unique_pix=unique_pix, idx=None, idx_rep=idx_rep + ) # 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)) @@ -294,8 +315,6 @@ def test_get_n_gal_map(cv, cat_and_params): # 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( @@ -519,7 +538,9 @@ def test_get_pseudo_cls_map_with_tomo(cv, cat_and_params): # =========================================================================== 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) + ell_eff, cl_all, wsp = cv.get_pseudo_cls_catalog( + catalog=cat_gal, params=params, tomo_bin_a="all", tomo_bin_b="all" + ) assert cl_all.shape == (4, N_ELL_BINS) # Effective ells share the binning math with the map path: bitwise-stable. @@ -597,7 +618,7 @@ def test_apply_random_rotation_preserves_magnitude(cv, 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) + e1_rot, e2_rot = apply_random_rotation(e1, e2) assert e1_rot.shape == e1.shape assert e2_rot.shape == e2.shape @@ -617,16 +638,16 @@ def test_apply_random_rotation_reproducible_with_seed(cv, 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)) + a1, a2 = apply_random_rotation(e1, e2, np.random.default_rng(42)) + b1, b2 = 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)) + c1, _ = 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) + d1, _ = apply_random_rotation(e1, e2) + f1, _ = apply_random_rotation(e1, e2) assert not np.allclose(d1, f1) @@ -641,9 +662,9 @@ def test_calculate_pseudo_cl_catalog_end_to_end(cv, tmp_path): (it drops the BE row); we pin the round-tripped table. """ ver = cv._test_version - cv._pseudo_cls = {ver: {}} + cv._pseudo_cls = {ver: {"tomo_bin_all_tomo_bin_all": {}}} out_path = cv._output_path(f"pseudo_cl_cat_{ver}.fits") - cv.calculate_pseudo_cl_catalog(ver, out_path) + cv.calculate_pseudo_cl_catalog(ver, out_path, tomo_bin_a="all", tomo_bin_b="all") assert os.path.exists(out_path) d = fits.getdata(out_path) @@ -714,7 +735,9 @@ def test_calculate_pseudo_cl_catalog_end_to_end(cv, tmp_path): # (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) + _, cl_prim, _ = cv.get_pseudo_cls_catalog( + catalog=cat_gal, params=params, tomo_bin_a="all", tomo_bin_b="all" + ) npt.assert_allclose(ee, cl_prim[0], rtol=RTOL_CAT, atol=ATOL_CAT) From 3146790653076f9da1cb5fb74b74a069c4b1e8e5 Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Fri, 17 Jul 2026 09:30:28 +0200 Subject: [PATCH 7/8] Fix Ruff errors --- src/sp_validation/cosmo_val/pseudo_cl.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index e2537630..13f747ad 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -1085,7 +1085,11 @@ def _merge_iNKA_covariance(self, ver, tomography): return covar + # ---------------- Plotting functions for pseudo-Cl's ---------------- # def plot_pseudo_cl( + self, + pol_list, + versions=None, ell_factor="ell", cov_type="iNKA", offset=0.15, From 7eca12e874b06b0f6045036cb285c0fdcb1d1e6a Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Fri, 17 Jul 2026 09:35:19 +0200 Subject: [PATCH 8/8] Add get_cosmo import in core CosmoValidation that disappeared in the merging process. --- src/sp_validation/cosmo_val/core.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 066bf4a1..c6204c65 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -9,6 +9,7 @@ import numpy as np import yaml from astropy.io import fits +from cs_util.cosmo import get_cosmo from shear_psf_leakage import run_object, run_scale from ..b_modes import (