nnsubspace_mnist.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This notebook compute the active subspace of the Neural Network trained on MNIST dataset.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\lib\\site-packages\\h5py\\__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "from __future__ import absolute_import\n",
    "from __future__ import division\n",
    "from __future__ import print_function\n",
    "\n",
    "import numpy as np\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "from nnsubspace.nndataset.dataset import Dataset\n",
    "from nnsubspace.nnmodel.model import NNModel\n",
    "from nnsubspace.nnsubspace.subspace import NNSubspace\n",
    "import nnsubspace.visual.subspaceplot as subspaceplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       ".output_png {\n",
       "    display: table-cell;\n",
       "    text-align: center;\n",
       "    vertical-align: middle;\n",
       "}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# center picture\n",
    "from IPython.core.display import HTML\n",
    "HTML(\"\"\"\n",
    "<style>\n",
    ".output_png {\n",
    "    display: table-cell;\n",
    "    text-align: center;\n",
    "    vertical-align: middle;\n",
    "}\n",
    "</style>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_train shape: (60000, 28, 28, 1)\n",
      "60000 train samples\n",
      "10000 test samples\n"
     ]
    }
   ],
   "source": [
    "# The package supports mnist/cifar10/cifar100/ImageNet, you can also add other datasets to the package\n",
    "dataset_name = 'mnist'\n",
    "\n",
    "dataset = Dataset(dataset_name=dataset_name)\n",
    "model = NNModel(dataset_name=dataset_name, model_id='0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Test single sample**\n",
    "\n",
    "* sample number of 259 from the test dataset is presented in the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample 259 | model prediction | ground truth\n",
      "label: 6, score: 0.69\n",
      "label: 6, score: 1.00\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKUAAAClCAYAAAA9Kz3aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAABdtJREFUeJzt3W1InWUcx/H76Jyz6dbIySo2rGS1iGFzlXW2ni0qAlm6WG3Bii1alLIa+KKX0SiC1sOq2ZYQQq2FWDDGHAY97AHTMptjiT2Qklarqa05H46n9/4v4fZ4POd3zv39vPxx3+dcwo8L/lzH+w5Fo1EPUJKR7AUAk1FKyKGUkEMpIYdSQg6lhBxKCTmUEnIoJeTMmc7FZRmVHP8gZkcmDoT8XMdOCTmUEnIoJeRQSsihlJBDKSGHUkIOpYQcSgk5lBJyKCXkUErIoZSQQykhZ1o/XQuCaLjYZDvra01Wkj3XZHdtesJkc5rb4rOwAGGnhBxKCTmUEnIoJeRQSsgJ7PT96/6Vzrwt/K7JskNZJts3uMRe1zNgskgMaws6dkrIoZSQQykhh1JCTmAHnTdLPnTmrqGmbmipyRor1pos0vXjzBcGdkrooZSQQykhh1JCTiAGne766022Zl6L89p9g1eY7LPKNSaLdDLUzBZ2SsihlJBDKSGHUkJOIAadj8N7TJYVynRe27ClzGShzva4rwlTY6eEHEoJOZQScigl5KTdoNP33C0mW5FlT29Kv93gvH9xyymT+X2jVWaRPQ36/b5Lfd172aE+Zx7p/sXnt6cPdkrIoZSQQykhh1JCTtoNOmN5NnOd3pxvy3feHx3rMtnQI6Umu+bZTpNV5B8y2T05/zm/Z7KmqvnOvGf0EpN9d26ZyVr32kcYFhw7a7KJk6d9rSeZ2Ckhh1JCDqWEHEoJOaFo1O95heeVZVT6vzhJHj3da7INeX+YrHztQ877z4TtCczBl1412cKMeTGsLrF2D1xlsqb1N5ksUf9vdGTiQMjPdeyUkEMpIYdSQg6lhBxKCTkpfczYX21/O1mRu8txpT1mPFWz2PmZO2/7yGSuSXtb760m+/5t+xz1vJ5R5/fMxM/r7d+z5+46kz198U8me73K/mPc8q3xWVe8sFNCDqWEHEoJOZQSclJ60BkusKeeUz35YrKuB+xLnKZyfMR+Zu8W+xz0RR3HfX/mTCz/3GYvlm822R273zFZ+/1vmOzOrdtNll+bmL/FhZ0Scigl5FBKyKGUkJPSg85saBmxP/l7Ybs98sjpcD8zPVlyT/5lMteAdnP2XJP9s3rcZPm18VlXLNgpIYdSQg6lhBxKCTmBHXTaRtz5Y59uM1lR44lZXs3MuR4Z6BrQmh2nPIfvtT/3e8YLx2dhMWCnhBxKCTmUEnIoJeQEdtDZ8bwdaDzP84oa9Icav3I7zyR7CTFhp4QcSgk5lBJyKCXkBHbQWdBuHw/oeZ5nf8SVusaWLEz2EmLCTgk5lBJyKCXkUErICeygk25C2dkmG6w55+veB088ZbJCr2PGa4oVOyXkUErIoZSQQykhh1JCTkpP34UHh012dtMFky1yPLO8+3H7ZjHP87yit+xnjve7jySVRG681mRHi/ea7Ldx+/cVvjYrS4oZOyXkUErIoZSQQykhJ+1erdzfuMJkrTfU+76/bsg+y3xXfbnJlr3SarLoWPxf5BTKso/u84qvNtG6D5pNtnlBj8mu+9o+G73w4cQcKfJqZaQsSgk5lBJyKCXkpN2gE7l9lcn+rLanGIdXvee8Pz8zx9f3lHyz0WT//j3f170uBV9kOfPhdQMma3MMbh2jEZNtrKs22dLm8yYLHW33s8QZY9BByqKUkEMpIYdSQk7aDTp+RcPFzrx/h30Y+ler3zfZRSHHSUuCfHnBfndV7ZMmu/zlY4lYjm8MOkhZlBJyKCXkUErICeygMy2lK03UVzNmMtdJi19Nw+7ToOoG+1OzKz9xPPmi5YeYvztRGHSQsigl5FBKyKGUkMOgg4Rh0EHKopSQQykhh1JCDqWEHEoJOZQScigl5FBKyKGUkEMpIYdSQg6lhBxKCTmUEnIoJeRQSsihlJBDKSGHUkIOpYQcSgk5lBJyKCXkUErImdYTMoBEYKeEHEoJOZQScigl5FBKyKGUkEMpIYdSQg6lhBxKCTn/A8dEGQA45s4hAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 144x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEECAYAAABdvxc1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xlc1HX+wPHXh2GA4ZDDGxQwr0ypPCp33cpOtlw3O93N1qu081fWdli5iVupZcfWbttuh2luh2Utm1LSlrWlpeZ9ZFopoCAiyeA1wDDz+f0xDHHMAANzMbyfjweP4Pv98v2+ofow7/l8Pu+30lojhBBCCCGEECLwwgIdgBBCCCGEEEIIB0nQhBBCCCGEECJISIImhBBCCCGEEEFCEjQhhBBCCCGECBKSoAkhhBBCCCFEkJAETQghhBBCCCGChCRoQgghhBBCCBEkJEETQgghWkEplaWU0g0+ihtcc5tSap9SqkIptVEpdW6g4hVCCNE+SIImhBBCtN5uoGedjwznCaXUeOA5YC4wFPgK+EgplRqAOIUQQrQTkqAJIYQIOUqps5VS/1VKHXYxy9XXi4+q1loX1/k4XOfcPcAirfXLWutdWuv/Aw4Ct3rx+UIIIUKMJGhCCCFCilJqCPA5sAsYDVwIFAPrgRuAvXWufUgpdbyZj6aWJZ6ilCqsWcb4tlLqlJr7RgDDgY8bXP8x8Esv/ahCCCFCkNJaBzoGIYQQwmuUUp8CZq311XWOzQOu0Vr3b3BtEpDUzC0LtdYWF8+5DIgDvgO6AbOAU4HBQCRQCJyvtf6izvc8AkzQWg9szc8mhBAi9IUHOgAhhBDCW5RSXYDzgcwGp04Ajd6R1FofAY605lla648aPHstjtm5ScDbzssahugqDiGEEMJJljgKIYQIJcMBA7C1wfERwDcNL/bCEsdaWuvjwE6gP1AK2IAeDS7rBhzy8GcSQgjRgcgMmhBCiFBiqPmnyXlAKdUPx4zalS6u/wfwTjP3LGzJg5VSUTiWOH6mta5SSm0ELgHerXPZJcB7LbmfEEKIjkkSNCGEEKFkHXASeFIp9WcgFXgeeFtrvbLhxW1Z4qiUegpYDhTgmBn7ExADLK655BlgiVJqPbAGuAVIxpEUCiGEEC5JgiaEECJkaK1/UkpdiyM52opj9utVYL4PHtcLeAvoAhwG1gIjtdb5NbEsVUp1xlE8pCewA7jceV4IIYRwRao4CiGEEEIIIUSQkCIhQgghhBBCCBEkJEETQgghhBBCiCAhCZoQQgghhBBCBAlJ0IQQQgghhBAiSEiCJoQQQgghhBBBIqTL7Hfp0kWnp6cHOgwhRI2NGzeWaq27BjoOf5ExSIjgImOQECKQWjoGhXSClp6ezoYNGwIdhhCihlKqQ/V/kjFIiOAiY5AQIpBaOgbJEkchhBBCCCGECBKSoAkhhBBCCCFEkJAETQghhBBCCCGCREjvQWuJ7M2FLMjdTZHZQnKCifsyBzJuaEqgwxJCCCGEEEJ0QB06QcveXMiD72/HYrUBUGi28OD72wEkSRNCCCFaSd78FEKI1uvQSxwX5O6uTc6cLFYbC3J3BygiIYQQon1zvvlZaLag+fnNz+zNhYEOTQgh2oWQTNCUUmOVUi+Vl5c3eV2R2eLRcSGEEEI0zV9vfhYcLfDq/ULR4ZOHAx2CEKIVQjJB01ov11pPj4+Pb/K65ASTR8eFEEII0bRCP7z5ubpwNVf85wpe3PoiWmuv3TfUlFhKAh2CEKIVQjJBa6n7MgdiMhrqHTMZDdyXOTBAEQkhhBDtV/bmQpSbc9568/Ob4m+Y8dkMqu3VHDv5E6x8EKpOeuXeQggRDDp0kRDnhuU/Ze/gWGU1sZHhPDZuiGxkFkIIIVphzvKduJrPUuCVNz+3lGzh9k9vp9JWydWnjOW+HatQhRuhohyufLHN9xdCiGDQoWfQwJGk/fX6oQAM7BEnyZkQQgjhoezNhQz988eUnbS6PK9pe3Xkb3/6lts+uQ1LtYXfpF3Kn3Z95UjOElJh9Mw23VsIIYJJh55BcxqamohSsP1AORVWG1ENlj0KIYQQwrWGLWtcSWnj8sYfyn7g5v/ezDHrMS7pdT6P7tmAoWgLJKTB5BxI6N2m+wshRDCRBA2INxkZ0C2O3YeOsbOonOFpSYEOSQghhGgXXFVtbKityxtNRhMxxhhOTxrEEz9sJfzgNk5E9+YGyyy2zN9GcsL30mtNBNyo+avcFsnxhZQEE2tmXujR90yePJnFixdz44038sorr9Q7d//997NgwQLGjBnDihUrADh27BhPPvkk7733Hvv27aNTp04MGjSIm2++mfHjxxMW1uEX4/mEJGg1hqUlsvvQMTbml0mCJoQQQrRQc9UZFXD30i0syN1dm0R52sg6JTaFJaP/Sqd3J2Ms3sHxmN785uiD5FljgJ97rUHbl1IK0VqFZgt588f47XnpM3Na9X29e/dm6dKlPPfcc8TEOP4fqq6uZsmSJaSmptZeZzab+dWvfkVZWRmPPfYYZ599NhEREaxevZpHH32UX/ziF6Snp3vjRxENSNpbY3haIgAb8soCHIkQQgjRfjRXnVHXfBSaLdy9dAvpM3O4e+mWZhtZF58o5s1dbzq+OPETXd+ZQmTxDkjqy4Tq2eRZE+pd74tea0KEotNPP53+/fvzzjvv1B7LyckhKiqK0aNH1x576KGH2LdvH+vWrWPKlCkMHjyY/v37M2XKFDZt2kSPHj0CEH3HIAlajRE1CdqmgjLpqSKEEEI0kL25kFHzV9FnZg6j5q+qTahctaxxRzf4p1PD5KrUUsq0j6cxb/083tm+EBaPhUPboXM/mJzDtvJol/f3Zq81IULZjTfeyMKFC2u/XrhwIVOmTEEpR6MMu93O22+/zYQJE+jVq1ej74+KiiIqKspv8XY0kqDVSOscTeeYCEqPV1FwRPqpCCGEEE7OQiDuZr0iw9v+csKZXJVVlDHt42nkHc3j1Ph+ZH75EpTshM79HQVBOvV0O2vnrV5rQoS666+/ng0bNvD9999TXFzMypUrmTx5cu350tJSysrKGDRoUOCC7MBkD1oNpRTD0xL5+NtDbMgrI61zTKBDEkIIIYKCq0IgFquNOct3UmG1N1skpCWSE0wcrTrKzf+9mR/MP3BKXCr/LCwi/vB30GUgTFoOcd0Bx6xdw8qRJqPBK73WhOgIEhMTufLKK1m4cCEJCQmMHj263v4zWU0WWDKDVodzH9rGAtmHJoQQQji5WzpYdtLqleQM4M6Le3PrJ7ey68guUmNSeOVgCUmHv4Oup8LkFbXJGTgKgcy7KoOUBBMKRzW7eVdlSIEQITwwdepUXn/9dRYuXMjUqVPrnevatSuJiYns2rUrQNF1bDKDVkdtgiaFQoQQQohayQkmn5YPVwpmr3kUY8I2EsI788qhUroe/h66DnLMnMV2bfQ944amSEImRBtcdNFFREREUFpayrhx4+qdCwsLY/z48bz++us88sgjjfahVVRUAMg+NB+RGbQ6hqTEE2EIY0/JMcot1kCHI4QQQgSF+zIHonx4f62h8nAmhhO9WbC3iJ6Hv+c7e2/GHH2A7O+rfPjk0Ddq/qpAhyCClFKKbdu2sW/fPiIjIxudnzt3LqmpqZxzzjm89tpr7Ny5kx9++IElS5YwfPhwiouLAxB1x9BuZtCUUr8AngGqgCJgotbaq1lUlNHAkJRObCows2W/mfMHNH7HTgghhOhoxg1NYcbSLV67n0Ep7FoTpjQ27Uj9etsqefPQbnqpUrbZ+zCxaibmqgjpb9ZG/myc3NGlJJha3Zustc9rq7i4OLfnEhMTWbt2LU8++SRPPPEEeXl5dOrUidNOO40//elP9fasCe9qNwkakA9cqLW2KKUeB64Alnn7IcPTEtlUYGZj3hFJ0IQQQggcVRwNSmHzQuEAo0ERHqawWK0YU97AUJFMr58G8WbEPHqoMjba+zOl6n6O4ijW5SzBLwma5yI0xPSdD/iveXJHtmbmhYEOoVmLFi3y6Hx8fDyPP/44jz/+uO+CEo20myWOWusirbXzbaBqwO6L5wxPSwKkUIgQQggBP5fY90ZyFm0Mw2bTWKxWolLewhi3i6ik1fzd9Dg9VBlr7YOYWDWzNjlzkv5mrdO/qoqwCHOgwxBCeMjvCZpS6g6l1AalVKVSalGDc0lKqX8rpU4opfKVUte7+P4+wGXACl/ENywtAYAtBWaqbT7JAYUQQoh2w1WJ/daqrNbYsTmSs047CbNF8NLBEk61l/OFLYPJVfdzgsbLtqS/mRCiIwnEDFoR8Biw0MW5F3DsMesOTABeVEoNdp5USnUCFgN/0Fr7ZNdwt7goUpOiOVFl47viY754hBAiSLTkTaE61w5TSn2hlDqulDqklLrLn7EK0RbZmwsZNX8VfWbmMGr+qtoG0y3hzdkrm64mKuVtjJ12EGaL4OWDhxhpLee/tmFMs/6RChoXKgjl/mYyBgkhXPH7HjSt9fsASqkRQG3NTqVUDHA1MERrfRxYrZT6APgDMFMpFQ68BWRprXf7MsYRaYkUHDnJpoIyhqTE+/JRQojAqvum0JlAjlJqq9Z6Z92LlFJdgJXA3Tj2vkZQZ/wSIpg5lyg6Z8EKzRaPCm9ERxg4UeWNGTQbUcnvYOy0nTCbkZeLiznbepwPbWdzl/UOrHVekqQkmCgyW0hOMHFf5sB6cWZvLmRB7m6359sZ345BhghvxyuE8INgKhIyALBprffUObYVOL/m898D5wCPKKUeAV7UWi9teBOl1HRgOtDq6jLD0hJ5f3MhG/PLmPiL9FbdQwgR3Jp7U6jB5fcAuVrrN2q+rgSke6doF1wtUXQW3nCed5XsZG8uZM7ynV5KzkAZLBhMBwizhTuSs6oT/Ns2inutt2DDUHtdSoKpttiCMxm7e+kWkhNMXHBqV97bWNjqZDOY+GUM6j4YqGTU/FXtooCFEMIhmIqExALlDY6VA3EAWuslWusuWuvRNR+NkrOa617SWo/QWo/o2rV1VRidDas3SMNqIUKZuzeFBru4diRwRCn1lVKqRCm1XCkl9YVFu+BuiaIzuSk0W9B1vs7eXFg761Z20nvdbLQtlqH7f8HCmuTs7erR/NF6a73kDKhdzuiMoW58b6wtaDLZbGf8NgZJqX0h2pdgStCOA50aHOsE+H0j2IDuccRFhlNotlBcXuHvxwsh/KPJN4Ua6AVMAu4CUoF9OJZcN6KUml5TCGnD4cOHvRiuEK3TVIENV8nO3e9sYc7ynV4qDGIjPG4roMkMW89C9U+GV51kcfUlPFh9E3YXL0OcM2GuZv7c1ZFsp1UeZQwSQrgUTAnaHiBcKdW/zrEzgJ1urvcZQ5jizFRHNcdNUm5fiFDlyZtCFuDfWutvtNYVwBzgl0qpRptUvTGLL4Q33Zc5EJPR0PyFNbTGSzNnNqJSlmLq9Rand1nEC8bniVA2/lk9htnVk9EuXoIkRhtrP/ck6WqnVR5lDBJCuBSIMvvhSqkowAAYlFJRSqlwrfUJ4H3gz0qpGKXUKBzNqJe04hljlVIvlZc3fGOq5UbU9EOTZY5ChCxP3hTaRv03752fKx/FJoTXjBuawryrMjAof/7n6iylv41wm4E5VV8Qruw8Xz2OedXX4+5/nYo6M2bukq6G39mOqzzKGCSEcCkQM2izcLwTNBO4oebzWTXnbgNMQAmOqftbG1Yyagmt9XKt9fT4+NZXYHTuQ5OG1UKEJg/fFHoNuFIpdaZSygj8CVittZYOsKJdGDc0BbsXGk23TDVRKW9i7LQDo83A4kOFDK2qZJ719zxTfR1N5RQW68/9R13N/JmMBiaMTCUlwYTCUVBk3lUZ7a5ACPhvDEq2VhM3aCaZyzK9Gb4QwocCUWY/C8hyc+4IMM6f8bhzZmoCYQp2FpZTYbUR5cHyECFEu3Ebjp6MJcBP1LwppJQ6F/hIax0LoLVepZR6CMgBooHVgNt+RUIEo3iTEbPFe0U/XFLVRKW8gTFuFxE2A68XFzKo0srD1VN5w3YxAKP6JrHmxyPN3qruXrQQKanvis/HoNwDRYyK+jdFPe+Sao6+9mwGlBf473nxqXD3do++ZfLkySxevBgAg8FAcnIyY8aMYe7cuSQmOiYntm7dyiOPPMK6deswm81069aNs846i2eeeYa0tDSv/xiisWAqs+81SqmxwNh+/fq1+h6xkeGc2qMT3x48ytb9Zs45pbP3AhRCBAV3bwpprb/EsYG/7rEXgRf9FJoQXuePFY6R3T7EGLeLKFsYi4sPMKDSxgzrbXxgH1V7TUuSM6dxQ1NCLSGrx19j0JqZF5KxWKo5+lx5AWS1fnuNx7Jat1Ls4osvZsmSJVRXV/Ptt98ydepUzGYzb731FocPH+aiiy4iMzOTnJwcOnfuTH5+Pjk5ORw9etTLP4BwJyQTNK31cmD5iBEjprXlPsPTEvn24FE2FpRJgiaEEKJdadjQ2Zsl892x/fQreps28PSRfE6phJutd/OpfXij6xSuKzImmIwujoo2iU91vJDvI51BhENkZCQ9evQAoFevXowfP55FixYBsGbNGsrKynjttdeIiHA0Ok9PT+f88893dzvhA8FUxTHo1O5Dk0IhQggh2hFXPcR8pxrQRFLF39Vr5BT/QO9KA5Ot97tMzsCRnDV8AWIMU2T91lULMNEmd2+vndVJSTAxav6qAAckgsnevXtZuXIlRqPjzZEePXpgt9tZtmwZ2m/7VkVDkqA1oW6hEPmPVAghRHvxwHvbvNTHrBmqClPvRcR0zeFV45NcathIuY7hhqqH+NredLJlMCgSTMbaYh8Lrj0jpJczBoM1My+UZY6ClStXEhsbi8lkom/fvnz77bc88MADAIwcOZKHHnqISZMmkZSUxKWXXsrcuXPJz88PcNQdS0gmaN4osw/QK9FEt7hIzCet7C094aXohBBCCN+Z8PLXVFbbm7+wrcIqMfVeRHjsDyQkrOZU43eU6ASuq3qELbr5PeBWmyYmMpx988ewZuaFkpz5WLINMhZn0Lff/YEORQTYeeedx5YtW1i/fj3/93//x+WXX86dd95Ze/7xxx+nuLiYl156iYyMDF599VVOO+00Pv300wBG3bGEZILmjTL7AEopWeYohBCiXfGkCEerhVmI7v0q4TF7SazWLCkuosKWxLVVj7BH927xbTxpRi3aJnfqdrZP2k6JMSRf+gkPREdH069fPzIyMnj++ec5efIkjz76aL1rOnfuzLXXXsvTTz/Nrl27SE9Pb3SN8J2QLBLiTcPTEvloRzEb88u47qyW/9ERQgghvKlh0Q9nyfm6xxOifV9kQxlOYEp9FUNUEV2tmkXFB6mq6sY1VQ9SjGcFtdw1oxa+lT4zB3AsLZWy+2L27NlcdtllTJ8+neTk5EbnIyIi6Nu3L0VFRQGIrmOSBK0Z0rBaCCFEoDmLfjj3lRWaLTz4/nY25B/hvY2Ftcd9XalRhR/FlPoKhsgSkq12Fh08yE/W3kysmskROnl0L5PRwH2ZA30UqXAn2QZFg2YCYK5KACRB6+hGjx7N4MGDeeyxx7j88st5++23+d3vfseAAQPQWrN8+XI+/PBD5syZE+hQOwxJ0JoxODmeyPAwfig5jvlkFQnREYEOSQghRAezIHd3o6IfFquNt9btx+bPIlZaEUUlKVXVLCwuZq91ADdV3csxoj26TYLJSNZvB8u+swDInbrdUXY/q5yMxRmBDie0OFsa+PN5XnLPPfcwZcoUxo8fT2xsLPfeey/79+8nPDycPn368NRTT3HXXXd57XmiaSGZoHmjUbVTRHgYZ/RKYH3eETYVlHHhqd3bHqAQQgjhAXd7tfyanAHn6x+YU/w9nahivXUYd1rvoBLP3ri8YWQqj42TxCBYpM/MkaWO3nL39kBH0Cxnv7OGrr/+eq6//noA6XkWBEJyp6i3ioQ4DXMuc8yXZY5CCCH8L5B7tcIiSojo8gnjwr7gFePTpNkryK06j9usd3mcnAEs/WY/2ZsLfRCp8FRyTDIDRzwvpfeFCDIhmaB524iaBG2DVHIUQggRAPdlDsQYpvz+3LDIIkxp/ySy6yecl7SEcGXnherfMrN6GjYMrbqn1aZZkLvby5EKj9QsxcvdsZaiE1L4QYhgE5JLHL3NOYO29YAZq82O0SB5rRBCCP+Ylb2dN9YW4N/FjBAWVUB06kKUoYJRJy1cfuIkj1pv4FXb5W2+t5TXD7C6S/EWZ5CSYKqt7AhS3VGIQJMErQWSYiI4pUsMe0tPsOvgUU7vlRDokIQQQnQAs7K386+1BX5/rsG0F1PvRShDFRedOMnjh8qYab2VbPuvvHJ/Ka8fXBomY3WTNSGE/7U6QVNKxQKdgSKttW/r+gaB4WmJ7C09wYa8MknQhBBC+Fz25sLAJGex32FK+RcqrJrLj59gVskxbrfew+f2oV65v9GgpLx+EEm20aiaY0zfBGBMYAISQni+B00p9Rul1CagHPgRyKg5/opS6novx9cqSqmxSqmXysvLvXZP6YcmhBDCX5x9z/zPTnSXXFRYNVcfPc79JRVMqnrIa8lZYrSRBdecIeX1g0ju1O1sn1TnY18BYRHmQIclRIfm0QyaUmoc8B7wKfAA8GSd0/uAScCbXouulbTWy4HlI0aMmOatezoTtE1SyVEIIYSPuep75g9dOcpfDh1kd6cyLjOH87uq2Xyve7XpnsYwxYJrJSlrbxouc5R9aUL4j6dLHGcDr2mtb1JKhVM/QdsB3Oa1yIJM366xxJuMHCyvoNBsIUXWzwshhPCC7M2FLMjdTZHZQnKCifsyB/q57LnGELub3ic68brxCVL1YXqU9eSaqpkU0rXNdw8LU4wbmuLy55SkLQjFp5JsraZo0Mx6h81VCYAkaEL4g6cJ2iDg/prPGxaUKsOxJy0khYUphqUm8Nnuw2zML5METQghRJs5lzI6Z8sKzRYefH87isZ/ZH3DTmS3D4novJpfl1lINR9mq/0UplTdzxE6eeUJldV2ZmVv572NhY1+TkCStGBz93ZyGx7LiiejTyCCEaJj8jRBOwp0cXMuHTjcpmiC3PC0REeClneE356RHOhwhBBCtHOuljL6b2mjjaie72NM2Ei41gytOs4q25ncbr0TC1FefdJb6/Zj0/VTTovVxoLc3ZKgtSNNVXeUJZAiUIYMGcI111xDVlYWAOnp6dxxxx3ce++9gQ2sDTxN0P4LPKiU+gg4VnNMK6UigTuAj7wZXLAZnpYESKEQIYQQ3uHfpYx1KCtRKW9hjPuWKLud50pKyT/2C2ZVT211A+qmNEzOnKQfWjvhZtljXUeBzGXJ5F7TaP6tw8hclunXxt/JMZ7/vidPnszixYsBMBgMJCcnM2bMGObOnUtiYk3f361beeSRR1i3bh1ms5lu3bpx1lln8cwzz5CWlub1n8PbvvnmG2JiYrx6z6ysLJYtW8aOHTu8el93PE3QHgbWA7uBD3GswJgJnA7EA+O8Gl2QOaN3PIYwxa6DxzhRWU1MpLSRE0II0bSm9l4ZlHKbvPhMWCWmXq8THvMjcTY7fz9UwmfHx/Kc7SpA+eSR7n5O6YfWTrha9giQFQ9ZjorZ6TNzmkzgOoKiE0Vsn+S/6qsN2yO01MUXX8ySJUuorq7m22+/ZerUqZjNZt566y0OHz7MRRddRGZmJjk5OXTu3Jn8/HxycnI4evSol3+Cn1mtVoxGo1fu1bVr2/fOBppHZfa11nnAMGAFcAlgA84D1gLnaK3997ZBE3xRZh8gOiKc03p2wmbXbN0vJWiFEEI0zbnHrNBsQfPz3qvszYWA+5klXzJ1/4DwmB/pUm3j1aIS3jw+iedsV+Or5MxoUPz+nN6YjPVn5kxGg/RDCyHOvfnpM3PcfoyavyrAUQqAyMhIevToQa9evbj00ksZP348H3/8MQBr1qyhrKyM1157jeHDh5Oens7555/Pk08+SUZGyxLCdevWMWzYMKKiohg6dCgffvghSik+//xzAD7//HOUUnz44YecffbZREREkJuby48//sgVV1xBjx49iImJYdiwYaxYsaLevUtKSrjiiiswmUykpaWxcOHCRs9PT0/nqaeeqv26vLyc6dOn061bN+Li4jj//PPZsGFD7flFixYRGxvLp59+ypAhQ4iJieGCCy5g3759tefnzJnDzp07UUqhlGLRokWe/Mo95vEUkNb6AHCjD2LxGl+U2XcanpbI9sJyNuaX8ct+7rbjCSGEEO73mDn3XqUkmPy6zDEGC08c+Z6cCAv3lJ5gruUuPvNSjzOXz4sw8PiVGYwbmsKItCSp4hhq4lMds2jAGiDTmux2Fi3ZBrv3zPdjcKIl9u7dy8qVK2tnr3r06IHdbmfZsmX8/ve/RynP3rg5fvw4v/nNb7jkkktYsmQJRUVFzJgxw+W1DzzwAE8//TT9+vUjLi6OoqIiLrvsMh577DFMJhNLly7lqquuYtu2bZx66qmAY4lmfn4+n3zyCdHR0dx9993k5eW5jUdrzZgxY4iPj2fFihUkJSWxePFiLrzwQnbv3k3Pnj0BqKysZN68eSxcuJCoqCgmTZrELbfcQm5uLuPHj2fHjh2sWLGiNsmMj4/36PfiKVmj56HhaYks+ipP9qEJIYRolrs9Vs7j6Z39k6Ap4xE6WxWvRSwgQ+UxsrgTU6seZJvu69PnJkRHADBq/qraxOzZ8WdKYhYq7q6/nK+p3VCtXY4nvG/lypXExsZis9moqKgA4JlnngFg5MiRPPTQQ0yaNInbb7+ds846i9GjRzNhwoQW7T974403sNlsvPrqq5hMJgYPHszDDz/MhAkTGl2blZXFpZdeWvt1165dOeOMM2q/fvjhh1m+fDnLli1j1qxZ7Nmzh48++ojVq1czatQoABYvXswpp5ziNp7PPvuMLVu2cPjwYUwmxyzvo48+yvLly1myZAn33+8oTl9dXc0LL7zAwIGOWf17772XKVOmYLfbMZlMxMbGEh4eTo8ePZr9HXiDp42qG88j1qe11kE9u9ZWdRtW2+2asDDfLAkRQgjR/iW7mSFLTjCRvbmQNT8e8XkMBtM+YnovIvPYCYaUFbHP3p3J1gfI175/oeFc0inl9QU4lkE2VQmMOFdNAAAgAElEQVSy7nVSEdJ3zjvvPF566SUsFgsvv/wyP/74I3feeWft+ccff5x77rmHVatWsXbtWl599VUef/xxPvjgAy666KIm7/3dd98xZMiQ2mQI4JxzznF57YgRI+p9feLECebMmcOKFSs4ePAgVquViooKTj/9dAB27dpFWFgYZ599du33pKWlkZzsvrL6xo0bOXnyZKN9aRUVFfz444+1X0dGRtYmZwDJyclYrVbMZjNJSUlN/sy+4OkM2oU0bs2SBMQB5pqPkJacYCI5Poqi8gp+OHycAd3jAh2SEKKVlFJJwKvApUAp8KDW+k0X12XhKJJUWefw6Vrrvf6IU7Rf92UOrJegOJ2orOaB97b5/PnhcTuITn4LHWbjhLGKjfZTmO7FHmfNMSjV5BLPjq4jjUHJNijqeRdxPZu/VipC+lZ0dDT9+vUD4Pnnn+eCCy7g0UcfrS1TD9C5c2euvfZarr32WubNm8fQoUN59NFHm03QtNYtXhbZsNLivffey8qVK3nqqafo378/0dHRTJw4kaqqqtp7e8put9O9e3e+/PLLRuc6dfp5HAwPr58SOX8Gu93u8TO9waMETWud7uq4Uuo84B9A4/nLEDQsLZGibQfZmF8mCZoQ7dsLQBXQHTgTyFFKbdVa73Rx7VKt9Q1+jU60e84kZM7ynZSdtNYeN1us7r7Fa4wJa4nqkY1WcN3RY4wsSeMPPuhx5o5Cyuu3QIcZg3KntqC6YU1VSKkI6V+zZ8/msssuY/r06S5noyIiIujbty9FRc3XAhw0aBCvv/46FouldhZt/fr1LYpj9erVTJw4kauvvhr4eZZrwIABtfe22+188803/PKXvwSgoKCgybiGDRvGoUOHCAsLa3IpZHMiIiKw2fzVo9LDKo7uaK2/AJ4F/uqN+wW7ETXLHDfkyT40IdorpVQMcDXwJ631ca31auAD4A+BjUyEmnFDU/BvsUZNRJePieqZDQruKDNzyqEzuMX6R78mZxNGptZW9mtIyuvLGCSCx+jRoxk8eDCPPfYYK1as4IYbbmDFihXs2bOH3bt389RTT/Hhhx9y5ZVXNnuvCRMmYDAYmDZtGt9++y2ffPIJc+fOBWh2Zm3AgAH8+9//ZtOmTWzfvp0bbrihdo8cwMCBA/n1r3/NzTffzNdff82WLVuYPHlyveWUDV188cWMGjWKK664go8++oh9+/bx9ddfM3v2bJezau6kp6eTn5/Ppk2bKC0tpbKysvlvagNvFgnZC/iuFFQQcTas3iSFQoQICKXURCBHa/2Ti3NJwG+01q83c5sBgE1rvafOsa3A+W6uH6uUOgIcBP6mtX6xFaGLEDcreztvrdtfb+YoJcHklxkzp8ikL4jouoowrXmk9Ajfl41lnm0sviqj31BKnQqNzjYDdZc5hkJ5fRmDfCslwYS5KqFNhUVa08S5rZJjkv1aDCU5xv3eK0/dc889TJkyhfHjxxMbG8u9997L/v37CQ8Pp0+fPjz11FPcddddzd4nNjaW5cuXc+uttzJ06FBOO+00srKyuOaaa4iKavoNomeeeYYbb7yRc889l8TERGbMmFEvQQNHyftp06Zx4YUX0qVLF2bPnk1JSYnbezrL+c+aNYtp06ZRUlJC9+7dGTVqFBMnTmzZLwe4+uqref/997noooswm8289tprTJ48ucXf7ynVmvWcjW6iVDjwGnC21jpoRt0RI0boun0OvMVqs3N61sdYrDY2zrqYzrGRXn+GEKFIKbVRaz2i+SubvY8N+IXWutG6CaXUcGC91trQ+DvrXXcu8K7WP1dKUEpNAyZorUc3uPY0HHtsDwHnAO8B92it33Jx3+nAdIDU1NTh+fn5Hv50or2alb2df60tCGgMMVh4MvIv/CulhKllJ8gun0qOfaTfnu+quENTjbr9TcagIPdsBpR74f+h+FQykvBr02jh3n/+8x+uvPJKSkpK6NKlY7eoaukY5GkVR1cdBiNwvAvUGbjFk/v5ilJqLDDWuQHS24yGMM7sncDXe39iY34Zlw72T8lNIUStpqYCYoDqFtzjODSqlNAJONbwQq31t3W+/Eop9RxwDdDoxZHW+iXgJXC8SdSCOEQ7Vjf5COi/7DAL3e3HWRTxNINUAWcXxnFz1b1s0gP8GsYFp3ZtdGzc0JRQLAgiY5Av3N36hCp9Zg5588c4vsiKh6RULwUlPOUsfd+7d2927NjBjBkzGDt2bIdPzjzh6RLHMBpXcTwGvA+8rbX+3BtBtZUvG1U7DU9LdCRoBZKgCeEPSqkzgWF1Do1VSg1pcJkJ+B3wfQtuuQcIV0r111o7rz8DcLU5vyGNv9aLiaDlavleIKjwMhJSX2J05SFO/amYH+09mWK9nwLd3e+xfPbdYb8/019kDApudUv459WspJOS/t43d+7c2j1lDZ177rl89NFHHDp0iNmzZ3Pw4EF69OjBmDFjeOKJJ/wcafvmaRXH0T6Ko92p2w9NCOEXVwCzaz7XOEpOu/IT0Gw/Rq31CaXU+8CflVI34aigdgXwy4bXKqWuAL7AscToLOBO4CFPfwDRfjlnygrNFgxKYdO69p+BFBZ1gPjer1AdXsEeBZ8zkBlVf6Sc2IDEE+LVGWUMCmL1kqxnU0m2VreoEuRRIGNxYPastUe33HIL1113nctzzmId999/f20DaNE63iwS0qEMS3UkaFsPlFNZbSMyvMml5kKItvsLsAjHu8Z7gauAzQ2uqQQO6ZZvrr0NWAiU4HhRdavWemfN3pCPtNbOV7m/q7kuEjgAPKG1XtyGn0W0Iw1nypxJWaCTM0PsLuJSllAdZuccSwUXHOzHzZU3U4UxYDGFeHVGGYPai7u30+JUq6a0vz+Le7RnSUlJAWnc3NE0m6DVVCpqsRZULQoJ8dFG+neL5fuS4+wsOlqbsAkhfENrXQ6UAyil+gAHtdZVbbznEWCci+Nfws9TEFrr37flOaJ9W5C7O+DLGBuKSPyKqO4fYFPw22PHSSo+j4erryWQq95CoTpjU2QMCl3pM3OIG9SyJZHuyFJJ4U0tmUFb5MH9NNAhEjRwLHP8vuQ4m/LLJEETwo+01u2oLJlo74Jt2Z4pfi3hPT5AA9PLjrK3ZDxv2N1VZvctheMPf0qAqzP6m4xBISQ+lTyuJ9Oa3Kbm2M6lkg3J0knRGi1J0Pr4PIp2anhaIm9/s58NeWXcdG6goxGi41BK2WlcsKie5kpcC9FSyQkmCoMkSevOEf5a8R+erazkivIqssvuYK39tIDE0tGSsrpkDAohNZUj25pC1asi6ZQVT4a8ihat0GyCJu8SuecsFLKxoAytdbMd0oUQXvNnGr846gxcimOPxiJ/ByRC132ZA7nnnS3YA7nlzHCCIbqIV43P0l2ZeexAd6ZbH+JHHZjkKDHa2NGXc8kYJITwGSkS0gZ9usSQFBPB4WOV7D9iIbVzdKBDEqJD0FpnuTqulDIAy6nZJyKEN2zIPxLQ5EwZS+mS+iKnVR6mW6mZtfZB3FI1AzNxAYvJfNIasGcHAxmDREN1y/w7rY7sQrK1utUFSJKt1eQeKIL41Db1iBPtj8cJmlIqE0dD6oFAVMPzWutTvBBXu6CUYlhqIp/sOsTGgiOSoAkRYFprm1Lq78DfcFRcE6JNsjcX8q+1BQF7fpgpj6Ter1BhqCZPh7PYfi6PV03DGuD3V0O8WmOryRjUcbmeUR7j8dLJukslMxZnQFa5o9Kk6FDCPLlYKXU58CEQDZwKfAcUAL0BO/A/bwcY7GqXOUo/NCGCRSQgNYBFmznL6weKKW4T8an/oNJQzbknLZxacAFZVbcEPDlTENLVGr1AxiAhRJt4Osr/CXgBuBuwArO01puUUgNw7K/8yMvxBb0R6Y4EbUOeJGhC+ItSKtXF4QhgCDAf2ODfiEQomrN8Z4DK62sSu+RQ3XU11cC4oxb2H5zI6/bhAYilPgVMGJnaIQuD1CVjkPCVukslY/omOGbR+qSCF/q0JVuryT0WLssl2wFPE7RTgUdwzJZp5/drrfcopbJwJHDveDPA1lBKjQXG9uvXz+fPykiJx2hQ7D50jGMVVuKiAtcgVIgOJA/XFdQU8CNwu1+jESEje3MhC3J3B7RqY2piDmVdV6O05qafrCwv/SPfa1f5gO8pIN5kpNxiJbkDV210IQ8Zg4QP1F8qOcbtde64rCZZI2NxBpQHbsm2aDlPEzQ7UK211kqpw0AqsL7mXBHQ15vBtZbWejmwfMSIEdN8/awoo4HByfFs2W9my34z5/bv6utHCiFgKo1fHFUA+cA3Wuvg6ios2oVZ2dt5Y21B07XTfezSsG+Yb3mP+y3xnGOO4x/mh/iJwO0/0UBMZDhZvx0siVl9MgYJIXzG0wRtN5Be8/kGYIZSag1QDfwRxztKHc6ItES27DezIa9MEjQh/EBrvSjQMYjQkr25MKDJWVhEMTfav+bh8HcB+HXhAB6yTqOKwK/KKDRbavfiSZLmIGOQCFauqkk6xfRNcPRla8NySWm87R+eJmhvAINqPp8NfAIcqPnaBlzvpbjaleFpibyyeh+bCmQfmhD+pJTqhGPPRwpQCOzQWh8NbFSiPVqQuztgyVlU7HZiUt7kJ8sJqg8pnqoezz9sY3GslgsOFquNBbm7JUFrQMYgEWya7k84Bp5t2zJHabztHx4laFrrF+p8vlEplQH8GkdVx0+01t96Ob52wVnJcXOBGZtdYwgLnj+qQoQqpdQjOGbuY/n5lewxpdQCrfVjgYtMtEdFAdlzpumclIu12+dUKVBacWv1//Ff28gAxNK8wPyOgpeMQaJd8qBAiMv9bF4oViKa16ZavVrrA8ArXoql3erWKYreSSb2H7Gwu/gYpyV3CnRIQoQ0pdQcHEWJXgHeBg4B3YHfA3OUUuHuGskK4UpygsnPhUGqSen5OkcT9gBwwxErn5bcwx6d5scYPCO9z34mY5DoCFwtl4wbhNsllK19RtOzfh2TRwmaUmoT8Drwltb6kG9Cap+Gpyay/4iFjQVlkqAJ4XvTgKe11vfVObYTWKWUKgemA1mBCEy0D85qjUVmCwnRRir9WU4/7CTpvf/KT9FlRNrt/KEkmtfKHsZMnP9i8JDJaJDeZ/XJGCRCnqvEKXPhTIoGzfTaM44Cma9Uk3ugCOJTpQVADU9n0A4BTwJPKqU+wZGsZWutK7weWTszPC2R7C1FbMov4w8jg/cdUCFCRDyO3ouurARu9WMsop1xNqB29jgrO2n127PDqWZ017+xLrqMztU2Linqx/MnbsKGwW8x1JU3fwyj5q9yOXtoUAq71lJe3zUZg0SHlDvVuwlU+swcR8KXVQ5ZgatYG2w83YN2mVKqG45iIDcAb+JYb70M+JfW+jMfxNguDE9LAmBD/pEARyJEh7AOOAtHoaKGzqo5LzqwujNkDROMBbm7A9KAujPl/D3iOYaWf8dj4Z35qXQML1f+2u9x1NVnZg4J0UaMYQqr/ecyKSajgXlXZdT+zrI3FzJq/iqXv88OSsYgIYTPeLwHTWtdAvwF+ItSahDwBxwJ22Sl1AGtg3gBvQ8N7BFHTISB/UcslBytoFunqECHJEQouxP4t1KqGniXn/d/XIejP9EVSqkw58Vaa3tAohQB0XCGrNBsYcbSLcxZvpPZYwcHoNiFndSkHF4++V8GUsohncC2orvYovv5OY7GNI4ZRKNBkeCmIbWr36eU3ZcxSAhvSEkwYa5KcDTR7pPaqiIkydZqFh2w0jPrBx9EGBhtLRKySyn1ZxzrrucDvbwSVTtkCFMMTU1k9Q+lbCoo49dDegY6JCFC2baaf86v+ahLAXXXYGjaONaJ9sXdDFnZSSszlm7xbzCqkj4p/6A07iAvnIDJxX25peoeSkj0bxzNsNo0MZHhbJl9aaNzrn6fUnZfxiAhvMGxz61tRUIyFmfQkyLvBBQkWj1gKKUuxDF7dhWOErPrgXleiqtdGp7mSNA25EmCJoSP/RkC1rZKBLlgKQdvMJbSq/ffKI2sIM5mp7u5D7+vup1KIgIdmkvufm+eHu8gZAwSQviMp1Uch+DYe3Y9jqaM+cBzwBKt9ffeD699cfZD2ygNq4XwKSlfLZri/5L5jXWK3o6p15scMWjSq6z0L/wVL1eMIxiaT6e4+f24K6Pv7vfZkcvuyxgkRPBIjkl2NND2UY+2ZGtNlUknP1Sb9HQGbRtQjmO99RKt9ZfeD6n9GpqagFKwo7CcCquNKGNgqnIJEeqUUguBR7XW+1ycSwNma62n+j8yEQzuyxxYb8+Uv/VOWs7Rbqs5qRQjT1ZRXvgH3q8eGpBYXHH1+2mqjL6n13cEMgYJETxyr3FXUNU1lw24m5CxOMNRZdLJD9Umw5q/pJ7xQA+t9XRJzhqLizIysHscVptme2F5898ghGityUBXN+e6AJP8F4oINuOGpjDvqgwSTEa/P/vKsC+5OOpj7Erx2zLFzvwHWB9EyZlz/m7eVRmkJJhQOGbU6lZsbMj5+2zp9R3EZGQMEkL4iKdl9t/1VSDNUUrFA/8FTgNGaq13BCqWpgxPS+S74mNszC/jrPSkQIcjRChzt/+jB9ChN8d0dM4S+2aL//qbGbEyK/xfTAr/L9YjYDzel8VHbwu6/WYaePD97cy7KsNlE1p3xg1N6egJmSsyBgnRDqUkmEifmdPomLsxMTkmmYzFGSTHJJO7v6ahto95XCREKXUKjjKyqUDDWvJaa32jNwJz4SQwBljgo/t7xYj0RN5YV8CGvDI4P9DRCBE6lFJXAlfWOTRHKVXa4DITcC6w0W+BiaDgTMoCsfesc/RmenV/h3GHCqi0hZNVPYm3Ki4kGPabuSIVGFtHxiAhQoOrRKxhwlaXcwllxuIMKC+ov9zRRzwtEnIFjv1nYUAJUNngEp9VNNJaW4HDSgXnHzyn4amOWbNNBWVorQn2eIVoR1JxvPABx1hzJo3HoErgK+BBP8YlAqxhny7/0ZzS+V1+6rqRvUrxz7hurCm5na1B0N+sOR28AmNryRgkhPALT2fQHgM+ByZorQ+35oFKqTtwrN3OAN7SWk+ucy4JeBW4FCgFHtRav9ma5wRK7yQTXWIjKT1eSd5PJ+nTJSbQIQkRErTWz+GoGotSah8wTmu9tS339HTMUUpF4CiWFKu17rB9H4ONu75nPhVm4dTkFyiMKwUUY44YefvQAxwJsv5m7nTkCoytJWOQEB2bq2qRydZqco+Fe72qo6cJ2inAH1ubnNUowpHoZeJYClDXC0AV0B3HO1M5SqmtWuudbXieXymlGJGWyMqdxWzIOyIJmhA+oLXu46VbeTrm3Idj9UCsl54vvMDfyxpNEfvp0fufFEZUE2ezc96hPrxbPh0b7aNyb0evwOgNMgYJEVpc7UtzHncuiXRVLbJ22aOXeZqgfQd0bssDtdbvAyilRgC17/4opWKAq4EhWuvjwGql1Ac4mmHPbMsz/W14TYK2qaCMa0f0DnQ4QoQcpdR5zV2jtf6imXt4NOYopfrg6AN5D/Bya+IW3pe9uRCF/zoG9zb+QEWflykNU/StspJ44DLerrzYT0/3nLPnmUEpbFqTkmDivsyBsv+sjWQMEiK0uCsQ0tTeNF/yNEG7H/iLUmqd1nqvl2MZANi01nvqHNtKnVIbSqkPcbzDNFAp9U+t9aKGN1FKTQemA6Sm+r7KiivDnA2r86VhtRA+8jnNvyZvbjqj2TGngb8CD9FMdbZgGIM6gkAUBckMW8+CsH/y4rFI8sNi2V10B1vs3ppI8b6mqpKJNvscGYOEED7iaYKWhWMGbZdS6nvgSIPzWmvd2tqFsTiaYNdVDsTVufnlzd1Ea/0S8BLAiBEj/PWmaj1DUjoRER7GnkPHKT9pJT7a/714hAhxF7g41hn4DY4XN3e04B7NjjlONdXbwrXW/1ZKjW7qpsEwBoU6fxcFMRrKuClyGQ/YPwMg4/Bgllhv5jjBu4RdljH6nIxBQoSIUfNXuX2zL8XNft3MZZkUnSgi2Vrtk7L7niZoNmC316NwOA50anCsE3DMR8/zmchwA6enxLMhv4xN+8u4YGC3QIckREjRWv/Pzan3lVLPAmOBj5q5TYvGnJplSE8Czb5BJHzHOWNWZLagFNj99LKzW/QmIlLeYZWuZmqhkRerfsdC268J1hL6gCxj9AMZg4QIHYVmC3nzx3j0PUUniti+z3cl9z1tVD3aJ1E47AHClVL9tdbf1xw7A2g3BULqGp6W6EjQ8iVBE8LPcoC3gduaua6lY05/IB34sqZtRgQQr5QqBkZqrfO8FLdwo+GMmfZLcmZnYJe3KO6yDYtSpFngRus9bLOd4Y+Ht5rC/V4K4TcyBgkh2sTjRtVtpZQKr3muATAopaKAaq31CaXU+8CflVI34dhrdgXwy1Y8Yywwtl+/wPWiGV6zD21DnuxDE8LPBgL25i7yYMzZAdSt9vNL4G/AMKAtFW1FC2V9sNOvZfQNhqMMSPkbB2KOojRkHonik0MzMJPgtxhaS8rnBwUZg4QIcnWXNdZdxuhcutgcXy1tdPI4QVNKpQB/BM7Dsd56rNZ6h1JqBvC11npdM7eYBcyu8/UNwBwc+9tuAxbiKCH7E3Bra0rsa62XA8tHjBgxzdPv9RZnoZAt+81U2+yEG8ICFYoQIUcpNdHF4QhgCHAj8H4Lb+VyzFFKnQt8pLWO1VpXA8V1nn0EsGuti13eUXhV9uZCzBar357X2bSD6F5vcCBck2izMezgYN47NhFN8I/hsu/Mf2QMEqJ9c7essehEEdsn1fQ0y4r32RLG5niUoCmlBgNf4tiL9jUwFMeABJAGnA1c39Q9tNZZOJIxV+eOAOM8iSlYdYmNpEtsBKXHq+j/8Ecky54AIbxpkZvjlcBS4K6W3MTdmKO1/hI3fYa01p9Tp0WI8J66+8ycY+aCXF9te25sZNi3XBf1ErPDoxlSUY218Dqyq8722/PbIjHayOyxg+VvjP8scnNcxiAhRJt5OoP2NLALR5PpChzNFZ2+Ap7wUlxtEgxLHLM3F1J2wvGur8aRqT/4viMjlz+gQrSZq9rmFVrrQ36PRHhFw31mhWYL9727FasfqoEorNxqyOGP4e9isGgOH+zPy+Z7KaWLz5/tDX8Zf6b8XfE/GYOECEbPOhpHZ/ZKpsjoPs2JGwQZixu3WU62VjtmzsCnSxib42mC9ivg91rr40qphv09DgE9vBNW2wTDEscFubuxNdjJbrHaWJC7W/6QCtF23YAL+XlfRgGwCsc4JNqhBbm7G+0z80dylhD1HUkpr/Pr0iIMlZq/VV/Bs+ZrsDXbwio4KORNvwCRMUiIAHNVHj8vqoD0ijeJM87k2K75br832PtEepqgNbXptQvNNE/sSIrc9FNwd1wI0byaPbCvA6NpXOdcK6X+B0zUWh/wd2yibfzZcNrBTr/OyyjtupHDSvFypySOFf6Bz+1D/RxH20hREP+SMUiI4OFyH1kW5M0fQ8bimR6Xzg8mniZo64EpwHIX564D1rQ5ohCRnGBy+YJD/pgK0TpKqQTgcyAJmAl8AOTVnE7HUf3sfuAzpdRZWmuz/6MUrZG9uRCFYzm4PxgMZvqn/J3CmKOA4mJzOF8dvIdDBF9LFJPRwNXDU1j6zX6stvq/IWOYkqIgfiRjkBC+l7kwg6IWLmBwuUyxTyosziA5Jtn7wfmRpwnao8AnSqmPgTdx/D29WCl1F3AljsqOARcMe9DuyxxYbz8FQJQxTP6YCtF6M4E4YJjWOr/Bud3Ak0qpd3EUMJpZ8yHagQW5u/2WnHWN2YAxeRmF4ZBos3HWoYFkl08JqiWNidFGzCet9YpLjUhLYs7ynZSddOxtTjAZyfqtFAXxMxmDhGiCqyWHnoobRJNLE+sK9mWKbeFpo+r/KaXGAX/BURYWYD6Od5DGtaDEvl8Ewx405x/NBbm7a/9jvXBgN/ljKkTrXQnMd/HCqJbWep9S6gngFuTFUdDL3lxI1gc7/VZGf3T4Wr5PeZ+jhjDOtFRjKRzPe9az/PJsT0RHhLP5kUvrHRs3NEX+fgSejEFCNMFd6XpPtPelid7icR80rXUOkKOU6odjk+xPWmv/1UFuR5x/UL/6oZTrX1nHV3t/4nhlNbGRfu8PLkQoSAU2tuC6jTXXiiCWvbnQb1UaI7DyYPibTAnP5ZNSE5+G92ZlyQzKgrTxdKHZwqj5q+q1G5DkLCjIGCQ6nDYvOfRQsq35azqCVmcKWusfgB+8GEvI+kXfzpyVnsg3eWUs/iqP2y8I3NJLIdqxEzj2fjQnETjp41hEGy3I3e2X5Kxn/GdcY/yYKcfzqdIGvi6/irdtlwV142nFz0VTpEVLUJExSIQsd8sTZclhYHjaqHpiE6ftQDmwWaoX1aeU4q6LBnDDq+t45cu9TP5lOjEyiyaEp9YDfwD+08x1E2uuFUGmbiNqn6dmYZX07/EyxfEH+JfW/MLSnfkn72Cb7uvrJzcpMdrImNN78tl3hyk0WxoVR3FVLEVatAQNGYNEyHK3PFGWHAaGp1nCIn7+21G3vGzdY3al1FJgita6biNrvwmGIiENjerXmeFpiWzML+P1r/O5dXRgXyQI0Q79BViplHoKeKjh+KKUigDmAuOAywIQn2hCw0bUvhRr2k1S8hKKI6ox2e1cdjiRm0/cyTFifP5sd9w1k66btLqr/gvSoiVIyBgkglrmskyKThS16nubbNws/M7TBG0U8AaOMvvLcDRk7I6jxP5vgNuAIcAcIB94yGuReiAYioQ05JhF68/Ehet5+cu9TPxFmsyiCeEBrfXHSqlZOKrJTlRK/Zf6Ja4vAToDs7XWHwckSNGIMwHxT58zG2ld38XceTNlStG/0krPg6NZbPkNjVtW+Ze72a+GxT/cLTOSFi2BJ2OQCCauxoq4QUVs31fg3QfFy3bKQPA0Q7gXeFtrXTfx2gN8qZQ6BkzXWl+plOoETCBACVqwOigYwv8AACAASURBVLd/F87sncCW/WbeWJfP9PNkFk0IT2it5yqlvsbRa2gc4HzVagG+ABZorVcFKj5Rnz9nzQzYGN7zOXYllKA0/LYMvi65i032wL+48CQ1dNWixWQ0SIuWICFjkAgWrpYkZiyeCVnlAYpIeJOnu6QvAT51c24VcFHN518Asli+AaUUd13cH4CXvtiLpUpK1QjhKa31Z1rry4BOQA+gJ9BJa32ZvDAKjOzNhYyav4o+M3MYNX8V2ZsLAUchEH8kZ73VId6J+DPPHttM/8oqrijsxzvFj7E3CJIzcOwBcP5OmjNuaArzrsogJcGEwrHpft5VGbL/LIjIGCSE8DVPZ9CqgOG4TtKG15wHR+J3og1xhazRA7pyRq94th4o5411+dx07imBDkmIdklrbQNKAh1HR9dwlqzQbOHupVuYsXSLz5+tDEc5PeldXj/6FfHKwkFrEhF5U1liH+LzZ3vKkyIf0vOsfZAxSHjKk5L1zXG1Z0xK1IcOTxO0d4E5Sikbjj1oJTh6oV0LZPFz8+ozAemN5oJzFm3qog384397mXBOGqYIL/3fKoQQfuZqlsz3xfMhLnYT0cnL2Guw8xEGOpnP5kHrTZQT64ene06KfAjRcTn3i3lSsr45UtI+tHmaoN0DxAFP1nzU9Sbwx5rPdwBfty201gvGKo51XTCwG6f3imfbgXLeXF/Ajb/qE+iQhBCiVfyeeIRV0rv7vzAnfM9J4KyTlaw5ehXLrZkEuhBIU6TIhxAdl3O/mJSsFy3lUYKmtbYANyil/gycg2Pd9UFgndZ6T53rcrwapYeCsYpjXUop7rywPze9voF//O9HJpyTSpRRZtGEEO1PU6XhvS0ieg9dkv+F2VhFhF1z3REDOaX3U6B7+uX5LTGqbxKbCsqlyIcQ7dGzGVBeQGavZIqM3qu07VyOKEsQRUu16r++mmRsT7MXCrcuGtSNwcmd2Fl0lLfXFzB5lMyiCSHah7ql85WfJq1SY9dQ1ns5x4BTK6sYXHwGL5+cQHXr/ox5hckYhsVqBxwNqGePHcy4oSmNepvdlzlQ9pQJEeRGzV/FmooC0iveJM4402tLEUGWIwrPNfuXTSmVChzUWltrPm+S1trLDRhCk1KKOy/qz81LNvLi/37kd2fLLJoQIvg1LAqifbzhLAw7Uw3/396dx8dVl3sc/zzZ0zZLQzcaKMUWLEvLVmUTFFmKIGuvirLKFQSu4IJIUa/iDpd7VXblXjZFVMRSRAQUCiogSFsspVC6UEqb0KalSZo9k8lz/ziTMk0nyUwya/J9v17z6sw5Z848M505+T3n9/s95zG+HHqAi9urmNVqvFT3Oe7xfVP7wgMoLsijanRxzCRMRT5Eck9NQxuUoKGIkhXiOfW4Fjgc+CfBBRkH+nOsLCNOJ+w7kX12Lef1d7bxwKL1nHf41EyHJCLSr3SVzsc6GTfuj/yweRkneFBz6qT1+/DDrnNpIfPzubq7ffvQzpqGNq6Zvwzo+4LUIpI6cx6cQ21L7ZD2UbYPzGQK3DuTyaMnJykykcGJJ0G7EFgTdT8dBbpGBDPji8dO55L7lnD7M2v41Ad2p7hA+a2IZK90FAXJK32LCZN/QUtRK0+UtHBQXQVXhy5iYffBKX/teIW6d/xT2BYKJ1RKX0SSp7allmXnL4sMUzyDqe33J7wPDUOUbDJggubu90bdv6ev7cwsH6hITlgjxwn7TmLGpDJWbGzigUUbOPewPTIdkojIdtHzzfLNUnuGzkJUjn+E7qp/0mIwrbOTveuncELH56mnPJWvnBQqpS+SWdHDFEVyWTxz0LYCx7n7kshjAx4GvuTub0ZtOht4niwY4pjtZfaj5eUFc9Eu+9USbn96NZ+cvZt60UQkK/SebxZO4YSzvNJ1jNv1l7QVN5Pnztn1bazf/G/8KHwU2VY+f+yoQupbQzstVyl9kSEaZBXFyaEuuLaCt0qAigHLJYhkvXi+/ZW9tssDPk5wYeqslO1l9ns7cb9J7D1xDCs3NfPg4g2cfah60UQk89I136y8oAbb43baDKZ2hpi7aSy3N1/NJqpS/tqJ6qnWGJ24gkrpiyRF49twbSO1985k2fnLEn761HmP8ta16j2T3Je5+sSyXV6ecflH9+LyX7/MbU+v4ROH7E5RQV6mwxIZ9sysCrgTOAHYAlzj7jtNXjCzLwFXAOOAZuC3wFXu3pXGcNMuHUP2jsl7mevy/5f7GgHPo27LqVzbdTzZ1msGQRLWU0ofUCl9GTIdg/p25HULE77GYrV6sWWYUIKWJU6auSs3PrWK1XXNzF+ygbM+qC56kTS4FegEJgIHAo+a2VJ3X95ru0eAe9y9IdKgepCgsfTjtEabZqm6CLXltzBm4kOc3bKJeR2LATj63b25MnQpb/vEpL/eYBXmGWNKCmhoDamUvqTKsD8GJVRhcc/3qii+0dCmuWQyYilByxL5ecblH53OF3/zL255ejVzD9mNwnz1oomkipmNBuYC+7t7M/Csmf0BOBeYF72tu6+JfirQDWT/JNchumrO+/nSb/+VxD06BeVLKZs4n66CTv5RGqJ1QyE/7vokd4U/RjfZdcwLdTujigp4+VsnZDoUGYZGyjGotqWWZWuDoYuJmLro0RRFJJL94k3Qqs3sfZH7+VHLGqK22S15YY1MH581mZueWsWazS08tKSGT35g90yHJDKc7Q2E3X1l1LKlwIdjbWxmnwF+BpQRDEW6MuURpljvCo1hd6qjeopOP6iarz24lM7w0IuDWEEDZbv+Dh+zhi7gg23tnFFXwcmdX2Wt7zr0N5MiqswoKTSijkFT5yWWcGm4ooxk8SZoD8ZYtqDXY0PXSBuS/MhctC/9NuhFO+PgavWiiaTOGKD3Kd1GgsbPTiLzQu43s72A84BNsbYzs4uBiwGmTMneocp9VWjsfdHl3tf7Slw3hWNfYPSEPxHO66Is3M0VW5t4vf5ULg/PwbOs16w3VWaUFBpRxyANVxSJXzwJ2mdTHoVs9/FZwVy0tVtaWPByDZ+YrV40kRRphp0urlUONPX3JHdfZWbLgduAM2OsvwO4A2D27NlZe9KqvwqNPRddBhhqZf2qvC2Ujvsj2/K6Oa6llRPqxvGjjiuzaq5ZT+9h77OMqswoKZaTx6CE5pQRKYGv0vciCUnoQtW5Ipeug9ZbQX4eXzhmOlf+bim3Pr2aMw6qpkC9aCKpsBIoMLO93H1VZNkBQO/J+bEUANNSFlkaDDR0L7onLWF57eDGqbaYawvv4dUtXWyjiBe2zeWy8LFZ12vW7c5b1528fcinKjNKmuTkMai2pZbyd27kufYzmNq+U8HJnTRXlsK8j6YhMpHhY1gWCcm166D1dtqBk7l54SreereVPyyt5cyDNb1PJNncvcXM5gPfNbPPEVRQOw04ove2kfV/cPc6M9sXuAZ4Iq0BJ1k8FRoTvwaaU1C2jNKJf+Ajzd3c1LACAGven+u6LmKDjx9ktKnVM4xRlRklnXL5GFTT0AYlGrYokirDMkHLdQX5efzHMdO56sFXuGXhak47sJr8vOy7JpDIMHAZcBdQB7wLXOruy83sKOAxdx8T2e5I4AdmNgbYDPwO+M9MBJwsV815/04XWx4KK3yXkkkPUzAmqHewrbSdhvpSrus6m9+EjyEbr2sGQSl9DWOUDMrKY1B/wxgnh7p4ouQzGrYokkJK0LLUGQdVc/PC1by5pYVHltbqrK5ICrj7VuD0GMv/TjCBv+dxzs/FjTV8b+4h1fz6xfXbC4QMinVRVPU3isc9BXlhysLdfLm+nvKGvflY6ELeYZfkvYkUGFNSoOOrZEy2HoN6SuPHGsL4BnBkZSnPfVnDFkVSRQlaluqZi/a137/CTQtXccoBk9WLJiKD0rtiY01DG1c9uBScoSVnea2MnnobecVbADilqYXzt4a5sf3f+VP3oWRrr1m0htZQpkMQyVoawiiSGUrQstgZB1dz08JVvLm5hT++UstpB+osr4gkLlbFxlASrm12iK+nvGszW6yLb767lTebj+ATXZ9m23sn/jMu34xud/IilRp7Uxl9GakGGsaoIYwimaMELYsVRuaiXTN/GTcvXM0psyaTp140EUlQ0i62bCGKdvkrxU17cnX4Gc4peIqtW/Ko65rIt0Nf40XfJzmvkySlhfn86MyZnH5Q9U69iD3rNf9MRqrallqWnb8Mrq2Aa3e8HNvUeY+Ces9EMkYJWpabe/Bu3LJwNavrmvnTq+/w8VmTMx2SiOSYylGF1A9pKJ9TULac4gl/JK+ogf3Kujj7nVpCns99nadwa9fpdFCUtHiTpSc5A7b/qzL6IrEded3C7ZVdq9WzLJJRStCyXFFBHpcdM41vPPQqNz21ipP231W9aCKSkCFNMyuqo3jiIxSMCS7TtHdHJ1+tr+fl7unMC13ESt89SVEmV3Vl6U7Jl8roy0gXPaxxcqgr6D2LDGWsaWjTnDORLKEELQd84pDduXXhalZuaubx5Rs5aeaumQ5JRHJIY9sges/yOige9ySFVc9h1k1ZuJvL6xv42LYuftx1FveFj6c7yy443cNAQxdFYuipzth7SKOIZBclaDmgqCCPS4+Zzn8uCHrRTtxvknrRRGQHscro9/QWxXNR6t7MOike+w8gzNxtzVxR38jzodmcGDqXjVlcOt+Asw+bop4ykX5ED2fsoWGNItlDCVqO+OTs3bh14WpWbGziz69t5MT91YsmIoFvLljGr154m56RjDUNbVwzfxmL1m3l6RWb407O8ketIdw6lSpa+Hre/ZRs2cieoRCj26v4UtdV/LX7gNS9iSGorizVvDKRfsy5aya1+e9VZ6zZpOGMItlMCVqOKC7I59KPTOPbf1jOjU+tZs5+kzBTL5rISLfg5ZodkrMebaEw973wdlz7yCvaRPHERykYs5L9N83g9pbnqLQWOpoL+Fn4VG7rOjUri4BAkJw9N08XzBXpT20+QcXGHvMezVwwIjKgYZmgmdkpwCnTp0/PdChJ9akP7M5tz6zm9Xe28ZfXNnHCfpMyHZKIZNgNT7yxU3IWL8tvpmj8Xyis/CdmTmkYzsh/nkpr4W/hmXyr6wLe8uztrVeZfJH4TY1KyjScUSS7DcsEzd0fAR6ZPXv2RZmOJZlKCvO55MPT+M4jr3HjU6s4ft+J6kUTGaF65pwlOrcMCK5nVvUcRbs8jeV3YA6faGziPxoaCYXL+Y/QFTzafSjBjK7sYoATNDA1nFEkfhrSKJI7hmWCNpx9+oNTuO2ZNSyv3cZTr9dx3L4TMx2SiKRZrIsuJ6Kg/BWKJzwOwKGtIa7ZupmpnWHuCZ/IT7rm0syoZIY7aJWlhXz8gF15esVmzTETGaQ5d81kQnd3psMQkQQoQcsxPb1o3/vja1z6q8V0hV2NFpER5oYn3kgwOXOs8F08NA6AfbeVMbHM+GzLJo5oa2dR995cHrqA13xqSuJNlHrHRJKnNh+aVv5XpsMQkQQoQctBZcXBf1soHMw86anYBqhBIzLMLXi5JqFhjXklb1M84THyS9dTvOYSruZRPlnwV9gCm7ySL4Yu5OHuI8n0cMZzDpvC90+fmdEYREREsoEStBx041OrdlrWFgrz9YeWsXZLC5MqSphYXsyEshImVZRQNaoooeum9Xc9JRHJnJ6hjX3JM+iOVAzJK6qjaMLjFJa9BkBRuIDrx3yf49q30en5/F/4ZG7tOo0WsqNYwNMrNmc6BBERkaygBC0H1fZx9ry1MxwzeSvIMyaUFTOhvIRJ5ZHkbfv94PHEihLKigt4+F+1O8xtUe+cSPYYaGhjt4MVNFA0/kkKKxZj5lh3Pmc2hriycS1l7jwZPojvd52TddUZ+zquiYiIjDRK0HLQ5MrSmEOcKkoLOfewPdi0rZ1NTR1samxnU1M7Da0hahvbqW1s73e/pYX5dIa7CXfvWLS7LRTmO48sp3JUIRWlhVSOKqKytJDy0kLyE+iZA/XOiQxFPElMyaQFFJStADcObyzhB42rGR/u5s3uSVzedS7PdB+UhkgTN1llv0VSRmX1RXKLErQcdNWc9+9Uwa20MJ/vnLpfzGSnPRSmblsHm5ra2bStnY2N7dQ1dexwf2Nje79n5utbQ1xw90s7LS8rKXgvcSstomL7/cLtyytKi6gcVcjidfXc/NQq2ruCalLqnROJra8TGbFOzlh+M+R14KFdgsdbjuZ9bOS/61eyV7idZi/hR11ncFf4Y4RiHPLzgFTUd+urAuOidVt3urC2rmcmkhqrtixnduh9PKWLuYvkFCVoOagnmYm3J6qkMJ8pu4xiyi59l852d5o6ujjuf/5KXVPHTutLC/OYPbWKxrYQDa0hGlo7aerooqk9uK1ncMOT2kJh5s1/hUXrtjKhrITxZcWMH1PMhPJixpcVM25MMYX5eXHtS71zMhz0LqEffSIj+uSM5TdTWPV3iqqeJ9y6Jx3rL+CT+c/wFX+QCVsaAPh9+ENcH/o0dYyN+VqpSM4MOLuPgh8LXq7h94trdkjODJh7SLV+qyIp0GmwZrUqOIrkGiVoOer0g5LboDEzyksK+fpJ+8TsnfvRmTN3er1wt9PUHiRsjW0hGtqCfxtbO2MsC/HPt7bGfO32UDf3vfB2n7FVjS5i/JggYZtQFvwbfZtQVsw/127le398jbaQeuckt8WaZ9YWCnPDE2/w3LyP0tLVyHf+dhtFVf/A8joBqKKZ24uvYZatB+Bf3dP4buhclvje/b5WMpKz0sI8qkYXx3ViJNZ7c1QgRCSVNLxRJPcoQZMdJNI7l59nwXy0UUVx7fvI6xbGnDtXNaqQK47di83NHWxu6qCuKfh3c1MHW5o72NrSydaWTt7Y1JTQe2kLhbnu8RVK0CSn9DXPrHbbu/x08U+5f+X9FI8Ltilp3oMr6rdwbtcLYLDex3N96Cwe88MJe8zdJMSAylGFuENDWyjmNu2hbp6Lc/hUn+9NBUJEUibe36eIZA8laLKTZPfO9ehr7ty3Tok9dw6CXrqtLZ2RxK09SNyaO6jb1rE9odvc1MHaLS0xn7+xsZ2P3PA0B+8xlkMit70nlCV02QGRdOqrCNCkygLue/0+OsIdzBh1AB9auZnLO58lz5xtPoqf+ZnMOPWr3DJ7TwCmXfMnwh47SystzKe4IK/PpKvH2utO3n6/rxMsiRT36Ou9qUCIiIjIe5SgSdokOncOgl66nqGM+1Le53Z9NR4NeOvdVt56t5X5S2qAoLDJQVPGcsiUIGE7cEolY4r1U5DssP1ERvc2iioW0bn1Q5QWFnH18R+ku/grTH97EbMW/xpCrXRZPvd0Hc/vRn+az5/4AU6N+i31lZxB0Lts9N/F1ntYVF8nWBIp7pGMfYiIiAx3apVKWqW7d+4Hp+/H3pPKWbyunkXr6lmyrp6ahjb+tnIzf1sZzHvJM5gxqXx7D9she4xlt7GlmL3Xy6YCJJIuH5pRzEcO/yfPbvoD5HVQUVzFN476DKd3/wUe/QG01AUb7nMqBcddy4W7TOPCGPup7qO3qkdrqO8ZaLGSpsGcYOktGfsQEREZ7sz7OcuabczsB8DRwCbgPHdv7W/72bNn+6JFi9ISm2RevEnUO41tLFnXwKJ1W1myrp7ltdvo6nXttwllxduTteaOLn721zW0RzVo+yqcIv0zs8XuPjvTcaRLIsegmuYa7n71bh5a9RCd3UHxjyN2PYIrxh/Kfv/4X6h7LdiwejbM+QFMOazf/fWuBpmIn37qQH23ZVgaaceg0j1LvW2t5niKZIt4j0E504NmZvsD09z9KDO7BLgQuCXDYUkWibd3bteKUk6eVcrJs3YFoK0zzCsbGrb3sC1+u566pg4ee3Ujj726MeY+ei7ePbq4gLKS4FZeUkhZSQFjigsoiPPSAKDeuUwysyrgTuAEYAtwjbvfH2O7q4DzgT0i293m7jckK44bl9zI3a/eTdiDZOrYKcfyuckfYf8X74HnvxxsVDkFjrsW9jsTbOA5lL17q+I9FVddWarvn0iaZMsxSESyS84kaMBRwGOR+48B16METZKgtCifQ9+3C4e+L7jQr7vz5pYWFr9Vz+J19fx20fqYz6tvDXHRL2L3jpQW5m9P3MoiiVtPAhckccH9VXVN/H7xBjojJfdqGtqYN/8V3J0zDt4tNW9Yot0KdAITgQOBR81sqbsv77WdAecBrwDTgD+b2Xp3/00ygpg0ahIAp047lQunzGHaol/BXy8AHEoq4Oir4IMXQ0FxQvuNPmnR1zzNaJoPJpJ2KT0GTehKvAddRDIv7UMczewLwAXATODX7n5B1Lo+zySZ2deB19x9gZlVAg+4+wn9vZaGOEoy9NWwLS3M5/Bpu9DUHtp+we6m9hDNHV10J+FnNam8hLGji6gaXUjlqCKqRhUxdnQRY0cVUjW6iLGjiqgaXURl5HFpYf4O8+Z6ZFMPXTYNLzKz0UA9sL+7r4ws+yVQ4+7zBnjuTQTHz8v72y7eY1BHuIMtW1ZSveQ+eOlO6A5BXiF84HPw4a/BqKq431dfYg15LMwzxpQU0NAayvh3QyQdRtwxaHK+L6pVkiaSLbJ5iGMt8H1gDtC7tnJ/Z5LqgYrIdhVA7KseiyRZXwVI+pqD5u60dIajErcQ26ITuKj79/5jXZ+vu3FbOxu3tccdZ3FBHmMjSVxPUtfQ0smLa7dun2OnHrod7A2EexpGEUuBD/f3JAuy4KOAn/ex/mLgYoApU6YMHEX7Nor/cQvVz98CoRa6MRaEP8T9RedwzqSjOT0JyRnELtBxzIzxPL1iMw2t/ZfbF5GUSP0xqEKXlBHJRWlP0Nx9PoCZzQa2txAjZ5LmEpxJagaeNbM/AOcC84BngWuAewmSu+fSHLqMUIlWnjMzxhQHc9F2rYi5yXZPvl4X+7pQFSX87tIjqG/ppL41uFB3cD+0/XFDayhYHnnc0dUdV1LXHurmKw8s5dZn1jC5spTqylKqK0uoHlvK5IpSqseWMqm8JK55dNnUOzcIY4DGXssagbIBnnctkAfcHWulu98B3AFBD9pAQbz526t539pgysmT4YP4765PscKnQAiWz18GkLTPNHrIY+8etZqGNq5J8uuJSL9SfwyanJ87leBEZLtsmoPW75kkd19mZuvM7O9AHcFY7J0kfPZaJA7pvjzA106cEUmc4r+Ab1tnmK2tnTskdV/8zb9ibuvA6rpmVtc1x1yfZ8EQy8mVQcL2XiL33uMnX9uU6w38Ztjp4nrlQFNfT4gM0T4POMrdO4YawIKXa7hp9ZF8j6X8pGsui3zGDuvbQmFueOKNlHyeNzzxxk4VHlP5eiKyk5Qfgzb52CEFKCKZkU0J2oBnktz9moF2kujZa5FMSuZ1oUqL8qku2jGp+6/H34jZQ7drRQl3f/YD1NS3UdMQudW3URu5X9fUQW1jO7WN7SxaVx/z9cyg9xTWHGvgrwQKzGwvd18VWXYA0HtyPgBmdiFBb/7R7r4hGQHc8MQb1ITGcjbf6HOb2gEKewxWX/tN1euJyE5SfgyqQwmaSC7KpgQt4TNJIsNBqnrnoO8euqtPnMGMSeXMmNT7Jxfo7OpmY2M7GxpaqW1o3yF56/m3oyv2hY5zpYHv7i1mNh/4rpl9jmDe62nAEb23NbOzgR8Cx7j7m8mKIZ7PanICvaiJmNzHhaxT9XoisqN0HIMmlCVW+VVEskM2JWgJnUkSkYENtoeuqCCPKbuMYsouo2Kud3cO/9HCmPPdcqyBfxlwF8Gw6XeBS919uZkdBTzm7mMi230f2AV4KapS5n3ufslQXryvJKmHQcrK3veVvKvMvkhapfQYNLG8JDVRi0hKpT1BM7OCyOvmA/lmVgJ0JXImKY7XOAU4Zfr06UmMXCQ3paKHzsyY97EZOd/Ad/etwOkxlv+dYNh1z+M9U/H6sZKkHgacfdiUlPWuJnN4rYgMTqaPQSKSnTLRg/ZN4NtRj88BvkNQlSjmmaREX8DdHwEemT179kVDjlZEYlIDf+iiP8OahjbyzQi7U52mzzKVw2tFRERkcDJRZv9agmQs1rqYZ5JEJDupgT90+gxFREQk2sAXOspBZnaKmd3R2Ni7KKSIiIiIiEj2GpYJmrs/4u4XV1QMcJVgERERERGRLDIsEzQREREREZFcpARNREREREQkS5i7ZzqGlDGzzcC6ODcfB2xJYTiporjTK1fjhuyIfQ93H5/hGNImwWNQNsuG704u0ecVv3R/ViPtGNQEvJHpONJkJP3uRtJ7heH1fuM6Bg3rBC0RZrbI3WdnOo5EKe70ytW4Ibdjl8zSdycx+rzip88qtUbS56v3OnyNtPcLGuIoIiIiIiKSNZSgiYiIiIiIZAklaO+5I9MBDJLiTq9cjRtyO3bJLH13EqPPK376rFJrJH2+eq/D10h7v5qDJiIiIiIiki3UgyYiIiIiIpIllKCJiIiIiIhkiRGfoJlZlZk9ZGYtZrbOzD6T6ZgGYmbFZnZnJN4mM3vZzD6W6bgSYWZ7mVm7md2X6VjiZWZnmdnrke/KGjM7KtMxxcPMpprZn8ys3sw2mtktZlaQ6bgk88ysudctbGY3R60/1sxWmFmrmT1tZntErSs2s7vMbFvke/WVzLyL9Onvt2RmB5rZ4shntdjMDox6npnZ9Wb2buT2X2ZmmXsn6WFm+5jZQjNrNLPVZnZG1Dp9t1IoF9s2gzEc2kODkYttqMHK1bbXUI34BA24FegEJgJnA7eb2X6ZDWlABcB64MNABfCfwANmNjWDMSXqVuClTAcRLzM7Hrge+CxQBhwNvJnRoOJ3G1AH7AocSPC9uSyjEUlWcPcxPTeCY2Ab8DsAMxsHzCc4vlQBi4DfRj39WmAvYA/gGOBrZnZi+qLPiJi/JTMrAh4G7gPGAvcCD0eWA1wMnA4cAMwCPg58Pr2hp1ckcX0Y+CPB9+di4D4z21vfrbTIxbbNYAyH9tBg5FQbarByYV+FeAAAC69JREFUvO01JCO6SIiZjQbqgf3dfWVk2S+BGnefl9HgEmRmrwDfcfffZzqWgZjZWcCZwGvAdHc/J8MhDcjMngfudPc7Mx1LoszsdeBKd/9T5PENQLm7D+sGoiTGzM4Hvg1Mc3c3s4uBC9z9iMj60cAW4CB3X2FmNcBn3f3PkfXfA/Zy97My9BZSrq/fEvB74G5gN4/8UTWzt4GL3f3xyPHjHne/I7Lu34GL3P2wTLyPdDCz/YEXgLKoz+TPwIsEDWp9t1JkOLVtBiOX2kODkYttqMHK5bbXUI30HrS9gXDPASxiKZBTZ5nMbCLBe1me6VgGYmblwHeBKzMdS7zMLB+YDYyPDNPZEBnaVJrp2OJ0I3CWmY0ys2rgY8DjGY5Jss/5wC96GtMEx8GlPSvdvQVYA+xnZmOBydHrycFj5yD09VvaD3gl6rMDeIX3Po8dPktGxmcVawinAfuj71aqDYu2zWDkUntoMHKxDTVYw6DtNSQjPUEbAzT2WtZI0I2aE8ysEPgVcK+7r8h0PHH4HsHZkPWZDiQBE4FC4N+AowiGNh0EfDOTQSXgrwR/mLcBGwiGEy3IaESSVcxsCsEQoXujFvd3fBwT9bj3uuGsr9/SQH9Leq9vBMYM83loKwiGg15lZoVmdgLBd2wU+m6lWs63bQYjB9tDg5GLbajByvW215CM9AStmWB4SrRyoCkDsSTMzPKAXxKMM/9ChsMZUGTS/HHATzIdS4LaIv/e7O7vuPsW4MfASRmMKS6R78gTBPM9RgPjCObIXJ/JuCTrnAc86+5ro5b1d3xsjnrce92wNMBvaaC/Jb3XlwPNvXrchhV3DxHMuzsZ2Ehwxv8BgsRW363Uyum2zWDkWntoMHK4DTVYOdv2SoaRnqCtBArMbK+oZQeQA13jkTOvdxKcYZgb+WOY7T4CTAXeNrONwFeBuWa2JJNBDcTd6wkaFbnYmKoCdgducfcOd3+XYK7MiDjASdzOY8feMwiOgwf0PIjMa5kGLI/8Jt6JXk+OHDuHoL/f0nJgVq8esVm893ns8Fky/D8rANz9FXf/sLvv4u5zgPcB/0TfrVTL2bbNYORoe2gwPkIOtqEGK8fbXkM2ohO0yLj3+cB3zWy0mR0JnEZwFibb3Q7sA5zi7m0DbZwl7iD4I3xg5PYz4FFgTiaDitPdwOVmNiEyR+JLBNXJslrkjNNa4FIzKzCzSoK5Rkv7f6aMFGZ2BFBNpHpjlIeA/c1srpmVAN8imGfVM3ToF8A3zWysmc0ALgLuSVPYaTfAb+kZIAxcESn73XMGf2Hk318AXzGzajObTNCbdE86488EM5tlZiWROXtfJah+eQ/6bqVUjrdtBiMX20ODkcttqMHKybZXMozoBC3iMqCUYKz8r4FL3T2rzzJZcL2YzxP8QDfae9cwOjvDofXL3VvdfWPPjWAYRru7b850bHH4HkFJ25XA68DLwA8yGlH8zgROBDYDq4Eu4MsZjUiyyfnAfHffYfhT5Hc5l+B7Xg8cCkRX0fs2QWGHdQRzs25w9+FefCbmb8ndOwmG850HNAAXAqdHlgP8HHgEWAa8StCo+nl6Q8+Icwl6w+qAY4HjI72P+m6lXs61bQYjV9tDg5HjbajByuW215CM6DL7IiIiIiIi2UQ9aCIiIiIiIllCCZqIiIiIiEiWUIImIiIiIiKSJZSgiYiIiIiIZAklaCIiIiIiIllCCZqIiIiIiEiWUIImSWdmF5iZm9n0JO3vGTN7Jhn7EhEREcmEqPZRz63JzJaa2RfMrCCyzTNR68NmVm9m/zKzm81svzj2GX07Lv3vUpKhINMBiIiIZCMz+zRwdeThHkArwUWqAW5y97syEpiI5LpPABuA8sj9m4EJwLci618huAA3kW32By4ELjGzL7r7bf3sM9prSY5b0kQJmoiISAzu/mvg1wBm9jzwE3f/XWajEpFh4F/uvjpy/8+REUdf4r0ErcndX4ja/s9mdjPB8ehmM3vJ3V/qZ5+S4zTEUVIu0l3/rJkdZ2ZLzKzVzF41s9NjbHuWma0wsw4zW25mZ/Sxz3FmdruZ1US2XWFmF0etz4u87ltmVhG1fKaZtZnZDal5tyIy3JhZHjALeLmP9R/tY3hRTXojFZEc9RJQZmYT+trA3UPAZUAXcEW6ApPMUA+apMs04EbgR8AW4ErgQTOb0XPGJzJW+n7g0cj68ZHnFAJv9OzIzMqB54BS4FpgLTAHuN3Mit39ZnfvNrNzgKXAz4GzzKwU+A2wHPhGyt+xiAwX7wfCwJo+1r8OPACcADwI/ANYD2xMS3Qikuv2JDjGNPe3kbvXmdki4MgYq/N75rG9t7mHkxijpJESNEmXccDR7r4KwMyWAO8AnwR+GNnmO8AK4DR3745s9zrwAlEJGvBFgvkgM3v2BzxpZpXAt83sdnfvcvcNZvY5YL6ZPQEcHnnewe7emco3KyLDysEEw4e8j/VfBQ4A9nF3JWUiMpCeZKqMoB10JvCIu7ea2UDPfZvgmNTbil6PnwM+NNRAJTOUoEm6rIpKpnrOAtUBUwDMLB/4AHBdT3IW2e5FM3ur175OBF4E1vY6W/QE8DlgX4IJtrj7Q2b2c+B2oBi40N1XJvvNiUj2M7P/IzhBc3DUssUEjZ3d3X1DZNl/A3OB6ZEz0AcBS/rZ9QXAt5SciUicopOpbuBXBHPQ4mFArJNFZ7BjkZCmwYUm2UAJmqTL1hjLOoCSyP1xBEMZN8XYrveyCcB0INTHa+3S6/G9BNWQ6giGUIrIyFRPcMYaCOaOEQxfdGAssMHMRhFUS/te1PCgg4F7+tlvA/B5M1tPkMhtcff25IcvIsNETzLVBKxL8HixO8EIpN5eVZGQ4UNFQiRbbCFIuCbGWNd72bvA8wQ9brFui3o2jDS27gJeBSqA65IduIjkjB0SNIK5rncRzPsYG1l2HpAP3Bm13YH034N2EsHx62GCuWdtZjYpSTGLyPDzqrsvcvc3EknOIkVEZgPPpi40yQbqQZOs4O5hM3sJ+DczuzZqDtqhwFRgXdTmjwOXA2+7e90Au74RqCZoYH0c+KmZPeHujyf7PYhI1qsnuKYQZrYPQVGP9xPM/+hJ0L4A/J+7b+t5krtX9bVDM5tJUHyoCPg+QaXHOg13FJFkMrNC4DaCtvtNGQ5HUkwJmmSTbwN/BhZE5o2NJygc0ruh8xPgU8DfzewnBAVERgMzgKPc/TQAM5tLMCftXHd/E7jJzE4A7jGzWXEkdyIyvNQDpZE5r1cCD7v7m2bWCIyNVJKdAZwcz84iPfSPEiRln9KwRhFJkjIzO6znPjAT+CzBCaXL3H1xxiKTtFCCJlnD3Z80s7MJSufPB1YTTJr9Yq/tGs3sCIILOl5N0EPWQJCo/R7AzHYH/hf4lbvfF/X0zxIUELnHzE7upyqbiAw/9ZF/pwPnAMdEHm8j6EE7E5jv7utiPDeWkwjmgxyv5ExEkmgWweU6nGCe2lrgGeAsd1+ewbgkTUztUxERGQnM7IMEFWB/Aezl7kdElj9OMA/208CR7v5C1HO+DpwP7AWc6e4LotZ9GfgxQbXHvq6RJiIikhAVCRERkZGipwftbOB/opZvA84CXoxOziKeIugp+1uM/T0DdAGPmNklZnacme2V3JBFRGSkUYImIiIjRU+C9jawIGp5I0Hlxh/3foK7v9hX75i7v0xQaGQt8F2CAkYrI3NoRUREBkVDHEVERAZgZs8AP40e4hhjGyNI8s7vr/KjiIhIf9SDJiIikhzTgA8Cf8l0ICIikruUoImIiCTHAwTDHS/JdCAiIpK7VGZfREQkCdz94EzHICIiuU9z0ERERPpgZt8k6BEbT3A9onZgtrtvzGhgIiIybClBExERERERyRKagyYiIiIiIpIllKCJiIiIiIhkCSVoIiIiIiIiWUIJmoiIiIiISJZQgiYiIiIiIpIllKCJiIiIiIhkCSVoIiIiIiIiWUIJmoiIiIiISJZQgiYiIiIiIpIl/h/4d9251uCzzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x252 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "i_sample = 259\n",
    "x = dataset.x_test[i_sample]\n",
    "x = np.expand_dims(x, axis=0)\n",
    "y = model.model.predict(x)\n",
    "print('sample {} | model prediction | ground truth'.format(i_sample))\n",
    "dataset.decode_predictions(y)\n",
    "dataset.decode_predictions(dataset.y_test[i_sample])\n",
    "\n",
    "subspaceplot.imshow(np.squeeze(x + dataset.x_train_mean), figsize=(2, 2))\n",
    "\n",
    "AS = NNSubspace(model=model.model, x=x, x_train_mean=dataset.x_train_mean)\n",
    "\n",
    "for sigma in np.array([50]):\n",
    "\n",
    "  AS.sampling_setup(num_gradient_mc=1000,\n",
    "                   num_rs_mc=50000,\n",
    "                   seed=7,\n",
    "                   bool_clip=True,\n",
    "                   sigma=sigma / 255,\n",
    "                   num_eigenvalue=10)\n",
    "  AS.run()\n",
    "\n",
    "  AS.threecolumnplot(title='$\\sigma = $'+str(sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bootstrapping**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "N_bootstrap = 10000\n",
    "\n",
    "n_mc_list = np.logspace(1,np.log10(AS.num_rs_mc), num=50, base=10, dtype='int32')\n",
    "\n",
    "mean_list = np.empty_like(n_mc_list, dtype='float32')\n",
    "std_list = np.empty_like(n_mc_list, dtype='float32')\n",
    "\n",
    "mean_list_boostrap = np.zeros(N_bootstrap, dtype='float32')\n",
    "\n",
    "for index, n_mc in enumerate(n_mc_list):  \n",
    "  for i_bootstrap in range(N_bootstrap):\n",
    "    mean_list_boostrap[i_bootstrap] = np.mean(AS.y_samples_rs[np.random.choice(AS.num_rs_mc, n_mc, replace=True)])\n",
    "\n",
    "  mean_list[index] = np.mean(mean_list_boostrap)\n",
    "  std_list[index] = np.std(mean_list_boostrap)\n",
    "\n",
    "n_mc_list_mc = np.copy(n_mc_list)\n",
    "mean_list_mc = np.copy(mean_list)\n",
    "std_list_mc = np.copy(std_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEVCAYAAAAo63jjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYVNW57/Hv2900DTTdNDKoDCLatAM4jyjaQY9JNLmCJpo4Gw2JPjEnw0lO9F5zDNGTo7knngxmIIZoTNQ4YCTilKjlkGgcbgRFA45EQECZGxm73/vHqpKiqB6qu2rv6qrf53n2U73XXrXrrSX223uvYZu7IyIikg8VcQcgIiKlQ0lFRETyRklFRETyRklFRETyRklFRETyRklFRETyRklFRETyRklFRETyRklFRETypiruAKJWVzfI99hjb/r2jTuS0rZhwwYGDBgQdxhlQW0djXJv5xdeeOF9dx/aWb2ySypDhgzn9NOfZ5994KMfhYaGuCMqTYlEgubm5rjDKAtq62iUezub2aKu1Cvb21//+AfceCOsWxd3JCIipaNskwrAhg1wxx2wbVvckYiIlIayTioAixfDAw/EHYWISGko+6QC8MILYRMRkZ5RUkm6//5w1SIiIt1XdqO/2tPaGvpXpk2D2tq4oxGRQli3bh0rVqxg69atOb+3vr6eV199tQBRxa9Pnz4MGzaMurq6Hp8r8qRiZoOBXwEnAe8Dl7v7rR3UrwbmAbXuPjLL8fOBm4DPu/uNPYlt3Tq480447zyorOzJmUSk2Kxbt47ly5czYsQI+vXrh5nl9P7169czcODAAkUXH3dn48aNLFmyBKDHiSWO2183AFuA4cDZwM/MbP8O6n8DWJHtgJk1AJcD8/MV3KJF8PDD+TqbiBSLFStWMGLECPr3759zQillZkb//v0ZMWIEK1Zk/VWbk0iTipkNAE4HrnT3Fnd/CpgNnNtO/T2Bc4DvtXPK7wE/Ilzx5M3f/gZz5+bzjCISt61bt9KvX7+4wyha/fr169ZtwUxR3/4aB7S6+8K0srnA8e3U/zFwBbAx84CZHQEcBlwKnNHRh5rZNGAawJAhQ2lqSnQa6CuvwIoV0KdPp1Uli5aWFhKJRNxhlAW1ddfU19fT0tLS7fe3trayfv36PEZUfDZt2tTjf0tRJ5VaYG1G2VpgpxuVZjYVqHL3e8ysOeNYJfBT4DJ3b+vsUtbdZwAzAMaObfIFC5o7rJ+yfHnouO/fv0vVJU25L2kRJbV117z66qs96hMp1T6VdDU1NRx88ME9OkfUfSotQGYvUB2wQ/pP3ia7DrisnfNcCsxz96fzHmGaNWvgrrugra2QnyIiUjqivlJZCFSZWaO7v5YsO5CdO9obgTHAk8mrkGqg3syWAUcBJwDHm9nJyfqDgYPN7CB3/1I+A37zTXjoIfj4x/N5VhEpFldd1bV6mzdX52V1865+XroLLriAm2++mYsuuogbb9xxkOs3v/lNvv/973PKKadw3333AeGq6rrrruPuu+/mrbfeoq6ujn333ZcvfOELnHnmmVRUFO56ItKk4u4bzGwWMN3MLgYOAk4FJmZUfRkYlbY/EfgJcAjwHnABUJN2fBZwF2Goct797W+wyy5wxBGFOLuISOdGjRrF73//e374wx9+uAT/tm3buOWWWxg9evSH9dasWcOxxx7L6tWrufrqqzniiCOorq7mqaee4rvf/S5HH300Y8aMKViccUx+vBSYSRgmvBK4xN3nm9kk4AF3r3X3bcCy1BvMbBXQ5u6psjXpJzSzLcA6d8/sr8mbBx6AwYNh770L9QkiIu074IADWLp0KXfccQcXXnghAHPmzKGmpobjjjuOlStXAnDFFVfw1ltvsWDBAkaO3D61r7Gxkc9+9rMFjzPyeSruvsrdp7j7AHcfnZr46O5PunvWuezunsg28THteHNPJz52xj1MjMzDMG4RkW656KKLmDlz5of7M2fO5MILL/xw3k1bWxu33347Z5999g4JJaWmpoaampqdyvNJa3/lYPNmuPVW6MGoRBGRbjvrrLN4/vnnee2111i2bBkPPvggF1xwwYfH33//fVavXs2+++4bW4xa+ytHa9bA7bfD+edrDouIRKuhoYGpU6cyc+ZMBg0aRHNz8w79Ke4eY3SBkko3LF4M994Lp58OWu1BRKL0uc99jvPPP5/a2lqmT5++w7GhQ4fS0NAQ68KXuv3VTS+/DE88EXcUIlJuTjjhBKqrq3n//feZMmXKDscqKio488wz+d3vfsfiLM/y2LRpE5s2bSpofEoqPfDYYzA/b0tZioh0zsyYN28eb731Fn2zTJz5z//8T0aPHs2RRx7Jr3/9a+bPn8/rr7/OLbfcwqGHHsqyZcuynDV/dPurh+65BwYNghEj4o5ERLqjq5MR16/fwsCBeZj9mAcdLRfT0NDAM888w3XXXce1117L22+/TV1dHfvttx9XXnnlDn0whaCk0kPbtoWO+89/HvLwfBsRkZ3cdNNNOR2vr6/nmmuu4ZprrilcUO3Q7a88WL8ebrsNtmyJOxIRkXgpqeTJu+/CH/4QJkmKiJQrJZU8euWV0HkvIlKulFTy7IknYN68uKMQEYmHkkoBzJ4N77wTdxQikqkYZpwXq3y1jZJKAaRGhCUXDRWRItCnTx82btzpyeSStHHjRvrkYe0pJZUC2bABZs6EpUvjjkREAIYNG8aSJUv44IMPdMWSxt354IMPWLJkCcOGDevx+TRPpYA2bICbboIzz4S99oo7GpHyVpecSLZ06VK2bt2a8/s3bdpU8GXj49KnTx+GDx/+YRv1hJJKgW3ZEpbLnzIFJkyIOxqR8lZXV9ftX5yJRIKDDz44zxGVHt3+ikBrK9x9NzzzTNyRiIgUlpJKhB58EP70J02QFJHSpaQSsb/8JYwM27w57khERPJPSSUGCxbAr34Fq1fHHYmISH5FnlTMbLCZ3WNmG8xskZmd1Un9ajP7h5ktTisbYmZ/MbOVZrbGzJ42s2MKH33+rFgBM2bAW2/FHYmISP7EcaVyA7AFGA6cDfzMzPbvoP43gBUZZS3A54ChQANwLfBHM+tVo9k2boRbboHnnos7EhGR/Ig0qZjZAOB04Ep3b3H3p4DZwLnt1N8TOAf4Xnq5u29y9wXu3gYY0EpILoMLGX8htLXBnDlw331hlJiISG8W9V/244BWd1+YVjYXOL6d+j8GrgCyrq1gZvOAfYA+wI3unnlFk6o3DZgGMGTIUJqaEt0KvpDWrw9L5zc0QEUJ9HS1tLSQSCTiDqMsqK2joXbumqiTSi2wNqNsLbDTszHNbCpQ5e73mFlztpO5+wFmVgNMBarb+1B3nwHMABg7tskXLMh6uqIwaFCYgb/bbnFH0jOJRILm5ua4wygLautoqJ27Juq/iVuAzOmsdcD69ILkbbLrgMs6O2HyVthtwLfM7MB8BRqXNWvCyLAXX4w7EhGR3EWdVBYCVWbWmFZ2IDA/o14jMAZ40syWAbOA3cxsmZmNaefcfYCxeY02Jtu2hVth99+vfhYR6V0ivf3l7hvMbBYw3cwuBg4CTgUmZlR9GRiVtj8R+AlwCPCemR1FiP1ZoBL4MmE02d8K+w2i9eyz4THFZ5wBA3e6QSgiUnzi6BK+FOhHGCZ8G3CJu883s0lm1gLg7tvcfVlqA1YBbcn9VqAvYWjySmAJcDJwiruX3ELz77wDv/gF/POfcUciItK5yOd1uPsqYEqW8icJHfnZ3pMARqbtP064bVYWWlrCEvonnghHHw1mcUckIpJdCQxeLQ9tbfDww3DbbfDBB3FHIyKSnZJKL7NwYbgd9s47cUciIrIzJZVeaO1a+PWvw4rHWkZfRIqJkkov1dYWns2i22EiUkyUVHq5hQvh5z/X6DARKQ5KKiVg3bpwO+yxx8IVjIhIXJRUSoQ7PP54SC56+JeIxEVJpcS88064HTZvXtyRiEg5UlIpQZs3w6xZcPfdsGlT3NGISDnpVU9KlNy89FLowD/pJNhvP83EF5HC05VKiVu7Fu68Myzz8u67cUcjIqVOSaVMLFoEM2bAvfeGtcRERApBt7/KiDv8/e8wfz5MmhQWp6zSvwARySNdqZShLVvgkUfgJz8Jo8S01IuI5IuSShlbsyaMEpsxA954I+5oRKQUKKkI774Lt9wSNnXmi0hP6I66fOiNN8J2wAEweTIMGhR3RCLS2yipyE7mzYOXX4YJE+CYY2DYsLgjEpHeQklFsmprg7lzw9bYGJLLHntoAqWIdExJRTr12mthGzEiJJd99oEK9caJSBaR/2ows8Fmdo+ZbTCzRWZ2Vif1q83sH2a2OK1snJnda2bvmdkqM3vIzJoKH315W7IE7rgDbrgBFizQUGQR2Vkcf2/eAGwBhgNnAz8zs/07qP8NYEVG2SBgNtCUPM+zwL35D1WyWbkyPHHyt7+FFZn/ZUSkrEWaVMxsAHA6cKW7t7j7U4TkcG479fcEzgG+l17u7s+6+6/cfZW7bwWuB5rMbJfCfgNJ98Yb8LOfwZw5eqSxiARR96mMA1rdfWFa2Vzg+Hbq/xi4AtjYyXmPA5a5+8psB81sGjANYMiQoTQ1JXKJWTqxbh3Mng0DB8KAAaGspaWFRCIRa1zlQm0dDbVz10SdVGqBtRlla4GBmRXNbCpQ5e73mFlzeyc0s5GEW2pfa6+Ou88AZgCMHdvkCxa0ezrpoUGD4MgjwT1Bc3Nz3OGUhURCbR0FtXPXRN2n0gLUZZTVAevTC5K3ya4DLuvoZGY2FHgY+Km735bHOKWb1qyBhx6C5cvD65o1cUckIlGK+kplIVBlZo3u/lqy7EBgfka9RmAM8KSFiRHVQL2ZLQOOcve3zayBkFBmu/s1kUQvXeYOTz8NzzwD++4bVkQeOVLzXERKXaRJxd03mNksYLqZXQwcBJwKTMyo+jIwKm1/IvAT4BDgPTOrAx4C/uLu3yp85NJd7vDKK2HbfXc4/HAYPx769Ik7MhEphDiGFF8K9CMME74NuMTd55vZJDNrAXD3be6+LLUBq4C25H4rMBU4HLjQzFrSttExfB/poqVLw0PCfvADePhhWLUq7ohEJN8in1Hv7quAKVnKnyR05Gd7TwIYmbZ/M3BzgUKUAtu4Ef7617DtvXe4emls1Cx9kVKgZVokVq+/HraBA2H//cOtsREj1Pci0lspqUhRWL8+dOo/8ww0NIQVkseP1wrJIr2NkooUndWr4YknwjZ8eEgu48eHZCMixU1JRYra8uVhe+SRcFts/Phwm6wuc7aTiBQFJRXpNZYsCdvDD8Po0SHBjBsH9fVxRyYiKUoq0uu4w6JFYZszBwYPhj33DNuYMVCbdQyhiERBSUV6vVWrwvbCC2F/6NCQYPbaKySZvn1jDU+krHQ5qZjZTOC77v5WlmN7AP/h7p/LZ3Ai3fHee2F79lmorIRRo0KC2Wsv2G03DVcWKaRcrlQuAH4O7JRUgCHA+YCSihSV1lZ4++2wPfII9O8fkktjY5h42b9/3BGKlJZcb3+19wDZXen8mScisfvgA3jppbCZhUUuGxtDh//w4bqKEempDpNK8pkmU9OKvmNm72dU6wdMAl7Ic2wiBeUO77wTtkcfDbP699orDF3effeQZKrU6yiSk87+lxlNSBgQrlIOAjZn1NkM/BW4PL+hiURr/Xp48cWwQeiPGTYsJJjddw9XNcOG6WpGpCMdJhV3/yHwQwAzewuY4u5zowhMJG6trfDuu2FLjSyrrQ0jy8aODZvmyIjsqMsX9+6+ZyEDEekNWlq298kA7LJLSC6jRoXbZUOGhCsckXKVy5Di4zqr4+5P9Cwckd5l5cqwPfdc2K+sDPNkhg+HXXcNr8OGwYABum0m5SGXbsgE7Y/+StHfaFLWWlth2bKwzU27UVxTE65ihgwJSSf1c0ODniMjpSWXpPKRLGW7AJ8Ajge+lJeIRErQpk2weHHY0lVVhVtoQ4eGK5qhQ8M2eLCSjfROufSpPN7OoVlmdj3wSeCBvEQlUia2bdu+EnO6yspwJTNs2PZkM2xYuLLRbTQpZvkahT8HuJ3w/HkR6aHW1uzJpk+fkGwGD96+bdkC69aFeTZKOBK3fCWVJqAtT+cSkXZs3bp9mHNKUxP84Ach4TQ0hESTek1t9fW6nSbRyGX013lZiquB8cBFwKx8BSUiudu6FVasCFumigoYNGj71tCw488anSb5ksuVyk3tlG8Gfg/8a1dOYmaDgV8BJwHvA5e7+60d1K8G5gG17j4yrXwGYYBAI/A5d28vPpGy19a2/REB2VRUhImd2baBA8OTNuvqQvLRFY90JJekkm3y4yZ3X56lvCM3AFuA4YRlX+aY2Vx3n99O/W8AK4DMRy/NJSSza3P8fBHJ0NYW+mXWreu4Xir51NWFZFNbG1Z6HjBg+2vq5/79lYDKUS6jvxb19MPMbABwOjDe3VuAp8xsNnAu8K0s9fcEzgG+BvwyI54bknU29TQuEemariaflH79dkw47SWg1KtWI+j9cu6oN7OPAEcDI4AlwNPu/lgX3z4OaHX3hWllcwm3sbL5MXAFPVxW38ymAdMAhgwZSlNToienky7o27dF7RyR3tLWbW1hmZuWlvbrmIWrm9RmFhJNernZ9v3010L3CbW0tJBIJAr7ISUgl476wcCdhEmQbcBqoCEcsgTwaXdv547th2qBtRlla4GBWT5vKlDl7veYWXNX48zG3WcAMwDGjm3yBQt6dDrpgqamBGrnaKitA7Pw6OjMraZmxy29rLo67Ke/VlVlT1CJRILm5ubIv1dvk8uVyo+Awwm3o+50961m1gc4A/gpYTXjczs5RwtQl1FWB6xPL0jeJrsOODmH+ESkjLmHlQs29fCGeEVFSC59+mx/TQ3Xvv327fvd2aqqtm+lOtoul6TySTJGarn7VuB3yauYq7twjoVAlZk1uvtrybIDgcxO+kZgDPCkhZavBurNbBlwlLu/nUPcIiJd1taWPTn17w8LFuTvc9ITTLbkky0Rpe9XVoatoiL7a3q99HOk6qTXz2eCyyWptAKvtXNsQfJ4h9x9g5nNAqab2cWE0V+nAhMzqr4MjErbnwj8BDgEeA8+HGpcARjQx8xqgC3urkmYIlL0tm0LWzFI9V2lEk3mtu++XT9XLgP+7gXObOfYZ4A/dPE8lxIeQbwCuA24xN3nm9kkM2sBcPdt7r4stQGrgLbkfip5PUzowJ9I6C/ZCHS6PL+IiOzIPSS4LVtg48YwmGLt2jCv6b33uj7aD3K7UvkjcL2ZzSF02C8nzDU5A9gf+Fczm7w9SH80e/C+CpiSpfxJdp6LkjqWAEZmlDXnELuIiEQgl6RyV/J1FPDxLMfvTr4a4bkrGnEuIlIAqX6S1O2p9J/T9zP7VDJva2X2w7S3DRrU9dhySSqT6fwhXSIiZaWqasfRYtk63jM72DM7zrN1zmd20qfXL+aRY7nMqE8UMA4RkYKqrAwz/NPnq2QbOpzaMo+98QZMnrxzfS1Fs6NcJj++CUx197lZjo0HZrv72HwGJyKSkpqg2NFkxo7Kejo35J//hOHD8/d9SlUut7/GAH3bOVYD7NHjaESk5FVVhSuG1FVD5s+pxSjTt379wvuk+OX6n6m9PpXDgDU9jEVEeql+/XZeMj+1YGTmApJ9+hR3n4D0TIdJxcy+Cnw1uevAH81sS0a1fsBgwuOERaREmIXl7VNb5jNWBgzY/qqrCEnp7J/Cm8AjyZ/PB54nOaM9zWbgFeDG/IYmIoVUXb09QdTX7/hkyEGDQpmWopdcdZhU3P1ewkx6kmtwTXf3tyKIS0R6oLp6x0RRXx+21BXHwIGhM1sk33IZUnxhIQMRka4zC0mioSH0V5x4Yvg59ez5fv3UbyHxyGVI8bc7qeLu/t0exiMiaSorYcgQGDoUhg0Lr0OHhuSRujWVSMCxx8YapsiHculeu6qDY6lRYUoqIt1UWwu77Qa77hq2YcNg8GD1a0jvksvtr53mjZpZA+E5K18nyyKRIrKjykqoqwtbfX1IHKkkUpt1OVWR3qVHAwHdfTXwGzPbBbgBPalRBAgd5aNGwejR4XZVqqN8wAD1dUhpy9fo8rno1peUsYEDQwJJbcOHa00oKU/5SiqfYOf5KyIla8AA2HPPsI0ZE/o+dAUiktvor5lZiquB8cAE4D/yFZRIsamthZEjtyeSoUOVRESy6enzVDYBi4D/AW7OV1AicerbF3bfPWwjRoStrk5JRKQrchn9NaaAcYjEpqIi9IOMGwd77x3mhag/RKR7cupTMbPdCMOHjycsIrkSSAA/cPdleY9OpEAGDIDGxrDttVdYdl1Eei6XPpVxwFPAIOAvwOvArsC/AueZ2SR3f60L5xkM/Ao4CXgfuNzdb+2gfjUwD6h195Fp5Qclz7Mv8Cpwkbu/2NXvI+XDLPSBpG5pjRwZJhnqdpZI/uVypXItsBY4wt3fThWa2R7Aw8njp3XhPDcAW4DhwEHAHDOb6+7z26n/DWAF8OHUsGSiuZfQl/NT4AvAvWbW6O6ZS/NLmenbN1x9jBoV+kN23TXMGxGRwsslqXwE+GJ6QgFw90VmdhXhl3uHzGwAcDow3t1bgKfMbDZwLvCtLPX3BM4Bvgb8Mu1QczL2/3F3B35kZv9GGEzwYA7fSUrEkCGhT2TcuJBMtLSJSDxySSrVwPp2jq1PHu/MOKDV3Remlc0l9NFk82PgCmBjRvn+wLxkQkmZlyzfKamY2TRgGsCQIUNpakp0IVTpib59Wwrezn36bH8MbSqJvP122MpJS0sLiUQi7jBKntq5a3JJKi8Cl5nZA+7eliq08KCVS5PHO1NLuIWWbi0wMLOimU0Fqtz9HjNr7u55ANx9BjADYOzYJl+wIPN0km9NTQkK0c7V1XDAAXDYYeG2lkAikaC5uTnuMEqe2rlrckkq04H7gFfN7PfAu4SO+k8DjcApXThHC1CXUVZHxhVQ8jbZdbS/lliXziOlY9gwOPzwkFD0cCmR4pXLPJUHzewTwNXA/waMMBnyBeAT7v5wF06zEKhKdqinRoodCGR20jcCY4Ank0+crAbqzWwZcFSy/tfNzNJugR1AGAQgJaKmBvbdFw46KMwj0WgtkeKX0zwVd38QeNDM+gMNwGp3/yCH928ws1nAdDO7mDD661RgYkbVl4FRafsTgZ8AhxDWGFsKtAJfNrOfA59P1ns0l+8jxaeqCpqaYMKEMBGxKl+r04lIJLr1v2wykXQ5mWS4FJhJGCa8ErjE3eeb2STgAXevdfdtwIeTKc1sFdCWNsGy1cymADcC/0WYpzJFw4l7p4oKGDs2JJJ99tHtLZHeLPK/A919FVke6OXuT5I2FyXjWAIYmVH2d+DQAoQoEUglkv33D1cm/fvHHZGI5INuLkhkKivDpMT99guJpF+/uCMSkXxTUpGCq62Fo4+GQw5RIhEpdUoqUjD19XDMMXDwwWGiooiUPiUVybshQ2DQIDjjDC2XIlJulFQkbwYPhhNOCHNLnnhCCUWkHCmpSI/V1MDxx8MRRyiRiJQ7JRXptoqKsHTK8cdrSLCIBEoq0i1NTfAv/xL6T0REUpRUJCeNjWFE15gxcUciIsVISUU6VVkZVgc++uiwWrCISHuUVKRdNTWhz+SII2Bg1ifViIjsSElFdtK3Lxx3XEgoera7iORCSUU+ZBZmv0+eHJZWERHJlZKKADBqFHz847D77nFHIiK9mZJKmRs4MAwNnjBBT1YUkZ5TUilTVVVw1FEwaZIeiiUi+aOkUoYmTAhrdA0aFHckIlJqlFTKyOjRcNJJMHJk53VFRLpDSaUMNDSEfpN991W/iYgUlpJKCauogOZmmDgx9KGIiBRaRdQfaGaDzeweM9tgZovM7Kx26n3FzN40s3VmttTMrjezqrTjE83sWTNbb2bzzOzY6L5F8Rs6FD7/+TCJUQlFRKISeVIBbgC2AMOBs4Gfmdn+Wer9ETjE3euA8cCBwJchJCZgNvB9YBBwHfBHM2sofPjF78gjYdo02G23uCMRkXITaVIxswHA6cCV7t7i7k8RksO5mXXd/Q13X5N6K9AG7J3cnwgsd/c73b3V3X8LvAecVvAvUcQGDoRzzgmTGPVMeBGJQ9Q3RsYBre6+MK1sLnB8tsrJW2M/BwYC7wNfTx1KbjtUJ1zRZDvPNGAawJAhQ2lqSnQz/OJVUwP19bB4cdji1tLSQiKRiDuMsqC2jobauWuiTiq1wNqMsrWEpLETd78VuNXMGoHzgOXJQ38FdjezzwJ3AWcBewFZnz/o7jOAGQBjxzb5ggXNPfsWRaS6Gk4+GQ48sLhGdiUSCZqbm+MOoyyoraOhdu6aqPtUWoC6jLI6YH1Hb3L314D5wE+T+yuBU4GvERLNx4A/A0XwN3p0Ro2CSy6Bgw4qroQiIuUr6iuVhUCVmTUmEwWEDvj5XXhvFeFqBAB3fxw4HCA5KuwN4L/zG25xqqgIz4WfNCn8LCJSLCL9leTuG4BZwHQzG2BmxxCuOG7JrGtmF5vZsOTP+wGXA4+kHT/YzPqYWR3wf4HF7v5QFN8jToMHw0UXhaSihCIixSaOX0uXAv2AFcBtwCXuPt/MJplZS1q9Y4CXzGwDcH9yuyLt+DcJnffvALsBU6MIPk6HHgpf/CKMGBF3JCIi2UU+Lc7dVwFTspQ/SejIT+1f2Ml5Ppv/6IqTGXzyk3DIIXFHIiLSMc21LnJmcOqpoTNeRKTYKakUMTOYMiUMFxYR6Q2UVIqUGZx2Wnj2iYhIb6GkUoQqKkJCGZ91fQARkeKlpFJkKirgU5+C/faLOxIRkdwpqRSR6upwhbLPPnFHIiLSPUoqRWLoUDjjjPAqItJbKakUgQkTwjyU6uq4IxER6RkllRhVVsLHPgaHHaYFIUWkNCipxKS+Ptzu0pIrIlJKlFRi0NgIU6dC/6xPfxER6b2UVCI2aRJMnqzbXSJSmpRUItTcHDYRkVKlpBKRj3wkPANFRKSU6TFPEVBCEZFyoSuVAps8GY47Lu4oRESioSuVAlJCEZFyo6RSIMcco4QiIuVHSaUADjgATjwx7ihERKKnpJJne+0VHv+reSgiUo4iTypmNtjM7jGzDWa2yMzOaqfeV8ybBoqBAAALiElEQVTsTTNbZ2ZLzex6M6tKO36QmT1pZmvNbLGZfTu6b5HdbruFpVcqK+OOREQkHnFcqdwAbAGGA2cDPzOz/bPU+yNwiLvXAeOBA4Evpx2/FXgCGAwcD1xiZv+rkIF3pKEBzj4b+vaNKwIRkfhFmlTMbABwOnClu7e4+1PAbODczLru/oa7r0m9FWgD9k6rMgb4nbu3uvsbwFNAtuRUcP37wznnQG1tHJ8uIlI8op6nMg5odfeFaWVzCVcaO0neGvs5MBB4H/h62uH/Ac4zsyuBscDRwHXtnGcaMA1gyJChNDUlevYtdjg37LILvPRS3k5ZElpaWkgkEnGHURbU1tFQO3dN1EmlFlibUbaWkDR24u63AreaWSNwHrA87fB9wG+AfwMqgenu/lw755kBzAAYO7bJFyxo7sFX2K6iAj7zGRg3Li+nKymJRIJmLXQWCbV1NNTOXRN1n0oLUJdRVges7+hN7v4aMB/4KYTOfuBBYDpQA4wCPmpml+Y74I588pNKKCIi6aJOKguBquSVR8qBhITRmSpgr+TPYwm30X7j7tvcfTFwO3ByXqPtwOTJcPDBUX2aiEjvEGlScfcNwCxgupkNMLNjgFOBWzLrmtnFZjYs+fN+wOXAI8nDC0OxnWVmFWa2K3AmoX+m4I49NjwXRUREdhTHkOJLgX7ACuA24BJ3n29mk8ysJa3eMcBLZrYBuD+5XQHg7uuA04CvAquBF4GXgWsKHfxRR8EJJ2hyo4hINpGvUuzuq4ApWcqfJHTkp/Yv7OQ8jwKH5z3ADhx+OHz0o0ooIiLt0TItXXTIIXDyyUooIiIdUVLpggMPDCO9lFBERDqmpNKJPffUApEiIl2lpNKB2lo4/fQwyVFERDqnX5ftMINPfUrreYmI5EJJpR0f+QiMGRN3FCIivYuSShZ7763JjSIi3aGkkqGuDk47TR3zIiLdoaSSpqIi9KP07x93JCIivZOSSpoTT4TRo+OOQkSk91JSSWpqgqOPjjsKEZHeTUkFGDQIpkxRP4qISE+VfVKprIRPfxr69Ys7EhGR3q/sk8pJJ8GIEXFHISJSGso6qUyYAEccEXcUIiKlo2yTyoQJ6kcREcm3yB/SVQwOPRROOUULRYqI5FvZJZWKCvjEJ3SFIiJSCObucccQKTNbDyyI6OPqgbURvb8rdTuq096xrpZn7g8B3u8knnwq1bbuSlmUbV2q7ZytXP+md9To7vWdnsndy2oDno/ws2ZE9f6u1O2oTnvHulqeZT+ydi7ltu5Kmf5N97yd22lX/ZvuxvnUq1BYf4zw/V2p21Gd9o51tbyn37WnSrWtu1oWlVJt52zl+jfdjfOV4+2v5939sLjjKHVq5+ioraOhdu6acrxSmRF3AGVC7RwdtXU01M5dUHZXKiIiUjjleKUiIiIFoqQiIiJ5o6QCmFm9mT1rZi1mNj7ueEqVmR1tZk+b2eNmdpuZ9Yk7plJlZsPN7K/Jtn7UzHaLO6ZSZmafNbP34o6jGCipBB8ApwB3xR1IiVsETHb344E3gVNjjqeUvQ8cm2zr3wAXxRxPyTKzCuBTwDtxx1IMlFQAd9/q7voro8Dcfam7b0zubgPa4oynlLl7q7un2ncgMD/OeErcWYQ/SPXvmRJLKmb2JTN73sw2m9lNGccGm9k9ZrbBzBaZ2Vkxhdnr9bSdzWxP4OPAfRGF3Gv1pK3N7CAz+xvwJeD/RRh2r9PddjazSuAM4PcRh1y0Sm1ByaXA1cBHgcxnOd4AbAGGAwcBc8xsrrvrL7jcdbudzawOuBk41923RBdyr9Xttnb3F4EjzewM4HLgi5FF3ft0t53PAe5w9zbTKrVAiV2puPssd/8DsDK93MwGAKcDV7p7i7s/BcwGzo0hzF6vu+1sZlXAbcBV7h7Vop69Wg/aum9a9bWEfkNpRw9+d+wHnGdmDwKNZvajKOMuRiWVVDowDmh194VpZXOB/VM7ZnY/cBLwSzO7INrwSkZn7fxZ4Ejg22aWMLMzow6whHTW1oeY2RNm9hjwFeD7UQdYIjpsZ3f/d3c/yd0/Brzm7l+OI8hiUmq3v9pTy87LOq8ldGAC4O4nRxpRaeqwnd39FuCWqIMqUZ219dPAcVEHVYI6/d2RonXBgnK5UmkB6jLK6oD1McRSytTO0VFbR0PtnKNySSoLgSoza0wrOxANs8w3tXN01NbRUDvnqKSSiplVmVkNUAlUmlmNmVW5+wZgFjDdzAaY2TGEiXe6FdMNaufoqK2joXbOoyifZFboDbgK8IztquSxwcAfgA3AP4Gz4o63t25qZ7V1qW1q5/xtWvpeRETypqRuf4mISLyUVEREJG+UVEREJG+UVEREJG+UVEREJG+UVEREJG+UVEREJG+UVEREJG+UVEQEADNrNjM3s+a4Y5HeS0lFRETyRklFRETyRklFeiUzuyp5q6bRzOaYWYuZLTKzb5tZRVq9C5L1xmR7f0aZm9nVZvb15Lk2JM89LLndYWZrzewdM/v3LsRYa2Y/NrN/mtlmM1tuZn82s33S6nzJzJ42s1VmtsbMnjGzUzLOMyYZ2xfN7HtmtszM1pvZb82sv5ntbWYPJdvgdTM7v522mmBmj5nZB2b2rplNT2+rDr7Hacm4PkjGeKeZjc6oc5aZ/T0Zw1oze8nMvtDZuaX0KKlIb3cP8CgwhbCS7HeA8zt8R8fOBSYDlwKXAZOA3yQ/Zx7heeX3A/9lZp09LfR64IxkTP8CfBF4ERiUVmcMcCPwaeBM4HngPjP7eJbzXQ7sTvh+307W/3kytjnA1GSMvzaz/bO8/w/AnwltdStwZfI87TKzLwJ3A68AnwK+AIwHHjezgck6xwK/BR5PnvvTwC8zvqeUi7iXSdamrTsb25cqvzCj/CXg4bT9C5L1xmR7f0aZk3woU1rZD5Ll/yetrApYAfy6kxhfBn6Qw3eqSJ77YeDetPIxyRgezag/K1l+TlpZA7AN+I8sbfWtjPf/kvAEw0HJ/eZkvebkfupRujMz3jcG2AJ8Jbn/b8CquP9NaCuOTVcq0tvNydh/GRidrWIX/cndt6Xt/yP5+lCqIHn8dWBUJ+d6DrjAzK4ws8PMrDKzgpkdamb3mdlyQjLYSriqacpyvgcy9rPFtpqQ8LLFdkfG/u2ExDG+nfiPJjw693fJh1hVmVkVsDj52celfc+G5O24T5iZrlDKmJKK9HarMvY3AzU9ON/qjP0tHZR39jmXAb8APkf4xbvCzK43s/4AZjYKeITwEKjLgInA4cCD7Zy7p7Etb2d/RDvxD0u+/pmQ7NK3CcAuAO7+OOGW1yjCrbj3kn1HB7RzXilhVXEHIFJgm5Kv1RnluxT6g929hdAPcrmZ7UHok/gvwi/9fwc+BtQDZ7j74tT7UkmnAIYDb2bsAyxpp/7K5OsFZH8m+/rUD+5+F3CXmdUSbqNdCzxoZiPdva0HMUsvoysVKXWLkq8f3uJJ3sI5Kcog3H2Ru/83oc8nFUsqeWxNi20ccEyBwjgjY/8zQAvhlmE2fyUkjr3d/fks24LMN7h7i7vfR7hC240IkrcUF12pSKl7DngD+H5y+OxmwsiuvoX+YDN7GphNSCQtwPHAgcDNySp/JvSj/MbM/pvwS/g7hOegF+IPvs8n2+A54KPAxYTnsK/JVtnd15nZN4AbzGwooU9nLeF22fFAwt1vNbPphKuex4ClwEjgy8CL7v5eAb6HFDFdqUhJS3aqnwq8A9wE3AD8KflzoT1BuDr4HWFAwaeAr7r7D5OxzQfOBvYgJJ9vAt9Kvq8QTiUMApgNnANcDXy3oze4+y+A/0UYOHALIbF8h/AH6YvJan8jjAi7ntC21xKGF5+ClB1z985riUivZWZXAf8B9MkY2SaSd7pSERGRvFFSERGRvNHtLxERyRtdqYiISN4oqYiISN4oqYiISN4oqYiISN4oqYiISN78f2AogfW/ADCkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot it!\n",
    "fig, ax = plt.subplots(1)\n",
    "\n",
    "# ax.plot(n_mc_list, mean_list, lw=2, label='MC', color='blue')\n",
    "ax.fill_between(n_mc_list, mean_list+std_list*2, mean_list-std_list*2, label='MC', facecolor='blue',alpha=0.5)\n",
    "\n",
    "ax.set_xlim(n_mc_list.min(), n_mc_list.max())\n",
    "ax.set_ylim(np.min(mean_list-std_list*2), np.max(mean_list+std_list*2) , AS.num_rs_mc)\n",
    "\n",
    "ax.set_xscale('log')\n",
    "\n",
    "ax.set_xlabel('num samples')\n",
    "ax.set_ylabel('output')\n",
    "\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Runing bootstrap for n_mc = 12\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c2680afb944c492cb713114de59ee5b8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Runing bootstrap for n_mc = 12\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3fd576e7775c42b3b7993faed0c6f1de",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Runing bootstrap for n_mc = 12\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "95f820712b604d8fa27c0b03bc515a1b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Runing bootstrap for n_mc = 12\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8ce8bf845394408895739ba649171ed4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Runing bootstrap for n_mc = 12\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9a10497ead0549728f2748a08aae5740",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "N_bootstrap = 500\n",
    "\n",
    "n_mc_list = np.linspace(12, AS.num_gradient_mc, num=10, dtype='int32')\n",
    "\n",
    "mean_list = np.empty_like(n_mc_list, dtype='float32')\n",
    "std_list = np.empty_like(n_mc_list, dtype='float32')\n",
    "\n",
    "mean_list_boostrap = np.zeros(N_bootstrap, dtype='float32')\n",
    "\n",
    "# re-evaluate the gradients\n",
    "# re-evaluate the response surface\n",
    "\n",
    "for index, n_mc in enumerate(n_mc_list):\n",
    "  print('Runing bootstrap for n_mc = {}'.format(n_mc))\n",
    "  for i_bootstrap in tqdm(range(N_bootstrap)):\n",
    "    AS.num_gradient_mc = n_mc\n",
    "    \n",
    "    gradient_index = np.random.choice(AS.num_rs_mc, AS.num_gradient_mc, replace=True)\n",
    "    \n",
    "    AS.x_samples_gradient = AS.x_samples_rs[gradient_index]\n",
    "    AS.x_samples_gradient_flat = AS.x_samples_gradient.reshape(AS.num_gradient_mc, -1)\n",
    "    \n",
    "    AS.evaluate_gradient()\n",
    "    AS.svd()\n",
    "    \n",
    "    AS.poly1d_gradient = np.poly1d(np.polyfit(AS.xv[gradient_index, 0], AS.y_samples_rs[gradient_index], 2))\n",
    "    \n",
    "    mean_list_boostrap[i_bootstrap] = np.mean(AS.poly1d_gradient(AS.xv[:, 0]))\n",
    "\n",
    "  mean_list[index] = np.mean(mean_list_boostrap)\n",
    "  std_list[index] = np.std(mean_list_boostrap)\n",
    "\n",
    "n_mc_list_rs = np.copy(n_mc_list)\n",
    "mean_list_rs = np.copy(mean_list)  \n",
    "std_list_rs = np.copy(std_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1000"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AS.num_gradient_mc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot it!\n",
    "fig, ax = plt.subplots(1)\n",
    "\n",
    "# ax.plot(n_mc_list, mean_list, lw=2, label='MC', color='blue')\n",
    "ax.fill_between(n_mc_list, mean_list+std_list*2, mean_list-std_list*2, label='AS+RS', facecolor='blue',alpha=0.5)\n",
    "\n",
    "ax.set_xlim(n_mc_list.min(), n_mc_list.max())\n",
    "ax.set_ylim(np.min(mean_list-std_list*2), np.max(mean_list+std_list*2) , AS.num_rs_mc)\n",
    "\n",
    "ax.set_xscale('log')\n",
    "\n",
    "ax.set_xlabel('num samples')\n",
    "ax.set_ylabel('output')\n",
    "\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot it!\n",
    "fig, ax = plt.subplots(1)\n",
    "\n",
    "# ax.plot(n_mc_list, mean_list, lw=2, label='MC', color='blue')\n",
    "ax.fill_between(n_mc_list_mc, mean_list_mc+std_list_mc*2, mean_list_mc-std_list_mc*2, label='MC', facecolor='blue',alpha=0.5)\n",
    "ax.fill_between(n_mc_list_rs, mean_list_rs+std_list_rs*2, mean_list_rs-std_list_rs*2, label='AS+RS', facecolor='red',alpha=0.5)\n",
    "\n",
    "ax.set_xlim(n_mc_list_mc.min(), n_mc_list_mc.max())\n",
    "ax.set_ylim(np.min(mean_list_mc-std_list_mc*2), np.max(mean_list_mc+std_list_mc*2) , AS.num_rs_mc)\n",
    "\n",
    "ax.set_xscale('log')\n",
    "\n",
    "ax.set_xlabel('num samples')\n",
    "ax.set_ylabel('output')\n",
    "\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}