CMB_School_Part_05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/jeffmcm1977/CMBAnalysis_SummerSchool/blob/master/CMB_School_Part_05.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matched Filter Techniques\n",
    "\n",
    "### Jeff McMahon and Renée Hložek "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have seen how to extract an estimate for the power spectrum from a CMB map.  In this exercise we extract the point sources from our simulated maps.  We do this using the matched filter technqiue. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First we read in the necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python -c \"import cmb_modules\" || ( \\\n",
    "    wget https://github.com/jeffmcm1977/CMBAnalysis_SummerSchool/raw/master/cmb_school.tar.gz && \\\n",
    "    tar xzvf cmb_school.tar.gz \\\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import sys\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.mlab as mlab\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.io.fits as fits\n",
    "\n",
    "import constants as cs  # the constants module\n",
    "import cmb_modules  # the module of functions\n",
    "\n",
    "N = cs.N\n",
    "N_iterations = cs.N_iterations\n",
    "c_min = cs.c_min\n",
    "c_max = cs.c_max\n",
    "X_width = cs.X_width\n",
    "Y_width = cs.Y_width\n",
    "beam_size_fwhp = cs.beam_size_fwhp\n",
    "\n",
    "pix_size = cs.pix_size\n",
    "\n",
    "Number_of_Sources = cs.Number_of_Sources\n",
    "Amplitude_of_Sources = cs.Amplitude_of_Sources\n",
    "Number_of_Sources_EX = cs.Number_of_Sources_EX\n",
    "Amplitude_of_Sources_EX = cs.Amplitude_of_Sources_EX\n",
    "\n",
    "Number_of_SZ_Clusters = cs.Number_of_SZ_Clusters\n",
    "Mean_Amplitude_of_SZ_Clusters = cs.Mean_Amplitude_of_SZ_Clusters\n",
    "SZ_beta = cs.SZ_beta\n",
    "SZ_Theta_core = cs.SZ_Theta_core\n",
    "\n",
    "white_noise_level = cs.white_noise_level\n",
    "atmospheric_noise_level = 0  # cs.atmospheric_noise_level\n",
    "one_over_f_noise_level = 0  # cs.one_over_f_noise_level"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rebuilding the simulated map from previous exercises"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Make a CMB map\n",
    "ell, DlTT = np.loadtxt(\"CAMB_fiducial_cosmo_scalCls.dat\", usecols=(0, 1), unpack=True)\n",
    "CMB_T = cmb_modules.make_CMB_T_map(N, pix_size, ell, DlTT)\n",
    "\n",
    "N = int(N)\n",
    "## make a point source map\n",
    "PSMap = cmb_modules.Poisson_source_component(\n",
    "    N, pix_size, Number_of_Sources, Amplitude_of_Sources\n",
    ")\n",
    "PSMap += cmb_modules.Exponential_source_component(\n",
    "    N, pix_size, Number_of_Sources_EX, Amplitude_of_Sources_EX\n",
    ")\n",
    "\n",
    "## make an SZ map\n",
    "SZMap, SZCat = cmb_modules.SZ_source_component(\n",
    "    N,\n",
    "    pix_size,\n",
    "    Number_of_SZ_Clusters,\n",
    "    Mean_Amplitude_of_SZ_Clusters,\n",
    "    SZ_beta,\n",
    "    SZ_Theta_core,\n",
    "    False,\n",
    ")\n",
    "\n",
    "## add them all together to get the sky map at a single freuqency\n",
    "total_map = CMB_T + PSMap + SZMap\n",
    "\n",
    "## incorporate the impact of the instrument\n",
    "## beam\n",
    "CMB_T_convolved = cmb_modules.convolve_map_with_gaussian_beam(\n",
    "    N, pix_size, beam_size_fwhp, total_map\n",
    ")\n",
    "## noise\n",
    "Noise = cmb_modules.make_noise_map(\n",
    "    N, pix_size, white_noise_level, atmospheric_noise_level, one_over_f_noise_level\n",
    ")\n",
    "\n",
    "total_map_plus_noise = CMB_T_convolved + Noise\n",
    "\n",
    "total_map_plus_noise_original = np.copy(total_map_plus_noise)\n",
    "SZCat_original = np.copy(SZCat)\n",
    "## plot the result\n",
    "p = cmb_modules.Plot_CMB_Map(total_map_plus_noise, c_min, c_max, X_width, Y_width)\n",
    "\n",
    "## we will need a window funciton below, so we creat that here\n",
    "window = cmb_modules.cosine_window(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few point sources and SZ sources are visible in our simulated map, however there are many that are faint and difficult to identify.  There are sophisticated techniques for isolating these sources. We now discuss the matched filter approach.\n",
    "\n",
    "### Constructing a matcher filter\n",
    "\n",
    "We construct a single frequency matched filter following Melin etc...  the punch line is that we build the following filter in Fourier space:\n",
    "\n",
    "$$\\psi = \\frac{B_\\ell S_\\ell }{B_\\ell^2 N^2_{ap,\\ell} +N^2_{ins,\\ell} } $$\n",
    "\n",
    "Both $\\psi$ and all the other quantities are two dimensional in the sense that they depend on $k_x$ and $k_y$. We refer to these coordinates with a $\\ell$ as shorthand. In this equation $B_\\ell$ is the Fourier transform of the beam pattern and filtering; $S_\\ell$ is the Fourier transform of the signal template; $N^2_{ap,\\ell}$ is the power spectrum (e.g., absolute value squared of the 2d FFT) of the astrophysical noise from sources other than what is desired to be found, and $N^2_{ins,\\ell}$ is the power spectrum of the instrumental response.\n",
    "\n",
    "\n",
    "In the next box we define a function to calculate the matched filter and a function to plot the resulting map. In subsequent boxes we do the work of estimating the terms in the filter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matched_filter(input_map, beam_and_filt, signal_profile, FT_noise_covar):\n",
    "    ## input_map: the map we are processing\n",
    "    ## beam_and_filt: the beam convolved with any map filtering, in real space\n",
    "    ## signal_profile: the shape of the signal we are looking for, in real spcae\n",
    "    ## FT_noise_covar: the B_N_{ap}^2 + N_{ins}^2 in Fourier space\n",
    "    ## calculating FT_npoise_covar is expensive so it is done externally\n",
    "\n",
    "    FT_beam_and_filt = np.fft.fft2(\n",
    "        np.fft.fftshift(beam_and_filt)\n",
    "    )  ## transform beam_and_filt to Fourier space\n",
    "    FT_signal = np.fft.fft2(\n",
    "        np.fft.fftshift(signal_profile)\n",
    "    )  ## transform cluster_profile to Fourier space\n",
    "\n",
    "    psi = (\n",
    "        FT_beam_and_filt * FT_signal / FT_noise_covar\n",
    "    )  ## define the matchedfilter function\n",
    "\n",
    "    filtered = psi * np.fft.fft2(np.fft.fftshift(input_map))  ## filter the map\n",
    "    filtered = np.fft.fftshift(np.fft.ifft2(filtered))  ## center the filter\n",
    "    filtered = np.real(filtered)  ## change the data type to real\n",
    "    return filtered\n",
    "\n",
    "\n",
    "def Plot_Matched_Filtered_Map(Map_to_Plot, X_width, Y_width):\n",
    "    from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "    print(\"map mean:\", np.mean(Map_to_Plot), \"map rms:\", np.std(Map_to_Plot))\n",
    "    plt.figure(figsize=[10, 10])\n",
    "    im = plt.imshow(\n",
    "        Map_to_Plot, interpolation=\"bilinear\", origin=\"lower\", cmap=cm.RdBu_r\n",
    "    )\n",
    "    # im.set_clim(c_min,c_max)\n",
    "    ax = plt.gca()\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "    cbar = plt.colorbar(im, cax=cax)\n",
    "    im.set_extent([0, X_width, 0, Y_width])\n",
    "    plt.ylabel(r\"angle $[^\\circ]$\")\n",
    "    plt.xlabel(r\"angle $[^\\circ]$\")\n",
    "    cbar.set_label(r\"matched_filter [S/N]\", rotation=270)\n",
    "    plt.show()\n",
    "    return 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following what we did above to compute the power spectrum, we will calculate the noice covariance (the denominator in the matched filter definition) using Monte Carlo techniques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## construct the 2d noise covariance in Fourier space\n",
    "FT_noise_covar = np.zeros((N, N))  ## a 2d array to hold the result\n",
    "\n",
    "N_iterations = 16\n",
    "\n",
    "## make a series of simulated maps, find the power spectrum, and average these to estimate the noise covariance\n",
    "i = 0\n",
    "while i < N_iterations:\n",
    "    ## simulate the astrophysical map\n",
    "    CMB_T = cmb_modules.make_CMB_T_map(N, pix_size, ell, DlTT)\n",
    "    PSMap = cmb_modules.Poisson_source_component(\n",
    "        N, pix_size, Number_of_Sources, Amplitude_of_Sources\n",
    "    )\n",
    "    PSMap += cmb_modules.Exponential_source_component(\n",
    "        N, pix_size, Number_of_Sources_EX, Amplitude_of_Sources_EX\n",
    "    )\n",
    "    SZMap, trash = cmb_modules.SZ_source_component(\n",
    "        N,\n",
    "        pix_size,\n",
    "        Number_of_SZ_Clusters,\n",
    "        Mean_Amplitude_of_SZ_Clusters,\n",
    "        SZ_beta,\n",
    "        SZ_Theta_core,\n",
    "        False,\n",
    "    )\n",
    "    CMB_T = CMB_T + PSMap + SZMap  ## the astrophysical map\n",
    "\n",
    "    ## fold in the instrument response\n",
    "    CMB_T_convolved = cmb_modules.convolve_map_with_gaussian_beam(\n",
    "        N, pix_size, beam_size_fwhp, CMB_T\n",
    "    )\n",
    "    Noise = cmb_modules.make_noise_map(\n",
    "        N, pix_size, white_noise_level, atmospheric_noise_level, one_over_f_noise_level\n",
    "    )\n",
    "\n",
    "    ## Fourier transform the map\n",
    "    temp = np.fft.fft2(\n",
    "        np.fft.fftshift(window * (CMB_T_convolved + Noise))\n",
    "    )  ## these are the two terms in the denominator\n",
    "\n",
    "    ## now average\n",
    "    FT_noise_covar += np.real(np.conj(temp) * temp / (N_iterations * 1.0))\n",
    "    ## note the progress\n",
    "    sys.stdout.write(\n",
    "        \"\\r matched filter noise realization, iterations complete: {} of {}\".format(\n",
    "            (i + 1), N_iterations\n",
    "        )\n",
    "    )\n",
    "    sys.stdout.flush()\n",
    "    ## iterate\n",
    "    i = i + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the noise covariance and a function to compute the matched filter, we filter our map. To interpret this map we will divide by the standard deviation of the filtered map to put our map in signal to noise units.   \n",
    "\n",
    "You can play around with the number of iterations you use in the Monte Carlo of the noise term. If it is too small, you will get high frequency ringing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## construct the beam and cluster profile for the numerator of the matched filter\n",
    "beam_and_filt = cmb_modules.make_2d_gaussian_beam(\n",
    "    N, pix_size, beam_size_fwhp\n",
    ")  ## this is the filtering we did on the map\n",
    "cluster_profile = cmb_modules.beta_function(\n",
    "    N, pix_size, SZ_beta, SZ_Theta_core\n",
    ")  ## this is the signal we are looking for\n",
    "\n",
    "## Apply the matched filter to our map\n",
    "filtered_map = matched_filter(\n",
    "    total_map_plus_noise_original, beam_and_filt, cluster_profile, FT_noise_covar\n",
    ")\n",
    "\n",
    "## make a S/N map\n",
    "SN_map = filtered_map / np.std(filtered_map)\n",
    "\n",
    "## make a few plots\n",
    "p = cmb_modules.Plot_CMB_Map(total_map_plus_noise, c_min, c_max, X_width, Y_width)\n",
    "p = Plot_Matched_Filtered_Map(SN_map, X_width, Y_width)\n",
    "\n",
    "cluster_hist, cluster_bin_edges = np.histogram(\n",
    "    SN_map, bins=100, range=[SN_map.min(), SN_map.max()]\n",
    ")\n",
    "plt.semilogy(cluster_bin_edges[0:-1], cluster_hist)\n",
    "plt.ylabel(\"number of pixels\")\n",
    "plt.xlabel(\"matched filter signal to noise\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we show three maps: (top) the raw simulated map, (middle) the matched filtered map, and (bottom) a histogram of the S/N of each pixel in the matched filter map. You could identify SZ sources by doing a S/N cut and taking all the pixels below say -5. To find point sources we would replace the beta function in the definition of the filter with a delta function and repeat the exercise looking for S/N greater than say 5.\n",
    "\n",
    "NOTE: The signal to noise cut at $5 \\sigma$ is not sacred. It is chosen so that the probability of a spurious cluster (or point source) detection is controlled to be irrelevant on a ~100 degree patch with arcminute resolution. Depending on the application, this cut should be adjusted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  Re-run this code with different matched filter cluster radii on your map from Notebook 02 that has multiple cluster radii and see how the S/N varies for different clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code and plots go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "comment on your results here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  Make a matched filter for point sources. Compare the S/N for SZ clusters and for point sources between this filtered map and the SZ filtered map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code and plots go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "comment on your results here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  Extract a source list by (1) finding the highest S/N cluster detection, (2) set the pixels to zero within some radius of that cluster, and (3) iterate until you reach some S/N threshold.   Compare the source list you extract with the SZCat catalogue of true sources in the simulation.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code and plots go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "comment on your results here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stacking\n",
    "\n",
    "One often wants to understand signals at low mass and therefore the low signal to noise end of things. If an external catalogue exists, one can co-add the signals from objects in the external catalogue to boost the signal to noise. In this example we stack the simulated observed map and the matched filtered map at the positions of the clusters we put into the simulation. This highlights the opportunities that exist with stacking.  \n",
    "\n",
    "Interpreting these results requires careful analysis. Monte Carlo simulations including selection effects and systematic effects are a particularly appealing approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## take SZCat and stack total_map_plus_noise on the SZ positions, do this in a mass bin\n",
    "\n",
    "\n",
    "def Stack_on_Positions(total_map, N, cat, N_objects, bin_min, bin_max, Radius):\n",
    "    Radius = int(Radius)\n",
    "    stack = np.zeros([Radius * 2, Radius * 2])\n",
    "    counter = 0\n",
    "    i = 0\n",
    "    while i < N_objects:\n",
    "        ampl = cat[2, i]\n",
    "        if (ampl > bin_min) and (ampl <= bin_max):\n",
    "            xc = cat[0, i]\n",
    "            yc = cat[1, i]\n",
    "            if (xc > Radius) and (xc < N - Radius):\n",
    "                if (yc > Radius) and (yc < N - Radius):\n",
    "\n",
    "                    stack += total_map[\n",
    "                        int(xc - Radius) : int(xc + Radius),\n",
    "                        int(yc - Radius) : int(yc + Radius),\n",
    "                    ]\n",
    "                    counter += 1\n",
    "        i = i + 1\n",
    "    return stack / counter\n",
    "\n",
    "\n",
    "stack = Stack_on_Positions(\n",
    "    total_map_plus_noise, N, SZCat, Number_of_SZ_Clusters, -100000, 100000, 50\n",
    ")\n",
    "stack_SN = Stack_on_Positions(\n",
    "    SN_map, N, SZCat, Number_of_SZ_Clusters, -100000, 100000, 50\n",
    ")\n",
    "\n",
    "p = cmb_modules.Plot_CMB_Map(\n",
    "    stack, c_min / 4.0, c_max / 4.0, X_width * 50 * 2 / N, Y_width * 50 * 2 / N\n",
    ")\n",
    "p2 = Plot_Matched_Filtered_Map(stack_SN, X_width * 50 * 2 / N, Y_width * 50 * 2 / N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The upper plot shows the CMB temperature map stacked on the locations of SZ clusters.  The lower plot shows the same stacking exercise repeated with the matched filtered map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  In this exercise you will calculate the SZ signal vs. input amplitude below the S/N 5 threshold in the maps. This is similar to calculating an SZ scaling relation from optical or x-ray selected clusters.\n",
    "\n",
    "To do this, divide the cluster into \"SZ-amplitude\" bins (the third column of the SZCat), calculate the peak of the stacked map (CMB and S/N) for this amplitude bin. Find a method for estimating the noise of this stack. Plot the recovered SZ signal vs. the true amplitude.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code and plots go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "comment on your results here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### centering errors\n",
    "\n",
    "In this example we introduce errors in the positions of the clusters in the catalogue.  It is obvious that the cluster signal broadens and washes out. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "centering_errors_x = np.random.normal(0, 2, Number_of_SZ_Clusters)\n",
    "centering_errors_y = np.random.normal(0, 2, Number_of_SZ_Clusters)\n",
    "SZCat_centering_errs = SZCat\n",
    "SZCat_centering_errs[0, :] += centering_errors_x\n",
    "SZCat_centering_errs[1, :] += centering_errors_y\n",
    "\n",
    "stack = Stack_on_Positions(\n",
    "    total_map_plus_noise,\n",
    "    N,\n",
    "    SZCat_centering_errs,\n",
    "    Number_of_SZ_Clusters,\n",
    "    -100000,\n",
    "    100000,\n",
    "    50,\n",
    ")\n",
    "p = cmb_modules.Plot_CMB_Map(\n",
    "    stack, c_min / 4.0, c_max / 4.0, X_width * 50 * 2 / N, Y_width * 50 * 2 / N\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  Repeat exercise 18 with this centering error turned on. Compare the results. This is a prime example of systematics that must be accounted for in extracting SZ scaling relations through stacking."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code and plots go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "comment on your results here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## cross-correlation of optical and CMB maps\n",
    "\n",
    "An alternative to stacking is cross correlation. For example, one can take an optical catalogue, make a map in the same pixelization as the CMB map or the density of sources (or sum of the masses or other property), smooth that map and take the power spectrum using the techniques we used for CMB analysis.\n",
    "\n",
    "<font color='red'>EXERCISE: </font>  Do this.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "comment on your results here"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}