diff --git a/cosmo_val/cat_config.yaml b/cosmo_val/cat_config.yaml index ea36a81d..7fc6122f 100644 --- a/cosmo_val/cat_config.yaml +++ b/cosmo_val/cat_config.yaml @@ -1272,7 +1272,7 @@ GLASS_mock_validation: 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 + mask: /n09data/guerrini/glass_mock_test/mask_nside_1024.fits psf: PSF_flag: FLAG_PSF_HSM PSF_size: SIGMA_PSF_HSM @@ -1280,7 +1280,7 @@ GLASS_mock_validation: star_flag: FLAG_STAR_HSM star_size: SIGMA_STAR_HSM hdu: 1 - path: unions_shapepipe_psf_2024_v1.6.a.fits + path: ../../../../n17data/UNIONS/WL/v1.6.x/unions_shapepipe_psf_2024_v1.6.a.fits ra_col: RA dec_col: Dec e1_PSF_col: E1_PSF_HSM @@ -1303,5 +1303,5 @@ GLASS_mock_validation: dec_col: Dec e1_col: e1 e2_col: e2 - path: unions_shapepipe_star_2024_v1.6.a.fits + path: ../../../../n17data/UNIONS/WL/v1.6.x/unions_shapepipe_star_2024_v1.6.a.fits diff --git a/src/sp_validation/cosmo_val/catalog_characterization.py b/src/sp_validation/cosmo_val/catalog_characterization.py index bb5f503d..e8117011 100644 --- a/src/sp_validation/cosmo_val/catalog_characterization.py +++ b/src/sp_validation/cosmo_val/catalog_characterization.py @@ -150,13 +150,17 @@ def area(self): @property def n_eff_gal(self): if not hasattr(self, "_n_eff_gal"): - self.calculate_n_eff_gal() + self.calculate_n_eff_gal(tomography=False) + if self.compute_tomography: + self.calculate_n_eff_gal(tomography=True) return self._n_eff_gal @property def ellipticity_dispersion(self): if not hasattr(self, "_ellipticity_dispersion"): - self.calculate_ellipticity_dispersion() + self.calculate_ellipticity_dispersion(tomography=False) + if self.compute_tomography: + self.calculate_ellipticity_dispersion(tomography=True) return self._ellipticity_dispersion def _get_binned_catalog_mask(self, ver): @@ -202,29 +206,81 @@ def calculate_area_from_binned_catalog(self, ver): return area - def calculate_n_eff_gal(self): + def calculate_n_eff_gal(self, tomography=False): self.print_start("Calculating effective number of galaxy") - n_eff_gal = {} + if not hasattr(self, "_n_eff_gal"): + self._n_eff_gal = {} for ver in self.versions: self.print_magenta(ver) + if ver not in self._n_eff_gal: + self._n_eff_gal[ver] = {} + + if tomography: + tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) + + if tomo_bin_ids is None or tomo_bin_pairs is None: + raise ValueError( + f"Version {ver} does not have tomography information." + ) + + else: + tomo_bin_ids, tomo_bin_pairs = ["all"], [("all", "all")] + with self.results[ver].temporarily_read_data(): w = self._read_shear_cols(ver, "w_col") - n_eff_gal[ver] = n_eff_density(w, self.area[ver]) - print(f"n_eff_gal = {n_eff_gal[ver]:.2f} gal./arcmin^-2") + for tomo_bin_id in tomo_bin_ids: + if tomo_bin_id == "all": + self._n_eff_gal[ver][f"tomo_bin_{tomo_bin_id}"] = n_eff_density( + w, self.area[ver] + ) + else: + tomo_bin = self._read_shear_cols(ver, "tomo_bin_col") + mask = tomo_bin == tomo_bin_id + self._n_eff_gal[ver][f"tomo_bin_{tomo_bin_id}"] = n_eff_density( + w[mask], self.area[ver] + ) + print( + f"n_eff_gal for tomo_bin_{tomo_bin_id} = {self._n_eff_gal[ver][f'tomo_bin_{tomo_bin_id}']:.2f} gal./arcmin^-2" + ) - self._n_eff_gal = n_eff_gal self.print_done("Effective number of galaxy calculation finished") - def calculate_ellipticity_dispersion(self): + def calculate_ellipticity_dispersion(self, tomography=False): self.print_start("Calculating ellipticity dispersion") - ellipticity_dispersion = {} + if not hasattr(self, "_ellipticity_dispersion"): + self._ellipticity_dispersion = {} for ver in self.versions: self.print_magenta(ver) + if ver not in self._ellipticity_dispersion: + self._ellipticity_dispersion[ver] = {} + + if tomography: + tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) + + if tomo_bin_ids is None or tomo_bin_pairs is None: + raise ValueError( + f"Version {ver} does not have tomography information." + ) + + else: + tomo_bin_ids, tomo_bin_pairs = ["all"], [("all", "all")] + with self.results[ver].temporarily_read_data(): e1, e2, w = self._read_shear_cols(ver, "e1_col", "e2_col", "w_col") - ellipticity_dispersion[ver] = ellipticity_dispersion_stat(e1, e2, w) - print(f"Ellipticity dispersion = {ellipticity_dispersion[ver]:.4f}") - self._ellipticity_dispersion = ellipticity_dispersion + for tomo_bin_id in tomo_bin_ids: + if tomo_bin_id == "all": + self._ellipticity_dispersion[ver][f"tomo_bin_{tomo_bin_id}"] = ( + ellipticity_dispersion_stat(e1, e2, w) + ) + else: + tomo_bin = self._read_shear_cols(ver, "tomo_bin_col") + mask = tomo_bin == tomo_bin_id + self._ellipticity_dispersion[ver][f"tomo_bin_{tomo_bin_id}"] = ( + ellipticity_dispersion_stat(e1[mask], e2[mask], w[mask]) + ) + print( + f"Ellipticity dispersion for tomo_bin_{tomo_bin_id} = {self._ellipticity_dispersion[ver][f'tomo_bin_{tomo_bin_id}']:.4f}" + ) def plot_footprints(self): self.print_start("Plotting footprints:") diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index b0351d07..a0ceae5d 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -43,11 +43,199 @@ def pseudo_cls(self): @property def pseudo_cls_onecov(self): if not hasattr(self, "_pseudo_cls_onecov"): - self.calculate_pseudo_cl_onecovariance() + self.calculate_pseudo_cl_onecovariance(tomography=False) + if self.compute_tomography: + self.calculate_pseudo_cl_onecovariance(tomography=True) return self._pseudo_cls_onecov # ---------------- Pseudo-Cl calculation methods ---------------- # - # TODO: some cleaning to clearly separate DV, covariance, and utility functions. + def calculate_pseudo_cl(self, compute_tomography=True): + """ + Compute the pseudo-Cl of a `CosmologyValidation` inputs with tomography. + """ + out_dir = self._output_path("pseudo_cl") + os.makedirs(out_dir, exist_ok=True) + + if compute_tomography: + self.print_start("Computing tomographic pseudo-Cl's") + else: + self.print_start("Computing non-tomographic pseudo-Cl's") + + self._pseudo_cls = getattr(self, "_pseudo_cls", {}) + + for ver in self.versions: + self.print_magenta(ver) + + if ver not in self.pseudo_cls.keys(): + self._pseudo_cls[ver] = {} + + if compute_tomography: + tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) + + if tomo_bin_ids is None or tomo_bin_pairs is None: + raise ValueError( + f"Version {ver} does not have tomography information." + ) + + else: + tomo_bin_pairs = [("all", "all")] + + # Loop on the different tomographic bin pairs + for bin_key1, bin_key2 in tomo_bin_pairs: + self.print_cyan(f"Tomo Bin Pair: ({bin_key1}, {bin_key2})") + + if ( + f"tomo_bin_{bin_key1}_tomo_bin_{bin_key2}" + not in self._pseudo_cls[ver].keys() + ): + self._pseudo_cls[ver][ + f"tomo_bin_{bin_key1}_tomo_bin_{bin_key2}" + ] = {} + + out_path = self._output_path_pseudo_cl( + ver, tomo_bin_pair=(bin_key1, bin_key2) + ) + if os.path.exists(out_path) and not self.force_run: + self.print_done( + f"Skipping Pseudo-Cl's calculation, {out_path} exists" + ) + cl_shear = fits.getdata(out_path) + self._pseudo_cls[ver][f"tomo_bin_{bin_key1}_tomo_bin_{bin_key2}"][ + "pseudo_cl" + ] = cl_shear + continue + + if self.cell_method == "map": + self.calculate_pseudo_cl_map( + ver, self.nside, out_path, bin_key1, bin_key2 + ) + elif self.cell_method == "catalog": + self.calculate_pseudo_cl_catalog(ver, out_path, bin_key1, bin_key2) + else: + raise ValueError(f"Unknown cell method: {self.cell_method}") + + def calculate_pseudo_cl_map(self, ver, nside, out_path, tomo_bin_a, tomo_bin_b): + assert (tomo_bin_a == "all" and tomo_bin_b == "all") or ( + isinstance(tomo_bin_a, (int, np.integer)) + and isinstance(tomo_bin_b, (int, np.integer)) + ), "tomo_bin_a and tomo_bin_b must be either both 'all' or both integers." + + params = get_params_rho_tau(self.cc[ver]) + + self.print_cyan( + f"Computing pseudo-Cl's for tomographic bins {tomo_bin_a} and {tomo_bin_b}..." + ) + + # Load data and create shear and noise maps + cat_gal = fits.getdata(self.cc[ver]["shear"]["path"]) + + # Get the tomographic bin + cat_gal_a = self._get_tomographic_bin(params, cat_gal, tomo_bin_a) + cat_gal_b = self._get_tomographic_bin(params, cat_gal, tomo_bin_b) + + del cat_gal + + self.print_cyan("Creating maps and computing Cl's...") + # Get the pixels and indices for the catalogs + unique_pix_a, idx_a, idx_rep_a = self.get_pixels(params, nside, cat_gal_a) + unique_pix_b, idx_b, idx_rep_b = self.get_pixels(params, nside, cat_gal_b) + + # Create number density maps for each tomographic bin + n_gal_map_a = self.get_n_gal_map( + params, + nside, + cat_gal_a, + unique_pix=unique_pix_a, + idx=idx_a, + idx_rep=idx_rep_a, + ) + n_gal_map_b = self.get_n_gal_map( + params, + nside, + cat_gal_b, + unique_pix=unique_pix_b, + idx=idx_b, + idx_rep=idx_rep_b, + ) + + # Create shear maps for each tomographic bin + shear_map_a_e1, shear_map_a_e2 = self.get_shear_map( + params, + nside, + cat_gal_a, + unique_pix=unique_pix_a, + idx=idx_a, + idx_rep=idx_rep_a, + n_gal_map=n_gal_map_a, + ) + shear_map_a = shear_map_a_e1 + 1j * shear_map_a_e2 + del shear_map_a_e1, shear_map_a_e2 + + shear_map_b_e1, shear_map_b_e2 = self.get_shear_map( + params, + nside, + cat_gal_b, + unique_pix=unique_pix_b, + idx=idx_b, + idx_rep=idx_rep_b, + n_gal_map=n_gal_map_b, + ) + shear_map_b = shear_map_b_e1 + 1j * shear_map_b_e2 + del shear_map_b_e1, shear_map_b_e2 + + # Compute the pseudo-Cl's + ell_eff, cl_shear, wsp = self.get_pseudo_cls_map( + shear_map_a, n_gal_map_a, shear_map_b=shear_map_b, mask_b=n_gal_map_b + ) + + # Remove the noise bias for auto-correlations. + if tomo_bin_a == tomo_bin_b: + # Compute the noise bias using noise_bias_method + cl_noise = self.get_noise_bias_from_gaussian_real( + params, + nside, + cat_gal_a, + unique_pix=unique_pix_a, + idx=idx_a, + idx_rep=idx_rep_a, + n_gal_map=n_gal_map_a, + wsp=wsp, + ) + + # Subtract the noise bias from the pseudo-Cl's + cl_shear = cl_shear - cl_noise + + self.print_cyan("Saving pseudo-Cl's...") + self.save_pseudo_cl(ell_eff, cl_shear, out_path) + + cl_shear = fits.getdata(out_path) + self._pseudo_cls[ver][f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ + "pseudo_cl" + ] = cl_shear + + def calculate_pseudo_cl_catalog(self, ver, out_path, tomo_bin_a, tomo_bin_b): + assert (tomo_bin_a == "all" and tomo_bin_b == "all") or ( + isinstance(tomo_bin_a, (int, np.integer)) + and isinstance(tomo_bin_b, (int, np.integer)) + ), "tomo_bin_a and tomo_bin_b must be either both 'all' or both integers." + + params = get_params_rho_tau(self.cc[ver]) + + # Load data and create shear and noise maps + cat_gal = fits.getdata(self.cc[ver]["shear"]["path"]) + + ell_eff, cl_shear, wsp = self.get_pseudo_cls_catalog( + catalog=cat_gal, params=params, tomo_bin_a=tomo_bin_a, tomo_bin_b=tomo_bin_b + ) + + self.print_cyan("Saving pseudo-Cl's...") + self.save_pseudo_cl(ell_eff, cl_shear, out_path) + + cl_shear = fits.getdata(out_path) + self._pseudo_cls[ver][f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ + "pseudo_cl" + ] = cl_shear + def calculate_pseudo_cl_inka_cov( self, compute_tomography=True, load_all_block=False ): @@ -71,6 +259,7 @@ def calculate_pseudo_cl_inka_cov( out_path_merged = self._output_path_pseudo_cl_cov( ver, "iNKA", tomography=compute_tomography ) + tomo_str = "tomo" if compute_tomography else "non_tomo" if compute_tomography: tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) @@ -87,7 +276,9 @@ def calculate_pseudo_cl_inka_cov( self.print_done( f"Skipping Pseudo-Cl iNKA covariance calculation, {out_path_merged} exists" ) - self._pseudo_cls[ver]["cov_iNKA"] = fits.open(out_path_merged) + self._pseudo_cls[ver][f"cov_iNKA_{tomo_str}"] = fits.open( + out_path_merged + ) if load_all_block: self.print_done("Loading all the iNKA covariance blocks") @@ -381,13 +572,13 @@ def calculate_pseudo_cl_inka_cov( ] # Merge the covariance blocks - self._pseudo_cls[ver]["cov_iNKA"] = self._merge_iNKA_covariance( + self._pseudo_cls[ver][f"cov_iNKA_{tomo_str}"] = self._merge_iNKA_covariance( ver, tomography=compute_tomography ) self.print_done(f"Done Pseudo-Cl covariance calculation for {ver}") self.print_done("Done Pseudo-Cl covariance") - def calculate_pseudo_cl_onecovariance(self): + def calculate_pseudo_cl_onecovariance(self, tomography=False): """ Compute the pseudo-Cl covariance using OneCovariance. """ @@ -407,11 +598,14 @@ def calculate_pseudo_cl_onecovariance(self): if not os.path.exists(template_config): raise ValueError(f"Template config file {template_config} does not exist") - self._pseudo_cls_onecov = {} + if not hasattr(self, "_pseudo_cls_onecov"): + self._pseudo_cls_onecov = {} for ver in self.versions: self.print_magenta(ver) - - out_dir = self._output_path(f"pseudo_cl_cov_onecov_{ver}/") + self._pseudo_cls_onecov.setdefault(ver, {}) + out_dir = self._output_path( + "pseudo_cl/", f"pseudo_cl_cov_onecov_{ver}_tomography_{tomography}" + ) os.makedirs(out_dir, exist_ok=True) if ( @@ -421,7 +615,7 @@ def calculate_pseudo_cl_onecovariance(self): and not self.force_run ): self.print_done(f"Skipping OneCovariance calculation, {out_dir} exists") - self._load_onecovariance_cov(out_dir, ver) + self._load_onecovariance_cov(out_dir, ver, tomography) else: mask_path = self.cc[ver]["mask"] if not os.path.exists(mask_path): @@ -434,7 +628,9 @@ def calculate_pseudo_cl_onecovariance(self): self.cc[ver]["shear"]["redshift_path"] ) - config_path = os.path.join(out_dir, f"config_onecov_{ver}.ini") + config_path = os.path.join( + out_dir, f"config_onecov_{ver}_tomography_{tomography}.ini" + ) self.print_cyan( f"Modifying OneCovariance config file and saving it to {config_path}" @@ -446,6 +642,7 @@ def calculate_pseudo_cl_onecovariance(self): mask_path, redshift_distr_path, ver, + tomography, ) self.print_cyan("Running OneCovariance...") @@ -457,70 +654,11 @@ def calculate_pseudo_cl_onecovariance(self): f"OneCovariance command failed with return code {ret}" ) self.print_cyan("OneCovariance completed successfully.") - self._load_onecovariance_cov(out_dir, ver) + self._load_onecovariance_cov(out_dir, ver, tomography) self.print_done("Done Pseudo-Cl covariance with OneCovariance") - def _modify_onecov_config( - self, template_config, config_path, out_dir, mask_path, redshift_distr_path, ver - ): - """ - Modify OneCovariance configuration file with correct mask, redshift distribution, - and ellipticity dispersion parameters. - - Parameters - ---------- - template_config : str - Path to the template configuration file - config_path : str - Path where the modified configuration will be saved - mask_path : str - Path to the mask file - redshift_distr_path : str - Path to the redshift distribution file - """ - config = configparser.ConfigParser() - # Load the template configuration - config.read(template_config) - - # Update mask path - mask_base = os.path.basename(os.path.abspath(mask_path)) - mask_folder = os.path.dirname(os.path.abspath(mask_path)) - config["survey specs"]["mask_directory"] = mask_folder - config["survey specs"]["mask_file_lensing"] = mask_base - config["survey specs"]["survey_area_lensing_in_deg2"] = str(self.area[ver]) - config["survey specs"]["ellipticity_dispersion"] = str( - self.ellipticity_dispersion[ver] - ) - config["survey specs"]["n_eff_lensing"] = str(self.n_eff_gal[ver]) - - # Update redshift distribution path - redshift_distr_base = os.path.basename(os.path.abspath(redshift_distr_path)) - redshift_distr_folder = os.path.dirname(os.path.abspath(redshift_distr_path)) - config["redshift"]["z_directory"] = redshift_distr_folder - config["redshift"]["zlens_file"] = redshift_distr_base - - # Update output directory - config["output settings"]["directory"] = out_dir - - # Save the modified configuration - with open(config_path, "w") as f: - config.write(f) - - def _load_onecovariance_cov(self, out_dir, ver): - self.print_cyan(f"Loading OneCovariance results from {out_dir}") - cov_one_cov = np.genfromtxt( - os.path.join(out_dir, "covariance_list_3x2pt_pure_Cell.dat") - ) - gaussian_one_cov = cov_from_one_covariance(cov_one_cov, gaussian=True) - all_one_cov = cov_from_one_covariance(cov_one_cov, gaussian=False) - - self._pseudo_cls_onecov[ver] = { - "gaussian_cov": gaussian_one_cov, - "all_cov": all_one_cov, - } - - def calculate_pseudo_cl_g_ng_cov(self, gaussian_part="iNKA"): + def calculate_pseudo_cl_g_ng_cov(self, tomography=False, gaussian_part="iNKA"): assert gaussian_part in ["iNKA", "OneCovariance"], ( "gaussian_part must be 'iNKA' or 'OneCovariance'" ) @@ -528,34 +666,40 @@ def calculate_pseudo_cl_g_ng_cov(self, gaussian_part="iNKA"): f"Gaussian and Non-Gaussian covariance of the Pseudo-Cl's using {gaussian_part} for the Gaussian part" ) - self._pseudo_cls_cov_g_ng = {} + if not hasattr(self, "_pseudo_cls_cov_g_ng"): + self._pseudo_cls_cov_g_ng = {} for ver in self.versions: self.print_magenta(ver) + self._pseudo_cls_cov_g_ng.setdefault(ver, {}) + key_to_use = "tomo" if tomography else "non_tomo" out_file = self._output_path( - f"pseudo_cl_cov_g_ng_{gaussian_part}_{ver}.fits" + "pseudo_cl", + f"pseudo_cl_cov_g_ng_{gaussian_part}_{ver}_tomography_{tomography}.fits", ) if os.path.exists(out_file) and not self.force_run: self.print_done( f"Skipping Gaussian and Non-Gaussian covariance calculation, {out_file} exists" ) cov_hdu = fits.open(out_file) - self._pseudo_cls_cov_g_ng[ver] = cov_hdu + self._pseudo_cls_cov_g_ng[ver][key_to_use] = cov_hdu continue if gaussian_part == "iNKA": - gaussian_cov = self.pseudo_cls[ver]["cov"]["COVAR_EE_EE"].data + gaussian_cov = self.pseudo_cls[ver][f"cov_iNKA_{key_to_use}"][ + "COVAR_EE_EE" + ].data non_gaussian_cov = ( - self.pseudo_cls_onecov[ver]["all_cov"] - - self.pseudo_cls_onecov[ver]["gaussian_cov"] + self.pseudo_cls_onecov[ver][key_to_use]["all_cov"] + - self.pseudo_cls_onecov[ver][key_to_use]["gaussian_cov"] ) full_cov = gaussian_cov + non_gaussian_cov elif gaussian_part == "OneCovariance": - gaussian_cov = self.pseudo_cls_onecov[ver]["gaussian_cov"] + gaussian_cov = self.pseudo_cls_onecov[ver][key_to_use]["gaussian_cov"] non_gaussian_cov = ( - self.pseudo_cls_onecov[ver]["all_cov"] - - self.pseudo_cls_onecov[ver]["gaussian_cov"] + self.pseudo_cls_onecov[ver][key_to_use]["all_cov"] + - self.pseudo_cls_onecov[ver][key_to_use]["gaussian_cov"] ) - full_cov = self.pseudo_cls_onecov[ver]["all_cov"] + full_cov = self.pseudo_cls_onecov[ver][key_to_use]["all_cov"] else: raise ValueError(f"Unknown gaussian_part: {gaussian_part}") self.print_cyan("Saving Gaussian and Non-Gaussian covariance...") @@ -564,198 +708,11 @@ def calculate_pseudo_cl_g_ng_cov(self, gaussian_part="iNKA"): hdu.append(fits.ImageHDU(non_gaussian_cov, name="COVAR_NON_GAUSSIAN")) hdu.append(fits.ImageHDU(full_cov, name="COVAR_FULL")) hdu.writeto(out_file, overwrite=True) - self._pseudo_cls_cov_g_ng[ver] = hdu + self._pseudo_cls_cov_g_ng[ver][key_to_use] = hdu self.print_done( f"Done Gaussian and Non-Gaussian covariance of the Pseudo-Cl's using {gaussian_part} for the Gaussian part" ) - def calculate_pseudo_cl(self, compute_tomography=True): - """ - Compute the pseudo-Cl of a `CosmologyValidation` inputs with tomography. - """ - out_dir = self._output_path("pseudo_cl") - os.makedirs(out_dir, exist_ok=True) - - if compute_tomography: - self.print_start("Computing tomographic pseudo-Cl's") - else: - self.print_start("Computing non-tomographic pseudo-Cl's") - - self._pseudo_cls = getattr(self, "_pseudo_cls", {}) - - for ver in self.versions: - self.print_magenta(ver) - - if ver not in self.pseudo_cls.keys(): - self._pseudo_cls[ver] = {} - - if compute_tomography: - tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) - - if tomo_bin_ids is None or tomo_bin_pairs is None: - raise ValueError( - f"Version {ver} does not have tomography information." - ) - - else: - tomo_bin_pairs = [("all", "all")] - - # Loop on the different tomographic bin pairs - for bin_key1, bin_key2 in tomo_bin_pairs: - self.print_cyan(f"Tomo Bin Pair: ({bin_key1}, {bin_key2})") - - if ( - f"tomo_bin_{bin_key1}_tomo_bin_{bin_key2}" - not in self._pseudo_cls[ver].keys() - ): - self._pseudo_cls[ver][ - f"tomo_bin_{bin_key1}_tomo_bin_{bin_key2}" - ] = {} - - out_path = self._output_path_pseudo_cl( - ver, tomo_bin_pair=(bin_key1, bin_key2) - ) - if os.path.exists(out_path) and not self.force_run: - self.print_done( - f"Skipping Pseudo-Cl's calculation, {out_path} exists" - ) - cl_shear = fits.getdata(out_path) - self._pseudo_cls[ver][f"tomo_bin_{bin_key1}_tomo_bin_{bin_key2}"][ - "pseudo_cl" - ] = cl_shear - continue - - if self.cell_method == "map": - self.calculate_pseudo_cl_map( - ver, self.nside, out_path, bin_key1, bin_key2 - ) - elif self.cell_method == "catalog": - self.calculate_pseudo_cl_catalog(ver, out_path, bin_key1, bin_key2) - else: - raise ValueError(f"Unknown cell method: {self.cell_method}") - - def calculate_pseudo_cl_map(self, ver, nside, out_path, tomo_bin_a, tomo_bin_b): - assert (tomo_bin_a == "all" and tomo_bin_b == "all") or ( - isinstance(tomo_bin_a, (int, np.integer)) - and isinstance(tomo_bin_b, (int, np.integer)) - ), "tomo_bin_a and tomo_bin_b must be either both 'all' or both integers." - - params = get_params_rho_tau(self.cc[ver]) - - self.print_cyan( - f"Computing pseudo-Cl's for tomographic bins {tomo_bin_a} and {tomo_bin_b}..." - ) - - # Load data and create shear and noise maps - cat_gal = fits.getdata(self.cc[ver]["shear"]["path"]) - - # Get the tomographic bin - cat_gal_a = self._get_tomographic_bin(params, cat_gal, tomo_bin_a) - cat_gal_b = self._get_tomographic_bin(params, cat_gal, tomo_bin_b) - - del cat_gal - - self.print_cyan("Creating maps and computing Cl's...") - # Get the pixels and indices for the catalogs - unique_pix_a, idx_a, idx_rep_a = self.get_pixels(params, nside, cat_gal_a) - unique_pix_b, idx_b, idx_rep_b = self.get_pixels(params, nside, cat_gal_b) - - # Create number density maps for each tomographic bin - n_gal_map_a = self.get_n_gal_map( - params, - nside, - cat_gal_a, - unique_pix=unique_pix_a, - idx=idx_a, - idx_rep=idx_rep_a, - ) - n_gal_map_b = self.get_n_gal_map( - params, - nside, - cat_gal_b, - unique_pix=unique_pix_b, - idx=idx_b, - idx_rep=idx_rep_b, - ) - - # Create shear maps for each tomographic bin - shear_map_a_e1, shear_map_a_e2 = self.get_shear_map( - params, - nside, - cat_gal_a, - unique_pix=unique_pix_a, - idx=idx_a, - idx_rep=idx_rep_a, - n_gal_map=n_gal_map_a, - ) - shear_map_a = shear_map_a_e1 + 1j * shear_map_a_e2 - del shear_map_a_e1, shear_map_a_e2 - - shear_map_b_e1, shear_map_b_e2 = self.get_shear_map( - params, - nside, - cat_gal_b, - unique_pix=unique_pix_b, - idx=idx_b, - idx_rep=idx_rep_b, - n_gal_map=n_gal_map_b, - ) - shear_map_b = shear_map_b_e1 + 1j * shear_map_b_e2 - del shear_map_b_e1, shear_map_b_e2 - - # Compute the pseudo-Cl's - ell_eff, cl_shear, wsp = self.get_pseudo_cls_map( - shear_map_a, n_gal_map_a, shear_map_b=shear_map_b, mask_b=n_gal_map_b - ) - - # Remove the noise bias for auto-correlations. - if tomo_bin_a == tomo_bin_b: - # Compute the noise bias using noise_bias_method - cl_noise = self.get_noise_bias_from_gaussian_real( - params, - nside, - cat_gal_a, - unique_pix=unique_pix_a, - idx=idx_a, - idx_rep=idx_rep_a, - n_gal_map=n_gal_map_a, - wsp=wsp, - ) - - # Subtract the noise bias from the pseudo-Cl's - cl_shear = cl_shear - cl_noise - - self.print_cyan("Saving pseudo-Cl's...") - self.save_pseudo_cl(ell_eff, cl_shear, out_path) - - cl_shear = fits.getdata(out_path) - self._pseudo_cls[ver][f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ - "pseudo_cl" - ] = cl_shear - - def calculate_pseudo_cl_catalog(self, ver, out_path, tomo_bin_a, tomo_bin_b): - assert (tomo_bin_a == "all" and tomo_bin_b == "all") or ( - isinstance(tomo_bin_a, (int, np.integer)) - and isinstance(tomo_bin_b, (int, np.integer)) - ), "tomo_bin_a and tomo_bin_b must be either both 'all' or both integers." - - params = get_params_rho_tau(self.cc[ver]) - - # Load data and create shear and noise maps - cat_gal = fits.getdata(self.cc[ver]["shear"]["path"]) - - ell_eff, cl_shear, wsp = self.get_pseudo_cls_catalog( - catalog=cat_gal, params=params, tomo_bin_a=tomo_bin_a, tomo_bin_b=tomo_bin_b - ) - - self.print_cyan("Saving pseudo-Cl's...") - self.save_pseudo_cl(ell_eff, cl_shear, out_path) - - cl_shear = fits.getdata(out_path) - self._pseudo_cls[ver][f"tomo_bin_{tomo_bin_a}_tomo_bin_{tomo_bin_b}"][ - "pseudo_cl" - ] = cl_shear - # ---------------- Utility functions for pseudo-Cl calculations ---------------- # def get_namaster_bin(self, lmin, lmax, b_lmax): """Build NaMaster binning object (thin wrapper, state -> primitive).""" @@ -970,6 +927,156 @@ def get_fiducial_cl(self, ver, is_tomography): return fiducial_cl + def _modify_onecov_config( + self, + template_config, + config_path, + out_dir, + mask_path, + redshift_distr_path, + ver, + tomography, + ): + """ + Modify OneCovariance configuration file with correct mask, redshift distribution, + and ellipticity dispersion parameters. + + Parameters + ---------- + template_config : str + Path to the template configuration file + config_path : str + Path where the modified configuration will be saved + mask_path : str + Path to the mask file + redshift_distr_path : str + Path to the redshift distribution file + ver : str + Version identifier for the current analysis + tomography : bool + Whether to compute tomography or not + """ + config = configparser.ConfigParser() + # Load the template configuration + config.read(template_config) + + # Update mask path + mask_base = os.path.basename(os.path.abspath(mask_path)) + mask_folder = os.path.dirname(os.path.abspath(mask_path)) + config["survey specs"]["mask_directory"] = mask_folder + config["survey specs"]["mask_file_lensing"] = mask_base + config["survey specs"]["survey_area_lensing_in_deg2"] = str(self.area[ver]) + + # Update ellipticity dispersion and effective number density + # Account for tomography if needed + if tomography: + tomo_bin_ids, _ = self._get_tomo_bins(ver) + else: + tomo_bin_ids = ["all"] + input_ellipticity_distribution = ", ".join( + str(self.ellipticity_dispersion[ver][f"tomo_bin_{bin_id}"]) + for bin_id in tomo_bin_ids + ) + input_n_eff_gal = ", ".join( + str(self.n_eff_gal[ver][f"tomo_bin_{bin_id}"]) for bin_id in tomo_bin_ids + ) + config["survey specs"]["ellipticity_dispersion"] = ( + input_ellipticity_distribution + ) + config["survey specs"]["n_eff_lensing"] = input_n_eff_gal + + # Handle the case where the redshift distribution file is tomographic + # Converts it to a non-tomographic distribution if needed + if not tomography: + # Load the redshift distribution + z, dndz = self.read_redshift_distribution(ver, is_tomography=True) + # Sum over tomographic bins to get the non-tomographic distribution + dndz_non_tomo = np.sum(dndz, axis=1) + # Save the non-tomographic distribution to a new file + non_tomo_redshift_distr_path = os.path.join( + out_dir, f"redshift_distribution_{ver}_non_tomo.txt" + ) + np.savetxt( + non_tomo_redshift_distr_path, + np.column_stack((z, dndz_non_tomo)), + header="z dN/dz", + ) + redshift_distr_path = non_tomo_redshift_distr_path + + # Update redshift distribution path + redshift_distr_base = os.path.basename(os.path.abspath(redshift_distr_path)) + redshift_distr_folder = os.path.dirname(os.path.abspath(redshift_distr_path)) + config["redshift"]["z_directory"] = redshift_distr_folder + config["redshift"]["zlens_file"] = redshift_distr_base + + self._update_onecov_cosmo_params(config, self.cosmo) + + # Update output directory + config["output settings"]["directory"] = out_dir + + # Save the modified configuration + with open(config_path, "w") as f: + config.write(f) + + def _update_onecov_cosmo_params(self, config, cosmo): + """ + Update the cosmological parameters in the OneCovariance configuration. + + Parameters + ---------- + config : configparser.ConfigParser + The configuration object to update. + cosmo + The cosmology object containing the parameters to set. + """ + # Update cosmo params section + config["cosmo"]["h"] = str(cosmo["H0"] / 100.0) + config["cosmo"]["omega_m"] = str(cosmo["Omega_m"]) + config["cosmo"]["omega_b"] = str(cosmo["Omega_b"]) + config["cosmo"]["omega_de"] = str(cosmo["Omega_l"]) + config["cosmo"]["sigma8"] = str(cosmo.sigma8()) + config["cosmo"]["ns"] = str(cosmo["n_s"]) + config["cosmo"]["w0"] = str(cosmo["w0"]) + config["cosmo"]["wa"] = str(cosmo["wa"]) + config["cosmo"]["neff"] = str(cosmo["Neff"]) + config["cosmo"]["m_nu"] = str( + np.sum( + cosmo["m_nu"] + ) # cosmo["m_nu"] is an array of neutrino masses in eV. OneCovariance expects the sum of neutrino masses in eV. + ) + + # Update powspec evaluation section + cosmo_dict = cosmo.to_dict() + + config["powspec evaluation"]["non_linear_model"] = str( + cosmo_dict.get("extra_parameters", {}) + .get("camb", {}) + .get("halofit_version", "mead2020_feedback") # Default to mead2020 feedback + ) # TODO: check what would be the default if not provided as input to align with this + config["powspec evaluation"]["HMCode_logT_AGN"] = str( + cosmo_dict.get("extra_parameters", {}) + .get("camb", {}) + .get("HMCode_logT_AGN", 7.8) # OneCovariance default + ) + + def _load_onecovariance_cov(self, out_dir, ver, tomography): + self.print_cyan(f"Loading OneCovariance results from {out_dir}") + cov_one_cov = np.genfromtxt( + os.path.join(out_dir, "covariance_list_3x2pt_pure_Cell.dat") + ) + gaussian_one_cov = cov_from_one_covariance(cov_one_cov, gaussian=True) + all_one_cov = cov_from_one_covariance(cov_one_cov, gaussian=False) + + key_to_update = "tomo" if tomography else "non_tomo" + self._pseudo_cls_onecov.setdefault(ver, {}).update( + { + key_to_update: { + "gaussian_cov": gaussian_one_cov, + "all_cov": all_one_cov, + } + } + ) + def _get_tomographic_bin(self, params, cat_gal, tomo_bin): """Extract tomographic bin from a given catalogue""" if tomo_bin == "all": diff --git a/src/sp_validation/statistics.py b/src/sp_validation/statistics.py index b8c6d1d5..b1df76b5 100644 --- a/src/sp_validation/statistics.py +++ b/src/sp_validation/statistics.py @@ -7,6 +7,8 @@ Extracted verbatim from the former basic.py. """ +import itertools + import numpy as np from scipy import stats @@ -119,10 +121,60 @@ def cov_from_one_covariance(cov_one_cov, gaussian=True): Square covariance matrix. """ - n_bins = np.sqrt(cov_one_cov.shape[0]).astype(int) - cov = np.zeros((n_bins, n_bins)) - index_value = 10 if gaussian else 9 - for i in range(n_bins): - for j in range(n_bins): - cov[i, j] = cov_one_cov[i * n_bins + j, index_value] + # Get the ell_bins and tomo_bins for each covariance entry + ell1 = cov_one_cov[:, 1] + ell2 = cov_one_cov[:, 2] + tomoi = cov_one_cov[:, 5].astype(int) + tomoj = cov_one_cov[:, 6].astype(int) + tomok = cov_one_cov[:, 7].astype(int) + tomol = cov_one_cov[:, 8].astype(int) + + # Get the values to save in the covariance + cov_col = 10 if gaussian else 9 + values = cov_one_cov[:, cov_col] + + # Map the ell bins to and index + ell_bins = np.unique(ell1) + n_ell_bins = len(ell_bins) + ell_to_idx = {ell: idx for idx, ell in enumerate(ell_bins)} + + # Map the tomo bin pairs to an index + tomo_bin_ids = np.unique(tomoi) + tomo_bin_pairs = list(itertools.combinations_with_replacement(tomo_bin_ids, 2)) + n_spectra = len(tomo_bin_pairs) + pair_to_idx = {pair: idx for idx, pair in enumerate(tomo_bin_pairs)} + + # Initialize the covariance matrix + cov_size = n_ell_bins * n_spectra + cov = np.zeros((cov_size, cov_size)) + + # Get the tomo bin pair indices + # the ordering of OneCovariance is the same than itertools + a_idx = np.fromiter( + (pair_to_idx[(bin_i, bin_j)] for bin_i, bin_j in zip(tomoi, tomoj)), + dtype=int, + count=len(tomoi), + ) + b_idx = np.fromiter( + (pair_to_idx[(bin_k, bin_l)] for bin_k, bin_l in zip(tomok, tomol)), + dtype=int, + count=len(tomok), + ) + + # Get the ell bin indices + ell1_idx = np.fromiter( + (ell_to_idx[ell] for ell in ell1), dtype=int, count=len(ell1) + ) + ell2_idx = np.fromiter( + (ell_to_idx[ell] for ell in ell2), dtype=int, count=len(ell2) + ) + + # Get the row and col + row = a_idx * n_ell_bins + ell1_idx + col = b_idx * n_ell_bins + ell2_idx + + # Assign the values + cov[row, col] = values + cov[col, row] = values # Symmetrize + return cov diff --git a/src/sp_validation/tests/test_statistics.py b/src/sp_validation/tests/test_statistics.py index 8e56d040..1b52d87b 100644 --- a/src/sp_validation/tests/test_statistics.py +++ b/src/sp_validation/tests/test_statistics.py @@ -150,47 +150,129 @@ def test_chi2_and_pte_diagonal_reduces_to_sum_of_squares(): def test_cov_from_one_covariance_selects_gaussian_column(): - """Pin the reshaped matrix and prove the gaussian column selection. - - WHAT IS PINNED: ``cov_from_one_covariance`` reads a flat OneCovariance - table with one row per ``(i, j)`` pair (row-major, ``k = i*n_bins + j``) - and lays the chosen column into a square matrix. Column 10 carries the - Gaussian-only term (``gaussian=True``); column 9 the Gaussian+non-Gaussian - term (``gaussian=False``). With a deterministic table whose entries encode - their row and column, both reshaped matrices are pinned as literals. - - WHY TEETH: the only difference between the two calls is the column index - (10 vs 9), so the two pinned matrices differ by exactly 1 in every entry, - proving the gaussian flag selects the right column. A companion check - places a unique tag ``10*i + j`` in column 10 and asserts the reshape is - row-major (``cov[i, j]`` lands at row ``i*n_bins + j``), so a transposed - or column-major refactor would change the recovered matrix. - """ - # Two-bin table -> 4 rows; entry (row k, col c) = 10*k + c, so the - # chosen-column values are distinct and self-documenting. Row k carries - # the (i, j) pair with k = i*n_bins + j, so rows 0..3 -> (0,0),(0,1), - # (1,0),(1,1). Column 10 holds values 10, 20, 30, 40; column 9 holds - # 9, 19, 29, 39. - one_cov = np.array([np.arange(11.0) + 10.0 * k for k in range(4)]) + """Pin the reshaped matrix and prove the gaussian column selection.""" + # obs, ell1, ell2, s1, s2, tomoi, tomoj, tomok, tomol, cov, covg, covng, covssc + one_cov = np.array( + [ + [0.0, 10.0, 10.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 99.0, 100.0, 0.0, 0.0], + [0.0, 10.0, 20.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 100.0, 101.0, 0.0, 0.0], + [0.0, 20.0, 20.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 110.0, 111.0, 0.0, 0.0], + ] + ) cov_gauss = cov_from_one_covariance(one_cov, gaussian=True) cov_nongauss = cov_from_one_covariance(one_cov, gaussian=False) - npt.assert_allclose(cov_gauss, [[10.0, 20.0], [30.0, 40.0]], rtol=1e-12) - npt.assert_allclose(cov_nongauss, [[9.0, 19.0], [29.0, 39.0]], rtol=1e-12) + npt.assert_allclose(cov_gauss, [[100.0, 101.0], [101.0, 111.0]], rtol=1e-12) + npt.assert_allclose(cov_nongauss, [[99.0, 100.0], [100.0, 110.0]], rtol=1e-12) - # TEETH: the gaussian flag shifts the column by one, so every entry of the - # gaussian matrix exceeds its non-gaussian counterpart by exactly 1. + # TEETH: the gaussian flag shifts the column by one, so every entry of + # the gaussian matrix exceeds its non-gaussian counterpart by exactly 1. npt.assert_allclose(cov_gauss - cov_nongauss, 1.0, rtol=1e-12) - # TEETH: the (i, j) -> row k = i*n_bins + j layout is row-major. Tag - # column 10 with 10*i + j and check it lands at cov[i, j], not cov[j, i]. - tagged = np.zeros((4, 11)) - for i in range(2): - for j in range(2): - tagged[i * 2 + j, 10] = 10.0 * i + j + +def test_cov_from_one_covariance_orders_tomo_blocks(): + """Pin the tomo-block ordering against ``combinations_with_replacement``.""" + # obs, ell1, ell2, s1, s2, tomoi, tomoj, tomok, tomol, cov, covg, covng, covssc + one_cov = np.array( + [ + [ + 0.0, + 10.0, + 10.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + -1.0, + 0.0, + 0.0, + 0.0, + ], # (a,b)=(0,0) + [ + 0.0, + 10.0, + 10.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 2.0, + 0.0, + 1.0, + 0.0, + 0.0, + ], # (0,1) + [ + 0.0, + 10.0, + 10.0, + 1.0, + 1.0, + 1.0, + 1.0, + 2.0, + 2.0, + 1.0, + 2.0, + 0.0, + 0.0, + ], # (0,2) + [ + 0.0, + 10.0, + 10.0, + 1.0, + 1.0, + 1.0, + 2.0, + 1.0, + 2.0, + 10.0, + 11.0, + 0.0, + 0.0, + ], # (1,1) + [ + 0.0, + 10.0, + 10.0, + 1.0, + 1.0, + 1.0, + 2.0, + 2.0, + 2.0, + 11.0, + 12.0, + 0.0, + 0.0, + ], # (1,2) + [ + 0.0, + 10.0, + 10.0, + 1.0, + 1.0, + 2.0, + 2.0, + 2.0, + 2.0, + 21.0, + 22.0, + 0.0, + 0.0, + ], # (2,2) + ] + ) + + cov_gauss = cov_from_one_covariance(one_cov, gaussian=True) + npt.assert_allclose( - cov_from_one_covariance(tagged, gaussian=True), - [[0.0, 1.0], [10.0, 11.0]], + cov_gauss, + [[0.0, 1.0, 2.0], [1.0, 11.0, 12.0], [2.0, 12.0, 22.0]], rtol=1e-12, )