From 1fdd6443528435a713060b2c7bd79f340e3a5985 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 12 Nov 2025 11:55:28 +0100 Subject: [PATCH 1/7] added script to plot r-band magnitude histograms for different selections; cleaned up plotting functions --- .../cosmo_val/catalog_paper_plot/hist_mag.py | 384 ++++++++++++++++++ src/sp_validation/__init__.py | 22 +- src/sp_validation/basic.py | 2 - src/sp_validation/plots.py | 83 ++++ src/sp_validation/run_joint_cat.py | 102 +---- 5 files changed, 490 insertions(+), 103 deletions(-) create mode 100644 notebooks/cosmo_val/catalog_paper_plot/hist_mag.py diff --git a/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py b/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py new file mode 100644 index 0000000..3ee2c88 --- /dev/null +++ b/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py @@ -0,0 +1,384 @@ +# %% +# hist_mag.py +# +# Plot magnitude histogram for various cuts and selection criteria + +# %% +import matplotlib +import matplotlib.pylab as plt + +# enable autoreload for interactive sessions +from IPython import get_ipython +ipython = get_ipython() +if ipython is not None: + ipython.run_line_magic("matplotlib", "inline") + ipython.run_line_magic("reload_ext", "autoreload") + ipython.run_line_magic("autoreload", "2") + ipython.run_line_magic("reload_ext", "log_cell_time") + +# %% +import sys +import os +import re +import numpy as np +from astropy.io import fits +from io import StringIO + +from sp_validation import run_joint_cat as sp_joint +from sp_validation import util +from sp_validation.basic import metacal +from sp_validation import calibration +import sp_validation.cat as cat + +# %% +# Initialize calibration class instance +obj = sp_joint.CalibrateCat() + +config = obj.read_config_set_params("config_mask.yaml") + +hist_data_path = "magnitude_histograms_data.npz" + +test_only = False + +# %% +def get_data(obj, test_only=False): + + + # Get data. Set load_into_memory to False for very large files + dat, dat_ext = obj.read_cat(load_into_memory=False) + + if test_only: + n_max = 1_000_000 + print(f"MKDEBUG testing only first {n_max} objects") + dat = dat[:n_max] + dat_ext = dat_ext[:n_max] + + return dat, dat_ext + + +def read_hist_data(hist_data_path): + """ + Read histogram data from npz file. + + Parameters + ---------- + hist_data_path : str + Path to the npz file containing histogram data + + Returns + ------- + hist_data : dict + Dictionary with keys for each selection criterion containing: + - 'counts': histogram counts + - 'bins': bin edges + - 'label': label for the histogram + """ + loaded = np.load(hist_data_path, allow_pickle=True) + hist_data = {} + + for key in loaded.files: + data = loaded[key] + hist_data[key] = { + 'counts': data[0], + 'bins': data[1], + 'label': str(data[2]) + } + + return hist_data + + + +def get_mask(masks, col_name): + + # Get mask fomr masks with col_name = col_name + for mask in masks: + if mask._col_name == col_name: + return mask + + +def compute_hist(masks, col_name, mask_cumul, mag, bins): + """ + Compute histogram for given mask and magnitude data. + + Parameters + ---------- + masks : list + List of mask objects + col_name : str + Column name to identify the mask + mask_cumul : array or None + Cumulative mask array + mag : array + Magnitude data + bins : array + Bin edges for histogram + + Returns + ------- + counts : array + Histogram counts + bins : array + Bin edges + label : str + Label for the histogram + n_valid : int + Number of valid data points after masking + mask_cumul : array + Updated cumulative mask array + """ + this_mask = get_mask(masks, col_name) + + # First time: + if mask_cumul is None: + print("Init mask_cumul") + mask_cumul = this_mask._mask + else: + mask_cumul &= this_mask._mask + + # Data values + my_mag = mag[mask_cumul] + + # Data count + n_valid = np.sum(mask_cumul) + + # Label + string_buffer = StringIO() + this_mask.print_condition(string_buffer) + label = string_buffer.getvalue().strip() + if label == "": + label = col_name + + counts, bin_edges = np.histogram(my_mag, bins=bins) + + return counts, bin_edges, label, n_valid, mask_cumul + + +def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): + """ + Plot histogram from counts and bins using bar plot. + + Parameters + ---------- + counts : array + Histogram counts + bins : array + Bin edges + label : str + Label for the histogram + alpha : float + Transparency for plot + ax : matplotlib axis + Axis to plot on + color : str or None + Color for the histogram + + """ + bin_centers = 0.5 * (bins[:-1] + bins[1:]) + ax.bar( + bin_centers, + counts, + width=np.diff(bins), + alpha=alpha, + label=label, + align='center', + color=color + ) + +# %% + +if os.path.exists(hist_data_path): + print(f"Histogram data file {hist_data_path} found.") + print("Reading and plotting.") + + dat = dat_ext = None + hist_data = read_hist_data(hist_data_path) + +else: + print(f"Histogram data file {hist_data_path} not found.") + print("Reading UNIONS cat and computing.") + + dat, dat_ext = get_data(obj, test_only=test_only) + hist_data = None + + +# ## Masking +# %% +# Get all masks, with or without dat, dat_ext +masks, labels = sp_joint.get_masks_from_config( + config, + dat, + dat_ext, + verbose=obj._params["verbose"], +) + +# %% +# List of basic masks to apply to all cases +masks_labels_basic = [ + "FLAGS", + "overlap", + "IMAFLAGS_ISO", + "mag", + "NGMIX_MOM_FAIL", + "NGMIX_ELL_PSFo_NOSHEAR_0", + "NGMIX_ELL_PSFo_NOSHEAR_1", + "4_Stars", + "8_Manual", + "64_r", + "1024_Maximask", +] + +# %% +masks_basic = [] +for mask in masks: + if mask._col_name in masks_labels_basic: + masks_basic.append(mask) + +if dat is not None: + print("Creating combined basic mask") + masks_basic_combined = sp_joint.Mask.from_list( + masks_basic, + label="combined", + verbose=obj._params["verbose"], + ) +else: + # Create dummy combined mask (for plot label) + print("Creating dummy combined basic mask") + masks_basic_combined = sp_joint.Mask( + "combined", + "combined", + kind="none", + ) + +masks.append(masks_basic_combined) + +# %% +# Metacal mask (cuts) +mask_tmp = sp_joint.Mask( + "metacal", + "metacal", + kind="none", +) + +# %% +if dat is not None: + cm = config["metacal"] + gal_metacal = metacal( + dat, + masks_basic_combined._mask, + snr_min=cm["gal_snr_min"], + snr_max=cm["gal_snr_max"], + rel_size_min=cm["gal_rel_size_min"], + rel_size_max=cm["gal_rel_size_max"], + size_corr_ell=cm["gal_size_corr_ell"], + sigma_eps=cm["sigma_eps_prior"], + global_R_weight=cm["global_R_weight"], + col_2d=False, + verbose=True, + ) + + g_corr_mc, g_uncorr, w, mask_metacal, c, c_err = ( + calibration.get_calibrated_m_c(gal_metacal) + ) + + # Convert index array to boolean mask + mask_metacal_bool = np.zeros(len(dat), dtype=bool) + mask_metacal_bool[mask_metacal] = True + + mask_tmp._mask = mask_metacal_bool + +masks.append(mask_tmp) + + +# %% +# Plot magnitude histograms for various selection criteria + +# Define magnitude bins +mag_bins = np.arange(15, 30, 0.05) +mag_centers = 0.5 * (mag_bins[:-1] + mag_bins[1:]) + + +# %% +# Create figure with multiple subplots +figsize = 10 +alpha = 0.5 + +col_names = ["combined", "N_EPOCH", "npoint3", "metacal"] + +# Define explicit colors for each histogram +colors = ['C0', 'C1', 'C2', 'C3'] # Use matplotlib default color cycle +color_map = dict(zip(col_names, colors)) + + +# %% +# If hist_data not loaded, compute it +if hist_data is None: + hist_data = {} + +if dat is not None: + # Get magnitude column + mag = dat['mag'] + + mask_cumul = None + for col_name in col_names: + counts, bins, label, n_valid, mask_cumul = compute_hist( + masks=masks, + col_name=col_name, + mask_cumul=mask_cumul, + mag=mag, + bins=mag_bins + ) + hist_data[col_name] = { + 'counts': counts, + 'bins': bins, + 'label': label, + 'n_valid': n_valid, + } + +# Plot histogram data +plt.figure() +fig, (ax) = plt.subplots(1, 1, figsize=(figsize, figsize)) + +for col_name in col_names: + if col_name in hist_data: + data = hist_data[col_name] + plot_hist( + data['counts'], + data['bins'], + data['label'], + alpha=alpha, + ax=ax, + color=color_map[col_name] + ) + #print(f"{col_name}: n_valid = {data['n_valid']}") + +ax.set_xlabel('$r$') +ax.set_ylabel('Number') +ax.set_xlim(17.5, 26.5) +ax.legend() + +plt.tight_layout() +plt.savefig('magnitude_histograms.png', dpi=150, bbox_inches='tight') + +# Save histogram data to file (only if we computed it) +if dat is not None: + np.savez( + hist_data_path, + **{ + key: np.array( + [ + val['counts'], + val['bins'], + val['label'], + val["n_valid"], + ], + dtype=object + ) + for key, val in hist_data.items() + } + ) + print(f"Histogram data saved to {hist_data_path}") + +# %% +if dat is not None: + obj.close_hd5() +# %% diff --git a/src/sp_validation/__init__.py b/src/sp_validation/__init__.py index b31411d..62e2307 100644 --- a/src/sp_validation/__init__.py +++ b/src/sp_validation/__init__.py @@ -16,14 +16,14 @@ ] # Explicit imports to avoid circular issues -from . import util -from . import io -from . import basic -from . import galaxy -from . import cosmology -from . import calibration -from . import cat -from . import plot_style -from . import plots -from . import run_joint_cat -from . import survey \ No newline at end of file +#from . import util +#from . import io +#from . import basic +#from . import galaxy +#from . import cosmology +#from . import calibration +#from . import cat +#from . import plot_style +#from . import plots +#from . import run_joint_cat +#from . import survey diff --git a/src/sp_validation/basic.py b/src/sp_validation/basic.py index 9c0972e..3720004 100644 --- a/src/sp_validation/basic.py +++ b/src/sp_validation/basic.py @@ -15,8 +15,6 @@ from scipy.spatial import cKDTree from scipy.special import gamma -import matplotlib.pyplot as plt - from tqdm import tqdm import operator as op import itertools as itools diff --git a/src/sp_validation/plots.py b/src/sp_validation/plots.py index d9c5d0e..d565048 100644 --- a/src/sp_validation/plots.py +++ b/src/sp_validation/plots.py @@ -669,3 +669,86 @@ def hsp_map_logical_or(maps, verbose=False): ) return map_comb + + +def plot_area_mask(ra, dec, zoom, mask=None): + """Plot Area Mask. + + Create sky plot of objects. + + Parameters + ---------- + ra : list + R.A. coordinates + dec : list + Dec. coordinates + zoom : TBD + mask: TBD, optional + + """ + if mask is None: + mask == np.ones_like(ra) + + fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(30,15)) + axes[0].hexbin(ra[mask], dec[mask], gridsize=100) + axes[1].hexbin(ra[mask & zoom], dec[mask & zoom], gridsize=200) + for idx in (0, 1): + axes[idx].set_xlabel("R.A. [deg]") + axes[idx].set_ylabel("Dec [deg]") + + +def sky_plots(dat, masks, labels, zoom_ra, zoom_dec): + """Sky Plots. + + Plot sky regions with different masks. + + Parameters + ---------- + masks : list + masks to be applied + labels : dict + labels for masks + zoom_ra : list + min and max R.A. for zoom-in plot + zoom_dec : list + min and max Dec. for zoom-in plot + + """ + ra = dat["RA"][:] + dec = dat["Dec"][:] + + zoom_ra = (room_ra[0] < dat["RA"]) & (dat["RA"] < zoom_ra[1]) + zoom_dec = (zoom_dec[0] < dat["Dec"]) & (dat["Dec"] < zoom_dec[1]) + zoom = zoom_ra & zoom_dec + + # No mask + plot_area_mask(ra, dec, zoom) + + # SExtractor and SP flags + m_flags = masks[labels["FLAGS"]]._mask & masks[labels["IMAFLAGS_ISO"]]._mask + plot_area_mask(ra, dec, zoom, mask=m_flags) + + # Overlap regions + m_over = masks[labels["overlap"]]._mask & m_flags + plot_area_mask(ra, dec, zoom, mask=m_over) + + # Coverage mask + m_point = masks[labels["npoint3"]]._mask & m_over + plot_area_mask(ra, dec, zoom, mask=m_point) + + # Maximask + m_maxi = masks[labels["1024_Maximask"]]._mask & m_point + plot_area_mask(ra, dec, zoom, mask=m_maxi) + + m_comb = mask_combined._mask + plot_area_mask(ra, dec, zoom, mask=m_comb) + + m_man = m_maxi & masks[labels["8_Manual"]]._mask + plot_area_mask(ra, dec, zoom, mask=m_man) + + m_halos = ( + m_maxi + & masks[labels['1_Faint_star_halos']]._mask + & masks[labels['2_Bright_star_halos']]._mask + ) + plot_area_mask(ra, dec, zoom, mask=m_halos) diff --git a/src/sp_validation/run_joint_cat.py b/src/sp_validation/run_joint_cat.py index 9177eeb..c9d5965 100644 --- a/src/sp_validation/run_joint_cat.py +++ b/src/sp_validation/run_joint_cat.py @@ -1181,89 +1181,6 @@ def run(self): """ -def sky_plots(dat, masks, labels, zoom_ra, zoom_dec): - """Sky Plots. - - Plot sky regions with different masks. - - Parameters - ---------- - masks : list - masks to be applied - labels : dict - labels for masks - zoom_ra : list - min and max R.A. for zoom-in plot - zoom_dec : list - min and max Dec. for zoom-in plot - - """ - ra = dat["RA"][:] - dec = dat["Dec"][:] - - zoom_ra = (room_ra[0] < dat["RA"]) & (dat["RA"] < zoom_ra[1]) - zoom_dec = (zoom_dec[0] < dat["Dec"]) & (dat["Dec"] < zoom_dec[1]) - zoom = zoom_ra & zoom_dec - - # No mask - plot_area_mask(ra, dec, zoom) - - # SExtractor and SP flags - m_flags = masks[labels["FLAGS"]]._mask & masks[labels["IMAFLAGS_ISO"]]._mask - plot_area_mask(ra, dec, zoom, mask=m_flags) - - # Overlap regions - m_over = masks[labels["overlap"]]._mask & m_flags - plot_area_mask(ra, dec, zoom, mask=m_over) - - # Coverage mask - m_point = masks[labels["npoint3"]]._mask & m_over - plot_area_mask(ra, dec, zoom, mask=m_point) - - # Maximask - m_maxi = masks[labels["1024_Maximask"]]._mask & m_point - plot_area_mask(ra, dec, zoom, mask=m_maxi) - - m_comb = mask_combined._mask - plot_area_mask(ra, dec, zoom, mask=m_comb) - - m_man = m_maxi & masks[labels["8_Manual"]]._mask - plot_area_mask(ra, dec, zoom, mask=m_man) - - m_halos = ( - m_maxi - & masks[labels['1_Faint_star_halos']]._mask - & masks[labels['2_Bright_star_halos']]._mask - ) - plot_area_mask(ra, dec, zoom, mask=m_halos) - - -def plot_area_mask(ra, dec, zoom, mask=None): - """Plot Area Mask. - - Create sky plot of objects. - - Parameters - ---------- - ra : list - R.A. coordinates - dec : list - Dec. coordinates - zoom : TBD - mask: TBD, optional - - """ - if mask is None: - mask == np.ones_like(ra) - - fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(30,15)) - axes[0].hexbin(ra[mask], dec[mask], gridsize=100) - axes[1].hexbin(ra[mask & zoom], dec[mask & zoom], gridsize=200) - for idx in (0, 1): - axes[idx].set_xlabel("R.A. [deg]") - axes[idx].set_ylabel("Dec [deg]") - - def confusion_matrix(mask, confidence_level=0.9): n_key = len(mask) @@ -1361,7 +1278,7 @@ class Mask(): """ - def __init__(self, col_name, label, kind="equal", value=0, dat=None, verbose=False): + def __init__(self, col_name, label, kind=None, value=0, dat=None, verbose=False): self._col_name = col_name self._label = label @@ -1439,8 +1356,7 @@ def print_strings(cls, coln, lab, num, fnum, f_out=None): print(msg) if f_out: print(msg, file=f_out) - - + def print_stats(self, num_obj, f_out=None): if self._num_ok is None: self._num_ok = sum(self._mask) @@ -1461,10 +1377,12 @@ def get_sign(self): elif self._kind =="smaller_equal": sign = "<=" return sign - - def print_summary(self, f_out): - print(f"[{self._label}]\t\t\t", file=f_out, end="") - + + def print_condition(self, f_out): + + if self._value is None: + return "" + sign = self.get_sign() if sign is not None: @@ -1472,6 +1390,10 @@ def print_summary(self, f_out): if self._kind == "range": print(f"{self._value[0]} <= {self._col_name} <= {self._value[1]}", file=f_out) + + def print_summary(self, f_out): + print(f"[{self._label}]\t\t\t", file=f_out, end="") + self.print_condition(f_out) def create_descr(self): """Create Descr. From 9a39d94a83dbfe5f9fa5259fe5da3d3d0cace56c Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Fri, 14 Nov 2025 11:29:43 +0100 Subject: [PATCH 2/7] Improved mask labels for plots --- config/calibration/mask_v1.X.6.yaml | 6 +- .../cosmo_val/catalog_paper_plot/hist_mag.py | 283 ++++++++++++++---- src/sp_validation/run_joint_cat.py | 32 +- 3 files changed, 244 insertions(+), 77 deletions(-) diff --git a/config/calibration/mask_v1.X.6.yaml b/config/calibration/mask_v1.X.6.yaml index c3c7b25..9dee194 100644 --- a/config/calibration/mask_v1.X.6.yaml +++ b/config/calibration/mask_v1.X.6.yaml @@ -31,13 +31,13 @@ dat: # Number of epochs - col_name: N_EPOCH - label: r"$n_{\rm epoch}$" + label: $n_{\rm epoch}$ kind: greater_equal value: 2 # Magnitude range - col_name: mag - label: mag range + label: r kind: range value: [15, 30] @@ -86,7 +86,7 @@ dat_ext: # Rough pointing coverage - col_name: npoint3 - label: r"$n_{\rm pointing}$" + label: $n_{\rm pointing}$ kind: greater_equal value: 3 diff --git a/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py b/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py index 3ee2c88..4ced303 100644 --- a/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py +++ b/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py @@ -36,14 +36,24 @@ config = obj.read_config_set_params("config_mask.yaml") -hist_data_path = "magnitude_histograms_data.npz" - -test_only = False +test_only = True # %% +# Funcitons def get_data(obj, test_only=False): - - + """Get Data. + + Returns catalogue. + + Parameters + ---------- + obj : CalibrateCat instance + Instance of CalibrateCat class + test_only : bool, optional + If True, only load a subset of data for testing; + default is False. + + """ # Get data. Set load_into_memory to False for very large files dat, dat_ext = obj.read_cat(load_into_memory=False) @@ -58,6 +68,8 @@ def get_data(obj, test_only=False): def read_hist_data(hist_data_path): """ + Read Hist Data. + Read histogram data from npz file. Parameters @@ -72,6 +84,7 @@ def read_hist_data(hist_data_path): - 'counts': histogram counts - 'bins': bin edges - 'label': label for the histogram + """ loaded = np.load(hist_data_path, allow_pickle=True) hist_data = {} @@ -87,13 +100,29 @@ def read_hist_data(hist_data_path): return hist_data - def get_mask(masks, col_name): - + """Get Mask. + + Returns mask corresponding to col_name. + + Parameters + ---------- + masks : list + List of mask objects + col_name : str + Column name to identify the mask + Returns + ------- + list + Mask object + integer + Mask position in list + + """ # Get mask fomr masks with col_name = col_name - for mask in masks: + for idx, mask in enumerate(masks): if mask._col_name == col_name: - return mask + return mask, idx def compute_hist(masks, col_name, mask_cumul, mag, bins): @@ -125,12 +154,12 @@ def compute_hist(masks, col_name, mask_cumul, mag, bins): Number of valid data points after masking mask_cumul : array Updated cumulative mask array + """ - this_mask = get_mask(masks, col_name) + this_mask, _ = get_mask(masks, col_name) # First time: if mask_cumul is None: - print("Init mask_cumul") mask_cumul = this_mask._mask else: mask_cumul &= this_mask._mask @@ -143,14 +172,15 @@ def compute_hist(masks, col_name, mask_cumul, mag, bins): # Label string_buffer = StringIO() - this_mask.print_condition(string_buffer) + this_mask.print_condition(string_buffer, latex=True) label = string_buffer.getvalue().strip() if label == "": label = col_name + print("MKDEBUG", label) counts, bin_edges = np.histogram(my_mag, bins=bins) - return counts, bin_edges, label, n_valid, mask_cumul + return counts, bin_edges, rf"{label}", n_valid, mask_cumul def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): @@ -184,7 +214,68 @@ def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): color=color ) + +def plot_all_hists( + hist_data, + col_names, + figsize=10, + alpha=0.5, + color_map=None, + fraction=False, + ax=None, + out_path=None, +): + + if ax is None: + plt.figure() + fig, (ax) = plt.subplots( + 1, + 1, + figsize=(figsize, figsize) + ) + + counts0 = None + for col_name in col_names: + if col_name in hist_data: + + data = hist_data[col_name] + if fraction: + if counts0 is None: + counts0 = data["counts"] + counts = data["counts"] / counts0 + else: + counts = data["counts"] + + plot_hist( + counts, + data['bins'], + data['label'], + alpha=alpha, + ax=ax, + color=color_map[col_name] + ) + #print(f"{col_name}: n_valid = {data['n_valid']}") + + ax.set_xlabel('$r$') + ylabel = "fraction" if fraction else "number" + ax.set_ylabel(ylabel) + ax.set_xlim(17.5, 26.5) + if not fraction: + ax.legend() + + if out_path: + plt.tight_layout() + plt.savefig( + out_path, + dpi=150, + bbox_inches='tight' + ) + + # %% +# Main program +scenario = 1 +hist_data_path = f"magnitude_histograms_data_scenario-{scenario}.npz" if os.path.exists(hist_data_path): print(f"Histogram data file {hist_data_path} found.") @@ -201,8 +292,8 @@ def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): hist_data = None -# ## Masking # %% +# Masking # Get all masks, with or without dat, dat_ext masks, labels = sp_joint.get_masks_from_config( config, @@ -212,22 +303,89 @@ def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): ) # %% +# Combine mask according to scenario # List of basic masks to apply to all cases -masks_labels_basic = [ - "FLAGS", - "overlap", - "IMAFLAGS_ISO", - "mag", - "NGMIX_MOM_FAIL", - "NGMIX_ELL_PSFo_NOSHEAR_0", - "NGMIX_ELL_PSFo_NOSHEAR_1", - "4_Stars", - "8_Manual", - "64_r", - "1024_Maximask", -] + +masks_labels_basic = ["overlap", "mag", "64_r"] +col_names = ["basic masks"] + +if scenario == 0: + masks_labels_basic.extend([ + "FLAGS", + "IMAFLAGS_ISO", + "NGMIX_MOM_FAIL", + "NGMIX_ELL_PSFo_NOSHEAR_0", + "NGMIX_ELL_PSFo_NOSHEAR_1", + "4_Stars", + "8_Manual", + "1024_Maximask", + ]) + + col_names.extend(["N_EPOCH", "npoint3", "metacal"]) + +elif scenario == 1: + + col_names.extend([ + "IMAFLAGS_ISO", + "FLAGS", + "NGMIX_MOM_FAIL", + "NGMIX_ELL_PSFo_NOSHEAR_0", + "NGMIX_ELL_PSFo_NOSHEAR_1", + "4_Stars", + "8_Manual", + "1024_Maximask", + "N_EPOCH", + "npoint3", + "metacal", + ]) + + combine_cols = { + "ngmix failures": [ + "NGMIX_MOM_FAIL", + "NGMIX_ELL_PSFo_NOSHEAR_0", + "NGMIX_ELL_PSFo_NOSHEAR_1", + ] + } + +# %% +# Combine columns if specified. +# Remove old columns after combining. +if combine_cols is not None: + for new_col, old_cols in combine_cols.items(): + print(f"Combining columns {old_cols} into {new_col}") + # Create combined mask + old_masks = [] + idx_first = None + for col in old_cols: + mask, idx = get_mask(masks, col) + old_masks.append(mask) + if idx_first is None: + idx_first = idx + if dat is not None: + print(f"Creating combined mask for {new_col}") + masks_combined = sp_joint.Mask.from_list( + old_masks, + label=new_col, + verbose=obj._params["verbose"], + ) + else: + print(f"Creating dummy mask for {new_col} (for plot label)") + masks_combined = sp_joint.Mask( + new_col, + new_col, + kind="none", + ) + masks.insert(idx, masks_combined) + col_names.insert(idx, new_col) + + for old_mask, old_col in zip(old_masks, old_cols): + masks.remove(old_mask) + col_names.remove(old_col) + + print("After combining: masks =", [mask._col_name for mask in masks]) # %% +# Createe list of masks masks_basic = [] for mask in masks: if mask._col_name in masks_labels_basic: @@ -237,15 +395,15 @@ def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): print("Creating combined basic mask") masks_basic_combined = sp_joint.Mask.from_list( masks_basic, - label="combined", + label="basic masks", verbose=obj._params["verbose"], ) else: # Create dummy combined mask (for plot label) print("Creating dummy combined basic mask") masks_basic_combined = sp_joint.Mask( - "combined", - "combined", + "basic masks", + "basic masks", kind="none", ) @@ -260,6 +418,7 @@ def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): ) # %% +# Call metacal if data is available if dat is not None: cm = config["metacal"] gal_metacal = metacal( @@ -290,22 +449,17 @@ def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): # %% -# Plot magnitude histograms for various selection criteria - # Define magnitude bins mag_bins = np.arange(15, 30, 0.05) mag_centers = 0.5 * (mag_bins[:-1] + mag_bins[1:]) - # %% # Create figure with multiple subplots figsize = 10 alpha = 0.5 -col_names = ["combined", "N_EPOCH", "npoint3", "metacal"] - # Define explicit colors for each histogram -colors = ['C0', 'C1', 'C2', 'C3'] # Use matplotlib default color cycle +colors = [f'C{i}' for i in range(len(col_names))] # Use matplotlib default color cycle color_map = dict(zip(col_names, colors)) @@ -333,32 +487,38 @@ def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): 'label': label, 'n_valid': n_valid, } - +# %% +# Create plots +fig, axes = plt.subplots( + 1, + 2, + figsize=(2 * figsize, figsize) +) # Plot histogram data -plt.figure() -fig, (ax) = plt.subplots(1, 1, figsize=(figsize, figsize)) - -for col_name in col_names: - if col_name in hist_data: - data = hist_data[col_name] - plot_hist( - data['counts'], - data['bins'], - data['label'], - alpha=alpha, - ax=ax, - color=color_map[col_name] - ) - #print(f"{col_name}: n_valid = {data['n_valid']}") - -ax.set_xlabel('$r$') -ax.set_ylabel('Number') -ax.set_xlim(17.5, 26.5) -ax.legend() - +plot_all_hists( + hist_data, + col_names, + alpha=alpha, + color_map=color_map, + ax=axes[0], +) +plot_all_hists( + hist_data, + col_names, + alpha=alpha, + color_map=color_map, + fraction=True, + ax=axes[1], +) plt.tight_layout() -plt.savefig('magnitude_histograms.png', dpi=150, bbox_inches='tight') +out_path = f"magnitude_histograms_scenario-{scenario}.png" +plt.savefig( + out_path, + dpi=150, + bbox_inches='tight' +) +# %% # Save histogram data to file (only if we computed it) if dat is not None: np.savez( @@ -381,4 +541,9 @@ def plot_hist(counts, bins, label, alpha=1, ax=None, color=None): # %% if dat is not None: obj.close_hd5() + +# %% +for mask in masks: + mask.print_condition(sys.stdout, latex=True) + # %% diff --git a/src/sp_validation/run_joint_cat.py b/src/sp_validation/run_joint_cat.py index c9d5965..e30bc43 100644 --- a/src/sp_validation/run_joint_cat.py +++ b/src/sp_validation/run_joint_cat.py @@ -1365,32 +1365,34 @@ def print_stats(self, num_obj, f_out=None): sf = f"{self._num_ok/num_obj:10.2%}" self.print_strings(self._col_name, self._label, si, sf, f_out=f_out) - def get_sign(self): + def get_sign(self, latex=False): sign = None - if self._kind =="equal": - sign = "=" - elif self._kind =="not_equal": - sign = "!=" - elif self._kind =="greater_equal": - sign = ">=" - elif self._kind =="smaller_equal": - sign = "<=" + if self._kind == "equal": + sign = "$=$" if latex else "=" + elif self._kind == "not_equal": + sign = "$\ne$" if latex else "!=" + elif self._kind in ("greater_equal", "range"): + sign = "$\leq$" if latex else ">=" + elif self._kind == "smaller_equal": + sign = "$\geq$" if latex else "<=" return sign - def print_condition(self, f_out): + def print_condition(self, f_out, latex=False): if self._value is None: return "" - sign = self.get_sign() + sign = self.get_sign(latex=latex) + + name = self._label if latex else self._col_name if sign is not None: - print(f"{self._col_name} {sign} {self._value}", file=f_out) - + print(f"{name} {sign} {self._value}", file=f_out) + if self._kind == "range": - print(f"{self._value[0]} <= {self._col_name} <= {self._value[1]}", file=f_out) - + print(f"{self._value[0]} {sign} {name} {sign} {self._value[1]}", file=f_out) + def print_summary(self, f_out): print(f"[{self._label}]\t\t\t", file=f_out, end="") self.print_condition(f_out) From f4f0bb99df63569837bfa19cd3bdad4a76791791 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Mon, 17 Nov 2025 13:13:17 +0100 Subject: [PATCH 3/7] hist plot, basic: mask size and SNR separately --- .../cosmo_val/catalog_paper_plot/hist_mag.py | 23 ++++++++++++++++ src/sp_validation/basic.py | 26 +++++++++++++++++++ 2 files changed, 49 insertions(+) diff --git a/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py b/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py index 4ced303..7d202d7 100644 --- a/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py +++ b/notebooks/cosmo_val/catalog_paper_plot/hist_mag.py @@ -417,6 +417,25 @@ def plot_all_hists( kind="none", ) +# %% +def get_info_for_metacal_masking(dat, mask, prefix = "NGMIX", name_shear = "NOSHEAR"): + + res = {} + + res["flag"] = dat[mask][f"{prefix}_FLAGS_{name_shear}"] + + for key in ("flux", "flux_err", "T"): + res[key] = dat[mask][f"{prefix}_{key.upper()}_{name_shear}"] + res["Tpsf"] = dat[mask][f"{prefix}_Tpsf_{name_shear}"] + + return res + +# %% +if dat is not None: + cm = config["metacal"] + + + # %% # Call metacal if data is available if dat is not None: @@ -547,3 +566,7 @@ def plot_all_hists( mask.print_condition(sys.stdout, latex=True) # %% +# print number of valid objects and name +for data in hist_data + +# %% diff --git a/src/sp_validation/basic.py b/src/sp_validation/basic.py index 3720004..fc96230 100644 --- a/src/sp_validation/basic.py +++ b/src/sp_validation/basic.py @@ -597,3 +597,29 @@ def jackknif_weighted_average2( all_est = np.array(all_est) return np.mean(all_est), np.std(all_est) + + +def mask_gal_size(T, Tpsf, rel_size_min, rel_size_max, size_corr_ell=False, g1=None, g2=None): + + Tr_tmp = T + if size_corr_ell: + Tr_tmp *= ( + (1 - g1 **2 + g2 ** 2) / (1 + g1 ** 2 + g2 **2) + ) + + mask = ( + (Tr_tmp / Tpsf > rel_size_min) + & (Tr_tmp / Tpsf < rel_size_max) + ) + + return mask + + +def mask_gal_SNR(SNR, snr_min, snr_max): + + mask = ( + (SNR > snr_min) + & (SNR < snr_max) + ) + + return mask From cc43bd1fed5a139df84bee84f3ffec649d77007f Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Fri, 5 Dec 2025 11:37:25 +0100 Subject: [PATCH 4/7] Added scripts to plot star-gal correlations --- .../catalog_paper_plot/check_gaia.ipynb | 73 +++++++++ .../download_gaia_catalogues.py | 100 ++++++++++++ .../catalog_paper_plot/plot_gaia_sky.py | 154 ++++++++++++++++++ .../catalog_paper_plot/run_cosmo_val_GAIA.py | 42 +++++ .../run_cosmo_val_gammat.py | 74 +++++++++ .../run_cosmo_val_lambda.py | 82 ++++++++++ 6 files changed, 525 insertions(+) create mode 100644 notebooks/cosmo_val/catalog_paper_plot/check_gaia.ipynb create mode 100644 notebooks/cosmo_val/catalog_paper_plot/download_gaia_catalogues.py create mode 100644 notebooks/cosmo_val/catalog_paper_plot/plot_gaia_sky.py create mode 100644 notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_GAIA.py create mode 100644 notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_gammat.py create mode 100644 notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_lambda.py diff --git a/notebooks/cosmo_val/catalog_paper_plot/check_gaia.ipynb b/notebooks/cosmo_val/catalog_paper_plot/check_gaia.ipynb new file mode 100644 index 0000000..0de529a --- /dev/null +++ b/notebooks/cosmo_val/catalog_paper_plot/check_gaia.ipynb @@ -0,0 +1,73 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "8049856e", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d7f7568b", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.io import fits" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "62854cbc", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/automnt/n17data/mkilbing/astro/Runs/flexion/simuls_MAB\n" + ] + } + ], + "source": [ + "!pwd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c311b19", + "metadata": {}, + "outputs": [], + "source": [ + "cat = fits.getdata(\"gaia2.fits\", 1)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "base", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/cosmo_val/catalog_paper_plot/download_gaia_catalogues.py b/notebooks/cosmo_val/catalog_paper_plot/download_gaia_catalogues.py new file mode 100644 index 0000000..746a641 --- /dev/null +++ b/notebooks/cosmo_val/catalog_paper_plot/download_gaia_catalogues.py @@ -0,0 +1,100 @@ +# %% +# Calibrate comprehensive catalogue + +# %% +from IPython import get_ipython + +# %% +# enable autoreload for interactive sessions +ipython = get_ipython() +if ipython is not None: + ipython.run_line_magic("reload_ext", "autoreload") + ipython.run_line_magic("autoreload", "2") + ipython.run_line_magic("reload_ext", "log_cell_time") + +import numpy as np +import time +from astroquery.gaia import Gaia +import astropy.units +from astropy.coordinates import SkyCoord +from textwrap import dedent + + +def get_bins(results, key="phot_g_mean_mag", n_bins=3): + + quantiles_raw = np.percentile(results[key], np.linspace(0, 100, n_bins + 1)) + # round + quantiles = np.round(quantiles_raw, 1) + + print("Magnitude bin edges:") + subsets = [] + for i in range(n_bins): + this_subset = results[ + (results[key] >= quantiles[i]) & (results[key] < quantiles[i + 1]) + ] + subsets.append(this_subset) + + print( + f" Bin {i+1}: {quantiles[i]:.2f}, {quantiles[i+1]:.2f}, N={len(this_subset)}" + ) + + return subsets, quantiles + + +def do_query(do_async=True, nmax=3_000_000): + + # Define query + query = dedent( + f""" + SELECT TOP {nmax} ra, dec, phot_g_mean_mag, ruwe + FROM gaiadr3.gaia_source + WHERE dec > 30 + AND dec < 75 + AND phot_g_mean_mag < 20 + AND phot_g_mean_mag IS NOT NULL + AND ruwe < 1.4 + ORDER BY random_index + """ + ).strip() + + + # Launch the query + if not do_async: + print("Retrieveing GAIA data (sync)") + job = Gaia.launch_job(query) + else: + print("Retrieveing GAIA data (async)") + job = Gaia.launch_job_async(query) + + while not job.is_finished(): + print(f"Job status: {job.get_phase()}") + time.sleep(10) + + results = job.get_results() + + print(f"Retrieved {len(results)} sources") + print(results[:10]) # Show first 10 results + + return results + + +def write_subsets(subsets, quantiles): + + for idx, subset in enumerate(subsets): + + if idx == 0: + out_path = f"gaia_stars_g_smaller_{quantiles[idx + 1]}.fits" + elif idx == len(subsets) - 1: + out_path = f"gaia_stars_g_larger_{quantiles[idx]}.fits" + else: + out_path = f"gaia_stars_g_in_{quantiles[idx]}_{quantiles[idx + 1]}.fits" + print(f"Writing {len(subset)} objects to {out_path}") + subset.write(out_path, format="fits", overwrite=True) + + +# %% Main program +results = do_query(do_async=True) + +subsets, quantiles = get_bins(results) + +write_subsets(subsets, quantiles) diff --git a/notebooks/cosmo_val/catalog_paper_plot/plot_gaia_sky.py b/notebooks/cosmo_val/catalog_paper_plot/plot_gaia_sky.py new file mode 100644 index 0000000..ebe3199 --- /dev/null +++ b/notebooks/cosmo_val/catalog_paper_plot/plot_gaia_sky.py @@ -0,0 +1,154 @@ +"""Plot GAIA stars on sky maps by magnitude bins.""" + +import glob +import os +import re +import numpy as np +from astropy.table import Table +import matplotlib.pyplot as plt +from sp_validation.plots import FootprintPlotter + + +def load_gaia_fits(file_path): + """Load GAIA FITS file and return RA, DEC arrays. + + Parameters + ---------- + file_path : str + Path to FITS file + + Returns + ------- + ra : np.ndarray + Right ascension values + dec : np.ndarray + Declination values + n_objects : int + Number of objects in file + """ + table = Table.read(file_path, format='fits') + print(f"Loaded {len(table)} objects from {file_path}") + return table['ra'], table['dec'], len(table) + + +def extract_label_from_filename(filename): + """Extract magnitude information from filename to create label. + + Parameters + ---------- + filename : str + FITS filename + + Returns + ------- + str + Label describing the magnitude bin + """ + basename = os.path.basename(filename) + + # Match pattern: g_smaller_XX.X + match = re.search(r'g_smaller_([0-9.]+)', basename) + if match: + mag = match.group(1) + return f"G < {mag} (bright)" + + # Match pattern: g_in_XX.X_YY.Y + match = re.search(r'g_in_([0-9.]+)_([0-9.]+)', basename) + if match: + mag1, mag2 = match.group(1), match.group(2) + return f"{mag1} ≤ G < {mag2} (medium)" + + # Match pattern: g_larger_XX.X + match = re.search(r'g_larger_([0-9.]+)', basename) + if match: + mag = match.group(1) + return f"G ≥ {mag} (faint)" + + # Fallback + return basename.replace('.fits', '') + + +def main(): + """Create sky plots for GAIA stars in three magnitude bins.""" + + # Auto-detect GAIA files using patterns + pattern_smaller = "gaia_stars_g_smaller_*.fits" + pattern_in = "gaia_stars_g_in_*.fits" + pattern_larger = "gaia_stars_g_larger_*.fits" + + files = [] + files.extend(sorted(glob.glob(pattern_smaller))) + files.extend(sorted(glob.glob(pattern_in))) + files.extend(sorted(glob.glob(pattern_larger))) + + if not files: + print("ERROR: No GAIA FITS files found matching patterns:") + print(f" - {pattern_smaller}") + print(f" - {pattern_in}") + print(f" - {pattern_larger}") + return + + print(f"Found {len(files)} GAIA FITS files:") + for f in files: + print(f" - {f}") + print() + + # Generate labels from filenames + labels = [extract_label_from_filename(f) for f in files] + + # Initialize the footprint plotter + plotter = FootprintPlotter(nside_coverage=32, nside_map=2048) + + # Process each magnitude bin + hsp_maps = [] + for file_path, label in zip(files, labels): + print(f"Processing: {label}") + + # Load data + ra, dec, n_objects = load_gaia_fits(file_path) + + # Create healsparse map + hsp_map = plotter.create_hsp_map(ra, dec) + hsp_maps.append(hsp_map) + + # Plot all regions (NGC, SGC, fullsky) + plotter.plot_all_regions( + hsp_map, + outbase=f"gaia_sky_{file_path.replace('.fits', '')}" + ) + print(f"Created plots for {label}") + + # Create combined plot showing all bins + n_files = len(files) + if n_files > 0: + print("Creating combined plot") + + fig, axes = plt.subplots(1, n_files, figsize=(10 * n_files, 10)) + + # Handle case of single file + if n_files == 1: + axes = [axes] + + for idx, (hsp_map, label) in enumerate(zip(hsp_maps, labels)): + ax = axes[idx] + + # Use fullsky region parameters + region = plotter._regions['fullsky'] + + projection = plotter.plot_area( + hsp_map, + ra_0=region['ra_0'], + extend=region['extend'], + vmax=region['vmax'], + title=label + ) + + plt.tight_layout() + plt.savefig('gaia_sky_combined_all_bins.png', dpi=150, bbox_inches='tight') + print("Saved combined plot: gaia_sky_combined_all_bins.png") + + print("All plots created successfully!") + + +if __name__ == "__main__": + main() diff --git a/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_GAIA.py b/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_GAIA.py new file mode 100644 index 0000000..af58f0e --- /dev/null +++ b/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_GAIA.py @@ -0,0 +1,42 @@ +# %% +from IPython import get_ipython + +ipython = get_ipython() + +# enable autoreload for interactive sessions +if ipython is not None: + ipython.run_line_magic("load_ext", "autoreload") + ipython.run_line_magic("autoreload", "2") + ipython.run_line_magic("load_ext", "log_cell_time") + +import matplotlib.pyplot as plt # noqa: E402, F401 +import numpy as np # noqa: E402, F401 +import sys # noqa: E402 +sys.path.insert(0, '/home/mkilbing/astro/repositories/gitlab.euclid-sgs/FDQA/rho_tau_stats') +from cosmo_val import CosmologyValidation # noqa: E402 + +# enable inline plotting for interactive sessions +# (must be done *after* importing package that sets agg backend) +if ipython is not None: + ipython.run_line_magic("matplotlib", "inline") + + +# %% +cv = CosmologyValidation( + versions=["SP_v1.4.6_GAIA_bright", "SP_v1.4.6_GAIA_medium", "SP_v1.4.6_GAIA_faint"], + data_base_dir="/n17data/mkilbing/astro/data/", + catalog_config="./cat_config.yaml", + output_dir="./output", + theta_min=1.0, + theta_max=250.0, + nbins=20, + cov_estimate_method='jk', + theta_min_plot=0.8, + theta_max_plot=260.0, + rho_tau_method='emcee', + n_cov=1, +) + +# %% +cv.plot_gammat_around_stars(offset=0.05, gammax=True) +cv.plot_gammat_around_stars(offset=0.05, gammax=True, logy=True) diff --git a/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_gammat.py b/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_gammat.py new file mode 100644 index 0000000..4ddb1c3 --- /dev/null +++ b/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_gammat.py @@ -0,0 +1,74 @@ +# %% +from IPython import get_ipython + +ipython = get_ipython() + +# enable autoreload for interactive sessions +if ipython is not None: + ipython.run_line_magic("load_ext", "autoreload") + ipython.run_line_magic("autoreload", "2") + ipython.run_line_magic("reload_ext", "log_cell_time") + +import matplotlib.pyplot as plt +import numpy as np +import sys +import os +sys.path.insert(0, '/home/mkilbing/astro/repositories/gitlab.euclid-sgs/FDQA/rho_tau_stats') +from cosmo_val import CosmologyValidation + + +def rename_output(output_bases, output_dir, suff, version_string, rand_str): + + for base in output_bases: + old_path = f"{output_dir}/plots/{base}{suff}" + new_path = f"{output_dir}/plots/{base}_{version_string}{rand_str}{suff}" + print(f"mv {old_path} -> {new_path}") + os.rename(old_path, new_path) + + +# enable inline plotting for interactive sessions +# (must be done *after* importing package that sets agg backend) +if ipython is not None: + ipython.run_line_magic("matplotlib", "inline") + +versions_146 = ["SP_v1.4.6", "SP_v1.4.6_leak_corr"] +versions_var = ["SP_v1.4.5", "SP_v1.4.6", "SP_v1.4.7", "SP_v1.4.8"] + +versions = versions_146 + +output_dir = "./output" + +# %% +cv = CosmologyValidation( + versions=versions_146, + data_base_dir="/n17data/mkilbing/astro/data/", + catalog_config="./cat_config.yaml", + output_dir=output_dir, + theta_min=1.0, + theta_max=250.0, + nbins=20, + cov_estimate_method='jk', + theta_min_plot=0.8, + theta_max_plot=260.0, + rho_tau_method='emcee', + n_cov=1, + star_weight_type="uniform", +) + + +output_bases = [ + "gammat_around_stars_lin_non_tomographic", + "gammat_around_stars_log_non_tomographic", +] +suff = ".png" +version_string = "_".join(versions) + +# %% +cv.plot_gammat_around_stars(offset=0.05, gammax=True, wo_rand_subtr=True) +cv.plot_gammat_around_stars(offset=0.05, gammax=True, logy=True, wo_rand_subtr=True) +rename_output(output_bases, output_dir, suff, version_string, "_wo_rand_subtr") + +# %% +cv.plot_gammat_around_stars(offset=0.05, gammax=True) +cv.plot_gammat_around_stars(offset=0.05, gammax=True, logy=True) +rename_output(output_bases, output_dir, suff, version_string, "") diff --git a/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_lambda.py b/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_lambda.py new file mode 100644 index 0000000..7f7d566 --- /dev/null +++ b/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_lambda.py @@ -0,0 +1,82 @@ +# %% +from IPython import get_ipython + +ipython = get_ipython() + +# enable autoreload for interactive sessions +if ipython is not None: + ipython.run_line_magic("load_ext", "autoreload") + ipython.run_line_magic("autoreload", "2") + ipython.run_line_magic("load_ext", "log_cell_time") + +import sys +import os +import matplotlib.pyplot as plt +import numpy as np + +sys.path.insert(0, '/home/mkilbing/astro/repositories/gitlab.euclid-sgs/FDQA/rho_tau_stats') +from cosmo_val import CosmologyValidation # noqa: E402 + + +def rename_output(output_bases, output_dir, suff, version_string, rand_str): + + for base in output_bases: + old_path = f"{output_dir}/plots/{base}{suff}" + new_path = f"{output_dir}/plots/{base}_{version_string}{rand_str}{suff}" + print(f"mv {old_path} -> {new_path}") + os.rename(old_path, new_path) + + +# enable inline plotting for interactive sessions +# (must be done *after* importing package that sets agg backend) +if ipython is not None: + ipython.run_line_magic("matplotlib", "inline") + +# Different options +versions_146 = ["SP_v1.4.6", "SP_v1.4.6_leak_corr"] +versions_var = ["SP_v1.4.5", "SP_v1.4.6", "SP_v1.4.7", "SP_v1.4.8"] + +# We use this combination of versions +versions = versions_146 + +output_dir = "./output" + +theta_min = 1.0 +theta_max = 300.0 +plot_range_fac = 1.1 + +# %% +cv = CosmologyValidation( + versions=versions, + data_base_dir="/n17data/mkilbing/astro/data/", + catalog_config="./cat_config.yaml", + output_dir=output_dir, + theta_min=theta_min, + theta_max=theta_max, + nbins=15, + cov_estimate_method='jk', + theta_min_plot=theta_min / plot_range_fac, + theta_max_plot=theta_max * plot_range_fac, + rho_tau_method='emcee', + n_cov=1, + star_weight_type="uniform", + random_multiple=5, +) + +output_bases = [ + "gammat_stars_around_galaxies_lin_non_tomographic", + "gammat_stars_around_galaxies_log_non_tomographic", +] +suff = ".png" +version_string = "_".join(versions) + + +cv.plot_gammat_stars_around_galaxies(offset=0.025, gammax=True, wo_rand_subtr=True) +cv.plot_gammat_stars_around_galaxies(offset=0.025, gammax=True, logy=True, wo_rand_subtr=True) +rename_output(output_bases, output_dir, suff, version_string, "_wo_rand_subtr") + +# %% +cv.plot_gammat_stars_around_galaxies(offset=0.025, gammax=True) +cv.plot_gammat_stars_around_galaxies(offset=0.025, gammax=True, logy=True) +rename_output(output_bases, output_dir, suff, version_string, "") + From 52ca866cc02414609ddc08fe0bc219bafa2d60fb Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Fri, 5 Dec 2025 11:37:57 +0100 Subject: [PATCH 5/7] cat yaml config file --- notebooks/cosmo_val/cat_config.yaml | 192 +++++++++++++++++++++++++++- 1 file changed, 191 insertions(+), 1 deletion(-) diff --git a/notebooks/cosmo_val/cat_config.yaml b/notebooks/cosmo_val/cat_config.yaml index 67a0cad..47e26d6 100644 --- a/notebooks/cosmo_val/cat_config.yaml +++ b/notebooks/cosmo_val/cat_config.yaml @@ -694,6 +694,8 @@ SP_v1.4.5: shear: path: ../../../../../../sguerrini/unions_shapepipe_cut_struc_2024_v1.4.5_rerun.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec e1_col: e1 e2_col: e2 e1_PSF_col: e1_PSF @@ -722,6 +724,7 @@ SP_v1.4.5: PSF_flag: "FLAG_PSF_HSM" star_flag: "FLAG_STAR_HSM" square_size: True + patch_number: 100 SP_v1.4.6: pipeline: SP @@ -738,6 +741,8 @@ SP_v1.4.6: shear: path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec e1_col: e1 e2_col: e2 e1_PSF_col: e1_PSF @@ -766,6 +771,179 @@ SP_v1.4.6: PSF_flag: "FLAG_PSF_HSM" star_flag: "FLAG_STAR_HSM" square_size: True + patch_number: 100 + +SP_v1.4.7: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + ls: dashed + colour: darkorchid + marker: d + getdist_colour: 0.0, 0.5, 1.0 + cov_th: + A: 2420.2651014497287 #deg^2 + n_e: 7.040818382014773 #arcmin-2 + n_psf: 0.752316232272063 #arcmin-2 + sigma_e: 0.30961528707207325 + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.7/unions_shapepipe_cut_struc_2024_v1.4.7.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + e1_PSF_col: e1_PSF + e2_PSF_col: e2_PSF + w_col: w_des + R: 1.0 + redshift_distr: CFIS/v1.0/nz/dndz_SP_A.txt + psf: + path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + ra_col: RA + dec_col: Dec + e1_PSF_col: E1_PSF_HSM + e2_PSF_col: E2_PSF_HSM + e1_star_col: E1_STAR_HSM + e2_star_col: E2_STAR_HSM + PSF_size: SIGMA_PSF_HSM + star_size: SIGMA_STAR_HSM + PSF_flag: "FLAG_PSF_HSM" + star_flag: "FLAG_STAR_HSM" + square_size: True + patch_number: 100 + +SP_v1.4.6_GAIA_bright: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + ls: dashed + colour: green + marker: s + getdist_colour: 0.0, 0.5, 1.0 + cov_th: + A: 2420.2651014497287 #deg^2 + n_e: 7.040818382014773 #arcmin-2 + n_psf: 0.752316232272063 #arcmin-2 + sigma_e: 0.30961528707207325 + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + e1_PSF_col: e1_PSF + e2_PSF_col: e2_PSF + w_col: w_des + R: 1.0 + mask: /n09data/guerrini/glass_mock/mask_v1.4.5_nside8192.fits + redshift_distr: ../../../sguerrini/UNIONS/WL/nz/v1.4.6/nz_SP_v1.4.6_A.txt + psf: + path: ../../../GAIA/gaia_stars_g_bright.fits + hdu: 1 + ra_col: RA + dec_col: Dec + square_size: True + patch_number: 100 + +SP_v1.4.6_GAIA_medium: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + ls: dashed + colour: darkgreen + marker: o + getdist_colour: 0.0, 0.5, 1.0 + cov_th: + A: 2420.2651014497287 #deg^2 + n_e: 7.040818382014773 #arcmin-2 + n_psf: 0.752316232272063 #arcmin-2 + sigma_e: 0.30961528707207325 + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + e1_PSF_col: e1_PSF + e2_PSF_col: e2_PSF + w_col: w_des + R: 1.0 + mask: /n09data/guerrini/glass_mock/mask_v1.4.5_nside8192.fits + redshift_distr: ../../../sguerrini/UNIONS/WL/nz/v1.4.6/nz_SP_v1.4.6_A.txt + psf: + path: ../../../GAIA/gaia_stars_g_medium.fits + hdu: 1 + ra_col: RA + dec_col: Dec + square_size: True + patch_number: 100 + +SP_v1.4.6_GAIA_faint: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + ls: dashed + colour: lightgreen + marker: d + getdist_colour: 0.0, 0.5, 1.0 + cov_th: + A: 2420.2651014497287 #deg^2 + n_e: 7.040818382014773 #arcmin-2 + n_psf: 0.752316232272063 #arcmin-2 + sigma_e: 0.30961528707207325 + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + e1_PSF_col: e1_PSF + e2_PSF_col: e2_PSF + w_col: w_des + R: 1.0 + mask: /n09data/guerrini/glass_mock/mask_v1.4.5_nside8192.fits + redshift_distr: ../../../sguerrini/UNIONS/WL/nz/v1.4.6/nz_SP_v1.4.6_A.txt + psf: + path: ../../../GAIA/gaia_stars_g_faint.fits + hdu: 1 + ra_col: RA + dec_col: Dec + square_size: True + patch_number: 100 + + +SP_v1.4.7_GAIA_bright: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + ls: dashed + colour: darkorchid + marker: d + getdist_colour: 0.0, 0.5, 1.0 + cov_th: + A: 2420.2651014497287 #deg^2 + n_e: 7.040818382014773 #arcmin-2 + n_psf: 0.752316232272063 #arcmin-2 + sigma_e: 0.30961528707207325 + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.7/unions_shapepipe_cut_struc_2024_v1.4.7.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1 + e2_col: e2 + e1_PSF_col: e1_PSF + e2_PSF_col: e2_PSF + w_col: w_des + R: 1.0 + redshift_distr: CFIS/v1.0/nz/dndz_SP_A.txt + psf: + path: ../../../GAIA/gaia_stars_g_smaller_16.0.fits + hdu: 1 + ra_col: RA + dec_col: dec + square_size: True + patch_number: 100 SP_v1.4.8: pipeline: SP @@ -782,6 +960,8 @@ SP_v1.4.8: shear: path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.8/unions_shapepipe_cut_struc_2024_v1.4.8.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec e1_col: e1 e2_col: e2 e1_PSF_col: e1_PSF @@ -791,7 +971,7 @@ SP_v1.4.8: mask: /n09data/guerrini/glass_mock/mask_v1.4.5_nside8192.fits redshift_distr: ../../../sguerrini/UNIONS/WL/nz/v1.4.8/nz_SP_v1.4.8_A.txt star: - path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits + path: ../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits ra_col: RA dec_col: Dec e1_col: e1 @@ -810,6 +990,7 @@ SP_v1.4.8: PSF_flag: "FLAG_PSF_HSM" star_flag: "FLAG_STAR_HSM" square_size: True + patch_number: 100 SP_v1.4.5.A: pipeline: SP @@ -998,6 +1179,8 @@ SP_v1.4.6_leak_corr: shear: path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec e1_col: e1_leak_corrected e2_col: e2_leak_corrected e1_PSF_col: e1_PSF @@ -1026,6 +1209,7 @@ SP_v1.4.6_leak_corr: PSF_flag: "FLAG_PSF_HSM" star_flag: "FLAG_STAR_HSM" square_size: True + patch_number: 100 SP_v1.4.7_leak_corr: pipeline: SP @@ -1042,6 +1226,8 @@ SP_v1.4.7_leak_corr: shear: path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.7/unions_shapepipe_cut_struc_2024_v1.4.7.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec e1_col: e1_leak_corrected e2_col: e2_leak_corrected e1_PSF_col: e1_PSF @@ -1069,6 +1255,7 @@ SP_v1.4.7_leak_corr: PSF_flag: "FLAG_PSF_HSM" star_flag: "FLAG_STAR_HSM" square_size: True + patch_number: 100 SP_v1.4.8_leak_corr: pipeline: SP @@ -1388,6 +1575,8 @@ SP_v1.4.5_leak_corr: shear: path: ../../../../../../sguerrini/unions_shapepipe_cut_struc_2024_v1.4.5_rerun.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec e1_col: e1_leak_corrected e2_col: e2_leak_corrected e1_PSF_col: e1_PSF @@ -1416,6 +1605,7 @@ SP_v1.4.5_leak_corr: PSF_flag: "FLAG_PSF_HSM" star_flag: "FLAG_STAR_HSM" square_size: True + patch_number: 100 SP_v1.5.4: pipeline: SP From dd962f34d4324190923e219d17ea385ac7a9b6a1 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 17 Dec 2025 15:44:10 +0100 Subject: [PATCH 6/7] cosmo val config file: added v1.4.6 test cases --- notebooks/cosmo_val/cat_config.yaml | 104 ++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) diff --git a/notebooks/cosmo_val/cat_config.yaml b/notebooks/cosmo_val/cat_config.yaml index 47e26d6..686dd5e 100644 --- a/notebooks/cosmo_val/cat_config.yaml +++ b/notebooks/cosmo_val/cat_config.yaml @@ -992,6 +992,110 @@ SP_v1.4.8: square_size: True patch_number: 100 +# For additive bias only +SP_v1.4.6_uncal: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: w_des + R: 1.0 + psf: + path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 + +SP_v1.4.6_uncal_w_iv: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: w_iv + R: 1.0 + psf: + path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 + +SP_v1.4.6_uncal_w_1: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: one + R: 1.0 + psf: + path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 + +SP_v1.4.5_uncal: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.5/unions_shapepipe_cut_struc_2024_v1.4.5.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: w_des + R: 1.0 + psf: + path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 + +SP_v1.4.7_uncal: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.7/unions_shapepipe_cut_struc_2024_v1.4.7.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: w_des + R: 1.0 + psf: + path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 + +SP_v1.4.8_uncal: + pipeline: SP + subdir: CFIS/v1.0/ShapePipe/ + shear: + path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.8/unions_shapepipe_cut_struc_2024_v1.4.8.fits + covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt + ra_col: RA + dec_col: Dec + e1_col: e1_uncal + e2_col: e2_uncal + w_col: w_des + R: 1.0 + psf: + path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + hdu: 1 + patch_number: 100 + + SP_v1.4.5.A: pipeline: SP subdir: CFIS/v1.0/ShapePipe/ From 6d3f5daa09bc08070fd5c9ffeb361bbc9ff14557 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 17 Dec 2025 20:05:00 +0100 Subject: [PATCH 7/7] updated cosmo_val config to new path structure --- notebooks/cosmo_val/cat_config.yaml | 43 +++++++++++-------- .../run_cosmo_val_lambda.py | 42 ++++++++++++------ 2 files changed, 54 insertions(+), 31 deletions(-) diff --git a/notebooks/cosmo_val/cat_config.yaml b/notebooks/cosmo_val/cat_config.yaml index f3ca111..a1a6b75 100644 --- a/notebooks/cosmo_val/cat_config.yaml +++ b/notebooks/cosmo_val/cat_config.yaml @@ -616,6 +616,8 @@ SP_v1.4.6: covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt path: v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + ra_col: RA + dec_col: Dec w_col: w_des e1_col: e1 e1_col_corrected: e1_leak_corrected @@ -629,6 +631,7 @@ SP_v1.4.6: 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 pipeline: SP @@ -661,6 +664,8 @@ SP_v1.4.7: covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt path: v1.4.7/unions_shapepipe_cut_struc_2024_v1.4.7.fits redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + ra_col: RA + dec_col: Dec w_col: w_des e1_col: e1 e1_col_corrected: e1_leak_corrected @@ -706,6 +711,8 @@ SP_v1.4.8: covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt path: v1.4.8/unions_shapepipe_cut_struc_2024_v1.4.8.fits redshift_path: /n17data/mkilbing/astro/data/CFIS/v1.0/nz/dndz_SP_A.txt + ra_col: RA + dec_col: Dec w_col: w_des e1_col: e1 e1_col_corrected: e1_leak_corrected @@ -723,9 +730,9 @@ SP_v1.4.8: # For additive bias only SP_v1.4.6_uncal: pipeline: SP - subdir: CFIS/v1.0/ShapePipe/ + subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + path: v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -734,15 +741,15 @@ SP_v1.4.6_uncal: w_col: w_des R: 1.0 psf: - path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 SP_v1.4.6_uncal_w_iv: pipeline: SP - subdir: CFIS/v1.0/ShapePipe/ + subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + path: v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -751,15 +758,15 @@ SP_v1.4.6_uncal_w_iv: w_col: w_iv R: 1.0 psf: - path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 SP_v1.4.6_uncal_w_1: pipeline: SP - subdir: CFIS/v1.0/ShapePipe/ + subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits + path: v1.4.6/unions_shapepipe_cut_struc_2024_v1.4.6.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -768,15 +775,15 @@ SP_v1.4.6_uncal_w_1: w_col: one R: 1.0 psf: - path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 SP_v1.4.5_uncal: pipeline: SP - subdir: CFIS/v1.0/ShapePipe/ + subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.5/unions_shapepipe_cut_struc_2024_v1.4.5.fits + path: v1.4.5/unions_shapepipe_cut_struc_2024_v1.4.5.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -785,15 +792,15 @@ SP_v1.4.5_uncal: w_col: w_des R: 1.0 psf: - path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 SP_v1.4.7_uncal: pipeline: SP - subdir: CFIS/v1.0/ShapePipe/ + subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.7/unions_shapepipe_cut_struc_2024_v1.4.7.fits + path: v1.4.7/unions_shapepipe_cut_struc_2024_v1.4.7.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -802,15 +809,15 @@ SP_v1.4.7_uncal: w_col: w_des R: 1.0 psf: - path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 SP_v1.4.8_uncal: pipeline: SP - subdir: CFIS/v1.0/ShapePipe/ + subdir: /n17data/UNIONS/WL/v1.4.x shear: - path: ../../../../../../UNIONS/WL/v1.4.x/v1.4.8/unions_shapepipe_cut_struc_2024_v1.4.8.fits + path: v1.4.8/unions_shapepipe_cut_struc_2024_v1.4.8.fits covmat_file: ./covs/shapepipe_A/cov_shapepipe_A.txt ra_col: RA dec_col: Dec @@ -819,7 +826,7 @@ SP_v1.4.8_uncal: w_col: w_des R: 1.0 psf: - path: ../../../../../../UNIONS/WL/v1.4.x/unions_shapepipe_psf_2024_v1.4.a.fits + path: unions_shapepipe_psf_2024_v1.4.a.fits hdu: 1 patch_number: 100 diff --git a/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_lambda.py b/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_lambda.py index 7f7d566..9ce7322 100644 --- a/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_lambda.py +++ b/notebooks/cosmo_val/catalog_paper_plot/run_cosmo_val_lambda.py @@ -15,16 +15,7 @@ import numpy as np sys.path.insert(0, '/home/mkilbing/astro/repositories/gitlab.euclid-sgs/FDQA/rho_tau_stats') -from cosmo_val import CosmologyValidation # noqa: E402 - - -def rename_output(output_bases, output_dir, suff, version_string, rand_str): - - for base in output_bases: - old_path = f"{output_dir}/plots/{base}{suff}" - new_path = f"{output_dir}/plots/{base}_{version_string}{rand_str}{suff}" - print(f"mv {old_path} -> {new_path}") - os.rename(old_path, new_path) +from cosmo_val import CosmologyValidation, rename_output # noqa: E402 # enable inline plotting for interactive sessions @@ -33,7 +24,7 @@ def rename_output(output_bases, output_dir, suff, version_string, rand_str): ipython.run_line_magic("matplotlib", "inline") # Different options -versions_146 = ["SP_v1.4.6", "SP_v1.4.6_leak_corr"] +versions_146 = ["SP_v1.4.6"] #, "SP_v1.4.6_leak_corr"] versions_var = ["SP_v1.4.5", "SP_v1.4.6", "SP_v1.4.7", "SP_v1.4.8"] # We use this combination of versions @@ -45,11 +36,12 @@ def rename_output(output_bases, output_dir, suff, version_string, rand_str): theta_max = 300.0 plot_range_fac = 1.1 + # %% cv = CosmologyValidation( versions=versions, - data_base_dir="/n17data/mkilbing/astro/data/", catalog_config="./cat_config.yaml", + data_base_dir="/", output_dir=output_dir, theta_min=theta_min, theta_max=theta_max, @@ -66,17 +58,41 @@ def rename_output(output_bases, output_dir, suff, version_string, rand_str): output_bases = [ "gammat_stars_around_galaxies_lin_non_tomographic", "gammat_stars_around_galaxies_log_non_tomographic", + "dgammat_stars_around_galaxies_lin_non_tomographic", + "dsize_stars_around_galaxies_lin_non_tomographic", ] suff = ".png" version_string = "_".join(versions) - +# lambda_1 +print("Computing λ_1...") cv.plot_gammat_stars_around_galaxies(offset=0.025, gammax=True, wo_rand_subtr=True) cv.plot_gammat_stars_around_galaxies(offset=0.025, gammax=True, logy=True, wo_rand_subtr=True) + +# lambda_2 +print("Computing λ_2...") +cv.plot_dgammat_stars_around_galaxies(offset=0.025, gammax=True, wo_rand_subtr=True) + +# lambda_3 +print("Computing λ_3...") +cv.plot_dsize_stars_around_galaxies(offset=0.025, gammax=True, wo_rand_subtr=True) + + rename_output(output_bases, output_dir, suff, version_string, "_wo_rand_subtr") # %% +# lambda_1 +print("Computing λ_1...") cv.plot_gammat_stars_around_galaxies(offset=0.025, gammax=True) cv.plot_gammat_stars_around_galaxies(offset=0.025, gammax=True, logy=True) + +# lambda_2 +print("Computing λ_2...") +cv.plot_dgammat_stars_around_galaxies(offset=0.025, gammax=True) + +# lambda_3 +print("Computing λ_3...") +cv.plot_dsize_stars_around_galaxies(offset=0.025, gammax=True) + rename_output(output_bases, output_dir, suff, version_string, "")