{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to use a FIGARO reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FIGARO, despite being a stochastic sampler, does not produce a set of samples (numbers) directly, but rather a number of `figaro.mixture.mixture` objects stored as `.json` or `.pkl` files. Each of these objects represents a probability density drawn around the probability density that generated the available data. If you need to use these realisations in a `python` script (to evaluate their pdf, for example, or to draw realisations from them), here you'll find the (few!) steps needed for doing so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from figaro import plot_settings\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic usage "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us assume to have our realisations stored in a file called `draws_file.json` (they can be the product of a hierarchical inference or a DPGMM reconstruction, it does not make any difference). The first thing to do is to load the realisations via the dedicated method, `figaro.load.load_density`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from figaro.load import load_density\n",
    "draws = load_density('./draws_example.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`draws` will be a list of `figaro.mixture.mixture` objects, each of them representing a probability distribution. The methods of this class are modelled after the `scipy.stats` methods to facilitate users that are already familiar with the SciPy package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.03470811]),\n",
       " array([-3.36078192]),\n",
       " array([0.98599125]),\n",
       " array([[-0.214056  ],\n",
       "        [ 0.09627442],\n",
       "        [-1.66075543]]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = draws[0]\n",
    "X = 2.\n",
    "(d.pdf(X), d.logpdf(X), d.cdf(X), d.rvs(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced use"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you need to draw samples from the median distribution, there is a dedicated method in the `figaro.utils` module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiEAAAGbCAYAAAASrkAJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfIUlEQVR4nO3de3wcd33v/9fMWnIkWdJKsuWLfJElOXFsixibcAkhScHiXtpSJ4EWejkFh3LaQ9q0VkMPv0LPOU1lDoUWerGBXikltmnhNNAWi9KEcHcUJ7IVObpYtiVLliXtru63mfn9sV5Fim47392Z2Zn9PB8PPzLe23znvZ/IH33nplmWZSGEEEII4TLd6wEIIYQQIjtJEyKEEEIIT0gTIoQQQghPSBMihBBCCE9IEyKEEEIIT0gTIoQQQghPSBMihBBCCE+s8XoAyzFNk2vXrlFYWIimaV4PRwghhBBJsCyLkZERtmzZgq6vPNeRsU3ItWvX2LZtm9fDEEIIIYSCq1evsnXr1hVfk7FNSGFhIRDfiKKiIo9H4z3TNLl69Srbtm1btbMU6iRnd0jO7pCc3SE5LzQ8PMy2bdvm/h1fScY2IYldMEVFRdKEEC/yjRs3UlRUJEXuIMnZHZKzOyRnd0jOS0vmUAotU+8dMzw8THFxMbFYTJoQIYQQwifs/PstLZtPGIZBR0cHhmF4PZRAk5zdITm7Q3J2h+SsTpoQn9A0jXA4LGcKOUxydofk7A7J2R2Ss7qMPSZELKTrOmVlZV4PI/AkZ3dIzu6QnN0hOauTmRCfMAyDtrY2me5zmOTsDsnZHZKzOyRnddKE+ISu62zYsEGOvHaY5OwOydkdkrM7JGd1sjvGJxL7HIWzJGd3SM7ukJzdITmrk7bNJwzDoLW1Vab7HCY5u0Nydofk7A7JWZ00IT6h63pS1+EXqZGc3SE5u0NydofkrE52x/iEpmly0TYXSM7ukJzdITm7Q3JWJ22bTxiGwYULF2S6z2GSszskZ3dIzu6QnNVJE+ITuq5TWVkp030Ok5zdITm7Q3J2h+SsTnbH+ISmaRQUFHg9jMCTnN0hObtDcnaH5KxO2jafMAyD5uZmme5zmOTsDsnZHZKzOyRndTIT4hO6rlNTUyPTfQ6TnN3x8px7ohNExqZXfE9JQS4V4Tw3hhcYUs/ukJzVSRPiE5qmkZcnP4CdJjm7Y37OPdEJDn3qSSZmVv4tMi8nROMj90ojYoPUszskZ3XShPhEYrqvtraWUCjk9XACS3J2x/ycI2PTTMwYfObB/dSUr1vy9e39ozz8+DkiY9PShNgg9ewOyVmdNCE+oes6e/bskek+h0nO7lgq55rydeyrKPZwVMEj9ewOyVmdJOYj0mG7Q3J2h+TsDsnZHZKzGmlCfMI0TZqbmzFN0+uhBJrk7A7J2R2SszskZ3XShPiEruvU1tbKdJ/DJGd3SM7ukJzdITmrk8R8RM5Bd4fk7A7J2R2SszskZzXShPiEaZq0tLTIdJ/DJGd3SM7ukJzdITmrk7NjfCIUCrF//36vhxF4krM7JGd3SM7ukJzVyUyIT1iWxcTEBJZleT2UQJOc3SE5u0NydofkrE6aEJ8wTZP29naZ7nOY5OwOydkdkrM7JGd1sjvGJ0KhELW1tV4PI/AkZ3dIzu6QnN0hOauTmRCfsCyLsbExme5zmOTsDsnZHZKzOyRnddKE+IRpmnR1dcl0n8MkZ3dIzu6QnN0hOauT3TE+EQqF2Lt3r9fDCDzJ2R2SszskZ3dIzupkJsQnLMtieHhYpvscJjm7Q3J2h+TsDslZnTQhPmGaJteuXZPpPodJzu6QnN0hObtDclYnu2N8IhQKsXv3bq+HEXiSszskZ3dIzu6QnNXJTIhPWJZFNBqV6T6HSc7ukJzdITm7Q3JWJ02IT5imyY0bN2S6z2GSszskZ3dIzu6QnNXJ7hifCIVC7Nq1y+thBJ7k7A7J2R2SszskZ3UyE+ITpmkyODgonbbDJGd3SM7ukJzdITmrkybEJ2SfozskZ3dIzu6QnN0hOauT3TE+EQqFqK6u9noYgSc5u0Nydofk7A7JWZ3MhPiEaZr09/fLdJ/DJOc0Gb4G/3g/jN5Y8mnJ2R2SszskZ3XShPjI+Pi410PICpJzioYuwfF7oO1b0PHtZV8mObtDcnaH5KxGdsf4hK7rVFZWej2MwJOc0+CFf4Xa+6F8D1x6Cu54z6KXSM7ukJzdITmrk5kQnzBNk76+Ppnuc5jknAYv/gfsejPsvCfehCxxsN6yOUevwvSYSwMNPqlnd0jO6qQJ8ZHp6Wmvh5AVJOcUTETh+nnY8Xoo2QHmLIz2L/nSRTmbs/Cnd8Bfv3XJxkWokXp2h+SsRpoQn9B1ne3bt6Pr8pU5SXJOUce3YecbYE1u/O9lNTDYvuhlL89Zw2TL0x+FW98Ks1PQ1+zmqANL6tkdkrM6OSbEJ0zTpLe3l82bN0uhO0hyTtGL34Jdb3np72U1RLpfoCdn34KXmabJwMAA69ev53LfIB1r38/48Gvhv/0L/MfvQ/ePYfMrFn18e//osqsuKcilIpyXtk0JAqlnd0jO6qQJEUKkh2VBeyPUfWLuoWhhNU986wwfm968zJte5D79HLtzKij42b9nXU4e3PqWeCNy4JchlAPEG4y8nBAPP35u2dXn5YRofOReaUSE8BFpQnxC13UqKiq8HkbgSc4pGB8EfQ0Ubpp7aKDkleznC3zmwf3UlK9b8m3lz/yYW9b+IkUby+MP3PoW+M7/gesXYMt+ACrCeTQ+ci+RsaX3u7f3j/Lw4+eIjE1LEzKP1LM7JGd10oT4hGmadHd3s3XrVpnuc5DknILoFQhvW/DQZOntbNKGyCkYZfe8H9ILcp7tg511Cz9r653Q/ZO5JgTijYg0GPZIPbtDclYnaflIbm6u10PICpKzolg3FG9d+Ji+hm8aryHc/rVFL5/LOdIF4R0LnyzfAzdaHRlmtpF6dofkrEaaEJ/QdZ1NmzZJl+0wyTkFsatQvG3Rw/9ivIHi9n9Z8NhczpoGQx1QUrnwTWXVMNjh4GCzg9SzOyRndZKYT5imSVdXl1wMx2GScwpi3Us2IeesavTZyfgxHjfN5dx9FvJKYN2GhW8qq5EmJA2knt0hOauTJsRH8vPzvR5CVpCcFUWvLN4dA4DGyPY3QueTCx7Nz8tD+9qH4DUPLX5LUQXMjMHYgDNjzSJSz+6QnNVIE+ITuq5TXl4u030Ok5xTEOtedGBqwnj5K6Hl6zAZg9kp9Ka/pfwLr0SLXoU7P7D4DZoG214DV37o8KCDTerZHZKzOknMJwzDoKOjA8MwvB5KoEnOKYhdXWYmBEa2H4K1hfDH2+G7n4InfovRza/B+M2m5T9v++vgyg8cGmx2kHp2h+SsTpoQn9A0jXA4jKZpXg8l0CRnRdPjMDsNt4SXfNrKyYf3fgV2vxOebMB82zGm3vE5tKIty3/mjrtkJiRFUs/ukJzVSRPiE7quU1ZWJtN9DpOcFSVOz13ph3BoDdz5awDo+36eso0VK+dcfjsMvCg3s0uB1LM7JGd1kphPGIZBW1ubTPc5THJWtMKumAWq3whHL2HcUrJ6zrkFsLYIRq+nb5xZRurZHZKzOmlCfELXdTZs2CCdtsMkZ0XD12ClXSvz5Zcmn3NZNQy0pT6+LCX17A7JWZ0k5hOyz9EdkrOi0esL7hmzmqRzLquBwfYUB5e9pJ7dITmrkybEJwzDoLW1Vab7HCY5Kxq7AevKk3550jmv3yVNSAqknt0hOauTJsQndF1ny5YtMt3nMMlZ0eh1KEi+CUk6Z5kJSYnUszskZ3W276Lb2dlJQ0MD1dXVAITDYY4cOWJ7xXV1dZw5c8b2+7KVpmkUFRV5PYzAk5wVjd6AdRuTfnnSOZdVSxOSAqlnd0jO6my1bZ2dnRw8eJCGhgaOHj3K0aNH6ejo4NixY7ZWeuzYMRobG229J9sZhsGFCxdkus9hkrOi0eu2d8cklXN4Bwz3gjGT4gCzk9SzOyRndbaakIaGBo4cOUI4HJ577NFHH6W+vj7pz+js7OQnP/mJndUK4tN9lZWVMt3nMMlZ0Wi/rSYk6Zz1UPzU38jlFAeYnaSe3SE5q7OV2MmTJ+d2wyQkGpJkZzZOnz7Ngw8+aGe1gvh0X0FBgRx97TDJWcHMJFhG/LoeSbKVsxwXokzq2R2Ss7qkm5BoNEo0GqWqqmrRc+FwmKamFe4BcdPp06c5fPiwvREKID7d19zcLNN9DpOcFYzZmwUBmzmvr4FBuVaICqlnd0jO6pJuQjo7O5d9rrS0lMHBwRXfH41GGRoaWrKJWe19w8PDc03Q1NQUpmlimibAgmXDMBYsWzcv92xn2bKsRcuA7WXTNBcsLzXeZJYT26TrOlVVVXOddhC2KRO/J4CampoFf/f7Njn+PY3Ez4xZbpsS45s/Rk3T2LlzJ7qur75NZTVYA22rblMy31m2fU+aplFTU7Pgu/D7NmXi9wRQVVWFruuB2aZUv6dkpW0HVjQaXfH5EydOKJ1Fs2PHDoqLiykpKaGkpITHHnuM7u5u+vv7Abhy5QoDAwMAdHV1EYlEgHjTFIvFAGhra2NkZASA1tZWxsfHAWhpaWFychKA5uZmZmZmME2T5uZmTNNkZmaG5uZmACYnJ2lpaQFgfHyc1tZWAEZGRmhri/+WFovF5pq1SCRCV1cXAAMDA1y5cgWA/v5+uru7Aejt7aW3txdg1W3SNI1r164xPDwcmG3KxO/p6tWr5OXlcePGjcBsk9Pf01S0F/LLlt0mYNE2TU1NcenSJTRNW32bymqYvd666jYlnpPv6aVtGh0dJS8vj4sXLwZmmzLxe+rr6yMajaJpWmC2KZXvqaOjg6RZSero6LAA68yZM4ueC4fD1pEjR5Z975kzZ6xnnnlm7u+nTp2yVlt1LBazAOvy5ctWLBazIpGIFYlErMnJScswDMswDMuyrAXLs7OzC5ZN07S9bJrmomXLsmwvG4axYHmp8SaznNim2dlZq6mpyZqZmQnMNmXi9zQ9PW09++yz1vT0dGC2yfHv6bmTlvXPDy35muevRqwd9U9Yz1+NLBjjzMyM1dTUNPe6FbdppN8yP3nritvU3B21dtQ/YT13ZUi+p3njnZmZsZ599llramoqMNuUid/T9PT0XD0HZZtS+Z6GhoYswIrFYtZqkr5OSGlpKbD0jEc0Gl1wxszLNTU1cfTo0eQ7o3nC4fCK51/PPxo5FAo5tqxpmq3l+eNKZTnxeZZlsXfv3rm/B2GbnFhOdZs0TWPPnj2sWbNmbteX37dJdTnpMU4Nw9qiZcee+O/8MYZCIfbu3TuX+YrbVLAebWaC0Ow4hAqXHe/8v8v3FF+2LIs9e/aQk5Mz9134fZtSXXZim9asWbOgnoOwTen6nlaTdBMSDocJh8MMDQ0t+XxdXd2Sj584cYKOjo4Fp/EmDmKtr6+nrKxMuUHJNna+WKFOcrZpahhusX+hpqRz1rSXLlq25ZW215PtpJ7dITmrsXXF1AceeGDRvp7EPqZDhw4t+Z6ljgM5ceIEjY2NNDQ02Fl9VkvsC6ytrZVid5DkbN9IbIhRvZDBntii59r7R5d8j+2c1++CwQ5pQmySenaH5KzOVhNSX19PXV3dgubh+PHjHD9+fO7v0WiU+++/n4aGBg4cOLDk56x2EKtYTNd1amtr5WI4DpOc7emJTvDkj1p53tjOV558esnX5OWEKCnIXfCY7ZzLamBATtO1S+rZHZKzOltNSFVVFadOnaK+vp4777yTzs5OysrKFsx2DA0Ncfbs2SV323R2dnL8+HFOnz4NwP33309dXZ3SWTPZyDAMKXIXSM7Ji4xNk2+N8u7X3c77Xnn3kq8pKcilIpy36HFbOZffDs99JZWhZi2pZ3dIzmps38DuwIEDy85wQLxRSZyys9RzDQ0NshtGgWmatLS0yHSfwyRn+wqZoHx9OZUVxUm/x3bOm14B//HRFEaZnaSe3SE5q5O2zSdCoRD79++XAneY5GxfWBvFuCVs6z22cw5vh5wCePE/7A8wi0k9u0NyVidNiE9YlsXExMTcFemEMyRn+8KMMrs2bOs9tnPWNHj1B6DtW/YHmMWknt0hOauTJsQnTNOkvb19wWWCRfpJzvaFtVEMm02IUs4ba6Gv2d7gspzUszskZ3W2jwkR3giFQtTW1no9jMCTnG2yTNYxSW+uveuEKOW8cQ/0t4JlxWdGxKqknt0hOauTmRCfsCyLsbExme5zmORsjz49wgh5tpsCpZzXFkLOLTA2YHOU2Uvq2R2SszppQnzCNE26urpkus9hkrM9a6YixKwC2+9TzrmsBoZs3Bwry0k9u0NyVidNiE8k7rUhR187S3K2JzQZJUKh/fep5lxaFb98u0iK1LM7JGd10oT4hGVZDA8Py3SfwyRne0JTUaIKMyHKOZfVxC/fLpIi9ewOyVmdNCE+YZom165dk+k+h0nO9oSmIkQVZkKUc07cyE4kRerZHZKzOjk7xidCoRC7d+/2ehiBJznbE5qMErHW2X+fas5lNTDUaf99WUrq2R2SszqZCfEJy7KIRqMy3ecwydmeNVMRogpNiHLOJTth6BLIb5xJkXp2h+SsTmZCfMI0TW7cuEFhYaEc/OQgydme0GSUKPabEOWcc26B/DIY6YXiikVPt/ePrvj25W6mF1RSz+6QnNVJE+IToVCIXbt2eT2MwJOc7QlNRYlY2+2/L5Wc19fAjdYFTUhJQS55OSEefvzcim/NywnR+Mi9WdOISD27Q3JWJ02IT5imSSQSoaSkRG4X7SDJ2Z74gal7bL8vpZy3vQau/ghq3jT3UEU4j8ZH7iUyNr3s29r7R3n48XNExqazpgmRenaH5KxOmhCfSOxzDIfDXg8l0CRne0KTqR0TopTzttfA059e9HBFOC9rmotkST27Q3JWJ02IT4RCIaqrq70eRuBJzvasmYoqNSEp5bz5Drh+Xu4hkwSpZ3dIzupk3sgnTNOkv79fzkN3mORsT2gqSkTxwFTlnPNLYc0tMNxj/71ZRurZHZKzOmlCfGR8fNzrIWQFyTlJxiyaMcUEa5XenlLOlW+A9m+rvz+LSD27Q3JWI02IT+i6TmVlpRz05DDJ2YbJKMbaYsD+LpGUc655E3Q9rfbeLCL17A7JWZ0k5hOmadLX1yfTfQ6TnG0YH8JYW6L01pRzLquGyCW192YRqWd3SM7qpAnxkenp5U8/FOkjOSdpIoJxi1oTAinmnLhyqliV1LM7JGc1cnaMT+i6zvbt9i8KJeyRnG2YGGJ2bVjprSnnnF8KlgETUchTG0M2kHp2h+SsTmZCfMI0TXp6emS6z2GSsw0TEQzFJiQtOW95JfQ8o/7+LCD17A7JWZ00IUIINeNDyk1IWmx7jRycKoTPSRPiE7quU1FRIUdfO0xytmFiSPmYkLTkXHs/PPcVMGbUPyPgpJ7dITmrk8R8wjRNrly5ItN9DpOcbZiIpHR2TMo5l1VD+W64+G/qnxFwUs/ukJzVSRPiI7m5uV4PIStIzkkaVz8wFdKU821vh87/Sv1zAkzq2R2SsxppQnxC13U2bdok030Ok5xtSHF3TFpy3v66eBMyO5Xa5wSU1LM7JGd1kphPmKZJV1eXTPc5THK2IcWzY9KS88a9ULQF/uoNqX1OQEk9u0NyVidNiI/k5+d7PYSsIDknaTyCcUtY+e1pyVnT4H1fjd/MblYuFrUUqWd3SM5qpAnxCV3XKS8vl+k+h0nONqRwYGpac16zFkqroP9C6p8VMFLP7pCc1UliPmEYBh0dHRiG4fVQAk1yTtLsFGgaVkjtYLy057zzHjlAdQlSz+6QnNVJE+ITmqYRDofRNPt3LBXJk5yTND4Eeer3jUl7znt+Bpr+AaZG0vN5ASH17A7JWZ00IT6h6zplZWUy3ecwyTlJE5GUmpC057zt1XD7O+FP9sLYQHo+MwCknt0hOauTxHzCMAza2tpkus9hknOSJlKbCXEk57o/hH3vhh9/Pn2f6XNSz+6QnNVJE+ITuq6zYcMG6bQdJjknaSISv5OtIsdyvu/34Jm/gb7z6f1cn5J6dofkrE4S8wnZ5+gOyTlJ40OQp96EOJZz4SZ4yx/BEw+n93N9SurZHZKzOmlCfMIwDFpbW2W6z2GSc5LSsDvGsZz3/Xy8SZLZEKlnl0jO6qQJ8Qld19myZYtM9zlMck5SGnbHOJazpsGuOrj8vfR/ts9IPbtDclYnifmEpmkUFRXJdJ/DJOckpWF3jKM5bzkAPU3OfLaPSD27Q3JWJ02ITxiGwYULF2S6z2GSc5JSPEXX8ZwrDsA1aUKknt0hOauTJsQndF2nsrJSpvscJjknKQ27YxzNubQaRq7DZMyZz/cJqWd3SM7qJDGf0DSNgoICme5zmOScpDRcMdXRnHUdNtXC9ey+n4zUszskZ3XShPiEYRg0NzfLdJ/DJOckTURSOibElZw33Ao3Ljr3+T4g9ewOyVmdNCE+oes6NTU1Mt3nMMk5CZYFk1HICyt/hCs5r78NBl507vN9QOrZHZKzOknMJzRNIy8vT6b7HCY5J2F6DNasBT2k/BGu5CwzIVLPLpGc1UkT4hOGYXDu3DmZ7nOY5JyEFHfFgEs5r78NBtqc+3wfkHp2h+SsTpoQn9B1nT179sh0n8Mk5yRMDKV0Zgy4lHPRlvjZMdNjzq0jw0k9u0NyVieJ+UgopD79LZInOa8ixWuEJDies6bB+pqs3yUj9ewOyVmNNCE+YZomzc3NmKbp9VACTXJOwtgA5Jel9BGu5bzhdrjR6uw6MpjUszskZ3XShPiEruvU1tbKdJ/DJOckjPRB4eaUPsK1nMt3Q/8Lzq4jg0k9u0NyVieJ+Ygc9OQOyXkVI70pNyHgUs5ZPhMCUs9ukZzVSBPiE6Zp0tLSItN9DpOckzDSC4WbUvoI13Iu353VTYjUszskZ3VrvB6ASE4oFGL//v1eDyPwJOckpGF3jGs5F1XARBSmRp1fVwaSenaH5KxOmhCfsCyLyclJbrnlFrkgjoMk54V6ohNExqYXPLYr0kPXeAEzPTHa+9X+cXctZ02D9bfCwEWgxrn1ZCipZ3dIzuqkCfEJ0zRpb29nz549ciqYgyTnl/REJzj0qSeZmJm/r9uiZW0vP/237UxzGYC8nBAlBbm2PtvVnMt3Q38rlGdfEyL17A7JWZ00IT4RCoWora31ehiBJzm/JDI2zcSMwWce3E9N+ToA9KkYuScL+OcjPzX3upKCXCrCebY+29WcN9wON16A8ne6s74MIvXsDslZnTQhPmFZFuPj4+Tn58t0n4Mk58Vqytexr6I4/pf+Xije8tLfFbmac/lu6Pwv2OvsajKR1LM7JGd1cnaMT5imSVdXlxx97TDJeRWxq/GDPVPkas5ZfJqu1LM7JGd1MhPiE6FQiL17s/BXOZdJzquIdEHpzpQ/xtWci7bARARtdtKd9WUQqWd3SM7qZCbEJyzLYnh4GMuyvB5KoEnOqxi6BCWpNyGu5qxpUFJJ7vBl59eVYaSe3SE5q7M9E9LZ2UlDQwPV1dUAhMNhjhw5sur7GhsbaWpqAqCjo4Pq6mqOHj1qd/VZyzRNrl27xq5du+ToawdJzqsY6oCq+1L+GNdzLq0iN3YJWOf8ujKI1LM7JGd1tpqQzs5ODh48yKVLlwiHwwDU19dz7NixFRuKpqYmotHogtdUV1fT0dHB8ePH1UaeZUKhELt37/Z6GIEnOa/iegts3JPyx7iec1k1a4cvAaufwbDUtVHmUzkbyCtSz+6QnNXZakIaGho4cuTIXAMC8Oijj1JSUrJiE3L8+HEaGxs5fPjw3GOHDh3ixIkT0oQkybIsYrEYxcXFcvS1gyTnFUzGYGo4LQemup5zaTW5F59mtSZk6WujLJSXE6LxkXt90YhIPbtDclZnqwk5efIkDQ0NCx5LNCSNjY0cOnRoyffV1dWpjU7MMU2TGzduUFhYKNN9DpKcV3D9AmzcFz/GIkWu51xWzdrhL636sqWujTJfe/8oDz9+jsjYtC+aEKlnd0jO6pJuQqLRKNFolKqqqkXPhcNhmpqalm1CDh8+vGAWBJZuaMTyQqEQu3bt8noYgSc5r+D6BdiYnjMAXM+5rObmMSHJWXBtFB+TenaH5Kwu6bNjOjs7l32utLSUwcHBVT/j9OnT1NfXU1dXx6lTp5I6MDUajTI8PDzXBE1NTWGa5tz52POXDcNYsJw4UtnOsmVZi5YB28umaS5YXmq8ySwntinRaSc+MwjblInf0+zsLIODg8zOzgZmm1L5nhIMw8DqOw8b96ZlmwzD4MaNG5im6c425ZWhz4yTx+Sq31Ni3Xa2z+vvabltMgyDwcFBZmZmfFd7fvr/aXZ2dq6eg7JNqX5PyUrbKbrRaHTV1xw+fJiGhgYeeugh6uvrV2xsEnbs2EFxcTElJSWUlJTw2GOP0d3dTX9/PwBXrlxhYGAAgK6uLiKRCBBvmmKxGABtbW2MjIwA0Nrayvj4OAAtLS1MTsavHdDc3MzMzAymadLc3IxpmszMzNDc3AzA5OQkLS0tAIyPj9PaGr/40cjICG1tbQDEYrG5bYpEInR1dQEwMDDAlStXAOjv76e7uxuA3t5eent7AVbdJsuy6OjomMs5CNuUqd9TNBoN3DapfE/Xbj6f2Cbj2nOwcV/atunixYtzV5p0fJsGB5ku3kGldp3hm2Nf7ntKjG+pbZqdnQWgvb09Y76n1WovGo1y8eJFX9WeH/9/unr1KpZlBWqbVL+njo4OkmYlqaOjwwKsM2fOLHouHA5bR44cSfajLMuyrKNHj1rhcNiKRCJLPh+LxSzAunz5shWLxaxIJGJFIhFrcnLSMgzDMgzDsixrwfLs7OyCZdM0bS+bprlo2bIs28uGYSxYXmq8ySzLNsk2ebVN5y4PWjvqn7Cau6PW7My0Zf7RVsuaGvXtNkX+7hetDz36Meu5K0PLjre5O2rtqH/Cev5qZMnXPH81Yu2of8I6d3kwI7YpMa6g1Z5sk7+3aWhoyAKsWCxmrSbpmZDS0lJg6RmPaDS64IyZZNTV1RGNRjlx4sSKrwuHwxQVFREOhwmHw6xduxZd19H1+NDnL4dCoQXLiaOU7SxrmrZoGbC9rOv6guWlxpvMcmKbTNNkcHBwbrorCNuUid8TMNfxB2WbUvmeEkLDV9EKNkBuQVq2ybIsBgcHMU3TtW2aLt7JTq1v1e8psW472+f197TcNlmWRX9//9xY7X5PmbhNmfj/E8RnFEzTDMw2pfo9JSvpJiTRBAwNDS35/EpnwJSUlHDs2LEFjyWaGlvTNlkuMU0nnCU5LyGNB6UmuJ3zdFEllVqfq+vMBFLP7pCc1dg6RfeBBx5Y1DQk9jEtd2ZMYubk5WfVJN538OBBO0PIWrquU1lZ6fUwAk9yXkbi9Nw08SLnqeKdVOrZ1YRIPbtDclZn68DU+vp6Tp8+veCx48ePL7jgWDQapa6ubu4S7YnLuh84cGDR+w4cOJDUJd9F/Ijkvr6+uSORhTMk52VcP5/WmRAvck7sjskmUs/ukJzV2ZoJqaqq4tSpU9TX13PnnXfS2dlJWVnZgkZiaGiIs2fPLtht09DQwIkTJ3jmmWcIh8N0dnZy4MABuU6ITdPTy19KWqRPEHJe7dLjkPzlx7XZCbh2Dur+V5pGF+d2zsbaEnKYRZ8eAfx/DZBkBaGe/UByVmP7BnYHDhxYNKsxX1VV1dwpO/PJjEdqdF1n+/btXg8j8IKQczKXHofkLz9eduFvYfMdUFKZtjF6krOm0WVtZN1wF7DV3XV7JAj17AeSszrbTYjwhmma9Pb2snnz5gVHZIv0CkLOq116HJa5/PhkDAbboSJ+nJY2O8mb9Z9QeOX7UPdoWi7XnuBVzl3WJg7ELgF3u7ZOLwWhnv1AclYnTYgQAWXr0uOWBf/0C3D5afjpP4PqN1J49WlO5H4a+oAddzk6Vrd0WZt4XazL62EIIW6SJsQndF2noiL1u5eKlQU2Z9OMz2RMxiAvvPC54V44+X7Qc+BD34N/+RB854/YPtrHC+Z2tm/eQMGatWkdjlc5XzI3kTvc5fp6vRLYes4wkrM6aUJ8wjRNuru72bp1q0z3OShwOc9MwBO/Dc99GcI7IHoZ3vdVCvtjfCP3E2w9czswBpVvgPt+D9ashV9/GoCOc0/x3q9c40vvehvpOzk3zqucu60N5Iw+69r6vBa4es5QkrM6aUJ8JDc31+shZIVA5Xz2b2C4B977FbjwNdh5D3zp56lYW8I/mQe5rfcHsO/n4I0fg5f98JzYcAdRRpb+3DTwIuduawO5I92ur9dLgarnDCY5q5EmxCd0XWfTpk1eDyPwApWzacD3/wx+6euw4Ta47W3xmZGSSlqrfo2Pfe77HHzPZ9izc5vrQ/Mq537CrJkcAGMWQsH/8Reoes5gkrM6mTfyCdM06erqkovhOCxIOd8y1AIF6+MNSEJOHtzzO6DF/9c3c4s8GZtXOVvozORvhJFrrq7XK0Gq50wmOauTJsRH8vPzvR5CVghKzgW9P4IdmXsqqlc5zxRuhegVT9bthaDUc6aTnNVIE+ITuq5TXl4uBz05LEg5F/T+ECpf7/UwluRlztPrsqcJCVI9ZzLJWZ0k5hOGYdDR0YFhrHwVTJGaoOS8T+skv78JdmRmE+JlzjOF2yDS5fp6vRCUes50krM6aUJ8QtM0wuEwWhqvWikWC0rOH1jzTa4f+C3IL/V6KEvyMuep4p0w2LH6CwMgKPWc6SRnddKE+ISu65SVlcl0n8MCkbNlcbd+nuGdb/N6JMvyMufp4ioYbHN9vV4IRD37gOSsThLzCcMwaGtrk+k+hwUh55yRK4xbazHy1ns9lGV5mXN8JqQzfhXZgAtCPfuB5KxOmhCf0HWdDRs2SKftsCDknN9/jnNWjdfDWJGXOVtr8qBwEwx1ur5utwWhnv1AclYnifmE7HN0RxByzrvxLM+amd2EeJ7z5jug95w363aR5zlnCclZnTQhPmEYBq2trTLd57Ag5Jzff45zGd6EeJ7z5jug73lv1u0iz3POEpKzOmlCfELXdbZs2SLTfQ7zfc6z06yNvEiLtcPrkazI85w33wG9z3mzbhd5nnOWkJzVSWI+oWkaRUVFMt3nMN/nfL2ZqXA1U2T2zbQ8z3lTLfQ+D5blzfpd4nnOWUJyVidNiE8YhsGFCxdkus9hvs+5+xkmNrzS61GsyvOc80shtwCG1e4h094/yvme2LJ/eqITaR6wGs9zzhKSs7rg30YyIHRdp7KyUqb7HOb7nHvOMl7+Oq9HsaqMyHlTbfzg1OKKpN9SUpBLXk6Ihx8/t+Lr8nJCND5yLxXhvNTGmKKMyDkLSM7qpAnxCU3TKCgo8HoYgef7nLt/wsTuDwM9Xo9kRRmR847Xw6WnYPc7kn5LRTiPxkfuJTI2vexr2vtHefjxc0TGpj1vQjIi5ywgOauTts0nDMOgublZpvsc5uucJ6IwPsR0UaXXI1lVRuRceTdc+YHtt1WE89hXUbzsn5rydQ4MVk1G5JwFJGd10oT4hK7r1NTUyHSfw3yd88CLsGE3+ODguIzIecNtgb9yakbknAUkZ3WSmE9omkZeXp4cfe0wX+d8ozX+D6sPZETOOXnxA1SHu70bg8MyIucsIDmrkybEJwzD4Ny5czLd5zBf53zjom+akIzJecNtcONFb8fgoIzJOeAkZ3VyYKpP6LrOnj17ZLrPYb7O+cZFqP4pW29p7x9Vei5VGZPz+lth4CLsOuTtOBySMTkHnOSsTpoQHwmFQl4PISv4NueBi7D+Nhhb/aV2TjUtKXDmwmcZkfOG26D7J16PwlEZkXMWkJzVSBPiE6Zp0tzcTG1trRS7g3yb8/QYjEegeCuMDa/68mRONYV4s+LEaaYZk/P2u+CpT8avnBrA/fkZk3PASc7qpAnxCV3Xqa2tlek+h/k254E2WF9j6x/SinCeZ9exyJic19dA7rr4fWS27Pd2LA7ImJwDTnJWJ4n5iBz05A5f5nzjYvz0XB/JmJxvfxec/aLXo3BMxuQccJKzGmlCfMI0TVpaWjADfE2DTODbnAcuxg+y9ImMyvl1H4b2b8NAu9cjSbuMyjnAJGd10oT4RCgUYv/+/bK/0WG+zdlnMyEZlfMtxfCKB+HpP/F6JGmXUTkHmOSsTpoQn7Asi4mJCayA33rca77N2UfXCIEMzPneo9D5JLcMNC983DThu5+Cx7bDX94Nl+1f5t1LGZdzQEnO6qQJ8QnTNGlvb5fpPof5MufZ6fgt6cM7vB5J0jIu55w8ePMfsu3b/x24+Q/J0CX4yi9A+39C3Sfgvnr4m7fCfzV4OlQ7Mi7ngJKc1cnZMT4RCoWora31ehiB58uchzqgpBJC/vnfOSNz3vtu9G/U87Xcj7H3C5chvwxe+yG46yMvZfu7nfC5V8GdvwYF670dbxIyMucAkpzVyUyIT1iWxdjYmEz3OcyXOfvonjEJGZmzptF71yeo0Aa59vr/BR95Dt7wyMLmrqAM9v4cNP29d+O0ISNzDiDJWZ00IT5hmiZdXV0y3ecwX+bss4NSIXNzHq56J3dO/SWR298HuflLv+jAL0HzKXcHpihTcw4ayVmdNCE+EQqF2Lt3rxx97TBf5nzjImzwz+m54NOcEza9Akb7YfSG1yNZla9z9hHJWZ00IT5hWRbDw8My3ecwX+bsw5kQX+acoOtQeTd0fdfrkazK1zn7iOSsTpoQnzBNk2vXrsl0n8N8l7MxC5EuKK32eiS2+C7nl9t5D1x6yutRrMr3OfuE5KxOmhCfCIVC7N69W6b7HOa7nKOXoWgLrHHmTrdO8V3OL7fzXl80Ib7P2SckZ3X+Oacvy1mWRSwWo7i4GC2Ad/vMFL7L2WcXKUvwMuf2/lGl5xYoq4aZCYh1x+9cnKF8V88+JTmrkybEJ0zT5MaNGxQWFkq37SDf5ezD03PBm5xLCnLJywnx8OPnVnxdXk6IkoJVZpY07eYume/C/vemb5Bp5rt69inJWZ00IT4RCoXYtWuX18MIPN/lPPAiVP2U16OwzYucK8J5ND5yL5Gx6RVfV1KQS0U4b/UPTBwXksFNiO/q2ackZ3XShPiEaZpEIhFKSkrQdTmUxym+y/lGK7zmQ16Pwjavcq4I5yXXYCRj5xvgO/8HLCs+M5KBfFfPPiU5q5O0fMKyLKLRqJwC5jBf5WxZ8dvPr/ffb2C+ynk54e3x72Ckz+uRLCsQOfuA5KxOZkJ8IhQKUV3tr9Mw/chXOce6Ib8Ecgu8Holtvsp5JZtq4fp5KNrs9UiWFJicM5zkrE5mQnzCNE36+/vlPHSH+SrngYuw3n8HpYLPcl7Jpn3Q97zXo1hWYHLOcJKzOmlCfGR8fNzrIWQF3+Ts09NzE3yT80o27oO+816PYkWByNkHJGc10oT4hK7rVFZWykFPDvNVzj5uQnyV80o2vSK+OyZDBSbnDCc5q5PEfMI0Tfr6+mS6z2G+ytmH94xJ8FXOKyndCcO9MJ2ZvwUHJucMJzmrkybER6anV76+gUgPX+RsWTePCfHX3XPn80XOq9FDUL4b+l/weiTLCkTOPiA5q5EmxCd0XWf79u0y3ecw3+Q8dgNCayEv7PVIlPgm52Rs3Juxu2QClXMGk5zVSWI+YZomPT09Mt3nMN/k7OPjQcBHOSdj/W3xK9dmoEDlnMEkZ3XShAjhRz69Z0wgbbg1Y5sQITKdXKzMJ3Rdp6KiwuthBJ5vco5egZJKr0ehzDc5J2N95jYhgco5g0nO6mQmxCdM0+TKlSsy3ecw3+Q83ANFW7wehTLf5JyMoq0wNgAzE16PZJFA5ZzBJGd10oT4SG7uKrcXF2nhi5yHr8X/8fMxX+ScDF2H0ioY7PB6JEsKTM4ZTnJWI02IT+i6zqZNm+Toa4f5JueYv2dCfJNzsjbcFj9lOsMELucMJTmrk8R8wjRNurq6ZLrPYb7IeXocJiJQmJk3TUuGL3K2Y/2tMNDm9SgWCVzOGUpyVidNiI/k5+d7PYSskPE5D1yE9bviuwF8LONztmP9rfHTpjNQoHLOYJKzGttnx3R2dtLQ0DB32+JwOMyRI0dWfV9jYyNnzpwhGo3S2dnJ/fffn9T7RJyu65SXl3s9jMDzRc79rVC+x+tRpMQXOdtRVgNDmXdMSOByzlCSszpbTUhnZycHDx7k0qVLhMNhAOrr6zl27BhHjx5d9n2NjY00NTXR0NAAQDQa5eDBgzzzzDMcP35cffRZxDAMurq6qKysJBQKeT2cwPJFzv0tUH6716NIiS9ytqN0Jwxdil9OP4MELucMJTmrszWf29DQwJEjR+YaEIBHH32U+vr6Fd93/PjxBU1KOBymvr6eEydO0NnZaW/EWUrTNMLhMJqmeT2UQPNFzv0v+L4J8UXOduQWQE4+jA96PZIFApdzhpKc1dlqQk6ePDm3GyYh0ZA0NjYu+77Tp08valRe9apXrfo+8RJd1ykrK5Ojrx3mi5yvXwjE7piMz9mu0ioYyqxfqgKZcwaSnNUlnVg0GiUajVJVVbXouXA4TFNT07LvPXz48KLmxc56h4eH59Y/NTWFaZpzRyHPXzYMY8GydXNq1M6yZVmLlgHby6ZpLlhearzJLCe2yTAMLl68yOzsbGC2KRO/p5mZGdra2piZmcnMbRobxDKmoHDTituUkKnf0+zsLBcvXpx7LhC1V7oTc6B9LvtM2KbZ2Vna2tqYnp6WnxEObtPMzAwvvvji3LiDsE2pfk/JSroJWWm3SWlpKYODy09Dnjp1atFBqGfPngXg0KFDK653x44dFBcXU1JSQklJCY899hjd3d309/cDcOXKFQYGBgDo6uoiEonMjTcWiwHQ1tbGyMgIAK2trYyPjwPQ0tLC5OQkAM3NzXP/8DQ3N8/9g9Tc3AzA5OQkLS0tAIyPj9Pa2grAyMgIbW3xU/NisdhcTpFIhK6uLgAGBga4cuUKAP39/XR3dwPQ29tLb28vwKrbpOs6U1NTc9sRhG3KxO+pu7ubDRs2zC1n2jZNXj7L2LqdoGnLbtPwzfVl8vc0MzPDyMgIuq4HpvZmi7Zz48Ufz2WfGK+X2zQ2NsaGDRtoa2uTnxEOblN/fz+apqHremC2KZXvqaMj+YO0NSvRuqyiqamJgwcPcubMmUWNQ3V1NYcOHbJ1kGl1dTUPPfTQsge0Dg8PU1xczOXLlwmHw3PdVl5eHjk5OUB8CizxuK7rGIYxVwiGYaDrOpqm2VqGeKc3fzkUCmFZlq1l0zSxLGtueanxJrMs2yTbtGDs3/8c1nAP+lsfW3abnr8a4V1//n2e+M272bO5MPO3KSjf0/l/xmr9Bi2v+xPe+dmn+fqHX8cd20v9vU1B/J5kmxzfpkgkQmlpKbFYjKKiIlaS9Nkx8w9GfbmhoaFkPwaA+++/n0OHDq14Rs389a60EfP3wc0/Kjndy5qm2VqeP65UlhOfZxgGbW1t7Nq1i1AoFIhtcmI51W0yDIPW1ta5nDNum65fQNt5j61tUl12cptM01yynv28TVppFVrk0uLHPdym+T83EgdNys+I9G+TZVkL6jkI25Su72k1Se+OKS0tBeLHaLxcNBpdsUmZ78SJE5SWlsqpuTbpus6WLVsWFIBIv4zP+XozbNrn9ShSlvE5qyjdmZEHpgYu5wwkOatLOrFwOEw4HF521qOurm7Vzzh9+jTRaHRBA7JUUyMW0zSNoqIiOQXMYRmdszETvxbF+tu8HknKMjpnVbcUgxYiNBn1eiRzAplzBpKc1dlq2x544IFFB5wkDnRZ7QDTpqYmhoaGFuyCiUajcopukgzD4MKFC7aOOhb2ZXTOAy9CSSWs8f/dOjM651SUVpE7ctnrUcwJbM4ZRnJWZ6sJqa+v5/Tp0wseO378+KKZjbq6ugWn7HZ2dvLYY49RWlrK6dOn5/7U19cvecqvWEzXdSorK2W6z2EZnXPfedjo/10xkOE5p6J0J7mxLq9HMSewOWcYyVmdrcu2V1VVcerUKerr67nzzjvp7OykrKxswem3Q0NDnD17dsFum4MHDxKNRhc1MIAcG5IkTdMoKCjwehiBl9E5B+R4EMjwnFNRWkXucBdQ5vVIgADnnGEkZ3W2b2B34MABDhw4sOzzVVVVc+cNJ7z878I+wzBoaWlhz549to48FvZkdM595+HulXd7+kVG55yK0ipye/4dOOj1SIAA55xhJGd1MnfkE7quU1NTI9N9DsvonK+fh421Xo8iLTI651TMzYRkhsDmnGEkZ3W2Z0KENzRNIy8vz+thBF6qOfdEJ4iMTS/7fElBLhVhhc8fuQ76GijIjGn+VAW2nkurWDucOQemBjbnDCM5q5MmxCcMw6C5uZna2lqZ7nNQKjn3RCc49KknmZhZ/gj5vJwQjY/ca78Rud4cmINSIcD1nFeCZkyzjnGvRwIEOOcMIzmrkybEJ3RdZ8+ePTLd57BUco6MTTMxY/CZB/dTU75u0fPt/aM8/Pg5ImPT9puQvvOBOSgVAlzPmsZUUSU7Rvq9HgkQ4JwzjOSsTpoQH5EO2x2p5lxTvo59FcVpGs1N18/DrW9N72d6LKj1PF1cyY5rfV4PY05Qc840krMaadt8Yv7dE4VzMjbnvvOw6aWDUnuiE5zviS35p71/1MOBJidjc06D6aJKKrXrXg8DCHbOmURyViczIT6h6zq1tbUy3eewjMx5uBfG+qG0Gkj+2JOSgsy9smpG5pwm00U72KE97/UwgGDnnEkkZ3XShPjI/Ns5C+dkXM7nT8Ped0Mo/r/raseeQApn4bgo43JOk+miSir1zNkdE9ScM43krEYS8wnTNGlpaZHpPodlZM6dT8KuxTeITBx7stSfTG9AMjLnNJkqqmRHBu2OCWrOmURyVidNiE+EQiH2798vBz85LONynhqB7p/Ajru8HklaZVzOaWTkraeASbTZCa+HEuicM4nkrE6aEJ+wLIuJiQksy/J6KIHmdM63alcpf+ZP4N9+D3qTOG7g/Feh6l5YW+jIeLwS6HrWNK5Y5eQOX/F6JMHOOYNIzuqkCfEJ0zRpb2+X6T6HOZbz+BDhiyf569xPok8Pw9p18Pj7YHqFi1oNX4P/aoC7fyu9Y8kAQa/nLmtjRly+Peg5ZwrJWZ0cmOoToVCI2tpg3DckkzmW84/+ik0/+SK/PvNLfOR1D7O+ohgmIvC1D8EDf7/gpYlLv1f81/9k+rZf5IZVBT2xuef9cAruaoJez5etTRzIgCYk6DlnCslZnTQhPmFZFuPj4+Tn56NpmtfDCSxHcp6dgme/ROc7T3LmH/r4SOLxt/wRfP6N8LUPw8/8OWja3Om35swEP1r779RNfZrhp59e9JGZfgruaoJez13WRtbGurweRuBzzhSSszppQnzCNE26urrYvXu3HPzkIEdyfuFfYeNepsM1wLxTN9eshQ/+J/ztO+HpP4E3PDJ3+u2/vOYyWuwNfLlu6auk+uEU3JUEvZ4vm5vIHf5Pr4cR+JwzheSsTpoQnwiFQuzdu9frYQSeIzk/+yW489eWfm7NWnjPP8LfvB3WFpGXcyv7tE72vfg5cn75XyjelObLv2eIoNdz/JgQ7++mG/ScM4XkrE6aEJ+wLIuRkREKCwtlus9Bac851h2/78uut8D1ZU7ZXFcOP/95OHEf1cATa+Hy3Z9nx6bg7mMOej33EyY0ORTfFbdmrWfjCHrOmUJyVidnx/iEaZpcu3ZNjr52WNpzfu4rsO8wrFnl+I0tr4T/L8L5D1xm9+TfMFL5lvSsP0MFvZ4tdGYKt0LE29mQoOecKSRnddKE+EQoFJL9jS5Ia86WBef+Efb/QnKv13XQNCbx7jdnt2RDPU8XVcJQp6djyIacM4HkrE6aEJ+wLItoNCoXw3FYWnO++iPILYDNr0j9swImG+o5E5qQbMg5E0jO6uSYEJ8wTZMbN25QWFgo3baD0przha9B7f1pGVfQZEM9d2ub0K620rs9tuTzbpzhlA05ZwLJWZ00IT4RCoXYtWuX18MIvLTmfOkpePfx9HxWwAS5nksKcsnLCfHJszMcCT3LLzctvs4LxK/10vjIvY42IkHOOZNIzuqkCfEJ0zSJRCKUlJTI7aIdlLacR2/A6HUol9P2lhLkeq4I59H4yL2M9m2l8ptf5okH7170mvb+UR5+/Fz8yrgONiFBzjmTSM7qpAnxicQ+x3A47PVQAi1tOXc9BZV3xw82FYsEvZ4rwnlQdDt89Qb7NuVDKMeTcQQ950whOauTJsQnQqEQ1dXVXg8j8NKWc+eTsPOeJZ9a6d4vQbgvTDKyop71EBRVQPQKlHmzrVmRcwaQnNVJE+ITpmkyMDDA+vXrZbrPQWnL+dJT8PqPLHgocazAw4+fW/Gtfr8vTDKypp7LamCww7MmJGty9pjkrE6aEB8ZH1/htu8ibVLNOWekG4xpKK1a8HjiWIHI2PSK7/f7fWGSlRX1XFYNg+3Amz0bQlbknAEkZzXShPiErutUVlZ6PYzAS0fOBde+DzvvhSUu31wRzsuKBmM1WVPPZTXQ1+zZ6rMmZ49Jzupk3sgnTNOkr69PLgvssHTkvO7a96Dq3jSOKniypp7Lam7OhHgja3L2mOSsTpoQH5meXnkaX6RHqjnn9/0YdtyVptEEV1bUc1mN51dNzYqcM4DkrEZ2x/iEruts377d62EEXqo5b2YQzTIhLN/VSrKmngs3wWQMpkZh7TrXV581OXtMclYnMyE+YZomPT09Mt3nsFRzfpV+kbGNr0rzqIIna+pZ02DTK6D3nCerz5qcPSY5q5MmRIg0Oqi/yPimO70ehsgkWw9C91mvRyFERpImxCd0XaeiokLOQXdYqjm/Sn+RcZkJWVVW1XP5Hrhx0ZNVZ1XOHpKc1UliPmGaJleuXJHpPoelkrM+Pcp2rZ/J0t0OjCxYsqqe1++CwTZPVp1VOXtIclYnTYiP5OYG+yqamUI157wb53je3Am6HO+djKyp57JdMPAiWJYnq8+anD0mOauRJsQndF1n06ZNMt3nsFRyzr/+DM9YtzowquDJqnq+pQjW5MHYgOurzqqcPSQ5q5PEfMI0Tbq6umS6z2Gp5Jzf38Sz5i4HRhU8WVfP62/Ohrgs63L2iOSsTpoQH8nPz/d6CFlBKWfLIq//WZ41a9I/oIDKqnpefyvceMGTVWdVzh6SnNVIE+ITuq5TXl4u030OU8451o2ZU0gM9y9I5UdZV88VB6GnyfXVZl3OHpGc1UliPmEYBh0dHRiG4fVQAk055xutTJbIrphkZV09b70Tun/i+mqzLmePSM7qpAnxCU3TCIfDaEvcmVWkj3LO/S1MlchBqcnKunouq4HRfpiIuLrarMvZI5KzOmlCfELXdcrKymS6z2HKOfedZ7JsrzODCqCsq2ddv7lL5hmXV5tlOXtEclYnifmEYRi0tbXJdJ/DlHPue57Jsj3ODCqAsrKet97p+uXbszJnD0jO6qQJ8Qld19mwYYN02g5TynlmAmI9TBVXOTewgMnKet7m/nEhWZmzByRndXJpR59I7HMUzlLK+XoLbLgN9JAjYwqirKznxO4YF68lkZU5e0ByVidtm08YhkFra6tM9zlMKee+52FTrXODCqCsrOe8EijY4Op9ZLIyZw9IzuqkCfEJXdfZsmWLTPc5TCnnvudh8yucG1QAZW09V74BLj3l2uqyNmeXSc7qJDGf0DSNoqIiOQXMYUo59zXDJmlC7Mjaeq66Fy496drqsjZnl0nO6qQJ8QnDMLhw4YJM9znMds6mAf2tUC5nxtiRtfW8eT/0Pufa6rI2Z5dJzuqkCfEJXdeprKyU6T6H2c55sAOKtkCu3DfCjqyt5/B2mBohNOnORcuyNmeXSc7qJDGf0DSNgoICme5zmO2c5XgQJVlbz5oGlXdT0PNdl1aXpTm7THJWJ02ITxiGQXNzs0z3Ocx2znJmjJKsrueaOgqv/pcrq8rqnF0kOauTJsQndF2npqZGpvscZjvnvmZpQhRkdT3vqmNd95NoOH+9kKzO2UWSszpJzCc0TSMvL0+m+xxmK2fLkjNjFGV1PRdtYTZvPXu1LsdXldU5u0hyVidNiE8YhsG5c+dkus9htnIe6QM9BwrWOz+wgMn2eh7ddh/36c6fJZPtObtFclYnTYhP6LrOnj17ZLrPYbZy7muWg1IVZXs9j2x7Iz8VOuf4erI9Z7dIzuokMR8JheTeJG5IOue+5+R4kBRkcz2PbzxAldZLaDLq+LqyOWc3Sc5qpAnxCdM0aW5uxnTx5lfZyFbOvc/L8SCKsr6e9Ry+b+5hXfd3HF1N1ufsEslZnTQhPqHrOrW1tTLd5zBbOV87B1te6fiYgkjqGf7deDXFnd90dB2SszskZ3WSmI/IQU/uSCrnsQGYGYfirc4PKKCyvZ6/bR4gv++HMDns6HqyPWe3SM5qbDchnZ2dPPTQQxw7doxjx45x4sQJW+8/ceIE9fX1dleb9UzTpKWlRab7HJZ0zolZEDklT4nUM4xzC6Nb74NzX3ZsHZKzOyRndWvsvLizs5ODBw9y6dIlwuEwAPX19Rw7doyjR4+u+L6GhgYATp48yZEjR9RHnKVCoRD79+/3ehiBl3TO156VXTEpkHqOu7H/Nwh/65fh4K9Azi1p/3zJ2R2SszpbMyENDQ0cOXJkrgEBePTRR1ed2aiqquL48eMcP36cqqoqpYFmO8uymJiYwLIsr4cSaEnnLE1ISqSe46ZKb4Ptr4Fn/taRz5ec3SE5q7PVhJw8eZLq6uoFjyUaksbGxrQNSixmmibt7e0y3eewpHPuPQdb9rsxpECSep7n3nr43p/CzETaP1pydofkrC7pJiQajRKNRpecyQiHwzQ1NaV1YPPXOzw8PLf+qakpTNOc+7LnLxuGsWA50ZXaWbYsa9EyYHvZNM0Fy0uNN5nlxDaFQqEFF8MJwjZl4vekaRq1tbVomrb8doz0gzGNWbDRF9uUid9T4uJOoVAoMNtkdznxmZTfjrX9NZjP/F3atylx1kZi3U5vUxC/p2S2SdM09u7dSygUCsw2pfo9JSvpJqSzs3PZ50pLSxkcHEx6pXbs2LGD4uJiSkpKKCkp4bHHHqO7u5v+/n4Arly5wsDAAABdXV1EIpG58cZiMQDa2toYGRkBoLW1lfHxcQBaWlqYnJwEoLm5mZmZmQXne8/MzNDc3AzA5OQkLS0tAIyPj9Pa2grAyMgIbW1tAMRisbmcIpEIXV1dAAwMDHDlyhUA+vv76e7uBqC3t5fe3l6AVbfJsixaW1uJRqOB2aZM/J4uX77M2NgY169fX7RNPdEJvv1sG80/amQkfDuNTS/y9IXLnO+J8a2zrTzXFd/Wnp6ejNqmTP2ennvuOSzLCtQ22fmeAK7dHG/k1vcw+8PPw+xUWrdpeHiYsbExqT0XtqmrqwvLsgK1TarfU0dHB8nSrCR3YjU1NXHw4EHOnDnDoUOHFjxXXV3NoUOHOH78+Kqfc/DgQQ4dOjR3oOpyhoeHKS4u5vLly4TD4bluKy8vj5ycHCD+21TicV3XMQwDTdPmlnVdR9M0W8sQ7/TmLyd+W7OzbJomlmXNLS813mSWE9tkWRYvvPACu3fvZs2aNYHYpkz8ngzD4MUXX+TWW28lFArNjb0nOsGbP/1dJmYMPhh6gjJthD+efe+ius3LCfGth+9ma2lBxmxTJn5Ps7OzvPDCC3Oze0HYJjvf0wt9o7zzs0/z9Q+/jju2l2IaBtrJ96Nt3IN530fTtk2WZXHx4kV27dpFTk6O1J5D2zQzM8PFixe5/fbb525i5/dtSuV7ikQilJaWEovFKCoqYiVJnx0z/2DUlxsaGkr2Y2wLh8MrbsT8i8PMv2xuupc1TbO1PH9cqSzPH8u+ffsCt02Z9j3pus7evXuZT9d1YhOzTMwYfObB/dx34RQj297M3TV383IlBblUhPMyaptUl538ntasWTO3myCZbfXDNtn9nua/Vw+F4K1/BF+oQ7/zg1C4MW3b9PJ6zvbac2KbcnJyFvx8DsI2pet7Wk3Su2NKS0sB5nYHzBeNRldsUkTqLMtieHh4bp+bcMZqOdeUryMce4Ftt7+afRXFi/7Mb0DE8qSel1BSCa/+IHzjtyFNuUjO7pCc1SXdhITDYcLh8LKzHnV1dWkblFjMNE2uXbs2NxUmnLFazvpUDEZ6Yf0ul0cWLFLPy7j7tyB6Bc5/NS0fJzm7Q3JWZ+sU3QceeGDRASeJA11efpyISK9QKMTu3bttTXMJ+1bLOb//HFQcBF2+h1RIPS8jlAPv+jP49ifAmE394yRnV0jO6mw1IfX19Zw+fXrBY4mLkCVEo1Hq6uqWPWU3caqtsMeyLKLRqEz3OWy1nG8ZbJaLlKWB1PMKtrwSSquh5Wspf5Tk7A7JWZ2ty7ZXVVVx6tQp6uvrufPOO+ns7KSsrGzBZdiHhoY4e/bsgt020WiUxx57jGg0SmdnJydPngTiZ9WsdLl38RLTNLlx4waFhYXSbTtotZzzBi7Aqx7wYGTBIvW8ivt+D776Qbjt7ZCbr/wxkrM7JGd1SZ+i67bEKbrJnOIjhNPO98R452ef5mL577P2l78KZdVeD0n4WKKePvPgfmrK1y35mq3f+Qjapr0UH/pdl0cnRGrs/PttayZEeMc0TSKRCCUlJQtOkRLptVLOBUywZnIQSnZ6NLrgyPZ6LinIJS8nxMOPn1v2Ndu1e/jnto/Tv+1uym97ndJ6sj1nt0jO6qQJ8YnEPkc5FdpZK+W8W7vCZMltFMgPmZRlez1XhPNofOReImPTy76mvX+U3z11jb/6xoeg5mz8oFWbsj1nt0jO6qQJ8YlQKLTo5oEi/VbKea/exWTZXgpcHlMQST3HG5HVrivzHfOVTJY2s/a7fwL3rXy38qVIzu6QnNXJr3Q+YZom/f39ch66w1bKea92mcmyPR6MKniknpP3vd2/z+wPj/PihWc53xNb8KcnuvKddyVnd0jO6mQmxEcSNzYSzloyZ8viLv0CY5s/4f6AAkrqeWWJ40Y+/LWr/Jz+AL/x+Ht57/T/YoSXzpbJywnR+Mi9K86oSM7ukJzVSBPiE7quU1lZ6fUwAm+5nHNHLgMwXbT4OWGf1PPqFh43cjfh70/w1MwZeu79JBA/ZuThx88RGZtetgmRnN0hOauT3TE+YZomfX19Mt3nsOVyzrvexE+s2+DmHTJFaqSek1MRzpu7L1HZzzxGyeAz7Jv4Cfsqipc9tXc+ydkdkrM6mQnxkenp5Y+kF+mzVM75/c/yrFnDrR6MJ6iknm3KzYdDH4enPwPVb5p7uL1/dNm3mKbJ1PAo5eXlzo8vy0k9q5EmxCd0XWf79u1eDyPwlss5v/8Zmsz38qAHYwoiqWdFt74NnvoktHyNkoq3rnqtEUgcN1Ihd3h2kNSzOmlCfMI0TXp7e9m8ebNcDMdBS+Y8PU5u7DIXrW3eDi5ApJ4VhdbAm/83fOMRKh5626rXGnnx+jC/ffJ5BkcmpQlxkNSzOmlChFhN7zkmy25ndkT+dxEZYOc9UL4Hzn6Ritf99xWbCzlGQWQ6adl8Qtd1KioqpMt22JI5X/kh4+UHvBtUAEk9p+inPgrf/yzMrHydkES+krOzpJ7VSWI+YZomV65ckd9sHLZkzh3/yejWN3g3qACSek7Rhttg++vgmb9b8WWJfCVnZ0k9q5MmxEdyc3O9HkJWWJDz5DD0Pc/4pld7N6CAknpO0T2/C9/7U5iZ9HokAqlnVdKE+ISu62zatEmm+xy2KOfO/4Ltd2GF1no6rqCRek6DjXtg66vg2X9Y9iWyO8YdUs/qJDGfME2Trq4ume5z2KKc28/ArjpvBxVAUs9pcu/R+HVDZqeWfFp2x7hD6lmdNCE+kp+fv/qLRMrmcrYsaGuUJsQhUs9psKkWtuyHs3/j9UiyntSzGmlCfELXdcrLy2W6z2ELcr5+AdYWQlguQpRuUs9pdOgT8N3/C1Mji56S3THukHpWJ4n5hGEYdHR0YBiG10MJtAU5y64Yx0g9p9H6Gqipg598cdFTiXwlZ2dJPauTJsQnNE0jHA6jyQ3UHLUg59Zvwq43ez2kQJJ6TrM3PAI/+HMYH1rwcCJfydlZUs/qpAnxCV3XKSsrk+k+h83lPNgOw9eg8m6vhxRIUs9ptr4Gdr8Dnll4bIjsjnGH1LM6ScwnDMOgra1NpvsclsjZbPp7eOUvgh7yekiBJPXsgNf+enyXjDEz95DsjnGH1LM6aUJ8Qtd1NmzYIJ22w3RdZ0NpGK35JOz/Ra+HE1hSzw7YcBts3AfP/dPcQzIT4g6pZ3WSmE/IPkd3aJpGuP+HaOW3Q8kOr4cTWFLPDvmpj8KTn5y7iqocE+IOqWd10oT4hGEYtLa2ynSfwwzDYPS7f4m5/31eDyXQpJ4dsmV//CqqPz4OyO4Yt0g9q5MmxCd0XWfLli0y3ecwPXaF/Fgb2u0/7fVQAk3q2UGH/iB+h93xIdkd4xKpZ3WSmE9omkZRUZFM9zlM+/5n0V99BC3nFq+HEmhSzw4qqYRXPAhPHpPdMS6RelYnTYhPGIbBhQsXZLrPSS/8K1bbf9BadLfk7DCpZ4fd8ztw/quEhtoA2R3jNKlnddKE+ISu61RWVsp0n1Mmh+EbvwPv+TLbbn2F5OwwqWeH5ZXAffVs+8HHAEtydpjUszpJzCc0TaOgoECm+5xgGvD4L8IdD6JtvkNydoHUswsO/jc0Y4qf178rOTtM6lmdNCE+YRgGzc3NMt3nhB98DjQd3vRxydklkrMLdJ3uu/4P9Tlf4fLlS5zviS37pyc64fVofU3qWd0arwcgkqPrOjU1NTLdl06WFb+w0w/+Ao58B3Qd3bIkZxdIPbujYPsdnLTeROW/HeWdMx9Z9nV5OSEaH7mXinCei6MLDqlnddKE+ISmaeTlyQ+ItPrx56Hp7+D9/wxFWwDJ2S2Sszu2luTzcx/5NKX/+BaeemWU4ap3LnpNe/8oDz9+jsjYtDQhiqSe1Unb5hOGYXDu3DmZ7ksHy4IfHYenPw3v+TJs3Dv3lOTsDsnZHYZhcOPqJXLf/Rds/9Efsq9gmH0VxQv+1JSv83qYvif1rE5mQnxC13X27Nkj032pmp2GJ34L+i/AB85A8VYAeqITRMamsSwLK7yVlt6RBQeZtfePejXiQJJ6dsdczjk58IZH4OT74Vf/HeQ6OGkl9axOmhAfCYXkjq4pMWbiZ8GsLVzwg7gnOsGhTz3JxMzKv8Xk5YQoKch1Y6RZQerZHXM5v/qDcK0JvvkIvOtzIGdypJXUsxppQnzCNE2am5upra2VYldx6Slo/DgUb4N3fx70lzKMjE0zMWPwmQf3s7Msj/b2dmpqahblXFKQK/vM00Tq2R2Lcn7np+Gv3xo/Huo1R7weXmBIPauTJsQndF2ntrZWpvvsMmbhzMeg5evw9k/CbW9f9jfAmvJ17N1SxL6Kg+i6Luf8O0jq2R2Lcs7Jg/f8I3z+TbDjdbCp1tsBBoTUszpJzEfkoCebZqfg1C9DpAs+9DTsfkdSU9CSszskZ3csyrl4K7z1MTj1KzAlxzqli9SzGmlCfMI0TVpaWjBN0+uh+MNoP3zhEOQWwAN/D/mlSb1NcnaH5OyOZXPe926ofiOc+hW02UlvBhcgUs/qpAnxiVAoxP79+2V/42osC154Ar5YB3t/Fn7uOIRykn675OwOydkdK+b81j+Ggg3s+NavcQtT7g8uQKSe1UkT4hOWZTExMYFlWV4PJXNZFjzZAI1/AO/4FNz927bPAJCc3SE5u2PFnPUQ/MyfM124lb/LbUCfiro+vqCQelYnTYhPmKZJe3u7TPctoSc6wfnuKDe+/vtMnH+C1rd+hfN5d3L+2rDt+2JIzu6QnN2xas66zrW7/5iz5q1U/ethiHW7O8CAkHpWJ2fH+EQoFKK2Vo5kf7nENT4+YJ7iLaGzvGn69xn+4ovAi3OvsXNfDMnZHZKzO5LKWdP45Ox7+Lk9nWz54lvgFx6HTfvcGWBASD2rkybEJyzLYnx8nPz8fDl1dJ7I2DTvM7/OB0qf49rPfo0v37LwAFS798WQnN0hObvDTs5De36JLduq4B9+Fu75XbjzgyCnnCZF6lmdVJhPmKZJV1eXTPe9TOmFv+U9oe9w7V1f4fbqnSnfF0Nydofk7A7bOe9+B3ygES58DU7/Kswkvyszm0k9q5OZEJ8IhULs3bt39Rdmk/NfZX3zF7hn+qN8IX/Dii9d6d4v85+TnN0hObtDKeeSSvilr8M3fgv+5u3w3n+Cwk2OjC8opJ7VSRPiE5ZlMTIyQmFhoUz3AVx7Fv7t97j8tn+i70vXl31ZSUEueTkhHn783Iofl7gvjOTsDsnZHco5r8mN31/mB5+LX2/nwS/Blv2OjdPvpJ7VSRPiE6Zpcu3aNXbt2iXnoo9ch8ffDz/zOaYKdgHLNyEV4TwaH7mXyNj0ih+ZuC+MYRiSswuknt2RUs6aBnf9Jqy/Ff7x/vhtD/b+rCPj9DupZ3XShPhEKBRi9+7dXg/De7NT8Pj74M4PwK1vgZ7Yqm+pCOclfeM5ydkdkrM77OS87C7Lgtey4d0n2fjEr0DPM3Df78WvRCzmSD2rkybEJyzLIhaLUVxcnL3TfZYFT/xWfJ/16z/i0CokZzdIzu5IJudkdlnm5YT4zw9/k80//iP489fA2xriB7EKQOo5FdKE+IRpmty4cYPCwsLsne774V9C/wvwq9+0fSXUZEnO7pCc3ZFMzqvtskyc5j5ormPzz/w5XP0x/L/fhPP/DG87BgVlTm6CL0g9q5MmxCdCoRC7du3yehjeaf92/CC5DzTGb0fukKzP2SWSszuSzdnOLku2vRoeegqePAZ/9fr4Aay7DqU4Un+TelYnTYhPmKZJJBKhpKQEPYsuINQTnWC0t52q//cQV978BSZGCmDkpeNAVjr1VkW25uw2ydkdjuW8Zi286WNw29vi1xOpfQDuexRCK/+T0hOdSPogcT+RelYnTYhPWJZFNBolHA57PRTX9EQn+OlP/Qd/p/0Bf2C8iy8/PgY8veh1idNr0yEbc/aC5OwOx3Pe+io48iR8/b/DX90Nhz4eP2B8id2liVssTMwYK36kndssZAqpZ3XShPhEKBSiurra62G4KjIyQQN/xrqqV/MLh/4/fmGFA+vS9QMrG3P2guTsDldyzi+NX9Cs4z+h8ePww7+An/0LKN664GWRsWkmZgw+8+D+Za9kbPc2C5lC6lmdNCE+YZomAwMDrF+/Pmum+zb98A/px2DsTX/Mvq1hV9aZjTl7QXJ2h6s5V78Rdt4Xb0K+cAje+sew52cWzYrUlK9jX0Wxs2NxmdSzOmlCfGR8fNzrIbjn+5+loO8n/MbMb3NSd7dMsypnD0nO7khXzqsdfzU3I3nXb8D218K/HYXvfxbq/hAqX5+WMWQyqWc10oT4hK7rVFZWej0M51kWPPV/ofkUl9/6Jca/2Obq6rMmZ49Jzu5IR852bn0wdyzH1lfBB74NL/y/+Om8ZdXk3/YBwEppLJlK6lmdNCE+YZom/f39lJeXB3e6z7LgW/8Tup6GX/0ms9EcwN0mJCtyzgCSszvSkXMytz5Y8lgOTYvvjrnt7fDcV9j89Mf5Zu4Ypc/9AoTeCRv3Ona9H7dJPauz3YR0dnbS0NAwdxBOOBzmyJEjjr1PvGR6euVT23xtegy+8TsQvQy//K9wSxFEV78kuyNDCXLOGURydkc6crZ1HZGXC+XAgffTUf7TfOLPv8hfTvTBVz8Ak1GofhPUvAmq7osf4OpjUs9qbDUhnZ2dHDx4kEuXLs2dilRfX8+xY8c4evRo2t8nXjIzM8Nf//Vf8+ijj7J27Vqvh5NeXd+Ln+K38x74xdOQm+/ZUAKdcwaRnN3hds4rHTfSfmOMn1i76XvtB1hfUQyx7vhFCFu+Bt/4bSiroXzDXbxdt1g7tB7K74CcW5TG4fb1SKSe1WmWZSW9k+6hhx4iHA7T0NAw91g0GqWkpISVPkblfcPDwxQXFxOLxSgqKkp2iIGVyCsSiQTnXPTBDvjP/w3dZ+GnPw01C6+6eL4nxjs/+zRP/Obdrh1NH8icM5Dk7A63ck75GiDGLPSc5fpz3+JHP/4+deXD5I1choL18bv4rr8V1u+K/7esBvLLQF/68uheXI9E6nkhO/9+25oJOXny5IJGApgLvLGxkUOHlr50r+r7RADNTkPXU/DsP8Ll78Fd/wN+9i+Vf+MRQngvmeNGYIXZh9Aa2P5aboRu53987yBPvPtu9m0uhOFuGHgRBtrg2jl4/hQMdcD4UPxOvgXrYd1GWFce/29+GdrEGt5u9HL4vjvYsGkbs3kbMPLKsEIvXdDQr9cjCaKkm5BoNEo0GqWqqmrRc+FwmKampiWbCdX3iYUSBzv56qCn6TGI9cBQJ/Q9D73PweXvQ/keeMUD8K7PwtqlL1rkFV/m7EOSszvczDml40aWousQ3h7/87JZUiwLpkZgfABG+2H0evy/YwPkjl7nVXobe3s7KeqOxJ8bGwA9J36zvfz1bA8V88drZtj446dhY0W8mclfH38+rwTWFsHawvjxLEkNVepZVdJNSGdn57LPlZaWMjg4mNb3JXbTXL16leLiYkzTBCAvL4+cnHhh6Lo+97iu6xiGgaZpc8u6rqNp2uLlH3wObXwA07LQsNA07ebyS+ueW8ZCJ35imWVZ6Fq8/i0sdE3Dsqybj99chpdeM+9xiB8IPrc8bxsXPq5hYYGVGJeJBsxMz/BX78hj6qu/znDu2pfGzs2xa/OWb35+/PF52xRf0Uvj4mVjX+1xy4KbeS3YDoDZabSZMazpUZgeQ5uIYJmzULQFrXgH1sY9WFXvQr/3f2PcEo5/T1Mmxnhk2e9pZHgEc2qcWDRCbB0LXgPxI9LnL4dCISzLsrVsmiaWZc0tR6NRACKRyNznz6+x5ZaTrr0Vlp3aptXG7sU2RSIR1qxZQzQaDcw2ZeL3NDw8DMDQ0NDcz4BM3qbhWAxzapxz7T0Mx6JJfE9r0LQKDGsTeqGOVqTRro/w0bELVL72NezeVPjSNk2OEJqMYk0M0XXlCt89f5YDRh4516+iTz6HNTYA44NokyNY0yMwNYpmGVi562DtOrTcdVhr10FOHpqWg6WHQA+h6SGmZgz+4h23MPW1jxBbewtoOpqmx3+W3/y5nvh5/NKyNu/n9M3lmz9T5/+btOjfp5uvSSzrN9djWcz7+b3Csq7f/LfKQtd0zPs+CqHctNZeJBKJjy+Joz3Sdopu4od3ut43MjICwL59+xRHFEwf+saXvR6CTb3AM8A/K3/C3Z9J11iSJ+f8u2PHjh1eDyEr7Ny50+sh2PL+z6T+GXcl8Rl/kdQnDSS9zg9/42+Tfm3m+FPHPnlkZITi4pWP50u6CVnpYJuhoaG0v2/Lli10dHSQk5MT/w3+prVr18rRx0IIIUSGsiyLkZERtmzZsuprk25CSkvj53AvNXOx0t0DVd+n6/qSx5EIIYQQIrOtNgOSkPRRNOFwmHA4vOzsRV1dXVrfJ4QQQohgs3Uo7wMPPEBHR8eCxxIHnq50hovq+4QQQggRXLaakPr6ek6fPr3gsePHj3P8+PG5v0ejUerq6mhqarL1PiGEEEJkF1tXTAVoamri8ccf584775ybzZh/6fXEJdpPnTq1YJZjtfeJ1NTV1XHmzBmvhxFIjY2NnDlzhmg0SmdnJ/fff7/c90iR3EPKPVK33pCfxfbYbkJE5jl27Bj19fVJnZMt7GlsbKSpqWmuYY5Goxw8eJBDhw7JTJ5Ny91DqqysTH4hSTOpW2/Iz2IFlvC1jo4O6/Dhw5Z8lc44fPjwoseOHz9uAVZHR4cHI/KvI0eOWEePHl3wWCQSkdp1gNSt++RnsRq5xqzPnT59mgcffNDrYQTW6dOnqa+vX/DYq171KiD+26ZI3smTJ+d2wyTMv4eUSB+pW/fJz2I10oT42OnTpzl8+LDXwwi0w4cPL/qHU9iXzD2kRPpI3bpLfharS9tl24W7otEoQ0NDVFVVyQ9wB506dWrRY2fPngXk9HI7VO8hJdRI3bpHfhanRmZCfOrEiRNypLtHGhoaaGhokCv6ppHqvadE8qRunSE/i1MjTYgPNTY2ym8zHrn//vs5dOiQnM1hk+o9pER6SN06Q34Wp052x3ikqamJD37wg0m//vOf/zwHDhyYe6/8MEleKlnPd+LECUpLS+UURwWq95ASqZO6dY78LE6dNCEeOXDgAM8884zt9504cYKOjo4FR74n9kPKNReWppr1fKdPnyYajS66OrD845kcuYeUN6RunSM/i9NDLlYWACdOnOChhx6SC+Q4pKmpibNnzy7Y7xuNRmlsbJQj4m146KGHCIfDNDQ0zD3W2dlJdXW11K4DpG7dJz+L7ZOZkACQg/qc09nZyWOPPcaDDz644P5HZ86c4aGHHvJwZP5TX19PXV3dgiZE7iHlDKlbb8jPYvtkJsTHOjs7OX78OKdPn6azs5PDhw9TV1cnR2qnUUlJybI/WOR/HfvkHlLukLp1l/wsVidNiBBCCCE8IafoCiGEEMIT0oQIIYQQwhPShAghhBDCE9KECCGEEMIT0oQIIYQQwhPShAghhBDCE9KECCGEEMIT0oQIIYQQwhPShAghhBDCE9KECCGEEMIT0oQIIYQQwhPShAghhBDCE/8/VZDnwgF2EdoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from figaro.utils import rvs_median\n",
    "\n",
    "x = np.linspace(*draws[0].bounds[0], 1000)\n",
    "med = np.median([d.pdf(x) for d in draws], axis = 0)\n",
    "\n",
    "samples = rvs_median(draws, size = 1000)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "_ = ax.hist(samples, histtype = 'step', density = True)\n",
    "_ = ax.plot(x, med)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gradient of the recovered distribution can be evaluated using the `gradient()` method of individual draws:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.12637676]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d.gradient(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to evaluate the gradient of the median distribution in a numerically stable way using the dedicated function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.09038139]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from figaro.utils import gradient_median\n",
    "gradient_median(X, draws)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note: these methods are painfully slow and we haven't really found a way of optimising them. If you manage to improve it, please open a pull request!\n",
    "\n",
    "For multivariate distributions, it might happen that one needs to evaluate the conditional distribution or the marginal distribution.\n",
    "Making use of the properties of the multivariate Gaussian distribution, we can obtain the conditional and/or marginal distribution analytically both via the methods included in the `figaro.mixture.mixture` class or via the ones in the `figaro.marginal` module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "draws = load_density('./draws_example_2d.json')\n",
    "\n",
    "d = draws[0]\n",
    "\n",
    "# Single draw\n",
    "d_marg = d.marginalise([1]) # Marginalisation over the second dimension\n",
    "d_cond = d.condition([1.2], [0]) # Condition on a specific value along the first dimension\n",
    "\n",
    "from figaro.marginal import marginalise, condition\n",
    "\n",
    "# Vectorised\n",
    "draws_marg = marginalise(draws, [1])\n",
    "draws_cond = condition(draws,[1.2], [0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that in both cases the original `draws` are preserved.\n",
    "\n",
    "A set of reconstructions can be saved using the `figaro.load.save_density` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from figaro.load import save_density\n",
    "\n",
    "save_density(draws_marg, folder = '.', name = 'marginalised')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots produced by the CLI can be easily reproduced using the methods included in the `figaro.plot` module, mainly `figaro.plot.plot_median_cr` and `figaro.plot.plot_multidim`. Please refer to the relevant documentation page for the details. "
   ]
  }
 ],
 "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}