GC-elFish-EFT-zeq0p5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EFTofLSS Fisher forecasts on GC multipoles at z=0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Name this run - this will determine the name for the file into which the covariance matrix will be saved\n",
    "runname = 'el024_mcmcfids_EFT'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the science packages/libraries\n",
    "from __future__ import division\n",
    "from __future__ import print_function\n",
    "import numpy as np\n",
    "import scipy, scipy.interpolate, scipy.misc, scipy.integrate, scipy.stats\n",
    "import elgc\n",
    "import getdist, getdist.plots, IPython\n",
    "from matplotlib.patches import Patch\n",
    "from matplotlib.lines import Line2D\n",
    "\n",
    "import os.path as op\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set plotting and printing\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "font = {'size'   : 16, 'family':'STIXGeneral'}\n",
    "axislabelfontsize='x-large'\n",
    "plt.rc('font', **font)\n",
    "plt.rcParams['text.usetex'] = True\n",
    "\n",
    "np.set_printoptions(precision=3,linewidth=200,suppress=False)\n",
    "\n",
    "from matplotlib.patches import Ellipse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['f', 'b', 'b2', 'N', 'cs0', 'cs2', 'cs4']\n"
     ]
    }
   ],
   "source": [
    "# Main question for the Fisher matrix calculation: which parameters are we varying?\n",
    "# Note that this list determines the order in the output FM as well\n",
    "\n",
    "#NB: the cs0, cs2, cs4 params here are actually the cs squared (like in the MCMC)\n",
    "varypars = ['f','b','b2','N','cs0','cs2','cs4']\n",
    "\n",
    "print(varypars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the fiducial model for the Fisher matrix calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fiducial (GR-COLA) cosmology \n",
    "#from cos.rtf file (changed some params names here)\n",
    "\n",
    "zbin = 0.5\n",
    "\n",
    "c=3e5\n",
    "hubble = 0.697\n",
    "H00 = 100*hubble\n",
    "Tcmb = 2.726\n",
    "ns = 0.971\n",
    "om0 = 0.281\n",
    "omb = 0.046\n",
    "sigma80_fid = 0.844 #this is sigma8 at z=0\n",
    "\n",
    "#Mean MCMC values for EFT at z=0.5 \n",
    "fg_fid = 0.7042\n",
    "b_fid = 1.468\n",
    "b2_fid = -0.2756\n",
    "N_fid = 816.1 \n",
    "#NB: the cs0, cs2, cs4 params here are actually the cs squared (like in the MCMC)\n",
    "cs0_fid = 7.650\n",
    "cs2_fid = 2.978\n",
    "cs4_fid = 24.40\n",
    "\n",
    "#kmax at z=0.5 for multipoles: 0,2,4 \n",
    "kmax = [0.245, 0.245,0.129]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ps column names: ('k', 'pl', 'pdd', 'pdt', 'ptt', 'sv', 'plnw', 'pddres', 'pdtres', 'pttres', 'ANW')\n",
      "abc column names: ('k', 'bu2f2nw', 'bu2f3nw', 'bu4f2nw', 'bu4f3nw', 'bu6f3nw', 'b2u2f1nw', 'b2u2f2nw', 'b2u2f3nw', 'b2u2f4nw', 'b2u4f2nw', 'b2u4f3nw', 'b2u4f4nw', 'b2u6f3nw', 'b2u6f4nw', 'b2u8f4nw', 'u2f4nw', 'u4f3nw', 'u4f4nw', 'u6f3nw', 'u6f4nw', 'u8f4nw')\n",
      "bias column names: ('k', 'ddb0b2nw', 'ddb02nw', 'ddb0nw', 'ddb22nw', 'ddb2nw', 'ddnonenw', 'dtb2nw', 'dtb0nw', 'dtnonenw')\n"
     ]
    }
   ],
   "source": [
    "#k values are the same in all files ### Paths to input files\n",
    "PATHTOINPUTFILES = op.abspath('./inputs/')\n",
    "\n",
    "# Function to read the output files from MG-COPTER\n",
    "def read_copter(path,suffix=''):\n",
    "    \n",
    "    #naming files (following the Mathematica notebook)\n",
    "    psfile = ('eft_ps_components_nw'+suffix+'.dat')\n",
    "    abcfile = ('eft_abc_components'+suffix+'.dat')\n",
    "    biasfile = ('eft_bias_components'+suffix+'.dat')\n",
    "    \n",
    "    # set the names of the columns (following the Mathematica notebook)\n",
    "    pscolnames = 'k', 'pl', 'pdd', 'pdt', 'ptt', 'sv', 'plnw', 'pddres', 'pdtres', 'pttres','ANW'\n",
    "    abccolnames = 'k', 'bu2f2nw', 'bu2f3nw', 'bu4f2nw', 'bu4f3nw', 'bu6f3nw', 'b2u2f1nw', 'b2u2f2nw', 'b2u2f3nw',\\\n",
    "    'b2u2f4nw', 'b2u4f2nw', 'b2u4f3nw', 'b2u4f4nw', 'b2u6f3nw', 'b2u6f4nw', 'b2u8f4nw', 'u2f4nw', 'u4f3nw', 'u4f4nw',\\\n",
    "    'u6f3nw', 'u6f4nw', 'u8f4nw' \n",
    "    biascolnames = 'k', 'ddb0b2nw', 'ddb02nw', 'ddb0nw', 'ddb22nw', 'ddb2nw', 'ddnonenw', 'dtb2nw', 'dtb0nw', 'dtnonenw'\n",
    "    \n",
    "    return np.genfromtxt(op.join(path,psfile), names=pscolnames),\\\n",
    "           np.genfromtxt(op.join(path,abcfile), names=abccolnames),\\\n",
    "           np.genfromtxt(op.join(path,biasfile), names=biascolnames)\n",
    "\n",
    "# use numpy to read the whole files, the options will make the\n",
    "# resulting object act a bit like a dictionary (a structured array)\n",
    "# so we can query columns by name\n",
    "ps_fid, abc_fid, bias_fid = read_copter(op.join(PATHTOINPUTFILES,'data_fid'),'_fid')\n",
    "print ('ps column names:', ps_fid.dtype.names )# let's print the names to check\n",
    "print ('abc column names:', abc_fid.dtype.names) # let's print the names to check\n",
    "print ('bias column names:', bias_fid.dtype.names) # let's print the names to check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sv_fid = ps_fid['sv'][0] \n",
    "ANW = ps_fid['ANW'][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the survey properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Set up survey with volume 4 Gpc^3/h^3\n",
    "\n",
    "Area = 15000.0*1.06 #deg^2 \n",
    "Dzbin = 0.2\n",
    "z_vals = np.array([zbin]) \n",
    "Nzbin = len(z_vals)\n",
    "\n",
    "#set kmin \n",
    "kmin = 0.006"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "401\n"
     ]
    }
   ],
   "source": [
    "#construct array of mu values between -1 and 1 (inclusive)\n",
    "mu_step = 0.005\n",
    "mu = mu_step*np.arange((2/mu_step)+1)-1.\n",
    "#print ('mu array', mu)\n",
    "print (len(mu)) # Note that this has to be large enough to convert to multipoles!\n",
    "\n",
    "# Select multipoles\n",
    "MLPS = [0,2,4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define cosmology functions: background and perturbations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.733197112933542 0.7784487714884997\n"
     ]
    }
   ],
   "source": [
    "# Functions for DE later\n",
    "def w_integrand(z,w0=-1.,wa=0.):\n",
    "    return (1 + w0+wa*z/(z+1)) / (1+z)\n",
    "def DE_evo(zc,w0=-1.,wa=0.): \n",
    "    return np.exp(3*scipy.integrate.romberg(lambda z: w_integrand(z,w0,wa), 0, zc))\n",
    "\n",
    "# Define E(z) = H(z)/H0\n",
    "def Ez(zc,w0=-1.,wa=0.):\n",
    "    return np.sqrt((1-om0)*DE_evo(zc,w0,wa) + om0*pow(1+zc,3))\n",
    "\n",
    "# LCDM growth rate and growth factor\n",
    "def fg(zz,w0=-1.,wa=0.,gamma=0.55):\n",
    "    omz=om0*pow(1+zz,3)/pow(Ez(zz,w0,wa),2)\n",
    "    return pow(omz,gamma)\n",
    "\n",
    "def Dg_dz(zz,w0=-1.,wa=0.,gamma=0.55):\n",
    "    return -fg(zz,w0,wa)/(1+zz)\n",
    "\n",
    "def Dgz(zz,w0=-1.,wa=0.,gamma=0.55):\n",
    "    start_z = 0.0\n",
    "    ans = scipy.integrate.romberg(lambda z: Dg_dz(z,w0,wa), start_z, zz)\n",
    "    return np.exp(ans)\n",
    "\n",
    "nz = 100\n",
    "fg_fid_interp = np.zeros(nz)\n",
    "Dg_fid_interp = np.zeros(nz)\n",
    "z_for_interp = np.linspace(0,2,nz)\n",
    "for i,z in enumerate(z_for_interp):\n",
    "    fg_fid_interp[i] = fg(z)\n",
    "    Dg_fid_interp[i] = Dgz(z)\n",
    "\n",
    "fg_fid_i = scipy.interpolate.interp1d(z_for_interp,fg_fid_interp,kind='cubic')\n",
    "Dg_fid_i = scipy.interpolate.interp1d(z_for_interp,Dg_fid_interp,kind='cubic')\n",
    "del z_for_interp, Dg_fid_interp, fg_fid_interp, i, z\n",
    "print (fg_fid_i(zbin), Dg_fid_i(zbin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.994\n"
     ]
    }
   ],
   "source": [
    "#Define E(z)=H(z)/H0\n",
    "def Ez(zc):\n",
    "    return np.sqrt(1-om0+om0*pow(1+zc,3))\n",
    "\n",
    "#Define the comoving distance\n",
    "def drdz(zp):\n",
    "    return (c/H00)/Ez(zp)\n",
    "def rcom(zc):\n",
    "    return scipy.integrate.romberg(drdz,0,zc)\n",
    "def DA(zc):\n",
    "    return rcom(zc)/(1+zc)\n",
    "\n",
    "def Vsur(zc): #  (the hubble**3 factor is for units consistency)\n",
    "    z_min = zc-Dzbin/2\n",
    "    DA_min = DA(z_min)\n",
    "    z_max = zc+Dzbin/2\n",
    "    DA_max = DA(z_max)\n",
    "    \n",
    "    F = (4.*np.pi/3.) / (4.*np.pi*(180./np.pi)**2)   \n",
    "    return F * (DA_max**3 *(1+z_max)**3 - DA_min**3 *(1+z_min)**3) * hubble**3 * Area\n",
    "\n",
    "print(round(Vsur(zbin)/1e9,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#construct the P_lin in 2d (i.e. in k-mu space) needed for the linear covariance used in the MCMC\n",
    "def P_lin(mu_array,f,sig8_z=sigma80_fid*Dg_fid_i(zbin),ps=ps_fid):\n",
    "    muones = mu_array**0.\n",
    "    k_array = ps['k']\n",
    "    \n",
    "    D = sig8_z/sigma80_fid #every term that contains D, contains the same power of sigma80_fid\n",
    "    \n",
    "    blin = 1.46 \n",
    "    \n",
    "    res = blin**2*D**2*np.outer(ps['pl'],muones) + (f/blin)**2*blin**2*D**2*np.outer(ps['pl'],mu_array**4)\\\n",
    "    + 2.*(f/blin)*blin**2*D**2*np.outer(ps['pl'],mu_array**2)\n",
    "                \n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "900 300\n",
      "k are preserved: True\n"
     ]
    }
   ],
   "source": [
    "#construct the P_EFT in 2d (i.e. in k-mu space)\n",
    "def Peft_2d(mu_array,f,\n",
    "         b=b_fid,b2=b2_fid,N=N_fid,\n",
    "         cs0=cs0_fid,cs2=cs2_fid,cs4=cs4_fid,sv=sv_fid,\n",
    "         sig8_z=sigma80_fid*Dg_fid_i(zbin),\n",
    "         ps=ps_fid,abc=abc_fid,bias=bias_fid):\n",
    "    muones = mu_array**0.\n",
    "    k_array = ps['k']\n",
    "    \n",
    "    D = sig8_z/sigma80_fid #every term that contains D, contains the same power of sigma80_fid\n",
    "    \n",
    "    #resummation part\n",
    "    pddresum = D**2 * ps['plnw'] + D**4 * ps['pddres']\\\n",
    "    + np.exp(-0.5*ps['k']**2*D**2*ANW) * ( D**2 * ps['pl'] + D**4 * ps['pdd'] - (D**2*ps['plnw'] + D**4*ps['pddres'])\\\n",
    "                                         +0.5*ps['k']**2*D**2*ANW*(D**2*(ps['pl']-ps['plnw'])))\n",
    "   \n",
    "    pdtresum = D**2 * f * ps['plnw'] + D**4 * f * ps['pdtres']\\\n",
    "    + np.exp(-0.5*ps['k']**2*D**2*ANW) * ( D**2 * f * ps['pl'] + D**4 * f * ps['pdt'] - (D**2*f*ps['plnw'] + D**4*f*ps['pdtres'])\\\n",
    "                                         +0.5*ps['k']**2*D**2*ANW*(f*D**2*(ps['pl']-ps['plnw'])))\n",
    "    \n",
    "    pttresum = D**2 * f**2 * ps['plnw'] + D**4 * f**2 * ps['pttres']\\\n",
    "    + np.exp(-0.5*ps['k']**2*D**2*ANW) * ( D**2 * f**2 * ps['pl'] + D**4 * f**2 * ps['ptt'] - (D**2*f**2*ps['plnw'] + D**4*f**2*ps['pttres'])\\\n",
    "                                         +0.5*ps['k']**2*D**2*ANW*(f**2*D**2*(ps['pl']-ps['plnw'])))\n",
    "    \n",
    "    #total power spectrum\n",
    "    FoG = (1 - f**2*np.outer(k_array**2,mu_array**2)*D**2*sv_fid)\\\n",
    "             \n",
    "    res = FoG\\\n",
    "    *((b**2*np.outer(pddresum,muones)\\\n",
    "    + D**4*(b*b2*np.outer(bias['ddb0b2nw'],muones) + b**2*np.outer(bias['ddb02nw'],muones)\\\n",
    "    + b*np.outer(bias['ddb0nw'],muones) + b2**2*np.outer(bias['ddb22nw'],muones)\\\n",
    "    + b2*np.outer(bias['ddb2nw'],muones) + np.outer(bias['ddnonenw'],muones)) + N \n",
    "    + 2*(b*np.outer(pdtresum,mu_array**2))\\\n",
    "    - 2*D**4*f*(b2*np.outer(bias['dtb2nw'],mu_array**2) + b*np.outer(bias['dtb0nw'],mu_array**2) + np.outer(bias['dtnonenw'],mu_array**2)))\\\n",
    "    + np.outer(pttresum,mu_array**4))\\\n",
    "    + D**4*( b*f**2*np.outer(abc['bu2f2nw'],mu_array**2) + b*f**3*np.outer(abc['bu2f3nw'],mu_array**2)\\\n",
    "    + b*f**2*np.outer(abc['bu4f2nw'],mu_array**4)+b*f**3*np.outer(abc['bu4f3nw'],mu_array**4)\\\n",
    "    + b*f**3*np.outer(abc['bu6f3nw'],mu_array**6)+b**2*f*np.outer(abc['b2u2f1nw'],mu_array**2)\\\n",
    "    + b**2*f**2*np.outer(abc['b2u2f2nw'],mu_array**2)+b**2*f**3*np.outer(abc['b2u2f3nw'],mu_array**2)\\\n",
    "    + b**2*f**4*np.outer(abc['b2u2f4nw'],mu_array**2) + b**2*f**2*np.outer(abc['b2u4f2nw'],mu_array**4)\\\n",
    "    + b**2*f**3*np.outer(abc['b2u4f3nw'],mu_array**4) + b**2*f**4*np.outer(abc['b2u4f4nw'],mu_array**4)\\\n",
    "    + b**2*f**3*np.outer(abc['b2u6f3nw'],mu_array**6) + b**2*f**4*np.outer(abc['b2u6f4nw'],mu_array**6)\\\n",
    "    + b**2*f**4*np.outer(abc['b2u8f4nw'],mu_array**8) + f**4*np.outer(abc['u2f4nw'],mu_array**2)\\\n",
    "    + f**3*np.outer(abc['u4f3nw'],mu_array**4) + f**4*np.outer(abc['u4f4nw'],mu_array**4)\\\n",
    "    + f**3*np.outer(abc['u6f3nw'],mu_array**6) + f**4*np.outer(abc['u6f4nw'],mu_array**6)\\\n",
    "    + f**4*np.outer(abc['u8f4nw'],mu_array**8) )\\\n",
    "    - 2.*D**2*( np.outer(k_array**2*ps['plnw'],muones)*cs0 + np.outer(k_array**2*ps['plnw'],cs2*mu_array**2)\\\n",
    "    + np.outer(k_array**2*ps['plnw'],cs4*mu_array**4) + np.outer(k_array**2*ps['plnw'],muones)\\\n",
    "    *(f**3*cs0 - f**2*cs2 + f*cs4)*mu_array**6)\n",
    "        \n",
    "    return res\n",
    "\n",
    "#print(Peft_2d(mu,fg_fid))\n",
    "\n",
    "# Now convert to multipoles\n",
    "def Peft(f,b=b_fid,b2=b2_fid,N=N_fid,\n",
    "         cs0=cs0_fid,cs2=cs2_fid,cs4=cs4_fid,sv=sv_fid,\n",
    "         sig8_z=sigma80_fid*Dg_fid_i(zbin),\n",
    "         ps=ps_fid,abc=abc_fid,bias=bias_fid):\n",
    "    \n",
    "    pk2d = Peft_2d(mu,f,b,b2,N,cs0,cs2,cs4,sv,sig8_z,ps,abc,bias)\n",
    "    \n",
    "    return elgc.pkto1d(ps['k'],mu,pk2d,MLPS)\n",
    "\n",
    "k_test, plk_test = Peft(fg_fid)\n",
    "print (len(k_test),len(ps_fid['k']))\n",
    "print ('k are preserved:',(elgc.get_k_once(k_test)==ps_fid['k']).all())\n",
    "#plt.plot(k_test,plk_test,'.')\n",
    "#plt.show()\n",
    "#del k_test, plk_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(900, 900)\n",
      "[[1106.96     0.       0.    ...    0.       0.       0.   ]\n",
      " [   0.      27.222    0.    ...    0.       0.       0.   ]\n",
      " [   0.       0.       9.367 ...    0.       0.       0.   ]\n",
      " ...\n",
      " [   0.       0.       0.    ...  104.19     0.       0.   ]\n",
      " [   0.       0.       0.    ...    0.     104.973    0.   ]\n",
      " [   0.       0.       0.    ...    0.       0.     105.765]]\n"
     ]
    }
   ],
   "source": [
    "# Save the k-array nicely for later\n",
    "k_once = ps_fid['k']\n",
    "k_long = []\n",
    "for l in MLPS:\n",
    "    k_long = np.concatenate((k_long,k_once))\n",
    "\n",
    "# Effective volume: basically Vsur*Cov^(-1)\n",
    "nbar = 1e-3 #h^3/Mpc^3\n",
    "def Veff(zc):\n",
    "    n = nbar    \n",
    "    Pkmu = P_lin(mu,fg_fid) #linear covariance \n",
    "    C = elgc.covll_reduced(k_once,mu,Pkmu,n,MLPS)\n",
    "    return Vsur(zc) * np.linalg.inv(C)\n",
    "\n",
    "veff = Veff(z_vals[0]) # save for later (we really only need to calculate this once per redshift)\n",
    "print (veff.shape)\n",
    "\n",
    "# This looks like it is twice the size of the 2d case \n",
    "# (but this is compensated for by an extra factor of 1/2 in the Fisher matrix)\n",
    "print (veff) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define how to calculate the derivatives with respect to parameters we wish to vary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Construct the derivatives wrt the bias and EFT parameters\n",
    "\n",
    "#choose the step (convergent derivatives)\n",
    "params_step = 1e-5\n",
    "\n",
    "#Derivative wrt b\n",
    "def dPeft_db():\n",
    "    \n",
    "    P_bp = Peft(fg_fid,b_fid*(1+params_step),b2_fid,N_fid,\n",
    "                 cs0_fid,cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_bm = Peft(fg_fid,b_fid*(1-params_step),b2_fid,N_fid,\n",
    "                 cs0_fid,cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    \n",
    "    deriv_b = (P_bp - P_bm)/(2*b_fid*params_step)\n",
    "    \n",
    "    return deriv_b\n",
    "    \n",
    "    \n",
    "#Derivative wrt b2\n",
    "def dPeft_db2():\n",
    "   \n",
    "    P_b2p = Peft(fg_fid,b_fid,b2_fid*(1+params_step),N_fid,\n",
    "                 cs0_fid,cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1] \n",
    "    P_b2m = Peft(fg_fid,b_fid,b2_fid*(1-params_step),N_fid,\n",
    "                 cs0_fid,cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1] \n",
    "    \n",
    "    deriv_b2 = (P_b2p - P_b2m)/(2*b2_fid*params_step)\n",
    "    \n",
    "    \n",
    "    return deriv_b2\n",
    "\n",
    "#Derivative wrt N\n",
    "def dPeft_dN():\n",
    "    P_Np = Peft(fg_fid,b_fid,b2_fid,N_fid*(1+params_step),\n",
    "                 cs0_fid,cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_Nm = Peft(fg_fid,b_fid,b2_fid,N_fid*(1-params_step),\n",
    "                 cs0_fid,cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    \n",
    "    deriv_N = (P_Np - P_Nm)/(2*N_fid*params_step)\n",
    "            \n",
    "    return deriv_N\n",
    "\n",
    "#Derivative wrt cs0 \n",
    "def dPeft_dcs0(): \n",
    "    P_cs0p = Peft(fg_fid,b_fid,b2_fid,N_fid,\n",
    "                 cs0_fid*(1+params_step),cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_cs0m = Peft(fg_fid,b_fid,b2_fid,N_fid,\n",
    "                 cs0_fid*(1-params_step),cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    \n",
    "    deriv_cs0 = (P_cs0p - P_cs0m)/(2*cs0_fid*params_step)\n",
    "    \n",
    "    return deriv_cs0\n",
    "\n",
    "#Derivative wrt cs2\n",
    "def dPeft_dcs2(): \n",
    "    P_cs2p = Peft(fg_fid,b_fid,b2_fid,N_fid,\n",
    "                 cs0_fid,cs2_fid*(1+params_step),cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_cs2m = Peft(fg_fid,b_fid,b2_fid,N_fid,\n",
    "                 cs0_fid,cs2_fid*(1-params_step),cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    \n",
    "    deriv_cs2 = (P_cs2p - P_cs2m)/(2*cs2_fid*params_step)\n",
    "    \n",
    "    return deriv_cs2\n",
    "\n",
    "#Derivative wrt cs4\n",
    "def dPeft_dcs4(): \n",
    "    P_cs4p = Peft(fg_fid,b_fid,b2_fid,N_fid,\n",
    "                 cs0_fid,cs2_fid,cs4_fid*(1+params_step),sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_cs4m = Peft(fg_fid,b_fid,b2_fid,N_fid,\n",
    "                 cs0_fid,cs2_fid,cs4_fid*(1-params_step),sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    \n",
    "    deriv_cs4 = (P_cs4p - P_cs4m)/(2*cs4_fid*params_step)\n",
    "    \n",
    "    return deriv_cs4\n",
    "\n",
    "#Derivative wrt f\n",
    "def dPeft_df():\n",
    "    P_fp = Peft(fg_fid*(1+params_step),b_fid,b2_fid,N_fid,\n",
    "                 cs0_fid,cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_fm = Peft(fg_fid*(1-params_step),b_fid,b2_fid,N_fid,\n",
    "                 cs0_fid,cs2_fid,cs4_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    \n",
    "    deriv_f = (P_fp - P_fm)/(2*fg_fid*params_step)\n",
    "    \n",
    "    return deriv_f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the Fisher matrix calculation and run to loop to calculate its elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.245, 0.245, 0.129]\n",
      "900 -> 450 or 450\n"
     ]
    }
   ],
   "source": [
    "k_step = k_once[6] - k_once[5] \n",
    "\n",
    "# k values range to include in the constraints calculation\n",
    "k_index = np.zeros(len(k_long),dtype=int)\n",
    "for il,l in enumerate(MLPS):\n",
    "    k_index+= (k_long<kmax[il]) * elgc.get_l_boolindex(il,k_long)# kmax for monopole\n",
    "k_index*= k_long>kmin   # And all have to be larger than minimum\n",
    "k_index = k_index.astype(bool)\n",
    "\n",
    "#print (k_index)\n",
    "print (kmax)\n",
    "print (len(k_long),'->',np.sum(k_index),'or',len(k_long[k_index]))\n",
    "\n",
    "#plt.plot(k_long[k_index],plk_test[k_index],'.')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n"
     ]
    }
   ],
   "source": [
    "# This dictionary contains all the parameters that can be varied in this code.\n",
    "derivs = {\n",
    "    'b'  : lambda z: dPeft_db(),\n",
    "    'b2' : lambda z: dPeft_db2(),\n",
    "    'N'  : lambda z: dPeft_dN(),\n",
    "    'cs0': lambda z: dPeft_dcs0(),\n",
    "    'cs2': lambda z: dPeft_dcs2(),\n",
    "    'cs4': lambda z: dPeft_dcs4(),\n",
    "    'f'  : lambda z: dPeft_df()\n",
    "}\n",
    "            \n",
    "# Which should depend on redshift\n",
    "zdep = 'b','b2','N','cs0','cs2','cs4','f'\n",
    "\n",
    "# Which ones are independent of redshift and should appear in the FM once (summ over z contributions)\n",
    "zidep = [i for i in derivs.keys() if i not in zdep]\n",
    "print (zidep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(dPeft_df())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Auxiliary bits needed for Fisher matrix calculation\n",
    "\n",
    "# Create array of zeros for the final FM\n",
    "Nzdep = len(set(zdep)&set(varypars))\n",
    "Npar = len(varypars)\n",
    "s = Nzdep*Nzbin + (Npar-Nzdep)\n",
    "Fishermat = np.zeros([s,s])\n",
    "\n",
    "# Fisher matrix integrand - to be integrated over k and summed over l\n",
    "# Do the whole sum together in one big matrix operation\n",
    "def dF(zc,veff):\n",
    "    di = deriv_i(zc)[k_index]*k_long[k_index]\n",
    "    dj = deriv_j(zc)[k_index]*k_long[k_index]\n",
    "    tmp = np.dot(di,veff[k_index,:][:,k_index])\n",
    "    tmp = np.dot(tmp,dj)\n",
    "    return tmp*k_step/(4*np.pi**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%time \n",
    "# We could optimize this still by precalculating the derivatives\n",
    "\n",
    "# Save the indexing for posterity\n",
    "fm_index = {'b' :[],'b2':[],'N' :[],'cs0':[],'cs2':[],'cs4':[],'f' :[]}\n",
    "\n",
    "fm_head = []\n",
    "\n",
    "# Loop over redshift bins\n",
    "for iz, zc in enumerate(z_vals):\n",
    "    \n",
    "    veff = Veff(zc) # Let's not calculate this for every parameter - it stays the same at each redshift!\n",
    "    \n",
    "    # First loop over parameters (rows of the FM)\n",
    "    for i, pi in enumerate(varypars): \n",
    "        \n",
    "        # Calculate row index for the Fisher matrix\n",
    "        if pi in zdep:\n",
    "            zind = iz*Nzdep + i\n",
    "        else:\n",
    "            zind = Nzbin*Nzdep + (i-Nzdep)\n",
    "                \n",
    "        # Set the i derivative in dF to correspond to the row parameter\n",
    "        def deriv_i(zc): return derivs[pi](zc)\n",
    "\n",
    "        # Second loop over parameters (columns of the FM)\n",
    "        for j, pj in enumerate(varypars):\n",
    "            \n",
    "            # Calculate column index for the Fisher matrix\n",
    "            if pj in zdep:\n",
    "                zjnd = iz*Nzdep + j\n",
    "            else:\n",
    "                zjnd = Nzbin*Nzdep + (j-Nzdep)\n",
    "            \n",
    "            # If not done yet, set the j derivative in dF to correspond to the column parameter\n",
    "            if j>=i:\n",
    "                \n",
    "                def deriv_j(zc): return derivs[pj](zc)\n",
    "                \n",
    "                Fishermat[zind,zjnd] += dF(zc,veff)\n",
    "            \n",
    "            # Otherwise just take the symmetric value\n",
    "            else: \n",
    "                Fishermat[zind,zjnd] = Fishermat[zjnd,zind]\n",
    "        \n",
    "        # Finally, save the index in the FM (same redshift order as z_vals)\n",
    "        fm_index[pi].append(zind) # Only doing i, since it is symmetric\n",
    "        \n",
    "        # Save also the order of the colmns/rows\n",
    "        headstr = pi\n",
    "        if pi in zdep: headstr+='_'+str(zc)\n",
    "        fm_head.append(headstr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Print out results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
       "    return false;\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
    "    return false;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a fiducial value dictionary\n",
    "fids = {\n",
    "    'b'  : b_fid*np.ones(Nzbin),   \n",
    "    'b2' : b2_fid*np.ones(Nzbin),  \n",
    "    'N'  : N_fid*np.ones(Nzbin),   \n",
    "    'cs0': cs0_fid*np.ones(Nzbin), \n",
    "    'cs2': cs2_fid*np.ones(Nzbin), \n",
    "    'cs4': cs4_fid*np.ones(Nzbin), \n",
    "    'f'  : fg_fid*np.ones(Nzbin)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Order of columns/rows:\n",
      "['f_0.5', 'b_0.5', 'b2_0.5', 'N_0.5', 'cs0_0.5', 'cs2_0.5', 'cs4_0.5']\n",
      "\n",
      "Correlation matrix for parameters:\n",
      "\n",
      "[[ 1.    -0.451 -0.506  0.344 -0.567  0.131  0.211]\n",
      " [-0.451  1.     0.967 -0.987 -0.253  0.412  0.391]\n",
      " [-0.506  0.967  1.    -0.965 -0.082  0.383  0.287]\n",
      " [ 0.344 -0.987 -0.965  1.     0.329 -0.456 -0.426]\n",
      " [-0.567 -0.253 -0.082  0.329  1.    -0.564 -0.401]\n",
      " [ 0.131  0.412  0.383 -0.456 -0.564  1.    -0.328]\n",
      " [ 0.211  0.391  0.287 -0.426 -0.401 -0.328  1.   ]]\n",
      "\n",
      "Correlation matrix diagonal:\n",
      "\n",
      "[5.741e-04 1.147e-03 4.537e-01 1.075e+06 6.652e+00 3.441e+01 8.465e+00]\n"
     ]
    }
   ],
   "source": [
    "print ('Order of columns/rows:')\n",
    "print (fm_head)\n",
    "#print 'Fisher matrix:\\n',\n",
    "#print Fishermat\n",
    "\n",
    "# Invert the Fisher matrix to get the parameter covariance\n",
    "CovMat = np.linalg.inv(Fishermat)\n",
    "print ('\\nCorrelation matrix for parameters:\\n', )\n",
    "print (CovMat/np.sqrt(np.outer(np.diag(CovMat),np.diag(CovMat))))\n",
    "\n",
    "print ('\\nCorrelation matrix diagonal:\\n', )\n",
    "print( np.diag(CovMat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "varying zdep = \n",
      "\t ['f', 'b', 'b2', 'N', 'cs0', 'cs2', 'cs4']\n",
      "at z = 0.5 : zdep fiducials = \n",
      "\t [ 7.042e-01  1.468e+00 -2.756e-01  8.161e+02  7.650e+00  2.978e+00  2.440e+01]\n",
      "at z = 0.5 : marged err = \n",
      "\t [2.396e-02 3.387e-02 6.736e-01 1.037e+03 2.579e+00 5.866e+00 2.909e+00]\n",
      "at z = 0.5 : unmarged err = \n",
      "\t [2.771e-03 1.092e-03 5.724e-03 7.241e+00 6.360e-02 2.516e-01 1.689e-01]\n",
      "at z = 0.5 : marged % err = \n",
      "\t [  3.403   2.307 244.396 127.016  33.713 196.991  11.924]\n",
      "at z = 0.5 : unmarged % err = \n",
      "\t [0.394 0.074 2.077 0.887 0.831 8.447 0.692]\n"
     ]
    }
   ],
   "source": [
    "# Print the errors nicely\n",
    "Nvar  = len(varypars)\n",
    "\n",
    "varyshape = list(set(zidep)&set(varypars))\n",
    "if len(varyshape)>0:\n",
    "    shapeind = [fm_index[p][0] for p in varyshape]\n",
    "    shapefid = np.array([fids[p] for p in varyshape])\n",
    "    print ('\\nvarying shape = \\n\\t', varyshape)\n",
    "    print ('shape fiducials = \\n\\t', shapefid)\n",
    "    print ('unmarged err = \\n\\t', np.sqrt(1./np.diag(Fishermat))[shapeind])\n",
    "    print ('marged err = \\n\\t', np.sqrt(np.diag(CovMat))[shapeind])\n",
    "    print ('unmarged % err = \\n\\t', np.sqrt(1./np.diag(Fishermat))[shapeind]/np.abs(shapefid)*100)\n",
    "    print ('marged % err = \\n\\t', np.sqrt(np.diag(CovMat))[shapeind]/np.abs(shapefid)*100)\n",
    "\n",
    "#varyzdep = list(set(varypars)&set(zdep))\n",
    "varyzdep = varypars\n",
    "\n",
    "zdepfid = np.array([fids[p] for p in varyzdep])\n",
    "print ('\\nvarying zdep = \\n\\t', varyzdep)\n",
    "if len(varyzdep)>0:\n",
    "    for iz, z in enumerate(z_vals):\n",
    "        zdepind = np.array([fm_index[p] for p in varyzdep])\n",
    "        print ('at z =',z,': zdep fiducials = \\n\\t', zdepfid[:,iz])\n",
    "        merrs = np.sqrt(np.diag(CovMat))[zdepind[:,iz]]\n",
    "        print ('at z =',z,': marged err = \\n\\t', merrs)\n",
    "        uerrs = np.sqrt(1./np.diag(Fishermat))[zdepind[:,iz]]\n",
    "        print ('at z =',z,': unmarged err = \\n\\t', uerrs)\n",
    "        print ('at z =',z,': marged % err = \\n\\t', merrs/np.abs(zdepfid[:,iz])*100.)\n",
    "        print ('at z =',z,': unmarged % err = \\n\\t', uerrs/np.abs(zdepfid[:,iz])*100.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEaCAYAAAAYOoCaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3U9s22aeN/CvnM7MYsaifBhgsdMwp0WnERVg3kEzA9E9NWks+xYPYuXUtEiUAEX9B1Pr8k4io25OlQeIiy6KRmnj3MwYzS2WnKYvdoGKwjrAHiakMj0tanoW2FMtyoedf+Z78Cu+liXZ1B9KlPj9DARYFCn+7E6eH5/f8/BhwLIsC0RE5DtDvQ6AiIh6gwmAiMinmACIiHyKCYCIyKeYAIiIfIoJgIjIp17qdQBULZVKIZvNwjRNxONxBINBAEC5XIZhGFBVFbIs4/79+8fur2kadF3H/Pw8EokEDMPA+fPnIUkSIpEIgsEgisUiVFW1jy2Xy1BVFaVSCc+ePXMUs2ma+Oyzz3Dq1ClsbW2hXC7bMY2NjWFiYgKxWKzh8blcDs+fP8fIyAgAQBRFSJIEXderjnO6X79QVRWGYUAURft3EEXxyGNM00Q2m4Usy/ZxiqJgfn4egiB0KXIaGBZ5zkcffWRdvHix7mdbW1s1nx21/927d62PPvrIsizLyufz9s8V+XzeeuWVV2qOa/R99bz99tvW1taWtbW1Zd26devYmA5aXV217t69WxPTxYsXrdXV1ab36xdbW1vW9PR01ba333772OM0TbNeeeUV+3Xu3DlL0zS3wqQBxx5AnxFFEdFo1PH+iUQC6XQawP7V4+XLlx0dF4/HYZqmo6tKTdPsK9fFxUXHsQFAJpPB06dPq7bJsoxIJNLSfv0ik8nU/LcIh8PI5XLH9maePn2KUCgEALzqp7ZwDKAPnTlzpqn9T506BQAolUrHlhgqRFFEqVRqOrZmmKbZ8LN4PN70fv0km83WJK8zZ85AVVVHxwuCwMaf2sYE0CcURbF/lmX52P0Nw7B/rlwtNnO1HIlE7OPcIggCDMOArus1n1XGKZrZr18YhlG3dyUIguMEQNQJTAB9Ip/P2z87ufI72JBUSgqSJDk+X+UKM5fLQVEU5HI5ZDKZqsSiqirS6TRM00Q6nbZLTc24du0a3n77bSiKUnOlfzBep/s5lU6ncfbsWaRSKQD7vYyDSdZNjXo0oVCo6u/biK7ryOVy9n+TeomRyAmOAXiUYRh2g7qxsdF0OSafz7ddHqk0tolEwt42MzODGzduQJIkyLIMWZbx8OFDJJPJls5ROS6VSiGVSkGSJIyPj1eds5n9nFBVFeVyGR9++KH9dy4UClhZWTnyuHfeecfxOcLhcMO/STultUryPzhOMDk5ieXlZcflPaIKJgCPEkXRbkCSySRmZmaO3L/SkFWmcbZbHzZNE0tLSzVTQS9fvoxbt27h0aNHbX3/QclkEslkEqqqIpfL4e7du8hmszXncLrfcSKRSFUZLZVKYXl5+di/WWXqbS+JoljT0EejUWQymaYH4IlYAuoTo6OjVe8PlwoqCWNxcRFPnz5tuy5+cGbPQbIsd6TkUK8MIssyFhcX8fXXXwNAw5KS0/0aOdjQp9NpR/PvO6nR2EqrPYNTp05x7IBawgTQJ8bHx6veH/cP3slA8VGOa+TbSQKmaULTNJimWfd7BEHA9evXUSgU7P2d7NesdDqNiYmJtv9Wzaokm8NJ0Mm020wmU7MtFAq5PmOLBhNLQH3iYMNgmuaxM3TavTO20RVxpdFq54rZMAy70VJVte4grizLuHv3LgA43q8Z6XQao6OjTQ0gd2oMoDLAXiqVqv67lkqlI+/xMAwDS0tLiMfjNcex/k+tYALoQ/XmkHeaLMu4detWzXZN0yBJUltjDOvr6/ZNUNlstu5ArmEYVb+j0/2cqDT+B6/8M5nMsQPKnRwDiEaj0HW9quHWdb2m1HeQKIpYXFys+dvn8/mmbg4kqmAJqM8YhoFMJuP6FV+lvHK45JDJZPDhhx/a74+6SaseVVVx7969qjJIvRr+0tIS5ufnq85z3H6GYeDs2bNHlsdyuRyKxSLy+TwymQxyudyxA+xuSCaTNdNONU2rmrlVWbvp8O9zcPxH13UUi8WWZ2GRvwUsi88E9hIni8EBwLffflt3f1mWHZV/KouIVRYkk2UZ4XAYN27cqLrCzOVyMAwDgiDANE3IsmyXTSozjzY2NhCPxxGLxSDLcsPfobLwXCX+yliAKIpYXV21z1kul5FIJKqShJP9DMPA5OQkpqamGjaIlSt90zQxOzsLTdNw/fr1lqaTtktVVbsX8Pz5c1y+fLkqsVd+n+Xl5areSmV67s7ODsrlMheCo5YxAdDAURSlb5eIIOomloCIiHyKCYAGSuV5CUR0PCYAGiimaXJKJJFDHAMgIvIp9gCIiHzKszeCvZr4Wa9DqJF/Y73XITT01+93ex1CXf85+W2vQ6jr5bv/1OsQGvpJ8nSvQ6jrf/a+73UIDb38k1+2dXwz7c0fM//V1rm8xLMJgIioWwI+rYUwARCR7w0NBXodQk8wARCR7wX82f4zARARDbEERETkTyde8mcXgAmAiHyPg8BERD415NNBACYAIvI99gCIiHyKg8BERD7l1x6Aq7/22tqam19PRNQRJ04EHL8GSVs9gHv37qFUKtX9zLIsFAoFXLp0qZ1TEBG5jiWgFsiyjNXVVYyOjtb9/ODDq4mIvCrApSCaFw6HMTo6irGxsbqf80HVRNQP2ANoUaPGHwCi0Wi7X09E5Dq/DgK3lAAMw8Ds7CxevHgBQRAQiUSQTCbx6quvdjo+IiLXuTm4qygKQqEQAKBUKiEejzs6Bth/xCkAJBIJV2JrKQEoioLbt28jGAyiWCwin8/jrbfewuXLl/Hb3/620zESEbnKrRJQpfGPxWIAgFwuB0VRjkwCmUymqsE3TRPpdBrJZLLj8bX0a585cwbhcBiiKGJsbAyLi4vY3NzE3t4ep34SUd8JDAUcv5qhKIrd+ANALBazr+4bOTx5RhAElMvlps7r1LEJIJVK4Te/+Q0WFhbw5MkT7O7uIhAIYHe39hGE8/PzdpeFiKhfDA05fzllmmbdmZC6rh/ZTmqa1rUZlMf+OouLi7AsC5Ik4fHjx3jjjTdw8+ZNnDt3Dmtra9je3q7a//vvvfvcUCKiegIB5y+nNE2za/8HCYIATdMaHjc/P4/JyUm7p5DL5TA/P9/07+SEozGA27dvIxwOY2pqCsB+F0VVVXzzzTdIp9Mol8sQRdEOnoionzTzSEjTNOtewQuCUDX13TTNulPh6yWFg2RZxsrKCiYnJ5HJZHD//n3XptQ7SgDhcLjqvSiKiMfj9kBGpbtSSQJERP2kmQfCPHjwAJ988knN9vfeew/T09NV21opiVcusJ89e4alpSWcP38e9+/fhyzLTX/XcTqyGBwbfiLqZ808D+DKW1dw8eLFmu2Hr9IbXbU3Wj6nIpPJYHFxEcB+CT4Wi2F2dhZff/11x3sCXA2UiHyvmRLQ4VJPI5FIpG5jb5omIpFI3WNUVa250pdlGdevX4eqqlUzijqBCYCIfK+ZBOCUIAgIhUI1YwGiKDZMII1KRpIkdTw+wOXloImI+sFQYMjxqxnxeBzZbNZ+ryhK1U1euq5jZmbGfh+LxbC+vl7zPfl83rtjAERE/cyNHgCwv4SDoihQVdW+uj94F3CpVEKhUKjqJSSTSaTTaYyMjEAQBJimiRs3brgSHxMAEfneSyfcK4YcteyDLMt49uxZ1TZRFF1Z9qEeJgAi8r0hn64H7dkEkH+jtg7mBaP/Z6LXIdT1wx9584EWV6c7X7fshHf/8C+9DqGhr/77s16HUNdLQz/odQgNvfyTX7Z1vFslIK/zbALwIq82/kTUnmbuAxgkTABE5HvsARAR+dSJEyd6HUJPMAEQke+xB0BE5FNMAEREPtXsHb6DggmAiHyPPQAiIp/iNFAiIp966SXOAiIi8iW/9gA6MvKxvb2NJ0+e4MWLF/a2crmMQqFQ89B4IiKvGRoKOH4NkrYTwMbGBiYnJ/H48WPMzMxgbm4OABAMBiFJEt588822gyQicpNfE0DbJSBVVbG5uWm/13UdCwsLSCaTEAQBlmW1ewoiIlf5dRpo27/14UeVSZKEDz74AIqioFwuI+DT2hoR9Q+/9gDaTgChUAiFQgFzc3PY3d21t1+9ehWaprEHQESe99KJHzh+DZK2S0BjY2PY3t7GxMQEhoeHqz6LRqN49OhRu6cgInLViQCngTbFMAzMzs7ixYsXEAQBkUgEp06dwquvvlq1XzgcbjtIIiI3DQ0xATRFURTcvn0bwWAQxWIR+Xweb731Fi5fvozf/va3nYyRiMhVQz7tAbQ8BnDmzBmEw2GIooixsTEsLi5ic3MTe3t7WFtb62SMRESuGhoacvwaJI56AKlUCrquIxKJYHR0FLIsIxAIYHd3t6buPz8/j88//9yVYImI3PCDEz/sdQg94SidLS4uwrIsSJKEx48f44033sDNmzdx7tw5rK2t1dzt+/3337sSLBGRG4YCQ45fg8TxGMDt27cRDocxNTUFYH8QWFVVfPPNN0in0yiXyxBFEcB+L4CIqF9wEPgYh2fziKKIeDyOeDwOYD8hVLYTEfUTTgNtExt+IupXg1bacYrLQROR77EERETkU36dBcQEQES+xxIQEZFPsQRERORTbi4FoSgKQqEQAKBUKtkzJ51KpVJYXFx0I7TOPBKSiKifubUURKXxj8ViiMViCIVCUBTF8fG5XA6qqjb76zjm2R7AX7/fPX6nLvvX//UQF4rNZe9u+cufvfnchR95dP1062c/6XUIDf1w9x96HUJd//jjV3odgmvcug9AUZSqJfFjsRgmJycd9QJM03QlpoM8mwC8yKuNPxG15yUXZgGZpmnfIHuQruswTROCIBx5vKqqkGW543EdxBIQEfmeG2sBaZpm1/4PEgQBmqYdeayu6643/gB7AEREONHELCDTNOuWZwRBqLqqb3SVXy8pHFYqlSBJkutlICYAIvK9ZmYBPXjwAJ988knN9vfeew/T09NV21ppwHO5HGKxWNPHtYIJgIh8r5nSzpUrV3Dx4sWa7Yev9hvV+EulUsPvNgwDkiQ5jqVdTABE5HsvDTkfBD5c6mkkEonUbexN00QkEql7jKqqVb2GnZ0dlEolZDIZSJLU8XEBJgAi8r2AC9NABUFAKBSqGQsQRbFhAjk8PdQwDGxsbCCRSHQ8PoCzgIiIMBQ44fjVjHg8jmw2a79XFKWqMdd1HTMzM0d+x1Elo3YxARCR7w3hhONXMyqNvaqqyOVyAKqv8kulEgqFQt3B4lwuh1QqBdM0kU6noet6G79hfS2VgD7//HNcvXq107EQEfWEm2sBHXXXryzLePbsWd3PKstHuKmlHoCiKNjd9d5SDURErXCrBOR1LfUAdnZ2oCgKgsEgAoEAIpEITp8+3enYiIi64kSAD4RxbH5+HlNTU/b7YrGIhw8fwjRNhEIhnDx5EtFotGNBEhG5iQ+EacLBxh8AwuEwwuGw/b5QKODChQsIBoOYmJjgeAERedqglXac6uh9AIVCAYqiIJfLQRAEXLhwARcuXOjkKYiIOo4JoEW7u7tQFAWKomBrawvhcBjLy8sYGxvrRHxERK5z40awftByAqhc7W9sbCAYDGJqagrxeByiKHYyPiIi173EQWDnLly4AMMwEI1GcefOHV7tE1FfYwmoCVtbW7h27RomJiaqBn+JiPoRE0ATpqamMD8/b0//BNDwfoBCocApoUTkaRwDaEJlfYvD0z8PJwRJkrC0tIQvv/yyA6ESEbmj2TV+BkVLCaDRQO/hhLCxsYHt7e3WIiMi6hLeCOaCsbExrK+vu3kKIqK2NfNAmEHi+gNhbt++7fYpiIjaEmAJyB3BYNDtUxARtYWzgIiIfIoJgIjIpwI+fTgiEwAREQK9DqAnmACIyPd4HwARkW+xBERE5EsBloC85T8nv+11CDWuTsu9DqGhH534Qa9DqOvOlX/tdQh1/eIva70OoaFf/Nu5XodQ14kf/6jXITQ22d7hAd4JTETkV+wBEBH5EgeBiYh8iyUgIiJf4iAwEZFP8U5gIiLfYg+gLYVCAcVi0X4fDof5KEgi6gvsAbRoY2MD6+vrEEURIyMj9vZ8Pg9FURCPx5kIiMjT+DyANiwvLzf87MmTJ504BRGRa9gDaFEg4M/aGRENDjdnASmKglAoBAAolUqIx+OuHNOKthOAZVmYnZ3FyMgIBEGwy0A7OzswDMO1wImIOsalpSAqDXksFgMA5HI5uzTeyWNa1fZvPTY2huXlZciyjFAoBMuyYFkWIpEI7ty5w/o/EXleoIn/NUNRFLshB4BYLAZFUTp+TKs6NgtobGys7vYXL14A2H828MmTJzt1OiKijnFjDMA0TRiGUbNd13WYpglBEDpyTDtcHfkwDAO/+93vcPr0aViWhUKh4ObpiIhaEsAJxy+nNE2z6/gHCYIATdM6dkw7XL0RTBRFPHr0yP5ZFEU3T0dE1JJmSjumacI0zZrtgiBUXaE3umKv18C3c0w7GiaA7e1tFItFiKKI06dPAwDK5TI0TYMoiiznENEAcZ4AHjx4gE8++aRm+3vvvYfp6emqbfUSxXFaOaZVdRPAxsYGbt26hWg0imKxCEmScOfOHQSDQUiShF//+td2bf/wcZqm4caNG/jd736HU6dO4caNGxgeHnb9FyEiapnlfNcrV67g4sWLNdsPX7k3qteXSqWG393KMe2omwBUVcXm5qb9Xtd1LCwsIJlMQhAEWFb9v5YgCHj//fexsLCA119/HZcuXcLa2houXbrkSvBERJ0QaNCm1XO41NNIJBKp23CbpolIJNKxY9pRdxBYkqSa9x988AEURUG5XG5485dhGNjY2EA2m8X4+DgAsFRERN6318TLIUEQEAqFako6oig2TCCtHNOOugkgFAqhUChgbm4Ou7u79varV69C07SGPYCpqSmIooiVlRUMDw9jaWmp4wETEXWcZTl/NSEejyObzdrvFUVBIpGw3+u6jpmZmaaO6aS6CWBsbAwnT57E+Ph4Tf0+Go3i008/rdp2MEmEw2GEw2EAwOjoKM6cOdPpmImIOipgOX81o9Jwq6qKXC4HAFV39JZKJRQKhaor/uOO6aSGs4AaTdt8+PAhFhYW8OjRI3t20PDwcN1aP+8CJqK+0GTD3oyjGm9ZlvHs2bOmjumkpm8EC4VC2NzctBv/ilgsxpU/iag/uVQC8rqmbwQzTRPBYLBmezAYxM7OTkeCIiLqqsFq1x1rugdgWVbDK30uDU1E/ShgWY5fg6TpBDA1NYXV1VXMzc3hq6++srfv7u7iu+++62hwRERdYTXxGiAtrQX0xRdfIJPJYHp6GoFAwJ6furKy0snYiIi6Y8AadqdaXgwukUggkUjAMAyUy2V76icRUd8ZsNKOU22vBsoVPomo3zU7v39QuLocNBFRX2AC8JaX7/5Tr0Oo8e4f/qXXITRk/ewnvQ6hrl/8Za3XIdT19s3/3esQGnoi3ut1CHWNvPbPvQ7BPXv+zACeTQBERN3i1xKQq4+EJCIi72IPgIiIs4CIiHzKn+0/EwARUYCDwEREPuXP9r8zg8Db29t48uRJ1YPiy+UyCoUCtre3O3EKIiL37FnOXwOk7QSwsbGByclJPH78GDMzM5ibmwOwvzy0JEl488032w6SiMhNlmU5fg2StktAqqpic3PTfq/rOhYWFpBMJiEIwsD9wYhoADXxsPdB0nYPQJKkmvcffPABFEVBuVzmMwKIyPOsPcvxa5C0nQBCoRAKhQLm5uaqHg5/9epVaJrGHgAReZ719z3Hr0HSdglobGwMhmFgfHwcw8PDVZ9Fo1E8evSo3VMQEblrwK7snerINFBRFBsuCx0IBPDixQsEg0GcPHmyE6cjIuoov1YqXF0LyDAM3Lx5E6dPn4ZlWSgUCm6ejoioNXtNvAaIqzeCiaKIL7/80v6ZD48hIi/yaw+gYQLY3t5GsViEKIo4ffo0gP2buzRNgyiKLOcQ0cCw/jZgl/YO1S0BtXpz18bGBn7/+99jd3cXs7Oz9s9ERJ5mWc5fA6RuD6DVm7sEQcD777+PhYUFvP7667h06RLW1tZw6dIld6InIuqAQZvf71TdHkCrN3cZhoGNjQ1ks1mMj48DAEtFROR9Ph0ErpsAWr25a2pqCqIoYmVlBcPDw1haWnInaiKiDuJaQAccd3PXp59+WrVtd3fX3i8cDtvbR0dHcebMmU7HTETUWSwBVRNFEWNjYzXbHz58iHfffbdq6efh4WGsra3V7BuNRmsSCBGR1/h1KYimbwQLhULY3Ny0p4ZWxGIxPHnypGOBERF1jU/HAJq+Ecw0TQSDwZrtwWAQOzs7HQmKiKibvDYLSFEUhEIhAECpVEI8Hm/q+FQqhcXFxWP3a7oHYFlWwyt9Lv1MRH3JQ/cBVBr/WCyGWCyGUCgERVEcH5/L5aCqqqN9m04AU1NTWF1dxdzcHL766it7++7uLr777rtmv46IqOe89DwARVEQi8Xs97FYzHECME2zqXO1tBbQF198gUwmg+npaQQCAQiCAABYWVlp5euIiHqriYbdNM26Da0gCHZb2CrTNGEYRs12Xddhmuax36+qKmRZdny+lheDSyQSSCQSMAwD5XK5avonEVE/aWZ2z4MHD/DJJ5/UbH/vvfcwPT3dVhyaptm1/4MEQYCmaUc27rquN9X4Ax1YDZQrfBJRv2smAVy5cgUXL16s2d7u1T+Ahlf59ZLCYaVSCZIkNVUGcnU5aCKiftBMbT/UgVLPUZqt4wP7A78Hxw2c8mwC+Eny9PE7ddlX//1Zr0No6Ie7/9DrEOr6xb+d63UIdT0R7/U6hIYuGNd6HUJd//Hmv/c6BNe4cYOXqqrIZDLH7hcMBvHxxx8DaNyLKJVKDY83DKNm/TanPJsAiIi6Zq/zCUCW5aZr8pFIpG5jb5omIpFI3WNUVa3qNezs7KBUKiGTyUCSpCNjYAIgIt/b++vfex0CgP0eQCgUqhkLEEWxYe/g8E1ilVWZE4nEsedz9ZnARET9wEtrAcXjcWSzWfu9oihVjbmu65iZmTnyO44qGR3EBEBEtLfn/OWySmOvqipyuRyA6qv8UqmEQqFQd7A4l8shlUrBNE2k02noun7kuVgCIiLf89paQEet/SPLMp49e1b3s8ryEU4xARCR7w3aMs9OuVYCKhQKAFD13AAiIi+y9vYcvwZJ2z2Ahw8folwu12x//vw5isUiVFXF559/3u5piIhcY3lkFlC3td0DiEQi+MMf/gBBEBAMBu3XoD5Dk4gGj5dmAXVT2z2AcDiM5eVlFAoFjIyMVD0pbGpqqukbIYiIum3QSjtOdWwQOBqNolwu48mTJwiHw/bDYbhKKBF53t/9Wano6CBwMBjEhQsXYFmWq4slERF1EgeBO0gURXuZ6MosoGAwiJMnT7pxOiKitgxabd8pV+8ENgwDN2/exOnTp2FZlj01lIjIS6y//t3xa5C4eiOYKIr48ssv7Z/58Bgi8qJBK+041TABbG9vo1gsQhRFe2ZPuVyGpmkQRZHlHCIaGCwBHbCxsYHJyUk8fvwYMzMzmJubA7Bfx5ckCW+++WbdL9vY2MDvf/977O7uYnZ21v6ZiMjLOAh8gKqq2NzctN/ruo6FhQUkk0kIgtDw5i5BEPD+++9jYWEBr7/+Oi5duoS1tTVcunTJneiJiDqB00D/v8OPF5MkCR988AEURUG5XLbn+B9WeRBBNpvF+Pg4ALBURESet/eXvzl+DZK6CSAUCqFQKGBubq6qhHP16lVomtawBzA1NQVRFLGysoLh4WEsLS25EzURUQdxKYgDxsbGYBgGxsfHMTw8XPVZNBrFp59+WrVtd3fX3u/gnb+jo6M4c+ZMp2MmIuoo62+D1bA71fA+AFEUMTY2VrP94cOHePfdd6uWeR4eHsba2lrNvtFotCaBEBF5jV97AE3fCBYKhbC5uVm16Buw/ySaJ0+edCwwIqJusf625/g1SJq+Ecw0TQSDwZrtwWAQOzs7HQmKiKibBq1hd6rpHoBlWQ2v9BvNDiIi8rK9P//d8WuQNJ0ApqamsLq6irm5OXz11Vf29t3dXXz33XcdDY6IqBv8OgbQ0lpAX3zxBTKZDKanpxEIBOyln1dWVjoZGxFRV/i1BNTyYnCJRAKJRAKGYaBcLvPBL0TUt5gAWsQVPomo3w1aaccpV5eDJiLqB+wBeMz/7H3f6xBqvDT0g16H0NA//viVXodQ14kf/6jXIdQ18to/9zqEhv7jzX/vdQh1/fKLX/c6hIb+GP2vto7f+/NgrfHjlGcTABFRt7AHQETkUxwDICLyKfYAiIh8igmAiMinvLbEg6IoCIVCAIBSqYR4PO7oGGB/vTZg/16t4zABEJHveWkMoNL4x2IxAEAul4OiKEcmgUwmU9Xgm6aJdDqNZDJ55LmaXguIiGjQeGk5aEVR7MYf2F9qv3J134hhGFXvBUFAuVw+9lxNJ4BCoVD1MBgion7nlQRgmmZNYw4Auq7bpZ16NE2re9xxmk4A2WwWn332mf3+4DODiYj6UTOrgZqmie3t7ZrXUQ20U5qm2bX/gwRBgKZpDY+bn5/H5OSk3VPI5XKYn58/9nxNJ4BTp07hzp079vuDyYCIqB81kwAePHiAc+fO1bwePHjQdhymadqrKx9ULykcJMsyVlZWkEqlcP78eUiSVPd7Dmt6EFiSJPzqV7+CKIqQZRmGYeBPf/oTXn755Wa/iojIE/6253wW0JUrV3Dx4sWa7U4aXCda6UkYhgFVVfHs2TMsLS3h/PnzuH//PmRZPvK4phNANBrF5uYmCoUC8vk8VFXFuXPnEAgEEA6HEYlE7JcoinwoPBF53t8t57V9QRAcNfaqqiKTyRy7XzAYxMcff2x/dz2lUunI78hkMlhcXAQALC4uIhaLYXZ2Fl9//fWRsbY8DTQajSIajQLYrz8ZhoFisYjnz5/jm2++wd27d7G9vY1wOIyJiQnE43EmAyLypL0mEoBTsiwfewV+WCQSqdtWZig6AAAIBUlEQVTYm6aJSCRS9xhVVWvOI8syrl+/DlVVq2YUHdb2fQCVuamiKEIURYyNjVV9XiwWoaoqbt68iYmJCVy4cKHdUxIRddSeZfU6BAD7PYBQKFQzFiCKYsMr+UYlI0mSjj1f2/cBHPdAmHA4jGvXriEejx/bjSEi6oU9a8/xy23xeBzZbNZ+ryhK1U1euq5jZmbGfh+LxbC+vl7zPfl8vvNjAK0ol8t45513HN2aTETUbc0MArstkUhAURSoqmpf3R+8C7hUKqFQKFT1EpLJJNLpNEZGRiAIAkzTxI0bN449V1cSQDAYxB//+MdunIqIqGnduLJvxlHLPsiyjGfPnlVtE0Xx2GUf6uFaQETke83MAhokTABE5HteGQTuNiYAIvI9r5WAuoUJgIh8jwmAiMinvDQLqJuYAIjI9zgITETkUywBERH5FGcBERH5FHsAREQ+xQRARORTf+UsICIif2IPgIjIp/w6CBywLJ/+5kREPtf2A2GIiKg/MQEQEfkUEwARkU8xARAR+dRAJoBUKtXrEIiIPG+gEoCu65icnISmab0OpS4mJiLykoG6D0CSJMTjcei63utQqui6jlu3bvU6DCKiKgOVAAAgn89jYmKi12FU8Wpi6ge5XA6lUgm6rkOWZcRisV6HZPNqbF6Ni7xn4BJAoVDA7du3ex1GDS8mJq8zDAOCINgN2NmzZyFJEkRR7HFk3o3Nq3GRNw3UGIBhGAAAVVWRy+WQTqftbb1WKBQgy3Kvw6jLMAykUin8/Oc/x/nz52s+VxQFZ8+exfnz55HJZLoWl67rWF1dtd9Ho1HP9KK8GptX42pEVVVkMhmkUimoqtrrcPzHGiCrq6vWa6+9ZpVKJcuyLEvTNOvixYs9jsqytra2rNdee83KZrNWNpu1PvroI2tra6vXYdX46KOPrFdeecXK5/M1n929e9f+u3bTwXOeO3fO0jSt6zE04tXYvBpXPa+99pplWZY1PT3t6TgH1UD1APL5PObn5yEIAgBAEISqqx9d17t6BVtRubKp1GMnJiYwOzvb9TiOous6JiYmIMty1RVkhSAI9t+1myrn1HUd4XAYkiR1PYZGvBqbV+M6TFVVuzT18ccfezbOQTZQCeBwmUXXdfsfQy6Xg6Io2NnZ6XpcRyWmSlypVAq5XK7rsVVommYPVm9sbMA0zZ7FUo+iKPj44497HUZdXo3Nq3FV5PN5RKPRXofhawOTAAzDgGmaVYNdqqoiEokAAGKxWM9q8I0SU2XALh6PY3FxEbdu3er5mEUsFoMgCFAUxd6m67r9d+yFTCaD+fl5APBcYvJqbF6NC9j/t5pOp/Hw4UP7Z6/F6BcDNQvo8EyHbDaLlZWV3gTz/xyVmHRdx/r6up0cKgN2vZ6xMTU1BUVRkEgkAOz3DuLxuOvnrSSdUCiEUqmEeDwOVVURj8chCAJM04SmaT1J5F6NzatxHUUURSSTSdy7dw/JZLLn/3/3s4FJAIf/T5ROpzE1NeWJumKjxCRJUtU/zGKxiBs3bnQ7vJor/MuXL+PevXtQVRWyLHfl6uydd95BPB5HLBaDaZo4e/YsfvrTn+Ldd9+t2u/p06eux9IvsXk1LicqPV02/r01MAkAAObn55FOpzEyMoJTp0515ar1OMclJi8M2B2+whdFEZIkYXV1FbIsuz74m8lkEAwG7bnrgiDg0aNHkCQJ3377ravn7tfYvBqXU7que+LizO8GKgHEYjFP3vXoJDF5bcDu+vXrmJ2drZqp4ZalpSUsLy9XbfNK4+DV2Lwal1PPnz/v6bgS7RuoBOBVxyWmwwN2vZhueVhlMDiVSrlaQqiUAtpJ3IdnT3XqIqBTsXV6WYZOxFWh63rXxngOMgwDo6OjXT0n1RqYWUDHUVUV6+vrKBQKVTNceq3egF035XK5hldiU1NTriejoxKe07GH9fV1O8mur697Jja3Znl14m9W8dlnn/VkBk6xWGQPwAN80wOQZdlTMyGA/auvd955p2pbtwbsTNPEzZs3USgUIIoi5ufna/4+ly9fxsjIiKtxVMoWBxs10zSRzWYdXZXmcrmqGEVRtAevex2bW7O82o2rQlVVjI6O9iQBGIbBAWAPCFiWZfU6CPK3SkN56tQpAPtTGg+WNxRFgSiKME2zpuyRyWSws7ODZDJpvwdgT2HtZWxAdSN9/vx5LC8vd6RW34m4DMOApmkwTbNjfy+nsd+6dQuPHj3q2jmpPt/0AMi7JElq2Cjquo5QKARZlpFOp+s2Zm72UtqNza1ZXu3GpaoqYrFYV0uOlWSoqirGx8e7dl5qjAmAPO/WrVtQVdUeKPcSp7F1e5bXUXEZhtGTGUNnz57F/fv38fz5c08u2e5HvhkEpv4kSRK+/vprSJKEmzdv1nwuimLV+k47Oztda9yOi62i28syHBeXrutQVRWKoiCfz+P58+ddWYr52rVrMAwDExMTnpjpRuwBkIepqop8Po9kMol4PI6trS37s8og4uGZP4ZhdGWw30lslf26uSyD079ZRSUhdeNvVhmnIe/gIDB51sFpsYZhYHx83L5ynJycxIcffghJkqCqqt2QCYLQlcbMSWyVnw96+vSpq7NfnP7NgP2ewNLSEoD9mxX76UYy6gwmAOpbXl5OwKuxeTUu6g2OAVBfqtxk5UVejc2rcVHvMAFQX/LyjURejc2rcVHvsARERORT7AEQEfkUEwARkU8xARAR+RQTABGRTzEBEBH5FBMAEZFP/V9zxXGaV7x5agAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns; sns.set()\n",
    "sns.set(font_scale=1.6)\n",
    "\n",
    "fig = plt.figure(figsize=(6,4))\n",
    "df = pd.DataFrame(CovMat/np.sqrt(np.outer(np.diag(CovMat),np.diag(CovMat))))\n",
    "df.rename(columns={0:'$b_1$',1:'$b_2$',2:'$N$',3:'$c^2_{s,0}$',4:'$c^2_{s,2}$',5:'$c^2_{s,4}$',6:'$f$',7:'$\\sigma_8$',8:'$D_A$',9:'$H$'}, inplace=True)\n",
    "df.rename(index={0:'$b_1$',1:'$b_2$',2:'$N$',3:'$c^2_{s,0}$',4:'$c^2_{s,2}$',5:'$c^2_{s,4}$',6:'$f$',7:'$\\sigma_8$',8:'$D_A$',9:'$H$'}, inplace=True)\n",
    "ax = sns.heatmap(df, cmap=\"PiYG\", vmin=-1, vmax=1)\n",
    "\n",
    "ax.set_title('$\\mathrm{EFTofLSS}, z=0.5$')\n",
    "#plt.savefig('corr-mat-eft-zeq0p5.pdf', format=\"pdf\", bbox_inches = 'tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>$b_1$</th>\n",
       "      <th>$b_2$</th>\n",
       "      <th>$N$</th>\n",
       "      <th>$c^2_{s,0}$</th>\n",
       "      <th>$c^2_{s,2}$</th>\n",
       "      <th>$c^2_{s,4}$</th>\n",
       "      <th>$f$</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>$b_1$</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.451485</td>\n",
       "      <td>-0.506048</td>\n",
       "      <td>0.343678</td>\n",
       "      <td>-0.566955</td>\n",
       "      <td>0.130634</td>\n",
       "      <td>0.210926</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$b_2$</th>\n",
       "      <td>-0.451485</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.967421</td>\n",
       "      <td>-0.986652</td>\n",
       "      <td>-0.253217</td>\n",
       "      <td>0.411791</td>\n",
       "      <td>0.391436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$N$</th>\n",
       "      <td>-0.506048</td>\n",
       "      <td>0.967421</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.965306</td>\n",
       "      <td>-0.081515</td>\n",
       "      <td>0.383110</td>\n",
       "      <td>0.286686</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$c^2_{s,0}$</th>\n",
       "      <td>0.343678</td>\n",
       "      <td>-0.986652</td>\n",
       "      <td>-0.965306</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.329496</td>\n",
       "      <td>-0.455758</td>\n",
       "      <td>-0.425693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$c^2_{s,2}$</th>\n",
       "      <td>-0.566955</td>\n",
       "      <td>-0.253217</td>\n",
       "      <td>-0.081515</td>\n",
       "      <td>0.329496</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.563878</td>\n",
       "      <td>-0.401088</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$c^2_{s,4}$</th>\n",
       "      <td>0.130634</td>\n",
       "      <td>0.411791</td>\n",
       "      <td>0.383110</td>\n",
       "      <td>-0.455758</td>\n",
       "      <td>-0.563878</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.327548</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$f$</th>\n",
       "      <td>0.210926</td>\n",
       "      <td>0.391436</td>\n",
       "      <td>0.286686</td>\n",
       "      <td>-0.425693</td>\n",
       "      <td>-0.401088</td>\n",
       "      <td>-0.327548</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                $b_1$     $b_2$       $N$  $c^2_{s,0}$  $c^2_{s,2}$  \\\n",
       "$b_1$        1.000000 -0.451485 -0.506048     0.343678    -0.566955   \n",
       "$b_2$       -0.451485  1.000000  0.967421    -0.986652    -0.253217   \n",
       "$N$         -0.506048  0.967421  1.000000    -0.965306    -0.081515   \n",
       "$c^2_{s,0}$  0.343678 -0.986652 -0.965306     1.000000     0.329496   \n",
       "$c^2_{s,2}$ -0.566955 -0.253217 -0.081515     0.329496     1.000000   \n",
       "$c^2_{s,4}$  0.130634  0.411791  0.383110    -0.455758    -0.563878   \n",
       "$f$          0.210926  0.391436  0.286686    -0.425693    -0.401088   \n",
       "\n",
       "             $c^2_{s,4}$       $f$  \n",
       "$b_1$           0.130634  0.210926  \n",
       "$b_2$           0.411791  0.391436  \n",
       "$N$             0.383110  0.286686  \n",
       "$c^2_{s,0}$    -0.455758 -0.425693  \n",
       "$c^2_{s,2}$    -0.563878 -0.401088  \n",
       "$c^2_{s,4}$     1.000000 -0.327548  \n",
       "$f$            -0.327548  1.000000  "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add Gaussian priors on parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b2 index 2\n",
      "cs0 index 4\n",
      "cs4 index 6\n",
      "cs2 index 5\n",
      "combined marged [1.650e-02 1.111e-02 1.831e-01 3.034e+02 1.572e+00 1.909e+00 1.565e+00]\n",
      "\n",
      "Correlation coefficients after adding EFT parameters priors:\n",
      "[[ 1.    -0.383 -0.318  0.078 -0.444  0.251  0.501]\n",
      " [-0.383  1.     0.753 -0.891 -0.285  0.139  0.365]\n",
      " [-0.318  0.753  1.    -0.809  0.145  0.065  0.106]\n",
      " [ 0.078 -0.891 -0.809  1.     0.452 -0.235 -0.542]\n",
      " [-0.444 -0.285  0.145  0.452  1.    -0.437 -0.699]\n",
      " [ 0.251  0.139  0.065 -0.235 -0.437  1.    -0.157]\n",
      " [ 0.501  0.365  0.106 -0.542 -0.699 -0.157  1.   ]]\n"
     ]
    }
   ],
   "source": [
    "EFTmat = Fishermat*0.\n",
    "EFT_errors = {\n",
    "    #'b': (b_fid/100)*10,\n",
    "    'b2': abs(b2_fid/100)*70,\n",
    "    #'N': (N_fid/100)*10,\n",
    "    'cs0': (cs0_fid/100)*70,\n",
    "    'cs4': (cs2_fid/100)*70,\n",
    "    'cs2': (cs2_fid/100)*70\n",
    "}\n",
    "for p in EFT_errors.keys():\n",
    "    if p not in varyzdep: continue\n",
    "    print (p, 'index', fm_index[p][0])\n",
    "    EFTmat[fm_index[p][0],fm_index[p][0]] = 1./EFT_errors[p]**2\n",
    "\n",
    "priorfactor = 1.0 #=0.0 if no priors   \n",
    "    \n",
    "Fishermat = Fishermat+priorfactor*EFTmat\n",
    "CovMat = scipy.linalg.pinv(Fishermat)\n",
    "print ('combined marged',  np.sqrt(np.diag(CovMat)))\n",
    "\n",
    "print ('\\nCorrelation coefficients after adding EFT parameters priors:')\n",
    "print (CovMat/np.sqrt(np.outer(np.diag(CovMat),np.diag(CovMat))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEaCAYAAAAYOoCaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3V1sG2d6L/A/7eyiwIpDt9iL4tTjyzTm0C3OIu6Co+CgWDsWpQO0sFqLvoqT49AGgugDiYi9WJtCFF+FWiAKAgQxk1i+KKCxToyigEXKcYBtEQ5xpIsC9YycoOhFNGrRokUjDtWLdrOac6HDOaJISsOP4Qw5/99iAHE4w3nIrN9n3s8JWZZlgYiIAueE1wEQEZE3mACIiAKKCYCIKKCYAIiIAooJgIgooJgAiIgC6jmvA6BamUwG+XwepmkimUwiHA4DACqVCgzDgKqqkGUZ9+7dO/Z4TdOg6zpmZ2eRSqVgGAYuXrwISZIQi8UQDoexubkJVVXtcyuVClRVRblcxsbGhqOYTdPExx9/jDNnzmBrawuVSsWOaWRkBGNjY0gkEk3PLxQKePr0KU6dOgUAEEURkiRB1/Wa85we1y9UVYVhGBBF0f4OoigeeY5pmsjn85Bl2T5PURTMzs5CEIQeRU4DwyLfee+996zLly83fG9ra6vuvaOOv3v3rvXee+9ZlmVZxWLR/ruqWCxazz//fN15zT6vkVdffdXa2tqytra2rNu3bx8b00HLy8vW3bt362K6fPmytby83PJx/WJra8uanJys2ffqq68ee56madbzzz9vbxcuXLA0TXMrTBpwrAH0GVEUEY/HHR+fSqWQzWYB7N89Xr161dF5yWQSpmk6uqvUNM2+c52fn3ccGwDkcjk8efKkZp8sy4jFYm0d1y9yuVzdf4toNIpCoXBsbebJkyeIRCIAwLt+6gj7APrQuXPnWjr+zJkzAIByuXxsE0OVKIool8stx9YK0zSbvpdMJls+rp/k8/m65HXu3DmoqurofEEQWPhTx5gA+oSiKPbfsiwfe7xhGPbf1bvFVu6WY7GYfZ5bBEGAYRjQdb3uvWo/RSvH9QvDMBrWrgRBcJwAiLqBCaBPFItF+28nd34HC5Jqk4IkSY6vV73DLBQKUBQFhUIBuVyuJrGoqopsNgvTNJHNZu2mpla8/vrrePXVV6EoSt2d/sF4nR7nVDabxfnz55HJZADs1zIOJlk3NavRRCKRmt+3GV3XUSgU7P8mjRIjkRPsA/ApwzDsAnVtba3l5phisdhx80i1sE2lUva+qakp3Lx5E5IkQZZlyLKMBw8eIJ1Ot3WN6nmZTAaZTAaSJGF0dLTmmq0c54SqqqhUKnj33Xft37lUKmFpaenI81577TXH14hGo01/k06a1qrJ/2A/wfj4OBYXFx037xFVMQH4lCiKdgGSTqcxNTV15PHVgqw6jLPT9mHTNLGwsFA3FPTq1au4ffs2Hj582NHnH5ROp5FOp6GqKgqFAu7evYt8Pl93DafHHScWi9U0o2UyGSwuLh77m1WH3npJFMW6gj4ejyOXy7XcAU/EJqA+MTw8XPP6cFNBNWHMz8/jyZMnHbeLHxzZc5Asy11pcmjUDCLLMubn5/Hll18CQNMmJafHNXOwoM9ms47G33dTs76VdmsGZ86cYd8BtYUJoE+Mjo7WvD7uH7yTjuKjHFfId5IETNOEpmkwTbPh5wiCgBs3bqBUKtnHOzmuVdlsFmNjYx3/Vq2qJpvDSdDJsNtcLle3LxKJuD5iiwYTm4D6xMGCwTTNY0fodDozttkdcbXQ6uSO2TAMu9BSVbVhJ64sy7h79y4AOD6uFdlsFsPDwy11IHerD6DawV4ul2v+u5bL5SPneBiGgYWFBSSTybrz2P5P7WAC6EONxpB3myzLuH37dt1+TdMgSVJHfQyrq6v2JKh8Pt+wI9cwjJrv6PQ4J6qF/8E7/1wud2yHcjf7AOLxOHRdrym4dV2va+o7SBRFzM/P1/32xWKxpcmBRFVsAuozhmEgl8u5fsdXbV453OSQy+Xw7rvv2q+PmqTViKqq+OSTT2qaQRq14S8sLGB2drbmOscdZxgGzp8/f2TzWKFQwObmJorFInK5HAqFwrEd7G5Ip9N1w041TasZuVVdu+nw9znY/6PrOjY3N9sehUXBFrIsPhPYT5wsBgcA33zzTcPjZVl21PxTXUSsuiCZLMuIRqO4efNmzR1moVCAYRgQBAGmaUKWZbvZpDryaG1tDclkEolEArIsN/0O1YXnqvFX+wJEUcTy8rJ9zUqlglQqVZMknBxnGAbGx8cxMTHRtECs3umbponp6WlomoYbN260NZy0U6qq2rWAp0+f4urVqzWJvfp9FhcXa2or1eG5Ozs7qFQqXAiO2sYEQANHUZS+XSKCqJfYBEREFFBMADRQqs9LIKLjMQHQQDFNk0MiiRxiHwARUUCxBkBEFFC+nQj2Quq/eR1CnY0//xuvQ2jquchveR1CQ7/+wx94HUJD3//Vv3gdQlP/lvhXr0No6K//4S+8DqGp1E8+6+j8Vsqbr3P/1NG1/MS3CYCIqFdCAW0LYQIgosA7cSLkdQieYAIgosALBbP8ZwIgIjrBJiAiomA6+VwwqwBMAEQUeOwEJiIKqBMB7QRgAiCiwGMNgIgooNgJTEQUUEGtAbj6tVdWVtz8eCKirjh5MuR4GyQd1QA++eQTlMvlhu9ZloVSqYQrV650cgkiItexCagNsixjeXkZw8PDDd8/+PBqIiK/CnEpiNZFo1EMDw9jZGSk4ft8UDUR9QPWANrUrPAHgHg83unHExG5LqidwG0lAMMwMD09jWfPnkEQBMRiMaTTabzwwgvdjo+IyHVudu4qioJIJAIAKJfLSCaTjs4B9h9xCgCpVMqV2NpKAIqi4M6dOwiHw9jc3ESxWMQrr7yCq1ev4q233up2jERErnKrCaha+CcSCQBAoVCAoihHJoFcLldT4JumiWw2i3Q63fX42vra586dQzQahSiKGBkZwfz8PNbX17G3t8ehn0TUd0InQo63ViiKYhf+AJBIJOy7+2YOD54RBAGVSqWl6zp1bALIZDL4sz/7M8zNzeHx48fY3d1FKBTC7u5u3bGzs7N2lYWIqF+cOOF8c8o0zYYjIXVdP7Kc1DStZyMoj/068/PzsCwLkiTh0aNH+NnPfoZbt27hwoULWFlZwfb2ds3x3333nWvBEhG5IRRyvjmlaZrd9n+QIAjQNK3pebOzsxgfH7drCoVCAbOzsy1/Jycc9QHcuXMH0WgUExMTAParKKqq4quvvkI2m0WlUoEoinbwRET9pJVHQpqm2fAOXhCEmqHvpmk2HArfKCkcJMsylpaWMD4+jlwuh3v37rk2pN5RAohGozWvRVFEMpm0OzKq1ZVqEiAi6ietPBDm/v37+PDDD+v2v/nmm5icnKzZ106TePUGe2NjAwsLC7h48SLu3bsHWZZb/qzjdGUxOBb8RNTPWnkewLVXruHy5ct1+w/fpTe7a2+2fE5VLpfD/Pw8gP0m+EQigenpaXz55ZddrwlwNVAiCrxWmoAON/U0E4vFGhb2pmkiFos1PEdV1bo7fVmWcePGDaiqWjOiqBuYAIgo8FpJAE4JgoBIJFLXFyCKYtME0qzJSJKkrscHuLwcNBFRPzgROuF4a0UymUQ+n7dfK4pSM8lL13VMTU3ZrxOJBFZXV+s+p1gs+rcPgIion7lRAwD2l3BQFAWqqtp39wdnAZfLZZRKpZpaQjqdRjabxalTpyAIAkzTxM2bN12JjwmAiALvuZPuNYYcteyDLMvY2Nio2SeKoivLPjTCBEBEgXcioOtB+zYBbPz533gdQkPn//f/8DqEhv77T4+eXOKV8ZE/9jqEhv7kVxmvQ2hq74c/8jqEhp7/3Z94HYJr3GoC8jvfJgA/8mvhT0SdaWUewCBhAiCiwGMNgIgooE6ePOl1CJ5gAiCiwGMNgIgooJgAiIgCqtUZvoOCCYCIAo81ACKigOIwUCKigHruOY4CIiIKpKDWALrS87G9vY3Hjx/j2bNn9r5KpYJSqVT30HgiIr85cSLkeBskHSeAtbU1jI+P49GjR5iamsLMzAwAIBwOQ5IkvPzyyx0HSUTkpqAmgI6bgFRVxfr6uv1a13XMzc0hnU5DEARYltXpJYiIXBXUYaAdf+vDjyqTJAnvvPMOFEVBpVJBKKBta0TUP4JaA+g4AUQiEZRKJczMzGB3d9fef/36dWiaxhoAEfnecyd/4HgbJB03AY2MjGB7extjY2MYGhqqeS8ej+Phw4edXoKIyFUnQxwG2hLDMDA9PY1nz55BEATEYjGcOXMGL7zwQs1x0Wi04yCJiNx04gQTQEsURcGdO3cQDoexubmJYrGIV155BVevXsVbb73VzRiJiFx1IqA1gLb7AM6dO4doNApRFDEyMoL5+Xmsr69jb28PKysr3YyRiMhVJ06ccLwNEkc1gEwmA13XEYvFMDw8DFmWEQqFsLu7W9fuPzs7i08//dSVYImI3PCDkz/0OgRPOEpn8/PzsCwLkiTh0aNH+NnPfoZbt27hwoULWFlZqZvt+91337kSLBGRG06ETjjeBonjPoA7d+4gGo1iYmICwH4nsKqq+Oqrr5DNZlGpVCCKIoD9WgARUb9gJ/AxDo/mEUURyWQSyWQSwH5CqO4nIuonHAbaIRb8RNSvBq1pxykuB01EgccmICKigArqKCAmACIKPDYBEREFFJuAiIgCys2lIBRFQSQSAQCUy2V75KRTmUwG8/PzboTWnUdCEhH1M7eWgqgW/olEAolEApFIBIqiOD6/UChAVdVWv45jvq0BPBf5La9DqPO319fxvzR/PuLyb/9P2esQGnrjR7/jdQgNVV74D69DaOqHiHgdQkN/8Nt/6nUIrnFrHoCiKDVL4icSCYyPjzuqBZim6UpMB/k2AfiRXwt/IurMcy6MAjJN054ge5Cu6zBNE4IgHHm+qqqQZbnrcR3EJiAiCjw31gLSNM1u+z9IEARomnbkubquu174A6wBEBHhZAujgEzTbNg8IwhCzV19s7v8RknhsHK5DEmSXG8GYgIgosBrZRTQ/fv38eGHH9btf/PNNzE5OVmzr50CvFAoIJFItHxeO5gAiCjwWmnauXbtGi5fvly3//DdfrM2/nK5+YANwzAgSZLjWDrFBEBEgffcCeedwIebepqJxWINC3vTNBGLxRqeo6pqTa1hZ2cH5XIZuVwOkiR1vV+ACYCIAi/kwjBQQRAQiUTq+gJEUWyaQA4PDzUMA2tra0ilUl2PD+AoICIinAiddLy1IplMIp/P268VRakpzHVdx9TU1JGfcVSTUaeYAIgo8E7gpOOtFdXCXlVVFAoFALV3+eVyGaVSqWFncaFQQCaTgWmayGaz0HW9g2/YWFtNQJ9++imuX7/e7ViIiDzh5lpAR836lWUZGxsbDd+rLh/hprZqAIqiYHd3t9uxEBF5wq0mIL9rqwaws7MDRVEQDocRCoUQi8Vw9uzZbsdGRNQTJ0N8IIxjs7OzmJiYsF9vbm7iwYMHME0TkUgEp0+fRjwe71qQRERu4gNhWnCw8AeAaDSKaDRqvy6VSrh06RLC4TDGxsbYX0BEvjZoTTtOdXUeQKlUgqIoKBQKEAQBly5dwqVLl7p5CSKirmMCaNPu7i4URYGiKNja2kI0GsXi4iJGRka6ER8RkevcmAjWD9pOANW7/bW1NYTDYUxMTCCZTEIUxW7GR0TkuufYCezcpUuXYBgG4vE43n//fd7tE1FfYxNQC7a2tvD6669jbGyspvOXiKgfMQG0YGJiArOzs/bwTwBN5wOUSiUOCSUiX2MfQAuq61scHv55OCFIkoSFhQV8/vnnXQiViMgdra7xMyjaSgDNOnoPJ4S1tTVsb2+3FxkRUY9wIpgLRkZGsLq66uYliIg61soDYQaJ6w+EuXPnjtuXICLqSIhNQO4Ih8NuX4KIqCMcBUREFFBMAEREARUK6MMRmQCIiBDyOgBPMAEQUeBxHgARUWCxCYiIKJBCbALyl1//4Q+8DqHO+Mgfex1CU2/86He8DqGh1P+853UIDc1v/ZfXITQ1jp97HUJD//yXX3sdQlO//fPOlqEPcSYwEVFQsQZARBRI7AQmIgosNgEREQUSO4GJiAKKM4GJiAKLNYCOlEolbG5u2q+j0SgfBUlEfYE1gDatra1hdXUVoiji1KlT9v5isQhFUZBMJpkIiMjX+DyADiwuLjZ97/Hjx924BBGRa1gDaFMoFMy2MyIaHG6OAlIUBZFIBABQLpeRTCZdOacdHScAy7IwPT2NU6dOQRAEuxloZ2cHhmG4FjgRUde4tBREtSBPJBIAgEKhYDeNd/OcdnX8rUdGRrC4uAhZlhGJRGBZFizLQiwWw/vvv8/2fyLyvVAL/2uFoih2QQ4AiUQCiqJ0/Zx2dW0U0MjISMP9z549A7D/bODTp09363JERF3jRh+AaZowDKNuv67rME0TgiB05ZxOuNrzYRgGfvGLX+Ds2bOwLAulUsnNyxERtSWEk443pzRNs9vxDxIEAZqmde2cTrg6EUwURTx8+ND+WxQ7W7KViMgNrTTtmKYJ0zTr9guCUHOH3uyOvVEB38k5nWiaALa3t7G5uQlRFHH27FkAQKVSgaZpEEWRzTlENECcJ4D79+/jww8/rNv/5ptvYnJysmZfo0RxnHbOaVfDBLC2tobbt28jHo9jc3MTkiTh/fffRzgchiRJ+OlPf2q37R8+T9M03Lx5E7/4xS9w5swZ3Lx5E0NDQ65/ESKitlnOD7127RouX75ct//wnXuz9vpyudz0s9s5pxMNE4CqqlhfX7df67qOubk5pNNpCIIAy2r8awmCgLfffhtzc3N46aWXcOXKFaysrODKlSuuBE9E1A2hJmVaI4ebepqJxWINC27TNBGLxbp2TicadgJLklT3+p133oGiKKhUKk0nfxmGgbW1NeTzeYyOjgIAm4qIyP/2WtgcEgQBkUikrklHFMWmCaSdczrRMAFEIhGUSiXMzMxgd3fX3n/9+nVomta0BjAxMQFRFLG0tIShoSEsLCx0PWAioq6zLOdbC5LJJPL5vP1aURSkUin7ta7rmJqaaumcbmqYAEZGRnD69GmMjo7Wtd/H43F89NFHNfsOJoloNIpoNAoAGB4exrlz57odMxFRV4Us51srqgW3qqooFAoAUDOjt1wuo1Qq1dzxH3dONzUdBdRs2OaDBw8wNzeHhw8f2qODhoaGGrb1cxYwEfWFFgv2VhxVeMuyjI2NjZbO6aaWJ4JFIhGsr6/bhX9VIpHgyp9E1J9cagLyu5YngpmmiXA4XLc/HA5jZ2enK0EREfXUYJXrjrVcA7Asq+mdPpeGJqJ+FLIsx9sgaTkBTExMYHl5GTMzM/jiiy/s/bu7u/j222+7GhwRUU9YLWwDpK21gD777DPkcjlMTk4iFArZ41OXlpa6GRsRUW8MWMHuVNuLwaVSKaRSKRiGgUqlYg/9JCLqOwPWtONUx6uBcoVPIup3rY7vHxSuLgdNRNQXmAD85fu/+hevQ6jzJ7/KeB1CU5UX/sPrEBqa3/ovr0NoKLPwF16H0NTLI694HUJD//6rv/c6hKZO//zlzj5gL5gZwLcJgIioV4LaBOTqIyGJiMi/WAMgIuIoICKigApm+c8EQEQUYicwEVFABbP8704n8Pb2Nh4/flzzoPhKpYJSqYTt7e1uXIKIyD17lvNtgHScANbW1jA+Po5Hjx5hamoKMzMzAPaXh5YkCS+/3OH4XCIil1mW5XgbJB03AamqivX1dfu1ruuYm5tDOp2GIAgD94MR0QBq4WHvg6TjGoAkSXWv33nnHSiKgkqlwmcEEJHvWXuW422QdJwAIpEISqUSZmZmah4Of/36dWiaxhoAEfme9Zs9x9sg6bgJaGRkBIZhYHR0FENDQzXvxeNxPHz4sNNLEBG5a8Du7J3qyjBQURSbLgsdCoXw7NkzhMNhnD59uhuXIyLqqqC2VLi6FpBhGLh16xbOnj0Ly7JQKpXcvBwRUXv2WtgGiKsTwURRxOeff27/zYfHEJEfBbUG0DQBbG9vY3NzE6Io4uzZswD2J3dpmgZRFNmcQ0QDw/p+wG7tHWrYBNTu5K61tTX88pe/xO7uLqanp+2/iYh8zbKcbwOkYQ2g3cldgiDg7bffxtzcHF566SVcuXIFKysruHLlijvRExF1waCN73eqYQ2g3cldhmFgbW0N+Xweo6OjAMCmIiLyv4B2AjdMAO1O7pqYmIAoilhaWsLQ0BAWFhbciZqIqIu4FtABx03u+uijj2r27e7u2sdFo1F7//DwMM6dO9ftmImIuotNQLVEUcTIyEjd/gcPHuCNN96oWfp5aGgIKysrdcfG4/G6BEJE5DdBXQqi5YlgkUgE6+vr9tDQqkQigcePH3ctMCKingloH0DLE8FM00Q4HK7bHw6HsbOz05WgiIh6yW+jgBRFQSQSAQCUy2Ukk8mWzs9kMpifnz/2uJZrAJZlNb3T59LPRNSXfDQPoFr4JxIJJBIJRCIRKIri+PxCoQBVVR0d23ICmJiYwPLyMmZmZvDFF1/Y+3d3d/Htt9+2+nFERJ7z0/MAFEVBIpGwXycSCccJwDTNlq7V1lpAn332GXK5HCYnJxEKhSAIAgBgaWmpnY8jIvJWCwW7aZoNC1pBEOyysF2macIwjLr9uq7DNM1jP19VVciy7Ph6bS8Gl0qlkEqlYBgGKpVKzfBPIqJ+0sronvv37+PDDz+s2//mm29icnKyozg0TbPb/g8SBAGaph1ZuOu63lLhD3RhNVCu8ElE/a6VBHDt2jVcvny5bn+nd/8Amt7lN0oKh5XLZUiS1FIzkKvLQRMR9YNW2vYjXWjqOUqr7fjAfsfvwX4Dp3ybAP4t8a9eh1Bn74c/8jqEpn6I4+8QvDCOn3sdQkMvj7zidQhNxdfqJ2D6gf7p116H4Bo3JnipqopcLnfsceFwGB988AGA5rWIcrnc9HzDMOrWb3PKtwmAiKhn9rqfAGRZbrlNPhaLNSzsTdNELBZreI6qqjW1hp2dHZTLZeRyOUiSdGQMTABEFHh7v/6N1yEA2K8BRCKRur4AURSb1g4OTxKrrsqcSqWOvZ6rzwQmIuoHfloLKJlMIp/P268VRakpzHVdx9TU1JGfcVST0UFMAEREe3vON5dVC3tVVVEoFADU3uWXy2WUSqWGncWFQgGZTAamaSKbzULX9SOvxSYgIgo8v60FdNTaP7IsY2Njo+F71eUjnGICIKLAG7Rlnp1yrQmoVCoBQM1zA4iI/Mja23O8DZKOawAPHjxApVKp2//06VNsbm5CVVV8+umnnV6GiMg1lk9GAfVaxzWAWCyGv/u7v4MgCAiHw/Y2qM/QJKLB46dRQL3UcQ0gGo1icXERpVIJp06dqnlS2MTERMsTIYiIem3Qmnac6loncDweR6VSwePHjxGNRu2Hw3CVUCLyvd8Es6Wiq53A4XAYly5dgmVZri6WRETUTewE7iJRFO1loqujgMLhME6fPu3G5YiIOjJobftOuToT2DAM3Lp1C2fPnoVlWfbQUCIiP7F+/RvH2yBxdSKYKIr4/PPP7b/58Bgi8qNBa9pxqmkC2N7exubmJkRRtEf2VCoVaJoGURTZnENEA4NNQAesra1hfHwcjx49wtTUFGZmZgDst+NLkoSXX3654Yetra3hl7/8JXZ3dzE9PW3/TUTkZ+wEPkBVVayvr9uvdV3H3Nwc0uk0BEFoOrlLEAS8/fbbmJubw0svvYQrV65gZWUFV65ccSd6IqJu4DDQ/+/w48UkScI777wDRVFQqVTsMf6HVR9EkM/nMTo6CgBsKiIi39v7r+8db4OkYQKIRCIolUqYmZmpacK5fv06NE1rWgOYmJiAKIpYWlrC0NAQFhYW3ImaiKiLuBTEASMjIzAMA6OjoxgaGqp5Lx6P46OPPqrZt7u7ax93cObv8PAwzp071+2YiYi6yvp+sAp2p5rOAxBFESMjI3X7Hzx4gDfeeKNmmeehoSGsrKzUHRuPx+sSCBGR3wS1BtDyRLBIJIL19fWaRd+A/SfRPH78uGuBERH1ivX9nuNtkLQ8Ecw0TYTD4br94XAYOzs7XQmKiKiXBq1gd6rlGoBlWU3v9JuNDiIi8rO9//yN422QtJwAJiYmsLy8jJmZGXzxxRf2/t3dXXz77bddDY6IqBeC2gfQ1lpAn332GXK5HCYnJxEKheyln5eWlroZGxFRTwS1CajtxeBSqRRSqRQMw0ClUuGDX4iobzEBtIkrfBJRvxu0ph2nXF0OmoioH7AG4DN//Q9/4XUIdZ7/3Z94HUJTf/Dbf+p1CA39819+7XUIDf37r/7e6xCa0j/1528mzb3gdQhNfZ37p47O3/vPwVrjxynfJgAiol5hDYCIKKDYB0BEFFCsARARBRQTABFRQPltiQdFURCJRAAA5XIZyWTS0TnA/nptwP5creMwARBR4PmpD6Ba+CcSCQBAoVCAoihHJoFcLldT4JumiWw2i3Q6feS1Wl4LiIho0PhpOWhFUezCH9hfar96d9+MYRg1rwVBQKVSOfZaLSeAUqlU8zAYIqJ+55cEYJpmXWEOALqu2007jWia1vC847ScAPL5PD7++GP79cFnBhMR9aNWVgM1TRPb29t121EFtFOaptlt/wcJggBN05qeNzs7i/HxcbumUCgUMDs7e+z1Wk4AZ86cwfvvv2+/PpgMiIj6USsJ4P79+7hw4ULddv/+/Y7jME3TXl35oEZJ4SBZlrG0tIRMJoOLFy9CkqSGn3NYy53AkiThj/7ojyCKImRZhmEY+Md//Ef83u/9XqsfRUTkC9/vOR8FdO3aNVy+fLluv5MC14l2ahKGYUBVVWxsbGBhYQEXL17EvXv3IMvykee1nADi8TjW19dRKpVQLBahqiouXLiAUCiEaDSKWCxmb6Io8qHwROR7v7Gct+0LguCosFdVFblc7tjjwuEwPvjgA/uzGymXy0d+Ri6Xw/z8PABgfn4eiUQC09PT+PLLL4+Mte1hoPF4HPF4HMB++5NhGNjc3MTTp0/x1Vdf4e7du9je3kY0GsXY2BiSySSTARFQaOYqAAAIE0lEQVT50l4LCcApWZaPvQM/LBaLNSzsTdNELBZreI6qqnXXkWUZN27cgKqqNSOKDut4HkB1bKooihBFESMjIzXvb25uQlVV3Lp1C2NjY7h06VKnlyQi6qo9y/I6BAD7NYBIJFLXFyCKYtM7+WZNRpIkHXu9jucBHPdAmGg0itdffx3JZPLYagwRkRf2rD3Hm9uSySTy+bz9WlGUmkleuq5jamrKfp1IJLC6ulr3OcVisft9AO2oVCp47bXXHE1NJiLqtVY6gd2WSqWgKApUVbXv7g/OAi6XyyiVSjW1hHQ6jWw2i1OnTkEQBJimiZs3bx57rZ4kgHA4jK+/9udDLoiIenFn34qjln2QZRkbGxs1+0RRPHbZh0a4FhARBV4ro4AGCRMAEQWeXzqBe40JgIgCz29NQL3CBEBEgccEQEQUUH4aBdRLTABEFHjsBCYiCig2ARERBRRHARERBRRrAEREAcUEQEQUUL/mKCAiomBiDYCIKKCC2gkcsqyAfnMiooDr+IEwRETUn5gAiIgCigmAiCigmACIiAJqIBNAJpPxOgQiIt8bqASg6zrGx8ehaZrXoTTExEREfjJQ8wAkSUIymYSu616HUkPXddy+fdvrMIiIagxUAgCAYrGIsbExr8Oo4dfE1A8KhQLK5TJ0XYcsy0gkEl6HZPNrbH6Ni/xn4BJAqVTCnTt3vA6jjh8Tk98ZhgFBEOwC7Pz585AkCaIoehyZf2Pza1zkTwPVB2AYBgBAVVUUCgVks1l7n9dKpRJkWfY6jIYMw0Amk8Hv//7v4+LFi3XvK4qC8+fP4+LFi8jlcj2LS9d1LC8v26/j8bhvalF+jc2vcTWjqipyuRwymQxUVfU6nOCxBsjy8rL14osvWuVy2bIsy9I0zbp8+bLHUVnW1taW9eKLL1r5fN7K5/PWe++9Z21tbXkdVp333nvPev75561isVj33t27d+3ftZcOXvPChQuWpmk9j6EZv8bm17gaefHFFy3LsqzJyUlfxzmoBqoGUCwWMTs7C0EQAACCINTc/ei63tM72KrqnU21PXZsbAzT09M9j+Mouq5jbGwMsizX3EFWCYJg/669VL2mruuIRqOQJKnnMTTj19j8GtdhqqraTVMffPCBb+McZAOVAA43s+i6bv9jKBQKUBQFOzs7PY/rqMRUjSuTyaBQKPQ8tipN0+zO6rW1NZim6VksjSiKgg8++MDrMBrya2x+jauqWCwiHo97HUagDUwCMAwDpmnWdHapqopYLAYASCQSnrXBN0tM1Q67ZDKJ+fl53L592/M+i0QiAUEQoCiKvU/Xdft39EIul8Ps7CwA+C4x+TU2v8YF7P9bzWazePDggf2332IMioEaBXR4pEM+n8fS0pI3wfw/RyUmXdexurpqJ4dqh53XIzYmJiagKApSqRSA/dpBMpl0/brVpBOJRFAul5FMJqGqKpLJJARBgGma0DTNk0Tu19j8GtdRRFFEOp3GJ598gnQ67fn/34NsYBLA4f8TZbNZTExM+KJdsVlikiSp5h/m5uYmbt682evw6u7wr169ik8++QSqqkKW5Z7cnb322mtIJpNIJBIwTRPnz5/Hj3/8Y7zxxhs1xz158sT1WPolNr/G5US1psvC31sDkwAAYHZ2FtlsFqdOncKZM2d6ctd6nOMSkx867A7f4YuiCEmSsLy8DFmWXe/8zeVyCIfD9th1QRDw8OFDSJKEb775xtVr92tsfo3LKV3XfXFzFnQDlQASiYQvZz06SUx+67C7ceMGpqena0ZquGVhYQGLi4s1+/xSOPg1Nr/G5dTTp0897VeifQOVAPzquMR0uMPOi+GWh1U7gzOZjKtNCNWmgE4S9+HRU926CehWbN1elqEbcVXput6zPp6DDMPA8PBwT69J9QZmFNBxVFXF6uoqSqVSzQgXrzXqsOulQqHQ9E5sYmLC9WR0VMJz2vewurpqJ9nV1VXfxObWKK9u/GZVH3/8sScjcDY3N1kD8IHA1ABkWfbVSAhg/+7rtddeq9nXqw470zRx69YtlEoliKKI2dnZut/n6tWrOHXqlKtxVJstDhZqpmkin887uistFAo1MYqiaHdeex2bW6O8Oo2rSlVVDA8Pe5IADMNgB7APhCzLsrwOgoKtWlCeOXMGwP6QxoPNG4qiQBRFmKZZ1+yRy+Wws7ODdDptvwZgD2H1MjagtpC+ePEiFhcXu9JW3424DMOApmkwTbNrv5fT2G/fvo2HDx/27JrUWGBqAORfkiQ1LRR1XUckEoEsy8hmsw0LMzdrKZ3G5tYor07jUlUViUSip02O1WSoqipGR0d7dl1qjgmAfO/27dtQVdXuKPcTp7H1epTXUXEZhuHJiKHz58/j3r17ePr0qS+XbA+iwHQCU3+SJAlffvklJEnCrVu36t4XRbFmfaednZ2eFW7HxVbV62UZjotL13WoqgpFUVAsFvH06dOeLMX8+uuvwzAMjI2N+WKkG7EGQD6mqiqKxSLS6TSSySS2trbs96qdiIdH/hiG0ZPOfiexVY/r5bIMTn+zqmpC6sVvVu2nIf9gJzD51sFhsYZhYHR01L5zHB8fx7vvvgtJkqCqql2QCYLQk8LMSWzVvw968uSJq6NfnP5mwH5NYGFhAcD+ZMV+mkhG3cEEQH3Lz8sJ+DU2v8ZF3mAfAPWl6iQrP/JrbH6Ni7zDBEB9yc8Tifwam1/jIu+wCYiIKKBYAyAiCigmACKigGICICIKKCYAIqKAYgIgIgooJgAiooD6v9hPeYlPQas3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns; sns.set()\n",
    "sns.set(font_scale=1.6)\n",
    "\n",
    "fig = plt.figure(figsize=(6,4))\n",
    "df = pd.DataFrame(CovMat/np.sqrt(np.outer(np.diag(CovMat),np.diag(CovMat))))\n",
    "df.rename(columns={0:'$b_1$',1:'$b_2$',2:'$N$',3:'$c^2_{s,0}$',4:'$c^2_{s,2}$',5:'$c^2_{s,4}$',6:'$f$',7:'$\\sigma_8$',8:'$D_A$',9:'$H$'}, inplace=True)\n",
    "df.rename(index={0:'$b_1$',1:'$b_2$',2:'$N$',3:'$c^2_{s,0}$',4:'$c^2_{s,2}$',5:'$c^2_{s,4}$',6:'$f$',7:'$\\sigma_8$',8:'$D_A$',9:'$H$'}, inplace=True)\n",
    "ax = sns.heatmap(df, cmap=\"PiYG\", vmin=-1, vmax=1)\n",
    "\n",
    "ax.set_title('$\\mathrm{EFTofLSS}, z=0.5$')\n",
    "#plt.savefig('corr-mat-eft.pdf', format=\"pdf\", bbox_inches = 'tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>$b_1$</th>\n",
       "      <th>$b_2$</th>\n",
       "      <th>$N$</th>\n",
       "      <th>$c^2_{s,0}$</th>\n",
       "      <th>$c^2_{s,2}$</th>\n",
       "      <th>$c^2_{s,4}$</th>\n",
       "      <th>$f$</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>$b_1$</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.383421</td>\n",
       "      <td>-0.317657</td>\n",
       "      <td>0.077592</td>\n",
       "      <td>-0.443676</td>\n",
       "      <td>0.251401</td>\n",
       "      <td>0.501286</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$b_2$</th>\n",
       "      <td>-0.383421</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.752971</td>\n",
       "      <td>-0.890965</td>\n",
       "      <td>-0.284757</td>\n",
       "      <td>0.138569</td>\n",
       "      <td>0.364981</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$N$</th>\n",
       "      <td>-0.317657</td>\n",
       "      <td>0.752971</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.808628</td>\n",
       "      <td>0.144903</td>\n",
       "      <td>0.064791</td>\n",
       "      <td>0.105794</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$c^2_{s,0}$</th>\n",
       "      <td>0.077592</td>\n",
       "      <td>-0.890965</td>\n",
       "      <td>-0.808628</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.451515</td>\n",
       "      <td>-0.234582</td>\n",
       "      <td>-0.541537</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$c^2_{s,2}$</th>\n",
       "      <td>-0.443676</td>\n",
       "      <td>-0.284757</td>\n",
       "      <td>0.144903</td>\n",
       "      <td>0.451515</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.437358</td>\n",
       "      <td>-0.698726</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$c^2_{s,4}$</th>\n",
       "      <td>0.251401</td>\n",
       "      <td>0.138569</td>\n",
       "      <td>0.064791</td>\n",
       "      <td>-0.234582</td>\n",
       "      <td>-0.437358</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.156943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$f$</th>\n",
       "      <td>0.501286</td>\n",
       "      <td>0.364981</td>\n",
       "      <td>0.105794</td>\n",
       "      <td>-0.541537</td>\n",
       "      <td>-0.698726</td>\n",
       "      <td>-0.156943</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                $b_1$     $b_2$       $N$  $c^2_{s,0}$  $c^2_{s,2}$  \\\n",
       "$b_1$        1.000000 -0.383421 -0.317657     0.077592    -0.443676   \n",
       "$b_2$       -0.383421  1.000000  0.752971    -0.890965    -0.284757   \n",
       "$N$         -0.317657  0.752971  1.000000    -0.808628     0.144903   \n",
       "$c^2_{s,0}$  0.077592 -0.890965 -0.808628     1.000000     0.451515   \n",
       "$c^2_{s,2}$ -0.443676 -0.284757  0.144903     0.451515     1.000000   \n",
       "$c^2_{s,4}$  0.251401  0.138569  0.064791    -0.234582    -0.437358   \n",
       "$f$          0.501286  0.364981  0.105794    -0.541537    -0.698726   \n",
       "\n",
       "             $c^2_{s,4}$       $f$  \n",
       "$b_1$           0.251401  0.501286  \n",
       "$b_2$           0.138569  0.364981  \n",
       "$N$             0.064791  0.105794  \n",
       "$c^2_{s,0}$    -0.234582 -0.541537  \n",
       "$c^2_{s,2}$    -0.437358 -0.698726  \n",
       "$c^2_{s,4}$     1.000000 -0.156943  \n",
       "$f$            -0.156943  1.000000  "
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "varying zdep = \n",
      "\t ['f', 'b', 'b2', 'N', 'cs0', 'cs2', 'cs4']\n",
      "at z = 0.5 : zdep fiducials = \n",
      "\t [ 7.042e-01  1.468e+00 -2.756e-01  8.161e+02  7.650e+00  2.978e+00  2.440e+01]\n",
      "at z = 0.5 : marged err = \n",
      "\t [1.650e-02 1.111e-02 1.831e-01 3.034e+02 1.572e+00 1.909e+00 1.565e+00]\n",
      "at z = 0.5 : unmarged err = \n",
      "\t [2.771e-03 1.092e-03 5.722e-03 7.241e+00 6.359e-02 2.498e-01 1.684e-01]\n",
      "at z = 0.5 : marged % err = \n",
      "\t [ 2.343  0.757 66.428 37.172 20.546 64.098  6.415]\n",
      "at z = 0.5 : unmarged % err = \n",
      "\t [0.394 0.074 2.076 0.887 0.831 8.387 0.69 ]\n"
     ]
    }
   ],
   "source": [
    "# Print the errors nicely\n",
    "Nvar  = len(varypars)\n",
    "\n",
    "#varyzdep = list(set(varypars)&set(zdep))\n",
    "varyzdep = varypars\n",
    "\n",
    "zdepfid = np.array([fids[p] for p in varyzdep])\n",
    "print ('\\nvarying zdep = \\n\\t', varyzdep)\n",
    "if len(varyzdep)>0:\n",
    "    for iz, z in enumerate(z_vals):\n",
    "        zdepind = np.array([fm_index[p] for p in varyzdep])\n",
    "        print ('at z =',z,': zdep fiducials = \\n\\t', zdepfid[:,iz])\n",
    "        merrs = np.sqrt(np.diag(CovMat))[zdepind[:,iz]]\n",
    "        print ('at z =',z,': marged err = \\n\\t', merrs)\n",
    "        uerrs = np.sqrt(1./np.diag(Fishermat))[zdepind[:,iz]]\n",
    "        print ('at z =',z,': unmarged err = \\n\\t', uerrs)\n",
    "        print ('at z =',z,': marged % err = \\n\\t', merrs/np.abs(zdepfid[:,iz])*100.)\n",
    "        print ('at z =',z,': unmarged % err = \\n\\t', uerrs/np.abs(zdepfid[:,iz])*100.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##### Alkistis & Dida, 2019"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}