diff --git a/.gitignore b/.gitignore index 8f583d18..553a1d09 100644 --- a/.gitignore +++ b/.gitignore @@ -174,6 +174,7 @@ notebooks/demo_add_bands_to_empty.ipynb notebooks/demo_comprehensive_to_minimal_cat.ipynb notebooks/demo_create_footprint_mask.ipynb notebooks/demo_check_footprint.ipynb +notebooks/demo_comprehensive_to_minimal_cat.ipynb notebooks/demo_calibrate_minimal_cat.ipynb notebooks/leakage_minimal.ipynb notebooks/demo_add_bands.ipynb diff --git a/config/calibration/mask_v1.X.10.yaml b/config/calibration/mask_v1.X.10.yaml new file mode 100644 index 00000000..4729bfdc --- /dev/null +++ b/config/calibration/mask_v1.X.10.yaml @@ -0,0 +1,110 @@ +# Config file for masking and calibration. +# Standard cuts without halo masks, type v1.X.10 + +# General parameters (can also given on command line) +params: + input_path: unions_shapepipe_comprehensive_struc_2024_v1.X.c.hdf5 + cmatrices: False + sky_regions: False + verbose: True + +# Masks +## Using columns in 'dat' group (ShapePipe flags) +dat: + # SExtractor flags + - col_name: FLAGS + label: SE FLAGS + kind: smaller_equal + value: 2 + + # Duplicate objects + - col_name: overlap + label: tile overlap + kind: equal + value: True + + # ShapePipe flags + - col_name: IMAFLAGS_ISO + label: SP mask + kind: equal + value: 0 + + # Number of epochs + - col_name: N_EPOCH + label: r"$n_{\rm epoch}$" + kind: greater_equal + value: 2 + + # Magnitude range + - col_name: mag + label: mag range + kind: range + value: [15, 30] + + # ngmix flags + - col_name: NGMIX_MOM_FAIL + label: "ngmix moments failure" + kind: equal + value: 0 + + # invalid PSF ellipticities + - col_name: NGMIX_ELL_PSFo_NOSHEAR_0 + label: "bad PSF ellipticity comp 1" + kind: not_equal + value: -10 + - col_name: NGMIX_ELL_PSFo_NOSHEAR_1 + label: "bad PSF ellipticity comp 2" + kind: not_equal + value: -10 + +## Using columns in 'dat_ext' group (post-processing flags) +dat_ext: + + # Stars + - col_name: 4_Stars + label: "Stars" + kind: equal + value: False + + # Manual mask + - col_name: 8_Manual + label: "manual mask" + kind: equal + value: False + + # r-band footprint + - col_name: 64_r + label: "r-band imaging" + kind: equal + value: False + + # Maximask + - col_name: 1024_Maximask + label: "maximask" + kind: equal + value: False + + # Rough pointing coverage + - col_name: npoint3 + label: r"$n_{\rm pointing}$" + kind: greater_equal + value: 3 + +# Metacal parameters +metacal: + # Ellipticity dispersion + sigma_eps_prior: 0.34 + + # Signal-to-noise range + gal_snr_min: 10 + gal_snr_max: 500 + + # Relative-size (hlr / hlr_psf) range + gal_rel_size_min: 0.707 + gal_rel_size_max: 3 + + # Correct relative size for ellipticity? + gal_size_corr_ell: False + + # Weight for global response matrix, None for unweighted mean + global_R_weight: w diff --git a/notebooks/cosmo_val/cat_config.yaml b/notebooks/cosmo_val/cat_config.yaml index 6c80062a..44e2cea1 100644 --- a/notebooks/cosmo_val/cat_config.yaml +++ b/notebooks/cosmo_val/cat_config.yaml @@ -88,48 +88,6 @@ SP_axel_v0.0: e1_col: E1_STAR_HSM e2_col: E2_STAR_HSM path: star_cat.fits -SP_test: - subdir: /n17data/mkilbing/astro/data/CFIS/v1.0/ShapePipe - pipeline: SP - colour: black - getdist_colour: 0.0, 0.5, 1.0 - ls: dashed - marker: d - cov_th: - A: 2420.2651014497287 - n_e: 7.040818382014773 - n_psf: 0.752316232272063 - sigma_e: 0.30961528707207325 - 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.4.1.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 - covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt - path: v1.4.x/v1.4.5/unions_shapepipe_cut_struc_2024_v1.4.5.fits - redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt - w_col: w_iv - e1_col: e1 - e1_PSF_col: e1_PSF - e2_col: e2 - e2_PSF_col: e2_PSF - star: - ra_col: RA - dec_col: Dec - e1_col: e1 - e2_col: e2 - path: unions_shapepipe_star_2024_v1.4.1.fits SP_v0.1.1: subdir: /n17data/mkilbing/astro/data/CFIS/v0.0 pipeline: SP @@ -200,104 +158,19 @@ SP_v1.3: e1_col: e1 e2_col: e2 path: unions_shapepipe_star_2022_v1.0.3.fits -SP_v1.4-P3: - subdir: /n17data/mkilbing/astro/data/CFIS/v1.0/ShapePipe/P3 +SP_v1.3.6: + subdir: /n17data/UNIONS/WL/v1.3.x pipeline: SP - colour: brown + colour: violet getdist_colour: 0.0, 0.5, 1.0 ls: dashed - marker: d + marker: h cov_th: - A: 337 - n_e: 7.45 - n_psf: 0.34 - sigma_e: 0.32 - path_rho: rho_stats_SP_v1.4-P3.fits - path_tau: tau_stats_SP_v1.4-P3.fits - path_xi_plus: xi_plus_SP_v1.4-P3.fits - path_xi_minus: xi_minus_SP_v1.4-P3.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: 2 - path: unions_shapepipe_psf_2022_P3_v1.4.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 - covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt - path: unions_shapepipe_extended_2022_P3_v1.4.fits - w_col: w - e1_col: e1 - e1_PSF_col: e1_PSF - e2_col: e2 - e2_PSF_col: e2_PSF - star: - ra_col: RA - dec_col: Dec - e1_col: e1 - e2_col: e2 - path: unions_shapepipe_star_2022_P3_v1.4.fits -SP_v1.4.1_noleakage: - subdir: /n17data/mkilbing/astro/data/CFIS/v1.0/ShapePipe - pipeline: SP - colour: black - getdist_colour: 0.0, 0.5, 1.0 - ls: dashed - marker: d - cov_th: - A: 3191.1827474082734 - n_e: 7.0582373812033286 - n_psf: 0.5705736293858423 - sigma_e: 0.317469332009887 - 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.4.1.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 - covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt - path: unions_shapepipe_2024_v1.4.1.fits - w_col: w_iv - e1_col: e1 - e1_PSF_col: e1_PSF - e2_col: e2 - e2_PSF_col: e2_PSF - star: - ra_col: RA - dec_col: Dec - e1_col: e1 - e2_col: e2 - path: unions_shapepipe_star_2024_v1.4.1.fits -SP_v1.4.2: - subdir: /n17data/mkilbing/astro/data/UNIONS/v1.x/ShapePipe/v1.4.x/v1.4.2 - pipeline: SP - colour: black - getdist_colour: 0.0, 0.5, 1.0 - ls: dashed - marker: d - cov_th: - A: 1826 - n_e: 7.75 - n_psf: 0.35 - sigma_e: 0.3133 + 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_v1.4.6_nside_8192.fits psf: PSF_flag: FLAG_PSF_HSM PSF_size: SIGMA_PSF_HSM @@ -305,7 +178,7 @@ SP_v1.4.2: star_flag: FLAG_STAR_HSM star_size: SIGMA_STAR_HSM hdu: 1 - path: unions_shapepipe_psf_2024_v1.4.a.fits + path: unions_shapepipe_psf_2022_v1.3.a.fits ra_col: RA dec_col: Dec e1_PSF_col: E1_PSF_HSM @@ -315,18 +188,21 @@ SP_v1.4.2: shear: R: 1.0 covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt - path: unions_shapepipe_cut_struc_2024_v1.4.2.fits - w_col: w_iv + path: v1.3.6/unions_shapepipe_cut_struc_2022_v1.3.6.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des e1_col: e1 + e1_col_corrected: e1_leak_corrected e1_PSF_col: e1_PSF e2_col: e2 + e2_col_corrected: e2_leak_corrected e2_PSF_col: e2_PSF star: ra_col: RA dec_col: Dec e1_col: e1 e2_col: e2 - path: unions_shapepipe_star_2024_v1.4.a.fits + path: unions_shapepipe_star_2022_v1.0.3.fits SP_v1.4.5: subdir: /n17data/UNIONS/WL/v1.4.x pipeline: SP @@ -629,6 +505,51 @@ SP_v1.4.6: e1_col: e1 e2_col: e2 path: unions_shapepipe_star_2024_v1.4.a.fits +SP_v1.4.6.1: + subdir: /n17data/UNIONS/WL/v1.4.x + pipeline: SP + colour: red + getdist_colour: 0.0, 0.5, 1.0 + ls: dashed + marker: p + 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_v1.4.6_nside_8192.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.4.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 + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + path: v1.4.6.1/unions_shapepipe_cut_struc_2024_v1.4.6.1.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des + e1_col: e1 + e1_col_corrected: e1_leak_corrected + e1_PSF_col: e1_PSF + e2_col: e2 + e2_col_corrected: e2_leak_corrected + e2_PSF_col: e2_PSF + star: + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + path: unions_shapepipe_star_2024_v1.4.a.fits patch_number: 100 SP_v1.4.7: subdir: /n17data/UNIONS/WL/v1.4.x @@ -722,7 +643,51 @@ SP_v1.4.8: e1_col: e1 e2_col: e2 path: unions_shapepipe_star_2024_v1.4.a.fits - +SP_v1.4.10.1: + subdir: /n17data/UNIONS/WL/v1.4.x + pipeline: SP + colour: blue + getdist_colour: 0.0, 0.5, 1.0 + ls: dashed + marker: s + 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_v1.4.6_nside_8192.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.4.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 + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + path: v1.4.10.1/unions_shapepipe_cut_struc_2024_v1.4.10.1.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des + e1_col: e1 + e1_col_corrected: e1_leak_corrected + e1_PSF_col: e1_PSF + e2_col: e2 + e2_col_corrected: e2_leak_corrected + e2_PSF_col: e2_PSF + star: + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + path: unions_shapepipe_star_2024_v1.4.a.fits # For additive bias only SP_v1.4.6_uncal: pipeline: SP @@ -827,6 +792,7 @@ SP_v1.4.8_uncal: patch_number: 100 SP_v1.4_LFmask_8k: +>>>>>>> upstream/develop subdir: /n17data/mkilbing/astro/data/CFIS/v1.0/SP_LFmask pipeline: SP colour: black @@ -966,3 +932,49 @@ nz: path: dndz paths: output: ./output +SP_v1.6.6: + subdir: /n17data/UNIONS/WL/v1.6.x + 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_v1.4.6_nside_8192.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 + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + path: v1.6.6/unions_shapepipe_cut_struc_2024_v1.6.6.fits + redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + w_col: w_des + e1_col: e1 + e1_col_corrected: e1_leak_corrected + e1_PSF_col: e1_PSF + e2_col: e2 + e2_col_corrected: e2_leak_corrected + e2_PSF_col: e2_PSF + cols: RA,Dec + 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/notebooks/cosmo_val/run_cosmo_val.py b/notebooks/cosmo_val/run_cosmo_val.py index f08f7c0a..2fae1705 100644 --- a/notebooks/cosmo_val/run_cosmo_val.py +++ b/notebooks/cosmo_val/run_cosmo_val.py @@ -57,6 +57,9 @@ # %% cv.plot_objectwise_leakage() +# %% +cv.plot_objectwise_leakage_aux() + # %% #cv.plot_ellipticity() @@ -64,11 +67,9 @@ cv.plot_weights() # %% -cv.plot_separation() +cv.calculate_additive_bias() # %% -cv.npatch = 1 -cv.treecorr_config["var_method"] = "shot" cv.plot_2pcf() cv.treecorr_config["var_method"] = "jackknife" cv.npatch = 100 diff --git a/notebooks/demo_apply_hsp_masks.py b/notebooks/demo_apply_hsp_masks.py index e9e718d3..1760cd7f 100644 --- a/notebooks/demo_apply_hsp_masks.py +++ b/notebooks/demo_apply_hsp_masks.py @@ -22,6 +22,7 @@ # + import os +import re import numpy as np import tracemalloc import healsparse as hsp @@ -57,16 +58,23 @@ # Set parameters base = "unions_shapepipe_comprehensive" year = 2024 -ver = "v1.5.c" +if ver == "v1.3.c": + year = 2022 +else: + year = 2024 +ver_maj = "v1.5" + +obj._params["input_path"] = f"{base}_{year}_{ver_maj}.c.hdf5" +obj._params["output_path"] = f"{base}_struc_{year}_{ver_maj}.c.hdf5" +obj._params["mask_dir"] = f"{os.environ['HOME']}/{ver_maj}.x/masks" -obj._params["input_path"] = f"{base}_{year}_{ver}.hdf5" -obj._params["output_path"] = f"{base}_struc_{year}_{ver}.hdf5" -obj._params["mask_dir"] = "/n17data/UNIONS/WL/masks" obj._params["nside"] = 131072 obj._params["file_base"] = "mask_r_" obj._params["bits"] = bits -obj._params["aux_mask_files"] = f"{obj._params['mask_dir']}/coverage_v1.5.x.hsp" +obj._params["aux_mask_files"] = f"{obj._params['mask_dir']}/coverage_{ver_maj}.x.hsp" +ver_x = re.sub(r"c$", "x", ver) +obj._params["aux_mask_files"] = f"{obj._params['mask_dir']}/coverage_{ver_x}.hsp" obj._params["aux_mask_labels"] = "npoint3" obj._params["verbose"] = True # - diff --git a/notebooks/extract_info.py b/notebooks/extract_info.py index 94d36db6..0bddceab 100644 --- a/notebooks/extract_info.py +++ b/notebooks/extract_info.py @@ -163,6 +163,19 @@ ) # - +# MKDEBUG: Moved from end of this script + +# ### Write PSF catalogue with multi-epoch shapes from shape measurement methods + +spv_cat.write_PSF_cat( + f'{output_PSF_cat_base}_{shape}.fits', + ra_star, + dec_star, + g_star_psf[0], + g_star_psf[1], +) + + #### Refine: Match to SPREAD_CLASS samples if "SPREAD_CLASS" in dd.dtype.names: spv_cat.match_spread_class(dd, ind_star, m_star, stats_file, len(ra_star), verbose=verbose) @@ -250,6 +263,7 @@ ) # - +do_selection_calibration = False if not do_selection_calibration: if verbose: print("No selection and calibration; exiting here") @@ -1007,13 +1021,4 @@ dec = dd['DEC'][cut_overlap] tile_id = dd['TILE_ID'][cut_overlap] write_galaxy_cat(f'{output_shape_cat_base}.fits', ra, dec, tile_id) - -# ### Write PSF catalogue with multi-epoch shapes from shape measurement methods - -write_PSF_cat( - f'{output_PSF_cat_base}_{shape}.fits', - ra_star, - dec_star, - g_star_psf[0], - g_star_psf[1], -) + diff --git a/pyproject.toml b/pyproject.toml index d0ee1d58..5880f8fa 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ authors = [ ] license = {text = "MIT"} readme = "README.md" -requires-python = ">=3.6,<3.12" +requires-python = ">=3.6,<3.13" classifiers = [ "License :: OSI Approved :: MIT License", "Programming Language :: Python :: 3", @@ -20,8 +20,6 @@ dependencies = [ "camb", "clmm", "colorama", - "cosmo-numba @ git+https://github.com/aguinot/cosmo-numba", - "cs_util>=0.1.5", "emcee", "getdist @ git+https://github.com/benabed/getdist.git@upper_triangle_whisker", "h5py", @@ -39,18 +37,17 @@ dependencies = [ "opencv-python-headless", "pyccl", "pyarrow", - "pymaster", "regions", "reproject", "scipy", "seaborn", "skyproj", - "shear_psf_leakage @ git+https://github.com/CosmoStat/shear_psf_leakage.git@develop", "statsmodels", "treecorr>=5.0", "tqdm", "uncertainties" ] +# "cosmo-numba @ git+https://github.com/aguinot/cosmo-numba", [project.urls] Homepage = "https://github.com/CosmoStat/sp_validation" diff --git a/scripts/merge_psf_cat.py b/scripts/merge_psf_cat.py new file mode 100644 index 00000000..43f680d5 --- /dev/null +++ b/scripts/merge_psf_cat.py @@ -0,0 +1,217 @@ +#!/usr/bin/env python3 +"""MERGE PSF CAT. + +Merge PSF catalogues (psf_catalog_ngmix.fits) from different patches +into a single FITS file. + +:Author: Martin Kilbinger +""" + +import sys + +import numpy as np +from astropy.io import fits + +from cs_util import logging +from cs_util import cat +from cs_util import args as cs_args + + +class MergePsfCat: + """Merge Psf Cat. + + Class to merge PSF catalogues from multiple patches. + + """ + + def __init__(self): + self.params_default() + + def params_default(self): + """Params Default. + + Set default parameter values. + + """ + self._params = { + "patches": "v1", + "sh": "ngmix", + "survey": "unions", + "year": "2024", + "version": "1.4.2", + "pipeline": "shapepipe", + "base_path": ".", + "hdu": 1, + "verbose": False, + } + self._short_options = { + "patches": "-p", + "sh": "-g", + "survey": "-s", + "year": "-y", + "version": "-V", + "base_path": "-b", + "hdu": "-H", + } + self._types = { + "hdu": "int", + } + self._help_strings = { + "patches": ( + "list of patches separated by '+', or shortcut " + "(allowed are 'v1', 'v1.5'), default={}" + ), + "sh": "shape measurement method, default={}", + "survey": "survey name, default={}", + "year": "year of processing, default={}", + "version": "catalogue version, default={}", + "base_path": "base path containing patch directories, default={}", + "hdu": "HDU number to read from input FITS files, default={}", + } + + def set_params_from_command_line(self, args): + """Set Params From Command Line. + + Only use when calling using python from command line. + Does not work from ipython or jupyter. + + """ + options = cs_args.parse_options( + self._params, + self._short_options, + self._types, + self._help_strings, + ) + + self._params.update(options) + + logging.log_command(args) + + def get_patches(self): + """Get Patches. + + Return list of patches according to option parameter value. + + Returns + ------- + list + patches, list of str + + """ + if self._params["patches"] == "v1": + n_patch = 7 + patches = [f"P{x}" for x in np.arange(n_patch) + 1] + elif self._params["patches"] == "v1.5": + n_patch = 8 + patches = [f"P{x}" for x in np.arange(n_patch) + 1] + elif self._params["patches"] == "v1.6": + n_patch = 9 + patches = [f"P{x}" for x in np.arange(n_patch) + 1] + else: + patches = self._params["patches"].split("+") + + return patches + + def merge_catalogues(self, patches): + """Merge Catalogues. + + Merge PSF catalogues from sub-patches into one FITS file. + + Parameters + ---------- + patches : list of str + list of patches/sub-directories + + """ + base_path = self._params["base_path"] + sh = self._params["sh"] + hdu_in = self._params["hdu"] + verbose = self._params["verbose"] + + input_sub_path = f"sp_output/psf_catalog_{sh}.fits" + output_path = ( + f"{self._params['survey']}_{self._params['pipeline']}_star_" + f"{self._params['year']}_v{self._params['version']}.fits" + ) + + dat_all = {} + for idx, patch in enumerate(patches): + + if verbose: + print(f" {patch}") + + input_path = f"{base_path}/{patch}/{input_sub_path}" + try: + dat = fits.getdata(input_path, hdu_in) + except Exception: + print(f"No data found in file {input_path} at HDU #{hdu_in}") + print(f"Trying at HDU #{hdu_in - 1}") + dat = fits.getdata(input_path, hdu_in - 1) + + if idx == 0: + col_names = dat.dtype.names + for name in col_names: + if name != "w": + dat_all[name] = [] + else: + dat_all[name + "_iv"] = [] + dat_all["patch"] = [] + + for name in col_names: + if name != "w": + dat_all[name] = np.append(dat_all[name], dat[name]) + else: + dat_all[name + "_iv"] = np.append( + dat_all[name + "_iv"], dat[name] + ) + + dat_all["patch"] = np.append(dat_all["patch"], [idx + 1] * len(dat)) + + col_names = col_names + ("patch",) + + column_all = [] + for name in col_names: + if name != "patch": + my_format = "D" + else: + my_format = "I" + if name == "w": + name = "w_iv" + column = fits.Column(name=name, array=dat_all[name], format=my_format) + column_all.append(column) + + if verbose: + print(f"Writing file {output_path}") + cat.write_fits_BinTable_file(column_all, output_path) + + def run(self): + """Run. + + Main processing function. + + """ + patches = self.get_patches() + + if self._params["verbose"]: + print("Merging PSF catalogues from patches:", patches) + + self.merge_catalogues(patches) + + +def main(argv=None): + """Main. + + Main program. + + """ + obj = MergePsfCat() + + obj.set_params_from_command_line(argv) + + obj.run() + + return 0 + + +if __name__ == "__main__": + sys.exit(main(sys.argv)) diff --git a/scripts/print_matrix_from_fits.py b/scripts/print_matrix_from_fits.py new file mode 100755 index 00000000..0b284b19 --- /dev/null +++ b/scripts/print_matrix_from_fits.py @@ -0,0 +1,197 @@ +#!/usr/bin/env python +"""Read matrix entries from FITS header and print in tex or pdf format.""" + +from astropy.io import fits +import numpy as np +import subprocess +import tempfile +import os + +from cs_util.args import parse_options + + +def format_scientific_latex(value, prec, exp_thresh=3): + """Format number in scientific notation for LaTeX. + + Parameters + ---------- + value : float + input value + prec : int + number of significant digits + exp_thresh : int, optional + exponent threshold for switching to scientific notation; + use scientific notation if |exponent| >= exp_thresh; + default is 3 + + Returns + ------- + str + LaTeX-formatted string, e.g. "1.23 \\times 10^{-4}" + + """ + if value == 0: + return "0" + + # Get the exponent + exp = int(np.floor(np.log10(np.abs(value)))) + + # Use regular notation if exponent is small + if abs(exp) < exp_thresh: + return f"{value:.{prec}g}" + + # Get the mantissa + mantissa = value / 10**exp + + # Format mantissa with specified precision + mantissa_str = f"{mantissa:.{prec-1}f}" + + return f"{mantissa_str} \\times 10^{{{exp}}}" + + +def get_matrix_from_header(header, prefix): + """Extract 2x2 matrix from FITS header using prefix. + + Looks for keys: prefix11, prefix12, prefix21, prefix22 + """ + matrix = np.zeros((2, 2)) + for i in range(1, 3): + for j in range(1, 3): + key = f"{prefix}{i}{j}" + if key not in header: + raise KeyError(f"Key '{key}' not found in FITS header") + matrix[i-1, j-1] = header[key] + return matrix + + +def matrix_to_tex(matrix, prec, exp_thresh): + """Convert matrix to LaTeX format.""" + rows = [] + for i in range(2): + row_vals = [format_scientific_latex(matrix[i, j], prec, exp_thresh) for j in range(2)] + rows.append(" & ".join(row_vals)) + + tex = r"\begin{pmatrix}" + "\n" + tex += (" " + r" \\" + "\n").join(rows) + "\n" + tex += r"\end{pmatrix}" + return tex + + +def generate_pdf(tex_content, output_path): + """Generate PDF from LaTeX content.""" + full_tex = r"""\documentclass[preview,border=2pt]{article} +\usepackage{amsmath} +\begin{document} +$""" + tex_content + r"""$ +\end{document} +""" + + with tempfile.TemporaryDirectory() as tmpdir: + tex_file = os.path.join(tmpdir, "matrix.tex") + with open(tex_file, "w") as f: + f.write(full_tex) + + # Run pdflatex + result = subprocess.run( + ["pdflatex", "-interaction=nonstopmode", "-output-directory", tmpdir, tex_file], + capture_output=True, + text=True + ) + + if result.returncode != 0: + raise RuntimeError(f"pdflatex failed: {result.stderr}") + + pdf_file = os.path.join(tmpdir, "matrix.pdf") + if output_path: + import shutil + shutil.copy(pdf_file, output_path) + print(f"PDF written to {output_path}") + else: + print(f"PDF generated at {pdf_file}") + + +def main(): + # Default parameter values + p_def = { + "input": None, + "prefix": "R_S", + "format": "tex", + "prec": 3, + "exp_thresh": 3, + "output": None, + "ext": 0, + } + + # Short options + short_options = { + "input": "-i", + "prefix": "-r", + "prec": "-n", + "exp_thresh": "-t", + "output": "-o", + "ext": "-e", + } + + # Types + types = { + "input": "string", + "prefix": "string", + "format": "string", + "prec": "int", + "exp_thresh": "int", + "output": "string", + "ext": "int", + } + + # Help strings + help_strings = { + "input": "Input FITS file", + "prefix": "Prefix for matrix entries (default: {})", + "format": "Output format: tex or pdf (default: {})", + "prec": "Number of significant digits (default: {})", + "exp_thresh": "Exponent threshold for scientific notation (default: {})", + "output": "Output file path (extension added based on format)", + "ext": "FITS extension to read header from (default: {})", + } + + options = parse_options(p_def, short_options, types, help_strings) + + if options["input"] is None: + raise ValueError("Input FITS file (-i) is required") + + if options["format"] not in ["tex", "pdf"]: + raise ValueError("Format must be 'tex' or 'pdf'") + + # Read FITS header + with fits.open(options["input"]) as hdu: + header = hdu[options["ext"]].header + + # Extract matrix + matrix = get_matrix_from_header(header, options["prefix"]) + + # Generate output + tex_content = matrix_to_tex(matrix, options["prec"], options["exp_thresh"]) + + # Determine output path with correct extension + output_path = None + if options["output"]: + output_path = f"{options['output']}.{options['format']}" + + if options["format"] == "tex": + full_tex = r"""\documentclass[preview,border=2pt]{article} +\usepackage{amsmath} +\begin{document} +$""" + tex_content + r"""$ +\end{document} +""" + print(full_tex) + if output_path: + with open(output_path, "w") as f: + f.write(full_tex) + print(f"Written to {output_path}") + elif options["format"] == "pdf": + generate_pdf(tex_content, output_path) + + +if __name__ == "__main__": + main() diff --git a/src/sp_validation/basic.py b/src/sp_validation/basic.py index fc96230a..9767ad79 100644 --- a/src/sp_validation/basic.py +++ b/src/sp_validation/basic.py @@ -370,10 +370,19 @@ def _masking_gal(self): else: snr_flux = data['flux'] / data['flux_err'] + if name == 'ns': + # This is a FHP hack, the ns PSF measured in shapepipe is not correct, + # fortunately it is the same dilated PSF as the other branches, + # thus we can simply use p1 + print("FHP using p1 PSF for ns in cuts") + Tpsf = self.p1['Tpsf'] + else: + Tpsf = data['Tpsf'] + mask_tmp = ( (data['flag'] == 0) - & (Tr_tmp / data['Tpsf'] > self._rel_size_min) - & (Tr_tmp / data['Tpsf'] < self._rel_size_max) + & (Tr_tmp / Tpsf > self._rel_size_min) + & (Tr_tmp / Tpsf < self._rel_size_max) & (snr_flux > self._snr_min) & (snr_flux < self._snr_max) ) diff --git a/src/sp_validation/cosmo_val.py b/src/sp_validation/cosmo_val.py index db619c41..e9eda22d 100644 --- a/src/sp_validation/cosmo_val.py +++ b/src/sp_validation/cosmo_val.py @@ -1396,7 +1396,6 @@ def calculate_objectwise_leakage_aux(self): self.print_magenta(ver) results_obj = self.results_objectwise[ver] - print("MKDEBUG query done") results_obj.check_params() results_obj.update_params() results_obj.prepare_output() @@ -1514,26 +1513,6 @@ def plot_weights(self, nbins=200): cs_plots.show() self.print_done("Weight histograms saved to " + out_path) - def plot_separation(self, nbins=200): - self.print_start("Separation histograms") - if "SP_matched_MP_v1.0" in self.versions: - fig, axs = plt.subplots(1, 1, figsize=(10, 7)) - with self.results["SP_matched_MP_v1.0"].temporarily_read_data(): - sep = self.results["SP_matched_MP_v1.0"].dat_shear["Separation"] - axs.hist( - sep, - bins=nbins, - density=False, - histtype="step", - label="SP_matched_MP_v1.0", - color=self.cc["SP_matched_MP_v1.0"]["colour"], - ) - print("Max separation: %s arcsec" % max(sep)) - axs.set_xlabel(r"Separation $\theta$ [arcsec]") - axs.legend() - else: - self.print_done("SP_matched_MP_v1.0 not in versions") - def calculate_additive_bias(self): self.print_start("Calculating additive bias:") self._c1 = {}