CosmoAndes_PartFour.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Map Analysis\n",
    "## CosmoAndes 2018\n",
    "### Adapted from notebook developed by\n",
    "### Jeff McMahon and Renée Hložek for the 2016 CMB Data Analysis Summer School at the University of Michigan\n",
    "\n",
    "We will analyze these maps in several ways. First, in this part, we will use a Monte Carlo technique to recover the power spectrum, with exercises along the way. In part five we will, use a matched filter to isolate the SZ signal and point sources, and work through stacking a map on a cluster catalogue from a simulated external survey (e.g., optical or x-ray).\n",
    "\n",
    "## Power spectrum\n",
    " \n",
    "In this section we compute the power spectrum of a CMB map.  We will work in the flat sky approximation which is described above in the section where we generated the simulated CMB map.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\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",
    "%matplotlib inline\n",
    "import constants as cs # the constants module\n",
    "\n",
    "from cmb_modules import * # the module of functions\n",
    "\n",
    "N = cs.N\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 = cs.atmospheric_noise_level\n",
    "one_over_f_noise_level = cs.one_over_f_noise_level"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recalculating the results from the previous stages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.  0.  0. ...,  0.  0.  0.]\n",
      " [ 0.  0.  0. ...,  0.  0.  0.]\n",
      " [ 0.  0.  0. ...,  0.  0.  0.]\n",
      " ..., \n",
      " [ 0.  0.  0. ...,  0.  0.  0.]\n",
      " [ 0.  0.  0. ...,  0.  0.  0.]\n",
      " [ 0.  0.  0. ...,  0.  0.  0.]]\n",
      "('Number of sources required: ', 5000)\n",
      "393\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-8cf79cbaf87f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;31m## make an SZ map\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mSZMap\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mSZCat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSZ_source_component\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpix_size\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mNumber_of_SZ_Clusters\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mMean_Amplitude_of_SZ_Clusters\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mSZ_beta\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mSZ_Theta_core\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;31m## add them all together to get the sky map at a single freuqency\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/reneehlozek/Code/CosmoAndes/cmb_modules.py\u001b[0m in \u001b[0;36mSZ_source_component\u001b[0;34m(N, pix_size, Number_of_SZ_Clusters, Mean_Amplitude_of_SZ_Clusters, SZ_beta, SZ_Theta_core, do_plots)\u001b[0m\n\u001b[1;32m    121\u001b[0m         \u001b[0mpix_y\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    122\u001b[0m         \u001b[0mpix_amplitude\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexponential\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mMean_Amplitude_of_SZ_Clusters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1.\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 123\u001b[0;31m         \u001b[0mSZcat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpix_x\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    124\u001b[0m         \u001b[0mSZcat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpix_y\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    125\u001b[0m         \u001b[0mSZcat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpix_amplitude\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices"
     ]
    }
   ],
   "source": [
    "## Make a CMB map\n",
    "ell, DlTT = np.loadtxt(\"CAMB_fiducial_cosmo_scalCls.dat\", usecols=(0, 1), unpack=True) \n",
    "CMB_T = make_CMB_T_map(N,pix_size,ell,DlTT)\n",
    "\n",
    "## make a point source map\n",
    "PSMap = Poisson_source_component(N,pix_size,Number_of_Sources,Amplitude_of_Sources) \n",
    "PSMap +=  Exponential_source_component(N,pix_size,Number_of_Sources_EX,Amplitude_of_Sources_EX)\n",
    "\n",
    "## make an SZ map\n",
    "SZMap,SZCat = SZ_source_component(N,pix_size,Number_of_SZ_Clusters,Mean_Amplitude_of_SZ_Clusters,SZ_beta,SZ_Theta_core,False)\n",
    "\n",
    "## add them all together to get the sky map at a single freuqency\n",
    "total_map = CMB_T + PSMap + SZMap\n",
    "\n",
    "## incorperate the impact of the instrument\n",
    "    ## beam\n",
    "CMB_T_convolved =convolve_map_with_gaussian_beam(N,pix_size,beam_size_fwhp,total_map)\n",
    "    ## noise\n",
    "Noise = make_noise_map(N,pix_size,white_noise_level,atmospheric_noise_level,one_over_f_noise_level)\n",
    "\n",
    "total_map_plus_noise = CMB_T_convolved + Noise\n",
    "\n",
    "## plot the result\n",
    "p=Plot_CMB_Map(total_map_plus_noise,c_min,c_max,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Apodize  the map to eliminate edge effects\n",
    "\n",
    "Before taking a 2d FFT (eg the obvious thing to do for computing a power spectrum) we must apodize the maps to eliminate edges effects.  Edge effects come about because the Fourier transform treats a square array as having periodic boundaries.  Thus if we take the Fourier transform of a 2-dimensional map and the values on the left and right side (and also, top and bottom) of the map don't match, we end up generating spurious signals.   In this example we use a cosine window to smoothly roll off the signal to zero as we approach the edges of the map.  There are many choices of windows that trade sensitivity loss, coupling of adjacent modes, and ringing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def cosine_window(N):\n",
    "    \"makes a cosine window for apodizing to avoid edges effects in the 2d FFT\" \n",
    "    # make a 2d coordinate system\n",
    "    N=int(N) \n",
    "    ones = np.ones(N)\n",
    "    inds  = (np.arange(N)+.5 - N/2.)/N *np.pi ## eg runs from -pi/2 to pi/2\n",
    "    X = np.outer(ones,inds)\n",
    "    Y = np.transpose(X)\n",
    "  \n",
    "    # make a window map\n",
    "    window_map = np.cos(X) * np.cos(Y)\n",
    "   \n",
    "    # return the window map\n",
    "    return(window_map)\n",
    "  ###############################\n",
    "    \n",
    "window = (cosine_window(N))\n",
    "    \n",
    "appodized_map = window * total_map_plus_noise\n",
    "\n",
    "p=Plot_CMB_Map(appodized_map,c_min,c_max,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows our simulated map with a cosine window applied to eliminate edge effects.  It is obvious from this map that we are suppressing the signal here.\n",
    "<font color='red'>EXCERCISE: </font>  There are an huge number of well studied windows with various combinations of properties.   Some minimize mode coupling, others minimize signal loss, while others maximize some combination of the two.  Find the wikipedia article on Fourier transform windows, choose one of your favorites and implement it as an option.  Compare the impact of the new window compared to the simple cosine window on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## your code and plots go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your comments go here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Naive Powerspectrum\n",
    "\n",
    "Here we compute a naive power spectrum and compare it to the input power spectrum for our simulations.  The power spectrum is computed by: (1) applying a 2d FFT, (2) taking the absolute value squared of this map in Fourier space ($k_x$ and $k_y$), and (3) averaging the signal in annular bins of $k = \\sqrt{k_x^2 + k_y^2}$.  These bins are converted to bins in $\\ell$ with the scaling: $\\ell = k* 2 \\pi$ per the flat sky approximation.   NOTE: step 3 (averaging in radial bins) is how we convert our 2d Fourier map into a 1d power spectrum.\n",
    "\n",
    "Our spectrum code takes two maps as inputs to allow for cross spectra when we get to polarizaiton in Part Seven."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#### parameters for setting up the spectrum\n",
    "delta_ell = 50.\n",
    "ell_max = 5000.\n",
    "\n",
    "if max(ell)< ell_max: \n",
    "        print 'WARNING: Your theory curves end before the binned ell_max'\n",
    "\n",
    "def calculate_2d_spectrum(Map1,Map2,delta_ell,ell_max,pix_size,N):\n",
    "    \"calcualtes the power spectrum of a 2d map by FFTing, squaring, and azimuthally averaging\"\n",
    "    N=int(N)\n",
    "    # make a 2d ell coordinate system\n",
    "    ones = np.ones(N)\n",
    "    inds  = (np.arange(N)+.5 - N/2.) /(N-1.)\n",
    "    kX = np.outer(ones,inds) / (pix_size/60. * np.pi/180.)\n",
    "    kY = np.transpose(kX)\n",
    "    K = np.sqrt(kX**2. + kY**2.)\n",
    "    ell_scale_factor = 2. * np.pi \n",
    "    ell2d = K * ell_scale_factor\n",
    "    \n",
    "    # make an array to hold the power spectrum results\n",
    "    N_bins = int(ell_max/delta_ell)\n",
    "    ell_array = np.arange(N_bins)\n",
    "    CL_array = np.zeros(N_bins)\n",
    "    \n",
    "    # get the 2d fourier transform of the map\n",
    "    FMap1 = np.fft.ifft2(np.fft.fftshift(Map1))\n",
    "    FMap2 = np.fft.ifft2(np.fft.fftshift(Map2))\n",
    "    PSMap = np.fft.fftshift(np.real(np.conj(FMap1) * FMap2))\n",
    "    # fill out the spectra\n",
    "    i = 0\n",
    "    while (i < N_bins):\n",
    "        ell_array[i] = (i + 0.5) * delta_ell\n",
    "        inds_in_bin = ((ell2d >= (i* delta_ell)) * (ell2d < ((i+1)* delta_ell))).nonzero()\n",
    "        CL_array[i] = np.mean(PSMap[inds_in_bin])\n",
    "        #print i, ell_array[i], inds_in_bin, CL_array[i]\n",
    "        i = i + 1\n",
    " \n",
    "    # return the power spectrum and ell bins\n",
    "    return(ell_array,CL_array*np.sqrt(pix_size /60.* np.pi/180.)*2.)\n",
    "\n",
    "## make a power spectrum\n",
    "binned_ell, binned_spectrum = calculate_2d_spectrum(appodized_map,appodized_map,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)\n",
    "plt.semilogy(ell,DlTT)\n",
    "plt.ylabel('$D_{\\ell}$ [$\\mu$K$^2$]')\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot shows the input CMB power spectrum (green) and the naive power spectrum we estimated from our CMB map (blue).  The naive spectrum does not match the input due to a combination of: instrumental noise, SZ and point source signals; and suppression from the beam and the apodization.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXCERCISE: </font>  Comment on how the measured power sepctrum (blue) compares to the input spectrum (green).  Why are they different?  What features are common to both?  How does your alternative window funciton affect the measured spectrum compared to the default?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## your code and plots go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your comments go here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correcting the biases in the naive power spectrum\n",
    "\n",
    "To correct the naive spectrum and obtain an unbiased estimate of the underlying power spectrum in our simulated map we must correct additive bias and multiplicative bias.   We can relate our measured spectrum $\\hat D_\\ell$ to the true underlying spectrum $D_\\ell$ as follows:\n",
    "\n",
    "$$\\hat D_\\ell = T*D_\\ell + N. $$\n",
    "\n",
    "Here $N$ represents an additive noise term and $T$ represents a transfer function of the instrument (beam) and filtering (in this case the apodization, but other processing can enter, for example a 1-d high pass filter to eliminate striping).   All of these terms are functions of $\\ell$.\n",
    "\n",
    "The true power spectrum can be recovered by using monte carlo technqies.  To apply this technique we use simulations to calibrate $T$ and $N$ and then use algebra to correct the naive measurement.  This is compactly expressed in the following equation:\n",
    "\n",
    "$$ D_\\ell = \\hat (D_\\ell - N)/T. $$\n",
    "\n",
    "In the next two sections we use simulations to calibrate $T and $N$ to recover an unbiased estimate of the underlying power spectrum.  Subsequently we will used monte carlo simulations to estimate the error bars on this measurement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calibrating the transfer function\n",
    "\n",
    "The transfer function can be calibrated by: (1) generating sky simulations with a known power spectrum, modeling the transfer function of the instrument and the post-processing, and keeping the noise level to zero,; and (2) calculating the naive power spectrum from each simulation, (3) running many simulations to reduce numerical noise, and (4) dividing the true spectrum by the the average signal only spectrum to recover our estimate for the transfer function.   The accuracy depends on the number of realizations used.  Here  we used 64 realizations as a compromise between speed and accuracy.  \n",
    "\n",
    "Here we use a CMB only spectrum to estimate the transfer function.  We could improve the estimate of the transfer function by using an input power spectrum that is modified to follow the effect of point sources.  This could be added as an exercise.  In this exercise compare the transfer function to find out how much the choice of input spectrum matters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N_iterations = 16\n",
    "\n",
    "signal_only  = np.zeros([N_iterations,int(ell_max/delta_ell)])\n",
    "i = 0\n",
    "while (i <N_iterations):\n",
    "    CMB_T = make_CMB_T_map(N,pix_size,ell,DlTT)\n",
    "    CMB_T_convolved = convolve_map_with_gaussian_beam(N,pix_size,beam_size_fwhp,CMB_T)\n",
    "    binned_ell_cur, binned_spectrum_cur = calculate_2d_spectrum(CMB_T_convolved*window,CMB_T_convolved*window,delta_ell,ell_max,pix_size,N)\n",
    "    signal_only[i,:] = binned_spectrum_cur\n",
    "    sys.stdout.write(\"\\r signal only sims, iterations complete: %d of %d\" % ((i+1),N_iterations) )\n",
    "    sys.stdout.flush()\n",
    "    i = i + 1\n",
    "\n",
    "def average_N_spectra(spectra,N_spectra,N_ells):\n",
    "    avgSpectra = np.zeros(N_ells)\n",
    "    rmsSpectra = np.zeros(N_ells)\n",
    "    \n",
    "    # calcuate the average spectrum\n",
    "    i = 0\n",
    "    while (i < N_spectra):\n",
    "        avgSpectra = avgSpectra + spectra[i,:]\n",
    "        i = i + 1\n",
    "    avgSpectra = avgSpectra/(1. * N_spectra)\n",
    "    \n",
    "    #calculate the rms of the spectrum\n",
    "    i =0\n",
    "    while (i < N_spectra):\n",
    "        rmsSpectra = rmsSpectra +  (spectra[i,:] - avgSpectra)**2\n",
    "        i = i + 1\n",
    "    rmsSpectra = np.sqrt(rmsSpectra/(1. * N_spectra))\n",
    "    \n",
    "    return(avgSpectra,rmsSpectra)\n",
    "\n",
    "\n",
    "sig_only_mean_spectrum, rms_not_needed = average_N_spectra(signal_only,N_iterations,int(ell_max/delta_ell))\n",
    "    \n",
    "sub_sampled_CLs = DlTT[binned_ell] * 2. * np.pi / (binned_ell * (binned_ell+1.))\n",
    "\n",
    "Multiplicative_Bias_est =  sub_sampled_CLs / sig_only_mean_spectrum\n",
    "   \n",
    "## make some plots\n",
    "plt.semilogy(binned_ell,binned_spectrum* binned_ell * (binned_ell+1.)/2. / np.pi,color='b')\n",
    "plt.semilogy(binned_ell,(sig_only_mean_spectrum)* binned_ell * (binned_ell+1.)/2. / np.pi,color='g')\n",
    "plt.semilogy(binned_ell,(Multiplicative_Bias_est),color='b')\n",
    "plt.semilogy(binned_ell,(binned_spectrum)*Multiplicative_Bias_est* binned_ell * (binned_ell+1.)/2. / np.pi,color='y')\n",
    "plt.semilogy(ell,DlTT,color='r')\n",
    "plt.ylabel('$D_{\\ell}$ [$\\mu$K$^2$]')\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot shows the estimate of the CMB power spectrum after correcting for the multiplicative bias (transfer function) in yellow.   In addition we show (red) the input CMB power spectrum, (green) the average of the signal only simulations, (blue, lower) the transfer function, and (blue, upper) the naive power spectrum of our map.  Consider how all these curves relate to creating the yellow estimate.\n",
    "<font color='red'>EXCERCISE: </font>  Why does the green curve look so much like the red curved but surpressed?  What is the meaning of the lower (and smoother) blue curve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your comments go here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calibrating the noise bias\n",
    "\n",
    "The noise bias can be computed by running noise only simulations through the naive power spectrum estimator and computing the average power spectrum.\n",
    "\n",
    "NOTE: An alternative approach exists for dealing with the noise.  If you can subdivide your data into subsets with common signal but independent noise, one can compute \"cross-spectra\" between these subsets.  (You compute these by doing a 2d FFT on each subset and then multiplying one by the complex conjugate of the other.)   This results in some information loss (since you are throwing out the auto-correlation of each subset, but it completely eliminates potential measurement bias from an incorrect noise model.  \n",
    "\n",
    "As a two part exercise: (1) use the wrong noise model in analyzing the power spectrum to see what happens, and (2) implement a cross-spectrum estimator to see that this noise bias goes away with cross spectra.  Also note that the error bars grow with the cross spectrum. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "noise_only  = np.zeros([N_iterations,int(ell_max/delta_ell)])\n",
    "i = 0\n",
    "while (i <N_iterations):\n",
    "    Noise = make_noise_map(N,pix_size,white_noise_level,atmospheric_noise_level,one_over_f_noise_level)\n",
    "    binned_ell_cur, binned_spectrum_cur = calculate_2d_spectrum(Noise*window,Noise*window,delta_ell,ell_max,pix_size,N)\n",
    "    noise_only[i,:] = binned_spectrum_cur\n",
    "    sys.stdout.write(\"\\r noise only sims, iterations complete: %d of %d\" % ((i+1),N_iterations) )\n",
    "    sys.stdout.flush()\n",
    "    i = i + 1\n",
    "\n",
    "noise_only_mean_spectrum, rms_not_needed = average_N_spectra(noise_only,N_iterations,int(ell_max/delta_ell))\n",
    "\n",
    "Additive_Bias_est = noise_only_mean_spectrum\n",
    "\n",
    "\n",
    "#plt.semilogy(binned_ell,binned_spectrum* binned_ell * (binned_ell+1.)/2. / np.pi,color='b')\n",
    "#plt.semilogy(binned_ell,(sig_only_mean_spectrum)* binned_ell * (binned_ell+1.)/2. / np.pi,color='g')\n",
    "plt.semilogy(binned_ell,(noise_only_mean_spectrum)* binned_ell * (binned_ell+1.)/2. / np.pi,color='y')\n",
    "#plt.semilogy(binned_ell,(Multiplicative_Bias_est),color='b')\n",
    "plt.semilogy(binned_ell,(binned_spectrum)*Multiplicative_Bias_est*\\\n",
    "             binned_ell * (binned_ell+1.)/2. / np.pi,color='b')\n",
    "plt.semilogy(binned_ell,(binned_spectrum -noise_only_mean_spectrum)\\\n",
    "             *Multiplicative_Bias_est* binned_ell * (binned_ell+1.)/2. / np.pi,color='g')\n",
    "plt.semilogy(ell,DlTT,color='r')\n",
    "plt.ylabel('$D_{\\ell}$ [$\\mu$K$^2$]')\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The green curve in this plot shows our unbiased estimate for the spectrum.  This estimate includes corrections for the noise bias and the transfer function.  The yellow curve shows our estimate for the noise only additive bias.  The blue curve shows the spectrum accounting for only the multiplicative bias as was done above.  The red curve shows the underlying power spectrum used in generating our simulated map.\n",
    "\n",
    "At this point we have an unbiased estimate of the power spectrum (the green curve). Note that at high-$\\ell$  the spectrum is dominated by the SZ and point source components.  Next we need to put error bars on this measurement.  Again we do this with the monty carlo approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXCERCISE: </font>  What is the yellow curve?  How is it generated, how is it used to get the green curve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "your comments go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantifying the error bars\n",
    "\n",
    "The error bars are computed by generating simulations including signal and noise, computing the naive power spectrum, taking the RMS of these results and then subtracting the noise bias and accounting for the transfer function.  In effect we are simulating running the experiment a bunch of times and looking at the RMS of the results--- this approach is broadly applicable to other problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SplusN  = np.zeros([N_iterations,int(ell_max/delta_ell)])\n",
    "i = 0\n",
    "while (i <N_iterations):\n",
    "    CMB_T = make_CMB_T_map(N,pix_size,ell,DlTT)\n",
    "\n",
    "    PSMap = Poisson_source_component(N,pix_size,Number_of_Sources,Amplitude_of_Sources) \n",
    "    PSMap +=  Exponential_source_component(N,pix_size,Number_of_Sources_EX,Amplitude_of_Sources_EX)\n",
    "    SZMap,trash = SZ_source_component(N,pix_size,Number_of_SZ_Clusters,\\\n",
    "                                      Mean_Amplitude_of_SZ_Clusters,SZ_beta,SZ_Theta_core,False)\n",
    "    \n",
    "    CMB_T  = CMB_T + PSMap + SZMap\n",
    "\n",
    "    CMB_T_convolved = convolve_map_with_gaussian_beam(N,pix_size,beam_size_fwhp,CMB_T)\n",
    "    Noise = make_noise_map(N,pix_size,white_noise_level,atmospheric_noise_level,one_over_f_noise_level)\n",
    "    binned_ell_cur, binned_spectrum_cur = calculate_2d_spectrum((CMB_T_convolved+Noise)*window\\\n",
    "                                                                ,(CMB_T_convolved+Noise)*window\\\n",
    "                                                                ,delta_ell,ell_max,pix_size,N)\n",
    "    SplusN[i,:] = binned_spectrum_cur\n",
    "    sys.stdout.write(\"\\r signal and noise sims, iterations complete: %d of %d\" % ((i+1),N_iterations) )\n",
    "    sys.stdout.flush()\n",
    "    i = i + 1\n",
    "    \n",
    "mean_not_needed,rms_sig_plus_noise = average_N_spectra(SplusN,N_iterations,int(ell_max/delta_ell))\n",
    "\n",
    "\n",
    "\n",
    "plt.errorbar(binned_ell, (binned_spectrum -noise_only_mean_spectrum)*Multiplicative_Bias_est* \\\n",
    "             binned_ell * (binned_ell+1.)/2. / np.pi, \\\n",
    "             yerr=rms_sig_plus_noise*Multiplicative_Bias_est* binned_ell * (binned_ell+1.)/2. / np.pi)\n",
    "plt.semilogy(ell,DlTT,color='r')\n",
    "plt.ylabel('$D_{\\ell}$ [$\\mu$K$^2$]')\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.show()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And there you have it.  That is how you compute a CMB power spectrum and error bars.   If you want to fit cosmology to these data you can re-run CAMB varying cosmological parameters and compute the likelihood difference between these models and the data.  This is left as an exercise.\n",
    "<font color='red'>EXCERCISE: </font>  Why dosn't the red curve (input CMB spectrum) fall wihin the error bars of our simulated data set above $\\ell \\sim 2500$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "your comments and discussion go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXCERCISE: </font>  One can worry that the transfer funciton calcuation depends on the input spectrum used in its calculation.  Modify the spectrum used in calculating the transfer funciton to better match the spectrum in the simualtion by adding a constant in quadrature to the input spectrum.   Comment on weather this changes the result or not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "your comments go here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}