Testing2RelicForecast.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Test 1: gforecast Branch, g Forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T01:14:58.007185Z",
     "start_time": "2020-05-04T22:34:50.039003Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns \n",
    "import cosmicfish as cf \n",
    "import matplotlib.pyplot as plt\n",
    "import dill\n",
    "\n",
    "# Instruct pyplot to use seaborn \n",
    "sns.set()\n",
    "\n",
    "# Set project, data, CLASS directories \n",
    "datastore = cf.correct_path(\"/Volumes/SSD01/data.nosynctest1/\")\n",
    "classpath = cf.correct_path(\"/Users/nicholasdeporzio/Desktop/cfworkspace/class\")\n",
    "\n",
    "# Specify resolution of numerical integrals\n",
    "derivative_step = 0.008 # 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 \n",
    "\n",
    "# Generate output paths  \n",
    "ps18_resultsdir = cf.correct_path(\"/Volumes/SSD01/results/test1/\")\n",
    "cf.makedirectory(ps18_resultsdir)\n",
    "\n",
    "# Linda Fiducial Cosmology \n",
    "ps18_fid = {\n",
    "        \"A_s\" : 2.2321e-9, \n",
    "        \"n_s\" : 0.967,\n",
    "        \"omega_b\" : 0.02226,\n",
    "        \"omega_cdm\" : 0.1127,\n",
    "        \"tau_reio\" : 0.0598,\n",
    "        \"h\" : 0.701,\n",
    "        \"T_cmb\" : 2.726, # Units [K]\n",
    "        \"N_ncdm\" : 1., \n",
    "        \"deg_ncdm\" : 1,\n",
    "        \"T_ncdm\" : (0.91/2.726), # Units [T_cmb]. \n",
    "        \"m_ncdm\" : 0.01, # Units [eV]\n",
    "        \"b0\" : 1.0, \n",
    "        \"alphak2\" : 1.0,\n",
    "        \"sigma_fog_0\" : 250000, #Units [m s^-2]\n",
    "        \"N_eff\" : 3.046, #We allow relativistic neutrinos in addition to our DM relic\n",
    "        \"relic_vary\" : \"N_ncdm\" # Fix T_ncdm or m_ncdm \n",
    "        }\n",
    "\n",
    "# DESI Parameters\n",
    "z_table = np.array([0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85])\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\n",
    "\n",
    "# Run Fisher Forecast\n",
    "masses = np.geomspace(0.01, 10., 21)\n",
    "temps = np.array([0.91])\n",
    "\n",
    "omegacdm_set = np.array([ps18_fid['omega_cdm'] - ((masses/cf.NEUTRINO_SCALE_FACTOR)*np.power(tval / 1.95, 3.)) for tidx, tval in enumerate(temps)])                                  \n",
    "ps18_fiducialset = [[dict(ps18_fid, **{\n",
    "    'm_ncdm' : masses[midx], \n",
    "    'omega_cdm' : omegacdm_set[tidx, midx], \n",
    "    'T_ncdm' : temps[tidx]/2.726}) \n",
    "    for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)] \n",
    "ps18_forecastset = [[cf.forecast(\n",
    "    classpath, \n",
    "    datastore, \n",
    "    'relic', \n",
    "    fidval, \n",
    "    z_table, \n",
    "    \"DESI\",\n",
    "    dNdz, \n",
    "    fcoverage_deg=skycover, \n",
    "    dstep=derivative_step,\n",
    "    RSD=True,\n",
    "    FOG=True,\n",
    "    AP=True,\n",
    "    COV=True) \n",
    "    for fididx, fidval in enumerate(fidrowvals)] \n",
    "    for fidrowidx, fidrowvals in enumerate(ps18_fiducialset)]\n",
    "\n",
    "for frowidx, frowval in enumerate(ps18_forecastset):\n",
    "    for fidx, fcst in enumerate(frowval):\n",
    "        if type(fcst.fisher)==type(None):\n",
    "            fcst.gen_pm()\n",
    "            fcst.gen_fisher(\n",
    "                fisher_order=[\n",
    "                    'omega_b',\n",
    "                    'omega_cdm',\n",
    "                    'n_s',\n",
    "                    'A_s',\n",
    "                    'tau_reio',\n",
    "                    'h',\n",
    "                    'N_ncdm',\n",
    "                    'sigma_fog',\n",
    "                    'b0',\n",
    "                    'alpha_k2'],\n",
    "                mu_step=mu_integral_step,\n",
    "                skipgen=False)\n",
    "            print(\"Relic Forecast \", fidx, \" complete...\")\n",
    "            dill.dump_session(os.path.join(ps18_resultsdir, 'test1.db'))\n",
    "        else:\n",
    "            print('Fisher matrix already generated!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T11:45:22.236207Z",
     "start_time": "2020-05-05T11:44:30.414274Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import dill\n",
    "dill.load_session('/Volumes/SSD01/results/test1/test1.db')\n",
    "\n",
    "ps18_resultsdir = cf.correct_path(\"/Volumes/SSD01/results/test1/\")\n",
    "\n",
    "#cmbsurvey=\"CMBS4\"\n",
    "cmbsurvey=\"Planck\"\n",
    "\n",
    "prepath = cf.priors_directory()\n",
    "tempnames=['0.91']\n",
    "\n",
    "head = 'A_s \\t h \\t n_s \\t omega_b \\t omega_cdm \\t tau_reio \\t N_ncdm \\t sigma_fog \\t beta0 \\t alpha_k2'\n",
    "\n",
    "\n",
    "for tidx, tval in enumerate(ps18_forecastset):\n",
    "    for midx, mval in enumerate(tval): \n",
    "        mval.reorder_fisher(['A_s', 'h', 'n_s', 'omega_b', 'omega_cdm', 'tau_reio', 'N_ncdm', 'sigma_fog', 'b0', 'alpha_k2'])\n",
    "        mval.load_cmb_fisher(\n",
    "              fisher_order=[\n",
    "                'A_s',\n",
    "                'h',  \n",
    "                'n_s',\n",
    "                'omega_b',                                    \n",
    "                'omega_cdm',                                                     \n",
    "                'tau_reio',                                                                                    \n",
    "                'N_ncdm'],\n",
    "            fisherpath = cf.priors_directory()+'/Relic_Fine_Priors/Fisher_T'+tempnames[tidx]+'K/Fisher'+cmbsurvey+'_bin'+str(midx+1)+'_T'+tempnames[tidx]+'K.txt'\n",
    "        )\n",
    "        outdir=os.path.join(ps18_resultsdir, cmbsurvey, 'T'+tempnames[tidx], 'M'+str(midx+1))\n",
    "        cf.makedirectory(outdir)\n",
    "        mval.export_matrices(outdir)\n",
    "        np.savetxt(outdir + '_Full_Fisher.dat', mval.numpy_full_fisher, delimiter='\\t', header=head)\n",
    "        mval.pandas_lss_fisher.to_csv(outdir + \"_LSS_Fisher.dat\", sep='\\t', index=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T11:42:04.298551Z",
     "start_time": "2020-05-05T11:42:04.290653Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T11:59:49.783386Z",
     "start_time": "2020-05-05T11:59:49.453449Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T11:58:16.271109Z",
     "start_time": "2020-05-05T11:58:16.265012Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Test 2: gforecast Branch, M Forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns \n",
    "import cosmicfish as cf \n",
    "import matplotlib.pyplot as plt\n",
    "import dill\n",
    "\n",
    "# Instruct pyplot to use seaborn \n",
    "sns.set()\n",
    "\n",
    "# Set project, data, CLASS directories \n",
    "datastore = cf.correct_path(\"/Volumes/SSD01/data.nosynctest2/\")\n",
    "classpath = cf.correct_path(\"/Users/nicholasdeporzio/Desktop/cfworkspace/class2\")\n",
    "\n",
    "# Specify resolution of numerical integrals\n",
    "derivative_step = 0.008 # 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 \n",
    "\n",
    "# Generate output paths  \n",
    "ps2_resultsdir = cf.correct_path('/Volumes/SSD01/results/test2/')\n",
    "cf.makedirectory(ps2_resultsdir)\n",
    "\n",
    "#Fiducial Cosmology\n",
    "ps2_fid = {                                                                    \n",
    "        \"A_s\" : 2.2321e-9,                                                      \n",
    "        \"n_s\" : 0.967,                                                          \n",
    "        \"omega_b\" : 0.02226,                                                    \n",
    "        \"omega_cdm\" : 0.1127,                                                   \n",
    "        \"tau_reio\" : 0.0598,                                                    \n",
    "        \"h\" : 0.701,                                                            \n",
    "        \"T_cmb\" : 2.726, # Units [K]                                            \n",
    "        \"N_ncdm\" : 3.0,                                                          \n",
    "        \"deg_ncdm\" : 1.0,                                                         \n",
    "        \"T_ncdm\" : (1.95/2.726), # Units [T_cmb].                             \n",
    "        \"m_ncdm\" : 0.02, # 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), #We allow relativistic neutrinos in addition to our DM relic\n",
    "        \"relic_vary\" : None # Fix T_ncdm or m_ncdm                          \n",
    "        }\n",
    "\n",
    "# DESI survey parameters\n",
    "z_table = np.array([0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85])\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\n",
    "\n",
    "# Run Fisher Forecast for cosmology\n",
    "ps2_forecast = cf.forecast(\n",
    "    classpath, \n",
    "    datastore, \n",
    "    'neutrino', \n",
    "    ps2_fid, \n",
    "    z_table, \n",
    "    \"DESI\",\n",
    "    dNdz, \n",
    "    fcoverage_deg=skycover,\n",
    "    dstep=derivative_step,\n",
    "    RSD=True,\n",
    "    FOG=True,\n",
    "    AP=True,\n",
    "    COV=True)\n",
    "ps2_forecast.gen_pm()\n",
    "ps2_forecast.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",
    "        'sigma_fog',                                   \n",
    "        'b0',                                         \n",
    "        'alpha_k2'\n",
    "    ],\n",
    "    mu_step=mu_integral_step, \n",
    "    skipgen=False)\n",
    "dill.dump_session(os.path.join(ps2_resultsdir, 'test2.db'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-02T06:38:09.139150Z",
     "start_time": "2020-05-02T06:38:08.258574Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import dill\n",
    "dill.load_session('/Volumes/SSD01/results/test2/test2.db')\n",
    "\n",
    "ps2_resultsdir = cf.correct_path(\"/Volumes/SSD01/results/test2/\")\n",
    "\n",
    "#survey=\"CMBS4\"\n",
    "survey=\"Planck\"\n",
    "\n",
    "prepath = cf.priors_directory()\n",
    "head = 'omega_b \\t omega_cdm \\t n_s \\t A_s \\t tau_reio \\t h \\t M_ncdm \\t sigma_fog \\t beta0 \\t alpha_k2'\n",
    "\n",
    "# Add CMB information\n",
    "ps2_forecast.load_cmb_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",
    "        fisherpath=os.path.join(cf.priors_directory(), survey+\"_Fisher_Neutrinos.dat\"))\n",
    "\n",
    "outdir=os.path.join(ps2_resultsdir, survey)\n",
    "cf.makedirectory(outdir)\n",
    "ps2_forecast.export_matrices(outdir)\n",
    "np.savetxt(outdir + '/Full_Fisher.dat', ps2_forecast.numpy_full_fisher, delimiter='\\t', header=head)\n",
    "ps2_forecast.pandas_lss_fisher.to_csv(outdir + \"/LSS_Fisher.dat\", sep='\\t', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Test 3: 2relic Branch, gM Forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T09:55:05.298211Z",
     "start_time": "2020-05-08T00:38:26.808638Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import os                                                                       \n",
    "import shutil                                                                   \n",
    "import numpy as np                                                              \n",
    "import pandas as pd                                                             \n",
    "import seaborn as sns                                                           \n",
    "import cosmicfish as cf                                                         \n",
    "import matplotlib.pyplot as plt                                                 \n",
    "import dill                                                                     \n",
    "                                                                                \n",
    "                                                                                \n",
    "# Instruct pyplot to use seaborn                                                \n",
    "sns.set()                                                                       \n",
    "                                                                                \n",
    "# Set project, data, CLASS directories                                                                                    \n",
    "projectdir = '/Volumes/SSD01/results/test3/'                                    \n",
    "datastore = '/Volumes/SSD01/data.nosynctest3a/'                                  \n",
    "classpath = '/Users/nicholasdeporzio/Desktop/cfworkspace/class/'                \n",
    "                                                                                \n",
    "# Specify resolution of numerical integrals                                     \n",
    "derivative_step = 0.008 # 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 \n",
    "\n",
    "# Generate output paths                                                                \n",
    "ps23_resultsdir = projectdir                                                               \n",
    "cf.makedirectory(ps23_resultsdir)                                                                                        \n",
    "                                                                                \n",
    "# Linda Fiducial Cosmology                                                      \n",
    "ps23_fid = {                                                                    \n",
    "        \"A_s\" : 2.2321e-9,                                                      \n",
    "        \"n_s\" : 0.967,                                                          \n",
    "        \"omega_b\" : 0.02226,                                                    \n",
    "        \"omega_cdm\" : 0.1127,                                                   \n",
    "        \"tau_reio\" : 0.0598,                                                    \n",
    "        \"h\" : 0.701,                                                            \n",
    "        \"T_cmb\" : 2.726, # Units [K]                                            \n",
    "        \"N_ncdm\" : 4.,                                                          \n",
    "        \"deg_ncdm\" : 1.0,                                                       \n",
    "        \"T_ncdm\" : (0.79/2.726), # Units [T_cmb].                               \n",
    "        \"m_ncdm\" : 0.01, # Units [eV]                                           \n",
    "        \"b0\" : 1.0,                                                             \n",
    "        \"alphak2\" : 1.0,                                                        \n",
    "        \"sigma_fog_0\" : 250000, #Units [m s^-2]                                 \n",
    "        \"N_eff\" : 0.0064, #We allow relativistic neutrinos in addition to our DM relic\n",
    "        \"relic_vary\" : \"N_ncdm\",  # Fix T_ncdm or m_ncdm                        \n",
    "        \"m_nu\" : 0.02                                                           \n",
    "        }    \n",
    "\n",
    "# DESI Parameters                                                               \n",
    "z_table = np.array([0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85])\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                         \n",
    "                                                                                \n",
    "# Run Fisher Forecast                                                                                               \n",
    "masses = np.geomspace(0.01, 10., 21)                                                    \n",
    "temps = np.array([0.79, 0.91, 0.94, 1.08])                                                        \n",
    "                                                                                \n",
    "omegacdm_set = np.array([                                                       \n",
    "    ps23_fid['omega_cdm']                                                       \n",
    "    - ((masses/cf.NEUTRINO_SCALE_FACTOR)* np.power(tval / 1.95, 3.))                                            \n",
    "    for tidx, tval in enumerate(temps)])    \n",
    "\n",
    "ps23_fiducialset = [[                                                           \n",
    "    dict(ps23_fid, **{                                                          \n",
    "        'm_ncdm' : masses[midx],                                                \n",
    "        'omega_cdm' : omegacdm_set[tidx, midx],                                 \n",
    "        'T_ncdm' : temps[tidx]/2.726})                                          \n",
    "    for midx, mval in enumerate(masses)]                                        \n",
    "    for tidx, tval in enumerate(temps)]   \n",
    "\n",
    "ps23_forecastset = [[cf.forecast(                                               \n",
    "    classpath,                                                                  \n",
    "    datastore,                                                                  \n",
    "    '2relic',                                                                   \n",
    "    fidval,                                                                     \n",
    "    z_table,                                                                    \n",
    "    \"DESI\",                                                                     \n",
    "    dNdz,                                                                       \n",
    "    fcoverage_deg=skycover,                                                     \n",
    "    dstep=derivative_step,                                                      \n",
    "    RSD=True,                                                                   \n",
    "    FOG=True,                                                                   \n",
    "    AP=True,                                                                    \n",
    "    COV=True)                                                                   \n",
    "    for fididx, fidval in enumerate(fidrowvals)]                                \n",
    "    for fidrowidx, fidrowvals in enumerate(ps23_fiducialset)] \n",
    "\n",
    "#dill.load_session('')                                                          \n",
    "for frowidx, frowval in enumerate(ps23_forecastset):                            \n",
    "    for fidx, fcst in enumerate(frowval):                                       \n",
    "        if type(fcst.fisher)==type(None):                                       \n",
    "            fcst.gen_pm()                                                       \n",
    "            fcst.gen_fisher(                                                    \n",
    "                fisher_order=[                                                  \n",
    "                    'omega_b',                                                  \n",
    "                    'omega_cdm',                                                \n",
    "                    'n_s',                                                      \n",
    "                    'A_s',                                                      \n",
    "                    'tau_reio',                                                 \n",
    "                    'h',                                                        \n",
    "                    'N_ncdm',                                                   \n",
    "                    'M_ncdm',                                                   \n",
    "                    'sigma_fog',                                                \n",
    "                    'b0',                                                       \n",
    "                    'alpha_k2'],                                                \n",
    "                mu_step=mu_integral_step,                                       \n",
    "                skipgen=False)                                                  \n",
    "            print(\"Relic Forecast \", fidx, \" complete...\")                      \n",
    "            dill.dump_session(os.path.join(ps23_resultsdir, 'test3.db'))\n",
    "        else:                                                                   \n",
    "            print('Fisher matrix already generated!')  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T16:45:42.976393Z",
     "start_time": "2020-05-08T16:45:41.955312Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T16:45:46.993356Z",
     "start_time": "2020-05-08T16:45:46.985831Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T17:55:34.591360Z",
     "start_time": "2020-05-08T17:54:24.651310Z"
    },
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import dill\n",
    "dill.load_session('/Volumes/SSD01/results/test3/test3_1.08K.db')\n",
    "import cosmicfish as cf\n",
    "\n",
    "ps23_resultsdir = cf.correct_path(\"/Volumes/SSD01/results/test3/\")\n",
    "\n",
    "#cmbsurvey=\"CMBS4\"\n",
    "cmbsurvey=\"Planck\"\n",
    "\n",
    "prepath = cf.priors_directory()\n",
    "tempnames=['1.08']\n",
    "#head = 'A_s \\t h \\t n_s \\t omega_b \\t omega_cdm \\t tau_reio \\t M_ncdm \\t N_ncdm \\t sigma_fog \\t beta0 \\t alpha_k2'\n",
    "\n",
    "\n",
    "for tidx, tval in enumerate(ps23_forecastset):\n",
    "    for midx, mval in enumerate(tval): \n",
    "        mval.reorder_fisher(['A_s', 'h', 'n_s', 'omega_b', 'omega_cdm', 'tau_reio', 'M_ncdm', 'N_ncdm', 'sigma_fog', 'b0', 'alpha_k2'])\n",
    "        mval.load_cmb_fisher(\n",
    "              fisher_order=[\n",
    "                'A_s',\n",
    "                'h',  \n",
    "                'n_s',\n",
    "                'omega_b',                                    \n",
    "                'omega_cdm',                                                     \n",
    "                'tau_reio',                                                                                    \n",
    "                'M_ncdm', \n",
    "                'N_ncdm'],\n",
    "            fisherpath = cf.priors_directory()+'/2relic_Priors/Fisher_T'+tempnames[tidx]+'K/Fisher'+cmbsurvey+'_bin'+str(midx+1)+'_T'+tempnames[tidx]+'K.txt'\n",
    "        )\n",
    "        outdir=os.path.join(ps23_resultsdir, cmbsurvey, 'T'+tempnames[tidx], 'M'+str(midx+1))\n",
    "        cf.makedirectory(outdir)\n",
    "        mval.export_matrices(outdir)\n",
    "        #np.savetxt(outdir + '_Full_Fisher.dat', mval.numpy_full_fisher, delimiter='\\t', header=head)\n",
    "        #mval.pandas_lss_fisher.to_csv(outdir + \"_LSS_Fisher.dat\", sep='\\t', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-11T09:30:55.114277Z",
     "start_time": "2020-05-11T09:30:55.112054Z"
    },
    "heading_collapsed": true
   },
   "source": [
    "# Test 4: Testing High-Precision CLASS Output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Generate High Precision Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-11T12:51:42.792742Z",
     "start_time": "2020-05-11T10:30:19.015100Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import os                                                                       \n",
    "import shutil                                                                   \n",
    "import numpy as np                                                              \n",
    "import pandas as pd                                                             \n",
    "import seaborn as sns                                                           \n",
    "import cosmicfish as cf                                                         \n",
    "import matplotlib.pyplot as plt                                                 \n",
    "import dill                                                                     \n",
    "                                                                                \n",
    "                                                                                \n",
    "# Instruct pyplot to use seaborn                                                \n",
    "sns.set()                                                                       \n",
    "                                                                                \n",
    "# Set project, data, CLASS directories                                                                                    \n",
    "projectdir = '/Volumes/SSD01/results/test4/'                                    \n",
    "datastore = '/Volumes/SSD01/data.nosynctest4/'                                  \n",
    "classpath = '/Users/nicholasdeporzio/Desktop/cfworkspace/class/'                \n",
    "                                                                                \n",
    "# Specify resolution of numerical integrals                                     \n",
    "derivative_step = 0.008 # 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 \n",
    "mu_integral_step = 1.00\n",
    "\n",
    "# Generate output paths                                                                \n",
    "ps23_resultsdir = projectdir                                                               \n",
    "cf.makedirectory(ps23_resultsdir)                                                                                        \n",
    "                                                                                \n",
    "# Linda Fiducial Cosmology                                                      \n",
    "ps23_fid = {                                                                    \n",
    "        \"A_s\" : 2.2321e-9,                                                      \n",
    "        \"n_s\" : 0.967,                                                          \n",
    "        \"omega_b\" : 0.02226,                                                    \n",
    "        \"omega_cdm\" : 0.1127,                                                   \n",
    "        \"tau_reio\" : 0.0598,                                                    \n",
    "        \"h\" : 0.701,                                                            \n",
    "        \"T_cmb\" : 2.726, # Units [K]                                            \n",
    "        \"N_ncdm\" : 4.,                                                          \n",
    "        \"deg_ncdm\" : 1.0,                                                       \n",
    "        \"T_ncdm\" : (0.79/2.726), # Units [T_cmb].                               \n",
    "        \"m_ncdm\" : 0.01, # Units [eV]                                           \n",
    "        \"b0\" : 1.0,                                                             \n",
    "        \"alphak2\" : 1.0,                                                        \n",
    "        \"sigma_fog_0\" : 250000, #Units [m s^-2]                                 \n",
    "        \"N_eff\" : 0.0064, #We allow relativistic neutrinos in addition to our DM relic\n",
    "        \"relic_vary\" : \"N_ncdm\",  # Fix T_ncdm or m_ncdm                        \n",
    "        \"m_nu\" : 0.02                                                           \n",
    "        }    \n",
    "\n",
    "# DESI Parameters                                                               \n",
    "z_table = np.array([0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85])\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                         \n",
    "                                                                                \n",
    "# Run Fisher Forecast                                                                                               \n",
    "#masses = np.geomspace(0.01, 10., 21)                                                    \n",
    "#temps = np.array([0.79, 0.91, 0.94, 1.08])     \n",
    "masses = np.array([np.geomspace(0.01, 10., 21)[0], np.geomspace(0.01, 10., 21)[13], np.geomspace(0.01, 10., 21)[-1]])\n",
    "temps = np.array([0.79, 1.08])\n",
    "                                                                                \n",
    "omegacdm_set = np.array([                                                       \n",
    "    ps23_fid['omega_cdm']                                                       \n",
    "    - ((masses/cf.NEUTRINO_SCALE_FACTOR)* np.power(tval / 1.95, 3.))                                            \n",
    "    for tidx, tval in enumerate(temps)])    \n",
    "\n",
    "ps23_fiducialset = [[                                                           \n",
    "    dict(ps23_fid, **{                                                          \n",
    "        'm_ncdm' : masses[midx],                                                \n",
    "        'omega_cdm' : omegacdm_set[tidx, midx],                                 \n",
    "        'T_ncdm' : temps[tidx]/2.726})                                          \n",
    "    for midx, mval in enumerate(masses)]                                        \n",
    "    for tidx, tval in enumerate(temps)]   \n",
    "\n",
    "ps23_forecastset = [[cf.forecast(                                               \n",
    "    classpath,                                                                  \n",
    "    datastore,                                                                  \n",
    "    '2relic',                                                                   \n",
    "    fidval,                                                                     \n",
    "    z_table,                                                                    \n",
    "    \"DESI\",                                                                     \n",
    "    dNdz,                                                                       \n",
    "    fcoverage_deg=skycover,                                                     \n",
    "    dstep=derivative_step,                                                      \n",
    "    RSD=False,                                                                   \n",
    "    FOG=False,                                                                   \n",
    "    AP=False,                                                                    \n",
    "    COV=False)                                                                   \n",
    "    for fididx, fidval in enumerate(fidrowvals)]                                \n",
    "    for fidrowidx, fidrowvals in enumerate(ps23_fiducialset)] \n",
    "\n",
    "#dill.load_session('')                                                          \n",
    "for frowidx, frowval in enumerate(ps23_forecastset):                            \n",
    "    for fidx, fcst in enumerate(frowval):                                       \n",
    "        if type(fcst.fisher)==type(None):                                       \n",
    "            fcst.gen_pm()                                                       \n",
    "            fcst.gen_fisher(                                                    \n",
    "                fisher_order=[                                                  \n",
    "                    'omega_b',                                                  \n",
    "                    'omega_cdm',                                                \n",
    "                    'n_s',                                                      \n",
    "                    'A_s',                                                      \n",
    "                    'tau_reio',                                                 \n",
    "                    'h',                                                        \n",
    "                    'N_ncdm',                                                   \n",
    "                    'M_ncdm',                                                   \n",
    "                    'sigma_fog',                                                \n",
    "                    'b0',                                                       \n",
    "                    'alpha_k2'],                                                \n",
    "                mu_step=mu_integral_step,                                       \n",
    "                skipgen=False)                                                  \n",
    "            print(\"Relic Forecast \", fidx, \" complete...\")                      \n",
    "            dill.dump_session(os.path.join(ps23_resultsdir, 'test4_pk_ref.db'))\n",
    "        else:                                                                   \n",
    "            print('Fisher matrix already generated!')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Compare default precision profile with HP cl_permille.pre profile "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import dill \n",
    "dill.load_session('/Volumes/SSD01/results/test4/test4_permille.db')\n",
    "sns.set()\n",
    "sns.set_palette('bright')\n",
    "\n",
    "#mass = '0.010'\n",
    "mass = '0.891'\n",
    "temp = '0.79'\n",
    "param = 'N_ncdm'\n",
    "data = 'cl'\n",
    "midx = 13\n",
    "tidx = 0 \n",
    "red_midx=1\n",
    "red_tidx=0\n",
    "\n",
    "default_path_hi = '/Volumes/SSD01/Cls/T'+temp+'K/M'+mass+'eV/'+param+'/high/test_'+data+'.dat'\n",
    "default_path_lo = '/Volumes/SSD01/Cls/T'+temp+'K/M'+mass+'eV/'+param+'/low/test_'+data+'.dat'\n",
    "hp_path_hi = (eval(('ps23_forecastset['+str(red_tidx)+']['+str(red_midx)+'].'+param+'_high[0].datapath'))[0:-6]+data+'.dat')\n",
    "hp_path_lo = (eval(('ps23_forecastset['+str(red_tidx)+']['+str(red_midx)+'].'+param+'_low[0].datapath'))[0:-6]+data+'.dat')\n",
    "\n",
    "print(default_path_hi)\n",
    "print(hp_path_hi)\n",
    "\n",
    "default_l_hi = np.loadtxt(default_path_hi, skiprows=12, usecols={0})\n",
    "default_TT_hi = np.loadtxt(default_path_hi, skiprows=12, usecols={1})\n",
    "default_EE_hi = np.loadtxt(default_path_hi, skiprows=12, usecols={2})\n",
    "default_TE_hi = np.loadtxt(default_path_hi, skiprows=12, usecols={3})\n",
    "default_l_lo = np.loadtxt(default_path_lo, skiprows=12, usecols={0})\n",
    "default_TT_lo = np.loadtxt(default_path_lo, skiprows=12, usecols={1})\n",
    "default_EE_lo = np.loadtxt(default_path_lo, skiprows=12, usecols={2})\n",
    "default_TE_lo = np.loadtxt(default_path_lo, skiprows=12, usecols={3})\n",
    "\n",
    "hp_l_hi = np.loadtxt(hp_path_hi, skiprows=12, usecols={0})\n",
    "hp_TT_hi = np.loadtxt(hp_path_hi, skiprows=12, usecols={1})\n",
    "hp_EE_hi = np.loadtxt(hp_path_hi, skiprows=12, usecols={2})\n",
    "hp_TE_hi = np.loadtxt(hp_path_hi, skiprows=12, usecols={3})\n",
    "hp_l_lo = np.loadtxt(hp_path_lo, skiprows=12, usecols={0})\n",
    "hp_TT_lo = np.loadtxt(hp_path_lo, skiprows=12, usecols={1})\n",
    "hp_EE_lo = np.loadtxt(hp_path_lo, skiprows=12, usecols={2})\n",
    "hp_TE_lo = np.loadtxt(hp_path_lo, skiprows=12, usecols={3})\n",
    "\n",
    "scaleEE = 100\n",
    "scaleTE = 100\n",
    "\n",
    "plt.figure(figsize=(15,7.5))\n",
    "plt.semilogx(default_l_hi, default_TT_hi, linestyle='solid', linewidth=3, label='Default, TT')\n",
    "plt.semilogx(hp_l_hi, hp_TT_hi, linestyle='dashed', label='HP, TT')\n",
    "plt.semilogx(default_l_hi, scaleEE*default_EE_hi, linestyle='solid', linewidth=3, label='Default, EE')\n",
    "plt.semilogx(hp_l_hi, scaleEE*hp_EE_hi, linestyle='dashed', label='HP, EE')\n",
    "plt.semilogx(default_l_hi, scaleTE*default_TE_hi, linestyle='solid', linewidth=3, label='Default, TE')\n",
    "plt.semilogx(hp_l_hi, scaleTE*hp_TE_hi, linestyle='dashed', label='HP, TE')\n",
    "plt.xlabel(r\"l\", fontsize=24)                                       \n",
    "plt.ylabel(r\"TT\", fontsize=24)\n",
    "plt.title(r'Power Spectrum (High)', fontsize=24)\n",
    "plt.legend(fontsize=18, loc='upper right')                                       \n",
    "plt.tick_params(axis='x', which='minttor')   \n",
    "plt.xlim(100, 5000)\n",
    "plt.grid(True, which='minor')                                                   \n",
    "#plt.savefig(prepath+'Plots/T'+str(temp)+'K/'+'Errors_CMBS4_Only_'+str(temp)+'K.png')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(15,7.5))\n",
    "plt.semilogx(default_l_hi, default_TT_hi-default_TT_lo, linestyle='solid', linewidth=3, label='Default, TT')\n",
    "plt.semilogx(hp_l_hi, hp_TT_hi-hp_TT_lo, linestyle='dashed', label='HP, TT')\n",
    "plt.semilogx(default_l_hi, scaleEE*(default_EE_hi-default_EE_lo), linestyle='solid', linewidth=3, label='Default, EE')\n",
    "plt.semilogx(hp_l_hi, scaleEE*(hp_EE_hi-hp_EE_lo), linestyle='dashed', label='HP, EE')\n",
    "plt.semilogx(default_l_hi, scaleTE*(default_TE_hi-default_TE_lo), linestyle='solid', linewidth=3, label='Default, TE')\n",
    "plt.semilogx(hp_l_hi, scaleTE*(hp_TE_hi-hp_TE_lo), linestyle='dashed', label='HP, TE')\n",
    "plt.xlabel(r\"l\", fontsize=24)                                       \n",
    "plt.ylabel(r\"TT\", fontsize=24)\n",
    "plt.title(r'Delta Power Spectum ', fontsize=24)\n",
    "plt.legend(fontsize=18, loc='upper right')                                       \n",
    "plt.tick_params(axis='x', which='minttor')   \n",
    "plt.xlim(100, 5000)\n",
    "plt.grid(True, which='minor')                                                   \n",
    "#plt.savefig(prepath+'Plots/T'+str(temp)+'K/'+'Errors_CMBS4_Only_'+str(temp)+'K.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Test 5: Testing Convergence "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Generate Convergence Forecasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-11T18:12:21.276566Z",
     "start_time": "2020-05-11T17:51:37.017825Z"
    },
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os                                                                       \n",
    "import shutil                                                                   \n",
    "import numpy as np                                                              \n",
    "import pandas as pd                                                             \n",
    "import seaborn as sns                                                           \n",
    "import cosmicfish as cf                                                         \n",
    "import matplotlib.pyplot as plt                                                 \n",
    "import dill                                                                     \n",
    "                                                                                \n",
    "                                                                                \n",
    "# Instruct pyplot to use seaborn                                                \n",
    "sns.set()                                                                       \n",
    "                                                                                \n",
    "# Set project, data, CLASS directories                                                                                    \n",
    "projectdir = '/Volumes/SSD01/results/test5/'                                    \n",
    "datastore = '/Volumes/SSD01/data.nosynctest5/'                                  \n",
    "classpath = '/Users/nicholasdeporzio/Desktop/cfworkspace/class/'                \n",
    "                                                                                \n",
    "# Specify resolution of numerical integrals                                     \n",
    "derivative_steps = [0.05, 0.1, 0.15] # 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 \n",
    "mu_integral_step = 1.00\n",
    "\n",
    "# Generate output paths                                                                \n",
    "ps23_resultsdir = projectdir                                                               \n",
    "cf.makedirectory(ps23_resultsdir)                                                                                        \n",
    "                                                                                \n",
    "# Linda Fiducial Cosmology                                                      \n",
    "ps23_fid = {                                                                    \n",
    "        \"A_s\" : 2.2321e-9,                                                      \n",
    "        \"n_s\" : 0.967,                                                          \n",
    "        \"omega_b\" : 0.02226,                                                    \n",
    "        \"omega_cdm\" : 0.1127,                                                   \n",
    "        \"tau_reio\" : 0.0598,                                                    \n",
    "        \"h\" : 0.701,                                                            \n",
    "        \"T_cmb\" : 2.726, # Units [K]                                            \n",
    "        \"N_ncdm\" : 4.,                                                          \n",
    "        \"deg_ncdm\" : 1.0,                                                       \n",
    "        \"T_ncdm\" : (0.79/2.726), # Units [T_cmb].                               \n",
    "        \"m_ncdm\" : 0.01, # Units [eV]                                           \n",
    "        \"b0\" : 1.0,                                                             \n",
    "        \"alphak2\" : 1.0,                                                        \n",
    "        \"sigma_fog_0\" : 250000, #Units [m s^-2]                                 \n",
    "        \"N_eff\" : 0.0064, #We allow relativistic neutrinos in addition to our DM relic\n",
    "        \"relic_vary\" : \"N_ncdm\",  # Fix T_ncdm or m_ncdm                        \n",
    "        \"m_nu\" : 0.02                                                           \n",
    "        }    \n",
    "\n",
    "# DESI Parameters                                                               \n",
    "#z_table = np.array([0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85])\n",
    "#dNdz = np.array([309., 2269., 1923., 2094., 1441., 1353., 1337., 523., 466., 329., 126., 0., 0.])\n",
    "z_table = np.array([0.65, 0.75])\n",
    "dNdz = np.array([309., 2269.])\n",
    "skycover = 14000. # Sky coverage of survey in degrees^2                         \n",
    "                                                                                \n",
    "# Run Fisher Forecast                                                                                               \n",
    "#masses = np.geomspace(0.01, 10., 21)                                                    \n",
    "#temps = np.array([0.79, 0.91, 0.94, 1.08])    \n",
    "masses = np.array([\n",
    "    np.geomspace(0.01, 10., 21)[0],\n",
    "    np.geomspace(0.01, 10., 21)[1],\n",
    "    np.geomspace(0.01, 10., 21)[13],\n",
    "    np.geomspace(0.01, 10., 21)[-2],\n",
    "    np.geomspace(0.01, 10., 21)[-1]])\n",
    "temps = np.array([0.79, 1.08])\n",
    "                                                                                \n",
    "omegacdm_set = np.array([                                                       \n",
    "    ps23_fid['omega_cdm']                                                       \n",
    "    - ((masses/cf.NEUTRINO_SCALE_FACTOR)* np.power(tval / 1.95, 3.))                                            \n",
    "    for tidx, tval in enumerate(temps)])    \n",
    "\n",
    "ps23_fiducialset = [[[                                                           \n",
    "    dict(ps23_fid, **{                                                          \n",
    "        'm_ncdm' : masses[midx],                                                \n",
    "        'omega_cdm' : omegacdm_set[tidx, midx],                                 \n",
    "        'T_ncdm' : temps[tidx]/2.726})                                          \n",
    "    for midx, mval in enumerate(masses)]                                        \n",
    "    for tidx, tval in enumerate(temps)]  \n",
    "    for didx, dval in enumerate(derivative_steps)]\n",
    "\n",
    "ps23_forecastset = [[[cf.forecast(                                               \n",
    "    classpath,                                                                  \n",
    "    datastore,                                                                  \n",
    "    '2relic',                                                                   \n",
    "    fidval,                                                                     \n",
    "    z_table,                                                                    \n",
    "    \"DESI\",                                                                     \n",
    "    dNdz,                                                                       \n",
    "    fcoverage_deg=skycover,                                                     \n",
    "    dstep=0.008,\n",
    "    gstep=derivative_steps[didx],\n",
    "    RSD=False,                                                                   \n",
    "    FOG=False,                                                                   \n",
    "    AP=False,                                                                    \n",
    "    COV=False)                                                                   \n",
    "    for fididx, fidval in enumerate(fidrowvals)]                                \n",
    "    for fidrowidx, fidrowvals in enumerate(dval)] \n",
    "    for didx, dval in enumerate(ps23_fiducialset)]\n",
    "\n",
    "#dill.load_session('')     \n",
    "for didx, dval in enumerate(ps23_forecastset): \n",
    "    for frowidx, frowval in enumerate(dval):                            \n",
    "        for fidx, fcst in enumerate(frowval):                                       \n",
    "            if type(fcst.fisher)==type(None):                                       \n",
    "                fcst.gen_pm()                                                       \n",
    "                fcst.gen_fisher(                                                    \n",
    "                    fisher_order=[                                                  \n",
    "                        'omega_b',                                                  \n",
    "                        'omega_cdm',                                                \n",
    "                        'n_s',                                                      \n",
    "                        'A_s',                                                      \n",
    "                        'tau_reio',                                                 \n",
    "                        'h',                                                        \n",
    "                        'N_ncdm',                                                   \n",
    "                        'M_ncdm',                                                   \n",
    "                        'sigma_fog',                                                \n",
    "                        'b0',                                                       \n",
    "                        'alpha_k2'],                                                \n",
    "                    mu_step=mu_integral_step,                                       \n",
    "                    skipgen=False)                                                  \n",
    "                print(\"Relic Forecast \", fidx, \" complete...\")                      \n",
    "                dill.dump_session(os.path.join(ps23_resultsdir, 'test5'+'_'+str(didx)+'_'+str(frowidx)+'_'+str(fidx)+'.db'))\n",
    "            else:                                                                   \n",
    "                print('Fisher matrix already generated!')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-11T14:01:29.958228Z",
     "start_time": "2020-05-11T14:01:29.952443Z"
    },
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Make Convergence Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import dill \n",
    "dill.load_session('/Volumes/SSD01/results/test5/test5_2_1_4.db')\n",
    "sns.set()\n",
    "sns.set_palette('bright')\n",
    "\n",
    "dsteps = [0.05, 0.1, 0.15]\n",
    "masses = np.array([\n",
    "    np.geomspace(0.01, 10., 21)[0],\n",
    "    np.geomspace(0.01, 10., 21)[1],\n",
    "    np.geomspace(0.01, 10., 21)[13],\n",
    "    np.geomspace(0.01, 10., 21)[-2],\n",
    "    np.geomspace(0.01, 10., 21)[-1]])\n",
    "temps = np.array([0.79, 1.08])\n",
    "\n",
    "var = 'N_ncdm'\n",
    "varfid = 1.0\n",
    "\n",
    "midx=4\n",
    "tidx=0\n",
    "\n",
    "k = [eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_high[0].k_table') for didx, dval in enumerate(dsteps)]\n",
    "pm_hi = [eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_high[0].ps_table') for didx, dval in enumerate(dsteps)]\n",
    "pb_hi = [eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_high[0].b_interp_table') for didx, dval in enumerate(dsteps)]\n",
    "pc_hi = [eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_high[0].cdm_interp_table') for didx, dval in enumerate(dsteps)]\n",
    "pm_lo = [eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_low[0].ps_table') for didx, dval in enumerate(dsteps)]\n",
    "pb_lo = [eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_low[0].b_interp_table') for didx, dval in enumerate(dsteps)]\n",
    "pc_lo = [eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_low[0].cdm_interp_table') for didx, dval in enumerate(dsteps)]\n",
    "dpm = [(pm_hi[didx]-pm_lo[didx])/(2.*dval*varfid) for didx, dval in enumerate(dsteps)]\n",
    "dpb = [(pb_hi[didx]-pb_lo[didx])/(2.*dval*varfid) for didx, dval in enumerate(dsteps)]\n",
    "dpc = [(pc_hi[didx]-pc_lo[didx])/(2.*dval*varfid) for didx, dval in enumerate(dsteps)]\n",
    "\n",
    "cl_l = [np.loadtxt(eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_high[0].datadirectory')+'/test_cl.dat', skiprows=12, usecols={0}) for didx, dval in enumerate(dsteps)]\n",
    "TT_hi = [np.loadtxt(eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_high[0].datadirectory')+'/test_cl.dat', skiprows=12, usecols={1}) for didx, dval in enumerate(dsteps)]\n",
    "TT_lo = [np.loadtxt(eval('ps23_forecastset['+str(didx)+']['+str(tidx)+']['+str(midx)+'].'+var+'_low[0].datadirectory')+'/test_cl.dat', skiprows=12, usecols={1}) for didx, dval in enumerate(dsteps)]\n",
    "dTT = [(TT_hi[didx]-TT_lo[didx])/(2.*dval*varfid) for didx, dval in enumerate(dsteps)]\n",
    "\n",
    "\n",
    "####\n",
    "# Pm \n",
    "####\n",
    "plt.figure(figsize=(15,7.5))\n",
    "for didx, dval in enumerate(dsteps): \n",
    "    plt.semilogx(k[didx], dpm[didx], linestyle='solid', linewidth=3, label=r'$\\Delta\\theta=$'+str(100.*dval)+r'%')\n",
    "plt.xlabel(r\"k\", fontsize=24)                                       \n",
    "plt.ylabel(r\"$dP_m/d\\theta$\", fontsize=24)\n",
    "plt.title(r'T='+str(temps[tidx])+', M='+str(masses[midx]),  fontsize=24)\n",
    "plt.legend(fontsize=18, loc='upper right')                                       \n",
    "plt.tick_params(axis='x', which='minttor')   \n",
    "#plt.xlim(100, 5000)\n",
    "plt.grid(True, which='minor')                                                   \n",
    "plt.savefig('/Users/nicholasdeporzio/Desktop/test/p1.png')\n",
    "plt.show()\n",
    "\n",
    "####\n",
    "# Pb \n",
    "####\n",
    "#plt.figure(figsize=(15,7.5))\n",
    "#for didx, dval in enumerate(dsteps): \n",
    "#    plt.semilogx(k[didx], dpb[didx], linestyle='solid', linewidth=3, label=r'$\\Delta\\theta=$'+str(100.*dval)+r'%')\n",
    "#plt.xlabel(r\"k\", fontsize=24)                                       \n",
    "#plt.ylabel(r\"$dP_b/d\\theta$\", fontsize=24)\n",
    "#plt.title(r'T='+str(temps[tidx])+', M='+str(masses[midx]),  fontsize=24)\n",
    "#plt.legend(fontsize=18, loc='upper right')                                       \n",
    "#plt.tick_params(axis='x', which='minttor')   \n",
    "##plt.xlim(100, 5000)\n",
    "#plt.grid(True, which='minor')                                                   \n",
    "##plt.savefig(prepath+'Plots/T'+str(temp)+'K/'+'Errors_CMBS4_Only_'+str(temp)+'K.png')\n",
    "#plt.show()\n",
    "\n",
    "####\n",
    "# Pc \n",
    "####\n",
    "#plt.figure(figsize=(15,7.5))\n",
    "#for didx, dval in enumerate(dsteps): \n",
    "#    plt.semilogx(k[didx], dpc[didx], linestyle='solid', linewidth=3, label=r'$\\Delta\\theta=$'+str(100.*dval)+r'%')\n",
    "#plt.xlabel(r\"k\", fontsize=24)                                       \n",
    "#plt.ylabel(r\"$dP_c/d\\theta$\", fontsize=24)\n",
    "#plt.title(r'T='+str(temps[tidx])+', M='+str(masses[midx]),  fontsize=24)\n",
    "#plt.legend(fontsize=18, loc='upper right')                                       \n",
    "#plt.tick_params(axis='x', which='minttor')   \n",
    "##plt.xlim(100, 5000)\n",
    "#plt.grid(True, which='minor')                                                   \n",
    "##plt.savefig(prepath+'Plots/T'+str(temp)+'K/'+'Errors_CMBS4_Only_'+str(temp)+'K.png')\n",
    "#plt.show()\n",
    "\n",
    "####\n",
    "# TT\n",
    "####\n",
    "plt.figure(figsize=(15,7.5))\n",
    "for didx, dval in enumerate(dsteps): \n",
    "    plt.semilogx(cl_l[didx], dTT[didx], linestyle='solid', linewidth=3, label=r'$\\Delta\\theta=$'+str(100.*dval)+r'%')\n",
    "plt.xlabel(r\"l\", fontsize=24)                                       \n",
    "plt.ylabel(r\"$dTT/d\\theta$\", fontsize=24)\n",
    "plt.title(r'T='+str(temps[tidx])+', M='+str(masses[midx]),  fontsize=24)\n",
    "plt.legend(fontsize=18, loc='upper right')                                       \n",
    "plt.tick_params(axis='x', which='minttor')   \n",
    "#plt.xlim(100, 5000)\n",
    "plt.grid(True, which='minor')                                                   \n",
    "plt.savefig('/Users/nicholasdeporzio/Desktop/test/p2.png')\n",
    "plt.show()\n",
    "\n",
    "####\n",
    "# deltaTT\n",
    "####\n",
    "plt.figure(figsize=(15,7.5))\n",
    "for didx, dval in enumerate(dsteps): \n",
    "    plt.semilogx(cl_l[didx], (dTT[didx]-dTT[-2])/dTT[-2], linestyle='solid', linewidth=3, label=r'$\\Delta\\theta=$'+str(100.*dval)+r'%')\n",
    "plt.xlabel(r\"l\", fontsize=24)                                       \n",
    "plt.ylabel(r\"$\\Delta(dTT/d\\theta)/(dTT/d\\theta_{10\\%})$\", fontsize=24)\n",
    "plt.title(r'T='+str(temps[tidx])+', M='+str(masses[midx]),  fontsize=24)\n",
    "plt.legend(fontsize=18, loc='upper right')                                       \n",
    "plt.tick_params(axis='x', which='minttor')   \n",
    "plt.ylim(-5, 5)\n",
    "plt.grid(True, which='minor')                                                   \n",
    "plt.savefig('/Users/nicholasdeporzio/Desktop/test/p3.png')\n",
    "plt.show()\n",
    "\n",
    "####\n",
    "# Pm \n",
    "####\n",
    "plt.figure(figsize=(15,7.5))\n",
    "for didx, dval in enumerate(dsteps): \n",
    "    plt.semilogx(k[didx], (dpm[didx]-dpm[-2])/dpm[-2], linestyle='solid', linewidth=3, label=r'$\\Delta\\theta=$'+str(100.*dval)+r'%')\n",
    "plt.xlabel(r\"k\", fontsize=24)                                       \n",
    "plt.ylabel(r\"$\\Delta(dP_m/d\\theta)/(dP_m/d\\theta_{10\\%})$\", fontsize=24)\n",
    "plt.title(r'T='+str(temps[tidx])+', M='+str(masses[midx]),  fontsize=24)\n",
    "plt.legend(fontsize=18, loc='upper right')                                       \n",
    "plt.tick_params(axis='x', which='minttor')   \n",
    "#plt.xlim(100, 5000)\n",
    "plt.grid(True, which='minor')                                                   \n",
    "plt.savefig('/Users/nicholasdeporzio/Desktop/test/p4.png')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Plot Contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T14:32:36.506798Z",
     "start_time": "2020-05-05T14:31:59.821801Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Show plots inline, and load main getdist plot module and samples class\n",
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "import sys, os\n",
    "sys.path.insert(0,os.path.realpath(os.path.join(os.getcwd(),'..')))\n",
    "from getdist import plots, MCSamples, loadMCSamples\n",
    "from getdist.gaussian_mixtures import GaussianND\n",
    "import getdist\n",
    "# use this *after* importing getdist if you want to use interactive plots\n",
    "# %matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "import IPython\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "print('GetDist Version: %s, Matplotlib version: %s'%(getdist.__version__, plt.matplotlib.__version__))\n",
    "# matplotlib 2 may not work very well without usetex on, can uncomment\n",
    "# plt.rcParams['text.usetex']=True\n",
    "\n",
    "#dataset = 'Full'\n",
    "dataset = 'LSS'\n",
    "\n",
    "#survey = 'CMBS4'\n",
    "survey = 'Planck'\n",
    "\n",
    "#Import fishers\n",
    "fisher_1 = np.loadtxt('/Volumes/SSD01/results/test1/'+survey+'/T0.91/M1_'+dataset+'_Fisher.dat', skiprows=1)\n",
    "fisher_2 = np.loadtxt('/Volumes/SSD01/results/test2/'+survey+'/'+dataset+'_Fisher.dat', skiprows=1)\n",
    "fisher_3 = np.loadtxt('/Volumes/SSD01/results/test3/'+survey+'/T0.91/M1_'+dataset+'_Fisher.dat', skiprows=1)\n",
    "\n",
    "fisher_1 = np.loadtxt('/Volumes/SSD01/results/test1/'+survey+'/T0.91/M1_'+dataset+'_Fisher.dat', skiprows=1)\n",
    "fisher_2 = np.loadtxt('/Volumes/SSD01/results/test2/'+survey+'/'+dataset+'_Fisher.dat', skiprows=1)\n",
    "fisher_3 = np.loadtxt('/Volumes/SSD01/results/test3/'+survey+'/T0.91/M1_'+dataset+'_Fisher.dat', skiprows=1)\n",
    "\n",
    "#Delete tau and invert for covariances\n",
    "cov_1 = np.linalg.inv(np.delete(np.delete(fisher_1, 5, axis=0), 5, axis=1)) \n",
    "cov_2 = np.linalg.inv(np.delete(np.delete(fisher_2, 4, axis=0), 4, axis=1)) \n",
    "cov_3 = np.linalg.inv(np.delete(np.delete(fisher_3, 5, axis=0), 5, axis=1)) \n",
    "\n",
    "fid_1 = [\n",
    "    2.2321e-9, #A_s\n",
    "    0.701, #h\n",
    "    0.967, #n_s\n",
    "    0.02226, #omega_b\n",
    "    0.1127, #omega_cdm\n",
    "#    0.0598, #tau_reio\n",
    "    1.0, #N_ncdm \n",
    "    250., #sigma_fog\n",
    "    1.0, #beta0\n",
    "    1.0, #alphak2 \n",
    "]\n",
    "fid_2 = [\n",
    "    0.02226, #omega_b\n",
    "    0.1127, #omega_cdm\n",
    "    0.967, #n_s\n",
    "    2.2321e-9, #A_s\n",
    "#    0.0598, #tau_reio\n",
    "    0.701, #h\n",
    "    0.06, #M_ncdm \n",
    "    250., #sigma_fog\n",
    "    1.0, #beta0\n",
    "    1.0, #alphak2 \n",
    "]\n",
    "fid_3 = [\n",
    "    2.2321e-9, #A_s\n",
    "    0.701, #h\n",
    "    0.967, #n_s\n",
    "    0.02226, #omega_b\n",
    "    0.1127, #omega_cdm\n",
    "#    0.0598, #tau_reio\n",
    "    1.0, #N_ncdm\n",
    "    0.06, #M_ncdm\n",
    "    250., #sigma_fog\n",
    "    1.0, #beta0\n",
    "    1.0, #alphak2 \n",
    "]\n",
    "\n",
    "\n",
    "\n",
    "params_1 = [\n",
    "    r'$A_s$', \n",
    "    r'$h$', \n",
    "    r'$n_s$', \n",
    "    r'$\\omega_b$', \n",
    "    r'$\\omega_{cdm}$', \n",
    "#    r'$\\tau_{reio}$', \n",
    "    r'$g_\\chi$', \n",
    "    r'$\\sigma_{fog}$', \n",
    "    r'$\\beta_0$', \n",
    "    r'$\\alpha_{k2}$']\n",
    "params_2 = [\n",
    "    r'$\\omega_b$', \n",
    "    r'$\\omega_{cdm}$', \n",
    "    r'$n_s$', \n",
    "    r'$A_s$',\n",
    "#    r'$\\tau_{reio}$',\n",
    "    r'$h$',  \n",
    "    r'$M_\\nu$', \n",
    "    r'$\\sigma_{fog}$', \n",
    "    r'$\\beta_0$', \n",
    "    r'$\\alpha_{k2}$']\n",
    "params_3 = [\n",
    "    r'$A_s$', \n",
    "    r'$h$', \n",
    "    r'$n_s$', \n",
    "    r'$\\omega_b$', \n",
    "    r'$\\omega_{cdm}$', \n",
    "#    r'$\\tau_{reio}$', \n",
    "    r'$g_\\chi$', \n",
    "    r'$M_\\nu$', \n",
    "    r'$\\sigma_{fog}$', \n",
    "    r'$\\beta_0$', \n",
    "    r'$\\alpha_{k2}$']\n",
    "\n",
    "gauss_1 = GaussianND(fid_1, cov_1, names=params_1)\n",
    "gauss_2 = GaussianND(fid_2, cov_2, names=params_2)\n",
    "gauss_3 = GaussianND(fid_3, cov_3, names=params_3)\n",
    "\n",
    "h = plots.get_subplot_plotter()\n",
    "\n",
    "# Set a default plot style \n",
    "from getdist.styles.tab10 import style_name\n",
    "plots.set_active_style(style_name)\n",
    "\n",
    "# Adjust plot sizes, styles, etc using getdist.plots.GetDistPlotSettings class\n",
    "h.settings.figure_legend_frame = False #Box around legend\n",
    "h.settings.scaling=False #Defaults to font size for Tick Labels: 11, Labels: 14, Legend: 12, Text: 12\n",
    "h.settings.axes_fontsize=14\n",
    "h.settings.axes_labelsize=20\n",
    "h.settings.axis_tick_x_rotation=45\n",
    "h.settings.legend_fontsize=20\n",
    "#h.settings.axis_tick_y_rotation=90\n",
    "#h.settings.colorbar_tick_rotation=90\n",
    "if dataset=='Full': \n",
    "    if survey=='CMBS4':\n",
    "        legendposttext = ', DESI+CMBS4'\n",
    "    elif survey=='Planck':\n",
    "        legendposttext = ', DESI+Planck'\n",
    "elif dataset=='LSS': \n",
    "    legendposttext = ', DESI Only'\n",
    "else: \n",
    "    legendposttext = ''\n",
    "h.triangle_plot(\n",
    "    [\n",
    "        gauss_1, \n",
    "        #gauss_2, \n",
    "        gauss_3\n",
    "    ], \n",
    "    filled=True,\n",
    "    legend_labels=[\n",
    "        r'$g_\\chi$-Only'+legendposttext, \n",
    "        #r'$M_\\nu$ -Only'+legendposttext, \n",
    "        r'$g_\\chi$ & $M_\\nu$ Simultaneous' + legendposttext\n",
    "    ], \n",
    "    line_args=[\n",
    "        {'ls':'solid', 'lw':2, 'color':'#A51C30'},\n",
    "        #{'ls':'solid', 'lw':2, 'color':'#93a1ad'}, \n",
    "        {'ls':'solid', 'lw':2, 'color':'blue'}\n",
    "    ], \n",
    "    contour_args=[\n",
    "        {'ls':'dotted', 'lw':2, 'color':'#A51C30', 'filled' : True},\n",
    "        #{'ls':'dotted', 'lw':2, 'color':'#93a1ad', 'filled' : True}, \n",
    "        {'ls':'dotted', 'lw':2, 'color':'blue', 'filled' : False}\n",
    "    ])\n",
    "plt.savefig(\"/Users/nicholasdeporzio/Desktop/diagnostics/comparison.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Add Priors to Cluster Output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T13:58:27.329068Z",
     "start_time": "2020-06-04T13:48:57.952738Z"
    },
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import dill\n",
    "import numpy as np\n",
    "\n",
    "#dtype=['Fine', 'Grid']\n",
    "dtype=['Fine']\n",
    "#lsssurvey=['BOSS', 'EUCLID', 'DESI']\n",
    "lsssurvey=['DESI']\n",
    "cmbsurvey=['CMBS4', 'Planck']\n",
    "\n",
    "for idx1, val1 in enumerate(dtype):\n",
    "    for idx2, val2 in enumerate(lsssurvey):\n",
    "        for idx3, val3 in enumerate(cmbsurvey):\n",
    "            prepath_in='/Volumes/SSD01/TEST_RESULTS_1/'+val2+'/Forecast_Files/'+val2+'_'+val1+'Plot_'\n",
    "            prepath_out='/Volumes/SSD01/TEST_RESULTS_1/'+val2+'/Fisher_Matrices/'+val1+'/T'\n",
    "            \n",
    "            if val1=='Grid': \n",
    "                nmax=120\n",
    "                modnum=10\n",
    "                fishprefix='gp'\n",
    "            \n",
    "            elif val1=='Fine':\n",
    "                nmax=84\n",
    "                modnum=21\n",
    "                fishprefix='fp'\n",
    "                masses = np.geomspace(0.01, 10., 21)\n",
    "                temps = np.array([0.79, 0.91, 0.94, 1.08])\n",
    "            \n",
    "            \n",
    "            for idx in range(nmax): \n",
    "            \n",
    "                temp_index=(idx // modnum)\n",
    "                mass_index=(idx % modnum)\n",
    "            \n",
    "                dill.load_session(prepath_in+str(idx)+'/'+fishprefix+'_'+str(temp_index)+'_'+str(mass_index)+'.db')  \n",
    "                \n",
    "                if val1=='Grid':\n",
    "                    masses = np.geomspace(0.01, 10., 10)\n",
    "                    temps = np.linspace(0.79, 1.5, 12)\n",
    "                elif val1=='Fine':\n",
    "                    masses = np.geomspace(0.01, 10., 21)\n",
    "                    temps = np.array([0.79, 0.91, 0.94, 1.08])\n",
    "                \n",
    "                #print(eval(fishprefix+'_forecastset')[0][0].fisher)\n",
    "            \n",
    "                mass = masses[mass_index]\n",
    "                temp = temps[temp_index]\n",
    "                \n",
    "                outdir = prepath_out+f'{temp:.2f}'+'K/M'+f'{mass:.3f}'+'/'+val3+'/'\n",
    "                cf.makedirectory(outdir)\n",
    "            \n",
    "                for tidx, tval in enumerate(eval(fishprefix+'_forecastset')):\n",
    "                    for midx, mval in enumerate(tval): \n",
    "                        if val2=='DESI': \n",
    "                            mval.reorder_fisher(['A_s', 'h', 'n_s', 'omega_b', 'omega_cdm', 'tau_reio', 'M_ncdm', 'N_ncdm', 'sigma_fog', 'b0', 'alpha_k2'])\n",
    "                        else: \n",
    "                            mval.reorder_fisher(['A_s', 'h', 'n_s', 'omega_b', 'omega_cdm', 'tau_reio', 'M_ncdm', 'N_ncdm', 'sigma_fog', 'beta0', 'beta1', 'alpha_k2'])\n",
    "                            \n",
    "                        mval.load_cmb_fisher(\n",
    "                            fisher_order=[\n",
    "                                'A_s',\n",
    "                                'h',  \n",
    "                                'n_s',\n",
    "                                'omega_b',                                    \n",
    "                                'omega_cdm',                                                     \n",
    "                                'tau_reio',                                                                                    \n",
    "                                'M_ncdm', \n",
    "                                'N_ncdm'],\n",
    "                            fisherpath = cf.priors_directory()+'/2relic_Priors/'+val1+'/Fisher_T'+f'{temp:.2f}'+'K/Fisher'+val3+'_bin'+str(mass_index+1)+'_T'+f'{temp:.2f}'+'K.txt'\n",
    "                        )\n",
    "                        mval.export_matrices(outdir)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-15T13:00:57.082746Z",
     "start_time": "2020-05-15T13:00:57.054378Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-15T13:01:00.532728Z",
     "start_time": "2020-05-15T13:00:59.930509Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-15T13:02:07.001191Z",
     "start_time": "2020-05-15T13:02:06.995102Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-15T13:38:04.766298Z",
     "start_time": "2020-05-15T13:38:04.762485Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-05T12:41:47.295174Z",
     "start_time": "2020-06-05T12:41:47.287788Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "df2 = pd.DataFrame(np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]), columns=['a', 'b', 'c'])\n",
    "\n",
    "test = np.array([df2.values, 10*df2.values])\n",
    "np.shape(test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Plot Errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Forecast Class Definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-08T07:58:07.403449Z",
     "start_time": "2020-06-08T07:58:06.590529Z"
    },
    "code_folding": [],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np \n",
    "from os.path import join as jn\n",
    "\n",
    "class parameter:    \n",
    "    def __init__(self, name, temp_vals=0, mass_vals=0):\n",
    "        self.name = name\n",
    "        self.fixed = False \n",
    "\n",
    "class forecast: \n",
    "    '''Forecasts Fisher/Covariance matrices from fid cosmo.'''\n",
    "    \n",
    "    def __init__(self, \n",
    "                 lss_survey_name, \n",
    "                 fisher_matrices_dir, \n",
    "                 masses, \n",
    "                 temps):\n",
    "        \n",
    "        self.lss_survey_name = lss_survey_name\n",
    "        self.fisher_matrices_dir = fisher_matrices_dir\n",
    "        self.masses = masses\n",
    "        self.temps = temps\n",
    "        \n",
    "        self.prepath = jn(fisher_matrices_dir, 'Fine', 'T')\n",
    "        \n",
    "        self.paths_lss = [[jn(self.prepath+str(tval)+'K', 'M'+f'{mval:.3f}', 'Planck', 'inv_lssfisher.mat') \n",
    "            for midx, mval in enumerate(masses)] \n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.paths_Planck = [[jn(self.prepath+str(tval)+'K', 'M'+f'{mval:.3f}', 'Planck', 'inv_cmbfisher.mat') \n",
    "            for midx, mval in enumerate(masses)] \n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.paths_CMBS4 = [[jn(self.prepath+str(tval)+'K', 'M'+f'{mval:.3f}', 'CMBS4', 'inv_cmbfisher.mat') \n",
    "            for midx, mval in enumerate(masses)] \n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.paths_Full_Planck = [[jn(self.prepath+str(tval)+'K', 'M'+f'{mval:.3f}', 'Planck', 'inv_fullfisher.mat') \n",
    "            for midx, mval in enumerate(masses)] \n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.paths_Full_CMBS4 = [[jn(self.prepath+str(tval)+'K', 'M'+f'{mval:.3f}', 'CMBS4', 'inv_fullfisher.mat') \n",
    "            for midx, mval in enumerate(masses)] \n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        \n",
    "        self.covs_lss = [[pd.read_csv(self.paths_lss[tidx][midx], sep='\\t', header=0)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.covs_Planck = [[pd.read_csv(self.paths_Planck[tidx][midx], sep='\\t', header=0)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.covs_CMBS4 = [[pd.read_csv(self.paths_CMBS4[tidx][midx], sep='\\t', header=0)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.covs_Full_Planck = [[pd.read_csv(self.paths_Full_Planck[tidx][midx], sep='\\t', header=0)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.covs_Full_CMBS4 = [[pd.read_csv(self.paths_Full_CMBS4[tidx][midx], sep='\\t', header=0)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        \n",
    "        self.column_names = self.covs_Full_Planck[0][0].columns\n",
    "        self.parameters = [parameter(mval) for midx, mval in enumerate(list(self.column_names))]\n",
    "\n",
    "        self.fishers_lss = [[pd.DataFrame(np.linalg.inv(self.covs_lss[tidx][midx].values), columns=self.covs_lss[tidx][midx].columns)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.fishers_Planck = [[pd.DataFrame(np.linalg.inv(self.covs_Planck[tidx][midx].values), columns=self.covs_Planck[tidx][midx].columns)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.fishers_CMBS4 = [[pd.DataFrame(np.linalg.inv(self.covs_CMBS4[tidx][midx].values), columns=self.covs_CMBS4[tidx][midx].columns)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.fishers_Full_Planck = [[pd.DataFrame(np.linalg.inv(self.covs_Full_Planck[tidx][midx].values), columns=self.covs_Full_Planck[tidx][midx].columns)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.fishers_Full_CMBS4 = [[pd.DataFrame(np.linalg.inv(self.covs_Full_CMBS4[tidx][midx].values), columns=self.covs_Full_CMBS4[tidx][midx].columns)\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "\n",
    "        self.clear()\n",
    "                \n",
    "    def clear(self): \n",
    "        for pidx, pval in enumerate(self.parameters): \n",
    "            pval.fixed=False \n",
    "            \n",
    "        self.final_covs_lss = [[pd.DataFrame(self.covs_lss[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.final_covs_Planck = [[pd.DataFrame(self.covs_Planck[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.final_covs_CMBS4 = [[pd.DataFrame(self.covs_CMBS4[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.final_covs_Full_Planck = [[pd.DataFrame(self.covs_Full_Planck[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.final_covs_Full_CMBS4 = [[pd.DataFrame(self.covs_Full_CMBS4[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "\n",
    "        self.final_fishers_lss = [[pd.DataFrame(self.fishers_lss[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.final_fishers_Planck = [[pd.DataFrame(self.fishers_Planck[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.final_fishers_CMBS4 = [[pd.DataFrame(self.fishers_CMBS4[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.final_fishers_Full_Planck = [[pd.DataFrame(self.fishers_Full_Planck[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        self.final_fishers_Full_CMBS4 = [[pd.DataFrame(self.fishers_Full_CMBS4[tidx][midx])\n",
    "            for midx, mval in enumerate(masses)]\n",
    "            for tidx, tval in enumerate(temps)]\n",
    "        \n",
    "    def fix_param(self, name): \n",
    "        flag=False\n",
    "        for pidx, pval in enumerate(self.parameters):\n",
    "            if pval.name==name:\n",
    "                pval.fixed=True\n",
    "                flag=True\n",
    "        if flag==True: \n",
    "            self.process()\n",
    "            print(name, ' has been fixed...')\n",
    "        else: \n",
    "            print(name, ' is not a valid parameter name. Try again... ')\n",
    "\n",
    "        \n",
    "    def process(self):   \n",
    "        \n",
    "        for tidx, tval in enumerate(temps):\n",
    "            for midx, mval in enumerate(masses): \n",
    "                for pidx, pval in enumerate(self.parameters):               \n",
    "                    \n",
    "                    start_params = self.final_fishers_Full_CMBS4[tidx][midx].columns\n",
    "                    if (pval.fixed==True) and (pval.name in start_params): \n",
    "                        \n",
    "                        # Process Fishers\n",
    "                        if (pval.name in self.final_fishers_lss[tidx][midx].columns): \n",
    "                            lss_drop_index = self.final_fishers_lss[tidx][midx].columns.get_loc(pval.name)\n",
    "                            self.final_fishers_lss[tidx][midx].drop(columns=self.final_fishers_lss[tidx][midx].columns[lss_drop_index], inplace=True)\n",
    "                            self.final_fishers_lss[tidx][midx].drop(index=[lss_drop_index], inplace=True)\n",
    "                            self.final_fishers_lss[tidx][midx].reset_index(drop=True, inplace=True)\n",
    "\n",
    "                        if (pval.name in self.final_fishers_Planck[tidx][midx].columns):  \n",
    "                            cmb_drop_index = self.final_fishers_Planck[tidx][midx].columns.get_loc(pval.name)\n",
    "                            self.final_fishers_Planck[tidx][midx].drop(columns=self.final_fishers_Planck[tidx][midx].columns[cmb_drop_index], inplace=True)\n",
    "                            self.final_fishers_Planck[tidx][midx].drop(index=[cmb_drop_index], inplace=True)\n",
    "                            self.final_fishers_Planck[tidx][midx].reset_index(drop=True, inplace=True)\n",
    "                            \n",
    "                            self.final_fishers_CMBS4[tidx][midx].drop(columns=self.final_fishers_CMBS4[tidx][midx].columns[cmb_drop_index], inplace=True)\n",
    "                            self.final_fishers_CMBS4[tidx][midx].drop(index=[cmb_drop_index], inplace=True)\n",
    "                            self.final_fishers_CMBS4[tidx][midx].reset_index(drop=True, inplace=True)\n",
    "                            \n",
    "                        full_drop_index = self.final_fishers_Full_Planck[tidx][midx].columns.get_loc(pval.name)\n",
    "                        self.final_fishers_Full_Planck[tidx][midx].drop(columns=self.final_fishers_Full_Planck[tidx][midx].columns[full_drop_index], inplace=True)\n",
    "                        self.final_fishers_Full_Planck[tidx][midx].drop(index=[full_drop_index], inplace=True)\n",
    "                        self.final_fishers_Full_Planck[tidx][midx].reset_index(drop=True, inplace=True)\n",
    "                            \n",
    "                        self.final_fishers_Full_CMBS4[tidx][midx].drop(columns=self.final_fishers_Full_CMBS4[tidx][midx].columns[full_drop_index], inplace=True)\n",
    "                        self.final_fishers_Full_CMBS4[tidx][midx].drop(index=[full_drop_index], inplace=True)\n",
    "                        self.final_fishers_Full_CMBS4[tidx][midx].reset_index(drop=True, inplace=True)\n",
    "                        \n",
    "                        \n",
    "                        # Process Covariances  \n",
    "                        self.final_covs_lss[tidx][midx]=pd.DataFrame(np.linalg.inv(self.final_fishers_lss[tidx][midx].values), columns=self.final_fishers_lss[tidx][midx].columns) \n",
    "                        \n",
    "                        self.final_covs_Planck[tidx][midx]=pd.DataFrame(np.linalg.inv(self.final_fishers_Planck[tidx][midx].values), columns=self.final_fishers_Planck[tidx][midx].columns)\n",
    "                        \n",
    "                        self.final_covs_CMBS4[tidx][midx]=pd.DataFrame(np.linalg.inv(self.final_fishers_CMBS4[tidx][midx].values), columns=self.final_fishers_CMBS4[tidx][midx].columns)\n",
    "                        \n",
    "                        self.final_covs_Full_Planck[tidx][midx]=pd.DataFrame(np.linalg.inv(self.final_fishers_Full_Planck[tidx][midx].values), columns=self.final_fishers_Full_Planck[tidx][midx].columns)\n",
    "                        \n",
    "                        self.final_covs_Full_CMBS4[tidx][midx]=pd.DataFrame(np.linalg.inv(self.final_fishers_Full_CMBS4[tidx][midx].values), columns=self.final_fishers_Full_CMBS4[tidx][midx].columns)\n",
    "                        \n",
    "                                                         \n",
    "                                                         \n",
    "    def print_fixed(self): \n",
    "        print('The following parameters are currently fixed: ')\n",
    "        for pidx, pval in enumerate(self.parameters): \n",
    "            if pval.fixed==True: \n",
    "                print(pval.name, '\\n')\n",
    "        print('Execute `forecast.clear()` to marginalize... ')\n",
    "                                             "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Validate that DESI-As-Euclid Matches Euclid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-05T17:22:47.872160Z",
     "start_time": "2020-06-05T17:22:28.958842Z"
    },
    "hidden": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import shutil\n",
    "import seaborn as sns\n",
    "import cosmicfish as cf \n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# Instruct pyplot to use seaborn \n",
    "sns.set()\n",
    "\n",
    "masses = np.geomspace(0.01, 10., 21)\n",
    "temps = np.array([0.79, 0.91, 0.94, 1.08])\n",
    "param = 'N_ncdm'\n",
    "\n",
    "boss_forecast = forecast('BOSS', '/Volumes/SSD01/TEST_RESULTS_1/BOSS/Fisher_Matrices/', masses, temps)\n",
    "print('BOSS Loaded...')\n",
    "\n",
    "desi_forecast = forecast('DESI', '/Volumes/SSD01/TEST_RESULTS_1/DESI/Fisher_Matrices/', masses, temps)\n",
    "print('DESI Loaded...')\n",
    "\n",
    "euclid_forecast = forecast('EUCLID', '/Volumes/SSD01/TEST_RESULTS_1/EUCLID/Fisher_Matrices/', masses, temps)\n",
    "print('Euclid Loaded...')\n",
    "\n",
    "\n",
    "boss_forecast.clear()\n",
    "desi_forecast.clear()\n",
    "euclid_forecast.clear()\n",
    "\n",
    "boss_forecast.fix_param('tau_reio')\n",
    "desi_forecast.fix_param('tau_reio')\n",
    "euclid_forecast.fix_param('tau_reio')\n",
    "\n",
    "boss_forecast.fix_param('beta1')\n",
    "euclid_forecast.fix_param('beta1')\n",
    "\n",
    "p1 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p2 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p3 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p4 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p5 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    for midx, mval in enumerate(masses): \n",
    "        p1[tidx][midx]=(desi_forecast.final_fishers_lss[tidx][midx].values-euclid_forecast.final_fishers_lss[tidx][midx].values)/euclid_forecast.final_fishers_lss[tidx][midx].values\n",
    "        p2[tidx][midx]=(desi_forecast.final_covs_lss[tidx][midx].values-euclid_forecast.final_covs_lss[tidx][midx].values)/euclid_forecast.final_covs_lss[tidx][midx].values\n",
    "        p3[tidx][midx]=(desi_forecast.final_covs_lss[tidx][midx][param][desi_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)]-euclid_forecast.final_covs_lss[tidx][midx][param][euclid_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])/euclid_forecast.final_covs_lss[tidx][midx][param][euclid_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)]\n",
    "        p4[tidx][midx]=np.sqrt(desi_forecast.final_covs_lss[tidx][midx][param][desi_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "        p5[tidx][midx]=np.sqrt(euclid_forecast.final_covs_lss[tidx][midx][param][euclid_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "p1 = np.array(p1)\n",
    "p2 = np.array(p2)\n",
    "p3 = np.array(p3)\n",
    "p4 = np.array(p4)\n",
    "p5 = np.array(p5)\n",
    "\n",
    "        \n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    for idx1 in range(len(p1[0, 0, 0, :])):\n",
    "        for idx2 in range(len(p1[0, 0, 0, :])):\n",
    "            plt.plot(masses, p1[tidx, :, idx1, idx2], linestyle='solid')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$((F_{ij})_{DESI}-(F_{ij})_{Euclid})/(F_{ij})_{Euclid}$\", fontsize=24)\n",
    "    plt.title('Fisher Relative Error, T='+str(tval)+'K', fontsize=24)\n",
    "    #plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    #plt.savefig('/Users/nicholasdeporzio/Desktop/06.05.2020/F_RelErr_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    for idx1 in range(len(p2[0, 0, 0, :])):\n",
    "        for idx2 in range(len(p2[0, 0, 0, :])):\n",
    "            plt.plot(masses, p2[tidx, :, idx1, idx2], linestyle='solid')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$((C_{ij})_{DESI}-(C_{ij})_{Euclid})/(C_{ij})_{Euclid}$\", fontsize=24)\n",
    "    plt.title('Covariance Relative Error, T='+str(tval)+'K', fontsize=24)\n",
    "    #plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    #plt.savefig('/Users/nicholasdeporzio/Desktop/06.05.2020/C_RelErr_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.plot(masses, p3[tidx, :], linestyle='solid')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$((C_{g})_{DESI}-(C_{g})_{Euclid})/(C_{g})_{Euclid}$\", fontsize=24)\n",
    "    plt.title(r'$C(g_X)$ Relative Error, T='+str(tval)+'K', fontsize=24)\n",
    "    #plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    #plt.savefig('/Users/nicholasdeporzio/Desktop/06.05.2020/Cg_RelErr_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, p4[tidx, :], linestyle='solid', linewidth=4, label='DESI')\n",
    "    plt.semilogx(masses, p5[tidx, :], linestyle='dashed', linewidth=4, label='EUCLID')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\\sigma(g_X)$\", fontsize=24)\n",
    "    plt.title(r'Error, T='+str(tval)+'K', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    #plt.savefig('/Users/nicholasdeporzio/Desktop/06.05.2020/Cg_RelErr_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "desi_only_gErrs_tauFixed_beta1Fixed_MnuMarged = np.array(p4)\n",
    "euclid_only_gErrs_tauFixed_beta1Fixed_MnuMarged = np.array(p5)\n",
    "\n",
    "print('Beta1 fixed, M marginalized, Euclid covariances: ')\n",
    "#print(desi_forecast.final_covs_lss[0][0])\n",
    "print(euclid_forecast.final_covs_lss[0][0])\n",
    "\n",
    "####################\n",
    "\n",
    "boss_forecast.clear()\n",
    "desi_forecast.clear()\n",
    "euclid_forecast.clear()\n",
    "\n",
    "boss_forecast.fix_param('tau_reio')\n",
    "desi_forecast.fix_param('tau_reio')\n",
    "euclid_forecast.fix_param('tau_reio')\n",
    "\n",
    "boss_forecast.fix_param('beta1')\n",
    "euclid_forecast.fix_param('beta1')\n",
    "\n",
    "boss_forecast.fix_param('M_ncdm')\n",
    "desi_forecast.fix_param('M_ncdm')\n",
    "euclid_forecast.fix_param('M_ncdm')\n",
    "\n",
    "p6 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p7 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    for midx, mval in enumerate(masses): \n",
    "        p6[tidx][midx]=np.sqrt(desi_forecast.final_covs_lss[tidx][midx][param][desi_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "        p7[tidx][midx]=np.sqrt(euclid_forecast.final_covs_lss[tidx][midx][param][euclid_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "p6 = np.array(p6)\n",
    "p7 = np.array(p7)\n",
    "\n",
    "desi_only_gErrs_tauFixed_beta1Fixed_MnuFixed = np.array(p6)\n",
    "euclid_only_gErrs_tauFixed_beta1Fixed_MnuFixed = np.array(p7)\n",
    "\n",
    "print('Beta1 fixed, M Fixed, Euclid covariances: ')\n",
    "#print(desi_forecast.final_covs_lss[0][0])\n",
    "print(euclid_forecast.final_covs_lss[0][0])\n",
    "\n",
    "##############\n",
    "\n",
    "boss_forecast.clear()\n",
    "desi_forecast.clear()\n",
    "euclid_forecast.clear()\n",
    "\n",
    "boss_forecast.fix_param('tau_reio')\n",
    "desi_forecast.fix_param('tau_reio')\n",
    "euclid_forecast.fix_param('tau_reio')\n",
    "\n",
    "\n",
    "p8 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p9 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    for midx, mval in enumerate(masses): \n",
    "        p8[tidx][midx]=np.sqrt(desi_forecast.final_covs_lss[tidx][midx][param][desi_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "        p9[tidx][midx]=np.sqrt(euclid_forecast.final_covs_lss[tidx][midx][param][euclid_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "p8 = np.array(p8)\n",
    "p9 = np.array(p9)\n",
    "\n",
    "desi_only_gErrs_tauFixed_beta1Marged_MnuMarged = np.array(p8)\n",
    "euclid_only_gErrs_tauFixed_beta1Marged_MnuMarged = np.array(p9)\n",
    "\n",
    "print('Beta1 Marginalized, M Marginalized, Euclid covariances: ')\n",
    "#print(desi_forecast.final_covs_lss[0][0])\n",
    "print(euclid_forecast.final_covs_lss[0][0])\n",
    "\n",
    "##############\n",
    "\n",
    "boss_forecast.clear()\n",
    "desi_forecast.clear()\n",
    "euclid_forecast.clear()\n",
    "\n",
    "boss_forecast.fix_param('tau_reio')\n",
    "desi_forecast.fix_param('tau_reio')\n",
    "euclid_forecast.fix_param('tau_reio')\n",
    "\n",
    "boss_forecast.fix_param('M_ncdm')\n",
    "desi_forecast.fix_param('M_ncdm')\n",
    "euclid_forecast.fix_param('M_ncdm')\n",
    "\n",
    "\n",
    "p10 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p11 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    for midx, mval in enumerate(masses): \n",
    "        p10[tidx][midx]=np.sqrt(desi_forecast.final_covs_lss[tidx][midx][param][desi_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "        p11[tidx][midx]=np.sqrt(euclid_forecast.final_covs_lss[tidx][midx][param][euclid_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "p10 = np.array(p10)\n",
    "p11 = np.array(p11)\n",
    "\n",
    "desi_only_gErrs_tauFixed_beta1Marged_MnuFixed = np.array(p10)\n",
    "euclid_only_gErrs_tauFixed_beta1Marged_MnuFixed = np.array(p11)\n",
    "\n",
    "print('Beta1 Marginalized, M Fixed, Euclid covariances: ')\n",
    "#print(desi_forecast.final_covs_lss[0][0])\n",
    "print(euclid_forecast.final_covs_lss[0][0])\n",
    "\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, euclid_only_gErrs_tauFixed_beta1Fixed_MnuMarged[tidx], linestyle='solid', color='green', label='Beta1 Fixed, M Marged')\n",
    "    plt.semilogx(masses, euclid_only_gErrs_tauFixed_beta1Fixed_MnuFixed[tidx], linestyle='solid', color='red', label='Beta1 Fixed, M Fixed')\n",
    "    plt.loglog(masses, euclid_only_gErrs_tauFixed_beta1Marged_MnuMarged[tidx], linestyle='dashed', linewidth=2, color='green', label='Beta1 Marged, M Marged')\n",
    "    plt.loglog(masses, euclid_only_gErrs_tauFixed_beta1Marged_MnuFixed[tidx], linestyle='dashed', linewidth=2, color='red', label='Beta1 Marged, M Fixed')\n",
    "\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\", fontsize=24)\n",
    "    plt.title('', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    #plt.savefig('/Users/nicholasdeporzio/Desktop/06.05.2020/F_RelErr_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "\n",
    "print(euclid_only_gErrs_tauFixed_beta1Marged_MnuMarged[0][0])\n",
    "print(euclid_only_gErrs_tauFixed_beta1Marged_MnuFixed[0][0])\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-05T12:54:20.830887Z",
     "start_time": "2020-06-05T12:54:20.825827Z"
    },
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Generate Paper Error Table "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-08T07:58:55.117860Z",
     "start_time": "2020-06-08T07:58:10.718012Z"
    },
    "hidden": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import shutil\n",
    "import seaborn as sns\n",
    "import cosmicfish as cf \n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# Instruct pyplot to use seaborn \n",
    "sns.set()\n",
    "\n",
    "masses = np.geomspace(0.01, 10., 21)\n",
    "temps = np.array([0.79, 0.91, 0.94, 1.08])\n",
    "\n",
    "\n",
    "boss_forecast = forecast('BOSS', '/Volumes/SSD01/FINAL_RESULTS_NewBias/BOSS/Fisher_Matrices/', masses, temps)\n",
    "print('BOSS Loaded...')\n",
    "desi_forecast = forecast('DESI', '/Volumes/SSD01/FINAL_RESULTS_NewBias/DESI/Fisher_Matrices/', masses, temps)\n",
    "print('DESI Loaded...')\n",
    "euclid_forecast = forecast('EUCLID', '/Volumes/SSD01/FINAL_RESULTS_NewBias/EUCLID/Fisher_Matrices/', masses, temps)\n",
    "print('Euclid Loaded...')\n",
    "\n",
    "param = 'N_ncdm'\n",
    "\n",
    "# Get LSS-Only, M_nu marginalized Errors \n",
    "p1 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p2 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p3 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "\n",
    "boss_forecast.fix_param('tau_reio')\n",
    "desi_forecast.fix_param('tau_reio')\n",
    "euclid_forecast.fix_param('tau_reio')\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    for midx, mval in enumerate(masses): \n",
    "        p1[tidx][midx]=np.sqrt(boss_forecast.final_covs_lss[tidx][midx][param][boss_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "        p2[tidx][midx]=np.sqrt(desi_forecast.final_covs_lss[tidx][midx][param][desi_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "        p3[tidx][midx]=np.sqrt(euclid_forecast.final_covs_lss[tidx][midx][param][euclid_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "p1 = np.array(p1)\n",
    "p2 = np.array(p2)\n",
    "p3 = np.array(p3)\n",
    "\n",
    "boss_lss_only_Mmarg = np.array(p1)\n",
    "desi_lss_only_Mmarg = np.array(p2)\n",
    "euclid_lss_only_Mmarg = np.array(p3)\n",
    "\n",
    "# Get LSS-Only, M_nu fixed Errors \n",
    "p1 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p2 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p3 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "\n",
    "boss_forecast.fix_param('M_ncdm')\n",
    "desi_forecast.fix_param('M_ncdm')\n",
    "euclid_forecast.fix_param('M_ncdm')\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    for midx, mval in enumerate(masses): \n",
    "        p1[tidx][midx]=np.sqrt(boss_forecast.final_covs_lss[tidx][midx][param][boss_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "        p2[tidx][midx]=np.sqrt(desi_forecast.final_covs_lss[tidx][midx][param][desi_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "        p3[tidx][midx]=np.sqrt(euclid_forecast.final_covs_lss[tidx][midx][param][euclid_forecast.final_covs_lss[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "p1 = np.array(p1)\n",
    "p2 = np.array(p2)\n",
    "p3 = np.array(p3)\n",
    "\n",
    "boss_lss_only_Mfixed = np.array(p1)\n",
    "desi_lss_only_Mfixed = np.array(p2)\n",
    "euclid_lss_only_Mfixed = np.array(p3)\n",
    "\n",
    "# Get M_nu marginalized Errors \n",
    "p1 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p2 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p3 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p4 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p5 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p6 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p7 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p8 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p9 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p10 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p11 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p12 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "\n",
    "boss_forecast.clear()\n",
    "desi_forecast.clear()\n",
    "euclid_forecast.clear()\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    for midx, mval in enumerate(masses): \n",
    "        p1[tidx][midx]=np.sqrt(boss_forecast.final_covs_Planck[tidx][midx][param][boss_forecast.final_covs_Planck[tidx][midx].columns.get_loc(param)])\n",
    "        p2[tidx][midx]=np.sqrt(desi_forecast.final_covs_Planck[tidx][midx][param][desi_forecast.final_covs_Planck[tidx][midx].columns.get_loc(param)])\n",
    "        p3[tidx][midx]=np.sqrt(euclid_forecast.final_covs_Planck[tidx][midx][param][euclid_forecast.final_covs_Planck[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "        p4[tidx][midx]=np.sqrt(boss_forecast.final_covs_CMBS4[tidx][midx][param][boss_forecast.final_covs_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "        p5[tidx][midx]=np.sqrt(desi_forecast.final_covs_CMBS4[tidx][midx][param][desi_forecast.final_covs_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "        p6[tidx][midx]=np.sqrt(euclid_forecast.final_covs_CMBS4[tidx][midx][param][euclid_forecast.final_covs_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "        p7[tidx][midx]=np.sqrt(boss_forecast.final_covs_Full_Planck[tidx][midx][param][boss_forecast.final_covs_Full_Planck[tidx][midx].columns.get_loc(param)])\n",
    "        p8[tidx][midx]=np.sqrt(desi_forecast.final_covs_Full_Planck[tidx][midx][param][desi_forecast.final_covs_Full_Planck[tidx][midx].columns.get_loc(param)])\n",
    "        p9[tidx][midx]=np.sqrt(euclid_forecast.final_covs_Full_Planck[tidx][midx][param][euclid_forecast.final_covs_Full_Planck[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "        p10[tidx][midx]=np.sqrt(boss_forecast.final_covs_Full_CMBS4[tidx][midx][param][boss_forecast.final_covs_Full_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "        p11[tidx][midx]=np.sqrt(desi_forecast.final_covs_Full_CMBS4[tidx][midx][param][desi_forecast.final_covs_Full_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "        p12[tidx][midx]=np.sqrt(euclid_forecast.final_covs_Full_CMBS4[tidx][midx][param][euclid_forecast.final_covs_Full_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "p1 = np.array(p1)\n",
    "p2 = np.array(p2)\n",
    "p3 = np.array(p3)\n",
    "p4 = np.array(p4)\n",
    "p5 = np.array(p5)\n",
    "p6 = np.array(p6)\n",
    "p7 = np.array(p7)\n",
    "p8 = np.array(p8)\n",
    "p9 = np.array(p9)\n",
    "p10 = np.array(p10)\n",
    "p11 = np.array(p11)\n",
    "p12 = np.array(p12)\n",
    "\n",
    "boss_Planck_only_Mmarg = np.array(p1)\n",
    "desi_Planck_only_Mmarg = np.array(p2)\n",
    "euclid_Planck_only_Mmarg = np.array(p3)\n",
    "boss_CMBS4_only_Mmarg = np.array(p4)\n",
    "desi_CMBS4_only_Mmarg = np.array(p5)\n",
    "euclid_CMBS4_only_Mmarg = np.array(p6)\n",
    "\n",
    "boss_Full_Planck_Mmarg = np.array(p7)\n",
    "desi_Full_Planck_Mmarg = np.array(p8)\n",
    "euclid_Full_Planck_Mmarg = np.array(p9)\n",
    "boss_Full_CMBS4_Mmarg = np.array(p10)\n",
    "desi_Full_CMBS4_Mmarg = np.array(p11)\n",
    "euclid_Full_CMBS4_Mmarg = np.array(p12)\n",
    "\n",
    "# Get M_nu fixed Errors \n",
    "p1 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p2 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p3 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p4 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p5 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p6 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p7 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p8 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p9 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p10 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p11 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "p12 = [[0. for midx, mval in enumerate(masses)] for tidx, tval in enumerate(temps)]\n",
    "\n",
    "boss_forecast.fix_param('M_ncdm')\n",
    "desi_forecast.fix_param('M_ncdm')\n",
    "euclid_forecast.fix_param('M_ncdm')\n",
    "\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    for midx, mval in enumerate(masses): \n",
    "        p1[tidx][midx]=np.sqrt(boss_forecast.final_covs_Planck[tidx][midx][param][boss_forecast.final_covs_Planck[tidx][midx].columns.get_loc(param)])\n",
    "        p2[tidx][midx]=np.sqrt(desi_forecast.final_covs_Planck[tidx][midx][param][desi_forecast.final_covs_Planck[tidx][midx].columns.get_loc(param)])\n",
    "        p3[tidx][midx]=np.sqrt(euclid_forecast.final_covs_Planck[tidx][midx][param][euclid_forecast.final_covs_Planck[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "        p4[tidx][midx]=np.sqrt(boss_forecast.final_covs_CMBS4[tidx][midx][param][boss_forecast.final_covs_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "        p5[tidx][midx]=np.sqrt(desi_forecast.final_covs_CMBS4[tidx][midx][param][desi_forecast.final_covs_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "        p6[tidx][midx]=np.sqrt(euclid_forecast.final_covs_CMBS4[tidx][midx][param][euclid_forecast.final_covs_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "        p7[tidx][midx]=np.sqrt(boss_forecast.final_covs_Full_Planck[tidx][midx][param][boss_forecast.final_covs_Full_Planck[tidx][midx].columns.get_loc(param)])\n",
    "        p8[tidx][midx]=np.sqrt(desi_forecast.final_covs_Full_Planck[tidx][midx][param][desi_forecast.final_covs_Full_Planck[tidx][midx].columns.get_loc(param)])\n",
    "        p9[tidx][midx]=np.sqrt(euclid_forecast.final_covs_Full_Planck[tidx][midx][param][euclid_forecast.final_covs_Full_Planck[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "        p10[tidx][midx]=np.sqrt(boss_forecast.final_covs_Full_CMBS4[tidx][midx][param][boss_forecast.final_covs_Full_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "        p11[tidx][midx]=np.sqrt(desi_forecast.final_covs_Full_CMBS4[tidx][midx][param][desi_forecast.final_covs_Full_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "        p12[tidx][midx]=np.sqrt(euclid_forecast.final_covs_Full_CMBS4[tidx][midx][param][euclid_forecast.final_covs_Full_CMBS4[tidx][midx].columns.get_loc(param)])\n",
    "\n",
    "p1 = np.array(p1)\n",
    "p2 = np.array(p2)\n",
    "p3 = np.array(p3)\n",
    "p4 = np.array(p4)\n",
    "p5 = np.array(p5)\n",
    "p6 = np.array(p6)\n",
    "p7 = np.array(p7)\n",
    "p8 = np.array(p8)\n",
    "p9 = np.array(p9)\n",
    "p10 = np.array(p10)\n",
    "p11 = np.array(p11)\n",
    "p12 = np.array(p12)\n",
    "\n",
    "boss_Planck_only_Mfixed = np.array(p1)\n",
    "desi_Planck_only_Mfixed = np.array(p2)\n",
    "euclid_Planck_only_Mfixed = np.array(p3)\n",
    "boss_CMBS4_only_Mfixed = np.array(p4)\n",
    "desi_CMBS4_only_Mfixed = np.array(p5)\n",
    "euclid_CMBS4_only_Mfixed = np.array(p6)\n",
    "\n",
    "boss_Full_Planck_Mfixed = np.array(p7)\n",
    "desi_Full_Planck_Mfixed = np.array(p8)\n",
    "euclid_Full_Planck_Mfixed = np.array(p9)\n",
    "boss_Full_CMBS4_Mfixed = np.array(p10)\n",
    "desi_Full_CMBS4_Mfixed = np.array(p11)\n",
    "euclid_Full_CMBS4_Mfixed = np.array(p12)\n",
    "\n",
    "savepath='/Users/nicholasdeporzio/Desktop/06.05.2020/'\n",
    "\n",
    "np.savetxt(savepath+'boss_lss_only_Mnu_marginalized.txt', boss_lss_only_Mmarg)\n",
    "np.savetxt(savepath+'boss_lss_only_Mnu_fixed.txt', boss_lss_only_Mfixed)\n",
    "np.savetxt(savepath+'boss_Planck_only_Mnu_marginalized.txt', boss_Planck_only_Mmarg)\n",
    "np.savetxt(savepath+'boss_Planck_only_Mnu_fixed.txt', boss_Planck_only_Mfixed)\n",
    "np.savetxt(savepath+'boss_CMBS4_only_Mnu_marginalized.txt', boss_CMBS4_only_Mmarg)\n",
    "np.savetxt(savepath+'boss_CMBS4_only_Mnu_fixed.txt', boss_CMBS4_only_Mfixed)\n",
    "np.savetxt(savepath+'boss_Full_Planck_Mnu_marginalized.txt', boss_Full_Planck_Mmarg)\n",
    "np.savetxt(savepath+'boss_Full_Planck_Mnu_fixed.txt', boss_Full_Planck_Mfixed)\n",
    "np.savetxt(savepath+'boss_Full_CMBS4_Mnu_marginalized.txt', boss_Full_CMBS4_Mmarg)\n",
    "np.savetxt(savepath+'boss_Full_CMBS4_Mnu_fixed.txt', boss_Full_CMBS4_Mfixed)\n",
    "\n",
    "np.savetxt(savepath+'desi_lss_only_Mnu_marginalized.txt', desi_lss_only_Mmarg)\n",
    "np.savetxt(savepath+'desi_lss_only_Mnu_fixed.txt', desi_lss_only_Mfixed)\n",
    "np.savetxt(savepath+'desi_Planck_only_Mnu_marginalized.txt', desi_Planck_only_Mmarg)\n",
    "np.savetxt(savepath+'desi_Planck_only_Mnu_fixed.txt', desi_Planck_only_Mfixed)\n",
    "np.savetxt(savepath+'desi_CMBS4_only_Mnu_marginalized.txt', desi_CMBS4_only_Mmarg)\n",
    "np.savetxt(savepath+'desi_CMBS4_only_Mnu_fixed.txt', desi_CMBS4_only_Mfixed)\n",
    "np.savetxt(savepath+'desi_Full_Planck_Mnu_marginalized.txt', desi_Full_Planck_Mmarg)\n",
    "np.savetxt(savepath+'desi_Full_Planck_Mnu_fixed.txt', desi_Full_Planck_Mfixed)\n",
    "np.savetxt(savepath+'desi_Full_CMBS4_Mnu_marginalized.txt', desi_Full_CMBS4_Mmarg)\n",
    "np.savetxt(savepath+'desi_Full_CMBS4_Mnu_fixed.txt', desi_Full_CMBS4_Mfixed)\n",
    "\n",
    "np.savetxt(savepath+'euclid_lss_only_Mnu_marginalized.txt', euclid_lss_only_Mmarg)\n",
    "np.savetxt(savepath+'euclid_lss_only_Mnu_fixed.txt', euclid_lss_only_Mfixed)\n",
    "np.savetxt(savepath+'euclid_Planck_only_Mnu_marginalized.txt', euclid_Planck_only_Mmarg)\n",
    "np.savetxt(savepath+'euclid_Planck_only_Mnu_fixed.txt', euclid_Planck_only_Mfixed)\n",
    "np.savetxt(savepath+'euclid_CMBS4_only_Mnu_marginalized.txt', euclid_CMBS4_only_Mmarg)\n",
    "np.savetxt(savepath+'euclid_CMBS4_only_Mnu_fixed.txt', euclid_CMBS4_only_Mfixed)\n",
    "np.savetxt(savepath+'euclid_Full_Planck_Mnu_marginalized.txt', euclid_Full_Planck_Mmarg)\n",
    "np.savetxt(savepath+'euclid_Full_Planck_Mnu_fixed.txt', euclid_Full_Planck_Mfixed)\n",
    "np.savetxt(savepath+'euclid_Full_CMBS4_Mnu_marginalized.txt', euclid_Full_CMBS4_Mmarg)\n",
    "np.savetxt(savepath+'euclid_Full_CMBS4_Mnu_fixed.txt', euclid_Full_CMBS4_Mfixed)\n",
    "        \n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, boss_lss_only_Mmarg[tidx, :], linestyle='solid', color='blue', label='BOSS, M Marginalized')\n",
    "    plt.semilogx(masses, boss_lss_only_Mfixed[tidx, :], linestyle='dashed', color='blue', label='BOSS, M Fixed')\n",
    "    plt.semilogx(masses, desi_lss_only_Mmarg[tidx, :], linestyle='solid', color='orange', label='DESI, M Marginalized')\n",
    "    plt.semilogx(masses, desi_lss_only_Mfixed[tidx, :], linestyle='dashed', color='orange', label='DESI, M Fixed')\n",
    "    plt.semilogx(masses, euclid_lss_only_Mmarg[tidx, :], linestyle='solid', color='pink', label='Euclid, M Marginalized')\n",
    "    plt.semilogx(masses, euclid_lss_only_Mfixed[tidx, :], linestyle='dashed', color='pink', label='Euclid, M Fixed')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\\sigma(g_X)$\", fontsize=24)\n",
    "    plt.title('LSS Only, T='+str(tval)+'K', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    plt.savefig(savepath+'lss_only_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, boss_Planck_only_Mmarg[tidx, :], linestyle='solid', color='blue', label='Planck, M Marginalized')\n",
    "    plt.semilogx(masses, boss_Planck_only_Mfixed[tidx, :], linestyle='dashed', color='blue', label='Planck, M Fixed')\n",
    "    plt.semilogx(masses, boss_CMBS4_only_Mmarg[tidx, :], linestyle='solid', color='orange', label='CMBS4, M Marginalized')\n",
    "    plt.semilogx(masses, boss_CMBS4_only_Mfixed[tidx, :], linestyle='dashed', color='orange', label='CMBS4, M Fixed')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\\sigma(g_X)$\", fontsize=24)\n",
    "    plt.title('CMB Only, T='+str(tval)+'K', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    plt.savefig(savepath+'cmb_only_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, boss_Full_Planck_Mmarg[tidx, :], linestyle='solid', color='blue', label='BOSS, M Marginalized')\n",
    "    plt.semilogx(masses, boss_Full_Planck_Mfixed[tidx, :], linestyle='dashed', color='blue', label='BOSS, M Fixed')\n",
    "    plt.semilogx(masses, desi_Full_Planck_Mmarg[tidx, :], linestyle='solid', color='orange', label='DESI, M Marginalized')\n",
    "    plt.semilogx(masses, desi_Full_Planck_Mfixed[tidx, :], linestyle='dashed', color='orange', label='DESI, M Fixed')\n",
    "    plt.semilogx(masses, euclid_Full_Planck_Mmarg[tidx, :], linestyle='solid', color='pink', label='Euclid, M Marginalized')\n",
    "    plt.semilogx(masses, euclid_Full_Planck_Mfixed[tidx, :], linestyle='dashed', color='pink', label='Euclid, M Fixed')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\\sigma(g_X)$\", fontsize=24)\n",
    "    plt.title('LSS + Planck, T='+str(tval)+'K', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    plt.savefig(savepath+'Full_Planck_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, boss_Full_CMBS4_Mmarg[tidx, :], linestyle='solid', color='blue', label='BOSS, M Marginalized')\n",
    "    plt.semilogx(masses, boss_Full_CMBS4_Mfixed[tidx, :], linestyle='dashed', color='blue', label='BOSS, M Fixed')\n",
    "    plt.semilogx(masses, desi_Full_CMBS4_Mmarg[tidx, :], linestyle='solid', color='orange', label='DESI, M Marginalized')\n",
    "    plt.semilogx(masses, desi_Full_CMBS4_Mfixed[tidx, :], linestyle='dashed', color='orange', label='DESI, M Fixed')\n",
    "    plt.semilogx(masses, euclid_Full_CMBS4_Mmarg[tidx, :], linestyle='solid', color='pink', label='Euclid, M Marginalized')\n",
    "    plt.semilogx(masses, euclid_Full_CMBS4_Mfixed[tidx, :], linestyle='dashed', color='pink', label='Euclid, M Fixed')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\\sigma(g_X)$\", fontsize=24)\n",
    "    plt.title('LSS + CMBS4, T='+str(tval)+'K', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    plt.savefig(savepath+'Full_CMBS4_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "\n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, boss_lss_only_Mmarg[tidx, :], linestyle='solid', color='blue', label='BOSS Only, M Marginalized')\n",
    "    plt.semilogx(masses, boss_lss_only_Mfixed[tidx, :], linestyle='dashed', color='blue', label='BOSS Only, M Fixed')\n",
    "    plt.semilogx(masses, boss_Full_Planck_Mmarg[tidx, :], linestyle='solid', color='orange', label='BOSS+Planck, M Marginalized')\n",
    "    plt.semilogx(masses, boss_Full_Planck_Mfixed[tidx, :], linestyle='dashed', color='orange', label='BOSS+Planck, M Fixed')\n",
    "    plt.semilogx(masses, boss_Full_CMBS4_Mmarg[tidx, :], linestyle='solid', color='pink', label='BOSS+CMBS4, M Marginalized')\n",
    "    plt.semilogx(masses, boss_Full_CMBS4_Mfixed[tidx, :], linestyle='dashed', color='pink', label='BOSS+CMBS4, M Fixed')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\\sigma(g_X)$\", fontsize=24)\n",
    "    plt.title('LSS + CMBS4, T='+str(tval)+'K', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    plt.savefig(savepath+'boss_all_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, desi_lss_only_Mmarg[tidx, :], linestyle='solid', color='blue', label='DESI Only, M Marginalized')\n",
    "    plt.semilogx(masses, desi_lss_only_Mfixed[tidx, :], linestyle='dashed', color='blue', label='DESI Only, M Fixed')\n",
    "    plt.semilogx(masses, desi_Full_Planck_Mmarg[tidx, :], linestyle='solid', color='orange', label='DESI+Planck, M Marginalized')\n",
    "    plt.semilogx(masses, desi_Full_Planck_Mfixed[tidx, :], linestyle='dashed', color='orange', label='DESI+Planck, M Fixed')\n",
    "    plt.semilogx(masses, desi_Full_CMBS4_Mmarg[tidx, :], linestyle='solid', color='pink', label='DESI+CMBS4, M Marginalized')\n",
    "    plt.semilogx(masses, desi_Full_CMBS4_Mfixed[tidx, :], linestyle='dashed', color='pink', label='DESI+CMBS4, M Fixed')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\\sigma(g_X)$\", fontsize=24)\n",
    "    plt.title('LSS + CMBS4, T='+str(tval)+'K', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    plt.savefig(savepath+'desi_all_'+str(tval)+'K.png')\n",
    "    plt.show()\n",
    "    \n",
    "for tidx, tval in enumerate(temps): \n",
    "    plt.figure(figsize=(15,7.5))\n",
    "    plt.semilogx(masses, euclid_lss_only_Mmarg[tidx, :], linestyle='solid', color='blue', label='Euclid Only, M Marginalized')\n",
    "    plt.semilogx(masses, euclid_lss_only_Mfixed[tidx, :], linestyle='dashed', color='blue', label='Euclid Only, M Fixed')\n",
    "    plt.semilogx(masses, euclid_Full_Planck_Mmarg[tidx, :], linestyle='solid', color='orange', label='Euclid+Planck, M Marginalized')\n",
    "    plt.semilogx(masses, euclid_Full_Planck_Mfixed[tidx, :], linestyle='dashed', color='orange', label='Euclid+Planck, M Fixed')\n",
    "    plt.semilogx(masses, euclid_Full_CMBS4_Mmarg[tidx, :], linestyle='solid', color='pink', label='Euclid+CMBS4, M Marginalized')\n",
    "    plt.semilogx(masses, euclid_Full_CMBS4_Mfixed[tidx, :], linestyle='dashed', color='pink', label='Euclid+CMBS4, M Fixed')\n",
    "    plt.xlabel(r\"Relic Mass [eV]\", fontsize=24)                                       \n",
    "    plt.ylabel(r\"$\\sigma(g_X)$\", fontsize=24)\n",
    "    plt.title('LSS + CMBS4, T='+str(tval)+'K', fontsize=24)\n",
    "    plt.legend(fontsize=18, loc='upper right')                                       \n",
    "    plt.tick_params(axis='x', which='minttor')                                      \n",
    "    plt.grid(True, which='minor')                 \n",
    "    #plt.xlim(1, 5)\n",
    "    #plt.ylim(0, 1.5)\n",
    "    plt.savefig(savepath+'euclid_all_'+str(tval)+'K.png')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-18T14:00:36.070430Z",
     "start_time": "2020-05-18T14:00:36.067015Z"
    },
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Interpolate Error at 3-sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-11T19:03:55.655046Z",
     "start_time": "2020-06-11T19:03:55.420049Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import scipy.interpolate\n",
    "import numpy as np\n",
    "\n",
    "importpath = '/Users/nicholasdeporzio/Desktop/Final_Results_Hold_06.05.2020/N_errors/'\n",
    "\n",
    "boss_lss_only_Mnu_marginalized = np.loadtxt(importpath+'boss_lss_only_Mnu_marginalized.txt')\n",
    "boss_lss_only_Mnu_fixed = np.loadtxt(importpath+'boss_lss_only_Mnu_fixed.txt')\n",
    "boss_Planck_only_Mnu_marginalized = np.loadtxt(importpath+'boss_Planck_only_Mnu_marginalized.txt')\n",
    "boss_Planck_only_Mnu_fixed = np.loadtxt(importpath+'boss_Planck_only_Mnu_fixed.txt')\n",
    "boss_CMBS4_only_Mnu_marginalized = np.loadtxt(importpath+'boss_CMBS4_only_Mnu_marginalized.txt')\n",
    "boss_CMBS4_only_Mnu_fixed = np.loadtxt(importpath+'boss_CMBS4_only_Mnu_fixed.txt')\n",
    "boss_Full_Planck_Mnu_marginalized = np.loadtxt(importpath+'boss_Full_Planck_Mnu_marginalized.txt')\n",
    "boss_Full_Planck_Mnu_fixed = np.loadtxt(importpath+'boss_Full_Planck_Mnu_fixed.txt')\n",
    "boss_Full_CMBS4_Mnu_marginalized = np.loadtxt(importpath+'boss_Full_CMBS4_Mnu_marginalized.txt')\n",
    "boss_Full_CMBS4_Mnu_fixed = np.loadtxt(importpath+'boss_Full_CMBS4_Mnu_fixed.txt')\n",
    "\n",
    "desi_lss_only_Mnu_marginalized = np.loadtxt(importpath+'desi_lss_only_Mnu_marginalized.txt')\n",
    "desi_lss_only_Mnu_fixed = np.loadtxt(importpath+'desi_lss_only_Mnu_fixed.txt')\n",
    "desi_Planck_only_Mnu_marginalized = np.loadtxt(importpath+'desi_Planck_only_Mnu_marginalized.txt')\n",
    "desi_Planck_only_Mnu_fixed = np.loadtxt(importpath+'desi_Planck_only_Mnu_fixed.txt')\n",
    "desi_CMBS4_only_Mnu_marginalized = np.loadtxt(importpath+'desi_CMBS4_only_Mnu_marginalized.txt')\n",
    "desi_CMBS4_only_Mnu_fixed = np.loadtxt(importpath+'desi_CMBS4_only_Mnu_fixed.txt')\n",
    "desi_Full_Planck_Mnu_marginalized = np.loadtxt(importpath+'desi_Full_Planck_Mnu_marginalized.txt')\n",
    "desi_Full_Planck_Mnu_fixed = np.loadtxt(importpath+'desi_Full_Planck_Mnu_fixed.txt')\n",
    "desi_Full_CMBS4_Mnu_marginalized = np.loadtxt(importpath+'desi_Full_CMBS4_Mnu_marginalized.txt')\n",
    "desi_Full_CMBS4_Mnu_fixed = np.loadtxt(importpath+'desi_Full_CMBS4_Mnu_fixed.txt')\n",
    "\n",
    "euclid_lss_only_Mnu_marginalized = np.loadtxt(importpath+'euclid_lss_only_Mnu_marginalized.txt')\n",
    "euclid_lss_only_Mnu_fixed = np.loadtxt(importpath+'euclid_lss_only_Mnu_fixed.txt')\n",
    "euclid_Planck_only_Mnu_marginalized = np.loadtxt(importpath+'euclid_Planck_only_Mnu_marginalized.txt')\n",
    "euclid_Planck_only_Mnu_fixed = np.loadtxt(importpath+'euclid_Planck_only_Mnu_fixed.txt')\n",
    "euclid_CMBS4_only_Mnu_marginalized = np.loadtxt(importpath+'euclid_CMBS4_only_Mnu_marginalized.txt')\n",
    "euclid_CMBS4_only_Mnu_fixed = np.loadtxt(importpath+'euclid_CMBS4_only_Mnu_fixed.txt')\n",
    "euclid_Full_Planck_Mnu_marginalized = np.loadtxt(importpath+'euclid_Full_Planck_Mnu_marginalized.txt')\n",
    "euclid_Full_Planck_Mnu_fixed = np.loadtxt(importpath+'euclid_Full_Planck_Mnu_fixed.txt')\n",
    "euclid_Full_CMBS4_Mnu_marginalized = np.loadtxt(importpath+'euclid_Full_CMBS4_Mnu_marginalized.txt')\n",
    "euclid_Full_CMBS4_Mnu_fixed = np.loadtxt(importpath+'euclid_Full_CMBS4_Mnu_fixed.txt')\n",
    "\n",
    "temp_index=0\n",
    "masses = np.geomspace(0.01, 10., 21)\n",
    "temps=[0.79, 0.91, 0.94, 1.08] \n",
    "surveys=['boss', 'desi', 'euclid']\n",
    "\n",
    "if temp_index==0: \n",
    "    gamma=1\n",
    "    gX=1\n",
    "elif temp_index==1:\n",
    "    gamma=0\n",
    "    gX=2\n",
    "elif temp_index==2: \n",
    "    gamma=1\n",
    "    gX=2\n",
    "elif temp_index==3: \n",
    "    gamma=0\n",
    "    gX=4\n",
    "\n",
    "#sigma_threshold = gX/3.\n",
    "sigma_threshold = 1./3.\n",
    "\n",
    "def m_convert(m_eq): \n",
    "    return m_eq * np.power(2./gX, 1./4.) * np.power(7./8., gamma/4.) * np.power(6./7., gamma)\n",
    "\n",
    "print('Temperature: ', temps[temp_index])\n",
    "\n",
    "int_boss_lss_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_lss_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_boss_lss_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_lss_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_boss_Planck_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_Planck_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_boss_Planck_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_Planck_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_boss_CMBS4_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_CMBS4_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_boss_CMBS4_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_CMBS4_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_boss_Full_Planck_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_Full_Planck_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_boss_Full_Planck_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_Full_Planck_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_boss_Full_CMBS4_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_Full_CMBS4_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_boss_Full_CMBS4_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, boss_Full_CMBS4_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "\n",
    "int_desi_lss_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_lss_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_desi_lss_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_lss_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_desi_Planck_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_Planck_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_desi_Planck_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_Planck_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_desi_CMBS4_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_CMBS4_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_desi_CMBS4_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_CMBS4_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_desi_Full_Planck_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_Full_Planck_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_desi_Full_Planck_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_Full_Planck_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_desi_Full_CMBS4_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_Full_CMBS4_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_desi_Full_CMBS4_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, desi_Full_CMBS4_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "\n",
    "int_euclid_lss_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_lss_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_euclid_lss_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_lss_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_euclid_Planck_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_Planck_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_euclid_Planck_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_Planck_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_euclid_CMBS4_only_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_CMBS4_only_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_euclid_CMBS4_only_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_CMBS4_only_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_euclid_Full_Planck_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_Full_Planck_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_euclid_Full_Planck_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_Full_Planck_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "int_euclid_Full_CMBS4_Mnu_marginalized = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_Full_CMBS4_Mnu_marginalized[temp_index, :]-sigma_threshold)\n",
    "int_euclid_Full_CMBS4_Mnu_fixed = scipy.interpolate.InterpolatedUnivariateSpline(masses, euclid_Full_CMBS4_Mnu_fixed[temp_index, :]-sigma_threshold)\n",
    "\n",
    "print('**************************')\n",
    "print('BOSS Only (M_nu Marginalized): ', m_convert(int_boss_lss_only_Mnu_marginalized.roots()))\n",
    "print('BOSS Only (M_nu Fixed): ', m_convert(int_boss_lss_only_Mnu_fixed.roots()))\n",
    "print('Planck Only (M_nu Marginalized): ', m_convert(int_boss_Planck_only_Mnu_marginalized.roots()))\n",
    "print('Planck Only (M_nu Fixed): ', m_convert(int_boss_Planck_only_Mnu_fixed.roots()))\n",
    "print('CMBS4 Only (M_nu Marginalized): ', m_convert(int_boss_CMBS4_only_Mnu_marginalized.roots()))\n",
    "print('CMBS4 Only (M_nu Fixed): ', m_convert(int_boss_CMBS4_only_Mnu_fixed.roots()))\n",
    "print('BOSS + Planck (M_nu Marginalized): ', m_convert(int_boss_Full_Planck_Mnu_marginalized.roots()))\n",
    "print('BOSS + Planck (M_nu Fixed): ', m_convert(int_boss_Full_Planck_Mnu_fixed.roots()))\n",
    "print('BOSS + CMBS4 (M_nu Marginalized): ', m_convert(int_boss_Full_CMBS4_Mnu_marginalized.roots()))\n",
    "print('BOSS + CMBS4 (M_nu Fixed): ', m_convert(int_boss_Full_CMBS4_Mnu_fixed.roots()))\n",
    "print('**************************')\n",
    "print('DESI Only (M_nu Marginalized): ', m_convert(int_desi_lss_only_Mnu_marginalized.roots()))\n",
    "print('DESI Only (M_nu Fixed): ', m_convert(int_desi_lss_only_Mnu_fixed.roots()))\n",
    "print('Planck Only (M_nu Marginalized): ', m_convert(int_desi_Planck_only_Mnu_marginalized.roots()))\n",
    "print('Planck Only (M_nu Fixed): ', m_convert(int_desi_Planck_only_Mnu_fixed.roots()))\n",
    "print('CMBS4 Only (M_nu Marginalized): ', m_convert(int_desi_CMBS4_only_Mnu_marginalized.roots()))\n",
    "print('CMBS4 Only (M_nu Fixed): ', m_convert(int_desi_CMBS4_only_Mnu_fixed.roots()))\n",
    "print('DESI + Planck (M_nu Marginalized): ', m_convert(int_desi_Full_Planck_Mnu_marginalized.roots()))\n",
    "print('DESI + Planck (M_nu Fixed): ', m_convert(int_desi_Full_Planck_Mnu_fixed.roots()))\n",
    "print('DESI + CMBS4 (M_nu Marginalized): ', m_convert(int_desi_Full_CMBS4_Mnu_marginalized.roots()))\n",
    "print('DESI + CMBS4 (M_nu Fixed): ', m_convert(int_desi_Full_CMBS4_Mnu_fixed.roots()))\n",
    "print('**************************')\n",
    "print('Euclid Only (M_nu Marginalized): ', m_convert(int_euclid_lss_only_Mnu_marginalized.roots()))\n",
    "print('Euclid Only (M_nu Fixed): ', m_convert(int_euclid_lss_only_Mnu_fixed.roots()))\n",
    "print('Planck Only (M_nu Marginalized): ', m_convert(int_euclid_Planck_only_Mnu_marginalized.roots()))\n",
    "print('Planck Only (M_nu Fixed): ', m_convert(int_euclid_Planck_only_Mnu_fixed.roots()))\n",
    "print('CMBS4 Only (M_nu Marginalized): ', m_convert(int_euclid_CMBS4_only_Mnu_marginalized.roots()))\n",
    "print('CMBS4 Only (M_nu Fixed): ', m_convert(int_euclid_CMBS4_only_Mnu_fixed.roots()))\n",
    "print('Euclid + Planck (M_nu Marginalized): ', m_convert(int_euclid_Full_Planck_Mnu_marginalized.roots()))\n",
    "print('Euclid + Planck (M_nu Fixed): ', m_convert(int_euclid_Full_Planck_Mnu_fixed.roots()))\n",
    "print('Euclid + CMBS4 (M_nu Marginalized): ', m_convert(int_euclid_Full_CMBS4_Mnu_marginalized.roots()))\n",
    "print('Euclid + CMBS4 (M_nu Fixed): ', m_convert(int_euclid_Full_CMBS4_Mnu_fixed.roots()))\n",
    "print('**************************')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-26T03:36:29.680025Z",
     "start_time": "2020-05-26T03:36:29.676901Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Equivalent N_eff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-11T19:05:58.800744Z",
     "start_time": "2020-06-11T19:05:58.797448Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def sig_Neff(g_x, sig_g, T_x): \n",
    "    return g_x * (sig_g/2.) * np.power(T_x/1.95, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-11T19:05:59.123853Z",
     "start_time": "2020-06-11T19:05:59.107907Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "for tidx, tval in enumerate(temps): \n",
    "    if tidx==0: \n",
    "        gx = 1.\n",
    "    elif tidx==1: \n",
    "        gx = 2\n",
    "    elif tidx==2: \n",
    "        gx = 2\n",
    "    elif tidx==3: \n",
    "        gx = 3\n",
    "        \n",
    "    print('Error on N_eff for T=', tval, '...')\n",
    "    print('BOSS Only: ', sig_Neff(gx, boss_lss_only_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, boss_lss_only_Mnu_fixed[tidx][0], tval), ')')\n",
    "    print('DESI Only: ', sig_Neff(gx, desi_lss_only_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, desi_lss_only_Mnu_fixed[tidx][0], tval), ')')\n",
    "    print('Euclid Only: ', sig_Neff(gx, euclid_lss_only_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, euclid_lss_only_Mnu_fixed[tidx][0], tval), ')')\n",
    "\n",
    "    print('Planck Only: ', sig_Neff(gx, desi_Planck_only_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, boss_Planck_only_Mnu_fixed[tidx][0], tval), ')')\n",
    "    print('CMBS4 Only: ', sig_Neff(gx, desi_CMBS4_only_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, boss_CMBS4_only_Mnu_fixed[tidx][0], tval), ')')\n",
    "\n",
    "    print('BOSS + Planck: ', sig_Neff(gx, boss_Full_Planck_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, boss_Full_Planck_Mnu_fixed[tidx][0], tval), ')')\n",
    "    print('DESI + Planck: ', sig_Neff(gx, desi_Full_Planck_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, desi_Full_Planck_Mnu_fixed[tidx][0], tval), ')')\n",
    "    print('Euclid + Planck: ', sig_Neff(gx, euclid_Full_Planck_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, euclid_Full_Planck_Mnu_fixed[tidx][0], tval), ')')\n",
    "\n",
    "    print('BOSS + CMBS4: ', sig_Neff(gx, boss_Full_CMBS4_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, boss_Full_CMBS4_Mnu_fixed[tidx][0], tval), ')')\n",
    "    print('DESI + CMBS4: ', sig_Neff(gx, desi_Full_CMBS4_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, desi_Full_CMBS4_Mnu_fixed[tidx][0], tval), ')')\n",
    "    print('Euclid + CMBS4: ', sig_Neff(gx, euclid_Full_CMBS4_Mnu_marginalized[tidx][0], tval), ' (', sig_Neff(gx, euclid_Full_CMBS4_Mnu_fixed[tidx][0], tval), ')')\n",
    "    print('********************')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-08T07:20:29.861173Z",
     "start_time": "2020-06-08T07:20:29.848008Z"
    },
    "hidden": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for tidx, tval in enumerate(temps): \n",
    "    if tidx==0: \n",
    "        gx = 1.\n",
    "    elif tidx==1: \n",
    "        gx = 2\n",
    "    elif tidx==2: \n",
    "        gx = 2\n",
    "    elif tidx==3: \n",
    "        gx = 3\n",
    "        \n",
    "    print('Marginalization Shift in Error for T=', tval, '...')\n",
    "    print('BOSS Only: ', max((boss_lss_only_Mnu_marginalized[tidx]/boss_lss_only_Mnu_fixed[tidx])))\n",
    "    print('DESI Only: ', max(desi_lss_only_Mnu_marginalized[tidx]/desi_lss_only_Mnu_fixed[tidx]))\n",
    "    print('Euclid Only: ', max(euclid_lss_only_Mnu_marginalized[tidx]/euclid_lss_only_Mnu_fixed[tidx]))\n",
    "\n",
    "    print('Planck Only: ', max(desi_Planck_only_Mnu_marginalized[tidx]/desi_Planck_only_Mnu_fixed[tidx]))\n",
    "    print('CMBS4 Only: ', max(desi_CMBS4_only_Mnu_marginalized[tidx]/desi_CMBS4_only_Mnu_fixed[tidx]))\n",
    "\n",
    "    print('BOSS + Planck: ', max(boss_Full_Planck_Mnu_marginalized[tidx]/boss_Full_Planck_Mnu_fixed[tidx]))\n",
    "    print('DESI + Planck: ', max(desi_Full_Planck_Mnu_marginalized[tidx]/desi_Full_Planck_Mnu_fixed[tidx]))\n",
    "    print('Euclid + Planck: ', max(euclid_Full_Planck_Mnu_marginalized[tidx]/euclid_Full_Planck_Mnu_fixed[tidx]))\n",
    "\n",
    "    print('BOSS + CMBS4: ', max(boss_Full_CMBS4_Mnu_marginalized[tidx]/boss_Full_CMBS4_Mnu_fixed[tidx]))\n",
    "    print('DESI + CMBS4: ', max(desi_Full_CMBS4_Mnu_marginalized[tidx]/desi_Full_CMBS4_Mnu_fixed[tidx]))\n",
    "    print('Euclid + CMBS4: ', max(euclid_Full_CMBS4_Mnu_marginalized[tidx]/euclid_Full_CMBS4_Mnu_fixed[tidx]))\n",
    "    print('********************')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-08T07:20:11.352133Z",
     "start_time": "2020-06-08T07:20:11.348386Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T18:56:44.757888Z",
     "start_time": "2020-05-12T18:56:44.752385Z"
    },
    "heading_collapsed": true
   },
   "source": [
    "# Debugging: Plot Contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T19:29:25.197691Z",
     "start_time": "2020-05-12T19:29:04.891489Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Show plots inline, and load main getdist plot module and samples class\n",
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "import sys, os\n",
    "sys.path.insert(0,os.path.realpath(os.path.join(os.getcwd(),'..')))\n",
    "from getdist import plots, MCSamples, loadMCSamples\n",
    "from getdist.gaussian_mixtures import GaussianND\n",
    "import getdist\n",
    "# use this *after* importing getdist if you want to use interactive plots\n",
    "# %matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "import IPython\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "print('GetDist Version: %s, Matplotlib version: %s'%(getdist.__version__, plt.matplotlib.__version__))\n",
    "# matplotlib 2 may not work very well without usetex on, can uncomment\n",
    "# plt.rcParams['text.usetex']=True\n",
    "\n",
    "#dataset = 'Full'\n",
    "#dataset = 'LSS'\n",
    "\n",
    "#survey = 'CMBS4'\n",
    "survey = 'Planck'\n",
    "\n",
    "dtype='cmb'\n",
    "mass_index=13\n",
    "temp_index=0\n",
    "\n",
    "masses = np.geomspace(0.01, 10., 21)\n",
    "temps = np.array([0.79, 0.91, 0.94, 1.08])\n",
    "mass=masses[mass_index]\n",
    "temp=temps[temp_index]\n",
    "path='/Volumes/SSD01/FINAL_RESULTS/Fine_Plot/Fisher_Matrices/T'+str(temp)+'K/M'+f'{mass:.3f}'+'/'+survey+'/inv_'+dtype+'fisher.mat'\n",
    "\n",
    "cov=np.loadtxt(path, skiprows=1)\n",
    "\n",
    "fid = [\n",
    "    2.2321e-9, #A_s\n",
    "    0.701, #h\n",
    "    0.967, #n_s\n",
    "    0.02226, #omega_b\n",
    "    0.1127, #omega_cdm\n",
    "    0.0598, #tau_reio\n",
    "    0.06, #M_ncdm\n",
    "    1.0, #N_ncdm\n",
    "#   250., #sigma_fog\n",
    "#    1.0, #beta0\n",
    "#    1.0, #alphak2 \n",
    "]\n",
    "\n",
    "params = [\n",
    "    r'$A_s$', \n",
    "    r'$h$', \n",
    "    r'$n_s$', \n",
    "    r'$\\omega_b$', \n",
    "    r'$\\omega_{cdm}$', \n",
    "    r'$\\tau_{reio}$', \n",
    "    r'$M_\\nu$',\n",
    "    r'$g_\\chi$',  \n",
    "#    r'$\\sigma_{fog}$', \n",
    "#    r'$\\beta_0$', \n",
    "#    r'$\\alpha_{k2}$'\n",
    "]\n",
    "\n",
    "gauss = GaussianND(fid, cov, names=params)\n",
    "\n",
    "h = plots.get_subplot_plotter()\n",
    "\n",
    "# Set a default plot style \n",
    "from getdist.styles.tab10 import style_name\n",
    "plots.set_active_style(style_name)\n",
    "\n",
    "# Adjust plot sizes, styles, etc using getdist.plots.GetDistPlotSettings class\n",
    "h.settings.figure_legend_frame = False #Box around legend\n",
    "h.settings.scaling=False #Defaults to font size for Tick Labels: 11, Labels: 14, Legend: 12, Text: 12\n",
    "h.settings.axes_fontsize=14\n",
    "h.settings.axes_labelsize=20\n",
    "h.settings.axis_tick_x_rotation=45\n",
    "h.settings.legend_fontsize=20\n",
    "#h.settings.axis_tick_y_rotation=90\n",
    "#h.settings.colorbar_tick_rotation=90\n",
    "\n",
    "\n",
    "h.triangle_plot(\n",
    "    [\n",
    "        gauss\n",
    "    ], \n",
    "    filled=True,\n",
    "    legend_labels=[\n",
    "        \"Planck, T=0.79, M=\"+f'{mass:.3f}'\n",
    "    ], \n",
    "    line_args=[\n",
    "        {'ls':'solid', 'lw':2, 'color':'#A51C30'},\n",
    "    ], \n",
    "    contour_args=[\n",
    "        {'ls':'dotted', 'lw':2, 'color':'#A51C30', 'filled' : True},\n",
    "    ])\n",
    "plt.savefig(\"/Users/nicholasdeporzio/Desktop/test/\"+dtype+\"_M\"+f'{mass:.3f}'+\".png\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T19:15:27.448667Z",
     "start_time": "2020-05-12T19:15:27.444295Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.geomspace(0.01, 10., 21)[6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T18:54:34.097037Z",
     "start_time": "2020-05-08T18:54:34.091340Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-07T18:38:19.546453Z",
     "start_time": "2020-05-07T18:38:19.542657Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T13:59:33.008187Z",
     "start_time": "2020-05-05T13:59:30.866526Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T22:08:21.905096Z",
     "start_time": "2020-05-05T22:08:19.001117Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T21:55:26.170496Z",
     "start_time": "2020-05-05T21:55:24.273316Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T21:56:28.997468Z",
     "start_time": "2020-05-05T21:56:28.992393Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T21:59:10.280973Z",
     "start_time": "2020-05-05T21:59:10.275533Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3.8.0",
   "language": "python",
   "name": "3.8.0"
  },
  "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.0"
  },
  "notify_time": "10"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}