{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "edb2cf2a",
   "metadata": {},
   "source": [
    "## Speed of Sound EOS Inference pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bc633d3",
   "metadata": {},
   "source": [
    "This is an example notebook about how to use our tools to analysis a observation constraint on speed of sound star equation of state.\n",
    "\n",
    "Here in this notebook, we are using a speed of sound EOS.\n",
    "\n",
    "The following is the package that we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f6f20856",
   "metadata": {},
   "outputs": [],
   "source": [
    "import InferenceWorkflow.BayesianSampler as sampler\n",
    "import InferenceWorkflow.Likelihood as likelihood\n",
    "import InferenceWorkflow.prior as prior\n",
    "import numpy as np\n",
    "from EOSgenerators.SpeedofSound_EOS import compute_EOS\n",
    "from TOVsolver.unit import g_cm_3, dyn_cm_2, km, Msun\n",
    "import EOSgenerators.crust_EOS as crust\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edfa0f72",
   "metadata": {},
   "source": [
    "### Set up prior"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38136542",
   "metadata": {},
   "source": [
    "Next step, we need to set up the prior, first use parameters array to specify the variable name, this process should consistent with what you need to call them.\n",
    "\n",
    "Define a prior transform function to define prior. Cube are set of random number from 0 to 1. This prior setting is standard set-up of UltraNest package, since we are using UltraNest to do nest-sampling.\n",
    "\n",
    "We provided two options call from prior:\"normal_Prior\" and \"flat_prior\".\n",
    " \n",
    "Then the Parameters prior should all set.\n",
    "\n",
    "We note that since we are doing Equation of state Inference from mass-radius data of neutron star measurement. The center density of the star should be also sampled. Otherwise will be a partially-defined prior, did not span all parameters space, and proved to be different with full-scope inference.\n",
    "\n",
    "This request as randomly generate a density from a EoS range, however, this process is not that trivial, since we need to determine the upper limit of the central density of compact star --- different equation of state will predict different upper bound, so here we need to use the prior-setting EoS parameters computing the EOS by\n",
    "\n",
    "EOSgenerators.SpeedofSound_EOS.compute_EOS\n",
    "Compute out EOS, put into\n",
    "\n",
    "main.OutputMR\n",
    "\n",
    "find out Mass Radius of this equation of state, then find out the last stable point of this equation of state.(first mass points that give the direvative to be negative)\n",
    "\n",
    "and find out that index by len() function, then reset this max_d to be upper limit of this density range.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2b87dd3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = [\n",
    "    \"a1\",\n",
    "    \"a2\",\n",
    "    \"a3\",\n",
    "    \"a4\",\n",
    "    \"a5\",\n",
    "    \"d1\",\n",
    "]\n",
    "\n",
    "\n",
    "def prior_transform(cube6):\n",
    "    a1 = prior.flat_prior(0.1, 1.5, cube6[0])\n",
    "    a2 = prior.flat_prior(1.5, 12, cube6[1])\n",
    "    a3 = prior.flat_prior(0.05 * a2, 2 * a2, cube6[2])\n",
    "    a4 = prior.flat_prior(1.5, 37, cube6[3])\n",
    "    a5 = prior.flat_prior(0.1, 1, cube6[4])\n",
    "    d1 = prior.flat_prior(14.7, 15.9, cube6[5])\n",
    "    return (a1, a2, a3, a4, a5, d1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f8e4822",
   "metadata": {},
   "source": [
    "In the upper part, we define a flat (uniform) prior for the parameters in the CS equation of state, due to the lack of constraints from terrestrial experiments.\n",
    "\n",
    "Note that the above code is an example of Bayesian analysis for a given mass and radius observation measurement.\n",
    "For example, if you use the NICER data for the measurements of J0030, then you should define another parameter, except the CS EOS parameters, e.g. \"d1\" for the centre density of this measurement, in the meantime add \"params[2]\" to this code.\n",
    "\n",
    "If you further consider the adjoint analysis with J0030+J0740, then you should define the other two parameters, e.g. \"d1\" and \"d2\" for the centre density of these two measurements, in the meantime add \"params[3]\" to the above code."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53ab6d43",
   "metadata": {},
   "source": [
    "### Set up likehood"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bfd12ff",
   "metadata": {},
   "source": [
    "We need to set up a likelihood, Using standard definition way of UltraNest, that is below.\n",
    "\n",
    "Here the likelihood is generated from a simulated mass radius measurement, which is  𝑀=1.4𝑀⊙\n",
    "  and  𝑅=13\n",
    "  km, With a 5% Mass radius measurement uncertainty, \n",
    "  \n",
    " so here\n",
    " \n",
    "      likelihood.MRlikihood_Gaussian\n",
    "      \n",
    "function will be use for our likelihood, please check [likelihood.MRlikihood_Gaussian](https://github.com/ChunHuangPhy/CompactOject/blob/main/InferenceWorkflow/Likelihood.py) to see the original code, and more choice of likelihood. eg:\n",
    "\n",
    "1.If we have some real mass-radius measurements, say PSR J0030 or PSR J0740, come from NICER, a KDE kernel could be trained to feed into\n",
    "\n",
    "likelihood.MRlikihood_kernel(eps_total,pres_total,x,d1)\n",
    "\n",
    "set the KDE kernel as a input for this function\n",
    "\n",
    "2.If we gain measurement from radio-timing, say only measure the neutron star mass, then\n",
    "\n",
    "likelihood.Masslikihood_Gaussian(eps_total,pres_total,x,d1)\n",
    "\n",
    "Which will give the likelihood from single mass measurement, x is the parameters of that measurement, you should specify where this measurement mass is located and what is the sigma width of this mass measurement.\n",
    "\n",
    "3.If we have nuclear measurements, and want to constrain this RMF model by nuclear properties like K(The Incompressibility of nuclear matter),J ( the symmetry energy at saturation density) and L( the slope of symmetry energy at saturation density). You can choose:\n",
    "\n",
    "likelihood.Kliklihood(theta,K_low,K_up)\n",
    "likelihood.Jliklihood(theta,K_low,K_up)\n",
    "likelihood.Lliklihood(theta,K_low,K_up)\n",
    "\n",
    "We are defaulting a hard-cut flat constrain, so if you don't like this default hard cut, also could define the likelihood by youself with similiar style.\n",
    "\n",
    "4.If we have a Tidal measurements from Gravitational wave detector, we can use it to do constraint:\n",
    "\n",
    "likelihood.TidalLikihood_kernel(eps_total,pres_total,x,d1)\n",
    "\n",
    "Where x is sampled distribution from real measurements, the standard is\n",
    "\n",
    "kernel, chrip = x,\n",
    "\n",
    "where the kernel is a whole set sampling from GW event, that is [chrip mass, M2/M1, tidal of M1, tidal of M2] four quantities. Chrip is the single smapling that comes only the chrip mass sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0d3a8359",
   "metadata": {},
   "outputs": [],
   "source": [
    "Tolos_crust_out = np.loadtxt(\"Test_Case/Tolos_crust_out.txt\")\n",
    "eps_crust_T_out = Tolos_crust_out[:, 3] * g_cm_3\n",
    "pres_crust_T_out = Tolos_crust_out[:, 4] * dyn_cm_2\n",
    "eps_com, pres_com = crust.PolyInterpolate(eps_crust_T_out, pres_crust_T_out)\n",
    "x_last = eps_com[-1]\n",
    "y_last = pres_com[-1]\n",
    "dydx_last = (pres_com[-1] - pres_com[-2]) / (eps_com[-1] - eps_com[-2])\n",
    "CS_EOS = compute_EOS(x_last, y_last, dydx_last)\n",
    "core_e_min = x_last\n",
    "core_e_max = 2e16 * g_cm_3\n",
    "n_core_e = 1000\n",
    "core_e = np.geomspace(core_e_min, core_e_max, n_core_e)\n",
    "full_e = np.concatenate((eps_com, core_e[1:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e6a0c9c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIcAAANrCAYAAADcUPjsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAB7CAAAewgFu0HU+AAC0XklEQVR4nOzdd3hUVf7H8c+khyTU0BNqCIQaWgRpggIiRZRqQxTr2kCUFVTAuoqKoLuuiwio/BQUEIEgIiJVgYCkUAIEhSS0EEpIb3N/f7hkHWbombkp79fz+Dic77k333FZyHxy7jkWwzAMAQAAAAAAoFxyM7sBAAAAAAAAmIdwCAAAAAAAoBwjHAIAAAAAACjHCIcAAAAAAADKMcIhAAAAAACAcoxwCAAAAAAAoBwjHAIAAAAAACjHCIcAAAAAAADKMcIhAAAAAACAcoxwCAAAAAAAoBwjHAIAAAAAACjHCIcAAAAAAADKMcIhAAAAAACAcoxwCAAAAAAAoBwjHAIAAAAAACjHCIcAAAAAAADKMQ+zG0DZkJOTo7i4OElS9erV5eHBby0AAAAAAIpbQUGBTp48KUlq1aqVfHx8rvuefIJHsYiLi1NERITZbQAAAAAAUG5s27ZNHTt2vO778FgZAAAAAABAOcbKIRSL6tWrF73etm2bateubWI3AAAAAACUTceOHSt6cuevn8WvB+EQisVf9xiqXbu2goKCTOwGAAAAAICyr7j2++WxMgAAAAAAgHKMcAgAAAAAAKAcIxwCAAAAAAAoxwiHAAAAAAAAyjHCIQAAAAAAgHKMcAgAAAAAAKAcIxwCAAAAAAAoxwiHAAAAAAAAyjHCIQAAAAAAgHKMcAgAAAAAAKAcIxwCAAAAAAAoxwiHAAAAAAAAyjHCIQAAAAAAgHKMcAgAAAAAAKAcIxwCAAAAAAAoxwiHAAAAAAAAyjHCIQAAAAAAgHKMcAgAAAAAAKAcIxwCAAAAAAAoxwiHAAAAAAAAyjHCIQAAAAAAgHKMcAgAAAAAAKAcIxwCAAAAAAAoxwiHAAAAAAAAyjHCIQAAAAAAgHKMcAgAAAAAAKAcIxwCAAAAAAAoxwiHAAAAAAAAyjHCIQAAAAAAUC4ZhqGFUYn6v62HzW7FVB5mNwAAAAAAAOBqpzPz9MLiWK3ec0JeHm7q2KCqQmsGmN2WKVg5BAAAAAAAypUN+0+q74wNWr3nhCQpr8Cqp7/aqdyCQpM7MwfhEAAAAAAAKBfyC636x/d7NWrONp1Mz7WpxR9P1097U0zqzFw8VgYAAAAAAMq8pNNZeuqrnYpOOmtXq+rnpbfubKU+LWq5vrESgHAIAAAAAACUad/HHdOExbFKzymwq93UtLqmDW2tGgE+JnRWMhAOAQAAAACAMim3oFCvrdij+VsS7Wqe7hZN7BemB7o0kMViMaG7koNwCAAAAAAAlDnH0rL1+PzfHD5G1qBaBX14Vzu1Cqrk+sZKIMIhAAAAAABQpmz9/ZSe+PI3pWbk2dVuD6+j1we3VICPpwmdlUyEQwAAAAAAoEwwDEPzfjmkNyL3qsBq2NR8PN306qCWGtYhqNw/RnYhwiEAAAAAAFDq5RVYNXFJnBb/lmxXq1e1gv5zX3uF1a5oQmclH+EQAAAAAAAo1c5m5enRL3Zo6x+n7Wo9Qqtr5shwVa7gZUJnpQPhEAAAAAAAKLUOn8rUA3Oj9Htqpl3tqV4hGntLqNzdeIzsUgiHAAAAAABAqbTj8Gk9/PkOnc603Xja19NdM0aGq2+LWiZ1Vrq4md0ALu3s2bN6+umn1blzZ9WqVUve3t6qW7euevXqpcWLF8swjIte++2336p3796qVq2afH191bBhQ911111KSkpy4TsAAAAAAKD4rduXontmb7ULhmoEeOubxzoTDF0FVg6VcKmpqZozZ446deqkwYMHq2rVqkpJSdHy5cs1dOhQPfzww5o1a5bNNYZh6LHHHtOsWbPUuHFjjRw5UgEBATp69KjWr1+vw4cPKzg42KR3BAAAAADA9VkRe1TjFkYrv9B2wUSzWgGaM7qj6lT2Namz0olwqIRr2LChzp49Kw8P2/+p0tPT1alTJ33yySd65pln1KJFi6Lahx9+qFmzZumJJ57QzJkz5e7ubnNtQUGBS3oHAAAAAKC4fbUtUZO+jdOFD9Lc1LS6/nl3O/l7E3VcLR4rK+Hc3d3tgiFJCggIUN++fSVJCQkJRePZ2dl65ZVX1KhRI82YMcMuGJLk8H4AAAAAAJR0n/96SBOX2AdDg8Pr6JNRHQiGrhHh0CWkpKRoxYoVmjx5svr166fAwEBZLBZZLBaNHj36qu6VmJio5557TmFhYfLz81PVqlUVERGhd999V1lZWVfdW05OjtauXSuLxaLmzZsXjf/44486ffq0Bg8erMLCQi1ZskRvvfWWPv74Y5sQCQAAAACA0uTLrYma/N1uu/FRnetr+vBweboTcVwrIrVLqFmzZrHcJzIyUvfcc4/S0tKKxrKyshQVFaWoqCjNnj1bK1euVKNGjS56j7Nnz2rGjBmyWq1KSUnRypUrlZSUpClTpqhJkyZF87Zv3y7pz9VBbdq00b59+4pqbm5uGjdunN59991ieV8AAAAAALjCN9uTNOnbOLvxJ3uGaHyfUFksHFV/PYjVrlBwcLD69Olz1dfFxMRo+PDhSktLk7+/v9544w398ssv+umnn/Twww9Lkvbt26f+/fsrIyPjovc5e/asXnnlFb322mv6z3/+o+PHj+udd97RlClTbOalpKRIkt577z1VrFhR27ZtU3p6ujZs2KDQ0FC99957+ve//33V7wMAAAAAADMsjzmqCYtj7caf6xOq5/o2JRgqBoRDlzB58mQtX75cx48fV2Jiov7zn/9c9T3Gjh2rrKwseXh4aPXq1Zo0aZI6d+6sXr16adasWZo2bZokKT4+XtOnT7/ofRo0aCDDMFRQUKA//vhDr776ql588UUNGTLEZoNpq9UqSfLy8tLSpUvVsWNH+fv7q1u3blq0aJHc3Nz03nvvXfX7AAAAAADA1X49eErjv46x22Po6Zub6MleTRxfhKtGOHQJr7zyigYMGHDNj5dFRUVp3bp1kqQxY8aoc+fOdnPGjx+vsLAwSdKMGTOUn59/yXu6u7urQYMGeuGFF/T666/r22+/1SeffFJUr1SpkiSpQ4cOqlOnjs21LVq0UKNGjXTw4EGdPXv2mt4TAAAAAACusP9Euh75YrvyCq0244/1aKxxtxAMFSfCISdaunRp0esHHnjA4Rw3NzeNGjVKknTmzJmiMOlKnH/M7a/XNG3aVJJUuXJlh9ecH8/Ozr7irwMAAAAAgCudOJej0XO2KT2nwGb8vk719fdbeZSsuBEOOdHGjRslSX5+fmrfvv1F5/Xo0aPo9aZNm674/kePHpVkezR9z549JUl79+61m5+fn6+EhAT5+fmpevXqV/x1AAAAAABwldyCQj02f4eOpuXYjPdpXlNTB7UgGHICwiEnOh/QhISE2AQ4F2rWrJndNedFR0fbnHJ23unTpzVp0iRJUr9+/YrGGzdurD59+ighIUGzZ8+2ueatt97S2bNndccdd1yyHwAAAAAAzDJ12R7tTDxrM9auXmV9cFdbubsRDDkDCYGT5OTkKDU1VZIUFBR0yblVqlSRn5+fMjMzlZSUZFObN2+eZs+erZ49e6p+/fry8/PT4cOHFRkZqYyMDA0ZMkR33323zTUfffSRbrzxRj388MNaunSpmjVrpp07d2rt2rWqX7++3nnnnat+P8nJyZesHzt27KrvCQAAAADAX325NVFfbUu0GQuu6qvZ93eUj6e7SV2VfYRDTpKenl702t/f/7Lzz4dDFx5nP3ToUKWlpWnLli3asGGDsrKyVLVqVXXt2lWjRo3SyJEj7ZbUNW7cWNu3b9fkyZO1atUqrV69WrVq1dITTzyhyZMnq0aNGlf9foKDg6/6GgAAAAAArtSuI2maumy3zZivp7tm3ddBVf28TOqqfCAccpKcnP89G+nldfnfxN7e3pLsN4ru2rWrunbtetVfPzg4WHPnzr3q6wAAAAAAcLXsvEI9s2Cn3clkbw9trbDaFU3qqvwgHHISHx+fotd5eXmXnZ+bmytJ8vX1dVpP1+PCx90udOzYMUVERLioGwAAAABAWfJa5B4dPJlpM/ZQ14Ya1KaOSR2VL4RDThIQEFD0+sJHxRzJzPzz/wRX8giaGS63bxIAAAAAANdi9e7j+nKr7T5DrepW0oRbm13kChQ3TitzEh8fHwUGBkq6/GbOZ86cKQqH2NsHAAAAAFBenMvJ14tLd9mM+Xq6a+bIcHl5EFm4Cv+lnSgsLEySlJCQoIKCgovOi4+Pt7sGAAAAAICy7p1V+3QyPddmbOqg5mpUvWQ+VVNWEQ450fmNpDMzM7Vjx46Lzlu/fn3R6y5duji9LwAAAAAAzPZb4hnN33rYZuymptU1vANP1Lga4ZATDR48uOj1xU4Os1qt+vzzzyVJlStXVs+ePV3RGgAAAAAApskvtGrSkjgZxv/GfDzd9NrtLWWxWMxrrJwiHHKiiIgIdevWTZL06aef6tdff7Wb895772nv3r2SpGeeeUaenp4u7REAAAAAAFf7dNMfij+ebjM27pZQBVetYFJH5RunlV3Cpk2blJCQUPTr1NTUotcJCQmaN2+ezfzRo0fb3WPmzJnq0qWLsrOz1adPH02aNEk9e/ZUdna2FixYoFmzZkmSQkNDNX78eKe8DwAAAAAASoqk01masWa/zVizWgF6sGtDkzqCxTD+uogLfzV69Gh99tlnVzz/Yv8ply9frnvvvVfnzp1zWA8NDVVkZKRCQkKuqc+SIDk5ueiktaSkJAUFBZncEQAAAACgpDEMQ6PnRmn9/pNFYxaLtOTxG9W2XhUTOys9nPH5m8fKXGDgwIGKjY3VuHHjFBoaqgoVKqhy5crq0KGD3n77be3cubNUB0MAAAAAAFyJyLhjNsGQJN3XqT7BkMlYOYRiwcohAAAAAMClpGXn65bp622Orq8R4K0143uoog/7714pVg4BAAAAAIBS6a3v422CIUmaOqgFwVAJQDgEAAAAAACcasvvp/TVtkSbsV7Naqhfy1omdYS/IhwCAAAAAABOk5NfqIlL4mzGKni569XbW8hisZjUFf6KcAgAAAAAADjNK8v36I/UTJux5/s2VVCVCiZ1hAsRDgEAAAAAAKf4ZnuS3eNk4cGVNapzA3MagkOEQwAAAAAAoNgtizmqFxw8TvbusDZyd+NxspLEw+wGAAAAAACAc+UVWHXwZIaOp+UoPbdAeQVWSZJhGDIkyZAMGTIMyZD+++///Vr/nWcY/7vGMP68t/Hf+/x3mgoNQ3FH0hQZe8yuj9cHt1RIDX+nv19cHcIhAAAAAADKoISUdH0fd1xr4lO0+0iaCqyGqf080r2R7mwXZGoPcIxwCAAAAACAMsJqNbR6zwnN++UPbfn9tNntFLmvU31N7NfM7DZwEYRDAAAAAACUAZsTUvXW9/GKO5JmditFArw99MJtzXR3RD2OrS/BCIcAAAAAACjFUtJzNOW73fp+1/HLzq1Z0VsVfTzl5eEmi0WyyPLff0uyWGT581///fdff/3n4F9/bbH8+Vr636/13+t8PNwU0bCq7mwXpKp+Xk553yg+hEMAAAAAAJRChmFoafQRTV22R2nZ+Q7nBPp7aUDrOro5rIbaBFdWRR9PF3eJ0oBwCAAAAACAUiY7r1AvLd2lxb8lO6w3qu6np3qF6LZWteXt4e7i7lDaEA4BAAAAAFCKHDyZob/N/037TqTb1ar6een5vk01rH2QPNzdTOgOpRHhEAAAAAAApcQvCal6dP4OpecU2NUGtqmjqQObq5q/twmdoTQjHAIAAAAAoBRYvCNZLyyJVX6hYTPu7+2ht4a00oDWdUzqDKUd4RAAAAAAACWYYRj6cG2Cpv+4367WrFaAPrqnnRpV9zehM5QVhEMAAAAAAJRQhmHore/j9Z8Nv9vVbm1RS++PCJevFxtO4/oQDgEAAAAAUAJZrYamLt+tz389bFd7uFtDTewXJjc3iwmdoawhHAIAAAAAoIQptBqatCROC7cn2YxbLNIrg1poVOcG5jSGMolwCAAAAACAEsQwDL201D4YcrNI7w5rozvbBZnUGcoqwiEAAAAAAEoIwzD0RuRefbXNNhjycLNo5si26t+6tkmdoSwjHAIAAAAAoISY+dMBzd70h82Yp7tFH93TXr2b1zSpK5R1bmY3AAAAAAAApNkbf9eMNQdsxtzdLPrwrnYEQ3AqwiEAAAAAAEy25LdkvR6512bMYpHeHdZat7asZVJXKC8IhwAAAAAAMNHmhFRNWBRrN/7a7S11R1s2n4bzEQ4BAAAAAGCS+OPn9NgXO1RgNWzGX+jXTPd2qm9SVyhvCIcAAAAAADDB8bQcPTA3Sum5BTbjY7o21GM9GpvUFcojwiEAAAAAAFwsPSdfo+du07G0HJvxfi1r6cXbwkzqCuUV4RAAAAAAAC5UaDX09Fc7FX883Wa8ff0qen9EuNzcLCZ1hvKKcAgAAAAAABeatipeP+87aTPWMNBPs0d1kI+nu0ldoTwjHAIAAAAAwEWW/Jas/2z43Wasqp+X5j3QUVX8vEzqCuUd4RAAAAAAAC6wM/GMXlgSZzPm6W7Rx/e2V/1qfiZ1BRAOAQAAAADgdMfTcvToFzuUV2C1GX/t9paKaFjVpK6APxEOAQAAAADgRDn5hXrki+1KSc+1GR99YwONjKhnUlfA/xAOAQAAAADgRJO/26XY5DSbsS4h1fRSf46sR8lAOAQAAAAAgJMsjErU19uTbcbqV6ugf93dTh7ufCRHycDvRAAAAAAAnCAuOU0vf7fbZszPy12zR3VQ5QqcTIaSg3AIAAAAAIBidjYrT4//n/0G1G8Pba0mNQNM6gpwjHAIAAAAAIBiZLUaGrcwWslnsm3GH+zSUANa1zGpK+DiCIcAAAAAAChG//o5QT/vO2kz1qF+FU28rZlJHQGXRjgEAAAAAEAx2XjgpKav2W8zFujvpX/e3U6ebECNEorfmQAAAAAAFIOU9ByNWxgtw/jfmJtF+uCutqpVyce8xoDLIBwCAAAAAOA6Wa2Gxn8do9SMPJvx5/s2042NA03qCrgyhEMAAAAAAFyn/2z4XRsPpNqM9WxaXY92b2RSR8CVIxwCAAAAAOA67Dh8Ru+u3mczVrOit94d1kZubhaTugKuHOEQAAAAAADXKC07X09/tVOF1v9tNGSxSDNGtFU1f28TOwOuHOEQAAAAAADXwDAMTVoSpyNns23Gn+rVRJ0bVzOpK+DqEQ4BAAAAAHANvtqWpMi4YzZjEQ2q6uleISZ1BFwbwiEAAAAAAK7SwZMZenXFbpuxyhU8NWNkuDzc+aiN0oXfsQAAAAAAXIX8QqvGLYxWTr7VZvydoW1Up7KvSV0B145wCAAAAACAq/Dh2gTFJqfZjI3qXF+9m9c0qSPg+hAOAQAAAABwhX5LPKN//ZxgM9a4up8m9gszqSPg+hEOAQAAAABwBTJzC/TswmibY+s93Cx6f0S4fL3cTewMuD6EQwAAAAAAXIE3Vu7VoVNZNmPP3NxErYMqm9MQUEwIhwAAAAAAuIy18Sf05dZEm7F29Srr8Zsam9QRUHwIhwAAAAAAuITTmXmasCjOZqyCl7umD+fYepQN/C4GAAAAAOASpizbrdSMXJuxlwc0V4NAP5M6AooX4RAAAAAAABexatcxLY85ajN2S1gNjewYbFJHQPEjHAIAAAAAwIEzmXl6aekum7HKFTz15p2tZLFYTOoKKH6EQwAAAAAAOPDK8t1KzcizGZs6sIVqBPiY1BHgHIRDAAAAAABc4Mc9J7Q02v5xstvD65jUEeA8hEMAAAAAAPzF2aw8TfrW9nSyij4eeuMOHidD2UQ4BAAAAADAX7y6Yo9OptueTjZ5YAvVrMjjZCibCIcAAAAAAPivtfEntOS3IzZjNzWtriHt6prUEeB8hEMAAAAAAEjKzC3QS9/ank4W4O2hf3A6Gco4wiEAAAAAACTNWLNfR9NybMZeGhCm2pV8TeoIcA3CIQAAAABAubf7aJrmbD5kM9a5UTUN7xBsTkOACxEOAQAAAADKtUKroUnf7lKh1Sga83J30+t3tORxMpQLhEMAAAAAgHLty62HFZN01mbs8Zsaq3F1f3MaAlyMcAgAAAAAUG6dOJejaav22Yw1CvTT4zc1NqkjwPUIhwAAAAAA5darK/YoPbfAZuz1O1rKx9PdpI4A1yMcAgAAAACUSz/vS1Fk7DGbsTvb1dWNjQNN6ggwB+EQAAAAAKDcyc4r1MtLd9mMVa7gqRdvCzOpI8A8hEMAAAAAgHLng7UHlHwm22ZsYr9mqubvbVJHgHkIhwAAAAAA5cq+4+n6ZMPvNmMRDapqWPtgkzoCzEU4BAAAAAAoN6xWQ5O+jVOB1Sga83S36I07WsrNzWJiZ4B5CIcAAAAAAOXGwu1J2nH4jM3YI90bqUnNAJM6AsxHOAQAAAAAKBdOpufqHyv32ozVq1pBT/VqYlJHQMlAOAQAAAAAKBfeiNyjczkFNmOvD24pH093kzoCSgbCIQAAAABAmbfpQKqWRh+1GRvUpo66h1Y3qSOg5CAcAgAAAACUaTn5hXppaZzNWICPh14aEGZSR0DJQjgEAAAAACjTPvo5QYdOZdmM/f3WZqoR4GNSR0DJQjgEAAAAACizElLS9e/1B23G2tarrLsj6pnUEVDyEA4BAAAAAMokq9XQpCW7lF9oFI25u1n05h2t5OZmMbEzoGQhHAIAAAAAlElfb0/StkOnbcYe6tpQYbUrmtQRUDIRDgEAAAAAypyT6bl6c+Vem7GgKr565pYmJnUElFyEQwAAAACAMue1FXt0LqfAZuz1wS1VwcvDpI6AkotwqBQ4e/asnn76aXXu3Fm1atWSt7e36tatq169emnx4sUyDOO65gMAAABAWbJuX4qWxRy1GRvUpo5ualrDpI6Ako1wqBRITU3VnDlz5Ofnp8GDB2v8+PHq16+fdu/eraFDh+rRRx+9rvkAAAAAUFZk5hbopaW7bMYq+njo5QHNTeoIKPksBstISrzCwkIZhiEPD9vlj+np6erUqZP27NmjXbt2qUWLFtc0vzgkJycrODhYkpSUlKSgoKBiuzcAAAAAXKmXlsZp/pZEm7G37mylkRxdjzLCGZ+/WTlUCri7u9sFPZIUEBCgvn37SpISEhKueT4AAAAAlAUb9p+0C4YiGlbV8A7BJnUElA6EQ5eRkpKiFStWaPLkyerXr58CAwNlsVhksVg0evToq7pXYmKinnvuOYWFhcnPz09Vq1ZVRESE3n33XWVlZV11bzk5OVq7dq0sFouaN7/8EsmrnQ8AAAAApUVadr7+vjjWZszX013vDG0tNzeLSV0BpQPbtF9GzZo1i+U+kZGRuueee5SWllY0lpWVpaioKEVFRWn27NlauXKlGjVqdNF7nD17VjNmzJDValVKSopWrlyppKQkTZkyRU2a2B/HeLXzAQAAAKA0MgxDU77bpWNpOTbjk25rpvrV/EzqCig9CIeuQnBwsMLCwrR69eqrui4mJkbDhw9XVlaW/P39NXHiRPXs2VPZ2dlasGCBPvnkE+3bt0/9+/dXVFSU/P39Hd7n7NmzeuWVV4p+7enpqXfeeUfjx48vlvkAAAAAUBotiErS0mjb08m6hgTqnhvqm9QRULoQDl3G5MmT1bFjR3Xs2FE1a9bUoUOH1LBhw6u6x9ixY5WVlSUPDw+tXr1anTt3Lqr16tVLTZo00YQJExQfH6/p06dr8uTJDu/ToEEDGYahwsJCJSUlacGCBXrxxRf1yy+/6Ouvv7bbZ+hq5wMAAABAabP7aJqmLNttMxbg7aG3eZwMuGLsOXQZr7zyigYMGHDNj5dFRUVp3bp1kqQxY8bYBEPnjR8/XmFhYZKkGTNmKD8//5L3dHd3V4MGDfTCCy/o9ddf17fffqtPPvmk2OYDAAAAQGmQkp6jRz7fobwCq834tKGtVbeyr0ldAaUP4ZCTLV26tOj1Aw884HCOm5ubRo0aJUk6c+ZMUZh0Jfr06SNJV3zN1c4HAAAAgJIoK69AY+Zt15Gz2TbjD3RpoH6tapvUFVA6EQ452caNGyVJfn5+at++/UXn9ejRo+j1pk2brvj+R4/++VztlT4idrXzAQAAAKCkSc/J1+g5UYo7kmYz3rZeZU3sF2ZSV0DpRTjkZHv37pUkhYSEXDKQadasmd0150VHR9uccnbe6dOnNWnSJElSv379rnk+AAAAAJQWv5/M0F2fbNG2Q6dtxutVraBPRnWQlwcfc4GrxfIRJ8rJyVFqaqokKSgo6JJzq1SpIj8/P2VmZiopKcmmNm/ePM2ePVs9e/ZU/fr15efnp8OHDysyMlIZGRkaMmSI7r777muefyWSk5MvWT927NhV3Q8AAAAAroZhGJq/NVFvRO5RTr7tHkOVK3hq7gMdFejvbVJ3QOlGOORE6enpRa8vdjz9X50PhzIyMmzGhw4dqrS0NG3ZskUbNmxQVlaWqlatqq5du2rUqFEaOXKkLBbLNc+/EsHBwVc1HwAAAACKS8q5HD2/KFbr95+0qwX6e2n+QzeocfXLf+YC4BjhkBPl5OQUvfby8rrsfG/vP1Pu7GzbDdW6du2qrl27XvHXvdr5AAAAAFBSrYw7pknfxulslv2pzo0C/fTJ/R0IhoDrRDjkRD4+PkWv8/LyLjs/NzdXkuTrW/KOXLzwUbcLHTt2TBERES7qBgAAAEBZdy4nX1O/260lO484rI/qXF8T+4XJ18vdxZ0BZQ/hkBMFBAQUvb7wUTFHMjMzJV3ZI2iudrk9kwAAAACguPxyMFXPfR2jo2k5drUaAd6aNrS1bmpaw4TOgLKJcMiJfHx8FBgYqNTU1Mtu6HzmzJmicIj9fQAAAACURzn5hXr3h32avekPh/X+rWrr9cEtVcXv8tt2ALhyhENOFhYWpo0bNyohIUEFBQUXPc4+Pj7e5hoAAAAAKE92H03TuIXR2n/C/qmLAB8PvXZ7S90eXueqD9cBcHluZjdQ1p3fGDozM1M7duy46Lz169cXve7SpYvT+wIAAACAkqDQauijdQka/K/NDoOhzo2qadXY7hrcti7BEOAkhENONnjw4KLXc+fOdTjHarXq888/lyRVrlxZPXv2dEVrAAAAAGCqxFNZGvGfXzVt1T7lFxo2NS8PN73UP0z/99ANqlu55B3aA5QlhENOFhERoW7dukmSPv30U/366692c9577z3t3btXkvTMM8/I09PTpT0CAAAAgCsZhqGFUYnqN3ODth8+Y1dvXruiVjzVVQ91ayQ3N1YLAc7GnkOXsWnTJiUkJBT9OjU1teh1QkKC5s2bZzN/9OjRdveYOXOmunTpouzsbPXp00eTJk1Sz549lZ2drQULFmjWrFmSpNDQUI0fP94p7wMAAAAASoLUjFy9sDhOa/aesKu5WaTHejTW2FtC5eXBWgbAVSyGYRiXn1Z+jR49Wp999tkVz7/Yf87ly5fr3nvv1blz5xzWQ0NDFRkZqZCQkGvq02zJyclFp6wlJSUpKCjI5I4AAAAAlDQ/7jmhFxbH6lRmnl2tXtUKmj68jTo0qGpCZ0Dp4YzP36wccpGBAwcqNjZWM2fOVGRkpJKTk+Xl5aWQkBANGzZMTz75pCpUqGB2mwAAAABQ7DJyC/Ta8j1auD3JYX1kx2C9NKC5/L35iAqYgZVDKBasHAIAAADgSNSh03r262glnc62q1Xz89JbQ1qrd/OaJnQGlE6sHAIAAAAAlAp5BVa9v2a/Pl5/UI6WJNwSVlNvDWmlQH9v1zcHwAbhEAAAAACgWO07nq5xC6O155j9nqt+Xu6aMrCFhnUIksXCSWRASUA4BAAAAAAoFlaroTmb/9C0H/Ypr8BqV+9Qv4qmDw9XvWrstwqUJIRDAAAAAIDrduRstp77Oka//n7KrubpbtG43qF6tHtjubuxWggoaQiHAAAAAADXzDAMLY0+oslLdys9t8CuHlrTX++PCFeLOpVM6A7AlSAcAgAAAABck9OZeXppaZxWxh23q1ks0pguDfVc36by8XQ3oTsAV4pwCAAAAABw1dbGn9DfF8fpZHquXa1OJR+9O7yNbmwcaEJnAK4W4RAAAAAA4Ipl5Bbojcg9+mpbksP6nW3rasqgFqrk6+nizgBcK8IhAAAAAMAV2fbHaY3/JlpJp7PtapUreOqNwa3Uv3VtEzoDcD0IhwAAAAAAl5RbUKjpq/dr1sbfZRj29Z5Nq+vtIa1Vo6KP65sDcN0IhwAAAAAAF7X7aJqeXRijfSfS7WoVvNz18oDmGtkxWBYLR9QDpRXhEAAAAADATkGhVf/Z8LtmrNmv/EL75UIdG1TRu8PaqH41PxO6A1CcCIcAAAAAADYOpWbq2a+j9VviWbual7ubxvcJ1UPdGsndjdVCQFlAOAQAAAAAkCQZhqH/25qoNyL3Kju/0K4eVrui3h/RRs1qVTShOwDOQjgEAAAAANDxtBxNWByrDftP2tXcLNJjPRrrmVuayNvD3YTuADgT4RAAAAAAlHPLYo7q5aW7lJadb1erX62Cpg9vo/b1q5rQGQBXIBwCAAAAgHLqTGaeXv5ul1bEHnNYv7dTPU3sFyY/bz46AmUZ/w8HAAAAgHLo530p+vuiWKWk59rVagR4a9rQ1rqpaQ0TOgPgaoRDAAAAAFCOZOYW6I2Ve/Xl1kSH9YFt6ui121uocgUvF3cGwCyEQwAAAABQTmw/dFrPfh2jxNNZdrVKvp56bXBLDWpTx4TOAJiJcAgAAAAAyrjcgkLNWHNA/1l/UFbDvt4jtLqmDW2tmhV9XN8cANMRDgEAAABAGbb32DmNWxit+OPpdjVfT3e92D9M99xQTxaLxYTuAJQEhEMAAAAAUAYVWg3N2vC7pv+4T/mF9suF2tWrrOnDw9Ug0M+E7gCUJIRDAAAAAFDGHD6VqfFfx2j74TN2NU93i8b1DtWj3RvL3Y3VQgAIhwAAAACgzDAMQ19tS9LrkXuUlVdoV29WK0DTh4ereZ2KJnQHoKQiHAIAAACAMiDlXI4mLI7Vun0n7WoWi/Ro98Ya17uJvD3cTegOQElGOAQAAAAApdyK2KN6aekunc3Kt6vVq1pB7w1vo44NqprQGYDSgHAIAAAAAEqps1l5mvzdbi2LOeqwfvcN9fTibWHy8+ajH4CL408IAAAAACiFNuw/qecXxejEuVy7WvUAb00b0lo9m9UwoTMApQ3hEAAAAACUIll5BfrHynh9seWww3r/VrX1+uCWquLn5eLOAJRWhEMAAAAAUErsOHxG47+O1qFTWXa1ij4eem1wSw1qU0cWC0fUA7hyhEMAAAAAUMLlFVg186f9+ve6g7Ia9vVuTQL1ztA2qlXJx/XNASj1CIcAAAAAoATbdzxd4xZGa8+xc3Y1H083vXhbmO7tVJ/VQgCuGeEQAAAAAJRAhVZDn276Xe/+sF95hVa7ett6lTV9eLgaBvqZ0B2AsoRwCAAAAABKmMRTWXrumxhtO3TarubhZtG43qF6tHsjebi7mdAdgLKGcAgAAAAASgjDMPTltkS9EblXWXmFdvXQmv6aPjxcLetWMqE7AGUV4RAAAAAAlADH0rI1YVGsNh5ItatZLNLD3Rrp2d6h8vF0N6E7AGUZ4RAAAAAAmMgwDH2784imLNut9JwCu3pQFV+9N6yNbmhUzYTuAJQHhEMAAAAAYJLUjFy9+G2cfth9wmH97hvqadJtYfL35qMbAOfhTxgAAAAAMMGqXcc06dtdOp2ZZ1erWdFbbw9prZua1jChMwDlDeEQAAAAALhQWla+pizbpaXRRx3W72hbV1MHtlClCp4u7gxAeUU4BAAAAAAusm5fiv6+OFYnzuXa1ar5eemNO1rq1pa1TegMQHlGOAQAAAAATpaRW6A3Ivfqq22JDut9mtfUm3e2UqC/t4s7AwDCIQAAAABwqi2/n9Jz38Qo+Uy2XS3Ax0Ov3t5Cg8PrymKxmNAdABAOAQAAAIBT5OQX6p0f9mnO5j9kGPb17qHV9faQVqpdydf1zQHAXxAOAQAAAEAxi046q2e/jtbvJzPtahW83PVS/+a6KyKY1UIASgTCIQAAAAAoJnkFVn3w0wF9tC5BVgerhSIaVNW7w9qoXrUKrm8OAC6CcAgAAAAAisHeY+f07Ncx2nvsnF3Ny8NNE/o21YNdGsrNjdVCAEoWwiEAAAAAuA4FhVb9Z8PvmrFmv/IL7ZcLtQ6qpOnD2yikRoAJ3QHA5REOAQAAAMA1OngyQ+O/jlF00lm7moebRc/c3ESP39RYHu5urm8OAK4Q4RAAAAAAXCWr1dDcXw7pnR/ilZNvtas3rRmg94a3Ucu6lUzoDgCuDuEQAAAAAFyFQ6mZen5RjKIOnbGruVmkR3s01thbmsjbw92E7gDg6hEOAQAAAMAVsFoNzfvlkKZdZLVQw0A/vTusjdrXr2JCdwBw7QiHAAAAAOAyDqVmasKiWG07dNphffSNDfT3W5vJ14vVQgBKH8IhAAAAALgIq9XQZ78e0turHK8WCq7qq2lD2qhz42omdAcAxYNwCAAAAAAcOHwqU88vitW2PxyvFhrVub7+fmsz+XnzsQpA6cafYgAAAADwF1aroc9/PaS3V+1Tdn6hXT2oiq+mDW2tGxsHmtAdABQ/wiEAAAAA+K/EU1l6flGMtl5ktdB9nerrhX6sFgJQtvAnGgAAAIByz2o1NH/rYb31fbyy8i6yWmhIa90YwmohAGUP4RAAAACAci3p9J+rhbb87ni10L2d6umFfmHyZ7UQgDKKP90AAAAAlEuXWy1Ut/Kfewt1YbUQgDKOcAgAAABAuXP4VKb+vjj2oquF7r6hnibdxmohAOUDf9IBAAAAKDcKrYbmbv5D767ep5x8q129bmVfvT2ktbo2YbUQgPKDcAgAAABAubD/RLqeXxSrmKSzDuusFgJQXvGnHgAAAIAyLa/Aqn+vO6h//nxA+YWGXb1uZV+9NaSVujWpbkJ3AGA+wiEAAAAAZVZM0ln9fXGs4o+nO6zf37m+nr+1GauFAJRr/AkIAAAAoMzJyS/U+z/u1ycbf5fVfrGQGgX66e2hrdWxQVXXNwcAJQzhEAAAAIAyZevvp/T3xbE6dCrLrubuZtEj3RvpmZubyMfT3YTuAKDkIRwCAAAAUCak5+Tr7VXxmr8l0WE9rHZFvTO0tVrWreTizgCgZCMcAgAAAFDq/bT3hF5euktH03Lsal7ubnrmliZ6pHsjebq7mdAdAJRshEMAAAAASq0T53I0ddlufb/ruMN6u3qVNW1oa4XUCHBxZwBQehAOAQAAACh1rFZD/7f1sKat2qf03AK7uq+nu57v21T339hA7m4WEzoEgNKDcAgAAABAqRJ//JwmLonTzsSzDutdQwL1jztbKbhqBdc2BgClFOEQAAAAgFIhO69QM386oNkbf1eBg/Ppq/p56aX+YbqjbV1ZLKwWAoArRTgEAAAAoMTbsP+kXlq6S4mn7Y+nl6Sh7YM06bYwVfXzcnFnAFD6EQ4BAAAAKLFOpufq9cg9+i76qMN6w0A/vXFHS93YONDFnQFA2UE4BAAAAKDEKSi06v+2Jurd1fuUnmO/4bSnu0WP92isv/UMkY+nuwkdAkDZQTgEAAAAoETZcfiMXl66S3uOnXNY79igit68o5Wa1OR4egAoDoRDAAAAAEqEUxm5entVvL7enuywXtHHQxNvC9OIDsFy43h6ACg2hEMAAAAATFVoNfTltkS9sype5xw8QiZJd7Stq4m3NVONAB8XdwcAZR/hEAAAAADTRCed1ctLdynuSJrDemhNf712e0vd0KiaizsDgPKDcAgAAACAy504l6Npq/Zp8W+OHyHz83LXuN6huv/GBvJ0d3NxdwBQvhAOAQAAAHCZnPxCzd74uz5ad1BZeYUO5wxqU0cv9g9TzYo8QgYArkA4BAAAAMDpDMPQyrjjenPlXh05m+1wTuPqfnrt9pa6MSTQxd0BQPlGOAQAAADAqXYdSdOry/do26HTDut+Xu566uYmerBLQ3l58AgZALga4RAAAAAApzh6Nlvv/7hfi35LlmHY1y0WaXj7YI3vG8opZABgIsIhAAAAAMUqLStfH61P0LzNh5RbYHU4J6JhVU0e0Fwt61ZycXcAgAsRDgEAAAAoFjn5hfrsl0P6188JOpdT4HBOUBVfvXhbmG5tWUsWi8XFHQIAHOGB3hLu7Nmzevrpp9W5c2fVqlVL3t7eqlu3rnr16qXFixfLcLQ+V1JUVJRuu+02ValSRX5+foqIiNCXX37p4u4BAABQHhRaDX0dlaSe767TP76PdxgM+Xm5a8KtTbXm2R7q16o2wRAAlCCsHCrhUlNTNWfOHHXq1EmDBw9W1apVlZKSouXLl2vo0KF6+OGHNWvWLJtr1q1bp759+8rLy0sjR45UpUqVtGTJEt1zzz06dOiQJk2aZNK7AQAAQFlitRpaEXdMH/x0QAkpGQ7neLhZdPcN9fRUryaqHuDt4g4BAFfCYlxs6QlKhMLCQhmGIQ8P2xwvPT1dnTp10p49e7Rr1y61aNFCklRQUKBmzZopOTlZv/76q9q2bVs0v3Pnztq3b5/27NmjJk2aFGufycnJCg4OliQlJSUpKCioWO8PAACAksNqNbRy1zHNXHNABy4SCknSgNa19VyfpmoQ6OfC7gCgbHPG528eKyvh3N3d7YIhSQoICFDfvn0lSQkJCUXja9eu1cGDB3X33XcXBUPn57/88ssqKCjQ3Llznd84AAAAyhyr1dDKuGPqN3Ojnvxy50WDoRsbV9OyJ7von3e3IxgCgFKAx8ouISUlRdu2bdO2bdsUFRWlqKgonTp1SpJ0//33a968eVd8r8TERH3wwQeKjIxUYmKivL29FRISouHDh+tvf/ubKlSocFW95eTkaO3atbJYLGrevHnR+Lp16yRJffr0sbvm/Nj69euv6msBAACgfCsotGrlruP66OcExR9Pv+i85rUr6oV+zdStSSB7CgFAKUI4dAk1a9YslvtERkbqnnvuUVpaWtFYVlZWUeA0e/ZsrVy5Uo0aNbroPc6ePasZM2bIarUqJSVFK1euVFJSkqZMmWLziNiBAwckyeFjY1WqVFFgYGDRHAAAAOBScvIL9c32JM3a+LuSTmdfdF6zWgEae0sT9WleS25uhEIAUNoQDl2h4OBghYWFafXq1Vd1XUxMjIYPH66srCz5+/tr4sSJ6tmzp7Kzs7VgwQJ98skn2rdvn/r376+oqCj5+/s7vM/Zs2f1yiuvFP3a09NT77zzjsaPH28z73wAValSJYf3qVixopKTk6/qPQAAAKB8ScvK1xdbDmnu5kM6lZl30XlNawbomVua6NYWhEIAUJoRDl3C5MmT1bFjR3Xs2FE1a9bUoUOH1LBhw6u6x9ixY5WVlSUPDw+tXr1anTt3Lqr16tVLTZo00YQJExQfH6/p06dr8uTJDu/ToEEDGYahwsJCJSUlacGCBXrxxRf1yy+/6Ouvv3a4LxEAAABwNRJSMvT5r4e0eEeyMvMKLzqvSQ1/PXNLE93WsjahEACUAU5PFF599VVnfwlJumiocj3+ulLnWkRFRRXtATRmzBibYOi88ePHa+7cudq7d69mzJihiRMnytPT86L3dHd3V4MGDfTCCy/I3d1dEyZM0CeffKLHH39c0v9WDP31Eba/Onfu3EVXFQEAAKD8sVoN/bwvRfN+OaSNB1IvObdNUCU91qOx+rSoJXdCIQAoM5weDk2dOtUlm9E5Ixy6XkuXLi16/cADDzic4+bmplGjRmnixIk6c+aM1q1bp969e1/R/fv06aMJEyZo3bp1ReHQ+b2GDhw4oPbt29vMP3PmjFJTU3XjjTdew7sBAABAWXI6M09LfkvW578eVuLprEvO7R5aXY/1aKTOjaqx0TQAlEEuexbJMAyn3buk/gW1ceNGSZKfn59dUPNXPXr0KHq9adOmKw6Hjh49Kkk2j5T16NFD//jHP7R69WqNHDnSZv75/ZL++vUAAABQfhRaDW08cFLfbE/W6j3HlV948e/R3SxS/9Z19Gj3RmpZl5XnAFCWuSwc2rVrl82R68V1z9atWxfrPYvT3r17JUkhISGX3BOoWbNmdtecFx0drYYNG9o9Cnb69GlNmjRJktSvX7+i8ZtvvlmNGjXSl19+qaefflrh4eGSpPT0dL322mvy8PDQ6NGjr+dtAQAAoJQ5fCpTi3cka9GOZB1Ny7nk3CoVPDUyop7u7VRfdSv7uqhDAICZSvUuxiV1xZAk5eTkKDX1z2e2g4KCLjm3SpUq8vPzU2ZmppKSkmxq8+bN0+zZs9WzZ0/Vr19ffn5+Onz4sCIjI5WRkaEhQ4bo7rvvLprv4eGh2bNnq2/fvurWrZvuuusuVaxYUUuWLNEff/yh119/XaGhoVf9fi53wtmxY8eu+p4AAABwnuNpOVoRe1TLY48pJunsZec3r11Ro29soEHhdeTj6e78BgEAJUapDodKsvT09KLXFzue/q/Oh0MZGRk240OHDlVaWpq2bNmiDRs2KCsrS1WrVlXXrl01atQojRw50i4k69mzpzZt2qQpU6bo66+/Vl5enlq0aKHXXntN99xzzzW9n+Dg4Gu6DgAAAK5z4lyOftxzQstjjmrbodO63M4OXh5u6teylu7tVF8d6lcp0T98BQA4j9PDoZ9//lmSrvoI+CvRsGHDovuXNDk5/1uu6+Xlddn53t7ekqTs7Gyb8a5du6pr165X/fUjIiL0/fffX/V1AAAAKD2sVkNxR9K0Nj5Fa+NTFHfE8Ym1F2pZt6JGdAjWoDZ1VanCxU/KBQCUD04Ph5y5+XGFChVK7ObKPj4+Ra/z8vIuOz83N1eS5OtbMp/rvvBxtwsdO3ZMERERLuoGAACg/MrILdCmAyf/GwidVGpG7hVdV9XPSwNb19bwjsFqUYcNpgEA/8NjZU4SEBBQ9PrCR8UcyczMlHRlj6CZ4XL7JgEAAMA58gqs+nZnsjYnnNKBlAwlpKRf8pSxvwrw9lDflrU0sE0ddWlcTR7ubk7uFgBQGhEOOYmPj48CAwOVmpp62c2cz5w5UxQOsbcPAAAAztt/Il3jFkZr99FzV3xNgLeHejStroFt6qhHaHU2lwYAXBbhkBOFhYVp48aNSkhIUEFBwUWPs4+Pj7e5BgAAAOWb1WpozuY/NO2HfcorsF52fqPqfrq5WQ31alZTHRpUkScrhAAAV8Gp4VB+fr7i4uLk4eGhVq1aXfT0g9jYWEVHR2vUqFHObMflunbtqo0bNyozM1M7duzQDTfc4HDe+vXri1536dLFVe0BAACgBEo6naXnvonR1j9OX3SOp7tFNzSspl7NaqhXsxpqEOjnwg4BAGWN036ksGjRItWpU0cdO3ZU27ZtFRwcrC+//NLh3G+//VYPPPCAs1oxzeDBg4tez5071+Ecq9Wqzz//XJJUuXJl9ezZ0xWtAQAAoIQxDENfb09Sv5kbHQZDTWr467XbW2jBI520c3IfzX/oBj3YtSHBEADgujklHNq2bZtGjhypc+fOqXfv3rrtttt06tQp3XfffXr88ced8SVLpIiICHXr1k2S9Omnn+rXX3+1m/Pee+9p7969kqRnnnlGnp4cJQoAAFDepGbk6pEvdmjColhl5BbY1CwW6aGuDbX8qa66r3MDdWpUTf7e7A4BACg+TvlbZdq0aXJzc9PatWuLHpNKTEzUfffdp1mzZik7O1tz58696GNmJcWmTZuUkJBQ9OvU1NSi1wkJCZo3b57N/NGjR9vdY+bMmerSpYuys7PVp08fTZo0ST179lR2drYWLFigWbNmSZJCQ0M1fvx4p7wPAAAAlFyrdx/XxCVxOpWZZ1erW9lX7w5ro86Nq5nQGQCgvHBKOLR582YNHjzYZv+cevXq6aefftIDDzygzz//XIWFhfr8889LdEA0e/ZsffbZZw5rmzdv1ubNm23GHIVDbdu21cKFC3Xvvffq3LlzmjRpkt2c0NBQRUZGKiAgoFj6BgAAQMmXnpOvV5fv0Tc7HJ9sO7R9kKYMbK4AH1aWAwCcyynh0OnTp9WkSRP7L+bhoc8//1xeXl6aO3eurFarvvjiC2e0UKIMHDhQsbGxmjlzpiIjI5WcnCwvLy+FhIRo2LBhevLJJ1WhQgWz2wQAAIAL5BVY9cPu43rr+3gdOZttV6/m56U372ylvi1qmdAdAKA8cko4VKtWLaWkpDisWSwWffrppzIMQ/PmzZPValVISIgz2rhu8+bNs3t07FrVr19f06dP1/Tp04vlfgAAAChdMnML9PH6g/pqW6JSM+wfIZOk3s1r6h93tlKgv7eLuwMAlGdOCYeaNWtmczy7I59++qmkPwMYHqcCAABAWbbj8Gk9+3WMDp/Kclj39/bQlIHNNbR9UInedgEAUDY55bSyfv36KSEhwW5Pnr86v4Lo/vvvV3p6ujPaAAAAAEyVV2DV26viNezjXy8aDN3QsKq+f6abhnUIJhgCAJjCKSuHhg8frhMnTujkyZOXnGexWDRnzhzVr19fhw8fdkYrAAAAgCnij5/TuIUx2nvsnMN6RMOqurdTfQ1oVVtuboRCAADzWAzDMMxuAqVfcnKygoODJUlJSUkKCgoyuSMAAABzFFoNzd74u95bvV95hVa7evv6VfTmHa3UtBZbKwAArp4zPn87ZeUQAAAAUB4lnc7S+K9jtO3Qabuap7tF43qH6tHujeXOSiEAQAlCOAQAAABcp/ScfH21LVEz1xxQZl6hXb1pzQBNH9FGLepUMqE7AAAujXAIAAAAuEanMnI1a8Pv+nJrotJzC+zqFov0SLdGerZPqLw93E3oEACAyysR4dDp06eVlpamSpUqqWrVqma3AwAAAFzW6t3HNXFJnE5l5jmsB1Xx1fTh4YpoyPe3AICSzdRwaPHixXr55Ze1b9++orGKFSuqdevWCg8PL/qnZcuW8vT0NLFTAAAA4E/pOfl6ZfkeLdqRfNE5IzoE6+WBzeXvXSJ+FgsAwCWZ9rfVsmXLNGzYMFksFv31wLS0tDRt3LhRmzZtKhrz8PBQs2bNbAKjnj17mtE2AAAAyrFfD57Sc9/E6MjZbLuaxSL1bV5Lj/ZopLb1qpjQHQAA18a0cOjNN9+UJBmGoTvuuEO33HKLrFarEhISFBMTo5iYGJ05c0aSlJ+fr7i4OO3atUvz58+XxWJRQYH9M90AAACAM+TkF+rdH/bp081/6C8/1yzSu3lNTbotTA0D/VzfHAAA18m0cCg2NlYWi0V/+9vf9OGHHzqck5SUpOjoaEVHRysmJkbR0dH6448/bFYaAQAAAM6060iaxi2M1oGUDLuav7eHpgxsrqHtg2SxcDw9AKB0Mi0c8vPzU25urkaMGHHROcHBwQoODtbAgQOLxtLT0xUbG+uKFgEAAFBO5RVYFZ10Vmv2ntCcTX+owGr/w8lOjarq3WFtFFSlggkdAgBQfEwLh8LCwrR582a5u1/dkZ4BAQHq0qWLk7oCAABAeZZ0OktvrYrXz/EpysordDjHy8NNE/o21YNdGsrNjdVCAIDSz82sL3z//ffLMAz9/PPPZrUAAAAASPpzH8yFUYm6dcYGRcYeu2gw1KJORa14qqse6taIYAgAUGaYFg7dd999atOmjd5//30dP37crDYAAABQzp1Mz9XDn2/X3xfHKfMioZCbRXqqV4i+/VsXhdYMcHGHAAA4l2nhkJeXlxYuXCgPDw/dcsst2rNnj1mtAAAAoJxateuY+s7YoDV7UxzWm9euqIe6NlTk0900vk9TeXmY9u0zAABOY9qeQ5IUGhqqtWvXqkuXLurYsaOefPJJDR8+XO3btzezLQAAAJRxaVn5emXFbi357YjD+siOwXq+b1NV8/d2cWcAALiexTDxXPjXX39db775pnJzc2UYRtHxn9WrV1d4eLjCw8PVtm1bhYeHKzQ0lONBS7Dk5GQFBwdLkpKSkhQUFGRyRwAAAPYMw9CymKN6bcUepWbk2dUD/b301p2tdUvzmiZ0BwDA5Tnj87dpK4fmzp2ryZMn24ydz6lSUlL0448/6scffyyq+fr6qnXr1kWh0SOPPOLSfgEAAFC6paTn6IXFcVob7/gRsltb1NIbd7RktRAAoNwxLRz68MMPJUne3t4aP368br75Zrm7u+vAgQOKiYnRzp07FRMTo4yMDElSVlaWtmzZoq1bt8pisRAOAQAA4Ip9H3dMk76N05msfLtagLeHpg5qoTvb1WWlOgCgXDItHDpw4IAsFotefPFFvfTSS0Xj3bt3t5mXkJCgnTt3Kjo6uujfnG4GAACAK5GWna9Xlu3Wkp2O9xa6JaymXr29hepU9nVxZwAAlBymhUN+fn7KyspSnz59LjkvJCREISEhGjZsWNFYSorjpcAAAADAeZsTUvX8NzE6mpZjVwv099brg1vq1pa1TOgMAICSxbRwKDw8XD/++KOysrKu+toaNWo4oSMAAACUBTn5hXp7Vbzmbj7ksN6vZS29cUcrVfXzcm1jAACUUG5mfeEHH3xQhmFoxYoVZrUAAACAMiY2+az6f7DRYTAU4OOh90e00Uf3tCMYAgDgL0wLh4YPH64+ffro3//+t+Li4sxqAwAAAGVAQaFVM9cc0J0f/aKDJzPt6jc2rqYfxnbXHW2D2HQaAIALmBYOLVy4UC+//LJatWqlXr16admyZWa1AgAAgFKq0Gpow/6TGvLxr3p/zX4VWA2bureHm6YMbK75Y25g02kAAC7CtD2H7rrrrqKf2hiGoTvuuEN9+/bViBEjdNttt6l69epmtQYAAIAS7mxWnj7d9IcW7UjWMQcbTktSq7qV9P6INgqpEeDi7gAAKF1MC4ekP0Ohv77+4Ycf9MMPP0iS6tSpo7Zt2yo8PFxt27ZV27Zt1aBBA5M6BQAAQEmxbl+KJiyKVUp6rsO6u5tFT/YM0ZO9QuTpbtpCeQAASg3TwqEDBw4oOjpaMTExRf9OSkoqqh85ckRHjx5VZGRk0VilSpWKwqL33nvPjLYBAABgkqy8Ar0RuVf/tzXxonMaBfpp+ohwhQdXdl1jAACUchbjr8t3THbmzBmbsCg6Olp79+5VXl6ezTyLxaLCwkKTuoQjycnJCg4OliQlJSUpKCjI5I4AAEBZ8lviGT27MFqHTmU5rLeoU1HDOwRrRMdg+Xi6u7g7AABcxxmfv019rOxCVapU0U033aSbbrqpaKygoEB79uyxCY1iYmLMaxIAAAAuk1dg1Qc/HdBH6xJkdfAjzR6h1fX3W5upeZ2Krm8OAIAyokSFQ454eHiodevWat26te677z6z2wEAAICL7D+RrnELo7X76Dm7mq+nu17sH6Z7bqjH0fQAAFwn03boi4+P14033qinn35a+/fvN6sNAAAAlDBWq6HZG3/XgA83OQyG2tarrJXPdNO9neoTDAEAUAxMWzm0cOFCbdmyRYmJiZo2bZpZbQAAAKAEST6Tpee+idGW30/b1TzcLBrXO1SPdm8kD04hAwCg2JgWDq1evVoWi0WjR4+Wj4/PZecnJCRo5cqVCg8PV9euXeXmxjcEAAAAZYVhGFr82xG9smy30nML7OpNavjr/RHhalm3kgndAQBQtpkWDh06dEiS1L179yua36BBA7399ts6fvy4Vq5cqb59+zqxOwAAALjKqYxcTfo2Tj/sPmFXs1ikMV0a6rm+TTmFDAAAJzFt+c2pU6ckSTVr1ryi+R4eHhoyZIgMw9Dy5cud2RoAAABcZM2eE+o7Y4PDYKhuZV99+VAnvTSgOcEQAABOZFo4VLVqVUnSmTNnrviaHj16SJK2bNnilJ4AAADgGhm5Bfr7olg99Pl2pWbk2dWHtg/S92O7qXPjaiZ0BwBA+WJaONSsWTNJ0i+//HLF1zRo0ECSdPToUWe0BAAAABfY9sdp3TpjgxZuT7KrVfXz0sf3tte7w9qooo+nCd0BAFD+mBYO9enTR4ZhaNasWcrPz7+ia84fVXr+kTQAAACUHrkFhfrHyr0aMetXJZ/JtqvfElZDP4ztrltb1jKhOwAAyi/TwqExY8bI19dXSUlJeuSRR67omn379kmSAgICnNkaAAAAitmeo+d0+z836z8bfpdh2Nb8vNw1bUhrfTKqg6oHeJvTIAAA5Zhp4VD16tU1depUGYahzz//XAMGDFBqauolr/nkk08kSS1atHBFiwAAALhOBYVW/XvdQd3+r02KP55uV49oUFWrxnbX8I7BRavEAQCAa5kWDknS888/rzFjxsgwDH3//fdq3LixJk6cqB07dtjMS05O1v33369169bJYrFo8ODB5jQMAACAK1JQaNWS35LV+/0NentVvPILbZcLebm7adJtzfTVI50UXLWCSV0CAABJshjGhQt7Xe/NN9/UlClTVFhYWPQTI09PT9WsWVN5eXlKSUmRJBmGocaNGys6Olp+fn5mtowLJCcnKzg4WJKUlJSkoKAgkzsCAABmKCi0alnMUX24NkF/pGY6nBNWu6LeH9FGzWpVdHF3AACUfs74/G3qyqHzJk2apOjoaPXr10/SnyFQXl6ekpKSdOLECRmGIcMw1KxZM61cuZJgCAAAoATadSRN/T/YpGe/jnEYDLlZpMdvaqylT9xIMAQAQAniYXYD57Vo0UKRkZE6ceKEVqxYobi4OB0/flw5OTmqU6eOevbsqTvvvFPu7u5mtwoAAIC/KCi06j8bfteMNfvtHh87r1XdSpoysLk6NKjq4u4AAMDluCQc2rFjh9q3b39Fc2vWrKkxY8Y4uSMAAAAUh0OpmXr262j9lnjWYb11UCWNvaWJejatwYbTAACUUC4Jhzp27Kg6deqof//+GjhwoG655Rb5+Pi44ksDAADACQzD0JfbEvX6ir3Kzi+0qzerFaAJtzYlFAIAoBRw2WNlR48e1ezZszV79mz5+PioV69eGjhwoAYMGKA6deq4qg0AAABcp5RzOZqwOFbr9p20q7lZpL/dFKKnb24iL48Ssb0lAAC4DJf8jZ2cnKyPP/5Yt912m3x8fJSdna3IyEg9/vjjCg4OVvv27TV16lS7I+wBAABQskTGHlOfGRscBkMNqlXQN4/dqOf6NiUYAgCgFHH5UfbZ2dlas2aNVqxYocjISB09evTPRv673LhWrVo2j5/5+vq6sj1cI46yBwCgbEvLzteU73ZpafRRh/V7bqinSbeFyc+7xJx3AgBAmeSMz98uD4cutGPHDi1fvlwrVqzQb7/99mdT/w2KePys9CAcAgCg7NqckKrnvonRsbQcu1r1AG9NG9paPZvWMKEzAADKnzIZDv3V0aNHtWLFCi1fvlxr165Vdna2pP+FReHh4Ro4cKAGDhx4xaefwTUIhwAAKHty8gv11vfxmvfLIYf1/q1q6/XBLVXFz8u1jQEAUI6V+XDor3JycrRmzRotX778oo+fDRgwQH/729/Upk0bM1uFCIcAAChrYpPPatzCaB08mWlXC/Dx0Gu3t9Tt4XU4iQwAABdzxufvEvtQuI+PjwYMGKABAwZI+vPxs/Orinbu3Kljx45p9uzZqlu3LuEQAABAMSkotOpfPx/Uh2sPqMBq/zPELiHV9M7QNqpTmX0hAQAoK0psOHSh9u3bq3379poyZUrR42crVqxQhQoVzG4NAACgTPj9ZIbGfR2jmKSzdjVvDze90K+Z7u/cQG5urBYCAKAsKbGPlaF04bEyAABKL8Mw9MWWw3pz5V7l5Fvt6q3qVtL7I9oopEaACd0BAIC/KlePlQEAAMD5jqfl6PlFMdp4INWu5u5m0RM9Q/RUrxB5uruZ0B0AAHAFU8Oh7OxszZkzR7/88osyMjJUrVo1NW3aVK1bt1Z4eLhq165tZnsAAABl2rKYo3p56S6lZefb1RoF+mn6iHCFB1d2fWMAAMClTAuHUlNT1b17d+3bt++icwIDAxUeHm7zT7NmzTgVAwAA4DqczcrTy9/t1vKYow7r93eurxf6hcnXy93FnQEAADOYFg69+OKLio+PlyR5e3urSZMmys7O1qFDh1RYWChJOnnypNasWaM1a9YUXefj46NWrVppy5YtpvQNAABQmq3ff1ITFsXoxLlcu1rNit56Z2gbdQ+tbkJnAADALKaFQ5GRkbJYLGrevLl++OEH1alTR5KUl5enXbt2KTo6Wjt37tTOnTsVGxurjIwMSX8+ihYVFWVW2wAAAKXS4VOZ+tfPCfp6e7LD+qA2dfTa7S1VqYKnizsDAABmM/WxMkl6+eWXi4IhSfLy8lK7du3Url07m/kHDhzQzp07i0IjAAAAXN7eY+f073UHtSL2qKwOzqit5Oup1we31MA2deyLAACgXDAtHKpZs6aSk5PVpEmTK5rfpEkTNWnSRMOHD3dyZwAAAKVfVl6B3ly5V/O3JF50TrcmgXpnaBvVquTjws4AAEBJY9qZpJ07d5YknThxwqwWAAAAyqTfEs+o/webLhoM+Xi66bXbW+jzByMIhgAAgHnh0JgxY2QYhr777juzWgAAAChT8gutem/1Pg399y/6IzXTru7n5a5HujfS+ud76r7ODTgBFgAASDLxsbLevXvrzjvv1Ny5czVmzBh16NDBrFYAAABKvYSUdI1dGK1dR87Z1Xw93fVoj0YafWMDVa7gZUJ3AACgJDMtHFqyZIlefPFFHT9+XL1799Ynn3yioUOHmtUOAABAqWS1Gpr3yyG9vSpeuQVWu3p4cGW9PyJcDQP9TOgOAACUBqaFQ0OHDpXFYpG7u7sKCgo0YsQI9erVS8OHD9eAAQNUu3Zts1oDAAAoFY6ezdbzi2K0OeGUXc3DzaJnbm6ix29qLA9303YSAAAApYBp4ZAkGYahgoKCotdr167V2rVrJf15mlnbtm3Vtm1btWvXTm3btlXDhg3NbBcAAKBEMAxDS6OPaPJ3u5WeU2BXb1zdTzNGtFWroEomdAcAAEob08KhuLg4xcTE2Pzz15PLjh8/rlWrVmnVqlVFYxUrVlR4eLjatWun9957z4y2AQAATHUmM08vLd2lyLhjDusPdGmgv9/aTD6e7i7uDAAAlFYWwzAMs5s47+TJk4qOjrYJjOLj44tWF51nsVhUWFhoUpdwJDk5WcHBwZKkpKQkBQUFmdwRAABlz7p9KZqwKFYp6bl2tdqVfPTusDbqEhJoQmcAAMBVnPH529THyi5UvXp19e7dW7179y4ay8/P1+7du20Co9jYWBO7BAAAcK2svAK9uXKv5m9JdFi/o21dTR3UQpV8PV3cGQAAKAtMC4feeecdjR07Vp6el/4mxtPTU+Hh4QoPD3dNYwAAACXIb4lnNP7rGP2RmmlXq1zBU28MbqX+rTnIAwAAXDvTjq74+9//rrCwMC1ZssSsFgAAAEqs/EKrpq/ep6H//sVhMNQjtLp+GNudYAgAAFw308Ihd3d3/f777xo2bJh69Oih3377zaxWAAAASpSElHTd+dEv+mBtgqwX7A7p6+mu1we31LwHOqpmRR9zGgQAAGWKaeFQbGys+vbtK8MwtGnTJkVEROiBBx7Q0aNHzWoJAADAVFaroTmb/lD/DzYp7kiaXT08uLJWPtNN93aqL4vFYkKHAACgLDItHAoLC9P333+vFStWKDQ0VFarVZ9//rmaNm2q1157TdnZ2Wa1BgAA4HJHz2brvjlb9eqKPcotsNrUPNwsGt87VIse66yGgX4mdQgAAMoq08Kh82677Tbt2rVLM2bMUJUqVZSZmampU6cqNDRU8+fPN7s9AAAApzIMQ0t3HlHfGRu0OeGUXb1xdT99+7cueurmJvJwN/1bNwAAUAaViO8w3N3d9fTTT+vAgQN68skn5e7uriNHjuj+++9XRESENm/ebHaLAAAAxe5MZp6e/HKnxi6MVnpOgV39gS4NFPl0N7UKqmRCdwAAoLwoEeHQeVWqVNEHH3ygmJiYov2Itm/fru7du2vEiBE6dOiQ2S0CAAAUi3X7UtR3xgZFxh2zq9Wu5KP/e+gGTRnYQj6e7iZ0BwAAypMSFQ6dd34/opUrV6pZs2YyDEOLFi1SWFiYJk6cqPT0dLNbBAAAuCZZeQV6aWmcRs+NUkp6rl19cHgdrRrbXV1CAk3oDgAAlEclMhw679Zbb1VcXJzeeecdeXh4KDc3V9OmTVOTJk00a9YsGYZx+ZsAAACUEDsTz6j/B5s0f0uiXa1yBU/96+52mjGyrSr5eprQHQAAKK88zG7gQoZhaPfu3YqKitL27dsVFRWl2NhYFRQUyGKxyDAMpaSk6PHHH9dHH32kf/7zn+ratavZbQMAAFxUfqFVH/50QP9ad1CFVvsfbnUPra53hrZWzYo+JnQHAADKO9PDoQMHDhSFQFFRUYqOjlZWVlZR/a+rgywWi8LCwlSjRg2tX79esbGx6tGjhx577DFNnz5d3t7eZrwFAACAi0pISde4hTGKO5JmV/P1dNek/mG694Z6slgsJnQHAABgYjjUu3dvbd++XefOnSsau/AxsWrVqikiIkKdOnVSp06ddMMNN6hixYqSpKioKE2YMEHr16/Xxx9/rJ07d+rHH3+Un5+fS98HAACAI1aroc9+PaS3vo9XboHVrh4eXFnvjwhXw0C+dwEAAOYyLRz66aefbH7t4eGh1q1bFwVBnTp1UkhIyEWv79ixo37++WfNnj1bTzzxhLZu3ap3331XU6ZMcXbrAAAAl3T0bLaeXxSjzQmn7GoebhY9c3MTPX5TY3m4l+jtHwEAQDlhWjhUp04dmyCoQ4cO8vG5+ufsH3roIR09elRTp07VwoULCYcAAIBpDMPQd9FH9fJ3u5SeU2BXb1zdTzNGtFWroEomdAcAAOCYaT+uSk5O1qJFi/Tcc8+pa9eu1xQMnTdo0CBJ0qFDh4qpu5LlyJEjmjFjhvr06aN69erJy8tLtWrV0pAhQ7R161a7+fPmzZPFYrnkPzfffLMJ7wQAgLLrbFaenvxqp8YujHYYDI2+sYEin+5GMAQAAEoc0zekLg4BAQGSpNzcXJM7cY4PP/xQb7/9tho3bqzevXurRo0aOnDggJYuXaqlS5fqq6++0vDhw4vmh4eHX3QF1aJFi7R792717dvXVe0DAFDmrduXogmLYpWSbv+9SO1KPnpnaBt1bRJoQmcAAACXZzEu3AW6mN1zzz1q06aNwsPDFR4erho1ahT718jOztb8+fO1c+dOffTRR8V+f7MtWbJE1atXV7du3WzGN27cqJtvvlkBAQE6evToZU9ry8vLU506dZSWlqbk5GTVrFmz2HpMTk5WcHCwJCkpKUlBQUHFdm8AAEqinPxCLYs5qkXbk7Xt0GmHcwaH19Ert7dUJV9PF3cHAADKKmd8/nZ6OOTm5mZzNGvNmjUVHh5uExiFhoZyfOs16tu3r1avXq2oqCh16NDhknMXLlyokSNHavDgwfr222+LtQ/CIQBAebL90Gk9+3WMEk9nOaxX8vXUG3e01IDWdVzcGQAAKOuc8fnb6Y+VhYaGKiEhQVbrn0e4Hj9+XD/88IN++OGHojm+vr5q1aqVTWjUunVrVahQwdntXVZKSoq2bdumbdu2KSoqSlFRUTp16s+TR+6//37Nmzfviu+VmJioDz74QJGRkUpMTJS3t7dCQkI0fPhw/e1vf7um9+vp+edPIj08Lv8/5aeffirpz028AQDA1csrsGrGmv36eP1BWS/y47XuodX1ztDWqlnx2vdTBAAAcCWnrxySpJycHMXFxSkmJqbon9jYWJ07d+5/jVywcshisSgkJMRulVHt2rWd3a5dHxdzNeFQZGSk7rnnHqWlpTmsN23aVCtXrlSjRo2uuLfExESFhoaqSpUqSk5Olru7+0XnHj58WI0aNVLt2rV1+PDhS869FqwcAgCUdftPpGvsgmjtOXbOYb2qn5fG9Q7VvTfUY0U0AABwmlK5ckiSfHx81LFjR3Xs2NFm/ODBg3riiSe0evVqXZhRGYah/fv368CBA/rmm2+KxgMDA9W2bVutWrXKFa3bCA4OVlhYmFavXn1V18XExGj48OHKysqSv7+/Jk6cqJ49eyo7O1sLFizQJ598on379ql///6KioqSv7//Ze+Zn5+v++67T7m5uZo2bdplw565c+fKarXqgQceKPZgCACAssxqNTRn8x+a9sM+5RVY7erhwZX1aPdGujmsprw8TDsIFgAA4JqZdlqZ1WrVs88+qx9//FGNGzfW888/r1atWqlixYo6fvy4tm7dqmXLlmnbtm021508eVI//vijy/qcPHlyUbBVs2ZNHTp0SA0bNryqe4wdO1ZZWVny8PDQ6tWr1blz56Jar1691KRJE02YMEHx8fGaPn26Jk+efMn7Wa1WPfjgg9qwYYMefvhh3XfffZedP3fuXFksFj344INX1TsAAOXZkbPZev6bGP1y8JRdzcPNonG9Q/VYj8Zyd2OlEAAAKL1c8liZIzNnztS4ceMUERGh9evXX/SkraioKD399NPaunWr/P39dddddyk+Pl7r1693ccd/+ms4dCWPlUVFRSkiIkKS9Oijj+rjjz+2m2O1WtWyZUvt3btXVapU0YkTJ4r2ErqQYRh66KGHNGfOHN1777367LPP5OZ26Z9S/vDDD7r11lt18803a82aNVfwLq8ej5UBAMoSwzC0NPqIJn+3W+k5BXb1JjX89f6IcLWsW8mE7gAAQHnmjM/fpq19/vTTT2WxWDRx4sRLHsHesWNHbdq0SaNGjVJGRoYOHjxoWjB0LZYuXVr0+oEHHnA4x83NTaNGjZIknTlzRuvWrXM4z2q1asyYMZozZ47uuusuzZs377LBkMRG1AAAXI0zmXl68sudGrcwxmEwNKZrQy1/qivBEAAAKDNMC4cOHjwoSapXr95l57q7u2v27Nlq06aNfv75Z3311VfObq/YbNy4UZLk5+en9u3bX3Rejx49il5v2rTJrm61WvXQQw9p7ty5GjFihL744osr2jvo1KlT+u6771S1alXdcccd1/AOAAAoP9btS1HfGRsUGXfMrlanko++fOgGvTyguXw82b8PAACUHabtOeTv76+cnBwlJSWpbdu2l53v4eGhp556Sg899JC++uor3XXXXS7o8vrt3btXkhQSEnLJ4+abNWtmd81551cMzZs3T8OGDdP8+fOveFPpL774Qnl5ebr33nsvuULrcpKTky9ZP3bM/ptoAABKi6y8Av1jZby+2HLYYf2OtnU1dVALVfJ1/Ng3AABAaWZaONSxY0d9//33+vLLLzVo0KAruqZ169aSpN9++82ZrRWbnJwcpaamStJlnwGsUqWK/Pz8lJmZqaSkJJvaq6++qnnz5snf31+hoaF6/fXX7a4fPHiwwsPD7caL65Gy888zAgBQ1kQnndWzC6P1e2qmXa1yBU+9MbiV+reubUJnAAAArmFaODR69GitXLlS33zzjYYMGaJhw4Zd9pr09HRJfz4qVRqc71fSFR1Pfz4cysjIsBk/dOiQJCkjI0NvvPGGw2sbNGhgFw5t27ZNu3btUkREhFq1anV1zQMAUMblF1r1z7UJ+ufPCSq02p/P0T20ut4Z2lo1K/qY0B0AAIDrmBYODR06VLfccovWrFmju+++W/v379eECRMuekqXJM2dO1eSVKFCBVe1eV1ycnKKXnt5eV12/vnHvrKzs23G582bd9lT0RyJiIhQcR1Gd+FqpgsdO3as6FQ2AABKuoMnM/TswmjFJKfZ1Xw83fRi/+a694Z6slg4oh4AAJR9poVDkvTNN9/o5ptv1m+//abJkydr9uzZeuSRRzRo0CC1aNGiaN7Bgwf16quvav78+bJYLLrxxhtN7PrK+fj87yeNeXl5l52fm5srSfL19XVaT9eKo+kBAGWBYRj6Ysthvblyr3LyrXb1NkGVNH1EuBpXv/yKXwAAgLLC1HCoUqVK2rx5s5544gnNmTNHhw8f1ksvvaSXXnpJvr6+CgwM1OnTp5WZ+b89ADw8PPTiiy+a2PWVCwgIKHp94aNijpx/n1fyCBoAALg6Saez9OLSXdqw/6Rdzd3Noqd6heiJniHydDftMFcAAABTmP7dj7e3t2bPnq2NGzfq5ptvlmEYMgxDWVlZSkxMVEZGRtFY5cqVtWDBAnXq1Mnstq+Ij4+PAgMDJV3+tK8zZ84UhUNs/gwAQPE5mZ6rqct2q9d76xwGQ40C/bTk8Rs19pZQgiEAAFAumbpy6K+6dOmiH3/8UceOHdOqVasUExOj48ePKzs7WzVq1FCnTp00bNgwVaxY0exWr0pYWJg2btyohIQEFRQUXPQ4+/j4eJtrAADA9TEMQ4t/O6JXlu1Wem6Bwzn3d66vF/qFydfL3cXdAQAAlBwlJhw6r3bt2nrggQfMbqPYdO3aVRs3blRmZqZ27NihG264weG89evXF73u0qWLq9oDAKBMOpaWrVeW7dGq3ccd1msEeOudYW3UI7S6izsDAAAoeVg77WSDBw8uen3+tLULWa1Wff7555KkypUrq2fPnq5oDQCAMufXg6c0eu42dXlrrcNgyMvdTaNvbKDV47oTDAEAAPwX4ZCTRUREqFu3bpKkTz/9VL/++qvdnPfee0979+6VJD3zzDPy9PR0aY8AAJR22XmFmvLdLt31yRat23dSVsN+zu3hdbT2uR6aOqiFKlfwcn2TAAAAJZTTHytLTEyUJNWtW1fu7sX7PH9hYaGOHDkiSapXr16x3vu8TZs2KSEhoejXqampRa8TEhI0b948m/mjR4+2u8fMmTPVpUsXZWdnq0+fPpo0aZJ69uyp7OxsLViwQLNmzZIkhYaGavz48U55HwAAlFVxyWkau3CnDp7MdFiv6uelN+9opVtb1nJxZwAAAKWDxTAMBz9bKz5ubm5yc3NTbGysmjdvXqz33r17t1q1aiU3NzcVFDjeaPJ6jR49Wp999tkVz7/Yf87ly5fr3nvv1blz5xzWQ0NDFRkZqZCQkGvq02zJyclFp6wlJSUpKCjI5I4AAGVdQaFVH68/qBlrDqjAwVIhf28P3R5eR8/c0kQ1AnxM6BAAAKD4OePzt0s2pHZy/uT0+xeHgQMHKjY2VjNnzlRkZKSSk5Pl5eWlkJAQDRs2TE8++aQqVKhgdpsAAJQKiaeyNO7raO04fMau5uPppuf7NtNdEcGq4FXizt4AAAAocVyycshisWjXrl3FfkT7+ZVDFotFhYWFxXpvXB1WDgEAXMEwDH2zPVmvLN+tzDz7v/tbB1XS+yPC1bi6vwndAQAAOF+pXTkkSX369Cn2jZbz8/OL9X4AAKDkOpWRq4lL4rR6zwm7mptFerJniJ66uYk83TlvAwAA4Gq47LGy8xtHAwAAXK218Sc0YVGcUjNy7Wr1q1XQ9OHhal+/igmdAQAAlH5OD4fuv/9+Z38JAABQRmXlFeiNyL36v62JDusjOwbr5QHN5efN3kIAAADXyunfSc2dO9fZXwIAAJRB0UlnNW5htP5ItT+ivpqfl94a0lq9m9c0oTMAAICyhR+zAQCAEqWg0Kp//pygD9cmqNDBEfU3N6uht4a0VvUAbxO6AwAAKHsIhwAAQInxR2qmxi6MVkzSWbuar6e7Xh7QXHdFBMtisbi+OQAAgDKKcAgAAJjOMAx9uS1Rr6/Yq+x8+yPqw4Mr6/0R4WoY6GdCdwAAAGUb4RAAADDVyfRcvbA4Vj/Fp9jV3N0serpXEz3Rs7E8OKIeAADAKQiHAACAaX7cc0IvLI7Vqcw8u1rDQD+9PyJc4cGVXd8YAABAOUI4BAAAXC4zt0CvLt+jhduTHNbv7VRPk24LUwUvvlUBAABwNr7jAgAALrXj8BmNWxitxNNZdrVAf2+9M7S1ejarYUJnAAAA5RPhEAAAcIn8Qqs++OmA/vVzghycUK8+zWvqH3e2UjV/jqgHAABwJcIhAADgdAdPZmjcwmjFJqfZ1fy83DVlYAsN6xDEEfUAAAAmIBwCAABOYxiGvthyWG+u3KucfKtdvX39Knp/eLjqVatgQncAAACQCIcAAICTpJzL0fOLYrV+/0m7moebRWNvaaLHenBEPQAAgNkIhwAAQLFbteuYJi6J05msfLta4+p+mjGirVoFVTKhMwAAAFyIcAgAABSb9Jx8TV22R4t/S3ZYv79zfb3QL0y+Xu4u7gwAAAAXQzgEAACKxbY/TuvZr6OVfCbbrlYjwFvvDGujHqHVTegMAAAAl0I4BAAArktegVXvr9mvj9cflOHgiPp+LWvpzTtaqYqfl+ubAwAAwGURDgEAgGu2/0S6xi6I1p5j5+xq/t4eemVQC93Zri5H1AMAAJRghEMAAOCqWa2G5v1ySG+tildegf0R9RENquq94W0UXJUj6gEAAEo6wiEAAHBVjqfl6PlFMdp4INWu5ulu0bO9m+qR7o3k7sZqIQAAgNKAcAgAAFyxFbFH9eK3u5SWbX9EfZMa/np/RLha1uWIegAAgNKEcAgAAFxWWna+pi7brW93HnFYf7BLQ024tal8PDmiHgAAoLQhHAIAAJf068FTGv91tI6m5djValX00bvD2qhrk0ATOgMAAEBxIBwCAAAO5RYU6r3V+/XJxt8dHlE/oHVtvT64pSpX4Ih6AACA0oxwCAAA2Ik/fk5jF0Qr/ni6XS3Ax0OvD26pQW3qcEQ9AABAGUA4BAAAilithuZs/kPTVu1TXqH9EfWdGlXVe8PDVbeyrwndAQAAwBkIhwAAgAoKrVoRe0wfrUvQ/hMZdnUvdzc937epxnRtKDeOqAcAAChTCIcAACjnog6d1oRFsfojNdNhvVmtAL0/IlxhtSu6uDMAAAC4AuEQAADlVF6BVTN/2q9/rzsoq4MNpy0W6eFujfRs71COqAcAACjDXBoOHTx4UF9++aXi4+NVWFio2rVrq3v37rrtttvk7e3tylYAACjXElLSNXZhtHYdOeew3qSGv165vYVubMwR9QAAAGWdy8Kh//znP3r66adVUFBgM/7BBx8oODhYH3/8sW699VZXtQMAQLlkGIa+2HJYb0TuVW6B/YbTLetW1FO9mqh3WE32FgIAACgnXBIORUVF6YknnpDVav9NqCQlJiZq0KBBWrRokQYNGuSKlgAAKHdSzuXo+UWxWr//pF3N092icb1D9Wj3xnInFAIAAChXXBIOffjhh7JarbJYLOrfv7/GjBmjOnXq6MiRI1qxYoW++OIL5efn68EHH1R8fLwCA1nCDgBAcVq165gmLonTmax8u1pIDX/NGBGulnUrmdAZAAAAzOaScGjTpk2yWCzq16+fli1bVjTesWNHDR48WKNGjVLfvn115swZzZo1S5MmTXJFWwAAlHnpOfl6ZfkeLdqR7LA++sYGeqFfMzacBgAAKMfcXPFFjh07Jkl69NFHHda7d++u8ePHyzAMLV682BUtAQBQ5m0/dFq3fbDRYTBUI8Bbnz0YoamDWhAMAQAAlHMuCYdyc3MlSQ0aNLjonLvuukuSFBcXp7y8PFe0BQBAmZRXYNU7P8Rr+H9+VdLpbLt6v5a19MPY7uoRWt2E7gAAAFDSuPQoe3f3i/9kMiQkRJJUWFiokydPqm7duq5qCwCAMiMhJUNjF+50eES9v7eHpg5qoSHt6spiYdNpAAAA/Mml4dCleHt7F71OT083sRMAAEqf80fUv7lyr3Ly7U8H7VC/it4fEa7gqhVM6A4AAAAlmUvDoSv9KeXFjrwHAAD2LnVEvYfbn0fUP9aDI+oBAADgmEvDoa5du6pNmzYKDw8v+qd58+by8CgxC5gAAChVOKIeAAAA18tlqYxhGDpz5ozWr1+v9evXF417enqqefPmCg8PLxrLz7f/BhcAAPxPRm6BXlm2W99wRD0AAACuk0vCoY8++kjR0dGKjo7Wrl27lJWVVVTLy8tTTEyMYmJiih4769Chgxo3bqzWrVurdevWatOmjVq3bq369eu7ol0AAEq07YdOa9zX0Q5PIqsR4K13hrXhJDIAAABcMYthGIYrv6BhGNq3b19RWBQdHa2YmBidOHHCvrkL9igKCAhQq1at1KZNG/3zn/90Vcu4AsnJyQoODpYkJSUlKSgoyOSOAKDsyS+0auaaA/poXYKsDv72vrVFLf3jzlaq4ufl+uYAAADgEs74/O3ycOhiTpw4oZ07d9qERgkJCQ43p7ZYLCosLDShS1wM4RAAONfBkxkauyBacUfS7Gr+3h6aMrC5hrYP4oh6AACAMs4Zn79LzE7QNWvW1K233qpbb721aCw7O1sxMTE2gdGuXbuUnW2/jB4AgLLIMAwtjErSK8v3KDvf/gcjHFEPAACA61ViwiFHfH191alTJ3Xq1KlozDAM7d+/38SuAABwjbNZeXphcZxW7T5uV+OIegAAABSXEh0OOWKxWNS0aVOz2wAAwKl+PXhK4xZG6/i5HLta4+p+mjGirVoFcUQ9AAAArl+pC4cAACjL8gutev/H/fr3+oNytCvgPTfU00v9m8vXiyPqAQAAUDwIhwAAKCEOpWbqmQU7FZNsv+l05QqeentIa/VtUcuEzgAAAFCWEQ4BAGAywzC0aEeypizbraw8+02nu4RU0/Th4apZ0ceE7gAAAFDWEQ4BAGCitOx8Tfo2TpGxx+xqnu4WPdenqR7u1khubDoNAAAAJyEcAgDAJNv+OK1xC6N15Gy2Xa1RoJ9mjmTTaQAAADgf4RAAAC5WUGjVBz8d0D9/TpDVwabTIzsGa/LA5qrgxV/TAAAAcD6+6wQAwIWSTmfpmQU79VviWbtaJV9PvXVnK/VrVdv1jQEAAKDcIhwCAMBFvt2ZrJeX7lZGboFdrVOjqpo+PFx1Kvua0BkAAADKM8IhAACc7FxOviYv3aWl0Uftah5uFo3rHarHejSWO5tOAwAAwASEQwAAONGOw2f0zIKdSj5jv+l0/WoVNHNkW4UHV3Z9YwAAAMB/EQ4BAOAEhVZD//o5QTN/OqBCB7tOD20fpKmDWsjfm7+KAQAAYC6+IwUAoJgln8nSuIXRijp0xq4W4OOhN+9opYFt6pjQGQAAAGCPcAgAgGK0LOaoXvw2Tuk59ptOd2xQRe+PCFdQlQomdAYAAAA4RjgEAEAxyMgt0JTvdmvxb8l2NXc3i565uYn+dlNjebi7mdAdAAAAcHGEQwAAXKfopLN6ZsFOHT6VZVcLruqrGSPaqn39KiZ0BgAAAFwe4RAAANeo0Gro4/UH9f6P+1XgYNPpO9rW1au3t1CAj6cJ3QEAAABXhnAIAIBrcPRstsYtjNbWP07b1fy9PfT64JYa3LauCZ0BAAAAV4dwCACAq/R93DG9sCROadn5drW29Spr5oi2qleNTacBAABQOhAOAQBwhbLyCvTq8j1aEJVkV3OzSE/2DNHTNzdh02kAAACUKoRDAABcgbjkND2zYKd+T820q9Wt7Kv3R4QromFVEzoDAAAArg/hEAAAl2C1Gvpk4+96d/U+5Rfabzo9oHVtvXFHK1XyZdNpAAAAlE6EQwAAXMSJczl69utobU44ZVer4OWuV29vqSHt6spisZjQHQAAAFA8CIcAAHBg9e7j+vviWJ3Jst90uk1QJc0c2VYNAv1M6AwAAAAoXoRDAAD8RU5+of6xcq8++/WwXc1ikR7v0VjjeofKk02nAQAAUEYQDgEA8F9/pGbqyS9/0+6j5+xqtSv5aPrwcHVuXM2EzgAAAADnIRwCAEDSd9FHNGlJnDLzCu1q/VrW0j/ubKXKFbxM6AwAAABwLsIhAEC5lpVXoKnLduvr7cl2NW8PN00d1EIjOwaz6TQAAADKLMIhAEC5te94up788jcdSMmwqzWp4a9/3t1OTWsFmNAZAAAA4DqEQwCAcscwDC2IStLUZbuVW2C1qw/vEKSpg1qoghd/TQIAAKDs47teAEC5kp6Tr0nf7tLymKN2NT8vd715ZyvdHl7XhM4AAAAAcxAOAQDKjdjks3rqq506fCrLrtaiTkX98+52ahjoZ0JnAAAAgHkIhwAAZZ5hGJqz+ZDe+n6v8gsNu/roGxto4m3N5O3hbkJ3AAAAgLkIhwAAZdqZzDw9vyhWa/aesKtV9PHQtKFtdGvLWiZ0BgAAAJQMhEMAgDIr6tBpPf3VTh1Ly7Grta1XWR/e1VZBVSqY0BkAAABQchAOAQDKHKvV0L/XH9T0H/er0Gr/GNljPRprfJ9Qebq7mdAdAAAAULIQDgEAypSU9Bw9uzBGmxJS7WrV/Lz03vA2uqlpDRM6AwAAAEomwiEAQJmx6UCqxi6MVmpGrl2tc6NqmjEyXDUr+pjQGQAAAFByEQ4BAEq9gkKrZqw5oH+tS5BxwVNkbhZp7C2heqJniNzdLOY0CAAAAJRghEMAgFLt6NlsPf3VTm0/fMauVrOitz4Y2VY3NKpmQmcAAABA6cBOnCXckSNHNGPGDPXp00f16tWTl5eXatWqpSFDhmjr1q1XdI9p06bJYrHIYrFoy5YtTu4YAFznxz0ndNsHGx0GQ72a1dD3z3QnGAIAAAAug5VDJdyHH36ot99+W40bN1bv3r1Vo0YNHThwQEuXLtXSpUv11Vdfafjw4Re9fu/evZo8ebL8/PyUmZnpws4BwHlyCwr11vfxmrv5kF3Nw82iF/o104NdGsqNx8gAAACAyyIcKuEiIiK0YcMGdevWzWZ848aNuvnmm/X444/r9ttvl7e3t921hYWFuv/++9WmTRuFhoZq/vz5rmobAJzm8KlMPfnlTsUdSbOrBVf11Yd3tVN4cGXXNwYAAACUUjxWVsLdeeeddsGQJHXr1k09e/bU6dOnFRcX5/Dat99+WzExMZozZ47c3d2d3SoAON2ymKPq/8Emh8HQba1qacVT3QiGAAAAgKtEOHQJKSkpWrFihSZPnqx+/fopMDCwaO+e0aNHX9W9EhMT9dxzzyksLEx+fn6qWrWqIiIi9O677yorK+ua+vP09JQkeXjYLwDbtWuXXnnlFb300ktq0aLFNd0fAEqK7LxCTVwSq6e/2qmM3AKbmpeHm14f3FL/urudKvl6mtQhAAAAUHrxWNkl1KxZs1juExkZqXvuuUdpaf/7SXdWVpaioqIUFRWl2bNna+XKlWrUqNEV3zMxMVFr1qxRrVq11KpVK5taQUGBRo8erbCwML3wwgvF8h4AwCwHTqTriS9/0/4TGXa1RtX99M+72ql5nYomdAYAAACUDawcukLBwcHq06fPVV8XExOj4cOHKy0tTf7+/nrjjTf0yy+/6KefftLDDz8sSdq3b5/69++vjAz7Dz6O5Ofn67777lNubq6mTZtm98jYm2++WfQ42fnVRQBQ2hiGoa+jkjTwn5scBkND2gVp+ZNdCYYAAACA68TKoUuYPHmyOnbsqI4dO6pmzZo6dOiQGjZseFX3GDt2rLKysuTh4aHVq1erc+fORbVevXqpSZMmmjBhguLj4zV9+nRNnjz5kvezWq168MEHtWHDBj388MO67777bOoxMTF6/fXX9dxzz6ldu3ZX1SsAlBQZuQV68ds4fRd91K5Wwctdr93eUkPaB5nQGQAAAFD2sHLoEl555RUNGDDgmh8vi4qK0rp16yRJY8aMsQmGzhs/frzCwsIkSTNmzFB+fv5F72cYhh5++GHNnz9f9957rz7++GO7Offff78aN26sqVOnXlPPAGC2XUfSNOCDjQ6DoWa1ArT8qa4EQwAAAEAxIhxyoqVLlxa9fuCBBxzOcXNz06hRoyRJZ86cKQqTLmS1WjVmzBjNmTNHd911l+bNmyc3N/v/+WJiYhQfHy8fH5+izbMtFos+++wzSVLnzp1lsVhsegOAksAwDM3b/Ifu/OgXHTplv1H/fZ3qa+kTXdS4ur8J3QEAAABlF4+VOdHGjRslSX5+fmrfvv1F5/Xo0aPo9aZNm9S7d2+butVq1UMPPaS5c+dqxIgR+uKLLy56NP2YMWMcjm/YsEEHDhzQoEGDVL16dTVo0OAq3w0AOE9aVr6eXxSj1XtO2NUCfDw0bUhr9WtV24TOAAAAgLKPcMiJ9u7dK0kKCQlxeNz8ec2aNbO75rzzK4bmzZunYcOGaf78+RcNhiRp9uzZDsdHjx6tAwcOaOLEierUqdPVvA1JUnJy8iXrx44du+p7AoAk7Th8Rk9/tVNHzmbb1doEV9Y/72qr4KoVTOgMAAAAKB8Ih5wkJydHqampkqSgoEvvjVGlShX5+fkpMzNTSUlJNrVXX31V8+bNk7+/v0JDQ/X666/bXT948GCFh4cXW++OBAcHO/X+AMofq9XQfzb8rndX71Oh1bCrP9K9kZ7r01ReHjwBDQAAADgT4ZCTpKenF73297/8/hjnw6ELj7M/dOiQJCkjI0NvvPGGw2sbNGjg9HAIAIpTakaunv06Rhv2n7SrVangqenDw9WzWQ0TOgMAAADKH8IhJ8nJySl67eXlddn53t7ekqTsbNvHKubNm6d58+Zddz/Xe58LVzRd6NixY4qIiLjm+wMoP35JSNUzC6N1Mj3XrhbRsKo+GNlWtSr5mNAZAAAAUD4RDjmJj8//Ptjk5eVddn5u7p8fknx9fZ3W0/W43KNxAHA5BYVWffDTAX34c4KMC54is1ikp3o10dO9QuThzmNkAAAAgCsRDjlJQEBA0esLHxVzJDMzU9KVPYIGAKVN4qksPfdNjLYdOm1XqxHgrRkjw3Vj40ATOgMAAABAOOQkPj4+CgwMVGpq6mVP+jpz5kxROMTGzwDKEqvV0Pyth/XW9/HKyiu0q3cPra7pw9so0N/bhO4AAAAASIRDThUWFqaNGzcqISFBBQUFFz3OPj4+3uYaACgLEk9l6flFMdr6h/1qIQ83i57v21QPd2skNzeLCd0BAAAAOI+NHZyoa9eukv58ZGzHjh0Xnbd+/fqi1126dHF6XwDgTFaroc9+OaS+MzY4DIbqVvbVwkc769EejQmGAAAAgBKAcMiJBg8eXPR67ty5DudYrVZ9/vnnkqTKlSurZ8+ermgNAJzi8KlM3fXJFk1ZtlvZ+faPkd19Qz39MK672tevYkJ3AAAAABwhHHKiiIgIdevWTZL06aef6tdff7Wb895772nv3r2SpGeeeUaenp4u7REAioNhGJq/5bBunbHxoquF5o+5QW/e0Ur+3jzRDAAAAJQkfId+CZs2bVJCQkLRr1NTU4teJyQkaN68eTbzR48ebXePmTNnqkuXLsrOzlafPn00adIk9ezZU9nZ2VqwYIFmzZolSQoNDdX48eOd8j4AwJly8gs16ds4LfntiMP6PTfU08TbwgiFAAAAgBLKYhiGYXYTJdXo0aP12WefXfH8i/2nXL58ue69916dO3fOYT00NFSRkZEKCQm5pj5LguTk5KKT1pKSkhQUFGRyRwBcIel0lh6bv0O7j9r/+Va3sq+mDW2tLiEcUQ8AAAAUF2d8/ubHuC4wcOBAxcbGaubMmYqMjFRycrK8vLwUEhKiYcOG6cknn/z/9u47OsoyYf/4NSmkEjoJPVSJSE8iCogBRVBURJoKGEFewLWABV9ZF9FVAQs/dnVRMEAQ1BVRUcFlrQgoYugghLL0AIYACaSXeX5/+GaWYRJIcJJnJs/3cw7nTO77nskVuE8yuXiKgoODzY4JAOWybl+aHvlgs85mF7jMcbQQAAAA4D04cghuwZFDgHUYhqG5aw7olVXJsl/0EySkmq9eH9pR/a5pYE44AAAAoIrjyCEAgKky8wo1edk2fbnjpMtci7ohmjuyq1qHVzchGQAAAIArRTkEACiTA6cyNW7xJu1LzXSZu/nqcL0+tKPCArnjIgAAAOBtKIcAAJf19a7f9PiHW3U+r9Bp3GaTHr+pjf4U10o+PjaT0gEAAAD4IyiHAAClstsNzf52n/7+7T6XubBAP/3tns6Ku6q+CckAAAAAuAvlEACgRBk5BZr4zy36fs8pl7m2EdU1d2RXNasTYkIyAAAAAO5EOQQAcJF88pzGLd6kw6ezXebu6NhQM+5ur+Bq/AgBAAAAqgLe2QMAnHyx7bgmL9uunIIip3FfH5ue6d9WY3o0l83G9YUAAACAqoJyCAAgSSossmvmqmS9s/agy1ydkGp6894uuq5lHROSAQAAAKhIlEMAAJ3OzNMjH2zRT/857TLXsXENvTWiqxrWDDIhGQAAAICKRjkEABa3/Vi6xi/epOMZuS5zw2OaaNod7RTo72tCMgAAAACVgXIIACxs6cajenb5TuUX2p3G/X1tev6Oa3TvtU1NSgYAAACgslAOAYAF5Rfa9fwXv+q9DUdc5iLCAjVnRBd1aVrLhGQAAAAAKhvlEABYzG/ncjVhySZtPpLuMhcbWVv/uK+L6lUPqPxgAAAAAExBOQQAFpJ06Iweem+zTp3Pc5l7oHukptwaJX9fHxOSAQAAADAL5RAAWIBhGHp3/WH9dcUuFdoNp7lAfx9NH9Red3VubFI6AAAAAGaiHAKAKi63oEhTPt2hTzanuMw1rhWkuSO7ql3DGiYkAwAAAOAJKIcAoAo7eiZb45ds0q/Hz7nM9WxdV2/c01k1g6uZkAwAAACAp6AcAoAqat2+ND3ywWadzS5wmftTXEs9fvNV8vWxmZAMAAAAgCehHAKAKsYwDM1dc0CvrErWRZcXUkg1X70+tJP6XRNhTjgAAAAAHodyCACqkMy8Qk1etk1f7jjpMteiXojmjeyqVvWrm5AMAAAAgKeiHAKAKuLAqUyNW7xJ+1IzXeZuvjpcs4Z2VPVAfxOSAQAAAPBklEMAUAV8ves3Pf7hVp3PK3Qat9mkJ25uo4dubCUfri8EAAAAoASUQwDgxex2Q7O/3ae/f7vPZa5GkL/+NryTbryqvgnJAAAAAHgLyiEA8FIZ2QWa+OEWfb/nlMtc24jqmjcyWk3rBJuQDAAAAIA3oRwCAC+UfPKcxi3epMOns13m7uzUUNMHtVdwNb7FAwAAALg8fnMAAC/z+bbjenrZduUUFDmN+/rYNOXWKI3uHimbjesLAQAAACgbyiEA8BKFRXbNXJWsd9YedJmrE1JNb97bRde1rGNCMgAAAADejHIIALzA6cw8Pfz+Fq0/cNplrmOTmnp7RBc1qBFkQjIAAAAA3o5yCAA83PZj6Rq/eJOOZ+S6zA2PaaLn72ynAD9fE5IBAAAAqAoohwDAgy1NOqpnP9up/EK703g1Xx89f2c73RPb1KRkAAAAAKoKyiEA8ED5hXY9/8Wvem/DEZe5iLBAvTWiizo3rWVCMgAAAABVDeUQAHiY387lavySTdpyJN1lLrZ5bf3j3i6qVz2g8oMBAAAAqJIohwDAgyQdOqMJSzYrLTPPZW509+Z65ta28vf1MSEZAAAAgKqKcggAPIBhGHp3/WH9dcUuFdoNp7lAfx/NGNRBAzs3MikdAAAAgKqMcggATJZbUKQpn+zQJ1tSXOaa1A7S2yO6ql3DGiYkAwAAAGAFlEMAYKKjZ7I1fskm/Xr8nMvcDW3q6e/DO6lmcDUTkgEAAACwCsohADDJ2n2n9MgHW5SeXeAy93BcK026uY18fWwmJAMAAABgJZRDAFDJDMPQ2z8c0Kv/TtZFlxdSaICfXhvSUf2uiTAnHAAAAADLoRwCgEqUmVeopz7apn/tPOky16JeiOaNjFar+qEmJAMAAABgVZRDAFBJDpzK1LjFm7QvNdNlru/V4Xp9aEdVD/Q3IRkAAAAAK6McAoBK8PWu3/T4h1t1Pq/Qadxmk57se5Um9GopH64vBAAAAMAElEMAUIHsdkOzv92nv3+7z2WuRpC//n5PZ/VqU8+EZAAAAADwO8ohAKggGdkFmvjhFn2/55TLXFSDMM0d0VVN6wSbkAwAAAAA/otyCAAqQPLJcxq3eJMOn852mbuzU0PNGNRBQdV8TUgGAAAAAM4ohwDAzT7fdlxPL9uunIIip3FfH5v+fGuUHugeKZuN6wsBAAAA8AyUQwDgJoVFds1clax31h50masbWk1v3ttF3VrUMSEZAAAAAJSOcggA3OB0Zp4efn+L1h847TLXqUlNvTWiixrUCDIhGQAAAABcGuUQAPxBO45laNzijTqekesyd09sE027o50C/Li+EAAAAADPRDkEAH/Aqp0nNfHDLcotsDuNV/P10fN3ttM9sU1NSgYAAAAAZUM5BABXwDAMzV1zQDNXJcswnOciwgL11ogu6ty0ljnhAAAAAKAcKIcAoJzyC+36y/Kd+nDjUZe52Oa19Y97u6he9QATkgEAAABA+VEOAUA5ZGQXaPySTSVeeHpYdBP9deA1qubnY0IyAAAAALgylEMAUEaH0rI0elGSDpzKchq32aSn+7XVuBtayGazmZQOAAAAAK4M5RAAlMEvB8/ofxZvVHp2gdN4oL+PZg/rrH7XRJiUDAAAAAD+GMohALiMTzYf09Mfb1dBkfOVp+tXD1DC/dHq0LimOcEAAAAAwA0ohwCgFHa7of/3zV698d1+l7moBmGaf3+0GtYMMiEZAAAAALgP5RAAlCC3oEhPfLRNK7efcJnr07a+/n5PZ4UE8C0UAAAAgPfjNxsAuMip83ka++5GbT2a7jI3untz/fm2KPn6cOFpAAAAAFUD5RAAXGDPyfManZiklPQcp3FfH5um3dFOI7s1MykZAAAAAFQMyiEA+D8/7D2lP723WZl5hU7j1QP89OZ9XdSrTT2TkgEAAABAxaEcAgBJi9cf0rQvdqnI7nxHssa1grQgPkZtwqublAwAAAAAKhblEABLK7IbenHlLi388ZDLXOemNfXOqGjVDQ2o/GAAAAAAUEkohwBYVmZeoR79YIu+S051mRvQoYFeG9JRgf6+JiQDAAAAgMpDOQTAko6n52h0YpKST553mXu0dytNvKmNfLgjGQAAAAALoBwCYDnbjqbrwXc36tT5PKfxar4+mjm4ve7q3NikZAAAAABQ+SiHAFjKqp0nNPHDrcotsDuN1wr217xR0YqJrG1SMgAAAAAwB+UQAEswDENv/3BAM1clu8y1qBeihfExalYnxIRkAAAAAGAuyiEAVV5+oV3PLt+hpRuPucxd37KO3rqvq2oE+5uQDAAAAADMRzkEoEpLz87X+CWb9POBMy5zw2Oa6K8Dr5G/r48JyQAAAADAM1AOAaiyDqZlaUxikg6kZTmN22zSM/3bamzPFrLZuCMZAAAAAGujHAJQJW04cFrjlmxSenaB03igv49mD+usftdEmJQMAAAAADwL5RCAKmfZpmN65pPtKigynMbrVw/Q/Ptj1L5xDZOSAQAAAIDnoRwCUGXY7YZmfb1Xb36/32Xu6gZhmh8frQY1gkxIBgAAAACei3IIQJWQW1CkJ5Zu08odJ1zmboqqr78N76yQAL7lAQAAAMDF+E0JgNdLPZ+rse9u0raj6S5zY3o015Rbo+Trw4WnAQAAAKAklEMAvFryyXMak7hRKek5TuO+PjY9f0c7jejWzKRkAAAAAOAdKIcAeK3Ve1L18PtblJlX6DRePcBPc0Z0Uc/W9UxKBgAAAADeg3IIgFd6d/0hTfv8V9mdb0imxrWCtDA+Rq3Dq5sTDAAAAAC8DOUQAK9SWGTXiyt3K/GnQy5zXZrW1LxR0aobGlD5wQAAAADAS1EOAfAamXmFeuT9zfp+zymXuTs6NtQrgzso0N/XhGQAAAAA4L0ohwB4hZT0HI1JTFLyyfMuc4/1aa2JN7WWzcYdyQAAAACgvCiHAHi8rUfT9eCijUrLzHMar+bro1cGd9DAzo1MSgYAAAAA3o9yCIBH+3LHCU36cKvyCu1O47VDqmnuyK6KiaxtUjIAAAAAqBoohwB4JMMwNGf1f/Tqv/e4zLWsF6IF8TFqVifEhGQAAAAAULVQDgHwOPmFdk35dIeWbTrmMte9VR3Nua+ragT5m5AMAAAAAKoeyiEAHuVsVr7GL9mkDQfPuMzdE9tEL9x5jfx9fUxIBgAAAABVE+UQAI9xMC1LoxOTdDAty2ncZpOm9I/Sgz2bc0cyAAAAAHAzyiEAHuHnA6c1fskmpWcXOI0H+ftq9vBOuqVdhEnJAAAAAKBq49wML5CSkqLZs2erb9++atq0qapVq6aIiAjdfffd2rBhQ4nPWbJkicaNG6fo6GgFBATIZrMpMTGxcoMDZfTRxqMaOX+DSzEUHhagj8ZfRzEEAAAAABWII4e8wBtvvKGZM2eqZcuWuvnmm1W/fn3t27dPy5cv1/Lly/XBBx9o6NChTs959tlndfjwYdWtW1cNGjTQ4cOHTUoPlM5uN/TaV3s0Z/V/XObaNQzT/PtjFFEj0IRkAAAAAGAdlENeIDY2VmvWrFHPnj2dxteuXas+ffpowoQJuvPOOxUQEOCYS0hIUOvWrdWsWTPNmDFDzzzzTGXHBi4pJ79IT3y0VV/uOOkyd1NUuP42vJNCAvgWBQAAAAAVjd+8vMCgQYNKHO/Zs6fi4uL01VdfaceOHYqOjnbM3XTTTZUVDyi31PO5Grtoo7Ydy3CZe7BHcz1za5R8fbjwNAAAAABUBq45dBmpqalasWKFpk6dqv79+6tu3bqy2Wyy2WyKj48v12sdOXJETz75pKKiohQSEqLatWsrNjZWr732mrKzs68on7+/vyTJz4+eD95h94lzuusfP7kUQ74+Nr18V3s9O+BqiiEAAAAAqEQ0CpcRHh7ultdZuXKl7rvvPmVk/PcX4uzsbCUlJSkpKUkJCQn68ssv1aJFizK/5pEjR/TNN98oIiJC7du3d0tOoCJ9n5yqh9/frKz8Iqfx6oF+mnNfF/VsXc+kZAAAAABgXRw5VA5NmjRR3759y/28bdu2aejQocrIyFBoaKheeukl/fTTT/r22281duxYSdKePXt02223KTMzs0yvWVBQoJEjRyovL0+vvPKKfH19y50LqEyJPx7UmEVJLsVQk9pB+mTC9RRDAAAAAGASjhy6jKlTpyomJkYxMTEKDw/XoUOH1Lx583K9xsSJE5WdnS0/Pz999dVXuu666xxzvXv3VuvWrTV58mQlJydr1qxZmjp16iVfz263a/To0VqzZo3Gjh2rkSNHXtHXBlSGwiK7/rpilxatd71jXtdmtTRvZFfVCQ0o4ZkAAAAAgMrAkUOX8fzzz2vAgAFXfHpZUlKSVq9eLUkaM2aMUzFU7IknnlBUVJQkafbs2SooKCj19QzD0NixY7VkyRKNGDFCb7/99hXlAirD+dwCPfjuxhKLoTs7NdR7D15LMQQAAAAAJqMcqmDLly93PH7ggQdKXOPj46NRo0ZJks6ePesoky5mt9s1ZswYLViwQPfcc48SExPl48M/ITzTsbPZGvzWeq3ec8plbuJNrTV7WCcF+nM6JAAAAACYjWahgq1du1aSFBISoq5du5a6rlevXo7H69atc5m32+168MEHtXDhQg0bNkyLFy/mOkPwWFuOnNXAf/ykPb+ddxqv5uujvw3vpIk3tZHNxh3JAAAAAMATcM2hCrZ7925JUqtWrS55u/m2bdu6PKdY8RFDiYmJGjJkiJYsWVLpxdCxY8cuOX/ixIlKSgJPt3L7CT2+dKvyCu1O47VDqmneyK6KjqxtUjIAAAAAQEkohypQbm6u0tLSJEmNGze+5NpatWopJCREWVlZOnr0qNPcCy+8oMTERIWGhqpNmzZ68cUXXZ4/cOBAderUyfFxQkKC4wikHTt2OMaKT1kbOHCgBg4cWOavpUmTJmVeC2syDENzVv9Hr/57j8tcq/qhWnB/jJrWCTYhGQAAAADgUiiHKtD58/89pSY0NPSy64vLoYtvZ3/o0CFJUmZmpl566aUSnxsZGelUDq1bt06LFi1yWvPjjz/qxx9/dKwvTzkEXEpeYZGmfLJTH292PcKsR6u6+sd9XVQjyN+EZAAAAACAy6EcqkC5ubmOx9WqVbvs+oCA3+/alJOT4zSemJioxMTEcn3uK3nOpVx8NNPFTpw4odjYWLd9PniPs1n5Grdkk345eMZl7p7Ypnrhznby9+XyZgAAAADgqSiHKlBgYKDjcX5+/mXX5+XlSZKCgoIqLNOVutxpcbCmA6cyNToxSYdOZzuN22zSn2+N0pgezbnwNAAAAAB4OMqhClS9enXH44tPFStJVlaWpLKdggaYbf1/Tmv8kk3KyClwGg/y99Xf7+msm68ONykZAAAAAKA8KIcqUGBgoOrWrau0tLTL3u3r7NmzjnKIiz/D0y3deFRTPtmhQrvhNB4eFqD598fomkY1TEoGAAAAACgvLgRSwaKioiRJ+/fvV2FhYanrkpOTXZ4DeBq73dDMVcmavGy7SzHUrmGYPvtTD4ohAAAAAPAylEMVrEePHpJ+P2Vs06ZNpa774YcfHI+7d+9e4bmA8srJL9Kf3t+st1b/x2Xu5qvD9dH46xRRI7CEZwIAAAAAPBnlUAW78HbxCxcuLHGN3W7Xu+++K0mqWbOm4uLiKiMaUGap53I1bN56/WvnSZe5/7mhhd4e0VXB1ThLFQAAAAC8EeVQBYuNjVXPnj0lSfPnz9f69etd1rz++uvavXu3JOmxxx6Tv79/pWYELmXX8XMa+I8ftf1YhtO4r49NL9/VXlNujZKvD3ckAwAAAABvxX/1X8a6deu0f/9+x8dpaWmOx/v371diYqLT+vj4eJfX+Nvf/qbu3bsrJydHffv21ZQpUxQXF6ecnBz985//1Lx58yRJbdq00RNPPFEhXwdwJdbtS9P4JZuUmed8vazqgX56676u6tG6rknJAAAAAADuYjMMw7j8MuuKj4/XokWLyry+tL/OL774QiNGjNC5c+dKnG/Tpo1WrlypVq1aXVFOsx07dsxxl7WjR4+qcePGJifCH/X5tuN6YulWFRQ57+mmtYO1ID5arepXNykZAAAAAFhXRfz+zWllleT222/X9u3bNWnSJLVp00bBwcGqWbOmoqOjNXPmTG3ZssVriyFUPfPXHdSjH2xxKYaim9XSpw9dTzEEAAAAAFUIRw7BLThyqGoovlX93DUHXOZubR+hWUM7KdDf14RkAAAAAACpYn7/5ppDACRJBUV2Pb1suz7ZkuIyd/91zTT19nZceBoAAAAAqiDKIQDKyivUhPc2a83eUy5zT91ylR66saVsNoohAAAAAKiKKIcAi0vLzNPoxKQSb1U/Y1B7DYluYlIyAAAAAEBloBwCLOzI6WyNWrBBh05nO40H+vtozn1d1LttuEnJAAAAAACVhXIIsKidKRmKX5iktMw8p/Fawf5aEB+jzk1rmZQMAAAAAFCZKIcAC1q3L03jFm9UVn6R03ijmkF6d0ysWtYLNSkZAAAAAKCyUQ4BFvP5tuN6YulWFRQZTuNtI6pr0ehYhYcFmpQMAAAAAGAGyiHAQuavO6i/rtjlMt6tRW3NGxWtsEB/E1IBAAAAAMxEOQRYgN1uaOaqZM1dc8Bl7rb2DTRrWEcF+PmakAwAAAAAYDbKIaCKKyiy6+ll2/XJlhSXufuva6apt7eTr4/NhGQAAAAAAE9AOQRUYVl5hZrw3mat2XvKZe6pW67SQze2lM1GMQQAAAAAVkY5BFRRaZl5Gp2YpO3HMpzGfX1smjGovYZENzEpGQAAAADAk1AOAVXQkdPZGrVggw6dznYaD/L31Zz7uiiubX2TkgEAAAAAPA3lEFDF7EzJUPzCJKVl5jmN1wr214L4GHVuWsukZAAAAAAAT0Q5BFQh6/aladzijcrKL3Iab1QzSO+OiVXLeqEmJQMAAAAAeCrKIaCK+Gxrip78aJsKigyn8agGYUp8IEbhYYEmJQMAAAAAeDLKIaAKSFh7QC+u3O0y3q1Fbc0bFa2wQH8TUgEAAAAAvAHlEODF7HZDM1cla+6aAy5zt7VvoFnDOirAz9eEZAAAAAAAb0E5BHipgiK7nl62XZ9sSXGZi78+UlMHXC0fH5sJyQAAAAAA3oRyCPBCWXmFmvDeZq3Ze8pl7qlbrtJDN7aUzUYxBAAAAAC4PMohwMukZeZpdGKSth/LcBr39bFpxqD2GhLdxKRkAAAAAABvRDkEeJEjp7M1asEGHTqd7TQe5O+rOfd1UVzb+iYlAwAAAAB4K8ohwEvsTMlQ/MIkpWXmOY3XCvbXgvgYdW5ay6RkAAAAAABvRjkEeIF1+9I0bvFGZeUXOY03qhmkd8fEqmW9UJOSAQAAAAC8HeUQ4OE+25qiJz/apoIiw2k8qkGYEh+IUXhYoEnJAAAAAABVAeUQ4MES1h7Qiyt3u4x3a1Fb80ZFKyzQ34RUAAAAAICqhHII8EB2u6GZq5I1d80Bl7nb2jfQrGEdFeDna0IyAAAAAEBVQzkEeJiCIrsmL9uuT7ekuMzFXx+pqQOulo+PzYRkAAAAAICqiHII8CBZeYWa8N5mrdl7ymVucr+rNKFXS9lsFEMAAAAAAPehHAI8RFpmnkYnJmn7sQyncV8fm2YMaq8h0U1MSgYAAAAAqMoohwAPcOR0tkYt2KBDp7OdxoP8fTXnvi6Ka1vfpGQAAAAAgKqOcggw2c6UDMUvTFJaZp7TeK1gfy2Ij1HnprVMSgYAAAAAsALKIcBE6/aladzijcrKL3Iab1wrSItGx6plvVCTkgEAAAAArIJyCDDJZ1tT9ORH21RQZDiNRzUI06IHYlQ/LNCkZAAAAAAAK6EcAkyQsPaAXly522X8uhZ1NHdUV4UF+puQCgAAAABgRZRDQCWy2w3NXJWsuWsOuMzd1qGBZg3tqAA/XxOSAQAAAACsinIIqCQFRXZNXrZdn25JcZmLvz5SUwdcLR8fmwnJAAAAAABWRjkEVIKsvEJNeG+z1uw95TI3ud9VmtCrpWw2iiEAAAAAQOWjHAIqWFpmnkYnJmn7sQyncV8fm2YMaq8h0U1MSgYAAAAAAOUQUKGOnM7WqAUbdOh0ttN4kL+v5ozoorir6puUDAAAAACA31EOARVkZ0qG4hf+orTMfKfxWsH+WhAfo85Na5mUDAAAAACA/6IcAirAun1pGrd4o7Lyi5zGG9cK0qLRsWpZL9SkZAAAAAAAOKMcAtzss60pevKjbSooMpzGoxqEadEDMaofFmhSMgAAAAAAXFEOAW6UsPaAXly522X8uhZ1NHdUV4UF+puQCgAAAACA0lEOAW5gtxuauSpZc9cccJm7rUMDzRraUQF+viYkAwAAAADg0iiHgD+ooMiuycu269MtKS5z8ddHauqAq+XjYzMhGQAAAAAAl0c5BPwBWXmFmvDeZq3Ze8plbnK/qzShV0vZbBRDAAAAAADPRTkEXKG0zDyNTkzS9mMZTuO+PjbNGNReQ6KbmJQMAAAAAICyoxwCrsCR09katWCDDp3OdhoP8vfVnBFdFHdVfZOSAQAAAABQPpRDQDntTMlQ/MJflJaZ7zReK9hfC+Jj1LlpLZOSAQAAAABQfpRDQDms25emcYs3Kiu/yGm8ca0gvTs6Vi3qhZqUDAAAAACAK0M5BJTRZ1tT9ORH21RQZDiNRzUI06IHYlQ/LNCkZAAAAAAAXDnKIaAMEtYe0Isrd7uMX9eijuaO6qqwQH8TUgEAAAAA8MdRDgGXYLcbmrEqWfPWHHCZu61DA80a2lEBfr4mJAMAAAAAwD0oh4BSFBTZNXnZdn26JcVlLv76SE0dcLV8fGwmJAMAAAAAwH0oh4ASZOUVasJ7m7Vm7ymXuaf7tdX4Xi1ks1EMAQAAAAC8H+UQcJG0zDyNTkzS9mMZTuO+PjbNvLuDBndtbFIyAAAAAADcj3IIuMCR09katWCDDp3OdhoP8vfVnBFdFHdVfZOSAQAAAABQMSiHgP+zMyVD8Qt/UVpmvtN47ZBqWhAfo05NapoTDAAAAACACkQ5BEhaty9N4xZvVFZ+kdN441pBend0rFrUCzUpGQAAAAAAFYtyCJb32dYUPfnRNhUUGU7jUQ3CtOiBGNUPCzQpGQAAAAAAFY9yCJaWsPaAXly522X8uhZ1NHdUV4UF+puQCgAAAACAykM5BEuy2w3NWJWseWsOuMzd1qGBZg3tqAA/XxOSAQAAAABQuSiHYDl2u6EnPtqmT7ekuMzFXx+pqQOulo+PzYRkAAAAAABUPsohWI6Pj031qwe4jD/dr63G92ohm41iCAAAAABgHZRDsKSn+7VV6vk8fbolRb4+Ns28u4MGd21sdiwAAAAAACod5RAsyef/CqGc/CINi22iuKvqmx0JAAAAAABTUA7Bsqr5+ejtkV3NjgEAAAAAgKl8zA4AAAAAAAAA81AOAQAAAAAAWBjlEAAAAAAAgIVRDgEAAAAAAFgY5RAAAAAAAICFUQ4BAAAAAABYGOUQAAAAAACAhVEOAQAAAAAAWBjlEAAAAAAAgIVRDgEAAAAAAFgY5RAAAAAAAICFUQ4BAAAAAABYGOUQAAAAAACAhVEOAQAAAAAAWBjlEAAAAAAAgIVRDgEAAAAAAFgY5RAAAAAAAICFUQ4BAAAAAABYGOUQAAAAAACAhVEOAQAAAAAAWBjlEAAAAAAAgIVRDgEAAAAAAFgY5RAAAAAAAICFUQ4BAAAAAABYGOUQAAAAAACAhVEOAQAAAAAAWBjlEAAAAAAAgIVRDgEAAAAAAFiYn9kBUDUUFhY6Hp84ccLEJAAAAAAAVF0X/s594e/ifwTlENzi1KlTjsexsbEmJgEAAAAAwBpOnTqlyMjIP/w6nFYGAAAAAABgYTbDMAyzQ8D75ebmaseOHZKkevXqyc+Pg9K83YkTJxxHgf3yyy9q0KCByYngrdhLcBf2EtyFvQR3YB/BXdhLKK/CwkLH2Tvt27dXYGDgH35NfoOHWwQGBiomJsbsGKggDRo0UOPGjc2OgSqAvQR3YS/BXdhLcAf2EdyFvYSycsepZBfitDIAAAAAAAALoxwCAAAAAACwMMohAAAAAAAAC6McAgAAAAAAsDDKIQAAAAAAAAujHAIAAAAAALAwyiEAAAAAAAALsxmGYZgdAgAAAAAAAObgyCEAAAAAAAALoxwCAAAAAACwMMohAAAAAAAAC6McAgAAAAAAsDDKIQAAAAAAAAujHAIAAAAAALAwyiEAAAAAAAALoxwCUCEiIyNls9lK/DN+/Hiz48FLJCYmlrqPiv/06dPH7JjwEna7XW+++aa6dOmi4OBghYWFqVevXvr888/NjgYPtGTJEo0bN07R0dEKCAiQzWZTYmKi29bDOsqzN9LT0/Xoo4/quuuuU0REhAICAtSoUSP17t1bH3/8sQzDqNzw8Cjl/T7D+3GUh5/ZAQBUXTVq1NDEiRNdxqOjoys/DLxSp06d9Nxzz5U4t2zZMv3666+65ZZbKjkVvJFhGBo6dKg+/vhjtWzZUmPGjFFeXp4+++wz3XnnnXrjjTf08MMPmx0THuTZZ5/V4cOHVbduXTVo0ECHDx9263pYR3n2RlpamhYsWKBu3bpp4MCBql27tlJTU/XFF19o8ODBGjt2rObNm1eJ6eFJruT7DO/HUVY2g/oZQAWIjIyUJB06dMjUHKia8vPz1bBhQ2VkZOjYsWMKDw83OxI83LJlyzRkyBB1795dX3/9tYKCgiT9/otYdHS0Tp48qeTkZMf3LuCbb75R69at1axZM82YMUPPPPOMFi5cqPj4eLesh3WUZ28UFRXJMAz5+Tn/H/758+fVrVs37dq1Szt37lS7du0qKT08SXm/z/B+HOXBaWUAAK/z6aef6vTp0xowYADFEMpk+fLlkqQpU6Y4iiFJqlu3riZNmqS8vDwtXLjQpHTwRDfddJOaNWtWYethHeXZG76+vi7FkCRVr17dcaTs/v373ZoP3oPvM6hIlENAFZOamqoVK1Zo6tSp6t+/v+rWres4t7i8/3t55MgRPfnkk4qKilJISIhq166t2NhYvfbaa8rOzr7s8/Py8rRo0SK9/PLLeuutt7Rt27Yr/KpgBk/aSxebP3++JOnBBx8s93NR+TxhL/3222+SpObNm7vMFY9999135cqCiuEJ+wVVQ1XaS7m5ufruu+9ks9l09dVXV/jngzNv3ku8H0dZcc0hoIpx11EUK1eu1H333aeMjAzHWHZ2tpKSkpSUlKSEhAR9+eWXatGiRamvcfLkSZcfmP369dPixYtVt25dt+RExfGkvXShw4cP69tvv1WjRo3Ur18/t2RExfKEvVSvXj1J0sGDBxUVFeU0d/DgQUnS3r173ZITf4wn7BdUDd68l9LT0zV79mzZ7Xalpqbqyy+/1NGjR/Xcc8+pdevWbvs8KBtv3ku8H0dZceQQUIU1adJEffv2Lffztm3bpqFDhyojI0OhoaF66aWX9NNPP+nbb7/V2LFjJUl79uzRbbfdpszMzBJfY/To0Vq9erVOnTqlc+fO6eeff1b//v21atUq3XHHHdxtw8uYuZcutnDhQtntdj3wwAPy9fUtdyaYy6y91L9/f0nSjBkzlJub6xg/ffq0Zs+eLen3X8bgWTzpew+8m7ftpfT0dD3//PP661//qrlz5+rkyZN69dVXS71JAyqPN+0l3o+jXAwAVcrUqVONL774wjh58qRhGIZx8OBBQ5Ihybj//vvL9Bo33nijIcnw8/MzfvrpJ5f5V155xfGazz//fJmzFRUVGT169DAkGStWrCjz82AOT9xLRUVFRtOmTQ2bzWYcOHCgXF8PzOMJe6mgoMCIi4szJBmtWrUyHn74YWPcuHFGeHi40aFDB0OSERQU9Ie+TriHJ+yXi02fPt2QZCxcuLBMn7+861ExqsJeKiwsNA4ePGhMnz7dqFatmnHXXXcZBQUFZXou3Kcq7KVivB9HaSiHgCquvD+8fvnlF8f6cePGlbimqKjIiIqKMiQZtWrVMvLz88ucZ/78+YYk45lnninzc+AZPGEvrVq1ypBk9OnT50q+BHgIs/ZSbm6uMW3aNKNNmzZGtWrVjHr16hn/8z//Y+zdu9eQZDRt2vSPfmmoAJ7wvYdyqGrwxr10oeLyYM6cOeV+LtzL2/cS78dREk4rA+Ck+I4+kvTAAw+UuMbHx0ejRo2SJJ09e1arV68u8+sXn9vMhUCrvorYS1yI2prctZcCAgL03HPPac+ePcrLy1Nqaqrmzp2rlJQUSVJ0dLTbs6PyVfTPMViHp+2l4lOZ2K/ex9P2Eu/HURLKIQBO1q5dK0kKCQlR165dS13Xq1cvx+N169aV+fU3bNggSYqMjLyygPAa7t5Lp0+f1meffabatWvrrrvucl9QeLyK/r703nvvSZKGDx9+hQnhSSp6v8A6PG0vHT9+XJJKvNU9PJun7SXej6MklEMAnOzevVuS1KpVq0u++Wjbtq3Lc4rt2rWrxAu7rlu3TrNmzVJAQIAGDRrknsDwWO7YSxdavHix8vPzNWLECAUEBLgvKDyeu/bSuXPnXMaWLVumBQsWKCYmhu9LVYS7v/fAuszYS1u3bnW6k1WxM2fOaMqUKZL+e4F9eA8z9hLvx1Fe1M4AHHJzc5WWliZJaty48SXX1qpVSyEhIcrKytLRo0ed5pYuXapXXnlFffr0UWRkpAICArRz50599dVX8vHx0dtvv62mTZtW2NcB87lrL12IU8qsyZ176dprr1WTJk0UFRWlwMBA/fLLL1q9erVatGihjz76iLvfVQHu3C8JCQmO/7nfsWOHY6z4VI+BAwdq4MCBV7wens2svZSYmKiEhATFxcWpWbNmCgkJ0eHDh7Vy5UplZmbq7rvv1r333uumrxKVway9xPtxlBflEACH8+fPOx6HhoZedn3xD6+Lb7cZFxen3bt3a/Pmzfrhhx+Um5ur8PBwDRs2TJMmTVJsbKzbs8OzuGsvFfvll1+0c+dOxcbGqn379m7LCc/nzr00bNgwffLJJ/r5559VUFCg5s2b69lnn9VTTz2lsLAwt+aGOdy5X9atW6dFixY5jf3444/68ccfJf1+OsaFZU9518OzmbWXBg8erIyMDP38889as2aNsrOzVbt2bfXo0UOjRo3S8OHDZbPZ/sBXhspm1l7i/TjKi3IIgENubq7jcbVq1S67vvjUnpycHKfxXr16OZ0zDetx114qFhsbK8Mw3BMOXsWde2natGmaNm2a27LB87hzvyQmJioxMbHMn7u86+HZzNpLPXr0UI8ePcoWEl7BrL3E+3GUF9ccAuAQGBjoeJyfn3/Z9Xl5eZKkoKCgCssE78Regruwl1Ae7Be4C3sJ7sJegregHALgUL16dcfj0k7vuVBWVpaksh0iC2thL8Fd2EsoD/YL3IW9BHdhL8FbUA4BcAgMDFTdunUlSceOHbvk2rNnzzp+eDVp0qTCs8G7sJfgLuwllAf7Be7CXoK7sJfgLSiHADiJioqSJO3fv1+FhYWlrktOTnZ5DnAh9hLchb2E8mC/wF3YS3AX9hK8AeUQACfFF0HMysrSpk2bSl33ww8/OB537969wnPB+7CX4C7sJZQH+wXuwl6Cu7CX4A0ohwA4ufA2uwsXLixxjd1u17vvvitJqlmzpuLi4iojGrwMewnuwl5CebBf4C7sJbgLewnegHIIgJPY2Fj17NlTkjR//nytX7/eZc3rr7+u3bt3S5Iee+wx+fv7V2pGeAf2EtyFvYTyYL/AXdhLcBf2EryBzTAMw+wQANxn3bp12r9/v+PjtLQ0PfXUU5J+Pzz1wQcfdFofHx/v8hpbtmxR9+7dlZOTo9DQUE2ZMkVxcXHKycnRP//5T82bN0+S1KZNG23cuNHpLgyoOthLcBf2EsqD/QJ3YS/BXdhLsAQDQJVy//33G5LK/Kc0n3/+uREWFlbq89q0aWPs27evEr8yVDb2EtyFvYTyYL/AXdhLcBf2EqyA08oAlOj222/X9u3bNWnSJLVp00bBwcGqWbOmoqOjNXPmTG3ZskWtWrUyOya8AHsJ7sJeQnmwX+Au7CW4C3sJnozTygAAAAAAACyMI4cAAAAAAAAsjHIIAAAAAADAwiiHAAAAAAAALIxyCAAAAAAAwMIohwAAAAAAACyMcggAAAAAAMDCKIcAAAAAAAAsjHIIAAAAAADAwiiHAAAAAAAALIxyCAAAAAAAwMIohwAAAAAAACyMcggAAAAAAMDCKIcAAAAAAAAsjHIIAAAAAADAwiiHAAAAAAAALIxyCAAAAAAAwMIohwAAAAAAACyMcggAAAAAAMDCKIcAAAAAAAAsjHIIAAAAVUZ6eroeffRRXXfddYqIiFBAQIAaNWqk3r176+OPP5ZhGGZHBADA49gMfkICAACgiti/f786deqkbt26qVWrVqpdu7ZSU1P1xRdfKDU1VWPHjtW8efPMjgkAgEehHAIAAECVUVRUJMMw5Ofn5zR+/vx5devWTbt27dLOnTvVrl07kxICAOB5OK0MAADAg8THx8tms7n8OXTokNnRvIKvr69LMSRJ1atX1y233CLp96OLLrZ69eoS/96nTZtW0ZEBADAd5RAAAAAcUlJSHMVIUlKS2XHcJjc3V999951sNpuuvvpqs+MAAOBRXP9bBQAAAKZr2LCh/v3vfzs+btSoUaV83hUrVkiSIiIiFB0dXSmfsyKkp6dr9uzZstvtSk1N1ZdffqmjR4/queeeU+vWrV3Wx8TEaMeOHY6P27dvX5lxAQAwFeUQAACAB/L399c111xT6Z/3iy++kCQNGDBANput0j+/u6Snp+v55593fOzv769XX31VTzzxRInrQ0JCTPn7BgDAE3BaGQAAACRJ2dnZ+u677yRJt99+u8lp/pjIyEgZhqHCwkIdPHhQL7zwgv785z/r7rvvVmFhodnxAADwKJRDAAAAkCR98803ysnJUWBgoG666Saz47iFr6+vIiMj9b//+7968cUX9emnn+qdd94xOxYAAB6FcggAAFjO3r179cQTT6hLly6qUaOG/P39Va9ePbVr10733HOPZs2aZXbEcjEMQ++//75uueUW1atXTyEhIerYsaNmz56twsJC5eTkyMfHRzabTTNnziz1dYqvN9S7d28FBweXui41NVXTp09XXFycGjRooICAADVs2FDXX3+9pk2bpt27d7s8Z//+/Y4LXX/wwQeSpKVLl6pfv34KDw9XaGioOnfurAULFsgwDMfz8vLyNH/+fN1www2qW7eugoKCdO211+qzzz4r999T3759Jf1+ZzIAAPBfXHMIAABYyvTp0/WXv/xFRUVFTuNpaWlKS0vTrl27tGHDBj3++OMmJSyf1NRUDRw4UOvXr3ca3759uyZNmqSvvvpKL730kqNw6dixY4mvYxiGVq5cKenSp5S9+uqreu6555STk+M0fuLECZ04cULr16/XqlWr9PPPPzvNb9261fG4cePG6tevn9MFt4vXjBkzRseOHdPUqVO1d+9eDRs2zOm5kvTLL79o4MCBWrp0qYYMGVJq1osdP35ckkq81T0AAFbGT0YAAGAZc+fO1ZQpUyRJnTp10ujRo9WhQwdVr15dp06d0tGjR7VhwwbVqlXL5KRlk5mZqRtvvFG7d++WzWbTPffco+HDh6tRo0Y6ePCgXn75Zf3rX/9Sbm6u4zmdOnUq8bU2bdrkKE8GDBhQ4pqxY8cqISFBktSkSRONHTtW3bt3V61atfTbb79pzZo1SkxMVNeuXV2eu23bNsfjJ598Utu3b9fEiRN1++23Kzg42FFi5efn6+WXX9bgwYPVp08f5efn64UXXlDv3r1VVFSkDz/8UHPmzJEk/eUvf3Eph7Zu3armzZurRo0aTuNnzpxx/Nv379//Un+tAABYDuUQAACwjJdeekmSdPfdd2vp0qXy8XE9w/7BBx+s7FhXbPz48dq9e7f8/Pz0ySefOB3x06VLF/Xr109t27bV999/L0mqX7++IiIiSnyt4ruUderUSY0bN3aZf/nllx3F0MiRI/XOO+8oICDAaU2/fv307LPP6siRIy7Pv/Don/3792v9+vVORVW3bt2Ul5enl19+WXl5eerevbsiIiL0zTffqFGjRo51N9xwg06cOKFPP/1Ue/bsUVpamurWreuYT0xMVEJCguLi4tSsWTOFhITo8OHDWrlypTIzM3X33Xfr3nvvLe2vFAAAS+KaQwAAwBLOnTuno0ePSpJuueWWEoshb7J27Vq99957kqRp06aVeCpYSEiIHnvsMcfHpR01JP33ekMlvc6vv/6qqVOnSvq9AEpMTHQphooFBwerbdu2LuMXHjmUkJBQYpY+ffo4Hufm5mrp0qVOxVCxCy+WnZGR4TQ3ePBgDRkyRPv379fixYs1a9Ysff/99+rRo4fef/99ffTRR17/bw8AgLtx5BAAALCE0NBQ1alTR6dPn9YLL7yg2rVrq1evXqpTp45sNpvZ8crthRdekCQ1bdpUTz/9dKnr2rVr53hc2vWGUlJStHnzZkkll0OTJ09WUVGRQkNDtWjRonKXK2fOnHEq5u66664S1/n6+joe/+lPf1L79u1LXGe32x2PLz59rEePHurRo0e58gEAYHX8twkAALAEHx8fJSYmKjg4WMeOHdPgwYNVr149x128lixZYnbEMjtx4oS++eYbSdJDDz10yQssX1ielHbkUPFRQxEREYqOjnaaO3nypP71r39JkiZMmKD69euXO++FRw0NGzas1HX79u0r07r9+/dLksLCwpxOKQMAAFeGcggAAFhG37599fbbb5dYKHTo0MGERFemuKyRpFtvvfWSa8+cOeN4XNqRQ8XXG7rttttcjqJasWKF405n5bkz2IUuvN7QpfIWl0j16tVTTEzMZdd5078ZAACejNPKAACAJWzevFkjR47Url271KdPH91///1q3ry5atasKUmKiooyN2A5FJctAQEBpZ56VWznzp2SpMDAwBKvBZSTk6PvvvtOUsmnlO3YsUOS5O/v73JUUVkVlzkNGzZUeHh4qeu2bNki6dLXRrrw9S63DgAAlA3lEAAAqPIOHz6sPn36KD09XW+88YYefvhhsyP9ISdPnpT0+xE2l/P1119L+v3aQxde06fYN998o5ycHAUGBjpd6PnizxUeHn7F12YqLrM6d+5c6hrDMMpU+hw5ckRnz5697DoAAFB2nFYGAACqvOnTpys9PV3XXnut1xdDkpSXlyfp9zuwXUpycrLjNvalFSnFp5T17t1bISEhLvO5ubmSpPz8/CvKWlBQoN27d18yg/T7dYQyMzMlXbpEuvAUNcohAADcg3IIAABUeWvWrJEktWrVyuQk7lF8Uehz58457gJ2MbvdrkceecRxvaCSrjdkGIZWrlwpqeRTyqTfTwWTpNTUVB05cqTcWXft2uUoli5V+hSfUiZduvQpLof8/Px0zTXXlDsPAABwRTkEAACqvOzsbEnS999/r7S0tFLX5ebmOo7K8WTdunVzPJ4+fbrLfFFRkR555BHHHc2kkguXTZs26fjx45KkAQMGlPi54uLiHI///Oc/O8qmixUUFDjdbazYhXcqK8sRQUFBQWrTpk2p64pfr23btgoICCh1HQAAKDuuOQQAAKq866+/XocPH9bx48d17bXX6tFHH1V0dLRCQkL022+/6fjx4/rhhx+0fPlyHT161ONLhyFDhuipp57S2bNn9dZbb6mgoED33nuvQkNDtWvXLr355pvauHGjmjZt6jjap6Q7exXfwr5Tp05q3LhxiZ9r0KBBatu2rZKTk7VkyRIdO3ZMY8eO1VVXXSXDMHTo0CGtWbNGy5Yt0+TJkzVx4kSn5xeXPjVq1FDz5s1L/ZqKjxzq0KFDiddGuvj1OKUMAAD3oRwCAABV3vTp07V+/XodOnRIBw4ccCkwirVu3VrVq1ev3HBXICwsTAkJCRo2bJgKCwuVkJCghIQEx7yvr6+mTJmi7OxszZ49W1FRUapRo4bL6xRfb6i0U8qk30/f+uyzz9S/f38dOHBAq1ev1urVq0tcW1Jhc+FFpi91QeuylD7nzp3TwYMHL7sOAACUD+UQAACo8po1a6YdO3Zo3rx5WrFihX799VedPXtW1apVU3h4uFq0aKEePXqof//+Zkcts0GDBmnNmjV68cUX9dNPPyk3N1eNGjVSnz599NBDD6ljx47q2rWrJOmGG25weX5KSorjaJ1LlUOS1KZNG23dulXz5s3T8uXL9euvv+r8+fOqU6eOGjZsqG7duun2229Xjx49XJ5bljuQ/fbbb467ol1q3fbt2x2ntVEOAQDgPjajtBPHAQAAUOni4+O1aNEiNWvWTIcOHbri19m5c6fat28v6ffb1ffp08dpfu7cuRo/frwiIiJ0/PjxK75NfVVV/Pfx3HPPadq0aeaGAQCggnHkEAAAgAcqKCjQzp07HR9fddVV8vf3L9Nz7Xa7Jk2aJOn3o3569+7tsqb4ekO33XYbxZCkrKwsxylrAABYDeUQAACABzp+/LjjyB9JOnjwoCIjIyVJ//nPf9SyZcsSn5ebm6uHHnrIcaey1157rcTyp2fPnuratavuvPNO94f3QklJSU53ZgMAwEoohwAAALzMkCFDFBAQoOHDh6tTp06qUaOGzp49q59//lnvvPOO4wiYxx9/vNTrCU2ePLkyIwMAAA/GNYcAAAC8SGFhoUJDQ5WXl1fqGj8/P02bNk1TpkzhlDEAAHBZlEMAAABepLCwUMuWLdPnn3+uLVu26NSpU8rIyFBYWJhatWqlPn36aPz48WratKnZUQEAgJegHAIAAAAAALAwH7MDAAAAAAAAwDyUQwAAAAAAABZGOQQAAAAAAGBhlEMAAAAAAAAWRjkEAAAAAABgYZRDAAAAAAAAFkY5BAAAAAAAYGGUQwAAAAAAABZGOQQAAAAAAGBhlEMAAAAAAAAWRjkEAAAAAABgYZRDAAAAAAAAFkY5BAAAAAAAYGGUQwAAAAAAABZGOQQAAAAAAGBhlEMAAAAAAAAWRjkEAAAAAABgYZRDAAAAAAAAFkY5BAAAAAAAYGGUQwAAAAAAABZGOQQAAAAAAGBh/x95MQdU8vZ5DQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1280x960 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make sure we dont mess things up\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "cs_a = CS_EOS.gen_a((0.2, 0.2, 0.3, 0.4, 0.5))\n",
    "print(CS_EOS.check_a(cs_a))\n",
    "core_p = CS_EOS.cal_core_p(core_e, cs_a)\n",
    "\n",
    "test_e = np.concatenate((eps_com, core_e[1:]))\n",
    "test_p = np.concatenate((pres_com, core_p[1:]))\n",
    "plt.figure(dpi=200)\n",
    "plt.plot(test_e / g_cm_3, test_p / dyn_cm_2)\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(r\"$\\epsilon\\quad [g/cm^3]$\")\n",
    "plt.ylabel(r\"$P\\quad [dyn/cm^2]$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7d8d8d0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood_transform(paras):\n",
    "    a1, a2, a3, a4, a5, d1 = paras\n",
    "    a6 = CS_EOS.cal_a6(a1, a2, a3, a4, a5)\n",
    "    a = (a1, a2, a3, a4, a5, a6)\n",
    "    if not CS_EOS.check_a(a):\n",
    "        return -1e100\n",
    "    core_p = CS_EOS.cal_core_p(core_e, a)\n",
    "    full_p = np.concatenate((pres_com, core_p[1:]))\n",
    "\n",
    "    probMRgaussian = likelihood.MRlikihood_Gaussian(\n",
    "        full_e, full_p, (1.4, 13, 0.07, 0.65), d1\n",
    "    )\n",
    "\n",
    "    return probMRgaussian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a3802c4",
   "metadata": {},
   "source": [
    "### Set up sampler"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76bce7fb",
   "metadata": {},
   "source": [
    "Here next, we define sampler, there is two different sampler we provided for. \n",
    "\n",
    "Considering where you need resume file:\n",
    "\n",
    "sampler.UltranestSampler   and  sampler.UltranestSamplerResume\n",
    "\n",
    "Here since it is our first run, so we only use first one. Some of the sampler parameters is requested, first is step number, our choice for UltraNest sampler is slicesampler, which could easily be sliced up your total computation load, and parallelize, speed up sampling. So step as suggested by documentation of UltraNest, we use 2*len(parameters).\n",
    "\n",
    "live_point we set 2000, it will influence the sampling precision, We suggest for 7 dimension space, maybe 5000 is a better choice, however, since my computer only have limited resources, we set 2000.\n",
    "\n",
    "max_calls set 10000, it is how many iteration after it will stop, we suggest to set this number significantly higher, otherwise maybe will broken before the inference converging to a definite value. That result will be un-phyiscal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f05849a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating directory for new run output\\run11\n",
      "[ultranest] Sampling 1000 live points from prior ...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\mambaforge\\lib\\site-packages\\scipy\\integrate\\_ode.py:431: UserWarning: dopri5: larger nsteps is needed\n",
      "  self._y, self.t = mth(self.f, self.jac or (lambda: None),\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\InferenceWorkflow\\Likelihood.py:116: RuntimeWarning: divide by zero encountered in log\n",
      "  likelihood = np.log(fx)\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:79: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = -(eps + pres) * (m + 4.0 * pi * r**3.0 * pres)\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:80: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = dpdr / (r * (r - 2.0 * m))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ultranest] Widening roots to 1018 live points (have 1000 already) ...\n",
      "[ultranest] Sampling 18 live points from prior ...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9014ce033a1b4ab3bdf09bcc33be494e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HTML(value=''), GridspecLayout(children=(HTML(value=\"<div style='background-color:#6E6BF4;'>&nb…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-1e+100(0.00%) | Like=-1e+100..1.2 [-1e+101..-162.1] | it/evals=418/1049 eff=3.2258% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\InferenceWorkflow\\Likelihood.py:116: RuntimeWarning: divide by zero encountered in log\n",
      "  likelihood = np.log(fx)\n",
      "d:\\mambaforge\\lib\\site-packages\\scipy\\integrate\\_ode.py:431: UserWarning: dopri5: larger nsteps is needed\n",
      "  self._y, self.t = mth(self.f, self.jac or (lambda: None),\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-623.4(0.00%) | Like=-615.03..1.23 [-1e+101..-162.1] | it/evals=428/1261 eff=4.5267% N=601 85856815104.0(0.00%) | Like=-731.18..1.23 [-1e+101..-162.1] | it/evals=419/1069 eff=3.9216% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:79: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = -(eps + pres) * (m + 4.0 * pi * r**3.0 * pres)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-233.6(0.00%) | Like=-219.28..1.23 [-1e+101..-162.1] | it/evals=477/2136 eff=5.3667% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:80: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = dpdr / (r * (r - 2.0 * m))\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:80: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  dpdr = dpdr / (r * (r - 2.0 * m))\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\InferenceWorkflow\\Likelihood.py:115: RuntimeWarning: overflow encountered in power\n",
      "  fx = 1/(sigma_x*sigma_y*(np.sqrt(2*np.pi))**2)*np.exp(-np.power(MR[1]/km-Rvalue, 2.)/(2*np.power(sigma_x,2.))-np.power(MR[0]/Msun-Mvalue, 2.)/(2*np.power(sigma_y,2.)))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-158.6(0.00%) | Like=-152.26..1.23 [-161.8854..-57.6311] | it/evals=553/3469 eff=5.5488% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\mambaforge\\lib\\site-packages\\scipy\\integrate\\_ode.py:431: UserWarning: dopri5: larger nsteps is needed\n",
      "  self._y, self.t = mth(self.f, self.jac or (lambda: None),\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\InferenceWorkflow\\Likelihood.py:116: RuntimeWarning: divide by zero encountered in log\n",
      "  likelihood = np.log(fx)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-144.8(0.00%) | Like=-138.87..1.23 [-161.8854..-57.6311] | it/evals=591/4079 eff=5.6844% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:79: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = -(eps + pres) * (m + 4.0 * pi * r**3.0 * pres)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-142.6(0.00%) | Like=-136.39..1.23 [-161.8854..-57.6311] | it/evals=599/4198 eff=5.7233% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:80: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = dpdr / (r * (r - 2.0 * m))\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:80: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  dpdr = dpdr / (r * (r - 2.0 * m))\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\InferenceWorkflow\\Likelihood.py:115: RuntimeWarning: overflow encountered in power\n",
      "  fx = 1/(sigma_x*sigma_y*(np.sqrt(2*np.pi))**2)*np.exp(-np.power(MR[1]/km-Rvalue, 2.)/(2*np.power(sigma_x,2.))-np.power(MR[0]/Msun-Mvalue, 2.)/(2*np.power(sigma_y,2.)))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-112.4(0.00%) | Like=-105.42..1.23 [-161.8854..-57.6311] | it/evals=689/5673 eff=5.8432% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\mambaforge\\lib\\site-packages\\scipy\\integrate\\_ode.py:431: UserWarning: dopri5: larger nsteps is needed\n",
      "  self._y, self.t = mth(self.f, self.jac or (lambda: None),\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-111.9(0.00%) | Like=-105.32..1.23 [-161.8854..-57.6311] | it/evals=690/5696 eff=5.8358% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\InferenceWorkflow\\Likelihood.py:116: RuntimeWarning: divide by zero encountered in log\n",
      "  likelihood = np.log(fx)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-111.2(0.00%) | Like=-105.19..1.23 [-161.8854..-57.6311] | it/evals=692/5734 eff=5.8312% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:80: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = dpdr / (r * (r - 2.0 * m))\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:79: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = -(eps + pres) * (m + 4.0 * pi * r**3.0 * pres)\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:80: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  dpdr = dpdr / (r * (r - 2.0 * m))\n",
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\InferenceWorkflow\\Likelihood.py:115: RuntimeWarning: overflow encountered in power\n",
      "  fx = 1/(sigma_x*sigma_y*(np.sqrt(2*np.pi))**2)*np.exp(-np.power(MR[1]/km-Rvalue, 2.)/(2*np.power(sigma_x,2.))-np.power(MR[0]/Msun-Mvalue, 2.)/(2*np.power(sigma_y,2.)))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-79.5(0.00%) | Like=-73.21..1.23 [-161.8854..-57.6311] | it/evals=823/8225 eff=5.6334% N=601   \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\mambaforge\\lib\\site-packages\\scipy\\integrate\\_ode.py:431: UserWarning: dopri5: larger nsteps is needed\n",
      "  self._y, self.t = mth(self.f, self.jac or (lambda: None),\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-79.3(0.00%) | Like=-73.07..1.23 [-161.8854..-57.6311] | it/evals=824/8244 eff=5.6324% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\InferenceWorkflow\\Likelihood.py:116: RuntimeWarning: divide by zero encountered in log\n",
      "  likelihood = np.log(fx)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-78.5(0.00%) | Like=-72.40..1.23 [-161.8854..-57.6311] | it/evals=828/8324 eff=5.6255% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:79: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = -(eps + pres) * (m + 4.0 * pi * r**3.0 * pres)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=-70.3(0.00%) | Like=-64.40..1.23 [-161.8854..-57.6311] | it/evals=882/9389 eff=5.5549% N=601 \r"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\51970\\Desktop\\GitHub\\CompactOject\\TOVsolver\\solver_code.py:80: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  dpdr = dpdr / (r * (r - 2.0 * m))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ultranest] Explored until L=1  23 [-161.8854..-57.6311] | it/evals=915/10016 eff=5.5346% N=601 \n",
      "[ultranest] Likelihood function evaluations: 10016\n",
      "[ultranest] Writing samples and results to disk ...\n",
      "[ultranest] Writing samples and results to disk ... done\n",
      "[ultranest] Reached maximum number of likelihood calls (10016 > 10000)...\n",
      "[ultranest] done iterating.\n"
     ]
    }
   ],
   "source": [
    "step = 2 * len(parameters)\n",
    "live_point = 1000\n",
    "\n",
    "max_calls = 10000\n",
    "samples = sampler.UltranestSampler(\n",
    "    parameters, likelihood_transform, prior_transform, step, live_point, max_calls\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.10.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}