testing_mcmc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I am using pyccl version 3.0.2\n"
     ]
    }
   ],
   "source": [
    "import pyccl as ccl\n",
    "print(f\"I am using pyccl version {ccl.__version__}\")\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scripts.presets import Presets\n",
    "from scripts.data_vectors import DataVectors\n",
    "from scripts.data_vector_metrics import DataVectorMetrics\n",
    "import emcee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecast_year = \"1\"\n",
    "presets = Presets(forecast_year=forecast_year, should_save_data=False)\n",
    "data = DataVectors(presets)\n",
    "cl_gc = np.array(data.galaxy_clustering_cls(include_all_correlations=True))\n",
    "cl_ggl = np.array(data.galaxy_galaxy_lensing_cls(include_all_correlations=True))\n",
    "cl_cs = np.array(data.cosmic_shear_cls(include_all_correlations=True))\n",
    "metric = DataVectorMetrics(presets)\n",
    "cls_noise = metric.get_matrix(cl_gc, cl_ggl, cl_cs, True)\n",
    "\n",
    "def likelihood(params):\n",
    "    Om = params[0]\n",
    "    Ob = params[1]\n",
    "    h = params[2]\n",
    "    sigma8 = params[3]\n",
    "    n_s = params[4]\n",
    "    w0 = params[5]\n",
    "    wa = params[6]\n",
    "    if Ob >= Om or Ob <= 0 or Om <= 0 or sigma8 <= 0:\n",
    "        return -np.inf\n",
    "    try:\n",
    "        cosmo = ccl.Cosmology(Omega_c=(Om-Ob),\n",
    "                                Omega_b=Ob,\n",
    "                                h=h,\n",
    "                                sigma8=sigma8,\n",
    "                                n_s=n_s,\n",
    "                                w0=w0,\n",
    "                                wa=wa)\n",
    "        \n",
    "        presets_mc = Presets(cosmology=cosmo, forecast_year=forecast_year, should_save_data=False)\n",
    "        \n",
    "        #presets_mc = Presets(cosmology=cosmo,redshift_resolution=600, forecast_year=forecast_year, should_save_data=False)\n",
    "        data = DataVectors(presets_mc)\n",
    "        cl_gc_mc = np.array(data.galaxy_clustering_cls(include_all_correlations=True))\n",
    "        cl_ggl_mc = np.array(data.galaxy_galaxy_lensing_cls(include_all_correlations=True))\n",
    "        cl_cs_mc = np.array(data.cosmic_shear_cls(include_all_correlations=True))\n",
    "        metric = DataVectorMetrics(presets_mc)\n",
    "        cls = metric.get_matrix(cl_gc_mc, cl_ggl_mc, cl_cs_mc, True)\n",
    "        return metric.get_loglike(cls, cls_noise)\n",
    "    except:\n",
    "        return -np.inf\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x17795c3b0>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGsCAYAAAB968WXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOB0lEQVR4nO3deVxU5eIG8GcWGBYZEJB9EUFQEPctcMtMLUtzLVOvmZlmZmpamvXLNs2yslVNS9NS29Qsr2upmeK+IsgmCqKIgDCswzLv7w9sruTGKDNnluf7+czn05x5OTzzXq/zeOac98iEEAJEREREJiKXOgARERHZFpYPIiIiMimWDyIiIjIplg8iIiIyKZYPIiIiMimWDyIiIjIplg8iIiIyKZYPIiIiMimWDyIiIjIplg8iIiIyKasoH/n5+XjhhRcQEREBJycnBAUFYfLkySgsLLztz82bNw8dOnSAi4sLvLy88NhjjyEpKanWGJlMdtPHBx98UGtcXFwcevbsCWdnZ7i5uaFHjx4oKyurNWbTpk3o1KkTHB0d4enpiUGDBt2QacWKFWjZsiUcHBzg4+ODSZMm6V+bM2fOTbM4OzsbNF/Z2dkYNWoUfHx84OzsjLZt2+Lnn382aB9ERER3y6LKR48ePbBixYobtl+8eBEXL17EggULcOrUKaxYsQJbtmzB2LFjb7u/3bt34/nnn8f+/fuxfft2VFVVoXfv3igpKdGPuXTpUq3HN998A5lMhsGDB+vHxMXFoW/fvujduzcOHjyIQ4cOYdKkSZDL/ze9v/zyC0aNGoUxY8bgxIkT2Lt3L5588slaeT766CPMnj0bM2fOxOnTp/HHH3+gT58++tenT59+Q57IyEgMHTrUoHkcNWoUkpKSsHHjRpw6dQqDBg3C448/jmPHjhm0HyIiorsiLEj37t3F8uXL6zT2xx9/FPb29qKysrLO+8/JyREAxO7du285ZsCAAaJnz561tnXq1Em89tprt/yZyspK4e/vL5YtW3bLMfn5+cLR0VHs2LGjznmPHz8uAIi//vqr1vaNGzeKtm3bCpVKJUJCQsScOXNqzYOzs7NYuXJlrZ9xd3e/bT4iIqL6YlFHPgxRWFgItVoNpVJp0M8AgLu7+01fv3z5MjZt2lTriEpOTg4OHDgALy8vxMTEwNvbG927d8fff/+tH3P06FFkZWVBLpejTZs28PX1xUMPPYTTp0/rx2zfvh06nQ5ZWVlo3rw5AgICMGzYMGRmZt4y77JlyxAeHo6uXbvqt23duhUjR47E5MmTkZCQgCVLlmDFihV499139WO6dOmCH374Afn5+dDpdFi7di20Wi169OhR57kiIiK6a1K3H0PU9chHbm6uCAoKErNnz67zvnU6nXj00UdFly5dbjlm/vz5omHDhqKsrEy/LS4uTgAQ7u7u4ptvvhFHjx4VU6ZMEfb29iI5OVkIIcSaNWsEABEUFCR+/vlncfjwYTF8+HDh4eEh8vLyhBBCzJs3T9jZ2YmIiAixZcsWERcXJx544AEREREhtFrtDVnKy8tFw4YNxfz582tt79q1q5g7d26tbatWrRK+vr765wUFBaJPnz4CgFAqlUKtVott27bVea6IiIjuhVmXj3fffVc4OzvrH3K5XKhUqlrb/v2VQ2FhoejUqZPo27evqKioqPPvmjhxoggODhaZmZm3HBMRESEmTZpUa9vevXsFADFr1qxa26Ojo8XMmTOFEEJ8//33AoBYsmSJ/vXy8nLh6ekpFi9erH+vAMTWrVv1Y3JycoRcLhdbtmy5Icvq1auFUqkUly5dqrXdyclJODg41JojBwcHAUCUlJQIIYSYNGmS6Nixo9ixY4c4fvy4mDNnjnB1dRUnT56sy1QRERHdk7p/JyGBCRMmYNiwYfrnI0aMwODBg2tdJeLv76//76KiIvTt2xcNGjTA+vXrYWdnV6ff88ILL2Djxo3466+/EBAQcNMxe/bsQVJSEn744Yda2319fQEAkZGRtbY3b94cGRkZtxyjUqnQpEmT245p1KgRPD099WOut2zZMjzyyCPw8fGptV2n0+HNN9+86ZU0Dg4OSEtLw+eff474+HhERUUBAFq1aoU9e/bgiy++wOLFi2/6/omIiOqLWZcPd3f3WudfODo6wsvLC2FhYTeM1Wg06NOnD1QqFTZu3AgHB4c77l8IgRdeeAHr16/Hrl27EBIScsuxX3/9Ndq1a4dWrVrV2t64cWP4+fndcIlucnIyHnroIQBAu3btoFKpkJSUhC5dugAAKisrce7cOQQHBwMAYmNjAQBJSUn6ApSfn4/c3Fz9mH+kp6dj586d2Lhx4w0527Zti6SkpJvOEQCUlpYCQK0rcQBAoVBAp9Pd8v0TERHVG6kPvRjiVud8aDQa0alTJxEdHS1SU1PFpUuX9I+qqir9uJ49e4rPPvtM//y5554Trq6uYteuXbV+prS0tNb+CwsLhZOTk1i0aNFNc3388cdCrVaLn376SaSkpIjXXntNODg4iNTUVP2YF198Ufj7+4utW7eKM2fOiLFjxwovLy+Rn5+vHzNgwAARFRUl9u7dK06dOiUeeeQRERkZecPXR6+99prw8/Or9d7+sWXLFqFUKsUbb7wh4uPjRUJCgli7dq3+/JeKigoRFhYmunbtKg4cOCBSU1PFggULhEwmE5s2bbrN7BMREdUPqygfO3fuFABu+khPT9ePCw4OFm+88Yb++a1+5t+/Y8mSJcLR0VEUFBTcMtu8efNEQECAcHJyEvfdd5/Ys2dPrdcrKirESy+9JLy8vISLi4vo1auXiI+PrzWmsLBQPP3008LNzU24u7uLgQMHioyMjFpjqqurRUBAgHj11VdvmWXLli0iJiZGODo6CrVaLTp27Ci++uor/evJycli0KBBwsvLSzg5OYmWLVvecOktERGRsciEEEKaYy5ERERki6x2nQ8iIiIyTywfREREZFJmd7WLTqfDxYsX4eLiAplMJnUcIiIiqgMhBIqKiuDn53fDFZX/Znbl4+LFiwgMDJQ6BhEREd2FzMzMW66Z9Q+zKx8uLi4AasKr1WqJ0xAREVFdaDQaBAYG6j/Hb8fsysc/X7Wo1WqWDyIiIgtTl1MmeMIpERERmRTLBxEREZkUywcRERGZFMsHERERmRTLBxEREZkUywcRERGZFMsHERERmRTLBxEREZkUywcRERGZFMsHERERmRTLBxEREZkUywcRERGZlNndWI7I1CqqdCgqr4SmvAqaskoUlVdBU14JTVklNOU1z7VVOvioHRDo7oRAd0cENnSCs4r/9yEiuhv825OsQn5JBfKKtTUF4lpx+F+JqNKXiOsLxT//XV6pu6vf6eFsjwB3JwQ2dESguxOC3J0Q2LCmnPi5OcJOwQOLREQ3w/JBFqmgtAL7z+ZhX1oe9qbmIu1KyT3vs4FKCbWDEmpHO7g4KKF2sIPa0Q5qByWUCjkuFZYhM78MGfmlKCyrRF5JBfJKKnAis+CGfcllgK+ro/4oyfVHTILcndDIRVWn204TEVkjlg+yCKUVVTh07ir2peZib1ouTl/UQIjaY9ycrisNDnZQOyrhcsN/K68Vipqxrtf+u4GDEgp53cuAprwSmfmlyMwvw4WrpcjML0VGfikyr9Y8L6/UIaugDFkFZdiP/Bt+XqWUI+DaEZPAhk4I926Ans294e/meK9TRURk9mRC/PuvcGlpNBq4urqisLAQarVa6jgkkYoqHY5nFmBfWi72pebhWOZVVFbX/qMa5tUAMaEeiAn1ROcm7nBzspcobW1CCFwp1iIzv+xaQSlF5tWaopJ5tRSXCstRrbv5/+2i/V3RO9IbvaN8EO7dgEdHiMhiGPL5zfJBZqFaJ5B4SYO9qbnYm5aHQ+n5KKusrjXG382xpmyE1RQOb7WDRGnvTWW1DpcKyq8VkpojJofPXcXh8/m4vpMEezihd6Q3+kT5oE1QQ4OOzBARmRrLB5k9IQTSrpQgLi0Xe1PzEHc2D4VllbXGuDvb475QD8SGeiI2zANB7k5WfSQgt1iLPxNzsC0hG3+l5KKi6n8nwno2sEev5t7oHeWNmFBPONgpJExKRHQjo5aPrKwsvPLKK9i8eTPKysoQHh6Or7/+Gu3atQMAFBcXY+bMmdiwYQPy8vLQuHFjTJ48Gc8991y9hyfLcrGgDHtTc7EvLQ/70nJxWaOt9XoDlRKdQtwRE+aJmFAPRHi7QG6j/9ov0Vbhr+Qr2JZwGX8kXoamvEr/mrO9Aj0ivNA7yhs9Irzg6mgnYVIiohqGfH4bdMLp1atXERsbi/vvvx+bN2+Gl5cX0tLS4Obmph8zdepU7Ny5E9999x0aN26Mbdu2YeLEifDz88OAAQPu6g2R5aqs1mHj8YtYuucszmQX1XrNXilH++CG175K8URLf1coeXkqAMBZpcRD0b54KNoXldU6HEzPx9bT2dh2+jKyNeXYdOoSNp26BKVchvtCPdA70hsPRvrAx9Uyv4oiItti0JGPmTNnYu/evdizZ88tx7Ro0QKPP/44Xn/9df22du3a4eGHH8bbb799x9/BIx/WobSiCj8cysSyPenIKigDUHP5acsAN8ReO2ejXXBDfn1gICEETmUVYtvpy9iWkI3ky8W1Xm8V6HbtPBFvhHm5SJSSiGyR0b52iYyMRJ8+fXDhwgXs3r0b/v7+mDhxIsaNG6cfM2HCBBw5cgQbNmyAn58fdu3ahf79+2Pz5s3o0qXLDfvUarXQav93+F2j0SAwMJDlw0JdLanAt3Hn8O2+c7haWnMOh2cDFcZ2CcGTHYPg6sSvCOpTem4Jtp3OxraEyziacbXW5cdNGjmjd6QPekd5o3WAm81+hUVEpmG08uHgUHNId9q0aRg6dCgOHjyIKVOmYMmSJfjPf/4DAKioqMC4ceOwcuVKKJVKyOVyLFu2DKNGjbrpPufMmYM333zzhu0sH5Ylq6AMy/acxdqDmfqrVII9nDC+WygGtfXnEQ4TyCkqxx+JOdh6Ohv7UvNQUf2/E1ZDPJ0xsUcoHmvjz5VXicgojFY+7O3t0b59e+zbt0+/bfLkyTh06BDi4uIAAAsWLMDSpUuxYMECBAcH46+//sKsWbOwfv169OrV64Z98siHZUu+XITFu9Ow8fhFVF27TrSFvxoTuofioRa+vDxUIkXlldidfAXbTl/GzjM5KNLWnLAa0NARE3uEYUi7ANgrWUKIqP4YrXwEBwfjwQcfxLJly/TbFi1ahHfeeQdZWVkoKyuDq6sr1q9fj379+unHPPPMM7hw4QK2bNlSr+FJOkfO52PRrjTsSMzRb4sN88CE7qHoEuZp1ZfEWpoSbRW+238eS/ecRW5xBQDAz9UBE3qEYlj7QB6VIqJ6YbSrXWJjY5GUlFRrW3JyMoKDgwEAlZWVqKyshFxe+19UCoUCOt3d3byLzIcQAjuTcrBoVxoOnbsKAJDJgL5RPpjQPRStAt2kDUg35axSYnz3UPznvsZYczADi3en4WJhOf7v19P4/M9UjO8eiic7BsHRniWEiEzDoCMfhw4dQkxMDN58800MGzYMBw8exLhx4/DVV19hxIgRAIAePXogNzcXn3/+OYKDg7F7924899xz+Oijj+q01gePfJifymodfj95EUt2/+9yWXuFHIPa+uPZbk3QpFEDiROSIcorq/HT4Uws2lVTQoCaRcye6doEozoHw1nFWz4RkeGMusjY77//jlmzZiElJQUhISGYNm1aratdsrOzMWvWLGzbtg35+fkIDg7Gs88+i6lTp9bpUDzLh/koq6jGD4cysPS6y2UbqJQY0SkIT3cJsdjlzalGRZUOvxy9gC93pSIzv+Z/34ZOdhjbJQT/iWkMtQOvTCKiuuPy6nRPCkorsDLuPFbsO4f8kppzBDwb2GNMbAhGdg7mippWprJah1+PX8QXO1ORnlsCAFA7KDEmNgRPx4bw8mgiqhOWD7orFwvK8PXf6VhzMAOlFTWXywa5O+HZbk0wpF0AT0y0clXVOmw6dQmf/ZmK1JyaxcsaqJQYHROMsV2awN3ZPO4aTETmieWDDFJaUYX3tyThu/3n9ZfLRvqqMaFHKB5u4cMlz22MTiewOT4bn/2Zoj/Hx8legZGdgzGuaxM0clFJnJCIzBHLB9XZkfP5eOnHEziXVwoAuK+JByb0CEW3prxc1tbpdALbEy/jsz9TEJ+lAQColHI82SkI47uF8j4yRFQLywfdkbaqGh9vT8FXf6VBJwBfVwfMH9wS3cIbSR2NzIwQAruSruCTP1JwPLMAQM3VTsM6BOC5HmHwd3OUNiARmQWWD7qt0xcLMe2HE0i6XHNIfXDbAPzfo5E8kZRuSwiBv1Nz8dkfqTh4Lh8AYKeQYXDbALzYqyl8XVlCiGwZywfdVFW1Dot2peGTP1JQpRPwbGCPdwdGo0+Uj9TRyMLsP5uHT/9Iwb60PACAi0qJV/s1xxMdAvl1HZGNYvmgG6TmFOOln07gxLXD5n2jfPDuwBbwaMCTB+nuHTmfj3c2JeJYRgEAoEuYJ+YNikagu5O0wYjI5Fg+SE+nE1ix7xzmbzkDbZUOagcl3hrQAgNa+/FfqFQvqnUCy/em44OtSdBW6eBkr8DMh5phZKdgyHljQSKbwfJBAIDM/FLM+PkE9p+t+X6+W3gjzB8cze/mySjSc0vwys8n9eeDdApxx/zBLdHY01niZERkCiwfNk4IgR8PZ+Lt3xNRrK2Ck70Crz7cHCM6BfFoBxmVTiewav95zN9yBqUV1XCwk2N67wiMiQ2BgkdBiKway4cNy9GUY+a6U/jzTM2t7js0bogFQ1sh2IP/+iTTycwvxSu/nNSfkNo2yA3vD2mFMC/ehJDIWrF82KjfTlzE67/Go6C0EvYKOab3CcfYLk34L06ShBACaw5mYu5/a47A2SvlmNorHOO6hnDVXCIrxPJhY66WVOD1X+Px+8lLAIAW/mp8NKw1wr1dJE5GBGQVlGHWulP4K/kKAKBVgCveH9IKET7880lkTVg+bMifZy7jlV9O4UqRFgq5DJPuD8OknmGw478syYwIIfDzkQt46/cEFJVXwU4hw+SeTTGhRyj/rBJZCZYPG1BUXol3fk/ED4czAQBhXg3w0bBWaBngJm0wotu4rCnH7PWnsCOx5pykSF81PhjaElF+rhInI6J7xfJh5eLS8jD9pxPIKiiDTAaMjQ3B9D4RvOU9WQQhBDaeuIg3Np5GQWkllHIZJvYIxfM9w6BS8s8wkaVi+bBS5ZXVmL/lDJbvPQcACHR3xIIhrdCpiYe0wYjuwpUiLf7v13hsjs8GAIR7N8AHQ1qhVaCbtMGI6K6wfFihM9kaTPz+KM5eKQEADO8YhNn9mqOBSilxMqJ7899Tl/D6hnjklVRALgOe7RaKKb2a8kgekYVh+bAy+8/mYdy3h1GkrYKXiwrzh7TE/RFeUsciqjf5JRWYs/E0Np64CABo0sgZHwxpiXbB7hInI6K6YvmwIlviszF57TFUVOnQsbE7loxqh4bO9lLHIjKKbaezMXtDPK4UaSGTAU/HhmB67wg42vMoCJG5M+Tzm9e4mbHvD5zHxO+PoKJKh96R3lg5tiOLB1m13lE+2DG1Owa3DYAQwNd/p6PvJ38hPqtQ6mhEVI9YPsyQEAKf7EjB7PXx0AlgeMdAfDmiLb8DJ5vg6mSHD4e1wvIxHeDr6oDzeaUYujgOm64tokdElo/lw8xU6wT+79fT+HhHMgBgcs8wzB0YzeWoyebcH+GFLVO6oVt4I5RVVuP51Ufx0fZk6HRm9U0xEd0FfqKZEW1VNV5YcxSr9p+HTAa8NSAK03pH8E60ZLNcHe3wzej2eKZLCADg0z9SMPH7oyjRVkmcjIjuBcuHmSgqr8RT3xzCf09lw14hx+fD2+I/9zWWOhaR5JQKOV57JBIfDGkJe4UcW05nY/CifcjML5U6GhHdJZYPM5BTVI7Hl+xH3Nk8NFApsWJMB/Rr6St1LCKzMrR9INY82wmeDVQ4k12EAV/sxcH0fKljEdFdYPmQ2LncEgxZFIeESxp4NrDH2mc7IybMU+pYRGapXbA7Nk6KRQt/NfJLKjBi2X6sOZghdSwiMhDLh4TiswoxZPE+ZOSXIsjdCT9PiEELf95gi+h2/Nwc8dP4GPRr6YvKaoFZ607hjV/jUVWtkzoaEdURy4dE9qbm4vElccgtrkCkrxo/P3cfGns6Sx2LyCI42ivw+fA2mN47HADwbdx5jF5+EAWlFRInI6K6YPmQwO8nL2LM8kMoqajGfU088MP4zvBycZA6FpFFkclkmNSzKZaMagcnewX2puZhwBd7kXK5SOpoRHQHLB8mtjLuHF5YcwwV1To8HO2DFU93gIuDndSxiCxWnygfrJsYg4CGjjifV4qBX+7Dn2cuSx2LiG6D5cNEhBD4cFsS/u/X0xACGNU5GJ8NbwuVkquWEt2rZj5q/Pp8LDqGuKNYW4Wx3x7G4t1pMLNbVxHRNSwfJlBVrcOr60/hsz9TAQDTHgzHWwOioJBz8TCi+uLRQIXvxnbCk52CIATw3uYzmPbjCZRXVksdjYj+heXDyMorqzHx+6NYczATchkwd2A0Jj/QlKuWEhmBvVKOuQOj8fa1cr/+WBYe/2o/LmvKpY5GRNdh+TCiwrJK/Ofrg9iWcBn2Sjm+HNEWT3YKkjoWkdUbdV9jrHq6I9yc7HAiswD9P/8bJzILpI5FRNewfBjJZU05Hl8Sh4Pn8uGiUmLl0x3RtwVXLSUylZgwT/z6fCyaejXAZY0WQ5fE4dfjWVLHIiKwfBjF2SvFGPTlPpzJLkIjFxV+GH8fOjfxkDoWkc0J9nDGuokx6NXcCxVVOry49jjmbznDO+MSSYzlo56dyCzAkMVxyCooQ4inM9Y9F4NIP7XUsYhslouDHZaMao+JPUIBAIt2peHZVYdRVF4pcTIi22Vw+cjKysLIkSPh4eEBJycntG7dGkeOHKk1JjExEf3794erqytcXFzQuXNnZGRY//0X/kq+guFL9yO/pAItA1zx84T7EOjuJHUsIpunkMvwct9m+OSJ1lAp5diRmINBX+7D+bwSqaMR2SSDysfVq1cRGxsLOzs7bN68GQkJCfjwww/h5uamH5OWloYuXbqgWbNm2LVrF06cOIHXX38dDg7WvYLnrqQcjP32EEorqtG1qSdWj+sMjwYqqWMR0XUGtPbHj+Pvg7dahZScYgz4Yi/2peZKHYvI5siEAavwzJw5E3v37sWePXtuOeaJJ56AnZ0dVq1adVeBNBoNXF1dUVhYCLXaMr6uiM8qxONL4lBSUY1+0b74+PHWsFfyGy0ic5WjKce4VUdwIrMASrkMnz/ZhieEE90jQz6/DfqE3LhxI9q3b4+hQ4fCy8sLbdq0wdKlS/Wv63Q6bNq0CeHh4ejTpw+8vLzQqVMnbNiw4Zb71Gq10Gg0tR6WJKugDE+vqLlPS0yoB4sHkQXwUjvgh2c745GWvqjSCTy/+hh+O3FR6lhENsOgT8mzZ89i0aJFaNq0KbZu3YoJEyZg8uTJWLlyJQAgJycHxcXFeO+999C3b19s27YNAwcOxKBBg7B79+6b7nPevHlwdXXVPwIDA+/9XZlIYVklxiw/iJwiLSK8XbB4VDsWDyIL4WCnwCdPtMGgtv6o1gm8uPYYNhzjpbhEpmDQ1y729vZo37499u3bp982efJkHDp0CHFxcbh48SL8/f0xfPhwrF69Wj+mf//+cHZ2xpo1a27Yp1arhVar1T/XaDQIDAw0+69dtFXVGP3NQew/mw9vtQrrJ8bCz81R6lhEZCCdTmDWulP44XAmZDLggyGtMKRdgNSxiCyO0b528fX1RWRkZK1tzZs311/J4unpCaVSedsx/6ZSqaBWq2s9zJ0QAq/8fBL7z+ajgUqJ5U91ZPEgslByuQzzBkVjxLV7wsz4+QTWHLT+q/OIpGRQ+YiNjUVSUlKtbcnJyQgODgZQc2SkQ4cOtx1jDRZsS8KG4xehlMvw5Yi2XMeDyMLJ5TK881gLPBXTGEIAs9adwqq4c1LHIrJaSkMGT506FTExMZg7dy6GDRuGgwcP4quvvsJXX32lHzNjxgw8/vjj6NatG+6//35s2bIFv/32G3bt2lXf2SWx+kAGvtiZBgCYOyga3cIbSZyIiOqDTCbDG49GQimXYdnf6Xj919OorBZ4ukuI1NGIrI5B53wAwO+//45Zs2YhJSUFISEhmDZtGsaNG1drzDfffIN58+bhwoULiIiIwJtvvokBAwbUaf/mfKntzjM1a3noBPDiA00x9cFwqSMRUT0TQuD9rUlYtKvmHxmvPtwMz3YLlTgVkfkz5PPb4PJhbOZaPk5dKMTjX8WhtKIaQ9oF4IMhLSGTyaSORURGIITAx9uT8emfqQCAGX0i8Pz9YRKnIjJvRjvh1FZl5pdizIr/rV46b1A0iweRFZPJZJjWOwLTrh3d/GBrEj7ZkSJxKiLrwfJxB4WllXhq+UHkFmvRzMcFX45oCzsFp43IFkx+oCle7hsBAPh4RzI+3JYEMztYTGSR+Cl6G9qqaoxbdRhpV0rg6+qAFWM6wsXBTupYRGRCE3uE4bV+zQEAn/2Zive2nGEBIbpHLB+3oNMJTP/pJA6m58NFpcTyMR3g42rdN8cjopt7pmsTzHm0Zv2iJbvP4u3fE1lAiO4By8ctvL81Cb+dqFnLY/GodmjmYz4nvxKR6T0VG4J3HmsBAPhmbzre2HgaOh0LCNHdYPm4iVX7z2Px7prL7OYPbonYME+JExGRORjZORjzB0dDJgNWxp3H7A3xLCBEd4Hl4192JFzGG7/GAwBeejAcg3mPByK6zuMdgvDBkFaQyYA1BzPwyi8nUc0CQmQQlo/rnMgswAtrjkEngCc6BGJST17XT0Q3GtIuAAsfbw25DPjpyAVM/+kEqqp1UscishgsH9dk5JVi7LeHUFZZje7hjfD2Yy24lgcR3dKA1v74bHhbKOQyrD+Whak/soAQ1RXLB4CrJRV4asVB5BZXIMpPjS+4lgcR1UG/lr744sm2sFPI8NuJi3hhzTFUsoAQ3ZHNf8KWV1Zj3MrDOHulBP5ujvjmqQ5ooDLofntEZMP6tvDBohHtYK+QY3N8NiZ+fxTaqmqpYxGZNZsuHzqdwEs/nsDh81fh4lCzloe3mmt5EJFhekV646v/tIO9Uo7tCZfx3HdHUV7JAkJ0KzZdPuZtTsSmU5dgp5Bhyah2CPd2kToSEVmoHhFe+GZ0BzjYyfHnmRyMW3mYBYToFmy2fKzYm46le9IBAAuGtkJMKNfyIKJ706WpJ5Y/1RGOdgrsScnF0ysOobSiSupYRGbHJsvHttPZePP3BAA1t8oe0Npf4kREZC3uC/XAt093hLO9AvvS8vDMt4d5DgjRv9hc+TiWcRWT1x6DEMCTnYIwsUeo1JGIyMp0DHHHyrGd9AVkxk8nuRIq0XVsqnyczyvB2G8Po7xSh57NvPBW/yiu5UFERtEuuCEWjWwHpVyGjScuYv7WM1JHIjIbNlM+8ksq8NTyQ8gvqUC0vys+G94GSq7lQURG1C28Ed4b3BJAzd1wv913TtpARGbCZj59i8uroBMC/m6O+Pqp9nDmWh5EZAJD2gVgeu9wAMCc305jS3y2xImIpCcTQpjVF5EajQaurq4oLCyEWl2/t7HPK9aisKwSTRo1qNf9EhHdjhACr66Px5qDGVAp5Vg9rhPaBbtLHYuoXhny+W0zRz4AwKOBisWDiExOJpPh7QFR6NXcC9oqHcZ+exhpV4qljkUkGZsqH0REUlEq5Ph0eBu0CnRDQWklRn9zEDlF5VLHIpIEywcRkYk42Svx9ej2CPZwwoWrZXh6xSGUaLkIGdkelg8iIhPybKDCt2M6wt3ZHvFZGkz8/ijvhEs2h+WDiMjEGns64+vR7eFgJ8fu5CuYvf4UzOzcfyKjYvkgIpJAm6CG+Hx4W8hlwI+HL2DhjhSpIxGZDMsHEZFEekV64+3HWgAAPvkjBWsPZkiciMg0WD6IiCQ0olMwJt0fBgCYvSEeO8/kSJyIyPhYPoiIJPZS73AMauuPap3AxO+P4uSFAqkjERkVywcRkcRkMhneG9QSXZt6oqyyGk+vOISMvFKpYxEZDcsHEZEZsFfK8eWItoj0VSO3uAKjlx9EfkmF1LGIjILlg4jITLg42GH5mA7wd3NEem4Jxn57CGUV1VLHIqp3LB9ERGbEW+2Ab5/uALWDEscyCvDi2mOo1nENELIuLB9ERGYmzMsFy0Z3gL1Sjm0JlzFn42kuQkZWheWDiMgMdQxxx8LHW0MmA1btP4/Fu89KHYmo3rB8EBGZqYejffFav0gAwPwtZ7DhWJbEiYjqB8sHEZEZG9slBM90CQEAzPj5BPam5kqciOjesXwQEZm5Vx9ujn4tfVFZLTBh1REkXtJIHYnonhhcPrKysjBy5Eh4eHjAyckJrVu3xpEjR246dvz48ZDJZFi4cOG95iQisllyuQwfDm2FjiHuKNJW4anlB3GxoEzqWER3zaDycfXqVcTGxsLOzg6bN29GQkICPvzwQ7i5ud0wdsOGDThw4AD8/PzqKysRkc1ysFNg6aj2aOrVAJc1Wjy1/CAKyyqljkV0VwwqH/Pnz0dgYCCWL1+Ojh07onHjxnjggQcQGhpaa1xWVhYmTZqE77//HnZ2dvUamIjIVrk62WHF0x3hrVYh+XIxnl15GNoqLkJGlseg8rFx40a0b98eQ4cOhZeXF9q0aYOlS5fWGqPT6TBq1CjMmDEDUVFRd9ynVquFRqOp9SAiopvzd3PE8qc6ooFKiQPp+Xh1XTzXACGLY1D5OHv2LBYtWoSmTZti69atmDBhAiZPnoyVK1fqx8yfPx9KpRKTJ0+u0z7nzZsHV1dX/SMwMNCwd0BEZGMi/dRYNLIt5DLgl6MX8P2BDKkjERlEJgyozPb29mjfvj327dun3zZ58mQcOnQIcXFxOHLkCPr164ejR4/qz/Vo3LgxpkyZgilTptx0n1qtFlqtVv9co9EgMDAQhYWFUKvVd/m2iIis3+LdaXhv8xnYKWT4Yfx9aBvUUOpIZMM0Gg1cXV3r9Plt0JEPX19fREZG1trWvHlzZGTUtO49e/YgJycHQUFBUCqVUCqVOH/+PF566SU0btz4pvtUqVRQq9W1HkREdGfjuzXBQy18UFktMPG7o7hSpL3zDxGZAYPKR2xsLJKSkmptS05ORnBwMABg1KhROHnyJI4fP65/+Pn5YcaMGdi6dWv9pSYiIshkMnwwtBVCGzkjW1OOF9YcRVW1TupYRHdkUPmYOnUq9u/fj7lz5yI1NRWrV6/GV199heeffx4A4OHhgRYtWtR62NnZwcfHBxEREUZ5A0REtqyBSoklo9rD2V6B/WfzMX/LGakjEd2RQeWjQ4cOWL9+PdasWYMWLVrg7bffxsKFCzFixAhj5SMiojsI82qABUNbAQCW7knH7ycvSpyI6PYMOuHUFAw5YYWIiP5n3n8TseSvs3CyV+DX52PR1NtF6khkQ4x2wikREZmvGX0icF8TD5RWVGP8qiMoKucKqGSeWD6IiKyEUiHHZ0+2ga+rA87mlmD6Tye4ABmZJZYPIiIr4tlAhUUj28FeIcfW05exaHea1JGIbsDyQURkZVoHumFO/5rbWyzYmoS/U3IlTkRUG8sHEZEVGt4xEMPaB0AngBfWHEVWQZnUkYj0WD6IiKyQTCbDWwNaoIW/GldLK/Hcd0dQXsk74JJ5YPkgIrJSDnYKLBrRDm5Odjh5oRBzNp6WOhIRAJYPIiKrFujuhE+faAOZDFh7KBNrD/IOuCQ9lg8iIivXLbwRpveuucXF//16GicyC6QNRDaP5YOIyAY81z0UD0Z6o6Jah+e+O4L8kgqpI5ENY/kgIrIBcrkMHw5rhRBPZ1wsLMfkNcdQreMCZCQNlg8iIhuhdrDD4pHt4GinwN+puViwLUnqSGSjWD6IiGxIhI8L3h/SEgCwaFcatsRnS5yIbBHLBxGRjXm0lR/GdgkBAEz/6QRSc4olTkS2huWDiMgGzXyoGTqGuKNYW4UJ3x1BsbZK6khkQ1g+iIhskJ1Cji+ebAtvtQqpOcV45eeTvAMumQzLBxGRjWrkosKXI9rBTiHDplOXsGxPutSRyEawfBAR2bB2wQ3xf49EAgDe23IGcWl5EiciW8DyQURk40Z2DsagNv6o1glMWn0Ulwp5B1wyLpYPIiIbJ5PJ8O7AaDT3VSOvpALPfXcU2ireAZeMh+WDiIjgaK/AkpHt4Opoh+OZBXj79wSpI5EVY/kgIiIAQJCHExY+0RoyGfDd/gz8fOSC1JHISrF8EBGR3v0RXpjyQDgAYPb6U0i8pJE4EVkjlg8iIqrlhZ5h6NnMC9oqHV5cewzllTz/g+oXywcREdUil8vwwZCW8GygQvLlYry3+YzUkcjKsHwQEdENPBqo8MHQmhvQrdh3DruTr0iciKwJywcREd3U/RFeGH1fMICaG9DlFWslTkTWguWDiIhuadbDzdHUqwGuFGkxc90p3v+F6gXLBxER3ZKDnQILn2gNO4UM2xMuY+2hTKkjkRVg+SAiotuK8nPFjD4RAIC3fkvA2SvFEiciS8fyQUREd/RMlyaICfVAWWU1pv5wHJXVOqkjkQVj+SAiojuSy2X4cFgruDra4cSFQnyyI0XqSGTBWD6IiKhOfF0dMXdgNADgy12pOJieL3EislQsH0REVGf9WvpicNsA6AQw9Yfj0JRXSh2JLBDLBxERGWRO/0gEuTshq6AM/7chXuo4ZIFYPoiIyCAuDnb4+PHWUMhl2HD8In49niV1JLIwLB9ERGSwdsENMen+MADAaxviceFqqcSJyJKwfBAR0V15oWcY2gS5oai8CtN+PIFqHVc/pboxuHxkZWVh5MiR8PDwgJOTE1q3bo0jR44AACorK/HKK68gOjoazs7O8PPzw3/+8x9cvHix3oMTEZG0lAo5Fj7eGs72ChxMz8eSv9KkjkQWwqDycfXqVcTGxsLOzg6bN29GQkICPvzwQ7i5uQEASktLcfToUbz++us4evQo1q1bh+TkZPTv398Y2YmISGLBHs54o38UAOCjbck4daFQ4kRkCWTCgLsEzZw5E3v37sWePXvq/AsOHTqEjh074vz58wgKCrrjeI1GA1dXVxQWFkKtVtf59xARkTSEEJj4/VFsjs9Gk0bO2PRCVzjaK6SORSZmyOe3QUc+Nm7ciPbt22Po0KHw8vJCmzZtsHTp0tv+TGFhIWQymf7oyL9ptVpoNJpaDyIishwymQxzB0bDW63C2SsleGdTgtSRyMwZVD7Onj2LRYsWoWnTpti6dSsmTJiAyZMnY+XKlTcdX15ejpkzZ+LJJ5+8ZQuaN28eXF1d9Y/AwEDD3wUREUmqobM9PhzaGgDw/YEM7Ei4LG0gMmsGfe1ib2+P9u3bY9++ffptkydPxqFDhxAXF1drbGVlJYYOHYqMjAzs2rXrluVDq9VCq9Xqn2s0GgQGBvJrFyIiC/TO7wlY9nc6PJztsWVKNzRyUUkdiUzEaF+7+Pr6IjIysta25s2bIyMjo9a2yspKDBs2DOnp6di+ffttQ6hUKqjV6loPIiKyTDP6RqCZjwvySirw8s8nYMC/b8mGGFQ+YmNjkZSUVGtbcnIygoOD9c//KR4pKSnYsWMHPDw86icpERGZPZVSgU+Ht4G9Uo6dSVewav95qSORGTKofEydOhX79+/H3LlzkZqaitWrV+Orr77C888/DwCoqqrCkCFDcPjwYXz//feorq5GdnY2srOzUVFRYZQ3QERE5iXc2wWzHmoGAHh3UyJSLhdJnIjMjUHnfADA77//jlmzZiElJQUhISGYNm0axo0bBwA4d+4cQkJCbvpzO3fuRI8ePe64f15qS0Rk+YQQGL38EP5KvoLmvmpseD4GKiUvv7Vmhnx+G1w+jI3lg4jIOuRoytH3kz3IL6nAs92a4NWHm0sdiYzIaCecEhER1ZWX2gHvDYoGACzdcxb7UnMlTkTmguWDiIiMpneUD4Z3DIQQwLQfT6CglOf/EcsHEREZ2euPRKKJpzOyNeWYvT6el98SywcRERmXk70SC59oDaVchk2nLuGXo1lSRyKJsXwQEZHRtQxww9QHwwEAb/waj4y8UokTkZRYPoiIyCQmdA9Fx8buKKmoxpQfjqGqWid1JJIIywcREZmEQi7DR4+3gotKiaMZBfhiZ5rUkUgiLB9ERGQyAQ2d8PZjLQAAn/6ZgvisQokTkRRYPoiIyKQea+OPh6N9UK0TmLnuJL9+sUEsH0REZHJz+kfB1dEO8VkaLPs7Xeo4ZGIsH0REZHJeLg54rV/Ncusfb09Gem6JxInIlFg+iIhIEkPaBaBrU09oq3SY+ctJ6HRcfMxWsHwQEZEkZDIZ5g6MhqOdAgfS87H2UKbUkchEWD6IiEgyge5OmN4nAgAw77+JyC4slzgRmQLLBxERSeqpmMZoHeiGIm0VXttwivd+sQEsH0REJCmFXIb5g1vCTiHDjsQc/H7yktSRyMhYPoiISHIRPi6Y2CMMADBn42lcLamQOBEZE8sHERGZhYn3h6KpVwPklVTg7U0JUschI2L5ICIis6BSKjB/SEvIZMC6o1nYlZQjdSQyEpYPIiIyG22DGuKpmMYAgNnr41GirZI2EBkFywcREZmV6b0j4O/miKyCMnywNUnqOGQELB9ERGRWnFVKzBsUDQD4Nu4cjpy/KnEiqm8sH0REZHa6hTfC4LYBEAJ45ZeT0FZVSx2J6hHLBxERmaXXH2kOzwb2SM0pxhc706SOQ/WI5YOIiMySm5M93uzfAgDw5c5UnMnWSJyI6gvLBxERma2Ho33wYKQ3qnQCr/x8EtW8861VYPkgIiKzJZPJ8M5jLeDioMSJC4VYvjdd6khUD1g+iIjIrHmrHfDqw80BAAu2JSEjr1TiRHSvWD6IiMjsPdEhEJ2buKO8UodZ60/yzrcWjuWDiIjMnkwmw3uDWkKllGNvah5+OnxB6kh0D1g+iIjIIjT2dMa0B8MBAO9sSkCOplziRHS3WD6IiMhijO0Sgmh/V2jKq/DGxtNSx6G7xPJBREQWQ6mQY/7gllDKZdgcn40t8ZekjkR3geWDiIgsSqSfGuO7NwEAvP7raRSWVkqciAzF8kFERBbnhZ5N0aSRM64UaTH3v4lSxyEDsXwQEZHFcbBTYP7glgCAHw5nYm9qrsSJyBAsH0REZJE6NHbHqM7BAIBZ606htKJK4kRUVywfRERksV7uGwFfVwdk5Jfio23JUsehOjK4fGRlZWHkyJHw8PCAk5MTWrdujSNHjuhfF0Jgzpw58PPzg6OjI3r06IHTp3k5FBER1T8XBzu8O7Dmzrff7E3HicwCaQNRnRhUPq5evYrY2FjY2dlh8+bNSEhIwIcffgg3Nzf9mPfffx8fffQRPv/8cxw6dAg+Pj548MEHUVRUVN/ZiYiI0LOZNwa09oNOAK/8chIVVTqpI9EdyIQBC+TPnDkTe/fuxZ49e276uhACfn5+mDJlCl555RUAgFarhbe3N+bPn4/x48ff8XdoNBq4urqisLAQarW6rtGIiMiG5RVr8eDHfyG/pALTHgzH5AeaSh3J5hjy+W3QkY+NGzeiffv2GDp0KLy8vNCmTRssXbpU/3p6ejqys7PRu3dv/TaVSoXu3btj3759N92nVquFRqOp9SAiIjKERwMV3ng0EgDw+Z+pSM3h0XZzZlD5OHv2LBYtWoSmTZti69atmDBhAiZPnoyVK1cCALKzswEA3t7etX7O29tb/9q/zZs3D66urvpHYGDg3bwPIiKycf1b+eH+iEaoqNbh5Z9PolrHO9+aK4PKh06nQ9u2bTF37ly0adMG48ePx7hx47Bo0aJa42QyWa3nQogbtv1j1qxZKCws1D8yMzMNfAtEREQ1nz3vDoyGs70CRzMKsCrunNSR6BYMKh++vr6IjIysta158+bIyMgAAPj4+ADADUc5cnJybjga8g+VSgW1Wl3rQUREdDf83Bwx86FmAIAF25KRU8Q735ojg8pHbGwskpKSam1LTk5GcHDNIi8hISHw8fHB9u3b9a9XVFRg9+7diImJqYe4REREtzeiUzBaBbiiWFuF9/57Ruo4dBMGlY+pU6di//79mDt3LlJTU7F69Wp89dVXeP755wHUHPKaMmUK5s6di/Xr1yM+Ph5PPfUUnJyc8OSTTxrlDRAREV1PLpfhzQEtIJMB645l4dC5fKkj0b8YVD46dOiA9evXY82aNWjRogXefvttLFy4ECNGjNCPefnllzFlyhRMnDgR7du3R1ZWFrZt2wYXF5d6D09ERHQzrQPd8Hj7mgsYXt8Qj6pqrv1hTgxa58MUuM4HERHVh/ySCty/YBcKyyrxZv8ojI5pLHUkq2a0dT6IiIgshbuzPab3DgcALNiWhNxircSJ6B8sH0REZLWe7BSMKD81isqr8P4WnnxqLlg+iIjIainkMrw1IAoA8OPhCziacVXiRASwfBARkZVrF+yOwW0DAAD/92s8Vz41AywfRERk9WY+1AwuDkrEZ2mw9lCG1HFsHssHERFZvUYuKkx7sObk0w+2JuFqSYXEiWwbywcREdmEUZ2D0czHBQWllXh/a9Kdf4CMhuWDiIhsglIhx5v9a04+XXsoAycvFEgbyIaxfBARkc3o1MQDj7X2gxDA//16GjqefCoJlg8iIrIpsx5uDmd7BY5nFuCnI5lSx7FJLB9ERGRTvNUOmNKr5uTT+VuSUFhaKXEi28PyQURENuep2MYI82qA/JIKfLidJ5+aGssHERHZHDuFHG9dO/n0u/3ncfpiocSJbAvLBxER2aSYME/0a+kL3bWTT83sJu9WjeWDiIhs1mv9msPRToEj569i3dEsqePYDJYPIiKyWb6ujnjhgTAAwLzNZ6Ap58mnpsDyQURENu2ZLk3QxNMZucVaLNyeInUcm8DyQURENs1eKcecayeffht3DmeyNRInsn4sH0REZPO6hTdCnyhvVOsETz41AZYPIiIiAK8/EgkHOzkOpudj44mLUsexaiwfREREAAIaOuH5HjUnn879byKKtVUSJ7JeLB9ERETXjOvWBMEeTris0eLTP3jyqbGwfBAREV3jYKfAG49GAgC++TsdqTlFEieyTiwfRERE1+nZzBu9mnuhSifwxkaefGoMLB9ERET/8n+PRMFeKcfe1Dz891S21HGsDssHERHRvwR5OGFC91AAwDubElBawZNP6xPLBxER0U1M7BGKgIaOuFRYjs//TJU6jlVh+SAiIroJBzsFXn+k5uTTpXvO4uyVYokTWQ+WDyIiolvoHemN7uGNUFktMOe3BJ58Wk9YPoiIiG5BJpNhTv8o2Cvk+Cv5CrYlXJY6klVg+SAiIrqNEE9nPNM1BADw1m8JKKuoljiR5WP5ICIiuoNJPcPg5+qArIIyLNqdJnUci8fyQUREdAdO9kq8du3k08W703A+r0TiRJaN5YOIiKgOHmrhg9gwD1RU6fDWbwlSx7FoLB9ERER1IJPJ8Gb/KCjlMvxxJgd/nuHJp3eL5YOIiKiOwrxc8HSXmpNP392UiMpqncSJLBPLBxERkQEm9QxDQyc7pF0pwdpDmVLHsUgsH0RERAZQO9hhSq9wAMDC7ckoKq+UOJHlMah8zJkzBzKZrNbDx8dH/3pxcTEmTZqEgIAAODo6onnz5li0aFG9hyYiIpLSk52CEOLpjLySCizmpbcGM/jIR1RUFC5duqR/nDp1Sv/a1KlTsWXLFnz33XdITEzE1KlT8cILL+DXX3+t19BERERSslPIMfOhZgCAZXvScbGgTOJElsXg8qFUKuHj46N/NGrUSP9aXFwcRo8ejR49eqBx48Z49tln0apVKxw+fLheQxMREUmtd6Q3Ooa4Q1ulw4JtSVLHsSgGl4+UlBT4+fkhJCQETzzxBM6ePat/rUuXLti4cSOysrIghMDOnTuRnJyMPn363HJ/Wq0WGo2m1oOIiMjcyWQyzH64OQBg3dEsxGcVSpzIchhUPjp16oSVK1di69atWLp0KbKzsxETE4O8vDwAwKefforIyEgEBATA3t4effv2xZdffokuXbrccp/z5s2Dq6ur/hEYGHhv74iIiMhEWgW6YUBrPwDAO5t419u6kol7mKmSkhKEhobi5ZdfxrRp07BgwQIsXboUCxYsQHBwMP766y/MmjUL69evR69evW66D61WC61Wq3+u0WgQGBiIwsJCqNXqu41GRERkEpn5pXjgo92oqNLh69Ht8UBzb6kjSUKj0cDV1bVOn9/Ke/lFzs7OiI6ORkpKCsrKyvDqq69i/fr16NevHwCgZcuWOH78OBYsWHDL8qFSqaBSqe4lBhERkWQC3Z0wJrYxluw+i7n/TUT38EZQKriSxe3c0+xotVokJibC19cXlZWVqKyshFxee5cKhQI6HVeAIyIi6/X8/Vx4zBAGlY/p06dj9+7dSE9Px4EDBzBkyBBoNBqMHj0aarUa3bt3x4wZM7Br1y6kp6djxYoVWLlyJQYOHGis/ERERJK7fuGxj7nw2B0ZVD4uXLiA4cOHIyIiAoMGDYK9vT3279+P4OBgAMDatWvRoUMHjBgxApGRkXjvvffw7rvvYsKECUYJT0REZC648Fjd3dMJp8ZgyAkrRERE5mTr6WyMX3UEKqUcO6f3gJ+bo9SRTMaQz2+eEUNERFRPekd6o2NjLjx2JywfRERE9UQmk2F2v5qFx9Yf48Jjt8LyQUREVI/+WXhMCODdTYlceOwmWD6IiIjq2fTeEbBXyhF3Ng9/nsmROo7ZYfkgIiKqZ/8sPAYAc/+biKpqrnd1PZYPIiIiI5jYgwuP3QrLBxERkRG4Ov5v4bGFO7jw2PVYPoiIiIzkn4XHcosrsGT3WanjmA2WDyIiIiOxU8gx86FmAICle87iYkGZxInMA8sHERGREXHhsRuxfBARERkRFx67EcsHERGRkXHhsdpYPoiIiEzg+oXHdibZ9sJjLB9EREQmUHvhsTM2vfAYywcREZGJ/LPwWGpOsU0vPMbyQUREZCJceKwGywcREZEJceExlg8iIiKT+vfCY5cKbW/hMZYPIiIiE7t+4bEPttrewmMsH0RERCZm6wuPsXwQERFJwJYXHmP5ICIikoitLjzG8kFERCQRW114jOWDiIhIQtcvPPbDYdtYeIzlg4iISEKujnZ48YGmAICPt9vGwmMsH0RERBIb0TnYphYeY/kgIiKSmK0tPMbyQUREZAauX3hswdZkqeMYFcsHERGRGZDJZHj12sJj645dwOmL1rvwGMsHERGRmWgd6IZHWvpCCOCjbdZ79IPlg4iIyIxMezAcCrkMf5zJwdGMq1LHMQqWDyIiIjPSpFEDDG7rDwBYYKU3nWP5ICIiMjOTH2gKO4UM+9LysDc1V+o49Y7lg4iIyMwENHTCiE7BAIAPtiZZ3U3nWD6IiIjM0MT7Q+FgJ8fxzAL8kWhdN51j+SAiIjJDXi4OeComBACwYFsSdDrrOfrB8kFERGSmJnRvAheVEmeyi7Dp1CWp49Qblg8iIiIz5eZkj2e6NgFQc9O5qmqdxInqh0HlY86cOZDJZLUePj4+tcYkJiaif//+cHV1hYuLCzp37oyMjIx6DU1ERGQrxnYNgbuzPc7mlmDd0Syp49QLg498REVF4dKlS/rHqVOn9K+lpaWhS5cuaNasGXbt2oUTJ07g9ddfh4ODQ72GJiIishUNVEo81z0UAPDJHynQVlVLnOjeKQ3+AaXyhqMd/5g9ezYefvhhvP/++/ptTZo0uft0REREhFH3BWPZ32eRVVCGtQczMTqmsdSR7onBRz5SUlLg5+eHkJAQPPHEEzh79iwAQKfTYdOmTQgPD0efPn3g5eWFTp06YcOGDbfdn1arhUajqfUgIiKi/3GwU2BSz6YAgM93pqKswrKPfhhUPjp16oSVK1di69atWLp0KbKzsxETE4O8vDzk5OSguLgY7733Hvr27Ytt27Zh4MCBGDRoEHbv3n3Lfc6bNw+urq76R2Bg4D2/KSIiImvzePtABDR0xJUiLb6NOyd1nHsiE/ewbFpJSQlCQ0Px8ssv44knnoC/vz+GDx+O1atX68f0798fzs7OWLNmzU33odVqodVq9c81Gg0CAwNRWFgItVp9t9GIiIiszi9HLuCln07A1dEOe165H2oHO6kj6Wk0Gri6utbp8/ueLrV1dnZGdHQ0UlJS4OnpCaVSicjIyFpjmjdvfturXVQqFdRqda0HERER3eixNv4I82qAwrJKLNuTLnWcu3ZP5UOr1SIxMRG+vr6wt7dHhw4dkJRU+w58ycnJCA4OvqeQREREBCjkMkx7MBwA8PWes8gvqZA40d0xqHxMnz4du3fvRnp6Og4cOIAhQ4ZAo9Fg9OjRAIAZM2bghx9+wNKlS5GamorPP/8cv/32GyZOnGiU8ERERLamb5QPovzUKKmoxuLdaVLHuSsGlY8LFy5g+PDhiIiIwKBBg2Bvb4/9+/frj2wMHDgQixcvxvvvv4/o6GgsW7YMv/zyC7p06WKU8ERERLZGLpdheu8IAMC3+87hsqZc4kSGu6cTTo3BkBNWiIiIbJEQAkMXx+Hw+asY2TkI7zwWLXUk051wSkRERKYnk8kwvU/N0Y+1BzORmV8qcSLDsHwQERFZoM5NPNC1qSeqdAILd6RIHccgLB9EREQW6qVr536sP3YBqTlFEqepO5YPIiIiC9U60A0PRnpDJ4CPt1vO0Q+WDyIiIgv2Uu9wyGTAplOXEJ9VKHWcOmH5ICIismDNfNTo38oPAPDhtqQ7jDYPLB9EREQWbmqvcCjkMuxMuoIj5/OljnNHLB9EREQWrrGnM4a2CwAAfLA1CWa2hNcNWD6IiIiswAsPNIW9Qo79Z/OxNzVP6ji3xfJBRERkBfzdHDGicxAA4IOtZ8z66AfLBxERkZWY2CMMjnYKnLhQiO0Jl6WOc0ssH0RERFaikYsKY2IbAwA+2p4Mnc48j36wfBAREVmR8d1C4eKgxJnsIvx28qLUcW6K5YOIiMiKuDrZ4dmuTQAAH29PRmW1TuJEN2L5ICIisjJjuoTAw9ke5/JK8cuRC1LHuQHLBxERkZVpoFLiuR6hAIBP/0iBtqpa4kS1sXwQERFZoZGdg+GjdsDFwnKsPpAhdZxaWD6IiIiskIOdAi88EAYA+GJnKkorqiRO9D8sH0RERFZqWPtABLk7Ibe4Aiv2nZM6jh7LBxERkZWyU8gx9cGmAIDFu9JQWFYpcaIaLB9ERERWrH8rfzT1agBNeRWW7TkrdRwALB9ERERWTSGX4aXe4QCAb/5OR16xVuJELB9ERERWr0+UD6L9XVFSUY1Fu9KkjsPyQUREZO1kMhmm94kAAKzcfx6XCsskzcPyQUREZAO6NfVEx8buqKjS4bM/UyXNwvJBRERkA64/+vHjoUxk5pdKlkUp2W8mIiIik+oY4o4RnYLQMcQdfm6OkuVg+SAiIrIh7w6MljoCv3YhIiIi02L5ICIiIpNi+SAiIiKTYvkgIiIik2L5ICIiIpNi+SAiIiKTYvkgIiIik2L5ICIiIpNi+SAiIiKTYvkgIiIikzKofMyZMwcymazWw8fH56Zjx48fD5lMhoULF9ZHTiIiIrISBt/bJSoqCjt27NA/VygUN4zZsGEDDhw4AD8/v3tLR0RERFbH4PKhVCpvebQDALKysjBp0iRs3boV/fr1u6dwREREZH0MLh8pKSnw8/ODSqVCp06dMHfuXDRp0gQAoNPpMGrUKMyYMQNRUVF12p9Wq4VWq9U/LywsBABoNBpDoxEREZFE/vncFkLccaxB5aNTp05YuXIlwsPDcfnyZbzzzjuIiYnB6dOn4eHhgfnz50OpVGLy5Ml13ue8efPw5ptv3rA9MDDQkGhERERkBoqKiuDq6nrbMTJRl4pyCyUlJQgNDcXLL7+M7t27o1+/fjh69Kj+XI/GjRtjypQpmDJlyi338e8jHzqdDvn5+fDw8IBMJrvbaDel0WgQGBiIzMxMqNXqet03/Q/n2TQ4z6bDuTYNzrNpGGuehRAoKiqCn58f5PLbX89i8Ncu13N2dkZ0dDRSUlIgl8uRk5ODoKAg/evV1dV46aWXsHDhQpw7d+6m+1CpVFCpVLW2ubm53UusO1Kr1fyDbQKcZ9PgPJsO59o0OM+mYYx5vtMRj3/cU/nQarVITExE165dMWrUKPTq1avW63369MGoUaMwZsyYe/k1REREZEUMKh/Tp0/Ho48+iqCgIOTk5OCdd96BRqPB6NGj4eHhAQ8Pj1rj7ezs4OPjg4iIiHoNTURERJbLoPJx4cIFDB8+HLm5uWjUqBE6d+6M/fv3Izg42Fj56pVKpcIbb7xxw9c8VL84z6bBeTYdzrVpcJ5Nwxzm+Z5OOCUiIiIyFO/tQkRERCbF8kFEREQmxfJBREREJsXyQURERCZl0eXjyy+/REhICBwcHNCuXTvs2bPnlmP//vtvxMbGwsPDA46OjmjWrBk+/vjjG8b98ssviIyMhEqlQmRkJNavX2/Mt2Ax6nuuT58+jcGDB6Nx48aQyWRYuHChkd+BZajveV66dCm6du2Khg0bomHDhujVqxcOHjxo7Ldh9up7ntetW4f27dvDzc0Nzs7OaN26NVatWmXst2H2jPF39D/Wrl0LmUyGxx57zAjJLUt9z/OKFSsgk8lueJSXl9dfaGGh1q5dK+zs7MTSpUtFQkKCePHFF4Wzs7M4f/78TccfPXpUrF69WsTHx4v09HSxatUq4eTkJJYsWaIfs2/fPqFQKMTcuXNFYmKimDt3rlAqlWL//v2meltmyRhzffDgQTF9+nSxZs0a4ePjIz7++GMTvRvzZYx5fvLJJ8UXX3whjh07JhITE8WYMWOEq6uruHDhgqneltkxxjzv3LlTrFu3TiQkJIjU1FSxcOFCoVAoxJYtW0z1tsyOMeb5H+fOnRP+/v6ia9euYsCAAUZ+J+bNGPO8fPlyoVarxaVLl2o96pPFlo+OHTuKCRMm1NrWrFkzMXPmzDrvY+DAgWLkyJH658OGDRN9+/atNaZPnz7iiSeeuLewFs4Yc3294OBglg9h/HkWQoiqqirh4uIivv3227vOaelMMc9CCNGmTRvx2muv3VVGa2Csea6qqhKxsbFi2bJlYvTo0TZfPowxz8uXLxeurq71FfGmLPJrl4qKChw5cgS9e/eutb13797Yt29fnfZx7Ngx7Nu3D927d9dvi4uLu2Gfffr0qfM+rZGx5ppqM9U8l5aWorKyEu7u7veU11KZYp6FEPjjjz+QlJSEbt263XNmS2TMeX7rrbfQqFEjjB07tt7yWipjznNxcTGCg4MREBCARx55BMeOHau33MA93ttFKrm5uaiuroa3t3et7d7e3sjOzr7tzwYEBODKlSuoqqrCnDlz8Mwzz+hfy87Ovqt9WjNjzTXVZqp5njlzJvz9/W+4D5OtMOY8FxYWwt/fH1qtFgqFAl9++SUefPDBen8PlsBY87x37158/fXXOH78uDFiWxxjzXOzZs2wYsUKREdHQ6PR4JNPPkFsbCxOnDiBpk2b1kt2iywf/5DJZLWeCyFu2PZve/bsQXFxMfbv34+ZM2ciLCwMw4cPv6d92gJjzDXdyJjz/P7772PNmjXYtWsXHBwc6jW3pTHGPLu4uOD48eMoLi7GH3/8gWnTpqFJkybo0aOHMd6CRajPeS4qKsLIkSOxdOlSeHp6GjO2xanvP8+dO3dG586d9WNjY2PRtm1bfPbZZ/j000/rJbNFlg9PT08oFIobml1OTs4NDfDfQkJCAADR0dG4fPky5syZo59wHx+fu9qnNTPWXFNtxp7nBQsWYO7cudixYwdatmxZv+EtiDHnWS6XIywsDADQunVrJCYmYt68eTZZPowxz2lpaTh37hweffRR/VidTgcAUCqVSEpKQmhoaD2/E/Nmqr+f5XI5OnTogJSUlPoJDgu91Nbe3h7t2rXD9u3ba23fvn07YmJi6rwfIQS0Wq3++X333XfDPrdt22bQPq2NseaaajPmPH/wwQd4++23sWXLFrRv375e8loqU/55tuU/88aY52bNmuHUqVM4fvy4/tG/f3/cf//9OH78OAIDA+v1PVgCU/15FkLg+PHj8PX1veusN9upRfrn8qKvv/5aJCQkiClTpghnZ2dx7tw5IYQQM2fOFKNGjdKP//zzz8XGjRtFcnKySE5OFt98841Qq9Vi9uzZ+jF79+4VCoVCvPfeeyIxMVG89957vNRWGGeutVqtOHbsmDh27Jjw9fUV06dPF8eOHRMpKSkmf3/mwhjzPH/+fGFvby9+/vnnWpfMFRUVmfz9mQtjzPPcuXPFtm3bRFpamkhMTBQffvihUCqVYunSpSZ/f+bCGPP8b7zaxTjzPGfOHLFlyxaRlpYmjh07JsaMGSOUSqU4cOBAveW22PIhhBBffPGFCA4OFvb29qJt27Zi9+7d+tdGjx4tunfvrn/+6aefiqioKOHk5CTUarVo06aN+PLLL0V1dXWtff70008iIiJC2NnZiWbNmolffvnFVG/HrNX3XKenpwsANzyu348tqu95Dg4Ovuk8v/HGGyZ8V+anvud59uzZIiwsTDg4OIiGDRuK++67T6xdu9aUb8ksGePv6OuxfNSo73meMmWKCAoKEvb29qJRo0aid+/eYt++ffWaWSaEEPV3HIWIiIjo9izynA8iIiKyXCwfREREZFIsH0RERGRSLB9ERERkUiwfREREZFIsH0RERGRSLB9ERERkUiwfREREZFIsH0RERGRSLB9ERERkUiwfREREZFIsH0RERGRS/w+nUasWxRk8jQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "pos = [0.3156, 0.049,0.6727,0.831,0.9645,-1,0.0] \n",
    "\n",
    "omega_m = np.linspace(0.3,0.35,20)\n",
    "ll = np.zeros_like(omega_m)\n",
    "for i in range(len(omega_m)):\n",
    "    pos[0] = omega_m[i]\n",
    "    ll[i] = likelihood(pos)\n",
    "\n",
    "plt.plot(omega_m,ll)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "levinpower_env",
   "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.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}