Relic_Neutrino_Forecast_Matching.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Goal\n",
    "\n",
    "Show that for a given neutrino forecast, with some fiducial cosmology, that we can reproduce the same results using a forecast for an equivlent dark matter relic.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistics\n",
    "\n",
    "1) For a degenerate three-neutrino scenario, with the mass of a single neutrino given by $m_{\\nu}$, we must forecast dark matter relic with mass $M_{relic} = m_\\nu * 3^{1/4}$\n",
    "\n",
    "2) For a degenerate three-neutrino scenario, with the temperature of a single neutrino given by $T_{\\nu}$, we must forecast dark matter relic with temperature $T_{relic} = T_\\nu * 3^{1/4}$\n",
    "\n",
    "3) When we run the Fisher forecast for the relic, we must fix $T_{relic}$ and vary $M_{relic}$\n",
    "\n",
    "4) When comparing the results, we should compare forecast for $\\omega_{ncdm}$ and not for $M_{ncdm}$. This is because our decisions in (1) and (2) above are to produce equivalent affects on the abundance of the neutrinos/relics. It will not, for example, produce the same uncertainties for $M_\\nu$ and $M_{relic}$.\n",
    "\n",
    "5) In both cases, we wish to ignore contributions to $N_{eff}$ by relativistic neutrinos. So, we reduce $N_{eff}$ by $3 * 1.0132$ in both cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup Analysis Environment\n",
    "\n",
    "1) Choose a directory to house your project in: \n",
    "```\n",
    ".../<project-directory>\n",
    "```\n",
    "\n",
    "\n",
    "2) Create and activate a fresh Python3 virtual environment (we use Python 3.6.8) there if you'd like: \n",
    "```\n",
    "$ cd .../<project-directory>\n",
    "$ python -m virtualenv env \n",
    "$ source env/bin/activate\n",
    "```\n",
    "\n",
    "3) Download the `cosmicfish` package from Git: \n",
    "```\n",
    "$ git clone git@github.com:ndeporzio/cosmicfish.git\n",
    "```\n",
    "\n",
    "4) Install the `cosmicfish` package. Note that its dependencies will install automatically.\n",
    "```\n",
    "$ cd cosmicfish\n",
    "$ pip install . \n",
    "```\n",
    "\n",
    "5) Launch Jupyter and open `tutorial.ipynb` notebook using Jupyter browser\n",
    "```\n",
    "$ jupyter notebook\n",
    "```\n",
    "\n",
    "6) Create a data folder where the analysis can store spectrum data. This can be anywhere you'd like - you'll specify the path below. \n",
    "```\n",
    "$ mkdir <project-directory>/data\n",
    "```\n",
    "\n",
    "7) Install and build CLASS (if you don't already have a build). Note the `cosmicfish` package includes a method for downloading and installing CLASS for you:\n",
    "```\n",
    "$ python \n",
    ">>> import cosmicfish as cf\n",
    ">>> cf.install_class('<project-directory>/class')\n",
    ">>> exit()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import & Configure _cosmicfish_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the analysis package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cosmicfish as cf "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import relevant python packages... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns \n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other setup steps... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Instruct pyplot to use seaborn \n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify the paths from the setup of your analysis environment above.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set project, data, CLASS directories \n",
    "projectdir = cf.correct_path(\"~/Desktop/cfworkspace/\")\n",
    "datastore = cf.correct_path(\"/Volumes/SSD01/datastore/\")\n",
    "classpath = os.path.join(projectdir, \"class\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify the granularity of numerical derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "derivative_step = 0.01 #How much to vary parameter to calculate numerical derivative\n",
    "mu_integral_step = 0.05 #For calculating numerical integral wrt mu between -1 and 1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Specify Fiducial Cosmologies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neutrinofid = {\n",
    "        \"A_s\" : 2.22e-9, \n",
    "        \"n_s\" : 0.965,\n",
    "        \"omega_b\" : 0.02222,\n",
    "        \"omega_cdm\" : 0.1120,\n",
    "        \"tau_reio\" : 0.06,\n",
    "        \"h\" : 0.70,\n",
    "        \"T_cmb\" : 2.726, # Units [K]\n",
    "        \"N_ncdm\" : 3., \n",
    "        \"T_ncdm\" : 1.95/2.726, \n",
    "        \"m_ncdm\" : 0.03, # Units [eV]\n",
    "        \"b0\" : 1.0, \n",
    "        \"alphak2\" : 1.0,\n",
    "        \"sigma_fog_0\" : 250000, #Units [m s^-2]\n",
    "        \"N_eff\" : 3.046 - (3 * 1.0132), \n",
    "        \"relic_fix\" : None #Not used for neutrino forecasts\n",
    "        } \n",
    "\n",
    "# For matching the neutrino/relic forecasts, we must modify temp/mass according to equation 44 \n",
    "# of the Cora/Julian paper. That is: \n",
    "# m_chi = m_nu * 3^(1/4) ***where m_nu is the mass of a single neutrino in the degenerate hierarchy\n",
    "# T_chi = 1.95 * 3^(1/4)\n",
    "\n",
    "relicfid = {\n",
    "        \"A_s\" : 2.22e-9, \n",
    "        \"n_s\" : 0.965,\n",
    "        \"omega_b\" : 0.02222,\n",
    "        \"omega_cdm\" : 0.1120,\n",
    "        \"tau_reio\" : 0.06,\n",
    "        \"h\" : 0.70,\n",
    "        \"T_cmb\" : 2.726, # Units [K]\n",
    "        \"N_ncdm\" : 1., \n",
    "        \"T_ncdm\" : (1.95/2.726) * np.power(3., 1./4.), # Units [T_cmb]  Are Julian's matrices at 1.64 K? \n",
    "        \"m_ncdm\" : 0.03 * np.power(3., 1./4.), # Units [eV]\n",
    "        \"b0\" : 1.0, \n",
    "        \"alphak2\" : 1.0,\n",
    "        \"sigma_fog_0\" : 250000, #Units [m s^-2]\n",
    "        \"N_eff\" : 3.046 - (3 * 1.0132), \n",
    "        \"relic_fix\" : \"T_ncdm\" # Fix T_ncdm or M_ncdm \n",
    "        }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Specify Experiment Observational Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Specify redshift bins, noise per bin, and sky coverage\n",
    "z_table = np.arange(0.65, 1.85, 0.1)\n",
    "dNdz = np.array([309., 2269., 1923., 2094., 1441., 1353., 1337., 523., 466., 329., 126., 0., 0.])\n",
    "skycover = 14000. # Sky coverage of survey in degrees^2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstrate Convergence\n",
    "\n",
    "Before running the forecast, we want to ensure our cosmological parameters are well converged about the points we are interested in using to calculate Fisher matrices. To do so, we can use the `convergence` class of `cosmicfish`. \n",
    "\n",
    "We pass to `convergence` some fiducial cosmology, and then it will vary the parameters of that fiducial cosmology by step sizes specified by the user and compute the corresponding power spectrum derivatives. \n",
    "\n",
    "All 'varyfactors' are computed relative to the fiducial cosmology. That is:\n",
    "\n",
    "dtheta = varyfactor * theta_fiducial\n",
    "\n",
    "**WARNING: This process takes a considerable amount of time (~5 mins per varyfactor w/o prexisting data). Only run when you need to.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neutrinoconvergencetest = cf.convergence(\n",
    "    classpath, # Path to CLASS installation\n",
    "    datastore, # Path to directory holding CLASS output data\n",
    "    'neutrino', # 'relic' or 'neutrino' forecasting scheme \n",
    "    neutrinofid, # The fiducial cosmology \n",
    "    z_table, # Redshift steps in observation\n",
    "    dNdz, # Redshift noise in observation\n",
    "    fcoverage_deg=14000, # Sky coverage in observation\n",
    "    RSD=True, # Use RSD correction to Pm\n",
    "    FOG=True, # Use FOG correction to Pm\n",
    "    AP=True, # Use AP correction to PM\n",
    "    COV=True, #Use AP Change of Variables correction to PM\n",
    "    mu_step=mu_integral_step,\n",
    "    varyfactors=[0.003, 0.005, 0.007] # Relative factors used to compute convergence\n",
    "    )\n",
    "neutrinoconvergencetest.gen_all_plots() # Display convergence plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relicconvergencetest = cf.convergence(\n",
    "    classpath, # Path to CLASS installation\n",
    "    datastore, # Path to directory holding CLASS output data\n",
    "    'relic', # 'relic' or 'neutrino' forecasting scheme \n",
    "    relicfid, # The fiducial cosmology \n",
    "    z_table, # Redshift steps in observation\n",
    "    dNdz, # Redshift noise in observation\n",
    "    fcoverage_deg=14000, # Sky coverage in observation\n",
    "    RSD=True, # Use RSD correction to Pm\n",
    "    FOG=True, # Use FOG correction to Pm\n",
    "    AP=True, # Use AP correction to PM\n",
    "    COV=True, #Use AP Change of Variables correction to PM\n",
    "    mu_step=mu_integral_step,\n",
    "    varyfactors=[0.003, 0.005, 0.007] # Relative factors used to compute convergence\n",
    "    )\n",
    "relicconvergencetest.gen_all_plots() # Display convergence plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run Forecasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neutrinoforecast = cf.forecast(\n",
    "    classpath, \n",
    "    datastore, \n",
    "    'neutrino', \n",
    "    neutrinofid, \n",
    "    z_table, \n",
    "    dNdz, \n",
    "    fcoverage_deg=skycover, \n",
    "    dstep=derivative_step,\n",
    "    RSD=True,\n",
    "    FOG=True,\n",
    "    AP=True,\n",
    "    COV=True)\n",
    "\n",
    "neutrinoforecast.gen_pm()\n",
    "\n",
    "neutrinoforecast.gen_fisher(\n",
    "    fisher_order=[\n",
    "        'omega_b',                                    \n",
    "        'omega_cdm',                                  \n",
    "        'n_s',                                        \n",
    "        'A_s',                                        \n",
    "        'tau_reio',                                   \n",
    "        'h',                                          \n",
    "        #'M_ncdm',                                    \n",
    "        'omega_ncdm',                                 \n",
    "        'sigma_fog',                                   \n",
    "        'b0',                                         \n",
    "        'alpha_k2'],\n",
    "    mu_step=mu_integral_step, \n",
    "    skipgen=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relicforecast = cf.forecast(\n",
    "    classpath, \n",
    "    datastore, \n",
    "    'relic', \n",
    "    relicfid, \n",
    "    z_table, \n",
    "    dNdz, \n",
    "    fcoverage_deg=skycover, \n",
    "    dstep=derivative_step,\n",
    "    RSD=True,\n",
    "    FOG=True,\n",
    "    AP=True,\n",
    "    COV=True)\n",
    "\n",
    "relicforecast.gen_pm()\n",
    "\n",
    "relicforecast.gen_fisher(\n",
    "    fisher_order=[\n",
    "        'omega_b',                                    \n",
    "        'omega_cdm',                                  \n",
    "        'n_s',                                        \n",
    "        'A_s',                                        \n",
    "        'tau_reio',                                   \n",
    "        'h',                                          \n",
    "        #'M_ncdm',                                    \n",
    "        'omega_ncdm',\n",
    "        #'T_ncdm',\n",
    "        'sigma_fog',                                   \n",
    "        'b0',                                         \n",
    "        'alpha_k2'],\n",
    "    mu_step=mu_integral_step, \n",
    "    skipgen=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Add CMB Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neutrinoforecast.load_cmb_fisher(\n",
    "    fisher_order=[\n",
    "        'omega_b',                                    \n",
    "        'omega_cdm',                                  \n",
    "        'n_s',                                        \n",
    "        'A_s',                                        \n",
    "        'tau_reio',                                   \n",
    "        'h',                                                                             \n",
    "        'omega_ncdm'],\n",
    "    fisherpath=os.path.join(cf.priors_directory(), \"CMBS4_Fisher_Neutrinos.dat\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relicforecast.load_cmb_fisher(\n",
    "    fisher_order=[\n",
    "        'omega_b',                                    \n",
    "        'omega_cdm',                                  \n",
    "        'n_s',                                        \n",
    "        'A_s',                                        \n",
    "        'tau_reio',                                   \n",
    "        'h',                                                                             \n",
    "        'omega_ncdm'],\n",
    "    fisherpath=os.path.join(cf.priors_directory(), \"CMBS4_Fisher_Neutrinos.dat\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save Fisher Matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neutrinoforecast.export_matrices(\"~/Desktop\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relicforecast.export_matrices(\"~/Desktop\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare Results - Check for Equivalance "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(np.array((neutrinoforecast.numpy_lss_fisher - relicforecast.numpy_lss_fisher)/neutrinoforecast.numpy_lss_fisher), \n",
    "            columns=[\n",
    "        'omega_b',                                    \n",
    "        'omega_cdm',                                  \n",
    "        'n_s',                                        \n",
    "        'A_s',                                        \n",
    "        'tau_reio',                                   \n",
    "        'h',                                          \n",
    "        #'M_ncdm',                                    \n",
    "        'omega_ncdm',                                 \n",
    "        'sigma_fog',                                   \n",
    "        'b0',                                         \n",
    "        'alpha_k2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3.7.4",
   "language": "python",
   "name": "3.7.4"
  },
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}