EECS491 - A2 - E1 - tdm47.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EECS 491 - A2 - E1\n",
    "### Tristan Maidment - tdm47\n",
    "\n",
    "## Goal\n",
    "The goal for this exercise is to use the TCP data acquired from the wireshark dump made in A1.  By using the same model, we can create a simple example with real data that can be solved via sampling and variable elimination.\n",
    "\n",
    "## Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "import csv\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import math\n",
    "from pymc3.math import switch\n",
    "from pymc3 import Bernoulli\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data from A1 is imported, and a new parameter, `hops_retrans_all`, is used to determine all retransmission events, including ones that are not Fast Retransmits and Timeouts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Packets successful 104458\n",
      "Packets unsuccessful 204\n",
      "Packets timed_out 73\n",
      "Packets fast_retransmitted 43\n",
      "Percent unsuccessful 0.0019491314899390418\n",
      "Percent timed_out 0.0006974833272821081\n",
      "Percent fast_retransmitted 0.0004108463434675431\n",
      "Packets total 104662\n"
     ]
    }
   ],
   "source": [
    "hops = []\n",
    "hops_retrans = []\n",
    "hops_timeout = []\n",
    "hops_retrans_all = []\n",
    "with open('TTL3.csv', newline='') as csvfile:\n",
    "    ttls = csv.reader(csvfile, delimiter=',', quotechar='\"')\n",
    "    for row in ttls:\n",
    "        if row[1]:\n",
    "            ttl = 64 - int(row[1])\n",
    "            #I removed all TTL that are below 0.  \n",
    "            #These TTL are either set at a value higher than 64 by the sender,\n",
    "            #or packets being sent to local devices on my home network. (TTL = 0)\n",
    "            if ttl > 0:\n",
    "                if \"[TCP Retransmission]\" in row[2]:\n",
    "                    hops_timeout.append(ttl)\n",
    "                elif \"[TCP Fast Retransmission]\" in row[2]:\n",
    "                    hops_retrans.append(ttl)\n",
    "                elif \"retransmission\" in row[2]:\n",
    "                    hops_retrans_all.append(ttl)\n",
    "                else:\n",
    "                    hops.append(ttl)\n",
    "pkt_total = len(hops) + len(hops_timeout) + len(hops_retrans) + len(hops_retrans_all)\n",
    "pkt_succ = len(hops)\n",
    "pkt_time = len(hops_timeout)\n",
    "pkt_fast = len(hops_retrans)\n",
    "print(\"Packets successful\", pkt_succ)\n",
    "print(\"Packets unsuccessful\", pkt_total - pkt_succ)\n",
    "print(\"Packets timed_out\", pkt_time)\n",
    "print(\"Packets fast_retransmitted\", pkt_fast)\n",
    "print(\"Percent unsuccessful\", (pkt_total - pkt_succ)/pkt_total)\n",
    "print(\"Percent timed_out\", pkt_time/pkt_total)\n",
    "print(\"Percent fast_retransmitted\", pkt_fast/pkt_total)\n",
    "print(\"Packets total\", pkt_total)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "P(R) & = 0.0019491314899390418 \\\\\n",
    "P(T | R) & = 0.357843137254902 \\\\ \n",
    "P(F | R) & = 0.21078431372549 \\\\\n",
    "P(T | \\bar R) & = 0.0111083296707496513 \\\\ \n",
    "P(F | \\bar R) & = 0.0236974833272821081 \\\\\n",
    "P(T) &= 0.0006974833272821081 \\\\\n",
    "P(F) &= 0.0004108463434675431 \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Due to inaccuracies present in Wireshark, and the present of other modern TCP optimizations, we will assume that Wiresharks representation of P(T) and P(F) are independent.  For the purpose of this exercise, we will assume that we estimated the following:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "P(T = true| R = false, F = false) & = 0.3212359034102345301 \\\\\n",
    "P(T = true| R = false, F = true) & = 0.0111083296707496513 \\\\\n",
    "P(T = true| R = true, F = false) & = 0.1767348075841529761 \\\\\n",
    "P(T = true| R = true, F = true) & = 0.0236974833272821081 \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this data we can define the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    retransmit = Bernoulli('Retransmit', 0.0019491314899390418)\n",
    "    fast_retransmit = Bernoulli('Fast Retransmit', switch(retransmit, 0.21078431372549, 0.0236974833272821081))\n",
    "    timeout = Bernoulli('Timeout', switch(retransmit, 0.357843137254902, 0.0111083296707496513))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using that model, we can sample the network.  Due to the very low probabilities, a large sampling value must be used for accurate results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "BinaryGibbsMetropolis: [Retransmit, Fast Retransmit, Timeout]\n",
      "100%|██████████| 40500/40500 [00:15<00:00, 2655.16it/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGoCAYAAABbkkSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XmcZFV9///Xp6v37pmeFYbZBxjBAVlkREVN0JgvICDGJaJJXGLCN1+XxPjz65qvEo3RmEXjkhBERYgREDckKoLgFgScAQSGxRlmBph96+l9qeXz++Oe7qnpqaqu7q6qW1X9fj4e9eh7b93lc8691feee84919wdERERERERmbmGuAMQERERERGpFypgiYiIiIiIlIgKWCIiIiIiIiWiApaIiIiIiEiJqIAlIiIiIiJSIipgiYiIiIiIlIgKWCIyKTN7iZk9EXccIiIiY3RukmqlApbIDJnZdjMbMrN+M9tjZteaWWcRy51vZjsqEeNMufsv3P2UsfGQ5pfHGZOIiOSnc5NIfFTAEimNS929EzgLOBv4YClWamaNpViPiIjMSjo3icRABSyREnL3PcBtRCczzKzFzP7JzJ42s71mdpWZtZlZB/BDYGm4u9hvZkvN7Eozu9nM/tPMeoG3mNm5ZvYrMztsZrvN7Atm1jy2TTNzM/sLM9sc5vmimVn47mQz+5mZ9ZjZATO7ccJybw/L9ZnZx83sJDO728x6zeymse1k39E0s+uBlcD3Q9zvq1D2iojINOjcJFJZKmCJlJCZLQcuAraESZ8CnkV0UjsZWAZ8xN0Hwny73L0zfHaFZS4DbgbmAV8H0sBfA4uAFwK/B7x9wqYvAZ4HnAH8IXBBmP5x4MfAfGA58PkJy10AnAO8AHgfcDXwx8AK4HTgDRPT6O5/AjxNuDPq7p8uMntERCQGOjeJVJYKWCKl8V0z6wOeAfYBHw136q4A/trdD7l7H/D3wOWTrOtX7v5dd8+4+5C7b3T3e9w95e7bgf8AfnfCMp9y98Pu/jRwF+EuJZAEVgFL3X3Y3X85YblPu3uvu28CHgF+7O5b3b2H6C7m2dPJDBERqQo6N4nEQAUskdJ4lbvPAc4HTiW6o7cYaAc2huYRh4EfhemFPJM9YmbPMrNbw0PKvUQnwkUTltmTNTwIjD3I/D7AgPvMbJOZ/emE5fZmDQ/lGJ/0gWgREalaOjeJxEAFLJEScvefAdcC/wQcIDoRnObu88KnKzxwDOD5VjNh/N+Bx4G17j4X+BDRiamYePa4+5+7+1LgfwP/ZmYnTylRxcUoIiJVSucmkcpSAUuk9D4L/D7wHOBLwGfM7DgAM1tmZmNt0PcCC82sa5L1zQF6gX4zOxX4P8UGYmavC23vAbqJTj6ZolOS317gxBKsR0REKkPnJpEKUQFLpMTcfT9wHfAR4P1EDxXfE5pQ3AGcEuZ7HPgGsDU001iaZ5XvBd4I9BGdFG/MM18uzwPuNbN+4Bbgr9x969RTdYxPAn8T4n5vCdYnIiJlpHOTSOWYu2pTRURERERESkE1WCIiIiIiIiWiApaIiIiIiEiJqIAlIiIiIiJSIipgiYiIiIiIlEhj3AFU2qJFi3z16tVxhyEiIgVs3LjxgLtP9uLTuqDzkohIbSj23DTrClirV69mw4YNcYchIiIFmNlTccdQKToviYjUhmLPTWVvImhmCTN7wMxuDeNrzOxeM9tiZjeaWXOY3hLGt4TvV2et44Nh+hNZL8LDzC4M07aY2QfKnRYREZm9zOwrZrbPzB7J872Z2efCOekhM3tupWMUEZH4VeIZrL8CHssa/wfgM+5+MtHbu98Wpr8N6A7TPxPmw8zWAZcDpwEXAv8WCm0J4IvARcA64A1hXhERkXK4lug8lM9FwNrwuQL49wrEJCIiVaasTQTNbDlwMfAJ4D1mZsDLiN78DfA14Eqik9BlYRjgZuALYf7LgBvcfQTYZmZbgHPDfFvG3vxtZjeEeR8tZ5pERGR2cvefZ7euyOEy4Dp3d+AeM5tnZie4++6KBLjrQR7+7Ra2znshZyzrYs2eHwHwUPsL6M208uIladj+C273c1l1/Hx2dg9x8oGfsOLU57GveTm/2nqQF43ezfZ93Zx5wVtpamoaX/VdT+zjuDkt7O8bYVFnC6cv64JMGjZ9B1acC/NWwubbofM4MsefwZ33buQsnmDR8y+HhgSP7uqlb8cmnt+5H9a9Mnf8B7bw0Mb/oeXMV3PKkjlHpg8egifvhFMuYuv9d3LowG565p/BS+c8Q8Opr4DmjvFZn/r51xkYGmbdaWfB8nMAcHdufWg3ZyztYNWe26MZ110GiSaOEtLz6/SzaF64kjNXzBv/6t6tB2loMAZH0yxob2ZhZzO/3n6IlsYEF56+BIBf/s9PObHpEEvPfTUcfJLUMxt5cMdh1h7XycPz/xeL+x/jlMa90LEYTvxdANIZ59aHdgHQ0tjA3LYmWhoTnLNq/lGhPbqrlz29Q7xs7UIe+vG1LOrqYOnyEyE9Cmtewn3bDrG7ZwiA05d1cdLizvFlNz7VzUgqzXknLTpqnb/efohMxnn+iQsB2NE9yManurmkeSOJxafA4lPG5/3hw7tZ1dXA4EPfY15bEyefexHMidK9r3eYX209yEWnn0Dzjruhbw+seD7MW5F7PwOkU/Dod2HlC6BrOQAjqTQ/emQP5520iAee7obDT7NycBO7WtZwbvseOluiS9a+tZdx5+P7SGRGuaRpY7S+014NDQ30j6T4yWN7Of+U4+hqa+J7D+4E4LKzlgFw5+N7WTK3jXVL5x4T0gNPd/NM9xDuzkX8D82JBjj9NXzvN7uOmu/sFfNZsaCN/354N89Z1sWqhdHx98jOHvb3j/DSU45jf98Idz95gLXHzeHpQwN0JNKs2vcTVp5zEXQex/1PdzM8mua8kxcdE8f48d7cCU1tMLAfTjwfOqJ5tx0Y4NFdvVx8xglHLfbIzh6e3N8PwO+0bWN+Y5L7H3+Sxc9+CStOPLIvSSfh0e/BqvNg7lIYPMSujf/NY10vZmHjCKmtP2fX0t/nguVpWvY+GC1z6iWw/Rcw3AOL1sIJZ0b7efsvYcGJ0LsLUsOw8oXQtQy23AHti2DJc8L/iMLHw48e2cPJx3VysH+EtDv7+0Z44YkLOW5uK/z2Npi7DJacHta7kLt6l3H83FbWLZ3LLzbvZ9mhezmxfRiWng27HoBTL+b+3cP587gCyv0M1meB9wFj/ykXAofdPRXGdwDLwvAy4BkAd0+ZWU+YfxlwT9Y6s5d5ZsL05+cKwsyuILqbyJLFC9hw69UzSJJI9Vl/yRVxhyAiWeexYOx8dUwBK/u8tHLlytJs/eAWRvoPwTx4+JmDrAnlhwM7n6Rv7slw4EkA0oOHeGx3AoC9h/tYsftBNrdFhYnt+7oBGOjvZd78heOr7h1K0juUBKBnKBkVsJLRBT17HokKWMM9MNxDcvHpdPY8wc7MIItGB6B1Lpv39bFs38PQtiB//LsfZDQ5wrY9vUcXsA5ti/727+XQgSgru7ofItUxj+a+PbDwpOh7d/b3DkTD3dvGC1ipjJNx5/Gn97KqOaxzpA/aJ8QyGi2b3vMI21l4VAFrT+/w+PDhwVH290eXTyOp9Pj01u7fsgtYCrDrAQZHo+/29gxzsHmE1p2PwaoF0QVzMDiaGh8eSWXY3zcCcEwBa/O+vhBjP6PpDLsO9bF0zpFDbaxwBdGFdnYBa0f3ILnsOjx01Phju6NtJIeHSOx5+KgC1mg6w45nnmEBcHgoCfufGC9gbd4XXdT3DidZ1LcnWmDvpsIFrNH+MN+j4wWs7oHo+Nq6v5+hZJoTeh6jP5Nk7shv2Z9uprOlM6Q12hdNI4dgrIycHoGGNvaEfHjm0CBdy7qO2WzfcIq+4b6cBaynDx3Jp/7hFAs6msEzx8y3aVcPS+e1ks44D+/sGS9gjRVuALaEPBnbb9a3l319I6w8sBk6j+OZQ7n3SZQR247k0Vg+HXxyvID10I7DORfL3v6u7U8wf2kXGc+w64l7jy5gjYRjad9jUQGrexu7eoZw28XA0G5agKZkH8NPPUZLa/R/gsGD0e8b4MDmqIC1//Fo/NDWI+ve92hUwBo6HH0WnHhkeoHjYSSVZtOunqOmbd7XHxWwRvqibS05fXy9venj6B1Osm7pXA4NjNK2fxcnrloQFa4ABvbzzKFE3u1VQtmaCJrZJcA+d99Yrm0Uy92vdvf17r5+flfn5AuIiIiUUfZ5afHiWdFZoojIrFHOGqwXAa80s1cArcBc4F+BeWbWGGqxlgM7w/w7gRXADjNrBLqAg1nTx2Qvk2+6iIhIpRU6X4mIyCxRthosd/+guy9399VEnVTc6e5/BNwFvDbM9mbge2H4ljBO+P7O0I79FuDy0MvgGqKHh+8Dfg2sDb0SNodt3FKu9IiIiEziFuBNoTfBFwA9FXv+ajZwjzuCKlD/eZArhZOmuoLHRrXsAZ9xmgsv7wXGKqHWf+5xvAfr/cANZvZ3wAPAl8P0LwPXh04sDhEVmHD3TWZ2E1HnFSngHe6eBjCzdwK3AQngK+6+qaIpERGRWcPMvgGcDywysx3ARwlPgbj7VcAPgFcAW4BB4K3xRApuVv6NFNxGBbZfpOqJ5FhWif1UpIpGUkS6p3d9XZlUTH2/TWX+Sh8T1XMM1pOKFLDc/afAT8PwVo70Apg9zzDwujzLf4KoJ8KJ039AdEITEREpK3d/wyTfO/COCoUjIiJVqhLvwRIREZG6lbuuodpa+HjREZU28mrLh9I4NlW13qQrp7pM1ORsuunOt1yl87EK9psKWCIiIlK9StaCqXRNoUzNqkqu2FZ3tZDz5Wz5WU3NSiU/FbBERERERERKRAUsERERKZP4m+pINZslTQ1l1lEBS0REpK5N/Yq1uC6ga7ipUq70qelVYXmyJ/fkSpSSyre/in1eb+qFwenmywzyc0pBHslTr+XfdxVQAUtERGYdM2szs1PijqOy8l8wlfJSqtqeT4rtmZXqygaREiviAJ/FNy1UwBIRkVnFzC4FHgR+FMbPMjO9qF5EREpCBSwREZltriR6H+NhAHd/EFgTZ0AiRdEDSnmUK1+q9xmxQnEUFeNU0hFDoot/rUJ1UgFLRERmm6S790yYVttn82PE2zSnzjKzbKqpAVVlW3Mdu7FjpkwjoEqlYbLNFJG6AgtX01Eh09UYdwAiIiIVtsnM3ggkzGwt8JfA3THHJCIidUI1WCIiMtu8CzgNGAG+AfQC7441olpWLW2mSqXe0lMhtd6kS7JpX86UarBERGRWcfdB4MPhM+tYhQoQ097KMfHlaTLlueadwmpnYtI2YtndXdchFUKnp6bybTpNFfOlr9Lpjj+fVcASEZFZxczuIscZ2N1fFkM4UiF6r099mO2PKBmTFx8Kfj/bM7BCVMASEZHZ5r1Zw63Aa4BUTLHUt5q6Yy9Q4aZ+OY4PNTWUelBUAcvMnuPuD5c7GBERkXJz940TJv2Pmd0XSzDVqu4LRvWevphV+viZsL0Zb75QF+hZX1YslTNI0PTrq3xGS892xdZg/ZuZtQDXAl/P0b2tiIhITTCzBVmjDcA5QFdM4VROoaZB07mOyrO+amuBVHw4pQ28yrKhBPKkKLaE1k4Oe7X9KLIdFVupu5Ov4nSXWVEFLHd/SejK9k+BjeFO31fd/fayRiciIlJ6GzlyezYFbAPeFmtEIiJSN4p+BsvdN5vZ3wAbgM8BZ5uZAR9y92+XK0AREZFScvc1cccgIlI71KR2qop9BusM4K3AxcDtwKXufr+ZLQV+BaiAJSIiVc3MXl3o+9lzs7ACF0vlaBJV98+FxauirdhybazQ9uPa9dM85qwCmTnzLKnu39OMfu5V8L+i2BqszwPXENVWDY1NdPddoVZLRESk2l1a4DtHNwulwqr50Zz45cuc/JlWH9lZH6mY7YotYF0MDLl7GsDMGoBWdx909+vLFp2IiEiJuPtb446hPsV/t7iUrM7SUzHKNpFxxRaw7gBeDvSH8Xbgx8B55QhKRESknMzsYuA0ovdgAeDuH4svouripWhiM+1VTFgwTzVP1F128Rsp5fX/5HUM9V4LMaNX3ZYmgpqs/qtMKbQkv9/p5G++7VZBk71KayhyvlZ3HytcEYbbyxOSiIhI+ZjZVcDrgXcRXQm/DlgVa1AVUYsXpDM3fmk3O5MvdabYck/+Qlae1yvMwkJQORVbwBows+eOjZjZOcBQgflFRESq1Xnu/iag293/Fngh8KyYY6pTumgrqQpcBFd2j+XYmg4ZqQPFNhF8N/BNM9tFVPRdQnT3T0REpNaM3SAcDL3hHgROiDEemQabYpWUKrBE8itJk8uabLZZHsW+aPjXZnYqcEqY9IS7J8sXloiISNncambzgH8E7ie6Z/6leEMqo7LVesR7MeVTrOo4au5ceaKLw4I8X9Oyqii6+oSx8lWDZR86JXnWqbitVmTZ7H08raSVMD9qvSKz6BcNA88DVodlnmtmuPt1ZYlKRESkTNz942HwW2Z2K9Fzxj1xxhS3Ul4kV8PldjUo5l1I1VE4qU75CnQyuekUTqb+DFYR+2cW37Qo9kXD1wMnAQ8C6TDZARWwRESkppjZQ8ANwI3u/iQwEnNIUlVq/d65SImpA4wpK7YGaz2wzitXFyoiIlIulxI9R3yTmWWAG4Gb3P3peMOqlGJO5TM93Rd4m5T7NO9s6xKknCpbm3bstgptf6Z73t2LqlEs1ZYrUXETXZLPYEMVvaSffb/dYnsRfISoYwsREZGa5u5Pufun3f0c4I3AGcC2mMMSkSIUKrxMrxBVZeohDVJ0DdYi4FEzu4+sphTu/sqyRCUiIlJGZraKqBbr9URN398Xb0Q1TI1bhNlYR1G/9E6smSu2gHVlOYMQERGpFDO7F2gCbgJe5+5bYw5JCqrUHf3itzNpTYlqIaTW6RiekWK7af9ZuNu31t3vMLN2IFHe0ERERMriTe7+RNxBVLN4718Xt3V3aqL2rAZCnLrJElXjiS5b9JXKlyltph4LUvEff0U9g2Vmfw7cDPxHmLQM+G65ghIRESmXWVu4iuOOdFVdaNfDhWQF8rOifR8cu7FqOmKq0Yw7I6no/4F6+M1NT7GdXLwDeBHQC+Dum4HjyhWUiIiIiIhILSq2gDXi7qNjI2bWiG4yiIjILGNmF5rZE2a2xcw+kOP7t5jZfjN7MHz+LI44j4qpXKfrGntG4+g3zegSZsry7O+qOAwm1ISVs+I0liNnJjV9U8mMrJ05vXSWLndq/c1QxRawfmZmHwLazOz3gW8C3y9fWCIiIuVhZu1m9v/M7EthfK2ZXVLEcgngi8BFwDrgDWa2LsesN7r7WeFzTUmDL5OquEiehZTvBdRV5lQ6LRUonNTV/im9YgtYHwD2Aw8D/xv4AfA3hRYwsxVmdpeZPWpmm8zsr8L0BWZ2u5ltDn/nh+lmZp8LdwUfMrPnZq3rzWH+zWb25qzp55jZw2GZz1ldvABBRETK7KtErxx5YRjfCfxdEcudC2xx962hVccNwGXlCbGyavxmcc1SvhdQV5lT6bRU5E3H5d9GDSuqgOXuGXf/kru/zt1fG4Yny9kU8P+5+zrgBcA7wp2+DwA/cfe1wE/COER3BNeGzxXAv0NUIAM+Cjyf6OT20bFCWZjnz7OWu7CY9IiIyKx2krt/GkgCuPsgxV2RLAOeyRrfEaZN9Jpwo/BmM1uRa0VmdoWZbTCzDfv3759i+FPjs/hB8/yUJ3Vlwv31ct5uj+XImUmCprls3L+QWq8zKbYXwW1mtnXip9Ay7r7b3e8Pw33AY0QnosuAr4XZvga8KgxfBlznkXuAeWZ2AnABcLu7H3L3buB24MLw3Vx3vycU9q7LWpeIiEg+o2bWRritbGYnEdVolcL3gdXufgbR+epruWZy96vdfb27r1+8eHGJNp3btJ7BKutDLNVw57vCMVRFmuvYJM9glTL7q+UZrHIvG/cRW+vPYBX7ouH1WcOtwOuABcVuxMxWA2cD9wLHu/vu8NUe4PgwnO/OYKHpO3JMz7X9K4hqxViyuOiwRUSkPn0U+BGwwsy+TtRL7luKWG4nkF0jtTxMG+fuB7NGrwE+PaNIK6TGbxaXVgUv7Oop30uebfWUOTX3DNbE5XOsr672T+kV20TwYNZnp7t/Fri4mGXNrBP4FvBud++dsF6nAoXk7DuF87s6y705ERGpYu5+O/BqokLVN4D17v7TIhb9NbDWzNaYWTNwOXBL9gyhdcWYVxK13hARkVmkqBqs7A4niApl64tZ1syaiApXX3f3b4fJe83sBHffHU5E+8L0fHcGdwLnT5j+0zB9eY75RUREjjHhXAYw1ppipZmtHGvWno+7p8zsncBtQAL4irtvMrOPARvc/RbgL83slUTPIR+iuJqxyqpADY1jU99MuCNu+RasqSZDM+3uOgZTqpCYYapy1H4UqhCJbddPc8O1UbdTwUyt9A6sgv8VxTYR/Oes4RSwHfjDQguEHv2+DDzm7v+S9dUtwJuBT4W/38ua/k4zu4GoQ4ueUAi7Dfj7rI4t/hfwQXc/ZGa9ZvYCoqaHbwI+X2R6RERk9vnnAt858LLJVuDuPyDqSTd72keyhj8IfHC6AUo51cZlrxRWaC/Wxx4unIrytcyLv1BST4oqYLn7S6ex7hcBfwI8bGYPhmkfIipY3WRmbwOe4khB7QfAK4AtwCDw1rDtQ2b2caKmGQAfc/dDYfjtwLVAG/DD8BERETnGNM9lMildmMl06diR+lRsE8H3FPp+Qg3V2LRfkr8Y/ns55nfgHXnW/xXgKzmmbwBOLxSbiIhINjNrJbpB92KiK7xfAFe5+3CsgUlulXqYfgrbqY+aktKILy+0F0rPJgyrADxdU+lF8HkceZj3UuA+YHM5ghIRESmj64A+jjQrfyNwPVEPuUJVPMJQpJoJdJap7f3iZYu/dOst9BuNLfdr5x9H2RVbwFoOPDe8zwozuxL4b3f/43IFJiIiUianu/u6rPG7zOzR2KKpmPx3/FUXUHplydNKXMDGfI2sa/Rym8aRmXOnFLGeWdyVe1HdtBO9q2o0a3yUI++vEhERqSX3hw6SADCz5wMbYoynRs3ei6fZqZT7W8eO1Ldia7CuA+4zs++E8VeR5+30IiIiVe4c4G4zezqMrwSeMLOHiR4JPiO+0MphOlUC5axGmGzdlailqXQ1iaplyssLjJU29z2WKraZbHOSZfPUMsWTzqztx7r1mSu2F8FPmNkPgZeESW919wfKF5aIiEjZXBh3AFK9VLcyPaV/bql+9oRPpalcCZrVzXxP1HrxJn7F1mABtAO97v5VM1tsZmvcfVu5AhMRESkHd38qvFtxBVnnwcleNCwiUnf00FtZFNtN+0eJehI8Bfgq0AT8J9G7rkRERGpGeLfiW4AnOXKrtqgXDdcDq8Td6WndhZ9smRxxx3RxOGnysmeokQvYKe2xMiSp8PbjysPpbbcm+nao6HFZG7+BUiq2BusPgLOB+wHcfZeZzSlbVCIiIuXzh8BJ7j466ZwiIvWsJkqDtafYXgRHw4uAHcDMOsoXkoiISFk9AsyLO4i6UfCFPNVz53pKz8FIZUzj+LAC+3E27OJZkMS6UGwN1k1m9h/APDP7c+BPgS+VLywREZGy+STwgJk9AoyMTXT3V8YXkuRXqUtKXbpOR6ECTzmpwFx6RxV3zSZOkSkothfBfzKz3wd6iZ7D+oi7317WyERERMrja8A/AA8DmZhjqUql7xGu9Gohxlmrimoup6Vc4U83X3IsV/D4jy37a3y/l9CkBSwzSwB3uPtLARWqRESk1g26++fiDkJk6sp/ARv3JbIKztVI+2SqJi1guXvazDJm1uXuPZUISkREpIx+YWafBG7h6CaC9d1Ne8FnV0rd3ErNt+ovB0qYIjXvi0+p8r6o9cze/VzsM1j9wMNmdjswMDbR3f+yLFGJiIiUz9nh7wuyptVvN+3TaZZUziZek647z/cljanSd+RVA1BWE44Nzzk+ey72s1M/7dcyxH3Ixr39GSq2gPXt8BEREalpocm7SB41fmUXk1p/7EqydqB25owVLGCZ2Up3f9rdv1apgERERMrNzC4GTgNax6a5+8fii0hmI7WUy6++LvGnsqN1UNSDyd6D9d2xATP7VpljERERKTszuwp4PfAuoquZ1wGrYg1K8pc2qrAUYpNeBFdfzJOJq7v1atl+adVTWmQ6JitgZR8hJ5YzEBERkQo5z93fBHS7+98CLwSeFXNMFVSZuoHK9AZX/fUc1R/hdFQ2VfH1Ol6fey9SwbTVdT7mNlkBy/MMi4iI1Kqh8HfQzJYCSeCEGOOpcYUuD3TpIIXo+CiHWVieqTqTdXJxppn1EtVktYVhwri7+9yyRiciIlJ6t5rZPOAfgfuJrvK+FG9IUn5qtjUTnqcJX6Vb9lmOoSNT6n8flz6/s1aYtXKPLvVLvbFZo2ABy90TlQpERESkEtz942HwW2Z2K9Baz+95nHY3zWVTmniqJVXuXmfPD5VCteydyFSjKV/0013zscsVX0s1zdciTEe1VJ1VQRyTNREUERGpC2b2PDNbkjX+JuAm4ONmtiC+yESKVIELx7gvTavg2lgm0k6ZMhWwRERktvgPYBTAzH4H+BRwHdADXB1jXBWSv5alMhUws6yWZ5YlV2pFqQ5MHeCFFPuiYRERkVqXcPdDYfj1wNXu/i2ipoIPxhiXALpgm020r6W+qQZLRERmi4SZjd1Y/D3gzqzvdMMx21SaBJW4+ZDlXV2FnhUpQ3MomwVNrMqZwsmzzwuMSWG5OguJPwcr85qH8tEJRUREZotvAD8zswNEXbX/AsDMTiZqJigiIjJjKmCJiMis4O6fMLOfEL3z6sfu4/fFG4B3xReZFFZ9zckmfWatBnsVjDviY7Zfg3k45tjQp5CWmNI9G2pZK0kFLBERmTXc/Z4c034bRyyxqafrKF0UxmPW5Hsdp7Oi+7CO8zEPPYMlIiIi5TFrLsRlWnIcHzpiZq7Wn1+qBypgiYiIiMi0VLpB25EWdDk6Z6jdVoVTUM5EZq17dmRm2aiAJSIiIiIiUiIqYImIiMhRytrAqFTNBmfQVjzIAAAgAElEQVS0mtKl8Jjk5ExfqbZXI02/Ym4aOnHzUw2nbOFPd8UzCWi6i05ruWo5PuOPQwUsERGRIpnZhWb2hJltMbMP5Pi+xcxuDN/fa2arKx+lyCyjZ/2mT00By0IFLBERkSKYWQL4InARsA54g5mtmzDb24Budz8Z+AzwD5WNsoACF1KlvsTyHNvyybZSbxd6NZacuLM/7u2XUlUnpVQZXU87rAzUTbuIiEhxzgW2uPtWADO7AbgMeDRrnsuAK8PwzcAXzMyy3rlVcv293QA0JtPh7wCQYdij8Yb0CI3JAYaGBrBkmsb0UJgnMpxMkxrqpTE5Oj5tdLB/fL1H1jlhm5lU2GaaVG/3+PZH+g6TSEfrGhjoxTONNCYHMDLRtvp7oOHYy4/G4WEAEunBo7bdMDRAQzJNZrD/qPlHUxkY7CczPu+RLB5ORjEBZNxpTA7QkBll2KIY0wN9ON1HrY90ksZkmgZP0pgcKJj+9NAIjclUlMa+bizrYrO/r5vG4WFGkxkAUmmnMTU4HhcwHlsy7cesezx/s/MmzDMwcGQ72evKuX8mLJtvnWPTU4O9NCbTjKTTNDbYeIxj8ybSQ+PjQ8ODpMeWC8fOcH9iPKaxYyKv9Gh0vCSHxucbGk7RmBwgPTQc7Yf0keMxlc6Mr3tkoCfEMzg+LT3Qh6eM4YFhGpMDJIeS9PceSeNg32FGE5Y3L7LzAyAZtpfu78m9f/q6j1lX9vjoYC+NyZHx+cfybmh4iHTvsctmaxjsp2E8HyMelgPIDPfR6M5gfw+pxJHjYWKcw6loHUbmqO1YaphEMo2no3U2hN9VQ2aEhsxIiHeYJMMMJ6J6mMyEmFK93STC/5OjpIaP+l+Q6u8N+7nw8ZArj9Nu9Pc2H1lX1nob00fyb2zZ4axYMoP9NCZ9fJ5srW0dNDY1542lVKyM//Or0rq1q/y6z3w47jBESmr9JVfEHYJISZnZRndfH3cc2czstcCF7v5nYfxPgOe7+zuz5nkkzLMjjD8Z5jmQb73r16/3DRs2TDuujf99De6ZaS8vIjJbrDjr9zh++UnTXr7Yc5NqsERERCrMzK4ArgBYuXLljNa19PSX4A6WHmVkaJA+2pjX3kzzwC4sk2KocznJtNPelKBh+BD9iS6aEg2k0hmaUn00tXSQaWimd3iUrsY0Q927aDvu6AuQwdEUjQ0NpDMZEg1Gc2MCgIbhbrx5Dt7QSMNoP97QiDe20jc0QqcPYO0LgKgGYmR0lDk2TKZ1Xu48yaQY7j9MQ8cCmhJHP8GQGDpIum0hqdEh0j27GelYxnzrJ9228Kj50oPdJFMZ2tra8Kb28el9w0namxtpGu3GMmnS7YtzxtAw3M2gtWMNjbQ0JcanjyTTZBzASTQYTYkE3YMjNCUamNPaFOXR0DCtPkRD+3zIpGkY7WVgYJC2jnaGGuaS8CRt6T68sY1Mc+f4unuHkgynMrQ3N9DYEKW7NWvbY/k3nErT2dLE6MGnSHQsoLGxCcukyDR3MpJMM5LKMJRMs7Cjmcas/BtNZUhnMrQ1H33JN5pKk85AW3O0Lfcon+Y1DOGNrXjiyF3+/pEkrY0Jug/sYUFHM4n2+dBwZLne4VG62pqx1FC0rzqW4DlqKSfm9dixk50Xna2NDI2m8UyK5lQ/ww0ddCVGaEgPR3nX0sXB/hFamhroSh7AUkOk5p04vo6eoSiWsfWlMxnmd7RE+ygcx82Nxz4hM5rKkMpkGE1lmN+UIjF0kNTcFQyOphgYSdPenGA0lWZuWzOJBhs/phINNr6PRlIZOloax+OY09LEwNg2R7vH8y3fPhmTGDpIprENAEsNkWmdDxZqkzLOwGhq/Lgbk0pH+797YJSVcxM0eIqB4SGa2+eRaDx63obhbjLNc4/sw4ED9Ce6aElA98F9jDTPZ+W8FhpGezB30m0LaRjtx9IjZJo68MZW8Mz4PrTUEHgm639BH55oxhMtOffzRAMjKVoaG0hnnIxHNYidrY00mEXramjCG1tpSA7g1sBgpolE2I9Do2kSmRHaR/aTmruCxHA36baFefN47vzj8sZRSjVfwDKzC4F/BRLANe7+qZhDEhGR+rQTWJE1vjxMyzXPDjNrBLqAgxNX5O5XA1dDVIM1k6BOWHVKnm9OndJ6lo0PnTWDaCqhOuJbNvksRVk61QVW59vf5bfypGfnnF7pvCg039I8w5UW57an7sgxtXptjGHUkZru5KLIB45FRERK4dfAWjNbY2bNwOXALRPmuQV4cxh+LXBnOZ+/EhGR6lPTBSyyHjh291Fg7IFjERGRknL3FPBO4DbgMeAmd99kZh8zs1eG2b4MLDSzLcB7gGO6chcRkfpW051cFPPAcZg+3tYdOB14pKKBztwiIO8D0lWsFuOuxZihNuNWzJVTi3Gf4u5z4g6iEsxsP/DUDFdTC/tYMZZOLcRZCzFCbcSpGEtnpnGucvfcD3FmqflnsIqR3dbdzDZUW89Uk6nFmKE2467FmKE241bMlVOLcZvZ9LvVqzHFnKwnUwv7WDGWTi3EWQsxQm3EqRhLp1Jx1noTwWIeOBYREREREamIWi9gFfPAsYiIiIiISEXUdBNBd0+Z2dgDxwngK+6+aZLFri5/ZCVXizFDbcZdizFDbcatmCunFuOuxZjjVAv5pRhLpxbirIUYoTbiVIylU5E4a7qTCxERERERkWpS600ERUREREREqoYKWCIiIiIiIiVSdwUsM1tgZreb2ebwd36e+dJm9mD43JI1fY2Z3WtmW8zsxtB5RlXEbWZnmdmvzGyTmT1kZq/P+u5aM9uWlaazyhjrhWb2RMijY16iaWYtIe+2hLxcnfXdB8P0J8zsgnLFOM2432Nmj4a8/YmZrcr6LufxUgUxv8XM9mfF9mdZ3705HE+bzezNVRTzZ7Li/a2ZHc76LpZ8Dtv+ipntM7Oc78mzyOdCuh4ys+dmfRdXXk8W8x+FWB82s7vN7Mys77aH6Q9Wskv0ImI+38x6so6Dj2R9V/DYmo3izpNcx1G+81klf0O5jrNSxmVm54R0bwnLWolivNLMdmYd/6/I+i7n+TPfMWAluJ4xsxVmdpdF58ZNZvZXYXq15WW+OKsmP82s1czuM7PfhBj/ttB6bRrXUvliL0GMOa8x49rfYT0JM3vAzG4N41WTjwC4e119gE8DHwjDHwD+Ic98/Xmm3wRcHoavAv5PtcQNPAtYG4aXAruBeWH8WuC1FYgzATwJnAg0A78B1k2Y5+3AVWH4cuDGMLwuzN8CrAnrSVQof4uJ+6VAexj+P2NxFzpeqiDmtwBfyLHsAmBr+Ds/DM+vhpgnzP8uos5pYsvnrG3/DvBc4JE8378C+CFgwAuAe+PM6yJjPm8sFuCisZjD+HZgURXm8/nArTM9tmbDpxryJNdxRJ7zWSV/Q7mOs1LGBdwX5rWw7EUlivFK4L055s15/ix0DFCC6xngBOC5YXgO8NsQS7XlZb44qyY/Q/o6w3ATcG9Id871MsVrqUKxlyDGa8lxjRnX/g7reQ/wX4TzRTXlo7vXXw0WcBnwtTD8NeBVxS4YStEvA26ezvIzNGnc7v5bd98chncB+4AZv6Byis4Ftrj7VncfBW4gij1bdlpuBn4v5O1lwA3uPuLu24AtYX1VEbe73+Xug2H0HqL3qsWpmLzO5wLgdnc/5O7dwO3AhWWKM9tUY34D8I0KxDUpd/85cKjALJcB13nkHmCemZ1AfHk9aczufneICarjmC4mn/OZye+hXlVrnuQ7n1XsN5TnOCtJXOG7ue5+j0dXatcxjWuFKf4W8p0/cx4Dpbqecffd7n5/GO4DHgOWUX15mS/OfCqenyFP+sNoU/h4gfVO9Vpqxv8PCsSYTyz728yWAxcD14TxQvun4vkIddhEEDje3XeH4T3A8XnmazWzDWZ2j5mN7YSFwGF3T4XxHRT+gZZSsXEDYGbnEpWsn8ya/IlQRfsZM2spU5zLgGeyxnPl0fg8IS97iPK2mGXLZarbfhvRnZUxuY6Xcis25teE/X6zmY29eDuuvC56uxY1wVwD3Jk1OY58Lla+tMV5XE/FxGPagR+b2UYzuyKmmPJ5YWii8kMzOy1Mq5V8rqRqyJNcx1G+81ncv6FSxbUsDJcr3neG/+lfsSOPC0w1xpJfz4SmVWcT1WpUbV5OiBOqKD9Ds7YHiW6Q3050DZdvvVO9lirJ72hijO4+lo+5rjHj2t+fBd4HZMJ4of0TSz7WZAHLzO4ws0dyfCbWSDj5S96r3H098Ebgs2Z2Uo3ETbgDcD3wVncfO7g+CJwKPI+oSvb95UlF/TOzPwbWA/+YNbnix0uRvg+sdvcziP5Zf22S+avJ5cDN7p7Omlat+VzTzOylRAWs7P8LL3b35xI1HXyHmf1OLMEd636i4+BM4PPAd2OORworeBxNdj6LS7XGBfw7cBJwFtFjAP8cbzgRM+sEvgW82917s7+rprzMEWdV5ae7p939LKLWBOcSXbdVlYkxmtnpVNE1ppldAuxz941xxVCMmixgufvL3f30HJ/vAXtDAWSsILIvzzp2hr9bgZ8S3e04SFS9OfYC5uXAzmqK28zmAv8NfDhUx46te3eooh0Bvkr5mt7tBFZkjefKo/F5Ql52EeVtMcuWS1HbNrOXAx8GXhnyEsh7vJTbpDG7+8GsOK8Bzil22TKZynYvZ0LzwJjyuVj50hbncT0pMzuD6Ni4zN0Pjk3Pyut9wHeoXHPdgty9d6yJirv/AGgys0VUeT7HJPY8yXMc5Tufxf0bKlVcOzm6uW3J4nX3veECNwN8iSO/y6nGWLLrGTNrIiq0fN3dvx0mV11e5oqzGvMzxHUYuAt4YYH1TvVaqqS/o6wYLyxwjRnH/n4R8Eoz207UfO9lwL9Sbfno03iwrJo/RLUO2Q9efjrHPPOBljC8CNjMkYcYv8nRD8m9vYribgZ+QnRnZuJ3J4S/RlR1+qkyxdlI9LDiGo48/HfahHnewdEPFN4Uhk/j6AcKt1K5Ti6Kiftsour6tcUeL1UQ8wlZw38A3BOGFwDbQuzzw/CCaog5zHcq0cPxFnc+T4hrNfk7X7iYox/mvS/OvC4y5pVE7crPmzC9A5iTNXw30Um0GmJeMnZcEJ3Enw55XtSxNZs+cedJvuOIPOezSv+GJh5npYyLYx/Uf0WJYsz+n/7XRM+IQJ7zZ6FjgBJcz4T0XQd8dsL0qsrLAnFWTX4SPTM/1jFZG/AL4JJ862WK11KFYi9BjDmvMePa31nxns+RTi6qJh/dvS4LWAuJCiGbgTuyduh64JowfB7wcMi0h4G3ZS1/Ytj5W8LOaqmiuP8YSAIPZn3OCt/dGdLyCPCfhF5gyhTrK4h66HmSqCYN4GNEtT4ArSHvtoS8PDFr2Q+H5Z5gmj3HlDHuO4C9WXl7y2THSxXE/ElgU4jtLuDUrGX/NOyDLUTNSasi5jB+JRNuAsSZz2H73yBqQpIkanP9NuAvgL8I3xvwxZCuh4H1VZDXk8V8DdCddUxvCNNPDPn8m3D8fLiKYn5n1jF9D1mFw1zH1mz/xJkn+Y4j8p/PKvYbynOclSwuonPzI2GZL5B1s2iGMV4fYngIuIWjCwg5z5/5jgFKcD0DvJio+d9DWf9HXlGFeZkvzqrJT+AM4IEQyyPARwqtl2lcS+WLvQQx5rzGjGt/Z63rfI4UsKomH919/C6hiIiIiIiIzFBNPoMlIiIiIiJSjVTAEhERERERKREVsEREREREREpEBSwREREREZESUQFLRERERESkRFTAEhERERERKREVsEREREREREpEBSwREREREZESUQFLRERERESkRFTAEhERERERKREVsEREREREREpEBSwREREREZESUQFLRMrOzBJm1m9mK+OORUREROclKScVsEQCM9tuZkPhH+7YZ+kM1vdyM9s+yTz/aWajYVuHzOzHZvasItffaGZuZqunG2OluHva3Tvd/WkYT/eVMYclIlLVdF4qH52XpJxUwBI52qXhH+7YZ1cFtvn37t4JLAP2AV8q1YrNrLFU6xIRkVjovCRSY1TAEpmEmTWY2c1mtsfMDpvZT83s2VnfX2Jmj5lZn5ntMLO/NrMu4PvAyqy7jscV2o67DwE3AWdN2P6fmdnjZtZtZj80sxXhq5+Hv5vC+l8zdnfSzD5kZnuAL5nZQjP7gZntD+v4vpkty1r/L83sb83s7pCGH5nZgvBdu5n9l5kdDGm/z8wWZS33MTO7x8wGzOy7YVvfMLNeM7t3rOlF9l1NM3s78HrgQyHu78xg94iIzDo6L+m8JNVNBSyR4twKrAWWAI8A12d991Xgbe4+BzgD+Jm79wCXAk9n3XXcV2gDZtYJvAHYkjXtNcD/BS4DFgP3Av8Vvv6d8Pe0sP5vhfHlQCewEng70e/8S2F8FZAE/nXC5t8IvBk4HugA3hOmvxVoD+tcGNY3nLXc68Oyy4FTgbuBq4EFwJPA/5uYTnf/N+BGwh1Sd/+DQvkiIiI56byk85JUKRWwRI723XBH7LCZfRfA3TPufq2797n7MHAlcI6ZdYRlksA6M5vj7ofc/f4pbvMDZnYY6APOJTqhjPkLon/4T7h7Cvg74NzsO305pIAr3X3U3Yfcfb+7fycM9wJ/D/zuhGW+7O6b3X0Q+CZH7lYmgUXAyaG9+gZ3789a7ivuvtXdu4HbgN+6+10h1m8CZ08xL0RE5Gg6L+m8JDVGBSyRo73K3eeFz6tgvKehT5vZVjPr5cidvEXh7x8ArwSeDs00nj/FbX7K3ecBa4hOHGuzvlsFfHHs5AocADJEd+by2evuo2MjZtZpZteY2dMh/juzYh+zJ2t4kOhOI8C1wB3ATWa208w+ZUe3n9+bNTyUY7wTERGZCZ2XdF6SGqMClsjk3gS8AngZ0AWcHKYbgLvf6+6vBI4jarJxQ/jep7IRd98O/DXweTNrCZOfIWrmMS/r0+bu9xZY/8Tp/5foJHmuu88N6Sg2plF3v9Ldnw28mOik/UfFpyr/qkuwDhGR2UrnJZ2XpIqpgCUyuTnACHCQqN33J8a+MLM2M3ujmc119yRRc4pM+HovsMjM5hS7IXf/IdHdwD8Lk64CPjz28LKZzTOz14Z50yGmE4uIfxDoNrOFwEeKjcfMXmZmp5tZA9BLdCczM8lixdjL5HGLiEhuOi/pvCRVTAUskcl9FdgVPpuIHpjN9mbgqdDM4W3AHwO4+yPAt4DtoSlFwd6asvwj8H4za3b3bwL/AnwzrP8h4IKseT8K/FdY/6vzrO9fiO5wHgyx/7DIOACWAt8mOoltImqW8V8FlyjONcCZofeom0uwPhGR2UTnJZ2XpIqZu2pERURERERESkE1WCIiIiIiIiWiApaIiIiIiEiJqIAlIiIiIiJSImUvYIV3NTxgZreG8TVmdq+ZbTGzG82sOUxvCeNbwvers9bxwTD9CTO7IGv6hWHaFjP7QLnTIiIiIiIiUkjj5LPM2F8BjwFzw/g/AJ9x9xvM7Cqi3m3+PfztdveTzezyMN/rzWwdcDlwGlHPMXeY2bPCur4I/D6wA/i1md3i7o8WCmbRokW+evXqkiZQRERKa+PGjQfcfXHccVSCzksiIrWh2HNTWQtYZrYcuJjo/QzvMTMjepncG8MsXwOuJCpgXRaGAW4GvhDmvwy4wd1HgG1mtgU4N8y3xd23hm3dEOYtWMBavXo1GzZsKEn6RESkPMzsqbhjmMjMvgJcAuxz99NzfG/AvxK9AHYQeIu73z/ZenVeEhGpDcWem8rdRPCzwPs48gK4hcBhd0+F8R3AsjC8jOjt4ITve8L849MnLJNv+jHM7Aoz22BmG/bv3z/TNImIyOx0LXBhge8vAtaGzxVENw9FRGSWKVsNlpmN3eXbaGbnl2s7xXD3q4GrAdavX68Xf4mIyJS5+8+znw/O4TLgOo9eMHmPmc0zsxPcfXdZA+vbw28e/g27G5awrPcBFnc04u480wfLF3SwtddINs9lPn2k+vaxZmEHh4eSHBhJkEgPMb+ri+HBfnppp6fr2Sw6cB8AhpHsWELLyAHSqSTDLYtpHYluUh6edxrzDm8CIJVoZ0lXOwcOHQBg/rJTSBzexoGBUdwS9HSdSmNqAHNn1Vxj/+6naT9uDanefSSH+0g3tJBubKd5tJtUYyctmQG62ppoamhgC8uY0/fkeFJHWhZgmTRYA8NzVnJO4zb6hkbpaG7EWzrZunPP+LwdS9aytz/J2e0H2N0zTEPnInoP7uH4Oa0cGhxl2JsYXLCOOQceoL25kYY5S+ju7aEp2QdAZ0sjT3WeScfADlqH95JOtJFID2XFspDWkUO0NjWwtKuNvX3DzG9vpnc4Sdep57N9xw44sIWR1oW0De1lpGUBxzf0cXhoFIDB9uU0ZJI0pvppTA2wemEHj7Oa9HAfyePOZP62W3GiS5YFHc3sSXaw7LiFcPgZ9veP0GBGY4PR/qzfpXnH3bQ0JugdTtIzlARgzaIOhkbT9I+kONSynNbmJjpshL79TwPQmGgglc5w2opFbHrmABlrom3lWSxodloOPUZjcxtbdh/g+Dmt7O0bZs78xQwePkDanVRjJ5mGRtwSzM8cZtm5l5H+7e3sHWogkRqgdziFNXfQ37SQZb6XRnP2DRmJrhPw7u3HHMIjLYtoyIzSlOzl+DmtjKQyDCfTzO9oYnN6CZ390TKLFp/AbltM076HOO2k1Ty0fS8DHStpG9o1vt9aW1tpWHwKg9YBT9/N8XNa2ZNqZ7BjBankKAsGttLcMY91607nVxvvhwWrOa/lKbbs6+fw0Ch9c04m2dSJeYa2oT20Du8bj3O49XjWtvVxsH0Ngzuj439hRzMHB0bH52lvbqSrrZH9tpjjOcDOw0McnncazSPdJJvn0dXzKIbR2ABL57expaeB0ZYFzBvaQTKTGd/fc7sW8NSuvePHQCrRTmN6kME1F9C+7bboGJ1/PJ02zO7ufpKNnTQne6I8aEzQnDAGk2lSmSOXvSct7mRgJIU7PJXsoiU9yOicZbQcehyAZx03h9/u6+OErtboN2OGOzz7hDk8smeQTKKFNkZI4LTOP54TEv08tqeX1gXL2JlYQevgTjoGnsGtgVSinTVz0uzuGaavdRkLRncznErTO/cUhtqOZ373wzSPdtN0/CksGn6K3T3DIUoDHFvxfIb2baV1ZD9DbUtoG9rDQMdKOgaeZlFnC/3DKWzF89jRM8qSgccZbJzHUGIObUO7aUr2smJ+O890D0b/C9o72N7xHOYsWMqLTl5I1NigvMr2omEz+yTwJ0AKaCV6Bus7RG/7XuLuKTN7IXClu19gZreF4V+ZWSOwB1gMfADA3T8Z1nsbR5oSXunuF4TpH8yeL591a1f5dZ/5cEnTKhK39ZdcEXcIIiVlZhvdfX3ccUwUCli35mkieCvwKXf/ZRj/CfB+dz+m/Z+ZXUFUy8XKlSvPeeqp6beI7N9wA4/v6Z328lJ9BjpW0DHwzOQzyrQt6mzhQP8IAMvnt7MjXIxL/dq5/GLOXD6P1Ys6pr2OYs9NZWsi6O4fdPfl7r6aqJOKO939j4C7gNeG2d4MfC8M3xLGCd/fGe4C3gJcHnoZXEPU9OI+4NfA2tArYXPYxi3lSo+IiEipuPvV7r7e3dcvXjyzvjwyGTXMqDeJ9OjkM8mMpNKeNZwpMKfUk1SF/l9WohfBid4P3GBmfwc8AHw5TP8ycH3oxOIQUYEJd99kZjcRdV6RAt7h7mkAM3sncBuQAL7i7psqmhIREZEjdgIrssaXh2kiIjKLVKSA5e4/BX4ahrdypBfA7HmGgdflWf4TRD0RTpz+A+AHJQxVRERkum4B3hl6tX0+0FP2569ERKTqxFGDJSIiUnPM7BvA+cAiM9sBfBRoAnD3q4hu+L0C2ELUTftb44lURETipAKWiIhIEdz9DZN878A7KhSOiIhUqXK/B0tERERERGTWUAFLRERERESkRFTAEhERERGRWaAy3bSrgCUiIiIiIlIiKmCJiIiIiIiUiApYIiIitcziDkBERLKpgCUiIrOOmbWZ2SlxxyEiIvVHBSwREZlVzOxS4EHgR2H8LDO7Jd6oRESkXqiAJSIis82VwLnAYQB3fxBYE2dAIiJSP1TAEhGR2Sbp7j0TplWm714REal7jXEHICIiUmGbzOyNQMLM1gJ/Cdwdc0wiIlJmXqFbaarBEhGR2eZdwGnACPANoBd4d6wRiYhI3VANloiIzCruPgh8OHxERERKSgUsERGZVczsLnI8c+XuL4shHBERqTNFFbDM7Dnu/nC5gxEREamA92YNtwKvAVIxxSIiInWm2BqsfzOzFuBa4Os5el8SERGpCe6+ccKk/zGz+2IJRkRE6k5RBSx3f0noaelPgY3hRPRVd7+9rNGJiIiUmJktyBptAM4BumIKZ8Ys7gBEROQoRT+D5e6bzexvgA3A54CzzcyAD7n7t8sVoIiISIltJHoGy4iaBm4D3hZrRCIiUjeKfQbrDOCtwMXA7cCl7n6/mS0FfgWogCUiIjXB3dfEHYOIiNSvYmuwPg9cQ1RbNTQ20d13hVotERGRqmZmry70vVpjiIhIKRRbwLoYGHL3NICZNQCt7j7o7teXLToREZHSubTAd45aY4iISAkUW8C6A3g50B/G24EfA+eVIygREZFSc/e3xh2DiIjUv2ILWK3uPla4wt37zay9TDGJiIiUlZldDJxG9B4sANz9Y/FFJCIi5XbMG+bLpKHI+QbM7LljI2Z2DjBUYH4REZGqZGZXAa8H3kXUk+DrgFWxBiUiInWj2BqsdwPfNLNdRCejJUQnJxERkVpznrufYWYPufvfmtk/Az+MOygREakPxb5o+NdmdipwSpj0hLsnyxeWiIhI2Yy1wBgMrxs5CJwQYzwiIlJHin7RMPA8YHVY5rlmhrtfV5aoREREyudWM5sH/CNwP1Gz/C/FG9L0mVgfA3QAACAASURBVFncIYiISJZiXzR8PXAS8CCQDpMdUAFLRERqirt/PAx+y8xuJerIqSfOmEREpH4UW4O1Hljn7pXqfENERKQszOwh4AbgRnd/EhiJOSQRiZEqgaXUiu1F8BGiji1ERERq3aVACrjJzH5tZu81s5VxByUi8VD1gZRasQWsRcCjZnabmd0y9ilnYCIiIuXg7k+5+6fd/RzgjcAZwLaYwxIRkTpRbBPBK8sZhIiISCWZ2Sqi1428nujZ4vfFG5GIiNSLYrtp/1k4Ga119zvMrB1IlDc0ERGR0jOze4Em4Cbgde6+NeaQRESkjhTbi+CfA1cAC4h6E1wGXPX/s3fn8ZLV9Z3/X+9a7r3dsjRbkL2BoD6AIQgtruPPdYK4kIwa0UlEgyGLGDPOPDKgM2o0jksmozFxYpCgYByRoEk6RGNASeKMQWyQsOggLYvSIFs33fRyb92q8/n9cb63u/p21b11u6vqnKr7fj4e9bhnrXrXt751z/nW+Z5zgJcOLpqZmdlAvDki7io6hJmZjadez8F6O/B8YAtARNwN/MygQpmZmQ2KG1dmZsvTsC5o0msDayYiGnMjkmrk98EyMzNbNiSdLekuSeslXdxh/lskPSrp1vR4WxE5zcysOL02sP5J0ruBFZJeDvwl8LcLrSDpGEk3SPq+pDslvTNNP1jSdZLuTn8PStMl6ZNpo3WbpDPanuv8tPzdks5vm36mpNvTOp+Ub2dvZmYDIqkKfAp4BXAy8EZJJ3dY9EsRcXp6XDbUkGZmVrheG1gXA48CtwO/DnwV+K+LrNME/lNEnAw8B3h72hBdDHwjIk4CvpHGId9gnZQeFwJ/CnmDDHgf8GzgLOB9c42ytMyvta13do/vx8zMlilJKyX9N0mfSeMnSXpVD6ueBayPiHtSr46rgHMHmdXMBsUdsWxwempgRUQWEZ+JiNdHxOvS8II1MyIeiohb0vCTwA/IL45xLnBFWuwK4BfS8LnAlZG7EVgl6Qjg54HrImJjRGwCrgPOTvMOiIgbU5Yr257LzMysm88CM8Bz0/gG4Pd7WO8o4Cdt4w+kafO9NvXEuEbSMZ2eSNKFktZJWvfoo48uIbqZ9Yc7Pdng9NTAknSvpHvmP3p9EUmrgWcC3wEOj4iH0qyfAoen4W4broWmP9BheqfX37kh27R5a6+xzcxsPJ0YER8DZgEiYjv929v6W2B1RJxG/oPgFZ0WiohLI2JNRKw57LDD+vTSZmZWBr3eaHhN2/AU8HryS7YvStJ+wJeB34mILe2nSUVESBr4MdqIuBS4FODkk47zMWEzs+WtIWkFqY+QpBPJj2gtZgPQfkTq6DRtp4h4vG30MuBj+xbVzMxGTa9dBB9ve2yIiE8Ar1xsPUl18sbVFyLiK2nyw6l7H+nvI2l6tw3XQtOP7jDdzMxsIe8D/h44RtIXyM8H/t0e1vsucJKk4yVNAOcBa9sXmNu+Ja8h7x5vZmbLSK83Gj6jbbRCfkRrwXXTFf3+HPhBRPzPtllrgfOBj6S/f9M2/SJJV5Ff0GJzRDwk6evAf2+7sMW/Ay6JiI2Stkh6DnnXwzcDf9zL+zEzs+UrIq6TdAv5BZgEvDMiHuthvaaki4CvA1Xg8oi4U9IHgHURsRb4bUmvIb/Q00bgLYN6H2ZmVk69dhH8w7bhJnAf8EuLrPN84FeA2yXdmqa9m7xhdbWkC4D7257nq8A5wHpgO/BWgNSQ+iD5L4cAH4iIjWn4t4DPASuAr6WHmZnZHub9WAgwdz7wsZKOnbsw00Ii4qvk26v2ae9tG74EuGRfs5qZ2ejqqYEVES9e6hNHxP+h+0nDL+2wfABv7/JclwOXd5i+Djh1qdnMzGxZ+sMF5gXwkmEFMTOz8dVrF8F3LTR/XhdAMzOz0tmbHwvNbPzJV2y3PlvKVQSfxa6TeV8N3ATcPYhQZmZmgyJpiryL+QvIj1x9C/h0REwXGszMCrHwnV3Nlq7XBtbRwBnphsFIej/wdxHxy4MKZmZmNiBXAk+y68JIbwI+T34LkpHjH9/NzMql1wbW4UCjbbzBrhsEm5mZjZJTI+LktvEbJH2/sDRmZjZWem1gXQncJOmv0vgv0OXu9GZmZiV3i6TnRMSNAJKeDawrOJOZmY2JXq8i+CFJXwP+bZr01oj43uBimZmZDcyZwLcl/TiNHwvcJel28ovanlZcNDMzG5RgOCfc9XoEC2AlsCUiPivpMEnHR8S9gwpmZmY2IGcXHcDMzMZXr5dpfx/5lQSfDnwWqAN/QX4zYTMzs5EREfdLOgg4hrbtYC83Gjaz8ePLtFu/9XoE6xeBZwK3AETEg5L2H1gqMzOzAZH0QeAtwI9gZ38R32jYbJnyZdqt33ptYDUiIiQFgKSnDDCTmZnZIP0ScGJENBZd0szMbIkqPS53taQ/A1ZJ+jXgeuAzg4tlZmY2MHcAq4oOYWbl4C6C1m+9XkXwf0h6ObCF/Dys90bEdQNNZmZmNhgfBr4n6Q5gZm5iRLymuEhmVhR3EbR+W7SBJakKXB8RLwbcqDIzs1F3BfBR4HYgKziLmZmNmUUbWBHRkpRJOjAiNg8jlJmZ2QBtj4hPFh2iX9y7ycysXHq9yMVW4HZJ1wHb5iZGxG8PJJWZmdngfEvSh4G17N5F0JdpN1uGfA6W9VuvDayvpIeZmdmoe2b6+5y2ab5Mu5WITwoaJp+DZf22YANL0rER8eOIuGJYgczMzAYpnVNsZmY2EIsdwfpr4AwASV+OiNcOPpKZmdlgSXolcAowNTctIj5QXCIzK4q7CC4fwzpaudh9sNqr3AmDDGJmZjYMkj4NvAF4B/l27vXAcYWGMrPCuIug9dtiDazoMmxmZjaqnhcRbwY2RcTvAc8FnlZwJjMzGxOLdRH8OUlbyH/hW5GGSeMREQcMNJ2ZmVn/7Uh/t0s6EngcOKLAPGZmNkYWbGBFRHVYQczMzIbkWkmrgD8AbiHvofGZYiOZWVF8Dpb1W6+XaTczMxsLEfHBNPhlSdcCUxGxuchMZlYcn4Nl/bbYOVhmZmZjQdKzJD21bfzNwNXAByUdXFyyfeRf38eP9/jNRpobWGZmtlz8GdAAkPRC4CPAlcBm4NICc5lZgdxF0PrNXQTNzGy5qEbExjT8BuDSiPgyeVfBWwvMZWYF8gFD6zcfwTIzs+WiKmnuh8WXAt9sm+cfHM3MrC+8QTEzs+Xii8A/SXqM/FLt3wKQ9LPk3QTNbBlyF0HrNzewzMxsWYiID0n6Bvk9r/4hYmfHoArwjuKSmZnZOHEDy8zMlo2IuLHDtB8WkcXMysHnYC0fw/qsfQ6WmZmZmZlZn7iBZWZmZmbLls/Bsn5zA8vMzMzMli13EbR+cwPLzMzMrESE9/jNRpkbWGZmZj2SdLakuyStl3Rxh/mTkr6U5n9H0urBhxr4K5iNNXcRtH5zA8vMzKwHkqrAp4BXACcDb5R08rzFLgA2RcTPAh8HPjrclDYOwq1ms5HmBpaZmVlvzgLWR8Q9EdEArgLOnbfMucAVafga4KXSYH8fX//w1kE+vRVgaubRoiOMvSd2NHYOP7R5usAkNiyHPnoj1R2PD+W13MAyMzPrzVHAT9rGH0jTOi4TEU1gM3DI/CeSdKGkdZLWPfrovu1MT9W8KR83jfqBRUcYeysndt0Kdr9J3xZ2OVBkrFpZH8prjXyNknQ28EdAFbgsIj5ScCQzM7MFRcSlwKUAa9as2acrGjzt5b/al0xmZtYfI93AausP/3LyXxK/K2ltRHy/2GRmw7Xu2kuLjmC2HGwAjmkbPzpN67TMA5JqwIHAcPqkmJlZKYx6v4Je+sObmZn1w3eBkyQdL2kCOA9YO2+ZtcD5afh1wDcjfJcdM7PlZKSPYNG5P/yz5y8k6ULgwjQ686xX//odQ8jWT4cCjxUdYi+MYu5RzAyjmduZh2cUcz+96ADzRURT0kXA18m7pV8eEXdK+gCwLiLWAn8OfF7SemAjeSNsQTfffPNjku7fx3ij8Bk7Y/+MQs5RyAijkdMZ+2dfcx7Xy0Ia5R/WJL0OODsi3pbGfwV4dkRctMA66yJizbAy9sMoZobRzD2KmWE0czvz8Ixi7lHMXKRRKC9n7J9RyDkKGWE0cjpj/wwr56h3EeylP7yZmZmZmdlQjHoDq5f+8GZmZmZmZkMx0udgdesPv8hqo3i5tVHMDKOZexQzw2jmdubhGcXco5i5SKNQXs7YP6OQcxQywmjkdMb+GUrOkT4Hy8zMzMzMrExGvYugmZmZmZlZabiBZWZmZmZm1idj18CSdLCk6yTdnf4e1GW5lqRb02Nt2/TjJX1H0npJX0oXzyhFbkmnS/oXSXdKuk3SG9rmfU7SvW3v6fQBZj1b0l2pjC7uMH8yld36VJar2+ZdkqbfJennB5VxL3O/S9L3U9l+Q9JxbfM61pcSZH6LpEfbsr2tbd75qT7dLen8+esWmPnjbXl/KOmJtnmFlHN67cslPSKp433ylPtkel+3STqjbV5RZb1Y5v+Qst4u6duSfq5t3n1p+q2S1pUo84skbW6rB+9tm7dg3VqOii6TTvWo2/ZsmN+hTvWsn7kknZne9/q0rvqU8f2SNrTV/3Pa5nXcfnarA+rD/oykYyTdoHzbeKekd6bpZSvLbjlLU56SpiTdJOlfU8bfW+h5tRf7Ut2y9yFjx33Moj7v9DxVSd+TdG0aL005AhARY/UAPgZcnIYvBj7aZbmtXaZfDZyXhj8N/GZZcgNPA05Kw0cCDwGr0vjngNcNIWcV+BFwAjAB/Ctw8rxlfgv4dBo+D/hSGj45LT8JHJ+epzqk8u0l94uBlWn4N+dyL1RfSpD5LcCfdFj3YOCe9PegNHxQGTLPW/4d5BenKayc2177hcAZwB1d5p8DfA0Q8BzgO0WWdY+ZnzeXBXjFXOY0fh9waAnL+UXAtftat5bDowxl0qke0WV7NszvUKd61s9cwE1pWaV1X9GnjO8H/nOHZTtuPxeqA/RhfwY4AjgjDe8P/DBlKVtZdstZmvJM72+/NFwHvpPed8fnZYn7Ugtl70PGz9FhH7Oozzs9z7uA/03aXpSpHCNi/I5gAecCV6ThK4Bf6HXF1Ip+CXDN3qy/jxbNHRE/jIi70/CDwCPAYUPKN+csYH1E3BMRDeAq8uzt2t/LNcBLU9meC1wVETMRcS+wPj1fKXJHxA0RsT2N3kh+X7Ui9VLW3fw8cF1EbIyITcB1wNkDytluqZnfCHxxCLkWFRH/DGxcYJFzgSsjdyOwStIRFFfWi2aOiG+nTFCOOt1LOXezL9+HcVXWMum2PRvad6hLPetLrjTvgIi4MfI9tSvZi32FJX4Xum0/O9aBfu3PRMRDEXFLGn4S+AFwFOUry245uxl6eaYy2ZpG6+kRCzzvUvel9vn/wQIZuynk85Z0NPBK4LI0vtDnM/RyhDHsIggcHhEPpeGfAod3WW5K0jpJN0qa+xAOAZ6IiGYaf4CFv6D91GtuACSdRd6y/lHb5A+lQ7QflzQ5oJxHAT9pG+9URjuXSWW5mbxse1l3UJb62heQ/7Iyp1N9GbReM782fe7XSJq78XZRZd3z6yrvgnk88M22yUWUc6+6vbci6/VSzK/TAfyDpJslXVhQpm6em7qofE3SKWnaqJTzMJWhTDrVo27bs6K/Q/3KdVQaHlTei9L/9Mu163SBpWbs+/5M6lr1TPKjGqUty3k5oUTlmbq13Ur+A/l15Ptw3Z53qftSffkezc8YEXPl2Gkfs6jP+xPA7wJZGl/o8ymkHEeygSXpekl3dHjMPyIRdG95HxcRa4A3AZ+QdOKI5Cb9AvB54K0RMVe5LgGeATyL/JDsfxnMuxh/kn4ZWAP8QdvkodeXHv0tsDoiTiP/Z33FIsuXyXnANRHRaptW1nIeaZJeTN7Aav+/8IKIOIO86+DbJb2wkHB7uoW8Hvwc8MfAXxecxxa2YD1abHtWlLLmAv4UOBE4nfw0gD8sNk5O0n7Al4HfiYgt7fPKVJYdcpaqPCOiFRGnk/cmOIt8v61U5meUdCol2seU9CrgkYi4uagMvRjJBlZEvCwiTu3w+Bvg4dQAmWuIPNLlOTakv/cA/0j+a8fj5Ic3527AfDSwoUy5JR0A/B3wnnQ4du65H0qHaGeAzzK4rncbgGPaxjuV0c5lUlkeSF62vaw7KD29tqSXAe8BXpPKEuhaXwZt0cwR8XhbzsuAM3tdd0CW8rrnMa97YEHl3Ktu763Ier0oSaeR141zI+LxueltZf0I8FcMr7vugiJiy1wXlYj4KlCXdCglL+eCFF4mXepRt+1Z0d+hfuXawO7dbfuWNyIeTju4GfAZdn0vl5qxb/szkurkjZYvRMRX0uTSlWWnnGUsz5TrCeAG4LkLPO9S96X6+j1qy3j2AvuYRXzezwdeI+k+8u57LwH+iLKVY+zFiWVlfpAfdWg/8fJjHZY5CJhMw4cCd7PrJMa/ZPeT5H6rRLkngG+Q/zIzf94R6a/ID51+ZEA5a+QnKx7PrpP/Tpm3zNvZ/YTCq9PwKex+QuE9DO8iF73kfib54fqTeq0vJch8RNvwLwI3puGDgXtT9oPS8MFlyJyWewb5yfEqupzn5VpN94svvJLdT+a9qciy7jHzseT9yp83b/pTgP3bhr9NvhEtQ+anztUL8o34j1OZ91S3ltOj6DLpVo/osj0b9ndofj3rZy72PFH/nD5lbP+f/h/JzxGBLtvPheoAfdifSe/vSuAT86aXqiwXyFma8iQ/Z37uwmQrgG8Br+r2vCxxX2qh7H3I2HEfs6jPuy3vi9h1kYvSlGNEjGUD6xDyRsjdwPVtH+ga4LI0/Dzg9lRotwMXtK1/Qvrw16cPa7JEuX8ZmAVubXucnuZ9M72XO4C/IF0FZkBZzyG/Qs+PyI+kAXyA/KgPwFQqu/WpLE9oW/c9ab272Msrxwww9/XAw21lu3ax+lKCzB8G7kzZbgCe0bbur6bPYD15d9JSZE7j72fejwBFlnN6/S+SdyGZJe9zfQHwG8BvpPkCPpXe1+3AmhKU9WKZLwM2tdXpdWn6Camc/zXVn/eUKPNFbXX6Rtoah53q1nJ/FFkm3eoR3bdnQ/sOdalnfctFvm2+I63zJ7T9WLSPGT+fMtwGrGX3BkLH7We3OkAf9meAF5B3/7ut7f/IOSUsy245S1OewGnA91KWO4D3LvS87MW+VLfsfcjYcR+zqM+77blexK4GVmnKMSJ2/kpoZmZmZmZm+2gkz8EyMzMzMzMrIzewzMzMzMzM+sQNLDMzMzMzsz5xA8vMzMzMzKxP3MAyMzMzMzPrEzewzMzMzMzM+sQNLDMzMzMzsz5xA8vMzMzMzKxP3MAyMzMzMzPrEzewzMzMzMzM+sQNLDMzMzMzsz5xA8vMzMzMzKxP3MAyK5ikd0u6rOgcZmZm4O2S2b5SRBSdwWysSdraNroSmAFaafzXI+ILw0+1MEmrgXuBekQ0i01jZmb95O2S2WDVig5gNu4iYr+5YUn3AW+LiOuLS2RmZsuZt0tmg+UugmYFk/R+SX+RhldLCklvlfQTSZsk/YakZ0m6TdITkv5k3vq/KukHadmvSzqubd7zJH1X0ub093lt8+6T9LJOOYB/Tn+fkLRV0nMHVwJmZlYm3i6Z7Rs3sMzK6dnAScAbgE8A7wFeBpwC/JKk/w9A0rnAu4F/DxwGfAv4Ypp3MPB3wCeBQ4D/CfydpEN6eP0Xpr+rImK/iPiXPr0vMzMbTd4umfXIDSyzcvpgRExHxD8A24AvRsQjEbGBfGP1zLTcbwAfjogfpD7p/x04Pf1a+Erg7oj4fEQ0I+KLwP8DXj38t2NmZiPO2yWzHrmBZVZOD7cN7+gwPtd//jjgj1IXjSeAjYCAo4AjgfvnPe/9aZ6ZmdlSeLtk1iM3sMxG20/Ir/i0qu2xIiK+DTxIvqFrdyywIQ1vI7961Jyntg378qJmZrY3vF2yZc8NLLPR9mngEkmnAEg6UNLr07yvAk+T9CZJNUlvAE4Grk3zbwXOk1SXtAZ4XdvzPgpkwAlDeRdmZjYuvF2yZc8NLLMRFhF/BXwUuErSFuAO4BVp3uPAq4D/BDwO/C7wqoh4LK3+34ATgU3A7wH/u+15twMfAv5v6ubxnOG8IzMzG2XeLpn5RsNmZmZmZmZ94yNYZmZmZmZmfeIGlpmZmZmZWZ+4gWVmZmZmZtYnbmCZmZmZmZn1Sa3oAMN26KGHxurVq4uOYWZmC7j55psfi4jDis4xDN4umZmNhl63TQNvYEmqAuuADRHxKknHA1cBhwA3A78SEQ1Jk8CVwJnkl+58Q0Tcl57jEuACoAX8dkR8PU0/G/gjoApcFhEfWSzP6tWrWbduXZ/fpZmZ9ZOk+4vOMJ+ky8kvMf1IRJzaYb7It0nnANuBt0TELYs9r7dLZmajoddt0zC6CL4T+EHb+EeBj0fEz5Lf5+CCNP0CYFOa/vG0HJJOBs4DTgHOBv6XpGpquH2K/N4KJwNvTMuamZkNwufIt0PdvAI4KT0uBP50CJnMzKxkBnoES9LRwCvJbwz3rvTr3kuAN6VFrgDeT74ROjcNA1wD/Ela/lzgqoiYAe6VtB44Ky23PiLuSa91VVr2+4N8T2ZmtjxFxD9LWr3AIucCV0Z+g8kbJa2SdEREPDTIXFmzyczsLBO1CjuawQo1mdUElaxBlYxGZSVqTTPThHq9Sk3QosJUTTSaTSq1SVpZRo0Ws1GlGk0mplZAq0kjKmQzW5mNClOTk1SrNRrTW2lpgomJOs3p7VTqE9SqNVrNBkRGpVKlItjRaDJVrxBUqNZq7JieZkW9SisLkGiRz6uQoWoNtRqoWqeVQYWMRqvFVLVC1pyhVp+k1WrSarWoVKtEpZa/TjRpZjBLlYlsmkp1ggCq1RotxGxjmhotMtXJshYTVZGpRnP6SahUmaxVyLKgVd+PLILJSsZ0U1RjlkrWoDa5kkYriMggWkRlkipN1GpQq08REjOzLeoVQPl7aapGRUJZk5lmENGiUq0zoYwsa0GlmpdXQIZoNlsEQYWMifokAZDNUq1NQNYCichaTLcqTNSr0GqAKrQaO6jUJmlk0MyCyZqYrNWZac4yUauRtZq0WhnUJqnRolKp0GxlNFvNvOQFldoUWatBRNCiQq1apdKcJqvUyFI+IqhU61QVzLZaqFKhFi0kEaoy22pRq9WpkpHNbIWJp6CsRVN1slYjX65SpyrIWg0m6hM0djxJpTZBpVIjWg1UX8lsc5ZmK2Oy0qJaqTGb5fVxRU1kWUZTEyia1Cv551qtT6DIqERGK2sR1UkmqqLZaqH6iryetWZpzGxjauUBzDYaTNRrtEJMqAm1KYiMRqNBhpCCWrVGIGqCyJrMRoUaGVmlRquxjdmoUa+KerVKFuR1tF5HlSpZ1iKyFlXBbFTJEFP1KmRNpmebiKBazXe7K5HRrExQV4tGln+Ha9UKlUqVubvTtlpNJup1IstoRIWYnaZSrVCv1vI6TjBLhVZzlkq0mGnByqkpGs2MFZN1ZhvTtLJA1RoVgpoyWiFq1SpZljETtfwzy5oARBZM1vNyDyooLxUy8nUUGRUFs5lotTJW1IJWlj+yCKYmJsgiaGakel6j0ppGBM2oUp+cIms2yFTN618qs0ZUWFHJoDpBM92bt5LldbRaq0PWJMsyGs0mrQym6lWaVJio1YBge6NFtVqlqkCIakU0MphtTDMxuZJazFJJZZZlQbU2nLOjBv0qnyC/S/f+afwQ4ImIaKbxB4Cj0vBRwE8AIqIpaXNa/ijgxrbnbF/nJ/OmP7tTCEkXkv+ayLHHHrsPb8fMzKyrnduxZG57NdAG1i1/f/k+rb995dGs3P7AbtNOffEbmLrnOm67f+Nu0x879CwOfewmALYc8HQO2HIXADOTBzM5s/uyvchUoxJNntz/RPZ/8kcpz1Gs3L5ht+VOO3oVtz3wxB7r7zdZY+vM7usDrFl9KP+46Wc4cPPCv7keMFVjy3STxsRBPHrYszlqw9+zY+pwVkw/DMCpRx7IHQ9ubltDkHaBjz/0KTz65AxbZ5p7PvEiTn7qAXz/p1uYnjyMqZlHd07fseKpBBVW7niQM489iPx3Zrjl/k1sn/oZnjjo33DEQ9d3fd7D95/i4Send5YLwNb9VrPf1vs48sAVPLh5x27Lz0wewuTM44vmbVWnqLamd5s2Vasy3WwB8OCR/44THvnGzvGVEzW2ZFPUmlsBePjwF7HqiTuYnHms4/OHqihaC2Z48MiXc+SD1y24zKoVdZ7YMcuOFU+lMXFw18//hEP34+BVB9LYsY3bNuxZr04+8gA2bWvw0ObpDmtDrSKaWewcnzz0BGYeu2e3ZXaseCpH/szhHLrjR9z9yNY9nmN+nZ2TqQ5AJWapV0QW8OSKo/b4jp521KqO2QEOPmA/Nm7Z/TWfMlFjW6PJEQdO8dDmaTYf+AwO3Pz/Oq7fzQGTNbakejX/O/nUA6b46ZbO5dWLHSuOYOMhZzA5/ejO/zFrjjsYgFvu7/y/5dCnTPDYtsZu0+rVCrOtDICt+x3PodM/5tSjDuSeLRU2bXqMNa+6cK8zLsXAGliS5vqp3yzpRYN6nV5ExKXApQAnn3RcrLv20iLjmPXdsP5hmFl/lOmHv/k7bgAzW59gqsOyU9O7GgMrduxqN+5N4wqgkn5vbd/RnN+4Amg0s47rzzUi9thRjWy3fN1smc7Xn2hsopLlO/hzjSuA6dn5O/27dqo375jdq8YVwPb0vO2NK4AVO37acfkgWDH9ME+2fnbB533kyRmA3XLtt/U+ADZtb+yxfC+NK2CPxhWwszEFUMkau41vbzSpsWsHv9bc1rVxBSzauMoz7Jl/vid2zAJ5OVZbM12X27KjwcFP2bFb5nbTjYzHt3Z/vfbGFbBH42ouw2MPbKa+MvaYB7Dfk/d2nF6J2Z3Ds+l1On1Hu2UH9mhcAWxrq9qozAAAIABJREFU5HVi07ZdZbRUW9rq1ey87+TDW7qXdy/mvq9L+V8yv3EF7GxcAey39V7mau6mTd3r3yAM8hys5wOvkXQf+UUtXkJ+8u8qSXMNu6OBuf+kG4BjANL8A8kvdrFz+rx1uk03MzMrQs/bpYi4NCLWRMSaww5bFhdLNDNbNgbWwIqISyLi6IhYTX6Rim9GxH8AbgBelxY7H/ibNLw2jZPmfzP1Y18LnCdpMl2B8CTgJuC7wEmSjpc0kV5j7aDej5mZ2SLWAm9W7jnA5kGff2VmVoTOx+VsThH3wfovwFWSfh/4HvDnafqfA59PF7HYSN5gIiLulHQ1+cUrmsDbI/JjyZIuAr5Ofpn2yyPizqG+EzMzWzYkfRF4EXCopAeA9wF1gIj4NPBV8ku0rye/TPtbi0k6OIGKjmAjqVy74+VKU5AoZymMy/+YoTSwIuIfgX9Mw/ew6yqA7ctMA6/vsv6HyK9EOH/6V8k3aGZmZgMVEW9cZH4Abx9SnMKNy46Q9cMw6sLo1bf8unblbMjYYA3jPlhmZmZmpRajt/9uAzR3BUezveEGlpmZmc3jncui+biH2ehyA8vMzMx6M8a/6sstmmVnfGvz8IV/EtiNG1hmZmZmQ+SdUbPx5gaWmZmZmZmNvRjS1RPdwDIzM7PdLMfjK6NwkQuV9NLao2axUhznI4zj+87KxQ0sMzNbdiStkPT0onPY8iSf/WM21tzAMjOzZUXSq4Fbgb9P46dLWltsKjMzGxduYJmZ2XLzfvIb3j8BEBG3AscXGWh0FHPkZVjnTQzLOHdB62y5vd/y8w2QB8sNLDMzW25mI2LzvGne29hNDw2pkl+yPUqez5ZmFM6R24Pr4NKNSZnVig5gZmY2ZHdKehNQlXQS8NvAtwvOZGZmY8JHsMzMbLl5B3AKMAN8EdgC/E6hiWz0+JinmXXhI1hmZrasRMR24D3pYR10azuEr35nHY1mvdCYdEez8nEDy8zMlhVJN9ChDRERLykgjpmNM7fhliU3sMzMbLn5z23DU8BrgWZBWczMRs6oXtkzYjjX0eipgSXp30TE7YMOY2ZmNmgRcfO8Sf9X0k2FhCkrd50aO/5El2A02w62AKGh3h6h1yNY/0vSJPA54AsdLm9rZmY2EiQd3DZaAc4EDiwozgjzLrvZYvwtWZ56amBFxL9Nl7L9VeDm9EvfZyPiuoGmMzMz67+byX+jFnnXwHuBCwpNZGZmY6Pnc7Ai4m5J/xVYB3wSeKbyy6+8OyK+MqiAZmZm/RQRxxedYWQV1HXQPbZG3Iier1Nq+1yk/kwGqddzsE4D3gq8ErgOeHVE3CLpSOBfADewzMys1CT9+4Xm+8fCcePOWTnhnWmz4er1CNYfA5eRH63aMTcxIh5MR7XMzMzK7tULzAv8Y2EbN06sbEaxTg4w8ygWxzLSawPrlcCOiGgBSKoAUxGxPSI+P7B0ZmZmfRIRby06g42PYo4J+UiUlYNr4sJ6bWBdD7wM2JrGVwL/ADxvEKHMzMwGSdIrgVPI74MFQER8oLhE5TKq97gZJzGsG/bYsuLv9nBUelxuKiLmGlek4ZWDiWRmZjY4kj4NvAF4B3lHm9cDxxUayszGktvIy1OvDaxtks6YG5F0JrBjgeXNzMzK6nkR8WZgU0T8HvBc4GkFZzIzswEb1vG7XrsI/g7wl5IeJP+176nkv/6ZmZmNmrkfCLenq+E+DhxRYJ7y8c/uYyMQwt0NdzMOveTG4T2MsV5vNPxdSc8Anp4m3RURs4OLZWZmNjDXSloF/AFwC/muymeKjWTLineObdS5Di+o5xsNA88CVqd1zpBERFw5kFRmZmYDEhEfTINflnQt+XnGm4vMZGaL8RG43bg4Sq3XGw1/HjgRuBVopckBuIFlZmYjRdJtwFXAlyLiR8BMwZFsufHOsdlY6/UI1hrg5PC1Hc3MbPS9mvw84qslZcCXgKsj4sfFxrKu9mrvY7RbMaO0wxWLnN+lkr6b0a4h+6i0u/Tj8an0ehXBO8gvbGFmZjbSIuL+iPhYRJwJvAk4Dbi34FglM2o7OaOW18zGWa9HsA4Fvi/pJtq6UkTEawaSyszMbIAkHUd+FOsN5F3ff7fYRLZvhv9r/D4dACjrwYO9pNIeDbFBiXGrxH3WawPr/YMMYWZmNiySvgPUgauB10fEPQVHKp3uu04+UmTjbySaDnsZciTe2xjo9TLt/5R+7TspIq6XtBKoDjaamZnZQLw5Iu4qOoTZUIiB71Uvdg6W7ams56VZf/R0DpakXwOuAf4sTToK+OtBhTIzMxsUN67MbBgG2oRymxaVuBB6vcjF24HnA1sAIuJu4GcGFcrMzKyMJJ0t6S5J6yVd3GH+WyQ9KunW9HhbETnNzGxPw7ogeq/nYM1EREPpELCkGu7GaWZmy4ikKvAp4OXAA8B3Ja2NiO/PW/RLEXHR0AP2k7t8ldTe73ott090oW6Lw9yBLfNRluVEGu6V6Xs9gvVPkt4NrJD0cuAvgb9daAVJx0i6QdL3Jd0p6Z1p+sGSrpN0d/p7UJouSZ9MvwreJumMtuc6Py1/t6Tz26afKen2tM4nJW8RzMxsYZJWSvpvkj6Txk+S9KoeVj0LWB8R90REg/xmxecOMqv1ypv/cefzvGyU9NrAuhh4FLgd+HXgq8B/XWSdJvCfIuJk4DnA2yWdnJ7rGxFxEvCNNA7wCuCk9LgQ+FPIG2TA+4Bnk2/c3jfXKEvL/Frbemf3+H7MzGz5+iz5LUeem8Y3AL/fw3pHAT9pG38gTZvvtemHwmskHbNPSa1H7lQz7pZ6Kfixv3T8mL+9UddTAysisoj4TES8PiJel4YX/Ggj4qGIuCUNPwn8gHxDdC5wRVrsCuAX0vC5wJWRuxFYJekI4OeB6yJiY0RsAq4Dzk7zDoiIG1OWK9uey8zMrJsTI+JjwCxARGynf4dA/hZYHRGnkW+vrui0kKQLJa2TtO7RRx/t00ubWdn4flHLU69XEbxX0j3zH72+iKTVwDOB7wCHR8RDadZPgcPTcLdfBhea/kCH6Z1ef+eGbNPmrb3GNjOz8dSQtIL0G7CkE8mPaC1mA9B+ROroNG2niHg8Iuae6zLgzE5PFBGXRsSaiFhz2GGHLTW/mY0In4O1PPV6kYs1bcNTwOuBg3tZUdJ+wJeB34mILe2nSUVESBp40z4iLgUuBTj5pOP8U4KZ2fL2PuDvgWMkfYH8Krlv6WG97wInSTqevGF1HvCm9gUkHdH2I+JryHtvmBXIO/h7w6U2nobVCOj1RsOPz5v0CUk3A+9daD1JdfLG1Rci4itp8sNzG6DUze+RNL3bL4MbgBfNm/6PafrRHZY3MzPrKiKuk3QL+fnBAt4ZEY/1sF5T0kXA14EqcHlE3CnpA8C6iFgL/Lak15Cfh7yR3hpuJTT6u5cxBu/B9l4vF8VYaAn/Gl+McbmYSU8NrPYr+pF3K1yz2Lrpin5/DvwgIv5n26y1wPnAR9Lfv2mbfpGkq8gvaLE5NcK+Dvz3tgtb/DvgkojYKGmLpOeQdz18M/DHvbwfMzNbfuZtywDmjjQdK+nYufOGFxIRXyW/0FP7tPe2DV8CXLKvWW13y/I8lrF6y2P1ZsbCeDRjyqvXLoJ/2DbcBO4DfmmRdZ4P/Apwu6Rb07R3kzesrpZ0AXB/2/N8FTgHWA9sB94KkBpSHyTvmgHwgYjYmIZ/C/gcsAL4WnqYmZl18ocLzAvgJcMKYuPAjQZbxIJHY1x/xlmvXQRfvNQnjoj/Q/cG8ks7LB/A27s81+XA5R2mrwNOXWo2MzNbfvZmW2Zmti/GpMebLVGvXQTftdD8eV0AzczMSkvSFHkPiBeQ/4z8LeDTETFdaLASKfYWQsK/7vdPpOvYeUe/dzHOZ/D5qzUUS7mK4LPIz5MCeDVwE3D3IEKZmZkN0JXAk+w6b/dNwOfJr5BrZmaLKMN9nKVy5Oik1wbW0cAZ6YbBSHo/8HcR8cuDCmZmZjYgp0bEyW3jN0j6fmFpysiHO8zM9lpPNxomvxlwo228wa4bBJuZmY2SW9IVaAGQ9GxgXYF5Rsa4XELZ+s31wvplMHVp2DW01yNYVwI3SfqrNP4LwBWDiWRmZjZQZwLflvTjNH4scJek28mvuXRacdHMbNDK2q3MxkevVxH8kKSvAf82TXprRHxvcLHMzMwG5uyiA5jtDbllYCVRhppY5uOmvR7BAlgJbImIz0o6TNLxEXHvoIKZmZkNQkTcn25efwxt28FebjRsZmaja1i/UfR6mfb3kV9J8OnAZ4E68BfkNxM2MzMbGenm9W8BfsSuH2J9o+HdlPm34R6N+Pli/dsPLEE5lOFwRyclKJrilPVD6V2Z30GvR7B+EXgmcAtARDwoaf+BpTIzMxucXwJOjIjGokuadbEvO3dl3jE0s33X61UEGxERpP8Jkp4yuEhmZmYDdQewqugQZma298p8ALLXI1hXS/ozYJWkXwN+FfjM4GKZmZkNzIeB70m6A5iZmxgRryku0qgY/C5NIORjPGbWT2Koh457vYrg/5D0cmAL+XlY742I6waazMzMbDCuAD4K3A5kBWcppXADx5azsa7+fXpzY11G+27RBpakKnB9RLwYcKPKzMxG3faI+GTRIUptxC8QYVYO/h4tV4s2sCKiJSmTdGBEbB5GKDMzswH6lqQPA2vZvYugL9M+VkZ157bPfZkkH22wETKq39vd9XoO1lbgdknXAdvmJkbEbw8klZmZ2eA8M/19Tts0X6Z9ybzXbqPKddcGq9cG1lfSw8zMbKSlLu9m+2RYNyxdCp87Z1YOCzawJB0bET+OiCuGFcjMzGzQJL0SOAWYmpsWER8oLtEoGo+uPKPFDSgrh1LUxL3oTTusHyEWuw/WX88NSPrygLOYmZkNnKRPA28A3kG+iX49cFyhoUZSKXaxrI3c6F0aF5cNyGINrPaqd8Igg5iZmQ3J8yLizcCmiPg94LnA0wrOZLabvWm+lrWLoO9rVkJl7OM6RhZrYEWXYTMzs1G1I/3dLulIYBY4osA8I2q8fv7XSOzl7E2Zj9fnZDYKFrvIxc9J2kL+7VyRhknjEREHDDSdmZlZ/10raRXwB8At5D8gfqbYSKOifWd9QC2SPl+lfLzsfcG4mVUM31KuHIb9MSzYwIqI6rCCmJmZDUNEfDANflnStcCU7/M4jxs4A+XeWeXWz4+nbJ91yeKMrcW6CJqZmY0FSc+S9NS28TcDVwMflHRwccnMzEZLWc/3Kws3sMzMbLn4M6ABIOmFwEeAK4HNwKUF5rIBGNXdv3BnvlLoR9c+f5ZLF2PSp7LXGw2bmZmNumpEbEzDbwAujYgvk3cVvLXAXGY2psakvWBL5CNYZma2XFQlzf2w+FLgm23z/INjiZXtPJZcKUOVVD/Kqn/l3cszlbPO2ajwBsXMzJaLLwL/JOkx8ku1fwtA0s+SdxM0M7NelKABujcHB4fVcHYDy8zMloWI+JCkb5Df8+ofInZuaivAO4pLZjY4JdgPLi333rNBcQPLzMyWjYi4scO0HxaRxczKa9zPnZKb3gPlc7DMzMzMzMz6xA0sMzMz68lwLqE85ocOCuJSteVs2PXfDSwzMzMzM7M+cQPLzMzMduOzM6x8llorux+zWOxKcuN8ifYxfmul4gaWmZlZjySdLekuSeslXdxh/qSkL6X535G0evgpzdqM+9UaykzCnTOXJzewzMzMeiCpCnwKeAVwMvBGSSfPW+wCYFNE/CzwceCjw01pO41qw2JEY5vZLm5gmZmZ9eYsYH1E3BMRDeAq4Nx5y5wLXJGGrwFeKo3env4TD/+Ejdsae0yfnH5s5/BEYzD3Zla09pi2ZXp2Sc+xcVuDicYTS1pnavqne0x7crrZdfntjT1z9mqh51VkADyxvcHGbY3dPoep6UeX/FrV1jQA22b3Pu9iVuzYs+zatdebvbXU915foH7uaDTZuK3B5h2dP4etM0129KG86o0tbNmxtLrbqydn9u55W1neSbDa2rFPr795id/JXqzY/iCTM4/vHH9820zH/0NzmllvHR43bmswm5ad3v7kvoXske+DZWZm1pujgJ+0jT8APLvbMhHRlLQZOATYbQ9T0oXAhQDHHnvsoPLutS333cyWDtNXTD/cNja8szkeeXJmScvf89hWqkt8jYM23b7HtIefnO66fK87d508saP7TuOcH2/cvse0/bbeu9ev2dqHvItZ9cQdC85fuePBPrzGnUtavhLdG7GzWXDPY1u7zt+0ffHPpxciY8t01pfnmu+hzd3r5jAs1PDZWwdv/N5u4/c+tq0vz9v+We/Y/BirDlzVl+ddiBtYZmZmQxYRlwKXAqxZs2af9nyf9oLXMjO9g6nJOtsbGSuYZbY6RbU1g8iYra5ErWmmZ4N6vU5dGS2qTNQrzDYa1CZW0Gy1qNGiQZVqq0F9aiVkTZpRIZvdzmxLTE1OUqnVmJ3eRkt1Jup1Zqe3UZtYgSoimg0iy6jWalQqYsdMg8laFRCq1ZmZnmaqXiELAUGmGhmCrEW1VqeSNQlVyBCVaNFoZkzWq8TsNJWJSWJ2llZkVCpVqNRAFSq0aGViNipMMINqE0SrRaVWpxWiNTtNNTKo1pltzTJZrZBVJ8l2bIZKjXpVRAbZxEpaGdQrGdNNmIhZFE0qE1O0WtCKjGpkNCsTVKNJJWtQrU+SqcJsM6NGCyo1KsrI0q6VoslME6BJpVKnXoHImoCo1OpkWRCVKrOzs0CFCk1qE1P5hzo7TaU+BdHKuzpGxnRTTNRrqDlDVGpEYxuqr2A2g0arxWS1wkS9RmN2lnq9Dq0mWZaRVSepRQtVRZZBs9mgQn6jWdVXQKsBBM2oUK3VqDa3k1UmaLVahKr5DWkrNWoEM60WtVqdSmuGSrVKRoVmK8s/PwUxvQUm9kfZLK3KBNnsNKpUiUqdagVaszPUJyZp7thKpVpH1Wr++hNPoTk7S6OVsaKSUalWaWXBTAYrqxVakdFUnUq0qFcqNJvTVKt1iIxKBZrNFlQn8rqdBVGbyi9UkTWZnd7K5FMOpDEzzWStRksV6swS1an8fc82aAbUCKrVOlmlQpUMRYvZLB+OSo1mYwdNqtQV1Ko1kJhpNpms16FSJVpNIKiQ0cgqhKpM1quQNZlpNFFF1Cp5x7EKLWaZoFbJaGYiazao1/I6nf9wIZqtZv45ZhkzmaA1ndejWoXIgqqCFlWarSaVbJaZqLJyos5MM2Nqok42O00zQJUaFaCqJllUqFQrEMF0VqEm0Wo2kESWZUxOTDDbbOY5IqNWrTDbgnqtirIWqkAzg2YrWFELslZGC8giY2Jikmi1aIaIyAjVqWXTiIyZqDI5MUXWbNBE1CQq1bzMZrIqU5UmUZskCxERVKOVf//qkyibJSKYmW0SAZPVCk1VqdWqCJhutFC1Rl15Q1bVCs0WzDZ2UJ98CrVooGqdQESrycr9B9+4gjFoYEk6G/gjoApcFhEfKTiSmZmNpw3AMW3jR6dpnZZ5QFINOBB4nAE6YNUhO4f3H+QLzTnw4LaRQ7ouNpQse+ugQ7vOOmCIMZZiz/I8rIAUi+lerrtZ1eNyfXH4EF+rszJ+F8qYaSHd8pb1fYz0OVg9nnBsZmbWD98FTpJ0vKQJ4Dxg7bxl1gLnp+HXAd+MGOeLPpuZ2XyjfgRr5wnHAJLmTjj+fqGpzIZs3bWXFh3BbOylc6ouAr5O3mvi8oi4U9IHgHURsRb4c+DzktYDG8kbYWZmtoxolH9Yk/Q64OyIeFsa/xXg2RFx0bzldp5MDJwKLHw2ZvkcyrwTpEfEKOYexcwwmrmdeXhGMffTI6KsvT/6StKjwP37+DSj8Bk7Y/+MQs5RyAijkdMZ+2dfcx4XEYv2zx31I1g9aT+ZWNK6iFhTcKQlGcXMMJq5RzEzjGZuZx6eUcwtaV3RGYall431YkbhM3bG/hmFnKOQEUYjpzP2z7ByjvQ5WPR2wrGZmZmZmdlQjHoDq5cTjs3MzMzMzIZipLsIdjvheJHVRvFqAKOYGUYz9yhmhtHM7czDM4q5RzFzkUahvJyxf0Yh5yhkhNHI6Yz9M5ScI32RCzMzMzMzszIZ9S6CZmZmZmZmpeEGlpmZmZmZWZ+MXQNL0sGSrpN0d/p7UJflWpJuTY+1bdOPl/QdSeslfSldPKMUuSWdLulfJN0p6TZJb2ib9zlJ97a9p9MHmPVsSXelMrq4w/zJVHbrU1mubpt3SZp+l6SfH1TGvcz9LknfT2X7DUnHtc3rWF9KkPktkh5ty/a2tnnnp/p0t6TzS5T54215fyjpibZ5hZRzeu3LJT0iqeN98pT7ZHpft0k6o21eUWW9WOb/kLLeLunbkn6ubd59afqtw7wkeg+ZXyRpc1s9eG/bvAXr1nJUdJl0qkfdtmfD/A51qmf9zCXpzPS+16d11aeM75e0oa3+n9M2r+P2s1sdUB/2ZyQdI+kG5dvGOyW9M00vW1l2y1ma8pQ0JekmSf+aMv7eQs+rvdiX6pa9Dxk77mMW9Xmn56lK+p6ka9N4acoRgIgYqwfwMeDiNHwx8NEuy23tMv1q4Lw0/GngN8uSG3gacFIaPhJ4CFiVxj8HvG4IOavAj4ATgAngX4GT5y3zW8Cn0/B5wJfS8Mlp+Ung+PQ81SGVby+5XwysTMO/OZd7ofpSgsxvAf6kw7oHA/ekvwel4YPKkHne8u8gvzhNYeXc9tovBM4A7ugy/xzga4CA5wDfKbKse8z8vLkswCvmMqfx+4BDS1jOLwKu3de6tRweZSiTTvWILtuzYX6HOtWzfuYCbkrLKq37ij5lfD/wnzss23H7uVAdoA/7M8ARwBlpeH/ghylL2cqyW87SlGd6f/ul4TrwnfS+Oz4vS9yXWih7HzJ+jg77mEV93ul53gX8b9L2okzlGBHjdwQLOBe4Ig1fAfxCryumVvRLgGv2Zv19tGjuiPhhRNydhh8EHgH2+QaVS3QWsD4i7omIBnAVefZ27e/lGuClqWzPBa6KiJmIuBdYn56vFLkj4oaI2J5GbyS/r1qReinrbn4euC4iNkbEJuA64OwB5Wy31MxvBL44hFyLioh/BjYusMi5wJWRuxFYJekIiivrRTNHxLdTJihHne6lnLvZl+/DuCprmXTbng3tO9SlnvUlV5p3QETcGPme2pXsxb7CEr8L3bafHetAv/ZnIuKhiLglDT8J/AA4ivKVZbec3Qy9PFOZbE2j9fSIBZ53qftS+/z/YIGM3RTyeUs6GnglcFkaX+jzGXo5whh2EQQOj4iH0vBPgcO7LDclaZ2kGyXNfQiHAE9ERDONP8DCX9B+6jU3AJLOIm9Z/6ht8ofSIdqPS5ocUM6jgJ+0jXcqo53LpLLcTF62vaw7KEt97QvIf1mZ06m+DFqvmV+bPvdrJM3deLuosu75dZV3wTwe+Gbb5CLKuVfd3luR9Xop5tfpAP5B0s2SLiwoUzfPTV1UvibplDRtVMp5mMpQJp3qUbftWdHfoX7lOioNDyrvRel/+uXadbrAUjP2fX8mda16JvlRjdKW5bycUKLyTN3abiX/gfw68n24bs+71H2pvnyP5meMiLly7LSPWdTn/Qngd4EsjS/0+RRSjiPZwJJ0vaQ7OjzmH5EIure8j4uINcCbgE9IOnFEcpN+Afg88NaImKtclwDPAJ5Ffkj2vwzmXYw/Sb8MrAH+oG3y0OtLj/4WWB0Rp5H/s75ikeXL5DzgmohotU0razmPNEkvJm9gtf9feEFEnEHedfDtkl5YSLg93UJeD34O+GPgrwvOYwtbsB4ttj0rSllzAX8KnAicTn4awB8WGycnaT/gy8DvRMSW9nllKssOOUtVnhHRiojTyXsTnEW+31Yq8zNKOpUS7WNKehXwSETcXFSGXoxkAysiXhYRp3Z4/A3wcGqAzDVEHunyHBvS33uAfyT/teNx8sObczdgPhrYUKbckg4A/g54TzocO/fcD6VDtDPAZxlc17sNwDFt453KaOcyqSwPJC/bXtYdlJ5eW9LLgPcAr0llCXStL4O2aOaIeLwt52XAmb2uOyBLed3zmNc9sKBy7lW391ZkvV6UpNPI68a5EfH43PS2sn4E+CuG1113QRGxZa6LSkR8FahLOpSSl3NBCi+TLvWo2/as6O9Qv3JtYPfutn3LGxEPpx3cDPgMu76XS83Yt/0ZSXXyRssXIuIraXLpyrJTzjKWZ8r1BHAD8NwFnnep+1J9/R61ZTx7gX3MIj7v5wOvkXQfefe9lwB/RNnKMfbixLIyP8iPOrSfePmxDsscBEym4UOBu9l1EuNfsvtJcr9VotwTwDfIf5mZP++I9Ffkh04/MqCcNfKTFY9n18l/p8xb5u3sfkLh1Wn4FHY/ofAehneRi15yP5P8cP1JvdaXEmQ+om34F4Eb0/DBwL0p+0Fp+OAyZE7LPYP85HgVXc7zcq2m+8UXXsnuJ/PeVGRZ95j5WPJ+5c+bN/0pwP5tw98m34iWIfNT5+oF+Ub8x6nMe6pby+lRdJl0q0d02Z4N+zs0v571Mxd7nqh/Tp8ytv9P/4/k54hAl+3nQnWAPuzPpPd3JfCJedNLVZYL5CxNeZKfMz93YbIVwLeAV3V7Xpa4L7VQ9j5k7LiPWdTn3Zb3Rey6yEVpyjEixrKBdQh5I+Ru4Pq2D3QNcFkafh5weyq024EL2tY/IX3469OHNVmi3L8MzAK3tj1OT/O+md7LHcBfkK4CM6Cs55BfoedH5EfSAD5AftQHYCqV3fpUlie0rfuetN5d7OWVYwaY+3rg4bayXbtYfSlB5g8Dd6ZsNwDPaFv3V9NnsJ68O2kpMqfx9zPvR4Aiyzm9/hfJu5DMkve5vgD4DeA30nwBn0rv63ZgTQnKerHMlwGb2ur0ujT9hFTO/5rqz3tKlPmitjp9I22Nw051a7k/iiyTbvWI7tuzoX2HutSzvuUi3zbfkdb5E9ojsXsSAAAAoUlEQVR+LNrHjJ9PGW4D1rJ7A6Hj9rNbHaAP+zPAC8i7/93W9n/knBKWZbecpSlP4DTgeynLHcB7F3pe/v/27RgFYBCIoqBH8/6XWrvEImnMBw3MHEBEm32CC7PU294De3ycMXfd97RWb3dgHXOOVXW9EgIAAPDRL/9gAQAAnEhgAQAAhAgsAACAEIEFAAAQIrAAAABCBBYAAECIwAIAAAgZc6d3NSfpLKIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x117820ac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with model:\n",
    "        trace = pm.sample(40000, chains=2)\n",
    "        pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imported from the \"Wet Grass\" example to determine the conditional probabilities from the trace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "## From \"Wet Grass Example\"\n",
    "def calcCondProb(trace, event, cond):\n",
    "    # find all samples satisfy conditions\n",
    "    for k, v in cond.items():\n",
    "        trace = [smp for smp in trace if smp[k] == v]\n",
    "    # record quantity of all samples fulfill condition\n",
    "    nCondSample = len(trace)\n",
    "    # find all samples satisfy event\n",
    "    for k, v in event.items():\n",
    "        trace = [smp for smp in trace if smp[k] == v]\n",
    "    # calculate conditional probability\n",
    "    return len(trace) / nCondSample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(R = true | T = true) =  0.06595744680851064\n",
      "P(R = true | T = false) =  0.0010624841892233746\n",
      "P(R = true | F = true) =  0.017130620985010708\n",
      "P(R = true | F = false) =  0.0014590692673936416\n"
     ]
    }
   ],
   "source": [
    "print('P(R = true | T = true) = ', \n",
    "      calcCondProb(trace, {'Retransmit' : 1}, {'Timeout' : 1}))\n",
    "print('P(R = true | T = false) = ', \n",
    "      calcCondProb(trace, {'Retransmit' : 1}, {'Timeout' : 0}))\n",
    "print('P(R = true | F = true) = ', \n",
    "      calcCondProb(trace, {'Retransmit' : 1}, {'Fast Retransmit' : 1}))\n",
    "print('P(R = true | F = false) = ', \n",
    "      calcCondProb(trace, {'Retransmit' : 1}, {'Fast Retransmit' : 0}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now know `P(R = true | T = true)`.  Due to the small probabilities, the value of this fluctuates betwee 0.04 - 0.07.  When calculating the conditional probability using variable elimination, we get the value ~0.05.  This seems to be accurate to the sampling result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variable Elimination\n",
    "$$\n",
    "\\begin{align}\n",
    "P(R = true | T = true) & = \\sum_{F} P(R, F|T) \\\\\n",
    "& = \\sum_{F} \\dfrac{P(R,F,T)}{P(T)} \\\\\n",
    "& = \\dfrac{P(R) \\sum_{F} P(T|R,F)P(F)}{\\sum_{R} P(R) \\sum_{F} P(T|R,F) P(F)} \\\\\n",
    "P(R) \\sum_{F} P(T|R,F)P(F) & = 0.0019491314899390418 \\times (0.176734807584152976 \\times (1-0.0004108463434675431) + 0.0236974833272821081 \\times 0.0004108463434675431) \\approx 0.00034435682 \\\\\n",
    "\\dfrac{P(R) \\sum_{F} P(T|R,F)P(F)}{\\sum_{R} P(R) \\sum_{F} P(T|R,F) P(F)} & = \\dfrac{0.00034435682}{0.00034435682 + 0.00619382033} \\\\\n",
    "& \\approx 0.05266862798\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "From this exercise, it can be concluded that the sampling method provides an accurate alternative for variable elimination, which is computationally expensive and slow to do by hand.  In addition, it can be noted that networks with very small probabilities have a large amount of variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}