From f98a91fbbc5fa648eb31ca9af40974e371988252 Mon Sep 17 00:00:00 2001 From: Andreas Tersenov <“[andrewtersenov@gmail.com> Date: Thu, 10 Aug 2023 15:03:07 +0200 Subject: [PATCH] Organized the code in a different way & added more functions --- example/cosmoSLICS_bookkeeping.ipynb | 469 ++++++++++++++------------- 1 file changed, 236 insertions(+), 233 deletions(-) diff --git a/example/cosmoSLICS_bookkeeping.ipynb b/example/cosmoSLICS_bookkeeping.ipynb index 6813967..1665b15 100644 --- a/example/cosmoSLICS_bookkeeping.ipynb +++ b/example/cosmoSLICS_bookkeeping.ipynb @@ -200,7 +200,11 @@ "from scipy import ndimage as ndi\n", "from lenspack.geometry.projections.gnom import radec2xy\n", "from lenspack.utils import bin2d\n", - "from lenspack.image.inversion import ks93" + "from lenspack.image.inversion import ks93\n", + "import lenspack.peaks as peaks\n", + "from lenspack.starlet_l1norm import noise_coeff, get_l1norm_noisy\n", + "from lenspack.image.transforms import starlet2d\n", + "from astropy.stats import mad_std" ] }, { @@ -209,12 +213,14 @@ "metadata": {}, "outputs": [], "source": [ - "# Define the path to the catalog file\n", - "catalog_file = \"/n17data/tersenov/SLICS/Cosmo_DES/16_a/LOS4/DES_MocksCat_16_a_4_Bin3_LOS4_R4.dat\"\n", - "\n", - "# Load the catalog data as a numpy array\n", - "catalog_data = np.loadtxt(catalog_file)\n", - "catalog_data = catalog_data.T" + "# Constants and Parameters\n", + "CATALOG_FILE = \"/n17data/tersenov/SLICS/Cosmo_DES/16_a/LOS4/DES_MocksCat_16_a_4_Bin3_LOS4_R4.dat\"\n", + "N_GAL = 7\n", + "PIX_ARCMIN = 0.4\n", + "SHAPE_NOISE = 0.44\n", + "NSCALES = 5\n", + "NBINS = 40 \n", + "KAPPA_SNR = np.linspace(-2, 6, 31)" ] }, { @@ -229,15 +235,8 @@ " g1_sim = catalog_data[6]\n", " g2_sim = catalog_data[7]\n", " kappa_sim = catalog_data[8]\n", - " return ra, dec, g1_sim, g2_sim, kappa_sim" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + " return ra, dec, g1_sim, g2_sim, kappa_sim\n", + "\n", "def create_kappa_map(ra, dec, g1_sim, g2_sim, size_x_deg=10, size_y_deg=10, pixel_size_emap_amin=0.4):\n", " x, y = radec2xy(np.mean(ra), np.mean(dec), ra, dec) # Project (ra,dec) -> (x,y)\n", "\n", @@ -248,33 +247,177 @@ " emap = np.array([e1map,e2map]) # stack the two components into a single array\n", "\n", " kappaE, kappaB = ks93(e1map, -e2map) # make kappa map (the minus sign has to be here for our data conventions)\n", - " return kappaE, kappaB" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + " return kappaE, kappaB\n", + "\n", "def add_noise_to_kappa_map(kappa_map, shape_noise, n_gal, pix_arcmin):\n", " sigma_noise_CFIS = shape_noise / (np.sqrt(2 * n_gal * pix_arcmin**2))\n", " noise_map_CFIS_z05 = sigma_noise_CFIS * np.random.randn(kappa_map.shape[0], kappa_map.shape[1]) # generate noise map\n", " kappa_map_noisy = kappa_map + noise_map_CFIS_z05 # Add noise to the mass map\n", - " return kappa_map_noisy, noise_map_CFIS_z05" + " return kappa_map_noisy, noise_map_CFIS_z05\n", + "\n", + "def smooth_kappa_map(kappa_map, pixel_size_emap_amin):\n", + " # Set the standard deviation of the Gaussian filter based on the pixel size of the kappa map\n", + " precision_Peaks = 2 / pixel_size_emap_amin # pixel_size_emap_amin is the pixel size of the kappa map in arcminutes\n", + " kappa_map_smoothed = ndi.gaussian_filter(kappa_map, precision_Peaks)\n", + " return kappa_map_smoothed\n", + "\n", + "def convert_to_snr_map(kappa_map_smoothed, noise_map_smoothed):\n", + " snr_map = kappa_map_smoothed / np.std(noise_map_smoothed)\n", + " return snr_map\n", + "\n", + "def compute_multiscale_snr_maps(image, noise, nscales):\n", + " \"\"\"\n", + " Compute SNR maps for each wavelet scale of a noisy image.\n", + " \n", + " Parameters:\n", + " image (numpy.ndarray): The noiseless image.\n", + " noise (numpy.ndarray): The noise to be added to the image.\n", + " nscales (int): Number of wavelet scales for starlet decomposition.\n", + " \n", + " Returns:\n", + " snr_maps (list of numpy.ndarray): List of SNR maps for each scale.\n", + " \"\"\"\n", + " # Add noise to the noiseless image\n", + " image_noisy = image + noise\n", + " \n", + " # Perform starlet decomposition\n", + " image_starlet = starlet2d(image_noisy, nscales)\n", + " \n", + " # Estimate the noise level\n", + " noise_estimate = mad_std(image_noisy)\n", + " coeff_j = noise_coeff(image, nscales)\n", + " \n", + " snr_maps = []\n", + " for image_j, std_co in zip(image_starlet, coeff_j):\n", + " sigma_j = std_co * noise_estimate\n", + " \n", + " # Compute SNR map\n", + " snr_map = image_j / sigma_j\n", + " snr_maps.append(snr_map)\n", + " \n", + " return snr_maps\n", + "\n", + "def compute_single_scale_peak_counts(snr_map, kappa_snr):\n", + " \"\"\"\n", + " Compute peak counts for a single SNR map.\n", + "\n", + " Parameters:\n", + " snr_map (numpy.ndarray): SNR map.\n", + " kappa_snr (numpy.ndarray): Array of kappa values corresponding to the SNR map.\n", + "\n", + " Returns:\n", + " kappa_th_center_snr (numpy.ndarray): Array of kappa threshold centers for peak counts.\n", + " peak_counts (numpy.ndarray): Peak counts for the given SNR map.\n", + " \"\"\"\n", + " kappa_th_center_snr = 0.5 * (kappa_snr[:-1] + kappa_snr[1:])\n", + " peak_counts = peaks.peaks_histogram(snr_map, kappa_snr)[0]\n", + " return kappa_th_center_snr, peak_counts\n", + "\n", + "\n", + "def compute_multiscale_peak_counts(snr_maps, kappa_snr):\n", + " \"\"\"\n", + " Compute peak counts for each wavelet scale of SNR maps.\n", + "\n", + " Parameters:\n", + " snr_maps (list of numpy.ndarray): List of SNR maps for each scale.\n", + " kappa_snr (numpy.ndarray): Array of kappa values corresponding to SNR maps.\n", + "\n", + " Returns:\n", + " kappa_th_center_snr (numpy.ndarray): Array of kappa threshold centers for peak counts.\n", + " peak_counts (list of numpy.ndarray): List of peak counts for each scale.\n", + " \"\"\"\n", + " kappa_th_center_snr = 0.5 * (kappa_snr[:-1] + kappa_snr[1:])\n", + " \n", + " peak_counts = [peaks.peaks_histogram(snr_map, kappa_snr)[0] for snr_map in snr_maps]\n", + "\n", + " return kappa_th_center_snr, peak_counts\n", + "\n", + "def plot_map(map_data, title='', cmap='inferno', vmin=None, vmax=None):\n", + " \"\"\"\n", + " Plot a 2D map using a colormap.\n", + "\n", + " Parameters:\n", + " map_data (numpy.ndarray): The 2D map data.\n", + " title (str): Title of the plot.\n", + " cmap (str): Colormap name.\n", + " vmin (float): Minimum value for colormap scaling.\n", + " vmax (float): Maximum value for colormap scaling.\n", + " \"\"\"\n", + " plt.figure()\n", + " img = plt.imshow(map_data, cmap=cmap, vmin=vmin, vmax=vmax, origin='lower')\n", + " plt.title(title)\n", + " plt.colorbar(img)\n", + " plt.show() \n", + "\n", + "def plot_peak_count_histograms(kappa_th_center, peak_counts, title, xlabel, ylabel, log_scale=False):\n", + " \"\"\"\n", + " Plot histograms of peak counts.\n", + "\n", + " Parameters:\n", + " kappa_th_center (numpy.ndarray): Array of kappa threshold centers.\n", + " peak_counts (numpy.ndarray or list of numpy.ndarray): Peak counts for each scale.\n", + " title (str): Title of the plot.\n", + " xlabel (str): Label for the x-axis.\n", + " ylabel (str): Label for the y-axis.\n", + " log_scale (bool, optional): Whether to use a logarithmic scale for the y-axis. Default is False.\n", + " \"\"\"\n", + " plt.figure()\n", + " \n", + " if isinstance(peak_counts, list): # Multiscale case\n", + " for scale, peak_count in enumerate(peak_counts):\n", + " plt.plot(kappa_th_center, peak_count, label=f'Scale {scale + 1}')\n", + " else: # Single-scale case\n", + " plt.plot(kappa_th_center, peak_counts)\n", + " \n", + " plt.legend()\n", + " if log_scale:\n", + " plt.yscale('log')\n", + " plt.xlabel(xlabel)\n", + " plt.ylabel(ylabel)\n", + " plt.grid(True)\n", + " plt.title(title)\n", + " plt.show()\n", + "\n", + "def plot_l1norm_histograms(bins, l1norm_histogram, title, xlabel, ylabel, xlim=None, log_scale=False):\n", + " \"\"\"\n", + " Plot L1-norm histograms for each scale.\n", + "\n", + " Parameters:\n", + " bins (list of numpy.ndarray): List of bin edges for each scale.\n", + " l1norm_histogram (list of numpy.ndarray): List of L1-norm histograms for each scale.\n", + " title (str): Title of the plot.\n", + " xlabel (str): Label for the x-axis.\n", + " ylabel (str): Label for the y-axis.\n", + " xlim (tuple or list): Limits for the x-axis (e.g., xlim=(0, 6)).\n", + " log_scale (bool): Whether to use a logarithmic scale for the y-axis.\n", + " \"\"\"\n", + " plt.figure()\n", + "\n", + " for scale, l1norm_hist in enumerate(l1norm_histogram):\n", + " plt.plot(bins[scale], l1norm_hist, label=f'Scale {scale}')\n", + "\n", + " plt.legend()\n", + " plt.xticks()\n", + " plt.yticks()\n", + " plt.xlabel(xlabel)\n", + " plt.ylabel(ylabel)\n", + " plt.grid(True)\n", + " plt.title(title)\n", + " \n", + " if xlim:\n", + " plt.xlim(xlim)\n", + " \n", + " if log_scale:\n", + " plt.yscale('log')\n", + " \n", + " plt.show()\n" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "def smooth_kappa_map(kappa_map, pixel_size_emap_amin):\n", - " # Set the standard deviation of the Gaussian filter based on the pixel size of the kappa map\n", - " precision_Peaks = 2 / pixel_size_emap_amin # pixel_size_emap_amin is the pixel size of the kappa map in arcminutes\n", - " kappa_map_smoothed = ndi.gaussian_filter(kappa_map, precision_Peaks)\n", - " return kappa_map_smoothed" + "Loading catalog" ] }, { @@ -283,16 +426,17 @@ "metadata": {}, "outputs": [], "source": [ - "def convert_to_snr_map(kappa_map_smoothed, noise_map_smoothed):\n", - " snr_map = kappa_map_smoothed / np.std(noise_map_smoothed)\n", - " return snr_map" + "# Load the catalog data as a numpy array\n", + "catalog_data = np.loadtxt(CATALOG_FILE)\n", + "catalog_data = catalog_data.T\n", + "ra, dec, g1_sim, g2_sim, kappa_sim = read_catalog_data(catalog_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Make SNR map from gaussian-smoothed noisy kappa-map" + "Create kappa map and make an SNR map from it" ] }, { @@ -301,34 +445,24 @@ "metadata": {}, "outputs": [], "source": [ - "# Read catalog data\n", - "ra, dec, g1_sim, g2_sim, kappa_sim = read_catalog_data(catalog_data)\n", - "\n", "# Create kappa map\n", - "kappaE, kappaB = create_kappa_map(ra, dec, g1_sim, g2_sim)\n", + "kappaE, _ = create_kappa_map(ra, dec, g1_sim, g2_sim)\n", "\n", "# Add noise to the kappa map\n", - "n_gal = 7 # galaxy number density in gal/arcmin^2\n", - "pix_arcmin = 0.4 # pixel size in arcmin\n", - "shape_noise = 0.44 # intrinsic ellipticity dispersion\n", - "kappaE_noisy, noise_map_CFIS_z05 = add_noise_to_kappa_map(kappaE, shape_noise, n_gal, pix_arcmin)\n", + "kappaE_noisy, noise_map_CFIS_z05 = add_noise_to_kappa_map(kappaE, SHAPE_NOISE, N_GAL, PIX_ARCMIN)\n", "\n", "# Smooth the noisy kappa map\n", - "kappaE_noisy_smoothed = smooth_kappa_map(kappaE_noisy, pix_arcmin)\n", + "kappaE_noisy_smoothed = smooth_kappa_map(kappaE_noisy, PIX_ARCMIN) \n", "\n", "# Compute SNR map\n", - "snr = convert_to_snr_map(kappaE_noisy_smoothed, kappaE_noisy_smoothed)\n" + "snr = convert_to_snr_map(kappaE_noisy_smoothed, kappaE_noisy_smoothed)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Plot:\n", - "1. Noiseless mass map\n", - "2. Noisy mass map\n", - "3. Noisy smoothed mass map\n", - "4. SNR map" + "Peak counts on the SNR map" ] }, { @@ -337,47 +471,15 @@ "metadata": {}, "outputs": [], "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "# Plot all four maps in subplots\n", - "fig, axs = plt.subplots(2, 2, figsize=(10, 7))\n", - "\n", - "# Plot the noiseless mass map\n", - "cmap0 = axs[0, 0].imshow(kappaE, cmap='inferno', vmin=-0.01, vmax=0.01, origin='lower')\n", - "axs[0, 0].set_title(\"Noiseless Mass Map\")\n", - "cbar0 = plt.colorbar(cmap0, ax=axs[0, 0])\n", - "\n", - "# Plot the noisy mass map\n", - "cmap1 = axs[0, 1].imshow(kappaE_noisy, cmap='inferno', vmin=-0.1, vmax=0.1, origin='lower')\n", - "axs[0, 1].set_title(\"Noisy Mass Map\")\n", - "cbar1 = plt.colorbar(cmap1, ax=axs[0, 1])\n", - "\n", - "# Plot the noisy smoothed mass map\n", - "cmap2 = axs[1, 0].imshow(kappaE_noisy_smoothed, cmap='inferno', vmin=-0.1, vmax=0.1, origin='lower')\n", - "axs[1, 0].set_title(\"Noisy Smoothed Mass Map\")\n", - "cbar2 = plt.colorbar(cmap2, ax=axs[1, 0])\n", - "\n", - "# Plot the SNR map\n", - "cmap3 = axs[1, 1].imshow(snr, cmap='inferno', origin='lower')\n", - "axs[1, 1].set_title(\"SNR Map\")\n", - "cbar3 = plt.colorbar(cmap3, ax=axs[1, 1])\n", - "\n", - "# Reversing x-axis for all plots\n", - "for ax in axs.flat:\n", - " ax.set_xlim(ax.get_xlim()[::-1])\n", - "\n", - "# Adjust layout\n", - "plt.tight_layout()\n", - "\n", - "# Show the plot\n", - "plt.show()" + "# Compute peak counts \n", + "kappa_th_center_snr, peak_counts_single = compute_single_scale_peak_counts(snr, KAPPA_SNR)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Make the peak-count histogram from the SNR map and plot it" + "Multiscale SNR maps" ] }, { @@ -386,41 +488,25 @@ "metadata": {}, "outputs": [], "source": [ - "import lenspack.peaks as peaks\n", - "\n", - "kappa_snr = np.linspace(-2, 6, 31) # SNR values to compute peak counts for\n", - "kappa_th_center_snr = 0.5 * (kappa_snr[:-1] + kappa_snr[1:]) # Center of SNR bins\n", - "\n", - "# Compute peak counts for the SNR map\n", - "kappa_snr = np.linspace(-2, 6, 31) # Adjust the range as needed\n", - "peak_counts_first = peaks.peaks_histogram(snr, kappa_snr)[0]\n", - "\n", - "# Plot peak counts for the first SNR map\n", - "plt.plot(kappa_th_center_snr, peak_counts_first)\n", - "# plt.yscale('log')\n", - "plt.xlabel('SNR smooth')\n", - "plt.ylabel('Peak Counts')\n", - "plt.title('Peak Counts Histogram for Gaussian SNR Map')\n", - "# plt.ylim(1e-4, 1e4)\n", - "plt.show()" + "# Compute multiscale SNR maps\n", + "multiscale_snr_maps = compute_multiscale_snr_maps(kappaE, noise_map_CFIS_z05, NSCALES)" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "Make multiscale SNR maps with starlet transform" + "# Compute peak counts for each scale\n", + "kappa_th_center_snr, peak_counts_multi = compute_multiscale_peak_counts(multiscale_snr_maps, KAPPA_SNR)" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "from lenspack.starlet_l1norm import noise_coeff, get_l1norm_noisy\n", - "from lenspack.image.transforms import starlet2d\n", - "from astropy.stats import mad_std" + "$\\ell_1$-norm" ] }, { @@ -429,53 +515,26 @@ "metadata": {}, "outputs": [], "source": [ - "def compute_multiscale_snr_maps(image, noise, nscales):\n", - " \"\"\"\n", - " Compute SNR maps for each wavelet scale of a noisy image.\n", - " \n", - " Parameters:\n", - " image (numpy.ndarray): The noiseless image.\n", - " noise (numpy.ndarray): The noise to be added to the image.\n", - " nscales (int): Number of wavelet scales for starlet decomposition.\n", - " \n", - " Returns:\n", - " snr_maps (list of numpy.ndarray): List of SNR maps for each scale.\n", - " \"\"\"\n", - " # Add noise to the noiseless image\n", - " image_noisy = image + noise\n", - " \n", - " # Perform starlet decomposition\n", - " image_starlet = starlet2d(image_noisy, nscales)\n", - " \n", - " # Estimate the noise level\n", - " noise_estimate = mad_std(image_noisy)\n", - " coeff_j = noise_coeff(image, nscales)\n", - " \n", - " snr_maps = []\n", - " for image_j, std_co in zip(image_starlet, coeff_j):\n", - " sigma_j = std_co * noise_estimate\n", - " \n", - " # Compute SNR map\n", - " snr_map = image_j / sigma_j\n", - " snr_maps.append(snr_map)\n", - " \n", - " return snr_maps" + "# Call the function to compute the L1-norm histogram\n", + "bins_l1, l1norm_histogram = get_l1norm_noisy(kappaE, noise_map_CFIS_z05, NSCALES, NBINS*2)" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "multiscale_snr_maps = compute_multiscale_snr_maps(kappaE, noise_map_CFIS_z05, nscales=5)" + "# Plotting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Plot the SNR maps of the 5 different scales and the coarse map" + "Plot:\n", + "1. Noiseless mass map\n", + "2. Noisy mass map\n", + "3. Noisy smoothed mass map\n", + "4. SNR map" ] }, { @@ -484,35 +543,24 @@ "metadata": {}, "outputs": [], "source": [ - "# Plot SNR maps in subplots with 2 rows\n", - "fig, axs = plt.subplots(2, 3, figsize=(15, 8))\n", - "\n", - "# Define vmin and vmax values for each SNR map\n", - "vmin_values = [-0.5, -3, -4, -4, -4, -4] # Replace with your desired vmin values\n", - "vmax_values = [0.5, 3, 4, 4, 4, 4] # Replace with your desired vmax values\n", + "# Define the map data and titles\n", + "map_data = [kappaE, kappaE_noisy, kappaE_noisy_smoothed, snr]\n", + "titles = ['Noiseless Mass Map', 'Noisy Mass Map', 'Noisy Smoothed Mass Map', 'SNR Map']\n", "\n", - "scales_arcmin = [2**(i+1) * pix_arcmin for i in range(len(multiscale_snr_maps))]\n", + "# Define the vmin and vmax values for each map (adjust these values as needed)\n", + "vmin_values = [-0.01, -0.1, -0.1, None]\n", + "vmax_values = [0.01, 0.1, 0.1, None]\n", "\n", - "for i, (snr_map, scale_arcmin) in enumerate(zip(multiscale_snr_maps, scales_arcmin)):\n", - " row = i // 3\n", - " col = i % 3\n", - " # cmap = axs[row, col].imshow(snr_map, cmap='inferno', origin='lower')\n", - " cmap = axs[row, col].imshow(snr_map, cmap='inferno', origin='lower', vmin = vmin_values[i], vmax = vmax_values[i])\n", - " axs[row, col].set_title(f\"SNR Map (Scale {scale_arcmin:.1f} arcmin)\")\n", - " cbar = plt.colorbar(cmap, ax=axs[row, col])\n", - "\n", - "# Adjust layout\n", - "plt.tight_layout()\n", - "\n", - "# Show the plot\n", - "plt.show()\n" + "# Plot each map using the plot_map function\n", + "for data, title, vmin, vmax in zip(map_data, titles, vmin_values, vmax_values):\n", + " plot_map(data, title=title, vmin=vmin, vmax=vmax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Calculate the peak counts for each scale" + "Plot the SNR maps of the 5 different scales and the coarse map" ] }, { @@ -521,30 +569,15 @@ "metadata": {}, "outputs": [], "source": [ - "def compute_peak_counts_on_snr_maps(snr_maps, kappa_snr):\n", - " \"\"\"\n", - " Compute peak counts for each wavelet scale of SNR maps.\n", - "\n", - " Parameters:\n", - " snr_maps (list of numpy.ndarray): List of SNR maps for each scale.\n", - " kappa_snr (numpy.ndarray): Array of kappa values corresponding to SNR maps.\n", + "# Define vmin and vmax values for each SNR map\n", + "vmin_values = [-0.5, -3, -4, -4, -4, -4] # Replace with your desired vmin values\n", + "vmax_values = [0.5, 3, 4, 4, 4, 4] # Replace with your desired vmax values\n", "\n", - " Returns:\n", - " kappa_th_center_snr (numpy.ndarray): Array of kappa threshold centers for peak counts.\n", - " peak_counts (list of numpy.ndarray): List of peak counts for each scale.\n", - " \"\"\"\n", - " kappa_th_center_snr = 0.5 * (kappa_snr[:-1] + kappa_snr[1:])\n", - " \n", - " peak_counts = [peaks.peaks_histogram(snr_map, kappa_snr)[0] for snr_map in snr_maps]\n", + "scales_arcmin = [2**(i+1) * PIX_ARCMIN for i in range(len(multiscale_snr_maps))]\n", "\n", - " return kappa_th_center_snr, peak_counts\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Plot the peak counts for the different scales" + "for i, (snr_map, scale_arcmin) in enumerate(zip(multiscale_snr_maps, scales_arcmin)):\n", + " title = f\"SNR Map (Scale {scale_arcmin:.1f} arcmin)\"\n", + " plot_map(snr_map, title=title, cmap='inferno', vmin=vmin_values[i], vmax=vmax_values[i])\n" ] }, { @@ -553,35 +586,20 @@ "metadata": {}, "outputs": [], "source": [ - "kappa_snr = np.linspace(-2, 6, 31) # SNR values to compute peak counts for\n", - "\n", - "# Compute peak counts for each scale\n", - "kappa_th_center_snr, peak_counts = compute_peak_counts_on_snr_maps(multiscale_snr_maps, kappa_snr)\n", - "\n", - "# Plot peak counts for each scale\n", - "plt.figure()\n", - "\n", - "for scale, peak_count in enumerate(peak_counts):\n", - " plt.plot(kappa_th_center_snr, peak_count, label=f'Scale {scale+1}')\n", - "\n", - "plt.legend()\n", - "plt.yscale('log')\n", - "plt.xticks()\n", - "plt.yticks()\n", - "plt.xlabel('SNR smooth')\n", - "plt.ylabel('Peak Counts')\n", - "plt.grid(True)\n", - "plt.xlim(-1, 6)\n", - "plt.ylim(1e-1, 1e5)\n", - "plt.title('Peak Counts for Different Scales')\n", - "plt.show()\n" + "# Plot single-scale peak count histogram\n", + "plot_peak_count_histograms(kappa_th_center_snr, peak_counts_single, \n", + " 'Peak Counts Histogram for Gaussian SNR Map', 'SNR smooth', 'Peak Counts')" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "Calculate the l1-norm and plot its histogram" + "# Plot multiscale peak count histograms\n", + "plot_peak_count_histograms(kappa_th_center_snr, peak_counts, \n", + " 'Peak Counts for Different Scales', 'SNR smooth', 'Peak Counts', log_scale=True)" ] }, { @@ -590,23 +608,8 @@ "metadata": {}, "outputs": [], "source": [ - "# Call the function to compute the L1-norm histogram\n", - "bins_l1, l1norm_histogram = get_l1norm_noisy(kappaE, noise_map_CFIS_z05, nscales=5, nbins=50)\n", - "\n", - "# Plot L1-norm histograms for each scale\n", - "for scale, l1norm_hist in enumerate(l1norm_histogram):\n", - " plt.plot(bins_l1[scale], l1norm_hist, label=f'scale {scale}')\n", - "\n", - "plt.legend()\n", - "# plt.yscale('log')\n", - "plt.xticks()\n", - "plt.yticks()\n", - "plt.xlabel('L1-norm')\n", - "plt.ylabel('Frequency')\n", - "plt.grid(True)\n", - "# plt.ylim(1e-4, 1e4)\n", - "plt.xlim(0, 6)\n", - "plt.show()\n" + "# Plot l1-norm histograms for different scales\n", + "plot_l1norm_histograms(bins_l1, l1norm_histogram, 'L1-norm Histograms for Different Scales', 'L1-norm', 'Frequency')" ] } ],