From 1e79a3a6408a610ef31f330171bb454db53d28a7 Mon Sep 17 00:00:00 2001 From: payerne Date: Fri, 3 Apr 2026 23:34:40 +0200 Subject: [PATCH 1/4] add scripts and notebook for DC2 examples --- ...DC2_redMaPPer_clusterensemble_object.ipynb | 518 ++++++++++++++++++ examples/extract_dc2_data/_utils_cosmoDC2.py | 140 +++++ .../extract_cluster_catalog_redMaPPer.py | 56 ++ ..._compute_ind_profile_cosmoDC2_per_split.py | 209 +++++++ .../jobs_extract_sources_cosmoDC2.py | 23 + 5 files changed, 946 insertions(+) create mode 100644 examples/DC2_redMaPPer_clusterensemble_object.ipynb create mode 100644 examples/extract_dc2_data/_utils_cosmoDC2.py create mode 100644 examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py create mode 100644 examples/extract_dc2_data/extract_sources_and_compute_ind_profile_cosmoDC2_per_split.py create mode 100644 examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py diff --git a/examples/DC2_redMaPPer_clusterensemble_object.ipynb b/examples/DC2_redMaPPer_clusterensemble_object.ipynb new file mode 100644 index 0000000..8356dfd --- /dev/null +++ b/examples/DC2_redMaPPer_clusterensemble_object.ipynb @@ -0,0 +1,518 @@ +{ + "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", + "def load(filename, **kwargs):\n", + " with open(filename, 'rb') as fin:\n", + " return pickle.load(fin, **kwargs)\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(redmapper_cat['redshift'],redmapper_cat['richness'], cmin=1, bins = [10, np.logspace(np.log10(10), np.log10(300), 50)], cmap = 'gist_rainbow')\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(name_file,)\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", + " [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", + " names=(\"cluster_id\", \"ra\", \"dec\", \"z\", \"radius\", \"DS_t\", \"DS_x\", 'W_l'),\n", + " masked=True,)\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", + " [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", + " names=(\"cluster_id\", \"ra\", \"dec\", \"z\", \"radius\", \"DS_t\", \"DS_x\", 'W_l'),\n", + " masked=True,)\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(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", + "r_stack, gt_stack, gx_stack = (cluster_ensemble.stacked_data[c] for c in (\"radius\", \"DS_t\", \"DS_x\"))\n", + "\n", + "cluster_ensemble_small.make_stacked_radial_profile(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", + "r_stack_small, gt_stack_small, gx_stack_small = (cluster_ensemble_small.stacked_data[c] for c in (\"radius\", \"DS_t\", \"DS_x\"))\n", + "\n", + "#covariance matrices\n", + "cluster_ensemble.compute_sample_covariance(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", + "cluster_ensemble.compute_bootstrap_covariance(tan_component=\"DS_t\", cross_component=\"DS_x\", n_bootstrap=500)\n", + "cluster_ensemble.compute_jackknife_covariance(n_side=16, tan_component=\"DS_t\", cross_component=\"DS_x\")\n", + "\n", + "cluster_ensemble_small.compute_sample_covariance(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", + "cluster_ensemble_small.compute_bootstrap_covariance(tan_component=\"DS_t\", cross_component=\"DS_x\", n_bootstrap=500)\n", + "cluster_ensemble_small.compute_jackknife_covariance(n_side=16, tan_component=\"DS_t\", cross_component=\"DS_x\")\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", + "\n", + " alpha=.7,\n", + " label=r\"Stacked profile (full sample)\",)\n", + "axs[0].plot(\n", + " r_stack,\n", + " gt_stack / 1e13,\n", + " color=\"C0\",\n", + " zorder=100,\n", + "\n", + " alpha=.5,)\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", + "\n", + " alpha=.5,\n", + " label=r\"Stacked profile (small sample)\",)\n", + "axs[0].plot(\n", + " r_stack_small,\n", + " gt_stack_small / 1e13,\n", + " color=\"C1\",\n", + " zorder=100,\n", + "\n", + " alpha=.5,)\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", + "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", + "\n", + " alpha=.7,\n", + " label=r\"Stacked profile (full sample)\",)\n", + "axs[1].plot(\n", + " r_stack,\n", + " gx_stack / 1e13,\n", + " color=\"C0\",\n", + " zorder=100,\n", + "\n", + " alpha=.5,)\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", + "\n", + " alpha=.5,\n", + " label=r\"Stacked profile (small sample)\",)\n", + "axs[1].plot(\n", + " r_stack_small,\n", + " gx_stack_small / 1e13,\n", + " color=\"C1\",\n", + " zorder=100,\n", + "\n", + " alpha=.5,)\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 ,'-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", + " \"-\",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", + " \"--\",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(which=\"both\",)\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", + " [cluster_ensemble.cov[\"tan_sc\"], cluster_ensemble.cov[\"tan_bs\"], cluster_ensemble.cov[\"tan_jk\"]],\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(np.corrcoef(cluster_ensemble.cov[\"tan_bs\"]), cmap=\"coolwarm\", vmin=-1, vmax=1, origin=\"lower\")\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(labels = labels, ticks=ticks, )\n", + "plt.yticks(labels = labels, ticks=ticks, )\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..042cad5 --- /dev/null +++ b/examples/extract_dc2_data/_utils_cosmoDC2.py @@ -0,0 +1,140 @@ +import numpy as np +from astropy.table import QTable, Table, vstack, join +import pickle +import pandas as pd +import clmm +import cmath +import healpy +import GCRCatalogs +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)**.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. + zmin = lens_z + theta_max = (rmax/lens_distance) * (180./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 + + + + + \ No newline at end of file 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..7a99f91 --- /dev/null +++ b/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py @@ -0,0 +1,56 @@ +import GCRCatalogs +GCRCatalogs.set_root_dir_by_site('in2p3') +import matplotlib.pyplot as plt +import pickle +import sys +import numpy as np +from astropy.cosmology import FlatLambdaCDM +import astropy.units as u +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..96747cb --- /dev/null +++ b/examples/extract_dc2_data/extract_sources_and_compute_ind_profile_cosmoDC2_per_split.py @@ -0,0 +1,209 @@ +# -*- coding: utf-8 -*- +import numpy as np +import GCRCatalogs +#GCRCatalogs.set_root_dir_by_site('in2p3') +import healpy +import glob, sys +import astropy.units as u +from astropy.io import fits as fits +from astropy.coordinates import SkyCoord, match_coordinates_sky +import astropy +import clmm +import pandas as pd + +from astropy.table import QTable, Table, vstack, join, hstack +import pickle, sys + +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) +#connection with qserv +import mysql +from mysql.connector import Error +import argparse +import _utils_cosmoDC2 + +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..b479f1a --- /dev/null +++ b/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py @@ -0,0 +1,23 @@ +import os, 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) From 258751f929f48304aa076c5d4f8ea7a48ed521da Mon Sep 17 00:00:00 2001 From: Michel Aguena Date: Thu, 9 Apr 2026 14:53:29 +0200 Subject: [PATCH 2/4] run black --- ...DC2_redMaPPer_clusterensemble_object.ipynb | 263 ++++++++----- examples/extract_dc2_data/_utils_cosmoDC2.py | 165 ++++---- .../extract_cluster_catalog_redMaPPer.py | 81 ++-- ..._compute_ind_profile_cosmoDC2_per_split.py | 369 +++++++++++------- .../jobs_extract_sources_cosmoDC2.py | 30 +- 5 files changed, 542 insertions(+), 366 deletions(-) diff --git a/examples/DC2_redMaPPer_clusterensemble_object.ipynb b/examples/DC2_redMaPPer_clusterensemble_object.ipynb index 8356dfd..9a84f89 100644 --- a/examples/DC2_redMaPPer_clusterensemble_object.ipynb +++ b/examples/DC2_redMaPPer_clusterensemble_object.ipynb @@ -26,9 +26,12 @@ "import clmm\n", "import matplotlib.pyplot as plt\n", "\n", + "\n", "def load(filename, **kwargs):\n", - " with open(filename, 'rb') as fin:\n", + " with open(filename, \"rb\") as fin:\n", " return pickle.load(fin, **kwargs)\n", + "\n", + "\n", "np.random.seed(1)" ] }, @@ -52,8 +55,8 @@ "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_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" ] }, @@ -64,12 +67,18 @@ "metadata": {}, "outputs": [], "source": [ - "redmapper_cat = load(path_DC2_data + 'cosmoDC2_v1.1.4_redmapper_v0.8.1_catalog.pkl')\n", - "plt.hist2d(redmapper_cat['redshift'],redmapper_cat['richness'], cmin=1, bins = [10, np.logspace(np.log10(10), np.log10(300), 50)], cmap = 'gist_rainbow')\n", - "plt.colorbar(label='count')\n", - "plt.yscale('log')\n", - "plt.xlabel('redshift', fontsize=15)\n", - "plt.ylabel(r'$\\lambda$', fontsize=15)\n", + "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()" ] @@ -109,11 +118,13 @@ "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", + "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(name_file,)\n", - " if i==0: \n", + " file = load(\n", + " name_file,\n", + " )\n", + " if i == 0:\n", " file_stack = file\n", " continue\n", " else:\n", @@ -139,19 +150,22 @@ "metadata": {}, "outputs": [], "source": [ - "#ClusterEnsemble object for the full sample\n", - "cluster_ensemble = ClusterEnsemble('stack full sample')\n", + "# ClusterEnsemble object for the full sample\n", + "cluster_ensemble = ClusterEnsemble(\"stack full sample\")\n", "dat = clmm.GCData(\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", - " names=(\"cluster_id\", \"ra\", \"dec\", \"z\", \"radius\", \"DS_t\", \"DS_x\", 'W_l'),\n", - " masked=True,)\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" ] }, @@ -174,19 +188,22 @@ "metadata": {}, "outputs": [], "source": [ - "#ClusterEnsemble object for the small sample\n", - "cluster_ensemble_small = ClusterEnsemble('stack reduced sample (low-lambda)')\n", + "# ClusterEnsemble object for the small sample\n", + "cluster_ensemble_small = ClusterEnsemble(\"stack reduced sample (low-lambda)\")\n", "dat = clmm.GCData(\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", - " names=(\"cluster_id\", \"ra\", \"dec\", \"z\", \"radius\", \"DS_t\", \"DS_x\", 'W_l'),\n", - " masked=True,)\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" ] }, @@ -207,21 +224,39 @@ "metadata": {}, "outputs": [], "source": [ - "#stacked profile\n", - "cluster_ensemble.make_stacked_radial_profile(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", - "r_stack, gt_stack, gx_stack = (cluster_ensemble.stacked_data[c] for c in (\"radius\", \"DS_t\", \"DS_x\"))\n", + "# 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(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", - "r_stack_small, gt_stack_small, gx_stack_small = (cluster_ensemble_small.stacked_data[c] for c in (\"radius\", \"DS_t\", \"DS_x\"))\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", + "# covariance matrices\n", "cluster_ensemble.compute_sample_covariance(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", - "cluster_ensemble.compute_bootstrap_covariance(tan_component=\"DS_t\", cross_component=\"DS_x\", n_bootstrap=500)\n", - "cluster_ensemble.compute_jackknife_covariance(n_side=16, 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(tan_component=\"DS_t\", cross_component=\"DS_x\")\n", - "cluster_ensemble_small.compute_bootstrap_covariance(tan_component=\"DS_t\", cross_component=\"DS_x\", n_bootstrap=500)\n", - "cluster_ensemble_small.compute_jackknife_covariance(n_side=16, tan_component=\"DS_t\", cross_component=\"DS_x\")\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", @@ -251,41 +286,41 @@ "metadata": {}, "outputs": [], "source": [ - "fig, axs = plt.subplots(2, 1, figsize=(6, 6), gridspec_kw={'height_ratios': [2,1]})\n", + "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", + " gt_stack / 1e13 - err_gt,\n", + " gt_stack / 1e13 + err_gt,\n", " color=\"C0\",\n", " zorder=100,\n", - "\n", - " alpha=.7,\n", - " label=r\"Stacked profile (full sample)\",)\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", - "\n", - " alpha=.5,)\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", + " gt_stack_small / 1e13 - err_gt_small,\n", + " gt_stack_small / 1e13 + err_gt_small,\n", " color=\"C1\",\n", " zorder=100,\n", - "\n", - " alpha=.5,\n", - " label=r\"Stacked profile (small sample)\",)\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", - "\n", - " alpha=.5,)\n", + " alpha=0.5,\n", + ")\n", "\n", "# for i in range(400):\n", "# axs[0].plot(\n", @@ -300,59 +335,62 @@ "# \"C1\",\n", "# alpha=.1,\n", "# linewidth=1,)\n", - " \n", - "axs[0].loglog([], [], )\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", + " gx_stack / 1e13 - err_gx,\n", + " gx_stack / 1e13 + err_gx,\n", " color=\"C0\",\n", " zorder=100,\n", - "\n", - " alpha=.7,\n", - " label=r\"Stacked profile (full sample)\",)\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", - "\n", - " alpha=.5,)\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", + " gx_stack_small / 1e13 - err_gx_small,\n", + " gx_stack_small / 1e13 + err_gx_small,\n", " color=\"C1\",\n", " zorder=100,\n", - "\n", - " alpha=.5,\n", - " label=r\"Stacked profile (small sample)\",)\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", - "\n", - " alpha=.5,)\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].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].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)" + "# plt.savefig(f'redmapper_clusters_lensing_stacked_profile_DC2.png', bbox_inches='tight', dpi=200)" ] }, { @@ -370,33 +408,36 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(7,6))\n", + "plt.figure(figsize=(7, 6))\n", "plt.plot(\n", " cluster_ensemble.data[\"radius\"][0],\n", - " (cluster_ensemble.cov[\"tan_jk\"].diagonal() ** 0.5)/1e13 ,'-s',\n", + " (cluster_ensemble.cov[\"tan_jk\"].diagonal() ** 0.5) / 1e13,\n", + " \"-s\",\n", " label=\"Jackknife\",\n", - " color='g',\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", - " \"-\",color = 'r',\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", - " \"--\",color='C0',\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", + " (cluster_ensemble_small.cov[\"tan_bs\"].diagonal() ** 0.5) / 1e13,\n", " \"--\",\n", " c=\"r\",\n", " label=\"Bootstrap (small sample)\",\n", @@ -410,9 +451,11 @@ "plt.minorticks_on()\n", "plt.grid(lw=0.5)\n", "plt.ylim(0.01, 1)\n", - "plt.grid(which=\"both\",)\n", - "plt.xscale('log')\n", - "#plt.savefig(f'redmapper_clusters_lensing_stacked_profile_DC2_variance.png', bbox_inches='tight', dpi=200)" + "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)" ] }, { @@ -436,7 +479,11 @@ "maximum = cluster_ensemble.cov[\"tan_sc\"].max()\n", "for ax, cov, label in zip(\n", " axes,\n", - " [cluster_ensemble.cov[\"tan_sc\"], cluster_ensemble.cov[\"tan_bs\"], cluster_ensemble.cov[\"tan_jk\"]],\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", @@ -466,14 +513,26 @@ "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(np.corrcoef(cluster_ensemble.cov[\"tan_bs\"]), cmap=\"coolwarm\", vmin=-1, vmax=1, origin=\"lower\")\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(labels = labels, ticks=ticks, )\n", - "plt.yticks(labels = labels, ticks=ticks, )\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.savefig(f'redmapper_clusters_lensing_stacked_profile_DC2_correlation_matrix.png', bbox_inches='tight', dpi=200)\n", "plt.show()" ] }, diff --git a/examples/extract_dc2_data/_utils_cosmoDC2.py b/examples/extract_dc2_data/_utils_cosmoDC2.py index 042cad5..0f46722 100644 --- a/examples/extract_dc2_data/_utils_cosmoDC2.py +++ b/examples/extract_dc2_data/_utils_cosmoDC2.py @@ -1,49 +1,65 @@ import numpy as np from astropy.table import QTable, Table, vstack, join -import pickle +import pickle import pandas as pd import clmm import cmath import healpy import GCRCatalogs -GCRCatalogs.set_root_dir_by_site('in2p3') + +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 + + # 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 + 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 + + # 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 + 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)**.5 - complex_ellipticity_uncorr = ellipticity_uncorr_e1 + 1j*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 + 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): + +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: @@ -59,82 +75,97 @@ def qserv_query(lens_z, lens_distance, ra, dec, rmax = 10): rmax: float maximum radius """ - zmax = 3. + zmax = 3.0 zmin = lens_z - theta_max = (rmax/lens_distance) * (180./np.pi) + 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 += "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 += ";" + query += ";" return query - - conn_qserv = mysql.connector.connect(host='ccqserv201', user='qsmaster', port=30040) + + 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': - + + 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) + 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 + 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_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),] + "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()) + + 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) + 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: + print("chunk = " + str(j)) + try: dat_extract_chunk = Table(next(chunk)) - except: continue + 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 - - + 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") - - \ No newline at end of file + return tab diff --git a/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py b/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py index 7a99f91..8344116 100644 --- a/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py +++ b/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py @@ -1,5 +1,6 @@ import GCRCatalogs -GCRCatalogs.set_root_dir_by_site('in2p3') + +GCRCatalogs.set_root_dir_by_site("in2p3") import matplotlib.pyplot as plt import pickle import sys @@ -7,50 +8,66 @@ from astropy.cosmology import FlatLambdaCDM import astropy.units as u from astropy.table import Table + + def load(filename, **kwargs): """Loads GalaxyCluster object to filename using Pickle""" - with open(filename, 'rb') as fin: + 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'] + 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,} +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_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') +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 index 96747cb..865dbff 100644 --- 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 @@ -1,7 +1,8 @@ # -*- coding: utf-8 -*- import numpy as np import GCRCatalogs -#GCRCatalogs.set_root_dir_by_site('in2p3') + +# GCRCatalogs.set_root_dir_by_site('in2p3') import healpy import glob, sys import astropy.units as u @@ -14,196 +15,262 @@ from astropy.table import QTable, Table, vstack, join, hstack import pickle, sys + def load(filename, **kwargs): """Loads GalaxyCluster object to filename using Pickle""" - with open(filename, 'rb') as fin: + with open(filename, "rb") as fin: return pickle.load(fin, **kwargs) + def save_pickle(dat, filename, **kwargs): - file = open(filename,'wb') + 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) -#connection with qserv + +# cosmoDC2 cosmology +cosmo = Cosmology(H0=71.0, Omega_dm0=0.265 - 0.0448, Omega_b0=0.0448, Omega_k0=0.0) +# connection with qserv import mysql from mysql.connector import Error import argparse import _utils_cosmoDC2 -mag_i_max = 24.25 -mag_r_max = 28 +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 + 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 = 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,) + 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 + +# 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/' +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' +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) +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']) + 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 + 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,) + [ + 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) - - - + 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 index b479f1a..fa6764f 100644 --- a/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py +++ b/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py @@ -1,23 +1,25 @@ import os, 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} ' + "#!/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: + 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}') + f.write("\n") + # os.system(f'sbatch {name_job}') os.remove(name_job) From 2f5b2b448d26ffc40f1a99ac228cd66818273612 Mon Sep 17 00:00:00 2001 From: Michel Aguena Date: Thu, 9 Apr 2026 14:54:27 +0200 Subject: [PATCH 3/4] bump version --- examples/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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" From c99854e8d0b57993c6ea712778ccc034ddee38ad Mon Sep 17 00:00:00 2001 From: Michel Aguena Date: Thu, 9 Apr 2026 15:06:22 +0200 Subject: [PATCH 4/4] run isort --- examples/extract_dc2_data/_utils_cosmoDC2.py | 11 ++++---- .../extract_cluster_catalog_redMaPPer.py | 5 ++-- ..._compute_ind_profile_cosmoDC2_per_split.py | 27 ++++++++++--------- .../jobs_extract_sources_cosmoDC2.py | 4 ++- 4 files changed, 26 insertions(+), 21 deletions(-) diff --git a/examples/extract_dc2_data/_utils_cosmoDC2.py b/examples/extract_dc2_data/_utils_cosmoDC2.py index 0f46722..c67ba62 100644 --- a/examples/extract_dc2_data/_utils_cosmoDC2.py +++ b/examples/extract_dc2_data/_utils_cosmoDC2.py @@ -1,11 +1,12 @@ -import numpy as np -from astropy.table import QTable, Table, vstack, join +import cmath import pickle -import pandas as pd + import clmm -import cmath -import healpy 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 diff --git a/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py b/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py index 8344116..cc23071 100644 --- a/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py +++ b/examples/extract_dc2_data/extract_cluster_catalog_redMaPPer.py @@ -1,12 +1,13 @@ import GCRCatalogs GCRCatalogs.set_root_dir_by_site("in2p3") -import matplotlib.pyplot as plt import pickle import sys + +import astropy.units as u +import matplotlib.pyplot as plt import numpy as np from astropy.cosmology import FlatLambdaCDM -import astropy.units as u from astropy.table import Table 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 index 865dbff..0a2a00c 100644 --- 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 @@ -1,19 +1,19 @@ # -*- coding: utf-8 -*- -import numpy as np -import GCRCatalogs +import glob +import pickle +import sys -# GCRCatalogs.set_root_dir_by_site('in2p3') -import healpy -import glob, sys -import astropy.units as u -from astropy.io import fits as fits -from astropy.coordinates import SkyCoord, match_coordinates_sky 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.table import QTable, Table, vstack, join, hstack -import pickle, sys +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): @@ -33,11 +33,12 @@ def save_pickle(dat, filename, **kwargs): # 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 -import argparse -import _utils_cosmoDC2 mag_i_max = 24.25 mag_r_max = 28 diff --git a/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py b/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py index fa6764f..59eb82d 100644 --- a/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py +++ b/examples/extract_dc2_data/jobs_extract_sources_cosmoDC2.py @@ -1,4 +1,6 @@ -import os, sys +import os +import sys + import numpy as np number_of_splits = 20