diff --git a/examples/DC2_redMaPPer_clusterensemble_object.ipynb b/examples/DC2_redMaPPer_clusterensemble_object.ipynb new file mode 100644 index 0000000..9a84f89 --- /dev/null +++ b/examples/DC2_redMaPPer_clusterensemble_object.ipynb @@ -0,0 +1,577 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4d6efa8a-7255-40b8-9787-c9731bf78172", + "metadata": {}, + "source": [ + "## Compute stacked excess surface density profile arround redMaPPer clusters from the extragalactic cosmoDC2 catalog\n", + "\n", + "This notebook can be run at NERSC or CC-IN2P3 where the DESC DC2 products are stored. You need to be a DESC member to be able to access those. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a2ce472-8077-48ce-8e7c-8d5fe2338978", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from astropy.table import Table, hstack, vstack\n", + "import random\n", + "import pickle\n", + "from clmm import GalaxyCluster, ClusterEnsemble, GCData\n", + "import healpy, glob\n", + "import clmm\n", + "import matplotlib.pyplot as plt\n", + "\n", + "\n", + "def load(filename, **kwargs):\n", + " with open(filename, \"rb\") as fin:\n", + " return pickle.load(fin, **kwargs)\n", + "\n", + "\n", + "np.random.seed(1)" + ] + }, + { + "cell_type": "markdown", + "id": "059d0d80-3a32-4709-bacf-a10ec2eb070e", + "metadata": {}, + "source": [ + "### 1. Extract the catalog of redMaPPer clusters detected in the cosmoDC2 extragalactig catalog\n", + "- `python extract_cluster_catalog_redMaPPer.py` extract the catalog of redMaPPer-detected clusters\n", + "- The output catalog in saved at\n", + "-- `/global/cfs/projectdirs/lsst/groups/CL/CLMM/DC2_data/cosmoDC2_v1.1.4_redmapper_v0.8.1_catalog.pkl` (if you use NERSC).\n", + "-- `/sps/lsst/groups/clusters/CLMM/DC2_data/cosmoDC2_v1.1.4_redmapper_v0.8.1_catalog.pkl` (if you use NERSC)\n", + "- THe figure below displays the redMaPPer catalog" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a027780-fdf9-4203-b187-f54b8736891c", + "metadata": {}, + "outputs": [], + "source": [ + "path_DC2_data_ccin2p3 = \"/sps/lsst/groups/clusters/CLMM/DC2_data/\"\n", + "path_DC2_data_nersc = \"/global/cfs/projectdirs/lsst/groups/CL/CLMM/DC2_data/\"\n", + "path_DC2_data = path_DC2_data_ccin2p3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "53a53ba9-1451-42a3-90b4-3635b655dc61", + "metadata": {}, + "outputs": [], + "source": [ + "redmapper_cat = load(path_DC2_data + \"cosmoDC2_v1.1.4_redmapper_v0.8.1_catalog.pkl\")\n", + "plt.hist2d(\n", + " redmapper_cat[\"redshift\"],\n", + " redmapper_cat[\"richness\"],\n", + " cmin=1,\n", + " bins=[10, np.logspace(np.log10(10), np.log10(300), 50)],\n", + " cmap=\"gist_rainbow\",\n", + ")\n", + "plt.colorbar(label=\"count\")\n", + "plt.yscale(\"log\")\n", + "plt.xlabel(\"redshift\", fontsize=15)\n", + "plt.ylabel(r\"$\\lambda$\", fontsize=15)\n", + "plt.tick_params(axis=\"both\", which=\"major\", labelsize=15)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "cfa567d0-9c1b-4aad-a26a-d2fd6666f5a2", + "metadata": {}, + "source": [ + "### 2. Extract background galaxy sources for a sub-sample of the redMaPPer clusters and estimate individual lensing profiles\n", + "- We extract background source catalog for each redMaPPer-detected cluster. We consider a sub-sample of the redMaPPer cluster catalog, wiht richness $20 < \\lambda < 40$ and redshift $0.2 < z < 0.4$. This sub-sample consists in $n_{\\rm clusters}\\sim 450$ clusters. \n", + "- We consider the source selection mag-$i < 24.25$ and mag-$r < 28$, to reach a galaxy density $n_{\\rm gal} = 26$ gal/arcmin$^2$. \n", + "- The background galaxies are selected using the simple cut on the true redshift of the source galaxy $z_{s}$, such as $z_{s} > z_{\\rm cl}+0.2$. \n", + "- Then, we estimate the lensing profile for each cluster in 15 logarithm equally-spaced radial bins between $0.5$ Mpc and $30$ Mpc. " + ] + }, + { + "cell_type": "markdown", + "id": "eb82124c-6709-4490-ab4e-06ac41474902", + "metadata": {}, + "source": [ + "#### Extraction of galaxy catalogs from cosmoDC2\n", + "\n", + "- To speed up the extraction of cosmoDC2 source galaxies and lensing profile estimation, we subdivide the cluster sample in 20 independant splits. We then run the extraction for each independent split (containing $\\sim 25$ clusters).\n", + "- `python extract_sources_and_compute_ind_profile_cosmoDC2_per_split.py --which_split 0 --number_of_splits 30` is an example, which extracts background source catalogs and lensing profiles for clusters in the split with index `0`.\n", + "-- The extraction of cosmoDC2 data can be done using GCRCatalogs. Here, we use Qserv available at CC-IN2P3, that provides a fast interface to extract cosmoDC2 quantities (only available at CC-IN2P3). Qserv has design optimized for astronomical queries (parallel distributed SQL database).\n", + "-- Some useful functions to deal with cosmoDC2 data are stored in `_utils_cosmoDC2.py`.\n", + "- We then sublmit 20 parrallel jobs for each idependant split. `python jobs_extract_sources_cosmoDC2.py` writes and runs 20 parrallel jobs at CC-IN2P3.\n", + "- The individual profiles (stored by split from 0 to 19) are saved in the directory `./DC2_data/individual_redMaPPer_cluster_lensing_profiles/`.\n", + "- We group the individual lensing profiles in a single file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76d586a5-0990-472e-afff-d6fedfb1b6b2", + "metadata": {}, + "outputs": [], + "source": [ + "path_to_splits = path_DC2_data + \"individual_redMaPPer_cluster_lensing_profiles/\"\n", + "list_files = glob.glob(path_to_splits + \"individual_profiles_split*\")\n", + "for i, name_file in enumerate(list_files):\n", + " file = load(\n", + " name_file,\n", + " )\n", + " if i == 0:\n", + " file_stack = file\n", + " continue\n", + " else:\n", + " file_stack = vstack([file_stack, file])" + ] + }, + { + "cell_type": "markdown", + "id": "2453672d-5fb0-41b6-9760-053fd8b29fd4", + "metadata": {}, + "source": [ + "### 3. Creating `ClusterEnsemble object`\n", + "- We add the individual lensing profiles to the `ClusterEnsemble object`\n", + "- We create three ensemble object, respectively with:\n", + "-- clusters within the full range $20 < \\lambda < 40$ and $0.2 < z < 0.4$ (450 objects)\n", + "-- 100 randomly selected clusters within $20 < \\lambda < 40$ and $0.2 < z < 0.4$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "246ea896-ca6e-42a7-b906-0b53d8705833", + "metadata": {}, + "outputs": [], + "source": [ + "# ClusterEnsemble object for the full sample\n", + "cluster_ensemble = ClusterEnsemble(\"stack full sample\")\n", + "dat = clmm.GCData(\n", + " [\n", + " file_stack[\"cluster_id\"],\n", + " file_stack[\"cluster_ra\"],\n", + " file_stack[\"cluster_dec\"],\n", + " file_stack[\"cluster_redshift\"],\n", + " file_stack[\"radius\"],\n", + " file_stack[\"binned_DS_t\"],\n", + " file_stack[\"binned_DS_x\"],\n", + " file_stack[\"W_l\"],\n", + " ],\n", + " names=(\"cluster_id\", \"ra\", \"dec\", \"z\", \"radius\", \"DS_t\", \"DS_x\", \"W_l\"),\n", + " masked=True,\n", + ")\n", + "cluster_ensemble.data = dat" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3953a669-75b1-4983-b16a-7cb63e66fdf8", + "metadata": {}, + "outputs": [], + "source": [ + "index_cluster_full = np.arange(len(file_stack))\n", + "index_cluster_small = np.random.choice(index_cluster_full, 100, replace=False)\n", + "file_stack_small = file_stack[index_cluster_small]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c0ed7e9-5504-4fb6-8383-7da6aa4c3d29", + "metadata": {}, + "outputs": [], + "source": [ + "# ClusterEnsemble object for the small sample\n", + "cluster_ensemble_small = ClusterEnsemble(\"stack reduced sample (low-lambda)\")\n", + "dat = clmm.GCData(\n", + " [\n", + " file_stack_small[\"cluster_id\"],\n", + " file_stack_small[\"cluster_ra\"],\n", + " file_stack_small[\"cluster_dec\"],\n", + " file_stack_small[\"cluster_redshift\"],\n", + " file_stack_small[\"radius\"],\n", + " file_stack_small[\"binned_DS_t\"],\n", + " file_stack_small[\"binned_DS_x\"],\n", + " file_stack_small[\"W_l\"],\n", + " ],\n", + " names=(\"cluster_id\", \"ra\", \"dec\", \"z\", \"radius\", \"DS_t\", \"DS_x\", \"W_l\"),\n", + " masked=True,\n", + ")\n", + "cluster_ensemble_small.data = dat" + ] + }, + { + "cell_type": "markdown", + "id": "29ebfcd7-6419-4582-a4e3-07c558ca96dd", + "metadata": {}, + "source": [ + "### 4. Estimating the stacked profile and the covariance matrices\n", + "- We compute the stacked lensing profiles for the considered cluster sample. \n", + "- We estimate the stacked profile covariance matrix using three different data-based methods, namely the Sample, the Bootstrap and the Jackknife approaches as detailled in the [example notebook](https://github.com/LSSTDESC/CLMM/blob/main/examples/demo_mock_ensemble.ipynb) of the `ClusterEnsemble object`. Here, we use $N_{\\rm side}=16$ and $N_{\\rm bootstrap} = 500$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f970031e-e9b2-49ff-9197-51f67002f4b5", + "metadata": {}, + "outputs": [], + "source": [ + "# stacked profile\n", + "cluster_ensemble.make_stacked_radial_profile(\n", + " tan_component=\"DS_t\", cross_component=\"DS_x\"\n", + ")\n", + "r_stack, gt_stack, gx_stack = (\n", + " cluster_ensemble.stacked_data[c] for c in (\"radius\", \"DS_t\", \"DS_x\")\n", + ")\n", + "\n", + "cluster_ensemble_small.make_stacked_radial_profile(\n", + " tan_component=\"DS_t\", cross_component=\"DS_x\"\n", + ")\n", + "r_stack_small, gt_stack_small, gx_stack_small = (\n", + " cluster_ensemble_small.stacked_data[c] for c in (\"radius\", \"DS_t\", \"DS_x\")\n", + ")\n", + "\n", + "# covariance matrices\n", + "cluster_ensemble.compute_sample_covariance(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", + "cluster_ensemble.compute_bootstrap_covariance(\n", + " tan_component=\"DS_t\", cross_component=\"DS_x\", n_bootstrap=500\n", + ")\n", + "cluster_ensemble.compute_jackknife_covariance(\n", + " n_side=16, tan_component=\"DS_t\", cross_component=\"DS_x\"\n", + ")\n", + "\n", + "cluster_ensemble_small.compute_sample_covariance(\n", + " tan_component=\"DS_t\", cross_component=\"DS_x\"\n", + ")\n", + "cluster_ensemble_small.compute_bootstrap_covariance(\n", + " tan_component=\"DS_t\", cross_component=\"DS_x\", n_bootstrap=500\n", + ")\n", + "cluster_ensemble_small.compute_jackknife_covariance(\n", + " n_side=16, tan_component=\"DS_t\", cross_component=\"DS_x\"\n", + ")\n", + "\n", + "\n", + "err_gt = cluster_ensemble.cov[\"tan_sc\"].diagonal() ** 0.5 / 1e13\n", + "err_gx = cluster_ensemble.cov[\"cross_sc\"].diagonal() ** 0.5 / 1e13\n", + "\n", + "err_gt_small = cluster_ensemble_small.cov[\"tan_sc\"].diagonal() ** 0.5 / 1e13\n", + "err_gx_small = cluster_ensemble_small.cov[\"cross_sc\"].diagonal() ** 0.5 / 1e13" + ] + }, + { + "cell_type": "markdown", + "id": "095d6eda-75c9-4651-b67b-ab76111a2e3a", + "metadata": {}, + "source": [ + "In the Figure below, we plot:\n", + "- Upper panel: \n", + "-- the stacked excess surface density profile for the full sample in black, with error-bars coming from the Bootstrap covariance matrix\n", + "-- the stacked excess surface density profiles for the small cluster sample\n", + "- Lower panel:\n", + "-- the same as the upper plot, but considering the cross component of the source galaxy ellipticities." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33eed8bc-831e-446d-bcfa-ec903714c08d", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(2, 1, figsize=(6, 6), gridspec_kw={\"height_ratios\": [2, 1]})\n", + "\n", + "axs[0].fill_between(\n", + " r_stack,\n", + " gt_stack / 1e13 - err_gt,\n", + " gt_stack / 1e13 + err_gt,\n", + " color=\"C0\",\n", + " zorder=100,\n", + " alpha=0.7,\n", + " label=r\"Stacked profile (full sample)\",\n", + ")\n", + "axs[0].plot(\n", + " r_stack,\n", + " gt_stack / 1e13,\n", + " color=\"C0\",\n", + " zorder=100,\n", + " alpha=0.5,\n", + ")\n", + "\n", + "axs[0].fill_between(\n", + " r_stack,\n", + " gt_stack_small / 1e13 - err_gt_small,\n", + " gt_stack_small / 1e13 + err_gt_small,\n", + " color=\"C1\",\n", + " zorder=100,\n", + " alpha=0.5,\n", + " label=r\"Stacked profile (small sample)\",\n", + ")\n", + "axs[0].plot(\n", + " r_stack_small,\n", + " gt_stack_small / 1e13,\n", + " color=\"C1\",\n", + " zorder=100,\n", + " alpha=0.5,\n", + ")\n", + "\n", + "# for i in range(400):\n", + "# axs[0].plot(\n", + "# cluster_ensemble.data[\"radius\"][i],\n", + "# cluster_ensemble.data[\"DS_t\"][i] / 1e13,\n", + "# \"C1\",\n", + "# alpha=.1,\n", + "# linewidth=1,)\n", + "# axs[1].plot(\n", + "# cluster_ensemble.data[\"radius\"][i],\n", + "# cluster_ensemble.data[\"DS_x\"][i] / 1e13,\n", + "# \"C1\",\n", + "# alpha=.1,\n", + "# linewidth=1,)\n", + "\n", + "axs[0].loglog(\n", + " [],\n", + " [],\n", + ")\n", + "\n", + "axs[1].fill_between(\n", + " r_stack,\n", + " gx_stack / 1e13 - err_gx,\n", + " gx_stack / 1e13 + err_gx,\n", + " color=\"C0\",\n", + " zorder=100,\n", + " alpha=0.7,\n", + " label=r\"Stacked profile (full sample)\",\n", + ")\n", + "axs[1].plot(\n", + " r_stack,\n", + " gx_stack / 1e13,\n", + " color=\"C0\",\n", + " zorder=100,\n", + " alpha=0.5,\n", + ")\n", + "\n", + "axs[1].fill_between(\n", + " r_stack,\n", + " gx_stack_small / 1e13 - err_gx_small,\n", + " gx_stack_small / 1e13 + err_gx_small,\n", + " color=\"C1\",\n", + " zorder=100,\n", + " alpha=0.5,\n", + " label=r\"Stacked profile (small sample)\",\n", + ")\n", + "axs[1].plot(\n", + " r_stack_small,\n", + " gx_stack_small / 1e13,\n", + " color=\"C1\",\n", + " zorder=100,\n", + " alpha=0.5,\n", + ")\n", + "x = np.linspace(-10, 50, 100)\n", + "axs[1].plot(x, 0 * x, \"--\")\n", + "axs[1].set_xscale(\"log\")\n", + "axs[0].set_xscale(\"log\")\n", + "axs[1].set_xlabel(\"R [Mpc]\", fontsize=15)\n", + "axs[0].tick_params(axis=\"both\", which=\"major\", labelsize=15)\n", + "axs[1].tick_params(axis=\"both\", which=\"major\", labelsize=15)\n", + "axs[0].set_ylabel(r\"$\\Delta\\Sigma_+$ $[\\times 10^{13} M_\\odot /Mpc^2])$\", fontsize=15)\n", + "axs[1].set_ylabel(r\"$\\Delta\\Sigma_\\times \\times 10^{13}$ \", fontsize=15)\n", + "axs[0].grid(True, which=\"both\")\n", + "axs[1].grid(True, which=\"both\")\n", + "axs[0].set_ylim(0.03, 10)\n", + "axs[1].set_ylim(-0.7, 0.7)\n", + "axs[0].set_xlim(0.5, 30)\n", + "axs[1].set_xlim(0.5, 30)\n", + "axs[0].legend()\n", + "# plt.savefig(f'redmapper_clusters_lensing_stacked_profile_DC2.png', bbox_inches='tight', dpi=200)" + ] + }, + { + "cell_type": "markdown", + "id": "086dcc79-4438-44cf-b1b2-0ee020cf7618", + "metadata": {}, + "source": [ + "In the Figure below, we plot the diagonal elements of the three covariance matrices, estimated respectively using the sample, bootstrap and Jackknife approach." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d92eb003-7678-4815-87b3-a19ded959a8c", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(7, 6))\n", + "plt.plot(\n", + " cluster_ensemble.data[\"radius\"][0],\n", + " (cluster_ensemble.cov[\"tan_jk\"].diagonal() ** 0.5) / 1e13,\n", + " \"-s\",\n", + " label=\"Jackknife\",\n", + " color=\"g\",\n", + " linewidth=2,\n", + ")\n", + "\n", + "plt.plot(\n", + " cluster_ensemble.data[\"radius\"][0],\n", + " (cluster_ensemble.cov[\"tan_bs\"].diagonal() ** 0.5) / 1e13,\n", + " \"-\",\n", + " color=\"r\",\n", + " label=\"Bootstrap\",\n", + " linewidth=1.5,\n", + " markersize=4,\n", + ")\n", + "plt.plot(\n", + " cluster_ensemble.data[\"radius\"][0],\n", + " (cluster_ensemble.cov[\"tan_sc\"].diagonal() ** 0.5) / 1e13,\n", + " \"--\",\n", + " color=\"C0\",\n", + " label=\"Sample\",\n", + " linewidth=2,\n", + ")\n", + "plt.loglog(\n", + " cluster_ensemble_small.data[\"radius\"][0],\n", + " (cluster_ensemble_small.cov[\"tan_bs\"].diagonal() ** 0.5) / 1e13,\n", + " \"--\",\n", + " c=\"r\",\n", + " label=\"Bootstrap (small sample)\",\n", + " linewidth=1.5,\n", + " markersize=4,\n", + ")\n", + "plt.xlabel(\"R [Mpc]\", fontsize=20)\n", + "plt.ylabel(r\"$C_{ii}^{1/2}$ $[\\times 10^{13} M_\\odot /Mpc^2])$\", fontsize=15)\n", + "plt.tick_params(axis=\"both\", which=\"major\", labelsize=15)\n", + "plt.legend(frameon=True, fontsize=12)\n", + "plt.minorticks_on()\n", + "plt.grid(lw=0.5)\n", + "plt.ylim(0.01, 1)\n", + "plt.grid(\n", + " which=\"both\",\n", + ")\n", + "plt.xscale(\"log\")\n", + "# plt.savefig(f'redmapper_clusters_lensing_stacked_profile_DC2_variance.png', bbox_inches='tight', dpi=200)" + ] + }, + { + "cell_type": "markdown", + "id": "4a7c79c8-b14d-42c5-ad49-9b0d5da98e0e", + "metadata": {}, + "source": [ + "Similarly, we show the associated correlation matrices." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73641494-017a-4a6b-a110-d156e1ac9181", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(1, 3, figsize=(20, 5))\n", + "fig.subplots_adjust(wspace=0.15, hspace=0)\n", + "\n", + "maximum = cluster_ensemble.cov[\"tan_sc\"].max()\n", + "for ax, cov, label in zip(\n", + " axes,\n", + " [\n", + " cluster_ensemble.cov[\"tan_sc\"],\n", + " cluster_ensemble.cov[\"tan_bs\"],\n", + " cluster_ensemble.cov[\"tan_jk\"],\n", + " ],\n", + " [\"Stack : Sample\", \"Stack : Bootstrap\", \"Stack : Jackknife\"],\n", + "):\n", + " ax.set_title(label, fontsize=20)\n", + " ax.set_xlabel(\"radial bin index\", fontsize=18)\n", + " ax.set_ylabel(\"radial bin index\", fontsize=18)\n", + " ax.tick_params(axis=\"both\", which=\"major\", labelsize=12)\n", + " im = ax.imshow(np.corrcoef(cov), cmap=\"coolwarm\", vmin=-1, vmax=1, origin=\"lower\")\n", + " plt.colorbar(im, ax=ax)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "0a102c0e-48af-4a2c-afe4-04f142c4baa9", + "metadata": {}, + "source": [ + "(to save just one figure)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e10dd16d-9dc9-4d03-af4b-f37a1f61708d", + "metadata": {}, + "outputs": [], + "source": [ + "plt.xlabel(\"R [Mpc]\", fontsize=18)\n", + "plt.ylabel(\"R [Mpc]\", fontsize=18)\n", + "plt.tick_params(axis=\"both\", which=\"major\", labelsize=12)\n", + "im = plt.imshow(\n", + " np.corrcoef(cluster_ensemble.cov[\"tan_bs\"]),\n", + " cmap=\"coolwarm\",\n", + " vmin=-1,\n", + " vmax=1,\n", + " origin=\"lower\",\n", + ")\n", + "ticks = np.array([0, 2, 4, 6, 8, 10, 12, 14])\n", + "labels = [f\"{r_stack[t]:.1f}\" for t in ticks]\n", + "plt.xticks(\n", + " labels=labels,\n", + " ticks=ticks,\n", + ")\n", + "plt.yticks(\n", + " labels=labels,\n", + " ticks=ticks,\n", + ")\n", + "\n", + "plt.colorbar()\n", + "# plt.savefig(f'redmapper_clusters_lensing_stacked_profile_DC2_correlation_matrix.png', bbox_inches='tight', dpi=200)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "687790d0-95fb-4de4-866b-21bfe1105276", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a481498-00ab-4b70-9f77-1213feaea267", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "myenv", + "language": "python", + "name": "myenv" + }, + "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.10.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/extract_dc2_data/_utils_cosmoDC2.py b/examples/extract_dc2_data/_utils_cosmoDC2.py new file mode 100644 index 0000000..c67ba62 --- /dev/null +++ b/examples/extract_dc2_data/_utils_cosmoDC2.py @@ -0,0 +1,172 @@ +import cmath +import pickle + +import clmm +import GCRCatalogs +import healpy +import numpy as np +import pandas as pd +from astropy.table import QTable, Table, join, vstack + +GCRCatalogs.set_root_dir_by_site("in2p3") +import mysql +from mysql.connector import Error + + +def _fix_axis_ratio(q_bad): + + # from https://github.com/LSSTDESC/gcr-catalogs/blob/ellipticity_bug_fix/GCRCatalogs/cosmodc2.py + # back out incorrect computation of q using Johnsonb function + e_jb = np.sqrt((1 - q_bad**2) / (1 + q_bad**2)) + q_new = np.sqrt( + (1 - e_jb) / (1 + e_jb) + ) # use correct relationship to compute q from e_jb + return q_new + + +def _fix_ellipticity_disk_or_bulge(ellipticity): + + # from https://github.com/LSSTDESC/gcr-catalogs/blob/ellipticity_bug_fix/GCRCatalogs/cosmodc2.py + # back out incorrect computation of q using Johnsonb function + q_bad = (1 - ellipticity) / ( + 1 + ellipticity + ) # use default e definition to calculate q + # q_bad incorrectly computed from e_jb using q_bad = sqrt((1 - e_jb^2)/(1 + e_jb^2)) + q_new = _fix_axis_ratio(q_bad) + e_new = (1 - q_new) / ( + 1 + q_new + ) # recompute e using default (1-q)/(1+q) definition + return e_new + + +def correct_shear_ellipticity(ellipticity_uncorr_e1, ellipticity_uncorr_e2): + + # from https://github.com/LSSTDESC/gcr-catalogs/blob/ellipticity_bug_fix/GCRCatalogs/cosmodc2.py + ellipticity_uncorr_norm = ( + ellipticity_uncorr_e1**2 + ellipticity_uncorr_e2**2 + ) ** 0.5 + complex_ellipticity_uncorr = ellipticity_uncorr_e1 + 1j * ellipticity_uncorr_e2 + phi = np.array([cmath.phase(c) for c in complex_ellipticity_uncorr]) + ellipticity_corr_norm = _fix_ellipticity_disk_or_bulge(ellipticity_uncorr_norm) + ellipticity_corr = ellipticity_corr_norm * np.exp(1j * phi) + ellipticity_corr_e1, ellipticity_corr_e2 = ( + ellipticity_corr.real, + ellipticity_corr.imag, + ) + return ellipticity_corr_e1, ellipticity_corr_e2 + + +def extract_cosmoDC2_galaxy(lens_z, lens_distance, ra, dec, rmax=10, method="qserv"): + + if method == "qserv": + + def qserv_query(lens_z, lens_distance, ra, dec, rmax=10): + r""" + quantities wanted + cuts for qserv + Attributes: + ----------- + z: float + lens redshift + lens_distance: float + distance to the cluster + ra: float + lens right ascension + dec: float + lens declinaison + rmax: float + maximum radius + """ + zmax = 3.0 + zmin = lens_z + theta_max = (rmax / lens_distance) * (180.0 / np.pi) + query = "SELECT data.coord_ra as ra, data.coord_dec as dec, data.redshift as z, " + query += "data.galaxy_id as galaxy_id, data.halo_id as halo_id, " + query += "data.mag_i, data.mag_r, data.mag_y, " + query += "data.shear_1 as shear1, data.shear_2 as shear2, data.convergence as kappa, " + query += "data.ellipticity_1_true as e1_true_uncorr, data.ellipticity_2_true as e2_true_uncorr " + query += "FROM cosmoDC2_v1_1_4_image.data as data " + query += f"WHERE data.redshift >= {zmin} AND data.redshift < {zmax} " + query += f"AND scisql_s2PtInCircle(coord_ra, coord_dec, {ra}, {dec}, {theta_max}) = 1 " + query += f"AND data.mag_i <= 24.6 " + query += f"AND data.mag_r <= 28.0 " + query += ";" + return query + + conn_qserv = mysql.connector.connect( + host="ccqserv201", user="qsmaster", port=30040 + ) + cursor = conn_qserv.cursor(dictionary=True, buffered=True) + qserv_query_extract = qserv_query(z, lens_distance, ra, dec, rmax=rmax) + tab = pd.read_sql_query(qserv_query_extract, conn_qserv) + tab = QTable.from_pandas(tab) + return tab + + if method == "gcrcatalogs": + + def gcrcatalogs_query(): + + return [ + "galaxy_id", + "ra", + "dec", + "shear_1", + "shear_2", + "mag_i", + "mag_r", + "mag_y", + "redshift", + "convergence", + "halo_id", + "ellipticity_1_true", + "ellipticity_2_true", + ] + + theta_apperture_deg = (180 / np.pi) * (rmax / lens_distance) + n = 1000 + dx = np.random.random(n) * (-2) + 1 + dy = np.random.random(n) * (-2) + 1 + x_rand, y_rand = ra + theta_apperture_deg * dx, dec + theta_apperture_deg * dy + healpix = np.unique(healpy.ang2pix(32, x_rand, y_rand, nest=False, lonlat=True)) + healpix_dc2 = GCRCatalogs.load_catalog( + "cosmoDC2_v1.1.4_image" + ).get_catalog_info()["healpix_pixels"] + healpix_list = healpix[np.isin(healpix, healpix_dc2)] + + cosmoDC2 = GCRCatalogs.load_catalog("cosmoDC2_v1.1.4_image") + + coord_filters = [ + "ra >= {}".format(ra - theta_apperture_deg / 2), + "ra < {}".format(ra + theta_apperture_deg / 2), + "dec >= {}".format(dec - theta_apperture_deg / 2), + "dec < {}".format(dec + theta_apperture_deg / 2), + ] + + z_filters = ["redshift >= {}".format(lens_z)] + mag_filters = ["mag_i < {}".format(24.6)] + ["mag_r < {}".format(28.0)] + + tab = Table(names=gcrcatalogs_query()) + + for i, hp in enumerate(healpix_list): + print(f"-----> heapix pixel = " + str(hp)) + chunk = cosmoDC2.get_quantities( + gcrcatalogs_query(), + filters=(coord_filters + z_filters + mag_filters), + native_filters=[f"healpix_pixel=={hp}"], + return_iterator=True, + ) + for j in range(3): + print("chunk = " + str(j)) + try: + dat_extract_chunk = Table(next(chunk)) + except: + continue + tab = vstack([tab, dat_extract_chunk]) + + tab.rename_column("ellipticity_1_true", "e1_true_uncorr") + tab.rename_column("ellipticity_2_true", "e2_true_uncorr") + tab.rename_column("shear_1", "shear1") + tab.rename_column("shear_2", "shear2") + tab.rename_column("convergence", "kappa") + tab.rename_column("redshift", "z") + + return tab diff --git a/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py b/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py new file mode 100644 index 0000000..cc23071 --- /dev/null +++ b/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py @@ -0,0 +1,74 @@ +import GCRCatalogs + +GCRCatalogs.set_root_dir_by_site("in2p3") +import pickle +import sys + +import astropy.units as u +import matplotlib.pyplot as plt +import numpy as np +from astropy.cosmology import FlatLambdaCDM +from astropy.table import Table + + +def load(filename, **kwargs): + """Loads GalaxyCluster object to filename using Pickle""" + with open(filename, "rb") as fin: + return pickle.load(fin, **kwargs) + + +def save_pickle(dat, filename, **kwargs): + file = open(filename, "wb") + pickle.dump(dat, file) + file.close() + + +catalog = GCRCatalogs.load_catalog("cosmoDC2_v1.1.4_redmapper_v0.8.1") +quantity = [ + "cluster_id", + "ra", + "dec", + "redshift", + "redshift_err", + "richness", + "richness_err", + "ra_member", + "dec_member", + "id_member", + "cluster_id_member", +] +dat = catalog.get_quantities(quantity) + +line_ra_member = [] +line_dec_member = [] +line_id_member = [] +line_cluster_id_member = [] +for i in range(len(dat["cluster_id"])): + mask = dat["cluster_id_member"] == dat["cluster_id"][i] + line_ra_member.append(dat["ra_member"][mask]) + line_dec_member.append(dat["dec_member"][mask]) + line_id_member.append(dat["id_member"][mask]) + line_cluster_id_member.append(dat["cluster_id_member"][mask]) + +dictionnary = { + "cluster_id": dat["cluster_id"], + "ra": dat["ra"], + "dec": dat["dec"], + "redshift": dat["redshift"], + "redshift_err": dat["redshift_err"], + "richness": dat["richness"], + "richness_err": dat["richness_err"], + "ra_member": line_ra_member, + "dec_member": line_dec_member, + "id_member": line_id_member, + "cluster_id_member": line_cluster_id_member, +} + +t = Table(dictionnary) + +path_DC2_data_ccin2p3 = "/sps/lsst/groups/clusters/CLMM/" +path_DC2_data_nersc = "/global/cfs/projectdirs/lsst/groups/CL/CLMM/" + +path_data = path_DC2_data_nersc + +save_pickle(t, path_data + "./DC2_data/cosmoDC2_v1.1.4_redmapper_v0.8.1_catalog.pkl") diff --git a/examples/extract_dc2_data/extract_sources_and_compute_ind_profile_cosmoDC2_per_split.py b/examples/extract_dc2_data/extract_sources_and_compute_ind_profile_cosmoDC2_per_split.py new file mode 100644 index 0000000..0a2a00c --- /dev/null +++ b/examples/extract_dc2_data/extract_sources_and_compute_ind_profile_cosmoDC2_per_split.py @@ -0,0 +1,277 @@ +# -*- coding: utf-8 -*- +import glob +import pickle +import sys + +import astropy +import astropy.units as u +import clmm +import GCRCatalogs +# GCRCatalogs.set_root_dir_by_site('in2p3') +import healpy +import numpy as np +import pandas as pd +from astropy.coordinates import SkyCoord, match_coordinates_sky +from astropy.io import fits as fits +from astropy.table import QTable, Table, hstack, join, vstack + + +def load(filename, **kwargs): + """Loads GalaxyCluster object to filename using Pickle""" + with open(filename, "rb") as fin: + return pickle.load(fin, **kwargs) + + +def save_pickle(dat, filename, **kwargs): + file = open(filename, "wb") + pickle.dump(dat, file) + file.close() + + +from clmm import Cosmology +from scipy.integrate import simps + +# cosmoDC2 cosmology +cosmo = Cosmology(H0=71.0, Omega_dm0=0.265 - 0.0448, Omega_b0=0.0448, Omega_k0=0.0) +import argparse + +import _utils_cosmoDC2 +# connection with qserv +import mysql +from mysql.connector import Error + +mag_i_max = 24.25 +mag_r_max = 28 + + +def source_selection_magnitude(table): + mask = table["mag_i"] < mag_i_max + mask *= table["mag_r"] < mag_r_max + return table[mask] + + +def collect_argparser(): + parser = argparse.ArgumentParser( + description="Let's extract source catalogs and estimate individual lensing profiles" + ) + parser.add_argument("--which_split", type=int, required=True) + parser.add_argument("--number_of_splits", type=int, required=True) + parser.add_argument( + "--lens_catalog_name", + type=str, + required=False, + default="./DC2_data/cosmoDC2_v1.1.4_redmapper_v0.8.1_catalog.pkl", + ) + parser.add_argument( + "--lambda_min", + type=float, + required=False, + default=20, + ) + parser.add_argument( + "--lambda_max", + type=float, + required=False, + default=40, + ) + parser.add_argument( + "--redshift_min", + type=float, + required=False, + default=0.2, + ) + parser.add_argument( + "--redshift_max", + type=float, + required=False, + default=0.4, + ) + parser.add_argument( + "--method", + type=str, + required=False, + default="qserv", + ) + parser.add_argument( + "--rmax", + type=float, + required=False, + default=30, + ) + return parser.parse_args() + + +# select galaxy clusters +_config_extract_sources_in_cosmoDC2 = collect_argparser() +lens_catalog_name = _config_extract_sources_in_cosmoDC2.lens_catalog_name +lens_catalog = np.load(lens_catalog_name, allow_pickle=True) +mask_select = lens_catalog["richness"] > _config_extract_sources_in_cosmoDC2.lambda_min +mask_select *= lens_catalog["richness"] < _config_extract_sources_in_cosmoDC2.lambda_max +mask_select *= ( + lens_catalog["redshift"] > _config_extract_sources_in_cosmoDC2.redshift_min +) +mask_select *= ( + lens_catalog["redshift"] < _config_extract_sources_in_cosmoDC2.redshift_max +) +lens_catalog = lens_catalog[mask_select] +lens_catalog_truncated = lens_catalog + +n_cl = len(lens_catalog_truncated) +index_cl = np.arange(n_cl) +split_lists = np.array_split( + index_cl, int(_config_extract_sources_in_cosmoDC2.number_of_splits) +) +lens_catalog_truncated = lens_catalog_truncated[ + np.array(split_lists[int(_config_extract_sources_in_cosmoDC2.which_split)]) +] +n_cl_to_extract = len(lens_catalog_truncated) + +path_DC2_data_ccin2p3 = "/sps/lsst/groups/clusters/CLMM/" +path_DC2_data_nersc = "/global/cfs/projectdirs/lsst/groups/CL/CLMM/" +path_CLMM_data = path_DC2_data_nersc + +path_where_to_save_profiles = ( + path_CLMM_data + f"./DC2_data/individual_redMaPPer_cluster_lensing_profiles/" +) +name_save = path_where_to_save_profiles + f"individual_profiles" +name_save += f"_split={_config_extract_sources_in_cosmoDC2.which_split}" +name_save += f"_number_of_splits={_config_extract_sources_in_cosmoDC2.number_of_splits}_ncl={n_cl_to_extract}.pkl" + +print( + f"[cluster catalog]: {_config_extract_sources_in_cosmoDC2.lambda_min} < richness < {_config_extract_sources_in_cosmoDC2.lambda_max}" +) +print( + f"[cluster catalog]: {_config_extract_sources_in_cosmoDC2.redshift_min} < redshift < {_config_extract_sources_in_cosmoDC2.redshift_max}" +) +print(f"[cluster catalog]: number of clusters = {n_cl}") +print( + f"[cluster catalog]: number of clusters in split {_config_extract_sources_in_cosmoDC2.which_split} = {n_cl_to_extract}" +) +print(f"[individual lensing profiles]: individual profiles " + name_save) + +for n, lens in enumerate(lens_catalog_truncated): + + z, ra, dec = lens["redshift"], lens["ra"], lens["dec"] + cluster_id = lens["cluster_id"] + richness = lens["richness"] + print(f"[load background source catalog]: for redMaPPer cluster {cluster_id}") + # cluster member galaxies + id_member_galaxy = lens_catalog["id_member"][n] + ra_member_galaxy = lens_catalog["ra_member"][n] + dec_member_galaxy = lens_catalog["dec_member"][n] + lens_distance = cosmo.eval_da(z) + + tab = _utils_cosmoDC2.extract_cosmoDC2_galaxy( + z, + lens_distance, + ra, + dec, + rmax=_config_extract_sources_in_cosmoDC2.rmax, + method=_config_extract_sources_in_cosmoDC2.method, + ) + # compute reduced shear and ellipticities + tab["g1"], tab["g2"] = clmm.utils.convert_shapes_to_epsilon( + tab["shear1"], tab["shear2"], shape_definition="shear", kappa=tab["kappa"] + ) + ellipticity_uncorr_e1 = tab["e1_true_uncorr"] + ellipticity_uncorr_e2 = tab["e2_true_uncorr"] + ellipticity_corr_e1, ellipticity_corr_e2 = ( + _utils_cosmoDC2.correct_shear_ellipticity( + ellipticity_uncorr_e1, ellipticity_uncorr_e2 + ) + ) + tab["e1_true"] = ellipticity_corr_e1 + tab["e2_true"] = ellipticity_corr_e2 + tab["e1"], tab["e2"] = clmm.utils.compute_lensed_ellipticity( + tab["e1_true"], tab["e2_true"], tab["shear1"], tab["shear2"], tab["kappa"] + ) + + dat_extract_mag_cut = source_selection_magnitude(tab) + mask_z = dat_extract_mag_cut["z"] > z + 0.2 + dat_extract_mag_cut_z_cut = dat_extract_mag_cut[mask_z] + + dat = clmm.GCData( + [ + dat_extract_mag_cut_z_cut["ra"], + dat_extract_mag_cut_z_cut["dec"], + dat_extract_mag_cut_z_cut["e1"], + dat_extract_mag_cut_z_cut["e2"], + dat_extract_mag_cut_z_cut["z"], + dat_extract_mag_cut_z_cut["galaxy_id"], + ], + names=("ra", "dec", "e1", "e2", "z", "id"), + masked=True, + ) + + cl = clmm.GalaxyCluster("redmapper_cluster", ra, dec, z, dat) + + bin_edges = clmm.dataops.make_bins( + 0.5, _config_extract_sources_in_cosmoDC2.rmax, 15, method="evenlog10width" + ) + + cl.compute_tangential_and_cross_components( + shape_component1="e1", + shape_component2="e2", + tan_component="DS_t", + cross_component="DS_x", + cosmo=cosmo, + is_deltasigma=True, + use_pdz=False, + ) + + cl.compute_galaxy_weights( + use_pdz=False, + use_shape_noise=False, + shape_component1="e1", + shape_component2="e2", + use_shape_error=False, + shape_component1_err="e_err", + shape_component2_err="e_err", + weight_name="w_ls", + cosmo=cosmo, + is_deltasigma=True, + add=True, + ) + + cl.make_radial_profile( + "Mpc", + bins=bin_edges, + error_model="ste", + cosmo=cosmo, + tan_component_in="DS_t", + cross_component_in="DS_x", + tan_component_out="binned_DS_t", + cross_component_out="binned_DS_x", + tan_component_in_err=None, + cross_component_in_err=None, + include_empty_bins=False, + gal_ids_in_bins=False, + add=True, + table_name="profile", + overwrite=True, + use_weights=True, + weights_in="w_ls", + weights_out="W_l", + ) + + binned_DS_t = cl.profile["binned_DS_t"] + binned_DS_x = cl.profile["binned_DS_x"] + W_l = cl.profile["W_l"] + radius = cl.profile["radius"] + cluster_data = [cluster_id, ra, dec, z, lens["richness"]] + cluster_data_colnames = [ + "cluster_id", + "cluster_ra", + "cluster_dec", + "cluster_redshift", + "cluster_richness", + ] + if n == 0: + ind_profile = {k: [] for k in list(cl.profile.colnames) + cluster_data_colnames} + + for index, name in enumerate(cl.profile.colnames): + ind_profile[name].append(cl.profile[name]) + for index, name in enumerate(cluster_data_colnames): + ind_profile[name].append(cluster_data[index]) + +save_pickle(Table(ind_profile), name_save) diff --git a/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py b/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py new file mode 100644 index 0000000..59eb82d --- /dev/null +++ b/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py @@ -0,0 +1,27 @@ +import os +import sys + +import numpy as np + +number_of_splits = 20 +for which_split in np.arange(number_of_splits): + lines_base = [ + "#!/bin/sh", + "# SBATCH options:", + "#SBATCH --job-name=extract_cosmodc2_redmapper # Job name", + f"#SBATCH --output=./log/{which_split}_nsplits_{number_of_splits}.log", + "#SBATCH --partition=htc # Partition choice", + "#SBATCH --ntasks=1 # Run a single task (by default tasks == CPU)", + "#SBATCH --mem=8000 # Memory in MB per default", + "#SBATCH --time=0-6:00:00 # 7 days by default on htc partition", + "source /pbs/home/c/cpayerne/setup_mydesc.sh", + ] + cmd = f"python extract_sources_and_compute_ind_profile_cosmoDC2_per_split.py --which_split {which_split} --number_of_splits {number_of_splits} " + lines = lines_base + [cmd] + name_job = f"job_which_split={which_split}_number_of_splits={number_of_splits}.job" + with open(name_job, "w") as f: + for line in lines: + f.write(line) + f.write("\n") + # os.system(f'sbatch {name_job}') + os.remove(name_job) diff --git a/examples/version.py b/examples/version.py index f102a9c..3b93d0b 100644 --- a/examples/version.py +++ b/examples/version.py @@ -1 +1 @@ -__version__ = "0.0.1" +__version__ = "0.0.2"