{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
\n", "\n", "# Running Parallel Chains\n", "\n", "Author(s): Paul Miles | Date Created: August 31, 2018\n", "\n", "This demonstration was made using version 1.6.0 of [pymcmcstat](https://github.com/prmiles/pymcmcstat/wiki).\n", "\n", "Determining whether or not your chains have converged to the posterior density is a significant challenge when performing MCMC simulations. There are many diagnostics available for assessing chain convergence. The most robust approach is to use the Gelman-Rubin diagnostic, which requires several sets of chains for comparison. The Gelman-Rubin approach essentially performs an analysis of the variances within each chain set and between each chain set. For more details regarding this approach see\n", "\n", "- Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical science, 7(4), 457-472. [https://doi.org/10.1214/ss/1177011136](https://doi.org/10.1214/ss/1177011136)\n", "- Brooks, S. P., & Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of computational and graphical statistics, 7(4), 434-455. [https://doi.org/10.1080/10618600.1998.10474787](https://doi.org/10.1080/10618600.1998.10474787)\n", "\n", "While this diagnostic approach does not require computations to be run in parallel, it is certainly more convenient if you can generate a set of chains by running parallel MCMC simulations. In this tutorial we demonstrate how to use the `ParallelMCMC` class of the [pymcmcstat](https://github.com/prmiles/pymcmcstat/wiki) package." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.9.0\n" ] } ], "source": [ "# import required packages\n", "import numpy as np\n", "import os\n", "from pymcmcstat.MCMC import MCMC\n", "from pymcmcstat.ParallelMCMC import ParallelMCMC\n", "from datetime import datetime\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import pymcmcstat\n", "print(pymcmcstat.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define Model and Sum-of-Squares Function" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# define test model function\n", "def modelfun(xdata, theta):\n", " m = theta[0]\n", " b = theta[1]\n", " nrow, ncol = xdata.shape\n", " y = np.zeros([nrow,1])\n", " y[:, 0] = m*xdata.reshape(nrow,) + b\n", " return y\n", "\n", "def ssfun(theta, data):\n", " xdata = data.xdata[0]\n", " ydata = data.ydata[0]\n", " # eval model\n", " ymodel = modelfun(xdata, theta)\n", " # calc sos\n", " res = ymodel[:, 0] - ydata[:, 0]\n", " return (res**2).sum(axis = 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define Data Set - Plot" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAENCAYAAABNdxylAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd5hU5fXHP4etwwKLfa0RUKMxYkMUjV0iKsbYC1lLVNZoAmogaPxpDFYEe9QAagyrxl5JiMFuEEWMYotGZS2JLDa6y9bz++Od2b07O3V3Zu6d2fN5nnmGeW+ZM3eH95457znfI6qKYRiGYRhGJunjtwGGYRiGYRQe5mAYhmEYhpFxzMEwDMMwDCPjmINhGIZhGEbGMQfDMAzDMIyMYw6GYRiGYRgZxxwMwzAMwzAyjjkYhmEYhmFknGK/DTAMIz8QkW2A7YENAQW+At5R1Q99NcwwjEAipuRpGEY8RGQ74CzgGKAqMhx+jkweS4EHgOmq+u/cWmgYRlAxB8MwjC6IyBBgCnAk0AC8BMwHPga+wTkZ6wJbAXsAewMh4BFgkqou9sFswzAChDkYhmF0QUQagbeBG4FHVHVNkv0rcFGO8cAPVLU8+1YahhFkzMEwDKMLIvITVX2im8ceoaqPZ9omwzDyC3MwDMMwDMPIOFamahhGtxCRShE5X0S28tsWwzCChzkYhmF0l/WBqcAOfhtiGEbwMB0MwzBiIiI3JdmlEldNcoaI7A+oqo7PvmWGYeQDloNhGEZMRKQNp3UhCXbzbldVLcq6YYZh5AUWwTAMIx6fAOsBvwMejbF9S+A54BfA33NmlWEYeYE5GIZhxGN74DJcnsX+wC9V9fPIRhGJzB9fqeqnPthnGEaAsSRPwzBioqoNqjoB2B3YDHhPRCaIiC2DGIaRFHMwDMNIiKr+C9gN+D1wKfC6iIzw1SjDMAKPORiGYSRFVdtUdRquJHUprjfJ9XQ0PDMMw+iEORiGYaSMqtap6sHAqbilk0QVJoZh9GKsTNUwjG4hIqW4DqprVLXFb3sMwwgW5mAYhmEYhpFxbInEMIyUEZFiERngKVE1DMOIiTkYhmEkREROEJHZIrIUaASWAY0isjQ8fqLPJhqGEUBsicQwjJiISF/gCeAA4DvgTeB/wFqgHNgU2AmXh/E8cLiqfueLsYZhBA4LcxqGEY/JwD7AOGCmqjZG7yAiZcBY4FqcTsbEnFpoGEZgsQiGYRgxEZHPgAfCap7J9r0WOE5VN8++ZYZh5AOWg2EYRjw2AP6d4r7vAetn0RbDMPIMczAMw4jHJ8CoFPc9NLy/YRgGYA6GYRjxmQEcLSIPiMieIlLi3SgiJSKyl4g8CPw0vL9hGAZgORiGYcRBRAS4AfhleKgN+BpXqlqGWxKJ/Ei5BRivNqEYhhHGHAzDMBIiIj8ATgKGAZsAfXFlq18ArwH3q+o7/lloGEYQMQfDMAzDMIyMYzkYhmGkTVgu/E4R2dZvWwzDCCbmYBiG0R1CuJbtm/hsh2EYAcUcDMMwDMMwMo45GIZhdBdL4DIMIy7mYBiG0V3EbwMMwwguVkViGEbaiEgfYHOgPlYTNMMwDHMwDMNIi3AH1fWBr1S1yW97DMMIJrZEYhhGSojILiLyLLAK+Az4UXh8QxF5RkQO8tVAwzAChTkYhmEkRUR2Al4ChgCzvNtU9Utc2eopPphmGEZAMQfDMIxUmIyTBt8euICuCZ7PAMNzbZRhGMGl2G8DMomItOCcppV+22IYWWIA0Kaquf6/uzdwlaquDudgRPMZvVh0y+Yeo5eQ1vxTUA4G7j+4VFZWVvptiGFkgxUrVoA/kcdyYEWC7QNyZUhAsbnHKHjSnX8KzcFYWVlZWbl8+XK/7TCMrDBw4EBWrFjhx6/kj4FdE2w/AHgvR7YEEZt7jIIn3fnHcjAMw0iFe4HqqEoRBRCRXwOjgFo/DDMMI5gUWgTDMIzsMA0YCTwFvI9zLq4XkQ2AKmAucKt/5hmGETQsgmEYPlFfD/Pnu+egExbUGglMABqAtcA2wNfAb4DRqtrmn4VGutSvrmf+5/OpX50HX0AjLzEHwzB8oLYWBg+GUaPcc20eLC6oaouqXq+qw1S1QlX7quqOqnqtqrb4bZ+ROrWLahl842BG3TOKwTcOpnZRHnwBjbyjoKTCRWS5JVoZQae+3jkVDQ0dY6EQLF4MVVWJjw0nWa1Q1YHZtdJIh3yae+pX1zP4xsE0tHR8AUPFIRaPX0xVvyRfQKNXk+7842sEQ0T2CEsMrxaRFSLymIgM8dMmw8g2dXVQUtJ5rLTUjQcdESkWkQEiYvlbeUrdsjpKijp/AUuLSqlblgdfQCOv8M3BEJHdgBeAzYDfAZcDOwIvichGftllGNlm0CBobu481tTkxoOIiJwgIrNFZCnQCCwDGkVkaXj8RJ9NNNJg0DqDaG7t/AVsam1i0DoB/QIaeYufEYzJuKZJe4TXcKcCewL9gQt9tMswMkasRM6qKpg+3S2LVFa65+nTky+P5BoR6SsiT+NKVPcDPgQeAu4OP38YHr87HIns65OpRhpU9ati+ujphIpDVJZVEioOMX30dFseMTKObzkYIrISmK2qJ0WNzwZ2UdW0ZYfzaR3UKHxqa6Gmxi2HNDc7J6K6umN7fb1bFhk0KHXnIpc5GCIyDRgHnA/MVNXGGPuUAWOBa4EbVXVitu0KIvk499SvrqduWR2D1hlkzoWREunOP36uo5bhyt2i+Q7YWEQ2VtUlObbJMDJCfb1zLhoaOpI5zzgDdtwRhg51r6uqYjgWjz4KL74I11+fU3vjcBxwk6r+Id4OYafjZhHZEjgB6JUORj5S1a/KHAsjq/jpYHwAjBCRPpH6eREpBXYPb98E6ORgiEiynwfWB8AIBJFETm+lSFMTDB8OU6fCsGFRkYvWVrj4YrjqKvd6u+1g7Nic2x3FBsC/U9z3PWD9LNpiGEae4WcOxq3AdsBMEfmBiPwQmAVsHN4e8s0yw+gG3nyLWImcAI2NMG4cHHywR//i22/hsMM6nAtwEYxYJ8gtn+AkwFPh0PD+hmEYgI8Ohqr+EbgSqAbeBd4GhgDXhHdZHeOYgYkeJO72aBhZI1o4a+5cl3NRWhp7/1WrXHTjD2cuomXnYfDUUx0bDz4Y5s3rWsuae2YAR4vIAyKyp4h0MkhESkRkLxF5EPhpeH/DMAzAZx0MVb0I2AjYGxiqqruFbVJc90bDCDzefIuVK91zTQ2MHAmvvQZlZbGPO5F7ea5xBMWfdegPXFP8W+4+8a+w7ro5sj4hNwA3A0cDLwHficgSEflERJbg8qVeBI4CbgnvbxiGAQRAKlxVl6nqP1X17fDQQcACVV3lp12GkSqJhLOGDoWZM10pav/+blsxzVzPudzLGPqG85xX0p8jeYRJLVcw9hdFgehPoo7xwA7AVcAzwFdAS/j5GVwUckdVHaeFJAtsGEaPCZQan4gcD+wGmHCPkTckE86qrnbRjLo6ePfZpWxzyfHs0/ZC+74f9NmWI9oe5QO2BaBPH1i4EEaPztUnSIyqvgf8n992RBMukZ2MW2ZdB1gEXKSqzyQ57ijgeGA4LoL6GfAkcLmq2jKrYWQIP5U8DxCRp0XkNyJyuojcDtwD3KOq9/lll2GkS7RwVnk5TJjQdZ8RRQs447ZdOzkXa0f9lH1KX213LgDWrIFjj82PBmg+cxdwHk74azzQBswRkRFJjpuBSzCvxel8PBV+nici5Vmz1jB6G6rqywPYGvgHLtS6FpfkOQ7o04NzLq+srFTD8IMlS1Qvvli1vFx1wADVUEh11qzwxhkzVEtLVcE9RFSvuEK1tVVnzXLHRDZFHqGQO6eXyspKBZZr7v6fro9r0X4FsJdnfBLwKU6N9x/ADrmyKfz+w3G5Wud6xsqBj4AXkxy7X4yxk8PnO7Wb9tjck4Qlq5boy5+9rEtWLUm+sxFI0p1/rJuqYWSIWF1SK8sb+eKYcfS921Ngsc46cO+9ruQkzOzZcMIJLnrRfmwlzJkDIzy/x3Os5FkFLMRp0oC7Af8M53RMBd7AlZP/EFfBtaOq/jfbdoVtuwY4F1hXVVd7xi/EOUObahpCfSLSH1gJXKOqk7phj809CahdVEvN7BpKikpobm1m+ujpVO9YnfxAI1DkVTdVwygkopM9N+W/zG3et7NzMXSoS7AY1VleYtgwaGvrfL4ANECbiBOvOw7YA5fjcDlwKrCbqo5Q1Z2AkUA/XKQjV+wMvO91LsIsAATYKc3zRSTPvu6pYUZn6lfXUzO7hoaWBlY2rqShpYGa2TXUrw5AJrORVczBMIwM4U323JsXeZ1d2a311Y4dTjqJ+kfnM3/p4C5VIgFtgHYocKeqPqSqC4ALgEHAE9pR9YWqPgfcB/w4h7ZtTJTSb5jIWLq9jCYBrcAjsTaKyPJED0xFOC7WHr73Yg6GYWSIqiqY/kfl/JKbeIYD2Ygv3YaiIrjuOmoPvpvBP+zbLsYVncRZXQ2LF7tlkcWLOzdG84ktcLlREd4NP78RY9/Xge9l3aIOQrjW8dGs9WxPCRE5CTgdtzySd/o79avrmf/5/MBGBKw9fO8lUGWqhpHXfPcd1XNrqG6+u2Nsgw3ggQeo33Y/agZ3bn4WEePyRiliNkDzj2Y6zxGRm3cXld3wNsm6RR004BomRlPu2Z4UEdkbuAP4K3BxvP2SrTn7FcXIh9yGSHv4mtk1lBaV0tTaZO3hewnmYBhGD4i0XN+qqI4Nao6CN9/s2LjbbvDww7D55tTN79r8LCLGFSCHIpovgM08r1cDvwLej7Hv93AVYbliCR19i7xExr5IdgIR2RF4AngLOF5VWzNnXvbx5jY0tLgvVs3sGkYOGRm4m3f1jtWMHDLS2sP3MmyJxDC6SaT/yFUHzKXP7sM6ORffnfhzXrnmRepLNgeSi3EFlH/hkjsB15pdVW9R1Vg37wNwSaC54k1gWxHpFzUe6cac0BYRGQL8HfgSOExV1yTaP4jkW25DVb8qRmw+wpyLXoQ5GIaRBpGOqW+9BTVjlXENV/Po2lGsx7cAtBaV8Mpp01n/0ds5+Ijy9lyLgCZxJmMyTgo8ISKyES6iMCvrFnXwEFACnOGxoww4DZgXcYJEZAsR2dZ7YLj89h84Ya6DVTUvK0d6Y25D0PNNjM6YDoZhpEhtrcubKCmB4oZV3N56Gke2Pdy+/X9swgnFD7OweA/Wru04LhRySZtVVR1LKoMGdc+5yKUORtARkQdwXVyvxzVHPAXXamB/VZ0X3ud5YF9VFc9xbwI74jo3vx112o9VdX43bPFl7onkYHhzG4KWg5Ep8iHfpNBJd/4xB8MwUsArorU1/+FRjmR73mvf/iJ7cxwPsKpvFSLJBbO6izkYHYRlvS/DiX+tg8ul+K2qPu3Z53m6OhiJJr0/q+qp3bDFt7mnfnV9wec21K+uZ/CNg9tzTQBCxSEWj1/c5TP3huvhF+nOP5bkaRgpUFcHxcVwOE9QSzWVrGzfdjO/5Hyuo4USytu6Hpss16KnUY3eiqquxYmBTUywz34xxnJZ7ZJ1qvpVFfyNNJJv4nUwIvkm3s9uUY5gYTkYhpECCxe0MWHVJTzBEe3ORQPlnMyfGcfNlPcrIRSCGTPcI9Vci0iiaDxtjKAhIn1F5EIR+aeIvC8iL4jIr8P5D4aRFVLJNzHF0M4EIV/FIhiGkYSlHyxnq/N/xiH8tX3sE77HUTzCG+xC//4wbRr85CcdjkSkPXuiqER9vcvpSKaN4RcishI4XVUfDL8eALwIDAWacKWgI4AfAUeJyP6q2uSXvb2dQl4aSEVLI9UoR28gKJEcczAMIxHvvEO/Hx/JIW0ftQ/N5SBO5C98w/oAtLR0di4gNcGsSO+SAGtj9MNVakS4HOdcXARMVdWWcOTicuDXuNbpU3JupZGzG4qfTkwyLY2gVNX4cY287wkERh/FlkgMIx4PPEDLrrtTsaTDubiaSYzi73zD+vTv37Ny0zzUxjgeuE9Vr1LVFmjXxpgIPBfebuSISAj8raVv5WRpoHZRLYNvHMyoe0Yx+MbB1C7K/XpeIi2NSJQjVByisqySUHEo54qhflyj6Pe89bVbA6OPYlUkRq8lbnJlSwtccAFce2370GoqOJW7eJhjKC+Ha65xHVB7mpgZKX0tLXXOxfTpiXuQ5LhdexvwM1W9V0QqgFXAyap6d4x9fw1cqqr9s21XEMn13OONWKxtWUuRFHVaGqgsq2TOmDmM2DwDpUukV8XhN35FWfy4RrHes7yoHATWtnTUykfb0d1rZFUkhkHyygyvpkVzs+fG/tVXcMIJ8Oyz7fv+h605kkd5j+2pqID77oPRozNjZ3V1avkaAaAJ1210ZZztq4Gi3JnTe4klER5NppcG8im/wa+qGj+uUaz3LCsuY9zwcUybPy1mvkou8zPMwTAKjrjOQ5h4yZWHbPg66489Cj77rH3fv/UZzUlttazAOextbS5ykUkC1uAsmrEiclD43w3AkDj7bQF8kxuTejexbiqhohBttFFeXJ6VZmJByW8IMn5co3jvefbwszl7+NldohS57l9jORhGQeF1HlaudM9nnOGkvSNEkiu9nMpdrPuTvdqdizaEy4p/zxOnP05TaGA+yXtnmn2AU8OPfsBxcfbbH3gnNyb1bmLdVBBYcOYC5oyZw+LxizP+izQI+Q1Bx49rlOg9Y+Wr5Lp/jUUwjIIiVmVGUxMMHw4zZ7pIRkUFNDa6bSU0cT3ncU7Dre37L6eSMdzD31oOo/TPToUzFAr8EkbGUdWUfoCIyLrAK8DzWTXIAOKXbA7daGhW37c3dETtaf6GH9conffMdZTFkjyNgsIr6R1NKARTpsCkSaAKA9cu4RE5hhH6cvs+7/XZniPaHuUjtm4fKyvrcE78xqTCg4kfc4+fJaOZfO+g6HcERTsi2/Skf431IjEHo9dTW+uWRZqiJJ/693djjY0wgpd5iGPYhCXt2z/c+Tj2/PcdfL02ugN454ZlfmIORjDpTXNPJm/EQbmp51OVTCbIVRWJ5WAYBUd1Nbz2mos8eGlshNIS5RfcyvPs1+5ctNKHiVzDD9+5j3P/rx+lpV3PGRHA6u2YVHjvJpNy3EGS9s51boLfJNITySTmYBgFydChblnD2xNk8oUN/OG7n3Mr51CKW4f8mvX4Mf9gGhNpahYuuwwuuaSrcxJwAaysICIrReRYz+sBwMvAFcAwoBQnFX4N8KyIxHDNjEIikzfiIN3UrUomO5iDYRQs1dVuWWPOHLh10qeMnLw3J7fd1b79dXZhV17nWQ5sH2tshCuugKlTU29YVsAkkgrvp6qDgf7AdThH47ycW2j0iHQbYmXyRpyJc2WqoZdVyWQHqyIxCpqqKij957NsdenxbMDX7eOz+pzCkotvo/6qkJOQ8tCnj4tWLF6cFwJYuaRdKjwyoKqNwEQR2SW83XqR5AndyX9IpelYqvT0XJnO3+gNVTK5xpI8jcJFFa69Fp00CWlrA6CZYs7lBu4ZcDZz/i5UVLgS1kjZaoTyctd2PQiVI15MKjyY5Nvc09OkRr+rSHpbUmZPyOTfyqTCjV5BMilwVq92pST334+Eh5ZQxTE8xMvsRai549iZM2HsWFjbId3P2rXBap0eAEwqvIBIJGsd2Z7ohpRJOe7unCufpMv9xO8qHcvBMPKO2lqndTFqlHu++WaYP985HQB89BGMGAH3399+zJdb78mPyl/n3cq9uuRUVFfDgw86AS4vVjkCOKnwO4HpmFR4wRAv/2HhFwt975iaCvmQlJmp/JCevL/fVTrmYBh5RSwp8HHj4OCDXUTi7hP/Stuuw+Adj2r1WWex4TvPMa9uE+bMcbkV0Usfw4a5PiNeemPlSAxMKrwAiZXUOOWgKUx6elIgykaTEfSkzCC0tg9ClY4tkRh5xcKFLgkzmtWr2vg/Luek+y6lD+G8orIyuO02OO00IHFTsaoqF9WIbp3em5dHTCq8sIlOasy3ZYegJmXmuqFYPIIQ5TEHwwgc8fIramu75koADGAFsziZI3iifexz2ZzQY4+w/qjUW5/mUev0QKGq3wLn+22HkT7R+Q9+35DSxY/W7MmSJoPiqGWy4qe7mINhBIp4rdYjSyPRzsV2vMdj/JRt+LB97Fn258x+93N35Qasn+b7B7x1umFkje7ckILSRyRVempvKkmTQYgcRPA7ymNlqkZgiNWoLNIDpK7OJXWu9NQwnFT2MHe2nUJZ85r2sWn8mgu4mtJQcafeIUmrTvIE60USTApp7kn1Jux3hUK69NTedEpje9JQLMhYmaqRt8RqtR6p5Bg0yEU0APrQyhVcxAWNHZpOTSV9OYM7eKLvCZRG5U/Ei4oYmUFENgf+DKiqHphsfyPYpLLs0JM8Az+iHpnIi0hn6cPvyEFQsCoSIzB4nYgIkUqOSBLmpuXfMLf4EC7wCkYOGULp669wzWcndKkSiVV1UlPjKWk1MkFfYL/ww+gFdLdCwa/qikxUVKS79JGrhmJBxhwMIzBEnIh4PUCqf/gGn2wwjANa5nYcdOihrnXqDjtQVeXkL7xLIJGoiBfTt8g4HwODgMF+G2Lkhu7kGfipyxDL3oaWBipKK+Ic0ZVMlsb6rZGRK8zBMAJDfT1stRW88gpd9Spqa2HPPSn+/JOOAy6+GJ58EtZZJ+45E0VFjMygqi2q+qmqfuq3LUZu6M7N1k9dBq+9oaIQAEVSxB6378HNr96c0s2+fnU9W627Fa+c8Qpzxsxh8fjF3cqrCIJGRq6wJE8jENx8M0yc6KILLS2ePInmZpgwAW66qWPn/v2dw3HEESmdO5KD4dW3yNccDEvyDCbZnnuCWq2Rjl1B6B/y1tK3GD5zOI2tnZsP9S/tT0tbS9xkzEwltAbhGvSEdOcfi2AYvnPzzU6Ns7ERVq3qyJP48u2lcNBBnZ2L7bZzSyIpOhfQuW17LBVPIz4iUiYiZ4vIrSLyOxHZOs5+B4nIs7m2rzcQ5F+86eQZdHeJIZPLCWua1lBWXNZlfFXTKhpaGjjjyTN4a+lbXd4/U0s72Y7iBG3pxapIDF+pr3cBimhGyCsM2O9o+PaLjsGjjoK77nIRjDQxfYv0EZG+wDxgKLT3jPutiFysqtdE7b4RsG8u7esNBEUVMlOkW12R6VLYWLkYXppamxg+czgzD5/Z/j6ZFM7KpkZGEMuGLYJhZIX6+qgGZFHbZs92j4UL3dKFlzOYyZzv9qE87Fy0Ibxx7JXw0ENJnYtE72ukzbnAjsCVOCfjMOB14CoRudVPw3oLQegnkWlSjXpkIynUG0XpXxp7Lmlsbez0Ppl0CrLVQyXVa5XrCIdFMIyMk0h3orYWTj+9I/GyuBgk/Nu4lEb+wC85k9vbz/Ut63Aif+Gl2QezeGniKITpXWScY4H7VfXi8Ot3ROQp4EbgHBEpUdUz/TOv8AmSKmSuyZbktjeKsvCLhUyYO4Gm1qZO+3jfJ9OS29nQyEjlWvkR4bAkTyOjJFLjBFe9ES33XVwM3yv6L/c1H82wtgXt42+yI0fyKJ8wiP794amnXBlquu9bSEsjuUzyFJGVwARVnRFj22XARcCfVPV0ERkDzFLVomzbFUSyOfcUqipkMnKVEBkr8TPW+wQ10RaSX6tMXUtL8jR8JZHuRF1d7E6oBxa/wLvlu3RyLv5SNIY9eZlPcL/UVq1yyylevMshpneRFdYCJbE2hKMak4HTRORP2FySNap3rGbx+MU9Ko3MR3LVkn3oRkOZefjMpO8TZOGsZNfKr6U2i2AYGSW9CIYynhuZxgSKaXVDRUUwbRo39xnPuPHiPXWniET0csiUKTBpkkUwMomIzAc+UtW4dzQRuRS4BPgU2MIiGDb3ZJpUIwc9jTAEOUKRKvE+g18RDFS1YB7A8srKSjX8ZdYs1VBItbLSPc+a1XlbSYlqiDVayxhV6HhsuKHq88+rqurLL6v27995c0WF6pNPqi5Z4s7r3RYKqd50U/z3LRQqKysVWK65+f80GVgB9E+y3yVAG9CaC7uC+LC5x19mvTlLQ5eHdMBVAzR0eUhnvZmb//xLVi3Rlz97WZesWpKT9+sJkWtUeVVlt69RuvOPRTCMtEmlM2mifb56dTGlJx5FZd2ijsHdd3dVIptt1n58dCQEoLzcCXLdeGPnzqqVlU7nYtCgwuiaGo8cRzC2Bn4OPKiq/0qy7znAMFU9Ldt2BRGbe/zDL/GqIJaFJqOnUZp05x9zMIy06EmlRn09fPuXp9h28on0Wb6sY8OZZzq1rbLOAji1tTB2bNek0PJy9+wdL8TlkFiYkmcwsbnHP+Z/Pp9R94xiZWPHL47KskrmjJnDiM3jZIX3kHxX5OwuluRpZI14nUnfeiu59kTtLOW2za9k2/MP6XAuSkthxgz3KOuqrlddDQ8+CBVR/YjKylwUI15TNMMweg9+lPJ2N2myJzoUQVPpTAXTwTBSJlKpEb1sMXy4u+nHi2gs/XAlA047ld+3Pdo+9j82pfzRh1nv0N0TvuewYdDW1nmsqQnOPts9Cnk5xDDSJRICryitYE3TGt8TFnORoJlpnYpU6I5T05MllVjHZlpLIxvYEomRMvHyIryUl7uow7Bh4Zv+++/z3agj6fvp++37vMA+nNH/AWY9tVFKOROF1Kysp9gSSTAJwtwTuQmpKmtb1xIqdl1D/coNSPWGmslGYrm84aajT9KTJZVYx5ZICUVFRZQWleY0/8NyMMzByCrem31Dg6sqjXY4Kipc1OHBMY9x8L0nU/zdqvZtNzCeiUylJFTSXlqaSj5HKomlvQFzMIKJ33NPrJtQBD9yA1K9oeZ7LkOqTk1P8kRiHRtNrq6Z5WAYWcXbmfS112Lv07CmlQsbLuaw249sdy6aikL8vKSWSytv6ORcROdzxMvjqKpyKp692bkwjHjEygmI4EfvklRzFPK910qq4ls9yRNJ1qANgnvNfHUwRGRrEblfRP4rImtE5D0RuUBEumb8GYEhcrMfOtRFHUKhjkTMdfiW2YzmYi5v37+OLdmn+GXOXfiz9pbpw4aZ8qbRM8Kt5KeIyBci0iAir4jIgSkeu6mIPCAiy0VkpYg8JiJ52+Aj0U3Ij94lFaUVNLY0dhqLZUdv6bXSEzN9R+sAACAASURBVFXS6GPLi8op6dN58gzqNfMtyVNENgUW4IR8/gB8C+wNXAVsD/TSVfb8oroaRo6EuXPhhp+/xQMtRzKExe3bn+LHnMS9tJavx5o1nXuJNEfNh01NbgnEMFLkLuBo4AbgI+BUYI6I7Kuq8+MdJCL9gOeA/sAVQAtwHvC8iOykqsviHRtUvImOKDS0NnTKwcjlcsPNr97MxLkT218nsiNTCZr5oMLZkyZn0cfO/XhuTpNau00yJS7cf94BqSp3pfoAJgEKbB81/hDQDJR045ympucDs2apnlxyr66hs7zmlVygfWhpV9pcsqTrcYWuvJlpcqnkGf0ABNgCKA2/7uN9nWNbhofnj3M9Y+U4R+PFJMf+Bqc8urNnbFucozG5m/bkdO6JpyAZGV9Uv8gXhcmbXrlJuZROj9LLSnVR/aKEx6WriOnd3y8Vz3j25Or8fqiIpjv/pBLB+AVwkohcDMxQ1UxlhQ4IPy+NGq/HORitGXofIwPESrKsr4fXX21h2Wm/4c+t17fvu4p+PHXiXVz22NH0L4XGRpgwoes5I9EPS97MG9YF6oCRwLPABlGvc8kxuHni9siAqq4VkTuAK0RkY1VdkuDYV1T1Dc+x74vIM8BxOOnznNCdX96Jqi4i7cX9oH51PRPmdv2PXlZUxpqmNQmPTcdu7+dvammiVVtpbmtuTxStmV3DyCEjc3Ydsq3oGe/8fv6tUyWVHIyhwELgNuANEdkvQ+/9Qvj5DhHZUUQ2D7d8PhWYoqpt8Q81ckltrStPHTXKPdfWusfug76k31EjGedxLt7n+xzYbwGb/upoFi+GcePc+I03dhzrxZI38w5J8jpX7Ay8r6qro8YX4GzaKdZBItKHjjktmgXANiLSN5OGxqN2US2DbxzMqHtGMfjGwdQuqk16TP3qempm19DQ0sDKxpU0tDRQM7smEOJLdcvqKC0q7TKeyfyA6M+/tnUtzW2d11pjJTxmS6Qq23+PIP+9UyLVUAdwGPABLrLwEPC9VI9NcM7/A77DhTojj4sT7L88yUNtiSQxS5a4RmLRyxWJ9o9uLFZerrpX6QL9lM07bXiUI3QAy9uXQ+I1JUv1vY2u+LxEsh5uaeGA8OuNvK9zbMs7wFMxxn8QnkdOj3Pc+uHtF8bYdnZ425AY23o890TNc/awR2AfmZp/Uq4iUdW/4pIvJwEHAf8WkStEpCLxkQmpA54HxuKSte4Efi8iZ/XgnEYcYkUiwC11xJP6jqh3ejm17U6ebtqbLfgcgDaEi5lMdd9HaA5Vtst2xzrWKkWMDBECGmOMr/Vsj3cc3TzWMIw0SKuKRFVbgGki8mdctcck4FQRuUBVk8f3PIjICcB0YBtV/SI8/Eg4hDlNRO7XqGxuTSLuISLLgcp07OgtePuIRISxampg+fLEYleDBnVUe5TSyI2M56ym6e3blzGQMdzDc+WH8uD9HgXPqGMjWKWIkSEagFjl7OWe7fGOI91jbe4xjPTpbpnqOrjIw9a40tK7wu2af6WqceSXunA28LrHuYjwBC4PY8fwexgZIFYfkeJi1zSssbGz0zFyZIeTUFXlnI5Lx37Bvc3HsHtrR/Xf27IDJ1c8wgetWzFjOowe3fk9I8dGy3xbvoWRAZYAG8cYj4xFzysRvsVFL+Idq+FzZxyNyo/vqdR0aVEpr535GkM3GppxW2sX1TL2ybGsbe3cyjhbipHpKnqm0+MkXaXQ7iRtxrMnU51eI+df+MVCJj09qUuCa6qfLdckdTBEpApXEhZ5DKPDU1fcWugCYD9gvohcC1yg0f+burIR8FWM8UhQ3RqxZZBY0YTGRtekrNETLI4sYXidgOpB/+SkAcdS9GXHGsqDRcfz9dV3cOteFQkrQKxSxMgSbwLjRaSfdk70jHTPWxTrIFVtE5G3cfNYNLsDH6rqd5k1NTbp6CJEFC+jb5TJqjO6QySxMNq5gI4EykzfwGJ9vuj3ir6Jp9LDo25ZHVMOmsKkpyelpBnhTapMpyolnj09FRLzfuZB6wziwFkHdrItmmz9fbpLKjfxL3COhOC8/5eBV4D5wAJVXQUgIsXARGByeP8Lkpz3P8BIERmiqh97xk/EJZK+lcbnMJIQK5oQkev20mkJQxVuuQXOO4+ilhYAWihiElO4rvV8QpcIixcndxqqqsyxMDLOQ8AE4AycVg9hBeDTgHmRyKiIbAH0VdX3o469SkR21nCpqoh8HzgAuDp3HyH18sxcKl7Gutln+z2Tfb50owrR+085aArDNhnWLUeuJzftngiJRX+GCXtOiPt3iRA0Rc9UHIwZhJ0KVf1PvJ3C+RlXicgA3BJHMgdjKnAIME9EIkqeo8Njf1TVL1OwzUiDWNGEgQM7nI5OehUNDXDWWTBrVvvxX8v6HKf38xwHALGjHYaRC1T1VRF5ELhGRDYGPgZOAb6Hm38izAL2pXM57a3AmcDfwhHXFuB83NLI9QSQXLYkjyc7Xl5cnrX3TPT50o0qxNp/0tOTUlo6yIYjlyhSFW9pJdZnmDpvapei8EhX1bKiskAqeiZ1MFQ13YqORbjlj2TnfVFE9gQuBc7BlcDVARfinA8jC0SiCZHKkZEjXW+QW2+FqVOdXsX9Uz7hmYFHsdmX7TpEvC7DOK7oYRa3bNE+Zgmbhs+cDFwWfl4HF/U8VFXnJTpIVVeF9XyuBy7G6QE9h1MF/SarFveAnkhNp0P0zb6xpZGJe03k7N3OzurNy/v5Kkor+Gz5Z8z+z2yAtKIKPYlCZMuRixWpShSVifUZyorLGDd8HNPmT+tkWy6+E90l4+3aRWQd4GBVvS+jJ07tva1dewpEWq5HKke8nU0P5Gnu4wTWp2OevZPTOJtbaSspp6jI5W1EEjbjtVc3soOf7dpFZD1c3tRBqvqsiGyE+9V/kKrmWskzUBTi3ONXf4/aRbWc/vjpNKuLJBRLMSKScjJjJlrAZ/uzJ7Mx0XbAN4ci3fkn44mU4dLSnDsXRmrEKledMAHKSpWJTOUqLqQIJ6LaRAnjuZE/chYgVPaFu++G9dazhM1eyrfAIJycPzhnw/vaSJF8aM7l/dWdK3vrV9cz9smx7c4FQIu2UEwx5cXlKS0FZCIKkW0Z7kRRlsj2RMmp2ajiycbf1yo1ehmxylXXLV3Nzd/9nGN4sH3sCzbmaB7mFTpKqZqaOutcGL2LcGXYp57Xbd7XRmpku3dFpsmlvXXL6ujTp0+XTlRlxWXcd/R9rNd3vZRugrlaTuou8XI9Fn6xkANnHZh2cmpPyObfN+NLJH5SiGHKTFNf71Q8Iw7GVnzIY3Ik2+u77fv8k704lgep90gFlJfDjBm2JOI3fi6RGPFJde7JRPg+l+Ta3vrV9Qy6YVCXEtny4nLqxgen/DITROugRCIWufxupPv3TXf+SVkq3CgMIuWqoRAc13c2r7FbJ+fi1WHn8OOiZ9udi5ISuPhiF/kw58IwekbdsjqK+3QOHMdqzpVr4jUDi4TyvWTT3qp+Vcw4fAYl0vGepX1KmTF6RkE5F+CiLIvHL2bOmDksHr+YYZsMy+m1huz/fW2JpBdSPaaNIxdNpt+1v+8YLCuD6dPZ/ZRTWFwPC8O9Jm1JxDAyx8IvFrKqaVWnMb+1CxKFyHOpvxEhsryx8H8LQWDYJsMKzrmIEJ3rketrne2/r0UwegGdmpktXw5HHNHZudhiC5g3D045BXAOxejR7mHOhWFkhvrV9Ux6elKX8SkHTfHtBpqsHXgkYTJUHKKyrJJQcSgnWgtV/aoY/f3RjN5mdME6F9H4ca2z/Z4WwchD6uvdkkVFBaxZk7iiw1uSulXjuzw38EgGLP2wY4cDD4S//AU22CA3xhtGLyVW5UD/0v4M2ySWannmSFQhkIpmRNATJgsJP651Nt/THIw8I+IwqMLatS6XAmJrUnhLUg9teIi7OJV+Sz29CyZOhCuvdF3PDMPIKrHC0S1tLVkNgSerEEg1RJ7tsk2jAz+udbbe05ZI8givw7A2nGQd0bMYOza8BOJh4UIollauZhIPcSz9cM5Fa6gC7r8frrnGnAsjJUTkWhEZIyLbiYjE2L6piIT8sC1fyHUIPNnyhx82Gb0Lu7vkEbE0LCKsXevkvidPdo7GrbfCHVO+5uGmExnJ0+37fSRbMfBvj7L+fj/MoeVGAXAerokhQIOIvAX8C3gj/DgB2B/YzR/z8oNchsBTlcy2JRB/yQfRte5iDoaPRHIpUlXFjNVy3cvUqS6V4je/ge3W/ot5HMWWHh2kv/U5jJW33M0J+5mEgpE26wK7RD1+gWu/FHE8ctLmPN/JVQg8nQoBWwLxh2yJXAXFabElEp+orXWCV6NGuefa2uTHeDUsSku7bi8pcWkVx66dxTz26uRcXFVyCW2PPsEJZ5lzYaSPqi5X1WdVdZqqnqSq2wIDgROB/+BkxM/01UijE91Z/oinh2FknlSWsLpD7aJaBt84mFH3jGLwjYOpXZTCzSVLmJKnD0SraYJzGhYvTi2SUR/WqTjmGNdiPULfkmau43xqmv/QPraCAfyMu3kmdHjM86cbRTH8JYhKniLSD3gHmKaqf0i2fyES5Lkn1V+z+SZhnguyGQmY//l8Rt0zipWNK9vHKssqmTNmDiM2H5HgyPhku0maKXnmAZFcCi+lpW48FSI6FTNnOsekshK2LK/nw80P6ORcvMsP2I3XeLr8cKZP7+pAdCeKYhjRqOpqYBYuT8MIGFX9qhix+YikkYts/JrOZ7IdCUi0hNXdSFI8Zc5bF9zqS1TDHAwfiJVL0dTkxtOhutpFPV66Zj4fVe7CJov/2b7tseJj2LfsVU64eJuYMt/eipSVK91zTU3XShTDSJGvwNO8xuhCkJcfci0JHnRy4XDFW8Ka+/HcbjsDsZyWxpZGps6f6ovzaA6GD3hzKSor3XN0hKGT+mY8VKl6fDo7/HJfipYucWN9+rDqt1ex0QsP8M4n/Zg8ueO83nP2NIpi9C5E5CsR+YeIXC0ix4nIVp5tAhwGLPTPwmATpHXxWPghCR5kcuVwRfcjGTlkZI8cm1hOy8S9JlJa1DlpL1fOo1WR+ER1NYwcGTv/wau+2dwcW0SLtWvhnHPgzjs7xtZdF+67j/4jRxK9ghd9zilTMhNFMXoN84CdgIPCr1VEVgHvAxuEH6eISImqJqh16n14fw1H1sZrZtcwcsjInGb4J8oniNyYvN09e7MeRi4dLm8Fz/zP56dUWpyI6LJjgGkvT+u0T6zPko18E3MwfKSqKnbSZWTpIpIEWlPjnJH2fT//HI4+Gl57rePAnXeGRx6BLbfs8j6xzjlpknMyJk1ykYumpq5RFMOIoKo/BRCRgThHY+fw807A5ri55CGgRUT+A7wNLFLVKf5YHBxS1aPIJqkkcJoeRgd+OVyZcmyiy46TfZZsJfhaFUkA8FZy1NW5pMuVHYnFVFbCnDkwYgTw3HNw/PHw1VcdO1RXd6y5xGD+/PjnjLynVZHkBwGtIikFfkiHw7EzMBTop6pFftqWKxLNPYky+3NxE4/1/qVFpbx25msM3Who1t8/n/FDTyJys/c6A9nUxkjn+5nu/GMRDJ9JeeliS4XrrncqWq2tbkNxMVx3HfzylxBWb45VdpooqTRWFMUw0kFVm3Cqnv/yjovIEH8sChZ+Lz/EiqA0tTYxfOZwZh4+s9eWoqbiPPghQJatSFK8z5LNCJslefpIrEqOyNKFNwH0jpvWUHnOGPj1rzuciw03hGeegV/9qt25iFd2mkpSqWF4EZEDe3DsQQCq+nHmLMpvopP5cnlTjxV2B2hsbSz4UtR4lTtBT7pNpbQ4U2Qz38QcDB+JV8kxbJgrP50zB/714Mfsf9GehB79S/s+r8gePPTbf8E++7SPJSs7jZS0zpnjnrskjRpGZ/4uIs+KyGgRSbrMISIlInKkiLwI/C0H9uUd6d40MlXWGomgRFcSQGGXosZzIkzzozPZbHhnSyQ+kmzp4t1pc9j22pNYh4513ds4i3P1BvpcUMaPju+IQsRqhBYpO43sY8shRhrsDFwHPAF8JSJPAwuAj3Gy4ILrT7I1sAdwIE46/B+4PAyjB2Q66a56x2p2rNqR4TOH09jaIf9bqKWoiSp3gpB0GzSytSxjEQwfibt0sWEbqyZdzv7XHtbuXDRSys+5g7O5jSbKKCrqrFlRUdFZNhys7NToPqr6jqr+GNgL5zQcDlwPPA68BLwIPAZMA34cHt9DVQ9R1ff8sbowyNYv7KEbDWXm4TN7RWv2RDoWpvkRm2wsy1gEw2e66GH0XQlHnUz/xx9v3+dzNuNoHuY1hrePtbZ2OA+RRNFwKkZ7MYnlWRg9RVXnA/PDyyS7Aj/AaV4oTr3zHeANVW3zz8rCIpu/sHtLKWoiJ8LvpNvehDkYAaB96eLf/4b9joQPPmjf9jz7chwP8BUbto+VlsKMGe4Yb+5FhNZWJ5Ex1CrQjAyhqq24JZIFfttS6GT7F3a2KiOC0iIcklfu9BZHy2/MwcgQkQ6n4JI0044cPPoonHwyrF7dPvTewefxkxem0KeshPJG+MUv4IADOp8/Vu5FKARr1vTs8xiGkTuib8759gs7iJ1YkzkRfpSg9jZMaCsD1NbC6ad3JGyWlMAdd6RYqdHaCpdcAlde2TEWCsHtt8NJJyVtp97T1u9GfhFEoS2jZ3NPvJtzkCICifBbSCzo5MvfMRWsXXuOqa+HsWM7V4M0N7uxpJ1Jv/0WDjusk3OxZqNBvDBlPvUHnAQ4J2HEiPjOgmlcGEb+kiihM5daCD3BOrHGJx29jSB32+0u5mD0kIVx+kdGV3l0YdEit9bx1FPtQ0/JwWy+dCH7jduRLbboEMpKhmlcGEZ+Ugg351RyRgrx5pmMdKqBgi781V3MwegBtbVw7LGusWk03iqPLtxzjwtLeDyQq4su4lD9K8tYF0gjChImWaTDMIzgUQglk8mEmgr15pmMVJ3HQhb+Mgejm0SqN2I5F94qj040N8N558HPftaRNNG/Px9c9QiXlV1OG50FE5NGQQzDyGuyqaKYS+JJoRfyzTMZqTqPhRDFiodVkXSTeNUbkyc7/yFSQtqeoClLXRfUF17oOOD734fHHqNy4La0/b7reySMghiGURAUSslkrKqM3qyamWo1UCFEseJhDkY3iSXzDR3OhbdL6o6NC5hTcTQV3/63Y8cjjoBZs2DAAKpwEQ9vJUrcKIhhGAVHoZZMFsLNsydVIKk4j/lYlpwqVqbaAyJORGmpk+WePt0lWL71Fgwf7qS7T+d2buEcymiKGAmXXQYXXgh9Oq9Q9VhLwyh4glqmKiKtwP+A/1PVWX7bk2v8KpEPAsluwJEyXO/N02+NjFTJpb5HPpSzpjv/mIPRQ6J1Kmpr4YwzgKZGbuZXjGVm+76NfQeyZsa9rDvmkKT6FoYRiwA7GJ8AFcB6wJuquou/FuWW3upgpHoDjnXzDPoN1fQ9umI6GDnGW70RSfxcv+l/vMC+nZyLt9iB4bKQTc84hJ/+1DkWo0Y5kaxUy1ENI6io6paqugGuk+q9fttjOLJZHppOAme0pkc+VJYUcvJlrjAHI4PU1cE+8hKvsyt78Gr7+L2cyAjm89aaIaxdC48/7qpPVq50SaI1NamXoxpGkFHVt1R1mt92GNm/iXf3BpwvlSWFkD/iN+ZgZApVtn/mJp787gCqWApAC0VM6HMdNRX38B0VcQ8tLbVyVMMwMkcubuLdvQHnS2SgUEqI/cSqSNIgbt7Ed99BTQ0D7r67fegr2YCflTzA6Gn70Top8Xmbmqwc1QgWInJyd47rjQmeQSQX5aHdrX7Ip8hAoiqQoOeQBAFzMFLEW3ba3NxRMUJdHRx1FLz5Zvu+TTvtxme/e5g/77E5VVUwcKBT5YwlylVebr1DjEByF6CApHGMAuZgBIBc3cS7o+GRb2WZsUqIg9g9NohYFUkKxOpYWl4Oz134D4bfcCJ9ln3bseH00+EPf3A7RJ3j1lth6lQoK3MlrBMnwtlnm3NhpE6uqkhEZN/uHKeqLyTfq/AIYhVJ0MtD8zUC0JurS9KdfyyCkQJdVTuVcWuvYbff/ZY+tLmhkhK4+WYXqpCuP/qqqpzK59lnW3mqEXx6q6MQBDJ14w26Qmi+iov1ZnXSdDEHIwW8qp39WMWfOI1jeLh9+xdsQtkjD7He6BFJz1VVZY6FYRixyXToPV9v4kEmn3JI/MaqSFKgqsrlSexQ9h9elT06ORcv8SP26/86/1kvvnNRXw/z51spqpH/iMgwETlHRP5PRC6JelycY1sGisgMEflKRNaIyLMislMKx/URkdNE5EkR+Tx87Dsi8lsRKcuF7bHIl/LN3k6i6pLe2JY+ERbBSJHqyic4qaSaosaV7WM380t+zbUUNZfGrQKJmxxqGHmEiISAR4Af4xI/vQmg6hm7LEf29AH+CuwATAO+Ac4GnheRXVX14wSH9wXuBF4B/gh8CYzA2X4AcFAWTY+Lhd7zh1jLT5b42RWLYCSjrQ0uuQSOOIKi1c65aKCck/kz47iZZkqZODH2skdE2bOhwUS1jLznEpxzcQWwP86hOAU4BHgJeA34QQ7tOQbYEzhZVSer6i3Afjgn53dJjm0C9lLVEap6harOVNWfA78HDhSR/bJod1wqSitobGnsNGah9+DiVSe16FNszMFIxPLlcPjhrjlZmE/le+zFPGpxMgHl5S5xMxaR5FAvJqpl5CnHAA+q6iXAO+Gx/6nqU7hf/KXAqTm25wvg8ciAqn4FPAD8VERK4h2oqk2q+nKMTY+Gn7fLpKGpULuolj1u3wMJB4VCxSFCxSGmHDSFumV1vf5GFXTyRTws15iDEY+333YtTf/2t46xgw7i1T8s5P3QLlRWQiiUuKV6rJbuJqpl5CmbA5HKktbwcymAqrYAfwFOyKE9OwOva9c6+wVAf2Crbpwz8j/5654Yli7eX79rW51YTqu2ctHeFzHp6UmB7tdhOCzxMzbmYMTigQdgjz3gY88y7m9+A3PmcNzZ67N4McyZA4sXJ86niCSHhkK0OyQmqmXkKavoyNlaBbQBm3i2r6DjBp0LNgaWxBiPjG0SY1syfoP7HP+I3iAiyxM9gMpuvB8Q+9dvWVEZl714mYXc8wSTFY+NJXl6aWmBCy+EaZ5eTRUV8Kc/wbHHtg+lU2paXQ0jR5r2hZH3fAxsA6CqrSLyLm6Z4k4REeAo4PPunDicsFmayr6qGtHDDQGNMXbxbk/Hht/ilnpqVHVFOsf2lFi/fhtbGykrKqOxteMjWsJnsAm67ogfWAQjwtdfu/7pXudi663h1Vc7ORfdwdvS3TDylKeBo0WkKPx6OjBKRD4GPsTdnO/o5rn3ARpSeYjI+uFjGoBYJaXlnu0pISLHA5cD01V1Rqx9VHVgogcu8tEtYv36nTZyGi1tLZ32s5B78IluS9/bsQgG8PVTr9P/lKMoW/pZx+Do0a7GdGBWFZkNI1+4GqglXJqqqreKSDnwM1xOxkzgmm6e+33gtBT3XRV+XoJbJokmMvZFKicTkZG4/ilPAuekaEPGifXrd2D5wLzp12EYsfC1F4mI3IUrdYvHZqr6vzTOl3Y/gHlj/8yuM2so90ZbL70ULr4Y+liAxwgWuepFEnRE5EFcmepm3kRPEZkBnAisq6rN8Y4P77s78AzwJjBSVVOOesQ4V3b6IOVpvw6jMMm3XiTTcaFXL4ITv/kkHecibZqaWHPW+ez1p1vah5ZTyemld3NLzWiqwr5F3BbthmH4yUO4HJAjgMcAwssnxwKPe50LERkC4BXfEpHtcEJdnwCH98S5yCYm9W3kM746GKo6H5jvHRORH+GU9u7J2hsvWQLHHkvFvHntQ++wPUfyKF+FtqauzjkTpsJpGA4ReTaF3VRVD8y6MY6HcEqcs0RkGq609GxcXtmlUfs+E37eEkBE+gNPAesAU4HDpHODwrdU9a1sGW4YvQW/IxixOAmnxndvVs4+bx4cc0wnOc37OY7TuYM19CMU1qnwqnBGuqjW1LiKEItkGL2Qwbj/l16KcTkPfXA3+DW5MiZcyXIozkEYh6saWYBT9vwoyeHr4XQ9wOWWRPN7wBwMw+ghgXIwwup7xwEvq+onGT25Ktx2G4wf78pRgVb6cAFXM40JgFBe3qFTMX9+dIv2DhVOczCM3oaqbhlrPNwc7Hxckua+ObZpGXBG+JFovy2jXn9CRx8VwzCyRKAcDOBg3K+LmMsjYUGbRMQXu3nxRTinI0n8G9bjOO7nWVxEt6IC7rvPFY9AbBXOxkZT4TQML6raCFwlIj8ArsMlWBqGYQROB+MkoBnXTyCz7LsvjB0LwL9kF3bh9XbnAqA1LH4cWTmJqHB6e4m0tsLcuRm3zDAKgX/ifiAYhmEAAXIwRKQfLiP8KVX9JtY+PRW7qf/tTVxafDl76T/5jO+1jxcXO+dhzBgYPNgld9bXwzrrdK5UbW62bqiGEYdBpKjGaRhG7yBISyQ/JcvVI3VflHFd6CLWruoYC4Wcc9HU1LEkcvrpUFTknIvGKDFiy8MweiMiskWcTeviVDzHAc/nzCDDMAJPkByMMcBq4IlsvcHChbBqVeex1lYoK3MORoTm5q75FxGsG6rRS/mErlUkEQT4AOdkGIaRJ2RbyC0QDoaIbID7FfQXVf0uG+9RXw+TJnUdv+QSuOKK5MdXVEBbm3VDNXotk+nqYCjwLfAf4GlVbcu5VYZhdIvaRbXUzK6hpKiE5tZmpo+eTvWOmRV6CoSDARyPsyV7yyN1XctO+/eHAw6ALbZw+Z9r18Y+trzcVZgMG2bOhdE7UdVL/bbBMIzMUL+6nprZNTS0NNDQ4m6KNbNrGDlkZEYjGUFJ8hwDfElX2fCMEavstKXFjVdXw4MPuihFNGVlMGOGK18158IwDMPId+qW1VFSVNJprLSolLpldRl9n0A4GKo6QlU3UtXWbL1HpOw0FILKSvfsN92IrQAACvZJREFUXe4YNswtgXgpLYUFC0we3DAMwygcBq0ziObWzr+4m1qbGLROZhMMg7JEkhOqq53Ud6zmZREHpKbGORZNTe710KH+2WsYfpFi75FoctmLxDCMblLVr4rpo6dTM7uG0qJSmlqbmD56esYTPX1t155pMtEy2bqnGkEmV+3aReQTuiZ1VgDrh/8d+U8WseNrYLWqDs6mXUElW+3aDSObpFtFkm/t2gNHVZU5FoYR3b9DRAYDzwE3AlNUtT48XgVcgNOxseiFYeQRVf2qslKeGiEQORiGYQSe63FNCM+LOBcAqlqvqufiWqdf75t1hmEEjkJbImkDpLIyfs8zw8hnVqxYAS7XIac/DkRkBTBJVf8YZ/tZwNXZXroJKjb3GL2BdOefQlsiaQP6rFixYmWCfSIzQMK+JUbK2PXMLMmu5wDc9zzXKLBdgu3b58qQgGJzT+6x65lZUrmeac0/BeVgqGrSzxNp+d5bf2llGruemSXA1/MfwC9E5HWgVsOhTxER4GSgBnjMR/t8xeae3GPXM7Nk43oW1BJJKtiXMrPY9cwsQb2eIrIZ8BKwBbAU+DC8aWtgI+Bz4Eeq+l9/LAw+Qf3b5it2PTNLNq6nJXkahpGUsOOwEzAFWAYMDz+Whcd2MufCMAwvFsEweoRdz8xi17Nwsb9tZrHrmVksgmEYhmEYRl5QUEmehmFkBhE5OfzPWlVVz+uEqOqsLJplGEYeYUskRo+w65lZgnI9w7oOCoRUtcnzWhIcpqpalBMD85Cg/G0LBbuemSUb19MiGIZhxGJ/AFVt8r42DMNIlV4XwTAMwzAMI/tYBMMwjG4hIsXAEcC6wJPeHiWGYRgWwTAMIykicg2wv6ruFn4tuO6qe+PyMr4B9lDVj/2z0jCMIGFlqoZhpMIonJJnhMOBfYCpwEnhsQtybZRhGMHFlkgMw0iFzemQBwfnYNSp6gUAIrI9MMYPwwzDCCYFEcEQkd1E5BYReU9E1ojIZyJyn4hsleLxm4rIAyKyXERWishjIjIo23YHlZ5cTxG5VEQ0xqPXrs+LyDAReVREPhWRBhGpF5G/i8ieKR4fhO9nKdDieb0/8LTn9WJg45xaFABs7sksNvdkFr/nnkKJYEwC9gIeBN4CqoBfAm+IyHBV/Xe8A0WkH24tuT9wBW4SPQ94XkR2UtVl2TY+gHT7enqoAb7zvG7IuJX5wxDc/7WZwBJgIO7X/osicoiqzo13YIC+n58DI4CZ4WjFYOASz/YNgdU5siVI2NyTWWzuySz+zj2qmvcPYE+gNGpsa2AtcFeSY3+D62+/s2ds2/DFnOz3Z8vD63kpTpBpoN+fI8gPoC9QD8xOsl8gvp/hv2srMBuoA5Z7/8bAfcArfl9XH/6ONvcE53ra3JPaNc7Z3FMQSySq+rJ2CAJFxj4E3gW2S3L4MbiJ8Q3Pse8DzwDHZdrWfKCH1zOCiMiAcLWBEYWqfgd8hftFkYigfD+vAu7CRTEUOFlVlwOISCXwk7BNvQqbezKLzT3ZJ5dzT0E4GLEIf7k2Ar5OsE8fYCiwMMbmBcA2ItI3OxbmF6lczyg+A1YAK0TkThFZN2vG5Qki0l9E1heR74vIlcAPSXBTDtL3U1UbVfV0VV1PVQer6hOezatw+ReX5sKWoGNzT2axuafn+DX3FEoORizGAJsCFyXYZ12gDLc2Fc0SXH3/xoDV9qd2PQGWATcDrwBNwAG4NdFdRGR3VW3MqpXB5k/A0eF/NwF/BK5MsH8gv58iUgasD3ylqk2q2oab0A2HzT2ZxeaenuPL3FOQDoaIbAvcAvwTqE2wayj8HOuLtzZqn15LGtcTVb0xaughEXknfPzJuGSj3srvgenAZkA17j9wCbG/fxCw76eI7AJMA34EFAEjgWdFZEPgL8BVqvp0glMUPDb3ZBabezKGL3NPwS2RiEgV8FecN3ts+NdVPCLZxWUxtpVH7dMrSfN6xuOPuKzuAzNpW76hqm+r6lxV/RNwMLArLq8hHoH5forITjihrSFAp5bsqvolbrI5JRe2BBWbezKLzT2Zw6+5p6AcjHCy2RygEjhYk/dG+BbnocWq398Yl8wWK0TUK+jG9YxJeGL4Hy7sZgCq2gw8DhwlIvF+CQTp+zkZ+ALYHqfYGZ1A9wwwPEe2BA6bezKLzT3ZI5dzT8E4GCJSDjwJbAOMVtUPkh0T/vK9DQyLsXl34MNwxm2vozvXM8G5SnBKkF9lyLxCIYS7UfePtTFg38+9gZmquho3uUTzGbBJjmwJFDb3ZBabe3JCTuaegnAwRKQIuB9XQnesqr4SZ78twmt6Xh4C9hCRnT37fR+XIPRglkwOND25niKyQYxdJ+LCak9l2tZ8INY1EZEBwLHA5+ElhqB/P8tJnMg5IFeGBAmbezKLzT2Zxe+5pyC6qYrIDcB4nNf7QNTm1ar6WHi/54F9VVU8x/YH3gAqgGtxIiLn47y7nVT1m6x/gIDRw+v5HU506R1ciG1/XPbyP3HdOFvoZYjIs7jEqJdxAjebA6fhEq5OUNUHwvs9T0C/n+FkuZdVdayIrIf7RXiQqj4b3v5XYH1V3T0X9gQFm3syi809mcX3uScVNa6gP4DncWHbWI9PoveLcfxmOI9sBa6m/wlgsN+fKx+vJy5T+73wdWwEPsCt34f8/lw+Xs+fh6/Vl0Az7ub8ZPg/dJfrHuN437+fwG9xSV0HAevhFP72D2/7NU7l85d+X2sf/rY29wTketrcE/N6+jr3FEQEwzCM7CIipbgw8z7A+zjJ4LeBDXD9IuYCh2r3Mv0NwyhACiIHwzCM7KJOvnkkMAEXyViLS8L7GtezYLQ5F4ZheLEIhmEYhmEYGcciGIZh9BgR2UtEel2zM8Mw4lOQUuGGYWSOcNXIEOBbVf0oatseuES6A3GJn4ZhGIBFMAzDiIOIFInIH4GlwHzgAxF5WUQ2DLfDvheYhysHvBfYwUdzDcMIGJaDYRhGTETkXOA64L+4DpVbATsBD+PK14bjGlBdpqrW9dMwjE6Yg2EYRkxE5HXcMuoIDcsCi8gtwC+Ab4CfqOp8H000DCPA2BKJYRjx2AaYpZ17DtwWfp5izoVhGIkwB8MwjHhU4OSFvURev51jWwzDyDPMwTAMIxHRa6iR1825NsQwjPzCylQNw0jEoSJS5XndF+dkHCsiO0Xtq6p6fe5MMwwjyFiSp5EWIhICPsRpHmytqo2ebbfjOvWNUdX7fDLRyBAikq6uhapqUVaMMXo9NvfkHxbBMNJCVRtE5HfA7cDZwPUAInIVcDpwjv0HLxj299sAw4hgc0/+YREMI21EpAhYBGwIDAbOwP1n/52qTv7/9u4QJYMoigLwudswuABBxGSwC4KgRQQX4Wqsdrtg02YxaBTEDZhNpmv4DUb9eTjM8H1tmHLDzOHMe8zMlLMByyV75kXBYC1VdZTkJsl9Vk+6l919Me1UwNLJnvlQMFhbVT0l2U1yneS8XUzAP5A98+A1VdZSVWdJdr4PP9zgwH+QPfNhBYM/q6qDrJYob7L6HsJpku3ufpl0MGDRZM+8KBj8SVXtJblL8pjkMKufXr0kue3ukylnA5ZL9syPLRJ+raq2ktwmeU1y0t2f33/RvEpyXFX7kw4ILJLsmScrGPxKVW0meUjymWS/u99/nNtI8pbkubvd6MAwsme+FAwAYDhbJADAcAoGADCcggEADKdgAADDKRgAwHAKBgAwnIIBAAynYAAAwykYAMBwCgYAMNwXujMdB9tKzBMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nds = 100\n", "m = 2.0\n", "b = 3.0\n", "x = np.linspace(2, 3, num=nds).reshape(nds, 1)\n", "y = m*x + b + 0.1*np.random.standard_normal(x.shape)\n", "res = y - modelfun(x, [m, b])\n", "\n", "sns.set_context('talk')\n", "plt.figure(figsize=(8,4))\n", "plt.subplot(1,2,1)\n", "plt.plot(x, y, '.b');\n", "plt.plot(x, modelfun(x, [m, b]), '-r', linewidth=3)\n", "plt.xlabel('$x$')\n", "plt.ylabel('$y$')\n", "plt.subplot(1,2,2)\n", "plt.plot(x, res, '.g')\n", "mr = res.mean()\n", "plt.plot([x[0], x[-1]], [mr, mr], '-k', linewidth=3)\n", "plt.xlabel('$x$')\n", "plt.ylabel(str('Residual, ($\\\\mu$ = {:5.4e})'.format(mr)))\n", "plt.tight_layout(rect=[0, 0.03, 1, 0.95])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Initialize MCMC Object and Setup Simulation\n", "- In the simulation options, we turn off the waitbar and set the output verbosity to 0. The results will be saved to a json file.\n", "- We define a directory name based on the date/time to ensure we do not overwrite any results.\n", "- We add the model parameters, but their initial values will actually be determine later." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "mcset = MCMC()\n", "# Add data\n", "mcset.data.add_data_set(x, y)\n", "datestr = datetime.now().strftime('%Y%m%d_%H%M%S')\n", "savedir = 'resources' + os.sep + str('{}_{}'.format(datestr, 'parallel_chains'))\n", "mcset.simulation_options.define_simulation_options(\n", " nsimu=5.0e3, updatesigma=True, method='dram',\n", " savedir=savedir, savesize=1000, save_to_json=True,\n", " verbosity=0, waitbar=False, save_lightly=True, save_to_bin=True)\n", "mcset.model_settings.define_model_settings(sos_function=ssfun)\n", "mcset.parameters.add_model_parameter(name='m',\n", " theta0=2.,\n", " minimum=-10,\n", " maximum=200,\n", " sample=1)\n", "mcset.parameters.add_model_parameter(name='b',\n", " theta0=2.75,\n", " minimum=-10,\n", " maximum=100,\n", " sample=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Setup Parallel Simulation and Define Initial Values\n", "- You can specify the number of chains to be generated (`num_chain`) and the number of cores to use (`num_cores`).\n", "- Note, the initial values are defined in a 3x2 array. The expected size is [num_chain x num_par], which is 3x2 in this case. If you don't specify initial values, then a random set will be generated within the parameter space." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# setup parallel MCMC\n", "parMC = ParallelMCMC()\n", "initial_values = np.array([[2.5, 2.5], [1.8, 3.8], [2.05, 3.42]])\n", "parMC.setup_parallel_simulation(mcset=mcset,\n", " initial_values=initial_values,\n", " num_chain=3,\n", " num_cores=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Run Parallel Simulation and Display Chain Statistics" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Processing: resources/20190719_131740_parallel_chains/chain_0\n", "Processing: resources/20190719_131740_parallel_chains/chain_1\n", "Processing: resources/20190719_131740_parallel_chains/chain_2\n", "Parallel simulation run time: 1.7680649757385254 sec\n", "\n", "****************************************\n", "Displaying results for chain 0\n", "Files: resources/20190719_131740_parallel_chains/chain_0\n", "\n", "\n", "------------------------------\n", "name : mean std MC_err tau geweke\n", "m : 2.0206 0.0372 0.0017 13.7387 0.9941\n", "b : 2.9497 0.0903 0.0042 14.3948 0.9910\n", "------------------------------\n", "==============================\n", "Acceptance rate information\n", "---------------\n", "Results dictionary:\n", "Stage 1: 17.16%\n", "Stage 2: 52.70%\n", "Net : 69.86% -> 3493/5000\n", "---------------\n", "Chain provided:\n", "Net : 69.86% -> 3493/5000\n", "---------------\n", "Note, the net acceptance rate from the results dictionary\n", "may be different if you only provided a subset of the chain,\n", "e.g., removed the first part for burnin-in.\n", "------------------------------\n", "\n", "****************************************\n", "Displaying results for chain 1\n", "Files: resources/20190719_131740_parallel_chains/chain_1\n", "\n", "\n", "------------------------------\n", "name : mean std MC_err tau geweke\n", "m : 2.0180 0.0424 0.0032 22.3207 0.9806\n", "b : 2.9565 0.1074 0.0081 21.9919 0.9662\n", "------------------------------\n", "==============================\n", "Acceptance rate information\n", "---------------\n", "Results dictionary:\n", "Stage 1: 22.72%\n", "Stage 2: 54.94%\n", "Net : 77.66% -> 3883/5000\n", "---------------\n", "Chain provided:\n", "Net : 77.66% -> 3883/5000\n", "---------------\n", "Note, the net acceptance rate from the results dictionary\n", "may be different if you only provided a subset of the chain,\n", "e.g., removed the first part for burnin-in.\n", "------------------------------\n", "\n", "****************************************\n", "Displaying results for chain 2\n", "Files: resources/20190719_131740_parallel_chains/chain_2\n", "\n", "\n", "------------------------------\n", "name : mean std MC_err tau geweke\n", "m : 2.0158 0.0512 0.0049 40.5945 0.9675\n", "b : 2.9616 0.1299 0.0127 42.6359 0.9448\n", "------------------------------\n", "==============================\n", "Acceptance rate information\n", "---------------\n", "Results dictionary:\n", "Stage 1: 15.16%\n", "Stage 2: 53.60%\n", "Net : 68.76% -> 3438/5000\n", "---------------\n", "Chain provided:\n", "Net : 68.76% -> 3438/5000\n", "---------------\n", "Note, the net acceptance rate from the results dictionary\n", "may be different if you only provided a subset of the chain,\n", "e.g., removed the first part for burnin-in.\n", "------------------------------\n" ] } ], "source": [ "parMC.run_parallel_simulation()\n", "parMC.display_individual_chain_statistics()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Processing Chains\n", "For this example I assume that the MCMC simulations were run, then the user went to process the saved results at a later time. We load the results from the simulation and demonstrate how to plot the parameter distributions and pairwise correlation using the [mcmcplot](https://prmiles.wordpress.ncsu.edu/codes/python-packages/mcmcplot/) package.\n", "\n", "To load the results we must specify the name of the top level directory where the parallel results are stored. We previously defined the `savedir` variable, so it still exists in memory. If it wasn't, then we would need to define it as the string matching the directory name." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'resources/20190719_131740_parallel_chains'" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "savedir" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAH8CAYAAABcl+LaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXRUdZrw8e+vqlJVSWXfQ3aWJCwiO4iKEBBcQMGtV51e35l2FLG7bV/Hnp7p1XkduwW3XqbH7pbW1nZBRWlQCIgIQtiFQAjZQ/Z9rf2+f1RSSSCELUkl4fmck5N769a997nnEOqp3/L8lKZpCCGEEEKIkU/n6wCEEEIIIcTAkMROCCGEEGKUkMROCCGEEGKUkMROCCGEEGKUkMROCCGEEGKUMPg6gOFEKeXEk+w2+zoWIXwoGHBrmjbi/n+Qv2EhgBH8NyyunJJyJ92UUm5AhYSE+DoUIXymqakJQNM0bcS16MvfsBAj+29YXDnJ5ntrDgkJCWlsbPR1HEL4TGhoKE1NTSO1xUv+hsVVb4T/DYsrJNm8EEIIIcQoIYmdEEIIIcQoIYmdEEIIIcQoIYmdEEIIIcQoIYmdEEIIIcQoIYmdGD6ay6G93tdRCCGEECOWlDsRvtdeD1v+DY78DZQexmXC0l9AdEbv9zWWQt5HoPeDpPkQOd438QohhBDDlCR2wrc0Dd74OhR/1rnvgtMfQ9GnsOhJmP1tqC+ET5+BnPdAc3eeqOCW/4J5/+Kz0MUo57R6vnTkboKKLyDuGki/DQLCwWD2dXRCCNEnWXmiB6VUoxQ3HWKnt8Jf77788+94AWbcP3DxiK7ipk2apoX6OpZLNWB/w04rlOyF1+4Fp637dYMJvvomJM2V5E4MWyP5b1hcORljJ3xr56+7t6d9Df5lF0RPPv/7k66DoDHd+1v/E+xtgxaeuEq115+b1Cmd59/fF38HazO47L6LTwghzkMSO+E7NblQstuzrTPAwicg9hr450/Im/5vFGqxALRo/nzomsOdtp/xbOJz8FA2BCd4zmuvhX3/46MHEKNW7qbeSd3kVfDtjyD1JmxjplPlaOGNU2/z090/5Y2Tb1DVVoWt5/uFEMJHZIyd8J2CHd3babdAaCIAbmXgn/PmUmCbghEHTvS4O7+DHNmWxzXxISxZ8AP44FHPuXtehOseAr38cxYDpOKL7u3Jq2DSSvjz7djSb+PwjC/x4Af3YXd3t9g9nf00Ly15iWlR0zAZTD4IWAghPKTFTvhO4c7u7bELvZvbc6spqPV0r5pM/mx/LJMbJ0R6jz/+9lFqJ9wLgTGeF9qquydfCDEQ4q7x/FY6uO5fYcP/AZeDxhtW8+DuH/dK6gDsbjsPbn2QRpuMzxVC+JYkdsI33C7PzNcuqQu8m/+7q9C7/ZW5SSRHWHjuy9OJCfa0hNS12XlmayFMurP7/Jx3Bz1kcRVJv80zUSJ1ARR+6umWTV3Ajqr95yR1XexuOzvKdgxtnEIIcRZJ7IRvVB4Fa5NnOzAWItMAqGmxsTu/DgCdgn+anwJAmMXIf9011Xv6G/tLyY+6uft6JzZ6kkUhBkJAuGf2a0gCNHR+0QiM4WTbmX5Py63PHYLghBDi/CSxE75RvKd7O/VGUAqAXadrvC/PTA4jPtTfu78oI5rMjGjAU/7uiWx/NG93bA2U7hv8uMXVwWD2lDRZ/B8Q1Vkou7WKDEt8v6elh6cPQXBCCHF+ktgJ36g+3r0dP9O7ufNUrXd7wYSoc0778e0T8dN7ksB9Jc2UhF/ffbBn164QV8pgBnMwZCz3dMsW7mRhzCyMOmOfbzfqjCxMWDi0MQohxFkksRO+UZXTvR09CQC3W+PTvO4WuwVp5yZ2Y6MC+eb1qd79P5cndh/sORlDiIFgMHu6Zb/yBuj9CN31HC/N/8U5yZ1RZ+SlJS8RapJ6sEII35L6EGLoud1Qc7J7P8ZTkPhkZQu1rZ6B6WEBfkyJD+nz9Iczx/POwTJqW+1sap3Af3QtAFC6DxxW8JMVAcQAMgVB8nWw+jCm3H8wraGCTcvfZEfVPnIbTpEens7ChIWEmcMw6vtuzRNCiKEiiZ0Yeg2F4Gj3bFuiweIpZXKwpMH7luvGRaDXKe/+7vLdrD2wlrLWMsZYxrBs9t28uj2YKsIpJo5kKsBlg7Jsz5g9IQaSwQzBY2D2tzEBMcCXwsaCowObvZVGezMbij7iZFM+6WFpLEhcyMn6k+ws28m40HEsTlpMoF8gwaZgXz+JEGKUk65YMfSqe3TDxkzybh4t664BNi2xu0vrnbx3eHDrg5yoP0GLvYXchlzer/wVkWM+B+Az58Tu60k9OzFUbC3Yak5wuOoAt31wH7/Y/zRv5b3NL/c9xYoNK3C6nbTaW3k6+2lWbFjB0dqjNNuafR21EGKUk8RODL1e4+u614U9Utrk3Z6a4EnsipuL+fnnP8elnVvKxBbyLvqAfLLdPWYinjkw8PEK0ZeORhp1igc/e7LPgsVPfPoE90++H53SYXfbeSTrEVodrT4KVghxtZDETgy9Plrs2mxO8qpbAE/lk67xdU9nP43T7QRgfOh4/nzLn7k26lrv6UEJ73KIlO7rle331EIRYrCVH2JHZfY5SZ1O6ZgXN4+lKUspbylnTswcwJPsZZVk0W5twtZaBU6rL6IWQoxyktiJoVef370dMQGAY2eacHfmYxOiAwk0GdhXsY+dZZ6ZrgrFL274BTNjZvKbhb8hyC8IAJe+hvKwPBo1i+fkjvrugrJCDCZr4zkFi5clL+OVW19hbtxcjHoj1R3VPDb7MZYlLwMgvzGfLxpOcrj6ELaaXEnuhBADTiZPiKGlaVBf1L0f7ildcrTs3G7Y13Nf9762cvxKJkd4um2jA6JZM3MNP//85wCYI3aT3TKWm1Xnwu1lByB87CA+hBCAOZQMXXfB4mXJy1iaspRvbf5Wr1Y8o87IUzc+xYSwCUyJnMIYyxieOvYXMub9OyaHzTMxQwghBoi02Imh1V4Hdk+XK34WsHhq1Z2o6B5Ufk18CLUdtWwv2e597YFJD/S6zKoJq4i1xAKg6Vt5M6hHaZQz+wcpeCF6GDOdhbGzMeqM6JSO+yffzxOfPnFO1+yixEXEWmIJMgaxtXgreY15/Of8/+Rgw0leK/yA8tZy2uxtPnoIIcRoI4mdGFr1PbpJw1O9S4mdrGzxvpwRG8R7p9/DqXnG1k2Pns74sPG9LuOn8+Mbk7/h3T8YVoe7c7u9KHtQQheiF/9QQt0aL13/S+bHzSe7j/F2Xa1439z8TQ5WHWTlhJUUNBXw4pEXOdN6hpmxs1h3YB2Haw5LcieEGBDSFSuGVs/xb2EpADhdbk7XdM8WTI8N4r8+/sC7f0/aPX1e6q4Jd/HS4Zdotjdj82tjn9nEPKsNQ02OpwiyTr63iEFkCsIUMYFptjj+K242zx56vtfhrla8b23+FosSF/XbTfv+6fcZG5KKxWgZ6qcQQowy8sknhtbZLXZAUV07dqenvS0m2ESjo5zTjacBMOlNLEla0uel/A3+LB+73Lv/WlAYAEZ3B+56mUAhhoApCFNwHCH+EWSEZ/Q6NCd2DtmV2Tg153m7abvKonx10lfJKtnuWTlFCCGugCR2Ymj1arHzJHanqrq7YdNjg9lastW7f/2Y6wnwCzjv5e6acJd3+1OLkcbOVrrinM8HKmIhLspNCQt7rSEb6R9JaUupN8E7O6nrYnfbOVR1CLfmBpdjqMIVQoxSktiJodVHi13P8XXpMYFsLe5O7JYk991a531/eDpTIqYA4FTwscUfgKpTMs5ODK0QUygvLlrrTe5qO2pJCU7xJnj9KW4p5taxt2Jrrx2KUIUQo5gkdmJo9dVi1yOxi4uwcbzuOAAGZWBBwoILXnL5uO7u2M0WzxglVXVsIKIV4qIF+JmZFj6Zt1a8xY9m/4iEwAQWJy2m0dpIYlDiec/TKR1zYuawMX8j+e52bLZWcPXduieEEBciiZ0YOk4btFZ5tpUOQhIAyO3RFdvcVYsOmBk7kxBTCBdyc/LNKDyza/ebTdTqdSTa8ylv7BjA4IW4MLMpkP3lnzN/zHwcbgcHqw7ylYyvMCd2Tq9u2i7Lkpex/tb1NNgaKGouYnfFHuoczRyuPUZVWyU2p80HTyGEGMkksRNDp6WiezswBvR+WB0uiuo8ZR6UgtOt3V2oC+Iv3FoHnoLFs2JnAeBWii2WAMaoerJzTg9c7EJcDL2RBcmL+MORP7AgYQG/3PtLkoOTCTAE8JuFv+mV3PUshfLUvqd4J+8d1h1cx4oNK2hxtNFibyW/MV+SOyHEJfFpYqeUmq2UelEplaOUalNKlSilXldKjb+Ic+copV5SSh1QStmVUrJA6HDXXN69HTwGgLyqVu/SrskRRvZV7vW+5WK6YbvcknKLd3tbgGeyRUWujLMTQy/UFMZdaXexrXgb629dj9Vlpbq9mpN1J3nrjrd4fPbj3DPhHh6c9uB5Z8qu2b4GN25sLhtN9qbz3EkIIc7l6xa7x4G7gK3AI8AfgIXAIaXUxAucexvw3c7t/P7eKIaJPhK7nt2wsdEVdDg93aeJQYkkBydf9KUzkzK92wfNJpp0OpzlR68wYCEunclgYlrUNL4/9Z/5ovooQcYgvqj5gptTbuZk3UmiAqKID4rnk7JP+p0pu6d8DxVtFTjdTtrsLX2+TwghzubrxO43QLKmaas1Tfujpmm/AG4E/PAkff35LRCsadpMYMsgxykGQnOPBdODPWts5lZ2LyWmBeR4t2+IvwHVuSrFxYj0j2Rq5FQAXErxqb+ZOOtpKppknJ0YeiaDiRi9mZVO+KT0E146+hItjhbez3+f1OBUipuKKWou6vcaRU1F7DqzixUbVnC45iiN1vqhCV4IMaL5NLHTNG23pmn2s17LA44D/bbYaZpWpWmafGqPJH222HWtOKFR6TjoPXxTwk2XfPmFiQu92zsC/JmkitlbIB+GwkcCwvGLuYZTDbm4NTfrj69n1YRV/O8X/8ud4+/sd6YswNjQsTw0/SF+MOsH/DXnr7Q7rbS21YBTihgLIc7P1y1251CeZpoYYMALOimlGvv7AS48BVNcvn5a7HTGGurtnskVAYYAZsfOvuTL90zsdgX4k6TOcLCg6vLjFcPOiPobNpghYpx3RYotxVv4qOgjvjrpqzRYG7g19dY+Z8qCZ6mxRYmL0Cs9y8cu52fX/wyby4bbYEKrzpXkTghxXsMusQO+BsQDf/d1IGKA9Wqxi6ex3U5Vs2fGnynkpPfQ/DHzMer7/sDrz/jQ8YyxeFoC23Q6cvwNNJZ8cYGzhBhEeiMLExd5E7gtxVt44B8P8EbuG2wv2c6vF/76nOTOqDOy6a5N6JSOj4o/4jcHfsPmos2Y9CasTiu1IXFY2xskuRNC9Mng6wB6UkplAC8Cu4D1A319TdNCL3D/4fWNf7Q5qyu254oTltBTdBV1uJTZsD0ppZgfP5+3Tr0FwB5/M341J7E6XJj99JcbtRhGRuLfcKgplBcXv8i/bvtX7G47bs3N3sq9HKo+xO9v/j1vrniTPRV7yG/MJzUklVtSbiG3IZdHsh7pNbli7YG1rMtcR0ZYBgX2asbaW9AHJeBnOv+Se0KIq8+wabFTSsUCHwINwL2aprl9HJIYSC4HtFR27wfFcbKic+KEvg2bwTOxWaEuO7EDz9qyXT7zNzOWMnIqmvs5Q4jBZTKYmB49nQ9XvseTc5/knrR7WDNjDS/f8jKvn3yd3eW7mRw+GbvLzu7S3TjcjnOSOvDMlO16Pdw/kib/EGzOVmxOl4+eTAgxHA2LxE4pFQL8A8837WWaplVe4BQx0rRWAZ0F6yzRYDByosLTYmewnPIeuybqGiL8Iy77NnPj5qLvXIXihNFInL6UI6WNVxK5EFfMZDAR69L4cuAErom4hj3le3jgHw/wccnHzIqZxbc/+jbv5b/H96Z/j20l2/otg5JVkoVJZ6KmowarUjjd0iUrhOjm88ROKWUGNgJpwHJN03J9HJIYDH3MiD3ZOXHCEHjCe2hhwsIruk2QMYipoRMA0JSi2VLF0TIp8CqGAf9QUApHRx17K/fi1tzMiZ3DrvJd3kQuxhJDfmP/ZTnzG/OxuW24NTdWl40PCj6gsKmQ6vZqWaVCCOHzlSf0wBvAdXi6Xz8/z/uSOsffiZHqrBmxLrfWWZzYhSHwlPfQlXTDdrkusbtUSpG/lZyS6iu+phBXzBQEkRNYOO5274SJSP9ISltKvW+paqtiXOi4fi8zLnQcla2VZFdmk1ufS2ZSJh8VfsTh6sOyBJkQwuctdr8G7sDTDRuulPp6j5+VPd73CnCi54lKqWSl1I+VUj8G5nS+9uPOnxVD9QDiIp3VYldc14bV4UZvyUfpPV1JcZY40sLSrvhWs8Zc590+YDahbzhNm815xdcV4oqZggg1R/DSoucw6ozUdtT2qmf33a3fZXHS4n7LoGQmZfLP2/6Z0pZSmmxNdDg7yIjIYGzIWIx6I422hqF6GiHEMOTrxG5a5+8VeGbB9vxZe4FzU4Gfd/50fZJ37d894JGKK3OeGbF+wYe9Ly9OWnxJq02cz9SoqRg7x9kVGf1I1OdzQiZQiGHCZDAxLXQ8m5b8kZsTFnFTwk3eRK7d0U6rrZV1mev6LIOyLnMdrbZW2h3tJAYlYjaYOVF3gjdy36Dd2U5VexUVbZV02Nt88WhCiGHA1ytPLNQ0TZ3nJ+Xs95117o5+zv3GUD+LuICzumJPVjSDsmMIOuZ9efm45QNyK5PexDXG7gkYZsspjpdLYieGD5M5hJiOZu479A5JfsGsW9SdyN39wd1Em6PZuGojj89+nHsm3MPjsx9n46qNRJujufuDuzHqjMyOnU18UDx7K/fy7/P+naq2KtbnrCc6IJoGWyNtktwJcVUaVnXsxCh2VovdiSMtGIJOoPSeQeMpwSlMCp80YLebHZbBgapdAHQEVHG8XCZQiGHEYIakuRA5AWP5YdLGTOVPt/yJnLocTjWcYn/1fubEzmHFuBV0ODuobK1k1furaHe0Y9QZeerGp7C77Lx+8nUmRkzkp3t+ylcnfpXr46/n4+KPmREzA4PeD4tO77mXEOKq4euuWHG1ODuxq2jGLyTb+9LyscsHpBu2y6yEG7zbpf7t0mInhh+DGYLHoNKWEeJnoc3RxpKkJdw9/m5uTbmVdmc7+yr24dbcnKg/wW0pt/HDWT/kzRVvkhqSyusnX2d76XZmxc7yLFumwV3j78KgDLyV+xZOtxNXYxmNzc1S606Iq4gkdmLwuV3QUuHdbTZGcaatBEPgaQB0Sscd4+4Y0FtOHXsLBs1TG6/MqCisLcXulJrXYngymYKZETWNzUWb+fKmL/N+/vu0Odp4O+9tskqySA9PJzEoEYfbwVN7n+KejfewvXQ7r972KjH+Mdw5/k4Kmwv57/3/jU7p+O7U71LZVglhSQQZXOwvapDkToirhHTFisHXVgPuzlmp/uGcqnNiDNvjPXxTwk3EBcYN6C39AyKY6lQc9PPsR5gPkld9O5PHDKvVpoTwMvn5c7rR82XnmQPPcHvq7fzHdf/BJ2WfkNeQx5LkJWwv3U58UDyrx6zm5uSbabI1cbT2KE98+kSvosbG/Uaey3yO6o46TtXnkhaTRlNHB9FBgb56PCHEEJEWOzH4enXDxnOkvBq/kAPel76c8eVBue3MHhMoAgJkAoUY/saHjvduf1j4ITm1OSxIWMCUyClUtFUQ5R+FzWVjb/ledEqHhnZOUgeeFSpWZ62mtqMWm9vO2oPrONX0hdS4E+IqIImdGHxnja/LKt2M0ns+YEINY5gXN29QbjsrLN273RFQSY4kdmKYy0zK7FXmZM0nazhRe4Iwcxj+Bn8SghL4qOgjMiIyOFF/gn2V+/pdfuxQ9SHeyXuHr076Ko/uWE29VbpkhRjtJLETg69HYqcFxXGibbN3f0nCKnRqcP4ZThszzzvOrsHUwZHyMxc4QwjfCjGG8Fzmc+ckdy8eehGDMmB32XnqxqdICkqixd7Sa9WKvhQ2FxJmDmN/5X6ujbqWnWU7MNiacNjaB/tRhBA+ImPsxODrUcNuv9GAQ+/Z19x+fHvafYN224CYqUzab+eo2QTAqcZjuN03o9MN3OxbIQaSxWhhWtQ0Plj1AdtKtpHfmM+40HEsSV5CfUc9CUEJxFpi8Tf4c7D6YK9VK/qSGpzKzrKdmPQmogKiyG04hT56Dro2B4SM8SxzJoQYVaTFTgy+Hi12r/VoYfC3zSYhJHzw7hs5gem27jFFLr98iuqkaKsY3ixGC3GBcaSZb2HN9MdYMW4Fv9jzC0paSqjtqMVisBBgCGByxGTmxM7pd/mxadHTyK7KJiM8gyj/KJYmLwU/C8o/BBxWcPXdjSuEGLkksRODrzOxa1eKne2nvS9fE7J0cO9rDmaaCvDu+gfkywQKMWLEhwXwP5+UsDF/I5+Wf8rYkLEEGAI4UnMEt+YmryEPheKpG5/qc/mxZxc9yyvHX8GgDCxIWMC96fdS3FzM38qyKHR38EbZNt7Ie4eqtiqZVCHEKCKJnRh8nV2xWywB2HEA4LJFc2Pi9EG/9bXBKd5tl7mCL87UD/o9hRgIIQF+fGVeMqcbT+PW3Gwr2UarvZVxYeP4/MznlLeWc6DqAFMjp/LWHW/x2KzHuCftHp6c+yRvrniT90+/T4ezgw13biDYGIy/zp8b4m8gPSyd+zbeR5g5nOzK/dz2zm0cqj4kyZ0Qo4QkdmJwaZq3xW5joMX7sqNxFtOSwgb99lGRk4h3eGroaTo3ByqOXeAMIYaHQJMfkRajtwTK747+jhBTCJqmEWAK4MaEG3n+0PPcuuFW3jn1DjNjZvLDmT8kxj+G7MpsHp7xMDfE38DLx17m7by3aXG2oKHRbG/mnTvfIaski/sn349Tc/Kv2/6VBlujj59YCDEQJLETg6u9Hlw2GnU6DnROYgBQrdOZGBc8+PePTGdaj3F2+S3H0Tpnygox3Jn89N4SKG7NzabCTRQ2FRJniUOndDy/+Hl06Phzzp/5xpZv0GRv4lTDKeIC47jrvbv41b5f8Xbe2zyz/xnuef8ecupyCDQG8sLBF7hj3B04XU5mx8zG7razo3QHdpeUQhFipJPETgyuzm7YnQH+uDvXgnV1JHLtmCTMfvrBv3/kBKZZuxM7q76AqmbpchIjR88SKL87+jvGhY5jd/luXjz8omfM3O1/48dzf8yrt71KVkkW18Vfx6PbH/XWt9MpHfPi5rE0ZSkb8jZg1Bv5l2n/wv7K/UQFRJGZmAlAbn0udocsuyfESCeJnRhcnd2wOwL8vS85WyYxK2UQZ8P2FJnWq8VO71/M8fKmobm3EAOgZwmUx2Y9RllLGacaTlHdXk2bo417P7iXBmsDSUFJuDU3B6oOeJO6ZcnLeOXWV5gbNxej3sjM2JkE+gVS11HHxMiJ6HV6lqUu4y/L/kJGWDpFdVLfToiRTurYicHVfAabgl3+Zu9LzpZJzBmqxC54DOMxEeB2067TofNrYm9JAYsnxgzN/YUYABajBYvRwtfHraTJ7aSwqZANpzfw0PSHMOqMXBd/HdmV2YwNHcvW4q2AJ6lbmrKUb23+Vq/VKX57+LesXbSW6rZqcmpzuC/9PlJDUokLjMNmk/InQox00mInBldzOUdMJjp0nn9qbnsEmiOaGcmDP3ECAKUwRE7gmh4fWAeqDg7NvYUYaHoDVkcz8+LmYVAG1h9fz2+X/Jbsymy+/8n3yQjPIDk4GZ3Scf/k+8+7juya7WvITMokwj+C3x/9PR8WfoiGhsHPRltbDTitPnpAIcSVksRODK7mcvabe7TWtY0nPSaYEH+/oYshMq3XOLui1pyhu7cQA8lgJtQcToeznadufIrtpds51XCKspYyrE4rJU0lLE5azPy4+WRXZve7juymwk3sKN3B23lv8/+y/x8rNqygqKWIdgVUn5TkTogRShI7Mbiaz5DdYzasqz2V68dHDm0MkRO4tsc4uw5dAU3tjqGNQYgBYjIGMi50PNOip/HWHW8RGxBLakgqAP+05Z+wOqx8f9b3L7iObGlLKVEBUd59u9vOI1mP4HA7ICwVzdY6qM8hhBgcktiJQWVrPsNRU8/EbizXj48Y2iAi05naI7HTmcs5WFo5tDEIMYBMBhPRAdGkhqSyIGEhy1KWeVefuPuDu9lVtouM8Ix+r5EYlEhNe02v1+xuO1klWVj9TDg0N1pLpbTcCTHCSGInBo+mcdRag13nKXOCLRy9O4Q5qUOd2KUR4tYYZ/d0SynlZkfRoaGNQYhBYtDrCDIE8MLiF7zJXZAxiHlx8/pdR3ZW7Cyyq7LPOZbfmI/L7aIBjSpcWJsraWxuxuaUGndCjASS2InBY23igF93MWB7+zimJYYSaBriydjhqaD0TOsxgeJIzeGhjUGIQWQxBTM+ZBzvrXyPJ+c+yYyYGeyr2Oetf9eTUWfkqRufYv3x9bi1c+vWjQ0dyx+++AMFzQXYXE4OWSvQudvZX9QgyZ0QI4AkdmLwNJfzRc9u2I6UoR9fB2AwQVgK1/aYQFHWcWLo4xBiEGUVbWHdgXUsSVpCclAy4f7htDnaeG/lezw26zHuTbuXJ+c+ycZVG9lWvI0txVvOuYZRZ2RR4iJeyXmFh7Y9hNFg5N92/RvtBo2XPyugoU3KoQgx3EkdOzFotKYzHDN1txa4OxJZlBHtm2Ci0plWUOzdteoLaGy3ExrQd1eVECPNyaZ8NhdvxmgwckvKLcyKmcXag2vJqc1hatRUDDoD40PHU9xczIrxK2i2NxNqDqW2o5Z9lfswKAPrFq1ja/FWnG7P+sqflH3CmyvepKyljOfvm8g7R6r5+rxkHz+pEKI/0mInBk1F3Qnq9Z5lw0wuRajfGKbGh/gmmMgJpDichHSuhan07WzO/cI3sQgxCLomS1S3V9NibwEgLSyNv5z4CxkRGfx6/695/eTrpASnkByczIKEBZj0JhYkLGDjyo1suXsL+Y35/PrAr73XPNVwig9Of0CrvZVWrYObJ/j3eW8hxPAhiZ0YNMfqu7s7I2yBLEyPQdc1kWKoRaahoNc4u+1F+3wTixCDYGHiIow6Iwitva0AACAASURBVPsq9xEfFE9dRx2LEhehQ8eR6iM8n/k8eqUnvymfle+u5Ff7fsXbeW/zdPbTrHpvFcfrjjM/fj7Lkpd5rzk2ZCzLxy+nw9GBTumwGjqoa2uTsXZCDGOS2IlBc6ylyLsdYA0n01fdsACR6QBc06PsSU79cV9FI8SACzWF8tKSl7wrUlhdVpptzbyw+AV2l+/21Kmb+QiPZD3S52oU39/xfTqcHSxNWcqyZE/5lBvjb6S2o5ZqazUvHHqBXWd2YdOaKGuUiRRCDFcyxk4MmqPWWuhsoFPWeG6cENX/CYMpcjxAr6XF6p2nsTldmAx6X0UlxIAxGUxMi5rGprs2caL+BFVtVSQFJ5EYlMjqGas5UnOETYWb+l2NIrsym+zKbB6c9iAPTHqA4pZiHt3+aK9zfnPgN6zLXEeIbQomQ+hQPZ4Q4iJJi50YFG7NzUmtu7BpfPC1Q7uM2Nn8w8ASzeQeLXbKVMHnBVW+i0mIAWYymIixxLAwcSG7zuzino338NPdP6XOWses2FkXtRpFmDmM47XHCTOHnZPUQfcKFXatnQ5by2A+jhDiMkhiJwbFmdYztHX+6wp1uZg5YbZvAwJvoeIUu2c5MaXcbDx5boFWIUaDjPAM3Jqbzys/Z1PBJooai0gLSzvv+3VKx9zYuYwLGYfD7SC/Kb/f1r2skiycgMPWPkhPIIS4HJLYiUFxvOKgd3uC3cn1U9N9GE2nKM+HWs/u2M/LpVCxGJ0WJi70Fidef2I9kQGR3BB/Q5+rUSxLXsb6W9fTYGugtLUUpRRJQUm9JlKcraCpAKUUerftvO8RQgw9GWMnBsVnp/d4txPsBsbHBPkwmk6RXYmdjY1BFgBq7XlUNVuJCTb7MjIhBlzXZIoHtz6I3W1na/FWZsbM5IXFL/DQtoe8rXHLkpexNGUp39z8zV4tdF0rVAB9FjOOD4ynsLmQ+MB4aG3E7jIRZjHKmFUhfExa7MSgOF7bXeokSYWilI/KnPQUOQHo3WKn9y/lk9ya850hxIjVczLFj+f9mJKWEjTNxYSwCXyw6gOemPME90y4hwenPcgTnz7R51i6Jz59gvsn349O9f6oMOqMzI6dzbc2fwur00qF9QyBhhbq62W2rBC+Ji12YlBUOCq8/7omBif5NpgunSVP0u12/DQNh1LojPW8f+wU981O9HFwQgw8k8FEjCGGL6V/CYC29nb0jiaa9IrbU2+HsbAxf2O/Y+kOVR1iTuwcPq/4HPAkdesy1xETEMPChIVklWRx1/i7aHe2E2HswNFhwBQks2WF8BVJ7MSAO15RSauhAwA/TWNGwiQfR9QpOB78AvBztDPRZueo2bOO7b7yIzS23yjLi4lRz2A00WD1J9BezRf2QmIDYjndeLrfc4pbivnl9b9ke+l2NDTmxM6hsKmQUw2n+OHsH5Jbn4tRb8SgN4AxGM1tw+6yY9TL35MQviBdsWLAvZvTPdN0nN2Bf3iqD6PpQaeDiHPr2WEq4aPjUvZEjH4mg56w4CCcxlimh0wg0BjoXYrsfBICE/jJ7p+QGJRImDGMFkcLhc2FbC3ZysaCjaSEpNBga+DvuX/nV/t+xcb8DzxFjdursTllYoUQQ00SOzHg9pR1r8GabrdDSIIPozlL5wSKKT3q2en9S3nvyBlfRSTEkDIZ9IQGB2O2ROJyO3vNnj2bUWdkVuws9lTs4eGsh8mIzOA7W77DuoPreD//fZptzeyr2Eebo42j1Ud5O+9tfrXvV6zYsIK8hjxON+bTbJVyKEIMJUnsxIByuNwUNXd37aTbHRAyjMavRXnG2U3tOYHCXMbu/FrONHb4KiohfCLMHI4OHS8ufvGc5K5rVuz64+txa27vzNpro67lR7N+xIerPiQ6IJpjdcf49MynrJm1hldve5VvTv4mbtyszlqNW3PR6mym1ebw0RMKcfWRMXZiQB0sbsBtPENXwYMMux1C4n0aUy+dM2MTnU5C0NGEG2VoB0M9Gw6W8VDmBB8HKMTQMRlMBGrBjA0dy4d3fciWoi0UNhWSGJTIrNhZrD++vlepk9KWUn42/2cUNBewYsOKXpMu1h5Yy7rMddwx/g7un3Q/mwo2caz2GJH+sUyPuA6XWyPEX8bdCTHYpMVODKgdpyrQG7vHq6Xpg8DP34cRnaWzK1YBkx1u78t6/zLeOlCGpmk+CkwI3wjwM6NQ+On8SAxMxOaysad8Dw/844Fz6td9Y9I3QOFZUuw8S435G/zZVLCJRUmLmBwxmekx13C0tJ7DJY00dfQ9+1YIMXAksRMDKiv/OOg8dazinE5CgofR+DqAiAmg86xZO7mtyfuy3lxKUV07+4sbfBWZED4TbAwGTWNixEQ+KvqIvZV7cWvuc94T5h/GtpJt/ZZH2V6ynQC/AO5+/27anG0UNxVTrT4hLtxGh0MSOyEGmyR2YsDUtdooaD7l3U+3DbOJEwAGI0R7ZgFO6THOTudfBsBb+8t8EpYQvmQymLCcOYhO6Xgu87k+x9t9sPIDKtsqyW/M7/da+Y35TImcwh+X/ZF3894lLjCOp7L/i69sXsmpxqMyU1aIQSaJnRgwu07XojNVePc9EyeGWWIHEDsV6J3Y6c3lgJsPjpbTbnf6KDAhfMccM4XAY+8QYgrm1dtf9a5MsWbGGv52+98w+5kpaS5hXOi4fq8zLnQcHxZ8yLc2f4vFyYvJa8jjF/N/gVNz8ugnD9NoaxyiJxLi6nTZiZ3ySFVKTe38PQzWjBK+tPNULTpzd2KXMdxKnXSJvQaAaJeLaOVpmVA6OzpTNW12F9tOVPsyOiF8IyAcS9QkUk9/SrBfELel3saUyCnsKd/DM/ufobq9mur2ahYlLuq3PEpmUiavnnzVuyTZmMAxLElZwjt3vENmYibbS7cP8YMJcXW55MROKWVUSq0FGoB84HDn7wal1LNKKdMAxyhGiN35NSOkxe4a7+ZkZ/c4Ip3Z0w275XjlkIckhM8ZzJA0F8uklYw5+REmnQGry8reyr3srdxLo62RRUmLyCrOYt2idX121z6X+RzbirfhdHtave1uO3sr9mJ1WDHpTDw661HGWMbQbGumzd7mi6cUYtS7nBa7F4E5wH1ANGDs/H1f5+vPX8xFlFKzlVIvKqVylFJtSqkSpdTrSqnxF3l+vFLq70qpRqVUs1LqXaXUMFni4OpT1tBOZVs1OoPnP2uL20280zm8ath1iZns3ZzS3D1ZQt85zm77yWqsDlnIXFyFDGYIHgMzvoZ/azWLEzMx6oy4NTfrj69Hc2uMDRtLQWMB7618jx/N/hH3TLiHH83+ERtXbSTBksDT+5/udcmCpgIcmgOd0nGs9hjh/uHY3XYO1xyW5E6IQXA5id3dwApN0z7SNK1W0zRn5++PgJXAvRd5nceBu4CtwCPAH4CFwCGl1MT+TlRKBQLbgRuBXwL/AcwAdiilwi7jmcQV2ldY36u1Ls1u9/zjGo4tdv5hEJoEwBRbd1Fi/8ByANrsLnbn1/okNCGGBb0RgscQ6hfEs4uexagzsqV4C88efJa0sDSWpCxhZ+lOYgJi+ErGV7hj3B1Utlay/L3l51wqLSwNf4M/ep0ePXoCDAHYXDYq2ipotjf74OGEGN0uJ7HTAL/zHDN0Hr8YvwGSNU1brWnaHzVN+wWeRM0PT9LXnweB8cBtmqb9t6ZpzwJLgXjg0Yu8vxhA2UX16M1ndcPq/MAS7cOo+hE/C4DJPSZQaH5nQHm6kLJOyjg7cZUzmDGbAnG5Xbx2+2s8OfdJgkxBfFT4ESHGEFZNWMW40HHEWeK47Z3b+Kct/9TrdJ3SccOYG1iSvIRWeys5dTlMjJxIcXMxLreLmdEz2VayDRxWHz2gEKPT5SR2fwU2K6XuVUqlKaWilVITlFL3Ah8Cr1zMRTRN261pmv2s1/KA40C/LXbAPcDnmqYd6nHuSWAbni5hMcT2Ftb3njhh61xxQjdMJ14nzAYgxO0msXNYqBuXt9VxT36dz0ITYjhJD0/n65u+zo6yHfxg5g+4O+1uWuwttDhaqGqrYkPeBp5e8HSvMXfLkpex/tb13JBwAy8ceoEPCz9kbOhYthZtRa/T87vDv6OqvYp4Szx2pxVsLT58QiFGl8tZUuz7wJPAM0AinhY6BZQCfwR+dbnBdM6sjQGO9PMeHTAVT9ft2fYBNyulAjRNO2flaaXUhebZh1xCuKJTXauNgpo2LGPLva+l2x0QOQzH13VJnOPdnGKzU9r5mWQMOIPVmkh+TRvVLVaig8w+ClD0Rf6Gh16IMYS1i9ayOms1b516izvG34FTcxKgAkgKTiI5JJnnDjzHn275EwerDuJv8CfOEsc3N3+zVyHj3x7+Lc9lPsfG/I2snrGabUXbWJyymAa3DYtbEWgK8uFTCjF6XHJziqZpLk3TfqZpWjIQBiQDYZqmJWua9nNN065k1PnX8HSn/r2f94QDJqCij2MVeJLMuCuIQVyiI2WNoGwoo6eVS6dpjHcM0xmxXWKv8YwjAia31Htfjozo7oL9vKD+nNOEuNpYjBamRU1j46oPMOgMmHQm6jrq6HB1YDaY2Vm6k01Fm7j/H/fz2ZnPmB49nUd3PNrnkmOrs1bznanfIa8hj+sTrqe2vRajzkhxRw3WlipoLgendM0KcSUuqsVOKRWDp4szBWgFDgIfa5rWBDT1c+pFU0pl4JlxuwtY389buxYe7at8ufWs9/SiaVroBWJoRL7xX7LDpU3ozeUo5RleOdbhwKxpEBzv48j6YTBB3LVQlt2rULEyl3q3Py+o445rx/giOnEe8jfsGxajBYvRwr0TvoJB56a6o5pDNYdIDEwktyEXwLMEmYKdZ3aed8kxFGiaRlpYGi2OFiIDInHhItI/klq3Az97K6FtNZii0j0zdIUQl+yCiZ1S6kZgExCApzWsS51S6ueapj13pUEopWLxjM9rAO7VtLMWKeytaxpjX/XyzGe9RwyBw6WN3iW5oMeEhLBkH0V0kRLnQlk2Eztn8LqBJucZUDbQTOwrlBY7IXoy+enB6WBWzCx2ndmF0+1kfGh3hapI/0hKW0r7PHftTWuZGDmRbSXbyG/MZ1zoOBYnLabD0UGTvQmz3syHZZ8wLWw806yxmAIlsRPiclxMV+wznb+/BSThabX7Mp6ixGuVUq9eSQBKqRDgH3i+ZS/TNO1C1WHr8bTW9dXdGodnzF9f3bRiEGiaxpHSRvTm7sTO2wIWluKboC7W2EUABGgaY12e7ywabgz+ZwDIr2ml1SbLiwnRi8GMn86PtLA09lXsY/6Y+d6JE7UdtSQGnTu2du1NazH5mVixYQVPZz/N23lv83T206zYsILilmKizFH46f1YlbaKf9v3Kxpd0h0rxOW6mMRuCvCspml/0TStTNO0Ek3T/q5p2jzgu8CXlVIPXc7NlVJmYCOQBizXNC33Qud0tuZ9Aczq4/BcIK+viRNicBTVtdPU4fAW94WeLXYpvgnqYiXP946zm9LePSsvJqoGAE2DnHKpsyXE2QKNgWQmZnKo+hAuzcVTNz6FUWdkX+U+ZsfO7jVD1qg3MjFyIo9kPdLnuLtHsh7BhYuc2hxcbhdvLH8DZfDH7jpPd64Qol8Xk9i1ACV9HdA07X+B14F/udQbK6X0wBvAdXi6Xz8/z/uSOsff9fQWME8pNb3H+9KBTODNS41FXL6jZY2ga0fXOXHCoGmk2+2gMwzvMXYAxgBPdyz0GmcXENTd4PvFmQEZQirEqBOu/Hhh0Tr+ePSPhJvD+fuKv/N/Z/9fTtad9BY1Bnhi9hNsK9l23nF3dredrJIs5sXN43TjaRSKVmcrzbZmbM6+hlILIfpzMYndduC2Cxwfdxn3/jVwB55u2HCl1Nd7/Kzs8b5XgBNnnfsSUABsUkr9UCm1BvgYTxfss5cRi7hMx8ub0ZvPePcn2B0YwbOUmE7vs7gu2jhPd2zPxK5DV+Td/qLsQtU1hLg6mUyBXBs0lkdmPkKDtQG35ubmlJvRNI0pEVN46463WDNjDTNiZpDfmN/vtfIb8wkyBhFniSOnLofPznyG3W2nxd4ss2SFuEQXk9j9DzBfKfXIeY6nAOXnOdafaZ2/V+CZBdvzZ21/J2qa1oJn+bFdwL8DPwcOAzdpmiaVZYfQ8fKmXt2wU2yd37CHezdsl7RbPL/sdvw0z6zeJkcl6D1rWEqLnRDnYTBjNlkItrYSZ4ml0dpIq72VhOAE1h5cy1+O/4U95Xuo7ahlXGj/3/3HhY7jUM0h9lbuJSk4iSPVR1ixYQW5Dadoba2S5E6IS3Ax5U62Ak7gN0qpVXiKEB/ofO0mPOu8XmgJsHNomrbwSt6naVoZF78urRgEmqaRU96MLmwEjq/rEj0JIibgV5dHhs3OF2bPZGu/gDIcLekU1LbRanMSaLqcWt5CjHKmIAL1fqTammnWHGgoChoKeHj6w2wu2szbeW9jUAZ+Mv8nrD2wts/uWKPOSGZSJiveXYHdZceoM/LMTZ45e6uzVrNx1Ub8rI2YAmOH+umEGJEupsXuZ3jKnZQCC/B0jR4DTgK/B3KBRqXUFKWUfPpdRSqarDS09544McU+whI7pWCyp+f/Wlv3eJ6ICE/3sqZBXpUsdyTEeRnMWCzRxAXGE2gIYGzYWO549w4yEzMx6oyE+YdR3V7Nusx1vSZVgCepW5e5jhO1J7yTJexuOz/85Ic8MPkBnJqTbSXbqHS0YW2vk5Y7IS7CBRM7TdP+U9O0lZqmpeBZ9WEx8EM8a8Yew7O811/xLAPWqpQ6opTqr8CwGCVyyptR+hZ0fp7uShOKsXaH5+BISewAJt0JwHRrd2JnCCj2budVtQ55SEKMREeLWpkaksbbd7xNg62BFxa/QKO1kSM1R7A5bGxctZHHZz/OPRPu4fHZj7Nx1UZsDhtrPlnT6zp2t50jNUeYHTObgsYCilrLaNIcaFUnJbkT4gIuqYVN07RGPJMltne9ppQy4imJMg2Y3vmzYgBjFMPU8fJmdP7dEyfSXeDXtTPcixP3FDMFoicxve6k96VWrRDPaAMDp6TFToiLMj4ugk2HC/lSXBVRURlUGIP5yXU/weF2sOq9VaDg8VmP87VJX6PB2uDtfu1LYVMhD09/mGO1xxgbMpZ2ZwfNkakYHR3o8MNkGAGTs4TwgUteK/ZsmqbZNU07qGnay5qmPaxp2g0XWvZHjA4nK5t7FyZu79GyFX45E6V9RCmY/nWiXG4SHJ4WRxd2dGbPnKC8ammxE+JihFmMJMZG80ZFDHqnnVBTKDaXDaPeyAuLXwANfr7357x/+n1O1p/st1ZdakgqQaYglqYsxag3crzmODl1OeS2FFNnrabNLi13QvTlihM7cfU6VdXSK7GbbO38j9YSDeZgH0V1maZ+CXQGplu7P2j0AUWAjLET4mKZDHpmpYSx8JpUdlfr0bnsGPVG3JqbjPAM3r3zXdbMWENKSArz4uadM+aui1FnZF7cPP70xZ9wuB3YXDbmx89ncsRk4ixxPHdgHUdrDktyJ0QfJLETl8XmdFFU14rOv7t2tbfUScQIaq3rYomEtFuYbuv+oPCzFAJQ3mSlxerwVWRCjCgmg57YEH9uyEggMCASs8FMnbWO7MpsTAYTy8cuZ0LoBEpbSr0rVvRk1Bl56sanKGkqISUkheq2aqL8o/iw4EM25m/Epbn4wewfsD7nFZpsDdicLh89qRDDk8xiFZelsLYNzVCDzuBZvS1UbybV0bmu6kjqhu1p+v3MPb3Zu2uwFAAuQM/p6lamJ4X5LDQhRqpwczg6DaL9o3G73dRYa7D4WQgxhfDioRd5+ZaXOVB5AIPOQHJwMglBCfzPkf/hoekPkRSSxOcVn/Nu/rukh6UzI2YG6w6sY8X4Fayevhqb24rN3obJMMJ6CIQYRJLYicuSV9Xq7aoEmKYPQXXtRIz1RUhXbvwSEs0RxDmdVBgMaMqGzr8Md0cyeZLYCXHZQv0CqHe0osNAclAyrY5WKtoqWDVhFX878Te+Mfkb7K3cS1ZpFolBiayeuZrTDad5dMejvWrfddW423h6I4/MfITK1koq2yqZapiKxWjx4RMKMXxIV6y4LHlVLej9i7z7M5zu7oMjtcVOb0BNvY95Hd3dsQZLHgBFtW2+ikqIkc9gxqI30+Jo4Y0Tb6ChsbNsJ+nh6dw5/k6+tulrPLP/Gd7Je4fnDz1PbUftOUkddNe4+9qkr7G9ZDuxgbGYdCbqrfWyrqwQnSSxE5flVFUr+h613qY313cfjBjvg4gGyOS7mNsjsfMLOA1AcV27ryISYlQwGQOpa/ZjSvQUsoqymBc3jxZ7Cw9te6hXAjcndg77Kvf1uUoFeJK7w9WHcWtummxNxAbGcqr+FO0O+fIlBEhiJy7TiZoydEbPsrx+Oj8m1ZV2HwwfoV2xAGOmM9cv0rurCygGXQeF0mInxBWLDQrh7c8MLE1ZQox/DJ+d+eycBC7SP5LSltLzXMGjsLmQZanLiPSPZMWGFfj7+VPRVgF2+QImhCR24pLZnC4qrCe8+5ODx2N0d06cCI4HY4CPIhsAShE58Q4mdc3wVRqGwFyK69rQNM23sQkxwoVZjKyansqOww0kGQLIbcg95z21HbUkBiX2e52U4BReO/Ea+Y35PDL9ER7Oehh/QwBaaxXYpDyRuLpJYicuWUFNG6rH+LpZ/jHdByPThj6ggTZxBYvaO7y7hsAc2uwualplDI8QV6Krzt1N16TS5DCSEZZ+znv2Ve5jduzsfmvczR8znyBjEAF+ASxLXcYLi19gT/ke7JZItI7GwX4MIYY1SezEJTtV1YLev8f4Oq3H5Oqoc/+jHnHiZ7LI3v2nYQw8CcpJUa108whxpbrq3EWEh7MwcdE5CZxbc7P++Pp+a9z9/sjvWXtwLd/c/E1yG3JJDEpkaepSGp3ttKFo7rBLfTtx1ZLETlyyE5U13uW2QHFta3P3wdHQYqf3Iy1+HvGddfk0vR19QAFFdTLOToiBFGAI5vnMF89J4LaXbifWEsu7d77LY7Me4960e1kzYw0v3/IyHxV9xJbiLYBnIsWj2x8FwOV24dJcOM2BtFid7C+S4sXi6iR17MQlO1R9FKU85U1izMmE1OV3HxwNLXaAGp/Jor2f89cQT+FTQ1COlDwRYoAFmwOYEjGVD1dtZEfZTnLrc0kPT2dyxGQKmgrISM1gauRUkoKT+GvOX3nu0HO4NXeva9jddrYWb2VpylKq2qow6o3Y3XZSY2NoszdhMoT76OmE8A1psROXrKj1mHd7SsS1UJvXfXA0tNgBjF1I5lnj7IpqW30XjxCjVLA5gFhzOHfGXc93rvkO6WHpVLRV8LM9P+PJT58kKTiJ7aXb2Vu595ykrktRcxHFTcUcqTmCS3MRGxDLP/L/QU5DDs1WGXMnri6S2IlLYnW4aKE7kbspagLYOxMecyhYonwU2QCLTGO6LogQl6crR+fXTG7jiQucJIS4LAYz/pYoQpQfMZYYthVv4+VbXiYjIoN6az1pYf1/YUwMSuR042lWjl+J3WUnpy6H28bdRkJgAh0umfQkri6S2IlLklfVhM5c4t2fa+ixjE9UOijVx1kjkFIYEmZzU49WuyrnASl5IsRgMZgJDIwm2BjMygkr+c6W77CnfA9/OvYnMhMz+50lOyt2Fjqlw6k5aXe2kxGeQVZJFqcaTlHTUUOblEARVxFJ7MQl2VVyFKX3FBQ1EkZcU2X3wVEyvs4rcXavsidu/2NS8kSIQRZoDGR6cyMbl7/JkuQlGPVGFIrnM58/7yzZ13JeY07cHGraa3jtxGvUWeuIDYxlW/E2DDoD9bZGbO3157mjEKOLTJ4Ql2R/5UHvdpxpEqq6R/dk9GQfRDSIEuYwv8OK0a1h1yn05iqyy05z+8Qpvo5MiFHNFDuF2Bdns/y721mYuJAWewsBfgG8c+c7bC3eSmlLKYlBicyKncVrOa+xOHkxpS2luNwuFiQsILsym+zKbL437Xvsr9xPZlImbbgxtdWBJcLXjyfEoJIWO3FJ8lt6TpyYCtXHuw/GTPJBRIMofgYB6Jhn7V47Nqtkuw8DEuLqoPxD0b70VwL/cBNh1lYMOgNR/lH89tBvWZi4kIWJCzHoDByoPMBXJ32VbcXbSAtLo6yljLSwNCx+FsLMYRyuPoxLc/HS4ZfIaThJOxo4rRcOQIgRTBI7cdE0TaPB1b0E0A3xM6H6ZPcbRluLndECMZNZ2N5dmPhI/Wc+DEiIq4QpCJU4F+2hA5iqcvm8/HM+Lv6Y5eOX8+UPvsxfc/7KyfqT7C7fzXe2fIcV41aQVZzFjNgZfG/r95gTO4fa9lqKmouYFTOLyrZKHsl6hGacYJPZ7WJ0k8ROXLSCxlLcek8xYs1lYlFYMDg7x6AFxozOLo4x01nYcwKF/QRNtiYfBiTEVcIU5Gm5S5jBqYZT/PrAr7E5bLy38j1uSrwJs97MTYk3seHODRQ0FhAbGMv64+uxuqzsrdgLChICE4ixxHCg+gB2t50dpTto0OmwOWWsrBi9JLETF+2j/D3ebaMzFUtDj/p10aOsG7bLmGlEudxMtXZ9ELjZWbbTpyEJcdUwBaHMIWSEeyZmrflkDdtLtnNLyi3cl34f0QHRfFz0MdfGXNtrRYpTDaeICYhhVuwsthVv4x+r/sG8uHnk1udS1nqG/MY8bO210i0rRiWZPCEu2r6KA97tGNNEqM7pPhgzyrphu4yZDsDC9g6Omk0A7CjdwYpxK3wZlRBXD2MACxMX8XT2f3ta3cp2YHVZ2VO+h6iAKGraa1h7aG2v4sUJQQlcG3Ut64+vJ8gYhJ/Oj7lxc8lMzOTjoo+5Lv46Gq2NxFhbwBIJpiAfPqAQA0ta7MRFy2v6wrs9MexaqOzeH7UtdtGT0HR+vcqefHpmF3bX/2fvvqPjLK/Ej3/vzEijXpEsucq9gm1wIVQ3IKZsqIEAJoQkhEAShGG6RQAAIABJREFUsim7SSCNQPa3JCE2JGwCCRBMQoAkQGyKKcY09wbG3XLHapatLk19fn+8o5mRLKtYGr0q93OOz7x97hx7rKun3MdrY1BK9S8Z7gx+N/d3xDviWVu8lml509hUuomle5eesCJFvCOeSwou4bkdz7HswDIGpw7m1X2vsmjjIq5bch3DM4ZTVlfGtuM74dXvYQ6v15Y71adoYqfapdJTSWXgMADGODhnyFQo+ihyQf5kmyKLMZcbGTCBkT4fQ3w+AOr9dawtXmtzYEr1H26Xm4nZE3nxcy/yranfYnv5dh6a9VCLde0evPBBFm5YyOsHXg8XL15Xsg6w1pX94fs/ZGjaUOLj02Dspcjq32O0xp3qQ7QrVrXL5tLN4e1gwyCmZLmg8pB1wOnue8WJo+VPQYo+YnZdPU+nxwHw4acfct6g82wOTKn+I82dhlOcfG7U59h5bCcN/gb+8R//YE3RGnYd38WglEFMy5vG4q2LWXZgWbh48eKti5u06HmDXtYWr2V63nR+7nuLcRMvZlaggQy/B7fLbeMnVKpraIudapdVR9aHt019AQXePZGTAyaAM86GqLpJaJzdOfWR7prVRavtikapfis5PpnUuExGpExmSs5USmtLKa0p5e6pdzN7yGw2FG8gNT6V7037Hk989okmEyqi7anYQ2pcKtXeWu7/6BEuXXodm8s262xZ1SdoYqfaZc2RphMnXKVR4+v6ajdso4FTAJja4MEVWit2T8UeyurK7IxKqX7J7XKSlZxIIBDPkLQhLN6xGF/Qx5HaI6wtXsug1EHkJOZwy2u3tJjUgVUG5eerfs7FBRdzybBL8Aa93PnWnVR4jnfzp1Gq62lip9rkCXjYVx1ZOmxS1hlNx9flnWFDVN0odwJBRxxJxjC1IfIbvbbaKWUPt8tJbmoKGa5UHpm9kN+s+w3jMsfxk8/8hERnIqMzR+OSlkcaNY67W1O8hhd3v8h3zvoO5w48F7/xs+LgctCJUaqX08ROtWnr0a0EjB+AoOc0Ts8fDEciY+76fIudy00wx5r1e3ZUd+yaojV2RaSUApISUpmUdTrfPuubVDQcY2/lXoamDiXTnRmeRRutcdzdppJNPPXZpzgr7ywe2/IYFw6+kOcvfx4Q8Nba82GU6iI6eUK1aVPppvB2oH4Yk7ICcKzQOuCMh7zTbYqs+zgHTYWSjzi7oYFHQsd0ZqxS9ktJTCOFNAYA+d58Kr3V7Dy+k8PVh1l61VLeOPAG+yr3MSR1iFUmpWQTA1MGctvrt+E3fmbkzcAT8PDe4fe4afxNNBzdTcKQGXZ/LKVOmbbYqTZtKNkY3g7UFTAxGLXiRN7p0A9mkklonN0Ej5f4oABQVFtEUU2RnWEppaIkx6fiFAcltSWMzhzNJ+WfkJOYgzfgpaKhghd2vsD5g8/nng/uYfaQ2Tw9/2lm5s8k3hnPWXlnkeZOoyp9INQf19p2qtfSFjvVqkAwwIaSyMSJZDOa9PKobtjB022IygahxM4FTGjwsznJCcCG0g1cnnK5jYEppaJluDMZlDqIoWlDWVK4hBvH38jg1MGsLV5Lla+KlUdW8sylz1Dtrea212/DG4yMqYt3xPPb2b8lxQdJxsDQmeBKsPHTKNVx2mKnWrXr+C7q/NaYk6AvlXGnjUA+jSR6DJpmU2TdLHcCfrFKupzjqQ4f3lSy6WR3KKVs4Ha5mZIzBbfTzfmDz+dw9WG+9PqXWLRxEf/a/S8eXPcgN75yI8cajjF7yOwm9/qNn2e3P4s3/3TY/z5o4WLVC2lip1q1viRSvy5QN4LxA1Lg08gxBveTxM7lpiptDABnRs2M3Vi68WR3KKVs4na5SY1PZVzWOL737veatMpBZAWKBRMX4BDrx+Alwy7h6flPc1beWSzcuIi/5Q/niMOwq7SM4sp6PP6AHR9FqQ7TrljVqvXF0YndcGamlkJDpXUg6TTILLAnMBsEBpwBlVs53ePFYYSgGPZU7KGioYKMhAy7w1NKNfPe4fdOSOoaeYNe1hevZ/qA6WS4M7i44OITumZ/s/4hFs1ZxK6SQWQmpjKtIBO3y9ld4St1SrTFTp1U0ATZUBrpdg3UDWeS75PIBcM+AyI2RGaPhGFnAdbYm0GeyO9Em8s2n+wWpZSNdh7b2er5Q9WHyE3KZcHEBfzw/R+22LJ39/K7OXO4m3HJNdTXaSkU1fNpYqdOqrCikEqP1ToX9CdjfLnkHY8aX1dwvk2R2SNleKTbeVpDfXh7Y4l2xyrVE43LGtfq+SGpQ8hOzGZd8bpWW/aWH1xO2tG1pJZuxOepi0WoSnUZTezUSTUdXzecUael4Dq0MnLBsHNtiMo+kjsRX2j0wqyopYeiWzWVUj3HrCGzTihS3CjeEc/0vOmU15dzqPpQq88prCjEN2I2zuduIFh7TMfbqR5NEzt1Us3H1100oBJqQ+ujJmZC7gSbIrOJK54i9wiAJkuLbSvfRr2//mR3KaVskuHO4NF5j7a4AsUjcx8hy53FrCGzGJE+otXnjMgYwa7KQhg8A8/WV6ht8McybKU6RRM71SJjTJP6dYG64cxybY1cMPQccPS/fz7HMyYCkBkMMsCkAOAP+vnk6Cet3aaUskFj6ZNXr36Ve8++l+vGXMe9Z9/L0qteISM+nQfWPMCh6kOcP+j8Vlv2Zg+ZzdJ9yyB1AEnHt7O9uLrFa5XqCfrfT2bVLvur9lPeUA6ACSQS9OQxriZq0ftRc2yKzF4mf0p4e1xDZOJIdBKslOo53C43A5IHcP3Y6/nJZ37C9WOvJzUuG19DLvee/WPAWh7woVkPtdiyt2j2It468BajUwZCdQm1mePZVVRlx0dRql1sLXciIvnA3cBMYBqQAsw2xqxox70CfAe4HSgAioGngPuNMb7YRNx/RI+v89cVkOYMkFoctej9qItsiMp+aSOmQ2gS7MyaY7ybaBUtjl5PVynVs6W44xiZk0GNx8P03LN5fvezHKs/xpOffZJ1xes4VH2IMZljuHDIhby5/00e2fQIr857HA7/GM8Vj3FxfJBntj3Dnoo9jMoYxZyhc0iPTyc5Ptnuj6aU7XXsxgL/DewBPgbO6cC9DwHfBl4AfgtMAH4EDAFu69ow+59VR1aFtwN1I7g2ex9SFVo7MWccZA6zKTJ7DRxzJh7jwi1+5jSU8iCDANhcuhl/0I/LYfdXSinVHinuOFLccWQmjeKbKd9kxeEVvLznZWbkzeDS4Zeytmgtl794OQ4cPHrO/WSsepRjX36PXXV7uPuVu5vMov3tht/y8JyHmZIzRZM7ZTu7fwptAE4zxpSLyJXAi+25SUQGAd8EnjTG3BZ1fBfwiIg8bIzR4mKnKBAMsKYo0joXqBnD/AHvRS4Y3T9b6wASEhLZ4SxgXHAPA/0BTnOlcdRfRZ2/jl3HdzEhu59NKFGql3O73AxwWV21voCPSm81y/a+QWHVLn5w5neYlTcDd8k+6s75MQ0pcdz94hdarHf3reXfYulVSzWxU7azdYydMabaGFN+CrfOBJzA35sdb9y/vlOB9XPbyrdR5bXGkAR9aYg3mzOq3o1cMGa+TZH1DIdSJgMgwLhAYvi41rNTqneLc8ax7KMqXllZwDBzC+cOvIYtpelUZM0gLT2D5QffbrXe3dsH3+7miJU6kd0tdqfKHXptXmOisXLkmS3dJCIVbTw3vTNB9RWriqK6YWtHcbZjO25PaDHs1HwY+hmbIusZqvNmQtU/AZhSVcEHqdbxjaUbuXnCzTZG1vfpd1jF2idHqlhZWM7KwnLW7T/Gzy4dRXrZBpyfHmSPZ1+r9xZWFHZTlEqdXG+dFdu4TkzzCrmNSyEM7MZY+pyVRyJFiP21o7kxaV3k5IQr+2WZk2jxI88jaKwZsRdUHgwf31iyEWOMXWEppbrAxIFpAFx2ej7zJ+WTbqqIe+56OFbIqIyRrd47so3zSnWHXvkT2hizEVgD/EhEbhWRAhGZD/wR8AGJJ7kvo7U/QGX3fYqeqdZXy0dlH4X342qHMCcYacFj0tU2RNWzDB8ymG3Gmjwy1uch1ZkAQHlD2xXsVefod1jF2rzxA0iMc/CV84ez9UglwR2vgd8Dq37PnDZWspg7dG43R6vUiXplYhdyDdZM2ieBfcAS4HlgE1BjY1y92vri9fiDVlX1QEM+V7CFxGBo4euskTB4uo3R9Qyjc1NZa6xJEg7gDEdq+JzWs1Oqd8tMjueFO85hVWE5gzISiSsLFR8P+knf8hIPz17YYr27h+c8TFp8mg0RK9VUbx1jhzHmU+A8ERkN5AG7jTHFInIE+NDe6HqvpuPrRnOjM2ow8LQvgUgLd/Uv8S4H+zNmQvWrAEypOMqHoYlwm0o3cdXoq2yMTinVGW6Xk7F5qfx19QEqEl14cyaREDqX/Ma9TLn4fpZe+TJvH3qHwoq9jMwYyZyhc0mJS8XlcFFSW8KKQyvYcWwH47LGMWvILDLcGbhd7lbfV6mu0ptb7AAwxuw2xrwfSuomAPmATk06RdHj6wpqE5ni2GvtON0w+Uaboup5gkPPpc5Y/1HPqCwOH9dCxUr1fnFOBxMHpfHEB/twjJ0PUUlZ8hv3kv/bydx8aCc/jR/GDcMuZ/6vN3OwvIFNpZu49F+Xcv+a+/nH7n9w/5r7ufRfl7K5bDMev6eVd1Sq6/SKxE5ERopIq6NSRcQBPAiUAn/rlsD6mOLaYvZVWrO+TNDFt3yR1Sc4/VpIzrYpsp5n7OAcPghOAmCSx0O8OAFrKbaj9UftDE0p1QXmjR+A0+Hg2U9q8F7/XJPkjqAf1j2GN30Yf91QTLUnSIBa7nr7rhZr3N351p1UeNqa0K1U17C9K1ZE7g1tjg+9LhCR84AKY8zvQscaW+AKou77PVb8m4F44EZgKnClMUYX8jsFKw6tCG8n1OdxmaOx9U7g3LttiamnmjgoneeCZ3KxcwPxwKSgk40SAKzZsRcXXGxvgEqpTslMjueJW6dz21Pr4OLhfOGujZidrxFXuhVfzkRk7Hye/aSGn75ayLmjsllb+n6rNe5WHF7B9WO1xKqKPdsTO+AXzfYbV5I4APyOk9uAtc7szVgzYVcC5xtj1rRyj2rF8oPLw9tX1B6PnBh3GeSMtSGinmt8XhrvBKcSMIJTDGdVlrExwyqhtrZ4rSZ2SvVybpeTaQWZvPv9Wby1vZSX9xiSki7jWM7FHDxWy5O//gh/0Lo2NzWBw7V7Wn3ezmM7Wz2vVFexPbEzxrQ5Gt8YU9DCsSeAJ2IRU39U5a1iXXGkXt3XGnaFtgRm/cCeoHqwxHgnOflDWV02gXOdW5lZ38DjocQuejk2pVTv5XY5yUtP5OazrfJGz6zez49f3nrCdaXVDUxIHtXqs8Zm6S/Hqnv0ijF2KvbeP/w+fmOVORneAHkBq1uRyV+AvNNtjKznml6Qxb+D5wAwxeMhPlSbeH/VfkpqS2yMTCkVC/PGD8DtOvHH5qrCcmbknt9qjbtZg2fFODqlLJrYKaBpN+zlddYgX+N0w5x77Aqpx5tWkMlrgel4jRO3gakNDeFza4vX2hiZUioWGsfdNU/u4pwOPJ4EfnX+wy3WuHt03qNkuDO6M1TVj9neFavs5w14+eDTD8L7c+qsJXjl7K9D+mC7wurxphdkUUUKbwXP4lLnWmbWN7Am0ap4tbpoNVeMvMLmCJVSXSl63N1rnxSzs7iaYdlJnD0imz+9vw8kjj/P/RfrSz/g07pCxmWPZdZgrWOnupcmdoo1RWuo89cBMMTnY6TPR70rncTz/tPmyHq2AWkJDM1K4tmKOVzqXMvZ9Q08HDq38shKgiaIQ7RRXKm+pHHc3ecmD+R/juzg/d1H+dWynQRDQzFe21LG+aNHs+iG68hIarlrVqlY0p86iuWHIt2wc+rqEaBqxrchUbsO2jJzeBYfBCdxKJjDBK+XzNDYxKP1R9l+bLvN0SmlYiU5wcWVUwex4cDxcFIHVrfs7ReMIDHeaV9wql/TxK6f8wf9vHPwnfD+nNp6jkguubPvtDGq3mP2uFwMDp4JzMMJnFsfGWf3/uH37QtMKRVT0d2y9185iZtmDOX+Kyfx7vdnMa0gE7dLEztlD+2K7efWFK2hvKEcgGx/gMkeDy8X3MHVcQlt3KkAzht9Gi6H8NfAXO5yvcz5dfUsTbEWjn3/0/e5Y/IdNkeolIqV5uVQonn8AY7Xenlrewlbj1QxcWAa88YPIDM5HvBT4anQNWVVTGhi188t2bskvH1pbS2FwUEMu/AWGyPqXdIS4phekMWqvYZnAvO4qX4pDmMIirClbAtH649yWuJpdoeplOpGHn+A9fuPc9tT6/A0VjEGfrF0O/+8czo17Dlh+bEH1z3Io/MeZUrOFE3uVKdoV2w/VuerY/mBt8L7l9fU8kT8TUwdpolIR8wdnwvAk/7Pkhh0MsVjLfZtME3KyCil+ofjtd4TkjoAX0DXlFWxp4ldP/b2wbepD1hJyEivF2/DYJLOuBKHo83FQFSUy87IRwTKyOCFwAVcXFsXPvfG/mU2RqaUssNb20tOSOoAPjOyfWvKKtUZmtj1Y0t2/Su8fXlNHb/1f57LJg+0MaLeKT89kXNHWq2cf/RfzpxaD2KsaXLritdRXl9uZ3hKqW629UhVi8d1TVnVHTSx66dK60pZU7ohvD+qKpPdKTOYOkRLnJyKq88cBMBBM4B1Moepoe7YIIa39r9hZ2hKqW42cWBai8dLqxsYrGvKqhjTxK6fem33SwSxWpWm1zfwgudzzD9joHbDnqLPTsojLcGai/SLms9xUZ0/fO7FLU/YFZZSyga6pqyykyZ2/ZAxhn9teya8/5lqJ68Hp3PZGfk2RtW7JcW7uHGmVfKgjAz8gYuID1Ut3VpfzI6STXaGp5TqRq2tKeskmd/P/b2uKatiRsud9ENrj6xkr/c4AEnBIHsq5zEwK4Upg/U/lM649ZwC/vzBXnwBw4NHL+HSlPdZlmh9xf75wX3cc82LNkeolOoO0cWL39peyvYjVYwfmMa88bmhOnZTefXqV1lxeAU7j+1kbJauKau6jiZ2/dBz6xeFty+u9vKCfw63Txmk3bCdlJeewLVnDebZtYeoIwFH/QWQuBKAf1ft4htFm0nPn2JzlEqp7tBa8WJwMsA1gOvHXt/tcam+T7ti+5lPqw+z/HjUGqbHp9GAm6vOHGxfUH3If84bQ2KctZTQCyWXUhCwvmJ1Dgd/ffPbYExrtyullFKdooldP/PkyvsJhBrmptV7WFL/OaYMyWD4acn2BtZH5KYlcPsFIwAwuKBiTvjcX/2lVG/918luVUoppTpNE7t+5Gj9UV4s+jC8X1A+igpSw6U6VNe4/YIR5KRa42S2lM4mN2htVzmd/N+HPwNfvY3RKaWU6ss0setHHv/wPryh1roJDV5er/48Lodw+RlalLgrJbtdfOeiMaE9J8dLLg2f+5vbsHP5T+wJTCmlVJ+niV0/sbdyL899+k54f+KxIRSRw6yxuWQlt1xTSZ26z08bwhmD0wE4VjGDwT5rZYqACD84uJS6oo/tDE8ppVQfpYldP2CM4f+99yMCof1p9Q28VXUDAJ+fppMmYsHpEH551elYE42FnQdvJj40b2JPvIufvP5lAv6W14tUSimlTpUmdv3AszueZdWxrQCIMUwvG8ZBk8+gjETmjh9gc3R916RB6dx27nAAgt48ssvnhc8tczTws5c/jz/oP9ntSimlVIdpYtfHbS7dzEPrfx3eX1BVzb/qrgPgxplDcWrtupj6z4vGMDA9AYBdZfM41xdpIX2pppA7ltzA8YbjdoWnlFKqj9HErg/beWwnd719F56gD4AxHi/DyiayxwwmMc7JF2YMtTnCvi/Z7eK+z00K77+552tc7ksI76+p2Mn1L1/N1vKtdoSnlFKqj9HEro96//D7fPH1L1LlrQIgKxDgf0qr+bXvRgBu+cwwnTTRTeZNGMBnJ+YBECCOraXf5es1vvD5ooaj3PzKTfx5y58JBAMne4xSSinVJk3s+phaXy33rbqPO9++k1pfLQDJwSCPFpfxfP2VlJJJQpyDr5w/wuZI+5f7PjeRjKQ4AD6uSqUq4Wc8Ul5NaiAIgN8EWLhxIbctu41Paz61M1SllFK9mCZ2fYQxhuUHl3P1y1fzwq4Xwsdz/X7+cqSEKs8Y/hyYD8A354wOF9BV3SM3LYH/uer08P4fdqdSMeoh/l5eyySPJ3x8Y+lGrnnpSl7e8zJGlx9TSinVQZrY9QGfHP2ELy37Ene/czdHao+Ej8+treO5T4vJ8KfzLc/XMTgYMyCFr2prnS3mn57PgqgFwb+90s2OGU/ztD+Lrx+vxBlK5GoDDdz74b18942v6cQKpZRSHaKJXS8VNEHeO/weX3njK3zhlS+woWRD+FxaIMD/lh7lt6VHcQRT+ELDDygjk6R4J7+9fgrxLv1rt8uPL5/AtGGZ4f07Xq/m2UlP8vVJX+HpkmMM80XG3r1ZvIqrn5/L2xv/iAkG7QhXKaVULyPa3RMhIhXp6enpFRUVtsVgjKGq3s+RynqKKxsoqmygrNpDgz9Ara+W4oYdfOrdQKl/Ew0cbXKvyxg+X1XDHRWVZAaD7AwO5iu+73LIDMDlEP644CytW9cDVNR5ueGx1eworg4fmzsul5+dn0T22vv5zdFVPJ+W2uSeyQEHX8w+i1nDLyEufwqkD4b45JjEl5GRQWVlZaUxJiMmbxBDPeE7rJTdevN3WHWeJnZRuvKHgjGGel+AGo+fqnoPFfV1VDTUU9Xgoc7ro8bjo9bro7rBT1l1AyXVdZTWVFNWW4HP1OJy1pLgOkZc3HHEVUHAXY43vqbF93Iaw0W1dXzjeCXD/H4CRngiMJ+H/NdSTwI5qW5+f+OZzBie1enPpbrG0RoPX316PZsORv6ticD5o3O4fkARjiML+V9HCeUuZ5P7UgNBzmpoYILXy1hJYFBiDjmpQ8jIGI4jYyhkFkDmMMgYBu6UU4qtN/9Q0MROqd79HVadp4ldlLZ+KPxy1a9ZsmMtQQIYE8AQJEgwvG3t+zD4MOIH8YEEEIlNN1pqIMhVNTXcVFXNQH8Aj4nj+cCFPBa4jENmACNzkrnmrMHcfPYw0hLiYhKDOnUNvgAPvLKdxasPtHh+sutjcnKXsiGtEr+0XkhajCHJGJKDQZKChkRjcIoDhzhxihOHw2ntI3xV0pmRkA/XL27xWb35h4Imdkr1/O+wiDwFTDPGTGrr2jaeMwXYBMw2xqzogtD6BE3sooiVgUl6enqL53PuyiNhjH0JktMYhvt8nF3fwIV19ZzV4MFhYGVZMksOpvDvAwmU1hoI+DANNRhPyy18qmdx5o4iYdpVxA2d0uL5FFcxQ7JepzJtL9VxnV9f9lelRzmztI4xj7ecLFZWVgIYY0yvG4zZ1ndYqf6gp3+HRWQkkGyM+biTz9HErgWa2EURET/WhJIqu2M5icafVpW2RtG1+uJngt79udKAoDHGZXcgHdULvsOt6c3/ZlrSlz5Pb/ssvfY73BGa2LWsT/+ld1RP/xKISAVAT21ePxV98TNB3/1cPV1P/w63pq/9m+lLn6cvfZaeILorVkRuBZ4EzgR+CVwAHAF+YYx5utl99wLfAFKAN4A/tPBsAb4L3A4MAz4FHjHG/DZ03gWsBvzAucaYQOj4D4CfA9M725Jotx7ZTKuUUkqpfuWvWMnalVitcE+JyPjGkyLyDeAXwGLgGmAv8OcWnrMIuA/4C3AZ8BTwvyJyB4Axxg8sACYDPwo9ezJWUveT3p7UgbbYKaWUUsp+vzPGPAogIiuxkrJrgPtFxAn8EFhsjPl+6PplIpKLlaQRum8kVoveHcaYx0KH3xKRJOCnIvKYMSZojNkuIj/CSvjexmr5Wwv8qhs+Z8xpi51SSiml7PZG44YxphY4AAwOHRoMDARebHbPP5rtzwu9/lNEXI1/gLeAPGBI1LULgQ+B5cBw4IvGmD5RCV4TO6WUUkrZrXmNIi+QENrOD72WNrumpNn+aYAARwFf1J83Q+fDiZ2xZo7+HXADrxlj9nYm+J5Eu2KVUkop1ZMVhV5zmx1vvpTSMcAA52Elhs3tbNwQkYFYkzU2AdeKyBxjzPKuCddemtgppZRSqic7jJXcXUXT7thrm133dug12xizpI1n/hkrETwfeAZ4UkRON8b0xlJJTWhip5RSSqkeyxgTEJH/BywSkRKsrtWLgdnNrtslIr8HFovIr4A1QBwwBqvW3ZUAoRmyFwEXGGNqReRrwCfAw8Ct3fSxYkYTu16kL9ZQ6oufCfru51Kx09f+zfSlz9OXPksv9giQAdwF3Ik1IeIrwOvNrvsWVpfr14CfADWh/RcgPHP218CvjDErAYwxpSJyO/CiiLxkjHkp9h8ndnTlCaWUUkqpPkJnxSqllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFJK9RGa2CmllFKqzxARIyILu/B5t4aeWdBVz4wlTeyUUkop1eOJyCgReUxE9omIR0QqReQ9Efm6iMTbHV9Hich/iMhGEWkQkYMi8lMRcXX2uZ1+gFJKKaVULInIFcDzQB3wNLAVSAQuAB4BhgP/FaO3Xwz8HfB01QNFZD7wErAc+CZwOvAT4LTQ/inTxE4ppZRSPZaIjASeBfYCs40xpVGnHxGRCVgJXkwYYwJAoIsf+2tgE3BJ6PmISBXwQxF52Biz+1QfrF2xSimllDoZJzAPWBB6ddoQw/eBZODLzZI6AIwx24wxf2h+XESuEZGtoW7brSLy2Wbnh4nIoyKyU0TqRaRcRF5oPpaupTF2IrJfRF4SkQtFZF2oO3WviNzS1ocJJaITgD82JnUhj2LlZde09YzWaIudUkoppU51+VrtAAAgAElEQVRwrNa7oN4b+OWyrcUpu0qq08YMSK26ZGJeTWK880dZyfGLuzGUK4BCY8zqDtxzIXAdVrJUA3wL+KeIDDXGlIeumQ6cg9XNehgoAL4OrBCRCcaYujbeYwxW9/DjwFPAl4GnRGSDMWZrK/dNDb2ujz5ojDkiIoejzp8STeyUUkop1cSxWu+CLYcrFt6+eEOWxx9sPJzxv6/vyHhswVkLTx+cQXckdyKSBgwEXu7greOB8caYfaHnvAN8BHwB+F3omleMMf9o9n5LgFVYrWZtfb7xwLnGmJWhe58HDgFfAr7Xyn35odeiFs4VYX3eU6ZdsUoppZSK5qz3Bn7ZLKkDwOMPcvviDVn13sADdE+3bFrotbqD9y1rTOoAjDEfA1XAiKhj9Y3bIhInItnAHqACOLMd7/FxY1IXel4ZsDP6PU4iMfTa0mSMhqjzp0QTO6WUUkpFm71sa3FK86SukccfZNm24jRgVjfEUhV6Te3gfQdbOHYcyGzcEZFEEblPRA5hJVlHgTIgA0jvivc4icaE0t3CuYSo86dEu2KVUkopFS1/V0l1WmsX7C6pTiHSpRgzxpgqESkCJnXw1pPNYpWo7Uewuk0XYnW/VgIGa8xdexq+2vMeLWnsgs3nxO7YfGAlnaCJnVJKKaWiFY0ZkFqF1XLVotEDUmtoeYxYLCwFvioiM40xa7rwudcCfzHGfLfxgIgk0Mrn7iKbQ6/TgI1R7z0QGBx1/pRoV6xSSimlor1zycS8Grer5RTB7XJwyYS8KmBFN8XzIFZh4j+JSE7zkyIyXkS+dgrPDXBi69o3ifHYwdCM2R3A7SIS/V5fB4LAPzvzfE3slFJKKRUtkBjv/NFjC8461jy5c7scPLZg2rHEeOc9dH3R3hYZY/YANwGjge0i8pCIfFlEviEifwc+pu0JCy1ZCiwQkYUicruIPIlVFqW8jfu6wvexJmgsE5Gvisgi4EdYte12debB2hWrlFJKqSaykuMXnz44g+XfnfXAsm3FabtLqlNGD0ituWRCXlVivPOebq5jhzHmJRGZjJUQXQ3chTWDdHNo+8lTeOzdWMnpTViTFj7EKsK8rCtibo0xZqmIXA38FGusXxlwP/CLzj5bjDGdfcapvbHINOAerIw1F2vQ4mbgvujpw+181qvAfGCRMebbnYjJj9WKWdXWtUr1YWlA0BjT637x0++wUkDXfoedWLNfGwf6r6CbWurUqbHzP+6Rofd/HOsfSwZW1vyeiMw3xrzZnoeIyGV03RpxDkDS09PbM81ZqT6psrISeu8wDf0Oq36vi7/DAeDtLnqW6ga2tdi1RESSsBb5XW+Mubwd18cDnwDPAD+n8y12Fenp6ekVFRWn+giler2MjAwqKysrjTGxnhnW5fQ7rFTv/g6rzutRv5WH1mVrLA7YHndjVWj+dcyCUkoppZTqJWwfQyMiqVjVl7OBL2IVIbyvHfflAT8G7jLG1Im0VQ/Q+m2+jUu0+0apHky/w0op1TrbEzusmSzXhLa9wB+AX7bjvv/BWpPtmRjFpZRSSinVq/SExO7nwB+xqi0vwGq9i6PlxXEBEJEZwC3AhaYDgwTbGm8Qag3Q3/iV6qH0O6yUUq2zPbEzxmwBtgCIyDPAeuAprKU+TiBWn+si4J/GmA+6KUyllFJKqR6vp02e8AEvA1eLSOJJLrsKmAH8n4gUNP4JnUsL7Z/sXqWUUkqpPsv2FrsWJGKt3ZYK1LdwfihWQrq8hXNfCv2ZD7weqwCVUqpfKtkK25fC0Z3gSoDMAsidAAMmQEYBOHpUW4FS/ZJtiZ2I5BhjypodSwOuAw4ZY0pDx4YCScaYHaHLlgD7W3jki1jrvv0Z2BiruJVSqt8J+OGd++GDhcBJhjW702Dc5XDetyFnbLeGp5SKsLPF7jkRaQBWAsXAEKzWtsHADVHXPQ1ciNWKhzGmEChs/rBQuZNCY8xLsQ1bKaX6mde+D+ufaP0aTxV89DfY/m/4/F9g1LzuiU0p1YSdid0zWDNbvwVkAhXAamCBMeZdG+NSSinVaOPTTZK6nQln8ET1TJwEGSWfMlYOcUb8p6QGQiUGvTXw7I1wxweQM8amoFV/JiKGTq5E1ex5t2KVZhtujNnfFc+MJdsSO2PME0AbvwKCMWZWO5/XdoVipZRS7Vd/HN74cXj30MD5fHbvTZjm8+58hhcvg6nrfwCVhyDggTd/DDc+180Bq75MREYB/wVcBAwEGoCPgGeBPxtjvDaG1yEicgcwB5iJNXfgL8aYW7vi2TrSVSmlVMve/w00WC1xwYzhfKF0QTipu/WcAq49a3DoQuHxg/lww18JjZqBXa/Dfq1IpbqGiFyBVRrtGuAl4C7gXqAEeAS4P4ZvvxhrYueBLnzmD4B5wHasxRm6jCZ2SimlTlR3DNY+Ht59Pe9rHK6xtgekufn+JWO5/YIR4fPv7izDkzMJptwYecamv3ZXtKoPE5GRWK1ye4Hxxpj/NMb8yRjziDHmOuCM0LmYMMYEjDENHVkQoR0uBLKNMZ+l5Qogp0wTO6WUUif66FnwNwAQyD2dH+6IJHHfv2QcyW4Xo3NTGJadBECtN8DKwnKY/pXIM3a8Av6TLiKkegcnVsvSgtCr04YYvg8kA19urJgRzRizzRjzh+bHReQaEdkqIp7Q62ebnR8mIo+KyE4RqReRchF5Iao2buN1t4qIiT4uIvtF5CURuVBE1olIg4jsFZFb2vOBjDEHujhRDNPETimlVFPGNJkw8W7G56hs8AMwLDuJK6cMBKxqBBdPGBC+7s1tJTBwKmQMsw54KqHwne6LW3Wt2vIFVBzaz+pHX+Df33qK1Y++QMWh/dSWL+jmSK7AqnqxugP3XAg8DPwNa1xeAvBPEcmOumY6cA7wd6yJnH8A5gIrRCSpHe8xBngeWAZ8F2sS6FMiMrEDcXa5nligWCmllJ0OrobyPQAYdyr37I7Mbr1r9ihczkibwLzxA3j8/X0ArN5bDiIw8Ur4cJF1wfYlMLZJQ4nqDWrLF1C0aSF/vzErqtU1g7d+lsENf1tI/lRIzl4c6zBC9W0HYq1K1RHjsbpt94We8w7WRIsvAL8LXfOKMeYfzd5vCbAKayxfW59vPHCuMWZl6N7ngUNYpdu+18F4u4y22CmllGpq+5LIZtY8iuqtNoAhWYlcNXVQk0snD8nA5bAmTOwtq6Wy3gdjL4tccODD2MerupoTX90vmyV1Fr8H/n5jFr66B+iebtm00Gt1B+9b1pjUARhjPgaqgBFRx8Jj20QkLtSatwer5e3MdrzHx41JXeh5ZcDO6PewgyZ2SimlIoxpktg9WhLpVbpr1ijinE1/bCTEORmXnxre33K4EgZOAafbOnB8H9ScMCxK9Wyz2bEk5aTjI/0e2LE0DZjVDbFUhV5TW73qRAdbOHYcq24uACKSKCL3icghwAMcBcqADCC9K97DDprYKaWUiij6CCqtn1deVyrL6qxu2Pz0BK4+c3CLt0wenBHe/uhwBbjcMCiqwePQmtjFq2Ihn9Idaa1eUbYjBciPdSDGmCqgCJjUwVsDJzkeXfP2EeAerHFynwcuxqqRV0778qP2vEe308ROKaVUxK5l4c0VnIUvNBT7y+cNJ97V8o+MyUMiid3mQ6EVKIbMiFygiV1vU0TuuKpWr8gZV4OVcHWHpcBIEZnZxc+9Fqsw8HeNMf8wxrwJfIDVYtdraWKnlFIqYl9kRceldVYjSWqCixtmDD3pLVOiEruPwond2ZELDmpi18u8w7granC5Wz7rcsO4y6uAFd0Uz4NAHfAnEclpflJExovI107huQFObF37JvaUdOkyOitWKaWUxVsLh9aGd1cGrfF1N589jBT3yX9cjMxJITHOSb0vQGm1h/IaD9nRLXbFH0MwAI5e/fOyPwkQl/Qjbvhb81mxVlJ3w9+OEZd0DyfviuxSxpg9InITVlmS7SLyNLAVazWI87BmsD50Co9eCiwQkUpgG/AZrFp95V0SeCtCK2lMDu26gTNE5N7Q/mJjzCmvcqGJnVJKKcvB1RD0AbAzOJijpBPvdPClcwpavc3pEEblprDl00oAdpfWkD3iNEjJg5piq9Dx8f2QPTLGH0B1meTsxeRPhW9seIAdS9Mo25FCzrgaxl1eRVzSPd1R6iSaMeYlEZmMVaz4aqwlxRqAzaHtJ0/hsXdjJac3YdW5+xArsVvW2k1d5Brgi1H7U0N/wOoO1sROKaVUJ0V1wza21l05dSC5aQlt3jq6WWJ39ohsyB1vJXYApds0settkrMXQ/bfOPvrs7AmShRhdb92S0tdc8aYncBX2nFdi5MXjDEFzfYrgNtauLT5dU8BT7X2rKjjs9qKL3TdrcCt7bm2ozSxU0opZTkYKey/KjgBgJtmDsMYw0dlH/HWgbfYcWwHld5KJudM5tox1zIuaxwAowakhO/dUxIqOZY7AfaGVp4o2Qbjr+iez6G6UgB42+4gVPtpYqeUUgr8XjiyOby7ITiG0QP9rDr2d36w9mUO1xxucvmOYzv4x65/8OAFD3JxwcWMyY2UGdtVUmNt5I6P3FC6LabhK6UsmtgppZSCki0QsAbJ7zK5VOeswp/+Ab/f7D/pLQET4L/f+2/S3emMHhApM7a7NJTYDZgQuVgTO6W6hZY7UUopBYfXA+AFvpeXjvu0FQSJJHWpcalcM/oaFs1exMLZCxmcYhUr9hs/D6x5gLz0eNyhOndHazwcr/VCzrjI88sLwdfQbR9Hqf5KW+yUUkrB4XUY4N6cbPYlRcpbnJFzBjePv5k5Q+fgdkbqmk3MnsiVL19Jra+WfZX7eHXfUkblnsbWI1Zd292lNcwYngUZw6DiAJgAHCuEARObv7NSqgtpi51SSik4vI6ViQm8lpIcPnTn5DtZPH8x84fPb5LUAeQl53HrxFvD+//30f8xKicpvL+7NDSBIntU5KZje2MSulIqQhM7pZTq7+orMMf3sygzsoLEJUMv447Jd+CQk/+YuGXCLWS6rfXOi2qLSEgvDJ/b3TiBIrrEiSZ2SsWcJnZKKdXflXzCiqREtrvjARATx3/N/A4ira9lnhSXxBUjIyVMPvW/F97e0ziBImtE5AZN7JSKOU3slFKqvyv6mJejumBPT51PblJuu269ctSV4e2tlSvBUQfArsZadtGJXXkhSqnY0sROKaX6ucqijbyfmBje/+LpN7T73tGZo5mYbU2I8Ad9JGRsBaC02kNlnQ+yorti93VNwEqpk9LETiml+rk3SzbhdVjdrknebC4aPamNO5q6bMRl4e2UrO3h7T1l1ZAxFBrH6VUdBl995wNWSp2UJnZKKdWf+T28TVV4d0jC+W2OrWtu7tC54W1v3E5wWPXqdpXUgCse0odELj6+v1PhKtUWETEisrALn3dr6JkFXfXMWNLETiml+jFvySesT4gP7180ouPruQ5MGciEbGuVCUMAV8oO4CQzY3WcnTpFIjJKRB4TkX0i4hGRShF5T0S+LiLxbT+hZxCRbBH5voi8LyJlIlIhIqtE5LqueL4mdkop1Y9t2vc2DQ7rR0G2z8FlEzrWDdto3tB54W1X6idAVC27zOGRCysOnlqgql8TkSuALcA1wEvAXcC9QAnwCHB/DN9+MZAIHOii530GeAAox4r7HqAeeF5EftzZh+vKE0op1Y+9e3BNeHugN5vBmUmtXH1yc4fN5eFNDwPgStkJ4o2UPMmI6orVxE51kIiMBJ4F9gKzjTGlUacfEZEJwAWxen9jTAAIdOEjtwKjjTHhRFFEHgXeAn4oIr82xpzyYFRtsVNKqX5sc32kEaLAPa6VK1s3In0EI9Kt0ibi8OFK3k1RZQPVDT5rAkUjTex6GycwD1gQenXaEMP3gWTgy82SOgCMMduMMX9oflxErhGRraFu260i8tlm54eJyKMislNE6kWkXEReaD6WrqUxdiKyX0ReEpELRWSdiDSIyF4RuaWtD2OM2Red1IWOGayWyESgoKX72ktb7JRSqp/yBDzscDYA1mSJKVHdqadi7tC57N1iFSF2pX6Cv2Yie0prmJoxLHJRpSZ2vcXxhuML6v31v3z74NsphRWFaSMzRlbNHTq3JtGV+KPMhMzF3RjKFUChMWZ1B+65ELgOeBSoAb4F/FNEhhpjykPXTAfOAf4OHMZKqL4OrBCRCcaYujbeYwzwPPA48BTwZeApEdlgjNnagVgb5YVej57CvWGa2CmlVD+15cg6fKEZsMO8PqafMbtTz5s3bB6Pb3kcAFfqdijys7u0hqljtSu2tznecHzB1vKtC+9efneWN+htPJyxcMPCjEVzFi2cmD2R7kjuRCQNGAi83MFbxwPjjTH7Qs95B/gI+ALwu9A1rxhj/tHs/ZYAq7DG8rX1+cYD5xpjVobufR44BHwJ+F5HghWRLOArwApjTFlH7m1Ou2KVUqqf+mDHm+HtkR4XBbkZrVzdtvFZ4xmYPBAAcTbgTN5rjbNLyQVXgnVRQ6X1R/Vkznp//S+bJXUAeINe7l5+d1a9v/4BuqdbNi30Wt3B+5Y1JnUAxpiPgSpgRNSx8Dg2EYkTkWxgD1ABnNmO9/i4MakLPa8M2Bn9Hu0hIg7gr0A6Vstip2hip5RS/dTHJZvD2wPJ7nD9uuZEhLnDIjXtXKlb2VFcDSJNa9lVHOrU+6iYm/32wbdTmid1jbxBL8sPLk8DZnVDLI1FFlM7eF9LTcPHgczGHRFJFJH7ROQQ4MHqAi0DMrCSrE6/Rzs9AlwCfMkYs6WD957AtsRORKaJyIsiciA0aLFYRF4XkXPace/VIvJcqJZNnYjsEJFfiUh7/iKUUkoB+7xHwtsjU8d2yTOblj3ZyvaiCmtHJ1D0JvmFFYVprV1QWFGYAuTHOhBjTBVQBHS0Ds/JZrFG//byCFapkeeBzwMXAxdhlSFpT37UnvdolYj8FLgT+C9jzLPtva81do6xGxl6/8ex/tIygJuA90RkvjHmzVbufQw4gtX/fRA4Hav5cr6ITDPGNMQ0cqWU6uUqPZUcdVr/VbqMYWrBzC557uScyWQnZFPeUI7DVcMx/27Kay4kW0ue9CZFIzNGVmH9XG7RyIyRNVg/u7vDUuCrIjLTGLOmzavb71rgL8aY7zYeEJEEWvncXUlE7gJ+BvzWGPPrrnqubS12xpjnjDFXGGPuN8b82RjzG+A8rKbQu9u4/VpjzBnGmJ8YY/5kjLkb+CowEWj/6tVKKdVPrTiwIbw91utl2Kj2DClqm9PhZM7QOeF9V+onbC+q1ha73uWduUPn1sQ7Wl7MId4Rz5yhc6qAFd0Uz4NAHfAnEclpflJExovI107huQFObF37Jt0wdlBErgcexhpb9902Lu+QHjXGLjS1uLF/u7XrVrRw+MXQ6/guDksppfqcFXvXhbcneby4ckZ32bObr0Kx7Uhl0zF2VZ922XupmAgkuhJ/tGjOomPNk7t4RzyL5iw6luhKvIeuLdp7UsaYPVg9eqOB7SLykIh8WUS+ISJ/Bz6mgxMWQpYCC0RkoYjcLiJPYvX+lbdxX6eIyAzg6dD7vA3cJCI3R/0Z0Jnn217uRERSATeQDXwRqx/9vlN4VJfUf1FKqf6gsOyj8PbwQAK4U7rs2dPzpuN2JOMJ1uKIr2DtkS3cPmxg5AJN7Hq8zITMxROzJ7LkqiUPLD+4PK2wojBlZMbImjlD51QluhLv6eY6dhhjXhKRyVjFiq/GWlKsAdgc2n7yFB57N1ZyehOQAHyIVYR5WVfE3IoJQDyQAzzRwvnZWEulnRLbEzusv4xrQtte4A/AL0/hOf+N9Rf0r5NdICIVbTxDJ18o1YPpd7jrlHkPhH8CjEgc3KXPjnPGceZp57Kq9A0AtlV+CGl3Ri6oOnKSO1VPkpmQuTiTzL/dPOHmWVgTJYqwul+7paWuOWPMTqxab21d1+LkBWNMQbP9CuC2Fi5tft1TWAWIT/qsqOOz2hHfCc/rSj2hK/bnWDNRbsPKlt1AXEceICI3YlV8ftAYU9jlESqlVB9SUl1FrdOqIuEwhgn5HZ1w2Lb/GH1JeLvSsZGGxNzIyepiCPi7/D1VTASwugufCb3aktSp9rO9xS5Us2ULgIg8A6zHymSvbc/9InI+8GfgFeDHbbxXq2P3Qq0B+hu/Uj2Ufoe7xmu7NtPYpjHU5yc9b2KXv8fcgvPhg3gQLw53KW/u28UVSadB3VEwAagthbSBbT9IKdUhPaHFLswY48NaNuRqEUls6/pQf/u/sQZOXm+M0d8klFKqDSsPRmqgjvF6IXtUl79HoiuRXOfk8P7SwjeaJnLaHatUTPSoxC4kEWv6catVpkVkJPA6UApcZoyp7YbYlFKq19txbFd4e4zPB6eNicn7TM+9ILy95fgHkDYocrLycEzeU6n+zs6VJ1qqRZMGXAccMsaUho4NFZFxza7LA94AgsAlxhidCauUUu3Q4AtQ6Y0MRR4dcMSsS/RzY+ZhjFUSrNrs53BK1EpL2mKnVEzYOcbuORFpAFYCxcAQ4EvAYJoWGX4auJCmRQRfx6pZ8yBwnoicF3Wu0BizKpaBK6VUb7X54HFwF4f3x6UMttZyjYEzh+QTrB2NM2UHAK8E6whXkdWSJ0rFhJ2J3TPALVjFADOBCmA1sMAY824b9zYO3PivFs79BdDETimlWrCicA9BpweAlGCQ/Oxxbdxx6twuJ/lxMyjFSuxeqz0cldhpi51SsWBbYmeMeYKWC/M1v25WC8di8+ulUkr1cSsPbQn/zz/G60Vyxsb0/S4YNJsXyv6KSIBCfylHXE4G+gPaYqdUjPTEyRNKKaViwBcIsqdid3h/tNcXkxmx0S4cNZRAbeQ93kpKsja0xU6pmNDETiml+omPD1cSjIskVGO8sZsR22haQRaBmkidvPeTEqyN6iIIaoUqpbqaJnZKKdVPrN13DKe7KLw/1uuD7JExfc8Ut4vRKdPD++sTEqgVgaAfasti+t5K9Uea2CmlVD+xal8xDnckmRqdOADi2qwF32nnDR9JoCEfAL8IqxNDrXY6zk6pLqeJnVJK9QP+QJCNR3aCGACG+HwkxbgbttHMEVn4ayKzb99LCiWTlZrYKdXVNLFTSql+YFtRFQ2OyGoP3TG+rtG0giyCtZHEbnVCY4udTqBQqqtpYqeUUv3Ah3vKm4yvG9MNM2IbpSXEMS5zIiYQD8CROBefupzaFatUDNhZoFj1YEv2L2myf0XBFTZFopTqCu/uKsURteLEWK+321rsAM4enkPhgQJcKdY6tesSEhikLXZKdTltsVMnaJ7UNR5r6bhSquer8fhZv/84joRIYtedXbEAZ4/IJlA3Iry/LsGtLXZKxYAmdqpDNMFTqvf5cM9RAlKFw1UDQGIwyCBnMqTkdlsM04dnNUns1ickaGKnVAxoYqdOiSZ4SvUe7+4qa9JaN9rrw5EzFqT7VmdMT4xjTMb4JuPsjtSWQjDYbTEo1R9oYqc6RZM7pXo2Ywzv7izD0WTihBdivEZsSz4zIpdA/dDw/sfxDqg72u1xKNWXaWKnOk1b75TquQrLavm0ov7EGbE541q5KzZmjsgi0BBJ7D5yx2vJE6W6mCZ2qstocqdUz/PuLmulCUdCJLEb5/Xak9gNzyIY3WLndmtip1QX08ROdSlN7pTqWd7dVQbiw+EuDR+zWuy6vys2Iyme4amRhHK7Ox5v5cFuj0OpvkwTO3WCztasW7J/Ca9tK237QqVUTNV7A6zeW47DXYqINUlhqM9HclwypA+2JaZzhxeQ5nUD4BNhR/kOW+JQqq/SxE7FhD9pjSZ3Stls9b5yvP4gDneku3NsY/26bpwRG+3sEVmkNmSG9z+u2mtLHEr1VZrYqRa11Wq389jOJn9aosmdUvZ6d6c1vs4ZNb5urE3j6xrNGJ6Noz4/vL+1ocy2WJTqizSxUx3WUiJ3sgTPn7RGx90pZZP3GidORM2IHeexp9RJo6zkeNLjIyte7AjU2haLUn2RJnbqpE5lrN3JWu80uVOqex0+Xsfeo7WAadZiZ0+pk2ij82eGt/c5gzT4G2yMRqm+RRM71WFjs1r/bf9krXea3CnVfT7cYxX+lbjjiNNKnDICAQYEAra22AFMHjmOoV4/AAERdpdstjUepfoSTexUqzozQ1aTO6Xs8+GecgCcCU0nTogrETKGnuy2bjF9eDbDPJHJG9uOrLMxGqX6Fk3sVJs6m9w1T/A0uVMqtowxrCy0WuyajK/zeuG00eBw2hUaAEOyEsn3JYb31x/RFjuluoomduqUtNUd21xLyZ0meErFxs6Sao7WeAFISC4OH7d7RmwjEWGI47Tw/vaqfTZGo1TfoomdapfOFi0Gbb1Tqrt8sPtoeDs+qSS8PdZjz4oTLRmdVhDe/jRQji/osy8YpfoQTexUuzVP7jraatdIkzulYmtloTW+DkcdHqwkL84Yhvt8kD/ZxsgiCvLHku+3JlD4Jci+Sm21U6oraGKnOuVkyd2uY7tO+BNNkzulYsMXCLJmb+PEicj4ulFeH3HQYxK7vEHDrZp6IasOfWRjNEr1HZrYqQ7pTJdsW8mdUqrzPjpUQa03AEBmZqQbdoLXC6n5kJJrV2hNODMGMd4bSew+PPSxjdEo1XdoYqc6rDNdsq0ld0v2L4FtL3UuOKX6uXA3LJCZGZk4MdHjgbwz7AipZWkDmRDVYrfr+A4bg1Gq79DETnW75l2zTZK7JKcmd0p1woYDx8Pb9Y794e1JHm+P6YYFIGUAY0NFigGO+fcRNEEbA1Kqb9DETp2SrphI0Wpyp5TqsGDQsPlQBQDirKHSVwqAOxhklNcH+T2oxc4ZR05iNlkBq9vYiIftR/faHJRSvZ8mduqUdUUJlJMmdzqZQqkO23u0lsp6q2xIekZ0/brQxInB0+0J7CScaQMZH9Ud+8buTTZGo1TfoImd6pTo5C661W5M1ph2P0OTO6W6xsaDkW7YvJyy8PYkj9ZG7ZsAACAASURBVBcyCyA1z4aoWpE60FoNI2TdkS02BqNU36CJneoRNLlTqvM2RSV2zsTD4e1JXg8MOduOkFqXNpBx3khh4n3Vu1q5WCnVHl2a2IlIvIi0a5CEiEwTkRdF5ICI1ItIsYi8LiLntPP+QSLyvPz/9u48TK6yzPv4967eO71l3xdJIIshJCEkyCIIqIOAIoqCqKOi77jryMx4jTqjzjDMjBs4ivu+DeAG7gMqAZEdAiSEQMi+dZJO0p1eq7ur7vePU11Vaaq7q9dTVf37XFdddbancp9Uneq7nu2YNZrZcTO73cxeNLwzkKEYiVq73pTciQzek3uaEkvOsVjqq/jF0U6Yty6coPpTM/OEueya47voimkAhchwjHSNnQELsjx2IVAMfBN4P/BZYBpwr5m9vN9/xKwKuBs4F/gP4JPAamC9mU0cUuQyLCPVJKuaO5GhiXbHeO5gMwBW3ERzV1B7VxWPs6CrOzdr7KpnMa+7m8p4kMxZcQv379QACpHhGHRiZ2bb+3oAzwKezeu4+63ufpm7X+/u33b3zwPnAA3AhwYo/l5gEfAqd/+su98IvAKYDfz9YM9JRt5QbzcGapYVGYqtB1vojgdfvzOmHUpuXxbtJFIxCaYuCSu0vtXMIgIsTutnd9fzGkAhMhxDqbGbCtxIkED1fvzjcIJx9zbgMFA3wKGvBx509+Q3gLtvAf4EvGE4McjQ9TVKdihNskruRAZn076m5HJdXfrExJ2w6EKI5GCX6ppZQDBqt8cThzaFFY1IQRjKlf4EcNDd7+j9AH5F0BybNTOrNrMpZrbYzG4AlhMkaH0dHwFWAI9m2P0wcIqZVfZRtrG/B1A7mNjlhYbSJJvpXrI923souRPQNdyfTftTiZ2X7kkuL49GYVG/vVvCUz0T4IQpT/a0bMM9q4YfEclgKInd/wBH+tjXBbx9kK/3XYJaui3AdcDXgBv6OX4SUAYcyLDvAEFiOXOQMcgIGmp/u0wJnpI7kexs2nc8sRTnSNe25PblnZ2w8IJwghpIWRWU1Z4w5Ul38V72HG0PMSiR/DboxM7df+ruGWvU3D3u7t8f5Et+mqB/3DuAvxIkbSX9HF+ReI5m2NfR65je8dX19wCaMpWT4RnsYIreCZ6SO+mhazizWNzZUh8kdpGyejribQBM6Y4xc9qpUDU1zPD6VzOTRZ1dFCdq6SKlR/jLtj0DFBKRvmSV2JnZdDP7gJl93sw+bWav6au5c7DcfaO73+Xu3wVeCZwOfK+fIj0/5coy7CvvdYyEpL9bjmXb507JnUh2dh9to6MrGFlaO3Ffcvvqjg5s2eVhhZWdmlmUAiel9bO7Z+dT4cUjkucGTOzM7FzgeeAmggES/wL8EthlZh8cyWDcvQu4A7jCzDLWugFHCWrrMjW3ziQYlZupmVbGWLb3k822/52SO5HMnq1vTi7X1qSmC1kdjUKuJ3bVwQCK9ObYTQ2bw4pGJO9lU2P3ucTzO4B5BPPUXQVsA24ysx+PcEwVBP3kqjPtdPc4sBFYk2H3OmBrYnSt5ICRGimr5E6kb6nEzumMbEluX1W1ACbl+LztNYkBFGk1dke7dtDY1tlXCRHpRzaJ3XLgRnf/vrvvdffd7n6bu58JvAu4yszeP9h/2Mxe0OnDzGqAK4E97n4osW2emfWegOlnwJlmtiqt7GLgAuCng41FRtdlbbHkcl9Nstn2vUt/BiV3IkBqYuKSYzRb0BNlQjzOKUtfF2ZY2UlOeZJK5CLl+3ls17G+SohIP7JJ7JqB3Zl2uPu3gVuAdw/h377VzH5nZp8ws3ea2aeBTQS1gv+QdtwPgGd6lf0KsB34nZn9g5l9GLiLoAn2xiHEIqNp2eVZJXfZyJTciYx3zyYSu+qKVG3daR1RipfnQWKXaIpNn/IkUnaQB3fW91VCRPqRTWJ3N/CqAfYvHMK//SOgEvgg8FWCu0k8CbzM3W/rr6C7NwPnA/cR9Pn7d4L59c5z976mYpEw9ZPc9RhqkqdaOxnPot0xdjS0AjBnwuPJ7avKpsDE+WGFlb3aOQBUuTOvO9hkFueBPRpAITIU2SR23wTOMrO+bvO1ANg/2H/Y3b/j7ue7+zR3L3H3qYlbjN3T67jz3f0Fkx4nmoWvdPdad69291e7u24ymMt6deLuSe6GmtClLyu5k/Fq++FWYolbicUqU1/Fp8+/MKyQBieR2AGsiHYkl7c1bSHaHctUQkT6kU1i90eCSYG/YGbrzezNZrbUzE42s3cS3Nf186MapRSMgQZTDGcqFCV3Mh5tPxzU1k0uOsD+0mDKk2J3lq+6NsywsldRB2U1ACzvSI1787I9J9wmTUSyk01i92/A74A9wEsJ+rxtIrhTxNeBZ4FGM1tuZsWjFagUjr7uTNFjKMldDyV3Mt7saGgBYOWEu5PbllFGxaSTwgpp8BK1di9O72dXvpdHd2oAhchgDZjYufun3P1yd19AUHN3IcHghh8RJHgrEstPAi1m9qSZ/XD0QpZCkCm5G+woWdBIWZHtif51EytTY8xWT14eVjhDk0jslnR2EUn8WSoqO8wDO/f1V0pEMhjULcXcvdHd73b3G939b939NKCKYE65dxH0x2sGMre3iaQpbluXXB7p5C6dkjspZDsaWqmijV2Vrcltqxbn+KTEvSUSu3J3TiquS27ecHAjnrjVmIhkZ9D3iu3N3Tvd/fHEYIgPuPs5iXs2ivTr4mXTBkzuBitTfzuRQrajoZXTSzawpawUgGKHtQteHnJUg1Q7N7l4WklVcrnddrLtcGumEiLSh2EndiLDcfGyaRm3azCFyMCOtXbS2NbFjKpHkttWlE2mqrSqn1I5KC2xW96VqqEL+tkdDSMikbylxE5C19dgCiV3Iv3r6V/XOiHVF+2cmWeFFc7QpU15srzteHK5qGIvj+oOFCKDosROcsJoJHc9lNxJodrR0Mp0GthYmdp21tI3hBfQUKUldguPHaAkEjQrR0oaeeqABlCIDIYSO8kZIz0NikbKSqHb0dDC4opHOVpUBMAkj7B02oqQoxqCmlkQCWbLKmlrYHFd6lrf1fIMnd3xsCITyTtK7CSnaBoUkeztaGilompzcv3MytlELA+/1iNFJw6gqE4te+keth5qDiMqkbyUh98AUuhGKrnroWlQpFBtP9zK0QmHk+tnz8rD/nU9Ji5ILi4vqU0uF1Xs5un9xzMUEJFMlNhJThqJ5G6gwRQi+Swed/YeOcjW8lQzZV72r+sxcX5ycWUsdROjoordPL2vMYyIRPKSEjvJCyOR3PVQk6wUggPHO5hVfj8xMwAWx2DK5KHPARm6tBq72S1HqSmZDIAVRdlwcEtIQYnkHyV2krPSa+3SDbVZVv3tpJDsONxKWVXqNmJnlc8MMZoRkJbYWeMuVkw9Lbm+4/gm4nHdgUIkG0rsJKcNNA1K7+VMMjXJ9qbkTvLNjoYWGisPJtfPnrY6xGhGQF2qKZbGXZw5K3U+XSU72HusPYSgRPKPEjvJeSOR3GWi/naSz546uI3G0m4AKuJxVi24MOSIhimtxo5jO1mZVmNXVLmbzQeaxj4mkTykxE7ywkgmd2qSlUKw+dhDyeW1HVFKZ50eYjQjoGIilNUEy11tLC2bQoQSACKlR3hsz+4QgxPJH0rsJG9km9wNdlCFkjvJR02dDyeX18VKoXJSiNGMALMTau1Km/Yyp/Lk5PpjB58IISiR/KPETvJKNsldz3r6YzCU3Emua41GaSndmVw/uy6PR8Omm7wotXzkeVZNX5Vc3dm6OUMBEelNiZ3knWyTu4H0VWsnkuvu2v4QsUgMgNld3bxo5pqQIxohU1I1dDRs5bx5qfNqj2zjaGtnCEGJ5BcldpKX+kvuhjKQAtQkK/nj7t1/SS6f3d6OTV8WYjQjqFeN3eq0Grui8n1s3NcQQlAi+UWJneSt3sndUGrv+pr+BJTcSe7aeDTVv+6s9g6YuiTEaEZQr8RucsVkKm06ABbp5t5dT4YUmEj+UGInea33JMbDaZqFFzbJ/n7zoaEFJjJKjrQf4XDndgCK3TmjowsmLww5qhGSntgd3Q7xGAuqUrWRTxzSAAqRgSixk7w3UHI3mAmMe+uufEjJneSURw4+klw+NRqluGoeFJeFGNEIKq+BqqCGjlgnNO46oTl2T9szfRQUkR5K7KQg9JfcweBq73rX2nVXPtTHkSJj7+EDqWbYte1RmFIgI2J7TE4fQPE8Lz9pXXK1LfI8nd2xEIISyR9K7KRgZJPcZUrwskr6Nt8+rNhERsqD+1M/NNZ1dFA+68UhRjMK0kfGHt7CadMXQ7wcACtu5qE920IKTCQ/KLGTgjJQcgcnJnh9JXW9a+1+XVmk5E5CV99az56W4A4MZfE4p3VEiUwrkIETPaanJaqHNlMUKaLGUn0I/7xDNegi/VFiJwUnm+QOBj+44teVRepvJ6F6uD7VDLsqGqUUYGrmz3fempY2dcvBTQC8qDqV7D3ZoAEUIv1RYicFKVNy11eC15dMkxarv52E6aEDac2w7VEcO7FPWiFIn5Pv8LMQ62LN9NXJTXvadAcKkf4osZOC1Tu5g75r7wZD89tJGNz9hBq7tR0dtFbMgtLKEKMaBRUToWZOsBzrhCPbuGjhGtyDP1cd7KO5sznEAEVymxI7KWiXLbgs66bZTHSrMckVe5r3UN9aD8CEeJxl0U66JxdYM2yP9H52BzexdPo0iM4M1s35y+5Hw4lLJA8osZNxYTjJXSaqtZOx9lB9qhl2TUeUYqB05tLwAhpNvRK7oohRG0n1iV2/65EMhUQElNjJODLUfnd91dopuZOxlD5/3br2DgAqZhXIPWJ7m7E8tbw/GCyxqCa17amGDWMdkUjeCC2xM7MzzOxmM9tsZq1mttvMbjGzRQOXBjO7yMzWm9kRMztmZg+Y2RtGO27Jb331uxsowVNyJ2Fydx4/+Hhy/YyOILGzqQVaYzd7TWp532MQj3PGzNOTmw50PEdXvCuEwERyX5g1dh8FrgD+CHwI+AZwPrDBzPr9tjKzS4E7gWLgk8C/ADHgVjO7dhRjlgKQKbmDkRlYITIa6lvrOdQeTLVTGY9zcmciqZlSYCNie9TNgwnTguXocWh4jrVzX0S8cxIAcTrZcmRLiAGK5K4wE7svAPPd/YPu/i13vx44FyghSPr68z7gAHChu3/Z3b8MXAjsB946mkFLYegvuRtsgqdaOxltGw6lmh5XRKMUAc1l04N7qxYiM5hzRmp97yMsmVlNrH1+ctMj9Y+FEJhI7gstsXP3+929s9e2rcDTwEDtCzXAMXePppWNAseA9pGOVQpTX8kdDC3BExktTxxOTcq7siP42myvLdDauh5z0ptjH6WmvIRqUuf8170aQCGSSU4NnjAzA6YDDQMceg/wYjP7dzNbmHj8O3AK8PnRjlMKR3/JHaQSvIGSPNXayWh64lBaYhcNfs8W3K3EektP7HYHI4JPrj01uenpo0/h7mMdlUjOKw47gF6uAWYDHx/guP8AFiaO+0RiWwvwane/q69CZtY4wOvWZhmnFJDLFlzG7zcf0l0l8sB4vIbbutp47thzAJjDio4gsauau7y/Yvlv9ulQVBpMUnz4GWg+yOkzFvPUvnKsqIPW7kZ2N+9mfs38gV9LZBzJmRo7M1sC3AzcB/xwgMOjwHPAT4GrgTcDjwO3mdkZ/RUUyeTiZdMobls3rNdQrZ2Mho0NG4l5DIC5nVCdqKUqn1mgU530KJ0Ac9Ouye3rWTarjljbguSm9JHCIhLIicTOzGYAvyXoI3elu8cHKPIl4FXA1e5+i7v/GLgIqAdu6quQu9f19wCaRuaMJB+NRHIno2s8XsPpzbBroq2pHVNPyXB0gTnpvNTy9rtZNqvmhAEU6YNKRCQQemJnZrXA7wmaUF7p7vUDHF8KvBP4TXoC6O5diddZa2a51sQseeLiZdMG7HcnMpbSB06cHg3mrztePDm4p2qhO+mC1PK2u5k3sYKKeGqq04cPqMZOpLdQEzszKwd+TTDo4VJ3z+bGnJMJ+gYWZdhXkthnIxakjEtK7iQXxD3Ok4efTK6vTPSva6nJah73/DdrJZQnuk221GP1T7Ji6nLcg6//fa27ONpxNMQARXJPmHeeKAJuBV5C0Pz6YB/HzUv0v+txCGgErjCzkrTjqoDLgE2J2juRYVFyJ2Hb0bSD5s5mACrjJczt7g52TB0nU/FEiuCUi1Prm3/F6rnTiLfPTm5Kb6oWkXBr7D4PvJqg+XSSmb057XF52nE/AJ7pWXH3GPA5grnuHjCzD5vZdcDDwBzg+jE7Ayl4Su4kTOl9yBa0FyebIqrmFPiI2HTLXp1a3nwHK+fUEmtfkNykfnYiJwqzL9rKxPNliUe6XcDtfRV09/8wsx0EtyL7JFAGPAVc4e6/HIVYZRzrSe406lXGWnpt1GmJ+8MCVM99cRjhhGPhBVAyAbpa4eg2VpXtJdY2P+iUAzymkbEiJwgtsXP384dznLv/BPjJCIYk0q+BEjzV7slIS+9fd0FHat52mzrQzXkKSEkFLP4b2PRzAOq2/oIZZeuSw583H9lMR3cH5cXl4cUokkM0elRkkJTAyVg41nGMncd3AlBEMau6gqlOWorqqJowOcTIQnDam5KJHU/ewulzLueulqkUlR0m5t1satjEmhlr+n8NkXEi9OlORETkhdKbYaf5FMoSd89qrl4YUkQhWvgyqJ4VLLc1cFn5U5rPTqQPSuxERHJQ+vx1c9pLUzumFvg9YjOJFMHKq5Oraw///IQ7UCixE0lRYicikoPSa+yWHE/dcWLC/JWZDi98p78NLPiTVXPgryyIpvrUbTi0gVg8FlJgIrlFiZ2ISI6JxqJsbNiYXL+gI3VDnuoFq8MIKXx182DJpcnV68oeJd5dDUBLVwtbjm4JKzKRnKLETkQkx2xq2ERXPJhnfVrZbFbHg8QuRgSbtizM0MJ15nuSixd1rqeoNdXP7oEDD4QRkUjOUWInIpJjHjv4WHJ5HtOJWDByoqFsHpRWhhVW+Oa9BGasAKDEO3lJe3ty14MHMt68SGTcUWInIpJjHk+bdHdeS+prunniOK6tAzA7odbufdGnkssbDm6go7sjUymRcUWJnYhIDumOd58wInbpsdRN7iMzV4QRUm5Z/jqYMA2AFfEjVEarAOiMd2p0rAhK7EREcsqzx56lNTEZ8bTKaaxueT65b+LJa8MKK3cUl8HadyVXz+1oSS6rOVZEiZ2ISE5Jb4ZdXLWURewDoJsIExedGVZYuWXNtVBcAcCrOg4nNyuxE1FiJyKSU9ITuzkdFcmBE3tLF0LphLDCyi0TJsPKNwGwpqODxH8Rzxx5hqZoUz8FRQqfEjsRkRzh7jx+KJXYvehwKkk5NnlVGCHlrpe8D8eoiTunRqMAOM7D9Q+HHJhIuJTYiYjkiB3Hd3C0IxgsUVNaw5KG55L7iuetCyus3DR5IbbkEgDWdaRGw96///6wIhLJCUrsRERyRPr8daumnsbi6Obk+vQXnxtGSLntrA8ET+2pxO6ePfcQ93hYEYmETomdiEiOeGB/6u4JJ8dqqbSgiXE3M5g695Swwspdc9fRPm0VKzui1MWCe8Uebj98wu3YRMYbJXYiIjmgO97Ng/tTozqX1KdGe26rWYeZhRFWbjOj7KUfohg4vy11F4o/7f5TeDGJhEyJnYhIDtjUsInmrmYgmL9u+f5Hk/uiC14WVlg5L7L0Mg4Xz+TC1lRi9+fdf8bdQ4xKJDxK7EREcsB9++5LLp8zZSWzozsA6PQipq94eVhh5b6iYra+6M28pKOdinjQt27X8V1sb9oecmAi4VBiJyKSA/6y7y/J5TXNqdqnB/xUli2YGUZIeaPqJW+nI17JOWmDKNQcK+OVEjsRkZDtb9nP5iPBCNjiSDGnb01N2fHExFdQVlwUVmh5Yen8mdzGy7mwtS25TYmdjFdK7EREQnbXrruSy2dOXMqs5qAZtsXLKVp2WVhh5Y2SogjPzL2aM1s7KU70rdt8ZDMHWg6EHJnI2FNiJyISsvTE7qLGo8nlP8TXcs6yeWGElHeWL1nC3bGzWJfWHPuHnX8IMSKRcCixExEJ0YGWAzx5+EkAiizCBdsfSu67teRyVsyuDSu0vHLOyVP4dvfF/E1ac+yvnvuZRsfKuKPETkQkRLc/f3tyeW28hImJkZ3/F1vD7FNWE4lo/rpsnDytiiNVp1B1fH5ydOzzzbuTfRdFxgsldiIiIYnFY/zi+V8k1684vA+Abo9wU/frePmyGWGFlnfMjHMWTeFH3Zfy8rRau9ufvS3EqETGnhI7EZGQ3L//fupb6wGoi8W4IJGQfD12KfvKFnLh0mlhhpd3zl40hXviK1hzvCK57Xfbf0tnrDPEqETGlhI7EZGQfP/p7yeXX93SSimwMb6A/+m+gstOm0V5iaY5GYyzF03BifBA68XM7uoG4Hg8yvpdmvpExg8ldiIiIXji0BM8VB8MlChy56rjzRzwSVzb+Y9EKeX1p88JOcL8M6O2nEXTqrgjdg4XtXQnt9/+5DdDjEpkbCmxExEZY+7OVx+4Prl+SUsr1Uzi6s6Pc4iJnHnSJFbNmxhihPnrnEVTiFJKW+MZyW1/bdrKgeb9IUYlMnaU2ImIjLF7HvkS9zc+C4C589bOci5r/Tg7Pbh12IcvOiXM8PLaRUunA3BHx6WsbY8CEDf430c+H2ZYImNGiZ2IyBhq2/8E/73xa8n110XhXzo/yR4PBkpcsGQaZ540Oazw8t66kyZRXV7MEWqZd2x+cvvP9/yJtq62fkqKFAYldiIiYyXawmd/+7fsTdz7tTbuHPSP83BjNQDV5cXc8NpTw4ww75UURXjZ4iBJvvf4Fczt6gLgODF+ueGrYYYmMiaU2ImIjJE/3PF2flYaT66vLH4tv9tZk1z/rytWMKO2PIzQCsorXhw0x27zeZzVPCm5/TvP/JhoLBpWWCJjItTEzszOMLObzWyzmbWa2W4zu8XMFg3iNd5kZg8nyh81s3vMbO1oxi0iMlibHriJT7Sm7oKwuugkfvP0uuT6R15+CpesmBlGaAXnwiXTqSorBuCvh9/I1O5ghOwhuvjZIzeFGZrIqAu7xu6jwBXAH4EPAd8Azgc2mNnSgQqb2fXA94FNifKfBrYBmq5dRHLG9l338r5nvkk0EnzlzvEy7n3mLUBwu7DXrJzFBy7I+vesDKCitIhXnRr8GXgmvohLoqmJnr+65Sc0th8LKzSRURd2YvcFYL67f9Ddv+Xu1wPnAiUESV+fzOws4GPAG9z9HYnyX0ws/2r0QxcRGdjOI8/yrrvfz9Gi4Ou2Jg4Ne9+Dx4K7I6yeV8d/v24FZron7Ei6YnVqHsA/H3ozsxK1dk0W56Y/fySssERGXaiJnbvf7+6dvbZtBZ4GBqqx+xDwiLv/0swiZlY1WnGKiAzFc8ee422/fROHzAGoiMc5ufkqDrcEtUmz6yr4xlvX6A4To2DtgkmcNHUCAE93zOJqX57c9/OGR7l32+/CCk1kVIVdY/cCFvxsnQ40DHDohcAjZnYD0AQ0m9lOM7umn9du7O8B1I7YiYjIiMuna3hTwybe/ttrOJL47VoRj3NV7HzW718JQGlRhK9cs5opVWVhhlmwIhHj7156UnL9+/XXcH40llz/2H0fY8/xPWGEJjKqci6xA64BZgO39XWAmU0EJgNXAdcSNNteDewBfmRmrx2DOEVEMnrs4GO88w9v53isA4CqeJx/6z6J/3n+4uQxn7h0KafNrQsrxHHh8lWzmVYdJM47Wko4c8I7mdbTJEuMa3/9Rva17AszRJERl1OJnZktAW4G7gN+2M+hPc2uk4FXu/tX3P0W4CJgL/CvmQq5e11/D4KaPxHJUflwDd+5807efee7aE0kdbWxGDc3l/LPO99Gz2CJS1bM5C1nzu/7RWRElBUXcd0rUnfx+NTmJfxr2UrK4sGUMwe6m3nTHa/j4QMPhxWiyIjLmcTOzGYAvwWOAVe6e7yfw9sTzzvc/aGeje4eBX4GnKY+dyIyltydbzz5Na675zo64sGkuFO6Y3y7Kc6/Hv8njsWC+emWzKjmMxosMWauPH0up88P7rvbFXM+tuPNfCZaTYkH/R6Pdrfyzjuv5foHr+dI+5EwQxUZETmR2JlZLfB7gv4xr3T3+gGKHAWiwMEM+w4S/CzOmb42IlLYjnce5x//cC1feuLm5Lb5XV18pznCDd2fZFNLMAlxXWUJ33zrGiYk5liT0ReJGJ+78jRqK0oA2NcKnzrwEW48ZkyKBX3uHLj12Vu5+Od/w42P3cixDk2HIvnLPPGrJbQAzMqBO4HTgQvd/cEsyz0IzHL3eb22fxl4N1Dt7u0ZC/f9mo21tbW1jY2NgykmUlDq6upoampqSjRt5pWxvobjHU3838M38bkdv+QQqY75a9s7+Lf4HN7T9D6eagqmNSmKGD98x1rOWjRlTGKTEz20/Qhv/c7DRLuDxqDpHOWLE/6bH07p4N7KihOOLbES1kw9n7OmX8KCCcuJxaErFqcr5pSXFDF/ciXzJlXm7GjmfL6GZfhCTezMrAj4BfAq4DXunnH8uZnNAyrdfUvatuuAzwGvcPe7EttqgOeBLe7+0iHEo8ROxr18/qMw4tdwZysc2xk8mvZBSz0013Pw+F7+2LGPn1or20pLTijy+uY21tglfHjvhXR5UDNXFDFueuNKLjtt1sjEJUOyaV8T7/nxY+w5Gvzmr6aN/yz5OpXVm7i5rpbnykpfUKa4u5RI2zza2hbR1T2FeFcNHq/AvJS5dXUsnT6ZF8+cxNKZNSyZWc3M2gqKIuE2s+fzNSzDF3ZidxPBfHS/5oWjYFvc/fbEceuB89zd0spWAo8RjKC9kaBv3rXAYuBl7v7XIcSjxE7GvXz+o9DvNewO9Rsh3gWxLoh1QqyTg+0NNLQdItZ6hFjbYbrbYVhj3AAADMtJREFUGoi1HCLWfIDOaBPNkQjHIxHqi4vYWVLCtpJi9paUvODlJ8ViXHZkGr85+mZ2eurWYDXlxdx01UouWDJ9NE9dstTeGeOHD+7kBw/sYu+xdsB5deQB/qnkJzwzoZ3v1NWwqWxwU9CYG5F4EeYlROLFlFBCiZVSYqWUWjlFRVVYyUTKiksoKyqhvLiU8pJSyopKKS0qpqy4lLKiEsqKSyiJlFAcKaE0UsKMmmomVVZQWlRKSaQEMyNCBDOjtqyWudVzM8aTz9ewDF/Yid164Lw+du9y9wXpx6UndontM4DPApcAFQSJ3sfc/d4hxhMHrLZW3fNk/GpqagJwd8+JPriD0f817DR++IVbvzCxju/W1Qz53yyPw8mNM9h5+PXsj885YV/Xrg203/c94s0DTcspYYhMnkfx1BdRNGUBlTV1vH7iVq6s2UhFVT2/qJnAHydUcrQoN5tbWx9r5egPjmbcl8/XsAxf6H3scomZdRMMKDkedix96PlrFfqUDiOoEM8J8vu8aoC4u+ddD/88uIb7k8+fmUwK6Xzy7Vzy9hqW4VNil0cSM+tTSNXrhXhOULjnJaOn0D4zhXQ+hXQuUvhUTSsiIiJSIJTYiYiIiBQIJXYiIiIiBUKJnYiIiEiBUGInIiIiUiCU2ImIiIgUCCV2IiIiIgVC89iJiIiIFAjV2ImIiIgUCCV2IiIiIgVCiZ2IiIhIgVBiJyIiIlIglNiFzMxmmtl/mdndZtZsZm5m52dZ1szsOjN71syiZrbLzD5tZiWjHPZAcZ1hZjeb2WYzazWz3WZ2i5ktyrL8bDO7zcwazey4md1uZi8a7biziGvI52Vma83sK2b2mJl1mplGLRWYYX4+diau/UyPrb2O7eu4d4/guawxs18mvlPazazezP5gZmdlWT7ra9jMrjWzZ8ysw8yeM7P3jdR5pP0bQz4fM7vCzG41sx1m1mZmW8zss2ZWm+HYUX9vRAZSHHYAwmLgo8DzwFNAVl+cCV8APgz8FLgRWAZ8DJgLvGNkwxyUjwJnJ+J6CpgBvB/YYGZr3f2ZvgqaWRVwN1AN/AfQDfw9sN7MVrr7sdEOvh9DPi/gVcC7EuW2AUtGOVYZe8P5fHwYqOq1bT5wPXBnhuP/D/hRr20PDSXoPiwk+PvwTeAAUAdcA9xrZhe7+119FRzMNWxmfwd8jeD/7AvAucCXzazc3T+fC+cDfAPYD/wQ2A2cCnwQuNjM1rh7R6/jR/u9Eemfu+sR4oPgy29yYvlywIHzsyg3m+AL8zu9tr8/8RorQzyns4DSXttOBjqA7w1Q9p+AOLAqbduSxLn+W8jv1XDOazpQkVi+Kbj0wv/86ZEbn48+Xu8TiWv5rF7bHbgphPOrBOqB3wxwXFbXMFABNAC39yr/I+A4UJsj53N+hm1vTbwPb8uF90YPPdIfaooNmbs3u/uRIRRdBxQBt/Ta3rP+xmEFNgzufr+7d/bathV4Glg6QPHXAw+6+4a0sluAPwFvGOlYB2M45+XuB929fTTjk3AN83OfyZuAHe5+f6adZlZhZuVDeN0hcfc24DBBbVd/sr2GXwZMBr7Sq/zNBD94Lx5uzP3J9nzcfX2Gzb9MPGd8X8f6vRFJp8Quf5UlnnsnC22J59VjGMuAzMwIaq0a+jkmAqwAHs2w+2HgFDOrHJ0Ihyab85Lxa6ifDzNbRZA0/KSPQ94JtALtZvaUmb12WIH2HUe1mU0xs8VmdgOwnCBB6+v4wVzDqxLPvY99jESN37CCzxzfoM6nHzMSz5ne1zF5b0T6osQufz2beD671/ZzE8+zxjCWbFxD0Hx8Wz/HTCJIWA9k2HcAMGDmyIc2LNmcl4xfQ/18XJN4/nGGffcT9KV9DfA+gmvmF2Z29VCD7Md3CWq1tgDXEfSHu6Gf4wdzDc8Eou5+NP2gRK3nEUbnO2yw59OXjwIx4Be9to/leyOSkQZP5Cl3f9zMHgI+Zmb1wHqCX/hfBboI+q/kBDNbQtC8ch9BB+S+9MQczbCvo9cxoRvEeck4NNTPR6LW6ypgg2cYcOHuZ/c6/vvAJuAzZnaLu4/kiOtPA18H5gBvIUhUSsh8jcLgruEKoDPDcT3Hjsa1PtjzeQEzexNwLfCf7r4tfd8YvzciGanGLr+9jmD03XeBHcCvCWoGNgAtIcaVZGYzgN8Cx4Ar3T3ez+E9zcplGfaV9zomVIM8Lxlnhvn5OI+gli9Tbd0LuHsrQc3THIJR9iPG3Te6+13u/l3glcDpwPf6KTKYa7i9j+N6jh3xa30I53MCMzsX+DbBe/svWfx7o/beiPRFNXZ5zN33AeeY2ckEfT62unu9me0H/hpudJCY5+n3QC1wtrvXD1DkKMEv50zNrTMJRpxlauIZU0M4LxlHRuDzcQ1BH7P/HUSZPYnnSYP8t7Lm7l1mdgfwCTOr6GMw0GCu4QNAqZlNSm+ONbNSgkEV+0f0BHrJ8nySzOw04FcEP6bf6O6xLP+pUX9vRNKpxq4AuPtWd/9LIqlbRvAFOpQOwSMmMSLs18ApwKXu/uwARUjUamwE1mTYvY4gcW3LsG/MDOW8ZPwY7ufDzMoIauLXu/tgEpuTEs+HB/PvDUEFQT+56kw7B3kNP5F47n3sGoK/TU8w+vo9nx5mthD4A3AIuCRRE5etsXpvRAAldnnDzBYmvlz6OyYCfIbgy6ev0XSjzsyKgFuBlxA0Qz3Yx3HzEv2Q0v0MODMxKrDnuMXABQSTmIZmmOclBW6EPh+vIph+I2MzrJlNybBtMvBegqlRtr6w1OCZ2dQM22qAK4E97n4osW041/CfCWr43tur/HsIupL8frjnkfbvD/l8Es3qdxLUor7S3TOOcB6r90ZkIKa+nOEzs08kFpcSzF31HYI+c43u/uXEMTsB3H1BWrmbCZrTnwBKE2VXAZe7++/GKPwXMLObgA+R6vOXrsXdb08ctx44z90trWw1QR/BCcDnCSY1/QjBr+qVQ5zzb0QM87zmE3TWhuCP90tI9dF50t1/PYqhyxgYzucj7TV+BlwKTHf3pgz7P0Uw4vI3BHdBmA38P2AawXX/mxE6lz8TDGC4n2AS37nA2wn6il3l7rf1dS6DuYbN7L0EA0x+SpA8nUsw+e9H3f0zI3EuI3A+TwCnEfxo3tjrpbe5+wOJ4z7FGLw3IgMa7RmQ9Rj4QdDvJNNjZ9oxO9PXE9veATxJMGdSI/A7YF0OnM/6LM9pPRnuwEDwZftToAloJujXclI+nxdwfj9lvxf2uekR7ucjsb2GYMDAz/v5N14B3EWQnHQS1Hj9lqAv30ieyzsScR4iGGV/mCBhPS/TOWcon/U1THCrvS0EffOeBz44Cu/NkM+nn/f0hGt3rN4bPfQY6KEaOxEREZECoT52IiIiIgVCiZ2IiIhIgVBiJyIiIlIglNiJiIiIFAgldiIiIiIFQomdiIiISIFQYiciIiJSIJTYiYiIiBQIJXYiIiIiBUKJnYTCzD5lZm5mF5rZ/5rZQTNrM7OHzeylYccnIv3TNSySm3RLMQmFmd0BXEJwP8X7Ce5zOw/4MMHNwhe5+4HwIhSR/ugaFslNxWEHIOPWSqAI+C93/0LPRjN7Hvgu8AbgiyHFJiID0zUskoPUFCtjzswmEvyyvy/9D0LCnxLPC8Y0KBHJmq5hkdylxE7CsCrx/K0M+3o+ky0AZlZsZl80s6Nm1mhm3zaz8jGJUkT6Mphr+A1mdp+ZtZjZzrEITmQ8U2InYViZeH40w751iecNieePAS8DTgVOBpYBnxnV6ERkIIO5ho8BXwY+PtpBiYgSOwlHzx+F7gz7PkLQGfvOxPo7gRvcfZ+7HwY+BbzNzIpGPUoR6UvW17C73+XutwC7xig2kXFNiZ2EoeePwnnpG83sWoJf+9e7e4uZ1QFzgSfSDnscqEb9d0TClNU1POZRiYhGxcrYMrMyYClBM82NZjYf2AmcD1wN3ArclDi8OvHcmPYSjb32icgYGuQ1LCJjTImdjLXlBJ+7LwB1wHXALGAb8PfAlzw1uWJz4rkWqE8s1/XaJyJjazDXsIiMMSV2MtZ6mnCecvenCDpVZ+TujWa2J1Hm2cTmVQRJ3c7RDFJE+pT1NSwiY0997GSsrQK6gC1ZHv8t4J/NbJaZTSUYPPE9d4+NUnwi0r9BXcNmVpSYoqgkWLXyRHOuiIwC1djJWFsJPOvunVkefwMwBXia4IfIz4CPjlJsIjKwwV7DbyG4E0WPdoIRsgtGOC4RQfeKlTFkZgY0Ab9x9zeFHY+IDI6uYZHcp8ROREREpECoj52IiIhIgVBiJyIiIlIglNiJiIiIFAgldiIiIiIFQomdiIiISIFQYiciIiJSIJTYiYiIiBQIJXYiIiIiBeL/A/fGSrjQh2+QAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from pymcmcstat import mcmcplot as mcp\n", "from pymcmcstat.chain import ChainProcessing as CP\n", "\n", "pres = CP.load_parallel_simulation_results(savedir)\n", "combined_chains, index = CP.generate_combined_chain_with_index(pres)\n", "\n", "settings = dict(\n", " pairgrid=dict(height=3.75, hue='index', despine=False),\n", " ld_type=sns.kdeplot,\n", " ld=dict(n_levels=5, shade=True, shade_lowest=False),\n", " md=dict(lw=3))\n", "sns.set_context('talk')\n", "fpg = mcp.plot_paired_density_matrix(\n", " chains=combined_chains,\n", " settings=settings,\n", " index=index)\n", "tmp = fpg.add_legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Qualitatively, it appears as though the chains have converged to about the same distribution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Gelman-Rubin Diagnostic\n", "In the previous section we generated the list of chains which are required for the Gelman-Rubin diagnostic. Note, this implementation assumes that all chains have the same dimension. The output will be a dictionary where the keys are the parameter names. Each parameter references another dictionary which contains\n", "- `R`, Potential Scale Reduction Factor (PSRF)\n", "- `B`, Between Sequence Variance\n", "- `W`, Within Sequence Variance\n", "- `V`, Mixture-of-Sequences Variance\n", "- `neff`, Effective Number of Samples" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parameter: $p_{0}$\n", " R = 0.9999050250834056\n", " B = 0.0005446504072835234\n", " W = 0.0013828493830858575\n", " V = 0.0013825867235498412\n", " neff = 7500\n", "Parameter: $p_{1}$\n", " R = 0.9998927399338449\n", " B = 0.0030520757695050565\n", " W = 0.008775468759398558\n", " V = 0.008773586345638533\n", " neff = 7500\n" ] } ], "source": [ "from pymcmcstat.chain import ChainStatistics as CS\n", "chains = CP.generate_chain_list(pres)\n", "psrf = CS.gelman_rubin(chains=chains, display=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In general, `R` closer to 1 indicates your chains have converged." ] } ], "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.8" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": true, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }