From 4d85a4f3fa1980e9428636674f4cbd47deb73979 Mon Sep 17 00:00:00 2001 From: Martin Kilbinger Date: Wed, 4 Sep 2024 16:19:04 +0200 Subject: [PATCH 01/10] nb leakage autocorr (alpha term; Baceon 2003) --- .gitignore | 1 + notebooks/leakage_PSF_auto_spin.py | 10 ++++++++++ 2 files changed, 11 insertions(+) create mode 100644 notebooks/leakage_PSF_auto_spin.py diff --git a/.gitignore b/.gitignore index 5fb506c..16031b7 100644 --- a/.gitignore +++ b/.gitignore @@ -6,5 +6,6 @@ notebooks/unions_shapepipe_extended_2022_W3_v1.0.3.fits __pycache__ __pycache__/* shear_psf_leakage/__pycache__/* +shear_psf_leakage/__pycache__ */__pycache__/* build diff --git a/notebooks/leakage_PSF_auto_spin.py b/notebooks/leakage_PSF_auto_spin.py new file mode 100644 index 0000000..7d01ae4 --- /dev/null +++ b/notebooks/leakage_PSF_auto_spin.py @@ -0,0 +1,10 @@ +# + +# Plot PSF auto-correlation terms + +Xi_11 = [] +Xi_22 = [] +for ndx in range(len(theta)): + Xi_11.append(obj.Xi_pp_ufloat[ndx][0, 0]) + Xi_22.append(obj.Xi_pp_ufloat[ndx][1, 1]) + +t_1 = obj.get_alpha_ufloat(0, 0) ** 2* Xi_11 From b8d65f6701d1045d315a5acaec2f4fac8c4004c3 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 4 Sep 2024 21:30:43 +0200 Subject: [PATCH 02/10] testing Bacon xi_sys --- notebooks/leakage_autocorr_psf.py | 206 ------------------------------ notebooks/leakage_xi_sys.py | 175 +++++++++++++++++++++++++ 2 files changed, 175 insertions(+), 206 deletions(-) delete mode 100644 notebooks/leakage_autocorr_psf.py create mode 100644 notebooks/leakage_xi_sys.py diff --git a/notebooks/leakage_autocorr_psf.py b/notebooks/leakage_autocorr_psf.py deleted file mode 100644 index 6938690..0000000 --- a/notebooks/leakage_autocorr_psf.py +++ /dev/null @@ -1,206 +0,0 @@ -# --- -# jupyter: -# jupytext: -# formats: ipynb,py:light -# text_representation: -# extension: .py -# format_name: light -# format_version: '1.5' -# jupytext_version: 1.16.1 -# kernelspec: -# display_name: shear_psf_leakage -# language: python -# name: shear_psf_leakage -# --- - -# # PSF leakage: PSF auto-correlation contributions -# -# Martin Kilbinger - -# %reload_ext autoreload -# %autoreload 2 - -# + -import os -import matplotlib.pylab as plt -import numpy as np - -from astropy import units - -from uncertainties import ufloat -from uncertainties import unumpy - -from cs_util import canfar -from cs_util import plots as cs_plots - -import shear_psf_leakage.run_object as run_object -import shear_psf_leakage.run_scale as run_scale -from shear_psf_leakage import leakage -# - - -# ## Set input parameters - -params_in = {} - -# ### Paths - -# + -# Input galaxy shear catalogue -params_in["input_path_shear"] = "unions_shapepipe_extended_2022_v1.0.fits" - -# Output directory -params_in["output_dir"] = "leakage_PSF_autocorr" -# - - -# ### Job control - -# Compute (spin-preserving) PSF ellipticity leakage -params_in["PSF_leakage"] = True - -# ### Other parameters - -# + -# PSF ellipticity column names -params_in["e1_PSF_col"] = "e1_PSF" -params_in["e2_PSF_col"] = "e2_PSF" - -# Set verbose output -params_in["verbose"] = True -# - - -# ### Retrieve test catalogue from VOspace if not yet downloade - -vos_dir = "vos:cfis/XXXX/" -canfar.download( - f"{vos_dir}/{params_in['input_path_shear']}", - params_in["input_path_shear"], - verbose=params_in["verbose"], -) - -# ## Compute leakage - -# Create leakage instance -obj_object = run_object.LeakageObject() - -# Set instance parameters, copy from above -for key in params_in: - obj_object._params[key] = params_in[key] - -# Run commands as given in LeakageObject.run() - -# + -# Check parameter validity -obj_object.check_params() - -# Update parameters (here: strings to list) -obj_object.update_params() - -# Prepare output directory -obj_object.prepare_output() -# - - -# Read input catalogue -obj_object.read_data() - -# Object-by-object spin-consistent PSF leakage -mix = True -order = "lin" -obj_object.PSF_leakage(mix=mix, order=order) - -p_dp = {} -for p in obj_object.par_best_fit: - p_dp[p] = ufloat(obj_object.par_best_fit[p].value, obj_object.par_best_fit[p].stderr) -print(p_dp) - -s_ds = leakage.param_order2spin(p_dp, order, mix) -print(s_ds) - -# Get spin elements -x0 = s_ds["x0"] -x4 = s_ds["x4"] -y4 = s_ds["y4"] -y0 = s_ds["y0"] - - -# Create scale-dependent leakage instance -obj_scale = run_scale.LeakageScale() - -# Set instance parameters, copy from above -for key in params_in: - obj_scale._params[key] = params_in[key] - - -obj_scale._params["input_path_PSF"] = "unions_shapepipe_psf_2022_v1.0.2.fits" -obj_scale._params["dndz_path"] = "dndz_SP_A.txt" - -# ### Run - -# Check parameter validity -obj_scale.check_params() - -# Prepare output directory and stats file -obj_scale.prepare_output() - -# Prepare input -obj_scale.read_data() - -# #### Compute correlation function and alpha matrices -# The following command calls `treecorr` to compute auto- and cross-correlation functions. -# This can take a few minutes. - -obj_scale.compute_corr_gp_pp_alpha_matrix() - -obj_scale.alpha_matrix() - -#### For comparison: scalar alpha leakage -obj_scale.compute_corr_gp_pp_alpha() -obj_scale.do_alpha() - - -# Plot terms -theta_arcmin = obj_scale.get_theta() - -term_1 = ( - (x0 ** 2 + x4 **2 + y0 **2 + y4 ** 2) - * (obj_scale.xi_pp_m[0][0] + obj_scale.xi_pp_m[1][1]) -) -term_2 = 4 * (x0 * y4 - x4 * y0) * obj_scale.xi_pp_m[0][1] - -theta = theta_arcmin * units.arcmin -xi_theo_p, xi_theo_m = run_scale.get_theo_xi(theta, obj_scale._params["dndz_path"]) - -# + -y = [ - unumpy.nominal_values(term_1), - -unumpy.nominal_values(term_2), - xi_theo_p, -] -dy = [ - unumpy.std_devs(term_1), - unumpy.std_devs(term_2), - [np.nan] * len(xi_theo_p), -] -x = [theta_arcmin] * len(y) - -title = "PSF auto-correlation additive terms" -xlabel = r"$\theta$ [arcmin]" -ylabel = "terms" -out_path = f"{obj_scale._params['output_dir']}/auto_corr_terms.png" -labels = ["t1" , "-t2", r"$\xi_p$"] - -cs_plots.plot_data_1d( - x, - y, - dy, - title, - xlabel, - ylabel, - out_path, - labels=labels, - shift_x=True, - xlog=True, - ylog=True, -) -# - - - diff --git a/notebooks/leakage_xi_sys.py b/notebooks/leakage_xi_sys.py new file mode 100644 index 0000000..f6f767a --- /dev/null +++ b/notebooks/leakage_xi_sys.py @@ -0,0 +1,175 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:light +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.16.1 +# kernelspec: +# display_name: shear_psf_leakage +# language: python +# name: shear_psf_leakage +# --- + +# # PSF leakage: PSF auto-correlation contributions +# +# Martin Kilbinger + +# %reload_ext autoreload +# %autoreload 2 + +# + +import os +import matplotlib.pylab as plt +import numpy as np + +from astropy import units + +from uncertainties import ufloat +from uncertainties import unumpy + +from cs_util import canfar +from cs_util import plots as cs_plots +from cs_util import args + +from shear_psf_leakage import run_object +from shear_psf_leakage import run_scale +from shear_psf_leakage import leakage +# - + +# ## Compute leakage + +# Create leakage instance +obj_scale = run_scale.LeakageScale() + +# ## Set input parameters + +# Read python parameter file or get user input +params_upd = args.read_param_script("params_xi_sys.py", obj_scale._params, verbose=True) +for key in params_upd: + obj_scale._params[key] = params_upd[key] + + +# + +# Check parameter validity +obj_scale.check_params() + +# Prepare output directory +obj_scale.prepare_output() +# - + +# Read input catalogue +obj_scale.read_data() + +# #### Compute correlation function and alpha matrices +# The following command calls `treecorr` to compute auto- and cross-correlation functions. +# This can take a few minutes. + +obj_scale.compute_corr_gp_pp_alpha_matrix() + +obj_scale.alpha_matrix() + +# Spin leakage elements +alpha_0_r = ( + 0.5 * ( + obj_scale.get_alpha_ufloat(0, 0) + + obj_scale.get_alpha_ufloat(1, 1) + ) +) +alpha_0_i = ( + 0.5 * ( + -obj_scale.get_alpha_ufloat(0, 1) + + obj_scale.get_alpha_ufloat(1, 0) + ) +) +alpha_4_r = ( + 0.5 * ( + obj_scale.get_alpha_ufloat(0, 0) + - obj_scale.get_alpha_ufloat(1, 1) + ) +) +alpha_4_i = ( + 0.5 * ( + obj_scale.get_alpha_ufloat(0, 1) + + obj_scale.get_alpha_ufloat(1, 0) + ) +) + + +def get_rho_0(self, idx, jdx): + """Get Alpha Ufloat. + + Return alpha leakage matrix element as array over scales. + + Parameters + ---------- + idx : int + line index, allowed are 0 or 1 + jdx : int + column index, allowed are 0 or 1 + + Returns + ------- + numpy.ndarray + matrix element as array over scales, each entry is + of type ufloat + + """ + mat = [] + n_theta = self._params["n_theta"] + for ndx in range(n_theta): + mat.append(self.Xi_pp_ufloat[ndx][idx, jdx]) + return np.array(mat) + + +# xi_sys terms +print("xi_sys terms") +xi_sys_term_p = ( + alpha_0_r ** 2 + alpha_0_i ** 2 + alpha_4_r ** 2 + alpha_4_i ** 2 +) * (get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1)) +print("xi_sys done") + + +#### For comparison: scalar alpha leakage +obj_scale.compute_corr_gp_pp_alpha() +obj_scale.do_alpha() + + +# Plot terms +theta_arcmin = obj_scale.get_theta() + +# + +y = [ + unumpy.nominal_values(xi_sys_term_p), + #xi_theo_p, +] +dy = [ + unumpy.std_devs(xi_sys_term_p), + #[np.nan] * len(xi_theo_p), +] +x = [theta_arcmin] * len(y) + +title = r"Bacon et al. (2003) $\xi_{sys}$" +xlabel = r"$\theta$ [arcmin]" +ylabel = "terms" +out_path = f"{obj_scale._params['output_dir']}/xi_sys.png" +labels = ["$t_+$"] + +cs_plots.plot_data_1d( + x, + y, + dy, + title, + xlabel, + ylabel, + out_path, + labels=labels, + shift_x=True, + xlog=True, + ylog=True, +) +# - + + From 66fd647d73db777e0d7d5d465e81fa055efd4d10 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Fri, 18 Oct 2024 14:49:57 +0200 Subject: [PATCH 03/10] xi_sys run and plot script updated to new eqs. --- notebooks/leakage_xi_sys.py | 60 +++++++++++++++++++++++++++---------- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/notebooks/leakage_xi_sys.py b/notebooks/leakage_xi_sys.py index f6f767a..569a7e7 100644 --- a/notebooks/leakage_xi_sys.py +++ b/notebooks/leakage_xi_sys.py @@ -6,14 +6,14 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.1 +# jupytext_version: 1.15.1 # kernelspec: -# display_name: shear_psf_leakage +# display_name: sp_validation # language: python -# name: shear_psf_leakage +# name: python3 # --- -# # PSF leakage: PSF auto-correlation contributions +# # PSF contamination: PSF leakage diagnostic # # Martin Kilbinger @@ -99,9 +99,9 @@ def get_rho_0(self, idx, jdx): - """Get Alpha Ufloat. + """Get Rho 0. - Return alpha leakage matrix element as array over scales. + Return (i, j) matrix element of rho_0. Parameters ---------- @@ -127,14 +127,25 @@ def get_rho_0(self, idx, jdx): # xi_sys terms print("xi_sys terms") xi_sys_term_p = ( - alpha_0_r ** 2 + alpha_0_i ** 2 + alpha_4_r ** 2 + alpha_4_i ** 2 -) * (get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1)) + 2 * ( + alpha_0_r ** 2 + alpha_0_i ** 2 + alpha_4_r ** 2 + alpha_4_i ** 2 + ) * (get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1)) +) +xi_sys_term_m = 2 * alpha_0_i * alpha_4_i * ( + (get_rho_0(obj_scale, 0, 0) - get_rho_0(obj_scale, 1, 1)) +) +xi_sys_term_mixed = ( + 4 * ( + alpha_0_r * alpha_4_i - alpha_4_r * alpha_0_i + ) * get_rho_0(obj_scale, 0, 1) +) print("xi_sys done") - -#### For comparison: scalar alpha leakage + +#### For comparison: scalar xi sys obj_scale.compute_corr_gp_pp_alpha() obj_scale.do_alpha() +obj_scale.compute_xi_sys() # Plot terms @@ -143,19 +154,23 @@ def get_rho_0(self, idx, jdx): # + y = [ unumpy.nominal_values(xi_sys_term_p), - #xi_theo_p, + unumpy.nominal_values(xi_sys_term_m), + unumpy.nominal_values(xi_sys_term_mixed), + obj_scale.C_sys_p, ] dy = [ unumpy.std_devs(xi_sys_term_p), - #[np.nan] * len(xi_theo_p), + unumpy.std_devs(xi_sys_term_m), + unumpy.std_devs(xi_sys_term_mixed), + obj_scale.C_sys_std_p, ] x = [theta_arcmin] * len(y) title = r"Bacon et al. (2003) $\xi_{sys}$" xlabel = r"$\theta$ [arcmin]" ylabel = "terms" -out_path = f"{obj_scale._params['output_dir']}/xi_sys.png" -labels = ["$t_+$"] +out_path = f"{obj_scale._params['output_dir']}/xi_sys_terms.png" +labels = ["$t_+$", "$t_-$", r"$t_{\rm mixed}$", "scalar"] cs_plots.plot_data_1d( x, @@ -168,8 +183,23 @@ def get_rho_0(self, idx, jdx): labels=labels, shift_x=True, xlog=True, - ylog=True, + ylog=False, ) # - +obj_scale.C_sys_p + +print(unumpy.nominal_values( + obj_scale.alpha_leak ** 2 * ( + get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1) + ) +)) +print() +print(obj_scale.alpha_leak ** 2 * obj_scale.r_corr_pp.xip) +print() +print(obj_scale.C_sys_p) + +obj_scale.alpha_leak + +alpha_4_i From 0cb4093c51aa3c74b2ead5206b436c0b7c28ab56 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Fri, 18 Oct 2024 14:50:13 +0200 Subject: [PATCH 04/10] x-limit for log plot changed (fixed?) --- shear_psf_leakage/run_scale.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/shear_psf_leakage/run_scale.py b/shear_psf_leakage/run_scale.py index d5783e6..006360a 100644 --- a/shear_psf_leakage/run_scale.py +++ b/shear_psf_leakage/run_scale.py @@ -771,15 +771,16 @@ def plot_alpha_leakage(self, close_fig=True, xlinlog="log"): """ # Set some x-values and limits for plot + x0 = self._params["theta_min_amin"] if xlinlog == "log": - x0 = self._params["theta_min_amin"] + x0 = x0 * 0.9 x_affine = np.geomspace(x0, self._params["theta_max_amin"]) xlim = [x0, self._params["theta_max_amin"]] xlog = True else: - x0 = -10 - x_affine = np.linspace(0, self._params["theta_max_amin"]) - xlim = [x0 * 2, self._params["theta_max_amin"]] + x0 = x0 ** 0.9 + x_affine = np.linspace(x0, self._params["theta_max_amin"]) + xlim = [x0, self._params["theta_max_amin"]] xlog = False # alpha(theta) data points From 84ca1c81877121075ffbe353c4438be648bd7f64 Mon Sep 17 00:00:00 2001 From: Martin Kilbinger Date: Sat, 26 Oct 2024 10:23:07 +0200 Subject: [PATCH 05/10] .gitignore --- citation.cff | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 citation.cff diff --git a/citation.cff b/citation.cff new file mode 100644 index 0000000..8143e34 --- /dev/null +++ b/citation.cff @@ -0,0 +1,16 @@ +cff-version: 1.2.0 +message: "If you use this software, please cite it as below." +authors: +- family-names: "Kilbinger" + given-names: "Martin" + orcid: "https://orcid.org/0000-0001-9513-7138" +- family-names: "Guerrini" + given-names: "Sacha" + orcid: "https://orcid.org/0009-0004-3655-4870 +- family-names: "Bonini" + given-names: "Clara" +title: "Shear PSF leakage" +version: 0.1.6 +doi: 00.0000/zenodo.0000 +date-released: 2024-12-12 +url: "https://github.com/CosmoStat/shear_psf_leakage" From 0317b11f3aff2fe0e066119e6e0633fbd639cd7a Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Tue, 29 Oct 2024 11:47:37 +0100 Subject: [PATCH 06/10] xi sys plots --- notebooks/leakage_xi_sys.py | 158 ++++++++++++++++++++++++++++-------- 1 file changed, 124 insertions(+), 34 deletions(-) diff --git a/notebooks/leakage_xi_sys.py b/notebooks/leakage_xi_sys.py index 569a7e7..132c365 100644 --- a/notebooks/leakage_xi_sys.py +++ b/notebooks/leakage_xi_sys.py @@ -13,12 +13,53 @@ # name: python3 # --- -# # PSF contamination: PSF leakage diagnostic +# # PSF contamination: Star-galaxy correlation function $\xi_{\rm{sys}}$ # # Martin Kilbinger +# The PSF contamination diagnostic function $\xi_{\rm{sys}} was introduced in Bacon et al. (2003) in the scalar PSF approximation. +# It is defined as +# \begin{equation} +# \xi_{\textrm{sys}, +}(\theta) +# = \frac{\left( \xi_+^\textrm{e, p} \right)^2(\theta)} +# {\xi_+^\textrm{p, p}(\theta)} +# = \frac{\tau_{+, 0}^2(\theta)}{\rho_{0, +}(\theta)} +# = \alpha^2(\theta) \rho_{0, +}(\theta). +# \end{equation} +# +# The generalisation of this function to spin-consistent notation is the matrix +# \begin{equation} +# \xi_{\textrm{sys}, ij}(\theta) +# = \sum_{k,l=1}^2 \alpha_{ik} \alpha_{jl} \, \rho_{0, kl}(\theta). +# \end{equation} +# The ‘+’ components can be computed as +# \begin{align} +# \xi_{\textrm{sys}, +}(\theta) +# = & \sum_{k,l=1}^2 \sum_{i=1}^2 \alpha_{ik} \alpha_{jl} \rho_{0, kl}(\theta) +# \nonumber \\ +# = & \, +# \left( +# \alpha_{11}^2 + \alpha_{22}^2 +# \right) +# \rho_{0, +}(\theta) +# + \alpha_{21}^2 \, \rho_{0, 11}(\theta) + \alpha_{12}^2 \, \rho_{0, 22}(\theta) +# \nonumber \\= & \, +# 2 \left[ \left(\alpha_0^\Re \right)^2 + \left(\alpha_4^\Re\right)^2 \right] +# \rho_{0, +}(\theta) +# \nonumber \\ +# & +# + \left( \alpha^\Im_0 + \alpha^\Im_4 \right)^2 +# \left[ \rho_{0, 11}(\theta) - \rho_{0, 22}(\theta) \right] +# \nonumber \\ +# & +# + 4 \left( \alpha^\Re_0 \alpha^\Im_4 + \alpha^\Re_4 \alpha^\Im_0 \right) +# \rho_{0, 12}(\theta) +# , +# \end{align} + # %reload_ext autoreload # %autoreload 2 +# %matplotlib inline # + import os @@ -39,13 +80,11 @@ from shear_psf_leakage import leakage # - -# ## Compute leakage +# ## Set up -# Create leakage instance +# Create instance of scale-dependent leakage object obj_scale = run_scale.LeakageScale() -# ## Set input parameters - # Read python parameter file or get user input params_upd = args.read_param_script("params_xi_sys.py", obj_scale._params, verbose=True) for key in params_upd: @@ -60,7 +99,9 @@ obj_scale.prepare_output() # - -# Read input catalogue +# ### Read input catalogues + +# Read input galaxy and star catalogue obj_scale.read_data() # #### Compute correlation function and alpha matrices @@ -71,7 +112,9 @@ obj_scale.alpha_matrix() -# Spin leakage elements +# + +# Compute spin coefficients from matrix elements + alpha_0_r = ( 0.5 * ( obj_scale.get_alpha_ufloat(0, 0) @@ -97,6 +140,10 @@ ) ) +print(alpha_0_r[0], alpha_0_i[0], alpha_4_r[0], alpha_4_i[0]) + + +# - def get_rho_0(self, idx, jdx): """Get Rho 0. @@ -116,7 +163,7 @@ def get_rho_0(self, idx, jdx): matrix element as array over scales, each entry is of type ufloat - """ + """ mat = [] n_theta = self._params["n_theta"] for ndx in range(n_theta): @@ -124,45 +171,72 @@ def get_rho_0(self, idx, jdx): return np.array(mat) +# + # xi_sys terms -print("xi_sys terms") xi_sys_term_p = ( - 2 * ( - alpha_0_r ** 2 + alpha_0_i ** 2 + alpha_4_r ** 2 + alpha_4_i ** 2 - ) * (get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1)) -) -xi_sys_term_m = 2 * alpha_0_i * alpha_4_i * ( - (get_rho_0(obj_scale, 0, 0) - get_rho_0(obj_scale, 1, 1)) + 2 * (alpha_0_r ** 2 + alpha_4_r ** 2) + * (get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1)) ) +xi_sys_term_m = (alpha_0_i + alpha_4_i) * (alpha_0_i + alpha_4_i) * (get_rho_0(obj_scale, 0, 0) - get_rho_0(obj_scale, 1, 1)) xi_sys_term_mixed = ( 4 * ( alpha_0_r * alpha_4_i - alpha_4_r * alpha_0_i ) * get_rho_0(obj_scale, 0, 1) ) -print("xi_sys done") + +xi_sys_total = xi_sys_term_p + xi_sys_term_m + xi_sys_term_mixed +# - -#### For comparison: scalar xi sys +#### For comparison: scalar xi sys. obj_scale.compute_corr_gp_pp_alpha() obj_scale.do_alpha() obj_scale.compute_xi_sys() -# Plot terms +# ## Plotting + +# angular scales in arcmin theta_arcmin = obj_scale.get_theta() # + +# Mean ellipticities for centered correlation functions + +e1_gal = obj_scale.dat_shear["e1"] +e2_gal = obj_scale.dat_shear["e2"] +weights_gal = obj_scale.dat_shear["w"] + +complex_gal = ( + np.average(e1_gal, weights=weights_gal) + + np.average(e2_gal, weights=weights_gal) * 1j +) + +e1_star = obj_scale.dat_PSF["E1_STAR_HSM"] +e2_star = obj_scale.dat_PSF["E2_STAR_HSM"] +complex_psf = np.mean(e1_star) + np.mean(e2_star) * 1j + +# + +# At the moment we compute the scalar xi_sys here using centered correlation functions. +# TODO: Implement consistent centering in shear_psf_leakage classes. + +# xi_sys = tau_0^2/rho_0 +xi_sys_scalar = (obj_scale.r_corr_gp.xip - np.real(np.conj(complex_gal) * complex_psf)) ** 2 / (obj_scale.r_corr_pp.xip - np.abs((complex_psf) ** 2)) +std_xi_sys_scalar = obj_scale.C_sys_std_p + y = [ unumpy.nominal_values(xi_sys_term_p), unumpy.nominal_values(xi_sys_term_m), unumpy.nominal_values(xi_sys_term_mixed), - obj_scale.C_sys_p, + unumpy.nominal_values(xi_sys_total), + xi_sys_scalar, + ] dy = [ unumpy.std_devs(xi_sys_term_p), unumpy.std_devs(xi_sys_term_m), unumpy.std_devs(xi_sys_term_mixed), - obj_scale.C_sys_std_p, + unumpy.std_devs(xi_sys_total), + std_xi_sys_scalar, ] x = [theta_arcmin] * len(y) @@ -170,7 +244,9 @@ def get_rho_0(self, idx, jdx): xlabel = r"$\theta$ [arcmin]" ylabel = "terms" out_path = f"{obj_scale._params['output_dir']}/xi_sys_terms.png" -labels = ["$t_+$", "$t_-$", r"$t_{\rm mixed}$", "scalar"] +labels = ["$t_+$", "$t_-$", r"$t_{\rm mixed}$", r"$\sum t_i$", "scalar"] + +ylim = [-1e-6, 2e-6] cs_plots.plot_data_1d( x, @@ -184,22 +260,36 @@ def get_rho_0(self, idx, jdx): shift_x=True, xlog=True, ylog=False, + ylim=ylim, ) -# - -obj_scale.C_sys_p +# + +# Test: compare different estimates of xi_sys. Some use uncentered, some centered correlation functions. -print(unumpy.nominal_values( - obj_scale.alpha_leak ** 2 * ( - get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1) - ) -)) -print() -print(obj_scale.alpha_leak ** 2 * obj_scale.r_corr_pp.xip) -print() -print(obj_scale.C_sys_p) +plt.clf() +plt.loglog(theta_arcmin, obj_scale.C_sys_p, "p:", label="scalar uncentered") + +# spin t_+ +plt.loglog( + theta_arcmin, + unumpy.nominal_values(xi_sys_term_p), + "d-", + label="spin $t_+$ centered", +) -obj_scale.alpha_leak +# tau^2 / rho +plt.loglog(theta_arcmin, obj_scale.r_corr_gp.xip ** 2 / obj_scale.r_corr_pp.xip, "v:", label="scalar centered") +plt.loglog( + theta_arcmin, + (obj_scale.r_corr_gp.xip - np.real(np.conj(complex_gal) * complex_psf)) ** 2 / (obj_scale.r_corr_pp.xip - np.abs((complex_psf) ** 2)), + "v-", + label="3c" +) + +plt.legend() +plt.show() +plt.savefig(f"{obj_scale._params['output_dir']}/xi_sys_test.png") + +# - -alpha_4_i From 275e9106dc7aee630a83a3b5203008f52352382b Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Tue, 29 Oct 2024 15:59:21 +0100 Subject: [PATCH 07/10] terms in agreement --- notebooks/leakage_xi_sys.py | 187 +++++++++++++++++++++++++++++++++--- 1 file changed, 173 insertions(+), 14 deletions(-) diff --git a/notebooks/leakage_xi_sys.py b/notebooks/leakage_xi_sys.py index 132c365..dc5691b 100644 --- a/notebooks/leakage_xi_sys.py +++ b/notebooks/leakage_xi_sys.py @@ -35,24 +35,42 @@ # The ‘+’ components can be computed as # \begin{align} # \xi_{\textrm{sys}, +}(\theta) -# = & \sum_{k,l=1}^2 \sum_{i=1}^2 \alpha_{ik} \alpha_{jl} \rho_{0, kl}(\theta) +# = & \sum_{k,l=1}^2 \sum_{i=1}^2 \alpha_{ik} \alpha_{il} \rho_{0, kl}(\theta) # \nonumber \\ # = & \, # \left( -# \alpha_{11}^2 + \alpha_{22}^2 -# \right) -# \rho_{0, +}(\theta) -# + \alpha_{21}^2 \, \rho_{0, 11}(\theta) + \alpha_{12}^2 \, \rho_{0, 22}(\theta) -# \nonumber \\= & \, -# 2 \left[ \left(\alpha_0^\Re \right)^2 + \left(\alpha_4^\Re\right)^2 \right] +# \alpha_{11}^2 + \alpha_{21}^2 +# \right) \rho_{0, 11}(\theta) +# + \left( +# \alpha_{12}^2 + \alpha_{22} +# \right)^2 \rho_{0, 22}(\theta) +# \nonumber \\ +# & \, +# + 2 \left( +# \alpha_{11} \alpha_{12} + \alpha_{21} \alpha_{22} +# \right) \rho_{0, 12}(\theta) +# \nonumber \\ +# = & \, +# \left| \alpha_0 + \alpha_4 \right|^2 \rho_{0, 11}(\theta) +# + \left| \alpha_0 - \alpha_4 \right|^2 \rho_{0, 22}(\theta) +# \nonumber \\ +# & \, +# + 4 \left( \alpha_0^\Re \alpha_4^\Im + \alpha_4^\Re \alpha_0^\Im \right) +# \rho_{0, 12}(\theta) +# \nonumber \\ +# = & \, +# 2 \left[ +# \left(\alpha_0^\Re \right)^2 + \left(\alpha_4^\Im\right)^2 +# + \left(\alpha_0^\Re \right)^2 + \left(\alpha_4^\Im\right)^2 +# \right] # \rho_{0, +}(\theta) # \nonumber \\ # & -# + \left( \alpha^\Im_0 + \alpha^\Im_4 \right)^2 +# + 2 \left( \alpha_0^\Re \alpha_4^\Re + \alpha_0^\Im \alpha_4^\Im \right) # \left[ \rho_{0, 11}(\theta) - \rho_{0, 22}(\theta) \right] # \nonumber \\ # & -# + 4 \left( \alpha^\Re_0 \alpha^\Im_4 + \alpha^\Re_4 \alpha^\Im_0 \right) +# + 4 \left( \alpha_0^\Re \alpha_4^\Im + \alpha_4^\Re \alpha_0^\Im \right) # \rho_{0, 12}(\theta) # , # \end{align} @@ -173,18 +191,25 @@ def get_rho_0(self, idx, jdx): # + # xi_sys terms + + +# Split into +, '-', and mixed parts + xi_sys_term_p = ( - 2 * (alpha_0_r ** 2 + alpha_4_r ** 2) + (alpha_0_r ** 2 + alpha_0_i ** 2 + alpha_4_r ** 2 + alpha_4_i ** 2) * (get_rho_0(obj_scale, 0, 0) + get_rho_0(obj_scale, 1, 1)) ) -xi_sys_term_m = (alpha_0_i + alpha_4_i) * (alpha_0_i + alpha_4_i) * (get_rho_0(obj_scale, 0, 0) - get_rho_0(obj_scale, 1, 1)) +xi_sys_term_m = ( + 2 * (alpha_0_r * alpha_4_r + alpha_0_i * alpha_4_i) + * (get_rho_0(obj_scale, 0, 0) - get_rho_0(obj_scale, 1, 1)) +) xi_sys_term_mixed = ( 4 * ( alpha_0_r * alpha_4_i - alpha_4_r * alpha_0_i ) * get_rho_0(obj_scale, 0, 1) ) -xi_sys_total = xi_sys_term_p + xi_sys_term_m + xi_sys_term_mixed +xi_sys_tot = xi_sys_term_p + xi_sys_term_m + xi_sys_term_mixed # - @@ -227,7 +252,7 @@ def get_rho_0(self, idx, jdx): unumpy.nominal_values(xi_sys_term_p), unumpy.nominal_values(xi_sys_term_m), unumpy.nominal_values(xi_sys_term_mixed), - unumpy.nominal_values(xi_sys_total), + unumpy.nominal_values(xi_sys_tot), xi_sys_scalar, ] @@ -235,7 +260,7 @@ def get_rho_0(self, idx, jdx): unumpy.std_devs(xi_sys_term_p), unumpy.std_devs(xi_sys_term_m), unumpy.std_devs(xi_sys_term_mixed), - unumpy.std_devs(xi_sys_total), + unumpy.std_devs(xi_sys_tot), std_xi_sys_scalar, ] x = [theta_arcmin] * len(y) @@ -245,6 +270,7 @@ def get_rho_0(self, idx, jdx): ylabel = "terms" out_path = f"{obj_scale._params['output_dir']}/xi_sys_terms.png" labels = ["$t_+$", "$t_-$", r"$t_{\rm mixed}$", r"$\sum t_i$", "scalar"] +markers = ["o", "d", "^", "x", "s"] ylim = [-1e-6, 2e-6] @@ -261,7 +287,140 @@ def get_rho_0(self, idx, jdx): xlog=True, ylog=False, ylim=ylim, + markers=markers, +) +# + +# Consistency test. Different computations of xi_sys terms + +# Using a_ij matrix +xi_sys_term_m_11 = ( + (obj_scale.get_alpha_ufloat(0, 0) ** 2 + obj_scale.get_alpha_ufloat(1, 0) ** 2) + * get_rho_0(obj_scale, 0, 0) ) +xi_sys_term_m_22 = ( + (obj_scale.get_alpha_ufloat(0, 1) ** 2 + obj_scale.get_alpha_ufloat(1, 1) ** 2) + * get_rho_0(obj_scale, 1, 1) +) +xi_sys_term_m_12 = ( + 2 * (obj_scale.get_alpha_ufloat(0, 0) * obj_scale.get_alpha_ufloat(0, 1) + + obj_scale.get_alpha_ufloat(1, 0) * obj_scale.get_alpha_ufloat(1, 1)) + * get_rho_0(obj_scale, 0, 1) +) + +# Sum of the three terms +xi_sys_m_tot = xi_sys_term_m_11 + xi_sys_term_m_22 + xi_sys_term_m_12 + +# Matrix sum +xi_sys_m2_tot = 0 +for idx in (0, 1): + for kdx in (0, 1): + for ldx in (0, 1): + xi_sys_m2_tot += ( + obj_scale.get_alpha_ufloat(idx, kdx) + * obj_scale.get_alpha_ufloat(idx, ldx) + * get_rho_0(obj_scale, kdx, ldx) + ) + +y = [ + unumpy.nominal_values(xi_sys_term_m_11), + unumpy.nominal_values(xi_sys_term_m_22), + unumpy.nominal_values(xi_sys_term_m_12), + unumpy.nominal_values(xi_sys_m_tot), + unumpy.nominal_values(xi_sys_m2_tot), + unumpy.nominal_values(xi_sys_tot), + +] +dy = [ + unumpy.std_devs(xi_sys_term_m_11), + unumpy.std_devs(xi_sys_term_m_22), + unumpy.std_devs(xi_sys_term_m_12), + unumpy.std_devs(xi_sys_m_tot), + unumpy.std_devs(xi_sys_m2_tot), + unumpy.std_devs(xi_sys_tot), +] +x = [theta_arcmin] * len(y) + +title = r"Bacon et al. (2003) $\xi_{sys}$" +xlabel = r"$\theta$ [arcmin]" +ylabel = "terms with matrix coeffs" +out_path = f"{obj_scale._params['output_dir']}/xi_sys_terms_m.png" +labels = ["$t_{11}$", "$t_{22}$", r"$t_{12}$", r"$\sum t_i$ (matrix terms)", r"$\sum t_i$ (matrix terms 2)", "$\sum t_i$ (spin terms)"] +markers = ["o", "d", "^", "x", "h", "s"] +linestyles = [":"] * 3 +linestyles.extend(["-", "-.", "--"]) + +ylim = [-1e-6, 2e-6] + +cs_plots.plot_data_1d( + x, + y, + dy, + title, + xlabel, + ylabel, + out_path, + labels=labels, + shift_x=False, + xlog=True, + ylog=False, + ylim=ylim, + markers=markers, + linestyles=linestyles, +) + + + +# + +# Plot difference, should be consistent with zero + +d1 = xi_sys_m_tot - xi_sys_tot +d2 = xi_sys_m2_tot - xi_sys_tot +d3 = xi_sys_m_tot - xi_sys_m2_tot + +d4 = xi_sys_term_m_11 + xi_sys_term_m_22 - (xi_sys_term_p + xi_sys_term_m) +d5 = xi_sys_term_m_12 - xi_sys_term_mixed + +y = [ + unumpy.nominal_values(d1), + unumpy.nominal_values(d2), + unumpy.nominal_values(d3), + unumpy.nominal_values(d4), + unumpy.nominal_values(d5), +] +dy = [ + unumpy.std_devs(d1), + unumpy.std_devs(d2), + unumpy.std_devs(d3), + #unumpy.std_devs(d4), + #unumpy.std_devs(d5), +] +x = [theta_arcmin] * len(y) + +title = r"Bacon et al. (2003) $\xi_{sys}$" +xlabel = r"$\theta$ [arcmin]" +ylabel = "difference" +out_path = f"{obj_scale._params['output_dir']}/xi_sys_diff.png" +labels = [r"(matrix - spin) tot", r"(matrix 2 - spin) tot", "(matrix - matrix 2) tot", "11 22", "mixed"] +markers = ["o", "d", "s", "x", "v"] +linestyles = ["-", "--", ":", ":", ":"] + + +cs_plots.plot_data_1d( + x[:3], + y[:3], + dy[:3], + title, + xlabel, + ylabel, + out_path, + labels=labels, + shift_x=True, + xlog=True, + ylog=False, + markers=markers, + linestyles=linestyles, +) + # + # Test: compare different estimates of xi_sys. Some use uncentered, some centered correlation functions. From 2b0005999193c977dc511078bf93df61655d9be2 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Tue, 29 Oct 2024 17:17:00 +0100 Subject: [PATCH 08/10] deleted obsolete notebook script --- notebooks/leakage_PSF_auto_spin.py | 10 ---------- 1 file changed, 10 deletions(-) delete mode 100644 notebooks/leakage_PSF_auto_spin.py diff --git a/notebooks/leakage_PSF_auto_spin.py b/notebooks/leakage_PSF_auto_spin.py deleted file mode 100644 index 7d01ae4..0000000 --- a/notebooks/leakage_PSF_auto_spin.py +++ /dev/null @@ -1,10 +0,0 @@ -# + -# Plot PSF auto-correlation terms - -Xi_11 = [] -Xi_22 = [] -for ndx in range(len(theta)): - Xi_11.append(obj.Xi_pp_ufloat[ndx][0, 0]) - Xi_22.append(obj.Xi_pp_ufloat[ndx][1, 1]) - -t_1 = obj.get_alpha_ufloat(0, 0) ** 2* Xi_11 From 8f0b9bf050210f1ab592aef4790106aa958c077a Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 30 Oct 2024 11:21:32 +0100 Subject: [PATCH 09/10] cleand up plot alpha --- shear_psf_leakage/run_scale.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/shear_psf_leakage/run_scale.py b/shear_psf_leakage/run_scale.py index 006360a..f78b087 100644 --- a/shear_psf_leakage/run_scale.py +++ b/shear_psf_leakage/run_scale.py @@ -771,14 +771,15 @@ def plot_alpha_leakage(self, close_fig=True, xlinlog="log"): """ # Set some x-values and limits for plot + factor = 0.9 x0 = self._params["theta_min_amin"] if xlinlog == "log": - x0 = x0 * 0.9 + x0 = x0 * factor x_affine = np.geomspace(x0, self._params["theta_max_amin"]) xlim = [x0, self._params["theta_max_amin"]] xlog = True else: - x0 = x0 ** 0.9 + x0 = x0 ** factor x_affine = np.linspace(x0, self._params["theta_max_amin"]) xlim = [x0, self._params["theta_max_amin"]] xlog = False From 877ba610cc5bdb5763dcfa536caf3ba401856f83 Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Wed, 30 Oct 2024 17:57:51 +0100 Subject: [PATCH 10/10] leakage spin coeff added to class --- notebooks/leakage_scale_spin.py | 228 ++++++++++++++++---------------- notebooks/leakage_xi_sys.py | 15 ++- shear_psf_leakage/run_scale.py | 32 +++++ 3 files changed, 156 insertions(+), 119 deletions(-) diff --git a/notebooks/leakage_scale_spin.py b/notebooks/leakage_scale_spin.py index ed20567..456354c 100644 --- a/notebooks/leakage_scale_spin.py +++ b/notebooks/leakage_scale_spin.py @@ -7,23 +7,22 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.16.1 +# jupytext_version: 1.15.1 # kernelspec: -# display_name: shear-psf-leakage-cHl61fDi-py3.9 +# display_name: sp_validation # language: python # name: python3 # --- -# # Scale-dependent PSF leakage +# # PSF contamination: Scale-dependent PSF leakage # -# ## Spin-consistent calculations -# -# Computes and plots scale-dependent spin-consistent PSF leakage functions. +# ## Spin-consistent formalism # # Martin Kilbinger # %reload_ext autoreload # %autoreload 2 +# %matplotlib inline # + import os @@ -35,45 +34,53 @@ from cs_util import plots as cs_plots from cs_util import args -import shear_psf_leakage.run_scale as run +from shear_psf_leakage import run_scale from shear_psf_leakage.leakage import * # - -# ## Compute leakage +# ## Set up -# Create leakage instance -obj = run.LeakageScale() +# Create instance of scale-dependent leakage object +obj_scale = run_scale.LeakageScale() # ### Set input parameters # Read python parameter file or get user input -params_upd = args.read_param_script("params_leakage_scale.py", obj._params, verbose=True) +params_upd = args.read_param_script( + "params_leakage_scale.py", + obj_scale._params, + verbose=True +) for key in params_upd: - obj._params[key] = params_upd[key] + obj_scale._params[key] = params_upd[key] # ### Run +# + # Check parameter validity -obj.check_params() +obj_scale.check_params() # Prepare output directory and stats file -obj.prepare_output() +obj_scale.prepare_output() +# - + +# ### Read input catalogues -# Prepare input -obj.read_data() +# Read input galaxy and star catalogue +obj_scale.read_data() # #### Compute correlation function and alpha matrices # The following command calls `treecorr` to compute auto- and cross-correlation functions. # This can take a few minutes. -obj.compute_corr_gp_pp_alpha_matrix() +obj_scale.compute_corr_gp_pp_alpha_matrix() -obj.alpha_matrix() +obj_scale.alpha_matrix() #### For comparison: scalar alpha leakage -obj.compute_corr_gp_pp_alpha() -obj.do_alpha() +obj_scale.compute_corr_gp_pp_alpha() +obj_scale.do_alpha() # #### PSF auto-correlation function correlation matrix # @@ -85,25 +92,26 @@ # $$ # Check symmetry of PSF auto-correlation matrix -diff = obj.Xi_pp_m[0][1] - obj.Xi_pp_m[1][0] +diff = obj_scale.Xi_pp_m[0][1] - obj_scale.Xi_pp_m[1][0] print( "Is r symmetrical? max abs (rel) diff =" - +f" {max(np.abs(diff)):.3e} ({max(np.abs(diff / obj.Xi_pp_m[0][1])):.3e})", + + f" {max(np.abs(diff)):.3e}" + + f" ({max(np.abs(diff / obj_scale.Xi_pp_m[0][1])):.3e})", ) # + # Plot r and ratios of r -theta = obj.get_theta() +theta = obj_scale.get_theta() -# Exact calculation +# Exact: Using centered correlation functions r = [] r_ratio_1 = [] r_ratio_2 = [] for ndx in range(len(theta)): my_r = ( - obj.Xi_pp_ufloat[ndx][0, 1] ** 2 - / (obj.Xi_pp_ufloat[ndx][0, 0] * obj.Xi_pp_ufloat[ndx][1, 1]) + obj_scale.Xi_pp_ufloat[ndx][0, 1] ** 2 + / (obj_scale.Xi_pp_ufloat[ndx][0, 0] * obj_scale.Xi_pp_ufloat[ndx][1, 1]) ) r.append(my_r) r_ratio_1.append(1 / (1 - my_r)) @@ -111,8 +119,8 @@ print("min max mean r = ", np.min(r), np.max(r), np.mean(r)) -# Approximate -r_fast = obj.xi_pp_m[0][1] ** 2 / (obj.xi_pp_m[0][0] * obj.xi_pp_m[1][1]) +# Approximate: Using uncentered correlation functions +r_fast = obj_scale.xi_pp_m[0][1] ** 2 / (obj_scale.xi_pp_m[0][0] * obj_scale.xi_pp_m[1][1]) n = 6 theta_arr = [theta] * n @@ -140,14 +148,14 @@ xlabel = r"$\theta$ [arcmin]" ylabel = r"functions of $r(\theta)$" -fac = 1.1 +fac = 0.9 xlim = [ - obj._params["theta_min_amin"] / fac, - obj._params["theta_max_amin"] * fac + obj_scale._params["theta_min_amin"] * fac, + obj_scale._params["theta_max_amin"], ] ylim = (-0.5, 2) -out_path = f"{obj._params['output_dir']}/r.png" +out_path = f"{obj_scale._params['output_dir']}/r.png" title = "" @@ -173,19 +181,17 @@ # + # Plot alpha matrix elements -theta = obj.get_theta() +ylim = obj_scale._params["leakage_alpha_ylim"] -xlim = [obj._params["theta_min_amin"] / fac, obj._params["theta_max_amin"] * fac] -ylim = obj._params["leakage_alpha_ylim"] - -theta_arr = [theta] * 4 +n = 4 +theta_arr = [theta] * n alpha = [] yerr = [] labels = [] for idx in (0, 1): for jdx in (0, 1): - alpha_ufloat = obj.get_alpha_ufloat(idx, jdx) + alpha_ufloat = obj_scale.get_alpha_ufloat(idx, jdx) alpha.append(unumpy.nominal_values(alpha_ufloat)) yerr.append(unumpy.std_devs(alpha_ufloat)) labels.append(rf"$\alpha_{{{idx+1}{jdx+1}}}$") @@ -197,7 +203,7 @@ xlabel = r"$\theta$ [arcmin]" ylabel = r"$\alpha_{ij}(\theta)$" title = "" -out_path = f"{obj._params['output_dir']}/alpha_ij.png" +out_path = f"{obj_scale._params['output_dir']}/alpha_ij.png" cs_plots.plot_data_1d( theta_arr, @@ -219,49 +225,39 @@ ) # + -# Plot spin elements: Spin-0 and spin-4 - -x0 = 0.5 * ( - obj.get_alpha_ufloat(0, 0) - + obj.get_alpha_ufloat(1, 1) -) -x4 = 0.5 * ( - obj.get_alpha_ufloat(0, 0) - - obj.get_alpha_ufloat(1, 1) -) -y4 = 0.5 * (obj.get_alpha_ufloat(0, 1) + obj.get_alpha_ufloat(1, 0)) -y0 = 0.5 * (obj.get_alpha_ufloat(1, 0) - obj.get_alpha_ufloat(0, 1)) +# Plot spin coefficients +obj_scale.compute_alpha_spin_coeffs() -ylim = obj._params["leakage_alpha_ylim"] +n = 4 +theta_arr = [theta] * n y_arr = [ - unumpy.nominal_values(x0), - unumpy.nominal_values(x4), - unumpy.nominal_values(y4), - unumpy.nominal_values(y0), + unumpy.nominal_values(obj_scale._alpha_0_r), + unumpy.nominal_values(obj_scale._alpha_0_i), + unumpy.nominal_values(obj_scale._alpha_4_r), + unumpy.nominal_values(obj_scale._alpha_4_i), ] dy_arr = [ - unumpy.std_devs(x0), - unumpy.std_devs(x4), - unumpy.std_devs(y4), - unumpy.std_devs(y0), + unumpy.std_devs(obj_scale._alpha_0_r), + unumpy.std_devs(obj_scale._alpha_0_i), + unumpy.std_devs(obj_scale._alpha_4_r), + unumpy.std_devs(obj_scale._alpha_4_i), ] -theta_arr = [obj.get_theta()] * len(y_arr) labels = [ - r"$\alpha^\Re_0 = (\alpha_{11} + \alpha_{22})/2$", - r"$\alpha^\Re_4 = (\alpha_{11} - \alpha_{22})/2$", - r"$\alpha^\Im_4 = (\alpha_{12} + \alpha_{21})/2$", - r"$\alpha^\Im_0 = (-\alpha_{12} + \alpha_{21})/2$" + r"$\alpha^\Re_0$", + r"$\alpha^\Im_0$", + r"$\alpha^\Re_r$", + r"$\alpha^\Im_r$" ] colors = ["blue", "orange", "green", "magenta"] -markers = ["o", "s", "h", "^"] +markers = ["o", "s", "^", "v"] linestyles = ["-"] * 4 xlabel = r"$\theta$ [arcmin]" ylabel = r"Components of leakage matrix" title = "" -out_path = f"{obj._params['output_dir']}/alpha_leakage_m_s0_s4.png" +out_path = f"{obj_scale._params['output_dir']}/alpha_leakage_m_s0_s4.png" cs_plots.plot_data_1d( theta_arr, @@ -283,13 +279,13 @@ # Including scalar leakage for comparison theta_arr.append(theta) -y_arr.append(obj.alpha_leak) -dy_arr.append(obj.sig_alpha_leak) -labels.append(r"$\alpha$ (scalar approx.)") -colors.append("blue") -markers.append("p") +y_arr.append(obj_scale.alpha_leak) +dy_arr.append(obj_scale.sig_alpha_leak) +labels.append(r"$\alpha^{\rm s}_{+}$") +colors.append("cyan") +markers.append("x") linestyles.append("--") -out_path = f"{obj._params['output_dir']}/alpha_leakage_m_s0_s4_as.png" +out_path = f"{obj_scale._params['output_dir']}/alpha_leakage_m_s0_s4_as.png" cs_plots.plot_data_1d( theta_arr, @@ -309,47 +305,49 @@ close_fig=False, shift_x=True, ) - +# + \xi_ {"incorrectly_encoded_metadata": "{22}^\\textrm{g,p} \\, \\xi_{11}^\\textrm{p,p}"} # #### Consistency relations for scalar leakage # If the leakage is a scalar function, it can be expressed in three different ways. # + -xlim = [ - obj._params["theta_min_amin"] / fac, - obj._params["theta_max_amin"] * fac -] -ylim = obj._params["leakage_alpha_ylim"] - alpha_1 = [] alpha_2 = [] for ndx in range(len(theta)): - my_a1 = obj.Xi_gp_ufloat[ndx][0, 0] / obj.Xi_pp_ufloat[ndx][0, 0] + my_a1 = obj_scale.Xi_gp_ufloat[ndx][0, 0] / obj_scale.Xi_pp_ufloat[ndx][0, 0] alpha_1.append(my_a1) - my_a2 = obj.Xi_gp_ufloat[ndx][1, 1] / obj.Xi_pp_ufloat[ndx][1, 1] + my_a2 = obj_scale.Xi_gp_ufloat[ndx][1, 1] / obj_scale.Xi_pp_ufloat[ndx][1, 1] alpha_2.append(my_a2) +# TODO: Use centered functions for all cases + y = [ - unumpy.nominal_values(x0), - obj.alpha_leak, + obj_scale.alpha_leak, unumpy.nominal_values(alpha_1), unumpy.nominal_values(alpha_2), + unumpy.nominal_values(obj_scale._alpha_0_r), ] dy = [ - unumpy.std_devs(x0), - obj.sig_alpha_leak, + obj_scale.sig_alpha_leak, unumpy.std_devs(alpha_1), unumpy.std_devs(alpha_2), + unumpy.std_devs(obj_scale._alpha_0_r), ] -theta_arr = [obj.get_theta()] * len(y) +theta_arr = [theta] * len(y) - -labels = [r"$\alpha_0^\Re$", r"$\alpha_+$", r"$\alpha_1$", r"$\alpha_2$"] +labels = [ + r"$\alpha^{\rm s}_+$", + r"$\alpha^{\rm s}_1$", + r"$\alpha^{\rm s}_2$", + r"$\alpha^\Re_0$", +] xlabel = r"$\theta$ [arcmin]" ylabel = r"$\alpha(\theta)$" -linestyles = ["-", "--", "--", "--"] title = "" -out_path = f"{obj._params['output_dir']}/alpha_leakage_scalar_consistency.png" +out_path = f"{obj_scale._params['output_dir']}/alpha_leakage_scalar_consistency.png" +colors = ["cyan", "black", "grey", "blue"] +markers = ["x", "h", "o", "o"] +linestyles = ["--", "--", "--", "-"] cs_plots.plot_data_1d( theta_arr, @@ -364,25 +362,23 @@ xlim=xlim, ylim=ylim, close_fig=False, - linestyles=linestyles, shift_x=True, + colors=colors, + markers=markers, + linestyles=linestyles, ) # - # If alpha is a scalar, the mixed-component centered cross-correlation functions should be identical. # + -xlim = [obj._params["theta_min_amin"], obj._params["theta_max_amin"]] - -theta_arr = [theta] * 3 - Xi_12 = [] Xi_21 = [] Xi_tr = [] for ndx in range(len(theta)): - Xi_12.append(obj.Xi_gp_ufloat[ndx][0, 1]) - Xi_21.append(obj.Xi_gp_ufloat[ndx][1, 0]) - Xi_tr.append(obj.Xi_gp_ufloat[ndx][0, 0] + obj.Xi_gp_ufloat[ndx][0, 0]) + Xi_12.append(obj_scale.Xi_gp_ufloat[ndx][0, 1]) + Xi_21.append(obj_scale.Xi_gp_ufloat[ndx][1, 0]) + Xi_tr.append(obj_scale.Xi_gp_ufloat[ndx][0, 0] + obj_scale.Xi_gp_ufloat[ndx][0, 0]) y = [ unumpy.nominal_values(Xi_12), @@ -394,12 +390,15 @@ unumpy.std_devs(Xi_21), unumpy.std_devs(Xi_tr), ] +theta_arr = [theta] * len(y) -labels = [r"$\tau_{12,0}$", r"$\tau_{21,0}$", r"tr$\tau_0$"] +labels = [r"$\tau_{0, 12}$", r"$\tau_{0,21}$", r"tr$\tau_{0}$"] xlabel = r"$\theta$ [arcmin]" -ylabel = r"correlation functions $\tau_{ij,0}$" +ylabel = r"correlation functions $\tau_{0,ij}$" title = "" -out_path = f"{obj._params['output_dir']}/Xi_mixed_consistency.png" +out_path = f"{obj_scale._params['output_dir']}/Xi_mixed_consistency.png" +markers = ["o", "s", "d"] +linestyles = ["-", "--", ":"] cs_plots.plot_data_1d( theta_arr, @@ -414,22 +413,20 @@ xlim=xlim, close_fig=False, shift_x=True, + markers=markers, + linestyles=linestyles, ) # + # For comparison, plot the same for the PSF - PSF correlations -xlim = [obj._params["theta_min_amin"], obj._params["theta_max_amin"]] - -theta_arr = [theta] * 3 - Xi_12 = [] Xi_21 = [] Xi_tr = [] for ndx in range(len(theta)): - Xi_12.append(obj.Xi_pp_ufloat[ndx][0, 1]) - Xi_21.append(obj.Xi_pp_ufloat[ndx][1, 0]) - Xi_tr.append(obj.Xi_pp_ufloat[ndx][0, 0] + obj.Xi_pp_ufloat[ndx][0, 0]) + Xi_12.append(obj_scale.Xi_pp_ufloat[ndx][0, 1]) + Xi_21.append(obj_scale.Xi_pp_ufloat[ndx][1, 0]) + Xi_tr.append(obj_scale.Xi_pp_ufloat[ndx][0, 0] + obj_scale.Xi_pp_ufloat[ndx][0, 0]) y = [ unumpy.nominal_values(Xi_12), @@ -441,13 +438,15 @@ unumpy.std_devs(Xi_21), unumpy.std_devs(Xi_tr), ] +theta_arr = [theta] * len(y) - -labels = [r"$\rho_{12,0}$", r"$\rho_{21,0}$", r"tr$\rho_0$"] +labels = [r"$\rho_{0,12}$", r"$\rho_{0,21}$", r"tr$\rho_0$"] xlabel = r"$\theta$ [arcmin]" ylabel = r"Centered correlation functions" title = "" -out_path = f"{obj._params['output_dir']}/Xi_pp_mixed_consistency.png" +out_path = f"{obj_scale._params['output_dir']}/Xi_pp_mixed_consistency.png" +markers = ["o", "s", "d"] +linestyles = ["-", "--", ":"] cs_plots.plot_data_1d( theta_arr, @@ -462,4 +461,9 @@ xlim=xlim, close_fig=False, shift_x=True, + markers=markers, + linestyles=linestyles, ) +# - + + diff --git a/notebooks/leakage_xi_sys.py b/notebooks/leakage_xi_sys.py index dc5691b..1e25819 100644 --- a/notebooks/leakage_xi_sys.py +++ b/notebooks/leakage_xi_sys.py @@ -70,7 +70,7 @@ # \left[ \rho_{0, 11}(\theta) - \rho_{0, 22}(\theta) \right] # \nonumber \\ # & -# + 4 \left( \alpha_0^\Re \alpha_4^\Im + \alpha_4^\Re \alpha_0^\Im \right) +# + 4 \left( \alpha_0^\Re \alpha_4^\Im - \alpha_4^\Re \alpha_0^\Im \right) # \rho_{0, 12}(\theta) # , # \end{align} @@ -158,8 +158,6 @@ ) ) -print(alpha_0_r[0], alpha_0_i[0], alpha_4_r[0], alpha_4_i[0]) - # - @@ -272,7 +270,10 @@ def get_rho_0(self, idx, jdx): labels = ["$t_+$", "$t_-$", r"$t_{\rm mixed}$", r"$\sum t_i$", "scalar"] markers = ["o", "d", "^", "x", "s"] -ylim = [-1e-6, 2e-6] +factor = 0.9 +xlim = [obj_scale._params["theta_min_amin"] ** factor, obj_scale._params["theta_max_amin"]] + +#ylim = [-1e-6, 2e-6] cs_plots.plot_data_1d( x, @@ -286,7 +287,7 @@ def get_rho_0(self, idx, jdx): shift_x=True, xlog=True, ylog=False, - ylim=ylim, + xlim=xlim, markers=markers, ) # + @@ -363,13 +364,12 @@ def get_rho_0(self, idx, jdx): shift_x=False, xlog=True, ylog=False, + xlim=xlim, ylim=ylim, markers=markers, linestyles=linestyles, ) - - # + # Plot difference, should be consistent with zero @@ -417,6 +417,7 @@ def get_rho_0(self, idx, jdx): shift_x=True, xlog=True, ylog=False, + xlim=xlim, markers=markers, linestyles=linestyles, ) diff --git a/shear_psf_leakage/run_scale.py b/shear_psf_leakage/run_scale.py index f78b087..ca0f5c0 100644 --- a/shear_psf_leakage/run_scale.py +++ b/shear_psf_leakage/run_scale.py @@ -1157,6 +1157,38 @@ def get_alpha_ufloat(self, idx, jdx): mat.append(self.alpha_leak_ufloat[ndx][idx, jdx]) return np.array(mat) + + def compute_alpha_spin_coeffs(self): + """Compute Alpha Spin Coefficients. + + Compute the spin coefficients of the PSF leakage alpha(theta) from + the matrix elements. + + """ + self._alpha_0_r = ( + 0.5 * ( + self.get_alpha_ufloat(0, 0) + + self.get_alpha_ufloat(1, 1) + ) + ) + self._alpha_0_i = ( + 0.5 * ( + -self.get_alpha_ufloat(0, 1) + + self.get_alpha_ufloat(1, 0) + ) + ) + self._alpha_4_r = ( + 0.5 * ( + self.get_alpha_ufloat(0, 0) + - self.get_alpha_ufloat(1, 1) + ) + ) + self._alpha_4_i = ( + 0.5 * ( + self.get_alpha_ufloat(0, 1) + + self.get_alpha_ufloat(1, 0) + ) + ) def do_alpha_matrix(self): """Do Alpha Matrix.