v1.0

.gitignore

@@ -5,5 +5,6 @@ notebooks/leakage_object
 notebooks/unions_shapepipe_extended_2022_W3_v1.0.3.fits
 __pycache__
 __pycache__/*
+shear_psf_leakage/__pycache__/*
 */__pycache__/*
 build

README.md

@@ -19,12 +19,14 @@ pypi.org,
 pip install shear_psf_leakage
 ```
 
-Alternatively, you can use
+Alternatively, you can use TBD
 
 ## Authors
 
 Martin Kilbinger <martin.kilbinger@cea.fr>
 
-Clara Bonini <mail>
+Clara Bonini
 
-Axel Guinot <mail>
+Sacha Guerrini
+
+Axel Guinot

notebooks/leakage_autocorr_psf.py

@@ -0,0 +1,206 @@
+# ---
+# 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 <martin.kilbinger@cea.fr>
+
+# %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,
+)
+# -
+
+

notebooks/leakage_object.py

@@ -15,81 +15,48 @@
 
 # # PSF leakage: object-wise estimation
 #
-# Demonstration notebook of shear_psf_leakage.run_object.LeakageObject() class.
+# Demonstration notebook of shear_psf_leakage.run_object.LeakageObject class.
 #
 # Martin Kilbinger <martin.kilbinger@cea.fr>
 
 # +
-import os
+import matplotlib
+matplotlib.use("agg")
+
+import os, sys
 import matplotlib.pylab as plt
 
 from cs_util import canfar
+from cs_util import args
 
 import shear_psf_leakage.run_object as run
+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_W3_v1.0.3.fits"
-
-# Output directory
-params_in["output_dir"] = "leakage_object"
-# -
-
-# ### Job control
-
-# +
-# Compute (spin-preserving) PSF ellipticity leakage
-params_in["PSF_leakage"] = True
-
-# Compute leakage with other parameters
-params_in["obs_leakage"] = True
-
-# Other input parameters
-params_in["cols"] = "RA Dec e1_PSF fwhm_PSF w mag snr"
-
-# Ratio between two input columns
-params_in["cols_ratio"] = "mag_snr"
-# -
-
-# ### 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 = run.LeakageObject()
 
 # Set instance parameters, copy from above
-for key in params_in:
-    obj._params[key] = params_in[key]
+params_upd = args.read_param_script("params_object.py", obj._params, verbose=True)
+for key in params_upd:
+    obj._params[key] = params_upd[key]
+
+#for key in obj._params:
+    #print(key, obj._params[key])
+
+# +
+#vos_dir = f"vos:cfis/weak_lensing/DataReleases/v1.0/ShapePipe/{patch}"
+#canfar.download(
+    #f"{vos_dir}/{obj_params['input_path_shear']}",
+    #obj._params["input_path_shear"],
+    #verbose=obj._params["verbose"],
+#)
+# -
 
 # ### Run
-# There are two options to run the leakage compuations:
+# There are two options to run the leakage computations:
 # 1. Run all at once with single class routine.
 # 2. Execute individual steps.
 
@@ -116,7 +83,8 @@
 
 if obj._params["PSF_leakage"]:
     # Object-by-object spin-consistent PSF leakage
-    obj.PSF_leakage()
+    for order in ("lin", "quad"):
+            obj.PSF_leakage(mix=True, order=order)
 
 if obj._params["obs_leakage"]:
     # Object-by-object spin-consistent PSF leakage