CMB_School_Part_12_v3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/jeffmcm1977/CMBAnalysis_SummerSchool/blob/master/CMB_School_Part_12_v3.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lensing reconstruction tutorial\n",
    "\n",
    "Primary Authors: Mathew Madhavacheril, Alexander van Engelen\n",
    "\n",
    "V2 additions: Jeff McMahon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python -c \"import cmb_modules\" || ( \\\n",
    "    wget https:/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 matplotlib.pyplot as plt\n",
    "# We use some stuff we learned before\n",
    "import cmb_modules as cm\n",
    "import lens_modules as lm\n",
    "np.random.seed(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will learn how to obtain the underlying lensing convergence -- the projected (dark) matter density -- from an observed map of the CMB. \n",
    "\n",
    "We will first simulate a lensed CMB map as we learnt in an earlier tutorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patch_deg_width = 40. # patch width in degrees\n",
    "pix_size = 1.5 # pixel size in arcminutes\n",
    "ells,ucltt,lcltt,clkk = lm.get_theory()\n",
    "N,lensed,kappa,ly,lx,modlmap = lm.get_lensed(patch_deg_width,pix_size,ells,ucltt,clkk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now convolve this map with a beam and add noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beam\n",
    "beam_arcmin = 1.4\n",
    "def gauss_beam(ell,fwhm):\n",
    "    # A gaussian beam transfer function (map-space, i.e. one power)\n",
    "    tht_fwhm = np.deg2rad(fwhm / 60.)\n",
    "    return np.exp(-(tht_fwhm**2.)*(ell**2.) / (16.*np.log(2.)))\n",
    "# Evaluate the beam on an isotropic Fourier grid\n",
    "kbeam2d = gauss_beam(modlmap,beam_arcmin)\n",
    "# Filter the map to convolve it with a beam\n",
    "beamed = lm.filter_map(lensed,kbeam2d)\n",
    "\n",
    "# Noise\n",
    "noise_uk_arcmin = 3. \n",
    "# White noise\n",
    "Clnoise = (noise_uk_arcmin*np.pi/180./60.)**2.\n",
    "Dlnoise = Clnoise*ells*(ells+1.)/2./np.pi\n",
    "# Make a GRF noise map\n",
    "noise_map = cm.make_CMB_T_map(N,pix_size,ells,Dlnoise)\n",
    "\n",
    "# The observed map\n",
    "observed = beamed + noise_map\n",
    "\n",
    "plt.title(\"simulated CMB map with lensing\")\n",
    "cm.Plot_CMB_Map(observed,-200,200,patch_deg_width,patch_deg_width)\n",
    "             \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lensing reconstruction\n",
    "\n",
    "Let's expand the lens equation\n",
    "\n",
    "$\\newcommand{\\al}{\\boldsymbol{\\alpha}}$\n",
    "$\\newcommand{\\x}{\\textbf{x}}$\n",
    "$\\newcommand{\\nab}{\\boldsymbol{\\nabla}}$\n",
    "$\\newcommand{\\uT}{{\\tilde{T}}}$\n",
    "$\\newcommand{\\uC}{{\\tilde{C}}}$\n",
    "$\\newcommand{\\dphi}{{\\nab\\phi}}$\n",
    "$\\newcommand{\\O}{{\\mathcal{O}}}$\n",
    "$\\newcommand{\\l}{{\\boldsymbol{\\ell}}}$\n",
    "$\\newcommand{\\L}{{\\textbf{L}}}$\n",
    "$$T(\\x) = \\uT(\\x+\\al)$$\n",
    "\n",
    "in a Taylor series (and substitute $\\al=\\dphi$) in the weak-lensing limit of small $\\phi$,\n",
    "\n",
    "$$\n",
    "T(\\x) \\approx \\uT(\\x) + \\nabla \\uT \\cdot \\dphi + \\O(\\dphi^2)\n",
    "$$ \n",
    "\n",
    " Homework exercise: Using 2d Fourier transforms,\n",
    "\n",
    "$$\n",
    "X(\\l) = \\int d^2\\x ~e^{-i\\l\\cdot\\x}X(\\x)\n",
    "$$\n",
    "$$\n",
    "X(\\x) = \\int \\frac{d^2\\l}{(2\\pi)^2} ~e^{i\\l\\cdot\\x}X(\\l)\n",
    "$$\n",
    "$$\n",
    "\\int d^2\\x e^{i\\l\\cdot\\x} = (2\\pi)^2 \\delta(\\l)\n",
    "$$\n",
    "\n",
    "show that to leading order in the lensing potential, lensing induces coupling between modes of the CMB map:\n",
    "\n",
    "\n",
    "$$\n",
    "\\langle T(\\l_1)T(\\l_2) \\rangle_{\\rm CMB} = \\phi(\\L)\\left[\\uC^{TT}_{\\ell_1}(\\L \\cdot \\l_1) + \\uC^{TT}_{\\ell_2}(\\L \\cdot \\l_2)\\right]\n",
    "$$\n",
    "\n",
    "where the averaging is over realizations of the unlensed CMB field $\\uT$.  Note: $L = \\ell_2 - \\ell_1$ represtnts the multiplole moment of the lensing reconstruciton.  \n",
    "\n",
    "The correlations present in the above result motivates extracting modes of the lensing potential $\\phi(\\L)$ by taking a weighted average of products of CMB modes. The weights can be derived such that the noise in the reconstruction is minimized. The resulting estimator is:\n",
    "\n",
    "$$\n",
    "\\hat{\\kappa}(\\L) = -A(\\L) \\int d^2\\l_1 \\left[\\uC^{TT}_{\\ell_1}(\\L \\cdot \\l_1) + \\uC^{TT}_{\\ell_2}(\\L \\cdot \\l_2)\\right]\\frac{T(\\l_1)T(\\l_2)}{C^{TT}_{\\ell_1}C^{TT}_{\\ell_2}}\n",
    "$$\n",
    "\n",
    "Homework exercise: Show that the above estimator can be written as a real-space multiplication\n",
    "\n",
    "$$\n",
    "\\hat{\\kappa} \\propto -\\nabla \\cdot \\left[\\left(\\nabla T_W\\right) T_H \\right]\n",
    "$$\n",
    "\n",
    "where $T_W$ is the temperature map Wiener filtered by $\\frac{\\uC^{TT}_{\\ell}}{C^{TT}_{\\ell}}$ and $T_H$ is the temperature map inverse variance filtered by $\\frac{1}{C^{TT}_{\\ell}}$.\n",
    "\n",
    "In what follows we use this real space expression to recontruct an estimate for the lensing map and explore some of the issues that must be taclked to extraxt science from these maps.   \n",
    "\n",
    "We refer you to the classic paper by [Hu and Okamoto](https://arxiv.org/pdf/astro-ph/0111606), and to a recent paper [Maniyar et al](https://arxiv.org/pdf/2101.12193) for more details of these quadratic estimators.  Maximum Liklihood techniques for lensign reconstruciton are know to produce higher signal to noise for the same experiment.  These are beyond the scope of this notebook.  If you are interested in details we recomend the classic paper by [Hirata and Seljak](https://arxiv.org/pdf/astro-ph/0209489) and the more recent paper by [Carron et al](https://arxiv.org/abs/1704.08230)\n",
    "\n",
    "### Exercise:\n",
    "Code up the real space estimate given above to complete the lensing estimator that has been started for you below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we attempt to reconstruct a lensing map from the observed map\n",
    "def qe_reconstruct(tmap,unlensed_cmb_power_2d,total_cmb_power_2d,ellmin,ellmax,modlmap,ly,lx):\n",
    "    \"\"\"\n",
    "    The simplest form of the quadratic estimator is\n",
    "\n",
    "    kappa_recon = - div ( grad(T) T )\n",
    "    where grad(T) is Wiener filtered with (unlensed_cmb_power/total_cmb_power)\n",
    "    and T is filtered with (1/total_cmb_power)\n",
    "    where the T map is beam deconvolved and\n",
    "    where total_cmb_power = lensed_cmb_power + noise_power / beam^2\n",
    "\n",
    "    The reconstruction is limited to ellmin < ells < ellmax of the CMB temperature map.\n",
    "\n",
    "    The grad and div operations require knowledge of:\n",
    "    modlmap: absolute wavenumbers of Fourier pixels\n",
    "    ly: map of ly Fourier coordinate of each pixel\n",
    "    lx: map of lx Fourier coordinate of each pixel\n",
    "    (note modlmap = sqrt(ly**2 + lx**2)\n",
    "    \"\"\"\n",
    "    inv_noise_filter = lm.kmask((1./total_cmb_power_2d),modlmap,ellmin,ellmax)\n",
    "    grad_filter = lm.kmask((unlensed_cmb_power_2d/total_cmb_power_2d),modlmap,ellmin,ellmax)\n",
    "\n",
    "    gradTy,gradTx = lm.gradient(tmap,ly,lx)\n",
    "\n",
    "    # The Wiener filtered gradient Grad(T_G)\n",
    "    filtered_gradTy = lm.filter_map(gradTy,grad_filter)\n",
    "    filtered_gradTx = lm.filter_map(gradTx,grad_filter)\n",
    "    # The inverse variance filtered T_H\n",
    "    filtered_T = lm.filter_map(tmap,inv_noise_filter)\n",
    "    \n",
    "    ############## COMPLETE THIS!!!\n",
    "    # ADD A LINE BELOW TO GET ukappa = divergence(Grad(T_G) * T_H)\n",
    "    \n",
    "    ############## COMPLETE THIS!!!\n",
    "    \n",
    "    # We filter with 1/L^2 at the end to make the normalization white\n",
    "    return -lm.filter_map(ukappa,lm.kmask(1/modlmap**2,modlmap,ellmin=2))\n",
    "\n",
    "# For this we also need to know how to calculate a divergence\n",
    "def div(imapy,imapx,ly,lx):\n",
    "    # Return divergence grad(Y)_y + grad(X)_x\n",
    "    ############## COMPLETE THIS!!!\n",
    "    # COMPLETE THIS FUNCTION THAT RETURNS THE DIVERGENCE\n",
    "    ############## COMPLETE THIS!!!\n",
    "    return divergence\n",
    "\n",
    "\n",
    "def interp(x,y,bounds_error=False,fill_value=0.,**kwargs):\n",
    "    # Just a simple interpolator that fills with zeros by default\n",
    "    from scipy.interpolate import interp1d\n",
    "    return interp1d(x,y,bounds_error=bounds_error,fill_value=fill_value,**kwargs)\n",
    "\n",
    "\n",
    "# To use this we need unlensed and total spectra interpolated on to the 2D Fourier grid\n",
    "### Prepare 2D Fourier space interpolations of theory spectra\n",
    "unlensed_cmb_power_2d = interp(ells,ucltt)(modlmap)\n",
    "total_cmb_power_2d = interp(ells,lcltt)(modlmap) + Clnoise/kbeam2d**2.\n",
    "\n",
    "\n",
    "### The noise was specified for a beam deconvolved map so we deconvolve the beam from our map\n",
    "tmap = lm.filter_map(observed,1/kbeam2d)\n",
    "\n",
    "### Next, we get the unnormalized reconstruction after choosing what multipoles to involve\n",
    "ellmin = 100\n",
    "ellmax = 3000\n",
    "ukappa = qe_reconstruct(tmap,unlensed_cmb_power_2d,total_cmb_power_2d,ellmin,ellmax,modlmap,ly,lx)\n",
    "# Calculating the normalization is a bit involved, so for now we focus on comparing\n",
    "# the largest scales of the reconstruction with that of the input, where the normalization is\n",
    "# expected to be somewhat constant\n",
    "\n",
    "### Filter the reconstruction and the input kappa to the same ellrange (where the normalization\n",
    "### is expected to be relatively white)\n",
    "kellmin = 10\n",
    "kellmax = 100\n",
    "ukappa_f = lm.filter_map(ukappa,lm.kmask(modlmap*0.+1.,modlmap,kellmin,kellmax))\n",
    "ikappa_f = lm.filter_map(kappa,lm.kmask(modlmap*0.+1.,modlmap,kellmin,kellmax))\n",
    "\n",
    "# Plot the filtered maps for comparison\n",
    "\n",
    "plt.title(\"filtered reconstructed kappa map\")\n",
    "cm.Plot_CMB_Lensing_Map(ukappa_f,patch_deg_width,patch_deg_width)\n",
    "             \n",
    "plt.title(\"filtered input kappa map\")\n",
    "cm.Plot_CMB_Lensing_Map(ikappa_f,patch_deg_width,patch_deg_width)\n",
    "         \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should be able to see a correspondence of peaks and valleys by eye in the above reconstruction if you have implemented it correctly.  Note- we have not computed the nomrmalization yet.   We will compare the input to recostructed map in two ways.  First using linear algebra to project out the common component, and secondly we will cross-correlate the reconstruction with the input and compare this to the input power spectrum. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## compare the input and reconstructed maps\n",
    "\n",
    "Since the normalization of the reconstructed lensing map is arbitrary, a simple difference between these maps will not provide a reasonable estimate of the residuals.  One must estimate the relative normalization to make this comparison.   Here we show a widely applicable trick from linear algera that is handy in a wide range of related circumstances.   \n",
    "\n",
    "If we treat the input map as a vector we can use methods from linear algebra to estimate the normalization and then make the comparison. This amounts to the following steps:\n",
    "\n",
    "1. Compute a unit vector aligned with the input map $\\vec I$: $ \\hat I = \\frac {\\vec I} {\\sqrt{|| \\vec I}||}$\n",
    "1. Compute the coefficent that fully removes the input map from the reconstructed map $\\eta= \\hat I \\cdot \\vec R$, where $\\vec R$ is the reconstructed map\n",
    "1. Use this coefficient to remove the input map from the reconstructed map to produce a noise map $\\vec N = \\vec R - \\eta \\hat I$\n",
    "\n",
    "Complete the code to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## remove the input map from the reconstructed map -- put your code here.  The first step is done for you.\n",
    "\n",
    "\n",
    "unit_i = ikappa_f / (np.sum(ikappa_f*ikappa_f))**.5\n",
    "eta = \n",
    "independent_part = \n",
    "\n",
    "plt.title(\"reconstructed minus input map\")\n",
    "cm.Plot_CMB_Lensing_Map(independent_part,patch_deg_width,patch_deg_width)\n",
    "         \n",
    "## "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise / Discussion Questions\n",
    "\n",
    "Discuss how well this method removed these signals.  \n",
    "1. How well did the reconstruction reproduce the input?\n",
    "1. What might cause these residuals?  How can you check?\n",
    "1. Did we make any choices in this analysis? If so, which could affect this?\n",
    "1. What real-wrold features of a true observing campaign could lead to additional biases?  If we have time, try implementing one of more of these to see what they do.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## cross-spectra\n",
    "\n",
    "EXERCISE: This looks good by eye! But are we sure it's not just a coincidence? The systematic way of checking that your lensing pipeline is working is to cross-correlate your reconstruction with the input kappa you provided. Calculate the binned 1d cross power and confirm it is non-zero and compare it to theory.  \n",
    "\n",
    "Below is code that can be addapted to carry out this cross-correlation.  It was pasted from the CMB power spectrum module. Complete the modificaitons needed to make this work. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Calculate binned 1d cross-power between reconstruction and input\n",
    "ell_max = 1000\n",
    "delta_ell = 10\n",
    "\n",
    "##\n",
    "Normalize_by_hand = 1\n",
    "\n",
    "## make a cross-power spectrum\n",
    "binned_ell, binned_spectrum = cm.calculate_2d_spectrum(,,delta_ell,ell_max,pix_size,N)\n",
    "#print binned_ell\n",
    "plt.semilogy(binned_ell,binned_spectrum* binned_ell * (binned_ell+1.)/2. / np.pi / Normalize_by_hand)\n",
    "plt.ylabel('$D_{\\ell}$ [$\\mu$K$^2$]')\n",
    "plt.xlabel('$L$')\n",
    "\n",
    "## reload the thoery\n",
    "ells,ucltt,lcltt,clkk = lm.get_theory()\n",
    "DlKK = clkk*ells*(ells+1.)/2./np.pi\n",
    "\n",
    "plt.semilogy(ells,DlKK)\n",
    "plt.xlim(0,ell_max)\n",
    "plt.ylim(1e-6,1)\n",
    "plt.show()\n",
    "\n",
    " \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise and Discussion\n",
    "\n",
    "1. adjust the nomarlization to make the curves line up.\n",
    "1. how well does the shape of the two curves match up?\n",
    "1. What expalins the abrupt cuttoff in the cross-correlation at $\\ell - 100$?\n",
    "1. Adjust the filtering so that the reconstruiton extends to $\\ell = 1000$. How well do these curves agree? You may have to adjust the normalization.  How do the maps look when you go up to this fine scale?\n",
    "1. Play wiht the high and low $\\ell$ limits for the reconstruciton and see hwo the angular scales and disucss what happens. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unforuntatly in real life, we do not know the true lensing signal to do this cross-correlation with an input and keep tuning it. However, we do understand the theory behind the quadratic estimator which provides an analytic expression for the response of a map to lensing.  This normalization is complicated, but possible to calculate analytically for the ideal case.  Small 10-20% corrections to this normalization can be estimated using well-understood simulations. We do not do that in this tutorial, but we may discuss it on the the third day of this summer school.\n",
    "\n",
    "We can compute cross-correlations between CMB lensing estimators and other tracers of large scale strucure like galaxy weak lensing, or the galaxy density estimated from visable or infrared photometric surveys.  However, these tracers include their own biases.   This provides measuermnts of a broad range of astrophysical paramters and cross-checks on cosmology and thus represnts a major research effort for the field. See the [DES 6x2 Analysis](https://arxiv.org/pdf/2206.10824) for example.\n",
    "\n",
    "\n",
    "## auto spectra \n",
    "Harvesting the full information content of CMB lensing requries dealing with these biases, but enables use of the CMB lensing autocorrelaiton.     In what follows we explore some of the issues that must be tackled to go beyond cross-correlation measumrnts.  \n",
    "\n",
    "### Exercise: \n",
    "  Modify the code you wrote above to compute the auto-correaltion.   Then ansewere the discusison quesitons below.  You should adjust the filting so you can see modes up to $L = 1000$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Normalize_by_hand = 2.5e6\n",
    "\n",
    "## make a auto-power spectrum\n",
    "binned_ell, binned_spectrum = cm.calculate_2d_spectrum(,,delta_ell,ell_max,pix_size,N)\n",
    "#print binned_ell\n",
    "plt.semilogy(binned_ell,binned_spectrum* binned_ell * (binned_ell+1.)/2. / np.pi / Normalize_by_hand**2)\n",
    "plt.ylabel('$D_{\\ell}$ [$\\mu$K$^2$]')\n",
    "plt.xlabel('$L$')\n",
    "\n",
    "## reload the thoery\n",
    "ells,ucltt,lcltt,clkk = lm.get_theory()\n",
    "DlKK = clkk*ells*(ells+1.)/2./np.pi\n",
    "\n",
    "plt.semilogy(ells,DlKK)\n",
    "plt.xlim(0,ell_max)\n",
    "plt.ylim(1e-6,1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise / Discussion questions\n",
    "\n",
    "1. Is there bias in this estimator?  How could you differentiate noise vs bias?\n",
    "1. What could affect the bias you see here?   Can you test abd confirm this?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bias In Lensing Estimators \n",
    "\n",
    "Dealing with bias in estimtes of the lensign power spectrum represts a large body of continuing research.   Major types of Biases include:\n",
    "\n",
    "1.   $N^{(0)}$ bias is the expectation value of the lensing power spectrum estimator that would be obtained if the CMB were an unlensed Gaussian field that had the power spectrum of the observed lensed CMB. \n",
    "1. $N^{(1)}$ and $N^{(2)}$ are smaller biases that arise from higher-order considerations\n",
    "2. Mean-field biases from statistical anisotropy due to a mask\n",
    "\n",
    "These effects can be mitigated through carful choice of estimators (bias hardened estimators, and cross-estimators) or removed through monte-carlo simulations.   These biases can be large, so dealing with these represntst the primary area of effort for prodcing cosmology from CMB lensing.  As the data go deeper, the required percisions for mitigating these effects must necissarly grow.   \n",
    "\n",
    "### Exercise:\n",
    "\n",
    "Reproduce the existence of the mean-field bias by applying apodization to the edge of the map before reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "## make a window\n",
    "window = cm.cosine_window(N)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}