examples/car_simulation_to_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From CMB simulations to power spectra analysis\n",
    "\n",
    "In this tutorial, we will show how to generate CMB simulations from theoritical $C_\\ell$, produce the corresponding CMB map tiles and finally analyze them using `psplay`. \n",
    "\n",
    "## CMB simulation\n",
    "\n",
    "Using the attached $C_\\ell$ [file](bode_almost_wmap5_lmax_1e4_lensedCls_startAt2.dat), we will first generate simulations with `pspy` in `CAR` pixellisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pspy import so_map\n",
    "\n",
    "ncomp= 3\n",
    "ra0, ra1, dec0, dec1 = -30, 30, -30, 30\n",
    "res = 0.5\n",
    "\n",
    "template = so_map.car_template(ncomp, ra0, ra1, dec0, dec1, res)\n",
    "\n",
    "cl_file = \"bode_almost_wmap5_lmax_1e4_lensedCls_startAt2.dat\"\n",
    "cmb = template.synfast(cl_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we make 2 splits out of it, each with $5$ µK.arcmin rms in temperature and $5\\times\\sqrt{2}$ µK.arcmin in polarisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "nsplits = 2\n",
    "splits = [cmb.copy() for i in range(nsplits)]\n",
    "for i in range(nsplits):\n",
    "    noise = so_map.white_noise(cmb, rms_uKarcmin_T=5, rms_uKarcmin_pol=np.sqrt(2)*5)\n",
    "    splits[i].data += noise.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We finally write on disk the two corresponding `fits` files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(nsplits):\n",
    "    splits[i].write_map(\"split{}_IQU_CAR.fits\".format(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generation of tile files\n",
    "\n",
    "We now have a total of 6 maps in `CAR` pixellisation we want to plot them into tile files. A tile file corresponds to a static `PNG` file representing a part of the sky at a given zoom level. To convert the maps, we will use a set of tools from `psplay`. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conversion to tile files\n",
    "You should get two files of ~ 1 Gb size. We should keep these files since we will use them later for the power spectra computation. Last step in the conversion process, we have to generate the different tiles "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from psplay import tools\n",
    "\n",
    "for i in range(nsplits):\n",
    "    tools.car2tiles(input_file=\"split{}_IQU_car.fits\".format(i), \n",
    "                    output_dir=\"tiles/split{}_IQU_car.fits\".format(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the end of the process, we get a new directory `tiles` with two sub-directories `split0_IQU_car.fits` and `split1_IQU_car_fits`. Within these two directories we have 6 directories which names refer to the zoom level : the directory `-5` correspond to the smallest zoom whereas the directory `0` corresponds to the most precise tiles. \n",
    "\n",
    "## Using `psplay` to visualize the maps\n",
    "\n",
    "Now that we have generated the different tiles, we can interactively see the maps with `psplay`. The configuration of `psplay` can be either done using a dictionary or more easily with a `yaml`file. Basically, the configuration file needs to know where the original `FITS` files and the tile files are located. Other options can be set but we won't enter into too much details for this tutorial. \n",
    "\n",
    "For the purpose of this tutorial we have already created the `yaml` [file](simulation_to_analysis.yml) associated to the files we have created so far. Here is a copy-paste of it\n",
    "\n",
    "```yaml\n",
    "map:\n",
    "  layers:\n",
    "    - cmb:\n",
    "        tags:\n",
    "          splits:\n",
    "            values: [0, 1]\n",
    "            keybindings: [j, k]\n",
    "\n",
    "          components:\n",
    "            values: [0, 1, 2]\n",
    "            keybindings: [c, v]\n",
    "            substitutes: [\"T\", \"Q\", \"U\"]\n",
    "\n",
    "        tile: files/tiles/split{splits}_IQU_car.fits/{z}/tile_{y}_{x}_{components}.png\n",
    "        name: CMB simulation - split {splits} - {components}\n",
    "\n",
    "data:\n",
    "  maps:\n",
    "    - id: split0\n",
    "      file: split0_IQU_car.fits\n",
    "    - id: split1\n",
    "      file: split1_IQU_car.fits\n",
    "\n",
    "  theory_file: bode_almost_wmap5_lmax_1e4_lensedCls_startAt2.dat\n",
    "  \n",
    "plot:\n",
    "    lmax: 2000\n",
    "```\n",
    "\n",
    "There are 3 sections related to the 3 main steps namely the map visualization, the spectra computation and the graphical representation of the spectra. The two first sections are mandatory. The `map` section corresponds to the tile files generated so far and can be dynamically expanded given different `tags`. Here for instance, we will built all the combination of split and component values. The tile and the name fields will be generated for each combination given the tag values. Dedicated keybindings can also be defined in order to switch between the different split and/or components.\n",
    "\n",
    "**The tricker part of this configuration is to set the path to tiles relatively to where your notebook/JupyterLab instance has been started**. We can't set an absolute path and so you have to make sure that your notebook has been initiated from the `examples` directory. Otherwise, you should change the path to tile files given that you have access to them from your JupyterLab instance. So make sure to initiate your JupyterLab session from a \"top\" directory.\n",
    "\n",
    "We can now create an instance of `psplay` application and show the different maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from psplay import App\n",
    "my_app = App(\"simulation_to_analysis.yml\")\n",
    "my_app.show_map()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we unzoom enough, we will see the $\\pm$ 30° patch size. We can also switch between the different I, Q, U and split layers. There are other options like the colormap and color scale which one can play and see how things change.\n",
    "\n",
    "## Selecting sub-patches and computing the corresponding power spectra\n",
    "\n",
    "Given the different map, we can now select patches by clicking on the square or the disk icons located just below the +/- zoom button. For instance, if we select a rectangle and another disk whose size are more or less the total size of our patch, we will get two surfaces of almost 3000 square degrees. Now we can ask `psplay` to compute the power spectra of both regions. Let's initiate the plot application"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_app.show_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now clik on the `Compute spectra` button and watch the log output by clicking on the `Logs` tab. Depending on your machine (mainly memory capacity), the process can be pretty long especially the transformation into spherical harmonics. It also depends of the size of your patches but within 1 or 2 minutes, we should get the final spectra for the different cross-correlation. For instance, we can switch between different combinations of spectra (TT, TE, EE,...). \n",
    "\n",
    "Finally, the `Configuration` tab offer different options like changing the $\\ell_\\mathrm{max}$ value of the method use for the computation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}