{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "849decd0-5b2c-4ec8-8f93-9b367aecce65",
   "metadata": {},
   "source": [
    "# Custom priors\n",
    "\n",
    "In this notebook, we demonstrate how to define custom parameter priors."
   ]
  },
  {
   "cell_type": "raw",
   "id": "0b3b0120-dc05-4e5f-a6ed-21ef347ed0ea",
   "metadata": {
    "raw_mimetype": "text/restructuredtext",
    "tags": []
   },
   "source": [
    "In this notebook:\n",
    "\n",
    "* :class:`RV <pyabc.random_variables.RV>`\n",
    "* :class:`DistributionBase <pyabc.random_variables.DistributionBase>`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "91fda0a7-b8a1-4f4a-8783-afc026213edf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from pyabc import *\n",
    "\n",
    "rng = np.random.default_rng()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e389904-a723-4003-b69b-43cc7731d48a",
   "metadata": {},
   "source": [
    "We consider a simple 2-dimensional test problem with four posterior modes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f5992100-8305-4328-822f-96378c2624c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# noise standard deviation\n",
    "std = 0.2\n",
    "\n",
    "\n",
    "def model(p):\n",
    "    \"\"\"Quadratic model with two in- and outputs.\"\"\"\n",
    "    return {\n",
    "        \"y0\": p[\"p0\"] ** 2 + std * rng.normal(),\n",
    "        \"y1\": p[\"p1\"] ** 2 + std * rng.normal(),\n",
    "    }\n",
    "\n",
    "\n",
    "# ABC distance function\n",
    "distance = PNormDistance(p=2)\n",
    "\n",
    "# ground truth parameters\n",
    "gt_par = {\"p0\": 1, \"p1\": 1.2}\n",
    "\n",
    "# observed data\n",
    "obs = model(gt_par)\n",
    "\n",
    "# ABC population size and maximum evaluations\n",
    "pop_size = 1000\n",
    "max_total_sim = 50 * pop_size\n",
    "\n",
    "# parameter boundaries\n",
    "prior_bounds = {\"p0\": (-2, 2), \"p1\": (-2, 2)}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f21be6c4-d0c4-433c-b458-1a66c3612bf9",
   "metadata": {},
   "source": [
    "In most pyABC examples and applications, we have independent parameter priors, which can be expressed in pyABC simply via a `Distribution` which assumes independency of passed one-dimensional priors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1af01956-f2d4-4135-b93c-d3baca52667f",
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = Distribution(\n",
    "    p0=RV(\"uniform\", -2, 4),\n",
    "    p1=RV(\"uniform\", -2, 4),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bbb74fd-8863-44bd-a4c4-d81d632f5180",
   "metadata": {},
   "source": [
    "If we use this prior, we find that the posterior exhibits four distinct modes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "211dbea5-8fb4-4961-9232-172964af0072",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ABC.Sampler INFO: Parallelize sampling on 4 processes.\n",
      "ABC.History INFO: Start <ABCSMC id=4, start_time=2022-04-04 16:54:41>\n",
      "ABC INFO: Calibration sample t = -1.\n",
      "ABC INFO: t: 0, eps: 1.49234810e+00.\n",
      "ABC INFO: Accepted: 1000 / 2007 = 4.9826e-01, ESS: 1.0000e+03.\n",
      "ABC INFO: t: 1, eps: 1.10051564e+00.\n",
      "ABC INFO: Accepted: 1000 / 2618 = 3.8197e-01, ESS: 9.6474e+02.\n",
      "ABC INFO: t: 2, eps: 8.41417921e-01.\n",
      "ABC INFO: Accepted: 1000 / 3349 = 2.9860e-01, ESS: 9.5689e+02.\n",
      "ABC INFO: t: 3, eps: 6.29025180e-01.\n",
      "ABC INFO: Accepted: 1000 / 4853 = 2.0606e-01, ESS: 9.4498e+02.\n",
      "ABC INFO: t: 4, eps: 4.65501309e-01.\n",
      "ABC INFO: Accepted: 1000 / 7645 = 1.3080e-01, ESS: 9.7925e+02.\n",
      "ABC INFO: t: 5, eps: 3.31555988e-01.\n",
      "ABC INFO: Accepted: 1000 / 12925 = 7.7369e-02, ESS: 9.8705e+02.\n",
      "ABC INFO: t: 6, eps: 2.39407788e-01.\n",
      "ABC INFO: Accepted: 1000 / 24502 = 4.0813e-02, ESS: 9.8812e+02.\n",
      "ABC INFO: Stop: Total simulations budget.\n",
      "ABC.History INFO: Done <ABCSMC id=4, duration=0:00:36.329396, end_time=2022-04-04 16:55:17>\n"
     ]
    }
   ],
   "source": [
    "# standard pyABC workflow\n",
    "abc = ABCSMC(model, prior, distance, population_size=pop_size)\n",
    "abc.new(create_sqlite_db_id(), obs)\n",
    "h = abc.run(max_total_nr_simulations=max_total_sim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d8e0debc-3c11-493d-a88f-cf07785587f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ABC.Transition INFO: Best params: {'scaling': 0.05}\n",
      "ABC.Transition INFO: Best params: {'scaling': 0.05}\n",
      "ABC.Transition INFO: Best params: {'scaling': 0.05}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[<AxesSubplot:ylabel='p0'>, <AxesSubplot:>],\n",
       "       [<AxesSubplot:xlabel='p0', ylabel='p1'>,\n",
       "        <AxesSubplot:xlabel='p1'>]], dtype=object)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABM6klEQVR4nO29eZxkZXXw/z3VXd3Ve8/07FsPwoDAuDKCW+ICeRlBwagoMKCCyRjQnxjjmzeGvCSv70s2NYkGQSdKEBlxN2LEDSSCUQxIEId9m4HZ1963Ws7vj+fe6uruWrvr1q3lfD+f51NVdz1VdevUuec5i6gqhmEYRuWJhC2AYRhGo2IK2DAMIyRMARuGYYSEKWDDMIyQMAVsGIYREqaADcMwQsIUsGEUiYisFZG7ROQREXlYRK4KWyajthGLAzaM4hCRlcBKVX1ARLqAXwNvVdVHQhbNqFHMAjaMIlHVfar6gPd8GHgUWB2uVEYtU3MKePPmzQrYCGP87O/dCE+GqkFE1gMvA341a/lWEblfRO4/9dRTw//ObFT1tVpzCvjw4cNhi9C4HH7SjQZHRDqBbwEfVtWhzHWquk1VN6nqpra2tnAENGqG5rAFMGqIt/9L2BKEjohEccp3u6p+O2x5jNqm5ixgwwgLERHgi8CjqvoPYctj1D6mgI3i+em1bjQurwEuBd4oIg9645ywhTJqF3NBFImqMjaVpKO1gT+yoT1hSxAqqvpzQMKWw6gfzAIuktt+s5fTr72D4Yl42KKEx1uvd8MwjLJgCrhI7tt5lNGpJLuPjYctimEYdYIp4CJ5Yv8IAPuHJkKWJETu+Cs3DMMoCw3s0CweVeXxA8MAHBhsYAU8djRsCQyjrghMAYvIjcCbgYOqujHL+tcD3wWe9RZ9W1U/HpQ8C+Hg8CSD487329AW8HmfCVsCw6grgrSAbwKuA27Os809qvrmAGUoC0941i/AgUZWwIZhlJXAfMCqejdQF/esj+93Cnh1bxv7G9kF8aOr3TAMoyyEPQn3KhH5jYj8QEROzbVRZoGTQ4cOVVI+wFnASzpbOXllF/uHJit+/qohMeGGYRhlIcxJuAeAflUd8bKJ/g3YkG1DVd0GbAPYtGlT0ZWGysXjB0Y4aUUny7tjPPDcQKVPXz2c+6mwJTCMuiI0C1hVh1R1xHt+OxAVkSVhyZOLVEp58sAwJy7vYkV3jKOjU0wmkmGLZRhGHRCaAhaRFV5xE0TkdE+WI2HJk4s9A+OMTSU5aXkXy3tiABxsVDfED/7MDcMwykKQYWi3Aq8HlojIbuAvgSiAqn4OeAdwhYgkgHHgQq3C/kj+BNyJK7oYmUgALhRt7eL2MMUKBUVJKTSFLYhh1AmBKWBVvajA+utwYWpVjZ+AsWFZJ/u8CIhGjYT4XNtWvnDPM9z04kFetKYnbHEMo+YJOwqi6nniwDCre9voikVZ3u1cEI0aC/z0oRGOjE5x8b/cy/076yLC0DBCxRRwAR7fP8xJK7oA6I410xZtalgL+JznPsmnOm5maVcrl37xv3jw+YGwRTIy2L59O+vXrycSibB+/Xq2b98etkhGAUwB5yGVUp45NMqGZZ0AiAgremLsa1ALeCTZTGusg6+9/1UoyncfbOz6wNXE9u3b2bp1K7t27UJV2bVrF1u3bm1YJVwrf0amgPMwMpVgKpliaVdretny7taGLcjzz83v5d9XfIClXa0sam9h2JuUNMLDVzSXXHIJY2NjM9aNjY1xySWXVLUCCoJsf0aXXnopIkIkEkFEEBGWLFkS+udi1dDy4Ec9dMWmP6bl3TF+vetYWCKFysB4nN72KOA+k4YuTh8yZ511FnfeeWdR2/rWMMCWLVuCFKsquPrqq+f8GfkBVpmBVkeOHOHyyy8HwvtczALOg2/hdbZG08tWdMc4ODRJFUbMBYqq8tHJz/Ku/Z8EoCsWNQs4JEpRvj5jY2NcfXVj1PHYtWtX0dtOTU3x/ve/P0Bp8mMKOA8jk87Cm20BTyVTHB2dCkusUBiPJzmS6kTbFgO+BWwKuBL4bgYRobm5uWTl61OKYqplmppKi1QfHR3lyiuvDEia/JgCzsOQbwFnKOAVXjZco9UFHhiL8/eJC3l840cA3wI2F0TQZPozAZLJ+afBi0joPs+gufLKK+f1Gd1www2hfDamgPPg+4C7Z1nA0HixwANjTtn2tmX6gM0CDpps/sz5oqq8+93vrlslfOWVV3LDDTfMe/9LLrmk4p+NKeA8ZPUB+xbwYGPVgxgYn+ITzZ/jtAedH9EUcGV47rnnynq8VCrF5ZdfXndKePv27QtSvj7vec97yiBN8ZgCzkM2H/CyrlZEGs8FMTgWZy99SM8aALpjUaaSKasMFzDr1q0r+zGnpqbqZkJu+/btLFmyhEsuuaQsx0smkxX1B5sCzsPwRIKIQHvLtFM/2hShr6PxYoEHxuP8Y+ICEq/7GDD9p2RWcLBce+21tLeXv/BTuS3rMPD940eOlLeI4uc///myHi8fpoDzMDyRoLO1Ga9qZpqlXa0cGW0wF0TaB9wCNKYCFpEbReSgiOyo1Dm3bNnCtm3byn7cICzrSlNO/3gmqVSqYskrpoDzMDyRoCsWnbO8Ef2fA+NTfLrlemLf+yMAujy/eINFQtwEbK70Sbds2cIVV1xR1mNee+21ZT1eGJRqxV+0sZlnr+okeU0Xz17VyUUbc+ehVSqV2xRwHkYm4zP8vz5drc2MTDaWAh4ci7OveQ2yxHWNakQLOMxGs695zWvCOG1VU4oVf9HGZv7lLW2s740QEWF9b4R/eUtbXiVcieQVU8B58F0Qs+mMNZ4CHhiL8+2ui+F1fwqQvjNoMAs4NMqtCOqhUE8p/vG/PjNGR8tMV2JHi/DXZ8by7he0r9wUcB5GJhNZLeDO1uZ0jHCjMDA+lfb/wrQFPNRgn0MhgurgXe4stnpITfb948Vkvq3rkZKWp9cH7Cs3BZyH4YkEnVl8wJ2xZoYb0AL+6PDfwjcuAxrTBVEMqrpNVTep6qalS5eW7bilptcWQz1EQmzZsoUvfelLBbd7bjB77ZZcywHa29sD95WbAs6Dm4TL7gOeSjRWDOzgeJzDHSfCihcBpF0z5oKoDAtJQc5FPURCQHGTlH9+5wSjUzOV7eiU8ud3Zg8n7evrY9u2bYFXSTMFnIfhiThdWXzAvv9zdLJxFPDAWJz/7r8MfsfVgmhuitDe0tRQFrDXaPaXwEkisltE3leJ8wbhq62EdVdJrr/+em655Rb6+vqyrr91R4I//N44OwdSpFTZOZDiD783zq07Zl6//f393HLLLRw+fLgiJSqtHnAOnIWbyukDBlcrYnFHy5z19cZEPMl4PElv+8z32mg1gQs1mg2Kcvtq+/v7ufbaa+uyNnBnZ2fOxIxbdyS4dcdI1nV9fX0cPnw4SNGyYgo4B36UQ64oCIChBlE+Q+Pufb7lsT+Fg+3wrlsAqwlcKcrpqxURdu7cWbbjVQt+Vtx8EzM+/elPl1mi4jAXRA6mu2FkScTwLeAGmYgb8BTwyNKXw5rT08u7GjAcLwzK6autF7/vbBaaFWcdMaoM37rtzOaCiE27IBoBPw356Ev+CF7zofTyrljUwtAqQLnqQdSb3zeTWo3oMAWcA9+yyzYJ19loFvCY6/7h94PzaTQfcFiUEu+aj0rM6ofFQiz7/v7+MkpSGqaAczCcxwXhW8CNEgvsW8DH3/EH8JUL08u7G7AmRlj48a7ztYT7+/vrVvnC/O8Swr4rMAWcA78WcDYXhF+IpmFcEOPOApYXvA5e8Lr0cmtLVFl8SzhXqJXP7Op9YSuZSuB/Nvms2ebmZq644gr6+/sREfr7+0O/KzAFnIPhLC3pfWLRCE0RSSvpemdgLE5TRGh97QfgldMB712tzUzEU8STqRClayy2bNnC4cOHcyrh/v5+vvzlL1eVkqkUW7ZsYefOnajqnJjgvr4+brrpJq6//np27txJKpVi586doX8upoBzMN2OaK4CFhEXAdAwFnCc3rboHMvK0pHD49Of/vScW27f0vUVUbUomTDw/6hUFVWtWGJFqZgCzsHwRIKWpgixaPaJj87WxqkHMTgWp6c9Cre83Q2PTquIFhqZt9yNZunWE4ElYojIjcCbgYOqujHLegE+DZwDjAHvVdUHgpKnVIYn4ln9vz6NVBHNVUKLwokza5GbBRwuW7ZsMYVb4wSZCXcTcB1wc471bwI2eOMM4AbvsSrIVYrSp5G6YgyMxVnRHYPT/3DG8q4Gywg0jHITmAuiiO4B5wM3q+NeoFdEVgYlT6nkKsbu09lAXTEGfBfELLrTLojG+BwMo9yE6QNeDTyf8Xq3t2wOQRW5zsdIjlKUPp2xaMMo4MHxuCvG/qXz3PDoarCMQMMoNzVRjEdVtwHbADZt2pS7gnIZGZqIs2ZR7sDuztbGcEHEkylGJhMuC27V22ass7ZEhrEwwlTAe4C1Ga/XeMuqgpHJBN0FfMCNEAc86BXi6W2PwmnvnbHOJuEMY2GE6YK4DXi3OF4JDKrqvhDlmYFrR5TfB9wISQh+GnJP21wfcLQpQiwaaZhwPMMoN0GGod0KvB5YIiK7gb8EogCq+jngdlwI2lO4MLTLgpKlVFS1qCgIgNHJxJxC5fXEMa8Qz6L2FvjXc93Cy76fXm/pyIYxfwJTwIW6B6iqAh8I6vwLYTyeJJlSOlvnWn0+0z3R6lsBHx11CnhxRwu89OI567tizVaS0jDmSU1MwlWakTx1IHzSEQB1fvvtl6Jc1NECL5sb9G9dMQxj/lgqchaGilDAvnVc78rn6KhzLyxub4Fk3I0Muq0msGHMG1PAWUgXY88bB+xbwPWtfI6NTRGLRmhraYKb3+pGBo2UEWgY5cZcEFnwLbpifcD1zNHRKTcBB/Dyd89Z7+Kh6/tPyDCCwhRwFkbylKL0aRQf8LFMBfySd81Zbz5gw5g/5oLIQr5i7D7pvnB1rnyOjU25CAiAqTE3MuiKNTM2lSRR5/HQhhEEpoCzMFyED7i9pYmINIAFPBZ3ERAA2y9wIwM/HXl0Mllp0Qyj5jEFnIVpH3BuBSwiDVEPwvmAPV/4Ky53I4NGK0kpIptF5HEReUpE/ixseYzaJq8PWESagfcBvw+s8hbvAb4LfFFV6/JXNzKRoC3aRHNT/v+nrjqviJZIphgcj0/7gDe+fc42fknKRlDAItIEfBb4PVz1vvtE5DZVfSRcyYxapdAk3JeBAeCvcBccuKI57wFuAebOytQBwwVKUfrUe1eMAa8QT9oHPDHoHmM96W262zwLeLx+P4cMTgeeUtVnAETkq7i61qaAjXlRSMucpqonzlq2G7hXRJ4ISKbQGRyP052l+MxsOmP1XZR9RhYcwK1eKnJGLYjuxipJma2GddV0cTFqj0IK+KiIXAB8S1VTACISAS4AjgUtXFgMel2AC9HZ2pxWUvWInwWX9gGf8f452/hV0qwehENEtgJbAdatWxeyNEa1U2gS7kLgHcABEXnCs3r3A2/z1tUlg+PxrOUXZ9MZq+/OyH4hnrQP+JTz3Mgg7QMebwgLuGANa1XdpqqbVHXT0qVLKyqcUXvkVcCqulNV3wX0A58HHgd+DtwLVE3t3nJTrALuqnMfsF+KMu0DHj3iRgadjRUFcR+wQUSOE5EWnBFyW8gyGTVMsZlwNwFDwD94ry/GTdBdkGuHWmaoWB9wnTfmnGMBf91LRc7wATdFXDheI0zCqWpCRD4I/AhoAm5U1YdDFsuoYYpVwBtV9ZSM13eJSF3O/CZTyvBkomgXxNiUqx3cFJEKSFdZBsamaIs2uUI8AK/+YNbtumPNjWIBo6q345oJGMaCKVYBPyAir/TaxyMiZwD3BydWePi+zKJcEJ7/c6RIhV1rHB2NT0/AAZz0pqzbdbdFG8UHbBhlpVgFfBrwCxF5znu9DnhcRH6La27x4kCkC4HBUhRw63RBnnpUwMfGpqZD0ACGD7jHruUztuuORRvGAjaMclKsAt4cqBRVRCkKOF0TuE4n4o6OZhTiAfiml4ac4QMGl4yxb3CigpIZRn1QlAJW1V1BC1ItpBVwe3GTcFC/RdmPjU2xbnH79ILX/nHW7bpiUR4/MFwhqQyjfrB6wLMYKMkH7D6+wTr1fx7LLMQDsOGsrNt1xxojCsIwyo1VQ5tFKS6Ivo5WYDpjrJ6IJ1MMTSRm+oAHd7sxi+4215o+ldIKSmgYtY8p4FmUEgWxqMNtc3R0MlCZwmBgbFYhHoBvv9+NWXTHoqQURqfMCjaMUjAXxCwGx+O0NkeIRZsKbtvZ2kxLU4Qjo/VXD8LPgksnYQD87kezbpuuiDaRSIfmGYZRGFPAsxgcKy4NGVxR9sUdLRyrRwU8OisNGeD4N2TdNrMexOretsBlM4x6wVwQsyi2DoTPoo6WdMpuPeFbwL2Zk3BHn3VjFl3pkpTmgjCMUjALeBalKuC+jpa6dEH4E4szLODveqnIWeKAoWEqohlG2TAFPIvB8Tgre2JFb7+4o4Xdx8YKb1hjZPUBv+FjWbdtpLZEhlFOzAUxi1It4MV1awFP0d7SNHMycv1r3ZiFXznOLGAjLLZv38769euJRCKsX7+e7du3hy1SUZgFPItiS1H6LO5oYXgiQTyZIlqgiWctcWxsaqb1C3D4Sfe4ZMOMxdOdkc0HbFSe7du3s3XrVsbG3J3orl272Lp1KwBbtmwJU7SCBKoxCrXwFpH3isghEXnQG38QpDyFKKUUpY/vI623SIhjo1PpOOc03/uwG7OINkVob2kyC9gIhauvvjqtfH3Gxsa4+uqrQ5KoeAJTwBktvN8EnAJcJCKnZNn0a6r6Um98ISh5iqGUJAwfXwHXmxvi0MhkOtMvzZnXuJEFq4hWOWr1djsonnvuuYLLq/UzC9ICTrfwVtUpwG/hXbWUkobsU68W8L6BCVb1zpqMXHeGG1nosnoQFcG/3d61axeqyq5du7jkkktoampCRKpKuVSKXM1P/eXZPrOtW7dWxecUpALO1sJ7dZbt3i4iD4nIN0VkbZb1iMhWEblfRO4/dOhQELIC0wq4t4hKaD59dWgBT8STHBmdYmXPrKSKA4+4kYXutijDdVoVrprIdrsNkEqlAOf/vPzyy6tCuVSKa6+9lvb29hnL2tvbufbaa4HqdlGEPWv0PWC9V9D9J8CXsm1UqU6z87GA/WI19ZSMcXDI1bZYMTsc7/b/6UYWrCJaZdi1q3Bl2KmpKf7oj/6oAtJUB1u2bGHbtm309/cjIvT19dHW1sall17K+vXrc35mu3btCv2OIUgFXEwL7yOq6ley+QKu80ZolFKK0mdRewsi9aWA9w2OA8yNh/4fH3cjC91t5gOuBE1NhWuUAIyMjHDllVcGLE31sGXLFnbu3MmXv/xlxsfHOXLkSNrdkI+w3RFBKuCCLbxFZGXGy/OARwOUpyDzsYCbIkJvW7SuFPD+IdfdYo4CXn2aG1nojllfuEqQTCaL3vaGG25oKFcEwFVXXZXVRZOPsbEx3vOe94TyWQWmgFU1AfgtvB8Fvq6qD4vIx0XkPG+zD4nIwyLyG+BDwHuDkqcYfAVSShww1F89CL+90IrZPuB9D7mRhe62ZoYmEqhaTeAg6ejoKGn7Sy+9tOpm/oNi+/btHDlyZF77JpPJUCzhQBMxsrXwVtVrMp5/DMie3xoCpZSizKSvzhTw/sEJumLN6ZZLaX7ofVWzakGAK8iTTCljU0k6Zu9nlIUrr7yS0dHRkvbx/xB37drFZZddBlR/csJ8Weikmj8xV8nPJ+xJuKqilFKUmSyuMwW8b3A8ez2MzX/jRhasHkRwnHXWWYgIN9xww4KOE4/Hueqqq8okVfXgx/gWM0FZiFwxxUFhpkoGpdaB8Fnc0cKvdw2UX6CQ2Dc4Mdf9ALDyxTn38SuiDU8kWNkTlGSNx1lnncWdd95ZtuPN9xa9WpmdhrxQcsUUB4VZwBksRAEfG5uqG//nvsEJVnZnsYD3/NqNLGQWZTfKRzmVbz2SKy56vpxwwgllO1YxmALOYP4KuJVkSusiDnYqkeLwyOTcGGCAH1/jRhbSFdHq1AUhIhd4E8YpEdkUtjzzpa+vL2wRykq5XQZ33XVXWY9XCFPAGcxfAXvNOcdq3w98cHgCVeamIQOc8wk3stDtV0Srgz+hHOwA3gbcHbYgC+Gd73xn2CKUlXK7DPyMwkphCjiDofE4PSWkIfssTrenr/3uyPtzhaABLD/FjSx01fkknKo+qqqPV/q8Z555ZlmP94UvfKGuwtHOOeecsh+zkp+PKWCP+ZSi9EnXgxipfQvYjwHOGgXx3K/cyEK6JnCD+4DLXbfkjjvuoLm5fHPl9RYJcfvttxfeqEQqGQ9sCthjPqUoffx6EMfqwAUxbQFnUcB3ftyNLMSiTbQ2R2q6KLuI3CEiO7KMoqv4BVG3pLW1tfBGJVBPkRBBhI1VslCPhaF5zCcN2aeeKqLtG5ygo6WJrmzJFG/5p7z7drfVdjqyqp4Vtgyz2b59e8nJF43EunXryhL/O5tKxQObBezhK885bXiKIBZtor2liaN14ILYPzTOyt42RGTuyiUb5rQjyqQ71myt6ctMUJZYvfiBs5WiLAeVigc2BezhdzZeuzjL5FMRLGpvqYsoiL0DE7m7Qu/8uRs56KmzokSZiMjvi8hu4FXA90XkR5U4bxDWHQSn2CuNX4qynGTWEg4aU8Aeu444Bbxm0fz+Tfs66yMdef/gBCuyJWEA3PU3buRg7eJ2njtavqD4akJVv6Oqa1S1VVWXq+rZQZ9z+/bt2e9EykClU26DZMuWLfT395flWCLCtm3bKlYPwhSwx64jY6zojpVciMenHupBJJIpDg7nsYDPv86NHKzv62Dv4DgT8eJLJhq5ufrqqwPLrqx0ym3QlMtiXbx4sRXjCYPnj46xbvH8fUlLOlvTEQS1yqGRSVKaIwYYYPFxbuTguCUdqFK3VnClCcpKreQtdqUol9I8evRoWY5TLKaAPXYdHWVd3/wV8AtXdHFweJLDI7WbjJE3Bhjg6bvcyMH6Ja5W7bOHbda+HARhpfb391f0FruSlCNpxYrxhMBEPMmBoUn6F2ABn7rKlQB7eO9QucSqOHuOea2IsqUhA9z9STdycFyfU8A7TQGXhXLP8IsIO3furEvlCy5pZSFKuKWlpeJ3BqaAce4HYEEW8CmrugHYsWewLDKFwSP7hog2CcctydF14W2fdyMHPe1RFrVH2XnEFHA5yGw2WQ7qze+bjTvuuANVLVkR9/X1ceONN1b8z8kUMNMREAvxAfe0RVm3uJ2H99auAt6xZ5ATl3fR2pxjIrJnjRt5WL+kw1wQZcRvNrnQybh69Pvm44477uCKK64gEsmv4trb27nllls4fPhwKHcGpoCBXZ4F3N9XWr+t2Wxc3V2zLghVZceeQTauylNN/ck73MjDcX0d7Dxsk3BBUKwl3N/fzxVXXJFu017Pft98XH/99SSTSW655ZYZZTh9pVwNn4ulIuNcEJ2tzSyaRyW0TE5d1cPtv90/77KWYbJ3cIJjY3E2rsmjgH/+j+5xQ+6M3fVLOvj2f+9hfCpJW8v8QvqM7Fx77bVzuj+0t7eHrkSqnS1btlTt52MWMLDryCjrFrcvOOh942qnvB6pQSvY911v9HzZWXnHjW7kwY+EMD9w+cn0CTeyZVtPmAWMc0GcuKxrwcc51VNeD+8d5FXH11bngR17BmmKCCevzKOAu5YXPE5mJETeYxnzopqtOaN0Gt4CTqWU3UfH6V9ABITPks5WVnTHatIPvGPPIBuWdebPBHz8B27kYf0S9zk+axawYRSk4S3g/UMTTCVTCwpBy2Tj6u6aC0VTVX67Z4jXnVigfu0vvDTkk96Uc5OuWJQlnS0WC2wYRdDwCrgcIWiZnLKqh58+drCmJqH8DL4XrS7gMnjnzUUdb71FQhhGUTS8C8JPwuhfvLAQNJ+Nq7pJKTy6v3bcEL/d7U3Arc4TAQHQ0edGAdYv6TAXhGEUQcMr4F1HR2mKSPYuwPPAV2L3PVvZoh4LYcfeQUSms/ly8shtbhTguCUdHBqeZGTSirMbRj5MAR8ZY3VvG81N5fkoVvbEOOO4xdzws6c5ViPlKXfsGeT4pZ20txTwSP3q824UwE9lNj+wYeSnoRVwIpniod2DZYmA8BER/s/5pzI8keATP654F/OSmYgnefD5gfzxvz4XfcWNApywrBOAXzx9eKHiGUZd09AK+Nb/eo7njo6x5YzyFDvxeeGKbt79qn5u/a/neGj3QFmPXW4+97OnOTwyxQWb1hbeONbjRgE2LOvkdzYs4bN3Pc1AHbRpMoygCFQBi8hmEXlcRJ4SkT/Lsr5VRL7mrf+ViKwPUp5MBsfj/MNPnuCM4xZz9qmFEwxK5Y9/70T6Olr5i3/bwfBEdXYK3nl4lOv/42ne8pJVvOaEJYV32PEtNwogIlx97skMT8T5zJ1PlUFSw6hPAlPAItIEfBZ4E3AKcJGInDJrs/cBx1T1BOAfgb8LSp7ZXPfTJxkYj/O/33xKIH23umNR/s95p7JjzyCb/+ke/qvKJuVUlWtue5iWpgh/ce7Jxe10341uFMELV3Tzrles5cv37rTqaIaRgyDjgE8HnlLVZwBE5KvA+cAjGducD/yV9/ybwHUiIpqn9t5UIsVzR0qLMVUUVUiqsufYOE8cGOamX+zkHS9fUzj0agGc++KVrOh5NR/5+oO8a9svOfuUFZx+3GJesraXrlgzzRGhKSIIwTRezERRUgqTiSSP7RvmF08f5u4nDnHNm09hea4mnLPZ8o2SzvnHv3citz24l498/UEuOG0tJ6/soqctSkQEEUp+3+VKljGMaiFIBbwaeD7j9W7gjFzbqGpCRAaBPiDn7M3jB4b53U/kbotTLP197Xz07JMWfJxCnNa/iNs/9Dt86sdP8KOH9/PDh/cHfs5iaIs28ZaXrOLdryrB/91SmgJc1hXjf7/5FK79/qP8+Xd+W6KEc9n5t+cu+BiGUU3URCaciGwFtgIsXbOeT13wkpKPEYk4i2t5d4wTlnWypLMlsJbfs+lobeaat5zCNW85hf2DE+zYM8hEIkkiqSRTwXS9zUaTZ3FvWN7JCUs7Sw+9+83X3ONL3lX0Lheevo53vWItu4+N8+i+Icamks4aT5V2asOoR4JUwHuAzKn1Nd6ybNvsFpFmoAc4MvtAqroN2AawadMmfftp+bsyVDMremKsyNX0stp5wEtFLkEBg5uUW7u4nbVlSvc2jHohSAV8H7BBRI7DKdoLgYtnbXMb8B7gl8A7gJ/m8/8aIfPufwtbAsOoKwKLglDVBPBB4EfAo8DXVfVhEfm4iJznbfZFoE9EngI+AswJVTOqiKaoGw2IiHxCRB4TkYdE5Dsi0hu2TEbtE6gPWFVvB26fteyajOcTwAVBymCUkf/e7h5f1pAFwX8CfMybLP474GPA/wpZJqPGaehMOKNEHvyKGw2Iqv7Yu6sDuBc3p2EYC0JqzeUqIsNAmEUWlpAnTM7OHygxVd0Y0rnTiMj3gK+p6i1Z1qUjdoCNwI5KyjaLRr5Wwj5/UddqLSrg+1V1k52/8c4f9LlF5A5gRZZVV6vqd71trgY2AW8rNGHcyN9Vo5+/2HPXRBywYVQCVT0r33oReS/wZuBMi9YxyoEpYMMoAhHZDPwp8DpVtX5LRlmoxUm4bXb+hj1/mOe+DugCfiIiD4rI54rYp5G/q0Y/f1HnrjkfsGEYRr1QixawYRhGXWAK2DAMIyRqUgGHnRYqIheIyMMikhKRioS5FOouEvC5bxSRgyISSkyriKwVkbtE5BHvc78qDDnmg12rlb1WvfOHdr2Weq3WpALGpYVuVNUXA0/g0kIryQ7gbcDdlThZkd1FguQmYHMFzzebBPAnqnoK8ErgAxV+/wvBrtXKf1c3Ed71WtK1WpMKOOy0UFV9VFUrmY2X7i6iqlOA312kIqjq3UBoPZVUdZ+qPuA9H8YVd1odljylYNdqZa9VCPd6LfVarUkFPIvLgR+ELUTAZOsuUhMKqNx4jVtfBvwqZFHmg12rDUQx12rVJmKUkBaaALaHcX6jsohIJ/At4MOqOhS2PD52rRqzKfZarVoFHHZaaKHzV5hiuovUNSISxV3Q21X122HLk4ldqzOwa7WEa7UmXRAZaaHnNUhaaLq7iIi04LqL3BayTBVDXPO+LwKPquo/hC1PKdi1atdqXlS15gbwFM7P9KA3Plfh8/8+zrc1CRwAflSBc56Dm0V/GndrWcn3eyuwD4h77/t9FT7/awEFHsr4zs+ppAwLkN2u1cp/5qFdr6Veq5aKbBiGERI16YIwDMOoB0wBG4ZhhIQpYMMwjJAwBWwYhhESpoANwzBCwhSwYRhGSJgCNgzDCAlTwIZhGCFhCtgwDCMkTAEbhmGEhClgwzCMkDAFbBiGERKhK+BabrhoGIaxEEKvhiYiK4GVqvqAiHQBvwbeqqqPZNt+iazQl8nvVFRGozA/SX1Dsi0/+w0deuRoMus+v35o8keqGmazz0DZvHmz/vCHPwxbDCOTn/29e3zdnwZ9pqy/h9mE3hFDVffhaneiqsMi4jexy6qAp5iqoHTGQjl8NMEvfpi9JVhs1bNLKixORTl8+HDYIhizOfxk2BLMIHQFnEmuJnYishXYChCjvfKCGfNGgRRWc9qoEt7+L2FLMIOqUcD5mtip6jZgG0C3LLZfcw2hKHHN7oIwjEanKhRwNTdcNBaGAnFSYYthGI6fXuse33h1uHJ4hK6Aa7nholEYBeJqCjgokillPJ6kszX0n3JtMFRdDZpDD0MDXgNcCrxRRB70xjlhC2WUB0WJ5xjGwvnX/3yW13/iP0im7PMsirde70aVEPrfpqr+nCJDNozaQxXiphsCY9eRMQ6PTLLzyCjHL+0MWxyjRKrBAjbqGEWIa/ZhLJzhiTgAj+4bKrClAcAdf+VGlWAK2AgUBaaIZB3GwhmaSADw2L7hkCWpEcaOulElhO6CMOqflFm7gTE0bhZwSZz3mbAlmIEpYCNQUghTNIUtRt0y7FvA+80CrkVMARuB4sLQzN0QFEOeD3jPwDiD43F62qIhS1Tl/MiL/z372nDl8LBfhhEobhKuOeswFs7wRILjl3YA8Ji5IQqTmHCjSrBfQT0j3v9riIkQqsKUmgsiCBLJFCOTCd5y3CqePjTKY/uHOeMFfWGLVd2c+6mwJZiBKWAjUFwqsingIBiZdP7fE5Z1sqg9ymP7zQKuNUwBG4HiuyCM8uNPwHXHmnnhim4esVC0wvzgz9zjm/42XDk87JcB6Vt1icwMl1I/vbPGahlIs5uIkaj39WakqWrCTdposjIVyhRzQQTFoBeC1t0W5YUru/jqfz1PMqU0RSzsr1YwBWwETsqiIAJh2gKOcvLKbsbjSXYdGeUFlpKcmyqxfH0aWgFHWloAkLY29xhr9VY4haGTk+5xYnLm6wpZj6XiW75NPd1uwaIeb0WGRTTkblNTA85fmJoMdkY4VQMWsIjcCLwZOKiqG7Osfz3wXeBZb9G3VfXjFRMwB34IWlesmZNXuO/8sf3DpoBrCDNNjEBxxXiqPgztJqBQb7p7VPWl3ghd+cJ0FlxPW5QNyzuJiIWiFeT7f+JGlWAK2AgUNwnXlHUUopiO2eL4jIg8JSIPicjLS5ZR9W6gegoEFEmmCyIWbWJJZysHhydDlqp6mUqkeHYghTbHwhYlTVWZIZXAv00HiCzqBSC5djkA46ucKyLV7G7ZW4acq6F1j7Mq5MAht35oBJie0AoNb/LQd6VEet1taOK4lQAMnOT65yVbp10QPc8uAiD2yF634JBrHJmaCqbZqa+A50kC+JPMjtki8pNZHbPfBGzwxhnADd5juXmViPwG2At8VFUfDuAcJeG7IDpj7mfcFWtm2AtNM+by+Z89zad+exbfufLVvCxsYTwaTgEblcWlIs/vMiuyY/b5wM2qqsC9ItIrIiu9fcvFA0C/qo54zQL+Dafw55DZQHbdunVlFGEuQ+MJOlub01EPnbFo2io2ZnJkZJLP3/0MAIeq6C7BXBBGoCzEBZFJro7ZOIX8fMbr3d6ysqGqQ6o64j2/HYiKyJIc225T1U2qumnp0qXlFGMOwxNxumLTf27dsWZGJkK+K6tS/vmnTzEymeCvm/+F4++tjn5w0IAWcMSPdABSq90P5ODpXQAMvNhZD5Eudzsu+5yvaOkDLr2z9zfeMeJuu9SYiw8OKypCmpwSiyzuBSB+gnM97D+jjXOaHuKqzjtYGR9iX7SLf1jxO/z7opOZ+qV7r6uOufcU8aIi8N5TuWOeVfO6IJaIyP0Zr7d5HbBnkK9jdiUQkRXAAVVVETkdZ7gcqbQcsxmaiNMdm3apdcWa2T9YPXUOqoVdR0bZ/qtdnP/SVQzs6GJQusIWKU3DKeBG4Jymh/g/Ld+jLe6sodXxYf7f7h8D8CNOr6gszgWRUwEfVtVN+fYvomP2HmBtxus13rKiEZFbgdfj/hB2A38JRAFU9XPAO4ArRCQBjAMXei6PUBmeSNDdNv0T7mxtNhdEFj714ydojkS4+pyTee2Oi7ls5XpKnqkNCFPAdciHo3fSJjNvRds0wUf231NxBQxCcp6JGEV2zL4N+KCIfBU3+TZYqv9XVS8qsP464LpSjlkJhibiLOuantHvikXT9SGMaX7x9GHOffFKlnXHWNQe5dhoMBPO86FxFLCfbtw67YKYWOqiBIZOdLfdb9n0IAC/2/MYADfvfTUATw6+AICOvS7AveWwi5Zg3L/dC8kF0eJuP7XPJVwMnOB+jCsig1m3XxUfZsLzXCbb3VcfiQZ7CRSwgAvhd8z+rYg86C37c2AdpK3T24FzgKeAMeCyBYhbUwyNJzhh6UwLeGQyYenIGagqA2NxlnW53/3/5Xo6djYDXw9XMI/GUcANxN7mbtYk5rpKj3oZf5VEERLzVMDFdMz2XAEfmNcJapzhiTjdbTN9wACjU4kZvuFGZnQqSSKl9La7z2O4ZTmD4XuP0jRcFIRqKj1QdSMpkBQG4m0MxNs4mujkaKKTsUSUsUSUSBIiSZBECkmkIOWNsJAISARpaUFaWtDWqBtNgjYJf9/9BsZk5n/rRKSJfz3uNJrHoXkcmiaTNE0mXaGeVHAXpCrEU5Gsw5g/qsrQRGJGFIT/3PzA0/juht42Fyv/05V/yOciF4Yp0gzMAq5DvtvhyhlcM/5j+sbHORTr4OYXvIyfrXzBzAjaCrDARAwjB+PxJMmUzoqCcM9HTAGn8SvG+RZwb3uUgbHqCdULXQEXKoRizI/vdmxk+FVO8R2eCK84iyIkUqaAy83QuFOyXRkKuLPVt4CrR8GEja9se9udBbxlz//jFVPDpFJnEakCP3noChhXCOU64OZAz+LFt/qVzQBaD4wB0Pt4LwA/bz4ZgHt6XJJT8z73pS170u0bPTQy4xh+/WANq/WP58uKjLjJwO6dbhIu2erkvmvC/Z9JhlhLd7p9mgbcew8qBTlTRGvKWX78NOTMMLS0C8IiIdIMjHsuCM8CHu9+AU/vO8jwRIKe9vD95KH/Mmq1EIpRHL4FnG0Y88e3cmcnYrh1poB9jvkWsDdZuWvjB/jn5Ns4OlYdoWihK+BiEJGtInK/iNwfp3ryuI3iSCFZhzF/pl0QmRawN9NvLog0g56i9a3dRR3u7vBYlSjganBBFMRLT90G0C2LFzRl7xdVB2jaexCAJQ84a6zjgIsLjne4L6v1qLuQ25895vY94CqH6ZR3gXupwOJFEWiiQi4I350yPu5eH3JZsbER51ZYftRVPOvaPdf3277Thaept0/6vQTkPlEwazcApl0Qc33ANgk3zcBYnPaWJlqb3TX40ns/wj9HD3FsNG8CZsWoCQvYqF1UhYRGsg5j/gxNzLWA21uaaIqIuSAyGBiPp90PALLiRTySWp92TYSN/QqMQHEWcCTrMOaP3w0j0wcsIulsOMMxMDZFjxcBARB53Ue4IXle1aQjh+6CyFYIRVW/GNT5MiuXJY+5lN2Idxveua/DrfArpnmpxur3UfNv+X3ZPRdEWF2TU14FM/EKxMvI6IzHzgPe+8noCafDbtvUiBfREXAlNzcJZ8q23AxPJGhpjhCLznTvdLY2p90ThnNBLMqIduhqbaY5IuYD9ilUCMWocRRzNwSAK0U59+fbFWs2H3AGA+NxTlw+PRciX7+Uba0H+cnYJ0OUaprQFXCY+C2Fkp6FK6NuEsu3bH3rMG0lzrJ01U/hDckCTk/GeZN/6dlJzzKOzLLYYdpqrpTMvgvCKC9D4/Gs9R66YlaSMpOBsTg9bdMuCNaczpPPPMOx0eq4S2hoBWwEjyIkTQGXneGJBF1t2RRwtKpa7oSJqjI4PpVOwgDgNR/izh2/hCpxQdgvwwgciwMuP7lcEK4oe3VYd2EzOpUkntQZPmCARe1RBqpEAZsFDHNv5YvtdhyW66EQnlxBpxkXgypmAQfA8ESCVT1zy4uaC2IaX8n2ZrogvnIhVx0c4t3jHwlJqpmYAjYCxlwQQTA0PrMhp0+ntaZP4xfimVHz4QWvY+/kAQaenEJVEQn3Tsx+GUagOAtYsg5j/gzNKsbu0x2LMpVIMZkIp0tLNZEuRZn5Ob3yCp45/lISKa2KPypTwEagKJDUSNZhzI+pRIqJeCqnDxgsHRmm6z349R98FnmJGQNVEAlhvwIjYLJbv2YBzx9/kq0rRxia28YU8MBYFgv4lrfzhvuvAKqjII/5gI1AUYWU+YDLiq9cfWs3k3RXjCq4vQ4b3wUxw1Vz4mbGjo7DTqqiJKX9MozAma8FLCI3ishBEdmRY/3rRWRQRB70xjVlF74K8ZVrZx4XhKUjuyiItmjTzHTt0/+QxGmXp9eHjVnARuCoztvdcBOFu6Xco6pvnu8JapHRyXwWsPmAfY7NqgPhs9jzCR+tAh+wKWAjUBQhNU8FrKp3i8j68kpU+4xOOeXakUcBmw/YS0NunzkBx5fOoweIyPvNAjYaAAXN7W5YIiL3Z7ze5hXfL4VXichvgL3AR1X14fmIWUuMTLoQs87WuYXurSvGNIPjUzMn4AA2vg0BenZFOVoFJSlNARuBk8qtgA+r6kJaEzwA9KvqiIicA/wbsGEBx6sJfBdENgs4HYZmk3AMjMU5YdmsrjCnvReARf/xH1XRnt4m4YxAUQVNRbKOhR9bh1R1xHt+OxAVkSXzOVYRE34iIp8RkadE5CERefkCRF8Q+RRwS3OE1uaIuSBwPuDe2S4Ij0XtLVVhAZsCNgJHU9nHQhGRFeLlkorI6bjr+cg8D3cTsDnP+jfhrOsNwFbghnmeZ8H41m1HS/Yb2C5LR85eCQ3gX8+Ffz2XRe1RiwM2GgHJ5wPOv2eWbilAFEBVPwe8A7hCRBLAOHChqs6raWsRE37nAzd7x79XRHpFZKWq7pvP+RbC6GSCtqjr/5aNrli04S3gMa8S2hwf8EsvBqD7iSjDE8MhSDYTU8BGsOSfhMu/a4FuKap6HS5MrRKsBp7PeL3bW1ZxBTwymczqfvDpbG1mpMEn4Qb8OhCzLeCXbQGg67kdVTFRaS4IowJIjlF/iMhWEblfRO4/dOhQIOcYnUxkjYDwsZKUpJtuzvEBJ+OQjNMVizIymWCeN0xlwxSwETypHKO22AOszXi9xls2A1XdpqqbVHXT0qVLAxFkdDKR1wLuilln5KyV0ABufivc/Fa6Ys2k1BVtDxNzQRjBsgAXRJVxG/BBEfkqcAYwGIb/F9wkXH4XhPmA04V4ZlvAL383AF0T0/HS2TIKK4UpYCN4akABFzHhdztwDvAUMAZcFo6kLhNuWVcs5/qumLWmHxj3XRCzLOCXvAuArt/sBVzG4Mqeioo2A1PARrAoSA24G4qY8FPgAxUSJy+jk0k6lhR2QVRDx4ewSHfDmO2CmHKdz6dTtsP9ozIFbASM1IQFXEuMFDEJp55/M8zb6zAZHI/T2hyZWQkNYPsFAHS98RYAhkJ21VTFtyMim4FPA03AF1T1b0MWySgnNWAB1xKjk4mcSRjgfMDgKqI1rAIei891PwC8wpWi7K6SokULioIQkVILp2Q7RhPwWVym0SnARSJyykKPa1QJCpKSrMMonVRKGZvKHwfs3143sh94cDw+1/0AsPHtsPHtVVO0qODfo4gszrUKNymxUE4HnlLVZ7zzfRWXdfRIGY5tVAPhhlrWFX4pynyWrW/5+aFYjcjA+FR2BTwxCEBXrAMI3wIu5v7kELCLmZHz6r1eVgYZsmUYnVGG4xpVgpgCLhujXinKfBZwb5sLvTpWBcVmwmJwPMHq3iyRIre6VOT29/47TRGpfgsYeAY4U1Wfm71CRJ7Psn3ZEZGtuAIoxGivxCmNcqHYJFwZSRfiyTMJ51vAAw1sAQ+Nxzl5ZdfcFWe8HwARobM1/IzBYhTwPwGLgDkKGPj7MshQMMPIK9K9DaBbFps9VWPUQhharZCvHZGP34a9Gjo+hEVOH/Ap56WfVkPKdkEFrKqfBRCRGHAl8FqcXfNzylOS7z5gg4gch1O8FwIXl+G4RrVgCrhs5KsF7NPR0kRzRKqi4HgYJJIpRiYTaVfMDEa9aqUdfV7VuOp3QfjcDAwD/+y9vthb9s6FCKCqCRH5IPAjXBjajY3QVqZREC8KwigPI0VYwCJCb3sLxxpUAfuxvT1tWT6jr7tUZC77vpcxWOUWcAYbVTUzPOwuESlLpILXzeD2chzLqD7MBVE+/CiI9pbcPmBwfuDB8cZ0QfjRHz3Z4oBf/cH00+5YM3sGJiolVlZKUcAPiMgrVfVeABE5A7i/wD5Go1Mjqci1wnRDzvw/3UXtUY5VQdv1MPB931l9wCe9Kf3UuSDCLcpeigI+DfiFiPiTceuAx0Xkt7hU+ReXXTqjPrBp07JRjA8YoKethT0D45UQqepIW8DZFPDwAffYtbw2JuEyyNcvyzByYhZw+RidTCBS2AWxqD3Kw3sHKyRVdTGtgLNMwn3TpSL7PuCwixYVrYBVdVeQghh1jFnAZWPEqwNRSGH0tkcbNgpiKJ8F/No/Tj/tikVJFpHaHSTWEcMIFs8HnG0UopZaxVcK1w0jv/ULrhD5eDzJRDzcjg9hkNcFseEsN8gsSRmeG8IUsBEowvwVMDXUKr5SjBZoyOmTzoZrQCt4cDxOW7SJluYs6m1wtxtQFQV5TAEbwbIAC1hV7waO5tkk3Srei87pFZGV5RG8OnG1gAsr4EVeK56BBgxFGxjLkQUH8O33u0Fm1bjwLODGLBZqVJbcynaJiGSGMm7z0s6LpWpaxVeKQrWAffxmlI0YipYzDRngdz+aftpdBV0xTAEbgZPH2j2sqpsqKErNMzKZYM2iwgWp/GaUjZiMMTgez56EAXD8G9JPp10Q5gM26hXNMxZOUa3i64nRqfztiHx8H3AjpiPntYCPPusGNglnNAgLmIQrxG3Au71oiFcSYqv4SlHsJFzaB9yACngonwL+7gfdoDom4cwFYQTOfJVtLbWKrxTFTsLFohFamiMNWZIyrwX8ho+ln3a0NBGRcC1gU8BGoIjOvyNGLbWKrwTxZIqpRKooC1hEWNSAyRjxZIrRqWRuBbz+temn00XZzQI26hhLRS4PxdaB8Olta+FYg1nAeZMwAA4/6R6XbAD8gjxmARv1jCngsjBdC7jwJBx46cgN1pbIV8BZW9IDfO/D7vGy7wOEXhPYFLARLFaOsmwU05Azk0XtLTxzeCRIkaoOXwF357KAz7xmxsvukLtimAI2AscUcHkYKdUF0YA+4IIuiHUzG653xZrZNxheUXYLQzMCx5+Imz2qCRHZLCKPe4V9/izL+veKyCERedAbf1BpGYtpyJlJb3sLA2Nx3FxlY5C3EhrAgUfc8OiKNTM8aRawUa8oVe8DFpEm4LPA7+HSme8TkdtUdXbLra+p6gfnHKBCpCfhikhFBmcBTyVTjMeTtBe5T63jW/w5FfDt/9M9pn3AUUbMBxwy4m4EpGnW5IamvAed8bpqmfU+5rwfQJPJGY9Bvye/GlqVczrwlKo+AyAiX8UV+ilLz8NyUUxDzkwWZWTDNYoCLuiC+B8fn/HS74oRVlF2c0EYwaIgKc06qohcRX1m83av7vA3RWRtlvWBMh2GVlwUhN8RopGSMQbH43S0NBFtyqHaVp/mhkdXLEoipUzEw7ESTAEbgRNgKnIl+R6w3ut9+BPgS9k2EpGtInK/iNx/6NChsgowOlVqFETj1QTOmwUHsO8hNzy6Qq6I1hj3JTmQZvdFRdpi7nVXp1vR7FkYE5MA6Jhrbpgad4/p2/ewme1yaPHeT4dXLau9be4+/nsZGXWP3nsM0hVRA8q2YFEfVT2S8fILwN9nO5BXTnMbwKZNm8pq5o9MJmiOCK3ZCo1nobcB60EMjsdzh6AB/NBLRc6IAwZXE3hZd9DSzaWhFbBRAWojDvg+YIOIHIdTvBcCF2duICIrMwr9nAc8WlkR/XZEhfvB+Uz7gBvLBZHXAt78NzNedodckMcUsBEobhKuqvy9c1DVhIh8EPgR0ATcqKoPi8jHgftV9TbgQyJyHpDAdel4b6XlLLYQj09P2gXRQAp4LE5/X556yStfPONl2CUpQ1XAInIB8FfAycDpqnp//j3KcM7m6X/HSLfncli1DIDhE3oBmOpyt3hth9yX0v7sMbf93gMAJL3b97CjInzXQ6Szw73u7QEg6d1LTSyf64Jo2+Nkj+zx/JOeckxNeT/Scr+n2rCAUdXbcdXVMpddk/H8Y8DHZu9XSYptyOnT2txEe0tTw7kg8lrAe37tHr2JuLCLsoc9CbcDeBtwd8hyGAFSJ5NwoZO311kOetuiDVWUfXA8nrsOBMCPr3HDo7vN9wE3oAtCVR8FKhN/501YRWKt04uWLQHgyMsWA3Dgd92/YM+yIbf8yV4AVt3TB0CnlzET8SbhUuNeCmOlLWF/8s2bdPMt36l+J+exF7pJxcETnHUryenPd9GjzjruSziZxZuEk4R7bxrA/KIp2/Kwf2iCl6zpLWmf3vaWhmlLNJVwSSd5/6TO+cSMl30drYjAgaFw0pFrwgcsIltxbceJUbgfllFFaPX7gGsBVWXf4ASbT42VtN/ijhYOjzSGAj484oyJxR2tuTdafsqMly3NEZZ0trJvoE4VsIjcAazIsupqVf1uMcfIDO3plsX2a64haiQTruo5OjrFVCLFip7SFPDq3jbufOxgQFJVF3sGXIjl6kVZwi99nvuVe8woyrOqJ8bewfEgRctJ4ApYVc8K+hylIK3T/45TK7oAOPIi9/p//873AHhn5y4A/r/lvwfA/Qc3AhA72AtAsx8PPOXdticqq2Ek4lwKEnW3WqkeNwk3sMH9OI++ylk857/oQbc8Pn3XcE/LqQB07HP7tO13saIEVbVQqy7rrSbxK3atLFEBr13cxuGRScankrS1FD+BV4vsPjYGwJp8CvhOLxXZiwMGWNnTxlOHwinbWRMuCKO2MQt44exPK+A8yiULaxe7P9/dx8bYsLyr7HJVE3uOeRZwb57P6C3/NGfRyt4Y9zx5KJR6EKFGQYjI73vNFl8FfF9EfhSmPEYAKJDU7KMIaqFMZCXYNzQ/C3jNIqeAn/esw3pm97FxlnS2EovmsfSXbEi3I/JZ1dPG6FQylM4YYUdBfAf4TiXO5d+2E5n+h0t5KZ2pmFMGq6Iu3rcz4i7ydW1HAfiFCzIg6QVtN3u3/n4criYqHMIyKwoi2eHcCOMuqIM1q1zW7Ou6HwPgqclpF/zPYicDoP7nkJxV8S0Iced57FopE1kJ9g2M0xwR+jrzTDBlYe1iZw0+fzQcH2cl2X1sPL/7AWDnz91jRnPOlb3u975vcLzkML+FEnYcsNEALKAge7pMpKpOAX6ZyIZj/+AEy7tjNEVKu0Ve2tlKLBrh+aONYAGPFVbAd/2NGxn4bp0wOmOYD9gIFMkfhrZERDKzH7d5ES8+2cpEzuwp43i7iPwu8ATwx6r6fJZtapp9gxMlR0CAi7Ffs6i97l0QqZSyd2CCszdmC7jK4Pzr5ixa5VvAIYSiNZwC9iMXAFoOuYuy57FeAK7sucS9XuTSdQcOukmLvue87Y94F/HkZAUkzYNfKD7ufFaRURf10H7Q/ZPv3eEuwj89+nYAEpPTX3PXE85tEjvk3qNOVqAaWm5/72FV3bTAw38PuFVVJ0Xk/bgykW9c4DGrjv1DE5yyan7lutYuaqt7F8ShkUmmkqm0zzsni4+bs2hZl7uz2BdCKJq5IIxgUXX1JrKNwhRVJlJV/X/ELwCnUWe4JIxxVs3DAgYXCVHvFnBRIWgAT9/lRgZNEWFZVyt7zQIODn+SSaems4Iih9yk27L/clZh9/Puy4u39wLQGXf7dDw3DIDscQHtfi3dik++eaTbCo25iy6y/zAAfV6dka7nnBWQaPcmbDJCa2IHvGI8u917SU4Eb80vIA64JspEBs3AWJyJeIoVJYag+axd1M7wRILBsXi6Qlq9sdsLQVtbSAHf/Un3ePwbZixe2RMLxQJuGAVshITmdUHk37VGykQGzXyTMHzSkRDHxuhp7ymbXNWEr4BX5YsBBnjb57MuXtnbxiN7h8otVkFMARuBs5BMuFooExk0+4eccpnPJBxkxAIfHWPj6vpVwH0dLYWbj/asybp4VU+MOx45UPFkjMZRwP7EVYYLIjUwCEDTpFvW8bwXwO3F+foxw+q5HJJDLl0xLNfDbNI1fI8NACCjziXRust9rTG/K3JG6JI/cZfy3BeBvxc/EcOYNwu3gOs/GaOoEDSAJ+9wjxtmVkhY2dPGZCLFsbE4i724+krQOArYCAVBETUFvBD2D04QERfTOx962qJ0x5rrOhJiz8A4L1xRRKr1z//RPc5SwH4o2t6BcVPARp2RsmIQC2HvgEvCaM7Var0I6jkSQlXZc2ycs05eXnjjd9yYdXFmMkYl3TQNp4BndDT2YmCTU9lvw/305XTUQbV0Q55F2hURLyGXvVJF5BcwCWc49g+Nz9v/67N2UTtPHhwuk0TVxaGRSSYTqeJcEF3ZlXRmOnIlsThgI2DUWcDZhlEU+wYn5u3/9Vm7uI3dx8bROnQH+REQRSngx3/gxiyWdLQSbZKKxwI3nAWcybRFm92yrblLNeQmoVmxSbgFoarsH5zg9ScuW9Bx1i5uZzKR4tDwJMu6F6bMq43pMpRFdMv5hZeKfNKbZiyORIQVIcQCN7QCNiqDmLU7b4YmEoxNJRduAWeUpaw3BexbwHk7Yfi88+acq1b2tFW8HoS5IIxgUXUlL7MNoyC+Reb7KOdLf59TwI/sqz8/8O5jYyxqj9LZWoQ92dHnRhbCaE1kCtgIHvMBz5udh13kwkIt4OOWdHDCsk6+88DucohVNaRSyt1PHiq+UNEjt7mRhVW9bewfnGB0snKF2U0BG8GiLKQYT8Pzjfufp6+jhVNXLSw0SkR456Y1PPDcAE8dDKf/WRDc89Rhnj86zrtesa64HX71eTeycObJy0mklG9V8E/KFLARMAqpZPZh5GXn4VF++vhBtpyxLn+bnSJ568tW0xQRvvnr+rGCv/KrXSzuaOHsU4uIAQa46CtuZOG0/kW8dG0v//qfO0lVyEAwBWwEi2I+4Hly0y920hwRLnllf1mOt6wrxhtOWsa3HthNog4+/wNDE9zx6EEuOG0Nrc1F/kHFetzIwfteexzPHh7lp48dLJOU+TEFbARMbcQBF9H8s1VEvuat/5WIrA9SnuGJON/89W7OfdHKskYtXLBpDYeGJ7n7yUNlO2ZYfP2+50mmlItOL9L9ALDjW27kYPPGFazsifHFnz9bBgkLYwrYCBYFksnso0rIaP75JuAU4CIROWXWZu8DjqnqCcA/An8XpEzfuH83I5MJLnvN3A4OC+GNL1zGks4Wbrn3uZq2gpMp5av3Pc9rTuhj/ZKO4ne870Y3chBtivCeV6/nl88c4eG9g2WQND8WB2wEjNaCuyHd/BNARPzmn5ndl88H/sp7/k3gOhERzZNaNplI8fShwhNermmIkkgqzx8b48cPH+CHO/bx8nW9vGRt7/zeUQ6iTREuPn0dn/npU7z+k//BZa85jtP6F9HSFKGlOUIFKzGWhKoSTyoT8ST3PHmYb/56N3sGxvmLc08u7UBbvlFwk4tesY7P3PkkF37+XjZvXMGbXrSCvg7X7r4pIkV9Rscv7SxKHFPARrBo9dbQyKCY5p/pbbxC8YNAH3A410GfODDMmZ/6WcnCdMeaOfvUFXzgjSeUvG8xfPisEzl1dQ9fvOdZ/u+/P1J4hyrk1cf38dGzT2JzoSacs2kpnC3X0x7lq1tfyc2/3MUPduznG/OYtNz5t+cWtZ0pYCNYVKvK3RA0IrIV2AqwbM16Pn3hS4varzkSoSkCve0tnNa/iOgCKp8VIhIRzj51BWefuoLH9g+xb2CCyUSSyUT13qmICC1NQnMkwgtXdhVuvpmL33zNPb7kXXk3e/GaXj55QS//760b+e/nBhibSjCZSBEv891cqApYRD4BvAWYAp4GLlPVgTBlMgKg+gvAFGz+mbHNbhFpBnqAI7MPpKrbgG0AmzZt0vNfujoQgcvFC1d088IV8+u2XJM84KUiF1DAPrFoE686PnvmXDkIexLuJ8BGVX0x8AR13lqmMVE0mcw6qoh0808RacE1/5ydLnUb8B7v+TuAn+bz/xpVyrv/zY0qIVQFrKo/VlU/7+9enOVh1BMLjIKoRHiYdw36zT8fBb7uN//0Gn4CfBHoE5GngI8Ac2QxaoCmqBtVQjX5gC8HvpZtRaZfLcY8fT9GKKjqvK3djPCw38NNjN0nIrepaubMUTo8TEQuxIWHFXd/OVPOQs0/J4ALSn8XRlXx39vd48u2hCuHR+AWsIjcISI7sozzM7a5GtdWfHu2Y6jqNlXdpKqbosyvL5YRHgtwQaTDw1R1CvDDwzI5H/iS9/ybwJlSyba2Rm3x4FfcqBIkbDeWiLwXeD9wpqoWbFolIoeAUfKE/9QQS6if9/GYqm6evUJEfuitz0YMyCzAus2bxPL3fQewWVX/wHt9KXCGqn4wY5sd3ja7vddPe9uE/rmKyDDweIgihH19NfL5Y6q6sdBGYUdBbAb+FHhdMcoXQFWXisj9qropWOmCp87exxzlC5BreYPweJjfb9jXVyOfX0TuL2a7sKMgrgO6gJ+IyIMi8rmQ5TGqi1LCw8gXHmYY1UioFrCXV28YuUiHh+EU7YXAxbO28cPDfomFhxk1RjVFQZTCtsKb1AT2PvLgpfz64WFNwI1+eBhwv6rehgsP+7IXHnYUp6SrhbC/Xzt/lZ879Ek4wzCMRiVsH7BhGEbDYgrYMAwjJGpWAYvIJ0TkMRF5SES+IyK9YctUCoVSbKsdEVkrIneJyCMi8rCIXBW2TNVI2NepiFzgfT8pEalISFbY17aI3CgiB70Y8Uqfu6TfRc0qYGq4kE+RHRiqnQTwJ6p6CvBK4AM1+B4qQdjX6Q7gbcDdlThZlVzbNwFhxZ+X9LuoWQVc44V8ikmxrWpUdZ+qPuA9H8YVsanu2oshEPZ1qqqPqmols/FCv7ZV9W5cREzFKfV3UbMKeBaXAz8IW4gSyNaBoWaVl1eB7GXAr0IWpdqptet0PtTVtb0QivldVHUcsIjcAWTrOXK1qn7X2yZvIR8jWESkE/gW8GFVHQpbnjAI+zot5vxGZSn2d1HVClhVz8q33ivk82ZcIZ9aCmguJsW26hGRKO4i266q3w5bnrAI+zotdP4KUxfX9kIo5XdRsy6IjEI+5xVbyKeKKKYDQ1XjlXz8IvCoqv5D2PJUKzV+nc6Hmr+2F0Kpv4uazYTzUk9bmS68cq+q/lGIIpWEiJwD/BPTKbbXhitRaYjIa4F7gN8CfqfCP/cKmxseYV+nIvL7wD8DS4EB4EFVPTvgc4Z6bYvIrcDrceUoDwB/qapfrNC5S/pd1KwCNgzDqHVq1gVhGIZR65gCNgzDCAlTwIZhGCFhCtgwDCMkTAEbhmGEhCngkPHiJX/lVY76mhc7aRg1j4h80LuuVURydcZuaEwBh8/fAf/o9cc7BrwvZHkMo1z8J3AWsCtsQaoVU8AVQkTWe3Vht4vIoyLyTRHpAN4IfNPb7EvAW0MT0jDmQY5ru11V/1tVd4YtXzVjCriynARcr6onA0PAFcBARrnChq0cZdQ8s6/tK0OWpyYwBVxZnlfV//Se34JLlzSMemD2tf3aMIWpFUwBV5bZed9xoFdE/Kp0DVc5yqgbZl/bVuOgCEwBV5Z1IvIq7/nFwM+Bu4B3eMveA1j9VqMWyXZtGwUwBVxZHsf1iHoUWATcAPwv4CNe1aw+XCk7w6g15lzbIvIhEdmNu7N7SES+EKqEVYhVQ6sQXnuSf1fVjWHLYhjlxK7t+WMWsGEYRkiYBWwYhhESZgEbhmGEhClgwzCMkDAFbBiGERKmgA3DMELCFLBhGEZI/P8xuJ/GVXg5ZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 1+2 dim marginal matrix plot via kernel density estimates\n",
    "visualization.plot_kde_matrix_highlevel(\n",
    "    h,\n",
    "    limits=prior_bounds,\n",
    "    refval=gt_par,\n",
    "    kde=GridSearchCV(),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6798c74c-8511-4496-ac9c-b48ec60d690a",
   "metadata": {},
   "source": [
    "Now let us assume that we know that the sum of parameters should actually be greater than 1.0. We can encode this by formulating a custom prior derived from the `DistributionBase` class. We need to implement a `rvs()` method that generates (pseudo) random samples from the distribution (here via rejection sampling), and a `pdf()` method that evaluates the density for a given parameter. The density need not be normalized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7c5152cc-a39a-4db2-b80b-7aa3d99ea96d",
   "metadata": {},
   "outputs": [],
   "source": [
    "class ConstrainedPrior(DistributionBase):\n",
    "    def __init__(self):\n",
    "        self.p0 = RV(\"uniform\", -2, 4)\n",
    "        self.p1 = RV(\"uniform\", -2, 4)\n",
    "        self.min_sum: float = 1.0\n",
    "\n",
    "    def rvs(self, *args, **kwargs):\n",
    "        while True:\n",
    "            p0, p1 = self.p0.rvs(), self.p1.rvs()\n",
    "            if p0 + p1 > self.min_sum:\n",
    "                return Parameter(p0=p0, p1=p1)\n",
    "\n",
    "    def pdf(self, x):\n",
    "        p0, p1 = x[\"p0\"], x[\"p1\"]\n",
    "        if p0 + p1 <= self.min_sum:\n",
    "            return 0.0\n",
    "        return self.p0.pdf(p0) * self.p1.pdf(p1)\n",
    "\n",
    "\n",
    "constrained_prior = ConstrainedPrior()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4718a699-57c8-4ca7-aa94-ee3fd6b00976",
   "metadata": {},
   "source": [
    "As expected, now only the upper right posterior mode is sampled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6014cd2f-98f3-4474-97ac-c83c9e8bd299",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ABC.Sampler INFO: Parallelize sampling on 4 processes.\n",
      "ABC.History INFO: Start <ABCSMC id=5, start_time=2022-04-04 16:55:21>\n",
      "ABC INFO: Calibration sample t = -1.\n",
      "ABC INFO: t: 0, eps: 1.57795072e+00.\n",
      "ABC INFO: Accepted: 1000 / 1928 = 5.1867e-01, ESS: 1.0000e+03.\n",
      "ABC INFO: t: 1, eps: 1.14261527e+00.\n",
      "ABC INFO: Accepted: 1000 / 2073 = 4.8239e-01, ESS: 9.6627e+02.\n",
      "ABC INFO: t: 2, eps: 8.62637704e-01.\n",
      "ABC INFO: Accepted: 1000 / 2437 = 4.1034e-01, ESS: 9.4258e+02.\n",
      "ABC INFO: t: 3, eps: 6.21836300e-01.\n",
      "ABC INFO: Accepted: 1000 / 2536 = 3.9432e-01, ESS: 8.8251e+02.\n",
      "ABC INFO: t: 4, eps: 4.60224825e-01.\n",
      "ABC INFO: Accepted: 1000 / 2651 = 3.7722e-01, ESS: 5.7476e+02.\n",
      "ABC INFO: t: 5, eps: 3.33140306e-01.\n",
      "ABC INFO: Accepted: 1000 / 3535 = 2.8289e-01, ESS: 7.6206e+02.\n",
      "ABC INFO: t: 6, eps: 2.35619085e-01.\n",
      "ABC INFO: Accepted: 1000 / 4980 = 2.0080e-01, ESS: 6.2725e+02.\n",
      "ABC INFO: t: 7, eps: 1.70215801e-01.\n",
      "ABC INFO: Accepted: 1000 / 7465 = 1.3396e-01, ESS: 5.8647e+02.\n",
      "ABC INFO: t: 8, eps: 1.21842793e-01.\n",
      "ABC INFO: Accepted: 1000 / 13116 = 7.6243e-02, ESS: 5.4529e+02.\n",
      "ABC INFO: t: 9, eps: 8.46962594e-02.\n",
      "ABC INFO: Accepted: 1000 / 25699 = 3.8912e-02, ESS: 3.8159e+02.\n",
      "ABC INFO: Stop: Total simulations budget.\n",
      "ABC.History INFO: Done <ABCSMC id=5, duration=0:00:48.108525, end_time=2022-04-04 16:56:09>\n"
     ]
    }
   ],
   "source": [
    "# standard pyABC workflow\n",
    "abc = ABCSMC(model, constrained_prior, distance, population_size=pop_size)\n",
    "abc.new(create_sqlite_db_id(), obs)\n",
    "h = abc.run(max_total_nr_simulations=max_total_sim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dd0cf6b1-8a87-42a5-b39a-7faccb62bbe3",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ABC.Transition INFO: Best params: {'scaling': 1.0}\n",
      "ABC.Transition INFO: Best params: {'scaling': 0.2875}\n",
      "ABC.Transition INFO: Best params: {'scaling': 1.0}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe25974e310>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8b0lEQVR4nO3deZycZZXo8d+p6n1JOkmH7AtLQGIgLBlQwAXDEhBBEUYgCAhjEOWKVx3vCJ/LON6L1xlnHHUUMEgEJCAIIojInhGiBAkQIAtryB6zdyfp6u7azv3jfatS3V1VXb1UPbWc7+fzfLqr6q16T5O3D0+f91lEVTHGGFN4AdcBGGNMpbIEbIwxjlgCNsYYRywBG2OMI5aAjTHGEUvAxhjjiCVgY7IQkSkiskREVovIKhG5znVMpnyIjQM2JjMRmQBMUNVXRKQZeBn4tKqudhyaKQPWAzYmC1Xdqqqv+N/vA9YAk9xGZcpFSSXgefPmKWDNVfvTv3nNbRzOiMh04FjgxTSvLRCR5SKy/IMf/KDT/0Z//vOfdcmSJa7/nSq95aSkEvDOnTtdh1DZdr7jtQokIk3Ag8DXVHVv79dVdaGqzlHVOfX19YUP8EAc7Nq1i927d2PlxeJX5ToAU0I+e5vrCJwQkWq85LtYVX/rOp5sRIRzzjkn+X08HicQKKl+VkWxfxljshARAW4H1qjqD13Hk4tAIEAgEKCrq4s77riD1157zXVIJgNLwCZ3z97ktcpyMvB54BMissJvZ7sOKhdVVVXU1tZSW1vrOhSTgZUgTO72biYSi1OlitcxLH+quhQoyR+2qqqKSy65JPlvFQ6HqampcRyVSWU9YJOz9jN+zNErPs2Tq7e5DsXkKJF8169fz49//GM2btzoOCKTyhKwydnWvZ10RmK8vqnNdShmgFpbWzn44IMZPXq061BMCkvAJmfNS2/iW1W/Zv2ukOtQzAA1NjZywQUX0NjYiKqyd2+fkXTGAUvAJmfxjt20sI+Nuy0Bl7JnnnmG2267jY6ODtehVDy7CWdy9peZ/5vr17xBiyXgkjZ79mzq6upoaGhwHUrFswRsctYWiiS/tndGGFlf7TgiMxhjx45l7NixAOzdu5dAIEBTU5PjqCqTlSBMzo5989+5vmoxgJUhykA8Hmfx4sX85je/sWnLjlgP2OQsFu6kjjAA63eFmDVppOOIzFAEAgHOPPNM6urqKmZcd7GxBGxy9qvR/4NX97dBexcbrAdcFg455JDk9++++y7Tpk2jutpKS4ViJQiTs7ZQhEkt9YxprGHDbruDXk727NnDPffcw/PPP+86lIpiPWCTs7/f+TMaa6u4ecw/WA+4zIwaNYqLL76Y6dOnuw6lolgP2OSsOxqjtirA1NENNhmjDM2YMYPq6mqi0Sgvvvii3ZgrAEvAJmf/ErmMPx3yDaaNbmBLWyeRWNx1SCYPVq1axeOPP86GDRtch1L2rARhctIdjREKx2ipr2b8yDriCpv3dDK9tdF1aGaYHX300YwdO5aJEye6DqXsWQ/Y5KS9M8J3q37JGet+wLQxXtJdb3XgsiQiyeS7detWli1b5jii8mUJ2OSkPRShixqCtQ1MHe1NYbUbceXv5ZdfZtmyZXR3d7sOpSxZCcLkpK0zwvei8/nACSdwSnMttVUBNuyyoWjl7uyzz6ajo8N21cgT6wGbnCTWgWhpqCYQEKaMbrAecAUIBAI0NzejqvzpT39i1apVrkMatMWLF9Pa2oqIJFtrayuLFy92FpMlYJOTtlCY71XdxiEvXA/ANBuKVlHi8Thr165l7dq1rkPp1+LFi5k+fTqBQIDW1laqqqoQES699FJ27drV49hdu3Zx6aWXIiJMnz694MnYShAmJ+2dEcI0U93cCsCU0Q0sW7sLraD94SpZMBhk/vz5yWnKxfrvvnjxYi677DLicW+IZO+Em8369etZsGABAPPnz89LfL1ZD9jkpC0U4T/iF1Nz5r8AMGFkHR3hGB3hmOPITKHU1NQgInR0dHDnnXeyadMm1yH1cfXVVyeT72CEQiFuuOGGYYwoO0vAJidtnWFG1lcnez2Ntd4fT6HuqMuwjAPxeJzOzk66urpch9LHcOzyUcgJKM4SsIjUichfReQ1EVklIv/iKhbTv7ZQhO/JzfC7LwPQWBsEsB5wBWpububqq6/msMMOAyAWK45r4Mtf/vKwfM7UqVOH5XNy4bIH3A18QlVnA8cA80TkQw7jMVm0hSLsrRkHIyYB0FDj9YA7rAdckQIBL3W888473HzzzbS1tbkNCLjllluG/BnV1dXcdNNNwxBNbpwlYPXs9x9W+81W/yhSbZ1h/jj2SviEVx9r9BNwyHrAFa25uZlRo0Y5Hyc8XL3fESNGFOwGHDiuAYtIUERWANuBp1T1xTTHLBCR5SKyfMeOHQWP0XjaQhFaGmqSjxuSJYjy7wGLyCIR2S4iK13HUmzGjx/PpZdeSn19PfF43Fld+NZbbx2Wz9m9e/ewfE6unCZgVY2p6jHAZOAEEZmV5piFqjpHVeckNhI0hdceivCFv90ED34RSOkBd1dED/gOYJ7rIIrdo48+yl133UU0Wvj/KeeydObFs6p4/7omYjc28/51TVw8q+8o3NGjR+cjvIyKYhywqraJyBK8i9x6GUUmEouzrzvK/uaDodUbB9xQUzk9YFV9TkSmu46j2B155JGMHj2aqqrCppVcJk9cPKuK2z5VT2ONN4pneotw26fqgU7uXenuGnY5CmKsiLT439cDpwNvuorHZLa305uG/PYR18DHvgXYMDTT14wZMzjllFMAaG9vL1hPOJdxu9+bW5dMvgmNNcL35tb1eG4gEzeGg8sSxARgiYi8DryEVwN+1GE8JoO2zsQ6ECk14BobhpbK7lUcEA6HWbRoEY8+Wphf51zG7U4dmX7WXu/ng8HgsMSUK2clCFV9HTjW1flN7hIL8Xz41X+Ed2vhwl9SWxUgGBBCFVCCyIWqLgQWAsyZM6eiR/PU1NTwsY99jEmTJhXkfFOnTmX9+vVZj9nQrkxv6ZuEN7T3/Kcq9Jhmmwln+tXeGQYgdtAsGH8U4C3a3VATpKMybsKZATruuOMYN24c4PVQ87m/3E033URNTU3WY65/pouOcM8YOsLK9c/0HLUxbdq0YY8vG0vApl+JHnD4Q9fBR76efL6xpqoiesAici/wAnCEiGwSkatcx1QqNmzYwC9/+UteffXVvJ1j/vz5LFq0KOsx966M8sXfd7KuLU5clXVtcb74+5434BoaGgo6CQOKZBSEKW6pawGnaqwNVkQNWFUvdh1DqZoyZQrnnnsuRx11VF7PM3/+fC6//PKsJYR7V0a5d+X+tK81Njby85//vKCTMMB6wCYHbZ0RRGDkI1fCfZcmn2+srbJRECYrEeHYY4+lqqqKcDjM6tWr83auxFKSA1VfX8/+/fsLnnzBErDJQXsozIi6amTKCTD5hOTzVgM2A/HCCy/wwAMP5G2o180338zcuXMH/L5QyN3GAlaCMP3aE4p45YeTv9rj+caaKv62t/iWJDTF6ZRTTmHatGmMGTMmb+d4+umnaW1tzTnJ33333XmLJRfWAzb9auuM0FJf3ef5htoqW4zH5CwYDDJ9+nQANm3axFtvvZWX8+SynoOIcM011zgpO6SyBGz6tb8rQlNdFdxzkdd8jTVBW47SDJiq8uyzz/LMM88MafeKTPpbz3fatGn86le/4uabbx72cw+UlSBMvzojcUY31sIhH+vxfEON9YDNwIkIF154IdFoNLmu8HC66aabWLBgQY/arojwpS99qSiSbirrAZt+dYaj1NcE4UPXeM3nDUOL5nWQvSlP9fX1ye3ulyxZwtatW4fts+fPn8/ChQuZNm0aIlJUPd7erAds+hUKx2io7jtHvqGmClXoisS9BG3MAHV2dvLaa68Ri8WYMGHCsH3u/Pnzndd3c2EJ2PSrMxLzEuzdn/WeuPRBIHVfuKglYDMoDQ0NfPGLX6ShocF1KE5YCcL0qzPsJ+DD53nN11BZi7KbPGlsbERE2LdvH/fccw979+51HVLBWA/YZBWJxYnG1StBnPDFHq81VtCi7Cb/9u/fz9/+9jfa2toYMWKE63AKwhKwySoxyiFdiaEhsSi7JWAzDCZMmMBXv/rV5I4aqopI+nV8y4WVIExWXZGUBHznuV7zJXvAVoIwwySRfFetWsUdd9xBd3e344jyy3rAJqtED7ihJgizzu/xWrIGbD1gM8wCgQCBQKDshzhaAjZZJZJrfXUQZl3R47UmvwRhPWAz3I488kg+8IEPICLE43FUteDbBRWClSBMVgdKEH3/X93gD0OzHrDJBxFBVXnggQd46KGHyrI3bD1gk1XyJlx1EH75Se/JL/wB8FZDA9hvPWCTJyLC5MmTCQQCZXlDzhKwyaoztQZ8zCU9XqurDiBiPWCTXyeddFLy+/379yfHDZcDK0GYrDpTR0EcO99rPhGhsabKasCmIPbu3cutt97K0qVLXYcybCwBm6x6lCBiEa+laKgJWg/YFERzczPHHnssRx55pOtQho0lYJNVjxLEXZ/2WorG2qqK2JjTuCcizJ07l9bWVgC2b9/uOKKhc5aARWSKiCwRkdUiskpErnMVi8ksUYKoqw7CcZd5LUVDTdA25jQFt3r1am655RbWrl3rOpQhcXkTLgp8Q1VfEZFm4GUReUpV87dtqhmwznCMgEBtVQBmf67P6401VbYWhCm4ww8/nNNOO41p06a5DmVInPWAVXWrqr7if78PWANMchWPSS8UjtFQU+XddQ6HvJaioTZou2KYgquqquLkk08mGAzS3d3Nhg0bXIc0KEVRAxaR6cCxwItpXlsgIstFZPmOHTsKHlul64zEvPIDwOILvZbCGwVhPWDjzhNPPMHixYvp7Ox0HcqAOR8HLCJNwIPA11S1z0KgqroQWAgwZ86c8psKU+Q6w1HvBhzA313Z53VvFER594BFZB7wYyAI/EJVv+84JJNi7ty5fPCDH6S+vt51KAM2pAQsIlXAVcBngIn+05uBh4HbVTWS6b3++6vxku9iVf3tUGIx+REKx7whaACzPtvn9cba8u4Bi0gQ+BlwOrAJeElEHrF7FcWjsbGRQw89FID169dTV1fHuHHjHEeVm6GWIH4FHAN8Bzjbb/8CzAbuzvZG8aay3A6sUdUfDjEOkyfJ7YgAutq9liLRAy7Hefq+E4B3VXWtqoaBXwPnOY7JpBGPx3nkkUd4/PHHXYeSs6GWII5X1cN7PbcJWCYib/fz3pOBzwNviMgK/7nrVfWxIcZkhlFnOHagBHGvPxXZXwsCvB5wNK6EY3Fqq8pvtSq8G8MbUx5vAk50FIvJIhAIcMkll1BXV+c6lJwNNQHvFpELgQdVNQ4gIgHgQmBPtjeq6lKgPCZ0l7HOSIyR9dXegxOv7vN6IjmHumPlmoBzIiILgAUAU6dOdRxN5RozZgzg7abx/PPPc9xxx9HU1OQ4qsyGWoK4CLgA2CYib/u93r8B5/uvmRKX3JATYOa5XkvRWJtYEa1s68CbgSkpjyf7z/WgqgtVdY6qzhk7dmzBgjPp7d69m6VLl/LGG2+4DiWrIfWAVXUd8DkRaQCuAT4ORIBlwNahBmfc63ETrmOX97VxTPL1xuSuGGU7EuIlYIaIHIyXeC8CLsn+FuPamDFjuOaaa2hpaXEdSlbDNQ74DuBI4IfAfwEz8W7QmRLXGUmpAd9/mddSJBZlL9fZcKoaBa4FnsCbLHS/qq5yG5XJxahRoxAR2tvbefjhh4lEsg7KcmK4xgHPUtWZKY+XiIgN0ykDneEYdYkEfNK1fV5P9oDLeElK/8aw3RwuUVu2bOHNN9/kxBNPZPz48a7D6WG4EvArIvIhVV0GICInAsuH6bONI9FYnHAsTkO1f5kccVafYxK943LtAZvSd+SRR3LwwQcX5eiI4SpBHA/8RUTWicg64AXg70TkDRF5fZjOYQossRJasgSxb5vXUiRuwtmawKaYJZLva6+9xsMPP1w049aHqwc8b5g+xxSR5FKUiQT8gD8VOXUccKIHXMYlCFM+2tvbaW9vJxqNUl1d7Tqc4UnAqrp+OD7HFJfkYuyJURCn/M8+xzRYD9iUkI985CPJVdTi8Tgi4nR/OeeL8ZjildyOKNEDnnFan2MSQ9SsB2xKgYgQDAaJRqPcf//9TJ48mY9+9KPO4imK5ShNceqxISdA+yavpQgGhPpq2xfOlJZgMEhDQwONjY1O47AesMmoTwnit/5U5JQaMEBjbdD2hTMlRUQ477zzkuWH7u5uamtrCx6H9YBNRp29SxAf/abXemmoqbJ94UzJSSTfHTt28JOf/ITVqws/dcF6wCajUO9haIeemva4xtoq9nVZAjalqaWlhcMOO8zJJA1LwCajTr+um9ySaPf73tfRB/c4rqW+mvbO4pvmaUwuqqur+cxnPpN83N7ezsiRIwtybitBmIySNWB/ujEPX+u1XkZaAjZl4q9//Ss/+9nPKNT+k9YDNhn1KUGc+u20x7U0VNNmCdiUgZkzZ9LR0ZFcVzjfLAGbjLrCMUSgtsr/Q2n6KWmPsx6wKRdNTU2ceqp3r6Ozs5NQKJTXZGwlCJNRYi3g5Eyhne94rZcR9dWEo3G6IjYUzZSPBx98kMWLFxOL5e+6th6wyagzkrIYO8Dvv+Z97TUOuKXBm1PfFoowfmTlbktkysvpp5/O/v37CQbzd01bAjYZ9diOCGDujWmPS+wZ194ZYfzI4lvyz5jBGDduXHJ7+3Xr1jF27NhhnzlnJQiTUSh1R2SAqSd6rZeW+hoA2kLhQoVmTMF0dXVx33338eSTTw77Z1sP2GTUpwSxzZ8pNG5mj+NSe8DGlJu6ujouuugiDjrooGH/bOsBm4z6lCAe+0ev9WIJ2JS7adOmUV9fTzweZ9myZUSjwzPz03rAJqPOSIzWppoDT5zx3bTHjWywBGwqw4YNG3jiiSdobGzkqKOOGvLnOU3AIrIIOAfYrqqzXMZi+gqFozTUNBx4YtLxaY9rrq1CxBKwKX/Tp09nwYIFTJgwYVg+z3UJ4g5sO6Oi1acEsfV1r/USCAgj66tpC1kCNuUvkXz37NnDs88+O6T95ZwmYFV9DtjtMgaTWZ+bcI9/22tp2Gw4U2lWrVrF8uXLaW9vH/RnWA3YZNRnGNq8/5fx2JZ6Ww/CVJaTTz6Z2bNn09zcPOjPKPoELCILgAUAU6dOdRxN5YjHle5o/MBSlAATjs54/AjrAZsKIyLJ5Pvyyy/T3d3NSSedNKDPcF0D7peqLlTVOao6Z+zYsa7DqRidvVdCA9j8stfSGFlfzV5LwKYCqSrvv/8+69atG3A9uOh7wMaNUDhNAn7Sn4rcay0I8JekLLOZcCJyIfAd4EjgBFVd7jYiU4xEJLmgu4igqjlvde+0Bywi9wIvAEeIyCYRucplPOaAxMpmPUoQZ//Aa2kkbsLF44O/I1yEVgLnA8+5DsQUt2AwSDAYJBwOs3jx4pzf57QHrKoXuzy/ySzUezcM6DMFOVVLfQ1xhf3hKCPqqvMdXkGo6hog596MMQCRSO6luKKvARs3EjXg+pqUS2TDi15LIzkduULHAovIAhFZLiLLC7WdjSk+NTU1XHHFFTkfbwnYpBXyN+Ssr07pAT/zXa+lUarTkUXkaRFZmaadN5DPsZvFJmEgfzHZTTiTVleyB5xSA/7UjzIeX6oL8qjqaa5jMJXLErBJK+0oiNYZGY8v1QRsjEtWgjBpdXQnShApCXjdUq+lkbotUbkQkc+IyCbgw8AfROQJ1zGZ8mI9YJPW29v2U1sV6LnF0BJ/KnKaccDl2ANW1YeAh1zHYcqXJWCT1oqNbRw1aSTVwZQ/ks77acbj66uD1AQDZZWAjck3K0GYPiKxOCs3t3PMlJaeL4w+2GtpiIi/HkR5zYYzJp8sAZs+3ty6j+5onNm9E/B7S7yWwcj6KusBGzMAVoIwfazYuAegbw/4uX/3vh56atr3tTTUlNVNOGPyzRKw6ePVjW20NtUweVR9zxfO/3nW942sr2bb3q48RmZMebEShOljxcY2jpnS0ndGz8jJXsugxdYENmZALAGbHtpDEdbu6OhbfgB452mvZTCivrpi14IwZjCsBGF6eH1zG0DfG3AAS//T+zoj/ezdloZq9nVHicbiVAXt/+3G9McSsOlhxYY2AI6e3NL3xQsWZX1vYjLG3q4ooxtrhjkyY8qPJWDTw4qNbRw6tjGZTHtoHpf1vamz4SwBG9M/+zvRJKmqfwNuVPoD3vqj1zI4sB6ETcYwJhfWAzZJ7+3oYFdHmOOnZUjAf/GnIh9xVtqXEz1g257emNxYAjZJf3lvJwAnHzYm/QF/f1fW908b0wjA6i17OfWIg4Y1NmPKkZUgTNLSd3YyeVQ9U0c3pD+gcYzXMmhtquWDE0fwp7dtSx5jcmEJ2AAQjcV5Ye0uTjmsNfOWKqsf8VoWH5kxllfW72Ffl5UhjOmPJWADwBub29nXFeXkw1ozH/Tiz72WxUcPbyUaV5at3T3MERpTfqwGbAD4y3u7ADjp0MwlBi6+p9/PmTNtNA01QZ57ewenz8w+bM2YSmcJ2ABe/ffICSMY01Sb+aC6kf1+Tk1VgA8fMobn3rE6sDH9cVqCEJF5IvKWiLwrIv/kMpZK1hmO8fL6PZySafRDwsoHvdaPjx4+lvW7Qqzf1TFMERpTnpwlYBEJAj8DzgJmAheLyExX8VSyl9btJhyLZ6//Ary0yGv9+OjhYwF4zkZDGJOVyxLECcC7qroWQER+DZwHrM70hnA0zoZdoQKFV/5CkShvbt3Hb17eSHVQOOHg0dnfMP83OX3u9DENTB5Vz5Ort/HhQ8cQDAQIZhpZMUBTx2QYImdMCXKZgCcBG1MebwJOzPaGt7bt46M/yLwljhmcmqoAn/u7KTTU9HM51OSW/ESEU484iF8tW89pP3xuGCI8YN33Pzmsn2eMS0V/E05EFgALAMZOns5/XDjbcUTlo7oqwOHjmjh0bFPP3Y8zee0+7+vsz/V76DfPPIKTD2slEosTjceJx4cYrDFlyGUC3gxMSXk82X+uB1VdCCwEmDNnjn72+Mw7Mpg8e8WfipxDAh5ZX828WePzHJAxpc1lAn4JmCEiB+Ml3ouASxzGY/pz2e9cR2BMWXE2CkJVo8C1wBPAGuB+VV3lKh6Tg2C11yqEiPxARN4UkddF5CERaXEdkykvTscBq+pjqnq4qh6qqje5jMXk4NXFXqscTwGzVPVo4G3g247jMWXG1oIwuVtxj9cqhKo+6f+lBrAM7z6FMcNGVNV1DDkTkX3AWw5DaAV2Ojx/McTg+vx1qjqr0CcVkd8D96nq3RleT47WAWYBKwsVWxqu/42KIQbX58/pOi21BLxcVedU6vmLIYZyO7+IPA2kG65xg6o+7B9zAzAHOF9z+IUpt/9GpRhDqZy/6McBG5NPqnpattdF5ArgHGBuLsnXmIGwBGxMBiIyD/gW8DFVtTnwZtiV2k24hRV+fnAfQyWd/6dAM/CUiKwQkVtzfF8l/TfKxHUMJXH+kqoBG2NMOSm1HrAxxpQNS8DGGONIySVg19NDReRCEVklInERKdgwF9e7h4jIIhHZLiJOxreKyBQRWSIiq/3//te5iCNXdp3adZrLdVpyCRj300NXAucDw7vQbRZFsnvIHcC8Ap8zVRT4hqrOBD4EfKXId1Cx69Su036v05JLwK6nh6rqGlUt9Gy85O4hqhoGEruHFIyqPgc422teVbeq6iv+9/vwFnCa5Cqe/th1atdpLtdpySXgXq4E/ug6iAJIt3tI0SaffBOR6cCxwIuOQ8mVXacVKJfrtCgnYgxgemgUGPbluXI5v3FDRJqAB4Gvqepex7HYdWrSyvU6LcoE7Hp6aH/ndyCn3UPKnYhU413Ui1X1t67jseu0D7tOGdh1WnIliJTpoedW0PTQ5O4hIlKDt3vII45jKigREeB2YI2q/tB1PP2x69Su05yuU1UtqQa8i1dnWuG3Wwt8/s/g1ba6gW3AEwU679l4d9Pfw/sTs9D/3e8FtgIR/+e/qsDnPwVQ4PWUf/uzC/3fYQDx2nVq12m/16lNRTbGGEdKrgRhjDHlwhKwMcY4YgnYGGMcsQRsjDGOWAI2xhhHLAEbY4wjloCNMcYRS8DGGOOIJWBjjHHEErAxxjhiCdgYYxyxBGyMMY44S8CltsmiMcYMN2eroYnIBGCCqr4iIs3Ay8CnVXV1pve0yng9Vj5SsBhN7p6K/0YyvXbmqY26a3cs7Wsvv979hKq63ERx2M2bN08ff/xx12GYdP70b97Xj30r32fK+PuQytmOGKq6FW/dTlR1n4gkNq/LmIDDhAsUnRlOO3dH+cvj6bcGq5v4fmuBw8m7nTt3ug7BZLLzHdcR9FAUWxJl27xORBYACwDqaChsYGZYKBDH1p02ReCzt7mOoAfnCbi/zetUdSGwEGCEjLbf4hKkKBFNX4IwppI5TcDFtsmiyQ8FIsRdh2EMPHuT9/UTN7iNw+csAZfaJotm8BSIqCVg41ZXJEZV20aqAsUz+tZlD/hk4PPAGyKywn/uelV9zF1IJh8UJWI1YOPY1+9fgeoXueXS412HkuRyFMRSchyqYUqbKkQs/xrHVm/Zy4j6atdh9FA8fXFTthQhoumbMYWgqmxt7+KSfYvg6e+4DifJ+SgIU/4UCNv/641DuzvCdEfjNET3Qmi363CS7LfCFERcJW0zphC2tncB8F35Epz7E8fRHGA9YJN3cYQwQddhmAq2pa0TgO5IcY1Htx6wyTtvGFogbTOmEBI94Ovid8ITxTEGGKwHbArAuwlnl5pxZ0u71wOu0TCxSGfR/D1mvxUm71SFsBbLJW8q0dY2rwd8Y/QLfOa0M2h2HE+CJWCTd95UZEvAxp2tfg8YoDMSo7muOMYDWxHO5F2iBJGuGVMIW9q6CAaEG6vuou5pqwGbCqJYCcK4E4sr2/Z2MWVUPbR7j4uF9YBNQcQ1kLb1R0QWich2EVmZ8twPRORNEXldRB4SkZYM710nIm+IyAoRWT6U+EUkKCKvisijQ/kcU3g793cTjSuHjG3iu9HL2HDCja5DSrIEXI4k0LclXgoGkWCwz/P5FPd7wOlaDu4Aem9Z9BQwS1WPBt4Gvp3l/aeq6jGqOmdQwR9wHbBmiJ9hHEiMAT64tRHwasDFwhKwyTtvMZ7B1YBV9Tlgd6/nnlTVqP9wGTB5+KM+QEQmA58EfpHP85j8SIwBPmRsI9+t+iVTXvjfjiM6wBKwyTvvJlwwbQNaRWR5SlswwI+/EvhjxlPDkyLy8iA+N9WPgG9B5lXlRWRB4mfYsWPHEE5lhluiB3xIaxNd1BCmxnFEB9hNuDIiVd7Qmo5zj0s+17RuPwCBrbu8JxrqAdDt3saRsX378h5XIgFnsHOw5QERuQGIAoszHHKKqm4WkYOAp0TkTb9HPZBznANsV9WXReTjmY5L3Tprzpw5xXOXx7ClrYv66iDjRtRycXQ+Y2fO5mDXQfksAZu886YiD++lJiJXAOcAc1U1bcJT1c3+1+0i8hBwAjCgBIy3ccC5InI2UAeMEJG7VfXSQQdvCmpreycTWupoqPGuwa5I8ezOYiUIk3f9lCAGTETm4ZUEzlXVUIZjGkWkOfE9cAawMt2xWWNX/baqTlbV6cBFwLOWfEvLlvYuJo6sp646wPeqbuO41/7ZdUhJ1gMuI8HxBwGwb0qQs6pf57q6Zxg/q52tjOSHMpdHA0czflkEgIbusPemQpQgNGsJIisRuRf4OF6teBPwz3ijHmrxygoAy1T1SyIyEfiFqp4NjAMe8l+vAu5R1ceH+rOY0rO1rZPDDx9LXXWQNprZHxjhOqQkS8Bl6Kzq1/lO/e+pFy/ZTqKd/6u/hzgs58iCx+OVIAaXgFX14jRP357h2C3A2f73a4HZgzpp5lj+G/jv4fxMk1/haJwd+7uZ0FJPbVWAH8Qu4n9MOYxi2RXOShBl6Lq6Z5LJN6GeCF/XZxxFJMQ0kLYZk0/b9nahChNH1iEi1FUFi2ocsPWAy0Cw2VvbKXTURADGB9rTHjeRdrad4I2UOHhV4XajGEoP2JihSIwBntDijf7516pbmfxePXCvw6gOsARchrbXNjG+e3+f57saA5D2llV+KULUErBxILEK2sSRdQDsDLRSH2xwGVIPloBLWWIq8dQJAFTv88oO32ucy/fDf6AhOVkMugJBbj34RBqf9B7r7raChakKkbiVG0zh7e30fidaGrzJF3c3XMqRo0ZwusugUlgCLkMPN80C4FvtS5gU20tXo5d8l050M/y8n4kYxuRNR9ir9zbWetdffXWwqPaFc5qARWQR3mD67ao6y2Us5ebhplk83DSLhafeAcDSre7m/ihCNG4J2BReKBxDBOqqvOvvW6H/INglwENuA/O57gHfAfwUuMtxHCVJgt5FJfu8wm6VP8FnyqOjksdct+GLANTu8R6Pf7kNAI0VrhfgLcZjJQhTeKHuKPXVQQIB76bzjpopRGLFMxPOaQJW1edEZLrLGEz+WQ/YuNIRjiWnIAM8PuZytrR1km5wuQuue8D98lexWgBQR/HcvTQDE6dww96MSQiFo8n6L0BddYAuqwHnLnWVqREy2laZSqFR7w5vdMMmAAI13p3e5u27ksc0bvRGSAT2eMPSou9v8N9cuD/DFKwHbJwI9eoB/8Pf/g87Qt14s9vds8KcyTtVIaqBtM2YfAqFozTUHPif/47Gw1mj0xxG1JP9Bpi883rAgbTNmHzq6I71SMAvTbmCm6PnOoyoJ9fD0PqsdKWqaRdaMf2Lh70VzuK7DpQg2O0Nf4gXsOTQm3cTzpKtKbzOcIzxI+qSj+uqg3RF4qgq/kp5TrkeBVEsNyNNPilWbjBOdPQqQXz6rf/FkdX76I7Oo67a/X0J+60odxov6A23tCEw+BJEhm3pR4vIUyLyjv91VIb3Xu4f846IXD58P5EpFaFwjIaUURC7Rs3mlfgMOsPFMRLCErDJO0WIxQNpWw7uoO+29P8EPKOqM4Bn/Mc9iMhovMXbT8TbiuifMyVqU746uqM0poyCeG/GldwWO6dolqS0BGwKIo6kbf1Jty09cB5wp//9ncCn07z1TOApVd2tqnuAp+ibyE0Zi8WV7mi8xzC0er8cUSxjgYt+HLApfapk6+22isjylMcL/bHf2YxT1a3+93/D236ot0nAxpTHm/znTIUIhb3VAFNrwB9+8SvcVt1OZ+QUV2H1YAnYFIBkS8CD3pYeQFVVRPI6QUdE6vB2U67F+515QFWLZ2dHk1bIr/Om1oD3TjiZv2xYzzlF0gO2EoTJO68HLGnbIG0TkQkA/tftaY7ZDExJeTzZf24wuoFPqOps4Bhgnoh8aJCfZQokkYBTa8C7j7qKX8bOKpqt6S0Bm7xTGO494R4BEqMaLgceTnPME8AZIjLKv/l2hv/cgKknscVItd9sWnyR6+j2ShD1NT3XggBsFISpJOl7v7n0gP3JOi8AR4jIJhG5Cvg+cLqIvAOc5j9GROaIyC8AVHU38H+Al/z2Xf+5wf0EIkERWYHX235KVV8c7GeZwkjXAz70icu5o/pfi2YUhNWATd6pQnyQM+GyTNaZm+bY5cA/pDxeBCwa1In7fnYMOEZEWoCHRGSWqqaOTU6u2jd16tThOKUZouRNuJQacPjQM3l67dvMLpIEbD1gUxDDXAN2RlXbgCX0GtKmqgtVdY6qzhk7dqyT2ExP6XrA0eOv4u7Y6UUzDM0SsCkIVUnbSoGIjPV7vohIPXA68KbToEy/EjXg1GFoB8YBF8dNOCtBmLxThHiJJNsMJgB3ikgQr9Nyv6o+6jgm04/kMLSUBNzw6/O5u3onr0TuzPS2grIEbPJPQUuw3JCgqq8Dx7qOwwxMsgRReyDNBWadzx/feYMRVoIwlSQel7TNmHwJhaOIQG1VSpo7/goeqTqjaGrA1gM2eacKausBmwLr6I7RWFPVZ93f+upg0SRg+60wBZFYFbN3MyZfOiM91wIG4Jef5OexG4tmIob1gE0BSEnXgE1p6uiO9aj/AnDMJTyz/a2iGQVhPWCTf/5NuHTNmHwJhaPU99714tj5PN90ps2EM5XGkq0pLK8H3CsBxyI0BuNFk4CtB2wKI56hGZMnoUisx2LsANz1aW7ccz3dloBNxbAShHEg1J3mJtxxl/FCyyetB2wqTFzSN2PyJBRO0wOe/TleHz2vaBKw1YBN/imIlRtMgYXC0b414HCI5mDERkGYSpKh95vbesBHiMiKlLZXRL7W65iPi0h7yjE35usnMaWjI10PePGFXLnuH+myccAgIvOAHwNB4Beq+n2X8Zg8GmSHQ1XfwtsGCH8xnM3AQ2kOfV5VzxlkdKbMRGNxwtF43xrw313J6yu20Lm7OBLwoHvAItLfzrX9vT8I/Aw4C5gJXCwiM4fymaZIKUhc0rYBmgu8p6rr8xClKSOhSN+V0ACY9VnWT5hHNK5EYu7LEFl7wCIyOtNLwNlDPPcJwLuqutY/16+B84DVQ/xcU4wy76A2kG3pLwLuzfDah0XkNWAL8E1VXTWoOE1ZCHX3XQkNgK52miXkfRuJUR10W4XtrwSxA1hPz1H06j8+aIjnngRsTHm8CThxiJ9pilSWjeNz2pZeRGqAc4Fvp3n5FWCaqu4XkbOB3wEzBhepKQcd4b6LsQNw7yXM29vFd/g6nZEYzXXVDqI7oL8EvBaYq6ober8gIhvTHD/sUvfaqqOhEKc0w00ZjiFnZwGvqOq2Ph+vujfl+8dE5GYRaVXVnUM9qSlNiR5wn5twJ17N2rW7YCt0F8FIiP763z8CRmV47d+GeO7NwJSUx5P953pI3WurmtohntK4IvH0bQAuJkP5QUTGi7/moIicgHdd7xpqzKZ0JTbkbOzdA555LrunnglQFGOBs/aAVfVnACJSB3wZOAWvP7MUuGWI534JmCEiB+Ml3ouAS4b4maZYDaGzISKNePuwXZ3y3JcAVPVW4ALgGhGJAp3ARaqauehhyl5yO6LeNeCOXYyItwMUxZKUuQ5DuwvYB/yX//gS/7m/H+yJVTUqItcCT+ANQ1tkN07Kk/ijIAZLVTuAMb2euzXl+58CPx30CUzZyVgDvv8yZndGgOuKYlH2XBPwLFVNHSK2RESGPFpBVR8DHhvq55jiV6oz4URkCl5nYxzeX38LVfXHbqMy/Um3IScAJ13Lju37YUNxlCByHYPxioh8KPFARE4Elmc53pgDdFhqwK5EgW/4HZAPAV+x8erFL9SdqAH36mMecRbhQ70acCn1gI8H/iIiidEQU4G3ROQNQFX16LxEZ8pHiVZkVXUrsNX/fp+IrMEbQmnj1YtYh98Dru/dA963jaZIYhyw+x5Argl4Xl6jMGWvRHq7WYnIdLzt6V9M81pyuOTUqVMLG5jpIxSOEgxIzx2RAR64kvGxOHBtsk7sUk4J2KZ+miEr0R5wgog0AQ8CX0sdd5zgz95bCDBnzpwS/2lLn7cUZbDPjsic8j8JxOME3ouwbW+3m+BS2HKUJv9KfDlKEanGS76LVfW3ruMx/Qv5W9L3MeM0qoBxI55hS1tnwePqzRKwyTuhdBOwP8HjdmCNqv7QdTwmNx3hNLthALRvAmBiS31RJGBbD9jkX2mPgjgZ+DzwiZT1hoe6EJXJs85wjIbei7ED/PZq+O3VTGypZ3MRJGDrAZvCKI1k24eqLsW2dC45Xg84TXr76DcBmPh2HU+s7CIeVwIBd/+8loBNQZRIb9eUiVA4xujGmr4vHHoqAJO2ryMci7Ozo5uDmusKHN0BVoIw+adZmjF5EApnuAm3+33Y/T6TWuoB2NLWVeDIerIEbAqihGvApgSl3ZIe4OFr4eFrmZhMwG7rwFaCMAVhydYUUkc41nc3DIBTvfX8Ewl48x5LwKbMiWbdEcOYYRcKR/tOQwaYfgoAI1Rpqq1yPhLCErApiKH0gEVkHd5yqDEg2nsLI3+s7o/x9ikMAVeo6iuDP6MpZeFonEhM+y7GDrDzHQCkdQYTW+qsBGEqxNBLEKdm2WLoLLw94Gbg7St4C7a/YMVqC4UBGNmQZhTE77/mff3CH7zJGO2WgE25y/9U5POAu/xdMJaJSIuITPBXMjMVZk8oAsDodAl47o3Jbye11PP6pvZChZWWJWBTEENMwAo8KSIK/DzNtvXpdtiehL+MpKkse/we8KiGNDseTz3wh9HElnp2d4TpDMfS14sLwBKwKYgsN+FaRSR1cf+FaRLsKaq6WUQOAp4SkTdV9bl8xGlK354OPwGnm4ixzV/GedzM5FjgzW2dHHZQU6HC68ESsMk/JVsNeGfvm2p93q662f+6XUQeAk4AUhNwTjtsm8qQKEGMSleCeOwfva9+DRi8scCWgE3ZGspqaP6OyAF/N4pG4Azgu70OewS4VkR+jXfzrd3qv5UrUYJoSVeCOOPApTOxxZuC7HIkhCVgk38KEh/0QOBxwEP+wtpVwD2q+nivbekfwxuC9i7eMLQvDDlmU7L2dIRpqAlSV52mrjvp+OS340bUERBLwKYCDLYHrKprgdlpnk/dll6Brww2NlNedofC6csPAFtf975OOJrqYIDxI+rY7HA9CEvApiBsKrIplLZQhFGNacoPAI97U5H5wh8A/HWBQwWKrC9LwCb/SnxLIlNa9mTrAc/7fz0eTmypZ8XGtvwHlYGthmbyzrsJp2mbMcNtT0eYlkwJeMLRXvNNbKlna3sncUfXopMELCIXisgqEYmLSNYhSKYMlPaWRKbE7AlFGJ1uBATA5pe95hs/opZITJMjJwrNVQ94JXA+PcdymjJWyglYRBaJyHYRWek6FpNdNBanvTOSuQf85I1e841uqgVgd4ebBOykBqyqawD8oUWmApRKss3gDuCnwF2O4zD9aO/014FINwsO4Owf9HjY6h+3qyPMjLxGll7R34QTkQXAAoA6GhxHYwZlaOOAnVPV50Rkuus4TP+yTsIAGDezx8PRTX4C3l9mPWAReRoYn+alG1T14Vw/x18XYCHACBldur/FFWwoM+FKRWpHYerUqY6jqVzJldAy9YA3vOh99RflSRy3u6M777Glk7cErKqn5euzTYnR8h/xkNpRmDNnTnn/sEUsuRBPphrwM/5UZH8ccGLJyp3l1gM2JlW594BNcei3BPGpH/V4WBUM0NJQ7ewmnKthaJ8RkU3Ah4E/iMgTLuIwBaJATNM3Y4ZRvyWI1hleSzGmsaayErCqPqSqk1W1VlXHqeqZLuIwhVPKEzFE5F7gBeAIEdkkIle5jsmkt6cjTE1VgPp0C/EArFvqtRRjGmvZub/MasDGpCrlXZFV9WLXMZjc7AmFGd1Qk3mI6xJ/KrJfAwavt/zejv0FiK4vS8Am76TEh6GZ0rG7I5K5/gtw3k/7PDWmqYa/rqugEoSpPBLTtC3re0SmiMgSEVntT12/Ls0xHxeRdhFZ4bcb032WqQxt2RbiARh9sNdSjGmsYU8oTMxBJ8F6wCb/VGFwF3cU+IaqviIizcDLIvKUqq7uddzzqnrOkOM0JW9PKMwHxo/IfMB7S7yvh56afGpMUy2q3ntb/anJhWIJ2BTEYEoQ/rZCW/3v94nIGrzdjnsnYGMAbxRE1hLEc//ufU1JwAcmY1gCNuVI6bfc0B9/KvCxwItpXv6wiLwGbAG+qaqrhnQyU5LicaUtFM48BA3g/J/3eWpM6nTkcfmKLj1LwKYgsvSA+92WXkSagAeBr6nq3l7vfwWYpqr7ReRs4HfgZF0V49jerghxJfNKaAAjJ/d5akyj1+vd5WA6siVgk3+JiRjpZd2WXkSq8ZLvYlX9bZ+PTknIqvqYiNwsIq2qunOIUZsSc2ASRpYSxDtPe19nHFgpIbUEUWiWgE3eCYrowEsQ4g3mvB1Yo6o/zHDMeGCbqqqInIA3smfXUOI1penANOQsPeCl/+l9TUnAoxqqEXGzHoQlYFMY8UEtBnEy8HngDRFZ4T93PTAVkjsjXwBcIyJRoBO4yN8l2VSYfhfiAbhgUZ+nqoIBWuqrnayIZgnY5N8gb8Kp6lK81SyzHfNTvMXSTYVLliCyJeDm9HfZxjTVOlkT2BKwKQAdbA/YmJwlesAt2WrAb/3R+3rEWT2eHt1Ywy6rAZuylP0mnDHDYk8oTFVAaK7Nktb+4v+x1CsBtzbV8Pa2wq8HYQnYFIRYD9jk2Z6Qtx191r0m/z79tn6jG2vY5WBFNEvAJv9UIWYJ2OTX2h0djB/Zz0y2xjFpnx7dWEtbZ4RoLE5VsHBL5NhiPKYw4vH0zZhhsLsjzEvrdnPqEQdlP3D1I17rpbWpxl8PIpKnCNOzHrDJP2Wwi/EYk5OnV28jrnDmB9PtA5ziRX8q8sxzezydOhljbHPh1oOwBGwKQCEecx2EKWN/XLmVyaPq+eDELCuhAVx8T9qnk9OR93cDzcMcXWaWgE3+KVYDNnmztyvCn9/dxeUnTct+Aw6gbmTap5ML8hR4KJrVgE0BaEnXgEVknoi8JSLvisg/uY7H9LTkze2EY3Hmzeqn/ACw8kGv9TLG0XoQ1gM2+adArDRLECISBH4GnA5sAl4SkUfSLApvHHl85d8Y21zLsVNG9X/wS/5U5Fmf7fG0N3yNgg9FswRsCqCkh6GdALyrqmsBROTXwHlkWRS+Oxp3tsljpYnGlP9+awefPX4SgUA/5QeA+b9J+3QwIIxqqGHdrtCw/NsdOrYpp+MsAZv8U9AS7QHj7cCxMeXxJuDEbG94e9s+5v7Hn/IalOnprFkTcjuwpiHjS+NH1PHIa1t45LUtQ45n3fc/mdNxloBN/qmWbAkiVyKyAFgAcNDk6fz4omPcBlRBmuuqOOnQ9BMs+njtPu/r7M/1eeknFx/Lqi3twxhZ/5wkYBH5AfApIAy8B3xBVdtcxGIKpHRXiNwMTEl5PNl/rgd/F4+FAHPmzNHzjplUmOjMwLziT0VOk4APO6iJww7KrXQwXFyNgngKmKWqRwNvA992FIcpCEVjsbStBLwEzBCRg0WkBrgI6DuVypSGy37ntSLhJAGr6pOqGvUfLsPrVZhylRgFka71o78hYCJSKyL3+a+/6G/eOXyhe9fptcATwBrgftv0s4QFq71WJIqhBnwlcF+mF1Nra3VkLqCb4qWqg+rt5jgE7Cpgj6oeJiIXAf8K9P37cghU9THgseH8TOPIq4u9r8fOdxuHL289YBF5WkRWpmnnpRxzAxAFFmf6HFVdqKpzVHVONYWbo22G1yBLEMkhYKoaBhJDwFKdB9zpf/8AMFf6nQ5lKtaKe7xWJMTV9lkicgVwNTBXVUM5vmcH0AGUy463rZTHz9IKvKmq89K9KCKP+8ekUwd0pTxObksvIhcA81T1H/zHnwdOVNVrUz57pX/MJv/xe/4xzv67isg+4C1X56c4rivXMbg+f52qzurvIFejIOYB3wI+lmvyBVDVsSKyPNs25qWkXH4W/+dIm3wBsr1Wpt5y+e9aDNeV6xiK4fy5HOdqFMRP8ZYcekpEVojIrY7iMMUtlyFgyWNEpAoYiW1Lb0qEkx6wqh7m4rym5CSHgOEl2ouAS3od8whwOfAC3hb1z9q29KZUFMMoiIFa6DqAYVQuP0tefg5VjYpIYghYEFikqqtE5LvAclV9BLgd+JWIvAvsxkvSrrn+d3V9fnAfQ0mc39lNOGOMqXS2HrAxxjhiCdgYYxwpyQQsIj8QkTdF5HUReUhEWlzHNBDlssOCiEwRkSUislpEVonIda5jKhaur1ERudD/N4mLSMGGY7m+tkVkkYhs98eHF9xAfydKMgFTwov5pEyvPQuYCVwsIjPdRjVoUeAbqjoT+BDwlRL+WYab62t0JXA+8FyhTlgk1/YdgMtx5wP6nSjJBFzii/nkMr22JKjqVlV9xf9+H95iNbYOI+6vUVVdo6qFno3n/NpW1efwRsM4MdDfiZJMwL1cCfzRdRADkG6HhZJPWv4qZMcCLzoOpRiV2jU6WGV5bQ9WLr8TRTsOWESeBtJtc3qDqj7sH9PvYj4m/0SkCXgQ+Jqq7nUdT6G4vkZzOb9xI9ffiaJNwKp6WrbX/cV8zsFbzKeUBjPntMNCqRCRarwLbbGq/tZ1PIXk+hrt7/wOlNW1PVgD+Z0oyRJEymI+5w5kMZ8iUTY7LPjLPt4OrFHVH7qOp5iU+DU6WGVzbQ/WQH8nSnImnD/ttJYDi64sU9UvOQxpQETkbOBHHJhee5PbiAZHRE4BngfeABL7zl/vL2Be0VxfoyLyGeC/gLFAG7BCVc8swHmdXtsici/wcbzlKLcB/6yqtxfw/AP6nSjJBGyMMeWgJEsQxhhTDiwBG2OMI5aAjTHGEUvAxhjjiCVgY4xxxBKwQ/54yRf9laPu88dOGlPyRORa/7pWEcm0I3bFswTs1r8C/+nvkbcHuMpxPMYMlz8DpwHrXQdSzCwBF4CITPfXhl0sImtE5AERaQQ+ATzgH3Yn8GlnQRozCBmu7QZVfVVV17mOr9hZAi6cI4CbVfVIYC9wDdCWsmRhRa8cZUpa72v7y47jKRmWgAtno6r+2f/+brzpksaUg97X9ikugyklloALp/ec7wjQIiKJFekqcuUoUxZ6X9u2vkGOLAEXzlQR+bD//SXAUmAJcIH/3OWAreFqSlG6a9vkwBJw4byFtz/UGmAUcAvwv4Cv+ytnjcFbxs6YUtPn2haRr4rIJry/7F4XkV84jbBI2WpoBeBvTfKoqs5yHYsxw8mu7aGxHrAxxjhiPWBjjHHEesDGGOOIJWBjjHHEErAxxjhiCdgYYxyxBGyMMY78f60vNTqHo6oCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 1+2 dim marginal matrix plot via kernel density estimates\n",
    "arr_ax = visualization.plot_kde_matrix_highlevel(\n",
    "    h,\n",
    "    limits=prior_bounds,\n",
    "    refval=gt_par,\n",
    "    kde=GridSearchCV(),\n",
    ")\n",
    "arr_ax[0][1].plot([-2, 2], [3, -1], linestyle=\"dotted\", color=\"grey\")"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}