{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4247d008-55ad-4811-ac4a-9db28e40abc9",
   "metadata": {},
   "source": [
    "# Example of likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9f79d871-874a-425e-83c9-ab33a4515d8e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using Normal as an approximation of Binomial\n",
      "Using aptext package from https://github.com/XENONnT/applefiles\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import jax.numpy as jnp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import appletree as apt\n",
    "from appletree.utils import get_file_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f28fac9a-8714-4022-99a7-c5a6363e1124",
   "metadata": {},
   "outputs": [],
   "source": [
    "apt.set_gpu_memory_usage(0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a374360-de4f-49dc-83ea-e0147f38c4bd",
   "metadata": {},
   "source": [
    "## Likelihood configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d91318c8-8130-4cfe-99e0-ecf1a24de3d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "config = dict(\n",
    "    data_file_name=get_file_path(\"data_Rn220.csv\"),\n",
    "    bins_type=\"equiprob\",\n",
    "    bins_on=[\"cs1\", \"cs2\"],\n",
    "    bins=[15, 15],\n",
    "    x_clip=[0, 100],\n",
    "    y_clip=[2e2, 1e4],\n",
    ")\n",
    "\n",
    "llh = apt.Likelihood(**config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb5202a-1e75-451e-8869-24fc36713428",
   "metadata": {},
   "source": [
    "## Register components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "44cefc08-01ac-4d5b-8250-dc7d4fb4a2e7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "----------------------------------------\n",
      "BINNING\n",
      "\n",
      "    bins_type: equiprob\n",
      "    bins_on: ['cs1', 'cs2']\n",
      "\n",
      "----------------------------------------\n",
      "DATA\n",
      "\n",
      "    file_name: /home/xudc/appletree/appletree/data/data_Rn220.csv\n",
      "    data_rate: 2000.0\n",
      "\n",
      "----------------------------------------\n",
      "MODEL\n",
      "\n",
      "    COMPONENT 0: rn220_er\n",
      "        type: simulation\n",
      "        rate_par: rn220_er_rate\n",
      "        pars: {'py1', 'rf1', 'drift_velocity', 'p_dpe', 'rf0', 'py4', 'py0', 'nex_ni_ratio', 'g2', 'g1', 's2_cut_acc_sigma', 'w', 's1_cut_acc_sigma', 's2_threshold', 's1_eff_3f_sigma', 'field', 'rn220_er_rate', 'py2', 'gas_gain', 'fano', 'py3', 'elife_sigma'}\n",
      "\n",
      "    COMPONENT 1: rn220_ac\n",
      "        type: fixed\n",
      "        file_name: AC_Rn220.pkl\n",
      "        rate_par: rn220_ac_rate\n",
      "        pars: {'rn220_ac_rate'}\n",
      "\n",
      "----------------------------------------\n"
     ]
    }
   ],
   "source": [
    "# Register components\n",
    "llh.register_component(apt.ERBand, \"rn220_er\")\n",
    "llh.register_component(apt.AC, \"rn220_ac\", file_name=\"AC_Rn220.pkl\")\n",
    "\n",
    "# To see all the components\n",
    "llh.print_likelihood_summary(short=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3c54a04-f500-4df1-a41a-4cbfd327fe34",
   "metadata": {},
   "source": [
    "## Load parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c62de767-55cf-48ba-a588-805c685a770d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load parameters(and their priors) in simulation. Note: these parameters will be shared among components\n",
    "\n",
    "par_manager = apt.Parameter(get_file_path(\"er.json\"))\n",
    "\n",
    "par_manager.sample_init()\n",
    "\n",
    "parameters = par_manager.get_all_parameter()\n",
    "\n",
    "# Have to specify the normalization factor of each component\n",
    "parameters[\"rn220_ac_rate\"] = parameters[\"ac_rate\"]\n",
    "parameters[\"rn220_er_rate\"] = parameters[\"er_rate\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ec69b81f-6d71-429d-a05f-8f154a1691e0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'w': 0.014274272336585105,\n",
       " 'fano': 0.059,\n",
       " 'gas_gain': 31.3,\n",
       " 'drift_velocity': 0.0677,\n",
       " 's2_threshold': 500.0,\n",
       " 'field': 23.0,\n",
       " 'nex_ni_ratio': 0.10151185848273721,\n",
       " 'g1': 0.15351120578938598,\n",
       " 'g2': 16.44133119465396,\n",
       " 'p_dpe': 0.22178487445536998,\n",
       " 'elife_sigma': -0.008225975273542147,\n",
       " 's1_eff_3f_sigma': 0.26542639778743793,\n",
       " 's1_cut_acc_sigma': 0.13441954291258612,\n",
       " 's2_cut_acc_sigma': 0.9857380490961171,\n",
       " 'py0': 0.09164779636874743,\n",
       " 'py1': 43.42783454398944,\n",
       " 'py2': -0.2609731511199965,\n",
       " 'py3': 0.8929235275852672,\n",
       " 'py4': 0.669344602200388,\n",
       " 'rf0': 0.024649716781456686,\n",
       " 'rf1': 0.1436405417797817,\n",
       " 'ac_rate': 8.781138062834259,\n",
       " 'er_rate': 1988.3417164358577,\n",
       " 'rn220_ac_rate': 8.781138062834259,\n",
       " 'rn220_er_rate': 1988.3417164358577}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b37b66a7-d73e-4b38-9976-4b2637ccae1f",
   "metadata": {},
   "source": [
    "## Generate histogram under parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0c9b5a96-2fad-4940-9ab3-ec39073865e1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGiCAYAAAD5t/y6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfMklEQVR4nO3deVxU5f4H8M857IKAiGwCSmqiiUtqSpqlkmjdfnn1VdeltJs3rwaWUrmUmplG2a/UyrTN7Xcl05uWaamoKS6AgqKCigsgCAyoCAjKOuf3BzIxiorMeYZh+Lxfr3nVnHnm+zznuPD1WSVFURQQERERNSJyQzeAiIiI6EExgSEiIqJGhwkMERERNTpMYIiIiKjRYQJDREREjQ4TGCIiImp0mMAQERFRo8MEhoiIiBodJjBERETU6DCBISIiokbngROYqKgoPPfcc/Dy8oIkSfjll1/0PlcUBXPnzoWnpyfs7OwQFBSEc+fO6ZXJy8vD2LFj4ejoCGdnZ0yYMAFFRUV6ZU6cOIEnnngCtra28PHxwaJFi+5oy8aNG+Hv7w9bW1sEBATg999/f9DbISIiokbogROY4uJidOvWDcuWLav180WLFuGLL77AihUrEBsbC3t7ewQHB6OkpERXZuzYsUhKSkJkZCS2bt2KqKgoTJw4Ufd5YWEhhgwZgjZt2iA+Ph6ffvop5s2bh2+//VZX5tChQxg9ejQmTJiAY8eOYfjw4Rg+fDgSExMf9JaIiIiosVEMAEDZvHmz7r1Wq1U8PDyUTz/9VHctPz9fsbGxUX788UdFURTl1KlTCgDlyJEjujJ//PGHIkmSkpmZqSiKonz99ddKixYtlNLSUl2ZGTNmKB07dtS9f/HFF5Vnn31Wrz19+vRR/v3vfxtyS0RERNQIWKqZDKWmpkKj0SAoKEh3zcnJCX369EF0dDRGjRqF6OhoODs7o1evXroyQUFBkGUZsbGx+Pvf/47o6GgMGDAA1tbWujLBwcH45JNPcO3aNbRo0QLR0dEICwvTqz84OPiOIa2aSktLUVpaqnuv1WqRl5eHli1bQpIkFZ4AERGZI0VRcP36dXh5eUGWxU0fLSkpQVlZmSqxrK2tYWtrq0osU6RqAqPRaAAA7u7uetfd3d11n2k0Gri5uek3wtISLi4uemX8/PzuiFH9WYsWLaDRaO5ZT23Cw8PxwQcf1OPOiIiIgIyMDHh7ewuJXVJSgpZ2DriBSlXieXh4IDU11WyTGFUTGFM3a9YsvV6bgoIC+Pr6YszbH8D34c4N2LK7S4qNwm8/fIkvvlqGDh06CK1rz549+PSTjzFmZjjc2zykauzTsVHYvnoZli9ehIfbt1M1dk2Re6Pw8WdL8dHnX+Chduo8r/379mDZ559ixScfoEM7v/t/4QHtijqI8C++wfLPPsbD7dV97jVF7o3Cx4u/wgC0gBOshNUDAJkowVEUYqb/w2hj30xoXdVir+ZhdVo6vgmfgw4PtVE9/q6oGHy07Ht80CcAbR0dVI9/u0PZufgm8QLe790FbZ3shdQRnXUZ355KwSd/6492LZ1Ujx+Vkokv9ifg67B/4mFvD9Xj17QrPhEfr/sNnyz5Cu3aq/93ZVLiCbw/8x00b95c9djVysrKcAOVGIfWsDZwkXAZtFiryURZWRkTmLrw8Kj6DZqTkwNPT0/d9ZycHHTv3l1XJjc3V+97FRUVyMvL033fw8MDOTk5emWq39+vTPXntbGxsYGNjc0d130f7gz/R/vU5RaNLi8nCwDQvXt3dO/RQ2hdly5dAgD4+neBb8cuqsa+dus+enTrih7dAlSNXVNGZlU9Xbp2R5du3VWJmZ1V9Vx6dH0Ejwaon+heyq7qNXy0Wxf06Cru2VzKqno2rrBBK9z550BNRbf+Benv2Bz+juL+wq8pp6RqeLhHl054tEsn1eNfyqr6+6ZzS2d0FvDD/naaGyW36nNCZxcx9WmKq+oI8HRFgKer6vGzCosBAD06tEWPh9uqHr+mS5fzAABdu3VHgEp/9mtjjOkG1pBhLRk4TKWo0xZTpupAnp+fHzw8PLB7927dtcLCQsTGxiIwMBAAEBgYiPz8fMTHx+vK7NmzB1qtFn369NGViYqKQnl5ua5MZGQkOnbsiBYtWujK1Kynukx1PURERI2RhSSp8jJ3D5zAFBUVISEhAQkJCQCqJu4mJCQgPT0dkiRh6tSpWLBgAbZs2YKTJ09i3Lhx8PLywvDhwwEAnTp1wtChQ/Haa6/h8OHDOHjwIEJDQzFq1Ch4eXkBAMaMGQNra2tMmDABSUlJ+Omnn7B06VK94Z8333wT27dvx2effYYzZ85g3rx5iIuLQ2hoqOFPhYiIqIHIEmBh4Es2//zlwYeQ4uLiMHDgQN376qRi/PjxWL16NaZPn47i4mJMnDgR+fn56N+/P7Zv3643Brdu3TqEhoZi8ODBkGUZI0eOxBdffKH73MnJCTt37kRISAh69uwJV1dXzJ07V2+vmMcffxwRERGYPXs23n33XXTo0AG//PILunRRd+iDiIjImNToQbGA+WcwD5zAPPXUU1CUuw+uSZKE+fPnY/78+Xct4+LigoiIiHvW07VrV+zfv/+eZV544QW88MIL924wERERmZ0mtQqJiIjI1FUPAxkUQ52mmDQmMERERCaEQ0h1w9OoiYiIqNFhDwwREZEJ4RBS3TCBISIiMiEcQqobDiERERFRo8MEhoiIyIRIqPrhbMjrQfpfli9fjq5du8LR0RGOjo4IDAzEH3/8ofu8pKQEISEhaNmyJRwcHDBy5Mg7jvJpCExgiIiITIixjxLw9vbGxx9/jPj4eMTFxWHQoEF4/vnnkZSUBACYNm0afvvtN2zcuBH79u1DVlYWRowYIer264xzYIiIiMxUYWGh3vvaDjV+7rnn9N4vXLgQy5cvR0xMDLy9vfHDDz8gIiICgwYNAgCsWrUKnTp1QkxMDPr27Sv2Bu6BPTBEREQmxNBzkGquYvLx8YGTk5PuFR4efs+6KysrsX79ehQXFyMwMBDx8fEoLy9HUFCQroy/vz98fX0RHR0t8jHcF3tgiIiITEhVAmLoKqQqGRkZcHR01F2/vfel2smTJxEYGIiSkhI4ODhg8+bN6Ny5MxISEmBtbQ1nZ2e98u7u7tBoNAa10VBMYIiIiEyImvvAVE/MvZ+OHTsiISEBBQUF+O9//4vx48dj3759hjVCMCYwRERETZy1tTXat28PAOjZsyeOHDmCpUuX4h//+AfKysqQn5+v1wuTk5MDDw+PBmptFc6BISIiMiHGXoVUG61Wi9LSUvTs2RNWVlbYvXu37rPk5GSkp6cjMDDQ0Fs1CHtgiIiITIiswhDSg/ROzJo1C8OGDYOvry+uX7+OiIgI7N27Fzt27ICTkxMmTJiAsLAwuLi4wNHREVOmTEFgYGCDrkACmMAQERE1abm5uRg3bhyys7Ph5OSErl27YseOHXj66acBAIsXL4Ysyxg5ciRKS0sRHByMr7/+uoFbzQSGiIjIpBj7LKQffvjhnp/b2tpi2bJlWLZsmUFtUhsTGCIiIhPC06jrhpN4iYiIqNFhDwwREZEJYQ9M3TCBISIiMiHGngPTWHEIiYiIiBod9sAQERGZEAuoMISkqNIUk8YEhoiIyITIKgwhyQZ+vzFgAkNERGRCVJnEa/75C+fAEBERUePDHhgiIiITosoqJA4hERERkTFxCKluOIREREREjQ57YIiIiEwIh5DqhgkMERGRCZElyeBl0E1hGTWHkIiIiKjRYQ8MERGRCZEsJEiyYT0oUhPogWECQ0REZEJkCwmygQkMh5CIiIiITBB7YIiIiEyJhQxJNrB/QTL/0xyZwBAREZkQSZYgGbgTnQTzH0JiAkNERGRCZAsJsoEJjNwEEhjOgSEiIqJGhz0wREREJkSSDZ8DIymcA0NERERGxCGkuuEQEhERETU67IEhIiIyIZIFVyHVBRMYIiIiE1KVwBg4BwZalVpjujiERERERI0Oe2CIiIhMCCfx1g0TGCIiIhMiSSqcRq01/wSGQ0hERETU6LAHhoiIyITIFjJkAyfxyor5908wgSEiIjIhqiyjVsx/CIkJDBERkQlhAlM35t/HRERERGaHPTBEREQmhHNg6oYJDBERkSlRYQgJHEIiIiIiMj3sgSEiIjIhsiRBNnAjO1ky/x4YJjBEREQmRLKQDT/MUWv+Ayzmf4dERERkdtgDQ0REZEJUOcyxCZyFxASGiIjIhKiykV0TSGA4hERERESNDntgiIiITAgn8dYNExgiIiITIltAhTkwKjXGhDGBISIiMiGSLEEycB8YQ7/fGKjex1RZWYk5c+bAz88PdnZ2aNeuHT788EMoiqIroygK5s6dC09PT9jZ2SEoKAjnzp3Ti5OXl4exY8fC0dERzs7OmDBhAoqKivTKnDhxAk888QRsbW3h4+ODRYsWqX07REREZIJUT2A++eQTLF++HF999RVOnz6NTz75BIsWLcKXX36pK7No0SJ88cUXWLFiBWJjY2Fvb4/g4GCUlJToyowdOxZJSUmIjIzE1q1bERUVhYkTJ+o+LywsxJAhQ9CmTRvEx8fj008/xbx58/Dtt9+qfUtERERGI8uy7kDHer9kzoF5YIcOHcLzzz+PZ599FgDQtm1b/Pjjjzh8+DCAqt6XJUuWYPbs2Xj++ecBAGvXroW7uzt++eUXjBo1CqdPn8b27dtx5MgR9OrVCwDw5Zdf4plnnsH//u//wsvLC+vWrUNZWRlWrlwJa2trPPLII0hISMDnn3+ul+gQERE1Jqosozb0MMhGQPUU7fHHH8fu3btx9uxZAMDx48dx4MABDBs2DACQmpoKjUaDoKAg3XecnJzQp08fREdHAwCio6Ph7OysS14AICgoCLIsIzY2VldmwIABsLa21pUJDg5GcnIyrl27VmvbSktLUVhYqPciIiJqysLDw9G7d280b94cbm5uGD58OJKTk/XKPPXUU5AkSe81adKkBmpxFdV7YGbOnInCwkL4+/vDwsIClZWVWLhwIcaOHQsA0Gg0AAB3d3e977m7u+s+02g0cHNz02+opSVcXFz0yvj5+d0Ro/qzFi1a3NG28PBwfPDBByrcJRERkRiqLKN+gO/v27cPISEh6N27NyoqKvDuu+9iyJAhOHXqFOzt7XXlXnvtNcyfP1/3vlmzZga10VCqJzAbNmzAunXrEBERoRvWmTp1Kry8vDB+/Hi1q3sgs2bNQlhYmO59YWEhfHx8GrBFRERE+iRZhmTgHJbq798+0mBjYwMbGxu9a9u3b9d7v3r1ari5uSE+Ph4DBgzQXW/WrBk8PDwMapeaVB9CeueddzBz5kyMGjUKAQEBePnllzFt2jSEh4cDgO7mc3Jy9L6Xk5Oj+8zDwwO5ubl6n1dUVCAvL0+vTG0xatZxOxsbGzg6Ouq9iIiIzJWPjw+cnJx0r+qfxfdSUFAAAHBxcdG7vm7dOri6uqJLly6YNWsWbty4IaTNdaV6D8yNGzfumP1sYWEBrbZqVx0/Pz94eHhg9+7d6N69O4CqDDE2NhaTJ08GAAQGBiI/Px/x8fHo2bMnAGDPnj3QarXo06ePrsx7772H8vJyWFlZAQAiIyPRsWPHWoePiIiIGoPqlUSGxgCAjIwMvX+s3977cjutVoupU6eiX79+6NKli+76mDFj0KZNG3h5eeHEiROYMWMGkpOTsWnTJoPaaQjVE5jnnnsOCxcuhK+vLx555BEcO3YMn3/+OV599VUAgCRJmDp1KhYsWIAOHTrAz88Pc+bMgZeXF4YPHw4A6NSpE4YOHYrXXnsNK1asQHl5OUJDQzFq1Ch4eXkBqHqYH3zwASZMmIAZM2YgMTERS5cuxeLFi9W+JSIiIuNRYQ4Mbn3/QUcbQkJCkJiYiAMHDuhdr7m6NyAgAJ6enhg8eDAuXLiAdu3aGdbWelI9gfnyyy8xZ84cvP7668jNzYWXlxf+/e9/Y+7cuboy06dPR3FxMSZOnIj8/Hz0798f27dvh62tra7MunXrEBoaisGDB0OWZYwcORJffPGF7nMnJyfs3LkTISEh6NmzJ1xdXTF37lwuoSYiIqqH0NBQ3b5r3t7e9yxbPRpy/vx580lgmjdvjiVLlmDJkiV3LSNJEubPn683m/l2Li4uiIiIuGddXbt2xf79++vbVCIiIpMjySqsQnqAScCKomDKlCnYvHkz9u7de8cK39okJCQAADw9PevbRIPxLCQiIiITouYqpLoICQlBREQEfv31VzRv3ly3XYmTkxPs7Oxw4cIFRERE4JlnnkHLli1x4sQJTJs2DQMGDEDXrl0NaqchmMAQERGZkKp9YCwMjFFZ57LLly8HULVZXU2rVq3CK6+8Amtra+zatQtLlixBcXExfHx8MHLkSMyePdugNhqKCQwREVETVvOw5dr4+Phg3759RmpN3TGBISIiMiHG3om3sWICQ0REZEJk2fDTpJvCadTmf4dERERkdtgDQ0REZEI4hFQ3TGCIiIhMCBOYujH/OyQiIiKzwx4YIiIiEyJJKmxkJ5l//wQTGCIiIhPCIaS6Mf87JCIiIrPDHhgiIiITwh6YumECQ0REZEJkCxmygQmIod9vDJjAEBERmRBJllQ4jVpSqTWmy/xTNCIiIjI77IEhIiIyIZwDUzdMYIiIiEwIE5i6Mf87JCIiIrPDHhgiIiITwp1464YJDBERkQmRLCwgW1gYHMPcmX+KRkRERGaHPTBEREQmhJN464YJDBERkQlhAlM35n+HREREZHbYA0NERGRCJFmFVUgGfr8xYAJDRERkQjiEVDdMYIiIiEyIJEuGJzA8zJGIiIjI9LAHhoiIyIRwDkzdMIEhIiIyIZJsAUk2cCdeA7/fGJh/ikZERERmhz0wREREpkS2qHoZGsPMMYEhIiIyJbJc9TI0hpkz/zskIiIis8MeGCIiIhMiWVhAsjBwEq+B328MmMAQERGZEs6BqRMOIREREVGjwx4YIiIiUyLLKvTAmH//BBMYIiIiE8KdeOuGCQwREZEpkVSYAyNxDgwRERGRyWEPDBERkSnhKqQ6YQJDRERkQjgHpm7M/w6JiIjI7LAHhoiIyJRwCKlOmMAQERGZEu4DUyfmf4dERERkdtgDQ0REZEJ4mGPdMIEhIiIyJbJs+BAQh5CIiIiITA97YIiIiEwJVyHVCRMYIiIiEyLJFpAMTEAM/X5jwCEkIiIiUyLJf82Dqe9LqvuP9/DwcPTu3RvNmzeHm5sbhg8fjuTkZL0yJSUlCAkJQcuWLeHg4ICRI0ciJydH7Tt/IExgiIiImrB9+/YhJCQEMTExiIyMRHl5OYYMGYLi4mJdmWnTpuG3337Dxo0bsW/fPmRlZWHEiBEN2GoOIREREZkUYw8hbd++Xe/96tWr4ebmhvj4eAwYMAAFBQX44YcfEBERgUGDBgEAVq1ahU6dOiEmJgZ9+/Y1qK31xQSGiIjIlKi4E29hYaHeZRsbG9jY2NzzqwUFBQAAFxcXAEB8fDzKy8sRFBSkK+Pv7w9fX19ER0c3WALDISQiIiIz5ePjAycnJ90rPDz8nuW1Wi2mTp2Kfv36oUuXLgAAjUYDa2trODs765V1d3eHRqMR1fT7Yg8MERGRKVFxI7uMjAw4OjrqLt+v9yUkJASJiYk4cOCAYfUbARMYIiIiE6LmUQKOjo56Ccy9hIaGYuvWrYiKioK3t7fuuoeHB8rKypCfn6/XC5OTkwMPDw+D2mkIDiERERE1YYqiIDQ0FJs3b8aePXvg5+en93nPnj1hZWWF3bt3664lJycjPT0dgYGBxm6uDntgiIiITImRd+INCQlBREQEfv31VzRv3lw3r8XJyQl2dnZwcnLChAkTEBYWBhcXFzg6OmLKlCkIDAxssAm8ABMYIiIi02LkBGb58uUAgKeeekrv+qpVq/DKK68AABYvXgxZljFy5EiUlpYiODgYX3/9tWFtNBATGCIiIhMiyTIkAyfxPsj3FUW5bxlbW1ssW7YMy5YtM6RZquIcGCIiImp0hCQwmZmZeOmll9CyZUvY2dkhICAAcXFxus8VRcHcuXPh6ekJOzs7BAUF4dy5c3ox8vLyMHbsWDg6OsLZ2RkTJkxAUVGRXpkTJ07giSeegK2tLXx8fLBo0SIRt0NERGQ8ksVfw0j1fUk8zPGBXbt2Df369YOVlRX++OMPnDp1Cp999hlatGihK7No0SJ88cUXWLFiBWJjY2Fvb4/g4GCUlJToyowdOxZJSUmIjIzULeuaOHGi7vPCwkIMGTIEbdq0QXx8PD799FPMmzcP3377rdq3REREZDySVHUYo0EvqaHvQjjV58B88skn8PHxwapVq3TXai7JUhQFS5YswezZs/H8888DANauXQt3d3f88ssvGDVqFE6fPo3t27fjyJEj6NWrFwDgyy+/xDPPPIP//d//hZeXF9atW4eysjKsXLkS1tbWeOSRR5CQkIDPP/9cL9EhIiIi86N6D8yWLVvQq1cvvPDCC3Bzc0OPHj3w3Xff6T5PTU2FRqPRO1PByckJffr0QXR0NAAgOjoazs7OuuQFAIKCgiDLMmJjY3VlBgwYAGtra12Z4OBgJCcn49q1a7W2rbS0FIWFhXovIiIik2Jw78utl5lT/Q5TUlKwfPlydOjQATt27MDkyZPxxhtvYM2aNQCgW1/u7u6u972aZypoNBq4ubnpfW5paQkXFxe9MrXFqFnH7cLDw/XOhPDx8THwbomIiNSlSLIqL3On+h1qtVo8+uij+Oijj9CjRw9MnDgRr732GlasWKF2VQ9s1qxZKCgo0L0yMjIauklERERUD6onMJ6enujcubPetU6dOiE9PR0AdOcm5OTk6JWpeaaCh4cHcnNz9T6vqKhAXl6eXpnaYtSs43Y2Nja6cyEe5HwIIiIio+EQUp2ofof9+vVDcnKy3rWzZ8+iTZs2AKom9Hp4eOidqVBYWIjY2FjdmQqBgYHIz89HfHy8rsyePXug1WrRp08fXZmoqCiUl5frykRGRqJjx456K56IiIgaFUlS52XmVE9gpk2bhpiYGHz00Uc4f/48IiIi8O233yIkJAQAIEkSpk6digULFmDLli04efIkxo0bBy8vLwwfPhxAVY/N0KFD8dprr+Hw4cM4ePAgQkNDMWrUKHh5eQEAxowZA2tra0yYMAFJSUn46aefsHTpUoSFhal9S0RERKSy/Px8g76vegLTu3dvbN68GT/++CO6dOmCDz/8EEuWLMHYsWN1ZaZPn44pU6Zg4sSJ6N27N4qKirB9+3bY2trqyqxbtw7+/v4YPHgwnnnmGfTv319vjxcnJyfs3LkTqamp6NmzJ9566y3MnTuXS6iJiKhxk2V1Xibkk08+wU8//aR7/+KLL6Jly5Zo3bo1jh8/Xq+YQs5C+tvf/oa//e1vd/1ckiTMnz8f8+fPv2sZFxcXRERE3LOerl27Yv/+/fVuJxERkalRYxWRqa1CWrFiBdatWwegarpHZGQk/vjjD2zYsAHvvPMOdu7c+cAxeZgjERGRKVFjEq6JJTAajUa3dcnWrVvx4osvYsiQIWjbtq1ubuuDMq07JCIiIrPTokUL3dYl27dv121mqygKKisr6xWTPTBERESmxAx7YEaMGIExY8agQ4cOuHr1KoYNGwYAOHbsGNq3b1+vmExgiIiITIkZJjCLFy9G27ZtkZGRgUWLFsHBwQEAkJ2djddff71eMZnAEBERkVDR0dGYOnUqLC31044pU6bg0KFD9YppWikaERFRE6dIkgpnIZnWRnYDBw5EXl7eHdcLCgowcODAesVkDwwREZEpMcMhJEVRINWSVF29ehX29vb1iskEhoiIiIQYMWIEgKr931555RXY2NjoPqusrMSJEyfw+OOP1ys2ExgiIiJTosZZRiYyhOTk5ASgqgemefPmsLOz031mbW2Nvn374rXXXqtXbCYwREREpsSMhpBWrVoFAGjbti3efvvteg8X1YYJDBEREQn1/vvvqx7TNFI0IiIiAgAVViAZfpaS2nJycvDyyy/Dy8sLlpaWsLCw0HvVB3tgiIiITImkwmnSJpbAvPLKK0hPT8ecOXPg6elZ64qkB8UEhoiIyJSY0RyYagcOHMD+/fvRvXt31WKa1h0SERGR2fHx8YGiKKrGZAJDRERkSqp7YAx9mZAlS5Zg5syZSEtLUy0mh5CIiIhMiRkOIf3jH//AjRs30K5dOzRr1gxWVlZ6n9d2zMD9MIEhIiIioZYsWaJ6TCYwREREJqT6MEdDY5iS8ePHqx7TtPqYiIiImjoznAMDABcuXMDs2bMxevRo5ObmAgD++OMPJCUl1Sue6d0hERERmZV9+/YhICAAsbGx2LRpE4qKigAAx48fr/cuvUxgiIiITEn1YY6GvkzIzJkzsWDBAkRGRsLa2lp3fdCgQYiJialXTM6BISIiMiVmuArp5MmTiIiIuOO6m5sbrly5Uq+YpnWHREREZHacnZ2RnZ19x/Vjx46hdevW9YrJBIaIiMiEmONhjqNGjcKMGTOg0WggSRK0Wi0OHjyIt99+G+PGjatXTNO6QyIioqbODFchffTRR/D394ePjw+KiorQuXNnDBgwAI8//jhmz55dr5icA0NERGRCqvaBMWwSrqntA2NtbY3vvvsOc+bMQWJiIoqKitCjRw906NCh3jGZwBAREZFQBw4cQP/+/eHr6wtfX19VYjKBISIis3HubLKQuBfOnxMStzaKUvUyNIYpGTRoEFq3bo3Ro0fjpZdeQufOnQ2OyQSGiIgavbKKCsiyjDf+PaGhm2IwraJAa2AGYuj31ZaVlYX169fjxx9/xMcff4yuXbti7NixGD16NLy9vesVkwkMERE1etaWltBqtfjg9XFo6+WhevwzaekI/3696nGbCldXV4SGhiI0NBSpqamIiIjAmjVrMGvWLAwYMAB79ux54JhMYIiIyGwM7dcbj3Zqr3rcqPiTRktglFsvQ2OYKj8/P8ycORPdunXDnDlzsG/fvnrFMa11VkRERE2cVlHnZYoOHjyI119/HZ6enhgzZgy6dOmCbdu21SsWe2CIiIhIqFmzZmH9+vXIysrC008/jaVLl+L5559Hs2bN6h2TCQwREZEJURQFioGTcA39vtqioqLwzjvv4MUXX4Srq6sqMZnAEBERmRA1hoBMbQjp4MGDqsfkHBgiIiIS7v/+7//Qr18/eHl54eLFiwCAJUuW4Ndff61XPCYwREREJkYx8GVqli9fjrCwMDzzzDPIz89HZWUlgKpTqpcsWVKvmBxCIiIioc6kZwmvI1VzWXgdxtIQQ0hRUVH49NNPER8fj+zsbGzevBnDhw/Xff7KK69gzZo1et8JDg7G9u3b6xT/yy+/xHfffYfhw4fj448/1l3v1asX3n777Qdr7C1MYIiISIjyykrIsoR/hn/T0E1pVBpiEm9xcTG6deuGV199FSNGjKi1zNChQ7Fq1SrdexsbmzrHT01NRY8ePe64bmNjg+Li4gdqazUmMEREJISVhQW0WgXzJo6Bn5eb0LoOnjiDbzfVrTeA7jRs2DAMGzbsnmVsbGzg4VG/XY79/PyQkJCANm3a6F3fvn07OnXqVK+YTGCIiEioYY8/ih7+7YTXYy4JjPbWy9AYAFBYWKh33cbG5oF6Tmrau3cv3Nzc0KJFCwwaNAgLFixAy5Yt6/TdsLAwhISEoKSkBIqi4PDhw/jxxx8RHh6O77//vl7tYQJDRERkQtQ8jdrHx0fv+vvvv4958+Y9cLyhQ4dixIgR8PPzw4ULF/Duu+9i2LBhiI6OhoWFxX2//69//Qt2dnaYPXs2bty4gTFjxsDLywtLly7FqFGjHrg9ABMYIiIis5WRkQFHR0fd+/r2vtRMMgICAtC1a1e0a9cOe/fuxeDBg+sUY+zYsRg7dixu3LiBoqIiuLkZNqzIBIaIiMzGmdR0IXHPXbwkJG5t1FyF5OjoqJfAqOWhhx6Cq6srzp8/X+cEplqzZs0MOkKgGhMYIiJq9MorKiDLMsbP/rShm2KwxnCUwKVLl3D16lV4enoKredemMAQEVGjZ2VpCa1Wi4ARk+DQykv1+AVZaTj120rV45qKoqIinD9/Xvc+NTUVCQkJcHFxgYuLCz744AOMHDkSHh4euHDhAqZPn4727dsjODi4wdrMBIaIiMyGV7fH4dLWX/W4uWeOGS2BUXMVUl3FxcVh4MCBuvdhYWEAgPHjx2P58uU4ceIE1qxZg/z8fHh5eWHIkCH48MMP6z2nRg1MYIiImrjzV/KFxM3Ivy4k7r0UZqWJiau5KCRubRSosArpAcs/9dRT9xx22rFjh2ENEoAJDBFRE1Wh1UKWJLyxeW9DN8Vg5RUVgCQj+pu5Dd0UuovQ0FDMnz8fLi4uqsRjAkNE1ERZyjK0ioLWg8bDxrl+O6zey/WMJFw+slX1uLWxsrQEFC16j3odzd28VY+fn5mCoz/Xb8O1B6VVFGgN7IIx9PtquXTpEry9q349IiIiMH36dLi4uCAgIAC///77HfvUPAgmMERETZxzh8dg79VBSOzLR7bidJr4JcipWbkAAN8eT6DVQ/Xbmv5esk7FGy2BUeNEadNIXwB/f3+0bNkS/fr1Q0lJCTIyMuDr64u0tDSUl5cbFJsJDBFRE3fzipi9U0ryMiHLEl6Zt0RIfHPVEKdRi5Kfn4+jR49i//792LRpE5555hm4u7ujtLQUO3bswIgRI+Du7l6v2ExgiIiaqAqtFpBkpPz8idB6psyYjdbebe5f0ADH4mLw05ofcC0zRUj8fEGTg81deXk5HnvsMTz22GNYsGAB4uPjkZ2djaCgIKxcuRJvvfUWfHx8kJyc/MCxmcAQEZmwlIIiYbFzbpQAiha9Xnwdzd1aqx5fk5yA05Eb8cTAp9G5a3fV49dUUV6ODdIq7PniPaH1GIUKZyGZyhiSs7Mzunfvjn79+qGsrAw3b95Ev379YGlpiZ9++gmtW7fGkSNH6hWbCQwRkQkq11ZCloBZ0ceF1+XToz9cBcwbAYDTkRuFxL2dpZUVFEUrrLcn9cJZfLPEOLv8aqFAa2AGYuj31ZKZmYno6GgcOnQIFRUV6NmzJ3r37o2ysjIcPXoU3t7e6N+/f71iM4EhIjJBVrIFtAoQNms2vH3FDL/EH47BulU/CIndUET19sRFHzRaAmNOXF1d8dxzz+G5557DihUrEBUVhdOnT2PcuHF4++238fLLL+Oxxx7Dvn37Hjg2ExgiIhP25OAh6NKtu7D461b9gPzMVCGxr+dmColr7hQVhpBMZBX1HZycnPDiiy9iwoQJ2LNnD5o1a1av5AVgAkNE1GSVl5VDkmT8+ZXYeSMp588KjQ8AmZeMt1OuaOa0CqmmEydOoHXrqrlWbdq0gZWVFTw8PPCPf/yjXvGYwBARNVFW1lXzRvqODoWjgEm8mUlxOL1nM2aGvqZ6bGp8am5al5iYaHA8JjBERA8gpcA45/tkFhUbpR4AaPtof7i16ywkdtKunzH57ffg5SN2GfXxuBj89/9WCuvtSUs5f/9CKjHnISQ1MYEhIqqD8vJyyBIw42BCQzel0ek38Gn4B3QTWkdFeTk2rVttFr095rQKSSQmMEREdWBlZQWtArwydSY8BG/KBgBJR2PxW8Rq4fUAQN4lMZN4C404idfSygparRbdRk6CvYDhsMLMVCRuWal6XKo/JjBERA/gsSeD0OGRrkap67eI1bhw7sF3KK2rtJQLkCQZO5fOElYHAKSeF3cP1bIyqibxenXvh5Zt/VWPrzlzzGgJDIeQ6kZ4AvPxxx9j1qxZePPNN7FkyRIAQElJCd566y2sX78epaWlCA4Oxtdff613HkJ6ejomT56MP//8Ew4ODhg/fjzCw8NhaflXk/fu3YuwsDAkJSXBx8cHs2fPxiuvvCL6loiIhKsoL4MsywibLH5I5LFRIXAScILzpaQjSP7zF8x589+qx76bQlFLwrPThMStjTmdRi2S0ATmyJEj+Oabb9C1q/6/VqZNm4Zt27Zh48aNcHJyQmhoKEaMGIGDBw8CACorK/Hss8/Cw8MDhw4dQnZ2NsaNGwcrKyt89NFHAIDU1FQ8++yzmDRpEtatW4fdu3fjX//6Fzw9PREcHCzytojIhJ05L+YHWFqGcfc0sbSyhlarxbx3p6Otr8/9v1APhw4fwbcr16Lto2JOcAaA07s3wXPgONg4ewiJX+36xRO4emwnDq6YK7QeY6jUVr0MjWHuhCUwRUVFGDt2LL777jssWLBAd72goAA//PADIiIiMGjQIADAqlWr0KlTJ8TExKBv377YuXMnTp06hV27dsHd3R3du3fHhx9+iBkzZmDevHmwtrbGihUr4Ofnh88++wwA0KlTJxw4cACLFy9mAkPUBFVotZBlGeOmmcFZODUMDRqEHt0ChMX/duVa5F0ScwBiQe4lAICNswfsXMUkYdVK8zWAosXo0Blwa61+b1Jm6gX899slqsel+hOWwISEhODZZ59FUFCQXgITHx+P8vJyBAUF6a75+/vD19cX0dHR6Nu3L6KjoxEQEKA3pBQcHIzJkycjKSkJPXr0QHR0tF6M6jJTp069a5tKS0tRWlqqe19YWKjCnRKRKbCUZWi1WvQe9TqaCxgO0Zw5hlM7NyL9gvhN2QBAY4SN2crLqzay2/XFu+IqkWSkbV4kLv5tHn1iINp1Vn+OUlJctNESGA4h1Y2QBGb9+vU4evRorSdMajQaWFtbw9nZWe+6u7s7NBqNrkzN5KX68+rP7lWmsLAQN2/ehJ2d3R11h4eH44MPPqj3fRGR6fPtIWY4RFtejjPyz/j47ddVj91QrG4dgDjuzRlCVlYlHY3Fth/XYNwbM+Du7at6/JpOHT2MbevXCK3DWLSKgkomMPelegKTkZGBN998E5GRkbC1tVU7vEFmzZqFsLAw3fvCwkK9nQGJiO5GvrVM1xg/jAHj/kDuPUDcyqptP65BrwGD0d4IK7fMJYGhulE9gYmPj0dubi4effRR3bXKykpERUXhq6++wo4dO1BWVob8/Hy9XpicnBx4eFRN8vLw8MDhw4f14ubk5Og+q/5v9bWaZRwdHWvtfQEAGxsb2NjYGHyPRGS6rmWKmc9x/dZ8DmP9MAb4A7mpqjoLydAeGJUaY8JUT2AGDx6MkydP6l375z//CX9/f8yYMQM+Pj6wsrLC7t27MXLkSABAcnIy0tPTERgYCAAIDAzEwoULkZubCzc3NwBAZGQkHB0d0blzZ12Z33//Xa+eyMhIXQwialoqtFpIkow9X5jXJN4zZ88Ji52WniEsdk0ZKeLuoVrOpXQAwCVBW/5npolJjGvDVUh1o3oC07x5c3Tp0kXvmr29PVq2bKm7PmHCBISFhcHFxQWOjo6YMmUKAgMD0bdvXwDAkCFD0LlzZ7z88stYtGgRNBoNZs+ejZCQEF0PyqRJk/DVV19h+vTpePXVV7Fnzx5s2LAB27ZtU/uWiKgRsJRlKIoWQydMg4un+pN4U0/GI/rXCKP8MAaqfmDKsoxXJk0RXldGipiJyZlpFyDJMj6dHiIk/u0kWcbSWaFGqYsaXoPsxLt48WLIsoyRI0fqbWRXzcLCAlu3bsXkyZMRGBgIe3t7jB8/HvPnz9eV8fPzw7Zt2zBt2jQsXboU3t7e+P7777mEmqiJ69T3SXh37HL/gg+osqIcsb+tN9oP42oDx70BZw/1EzIAuHjiCBJ2/oxF74i9p9aDxovfByYjCZePbBVW180r6ciO+lH1uLXhKqS6MUoCs3fvXr33tra2WLZsGZYtW3bX77Rp0+aOIaLbPfXUUzh27JgaTSQiAdKKbxitrqybN4XGt7CsmsT7/pz34NdG/FlIh2Ji8O33K9G+9wB4dnhEWD1Ht2/EO+/Nha+v+vd0JDYGa1d+B5sWHrBzFTvxuTS/aoWqc4fHYO/VQfX4hWknjJbAVKqwCsnQ7zcGPAuJiFSnhQIZwLyk0w3dFNUNHfI0enTvbpS6vv1+Ja5kiJt7ka+pmpg8KGgIArp1Vz1+eXkZ1q6SkfLzJ6rHvpubV9KFxC25eklIXKo/JjBEpDoZErQA3n/5ObR1b2mUOqNPXcC3v+9HzsULQuLnZRv3B1hZeTkkWcbmRdONWq+arKysAUXcCdE15Z6Ox/l9W4yaLImiheGriJrAHF4mMEQkTnCvR9Cjvfg9UwCgrKIC0vaDiFjwllHqE83aygqKVoshr06Fi6eY/arSTsYjZkuEkNg1iToh+nbn9/4Ch46DYNGsheqxK65fRvH5KNXj1qZSq6DSwAzG0O83BkxgiMgsWFtaQtFq8caM2fAWMJ/j6OEYrF/zg+px78e/71PwfljcHBhjJDDGoK0oByQJRcl7GropBlNUmMSrcA4MEVHjMmDw03ika3chsdev+QFnko1zFlLqxaqzkHIvitnXBADysqv2gTl3NllI/PT0qnsozBRzQnhNN/NyAUVB++deQzNXL9XjF2WnIWU7NxY0JUxgiIjqoLy8HLIs45//mmi0OiVZxvqFbwutQ5ZlvPHvCeIqkGQcXDFXXPzbuHUJhJNvR9XjXj2XYLQEplKpehkaw9wxgSEiqgOrW2chifoX/u3yzp/Apf2b4fP0q7B1EbOHSmFaInJit6BVv5dg7eimevzirNPIP/4H7Pz6wcLOSfX4NZXnZ6I067jQOoyF+8DUDRMYIjIrKefEDPFcujUc0szVEw4e4veBuXElCwBg6+KBZm5iJkKX5FXtnWLt5A4bF/U3yyvNzwYg4WbqQdVjEzGBISKzUFZRAVmWMT3kNXGVSDJOrJp//3Iq1nfup48E1yEh8/fPhFbhMWAcrJ3dhdZRnHEKV4+Zx1EyXIVUN0xgiEiYMxkao9WVdbUAWq0WASMmwaGV+kM8l88ex/k/fzb6EJLfMxNg19JTSB35KSeRdfBXuAaOgZWTgO33M08h/+R2WDu7w9ZVzFLwamX5OQCAIk2akPjFOWI2yKsNh5DqhgkMUROTj3LhdRSiHLIk4dX/XS28rtt5dXscLoL2HDn/589GH0Jy7dwXzX3Un5RaLevgr7ASNIRUdi0TkCSkb/lU9di1MnYPGTUoJjBEJuL0WXHLZQHg3IU0SAD24KrQenQU4LmJb6GVl9h/eVc7dyIO+zf9B4VZaULiX9ekN8gPyOKci8Ji37h8CZAkZG9fLKwOALBt0xeybXOhdVQUaFCmOQnXvmNg6aj+cFVZ3iXkxW1UPW5tuAqpbpjAEDWw6uW5r7w+1Sj1vf/PEfDzaCW0jkOJZ/Htb3/C1csHHm3bC62rWk56GiDJiP5G7JLdvqND4Sh4W3wAyEyKQ9LuzTj1fwuE19V/7BQ4uat/T+mJcTi582eUXIxRPfbd2LftCVu3dqrHvZGZZLQEhkNIdcMEhqiBVS/Pff2d99BawA6y1RKOxGDj2pUY+lg39Hi4rbB6gKoJtd9t24dV86YKrac2YbPE7MQbfzgG61b9AEe31nDx9lM9/u0KczMBRYuHnv0XbAXNgSlIOYnMA7/goV5PwL1dZyF1nNz5Xzj1eAGWDmKT5tLcsyhK3iW0DjItTGCITET/QU+jU0A3oXVsXLtSaPxq1dv6N8QQ0pODh6CLiJOVy8oRIa3CzqWzVI99Ly0f6QtHgXNgMg/8Iix2NTvvbrB2FZ/0FSXvQlmemEM3y65lColbG61WgdbAVUSGfr8xYAJDRMIEPD4Qvv5djFbf/k3/ERbbytoKiqKF//MTjbYKKW3fJtzQiJsDU3I1W1jsmsoLssTXUaipms+zU+x8HmPQqjAHpgnkL0xgiEic7DSxE5NrupyVYZR6mrl6oblnW+H1FOVUTRpOWvuh8LquZqQIiXst6yIgSbga9bWQ+LVpHjAcFg6uqsetKMhG0Snj7DPDOTB1wwSGiFRXVlEBSZYbZA7MhXNiDiZMS7kASDKO/jBPSPy7cXh4IORmLYTELruahtJLCdj2+Uwh8as9/D8TYe8qZh5PtavnTyI9ahNsWwfA2kX9OVAluWeNlsBQ3TCBISLVVc+B6TZyEuyNsGIHAHJPxyMlagvCJgvciRfAiyHT4eYtfl7PqSMx2LM5AkVn/xRel6hei7LL53Hj/F64B4g5YPF26VGbhNdhDJWKgkoDe1AM/X5jwASGyESkCuo5qJZ56yyfM+ni5yOkai4DABxaecHRq63w+gCgODcTWq0Ww14LQ0tP9TdlSzkZj0Ob16HHE4PwUKcA1ePXZs+mdbBu3QOyjZg9VCqu56DicjIsHVxh6aR+D0ll0RXVYzYFnMRbN0xgiBpY9T4w773xb+F1ybKEf4Z/I7weAIAk4+AKsXuy1KZz3yfh01HMxOFDm9cJiVubyvJyABLKMo8JrknCtejvBddBpD4mMER1dObcBSFxM7NzoNVqhQ9NJB87gsgNa4WerVOt+oydJ156A84CNkirzaVTx3Dsj/VGqcsYLKysACjwefqfsGmh/jlFAHD9YiJyYn+DY7cRsLAXMYR0DsXnxA+B1VRRIGZlVWWB8c71qoQKO/Gq0hLTxgSG6D7Ky6p6SMa/MUNoPcYYmojcsBZ2LT1h7y72LJ+bt5bnOru3Rkufh4TWVS0/p2qfjpyLYhLNq9lV+4tkppwTEv92uZeqVlW16PgYHFo/LKyenNjfYOsVAOuWbYXEN1YCU1lZAXPpTeIqpLphAkN0H1bWVTvlfvDOG2jrq/7cikNHjuKbteJ7DirLywFJNsrW9AAgSTJ++0xs0ndHnbKM/8wPExr/q3enCItfm5u54k5BLr1W1atQUSio16LYeHNgLCwsASiw8OoNydpB9fjakgJoNUdVj0v1xwSGqI6GDhqARwPEbLf+zdr1wv9lfzVHAyhavDptFjy9fYXWdTL+MLZErMKM9+bCu434k5sBIC4mBmtWfids2XF5Xjpupseh/d9eEz4EBwB55xOQeWgrzm0IF1yThLyD3wqt4Xp2mtD4AFB8pSoJk518INsLOLagMMtoCUxDrEKKiorCp59+ivj4eGRnZ2Pz5s0YPny47nNFUfD+++/ju+++Q35+Pvr164fly5ejQ4cOBrXTEExgiBpYeXm5Uf9l3/fJwXi4i9gjCwBgS8QqDHx6CLoK2Na/NuXlZcAqSfiyY9dH+sLRCEuCASDz4BZh81OAv+ao2LYfDItmLqrHL7+airKso0hY+YHqse9GKcmHVkTc0gIBUWun1SqoNPIqpOLiYnTr1g2vvvoqRowYccfnixYtwhdffIE1a9bAz88Pc+bMQXBwME6dOgVbW1uD2lpfTGCIGpiVlRUUrRaBY0KFnAhcLfP0MZzcvgEXL4ifw5F9Sdywx91YWVkDioLWg8bDxln9Sa/XM5Jw+chWFAvc2r+m6nlEIuenAFVzVKzdOsLSSczxCGWZccKGdWqqvJ4N5WoyKlL3CK3HXA0bNgzDhg2r9TNFUbBkyRLMnj0bzz//PABg7dq1cHd3xy+//IJRo0YZs6k6TGDILJ1XcU+VjHTj/MDy6ynuRGAAqKyoQNLO/2LhW5OF1XG788li97ap6dLFql8nmxYesHNVf4isJC8TkGScXDNf9dj3Imp+CvDXHJXKolwx8W/kAQAkm+aQbJ2F1FFNLitCJRTIHj2EzYFRck+oHrc2lSr0wFR/v7CwUO+6jY0NbGxsHihWamoqNBoNgoKCdNecnJzQp08fREdHM4EhUkP1niqid2MVIe+SmPNoqhXmZkGrFdG5XjtZlhH67wlGqw8AIMlI+fkToVU49XgBlg4C5ljcpkRzBsXn/hQ+PwWQUHziv0LjG7NXRHb0htRMwJDb9WxUNsIExsdHf2uG999/H/PmzXugWBpN1WRvd3d3vevu7u66zxoCExgyK1ZWVSuGXn/nPbT2VWfyaMKRGGxcu1KVWLUpLy+HJMnYvniWsDpqeuH1d9DKS+xW+GePx2HXxrXo8cJkOLQSf3IzAOSePY7kXf+Fg/8gWNqpP4m39Fo6bqYdgWXzVrASNNxSU0XRZQAKmnUaClnA/BQAqLiaipK0aGFJWWnuWRQl74JDl+dgad9S9fh6dV25gJsX9gutw1gqtVAhgan6b0ZGBhwdHXXXH7T3xZQxgSGz1H/Q0+gUoN5E1Y1rVwrdyE5RtHhzxmyhK3aOxsbgxzU/oFu/gfAzwlb4uzauhUMrLzh5+QmvCwAKNOmAJKHojMB/7Rv5ZGUAsHbvBCtn9ZfvA8BNAEiLhp13N1i7ivl1KkreBVvPR2DVQuzKNwC4eWE/lBIxk22V0sL7FzJBjo6OeglMfXh4VM0py8nJgafnXyvwcnJy0L17d4NiG4IJDNF9VA9Lid7IbkDQ03ika3ehdfy45gdkpZ4XWgcAaNJTAEnG/q/nCK/rdqLODqo+N8g1cAysnMTsjFvTzcxTyD+5HZXXc4TVob01R8UcKLc2sqtMj2rophhMzSEkNfj5+cHDwwO7d+/WJSyFhYWIjY3F5MnGm1N3OyYwRPdRPSzl8/SrsHVR/wdXYVoicmK3qB73dtXLtb+e/Ybwuqq1CBgCSwcxwx+3K8m5gOspccLPDnLw6wVbt3ZC6wBu/UA+uQPX4yOE12UOpFsb2XES718xHkRRURHOn//rHzepqalISEiAi4sLfH19MXXqVCxYsAAdOnTQLaP28vLS2yvG2JjAENWRrYsHmrmJWN1SNQku5exZ1WPXlJOVBUWrRf+xU4Qu1waA9MQ4nIz8GddO7hRaT23aDnsVti7qbzRXkHoS2Ye2oDQvQ/XYtakougJAgXPvUbBycBNSx82cMyg6ZfxfI5HMYRJvQ4iLi8PAgQN178PCqna0Hj9+PFavXo3p06ejuLgYEydORH5+Pvr374/t27c32B4wABMYovsqv7UF/7mfPhJWhyzLeCfEOCunnIxwPlFBTiagKGgzVEwyUWudqSehid4Cl0590dxb/bODtJUVyI7+DdnbF6se+17sfR6FTStxv17mlsCYg4bYyO6pp56Cco/deyVJwvz58zF/vnG3EbgXJjBkllJV3AcmNzsLULT4+6S34Npa/dU75xPi8OfP/8GEMLFb/CccPoRtP63Dts9nCqvjdqJ6rWpTkiduvxQAkC0sAUWBbbuBkOychdYFAJX5GSi7FCe8HjI9lYoKQ0g8zJFIXWfOit0F9nxKCmRZxntv/Fv12AH9BqKNv5jVO3/+/B/0fTIIHQVv8b91/f/h6VenooWHmFUt1VKOH0bc7/9FcoS4XquGIts5Qxa0tX9Nys184XUQNWZMYMgoKiqqVvK8Msk45/08Oe4NOLmr80P60qmjOLptPbIFrd65klk1p+LiBbFzYKq39+/42JNo/fAjQusCgLhtG+Ad9E/YthC/YgcACi8mIvfwb7iRI2bn5BuXLwGQcCNxs5D4RNVMbRWSqWICQ0ZhaVm1kmfIq1Ph4iluE7W0k/GI2RKBdr0GwLO9OtvyV5aXQ5I34Lu5U1WJVxtZlrEgzDjLEXPTxexnU9M1zSUAQIuHH4N9a+OcVqutrEDukW04s26h0HqMsSkbYNyN2coLsoTErdqMjx4UE5i6YQJDRuXf9yl4C/7Xf8yWCFzNUO+HdFFeDhStFk+89AacBazeuXTqGI79sR4tewTDSuCS4xuaC7h+IQ4bPnpbWB23u3nZeIc6lhVeARQtWvYZAytH9Vft3Mw+g4LE7UbblA2A8ARG0VY0yOZ8RGpgAkNmpbKiaq+TXz9Vf9M5Z0Grd/Ky0gBJxtVjO1SPXRu7to/Dws5JaB1l1y6hLCcJ5zeGC62nNvZtHhW2T0tB4nYhce+lLP+SsNgVN/IAReDeKcW5UK4a70BPc1GhVWBhYA9KBXtgiNSVe1HsLrAFlzVQtFq89MYMuKm0oifpSDR2/jcCv30mdifegBGThJ4blHPmKFL2/YqbaYeE1XE7n6f/CRsjzYG5fjERObG/GaUuY1Aqq3pHLu/5QnhdwvZOAVDJBOaBcQipbpjA0B00aeonGbnpaZBkGesXGmf4oueAwWjfuatq8XZs/A+GvRaGlp7qr95JORmPQ5vXwaGVFxy92qoev1rR5SwACvyemQC7lmL3ZslPOYmsg7+iRcfH4NBa/T1Z7sacEhjp1rJtS+++kGzVPxoBALSFGlTmnhQSm+qvIfaBaYyYwJBOeVkZJFnG6nnThNUxdMI0uAhIAqqlnoxH9K9itl7v3PdJ+HTsonrcyvJyHJJ+RPQ3c1WPXRu7lp6wdxd3aCQA3LxatSfLzVzjzYEpvVa1o3HZNTFDLuWFuQCAikKNkPi3qyi+CgCwaNEGsr2YnXgrACYw1GgxgSEdK2trKFotXnz9HbRSeUO1s8eOIHLjWnTq+yS8BSQBNUX/GoFLF9Tbbyb31vLjnItiVu/kX8kBFHFnLVUruHAcufF/4NT/LRBWhx5JxrkNRp4DI0nQRC4RWQHyY1cLjE90ayM7Azei40Z21CR1f2IQ/Dqpv2Fb5Ma1wpKAapcvVQ1VfTYjRNW4kizjP/PDVI15O9G71pbkaQBFixa9RwtZpVPTzewzuH5qR4McJWDnPwQWduqv5iq/lobStBjYPfw0LOxbqB7/jvquXkRpeqzwesj0cA5M3TCBIaOovHUScsSCt4xSX6e/T4S9qzoTYq+cP4G0PzehzRPPwdZJ/YmO+elnkXMyWuhZSzVZObrBuoXYnXirh1tc/PsIOZfobjTRW2Dj5g9LZzH3V5oWA2v3jrB0EnsYpq4+JjBEd8UEhozCwsoKilYrfKXN5bPHcf7Pn2Hv6oXmnm1ViXk9Jx2QZFzcL3aCqP/zE9FMpaSrNleSjyH94Bbk7l4qrA4iMhx7YOqGCQwZleiVNtc1VclG3HfzVI8tagJy9cTjZiomXbW5cSULUBQ4PfoiLJuLHUIqzUlG0ZlIoXUQmSsmMHXDBIbukJmi/oGL2RcvAJJstJU2Ho89AytHdbZ7v5GdgqtJB4VNQK6sKEf0lvU4+sM81WPXxrK5G6ycxQ6BVFzPFRr/nnUXiam78maekLhEVD9MYEinehn1svfeEFaH/6C/o1kLcSf5XklLRvrRA9Ac/l312KImIBdcrlqFZOfXT+gOuWXXMlCmScTVfV8Jq6MhaSsqAEgoOrq+oZtCZJBKRYtKrdbgGOaOCQzpVC+jDhwTCieVz/y5lBiHpF2bcWaPcU7ydXrkaViqtFKkJDcFN9KPCp6ALOFm6kGB8f/i3HsUrBwEr0LKOYOiUztxI1fMydC1KSu8DECBVZt+kG0dVY9fWZCFCs1x1eMS3Y4b2dUNExi6g1/PJ+DeTp2TnGtKjPwZ88Mmwc9b3PDFvsPx+P6nX1GQpP78iy4j/q3ayqaarpw7jgt/boJr3zGwdHRXPX61kuzTyD+5HfY+j8KmlfpnOtWkaCtQJO1CcoRxVlbVZOniBwsHMc+RCQyR6WACQ3fIu5SiesyCnEwAwLAn++HRLv6qx6/p+/W/wK33MFg7qDMHpjgnBddOHYJnwONo0VZM2y/8uQlWju6wdhG3vLmiMAeA2MMBdXXdyAMULWTPnpBsxGyDfzttUQ6UK6eNUheRSJVaBTIn8d4XExjSKS8rgyTJ2L54VkM3pd7KyssBSUbukT9Uj12YnaZ6TKB6mbaE7J2LhcTXY6TDAXXV2TSHZOtsnLpKr8P8/8qmpqBCC0gGn0atUmNMGBOYRiI5WfyJrlnZ2VAULdoOU3/31ILUk8g+tEXVmLWxtrICFC3sOzwFi2bOqsQsvXoRpZnHEfvt+6rEuxvRO+RW747r3GuU8GXUJVlJKEreg8q0vULrITJH7IGpGyYwJq6iohyyLOO1Ca8arU5bF080U/mwv5K8qsP9Tl9IVTXu7c6mXQQkCcXn9qoe23eImLOKCtMSoYnZAvs24uemXD+1A3Y+PWDjKrYeAChK3g3b9oMg2zkLrwsAKq5loOzSEaPURUQNjwmMibO0tIJWq8XAcW/A2UPs9u8XTxzB0R0/48y6hULiy7KMcWHG2QfGO+ifsG2hTrJReDERuYd/E3ZWUUme2FOUq1Vv729MVq06wNJR3O7Ct2MCQ+aAPTB1wwSmkWjfewA8OzwivJ6j2zdidOgMuLVWN1k6cywOOzaswUtvzICbyidd13T66GH8sX4NbFt4wK6VOvXcvHoJkGScXS9wRY0kcYt/IgLABKaumMDQHR59YiDade6qetwdG9bAvbUvvB/qoHrsatmpKYAk4/zGcNVjewwYB2tn9ZfnFmecwtVj2+AaOAZWTuoPUVW7mXkK+Se3C4tPRGRMTGAaiSsZ6i9tvl2+pmoI41LKedVjZ6WnQpJlfDYjRPXYtek2chLs3dTZb+ZycgLO7fkZju17oZlHe1Vi3u7qsW1w8OsFW7d2QuJXYwJDZPq4kV3dMIExcRUV5ZBkGZsXTTdKfZIsY+msUGHxnToHwcLeRVj8sitpKEo9rOqhkcW5VXvYlFzJUCXe7cryc4TEvZvy/EzhdTTkWUhEjV2lVjF4GTWHkKjBWVpaQdFq8dzEMLT08hFa14UTcdi/aR0mhM2Cp8rzVBIOH8LWn/6DglO7VI1bK0nGwRUqTxaWJKRv+VTdmEamVFYAkoQrfxpvHxgiIlGYwDQSjzw+EL4CTkK+3f5N6+Dp7Ys27R5WNW72pXRAUTBy8jtoJTARO3s8Dnv+uxaO3UbAwl6dQyPLLp9D8bk/hZ0hVH1uUGmemB6eahVFVwBFgU3bJyAJXtpcmX8J5dnHhNZBZK4URYFiYA+KorAHhkyEJk39eSm3y01PgyzLWBA2WVgdrbx84OUnZh4JAFzOqkoCLB1cYemozmZ8lcVXAEDYGUJV5wZFInu7EXbiBWDp2h4WzcVNFgaAMoAJDFE9abWKwXNYOAemHsLDw7Fp0yacOXMGdnZ2ePzxx/HJJ5+gY8eOujIlJSV46623sH79epSWliI4OBhff/013N3/WuGRnp6OyZMn488//4SDgwPGjx+P8PBwWFr+1eS9e/ciLCwMSUlJ8PHxwezZs/HKK6+ofUsNqnoOzOp504xWp4hl1Elxsdi9OQIr5ryhatzaScg7+K3qUUWdIVR1bpCiaq9Rbap7koiIzIHqCcy+ffsQEhKC3r17o6KiAu+++y6GDBmCU6dOwd7eHgAwbdo0bNu2DRs3boSTkxNCQ0MxYsQIHDx4EABQWVmJZ599Fh4eHjh06BCys7Mxbtw4WFlZ4aOPqvbiSE1NxbPPPotJkyZh3bp12L17N/71r3/B09MTwcHBD9RmzcUU2DZrpu6DUEl+rgaKVosXXhc79AJUDb/s2rhW2DLqXT//R8j8mppOxh/Gr+tWoVW/l2Ct0rb8RZeSUJC4U/gZQrZeAbBu2VZoHUxgiEyfoigGDwFxCKketm/XX6a5evVquLm5IT4+HgMGDEBBQQF++OEHREREYNCgQQCAVatWoVOnToiJiUHfvn2xc+dOnDp1Crt27YK7uzu6d++ODz/8EDNmzMC8efNgbW2NFStWwM/PD5999hkAoFOnTjhw4AAWL1581wSmtLQUpaWluveFhYVVbVw4Q+3HoDrRQy/AX8MvIpZR52ZW9V70fTIIHbt0Uz1+Tb+uWwVrJ3fYqHSyc1lhLqAosG03UMjckcr8DJRdilM9LhE1TopWhTkwHEIyXEFBAQDAxaVq6Wx8fDzKy8sRFBSkK+Pv7w9fX19ER0ejb9++iI6ORkBAgN6QUnBwMCZPnoykpCT06NED0dHRejGqy0ydOvWubQkPD8cHH3xwx/WZY/4Gf1/jbXf+IKJOnMGq7Qfw9WxjDL2IX0Z98cJZYbEB4FJaCiBJyPz9M9Vjy3bOkAUM8Sg381WPSUSNF+fA1I3QBEar1WLq1Kno168funSpWkGj0WhgbW0NZ2dnvbLu7u7QaDS6MjWTl+rPqz+7V5nCwkLcvHkTdnZ2d7Rn1qxZCAsL070vLCyEj48Pgnp2Qf+u6q66UdPK36MQOuM9tPZW94DF2x2Li8GGNSvRa9gLaN6ylaqxs86fwbnDe4VOEK5J9SGkkztxI3GzKvGIiMhwQhOYkJAQJCYm4sCBAyKrqTMbGxvY2Ng0dDPqpbVPG/i1F5tkpadegCTLiPtjo7A6Hv6fibB3VWd1UG2unj+J9KhNaO7XE3bu6g25FZzcwSEkIjIKRVv1MjSGuROWwISGhmLr1q2IioqCt/dfcxE8PDxQVlaG/Px8vV6YnJwceHh46MocPnxYL15OTo7us+r/Vl+rWcbR0bHW3pfGqqyiArIsY1boRKPV6dZrGKybq7tbbrEmFddOH0IzV084eLZVNbZePVeyhcW2dG2v2tLsmkoBgAkMEd3CSbx1o3oCoygKpkyZgs2bN2Pv3r3w8/PT+7xnz56wsrLC7t27MXLkSABAcnIy0tPTERgYCAAIDAzEwoULkZubCze3qmGAyMhIODo6onPnzroyv//+u17syMhIXQxzYW1pCa1WC9ceQ2HpIG4LfgC4qbmAgvNHkBv3h5gKJBkJK++cg0RERPSgVE9gQkJCEBERgV9//RXNmzfXzVlxcnKCnZ0dnJycMGHCBISFhcHFxQWOjo6YMmUKAgMD0bdvXwDAkCFD0LlzZ7z88stYtGgRNBoNZs+ejZCQEN0Q0KRJk/DVV19h+vTpePXVV7Fnzx5s2LAB27ZtU/uWGlRZRQUgybhyzHiH8AWMmASHVupOar589jjO//mzqnNTalOcdRr5x/9AaZ56e7aUFfJcHyIyHk7irRvVE5jly5cDAJ566im966tWrdJtMrd48WLIsoyRI0fqbWRXzcLCAlu3bsXkyZMRGBgIe3t7jB8/HvPnz9eV8fPzw7Zt2zBt2jQsXboU3t7e+P777x94DxhTZ21pCSha4funAH/toaLmQYjVii5nAYCqy5trU5qfLWwVkrb4CipUj8pVSESkj8uo60bIENL92NraYtmyZVi2bNldy7Rp0+aOIaLbPfXUUzh2zPDtys9eyoa9nWlO7k3VXAYAIecT3a5qCbKM6G9UPgixmqDEojZqnlt0MzsJRWf2cBUSEZEJ4VlIAF5fvKahm3BPos8nut2g8W/C2UPdXpKMpKM4svVHeD45DjbO4s7hKcpIwpWj21Q/t6jozG5YePWGZO2gWsxq2qIcaK+cUj0uETVSKvTAgD0wTcO7UybCv73f/Qs2gKjYOHwfscmodXZ4bAC8OjyietwjW3+EU/veaOYpdkfhK0fFzIOSnXwg26u7P041JjBEVE2rKJAMXEWk5SqkpiFoQF8M6NOroZtxV99H/AyPPs/C2rGl0HqKs1NwNfEArqSnqB47XyPmIEQiImqamMCYuPLyckCSoYk1zuoqSZbx8yfvCItfciVdWGwAKM3XCI1PRCSaoqgwiZc9ME3DuZSLcDDR06gzNblGX4XUsnswrFTec+ZGzgUUnj+CtF8/VTUuEZG54SqkumECA+DfM0x/czXjrUKScDVhh7A6bNsPgixgO/5qFdcyUHbpiLD4RESiabWAZPA+MCo1xoQxgQEwM+wN+HcQO7G0vqIOxWD1uvVGXYXk8thoWKq82VxJ9hkUJu2AbNcCFg7qn+hcTcs9VYiImgQmMADaeLc22QQmLT0DWq0WIdPfg7eP4NOoj8Rgw9qVsHJ0g1ULdZdRl+dnAZBw4+TPqsYlIjI3PAupbpjAAJgcNqOhm3BfrX3awK+D2CGkiykXAElCzq6lwupwCPgfWNqL64EpvXIeN89HCYtPRCQaT6OuGyYwjYAsy3h3ivFOoxY5hGTrGQBrF7GTkZnAEBGZPyYwAKa+PQPtHhbbu1FfcTExWLPyO6OcRl2Sm4r85Gg0a/MobFu1Uz1+YZK4ycFEROZCq1VUmMRb9+/PmzcPH3ygv5ilY8eOOHPmjEFtEI0JDIDW3j7o8HDHhm5GrVIvnDf6adTl19TfdK6CJzoTEdVJQyyjfuSRR7Br1y7de0tL008PTL+FRvDO1NCGbsJ9BY4JhZN7a6F1XEqMQ+KuTULnwBARkfEUFhbqvbexsYGNzZ2HF1taWsLDQ9w5dSIwgWkk/Ho+Afd2nYXXkxj5M1wDx8DKSd3fyDczTyH/pPF6kYiIGis1e2B8fHz0rr///vuYN2/eHeXPnTsHLy8v2NraIjAwEOHh4fD1FTtf0VBMYACMDwlD2/YdGroZtToRdxib161C3iX1zye6XUFOJgDAwa8XbN3UnwPDBIaI6P7UPMwxIyMDjo6Ouuu19b706dMHq1evRseOHZGdnY0PPvgATzzxBBITE9G8eXOD2iESExgAHt7eaNPeNCfxpqelQJJkbF88q6GbQkREjYyjo6NeAlObYcOG6f6/a9eu6NOnD9q0aYMNGzZgwoQJoptYb0xgAHwyK6yhm3BfLt2GqH4+0e1u5qah8FyM0DqIiOjeGvosJGdnZzz88MM4f/68QW0QjQkMgKDhL8LNS92dZ9WSciYRMX/uQt7xnQ3dFCIiMoKGPo26qKgIFy5cwMsvv2xQG0RjAgNg1y8bGroJ9/XC6++glZfP/Qsa4OzxOOzauFZoHUREdG+KVnmgfVzuFqOu3n77bTz33HNo06YNsrKy8P7778PCwgKjR482qA2iMYEBMOqlcfDxFXvOUH2dSjyJbVt+Qbd+A+HXKUB4fUxgiIialkuXLmH06NG4evUqWrVqhf79+yMmJgatWrVq6KbdExMYAD169kZAt+4N3Yy72rbll4ZuAhERGYmxD3Ncv369QXU1FCYwAGZMm9LQTSAiIgLQ8JN4GwsmMABCw6ajveCTnusr7nAM/rPq+4ZuBhERkUlhAgPA29sH7U30LKSM9IsN3QQiIjIirVYBjHiYY2PFBAbAzDAOIRERkWlQtJVQtJUGxzB3ckM3gIiIiOhBsQcGwL/eeBt+JjoHJuFIDDauXdnQzSAiIiNhD0zdMIEB0HfAQPTs+3hDN+OumMAQETUdilarQgKjVak1posJDICLF86hmb19QzejVpmcxEtERHQHJjAAPpwxtaGbQEREBABQKiuhVBrYA2Pg9xsDJjAA+kpOaClbNXQzapWpLUGCUtTQzSAiIiNRFBXmwChMYJqEtrIdfGW7hm7GXSVUMoEhImoqOIm3briMmoiIiBod9sAQERGZEPbA1A0TGCIiIhPCBKZuOIREREREjQ57YIiIiEwIN7KrGyYwREREJkSrrQQMTGC0HEIiIiIiMj3sgSEiIjIhnMRbN0xgiIiITAgTmLrhEBIRERE1OuyBISIiMiWVlVBkA3tQeJgjERERGZOiGL4KiYc5EhERkVEpWq3hCUwT2AeGc2CIiIio0WEPDBERkQlRVNjIrimsQmICQ0REZEKqhpAMGwLiEBIRERGRCWIPDBERkQnhEFLdMIEhIiIyIUxg6oZDSERERNTosAeGiIjIhGi1lZDYA3NfTGCIiIhMiFKpBSQDE5hKrkIiIiIiMjnsgSEiIjIhPAupbpjAEBERmRBFW2n4EBLnwBAREZExMYGpG86BISIiokanSffAKIoCAMjVljZwS+7uqrYcAJB2+iRKbtwQWldW2nkAwM3c89CWl6gau/RqBgCgPO8itBXinnd5YTYAoOTKBdXuoawgEwCgvXEFSmW5KjFr0pYUALj1bMrFPZuKW8+msjAbSmWZsHoAQFt8parOgmwoFWLrqlZ5q87K6zlC7q/yRh4AoCI/0yj3VFF4uare4stCft8BgPbGtVv/vQpJ6O/tdKF/7gGgvFBTVaeoe7lZ9etf/XNDJKW8xPAeFEG/Z0yJpBjjV8NEpaSkoF27dg3dDCIiaiQuXLiAhx56SEjskpIS+Pn5QaPRqBLPw8MDqampsLW1VSWeqWnSCUx+fj5atGiB9PR0ODk5NXRzGpXCwkL4+PggIyMDjo6ODd2cRoXPzjB8fvXHZ1d/BQUF8PX1xbVr1+Ds7CysnpKSEpSVqdPDZ21tbbbJC9DEh5BkuWoKkJOTE/8w15OjoyOfXT3x2RmGz6/++Ozqr/rnhii2trZmnXSoiZN4iYiIqNFhAkNERESNTpNOYGxsbPD+++/DxsamoZvS6PDZ1R+fnWH4/OqPz67++OxMT5OexEtERESNU5PugSEiIqLGiQkMERERNTpMYIiIiKjRYQJDREREjQ4TGCIiImp0mmwCs2zZMrRt2xa2trbo06cPDh8+3NBNMjnh4eHo3bs3mjdvDjc3NwwfPhzJycl6ZUpKShASEoKWLVvCwcEBI0eORE5OTgO12HR9/PHHkCQJU6dO1V3js7u3zMxMvPTSS2jZsiXs7OwQEBCAuLg43eeKomDu3Lnw9PSEnZ0dgoKCcO7cuQZssWmorKzEnDlz4OfnBzs7O7Rr1w4ffvih3iGEfHZ/iYqKwnPPPQcvLy9IkoRffvlF7/O6PKu8vDyMHTsWjo6OcHZ2xoQJE1BUVGTEu2iilCZo/fr1irW1tbJy5UolKSlJee211xRnZ2clJyenoZtmUoKDg5VVq1YpiYmJSkJCgvLMM88ovr6+SlFRka7MpEmTFB8fH2X37t1KXFyc0rdvX+Xxxx9vwFabnsOHDytt27ZVunbtqrz55pu663x2d5eXl6e0adNGeeWVV5TY2FglJSVF2bFjh3L+/HldmY8//lhxcnJSfvnlF+X48ePK//zP/yh+fn7KzZs3G7DlDW/hwoVKy5Ytla1btyqpqanKxo0bFQcHB2Xp0qW6Mnx2f/n999+V9957T9m0aZMCQNm8ebPe53V5VkOHDlW6deumxMTEKPv371fat2+vjB492sh30vQ0yQTmscceU0JCQnTvKysrFS8vLyU8PLwBW2X6cnNzFQDKvn37FEVRlPz8fMXKykrZuHGjrszp06cVAEp0dHRDNdOkXL9+XenQoYMSGRmpPPnkk7oEhs/u3mbMmKH079//rp9rtVrFw8ND+fTTT3XX8vPzFRsbG+XHH380RhNN1rPPPqu8+uqretdGjBihjB07VlEUPrt7uT2BqcuzOnXqlAJAOXLkiK7MH3/8oUiSpGRmZhqt7U1RkxtCKisrQ3x8PIKCgnTXZFlGUFAQoqOjG7Blpq+goAAA4OLiAgCIj49HeXm53rP09/eHr68vn+UtISEhePbZZ/WeEcBndz9btmxBr1698MILL8DNzQ09evTAd999p/s8NTUVGo1G7/k5OTmhT58+Tf75Pf7449i9ezfOnj0LADh+/DgOHDiAYcOGAeCzexB1eVbR0dFwdnZGr169dGWCgoIgyzJiY2ON3uampMmdRn3lyhVUVlbC3d1d77q7uzvOnDnTQK0yfVqtFlOnTkW/fv3QpUsXAIBGo4G1tfUdR8u7u7tDo9E0QCtNy/r163H06FEcOXLkjs/47O4tJSUFy5cvR1hYGN59910cOXIEb7zxBqytrTF+/HjdM6rtz3FTf34zZ85EYWEh/P39YWFhgcrKSixcuBBjx44FAD67B1CXZ6XRaODm5qb3uaWlJVxcXPg8BWtyCQzVT0hICBITE3HgwIGGbkqjkJGRgTfffBORkZGwtbVt6OY0OlqtFr169cJHH30EAOjRowcSExOxYsUKjB8/voFbZ9o2bNiAdevWISIiAo888ggSEhIwdepUeHl58dmRWWlyQ0iurq6wsLC4Y7VHTk4OPDw8GqhVpi00NBRbt27Fn3/+CW9vb911Dw8PlJWVIT8/X688n2XVEFFubi4effRRWFpawtLSEvv27cMXX3wBS0tLuLu789ndg6enJzp37qx3rVOnTkhPTwcA3TPin+M7vfPOO5g5cyZGjRqFgIAAvPzyy5g2bRrCw8MB8Nk9iLo8Kw8PD+Tm5up9XlFRgby8PD5PwZpcAmNtbY2ePXti9+7dumtarRa7d+9GYGBgA7bM9CiKgtDQUGzevBl79uyBn5+f3uc9e/aElZWV3rNMTk5Genp6k3+WgwcPxsmTJ5GQkKB79erVC2PHjtX9P5/d3fXr1++OJftnz55FmzZtAAB+fn7w8PDQe36FhYWIjY1t8s/vxo0bkGX9v9otLCyg1WoB8Nk9iLo8q8DAQOTn5yM+Pl5XZs+ePdBqtejTp4/R29ykNPQs4oawfv16xcbGRlm9erVy6tQpZeLEiYqzs7Oi0WgaumkmZfLkyYqTk5Oyd+9eJTs7W/e6ceOGrsykSZMUX19fZc+ePUpcXJwSGBioBAYGNmCrTVfNVUiKwmd3L4cPH1YsLS2VhQsXKufOnVPWrVunNGvWTPnPf/6jK/Pxxx8rzs7Oyq+//qqcOHFCef7555vsUuCaxo8fr7Ru3Vq3jHrTpk2Kq6urMn36dF0ZPru/XL9+XTl27Jhy7NgxBYDy+eefK8eOHVMuXryoKErdntXQoUOVHj16KLGxscqBAweUDh06cBm1ETTJBEZRFOXLL79UfH19FWtra+Wxxx5TYmJiGrpJJgdAra9Vq1bpyty8eVN5/fXXlRYtWijNmjVT/v73vyvZ2dkN12gTdnsCw2d3b7/99pvSpUsXxcbGRvH391e+/fZbvc+1Wq0yZ84cxd3dXbGxsVEGDx6sJCcnN1BrTUdhYaHy5ptvKr6+voqtra3y0EMPKe+9955SWlqqK8Nn95c///yz1r/nxo8fryhK3Z7V1atXldGjRysODg6Ko6Oj8s9//lO5fv16A9xN0yIpSo3tGYmIiIgagSY3B4aIiIgaPyYwRERE1OgwgSEiIqJGhwkMERERNTpMYIiIiKjRYQJDREREjQ4TGCIiImp0mMAQERFRo8MEhoiIiBodJjBERETU6DCBISIiokbn/wGADujbnK9njwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "key = apt.get_key()\n",
    "\n",
    "key, hist = llh._simulate_model_hist(key, int(1e6), parameters)\n",
    "apt.plot_irreg_histogram_2d(*llh._bins, hist)\n",
    "plt.show()\n",
    "\n",
    "key, log_llh = llh.get_log_likelihood(key, int(1e6), parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ebfaec0b-b531-4ba8-94ab-1174fd3d263a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The log (posterior) LLH is 1978.69 now.\n"
     ]
    }
   ],
   "source": [
    "print(f\"The log (posterior) LLH is {log_llh:.2f} now.\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}