notebooks/datastruct.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import yaml\n",
    "from classy import Class\n",
    "from iminuit import minimize as minuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'size': 14}\n",
    "axislabelfontsize='large'\n",
    "matplotlib.rc('font', **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pybird.correlator import Correlator\n",
    "from pybird.io_pb import ReadWrite\n",
    "from pybird.likelihood import Likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data format for EFT likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best way to understand what is the format of the data read by the EFT likelihood is to construct one. \n",
    "\n",
    "***Let us construct a fake data from an EFT prediction.***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = os.path.join(os.path.expanduser('~'), 'Documents', 'github', 'pybird', 'data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting the EFT prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Redshift?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = 0.57"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we create a CLASS engine `M`, that will provide linear cosmology stuff to the PyBird correlator engine. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = Class()\n",
    "cosmo = {'omega_b': 0.02235, 'omega_cdm': 0.120, 'h': 0.675, 'ln10^{10}A_s': 3.044, 'n_s': 0.965}\n",
    "M.set(cosmo)\n",
    "M.set({'output': 'mPk', 'P_k_max_h/Mpc': 10, 'z_max_pk': z})\n",
    "M.compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d\n",
    "kk = np.logspace(-5, 0, 200)\n",
    "pk_lin = np.array([M.pk_lin(k*M.h(), z)*M.h()**3 for k in kk]) \n",
    "ipk_lin = interp1d(kk, pk_lin, kind='cubic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = M.scale_independent_growth_factor_f(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need some distortion (AP) parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_ap, D_ap = M.Hubble(z) / M.Hubble(0.), M.angular_distance(z) * M.Hubble(0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Omega_m_fid: 0.310\n"
     ]
    }
   ],
   "source": [
    "M_fid = Class()\n",
    "cosmo_ap = {'omega_b': 0.02235, 'omega_cdm': 0.120, 'h': 0.6776, 'ln10^{10}A_s': 3.044, 'n_s': 0.965}\n",
    "M_fid.set(cosmo_ap)\n",
    "M_fid.set({'output': 'mPk', 'P_k_max_h/Mpc': 10, 'z_max_pk': z})\n",
    "M_fid.compute()\n",
    "\n",
    "Omega_m_fid = M_fid.Omega_m()\n",
    "H_fid, D_fid = M_fid.Hubble(z) / M_fid.Hubble(0.), M_fid.angular_distance(z) * M_fid.Hubble(0.)\n",
    "\n",
    "print ('Omega_m_fid: %.3f' % Omega_m_fid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create a PyBird Correlator engine `N`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "dk = 0.01\n",
    "kd = np.arange(0.005, 0.3, dk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = Correlator()\n",
    "N.set({'output': 'bPk', 'multipole': 3, 'z': z, 'km': 0.7, 'kr': 0.35, 'nd': 4e-4, 'kmax': 0.3, 'xdata': kd, \n",
    "       'eft_basis': 'westcoast', 'with_stoch': True, \n",
    "       'with_ap': True, 'H_fid': H_fid, 'D_fid': D_fid,\n",
    "       'with_binning': True, 'binsize': dk,\n",
    "       'with_bias': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us compute an EFT prediction on some fiducial EFT parameters. \n",
    "\n",
    "We load the fiducial EFT parameters from a previous fit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'b1': 1.9535, 'c2': 0.5863, 'c4': 0.0, 'b3': -0.3948, 'cct': 0.1839, 'cr1': -0.8414, 'cr2': -0.8084, 'ce0': 1.5045, 'ce1': 0.0, 'ce2': -1.6803, 'b2': 0.4146, 'b4': 0.4146}\n"
     ]
    }
   ],
   "source": [
    "outdir = os.path.join(datadir, 'eftboss', 'out')\n",
    "with open(os.path.join(outdir, 'fit_boss_onesky_pk_wc_cmass_ngc_l0.dat')) as f: data_file = f.read()\n",
    "eft_params_str = data_file.split(', \\n')[1].replace(\"# \", \"\")\n",
    "eft_truth = {key: float(value) for key, value in (pair.split(': ') for pair in eft_params_str.split(', '))}\n",
    "print (eft_truth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We give the linear power spectrum together with a bunch of cosmological parameters and EFT parameters to the PyBird Correlator engine: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.compute({'kk': kk, 'pk_lin': pk_lin, 'f': f1, 'H': H_ap, 'DA': D_ap, 'bias': eft_truth})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the multipoles of the galaxy power spectrum: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "bpk = N.get()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we need a covariance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import legendre\n",
    "\n",
    "def get_cov(kk, ipklin, b1, f1, Vs=3.e9, nbar=3.e-4): # BOSS is about Vs ~ 3e9\n",
    "    dk = np.concatenate((kk[1:]-kk[:-1], [kk[-1]-kk[-2]])) # this is true for k >> kf\n",
    "    Nmode = 4 * np.pi * kk**2 * dk * (Vs / (2*np.pi)**3)\n",
    "    mu_arr = np.linspace(0., 1., 200)\n",
    "    k_mesh, mu_mesh = np.meshgrid(kk, mu_arr, indexing='ij')\n",
    "    legendre_mesh = np.array([legendre(2*l)(mu_mesh) for l in range(3)])\n",
    "    legendre_ell_mesh = np.array([(2*(2*l)+1)*legendre(2*l)(mu_mesh) for l in range(3)])\n",
    "    pkmu_mesh = (b1 + f1 * mu_mesh**2)**2 * ipklin(k_mesh)\n",
    "    integrand_mu_mesh = np.einsum('k,km,lkm,pkm->lpkm', 1./Nmode, (pkmu_mesh + 1/nbar)**2, \n",
    "                                  legendre_ell_mesh, legendre_ell_mesh)\n",
    "    cov_diagonal = 2 * np.trapz(integrand_mu_mesh, x=mu_arr, axis=-1)\n",
    "    return np.block([[np.diag(cov_diagonal[i,j]) for i in range(3)] for j in range(3)]) \n",
    "\n",
    "cov = get_cov(kd, ipk_lin, eft_truth['b1'], f1)\n",
    "err = np.sqrt(np.diag(cov)).reshape(3,-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us plot our fake data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdwAAAEWCAYAAADM0CYnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu8VmWd9/HPDzZHkQA5ihzkCdERFYetgUbNWGpMY/XgscjERpHsYNLUy+yRnJk8lGlaE6mIkngcrWmITKQMI9lkew9MksZgCijKBgwRkNOG3/PHWhtv7n2617rXve7T9/16rdfe91rXWuu6WJv929e1roO5OyIiIlJYnYqdARERkWqggCsiIpICBVwREZEUKOCKiIikQAFXREQkBQq4IiIiKVDAFRERSYECroiISAoUcEVERFJQU+wMlJP+/fv7yJEji50NEREpIQ0NDVvcfUBH6RRwIxg5ciT19fXFzoaIiJQQM1uXSzo1KYuIiKRAAVdERCQFCrgiIiIpUMAVERFJgQKuiIhIChRwRUREUqCAKyIikoKSCbhm9nUz+4OZvW1mm83s52Y2NiuNmdn1Zva6me0ysyVmdnxWmr5mNt/MtoXbfDPrk5XmBDN7JrzGBjObZWaWRjlFRKQ6lUzABf4OmA2cBpwBNAG/MrN+GWm+BnwF+CJwCrAJWGxmh2ekeQj4W2Ay8JHw+/nNB82sN7AYaAyv8SXgq8DMQhRK0nHhXXVceFddsbMhItKmkplpyt3PzvxsZhcD24DTgZ+HNdAvAze7+0/CNJcQBN1PAXeZ2XEEQfb97r4sTHMFsNTMxrj7amAq0BO4xN13AavC82aa2W3u7mmUV9rWHDgfvWJizuds372Pt3c30bBuK+NH9C1U1kREYiulGm62wwnytzX8fDQwGHiqOUEYMH9LUCsGmAjsAJZlXOdZYGdWmqXhuc0WAUcCIxMtgcSyffc+Nry1i4Z1WztODDSs28qfN27nta27mHrP8pzPExFJUykH3DuAlUBzO+Hg8GtjVrrGjGODgc2ZtdTw+01ZaVq7RuY9DjKz6WZWb2b1mzdvjlMOiSBO8Fz+8pscCJ/4vqYDLH/5zQLnUkQkupIMuGZ2G/B+4Fx33591OLvJ17L2tdYk3FEaa2M/7n63u9e6e+2AAR0uBiF5ihM8J4w6gk7hE+xS04kJo44oYA5FROIpuYBrZt8DPgmc4e4vZxzaGH7NroUO5N0a6kZgYGaP4/D7AVlpWrsGtKz5SsriBM/xI/py7ODDOapvDx68bILe4YpISSqZTlMAZnYHcBHwd+7+56zDrxAEyzOBP4TpuwOTCHoZQ9D83IvgPW3ze9yJwGEZn+uAb5tZd3ffHe47E3gdWJtwkSSi5uD59u4m7rjo5JyD5xNXfaDAORMRyU/J1HDN7IfApQS1261mNjjcesHBd7G3A9eY2ZRwjO48gk5SD4VpXgSeJOixPMHMJgJ3AQvDHsqEad8B5pnZWDObAlwDqIdyiTi8exeG9umhmqqIVJRSquFeGX79ddb+fwGuD7//DtAD+CHQF/g9cJa7b89IPxX4Pu/2Zl4AfKH5oLtvM7Mzw2vUE/SCvhW4LamCSH6iDAeKK87QIxGRfJRMwHX3Dmd6Cmug1/NuAG4tzV+BT3dwnecBtUGmpBSDm8btikjaSqZJWSQtGrcrIsWggCsFF3Uii0LTuF0RKQYFXCmoUqxNatyuiBRDybzDlcrUWm2y2O9Mx4/oy2MzTmP5y28yYdQRRc+PiFQHBVwpqOba5AEvrdrk+BF9FWhFJFUKuBJZlF7HcSeyEBGpNAq4UnCaBUpERJ2mREREUqGAK1x4V93BZuJclNowHxGRcqCAK5GU4jAfEZFyoIArkVTrpBH/cMdvef+3n9YfGCISmwKuRGoirsZJI1SrF5EkKOBWuajBpBoXe6/WWr2IJEsBt8rFCSbVtl5tNdbqRSR5Godb5eLMBFVKy+ylQVNBikgSFHCrnGaCyo2mghSRfCngCod378Lh3bsooCSoYd1W1YhF5BAKuCIJa1i3lfPvXMYBh+5dOlVN5zIRaZ86TYkkLE5HtDjjfDU2WKS8KOAKj14xseo6QhVS1F7Nccb5xj3nh795SQFapEjUpCySsKi9mlurESd9jpq5RYpPAVekAKL0ap4w6gi6d+nEvqYDOQ/NinpOnKAuIslSwK0wURaHl9IwfkRfHrxsQqRezVHPiTPeWkSSpYArUgLijPONco4m7xApPgVckSoRNahrLLFIshRwRaQFdbISSZ6GBVWYKEvtibRFKySJJE8Bt4Jo3VZJilZIEkmempQriIZ+SFLUyUokeQq4FURDPyRJWiFJJFkKuBVES+1JMalXs0j7FHArjJbak2JQr2aRjpVMpykz+4CZLTCzDWbmZjYt6/i8cH/mtjwrTTcz+4GZbTGzneH1jspKM9zMfh4e32Jm3zezrikUUaRiqVezSMdKJuACvYBVwFXArjbS/AoYkrH9Q9bx24FzgU8Ck4DewEIz6wwQfv0FcHh4/JPAecCtSRZEpNqoV7NIx0qmSdndnwCegKA220ayPe6+sbUDZvYe4J+AS919cbjvYmAd8GFgEXAWcDwwwt1fDdN8DbjHzL7h7m8nVyKR6qFezSIdK5mAm6P3m9km4C3gGeAb7r4pPDYe6AI81ZzY3V81sxeB0wgC7kTgxeZgG1oEdAvP/03hi1BYWrRAiiVOr2Z1tJJqklPANbPuQH93fy1r//Hu/qeC5KylJ4GfAq8AI4FvAU+b2Xh33wMMBvYDW7LOawyPEX5tzDq+JTxvMCKSGnW0kmrT4TtcM/u/wP8CvzCzP5nZ+zIOzy9YzrK4+yPuvsDdn3f3nwOTgTHARzs41QDPvFRbt2j1ZLPpZlZvZvWbN2+OnG8RaV2cjlYN67byw9+8pFnUpCzl0mlqFjDe3U8CLgHuNbNPhcesYDnrgLu/DrwGjA53bQQ6A/2zkg7k3VrtRlrWZPuH52XXfJvvc7e717p77YABA5LIuogQdLTq3qUTnS23jlbNNeJbFq3W1KVSlnJpUu7q7psB3L3ezD4A/NTM3kvbtcWCM7P+wFDgjXBXA7APOBN4KExzFHAcsCxMUwf8PzM7KqN5/ExgT3i+iKRk/Ii+PHjZhJzf4caZulTviKWU5BJwN5nZie7+RwB3f9PMzgR+DJyYVEbMrBfw3vBjJ2C4mY0D/hpu1wM/IQiwI4GbgE3Af4b52mZmc4Fbwo5VbwK3AX8kGE4EQYeqPwH3m9lXgCOAW4A56qEskr4oHa2aa8T7mg5EqhFHeUesAC2FlEvAvRhoytzh7nuBT5rZvyeYl1oO7SX8L+H2Y+BzwAnAZ4A+BEH3N8AF7r4945yrw7w+CvQAfg18xt33h/neb2YfBWYDzxKM930I+OcEyyEiBVDoGrE6cUmhdRhws3smZx17NqmMuPsS2n8nfHYO19gNfDHc2kqzHvjHqPkTkeIrZI1Yq21JoeU8DtfMpsS4/i/dva1Zo0RECiZqjThqgG6mZmjJlbnn1u/JzA5EvLYDo9395ci5KlG1tbVeX19f7GyISIFEDZ5qhhYAM2tw99qO0kWdaWpwxsxOHWVge8epRERKR9TZstQMLVFEWbzgx7S9qEBrHgDU81dEKlbUscSgyTuqWc5NyqImZRFpKUoztJqgK1OhmpRFRCRDlGbozOkrNXlH9YkVcM3sgwSTSgwAXgJWZmwvuHtTO6eLiFSlCaOOoGtNtMk7pt6znL1NB+haoxpxuYtbw51LsFj89whmhzqZYCH3o4C9QPdEciciUkHiTN6xe18wQEQ14vIXN+AOBs7KHvJjZv0Igq+IiLSi0NNZqkZcuuIG3GeAo4FDAq67/5VgOkUREclTGjViUK04LVFmmloE/E+4/Ri4zsz+2LySkIiIJK+QNWLQIg9pilLDXQGMI1jMYFC4b7WZ/Zxg2bsVwB81laNIfg4cOMBrr73Gzp07i52VstWlSxcGDhxI7969i52VVEWtEUO8RR7UbB1PzgHX3a9p/t7MBhME3+btywSdp9zM1rj73ySdUZFqsWXLFsyMMWPG0KlTlLlpBMDd2bVrFxs2bACoyqAbJQDGWeRBHbniifQO18xuBX4GPOvuTwJPZhzrQRB8E1sjV6QavfXWW4wcOVLBNiYzo2fPngwdOpTXX3+96gJuVIVe5EGTfbwraqepnsDDQDcz+wXB4u9PufuusCm5LtxEJKb9+/fTpUuXYmej7PXo0YN9+/YVOxtlIUqtuNDrEkPl1ogjBVx3/xzwOTM7Ffg4cAPwkJn9mqDm+3N1ohLJn1l7S0NLLvRvWDiFHtoUp0ZcDkE61rAgd38OeA74hpm9lyD4TgN+ZGZ/IAi+D7v7hqQyKiIi5SdOjbhZlBpxOXTkynsuZXd/CbgVuNXMBgDnAB8LD3833+uLSDwX3hW83Xn0iolFzolUu6g14ijTX0L5dORKdPGCsDn53nATkSqzZcsWrr76anr27Mldd91V7OxIGYoztKlcOnLF6gZpZjeY2YxW9s8ws3/LP1siUo6uu+46hg0bxoMPPpj4tWfPns3RRx9N9+7dGT9+PEuXLk38HlIaxo/oy+f//r2RO3LNPGtMTsGztWbrNMQdd3AxwUQX2RqAz8TPjoiUq6amJh577DEmT57MkUcemei1H330Ua666iquvfZaVqxYwWmnncbkyZNZv359oveR8hUlSDc3W3c2cm62TkLcgDsQaK038pu8OwuViFSRuro6jjnmGBoaGjj99NMTvfZtt93GtGnTuPzyyznuuOP4wQ9+wJAhQ/jRj36U6H2kOkStESclbsBdD0xqZf8HgNfiZ0dEkrJ99z42vLWLhnVbU7nf0qVLmTBhAgsWLGDKlCktjt9444306tWr3a21ZuK9e/fS0NDAWWeddcj+s846i2XLlhWsPFLZojZbJyFup6m7gO+ZWVfg6XDfh4CbgG8nkTERia9h3Vb+vHE7Bxym3rM8lb/i6+vrGT16NGvWrGHy5Mktjs+YMYMLLrig3WsMHTq0xb4tW7awf/9+Bg06tPFs0KBB/OpXv8ov0yIpijsO91Yz6w98H+ga7t4L3OHu30kqc6KhHRJPnNl98rV27VpefPFFZs6cSU1Ny18t/fr1o1+/frGvnz2RhbtrcgspK7Ena3X3rwP9gQnARGBA5gIHIlI8E0YdQacwFqXVKaSxsZEdO3ZwxRVXtHo8bpNy//796dy5Mxs3bjxk/6ZNm1rUekVKWV7jcN19p5m9GH6/I5ksiUi+xo/oy7GDD+ft3U3ccdHJqbynqqmp4YYbbqBnz56tHo/bpNy1a1fGjx/P4sWLOf/88w/uX7x4Meeee25+mRZJUeyAa2ZfBmYCQ8PPrwO3Abe7uyeTPRGJ6/DuXTi8e5dUgu3ChQtZv3493bp144UXXmDFihVMnTr1kDT5NCnPnDmTiy++mFNPPZXTTz+dO++8k9dff50ZM1pMByBSsmIFXDP7DjAduIV3VweaCMwChgBfSyR3IlLy9u/fz7x583jggQeYNWsWQ4YMYe7cuYne48ILL+TNN9/kW9/6Fm+88QZjx47liSeeYMSIEYneR6SQ4tZwLwMuc/fHM/Y9bWarCXowK+AmZPvufby9u4mGdVtLcjJukc6dO/P448GvguxabZKuvPJKrrzyyoJdX6TQ8lnh+o9t7NOq2QlpHtrx2tZdTL1neWrjKUVEJHlxa7j3A58Hrsra/zlgfl45koOKMbRDKoeGkomUlrgBtxvwKTM7G1ge7nsfcCTwoJl9vzmhu38pvyxWr+ahHQc83fk+RUQkeXGbf48F/ht4AxgRbhvDfccBJ4Tb2FwvaGYfMLMFZrbBzNzMpmUdNzO73sxeN7NdZrbEzI7PStPXzOab2bZwm29mfbLSnGBmz4TX2GBms6xER883D+04qm+Pkl1QWUREchN3pqm/TzojQC9gFUFz9f2tHP8a8BVgGrCaoEf0YjMb4+7bwzQPAcOByYAD9xA0cZ8DYGa9gcXAb4FTgDHAPGAncGsBypS3NId2iIhI4SS6AH0+3P0J4AkAM5uXeSysgX4ZuNndfxLuuwTYBHwKuMvMjgM+Arzf3ZeFaa4AloZBeTUwFegJXOLuu4BV4Xkzzew2jR8WEZFCiRRwzWxBLunc/WPxstOmo4HBwFMZ99hlZr8FTiMYijQR2AFkLh/yLEHt9TSCWvFEYGkYbJstAv4NGAm8knC+RUREgOg13H8E1gFLks9KuwaHXxuz9jcSznQVptmcWUt1dzezTRnnD6bl8oGNGcdaBFwzm04wyQfDhw+Pm38REalyUQPud4FPE6x7ex8wz93TXP82u8nXsva11iTcURprY3+w0/1u4G6A2tpaNTmLiEgskXopu/vXgGHA1UAtsMbMfmlm55lZl0JkMNS8TMjgrP0DebeGuhEYmNnjOPx+QFaa1q4BLWvPJeHRKyZqPKXEc99Hg01ESkLkYUHuvt/dF7j7Jwjerf4G+Bawwcx6JZ3B0CsEwfLM5h1m1h2YxLvvbOsIejpnRqeJwGFZaSaF5zY7E3gdWFuIjItUky1btnDxxRe3uUSfSDXLdxrGw4A+BIFuB200y+bCzHqZ2TgzGxfma3j4eXj4XvZ24Bozm2JmYwmG8+wgGAqEu78IPEnQY3mCmU0k6Ey1MOyhTJj2HWCemY01synANYB6KIsk4LrrrmPYsGE8+OCDiV73pptu4pRTTqF3794MGDCAc845h1WrViV6D5FCixxwzayHmV0S9hB+nmDSi0vcfZS778wjL7XAinDrAfxL+P2/hse/Q7D83w+BeoJVic7KGIMLwbCf/yHozbwo/P7i5oPuvo2gRntkeI0fEoy/vS2PfIsI0NTUxGOPPcbkyZM58sgjE732kiVLuPLKK1m2bBlPP/00NTU1fPjDH+avf/1rovcRKaSow4LuBi4E1gBzgY+5+1tJZMTdl/BuB6bWjjtwfbi1leavBJ262rvP8wSdvkQkQXV1dRxzzDE0NDRw+umnJ3rtRYsWHfJ5/vz5vOc97+HZZ5/lnHPOSfReIoUStYZ7GbCVYErHycD94XSMh2yJ51JEotvzNmx7FV59LpXbLV26lAkTJrBgwQKmTJnS4viNN95Ir1692t2WLl2a0722b9/OgQMH6NtXM7BJ+Yg6LOh+8nhPKyWquSfrpb8obj4kOa8+B42rwA/Ajz8GlyyAYacW9Jb19fWMHj2aNWvWMHny5BbHZ8yYwQUXXNDuNYYOHdru8WZXXXUV48aNY+JE9eCX8hEp4Lr7tALlQ5Ki4CkAa5cGwRZg/97gc4ED7tq1a3nxxReZOXMmNTUtf7X069ePfv365X2fmTNn8rvf/Y7f/e53dO7cOe/riaQl5yZlMzvVzHL+6Taz8QUemysibRk5CSz87925a/C5wBobG9mxY0ebQ4KSaFK++uqrefjhh3n66acZNWpUIYohUjBRarh1hNMn5pj+N8A44OWomRKRPA07FQaNhd3b4Nx7Cl67BaipqeGGG26gZ8+erR7Pt0n5qquu4pFHHmHJkiUce+yxeeVVpBiiBFwDbjKzd3JM3zVGfqQcqNm6PHTrHWwpBNuFCxeyfv16unXrxgsvvMCKFSuYOnXqIWnyaVL+/Oc/z/z58/nZz35G37592bgxmHyuuWYsUg6iBNzfAv8nQvo6YFeHqaTyKUBXtP379zNv3jweeOABZs2axZAhQ5g7d26i95g9ezYAH/rQhw7Z/81vfpPrr78+0XuJFErOAdfd/66A+RB5lwJ0WencuTOPP/44QItabVI0EZxUgnyndhQREZEcRB2HKyLlQi0EIiVFNVypDFqKTkRKnAKuiIhIChRwJfU5d0VEqlHsgGtm3czsaDP7GzMbkGSmJEXNc+6+tS6YczeXoBs1QBc6vYhIGYgUcM3scDP7XLgW7jbgJWAVsNHMXjWzOWZ2SiEyKgXS2py77YkaoAudXkSkTESZS/lqYC3wWWAx8HGCqRuPASYC3yTo9bzYzJ40s9GJ51aSF3XO3agButDpm6lWLCIlLsqwoNOAD7r7qjaOPwfca2YzgH8CPkiwUL2Usqhz7jYHaD+QW4AudHqItxSdJtcQkZRFmWnq/ObvzWy4u69vI90eYHYCeZM49rwdBM9Xn8t9Dt0oc+5GDdCFTg9FWYquHFz65KUA3PeR+4qcExGB+J2mfmpm3Vo7YGbd88iP5COt95/desN7hkUL6IVMX4Sl6KR1W7Zs4eKLL25ziT6RahY34L4E3J2908yOBHJ86SY5iTKhQ9z3n+WuuVbcZ0RuzclSMNdddx3Dhg3jwQcfLNg9brzxRsyML3zhCwW7h0ghxA24nwXGm9kXm3eY2TiC97h/SSJjEkM11/Si1orVySpxTU1NPPbYY0yePJkjjzyyIPdYvnw5c+bM4cQTTyzI9UUKKVbAdfd3gHOBb5rZ+83sEwQ123vd/aIkMygRqKaXGw09Koi6ujqOOeYYGhoaOP300xO//rZt25g6dSpz586lb9++iV9fpNCiDAtaZGbfNrOLzOxY4H+B6cBC4AFgurvPKlA+JVdRa3rVKE7TexnO1bxj7w7e2PkGKzetTOV+S5cuZcKECSxYsIApU6a0OH7jjTceXDC+rW3p0rafxfTp0znvvPM444wzClkMkYKJMixoBcG4288Ag4B3gOeB/cB/AP9rZt3CXsoi8RV6qE6coUdlZuWmlazeuhrHufypy5lz1hzGDRxX0HvW19czevRo1qxZw+TJk1scnzFjBhdccEG71xg6dGir++fMmcNLL73E/PnzE8mrSDFEGRZ0TfP3ZjYIOJkgAI8DPkDwXveAma1x9+OTzqiUsTTGuka5R5yhR2WmvrEeJ1i0fd+BfdQ31hc84K5du5YXX3yRmTNnUlPT8ldLv3796NevX+Trrl69mmuvvZalS5fStWvXJLIqUhRx3+E2uvuT7n6zu1/k7scChxNMdvH9RHMoUggV3smqdlAthgHQpVMXagfVFvyejY2N7Nixo80hQXGblOvq6tiyZQtjx46lpqaGmpoannnmGWbPnk1NTQ179qhRTcpDYgvQu/suoC7cpJJV2+xMZTiT1biB4xjTdwzb923n5kk3F7x2C1BTU8MNN9xAz549Wz0et0n5E5/4BLW1h/7BcOmllzJ69GiuvfZa1XqlbOQccM3saHd/Jce0Bhzl7q/GzplIqSjTmax6de1Fr669Ugm2CxcuZP369XTr1o0XXniBFStWMHXq1EPSxG1S7tOnD3369Dlk32GHHUa/fv0YO3ZsXvkWSVOUGm6dmf0CuMfdW63Fmllf4CLgS8APgX/PP4tScNVWY4VoZY7TySrOFJtlav/+/cybN48HHniAWbNmMWTIEObOnVvsbImUnCgB91jgG8AvzGw/0AC8AewG+gJ/AxxHMPnFl919UcJ5FSmOqJ2syrAJOh+dO3fm8ccfB2hRqy2UJUuWpHIfkSTl3GnK3d9y968CQ4HPAX8G+gBHA03Aj4GT3f10BVupOFE6WVXrFJsi0q7InabCzlGPh5uIZMu3CZrDE8mGVgkSKS1x51I+yMyON7PLzOx9Wft7tbWiUB73ut7MPGvbmHHcwjSvm9kuM1tiZsdnXaOvmc03s23hNt/M+rS8m0iGS3+Re3Nv1Ck2s6eabMpxmMuWNcEmImUh74BLML3jqcBZZvZTM/uamfUEugGF6DmxGhiSsZ2QcexrwFeALwKnAJuAxWaWWWV4CPhbYDLwkfB7TV8jycqnCTrXgHtgf5B+7874+RSR1CQRcK8lmFd5AkET9ccJguK1QCEGyDW5+8aMbTMcHIr0ZeBmd/+Ju68CLiFon/tUmOY4giA73d2Xhb2trwD+0czGFCCvIh3LXuWpphvu3v45e3dC064g4G55SUG3FR3+G4qkLO+JL9x9J/Bd4Ltm1hUYAwwOt0L8xI8ysw3AXuD3wLXu/jJB563BwFMZedtlZr8FTgPuAiYCO4BlGdd7FtgZplldgPymqwx7uVakPKaa7LyrJ/v27Wt/Qoc92zM+HAg+dz2s/fs0Nz/3H5173srYrl276NKlS7GzIXJQ5BqumZ1kZtPM7B/CAHuQu+919+fdfbG7zydo4k3S74FpBM3BlxME2GVmdkT4PUBj1jmNGccGA5s940/f8PtNGWkOYWbTzazezOo3b96cVDlyV2ZTCkpMGU3Qffr0obGxkQMHDrSTPvMtSaesz22okiZod+edd95hw4YNDBw4sNjZETkoUg3XzKYDP4JwklZYY2ZnuPuG1tK7ezu/MaJz919m5Wc58DJB0/Hy5mRZp1nWvtZq3dlpMu95N3A3QG1tbbptVHHGc0p5yqgR9+/fn9dee43VqztocNn+VvCz0fMI2La+/bRNe2DHJsDBNsBhA6Gmgz6NOzYFX3uVX9Dq0qULgwYNonfv3sXOishBUZuUvwbMBv4NGAbcDnwb+HTC+cqJu+8wsz8Bo4GfhbsHA5lTSg7k3VrvRmCgmVlzLTd89zuAljXj4ivTKQUlP506dWL48OEdJ7zvn4OvuTRfL70Vfv2vwffWGc74Bkz6Svvn3Dm9oldUEklb1CblEcB33X2TuzcQNO+2XGk6JWbWnWAGrDeAVwgC6plZxyfx7jvbOqAXwbvcZhOBwzj0vW5pyO5MU4HrtkoeogxVivqzlD1UKZdXGvd99N0Zs0SkhagBtzOwq/mDu/8FwMyGJJmptpjZd83sg2Z2dDju93GCYPnjsMZ6O3CNmU0xs7HAPIJOUg+F+X0ReBK4y8wmmNlEgs5UC9299DpMRR3PKdKWqD9LcWbLUn8DkXbFGRY03czOMLPmZT/2Az0SzFN7jgIeJuhN/FNgDzDB3deFx78D3EawcEI9wTjds9w9s0vnVOB/CHozLwq/vziV3McRdd1WkbZE+VlKo0YMqhVLVYn6DncJMBP4V8DN7HWgO0EQ/jVQ7+5bk83iu9z9og6OO3B9uLWV5q8U6Z2zSFHlMVQpVo04l8BeRasqiUQKuO5+BoCZjQLGZ2yXEXSocjN72d2rY6CfSCWbEWHRhTjzR6sXvlSZWBNfhBNNvAw81rzPzEYCtQRTJYpINYlaI4Z4teI7J6ls+EieAAAQLklEQVTntJStvGeaaubua4G1aBUhkeoUpUYM0WvFqhFLmUss4IqIRDLsVPjsoqBmO3JSYd4TN3fI0pSnUgIUcEWkeIadmnstNd91hlUbliKLtVqQmeUwDY6ISILyXWdY44OlyOLWcH9qZqe7e4uFO82su7vvzjNfIiItRXlPrE5ZUmLirof7EuGE/pnM7EggYs+J6nLhXXVceFddsbMhUvnSmrxDJEdxA+5ngfFm9sXmHWY2DngO+EsSGRMRyUtzp6wPzSrcdJYQ1IpvP1EBWjoUdxzuO2Z2LvCsma0A+gPzge+5+6wkM1jtLrVgEaP7ipwPkbJU6E5ZGqokEeQccM1sEbASWBF+XQ1MBxaG17nc3R8uRCYryfbd+3h7dxMN67YyfkTfxK9/6ZOXAnDfRxSiRSKJOkwJ9J5YIonSpLwCOAn4HvACsB34KsHiBQ8B/2tmHaxoXd0a1m3lzxu389rWXUy9ZzkN6wo27bSIxDHs1GCd4Ki1YtB7YulQzjVcd7+m+XszGwScDIwLtw8QvNc9YGZr3P34pDNaCZa//CYHPPh+X9MBlr/8ZkFquSKSkjQm71CNuGLEfYfbSLCu7JPN+8ysB0HwPTGZrFWeCaOOoJPBAYcuNZ2YMOqIjk8afELhMyYi8RXyPbHeEVeUJOdS3gXUhZu0YvyIvnynx/08v384H/vstardilQb1YirmqZ2TNlxNRs4rmYDx5dQsFVHK5EUqUZctRRwU3b8kPcUOwsiUi7SqBFLahRwRURKWaHHEoOaoVOigCsiUinijCVWM3RqFHAlMr3zFSlhUWrEEK8Z+tXnogV1ARRwU6epGkWkpIycBDU9gmCba8ese88OgnRND9WII1DAlYJTjVikhA07NQia6phVcAq4UnIUoEVSVuiOWeqUBSjgpu7/7N7Fsbt3Bc0yVfyDJyJlKmrHLHXKOkgBN02vPsdXGzdQ4171P3hJU61YJEVRasRqgj4o7gL0EsfapdS40xmiLXAtibv0yUsPBmkRKaCoKypB0AR9+4kVt5qSarhpGjmJJjNwp3OUQelSdKpBi8SkJuiDFHDTNOxUbhk0lGN37+K8cx+pmB8iaUkBWiRDGk3QZTA2WAE3ZX/p3oO/dO/BeSX6AyHFESdARz2n1NKLtCrquGAom7HBCrgiUjSl9kdAtf7hU1J5ijouGMqmY5Y6TYlIInbs3cEbO99g5aaVBTun1NKXYp4qoQznr7iFj2xazMpuXXNKz8hJ/KF7N+5+T29W9uiZU634y4+cxVX3jmflqodyu0cCFHArTCX8Z6vGMpR7mVduWsnqravZsGMDlz91eUHOKbX0pZinqi1Dt65cPmQw/963D5cPGdhhoF656iF+t2sDSzrt4fI/3Jha0K3agGtmV5rZK2a228wazKwkuwzrl17ll6ESylzfWI/jAOw7sI/6xvoOyxD1nFJLX4p5quYy7Mdxg31+oON7vLyIfWYcMGOfBZ/TUJUB18wuBO4AbgROBpYBvzSz4YW+9w4O8AZNBfmlWs3/2UopT9VY5tpBtRgGQJdOXagdVNthGaKeU2rpSzFPKkOO9xh1Nl3c6exOFw8+p6FaO03NBOa5+5zw8xfN7CPA54CvF+qmKzetZDX7cODypy5nzllzGDdwXJvpW/ul11765h86xyP/oOZ6TqmlL8U8VWOZxw0cx5i+Y9i+bzs3T7q53Z/TuOeUWvrmc+6ffD/1jfXUDqrN6R6llL4U85RKGcZ+insIara1o85m3NhPdXiPRLh7VW1AV6AJOD9r/w+BZ9o7d/z48Z6POX+c42PvG+tj5431k358ks/545x2069oXOEnzDvBx84b67Xza31F44oO73Hef53nZz9+dk5p455TaulLMU/VWOZpv5zm0345Lae0+ZwjUmqAes8h/lRjDbc/0BlozNrfCHw4O7GZTQemAwwfnl+Lc1BrAHfo0rkwtYZeXXvRq2uvnNLGPafU0pdinqqxzHHG32rMrlSTagy4zTzrs7WyD3e/G7gboLa2tsXxKMYNHMeYfscWPICKtCZqcFMwFElWNQbcLcB+YHDW/oG0rPUmTgFUWqPaoUjlq7qA6+57zawBOBN4LOPQmcBPipMrqTQKhiKSreoCbug2YL6ZPQc8C8wAjgTuLGqupGQpgIpIvqoy4Lr7o2Z2BPD/gCHAKuAf3H1dcXMmaVEAFZG0VWXABXD32cDsYudDRESqQ9UGXKksqrGKSKlTwJWSo+ApIpVIAVcKTgFURKRKFy8QERFJm2q4EplqrCIi0amGKyIikgIFXBERkRSoSTllpdgcW4p5EhGpNKrhioiIpEABV0REJAUKuCIiIilQwBUREUmBOk1VGHWAEhEpTarhioiIpEABV0REJAVqUi5xaiIWEakMquGKiIikQAFXREQkBQq4IiIiKVDAFRERSYECroiISAoUcEVERFKggCsiIpICBVwREZEUKOCKiIikwNy92HkoG2a2GVgX4ZT+wJYCZadUqczVoRrLDNVZbpW5YyPcfUBHiRRwC8jM6t29ttj5SJPKXB2qscxQneVWmZOjJmUREZEUKOCKiIikQAG3sO4udgaKQGWuDtVYZqjOcqvMCdE7XBERkRSohisiIpICBVwREZEUKODmyMyuNLNXzGy3mTWY2aQO0n8wTLfbzF42sxn5XjNtSZfZzK43M8/aNha2FNFFKbeZDTGzh8zsz2a238zmtZHuXDN7wcz2hF//b8EKEEPSZTazaa08azez7gUtSAQRyzzFzJ4ys81mtt3Mfm9mH2slXSU95w7LXIHP+YNmtszM3jSzXeHP+D+3ki7ec3Z3bR1swIXAPuBy4DjgB8AOYHgb6Y8GdobpjgvP2wecG/eaFVLm64E/A4MztgHFLmue5R4JfB+YBiwD5rWSZiLQBHwjvOY3ws/vK3Z5C1jmaeHPQ+azHlzssuZR5juAa4BTgfcC3wT2A5Mq+DnnUuZKe87jgYuA48PfaZ8Oy3dlEs+56P8g5bABvwfmZO1bA9zURvpvA2uy9t0D1MW9ZoWU+XpgVbHLlmS5s9ItbCP4PAosztr3K+DhYpe3gGWeBuwodtkKUeaM9M8Bt1bDc26nzNXwnH+a+Qzzec5qUu6AmXUl+KvnqaxDTwGntXHaxFbSLwJqzaxLzGumphBlztg3ysw2hE08j5jZqEQynYACPpe2/m3K9VnnqoeZrTOz18xsoZmdnOf1EpFgmQ8HtmZ8robnnF1mqODnHJblNOCZjN2xn7MCbsf6A52Bxqz9jQTNJ60Z3Eb6mvB6ca6ZpkKUGYK/NqcBkwmaeAYDy8zsiPyznIhCPZe2/m3K9VnnYjXwWeDjwCeB3cCzZjY6j2smJe8ym9nngaOA+Rm7K/o5t1HminzO4R8Pe4B6YLa735lxOPZzrukogRyUPWDZWtnXUfrm/dZOmlIaGJ1kmXH3Xx5y0Gw58DJwCXBb/GwmrhDPpdKedfsXc68D6g5ezGwZsBL4IvCluNdNWKwym9m5wC3ARe6evZhJRT7ntspcwc95EtALmAB828xecffMPzRi/Tsq4HZsC0FHgey/XgbS8q+cZhvbSN8EvEnwcKJeM02FKHML7r7DzP4ElMJfwxCv3Llo69+mXJ91ZO6+38zqKY1nHbvMYeCZD3zG3RdkHa7I59xBmQ9RKc/Z3V8Jv33ezAYR9D9pDrixn7OalDvg7nuBBuDMrENnEvTQbE0d8OFW0te7+76Y10xNIcrc2gnh0IFjgTfi5zY5BXwudQW4ZiLS+lk0MwNOpASeddwym9kFwAPANHd/vJUkFfeccyhzdvqyf86t6AR0y/gc/zkXuxdZOWwEXcv3ApcRdAO/g6Br+Yjw+P3A/Rnpm4fI3B6mvyw8P3tYUJvXLPZWoDJ/F/hgmPZ9BD1c3y6VMscpd7hvXLj9FlgQfv83GcdPI6jpf53gD4yvEwxVKKXhIkmX+ZvA2cCo8Ni9YZlPLXZ5Y/58XxTm/yoOHQLTr1Kfc45lrrTn/EXgHwlq6KOBfwp/R92cxHMu+j9IuWzAlcBaYA/BX00fyDi2BFiSlf6DwH+H6V8BZkS5ZilsSZcZeAR4PfwPsAH4CRm/pEtli1Fub2Vbm5XmPIIxyHuBF4EpxS5nIcsMfA9YF15vE0EvzonFLmfcMoefWytz9r9LxTznXMpcgc/5y8CfCCoP2wh+n10JdEriOWvxAhERkRToHa6IiEgKFHBFRERSoIArIiKSAgVcERGRFCjgioiIpEABV0REJAUKuCIiIilQwBUREUmBAq5IGTOzb5vZ4hzTLjEzD7cJca9TLGY2LyP/5xU7PyJRKeCKlLdxBMuh5eo+YAjBFHfZ1/mf9k40szvN7Hvh983B755W0n0nPLYwQr5ycRVB3kXKkgKuSHk7CVgRIf077r7RW67gdBLtBNxwFZhzgP/K2P0qcKGZHZaRrga4GFgfIU85cfdt7r4x6euKpEUBV6RMmdlgYBBhDdfMDjOzR8zsv81sZIzr7DWzJ8xsp5n9xcz+PiPZKUB34HcZ+/4IrAEuyNj3UWA3waTwmfdYEtaQ7zCzreF2i5l1ykhjZvYVM1tjZnvM7DUzuynXcoiUOgVckfJ1MrALWG1mY4DnCJYNO93d10a8DsDnCVZ/OQlYBdyWkeYTwC/cvSnr3LnAZzM+f5ag2bq1VVGmEvzOmQhcAUwnWJ2l2Y3AdcBNwPHA+QS1aJGKoIArUr7GAc8TBMNlwBx3/7S774pxnW3ABe6+2N1fAh4HBmSk+Tjws1bOfQioNbPRYU35I8C8Nu7zBvAld/+zu/8HcAswE8DMegFXA9e4+73u/pK717n77IhlESlZNcXOgIjENo5gkex7gY+5+zN5XOfnWe9H3wu8BGBm7yVYYHxR9onuvtXM/pOgZvsWwdqi64NXvi0s90PXA60D/s3MehMs5N0N+HXMMoiUPNVwRcrXOOCnQBfgiDyvU5e172Te7f38CeDX7r6zjfPvBT5DEHTvjZmHViO0SCVRwBUpQ2bWk6AWehdwGXC/mf1tHtfJ7umcGXDbak5u9mtgL9C/g3Tvs0OrvhOA1939beAFYA/wodxzL1Je1KQsUp5OIuiYtMrd/2BmxwE/N7NT3X1DxOtA0OMYADM7AjgKWGlmAwgCY5sTTbi7m9mJgLn7nnbudSRwu5nNBk4Avgp8K7zGdjO7A7jJzPYAvyWotY939x9FKI9IyVLAFSlPJwFrMjpIzQLGAAvMbJK7vxPxOpnNxScD+whqnZ8G/uDuje1dxN2353CvB4HOwO8J/liYS9ArutnXga0EPZWPAhqB+3Mrhkjps0P7MIhIpTKzJQQ14i9EOOe/gGfd/Ttp37udazlwvrs/nu+1RNKkd7gi1WW6me0ws1NyTP8s8HAhM5SrcOKMHcXOh0hcquGKVAkzGwr0CD++2sH71qTvvYQ8a7hmNhDoHX58o51e0yIlSQFXREQkBWpSFhERSYECroiISAoUcEVERFKggCsiIpICBVwREZEUKOCKiIikQAFXREQkBQq4IiIiKfj/SNhSsljIRvwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7,4))\n",
    "for i in range(3): \n",
    "    plt.errorbar(kd, kd * bpk[i], yerr=kd*err[i], fmt='.', label=r'$\\ell=%s$' % (2*i))\n",
    "plt.xlabel(r'$k \\ [h/{\\rm Mpc}]$')\n",
    "plt.ylabel(r'$k \\ P_\\ell(k) \\ [{\\rm Mpc}/h]^2$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating the fake data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to create our fake data. For that, we can use the build-in function in the `ReadWrite` class. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "io = ReadWrite()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data structure is a dictionary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each sky corresponds to an entry in this dictionary, which is itself a dictionary containing the data. \n",
    "\n",
    "Here we have only one sky, and we name it 'sky_1'. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "d['sky_1'] = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now start filling this dictionary with the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['min', 'max', 'eff'])\n",
      "dict_keys(['Omega_m', 'H', 'D'])\n"
     ]
    }
   ],
   "source": [
    "zmin, zmax, zeff = 0.43, 0.7, z\n",
    "io.write_common(d['sky_1'], zmin, zmax, zeff, Omega_m_fid, H_fid, D_fid)\n",
    "print (d['sky_1']['z'].keys())\n",
    "print (d['sky_1']['fid'].keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['multipole', 'x', 'cov', 'nsims', 'survey_mask', 'binsize', 'l0', 'l2', 'l4'])\n"
     ]
    }
   ],
   "source": [
    "mult_pk = 3\n",
    "io.write_pk(d['sky_1'], mult_pk, kd, bpk, cov, nsims_cov_pk=0, \n",
    "            survey_mask_arr_p=None, survey_mask_mat_kp=None,\n",
    "            binsize=dk)\n",
    "print (d['sky_1']['bPk'].keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we are not providing a survey mask. One can compute the survey mask from random catalogs using https://github.com/pierrexyz/fkpwin , in a format compatible with PyBird: \n",
    "```\n",
    "survey_mask = np.load('output_fkpwin/Wlmkp_skyname.npz')\n",
    "survey_mask_arr_p = survey_mask['p']\n",
    "survey_mask_mat_kp = survey_mask['Wlm']\n",
    "```\n",
    "Also, there are two other build-in functions that we are not using here, feel free to explore!\n",
    "    \n",
    "```\n",
    "io.write_cf(d_sky, mult_cf, ss, cf, cov_cf, nsims_cov_cf=0, binsize=None)\n",
    "io.write_bao_rec(d_sky, bao_rec_rd_fid, bao_rec_H_fid, bao_rec_D_fid, alpha_par, alpha_per, \n",
    "                 cov_alpha, cov_cross_pk=None, cov_cross_cf=None)\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we save our fake data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(os.path.join(datadir, 'fake', 'fake_example_2pt.npy'), d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing the fake data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool. Now let us test the fake data by maximazing the likelihood: do we recover the truth?\n",
    "\n",
    "We first need a configuration file. Let us create one. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "configdir = os.path.join(datadir, 'eftboss', 'config')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity, let us load a configuration file and modify it accordingly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "config = yaml.full_load(open(os.path.join(configdir, 'boss_cmass_ngc_pk_wc.yaml'), 'r'))\n",
    "config[\"data_path\"] = os.path.join(datadir, 'fake')\n",
    "config[\"data_file\"] = 'fake_example_2pt.npy'\n",
    "config[\"multipole\"] = 3\n",
    "config[\"with_survey_mask\"] = False\n",
    "config[\"with_binning\"] = True\n",
    "\n",
    "config[\"sky\"]['sky_1'] = config[\"sky\"].pop('cmass_ngc')\n",
    "config[\"sky\"]['sky_1']['min'] = [0.01, 0.01, 0.01]\n",
    "config[\"sky\"]['sky_1']['max'] = [0.20, 0.20, 0.20]\n",
    "\n",
    "config[\"write\"][\"fake\"] = True\n",
    "config[\"write\"][\"out_path\"] = os.path.join(datadir, 'fake', 'out')\n",
    "config[\"write\"][\"out_name\"] = 'sky_1'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We center the prior on the truth for this check. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, truth in eft_truth.items(): \n",
    "    if name not in ['b1', 'b2', 'b4', 'c4']: \n",
    "        config[\"eft_prior\"][name]['mean'] = [truth]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data_path: /Users/pierre/Documents/github/pybird/data/fake\n",
      "data_file: fake_example_2pt.npy\n",
      "output: bPk\n",
      "multipole: 3\n",
      "with_ap: true\n",
      "with_survey_mask: false\n",
      "sky:\n",
      "    sky_1:\n",
      "        min:\n",
      "        - 0.01\n",
      "        - 0.01\n",
      "        - 0.01\n",
      "        max:\n",
      "        - 0.2\n",
      "        - 0.2\n",
      "        - 0.2\n",
      "write:\n",
      "    save: true\n",
      "    fake: true\n",
      "    plot: true\n",
      "    show: true\n",
      "    out_path: /Users/pierre/Documents/github/pybird/data/fake/out\n",
      "    out_name: sky_1\n",
      "km: 0.7\n",
      "kr: 0.35\n",
      "nd: 0.0004\n",
      "eft_basis: westcoast\n",
      "with_stoch: true\n",
      "eft_prior:\n",
      "    b1:\n",
      "        type: flat\n",
      "        range: null\n",
      "    c2:\n",
      "        type: gauss\n",
      "        mean:\n",
      "        - 0.5863\n",
      "        range:\n",
      "        - 2.0\n",
      "    c4:\n",
      "        type: unvaried\n",
      "    b3:\n",
      "        type: marg_gauss\n",
      "        mean:\n",
      "        - -0.3948\n",
      "        range:\n",
      "        - 2.0\n",
      "    cct:\n",
      "        type: marg_gauss\n",
      "        mean:\n",
      "        - 0.1839\n",
      "        range:\n",
      "        - 2.0\n",
      "    cr1:\n",
      "        type: marg_gauss\n",
      "        mean:\n",
      "        - -0.8414\n",
      "        range:\n",
      "        - 2.0\n",
      "    cr2:\n",
      "        type: marg_gauss\n",
      "        mean:\n",
      "        - -0.8084\n",
      "        range:\n",
      "        - 2.0\n",
      "    ce0:\n",
      "        type: marg_gauss\n",
      "        mean:\n",
      "        - 1.5045\n",
      "        range:\n",
      "        - 2.0\n",
      "    ce1:\n",
      "        type: unvaried\n",
      "        mean:\n",
      "        - 0.0\n",
      "    ce2:\n",
      "        type: marg_gauss\n",
      "        mean:\n",
      "        - -1.6803\n",
      "        range:\n",
      "        - 2.0\n",
      "with_binning: true\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print ( yaml.dump(config, indent=4, sort_keys=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us setup the likelihood with this configuration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading data file: /Users/pierre/Documents/github/pybird/data/fake/fake_example_2pt.npy\n",
      "-----------------------\n",
      "sky: sky_1\n",
      "output: bPk\n",
      "multipole: 3\n",
      "min bound (per multipole): [0.01, 0.01, 0.01]\n",
      "max bound (per multipole): [0.2, 0.2, 0.2]\n",
      "coordinate (AP) distortion: on\n",
      "binning: on\n",
      "-----------------------\n",
      "-----------------------\n",
      "get_maxlkl : False\n",
      "with_boss_correlated_skies_prior : False\n",
      "with_rs_marg : False\n",
      "-----------------------\n"
     ]
    }
   ],
   "source": [
    "L = Likelihood(config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We initialize the likelihood. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "eft_name = ['b1', 'c2'] \n",
    "eft_init = [2, 0.] \n",
    "init = L.loglkl(eft_init, eft_name, M, need_cosmo_update=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now maximize the likelihood. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min chi2: 0.00, p-value: 1.000\n"
     ]
    }
   ],
   "source": [
    "L.c[\"get_maxlkl\"] = True # minimization mode\n",
    "\n",
    "def chi2(eft_params): return -2. * L.loglkl(eft_params, eft_name, M, need_cosmo_update=False)\n",
    "minimum = minuit(chi2, eft_init)\n",
    "\n",
    "from scipy import stats\n",
    "def pvalue(minchi2, dof): return 1. - stats.chi2.cdf(minchi2, dof)\n",
    "dof = L.y_all.shape[0] - len(eft_name) - L.Ng # number of data points - EFT parameters\n",
    "print ('min chi2: %.2f, p-value: %.3f' % (minimum['fun'], pvalue(minimum['fun'], dof)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like we have recovered the truth ($\\chi^2 = 0$)! Let us see how well we have recovered it: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- ratio fit/truth -\n",
      "b1: 1.00\n",
      "c2: 1.00\n",
      "b3: 1.00\n",
      "cct: 1.00\n",
      "cr1: 1.01\n",
      "cr2: 0.99\n",
      "ce0: 1.00\n",
      "ce2: 1.00\n"
     ]
    }
   ],
   "source": [
    "print ('- ratio fit/truth -')\n",
    "for name, truth, fit in zip(eft_truth.keys(), eft_truth.values(), L.out[0]['eft_parameters'].values()):\n",
    "    if name not in ['b2', 'b4', 'c4', 'ce1']: # those are not in the westcoast prior\n",
    "        print ('%s: %.2f' % (name, fit/truth))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build-in fake data generator from best fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By the way, did you know that there is a build-in function that can create fake data directly from a fit?\n",
    "\n",
    "This happens if we set: ```config[\"write\"][\"fake\"] = True``` and run the following function: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fake data from best fit saved to /Users/pierre/Documents/github/pybird/data/fake.\n",
      "data files with best fit saved to /Users/pierre/Documents/github/pybird/data/fake/out.\n",
      "data plot with best fit saved to /Users/pierre/Documents/github/pybird/data/fake/out.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEmCAYAAADYym9WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4VFX6wPHvOyUF6YSEEAKB0IuiICJNRUAQUBQVEVDZVbD+UFzL2tZdFVjrqouKgGKoK6iIUVAQlRaEIFhBSmgmkBB6SWAyc35/3ElIQgLJZJJJwvt5nvvMzL3nnvuGMm/OueeeI8YYlFJKqfLIFugAlFJKqcJoklJKKVVuaZJSSilVbmmSUkopVW5pklJKKVVuaZJSSilVbmmSUqqcEZHvRGRToONQqjzQJKXUeU5EqorIP0XkSxHZJyJGRJ4IdFxKgSYppRSEAc8C7YD1AY5FqTwcgQ5AKRVwe4AoY0yKiMQA2wMbjlKnaUtKqTLm7V57RUS2i0imiKSKyFIRufIs53QTkcPeLrkQEVkpIj8VUnatiCQWNR5jzEljTIoPP4pSpU6TlFJl7x3gQeBT4H7gZWA/0L6gwiJyNfAV8A0wyBiTCUwDLhSRC/OVbQ50BKaXVvBKlSXt7lOq7A0AJhtjxp6roIj0B+ZhJbTbjTFZ3kMfAW8Cw4HHcp0yAsgCZvs1YqUCRFtSSpW9I0AnEYk6WyERuQkrOc0ChudKUBhjDgPzgdtExOYtL8Aw4GtjTFppBa9UWdIkpVTZexxrJN0uEUkUkRdEpFW+Mg2BOcAXwF3GGE8B9XwIRAFXeT93BRqjXX2qEtEkpVQZM8bMAZoA9wE7gf8DfhaR4bmKpQLLgT5Yyacgi7FG5mWfNxw4CnxWCmErFRCapJQKAGPMHmPMJGPMYKxW0zbgH7mKnAQGAj8DX4hIxwLqcAMzgMEiUhO4BZhnjMko9R9AqTKiSUqpMiQidhGpkXufMeYQ1rNJtfLtPwb0A5KAr0SkbQFVfghUAyZ5z48rjbiVChQd3adU2aoGJIvIx8BPWIMougJ9gYn5CxtjDolIH+B7YLGI9DDGbMl1/DcRWYfVitrtLVdsIvIAUNO7AVwlItnfD295B2ooVebEGBPoGJQ6b4hIEPAC0BvrvpQDqxX1AfCGMSZLRL4D6hljWuY6rz6wDAgCuhtjduY69iDWcPTxxpgnfYxrB9CokMONjTE7fKlXqZLSJKVUBScio4F3gdbGmI2Bjkcpf9IkpVQFJyJrAWOM6RToWJTyN70npVQFJCIXYI3+6441DdKtBZSpAYSeo6p93lGCSpVL2pJSqgLKNVv5IWCKMebRAspMA+44R1V6v0mVa5qklKqkRKQ1UP8cxVZ4J6xVqlzSJFUMYWFhJiYmJtBhKKVUhbJu3bp0Y0xdX87Ve1LFEBMTQ2JikZfpUUopBYjIznOXKpjOOKGUUqrc0iSllFKq3NIkpZRSqtzSJKWUUqrc0iSllFKq3NIkpZRSqtzSJKWUUqrc0iSllFKq3Co3SUpE/i4ia0XkiIjsE5HP869EKpbnRCRFRDJE5DsRaZOvTC0RmS4ih73bdO/S2rnLtBOR7711JIvIsyIiZfFzqvJjyKQEhkxKCHQYSqmzKDdJCrgSeBvoAvQEsoAlIlI7V5nHgEeAB4FLgTSs1Uqr5SozC7gEa9ntvt7307MPikh1YDGQ6q3j/4BHgbGl8UOp0uGPBHM000XyoQzW7Tzop6iUUv5WbqZFMsZck/uziIwADmMtrf25t6XzEDDBGPOxt8wdWInqNmCSiLTCSkzdjDGrvGVGA8tFpIUx5g9gGFAFuMMYkwH86j1vrIi8ZnQywwrhaKaLI5lZrNt5kA6NahVYxu12s3//ftLT00lPT2ffvn05r5v3u/jdcTHYHAybvJqZd3cutB6lVOCUmyRVgGpYLb3sX3MbA/WAr7MLGGMyRGQZVutrEnA5cAxYlauelcBxb5k/vGWWexNUtq+A54EYrOUPcojIKGAUQMOGDf3zk6kSWbfzIBv3HMEYGPLOCq7iF9i/PScJZSeigwcPUtjvHOFXDCekUwdEhIxTLl7+cD7vPTSY6tWrl/FPo5Q6m/KcpN4ANgDZfTr1vK+p+cqlAlG5yuzL3RoyxhgRSct1fj3gzwLqyD6WJ0kZY94D3gPo2LGjtrICyBjDkiVL+Nv7i/FEd0NsdlxuN/NW/Uzo9mWEhYVRt25dLrroopz3Bb2GhYXxy57j3PzuKjweg814+HTSSyx89SHuuecexowZQ/3651rhQilVFsplkhKR14BuWN12+VcNzZ8oJN++ghLJucpIIftVOWCM4YsvvuCFF17ghx9+IOqi7jiiuwGGkOAgPpn5Nh1iap+zntw6NAqiZb1qHMnM4o1bL0aGzuTll1/mlVde4fXXX2f48OH87W9/o3Xr1qXzQymliqQ8DZwAQEReB4YCPY0xSbkO7fW+1st3SjinW0J7gfDcI/W87+vmK1NQHXBmK00FkMfj4eOPP+aSSy5h4MCBpKam8u6777Lth8W0jqpJg1pVmHV352InqGxfjunBisd70qFRLS655BJmz57N1q1bGT16NHPmzKFNmzZcd911rFixotBuQ6VU6SpXSUpE3sAaBNHTGLMp3+HtWAmmd67yIUB3Tt+DSgCqYt13ynY5cEG+Mt2952brDaQAO/zyg6gScbvdzJo1i3bt2nHTTTdx/PhxPvjgAzZv3szo0aMJDg6mWoiTqJqhfh/s0LhxY9566y127drFc889R0JCAt27d6dLly58+umnuN35G/ZKqVJljCkXGzAROII1/Lxerq1qrjKPe8vcCLQF5mAll2q5yiwEfgE6YyWoX4DPcx2vgZXs5njruNFb5yPnirFDhw5GlZ5Tp06Z999/3zRr1swApk2bNmbWrFkmKysrYDEdP37cTJw40TRp0sQApnnz5mbSpEkmIyPDGGPMLe+uMre8uypg8SlVEQCJxtfc4OuJ/t6w7gcVtD2Xq4wAzwF7gEzge6BtvnpqAzO8ieeI933NfGXaAcu8dewB/gHIuWLUJOUf+b/YMzMzzbvvvmtiYmIMYC6++GLz8ccfG7fbHcAo88rKyjIfffSR8Q6eMREREebFF180vScsNF0nfGMSdxwIdIhKlVslSVJina+KomPHjkaXjy+57IdwPxhxEVOmTOGll14iOTmZyy67jGeeeYZrr72W8joBiDGG77//npdeeolvftpOxK0vYnM4CXE69FkrpQohIuuMMR19Obdc3ZNS54eDR46xYW0CsZ16MWbMGGJjY1m8eDEJCQn079+/3CYoABHhyiuv5Msvv+TxVyYjdieIjcxTLuZ+vyHQ4SlV6WiSUmVm48aNDH/oGTbtPUpm7aaEXPsYU+Yv5fvvv6dXr17lOjkV5PoubbDbbWAMxuPmP0/ex4MPPsiRI0cCHZpSlYYmKVWqTp06xUcffcRVV11F69at+XLtFsTuRGx2bHYnmdWjAx2izzo0qsXce7rwaN+WzPhrJ0YP7sPEiRNp3bo1n332WaDDU6pS0CSlSsWuXbt4+umnadiwIUOGDGHHjh2MHz+ejyaOs1ofgNNho3OTOgGOtGQ6NKrF/Vc1pXurKN544w0SEhKoU6cOgwYNYvDgwSQnJwc6RKUqNE1SqljONvu4x+Nh0aJFXH/99TRu3Jhx48Zx6aWX8sUXX7B161aeeOIJerVvQst61WhQK5SZd1W+gQaXXXYZiYmJTJgwgS+//JLWrVvz9ttv4/F4Ah2aUhWSju4rBh3dd3pk3v9Gn35eOj09nffff59JkyaRlJREeHg4d911F6NGjaJRo0aBCjXgtm3bxj333MOSJUu4/PLLmTx5Mm3atDn3iUpVMjq6T5U5YwyrVq1i+PDhREVF8fjjj9OgQQNmz57N7t27efHFF8/rBAUQGxvL119/TVxcHJs3b+biiy/mmWeeITMzM9ChKVVhaJI6j/hjocCDh4/wy/pE2l45kK5du/L5558zatQofv31V77//ntuvfVWgoKC/BRxxScijBgxgk2bNjF06FBeeOEFLrzwQr799ttAh6ZUhaBJShXKGMP27duZOXMm9913H22uGMCm1OMcr96I453+yrNvfkBycjJvvfWWdmOdQ1hYGB9++CGLFy/G4/HQs2dP/vKXv7B///5Ah6ZUuVYul+pQgXHq1CnWr1/PqlWrWLlyJatWrWLPnj0AVKtWjebXP8Bxh/Xwql3shLftRtWqVQMcdcXSq1cvfvnlF55//nlefvll4uPjaTDoUao0uYQ3hl5S6QaSKFVSmqTOI/mXXN+/fz8JCQk5CWnNmjU590tiYmLo2bMnXbp0oWvXrrRt25YNfx6xFgo0lWP4eKCEhoYybtw4br31Vu545J+kV4tFDhxnyDsrePfWNvRq3yTQISpVbujovmKoyKP71u08yE3vrMIYgw03juXvsCXhKwAcDgeXXHIJXbt2pUuXLnTp0qXQlWmvfWNZzkKB+lt/yf136WZe+WoziGDcWRxdNZteUR5Gjx7NVVddhc2mPfKq4ivJ6D5tSZ0nPl35Kx6PG7HZcXuEui06Mf66K+nSpQuXXnopoaGhRaqnWoiTaiFOTVB+cnlsXWy2LXgMhAQ5ubJzS+KnvcXcuXNp2rQpd999N3feeSfh4eHnrkypSkhbUsVQUVtSixYtYugDT1J90NOIPYiQIHulfJC2olq38yCrk/bTuUkdOjSqRWZmJh9//DHvvfcey5Ytw+l0csMNNzBq1ChtXakKqSQtKU1SxVARk9TkyZO59957adeuHbVv+icnHVW1q64C2bhxI5MnT+bDDz/kwIED2rpSFZI+zKvO4PF4ePLJJxk1ahS9e/dm2bJlhIXVLZUl11XpadWqFa+99hrJycnMmDGDyMjInAenhwwZwtKlS0ncsZ+J325l3c6DgQ5XKb/Te1KV0MmTJxk5ciSzZ89m1KhRTJw4EYdD/6orspCQEIYNG8awYcPYuHEj7733Hh9++CHzV/xMvaHjELsDp014ovMFdG0RSUREBDVr1qxwy58olZ929xVDRejuO3DgADfccAPLli1j/PjxPP744zlfVAXNu6cqrszMTB58J56v9wQhNjvGncWhFTM5snouAEFBQYSHhxMeHk5ERETOVtDnO2dv5OhJT56uYGMMbrcbl8uFy+UiKysr531Bn8fM/ZVMCWH8LR3o1T5W752pHHpPqoyU9ySVlJTEtddey/bt25k2bRpDhw4NdEiqlK3beTDn2bUgu/DIxQ6qntxHamoqaWlppKam5tnS0tJwuVx56giq35KIW19E7A5wZ3Hw03+Rsfs3srKyihxH7jqMO4v98/5BuO0Y0dHROVvDhg3zfK5du/YZLb38g0hU5aBJqoyU5yS1Zs0aBg4ciMvl4rPPPqN79+6BDkmVkeJ8sRtjOHToUJ6ktWBLJssP1wSxgfHQ1uzgQscenE4nDocDp9OZsxX2+dtUJ5/vEBBBMLRjJ9X/TGD37t3s2rWL5OTkM5JjaGhongQWHNWKhSdbgNh0BGolo0mqjJTXJDV//nxuu+026tWrx8KFC2nRokWgQ1IVyLqdBxk2ZTWuLA9Oh82n5HCuOjweD6mpqezevTsncWW/z96ON+xKje7Dra5Lj5tGR37hnh6N6du3L7Vr1/b3j63KkCapMhLIJFXY/aQ33niDhx9+mEsvvZTPP/9chyUrn/ijm62kdazZto8hk3/AGLDh5kT8BFJ/X43dbqdr164MGDCAAQMG0LJlSx0QUsFokioj5SlJud1uHnnkEd544w0GDRrEzJkzqVKlSkBiU8pfcie6i6NrsHbtWuLj44mPj2fDhg0ANGnShAEDBjBw4EB69OihS8NUAJqkykh5SVInTpxg+PDhfPrpp4wZM4ZXX30Vu90ekLiUKiu7d+/miy++ID4+nm+++YbMzEyqVatGnz59GDBgANdeey3h4eE6+KIc0iRVRspDknrrhliuu+461qxZw+uvv86YMWMCEo9SgXTixAmWLl2a08pKTk5GRGjf6wYOtL8dbHZCnDr4orzQCWbPA0czXezbl07n/k+z97ef+eSTTxg0aFCgw1IqIKpUqZJzj8oYw08//UR8fDxzfjqAQRCEzJMu/v3+x7w0shcxMTGBDln5SFtSxRColpS1zMZKPB4PeLJ4oWc4I67tVuZxKFXeZf9fMcYgHjcps/7OqZRNdO3aldtuu42bb76ZunXrBjrM847O3VfJrU7aj8djEJsduzOYI6H1Ah2SUuVSh0a1mHdvVx7t24p5D1zBHysXMm7cOA4dOsT9999PZGQk1157LTNnzuTYsWOBDlcVgSapCiC2qhvjdoHxEKQr4ip1Vh0a1eL+q5rSoVEtYmJi+Pvf/86vv/7Kzz//zKOPPspvv/3G8OHDCQ8PZ+jQoXz++eecOnUq0GGrQmiSqgB+Wjqf1DlPER6cpTeClfJRu3btGD9+PNu3b2fFihXceeedLF68mOuuu47IyEjuueceli1bhsfjYd3OgzqzfDmh96SKIRD3pIwxtG3blrSTDno++q5ODquUH7lcLhYvXsysWbOYP38+x48fp0H77th7P6IjBP1I70lVYj/++CO///47jS7rG+hQlKp0nE4n1157LTNmzCA1NZVZs2ZRv/2VGAQQMk66eOBfbzBx4kQ2bdqE/lJf9spNkhKRHiKyQESSRcSIyJ35jk/z7s+9rc5XJlhE3hKRdBE57q2vQb4yDUXkc+/xdBF5U0TK7SPrcXFxBAUFEd2hZ6BDUapSu+CCCxg6dChvP/cwNpsdjMEuhj0/LeOBBx6gVatWREdHc8cddzB9+nRSUlICHfJ5oTw9J1UV+BWI824FWQKMyPU5/93O/wDXA0OB/cBrQLyIdDDGuEXEDnzhPdYdqAN8CAjwoJ9+Dr9xuVzMnj2b6667DtsF1QMdjlLnBWuEYJfTs1b8+3qSkpL45ptvWLJkCV988QVxcdZXVKtWrejVqxe9evXiiiuuoEaNGgGOvvIpl/ekROQY8IAxZlqufdOAMGPMgELOqQHsA0YaY2Z690UDO4F+xpivRKQfVpJqZIzZ7S0zHJgChBtjjpwtrrK+JxUfH8/AgQNZsGABAwcOLLPrKqUK5/F4+Omnn3KS1rJly8jIyMBut3PppZfSq1cvrr76akIatGLd7qM6PROVcFqksySpQVitp0PA98BTxpg07/GewDdYyWZfrvN+A+YZY/4hIv8CBhtj2uQ6XhdIA3oaY749W1xlnaRuueUWvv32W1JSUnA6nWV2XaVU0Z08eZLVq1ezZMkSvvnmG9asWYM9olnOIpB2DP1CttClRSSxsbE0bdqUOnXqnFczuZf6tEgiEoLVivkz3/42xpjffLmwDxYBnwDbgRjgBWCptyvvJFAPcAPp+c5L9R7D+5qa73i697wCn5AVkVHAKICGDRuW+IcoqkOHDrFgwQJGjRqlCUqpciw4OJgrrriCK664gueff57Dhw/z1MxlxO8ExIbH42bm4jW8/fzcnHNq1KhB06ZNc5JW9hYbG0tkZGROAtPJcouQpETkBuAN4KCIOIC/GGN+8B6eDlxSivHlMMbMyfXxFxFZh9WV1x8reRVGgNzNxcKajgXuN8a8B7wHVkuqyAGX0Ny5czl58iS33357WV1SKeUHNWrUYGT/bizJXgQy2Mmi6W9RM+tfbN26NWfbtm0bP/74Ix9//DFutzvn/CpVqhAbG0u9tl3Y3KAfiA2nTXi4vZ1OsXUJCwsjLCyMCy644LxojRWlJfUs0MEYs09EOgIfisiLxphZWAkgIIwxKSLyJ9DMu2svYAfCsO5NZQsHluUq0zVfVWHe8/K3sAIqLi6OVq1a0aFDh0CHopQqpg6NajHzrs75WkERtGzZ8oyyLpeLXbt2sW3btjxJbOMxuzVZrtg4lZXF3/8Tx5HVp1tjwcHBOQmroK1uXSuhpVONHced9GwXzaWNw8rwT8E/ipKkgrLv8RhjEkWkB/CJiDSl8FZJqRORMCAK2OPdtQ5wAb2BWd4yDYBWwCpvmQTgaRFpkKvrsjdw0nt+ubBt2zZWrFjB+PHjz4vflJSqjDo0qlWkLjqn00lsbCyxsbH06dMnZ/+6nQcZ5m2NOYKcvPn0A9T2jCA9PZ309HT27duX8z49PZ0ff/yR9PR0Dh48PUtGUP2WOffGJn6/nazFr1LPkUFkZCT16tUjMjIyZ8v9OSQkJE+Mgex2LEqSShORC40xPwMYY/aLSG+sodsX+isQEakKNPV+tAENRaQ9cMC7PQd8jJWUYoDxWAMePvXGdVhEpgIvi0gap4eg/4w1dB3ga+A3IE5EHsEagv4yMPlcI/vK0owZMxARhg0bFuhQlFIBUnBr7NyysrI4cOAA6enpfLA6mbl/nAQRbDYbza64npCk79mzZw/r168nNTXVWl0hn5o1a+YkrQsateOXuj0DNgPHOUf3eVsjWcaYvQUc62qMWemXQESuBAoaXfchcC8wH7gYqImVqL4FnskeSu6tIwQr6dwGhGKN9rsvX5mGwNtATyADq9X1N+/gi7Mqi9F9xhiaNWtGo0aN+Oabb0r1Wkqpyi13a8zpsJ2RYNxuN/v27WPv3r3s2bOHPXv25Hm/Z88eUmq2w9PmWmsVBoGxfVpw/1VNz3LVM1W6IejlVVkkqVWrVtG1a1emTZvGHXfcUarXUkpVfiXtqjtXoiuKMlmZV0Ru9KH+hcaYDB/OO2/FxcVRpUoVbrzRlz9upZTKq6j3xs52vi/djv5SnGmR5hWzboM18i6pmOedtzIzM/nf//7HDTfcQLVq1QIdjlJKASVPdCVR3Ln76mXP8HAuInLUh3jOa/Hx8Rw6dEifjVJKKa/izIL+IdZAg6KaAZSbEXMVQVxcHJGRkVx99dWBDkUppcqFIicpY8xIY0yRW0fGmHuNMfmnKFKF2LdvHwsXLmT48OHY7fZAh6OUUuVCuVlP6nw3Z84csrKytKtPKaVy8Wk9KRG5AutB2brAVmBDru13Y0yW3yI8T8TFxdG+fXvatm0b6FCUUqrc8HXRw6lYCxS+jjVLxMXATUADrKU0Qgo/VeW3ceNGEhMTee211wIdilJKlSu+Jql6QB9jTJ7h5SJSGythqWKYPn06drudoUOHBjoUpZQqV3y9J/U90Dj/TmPMAWOMzuVTDB6Ph+nTp3PNNddQr16BS1oppdR5qzgzTnwF/OTdPgSeEZGfc6+Cq4rvu+++488//+Tll18OdChKnXeOHDlCWloaLpcr0KFUWE6nk/DwcKpXr14q9Renu2890B4YAUR49/0hIp9jLYGxHvhZp0Eqnri4OKpXr871118f6FCUOq8cOXKE1NRUoqKiCA0N1WVxfGCMISMjg+TkZIBSSVRFTlLGmCey34tIPayElb09hDWAwojIFmNMa38HWhkdP36cefPmceuttxIaGhrocJQ6r6SlpREVFUWVKlUCHUqFJSJUqVKFqKgoUlJSApukvAG9irVkxkpjzCJgUa5joVgJy29rTFV28+fP5/jx4/pslFIB4HK59JdDPwkNDS21LtPiju6rAswGgkXkC6wFB782xmR4u/kSvJsqgri4OGJiYujWrVugQ1HqvKRdfP5Rmn+OxRrd553qqAHQH0gGXgTSRWSBiPxFROqWRpCVUXJyMkuWLGH48OHYbDrxh1JKFcSnb0djzBpjzFPGmLbARVhD0u8E/hSRFSLyNxGJ8mOclc6sWbPweDyMGDEi0KEopUpoyKQEhkzSTqTSUOJf4Y0xW40xrxpjemDNOPE+0A3QJ1MLYYwhLi6Ozp0707x580CHo5SqgNLT0xkxYgSjR48OdCilyq/9TMaYfcaY940xg4wxr/iz7srkp59+4tdff9UBE0opnz3zzDNER0czc+ZMv9f99ttv07hxY0JCQujQoQPLly/3+zWKyqckJSIvisg9Bey/R0SeL3lYlVtcXBxOp5MhQ4YEOhSlVAWUlZXF3Llz6devH/Xr1/dr3f/73/8YM2YMTz75JOvXr6dLly7069ePXbt2+fU6ReVrS2oE1sO7+a0DtHlwFllZWcycOZMBAwZQu3btQIejlKqAEhISaN68OevWraNr165+rfu1117jzjvv5O6776ZVq1a89dZbREZG8s477/j1OkXla5IKBwqaDmk/p2ejUAX4+uuvSUtL064+pSqRo5kukg9lsG7nwTK53vLly+ncuTMLFizgxhtvPOP4uHHjqFq16lm3grrwTp06xbp16+jTp0+e/X369GHVqlWl9vOcja+zoO8CugNJ+fb3AP4sUUSV3PTp06lduzbXXnttoENRSvnBup0H2bT3KB4Dw6asZuZdnenQqFapXjMxMZFmzZqxZcsW+vXrd8bxe+65h1tuueWsdURFnTkAOz09HbfbTURE3rZGREQES5YsKVnQPvI1SU0CXheRIGCpd9/VwHjg3/4IrDLJHpr63q2tmT9/Pn/9618JCgoKcFRKKX9YnbQfj7Heu7I8rE7aX+pJaseOHWzcuJGxY8ficJz5NV67du0S3U7I/3CuMSZgDz77+pzUq1iJ6k1gs3d7A5hsjHnJf+FVLvPmzSMzM1OfjVKqEuncpA427/e302Gjc5M6pX7N1NRUjh07Vujwc1+7+8LCwrDb7ezduzfP/rS0tDNaV2XF15YUxpi/i8gLQGtAsJaNP+a3yCqhuLg4mjdvTqdOnQIdilLKTzo0qkXLetU4kpnFG7deXOqtKACHw8GLL75Y6OS4vnb3BQUF0aFDBxYvXszNN9+cs3/x4sUMHjy4ZEH7yOckBWCMOS4iG73vNUGdxfH0PSxbtowXXnhB5wtTqpKpFuKkWoizTBJUfHw8u3btIjg4mN9//53169czbNiwPGVK0t03duxYRowYQadOnejatSvvvvsuKSkp3HPPGU8dlQmfk5SIPASMBaK8n1OA14D/GGOMf8KrPHausSaMHz58eIAjUUpVVG63m2nTpjFjxgyeffZZIiMjmTp1ql+vMWTIEPbv388LL7zAnj17aNu2LV9++SWNGjXy63WKyqckJSIvAaOAlzk96/nlwLNAJPCYX6KrJI5kuNh73NB5wNCA/UUrpSo+u93OvHks0gS8AAAgAElEQVTzAM5oPfnTfffdx3333Vdq9ReHry2pu4C7jDHzcu1bKiJ/YA2o0CTltW7nQTbtOUyVSwayz25j3c6DZdIloJRSlUFJ7kn9XMg+XXcil9VJ+zGA2OwY72dNUkpVLv8bfXmgQ6i0fE0occD9Bey/F5jueziVT+cmdbzrRZkyG56qlFKVha9JKhi4U0Q2icg077YR+AvgEJE3s7eiVigiPbyLJyaLiBGRO/MdFxF5TkRSRCRDRL4TkTb5ytQSkekicti7TReRmvnKtBOR7711JIvIs1KKw+2yh6c2qFWlTJ5EV0qpysTX7r6WwI/e99kjAfZ6t1a5yhVnlF9V4FesVlpcAccfAx7BWlzxD6xBGotFpIUx5qi3zCygIdDPe+0pWC27gQAiUh1YDCwDLgVaANOA48CrxYi1WMpyeKpSSlUmPiUpY8xV/g7EGPMl8CWAiEzLfczb0nkImGCM+di77w4gDbgNmCQirYC+QDdjzCpvmdHAcm8i+wMYBlQB7jDGZAC/es8bKyKv6dB5pZQqXyrKIIfGQD3g6+wd3iSzDOji3XU5cAzIPVXvSqxWUu4yy73nZvsKqA/ElEbgSimlfFeslpSILChKOWPMdb6FU6h63tfUfPtT8T5M7C2zL3dryBhjRCQt1/n1OHOW9tRcx7bnv7CIjMJ6JoyGDRv6Gr9SSikfFLe7bwCwE/jO/6EUSf7uOMm3r6DuunOVkUL2WzuNeQ94D6Bjx47aHaiUUmWouEnqFWA41rpRHwDTjDFlsX5U9pS89YDdufaHc7oltBcIFxHJbk1572XVzVemHnmFe1/zt9L8Rp+hUKqS+6C/9Tryi8DGUQkV656UMeYxIBp4GOgIbBGRhSJyk4g4SyNAr+1YCaZ39g4RCcFaeDH7HlQC1gjB3BnhcuCCfGW6e8/N1htIAXaURuBKKVUa0tPTGTFiRKHLdVQWxR44YYxxG2MWGGMGYQ1o+BZ4AUgWkaq+BiIiVUWkvYi098bV0Pu5obdl9B/gCRG5UUTaYg0dP4Y17BxjzEZgEdZIv84icjnWFE3x3pF9eMueAKaJSFsRuRF4AtCRfUqpCuWZZ54hOjqamTNn+rXe8ePHc+mll1K9enXq1q3LwIED+fXXX/16jeIo6ei+C4CaWC2YYxTvuaj8OgLrvVso8E/v+395j7+ENcv6RCARayLbPrmekQJriPlPWKMAv/K+z1lh0BhzGKvlVN9bx0Ss56NeK0HcSilVprKyspg7dy79+vWjfv36fq37u+++47777mPVqlUsXboUh8NBr169OHDggF+vU1TFfk5KREKBW4C/YiWWT7GeO/qmJIEYY77j9CCGgo4b4DnvVliZA1j3zM52nV+w7qkppVSFlJCQQPPmzVm3bh1du3b1a91fffVVns/Tp0+nRo0arFy5koEDB/r1WkVRrJaUiLyHdW/oQWA2UN8YM6ykCUoppSq0k0fg8G7YvaZMLrd8+XI6d+7MggULuPHGG8847uvy8QU5evQoHo+HWrUCM2NOcVtSdwG7gD1YUw/1K2jau1J4Tkr5g45AUsr/dq+B1F/BeODD6+COBRDdqVQvmZiYSLNmzdiyZQv9+vU747ivy8cXZMyYMbRv357LLw/MKOXiJqk4SnbfSflKE4xS5dOO5VaCAnCfsj6XcpLasWMHGzduZOzYsTgcZ36Nl2T5+NzGjh3LihUrWLFiBXa7vcT1+aJYScoYc2cpxaGUUhVTTHcQm5Wo7EHW51KWmpqKzWYrdPj5uHHjGDdu3FnrWLhwId27Fx7rww8/zJw5c/j2229p0qRJieItiSInKRHpBKwzxriLWL4D8LMxxuVrcEopVe5Fd4KItpB5GAZPKfVWFIDD4eDFF1+kSpUqBR4vaXffmDFjmDNnDt999x0tW7YsUawlVZyWVALe+fGKWP5boD2QVNygVDmm3Y5KnSm4urWVQYKKj49n165dBAcH8/vvv7N+/XqGDRuWp0xJuvvuv/9+pk+fzvz586lVqxZ791oT/mQPuChrxUlSAowXkRNFLB/kQzzqfKCJTimfuN1upk2bxowZM3j22WeJjIxk6tSpfr3G22+/DcDVV1+dZ/8//vEPnnvuOb9eqyiKk6SWAbHFKJ8AZJyzlFK+0ESnzkN2u5158+YBnNF68pfyNvlOkZOUMebKUoxDKaWUOoOvy8crpZTKpi36UlNRVuZVyv8+6H+621ApVS5pklJKKVVuaZI6n5Tx/GJKKVVSmqTOE2bXD3j2/oI5uBPPtAHsTVzA4cOH8Xg8xavIH4muvNShlCr3fB44ISLBWOsyhQL7jDFFfchXlQGPx8Pvv//OsmXLWL58ORceXsKjHT04bELWqUzeePhmJqw4BUC1atWoUaMG1atXp0aNGjlb7s/Vq1enafAB+uz5BRsG835/VsQ+xr6Qxrjd7nNuWVlZuN1u6rl2MfTkz9gFsqb0Y8qpgWw7VYesrCxcLhdZWVl53hf02qLKYV5vvxWnDdxT+vJZjb9yQcuraNKkCTExMYSEhJzjT0cpVVEUK0mJSDWs9ZqGAp0AJ9ZDvkZEUrBWxn3PGLPW34Gqs3O5XKxfvz4nKa1YsSJnkbKoqCia9+mIIQGPMYgjiC5Dx/DqDbU4fPgwR44c4fDhwznv9+/fT1JSUs7+jAzrcbcnugXR56pgbDbBlXWSRZOezkl0RfVEtyDkqmBsIojbxZ6EuUxKFJxOJw6HI+c19/v8r5fUOYHTBnYbeNxZ/Pjpm0xY8QoAIkJUVBRNmjShSZMmxMbG5nkfFhZGzsz9J49YU9nsXlMmMwUopYqvOHP3PQw8jTXN0QLgRSAF64Hd2kBboDuwWERWAw8aY7b4PWIFQEZGBj/88APLly9n2bJlJCQkcPz4cQCaNWvGoEGD6NGjBz169CAmJsb6Yn63O2QexjZ4CgOL8aV86tQpjh49SubWZdgWjsQYg90Zwsh/vs1tdS/Cbrefc3M4HNbr3h+xz7gOYzw4g6vw/PQFPF/cBLF7Dbx/DRgPjuAqPPrfj7juRA2SkpLYtm1bzuvXX39NSkpKnlOrVq1KkyZN6NcujH81/hmHDeTDgcgdn/uWqPShYqVKVXFaUl2AK4wxhS12vwZ4X0TuwVq19wpAk5SfuI4fJunPPSx9ciQzvt/M2rVrcblciAgXXnghI0eOpEePHnTr1o3IyMiCK/FxfrGgoCDq1KkDdW6A9a9B5mFk8BSa+/KlHtu95JNx5prQUwZPoXZ0Jy6HAte7ycjIYPv27XkSWFJSEhEnNmATg00E18kMPnjydtJb3U7fvn1p3749NpverlVFN3LRSAA+6PtBgCOpfIoz48TN2e9FpKExZlch5U4Cb/shNuV1atsKPHt/ITbIEO3+mB+rNabbww/To0cPunTpUrYrZvpjIs0yrCM0NJTWrVvTunXrvAd2r8G8fw0ejwdjc7Jqj40PZzzFU089RXh4ONdccw19+/ald+/e1K1b1/c4lSol6enpPPzww1SpUoVJkyYFOpxS4+uvi594B06cQUT0rnV+JXho1BjDp68/gsNmcNiE0CAHk5+6nX//+9/0798/YEs6V3jRnZCItthqNSLoroVMW/I7e/fuJS4ujl69erFw4UKGDRtGREQEnTp14plnnmHlypVkZWUFOnKlAHjmmWeIjo5m5syZpXaNcePGISI88MADpXaNc/E1SW0F3su/U0TqA8tLFJHK46mnnuI/n63DYN3slzJaVO28EFwdakTntMYiIiIYMWIEM2fOJDU1lbVr1/Kvf/2LoKAgxo0bR7du3QgLC+Pmm29m6tSp/PnnnzoUXgVEVlYWc+fOpV+/ftSvX79UrrF69WomT57MhRdeWCr1F5WvSeovQAcReTB7h4i0x7ovtc0fgSmYNGkS48eP56L+d+GIbAc1G8EdC3QkWhmw2Wx07NiRp59+mhUrVpCens7cuXO56aabWL16NXfddRc3Xx5Dxp8/4T64E/cH/Tm1bUWgw1bniYSEBJo3b866devo2rWr3+s/fPgww4YNY+rUqQHvrfEpSRljTgCDgX+ISDcRGYTVgnrfGHOrPwM8X8XHx3PffffRv39//vvf/yIhNfL81q/KVq1atbjpppuYMmUKu3bt4tdff+WlewcSZBPsAh7XSV68uy8DBw5k4sSJJCUVY61PnUOwwjt26hh7ju9hQ9qGMrne8uXL6dy5MwsWLODGG2884/i4ceNyFiksbFu+vPBOr1GjRnHTTTfRs2fP0vwxiqQ4Q9C/AjYA672vfwCjgHhvPXcbY2aXRpDnm7Vr1zJkyBAuvvhi5syZg8Ohk9WfIYBDvkWENm3aQPW/w/vfYYwHmyOY6hf25PfPE4mPjwesRwH69u1L3759ufLKKwtd6ltVbBvSNvDHwT8wGO7++m4m95lM+/D2pXrNxMREmjVrxpYtW+jXr98Zx0uyfPzkyZPZunUr06dP90usJVWcb7/1WMvB3w5EACeAXwA38BGwWUSCvaP7lI+2b9/OgAEDCA8PJz4+PiDLNZe68vJMUUnj8A6Fl8zD2AdP4ZHoTjzyH9i6dSuLFi1i0aJFTJkyhbfeeovg4GB69OhBv3796Nu3Ly1btjz9ULGq0BJTEzFYCwW6PC4SUxNLPUnt2LGDjRs3Mnbs2AJ/ifV1+fg//viDJ598kuXLlxMUVD4WVy9yd58x5gljTF9jTCQQCdwEfAYsBnoAPwBHReS3Uon0PLB//3769euHy+Vi4cKF1KtXL9AhqXPJN/gCoGnTpjzwwAPEx8dz4MABvv76a+6//36Sk5MZO3YsrVu3JiYmhtGjRzN//nyOHzmggy8qsI4RHRHvwCanzUnHiI6lfs3U1FSOHTvG6NGjCzzua3dfQkIC6enptG3bNmfml++//563334bh8PByZNl3wbxqR/JGJOKNQXSoux9IhKK1dIK7FCQCiozM5Prr7+eHTt2sGTJElq2bBnokApWXlpBFURISAi9e/emd+/evPrqq+zcuZOvvvqKRYsWMXv2bH7+cip9bq9Clh08k/vyRdhoorsMpl27dgQHF/iUx5l01ouAah/enha1WnDUdZQJ3SeUeisKwOFw8OKLLxbahexrd9+gQYPo2DFvkh05ciTNmjXjySefDEjrym83O4wxGUCCd1PF4PF4uP3221m5ciUfffQR3bp1C3RIqpQ0atSIUaNGMWrUKFwuF3/OHkvwtjjsAi6PmzVzX+PGMRNwOp1ceOGFdOjQgY4dO9KxY0fatm2L0+kM9I+gClA1qCpVg6qWSYKKj49n165dBAcH8/vvv7N+/XqGDRuWp4yv3X01a9akZs2aefZdcMEF1K5dm7Zt25Yobl8VZ+BEY2PM9iKWFaCBMWa3z5GdRx599FHmzp3LK6+8ws0333zuE3ylv2n7Xwn+TJ1OJ42vGgHbZ1jzEAaF8n+vTafjLheJiYkkJiby0Ucf8d571iOJwcHBXHTRRTlJq2PHjrRq1QqHTpR73nC73UybNo0ZM2bw7LPPEhkZydSpUwMdVqkqTksqQUS+AKYYYwpsLYlILeBW4P+AicB/Sx5i5fbmm2/y2muv8eCDDzJ27NhAh6PKWr55CCOjOzH4Uhg8eDBgzTiSlJSUk7QSExOZPn06b79tzTx2RWwVFt3qwGkHz9Rr+an9v4i89Hrq169fvIEZ2mVYIdjtdubNmwdwRuuptHz33Xdlcp3CFCdJtQSeAr4QETewDtgDZAK1gNZAK6wHeh8yxnzl51grnU8//ZSHHnqIG264gddff11He52vzjIPoYgQGxtLbGwsQ4YMAazu4a1bt5KYmEiNn6fisCVaz2q5T/Hxqw8zYcX9VK1alRYtWtCyZUtatGiR875Zs2aEhoaW9U+olM+KM8HsIeBREXkW6A90AxphLXqYDnwIfHWWWdJVLgkJCdx2221cdtllzJw5E7vdHuiQVAVhs9lo3rw5zZs3h+5N4f1rMMaDPSiUG8dOIHoobNq0iT/++IPly5fnmdtNRGjUqFGeBNayZUva7ksnzJmBaJehT3T289JT7IET3gES87xbmRKR54B/5Nudaoyp5z0u3uOjsFp3PwD3G2N+y1VHLeBN4DrvrgVYa18dKt3oT9uyZQsDBw6kQYMGLFiwQH+zPd+VpIst17NaMngKl0Z34tJ8RU6cOMGWLVtyElf264oVKzh+/DidG9j55vYquO2Q9V4fJuzpgiOmS57WV5FWO9YuQ1UKSjy6T0TaAJcDvxhjfsi1vyrgKoWHe/8Arsz12Z3r/WPAI8Cd3nLPYi3C2MIYc9RbZhbQEOgHGGAKMB0Y6Oc4C5SWlka/fv0QERYuXKjLQKiSO8eyJVWqVOGiiy7ioosuyrPfGENycjKZS8YRsn0mNgGMIXTvOp6YujCnnIgQExNzRuurRYsW1KtXT1c6VqXKH0PQR2F1+UWKyOPAaqwBE8HAu1jLzftTljFmb/6d3lbUQ8AEY8zH3n13AGnAbcAkEWkF9AW6GWNWecuMBpZ7E9kffo41jxMnTjBw4EBSUlJYunQpTZs2Lc3LKXVWIkKDBg3g6jvh/dnWCENnKI+/s4AHardhy5YteVpemzZtYtmyZZw4cSKnjurVq9OiRQsGXBTO3yN/xi7AtAG4h32Ks4n/Jz71N2OM3gv2A2NMqdXtjyT1JHAvcBVWq+Z64EGsqZJK48mvJiKSDJzC6s570hiTBDQG6gFfZxc0xmSIyDKsVYUnYbX4jgGrctW3EjjuLVNqScrtMdx2222sXbuWTz75hM6dOxevAu1CUYXx0/ROuVdLvgBo37497dvnfe7H4/GQnJyck7SyE1hQyhok0rvS8alMnruzF/NSG9KqVaucrWXLlrRs2ZLq1asXHksZdhk6nU4yMjJ0TkU/yMjIKLVn+EqcpIwxx4FXgFdEJAhogZUs6gH+Tq8/YHXlbQLCgaeBVd4ux+w5hFLznZMKZD9aXQ/YZ3KlfWOMEZG0XOfnISKjsFqLNGzY0KegTeZhvvhhC6nr9vPmm28yaNAgn+pRqtQUcaVjm81GdHQ00dHR9OrV6/SB3WswU6/B4EEcQTTscTNtNx5i06ZNfPnll7hcrpyiUVFROUkrdxKLiIhAyrDLMDw8nOTkZKKioggNDdUWlQ+MMWRkZJCcnExERESpXEOK20wTkYuAi7G60ZYYY06dpazNGOMpWYhnjaUqkARMwOpmXAk0zP0QsYh8AEQaY/qKyJPAX40xsfnq2Q5MMsZMONv1OnbsaBITE4sX5O41uCb3QfDgEQdBdy3S/npVOb3bPU9rLJvL5SIpKYmNGzeyadMmNm7cmLMdO3Ysp1zvljVYMNjgtIHH5mBJ1P9RvXUvGjduTL169bDZijjVaDFaY0eOHCEtLS1PElXF43Q6CQ8PP2sLWUTWGWN8mtSwWC0pb6viHSD7V44tItLTGJNcUPnSTFDe+o95J7RtBsz37q4H5J7pIpzTrau9QLiISHZrynsvqy5ntsD8Y8dyHDaDIBixPmuSUpVSIa0xp9OZ86xWbtkDN7ITV8Odc3HafsJuA4/bxbK4cUxY8RxgzYEYExNDkyZNaNy48Rmveb4gi9Eaq169+tm7H1XAFbe77zHgbeB5IBr4D/Bv/D84okhEJATrIeNvge1YSag3sDbX8e7Ao95TEoCqWPemsu9LXQ5cQN77VP4T0x0RGxiPLv2uKrdi3kfKHrjRoEEDq+tw92Xw/jXWAI7gUO6ZMJUeR6qyfft2tm/fTlJSEtu3b2flypUcPnw4T121a9emSZMm9Gldi380/BmHDTzv92fz5a9Qs21vwsPDi7cumw6nLzeKm6QaAa8YY9KANBG5E2tNqTIhIq8AnwO7sFpIz2AlmA+995b+AzwlIpuAzVj3rI5hDTvHGLNRRBZhjfS7G6tFOAmIL7WRfQXclFZKFSDfFFGNojvRqJCiBw8ezElauRNY7cO/YBNrAIc76yTTnx/FhBWnsNlsREREUL9+faKioqhfv36eLXtfnTp1rHtTOpy+3ChukrIDGdkfjDHbRAQRiTTG7PFvaAVqAMwGwoB9WPehOhtjdnqPv4Q1HH4ipx/m7ZPrGSmAYVgP82aPAlwAPFCqURfxprRS570i/l+pVasWHTp0oEOHDnkP7F6DyZ6BwxlC//ueI2b4BSQnJ5OSkkJKSgo7d+5k1apVpKenn1FvUFAQ/S+sy6xrjuK0gXtKX+ZWuR1bw8uIioqiQYMG1K9fXx9uLkO+jO4bJSKrgA3GmANYw87LZMoEY8yt5zhugOe8W2FlDhCg7kml1Dn4YTi95GqNdYvuRGEL35w8eZK9e/eSkpKSJ4ldcuI7nLaN3ntjWfwa/w4TVryR59w6derQoEEDoqKicpJX7teoqChqZh62Ritqa6xEijW6T0SWYi1sWBNreHkK1kCFV4FvgERjzMFSiLNc8Gl0H+hvVEpVJLvX5NwbwxHK8Ztms9MTQXJyMsnJyfz55595XpOTk0lLS8tTRfZUU0F2cBs7bx/vi61RZ2JjY3MGe5xP06GV2eg+Y0xP7wWbAB1ybXdhDaowIpJkjGnmSzBKKRVw+e4jXxDdidZA69atCz3l5MmT7NmzJyd51U+aQ3DmUmsGDo+Hwz/F88/X5+Q5p379+jlJK/9r3bp1rXtjhQzrP58U+zmpQisSiQE6ApcYY570S6XljLaklFJFkq81Zm7/jPTQJiQlJbFt2za2bduW8z4pKYnk5LxP8VStWpVBHeszpdteHDYwNiebOr1EVOcbqFWrVoB+KN+VWUvqbIwxO4AdBGB2dKWUKleiO8FfvrKei4zpjkR3oi5Qt25dLrvssjOKZ2RksGPHjjzJ65IT32EXg10EV9YpZo67lwkr/kpERETOFFO5t4YNGxb8wHMF/yXZb0lKKaVULtGditxFFxoamjM9VA5va8xaK8waqVj7OnfOnIkfffQRBw+eHgIQEhKSM0N9nu3oIUI8RyvsAA5NUkopVR7lWyss/0hFYwzp6ek5SSt7W7t2LR999BHGmJwBHFnetcLeOtwbZ5OuNGvWjObNmxMTE1NqE8P6i09JSkQaGmN2+TsYpZRSudyzvNBDIkLdunWpW7cu3bvnnckmMzOTLVu2wPLXCEn9xForzGPI2LiYx946fUfG4XDQpEmTnKSVvTVr1oyoqKjT3YcBHMDha0vqExHpWtCChiISYozJLGFcSimlfBQSEkK7du2g5r3w/nxrqqmgUJ6Z+hn3VYll8+bNOduWLVvYvHkzS5cuJSMjZ64GQkNDadasGQMuCucfjX7GaQf58Dq4Y0GZJipfk9RW4D3gjtw7RaQ+8BmcsYK1UkqpslbAAI4wICwsjC5duuQpmr1WWHbSyk5gtY/+iE2sSbJxnyrzSbJ9TVJ/AdaIyIPGmLcARKQ9EA+s8FdwlcVIsSZY/yDAcSilzkNFHMCRe62wnj17nj6Qa6qpQEyS7VOSMsacEJHBwEoRWY81l9504HVjzLP+DFBZRi4aCcAHfTXVKaXKUHQnJFdrrNzekxKRr4ANwHrv6x9YK9bGe+u52xgzuzSCVEopFUDFGE7vb8VpSa3HmrfvdiACOIG1TIcb+AjYLCLBBQ2mUEoppXxR5CRljHki+72IRGAtId/eu/XAuk/lEZEtxpg2/g60QqvXLtARKKVUheTrPalUYJF3A0BEQrES1oX+CU0ppdT5zp9z92VgLc+e4K86lX/p4AulVEVTwGyESimlVPmgSUoppVS5pUlKKaVUuaVJShXLyEUjc+5tKaVUadMkpZRSqtzSJKXKnLbGlFJFpUlKVUia6JQ6P2iSUkopVW5pklLnLW2NKVX+aZJSqgQ00SlVujRJKRVgmuiUKpwmKaUqOE1yqjLTJKWU8kui0zpUadAkpZRSpaC8JMvyUoev/LZUh1JKVRbHTh3jqOsoG9I20D68fYWv48ipI6xPXU/b2m3xeDwFbm63u9D9qSmpuEJdJYrDV5qkKgCPx8PBIwc5euoo3/7xLW1q+bbw8f4D+zlhTrB823IuqnsRDocDu92es4nIOesoT//xtI7yE0Np1WGMwePxkJWVVeTtz6Q/ybRlErckjoaOhnmOuVyuM8rn37dH9rCpxiYQuPOLO+mS3IUax2vgdrsL3LK/yHNvR6sdZcdlO8AGt39+OxHLInDuc541KeT/nBWehfsmN2IThi8YzqkZp3DvznvN/K/59wXFBNHwkYaIQxj+2XC2v7SdjG0Zxfo7CY0NpfFjjbGdtHH313czuc/kMk1UYowps4uVJyJyH/AoEAn8BjxkjFl+tnM6duxoEhMTi32t62dfz8GjB/lLzF8IOxXGkSNHcrajR4/m+VzQ/qzwLBo/1hhxCCbLlOgf2tnqsNlseRJX/vdBMUFUubMKYhfwgGeOB3uqddxms+VJeAXts9lsuOq6SOmWAnawGRuN1zam2tFqhZ5T0L5DVQ6xLHIZCNixc/WBq4lwR+SUs9lsed4XtC+FFKa7poOAQxwMcw4jSqIAyP4/kfu1oH0ppDDLM8uqAwcDMwdS11WXrKysnC+M7Pf5X7PfHww9yIYWG8Bm/Xk0/bEpIQdCCv1CzL+dCjvFiQEnwA7iFuyf2LHttZ0z9jyvkQZuIefv1fzPQErx/o1Tnzx1ZM3Kwr3b+qLMTjTneu9o6CDsvjDELhi3Ydcruzi2+VixwijKv/NzCesfRsSNETlxpH2axqGvDp3x7/Js/2YdXR3Yutty/jzsq+04f3Tm+beY/99m/s+HWx5mf6v9iM2qI2JrBJG7Igv8t13Qq91uZ1udbWyss9G6sWOg3eF2tDvRLs81zxaPzWZjrW0tK1hh/X8TOw9c/AB3tburWH+mIrLOGNOxWCdln3s+JikRGQLMAO4DVnhfRwKtjTG7CjvPlyS1IW0DI+JHYDAF/qex2WxUr16d6tWrU61atZz3uT/vabiH32r9BgJihE6uTlziuqRYcTJpj8UAAAwtSURBVPzo/JEfnD/k1HHxiYtpc6zNWb9Q87/fXX83f8b8mfOfpu4fdam5uWaBv0nm35f92XWJC/fl1m+Hxm1wL3NzavmpQn8jLaieWn1rEX5DeM6XSOonqaR/kV6sP4/8X0RlUUfuxJ/9Wr13dapeUzWnDtd3LjwJnkK/EPNvGRdmcOTCIzl/JzV/q0mtzbUAclrGIpLnff7X/U33s6/Fvjx/r3W21SnWn8X+2P2ktUjL84UasSMCEcn5sjvX+12Ru9gWuc2qw0Cr9Fa0OtoKh8NR5G2VexVLMpdY/84R+tfoT/+w/meUczqdhe7beHgjdy+9G4MhxB7iU8thQ9oGbl94u9bhVZIkdb52940FphljJns/PygifYF74f/bu/tgueY7juPvT3JzE6RmSgSt0UpFqh6SVIjoqOlUpsoUJQRVzRjUoBRTGqNlph2PmWB0VIkwgnqsIhg0nVRLiHioqspQhCIhnSCekjDf/nHOTTeb3Zvds2ezZ3c/r5kzuXvO73zP73xz7v3e39mz98fUPA+0YMkCGJD8MBg4cCBnXH4GR3zliNUFaMMNN1znbbbSi2Rwz2BO2vekTBfa/Pvnr45x6qRTG7tYBw1h+mnTG4qxQe8GXH1Rg980vUO498p72XGTHde61dHf18+/+zxnP3c2EcHgQYOZcc4MRk0b1e8P8/J1C5cv5Mynz1ydj1un3cpO1+20ViHq6elhwIDKzyiV5+OGaTfUlY/y/5PLz7i85f+vQwYNYdop0xqL0TOEs446q+4YI98eyZz75yTX+cDBTN5jct0xxm88nlGfH8XyVcu5YM8LMt3aGjN8jGPkpOuKlKReYBdgWtmmB4E98j7euM3HIYkg6B3Yy4G7HMh2w7erK0ZRLrROirEt2zLrzVkNxRjBCGYumtnScylKPjspBsDQ3qEM7R3a0A/kPGLctv9tmfctWoysuq5IAcOAgcCSsvVLgL3LG0s6DjgOYOutt677YJ32TeMYxYtRhD50Woxr97k28755xrDuLFJ9yt+MU4V1RMRVwFWQvCeV5UB5fNOYdQMXByvXjUVqKfAZsEXZ+uGsPboy6wouDlZUXVekImKlpCeBiUDpjdaJwB2t6ZVZdi4O1sm6rkilpgOzJM0HHgGOJ/mkx5Ut7ZV1JRcZs+q6skhFxC2SNgXOJvkw73PAvhGxqLU9s3bjAmPWXF1ZpAAi4grgilb3w8zMquvaImXmUZBZ8blIWVtygTHrDi5Stt65wJhZrTzpoZmZFZZHUlYXj4LMbH3ySMrMzArLRcrMzArLt/u6iG/VmVm78UjKzMwKy0XKzMwKy0XKzMwKy0XKzMwKyw9OtAk/9GBm3cgjKTMzKywXKTMzKyzf7lsPfKvOzCwbj6TMzKywXKTMzKywXKTMzKywXKTMzKywXKTMzKywXKTMzKywXKTMzKywXKTMzKywXKTMzKywFBGt7kPbkPQOsKjV/WgDw4Clre5Eh3FO8+ec5qu/fH4pIjbLEtRFynInaUFEjGt1PzqJc5o/5zRfzcqnb/eZmVlhuUiZmVlhuUhZM1zV6g50IOc0f85pvpqST78nZWZmheWRlJmZFZaLlJmZFZaLlJmZFZaLlK2TpBMkvSLpE0lPStpzHe33Stt9IullSceXbT9XUpQti5t7FsVRTz4lbSnpJkkvSPpM0nVV2h0s6XlJK9J/v9+0EyigvHMqaUqFazQkDWnqiRRInTk9SNKDkt6RtFzS45L2r9Cu7uvURcr6JWkycBlwHjAWeBS4X9LWVdpvA9yXthsLnA9cLungsqYLgS1Llp2acgIFU28+gcEkn+K/AHi8SswJwC3AjcCY9N/bJI3Pt/fF1Iycpj5izWt0y4j4JK9+F1mGnO4F/BnYL21/H3BnaWHLep366T7rl6THgWcj4tiSdS8Ct0fE1ArtLwQOioiRJetmADtExIT09bnApIjYsdn9L5p681m272xgaURMKVt/C7BJREwsWfcn4J2IODzP/hdRk3I6BfhNRAzNv8fF10hOS9rPB/4aEaenrzNdpx5JWVWSeoFdgAfLNj0I7FFltwkV2j8AjJM0qGTdCElvpLcTbpY0IpdOF1jGfNaiWs4bidkWmphTgA0kLZL0H0mzJY1tMF5byDGnnwOWlbzOdJ26SFl/hgEDgSVl65cAW1TZZ4sq7XvSeJDcYpkCfBc4Nt3nUUmbNt7lQsuSz1pUy3kjMdtFs3K6EDgaOAA4HPgEeETSyH736gwN51TSicBWwKyS1Zmu055aDmhdr/yesCqsW1f71esj4v41NkqPAS8DPwKmZ+9m26g3n62K2U5yPf+ImAfMWx1MehR4BvgJcHLWuG0mU07T958vBg6LiPJZI+qO6ZGU9Wcp8Blr/6YznLV/I+qzuEr7T4H/VtohIj4A/gl0+m+pWfJZi2o5byRmu2hWTtcQEZ8BC+j8axQayGlaoGYBR0XE3WWbM12nLlJWVUSsBJ4EJpZtmkjytE8l84C9K7RfEBGrKu2QPtb7VeCt7L0tvoz5rMW8JsRsC03M6RokCdiZDr9GIXtOJR0K3ABMiYjbKzTJdp1GhBcvVRdgMrASOAbYnuSx1A9IJjEDuB64vqT9NsCHwKVp+2PS/Q8uaTON5JHVbYDxwGzg/b6YnbzUm8903Zh0eRi4O/36ayXb9yAZqU4lKfZTgVXA+Fafbxvn9BzgO8CIdNvMNKe7tfp8i5hT4LA0P6eQjJb6lk0avU5bngwvxV+AE4BXgRUkv2F9s2TbXGBuWfu9gKfS9q8Ax5dtvxl4M/0meAO4o/QHRKcvGfIZFZZXy9pMAl5Ic/ovko8BtPxc2zWnwCUks3CvAN4meQptQqvPs6g5TV9Xyml53uu+Tv05KTMzKyy/J2VmZoXlImVmZoXlImVmZoXlImVmZoXlImVmZoXlImVmZoXlImVmZoXlImVmZoXlImWWI0kXSnqoxrZzS6Yl3z1rnFaRdF1J/ye1uj/WmVykzPI1hmRKh1pdSzI1+ZMV4vy9vx0lXSnpkvTrvoIxo0K7i9Jts+voVy1OIem7WdO4SJnlazTwdB3tP4qIxbH2X4gfTT9FKv2r3N8D7ipZ/TowWdJGJe16gB8Cr9XRp5pExHsRsTjvuGalXKTMciJpC2Bz0pGUpI0k3SzpKUlfzhBnpaT7JH0o6d+SvlXSbFdgCPC3knXPAi8Ch5as249kVtm5ZceYm47ELpO0LF0uljSgpI0knS7pRUkr0mnUz6/1PMzy4CJllp+xwMfAQkmjgPkkUxN8IyJerTMOwIkkf417NPAca85afCBwb0R8WrbvNSTTnvc5muSWYqW/JP0Dkp8BE4AfA8cBPy3Zfh7wC+B8YAfgEJLRmtl64yJllp8xwD9ICsijwNURcWREfJwhznvAoRHxUES8BNwObFbS5gDgjxX2vQkYJ2lkOiLbB7iuynHeAk6OiBci4laSKb9PA5A0FDgV+HlEzIyIlyJiXkRcUee5mDWkp9UdMOsgY0imF58J7B8Rf2kgzj1l7/dsC7wEIGlbksn4HijfMSKWSbqTZAT1Lsl8Pq8lb2Gt5bFYc66eecCvJG1MMindYGBOxnMwy4VHUmb5GQP8ARgEbNpgnHll68by/6cGDwTmRMSHVfafCRxFUqhmZuxDxapmtr65SJnlQNKGJKOd35FMuX29pK83EKf8CcHSIlXtVl+fOSQznw5bR7vxWnOItTvwZkS8DzxPMiPrt2vvvVn+fLvPLB+jSR5OeC4inpC0PXCPpN0i4o0640DypB4AkjYFtgKekbQZSTGp+uHZiAhJOwOKiBX9HOsLwKWSrgB2An4G/DqNsVzSZcD5klYAD5OMDneJiN/WcT5mDXGRMsvHaODFkockfgmMAu6WtGdEfFRnnNJbeWOBVSSjmyOBJyJiSX9BImJ5Dce6ERgIPE5SYK8heZqwz1RgGckTflsBS4DrazsNs3xozfdNzWx9kTSXZOR1Uh373AU8EhEXre9j9xMrgEMi4vZGY5mV83tSZq11nKQPJO1aY/tHgN83s0O1Sj8M/EGr+2GdzSMpsxaR9EVgg/Tl6+t4/yjvY8+lwZGUpOHAxunLt/p52tAsMxcpMzMrLN/uMzOzwnKRMjOzwnKRMjOzwnKRMjOzwnKRMjOzwnKRMjOzwnKRMjOzwvofX/UmiEfltbwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "L.write()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see if this is real! This is the automatically generated file from the previous minimization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "config[\"data_file\"] = 'fake_sky_1.npy'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us try it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading data file: /Users/pierre/Documents/github/pybird/data/fake/fake_sky_1.npy\n",
      "-----------------------\n",
      "sky: sky_1\n",
      "output: bPk\n",
      "multipole: 3\n",
      "min bound (per multipole): [0.01, 0.01, 0.01]\n",
      "max bound (per multipole): [0.2, 0.2, 0.2]\n",
      "coordinate (AP) distortion: on\n",
      "binning: on\n",
      "-----------------------\n",
      "-----------------------\n",
      "get_maxlkl : True\n",
      "with_boss_correlated_skies_prior : False\n",
      "with_rs_marg : False\n",
      "-----------------------\n"
     ]
    }
   ],
   "source": [
    "L = Likelihood(config)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "init = L.loglkl(eft_init, eft_name, M, need_cosmo_update=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min chi2: 0.00, p-value: 1.000\n"
     ]
    }
   ],
   "source": [
    "L.c[\"get_maxlkl\"] = True # minimization mode\n",
    "\n",
    "def chi2(eft_params): return -2. * L.loglkl(eft_params, eft_name, M, need_cosmo_update=False)\n",
    "minimum = minuit(chi2, eft_init)\n",
    "\n",
    "from scipy import stats\n",
    "def pvalue(minchi2, dof): return 1. - stats.chi2.cdf(minchi2, dof)\n",
    "dof = L.y_all.shape[0] - len(eft_name) - L.Ng # number of data points - EFT parameters\n",
    "print ('min chi2: %.2f, p-value: %.3f' % (minimum['fun'], pvalue(minimum['fun'], dof)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}