Symbolic-PyMC Beta-Binomial Conjugate Example

beta-binomial-conjugate-example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Symbolic-PyMC Beta-Binomial Conjugate Example\n",
    "\n",
    "Using the example model from\n",
    "[here](https://github.com/zachwill/covid-19/blob/master/covid-19.ipynb),\n",
    "we'll show how auto-conjugation can save one from needlessly sampling a\n",
    "posterior known in closed form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "import pymc3 as pm\n",
    "\n",
    "import theano\n",
    "import theano.tensor as tt\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from operator import add\n",
    "\n",
    "from unification import var\n",
    "from etuples import etuple\n",
    "\n",
    "\n",
    "from kanren import run\n",
    "from kanren.core import eq, lall\n",
    "\n",
    "from symbolic_pymc.theano.meta import mt\n",
    "from symbolic_pymc.theano.pymc3 import model_graph, graph_model\n",
    "from symbolic_pymc.theano.utils import canonicalize\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "\n",
    "theano.config.cxx = \"\"\n",
    "theano.config.mode = \"FAST_COMPILE\"\n",
    "tt.config.compute_test_value = \"ignore\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(\n",
    "    [\n",
    "        [\"Singapore\", 72, 0],\n",
    "        [\"Italy\", 46, 21],\n",
    "        [\"Japan\", 32, 5],\n",
    "        [\"Hong Kong\", 30, 2],\n",
    "        [\"Thailand\", 28, 0],\n",
    "        [\"South Korea\", 24, 16],\n",
    "        [\"Malayasia\", 20, 0],\n",
    "        [\"Vietnam\", 16, 0],\n",
    "        [\"Germany\", 16, 0],\n",
    "        [\"Australia\", 15, 0],\n",
    "        [\"France\", 11, 2],\n",
    "        [\"UK\", 8, 0],\n",
    "        [\"USA\", 7, 0],\n",
    "        [\"Macau\", 6, 0],\n",
    "        [\"Taiwan\", 5, 1],\n",
    "        [\"UAE\", 5, 0],\n",
    "        [\"Canada\", 3, 0],\n",
    "        [\"Spain\", 2, 0],\n",
    "    ],\n",
    "    columns=[\"country\", \"recovered\", \"deaths\"],\n",
    ")\n",
    "\n",
    "df[\"combined\"] = df.recovered + df.deaths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    p = pm.Beta(\"p\", alpha=2, beta=2)\n",
    "    y = pm.Binomial(\"y\", n=df.combined.values, p=p, observed=df.deaths.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numeric Optimization\n",
    "\n",
    "The naive approach to estimating this model would involve some wasteful MCMC sampling, which we do below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [p]\n",
      "Sampling 4 chains, 0 divergences: 100%|██████████| 32000/32000 [00:48<00:00, 660.25draws/s]\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    brute_trace = pm.sample(draws=6000, tune=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Symbolic Optimization\n",
    "\n",
    "Now, we'll walk through the process of creating a rewrite rule--from scratch--that converts Beta-Binomial models into their closed-form posteriors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the PyMC3 graph into a symbolic-pymc graph\n",
    "fgraph = model_graph(model)\n",
    "\n",
    "# Convert the graph to a more consistent form \n",
    "# (this helps \"normalize\" the patterns that we want to match)\n",
    "fgraph = canonicalize(fgraph, in_place=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def betabin_conjugateo(x, y):\n",
    "    \"\"\"Replace an observed Beta-Binomial model with an unobserved posterior Beta-Binomial model.\"\"\"\n",
    "    obs_lv = var()\n",
    "\n",
    "    beta_size, beta_rng, beta_name_lv = var(), var(), var()\n",
    "    alpha_lv, beta_lv = var(), var()\n",
    "    # Match a generic Beta random variable\n",
    "    beta_rv_lv = mt.BetaRV(alpha_lv, beta_lv, size=beta_size, rng=beta_rng, name=beta_name_lv)\n",
    "\n",
    "    binom_size, binom_rng, binom_name_lv = var(), var(), var()\n",
    "    N_lv = var()\n",
    "    # Match a generic Binomial RV that uses the aforementioned Beta RV\n",
    "    binom_lv = mt.BinomialRV(N_lv, beta_rv_lv, size=binom_size, rng=binom_rng, name=binom_name_lv)\n",
    "\n",
    "    # Construct the posterior parameters from the prior parameters and observations\n",
    "    obs_sum = etuple(mt.sum, obs_lv)\n",
    "    alpha_new = etuple(mt.add, alpha_lv, obs_sum)\n",
    "    beta_new = etuple(mt.add, beta_lv, etuple(mt.sub, etuple(mt.sum, N_lv), obs_sum))\n",
    "\n",
    "    beta_post_rv_lv = etuple(\n",
    "        mt.BetaRV, alpha_new, beta_new, beta_size, beta_rng, name=etuple(add, beta_name_lv, \"_post\")\n",
    "    )\n",
    "    # Construct the posterior from the posterior parameters\n",
    "    binom_new_lv = etuple(\n",
    "        mt.BinomialRV,\n",
    "        N_lv,\n",
    "        beta_post_rv_lv,\n",
    "        binom_size,\n",
    "        binom_rng,\n",
    "        # Give it a descriptive name\n",
    "        name=etuple(add, binom_name_lv, \"_post\"),\n",
    "    )\n",
    "\n",
    "    # TODO: We could also transform non-observed conjugates.\n",
    "    return lall(eq(x, mt.observed(obs_lv, binom_lv)), eq(y, binom_new_lv))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: We could walk the entire graph and find all occurrences of Beta-Binomial\n",
    "# conjugates, but, since we're only looking for observed Beta-Binomials,\n",
    "# they'll always be the function graph outputs.\n",
    "q = var()\n",
    "res = run(1, q, betabin_conjugateo(fgraph.outputs[0], q))\n",
    "expr_graph = res[0].eval_obj\n",
    "fgraph_conj = expr_graph.reify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_conjugated = graph_model(fgraph_conj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For our conjugate model, we can just sample the posterior term directly\n",
    "conj_samples = model_conjugated.p_post.random(size=len(brute_trace[\"p\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/bwillard/apps/anaconda3/envs/symbolic-pymc/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/distplot.py:38: UserWarning: Argument backend_kwargs has not effect in matplotlib.plot_distSupplied value won't be used\n",
      "  \"Argument backend_kwargs has not effect in matplotlib.plot_dist\"\n",
      "/home/bwillard/apps/anaconda3/envs/symbolic-pymc/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/distplot.py:38: UserWarning: Argument backend_kwargs has not effect in matplotlib.plot_distSupplied value won't be used\n",
      "  \"Argument backend_kwargs has not effect in matplotlib.plot_dist\"\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7f0d137ed860>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f0d13063b38>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7f0d1386b9b0>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f0d1303f630>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEoCAYAAAAt0dJ4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3zT1f4/8Ncnu026J20po1CmUKYoS5a4ABeXofVyUdGLG0XFLSKi96f3irhFcCJ8xatecStS9i6zBQpltKV7Jm2SJjm/P9qEJP0k+STNbN/Px+M+rjTrfD75JDnvc97nfTjGGAMhhBBCCCGEkIATBboBhBBCCCGEEEJaUIBGCCGEEEIIIUGCAjRCCCGEEEIICRIUoBFCCCGEEEJIkKAAjRBCCCGEEEKCBAVohBBCCCGEEBIkKEAjxE+eeuopnDlzxqvPefjwYSxfvtyrz0kIIYR4E/3+EeIeCtAI8ZNdu3bB29sOFhQUoKKiwqvPSQghhHgT/f4R4h6ONqomncnu3bvxn//8B+Hh4SgpKUF6ejpee+01REVF8d7/rbfewpkzZ1BSUoKqqipMmTIFjz/+ODiOw7fffovVq1cDANLT0/Hiiy8iPj4emzZtwocffgiRSITo6Gj861//wrp16/Dee++ha9eu+PDDD5GWlsb7elqtFi+//DKOHj0Kg8GAadOmYcGCBVCr1XjsscdQVlYGxhhmzZqFq6++Grfccgs0Gg1mzZqFxx57zGfnjRBCSGij3z9CQggjpBPZtWsXGzBgAMvLy2OMMfbss8+yl156yeH9V65cya6++mpWV1fHGhsb2Y033sh+/fVXduLECTZ+/HhWXl7OGGPsvffeYwsXLmSMMTZx4kRWVFTEGGNs7dq1bPPmzYwxxiZMmMAKCgqctu+NN95g7777LmOMMZ1Ox+644w62efNm9t///pc9+eSTjDHGqqqq2KOPPspMJhPbuHEje/jhhz0/IYQQQjoF+v0jJHRIAh0gEuJvWVlZ6Nu3LwDglltuwdNPP+30/tdffz0iIyMBANdeey127tyJ8vJyjBs3DgkJCQCAOXPm4N133wVjDFdffTWys7MxceJETJw4EVdeeaXgtm3duhWNjY346aefAACNjY04efIkrrnmGrzxxhtYsGABxo4diyeffBIcx3ly+IQQQjop+v0jJDRQgEY6HbFYbPlvxpjNv13d3/xvk8lk8zfGGAwGAziOwxNPPIFbbrkFOTk5WLFiBcaPH49HH31UUNsYY3j55ZcxbNgwAEB1dTXCwsIQFhaGn3/+GTt27EBOTg7ef/99bNy4UdBzEkIIIQD9/hESKqhICOl0Dh06hAsXLgAAvvnmG4wbN87p/X///XdotVo0NTXhxx9/xLhx43D55Zdjy5YtlgXK69atw+WXXw6DwYDJkycjLCwM8+fPx/z583HixAkALT9sRqPR6WuNGjUKX375JUwmEzQaDW6//Xbs2LEDX3zxBV577TVMnjwZzz33HJRKJcrKygQ9JyGEEALQ7x8hoYJm0Eink5CQgBdeeAElJSXo378/nnzySaf3V6lUyM7ORn19PW666SaMHTsWAPDII49g/vz5MBqNSEtLw/LlyyGRSLBo0SIsWLAAcrkccrkczz//PABg0qRJuP/++7Fy5UpLiom9++67Dy+//DKmT5+O5uZmTJ8+HZMmTYJarcbixYtxww03QCqVYvLkyRg0aBCioqKwatUqvPjii5bXIYQQQvjQ7x8hoYGqOJJOZffu3Xj99dexYcMGQfd/6623oNPpqEIUIYSQkEa/f4SEDppBI53eo48+ioKCgjZ/79atG3r37u3X11u5cqXXX48QQgjhQ79/hAQnmkEjhBBCCCGEkCBBRUIIIYQQQgghJEhQiiMhhJAOwWQy4YUXXsCJEycgk8mwbNkydOvWzXL72rVrsWnTJgDA+PHjcf/990Or1WLx4sWoqqqCUqnEq6++itjY2EAdAiGEEEIzaIQQQjqG33//HXq9HuvXr8ejjz6KFStWWG67cOECvv/+e3z11VfYsGEDtm3bhvz8fKxbtw6ZmZn48ssvceONN+Kdd94J4BEQQgghHs6g5ebmQi6Xe/yiOp2uXY8PVh31uICOe2wd9bgAOrZQ1FGPC3B8bDqdDllZWV55jf3791vKgGdlZeHo0aOW25KTk/HRRx9ZNt41GAyQy+XYv38/7rrrLgDAuHHjBAVoBw4cQFhYmFfaTDr2dR8IdD69i86nd9H5vMTZ759HAZpcLke/fv08blBeXl67Hh+sOupxAR332DrqcQF0bKGoox4X4PjY8vLyvPYaarUaKpXK8m+xWAyDwQCJRAKpVIrY2FgwxvDaa6+hf//+6NGjB9RqNSIiIgAASqUSDQ0NLl+Hamt5F51P76Lz6V10Pr2LzuclzgJVWoNGCCGkQ1CpVNBoNJZ/m0wmSCSXfuZ0Oh2eeuopKJVKy8a21o/RaDSIjIx0+ToikajDBtKB0JEHJgKBzqd30fn0LjqflzgboKQ1aIQQQjqEoUOHIicnB0BLKn5mZqblNsYYFi5ciD59+mDp0qWWVMehQ4diy5YtAICcnBwMGzbM/w0nhBBCrNAMGun0jCaGwko18i42oKimCRUNOjQ1GwEAEQoJkiIV6JmgxKDUKMSpKG+akGA1ZcoUbN++HbNnzwZjDMuXL8eaNWuQnp4Ok8mEPXv2QK/XY+vWrQCARYsWYc6cOXjiiScwZ84cSKVSvP766wE+CuKJbacqkZmsQmKEwqvPqzMYwRigkIq9+rzuMpoYKtU6JEV69/gIIcGJAjTSKVWpddh05CJ+2F+K41+dB8cB/bpEIj02HEmRcsQqpQCABq0Bhy7U4psDRThZ1oCe8SpMHZCE6Vmp6JWocvEqhBB/EolEWLp0qc3fMjIyLP995MgR3setXLnSp+0KdYwxMAaIRFygm+JQlUaH4hqx1wO0zfnlaDYyTBuc4tXndVdhpQbHSuowIys1oO0ghPgHBWik09A2G/Hb8TJ8e7AYWwsqMapnHEakhuOVWSPQM14JjnPe+VDrDNhTWIUfDl/EjW9vx+CuUbh7bE+Mz0xw+VhCCAlV+87VoKxeixsGBTZICQSdwRToJrSiwgqEdCYUoJEOr0lvxJd7zuP9LacRr5LjlmFpeOWWy5AYoUBeXh4yEoTNhKnkEkzsm4SJfZPQoG3G+r0X8OTGI0iPC8fT1/XD4K7RPj4SQgjxvxqNHkYTBQiEEOIvFKCRDkvbbMRnO8/h/Zwz6BEfjjf+loXRveK8MtsVoZDirrE9cdvl3bBmRyFu/2g3bhqaisVT+yBCIfVC6wkhxL/qmppRo9FDJhEhKVIBsQcpjTUaPaob9YIHvgghhLRFVRxJh7T5RDmm/icHPx8rxco5WdhwzxUY0zve66mIYTIxFl7VC788Mg4ltVpct3IrDl2o9eprEEKIPxwrqcOholrsPVuN89WNHj5HPY4W13m5ZaGrokEX6CYEDW2zEcdL6gPdjKB2sa6J9gkjAChAIx2MRmfAE18fxmMbDuGhSb3x9b1X4MoM7wdm9lKiw/DhHcNwz7gM3PbRbqzeVkhfsoQQ0onVNuqx43RlwF6fMRZUqalFNY04Ve56I/jObE9hNQX1BAAFaKQDyb1Qi+tXbkV5gxY/PzwONw9N82vxDo7jcPuoblh/zyh8vuscFny2H3VNzX57fUKI/+w7Ww2DMVgKSJBgFOgxumMl9fjhcElgGxEge89WwxREwak7QrPVxNsoQCMhjzGG97acxu0f7cadY3rg43kjkBARuP3KBqRE4X8PjIFUzGHmeztQXNsUsLYQQnyjuLYJap2B9za9wYRGPf9t9rTNRnyXW+zNpnlFMGcA/Jlfhj2F1W3+HrwtDox6qwFCxhh0BmMAW+M/jDGU1DZZ9jMNNUH80SN+RAEaCWnNRhOe3HgEX+4+j/8uvBLZV3QPipL3KrkEq+YMxYS+ibjp7e20JoOQTmTXmSr8drxM0H2b9KHZiQykBq0B5Q3aQDcjpJwsU+Pno6WBbkaHxRjDd7nF0AfNtgzCBVMaLLmEAjQSshq0zZi/di9OljfgvwuvRO+kiEA3yYZIxGHJtf3wwKTemPPhLmw+UR7oJhFCvMjRSLc2REfuvaFK45/1Mx1xlsGX143QGd1A0TYbg3bW1p3Zx2YHac/HS+px4HyNN5vltrqmZmh4Zv1/OFyCusa2yzFqG/XIdVH0jDHWqb/vfIkCNBKSyuq1mPneTqjkEqy7exTiVIFLaXQle1Q3/GdWFh788mCnXQ9ACCG+Vq3Wh+y6wIoGHX455uUZrsAnk9hx3KBfjpXidIUa+85Wo1qj92ObXCuqb2737GNBhRoXrCqjMsZQ19SMoprGNmvVmcBkXZOJuVVQ5K8T5cg5WcF7m5YnAC2s1OBclcbpcxaUq71/3RIAFKCREHSxrgmz3t+JUT3j8PbcoVBIxYFukkuT+iXhw78Px5JvjuDbg8G33oQQ4nuMMRTV2Jav90ZGdnmDFgcDPDrvTLVG75fZEY3egIIKtc9fxxf0LgLL9p6+4JybstWkN6G4tskmkAkGWkP7z5799X++uhF/nSjH/nM1LmepHDlf3eh2lVCjlz+HjSGYou1o7XCwoQCNhJSyei3mfLALk/ol4flp/SHyYCPVQBnVMw5r5o3Ac98dxdf7iwLdHEKIjzhaB1uh1mH/Oe8HUoUVGo/3LXPEm924racqvN4+R0wBnkAL1jS9YPJdbjFOlfmv3H5714b54i01eGHdF11pnvkjryzoZmn5UIBGQka9thl//3gPxmUm4Jnr+wVFMRB3De8ei0/mj8RLPxzH+r3nA90cQgixcJQu5Wp2R4hgqENQUN6AKrXv1sjVNTXj+0OdO429QduMYyUtRbHqmpodzlZUqtt2kO1/0r1ZvMLVGrIGbTOMJtbu0vwMQEltk8vncRb0+TLG78jjB0eL6wSvhwuFVGgK0EhI0BmMWPDpPmQkqPD8tAEhGZyZDUmPwWd3jsTLm/I6/Y85IR1ZcW0TvsstxoXqRp+W0vd18OONjrK7M0s1Gr3TSpgmD3qax0rqkXfRdzM35s6hq2Pl67z7e+ZNZzBi79m2WxW019nKRhSUt6SZ/nWiHH95WBzrdIXar2u2/8wvxw+HS/Crk2uOr8AGn71nqzvE9jp8/awmvREFVpuNW1+1R4rqcKI0cBuRn65Qo7TOeXXXUKqySQEaCXpGE8Oi9YcAAK//bTDEIZTW6MigtGh8eMdwPP3NEfyZL6wcNyEktJg7A0LS+5r0RsGV4gxGk02HviOWnC9v0HlcebC0Tov80novt6gtxhhqG9vOBJXVO56la9Qb8D8ngYcn1QyLahrxPzcH+6o1epTwBBGNeoNX1zN6Gnc2aN17748W17UJfJuNJsvAiNBmmD+Df+SVobze9nPl6VoxfzCaWLvS9sz7MbqagSqs1OBYCf9n60yl2i+fu/Y4Uxk6a1QpQCNB76UfjuN0hRof3DE8JAqCCHV5zzisnDsED63Lxc7TVYFuDiFEoEZ9S+AltNNX6SStjmutbPfr8VJsO8W/4J8xhqPFdZaO+6YjF3G4yLd7KwY6Fcq+mIo7jl90PpLvrdmq0nottlhVxTM/rbkQA1+ih6sR/F+OlaKwUoOyeq3LAhD5pfXIL61HpVrPO6PIeVDGsaRW2671gkIrEHrb6Qo1au2qIbZntkStM6C03j8DH/bXoydn8EyFGltP8Vdo5Hvug+drbPod5uCzvqltuf3Ojm8gwx8oQCNB7as95/HT0YtY+4+RiFRIA90cr5vQJxErbhmEBZ/tw+Gi4B2dI4Rc4u1KaGZ863V2nq6CRm/E6Qq1TWGBeq37HamckxXYZ5fWVl6vFTxz50sF5Q02MwD250Lo2iBX2e+1TXpBqeX26Wx8KaoGo2+ug6ZmIy5UN7osoX6itCGgKWVHiuo8XtPni1UKewq9n7LpLnOw5ejzaTCaPPrsuuLqO8ncLvP/n69uDJmZ9yZfVIp0cv1Zz4ybTAx7z1b7pg0uUIBGgtaB8zVYtikP790+DMlRikA3x2euH9QFS67th/lr9wVdeWFCSGCVN2jRYNWha8+msDWNepS0rtHIvVALk4lh55kqhylL/tSyPsxxO/53uKTNflGesF5Pp9YZHAanv+eVoczNGRRHMy58HfJmowlHiy/NgrYn5ncV6/BtQuzo9dwJ1s9UqnGm0vk+WZ7QG00O99+qa2xuE6ybgw5nbffk/LozA2k/oFBQ3rIeyv69P36xHpvzbdfl2a/1sp5RM5oYvsstRl1Ts80gAd976oy7x2/9PWNdUMPdwHr3mSrBFTt/OFzSJk2zpLYJvx4vRbPR1O4CLoCw83ChuskyMx7IRAIK0EhQKq/X4t7P9uO5G/pjSHpMoJvjc3MvT8esEWmYt2YP75oGQggxmphXNoU1mRjOVWksm9P6O52RMWYpJOEOnV1w6ihVUejx/JFX5jS93N1CD46CCvsOOdCyxu60n/Zs++uk8EIdzt4XvcHk1WI3JhPjLbzRwBPQ7imsxpaTFfjrZDlO2bXR4Xn0cJbOk0GQP3ne492FVW1mrJtbgx2D1X4QzlJuzfezH6T562S5TdaN9VPwDWQI/YgzAFXqlk3Tja2vfchJOrV9283BpkZnwP5z1Sit16KwNYjPvVDrdMbVaGJtbte1pqn+eOQiDvhpr0dfZUi4iwI0EnT0BhP++cUBTB2QjL+N6Bro5vjNY1f3wcDUKCz4bH9QpBwRQnyvyY3OIF+/wXqE3tO1EiWt1SbtXs2j5xKiXmuwlGLn4yidyD54qNLoca5K02atkZCNaM3ru5rt0hSb9EbLubDPaPj+UInDIKw9mvRG1LoxI2Iwmngra/5wuMRl2lpJbRMKytUu0zPLG7RtZiz43hdXfVlnMy6nK9T4PU9YkayLdU2WwUuD3WZ31gVFbGZgPLyEvTEI4opG5/lvvP16Q8YYTlrNUv11otzpuldXzNtqFNW0fJ/YD4wAl84z33YJQMv7ZX681mBCUU0jzlVpeGdcyxsupVkzAFtOVlgCP+vLh69wTHsyCrypXtssuN/GGBP0XU0BGgk6y3/Mg4gDnr2hf6Cb4lccx+G1WwdBxAGP/d9hr0znE0KCl8nEsLvQdYGgKp5OULVGj/xS23QpIWXTGWPQ2FVH5Csw0d5BZG2zEX8I7HybNemN+OVYKXbZnRPzLEsFT6cz90ItzloFTfYdnyq1Dhqdoc16Lvt/7z1bje9yi12kyTGbjq918GE98i+kCIl1x/PX46WCq8vpDSb8kV+OnJNtC0IYTQyVDc4zMPaercaxkjocKXZeZGbn6SpLsZCKBp3NcXuj8jCHS7Mj3mRTKMP6JLOWz4x9ZUZnGnS2QbN3JlZcT+vxvc7RYsfpv3z390X/oVFvsASC5vMspNIqYwz7zzme/dp5uspyfHqDCbWNejQ1t1S1PWQ1S8iXKvzLsVIUlDdg/7m2331Hi+sc7nem0RndylZyldq5Ob8cewtrLMdgnjnVNhvbzBJXqHWCvqspQCNB5a8T5fg2txhvzRkKmaTzXZ5yiRjvZw/H8ZI6vPVnQaCbQwhx4lyVhjdlx9FvuafVWs0pXPYpd3zpaIwx5F6odRok8KVj2TvuZD2YEOZNivkKXTjq7Gj0BmibjW06l9azLNbBAd+sk/1s2raCSmw9VeGyIqI5sPN0z7dtBc6f3x3O1j6pdS3nqL2FJuxnofiYz8SO05XYbnV8QkrgO6pIyqe9xSoc7Ytq/xHYVlCJnWeqbNb+mUwMPx25yPv4iobW4F4jvOS//dXToDXwpmy6yzxwcLGuyaPqnI7Yz7SdLlc7fH/Vbm194LyNBeVqp2tKfztexjvLazIxy0yueZ/EYyX1ltm6SrUOx1vX1J6uUKPaQRB2uKjWpgKrI+7st2j+TNU06i3f1VtOVrSdJRb4lJ2vB0yCVrVGj8VfH8bymy7r0EVBXIkKk+Kjv4/Amh2FDn84CCGBd766EdsKKtuk1An5/S0ob8A5N4sCCdnnyNi6vsw8O1Gt0eN8VSN2eDGAsPfz0Ys2KVbWXAVGrtinGVp3Hs1pktZ9KL5z785MjSdr4+yDOmfvk5D+nlpncBgo8pVSP8uTdlmvbbYJhB093x95ZZY0MfvCEybGbP/mRlxQpdHZvK6zAYMD55xXMHaVxmb/3OY1U7vOXBoQYWCW+5kHPPafq0aVRm9J6eNjYqxNGqwzuua2z+WqiD7fqWnQthQFsQ+W8l1U7bQeMLIOgKwL3tRYXZ/b7b4XKtQ6wZVBPemfmI/1WEnLVhjO0pHNQZe1vNJ6/Hq8FI16A+/sXXFtE06VO25/jYBZM+tBnvLWASbGgH1nqy1rCPk0G1uCfevrztn9XZF4/EhCvIgxhiXfHMbY3vG47rIugW5OwPWIV+KtOUOw8PMDSI8Lx4CUqEA3iRDiwMXaJvROioC22QiFVCwo7cdR5cTcC7VIjlR4PEhl7uuV1DZBLhW3KVJgrdnguuNp3VljjPHOVugMJlQ06CB3kvVgfqWimkYcPC98SxFfFwawX7/jrLNu72JrRUz76pOFlRrEq+RO05jMI/N8qWgX6y51TNU6A44V1yFWKRPcLqBtYRJH8ZFaZ8Avx0oxbVBKm2IiZys1NrNNx62uWfM6PfuU0EKeNUbVGj22nqpw+DtmP5tnPUNUVq9tM9PiKsjlWydoP1uiN5hQVNMEkYvcNUcvVdvYEtgZTQxdosIsf3cWHACXgg7r1Dv72VDra75RwFo168dbF7j460Q5JvVLgkousQlMXO1zJ3TWyNlnxdkaU7OLdU0Ic7K3LV8VVXOwbp49c5f9jD7fdxrf7HRhpQbFtU2WwZ4rM+IsjzN/zvm+93lPpcCBDgrQSFD4v/1FOFZSj58eGhvopgSNsb0TsOjqTCz4dD++u3804lXyQDeJEMLj+MV6SMQiHC6qxYysVJf3d5byZE6b5AvQhGwia56xcLXGCBC2qbB1+tn3h0raHJ/1DIt5s1ug7T5iR4vrIOY45Jc2uOwAtmetj/WItbOy/dYqrTptrmYp661mNMxpkXwdVb4ZxXpts2Utjjkm4ZslsHaxtgml9Vqb13XG0/3F9vMEwvazG3ypiEJmdfkqNVpzllZa19TcJtXtdIUaPeKVUMpburCebKy9r3XNkqvHbrdK1azS6FBc24TU6DCbgM86yHC1d52Z0FRaIWsTna1P8zRl116Dthk7zzhO0bYU+XDzw2s9GAHYVuTkq6bYnsPhGwyp1xrQbDRBIRVDJZdAZzBagkCTiVnSOs2DBOaUUIOJQSpu+bAJKfhhNDGIRbYfTlefC0pxJAF3vqoRL/1wHG/8LQsRHXAz6vaYd2V3jOkVj39+vr/N2gpCSPDgW3B+3MEs2bkqz/Y7tC/u4WuFlRreDvh3ucX4LrcYjDHLrIv9Wha+APFQUa3NjIs5gKp0Y40PH0eBpqO0S3uNbmxCW9/UjKIaz96/4pqmNu/hwQvCZgiFzMo64yoo9rQCqBDmIME6VdHRnnF8+Aq3aJuN7doXT2jaon3wzTcj7aoSa3sGHOyPkYOwMvC7Wzft9tYGy/ZbG9j79Zj12lDHQbvBaIJG37YvI3Rgwd3r1Ho21tFrbC+oxO7WtYmb8yss33m7C6uRX1rP+9iaRr3L9EWdwWj53Kl1hjbvpasKphSgkYAyGE14ZEMuskd1w8gesYFuTtDhOA4v3TgQAPDst0fdHp0ihPgH3ydTSCEGPnqDCReqG/HXCeH7V3nCnKLnyFkXmxBXCZg9cca8vqy8nQGaq8qFzjQ1G90OfJ1VpDPzxqbagPuV4h0FC+5uuu0OV9samCvx2e9XZn+O+J7HUSf4dIXG558PR9zdBoevnXwVSYWoa2puUx2VL/g2/213YZXTvceE0BlNbdaCOno9oG06qbUKtQ57i9o+l5DsAHeYB4yKahqdVrTc3TorqNYZcLpCbfPeWs8Y28d2LZUn65wGjNZFoQ6er3H7eqUAjQTUu3+dhs5gxMOTMwPdlKAlk4jw7u3DsK2gEmt3nA10cwghTpS6CHoA1wvH9UYTzlZpXHby+dKXhLy+mdAZJkfsCwx4wmhiKFe7X2TFmrkgRaAIWSfkL6460r5gHjhs755U5iDSeg0QY/xVGu1T4/zp56O+3yfNHa4KX7hT5ISPt4Mna+YBBb7vOmM72m3+Hjxf3eh0Lag7+1Dac/Xc1sdk/u+TZQ2Cq/nSGjQSMIcu1OL9nDP49r4rO2VJfXfEq+T44I5hmP3+LvRNjsQVGXGBbhIhQcdkMuGFF17AiRMnIJPJsGzZMnTr1s3mPtXV1ZgzZw6+//57yOVyMMYwbtw4dO/eHQCQlZWFRx991O3XNg8gC9kgVsiaGSFre/hY7xsUCtwpYx1MrEfahQSIvx0vs/mdOy1w37NQYDK1LFUQmq5p1p7NlNvLnT2wOhpXa5+CRXtSuq1na83ptO35qnG0lYO73JnJpgCNBESj3oBH1udi8dQ+6JUYEejmhIQBKVFYdtNA3P/lAXz/wBikRoe5fhAhncjvv/8OvV6P9evXIzc3FytWrMC7775ruX3r1q14/fXXUVFxKQXn/PnzGDBgAN577z2vtMHd9KfOjIG/0xQKqdzuzqK0lAW/9G+hsxJCC51Ya++MibtOlTcgMsx368edzSTT2mxbjmZnrNdieZpeGcrqtc0orHCesu2Ml+IztwbeaNqCBMTLm/LQNTYcd1zRzfWdicWMrFTcNCQV9362v93pJIR0NPv378fYsS2VYLOysnD06FGb20UiEdasWYPo6GjL344dO4aysjJkZ2fj7rvvxpkzZzx6baHV20KJwUsV4BzRNht59/Ai7SOkiIQ3NeqNHgWSQpyuUDud7eIr6x+Kfjnm27TJ0vrApYQGg8355SH3XUMzaMTv/swvw49HLuKXh8d5bdq4M3ny2r644+M9eOq/R/D6zMF0DglppVaroVKpLP8Wi8UwGAyQSFp+6kaPHt3mMQkJCViwYAGuvfZa7Nu3D4sXL8bGjRudvg5jDCXFJby3GeokbdZUdSR7UIOSYu919kqKS9DcrLc5n3IJB52APdoIv6hYDiV+XItWLeGgDaECAWUAACAASURBVND7JWusREmNb9MV7a9PX8kT1aCk2Pupr3vVFSiq8906Mnf563x6k75WjEqN9wfFM3uqHN5GARrxq0q1Do9/fRiv3HwZEiM924i1s5OIRVg1dyimvbUNn+w4i3mjewS6SYQEBZVKBY3m0iipyWSyBGeODBw4EGJxyz5Gw4cPR3l5ucMNmc04jkNKahfe29JiwiBxsbdVKEtMi0YKvLvOraS4BCmpKV59zs4sz8/nk+O4gKWl9kyOhL7UN7N3Zv66PrtnJCNF5/2ZNBOAFMdxgN+F6uc9Jdr1fdzn+NqlFEfiN4wxPLnxCCb2TcQ1A/k7N0SYWKUM72cPw//79SR2Odk8kpDOZOjQocjJyQEA5ObmIjPTdXXYVatW4ZNPPgEA5Ofno0uXLjQrTYgbArlmUMhm66FidyH9lpNLaAaN+M1Xey/gZFkDfnxobKCb0iEMTI3Cshtbi4bcPwYpVDSEdHJTpkzB9u3bMXv2bDDGsHz5cqxZswbp6emYNGkS72MWLFiAxYsXY8uWLRCLxXjllVfa1YaiDjx7BgDFPtzQmBB3nfdw0/dg5K2980jHQAEa8YvCSg2Wb8rD2vkjoJLTZectNw5JxZHiOtz7+X5suOcKKKTiQDeJkIARiURYunSpzd8yMjLa3O/PP/+0/HdUVBQ++OADn7etowhkaXRC7LVnHytCghmlOBKfMxhNeGR9LuaN7o5h3WID3ZwOZ8m1faGUSfDMt0dDojw1IYQQQghxjAI04nOrNhfAxBgenNQ70E3pkFqKhgzBztNV+HTnuUA3hxBCCCGEtAMFaMSnDp6vwUdbC/HvWVmQiuly85U4lRzvZw/Dv345gd1UNIQQQgghJGRRj5n4jEZnwCPrc/HktX2RkRBENV47qIGpUXjpxgG478sDKKGF/IQQQgghIYkCNOIzyzYdR88EFW67PD3QTek0bhqShumDU/HPz/dDbzQFujmEEEIIIcRNFKARn/jteBl+PVaGV28ZRHsK+dmS6/oiTCbGql2VVDSEEEIIISTEUIBGvK68QYsnNx7GilsGISFCHujmdDpSsQhvzx2K3ItN+HwXFQ0hhBBCCAklFKARr2KM4YmvD+PqAUmY0j8p0M3ptOJUcjw7IRmv/nwCewqrA90cQgghhBAiEAVoxKu+2H0ehZUaPHN9/0A3pdPrHSfH0hkDsPCLA7hYR0VDCCGEEEJCAQVoxGtOV6ix4qd8vDErC0q5JNDNIQBuHpqGGwZ1wb2fH4C22Rjo5hBCCCGEEBcoQCNe0Ww04ZH1ubhzTA8MTY8JdHOIlaev7weFRITnvjtKRUMIIYQQQoIcBWjEK/7920mIRRzun9gr0E0hdqRiEd6+bSi2narE57vPB7o5hBBCCCHECQrQSLvtOF2Jz3adw5uzhkAqpksqGMWr5Hgvexhe/Skfe89S0RBCCCGEkGBFvWnSLjUaPRatP4SlMwYgPS480M0hTgxKi8aL0wfgn58fQGmdNtDNIYQQQgghPChAIx5jjGHJN0cwqmcsbhqSFujmEAFuGWYuGrIfOgMVDSGEEEIICTYUoBGPfbX3Ao5drMPSGwcGuinEDU9f3w8yiQjPfXuMioYQQgghhAQZCtCIRwrK1Vi+KQ//mTUEkQppoJtD3CAVi/DObUORc6oCX1DREEIIIYSQoEIBGnGbzmDEQ18dxN3jemJYNyqpH4riVXK8nz0MK37Kxz4qGkIIIYQQEjQoQCNue/WnEwiXiXHfBCqpH8oGpUXj+Wn98c8vDqCsnoqGEEIIIYQEAwrQiFt+PnoR3+YW483ZQyAWcYFuDmmnmcO74rqByVjw6T406aloCCGEEEJIoFGARgQ7V6XB418fxht/G4yU6LBAN4d4ybM39Ed0uAwPfXUQRhMVDSGEEEIICSQK0IggeqMJC784gL9f2R1X9UkMdHOIF0nEIrx921BcqGnCy5vyAt0cQgghhJBOjQI0Isj7e6oQqZDi4cmZgW4K8QGVXII180bgxyMXsWZ7YaCbQwghhBDSaVGARlz69mAxdl5oxJtzsmjdWQeWHKXAmn+MwBu/ncSvx0oD3RxCCCGEkE6JAjTi1KmyBjz77VE8OS4RiRGKQDeH+Fi/LpFYNXcoHt1wCIcu1Aa6OYQQQgghnQ4FaMQhtc6AhV8cwL1XZWBQMhUF6SzGZybgqev74c5P9uFclSbQzSGEEEII6VQoQCO8TCaGRzfkolucEv8cnxHo5hA/mzMyHXMvT8ftq3ejtI72SCOEEEII8RcK0AivlX+ewpkKDf49azBEtO6sU3pkcm9M6puE7NW7Ua3RB7o5hBBCCCGdAgVopI2fj5Zi7Y6z+PCO4YhQSAPdHBIgHMfhuRv6Y3DXaMxbswcN2uZAN4kQp0wmE5577jnMmjUL2dnZOHfuXJv7VFdXY+rUqdDpdAAArVaLBx54AHPnzsXdd9+N6upqfzebEEIIsUEBGrGRX1qPxV8fwsrZQ9A9Xhno5pAAE4k4rLj5MqREheHOT/ZBozMEukmEOPT7779Dr9dj/fr1ePTRR7FixQqb27du3Yr58+ejoqLC8rd169YhMzMTX375JW688Ua88847/m42IYQQYoMCNGJR0aDDXZ/sw4MTe2NcZkKgm0OChEQswptzshCpkOAfa/ZCTUEaCVL79+/H2LFjAQBZWVk4evSoze0ikQhr1qxBdHQ072PGjRuHnTt3+q/BdqLCOl7GQnIkVf8lpCOQSyhk8CdJoBtAgkOj3oA7P9mLcZkJuGtsj0A3hwQZuUSMd24bhvu/PIC/f7wHa/8xgtJfSdBRq9VQqVSWf4vFYhgMBkgkLT91o0eP5n1MREQEAECpVKKhocHl6zDGUFJc4qVWXxKfrEBJaccqyiPVSFFS6zw9urlZ75Pz2VnR+fQuOp8tFBIOWgNr9/PQ+bwks6fK4W0UoBEYTQwPrjuIOKUMS6cPAMdRURDSlkwiwtu3DcWD6w4ie/UefDJ/ZIcc8SehS6VSQaO5tDWEyWSyBGdCHqPRaBAZGenydTiOQ0pql3a1NSpMCr3BhKZmo+VvvXvGoVJcBQBIiQ5DSW1Tu14jGAzuHovms87X9ZUUlyAlNcVPLQq8GVmp+C632GfP39nOZ3slRMhR0aBzeDudzxZKmQQaffszaPx9PtNjw3G+utFvr+eeeoe30HxlJ8cYwwvfH0NpvRar5g6FREyXBHFMKhZh5ZwhSIsJw+wPdqG8vmON9pPQNnToUOTk5AAAcnNzkZmZKegxW7ZsAQDk5ORg2LBhPm2j2YjusW1moa2/f2Ve/i5OiwmNvSyzuka7vpMHJCLH57MHrbf2GqWsfeP+colY0P26x3nvPbsyI95rzxVMhJ7Ljk4WoqmZodlq4jXvbjmNP/PL8fHfR0AppwlV4ppULMKbs4dgRPcY3PTODpyuUAe6SYQAAKZMmQKZTIbZs2fjlVdewZIlS7BmzRr88ccfDh8zZ84cnDp1CnPmzMH69etx//33u/WaQ9Nj2ttsC7FV9kJGouPUF0/0Sozw6vP5Src4JUQ+yOJQSB13d5iTrK3U6NAIbINFe9+6CIWwfkjfLp5fz72dfBbCpMEZ1ChlEqeDDPw8T0cUizgkqOQ2fwvG5Kq+ya4zHoLR8O6xmJGV6vQ+1CPvxD7ZcRYfbyvEVwuuQCIt5CZuEIs4vDh9AJIiFbj13R346O8jMKyb9zqqhHhCJBJh6dKlNn/LyMhoc78///zT8t9hYWFYuXKl268Vq5ShWqMHxwEDUiJxrMRxqopQUeGXZtQ82X5yfGYCahubcaioFgAwuV8Sfs8rAwBIRByGpsfgwPmadrfTHZ706aYNTnGYAnhVn0Rom43YdaaqfQ0TKCpMiuIAppomRSpQJiBTQcRxMDmJNF2lzJqLuZQKeK3ECAXKG3yTPTGqZxx+OOx6fZL7wcolvRJVOFXOv9Y0OUqBwkoN722BxHFAvEom6P3xhhsGpWD/ueDfciQpUo78Uuf38cWAT3vFKWUu70MzaJ3U/+27gH//fhKfzr8cvbw8Uks6B47jcN+EXnj6+v6Y9/Ee/HjkYqCbRIjfjO19qdKtubMYE+76R9dXhqbHIDpcBoXVDIBcIrJUXpOIOUTyFPbxZuqj9Tkxi2xdp5oeG+6V14gKk0IsMHq1Pl5ns2TO5hn8lfYv4jib986cdif099nV/VzNBIbLJJAKTAXrFued95KP0PfWV7jWIYXocBkSIzwfuL6iZxz/87cjWBiYGuXWQKijFEeh57hXQgR6JaowKK0l7dgbcU73OCVGd/NvSnHvRBWGd491eb9Yu6DJfNzuiFPKXd8JsPmsO0IBWie06fBFLP3hONbMG4H+KaE5PUyCx63D0vBe9jA8/d8j+Ncv+TCa2l/liZBQwVhLh3Vc7wSMctApc/hYJ6GBs4CCT9fWACgpsqWDcHmPOEjEImQktHTc5RIxosKlbUply8RtOwoZCSr0jLft8Nuv+Zk6IFlQu1RyCWZkpWJIeozD9DGF1Pud8v5dIjGhb6Ll3550Lq/qk2g5n/6WECG3SYFytbarb3KEw6BAiGir2VshwfTEvom8AYD9ebaevRjZIxYjrDrKvgzGpg1KQb8ukUgSmB1kvt7Nn8mkSDkyEpwHEnwDHmaO1j2N6uk6UHBEKZcgLcbxe2N/jVyWFsV7P/vURUeiwqUYkBJlWaNpPwB1Wart8/fv4ro/2a9LJKRi4e97cqQC0walYELfRAzp2jY4FbI0RyIWOTzmhIiWv/O9l87Wpvpj3SoFaJ3MT0cu4smNh/FB9nAM8eLaCdK5je4Vj+/uG4M/8spx96f7UK91XlabkFBmnzLDcRxilDKHnTK+jqyv1rqYR+jD5WKbf5vZp4Yp5fztsO7c9YxXtemMCRkBFmJGVipGdeXv7HSLU2J8ZkKbdX5CgteuToKMASlRmD44xdI5Y61PeO3AS5U5r7+sC6LCpJYZFVfcnRmxv4ZcHZK5rc5eX0jAYw7Y7VmfL0fVeaNbO+gcB0QopIKC3uHdL713XaLCkGI1kzesWwzvzJ+zQPOqPokOb7MmEnHITIpAvIp/Vtt+/U9mssrtVLixvePd+hwrZRK3ZuXiBQZSfZIjMDA1CpP7J3n0eKHEIg7drAZqpB7MLjs7xXzncmSPWIhELbP/cp51pK7aYH5ORx8N87Xm7PvCFevvd+vXiVXK2jVjSgFaJ/LNgSI8vvEw3r9jGK7I8HykjRA+6XHh+GbhlQiTijFj1XYcK6kLdJMI8Qml7NJPp5BOsf2IsYjjIBLwuPasszE/e494JUb3ulSlboyATqV9n+KytCje9kaHy2yCJ+vZucn9kgSlFZnZpxNNG5SCwWlRlqAAaDtr56zYgzPd4sJtOk7mUXjrANuc2sg32m9fmOCqzETEW6VHCUkbnTZYWJlxIcGo+dyZZzhkEsfXlrPr1fxaYhFnSU01m5GVigFOMm7CW2dv7AcjukTxn4uoMCm6RIXZvL9mzqoPSuzaz/d4T/o3colY8HtiaYtY5FaQEu5gMIRPcqQCw7rFWN5bR5/ZGVmp6Jsc6TDwtherlCE1Jszj2UvrtVOuBg3cJZWIMCMr1Wa2S0iAY75WzdegeQZYLhFZPqsSschmcKDtc0gs13+CSo4bBgm7FmZkpeKqPpdSu4d2i7G55q8daPud5U7ARgFaJ/HF7nN48X/HsfYfIztsSVkSeOEyCVbNHYLbLk/HrPd3Yc32QsvoNCEdzdjeCegSZTsibh1UmTsH9p24weZgpPWj4ajjYD8j58msm1jE2YykK6RiqFxUyouQC9vfcHxmgs3Is3W6kVIucbruyT5FqEe80mZWQyTi2nRm7GftrNPy7Jk78ubnNN83Mkza5v3olaDiXT8HtHTs7Gfw7NNEo8KliG/trE4dkGzpKCa7UXzL1doivtkDM/O5tA6ix/VOsPmtt2+zM11jwnnvH6+SIyFCbgkEI6ze72HdYjCpXxL6JAsLmp11VJ2l/tpTymyviRlZqYJnqfp3iXSYzil05tRXaZpR4VIopGL0iFdi6oBkjOjheWqktbG9E5ymSTrSI16JbnFKm8+7/efR/O/xmfyfJaEu7xmHgan86Zl8JrTOqprfCfPSnYGpUUi3Wi9p/v7kq/wYLpNYrrqRPWJdvq/W3Rrr86CQih1+L/VJjoDKjSCdArRO4MOcM3jj15P44q7LqdIe8TmO43DX2J74asEofLrzHO78ZB+q1I43ASUkVPGlsFh3LrvFheOagcmWEV5zx9Z+YkzBM2MgFYswtne8TYc3pnX0OiNB5XINmNCRWnNrr2kd6c1IUNl0anylPXs0RVmdz2k8I90zslJtinuY1yMB/PvLiURcmwIBzvCd2sykCMzISoVCKrYEN9ZbJVx3mfONzVOjw5DoZEYi083ZwhilzGaGI04lxzUDk5EYIbcEkI6IRByGd2sJCOzXcF2ZEW/pkPa0mrURcS3rDX1ZMW9SvyRLcG1+lX4C1j1ZP95a76QIt5d6mAN+VWtwOlJA4DSYp9iEO0VzFFKxWzN1KrmkTZqoo7RVoQalRbt8jq6xLd93fLOaZtaBTUKEnDdVWizikJGgclqG3tn6P/NrOBsbvu6yLpjcL8nyPSnmOMtgsvV3p3Vq95UZ8Zbg0x/FbChA68CMppZNqD/eXoh1C0a5NSJBSHsNTI3CDw+MQaxShqv/nYNNh6nKI+n4THYjq3KJGH2SWjrX5pvMHQBHC9wjw6SQSUQts10894lRytq1Boyvs2cOmJx1PITuUeVr1m0UiThLMZBh3WJ4K2mKRJzVWhT3O1Z8DxnePRbjeifwBgjmWQ9r5nNuPZt5ZUY8hnSNsXREh6THOJxNFZIS64pcIkacSo4prWuVzIEa36bP5tlbvuuPj8QucGoPvo51cqQCKrkEstY0OJGIw9QByS6LRKhaZ4PjVfI26ZHOXrtngtLyee0d3zZwNl+D5s+hswGR7vFKJEcqbM6zfbAzvHssbyDnDF/RDKDl+8N+66Se8SqnKX5Ay3XQv4vwdEk+zgZehnWLsckK4MAJLjYE2Ka3usoCcPo8Yg5Sschy7Uzql4QwGX+7e9qdC/NjeiYoMcYqddz6nHV1MkOZGKEQPJtNAVoHpdEZsODTfdh7thrf3jcamUmhsUkp6ViUcgn+38zBePWWQVj6wzH88/P9Pts/h5BgwDdDY/3jf0VGHFJa0yIvS42yzFy5o72bJw9Oi7Z0LtwJWK7K5C/QYF9AxBFv7bcp4lo63OaOTqRCihlZqUiLCcc4B+lV5g60qz46X7W3lKgwDE6LtgQ2EQopUqPDEKOU8f62cnYl8x2+VoTc5WxlvEpmU/lQCKHvqHnGYHDXlsDA3QA8OlyK6HAZrhmYLDiQM7Nvo6sZpct5CoeYz3G3OCXEIg7TedaQJUcpMH1wCkb3inerkqdULLIEa6mRrmef7FMt7V3eM84SII3PTOANiq3xzaqbmQcknG2+bi1MKkZ6XDjvdWR9ziRiDr2TItwazHe0xpCPJ6mV1uJUckzul4RrBibznj+h729Pu/Rq87XraMJtcuvMq3VBJRHHIc7qu8I66HIU7AEtA9fXDHQ+m255DUH3IiGlrF6Lv72/ExwHbLjnCsFlZgnxlcn9k/DrI+MRqZDi6n/n4Ks952GicvykA7qqTwKyujoeCU+MUFwKFkRcu1L9nHE2WyCTiBCnkuO6y7o4vd+Evok2VfMczeL0iFcKqq7n7Ly4g+M4jO4V7/as0uhe8ZZgxBG+byWRiEP3eCXCZS1bBriTDmktOlzmdnDNcZzDmY/ECAVviqc9oQF070SVTSpmmFSM5CjH/QeFVIzxmQk217BSLvFoptU6wFMpJIJLwQMt19UNg1IczmK5U5jB+r2NDpdajiVWKWt3muCl55VBJOJsUuasL+WpA5Kd7jXnLL3PnlImcRpEtafKICAsqBdakl4soC1KuQRyibjNQEt0eEvxE2uOBkncPWZl61YhrlKD+VgPgPEN3jl9rNuvRoLarjNVuH7lNozsEYv3s4cL2iOCEH+ICpPi1VsH4e25Q/HB1jO46d0dOFxUG+hmEeJVCqnY4Q+50G6B0B/ysb0T0Cc5ok0ximsGJguawZGKRU6DnEiFVFCnlOM4tzuvyVEKwZu6eku8ynbNS0pUmF/TNsdnJqB7O/ZPGtc7wSaVKipMKihIdbT9Q5xKbtNZ5TjOJv316gHJiFfJna4F4nutiX1t13rxBVvml+XrlEvFIlzZKzDFzLrGhluOVyEVW45lTK94h8UvuscpeSuKugrkEyLkiA6XYVzvBJvPsEIqFhRE8A0mDEqLRl+rdauT+yfx7odmnUbqzvsLtHzvmOsZ9E2O4A3ArAtxCC1uJGQ9n1msUobLe1yaVR2fmYC+yZG4LDUKcokI1w7s4n6VyXaMGTuaNROLOMzISsWU/klu751GvfcOgjGG93PO4O0/C7DspoFuf+AI8ZfRveLx80Pj8PH2Qsz9cDemDU7Bo1dnen3PFkJCRe8kFWKULQHOxL6JDjvU9rNdsUoZbyfQnVm5lCiFwwqGfAanRUPixkaz9i7vEQeDyYSoMCnG9A5sReHu8co2AZM3a1yYn8q6I9keMUoZpBIRRByHPskRLtM1ASBOKXcYCCdHKfxSOCwxUoEKB4WikqMUGJ+ZgIoG/tuHdI3x2syVuTKjJ++xs4BpcNdoNBtNOHjB9u9ZXaNRrzW4fO4YD2dk+bRnA2WhlSutv3c4juNNk+6THAG1zuDwXLfM7NkOLLmzrpbjON7ZXfOaMWdbTThiXlfoyXLPtJjwNtVDrQcfPJl9owCtA6ht1OPxrw+joFyNjQuvpPVmJOjJJCLcOz4DM7JSsOKnfEz4119YMK4n7hrb02n+NiGhwJ0y4UBLx8S8PiPCLn0pJSoMpxRq1GubMTA1qs2i9fbiOPcqGLZnBgiA05Q5R5IiFe0qXOCOoekxaNQbvfJcsUoZrugZ57W1d0BLKmB/J/uRAUCY1dokR0HwlRnxDjdx9oZYpQwGk7kqnvP7RofLHAZoSrkYUU62U3CHTCLC+MwEjzZYdoXvECMU0jafZ2ecbaXga1f0jPNqoAhc2jriVFlDm9vsN9X2h3iVHGerNA5v7xKlQE2jXtDsJd9drAfWrspMhELWvveTArQQt+VkBR7/+hBGdI/Fd/ePduvLgJBA6xIVhjdnD8GRojos/zEPn+8+h4cnZ+LWYWk++RElJNSIRByu6pMAncEEqViEqLDO97mQSUR+q0KskIrbVSHTGsdxXg3OhLhmYDJOnahzeT9vbzJsb4yLFEX72RpvzZK54qwEvDvcLYoiRJeoMLeqGnqTv69TaxkJKugMJp+/TnKUwukG1L2TItDbSxMc3hhUoAAtRGl0Biz/MQ+bjlzE0hkDeasXERIqLkuLwpd3X47NJ8rx/345ibc3F+C+Cb1wy9A0h+lehIQapRublFoTWhWQELlE7Jc9mlxxtxBDYqQiZJZmXDMw2WZDej6eDjB2xM9519hwGJ1sShZKW0CJOQ4ysQhSF++/N1CAFoI255fj+e+PoUe8Er88PI6qNJIOgeM4TOybhAl9EvHb8TK8+ccprPqzAPdPpECNhJY4pRwD7NLQhHTqAi0qTNqm4Ii7JCIRDCbfj4Z3VNFhUpTUNgW6GV7lbvU6a8FY65dvjadYxCFSIUW9thmA4z0OvduO4Pg+6ZWocro0QSEV2xQNCWUiEYdrXWw67y0UoIWQoppGLP3fcRw4X4OnruuHm4aktrtEKiHBhuM4XD0gGVP6J+H3vHK8+cdJrPqzAHeN7YGZw7v6JLWEEG8Sizj0sqvq5qty+t4kpFS+K6N7xaFSrfdCazqn3kkRfltv5y9pMWFtUr46WteF41o2TP8ut9gvr3fDoJR2zZSGScXQ6F0XMBFCIRV3uGs2GFBPJwSodQas3lqID7eewa3D0vCvmYP9lq9NSKBwHIcp/ZMwuV8iNp8ox4c5hXjjt5OYOzIdoxMN6BfoBhJC2ogOl3ltnU9n5e7+bsGO4zi39u4y65sciRi6lni1N411bGY8aKI7uFGAFsS0zUZ8sfs83tlcgP4pkVh/zygMSAmdXF1CvMGc+jixbxKOFtfh422FuHt7CaaeMuAfo7sjq2s0zSQTQkgIEZIC2CeZKlL7SijM6Hd2FKAFoUa9Af+3rwgf5JxBQoQcb80dgiszArtfDCHBYGBqFN6YlYWbMsTYUSnF/LV7kRSpwJyR6bgxK9Vr5ZgJIYT4xrUDu7TZ048QYosCtCBS0aDDpzvP4rNd59A7UYUXpw/ApH6JNDtAiJ14pQRPDO+Lhyf3xq/HyrB+7wWs+Ckf1w5MxqwRXTGyRyx9bgghJAhRwSdCXKMALcBMJoZtBZVYv+8C/sgrw4Q+ifh43ggMTY8JdNMICXpyiRjTBqdg2uAUnK9qxPp95/HAuoOQikW47rJkXD8oBYPToihYI4QQQkjIoAAtABhjyC9twI9HLuKbAy0Vf/42vCueuq4fUqPDAtw6QkJTelw4Fk/ti0VT+mBPYTU2HSnBXZ/shVwixvWDuuDagckYlBYdFHsEEUIIIYQ4QgGanxhNDIeLavFHXjl+PHIRlWodpvRPxis3X4YxveI7XNUmQgJFLOJwRUYcrsiIw4vTB2J3YRU2Hb6IBZ/th9HEMLZ3PK7qk4DRGfFIpD0ECSGEEBJkKEDzEcYYimqasKewGjmnKpBzsgISsQgT+iTg2Wn9MTojnvKwCfExsYjDlRnxuDIjHstuHIjjF+vx14kKrNt9AU98fQRJUXKM6BaLYd1jMKJ7LHolqGiwhBBCiEfilHLQTwjxBgrQvERnMCK/QovtlWew/1wN9p2rgVprQFbXaIzpHY8F43qif5dIWgtDSIBwHIcBKVEYkBKF+yb0grbZiMNFddh3rhp/5pXjtZ9PwMQY+iZHoE9yBPokR6JPUgT6cpTr/QAAIABJREFUJEVQdUhCCCEujelNFbeJd1CA5iadwYjimiacLFPjRGkDTpY14ERZAworNYhRiDCqF8PIHrFYeFUv9O0SAamYZskICUYKqRgje8RiZI9YAC0Fe85UanCyrAH5pQ3YfqoSH28rxLkqDWKVMqTFhCMtJgxdY8PRtfW/EyLkiFfJEauU0do2QgghhHgFBWitGGOo1xpQpdahSqNHlVqHSrUepXVaFNU0oqimCUU1TShr0EIllyAzKQKZSREY0T0Gt4/qhswkFcovnEG/fv0CfSiEEA+IRBx6JarQK1GF6y7rYvl7k96I89WNKKppxIXqRlyoacKWk+UoqmlCRUPL94WJMcSGyxCvkiM+QoY4pRzR4VJEKCSIUEihkksQoZAgUiGFSiGx/F0pE0MhFUMmFlFqJSGEEEIABDhAO12hRoPWABNjYIyBMcDE0PrvlqCJoeXfJvO/rW43td7OGIPeyKA3mFr/Z4Te2PLfOsOl/9c2G6HWGVr+pzVY/lujM6BBa4DBxBAdLkWcUoY4lRzxKhkSIxQYmBqFay/rgrSYMKTFhCMqjD/dqdyvZ48Q4g9hMnFrymME7+0mE0NdUzMq1TpUtA7sVDboUNfUjAatARUNLd9zDVoDGnQGNGhb/q7WGtDUbLQ8j0wsgkTEEC4vhlwiglwqglwihkIqglwigkwihkTEQSziIBFxkIhFdv/mIBGJLP8W89xHxHEwZ1lzHAcOAMeh9f85y3+b73Tpds7qfi0bhg9IifLhWSeEEEI6L44xxtx9UG5uLuRyuS/aQwghpBPR6XTIysoKdDPcQr+BhBBC2svZ759HARohhBBCCCGEEO+jChaEEEIIIYQQEiQoQCOEEEIIIYSQIEEBGiGEEEIIIYQECQrQCCGEEEIIISRIUIBGCCGEEEIIIUGCAjRCCCGEEEIICRJeDdBMJhOee+45zJo1C9nZ2Th37pzN7Rs2bMDNN9+Mv/3tb9i8eTMAoLGxEY8//jjmzp2LmTNn4vDhw95sktd4cmwlJSW4/fbbcdttt2HhwoVoamoKRNOdcnVcAFBdXY2pU6dCp9MBALRaLR544AHMnTsXd999N6qrq/3dbEE8ObaGhgbce++9uP322zFr1iwcPHjQ380WxJNjMzt9+jSGDRvW5u/BwJPjMhqNWLZsGWbPno2bb77Z8vkLNp5ej3fddRfmzp2LefPmoaKiwt/NFsTVsa1duxYzZ87EzJkzsWrVKgCh8z1iJuT9I7ZuuukmZGdnIzs7G0uWLEFubi5mzpyJ2bNnW64DR+eV776d1aFDh5CdnQ0AOHfuHObMmYO5c+fi+eefh8lkAgCsWrUKt956K2bPnm3pR7lz387E+nweP34cY8eOtVynP/74IwA6n0I0Nzdj8eLFmDt3Lm699Vb88ccfdH16E/OiX375hT3xxBOMMcYOHjzI7r33Xstt5eXl7IYbbmA6nY7V19db/nvlypXsgw8+YIwxlpeXx/773/96s0le48mxvfzyy+zzzz9njDH2xhtvsE8//TQgbXfG2XExxlhOTg6bMWMGGzJkCNNqtYwxxj7++GO2cuVKxhhjP/zwA3vppZf822iBPDm2N998k61Zs4Yxxtjp06fZjTfe6Nc2C+XJsTHGWENDA7v77rvZqFGjbP4eLDw5ro0bN7Lnn3+eMcZYaWmp5f0LNp4c29q1a9mrr77KGGNs/fr17JVXXvFvowVydmznz59nN910EzMYDMxkMrFZs2axvLy8kPkeMXP1/hFbWq2WzZgxw+Zv06dPZ+fOnWMmk4nddddd7NixYw7PK999O6MPPviA3XDDDWzmzJmMMcbuuecetmvXLsYYY88++yz79ddf2dGjR1l2djYzmUysuLiY3XzzzW7ft7OwP58bNmxgq1evtrkPnU9hvv76a7Zs2TLGGGM1NTVs/PjxdH16kVdn0Pbv34+xY8cCALKysnD06FHLbYcPH8aQIUMgk8kQERGB9PR05OfnY9u2bZBKpbjzzjvxzjvvWB4fbDw5tn79+qG+vh4AoFarIZFIAtJ2Z5wdFwCIRCKsWbMG0dHRvI8ZN24cdu7c6b8Gu8GTY5s3bx5mz54NoGVmRi6X+6/BbvDk2BhjePbZZ7Fo0SKEhYX5tb1CeXJc27ZtQ1JSEhYsWIBnnnkGEydO9GubhfLk2DIzM6HRaAAE73cI4PzYkpOT8dFHH0EsFoPjOBgMBsjl8pD5HjFz9f4RW/n5+WhqasL8+fNxxx13YO/evdDr9UhPTwfHcRgzZgx27NjBe17VajXvfTuj9PR0vPXWW5Z/Hzt2DCNHjgTQ8rkxn8MxY8aA4zikpKTAaDSiurrarft2Fvbn8+jRo/jrr79w22234amnnoJarabzKdA111yDhx56CEBL/0IsFtP16UVeDdDUajVUKpXl32KxGAaDwXJbRESE5TalUgm1Wo2amhrU19dj9erVmDhxIl599VVvNslrPDm25ORkfPHFF7j++uuRk5ODa665xu/tdsXZcQHA6NGjERMT0+Yx5uNVKpVoaGjwT2Pd5MmxRUZGQqFQoKKiAosXL8aiRYv81l53eHJsq1atwvjx49G3b1+/tdNdnhxXTU0Nzp8/j/fffx933303lixZ4rf2usOTY4uJicH27dtx3XXXYfXq1bj11lv91l53ODs2qVSK2NhYMMbw6quvon///ujRo0fIfI+YuXr/iC2FQoE777wTq1evxosvvoglS5bYDAyZ33O+82r/t1C4Pnxl6tSpNgMzjDFwHAfA8Tk0/92d+3YW9udz0KBBePzxx/HFF1+ga9euePvtt+l8CqRUKqFSqaBWq/Hggw/i4YcfpuvTi7waoKlUKstoL9CSW27+INjfptFoEBERgejoaMuI94QJE4J2VNKTY3vttdfwyiuvYNOmTXj66afxxBNP+L3drjg7LiGP0Wg0iIyM9GkbPeXJsQHAiRMnMG/ePDzyyCOW0Z1g48mxff/999i4cSOys7NRUVGB+fPn+7qZbvPkuKKjo3HVVVeB4ziMHDkSZ8+e9XErPePJsa1atQp33XUXfvzxR6xevRoPPPCAr5vpEVfHptPp8Nhjj0Gj0eD5559v85hg/h4x8/T7pLPq0aMHpk+fDo7j0KNHD0RERKC2ttZyu/k95zuvfL+pwX59+ItIdKnb5ugcmvsg7ty3s5oyZQoGDhxo+e/jx4/T+XTDxYsXcccdd2DGjBmYNm0aXZ9e5NUAbejQocjJyQHQssA3MzPTctugQYOwf/9+6HQ6NDQ04PTp08jMzMSwYcOwZcsWAMDevXvRq1cvbzbJazw5tsjISMuFlZiYaEl3DCbOjsvZY8zvWU5ODoYNG+bTNnrKk2MrKCjAQw89hNdffx3jx4/3dRM95smx/fbbb/jss8/w2WefISEhAR9//LGvm+k2T47L+jskPz8fXbp08WkbPeXJsVl/h8TFxdn8cAUTZ8fGGMPChQvRp08fLF26FGKx2PKYUPgeMfPk/evMvv76a6xYsQIAUFZWhqamJoSHh+P8+fNgjGHbtm0YPnw473lVqVSQSqVt7kuA/v37Y/fu3QBaPjfmc7ht2zaYTCaUlJTAZDIhNjbWrft2VnfeeaelEMXOnTsxYMAAOp8CVVZWYv78+Vi8eLElu4OuT+/hGGPMW09mMpnwwgsv4OTJk2CMYfny5cjJyUF6ejomTZqEDRs2YP369WCM4Z577sHUqVNRW1uLZ555BhUVFZBIJHj11VeRlpbmrSZ5jSfHVlBQgKVLl8JkMoExhqeffhr9+/cP9KHYcHVcZhMnTsRPP/0EuVyOpqYmPPHEE6ioqIBUKsXrr7+OhISEAB4FP0+O7Z///CdOnDiB1NRUAC2j5u+++26gDsEhT47NmqO/B5onx6XX6/H888/j9OnTYIzhhRdewIABAwJ4FPw8ObaysjI888wzaGxshMFgwIMPPojRo0cH8Cj4OTs2k8mERYsWISsry3L/RYsWoW/fviHxPWLGd4wZGRmBblbQ0uv1WLJkCUpKSsBxHB577DGIRCIsX74cRqMRY8aMwSOPPOLwvObm5ra5b2dVVFSERYsWYcOGDSgsLMSzzz6L5uZm9OzZE8uWLYNYLMZbb72FnJwcmEwmLFmyBMOHD3frvp2J9fk8duwYXnrpJUilUsTHx+Oll16CSqWi8ynAsmXL8NNPP6Fnz56Wvz399NNYtmwZXZ9e4NUAjRBCCCGEEEKI52ijakIIIYQQQggJEhSgEUIIIYQQQkiQoACNEEIIIYQQQoIEBWiEEEIIIYQQEiQoQCOEEEIIIYSQIEEBGiGEEEIIIYQECQrQCCGEEEIIISRIUIBGCCGEEEIIIUGCAjRCCCGEEEIICRIUoBFCCCGEEEJIkKAAjRAfeuqpp3DmzBm/vFZZWRnmzZvnl9cihBBCnKHfP0I8RwEaIT60a9cuMMb88lpJSUlYu3atX16LEEIIcYZ+/wjxnCTQDSDE13bv3o3//Oc/CA8PR0lJCdLT0/Haa68hKiqK9/5vvfUWzpw5g5KSElRVVWHKlCl4/PHHwXEcvv32W6xevRoAkJ6ejhdffBHx8fHYtGkTPvzwQ4hEIkRHR+Nf//oX1q1bh/Lyctx///348MMPkZaWxvt6er0ey5cvx87/3959xzdd7f8Df2U1HekelAIFSilTKFv2EnGDXpGhRUXBcR1XXBcXDq6CX7n3dxEXijivFxS9ooIbKCAIFAq00EKBDlrapjtNsz/n90dJyPgk+SRNk7R9Px8PH9Lkk+R8Tj5JzvuM99m/HxKJBLfeeiuWLFmCkpISvPDCC6irq4NMJsNjjz2GyZMn4+uvv8auXbug0WhQWlqKjIwMvPHGG1AqlZg/fz727dvXbnVJCCGk46DfP0I6KEZIJ3fgwAE2ZMgQdurUKcYYY88//zx75ZVXnB6/bt06dvXVV7PGxkbW0tLC5s6dy37++WdWWFjIpk6dyqqrqxljjL377rvswQcfZIwxNmPGDHbhwgXGGGMfffQR27lzJ2OMsenTp7OioiKX5du0aRO77777mMFgYGq1ms2ZM4eVlZWxefPmsW+++YYxxtj58+fZhAkTWGVlJdu6dSubMGECq6+vZ0ajkd1yyy3sxx9/ZGVlZWzChAltqitCCCGdB/3+EdIx0RRH0iVkZmZi4MCBAIC//OUvOHDggMvjr7/+ekRFRSEsLAzXXnst9u/fj0OHDmHKlClITEwEACxcuBD79u0DYwxXX301srKysGrVKvTv3x/Tpk0TXLYDBw7gpptuglQqRXh4OP73v/8hNjYWhYWFmDNnDgCgT58+GD58OHJycgAAw4cPR0xMDCQSCdLT01FfX+9FrRBCCOns6PePkI6HAjTSJUgkEsu/GWM2f7s73vw3x3E2tzHGYDQaIRKJ8PTTT2PDhg1ISUnB6tWrsXbtWq/KBgDl5eW88/YZYzCZTAAAuVxuuV0kEvltnj8hhJCOhX7/COl4KEAjXcKxY8dQVlYGAPj6668xZcoUl8f/+uuv0Gq10Gg02L59O6ZMmYJx48Zh9+7dUCqVAIAvvvgC48aNg9FoxFVXXYWwsDAsWbIES5YsQWFhIYDWHx/zj4oz48aNw/bt28FxHDQaDe655x5cuHABGRkZ+PbbbwEAxcXFyMnJwahRo9paFYQQQroQ+v0jpOOhJCGkS0hMTMSLL76IiooKDB48GH//+99dHq9QKJCVlYWmpibcfPPNmDx5MgDgsccew5IlS2AymdCzZ0+8+uqrkEqlWL58OZYtWwa5XA65XI6VK1cCAGbOnImHHnoI69ats0wxsbdw4UKUlJRgzpw54DgOt99+OwYOHIg33ngDK1eutCzK/sc//oGUlBQf1gohhJDOjn7/COl4RIzGhkkn9+eff2Lt2rXYsmWLoOPffPNN6HQ6PPHEE+1cMt8qKCjAww8/jF9++SXQRSGEEBIE6PePkI6JRtBIl/T444+jqKjI4fbevXujf//+fn29devWtfn5Dxw4gMcffxx33nlnm5+LEEJI50W/f4QEPxpBI4QQQgghhJAgQUlCCCGEEEIIISRIUIBGCCGEEEIIIUHCqzVoubm5NvtQ6HQ6m787i856XkDnPbfOel5A5z23znpeQOc9N1+el06nQ2Zmpk+ey1+OHDmCsLCwQBej0+isn5NAofr0LapP36L6vMzV759XAZpcLsegQYMsf586dcrm786is54X0HnPrbOeF9B5z62znhfQec/Nl+d16tQpnzyPP4lEok75vgZKZ/2cBArVp29RffoW1edlrn7/KIsjIYSQToHjOLz44osoLCxESEgIVq1ahd69e1vu/+ijj/DDDz8AAKZOnYqHHnoIWq0WTz75JGpraxEREYE1a9YgLi4uUKdACCGE0Bo0QgghncOvv/4KvV6PzZs34/HHH8fq1ast95WVlWHbtm3473//iy1btmDv3r0oKCjAF198gYyMDPznP//B3Llz8fbbbwfwDAghhBAK0EgnxXEMSpUOhZUqFFQ2obhGDZ3RFOhiEULaUU5ODiZPngwAyMzMRF5enuW+5ORkfPDBB5BIJBCJRDAajZDL5TaPmTJlCvbv3x+QshPiqXq1HiaOdkoipDOiKY6kU9AZTfjzXB22HqxB8W97ceqiChxjiAkPgVgEaPQmqHRG9IoLw5jecZg6IBGzBndDeAh9BAjpLJqbm6FQKCx/SyQSGI1GSKVSyGQyxMXFgTGG119/HYMHD0bfvn3R3NyMyMhIAEBERARUKpXb1+E4rkOunQtWWq2W6tMLu841o09sCPrEhtjcTvXpW1SfvkX1KQy1TkmHVlStwmcHSvHN0XLERYRgRJIUf52ejit6RqNbZCjEYpHl2CatAYWVKhw8X4cP957Hiq9PYO6IHrh3Ul+kJSpcvAohpCNQKBRQq9WWvzmOg1R6+WdOp9PhmWeeQUREBFauXOnwGLVajaioKLevIxaLaZG7D1HSAO+c1pWjd1IkBqXYXrNUn75F9elbVJ+XuQpUaYoj6ZBOVjRh2SeHMWf9Phg5Dp/dMw6/Pz4VS8fE4+ohyegeHWYTnAFAVKgMY/rE4a/T0/HtQ5Pw3cOTIBGJcMObe/HEl8dQ2agN0NkQQnxh5MiRyM7OBtC6HUxGRoblPsYYHnzwQQwYMAAvv/wyJBKJ5TG7d+8GAGRnZ2PUqFH+LzghhBBihUbQSIdSrdJi9Y4C/Jxfhbsn9sHqvwxDXESI+wfy6JeowCtzh+LhGelY+/NpXPXP3XhsVgbumtAHErvgjhAS/GbNmoV9+/ZhwYIFYIzh1VdfxaZNm5CamgqO43Dw4EHo9Xrs2bMHALB8+XIsXLgQTz/9NBYuXAiZTIa1a9cG+CwIIYR0dRSgkQ6BMYavj5TjlR9O4urB3fD7E1ORFBnqk+dOigrFmluHYeG4VDz91XH8cLwCby4aiR4xtBEtIR2JWCzGyy+/bHNbv379LP8+ceIE7+PWrVvXruUihBBCPEFTHEnQq2jQ4O6PDuFfv57G+oUj8fqtw30WnFnL7BWD7x6ehFG9Y3HDuj34vaDK569BCCGEEEKIKxSgkaD2e0EVrv33HvSJj8BPf5uCSf0T2vX1QqRiPHv9YKz5yzA8tvkYXv+xAEYT166vSQghhBBCiBlNcSRBieMY1v1+Bh/uPY9/zc/EzEHd/Pr6Vw9JxsDkKPz1P0dworwRb90+ElGhMr+WgRBCCCGEdD00gkaCTpPWgGWfHsYPxy/i24cm+T04M0uND8eX94+HQi7Fbe/uR0WDJiDlIIQQQgghXQcFaCSoVDVpces7f0AmEeN/f52IvgkRAS1PqEyCtxaNxOT+Cbj57X3Ir2gMaHkIIYQQQtqDzmhCi94Y6GIQUIBGgkhxjRp/eecPTOiXgLcWjUSEPDhm4IrFIjx7/WD8dXo6Fmw4gH1FNYEuEiGkA6lq0qJJawh0MTolE8fw3bEKMMYCXRRCOrz9Z2vxy0lKkBYMKEAjQeFkRRNufXc/5o3qhZU3DnbYZDoYLB7fB/936zDc92kOfqUvMEKIQAfO1WJnQXWgi9Ep6Y0cOMbQ0eMzrcHk0fEc13rC52qase1YRXsUiXRBWkNwJ0VT+bCjK7+iMaizdVOARgIut6wBC98/gIdnpOPRq/pDJAq+4MzsmqHdsX7RCDy2ORff0Y8iIaSD8WUDx9cYY/g2txwmroNHWx6qVmnxU36l2+OadUZL3ewsbA34TRyj0cMgUVyjDtq16hq9Cd/mlge6GF5r1Big0hrwuw87uiobtVBpg3c6JwVoJKBOXWzCXZsO4tnrBuHOCX0CXRxBpg1Iwgd3jsYz35zA5kOlgS4OIYQI0tCi92kDx9fMcZmhi21tojcKO9/fTlXh1MUmAK3BWlemM5qw90xwLTc4dqEBOSX1gS4GL0+vlx+OX7T8u7ZZ5+vieERnNGFXYbXgz4m/lTdo2qVTiQI0EjDnlM3I2ngQj87sj9vG9Ap0cTwyLi0en90zDq/tKMB/D1KQRggJfgZT8I60tEfjS2sw+a1xeaZKhbPK5nZ/HWfBq85owvfHu86sjiaNEbXqwAYOnQXHMfx6surStNnW7wgj13qdVau02BvgdffmAeJg/fY6XFyHklq1z5+XAjQSEOUNGtzxwZ+4a0Jv3D2xb6CL45XhvWLwyZKxeG1HAbYcLgt0cQghpEMqbzJgR95Fm9sYYzZB28VGDSobtZa/fztVhUaN6+maOSX1fmtcnrzYhLzywGX5VetMnWpqaGWj1mYUx58MJg7f5pZb1vl5Q2c0obEleKcTWzNwHNR6I4w85xtMs2d9XRZfrqbh2qGiKEAjflev1iPrgz9xY2YK/jo9PdDFaZNhPWPw8ZKxWPX9SWzNuRDo4hBC/IAxBqVKB8b8t/6nI2eB1Bs5l+tfNDyJCc7XqG2CtoPn6/Dn+VrL3806I5Qq1yMovm40GU2cx8k8gl2jxoCGFr3NbRzHHG6z5uy6L2/Q+Kx+lCqdZRTH38yBbluun6OlDdh1OjDTibUGU6eaAuuLQMpo4qAz+ubaPFpaD43e9rnyK5rcfh95igI04lc6own3fZaDzF4x+Ps1A4M6IYhQmb1i8NGSsXjpu/wOvQiXECKMzsjwx9ka7D6tRLYH62B2n1Z6NcpRp9Z36CyQ2ksNI0++H1v0wRcIHSyuE5TMI9i4WtO3q7Aau08rbW4rrlU73GZt/7laS5ISa4eL63Cmqv2neXYEfJ/zqiZtm0blhNp/rha/nXKfnZDjmM2odLBjAiY5mrjWREP2HQgHztXhxzzffHZL61pQ2eRYb77uvKEAjfgNYwzPfJ0HMOC1v1zRKYIzs5GpsfjwrjF47pu8oE7bSghpO/NPP9/ogysNLXqvfsTNjb0dJwIz5atRY0BZXYtPnstVkBbs2QjVOtv3rrJRG7SJC6xtP3ER1TwNSmeq3YwEKFU6p9nvnDWiGWMebYDsqnlwpLTecq2odUaotAZUq9y/Fz/mXfR4lENnNKFR6/iZPXWxyePPxIFztSj3YZbHi40aXKh3LINB4DVZ3qCxGZUG2mdKY00b14GK4L6tyBhDeYPG0hlhHwfzjSgKeV4hr+vOt7nlNp+/vWdqUK92/7tBARrxm7d3ncXhkjq8mzUKcqkk0MXxudF94rBu0Qg8+kUuDp6vC3RxCCGdjL6N2Q096ROraNCguKZ14XtuWQOOlLrOTsfXAKps1LrNyKjSGnCh0bvpm/bn09Cix87CalSr/DMq8Of5WpyuUjm9f19RjUcBfHvS8HQM1PE0Es/XqFHlQTAnVGldi2UD5D+Katq0PqusrsUywvpbQTV+L6jG/rO270Vts86h8awzch4HaLmlDTha4RhUna5SIb/C8zWHvgyADp6vc5k1kuNcb1vhj/4Qo4nDvqIan0wvdJXkqFFjwOHiOhRUOv88+oJ1fZ9TCt+DsNbqs1ar1uGigJFLCtCIX+w4cRHv7zmHD+8ag7iIkEAXp91MH5CEVTcPxdJPDuNkRVOgi0MI8TOOYwFL1qA3coJ6dK2DBusAynqEI6ekHscuNAh+7d9OVeFio21D9s/ztSiqdj3lrcGHiRTyypvQpDEgp7geRdXNDg3QerXe7TTLi42aS1P12t7Qq2nW4UK9pk2jbLXNOkFBnvl9b2wxoK7FCK3B5HbPu3qe5z1nlYnyYqPGZyM+1nWgbNahotE3z+vset9bVIOyutbX2HasgjcYdfvcAEweRjHuOkHcJbbxJXPZ23PbityyBst2B3xBmLn22hIMmuv0cLH7ju+2ZFM0mDi3393WI5b272VuWYNPp7BSgEba3ekqFZ766jjWLxyJfomKQBen3c3J7IHlszKw+MOD7ZJ6lRASvA6X1OMHH09FFLL2AgB25F10GxABrWvhzIketp+4CK3BhIoGjWWEw57QUSCDkb+cvpzM3qgxOE2AYG6ImhhDfkWjw5Q6IYkTDp6vQ3mDBicvOnaweTMF86yyGT+2Yd3a3qIal+vBzMrqW4ORP87W4HilFn+er/NqzzvrOjp4vs5lo1ijN/kkkG0v5iQfjDHUtUNKfm+CjnM1jp/P01Uq5JYJ7wxpT/an5O6ze6G+BbVqHZQqnc/WeAlVWKly2vnRrDW6X/Nqd3I/5lXij7PeZ33lGIPOh1OeKUAj7UqlNeD+T3Pw1xnpmNQ/IdDF8Zs7J/TBHVem4o6Nf3o0958Q0jEdK2vAOWUz6lv0LhvyfFPNfMHcUFELTK4hEoksDUwTx9o0yuNNr/EfRTUob9DYTP0RYldhNXbxJKhwR6U1OB290BlN+ONsjVfrh4VMpzRfDw0tepyvUePPc7W8x5XUqvFtbrlXIx7mzkDzNFghQXWzk3VkQp2vUfMGsq6Y12DajzTllTf6bP22/UiuNc5N5lWOYzaJYFytU9KbOJev5YxSpbMJyE5XqZx25ja2GNyOhGr0JtSp9Zbpm+bvGOuSu+o5XoioAAAgAElEQVTkMV8HnkyB/qOoBmqrYN5dYCPkuRs1BpyvUeOPohpBn/GCyibegBcAVDrPRyo5xtCk8WSdpPOTulDfYhMg2h8qpNONAjTSbhhjeOLLY+jfTYH7pqQFujh+9+jM/pgxIAmLPzzo12kNhBD/K65V47RVBjtnozWu1ozwqW7SYv9Z/ga9tTwB62F8seZEa3Dcb+t0desoCl+jw9mUIWWzDmV1LV7NMjA/p9F+TcqlRhDfef5eUO10I+mGFoPLxBeuOHsM33nvPq3E8QsNvBnggMvJFPw1Rba4HWZ4OLvGzA3UU04CuouNWq/qn4+rNeBnlc3Y7yRABlxPadQbHbdZ8CbIPeMiILO363Q1dha6HkE9WFyHPWeULoMkV98hQkfvDKbLa/iUzTpB6/msq1NrMOHb3HIUKLWWAMxg4iwp63cVVuP4hQYom3Vu20zmbUcKnaw58/R71h2+NZOuAn3z69un4/cEBWik3byXfQ5nqprxxrzhnSpjo1AikQgrbxyC9CQFHvgsp0Nk+yKEuMe3oauZuQFncrGg3bph821uOfacudwAs2+YOGvM27MexcotaxA0nc+br+Wf8isdEoa4aoQ4C4pcybaa0ucq819BpbDRG39uf1LVpIVSpcP3x4UlD+gITBxz2hB2Jfu00uMA3FmW08YWA/afrXUa4Ll3+WKvab48usgYE5xZdW+R0ifbLCg9zGrobmptW/f7Mydccfc8Z6qa3Y6UqXVG1DTrHNo7P+ZVWkbcKlVGy/fcofN1+PmkZ3V6uLjOr9uOtOiNgve00xhMNnVk/f2n0Zs8yuJLARppF/vP1uKtnUV4L2sUIkNlgS5OwIjFIrwxbzj0Rg4rvj4R9GmkCSG+4WwKi9buBxywzaa3q7DaMpKi0Zu8WtNQUqvmTb/tKiDzJFhT64w232XWHXD2KaUBx95sZyNE5tEL62mgv5x0TD7ijLkUQr5lGWOC1uvZPEbAMQfO1bZpHQvQWj/Opo2Wepja3VzXuWUNHp+vWU2zzhIMe5KVsr5Fjwv1nk0BtC9jY0vr1NQLDS2oVmldZs30ZOTRfOxZpdpuSqNz9tss+Jsn2xRYE9rs+DGv0uX0WnMA5ypFfIvehH1FNbwZLvlGJy3BoQfvnX3imvIGjaB1mmZHS+txsVHD+17zbY7urGh8Aw/7impsOuCsv79/PlnpUWBJARrxuZpmHR7971GsmjsU/btFBro4ARcqk2DD4tE4UlqPdb8VBbo4hBAfc5XFTCjrdSZqnRHlDRr8fLISFU6y6DHGeBtK9p1AWoPJ5xllGzUGbDtWAZ3RBJ3RZOktNzdklM06m6DTflTQ2eiSdZIB6+l3QmYfGDnm0XSiJq3RqzTp5vIIHb0zc9U5Z27oaQ2t57mzsBp7itoW5JnllV8+R1ejmU1u1jkBrde3s4aws1FVldaAqiYt8r28Bnedbl1zKCTIKK5VC1q/xxjD98crsPu00qNrQEgIwTHbTH/1aj1Ka4UH1VqDibdzBWjtrBDyPpm118yl7DPugyGG1uB6R97lhEmupll+14YRZ40HgWtpbQtK61r8thWSOZg0r2U0rxEVskE4BWjEpziO4fEtxzBtQCLmZPYIdHGCRlxECDbdNQYf7y/G1pwLgS4OIcSHankyxJmnQQoNGn4vqLaMHuSWNTht1NWpWze7rmjUIvuM0pKkwr4tZp5OVN6gwZlqlUM5rAOGsvoWQQ1go13j12hi2FlQbRnhOn4pLX9RdbPNtE1veLoOK7+isc37xAl1sVFjM+XvmIA1PEL2Z7Jeg+ar/dOs9zRzNb1qZ0G10+DGPCLwUz5/Eo9qldayabOR49CkNVg6K3RGzuWGzkqV435lfOyDS2eP2O4kgyrfeim+Om7rHBfGmM2IcU5JPY6WOV8PVVyjtrnWCypVyCmpxy8nq3jfr4oGjUOg42kYllfeiF9PVlnK1R7zesrqWpwm8BDK1Wipt6zfC+uOI75OoG3HKlDeoPFoCqnQI5t1RvzhphOGAjTiU+/vOYfyBg1evGlIoIsSdPokROD9xaPw4rb8Nk+BIYQEt32Xfnw9WV9RZtVz7iw74J4zShwpqbescTtcbNv4M0+B010ajTG3LexH974/ftGS8bGwUmUJMquatCipVdukKDfji5naklZa6IbIQl7D+hh3Df4aD9cAuSMk2Ya309PaypOsoccv8I8mmYMj+3o1B+bWs8Iu1GsETeMyB7h/nK1xCJ5MHPPJxsbA5a0XhG5efqbK86BCrTM6HdlSX3rfnQWp1nsNFteoLWv2WvRGmwDC/PjS2hbB52KuQ7XVtac1mHBW2Wxzm7PPy7e55ZZOD45jlumIljLVtzh0pFhP7fZ0Oq4979cb2tp/ttbt9WSfKbO4Rg3GGGqbdbhwaT89Z7MZrNnXpavvIndrEaVuX40QgY6U1uPN34vw1QPjER5ClxafUb3jsPovw/DAZ0fw1f3jaQooIV2U1mCCXOpdH6n9Tz5fEgyVmyQhHGP47dTlEZGmS9noDthluNtz5nJnkpBMke1BSHY/T5b37iuqwbQBSW6P+za3HMN6xqBvQgQA1yNQ7pgDZrOyuhb0igsX/Hj70cu2yK9oRAnPCK2no5ZCEtEAjmuG7JXVaywBrLJZh7M+SrBy7EKD4M3WOa41tPBm5HJnYbVD3dl3AvCNoNqvrXQ1DdR+Cqkz1kFE+aX1f/usRmo8TXJinqp8VtnsUF5vNv9uD64++4wxQQFtk9V3TG5ZgyVQtn7uQwI2ynZaDi/GKWkEjfhEk9aAR744ihXXDcTA5KhAFyeoXT+sOx6Y1g9LPj4UNF9whBDfc7Ue5qf8Spyr8S7NuXWjwVlGSfPm02buFtE3OUlrXW/VYHU1VS1QzPtwedMAEoJvg2ODiwydztj3lh8prbc0AtUCAp1TF3033auoupn32vQ2iZW3dW8OeK3XXDm7DtvboeI6/HaqyuWohrP6ERLY8j3Wfh2UfYIK+84S4PIotnUiHuuRNusgolpAGnx3zl/6jhK63519R4QreiPHm5TDnjeZYM2EBrbHrQL5cqvENsW1ap/sXelNFlQK0IhPrPw2H0NSorBobGqgi9Ih3DclDVf2jcf9lH6fkE7L3RQdb/fIsV7z5qpR3dpYu3x/W0aAAsv5OZp79fkahq6mRh661Dhu0RttglA+HMdsAhpvk4vYM08pdHadWI/mGAQ0ZAOF7zfM1QbPZq72I/O3yiYtb0PcvLbT2e+00GDdVxvUm6fqCak7X0/lFcKTdWM78i66nb68s6Daq+DGzNMsooBjoMyXtOX4hQav9nD0BAVopM22HavAvqIavHbLsC6535k3RCIRVt08FIwxPP+/PEq/T0gn5Gq/NF4eHC501EKjv9zY8PXmra429fUlbxpZ7pjX4PxyssplA5Ax4IcTF50mn2gLd8kHrEcOrA/1VYBoz9sgQuhGx/av0ZbRsnNtGFURynoNlXU2Qmtnqn2fyIJc5knWSm+d93ImQ3ujAI20SUWDBs//Lw9vzBuOuIiQQBenQ5FLJXjnjlHYW1SDjXvPB7o4hHR4HMfhhRdewPz585GVlYWSkhKHY+rq6jB79mzodK29y4wxTJ48GVlZWcjKysLatWt9Vh538Yt9A92TxohW4FSitmZSc0XIonl/aa+GnKdZ3DzlalTTOjC1TnLg7X5mwOUtIfjYb4fQFlVO1v3wZVP0htD1b21hvUWBNV8lryDBwVcZU32NMjkQr3Ecw/ItubhlZA9MyUgMdHE6pASFHBvvGo157+5HvyQFpgtYuE4I4ffrr79Cr9dj8+bNyM3NxerVq/HOO+9Y7t+zZw/Wrl0LpfLyeqzS0lIMGTIE7777rs/L465h35ae2/PtGHgJ1ZapR77SOpWu485AOCJwVNNXgQ3flhDtwZO9vzqa9kj/3ln4I3C2J3R9nDNC9iQLBBpBI177YO851Kn1ePqagYEuSoc2MDkK/7otE49+cRRn6IufEK/l5ORg8uTJAIDMzEzk5eXZ3C8Wi7Fp0ybExMRYbsvPz0dVVRWysrKwdOlSnDt3zmflac8Rprakt+9MhCQZCGbuUm13VC0ddr0j6Wr8tX+ip2gEjXjlZEUT1v1WhC33jUeoTBLo4nR4Vw3uhgenp+PeTw7jfw9ORCxNFyXEY83NzVAoFJa/JRIJjEYjpNLWn7qJEyc6PCYxMRHLli3Dtddei8OHD+PJJ5/E1q1bXb4OYwwV5b5JBU4Ag0FP9elDVJ++U1FO9elrVJ+XZaQpnN5HARrxmNZgwqP/PYpHZ/bH4BRKqe8r901Jw+kqFR74PAefLBmHEC/3SCKkq1IoFFCrL08b5DjOEpw5M3ToUEgkrZ1Mo0ePRnV1NRhjLhMeiUQipPTo7ptCE1SUVyClR0qgi9FpUH36FtWnb1F9WnM+PZNagMRjq3cUIEEhxz2T+ga6KJ2KSCTCa7dcAYOJYeW2fMrsSIiHRo4ciezsbABAbm4uMjIy3D5m/fr1+PjjjwEABQUF6N69O2WjJYQQElA0gkY8svu0Ev/LLcf2RyZDLKZGjK/JpRK8lzUKc9bvw0d/FOPuiRQEEyLUrFmzsG/fPixYsACMMbz66qvYtGkTUlNTMXPmTN7HLFu2DE8++SR2794NiUSC1157zc+lJoQQQmxRgEYEq1Pr8eSXx7Bq7lCkxIQFujidVoJCjg/uHI3b3tuPtEQFplKGTEIEEYvFePnll21u69evn8Nxv//+u+Xf0dHR2LBhQ7uXjRBCCBGKpjgSQRhj+PvW45iUnoAbhtHc4fY2qHsU1s4bjke+ONqmPW8IIYQQQkjHQgEaEeTLwxdw8mITXpozJNBF6TKuHpKM+6am4d6PDwXtRoqEEEIIIcS3KEAjbpXUqvHK9yfxz9syERkqC3RxupQHpvbDyNRYPPDZERg5ShpCCCGEENLZUYBGXDKaOPxtcy4WT+iNsX3jAl2cLkckEuHVW66A3sThrQM1lNmREEIIIaSTowCNuLR+ZxFMHMPfrnKfrpq0j1BZa2bHoxc1+GDP+UAXhxBCCCGEtCMK0IhTOSX12LjnPP41PxMyCV0qgZSgkOOlGcl48/cz+OVkVaCLQwghhBBC2gm1ugkvldaAv20+ihXXDUK/REWgi0MA9I4NwbqFI7B8Sy7yKxoDXRxCCCGEENIOKEAjvFZuy8fA5CgsHNsr0EUhVqYNSMITVw/AvR8fRnWTNtDFIYQQQgghPkYBGnGw7VgF9p6pwZq/DINIJAp0cYidOyf0wazB3bD0k8PQ6E2BLg4hhBBCCPEhCtCIjQv1LXjumxNYe9twxEWEBLo4xIkXbhiMqDAZHv8yFxyl3yeEEEII6TQoQCMWJo5h+eZjmD+mFyb3Twx0cYgLUokYb90+EmeqmvGvX08HujiEEEIIIcRHKEAjFu/sKoJKZ8QTswcEuihEgKhQGTbeOQb/+bMU3xy9EOjiEEIIIYQQH6AAjQAAjpbW493d57BuQSbkUkmgi0MESo0Px7tZo/D8//JxuLgu0MUhhBBCCCFtRAEaQWOLAQ9/cRRPXzsQ/btFBro4xENj+sTh5TlDsOzTHJxTNge6OIQQQgghpA0oQOviGGN44qtjyOwVgzvGpQa6OMRLt4zsibsn9MGdmw6iWkXp9wkhhBBCOioK0Lq4jXvPo6i6Ga/dcgWl1O/gHpqRjqkZibh70yE064yBLg4hhBBCCPECBWhd2NHSevy/X89g/aIRiAyVBbo4pI1EIhFeumkoesaG4f5Pc6A3coEuEiGEEEII8RAFaF2USmfCQ/85imeuG4QhKdGBLg7xEYlYhH8vGAGtwYSnvjpGe6QRQgghhHQwFKB1QYwxrN2rxOg+sVg4tlegi0N8LFQmwQd3jkZ+RRPW/FgQ6OIQQgghhBAPUIDWBb2/5xzKmwx49WZad9ZZxYSH4OMlY/FtbgXezz4X6OIQQgghhBCBKEDrYv48V4s3fyvCM9O6IUIuDXRxSDtKiQnDp/eMxbu7z+LTAyWBLg4hhBBCCBGAArQupLxBgwc/P4JVNw9F39iQQBeH+EH/bpH4eMlYrP25EF/lXAh0cQghhBBCiBsUoHURGr0Jyz45jFtH9cSczB6BLg7xo6E9orHprjF46bt8fH+8ItDFIYQQQgghLlCA1gUwxvD01uOIV8jx1DUDA10cEgAjUmPx/uLRWLH1BHacuBjo4hBCCCGEECcoQOsCNmSfw/ELDXhzwQhIxJQUpKu6Mi0e72WNwlNfHcd3x2gkjRBCCCEkGFGA1sn9mFeJt3YWYcPi0YgOp82ou7oJ6Qn44M7ReOabE/jmKK1JI4QQQggJNhSgdWK5ZQ148stjeOv2kcjoFhno4pAgMS4tHh/dPQYrv83HlsNlgS4OIYQQQgixQgFaJ1VW14J7Pz6E524YhMn9EwNdHBJkRvWOw8dLxuLV7afwwR7aJ40QQgghJFhQgNYJNbYYcPdHh3Db6F6YPyY10MUhQWpEaiy23Dce7+85hzU/FoAxFugiEdImHMfhhRdewPz585GVlYWSEsf9/+rq6jB79mzodDoAgFarxcMPP4xFixZh6dKlqKur83exCSGEEBsUoHUyGr0JSz85jEHdo/DE1QMCXRwS5DK6RWLrAxPwU14lVnx9AkYTF+giEeK1X3/9FXq9Hps3b8bjjz+O1atX29y/Z88eLFmyBEql0nLbF198gYyMDPznP//B3Llz8fbbb/u72IQQErREIkouFwgUoHUiBhOHv/7nCOQyMdbOGw4xZWwkAvSMDceX949HfkUT7v8sB2qdMdBFIsQrOTk5mDx5MgAgMzMTeXl5NveLxWJs2rQJMTExvI+ZMmUK9u/f778CE0JIkEuICAl0EbokaaALQHyD4xie+PIYGlr0+OzecQiRUuxNhItXyPHfZVfi4S+O4rb39uPDu8agW1RooItFiEeam5uhUCgsf0skEhiNRkilrT91EydO5H1MZGRrEqWIiAioVCq3r8MYQ0U5bVXhKwaDnurTh6g+faur16dBIYWy2Xcdt129Pq1lpCmc3kcBWifAGMNL3+Wj4KIKm++7EuEh9LYSz0XIpXh/8Wi88v1JzH1rHzbeOQaDU6ICXSxCBFMoFFCr1Za/OY6zBGdCHqNWqxEV5f6aF4lESOnRvW2FJRYV5RVI6ZES6GJ0GlSfvtUR6jM6TIZGjaFdnntw9yicvNjks+frCPXpP87rlYZZOjjGGFbvKMDOQiU+vWcsYsJpKJp4TyIW4cWbhmDZlDTM37AfP+VXBrpIhAg2cuRIZGdnAwByc3ORkZEh6DG7d+8GAGRnZ2PUqFHtWsauZPrApEAXgZCAk/hhuUnP2HAo5MHVOd8rLjzQRejQKEDrwMzB2Y68Snyx7Eok0ZQ04iN3T+yLdQtH4KmvjuOfv5wGx1GGRxL8Zs2ahZCQECxYsACvvfYaVqxYgU2bNuG3335z+piFCxfizJkzWLhwITZv3oyHHnrIjyUOjPFp8ZCK+X/+Q2USn77WtAz3QZpMQuul+aQlXJ7+FOFiZoycljRYDEz2fNbHoO7CHzMyNdbj55/hRUdFsoftuUDl8UiKdF7ORIXcjyXpfIIr3CaC2QdnPWLCAl0k0slMH5CEb/86Ecs+PYyTFY345/xMRIXKAl0sQpwSi8V4+eWXbW7r16+fw3G///675d9hYWFYt25du5fNXlJkKKpVWoxPi8f+c7X+fe2oUAxIViC/wnF6zbQBifgxz/uR8+uu6I6Ciyqcq2kGAESGum9mtEd8cUWPaJwob/T9E/PI7BWD3LIGnz8vw+WOsczUGOwrquE9rnd8BE5XuV87mRoXjtK6Fp+Vr71IxCKYvOgUlEnESFAIm0Ukl4qhM3qetVjaxs4EoVMRO8pSlSvT4rDtGP96ssTIywFaVKgMTVrH8549JLnNM3XCZBJoDCab224cloLvjgf3OjeZxPUXH3W7dECMMbxGwRnxgz4JEfj6wYmQisWYu34fTvlwHjohnUVUWGvHhbPRp9lDkh1uG9pDWK+9dSPHt/gbmnJp20bQZBIxrugZbflbLBahZ2zrVCdPRwXaIi3R+eJ7ZyalJ3j1Wu2xhaT9+xAeYvu3uyROUzMSAQAhbhqBZkKP84ehKdG8t0eESDGlf+t5mf9vvazjuiu6I95u1Mb82bRnnTre26mBczJ7OL3PPKopk4htro9pA5Iw+VLZfa0fzzXfNyHC4TZno+euXDuUf82tsxT84SFSm+/DAcmRTh7vcVEcxPJkmbTOYj6qt/NRT29GXAHv6tBsfFo8AGB0H9ejscHziSSCmDiGZ77Jwy8nqyg4I36hkEvxzh0jsXBsKua9ux+f/1lCm1oTcsm1Q7tDeqkx4KytwRe4RZpHo900UKwb6tZT3twRB3jvosxeMYi0a/gO6XG54W0fENw4jD9pgKtGsKfMQQufnrHhiPMynbhMKnYIpLtHh2Fc33i3j3U2dfGaobZBfZjdNWR+d92V2bqXPvZSMMMXjF17RfskvXG2/qpnLH/b5YZhKejDE1QAwOSMBMRGhOCm4SmWRvmAbvwNfzMh0z9TYsJw47AUQe+XEGkJCmT2am18X8PTOWMfbPNRhEod3qfu0ZfrrHe8bR3FR4SgT0KEJUi7alA3AK0jyd68vr0QqdhloGNvbN84AHA760Z2KdBxFkj7grmDiI+zwNGdMX1ikdkrxuF2T64hd3VDAVoHojdyeOSLozhW1oAv7x9PwRnxG5FIhKVT0vDJPWPx9s6zePiLo1DxTFcgpKux3tJE6PotvjUpzkaXzA3c6500oPl6cuMj5FAImFpoz5ORhMFu1u30jo9w6F139fxC9u20Hg3w5Pdv5qBu6BMfgZjwEKfTioZ4kbHWHByJAEzoZzv6NrZvHJKj3Y8Y8iVSsB9x7Z8U6VCXGZcCkzCZxGEUYki3UN5RPfP0PPtRJj4pMWFItSubN6NsQkYarBvQrhJqmDsr7OsiTCZxvkaPtX4erI/lIxaLkBwdaqlXb5jXY0WHyWA+DbFYZAmIzCMnQvCNfFlXTcylgMYcgJpHEoUsQ4gOtz1mKE8Qx8f+O8r83pqDxagwGUJlElw1qBui7WYVOOtIEItFmJPZA5OtRq/5AiqxSMQ7ws03O8EZmUTscE0Dl0ceJ/RLwNWDk22O758UaSl7VJjMpgz2QTIAQZ95oShA6yBa9Ebc8/EhKFU6/Pe+K5FAiy9JAIxMjcX2RyZDb+Rw7b/34E8/r50hJBiZR8OuTIu39FzbiwkPsTRw+RqTUqvG76T0BIzvF49pA5J477c2c1ASesaG2YwUTOjnm5EAV1La2EEoEtk2nF0xjwZaN/bDQ6RugzRz8geFXIrhl3q7BzrpMfck0579sREuAk93DUi+BCnmRq05yDJvd2LdO28eQbIPzlJiwpAY4TrQ7psQYTNlzTxd0Dqhy5g+cZb6czW6ISQJjDMznXxWPDV9YBKmDbg8Ojq8p+3IxqT+CZbPXlvnfpjXhvF9hkekxuDGYSlIjbcNAuwDSqGD296M6NoHX7zHhMlsyt8vUWHTccE3VRJwPI+rBre+95JLt49Pi8fsIck2nwe5rPV5Q2USzMnsISigsl9HOHNQN1w9pBtvx0KoTOJuEoINvqnP5gA1MlSKsEvBdHyEHNdd0R2DU6Is59MzJkxQ54YzY/rEeXQ8BWgdQLVKi4UbDkAmEePjJWMpUQMJqOhwGd7LGoVHZvTHvZ8cxj9+OAmt3QJdQjqzMJnYZjRhWI9oXHdFd4RIxU4b61MzEjHsUpBgPVok4mlexCvkSIoMtfRCm/GtRwuVSTCqdxyusWpwi8Ui9IwJs/mtEBoMmQ3vGYNZg503oCPkUqfBqFCT+idgXN94DOnmutfZvsFtNrpPnCW4MLMexYwX2MBVyKVuF+wDjoGZeWTE/D7x1VeoTMLbKO0VF47Y8BAkRYUKnsaZHB3qUE6RSGSZYmvNHIjwjaSGh0hsRn7NwZ59495cl/bPbn3tuwsIxqfFYzxPh0F4iNQnaeEZGGQSsU0HBt8UyYxLgXmf+AgkRsrR3cuRjigXI9NSscjy2eZbyymxuo+vg2NqRiImpSdYrumxfeMso+09Y8NszsscNI9MjbXJLBkdJhN0PU3qbzsadZ3VCH2Mk/dUIhZhNE+QYY7bhMwgCJVJeDslnK1lA1o/n+b6zOgWaemosP9eEBL8RYfJHL4zAOCm4Sk25c9Mvfzc5uu0fxtGVwGrTi2BESUFaEHu1MUmzF2/D8N6xmBD1ihLdE9IIIlEItw2phe2PzIZxy804sY39+JoaX2gi0WIX0jFIptecLFYZNNw7hkbjn6JCgxIjnSY+mZteM8YhyDCVfY2T6bP9O8WabMPWWyEzKYRxPdcPazWBUklIreZ5CLkUswekixgJMRx3EIhl1nKobj0u5YUGWopo/XaGfsRCaD1fPhYr3uLV8htGp58pg1IsrxH1vXDN0JiHVDMyeyBOLug11l9OWu4TslI9DhIsV9TFRUqRVqCwmF9nfm4MX3ibEbL5mT2cDnix8dhqmqo1GHqov1oprlOo8NlDqnYrxma7DCSOcxJEO4r5jMwfyaH9YzBuL7xDoGUOXiIiwhxuC8qTOZQF6P7xFmmzVkH8NHhMpt6nznIdgRoRK8Yh/T+MeEhiFfIrUZHRZbR45GpsTYzp8xTHOMVco/2GwuTSZCgkHu9nYarANVTE71IyjOoe5Slc8EcsJqnd9qfU4+YMJtRSPNngi+piKsAUSjz1MxxfeNdTgE3d8q5WyfcMfJ4dlE7C6rxyH+PYvmsDNw1oY9PLiBCfKlXXDi+WHolPt5fjMUbD+LmkT3wxOwBNMpLOr3wECl0Rj3vfUIX0/P19Av9lpdJxDCYhKUJD5VJMCg5Cr3jI3CkpLUjhS8wEJrRzLrh6W1Dj29p0vh+8SisVKGgssnt7505YYL9YaP7xEFvVS/uRrNucccAABdOSURBVMbsRymB1qCGYww5JbadTpPSE/DDiYuWv5OjQlFnl+xiQHIkCisdU95PzUhESW0LimvV6BYVil4uEheY8U3Hm5ie4JCCXiwWWRqp5uqIkEstIynmwCFS7v33creoUHCMQanSgTGgfzcFmnVGp2VNjJTjBieJX/hGl+LCbRvNQ1KiUVipgpFzfo3HhIdYEp/Ym5SegL1OtiUwS44OxTXRtqMuaQkKJEWFQhEiBQPw/aVU7elJCt6EJD1iwtAjJgwjePZHsx6ltP+8SSViZHSLhFKlQ02zzmU5+ZjXbnnqagGjTK46DSJDZRjVOxY5JfWWID05unXLEKG6KaSIiQjhXarjar8/Z9KTFEhPcpy6aD3ad9Wgbh5NY7buCOmbEOF0u5DJ/ROx54wSwOXv/eToUJffzWJRazIcd+WhAC0IMcbwXvY5vLWzCP9vfqbP5mkT0h7EYhHuntgXs4ckY+W2fFy1djdW3jgE112RTJ0KpNPy5tKOCw/h/aEf1D0KydGhSEuIsKzZcGdqRiJMArOpRoRIIRaLoJBLMcVFJkMhrkyLt2l4WvMkCYIz6UkKJChC0KQ1OtwnErUuzLdu/FinWc/oFokQqdhp+VqfQ2R5XEMLf4AdGSq1CT6A1ml99usAw0Jap5cKERMeAolYBIlYxJuUYU5mD3ybW+72eewDYqGjYc4a8/bry6ZlJNlcVyNTYxEVKkN0uAyHiusstwtJpuFJg9heepICeiOHM9XO93dzlZXT27VCYrHIoYNRIm5NihMMv2d9EyKg92L/NiESFHLUNOss68V+OH4RRo5zSFjSMzbcJpFHgkKOGQP526l8NZYeL8cgJ1sNMFhluAX/ui1XWaRFIhHv/Z6MGtt/VmQSsU0GzTF94izTzb3N/Crks0EBWpBp1BjwxJfHcFbZjK0PTGhTRiFC/CklJgzvLx6Nn/Ir8dK2fHyyvxjP3zBYcIYoQjoSb3aaCAuR8DZkPPmen5iegDq13qMGR1s2101UyAFRa48yY67333I2WjUgOcphOqAzErEI8Qo5b4AGgDe1tXlEz1VgZs9Z+yg1LhwRIVK06FvX1ZoDJ6Ftc1ezByJDZS6/D8NDpE7LxefGYSk26xnH9o1DbHgIzhdVCX4O+zTs9mvKrKfPJSjkqGjQ8I7sxYbL0DM2DAfP19kEzfa82f4hLiKk3YISd/olKhChkwdFcAa0bSqo/VrA6DCZzWbd4/rG2YwQzxyUBMacJygSwtP9whhjlqBnYnoC/yibD9YutkVbEyQJRQFaEDlZ0YQHPs/B0JRobHtokk8W0BLib7OHJGNK/0S8v+ccFmw4gOuv6I7HZ2c4rEMghHguQSF3msV3VO9Y2M8Imz4wCaFt2Hx6ggfrRJy1YRXyywkhpg9Mws6CarfPxZdIwNk0I08CMzO+5CwALFPVkiLlLtcPOpMSE4abhvNP7XNn+oBEjwIB+60JrHv520PfhAgcv9DgEEQKnWoXEx7idNSLucivOCk9oV02AxdiaI9oSJpsg1apWMy7jskb0WEyr6Y4usPXWWL/G+xuw2xvpy9bG9Q9kncNqT37q37agCTe6cfmcvlyf0R/MAed49PiBY+6UQQQBEwcw4d7z+Pfv53BY7MysGQirTcjHVtYiASPzOyP20b3wus/FWD6/+3C4gl9cN+UNJe9q4QQ7/HtHyRkPahCbjulr1+iwqsOFSE97ULXp/aMDbdZX3TT8BSf/i6O7hNrGSXjIxKJeLNmAu73A/O2nPb152oql694WtLJ/RO92ujYG+ZqFIlEXk0pBtqeVp/PdVcI33vLnSEpUdAZTbhQr/HZcwKtnRbughj7aXbm69aXzU+pRIzoMM87UJwFZ8EoVCZxyGRtX4fmYDjJxQwEexSgBVhpbQue+PIYGjUG/HfZlTQdjHQqydGh+OdtmSiY0oR//nwak1/fiWWT03DXxD4288wJ6WgC1KHfLuIiQmwCNG9+h2YPSRbc4y4Ri1qnTrphPZXJ152WoTKJVyMEniYbaAt/TOXydApaW9bcuNr2gG9Ek2/7gGDgy2uxNfh0/nxtmZ7sKYlYhBkDk3iTuPiLkNE2V4b1iLaZtukPI1Nj8cdZ24Q0PWLCcLKiCRqDibfjTAgK0ALEaOLw+Z+leOOnQtx+ZW88Nqt/QD8UhLSngclR2LB4NI6VNWDtL6fx/p5zyBrfG3dP7BvoohHS5Q3vGWPZDNlbngQ7zrL7BSuZRGxZQ+fP9S/9kxTo3cYGqytXDermt/Nx955HhUkd9tjql6hwOorpqWDuEIyPCEFZXQvvfXKppM2jx/ERcsEJhQJdT6kebBnAhy8zbnvje2vMI/CldS2Cs/raowAtAA6er8PKbflgjOGjJWMEZ4EipKMb3isGnywZixMXGvHu7rOY8vpOTO8bgcfim3nT5BJC2p9YLIJc3DU6CL2ZOnXNkGSfTvsSSiQStWvHbaCTLZiJRK3nar+XmvX2AW2VGCkP2nVLveMj0DveeWDR1hG7ieltz67a3jryqh5v1sAKERyfzi6ivEGD138swK5CJZbPysDt41LblB2HkI7qip7ReOv2kTinbMYb3x3BnPV7kZkag8Xj+2DmwCT6XJCg1y1K7rDuoD0MTI5EgoLWbfpKXESIxw11+2QcxLe8yexIhKOcBo7kUglCpL6pl6hQ203JzXrEhKFRY/D6eSlA84OKBg3e2lmEr4+UY+6IHvj98ale79FBSGeSlqjAw+MTsXpROrbmXMCaHQVY+W0+5o7ogVtH9UB6Em0zQYLTwOQowRs7t0WoTOL1GgZCgt2k9ISAT6sjgScSiTC+XzzCvdio2hszByV5nCDHFb5RtKSoUI+SgtijAK0dnVM2Y9O+YnyVcwE3Du+On/42pc0LIAnpjKJCZbh7Yl/cOb4PDhXXYeuRC5j71h/olxiBG4enYPaQZJv9eAghhAQnmQcjjtRZTcz8uRWPsz0bgwkFaD7GcQy7Tlfj4z9KcKi4DnMyU7Dj0ckBWbhISEcjFoswLi0e49Li8dJNQ/HzyUpsP3ERb/xciPQkBWYPTsa0AUkYkhJF044IISTITBuQhFBZ8Dd+CQl2FKD5SGGlCv/LLce23AqIRMDi8b3x7wWZtOcTIV4KC2ndjHJOZg+06I3YXajET/mV+OiPYjC0To2ZlJ6Akb1jkJagoICNEEICrCPtX0VIMKMAzUscx3CivBE7C6vxY14lyus1uPaKZLx+6zBcmRbvt31SCOkKwkOkuPaK7rj2iu7gOIaTF5uQfUaJ745X4JUfTkIEIDM1FiNTYzAyNRaDU6KQQFNnCHEpUSGH3uTfPYMI8SdqiZGOigI0gRhjKKltwcHiOhw4V4vs00oAwNSMJPztqv6XhvW7RppiQgJJLBZhaI9oDO0RjQenpYPjGM7VNONISQOOltXj1e2nUFTdjOgwGfp3U2BAt0hkJEciPVGB3vERSIqU02gbIQAmpCeACdwfiZCOZlJ6QtBsJUCIp+jK5cEYg1Klw8ELLfi14gzyK5qQU1oPtc6IkamxGNs3Dnfe2QdX9Iimhh4hASYWi5CeFIn0pEjcNqYXAEBv5HC+Ro3CKhXOVKmwu1CJjXvP40K9BgDQMzYMvWLDkRrX+l+vuDB0u5RxKVEhb7d9TQgJNpSCm3RWlICEdGRdNkDT6E1QqnSoUmlRWtuCklo1zl/6f3GNGhqDCanRMoxKk2J8v3g8OL0fBnePov2ZCOkAQqRiDEiOxIBk2zT9HMdQrdKhrL4FpbUtKK1rwanKJvx8shJVTTpUq7TQGjjERYQgKVKOpKhQdIuUIylKjqTIUMQrQhATFoKYcBmiw2SIjQhBRIiEGrmEEEII8ZkOFaAxxqAzcpf+M0FnaP2/1nDpNoMJzTojmrRGqLQGNGku/V9rgEprRJ1aD6VKB6VKB5XOCIVcisRIOVLjwtEnPhwjesXg5hEp6B0fgV6x4Th7phCDBg0K9GkTQnxELBYhOToUydGhGNMnzuF+xhiatEYoVVpUN7V24FQ36VDVpMOh4jrUqfVoaDGgUWNAQ4sear0JUrEIMeEyRIXKECGXIjxEggi5FBFyKQwtKqQUMUTIJZb75FIxQqTi1o0yJWLIZeJL/2/9u/U+838ShFw6nta1EkIIIV1DwAI0g4nDnPX70KQ1gOMYOAaYGANjDKZLf3OMXb6PYzaLmWUSEeRSyeWGjKz13wq5FFFhMkSGShEZKkVUqAx9ExSICpMiJizkUk+4HImRcr9tiEcI6RhEIhGiw1pHx4Rskq0zmtCoMaDxUtCm1pvQojNCrTdBrTOi+IIW4SESNOuMqGrSQq03QW/koL/UyaS/1OGkN3LQmzjoDOb/m6A3cTCYLq8PkopFCJGKIZO0BmtikQgSMSARiSCRiCARiSAWt/7ffL9U0vp/X8d2o5PEoL4rQgghpH2ImBcrhHNzcyGX09xeQgghbaPT6ZCZmRnoYniEfgMJIYS0lavfP68CNEIIIYQQQgghvkcZLwghhBBCCCEkSFCARgghhBBCCCFBggI0QgghhBBCCAkSFKARQgghhBBCSJCgAI0QQgghhBBCgoTbAI3jOLzwwguYP38+srKyUFJSYnP/li1bcMstt+C2227Dzp07AQAtLS146qmnsGjRIsybNw/Hjx9vn9K3gTfnVVFRgTvuuAO33347HnzwQWg0mkAU3S135wYAdXV1mD17NnQ6HQBAq9Xi4YcfxqJFi7B06VLU1dX5u9hueXNeKpUK999/P+644w7Mnz8fR48e9Xex3fLmvMzOnj2LUaNGOdweLLw5N5PJhFWrVmHBggW45ZZbLJ+/YOLttXjvvfdi0aJFuOuuu6BUKv1dbEHcndtHH32EefPmYd68eVi/fj2AjvH90VZC3nNi6+abb0ZWVhaysrKwYsUK5ObmYt68eViwYIHl2nFWr3zHdlXHjh1DVlYWAKCkpAQLFy7EokWLsHLlSnBc6/6w69evx6233ooFCxZY2lyeHNuVWNfnyZMnMXnyZMt1un37dgBUn0IYDAY8+eSTWLRoEW699Vb89ttvdH36EnPjp59+Yk8//TRjjLGjR4+y+++/33JfdXU1u+GGG5hOp2NNTU2Wf69bt45t2LCBMcbYqVOn2DfffOPuZfzOm/P6xz/+wT777DPGGGP//Oc/2SeffBKQsrvj6twYYyw7O5vNmTOHjRgxgmm1WsYYYx9++CFbt24dY4yx77//nr3yyiv+LbQA3pzXv//9b7Zp0ybGGGNnz55lc+fO9WuZhfDmvBhjTKVSsaVLl7Irr7zS5vZg4s25bd26la1cuZIxxlhlZaXl/Qsm3pzXRx99xNasWcMYY2zz5s3stdde82+hBXJ1bqWlpezmm29mRqORcRzH5s+fz06dOtUhvj/ayt17TmxptVo2Z84cm9tuuukmVlJSwjiOY/feey/Lz893Wq98x3ZFGzZsYDfccAObN28eY4yx++67jx04cIAxxtjzzz/Pfv75Z5aXl8eysrIYx3GsvLyc3XLLLR4f21XY1+eWLVvYxo0bbY6h+hTmq6++YqtWrWKMMVZfX8+mTp1K16cPuR1By8nJweTJkwEAmZmZyMvLs9x3/PhxjBgxAiEhIYiMjERqaioKCgqwd+9eyGQy3HPPPXj77bctjw8m3pzXoEGD0NTUBABobm6GVCoNSNndcXVuACAWi7Fp0ybExMTwPmbKlCnYv3+//woskDfnddddd2HBggUAWkdmgnFzWW/OizGG559/HsuXL0dYWJhfy+sJb85t79696NatG5YtW4bnnnsOM2bM8GuZhfDmvDIyMqBWqwF03O+P5ORkfPDBB5BIJBCJRDAajZDL5R3i+6Ot3L3nxFZBQQE0Gg2WLFmCxYsX49ChQ9Dr9UhNTYVIJMKkSZPwxx9/8NZrc3Mz77FdUWpqKt58803L3/n5+Rg7diyA1s+auQ4nTZoEkUiElJQUmEwm1NXVeXRsV2Ffn3l5edi1axduv/12PPPMM2hubqb6FOiaa67Bo48+CqC1TSKRSOj69CG3AVpzczMUCoXlb4lEAqPRaLkvMjLScl9ERASam5tRX1+PpqYmbNy4ETNmzMCaNWvaoeht4815JScn4/PPP8f111+P7OxsXHPNNX4vtxCuzg0AJk6ciNjYWIfHmM85IiICKpXKP4X1gDfnFRUVhdDQUCiVSjz55JNYvny538orlDfntX79ekydOhUDBw70Wzm94c251dfXo7S0FO+99x6WLl2KFStW+K28QnlzXrGxsdi3bx+uu+46bNy4EbfeeqvfyusJV+cmk8kQFxcHxhjWrFmDwYMHo2/fvh3i+6Ot3L3nxFZoaCjuuecebNy4ES+99BJWrFhh05lkvk746tX+ts56TQkxe/Zsm84cxhhEIhEA53Vovt2TY7sK+/ocNmwYnnrqKXz++efo1asX3nrrLapPgSIiIqBQKNDc3IxHHnkEf/vb3+j69CG3AZpCobD0+gKt88XNF7f9fWq1GpGRkYiJibH0ek+fPj0oexq9Oa/XX38dr732Gn744Qc8++yzePrpp/1ebiFcnZuQx6jVakRFRbVrGb3hzXkBQGFhIe666y489thjlt6aYOLNeW3btg1bt25FVlYWlEollixZ0t7F9Io35xYTE4Np06ZBJBJh7NixKC4ubudSes6b81q/fj3uvfdebN++HRs3bsTDDz/c3sX0irtz0+l0eOKJJ6BWq7Fy5UqHxwTr90dbefv901X17dsXN910E0QiEfr27YvIyEg0NDRY7jdfJ3z1yvcb3BmvKW+IxZebbc7q0Nxm8eTYrmrWrFkYOnSo5d8nT56k+vTAxYsXsXjxYsyZMwc33ngjXZ8+5DZAGzlyJLKzswG0LtrNyMiw3Dds2DDk5ORAp9NBpVLh7NmzyMjIwKhRo7B7924AwKFDh5Cent5OxfeeN+cVFRVluVCSkpIs0x2Djatzc/UY83uWnZ2NUaNGtWsZveHNeRUVFeHRRx/F2rVrMXXq1PYuole8Oa9ffvkFn376KT799FMkJibiww8/bO9iesWbc7P+/igoKED37t3btYze8Oa8rL8/4uPjbX6Egomrc2OM4cEHH8SAAQPw8ssvQyKRWB4T7N8fbeXNe96VffXVV1i9ejUAoKqqChqNBuHh4SgtLQVjDHv37sXo0aN561WhUEAmkzkcS4DBgwfjzz//BND6WTPX4d69e8FxHCoqKsBxHOLi4jw6tqu65557LIko9u/fjyFDhlB9ClRTU4MlS5bgySeftMwIoevTd0SMMebqAI7j8OKLL+L06dNgjOHVV19FdnY2UlNTMXPmTGzZsgWbN28GYwz33XcfZs+ejYaGBjz33HNQKpWQSqVYs2YNevbs6a9zEsSb8yoqKsLLL78MjuPAGMOzzz6LwYMHB/pUHLg7N7MZM2Zgx44dkMvl0Gg0ePrpp6FUKiGTybB27VokJiYG8CwceXNeDzzwAAoLC9GjRw8Arb3g77zzTqBOgZc352XN2e3BwJtz0+v1WLlyJc6ePQvGGF588UUMGTIkgGfhyJvzqqqqwnPPPYeWlhYYjUY88sgjmDhxYgDPgp+rc+M4DsuXL0dmZqbl+OXLl2PgwIFB//3RVnz10q9fv0AXK2jp9XqsWLECFRUVEIlEeOKJJyAWi/Hqq6/CZDJh0qRJeOyxx5zWa25ursOxXdWFCxewfPlybNmyBefPn8fzzz8Pg8GAtLQ0rFq1ChKJBG+++Says7PBcRxWrFiB0aNHe3RsV2Jdn/n5+XjllVcgk8mQkJCAV155BQqFgupTgFWrVmHHjh1IS0uz3Pbss89i1apVdH36gNsAjRBCCCGEEEKIf9BG1YQQQgghhBASJChAI4QQQgghhJAgQQEaIYQQQgghhAQJCtAIIYQQQgghJEhQgEYIIYQQQgghQYICNEIIIYQQQggJEhSgEUIIIYQQQkiQoACNEEIIIYQQQoLE/wfb5H6QxB6PRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot({\"p_post_est\": brute_trace[\"p\"], \"p_post_conj\": conj_samples})\n",
    "plt.savefig(\"beta-binomial-samples.png\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}

beta-binomial-conjugate-example.py

"""Symbolic-PyMC Beta-Binomial Conjugate Example

Using the extremely naive model from
https://github.com/zachwill/covid-19/blob/master/covid-19.ipynb as an example,
we'll show how auto-conjugation could save one from needlessly sampling a
posterior known in closed form.
"""
import pandas as pd

import pymc3 as pm

import theano
import theano.tensor as tt

import matplotlib.pyplot as plt
import seaborn as sns

from operator import add

from unification import var
from etuples import etuple


from kanren import run
from kanren.core import eq, lall

from symbolic_pymc.theano.meta import mt
from symbolic_pymc.theano.pymc3 import model_graph, graph_model
from symbolic_pymc.theano.utils import canonicalize

sns.set_style("whitegrid")

theano.config.cxx = ""
theano.config.mode = "FAST_COMPILE"
tt.config.compute_test_value = "ignore"

df = pd.DataFrame(
    [
        ["Singapore", 72, 0],
        ["Italy", 46, 21],
        ["Japan", 32, 5],
        ["Hong Kong", 30, 2],
        ["Thailand", 28, 0],
        ["South Korea", 24, 16],
        ["Malayasia", 20, 0],
        ["Vietnam", 16, 0],
        ["Germany", 16, 0],
        ["Australia", 15, 0],
        ["France", 11, 2],
        ["UK", 8, 0],
        ["USA", 7, 0],
        ["Macau", 6, 0],
        ["Taiwan", 5, 1],
        ["UAE", 5, 0],
        ["Canada", 3, 0],
        ["Spain", 2, 0],
    ],
    columns=["country", "recovered", "deaths"],
)

df["combined"] = df.recovered + df.deaths


with pm.Model() as model:
    # No idea why one would use such a specific, "informative" prior. Why would
    # we want to put that much density on 1/2?  Why?!
    # For that matter, I'm not sure why one would use this model with this data
    # at all.
    p = pm.Beta("p", alpha=2, beta=2)
    y = pm.Binomial("y", n=df.combined.values, p=p, observed=df.deaths.values)


# Convert the PyMC3 graph into a symbolic-pymc graph
fgraph = model_graph(model)

# Perform a set of standard algebraic simplifications
fgraph = canonicalize(fgraph, in_place=False)


def betabin_conjugateo(x, y):
    """Replace an observed Beta-Binomial model with an unobserved posterior Beta-Binomial model."""
    obs_lv = var()

    beta_size, beta_rng, beta_name_lv = var(), var(), var()
    alpha_lv, beta_lv = var(), var()
    beta_rv_lv = mt.BetaRV(alpha_lv, beta_lv, size=beta_size, rng=beta_rng, name=beta_name_lv)

    binom_size, binom_rng, binom_name_lv = var(), var(), var()
    N_lv = var()
    binom_lv = mt.BinomialRV(N_lv, beta_rv_lv, size=binom_size, rng=binom_rng, name=binom_name_lv)

    obs_sum = etuple(mt.sum, obs_lv)
    alpha_new = etuple(mt.add, alpha_lv, obs_sum)
    beta_new = etuple(mt.add, beta_lv, etuple(mt.sub, etuple(mt.sum, N_lv), obs_sum))

    beta_post_rv_lv = etuple(
        mt.BetaRV, alpha_new, beta_new, beta_size, beta_rng, name=etuple(add, beta_name_lv, "_post")
    )
    binom_new_lv = etuple(
        mt.BinomialRV,
        N_lv,
        beta_post_rv_lv,
        binom_size,
        binom_rng,
        name=etuple(add, binom_name_lv, "_post"),
    )

    # TODO: We could also transform non-observed conjugates.
    return lall(eq(x, mt.observed(obs_lv, binom_lv)), eq(y, binom_new_lv))


# TODO: We could walk the entire graph and find all occurrences of Beta-Binomial
# conjugates, but, since we're only looking for observed Beta-Binomials,
# they'll always be the function graph outputs.
q = var()
res = run(1, q, betabin_conjugateo(fgraph.outputs[0], q))
expr_graph = res[0].eval_obj
fgraph_conj = expr_graph.reify()

# Convert the symbolic-pymc graph into a PyMC3 model
model_conjugated = graph_model(fgraph_conj)

# Do some wasteful MCMC sampling using the original model
with model:
    brute_trace = pm.sample(draws=6000, tune=2000)
    # Zzzz

# For our conjugate model, we can just sample the posterior term directly
conj_samples = model_conjugated.p_post.random(size=len(brute_trace["p"]))

# Wow, that was quick!

# Now, let's compare the two...
pm.traceplot({"p_post_est": brute_trace["p"], "p_post_conj": conj_samples})
plt.savefig("beta-binomial-samples.png")

# What do ya know; they're essentially the same!