GC-elFish-TNS-zeq0p5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TNS 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'"
   ]
  },
  {
   "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', 'sv']\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",
    "varypars = ['f','b','b2','N','sv'] \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",
    "\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",
    "\n",
    "#Mean MCMC values for TNS at z=0.5 (to compare with the MCMC)\n",
    "sv_fid = 9.544  \n",
    "b_fid = 1.498 \n",
    "b2_fid = 0.4233 \n",
    "N_fid = -468.4 \n",
    "fg_fid = 0.777\n",
    "\n",
    "#kmax at z=0.5 for multipoles: 0,2,4 \n",
    "kmax = [0.227, 0.227, 0.129]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ps column names: ('k', 'pl', 'pdd', 'pdt', 'ptt', 'P5')\n",
      "abc column names: ('k', 'bu2f2', 'bu2f3', 'bu4f2', 'bu4f3', 'bu6f3', 'b2u2f1', 'b2u2f2', 'b2u2f3', 'b2u2f4', 'b2u4f2', 'b2u4f3', 'b2u4f4', 'b2u6f3', 'b2u6f4', 'b2u8f4', 'u2f4', 'u4f3', 'u4f4', 'u6f3', 'u6f4', 'u8f4')\n",
      "bias column names: ('k', 'ddb0b2', 'ddb02', 'ddb0', 'ddb22', 'ddb2', 'ddnone', 'dtb2', 'dtb0', 'dtnone')\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', 'P5' #P5 not used \n",
    "    abccolnames = 'k', 'bu2f2', 'bu2f3', 'bu4f2', 'bu4f3', 'bu6f3', 'b2u2f1', 'b2u2f2', 'b2u2f3',\\\n",
    "    'b2u2f4', 'b2u4f2', 'b2u4f3', 'b2u4f4', 'b2u6f3', 'b2u6f4', 'b2u8f4', 'u2f4', 'u4f3', 'u4f4',\\\n",
    "    'u6f3', 'u6f4', 'u8f4' \n",
    "    biascolnames = 'k', 'ddb0b2', 'ddb02', 'ddb0', 'ddb22', 'ddb2', 'ddnone', 'dtb2', 'dtb0', 'dtnone'\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": "markdown",
   "metadata": {},
   "source": [
    "## Set up the survey properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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": 8,
   "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": 9,
   "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",
    "\n",
    "print (fg_fid_i(zbin), Dg_fid_i(zbin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "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": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "900 300\n",
      "k are preserved: True\n"
     ]
    }
   ],
   "source": [
    "#construct the P_TNS in 2d (i.e. in k-mu space)\n",
    "def Ptns_2d(mu_array,f,\n",
    "            b=b_fid,b2=b2_fid,N=N_fid,\n",
    "            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",
    "    FoG = 1. / ( 1. + np.outer(k_array**2,mu_array**2)*sv**2/2. )\\\n",
    "    \n",
    "    res = FoG\\\n",
    "    *((b**2*(D**2*np.outer(ps['pl'],muones) + D**4*np.outer(ps['pdd'],muones))\\\n",
    "    + D**4*(b*b2*np.outer(bias['ddb0b2'],muones) + b**2*np.outer(bias['ddb02'],muones)\\\n",
    "    + b*np.outer(bias['ddb0'],muones) + b2**2*np.outer(bias['ddb22'],muones)\\\n",
    "    + b2*np.outer(bias['ddb2'],muones) + np.outer(bias['ddnone'],muones)) + N\\\n",
    "    + 2*(b*(D**2*f*np.outer(ps['pl'],mu_array**2) + D**4*f*np.outer(ps['pdt'],mu_array**2))\\\n",
    "    - D**4*f*(b2*np.outer(bias['dtb2'],mu_array**2) + b*np.outer(bias['dtb0'],mu_array**2) + np.outer(bias['dtnone'],mu_array**2)))\\\n",
    "    + (D**2*f**2*np.outer(ps['pl'],mu_array**4) + D**4*f**2*np.outer(ps['ptt'],mu_array**4)))\\\n",
    "    + D**4*( b*f**2*np.outer(abc['bu2f2'],mu_array**2) + b*f**3*np.outer(abc['bu2f3'],mu_array**2)\\\n",
    "    + b*f**2*np.outer(abc['bu4f2'],mu_array**4)+b*f**3*np.outer(abc['bu4f3'],mu_array**4)\\\n",
    "    + b*f**3*np.outer(abc['bu6f3'],mu_array**6)+b**2*f*np.outer(abc['b2u2f1'],mu_array**2)\\\n",
    "    + b**2*f**2*np.outer(abc['b2u2f2'],mu_array**2)+b**2*f**3*np.outer(abc['b2u2f3'],mu_array**2)\\\n",
    "    + b**2*f**4*np.outer(abc['b2u2f4'],mu_array**2) + b**2*f**2*np.outer(abc['b2u4f2'],mu_array**4)\\\n",
    "    + b**2*f**3*np.outer(abc['b2u4f3'],mu_array**4) + b**2*f**4*np.outer(abc['b2u4f4'],mu_array**4)\\\n",
    "    + b**2*f**3*np.outer(abc['b2u6f3'],mu_array**6) + b**2*f**4*np.outer(abc['b2u6f4'],mu_array**6)\\\n",
    "    + b**2*f**4*np.outer(abc['b2u8f4'],mu_array**8) + f**4*np.outer(abc['u2f4'],mu_array**2)\\\n",
    "    + f**3*np.outer(abc['u4f3'],mu_array**4) + f**4*np.outer(abc['u4f4'],mu_array**4)\\\n",
    "    + f**3*np.outer(abc['u6f3'],mu_array**6) + f**4*np.outer(abc['u6f4'],mu_array**6)\\\n",
    "    + f**4*np.outer(abc['u8f4'],mu_array**8)))\n",
    "                \n",
    "    return res\n",
    "\n",
    "# Now convert to multipoles\n",
    "def Ptns(f,b=b_fid,b2=b2_fid,N=N_fid,\n",
    "         sv=sv_fid,sig8_z=sigma80_fid*Dg_fid_i(zbin),\n",
    "         ps=ps_fid,abc=abc_fid,bias=bias_fid):\n",
    "    \n",
    "    pk2d = Ptns_2d(mu,f,b,b2,N,sv,sig8_z,ps,abc,bias)\n",
    "    \n",
    "    return elgc.pkto1d(ps['k'],mu,pk2d,MLPS)\n",
    "\n",
    "k_test, plk_test = Ptns(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": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(900, 900)\n",
      "[[1083.117    0.       0.    ...    0.       0.       0.   ]\n",
      " [   0.      26.35     0.    ...    0.       0.       0.   ]\n",
      " [   0.       0.       9.061 ...    0.       0.       0.   ]\n",
      " ...\n",
      " [   0.       0.       0.    ...   99.981    0.       0.   ]\n",
      " [   0.       0.       0.    ...    0.     100.746    0.   ]\n",
      " [   0.       0.       0.    ...    0.       0.     101.52 ]]\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": 14,
   "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 sv \n",
    "def dPtns_dsv():\n",
    "    \n",
    "    P_svp = Ptns(fg_fid,b_fid,b2_fid,N_fid,sv_fid*(1+params_step),\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_svm = Ptns(fg_fid,b_fid,b2_fid,N_fid,sv_fid*(1-params_step),\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1] \n",
    "    \n",
    "    deriv_sv = (P_svp - P_svm)/(2*sv_fid*params_step)\n",
    "\n",
    "    return deriv_sv\n",
    "\n",
    "#Derivative wrt b\n",
    "def dPtns_db():\n",
    "    \n",
    "    P_bp = Ptns(fg_fid,b_fid*(1+params_step),b2_fid,N_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_bm = Ptns(fg_fid,b_fid*(1-params_step),b2_fid,N_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 dPtns_db2():#f,D,b,b2,N,sv,ps=ps_fid,abc=abc_fid,bias=bias_fid):\n",
    "    \n",
    "    P_b2p = Ptns(fg_fid,b_fid,b2_fid*(1+params_step),N_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1]\n",
    "    P_b2m = Ptns(fg_fid,b_fid,b2_fid*(1-params_step),N_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 dPtns_dN():\n",
    "    \n",
    "    P_Np = Ptns(fg_fid,b_fid,b2_fid,N_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_Nm = Ptns(fg_fid,b_fid,b2_fid,N_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_N = (P_Np - P_Nm)/(2*N_fid*params_step)\n",
    "            \n",
    "    return deriv_N\n",
    "\n",
    "#Derivative wrt f\n",
    "def dPtns_df():\n",
    "    \n",
    "    P_fp = Ptns(fg_fid*(1+params_step),b_fid,b2_fid,N_fid,sv_fid,\n",
    "                 sigma80_fid*Dg_fid_i(zbin),ps=ps_fid,abc=abc_fid,bias=bias_fid)[1] \n",
    "    P_fm = Ptns(fg_fid*(1-params_step),b_fid,b2_fid,N_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": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.227, 0.227, 0.129]\n",
      "900 -> 422 or 422\n"
     ]
    }
   ],
   "source": [
    "#for i in range(1,299):\n",
    "#    print(ps_fid['k'][i+1]-ps_fid['k'][i])\n",
    "\n",
    "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": 16,
   "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",
    "    'sv' : lambda z: dPtns_dsv(),\n",
    "    'b'  : lambda z: dPtns_db(),\n",
    "    'b2' : lambda z: dPtns_db2(),\n",
    "    'N'  : lambda z: dPtns_dN(),\n",
    "    'f'  : lambda z: dPtns_df()\n",
    "}\n",
    "            \n",
    "# Which should depend on redshift\n",
    "zdep = 'sv','b','b2','N','f'\n",
    "\n",
    "# Which ones are independent of redshift and should appear in the FM once (sum over z contributions)\n",
    "# We do not consider shape parameters or multiple z-bins in this work, but this would be useful for future extentions\n",
    "zidep = [i for i in derivs.keys() if i not in zdep]\n",
    "print (zidep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%time\n",
    "\n",
    "# Save the indexing for posterity\n",
    "fm_index = {'sv':[],'b' :[],'b2':[],'N' :[],'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": 19,
   "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": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a fiducial value dictionary\n",
    "fids = {\n",
    "    'sv' : sv_fid*np.ones(Nzbin),   \n",
    "    'b'  : b_fid*np.ones(Nzbin),   \n",
    "    'b2' : b2_fid*np.ones(Nzbin),  \n",
    "    'N'  : N_fid*np.ones(Nzbin),   \n",
    "    'f'  : fg_fid*np.ones(Nzbin) #fg_fid_i(z_vals)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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', 'sv_0.5']\n",
      "Fisher matrix:\n",
      " [[ 1.570e+05  1.858e+05  1.423e+04  1.747e+01 -1.252e+04]\n",
      " [ 1.858e+05  5.221e+05  5.014e+04  6.561e+01 -1.665e+04]\n",
      " [ 1.423e+04  5.014e+04  5.484e+03  7.171e+00 -1.429e+03]\n",
      " [ 1.747e+01  6.561e+01  7.171e+00  9.406e-03 -1.775e+00]\n",
      " [-1.252e+04 -1.665e+04 -1.429e+03 -1.775e+00  1.033e+03]]\n",
      "\n",
      "Correlation matrix for parameters:\n",
      "\n",
      "[[ 1.    -0.619  0.211 -0.06   0.974]\n",
      " [-0.619  1.     0.324 -0.48  -0.465]\n",
      " [ 0.211  0.324  1.    -0.984  0.385]\n",
      " [-0.06  -0.48  -0.984  1.    -0.246]\n",
      " [ 0.974 -0.465  0.385 -0.246  1.   ]]\n"
     ]
    }
   ],
   "source": [
    "print ('Order of columns/rows:')\n",
    "print (fm_head)\n",
    "print ('Fisher matrix:\\n',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 ('\\nCovariance matrix diagonal:\\n', )\n",
    "#print( np.diag(CovMat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "varying zdep = \n",
      "\t ['f', 'b', 'b2', 'N', 'sv']\n",
      "at z = 0.5 : zdep fiducials = \n",
      "\t [ 7.770e-01  1.498e+00  4.233e-01 -4.684e+02  9.544e+00]\n",
      "at z = 0.5 : marged err = \n",
      "\t [2.923e-02 8.063e-03 4.339e-01 3.387e+02 3.750e-01]\n",
      "at z = 0.5 : unmarged err = \n",
      "\t [2.524e-03 1.384e-03 1.350e-02 1.031e+01 3.111e-02]\n",
      "at z = 0.5 : marged % err = \n",
      "\t [  3.761   0.538 102.494  72.309   3.929]\n",
      "at z = 0.5 : unmarged % err = \n",
      "\t [0.325 0.092 3.19  2.201 0.326]\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": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the covariance matrix and the fiducial values into a binary file\n",
    "#np.savez(op.join(PATHTOINPUTFILES,'covmats/covmat_'+runname),covmat=CovMat,fmeans=zdepfid[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEUCAYAAAAlXv26AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAH9JJREFUeJzt3V+IG2eaLvBHbc/MwolKvjiwh52UWRY2iVUyB5Z4iMp7Ffd0V5q96Q7T5Ss7xpF9sXZ3NqNmWWjLjGJYNmrOxsZzICl7Op4rV8wYFibu6owNe2GVwH3ZKiUDB4ZxtXf3Ll3VPTDnTGydCyMhtf6VWqWW9On5DYJUSVX1yh4/9emtT6VIuVwug4iIhDUx6AKIiKi/GPRERIJj0BMRCY5BT0QkOAY9EZHgGPRERII7POgCqJHrupicnISiKEgkEohGoyiVSrBtG7quIxqNYmdnB7Ztw/M8bGxsAAByuRwKhQJc18WjR48gSVJ1n77vY3l5Gevr61BVFdlsFrIsAwAsy8Lm5iaOHDkCAJBlGYqiwHEcaJrWUF8mk8Ha2hp838f09DRmZmaavm7c2LYN13Uhy3L1z67yZ9yK7/tYW1uDqqrV7UzTRDqdrvv7I+pJmYZOPp8vf/zxxw3rXnvttYbXzs7O1i3fvXu3/N5775UvX77cdN9XrlxpeP1nn33WcKzZ2dny3bt3W9b48ccfNxx7nD19+rThz/y9997ruF2xWCy/9tpr1cepU6fKxWKxX2XSmOKIfgj5vo/Tp08Heq2u6/B9v270l81mMTk5CcdxoChK3ev3LhuGgYcPH9atU1UViURin9WPJ8MwGv7O4vE4LMvq+Gnn4cOHiMViAMBRPPUFe/RDyPO8jh/5K2RZhud5Det0XceVK1fabuv7fsvndF0PdHx6aW1treHkePz4cdi2HWh7SZIY8tQ3DPoh1M1oOpFIVEeDtdLpNBzHgWVZLbeVJAmu68JxnIbnKtcHqDPXdRs+VQEv/3yDBj1RPzHoh9De9ko7rUaCkiQhnU53HNW///77eO+992CaZsMIv5s62rEsC6ZpwrIsGIYB13WrzzmOg7m5OczNzcF1XViWhXPnzrX9tNFMLpfDiRMnkMlkALz8tGKaZij1d9Kq1lgsVvdeW6mckCt/Ps1OvES9YI9eYKlUCqZpwjAMpFKppq9ZWloC8HImTSaTgaIoeOedd1q+vluVE0jt/hYWFnDx4kUoigJFUfDRRx9hcXERrutC07TqbKKgrQzbtrGzs4OPPvoIrutWZx99/vnnbbc7d+5c4PcRj8erf1Z77W2ddaPyHmv7+HNzc7h+/Xrg9h1RR4O+GkzBtJp1s9femTKV7TzPa/r83tdeuXKl/Oabb3acURNk1o3neeU333yz6XFqt/U8r67Gbu3d7sqVK+WnT5/ua1/70ervpjKjplsff/xxw+wool6wdSM4VVWhqipWVlYCvTabzeLRo0cAXrZD9sv3fRSLxaajUlVVG9oTvVyMrN0ul8sFmr8epmbXSID9j/SPHj3K3j6FikE/BrLZLEzTbAhX3/eb9oMlScKFCxdQKBT2dbxKyHfqNdc+3yosu5HL5TAzMwNVVXveVzcqJ5W9vfpmF2j3MgyjYV0sFuupHUS0F3v0Y6B2umXttEnP82DbdtOLrqqq4rPPPtvX8VzXRSwWazmqrgRimKPuXC6HkydPdnUBOawefeXTyN7rCp7nIZlMttyn67pYWVmBrusN27E/T2Fi0I+JdDqNEydOIJ/P14X92tpa0wuvruvue3rlgwcPcPr0aaiq2nTWT7FYhKIooc0br4R87Ui+3QXoitXV1VCODwDJZBKO49QFtOM4OHnyZMttZFlGNptt+HPI5/NtTxBE3WLrRjCt2iWV6ZY7Ozt1633fb9qLX1lZQTqd7vr4tm3j1q1bkGW52gLa254wDAMfffRRdbldm8J1XZw4caJlz9qyLJRKJeTzeRiGAcuysLCw0HXdvVpaWmqYzlksFutOqpV7GO19L3unm5ZKpZafHoj2gyP6IVe5yVUlHM6dO4d4PI6LFy82jAQzmQxM04Trukin0w1tjFQqhe3t7epyLBar3tysNux3dnaajjRrj1O5qVkmk0E0GgWA6o3X9h6zMj9ckiT4vl9Xm+M4WFlZgeu6yGQy0HW9afsln8837b27rovV1VX4vo/FxUUUi0VcuHAhtOmhQcmyjFQqBcMwIMsyNjc3cf369YbX7T2p6bpenYK6vb2NnZ0d3L9//6DKpjERKZf54+A0/EzT5G0ZiPaJrRsiIsEx6Gno2bZ94FMmiUTCoKeh5/s+pxsS9YA9eiIiwQ181s0bqb8YdAmh+/Wf/eugS+iLH/zzXw26hND9+fNw7tA5bP7mn/560CX0RfF/P+t5H91kzjfGf/R8vGHA1g0RkeAGPqInIjpIhw5HBl3CgWPQE9FYiYxfzjPoiWi8TIxhw5pBT0RjJTIxfkN6Bj0RjRWO6ImIBBdh0BMRie3QIbZuiIiExtYNEZHg2LohIhLcBGfdEBGJjV+YIiISHHv0RESC471uiIgEx4uxRESCmxjDJv2+g/7WrVuIRCJQVRWyLOOVV14Jsy4ior7o14jeNE3EYjEAgOd50HU90DbAy5/LBIBUKtWX2vYd9IqiIJlMAgBKpRI8z8Px48cZ+EQ01PpxMbYS8pqmAQAsy4Jpmm3D3jCMumD3fR+5XA5LS0uh17fvt1wJeQCIx+NIJpOwbTuUooiI+uXQ4UjgR1CmaVZDHgA0TauO1ltxXbduWZIk7OzsdPdmAtr3iH5lZQVff/014vE4VFXF8ePHsb29HWZtREShC7tF7/t+Q2gDgOM48H0fkiQ13a5YLMJ1XciyHG5BTex7RH/8+HHcvn0bqqrCcRwsLy8fSMFERL2YmIgEfvi+j62trYZHpacOvAzsSm++liRJKBaLLetIp9OYm5urjvwty0I6nQ7/DaOHEb2qqvjqq68wNTVV18YhIhpm3VyMvXPnDm7evNmw/tKlS7h8+TIAtBy1Nwv/Wqqq4vPPP8fc3BwMw8Dq6mrL0X+vAgX97du3cf78+bp10WgUU1NTfSmKiKhfurnXzdkzZzE7O9uwfm8g147wg3JdF7ZtY2NjAysrK5icnMTq6ipUVe16X50ECvqVlZWGoCciGkXdzKOXJKnjKLvV857ntd3OMAxks1kAQDabhaZpWFxcxKNHj0If2QcK+nK5jGfPnsGyLEiShEQigWPHjoVaCBHRQTh8ONz5lYlEommo+76PRCLRdBvbthtG7qqq4sKFC7Btu24GTxgCv+O5uTn8/ve/x+PHj3H27FmcP38eW1tboRZDRNRvE5GJwI8gJElCLBZraN/IstxyZN6q1aMoSl/69IGDfnV1FdlsFtevX8eTJ09w/vx5XL16Fbu7u6EXRUTUL93MuglK13Wsra1Vl03TrPsylOM4WFhYqC5rmoYHDx407Cefzw+uRy9JEuLxeN26ytz5Tz/9FD/96U9DL4yIqB/68cMjqVQKpmnCtu3qaL32W7Ge56FQKNTN0FlaWkIul8ORI0cgSRJ838fFixdDrw0IGPTz8/MtZ95w7jwRjZJ+3dSs3e0OVFXFxsZG3TpZlvtyu4NmArVu0uk08vk8rl69imfPntU916+v7BIR9cPExETghygCf2HqF7/4BQzDwKlTp3D06FHIsozt7e2+fZOLiKgfDh8SJ8CD6uqbsalUCqlUqnq3ykQigWg02q/aiIhCxx8HD2jvhVkiolEhUksmKP7CFBGNlQh/YYqISGxs3RARCe7woUODLuHAMeiJaKxwRE9EJLh+fWFqmDHoiWiscERPRCQ4Tq8kIhIcWzdERIIL+4dHRsHAg/7Xf/avgy4hdH/3x38YdAl98cs//MugSwjdf+Ib/OWvXxt0GaH797f+bdAlDK2gPygikoEHPdEgiRjy1B4vxhIRCW4iwi9MEREJjbNuiIgE971D3x90CQeOQU9EY4WtGyIiwbF1Q0QkuEMc0RMRiW1igkFPRCQ0fmGKiEhwnHVDRCQ4juiJiATHHj0RkeA4j56ISHCHOKInIhIbe/RERII7zFk3RERi69ctEEzTRCwWAwB4ngdd17vaPpPJIJvN9qM0jN9nGCIaa4cihwI/gqqEvKZp0DQNsVgMpmkG3t6yLNi2vZ+3E0hoQX/v3r2wdkVE1DcTkUOBH0GZpglN06rLmqYFDnrf97t+D90K3Lq5desWPM9r+ly5XEahUMBPfvKT0AojIuqHsOfR+74P13Ub1juOA9/3IUlS2+1t24aqqqHWtFfgoFdVFXfv3sXJkyebPt/sjRIRDZvvHfpeqPsrFovV3nwtSZJQLBbbhrjjOH0PeaCLoI/H4zh58iSmp6ebPt/prEVENAy6acn4vt+0tSJJUjXzWo3am4X/Xp7nQVGUvrdvupp10yrkASCZTPZcDBFRv3Uzj/7OnTu4efNmw/pLly7h8uXL1eX9BLVlWXV9/X7qGPSu62JxcRFff/01JElCIpHA0tIS3njjjYOoj4goVN2M6M+ePYvZ2dmG9bUj+FbdjFbXNIGXuaooSuA6etUx6E3TxLVr1xCNRlEqlZDP53HmzBmcPn0aH3744UHUSEQUmm6CvrZF00oikWga6r7vI5FINN3Gtu26TwHb29vwPA+GYUBRlND79h2D/vjx44jH4wAAWZYxPT2NbDaLlZUV3Lt3jzNtiGikREK+qZkkSYjFYg29elmWW54k9n6ZynVdrK+vI5VKhVpbRV2zKpPJ4N1338XVq1fx1VdfYXd3F5FIBLu7uw0bptPpA5n/SUQUpsOR7wd+BKXrOtbW1qrLpmnWhbbjOFhYWGi7j3atnl7VBX02m0W5XIaiKPjyyy/x9ttvY3l5GadOncK9e/ewtbVVt/G3337bt8KIiPqhH1+YqoS6bduwLAtA/ajd8zwUCoWmg2PLspDJZOD7PnK5HBzH6fEdNoqUy+Vy7YpSqVRt1QAvP1LYtg3btlEoFLCzswNZlgG8HNVPTU31VMD/uRz8a8Kj4u/++A+DLqEvfnntXwZdQuj+8tevDbqEvpj4npi34v3vZ97seR+F/7od+LXJ/3G+5+MNg4YefW3IAy/7TLquV89OlS9GVcKeiGiUhN2jHwVd372SAU9Eo2wCDHoiIqHxh0eIiAR3eII/PEJEJDT26ImIBMcePRGR4LqZHy8KBj0RjRUGPRGR4A7xYiwRkdjYoyciEhzn0RMRCY49eiIiwTHoiYgEF2GP/uD94J//atAlhO6XfxDvdr4AcGb5HwddQuj+fvVvB11CXxz73YlBl9AXk+j9NsW8BQIRkeAi4MVYIiLBRQZdwIFj0BPRWOGInohIcBGO6ImIRMcRPRGR0DiPnohIeGzdEBEJjRdjiYgEx4uxRETC44ieiEhovB89EZHw2LohIhIaL8YSEQmOF2OJiITHET0RkdDYuglga2sLpVIJsizj2LFjAICdnR0Ui0XIsoxXX3019CKJiMIS6dMtEEzTRCwWAwB4ngdd1/uyzX50dWpbX1/H3NwcvvzySywsLOCDDz4AAESjUSiKgh//+Md9KZKIKCyRLv4XVCWwNU2DpmmIxWIwTTP0bfarq6C3bRtPnjzB9evX8Zvf/AapVApXr17F7u4uJElCuVzuS5FERGGJYCLwIyjTNKFpWnVZ07RAQd/tNvvVVdAritKw/LOf/QymaWJnZweRyPhdzSaiURPp4tGZ7/twXbdhveM48H0/tG160VXQx2IxFAoFfPDBB9jd3a2uP3/+PIrFIkf0RDT0uhnR+76Pra2thkdtGBeLxWqfvZYkSSgWi01r2M82vejqYuz09DS2trYwMzODV155pe65ZDKJ+/fvh1ocEVHYumnJ3LlzBzdv3mxYf+nSJVy+fBnAy9G5JEkNr2kW5BX72aYXgYLedV0sLi7i66+/hiRJSCQSOHr0KN54442618Xj8b4USUQUmi4aD2fPnsXs7GzD+r0hvZ92Sz9aNK0ECnrTNHHt2jVEo1GUSiXk83mcOXMGp0+fxocfftjvGomIQhPposUsSVLTkffe1zTjeV6o2/Qi0GeY48ePIx6PQ5ZlTE9PI5vN4smTJ3jx4gXu3bvXl8KIiPqi3MUjgEQi0TSgfd9HIpEIbZteNAR9JpPBu+++i6tXr+Krr77C7u4uIpFI3cXXinQ6faAfP4iIehZy0EuShFgs1pCFsiy3HLnvZ5teNAR9NptFuVyGoij48ssv8fbbb2N5eRmnTp3CvXv3sLW1Vff6b7/9NvSiiIj6plwO/ghI13Wsra1Vl03TRCqVqi47joOFhYWutglT0x79tWvXEI/HMT8/D+DlxVjbtvH48WPkcjns7OxAlmUAL0f1RESjIvIi/H2mUimYpgnbtquj9NrbGXieh0KhUDfbptM2YYqU9zH5vTLRvxL2vXB3N3rex7D5zz98M+gS+uLM8j8OuoTQ/f3q3w66hL449rsTgy6hLyblpZ738f+2/xD4td8/8t96Pt4w2NfdK8MIeCKigRjD73WO3/06iYjGDO9HT0RjpZt59KJg0BPReBm/nGfQE9GYeTF+Sc+gJ6KxEhm/nGfQE9GYYdATEQmOQU9EJDjOuiEiEht79EREouOInohIcOOX8wx6IhovbN0MwJ8/VwZdQuh+8Os+3Ad1CIh6p8efn3s86BJCt7n780GXMLzYuiEaLyKGPHUwfjnPoCei8VJ+LuYn7nYY9EQ0XsYv5xn0RDReyrypGRGR4HgxlohIbBzRExGJjkFPRCS28nMGPRGR0Mrs0RMRCY6tGyIisfFiLBGR6Ni6ISISW/k7Bj0RkdB4MZaISHRj2KOf6PSC27dvH0QdREQH40U5+EMQHYPeNE3s7u4eRC1ERH1XLpcDP0TRsXWzvb0N0zQRjUYRiUSQSCRw7Nixg6iNiCh8vE1xo3Q6jfn5+epyqVTCF198Ad/3EYvF8OqrryKZTPa1SCKisAzTD4+YpolYLAYA8DwPuq53tX0mk0E2m+34uo5BXxvyABCPxxGPx6vLhUIBU1NTiEajmJmZwfnz57sqlIjoIA3LF6YqIa9pGgDAsiyYphk47C3Lgm3bgV6771k3hUIBpmnCsixIkoSpqSlMTU3td3dERAeji6D3fR++7zeslyQJkiT1VIZpmrh//351WdM0zM3NBQr6ZjW101XQ7+7uwjRNmKaJp0+fIh6P4/r165ienu7qoEREg9LNRdY7d+7g5s2bDesvXbqEy5cv77sG3/fhum7Desdx4Pt+x5OIbdtQVTXw8QIFfWX0vr6+jmg0ivn5eei6DlmWAx+IiGgodDGiP3v2LGZnZxvW9zqaLxaL1d783v0Wi8W2Ie44TlchDwQI+qmpKbiui2QyiU8++YSjdyIaad1cjI2F0KJpptWovVn47+V5HhRF6ap90zHonz59ivfffx8zMzN1F2GJiEbRiz89H3QJALrvswMvL8BWLt52I9Csm3Q6XZ1WCaDlfPpCocCplkQ03F6EO73Stm0YhtHxddFoFDdu3ADQuvXjeV7L7V3XhaIo+6qxY9CnUikAjdMq9wa/oihYWVnBr371q30VQkR0EML+KUFVVbvumScSiaah7vs+EolE021s2677FLC9vQ3P82AYBhRFaVtDx6BvdcF1b/Cvr69ja2ur0+6IiAaqHPKIfj8kSUIsFmvo1cuy3HK0v3fapeu6WF9frw7G2+l4r5ugpqen8dZbb4W1OyKivig/fxH40U+6rmNtba26bJpmXWg7joOFhYW2+2jX6qkVWtADwLVr18LcHRFR6MrfPQ/86KdKqNu2DcuyANSP2j3PQ6FQaHrR1rIsZDIZ+L6PXC4Hx3HaHivU+9FHo9Ewd0dEFLphutdNu2/BqqqKjY2Nps9pmtbV7Bv+8AgRjZVh6NEfNAY9EY2XIRrRHxQGPRGNlWG5e+VBYtAT0VgZph79QWHQE9FYKQ/JLRAOEoOeiMYKL8YSEQmOrRsiIsFxRE9EJLqQb2o2Chj0RDRWXvzpu0GXcOAY9EQ0Vti6ISISHVs3B+9v/umvB11C6P79rX8bdAl9cex3JwZdQug2d38+6BL64vj/+p+DLqEvvjH+o+d9cERPRCS48ncMeiIioXEePRGR4F78X866ISISGls3RESCY9ATEQmuzOmVRERi44ieiEhwnHVDRCS4F3/krBsiIqFxRE9EJDj26ImIBMegJyISHFs3RESC4y0Q9igUCjhy5AiOHTt2UPUQEfXVOLZuJto9uba2hk8//bS6vLu72/eCiIj6qfzdi8APUbQN+qNHj+KTTz6pLteGPhHRKCo/Lwd+iKJt60ZRFPzoRz+CLMtQVRWu6+LZs2f44Q9/eFD1ERGFaphG6qZpIhaLAQA8z4Ou64G2AQDf9wEAqVSq4zZtgz6ZTOLJkycoFArI5/OwbRunTp1CJBJBPB5HIpGoPmRZxiuvvNLxgEREgzQss24qIa9pGgDAsiyYptk27A3DqAt23/eRy+WwtLTU9liBZt0kk0kkk0kAQDqdhuu6KJVK2NzcxOPHj/HZZ59ha2sL8XgcMzMz0HWdoU9EQ2lYboFgmibu379fXdY0DXNzc22D3nXdumVJkrCzs9PxWF1Nr6wUIMsyZFnG9PR03fOlUgm2bWN5eRkzMzOYmprqZvdERH3XzYje9/1qi6SWJEmQJGnfNfi+3xDaAOA4Dnzfb7nvYrEI13Uhy3JXx+sq6DvtPB6PIx6Po1AoYGtrq6tCiIgOQjdBf+fOHdy8ebNh/aVLl3D58uV911AsFqu9+VqSJKFYLEJV1abbpdNpzM3NIZ1OQ9d1WJaFdDrd8Xihf2FqZ2cH586dC3SBgIjooL0oBw/6s2fPYnZ2tmF9L6N5AC1H7c3Cv5aqqvj8888xNzcHwzCwuroaqJbQgz4ajeKbb74Je7dERKF4UQ4+bbLXFk07zVpCnbiuC9u2sbGxgZWVFUxOTmJ1dbXlJ4AK3gKBiMbK8y5G9EHYtg3DMDq+LhqN4saNGwBafyLwPK/tPgzDQDabBQBks1lomobFxUU8evSo7QmJQU9EY+W7F89D3Z+qqh1H1HslEommoe77PhKJRNNtbNtuOI6qqrhw4QJs265O02yGQU9EY6WbHn2/SJKEWCzW0KuXZbnlyLxVq0dRlI7Ha3sLBCIi0bwovwj86Cdd17G2tlZdNk2zbhKL4zhYWFioLmuahgcPHjTsJ5/Ps0dPRFSrm4ux/ZRKpWCaJmzbro7Wa78s5XkeCoVC3ah/aWkJuVwOR44cgSRJ8H0fFy9e7HgsBj0RjZVhaN1UtPsWrKqq2NjYqFsny3LH2x00w6AnorHyp5Avxo4CBj0RjZVhGtEfFAY9EY0VBj0RkeCG5WLsQWLQE9FY4YieiEhwYd8CYRQw6IlorIR9C4RRwKAnorHC1g0RkeAY9EREguOsGyIiwfFiLBGR4Ni6ISIS3DjOuomUy2PYsCIiGiP84REiIsEx6ImIBMegJyISHIOeiEhwDPoRlMlkBl0CEY0QBv0IcRwHc3NzKBaLgy6lL3gCI+oPzqMfIYqiQNd1OI4z6FJC5TgOrly5MugyiITFEf2IyefzUFV10GWEqnICSyQSgy6FSEgM+hFTKBSEC3pAzBMY0bBg0I8Q13UBALZtw7Is5HK56rpRJ+oJzHVdZDIZvP7665icnGx43jRNnDhxApOTkzAMYwAVhse2bRiGgUwmA9u2B10O1RC2R5/L5fDFF1/A9/269el0GqlUakBV9abyj0dVVUiSBFmWsbi4iPv37w+4st7UnsAAYHNzE6dPn4Ysy4MsKxSyLCObzSIajeLWrVuwbbvuhKbrOnzfh67rkCRpgJX2bnFxERsbG1hYWEAsFht0OVRDyBH9uXPncOTIEWxsbOD+/fuQJAkPHz7Eb3/725ENeeBleyOdTlcDQZKkuguzjuOM5Kiw9gSmaRpmZmawuLg44KrC4zgOZmZmoKoq7t692/C8JEkjH/K2bVdPzDdu3ICiKAOuiGoJF/SGYUCW5WqgK4qCd955B5ZlDbiy3u1tbziOUw0Iy7Jgmia2t7cHVd6+tTuBVd5XJpMZ2b/DYrFYveC8vr7e8ClTBPl8HslkctBlUAvCBf3KykrDqH17e3vkR0yu68L3/bp2hm3b1ZkqmqaNbI+71QnMdV1IkgRd15HNZnHlypWRviahaRokSYJpmtV1juOM9Gwj13WrbdLKf4t4Iht1QgW94ziQZbmht1soFEb6H1PF3ve1traGdDo9oGrC0e4E5jhOXasjmUyO/HcI5ufn64K+MtofVbIsY2lpCb7vY2lpCUtLSyM/qBKRUEHveV5DGFZCo/Yfk+/71VkrldcMe2977/vK5XKYn58f6ZCoaHUC0zQN165dq64vlUojd4F274j99OnTcF23el1ChNFv5VPWqP3djBOhgj6RSMDzvOqy7/swDAPZbLbudbZtQ9M0lEolABiZlkc6nUYul4NhGDh69CiWlpYGXVLPOp3AKqNDx3EQj8dH7sS2d8QuyzIURal+UhFh9Os4zsj9vYwboaZXSpKEdDoNwzAgSRKePn2K69evN/xj0jQNtm0jHo9X141C2GuaBk3TBl1G6ConsCNHjuDo0aPQdb3hNaZp4saNGwOoLnwXLlzA4uJi3UyVUba5uSlEa1RkQgU98DKwg4S24zg4efJk9b9F+Ac3qjqdwAzDqF6L8H1/5EfBlYuymUwGDx8+HHQ5PXNdt/pviYaTUK2bbtSOQjzPG/nwsG0bDx48QKFQqLvYN+ps265+mcj3/ZG6c6dlWS1HuvPz8yP//7mKUqnEEf2QG9sfB7dtG67rIhaLVb9pSsOlclvmWg8fPhz6T1++72N5eRmFQgGyLCOdTjd8ynRdF5ZljfQX+Cpef/11bGxs8N/QEBvboCei3lVuMT3qt+EQ3di2boho/yrTQm3bxjvvvDPgaqgTBj0Rde3EiROwbRubm5tNZ0nRcGHQE1HX3n//fbiui5mZGfbmRwB79EREguOInohIcAx6IiLBMeiJiATHoCciEhyDnohIcP8fexLmk2LAeEYAAAAASUVORK5CYII=\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:'$\\sigma_v$',1:'$b_1$',2:'$b_2$',3:'$N$',4:'$f$'}, inplace=True)\n",
    "df.rename(index={0:'$\\sigma_v$',1:'$b_1$',2:'$b_2$',3:'$N$',4:'$f$'}, inplace=True)\n",
    "ax = sns.heatmap(df, cmap=\"PiYG\", vmin=-1, vmax=1)\n",
    "\n",
    "ax.set_title('$\\mathrm{TNS \\, Lor}, z=0.5$')\n",
    "#plt.savefig('corr-mat-tns-zeq0p5.pdf', format=\"pdf\", bbox_inches = 'tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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>$\\sigma_v$</th>\n",
       "      <th>$b_1$</th>\n",
       "      <th>$b_2$</th>\n",
       "      <th>$N$</th>\n",
       "      <th>$f$</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>$\\sigma_v$</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.618788</td>\n",
       "      <td>0.211468</td>\n",
       "      <td>-0.060455</td>\n",
       "      <td>0.974429</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$b_1$</th>\n",
       "      <td>-0.618788</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.323727</td>\n",
       "      <td>-0.480270</td>\n",
       "      <td>-0.465335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$b_2$</th>\n",
       "      <td>0.211468</td>\n",
       "      <td>0.323727</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.983931</td>\n",
       "      <td>0.384652</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$N$</th>\n",
       "      <td>-0.060455</td>\n",
       "      <td>-0.480270</td>\n",
       "      <td>-0.983931</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.245635</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$f$</th>\n",
       "      <td>0.974429</td>\n",
       "      <td>-0.465335</td>\n",
       "      <td>0.384652</td>\n",
       "      <td>-0.245635</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            $\\sigma_v$     $b_1$     $b_2$       $N$       $f$\n",
       "$\\sigma_v$    1.000000 -0.618788  0.211468 -0.060455  0.974429\n",
       "$b_1$        -0.618788  1.000000  0.323727 -0.480270 -0.465335\n",
       "$b_2$         0.211468  0.323727  1.000000 -0.983931  0.384652\n",
       "$N$          -0.060455 -0.480270 -0.983931  1.000000 -0.245635\n",
       "$f$           0.974429 -0.465335  0.384652 -0.245635  1.000000"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add Gaussian priors on parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sv index 4\n",
      "b index 1\n",
      "N index 3\n",
      "combined marged [2.728e-02 6.873e-03 9.413e-02 4.639e+01 3.397e-01]\n",
      "\n",
      "Correlation coefficients after adding parameters priors:\n",
      "[[ 1.    -0.717  0.661 -0.004  0.99 ]\n",
      " [-0.717  1.    -0.693 -0.08  -0.656]\n",
      " [ 0.661 -0.693  1.    -0.618  0.654]\n",
      " [-0.004 -0.08  -0.618  1.    -0.033]\n",
      " [ 0.99  -0.656  0.654 -0.033  1.   ]]\n"
     ]
    }
   ],
   "source": [
    "Priormat = Fishermat*0.\n",
    "Prior_errors = {\n",
    "    'sv' : (sv_fid/100)*10,\n",
    "    'b': (b_fid/100)*10, \n",
    "    #'b2': (b2_fid/100)*10, \n",
    "    'N': (abs(N_fid)/100)*10\n",
    "}\n",
    "for p in Prior_errors.keys():\n",
    "    if p not in varyzdep: continue\n",
    "    print (p, 'index', fm_index[p][0])\n",
    "    Priormat[fm_index[p][0],fm_index[p][0]] = 1./Prior_errors[p]**2\n",
    "\n",
    "Fishermat_comb = Fishermat+Priormat\n",
    "CovMat_comb = np.linalg.inv(Fishermat+Priormat)\n",
    "print ('combined marged',  np.sqrt(np.diag(CovMat_comb)))\n",
    "\n",
    "print ('\\nCorrelation coefficients after adding parameters priors:')\n",
    "print (CovMat_comb/np.sqrt(np.outer(np.diag(CovMat_comb),np.diag(CovMat_comb))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEUCAYAAAAlXv26AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3U9sG+fdJ/Av5bTZF2sNtZfgxTbjw/umcUxSWLxI1K1Gp0aqNBFeYCEFFXXZ2IZD52JJaUoB70GmW1YFClN4GxvGvrBo13EOLzwx6myBRhyldrMHc4i1CgOFOUqLOIdqlPewPkQz0uJ1t4m5B4MDUfw3FIcacvj9GATMIWfmxz/6zsPneTgMFAqFAoiIyLd6vC6AiIhai0FPRORzDHoiIp9j0BMR+RyDnsgjlmVB13Wvy6Au8IzXBbSSYRgYGRlBOBxGJBJBb28v1tfXoWkaotEoent7sb29DU3TYJom1tbWAACpVAq5XA6GYeDOnTsQBMHepmVZWFhYwOrqKiRJQjKZhCiKAABVVfHgwQP09fUBAERRRDgchq7rkGW5rL5EIoFMJgPLsjA2Nobx8fGK9+s2mqbBMAyIomg/d8XnuBrLspDJZCBJkr2eoiiIx+Mlr18rKYoCRVFw69atuvfVdR1LS0vI5/P2++4gNFIj+UjBx7LZbOH8+fNly1588cWy+05MTJRcv3HjRuHEiROFmZmZits+e/Zs2f2Xl5fL9jUxMVG4ceNG1RrPnz9ftu9utrGxUfacnzhxou56+Xy+8OKLL9qX4eHhQj6fb1WZVWvY+x4oLq/ENM3CK6+80uqyymqpVCP5m69b9JZlYXp62tF9o9EoLMsqaf0lk0mMjIxA13WEw+GS+++9nk6ncfv27ZJlkiQhEonss/rulE6ny16zUCgEVVXrftq5ffs2gsEgABxYK363cDhc9r4Ann5CqbTcC9VqJH/zdR+9aZp1P/IXiaII0zTLlkWjUZw9e7bmupZlVb0tGo062j89lclkyg6O/f390DTN0fqCIHgS8rVkMhmvS6Au5+ugb6Q1HYlE7NbgbvF4HLquQ1XVqusKggDDMCoOrBXHB6g+wzDKPlUBT59fp0HvJV3XS+q0LAuJRAKGYXhYVam9NVJ38HXQN/IRtVpLUBAExOPxuq36N998EydOnICiKGUtfLc+KquqCkVRoKoq0ul0SYDouo7JyUlMTk7CMAyoqoqTJ0/W/LRRSSqVwsDAABKJBICnYaUoiiv111Ot1mAw6Cgsiwfk4vOz3xktqVQKR48excDAAHRdh2EYmJycxNGjR+3nBQBGRkYwMDAAVVVhGAYuX76Mubk5+/ZiS96yLKRSKaRSqYoha1mWXXe1A0Px9R0ZGbFfX0VRkEgkSh5nrfdBpRqLmnlvqaoKTdOgaZr9GKjNeD1IcNCqDcbutXcAdXh4uGQQq9IA6/nz5+3BwImJCUeDXk4HYysN9s7MzJQM9OXz+cLw8HAhm80WCoWnA8YbGxt1t12UzWYLZ8+eLWQymcLy8rJdm2maNdc7ceKE48vewfG9+6/02hQHWmvZ2NgoZDKZkmUTExMNPf7dzp8/XzbgvnfgNJvN2s91oVB5cHVjY6PqgGvx/ru3kc1mqw4+F5+HvfscHh4ueZy13geVamzmvZXJZMqedyeD53SwGPRV7A3y4nrF0Ks1k6YYmK+88krdEHcS9NVmZxRn9ey+3+4aG7V3vUYPFM1qJugrqRTWTlUK6OHh4ZJQ2xtwhUL5waBe0O99vWrNxKl22/LycslMpXrvg93baPa9tby8XPa3UOtvg7zh664bN0mSBEmSsLS05Oi+yWQSd+7cAfC0K2C/LMtCPp+vOKgsSVJZ90Qzg5G710ulUo7mr7up0hgJgLJBcqeOHDmy7/5oURQRDAbt51fTNMTjcaysrOxre7Xsfb0a7W6TJAmrq6tl23TyPmj2vSXLMpaWlnDy5Em725ITENoPg74ByWQSiqKU/QFU+4ajIAg4ffo0crncvvZXDPl6fc27b68Wlo1IpVIYHx+HJElNb6sRxcDZG3SVBmj3SqfTZcuCweC+DxIAMDY2Zge7ZVmQZdl+LS3LOtCDYC3F52b38+b0fdDse0sURdy5cweSJEFRFAwMDFR8Lchbvp5H77bd0y13t1pM06w6V1qSJCwvL+9rf4ZhIBgMVg2U4h+2m4GTSqUwNDTU0ADyyZMnHd83FAphfn6+4m3FFqNpmiXBbpomBgcHq27TMAwsLS0hGo2WrdfMczM+Po65ubmSegcHB6Fpmh38jdI0zfUDaPFAuJ9Pcs2+t4rfMYnFYojFYrAsC5OTkwf+aZBqY9A3KB6PY2BgANlstiTsM5kMYrFY2f0Nw9j39MqVlRVMT09DkqSKs37y+TzC4bBr88aLIb87iNLpdMXHtdu1a9dc2T/wNEh1XS8JCV3XMTQ0VHUdURSRTCbLnodsNlvzAFFPOByGaZpIp9P2az0+Po4bN25gfHzc0Tb2toKbnWpZqVtH07R9P85m31srKysQRdG+nyAIiMVi9iksqD2w66aKah9pi9Mtt7e3S5YXp9DttbS0hHg83vD+NU3DlStX7D+i06dPl30kTqfT+OlPf2pfr9VNYRgGBgYGqvZZq6qK9fV1ZLNZpNNpqKqK2dnZhutu1vz8fNl0znw+X3JQLZ7DaO9j2TslcH19vaQ1Xu85qGRqagqaptlBVuy+qRSAlmWVBXGlbpVm7Z2rrygKFhcX7WW13gd7a3TjvXX58uWyffC7I+3l0I9//OMfe13EQdB1HZcuXcL7778Py7Jw//59fPbZZ+jv78ezzz5bct9EIgFFUXD//n18+9vfxnPPPVdy+8svv4zPPvvMbmUGAgGEw2EMDQ1heXnZnlO8urqKeDxetv7u/fz617/G5uYmHj16hN///vd2wF+6dAkAMDMzY+/z0aNH+OSTT/Dw4UPkcjlEo1G7i0XXdfzsZz+Drut49OgRnnvuuZL9FgNBEISKreNPPvkEyWQS/f39SKfTUBQF4+PjdVvzbgsGg3juueegqioePXqEDz/8EPPz8yUtY8uy8P777+O1116zW42RSASqqiKXy+F3v/sd7t+/jwsXLpS8tvWeg0r6+vogiiJeeOEFe5lpmmWnaTAMA4uLi/j888/x6NEjfO9737Nvi0Qi+PDDD7G9vY2XX37ZHuRNpVL2/V944QVsbm6WLHv55ZdL6v/LX/4CRVEwPT2Nhw8fQtd1rKysYGFhwX6ta70PqtXYzHtrc3MTQ0NDuH//Ph4+fIiHDx+WPV/kvUChwN+M7SaKonT9rIhOfQ4sy8Lw8PCBnu2S/IFdN0REPseg7yKtmPHRaTr5OWhmqih1NwZ9F2mnud9e6dTnQNd1JBIJ+0Rpbg7ukv+xj56IyOfYoici8jnPvzD1Uuw/e12C6258/k9el9AS783+q9cluO7n8v/yuoSW+Ie5v/O6hJb49PK/Nb2NRjLnj+nm99cO2KInIvI5z1v0REQH6dAzAa9LOHAMeiLqKoHuy3kGPRF1l54u7LBm0BNRVwn0dF+TnkFPRF2FLXoiIp8LMOiJiPzt0CF23RAR+Rq7boiIfI5dN0REPtfDWTdERP7GL0wREfkc++iJiHyO57ohIvI5DsYSEflcTxd20u876K9cuYJAIABJkiCKIg4fPuxmXURELdGqFr2iKAgGgwCe/pB7NBp1tA4A+zeAY7FYS2rbd9CHw2EMDg4CANbX12GaJvr7+xn4RNTWWjEYWwx5WZYBAKqqQlGUmmGfTqdLgt2yLKRSKczPz7te374fcjHkASAUCmFwcBCaprlSFBFRqxx6JuD44pSiKHbIA4Asy3ZrvRrDMEquC4KA7e3txh6MQ/tu0S8tLeHTTz9FKBSCJEno7+/H1taWm7UREbnO7S56y7LKQhsAdF2HZVkQBKHievl8HoZhQBRFdwuqYN8t+v7+fly9ehWSJEHXdSwsLBxIwUREzejpCTi+WJaFzc3NskuxTx14GtjFvvndBEFAPp+vWkc8Hsfk5KTd8ldVFfF43P0HjCZa9JIk4eOPP8bo6GhJNw4RUTtrZDD2+vXruHTpUtnyM2fOYGZmBgCqttorhf9ukiThvffew+TkJNLpNK5du1a19d8sR0F/9epVnDp1qmRZb28vRkdHW1IUEVGrNHKum+NvHMfExETZ8r2BvLuF75RhGNA0DWtra1haWsLIyAiuXbsGSZIa3lY9joJ+aWmpLOiJiDpRI/PoBUGo28qudrtpmjXXS6fTSCaTAIBkMglZljE3N4c7d+643rJ3FPSFQgFffPEFVFWFIAiIRCI4duyYq4UQER2EZ55xd35lJBKpGOqWZSESiVRcR9O0spa7JEk4ffo0NE0rmcHjBsePeHJyEn/+859x9+5dHD9+HKdOncLm5qarxRARtVpPoMfxxQlBEBAMBsu6b0RRrNoyr9bVEw6HW9JP7zjor127hmQyiQsXLuDevXs4deoUzp07h52dHdeLIiJqlUZm3TgVjUaRyWTs64qilHwZStd1zM7O2tdlWcbKykrZdrLZrHd99IIgIBQKlSwrzp2/fPkyfvSjH7leGBFRK7Tih0disRgURYGmaXZrffe3Yk3TRC6XK5mhMz8/j1Qqhb6+PgiCAMuy8NZbb7leG+Aw6KempqrOvOHceSLqJK06qVmt0x1IkoS1tbWSZaIotuR0B5U46rqJx+PIZrM4d+4cvvjii5LbWvWVXSKiVujp6XF88QvHX5j65S9/iXQ6jeHhYRw5cgSiKGJra6tl3+QiImqFZw75J8CdauibsbFYDLFYzD5bZSQSQW9vb6tqIyJyHX8c3KG9A7NERJ3CT10yTvEXpoioqwT4C1NERP7GrhsiIp975tAhr0s4cAx6IuoqbNETEflcq74w1c4Y9ETUVdiiJyLyOU6vJCLyOXbdEBH5nNs/PNIJPA/6G5//k9cluG7673/udQktcfNqwusSXPfZ1fdw9P3/7nUZrvvwsT/fg25w+oMifuJ50BN5yY8hT7VxMJaIyOd6AvzCFBGRr3HWDRGRz33j0De9LuHAMeiJqKuw64aIyOfYdUNE5HOH2KInIvK3nh4GPRGRr/ELU0REPsdZN0REPscWPRGRz7GPnojI5ziPnojI5w6xRU9E5G/soyci8rlnOOuGiMjfWnUKBEVREAwGAQCmaSIajTa0fiKRQDKZbEVp6L7PMETU1Q4FDjm+OFUMeVmWIcsygsEgFEVxvL6qqtA0bT8PxxHXgv7mzZtubYqIqGV6AoccX5xSFAWyLNvXZVl2HPSWZTX8GBrluOvmypUrME2z4m2FQgG5XA4/+MEPXCuMiKgV3J5Hb1kWDMMoW67rOizLgiAINdfXNA2SJLla016Og16SJNy4cQNDQ0MVb6/0QImI2s03Dn3D1e3l83m7b343QRCQz+drhriu6y0PeaCBoA+FQhgaGsLY2FjF2+sdtYiI2kEjXTKWZVXsWhEEwc68aq32SuG/l2maCIfDLe++aWjWTbWQB4DBwcGmiyEiarVG5tFfv34dly5dKlt+5swZzMzM2Nf3E9Sqqpb067dS3aA3DANzc3P49NNPIQgCIpEI5ufn8dJLLx1EfURErmqkRX/8+HFMTEyULd/dgq/Wm1FtTBN4mqvhcNhxHc2qG/SKomBxcRG9vb1YX19HNpvFG2+8genpabzzzjsHUSMRkWsaCfrdXTTVRCKRiqFuWRYikUjFdTRNK/kUsLW1BdM0kU6nEQ6HXe+3rxv0/f39CIVCAABRFDE2NoZkMomlpSXcvHmTM22IqKMEXD6pmSAICAaDZX31oihWPUjs/TKVYRhYXV1FLBZztbaiks6qRCKB119/HefOncPHH3+MnZ0dBAIB7OzslK0Yj8cPZP4nEZGbngl80/HFqWg0ikwmY19XFKUktHVdx+zsbM1t1OrqaVZJ0CeTSRQKBYTDYXz00Ud49dVXsbCwgOHhYdy8eRObm5slK3/55ZctK4yIqBVa8YWpYqhrmgZVVQGUttpN00Qul6vYOFZVFYlEApZlIZVKQdf1Jh9hubKum8XFRYRCIUxNTQF4+pFC0zTcvXsXqVQK29vbEEURwNNWPRFRJ2nV+ehrndtGkiSsra1VvK142oRWKgv6Yn98kSiKiEaj9oMofjGqGPZERJ3E7T76TtDw2SsZ8ETUyXrAoCci8jX+8AgRkc8908MfHiEi8jX20RMR+Rz76ImIfK5V0yvbGYOeiLoKg56IyOcOcTCWiMjf2EdPRORznEdPRORz7KMnIvI5Bj0Rkc8F2Ed/8N6b/VevS3DdzasJr0toiR/8bdLrElz3m8R/8rqEllie/BevS2iJf8YbTW+Dp0AgIvK5ADgYS0TkcwGvCzhwDHoi6ips0RMR+VyALXoiIr9ji56IyNc4j56IyPfYdUNE5GscjCUi8jkOxhIR+R5b9EREvsbz0RMR+R67boiIfI2DsUREPsfBWCIi32OLnojI19h148Dm5ibW19chiiKOHTsGANje3kY+n4coinj++eddL5KIyC2BFp0CQVEUBINBAIBpmohGoy1ZZz8aOrStrq5icnISH330EWZnZ/H2228DAHp7exEOh/H973+/JUUSEbkl0MA/p4qBLcsyZFlGMBiEoiiur7NfDQW9pmm4d+8eLly4gN/+9reIxWI4d+4cdnZ2IAgCCoVCS4okInJLAD2OL04pigJZlu3rsiw7CvpG19mvhoI+HA6XXf/JT34CRVGwvb2NQKD7RrOJqNMEGrjUZ1kWDMMoW67rOizLcm2dZjQU9MFgELlcDm+//TZ2dnbs5adOnUI+n2eLnojaXiMtesuysLm5WXbZHcb5fN7uZ99NEATk8/mKNexnnWY0NBg7NjaGzc1NjI+P4/DhwyW3DQ4O4tatW64WR0Tktka6ZK5fv45Lly6VLT9z5gxmZmYAPG2dC4JQdp9KQV60n3Wa4SjoDcPA3NwcPv30UwiCgEgkgiNHjuCll14quV8oFGpJkURErmmg4+H48eOYmJgoW743pPfT3dKKLppqHAW9oihYXFxEb28v1tfXkc1m8cYbb2B6ehrvvPNOq2skInJNoIEuZkEQKra8996nEtM0XV2nGY4+w/T39yMUCkEURYyNjSGZTOLevXt48uQJbt682ZLCiIhaotDAxYFIJFIxoC3LQiQScW2dZpQFfSKRwOuvv45z587h448/xs7ODgKBQMnga1E8Hj/Qjx9ERE1zOegFQUAwGCzLQlEUq7bc97NOM8qCPplMolAoIBwO46OPPsKrr76KhYUFDA8P4+bNm9jc3Cy5/5dfful6UURELVMoOL84FI1Gkclk7OuKoiAWi9nXdV3H7OxsQ+u4qWIf/eLiIkKhEKampgA8HYzVNA13795FKpXC9vY2RFEE8LRVT0TUKQJP3N9mLBaDoijQNM1upe8+nYFpmsjlciWzbeqt46ZAYR+T34sT/Yth34y3/+d3m95Guzl19YTXJbTED/426XUJrvvNf/iF1yW0xP8YedfrElrin/9brult/L+t/+v4vt/s+49N768d7OvslW4EPBGRJ7rwe53dd75OIqIuw/PRE1FXaWQevV8w6Imou3RfzjPoiajLPOm+pGfQE1FXCXRfzjPoiajLMOiJiHyOQU9E5HOcdUNE5G/soyci8ju26ImIfK77cp5BT0TdpRu7bvZ19ko3PX782Mvdt0TPv3/tdQktsZH4jdcltMQ/Pv6h1yW47g+/+NzrElri2cN/0/Q2nmxsOb5vz5G+pvfXDtiip67mx5CnOrqwRc+gJ6KuUvi6Bb880uYY9ETUXbov5xn0RNRdCjypGRGRz3EePRGRv7FFT0Tkdwx6IiJ/K3zNoCci8jWPvyPqCQY9EXUXdt0QEfkbB2OJiPyOXTdERP5W+IpBT0TkaxyMJSLyuy7so++pd4erV68eRB1ERAfjScH5xSfqBr2iKNjZ2TmIWoiIWq5QKDi++EXdrputrS0oioLe3l4EAgFEIhEcO3bsIGojInIfT1NcLh6PY2pqyr6+vr6ODz74AJZlIRgM4vnnn8fg4GBLiyQicks7/fCIoigIBoMAANM0EY1GG1o/kUggmUzWvV/doN8d8gAQCoUQCoXs67lcDqOjo+jt7cX4+DhOnTrVUKFERAepXb4wVQx5WZYBAKqqQlEUx2Gvqio0TXN0333PusnlclAUBaqqQhAEjI6OYnR0dL+bIyI6GA0EvWVZsCyrbLkgCBAEoakyFEXBrVu37OuyLGNyctJR0FeqqZaGgn5nZweKokBRFGxsbCAUCuHChQsYGxtraKdERF5pZJD1+vXruHTpUtnyM2fOYGZmZt81WJYFwzDKluu6Dsuy6h5ENE2DJEmO9+co6Iut99XVVfT29mJqagrRaBSiKDreERFRW2igRX/8+HFMTEyULW+2NZ/P5+2++b3bzefzNUNc1/WGQh5wEPSjo6MwDAODg4N499132Xonoo7WyGBs0IUumkqqtdorhf9epmkiHA431H1TN+g3Njbw5ptvYnx8vGQQloioEz3569delwCg8X524OkAbHHwthGOZt3E43F7WiWAqvPpc7kcp1oSUXt74u70Sk3TkE6n696vt7cXFy9eBFC968c0zarrG4aBcDi8rxrrBn0sFgNQPq1yb/CHw2EsLS3hV7/61b4KISI6CG7/lKAkSQ33mUcikYqhblkWIpFIxXU0TSv5FLC1tQXTNJFOpxEOh2vWUDfoqw247g3+1dVVbG5u1tscEZGnCi636PdDEAQEg8GyvnpRFKu29vdOuzQMA6urq3ZjvJa657pxamxsDN/97nfd2hwRUUsUvn7i+NJK0WgUmUzGvq4oSklo67qO2dnZmtuo1dWzm2tBDwCLi4tubo6IyHWFr752fGmlYqhrmgZVVQGUttpN00Qul6s4aKuqKhKJBCzLQiqVgq7rNfcVKHh8irbHjx97ufuW6Pn39hjVd9tG4jdel+C6f3z8Q69LaIk//OJzr0toiWcP/03T2/g/V/+34/s+d+q/Nr2/dsAfHiGirtIOffQHjUFPRN2ljc5eeVAY9ETUVdrl7JUHiUFPRF2lnc5Hf1AY9ETUVQptcgqEg8SgJ6KuwsFYIiKfY9cNEZHPsUVPROR3Lp/UrBMw6Imoqzz561del3DgGPRE1FXYdUNE5Hfsujl4/zD3d16X4LoPH//c6xJaYnnyX7wuwXV/GPbnyb/+yw//3usSWuKP6X9rehts0RMR+VzhKwY9EZGvcR49EZHPPfkLZ90QEfkau26IiHyOQU9E5HMFTq8kIvI3tuiJiHyOs26IiHzuyWPOuiEi8jW26ImIfI599EREPsegJyLyOXbdEBH5HE+BsEcul0NfXx+OHTt2UPUQEbVUN3bd9NS6MZPJ4PLly/b1nZ2dlhdERNRKha+eOL74Rc2gP3LkCN599137+u7QJyLqRIWvC44vflGz6yYcDuM73/kORFGEJEkwDANffPEFvvWtbx1UfURErmqnlrqiKAgGgwAA0zQRjUYdrQMAlmUBAGKxWN11agb94OAg7t27h1wuh2w2C03TMDw8jEAggFAohEgkYl9EUcThw4fr7pCIyEvtMuumGPKyLAMAVFWFoig1wz6dTpcEu2VZSKVSmJ+fr7kvR7NuBgcHMTg4CACIx+MwDAPr6+t48OAB7t69i+XlZWxubiIUCmF8fBzRaJShT0RtqV1OgaAoCm7dumVfl2UZk5OTNYPeMIyS64IgYHt7u+6+GppeWSxAFEWIooixsbGS29fX16FpGhYWFjA+Po7R0dFGNk9E1HKNtOgty7K7SHYTBAGCIOy7BsuyykIbAHRdh2VZVbedz+dhGAZEUWxofw0Ffb2Nh0IhhEIh5HI5bG5uNlQIEdFBaCTor1+/jkuXLpUtP3PmDGZmZvZdQz6ft/vmdxMEAfl8HpIkVVwvHo9jcnIS8Xgc0WgUqqoiHo/X3Z/rX5ja3t7GyZMnHQ0QEBEdtCcF50F//PhxTExMlC1vpjUPoGqrvVL47yZJEt577z1MTk4inU7j2rVrjmpxPeh7e3vxxz/+0e3NEhG54knB+bTJZrtoaqnUJVSPYRjQNA1ra2tYWlrCyMgIrl27VvUTQBFPgUBEXeXrBlr0TmiahnQ6Xfd+vb29uHjxIoDqnwhM06y5jXQ6jWQyCQBIJpOQZRlzc3O4c+dOzQMSg56IuspXT752dXuSJNVtUe8ViUQqhrplWYhEIhXX0TStbD+SJOH06dPQNM2eplkJg56IukojffStIggCgsFgWV+9KIpVW+bVunrC4XDd/dU8BQIRkd88KTxxfGmlaDSKTCZjX1cUpWQSi67rmJ2dta/LsoyVlZWy7WSzWfbRExHt1shgbCvFYjEoigJN0+zW+u4vS5mmiVwuV9Lqn5+fRyqVQl9fHwRBgGVZeOutt+rui0FPRF2lHbpuimp9C1aSJKytrZUsE0Wx7ukOKmHQE1FX+avLg7GdgEFPRF2lnVr0B4VBT0RdhUFPRORz7TIYe5AY9ETUVdiiJyLyObdPgdAJGPRE1FXcPgVCJ2DQE1FXYdcNEZHPMeiJiHyOs26IiHyOg7FERD7HrhsiIp/rxlk3gUKhCzusiIi6CH94hIjI5xj0REQ+x6AnIvI5Bj0Rkc8x6DtQIpHwugQi6iAM+g6i6zomJyeRz+e9LqUleAAjag3Oo+8g4XAY0WgUuq57XYqrdF3H2bNnvS6DyLfYou8w2WwWkiR5XYarigewSCTidSlEvsSg7zC5XM53QQ/48wBG1C4Y9B3EMAwAgKZpUFUVqVTKXtbp/HoAMwwDiUQCR48excjISNntiqJgYGAAIyMjSKfTHlToHk3TkE6nkUgkoGma1+XQLr7to0+lUvjggw9gWVbJ8ng8jlgs5lFVzSn+8UiSBEEQIIoi5ubmcOvWLY8ra87uAxgAPHjwANPT0xBF0cuyXCGKIpLJJHp7e3HlyhVomlZyQItGo7AsC9FoFIIgeFhp8+bm5rC2tobZ2VkEg0Gvy6FdfNmiP3nyJPr6+rC2toZbt25BEATcvn0bf/rTnzo25IGn3RvxeNwOBEEQSgZmdV3vyFbh7gOYLMsYHx/H3Nycx1W5R9eLFlLKAAADMUlEQVR1jI+PQ5Ik3Lhxo+x2QRA6PuQ1TbMPzBcvXkQ4HPa4ItrNd0GfTqchiqId6OFwGK+99hpUVfW4subt7d7Qdd0OCFVVoSgKtra2vCpv32odwIqPK5FIdOxrmM/n7QHn1dXVsk+ZfpDNZjE4OOh1GVSF74J+aWmprNW+tbXV8S0mwzBgWVZJd4amafZMFVmWO7aPu9oBzDAMCIKAaDSKZDKJs2fPdvSYhCzLEAQBiqLYy3Rd7+jZRoZh2N2kxf/78UDW6XwV9LquQxTFsr7dXC7X0X9MRXsfVyaTQTwe96gad9Q6gOm6XtLVMTg42PHfIZiamioJ+mJrv1OJooj5+XlYloX5+XnMz893fKPKj3wV9KZploVhMTR2/zFZlmXPWinep937tvc+rlQqhampqY4OiaJqBzBZlrG4uGgvX19f77gB2r0t9unpaRiGYY9L+KH1W/yU1WmvTTfxVdBHIhGYpmlftywL6XQayWSy5H6apkGWZayvrwNAx3R5xONxpFIppNNpHDlyBPPz816X1LR6B7Bi61DXdYRCoY47sO1tsYuiiHA4bH9S8UPrV9f1jntduo2vplcKgoB4PI50Og1BELCxsYELFy6U/THJsgxN0xAKhexlnRD2sixDlmWvy3Bd8QDW19eHI0eOIBqNlt1HURRcvHjRg+rcd/r0aczNzZXMVOlkDx488EXXqJ/5KuiBp4HtJLR1XcfQ0JD9fz/8wXWqegewdDptj0VYltXxreDioGwikcDt27e9LqdphmHYf0vUnnzVddOI3a0Q0zQ7Pjw0TcPKygpyuVzJYF+n0zTN/jKRZVkddeZOVVWrtnSnpqY6/j1XtL6+zhZ9m+vaHwfXNA2GYSAYDNrfNKX2Ujwt8263b99u+09flmVhYWEBuVwOoigiHo+Xfco0DAOqqnb0F/iKjh49irW1Nf4NtbGuDXoial7xFNOdfhoOv+varhsi2r/itFBN0/Daa695XA3Vw6AnooYNDAxA0zQ8ePCg4iwpai8MeiJq2JtvvgnDMDA+Ps6++Q7APnoiIp9ji56IyOcY9EREPsegJyLyOQY9EZHPMeiJiHzu/wMteOSXY1bjXwAAAABJRU5ErkJggg==\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_comb/np.sqrt(np.outer(np.diag(CovMat_comb),np.diag(CovMat_comb))))\n",
    "df.rename(columns={0:'$\\sigma_v$',1:'$b_1$',2:'$b_2$',3:'$N$',4:'$f$'}, inplace=True)\n",
    "df.rename(index={0:'$\\sigma_v$',1:'$b_1$',2:'$b_2$',3:'$N$',4:'$f$'}, inplace=True)\n",
    "ax = sns.heatmap(df, cmap=\"PiYG\", vmin=-1, vmax=1)\n",
    "\n",
    "ax.set_title('$\\mathrm{TNS \\, Lor}, z=0.5, \\mathrm{\\, with \\, priors}$')\n",
    "#plt.savefig('corr-mat-tns-zeq0p5-priors.pdf', format=\"pdf\", bbox_inches = 'tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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>$\\sigma_v$</th>\n",
       "      <th>$b_1$</th>\n",
       "      <th>$b_2$</th>\n",
       "      <th>$N$</th>\n",
       "      <th>$f$</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>$\\sigma_v$</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.717498</td>\n",
       "      <td>0.660734</td>\n",
       "      <td>-0.004330</td>\n",
       "      <td>0.990139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$b_1$</th>\n",
       "      <td>-0.717498</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.692853</td>\n",
       "      <td>-0.080113</td>\n",
       "      <td>-0.655661</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$b_2$</th>\n",
       "      <td>0.660734</td>\n",
       "      <td>-0.692853</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.618372</td>\n",
       "      <td>0.654325</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$N$</th>\n",
       "      <td>-0.004330</td>\n",
       "      <td>-0.080113</td>\n",
       "      <td>-0.618372</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.032547</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$f$</th>\n",
       "      <td>0.990139</td>\n",
       "      <td>-0.655661</td>\n",
       "      <td>0.654325</td>\n",
       "      <td>-0.032547</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            $\\sigma_v$     $b_1$     $b_2$       $N$       $f$\n",
       "$\\sigma_v$    1.000000 -0.717498  0.660734 -0.004330  0.990139\n",
       "$b_1$        -0.717498  1.000000 -0.692853 -0.080113 -0.655661\n",
       "$b_2$         0.660734 -0.692853  1.000000 -0.618372  0.654325\n",
       "$N$          -0.004330 -0.080113 -0.618372  1.000000 -0.032547\n",
       "$f$           0.990139 -0.655661  0.654325 -0.032547  1.000000"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "varying zdep = \n",
      "\t ['f', 'b', 'b2', 'N', 'sv']\n",
      "at z = 0.5 : zdep fiducials = \n",
      "\t [ 7.770e-01  1.498e+00  4.233e-01 -4.684e+02  9.544e+00]\n",
      "at z = 0.5 : marged err = \n",
      "\t [2.728e-02 6.873e-03 9.413e-02 4.639e+01 3.397e-01]\n",
      "at z = 0.5 : unmarged err = \n",
      "\t [2.524e-03 1.384e-03 1.350e-02 1.007e+01 3.109e-02]\n",
      "at z = 0.5 : marged % err = \n",
      "\t [ 3.511  0.459 22.237  9.905  3.559]\n",
      "at z = 0.5 : unmarged % err = \n",
      "\t [0.325 0.092 3.19  2.15  0.326]\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_comb))[zdepind[:,iz]]\n",
    "        print ('at z =',z,': marged err = \\n\\t', merrs)\n",
    "        uerrs = np.sqrt(1./np.diag(Fishermat_comb))[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": {},
   "source": [
    "## Draw Ellipses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdgAAAHACAYAAAABXHlvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XtUVPfBLv5nLtxhz4CiKAyCgkZAk6ppDG2atjFvsGl7+ppG0mqrNjGXvicmTdSsdX69pPb8Ey+pSd/TxqAhNvrWSdI0vZxAGtP2bZoxqdEkyngDBRlQLgIzwx3m8vvDA3WEAfaevdl7D89nLdYKw3z3PGSEh337fg3BYDAIIiIikpVR7QBERETRiAVLRESkABYsERGRAliwRERECmDBEhERKYAFS0REpAAWLBERkQJYsERERApgwRIRESmABUtERKQAFiwREZECzGoH0Jq+vj5UVVUhPT0dJpNJ7ThERKQyv9+P1tZWFBUVIT4+fsLjWLDXqaqqwpo1a9SOQUREGnPw4EEsW7Zsws9nwV4nPT0dwNX/kRkZGSqnIVLX8uXL8cEHH6gdg0hVTU1NWLNmzXA/TBQL9jpDh4UzMjKQlZWlchoidQWDQf4cEP0/Yk8b8iInIiIiBbBgiSisiooKtSMQ6RYLloiISAEsWCIKa+XKlWpHINItFiwREZECWLBEREQKYMESUVhPPPGE2hGIdIsFS0RhPfnkk2pHINItFiwRhbVkyRK1IxDpFguWiMJqbm5WOwKRbrFgiYiIFMCCJaKwFi1apHYEIt1iwRJRWJWVlWpHINItFiwRhbV161a1IxDpFguWiMI6ePCg2hGIdIsFS0REpAAWLBERkQJYsEQU1rFjx9SOQKRbLFgiCuvEiRNqRyDSLRYsEYW1YcMGtSMQ6RYLloiISAEsWCIiIgWwYIkorGeeeUbtCES6xYIlorDWrl2rdgQi3WLBElFYmZmZakcg0i0WLBERkQJYsERERApgwRJRWCtWrFA7ApFuRVXBulwuOBwOtWMQRY39+/erHYFIt1Qv2MrKSvz4xz8e9Wt2ux2VlZWorKyE3W4fd1sOhwOPPfYYFixYgJtvvhmbNm2Cy+WSOzLRlLFu3Tq1IxDpllmtF3Y6nXjrrbfQ2dkJt9s94ut2ux0WiwUlJSUAMFyypaWlYbdpsVhw9OhReL1eCIKgWHaiqeLw4cNqRyDSLdX2YAsLC7FlyxYUFxeP+nW73T5crgBQUlIyob1YACxXIiJSneqHiEfj9XpHPbTrdDrh9XpVSERERCSOaoeIx1JVVQWLxTLicUEQUFVVFXavF0DIRU5OpxMbN25UJCPRVNDY2Kh2BCLd0mTBhjuHOlrpXquwsBAWi2V4rNfrxY4dO7Bly5ZRn+92u0ec/21qapKYmij6HDhwgNMlEkmkyUPEACQdCrbZbCHFXFJSgr1794Z9/u7du5Gbmxvycdttt0nKSxSNnnrqKbUjEOmWJgs23EVKHo9H0racTueoX3v88cdRW1sb8vHee++Jfg0iIqLrafIQcVFR0ahl6vV6UVRUNOoYp9OJ9evX4+jRoyO+Fq6wrVYrrFZryGNmsyb/lxARkc5odg/WYrGMOEx8/SHg68ds3rw55LGh8TabTZmgRFGuvLxc7QhEuqXJggWA0tJSVFRUDH9ut9tDrgh2Op3YtGnT8OejleiePXtGlC4RTdzixYvVjkA6FQgG4OnzoN5TjzNXzqCqpQrOFifq3HXo9/WrHW9SqHY81OVy4dChQzhy5AhcLhd27NiBRYsWDU8usXHjRtjtdjgcjuE90WtncfJ4PDhy5EjIFccrV65EWVkZBEFAfX09srOzx5z5iYjGtnTpUt6qQxMWCAbQ3tsOd58bnf2dCAQDI57T5+uDu8+NgvQCxJpiVUg5eQzBYDCodggtaWhowB133IF3330XWVlZaschUlVmZiYLlsblC/jQ0t2C1u5W+AK+CY3JSM5AppCpcDJ5SO0FXtFDRESS9Pv60dzdjLaetlH3Vscy4B9QKJV2sGCJKKw1a9aoHYE0qN/Xj0udl9DR1wGpB0H9Qb/MqbSHBUtEYW3fvl3tCKQh/oAfl7suo6W7RXKxDkkwJ8iUSrs0exUxEanv2hWtaOoKBoNo7W5FVUsVmruaIy5XAEiJS5EhmbZxD5aIwjp58qTaEUhl3n4vGrwN6B3slW2bBoMBybHJsm1Pq1iwREQ0Qr+vHy6vC54+8VPUjicxJhFGQ/QfQGXBElFYM2fOVDsCTbJgMIjm7mZc7rws+srgiZqWME2R7WoNC5aIwjp+/LjaEWgS9Qz2oM5dJ+vh4OuZjWZMS5waBRv9++hEJNmuXbvUjkCTIBAMoMHbgDNXzihargCQnpQ+JQ4PAyxYIhrDs88+q3YEUpi334tTradkuzp4LAaDAemJ6Yq+hpbwEDER0RTkC/jQ4G1AW0/bpL1mWkIaYkwxk/Z6amPBEhFNMZ4+D+rcdROeN1gOBoMBGckZk/Z6WsCCJaKwrl0ykvQvEAzA5XHhSs+VSX/t6YnTEW+On/TXVRMLlohoCuge6Eatu1aVtVhNRhNmp8ye9NdVGy9yIqKwVq5cqXYEilAwGMSlzks423ZWtYXOZyXPgtk49fbnpt53TEQ0RfT5+lDbUYuewR7VMsSZ4zAjaYZqr68mFiwRURRq6W5Bo7dRsdmYJipLyILBYFA1g1pYsEQU1hNPPKF2BBJp0D+IOncdvP1etaMgJS4F1nir2jFUw4IlorCefPJJtSOQCN5+L2o7aif19ptwDAYDsi3ZasdQFS9yIqKwlixZonYEmoBgMIhGbyOq26o1Ua7A1QubptptOdfjHiwRhdXc3Kx2BBrHoH8QFzouoGugS+0owxJjEqfcpBKjYcESEemUlg4JDzEYDJhjnTNlL2y6FguWiMJatGiR2hFoFEP3tjZ1NakdZYSM5AwkxiSqHUMTWLBEFFZlZaXaEeg6A/4B1HbUauqQ8JB4czxmJc9SO4Zm8CInIgpr69atakega3j6PDjdelqT5WowGJBjzeGh4WuwYIkorIMHD6odgfCvq4Rr2ms0db71WjOSZiApNkntGJrCQ8RERBqm5UPCQ+LN8VNyMv/xsGCJiDRKjXVbxTIYDJibOhdGAw+IXo8FS0RhHTt2TO0IU5KWrxK+XpaQhYSYBLVjaBL/5CCisE6cOKF2hClnwD+Ac23ndFGulnjLlF0pZyJYsEQU1oYNG9SOMKV4+72avUr4ejGmGORYc9SOoWk8RExEpAHNXc1o7GxEMBhUO8qE5FhzpuQi6mLw/w4RkYoCwQAuui+ivbdd7SgTNjN5JoQ4Qe0YmseCJaKwnnnmGbUjRLUB/wBq2mvQO9irdpQJS4xJRGZKptoxdIEFS0RhrV27Vu0IUauzvxMXOi5o+hac6xkNRuSm5nK2pgniRU5EFFZmJvdUlNDc1Yzqdu2s3TpRNottyq/xKgb3YImIJokez7cOmZY4DdMTp6sdQ1dYsEREk2DAP4Dz7efRM9ijdhTREmISkG3JVjuG7rBgiSisFStWqB0hKujxfOsQo8HIqRAlYsESUVj79+9XO4Lu6e3+1uvNsc7heVeJ+CcJEYW1bt06tSPoVjAYRJ27Dg3eBt2Wa3pSOtIS0tSOoVvcgyWisA4fPqx2BF0a9A/ifMd5dA90qx1FsqTYJNgEm9oxdI0FS0Qko57BHtS012DQP6h2FMliTDGYlzqP97tGiAVLRCST9t52XHRfRCAYUDuKZAaDAbnWXMSYYtSOonssWCIKq7GxUe0IutHobdTFEnPjyUzJREpcitoxogIvciKisA4cOKB2BM3zB/yoaa+JinJNTUjFzOSZaseIGixYIgrrqaeeUjuCpvX7+nHmyhl4+jxqR4lYQkwC13eVGQ8RExFJoOfJI65nMpo4mYQCWLBERCK1drfC5XXp9v7W6+VYcziZhAJYsEQUVnl5udoRNCUYDMLldaG1u1XtKLLJSM6ANd6qdoyoxIIlorAWL16sdgTN8AV8uNBxAZ39nWpHkY0QJyBT4JKESuEBdyIKa+nSpWpH0IQ+Xx/OXDkTVeUaZ45Dbmqu2jGiGvdgiYjG0DXQhfPt56PiYqYhJqMJeWl5MBtZAUri/10iojDae9tR566LmouZgKszNc1LnceLmiYBC5aIwlqzZo3aEVRzufMyLnVeUmz7voAP1W3VONt2Fl2DXbDGW3Hn3DsRZ4pT7DUBINuSzZmaJgkLlojC2r59u9oRJl0wGES9px5Xeq4osv3Gzka8e+FdOFud6Pf3h37N24iHlj6kyOsCwMzkmZieOF2x7VMoFiwRhVVSUoLKykq1Y0waf8CPCx0X4O33yr7tOncd3rnwDpytzrDPcXldsr/uEGu8FVlClmLbp5FYsEQU1smTJ9WOMGkG/AOoaa9B72CvrNutaa/BO+ffwdn2s+M+12QwyfraQxJjEnnFsApYsEQ05Smxhmt7bzteO/UaTl85PeEx81Lnyfb6Q2JMMchLy+M0iCpgwRJRWDNnRv/KKp4+Dy50XJBtDddgMIi/1/8db1W/NeIc63jy0/JlyTDEaDAiLy2Pa7uqhAVLRGEdP35c7QiKkntO4cudl2F32lHnqRM9NikmCUtmL5Elx5Dc1FwkxiTKuk2aOB4zIKKwdu3apXYExTR4G1DvqZelXH0BH96qfgu7PtglqVwB4PO2z8t6i06mkMk5hlXGgiWisJ599lm1I8guEAzgfPt5NHc1y7K9Cx0XsOP9HfjzhT9Lnu0pzhSHL+R8QZY8ADAtcRoykjNk2x5Jw0PERDRl+AI+1LTXoHugO+Jt9fn68Mdzf4TD5UAQke0FL89cjqSYpIgzAUBKXArmWObIsi2KDAuWiKaEPl8fatpr0O8Td+HRaKpaq/D6qdfh7nNHvC2z0Ywv5X4p4u0AVyfwn5s6FwaDQZbtUWRYsEQUVkVFhdoRZNHZ34nzHefhD/gj2o4v4MObZ97EP1z/kCkZ8KWcL8lyrpQT+GsP3wkiimpyTdjf0t2CX3/6azR0NsiUDMhIysBd8+6KeDucwF+bVC/YyspKOBwObNu2bcTX7HY7LBYLAMDj8aC0tHTc7UkZQ0SjW7lyJRobG9WOIZlcE/Z/dPkjvOZ8TfR9rWMxwID7iu6TZY+TE/hrk2pXETudTuzYsQMOhwNu98jzGENFWVJSgpKSElgsFtjt9jG3KWUMEUWfYDCIOnddxOU64B/AoapDOHDigKzlCgC3z7kdOdaciLfDCfy1S7WCLSwsxJYtW1BcXDzq1+12O0pKSoY/LykpmVDBih1DRNHFH/Cjur0abT1tEW2nqasJzx55Fh80fiBTsn+ZnjgdX8n/SsTb4QT+2qb6IeLReL1euFwjV5VwOp3wer0QBEGWMUQ0tieeeELtCKLINWH/Bw0f4I3Tb2AgMCBTsn8xwID7Cu9DrCk2ou1wAn/t02TBVlVVDZ9HvZYgCKiqqhp1r1fKGLfbPeLwdFNTUwTJiaLLk08+qXaECZNjwv5+fz9ec76Gjy5/JGOyUMW2YuSl5UW0DU7grw+aLNhwe5yjFWgkY3bv3o2f/vSnIY+ZzWbMnTtXRFqi6LVkyRJdzEfs7nOjtqM2ogn7GzobsP+T/WjtaZUxWajU+FR8bf7XItoGJ/DXD00WLHC1MJUe8/jjj2P9+vUhjzU1NWHdunWiX5soGjU3yzOdoJJaulvg8kS2UPk/6v+BN8++KXmqw4laXbA64ltpOIG/fmiyYMOdL/V4PLKOsVqtsFpDb/A2mzX5v4SIRuHyuNDS3SJ5fJ+vD4eqDuGT5k9kTDW6z87+LBamL4xoG5zAX1802SZFRUWjFqPX60VRUZFsY4hobIsWLVI7wqiCwSBq3bXo6O2QvI2W7hbsPb4XLT3SC3qiZiTOwKqFqyLaxvTE6ZzAX2c0WbCCIMBisYw4r2qz2cLuqUoZQ0Rjq6ysVDvCCP6AH+c7zqOzv1PyNqpaq3DgxAH0+fpkTDa6GGMM1t+0PqJDw9Z4K7It2TKmosmg2UvQSktLQ+ZBtdvt2Lhx4/DnTqcTmzZtEjWGiMTZunWr2hFC+AI+nGs7J7lcg8EgKmsqse/4vkkpVwC4Z+E9mJ0yW/L45Nhk5KbmcgJ/HVKtYF0uF3bs2IEXX3wRR44cwY4dO0L+Wh4qRofDMfz4tdMeejweHDlyJOTCpvHGEJE4Bw8eVDvCsH5fP85cOYOewR5J4/t8fdj3yT5Unq+MeHm5iVo2exmWZy2XPD4hJoG34+iYaoeIbTYbtmzZMuZzxirH4uJiHD16VNQYItKn3sFeVLdXS77HtaW7BS99/BKauifvPveMpAzcW3Cv5PGxpljkp+XDZDTJmIomkybPwRIRDeka6EJNe43kpeacrU4cOHEAvb7IZncSI9YYi/U3rUecKU7SeLPRjPxp+bzXVedYsEQU1rFjx1R9/UgnkDhcexj/99z/nbRDwkNWF62WfMXv0EQSXHpO/1iwRBTWiRMnkJGhzq0hV3quoN5TL2kd18HAIA5VHcKxy5P/B8IXsr+AZbOWSRprMBgwN3UukmKTZE5FamDBElFYGzZsUGU92KauJjR6pb2up9+Dfcf3od5bL3Oq8eVac/H1G74uefwcyxxY4sNP70r6woIlIk2JZHami56LeOnjl+DpDz+Dm1KEWAHrb1oPs0Har9VMIRPTEqfJnIrUxIIlIk0YWiS9vbdd0viPLn0Eu9OOwYD01XSkMhvNWH/TeljipO19zkyeyVmaohALlojCeuaZZybldQLBAM63n4e3X/wiH8FgEH8890f8pe4vCiSbmK/P/zrmpkpbhSstIY2LpkcpFiwRhbV27VrFX8MX8KGmvQbdA92ix/b7+/HrT38NZ6tTgWQTs2z2MnxhzhckjRXiBORYc+QNRJrB6UGIKKzMzExFtz/gH8DZK2cllaun34NffPgLVcs115qL+wrvkzQ2KTYJ89LmcQrEKMY9WCJSRZ+vD9Vt1RjwD4gee6nzEl48/iLcfW4Fkk3MtIRpuP8z98NsFP9rNM4cxykQpwAWLBFNuu6BbtS010ha4Pxs21mUf1I+aZP1jybBnICNSzYiOTZZ9NgYUwzmT5svqZhJX/gOE1FYK1askH2b3n4vzreflzQ70wcNH+C1U6/BH5Q2baIcjAYj1t24TtJVvyajCflp+Yg1xSqQjLSGBUtEYe3fv1/W7bX3tqPOXSd6dqZgMIi3at7COxfekTWPFKsWrsIN028QPc5gMGBe6jwkxCQokIq0iCcAiCisdevWybatlu4W1HbUii5XX8CHV068ooly/UL2F/B52+cljc215iIlLkXmRKRl3IMlorAOHz4sy3YavY1o6hK/VFz3YDde+vglnO84L0uOSBRML8A3bviGpLHZlmykJqTKnIi0jgVLRIoJBoOo99TjSs8V0WOv9FzBi8deREuPtGkT5TQ7eTa+e+N3JV31OytlFtKT0hVIRVrHgiUiRQSCAdR21Eq6labOXYe9H+9F10CXAsnESYlNwcYlGyUtH5eelI7ZKbMVSEV6wIIlorCkrqTjD/hR014jqSA/afoEB08eVGVO4evFGGNw/2ful3R41xpvhU2wKZCK9IIXORFRWAcOHBA9ZtA/iLNtZyWV619q/4L9n+7XRLkaYMCaRWskTWWYHJuM3NRcztI0xbFgiSisp556StTz+339ONt2Fr2DvaLGBYIBvH7qdfzh3B8QhPgF1pWwMm8lbsq4SfS4hJgEztJEAHiImIhk0jPYg5r2Ggz6xe199vv7sf+T/Th15ZRCycRbNnsZ/m3ev4keF2uKRX5aPkxGkwKpSG9YsEQUsc7+TpzvOA9/QNwMS55+D8qOlaGhs0GhZOJJncDfbDQjf1o+YkwxCqQiPWLBElFY5eXl4z7H0+fBhY4Loqc+bO5uxgsfvYCOvg6p8WSXnpiO+5eIn8DfZDQhf1q+pCuNKXqxYIkorMWLF4/59Y7eDtS6xc/OVOeuQ9nxMnQPil+mTimWOAseWfYIkmPETeBvNBiRl5aHxJhEhZKRXvEsPBGFtXTp0rBfa+tpk1Sup1tP41cf/UpT5RpvjseDSx9EWkKaqHEGgwHz0uZJWlWHoh/3YIlItNbuVtR76kWPO3b5GH5T9RtJy9QpxWw044HPPIDMFHGLyw9N3i/ECQolI71jwRKRKM1dzWjwir8o6b8v/jfePPOmZm7DAa7e6/qdxd9BXlqe6LG51lxY4i0KpKJowYIlorDWrFkT8vnlzsu41HlJ9Hb+VP0nHL4gz8IBcvpmwTdx48wbRY/LseZw8n4aFwuWiMLavn378H9LWREnEAzgVeer+KDxA7mjReyuuXfhc7bPiR6XbcnGtMRpCiSiaCOqYBsaGuByudDZ2QmXywW32w2r1QpBEJCVlQWr1YqFCxcqlZWIJllJSQkqKyvh8rjQ0i1uVZvBwCB+/emvcbLlpELppCvOKsbK/JWix2UKmVwZhyZs3IJtaGiA3W7H22+/DZfLFXLFoCBcPbnv9XqHHzMYDLjrrrtw9913484771QgMhFNlpMnT6LOXYe2njZR43oHe7Hv432o6ahRKJl0i2csxjcLvil63KyUWchIzlAgEUWrsAXb1dWF7du349VXX0VBQQFWr16NwsJCFBUVISUlZdQxp06dwsmTJ+FwOPDoo48iOzsb27Ztw/LlyxX7BohIGcFgEL4Un+hy9fZ7sefYHjR2SluJR0l5qXn4zo3fET1P8MzkmVx2jkQbtWBdLhfuueceFBcX44033kBBQcGENlZQUICCggKUlpYCAOx2OzZt2oSvfOUrePrpp2ULTUTKCgaDuNBxAamzxF3Ic6XnCl746AVc6RW/wLrSZifPxv2fuR8xRnFTGaYnpSNLyFIoFUWzUf+MO3LkCH72s59h9+7dEy7X0ZSWluLw4cPo6NDOVGhENLZAMICa9hq4+9w4ePDghMc1dDbguQ+f02S5psWn4aFlDyEhJkHcuIQ0ZFuyFUpF0W7Ugl29ejXuuusuWV5AEAQ899xzsmyLiJTlD/hR3VYNb//V6yomuh5sdXs1/vOf/4nOgU4l40mSHJuMh5c9DEucuHtWrfFWSWvBEg3hVIlEBADwBXyobq8OWSj9wMHxC/bT5k+x59ge9Pn6lIwnSZwpDhuXbMSMpBmixlniLZibOpcLplNExizYffv2oaura8TjDQ0NePzxx3HLLbdg4cKFuOWWW/D000+P+lwi0r5B/yDOtZ1D94C4+YEdLgde/uRlTU19OMRkMGH9jesxxzJH1LiUuBSWK8lizILduXMnXC5XyGOnTp3CihUr4HA4sHz5ctx7770oLCzEoUOHcMcdd7BkiXRmwD+Ac23n0DvYK2rc2+ffxqunXtXU1IdDDDDgW0XfwsJ0cfflJ8UmIS8tT/RVxkSjGfM+2NFWyfjhD3+Iu+66a8R5Va/Xi3vuuQc/+tGP8POf/1zelESkiH5fP861ncOAf2DUr//i+V+MeCwYDOJ3Z36Hv9f/Xel4kn19/texbPYyUWMSYxKRn5bPciXZiP6XdOrUKdx9990jHhcEAU8++SQcDocswYhIWX2+PpxtOxu2XEcTDAZhd9o1Xa53zb0LX8r9kqgx8eZ45E/Lh8loUigVTUXjFuz15yHuuusu2O32sM+1WLi6BJHW9Q724uyVsxj0D475vEc3PTr83/6gH6+ceEWT8woP+eKcL4qeAjHOHIf50+bDbOTU7CSvcf9Fbdq0KWTmJq/XC5fLhR/84AcjDgXv2bMHt956q/wpiUg23QPdqG6vhj/gn/AYX8CHlz99GVUtVQomi0xxVjG+ccM3RI2JNcVi/rT5iDGJm3yCaCLGLNjVq1eP+vitt96K7OzQm69XrVqFrq4ubNmyRb50RCSrroEu1LTXiCrXAf8AXvr4JZxpO6Ngssgsm70M9xbcK2pMjCkG+dPyEWuKVSgVTXVjFuy2bdsmvKHnnnsOVqsVycnJEYciIvl5+704334egWBgwmPu+/Z92HNsD853nFcwWWRunHEjvl30bVG31ZiNZuSn5SPeHK9gMprqZDvpYLPZ5NoUEcnM3efGhY4Lo94ZEE7PYA/a89pR31GvYLLILJy+EN+5Sdzk/SajCfnT8kVPm0gklqzXozc0NODIkSNybpKIItTe2y66XLsGuvCf//xP/P5vv1cwWWTy0/Lxvc98D2bDxPcTjAYj8tLykBiTqGAyoqtkLdgdO3bge9/7npybJKIIXOm5gtqOWlHl6u5z4/kPn8elrksYGJz4LTyTKceSgweWPCBqZRyDwYB5afOQHMvTWDQ5ZL0uffPmzbjvvvvk3CQRSdTS3QKXxzX+E6/R1tuGXx79Jdp6xa0BO5myUrLw0NKHEGeKm/AYg8GAeanzIMQJCiYjCiVrwdpsNp6LJdKApq4mNHrFLXje0t2CXx79Jdz97uHHkhKT5I4WkYykDDx888Oiz5/mWnNhiec9+jS5RBVsV1cX3G43OjuvLkmVkpLCK4eJNOZS5yVc7rwseswvP/plyEo6AFC0qEjOaBGZnjgdj9z8CJJjxP2+ybHmIDVB3MLxRHIYt2AbGhpgt9vx9ttvj5j4f4jNZkNxcTEeeOABZGVlyR6SiCamwduA5q5mUWMuei5iz7E96BnsGfG1ugt1yJmbI1M66azxVnx/2fdFr+mabcnGtMRpCqUiGtuYBXvq1CmsWrUK2dnZWL169YhDwB6PBy6XC/X19Xj11Vdht9vxu9/9DjfccIPiwYkolJRyPd9xHmXHytDnH30t1+bWZtULVogV8B83/wfSEtJEjbNZbEhPSlcoFdH4xizYnTt3jrpyzrWGpkbcvHkzNm3ahB/+8Id4/fXX5U1JRGNyeVxo6W4RNebMlTN46eOXMBDQ5pXCAJAUk4RHlj2C9ERxRWmz2EQvsk4ktzFv03E4HKOunBPO3XffDafTGXEoIpo4KeVa1VKFvR/v1XS5xpvj8dDShzArZZaocSxX0ooxC7agoAAVFRUT3lhFRQWKi4sjDkVEEyOlXI9fPo7yT8rhC/jGfe6SzyyRGi0iscZYPLTkIWRbssd/8jWyhCyWK2nGmIeIf/azn+Gee+6B0+nEfffdh4Uo3V7BAAAgAElEQVQLF8JqtSIlJQWdnZ1wu91oaGhAfX398EVQb7zxxmRlJ5rS6j31aO1uFTXmw8YPcajqEIKY2MQT3d3dsMZapcSTzGw04/4l9yM3NVfUuCwhCzOTZyqUiki8MQu2sLAQ77zzDsrKyvDCCy/A6/WGTKg9NDuMIAgoKSnBvn37eB8s0SSQUq4OlwOvnXptwuUKAGfPncUtt9wiNp5kZqMZ37vpe1gwbYGocZlCJsuVNGfc23RsNhu2bduGbdu2weVyobOzEx6PBwBgsVhgs9lC1oslImVddF/ElZ4rosa873ofr596XVS5Tjaz0YwNN21AQXqBqHGZQiYykjMUSkUknaiJJrh3SqQuKeX6D9c/8NtTv9V8ua6/aT0K0wtFjWO5kpZJmuz/Jz/5yZifE5H8pJTr3y/+PaI919xccedBpTAbzVh/43oUpYubNWp2ymyWK2mapLmIr5/RqaGhQZYwRDS6Oncd2nrETcD/3xf/G78787uIXnfGDGWvyDUbzVh34zoUzRBfrmJv3yGabJL2YK+90ImIlCWlXP9W97eIyxUAPvzww4i3EY7JYMJ3Fn8Hi2YsEjVuVsoslivpgqzrwRKRvKSW65tn31QokTyGyvXGmTeKGjcrZRZmp8xWKBWRvGRdro6I5BEMBlHnrkN7b7uocX+p/Qv+cO4PCqWSx1C53pRxk6hxLFfSGxYskcZILdfDtYfxp3N/kjVLqlXeZd6MBiPLlaaMqDtE7HK54HA4It7O8uXLkZmZiczMTJw4cQInTpwY/jwzMxO7du0CACxZsmT4sZKSEgDA1q1bQ57b1NSEP//5zyGPHThwAABCHlu3bh0AYN26dSGPA8CBAwdCHvvzn/+MpqamkMe2bt0KACgpKRl+bMmSq1Pd7dq1K+S5/J60+T3t3LkTGQsz8NkvfhYlK0tQXV2N6upqlKwsGf4YyrRmzZrhx771/30Lfzr3J9RdqMOHH344/DE4MAh3hzvksZaWq1MrXvvYubPnAADnzp4LeXz+gvloaWkJeczd4cbgwGDIY3UX6gAAVSerhh/7+PjHAK5eBPnhhx/inx/8E5/++lMkdSZN6Hv6n4/+TwDA3uf24uYbbtbU+xSN//b4PYX/npYvXw4pDMGh6ZhEuP/++7Fv376wn8uprKwMwNXZogCgtLR0zOfb7Xbs3LkTXq8XgiDg1ltvxZYtWyZ8D29DQwPuuOMOvPvuu1zbliZdbUet6D3Xv9b+Fb8/93tF8pw7ew7zF8yPeDs8LEx6JrUXNH2I+Mc//jGKi4uH/5ooKytDZWXl8OejsVgsOHr06HDBEumFlMPCf6v7m2LlCgAd7o6It2EymPDdG7/LC5poytHsIWKv14uKioqQMi0tLcXOnTsnNJ7lSnoSrVcLm41mSeU6O2U2y5V0T7MFW1VVBYvFEvKYIAhwuVzwer0qpSKSn9RJJPRQrutuXCepXHmfK0UDSYeIk5OTx/xcDteX6xBBEODxeMbcQ732Iien04mNGzfKno9IDhfdF0WX698v/l2WSSQmQupKOsPTH3KGJprCJBXsc889N+bncigsLBxetWeI1+uF1+uFy+UKe9FSYWEhLBbLcAF7vV7s2LEDW7ZsGfFct9sNt9sd8lhTU5NM3wHR2CRN3F//D7xxZvLWXG5paRE9XeLQqjicuJ+mOtGHiMebd/j06dOSw1xv8+bNw1cRA1cPGwuCEHbvFri64s+1e7clJSXYu3fvqM/dvXs3cnNzQz5uu+022fIThSN1VZzXT7+uUKLR1dbWino+y5XoX0QX7IoVK3D//fejq6trxNeOHDmCVatWyRIMuHpRU3FxMRwOBxwOB4qLi+H1ekUvmycIApxO54jHH3/8cdTW1oZ8vPfee3LFJxpVvadedLk6XA789tRvFUokD7PRjPtvul90uWYJWSxXikqSDhG///77uOOOO/C///f/xp133il3phCFhf/6YR269Sbc+Ven04n169fj6NGjI7422hir1Qqr1RrymNms6TuXSOfqPfVo7W4VNeafjf/Ea6de0/R6rrHGWNy/5H4smLZA1LgsIQszk2cqlIpIXZKuIi4vL8ctt9yCRx99FD/4wQ9G3ZuVQ2VlZcjnFRUVePDBB8M+XxAEbN68OeSxoSuOuVg8qU1KuX7S9AkOOQ+pVq4L5o9fmPHmeDy87GGWK9F1JBWs1WrF888/j927dw/vzX7wwQdyZ8OLL744XJBerxeVlZUhMzk5nU5s2rRp+PPRSnTPnj0jSpdosrk8LtHl6mx14sDJAwgEAwqlGl9SUtLYX49JwveXfR9zU+eK2q7NYmO5UtSL6HhoSUkJiouL8dhjj2HDhg0oKCiQKxeAqxc52e12AFev+H3uuedCDvV6PB4cOXIkZNamlStXoqysDIIgoL6+HtnZ2eNOr0ikJJfHhZbuFlFjzrWdQ/kn5fAFfAqlmpjjHx8Pe6tOSmwKHln2iOgJIWwWG2YkKbuQO5EWRHzCURAElJeXw2634yc/+Ymsi7EXFxejuLh4zK9ff75VEATe90qaIaVcaztqse/jfaqX61is8VZ8f9n3RRdltiUb6UnpCqUi0hbRBVteXo6FCxeOeLy0tBQrV67Eq6++KkswIr1r8DaILleX14U9x/eg39+vUKrITU+Yju/f/H2kJaSJGsdypalGdMHeeuutYb8mCAIeeOCBiAIRRYNLnZfQ3NUsaszlzst44aMX0OfrUyiVeDPTQ8+TzkyciUdufgTWeGuYEaPLseZgWuI0OaMRad6oFzm9+uqrOHLkiCwv0NXVhccff1yWbRHpQVNXEy53XhY1prWnFb/66FfoHuxWKJU0OXNzhv87MyUTj97yqKhyNRgMmJs6l+VKU9KoBVtUVITHHnsMTz/9NBobGyVv/M9//jPuuOOOMWdeIoomrd2taPSK+5lp723H/zn6f+Ad0N4iFlUnqwAAc4Q5+I+b/wPJsROfd3yoXFMTUpWKR6Rpox4iLigowOHDh/HDH/4QK1asQEFBAb7yla+goKAAixYtCju5f0NDA06dOoX3338fr776KgRBwLZt23DXXXcp+k0QaUFbTxvqPfWixnj6PfjlR7+Eu889/pNV0N3TjbzUPDyw5AHEm+MnPM5oMGJe2jwIcVw2kqausOdgBUHA888/D6fTibfeeguHDh2Cy+Uavkr42ttlrl0+LhgMori4GD/96U+xevVqBaMTaUdHbwcuei6KGtM12IVfHf2V6GkTJ5PZbcaDSx9ErCl2wmNMRhPy0vJE7e0SRaNxL3IqLCxEYWEhtmzZApfLBZfLhYaGBni93uGVaLKzs2GxWGCz2WS/F5ZI6zx9HtS6axEMTny2pd7BXrxw9AU0dWt39aaiGUVoaWkRVa5moxl5aXlIih17ggqiqUDUVcQ2m41TDhJdo7O/E+c7zosq135/P/Yc24OGzrFXplLTkowlWLN4DR44MPG7AmJMMchPy0dCTIKCyYj0Q/RUif/2b/+Gl156adxl64iiXddAF2raa0SV62BgEHuP70Wdp065YBFanrkc31n8HZgMJhw4cGBCY2JMMZg/bT7Llegaogv2hhtuwPbt23HnnXfinnvuYdnSlNQz2IOa9hpR8wT7gj689PFLqG6vVjBZZL6Q/QWUFpYOX2tx4OD4BRtnjsOCaQtEXQRFNBWILtjnn38eZ86cwc9//nNkZWXhhRdeGC7bffv2sWwp6vUO9qK6rRr+gH/CYwLBAF755BWcvnJawWSRWZG7AqsWrhI13Wm8OR4Lpi1AnDlOwWRE+iR5LuKSkhKUlJQAAE6dOgWHw4FDhw5h165dOHXqlGwBibSk39eP6vZqUfMEB4NB/FfVf+HTlk8VTBaZu/Pvxp1zxa3tnBCTgPnT5sNs5BrKRKOJ+Cfj9OnTqKiowNtvvw2XyxV2MXQivRvwD+Bc2zkM+gdFjfvt6d/io0sfKZQqMgYY8I0bvoHb59w+6td/8fwvRn08KTYJ+Wn5MBlNSsYj0jVJBXvkyBFUVlaisrISHo8HgiCgpKQEu3fv5m06FJUG/YM413YOA/4BUePePv82/uH6h0KpImM2mnFf0X1YNmuZqHHJscnIS8tjuRKNQ3TBfvazn0VnZydSUlJQUlKC++67j6VKUc0X8OFc2zn0+8StcHPEdQQVNRUKpYpMrDEWG27agIXpI1fGutajmx5FZUXl8OdCnIB5afNgNIi+fINoyhFdsCUlJVi5cuWYq+oQRQt/wI/qtmrRK9ycbDmJ106/plCqyCTFJGHjko3IseaIGmeNt2Ju6lxZ13wmimaiC3bbtm0jHjt9+jSqqqpQVFQ06lqxRHoUCAZQ3V6NnsEeUePOd5zHrz/9tahbeCaLNd6Kh5c+jIzkDFHj0hLSkGPNYbkSiSC6YB9//HF0dnZi3759AICdO3di3759CAaDMBgMeP7553HnneKuRiTSmkAwgJr2GnQPiFs+7nLnZew9vheDAXEXQk2GjKQMPLzsYVHLza1dsxbTE6djjnWOgsmIopPoEykOh2P49hyXy4W9e/figQcewJkzZ3Dvvfdix44dsockmkzBYBAXOi6gs79T1LiO3g7sObYHvb5ehZJJl2PJEb2WKwA88dATLFciiUQXrNfrRVZWFoCrVxMbDAaUlpYCAD73uc/B5XLJm5BoEg2Vq6fPI2pc92A3XvjoBbj7tbfs3MLpC/HIzY8gKUbcBPwZyRn4H1/6HwqlIop+ogvWZrPh9Omrs9EcOnQINpttuHBdLhcXAyBdq3PXiV6bdcA/gLJjZWjuaVYolXRLZy3F/UvuR5xJ3ExLs1NmI1PIRHOz9r4nIr0QXbBPPvkktm/fjoULF+LUqVPYvHnz8NfsdjuvLibdqvfUo723XdQYf9CPlz95WZOT99+efTvWLloLs0HcpRY2iw2zUmYplIpo6pB0m84777yDU6dOoaCgIGSPdePGjVi5cqWsAYkmQ4O3Aa3draLGBINBHKo6hFNXtDc16Ffzv4oVc1eIHjfHOgfTE6cPf75o0SI5YxFNKZJmcgq3Luzq1asjDkQ02S51XkJzl/hDoX8890ccvXRUgUTSGQ1GrC5YjeVZy0WNMxgMyLHmIC0hLeTxysrKMCOIaDycjoWmtKauJlzuvCx63N/q/oa/1P1FgUTSxRhjsOGmDZLKdW7q3BHlCgBbt26VKx7RlMOCpSmrtbsVjd5G0eOOXT6G35/9vQKJpEswJ+ChZQ9h0Qxxh3SNBiPmpc4Le/vOwYMH5YhHNCVxnSmaktp62lDvqRc97syVM/hN1W8QRFCBVNIIsQIeWvYQMlMyRY0zGozIS8tDSlyKQsmIpjYWLE05Hb0dqHPXiR5X76lH+SflotaCVdr0xOl4eNnDmJ4wffwnX8NkNCE/LR9JseLujSWiiWPB0pTi6fOg1l0retyVnit48fiL6PeLW1FHSVkpWXho2UNIiRW3B2o2mpE/LR+JMYnjPvfYsWNS4xFNeTwHS1NG10AXLnRcQDAo7vBuz2APXjz2IroGuhRKJt7C6Qvx6C2Pii7XWFMsFkxfMKFyBYATJ05IiUdEYMHSFNE72Iua9hrRK9z4gj6Uf1KOlp4WhZKJ9znb57BxyUbRszPFm+OxYPoCxJvjJzxmw4YNYuMR0f/DQ8QU9Qb8A6hur4Y/4Bc99nXn66hur1YglXgGGPC1+V/Dl3O/LHpsUmwS8tLyYDbyR55osvCnjaKaL+BDdVs1Bv3il4/7a+1f8UHjBwqkEi/GGIM1i9bgpoybRI8V4gTMS5sHo4EHrIgmEwuWotbQmq59vj7RY6taqvCHc39QIJV4ybHJeOAzDyDHmiN6bKQLpT/zzDOSxhERC5ai1NCyc2IXTAeAhs4GvHLiFU3c6zojcQYeXPpgyPzAEx6bNAM2S2SrW61duzai8URTGY8ZUVS66Lkoek1XAPD0e7D32F5N3I4zL3UeHlv+mKRyzRQyIy5XAMjMFDd5BRH9C/dgKeo0ehvR1tMmetyAfwBlx8s0sWj6kowl+Paib4u+KMlgMCDbki2plIlIXixYiirNXc1o6moSPS4YDOLAiQNo8DYokEqcO+feia/kfUX0eVOjwYjc1Nyw8woT0eRiwVLUaO9tl1yQf6r+E060qDupgslgwr0F94peDQe4OvVhXloekmOTZc20YoX4NWWJ6CoWLEUFb79X0vzCAPBBwwd4t/ZdeQOJFG+Ox/ob1+OG6TeIHhtjikF+Wj4SYhJkz7V//37Zt0k0VfAiJ9K9nsEenG8/L3oKRACoaa/B66dfVyDVxFnjrdj02U2SyjXOHIcF0xYoUq4AsG7dOkW2SzQVcA+WdK3f14/qtmrRUyACQEt3i+qr42SlZGHj0o2wxFlEj02MSUReWh5iTDEKJLvq8OHDim2bKNqxYEm3Bv2DqG6vllSQ3YPdKDtehu5B8ffJyqUwvRDfvfG7oucUBoCUuBTMS50Hk9GkQDIikgMLlnTJH/Cjur0a/T7x96sGggG8/MnLaO1pVSDZxHw++/NYdcMqSdMXpiakIteaK3l2JiKaHCxY0p1gMIjzHefRO9grafwfz/1RtQn8zUYz/v2Gf8fnbJ+TND49KR3ZlmyZU4XX2Ng4aa9FFG14kRPpTq27Fp39nZLGftL0Cf5a91eZE01McmwyHl76sORynZUya1LLFQAOHDgwqa9HFE1YsKQr9Z56dPR2SBrb1NWE31T9RuZEEzM7eTZ+sPwHyEvLkzTeZrFhdspsmVON76mnnpr01ySKFjxETLpxufMyWrulnTft8/XhpY9fUmWO4UUzFmHt4rWSLmYyGAzIseYgLSFNgWREpCQWLOnClZ4ruNR5SdLYYDCI/zr5X2jpaZE51fikTnsIXJ36cF7aPAhxggLJiEhpLFjSPHefG/Weesnj3617d9KnQYxkgXTg6sVQeWl5SIpNkjmZOOXl5aq+PpGesWBJ07oGulDbUStplibg6kxNb1W/JXOqsaXFp+F7S76HrJQsSeNjTbHIn5aPeHO8zMnEW7x4sdoRiHSLBUua1TvYi5r2GkmzNAFXy/mVE69IHi9Fflo+1t20Dskx0ibdT4hJQH5avqKzM4mxdOlS3qpDJBELljRpwD+A6vZq+AN+SeODwSB+/emv4ekXv+i6VF+c80V8bcHXYDJIm10pOTYZeWl5nJ2JKEqwYElzfAEfqtuqMegflLyNt8+/jXPt52RMFV6MMQalhaVYNnuZ5G1Y4i2YmzpX0sxORKRNLFjSlEAwgJr2GvT5+iRv42zbWbx9/m0ZU4WXGp+K+z9zP7IEaedbAWBa4jTMsczR5NSHa9asUTsCkW6xYEkzgsEgLnRcQPeA9An4Pf0evHLiFQQh7aIoMfLS8rD+xvURLXI+M3lmROWstO3bt6sdgUi3eDyKNKPB2wBPX2TnTH9z8jfoGuiSKdHoDDDgyzlfxiPLHomoXG0Wm6bLFQBKSkrUjkCkW9yDJU1o6W5BS3dkE0H8s/GfONN2RqZEo0uOTca3i76NgvQCydswGozITc2FNd4qYzJlnDx5Uu0IRLrFgiXVefo8aPA2RLaNfg/ePPumTIlGl5eah7WL10ZUjDGmGOSl5SExJlHGZESkRSxYUlXvYC9q3dInkhjy2qnX0DPYI1OqUAYYcOfcO1GSVxLRVb7x5njkT8tHrClWxnTKmjlzptoRiHSLBUuq8QV8qGmvkXyv65Djl4+jqqVKplShhFgBaxevxfxp8yPaTkpcCualztPdPa7Hjx9XOwKRbvEiJ1LF0O04A/6BiLbTNdCFN868IVOqUAvSFmDL57ZEXK7TEqchPy1fd+UKALt27VI7ApFusWBJFRfdFyO6HWfIb0//Vvarho0GI+7OvxsPL3sYKbEpEW1rdsps5FhzNHmP60Q8++yzakcg0i0eIqZJd6nzEtp72yPezsmWk/i46WMZEv2LNd6K7y7+Luamzo1oO1zHlYhYsDSp2nvbcbnzcsTb6R3sxeunXpch0b8Uphfi24u+jaSYyJaIMxlNmJc6Dylxke39EpG+sWBp0nQNdKHOXSfLtt4886ZsE/mbjWZ8Nf+ruH3O7REfyo0zxyEvLU8TS83JoaKiQu0IRLrFgqVJ0e/rx/n28xHfjgNcPcT84aUPZUgF2AQbSotKJa/deq2k2CTkpeXBbOSPFRHpoGDLysoAAIIgAABKS0vHHWO322GxWAAAHo9nQmNIOf6AHzXtNfAFfLJs793adyPeRqwxFivzVuL2nNtlWcHGGm9Fbmpu1K2Gs3LlSq4HSySRpn8b/PjHP4bNZsPGjRtRWloKr9eLysrKMccMlWtJSQlKSkpgsVhgt9snKTFdb2gC/0hWx7lWW29bxBc2zU+bj62f34ov5X5JlkKcmTwT89LmRV25ElFkNPsbwev1oqKiImSy8dLSUuzcuXPMcXa7PWRMSUkJC1ZFLq8L3n6vbNv7a+1fEQgGJI1NjEnEt4q+he/f/H1MT5gecRaDwYBsS7bmJ+wnInVo9hBxVVXV8GHeIYIgwOVywev1Dh8yvpbX64XL5RrxuNPpDDuGlNPS3YLW7lbZttc10IUPG6Wde71p5k24p+CeiO9rHWI0GDE3dS4s8Zbxn6xjTzzxhNoRiHRLswV7fbkOEQQBHo9n1LIcrZSHxlRVVaG4uDjkcbfbDbfbHfJYU1NTBKlpiKfPA5dn5B87kTjaeBSDgUFRY6xxVnyz4JsomlEkW46pNGH/k08+qXYEIt3SbMEWFhbC4wm9DcPr9Q7vpdpsthFjwu2lhivr3bt346c//WnIY2azGXPnRjbJwFTXO9iLCx0XZN9udXv1hJ9rNBhRnFWMr87/qqy3zCTEJCAvLU9XE/ZHYsmSJZyPmEgizZ6DBYDNmzcPX0UMXN1DFQQhbGECV0t2oh5//HHU1taGfLz33nsRZZ7qBv2DqGmvkXyeNJxgMIhad+24zzMbzVieuRz/6/P/C98s+Kas5SrECVgwbcGUKVcAaG5uVjsCkW5pdg8WuHpRk9PphMPhAAAUFxfD6/WOuvcKIOw51uv3hIdYrVZYraFre5rNmv5fommBYADnO85HPIH/aAwGw5iTQMQYY3Br1q34cu6XFVnIfHridGRbsnU7pzARTT7Nt0lhYeHwfw8dAg5XpEVFRaOWqdfrRVGRfOfgaHR17jpZJvAPJz81H5+2fBry2PSE6bgp4ybcnnO7bBcwXW92ymzMSpmlyLa1btGiRWpHINItTRdsZWVlyC03FRUVePDBB8M+f+jw8fXnYm02G68gVtilzkvo6O1Q9DW+e9N34XA54On3ICU2BQumLUBGcoZir8cJ+zHufedEFJ6mz8G++OKLw+dUhyaZuHZWJqfTiU2bNoWMKS0tDZk/1W63Y+PGjZMTeIpq62mTZQL/8ZgMJtyWfdvwvMFKlqvZaMb8afOndLkCwNatW9WOQKRbmt6D3bx58/AkEW63G88991zInqjH48GRI0dC9lg3btwIu90Oh8MxXM6cKlE5XQNduOi5qHYMWUXbhP2ROHjwILZv3652DCJd0nTBFhcXj7h39fqvHz16dMTjLNTJIecE/lrBCfuJSC78LUKSBIIBXOi4INsE/lqQlpCGOdY5nFOYiGTBgiVJ6j316BnsUTuGbDKFTEXP6erVsWPH1I5ApFv8U51Ea+1uRVtPm9oxZGEympCXlsdyDePEiRNqRyDSLRYsidI10AWXV945htUSZ47DDdNviPoJ+yOxYcMGtSMQ6RYPEdOEDfoHcaHjQlRc1JQSl4K5qXN5MRMRKYa/XWhChhZOH/SLW81Gi2YkzUCWkMVpD4lIUSxYmpDGzkZ0DXSpHSMiQwukT0+MfLH1qeKZZ55ROwKRbrFgaVztve1o7tL3qipmoxnz0uYhOTZZ7Si6snbtWrUjEOkWL3KiMfUO9uKiW98zNSXEJGBh+kKWqwSZmZlqRyDSLe7BUlj+gB/nO87LvrbrZLLGW5GbmsvJI4ho0rFgKaw6dx36ff1qx5BsVsoszE6ZrXYMIpqiWLA0quauZrj73GrHkMRkNCHHmqPIwutTzYoVK9SOQKRbLFgaoWugC42djWrHkCTeHI95afO4Eo5M9u/fr3YEIt3iiSkK4Qv4dDuZhCXeghum38ByldG6devUjkCkW9yDpRB6nUyC51uVcfjwYbUjEOkWC5aGXe68jM7+TrVjiGI0GJGbmsvzrUSkOSxYAnD1vOvlrstqxxAlzhyHeanzkBCToHYUIqIRWLAEf8CP2o5aXZ13tcRbkGvNhcloUjtKVGts1OfFbkRawIucCBc9FzHgH1A7xoRlJGcgLy2P5ToJDhw4oHYEIt1iwU5xV3quoKO3Q+0YE2I0GDE3dS4yBU7fN1meeuoptSMQ6RYPEU9hfb4+uDz6WDyd51uJSG9YsFNUMBhEbUetLuYZTktIQ7Ylm4eEiUhXWLBTVGNnI3oGe9SOMSajwYgsIQvpSelqR5myysvL1Y5ApFss2CnI2+/V/PquceY4zE2di8SYRLWjTGmLFy9WOwKRbvEipynGF/Chzl2ndowxpSWkYeH0hSxXDVi6dKnaEYh0i3uwU4zL49LsVIg8JExE0YQFO4W4+9xo721XO8ao4s3xmJs6l1cJE1HUYMFOEf6AH/WeerVjjIpXCWvXmjVr1I5ApFss2CmiwduguUPDRoMRNosN0xOnqx2Fwti+fbvaEYh0ixc5TQHefi+u9FxRO0aIeHM8bph+A8tV40pKStSOQKRb3IONcoFgABfdF9WOESItIQ1zrHNgNPDvO607efKk2hGIdIsFG+UavA2amcifh4SJaCphwUaxroEutHa3qh0DAK8S1quZM2eqHYFIt1iwUUwrVw2nJ6UjS8jiIWEdOn78uNoRiHSLv/Gi1JWeK+gd7FU1Q4wpBnlpeci2ZLNcdWrXrl1qRyDSLf7Wi0L+gB+N3kZVM1jjrShIL4Al3l0D4eUAAA6/SURBVKJqDorMs88+q3YEIt3iIeIo1NTVBF/Ap8pr80ImIqKrWLBRpt/Xj+ZudVbKSY5NRo41B3HmOFVen4hIS1iwUaaxsxHBYHBSX9NgMGB2ymzMTJoJg8Ewqa9NyqqoqFA7ApFusWCjSNdAFzp6Oyb1NePN8chNzeXSckRE1+FFTlFkMhdRNxgMmJUyCwXpBSzXKLZy5Uq1IxDpFvdgo0S/rx/uPvekvFZSbBLmWOZw0ggiojGwYKPEZFzYZDQYkSlkYkbSDMVfi4hI71iwUcAX8KGtp03R17DEW5BtyUasKVbR1yFteeKJJ9SOQKRbLNgocKXnCgLBgCLbNhvNsFlsSEtIU2T7pG1PPvmk2hGIdIsXOUUBJfZeDQYD0pPSUTijkOU6hS1ZskTtCES6xT1YnRv0D6LP1yfrNoU4ATaLDfHmeFm3S/rT3KzOpCVE0YAFq3NdA12ybSveHI8sIYvzBxMRyYAFq3OdA50Rb8NsNGN2ymxMT5zOmZgoxKJFi9SOQKRbLFidi+TwsMloQnpiOjKSM2AymmRMRdGisrJS7QhEusWLnHTObBT/N1KMKQZZQhYWzViETCGT5Uphbd26Ve0IRLrFgtU5MRcixZvjkWPNwaIZizAzeSaLlcZ18OBBtSMQ6RYPEetcemI6rvRcwaB/cNSvx5piYY23whJvgRAnTHI6IqKpiwWrczGmGCycvhCtPa3oHuhGEEGYjWYkxiTCEmfhfMFERCphwUaBGFMMZqfMVjsGRaFjx46pHYFIt3gOlojCOnHihNoRiHSLBUtEYW3YsEHtCES6xYIlIiJSAAuWiIhIASxYIgrrmWeeUTsCkW6xYIkorLVr16odgUi3WLBEFFZmZqbaEYh0iwVLRESkABYsERGRAliwRBTWihUr1I5ApFtRV7AulwsOh0PtGERRYf/+/WpHINItzc9FbLfbAQBerxcAsHHjxjGf73A4sHPnTni9XgiCgFtvvRVbtmyBzWZTPCtRtFm3bh1LlkgiTRdsWVlZSKF6vV7s2LEDW7ZsCTvGYrHg6NGjwwVLRNIdPnxY7QhEuqXpQ8Qulyvkc0EQ0NnZOaGxLFciIlKTpgu2qqpqRMkSERHpgaYLdvPmzVi1atXwedjKykps3rx53HEOh2P4o6ysTOmYRFGrsbFR7QhEuqXpc7DFxcV4+eWXsWrVKpSVlaG8vHzcQ7+FhYWwWCzDzxvrvK3b7Ybb7Q55rKmpSb5vgEjnDhw4wOkSiSTS9B7s0C03R48eRXFxMVasWDHuLTg2my2khEtKSrB3795Rn7t7927k5uaGfNx2222yfg9EevbUU0+pHYFItzRdsENXEQuCgG3btqG8vByPPfbY8C07EyUIApxO54jHH3/8cdTW1oZ8vPfee3LFJyKiKUyzBetwOFBcXBzyWHFxMR588MGwe7FOpxM333zzqF8b7dCy1WpFTk5OyEdWVlbk4YmIaMrT7DnYcHuphYWFYccIgjDiIqih7XCiCSLxysvL1Y5ApFua3YMtKSnBW2+9NeLx999/f3jP1ul0YtOmTcNfG61E9+zZM6Erj4lopMWLF6sdgUi3NLsHCwBbtmzBjh07YLVaIQgCvF4vHnrooeGvezweHDlyJGTWppUrV6KsrAyCIKC+vh7Z2dkoLS1V61sg0rWlS5fyVh0iiTRdsDabbcxpEYuLi3H06NGQxwRBGHe+YiIiIqVp9hAxERGRnrFgiSisNWvWqB2BSLdYsEQU1vbt29WOQKRbLFgiCqukpETtCES6xYIlorBOnjypdgQi3WLBEhERKYAFS0RhzZw5U+0IRLrFgiWisI4fP652BCLdYsESUVi7du1SOwKRbrFgiSisZ599Vu0IRLrFgiUiIlIAC5aIiEgBLFgiCquiokLtCES6xYIlIiJSAAuWiMJauXKl2hGIdIsFS0REpABNL7iuBr/fDwBoampSOQmR+gwGAxoaGtSOQaSqoT4Y6oeJYsFep7W1FQDXwSQCgKSkJNxxxx1qxyDShNbWVsyZM2fCzzcEg8Gggnl0p6+vD1VVVUhPT4fJZFItR0NDA2677Ta89957yMrKUi0HjY7vj/bxPdI+vbxHfr8fra2tKCoqQnx8/ITHcQ/2OvHx8Vi2bJnaMeDz+eDz+ZCRkaHpf3hTFd8f7eN7pH16eo/E7LkO4UVORERECmDBEhERKYAFS0REpADT008//bTaIWh08fHx+OIXvyjqpDpNHr4/2sf3SPui+T3iVcREREQK4CFiIiIiBbBgiRTicrngcDjUjkGkW3r/GeJ9sJPEbrfDYrEAADweD0pLS8d8/qpVq+B0Okc8XlxcjPLycknbpLHJ/R45HA7s3LkTXq8XgiDg1ltvxZYtW2Cz2RTJPxVI+Tdvt9sBAF6vFwCwcePGiLdJo5P7/dH9z1CQFHfo0KFgRUXF8OcVFRXBQ4cOjTlm+/btwaqqqmB9ff3wR0VFRbCqqkryNik8Jd6joe15PB7lgk8hUt6jF198MeRzj8cT3L59e0TbpNEp8f7o/WeIBTsJ/v3f/31Cjw3xeDzDv6Svde0/VrHbpLEp8R5d+8uGIifl3/yPfvSjMR/jz5F8lHh/9P4zxHOwCvN6vXC5XCMedzqdw4dEricIAgoLC0Mes9vtw4dbpGyTwlPiPSJ5Sf03X1VVNeq4SLZJIynx/kQDFqzCqqqqhs9JXEsQBFRVVU1oG16vN+ScgxzbpH9R4j0a4nA4hj/KysoizjpVSX2PNm/ejFWrVg2f56usrMTmzZsj2iaNpMT7M0TPP0MsWIUNnZy/3mj/GMOx2+0oLi6WdZv0L0q8RwBQWFiIoqIiFBcXo7i4GDabDTt27Ig471Qk9T0qLi7Gyy+/jB//+MdYsWIFCgsLh7fDnyP5KPH+APr/GWLBToJIDjd5vV643W5Zt0kjKfEe2Wy2kF8WJSUl2Lt3r+TXmeqkvEdDt3kcPXoUxcXFWLFiRchtH/w5ko8S74/ef4ZYsAob7a864Ool7BNht9uRnZ0t6zYplBLv0VivNdqtPTQ2qe9RWVkZNm7cCEEQsG3bNpSXl+Oxxx4Lu8c1kW3SSEq8P2O9ll5+hliwCisqKhr1H5nX60VRUdG44x0Ox4jDLJFuk0Ip8R45nU7cfPPNoz4/3C8jCk/Ke+RwOEYcti8uLsaDDz4Ih8PBnyMZKfH+RMPPEAtWYYIgwGKxjPiL7PpDH+E4HI4RF89Euk0KpdR7dP3FGkPb181N8hoi5T0Ktxc0dJ6PP0fyUer90fvPEAt2EpSWlqKiomL4c7vdHjJbidPpxKZNm8KOH+0f6HjbJHHkfo9G+wWwZ8+eEb8waOLEvkclJSV46623Rmzn/fffH95z4s+RfOR+f6LhZ4ir6UwSu90Om80Gr9c7Ygoxh8OBxx57DO++++6IX9QrVqzAG2+8MWrJjrVNEk/u98jr9cJut0MQBNTX1yM7O5vvUYTEvkcul+v/b+/+eRJm4jiAf5/kmU/ewLk2sToy4GwCON7CKgOytYtuwiBusODq4AoOHbFNcK1vQN4A9QXA8QJ8BtOG8l+lSnm+n4QBCvQaQr753V3v0Ol0kMlkIISA1hqlUin2W/F/tD3b/n3S/h9iwBIRESWAXcREREQJYMASERElgAFLRESUAAYsERFRAhiwRCk3GAyixdLTqtvt7vWuKvT/xIAlSjGtNS4uLlJz4/0yBwcHUEpxbWDaKwxYohTQWsOyLGSzWWSz2ajau7m5QbFYnFtyLimDwQBKKRiGgXK5vNGasPV6HYZhzD2mtx4rFArI5XJotVpJNp/oV/371w0govVs28Z4PEaj0QDwWfENBgN4nod+v/8rbQiCAEop5PN5XF5eotfrQSmFfr+/toKWUqLdbs+9Nq1arUIphUqlkvqKnAhgwBLtPK01fN+H4zgwTTN6vdfrwTTNXwujh4cHSClxf38P4LPqPDs7Q6fTwfX19crPCiFibV/ENE2YprnR9xGlAbuIiXZc2B08uxSj53koFou/1o7n5+e5ZepKpRI8z9vaOYrF4la/j+gvMWCJUigIAgRBsLYq3Ob5tNZz55NSRse2wTTN6NqI0o4BS7TDms0mlFIAPjcVMAwjFkCze236vr9wQlH4ue8KA3R239uwe3rdxtrhJC3DMKCUWjo5KrweBiztA47BEu2warWKk5MT2LaNdrsNKSWklPB9H8B8t/Hp6Skcx4mej8djlMvl2PZflmVhMplsdH4pJW5vb6MA/erenqEgCKKJTqsmR03vokKUdgxYoh0mhIhCSEoZddFqrZeG3XQ3rmVZEELEZvCGk5S+4zvBJ6VEPp+PTY7KZrNoNptL28IKlvYBA5YohUaj0dr3uK4Lz/PQbreXhvGmwq7h2a7gdZUtgIUbmOdyOby+vi58vxBio+sj2nUcgyXaQ1pr1Go15PN5FAqFH39fWEXPVpbh89mx2XUymczKajiTyXyxhUS7hxUsUQodHh6uDCjbtgEAd3d3c8e+MwYrhIAQAsPhMHZ8OBxCSvnlCjkIgpXjuT+tuIl2AQOWKIVWVYzdbhe+7+Px8XFhUH13DDa8R3V6EYinp6eV9+IGQQDbtmMTr8KFM2bvqZ3GlZxoH7CLmCiFlnXZaq1Rr9ejiU6+70ePn6pUKgiCAJZlwXVdWJYFrXVsjDU8Pt11rLWGUgqu68J1XSilIITA1dXV3DnC23e+2uVMtIsYsEQpFAbo7P2k4bZ1g8EA5XI59thkYf5VpJRwHAeTyQS2beP9/R2O48SqzSAI4HlebPUpx3FwfHyMVquFWq2Go6MjvLy8LKyu397eYtdHlGb/fHx8fPx1I4jo68rlcjRGui8sywLws1uJiHYFK1iilCoUCqnfaH2W53k4Pz//62YQbQUDliilSqUShBBwXfevm7IV3W4XUsqt3FZEtAsYsEQp1mg0UKvVUr+0oNYarVZrr7q7iTgGS5Ry9Xodo9Eo1eOWlmUhk8kwYGmvMGCJiIgSwC5iIiKiBDBgiYiIEsCAJSIiSgADloiIKAEMWCIiogQwYImIiBLwH1wsewsmjViVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set_style(\"ticks\")\n",
    "\n",
    "# Choose what to plot\n",
    "zi = 0; par0 = 'f'; par1 = 'sv';\n",
    "\n",
    "# Get the indices in the Covariance Matrix of the two parameters we want to plot\n",
    "p0i = varyzdep.index(par0); index0 = zdepind[:,zi] # change this to shapeind if using shape!\n",
    "p1i = varyzdep.index(par1); index1 = zdepind[:,zi] # change this to shapeind if using shape!\n",
    "\n",
    "# Get the fiducial values for the ellipse centre\n",
    "c0 = fids[par0][zi] # remove the [zi] if par0 is a shape parameter!\n",
    "c1 = fids[par1][zi] # remove the [zi] if par1 is a shape parameter!\n",
    "\n",
    "# Marginalise out everything else and make a new 2x2 Covariance Matrix\n",
    "marg_all = True\n",
    "if marg_all:\n",
    "    covmat_2D = np.zeros((2,2))\n",
    "    covmat_2D[0,0] = CovMat[index0[p0i],index0[p0i]]\n",
    "    covmat_2D[1,1] = CovMat[index1[p1i],index1[p1i]]\n",
    "    covmat_2D[0,1] = CovMat[index0[p0i],index1[p1i]]\n",
    "    covmat_2D[1,0] = covmat_2D[0,1]\n",
    "    \n",
    "    \n",
    "else:\n",
    "    fishmat_2D = np.zeros((2,2))\n",
    "    fishmat_2D[0,0] = Fishermat[index0[p0i],index0[p0i]]\n",
    "    fishmat_2D[1,1] = Fishermat[index1[p1i],index1[p1i]]\n",
    "    fishmat_2D[0,1] = Fishermat[index0[p0i],index1[p1i]]\n",
    "    fishmat_2D[1,0] = fishmat_2D[0,1]\n",
    "    covmat_2D = np.linalg.inv(fishmat_2D)\n",
    "    del fishmat_2D\n",
    "\n",
    "def get_ellipses(centre,CM_2D):\n",
    "    # We will need coefficients for the ellipse depending on the number of sigma level we want\n",
    "    CL_alpha = {\n",
    "        1: 1.52, \n",
    "        2: 2.48,\n",
    "        3: 3.44\n",
    "    }\n",
    "\n",
    "    # Get the errors\n",
    "    FM_2D = np.linalg.inv(CM_2D)\n",
    "    em0 = np.sqrt(CM_2D[0,0])\n",
    "    eu0 = np.sqrt(1./FM_2D[0,0])\n",
    "    em1 = np.sqrt(CM_2D[1,1])\n",
    "    eu1 = np.sqrt(1./FM_2D[1,1])\n",
    "\n",
    "    #print par0, eu0, em0\n",
    "    #print par1, eu1, em1\n",
    "\n",
    "    # Now calculate the ellipse parameters\n",
    "    a = np.sqrt(0.5*(CM_2D[0,0]+CM_2D[1,1]) + np.sqrt(0.25*(CM_2D[0,0]-CM_2D[1,1])**2 + CM_2D[0,1]**2))\n",
    "    b = np.sqrt(0.5*(CM_2D[0,0]+CM_2D[1,1]) - np.sqrt(0.25*(CM_2D[0,0]-CM_2D[1,1])**2 + CM_2D[0,1]**2))\n",
    "    theta = 180./np.pi * 0.5*np.arctan(2.*CM_2D[0,1] / (CM_2D[0,0]-CM_2D[1,1])) # should be in degrees\n",
    "\n",
    "    if CM_2D[0,0]<CM_2D[1,1]:\n",
    "        width=2.*b\n",
    "        height=2.*a\n",
    "    else:\n",
    "        width=2.*a\n",
    "        height=2.*b\n",
    "\n",
    "    # Create all the ellipses\n",
    "    ellipse1 = Ellipse([centre[0],centre[1]], width*CL_alpha[1], height*CL_alpha[1], theta, \n",
    "                       linewidth=0, fill=True, zorder=2, alpha=0.7, color='green')\n",
    "    ellipse2 = Ellipse([centre[0],centre[1]], width*CL_alpha[2], height*CL_alpha[2], theta, \n",
    "                       linewidth=0, fill=True, zorder=2, alpha=0.2, color='green')\n",
    "\n",
    "    return ellipse1, ellipse2, CL_alpha[1]*em0, CL_alpha[1]*em1, CL_alpha[1]*eu0, CL_alpha[1]*eu1\n",
    "\n",
    "# Make a test plot\n",
    "ell1, ell2, errm0, errm1, erru0, erru1 = get_ellipses([c0,c1],covmat_2D)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "\n",
    "ax.axhline(c1,ls='--',c='k',lw=1); ax.axvline(c0,ls='--',c='k',lw=1) # Centre lines\n",
    "#ax.axvline(c0-errm0,ls=':',c='g',lw=1); ax.axvline(c0+errm0,ls=':',c='g',lw=1) # Marginalised error lines par0\n",
    "#ax.axvline(c0-erru0,ls=':',c='c',lw=1); ax.axvline(c0+erru0,ls=':',c='c',lw=1) # Marginalised error lines par0\n",
    "#ax.axhline(c1-errm1,ls=':',c='g',lw=1); ax.axhline(c1+errm1,ls=':',c='g',lw=1) # Marginalised error lines par1\n",
    "#ax.axhline(c1-erru1,ls=':',c='c',lw=1); ax.axhline(c1+erru1,ls=':',c='c',lw=1) # Marginalised error lines par1\n",
    "\n",
    "ax.add_artist(ell1)\n",
    "ax.add_artist(ell2)\n",
    "ax.set_xlim(c0-errm0*2.1, c0+errm0*2.1)\n",
    "ax.set_ylim(c1-errm1*2.1, c1+errm1*2.1)\n",
    "ax.tick_params(color='black')\n",
    "ax.set_facecolor('none')\n",
    "#ax.axis('off')\n",
    "\n",
    "plt.xlabel(par0+'(z='+str(z_vals[zi])+')')\n",
    "plt.ylabel(par1+'(z='+str(z_vals[zi])+')')\n",
    "\n",
    "#plt.savefig('test.png', format=\"png\", bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdgAAAHACAYAAAABXHlvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Vl0E2eeP/xvafWmzTa2sS2vgMGyISSQxCJ7mMbQM+/7H7qD6cmckzAd6KuBnBnnvDcTupu7t02f0/RdQibMzfxfPJ2mu+c/CTZJmgSCBGHHlhfwLhnv1mpbe70XRgqySrZKLnn9fc7hHKyqeupB2Pq6qp7n9zAsy7IghBBCiKBEy90BQgghZC2igCWEEEKSgAKWEEIISQIKWEIIISQJKGAJIYSQJKCAJYQQQpKAApYQQghJAgpYQgghJAkoYAkhhJAkoIAlhBBCkoAClhBCCEkCyXJ3YKVxu91obW3Fhg0bIBaLl7s7hBBCllkgEMDY2BiqqqqQkpIS93EUsHO0trbi7bffXu5uEEIIWWH+8z//E7t27Yp7fwrYOTZs2ABg9o3My8tb5t4QQghZbsPDw3j77bfD+RAvCtg5QreF8/LyUFhYuMy9IYQQslLwfWxIg5wIIYSQJKCAJYQQQpKAApYQQghJAgpYQgghJAkoYAkhhJAkoIAlhBBCkoAClhBCCEkCClhCCCEkCShgCSGEkCSggCWEEEKSgAKWEEIISQIKWEIIISQJKGAJIYSQJKDVdJLg7/6/v1vuLoT9n5/9n+XugiAOHjyIwsJC/P73v1/urhBCSFwoYMmqUFdXB5VKtdzdIISQuFHAklWhrq5uubtACCG80DNYQgghJAnoCpZwOnLkCLRaLRQKBf7rv/4LKpUK9fX1qK2tDe/T0NCA5uZm2O127N+/H6dOnYqrXYVCAa1WGz722LFjOHr0aHh7ZWUlqqur8fHHH+PYsWOora0N9+fpc8x3/ljtEELIUqErWBJTY2MjioqKcOHCBVRWVuLDDz8Mbzt+/DiMRiPOnDmDCxcuwGw24+DBg3G129zcjKKiIpw5cwY1NTU4ffo0TCZTxPbTp0/j2LFj0Ov1nG3Ec/542iGEkGShK1gSk1arDT/7PHz4MJqbm+FwOGA2m9Hc3IybN29CqVQCAM6dO4fdu3ejsbFxweel+/btC+/z+9//Hnv37kVjY2P4CtRsNke0PZfJZIrr/Au1QwghyURXsCSmysrK8N9DI3jtdjtaW1uh1WqjgquqqiriSjReer0era2tEV/PF4rxnn+hdgghJJkoYElMarWa83WHwwG73Z6082q12nm3x3v+hdohhJBkooAlvOn1+vCt4qe1trZCp9Pxbs9gMKCqqmrZzk8IIclAAUt40+l00Ov1OHLkCAwGA0wmE44cOQKVShXXfNXm5mY0NTXBbDbj+PHjMJvN4VHES3F+QghZChSwJCHnzp3Dvn37cOLECbz77rtQKBS4cOFCXMfq9XoYDAYcPHgQbW1tOHfuHO/buYs5PyGELAWGZVl2uTuxklgsFrz55pv4+uuvUVhYuNzdWXNC82CppjAhZLVINBfoCpYQQghJAgpYQgghJAkoYAkhhJAkoEpOZEmdO3duubtACCFLgq5gCSGEkCSggCWEEEKSgAKWEEIISQIKWEIIISQJKGAJIYSQJKCAJeuGwWCAwWCAw+GIuT3WNkII4YsCliypI0eO4OTJk0t+3rNnz0Kr1UKv16OxsTFqJR6DwYCzZ8/S+rELaGxsxNmzZ6Pev6WwHN83sRgMBjQ0NCS0/jFZPyhgiSAOHjyI48ePJ/08TU1NMT9oGxsb0dTUhKamJjQ2NkZsMxgM4QUFQsvdha5oDQYDmpqaUF9fn9S+L9V7lCx79+6FSqXC0aNHl3StXZPJhIMHD6K1tVXwts1mM06ePImKigrs3bs3antjYyN2796NvXv34uzZs+HX9Xo9PvjgA5w+fTrqe42QECo0QQRRV1cHlUqVtPZNJhO++OILOJ1O2Gy2qO2NjY1QqVSora0FgHDI1tXVwWw2Q6FQhPfV6XQwGAzQ6/UAEN6e7LVkk/0eJZPJZILZbA6/v0tJp9Ohrq4uKVeLWq0Wp06dgkKhwCeffBLxfQHM/p85HA7U1dVx3t04evQoTpw4QcskEk4UsElw7v/5brm7EHbk/31pSc4z9wOmqakJH374IW7evClI+zqdDjqdDk1NTfjiiy+itjc2NkYsV1dbW4uDBw/G9cHX0NCwJKv70Idw4q5du4YDBw4kpW2TyYQDBw6gra0N58+fjwhYAFAqlTEfHahUKnpuT2KiW8Rk1XM4HJzPBE0mExwOB7RaLZxOZ/h1s9kcvsXZ1NSUtA9uIhyj0RgVfEJpbW0NXyU3NzdTYBLBUMASTsePH8eRI0fCX5tMJlRUVMBgMIRfa2howMGDBwFEDl46cuQITpw4AYfDgYqKClRUVEQFYENDQ/jZVlNT06L62traynnrValUhp/b6fX6cB9MJhO0Wm34OSzf255HjhzB8ePH0dDQgL1792L37t0Rz+eOHDmChoYGNDU14eDBg+F/H9cAr6fbeHpbrDbmE6utp88t5PvOdf7du3eH/89Df55+bxIR+n8LPStvaGhIyiCr2tpaKJXKiGeqJpMJVVVVgp+LrA8UsITTgQMHIsLUYDBAqVRGfCgbjUbU1NREHXvu3DmcOnUKSqUSnZ2d6OzsjBgU09jYiKKiIly4cAGVlZX48MMPF9VXh8PBeQvv6dA9evQoTCZTxDO206dPhwc2hbbF+8Hd3NyMoqIinDlzBjU1NTh9+nTEM8Lm5macPn0ax44di3nldfz4cRiNRpw5cwYXLlyA2WwO/8ISbxvxtgUs7n232+3zjrA+cuQI1Go1bt68iQsXLkCpVOKrr75CZ2cnjh49Gvd5uIS+D/V6PWpra3HgwAGcOHFiUW3GcujQoYiADV3dEpIICljCKfSBHvpwMxgMOHToUPhrh8MRfnbFZb7BPFqtFnV1ddBqtTh8+DAcDseib8vFc3xtbS30ej2USmX4Kjb099bW1vBVbjxt7du3D3V1ddDpdPj9738PrVYb8cFsNptx4cKF8FXRXCaTCc3NzfiP//gP6HQ6aLVanDt3DmazOdzOQm3waQtY3Pt+9uxZHDt2LOY2rVYbDlKdTof9+/cLdoV87do11NfXh9+D0P9ZiMlkSvgqee4V6uHDh2E2myO+z+dD07rIfNZUwD79g0EWR6lUQqfT4dq1awBmAzb04RO6tRrah6/Kysrw30NBbLfbF9VXLvO1+dFHH4UD4aOPPsL+/fsBzP5ikcj3kF6vj5hGEgryWFpbW8MB/7SqqqpweCzUBp+2gMTe96amJuzevRu1tbUxr0RPnz4dtc1mswkWPnOfv5pMpnDbodHiXCPL4zH3ClWr1UKn0+H8+fMAFg7Q0NzqhoaGhM5P1rZlD9hE5zVyMRgMOHHiBCoqKrB7924cP358WSbErxX79+9HW1sbTCZT+MpIq9WG545y3R6Oh1qtFrSfVVVVnEHhcDg4n5+dPXsWv/jFL8JfWywWwa9EFpon6nA4Fgy3eOeaxtMWkNj7Xltbi5s3b6KxsZHzZzB0J2BuX41GoyDPLkO/0D3dvsFgCLcduishpGPHjqG5uTli7vR8QncLDh48SJ83JMKyBazJZEJDQwMMBsOC8xpra2uhUqkWDFmVSoWbN2+G/4Ru3ZHEhK7mWltbw2G6b98+tLS0oLW1FXv27FnmHs5SKpWc0yW4rurMZjNsNpvgz9We/tCPR6jYxdwP5ESe+QnZVizHjh3jvA1rt9ujfsZC78XT53Y4HOEBSqF94r2tO7f9ixcvJrUoSOiW/MmTJ+MK74MHD2LPnj24cOECfd6QCMsWsDqdDh988EHMb+DGxsaI0Z21tbVxV0yh5yLC0Ol04VGVoWete/bsQXNzM0wmU/i2KpenR+k2NTUlvaRcXV0dLl68GP66sbGR85ZmQ0NDxNUrABQWFvJ+Btzc3IympiaYzebwnRI+g3l0Oh30ej2OHDkCg8EAk8mEI0eOQKVS8Z4vK2RbsWi1Ws6rs7l3DxwOB86ePYtTp05F7Bcard3W1gYAUbd8Y1W4mhtYDQ0NOHToUFy/OMzXLjB79yzWL0WHDh2K63PEZDLBZDLRHGfCaUUWmlhoXiMF6NKpqalBc3Nz+AMtNBCI6+rwaaHCECdOnEBVVVXUBy5fZrMZ58+fh9FohNlsRkNDA6qrq8O/hB09ehSNjY0RBfu5il8cOHAgqt+/+MUvcPHiRdTV1UVV8okldHX/4YcfQqVS4dy5c7yvXs6dO4eGhobwiNiampqIYhnL1RYfSqUS9fX14TrOAwMDOHPmTNR7XFtbC4PBEPEcOPQ+2+12GI3GmO99fX09GhoaoFarUVRUFHeYxWrX4XDg3/7t32A0GqHValFfXx913sOHDwv+KIOsPwzLsuxydiBUmefpSjoGgwEnT57EV199FbHv7t27cebMmZgfgE1NTRE/2CaTifcUAYvFgjfffBNff/01CgsLeR0bsh4rOa0GTU1NMee8mkym8O3OhYLyyJEjUCgUS1L9aaUI1QPu7OxMuI2zZ8+Gr7afHsUdMt//z3yamprQ0tKCDz74IOb2ZJV4FOJ9IStformwYq9gF5rXyEWn00GlUoWPdTgcaGhoiPmDZ7PZop7/Dg8PJ9hrstLN9yFLcx2Tr6WlJXz1abfbI97zp6trCSlZ7RISj2UfRRxLIvMi5/5GXFtbi08++STm/r/73e9QWloa8efll19OqL+ErGVCPJY5fPgwLl68yPnss6mpKaFfcgwGA7744gsYjUbOMRqJthuvhQpwkPVtRV7BJjKvcb62QtNM5nr//ffx7rvvRrw2PDyMd955h/d5nka3Zcla8/RauokO6In1aMfhcCQcgnq9Pintxuv8+fM4dOhQUs9BVq8VGbB85zUCs89C3n33Xc7VW2IFtlqtjhrIIJGsyLeErCDnzp1b7i4si3PnzoWn2hw+fFiwW69KpTIphfyT1S7ww3PfAwcOLMsSfmR1WJFp8vS8xqfDcb6Rq6HRjE8L3WamZzCECCM0L329o/chPizLYnx6HFa3Fb6ADwzDICs1CznpOWAYZrm7l3QrMmCBH+Y1hm5HzZ3XaDKZ8NFHH4VHcnKF6EcffZTUCemEEEIizfhmMOQawmPnY4y4RuANeKP2yUnPwWslr0EsEi9DD5fOsgXsYuc1hua4PX2Vu3///oj5eHzmzBFCCEnc6NQoOsc7MegcxEKzP0enRtE12YWK7Iol6t3yWPZ5sCuNEPNgCSFkPQiyQQzYB9Ax3gHrjJXXsWWaMrxQ+EKSeiasNTUPlhBCyMrlC/jwaPIRHk48xIxvJqE2ZvyJHbeaUMASQgiJS5ANomuyC6ZRE9x+96LaUsrX/vxhClhCCCELGrAP4MHIAzg9TkHay03PFaSdlYwClhBCSEwjrhHcG76HyZlJwdpkGAY56TmCtbdSUcASQgiJYnfbcXf4LoacQ4K3nZWaBalYKni7Kw0FLCGEkLBAMICW0RZ0jnciyAaTco7yzPKktLvSUMASQggBAFgcFtx+fBvTvumknSNNmoYSdUnS2l9JKGAJIWSdm/JO4fbQbQw6BpN+rq3ZWyFiVuxCboKigCWEkHUqyAbRMd6B1tFWBIKBpJ9PJpatm9vDAAUsIYSsS6NTo7g5eBMOD/+1t+diWYD1MECQgSgt9nPbLVlbIBGtn9hZP/9SQggh8AV8uDd8D12TXYtuK+hmEJiUIGCVgPXOro4jLfBCkuOP2lcuka/52sNzUcASQsg6MewaxveD32PKO7WodoIzDPyPZQg4olfDYX3cy9DtyN0BmVi2qPOuNhSwhBCyxgl11cp6GfiGpQhMSoAYy8SI1dHPcrPSstbVs9cQClhCCFnDhl3DuGG5saipN2wA8I9KERiVgA3GXiidkbEQpUc+g2UYBrvydyV87tWMApYQQtYgX8CHu8N30T3ZnXAbLAsEJiTwD0tj3vp9mjQ/enH1ck05MlMzE+7DakYBSwgha4wQV63BGQa+fjmCM/HNWRVlBCDWRN4elkvk2JG3I+E+rHYUsIQQskYEggHcGbqzqGetLAsERiTwj8gQd6VEBpAWRl+9rseBTU+jgCWEkDXAOmOFwWxY1LzWoPvJVes0v0pLkhwfRKmRo542Kjauy4FNT6OAJYSQVa59rB0PRh4kXJyfZYHAmAT+Iem8g5i4iNKDkGz0RbwmE8vwQsELCfVlLaGAJYSQVWrGNwOjxYgR10jCbQQ9DHwDMgRd0XNaF8KIWchKPGDmZPKu/F1IlaYm3Ke1ggKWEEJWIYvDgu8Hv4fH70m4Df+4BP5B/letIdIiLxhZ5K3hIlURitXFCfdpLaGAJYSQVUSQgUxeBj4zdyWmeEmy/VFFJVKlqdhdsDvhNtcaClhCCFklhBjI5J8Qwz8oAxtI7KoVAERpQUgKokcNP1/w/LoeNTwXBSwhhKwCHeMduD98P/GBTD7AZ5YjYE/8qhUAGBELabEHc5d0Lc8sR74if1FtrzUUsIQQsoLN+GZw3XIdw67hhNsI2MTwmWVg/YlftYZItV6IUiKfuyrlSjy78dlFt73WUMASQsgKtdiBTGwQ8D+Wwj8mFaQ/klwfxJmRz13FIjH2FO1ZV+u8xoveEUIIWWGEGMgU9DDw9fEvGhGLWOOHNN8X9fpzG5+DOkUtyDnWGgpYQghZQYQYyBSwieEbWNxApqeJMgKQFkUPaipRl6z7ak3zoYAlhJAV4tHEI9wZupP4QCaBbwkDgEgehKw0elCTUq6kKTkLoIAlhJBl5g/68f3g9+i39SfcButl4O0V7pYwADASFtJyD5g5SSEWifFS0Uv03HUB9O4QQsgysrvt+G7guxV1SxiYnY4jK/NAJGejtu3K3wVVikqwc61VFLCEELJM+mx9uDl4E/6gP6HjWfbJLeFR4W4Jh0iLvRClR9+qLtOUoUxTJvj51iIKWEIIWWJBNojbj28vutyh0LeEQ6T53qgyiACgSlFhV/4uwc+3VlHAEkLIEnJ5Xbg2cA2TM5MJtyFk4Yi5JNl+SHKjr6ilYileKnoJYtHiKkGtJxSwhBCyRAYdg7huuQ5vIHrKSzySeUsYAMTKACSF3H2rKayBUq5MynnXKgpYQghJMpZl8WDkAdrG2hJvw8vA2ydHcEr4W8IAIEoNQsqxtisAVOdWo0BZkJTzrmUUsIQQkkRuvxvXBq5hdGo04TYC9iejhJNwSxgAGOnsiGGG4+5vgbIAVTlVSTnvWkcBSwghSTI6NYprA9fg9rsTOj7Zt4SB0HQcd9TC6cBsMYmawpqknXuto4AlhJAk6BjvwL3he2DZ6OCKR7JvCQMAwwDSUg9EadF9lIlleKX4FUjFyQv3tY4ClhBCBOQP+nHDcgMD9oGE2wi6RPD2ycH6knNLOERa5IFYGT3XlWEY7CnaA4VckdTzr3UUsIQQIhCHx4HvBr6D3W1PuA3/mAT+QRkSvPCNm7TQG7X0XMjOvJ3Iy8hLbgfWAQpYQggRgMVhwXXLdfgC0Uu6xYMNAr4BGQLW5H8sSzf6INnAXT2qTFOGiuyKpPdhPaCAJYSQRRBiCo7Qa7fOR5LjgySP+5eA7LRsWiFHQBSwhBCSII/fA4PZgGHXcMJtBBwi+PrlSZuC8zRxlh/SAu5wTZOm4eXilyGauy4dSRgFLCGEJGByZhJX+69i2jedcBv+YQl8wzIgyc9bAUCsDkCq5a7SJBaJ8UrxK0iRpCS/I+sIBSwhhPDUPdmNW49vJb4wegDwDcgRsC1NXV+xIgBpMXeVJgB4sfBFaFI1S9KX9YQClhBC4hRkg7j1+Ba6J7sTb8PNwNcrR9C9NLdiRelBSEs9iHXntyqnCkWqoiXpy3pDAUsIIXGY9k3jav/Vxa2CYxfD1y/swujzEaUGZ6s0xbhQ1qq0qM6tXpK+rEcUsIQQsoCxqTF8N/BdwiUPAcA3LIV/aOmqIonkQcjK3WBifMpnpmbixcIXl6w/6xEFLCGEzKPH2oObgzcX97y1X46AfenWUWVkLGSbPGBi5HmGLAOvlrwKiYgiIJno3SWEEA4sy+Lu8F10jncm3MZSP28FnqyMs4m7eD8AyCVyvFbyGo0YXgIUsIQQMocv4MM18zUMOYcSbiPgEMHXJ+d83sqChd1jx+T0JKxuK3wBH+QSGbZkbUG6NCPhczKS2XAVybnDNTQdh2oMLw0KWEIIeYrT48SV/itweBwJt+EfkcA3FD2/1R/0YdD5GIMOC7xzSiq6/W60j7VjV35ilZQYMQtZuQeiFO5wZRgGe7R7kJ2WnVD7hD8KWEIIeWLENYLvBr6DN8BdkGEhbADwmaPrCXsCHpgdZgw5HyMQjP0sd8Y/k9B5GdGTcE2L3fau/F0oUBYk1D5JDAUsIYQAeDjxEHeG7iS8fmvQ8+R568wPz1tnfNPos/djdGo0rnbTJGm8z8uIWEjLPBClxw5XXY4OmzI38W6bLA4FLCFkXWNZFrce30LXZFfCbQScT563PqknHGSD6Lf3wewwIxiMP7DVPKspMQwgLfFCrIgdrqWaUmzP3c6rXSIMClhCyLrl8Xvw3cB3GJ0aTbgN/5gEvsEfnrdOzkzg4cSjhObM5qRviH9nBpAWeyBWca/pCgB5GXl4vuB53v0gwqCAJYSsS3a3HVf6r8DldSV0PMsCfrMM/onZj1FPwINHE48wPj2eUHuZqRooZMr4dmYAWbEHYk3scNWkavBS0Uu0Os4yooAlhKw7j52PYTAbEl8c3Q94e+UIusRgwcLisKDP1jvvAKaFlKhL4tuRAWRF84druiwdr5W8Bql46SpHkWgUsISQdaV9rB33R+4nPpjJzcDbkwLWw8DuseHhxCNMeacWPE7MiqGZ2QilJxsumRWj6f3hbZpUNZRy1cInD4VrZuxwlYllVEhihaCAJYSsC0E2iO8Hv0evtTfhNkLFI7w+P7qt3RhxDi+4lKuYFSPHVYqc6SIw7Ozt2iATGZDF8Vy9xhGuYpEYr5a8CqU8zlvNJKkoYAkha57b78bV/qsJPx8FAP+oBL7HMoxPjaFzohO+gH/e/RkwyJouxEZXGcTByFu1k6mPw3/PSsuCWq6e/+RxhCvDMNBr9VRIYgWhgCWErGnWGSuu9F/BtG86oeNDg5k84wy6Jjsw5Bxe8BiFJxNa51bI/NHzWj2SaTjkEwAAqViCiuyK+RtjAKl2/nAFgGc3PotCZeGCfSNLZ9kDtqmpCQaDAadOnYra1tjYCJVq9rmE3W5HXV3dgu0lcgwhZG0y280wWowIBOcPp1hCg5nsk1NoG2vDjG/+qTciVowC5xZkTceqmMSiX9US/qoiqwIykWzeNqVaDyRZ8/d/24Zt2JK1Zd59yNJbtoA1mUz44osv4HQ6YbPZoraHgrK2thbAbBA3NjbOG5iJHEMIWZtaR1vRMtKy8I4xBN0MPD0y9I9a0GfrW3BQVIZXgyJ7JWSB1Jj7DCm6MC11AgBy03ORnTb/vFdp0cLhWqwuxjN5z8y7D1keyzZBSqfT4YMPPoBer+fc3tjYGA5KAKitrUVjY+O8bSZyDCFkbQkEA/hu4LtFhWvAIYK9Dbjbfx+91t55w1XEilDorMCmyefmDVeXbBIjT0YOyyVybM6av3ShtMi7YLjmZuTSoukr2LLfIubicDhgNpujXjeZTHA4HFAqo0fIJXIMIWRtmfZN40r/FVhnrAm34R+VwPJwEo8mH8EfmD/g0n1KFNmrIOd41vq0gMiHPnUrAIABsDW7AhJR7DmqUq0Xkqz5B1GpU9R4uehlKiSxgq3IgG1tbQ0/R32aUqlEa2sr51VvIsfYbLao29PDwwsPYCCErDzj0+O42n81oRKFwOxgpuk+oONR14KlExkwyHOWI3eqJK62+1Um+EWzK/RsVOZDk5IZc1+p1gtJ9vzhmiZNo0ISq8CKDNhYV5xcAbqYY373u9/h17/+dcRrEokEZWVlPHpLCFluvdZefD/4PYJsYpWUWD8w1OZCZ3833H7PvPvKA2kotW5Hij++hdFH0/vgkM9OD0qVpmKTpjzmvvGEa6iQRKo09u1osjKsyIAFZgMz2ce8//77ePfddyNeGx4exjvvvMP73ISQpceyLO6P3Ef7WHvCbfhmgui4NYjB8aEFi0Zo3HnQ2rdBxIrjatueMorHitlVehgG2Ja9DSKG+1hp4cLhKmJEeLn4ZahS4qj6RJbdigzYWM9L7Xa7oMeo1Wqo1ZETvCWSFfmWEELm8AV8MFqMGHQMJtzG5JgLplt9cLnnnyPLsAy0zq3IjDn9Jtq01I5+VWv46yJVUcwKS9JCLyQb5g9XAKjR1iAnPSfuPpDltSLTpKqqijMYHQ4HqqqqBDuGELI6Tfum8W3ft7C5o6f4xavn0WN0tw8tOJBJHkhFqXVH3LeEAcArnkGP5h6CzOwta02qBiXqUs594w3X5/KfQ5GqKO4+kOW3IgNWqVRCpVJFPVfVarUxr1QTOYYQsvrY3DZ80/cNZnwzCR0fCARx904nxizOBffVuHOhtVfGfUsYAIKMHz2ae/CLZlfqSZHIUbmhEgyYqH3jDdcdeTuokMQqtGLHd9fV1eHixYvhrxsbG3H06NHw1yaTCcePH+d1DCFkdRt2DeOrnq8SDlfn9BS+/ebOguHKsAy0jm0otlXzCleARY/mPtyS2dV1RCIRqnKqIOWYkiMtiC9cKzdUonJDJY8+kJVi2a5gzWYzzp8/D6PRCLPZjIaGBlRXV4cLRRw9ehSNjY0wGAzhwUtPV2Sy2+0wGo0RV6wLHUMIWb26J7tx8/HNhJeZ6x8dRNvNAcA7f2Amcks4ZEDVBpfshzm4WzK3IEOmiNpPWuCFJGfhcN2StQU78nbw7gdZGRg20e/WNcpiseDNN9/E119/jcJCKpxNyErwYOQBTKOmhI5lweJW532MdsxAFJw/XBO5JRwynNGD4Yye8Nf5ynxsyYy+rRtvuJZpyvBC4Qu8+0GEl2gurMhnsIQQAsyu4XrDcgN9tr6Ejnf73fj29g0Eh1LmDU2GZVDorEDWdGK/VFtThyLCVZWixObM6FKI8YbaPGqvAAAgAElEQVSrVqXF8wXPJ9QXsnJQwBJCViRvwIur/VcXrKoUy4hrBNdvtEDmmH+QoywoR6l1B1J9iQ2GdMrHMaBqC38tl8hRlVMFZs4QF2l+fOGar8iHXqsHw0QPiiKrCwUsIWTFmfJO4Zu+b+Dw8C84AwCtQyY8vPcYqe75CzJkeDUotW2PWhA9XlMyG3rVD8A+KVEhFolQnVMF6Zwl6OIdLZyTnoOXil6i+sJrBAUsIWRFmZyZxLd93yZUUzjABvDNw6twPmSR6p8/XHOmipDvTHzqi1vqQo/mbniu62wR/61Rg5riKX8IAFlpWXil+BWIRfyf/5KViQKWELJiDDoGYTAb4A8uHEhzOb1OXGz5EqmWHMiDsRcxF7EiFNl1ULtzE+6nVzKNLs1tBJgfilQUqYuxIe2pKktMfIulA7Mr41Dx/rWHApYQsiI8mniE20O3E5qGY3aY8df73yF7omTewUyLmYIT4hd50KW5HS4kAQAb0rJRqi4Jf80wgLTYA7Fm4XBVyBV4vfR1yMSxfykgqxMFLCFk2d0duouO8Y6Ejr31+BZa27qR6yoD2NgDgxSeLJTaqiFiE//YC4h8eJR1G17xDyvuZMgysHXDNuBJpSZG9CRc1QuHa7osHW+UvoEUSUrCfSIrFwUsIWTZBIIBGC1GmO1m/seyAVzquoTJbg9yZ7jr/IbkThVjo3Nzot0EMFsCsUtzGx7xDwsDyCVyVOdWQ/xkhRxGxEJa6oFYufCyeSmSFLxe8jrSpPMv1k5WLwpYQsiy8Pg9uNJ/BePT47yPnfZN4y9t/w3xoAqZnvyY+zEsgyK7Dhp33mK6iiATQHfmHcxIXeHXJGIxduRuh1wsnz0Xj3CVS+R4o/QNKOTRVZ7I2kEBSwhZck6PE9/2fwunZ+GC+3ONTI3gf0xfIHOkeN5nqbPzW59Bqm9xIcYyQfRo7mFK+sOUIZGIQVVONdKk6QBmw1VW7oEoY+FwlYqleK3kNVrTdR2ggCWELKnx6XFc6b8Cj9+z8M5zPJx4iL+2f4v8ya2QBuQx90v3qlBq2wHJPKOJ48OiV30/or4wg9mF09Xy2bWkGfGTcE1fOFzFIjFeLX4VmamZi+wXWQ0WDNj29nZ8/vnnaGtrg9lsjqouolAoUFVVhdraWtTU1CSto4SQ1c9sN8NoMSIQXHgA0FwGiwEtjzqgtVfNO1I4a6YAWvtWgGN5OL561S1wyCciXivPLA9Px2EkT8I1beFwFTEivFL8Cjakb1h0v8jqEDNgm5ub0djYCJ1Ohz179uAXv/gFFIroWy1OpxM2mw1GoxHnz5+HWq3G0aNHqVA+ISRCx3gH7g7d5X2cP+jHxa6LGBtwoMilizlSmAGDQkfi9YTn6le3wJ4SWaaxUFWIQqV29nxSFrJyN0SpC08rYhgGe4r2IC9jcc+CyeoSFbBmsxmNjY348Y9/jE8//XTBBhQKBRQKBbRaLQ4dOgQA+OSTT6BSqfDWW28J32NCyKrCsixuD93Go4lHvI91ep34U9ufIBrJwMbp2KOAJUEpSm07kO5VL6arYf3qFlhTRiJey0nfgE2acgBPwnWTG6KU+Obsvlj4IgqVdNGx3kQErNlsRnt7O+rr6xfV6HvvvQen04lLly7hRz/60aLaIoSsXoFgANfM1zDoGOR97JBrCH9p+29oxrVQerJj7pfqU6DM+gykwdjPZPnoV7dGhas6RYWt2bNzXRnZk3CVxxeuuwt2o+SpIhRk/YgIWK1WC61WK0jDCoWCwpWQdcztd+Pbvm8xOTPJ+9i2sTZ8/fCvKLBum3cUsMadiyK7DgwrTHH8AXUrrCnDEa+lS9NRlVMFESMCI2ch3+QGI4svXHdu3IlNHMvWkfWBRhETQgTn8DjwTd83mPJO8TqOBYur/Vdxr78FJbbtkAZiVzjKd25CzlTJInv6gwGVCZNzwlUukWF73nZIRFKI5EHINnniDteqnCpszd4qWP/I6sMZsC6XCxkZ0fPLLBYLPv74Y1y/fh12ux2FhYXQ6/X413/916R3lBCyOoxOjeJq/1V4A15ex3mDXnz+8HMMD02g3LEToiD37/9iVowSWzUU89w25mtAZcJk6lDEaxKxGNufFJIQpQQh2+QGE2ct/orsClTnVgvWP7I6cd5X+c1vfoP29vaI19ra2nD8+HEolUq89957OHXqFOrq6jAwMIAf/ehHcLlcXE0RQtaRPlsfLvde5h2udo8d//vB/8bE4BSK7VUxw1UWSMGWieeTHq4iEYOqDVVIl2ZAlBqEbHP84VqeWY5nNz4rWP/I6sX5XVxXV4cTJ07g0qVL4deMRiMuXLgQte+hQ4dgMpnQ0NCAX//618nrKSFkRTONmvBg5AHv48wOM/7S+ReorBuROxV7DEi6T4ky6zMQL7p4xFPnVrVFhSsDYGvWVqhTNBClBSErd4OJ82FasboYu/N3C9Y/srpxftvodDoMDAxE3Cqeb16rTqfjnCNLCFn7WJbFzcc30T3ZzfvYlpEWfNXzNTbaNkPljl2AQeXOQYm9SrDBTMBsuE6kPo56vSyzDDnpuRClPwnXONc/L1AWoKawJqoYD1m/OAOWq2LTQt80Go1GuF4RQlaFQDAAg9kAi8PC67jQYKZb5tsotlUjzaeMuW/OVDHyF7kSzlyxwrVIpYVWWQRRRgCyMg+vcH2p6CUKVxKB89dBo9EIlmVhs9lgsVhgsVjAsmzELWMA4eeuZrOZrmAJWWd8AR++6fuGd7j6g378z8P/wd3+ByiffDZmuDIAtI6tSxau+cp8lGnKIVYGICvnH64iRrira7I2cF7B7t+/H1VVVXA4flg9QqvVRoToJ598gt/+9rc4deoUAISrOBFC1j63343LvZdhc9t4HTfjn8GfO/4M67gT5fadEAe5Rw6JWDFKbduh8GQJ0d0wi6qdM1zzMnKxJXMzxKoApCUexJuVFK5kPpwBq1AoUFlZOe+B7733Hvbt2ydYYQpCyOrg9DjxTd83cHn5zRywu+34rP0zBKwSlDi2x3yeKg3KUT65c96l6BJhUXZgPDW6otSGtGxUZG+FWB2cDdc47/JSuJKFLKrQBIUrIeuLdcaKb/q+gdvv5nXckGsIF9ovIN2WjXxXacz9Un0KlFt3CrDMXKRYt4Uz0zJRmVMJSVYAUq2XwpUISpBKTkajEQBouTpC1rAR1wiuDlyFL+DjddyjyUf4vPNzZDuKkTVdEHM/pScbpbbYV7aJYdGvNkWVPwRm6wtXbaiCNCcAaUH8/yYKVxIvQQL25MmTcDgcuHHjhhDNEUJWGLPdDIPZgCC78LqnT7s9dBvf9HyLQvvWeQv2Z09rUeioWGw352DRy7HkHAAo5ApU526HbGMA0o0UriQ5BAnYeJa1I4SsTl2TXbj1+BZYNr4avMDsNJzLvZdxb/A+SmzVSPOpYu5b6KhA9rSwj5tYhkWv+j4c8vGobemydGzP3Y6UwgAkOf6426RwJXwJErD0LJaQtal1tBUtIy28jglNw+kbM6PMthNyfxrnfiJWhBJbNZSe2AUmEsEyQfRo7sIps0ZtS5WmYkfuDqSVBCHJCsTdJoUrSQStpkMIicKyLO4M3cHDiYe8jpv2TeNCxwVYJ50ot8UerCQJylBu3TnvUnSJCDIBdGvuYEpmj9qWIpHjmY07kFHOQqymcCXJN2/AXrp0Cd999x2uX78Os9kcsU2pVKKmpgZ79uzBW2+9ldROEkKWTpANwmg2YsA+wOu4yZlJ/LHtj/A7GJTZn4GI5a7UkOLPQPnkTsEWSA8JMn50Z97BlNQRtU0mlmJH/g4oNgNiJYUrWRoxSyX+8pe/hE6nw549e7B//36oVD88Q7Hb7XA4HLDb7WhpacH58+dx5syZeesVE0JWPn/Qj6v9VzHsih51Ox+L04I/tf8JqU41ip2bY44EVngyUWrbDhEr7M2zgMiHLs1tzEij5+ZKxBLsyN8OdYUIooz4B2lRuJLF4vwuv3TpUtwDl0IVnH7729/SurCErGIevwff9H2DyZlJXsd1jHfgi64vkOUsRI6rJOZ+WTMF0Nq3YrYIonD8Ii+6Mm/DLYle3F0sEmNHfjUyKyUQpVG4kqXFGbALVXHiUlVVtejOEEKWx5R3Cpf7LsPpcfI67sbgDVztu4p852ZoZjbG3C/fuQk5UyWL7GU0n9iDrsxb8IhnoraJRCJUF1Qiu0oGUUr8I6ApXIlQOAPWYuFXvDvRYwghy8/utuNy32XM+KJDKpYggviq+ys8GG5BkV0Xs2YwwzIosVdD5c4RqrthXvEMujJvwyuOriolEjHQFVQgb3saGBmFK1kenN9FhYWFeP/992E0GsMr5szlcrlgsVhw6dIl/PznP0/oqpcQsrzGpsbwVc9XvMLVF/Thz+1/RutQG8qsO2KGqyQoxebJXUkJV49kGg+zbnKGK8Mw2FqwGQU7FBSuZFlxXsHW1NSgsLAQZ8+exYkTJ+B0ct82Co0k/tWvfkVzYQlZZQYdg7hmvoZAMP5RtW6/G39s/yPGrVaU23dC5k/l3C/Fn44y6zOQBbi3L4Zb4kJX5m34RdEVmBiGwdaichTtUMe93BxA4UqSI+ZQPq1Wi1OnTuHUqVMwm81wOp2w22fnlqlUqqjl6wghq0ePtQffD37PqzqT0+vEZ6bPMOXwoNwWe6m5DK8GZdYdgo8UBoAZqRNdmtsIiKIrMDEMg4qSMhRXZ1K4khUhrp8AujolZO1oH2vHveF7vI6ZnJnEH9r+ANYhRal9R8w5rpkzG1Fkr4TQI4UBYFpqR3fmHQSY6CtuhmFQUVaCkqqsuNdyBShcSXJRJSdC1pG7Q3fRMd7B65gh1xD+2P5HpDjUKHBuBlju8NzoKkfuPEvRLYZLZkWP5h6CMcJ1y6ZilOg2xL3cHEDhSpKPApaQdYBlWdwYvIFeay+v4/psffhz55+RaS/Ahqlizn0YlkGRXQeNO0+IrkZxyifQq76PIBM9j5VhGGzZokVZJb+BVBSuZClQwBKyxgWCAXw38B0eO6MXHJ9P+3g7Lj5qQp69HJoZ7vAUByUosz2DdK9aiK5GccjH0Kt5ABbRz4oZhsHmrYUo28ov2ClcyVKhgCVkDfMGvPi271uMT0cv2zaf20O38W33FWgd25DhyeTcRx5IRbl1J2QxVstZLFvKCPrVrbHDVZeP8s2xi1twoXAlS4kClpA1asY3g8t9l2F3R68sM5+rA1dxq/8OSu3bkRJjtZt0rwpltmdijiRerMnUxzCr2jiiFWAYYFN1HsrLC3i1SeFKlhqv77Tf/va3835NCFkZHB4HLnVf4hWuQQTR3N2MO333UW7dGTNcNe5cbLLuSlq4TqRZMBAjXCFisWlHHjaV85vZQOFKlgOvK9i2trZ5vyaELL+J6Ql82/8tPH5P3MeEFkkfHB5BuT32HNfcqWJsdG4WqqtRRtP78VjxiHMbKwpi8zN52FRcxKtNCleyXOgWMSFryLBrGFf7r8IfjC7EEIsn4MGfOv4E+8gMSh3bwcSY41rg3IINU/zCjY8hxSOMpPdzbmNFAWx+Ng+btSW82qRwJcuJApaQNcLisODawDUE2fiXZZvyTuGz9s8QGJOgyFXJOceVAYNiWxXU7lwhuxvBrGrDRCr3KOeg2I/Nz+VhcwG/ObZalRZ6rZ7ClSwbClhC1oA+Wx+uW67zKn1odVvxmekzyCfU2BhjjquIFaPMugMZXu6RxIvFMiz6VC2wp4xybvdLPKh4Lh9b8st4tVuiLsGLhS+C4VN5ghCBUcASssp1TXbh1uNbvMJ1ZGoEn5k+g2ayMOY6rpKgFOXWZ5EaY7DTYgWZAHo09+CSWTm3e6XT2LqrAFvyynm1W55Zjt35uylcybKjgCVkFesY78Ddobu8jhmwD+AvHf+N3IkyKDzZnPvMznF9NuZqOYsVEPnQpbmDGSn3Sl0zMid0u7WoyOE3oGpL1hY8l/+cEF0kZNEoYAlZpVpGWtA62srrmIcTD3GxswmF1m1I86o490n1KVBu3QlJUCZEN6P4xG50Zd6GR8y9Bu2U3Irq3aWo2LCFV7uVGyqxI2+HEF0kRBAUsISsQneG7qBzvJPXMa2jrfjq4V9RbK1Gij+dcx+FV4NS6zMxV8tZLI9kCl2Zd+ATcU8hcqSNYedzm7Elm1+4VudWoyqnSoguEiIYClhCVhGWZXHz8U10T3bzOu720G1ce3QdZbZnIA2kcO6jceei2FaFZCw1BzxZbk5zl3MtVwCwZQzhhV3VKNPwG9D0TN4z2LZhmxBdJERQFLCErBJBNojrluvot3HPFY3luuU6bnXfm7e04YYpLQqcFUJ0k5NTNolezT3OFXEAYEJpxiu7dqNIxW+e7a78XdiclbzCF4QsBq+AnTtKkc+oRUJI4gLBAK6Zr2HQMcjruCv9V2DqeTjvIun5zk3ImSoRoJfc5ivaD4bFmLoPbzy3B4WKwrjbZBgGzxc8z/tql5ClxCtgP/3003m/JoQIzx/040r/FYy4RuI+hgWLr3q+Ql/vEIqdVWDY6GILDBho7duQOZMvZHcjTKRZYFZyL/DOMkGMZvag9tk3kJcR/5JzDMOgprAGxWruubuErBR0i5iQFSyR5eaCCKLpURNG+u0odG7l3EfEilBq2x5zmo4QRjJ6MZTB/aw4KPJjNKsHP975I+Skx79YuogRQa/VQ6viV+yfkOVAAUvICuX2u3G59zJsblvcxwTYAP6n83/g6Pcjb5q7QIM4KEG5dSfSfNzTdIRgUXZiPM3M3UeRD6M53fi/d/wYWWlZcbcpFonxUtFLyFck74qbECFRwBKyAk37pnG59zIcHkfcx/iDfvy548/w9cuR7ea+wpMFUlBufRbyJC2SDrDoV5tgTRnm3OoTuzGe24u/3/F/QZOiibtViUiCV4pfQW5G8uohEyI0ClhCVhiX14W/9v4VU96puI/xBrz4k+nPYCwKqD3cV4Up/nRsmnwuaQUkWCaIXvV9OOQTnNvdkilY8/rx0+0HoZQr425XKpbi1eJXsSF9g1BdJWRJrLllJsxmMwwGw3J3g5CEODwOfNXzFa9wdfvd+Kzlj5AMaKCIEa7pPhW2TOxOWrgGRD50Zd6KGa7TMjscBWbU7XiLV7jKJXK8UfoGhStZlVb8FezZs2cBAErl7A9lXV3dvPsbDAacPn0aDocDSqUSNTU1+OCDD6DV0qAIsrLZ3DZc7r0Mt98d9zEz/hn88cGfkDG4EfIY1ZmUniyU2nZwjiQWgk/seVL6cJpzu1M+Dm/BJA7p3kKqJP7axqnSVLxe8jpUKcl7VkxIMsUdsH/4wx+wf/9+ZGRkJLM/EU6ePAm9Xo/a2loAs2Hb1NQU/pqLSqXCzZs3wwFLyGowOTOJy72X4Q144z7G5XXhT/f/G6qhwnmqM+Wh2KZDsqozuSVT6J6n9KE1dQiighkcqjwEmTj+q+d0WTreKH0DGbKl+7whRGhx/0q7b98+vPPOOxgcjJzo/sknn+Df//3fBe+Yw+HAxYsXI8K0rq4Op0+fjut4CleyWoxNjeGvvX/lFa4OjwMX7vw3NI+LY4Zr9rQ2qaUPp2Q2PMr6Pma4jqX3Q1bsw091P+UVrkq5EnvL9lK4klUv7oAN3W79+OOPYTQaw6+/9957ePDggeBFJ1pbW6FSRd4aUiqVMJvNcDjiH1lJyEo24hrBN33fwBfwxX2MzW3Dn29/jqzh0pilDze6ylHoSF7pQ7t8FF2a2wgwgeiNDIshxSOoimT4+61/D4ko/idRmlQN3ix7E2nSZI1yJmTpxP2d/8tf/hJFRUWor6/H6dOnoVarsW3bbIHtAwcO4OTJk/inf/onwTo2N1xDlEol7Hb7vFeoTw9yMplMOHr0KOd+NpsNNlvkHMPhYe7pBYQIbcg5hKsDVxEIcoRUDBMzE/ji1tfIniiN+UxV69iKrOn4yw7yNZE2CIuynavwIVgmCIuyA9qiXOzbtA8iHuMos9Oy8WrJq7yudglZyeL+7r948WJ4oFB9fT0+//zz8JWs2WyGTqcTtGM6nQ52uz3iNYfDAYfDAbOZewJ76Liqqiro9Xro9XpotVo0NDRw7vu73/0OpaWlEX9efvllQf8dhHCxOCy40n+FV7iOTY+h6fvLyJ4ojln6sMRWndRwHc7ogTlGuAZFfvSpW1BWUojaTbW8wjU3Ixevl75O4UrWlLh/At566y1cu3Yt/HV9fT2uXbsGo9EIu92OM2fOCN65+vr68ChiYPa2sVKpjHl1CwBarTbi6ra2thaffPIJ577vv/8+ent7I/5cvXpVuH8AIRz6bf34buA7BFnulWW4DLuG8dX1a8icLALY6GeqIlaEMuszULuTV4jBomrHcEYP5za/yIse9X1UlVZgb9leMDye+xYoC/Bq8au8biUTshrE/R39wQcfwGg0wuVyhUcSh24Xt7e3Q6FQCN65uro6mEym8C1fvV4Ph8PBe8qNUqmEyWSKuspWq9VQq9URr0kk9ENOkqfH2oPvB7/ntRLVkGsIV67fhtrJXSIw2aUPWYZFn6oF9pRRzu1e8Qz61A/wYvnzeKHgBV5tF6uLUVNYA4ZJzkAsQpYTr4lxNTU1UdN06uvrcejQoaSMJAZmb/mGbveGpt7Eev5qMpmwe/duzm00qpgst0cTj3DDcoNfuDqHYDA8gNLJXRBfGpRjy+TupIVrkPGjW3M7Zri6pU70Zt7DGxWv8Q7X8sxyCleypglyubZv3z4hmokyd87rxYsXcezYsZj7K5VK1NfXR7wWGnFMhSbIcuoY78Ddobu8jhlyDOG6oQ1pM9w1e+WBVGyafC7mNJ3F8os86Mq8C7fExbl9SmbDoKYDf7v1x9iUuYlX2xXZFXh247NCdJOQFWtFl0r8+OOPwwHpcDjQ1NQUUcnJZDLh+PHj4a+5QvSjjz6KCl1CllLraCv/cLUN4+a1h0iZ4X70kupTYMvE80kLV49kGg+zvo8ZrvaUMQxld+InVQd5h2tVThWFK1kXIq5gzWYzXC5XePrNYl26dAk/+tGPEj6+vr4ejY2NAGan1Jw5cybiVq/dbofRaIyo2rR//36cPXsWSqUSAwMDKCoqWrC8IiHJcn/4PtrG2ngdM2QdwW1DF6Re7rKCGV4NyqzPQMSKhehilGmpHd2auwiI/JzbJ9IG4cwcxmHdYWxI41cjeOfGndiazb1GLSFrDcPOeSBkNBrR1taGn//85wk32tbWBqPRiLq6uiUtrSgEi8WCN998E19//TUKC5M33YGsfXeG7qBzvJPXMUOTo7hr7ALj5X56o3LnoMReDYZjJLEQHPJx9KkfIMhwj3AezuhBIMuFt3RvQSWP/7kvwzDYlb+L99UuIStBorkQ9VNcU1ODwsJCnDx5EgzDoLa2FjU1NXF1oKmpCa2trairq1tUQBOy2iUUruOjuHe9B4yPO1yzZvKhtVcK0T1Ok6mPYVa1g+WY5TpbQKITKVkM6ir/gVelJREjQo22BkWqIiG7S8iKx/mTrNVqcerUKZjNZhiNRpw/fx5tbW1QKpXQarVQqVSw2+1wOBywWCxQKBTQarWoq6vDe++9t9T/BkJWlITCdXQM92/0An7uYREbprQocCav9OFIRi+GMro5twVFfgyo2pC1QYn/tfV/8SoGIRVL8XLRy7RQOlmX5h1FrNVqodVqcejQIQCA0+mE2WyG3W6HSqWCQqGAWq1OyhxYQlajhMJ1ZBwPvu8D6+e+7Zs7VYKNzmTdWmUxoGrDZOoQ51af2IN+VQuK8wrx4y0/hpiJ/7mvXCLHayWvITM1U6jOErKq8Jqmo1AoUFmZvFtUhKxmCYXr4wk8uNWLYIB7buxGVzlyXaVCdC9KkPGjR3MfLpmVc7tHMoU+dQt0Bdt4V2dKl6Xj9ZLXoZDTL99k/aKyRYQIIJFwHbZMouVO7HDNd25GzlSxEN2L4hXPoDvzbsxF0qdkNgyoTHih+Hns0e7h1bY6RY3XSl5DqjT+xdUJWYsWHbB/8zd/A71ejz17Zn8IFzMth5DVKKFwHbCi5V4vAgHu0bqFjq3ITlLR/mmpHT2ae/CLuJfIs6eMYVDZgTfK38DOvJ282t6QvgGvFL9CRfsJgQAB++WXX0Z83dbWBoZhBJtLS8hKllC49trR0tIDf4B7JZ0ieyUyZ7jrDi+WPWUU/arWmNNwxtMsGFP14W83/y0qsvgNqipQFmCPdg/EouTMzyVktRG8klNlZeW8y8kRslYkFK7dDrS0dHGGa2i5uWSF62h6H3rnmeM6pOjCpNqMn2z7Ce9wLdOU4eWilylcCXnKoq9gT548CaPRiKKiIuj1elRWVqKlpYVuFZM1jW+4siww8siJlo5HscPVWg2Vh7uo/+KwMKs6MJE6yL31ySLpPsUU6irrkJeex6v1yg2V2JG3Q4iOErKmLPoKds+ePfjyyy/xq1/9CizL4uLFizhw4IAQfSNkRUokXEc7p9DS8ZA7XFkGZdYdSQnX2dVw7sYM14DIj171A7AqD/6h+h94h+vOjTspXAmJYdFXsPv27QvXHKYiE2St4x2uQWC0fRoPujs4w3V2ofSdyPByr5izGD6RB93zrIbjE7vRp26BUpmOn1b+FBmy+MuaMgyDFwtfRIm6RKDeErL2xH0F+4c//CHmNrodTNYD3uEaACbaPTHDVcyKUW59NinhOiN1ojP7RsxwdUtc6NbcRV7WBhyuOswrXMUiMV4tfpXClZAFxB2wJ0+ehMvF/cNKyFrHO1z9gLXdj/vdbfAHolelEQcl2DS5C+letZDdBAA45GN4lHkTfpGXc7tLZkWP5h62btyCn1T+BCmS+Je8k4lleKP0DWxUbBSqu4SsWXHfIg4tuuNyuWCz2WilGbJu8A5XLwNbZxD3+kzwBqLnmkqCUmyafA4pfuFXmhpPN8OiiN1Xe8ooLMoO6Iv00Gv1vNpOk6bhtZLXoEqJfxUdQtYzXs9g3yvNTicAACAASURBVHnnHZhMJgCzz2D0ej0++OADbN1K6zuStYlvuAY9DFydDO4PtMDj90RtlwRl2Dz5HOT+dCG7CQCwKDsxnhZ7itx4uhljyn7UltWiKqeKV9tKuRKvl77OaxUdQtY7XqOIq6qq0NHRgY6ODly6dAkvvvgi/vmf/xmffvppsvpHyLLhHa5uBlOdItwbaMGMzx21XRZIwZbJ3YKHK8sE0aO5FztcGRZDikewqgfxk20/4R2uWWlZ2Fu2l8KVEJ54BezRo0fDf9dqtTh69Ci+/PJL9PX1wWg0Ct45QpYL73CdEmG6U4L75hZM+aaitssDqdg8uRsyv7D1ef0iLx5m3oRDPs65nWUCGFC2wadx4WdVP0Oxil9t442KjXij9A3IJXIhukvIuhJ3wFZWVsJisXBuO3XqFAwGg2CdImQ58Q3XgFOEmS4JHjxugdMbPRAwxZ+OzRO7IQ0IG1JuiQud2TcwI3Vy90vkQ6/6AVKyRHi7+m1sSNvAq/1idTFeKX4FEhGtCUJIIuIO2Pr6epw8eRKDg9wT1ouKigTrFCHLhXe42sXwdMvQOmSC3e2I2p7qy8DmyV2QBIUtfu+Uj+Nh1k34RNHPeYHZOa7dmrvIy83Gz6p/xmsaDgBUZFdAr9VDxAheTZWQdSPunx69Xo9/+Zd/wZtvvolPP/004mrW5XKhv78/KR0kZKnwDlerGJ5eGUwjbZiciV5TNc2nxObJXRAHpUJ2ExNpFvRo7iHIcC8W4JY60a25i8rCCvz9tr+HTMQv3Hfk7cCzG58VoquErGu87v3U1tbiyy+/xC9/+Uv85je/gUqlgkqlgkKhwJkzZ5LVR0KSjm+4+scl8Fmk6BjrwPh09PPPdK8K5dadELHC3l59rHiE0fTYv8y65JMwK9vxUukePF/wPK+2GYbB7vzdKM8sX2w3CSFIoFSiVqsNjxpua2uDQqGAVqsVvGOELBXe4Toige+xDA8nH2LENRK1XeHVoMy6Ewwr3O1VlgmiT9UKe8pozH2sqcMYVnXjwOb92JrNb+qcWCSGXqtHoZLmtxMilEX9el1ZWSlUPwhZFnzD1fdYCv+IFN3Wbjx2PI7arvBkocz2DBiWEayPfpEXvZp7mJJGP+MNGUvvh109jEPb3kKhgl9ISsVSvFL8CnLSk7GSDyHrFw0PJOsW73C1SOEfk6Lf3gezPXrOaTLC1SOZQrfmLrzi6Hm1wOyV7WPFI7BqN96ufBuZqZm82k+RpOD10tehThG+ZCMh6x0FLFmXeIerWQb/uAQWhxm91r6o7QpPpuDh6pJZ0au5h0CMwUwBkR9mVRsUWSk4uPVt3oUgFHIFXit5jfcIY0JIfChgybrDJ1xZFvANyBCYlGDYNYTuye6ofRTeTJTZdgoarpOpj2FWtYMFy7k9tNScNicff7fl73jPVdWkavBayWu8Cv0TQvihgCXrCu9w7ZchYJVgbHoUneOdUXE3O6BJ2CvXoYxujGT0xtw+I3WiX9WKHdpqvF7yOhjwO3duRi5eLnoZUrGw04cIIZEoYMm68WDkAb9w7ZMjYBPD7rGhfbw9RrgKN1qYZVgMqFphTYkemRzikI/jseYh9pa9iercat7nKFGX4IXCF6iABCFLgAKWrAttY20wjZri2pcNPglXuxjTvim0jLYiGIyM1wyBwzUg8qFHfQ9TMnvMfcbTLHBohvDWtp/yHikMALocHbbnbl9MNwkhPFDAkjXv0cQj3B++H9e+bBDw9cgRcIrhDXjwYORB1ILpGV4NygUMV49kCj2ae/CIZ7h3YFg8zuiCONuHf9z2j1DJ+a3HSgUkCFkeFLBkTeu19uLW41tx7csGAG+PHEGXGAHWjwcjLXDPWdNV6HB1yMfQp26NWfYwyAQwoGpD3sZM/HjLj3mXPZSKpdij3YONio1CdJcQwgMFLFmzzHYzbgzeiGtfNgB4u1MQnBKBRRCtoya45qyMk+FVCxquIxm9GMqIHpUc4hN70K9qxTOlVXip6CXeg5lSpal4reQ1muNKyDKhgCVr0pBzCAazASzLPc3laaz/SbhOiwCw6BjvhHVO8f7ZcH1WkHBlmSD6Va2wzVP20C11wqLpwN6KN7Atexvvc6hT1Hi15FVaJJ2QZUQBS9ac0alRXB24iiAbXHBf1g94u1IQnJkNzh5rb1R94XSvSrABTV7xDHo09+GWRK8bG+KUT2AyewBvVf4EeRl5vM+Rl5GHl4peomk4hCwzCliypkzOTOJK/xUEgtzPNJ/G+p6Eq3s2OAedgxiwD0Tsk+5Todz6LESseNF9c8ms6FXfR0Dkj7nPRNog2DwX/nHb2wlVWCrTlOH5gufBMMLNyyWEJIYClqwZdrcdl3svwxfwLbgv62Xg7ZIj6JkN1/HpMXRNPIrYJ92nRPmkMOE6nm7GoOJhzMpMYFgMZXRjQ5EC+zf9jHdlJgDYnrsduhzdIntKCBEKBSxZE5weJ/7a+1d4A94F92W9DDxdKWA9s1d5Do8dbXMKSaT4MwQJV5ZhYVa2YTJ1KOY+QSYAi6oD/3979xLcxJW3DfxpSZbkW0u+AcaWDQ4mYBsIMR5icUmqhnrLM/VVfRXeqvF8NZuZhZnVQBb2MtSE2cUswuwSpYos0SbL2ItkZhKwPG9IZipB8hBewGBxN9hS+y6p1d8CpFjW1S213bKeXxULt7pbB1rywzl9+n8Ote/HW81vrfs9DIIBR5uPYpd9Vx4tJaJCY8BS0VsML+Jvk3/DciT1ijOrRZcFhO5YoYRehutSeBE/PruBaPTn+7UWuQLtM2/mvVh6LsvMRYwhPKy9iV92vI322vZ1v4fZaMaJ1hNcao5IhxiwVNSWI8v42+TfsBhezLpvdFlA6LYVSvhluIbkEH54+kNCIQmzbEX7i24Yo+t73nStxbIg7tb8gIghfY96xbSAF9vu4b8P/F80VDSs+z0qzZV4Z9c7EC1iPk0lIo0wYKloheQQ/j75d8ytzGXdN7r0KlwjL8NVVmTcePpjQiEJU9SMPTPdMEUtebVrpvwR/OJ/oAjpHxGaN89CdgTx//b3q3qUpra8Fm/vepur4RDpGAOWilJYDuPvk39HYDmQdd/oogGhO5Z4uCpQ4Jv2YW5VIQlj1IT2mW6Y5fI8WqXgoXgL0xXJi7GvNlv+GHWvWXFqz3/DKKz/Hm+T2ASnw6lqIhQRbRx+Q6noyFEZ39z/BjNLM1n3jS68Clf558dW7szewcziz8caFSP2zHbDEqlU3yZDGJP2HzFvns2437R4D2907VNddL+9rh3djd18DIeoCDBgqahElSiuTl3Fs4X0VZDi+86/Ctfoz2H0ZP4xHgQfxH82KAa0zR5GebhadZuWyiRM1vyAkGEl7T6KIGO2/gH+642TqopHAMDhxsPYV79PbTOJaIMxYKloKIoCj9+Dx3PpH3mJkSUDwpOJ4SqtBHFr5lb8ZwECdgUOojKkvlZvLvdbZUMYEUcApw/+H5Sb1j8EbTQY0dvcC4fNobqdRLTxGLBUNP7n4f/AH8x8fxMA5KAR4XvmhHBdkVfgXbOua2ugC+JKvcrWKHhgu4nn5Q8z7rViWkT9fjOce34FA9ZfatFisuBk60nUV6htJxFtFgYsFYXvHn2HydnJrPvJASPC9y1YXYY4qsi48fQGQqsqPDmkfbAvb1fVllyebwWAFes8Orod2Ltt/c+3AkC1pRpvt76Naov64Wsi2jwMWNK9H578gP9dU8YwFTlgRPieBWsX0Ln5/GbC0nM75ttQt9isqi3z5lncs9/I+HwrAMi2BZx86w3UVtSqep/6inqcbD0Jiym/R4aIaPMwYEnXfM98mJieyLpfunD1B6fwbGE6/nPtUiN2zLepasvTqkk8qbqTrprwS4ICS7OM428eW/fi6DEOmwO9zb0wGvKvgUxEm4cBS7p168Ut/Pj0x6z7vbznmhyus8szuBu4G/+5KlSDlmDHutshG0K4Z/dizpz5sSDFGIWj044Dr6mf6buvfh8ONx5WfTwR6QcDlnTp7uxdfP/o+6z7pQvX5cgyJqYn4tutkUq0zR4CsL7nR3MdEjZagYM9r2FH/fpLHgKAIAjobuxGe526+7VEpD8MWNKdqeAUvn34bdb94uG6Zl31qBKF75kX4Vc1hk1RM9pmD6+7eH9OQ8IAqmos6H7rdVRY1ZUtNBlMcDqcaBKbVB1PRPrEgCVdeTT3COP+cShru6RrpAtXALg1cyteBtGgGPDa7Bswy7mHn2wI457tBuYsmYeEBQHY2VKHrsNtMKisrGQ1WfH2rrdRW65uMhQR6RcDlnTj+eJzXJu6hmiq1FxFlgxpw/Xx/CM8mXsS/7k10IXycO6rzSyYA5i0/5h1SNhsMuP1zmY0t6l/PtVuteNk60lUmtWXaCQi/WLAki7Mrczhm/vfQI7KGfd7WaHJmjJclyJLuD1zO/5z09xe2FZyXyf1aeU9PKm+nXVIuLaqBl3du1BZW5bzuddqFpvR6+hlwX6iLYzfbtp0y5Fl/OPeP7ASSV/LF8gcrgoU/Of5BORXC6fXLzajYaElp/eXDWHct3khWV5k3M9oMKBt+y60vdEAgzVbDKfX0dCBQzsOqT6eiIoDA5Y2VSQawdf3vk4oBJHKz7WFU7/uD05BWn65LmxFWESz9HpO7y9ZpjFlm0DEEM64X7WlGh2te1DzehkEk7pwNRqMONp0FK32VlXHE1FxYcDSplEUBdemrmVddi5V4f7V5kNzuBe8BwAwKEbsCh5AtsdxokIED8SfMFOeeeEAAUCLvQVtbc2wtESgdpW48rJynGw9yclMRCWEAUub5tuH32ZdGUeeyxyuUSWKm89vxov4O6R9MEcyr1gzZ57BlN2HcIbl5QDAarJgX8N+NLRVwLQtknHfTGrLa3Gy9STKy/JZzJ2Iig0DljbFjac3cHf2bsZ95DkDwnfThysATAYmMR9aAADULO9AzVJj2n2jgoxH4q2sK+AAwLbKbdi7rR0VbTKMovpwbbW34mjTUZY9JCpBDFjacLdnbsP7zJtxn1zCNbgSwINXy9dZ5HI4gulLFC6YA7hv8yJkXM74viajEe21e7GjdhvMbct5TWY6uP0gOrd1qj6eiIobA5Y21EPpIb579F3GfXIJV1mRcXP6ZvyRmtZAV8pKTVFBxpOqO3hWOZW1bTariP31Haiwl8G8ewmCym+HyWBCr6MXzaK6FXuIaGtgwNKGebH4AmP+sYxVmqLz2cMVAO7M3MZS5GVvtCpkR0XYtmYPBS8qHuJx1d2sRSMEQcBu+y44bC0oq5dhal5RPZmp0lyJk60nYbfa1Z2AiLYMBixtiLmVOXx9/+uMhSSi8waE7mQP1/nQfMLkqNqlnQmvS5ZpPBT/FyvGxaztEq3V2Fu7F1WWapTtDOU1mamhsgEnWk5wDVciAlAEAetyuQAAoviy3F1/f3/WY9xuN2y2lz2aYDCY0zGknVwKSeQargBwP3AvqdpS2LCCoPUZXpQ/wlLZXNZzmIwmvFbThsaqRghGoGzXMoxi5hKNmbTVtKGnqQcGwaD6HES0tej6t8H58+fhcDgwMDCA/v5+SJKE0dHRjMfEwrWvrw99fX2w2Wxwu90b1GJaK5dCEusJ14XwAp4vPk/YNmWbgG/bVTwQf8oargKAHdU7cLTpF2is2gmDRYF5r/pwFQQBbza+iaPNRxmuRJRAt78RJEnCyMgI+vr64tv6+/tx8eLFjMe53e6EY/r6+hiwm0RRFIxNjWUsJBFdMCCUwz3XmPuB+1lrBadTWVaJNxrfwL66fSgzmGGokmHeq36mcJmxDG+3vo3X63OrGkVEpUW3Q8Rerzc+zBsjiiL8fj8kSYoPGa8mSRL8fn/Sdp/Pl/YY0s71R9fxaO5R2tejC696rnJu4boYXsD0wrN1t8NoMGCXfReaRQeEVxWeTPURmJpDqiczVVuqcbL1JEQLP1NElJpuA3ZtuMaIoohgMJgyLFOFcuwYr9cLp9OZsD0QCCAQCCRse/LkCSh/3mde3Jm5k/b16KKwrnAFgKcLT9fde62vqEN7bTssppfrwQoCYGoKwdSgfjLT9qrtON5yHGajWfU5iGjr023AdnZ2IhgMJmyTJCneS3U4HEnHpOulpgvrjz76CB988EHCNpPJhLa2tjxaTndn7+LG0xtpX48uCwjdsa4rXAFgdnk2532rzFXYXbMbdeV18W1CmQLzrhUYqtRPZmqva0d3YzcEtV1fIioZur0HCwCDg4PxWcTAyx6qKIppAxN4GbK5eu+99zA5OZnw5+rVq3m1udQ9nnuMbx9+m/b16IqA0G0rlMh6A0rJuuIOANSU1+DQjoM4svNIQrgaKqKw7F1WHa4GwYCeph4c2XmE4UpEOdFtDxZ4OanJ5/PB4/EAAJxOJyRJStl7BZD2HuvannCM3W6H3Z5YEMBk0vU/ia7NLM3g2tS1tIUklNCrcA2rCSgBZoMZy9HkR30EAWio2IYWmwNV5uqk1421EZQ5QlA7yddisuB4y3Fsq8x98XYiIt2nSWfnz7VcY0PA6YK0q6srZZhKkoSuri7N2kgviz98fe9rRKKp720qYSB02wIlpL7391rta/jpxU+IyDKMBgPsVjtqy2tRX14fv8e6WiHut9qsNpxsPYkqc5XqcxBRadJ1wI6OjiY8cjMyMoIzZ86k3T82fLz2XqzD4eAMYg2tRFbwj3v/wHIkdSF9JQKE7lgRXcnvjkRDxTbUVzQgqsgwCkZkWvNVMCko27UCY7X6+61NYhN6m3tRZixTfQ4iKl26vgf7ySefxO+pxopMrK7K5PP5cPbs2YRj+vv7MTIyEv/Z7XZjYGBgYxpcguSojG/uf4O5ldQFHhT5VbguFeajJkCAUTAhU7gaKqKwvL6cV7jub9iPEy0nGK5EpJque7CDg4PxIhGBQACXLl1K6IkGg0GMj48n9FgHBgbgdrvh8Xji4cxSidpQFAVj/rGkykrx12UgdNeC6OLG/T/OWPPqfqvK5VcNggG/aPoFdtfsLmzDiKjk6DpgnU5n0rOra1+/fv160nYG6sb495N/46GUevFyJQqE71kQnd+ghcYFoKwxBNN29fdby8vKcbzlOOor6gvYMCIqVboOWNKvu7N38dPzn1K+pihA+L4FsrQx4SqYFJS1ruRVrH9b5TYcazkGa4rJUkREajBgad2mF6Zx/WHyyEFM2G+GHNiYcDVURGHevQLBrLZCMbCvfh/e2PEGn28looJiwNK6LIYXcW3qGqJK6t5i+GEZ5Bcb87HK9/lWk8GEo81H0WJrKWzDiIjAgKV1iM0YTvc4TuSJCZFn2s+6LcTzrdWWapxoOQGbNX1VMCKifDBgKWfjD8Yxu5S6HnBk2oTwY+2L3xeinnCz2Iy3mt/iIzhEpCkGLOXE+8wLfzB5KUAAkGeMCD/QPlwNlVGYd6m/3yoIAg5uP4iOho4Ct4yIKBkDlrLyB/1pV8eRg0aEpyyat8FUH4GpKb96wk6HEzuqdhS2YUREaTBgKaPAcgD/fPDPlK9F5w0I3zMjTW3/ghAMgKk5BFOd+vutteW1ON5yHJXmygK2jIgoMwYspbUSWcE3979JWcA/uiggdNcCJardoy2C+dX91kr191vbatpwZOcRGA0bVPCCiOgVBiylFFWiuDZ1DQuhhaTXlJCA0N31L5i+HkabjLKWFQgqP6EGwYDund3YU7unsA0jIsoRA5ZS+vfjf+PZwrOk7bH6wurWdM1OEABTniUPq8xVON5yHDXlNQVsGRHR+jBgKcnd2bu49eJW0nZFeVVfuEAr46wlmBWYW/N7BKdJbMJbzW/BbNR+VjMRUSYMWEowszSD7x59l/K1yAOzZvWFjaKMslb1Q8KCIODQ9kPY37C/sA0jIlKJAUtxy5FlXL1/FXJUTnot8tSEyPPCf1wKMSRcXlYOp8OJbZXbCtgyIqL8MGAJwMu1XT1+DxbDi0mvybNGhB8Vfsi1EFWZuAoOEekVA5YAvKzU9HT+adL26IIB4anCh2u+Q8IA0NHQgYPbD3IVHCLSJQYs4cn8E/imfUnboysaPOsqAGU7QjDtUD8kbDaa0evoxc7qnYVrFxFRgTFgS9xSeAnj/nEoa8oxKREgfMcCJVK4cC3EkDCrMhFRsWDAlrDYfde1y88pUSA0aUF0pXCP4xRiSLi9rh1vNr4Jg9qCxEREG4gBW8JuPLuRsphEeMqM6HyBHscp0JBwT1MPF0YnoqLCgC1Rj+cew/cs+b5r+FEZ5NnCfCwKNUu419GLirKKgrSJiGijMGBL0GJ4EeMPxpO2R16YEHlamEXIjdWvhoRVni62duv++v2cJUxERYkBW2Ji911XIisJ2+U5AyL+/B/HEQTAtCMM4/Yw1OZitaUaTocTteW1ebeHiGizMGBLjPeZF9ML0wnblDAQvm/Je11XgzWKstYQDBX5LS/XvbMbJgM/mkRU3PhbrITMLM2kfN41PJXn6jgCYNoWhmlHGGon+HIiExFtNQzYEiFH5ZTPu0aemvIq4G+wvOq15rEoOicyEdFWxIAtET88/QHSipSwLbpgQOSx+vuupoYwTI1hCCrzWRAEHNh2AB0NHZzIRERbDgO2BDxbeIafnv+UsE2JAKF76u67CmYFZS0rMFar77VWmavgdDhRV1Gn+hxERHrGgN3iwnIY/3zwz+TtUxYoofX3Gk11EZiaQqp7rQCwu2Y3uhu7UWYszCNBRER6xIDd4v71+F9YCC0kbItMmyAH15eQQtmrXquovtfKiUxEVEoYsFvYQ+kh7s7eTdimhAREHq2v52isiaCsOZRXHeGGygY4HU5OZCKiksGA3aLkqIzvH3+ftD380Jzz8nOCSUGZIwSjXVbdDqPBiEPbD2Fv3V5OZCKiksKA3aImpieShoblOQPkQPahYUEAjHURmBrz77UebTqKaku1+pMQERUpBuwWNB+ax8T0RNL2yMPsj+QY7TJMjSEYrOrLOrHXSkTEgN2Svn/0PaJK4mSk6KIB0aX0ZZYMlVGU7QzltfINwF4rEVEMA3aLeSg9xKO5R0nb0w0NGyxRmHaG87rPCrDXSkS0FgN2C5GjMv71+F+pX1w14isIgEGUYayNwGCTVa96E8NeKxFRMgbsFnJ75jbmQ/MpXzPtCEOwKhAMCgxVsup1WlcrM5bh0PZDaK9rz/9kRERbDAN2i4gqUdx8fjPt64LxZRWmQtll34XDjYdhNVkLdk4ioq2EAbtF3Avcw2J4UfP3ES0ijuw8gu1V2zV/LyKiYsaA3SL+M/0fTc9vNBjR2dCJ/Q37YVC76CsRUQlhwG4BD6QHSUvRFdLO6p3o3tmNKnOVZu9BRLTVMGC3gFsvbmly3ipzFQ43Hkaz2KzJ+YmItjIGbJGLRCOYXpgu6DktJgu6tnVhT+0eDgcTEanEgC1y0wvTSVWb1DIZTNhXvw/76vdxrVYiojwxYIvc04WneZ9DEAS8VvMaDmw/wMduiIgKhAFb5GaWZlQfKwgCWmwt6NrWBdEiFrBVRETEgC1yZmP2FXJSHbOndg/a69q5ADoRkUYYsEVuPT3Paks1Xq97HbtrdsNk4KUnItISf8sWub11ezE5O5m2ipNoEdFY3Yim6iZWXyIi2kAM2CJnNVnRt6cPk4FJvFh8AQUKzEYzaqw1aKxuZHEIIqJNwoDdAiwmC/bV79vsZhAR0SqsIkBERKQBBiwREZEGGLBEREQaYMASERFpgAFLRESkAQYsERGRBhiwREREGmDAEhERaYABS0REpAEGLBERkQa2XMD6/X54PJ7NbgYREZU43dcidrvdAABJkgAAAwMDGff3eDy4ePEiJEmCKIro7e3F0NAQHA6H5m0lIiKK0XXAulyuhECVJAnDw8MYGhpKe4zNZsP169fjAUtERLQZdD1E7Pf7E34WRRFzc3M5HctwJSKizaTrgPV6vUkhS0REVAx0PUQ8ODiI06dPY3BwEP39/RgdHcXg4GDW41ZPcvL5fGnv2wYCAQQCgYRtT548ya/RRERE0HnAOp1OfPbZZzh9+jRcLhcuX76cdei3s7MTNpstvl+m+7YfffQRPvjgg4RtJpMJbW1thftLEBFRSdL1EHHskZvr16/D6XTi1KlTWR/BcTgcCSHc19eHTz/9NOW+7733HiYnJxP+XL16taB/ByIiKk267sG6XC5cuHABAHDhwgX09fXh3Llz+Oqrr9Y1iUkURfh8PnR2diZst9vtsNvtCdtMJl3/kxARUZHQbQ/W4/HA6XQmbHM6nThz5kzaXqzP50NPT0/K1zirmIiINpJuu2uxwhJrre2FriaKYtIkqNh5WGiCiIg2km57sH19ffjiiy+Sto+NjcV7tj6fD2fPno2/lipEP/7445xmHhMRERWSbnuwADA0NITh4WHY7XaIoghJkvDHP/4x/nowGMT4+HhC1aZf/epXcLlcEEURU1NTaGlpQX9//2b9FYiIqETpOmAdDkfGsohOpxPXr19P2CaKYtZ6xURERFrT7RAxERFRMWPAEhERaYABS0REpAEGLBERkQYYsERERBpgwBIREWmAAUtERKQBBiwREZEGGLBEREQaYMASERFpgAFLRESkAQYsERGRBhiwREREGmDAEhERaUDXy9VtBlmWAQBPnjzZ5JYQEZEexPIglg+5YsCuMT09DQD43e9+t8ktISIiPZmenkZra2vO+wuKoigatqfoLC8vw+v1oqGhAUajcdPa8eDBA5w4cQJXr15Fc3PzprWDUuP10T9eI/0rlmskyzKmp6fR1dUFq9Wa83Hswa5htVpx5MiRzW4GIpEIIpEIduzYoesPXqni9dE/XiP9K6ZrtJ6eawwnOREREWmAAUtERKQBBiwREZEGjH/+85//vNmNoNSsViveeeeddd1Up43D66N/vEb6t5WvEWcRExERaYBDxERERBpgwBJpxO/3w+PxbHYziIpWsX+H+BzsBnG73bDZbACAYDCIqcb2ewAACCZJREFU/v7+jPufPn0aPp8vabvT6cTly5dVnZMyK/Q18ng8uHjxIiRJgiiK6O3txdDQEBwOhybtLwVqPvNutxsAIEkSAGBgYCDvc1Jqhb4+Rf8dUkhzV65cUUZGRuI/j4yMKFeuXMl4zIcffqh4vV5lamoq/mdkZETxer2qz0npaXGNYucLBoPaNbyEqLlGn3zyScLPwWBQ+fDDD/M6J6WmxfUp9u8QA3YDvPvuuzltiwkGg/Ff0qut/rCu95yUmRbXaPUvG8qfms/8+++/n3Ebv0eFo8X1KfbvEO/BakySJPj9/qTtPp8vPiSyliiK6OzsTNjmdrvjwy1qzknpaXGNqLDUfua9Xm/K4/I5JyXT4vpsBQxYjXm93vg9idVEUYTX683pHJIkJdxzKMQ56WdaXKMYj8cT/+NyufJua6lSe40GBwdx+vTp+H2+0dFRDA4O5nVOSqbF9Ykp5u8QA1ZjsZvza6X6MKbjdrvhdDoLek76mRbXCAA6OzvR1dUFp9MJp9MJh8OB4eHhvNtbitReI6fTic8++wznz5/HqVOn0NnZGT8Pv0eFo8X1AYr/O8SA3QD5DDdJkoRAIFDQc1IyLa6Rw+FI+GXR19eHTz/9VPX7lDo11yj2mMf169fhdDpx6tSphMc++D0qHC2uT7F/hxiwGkv1vzrg5RT2XLjdbrS0tBT0nJRIi2uU6b1SPdpDmam9Ri6XCwMDAxBFERcuXMDly5dx7ty5tD2uXM5JybS4Ppneq1i+QwxYjXV1daX8kEmShK6urqzHezyepGGWfM9JibS4Rj6fDz09PSn3T/fLiNJTc408Hk/SsL3T6cSZM2fg8Xj4PSogLa7PVvgOMWA1JooibDZb0v/I1g59pOPxeJImz+R7Tkqk1TVaO1kjdv6ieUheR9Rco3S9oNh9Pn6PCker61Ps3yEG7Abo7+/HyMhI/Ge3251QrcTn8+Hs2bNpj0/1Ac12TlqfQl+jVL8APv7446RfGJS79V6jvr4+fPHFF0nnGRsbi/ec+D0qnEJfn63wHeJqOhvE7XbD4XBAkqSkEmIejwfnzp3DV199lfSL+tSpU/j8889Thmymc9L6FfoaSZIEt9sNURQxNTWFlpYWXqM8rfca+f1+XLlyBXa7HaIoQpIk9Pf3J1wrfo8Kp9DXp9i/QwxYIiIiDXCImIiISAMMWCIiIg0wYImIiDTAgCUiItIAA5aoSG3lVUhSKbW/LxU/BixREfJ4PCVXR9fn8zFkqagwYIl0yOVyweVy4fz580mh4vf74fP5ktaj1QuPxwO32x1fXiyXUIw97xjb1+fz4fz58wn/iejr68OVK1c0azdRwW3ueu9EtNbY2Jjy7rvvKsFgUNm7d68yNjaW8Pr777+/SS3LbmpqSvnTn/6UsO33v/991uO8Xq+yd+/e+J9f/vKXitfrTbnfyMhIwdpLpCXTZgc8ESW6cuUKent7IYoivvzyy4SScXofInW5XPjtb3+bsK2jowOjo6Po6+vLeOyXX34ZXzQhXf3azs5OXLx4Meu5iPSAQ8REOjM+Po4DBw4ASK5p7PF4dB0uIyMjSaunHDhwIGGNz0xiRd4z6ejoKJrlyqi0MWCJdESSJEiSlLSMV8zqQvV64/f7U66zKopizgGbi2PHjqUsEk+kNwxYIp0YHh7GuXPnAAAXL17E+fPnk/aZmJjIeo6enp74sbHJQxsh3axmm82W09C2z+fD6OgoRkdH4XK50vZSu7q6sv47EOkB78ES6cTQ0BBcLhcA4MKFC0mvp+odrubxeDA3N4e//OUv8Pv9GB4exvj4OD777LOM7/uHP/wh5zZ2dHRgaGgo5WupFtzOVezvtXr4+/Tp07h06VLKtXb1fi+aCGDAEunKjRs30NHRkfI1v9+P5ubmtMd2dXUlDB+fP38ely5dynpP8/Lly+oaW0AOhyMpSHt7e+FyuVL+Z4OoGHCImEhHJiYm4hOc1goGg7Db7WmPXR2kw8PD6OvrS7lotVZiM4DXUtuzbWlpKei9W6KNxoAl0hG/3582FG02GwKBQNZzDA8P49e//vWGT4aKtXvtvdhsQ9sA4kPjq9lstryGnYk2G4eIiXQiNqknXYUmURQxNzeX8RzDw8M4duzYuqo8FeoebOwRm2AwmBCowWAQvb29ac/p9/tx8eJF9Pf3Jx23kT1wokJjwBLphNfrzRgo2WbjxsJ1dc/V5XJhYGAg4/sW8h5sb28vfD5fwt/D5/Ph2LFjaY9xOBy4cOFCUi93bGwsZTDn0iMm0gMOERPpxNTUVNoJTkDmIgyjo6OYmJjA2NgYXC4XRkdHcfbsWa2amtbQ0FDSY0Ferxf9/f3xn/1+P06dOpV0f3X1fx58Ph8mJiZS9paDwWBSMQsiPWIPlkgnJiYmslZpam5uTnmf1u/34/Lly5AkCefOnYPX68WZM2ey9l4LzeFwYGBgAC6XCw6HAzdu3MClS5eS9lt7b7W/vx9utxuSJCEQCGBubg6ff/55yvfQezUrohhBURRlsxtBREBPTw8+//zzjMPEHo8HPp9vw4NTT86ePYu//vWvm90MoqzYgyXaRD6fLz7sa7PZsk7qcTqdJb9kW6ZHlYj0hD1Yok3U09ODM2fOQBRF2Gy2nIY+PR4PbDabbteD1ZLL5drw53uJ1OIkJ6JN9Jvf/AaiKGJqairn+4pOp7MkCzDEZg8zXKlYsAdLVIQkSSq5yT65rClLpCcMWCIiIg1wiJiIiEgDDFgiIiINMGCJiIg0wIAlIiLSAAOWiIhIAwxYIiIiDfx/FBaiMD+peP8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set_style(\"ticks\")\n",
    "\n",
    "# Choose what to plot\n",
    "zi = 0; par0 = 'f'; par1 = 'sv';\n",
    "\n",
    "# Get the indices in the Covariance Matrix of the two parameters we want to plot\n",
    "p0i = varyzdep.index(par0); index0 = zdepind[:,zi] # change this to shapeind if using shape!\n",
    "p1i = varyzdep.index(par1); index1 = zdepind[:,zi] # change this to shapeind if using shape!\n",
    "\n",
    "# Get the fiducial values for the ellipse centre\n",
    "c0 = fids[par0][zi] # remove the [zi] if par0 is a shape parameter!\n",
    "c1 = fids[par1][zi] # remove the [zi] if par1 is a shape parameter!\n",
    "\n",
    "# Marginalise out everything else and make a new 2x2 Covariance Matrix\n",
    "marg_all = True\n",
    "if marg_all:\n",
    "    covmat_2D_comb = np.zeros((2,2))\n",
    "    covmat_2D_comb[0,0] = CovMat_comb[index0[p0i],index0[p0i]]\n",
    "    covmat_2D_comb[1,1] = CovMat_comb[index1[p1i],index1[p1i]]\n",
    "    covmat_2D_comb[0,1] = CovMat_comb[index0[p0i],index1[p1i]]\n",
    "    covmat_2D_comb[1,0] = covmat_2D_comb[0,1]\n",
    "    \n",
    "else:\n",
    "    fishmat_2D_comb = np.zeros((2,2))\n",
    "    fishmat_2D_comb[0,0] = Fishermat_comb[index0[p0i],index0[p0i]]\n",
    "    fishmat_2D_comb[1,1] = Fishermat_comb[index1[p1i],index1[p1i]]\n",
    "    fishmat_2D_comb[0,1] = Fishermat_comb[index0[p0i],index1[p1i]]\n",
    "    fishmat_2D_comb[1,0] = fishmat_2D_comb[0,1]\n",
    "    covmat_2D_comb = np.linalg.inv(fishmat_2D_comb)\n",
    "    del fishmat_2D_comb\n",
    "\n",
    "def get_ellipses(centre,CM_2D,sty={'w':0,'s':'-','f':True,'a':0.7,'c':'green'}):\n",
    "    # We will need coefficients for the ellipse depending on the number of sigma level we want\n",
    "    CL_alpha = {\n",
    "        1: 1.52, \n",
    "        2: 2.48,\n",
    "        3: 3.44\n",
    "    }\n",
    "\n",
    "    # Get the errors\n",
    "    FM_2D = np.linalg.inv(CM_2D)\n",
    "    em0 = np.sqrt(CM_2D[0,0])\n",
    "    eu0 = np.sqrt(1./FM_2D[0,0])\n",
    "    em1 = np.sqrt(CM_2D[1,1])\n",
    "    eu1 = np.sqrt(1./FM_2D[1,1])\n",
    "\n",
    "    #print par0, eu0, em0\n",
    "    #print par1, eu1, em1\n",
    "\n",
    "    # Now calculate the ellipse parameters\n",
    "    a = np.sqrt(0.5*(CM_2D[0,0]+CM_2D[1,1]) + np.sqrt(0.25*(CM_2D[0,0]-CM_2D[1,1])**2 + CM_2D[0,1]**2))\n",
    "    b = np.sqrt(0.5*(CM_2D[0,0]+CM_2D[1,1]) - np.sqrt(0.25*(CM_2D[0,0]-CM_2D[1,1])**2 + CM_2D[0,1]**2))\n",
    "    theta = 180./np.pi * 0.5*np.arctan(2.*CM_2D[0,1] / (CM_2D[0,0]-CM_2D[1,1])) # should be in degrees\n",
    "\n",
    "    if CM_2D[0,0]<CM_2D[1,1]:\n",
    "        width=2.*b\n",
    "        height=2.*a\n",
    "    else:\n",
    "        width=2.*a\n",
    "        height=2.*b\n",
    "\n",
    "    # Create all the ellipses   \n",
    "    ellipse1 = Ellipse([centre[0],centre[1]], width*CL_alpha[1], height*CL_alpha[1], theta, \n",
    "                       linewidth=sty['w'], ls = sty['s'], fill=sty['f'], zorder=2, alpha=sty['a'], color=sty['c'])\n",
    "    leg1 = ellipse1.set(clip_box=ax.bbox, alpha=sty['a'], facecolor=sty['c'])\n",
    "    ellipse2 = Ellipse([centre[0],centre[1]], width*CL_alpha[2], height*CL_alpha[2], theta, \n",
    "                       linewidth=sty['w'], ls = sty['s'], fill=sty['f'], zorder=2, alpha=sty['a']/2, color=sty['c'])\n",
    "    leg2 = ellipse2.set(clip_box=ax.bbox, alpha=sty['a']/2, facecolor=sty['c'])\n",
    "\n",
    "    return ellipse1, ellipse2, CL_alpha[1]*em0, CL_alpha[1]*em1, CL_alpha[1]*eu0, CL_alpha[1]*eu1, leg1, leg2\n",
    "\n",
    "# Make a test plot\n",
    "ell1, ell2, errm0, errm1, erru0, erru1, leg1, leg2 = get_ellipses([c0,c1],covmat_2D) #without priors\n",
    "ell1_p, ell2_p, errm0_p, errm1_p, erru0_p, erru1_p, leg1_p, leg2_p = get_ellipses([c0,c1],covmat_2D_comb,\n",
    "                                                    {'w':0,'s':'-','f':True,'a':0.9,'c':'tab:purple'}) #with priors\n",
    "ells = [ell1,ell2,ell1_p,ell2_p]\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "\n",
    "#ax.axhline(c1,ls='--',c='k',lw=1); ax.axvline(c0,ls='--',c='k',lw=1) # Centre lines\n",
    "\n",
    "ax.add_artist(ell1)\n",
    "ax.add_artist(ell2)\n",
    "ax.add_artist(ell1_p)\n",
    "ax.add_artist(ell2_p)\n",
    "ax.set_xlim(c0-errm0*2.1, c0+errm0*2.1)\n",
    "ax.set_ylim(c1-errm1*2.1, c1+errm1*2.3)\n",
    "ax.tick_params(color='black')\n",
    "ax.set_facecolor('none')\n",
    "#ax.axis('off')\n",
    "\n",
    "#plt.xlabel(par0+'(z='+str(z_vals[zi])+')')\n",
    "#plt.ylabel(par1+'(z='+str(z_vals[zi])+')')\n",
    "plt.xlabel('$f(z=0.5)$')\n",
    "plt.ylabel('$\\sigma_v(z=0.5)$')\n",
    "\n",
    "ax.legend([ell1,ell1_p], [r'no prior',r'with $10\\%$ prior on $\\{\\sigma_v,b_1,N\\}$'],loc='upper left',frameon=False)\n",
    "\n",
    "#plt.savefig('contour.pdf', format=\"pdf\", bbox_inches = 'tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "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
}