L96 Statistics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import itertools\n",
    "from SparseARD import*\n",
    "from Oracle_Training import *\n",
    "import json\n",
    "np.random.seed(0)\n",
    "\n",
    "retrain = False\n",
    "\n",
    "noise_percent = 1\n",
    "tol = 1e-8 # tolerance for ARD algorithm\n",
    "verbose = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 200)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAD8CAYAAAB6tolUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eYxl2X3f9zt3eXu9elXV+zLTQ84MSZEUR+KIpqIIomlBUZRFNKAgFhxDiWSMg4SCnJFiKnKEWEoMxKscBElg2rJFBHYiRYsVEDQQheIiSiRlLsMhh82Z4fR0c3qtrr3e/u69J390adK/7+/Hfo8908Pi9PcDEJxz+3fvPfec31ludd9PhRijEEIIIYQQQggh326Sb3cFCCGEEEIIIYQQEb6gEkIIIYQQQgg5JPAFlRBCCCGEEELIoYAvqIQQQgghhBBCDgV8QSWEEEIIIYQQcijgCyohhBBCCCGEkEPBwi+oIYQ0hPDFEMKHD8oPhRA+G0J4PoTwmyGE2r2rJiGEEEIIIYSQ1zvfyt+g/pyInL+t/HdE5NdijI+IyLaI/MyrWTFCCCGEEEIIIfcXC72ghhDOiMi/JyL/9KAcROS9IvLbByEfEpH33YsKEkIIIYQQQgi5P8gWjPtHIvI3RGTpoLwmIjsxxuKgfFlETnsnhhCeEJEnRERCrfbO2tFjL/9ZDAvcGWMWOifOv84i95L51zGnOPcJ3nUwZoH6LXKdufUTkQDtkyeliWmnE1VuJVMTM6zqqtwvdHk40mURkXxfl5PCPlOV61qXuQmRmGHZXifAscTJi7LUP6NJxrbFkpm9v7q382MePObleqjuXL51Ip7khOCxBeqz0HhYIN0MXv0WuTf2TeL0J17HqV+s4OKF05+FLofChNgpZIF+qFI4JbUxi8xn39Z5EWK8eSlJdKJmiU1cPBadh5qUehCXxfyflyapvRfWsSy9zsGg+ePcjMdF+mqRflgEr6vgmJsnYX7M3GXOG/cL5L+5t5P/i81D3/o6t1CMk/84pyTQ6V7+45ohY5u3sHy6Q2+R3MEYt89xrfGG0bz2upu53rnuIjm50Px/tyySp+beC9RlgX2VOeUunykucHXcCy6yfzT38XIS7u3F2KRc4GaLjPOFFr75l/EnPWivFPbAmd0D12Gj4O0fnRnFHKmgPm5zQYy3FSxhYJeVHehm7cN1D8tehdw55s7tJ+Ln+/iFaxsxxqPOFRdi7gtqCOHfF5H1GOPnQwjv+bPDTqjf7jF+UEQ+KCLSOHM2nn3/ky//WQUvEPjSISISoSGwLCIiZmPovaxAlzuTeIBNldsJ0FEpbJjwz0VsYidOTAr3XijGGx7QM96gqmV64J1s7ZuYd/cuqPL3Ni+amH8zekiV/3jrYVX+4lf0n4uInPqormBzw774Dk7pF9vBCZtu4yO6PFuzbxn1tZG+V92+ae7utFS58WzDxLSvQg5CI8865hQp2rqMLy8iItlQl/OBjcFNcuW8rOMLfNl0YuCxvOvMu/etE6EM46jKnR841LDsXLgBud2w/ZnX9THzMioi076+WbppH7RxU59X37HVyUYwIS/QD+Me5MWyjcHzvLywc54Tg3PnIvMizoEiEqC/Elios7pduDtNvftea9nEXa3rsVc4u+YXttdUeeemM5DgtPbyyITgvNjv2zFcDvXikm3bxQbzAsej1w+mP701bIGXMswv5+eBgj9H9HLHjHOnPrgk4A9oUufe88a9d29/XoQfXNQWeBP3wD2y9wMt/AFlbnM5h3mm3dAP7/3wZWtPrxnJc20T030RNsAjW78SfxBbcza3MHcW9me+Urag7MR4eXA7bpMv0Od4XW/+N3OV88NkwZhFfmB/lz9IwVyZNweK2L2gtxnHvVetZtcw3I8V+MMOsS80XlPgvjNLnTrPacNZYScQrE8xszEV/CAxej8QxBd6b3xCW1ROfRaaC8ZwnvND6Qh7jlpXj/PTq3YT8ODStip3s7GJmcHkXjiLBP4g1lsLC5jMx6UdsPtTPbB3xnajt7+vj5X7eoEKzg/TAvSft5eo2jq/akt2kWg4++uv/PivXjIHvwUW+RvUHxCR/zCE8GMi0hCRrtz6G9VeCCE7+FvUMyJy9ZVUhBBCCCGEEELI/c3cf1MVY/xvYoxnYoznROQvicgfxhj/soh8TER+4iDsp0Tk9+9ZLQkhhBBCCCGEvO55Jb8H9QMi8mQI4ety65vUX391qkQIIYQQQggh5H5kUUmSiIjEGD8uIh8/+O8LIvKuV79KhBBCCCGEEELuR76lF1RyuFjE3FY6Ipm9sbZZDKY1E7Mz1THjozZV3lC7qcqNVH8k/Y63XjTn7D6sr3t1b8nEjPbho+2p/fDcfDvviBcm+/qj8tnEPsPRNS2IavzQlon5xlUtc2k9q6+7fMGKNFrX9Efk6cBRAcO/X5iuWLnL3oP6I/f+GXuZ8rT+eL/bG5qY1YYW27hGY7DCjWbWBDQY61zBNvWECShIqDninWZdt1erZturU9PP0IWyiEgLDC/11IoqZiAk2Jq0TMylnRVV3rva1fd5yeZkY1OXHaeCTCHdUaYlInZge7bFBMaIM5MbCV1q/8EMCkuqVD/X1JGebOU6B7ZSa8NJ8wrKts/Nc01sm7Yv6mP1LTtX4ThqOnIqdFfkfRuDoiKU/Ex79pyiA9Kfup0LUDrhyffQkhhnjhl2F0RPw/kLQOVIiKo69vncy0gC4pHUuTdK3+qbJkSa6yBY8cQ7KH1zuhzz3ZVTmfy38xk6YLYgb5OelYEcgzVj+Qc3TMz243pOubbeNTH5Vb2OtK6bEGlu6GdoDx3RWQX9mdm+KevhjuVZy55TwLToSa9wB+kKcKcg/RnZTjcGbWe6MNI+z5+D0lJvXoQ0MOMz2jwJkP+eQBDvNfFEl/BbArznRLmd48sR9AmNnf2PndshwJNV4TzkyY1gTxkmzlw10jHZwLGmw9BCyZqIFX5l1pFn5p2Z3VLKBDyyJcyvl67rPZ6IyMXnT+j77NqJEu/t/EKMxZgjrhOxY8S7VwfFyJhLzlxawLj2xjmKlDzjftryrJqvjFfyT3wJIYQQQgghhJBXDb6gEkIIIYQQQgg5FPAFlRBCCCGEEELIoYAvqIQQQgghhBBCDgV8QSWEEEIIIYQQcigIMXoqtHtD48zZePb9T75cNpZJz8SXoonPqS+YqryYkIFhynk1D4mO8WyLaCRN0+qOfy5ijamJZ3eDey8SkzrKvDzVaq/gxKCddXvHqkTLoe6M4LRX+1l9ndYNfa+tt9lzBIyzlWd9hTp7/eC1D4J9c2J5z8SsNbWGDU3EIiLt1Jocb2d76lhg97QFdt2xOGbXtKKuecNeu74L+ebIgAs0MjoWuyncfrbkWHzbMEYcI6mx/KHdcLrAz7y8MYy2aS8v8DTXbgsHPFMhWIQTxyrcaOo+r4GFdor6TxEZbOo8aF62Nsjmui6nE1u/WVs/+8Sxx5ZNOODZM7H7HOkrtpc1/9pzcO70TI+xAXPVkk3cY6t6PD6wtGNi+jOtHTx/+YSJaTyjG2PlOZu32VgfG6/YB8N2R4OkawTFtcd2ubVnOkMEx2N13Bqqu8t6rqqlNm8r6OQYHaOrY3a/Hc8Oj2uYt64UpX6w4ahuYqa7+li+bfsBDaCeMdXYKe2tpEIDtTPvoOU4QJ42Wnbuz1Jcq518gxivr7BNx4VNsD60IZrpRURkoBsjHTuWXDDpGiOos65guy9iXC6a863R0TGD47oiTt4KrDVolhaxZlh8bhEx82DRgrHXnG/iToe2MRpgqK7t2ltje6HV/VZ9dNmztWYDXU4dYzyOiSnYzYsVq4pNuzoROh2rzV1p6mNLjk2/gGS5umf3P3sv6oVt+Tnbpq113RdF0/bn7hthvTzr7NfgtHRDd0R9y143xUd38h+tuJVnHIfzPHMzbjG98YhrvDce0YRsTNxdZzLtQB54e+sFXhO9ffrF/+SXPh9jfHz+2T78G1RCCCGEEEIIIYcCvqASQgghhBBCCDkU8AWVEEIIIYQQQsihgC+ohBBCCCGEEEIOBZQk3R5zyCVJeJ3MkTOgJKmV2w/GW5n+AntzbCU/V84fV+Vjf2o/Is9G+v5X3wMxPXvv5Lr+ct/7OH2G39N732xD7hTL9svz5lEtFRntNUxMuqk/lndFEC24NsqDvL6CvEhz21eYO55sAykrW8EK5CSVI0GJc8QoIrbOiSenghwsCn3v2Z6VeNTX9cDO9+29UUDjzQWI11coT3HFNvl8aUcFkp+koUUC9Za1GLQad5ZpiYjs7OqxFq7anGxd0X2FoiwRkZjoGFe2ge6zBSRJBciXvHGVdPVztjuO0KeprR3durV4dHN93nJuhRy9XI/hTmrvNQQzxcevPWJitv7NMX2dyyZEJtprJoMHQTjniJ5iH2QbN6z0B4Um0xUTIrNTuk2T3BEg7ejnrG0694ImdAVDkAdmzHiiJ5B/lC174bKt65w2rYQlQyGZM7dPhzAnF3ag43q+yFqd1219cEiM+3r+QpmKiEg6BhGVt2bA+oSiIBGRCIImqTmdhTIqRx4UF5ChBVh/ajB/eXNXp67HWj217Yfyp8xJOFwzPMweyTknkQVi4P5eDNax8oRMc+o3LKwN57mtI6q8d9Ha7eqbOlmCbVIrtlm1c0G2oufTTtvOi9g3BewdZoWdP+YJ1ESsRC3PbP3aNZ1Pq42hicF96JWBFSldfvGoKnfP2wV96RsgPSwcAWRNP/t0CcSSHUcmhwKkReZFu5ybudLI20TMmA2uACzcsSxi5Uooz/JkWijhMlI/sTnpSumcOe6Fv/HzlCQRQgghhBBCCPnOhy+ohBBCCCGEEEIOBXxBJYQQQgghhBByKOALKiGEEEIIIYSQQwFfUAkhhBBCCCGEHApo8b095jvM4hsWMN3Vc6uJO9nSGtW3d6+YmOVMW9c+fP3tJmbj/3pAlU/84boqVy1ruhud6ajy1X9rvq615lhfx2/Vusq2Y7Ebf1Vb9Jo37HVmujoyWXVsi5CnaHFMZtamZqy0nhEXLWyOQK9szs//iLmSOWMEzHHBsYSmNX2sVrO50wHb43JD90PDMT1ujLRO9tqlNRPTelE3GNroREQK6CvX4gvHPIsv9qdr8UWTHVqYHeNmgu3XsG2Btkw0LYqIDCZ63OzdRB2vSPOyVgp6uZ1O9XNNlm2CjU7omNmqrnPSsPVDPMuqTJM7l0UkcUyF5trJ/L5Cu2532RojJ1PdXrOLHRPTunzn8YhtJSJSHINEdYZedkP3p2cQru3BiU7TFHU46DQ7micLzyoJx6o6/rkzHsA4G9FkLmLs5p6JGA2gxkArIuUIrN9bdhDne2DSdcb5rD1/b1OD6zQ24D5Dew00gI7tdCazzgIWX5zL5w8HCaUNCrD+4PokIoLy68SxxyJ2DbMxZg5ewBTuSnPx2ALmX+86Jg+87SKu39AW3rpStCD/j9n9xtseuKrKb+qum5hE5pv6KxjYO7OmiTm/rX/LwpWLR0xM46puDDS4ev2JplrP6Gpi0EYtIhUYvbOeba+jvb4qn2jbjV4315UunM65sKsHoLe/aL+g26J9HfZDTreMezpPRsdtzPSYTp7cec5GXa8RnvUY9/K4jxcRSWEt9Pb/+B4xLfRzb+/b39hR3ND5VduybZyC5Nv97QjOsed/+cl7a/ENITRCCH8aQvhSCOGZEMKvHBz/jRDCiyGEpw7+99jdVoIQQgghhBBCCFngtw3KRETeG2PshxByEflUCOFfH/zZfx1j/O17Vz1CCCGEEEIIIfcLc19Q461/A/xnfxefH/zvtft3wYQQQgghhBBC7gsWkiSFENIQwlMisi4ifxBj/OzBH/3tEMLTIYRfCyHU73AJQgghhBBCCCHkjnxLkqQQQk9Efk9EflZENkXkuojUROSDIvJCjPFXnXOeEJEnRESy3so7z33gl1/+M0qSbrvOApIkjPEoK/1gVWVNAviR9qnOnon5oSPPqfIjdWtheWm6qsp//4/+XVVuXbQd2rypy0c/s2liJieXVPnFn7CdVVvRH883/8hKT44+pQU+VWavMz6iv/jfPWdjBufgo3YQAMi+85w39HVqtollBt+rT1fmj0VPyIQsJAZyREpmPDoiFJTmNFr66/leW7e5iMixVt8cQ85f1waC5Cu2PxuQKii4EhGZgU/InVOMkMOpEIxZI/rwhCE4VzlCnwhypdQRKaFcCSULIiI4be/tWPlBdlnbcFpXTYjU9u+cc1VuGwdFPCjNELE56MbAjzQ9IQe2c9FxcvKYngsePWnlJN/be0mVj+Tzc3IGNy+dREnhHxMlnm0DuDJZMcc+ee2Nqrz7ZSs96ehHcPN/eErf38tBA+atMx4CCK2w7J7nrGHYhEl3akKOH9GTZe6IxL5xWbdP9ykr5Otd0OMo6zvCLxjXI1gPtt/srBkP6TonjvSkHOrruHKjNsi9elbudayj87RXG5uYRqqvUzl5Oix0++xM9NzQn9i/X0BZ22Rs2xiFVjKy7WWe3ZUbzV+fcGglIzth1Lb1xbOBdy9dxrnJEyWGk3pdW16y69xooi80umbldrVtXWdPBojHMpsWUt/VdfTmcVyzph3dNtOu7YgCquxK1kDcWDbtnOdK1IAwQ6OhE4Tp5Ej7chhH3no5K3S7j7b1g9Wv2wWqDvuNzA49k8u4NnrHcN0TccRT3r9tXUCihoK7alW3xeoRK6I6Cfv/Vmbbz7zDOOscxoiI/Mt3/9N7K0m6nRjjjoh8XER+NMZ4Ld5iIiL/XETe9U3O+WCM8fEY4+Np2w5YQgghhBBCCCFEZDGL79GDvzmVEEJTRH5YRL4WQjh5cCyIyPtE5Cv3sqKEEEIIIYQQQl7fLGLxPSkiHwohpHLrhfa3YowfDiH8YQjhqNz6i+enROQ/v4f1JIQQQgghhBDyOmcRi+/TIvI9zvH33pMaEUIIIYQQQgi5L/mWJEmvlMaZs/Hs+598uUxJ0m3XgXtnqf0IOQ3zYyKYKlB0ICLS39cfiMfSNkYEGU+Y2pjv/x4tUnrPyrOq/I3JmjnnzU1tavl7X/sRE3Ps7+mvyJNPf9nEJC2Qwpw5aWLWf1CLNLbebtureUpLKM6s7JiY1boWIuAH4hXaCETkyqCrypefPW5i1r6kz2tuWgHArKVjxquO2KAJB5yP6fFD/cI6daSqzx9r5lFzaFNHjpC39Ef3R5atoAZFSlf3uyZm6zyIUV6wD5qAc6V0hARGWuCIDfAYSgu8ucrGOO0Hx3AOFBGJNX0s1G1eNNr6QVc61qTRynXMpLSV3h7oRBju6waLKEERkTDRSRAccVcy1cfSiQmx8g8nb8dHoL2O2wthW1SO5WfSh0ToW8FKAnNcAH+V44Aw46Gq2aCqpfuvvmyfYXVJN8assnPK5otartQ770hiQJZS1D3jkS6iiMoTfWBue20x0247GZ+zAqTTp7ZUuVe3shmU/DRSKxLrZLoNNybWb3FhQ68/4y2cKEXyTf1g7Stwn2t2PgulfvhZx/bV6Ih+huEJEyKzNT0vooBORCSBORj3KCJ2P4F7EhErRsRyzRFR1aHdl2o2b4809Ly9VrPzUBMmZdzHiIgkd/GbC/vO5P7l3VOqfP7ZMyamcUX3+fiors/yg3YPMIJ9VPmizbf2Vd3nngAJxWaebA/nxdqubZtsrI9NejYH994A553Wlp9G045PlGoWha1gMdHH4sAKhrIdHeM56QKknCcPKlownzmSPGnrPK137HO1G/pYq6bLmZOTASa50tnnjab62ftj+xDjoc6d2Lftlfb1tbOBI1XD9nLm6VkXJHlLMHd6H3XCep7u2z736oNgX4mIXPiFn3/tJEmEEEIIIYQQQsi9gi+ohBBCCCGEEEIOBXxBJYQQQgghhBByKOALKiGEEEIIIYSQQwFfUAkhhBBCCCGEHApo8b095hBZfBeJaWTWbojmTry3iMjV3WVVHn2tZ2Ka13W5/wbHDHhS6+ayL2pF3dEvWo3d9Xdrg9mP/QefNTFHQPn2O5ceMzH9L2pD45Gnbf0aW/r+ZcPaySo4VDRtYvRP6WP9c/pe6XFrojyxuqfKaw1rN7y0q62cuxdsP7Suzv8Z0ngNTHerNi8CStiGnqlN3wuNuCKO8RPMbWgsFRFJ2rof2h3HBtnRfX6qvWdipmCh/dK1UyamvKjtinnfsc9Bqrh2VjwNusER+lmj6wLzmTcPYRu7NuUa2PocG7CZv5wH9eYrHWAPLWQWxXs7t5lua5t447JtsOZNqI5NL2MzLK0kUSo45hkQcf3B/nTkwOYcz+JbNqB90KwoIssrA1V+YNmaRNGY+sLeERNz+WltM1960bF+g8Eb7bv5vjlFanBsZG8t0zfrebC3bOe8rU29RuSXrfUShZrT03YieuScXqAe6d40Mbszbe29PlwyMUgGufzStp2TJ1/XhvHuBdvG7es6UZPC5gXaf8crdjKY6KVaCiuPlbKBeyQbY5g/NZjreMZxPDGU843eoZhvBPXmM2OZ79r9xVLPrsXI8Ou6UVee0X/e2LbzWdnQdZ6256/L3hyD42Z8zLHS4vzgtGkK9unWVRtT37vz3F46hu8ChmPZMCFmDs7tMJcEHsH8pgERmejtj0xXbFtUTbiZJyXH9nHaC3MOz/HWFbOe5846DPVLG3Zuz2owFzhrbgm/SWPWt8mTbuk+r2/a58z1MmL6wVuHzTrn7FtmMO9M7S9ZkLJt++/iz/4CLb6EEEIIIYQQQr7z4QsqIYQQQgghhJBDAV9QCSGEEEIIIYQcCviCSgghhBBCCCHkUEBJ0u0xh0iSFDyhCXwTjeeIiLRrWihxur1rYs42t1X5xYE1Xjz1yUdV+dQnrZAg76OESHfg1R+0H3qf+oHLqtzK7HVvDrVIo3TsJAl87b03sl/zT6/qL7sbNxYQDnnSglWQdKD/xalfa2msyo8c2TAx59qbqpw7/TkE48v2tGVints+qsqbF1ZNTPsb+tnxY3oRK5JBeYqIyHRZP3zZAllP3T5DyBaQ6kBMre4IMJparuQJwDZ2dZ8X247lAYUJTpdbqYKNMedUkAfe1Apj3xPvJHDv1DqlJIFj3tyJYgpPYoBzMMoj0pEdM9he3nyLgqhqzfbnudNabPNQd9PEIOsjm5QXNrUwbXylY2JaV/Rz1Oy0aJ4LBSGeGGLahfHQcZKpoZPHrEVi1w1v7anVtPFiuTU2MdNSL4bbL6yYmOZ13RaTlfmStQwEHb1nnTlvQz+nl9tFUx8cHnekdCChizj/ikgs9HnZTWvGal7X98ocf84U0ml4BkSEZ7SYSkTkeFcbo/LUTg5FpetXOla1EmIqp8FwbSlxjhGRGZj+Ki+m0DHTiW6vcteu1fV1fU7DGZ55H/LWSf8y1/XxRE/YD9OeM6es6jlkec32zduOannWu5ZfNDFnc/0gT48eUOX//Zl3mXOWP67X3e5FO5/tvgHW6nfZvH3o7Lq+Ts1O7mOQAXp5gWvfpLALwI093ajjy3pe7Lxoc7J9Q3dgfdvmdqhgzmtaK9cUBGClI4yqamF+DIq63OvosicwNGsWLtXO+mnFek5O5vPfTyKclzQduSmII3stO1l1nFxB+lNtudoZ6k3AcN9K6WRbN2Bt2zYgrpeJTX+Xp//nJylJIoQQQgghhBDynQ9fUAkhhBBCCCGEHAr4gkoIIYQQQggh5FDAF1RCCCGEEEIIIYcCvqASQgghhBBCCDkU0OJ7e8xrZPHNUsdiaiy+tn7YVVXlGDbhXu26Nck92tP2zB9eecbEvLGmbXNfnZw2Mf/r139I1+cj2qZ59KmhOado607eP2Pti5MV/fCePXOG9swla0YLYFiTvrXNNcBUWN82IcYANwML4WTN9nmxBCd5ttYxmCj3bRDa0xrb9l6tm9q0l/WteS+mYMxrOlbJGpj36o4NsnlnI+PMClSNodGzyca2rnPSsP1Zb2h1XNMx/U5mOr8GO00TI2Pd59Gxqta6etzkua5P6o1zGMNeTAYxnokYmTnjPIXzGrlj2U6gTZ17zcD6Opxpo9/WnrVGV9d1m+Y7tn5FR9+rOmYNhFld168s7HUiWIXR9iwikubzn3M6ATPmntVBZru6LTIwXS9iLox2irFWSWftQetxWbfPgIbgpVP7Jua9Z59T5e/tXDIxOSipx6CrHDqqzBIq/dJkzcR88tobVXndsYkvXdAN1Llqx1421Meq3JmrYG4q7TIis46OGR2zMZMTelwnHehkZ3hWM+jkse30MAULeGHnUkxTzz5a1aB9WnZuz1q6zu2WHWvLTW18XqlrS2g3t0ZoHEf9whpA98AaivOJh/cbChabB/W1xzPb6WOY/4vSNmqe6TY80tYD/VRrz5zzUr+ny188ZWKOfV4/Q33H9tWsAzblJWeNbcP+Z9mEWPP2sr1X3tV50GnrctNZM9ACvuPM//Ganv8bG/YZ0Mg7Oerk7ZrOueisc9VNnV/tl2xMYwtN9CbE2MMLMLR78zbimvKheWZLjum3CWO44ZiRF3iPwD2k946QwHUyyHXMfRH/fQQZT/XDj0fOGjG0DfSNv/oBWnwJIYQQQgghhHznM/cFNYTQCCH8aQjhSyGEZ0IIv3Jw/KEQwmdDCM+HEH4zhOD8hiJCCCGEEEIIIWQxFvkb1ImIvDfG+A4ReUxEfjSE8G4R+Tsi8msxxkdEZFtEfubeVZMQQgghhBBCyOuduS+o8Rb9g2J+8L8oIu8Vkd8+OP4hEXnfPakhIYQQQgghhJD7goUkSSGEVEQ+LyIPi8j/IiJ/T0Q+E2N8+ODPz4rIv44xvs059wkReUJEJOutvPPcB3755T+jJOm2GLi3GwPGhtKxKgzG+l9azya2UfFD6mbDfiy/BuKAcWGFBO2a/ui+hI/cUbgiYj/CLx2JQVnpr7892Qx+7N2qWRkUiiA6uZVHjMGu8ezNoyZmdl5bCpZfAEHN1vyPzMuGIxLIna/cgaSAfJs4OQmHph3bpqOj+l5j6zgxYqeq6XzMn4PkB8tuboPEwOnPFPLf+5i/CaIizD8RkVqqz7u6Zw1bu5fBOuF1A9Qx2wLJjiOxwXOkmt+/Vd22V+jq8dho2dyejHXeVjt2rJyeZJ4AACAASURBVAUQDEUUrohIuqyvfWxFi3eONMEUJCIbI23GuvE1O2ZWv6Tv3dy0/VnWdZ6OV23eDo/r8vikvU73lJaaPLSyZetTt8+BVDCfJmBHQzGViEgZ9XP2Zw0Tc22oLWHXNnomptrQMhAUs4mIVMs6/+tLNv8j1Ge6ZevTuKxzuQHNhRI4EZHBAzDfnuqbmGNL+lg7t3lrxDvOGrE90uaR/b59hnJXn5dv2nWusaHLua2yaWfcg3jyJYxx5Ua4J3GugyKZ0jqIzDwTHUmY2e84+xaMCTAX1Nu2r9oNewzZH+i+mW3bh8i2dYN5/ZDArbz8x3kaxWIiTv85H54VbdhHwbjqHrFzxaNHtFjyeMOKlLanOm/Pbx63MTf0XJDuOQ8B64i31tj+dNZLmB9qsH4OR7aviqlO3MTJpQqFX9u2kRvrIIAcmRCZgTxx7IiUkhWdGJ4YqNzTgyvbsW2aggMMpUhGaikiZUfXx7xDiIjAfCsTR1wKx5KJfYgE2tTN/zm3FrHPhfNFheJQEYlwzNsnmHchb2vj1PnSf/aL916SFGMsY4yPicgZEXmXiLzFC/sm534wxvh4jPHxtO2sfIQQQgghhBBCiHyLFt8Y446IfFxE3i0ivRDCn/2o4oyIXH11q0YIIYQQQggh5H5iEYvv0RBC7+C/myLywyJyXkQ+JiI/cRD2UyLy+/eqkoQQQgghhBBCXv84/wDecFJEPnTwHWoiIr8VY/xwCOGrIvJ/hhD+BxH5ooj8+j2sJyGEEEIIIYSQ1zlzX1BjjE+LyPc4xy/Ire9RCSGEEEIIIYSQV8xCFt9Xi8aZs/Hs+598uXy/WnzRaiqymMUUDaX1tDAxk1I34o1NazGVdbAiOnbKdJDMC5G1x2+o8qM9bbr79DfO2frtgq3S60+nfQyoMHNOiWAx9cxjaVO34ZEVqxg8vbR7x6rsTKxlcjDVzzme2eSeFTpxCyiLiBQTfSwOrA4y29UxNSsYlNSx6CFVfueyiGOJS+785yIiVQ1Mco5ZsWrA+Gza3K61td12qTU2MUda2sBYS+w4eualk/o6n2mamOaWrvPOwzp5Zh2bcNloAWsvWiab8w2N2cBOVk099CSxIm4ZGQOubdME2hmtjXXH5Hm8q02/p9o24foznf/rw46JQaM3WspFrLURzxERGQz1vaa7jkl0T5+XDWxfBUgVzG1vPGD/VW3HfN26cxuLiJT7+uL1646VdlOXUyvxNSZF1wYM428EEubxcfsMguNz6tgqp/Pz31hoHSutOebYM7O6btNazeY25s5sZnNnNgY7N5S9Z0I7dkCrqYgkUJ1kNj/fgn0Es+665k60CntzMMQYm6dzXZxTsqENykB4681D8+pyqz53LossaE/G69RsflVgKY0Q49pa8dGdDRHuZdxxPoN91Y590PoWWF8dmXKhhcEy7TnmVThU39TXzfU0fuu64DEdOdb0pdP6xHOuNX2o7+WswxkMgDoOGudY5Wzccc87cpJnUOhJb2eq1/wru2D2F5FdMC57c3J9W5czuyUxe1N3fEKVXXs45vZdjCNvPJh3LM/0C3Oyu293uPTTH7j3Fl9CCCGEEEIIIeRewxdUQgghhBBCCCGHAr6gEkIIIYQQQgg5FPAFlRBCCCGEEELIoWCRXzNDXmU8L1VZwc8KHGECSkQajmDiTEcLfU51rMDkC9VZVV7+tJX8oOxg6/vtl/qbe/qL+s/9obayHHvOfhg/7ernHB614gWUeCzyUbknSUIJxWzZCQKBidc3L+6sqnJ/oCUsxdh+0R5BIoJiDRdPDoV97MiDSpAHDU84P3ca62Np3zYqimNc4QV+8I8iGU90hiIl7znxujP7DFMQRO07l8lBNvaW1esm5r/4nk+o8h8/+LCJefYj+tiZP9TSB08+U7T0w5cN+wxVBm1cOlK1CZpubO4MTup7bb7D1ic5CXUeWDtVekWPfRQmjE44YpSenmOWXDOE5qWiZ45t3dACt+ymHUco/PJkXy0Y5y0bIji9ehKWCtxKOA+VbdvpZQdEH11rLjrW01KRN3Y3TcwDTS0aaTiDb32qpR1f3z9qYi5srqny6HrbxLSu6LHfgiHSuezkba7PKR3RWQntN1uyMQVIpTz5HpqAjNBHRIqGzv9Z3ZHEZPOlhwnMr9my7j+UIHrHUHAoYkWIKEoUEQlwzOwBRGQy1c85GdqGjyDYQmmeiJW55Pu6jfOh0zaQgqHy5EH6OoX1k0nR0jGYJyIiJXjqvBhcN0LlyKlgm5L3vRjdPkaO5or+oH52yyRFW9evWHJEf0u6gtkZO1+MV/TN0qv2Zp2XdHntGVsf7L9JV9evf9a2zeisPqnWtfu+CQgfn7l80sRU0MaeVA2PeSKxBCRlnhgO5WKpI5XCvMhgHWn2bf4vjfCYzf+iCfLEtn2GGUzBKLgSsfOpJ4w1UllvDZsjoXMFsjU9AHJnj9mo67zwJK6lMx5fKfwbVEIIIYQQQgghhwK+oBJCCCGEEEIIORTwBZUQQgghhBBCyKGAL6iEEEIIIYQQQg4FfEElhBBCCCGEEHIoCNHTlt4jGmfOxrPvf/LlsrFSeQbQdL6FSkAY5pqq8JgXk8y3/gUwyaHRD/9cxFr+HCmnsfxVcb4Rq55b29Zaa6DK37VsLabnGtoi+XtXrQK0/9unVPnoZ3ZMDFr9ZqtaT3b9XdY+N31c129laWhiNrY6qpxdaJqY9lWsjAmRvpYVy+y0o4AD0nVrSWxs6IujNQ4ttSLWBOha2cDC5hlwq1q8Y1lEJKLBchGjpWOVrGZQ6b6tUL4PNmBoUq8t0NzptYWxMDtjr2rAM7Rt/neWtFH21PKuiXl775oq/7nO1+294EH+7vM/osrDTx8x56CJr2zaZygbcAwVtCISwfhcX7cN1r2gr9O+5mi/gemy1VPun9bP2T8HbbxslYgR8iTdsvVrrMMz2OlD0hnkhZM7ZQ2NrjYGx40Xg9bNWcfpm67ui6QD5sKGbeMM5n/P1oosMrd7Zlg0J9ZSmzuNbDY3pgBb7M2+nm93b+qyiEi+oRsVzaciIrNVqN+ao1yGR5/sWl1rvqk7NPNMrDgHL2B692KwK/C5PFs3PoOXt+bejnEf7cTRsRV78yBi7LbeOZ45XVXGyUk00U/sg6ZjsJJP5veVZ9zHvyrxxzBYj+vOhWrQYbmzFsKxNNednmXzx55nMcWxVsvsfNHO9Xzaze2eZKWm90i5k4TrE63I/vKNEyZmckFb0jsXdd80t5z9bQUm4rrtz7I+f072xto9A5rH5JuIpFPYt0Oz43OLiBSN+YZeXFfc32KwwPpk9oveXx1CFb05GJcfvG7p7B8rGFe4zxKx+0dvPknrtkIX/tJ/+/kY4+O2povBv0ElhBBCCCGEEHIo4AsqIYQQQgghhJBDAV9QCSGEEEIIIYQcCviCSgghhBBCCCHkUPCaS5LO/Oz/L0lCcYArMbgbSZLzkXvAV/FvoyQJxRoiIgmc53XLrNAPWpS2wVDSsdKxEqJ3rF1R5T+//DUTk4MJ6J9f/bdNzFe+dE6Ve+d1Izdv2uccr+mYnUftg9Ye7Ktyoz4zMTsg8mg9Z+VG3UvQ7iPnA/FMf/g+7doP4Udr+ljRhoD5zhNfpAEfz3siDfzAHoUY7nnOdVAekeT2g/Y0ny9Smk2gQtu63esb9kHzPVsdc90lKHcXEAw5AgxsZ1cYAse8tlha0oKXlaYuv3jxmDmne17bD4an7L07j26r8pmulTh1QKSxlI9NzKTU/fCF62dMzPT8sir3nrP1QbmSkWQ0bX/isSq3AwBzG8UaIiKzFsRYF5qU4NBBEZWImPHnSjLA15M5/p7ETjOKahERzwISJ3+dm38vuxY6MbimeuKdBsigamgGsqdUff1g9WuOGEu796RomRDpv1F3zsmzWybm9JIeE43UdkwCxpDKqXQBjTh1GnVY6OfaHulKb+/bh5ht66RM9511eKbr487/ICOpujZxe0f1Wvim1Zsm5kRDt1fdGQATSMKtqV7Etib2Ofdn+jmnhe1zFH55kjDc/3gCsDTMl43djYAM50kRkf2xfq7hSJdnI2cQoyBqZvMteKIpwKzfnkwI1uq06ciW2npN6LXshFZP9XkjyPW9kZVYDkd6gi2ntoIRntPb8yYwDyXOntc7bx7engRb3btqMdV5UAxAxLbjCSHhOR1hGkqSpstO0JLuh6zhzGe4J/EklpDbxcT2TbWL+zEdU9fbj1v3Ag+iN2+PjoGg7Kh9hnbP5uBX3/crlCQRQgghhBBCCPnOZ+4LagjhbAjhYyGE8yGEZ0IIP3dw/G+FEK6EEJ46+N+P3fvqEkIIIYQQQgh5veL9RkKkEJGfjzF+IYSwJCKfDyH8wcGf/VqM8e/fu+oRQgghhBBCCLlfmPuCGmO8JiLXDv57P4RwXkRO3+uKEUIIIYQQQgi5v/iWJEkhhHMi8kkReZuIPCki/6mI7InI5+TW37KaT3BDCE+IyBMiIllv5Z0P/uIvv/xnRujgvC7jR+WeSAalJ55IychSnH/cfDeSJJQb3a0kKc+0qKKWWnELfiTdH1tjyKivP/jHD9pFRFL4CB8/uBcR6Tb1saXaxMTUEpBtwIfdmfNVeQ0+3EfZhXcdT/rQTPWX3ZVjobjQX1Pl81dOmJjsBW1mWbpkQqS+p5+jaIKIoecIOtA/4PgTsMooUxGxY8KLqXIUo3gx+hli3RF11XV/5o6coVnX7Y65PJrYmw829Vf3jas2BgUr+OG+iJXCmDYW+4F/Vbcx2F5FyxF1HddysbeevK7KncyOh0/96VtU+YH/x8n/LX3edNmO4eEJ3cn904646wHdN71T1kS11hqYY8hgqhsIBSLjka1fOdL1CwMra0iHOrlTO8UYmVH0xgjmv5PbZUu3c9m2c2eCYiBHpGfuXeoKVYWzaJT6GJ4j4ss1bNACMTiXe9fFGE+qBvN/vqQH25FlmzdHW1rW40k8ikq3hbeGoZTo5qBjYvb6emAXY9vpsZjfYAH2CinKoESk3tDP3m3q8dnIrAykgIl7q2+tIsMdMH6NHRtOU9ensWTnFBQE7vftpFdu6TGb7dl71WB6QElYcMRiuEfyRE9mDfPEXdB97jqHostF7CjOVhC2JJLYJjVzEZ7jSsIWmIfwmBezkAwN1idcr25dZ75sCWPsHthrQNjfOvNkms0XgyKVMy8WIFyMU9vpYZZAeRERlXMQhYqePA4Z6QtljgwNx1Xq5JsZI47or4DpYtZ1BJAdWMMadtCijCrgcuC8ezRhjvEkZihonczsIC6d9fGF//iXXxtJUgihIyK/IyJ/Pca4JyL/m4i8UUQek1t/w/oPvPNijB+MMT4eY3w8baP+lBBCCCGEEEIIucVCL6ghhFxuvZz+ixjj74qIxBhvxBjLGGMlIv9ERN5176pJCCGEEEIIIeT1ziIW3yAivy4i52OM//C24ydvC/uLIvKVV796hBBCCCGEEELuFxax+P6AiPwVEflyCOGpg2O/JCI/GUJ4TG59CXBRRP7aPakhIYQQQgghhJD7gkUsvp8SX+HwkVe/OoQQQgghhBBC7lcW+RvUV5fbjV/42uuYAReyG34bCV6dgQh6ytIx6wawQWaOSasNJl0si4isw6UHV60lMW7pe/XFmgH78FjBChAlAZNiYoWH9t746M4/MkdrXVW3bYzHQtvevNHWhsa1Xt/EpN+3r8rb39U0MVtXtNyr/ZJ+iLpxV4ukYEmcdm1MsaTLvsX6zmURsWPEG0Zo0XOMeWh8npY2Bk1traZu49WOtt+KiDy4qhuoetjm//ZImzB3+rYfZpP50xUaBbPcJm63pcdNp27HEdql0T56prljzvmb/86/UuVnftD+Nq4/vv4GVd6y8l2JlbbzoV1cxE7ce3u2vXY2YOw7fT7P+rqIET12HdNpV1935tgXA5hNcT4RcSzuDdsWoaXbq+7YDdFMuMi8fTeknsUdzImeJTGF+uSOARfrXFa2PyeFzoyhY3ovYVzjmNnat1banQHYzptWy7w71DHV1+3a07wJBxwTceNujKmOGbOENaKq2fljWNcXGubw7I7F1DNhIhHmydwx69Yv6WOZY8wOMLQ6Tlugzdy1RkNalmA3rxx/Ja5HnmUb1x7X9AuP7tUP9xeu9RvawjO9oxncuxfmygzSFNtGxLOJe2bdO5dvHdTFZOr8loWhPpY5+0WsT9F09kg13MSB1XpsOyvf0/fK902I6YdFLOVeXuAYLu0SZtrL24eWkP+zrmMYhz7N4DcUtFp2D9CG31hQyzzVtaY/scmzu6fnlHLf0ztD0Vt3cY2o2/o8AHutIw1tZH+p3zPnvHTxqCo3L9t5MoNtXebkdnB+q8IrZWGLLyGEEEIIIYQQci/hCyohhBBCCCGEkEMBX1AJIYQQQgghhBwK+IJKCCGEEEIIIeRQ8JpKkmLQH9rjR/feR/jm2AIf6t9L5sk1FpFvVI7cYuYaCDQo1zjWstKfE8f1V+3ns+MmZvjUqirn9jJGHJCNbIzpPxT6eNm1SF+BFABFESIiCcTEqRVMjPogwKjZr/BTEKrUHMFK9yEtxJmc1tftb9rrNq7qh69vmhDJ4cNzlF2IiBQgr6jcNoW28MZRpg9WmSN5gG/348SRJA311/H7qW73vdTaNlCqgyIejziz9077+uGTiX1QlFdMO9aqsLmnRQY7u54UAIRRR3Ve7D5o+7xf6raYOp2VwPxQFY6sakvXL99z2gIkImG+v8GVdqCoAmO8fIsoMfOEIXDM5ICISA52jcJWsLYD8ri+0xYztLnYW5m5yWsLvAxKT5y2wHNmznUXmNqNFMO9F7azM4bN2EI5moikQxD47MP84fzYetbW192eLZmY9lVdbt209pSioe816dn6zWBoFdbZZAQw0ZvPzLxjr2POgVNQcCUi0uvpiftNq+sm5i2d66p8xLHN9MHu8sW9sybmMy88pMrtL9lFYukyzHlt26bjNV1GsQyWRUTKJojFHFmPNwcj2Dd4XRERaep5OjhyKiNk8uRBEz0A0307IDHfU/DjuGIlnD/uYkyLiMSOnqizJWt66jR1hZq53QDliW4vT7yGaw2WC2cfuj/RibA7tIkxHepFA4Vg3rF0YPuhcUM3YgMFamLn4OFJGzNZ020R6/OtTcVAT7B7Azvh7gW9l0GZnIiY+dZb50yeOnNyvg156wjAZkv62mOn3Z/f0JPnpZv6udpX7HXPmHnatt8M5hQsi4hM74F3kH+DSgghhBBCCCHkUMAXVEIIIYQQQgghhwK+oBJCCCGEEEIIORTwBZUQQgghhBBCyKGAL6iEEEIIIYQQQg4FIaK27h5SP3s2nvnr/9XLZWNJ9AxYeMyzJC5gCTUGrgViPGsXWrpSsM15NjU0+wbHAFehldbRxKFRsNscm5hz3S1VfrC1ZWI+u/GgKm985IyJ6VzV99p+k/1ZRgKmMTSLusZZoLLyXWsW9X6MAscqLy/AHohm0Vv3ghjHAJfUtSWu3tRWvXrNWvaKUldosGtteNk1bWttWhmkMQxi24iIFGC99IyMaAV1zatOGyKhAMMyNJc3hqs6WF8di2Noabthmlv7bjEBi++WTZ7Guq5ffceESDrR9Zm17FgbnIZ7H9N97E2b2Q6Ym7fsdbOBLgf7mAvl9jz7rnfMyx0QD0vVgL6qOUbEBtg9a/Yh8rruz9zpzwaMm9yZOyuYByeO6Xcy1Q9WOHbDcqrP88yTAgZQNJSi/VbEznmJYxxHvH7AY8WS0xY93abdlYGJOdndU+UjDRszLXWePnVFJ3vjM9bEvXpeP1ht12omi7Z+iJuP2fG5/2Z9naxtGwz7SvqOZRtszsnUM3pD/Ryjd1jW92919IS71IAJWERqme6HWWlzcog56cTkma7PWsv21ZGGNgavj2zfvPisVpsun3fuNQCTbj1A2Zwi46O6PD1l+7zW0cc86zHuf0rHHlvCeKxKG4P7puCZa2HPhvsqEZFirPsm2dblmjNv59A13t5m1tXl6Yozubd17rjWV5ybRrY/E8ewj1RgRq6t6snqWM+apU+29LFGascn2un7MzvOp7BooUFYRGQPjMHXrqyamPazum/a1+x1QGhs1jQRaw9HMzjuoUTs3tQzN5ewXpZHbXutHdNzcq9hfyXG9T2dPNPzyyZm5Wv6Xo1tZz4roH2gWNXsQwyPgdH4IRMixVmdO8srQxPTrtn54U9+5O9+Psb4uL3iYvBvUAkhhBBCCCGEHAr4gkoIIYQQQggh5FDAF1RCCCGEEEIIIYcCvqASQgghhBBCCDkUfHslSfCheXSkRCi6kQVESt6H5/ZV/C4lSXCdFEQH+JG+iJUhJM7H/ShFmjkykBLEAV79UJz06OpNE/Ou5RdV+Zn+aRPzJ7/5mCqf/vieidl905Iqb75VP0NmvwWXHL7L9z5or0DY4H2cbgVbznUwvzCX3OvMz8FYQ0mM7c8aiJQa9fn2lMHA2o2qm7oxGjcdUQu0syeViiAM8aQYsy4IJnq2zvUlLQ3pNHXZ+1C+nmoxRCuz121l+rzl3CbPGpgqOpmVhE3AjPLFnbMm5ktf1ZKwtc9Z40X3EjwHpEXZsI1cwDGc3zxcSRJQ5XYATNswXyyZEJmCZ2HWcQRWIDwSzHVn7IkjvJhHcMZekt9ZPiYi0mroflhEpFS6lYZ7O/N/ivM/PGfhiFv6Yz2QRvvOwALJT3DELbENQpOulfMstXS+19L5yeNJ+1bqemwdb+pJeVDYSfn85nFV3rnSNTHdr+nn7H29MDGhBKmIM47KGgp8bHtNuiDPsn4Vma7AfNay7YV7jgCCmmRk65cN9b1RFnjrwrroyQBRhIUSLBErwjq9vGtiGjC/Prdx1MRMntP9lfX1M4xP2bZpn+ir8mxmF9nphrbLZPtOew1gX2D9Kmae8eZOlF55a1jR0W0au7ZNUeyUwTiaTKzFrNyEdXjdtkV9W5dRcOjh7ltQbue0RQGurFnHxpRtaAvcpzjzEMrG0oHtz1ynhZv/uIcbH7PzUHZaJ8LZ1W0Tg3Kll7ZWTMz4mm4MHJ8idvxVIMP09oZJQ+dOp2v3G8stZ5MLoCAN801EZA1kaJ3cJs/ORI+1a307B+/savtTta0fvL5hE66xqcsoVBMRiQnMyQvITUVEvvyPnqQkiRBCCCGEEELIdz5zX1BDCGdDCB8LIZwPITwTQvi5g+OrIYQ/CCE8f/D/9kcbhBBCCCGEEELIgizyN6iFiPx8jPEtIvJuEfkvQwjfJSK/KCIfjTE+IiIfPSgTQgghhBBCCCF3xdwX1BjjtRjjFw7+e19EzovIaRH5cRH50EHYh0TkffeqkoQQQgghhBBCXv98S9+ghhDOicj3iMhnReR4jPGayK2XWBE59k3OeSKE8LkQwueqwcALIYQQQgghhBBCFrf4hhA6IvIJEfnbMcbfDSHsxBh7t/35dozxjt+hvioWX9fEuoDFN8w39ArIv9DieOs8LGsbWOLcG61dtdxavPJEH6uwMiIyHGt11tSxzSGNplWsvWFVa7vee+TZudf5nz72I+bYI/9Sm8aym2CDfPOaOWf9Md3paFoUEcn3A5SdCi1gSUSrmGfMw2No0LsVA2ZTuG7lWHwjWOKSuu1zY/qtLWD6HVp1YbGp7b9Z3/7cCZ+z6Nr6ZEs6V9LMPlcJlstyoi8cp7YBw2y+VdXYuXNn7MG4SZy5IIUYtMCKiPQa2rw3czp9q69teJOR7vTKsWwbMaw3taKu0mmbZKzb2DdjQtmRCSbz08kYn3FO9gySeMyz9+Gxqu4YhPE8z5pegFVy5LSXFXUa8DndOtfAMAt1xjEt4hu8DQOdK7Utmzs5SNK9vjNzlWcxBRF44Zib4xE9b7/twauq/JbudXPOjYk2RvYLe3O0KX9t3f7MunhWX6d11YQYi2TqtEVSwHruCI1tn9vcKTFPIcbLEzRYlk0bU8Axrx/QoB1zJ78y3NvMN8Y3HWP8fh/WiGf1/NZ90dm3jMFi3XRsyj19DPNPxNpsPQk4thcayEWcNcuxMifwGxO8vWAKbVivga21YQ2q3bo2uGbB9sMU1pH9iW2MnYF+0LFj/U52dNLVdmy75yBzzhxjMFbRmJK9v57CGGfpxjExszJZmaxBbi/bnMzBkps4uY0WX4/xUA/I5IZt0/qmfhB89smqMz7XsM6O6XdL37ux4RiEob1GZ2xbnDq7pcoPLO2YmP5M3+vm0KqbR9M7vxN4v0GkqnRjjJ33iumevne2YzcG+Z599md/9TWw+IYQchH5HRH5FzHG3z04fCOEcPLgz0+KyPrdVoIQQgghhBBCCFnE4htE5NdF5HyM8R/e9kf/t4j81MF//5SI/P6rXz1CCCGEEEIIIfcLC/wqefkBEfkrIvLlEMJTB8d+SUT+RxH5rRDCz4jIN0TkP7o3VSSEEEIIIYQQcj8w9wU1xvgpMf8y/WX+wqtbHUIIIYQQQggh9ysLS5JeDepnz8YzP3cHSZLzQbuRJC0Q40qSQMDhSZLwWPBey+d8tJ0418WPv2uZ/bi/kesPp+u5NX+U8DEzfnAvIjLCj+6dr9xTEPb0lq1d+buPaHvFo+0bJubjNx9V5Usfe1CVT39CiwVERPKNvipXbftB+/CUFjjsnbNSkTH4l1BwIiISSv3sdyONEbFSGMy36Ah9rHxmgbx1hCsJSH88cVGBYqIta4xqrOsK1fZMiG0fT5CAghwUtziyqhK6eBGpjjcXmKHn5HbZ0u3TecA+6NuOXVPlwpEkXdhdVeWtLS0kiH37EMlkvmECRXDRyVuUpQSnz13JG16nhCSc2uQOIGRKprrOOIZErHwjOJIiFKNgWcTmm+MdMXi5U+jpwkhsRBzBkNPuZQOOgYQlbdoHbbX0g620hybmaFPPr63UirtQirc3tYKV7Yme7zf2rCRjvAeDbeb0OTxHrPS9s2t2Tm7B9I9jOpbHGwAAIABJREFUWkRk8CAIyk70TUwKedvft89pRGuu0HCRfQvkrpPLOCZQUJYO7TkoJEsWyG0P3P9gHouIzECuVDpyu3xFV+D0qhWsvLmnO7AJ5qmnt0+Zc154Vh9bedrOk8sXdC4npe2X/mk9IHcesW06PQOCvoaz/wEhX9izk0E6xIXXhBhhm5FVOeuw2VMusDf0pD8pzO25s8/L4DxvT1nBmJ050r4CjqHgsJo5oj9cDwY2pgYyHBS8idj89wSVOIfM7HQms2Xom54dWI22zh1vZphs6Lmz+ZIefDi/iYjU+rAOOyK2oqHbYnDSJgbOi6FrnwHnvGzD/t0hip4Su4yY9bGEOWXWdta9NuRp2z5oWgNBpZPbpbPWvPiX/+a9lyQRQgghhBBCCCH3Gr6gEkIIIYQQQgg5FPAFlRBCCCGEEELIoYAvqIQQQgghhBBCDgWvqSSpcUZLkubJZ0QcWYrzsXXMUCoyX5LkfeSLH8IHR8TgipPmgNfxPnqv1/SH062a/ZC6U9MyhKKyP1+4trOsyuNNK1IS+MAepSwiVlKw1B2ZmFNd/XV8K9N1fmmvZ85Zv6DlM93nHfHCRf1Bdt63IoEq088+69jrTDv6Ocu67Twj/VlA4IMiIE8YYsRAdUd0gDKGhiNJgn5otuyX8b2W7hvMExGRaanbZ2tojRx7u/pY3HdEQCOQiIxRqmNOsTg/FkOpVOWIp8olffHmUSuk6bV1W+wOrYRldHlJlVuXbYVQIoUCH6/PCxhqnvSkbIKgw8sLnL8WmBeD96NGlL45sqV0joQrceZAnL+8mAIETSjsELFShXJsxRDprj5W33akNeB4wzEtIjLr6vJ0xUnUJT3WsrouJ464C589SZw2NvI9b12ZL+hL4drYxiIiu9ttVa5dsona1v47qe2DuMvJpcEpXaH+G+2c3FzVY89rC8yDKUr9RCTd0h2YzGxjFCD2iHV7L3wOV56I+wC8lSNWisX8n+ubseblDsR4AjyUofnjEa7j7C9SeM4mSBlXG3Yu7da05HBnavcSz107psrZs3bSW35B16ex7Yw9aOaibtu4BCENru8iIhPYcuC4F3EEaeg/Kux1E6yyk0qLrGEoxfNicL53hYFG2jRf+olt7M1Di0icarBm5KZxLCiBExFJoBFRDiUiUk9h/5PbfXGvpuedEw1rbTqSW2Hb7bw0XjHHvnDzrCqvX1w1MZ0Lej5rX7Ntmo2xz+39y3z+XnW6BGW7vZaidWeBpze323esBYSe3juDc+zCT1KSRAghhBBCCCHkdQBfUAkhhBBCCCGEHAr4gkoIIYQQQggh5FDAF1RCCCGEEEIIIYcCvqASQgghhBBCCDkUvPYW35998uWyMUx5tkp4hXZj0DDl2PCMxdcz5oG9Cs13t4KgCCa0GB39IlZ5AUNjvWbtad2Gtuqt1K1ZF82+z988amOugGnPqXNM7mwDExGJYKFN6mB3a9hnaDh2YmQ40prcYtuaWOvr2p5W37bXSXVzuXbPQksvZdaxMca8WlugbdCyuoBB0stJPJY5MXlW3rEsIlJLwdbq2kahPo6qMIfr1MDgh/fxYjwyiPHu3S+08fOrN46bmOo5rbrrvGTvlffRkmhjypoeE8Y47tnE4VjlxBgjtGP08649D8/Oh9f28hRJpvq5E2e44nVny07/du5sxBVxbMCOfRTn18nYNljc1HlRv2kbI9+H65oIkRLs3GhhRnu3iNPunkAemsc1XaNY1OnPogPz0Jo1ep84tqvKyw27RmwM9aS3cVWb35sv2TZurtv6IOM1KB9z1s8VXee8afNiOtD3r122pt/mdV1OnTwtYNmYLTkxLYyBOXDVGtHRpF5Vjk12pJ+hHM43oicTx3QKRlm0iXu4cwGaYXEecgzy0gKDase2xXJLL7K1zPbnaKpvtr1vTb/Fhu6s2radBGs7ULayVslHsO46j4XzMs71pd1umFzyYvBY2fD2qnjAXgfnB9cqDM3szSko6bWWYXsO1rny7Niw73N/a4ZjCDbguJk5iYsWbc8w6+2t5oAW2nrbzqWY20v1sYlBq3Z/aueq7b7O9/GuTR40l9d2bJ+jiBhzQMT5rRSL/JYKqDLmuohIhXmBv4FCxP0tFJd++gO0+BJCCCGEEEII+c6HL6iEEEIIIYQQQg4FfEElhBBCCCGEEHIo4AsqIYQQQgghhJBDgaONeQ3Bb5s9wZCxR9gI8z22ex0IcWIW+a4bpR34kbR3EbxX5dy7BLnRZGq7ph/018ypYwA40dI2kHecvGpiPj89q69zsWliYgJ1dAQmMdV1jomu8yy1VpGpkTU4YhSULy3ZD9hReDGZWqlCsqfrk+/Yn8dkQ12u7ZoQifvQfzmW7TlGUOOIBIy0yfnYHx0Bnl8F1RVuHoOQwJUqgIzB+wjfpNwifgJ4Bk9Whe2TzOwYaWzo8up1Ryo11Q82XbJ9Pl7T1570bH1my7o+JUh/UAgmIpLWQETlCK2siMqCc8FsbBusApFMtuNIRXb1c2YDR7YBchnM26n259yqX3cB2dI2XGhk5RHo55u0HSFHT1dweWVgYo6c0InRyuaL2ArHJDMtdRvOwKaC/SIiMoNzSk84B8dwDRGxcrFOzQppejUt6Vir901MN9MxdWcQVyu6PrNTkDuO0qKCTJ2UNicnMLATZ31qgs1oBSdgEVmGY7vfa6U6f7z5BlV+9hsnTIwU0F/eeARZSg0Ec9O+XcMml7uqXN+0fb4Ej+XJeqwA0saURnRmY4x4zZGshYB2R5hvHUFTnOkLj/ZtBUcJWAY9YQ2s56kzd9ZP6AZLTnrCTH1sWNrxiHKquG8bLB3o81LwiLlOvwXkVCkM2eDIs2ZLIJvp2bkqa+ljnjyuLGDemTid3of9z75+bmfoSW0bpVyesQ+Ozd9u+/sE3Be40qb5MTi2vL0NigcLEF9O9uyF1zM979xwctvIL52YBCWWHbufDV2dPKXTXiPI7bBr65zvgXgN0qu0y7CULd2AnmQNx547zp08faXwb1AJIYQQQgghhBwK5r6ghhD+WQhhPYTwlduO/a0QwpUQwlMH//uxe1tNQgghhBBCCCGvdxb5G9TfEJEfdY7/WozxsYP/feTVrRYhhBBCCCGEkPuNuS+oMcZPisjWa1AXQgghhBBCCCH3Ma/kG9T3hxCePvgnwCvfLCiE8EQI4XMhhM+VAyu4IIQQQgghhBBCRERCRJWiFxTCORH5cIzxbQfl4yKyIbfcXP+9iJyMMf70vOs0zpyNZ9//5MvlCo2unt0zxRinvgvEoG0rpFbLhjGehTOB89DKmS5gsiodu5tniDT3BvNes24NcEfb2uz4YHvbxKBt8RPPvMnENC9oQ1hlZYYGYyX0bGBgU/P6Co9VjulXGmBDrVtbJRoaK8f6V450W4SBtdaF2Z0NuPhMHqljSUSzXVV3bH1NyMGWVdQlC9hjA+SlsU+LSAK5nHljBM6bznT7jQdOomzrY80bth+a67rc2HbuDWq7WduxqrZ1O5dOdUwOOl2DhzBmEXOzN2bQomf6V0Sqtu7PrG3HeZ7rGOw7EccM7vRN2NGVbtzU56A5WUQkG6G53MaUdTAld23MFH6sOes5+sUmHPMU1TCfRrS3ikhAHbZzGZx3AoyrWsv2Q6upjYzdxtjENBewCvenOjF2htasPtzXMbFvkzAZgaHUmXc8y+XteJbtqjZ/jcW1GmWft87TeRqc66DltdG01ssMrMe7mx0TU780f9Ga9XR96pD/rRv2nGQG61PmrOdgH3XtmRizwFzljTWcm1yLOw7Zu7Gxe5eFoebNeWb99uZbWL49Ay5ep3LmTpwvvH2euTeujZ7tGY556yeul6lz7wqeq3D2JPOuKyLSyHWDtXI7Rpq5nndqrp4Y6gfJNHUsvqMC7M5TOw/1xzrhR/t2ACSw9tQ2bVvU4bcqJPYxTc5NHPP8DMzzOK7Qditif6uBN65w7Hl7wUX2s/jbLGLN2f/g3Om8a1Tw2yyyLT2ZNzbsuEqhTb25agbTa9Ga/wwiIhd/9hc+H2N0vPCLcVd/gxpjvBFjLGOMlYj8ExF5191WgBBCCCGEEEIIEbnLF9QQwsnbin9RRL7yzWIJIYQQQgghhJBFcP4xjyaE8H+IyHtE5EgI4bKI/Hci8p4QwmNy6x+GXBSRv3YP60gIIYQQQggh5D5g7gtqjPEnncO/fg/qQgghhBBCCCHkPmbuC+qhY5GP+ReQBETHjBJAGIUftIuIhASEHOY68+UzifMPq9MKBDCF/doaP2AfTezH6RuhrcqZsSGIPNrVRpoffvtXTcxH+29X5c4lW+mipcumKTxJEkoVUF4iIqbGnhgLxEWVk8pVMV+MhXKlmFuRQAkfnqd7+l7pyD7DbBUEKw8OTcxSUwtVPBkCyhhQyiUikoP8oJFaYVQLRC2tzNoG6ok+r3JsFiV0cg751cQv7p3rjhwbyIX+mio/e/2YiZFLOuFqO/PlJEX7LiVccOlkDPIZp8+NJMy5rhF9tG2+oQyhcqQ/FcQ0nLxtZBNVzrp9E5Of1ufNQIqx68h6trd1P6Sbdh6qwy8myxyBe+clXY5XnDkv08c8CQVKfdwYFFh5kqsa5oW+ULFjT9oNui32nPksmepjqfUoSTbSZZTGiIi0FxFywDHHcWIlYQtI3hKsz9Tm5AIOFtOfrhiuBXIvZ87rregGO/XQVRPzXP24Knc+bXP55Kd1padL+s/Hq/Y5p0sgurGXNfOQKzfCvvLkVChPcfL2rkSSKFjx9kOQy4kj3EoHun1qe86txyCbcbxFJTxXYZ1Xdsx6+7wxzBelTe5kouucwVzuCXNwbp85IpkJ9FU2sO2F82Lq5MW0p8vjrn3QPexjr8+hj00Xe9YfFG6h+EzsPiqt2YGPEr/eml0Aasdg7Dl9tb+vB1d10zZ88zqIzdZNiCRXdLnK9YPO9LZZRKwsyBufKHaNTm7jPio4fY7tHh3xGvZFsz0xMe067L9O6uLWnvOgV3QbN6/be7euQ/2cxC2a3oO9Ml7Jr5khhBBCCCGEEEJeNfiCSgghhBBCCCHkUMAXVEIIIYQQQgghhwK+oBJCCCGEEEIIORTwBZUQQgghhBBCyKHgO8/i67GASdee4xyC63iSs3k2YDTtilhTm2diRatecG6OZl/vXsOxNqRelyUTU0u1Pe0dy5dNzF/4/i+r8if67zAxaPZCm6FnykQLIdrxbp0H5tqRDQoF2JQXMKbGZUeNCUa6ampNcmgzbD6kVYVnejvmnO2Rtnuury+bmOklrSpE26GHb+6Mc2OMnc8xbqJt1DOJorWubID5tOFo7FpgNG5Z02+zoY91O1Z12n9I32t8xRrpmlfBMrxn23QMguCqZ7WNjY6uT1nqHJzdRE2nSH0TDZKeZRiM3gPbWdlAD5KyZdt0muuYYmLbYh/7zxlrJVpUl3RbYDuIiKwe3Vfl2gmbTN4ch8zA2jhxzOVjmM+mQ2dSGYJVe2AfNBvON+nm+5D/8Fg4D4jYNcIzlJq1ZgGjqzcvot7ctY3itT3D7BxrL84nIiLlkr556NibZ03InbozrjKdlLXMTjJ1sJB3a9ZWuVbTVtCHHXXnT5/+lCp//buPm5h//LkfUuUjn7hzDojYdQ7nRBGRAoajN4ZNOzvtjjGeVdWMa2fvkMA6F+CcNLMPmsCeJHUs8/CLD2Q2s40x7OsxnOzamNqubnfP+o3m9LJpE3na03WMXZtfYVXPaSD7l2rfzjG1dV3ntt0yScNsA2w/DI/pZxiccfYtNVgjJnYyyPdRGW/rY9Jgkd+AgTZZb3+GFnwn/0ew9xrWbe6EOu4L7HyxtKRt3Y2VfRMze1hXcq9vJ73Zuj7WuqYftLFpTpEcpPeerRuP4W+2EHF+a0DN6Qi0MDtjOM70cw53bYVGmVYPL3V1+z184qY5p3Faj5HNsX2IjT29V8X9kIidY14N+DeohBBCCCGEEEIOBXxBJYQQQgghhBByKOALKiGEEEIIIYSQQwFfUAkhhBBCCCGEHApeH5Kk1xCUIkX48hylAd/sGJKBVCR1bBsoTpo6QoKygg+pQTIiInJZtLAnc+71fb2Lqjz48+dNzFMfeYu+DnxUvsiH+3XrF5L6tm7jWdtaRVAqsveIvVlY1XKNMN9BJCvH7Uf4bz2ibVCX9ldU+eufe8Ccs/oVXT67Y9u4AMHQpGsrONPfpn+Tj/DRbGBjUKjiyT9QiFA6UgBIL3OdxBFaCUgeZns2b6epvll0pApZWz9E7SwmnMigq+VFzRet8KL7AojNLtVNzOi4PjY9re+dH7GWnUmq7915wUo8Vs7rtigats/7kE5OV0ljXV8737MxJXicJisLDMiR7puxIwyZlmiAWWBgOXIXlG14fR5yEIbkNiZd030RjjiSnwrmFEfygFKrWC3wXGgV8dpirO/lSt9mkJOOSKOC9vGkTfmevnZupzPJtDdDpuBvK3o24zKUIjkSjwkIrMa7ViQmswXaCw95YiDIg4823mRCULDyUG/LxPzwd+l17eONh1W580dWPlbf1eWdtzh9tazbK92w46i+o/vKFwbqsuNFtHgSLnTqwNRUePkGsh6U3YmIZCC2aTWtVO1YTydh45yV4aC4MXEWsQRyDssLx8C1ayDlaqa2fiU0/OVBz8S8cPWovu5zdgHtXtT37l6y9SuaurPGK7ZDp11d9sSICfqhcOp0BUjzr4tJWHl9JTiXOvMt7AumA7svmGR6DsH1QEQkqYFsqW7FWK0H9UZzekrfa3PL9lVtXQ++1Ka2rYsjrktA1Jjifk1sO4fKNjxe2+ub2bJun92JDtrbthvIkOA7jDOBgIgw27E3j45U85XCv0ElhBBCCCGEEHIo4AsqIYQQQgghhJBDAV9QCSGEEEIIIYQcCviCSgghhBBCCCHkUMAXVEIIIYQQQgghhwJafG8HZGSLGHmN1dcxYFUV2vrm2zSz1LFVwnmelXYy1V2KVl8Rkf5YG0ovhRUTkyXajPbO5Usm5soPaf3j5v97SpVnVoAoFUiFp989MDGDb2jTWO85e52NPwf2wNWRiSmf0ao7vLeIyMrbb6pyt2HtrH/yp29W5SNf0A1/yjH0jnu63Tffaq1noxNg5XSeodnW6rhmPl+VNnHszuORfvhy2zZGbVPXseaYYWtgsKxyMPpZWaWU2OeOTTYDE+uptV0T82BnW5WP1K3FdxkUpcn32Xv1S53/E1QXikgF2kscD3WjSBTJQWl88bE1E/NHz2pL6MpnrEH49CfBPl3YZxgd0426/agd56PTuj4xc9p9D2zAu2AcXMBc6DSfGWtl3ZmswBzqWXNjoZ+rnNmxhjZgceZXNKB7ptMAttgU7JCJY5PF63rzf9rTda45YxjP88awuW7imMEL3Z+DTWunbFzXMa1r+s+XLtl7x6CPeeMc88DLi4gxi1zHiYl5hBh7sx0wUH9pz7bF8TU9yb3z7GVVfvq9J805tU/odeXch+cPkv4p59hZXfbWS5h2XPu6kah6pt85Md51kxnsW/bthQswr+4Fa27eg9x2xLqmfp6hGusYCicGh4TTFhXmTh3Kjq047+g+7nbsPuHBE5v6wAl77xuPLany6KUlE9O5qCvdvmEbrHNFPyiaf0VExqv42xD0n3vW6LvB7U/sK8+Ibs6bb/SOqd1HVbk+NqrZCWNc1wtSa0n338Nv0L+pQUSk86b54xp/+4Znn0Yq5znHhR5H1wc2LzZf0ubo1kU75zVvwHpZ0+XJEVufyfH5e+nmmt5r7XStDTh93vm1D6+QuSkaQvhnIYT1EMJXbju2GkL4gxDC8wf/b99wCCGEEEIIIYSQb4FFfobyGyLyo3DsF0XkozHGR0TkowdlQgghhBBCCCHkrpn7ghpj/KSI4G+4/nER+dDBf39IRN73KteLEEIIIYQQQsh9xt3+K/TjMcZrIiIH/3/smwWGEJ4IIXwuhPC5cmC/NySEEEIIIYQQQkREQvRMQBgUwjkR+XCM8W0H5Z0YY++2P9+OMc79DrVx5mw8+/4nXy6j3AIFCiIiEcQUWBYREfhuOmZWHoECDFekgRIi514JyIsSOAfLIiIpnJOl9iP8GhzzJUn6mCdAGs30B+Ljif1gvAJpgVef5Zb+iPzNKzdMzKNtfewff/Y9qnz0j2yHbv2Ivu5f/e5PmZjNmf5A/MO/924Tc+5faWFOdIxRG4/rj8q33mPFBig5aX/efui9+rWZOXY7nqCgbMDH6q48SNe5tH4JKeBb9FnHxsy6Oi/ispWw1EC2VHdELUWp6zzatxVK1/WD5CDO8ERU0/+vvTOPseu67/v39/Zl3iwcLuJqUjLlyDJQ21WtBm5St6m3tIhSIy0sFLGaplADREECI4Dt9I8GAQo4adOiLYIWLiIgAWI56SJECFJHblM4gG25smRHiyWZpEiK5JAccvY3M28//WOegvktw3nhjDjvhd8PQJDn8Nx7zz3L75w78+7nHdTXKuzz/WDrY2VfANBaMUKhRvDzNSvMKfh5lIyAI7PmxQtWwOHiUMmfNz+p5UbHphdcmfvGtUhjPRgYL84cU+nu98ddmak3dH3Ks7ceowDQLfv7bI4bicKUEWsE4804pmIxVtmKSHx7wbZp5Miw4TQUrJgDg0tZ6Yq0A5lRwwiidHciEzWxOY1tGwDoVMy6EghgiuYzSkXvCEPbhKbV44H06oT+IfDxfX4MThb1/Gt19bi4Hgg6bi6YgXDD32h+yczHSDJohmAoScrbfYE/kdsH5IOLmbmfLfp1bmxMC0FOTur2et/4jDummXSl//j8A66MfFsLBPe/5AdPYVHH5G7Jx7x2TTdYtNZ0SmLKBDIjE8pdnB5ErLT9VjHGnjvy5VhpUzR27HFB+HeCrWA9cpIkMwajPeZA483EuHzZ93mlrINKJD3sGqHQwooX0vSu6WBQvuIbo2Q/92jar+WXFZitVyhSsnO4W/Jt0a2aIFz295kr6bxI+tlp6ov1Gn6OZMw+wImy4O+jV9axQDr+Rgs3jEAw+P2a3Z/ZWA/455pofFkxV3Xfmitz37TeO+wveUmklTvavX4mNFppIgHkeF6vGQfzgUEz4HMPfvWFlNJDAxUOuN3foF4XkcMA0P979nYrQAghhBBCCCGEALf/gPoMgMf6/34MwB/uTnUIIYQQQgghhNytDPI1M08B+BaA94jIZRH5WQBfBPBRETkD4KP9NCGEEEIIIYQQctts+23gKaVHt/ivH9vluhBCCCGEEEIIuYsZSJK0W9wxSVJQRqw4KfjdsRgJ0W5JkjIZK0nyb3HbPCtN2ihjRErB2+Dtnm6M1Za3BKwbcVIK7CT5nL7WwdqKK/O3Dryp0scK+q38f/+U//YhKze68TcCt9anbqrkz5z6livyuxe0OKn51QOuzOH/o88jS8Fb7lUtG+hMVV2R5fv0m/D1Y9uLZJx4IRJ9mLyed9h4kUxwHjdOo/Fvxle+6F+ELxZ0npV7AcBaQ4+n9oxum+qlQDYwwPv0Vm4RxQIr+uj4rnLCBis7AoCc9qJAfFOgZxwwVn7QsRIIADCig0wgZSmYdi8VvUijWtDylELWV3C9rfvhxqIfhOmqHtuVGd8W5Zv6vnINne7m/THdopF7BSISKwuy8o2NPHOtciC0MuIRyQfyONPOxaBNC3ldJoqvlroZ62uLXqCWv6bLVK7685TmTRuvB3G7qudN/bhv97UTehzkjJQL8KKRbtvPx56RoVlRXLbg2yZr1s9uIBXpLOlOzy0EQhMzlLvFSEhjDwqkIrbMACKlaD2HuS9779FYmqrqALK/7NcVKyM5v7jPlVm8qEVKlRm/AGSNK8W1DXystLEU8LIgGycj0Y1bawaQB0n0mbwBxCxWXBdJa6RlJGYNP0cyTSOFCWK7lQU5gU60D7V71bAfBljPB5D4ZcwakQvWEbs/s4JDAGjM63hVmtE3Vr7hq2dFcNHeZt1stTqTQSy19xWNAVvlQcxYzWCAtXRepuXHRWFR59l7L9T9tdtVfUz9qL90e9I8DwTXtuPUSi0BIH9YT/SDE36/bbm57DtnfUFvkrIrus/t/ADg5kMvEEB2x0xcnPBrz3TNx8HnPv7reyJJIoQQQgghhBBCdhU+oBJCCCGEEEIIGQr4gEoIIYQQQgghZCjgAyohhBBCCCGEkKGAD6iEEEIIIYQQQoYCWnw3lxnA4ivGKGhNp/b/ASBrzpuNyljLaiYwgBpzW1TGGhlbXd+oKw1tW2wYq29EpdxyeaentCX3M/d8Q6X/fqXhjnmuqQ11Z1r3uDJdo9Vb6ZVcmabR5HatigzA+TWtm/vGlVOuTOPsuEqPXfTnKS4Z46GxmDYCEXHTSBt7gaHUGuoyvomdhTYQNzsLYWT6tVbhXj4wMpb0eMqWvQJxfEIbLA/XtKK3kvPWy+WWHm+zdW+fW7qhNa+FGT8my7M6XVwM5lF7ewttc8L0335XBK0p3dCpYK4V9IPrq8AgbEWFtu+AwBhZDC5mzI7WZAsAufz2ptrWumnnBa37LCz4CjrTY2BTbk8bE2XND+68qV+0EqVeYB00uJic8e1lLb7lvB+nRWNLrhW0qXAs782F1qS+2PKm34uLOkAsvzXhyoyf1e1cvRbYPZ0B17dNu7K9YbxljMrOpjzux02mpturXAmClWF1ycft/IyOBcV5V8TFquge7LyJ4uJ2xwDR/sIUCPYAbt8SzM+MiaX5UqST1bRWvQ47M68bI7seWLXLxh477ftmYkobNseK2/dfo6P3Dq2230s0Wzqv0/bqWmvoDZZqZEw7Z4PYFc1rS9tcv7MW7G3WdJmMjdPBtxrYeB+NN5cXnCdl/vI24BSs1TDjK1fx8Wyiptdq+20NCws+cGev6DlbueYvXVje/pnBxqaOD4toTpr0EX8Ph44sqvTpKa8e3l+oq3QxUDcvtvW3Dbw0d0Slr531m4DaeWMH9tXD2mGdTvfgzHnZAAAao0lEQVR6k+2hKW3kbXV8p88tmCA362Nncd5YrNddEbe/sLG060/rvh2hZ/c6CPaU2y/LAIALv/DLtPgSQgghhBBCCBl9+IBKCCGEEEIIIWQo4AMqIYQQQgghhJChgA+ohBBCCCGEEEKGAkqSNpd5hyRJGfPmci4bSBXMtaMyeVOmkPMvg1txkr02AKy1tYxhYdW/wd5Y98IGd62ivv6JfQsqPVnwkqRXrmsp0vpV/6J+YVF3qASuF/vyd6fq2yuNG1FL2b/lbuUMnWYgeZjTYo/yVT14yv69feTXzbgNXirv6tOiUw4EGK7M9ufpBV3n5lEw17pV3dCFCS+FGR/Tb+ZnjNpmreUv3u3q9rJzBgAqBS3tiERiVtqxUvdv/HfnjYTlhu/PknZ7Ib8azHMrG8jpvrFtHuWFZayQIHB42LgYxSor13DHAEDetHPBt7uLi2YI9tb9QCnc1HmFJX9psUKfSM5gpr6VvQBeKiJWuAIvrwhlOANIdeyc6Nn2KkUyHH2jxZKPMZWSHtuFrA9o6209EJaX/UTvrukKSiCQcu0VrZe2Ley6F8w9e55IwFU192nnNADUm3pSLF7ywqixc3pwFxddEddX3SDmOUFIVMbk9dwYiM5rY+kAUpFB5mfQV6mtOys37+dj6aYeB/m6KxLEM5MO45A9iS9jzWaR9Me1cXCtrhGzpKCMEwxFoj8zZyWqs5UiBTI7h92rRvPKemRagWCuYcU7A4iUBtkXhyI90xYm1tfGvWXnwFgweAxzazpwL84Flry6qXTQV/lJvb+oVvx+wz6e1FeDNd+uUUHfTO/XoqL7jWzpaNkvYrlo42lY7eh49t35o67MldcOqfTEG35clOeNlDHaLxbklmnAzz8n0Ix+JWnyegPsN+Ix6fNe+9efpSSJEEIIIYQQQsjowwdUQgghhBBCCCFDAR9QCSGEEEIIIYQMBXxAJYQQQgghhBAyFPABlRBCCCGEEELIUBC4mMg7TS9QdNmfFKSgjD2uG+i2skZPmct4u9t40dt1LbNtPTTa896e1mprg9mZ6xWVLl3zWq+8lqmhFvyIxBo/OxVfxlpLA1kxkjHJdZr+Yta2GBksswe0XU5OaDvlWmDTXF/WN5G96bWExQV9XDbqFlOdrBfdIWOMqcnL+dAt6Wu1a4Gp0wyV9oLX0C5c1p1RWNLntUZVAOiZ/mtM+jG5Oq4PzJe8odraf0uBMTV3XDeinPD1aXf0uFxZ9arO3qLOs2bp3Ko/b6Z16zTg7baR+S4ZY3Bo2DR5meg8ojNT1hfqVXSbZo1Z8cAxbeYGgOn79c3nAiWujVWdQB9oy0Rx0ca45aaPQ/NzY7o+V/24Lc6b80aGZWMRFmNQxXoQPzI6xjTF169hrZwTftzWJvWkHQ8Mm6s5XenWih+3mVXdx9m1YI2w49KaWINdQaeqCzVqfn6mCd1/5by/z9NTWqFdOTDjylw8PaXS588fcmVKV3Ulo3nULW1v2+3Z4waxtVpTciSBtbb/ojeCjhkj+r7Kmiuzv6TzKjkfVBpdHQwu170Z+dqczkuz1nbux4m1c7txE+DMv/D25F4w97IDWNLbk7pNMxO+QnbelIJvOrA2+NU1fbFOyw8mEbtP8PXLmPXJpgGgZ4z20b4qu2JNv/5a2WVdprjgK9Tcp9Ote/XauBZ8U8O5OT/XHHaP1A1sxWavVZrxkyT3fTOHu4EN2Jy6Eoyvdk2nW5N+zs41JlX6W5d1WoK9oY2ddu8KeMN4acFf+3BL57WCvdfKUTMuxv21OhVrsQ7Glx26dmMcxaoBJNa7csxtsKMHVBG5AGAFQBdAZyc6YUIIIYQQQgghdze78RvUv5NSurl9MUIIIYQQQgghZGv4DiohhBBCCCGEkKFgpw+oCcCzIvKCiDy+GxUihBBCCCGEEHJ3stOP+H44pTQjIgcBfE1EXk8p/dnmAv0H18cBIDc5FZ2DIJIi+Repk303PZAkdTM6r5f8i9SFjBY2HKzUXZlKXgsILmDalcm+ql9qnzyjr9XLB1KKnE0HL88bV0R51hVBt6CPa0348zQO6uv3MsGb3UaCUah58UK5qPPWm1oukAJJUmFMH5Mq3nSwdtA0RiMwfbTNuQPRhyW77M9Tvq7PM/aWr3M2kENYrJynZV7mb/hhgtaUEYQEl5FFfeJe0wscUteIdwKhlfU+9Ip+/KOs65MpeIFJbr8WSnSnjWTKCnQApLaREnV9GzuZl+1fALlVK2cIxteSzfFlrFwskkeknG4fW+e5JS+uuDGnrRSpE7SFOU8okjHXjsRYlbKWNk2UvUls+riWNt3YN+bKLF7QUoza+UCKcePWcqpuJKuy4y0o0zHV6a76ZXd1QUtsCvO+wSpGyFELhGl2GFhBDQB0rITO1K9TDuJ2wcyjIAxZadOV1X2uzNWcuc+gz2umjyfvWXZlltr6PJXLPqjk1kx/BmuNlSvZ/rMyPsBLpHoVHz8yY+a+AovfynJZpZeXvA3wopUt5YL1vKCvVSn4Nez4AW0JE7M2Njp+4C7UdX0a9cgsZiR5QX/mizovE7RFzgiFioFgqNHUdWw1/DxaNm26EqwRsLLJph4EsuzPK019TArmeauqx0FuLJD4mb7KBPuN1DATNJAktd+tZVAnjwabJMMbFw6rdPFlL2gqm+5bOxy03zF97RMHbrgi9zygjUIdZ+8Brq7pdeTavDcDda/r/iwFMq+SecGwejlYbOTWUrXA4eekkVZwCAC5phF3tXx7SU/n5QOJZVo08yi4VqZlxm3J32e3aPa8ZspGoji3dQjmp2ufoIxEQrkdsqPfoKaUZvp/zwJ4GsCHgjJfSik9lFJ6KFsNLF2EEEIIIYQQQgh28IAqIlURqb39bwAfA/DKblWMEEIIIYQQQsjdxU4+4nsIwNOy8dmtHIAvp5S+uiu1IoQQQgghhBBy13HbD6gppTcB/LVdrAshhBBCCCGEkLsYfs0MIYQQQgghhJChYKcWX/IO0XNWX0CMCjMbqBSt2bcdKVMNFfGauCMVbU78wP2XXZnFU9ry9+zLD6p06aJX3TWOaj2ZNAMDqDFGZtb8PVijWe+4V6P1rOUvMI+NTenj1la8qTDzktZcFrWgDq1ATt2cNmbHwPTo7GmBDTi7ptsn2/Jl2jWjmzvmTafptLYF1ru+3Zur2sIpS77/Cgv6uKKWQ6J2wR0CedOYpQObZlfL+tApuyIur5cP5oixNmYCfWyvq8dFL+vHV9fa7szYkcjQ29y+r1xdin5Mto31uB2Mr0xDXytX9/2ZN3LuykxQn8umj0WnI7vhdubTKC8q0y3pe28X/XK0VNZjcrXq5+e+ca39vm9qzpUp7b+m0q+dOuTKLJ7VDW3nXnsimMOm+ypX/FiaOKfTRS+lhfT0uaM50q6a+lQji6OpXtDlORMqC2acpGDOpJy+r04xmDN2DnspLTpVHasatSC2Z3SZQ+Mrrsz+H9Lm5nPVg65M5XXdGNnAetyYNvHCWLXtWAeA3pRZL6NpflNfO78QTCQrRg52Y4NYhdfNcWuBTfOGOS4Zw77kA9u5JViHCws6L1cPJrppn2hctEwsyK0GFutrJh1IQ1ePmPOe8mvhPQe1An1fScePxaZffK5cM0G4Hph+7bcsLHmFdv6A3v986sE/d2X+du01lR7P+Hv4XuOESv+nlz7iytS+rqWkx6/oa6/e4w7BwoN6HFRO+LlXK+mJtNryMfmcyctm/Pgq5fQ8euDwdVemclyXyWV8DLaG4Hrbt/tSSxuLl9d1er3hj+mYbzUIbfXG5i/Bmp81a3V23ZfJGJlzcJv+mMD0K2YPKXVjBw6M0C7GBOZ3MfWxhmMAaE0NEEP+kvA3qIQQQgghhBBChgI+oBJCCCGEEEIIGQr4gEoIIYQQQgghZCjgAyohhBBCCCGEkKFAUgreNH+HKB07no4/8dm/SNsX/iNJQDIv/Ns0AFgPUFRGcuYF3uDRXMyL3BKdx0hYstneLf8fADLm7flMUMa+RD5ImWwg/cln9dvM0cvpmUCutB0TRf+i/l+ffEulP157WaXPtb244lf+9KdUeuyM7/T2D+sX86N7yH19QqUnz/k3xpffpc+9+F7/5rkVoUycCYwXprmW360zuke8fSNf1PXpBQIkS6Hg7+FgTRtMDlXqrkzBvC1f73hpwWLTSAECkUDXvmAfVNmO5VZXT77lur4OAHTmdV7xRiDbWNTpSGhiCcU75ra6gZCjXdX30Cv5+WBFXTDpXMXbBvZPaHHL4TFvw5ksaCFHMTAd5CMDwS6wHjTYXFOLNK7Wx1V6dr7mjklGABMJmpxIqRC0sY3tkZCpaOJr0c/hnJlr0TwqF3V/5YOYMjun7zV/TstSrCgCABpH9LVK+72szd55MxCxybKOVVGbZo1II2ovK56yaQBIhW3W1GDtcWthsPbYMplg/cyY9dKunwBQLOi+Giu2XJmDJg4WAqvIi5ePqXT5m1VXpnpdX3/5hG7U1fv8WLIytMol3xFF7eFBO4hDzWmdjvrKSpuiMWjHQTTX7LlTyewTyv4+i2XdD3YOAUC7oyf6yry/0eJbOiiPXXJFUJ7T9ZGev4f6Ebue+7EzflIvJBNlv2/JmvhaMHumXBB/r6zo/cbKa95cZyV0Db/9gZzW47aY9+2+sqDb0LYfAEye0e1Tu+DjjqVb2t6J2i3pwbQ+7dfqxpS+z7ZfItAzIa4bxX+zHKUg7rjYFDk/M9vEMwwQv6K9jt3bR/HMXDva/9v9az7nJ3Exp8dBNC7KRioVxbyeuZHZVS34nLs06Y4ZO6vHRWnO30Onos+78q4gxhzycfriZ77wQkrpIfcfA8LfoBJCCCGEEEIIGQr4gEoIIYQQQgghZCjgAyohhBBCCCGEkKGAD6iEEEIIIYQQQoYCSpI2l/krIElyL2Rn/YvUefNydbfnG+NGXb9cvRrID2Dax4o+Dn3LH1J76jmfaVj/yYdVev4zXgz0wMFZlf7ut9/typz8I/3SdvH8TX+xvK5zd8LfZ7eq3+ZvTun08kn/5n79uBHx7PcvkE/u01KdyZIXHSw2tKhl8cq4K1O6agUrrohzAliZEBCIDbzLBZ2qHl+9qhFMVALpTyCtsbRb+h46dS/0yazrcRoMf/SKJrMUWEVke7FBvqTrXKtq2UY574Uh9aZusOXlsivTW9YNn1nzcy/T0kKCwIWAZKQOvfz2Mpxeycs/xPRX3shSIqFDs2n66oa/z9Ksvi8r+AGAjplqnfL2Ig27ZgDwUp+gjL3PIwcXXZm/e88PVPqB8oxK58WP47ZZtGbaXkLx5voBlb66PuHKLLf02Gl2/WLY6ug8KygDgG5Xt3ur48t02jqv29Tp1PLHWDFQpuOtIk7gE3SVFfqkoK+SHcvlQCpS0xa1Y1O+P0+P31Dp8/VpV+bNb75LpY9808hKbnrJjpXNtCZ8X3Uq+kbbZd9ejWmdt3bEt0Vvn5F7lX3csXuQSDxl9wU5U8ZKWgCgkteTdrzgzXVW+jaZ317W89baPpf30pUjKp153Qutxs/r9smv+vZqjut2b+z3128cMCKxad2mE9N+AT05uaDSdr4CwMWXjqr0wed9/YoLZu9V9vG/WzDxP1g+M23df53gPPWjOm/lXn3t0j267wCgZqRSUYxZW9f33l4PNu4mpkgjkL41zH22gphihvJgkjBfxknCrNCt6a+dW9N5WR8KfN8M8jgVeThtXAxkUNZxGN1nu2aeR/brOTs54fu82db9t3rZW69q53QFx64F8lUfmvDN//7LlCQRQgghhBBCCBl9+IBKCCGEEEIIIWQo4AMqIYQQQgghhJChgA+ohBBCCCGEEEKGAj6gEkIIIYQQQggZCmjx3VxmxCy+mUAZljXHFQI7X9XY+SYL3rzX6elGff7su1yZe57VGrHJP5/TBS5qCyYAdN6vbbtv/aK3gf3WB7+sr5NbcWVWIoWZ4fn1e1X6ybM/7Mo0vquNgtlAQrh2Urfh5NFllS4FbbzW0sq1RstbaVur+h7yM/6eKqYJc0H9uiVzXi/6dcZUa4QDvFEzMmz2rF3OHROYYu08GmBepW4wQev64oVFr7rLm6ESWf9se7Vrvj7WViw9rd7LL/r6FbXoMbQv2vYL+8Ha+oIyvdwAZfLb92eyZl9jPS5WvX63WtJ5EuiUFxa1hTNzyZt+y9d0OmwvMyWaU75M87DWB1am/CSZrOq8yJi92tYXu3TNxIYZb+4sLBrTYzA/nYU5WnIDs6MrYs4TmZGlZ/o8409sDZFuLIXr8ABlnKHXl7H2cGvXBICutTlX/CTOGtv02Jhv+GPjSyr94MRVV2Yqr03qr6xoE+s33vB2+NpLJm775Qmrx3S6/W5fv0PTeh1JVs0NYG5Zz6P2gh+DWRMXBzGd2v4M91WmOkH1vDG7EBiETV9Vq16HOlnR7VPOeSVox3zbwFrbr5dNY7ru9gIzrMmy6/d0WY8JAHhg4rpKPzx21pU5mdf7n0ttH6y+vvxDKv3K4mFXZnZFW1RXlnzszCzogF+84dejgh2XpmvsngAAWkZCbq2wANx8lELwTRFFY8MOTP7Vog5gpaDPo7XFUjDfVFGItPeGllmIF9Z9Yyys6Lz2arDINm1gDC5mxmB0S25uRfPRrM2FseCbIUwcPDymY8yBkjdUV7Pa9FsMFuJ1E7h/sHLAlXnjoh/LFx/7PC2+hBBCCCGEEEJGnx09oIrIJ0TkDRE5KyKf361KEUIIIYQQQgi5+7jtB1QRyQL4LQCfBPBeAI+KyHt3q2KEEEIIIYQQQu4udvIb1A8BOJtSejOl1ALwFQCP7E61CCGEEEIIIYTcbdy2JElEfgrAJ1JK/7yf/mkAD6eUnjDlHgfweD/5PgCv3H51yR6zH8DNva4EuS3Yd6MN+290Yd+NNuy/0Yb9N7qw70ab96SUatsXiwk8ewMTOd3c025K6UsAvgQAIvKdnRidyN7C/htd2HejDftvdGHfjTbsv9GG/Te6sO9GGxH5zk6O38lHfC8DOL4pfQyA/14RQgghhBBCCCFkAHbygPo8gNMickpECgA+DeCZ3akWIYQQQgghhJC7jdv+iG9KqSMiTwD4EwBZAE+mlF7d5rAv3e71yFDA/htd2HejDftvdGHfjTbsv9GG/Te6sO9Gmx31321LkgghhBBCCCGEkN1kJx/xJYQQQgghhBBCdg0+oBJCCCGEEEIIGQruyAOqiHxCRN4QkbMi8vk7cU1y+4jIcRH5vyLymoi8KiK/2M//VRG5IiLf6//58b2uK4kRkQsi8nK/n77Tz9snIl8TkTP9v6f2up5EIyLv2TS/viciyyLyS5x7w4uIPCkisyLyyqa8cK7JBv+xvxa+JCIf3LuaE2DL/vs3IvJ6v4+eFpHJfv5JEVnfNA//y97VnGzRd1vGShH5Qn/uvSEiH9+bWpO32aL/fn9T310Qke/18zn3hohbPCfs2tr3jr+DKiJZAD8A8FFsfDXN8wAeTSl9/x29MLltROQwgMMppRdFpAbgBQA/CeAfA6inlP7tnlaQbIuIXADwUErp5qa83wAwn1L6Yv8HRVMppc/tVR3JrenHzisAHgbwM+DcG0pE5EcB1AH8bkrpff28cK71N8u/AODHsdGv/yGl9PBe1Z1s2X8fA/CnfRnkrwNAv/9OAvijt8uRvWWLvvtVBLFSRN4L4CkAHwJwBMD/BnB/Sql7RytN/oKo/8z//yaApZTSr3HuDRe3eE74p9ilte9O/Ab1QwDOppTeTCm1AHwFwCN34LrkNkkpXU0pvdj/9wqA1wAc3dtakV3gEQC/0//372AjmJDh5ccAnEspXdzripCtSSn9GYB5k73VXHsEG5uxlFJ6DsBkf6Ene0TUfymlZ1NKnX7yOWx8zzsZMraYe1vxCICvpJSaKaXzAM5iY39K9ohb9Z+ICDZ+KfLUHa0UGYhbPCfs2tp3Jx5QjwK4tCl9GXzYGRn6P7X6AIBv97Oe6P96/kl+RHSoSQCeFZEXROTxft6hlNJVYCO4ADi4Z7Ujg/Bp6MWZc2902GqucT0cPf4ZgP+1KX1KRL4rIl8XkR/Zq0qRWxLFSs690eJHAFxPKZ3ZlMe5N4SY54RdW/vuxAOqBHn8bpsRQETGAPwPAL+UUloG8J8B3Afg/QCuAvjNPaweuTUfTil9EMAnAfx8/6M0ZEQQkQKAnwDw3/pZnHt/NeB6OEKIyL8E0AHwe/2sqwBOpJQ+AOCzAL4sIuN7VT8SslWs5NwbLR6F/gEt594QEjwnbFk0yLvl/LsTD6iXARzflD4GYOYOXJfsABHJY2PQ/V5K6X8CQErpekqpm1LqAfiv4MdjhpaU0kz/71kAT2Ojr66//ZGK/t+ze1dDsg2fBPBiSuk6wLk3gmw117gejggi8hiAfwDgn6S+rKP/8dC5/r9fAHAOwP17V0tiuUWs5NwbEUQkB+BTAH7/7TzOveEjek7ALq59d+IB9XkAp0XkVP+3Ap8G8MwduC65Tfqf/f9tAK+llP7dpvzNnxf/hwBesceSvUdEqv2X1iEiVQAfw0ZfPQPgsX6xxwD84d7UkAyA+ukx597IsdVcewbAZ/pGw7+JDQHI1b2oINkaEfkEgM8B+ImU0tqm/AN9eRlE5F4ApwG8uTe1JBG3iJXPAPi0iBRF5BQ2+u7/3en6kYH4ewBeTyldfjuDc2+42Oo5Abu49uV2uc6OvgXvCQB/AiAL4MmU0qvv9HXJjvgwgJ8G8PLbim8AvwLgURF5PzZ+LX8BwL/Ym+qRbTgE4OmN+IEcgC+nlL4qIs8D+AMR+VkAbwH4R3tYR7IFIlLBhvV88/z6Dc694UREngLwEQD7ReQygH8F4IuI59ofY8NieBbAGjbszGQP2aL/vgCgCOBr/Tj6XErp5wD8KIBfE5EOgC6An0spDSrpIbvMFn33kShWppReFZE/APB9bHxs++dp8N1bov5LKf02vH8B4NwbNrZ6Tti1te8d/5oZQgghhBBCCCFkEO7ER3wJIYQQQgghhJBt4QMqIYQQQgghhJChgA+ohBBCCCGEEEKGAj6gEkIIIYQQQggZCviASgghhBBCCCFkKOADKiGEEEIIIYSQoYAPqIQQQgghhBBChoL/D4L5jk0afPQcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 40\n",
    "F = 16.\n",
    "\n",
    "def shift(x, s):\n",
    "    if s == 0: return x\n",
    "    else: return np.concatenate([x[s:], x[:s]])\n",
    "    \n",
    "def L96(x,t,params):\n",
    "    F = params\n",
    "    return (shift(x,1)-shift(x,-2))*shift(x,-1)-x+F\n",
    "\n",
    "x0 = np.exp(-(np.arange(n)-n/2)**2 / 16)\n",
    "\n",
    "m = 200\n",
    "dt = 0.05\n",
    "\n",
    "t = dt*np.arange(m)\n",
    "X = odeint(L96, x0, t, (F,))\n",
    "t = t\n",
    "\n",
    "figure(figsize = (16,4))\n",
    "pcolor(X.T)\n",
    "xlim(0,m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "861\n"
     ]
    }
   ],
   "source": [
    "p = 2\n",
    "\n",
    "measurement_noise = 0.01*noise_percent*np.std(X)*np.random.randn(*X.shape)\n",
    "X_obs = X + measurement_noise\n",
    "\n",
    "feature_maps, feature_descriptions = polynomial_feature_maps(n, p)\n",
    "d = len(feature_maps)\n",
    "print(d)\n",
    "    \n",
    "Theta = np.hstack([f(X_obs) for f in feature_maps])\n",
    "Theta_clean = np.hstack([f(X) for f in feature_maps])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "if retrain:\n",
    "\n",
    "    L1_err = {'ARD' : [], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "    L2_err = {'ARD' : [], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "    Added = {'ARD' : [], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "    Missed = {'ARD' : [], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "    Xi = {'ARD':[], 'ARDr' : [], 'ARDvi' : [], 'mSTSBL' : [], 'lSTSBL' : [], 'dSTSBL' : []}\n",
    "\n",
    "    for dim in range(40):\n",
    "\n",
    "        print('Dimension', dim+1, 'of', str(40)+'.')\n",
    "\n",
    "        # True predictor\n",
    "        e_dim = np.zeros(40); e_dim[dim]=1  # linear term in this dimension\n",
    "        e_dim_pm1 = np.zeros(40); e_dim_pm1[(dim+1)%n]=1; e_dim_pm1[(dim-1)%n]=1  # quadratic x_{dim +- 1}\n",
    "        e_dim_m2m1 = np.zeros(40); e_dim_m2m1[(dim-2)%n]=1; e_dim_m2m1[(dim-1)%n]=1  # quadratic x_{dim - 2}x_{dim - 1}\n",
    "\n",
    "        c_dim = np.where(np.array([np.all(feature_descriptions[j] == \\\n",
    "                                          e_dim.astype(int)) for j in range(d)]).astype(int) == 1)[0]\n",
    "        c_dim_pm1 = np.where(np.array([np.all(feature_descriptions[j] == \\\n",
    "                                          e_dim_pm1.astype(int)) for j in range(d)]).astype(int) == 1)[0]\n",
    "        c_dim_m2m1 = np.where(np.array([np.all(feature_descriptions[j] == \\\n",
    "                                          e_dim_m2m1.astype(int)) for j in range(d)]).astype(int) == 1)[0]\n",
    "\n",
    "        xi_true = np.zeros((len(feature_maps),1))\n",
    "        xi_true[0,0] = F\n",
    "        xi_true[c_dim,0] = -1\n",
    "        xi_true[c_dim_pm1,0] = 1 \n",
    "        xi_true[c_dim_m2m1,0] = -1\n",
    "        xi = xi_true\n",
    "\n",
    "        y = FiniteDiff(X_obs, dt, 1, 6)[:,dim].reshape(m,1)\n",
    "\n",
    "        # Test each method\n",
    "        gamma_ARD,sigma2_ARD,mu_xi_ARD,_ = SBL([Theta,y], \\\n",
    "                                           sigma2=1, \\\n",
    "                                           estimate_sigma=True, \\\n",
    "                                           maxit=500, \\\n",
    "                                           verbose=verbose, \\\n",
    "                                           tol=tol)\n",
    "        ARD_results = gamma_ARD,sigma2_ARD,mu_xi_ARD\n",
    "\n",
    "        lams = [np.linalg.norm(y)**2*l for l in [0,1e-3,1e-2,1e-1,1e0,1e1,1e2]]\n",
    "        gamma_ARDr, sigma2_ARDr, mu_xi_ARDr = Train_ARDr(Theta,y,xi_true,ARD_results=ARD_results,lams = lams)\n",
    "\n",
    "        alphas = [1,2,4,8,10,16,32,64,128,256,512]\n",
    "        gamma_ARDvi, sigma2_ARDvi, mu_xi_ARDvi = Train_ARDvi(Theta,y,xi_true,ARD_results=ARD_results, alphas=alphas)\n",
    "\n",
    "        taus = [0,1e-5,1e-4,1e-3,1e-2,1e-1,1]\n",
    "        gamma_STSBL, sigma2_STSBL, mu_xi_STSBL = Train_M_STSBL(Theta,y,xi_true,ARD_results=ARD_results, taus=taus)\n",
    "\n",
    "        taus = [np.inf,1e-8,1e-6,1e-4,1e-2,1e-1,1,2,5,10]\n",
    "        gamma_lSTSBL, sigma2_lSTSBL, mu_xi_lSTSBL = Train_L_STSBL(Theta,y,xi_true,ARD_results=ARD_results, taus=taus)\n",
    "\n",
    "        taus = [0,1e-2,1e-1,1,1e1,1e2,1e3,1e4]\n",
    "        gamma_dSTSBL, sigma2_dSTSBL, mu_xi_dSTSBL = Train_MAP_STSBL(Theta,y,xi_true,ARD_results=ARD_results, taus=taus)\n",
    "\n",
    "        # Record predictions\n",
    "        Xi['ARD'].append(list(mu_xi_ARD.flatten()))\n",
    "        Xi['ARDr'].append(list(mu_xi_ARDr.flatten()))\n",
    "        Xi['ARDvi'].append(list(mu_xi_ARDvi.flatten()))\n",
    "        Xi['mSTSBL'].append(list(mu_xi_STSBL.flatten()))\n",
    "        Xi['lSTSBL'].append(list(mu_xi_lSTSBL.flatten()))\n",
    "        Xi['dSTSBL'].append(list(mu_xi_dSTSBL.flatten()))\n",
    "\n",
    "        # Record errors\n",
    "        L1_err['ARD'].append(np.linalg.norm(xi - mu_xi_ARD, 1))\n",
    "        L1_err['ARDr'].append(np.linalg.norm(xi - mu_xi_ARDr, 1))\n",
    "        L1_err['ARDvi'].append(np.linalg.norm(xi - mu_xi_ARDvi, 1))\n",
    "        L1_err['mSTSBL'].append(np.linalg.norm(xi - mu_xi_STSBL, 1))\n",
    "        L1_err['lSTSBL'].append(np.linalg.norm(xi - mu_xi_lSTSBL, 1))\n",
    "        L1_err['dSTSBL'].append(np.linalg.norm(xi - mu_xi_dSTSBL, 1))\n",
    "\n",
    "        L2_err['ARD'].append(np.linalg.norm(xi - mu_xi_ARD, 2))\n",
    "        L2_err['ARDr'].append(np.linalg.norm(xi - mu_xi_ARDr, 2))\n",
    "        L2_err['ARDvi'].append(np.linalg.norm(xi - mu_xi_ARDvi, 2))\n",
    "        L2_err['mSTSBL'].append(np.linalg.norm(xi - mu_xi_STSBL, 2))\n",
    "        L2_err['lSTSBL'].append(np.linalg.norm(xi - mu_xi_lSTSBL, 2))\n",
    "        L2_err['dSTSBL'].append(np.linalg.norm(xi - mu_xi_dSTSBL, 2))\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_ARD,xi)\n",
    "        Added['ARD'].append(added)\n",
    "        Missed['ARD'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_ARDr,xi)\n",
    "        Added['ARDr'].append(added)\n",
    "        Missed['ARDr'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_ARDvi,xi)\n",
    "        Added['ARDvi'].append(added)\n",
    "        Missed['ARDvi'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_STSBL,xi)\n",
    "        Added['mSTSBL'].append(added)\n",
    "        Missed['mSTSBL'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_lSTSBL,xi)\n",
    "        Added['lSTSBL'].append(added)\n",
    "        Missed['lSTSBL'].append(missed)\n",
    "\n",
    "        _,added,missed = sparsity_err(mu_xi_dSTSBL,xi)\n",
    "        Added['dSTSBL'].append(added)\n",
    "        Missed['dSTSBL'].append(missed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "if retrain:\n",
    "    with open('./saved_results/L96_results.json','w') as f:\n",
    "        json.dump({'L2_err': L2_err, \\\n",
    "                   'L1_err': L1_err, \\\n",
    "                   'Added' : Added, \\\n",
    "                   'Missed' : Missed, \\\n",
    "                   'Xi' : Xi}, f)\n",
    "else:\n",
    "    with open('./saved_results/L96_results.json','r') as f:\n",
    "        results = json.load(f)\n",
    "    L2_err = results['L2_err']    \n",
    "    L1_err = results['L1_err']    \n",
    "    Added = results['Added']\n",
    "    Missed = results['Missed']\n",
    "    Xi = results['Xi']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([<matplotlib.axis.XTick at 0x7f39a132f7f0>,\n",
       "  <matplotlib.axis.XTick at 0x7f39a132f128>,\n",
       "  <matplotlib.axis.XTick at 0x7f39a1353550>,\n",
       "  <matplotlib.axis.XTick at 0x7f39a12c1550>,\n",
       "  <matplotlib.axis.XTick at 0x7f39a12c1a58>,\n",
       "  <matplotlib.axis.XTick at 0x7f39a12ca080>],\n",
       " <a list of 6 Text xticklabel objects>)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA24AAAEtCAYAAABwAjOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde7hcVX3/8feHBEES7kRUNEQooIZy81DLD5FYUCCVioIUDUKsbVSMFisIKpiAKGqLVguCURBQSosQQaqAWgm3UjQIEaOQKjcVkRPAmIQQMHx/f+w1ZJ9hZs6cM3tmz+Xzep55MrP3mrXWnjnfzF57XbYiAjMzMzMzM+teG5RdATMzMzMzM2vMDTczMzMzM7Mu54abmZmZmZlZl3PDzczMzMzMrMu54WZmZmZmZtbl3HAzMzMzMzPrcm64mZmZmZmZdTk33LqApI0lXSnpF5LulHSdpB3KrpdZP5L0UUn3SHpG0mFl18fMzMysGW64dY9zI+IVEbEHcDXw1bIrZNan/huYCdxYdkXMzMzMmuWGWwkkfVzSXZKOBIiIJyPiulyS/wVG7XGTNFtS1Hn8oV31N+sV1bEGEBG3RcSvxpiPY83MzMxK5YZbh0maCcwBbgHeUCfZ+4GrxpDtW4F9qh4HtlBNs57XZKyNlWPNrAmSvpoubHxuDO9ZJGlRE+nmS4qWKjgyv8qFmWlF5WnWLlUXEneusX9Gbv+BaduFku7veGUbSPWbX3Y9es3EsiswgP4B+AqwM/Ccq/6SPpL2HTCGPO+MiF+OpRKSNoqItWPd12reZh3UMNbGybFmNgpJzye7yAEwS9KHI+JPZdbJrA+tBN4BnFq1/Zi0b9Pctk8AX+hQvayN3OPWQZI2Al4PfBv4f8D/VO0/ATgcOCQiniiw3PnpysauaeGTVcBlTew7WNKtktZIWpEWUNml2bzNyjJarLWxXMeaGbwZ2Az4LvAC4OByq2PWlxYCR0tSZUO6aHI4cEU+YUT8KiLu6HD9rA3ccOusvwSeAZ4CtgJuqOyQ9E/A24DXR8RY58xMkDSx6lHru70qlfk3wOcb7ZN0MPAdYBXwt8B7gV2BmyVtN8a8zTqtbqy1yLFmNrpjgceB2cAash6AESQdJeluSWslLZX05loZSdpT0k2SnpT0W0mnAqqRbqKkj+TyfEjSWZI2rkq3g6TvSHpC0rCkLwAbFXDMZp32dWB74DW5bW8GJlDVcKseKpni5ROSfpVia7mkmyW9Jpfm7ZLukLQqXVC8S9K7q/LdX9J/S1opaXW6qLhrVZoJks6Q9LsUd4skTS/uYxgsHirZWfsCtwNHA5dGxBoASS8BzgLuBa5PF0/+FBFDTeZ7d41t3wHeWLXtixFRr6t8xD5Ji1N9DqkMcZF0K7AM+BDwT2PI26zTasYagKRTgPcAU4BdJZ0NDEXEw03k61gza0DSi8nmfS6IiGFJVwJvkbRlRDye0hwI/DtZ7HyILBa/AGwI3JPLaxvgh8DDZI3BtcCJwNQaRX8DOBT4DFkP+yvIhodNI+uBQNLzgO8DzwfeBzwCvBt4S2EfgFnnPEC2OvI7gJvStmOAb5FdCGzkJOCDwMeAO8l6yIfILnSSGnDfAL5IFnMbAC8HtqhkIOmvyS4kfofst7aS702SdouIX6dt84GPAp8DvpfK+fY4jtdww63TpgP3AX8HvLayMSJ+Q40riGPwZuA3Vdtq9dp9q0Eez+6TNAnYC/hUfl5CRNwn6RZg/zHmbdZpNWMNICLOAM4YZ76ONbPG3kF2kndxen0R2WiSvwXOS9tOI7sI8qaIeAZA0i/IVlS+J5fXB4FJwEER8WBK932yE9ZnSdov5X9sRFTK/YGkx4BvSNojIu4ka/ztAOwTEf+b3nsNcFdBx27WaRcDZ0n6ALAl2UWTQ5p43z7A96ouAl6de/6XwB8i4vjctu9V5fEF4IaIeFNlg6TryS5Efgg4XtKWZHG8ICJOqOQjaR3w6SbqaVU8VLKzXgy8jixY7hkt8Rj8LCIWVz1qLaDwuwZ55PdtSdaQrJX+YdIVmTHkbdZpjjWzchwD/F9E3Jpe/wB4KG1H0gRgb+DySqMNstt0APdX5bUP8L+VRltKt5qRJ5iQzaF7CrgiP4yZ9SealYs3+wC/rjTaUn7P4Lmi1ru+STbU91BgFtnvxn838b4fAzMlfVLSa1JvdPX+LSV9Q9IbJW2R3ylpJ2BH4JKqmHsCuJX1MffnZBdfqmPsP5o/RMtzw62ztiAbEvKRkspvtHxyft/j6fULa6R7IfDoGPM26zTHmlmHSdobeCWwUNIW6WRvU7JFFPZRtnT5NmRDIn9fI4vqbS9qMt0LgOeRDQ97Ovd4JO3feoz5mfWEiFgJXEnW030McEn+gkgDnwLmkc2Vvgl4VNLX0vBkIuIGspVhX0o2ymNY0g8k7Zbe/4L07/mMjLmnyaYO5GMOnhtjjrlxcsOtszYAvpCGRnatdEXzduCt6eooAJK2J1uhr6iFHszaxbFm1nnHpn9PIrsoUXnMTduPAZaTndxtW+P91dt+12S6R4EnyXryaj2+PMb8zHrJxcBfk/VuXTxKWgAi4umI+ExE/DlZ4+qDZHNBz8mluTwi9icbGfLmlO7atCBX5aLiR6gdc4em/ZURItUx5pgbJ89x6xBJx5IF1c/SCdrngbMjYlkB2e9RuUpSZXEL9845lWzC6X9J+hIwmWxewgqyhVTMupJjzazz0lCro4DbgJNrJPk86+859WPgCEnzc3PcXk22kEh+/tqtwImSXlpZ6CDNCz2Uka4layxuHhGNhondCrxT0l/m5rhtABw5lmM16zLfJxuK+IeIWDrWN6eFub4qaSbZisbV+1eR/T7tQDavbWuyuaj3A9MjotFctZ8Cq8li7Ie57UeNtZ6WccOtAyRtQtbl/EayVa5+BnyloBNJyMY41zKF7OrmmEXEtWnFoHlk/yE8BSwCPhwRD40nT7N2c6yZlaYyPOpDEbGoeqekLwPnAjPI/ta/B1yZtk8hu1hRvbLr54HjyBYzmM/6VSXX5BNFxCJJlwKXS/oc8COy24FMA2YCJ6X/Ay4ia1QulPRRsqGU7yFbUc+sJ0XEOrIFgJom6SpgCfATsl7xPcnmin457T+drFfserI5qi8BPgDcGRHDKc37gKvSRZvLyH4DtyUbLfJgRHwuIv4g6fPAxyStJIv7vYF3tXTQA0wRni5hZmZm45dOBF8HvDAinqixf3OyYVOXRcRsSW8jWyZ8GvBL4BTgHwEiYkbufXuRXeXfm2x41nlkF50/HhH5Gw9vALyfbCXZXcgaefcD1wGfjIgVKd0OwNlkDcjVZLcl+HnK92URcX/LH4ZZG0maDXwN2KnO4lhImkHW6Hp9RPxA0oXAjIiYlvZ/iOwi507AJsCDwKVksfJ0upj4AWA3skWyHiFrdJ2av6AoaR+yWwrsS3abjYfJVof9QmWBojTy5TTg74HNyXrljwOWAqdFxPwCPpaB4YabmZmZmZlZl/PiJGZmZmZmZl3ODTczMzMzM7Mu54abmZmZmZlZl3PDzczMzMzMrMuVcjuAbbbZJqZNm1ZG0WZtc/vtty+PiCll1yPPsWb9qNtizXFm/chxZtZ+Y42zUhpu06ZNY/HixWUUbdY2kh4YPVVnOdasH3VbrDnOrB85zszab6xx5qGSZmZmZmZmXc4NN7M2kjRX0mJJa9MNMCvbZ0lalXs8ISkkvapOPoskPZlLf0/HDsLMzMzMSueGm1l7PQScAVyQ3xgRl0TE5MoDOA64F/hJg7zm5t6zS/uqbGZmZmbdppQ5bmaDIiIWAkgaAl7SIOmxwMURER2pmJmZmZn1FPe4mZVM0vbAa4GLR0l6pqTlkm6RNKP9NTMzMzOzbtFUw63ePJ20bxNJX0onlCsk3diWmpr1r2OAmyLivgZpTgJ2ALYDFgBXS9qxVkJJc1K8Lh4eHi6+tmZmNnAkbSTpfEkPSFop6Q5JhzRI/0FJD6dzwwskbZTbN03S9Wl+992SDuzMUZj1tmZ73GrO00kWAFsBr0j/frCYqpkNjGOAixoliIjbImJlRKyNiIuAW4CZddIuiIihiBiaMqVrbsFjZma9bSLwa2B/YHPgVOAySdOqE0o6CDgZOACYRnbh8bRckkuBO4CtgY8Bl0vyD5bZKJpquEXEwoi4Eng0v13SLsDfAHMiYjgi1kXE7W2op1lfkrQv8GLg8jG+NQAVXyMzM7PniojVETE/Iu6PiGci4r+A+4BaqyEfC5wfEUsj4nHgE8BsAEk7A3sB8yJiTURcAdwFHN6RAzHrYa3OcXs18ABwWhoqeZckB571DEmjPlrMf6KkjYEJwARJG0vKLwp0LHBFRKxskMcWkg6qvFfSLLI5cde1VDmzDmpnnPUafxbWDyRtC+wMLK2xezqwJPd6CbCtpK3TvnurfveWpO3VZXjov3Wddp87NtJqw+0lwK7ACrJeg7nARZJeUZ3QwWfdKCJGPOpta8EpwBqyISNHp+enAKQG3ZHUGCYp6aOSrkkvNyQbqjwMLAfeDxwWEb6Xm/WM0WJtkLT5/xyztpO0IXAJcFFE3F0jyWSyc8OKyvNNa+yr7N+0OhMP/bdu1IFzx7pavR3AGuBp4IyI+BNwg6TrgTcAv8gnjIgFZPPhGBoa8i+TDYSImA/Mr7PvSWCLOvs+lXs+DOzdhuqZmZmNiaQNgK8DT5FdsK9lFbBZ7nXl+coa+yr76448MbNMqz1uPy2kFmZmZmbW1ZSNATsf2BY4PCKerpN0KbB77vXuwO8j4tG0bwdJm1btrzXk0sxymr0dQL15OjcCDwIfSWn2BWbguTdmZmZm/eZcslXED42INQ3SXQy8S9IrJW1JNkXgQoCIWAbcCcxL55NvBnYDrmhrzc36QLM9bjXn6aQrLW8iW5Z8BfAV4Jg6453NzMzMrAdJ2h54N7AH8LCkVekxS9LU9HwqQERcC3wWuJ5sEbsHgHm57I4ChoDHgU8DR6RpAWbWQFNz3EaZp7MU2Ke4KpmZmZlZN4mIB2h8G5rJVek/B3yuTl73k43QMrMxaHWOm5mZmZmZmbWZG25mZmZmZmZdzg03MzMzMzOzLueGm5mZmZmZWZdzw83MzMzMzKzLNbWqpJmZmY2dTmu0CB8wf/Q0MS+Kq5CZmfUs97iZmZmZmZl1Ofe4mZmZtdl4es1G7a0zM7OB4h43MzMzMzOzLueGm5mZmZmZWZdzw83MzMzMzKzLueFmZmZmZmbW5dxwMzMzMzMz63JuuJmZmZmZmXU5N9zMzMzMrCFJcyUtlrRW0oUN0p0naVXusVbSytz+RZKezO2/pyMHYNYHmmq4NROskuZJCkkHFlpDsx5WL3YkTUvxkv9xO7VBPtMkXS/pCUl3O87MzKzDHgLOAC5olCgi3hMRkysP4FLgm1XJ5ubS7NKm+pr1nWZvwF0J1oOA51fvlLQjcATwu+KqZtYXGsYOsEVE/KmJfC4FbgVmpsflknaKiOHCampm7fE1mHH9jBGbjjzySI477jieeOIJZs6c+Zy3zJ49uzN1M2tSRCwEkDQEvKSZ90iaBBwOvLGNVTMbGE31uEXEwoi4Eni0TpKzgZOAp4qqmFk/aCJ2RiVpZ2AvYF5ErImIK4C7yH4MzczMutXhwDBwY9X2MyUtl3SLpBmdr5ZZb2q2x60uSW8FnoqI70pqlG4OMAdg6tSprRZr1i8ekBTA94ETI2J5jTTTgXsjYmVu25K03cy63Tth0bxFNXdtsskmLFpUe987T3tn++pk1hnHAhdHROS2nQT8nOxi/1HA1ZL2iIhfVb/Z545mI7W0OImkycCngONHSxsRCyJiKCKGpkyZ0kqxZv1gObA3sD3wKmBT4JI6aScDK6q2rUjveQ5Jc9K8usXDwx5JaWZmnSfppcD+wMX57RFxW0SsjIi1EXERcAvZFIDn8Lmj2Uitrip5GvD1iLiviMqYDYqIWBURiyPiTxHxe2Au8AZJm9VIvgqo3r4ZsLJGWv/QmZlZNzgG+J+IuHeUdAHUH7JlZs9qteF2APABSQ9Lehh4KXCZpJNar5rZQKkMI6n147UU2EFSvodt97TdzMys7SRNlLQxMAGYIGljSY2m3BwDXFiVxxaSDqq8V9Is4LXAdW2ruFkfafZ2APWC9QBgV2CP9HgIeDdwTpvqa9ZT6sWOpFdL2kXSBpK2Br4ILIqI6iGRRMQy4E5gXnr/m4HdgCs6eSxmZjbQTgHWACcDR6fnp0iamm5p8+wkNEn7kK08WX0bgA3JVloeJpsy8H7gsIjwvdzMmtDs4iSnAPNyr48GTouI+flEktYBj0fEqmKqZ9bzasYOcA/Z/NAXAH8kW5zkbZVEks6D7H44adNRZFcuHwceBI7wrQDMzKxT0jnf/Dq7J1elvRWYVCOPYbL53WY2Dk013EYJ1ny6aa1Vx6y/jBI7lzZ433uqXt8PzCioWmZmZmbWY1qd42ZmZmZmZmZt5oabmZmZmZlZl3PDzczMzMzMrMu54WZmZmZmZtblml1V0szMzMZJp/n+wmZm1hr3uJmZmZmZmXU597iZmZm1ScyLhvslEdE4jZmZGbjHzczMzMzMrOu54WZmZmZmZtbl3HAzMzMzMzPrcm64mZmZmZmZdTk33MzMzMzMzLqcG25mZmZmZmZdzg03MzMzMzOzLueGm5mZmZk1JGmupMWS1kq6sEG62ZLWSVqVe8zI7Z8m6XpJT0i6W9KBnai/WT9oquFWL1gl/aWk70t6TNKwpG9KelHbamtmZmZmZXgIOAO4oIm0t0bE5NxjUW7fpcAdwNbAx4DLJU0pvLZmfajZHrd6wbolsACYBmwPrAS+VlTlzMzMOknS5pJ+lHoJdi27PmbdIiIWRsSVwKPjzUPSzsBewLyIWBMRVwB3AYcXVE2zvjaxmUQRsRBA0hDwktz2a/LpJJ0N3FBkBc3MzDroCeCvgX8uuyJmPWxPScuBx4CvA2dGxJ+A6cC9EbEyl3ZJ2m5moyh6jttrgaW1dkiak4ZbLh4eHi64WDMzs9ZFxNMR4R8ps/G7EdgVeAFZT9rbgBPTvsnAiqr0K4BNa2Xkc0ezkQpruEnaDfg464NzhIhYEBFDETE0ZYqHMttgKGp+qKRFkp7MTfS+pyMHYNajGi2kIGkrSd+StFrSA5LeXlI1zfpORNwbEfdFxDMRcRdwOnBE2r0K2KzqLZuRTbWplZfPHc1yCmm4Sfoz4BrgHyPipiLyNOsTRc4PnZub6L1L0RU16zONFlI4B3gK2BaYBZwryUO1zNojAKXnS4EdJOV72HanzmgtMxup5YabpO2BHwCfiIivt14ls/5RbzJ3RFwTEd+MiD9GxBPA2cC+pVTSrA/Viz1Jk8iGb50aEasi4mbg28A7SqimWc+QNFHSxsAEYIKkjSU9Z60ESYdI2jY9fzlwKnAVQEQsA+4E5qX3vxnYDbiiU8dh1suavR1AzWCVtB3wQ+CciDivnRU163N154fmnClpuaRb8vfEqeY5AWYN7QysSyeQFSMWR5D0XeANwFckza6ViePMBtApwBrgZODo9PwUSVPTEP6pKd0BwE8lrQa+CywEPpXL5yhgCHgc+DRwhOeVmjWnqVUlyYJ1Xu710cBpZN3fO5BdOXl2f0RMLqyGZn0uNz/0TQ2SnQT8nGx411HA1ZL2iIhfVSeMiAVkwzAZGhqK4mts1tNGXRwhImaOlonjzAZNRMwH5tfZPTmX7gTghAb53A/MKK5mZoOjqR63iJgfEap6zI+I09Lz/E0W3Wgza1Kz80Mj4raIWBkRayPiIuAWYNSTSzN7jjEtjmBmZtYtir4dgJk1qcX5ofnJ3mbWvGXAREk75bZ5cQQzM+t6briZtVER80MlbSHpoNx7Z5HNibuu/Udg1pvqxV5ErCabc3O6pEmS9iUbpuzFtczMrKu54WbWXjUncwN/z/r5oZV7s62qvEnSRyVdk15uSLas+TCwHHg/cFhE+F5uZvXViz2A44DnA48AlwLvjQj3uJmZWVdrdnESMxuHUSZzn9bgfZ/KPR8G9i60YmZ9rlHsRcRjwGGdrI+ZmVmr3ONmZmZmZmbW5dxwMzMzMzMz63JuuJmZmZmZmXU5N9zMzMzMzMy6nBtuZmZmZmZmXc4NNzMzMzMzsy7nhpuZmZmZmVmXc8PNzMzMzMysy7nhZmZmZmZm1uXccDMzMzMzM+tyTTXcJM2VtFjSWkkXVu07QNLdkp6QdL2k7dtSUzMzMzMrRaNzwap0x0q6XdIfJf1G0mclTcztXyTpSUmr0uOejhyAWR9otsftIeAM4IL8RknbAAuBU4GtgMXAfxZZQTMzMzMrXc1zwRo2AY4HtgFeDRwAnFCVZm5ETE6PXQqvqVmfmjh6EoiIhQCShoCX5Ha9BVgaEd9M++cDyyW9PCLuLriuZmZmZlaCBueC1enOzb38raRLgNe1uXpmA6HVOW7TgSWVFxGxGvhV2m5mZmZmg+21wNKqbWdKWi7pFkkz6r1R0pw0PHPx8PBwWytp1gtabbhNBlZUbVsBbFqd0MFnZmZmNjgkvRMYAv4lt/kkYAdgO2ABcLWkHWu9PyIWRMRQRAxNmTKl7fU163atNtxWAZtVbdsMWFmd0MFnZmZmNhgkHQZ8GjgkIpZXtkfEbRGxMiLWRsRFwC3AzLLqadZLWm24LQV2r7yQNAnYked2iZuZmZnZAJB0MPAV4NCIuGuU5AGo/bUy633N3g5goqSNgQnABEkbp6VdvwXsKunwtP/jwE+9MImZmZlZ/2hwLlid7q+AS4DDI+JHVfu2kHRQ5b2SZpHNgbuuE8dg1uua7XE7BVgDnAwcnZ6fEhHDwOHAJ4HHyZZ9PaoN9TQzMzOz8tQ8F5Q0Nd2PbWpKdyqwOfDd3L3arkn7NiS7pcAwsBx4P3BYRPhebmZNaPZ2APOB+XX2/QB4eXFVMusfkuYCs4E/By6NiNm5fQcA5wBTgduA2RHxQJ18pgFfI7s48iDZPXB+0Maqm5mZPavRuSDZYnWVdHWX/k8X/PcutGJmA6TVOW5m1lhRN6+/FLgD2Br4GHC5JK/yY2ZmZjYg3HAza6OIWBgRVwKPVu169ub1EfEk2VXM3SU9p/da0s7AXsC8iFgTEVcAd5ENUzYzMzOzAeCGm1k5xnLz+unAvRGRv83Gkjppfc9EMzMzsz7khptZOZq+ef0Y0/qeiWZmZmZ9yA03s3I0ffP6MaY1MzMzsz7khptZOcZy8/qlwA6S8j1su9dJa2ZmZmZ9yA03szYq4ub1EbEMuBOYl97/ZmA34IrOHYmZmZmZlckNN7P2GtfN6yWdJ+m8XD5HAUMp7aeBI1IeZmZmZjYAmroBt5mNz3hvXh8R76l6fT8wo9DKmZmZmVnPcI+bmZmZmZlZl3PDzczMzMzMrMu54WZmZmZmZtbl3HAzMzMzMzPrcm64mZmZmZmZdTk33MzMzMzMzLpcIQ03SdMkfVfS45IelnR2usmwmZmZmfU4SXMlLZa0VtKFo6T9YDofXCHpAkkb5fZNk3S9pCck3S3pwLZX3qxPFNXj9iXgEeBFwB7A/sBxBeVtZmZmZuV6CDgDuKBRIkkHAScDBwDTgB2A03JJLgXuALYGPgZcLmlKG+pr1neK6hV7GXB2RDwJPCzpWmB6QXmbmZmZWYkiYiGApCHgJQ2SHgucHxFLU/pPAJcAJ0vaGdgLeENErAGukHQ8cDhwXjvr3w5atKiQfGLGjELy6Wdlftajln399aOmKeo7LqrH7QvAUZI2kbQdcAhwbUF5m5mZmVlvmA4syb1eAmwraeu0796IWFm1v+bFfklz0vDMxcPDw2OqhKSGDyuOP+vOKarH7QbgH4A/AhOAi4Ar8wkkzQHmAEydOrWgYs3MzMysi0wGVuReV55vWmNfZf92tTKKiAXAAoChoaEYSyUiRiaX9JxtRRlvb0pRvUhlG5TPui29dWPUco+bpA2A64CFwCRgG2BL4DP5dBGxICKGImJoyhQPZTYzMzPrQ6uAzXKvK89X1thX2b8SMxtVEUMltwJeSjbHbW1EPAp8DZhZQN5mZmZ9pXoIkYcVWZ9ZCuyee7078Pt0frgU2EHSplX7l3awfmY9q+WGW0QsB+4D3itpoqQtyCamLmn8TjMzs8ETEQ0fZt0oneNtTDYlZoKkjevc+uli4F2SXilpS+AU4EKAiFgG3AnMS+9/M7AbcEVHDsKsxxW1OMlbgIOBYeCXwJ+ADxaUt5l1wGiTi90TYGY20E4B1pAt9X90en6KpKmSVkmaChAR1wKfBa4HHkiPebl8jgKGgMeBTwNHRMTYVh4xG1CFLE4SEXcCM4rIy8zK0cnJxWZm1lsiYj4wv87uyVVpPwd8rk4+9+NzRrNxKarHzczMzMzMzNrEDTczMzMzM7Mu54abWUnSnID8Y52kf6uTdnban08/o8NVNjMzM7OSFHUDbjMbo4h4dk6ApEnA74FvNnjLrRHxmrZXzMzMzMy6jnvczLrDEcAjwE1lV8TMzMzMuo8bbmbd4Vjg4mi8jOOekpZLWibp1Dr3zzEzMzOzPuQTP7OSpXvf7A+8q0GyG4Fdye6HMx34T7L7JZ5ZI785wByAqVOnFl1dMzMzMyuBe9zMyncMcHNE3FcvQUTcGxH3RcQzEXEXcDrZ8MpaaRdExFBEDE2ZMqVNVTYzMzOzTnLDzax8xwAXjfE9AagNdTEzMzOzLuSGm1mJJP0/YDsaryaJpEMkbZuevxw4Fbiq/TU0MzMzs27ghptZuY4FFkbEyvxGSVPTvdoqk9QOAH4qaTXwXWAh8KnOVtVsMEj6pKSbJF0uaZOy62NmZgZenMSsVBHx7jrbHwQm516fAJzQqXqZDSpJuwI7RsR+kt4D/B1wdsnVMjMzc4+bmZlZzn7ANen5NYBvem9mZl3BDTczM+s7kuZKWixpraQLq/ZtJelbklZLekDS23O7twRWpOcrgK06VGUzM7OG3HAzM7N+9BBwBnBBjX3nAE8B2wKzgHMlTU/7Hgc2T883Bx5rcz3NesYoFz3y6a5J87Qrj6ck3ZXbf7+kNbn93+vcUZj1rsIabpKOkvSLFMy/krRfUXmbmZmNRUQsjIgrgUfz2yVNAg4HTo2IVRFxM2RpPBQAACAASURBVPBt4B0pyc3AQen5QcAtHaqyWS9odNHjWRFxSERMrjyA/+G5qycfmkvzhrbX3KwPFLI4iaTXA58B/hb4EfCiIvI1MzMr2M7AuohYltu2BNgfICLuSj0JNwGPkN1n8TkkzQHmAEydOrVWkq4jNb71Y0R0qCadMdrxQvuOuazPup3HnLvosWtErAJullS56HFyg/dNI5s7+s5xFWxmzypqVcnTgNMj4n/T698WlK+ZmVmRJrN+DlvFCmDTyouI+MhomUTEAmABwNDQUE+0eKpP2CX1XWMtr8zjLavsNpfb8KJHA8cAN0XEfVXbL5G0AXAHcGJELCmqomb9quWhkpImAEPAFEm/lPQbSWdLen5Vujlpovji4eHhVos1MzMbj1XAZlXbNgNW1khrZuuNetGjjmOAC6u2zQKmAdsD1wPXSdqi+o0+dzQbqYg5btsCGwJHkHWF7wHsCZySTxQRCyJiKCKGpkyZUkCxZmZmY7YMmChpp9y23YGlJdXHrFeM+aKHpNcALwQuz2+PiFsiYk1EPBERZwJ/IDuHpCqdzx3NcopouK1J//5bRPwuIpYDnwNmFpC3mZnZmEmaKGljYAIwQdLGkiZGxGpgIXC6pEmS9gXeBHy9zPqa9YDxXPQ4FliY5sQ1EsDoE/TMBlzLDbeIeBz4DVnQmZmZdYNTyC4sngwcnZ5XRoIcBzyfbPGRS4H3RoR73MwaGOtFjzRl5q1UDZOUNFXSvpKely6onAhsg1dwNRtVUYuTfA14v6RrgaeB44H/KihvMzOzMYmI+cD8OvseAw7rZH3M+sRxZPdGfITsVhvvjYil6RZQ16Sl/ysOI5sDd31VHpsC5wI7Ak8CdwKHRMSjmFlDRTXcPkF2tWQZWRBeBnyyoLzNrI+VuWR3Wco65kH8rM2sOPUuekTETWSLl+S3XUrWo12ddimwW7vqaNbPCmm4RcTTZFdhjisiPzMbHIO2RDn07VLhZmZm1kZFLE5iZmZmZmZmbeSGm5mZmZmZWZcrao5bW402L8NDfczMzMzMrJ/1RMPN8zLMzMzMzGyQeaikmZmZmZlZl+uJHreyeOlsMzMzMzPrBu5xayAiRjzqbesnkho+yiq3nWWXSdIiSU9KWpUe99RJJ0mfkfRoenxW/fqhmJmZmdlzuMfNRvD9pUoxNyK+OkqaOWQ3Pd0dCOD7wL3AeW2um5mZmZl1Afe4mfWGY4GzIuI3EfFb4CxgdrlVMjMzM7NOccPNrHxnSlou6RZJM+qkmQ4syb1ekraNjzTqI5pJZ2ZmZmYd4YabWblOAnYAtgMWAFdL2rFGusnAitzrFcDkWvPcJM2RtFjS4uHh4XbU2czMzMw6zA03sxJFxG0RsTIi1kbERcAtwMwaSVcBm+VebwasihoTASNiQUQMRcTQlClTRqvA+B5mZmZm1lFuuJl1lwBqjUFcSrYwScXuaZuZmZmZDQA33MxKImkLSQdJ2ljSREmzgNcC19VIfjHwT5K2k/Ri4EPAhR2srpmZmZmVqNCGm6Sd0j2pvlFkvmZ9akPgDGAYWA68HzgsIu6RtJ+kVbm0XwauBu4CfgZ8J20zMzPrCElbSfqWpNWSHpD09jrp5kt6OneP0lWSdsjt30PS7ZKeSP/u0bmjMOtdRd/H7RzgxwXnadaXImIY2LvOvpvIFiSpvA7gw+lhZmZWhnOAp4BtgT2A70haEhG1hu7/Z0QcXb1R0vOAq4B/Bb4EvBu4StJOEfFU+6pu1vsK63GTdBTwB+C/i8rTzMzMzMonaRJwOHBqRKyKiJuBbwPvGGNWM8g6Dv41Lcz1RbK53X9VZH3N+lEhDTdJmwGnk827MTMzM7P+sjOwLiKW5bY1uqfooZIek7RU0ntz26cDP61aFfmnDfIxs6SoHrdPAOdHxK/rJfC9pczMzMx6VvX9REmvN62R9jLgFcAU4B+Aj0t621jz8bmj2UgtN9zShNIDgc83Sjeme0uZmZmZWTepvp8o6fXK6oQR8fOIeCgi1kXE/wBfAI4YRz4+dzTLKaLHbQYwDXhQ0sPACcDhkn5SQN5mZmZmVr5lwERJO+W2NXtP0fw9SpcCu0nK37N0tybzMRtoRTTcFgA7kq0utAdwHtlS5QcVkLeZmfUaadRHjJbOzLpKRKwGFgKnS5okaV/gTcDXq9NKepOkLZX5C+ADZCtJAiwC1gEfkLSRpLlp+w/bfhBmPa7lhltEPBERD1ceZF3gT6alzs3MzMysPxwHPB94BLgUeG9ELK1x79GjgF+SDX+8GPhMRFwEkJb8Pww4hmw18r8ju4epbwVgNoqi7+NGRMwvOk8z6yOj9KZEE2kYsRiZda3xfE/ubTPrWhHxGFmjq3p79b1H31adpir9HcCrCq+gWZ8r7D5uZmZmZmZm1h6F97iZmTXFvTFmZmZmTXOPm5mZmZmZWZdzw83MzMzMzKzLueFmZmZmZmbW5TzHzQbK6FOkwgsaDhCNtsKlv2wzMzPrEm64mdnAqm6YSXJjzczMzLpS1zXcmls0rnGvSK+fd43WCwDuCWiVFzQ0MzMzs17SdQ03cy+AmZmZmZmN1LUNt073iBTR0we939tnZmZmZmbdx6tKmpmZmZmZdTk33KpEjO9hNlaSNpJ0vqQHJK2UdIekQ+qknS1pnaRVuceMDlfZzMzMzErStUMlzQbARODXwP7Ag8BM4DJJfx4R99dIf2tEvKaD9TMzMzOzLuGGm1lJImI1MD+36b8k3Qe8Cri/jDqZmZmZWXfyUEmzLiFpW2BnYGmdJHtKWi5pmaRTJfnCi5mZmdmAaLnhNpZ5OmZWm6QNgUuAiyLi7hpJbgR2BV4AHA68DTixTl5zJC2WtHh4eLhdVTYzswEjaStJ35K0Op33vb1OuhMl/SydF94n6cSq/fdLWpObs/29zhyBWW8rosctP09nc+BUsnk60wrI26zvSdoA+DrwFDC3VpqIuDci7ouIZyLiLuB04Ig6aRdExFBEDE2ZMqVt9TYzs4FzDtlv1bbALOBcSdNrpBNwDLAlcDAwV9JRVWkOjYjJ6fGGdlbarF+0PNTK83TMxk+SgPPJfgRnRsTTTb41yH4YzczM2k7SJLIRH7tGxCrgZknfBt4BnJxPGxGfzb28R9JVwL7Af3SqvtYeWrRo9ETXX98wXcyYUVR1Bk7hc9yamKdjZuudC7yC7MrjmnqJJB2SYgtJLyfr2b6qM1W0niU1fEQTaczMkp2BdRGxLLdtCVCrx+1Z6QLlfjz3vPASScOSvidp9zrv9dB/s5xCFzdoNE9H0hxgDsDUqVOLLNasJ0naHng3sBZ4WOtPkt8N3AT8HHhlRDwIHABcKGky8HvgG8CnOl5pMzMbVJOBFVXbVgCbjvK++WQdBV/LbZsF/IRs5Mg/AtdJenlE/CH/xohYACwAGBoa8l1zu8h4es2a6q2zhgpruI02T8fBZzZSRDxA4+GOk3NpTwBOaHulrD/FOP7LdW+bmY20CtisattmwMp6b5A0l2yu234RsbayPSJuySU7U9KxZL1yVxdXXbP+U0jDrYV5OmZWkhkAVVfMjjzySI477jieeOIJZs6c+Zz3zJ49m9mzZ7Mc2KYDdTQzs66xDJgoaaeI+L+0bXfqTI2R9Hdkc99eGxG/GSVvz9s2a0JRc9yamqdjZmZmZr0nLUa3EDhd0iRJ+wJvIhttNYKkWWTD+V8fEfdW7ZsqaV9Jz5O0cbpVwDbALdX5mNlILfe4NZqnExGXtJq/mbXHIoA648032WQTFjUYi+7etjFoYsjhs4uENEzkEeZmVrrjgAuAR4BHgfdGxFJJ+wHXRERliP8ZwNbAj3Pnhd+IiPeQzYk7F9gReBK4EzgkIh7t3GGY9aYibgcw2jwdMzMzM+txEfEYcFiN7Tcxcl72yxrksRTYrS0VNOtzha4qadaM5tY8iIbp3PlgPWe8f7ReJMTMzMxww22guQFlZmZmZtYb3HCz0niFcuuoIuaa+UqFmZmZlcQNty4w+vlk414vaO180g0oMzMzM7Pu1qUNtxnVt5dq6v5SMLsDdTOznuYrFWZmZtaDurThNph8PmlmZmZmZrV0acNtUb3bSzW8v9Q739m2Cpn1J7f8zczMzHrCBmVXwMzMzMzMzBrr0h43M2urJsblSiK8iqKZmZlZV3CPm5mZmZmZWZdzj5uZDYwZQPWStc2sWDu7A3UzMzMza6RrG26Dt2aCb4FgZmZmZma1dW3DrRzjazxlDajlwDYdqKNZ75sBpfR8LQLqLVnbaMXaVpesnQHjO97Zs1v+n2XcZbdQZpnlmpmZ9auua7g1sxZCfy6aUNYtEMrq6Suzh9G9mwNt8LrzzczMrA90XcOtXONrPGVa723z+aQNikUAN9wwcuMNN8D73scmlf3VqtP3kEUwvp4+Wv+fZdxlt3hVqKxyWyFpK+B84A1kwyg+EhH/XlqFzLpMszEiScCngb9Pm84HTop01V3SHmnbK4BfAO+KiDvbfwRmva2Qhls//dgNXuNpUb3zZ6hzCl3M+XNZ5ZZd9khF/QhaE8q+BUKZ/7mUVXbv/Yd6DvAUsC2wB/AdSUsiYmm51TLrGs3GyBzgMGB3IIDvA/cC50l6HnAV8K/Al4B3A1dJ2ikinurMYZj1pqJuB5AP5FnAuZKmF5R334to/AA1kcZ6VLOxk/8R3A14I9mPXe8Z5Y9ZTaQxK5qkScDhwKkRsSoibga+Dbyj3JqZdYcxxsixwFkR8ZuI+C1wFuvnG8wg6zj414hYGxFfBAT8VZsPwazntdzjlgvkXSNiFXCzpEogn9xq/p0yiHPryjrmMj/r0bLs5Hc8xth59kcwvfcs4B+A8zpSWRu/Mnv6RsmzrHLbWvb47Qysi4hluW1LgP1Lqk/LNK1Gj+d04C/ILhdtXyPNHsCeEPPG/92MWu4lNd6UymU1xGfbVPZH6+zvxDGX9Vk3KHccn/VYYmR62pdPNz2376dVI0Z+mrZfO5YK5WmPPZ67ccYMOOwwePJJ2H3356Y5+GA4+GCieuL7WBx/PM/5FvLlnlzjVDiVy4oVsPnm4y561GNuUHZpx9yKEj/r0o65ShFDJZsKZElzyHoNmDp16pgKUI3hNvlt7frxH63cdpVdVrnNlO3PulBF/QiOMN5YG8TvwH/v/f9Zj8FkYEXVthXAprUSt/KbZtajxhIj1WlXAJPTsP+m83GcmY2kVn8sJe0HfDMiXpjb9g/ArIiYUes9Q0NDsXjx4pbKNes2km6PiKExpG86diStA6ZHxN3p9U7AMmCDRvPcHGvWj8Yaa03muSdwS0Rsktv2IWBGRBza6L2OM+tH1XE2lhiRtAJ4fUT8KL1+FbAoIjaV9MG0b2Yu/dVp/1n16uM4s3401t+zIua4rQI2q9q2GbCygLzN+tlYYqc67WbAKi9OYlaYZcDEdFGkYnfAC5OYZcYSI0vTvlrplgK7aWQ3/G518jGznCIabv6xMxufon4EzaxFEbEaWAicLmmSpH2BNwFfL7dmZt1hjDFyMfBPkraT9GLgQ8CFad8iYB3wAUkbSZqbtv+wnfU36wctN9z8Y2c2PgX+CJpZMY4Dng88AlwKvNe3AjAboWaMSNpP0qpcui8DVwN3AT8DvpO2kZb8Pww4BvgD8HfAYb4VgNnoiroB93HABWSB/Cj+sTNrVs3YSfPfromIySndl4EdyH4EAb6atplZQSLiMbITSjOroV6MRMRNZIuOVF4H8OH0qJXPHcCr2lRNs75VSMPNP3Zm41PUj6CZmZmZ9beibsBtZmZmZmZmbeKGm5mZmZmZWZdzw83MzMzMzKzLtXwD7nEVKg0DD7SQxTbA8oKq43K7s+xeLHf7iJhSZGVa1WKs9eJ30KtlD1q5rZbdVbHWo3FWZtmDVm6ZZTvO1uvF76AXyy2z7F4sd0xxVkrDrVWSFo/lLuMut/fKHrRyu9EgfgeDdsyD+Fl3m0H8Dgat3DLLdpytN2jfgf/e+7NcD5U0MzMzMzPrcm64mZmZmZmZdblebbgtcLl9X/aglduNBvE7GLRjHsTPutsM4ncwaOWWWbbjbL1B+w78996H5fbkHDczMzMzM7NB0qs9bmZmZmZmZgPDDTczMzMzM7Mu54ZbD5Oksutg1u8cZ2ad4Vgzaz/HWW9zw61FkiZUve5IQEjaBVgoaWInyqsq+28kTe50ud1G0l5l12FQOM4Gl+Osc8qKs1RWKbHmOFvPsdY5/k0bXK3GmRtuLZCkiFgnaUNJn5S0UXRutZd7gf+KiD91qDwAJF0LHBMRqzpY5n9Imtmp8poh6WrgE2XXYxA4zjpablfFmuOsc0qOMygh1hxn6znWOse/aR0rsy/jrGdXlUwt1p2A1cBDEfGTtF2dCABJEyJiXXp+ETAhIo5ud7mpvImVoJO0AfCvwEfbHRCSLgNeGhH7pNdt/6wlXQW8MCJe3c5yxkLSQmDLiHhd2XVpN8fZYMRZKqerYm2Q4gzKjbUy4yyV2fFYc5ytN0ix5t+0wfhN6+c463hXaREkfRXYAXgh8DzgGUkLIuJfIiI68UeRrpZsAJwH/J7Ugu5Q8OfzFzAL2EXSWyJidTsKlPQN4HBgRjvyr1PmxcDOEfGK9HrDiHi6U+XXqdMlwGHA9un1RhGxtsw6tYvjbDDiLJXbVbE2SHEG5cdayXEGHY41x9mIOg1MrJUdZ1B6rA3Eb1q/x1nPDZWUdCUwHXgnsBtwENmN7z4h6SMAHexy3jyVfzzZfwTQgc80Bb4kfSxduXkRMA24StKkostLVy5eAVwBvEvSgakeIbVnXLaky8mC/feSdk6bO9q1Xy19DnsCXwPOlrRTRKxN/wn3FcfZYMRZKrerYm2Q4gy6KtZKiTPobKw5ztYbpFjrojgD/6b53LGVOIuInnkAlwN31Nj+fOC9wG+BmW0sfy/gb4E3Aq9K214A/BT4QS7dhA58Fi8CVgBz0usNgbuBHwCTCiznK8BP0vOJZAF4KTAjl0YFH9s1wK3ALsCZwE3A3u0oawx1+mLuc9gFuBC4FpiWtm1QRr3adKyOs/Vl9G2cpTy7KtYGKc7S8ZQWa90UZ6mctsea42xEnQYm1vybNqIuffubNihx1vGDaOHgd0pfwjcrXwAwMbf/xcD1wD+1qfyvAj8Efg78Evg/4MS0bwpwTzsDsNaXCxwLfBt4WXq9IXAX8GNgkwLKnAJsV7Vtm/QdtCUAgf2BL+de7042DvvGsgIQeClwRP67BfYGLui3HzrH2WDEWcqrq2JtkOIsHUdpsVZ2nNX7HtsZa46zEXUamFgrM85S/v5NC587Vr7bIuKsl7rDfwWcStbleYWkLSPiT0pLmkbEQ8Aasi+rUHW62L8MnJ66nIeB/YBtJf041WddkXWIiGdSF/ef5TbfBWxB9sdBZGN4X0U2jnmbVsqTdAFwCbBlej1B2aTa5cB7yK6gvEfS/kWNzZb0DuBJYG56rYhYQjac4SfAWZL2jmjvkJaqOp1H9r3fmV5vGBHrIuLHwLnA74DzJE1L31EvxVQtjrM+j7NUTlfF2gDGGZQUa90QZynPjsWa42xEnQYt1vyb1ue/aQMXZ51sfY6zxToZ2CbXWn0dcBlZ1/eWafvz0r9fJ1tutMjym+li/+u07UVk3bTbt1hmzasCwOeA7wP/ljvmj5B1t29c4DFfRXaFaidg66p9G6R/tya7cvJd4DUFfc4/Bg6sHEv+cwBeSXb15IfAPh3627sSuAX4C2DT6s8gPd+bbEjArcDUTtSrTcfqOFu/vW/jLPdZd02sDVKcpWMpLdbKiLPqv6+q7W2PNcfZiDoNTKyVGWe579+/aVV/Y0XH2iDGWdsPoMWDPw+4DrgdmFU5cLIu0cvIxsxulbYfT9blPK3A8pvtYv9QbltL3dxVf3D5cdEvI5vQemD6I70JOJ1spZ4LKWh8NvD3wE1V2yYB2wIbVm2fksrersUyvwHcTLbKUd3Pj2yS61eB7wAb1/tPqqDP4aPAbVXbNkr/Pq9q+z7A2WMNvm55OM4GI85SXl0Va4MUZ+kYSou1MuIs5VFarDnORpQ1MLFWZpylPP2bFj53zG0rNM7aUvGCDv4q4Aay5TM/CawFXlv5A01/dJeRjRU9GXgU2LPgOmyQyrkUqNx/oToAvwtcVFB5+cDLj4v+VXr8Y27/3wPnkE0yfQZYUFAdTgD+vfJHBvwlcFv6j+1bwOvTvgmVz6jF8l4MfI90ZSxt246si3kW8GdV6V8ObNuBv79/Zv3k3Y3IhlF8L33f/052f5B8+sKuWnXy4TgbjDhLeXRdrA1KnKW6lxprnY6zynHlnnc81hxnI8oZiFgrO84qf0edjLWy4yzl63PH6EyctfUAWjjwLwO/yb1+EfA/wEG5bRsAryG7avEksFeB5Xe8i70q8K4kdZuTjQfeMQXFGuDjVe87OAXqrkWUDxyZ/uM7i+z+Ir8nuyLwFuC/ya4gFTZ5FtgZ+CPZFZEtgLcCjwM/IxuHvRjYsYS/wbPIJhLPBD5O9p/714AzUgBenP5z6tnJ246zwYmzVGbXxdogxFk6ztJirYw4S3mVFmuOs5r16vtYKzPOUt7+TfO5Y9vjrKMH1ORBb03WjX0R8Lq0bTtgZfoiPgkcDWyW9u1Ggd35lN/F3sy46DdW7ZtYYPkvSn9sNwJfAo7K7TuQbEz01kWVl/I9O/3HchOwBPhw2v5q4EfA20r4O9yC7OrIb8mums3O7TsZWNjpOhV8fI6zAYuzlHdXxVq/x1k6jtJirew4S/mWFmuOsxF16utYKzPOUn7+TfO5Y0fibCJdJiIelfQxsqsEsyTtAMwj+4LuI1up5lPACZJuj4h3FVV2ukneFsDnSUt2Svp1RNwo6UaybvbjgH+RtAw4ETgwIu4vqPydyMYD/7Ky2o6kiRHxp4hYk+p3JNmVhsrKORERhdxcMOX3O0mfJvuMN4qI1bkkewEPkV2lKtI/kl0l+iPw+4j4EUBE3Cap0p3fURHxB+DtkrYFVkfEqtzuDYGnJW0MrI0Ukb3EcTaQcQZdFmv9HmdQXqyVHWepDqXFmuNspH6PNf+mDeRv2mDGWadaoaM9yMYjb8D6buTpwPnAw8AlVWlfSNYF/WcFlt/xLnaqJkjS+XHRp+Y/Q+pcfSGb2Ho8WZfv7gV/73W7i8muEj1AwVeAR/s+yE2kJTeZNH0/7wMeA/68E3VqwzE6zgYwzirH3WBfx2JtEOIsHUdpsVZGnOW/26oyOhJrjrP630c/x1qZcZby9G+azx07GmdtP6AmD/pDZC3jy8i6s6em7dulADyf9t7VvuNd7Iwcl1zGuOi56TP/3+o/7hT470rPZ5CN2f0ZsEfBn/uE3PPXkG78CBxANh74kVb/g2uiDrPIrtoM1ajTacAB6fm+6bP/PwqeyNyph+NsMOOscry55x2PtUGKs3QcpcVaGXGW8ikt1hxnI+owMLFWZpylcvybNnKfzx2j/XHWtoMa4wfwV+kP4VqyG/etBD4G7El2M8DzyW6kd3gb6/BysomaXwXeBTwInAnMAU5Kr+8Ezi+grHzglTIummyJ2BvSH/gvWX8H928BS3PBvmUKhlaXba07ljr9J/MdsuEEk4F3kE1knd7mv7tvkU1ivRFYBxya23cRcH+l3sALyJbXfVk769Tm43Wc9Xmcpby6KtYGLc7ScZQaa52Ms1ReqbHmOHu23IGKtbLjLNXBv2nRvlhznNUov50HN8YP4ospwLYkG//772Tdih9OwXAF8BlgcoFldryLvSrwSlm2luzKzEbpsz005f9bYCj9Z1DznhMtlHcm2dWJjWrs+1T6T23DqvpNavPf27eB7+defwm4g+weH68FrqnUiaphCb38cJz1b5ylvLoq1gY1ztLxdDTWyoiz6u+tjFhznD1bxkDGWqfjLJXp3zSfO1ZedzzO2nZw4/gw3kd2RWDb9Hp2+mNcRnYH8l8AL/7/7Z1PqFVVFIe/9cxQQwKDwAzTyMhGaiGokBr+gQaF5kAIUimoxIE4iUSNcKCjEhsoCSWCYJREk0SUqEkIDZo2EQ3MCIWCIkvK1WBteder76V2997nnvP74AzuvYe79t7nfXe//X+A8WoPsdeYF72IG4dzNxE7/UwiegnOAo+mzyb+n1h9cVelZ7umX0Cix6x3Hnb2CgV4Pj37NX1lcZTRM0ZuOjSzDZc8a69n6fsa41qXPUt5KuZabc9SrKKuybMbYnbWtZKepe9Xnab/Hat6ljWTd1Eo3xFzYpcAF4jhxUnEQX4D27a156FXGWKnzrzowym/HxBnbIyk9/eSDmckzgA5ywAX7vbEfyb9kK7tFzB9/jQwudDf2SxgB7GF7Ir0N3Ye2FIifu1LnrXXs/T9jXCt656lMijiWk3PUvyirsmzm2J12rVSnqVYqtNc/zvW9Ox6y7AqZjbi7tfMbCOxW839wA53P5g57n6i4N8kTpOfTxxKuJcYWn6OmMO722/c0vNO44y4+7W+954gtq2F2M70beJH4Bwx5H99B5qBbFtrZouJ09uPEHOUjZgfPZ0Y1t6a7jtFbCu7wAe0VWxPGpYCh4ih+xPufiW9v504PHGVu18aZMxx0jITeJnoyZkLbHf3fSVi10KedcOz9P2NcK2LnkEd10p5lmJVdU2e3TItnXNNdVo36jR51kfJVuJttGZnEK3XY4XilR5ir71t7RJiQeVKYu73VuAy0Zuysue+hzOW+dJUvi+m19tSmS/M/Kxn0dcrk97bSfQUXZ8fPkKL5v+PURbybPTeVnqWvr+4a/LspvIo5lppz1KMmkcedNazFEeujeZbddrova2s0+RZT/wSQe6wkDalhzOnULxSU1mqr0FI8ZYnAdek1/OAdcB9jHMuxoDTsJSYH/0xcer9U5nj7Ull/yHwQt9ns4mh7zPA2hL5b8Ilz7Lnt7pnKW4x1+TZmOVSzLVSnqVY1V3romcpnly7uUxUp+XNb3XX5FmKCWCcjgAABA5JREFUXTLYbRbUI8Ap0tkUGeNcn6e7kWg1XwZezxiv+ra1PWlZRpwpsa7ic14G/Ezm82OIXYZ2pbJ/D/gdOAa80XPPg0Tvyff0bOva5kuedcOznnRkdU2ejVs22V0r7VmK1QjXuuRZiiPXbl0uqtPyl3F11+RZQ9a49WNmkz3NYS0QawYxtH3G3ddnjlVsDcJtpGUpcIBYbHrc3a/mjDdGGia5+58F4kwjDoHcQ/z47SDmJk8E9hHngPxN2jrY3c/nTlMTkGfd8CylI7tr8mxsSrlW0rMUrxGudcmzFEeu3QLVad2o0zrvWa1Wc5MuCg2xU2ENwn+kZwXwLTC19jPIkLed9MzxJha1fgFMSa83ED0pZ4gerPW0bJvkpl3yTJ7JsyLPpOTUzMa41mbPUv7kWsMu1Wntc20YPKteSE24KDTEnmIVW4Nwm+mZUrv8M+RpS49Ys9J7y4idluYTPSYX0o/uTGAz8GTtdLf9kmf1n8GA8yTPGniV9CzFa4xrbfQs5UuuNfBSnVb/GQw4T0PhWSOnStYg9xB7rW1ru4iZzSbOHplL9IiscPdzZvY+8BIxtL1LZV8eedYe5FlzKTFlTK6VQ641F9Vp7WFYPBupGbxJ5K7kfPQsjlPEosfTtR9+GzGzCcBFYu7xK8DXwDdJyKPE9IJ33f2gmd1TL6XdRJ61A3nWbHJ7lmLItQLItWajOq0dDJNnargVxt1/BN4BFpjZnNrpaQtmtsjMJrj7P+7+F3CJ2Dp3M3F45GngJ2Je9moAz3Agq2gG8iwP8kz0I9fyINdEL/IsD8PomRpudfgS+AH4pXZC2oCZHSYW6B4ws91pasFHxILS19x9A7F493PgE+BXM3vAzKxaokUJ5NkAkWdiHOTaAJFrYgzk2QAZVs+0xq0SJbetbTtmtpjoGTlCHIpowElgOnCvu29N950EHgcWuvulSskVBZFng0OeifGQa4NDromxkGeDY1g9U8NNtAIzWwIcBLYRu/1MJc7cmAasdvdT6b6H3P1itYQKMcTIMyHKINeEyM8weqaGm2gNZrYc2E/s+vOZmc0DHgNOAFd6FvkKIe4SeSZEGeSaEPkZNs/UcBOtwsyWAYeAt9z908rJEaKVyDMhyiDXhMjPMHmmrWNFq3D3r8zsVWKx6UTguLtfrZ0uIdqEPBOiDHJNiPwMk2cacROtxMxWAHuAZ939t9rpEaKNyDMhyiDXhMjPMHimhptoLWY2xd3/qJ0OIdqMPBOiDHJNiPw03TM13IQQQgghhBCi4egAbiGEEEIIIYRoOGq4CSGEEEIIIUTDUcNNCCGEEEIIIRqOGm5CCCGEEEII0XDUcBNCCCGEEEKIhqOGmxBCCCGEEEI0HDXchBBCCCGEEKLh/AtvxyVcc4p+bwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize = (15,4))\n",
    "\n",
    "lw = 2\n",
    "label_fs = 16\n",
    "\n",
    "whis = [0,100]\n",
    "showfliers = False\n",
    "\n",
    "subplot(1,4,1)\n",
    "title(r'$\\ell^2$ Error', fontsize = 16)\n",
    "bp1 = boxplot([L2_err[key] for key in L2_err.keys()], positions = np.arange(6), \\\n",
    "              patch_artist=True, whis=whis, showfliers=showfliers)\n",
    "for box, whisker,cap, median in zip(bp1['boxes'], bp1['whiskers'], bp1['caps'], bp1['medians']):\n",
    "    box.set(color='b', linewidth=lw)\n",
    "    box.set(facecolor = 'w' )\n",
    "    median.set(color='k', linestyle = '--', linewidth=1.5)\n",
    "yticks(fontsize = 12)\n",
    "xticks(np.arange(6), ['ARD','ARDr','ARDvi','M-ST','L-ST','MAP-ST'], fontsize = 12, rotation=45)\n",
    "    \n",
    "subplot(1,4,2)\n",
    "title(r'$\\ell^1$ Error', fontsize = 16)\n",
    "bp2 = boxplot([L1_err[key] for key in L1_err.keys()], positions = np.arange(6), patch_artist=True, whis=whis, showfliers=showfliers)\n",
    "for box, whisker, cap, median in zip(bp2['boxes'], bp2['whiskers'], bp2['caps'], bp2['medians']):\n",
    "    box.set(color='r', linewidth=lw)\n",
    "    box.set(facecolor = 'w' )\n",
    "    median.set(color='k', linestyle = '--', linewidth=1.5)\n",
    "yticks(fontsize = 12)\n",
    "xticks(np.arange(6), ['ARD','ARDr','ARDvi','M-ST','L-ST','MAP-ST'], fontsize = 12, rotation=45)\n",
    "\n",
    "subplot(1,4,3)\n",
    "title(r'Added', fontsize = 16)\n",
    "bp3 = boxplot([Added[key] for key in Added.keys()], positions = np.arange(6), patch_artist=True, whis=whis, showfliers=showfliers)\n",
    "for box, whisker, cap, median in zip(bp3['boxes'], bp3['whiskers'], bp3['caps'], bp3['medians']):\n",
    "    box.set(color='g', linewidth=lw)\n",
    "    box.set(facecolor = 'w' )\n",
    "    median.set(color='k', linestyle = '--', linewidth=1.5)\n",
    "yticks(fontsize = 12)\n",
    "xticks(np.arange(6), ['ARD','ARDr','ARDvi','M-ST','L-ST','MAP-ST'], fontsize = 12, rotation=45)\n",
    "yscale('symlog', linthreshy = 1)\n",
    "\n",
    "subplot(1,4,4)\n",
    "title(r'Missed', fontsize = 16)\n",
    "bp4 = boxplot([Missed[key] for key in Missed.keys()], positions = np.arange(6), patch_artist=True, whis=whis, showfliers=showfliers)\n",
    "for box, whisker, cap, median in zip(bp4['boxes'], bp4['whiskers'], bp4['caps'], bp4['medians']):\n",
    "    box.set(color='c', linewidth=lw)\n",
    "    box.set(facecolor = 'w' )\n",
    "    median.set(color='k', linestyle = '--', linewidth=1.5)\n",
    "yticks(fontsize = 12)\n",
    "xticks(np.arange(6), ['ARD','ARDr','ARDvi','M-ST','L-ST','MAP-ST'], fontsize = 12, rotation=45)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}