Merge pull request #133 from CosmoStat/develop

v0.1.2

Author: martinkilbinger

.gitignore

@@ -158,3 +158,4 @@ cosmo_inference/cosmosis_config/priors_*
 cosmo_inference/cosmosis_config/values_*
 
 notebooks/tests/star_response.ipynb
+sp_validation.code-workspace

README.md

@@ -11,153 +11,43 @@ Validation of weak-lensing catalogues (galaxy and star shapes and other paramete
 | [![coc](https://img.shields.io/badge/conduct-read-lightgrey)](https://github.com/martin.kilbinger/sp_validation/blob/master/CODE_OF_CONDUCT.md) | [![Updates](https://pyup.io/repos/github/martin.kilbinger/sp_validation/shield.svg)](https://pyup.io/repos/github/martin.kilbinger/sp_validation/) | |
 
 ---
-> Author: <a href="www.cosmostat.org/people/kilbinger" target="_blank" style="text-decoration:none; color: #F08080">Axel Guinot, Martin Kilbinger, Samuel Farrens, Emma Ayçoberry</a>  
+> Authors: <a href="www.cosmostat.org" target="_blank" style="text-decoration:none; color: #F08080">CosmoStat</a> lab at CEA Paris-Saclay, including:
+  Axel Guinot, Martin Kilbinger, Lucie Baumont, Sacha Guerrini, Fabian Hervas Peters, Samuel Farrens, Emma Ayçoberry.</a>
 > Email: <a href="mailto:[email protected]" style="text-decoration:none; color: #F08080">[email protected]</a>  
-> Year: 2021  
 ---
 
-See [pyraliddemo](https://github.com/sfarrens/pyraliddemo) for a demo package created with the Pyralid template.
+This package contains a library and several scripts and notebooks. The main
+tasks that can be performed by `sp_validation` are:
+- Shear validation, in particular for the output of the `shapepipe`
+  pipeline. This task takes on input a shear catalogue with metacal information,
+  performs the calibration and carries out various tests, e.g. PSF leakage.
+  A calibrated shear catalogue is then created on output.  
+- Post processing. A number of scripts allow further processing of the above
+  output calibrated shear catalogue.  
+- Cosmology validation. This task uses the calibrated shear catalogue from
+  above to run detailed diagnostics useful for further cosmology analysis,
+  e.g. rho- and tau-statistics, E-/B-mode decomposition. Several catalogues
+  can be compared and useful plots are created.
+- Cosmology inference. This task uses the calibrated shear catalogue from
+  a shear validation run and performes cosmology inference using the two-point
+  correlation function.
 
 
-## Run validation
+## Run shear validation
 
-### Set up
+See the [documentation](docs/source/run_validation.md) for instructions how to set up and run `sp_validation`.
 
-Edit the file `notebooks/params.py` according to your data.
 
-Make sure that all input files set in `params.py` are accessible from the run directory. These are
-the ASCII file containing the tile IDs (`path_tile_ID`), the FITS galaxy catalogue (`galaxy_cat_path`),
-and the FITS star catalogue (`star_cat_path`).
+## Post processing
 
-The file `param.py` needs to be in the directory where the validation is run.
+The output(s) of one or more [shear validation runs](#run-shear-validation) can
+be processed further with post-processing scripts. See
+[here](docs/source/post_processing.md) for details.
 
-### Run
+## Cosmology validation
 
-There are two possibilities to carry out the validation, by running the jupyter notebooks
-in a browser, or by running a python script in the command line via `ipython`.
+TBD.
 
-#### Running the jupyter notebooks
-
-1. In the run directory start JupyterLab:
-  ```bash
-  jupyer-lab
-  ```
-
-2. Load and run first notebook `main_set_up.ipynb`:
-   1. Double-click in file manager tab on the left.
-   2. Run notebook.
-
-3. Run all further notebooks:
-   1. Double-click as above.
-   2. Change kernel to `main_set_up.ipynb`:
-      Either via the menu `Kernel -> Change Kernel` or
-      by clicking on `Python 3` on the top left of the notebook.
-   3. Run notebook.
-
-Run the notebooks in the following order:
-   1. `main_set_up.ipynb` (main notebook)
-   2. `metacal_global.ipynb`
-   3. [`metacal_local.ipynb`] optional     
-   4. `psf_leakage.ipynb`
-   5. `write_cat.ipynb`
-   6. `maps.ipynb`
-   7. `cosmology.ipynb`
-
-#### Running the python script
-
-1. Create the python script from the jupyter notbooks. In `notebooks`:
-  ```bash
-  python config_convert.py
-  ```
-
-2. Run python script. In run directory:
-  ```bash
-   ipython /path/to/sp_validation/notebooks/validation.py
-   ```
-
-## Further post-processing
-
-After the validation is run, further processing steps can be carried out using python scripts, as follows.
-
-### Combine validation runs
-
-Summary statistics created by validation runs of sub-areas of a survey can be combined to create joint summary statistics.
-This is useful in cases where the galaxy catalogue of an entire survey is too large to process, and needs to be broken
-down in smaller patches. This step provides global summary statistics from those patches.
-
-Depending on the type of summary, their combination can be the sum (e.g. for number of objects), average, weighted average (e.g. for the additive bias),
-the weighted average of the square (e.g. the ellipticity dispersion), the weighted variance (to combine variance estimates), or the weighted variance of the mean
-(to combine mean variance estimates).
-
-In a directory containing the subpatches as subdirectories, and within each their own `sp_output` results of the validation runs, type
-```bash
-/path/to/sp_validation/scripts/combine_results.py
-```
-This script creates a number of output files, including `R.txt` and `c.txt` with the combined multiplicative and additive biases, respectively.
-
-### Create combined calibrated shear catalogue
-
-After creating the combined results described above, the global calibration outputs can be used to create a combined, globally calibrated shear catalogue.
-The calibration is obtained from the files `R.txt` and `c.txt` created above.
-
-In the same directory containing the subpatches as above, type
-```bash
-/path/to/sp_validation/scripts/create_joint_shape_cat.py
-```
-It creates the joint output catalogue `joint.fits`.
-
-### PSF - galaxy object-wise ellipticity leakage
-
-The validation notebooks, in particular `psf_leakage.ipynb`, compute the leakage from PSF to galaxy ellipticity, as part of the global validation.
-This can also be done with the stand-alone python script `scripts/leakage_object.py`.
-
-For example, for a given patch, run
-```bash
-leakage_object.py -i sp_output/shape_catalog_extended_ngmix.fits -o leakage --hdu_psf 2 -v
-```
-to output plots and text files for the object-wise and scale-dependent leakage for that patch.
-
-Leakage for the joint v1.0 catalogue can be computed via
-```bash
-leakage_object.py -i SP/unions_shapepipe_extended_2022_v1.0.fits -I SP/unions_shapepipe_psf_2022_v1.0.1.fits -o leakage -v
-```
-assuming `SP` is a link to the v1.0 ShapePipe data directory.
-If this call was done in a subdirectory from where in `..` are the patch runs, joint plots of the scale-dependent leakage can be produced by
-```bash
-plot_leakage.py leakage/alpha_leakage_ngmix.txt ../P[1234567]/leakage/alpha_leakage_ngmix.txt`
-```
-This will read in the text files produces by the previous calls of `leakage_object.py`.
-
-This script fits a consistent linear or quadratic model of the galaxy
-ellipticity as function of PSF ellipticity. The best-fit parameters are written
-in an output `json` file.
-
-A summary table of the object-wise leakage linear parameters will be created by
-`combine_results.py`, see above. This call will add the leakage from the joint
-catalogue, assuming the results are stored in `joint/leakage`.
-
-### PSF - galaxy scale-dependent ellipticity leakage
-
-Use the script `scripts/leakage_scale.py`.
-
-### Correct for PSF - galaxy ellipticity leakage
-
-Experimentally we can apply the correction alpha * eps_PSF to the galaxy
-elliptity.
-
-First, `leakage_object.py` (see above) needs to be run.
-
-Second, the best-fit result from that previous call of `leakage_object.py` is
-read, applied to the extended shear catalogue, and the corrected ellipticity is
-written in an output file.
-
-```bash
-apply_alpha.py`
-```
-
-Third, compute the shear correlation function and E-/B-mode statistics
-using `treecorr`, with the notebook (or ipython script)
-`notebooks/analt/xip_xim_v1_alpha_cor[.ipynb|.py]`.
-
-Fourth, use the notebook `TBD` to plot the results. 
+## Cosmology inference
 
+See the corresponding [documentation](cosmo_inference/README.md).

cosmo_inference/cosmocov_config/cosmocov.ini

@@ -16,22 +16,32 @@ h0 : 0.7
 # n_gal,lens_n_gal in gals/arcmin^2
 
 ; FOR SHAPEPIPE 3500
-; area : 3218.19
+; FOR SP_v1.3_LFmask_8k
+; area : 2138
 ; sourcephotoz : multihisto
 ; lensphotoz : multihisto
 ; source_tomobins : 1
 ; lens_tomobins : 1
-; sigma_e : 0.491712
-; source_n_gal : 8.42
+; sigma_e : 0.4384
+; source_n_gal : 7.6
 
-; FOR SHAPEPIPE 1500
-area : 1453
+; FOR SP_v1.4_P1+3
+; area : 1133
+; sourcephotoz : multihisto
+; lensphotoz : multihisto
+; source_tomobins : 1
+; lens_tomobins : 1
+; sigma_e : 0.4384
+; source_n_gal : 7.67
+
+; FOR SP_V0.0
+area : 1500
 sourcephotoz : multihisto
 lensphotoz : multihisto
 source_tomobins : 1
 lens_tomobins : 1
-sigma_e : 0.4808326112068524
-source_n_gal : 7.92
+sigma_e : 0.4243
+source_n_gal : 6.47
 
 # IA parameters
 IA : 1
@@ -42,8 +52,8 @@ oneplusz0_ia: 0.6
 # Covariance paramters
 #
 # tmin,tmax in arcminutes
-tmin : 1
-tmax : 200
+tmin : 0.1
+tmax : 250
 ntheta : 20
 
 # whether to calculate the non-Gaussian part as well

cosmo_inference/notebooks/MCMC.ipynb

@@ -247,7 +247,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "my_env",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -261,7 +261,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.8"
+   "version": "3.9.0"
   },
   "vscode": {
    "interpreter": {

cosmo_inference/scripts/cosmocov_process.py

@@ -7,19 +7,30 @@
 import matplotlib.pyplot as plt
 import subprocess
 
-def calc_cov(root, nz_file):
+def calc_cov(root, nz_file, cosmocov_path):
 	Path("data/%s/covs" %root).mkdir(parents=True, exist_ok=True)
 	# write to the ini file, paths to the nz file
-	with open('cosmocov_config/cosmocov_%s.ini' %root,'a') as f:
-		f.write('shear_REDSHIFT_FILE : %s' %nz_file + '\n')
-		f.write('clustering_REDSHIFT_FILE : %s' %nz_file + '\n')
-		f.write('outdir : data/%s/covs/' %root + '\n')
-	f.close()
+	print("Creating cosmocov ini file for %s..." %root)
+	if not os.path.exists('cosmocov_config/cosmocov_%s.ini' %root):
+		with open('cosmocov_config/cosmocov_%s.ini' %root,'a') as f:
+			with open('cosmocov_config/cosmocov.ini', 'r') as f2:
+				f.write(f2.read())
+			f.write('\n')
+			f.write('shear_REDSHIFT_FILE : %s' %nz_file + '\n')
+			f.write('clustering_REDSHIFT_FILE : %s' %nz_file + '\n')
+			f.write('outdir : data/%s/covs/' %root + '\n')
+		f.close()
+
+	else:
+		print("Cosmocov ini file for %s already exists!" %root)
 
 	print("Running CosmoCov...\n")
 
+	if not os.path.exists('%s/CosmoCov/covs/cov' %cosmocov_path):
+		raise Exception("CosmoCov executable not found! Please check the path to CosmoCov in the config file.")
+		
 	results = subprocess.run('for i in {1..3}; do \
-     						CosmoCov/covs/cov $i cosmocov_config/cosmocov_%s.ini; done' %root, 
+     						%s/CosmoCov/covs/cov $i cosmocov_config/cosmocov_%s.ini; done' %(cosmocov_path, root), 
 							capture_output = True, 
 							text = True,
 							shell=True)
@@ -55,8 +66,9 @@ def convert_cov(filename):
 
 	root = sys.argv[1]
 	nz_file = sys.argv[2]
+	cosmocov_path = sys.argv[3]
 
-	calc_cov(root,nz_file)
+	calc_cov(root,nz_file, cosmocov_path)
 
 	covfile = "data/%s/covs/cov_%s" %(root,root)
 

cosmo_inference/scripts/cosmosis_fitting.py

@@ -7,7 +7,7 @@
 import sys
 
 #transforms treecorr fits file of correlation functions into CosmoSIS-friendly 2pt FITS extension to be read by 2pt_likelihood
-def treecorr_to_fits(xipm_file,root):
+def treecorr_to_fits(filename1, filename2, root):
 
     xiplus_hdu = fits.open(filename1)
     ximinus_hdu = fits.open(filename2)
@@ -118,7 +118,7 @@ def rho_to_fits(filename):
 
     #create the required FITS extensions
     print("Creating 2PT fits extension...\n")
-    xip_hdu, xim_hdu = treecorr_to_fits(two_pt_file,root)
+    xip_hdu, xim_hdu = treecorr_to_fits(two_pt_file_xip, two_pt_file_xim, root)
     print("Creating CovMat fits extension...\n")
     cov_hdu = covdat_to_fits(cov_file)
     print("Creating n(z) fits extension...\n")

cosmo_inference/scripts/xi_sys_psf.py

@@ -6,6 +6,8 @@
 #This file should be added to your cosmosis_standard_library following the path shear/xi_sys/xi_sys_psf.py
 def setup(options):
     filename = options.get_string(option_section, 'data_file')
+    xi_p_name = options.get_string(option_section, 'xi_plus_name')
+    xi_m_name = options.get_string(option_section, 'xi_minus_name')
     data = fits.open(filename)
     rho_stats_name = options.get_string(option_section, 'rho_stats_name')
     samples_path = options.get_string(option_section,'samples')
@@ -16,31 +18,37 @@ def setup(options):
 
     rho_stats = data[rho_stats_name].data
 
-    return mean, cov, rho_stats
+    range_xi_p= options[option_section, 'angle_range_XI_PLUS_1_1']
+    range_xi_m= options[option_section, 'angle_range_XI_MINUS_1_1']
+
+    mask_xi_p = (data[xi_p_name].data['ANG'] > range_xi_p[0]) & (data[xi_p_name].data['ANG'] < range_xi_p[1])
+    mask_xi_m = (data[xi_m_name].data['ANG'] > range_xi_m[0]) & (data[xi_m_name].data['ANG'] < range_xi_m[1])
+
+    return mean, cov, rho_stats, mask_xi_p, mask_xi_m
 
 def execute(block, config):
 
-    mean, cov, rho_stats = config
+    mean, cov, rho_stats, mask_xi_p, mask_xi_m = config
 
     alpha, beta, eta = np.random.multivariate_normal(mean, cov)
     block['xi_sys', 'alpha'], block['xi_sys', 'beta'], block['xi_sys', 'eta'] = alpha, beta, eta
 
     xi_sys_p = (
-        alpha**2*rho_stats["rho_0_p"]
-        + beta**2*rho_stats["rho_1_p"]
-        + eta**2*rho_stats["rho_3_p"]
-        + 2*alpha*beta*rho_stats["rho_2_p"]
-        + 2*beta*eta*rho_stats["rho_4_p"]
-        + 2*alpha*eta*rho_stats["rho_5_p"]
+        alpha**2*rho_stats["rho_0_p"][mask_xi_p]
+        + beta**2*rho_stats["rho_1_p"][mask_xi_p]
+        + eta**2*rho_stats["rho_3_p"][mask_xi_p]
+        + 2*alpha*beta*rho_stats["rho_2_p"][mask_xi_p]
+        + 2*beta*eta*rho_stats["rho_4_p"][mask_xi_p]
+        + 2*alpha*eta*rho_stats["rho_5_p"][mask_xi_p]
     )
 
     xi_sys_m = (
-        alpha**2*rho_stats["rho_0_m"]
-        + beta**2*rho_stats["rho_1_m"]
-        + eta**2*rho_stats["rho_3_m"]
-        + 2*alpha*beta*rho_stats["rho_2_m"]
-        + 2*beta*eta*rho_stats["rho_4_m"]
-        + 2*alpha*eta*rho_stats["rho_5_m"]
+        alpha**2*rho_stats["rho_0_m"][mask_xi_m]
+        + beta**2*rho_stats["rho_1_m"][mask_xi_m]
+        + eta**2*rho_stats["rho_3_m"][mask_xi_m]
+        + 2*alpha*beta*rho_stats["rho_2_m"][mask_xi_m]
+        + 2*beta*eta*rho_stats["rho_4_m"][mask_xi_m]
+        + 2*alpha*eta*rho_stats["rho_5_m"][mask_xi_m]
     )
 
     block['xi_sys', 'xi_sys_vec'] = np.concatenate([xi_sys_p, xi_sys_m])