{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p56TXrchLlUN"
      },
      "source": [
        "# Part 1 - Bootstrapping"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "f0iIfAGvLlUP"
      },
      "source": [
        "In a medical study, doctors aimed to examine whether aspirin can reduce the incidence of heart attacks. To this aim, they recruited a large amount of subjects and divided them into an aspirin group and placebo group. After receiving the aspirin or the placebo pill for a period of time, doctors observed that people in the aspirin group ended up only having 104 heart attacks incidents, and people in the placebo group ended up having 189 incidents.\n",
        "\n",
        "Given the data, we can compute the percentages of the heart attack incident for both groups, and compute a ratio. The fact that the ratio is smaller than 1 suggests that aspirin therapy is effective in preventing heart attacks. But how sure are we of this estimate? Does the 95%\n",
        "confidence interval include 1?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "pWNv4OqsLlUQ",
        "outputId": "4525d48b-14c5-4ae7-a3ed-f40c43bb5528"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount(\"/content/drive\", force_remount=True).\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "from sklearn.linear_model import LinearRegression\n",
        "from random import choices\n",
        "from scipy.stats import norm, pearsonr\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "from scipy.optimize import minimize\n",
        "import sys\n",
        "import math\n",
        "from google.colab import drive\n",
        "from scipy.io import loadmat\n",
        "\n",
        "drive.mount('/content/drive')\n",
        "sys.path.append('/content/drive/My Drive/Colab Notebooks/')\n",
        "\n",
        "np.random.seed(0)\n",
        "cMAP1  = np.array([[50,205,50],[255,165,0]])/255"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BwhrXv0NLlUQ",
        "outputId": "2c11d49a-7d0c-4b55-e654-af3e863d3239"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The ratio computed from empirical data is  0.5501149812103875\n"
          ]
        }
      ],
      "source": [
        "#aspirin group\n",
        "aspirin_heart   = 104\n",
        "aspirin_total   = 11037\n",
        "aspirin_data    = np.concatenate((np.ones(aspirin_heart), np.zeros(aspirin_total-aspirin_heart)))\n",
        "\n",
        "#placebo group\n",
        "placebo_heart   = 189\n",
        "placebo_total   = 11034\n",
        "placebo_data    = np.concatenate((np.ones(placebo_heart), np.zeros(placebo_total-placebo_heart)))\n",
        "\n",
        "#Calculate statistic for original sample\n",
        "ratio_empirical = (aspirin_heart/aspirin_total)/(placebo_heart/placebo_total)\n",
        "print('The ratio computed from empirical data is ', str(ratio_empirical))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "ooPzcJRGLlUR"
      },
      "outputs": [],
      "source": [
        "#Let's say we bootstrap (sample with replacement) 10,000 times\n",
        "n_boot          = int(1e4)\n",
        "#ratio_boot is used to store the ratio for each bootstrapped dataset\n",
        "ratio_boot      = np.zeros(n_boot)\n",
        "\n",
        "\n",
        "# ----------------- YOUR CODE HERE ----------------- #\n",
        "\n",
        "\n",
        "# ----------------- YOUR CODE HERE ----------------- #\n",
        "\n",
        "\n",
        "#Find 95% confidence interval\n",
        "ratio_boot_sorted = np.sort(ratio_boot) # Sort ratio_boots from lowest to highest\n",
        "lower_bound_ind   = int(np.ceil(0.025 * n_boot)) # Find index of 2.5% value\n",
        "upper_bound_ind   = int(np.floor(0.975 * n_boot)) # Find index of 97.5% value\n",
        "\n",
        "#Use lower/higher index to find value corresponding to 2.5%/97.5% position\n",
        "lower_bound       = ratio_boot_sorted[lower_bound_ind]\n",
        "upper_bound       = ratio_boot_sorted[upper_bound_ind]\n",
        "#----------------------------------------------------------------------"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D-72p3tELlUS"
      },
      "source": [
        "Let's visualize the bootstrapped ratios and the 95% confidence interval."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 467
        },
        "id": "vsonfT0RLlUS",
        "outputId": "3dc1cefc-ca38-4ee1-fef5-8716b1110a98"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x78ace234fb50>"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfQklEQVR4nO3dd1gU1/s28Ht36YuAKNWgqFgAUbGvGks0YonR2I1RMJZEQWOMJfausUajqCkqaiyJGv0maiwhwViwoVgRG4pRig0QkAV25/3Dl/m5AgoIO+xyf65rrmt25uzsM+OO+3DOmXNkgiAIICIiIjJScqkDICIiIipJTHaIiIjIqDHZISIiIqPGZIeIiIiMGpMdIiIiMmpMdoiIiMioMdkhIiIio2YidQClgVarxYMHD1CuXDnIZDKpwyEiIqICEAQBz549g6urK+Ty/OtvmOwAePDgAdzc3KQOg4iIiIrg3r17eOedd/Ldz2QHQLly5QC8uFg2NjYSR0MFkZaWBldXVwAvklWlUilxRPRaaWnA///3woMHAP+9iKgYpKSkwM3NTfwdzw+THUBsurKxsWGyYyAUCoW4bmNjw2SntHvp3ws2Nkx2iKhYvakLCjsoExGR3mRkZKB3797o3bs3MjIypA6HyggmO0REpDcajQY7d+7Ezp07odFopA6Hygg2Y5FBMjExgb+/v7hORESUH/5KkEEyNzdHSEiI1GEQFZpGo0FWVpbUYUhGrVajSpUq4vrL/e+IXmVqalos3xEmO0REeiAIAuLj45GUlCR1KJLSarVYu3YtACAuLu61Y6MQAYCdnR2cnZ3fahw8JjtkkARBQHp6OgDAysqKg0FSqZeT6Dg6Opbp76xGo8Hz588BAO7u7qzZoXzl/D+fmJgIAHBxcSnysZjskEFKT0+HtbU1ACA1NZWPnlOpptFoxESnQoUKUocjqZc7JVtYWDDZodeytLQEACQmJsLR0bHI3xfWHxIRlbCcPjpWVlYSR0JkeHLum7fp68Zkh4hIT8pq0xXR2yiO+4bJDhERERk1JjtERERk1JjsEBERSSwkJAR2dnZ6+7ywsDDIZLIyMxQCkx0iIspXQEAAZDKZuFSoUAEdO3bExYsXi/Vz7ty5A5lMhsjIyAK/Z+bMmahfv36xxmGM2rRpgzFjxuhsa968OeLi4mBraytNUHrGZIcMkkKhQK9evdCrVy8+ukpUwjp27Ii4uDjExcUhNDQUJiYm+OCDD4p0LJlMhvLly6N8+fJ667BtrCNWv815mZmZvfVAfYaEyQ4ZJAsLC+zYsQM7duyAhYWF1OEQFZ4gAGlp0iyCUKhQzc3N4ezsDGdnZ9SvXx9ff/017t27h4cPH4plLl26hPfeew+WlpaoUKEChg8fjtTUVHG/VqvF7NmzUblyZXh5eaFnz544dOiQuL9q1aoAAF9fX8hkMrRp0wbAi+aWJk2aQKlUws7ODi1atMDdu3cREhKCWbNm4cKFC2KtU84UMjKZDGvWrMGHH34IpVKJefPmQaPRYMiQIahatSosLS1Rq1YtrFixQuc8AwIC0L17d8yaNQsODg6wsbHB559/jszMTLFMmzZtEBQUhKCgINja2qJixYqYNm0ahJeuqVqtxrhx41CpUiUolUo0bdoUYWFhOp8VEhKCypUrw8rKCh999BEeP3782n+DnJqvX375Ba1bt4aFhQW2bNmCx48fo3///qhUqRKsrKzg4+ODbdu26ZzTkSNHsGLFCvE63blzJ89mrF27dsHb2xvm5uZwd3fH0qVLdWJYvXo1atSoAQsLCzg5OaFXr16vjblUEUhITk4WAAjJyclSh0JkFObenauzLIyaLggvfmKFhVHThbl350odol49f/5cuHr1qvD8+fP/25iaKl4TvS+pqQWO3d/fX+jWrZv4+tmzZ8Jnn30meHh4CBqN5v+fSqrg4uIi9OjRQ7h06ZIQGhoqVK1aVfD39xfft2zZMsHGxkbYtm2bcO3aNWHChAmCqampcP36dUEQBOH06dMCAOGvv/4S4uLihMePHwtZWVmCra2tMG7cOOHmzZvC1atXhZCQEOHu3btCenq68NVXXwne3t5CXFycEBcXJ6SnpwuCIAgABEdHR2H9+vXCrVu3hLt37wqZmZnC9OnThTNnzgi3b98Wfv75Z8HKykr45ZdfdM7V2tpa6Nu3r3D58mVh7969goODgzB58mSxTOvWrQVra2vhiy++EK5duyYe54cffhDLDB06VGjevLnw77//Cjdv3hQWL14smJubi+d68uRJQS6XCwsXLhSio6OFFStWCHZ2doKtrW2+/w4xMTECAMHd3V3YtWuXcPv2beHBgwfCf//9JyxevFg4f/68cOvWLeG7774TFAqFcOrUKUEQBCEpKUlQqVTCsGHDxOuUnZ0t/PPPPwIA4enTp4IgCMLZs2cFuVwuzJ49W4iOjhY2bNggWFpaChs2bBAEQRDOnDkjKBQKYevWrcKdO3eEc+fOCStWrCjw9+ht5Hn//H8F/f1msiMw2SEqbkx2dBl6sqNQKASlUikolUoBgODi4iJERESIZX744QehfPnyQupLx923b58gl8uF+Ph4QRAEwdXVVZg3b57OsRs3biyMHDlSEIT/+zE/f/68uP/x48cCACEsLCzP2GbMmCHUq1cv13YAwpgxY954boGBgULPnj11ztXe3l5IS0sTt61Zs0awtrYWE7vWrVsLnp6eglarFctMnDhR8PT0FARBEO7evSsoFArh/v37Op/Vrl07YdKkSYIgCEL//v2Fzp076+zv27dvgZKd5cuXv/G8unTpInz11Vfi69atWwtffPGFTplXk52PP/5YeP/993XKjB8/XvDy8hIEQRB27dol2NjYCCkpKW/8/OJWHMkOp4sgg5SWlsbpIsiwWVkBLzXz6P2zC6Ft27ZYs2YNAODp06dYvXo1OnXqhNOnT6NKlSqIiopCvXr1dO7DFi1aQKvVIjo6GpaWlnjw4AFatGgBjUaD8+fPAwBUKhUuXbqU7+fa29sjICAAfn5+eP/999G+fXv06dOnQHMkNWrUKNe24OBgrF+/HrGxsXj+/DkyMzNzdXCuV6+ezkjXKpUKqampuHfvnjhbe7NmzXT6uqhUKixduhQajQaXLl2CRqNBzZo1dY6rVqvFqUKioqLw0Ucf6exXqVQ4cOBAoc9Lo9Fg/vz5+PXXX3H//n1kZmZCrVYXerTuqKgodOvWTWdbixYtsHz5cmg0Grz//vuoUqUKqlWrho4dO6Jjx4746KOPDGZUcCY7RERSkMkAA0nSlUolPDw8xNc//fQTbG1t8eOPP2Lu3Lkl+tkbNmzA6NGjceDAAfzyyy+YOnUqDh8+jGbNmr0x5pdt374d48aNw9KlS6FSqVCuXDksXrwYp06dKtZ4U1NToVAoEBERkevhiZw/0N7Gq+e1ePFirFixAsuXL4ePjw+USiXGjBmj08+oOJQrVw7nzp1DWFgYDh06hOnTp2PmzJk4c+aMXh+ZLyp2UCYiokKRyWSQy+Xi7OWenp64cOEC0tLSxDLHjx+HXC5HrVq1YGNjA1dXV3FbvXr1UK9ePYSHh8PLywvAi6eDAN2JQnP4+vpi0qRJOHHiBOrUqYOtW7eK78mrfF6OHz+O5s2bY+TIkfD19YWHhwdu3bqVq9yFCxfE8wKAkydPwtraGm5ubuK2VxOkkydPokaNGlAoFPD19YVGo0FiYiI8PDx0FmdnZ/F65XWMojh+/Di6deuGTz75BPXq1UO1atVw/fp1nTIFuU6enp44fvx4rmPXrFlTTNpMTEzQvn17LFq0CBcvXsSdO3fw999/FylufWOyQ0REr6VWqxEfH4/4+HhERUVh1KhRSE1NRdeuXQEAAwYMgIWFBfz9/XH58mX8888/GDVqFAYOHAgnJycAwPjx47Fw4UL8+uuvuH37NqZNm4bIyEh88cUXAABHR0dYWlriwIEDSEhIQHJyMmJiYjBp0iSEh4fj7t27OHToEG7cuAFPT08AgLu7O2JiYhAZGYlHjx5BrVbnew41atTA2bNncfDgQVy/fh3Tpk3DmTNncpXLzMzEkCFDcPXqVezfvx8zZsxAUFAQ5PL/+7mMjY3F2LFjER0djW3btmHlypXiedSsWRMDBgzAoEGD8NtvvyEmJganT5/GggULsG/fPgAQa6qWLFmCGzduYNWqVQVqwsrvvA4fPowTJ04gKioKn332GRISEnTKuLu749SpU7hz5w4ePXoErVab6zhfffUVQkNDMWfOHFy/fh0bN27EqlWrMG7cOADA3r178d133yEyMhJ3797Fpk2boNVqUatWrSLFrXcl1aHIkLCDsuFJTU0VAAgAdDpFUunADsq6XtfBsrTz9/cX7zUAQrly5YTGjRsLO3fu1Cl38eJFoW3btoKFhYVgb28vDBs2THj27Jm4X6PRCDNnzhQqVaokmJqaCvXq1RP+/PNPnWP8+OOPgpubmyCXy4XWrVsL8fHxQvfu3QUXFxfBzMxMqFKlijB9+nSxs3BGRobQs2dPwc7OTgAgPjkEQNi9e7fOsTMyMoSAgADB1tZWsLOzE0aMGCF8/fXXOh2cc548mz59ulChQgXB2tpaGDZsmJCRkSGWad26tTBy5Ejh888/F2xsbITy5csLkydP1umwnPPkl7u7u2Bqaiq4uLgIH330kXDx4kWxzLp164R33nlHsLS0FLp27SosWbKkQB2UX+7ALQgvOnF369ZNsLa2FhwdHYWpU6cKgwYN0nmCLjo6WmjWrJlgaWkpABBiYmJydVAWBEHYuXOn4OXlJZiamgqVK1cWFi9eLO47evSo0Lp1a6F8+fKCpaWlULduXZ0n2UpScXRQlglCIQdcMEIpKSmwtbVFcnIybGxspA6HCoAdlEu3ebHzdF6bpmdigudsAMCiqOnIsjLDlMpTpAhNEhkZGYiJiUHVqlXL/LhQWq0W9+7dAwC4ubnp1JhILSAgAElJSdizZ0++Zdq0aYP69etj+fLleourrHvd/VPQ3+/S8y0jIiKjJwgCHj58iIcPH4J/a5O+8GksMkgKhQKdO3cW14mIiPIjec3O/fv38cknn6BChQqwtLSEj48Pzp49K+4XBAHTp0+Hi4sLLC0t0b59e9y4cUPnGE+ePMGAAQNgY2MDOzs7DBkyRGeYcjI+FhYW2LdvH/bt21fmmwWIqHiEhIS8tgkLeDF9BZuwDI+kyc7Tp0/RokULmJqa4s8//8TVq1exdOlSlC9fXiyzaNEifPfdd1i7di1OnToFpVIJPz8/ZGRkiGUGDBiAK1eu4PDhw9i7dy/+/fdfDB8+XIpTIiIiolJG0mashQsXws3NDRs2bBC35UwGB7yo1Vm+fDmmTp0qjuy4adMmODk5Yc+ePejXrx+ioqJw4MABnDlzRhxZcuXKlejcuTOWLFkCV1dX/Z4UERERlSqS1uz8/vvvaNSoEXr37g1HR0f4+vrixx9/FPfHxMQgPj4e7du3F7fZ2tqiadOmCA8PBwCEh4fDzs5OZwjt9u3bQy6X5zsyplqtRkpKis5ChiUtLQ1KpRJKpVJnIDMiIqJXSZrs3L59G2vWrEGNGjVw8OBBjBgxAqNHj8bGjRsBAPHx8QAgDkqVw8nJSdwXHx8PR0dHnf0mJiawt7cXy7xqwYIFsLW1FZeXR8Ykw5Geno709HSpwyAiolJO0mRHq9WiQYMGmD9/Pnx9fTF8+HAMGzYMa9euLdHPnTRpEpKTk8UlZ8wHIiIiMj6SJjsuLi7ivCg5PD09ERsbCwDiPCKvDn2dkJAg7nN2dkZiYqLO/uzsbDx58kQs8ypzc3PY2NjoLEREVLoEBASge/fubywnk8ne+BRVYbi7u+v1ias2bdpgzJgxevu8skjSZKdFixaIjo7W2Xb9+nVUqVIFwIvOys7OzggNDRX3p6Sk4NSpU1CpVAAAlUqFpKQkREREiGX+/vtvaLVaNG3aVA9nQURkvAICAiCTyXItHTt2LPHPXrFiBUJCQt5YLi4uDp06dSrxeN5WWFgYZDIZkpKSdLb/9ttvmDNnjjRBlRGSPo315Zdfonnz5pg/fz769OmD06dP44cffsAPP/wA4EW2PmbMGMydOxc1atRA1apVMW3aNLi6uorZvqenJzp27Cg2f2VlZSEoKAj9+vXjk1hEpdirU0rkpSxNKVGadezYUeepWeBFDXlJs7W1fe3+zMxMmJmZ5VuLry85cRSVvb19MUZDeZG0Zqdx48bYvXs3tm3bhjp16mDOnDlYvnw5BgwYIJaZMGECRo0aheHDh6Nx48ZITU3FgQMHdAaS27JlC2rXro127dqhc+fOaNmypZgwERHR2zE3N4ezs7PO8vJ4aDKZDN9//z0++OADWFlZwdPTE+Hh4bh58ybatGkDpVKJ5s2b49atW+J7fvjhBzRs2BDff/893NzcYGVlhT59+iA5OVks82ozVps2bRAUFIQxY8agYsWK8PPzEz//5Was//77D/3794e9vT2USiUaNWokPp1769YtdOvWDU5OTrC2tkbjxo3x119/Fep65MQ1b948uLq6ijN/b968GY0aNUK5cuXg7OyMjz/+WOxmcefOHbRt2xYAUL58echkMgQEBIjn9XIz1tOnTzFo0CCUL18eVlZW6NSpU67BdKlwJJ8u4oMPPsAHH3yQ736ZTIbZs2dj9uzZ+Zaxt7fH1q1bSyI8KqXkcjlat24trpNxKgu1P68bOkGhUOj8Yfe6snK5HJaWlm8sW1KT5s6ZMwfLli3DsmXLMHHiRHz88ceoVq0aJk2ahMqVK+PTTz9FUFAQ9u3bh3LlysHMzAw3b97Er7/+ij/++AMpKSkYMmQIRo4ciS1btuT7ORs3bsSIESNw/PjxPPenpqaidevWqFSpEn7//Xc4Ozvj3Llz0Gq14v7OnTtj3rx5MDc3x6ZNm9C1a1dER0ejcuXKBT7f0NBQ2NjY4PDhw+K2rKwszJkzB7Vq1UJiYiLGjh2LgIAA7N+/H25ubti1axd69uyJ6Oho2NjY6Px7vSwgIAA3btzA77//DhsbG0ycOBGdO3fG1atXYWpqWuAY6f9InuwQFYWlpSXCwsKkDoPorVlbW+e7r3Pnzti3b5/42tHRMd/hFlq3bq1zT7i7u+PRo0e5yhVl8s29e/fminPy5MmYPHmy+Hrw4MHo06cPAGDixIlQqVSYNm2aWPvyxRdfYPDgwZDL5ahVqxYqVqyIjIwMbNq0CZUqVQLwYkDYLl26YOnSpfk2TdWoUQOLFi3KN9atW7fi4cOHOHPmjNg85OHhIe6vV68e6tWrJ76eM2cOdu/ejd9//x1BQUEFviZKpRI//fSTTvPVp59+Kq5Xq1YN3333ndgiYW1tLcbj6OgIOzu7PI+bk+QcP34czZs3B/Ci9cLNzQ179uxB7969Cxwj/R8mO0RE9Fpt27bFmjVrdLa92s+kbt264nrO2Gg+Pj462zIyMpCSkiI+AVu5cmUx0QFePHCi1WoRHR2db7LTsGHD18YaGRkJX1/ffPvBpKamYubMmdi3bx/i4uKQnZ2N58+fi08BF5SPj0+ufjoRERGYOXMmLly4gKdPn4q1SbGxsbmePM5PVFQUTExMdB6wqVChAmrVqoWoqKhCxUj/h8kOEZGEXjdpsUKh0Hn96jAbL3u1OffOnTtvFdfLlEqlTu1IXl5uXpHJZPluy0kA3iaW18mvaSjHuHHjcPjwYSxZsgQeHh6wtLREr169kJmZ+VZxpKWlwc/PD35+ftiyZQscHBwQGxsLPz+/Qh+bih+THTJIaWlpcHd3B/DiP/WS6odAVNIK890tqbL6pNFoEBkZifj4eMTGxuLBgwfik7MnT54Um7mKqm7duvjpp5/w5MmTPGt3jh8/joCAAHz00UcAXiSbxZEYXrt2DY8fP8Y333wjjsp/9uxZnTI5NUEajSbf43h6eiI7OxunTp0Sm7EeP36M6OjoAtcOUW7s2UkG69GjR3n2SSCi4qVWqxEfH6+zvM29l52dDa1WCwsLC/j7++PChQs4evQoRo8ejT59+rzVo+T9+/eHs7MzunfvjuPHj+P27dvYtWuXOJ9ijRo18NtvvyEyMhIXLlzAxx9//Na1TcCLJjkzMzOsXLkSt2/fxu+//55r7JwqVapAJpNh7969ePjwYZ61ejVq1EC3bt0wbNgwHDt2DBcuXMAnn3yCSpUqiRNiU+Ex2SEiotc6cOAAXFxcdJaWLVsW6VhyuRze3t5wcHCAh4cHevTogc6dO6NDhw6oW7cuVq9e/VaxmpmZ4dChQ3B0dETnzp3h4+ODb775RmwSXLZsGcqXL4/mzZuja9eu8PPzQ4MGDd7qMwHAwcEBISEh2LFjB7y8vPDNN99gyZIlOmUqVaqEWbNm4euvv4aTk1O+HaI3bNiAhg0b4oMPPoBKpYIgCNi/fz+fxHoLMqEoXfONTEpKCmxtbZGcnMypIwxEWlqa+HRIampqqa2yL6tefWTcND0TEzxfDB+xKGo6sqyKPgDbqwzh0fOMjAzExMSgatWqOo+Sl2UzZ87Enj17EBkZKXUoVMq97v4p6O83a3aIiIjIqLGDMhER6Y1Wq0V8fDyePXsmdShUhrBmh4iI9EYQBDx48AD9+/fXmcCZqCSxZocMklwuR6NGjcR1IiKi/DDZIYNkaWmJM2fOSB0GEREZAP5JTEREREaNyQ4REREZNSY7ZJDS09Ph7u4Od3f3fGeBJiIiAthnhwyUIAi4e/euuE5ERJQfJjtEJHp15OO8GMKIxUREL2MzFhERlUmtWrXC1q1bxdcymQx79uyRLqACCggIQPfu3aUOI19t2rTBmDFjxNfNmjXDrl27pAsITHaIiOg1nj17hjFjxqBKlSqwtLRE8+bNcw37EBAQAJlMprN07NhR3K9WqzFw4EDY2NjA09MTp06d0nn/4sWLMWrUqALFk5KSgilTpqB27dqwsLCAs7Mz2rdvj99++01s0n71xzYvv//+OxISEtCvX78CfS4V3dSpU/H1118Xy+zyRcVkh4iI8jV06FAcPnwYmzdvxqVLl9ChQwe0b98e9+/f1ynXsWNHxMXFicu2bdvEfT/88AMiIiIQHh6OoUOHYtq0aWJiEhMTgx9//BHz5r25CTUpKQnNmzfHpk2bMGnSJJw7dw7//vsv+vbtiwkTJiA5ObnA5/Xdd99h8ODBpXpQ0szMTKlDKBadOnXCs2fP8Oeff0oWQ+n9VyYiKgPS0tIKvWRnZ4vvz87ORlpaGp4/f16g4xbG8+fPsWvXLixatAitWrWCh4cHZs6cCQ8PD6xZs0anrLm5OZydncWlfPny4r6oqCh8+OGH8Pb2xsiRI/H06VMkJSUBAEaMGIGFCxe+dsbqHJMnT8adO3dw6tQp+Pv7w8vLCzVr1sSwYcMQGRkJa2vrAp3Xw4cP8ffff6Nr166vLXfp0iW89957sLS0RIUKFTB8+HCkpqYCAC5fvgy5XI6HDx8CAJ48eQK5XK5TUzR37ly0bNlSfH358mV06tQJ1tbWcHJywsCBA/Ho0SNxf5s2bRAUFIQxY8agYsWK8PPze218s2bNgoODA2xsbPD555/rJEdqtRqjR4+Go6MjLCws0LJlS50auZCQENjZ2ekcb8+ePZDJZOLrmTNnon79+ti8eTPc3d1ha2uLfv366cxrlpaWhkGDBsHa2houLi5YunRprjgVCgU6d+6M7du3v/Z8ShKTHTJIMpkMXl5e8PLy0rk5iQyNtbV1oZfdu3eL79+9ezesra3RqVMnneO6u7vn+d7CyM7OhkajgYWFhc52S0tLHDt2TGdbWFgYHB0dUatWLYwYMQKPHz8W99WrVw/Hjh3D8+fPcejQITg4OMDZ2Rlbt26FhYUFPvroozfGotVqsX37dgwYMACurq659ltbW8PEpGDP3Bw7dgxWVlbw9PTMt0xaWhr8/PxQvnx5nDlzBjt27MBff/2FoKAgAIC3tzcqVKiAI0eOAACOHj2q8xoAjhw5gjZt2gB4USv13nvvwdfXF2fPnsWBAweQkJCAPn366Hzuxo0bYWZmhuPHj2Pt2rX5xhcaGoqoqCiEhYVh27Zt+O233zBr1ixx/4QJE7Br1y5s3LgR586dg4eHB/z8/PDkyZMCXaMct27dwp49e7B3717s3bsXR44cwTfffCPuHz9+PI4cOYL//e9/OHToEMLCwnDu3Llcx2nSpAmOHj1aqM8uTkx2yCBZWVnhypUruHLlCqysrKQOh8golStXDiqVCnPmzMGDBw+g0Wjw888/Izw8HHFxcWK5jh07YtOmTQgNDcXChQtx5MgRdOrUCRqNBgDw6aefol69evDy8sKCBQvw22+/oVKlSpg5cyZWrlyJqVOnij/GrzaP5Xj06BGePn2K2rVrv/V53b17F05OTq9twtq6dSsyMjKwadMm1KlTB++99x5WrVqFzZs3IyEhATKZDK1atUJYWBiAF8ne4MGDoVarce3aNWRlZeHEiRNo3bo1AGDVqlXw9fXF/PnzUbt2bfj6+mL9+vX4559/cP36dfFza9SogUWLFqFWrVqoVatWvvGZmZlh/fr18Pb2RpcuXTB79mx899130Gq1SEtLw5o1a7B48WJ06tQJXl5e+PHHH2FpaYl169YV6lpptVqEhISgTp06ePfddzFw4ECEhoYCAFJTU7Fu3TosWbIE7dq1g4+PDzZu3KhT85jD1dUV9+7dk6zfDh89JyKSUE6zSGGYm5uL6x999BFSU1Nz/XDfuXPnbUMDAGzevBmffvopKlWqBIVCgQYNGuSasfzlphsfHx/UrVsX1atXR1hYGNq1awdTU1MEBwfrHHfw4MEYPXo0zp8/jz179uDChQtYtGgRRo8eneeTO8U5ntbz589z1Va9KioqCvXq1YNSqRS3tWjRAlqtFtHR0XByckLr1q3xww8/AHhRizN//nxcv34dYWFhePLkCbKystCiRQsAwIULF/DPP//kWbt269Yt1KxZEwDQsGHDAp1DvXr1dP7QU6lUSE1Nxb1795CcnKzz2QBgamqKJk2aICoqqkDHz+Hu7o5y5cqJr11cXJCYmCjGnZmZiaZNm4r77e3t80zSLC0todVqoVarYWlpWagYigOTHaIyoiBj6OjzOPTCyz+mRWFiYpJn883bHjdH9erVceTIEaSlpSElJQUuLi7o27cvqlWrlu97qlWrhooVK+LmzZto165drv3//PMPrly5gp9++gnjx49H586doVQq0adPH6xatSrPYzo4OMDOzg7Xrl1763OqWLEinj59+tbHyXnq68aNG7h69SpatmyJa9euISwsDE+fPkWjRo3EhCQ1NRVdu3bFwoULcx3HxcVFXC+uf7c3kcvluRLIrKysXOVMTU11XstksiLVzjx58gRKpVKSRAdgMxYZqPT0dHh7e8Pb25vTRRDpgVKphIuLC54+fYqDBw+iW7du+Zb977//8PjxY50f8RxpaWkYOnQoxo8fDwDQaDTij2xWVpbY9PWqnM6/W7ZswYMHD3LtT01NzbP5JC++vr6Ij49/bcLj6emJCxcu6HTqPn78OORyuVhz4ePjg/Lly2Pu3LmoX78+rK2t0aZNGxw5cgRhYWFifx0AaNCgAa5cuQJ3d3d4eHjoLEVJcC5cuKDTKf3kyZOwtraGm5sbqlevLvb7yZGVlYUzZ87Ay8sLwIvk8dmzZzrnFxkZWagYqlevDlNTU52hBJ4+farTLJfj8uXL8PX1LdTxixOTHTJIgiDg6tWruHr1KqeLICpBBw8exIEDBxATE4PDhw+jbdu2qF27NgYPHgzgRZIxfvx4nDx5Enfu3EFoaCi6desm9sF51dy5c6FSqVC1alUAL5qGfvvtN1y8eBGrVq3SaXp51bx58+Dm5oamTZti06ZNuHr1Km7cuIH169fD19e3wE2Cvr6+qFixok4y8KoBAwbAwsIC/v7+uHz5Mv755x+MGjUKAwcOhJOTEwCI/Xa2bNkiJjZ169aFWq1GaGio2F8HAAIDA/HkyRP0798fZ86cwa1bt3Dw4EEMHjw43wTvdTIzMzFkyBBcvXoV+/fvx4wZMxAUFAS5XA6lUokRI0Zg/PjxOHDgAK5evYphw4YhPT0dQ4YMAQA0bdoUVlZWmDx5Mm7duoWtW7ciJCSkUDFYW1tjyJAhGD9+PP7++29cvnwZAQEBefaFOnr0KDp06FDo8ywuTHaIiChfycnJCAwMRO3atTFo0CC0bNkSBw8eFJs3FAoFLl68iA8//BA1a9bEkCFD0LBhQxw9elSnbxHw4q/7nTt3YsmSJahVqxbkcjl69eqFLl264N1338XFixexYsWKfGOxt7fHyZMn8cknn2Du3Lnw9fXFu+++i23btmHx4sWwtbUt0DkpFAoMHjwYW7ZsybeMlZUVDh48iCdPnqBx48bo1asX2rVrl6uZrXXr1tBoNGKyI5fL0apVK8hkMp3EzdXVFcePH4dGo0GHDh3g4+ODMWPGwM7Orkhj/bRr1w41atRAq1at0LdvX3z44YeYOXOmuP+bb75Bz549MXDgQDRo0AA3b97EwYMHxSEB7O3t8fPPP2P//v3w8fHBtm3bdN5fUIsXL8a7776Lrl27on379mjZsmWufkf379/HiRMnxARZCjKBfxYjJSUFtra2SE5OLtBYDyS9tLQ0saNfamqq3tq5DZmUfW1M0zMxwXM2AGBR1HRkWZkV27ENYa6ujIwMxMTEoGrVqm/sGEv6ER8fD29vb5w7dw5VqlSROhyjNnHiRDx9+lTszF1Yr7t/Cvr7zZodIiIqc5ydnbFu3TrExsZKHYrRc3R0xJw5cySNgU9jERGR3mi1WnHU4IoVK0o6XUNpnkzTmHz11VdSh8Bkh4iI9EcQBLE2pUKFChJHQ2UFkx0ySDKZTGxn53QRRET0Okx2yCBZWVkV2wixRERk3NhBmYiIiIwakx0iIiIyakx2yCA9f/4cjRs3RuPGjXWGTCciInoVkx0ySFqtFmfPnsXZs2eLNCkdEdGr1q1bJ+mUBsVh5syZqF+/vtRh5CsgIEDnkf9+/fph6dKlJf657KBMRCQhfY9sXdgRp589e4Zp06Zh9+7dSExMhK+vL1asWIHGjRuLZQICArBx40ad9/n5+eHAgQP5HvfPP//ERx99hNTUVAwePBjLli0T9925cwcdOnTA2bNnCzSq/a5du7By5UqcP38eGo0G1apVQ69evRAUFAR7e3uEhIRgzJgxSEpKyvcYGRkZmDZtGnbs2PHGz6PiM3XqVLRq1QpDhw4t8HQfRcGaHSIiytfQoUNx+PBhbN68GZcuXUKHDh3Qvn173L9/X6dcx44dERcXJy7btm3L95hJSUmYN28eFi1ahEOHDuHnn3/G3r17xf0jR47EN998U6BEZ8qUKejbty8aN26MP//8E5cvX8bSpUtx4cIFbN68ucDnuXPnTtjY2Lx2ItLSJDMzU+oQikWdOnVQvXp1/PzzzyX6OUx2iMjozYud98aFcnv+/Dl27dqFRYsWoVWrVvDw8MDMmTPh4eGBNWvW6JQ1NzeHs7OzuORMOJmX+/fvQ6lUok+fPmjcuDHatm2LqKgoAMC2bdtgamqKHj16vDG+06dPY/78+Vi6dCkWL16M5s2bw93dHe+//z527doFf3//Ap/r9u3b0bVrV51tbdq0wZgxY3S2de/eHQEBAeJrd3d3zJkzB/3794dSqUSlSpUQHBys8x6ZTIY1a9agU6dOsLS0RLVq1bBz506dMvfu3UOfPn1gZ2cHe3t7dOvWTWd4jZzmn3nz5sHV1RW1atV67fl8//33cHNzg5WVFfr06YPk5GRxn1arxezZs/HOO+/A3Nwc9evX16mFCwsLg0wm06kJi4yMhEwmE2MKCQmBnZ0dDh48CE9PT1hbW4sJbw6NRoOxY8fCzs4OFSpUwIQJE5DXdJxdu3bF9u3bX3s+b4vJDhER5Sk7OxsajSbX5IuWlpY4duyYzrawsDA4OjqiVq1aGDFiBB4/fpzvcd3c3KBWq3H+/Hk8efIEZ86cQd26dfH06VNMmzYt18zi+dmyZQusra0xcuTIPPfb2dkV6DgAcOzYMTRq1KjA5V+2ePFi1KtXD+fPn8fXX3+NL774AocPH9YpM23aNPTs2RMXLlzAgAED0K9fPzHBy8rKgp+fH8qVK4ejR4/i+PHjYvLwcg1OaGgooqOjcfjwYZ2asFfdvHkTv/76K/744w8cOHAA58+f17lGK1aswNKlS7FkyRJcvHgRfn5++PDDD3Hjxo1CnXd6ejqWLFmCzZs3499//0VsbCzGjRsn7l+6dClCQkKwfv16HDt2DE+ePMHu3btzHadJkyY4ffo01Gp1oT6/MJjsEBFRnsqVKweVSoU5c+bgwYMH0Gg0+PnnnxEeHq7zF3zHjh2xadMmhIaGYuHChThy5Ag6deoEjUaT53FtbGwwY8YMDB48GE2aNMGgQYPg5+eHcePGISgoCDExMfD19UWdOnVy1YC87MaNG6hWrRpMTU3f6jyTkpKQnJwMV1fXIr2/RYsW+Prrr1GzZk2MGjUKvXr1wrfffqtTpnfv3hg6dChq1qyJOXPmoFGjRli5ciUA4JdffoFWq8VPP/0EHx8feHp6YsOGDYiNjUVYWJh4DKVSiZ9++gne3t7w9vbON56MjAxs2rQJ9evXR6tWrbBy5Ups374d8fHxAIAlS5Zg4sSJ6NevH2rVqoWFCxeifv36WL58eaHOOysrC2vXrkWjRo3QoEEDBAUFITQ0VNy/fPlyTJo0CT169ICnpyfWrl2bZ78cV1dXZGZmivGVBHZQJoNVsWJFqUOgUoBNUCVr8+bN+PTTT1GpUiUoFAo0aNAA/fv3R0REhFimX79+4rqPjw/q1q2L6tWrIywsDO3atct1TBMTE7z//vsYO3YsFAoFAODIkSO4ePEiVq5cCQ8PD2zbtg3Ozs5o0qQJWrVqBUdHx1zHyatJpChyhq94tQaroFQqVa7XryYOeZWJjIwEAFy4cAE3b95EuXLldMpkZGTg1q1b4msfHx+YmZm9MZ7KlSujUqVKOp+l1WoRHR0NKysrPHjwIFffpBYtWuDChQtvPPbLrKysUL16dfG1i4sLEhMTAQDJycmIi4tD06ZNxf0mJiZo1KhRrn83S0tLAC9qikoKkx0ySEqlEg8fPpQ6DCKjV716dRw5cgRpaWlISUmBi4sL+vbti2rVquX7nmrVqqFixYq4efNmrmRHoVDkejRarVZj5MiR2Lx5M27evIns7Gy0bt0aAFCzZk2cOnUqV3+anH3Hjh1DVlbWW9XuVKhQATKZDE+fPtXZLpfLc/0wZ2VlFflz8pOamoqGDRtiy5YtufY5ODiI60qlstg/Oy85M9G/fO55nfer11wmkxUpAX3y5AkA3XMtbmzGIiKiN1IqlXBxccHTp09x8OBBdOvWLd+y//33Hx4/fgwXF5cCHXvu3Lno2LEjGjRoAI1Gg+zsbHFfVlZWvs1hH3/8MVJTU7F69eo897/uUfOXmZmZwcvLC1evXtXZ7uDgkKvD7eXLl3O9/+TJk7lee3p6FrhMgwYNcOPGDTg6OsLDw0NnKcrj2LGxsXjw4IHOZ8nlctSqVQs2NjZwdXXF8ePHdd5z/PhxeHl5iecNQOfcc2qhCsrW1hYuLi44deqUuC07O1unRjDH5cuX8c4775Robb2kyc7MmTMhk8l0ltq1a4v7MzIyEBgYiAoVKsDa2ho9e/ZEQkKCzjFiY2PRpUsXWFlZwdHREePHj9e5UYiIqOgOHjyIAwcOICYmBocPH0bbtm1Ru3ZtDB48GMCLWonx48fj5MmTuHPnDkJDQ9GtWzd4eHjAz8/vjce/evUqfvnlF8yePRsAULt2bcjlcqxbtw779u3DtWvXdMb0eVnTpk0xYcIEfPXVV5gwYQLCw8Nx9+5dhIaGonfv3rnG/nkdPz+/XJ2u33vvPezbt0+MY8SIEXkmUMePH8eiRYtw/fp1BAcHY8eOHfjiiy90yuzYsQPr16/H9evXMWPGDJw+fRpBQUEAgAEDBqBixYro1q0bjh49ipiYGISFhWH06NH477//CnwOOSwsLODv748LFy7g6NGjGD16NPr06QNnZ2cAwPjx47Fw4UL88ssviI6Oxtdff43IyEgxZg8PD7i5uWHmzJm4ceMG9u3bV6SB/7744gt888032LNnD65du4aRI0fmef2OHj1a4oM5St6M5e3tjb/++kt8bWLyfyF9+eWX2LdvH3bs2AFbW1sEBQWhR48eYkaq0WjQpUsXODs748SJE4iLi8OgQYNgamqK+fPn6/1cSH+eP3+OTp06AXgxOFlOmy+RoSnsIH/6lpycjEmTJuG///6Dvb09evbsiXnz5olNGAqFAhcvXsTGjRuRlJQEV1dXdOjQAXPmzIG5uXmu42m1WvGpHw8PDwwfPhzLli0Tm2gsLS0REhKCwMBAqNVqrFq1Sqf/yasWLlyIhg0bIjg4GGvXroVWq0X16tXRq1evQj16PmTIEDRq1AjJyclibcqnn36KCxcuYNCgQTAxMcGXX36Jtm3b5nrvV199hbNnz2LWrFmwsbHBsmXLciV6s2bNwvbt2zFy5Ei4uLhg27ZtYk2KlZUV/v33X0ycOBE9evTAs2fPUKlSJbRr165AYw29ysPDAz169EDnzp3x5MkTfPDBBzq1X6NHj0ZycjK++uorJCYmwsvLC7///jtq1KgB4EXz1LZt2zBixAjUrVsXjRs3xty5c9G7d+9CxfHVV18hLi4O/v7+kMvl+PTTT/HRRx/pPAafkZGBPXv2vHYAyuIgE4qrh1cRzJw5E3v27Mmzeiw5ORkODg7YunUrevXqBQC4du0aPD09ER4ejmbNmuHPP//EBx98gAcPHsDJyQkAsHbtWkycOBEPHz7MtyOXWq3WecQtJSUFbm5uSE5OLtIXi/QvLS0N1tbWAF78ZamvtmxDJmVHXtP0TEzwfPGX+6Ko6ciyenMnS30ryaQjIyMDMTExqFq1apE7wRoLjUaD8+fPAwB8fX3FDsqlQe/evdGgQQNMmjSpwO9xd3fHmDFjco3H8zKZTIbdu3frTJNAL6xZswa7d+/GoUOH8i3zuvsnJSUFtra2b/z9lrxm58aNG3B1dYWFhQVUKhUWLFiAypUrIyIiAllZWWjfvr1Ytnbt2qhcubKY7ISHh8PHx0dMdIAXVZEjRozAlStX4Ovrm+dnLliwALNmzSrxcyPSFz6RRIZCLpeLnZtzOsKWFosXL8Yff/whdRhliqmpqfgIfkmS9JvWtGlThISE4MCBA1izZg1iYmLw7rvv4tmzZ4iPj4eZmVmuQaGcnJzEZ/Hj4+N1Ep2c/Tn78jNp0iQkJyeLy71794r3xIiIKE8ymQz29vawt7eHTCaTOhwd7u7uGDVqlNRhlClDhw5942jQxUHSmp2cPhcAULduXTRt2hRVqlTBr7/+WqJ9MMzNzfNsSyYiIiqMl6d0yI+EvUXo/ytVdYh2dnaoWbMmbt68CWdnZ2RmZubquZ2QkCD2KHd2ds71dFbO65wyRERUegiCgCdPnuDJkydMAkhvSlWyk5qailu3bsHFxQUNGzaEqampztDT0dHRiI2NFUeiVKlUuHTpkjhiIwAcPnwYNjY2Yi93IqLSgj/uL57Gun37Nm7fvg2tVit1OGQAiuO+kbQZa9y4cejatSuqVKmCBw8eYMaMGVAoFOjfvz9sbW0xZMgQjB07Fvb29rCxscGoUaOgUqnQrFkzAECHDh3g5eWFgQMHYtGiRYiPj8fUqVMRGBjIZqoywMrKSuoQiAok5zHt9PR0DpNAVEg500i8zSjZkiY7//33H/r374/Hjx/DwcEBLVu2xMmTJ8XRG7/99lvI5XL07NkTarUafn5+OmMFKBQK7N27FyNGjIBKpYJSqYS/v784OBUZL6VSibS0NKnDICoQhUIBOzs7sRbaysqq1HXO1ZeXR0POyMgoVY+eU+kiCALS09ORmJgIOzu7t/quSJrsbN++/bX7LSwsEBwcjODg4HzLVKlSBfv37y/u0IiIilVOP8KXm93LIq1Wi0ePHgF40bm3tD1+TqWPnZ3dW/fDlXycHSKiskAmk8HFxQWOjo4lMpmkoUhPT0eXLl0AAOfOnWNzNL2WqalpsdT+Mdkhg5SRkYGePXsCAHbt2lXmR6Ulw6FQKMp0041Go8Hdu3cBvBgGhPcu6QOTHTJIGo1GbL7Mb0ZkIiIioJQ9ek5ERERU3JjsEBERkVFjskNERERGjckOERERGTV2UCYqxebFzpM6BCIig8eaHSIiIjJqrNkhg6RUKjmpIpEB4r1LUmDNDhERERk1JjtERERk1JjskEHKyMhA79690bt3b2RkZEgdDhEVEO9dkgKTHTJIGo0GO3fuxM6dOzldBJEB4b1LUmAHZSIi0hszMzOsWrVKXCfSByY7RESkN6ampggMDJQ6DCpj2IxFRERERo01O0REpDcajQZHjx4FALz77rtQKBQSR0RlAZMdIiLSm4yMDLRt2xYAkJqaCqVSKXFEVBawGYuIiIiMGmt2yCBZWVkhNTVVXCciIsoPkx0ySDKZjNXfRERUIGzGIiIiIqPGZIcMklqtRkBAAAICAqBWq6UOh4iISjEmO2SQsrOzsXHjRmzcuBHZ2dlSh0NERKUYkx0iIiIyakx2iIiIyKgx2SEiIiKjxmSHiIiIjBqTHSIiIjJqHFSQiAjAvNh5bywzpfIUPURCRMWNyQ4ZJCsrKyQmJorrRGQYeO+SFJjskEGSyWRwcHCQOgwiKiTeuyQF9tkhIiIio8ZkhwySWq1GYGAgAgMDOV0EkQHhvUtSYLJDBik7OxurV6/G6tWrOV0EkQHhvUtSYJ8dIiLSG1NTU8yYMUNcJ9IHJjtERKQ3ZmZmmDlzptRhUBnDZIdIIgUZ14WIiN4ekx0iItIbrVaLqKgoAICnpyfkcnYdpZLHZIeIiPTm+fPnqFOnDgAgNTUVSqVS4oioLGBKTUREREaNNTtkkCwtLRETEyOuExER5YfJDhkkuVwOd3d3qcMgIiIDwGYsIiIiMmqlJtn55ptvIJPJMGbMGHFbRkYGAgMDUaFCBVhbW6Nnz55ISEjQeV9sbCy6dOkCKysrODo6Yvz48RyVswzIzMzE+PHjMX78eGRmZkodDhERlWKlItk5c+YMvv/+e9StW1dn+5dffok//vgDO3bswJEjR/DgwQP06NFD3K/RaNClSxdkZmbixIkT2LhxI0JCQjB9+nR9nwLpWVZWFpYsWYIlS5YgKytL6nCIiKgUkzzZSU1NxYABA/Djjz+ifPny4vbk5GSsW7cOy5Ytw3vvvYeGDRtiw4YNOHHiBE6ePAkAOHToEK5evYqff/4Z9evXR6dOnTBnzhwEBwe/9q99tVqNlJQUnYWIiIiMk+TJTmBgILp06YL27dvrbI+IiEBWVpbO9tq1a6Ny5coIDw8HAISHh8PHxwdOTk5iGT8/P6SkpODKlSv5fuaCBQtga2srLm5ubsV8VkRERFRaSJrsbN++HefOncOCBQty7YuPj4eZmRns7Ox0tjs5OSE+Pl4s83Kik7M/Z19+Jk2ahOTkZHG5d+/eW54JERERlVaSPXp+7949fPHFFzh8+DAsLCz0+tnm5uYwNzfX62cSERGRNCRLdiIiIpCYmIgGDRqI2zQaDf7991+sWrUKBw8eRGZmJpKSknRqdxISEuDs7AwAcHZ2xunTp3WOm/O0Vk4ZIqLiUpDJW6dUnqKHSIioMCRrxmrXrh0uXbqEyMhIcWnUqBEGDBggrpuamiI0NFR8T3R0NGJjY6FSqQAAKpUKly5dQmJioljm8OHDsLGxgZeXl97PiYiIiEofyWp2ypUrJ04Gl0OpVKJChQri9iFDhmDs2LGwt7eHjY0NRo0aBZVKhWbNmgEAOnToAC8vLwwcOBCLFi1CfHw8pk6disDAQDZTGTlLS0tcvnxZXCciw8B7l6RQqqeL+PbbbyGXy9GzZ0+o1Wr4+flh9erV4n6FQoG9e/dixIgRUKlUUCqV8Pf3x+zZsyWMmvRBLpfD29tb6jCIqJB475IUSlWyExYWpvPawsICwcHBCA4Ozvc9VapUwf79+0s4MiIiIjJUpSrZISqozMxMzJ8/HwAwefJkmJmZSRwRERUE712SApMdMkhZWVmYNWsWAGD8+PH8D5PIQPDeJSkw2SEiIr0xMTHByJEjxXUifeA3jYiI9Mbc3Py1/TCJSoLkc2MRERERlSTW7BARkd4IgoBHjx4BACpWrAiZTCZxRFQWMNkhIiK9SU9Ph6OjIwAgNTUVSqVS4oioLGAzFhERERk11uyQQbKwsBAngbWwsJA4GiIiKs2Y7JBBUigUaNy4sdRh5Ksgs2MTEZF+sBmLiIiIjBprdsggZWZmYsWKFQCAL774gqOwEhFRvopUs3Pu3DlcunRJfP2///0P3bt3x+TJk5GZmVlswRHlJysrCxMmTMCECROQlZUldThERFSKFSnZ+eyzz3D9+nUAwO3bt9GvXz9YWVlhx44dmDBhQrEGSERERPQ2ipTsXL9+HfXr1wcA7NixA61atcLWrVsREhKCXbt2FWd8RERERG+lSMmOIAjQarUAgL/++gudO3cGALi5uYkjYxIRERGVBkVKdho1aoS5c+di8+bNOHLkCLp06QIAiImJgZOTU7EGSERERPQ2ipTsfPvttzh37hyCgoIwZcoUeHh4AAB27tyJ5s2bF2uARERERG+jSI+e16tXT+dprByLFy+GiQmfZiciIqLSo0iZSbVq1XDmzBlUqFBBZ3tGRgYaNGiA27dvF0twRPmxsLDAP//8I64TkWHgvUtSKFKyc+fOHWg0mlzb1Wo1/vvvv7cOiuhNFAoF2rRpI3UYRFRIvHdJCoVKdn7//Xdx/eDBg7C1tRVfazQahIaGomrVqsUXHREREdFbKlSy0717dwCATCaDv7+/zj5TU1O4u7tj6dKlxRYcUX6ysrLwww8/AACGDx8OU1NTiSMiooLgvUtSKFSykzO2TtWqVXHmzBlUrFixRIIiepPMzEwEBQUBAAICAvgfJpGB4L1LUihSn52YmJjijoOIiMoAhUKBXr16ietE+lDk58RDQ0MRGhqKxMREscYnx/r16986MCIiMj4WFhbYsWOH1GFQGVOkZGfWrFmYPXs2GjVqBBcXF8hksuKOi4iIiKhYFCnZWbt2LUJCQjBw4MDijoeIyKDNi533xjJTKk/RQyRElKNI00VkZmZyWggiIiq0tLQ0yGQyyGQypKWlSR0OlRFFSnaGDh2KrVu3FncsRERERMWuSM1YGRkZ+OGHH/DXX3+hbt26uR4dXLZsWbEER5Qfc3Nz7N27V1wnIiLKT5GSnYsXL6J+/foAgMuXL+vsY2dl0gcTExN06dJF6jCIiMgAFCnZyZnEjYiIiKi0K/I4O0RSysrKwpYtWwAAAwYM4CisRESUryIlO23btn1tc9Xff/9d5ICICiIzMxODBw8GAPTu3ZvJDhER5atIyU5Of50cWVlZiIyMxOXLl3NNEEpEREQkpSIlO99++22e22fOnInU1NS3CoiIiIioOBVpnJ38fPLJJ5wXi4iIiEqVYk12wsPDYWFhUZyHJCIiInorRWrG6tGjh85rQRAQFxeHs2fPYtq0acUSGBEREVFxKFKyY2trq/NaLpejVq1amD17Njp06FAsgREREREVhyIlOxs2bCjuOIgKxdzcHL/++qu4TkSGgfcuSeGtBhWMiIhAVFQUAMDb2xu+vr7FEhTRm5iYmKB3795Sh0FEhcR7l6RQpGQnMTER/fr1Q1hYGOzs7AAASUlJaNu2LbZv3w4HB4fijJGIiIioyIqU7IwaNQrPnj3DlStX4OnpCQC4evUq/P39MXr0aGzbtq1YgyR6VXZ2Nnbv3g0A+Oijj2Bior+ZT+bFztPbZxEZGynvXSq7ivTo+YEDB7B69Wox0QEALy8vBAcH488//yzwcdasWYO6devCxsYGNjY2UKlUOu/PyMhAYGAgKlSoAGtra/Ts2RMJCQk6x4iNjUWXLl1gZWUFR0dHjB8/HtnZ2UU5LTIgarUaffr0QZ8+faBWq6UOh4gKiPcuSaFIKbVWq81zLiJTU1NotdoCH+edd97BN998gxo1akAQBGzcuBHdunXD+fPn4e3tjS+//BL79u3Djh07YGtri6CgIPTo0QPHjx8HAGg0GnTp0gXOzs44ceIE4uLiMGjQIJiammL+/PlFOTUiIipBcrkcrVu3FteJ9EEmCIJQ2Dd169YNSUlJ2LZtG1xdXQEA9+/fx4ABA1C+fHmxirIo7O3tsXjxYvTq1QsODg7YunUrevXqBQC4du0aPD09ER4ejmbNmuHPP//EBx98gAcPHsDJyQkAsHbtWkycOBEPHz6EmZlZgT4zJSUFtra2SE5Oho2NTZFjJ/1JS0uDtbU1ACA1NRVKpVJvn81mrMIzTc/EBM/ZAIBFUdORZVWwe9NYTak8ReoQiIxCQX+/i5RWr1q1CikpKXB3d0f16tVRvXp1VK1aFSkpKVi5cmWRAtZoNNi+fTvS0tKgUqkQERGBrKwstG/fXixTu3ZtVK5cGeHh4QBejNjs4+MjJjoA4Ofnh5SUFFy5ciXfz1Kr1UhJSdFZiIiIyDgVqRnLzc0N586dw19//YVr164BADw9PXUSk4K6dOkSVCoVMjIyYG1tjd27d8PLywuRkZEwMzMTn/bK4eTkhPj4eABAfHy8TqKTsz9nX34WLFiAWbNmFTpWIiIiMjyFqtn5+++/4eXlhZSUFMhkMrz//vsYNWoURo0ahcaNG8Pb2xtHjx4tVAC1atVCZGQkTp06hREjRsDf3x9Xr14t1DEKa9KkSUhOThaXe/fulejnERHRC2lpaXBwcICDgwPS0tKkDofKiELV7CxfvhzDhg3Ls13M1tYWn332GZYtW4Z33323wMc0MzODh4cHAKBhw4Y4c+YMVqxYgb59+yIzMxNJSUk6tTsJCQlwdnYGADg7O+P06dM6x8t5WiunTF7Mzc05cicRkUQePXokdQhUxhSqZufChQvo2LFjvvs7dOiAiIiItwpIq9VCrVajYcOGMDU1RWhoqLgvOjoasbGxUKlUAACVSoVLly4hMTFRLHP48GHY2NjAy8vrreKg0s3MzAwbNmzAhg0bCtwRnYiIyqZC1ewkJCTk+ci5eDATEzx8+LDAx5s0aRI6deqEypUr49mzZ9i6dSvCwsJw8OBB2NraYsiQIRg7dizs7e1hY2ODUaNGQaVSoVmzZgBeJFdeXl4YOHAgFi1ahPj4eEydOhWBgYGsuTFypqamCAgIkDoMIiIyAIVKdipVqoTLly+LzU6vunjxIlxcXAp8vMTERAwaNAhxcXGwtbVF3bp1cfDgQbz//vsAgG+//RZyuRw9e/aEWq2Gn58fVq9eLb5foVBg7969GDFiBFQqFZRKJfz9/TF79uzCnBYREREZsUIlO507d8a0adPQsWNHWFhY6Ox7/vw5ZsyYgQ8++KDAx1u3bt1r91tYWCA4OBjBwcH5lqlSpQr2799f4M8k45CdnY2DBw8CeDHcAIecJyKi/BTqF2Lq1Kn47bffULNmTQQFBaFWrVoAXgz2FxwcDI1GgylTOFgWlTy1Wi0m1qmpqUx2iIgoX4X6hXBycsKJEycwYsQITJo0CTmDL8tkMvj5+SE4ODjXuDdEREREUir0n8M5zUZPnz7FzZs3IQgCatSogfLly5dEfERERERvpch1/+XLl0fjxo2LMxYiIiKiYseODkREelaQyWQ5WShR8SnSRKBEREREhoLJDhERERk1NmORQTIzM8OqVavEdSIyDLx3SQpMdsggmZqaIjAwUOowiKiQeO+SFNiMRUREREaNNTtkkDQaDY4ePQoAePfdd6FQKCSOiIgKgvcuSYHJDhmkjIwMtG3bFsCL6SKUSqXEERFRQfDeJSkw2SEiIr2RyWTw8vIS14n0gckOERHpjZWVFa5cuSJ1GFTGsIMyERERGTUmO0RERGTUmOwQEZHepKenw9vbG97e3khPT5c6HCoj2GeHiIj0RhAEXL16VVwn0gcmO2SQTE1NsWjRInGdiIgoP0x2yCCZmZlh/PjxUodBREQGgH12iIiIyKixZocMkkajwblz5wAADRo04JDzRESULyY7ZJAyMjLQpEkTABxynoiIXo/NWERERGTUmOwQERGRUWOyQ0REREaNyQ4REREZNSY7REREZNSY7BAREZFR46PnZJBMTU0xY8YMcZ2IDAPvXZICkx0ySGZmZpg5c6bUYRBRIfHeJSmwGYuIiIiMGmt2yCBptVpERUUBADw9PSGXM28nMgS8d0kKTHbIID1//hx16tQBwOkiiAwJ712SApMdIiLSq4oVK0odApUxTHaIiEhvlEolHj58KHUYVMYw2SEiKoXmxc57Y5kplafoIRIiw8eeYURERGTUmOwQEZHePH/+HG3atEGbNm3w/PlzqcOhMoLNWEREpDdarRZHjhwR14n0gckOGSRTU1OMGzdOXC8uBeknQUREhoXJDhkkMzMzLF68WOowiIjIALDPDhERERk11uyQQdJqtYiNjQUAVK5cmUPOExFRvpjskEF6/vw5qlatCoBDzhMR0etJ+ufwggUL0LhxY5QrVw6Ojo7o3r07oqOjdcpkZGQgMDAQFSpUgLW1NXr27ImEhASdMrGxsejSpQusrKzg6OiI8ePHIzs7W5+nQkRERKWUpMnOkSNHEBgYiJMnT+Lw4cPIyspChw4dkJaWJpb58ssv8ccff2DHjh04cuQIHjx4gB49eoj7NRoNunTpgszMTJw4cQIbN25ESEgIpk+fLsUpERERUSkjaTPWgQMHdF6HhITA0dERERERaNWqFZKTk7Fu3Tps3boV7733HgBgw4YN8PT0xMmTJ9GsWTMcOnQIV69exV9//QUnJyfUr18fc+bMwcSJEzFz5kyYmZlJcWpERERUSpSqXp3JyckAAHt7ewBAREQEsrKy0L59e7FM7dq1UblyZYSHhwMAwsPD4ePjAycnJ7GMn58fUlJScOXKlTw/R61WIyUlRWchIiIi41Rqkh2tVosxY8agRYsWqFOnDgAgPj4eZmZmsLOz0ynr5OSE+Ph4sczLiU7O/px9eVmwYAFsbW3Fxc3NrZjPhoiIiEqLUvM0VmBgIC5fvoxjx46V+GdNmjQJY8eOFV+npKQw4SEig8OZ0YkKplQkO0FBQdi7dy/+/fdfvPPOO+J2Z2dnZGZmIikpSad2JyEhAc7OzmKZ06dP6xwv52mtnDKvMjc3h7m5eTGfBemTiYkJRo4cKa4TkWHgvUtSkPSbJggCRo0ahd27dyMsLEwcNyVHw4YNYWpqitDQUPTs2RMAEB0djdjYWKhUKgCASqXCvHnzkJiYCEdHRwDA4cOHYWNjAy8vL/2eEOmNubk5goODpQ6DiAqJ9y5JQdJkJzAwEFu3bsX//vc/lCtXTuxjY2trC0tLS9ja2mLIkCEYO3Ys7O3tYWNjg1GjRkGlUqFZs2YAgA4dOsDLywsDBw7EokWLEB8fj6lTpyIwMJC1N0RERCRtsrNmzRoAQJs2bXS2b9iwAQEBAQCAb7/9FnK5HD179oRarYafnx9Wr14tllUoFNi7dy9GjBgBlUoFpVIJf39/zJ49W1+nQRIQBAGPHj0CAFSsWBEymUziiIioIHjvkhQkb8Z6EwsLCwQHB7+22rNKlSrYv39/cYZGpVx6errYbMnpIogMB+9dkkKpefSciIiIqCSwKzwREemNUqksUK0+UXFizQ4REREZNSY7REREZNSY7BARkd5kZGSgd+/e6N27NzIyMqQOh8oIJjtERKQ3Go0GO3fuxM6dO6HRaKQOh8oIdlAmg2RiYgJ/f39xnYiIKD/8lSCDZG5ujpCQEKnDICIiA8BmLCIiIjJqrNkhgyQIAtLT0wEAVlZWHHKeiIjyxZodMkjp6emwtraGtbW1mPQQERHlhckOERERGTU2Y1GZMS92ntQhEBGRBFizQ0REREaNyQ4REREZNSY7REREZNSY7BAREZFRYwdlMkgKhQK9evUS19n5mMgwvHrvEukDkx0ySBYWFtixY4fUYRBRIfHeJSmwGYuIiIiMGmt2iIiMWEGaeKdUnqKHSIikw5odMkhpaWmQyWSQyWRIS0uTOhwiKiDeuyQFJjtERERk1NiMRUREemNlZYXExERxnUgfmOwQEZHeyGQyODg4SB0GlTFsxiIiIiKjxmSHiIj0Rq1WIzAwEIGBgVCr1VKHQ2UEkx0iItKb7OxsrF69GqtXr0Z2drbU4VAZwT47ZJAUCgU6d+4srhMREeWHyQ4ZJAsLC+zbt0/qMIiMAgceJGPHZiwiIiIyakx2iIiIyKgx2SGDlJaWBqVSCaVSySHniYjotdhnhwxWenq61CEQEZEBYM0OERERGTUmO0RERGTUmOwQERGRUWOyQ0REREaNyQ4REREZNT6NRQZJLpejdevW4joRGQbeuyQFJjtkkCwtLREWFiZ1GERUSLx3SQpMq4mIiMioMdkhIiIio8ZmLDJIaWlpcHd3BwDcuXNH0liIyoLimhn91XtXqVS+bWhEb8RkhwzWo0ePpA6BiIqA9y7pG5MdIiLSG0tLS1y+fFlcJ9IHSfvs/Pvvv+jatStcXV0hk8mwZ88enf2CIGD69OlwcXGBpaUl2rdvjxs3buiUefLkCQYMGAAbGxvY2dlhyJAhSE1N1eNZEBFRQcnlcnh7e8Pb25uPnpPeSPpNS0tLQ7169RAcHJzn/kWLFuG7777D2rVrcerUKSiVSvj5+SEjI0MsM2DAAFy5cgWHDx/G3r178e+//2L48OH6OgUiIiIq5SRtxurUqRM6deqU5z5BELB8+XJMnToV3bp1AwBs2rQJTk5O2LNnD/r164eoqCgcOHAAZ86cQaNGjQAAK1euROfOnbFkyRK4urrmeWy1Wg21Wi2+TklJKeYzo+KUV8fIzPRMcX3RvUUwszLTZ0hEVESZmZmYP38+AGDy5MkwM+O9SyWv1NYhxsTEID4+Hu3btxe32draomnTpggPDwcAhIeHw87OTkx0AKB9+/aQy+U4depUvsdesGABbG1txcXNza3kToSIiERZWVmYNWsWZs2ahaysLKnDoTKi1CY78fHxAAAnJyed7U5OTuK++Ph4ODo66uw3MTGBvb29WCYvkyZNQnJysrjcu3evmKOnkiaTy1CpbiVUqlsJMrlM6nCIiKgUK5NPY5mbm8Pc3FzqMOgtmFqYYsQfI6QOg4iIDECprdlxdnYGACQkJOhsT0hIEPc5OzsjMTFRZ392djaePHkiliEiIqKyrdQmO1WrVoWzszNCQ0PFbSkpKTh16hRUKhUAQKVSISkpCREREWKZv//+G1qtFk2bNtV7zERERFT6SNqMlZqaips3b4qvY2JiEBkZCXt7e1SuXBljxozB3LlzUaNGDVStWhXTpk2Dq6srunfvDgDw9PREx44dMWzYMKxduxZZWVkICgpCv3798n0Si4xD5vNMfNf+OwDA6L9Gw8yST3QQEVHeJE12zp49i7Zt24qvx44dCwDw9/dHSEgIJkyYgLS0NAwfPhxJSUlo2bIlDhw4AAsLC/E9W7ZsQVBQENq1awe5XI6ePXviu+++0/u5kJ4JQNJ/SeI6ERFRfiRNdtq0aQNByP+XSiaTYfbs2Zg9e3a+Zezt7bF169aSCI+IiAqhIJOFjqkwpuQDIXpFqe2zQ0RERFQcmOwQERGRUWOyQ0REREaNyQ4REREZtTI5gjIZARngWMNRXCciwyCTyeDl5SWuE+kDkx0ySGaWZhj912ipwyCiQrKyssKVK1ekDoPKGDZjERERkVFjzQ5JqiDjchAREb0NJjtkkDKfZ2Jt17UAgM//+JzTRRAZiBnRM954706pPEXfYZGRY7JDhkkAEm8kiutEZCB475IEmOwQEZHemJib4NPtn4rrRPrAbxoREemNXCFHNVU1qcOgMoZPYxEREZFRY80OERHpjSZLgzNbzwAAGn/cGApThcQRUVnAZIeIiPRGk6XB3ul7AQANejdgskN6wWSHDJMMsHvHTlwnIiLKD5MdMkhmlmYYd3yc1GEQEZEBYAdlIiIiMmpMdoiIiMiosRmLDFJWRhZ+6v0TAGDojqEwtTCVOCIiIiqtmOyQQRK0Au5fvC+uExER5YfNWERERGTUmOwQERGRUWOyQ0REREaNyQ4REREZNXZQphIzL3ae1CEQkQEqyP8dUypP0UMkZCyY7JDBsrK3kjoEIioC3rukb0x2yCCZWZlh8vnJUodBRIXEe5ekwGSHiIgMDpu6qDDYQZmIiIiMGmt2yCBlZWRho/9GAID/Rn9OF0FkIPR577L2h3Iw2SGDJGgF3Dl5R1wnIsPAe5ekwGSHioSPlRNRUSjMFOi3up+4TqQPTHaIiEhvFCYK1OlSR+owqIxhB2UiIiIyaqzZISIivdFkaxB1MAoA4OnnCYUJm7Ko5LFmh4iI9EaTqcH2kduxfeR2aDI1UodDZQRrdshgmVrycXMiInozJjuUiyE8aWVmZYYZ12ZIHQYRERkAJjtERFRmceDBsoHJDhER0WswITJ8THbIIGVlZGHb59sAAP3X9ud0EURElC8mO2SQBK2A6/9cF9eJiKTE2p/SjckOERGRHjAhkg6THSIiolKCCVHJMJpkJzg4GIsXL0Z8fDzq1auHlStXokmTJlKHpVeG8Mg4ERGRvhlFsvPLL79g7NixWLt2LZo2bYrly5fDz88P0dHRcHR0lDo8IiKiYlNcf9iWpRoio5guYtmyZRg2bBgGDx4MLy8vrF27FlZWVli/fr3UoREREZHEDL5mJzMzExEREZg0aZK4TS6Xo3379ggPD8/zPWq1Gmq1WnydnJwMAEhJSSn2+BbfW/zGMuPdxhfLccqSrOdZ4ro6VQ2tRithNPQmmudZyLm7MlLVyOK/V5nFe7f0mHZlmt4+qyC/c0WR87stCG94KlcwcPfv3xcACCdOnNDZPn78eKFJkyZ5vmfGjBkCAC5cuHDhwoWLESz37t17ba5g8DU7RTFp0iSMHTtWfK3VavHkyRNUqFABMplMwsgMW0pKCtzc3HDv3j3Y2NhIHY5R47XWD15n/eB11g9jvM6CIODZs2dwdXV9bTmDT3YqVqwIhUKBhIQEne0JCQlwdnbO8z3m5uYwNzfX2WZnZ1dSIZY5NjY2RnMjlXa81vrB66wfvM76YWzX2dbW9o1lDL6DspmZGRo2bIjQ0FBxm1arRWhoKFQqlYSRERERUWlg8DU7ADB27Fj4+/ujUaNGaNKkCZYvX460tDQMHjxY6tCIiIhIYkaR7PTt2xcPHz7E9OnTER8fj/r16+PAgQNwcnKSOrQyxdzcHDNmzMjVREjFj9daP3id9YPXWT/K8nWWCcKbntciIiIiMlwG32eHiIiI6HWY7BAREZFRY7JDRERERo3JDhERERk1JjtUKMHBwXB3d4eFhQWaNm2K06dPF+h927dvh0wmQ/fu3Us2QCNSmGsdEhICmUyms1hYWOgxWsNV2O90UlISAgMD4eLiAnNzc9SsWRP79+/XU7SGqzDXuU2bNrm+zzKZDF26dNFjxIapsN/n5cuXo1atWrC0tISbmxu+/PJLZGRk6ClaPSqeGaqoLNi+fbtgZmYmrF+/Xrhy5YowbNgwwc7OTkhISHjt+2JiYoRKlSoJ7777rtCtWzf9BGvgCnutN2zYINjY2AhxcXHiEh8fr+eoDU9hr7NarRYaNWokdO7cWTh27JgQExMjhIWFCZGRkXqO3LAU9jo/fvxY57t8+fJlQaFQCBs2bNBv4AamsNd5y5Ytgrm5ubBlyxYhJiZGOHjwoODi4iJ8+eWXeo685DHZoQJr0qSJEBgYKL7WaDSCq6ursGDBgnzfk52dLTRv3lz46aefBH9/fyY7BVTYa71hwwbB1tZWT9EZj8Je5zVr1gjVqlUTMjMz9RWiUSjK/x0v+/bbb4Vy5coJqampJRWiUSjsdQ4MDBTee+89nW1jx44VWrRoUaJxSoHNWFQgmZmZiIiIQPv27cVtcrkc7du3R3h4eL7vmz17NhwdHTFkyBB9hGkUinqtU1NTUaVKFbi5uaFbt264cuWKPsI1WEW5zr///jtUKhUCAwPh5OSEOnXqYP78+dBoNPoK2+AU9fv8snXr1qFfv35QKpUlFabBK8p1bt68OSIiIsSmrtu3b2P//v3o3LmzXmLWJ6MYQZlK3qNHj6DRaHKNSu3k5IRr167l+Z5jx45h3bp1iIyM1EOExqMo17pWrVpYv3496tati+TkZCxZsgTNmzfHlStX8M477+gjbINTlOt8+/Zt/P333xgwYAD279+PmzdvYuTIkcjKysKMGTP0EbbBKcp1ftnp06dx+fJlrFu3rqRCNApFuc4ff/wxHj16hJYtW0IQBGRnZ+Pzzz/H5MmT9RGyXrFmh0rEs2fPMHDgQPz444+oWLGi1OEYPZVKhUGDBqF+/fpo3bo1fvvtNzg4OOD777+XOjSjotVq4ejoiB9++AENGzZE3759MWXKFKxdu1bq0IzWunXr4OPjgyZNmkgditEJCwvD/PnzsXr1apw7dw6//fYb9u3bhzlz5kgdWrFjzQ4VSMWKFaFQKJCQkKCzPSEhAc7OzrnK37p1C3fu3EHXrl3FbVqtFgBgYmKC6OhoVK9evWSDNlCFvdZ5MTU1ha+vL27evFkSIRqFolxnFxcXmJqaQqFQiNs8PT0RHx+PzMxMmJmZlWjMhuhtvs9paWnYvn07Zs+eXZIhGoWiXOdp06Zh4MCBGDp0KADAx8cHaWlpGD58OKZMmQK53HjqQ4znTKhEmZmZoWHDhggNDRW3abVahIaGQqVS5Spfu3ZtXLp0CZGRkeLy4Ycfom3btoiMjISbm5s+wzcohb3WedFoNLh06RJcXFxKKkyDV5Tr3KJFC9y8eVNM3AHg+vXrcHFxYaKTj7f5Pu/YsQNqtRqffPJJSYdp8IpyndPT03MlNDmJvGBs02ZK3UOaDMf27dsFc3NzISQkRLh69aowfPhwwc7OTnzEeeDAgcLXX3+d7/v5NFbBFfZaz5o1Szh48KBw69YtISIiQujXr59gYWEhXLlyRapTMAiFvc6xsbFCuXLlhKCgICE6OlrYu3ev4OjoKMydO1eqUzAIRf2/o2XLlkLfvn31Ha7BKux1njFjhlCuXDlh27Ztwu3bt4VDhw4J1atXF/r06SPVKZQYNmNRgfXt2xcPHz7E9OnTER8fj/r16+PAgQNih7jY2FijqvaUUmGv9dOnTzFs2DDEx8ejfPnyaNiwIU6cOAEvLy+pTsEgFPY6u7m54eDBg/jyyy9Rt25dVKpUCV988QUmTpwo1SkYhKL83xEdHY1jx47h0KFDUoRskAp7nadOnQqZTIapU6fi/v37cHBwQNeuXTFv3jypTqHEyATB2OqqiIiIiP4P/wwnIiIio8Zkh4iIiIwakx0iIiIyakx2iIiIyKgx2SEiIiKjxmSHiIiIjBqTHSIiIjJqTHaIiIjIqDHZIaIyKSwsDDKZDElJSVKHQkQljMkOEZVqAQEBkMlkkMlkMDU1RdWqVTFhwgRkZGQU+Bht2rTBmDFjdLY1b94ccXFxsLW1LeaIiai04dxYRFTqdezYERs2bEBWVhYiIiLg7+8PmUyGhQsXFvmYZmZmcHZ2LsYoiai0Ys0OEZV65ubmcHZ2hpubG7p374727dvj8OHDAIDHjx+jf//+qFSpEqysrODj44Nt27aJ7w0ICMCRI0ewYsUKsYbozp07eTZj7dq1C97e3jA3N4e7uzuWLl2q71MlohLAZIeIDMrly5dx4sQJmJmZAQAyMjLQsGFD7Nu3D5cvX8bw4cMxcOBAnD59GgCwYsUKqFQqDBs2DHFxcYiLi4Obm1uu40ZERKBPnz7o168fLl26hJkzZ2LatGkICQnR5+kRUQlgMxYRlXp79+6FtbU1srOzoVarIZfLsWrVKgBApUqVMG7cOLHsqFGjcPDgQfz6669o0qQJbG1tYWZmBisrq9c2Wy1btgzt2rXDtGnTAAA1a9bE1atXsXjxYgQEBJTo+RFRyWKyQ0SlXtu2bbFmzRqkpaXh22+/hYmJCXr27AkA0Gg0mD9/Pn799Vfcv38fmZmZUKvVsLKyKtRnREVFoVu3bjrbWrRogeXLl0Oj0UChUBTb+RCRfrEZi4hKPaVSCQ8PD9SrVw/r16/HqVOnsG7dOgDA4sWLsWLFCkycOBH//PMPIiMj4efnh8zMTImjJqLSgskOERkUuVyOyZMnY+rUqXj+/DmOHz+Obt264ZNPPkG9evVQrVo1XL9+Xec9ZmZm0Gg0rz2up6cnjh8/rrPt+PHjqFmzJmt1iAwckx0iMji9e/eGQqFAcHAwatSogcOHD+PEiROIiorCZ599hoSEBJ3y7u7uOHXqFO7cuYNHjx5Bq9XmOuZXX32F0NBQzJkzB9evX8fGjRuxatUqnf5ARGSYmOwQkcExMTFBUFAQFi1ahK+++goNGjSAn58f2rRpA2dnZ3Tv3l2n/Lhx46BQKODl5QUHBwfExsbmOmaDBg3w66+/Yvv27ahTpw6mT5+O2bNns3MykRGQCYIgSB0EERERUUlhzQ4REREZNSY7REREZNSY7BAREZFRY7JDRERERo3JDhERERk1JjtERERk1JjsEBERkVFjskNERERGjckOERERGTUmO0RERGTUmOwQERGRUft/RFQyzd57qAQAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.figure()\n",
        "plt.hist(ratio_boot, bins=50, color=cMAP1[0,], alpha=0.6)\n",
        "plt.axvline(x=ratio_empirical, color='r', linestyle='-')\n",
        "plt.axvline(x=lower_bound, color='k', linestyle='--')\n",
        "plt.axvline(x=upper_bound, color='k', linestyle='-.')\n",
        "plt.xlabel('Ratio'); plt.ylabel('Counts');\n",
        "plt.legend(['Bootstrapped ratios', 'Empirical ratio', '95% CI (lower bound)',\n",
        "    '95 % CI (upper bound)'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xX-phmNyLlUS"
      },
      "source": [
        "The fact that the 95% confidence interval does not include 1 implies that we can be very confident that the aspirin therapy is effective in preventing heart attacks."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zKDOcBcULlUS"
      },
      "source": [
        "# Part 2 - Cross validation (leave-one-out)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GpzuTiKxLlUS"
      },
      "source": [
        "Cross-validation is a resampling method that uses different portions of the data to test and train a model on different iterations. The goal of  cross-validation is to test the model's ability to predict new data that was not used in estimating it, in order to flag problems like overfitting and to give an insight on how the model will generalize to a new (or unknown) dataset."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 393
        },
        "id": "JDJW2Q7VLlUS",
        "outputId": "9885426f-f686-4628-d4b7-1f270f86828d"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAF4CAYAAAD9rzvUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7u0lEQVR4nO3de3BbZX438O+5SEeyHcnxNfD6OBvDhrzEAdJw2bAdNlmyu9Ab3r5l+nqabaAMtDQwZbLTFtrpMrTdDTswCx3KsKFsgYFi9taYvrSUZimBtiwLBLJNnA1pEKzlkIsvseQo0pHOOc/7h1bCdmRZsiwfnaPvp6PZ2pGUR8bR+eq5/H6SEEKAiIiIqEpkpwdARERE3sawQURERFXFsEFERERVxbBBREREVcWwQURERFXFsEFERERVxbBBREREVcWwQURERFXFsEFERERVxbBBREREVeXasHH//fdDkiTcddddTg+FiIiIinBl2Hj77bexa9cuXHLJJU4PhYiIiObhurBx5swZ/M7v/A7+7u/+DsuXL3d6OERERDQP1ekBlGv79u341V/9VWzZsgV//dd/XfS+hmHAMIz817ZtY2JiAq2trZAkqdpDJSIi8gwhBKampnD++edDlsubq3BV2Hj++efx7rvv4u233y7p/jt37sR9991X5VERERHVj2g0iq6urrIe45qwEY1G8Ud/9EfYs2cPAoFASY+55557sGPHjvzXsVgM3d3diEajCIVC1RoqERGR58Tjcei6jmXLlpX9WEkIIaowpkU3ODiIL3/5y1AUJf89y7IgSRJkWYZhGDP+rJB4PI5wOIxYLMawQUREVIZKrqGumdm49tprceDAgRnfu/nmm7FmzRr86Z/+6bxBg4iIiJzhmrCxbNky9Pb2zvheY2MjWltbz/k+ERER1Q7XHX0lIiIid3HNzEYhe/fudXoIRERENA/ObBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVW5Jmw89thjuOSSSxAKhRAKhbBx40a89NJLTg+LiIiI5uGasNHV1YX7778f+/btwzvvvIPPf/7zuOGGGzA0NOT00IiIiKgISQghnB7EQrW0tOCBBx7ALbfcUtL94/E4wuEwYrEYQqFQlUdHRETkHZVcQ9UqjamqLMvC97//fSQSCWzcuHHO+xmGAcMw8l/H4/GlGB4RERFN45plFAA4cOAAmpqaoGka/uAP/gC7d+/GxRdfPOf9d+7ciXA4nL/pur6EoyUiIiLAZcso6XQaw8PDiMVi+MEPfoAnnngCr7322pyBo9DMhq7rXEYhIiIqUyXLKK4KG7Nt2bIFF1xwAXbt2lXS/blng4iIaGEquYa6ahllNtu2Z8xcEBERUe1xzQbRe+65B9dffz26u7sxNTWF5557Dnv37sXLL7/s9NCIiIioCNeEjVOnTuF3f/d3cfz4cYTDYVxyySV4+eWX8YUvfMHpoREREVERrgkb3/nOd5weAhERES2Aa8IGERHRYhFCICVSMG0TqqwiIAUgSZLTw/Ishg0iIqobhm3gSPII3jvzHqKpKCxYUKBAD+hY37Qeq4Orocma08P0HIYNIiKqC9FUFIPjgxgxRiBBQlgJwy/5YQoThxKHMJQYQpfWhb7WPugBFoFcTAwbRETkedFUFAOnBjBmjqHb3w2f7Jvx5y2+FmTsDIaNYQyMDqC/vZ+BYxG5us4GERHRfAzbwOD4IMbMMfRoPecEjRyf7EOP1oOxzBgGxwdh2KzjtFgYNoiIyNOOJI9gxBhBt7973k2gkiRB9+sYMUZwJHlkiUbofQwbRETkWUIIvHfmPUiQ5pzRmM0v+wEA7515Dy7u6FFTGDaIiMizUiKFaCqKsBIu+OdCCJjChGEbMIWZDxfNSjOiqShSIrWUw/UsbhAlIiLPMm0TFiz4Jf/M7wsT45lxnEifQMyKQQgBScqeUFnhXwG/5IcFC6Zt8mP5ImDYICIiz1JlFQoUmMLMfy9mxnA4eRhxMw5JkqBJGhRJgQ0bo+YoTmVOwS/5sTKwEqrMy+RiYF4jIiLPCkgB6AEdMSsGIBs0Dp49iLgZR0gJoVlpRlAOQpM1BOUgmpVmhJQQTluncTJzEqPpUYdfgTcwbBARkWdJkoT1TeshIJC0kjicPIyz1lk0K81QJGXOxzVIDfBLfrww/gKPwC4Chg0iIvK01cHV6NK68LPkz/IzGsWOwMatOMK+MNYE1/AI7CJh2CAiIk/TZA03tNyAjJ1Bwk7MeT9LWDhtnkaD0oA1DWsQVIIAeAR2MXDnCxEReV671o5OfydSIoW4FQcABOQAZMiwYSNlZ4+4htQQ1jSsyR+VnX4ENigFHRu/2zFsEBGR55m2iaASxLqGdcggkz3yasZgwoQECR3+Dqzwr0Cr2gpV+uTSqEoqkiLJI7AVYtggIiLPyx2BBYBOXyc6fB0whQlb2JAlGaqkQsK5+zhMYUKBwiOwFWJOIyIiz5t9BFaCBJ/kgyZr8Em+gkEDACatSegBHQEpsJTD9RyGDSIi8rzpR2Azdqakx6TtNABgfdP6eRu4UXEMG0REVBdyR2CH08Pzni4RQiCajqJL68Lq4OolGqF3MWwQEZEnCCGQtJOYMqeQtJPnBApN1tDX2oc2XxsiRiQ/czFb2k4jYkTQ5mtDX1sfNFlbiuF7Gne8EBGRqxm2gSPJI3jvzHuIpqKwYEGBAj2gY33TeqwOrs4HBj2go7+9H4PjgxgxRgBkj7eqkgpTmJi0JrP303T0tfVB13SnXpanSKKOKpXE43GEw2HEYjGEQiGnh0NERBWKpqL54CAh27U1FxxiVgwCAl1aF/pa+6AHPgkO5QQUyqrkGsqwQURErhRNRTFwagBj5hi6/d3wyb5z7pOxMxhOD6PN14b+9v4ZgQPILr2kRAqmbUKVVQSkADeDzqGSayj3bBARkesYtoHB8UGMmWPo0XoKBg0A8Mk+9Gg9GMuMYXB88JymapIkISgHsUxdhqAcZNCoEoYNIiJynSPJIxgxRtDt7543IEiSBN2vs6magxg2iIjIVYQQeO/Me9nCXHPMaMzml/0A2FTNKQwbRETkKimRQjQVzTdLK9X0pmq0tBg2iIjIVUzbhAVrRsO0UqiSCgtWtqkaLSmGDSIicpVcUzVTlBca2FTNOQwbRETkKrObqpWKTdWcw7BBRESuwqZq7sOwQURErsOmau7imrCxc+dOXHHFFVi2bBk6OjrQ19eH999/3+lhERGRA9hUzV1cEzZee+01bN++HW+++Sb27NmDTCaDL37xi0gkEk4PjYiIHJBrqqZrOkbSI4ikIpjITCBuxjGRmUAkFcFIegS6pqO/o59N1Rzk2t4oo6Oj6OjowGuvvYZrrrmmpMewNwoRkfewqdrSqOQa6trzP7FYdhdyS0vLnPcxDAOG8Ukd/Hg8XvVxERHR0tJkDesa16G3oZdN1WqUa5ZRprNtG3fddRc++9nPore3d8777dy5E+FwOH/TdU6hERF5FZuq1S5XLqPcfvvteOmll/Cf//mf6OrqmvN+hWY2dF3nMgoREVGZ6moZ5Y477sCLL76I119/vWjQAABN06BpXKcjIiJykmvChhACd955J3bv3o29e/di1apVTg+JiIiISuCasLF9+3Y899xzeOGFF7Bs2TKcOHECABAOhxEMBh0eHREREc3FNXs25tro8+STT+Kmm24q6Tl49JWIiGhh6mLPhksyEREREc3iyqOvRERE5B4MG0RERFRVDBtERERUVa7Zs1GPhBAsvUtERK7HsFGD2FSIiIi8hGGjxkRTUQyOD2LEGIEECWElDL/khylMHEocwlBiCF1aF/pa+6AH2OuFiIhqH8NGDRlODePZU89iLD0GXdMRVIKQ8MmySYuvBRk7g2FjGAOjA+hv72fgICKimscNojXAsA3sm9qHrw9/HW/G3sTx9HG8c+Yd/Hfiv3EycxKmMPP39ck+9Gg9GMuMYXB8EIZtFHlmIiIi53Fmw2G5ZZOhxBCGjWEsV5ZDlVTYsDGaHsWp9CmE1BDWBNcgrIYBZKup6n4dI8YIjiSPYF3jOodfBRER0dw4s+GgaCqKgVMDGE4NQwiBJqUJjUojNFlDUA6iWW1GWAkjbsZx8OxBxMxY/rF+2Q8AeO/Me6yuSkRENY1hwyGGbWBwfBBj5hhWaitxxj4DTTr3hIksyViuLsdZ6ywOJw/PWFJpVpoRTUWREqmlHDoREVFZGDYcciR5BCPGCLr93bBhQwgBuch/jpASQtyMY9wcz39PlVRYsGDa5pyPIyIichrDhgOEEHjvzHuQIMEn+6BICiRJgg17zscokgIAOJE+AYHssokpTChQoMrcekNERLWLYcMBKZFCNBVFWMlu+FSgIKyEYYjiJ0sCcgAxM5ZfSpm0JqEHdASkQNXHTEREtFAMGw4wbRMWLKhSdkZCkiSs8K+AEAKWsOZ8nAwZAgK2sJG20wCA9U3rWcKciIhqGsOGA1RZhQJlxmbPVl8rQmoIcSs+5+kSGzakX/xfNB1Fl9aF1cHVSzVsIiKiBWHYcEBACkAP6IhZnxxlVSUVa4Jr0KA0YNKaLDjDkbJTaFKaEE1H0eZrQ19bH3ukEBFRzWPYcIAkSVjftB4CAhk7k/9+WA2jt6E3P8MxaU0iaSdh2AYSVgIJOwEJEnRNR39HP3SNpcqJqL4JIZC0k5gyp5C0k6w7VKN4jMEhq4Or0aV1YdgYRo/Wk993EVbD2NC0AROZCRxPH0fMym4InbKmoGs6bjvvNqxtXMsZDSKqa+yO7S6SqKMYGI/HEQ6HEYvFEAqFnB5OtoLo6ADGMmPQ/Xq+KmhOLrEPG8No97dja8dWdAe6HRotEVFtKNQdW5VUmMJEzIpBQLA7dhVUcg1l2HDY9H80QLYqaO4fzaQ1CQDZfzRtfVw2IaK6l2vzMGaOodvfDZ/sO+c+GTuD4fQw2nxt7I69iBg2SlSLYQPgdCARUSkM28Djxx9H1IjOWH4uRAiBiBHJLz/zPbRylVxDuWejBmiyhnWN69Db0IuUSMG0TaiyioAUYA0NIqJfmN7mYb73RnbHri08jVJDJElCUA5imboMQTnIoEFE9Auz2zyUgt2xawfDBhER1bzZbR5Kxe7YtYFhg4iIat7sNg+lYnfs2sCwQURENa9Qm4dSsDt2bWDYICKimleozUMp2B27NjBsEBFRzZurzUMx7I5dOxg2iIjIFfJtHtLD854uEUKwO3YNYdggIiJX0GQNfa19aPO1IWJE8jMXs6XtNCJGhN2xawh3zBARkWvoAR397f3ztnnQNZ1tHmoIwwYREbmKHsiWIJ/e5iEpklCgYG3jWrZ5qEGuChuvv/46HnjgAezbtw/Hjx/H7t270dfX5/SwiIhoibHNg7u4as9GIpHApZdeikcffdTpoRARURUJIZC0k5gyp5C0k3NuCGWbB3dw1czG9ddfj+uvv97pYRARUZWwC7Y3uSpslMswDBiGkf86Ho87OBoiIiommormN35KkBBWwvBLfpjCxKHEIQwlhtCldaGvtQ96gBs/3cRVyyjl2rlzJ8LhcP6m6/zlJCKqRdFUFAOnBhA1otD9OlYFVqHF14KQGkKLrwWrAqug+3VEjSgGRgcQTUWdHjKVwdNh45577kEsFsvfolH+chIR1RrDNjA4Pogxcww9Ws+cLeR9sg89Wg/GMmMYHB+EYRsF70e1x9PLKJqmQdPqd21PCMFd2kRU844kj2DEGEG3v3ve9yhJkqD7dYwYIziSPIJ1jeuWaJRUCU+HjXrFDVZE5BZCCLx35j1IkOac0ZjNL/sBAO+deQ+9Db38EOUCrgobZ86cwdGjR/Nff/jhh9i/fz9aWlrQ3d3t4MhqBzdYEZGbpEQK0VQUYSVc1uOalWZEU1GkRApBKVil0dFicVXYeOedd7B58+b81zt27AAAbNu2DU899ZRDo6oduQ1WY+YYuv3d53xKaPG1IGNnMGwMY2B0AP3t/QwcROQo0zZhwYJf8pf1OFVSkRRJmLbp8d2H3uCqsLFp06Z5O/3Vq9kbrOaaVsxtsIoYEQyOD+K2827jkgoROUaVVShQYAqzrMeZwoQCBarsqstY3WIe9IhKNlgRETklIAWgB3TErFhZj5u0JqEHdASkQJVGRouJYcMDKt1gxdkiInKKJElY37QeAgIZO1PSY3Kt5dc3refmUJdg2PCAxdhgRUTklNXB1ejSujCcHp73w48QAtF0FF1aF1YHVy/RCKlSDBsekNtgpUrlrV2qkgoLVnaDFRGRQzRZQ19rH9p8bYgYkfzMxWxpO42IEUGbrw19bX3cb+Yi3FnjAdxgRURupwd09Lf354/uA9nZV1VSYQoTk9Zk9n6ajr62PugaT9K5Ca8yHpDbYHUocQgtvpaSHzdpTWJt41pusCKimqAHdNx23m0zihImRRIKFKxtXMuihC7GsOEBuQ1WQ4khZOxMSZtEy91gxdLnRLQUNFnDusZ16G3o5XuOhzBseER+g5UxXLTOBvDJBitd0+fdYMXS50TkBEmSspVBubPQExg2PCK3wWpgdAARIwLdr+ePt06XttOIpqMlbbBi6XMiIloMkqijIgvxeBzhcBixWAyhUMjp4VTF9IAAFN5g1aV1zbvBar7S5wCypc/Tw2jztbH0ORGRx1VyDWXY8KBiSx+XNV6GlYGVkCHPuQ5q2AYeP/44oka0pCWZiBGBruksfU5E5GGVXEO5jOJBhTZYWbAwnBrG/sR+/L/x/1d070Ulpc/XNa4DwA2lRET0CYYND8ttsIqmS9970aV1VVT6/NOBT+N/Uv/DDaVERJRX9jLKtm3bcMstt+Caa66p1piqpl6WUaYrd+/Fb7b+Jr47+l3IkMuq2TGRmUDciqPD34ET6RP5UJPbLxKzYhAQ3FBKRORSlVxDyz5UFIvFsGXLFnz605/GN77xDRw7dqzcp6AlMrvt/FwzFbm282OZMQyODyIjMjNKnwsIZEQGhm0gIzIQODefJu0k/if5P4imotD9OlYFVqHF14KQGkKLrwWrAqug+3VEjSgGRgcQTUWr9rqJyL2EEEjaSUyZU0jaSTaK9IgFbRAdHR3FM888g6effhqHDh3Cli1bcMstt+CGG26Az1fa1LsT6m1m40DiAJ49+Sx0v15yoa+fGz+HgECb2oaQGsK4OY4T6ROImdmZCQkSwmoYK/wr0Kq25mcu/iv2XzhjncEXln+h4JHbHG4oJaJCWNOn9jl6GuXdd9/Fk08+iSeeeAJNTU3YunUr/vAP/xCf/vSnK3naqqinsCGEwDOnnsGhxCGsCqwq+XEfJD+AIQwICJyxziBuxiFBgiZrkCHDhg3Dzv55SA1hTXANUnYKP576Mbq17mxFUhTfCJq20xhJj2Br59b8hlIiql+FavpwCbb2LOkyynTHjx/Hnj17sGfPHiiKgl/5lV/BgQMHcPHFF+Ohhx6q5KmpQgttO79cXY60ncZHqY8wmZlEWAmjWW1GUA5CkzUE5SCa1WaElTDiZhwHEgfwQfIDAMD/0v7XvEEDmLmhlFOkRPUtt68sanAJ1svKDhuZTAY//OEP8Wu/9mtYuXIlvv/97+Ouu+7Cxx9/jKeffho/+tGP8L3vfQ9/+Zd/WY3xUokW2nZegoSxzBhkyJAlec7wIEsylqvLccY8g4+Mj9CsNKNVbS3572lWmhFNRZESqbLGR0TesdB9ZYZtLPFIqVJlH30977zzYNs2+vv78dZbb+Gyyy475z6bN29Gc3PzIgyPFmqhbefHMmMwhIHLmy7H0dRRTFqTCCkhKJJyzn0tYUFIAraw0eZrKyvYqJKKpEjCtE32PiCqU4tR04fcoeyw8dBDD+HGG29EIDB3W/Lm5mZ8+OGHFQ2MKrOQtvO5Bm0NcgM6/B0IyAEcTh7O7tuQJGjStH0bwoAQAs1qM9J2GlPWVH4DaSlMYUKBAlVmqReieiSEqKimT29DLwsFukjZnym/8pWvFA0aVBtybecFBDJ2pqTHJO0kklYSXVpX/tTJhqYNuLTxUrSr7RAQMIUJAYF2tR2XNl6Ky5suR4evAxPmRFmzKJPWJPSAjoDE3yWierTQfWVcgnUnfqz0sHLbzg8bw2hQGmbsvVAlFR3+DrT72mHBgiUsKJICBUr++Tr9nRhPjsOwDPjU0o7YAsieXOEnE6K6lNtX5pfmPipfCJdg3Yn/qTws13a+zdeGiBHJX+RnS9tpRIwI2v3tcx6TlSQJqqRCkzWokjojJCxTlqFJbsLHmY/nPV2SW6rp0rqwOrh64S+OiFxtofvKuATrTgwbHqcH9Gz7d03HSHoEkVQkW1rcjGMiM4FIKoKR9Ah0TcfWjq34343/GzErVtbfMWVPYWN4Izp9nSWFmjZfG/ra+ligh6iO5faVlft+wyVYd2I0rAN6IFutc3p1vqRIQoGCtY1rZ1TnW2+tx1BiCBk7U3LVUQDY3LwZzUpzvjAPkF1bzRXmmbQms2PRdPS19UHXWJiHqJ7l9pUt5P2GS7Duw7BRJwq1nS/U+r3cfR7RdBS6pufDSqmhhoiokvcbcpeKy5W7ST2VK69ENJWt1DeWGYPu1wv2OknbaUTTUbT52tDf0X/OTIUQomioISICFuf9hpaGo71R3IRho3TTexUAhZdEurQuLokQUcX4fuMODBslYtgoT6ldGDmLQUSVYtfX2sewUSKGjYWZK0zwzYGIFhs/vNSuSq6h3CBK85IkCUEpOOOgdKGW0H7JD1OYOJQ4hKHEEFtCE1HZCr3fkPsxbFDZci2hx8wxdPu7zzmy1uJrQcbOYNgYxsDoQLbOBwMHEVHdcl12fPTRR/GpT30KgUAAV111Fd566y2nh1RX2BKaiIjK5aqw8d3vfhc7duzAvffei3fffReXXnopvvSlL+HUqVNOD61uVNISmoiI6pOrwsa3vvUt3Hrrrbj55ptx8cUX49vf/jYaGhrw93//904PrS5U2hK6jvYiExHRNK4JG+l0Gvv27cOWLVvy35NlGVu2bMGPf/zjgo8xDAPxeHzGjRaOLaGJiGghXBM2xsbGYFkWOjs7Z3y/s7MTJ06cKPiYnTt3IhwO52+6zk2Klci1hFal8vYVq5IKC1a2JTQREdUd14SNhbjnnnsQi8Xyt2g06vSQaoYQAkk7iSlzCkk7WdISB1tCExHRQrjm3b+trQ2KouDkyZMzvn/y5EmsWLGi4GM0TYOmsajUdJUU4sq1hD6UOIQWX0vJf+ekNYm1jWvZEpqIqE65ZmbD7/djw4YNeOWVV/Lfs20br7zyCjZu3OjgyNwjmori8eOP49mTz+JQ4hBkyAhKQciQcShxCM+efBaPH38c0VThGaBcS2gBgYydKenvZEtoIiJyzcwGAOzYsQPbtm3D5ZdfjiuvvBIPP/wwEokEbr75ZqeHVvMWqxAXW0ITEVG5XDOzAQC//du/jQcffBBf+9rXcNlll2H//v3413/913M2jdJMi1mIS5M19LX2oc3XhogRyc9czJa204gYEbT52tDX1sceKUREdYyN2OrAgcQBPHvyWeh+vaT6GGk7jZH0CLZ2bsW6xnUF78OW0ERE9YWN2GhOlRbi6m3oLbhUogd03HbebTM2myZFEgoUrG1cy66vRESUx7DhcYtRiCsoBQveR5M1rGtch96GXraEJiKiOTFseFyuEJdf8pf1OFVSkRTJbCGueXb2sCU0EREVw7DhcV4pxCWE4OwJEZFL1caVhKrG7YW4KilCRkREtYFhw+NyhbiGEkPI2JmST6MAzhfimn7iRYKEsBKGX/LDFCYOJQ5hKDGUPfHS2lewJggREdUGho064MZCXItVhIyIiJzHLX11wG2FuBazCBkRETmPMxt1Qg/o6G/vn7cQl67pjhfiOpI8ghFjBN3+7nmXcSRJgu7XMWKM4EjyyJxFyIiIyDkMG3XEDYW4qlWEjIiInMOwUWdqvRBXNYuQLQSP3BIRVY5ho07VaiGupShCVgoeuSUiWjwMG1RTaqEIGY/cEhEtrhr7XEv1LleELGbFynrcpDUJPaBXXIQsd+Q2akSh+3WsCqxCi68FITWEFl8LVgVWQffriBpRDIwOIJqKVvT3ERHVA4YNqim5ImQCAhk7U9JjFqsIGY/cEhFVB8MG1Zx8EbL0MIQQRe+bK0LWpXVVXISskiO3REQ0N4YNqjlOFCGr9MjtfKGIiKiecYMo1aSlLkJWa0duiYi8hGGDatZSFiGrlSO3RERexLBBNW2pipDVwpFbIiKv4mcxcgVJkhCUg1imLkNQDi56FU+nj9wS1RMhBJJ2ElPmFJJ2knue6gA/jhHhkyO3Q4khZOxMSZtEF+vILVG9YGXe+sWwQfQL+SO3xjB6tJ6iASJ35FbX9IqP3BLVA1bmrW9cRiH6BSeO3BLVA1bmJYYNomlyR251TcdIegSRVAQTmQnEzTgmMhOIpCIYSY9A13T0d/RXfOSWyOtYmZcALqMQnWMpj9wSeV0llXnXNa5bolFStTFsEBWwVEduibys0sq8vQ29/PfmEQwbREVIkpStDMoFR6KysTIv5fAtlIiIqiJXmVeVyvtcq0oqLFjZyrzkCQwbRERUFazMSzkMG0REVBWszEs5DBtERFQVucq8AgIZO1PSY1iZ15sYNoiWAHtBUL3KV+ZND8/7e5+rzNuldbEyr8e4ZkHs61//Ov75n/8Z+/fvh9/vx+TkpNNDIpoXe0FQvctV5h0YHUDEiED36/njrdOl7TSi6Sgr83qUa8JGOp3GjTfeiI0bN+I73/mO08Mhmhd7QRBl5Srz5v49ANnjraqkwhQmJq3J7P00HX1tfazM60GScNl87lNPPYW77rprQTMb8Xgc4XAYsVgMoVBo8QdHdUkIcU7hrxFjBAOnBjBmjqHb312woFHGzmA4PYw2X1u2RDoDB3kcZ/rcrZJrqGtmNhbCMAwYxif19ePxuIOjIa+Z643zPO08fJT6CGets7gwcOGcm9xyvSAiRgSD44O47bzb+EZLnsbKvPXL02Fj586duO+++5weBnlQsSWSt+Jv4efGz9Hh60CHrwNhde7qiewFQfWIlXnrj6P/qe+++25IklT0dvjw4QU//z333INYLJa/RaNsW0yVK9Yue7m6HDJkNMgNSFgJHDx7EDGzeI2B6b0gXLaqSURUEkdnNr761a/ipptuKnqfnp6eBT+/pmnQNE5L0+KZ3S579tSvBQsxK4agHERQDuK0eRqHk4exoWlD0ZLN7AVBRF7maNhob29He3u7k0MgKst87bItYUEIAUVSAAAhJYS4Gce4OY5OX+ecz6tKKpIime0FwallIvIY1+zZGB4exsTEBIaHh2FZFvbv3w8AuPDCC9HU1OTs4KgulNIuW5EUSJIEG3b+awA4kT6BDl8HJBTeBMdeEETkZa55Z/va176Gp59+Ov/1+vXrAQCvvvoqNm3a5NCoqJ6U0i5bgYKwEsaoOYogssshATmAmBmDKUz4pMIhZdKaxNrGtewFQUSe5JoJ26eeegpCiHNuDBq0VEpply1JElb4V0AIAUtYAAAZMgQEbGEXfAx7QRCR17kmbBA5rdR22a2+VoTUEOJWHEII2LAhQYIsnfvPjb0giKgeMGwQlajUdtmqpGJNcA0alAZMWpM4a51FWA2fMyOSttOIGBH2giAiz2PYICpROe2yw2oYvQ29aFKakBTZLq+nM6cRN+OYyEwgkopgJD0CXdPR39HPXhBE5Gmu2SBKVAvy7bKN4YJ1NqYLKSG0+lrRHejGp7RP4Xj6OJIiCQUK1jauXVAviEJ9WLjPg4hqHcMGURnKbZfd6e9Ef0c/uvxdFYUENrCiWsHASwvhuq6vlWDXV1os03ujAIXbZXdpXYvSLrtQH5bc3xWzYhAQnm9Vzwuc8xh4qZJrKMMG0QItxZtvrg9Lvbaq5wWuNjDwEsCwUTKGDaqGan3qNmwDjx9/HFEjOu/+ECFEdllH0z3Tqp4XOGfM/n0eTY/i+VPP123gpU9Ucg3lng2iClWrXfZ8fVhmj8FLrernm9Fp8bVkL3DGMAZGB3iBWwSFZpEggGPpY5Agobehd84y/T7Zhx6tBxEjgsHxQc8EXlo8PPpKVINK6cMyWymt6oUQSNpJTJlTSNrJmmxpP7uz7nwXuLHMGAbHB2HYxhKP1DuiqSgeP/44nj35LA4lDkGGjKAURNyMY9gYxqn0KbybeBcxc+4aM7MDL9F0nNkgqkGl9GEpZK5W9W7a+1DPMzpOmGsWSQiBD+0P0aQ0ISSHEDNjOHj2IHobehFWC/9eTg+8vQ293MRLeQwbRDUo14fFL517rLaYQq3qC+198Et+mMLEocQhDCWGambvQ6UzOrzAlWf2LNL0n50FCzErBk3SIEsylqvLcdo8jcPJw9jQtGHOHkFzBV6qbwwbRFVQ6abRUvuwzDa7Vb3b9j4s9owOFVdsFskSFoQQUCQl/72QEkLcjGPcHEenr7PgcxYKvEQMG0SLaLGWK3J9WA4lDqHF11Ly3z+9VX2xT63T1dLmvsWc0aHi5ptFUiQFkiTBhj3jewBwIn0CHb4OSDj3d2p24CUCGDaIFs1iLlfk+rAMJYaQsTMlLSnMblV/5Kz79j4s1oxOqeq5WFhuFqlRaYRhG5AlGaqk5gOEAgVhJYxRcxRBfDJbFJADiJkxmMKETzr393J64CXKYdggWgTVWK4opw9LrlW9rulYHVzt2r0PizGjUwo3bZitBsM28NMzP8X7yfeRtJLwyb5sQFbDWOFfgVa1FaqkYoV/BU5lTsESVn5WQ4YMEyZsYWP2xMbswEuUwwlHogpV66hmrg9Lm68NESOSfyOfrVCr+sXY++CEcjrr5pR7gZvrmKcMGYcSh/DsyWezxdRS0YpeS63Kvf4fjP4Ak5lJAICK7IzGaHoU/33mv7HvzD7EzBhafa0IqSHErXj+mLQNGxIkyNLMy0cu8HZpXVgdXL3UL4tqHMMGUYUqOao5Hz2gZ2dBNB0j6RFEUhFMZCbmbVWf2/sw14mBuaiSCgtWdu+DQ/IzOunheeuAlHuBy81ARY0odL+OVYFVaPG1IKSG0OJrwarAKuh+HVEjioHRAc8Fjumvf6V/Jbq0LkDKBtugHESz2oywEkbcjOPg2YNIWAmsCa5Bg9KASWsSlrCQslMIq+EZv1uFAi/RdAwbRBWoVvGt6fRAtgT51s6tWNu4FjZsJEUSNmysbVyLrZ1bcdt5t81o+LbUex8WUyUzOsXUe7Gw2a/fr/ixwr8CQghYwsrfL3fM9ax1FoeTh9GoNKK3oRchNYRJaxIJO4EGuQFT5lTRwEs0HfdsEFVgqY5qarKGdY3r0NvQW9KGxqXa+1AtuRmd+Trr6ppecmfdei8WVuj155dJzDialeYZP5fZx1x/qfGXcODsAQgIhJQQkiIJBQrWNq6ti30uVBmGDaIKLPVRzVL7sCzGaRan5WZ0pm/kXOgFzq0bZhfLXK9flVSsCa7BwbMHMWlNIqSE8htBpx9zbVaaMZIewQXBC/B/2/8v2v3tdXmChxaOYYOoArW8XFHJaZZaUe6MzlzqvVhYsdcfVsPobejF4eRhxM04JEnKVg2FDAkSRowRNMgNWBlYOXMWiYvwVAb+uhBVILdcEbPmblBVyKQ1CT2gV3W5olp7H5wgSRKCchDL1GUIysGyP0m7ecPsYpjv9YfVMDY0bcCljZeiXW2HgIApzHy9mBvbbzxnXxBROTizQVSBWl+uqMbeBzeq5RmopVDK61clFR3+DrT72mHBgiUsxMwYJEi4pOmSmgyh5B7u/hdEVANqfbliMfc+uJXbN8wWUk7103JevyRJUKFClVScsc/U7Osnd2HYIKpQbrliYHQAESMC3a/nNxdOl7bTiKajjixXLNbeB7eq9Rmociyk+qmXXj+5E8MG0SJwy3JFqadZvKjWZ6BKUUn/HS+8fnIvSZRSVcgj4vE4wuEwYrEYQqGQ08MhD6r3nhu1LprKVgYdy4yVNANVS0Wq5uu/AyDbfyc9nB17gf47bn795LxKrqEMG0RVUM/dRGvd9NkBoPAMVJfWVVMbZg3byPZrMaIlzUpEjAh0LbtXZ3a4dePrp9rAsFEihg0iAtw3A3UgcQDPnnwWul8veb/FSHoEWzu3Fqx+6rbXT7Whkmso92wQUd1x04bZalQ/ddPrJ29g2CCiuuWGDbPVrH7qhtdP3sBfMSJyNSEEknYSU+YUknaypE66blLv1U/JG1wxs/HRRx/hr/7qr/Dv//7vOHHiBM4//3xs3boVf/7nfw6/v7wGWETkDfWy76Deq5+SN7jit/Dw4cOwbRu7du3ChRdeiIMHD+LWW29FIpHAgw8+6PTwiGiJVVJvwm28WP2U6o8rwsZ1112H6667Lv91T08P3n//fTz22GMMG0R1Zr56Ey2+lmy9CWMYA6MDBetNuAmrf5IXuHbPRiwWQ0tL8ZRvGAbi8fiMGxG5l2EbGBwfxJg5hh6tZ84Lr0/2oUfrwVhmDIPjgzBsY4lHurjy1T/Tw/PuSclV/+zSulj9k2qGK8PG0aNH8cgjj+D3f//3i95v586dCIfD+Zuuu/fTDREBR5JHMGKMoNvfPe8ndkmSoPt1jBgjOJI8skQj/MRiblzN9d9p87UhYkTyMxezpe00IkbEkf47RMU4WtTr7rvvxje/+c2i9/nZz36GNWvW5L8+duwYPve5z2HTpk144oknij7WMAwYxiefaOLxOHRdZ1EvIhcSQuCZU8/gUOIQVgVWlfy4SCqCtY1r8ZWOryzJkkI1N66y+ic5ybUVREdHRzE+Pl70Pj09PfkTJx9//DE2bdqEz3zmM3jqqacgy+VNzLCCKJF7Je0kvhX9FmTIZW2UnMhMwIaNHfoOBOXC9SbKNVc5+kIbV3NhIGbFICAq3rhaL6dwqPa4toJoe3s72tvbS7rvsWPHsHnzZmzYsAFPPvlk2UGDiNwtV2/CL5V33F2BgqRI4nTmNOBDRVUyi13oz/efj7en3sZp83RVN66y+ie5kStOoxw7dgybNm3CypUr8eCDD2J0dDT/ZytWrHBwZES0VMqtN2EKE+OZcUSSEUzZU9h1fBf8kn/BMwDFjtseOHMAL6RegIDAVU1XzbtxNWJEMDg+WLBRWqlY/ZPcxBVhY8+ePTh69CiOHj2Krq6uGX/mtWqBRFRYOfUmYmYMh5OHETfjSNgJdPo60Sg3whLWgupwzHfcNiMyEBCAAIaSQ+iVehFWC5cXn71xtVCjNCKvcUUmvummmyCEKHgjovqQqzchIJCxM3PeL2bGcPDsQcTNOBrlRjTKjegJ9iCshtHia8GqwCrofh1RI4qB0QFEU9Gif+98x22FEDiRPgFFUtDqa8VZ6ywOJw8XnYGZ3iht+vuY10uvU/1yxcwGEREwrd6EMYwereecPQqmMHE4eRhnrbNoVpoxaU0ipIbQqrbOuF85yxnzHbe1YCFmxaBJ2ceHlBDiZhzj5jg6fZ1zvpbpjdJkIXPTJ3kawwaRy811MsKLcvUmBkYHEDEi0P16fpYAAMYz4/kZjUlrEg1KA9Y0rCnYxKyU5YxS2rtbwoIQAoqkAED+f0+kT6DD1wEJhf9bqJKKpEjio+RH+NHkj+qi9DrVL4YNIpeq1yOQekBHf3v/OfUmFCiIJCNI2AkAQEgNYU3DmqKt2acvZ/Q29J4T0kpp765ICiRJgg07/72AHEDMjMEUJnxS4ZBiChNJK4nd47sxaU7WRel1ql8MG0QuVE+NyArRAzpuO++2GWErKZKYsqfQ6etET7AHrWprSW3Zpy9nBKWZdThKOW6rQEFYCWPUHEUQ2cfLkGHChC1szDGxgXFzHFPWFCxYBZeEchbzBAuRUxg2iFym3hqRzWV2vYnTmdPYdXwXGuXGOU+CFJJbzjBt85wt86Uct5UkCSv8K3AqcwqWsKBICmzYkCBBlgrvwU/baZyxzkCCtKDS6zzBQm7jitMoRJRVr43IipEkCUE5iOW+5fBLfljCKuvxpjChQIEqn/vZK3fcNmbFij5Hq68VITWEuBXP7qGxUwir4YIzK0IIDBvDsGEjpIZK6uIKzH2ChcgNGDaIXMRNjciWWqnBYLZJaxJ6QEdACpzzZ6Uet1UlFWuCa9CgNGDCmoAtbKzwrzhnc2iuUdpy33K0qq1Yriwva6zTl3yI3IRhg8glSjkZMVs9fRouNRhMl+ueur5p/ZzhrdT27mE1jLXBtb8YDDBlTWEiM4G4GcdEZgKRVAQj6RHomo4vt30ZQSVY0p6S6VRJhQUru+RD5CLcs0HkEqWcjCik2AZIr5mvDsd0QghE01Homo7VwdVz3m++47Y5aTuNcWscVyy7AlcuuxLH0sfyG1cVKFjbuDZ/SsiGXVbp9ZxiSz5EtYy/sUQusdBGZMU2QHpNOcEgmo6izdeGvra+eU93zHXcdnZ7d13T8+3di9U/EUKUXHp9uklrEmsb1xZc8iGqZQwbRC5RbiOynHr7NLyQYFDq8xY6bjt71iIXXIo1Ssst+QwlhpCxMyUti5Wy5ENUq+rj3YfIA8ppRDZdPX4aLjcYlGox27tXY8mHqFYxbBC5BD8Nl2cxg8Fsi9HevVpLPkS1iGGDyEX4abh8ixEMKlFs70a1lnyIao0kvH4ebpp4PI5wOIxYLIZQKOT0cIgWJJrKtkYfy4yV9Gm4v6OfFykHlNO7pl773JC7VHINZdggcqHpvVGAwp+Gu7QufhquomIzFoV61+T++8SsGAREwd419dTBl9yHYaNEDBvkJfw07Iz5fu4NcgN+OPrDOXvXAMj2rkkPZ2eePNq7hryHYaNEDBvkRfw0vHTmm7EwhYmxzBiWKcuwtmHtvHtqIkYEuqazkyu5QiXXUG4QJXI5pzdA1otSuu0eM47hZOYkMiKDuBUv2n2WnVypnvDtiYhoHqV02xVCYCwzhia5CYZt4HDy8LwF2Oqpdw3VN4YNIqJ5lNJt14KFmBWDJmsIKSHEzTjGzfF5n5udXKkeMGwQERVRarddS1gQQkCGDEVSAAAn0icgUHzGgp1cqR4wbBARFVFqt11FUiBJEmzYAICAHEDMjM27lFJvvWuoPjFsEBEVkeu2q0rFw4ACBWElDEMYAAAZMgQEbGEXfdykNQk9oNdV7xqqPwwbRERFlNptV5IkrPCvgBAClrBgw4YECbI099tsPfeuofrCsEFEVESu227Mis1731ZfK0JqCHErjqSVRFgNzzkjkutd06V11XXvGqoPDBtEREXkuu0KCGTsTNH7qpKKNcE10GQNZ8VZtPpaIeHcGYu0nUbEiLCTK9UN7kgiIppHOd12Q0oIy9Xl8Ek+JK0kIqkIO7lS3WPYICKahyZr6Gvtw8DoQLbE+DzddrsD3fg/bf8HZ+2z+R4qSZGEAgVrG9eydw3VHYYNIqIS6AEd/e3983bbnT1j0dvQy941VPcYNoiISqQHsk3Tpnd9nW/Ggr1riBg2iIjKoska1jWu44wFURkYNoiIFoAzFkSlc80/k9/4jd9Ad3c3AoEAzjvvPHzlK1/Bxx9/7PSwiIiIaB6uCRubN2/G9773Pbz//vv44Q9/iA8++AC/9Vu/5fSwiIiIaB6SEKJ4S8Ia9U//9E/o6+uDYRjw+ebuxDhdPB5HOBxGLBZDKBSq8giJiIi8o5JrqCv3bExMTOAf/uEfcPXVVxcNGoZhwDCM/NexWLbccDwer/oYiYiIvCR37VzQHIVwkT/5kz8RDQ0NAoD4zGc+I8bGxore/9577xUAeOONN9544423Rbp98MEHZV+/HV1Gufvuu/HNb36z6H1+9rOfYc2aNQCAsbExTExM4Oc//znuu+8+hMNhvPjii3MeN5s9szE5OYmVK1dieHgY4XB48V4IzSkej0PXdUSjUS5dLRH+zJcef+ZLjz/zpReLxdDd3Y3Tp0+jubm5rMc6GjZGR0cxPj5e9D49PT3w+88tCzwyMgJd1/HGG29g48aNJf193LOx9PgzX3r8mS89/syXHn/mS8+1ezba29vR3t6+oMfatg0AM2YuiIiIqPa4YoPoT37yE7z99tv45V/+ZSxfvhwffPAB/uIv/gIXXHBBybMaRERE5AxX1NloaGjAP/7jP+Laa6/FRRddhFtuuQWXXHIJXnvtNWha6V0TNU3DvffeW9ZjqDL8mS89/syXHn/mS48/86VXyc/ctXU2iIiIyB1cMbNBRERE7sWwQURERFXFsEFERERVxbBBREREVVWXYeOjjz7CLbfcglWrViEYDOKCCy7Avffei3Q67fTQPO3rX/86rr76ajQ0NJRdfY5K8+ijj+JTn/oUAoEArrrqKrz11ltOD8nTXn/9dfz6r/86zj//fEiShMHBQaeH5Gk7d+7EFVdcgWXLlqGjowN9fX14//33nR6Wpz322GO45JJLEAqFEAqFsHHjRrz00ktlP09dho3Dhw/Dtm3s2rULQ0NDeOihh/Dtb38bf/Znf+b00DwtnU7jxhtvxO233+70UDzpu9/9Lnbs2IF7770X7777Li699FJ86UtfwqlTp5wemmclEglceumlePTRR50eSl147bXXsH37drz55pvYs2cPMpkMvvjFLyKRSDg9NM/q6urC/fffj3379uGdd97B5z//edxwww0YGhoq63l49PUXHnjgATz22GOIRCJOD8XznnrqKdx1112YnJx0eiiectVVV+GKK67A3/7t3wLIVtnVdR133nkn7r77bodH532SJGH37t3o6+tzeih1Y3R0FB0dHXjttddwzTXXOD2cutHS0oIHHngAt9xyS8mPqcuZjUJisRhaWlqcHgbRgqTTaezbtw9btmzJf0+WZWzZsgU//vGPHRwZUfXEYjEA4Hv3ErEsC88//zwSiUTZ1btdUa682o4ePYpHHnkEDz74oNNDIVqQsbExWJaFzs7OGd/v7OzE4cOHHRoVUfXYto277roLn/3sZ9Hb2+v0cDztwIED2LhxI1KpFJqamrB7925cfPHFZT2Hp2Y27r77bkiSVPQ2+4332LFjuO6663DjjTfi1ltvdWjk7rWQnzkRUaW2b9+OgwcP4vnnn3d6KJ530UUXYf/+/fjJT36C22+/Hdu2bcOhQ4fKeg5PzWx89atfxU033VT0Pj09Pfn//+OPP8bmzZtx9dVX4/HHH6/y6Lyp3J85VUdbWxsURcHJkydnfP/kyZNYsWKFQ6Miqo477rgDL774Il5//XV0dXU5PRzP8/v9uPDCCwEAGzZswNtvv42/+Zu/wa5du0p+Dk+FjXJa1h87dgybN2/Ghg0b8OSTT0KWPTXJs2TK+ZlT9fj9fmzYsAGvvPJKfoOibdt45ZVXcMcddzg7OKJFIoTAnXfeid27d2Pv3r1YtWqV00OqS7ZtwzCMsh7jqbBRqmPHjmHTpk1YuXIlHnzwQYyOjub/jJ8Cq2d4eBgTExMYHh6GZVnYv38/AODCCy9EU1OTs4PzgB07dmDbtm24/PLLceWVV+Lhhx9GIpHAzTff7PTQPOvMmTM4evRo/usPP/wQ+/fvR0tLC7q7ux0cmTdt374dzz33HF544QUsW7YMJ06cAACEw2EEg0GHR+dN99xzD66//np0d3djamoKzz33HPbu3YuXX365vCcSdejJJ58UAAreqHq2bdtW8Gf+6quvOj00z3jkkUdEd3e38Pv94sorrxRvvvmm00PytFdffbXg7/S2bducHponzfW+/eSTTzo9NM/6vd/7PbFy5Urh9/tFe3u7uPbaa8W//du/lf08rLNBREREVcWNCkRERFRVDBtERERUVQwbREREVFUMG0RERFRVDBtERERUVQwbREREVFUMG0RERFRVDBtERERUVQwbREREVFUMG0RERFRVDBtERERUVQwbROSo0dFRrFixAt/4xjfy33vjjTfg9/vxyiuvODgyIlosbMRGRI77l3/5F/T19eGNN97ARRddhMsuuww33HADvvWtbzk9NCJaBAwbRFQTtm/fjh/96Ee4/PLLceDAAbz99tvQNM3pYRHRImDYIKKakEwm0dvbi2g0in379mHdunVOD4mIFgn3bBBRTfjggw/w8ccfw7ZtfPTRR04Ph4gWEWc2iMhx6XQaV155JS677DJcdNFFePjhh3HgwAF0dHQ4PTQiWgQMG0TkuD/+4z/GD37wA/z0pz9FU1MTPve5zyEcDuPFF190emhEtAi4jEJEjtq7dy8efvhhPPPMMwiFQpBlGc888wz+4z/+A4899pjTwyOiRcCZDSIiIqoqzmwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVUxbBAREVFVMWwQERFRVTFsEBERUVX9f6Co9Je5ncPkAAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 600x400 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "dataStruct = loadmat('/content/drive/MyDrive/Colab Notebooks/regress1.mat')\n",
        "\n",
        "#let's first visualize the data\n",
        "x = dataStruct['x'].T\n",
        "y = dataStruct['y'].T\n",
        "\n",
        "fig, ax = plt.subplots(nrows=1, ncols=1,figsize= (6,4))\n",
        "plt.scatter(x, y, s = 100, color = cMAP1[0,], alpha = 0.5)\n",
        "plt.xlim([-2,3]); plt.ylim([-3,4]); plt.xlabel('x')\n",
        "plt.ylabel('y');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aKQAuHcRLlUT"
      },
      "source": [
        "Now, we will fit the data with polynomial linear models without using cross validation. This will serve as a comparison to see how cross validation alters our interpretation of the 'best' model."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "O0TXlWaJLlUT"
      },
      "outputs": [],
      "source": [
        "#Number of data points\n",
        "n_pts    = y.shape[1]\n",
        "x_zero   = np.ones([1,n_pts])\n",
        "#All regressors that we will use for each fit\n",
        "XX       = np.stack((x_zero[0], x[0], x[0]**2,x[0]**3,x[0]**4,x[0]**5,x[0]**6,\n",
        "                     x[0]**7,x[0]**8,x[0]**9,x[0]**10,x[0]**11), axis = 0)\n",
        "#Number of regressors (should be 6)\n",
        "n_models = 12\n",
        "order    = np.arange(0,n_models)\n",
        "#initialize matrices mse_vec (which stores mean squared error) and fit_mat (which stores\n",
        "#predicted y values)\n",
        "mse_vec  = []\n",
        "fit_mat  = []\n",
        "\n",
        "#------------------------YOUR CODE STARTS HERE-------------------------\n",
        "\n",
        "\n",
        "#----------------------------------------------------------------------"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UNdEl3RaLlUT"
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(nrows=3, ncols=4,figsize= (14,9))\n",
        "for i in range(n_models):\n",
        "    idx_row = i//4\n",
        "    idx_col = i%4\n",
        "    ax[idx_row, idx_col].scatter(x, y, s = 100, color = cMAP1[0,], alpha = 0.5)\n",
        "    ax[idx_row, idx_col].plot(x[0], fit_mat[i], color = cMAP1[1,], linewidth = 3)\n",
        "    ax[idx_row, idx_col].set_xlim([-2,3]); ax[idx_row, idx_col].set_ylim([-3,4]);\n",
        "    ax[idx_row, idx_col].set_xlabel('x'); ax[idx_row, idx_col].set_ylabel('y');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8FFK2hUSLlUT"
      },
      "source": [
        "Next let's re-do the model fitting with leave-one-out cross validation and then re-evaluate our models."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vygGeAXLLlUT"
      },
      "outputs": [],
      "source": [
        "#Let's fit this data with polynomials:\n",
        "\n",
        "#------------------------YOUR CODE STARTS HERE-----------------------------\n",
        "#Initialize a matrix to hold MSE values for each model for each cross validation\n",
        "mse_train = np.zeros([n_models, n_pts])\n",
        "mse_test  = np.zeros([n_models, n_pts])\n",
        "\n",
        "\n",
        "\n",
        "#find lowest MSE\n",
        "mean_mse_test = np.mean(mse_test,1)\n",
        "min_mse       = np.min(mean_mse_test)\n",
        "min_ind       = np.argmin(mean_mse_test)\n",
        "mean_mse_train = np.mean(mse_train,1)\n",
        "#----------------------------------------------------------------------------"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Ej8dKN79LlUT"
      },
      "outputs": [],
      "source": [
        "#Plot the mean MSE\n",
        "fig, ax = plt.subplots(nrows=1, ncols=1,figsize= (8,6))\n",
        "plt.scatter(min_ind, mean_mse_train[min_ind], s = 100, c = 'red')\n",
        "plt.plot(order, mean_mse_test, color = cMAP1[1],ls = '-', lw=3, ms=5)\n",
        "plt.plot(order, mean_mse_train, color = cMAP1[0],ls = '-', lw=3, ms=5)\n",
        "plt.xlim([-0.5,n_models-0.5]);  plt.xlabel('Order')\n",
        "plt.ylabel('Mean squared error (MSE)')\n",
        "plt.legend(['MSE of the ''best'' model','MSE (training set)', 'MSE (test set)'])\n",
        "plt.show()\n",
        "\n",
        "#Which model is best?\n",
        "print('The best model is when order = ' + str(order[min_ind]))\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Part 3 - Regularization"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this section, we will use Ridge (L2) regularization to reduce variance of an overparameterized model.\n",
        "\n",
        "First we need data generated by a groundtruth function. Let's say this function is a 7th degree polynomial with some added noise.\n",
        "\n",
        "Now, let's assume our model has 30 degrees of freedom (akin to a 29 degree polynomial). We can find the best fit coeffs by using np.polyfit (built-in function)\n",
        "\n",
        "The 29 degree polynomial gives us a very bad fit to the held-out test data. Let's regularize with ridge and lasso and observe how test error changes"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "np.random.seed(3003)\n",
        "# Generate train data\n",
        "N = 200\n",
        "x = np.linspace(-1, 1, N)\n",
        "\n",
        "true_y = (10 * x**7 - 8 * x**5 + 3 * x**3 - 2 * x  + 0.5)\n",
        "y = true_y + np.random.normal(0, 1., N)\n",
        "\n",
        "\n",
        "\n",
        "# Plot results\n",
        "plt.figure(figsize=(12, 6))\n",
        "\n",
        "plt.plot(x, y, marker='o', linestyle='', label='Train data')\n",
        "plt.plot(x, y_pred_poly, label=f'Polynomial (degree={degree})')\n",
        "plt.plot(x, y_pred_ridge, label=f'Ridge (alpha={ridge_alpha})')\n",
        "\n",
        "plt.xlabel('X')\n",
        "plt.ylabel('Y')\n",
        "plt.legend()\n",
        "plt.title('Comparison of models')\n",
        "plt.grid(True)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\n",
        "\n",
        "# Plot results\n",
        "plt.figure(figsize=(12, 6))\n",
        "\n",
        "plt.plot(x_test, y_test, marker='o', linestyle='', label='Test data')\n",
        "plt.plot(x_test, y_pred_poly, label=f'Polynomial (degree={degree})')\n",
        "plt.plot(x_test, y_pred_ridge, label=f'Ridge (alpha={ridge_alpha})')\n",
        "\n",
        "plt.xlabel('X')\n",
        "plt.ylabel('Y')\n",
        "plt.legend()\n",
        "# plt.xlim(-5, 15)\n",
        "plt.ylim(-100, 300)\n",
        "plt.title('Comparison of models')\n",
        "plt.grid(True)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Calculate mean squared error (MSE)\n",
        "mse_poly = np.mean((y_pred_poly - y_test)**2)\n",
        "mse_ridge = np.mean((y_pred_ridge - y_test)**2)\n",
        "\n",
        "# Plot MSE values\n",
        "plt.figure(figsize=(8, 6))\n",
        "plt.bar(['Polynomial', 'Ridge', 'Lasso'], [mse_poly, mse_ridge])\n",
        "plt.xlabel('Model')\n",
        "plt.ylabel('Mean Squared Error (MSE)')\n",
        "plt.yscale('log')\n",
        "plt.title('Comparison of MSE on held-out data')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Extract coefficients\n",
        "\n",
        "plt.figure(figsize=(12, 12))\n",
        "\n",
        "# Polynomial model\n",
        "plt.subplot(211)\n",
        "plt.stem(np.arange(degree + 1), coeffs_poly)\n",
        "plt.xlabel('Polynomial degree')\n",
        "plt.ylabel('Coefficient value')\n",
        "plt.title('Polynomial Coefficients')\n",
        "plt.grid(True)\n",
        "\n",
        "# Ridge model\n",
        "plt.subplot(212)\n",
        "plt.stem(np.arange(degree + 1), coeffs_ridge)\n",
        "plt.xlabel('Feature index')\n",
        "plt.ylabel('Coefficient value')\n",
        "plt.title('Ridge Coefficients (alpha={})'.format(ridge_alpha))\n",
        "plt.grid(True)\n",
        "\n",
        "\n",
        "plt.suptitle('Coefficients of Polynomial, Ridge, and Lasso Models')\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8qOa_iuRLlUV"
      },
      "source": [
        "# Part 4 - Permutation Testing\n",
        "\n",
        "In Master, Subramanian, & Zhao, 2005, Science, the authors ran the above task and computational models of behavior. After confirming accurate model recoverability with a generate & recover procedure, they found a significant correlation between Math Tools grade & attentional slope. They concluded that students who do well in Math Tools are better at allocating their attention to relevant information.\n",
        "\n",
        "Let's try and confirm their finding using a permutation test.\n",
        "The r-value they published in the paper is 0.23."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 73,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 470
        },
        "id": "m1Fu62octf9W",
        "outputId": "4733ea60-512b-405e-83ac-5478ef81fd2a"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, 'p = 0.117')"
            ]
          },
          "execution_count": 73,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGzCAYAAAAFROyYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQuUlEQVR4nO3de1yO9/8H8NfV3ekuFalUW+RQKVNOjZzK2BSGOYT1G5E55hSGbcjZzMYcYr5GmLMZhjFrymE5iyHnciwhSufDff/+aK651U3R7brS6/l4XI9dh8913e/7+oz77fpc1/sS1Gq1GkREREQyoid1AERERETPY4JCREREssMEhYiIiGSHCQoRERHJDhMUIiIikh0mKERERCQ7TFCIiIhIdpigEBERkewwQSEiIiLZYYJCREREssMEhYgk9/jxYwwYMADW1tYwNTVFq1atcOrUqWLte+zYMQwZMgQNGzaEgYEBBEHQ2nbJkiXo3r07qlatCkEQEBgYWGQ7Hx8fCIJQ5GRgYPAqX5GISkhf6gCIqHxTqVRo3749zpw5g7Fjx8LKygphYWHw8fHByZMn4eTk9ML9d+/ejeXLl8Pd3R01atTA5cuXtbb95ptv8OTJE7z//vtISEjQ2u6rr75C//79Ndalp6dj0KBB+Oijj0r2BYnolTBBISJJbdmyBX///Tc2b96Mbt26AQD8/f3h7OyMyZMnY926dS/cf/DgwRg3bhyUSiWCg4NfmKBERUWJV08qVKigtd2HH35YaN3PP/8MAAgICCjO1yKi18QhHqJyKDQ0FIIg4OLFi/D394e5uTkqV66MESNGICsr643GsmXLFlSpUgVdunQR11lbW8Pf3x/bt29Hdnb2C/evUqUKlEplsT6rWrVqLxwCepF169bB1NQUnTp1eqX9iahkeAWFqBzz9/eHo6MjZs2ahSNHjmDBggV49OgRVq9e/cL9MjIykJGR8dLjKxQKVKpU6YVtTp8+jQYNGkBPT/PfS++//z6WLVuGy5cvo27dui//Mjp0//597Nu3Dz169ICpqamksRCVF0xQiMqx6tWrY/v27QCAoUOHwtzcHGFhYRgzZgzc3d217jdnzhxMmTLlpcevVq0a4uPjX9gmISEBLVu2LLTezs4OAHD37l3JE5SNGzciLy+PwztEbxATFKJybOjQoRrLw4YNQ1hYGHbv3v3CBKV3795o3rz5S49fnKGXzMxMGBkZFVpvbGwsbpfaunXrYG1tXeS9KUSkG0xQiMqx55+QqVmzJvT09F561aNGjRqoUaNGqcSgVCqLvM/k6b0wxb2/RFeuX7+O6OhoBAcHQ1+ff2USvSn800ZEouLeQJqWloa0tLSXtlMoFLC2tn5hGzs7uyIf+X26zt7evlgx6crTp4g4vEP0ZjFBISrHrly5gurVq4vLV69ehUqlgqOj4wv3mzt3bqndg1KvXj0cPHgQKpVK40bZo0ePwsTEBM7Ozi/9HF1at24datasiSZNmkgaB1F5wwSFqBxbvHixRuGxhQsXAgD8/PxeuF9p3oPSrVs3bNmyBVu3bhXroDx48ACbN2/Gxx9/rHF/yrVr1wAUDEW9CadPn0ZsbCwmTpz4Rj6PiP7DBIWoHIuLi0PHjh3h6+uL6Oho/Pzzz/j000/h4eHxwv1K8x6Ubt26oUmTJujbty8uXLggVpLNz88vdJWmdevWAKBxVebGjRtYs2YNAODEiRMAgOnTpwMouILz2WefiW1/++03nDlzBgCQm5uLs2fPim07duxY6MbgtWvXAuDwDpEk1ERU7kyePFkNQH3hwgV1t27d1GZmZupKlSqpg4OD1ZmZmW88nuTkZHVQUJC6cuXKahMTE7W3t7f6+PHjhdpVq1ZNXa1aNY11+/fvVwMocvL29tZo26dPH61tV65cqdE2Pz9f/c4776gbNGhQyt+WiIpDUKvVakkyIyKSTGhoKKZMmYL79+/DyspK6nCIiAphqXsiIiKSHSYoREREJDtMUIiIiEh2eA8KERERyQ6voBAREZHsMEEhIiIi2SmThdpUKhXu3r0LMzOzYr87hIiIiKSlVqvx5MkT2Nvba7zaoihlMkG5e/cuHBwcpA6DiIiIXsGtW7fw7rvvvrBNmUxQzMzMABR8QXNzc4mjISIiouJITU2Fg4OD+Dv+ImUyQXk6rGNubs4EhYiIqIwpzu0ZvEmWiIiIZIcJChEREckOExQiIiKSnTJ5D0pxqNVq5OXlIT8/X+pQiMoshUIBfX19Ps5PRG/cW5mg5OTkICEhARkZGVKHQlTmmZiYwM7ODoaGhlKHQkTlyFuXoKhUKsTFxUGhUMDe3h6Ghob81x/RK1Cr1cjJycH9+/cRFxcHJyenlxZWIiIqLW9dgpKTkwOVSgUHBweYmJhIHQ5RmaZUKmFgYIAbN24gJycHxsbGUodEROXEW/vPIf5Lj6h08M8SEUmBf/MQERGR7DBBISIiItlhgvIWCgwMROfOncVlHx8fjBw5skTHCA0NRb169bQek4iISJeYoMhEYGAgBEHA7NmzNdZv27ZNFk8h/fDDDwgPDy9WWyYzRET0upigyIixsTG++eYbPHr0SOpQCrGwsEDFihWlDoOIiMqJ8pWgpKdrn7Kyit82M/PlbV9BmzZtYGtri1mzZmlt8/zQCwDMnz8fjo6Or/SZT82ePRtVqlSBmZkZgoKCkPXc+Xj+qsiWLVtQt25dKJVKVK5cGW3atEF6ejpCQ0OxatUqbN++HYIgQBAEREZGvlZsRERU/pSvBKVCBe1T166abW1stLf189Ns6+hYuM0rUCgUmDlzJhYuXIjbt2+/2nd8BZs2bUJoaChmzpyJEydOwM7ODmFhYVrbJyQkoFevXujXrx9iY2MRGRmJLl26QK1WY8yYMfD394evry8SEhKQkJCApk2bvrHvQkREb4e3rlBbWffJJ5+gXr16mDx5Mn766ac38pnz589HUFAQgoKCAADTp0/Hn3/+WegqylMJCQnIy8tDly5dUK1aNQBA3bp1xe1KpRLZ2dmwtbXVffBERPRWKl9XUNLStE+//KLZNilJe9vff9dsGx9fuM1r+Oabb7Bq1SrExsa+1nGKKzY2Fo0bN9ZY5+XlpbW9h4cHWrdujbp166J79+743//+J8v7ZohkIz0dEISC6RWHgMuD9PR0cWg4nedJMnLph/KVoJiaap+eL+H9orZK5cvbvoaWLVuibdu2mDBhQqFtenp6UKvVGutyc3Nf6/NKSqFQYN++ffj999/h5uaGhQsXwsXFBXFxcW80DiIienuVrwSlDJk9ezZ+++03REdHa6y3trZGYmKiRpISExPzWp/l6uqKo0ePaqw7cuTIC/cRBAHNmjXDlClTcPr0aRgaGuLXX38FABgaGiI/P/+1YiIiovKN96DIVN26dREQEIAFCxZorPfx8cH9+/cxZ84cdOvWDXv27MHvv/8Oc3PzV/6sESNGIDAwEI0aNUKzZs2wdu1anD9/HjVq1Ciy/dGjRxEREYGPPvoINjY2OHr0KO7fvw9XV1cAgKOjI/bu3YtLly6hcuXKsLCwgIGBwSvHR0RE5Q+voMjY1KlToVKpNNa5uroiLCwMixcvhoeHB44dO4YxY8a81uf06NEDEydOxBdffIGGDRvixo0bGDx4sNb25ubmOHDgANq1awdnZ2d8/fXX+O677+D379NNn3/+OVxcXNCoUSNYW1vj8OHDrxUfERGVP4L6+RsayoDU1FRYWFggJSWl0JWDrKwsxMXFoXr16nw1PFEp4J+pUpCe/l/5gbS0175P7W2Vnp6OCv+ep7S0NJjyPElCl/3wot/v5/EKChEREckO70EhItI1hQJo1+6/eSqSQqFAu3/Pk4LnSTJy6QcmKEREumZsDOzaJXUUsmdsbIxdPE+Sk0s/cIiHiIiIZIcJChEREckOExQiIl1LT/+vyjRLuGuVnp4OU1NTmJqastS9hOTSD7wHhYjoTcjIkDqCMiGD50kW5NAPvIJCREREssMEhYiIiGSHQzxvifj4eFSvXh2nT59GvXr1irVPeHg4Ro4cicePH0sah9RK4zw4Ojpi5MiRGDly5GvHc/b2q8ehC+q8HCQ9ykT/rZG48yQf8bPbSx0SEZUDvIIiI7du3UK/fv1gb28PQ0NDVKtWDSNGjMDDhw9fuq+DgwMSEhLw3nvvFfvzevTogcuXL79OyJIRBAHbtm0r8X6Ojo6YP3++xrqSnIfw8HBUrFix0Prjx49jwIABJY6HiIiKxgRFJq5fv45GjRrhypUrWL9+Pa5evYqlS5ciIiICXl5eSE5O1rpvTk4OFAoFbG1toa9f/ItiSqUSNjY2pRH+a8nJyZH080vjPFhbW8PExKSUIiIiovKRoKjVQE76m59K8B7GoUOHwtDQEH/88Qe8vb1RtWpV+Pn54c8//8SdO3fw1VdfiW0dHR0xbdo09O7dG+bm5hgwYADi4+MhCAJiYmLEdjt27ICTkxOMjY3RqlUrrFq1CoIgiEMZz18NCA0NRb169bBmzRo4OjrCwsICPXv2xJMnT8Q2e/bsQfPmzVGxYkVUrlwZHTp0wLVr10rUHUXF/7ycnBwEBwfDzs4OxsbGqFatGmbNmiXuDwCffPIJBEEQl69du4ZOnTqhSpUqqFChAjw9PfHnn3+Kx/Tx8cGNGzcwatQoCIIAQRCKPA9nzpxBq1atYGZmBnNzczRs2BAnTpxAZGQk+vbti5SUFHH/0NBQMaZnr8w8fvwYAwcORJUqVWBsbIz33nsPO3fuLNF5oreInh7g7V0w6ZWPv3ZfhZ6eHry9veHt7Q09nifJyKUfysc9KLkZwEz7N/+5X94FDF/+Fsjk5GTs3bsXM2bMgFKp1Nhma2uLgIAAbNy4EWFhYeKP6ty5czFp0iRMnjy5yGPGxcWhW7duGDFiBPr374/Tp09jzJgxL43l2rVr2LZtG3bu3IlHjx7B398fs2fPxowZMwAUPB8fEhICd3d3pKWlYdKkSfjkk08QExNTov+RXxb/ggULsGPHDmzatAlVq1bFrVu3cOvWLQAFwyk2NjZYuXIlfH19xXdFpKWloV27dpgxYwaMjIywevVqfPzxx7h06RKqVq2KrVu3wsPDAwMGDMDnn3+uNbaAgADUr18fS5YsgUKhQExMDAwMDNC0aVPMnz8fkyZNwqVLlwBAfOPns1QqFfz8/PDkyRP8/PPPqFmzJi5cuMB3i5RnSiUQGSl1FLKnVCoRyfMkObn0Q/lIUGTuypUrUKvVcHV1LXK7q6srHj16hPv374tDER988AFGjx4ttomPj9fY58cff4SLiwu+/fZbAICLiwvOnTsnJhraqFQqhIeHw8zMDADw2WefISIiQtyva9euGu1XrFgBa2trXLhwoUT3vzwf//Nu3rwJJycnNG/eHIIgoFq1auI2a2trAEDFihVha2srrvfw8ICHh4e4PG3aNPz666/YsWMHgoODYWlpCYVCATMzM439ivrssWPHonbt2gAAJycncZuFhQUEQXjh/n/++SeOHTuG2NhYODs7AwBq1KihtT0RERVWPhIUA5OCqxlSfG4JqEswJNSoUaMXbr906RI8PT011r3//vsvPa6jo6OYnACAnZ0dkpKSxOUrV65g0qRJOHr0KB48eACVSgWg4Ee9JAnKs/EPGjQIP//8s7iclpaGwMBAfPjhh3BxcYGvry86dOiAjz766IXHTEtLQ2hoKHbt2oWEhATk5eUhMzMTN2/eLHZcABASEoL+/ftjzZo1aNOmDbp3746aNWsWe/+YmBi8++67YnJCREQlVz4G+QShYKjlTU//Dse8TK1atSAIAmJjY4vcHhsbi0qVKolXDgDA1PTlQ0evwsDAQGNZEAQxCQGAjz/+GMnJyfjf//6Ho0eP4ujRowBKfqPrs/FPnToVMTEx4gQADRo0QFxcHKZNm4bMzEz4+/ujW7duLzzmmDFj8Ouvv2LmzJk4ePAgYmJiULdu3RLHFhoaivPnz6N9+/b466+/4Obmhl9//bXY+z8/TEeE9HTA2rpgYgl3rdLT02FtbQ1ra2uWupeQXPqhfCQoMle5cmV8+OGHCAsLQ2Zmpsa2xMRErF27Fj169BDvPykOFxcXnDhxQmPd8ePHXyvOhw8f4tKlS/j666/RunVrcejpddnY2KBWrVri9JS5uTl69OiB//3vf9i4cSN++eUX8WkmAwMD5Ofnaxzn8OHDCAwMxCeffIK6devC1ta20NCXoaFhof2K4uzsjFGjRuGPP/5Aly5dsHLlymLv7+7ujtu3b5fZR7hJRx48KJjohR48eIAHPE+Sk0M/MEGRiUWLFiE7Oxtt27bFgQMHcOvWLezZswcffvgh3nnnnZfeO/K8gQMH4uLFixg3bhwuX76MTZs2ITw8HABKlOg8q1KlSqhcuTKWLVuGq1ev4q+//kJISMgrHetlvv/+e6xfvx4XL17E5cuXsXnzZtja2opP2zg6OiIiIgKJiYlikuTk5IStW7ciJiYGZ86cwaeffqpx9efpfgcOHMCdO3eK/MOXmZmJ4OBgREZG4saNGzh8+DCOHz8u3h/k6OiItLQ0RERE4MGDB0W+r8Lb2xstW7ZE165dsW/fPsTFxeH333/Hnj17SvksERG9vZigyISTkxNOnDiBGjVqwN/fHzVr1sSAAQPQqlUrREdHw9LSskTHq169OrZs2YKtW7fC3d0dS5YsER9VNjIyeqUY9fT0sGHDBpw8eRLvvfceRo0aJd6EW9rMzMwwZ84cNGrUCJ6enoiPj8fu3bvFJ4W+++477Nu3Dw4ODqhfvz6AgqSmUqVKaNq0KT7++GO0bdsWDRo00Dju1KlTER8fj5o1a2oMmT2lUCjw8OFD9O7dG87OzvD394efnx+mTJkCAGjatCkGDRqEHj16wNraGnPmzCky/l9++QWenp7o1asX3Nzc8MUXXxTryg0RERUQ1CW5M1MmUlNTYWFhgZSUFJibm2tsy8rKQlxcHKpXrw5jY2OJIpSnGTNmYOnSpeLjuiRPsix1f/c2Qvcn4c4T+SdZsizFn54OPH0kPS0N0NE9ZGVdenq6+Oh+Wlqazu61oxfTZT+86Pf7eSW6gjJr1ix4enrCzMwMNjY26Ny5s1gP4qmsrCwMHToUlStXRoUKFdC1a1fcu3dPo83NmzfRvn17mJiYwMbGBmPHjkVeXl5JQqFiCAsLw/Hjx3H9+nWsWbMG3377Lfr06SN1WERERC9VogQlKioKQ4cOxZEjR7Bv3z7k5ubio48+0rjLd9SoUfjtt9+wefNmREVF4e7du+jSpYu4PT8/H+3bt0dOTg7+/vtvrFq1CuHh4Zg0aVLpfSsCUPBIcKdOneDm5oZp06Zh9OjRYuVTIiIiOXutIZ6nhcOioqLQsmVLpKSkwNraGuvWrRMfCb148SJcXV0RHR2NJk2a4Pfff0eHDh1w9+5dVKlSBQCwdOlSjBs3Dvfv34ehoeFLP5dDPPQ24xDP65HlEE9mJtCyZcH8gQMFlWWpkMzMTLT89zwdOHCAj+xLRJf9UJIhntcq1JaSkgIA4g2cJ0+eRG5uLtq0aSO2qV27NqpWrSomKNHR0ahbt66YnABA27ZtMXjwYJw/f1684fFZ2dnZyM7O1viCRERlhlIJvOZj/uWBUql87XII9Prk0g+v/BSPSqXCyJEj0axZM7GCaGJiIgwNDQu9jr5KlSpITEwU2zybnDzd/nRbUWbNmgULCwtxcnBweNWwiYiIqAx45QRl6NChOHfuHDZs2FCa8RRpwoQJSElJESc+hUJERPR2e6UEJTg4GDt37sT+/fvx7rvviuttbW2Rk5ODx48fa7S/d++e+HI1W1vbQk/1PF3W9gI2IyMjmJuba0xERGVGRgbg6FgwFVHcjwpkZGTA0dERjo6ORRZBpDdDLv1QogRFrVYjODgYv/76K/766y9Ur15dY3vDhg1hYGCAiIgIcd2lS5dw8+ZNeHl5AQC8vLzwzz//aLyAbt++fTA3N4ebm9vrfBciInlSq4EbNwqmsld66o1Rq9W4ceMGbty4UaKXp1Lpkks/lOgm2aFDh2LdunXYvn07zMzMxHtGLCwsoFQqYWFhgaCgIISEhMDS0hLm5uYYNmwYvLy80KRJEwDARx99BDc3N3z22WeYM2cOEhMT8fXXX2Po0KGvXOGUiIiI3i4luoKyZMkSpKSkwMfHB3Z2duK0ceNGsc28efPQoUMHdO3aFS1btoStrS22bt0qblcoFNi5cycUCgW8vLzwf//3f+jduzemTp1aet+KSEfUajUGDBgAS0tLCIIgvn35RQRBwLZt23Qe26u4c+smPBwq4eL5f6QOhYhIQ4mHeIqaAgMDxTbGxsZYvHgxkpOTkZ6ejq1btxa6t6RatWrYvXs3MjIycP/+fcydOxf6+q/1xHOZFxgYiM6dO0sdRok9efIEI0eORLVq1aBUKtG0adNCj6cFBgZCEASNydfXV9yenZ2Nzz77DObm5nB2dsaff/6psf+3336LYcOGvZHv8zJ79uxBeHg4du7ciYSEBPEJNiIiKl3lOyugYsvJySmyiF7//v1x7tw5rFmzBvb29vj555/Rpk0bXLhwAe+8847YztfXFytXrhSXnx3OW7ZsGU6ePIno6Gj8/vvv+PTTT3Hv3j0IgoC4uDj873//w4kTJ3T7BYvp2rVrsLOzQ9OmTaUOhYjorVau3macnp6udcrKyip228zMzJe2LW1RUVF4//33YWRkBDs7O4wfP158f9HOnTtRsWJF8W25MTExEAQB48ePF/fv378//u///k9cPnToEFq0aAGlUgkHBwcMHz5cI25HR0dMmzYNvXv3hrm5OQYMGFAopszMTPzyyy+YM2cOWrZsiVq1aiE0NBS1atXCkiVLNNoaGRnB1tZWnCpVqiRui42NRceOHVGnTh0MHToU9+/fx4MHDwAAgwcPxjfffFPsJ7dWrFiBOnXqiOcpODhY3Hbz5k106tQJFSpUgLm5Ofz9/TWeKAsNDUW9evWwZs0aODo6wsLCAj179sSTJ08AFFwJGjZsGG7evAlBEODo6FismAAgISEBfn5+UCqVqFGjBrZs2aKx/Z9//sEHH3wApVKJlnVrYOq4kchITxO3B3XvgDmhEzT2GRkUgImjhojLfl7uWL7wO0waHQyv2g5o2/g9bFkbrvk5p0/C37clPGvZole7Vrh47myxvwMR0ZtUrhKUChUqaJ26du2q0dbGxkZrWz8/P422jo6OhdqUpjt37qBdu3bw9PTEmTNnsGTJEvz000+YPn06AKBFixZ48uQJTp8+DaAgmbGyskJkZKR4jKioKPj4+AAouArg6+uLrl274uzZs9i4cSMOHTqk8WMOAHPnzoWHhwdOnz6NiRMnFoorLy8P+fn5hV4poFQqcejQIY11kZGRsLGxgYuLCwYPHoyHDx+K2zw8PHDo0CFkZmZi7969sLOzg5WVFdauXQtjY2N88sknxTpPS5YswdChQzFgwAD8888/2LFjB2rVqgWgoLBgp06dkJycjKioKOzbtw/Xr19Hjx49NI5x7do1bNu2DTt37sTOnTsRFRWF2bNnAwB++OEHTJ06Fe+++y4SEhJKVGlx4sSJ6Nq1K86cOYOAgAD07NkTsbGxAAoS3LZt26JSpUo4fvw4vl0ajiOHIjHr6y+KffynVi9bjDru9bDx9yj49w7CjC9HI/7aFQBARnoahvXtiZpOLli/az8GhYzD99ML9yvpgCAAbm4FkyBIHY1sCYIANzc3uLm5QeB5koxc+oFDPGVAWFgYHBwcsGjRIgiCgNq1a+Pu3bsYN24cJk2aBAsLC9SrVw+RkZFo1KgRIiMjMWrUKEyZMgVpaWlISUnB1atX4e3tDaCgMm9AQABGjhwJAHBycsKCBQvg7e2NJUuWiAnHBx98gNGjR2uNy8zMDF5eXpg2bRpcXV1RpUoVrF+/HtHR0WJiABQM73Tp0gXVq1fHtWvX8OWXX8LPzw/R0dFQKBTo168fzp49Czc3N1hZWWHTpk149OgRJk2ahMjISHz99dfYsGEDatasiRUrVmgMHT1r+vTpGD16NEaMGCGu8/T0BABERETgn3/+QVxcnFiJePXq1ahTpw6OHz8utlOpVAgPD4eZmRkA4LPPPkNERARmzJgBCwsLmJmZQaFQaK3Zo0337t3Rv39/AMC0adOwb98+LFy4EGFhYVi3bh2ysrKwevVqmJqaQlXxXUyYNgfD+/bCyC9DUdnaptif0/yDD9GjT8Hn9BsyEj8vX4Jjfx+EY00n7N62BWqVCqHfLoSRsTFqubjiXsJdzPhSex9TKTExAc6flzoK2TMxMcF5nifJyaUfylWCkpaWpnWbQqHQWH62Tsvz9PQ0LzzFx8e/VlwvExsbCy8vL41MtlmzZkhLS8Pt27dRtWpVeHt7IzIyEqNHj8bBgwcxa9YsbNq0CYcOHUJycjLs7e3h5OQEADhz5gzOnj2LtWvXisdTq9VQqVSIi4uDq6srAKBRo0YvjW3NmjXo168f3nnnHSgUCjRo0AC9evXCyZMnxTY9e/YU5+vWrQt3d3fUrFkTkZGRaN26NQwMDLB48WKN4/bt2xfDhw/H6dOnsW3bNpw5cwZz5szB8OHD8csvvxSKIykpCXfv3kXr1q21nkMHBweN1yS4ubmhYsWKiI2NFRMUR0dHMTkBADs7uxf+v1BcT+sAPbv89Amg2NhYeHh4wNTUVNxer1FjqFQqxF+7UqIExdm1jjgvCAKsrG2Q/LBguCzuymU4udaB0TNXvDwaer7K1yEi0rlylaA8+wMgVVtd8fHxwYoVK3DmzBkYGBigdu3a8PHxQWRkJB49eiRePQEKErWBAwdi+PDhhY5TtWpVcb4436tmzZqIiopCeno6UlNTYWdnhx49eqBGjRpa96lRowasrKxw9erVIhOK/fv34/z581i+fDnGjh2Ldu3awdTUFP7+/li0aFGRxyytt20aGBhoLAuCAJVKVSrHfh2Cnl6hgkl5ebmF2unrF45fLYP4iYhKqlzdg1JWubq6Ijo6WuMH6vDhwzAzMxNfNfD0PpR58+aJycjTBCUyMlK8/wQAGjRogAsXLqBWrVqFpqKe1CkOU1NT2NnZ4dGjR9i7dy86deqkte3t27fx8OFD2NnZFdqWlZWFoUOH4scff4RCoUB+fj5ycwt+iHNzc8UbgZ9nZmYGR0dHjSrGz3J1dcWtW7c03uN04cIFPH78+I1UMD5y5Eih5adXqlxdXXHmzBmNm5RjThyFnp4eHGsWXPWqZGmFB0n/3dCbn5+Pq5diSxRDdSdnXIk9j+xnbgg/e0oeT0e99TIygDp1CiaWcNcqIyMDderUQZ06dVjqXkJy6QcmKDKSkpKCmJgYjenWrVsYMmQIbt26hWHDhuHixYvYvn07Jk+ejJCQEHG4qVKlSnB3d8fatWvFZKRly5Y4deoULl++rHEFZdy4cfj7778RHByMmJgYXLlyBdu3by90k2xx7N27F3v27EFcXBz27duHVq1aoXbt2ujbty+Agqs1Y8eOxZEjRxAfH4+IiAh06tQJtWrVQtu2bQsdb9q0aWjXrh3q168PoGAoa+vWrTh79iwWLVqEZs2aaY0lNDQU3333HRYsWIArV67g1KlTWLhwIQCgTZs2qFu3LgICAnDq1CkcO3YMvXv3hre3d7GGsrRZtGiR1mGlZ23evBkrVqzA5cuXMXnyZBw7dkw83wEBATA2NkafPn1w7tw5HPv7IGZPHIcOXXqIwzvvN2uBAxF/4EDEXsRdvYwZX47Gk9SUEsXarnM3CIKAKeNG4Nrlizj41x9Y/WPRV6SolKnVwIULBRNLuGulVqtx4cIFXLhwgaXuJSSXfihXQzxyFxkZKf4wPxUUFITly5dj9+7dGDt2LDw8PGBpaYmgoCB8/fXXGm29vb0RExMjJiiWlpZwc3PDvXv34OLiIrZzd3dHVFQUvvrqK7Ro0QJqtRo1a9Ys9ERLcaSkpGDChAm4ffs2LC0t0bVrV8yYMUMcKlEoFDh79ixWrVqFx48fw97eHh999BGmTZtW6NUG586dw6ZNmzSqs3br1g2RkZFo0aIFXFxcsG7dOq2x9OnTB1lZWZg3bx7GjBkDKysrdOvWDUDBUMf27dsxbNgwtGzZEnp6evD19RUTmFf14MEDXLt27aXtpkyZgg0bNmDIkCGws7PD+vXrxSs3JiYm2Lt3L0aMGAFPT08YGSvRpl1HjJk0Xdy/c4//w+UL5/D1yMFQ6Ovj//oPhqdXixLFamJaAQtWrsf0CSHo4eeNGk4uGPFlKEYP6F2yL01E9AYI6jKYpqampsLCwgIpKSmF6mNkZWUhLi4O1atXL/T4K1FZcPb2Y6lD0KDOy0HS3dsI3Z+EO0+KHmKTk/jZ7aUOobD0dOBp+YG0NEAG963JUXp6ulimIS0tTRb395VHuuyHF/1+P49DPERERCQ7TFCIiIhIdpigEBERkezwJlkiIl0TBKBatf/mqUiCIKDav+eJpe6lI5d+YIJCRKRrJiaAjitOvw1MTEx0XpmbXk4u/cAhHiIiIpIdJihEREQkO0xQiIh0LTMT8PQsmDIzpY5GtjIzM+Hp6QlPT09k8jxJRi79wHtQiIh0TaUCTpz4b56KpFKpcOLf8ySHl3SWV3LpB15BISqhw4cPo27dujAwMEDnzp1f2j40NBT16tXTeVyvauKoIRgZFCB1GEREGpigyISPjw9GjhxZaH14eDgqVqz4xuMpiYiICDRt2hRmZmawtbXFuHHjkJeXJ26Pj4+HIAiFpmff8Ltv3z44OzvD3Nwcn332GXJycsRtKSkpcHZ2xo0bN97o99ImJCQE9erVQ1xcHMLDw6UOh4jorcQEhYolPz+/yEt9Z86cQbt27eDr64vTp09j48aN2LFjB8aPH1+o7Z9//omEhARxatiwIYCCS4iffvopBg0ahOjoaJw4cQLLli0T9xs/fjwGDRokPpcvtWvXruGDDz7Au+++K/vkkYiorCpXCUp6errWKSsrq9htn79pqKg2uhIYGIjOnTtjypQpsLa2hrm5OQYNGqRxxcHHxwfBwcEIDg6GhYUFrKysMHHiRI3XZmdnZ2PMmDF45513YGpqisaNGyMyMlLc/vTKzY4dO+Dm5gYjIyPcvHmzUDwbN26Eu7s7Jk2ahFq1asHb2xtz5szB4sWL8eTJE422lStXhq2trTg9fePxgwcP8ODBAwwZMgR16tRBx44dERsbCwD4+++/cfz4cYwYMaJY5+f27dvo1asXLC0tYWpqikaNGuHo0aPi9iVLlqBmzZowNDSEi4sL1qxZo7G/IAhYvnw5PvnkE5iYmMDJyQk7duwA8N+VoIcPH6Jfv34QBKFEV1B+/PFHODg4wMTEBP7+/khJSRG3qVQqTJ06Fe+++y4a1awC/7YtcHj/n+L249GH4OFQCanP7HPx/D/wcKiEO7cK+mX7pnVoXqcaDkdGoHOrxmji8i4G/1833L+XKO6Tn5+Pb6d8heZ1qqFl3RqYN2MS1Chz7wslonKgXCUoFSpU0Dp17dpVo62NjY3Wtn5+fhptHR0dC7XRpYiICMTGxiIyMhLr16/H1q1bMWXKFI02q1atgr6+Po4dO4YffvgB33//PZYvXy5uDw4ORnR0NDZs2ICzZ8+ie/fu8PX1xZUrV8Q2GRkZ+Oabb7B8+XKcP38eNjY2hWLJzs4u9NZopVKJrKwsnDx5UmN9x44dYWNjg+bNm4s/+gBgbW0NOzs7/PHHH8jIyMDBgwfh7u6O3NxcDB48GD/++CMUCsVLz0taWhq8vb1x584d7NixA2fOnMEXX3whXvn59ddfMWLECIwePRrnzp3DwIED0bdvX+zfv1/jOFOmTIG/vz/Onj2Ldu3aISAgAMnJyXBwcEBCQgLMzc0xf/58JCQkoEePHi+NCwCuXr2KTZs24bfffsOePXtw+vRpDBkyRNz+ww8/4LvvvsPcuXOx5Y9D8PL+AMODPsWNuGvFOv5TmZmZWL1sEWbMX4qVW3Yh8e5tfD99orh99bJF2LF5HULnLkL41t+R8vgx/tqzq0SfQUT0JpSrBOVtYWhoiBUrVqBOnTpo3749pk6digULFmgMwTg4OGDevHlwcXFBQEAAhg0bhnnz5gEAbt68iZUrV2Lz5s1o0aIFatasiTFjxqB58+ZYuXKleIzc3FyEhYWhadOmcHFxgYmJSaFY2rZti7///hvr169Hfn4+7ty5g6lTpwIAEhISABQkht999x02b96MXbt2oXnz5ujcubOYpAiCgE2bNmHatGmoU6cO6tevj379+mH27Nlo1aoVjI2N0axZM7i4uGDRokVaz8u6detw//59bNu2Dc2bN0etWrXg7+8PLy8vAMDcuXMRGBiIIUOGwNnZGSEhIejSpQvmzp2rcZzAwED06tULtWrVwsyZM5GWloZjx45BoVDA1tYWgiDAwsICtra2UCqVxeqzrKwsrF69GvXq1UPLli2xcOFCbNiwAYmJiWJs48aNQ8+ePeFY0wmjvpwCF7e6WLt8SbGO/1Rebi6+nvk96njUh2tdD/Ts0x9HDx8Qt69dvhRBwaPQxu9j1HBywdezvkcFM7MSfQa9IiurgoleyMrKClY8T5KTQz+Uq8eM09LStG57/l/oSUlJWtvq6WnmdW+6JLCHh4dGsuDl5YW0tDTcunVLvE+jSZMmGu9Q8PLywnfffYf8/Hz8888/yM/Ph7Ozs8Zxs7OzUblyZXHZ0NAQ7u7uL4zlo48+wrfffotBgwbhs88+g5GRESZOnIiDBw+K58nKygohISHiPp6enrh79y6+/fZbdOzYEQDQvHlzHD9+XGxz+fJlrF69GqdPn0bLli0xYsQI+Pn54b333kPLli2LjCsmJgb169eHpaVlkbHGxsZiwIABGuuaNWuGH374QWPds8c2NTWFubn5C/9/KI6qVavinXfeEZe9vLygUqlw6dIlmJiY4O7du2jWrJnGPvUbNcal2HMl+hxjpQkcHKuLy1Y2tkh+cB8A8CQ1BfeTElG3XkNxu76+Ptzc6wNqDvPolKkpcP++1FHInqmpKe7zPElOLv1QrhIUU1NTydtqY25urnFPwlOPHz+GhYXFax//WWlpaVAoFDh58mShxOzZ4SmlUlmsF0WFhIRg1KhRSEhIQKVKlRAfH48JEyagRo0aWvdp3Lgx9u3bp3X7wIED8d1330GlUuH06dPo3r07TExM4O3tjaioqCITlOJezXiZp/fGPCUIguQ1GfSe9sMziURebm6hdgYGmn+kBUHQuPeIiKis4BCPTLi4uODUqVOF1p86darQlY4zZ85o3Kh75MgRVKhQAQ4ODuK6Z28MfdrGyckJCoUC9evXR35+PpKSklCrVi2NydbW9pXiFwQB9vb2UCqVWL9+PRwcHNCgQQOt7WNiYmBnZ1fktp9++gmWlpbo2LEj8vPzARQMNz3979N1z3N3d0dMTAySk5OL3O7q6orDhw9rrDt8+DDc3Nxe+v1e182bN3H37l1x+ciRI9DT04OLiwvMzc1hb29fKLbTJ46ihpMLAKBS5YJLrfeT/rvh9dKFf0oUg5m5BaxtbPFPzH/3BuXl5SH2n5iSfh0iIp1jgiITgwcPxuXLlzF8+HCcPXsWly5dwvfff4/169dj9OjRGm1zcnIQFBSECxcuYPfu3Zg8eTKCg4M1hp5u3ryJkJAQXLp0CevXr8fChQvFJ2GcnZ0REBCA3r17Y+vWrYiLi8OxY8cwa9Ys7NpV8hsmv/32W/zzzz84f/48pk2bhtmzZ2PBggXi1ZlVq1Zh/fr1uHjxIi5evIiZM2dixYoVGDZsWKFjJSUlYfr06Vi4cCEAoFKlSnB1dcX8+fMRHR2NiIiIQkMhT/Xq1Qu2trbo3LkzDh8+jOvXr+OXX35BdHQ0AGDs2LEIDw/HkiVLcOXKFXz//ffYunUrxowZU+Lv/KzatWvj119/fWEbY2Nj9OnTB2fOnMHBgwcxfPhw+Pv7iwnh2LFj8c0332Djxo2Iv3YF82eF4tKFfxAQNAgA4OBYA7b272DpvG9wI+4aDkTsxepli0sc66dBA7Fi8Xz8tWcX4q5exsyvxuBJamrJvzSVTGYm4ONTMLGEu1aZmZnw8fGBj48PS91LSC79UK6GeOSsRo0aOHDgAL766iu0adMGOTk5qF27NjZv3gxfX1+Ntq1bt4aTkxNatmyJ7Oxs9OrVC6GhoRptevfujczMTLz//vtQKBQYMWKExv0XK1euxPTp0zF69GjcuXMHVlZWaNKkCTp06FDi2H///XfMmDED2dnZ8PDwwPbt2ws96TRt2jTcuHED+vr6qF27NjZu3Ihu3boVOtbTp2zs7e3FdeHh4ejTpw8WLFiAsWPHwtPTs8g4DA0N8ccff2D06NFo164d8vLy4ObmhsWLC37IO3fujB9++AFz587FiBEjUL16daxcuRI+Pj4l/s7PunTpUpHDc8+qVasWunTpgnbt2iE5ORkdOnRAWFiYuH348OFISUnB6NGjcS8pCTWdXLDgp3WoVr0mgIJhp9mLlmPGl6PR/cPmqONRH8Fjv8KYQYElirX3gGA8SLqHiSGDIejpobP//+ED3/ZIY5KiWyoVEBX13zwVSaVSIerf8yT1sGp5Jpd+ENRlcIA6NTUVFhYWSElJgbm5uca2rKwsxMXFoXr16oUef30bBAYG4vHjx9i2bZvWNj4+PqhXrx7mz5//xuKi0nP29mOpQ9CgzstB0t3bCN2fhDtPih5ek5P42e2lDqGw9HTg6f1daWkFN81SIenp6eJ9cGlpaaVyfx+VnC774UW/38/jEA8RERHJDhMUIiIikh3eg1LGFKe0+rMl64mIiMoiXkEhIiIi2Xlrr6CUwXt/ieRJrQaghop/pF5PEa+KoMKKeqUGvXly6Ie3LkF5WgU0IyOj1CqLEpVn6rwc5Oar8SiLj32+MlPTgid56IVMTU11+jZ4Kh659MNbl6AoFApUrFhRfHeKiYlJscq1E8mFOi9H6hAKqNVQ5+XgUfIDRFxPQ1YeL6EQ0Zvz1iUoAMTqnK/7gjciKSQ9kksFTTVy89WIuJ6GrbHS/2uKiMqXtzJBEQQBdnZ2sLGxEd/hQlRW9N8aKXUIAACVGniUpeKVk9KQlQV07Vow/8svwFtYRLI0ZGVloeu/5+mXX355K4ttlgVy6Ye3MkF5SqFQFHpbL5HclYVqrVRC+fnA7t3/zVOR8vPzsfvf86TtpaCke3Lph7c6QSEqiuP4kr8QkYiI3izWQSEiIiLZYYJCREREssMEhYiIiGSHCQoRERHJDhMUIiIikh0+xUNEpGumpv++04hexNTUlO9RkwG59AOvoBAREZHsMEEhIiIi2WGCQkSka1lZQPfuBVNWltTRyFZWVha6d++O7t27I4vnSTJy6QdBLYeBphJKTU2FhYUFUlJSYG5uLnU4VMawkuzbLX52e6lDKCw9HahQoWA+La3gnhQqJD09HRX+PU9paWkw5XmShC77oSS/37yCQkRERLLDBIWIiIhkhwkKERERyQ4TFCIiIpIdJihEREQkO0xQiIiISHZY6p6ISNdMTAoeL346T0UyMTFB2r/nyYTnSTJy6QcmKEREuiYIrH1SDIIgsPaJDMilH5igENFbRe6F+GRZSI5IhpigEBHpmGFeLmbuXVSwMKUNYGQkbUAylZ2djYEDBwIAfvzxRxjxPElCLv3AUvdU7sj9X9j09lHmZCF2XreCBZa614ql7uWBpe6JiIiItGCCQkRERLLDBIWIiIhkhwkKERERyQ4TFCIiIpIdJihEREQkOyVOUA4cOICPP/4Y9vb2EAQB27Zt09geGBgIQRA0Jl9fX402ycnJCAgIgLm5OSpWrIigoCCxrC4R0dsm08AIDYatRYNha1nq/gVMTEyQlJSEpKQklrqXkFz6ocSF2tLT0+Hh4YF+/fqhS5cuRbbx9fXFypUrxeXni7wEBAQgISEB+/btQ25uLvr27YsBAwZg3bp1JQ2HiEj+BAHJJhbiPBVNEARYW1tLHUa5J5d+KHGC4ufnBz8/vxe2MTIygq2tbZHbYmNjsWfPHhw/fhyNGjUCACxcuBDt2rXD3LlzYW9vX9KQiIiI6C2jk3tQIiMjYWNjAxcXFwwePBgPHz4Ut0VHR6NixYpicgIAbdq0gZ6eHo4ePVrk8bKzs5GamqoxERGVFYZ5uZj6xxJM/WMJkJ0tdTiylZ2djaFDh2Lo0KHI5nmSjFz6odQTFF9fX6xevRoRERH45ptvEBUVBT8/P+Tn5wMAEhMTYWNjo7GPvr4+LC0tkZiYWOQxZ82aBQsLC3FycHAo7bCJiHRGocpH79O70Pv0LiAvT+pwZCsvLw9hYWEICwtDHs+TZOTSD6X+ssCePXuK83Xr1oW7uztq1qyJyMhItG7d+pWOOWHCBISEhIjLqampTFKIiIjeYjp/zLhGjRqwsrLC1atXAQC2trZISkrSaJOXl4fk5GSt960YGRnB3NxcYyIiIqK3l84TlNu3b+Phw4ews7MDAHh5eeHx48c4efKk2Oavv/6CSqVC48aNdR0OERERlQElHuJJS0sTr4YAQFxcHGJiYmBpaQlLS0tMmTIFXbt2ha2tLa5du4YvvvgCtWrVQtu2bQEArq6u8PX1xeeff46lS5ciNzcXwcHB6NmzJ5/gISIiIgCvcAXlxIkTqF+/PurXrw8ACAkJQf369TFp0iQoFAqcPXsWHTt2hLOzM4KCgtCwYUMcPHhQoxbK2rVrUbt2bbRu3Rrt2rVD8+bNsWzZstL7VkRERFSmlfgKio+PD9Rqtdbte/fufekxLC0tWZSNiIiItCr1p3iIiEhTloEhmg/6CQBwSKmUOBr5UiqViIuLE+dJGnLpByYoREQ6phb0cNuiSsGCHt/Rqo2enh4cHR2lDqPck0s/8E8KERERyQ6voBAR6ZhBfi7GHFhTsJDzIWBoKG1AMpWTk4OvvvoKADBjxgwY8jxJQi79IKhfdMerTKWmpsLCwgIpKSks2kYl5jh+l9QhUDmjzMlC7LxuBQtpaYCpqbQByVR6ejoqVKgAoKCkhSnPkyR02Q8l+f3mEA8RERHJDhMUIiIikh0mKERERCQ7TFCIiIhIdpigEBERkewwQSEiIiLZYR0UIiIdyzIwxIf9FgMA9rGEu1ZKpRLnzp0T50kacukHJihERDqmFvRwxbpawQJL3Wulp6eHOnXqSB1GuSeXfuCfFCIiIpIdXkEhItIxg/xcDI3eVLDAUvda5eTkYObMmQCAL7/8kqXuJSKXfmCpeyp3WOqe3jSWui8elrqXB5a6JyIiItKCCQoRERHJDhMUIiIikh0mKERERCQ7TFCIiIhIdpigEBERkeywDgoRkY5l6xugY+/vAQA7jI0ljka+jI2NcezYMXGepCGXfmCCQkSkYyo9Bc7aORcsKBTSBiNjCoUCnp6eUodR7smlHzjEQ0RERLLDKyhERDpmkJ+Lvid2FCyw1L1WOTk5+OGHHwAAI0aMYKl7icilH1jqnsodlrqnN42l7ouHpe7lgaXuiYiIiLRggkJERESywwSFiIiIZIcJChEREckOExQiIiKSHSYoREREJDusg0JEpGPZ+gbo2WsmAGADS7hrZWxsjP3794vzJA259AMTFCIiHVPpKXCkqnvBAkvda6VQKODj4yN1GOWeXPqBQzxEREQkO7yCQkSkY/r5eeh1Zk/BQu5HgIGBtAHJVG5uLpYtWwYAGDBgAAx4niQhl35gqXsqd1jqnt40lrovHpa6lweWuiciIiLSggkKERERyQ4TFCIiIpIdJihEREQkO0xQiIiISHaYoBAREZHssA4KEZGO5egboG+3yQCAlUZGEkcjX0ZGRti5c6c4T9KQSz8wQSEi0rF8PQX21/QsWNDnX7va6Ovro3379lKHUe7JpR84xENERESyw1SeiEjH9PPz0PlCZMECS91rlZubi7Vr1wIAAgICWOpeInLpB5a6p3KHpe7pTWOp++JhqXt5YKl7IiIiIi2YoBAREZHsMEEhIiIi2WGCQkRERLLDBIWIiIhkhwkKERERyQ7roBAR6ViOvgGGdBoPAAhjCXetjIyMsGnTJnGepCGXfmCCQkSkY/l6Cuyu3bxggaXutdLX10f37t2lDqPck0s/cIiHiIiIZIepPBGRjilU+Wh7ObpgIa8tr6JokZeXh19//RUA8Mknn0Cf50kScukH9j4RkY4Z5uUibPvsgoXsr5mgaJGdnQ1/f38ABSXWmaBIQy79wCEeIiIikh0mKERERCQ7TFCIiIhIdpigEBERkewwQSEiIiLZYYJCREREssNnuIiIdCxXoY8x7UYCAOYaGkobjIwZGhpi5cqV4jxJQy79IKjVarVkn/6KUlNTYWFhgZSUFJibm0sdDpUxjuN3SR0ClWPxs9tLHQKRZEry+80hHiIiIpKdEicoBw4cwMcffwx7e3sIgoBt27ZpbFer1Zg0aRLs7OygVCrRpk0bXLlyRaNNcnIyAgICYG5ujooVKyIoKAhpaWmv9UWIiORKocpHq2vH0eracSAvT+pwZCsvLw+7du3Crl27kMfzJBm59EOJE5T09HR4eHhg8eLFRW6fM2cOFixYgKVLl+Lo0aMwNTVF27ZtkZWVJbYJCAjA+fPnsW/fPuzcuRMHDhzAgAEDXv1bEBHJmGFeLlZumYKVW6YA2dlShyNb2dnZ6NChAzp06IBsnifJyKUfSnyTrJ+fH/z8/IrcplarMX/+fHz99dfo1KkTAGD16tWoUqUKtm3bhp49eyI2NhZ79uzB8ePH0ahRIwDAwoUL0a5dO8ydOxf29vav8XWIiIjobVCq96DExcUhMTERbdq0EddZWFigcePGiI4ueJNndHQ0KlasKCYnANCmTRvo6enh6NGjRR43OzsbqampGhMRERG9vUo1QUlMTAQAVKlSRWN9lSpVxG2JiYmwsbHR2K6vrw9LS0uxzfNmzZoFCwsLcXJwcCjNsImIiEhmysRTPBMmTEBKSoo43bp1S+qQiIiISIdKNUGxtbUFANy7d09j/b1798Rttra2SEpK0tiel5eH5ORksc3zjIyMYG5urjERERHR26tUE5Tq1avD1tYWERER4rrU1FQcPXoUXl5eAAAvLy88fvwYJ0+eFNv89ddfUKlUaNy4cWmGQ0RERGVUiZ/iSUtLw9WrV8XluLg4xMTEwNLSElWrVsXIkSMxffp0ODk5oXr16pg4cSLs7e3RuXNnAICrqyt8fX3x+eefY+nSpcjNzUVwcDB69uzJJ3iI6K2Uq9DHxA8HAQCmsYS7VoaGhli0aJE4T9KQSz+UuNR9ZGQkWrVqVWh9nz59EB4eDrVajcmTJ2PZsmV4/PgxmjdvjrCwMDg7O4ttk5OTERwcjN9++w16enro2rUrFixYgAoVKhQrBpa6p9fBUvckJZa6p/KsJL/ffBcPlTtMUEhKTFCoPCvJ7zffZkxEpGN6qny8f/t8wUK+L6BQSBuQTOXn5+PgwYMAgBYtWkDB8yQJufQDExQiIh0zysvFhvVfFiz8bzhgaiptQDKVlZUl3kKQlpYGU54nScilH8pEHRQiIiIqX5igEBERkewwQSEiIiLZYYJCREREssMEhYiIiGSHT/FQqWOdESIiel1MUIiIdCxPocBMn74AgC8NDCSORr4MDAwwZ84ccZ6kIZd+YCVZKnW8gkKkHSvJUnlWkt9v3oNCREREssMhHiIiHdNT5eO9e9cKFljqXqv8/HycOnUKANCgQQOWupeIXPqBCQoRkY4Z5eVix+qQgoWwASx1r0VWVhbef/99ACx1LyW59AOHeIiIiEh2mKAQERGR7DBBISIiItlhgkJERESywwSFiIiIZIcJChEREckOHzMmItKxPIUC85v1AgCMZAl3rQwMDDB58mRxnqQhl35gqXsqdSx1T6QdS91TecZS90RERFSmcYiHiEjHBLUKtR7cKlhQqQA9/tuwKCqVCrGxsQAAV1dX6PE8SUIu/cAEhYhIx4xzc7BvxdCChQV9WOpei8zMTLz33nsAWOpeSnLpB6anREREJDtMUIiIiEh2mKAQERGR7DBBISIiItlhgkJERESywwSFiIiIZIePGRMR6VieQoEf3+8CABjIEu5aGRgYYMyYMeI8SUMu/cBS91TqWOqeSDuWuqfyjKXuiYiIqEzjEA8RkY4JahXeSb1fsMBS91qpVCrcvHkTAFC1alWWupeIXPqBCQoRkY4Z5+bg0NIgAICr0gKZhsYSR6SdlENQmZmZqF69OgCWupeSXPqB6SkRERHJDhMUIiIikh0mKERERCQ7TFCIiIhIdpigEBERkewwQSEiIiLZ4WPGREQ6lq+nwOr67cV5Kpq+vj6GDBkizpM05NIPLHVPpY6l7onKLpbiJ11iqXsiIiIq03gNjYhI19RqWGamAgCSleaAIEgckDyp1Wo8ePAAAGBlZQWB50kScukHJihERDqmzM3GqYUBAADXUVtkXepeShkZGbCxsQHAUvdSkks/cIiHiIiIZIcJChEREckOExQiIiKSHSYoREREJDtMUIiIiEh2mKAQERGR7PAxYyIiHcvXU2DLe63FeSqavr4++vTpI86TNOTSDyx1T6WOpe6Jyi6WuiddYql7IiIiKtN4DY2ISNfUaihzswEAmQZGLHWvhVqtRkZGBgDAxMSEpe4lIpd+4BUUIiIdU+ZmI3ZeN8TO6yYmKlRYRkYGKlSogAoVKog/kPTmyaUfmKAQERGR7DBBISIiItlhgkJERESywwSFiIiIZIcJChEREckOExQiIiKSHdZBISLSMZWeHna5NBPnqWgKhQLdunUT50kacukHlrqnUsdS90RlF0vdky6x1D0RERGVaUxQiIiISHaYoBAR6ZgyJwvx33RA/DcdoMzJkjoc2UpPT4cgCBAEAenp6VKHU27JpR+YoBAREZHslHqCEhoaKmZeT6fatWuL27OysjB06FBUrlwZFSpUQNeuXXHv3r3SDoOIiIjKMJ1cQalTpw4SEhLE6dChQ+K2UaNG4bfffsPmzZsRFRWFu3fvokuXLroIg4iIiMoondRB0dfXh62tbaH1KSkp+Omnn7Bu3Tp88MEHAICVK1fC1dUVR44cQZMmTYo8XnZ2NrKz/3tFeWpqqi7CJiIiIpnQyRWUK1euwN7eHjVq1EBAQABu3rwJADh58iRyc3PRpk0bsW3t2rVRtWpVREdHaz3erFmzYGFhIU4ODg66CJuIiIhkotQTlMaNGyM8PBx79uzBkiVLEBcXhxYtWuDJkydITEyEoaEhKlasqLFPlSpVkJiYqPWYEyZMQEpKijjdunWrtMMmIiIiGSn1IR4/Pz9x3t3dHY0bN0a1atWwadMmKJXKVzqmkZERjIyMSitEIqI3SqWnh79qNBLnqWgKhQLt2rUT50kacukHnb+Lp2LFinB2dsbVq1fx4YcfIicnB48fP9a4inLv3r0i71khInobZOsbol/3UKnDkD1jY2Ps2sVXZUhNLv2g81Q+LS0N165dg52dHRo2bAgDAwNERESI2y9duoSbN2/Cy8tL16EQERFRGVHqV1DGjBmDjz/+GNWqVcPdu3cxefJkKBQK9OrVCxYWFggKCkJISAgsLS1hbm6OYcOGwcvLS+sTPERERFT+lHqCcvv2bfTq1QsPHz6EtbU1mjdvjiNHjsDa2hoAMG/ePOjp6aFr167Izs5G27ZtERYWVtphEBHJhjInCycXBQAAGgavRaahscQRyVN6ejpsbGwAAElJSTA1NZU4ovJJLv1Q6gnKhg0bXrjd2NgYixcvxuLFi0v7o8sNx/HSjw0SUcmY5Ga/vBEhIyND6hAI8ugH3k5OREREssMEhYiIiGSHCQoRERHJDhMUIiIikh0mKERERCQ7Oq8kS0RU3qkEAUcc3hPnqWh6enrw9vYW50kacukHJihERDqWbWCEnp/OljoM2VMqlYiMjJQ6jHJPLv3AFJWIiIhkhwkKERERyQ4TFCIiHVPmZOHkgk9xcsGnUOZkSR2ObKWnp8Pa2hrW1tZIT0+XOpxySy79wHtQiIjegMqZqVKHUCY8ePBA6hAI8ugHXkEhIiIi2WGCQkRERLLDBIWIiIhkhwkKERERyQ4TFCIiIpIdPsVDRKRjKkHAGVsncZ6Kpqenh0aNGonzJA259AMTFCIiHcs2MEKnPvOkDkP2lEoljh8/LnUY5Z5c+oEJChERiRzH75I6hBeKn91e6hDoDeE1NCIiIpIdJihERDpmnJuFQ0v64dCSfjDOZal7bVS5Wbi9pB9uL+kHFc+TZDIyMuDo6AhHR0dkZGRIFgeHeIiIdExQA++mJonzpIUayP/3PIHnSTJqtRo3btwQ56XCKyhEREQkO0xQiIiISHaYoBAREZHsMEEhIiIi2WGCQkRERLLDp3iIiHRMLQCXK1cV50kLATD49zyB50kygiDAzc1NnJcKExQiIh3LMjDGR/3DpA5D9vQMjGHP8yQ5ExMTnD9/XuowOMRDRERE8sMEhYiIiGSHCQoRkY4Z52bhj+VD8MfyISx1/wKq3CzcXT4Ed5cPYal7CWVkZKBOnTqoU6cOS90TEb3NBDXg/PCmOE9aqIHcf88TS91LR61W48KFC+K8VHgFhYiIiGSHCQoRERHJDhMUIiIikh0mKERERCQ7TFCIiIhIdvgUDxGRjqkF4La5jThPWgiA4t/zxFL30hEEAdWqVRPnpcIEhYhIx7IMjNF88Aqpw5A9PQNjvMvzJDkTExPEx8dLHQaHeIiIiEh+mKAQERGR7DBBISLSMaPcbGxfNQrbV42CUW621OHIlio3GwmrRiFh1SioeJ4kk5mZCU9PT3h6eiIzM1OyOHgPChGRjump1fBIvCLOkxZqNXL+PU/geZKMSqXCiRMnxHmp8AoKERERyQ4TFCIiIpIdJihEREQkO0xQiIiISHaYoBAREZHs8CkeIqI34KHSXOoQygQ9nidZsLKykjoEJihERLqWaWiMhsPXSR2G7OkZGsOB50lypqamuH//vtRhcIiHiIiI5IcJChEREckOExQiIh0zys3GhnXjsWHdeJa6fwFVbjYS141H4rrxLHUvoczMTPj4+MDHx4el7omI3mZ6ajWa3DonzpMWajWy/z1PLHUvHZVKhaioKHFeKryCQkRERLLDBIWIiIhkhwkKERERyQ4TFCIiIpIdJihEREQkO3yKpwiO43dJHQIRvWUyDIykDqFMEHieZMHExETqEJigEBHpWqahMdxCfpE6DNnTMzRGVZ4nyZmamiI9PV3qMDjEQ0RERPLDKyhERFRmyH0IPn52e6lDeGvwCgoRkY4Z5eVgxeZQrNgcCqO8HKnDkS11Xg6SNociaXMo1DxPksnKykL79u3Rvn17ZGVlSRYHr6AQEemYnkqFD66fEOepaGqVCpn/nie1SgVB4njKq/z8fOzevVuclwqvoBAREZHsMEEhIiIi2ZE0QVm8eDEcHR1hbGyMxo0b49ixY1KGQ0RERDIhWYKyceNGhISEYPLkyTh16hQ8PDzQtm1bJCUlSRUSERERyYRkCcr333+Pzz//HH379oWbmxuWLl0KExMTrFixQqqQiIiISCYkeYonJycHJ0+exIQJE8R1enp6aNOmDaKjowu1z87ORnZ2trickpICAEhNTdVJfKrsDJ0cl4jKp/ycLDz92yo/OwMqNZ/kKYoq579HWlXZGUAZPE+6+l16k56tIpuamlqqT/I8PT9qtfqlbSVJUB48eID8/HxUqVJFY32VKlVw8eLFQu1nzZqFKVOmFFrv4OCgsxiJiEqTxdOZsN5ShlFm3Cmj58livtQRlC57e3udHPfJkyewsLB4YZsyUQdlwoQJCAkJEZdVKhWSk5NRuXJlCAKflCcqbampqXBwcMCtW7dgbm4udThE9JZQq9V48uRJsRIfSRIUKysrKBQK3Lt3T2P9vXv3YGtrW6i9kZERjIw033BZsWJFXYZIRADMzc2ZoBBRqXrZlZOnJLlJ1tDQEA0bNkRERIS4TqVSISIiAl5eXlKERERERDIi2RBPSEgI+vTpg0aNGuH999/H/PnzkZ6ejr59+0oVEhEREcmEZAlKjx49cP/+fUyaNAmJiYmoV68e9uzZU+jGWSJ684yMjDB58uRCQ6tERG+KoC7Osz5EREREbxDfxUNERESywwSFiIiIZIcJChEREckOExQiIiKSHSYoRGWIIAjYtm2b1GEgPDz8pcUSQ0NDUa9evRe2iY+PhyAIiImJKbXYiOjtwASFSGYCAwPRuXPnIrclJCTAz8+v1D5r/PjxqF27tsa6ixcvQhAEBAYGaqwPDw+HkZERMjMz0aNHD1y+fLlEn/Wi70VE9DwmKERliK2tbanWJmnVqhUuXbqExMREcd3+/fvh4OCAyMhIjbb79+9HkyZNoFQqoVQqYWNjU2pxEBE9jwkKURny7BDP0+GRrVu3olWrVjAxMYGHhweio6M19jl06BBatGgBpVIJBwcHDB8+XHydevPmzWFgYKCRjERGRmLo0KFITk5GfHy8xvpWrVoBKHqIZ/bs2ahSpQrMzMwQFBSErKwscVtoaChWrVqF7du3QxAECIKg8ZnXr19/4XcgovKHCQpRGffVV19hzJgxiImJgbOzM3r16oW8vDwAwLVr1+Dr64uuXbvi7Nmz2LhxIw4dOoTg4GAAgKmpKTw9PbF//37xeJGRkWjdujWaNWsmrr9+/Tpu3rwpJijP27RpE0JDQzFz5kycOHECdnZ2CAsLE7ePGTMG/v7+8PX1RUJCAhISEtC0adNifQciKp+YoBCVcWPGjEH79u3h7OyMKVOm4MaNG7h69SoAYNasWQgICMDIkSPh5OSEpk2bYsGCBVi9erV4haNVq1bi1YwLFy4gKysL9evXR8uWLcX1kZGRMDY2RpMmTYqMYf78+QgKCkJQUBBcXFwwffp0uLm5idsrVKgApVIJIyMj2NrawtbWFoaGhsX6DkRUPjFBISrj3N3dxXk7OzsAQFJSEgDgzJkzCA8PR4UKFcSpbdu2UKlUiIuLAwD4+Pjg8uXLSEhIQGRkJJo3bw6FQgFvb2+NBKVp06Za73+JjY1F48aNNdaV5M3kL/oORFQ+SfayQCIqHQYGBuK8IAgAAJVKBQBIS0vDwIEDMXz48EL7Va1aFQDQrFkzGBoaYv/+/di/fz+8vb0BAJ6ennjw4AGuX7+OyMhIDBw4UJLvQETlExMUordYgwYNcOHCBdSqVUtrG6VSicaNGyMyMhJRUVEYO3YsgIKkoUmTJvjpp59w69YtrfefAICrqyuOHj2K3r17i+uOHDmi0cbQ0BD5+fmv+Y2IqLxggkIkQykpKYWKl1WuXLnExxk3bhyaNGmC4OBg9O/fH6amprhw4QL27duHRYsWie1atWqFefPmAShIap7y9vbG3LlzxZtptRkxYgQCAwPRqFEjNGvWDGvXrsX58+dRo0YNsY2joyP27t2LS5cuoXLlyrCwsCjx9yGi8oP3oBDJUGRkJOrXr68xTZkypcTHcXd3R1RUFC5fvowWLVqgfv36mDRpEuzt7TXatWrVCk+ePEGzZs2gr//fv1u8vb3x5MkT8XFkbXr06IGJEyfiiy++QMOGDXHjxg0MHjxYo83nn38OFxcXNGrUCNbW1jh8+HCJvw8RlR+CWq1WSx0EERER0bN4BYWIiIhkhwkKERERyQ4TFCIiIpIdJihEREQkO0xQiIiISHaYoBAREZHsMEEhIiIi2WGCQkRERLLDBIWIiIhkhwkKERERyQ4TFCIiIpKd/wfyGQRToI4kugAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# First, let's load the original data from their paper.\n",
        "data = loadmat('/content/drive/MyDrive/Colab Notebooks/Master_et_al_2005.mat')\n",
        "nsubjs = 30\n",
        "r_stat = 0.23\n",
        "\n",
        "sigmas = data['sigmas']\n",
        "grades = data['grades']\n",
        "\n",
        "# At the core of a permutation test is the creation of a null distribution\n",
        "# of results. I.e. How many times do we achieve the r-value Master et al.\n",
        "# published, purely due to chance?\n",
        "\n",
        "# To do this we permute (randomize) our variables of interest such that we\n",
        "# no there is not any relationship between them.\n",
        "# Then we compare the published r-value to our null distribution to\n",
        "# compute a p-statistic (a percentage describing how many times our null\n",
        "# distribution achieves the value originally obtained).\n",
        "\n",
        "n_perms = 1000\n",
        "r_perms = np.zeros(n_perms)\n",
        "\n",
        "# --------------- YOUR CODE HERE ---------------- #\n",
        "\n",
        "\n",
        "\n",
        "# --------------- YOUR CODE HERE ---------------- #\n",
        "\n",
        "p_value = np.sum(r_perms >= r_stat)/n_perms\n",
        "\n",
        "# Find confidence intervals\n",
        "ratio_boot_sorted = np.sort(ratio_boot) # Sort ratio_boots from lowest to highest\n",
        "lower_bound_ind = int(np.ceil(0.025 * n_boot)) # Find index of 2.5% value\n",
        "upper_bound_ind = int(np.floor(0.975 * n_boot)) # Find index of 97.5% value\n",
        "\n",
        "lower_bound = ratio_boot_sorted[lower_bound_ind]\n",
        "upper_bound = ratio_boot_sorted[upper_bound_ind]\n",
        "\n",
        "# Let's see whether we can conclude with confidence that there is a\n",
        "# relationship between attentional focus & math tools grades.\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "plt.hist(r_perms)\n",
        "lowery,uppery = ax.get_ylim()  # return the current ylim\n",
        "plt.plot([r_stat,r_stat],[lowery,uppery],'--r','LineWidth',2)\n",
        "plt.plot([lower_bound,lower_bound],[lowery,uppery],'--k')\n",
        "plt.plot([upper_bound,upper_bound],[lowery,uppery],'--k')\n",
        "plt.legend(['Null dist','Original r-statistic','Lower 95% conf. bound','Upper 95% conf. bound'])\n",
        "plt.title('p = ' + str(p_value))"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "math_tools",
      "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.10.6"
    },
    "vscode": {
      "interpreter": {
        "hash": "ab2b88acfe38f77fc88fd37d062c2dfd9d065fa5edb69e3c34007729e94448e5"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}