{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8545ca88",
   "metadata": {},
   "source": [
    "### Exercise 43"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb891b6d",
   "metadata": {},
   "source": [
    "#### Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "accffd3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.linalg import det, inv\n",
    "from scipy.stats import norm\n",
    "from sklearn import datasets\n",
    "from sklearn.decomposition import PCA\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LinearSegmentedColormap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a13b8fa",
   "metadata": {},
   "source": [
    "#### Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b8a54af2",
   "metadata": {},
   "outputs": [],
   "source": [
    "iris = datasets.load_iris()\n",
    "\n",
    "X = iris.data\n",
    "y = iris.target\n",
    "\n",
    "X = X[y < 2]\n",
    "y = y[y < 2]\n",
    "\n",
    "n = len(y)\n",
    "X = np.column_stack((np.ones(n), X))\n",
    "p = len(X[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da58494a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "36d2b377",
   "metadata": {},
   "source": [
    "#### Principal component analysis (PCA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d40138ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.        ,  0.32327447, -0.17121498,  0.86913407,  0.33284383],\n",
       "       [-0.        ,  0.65835816,  0.74710188, -0.08837738, -0.02434645]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pca = PCA(n_components=2)\n",
    "pca_fit = pca.fit(X)\n",
    "\n",
    "M = pca.components_\n",
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08766f84",
   "metadata": {},
   "source": [
    "### Find the components directly via the spectral theorem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "03dc9f07",
   "metadata": {},
   "outputs": [],
   "source": [
    "SS = (1 / (n - 1)) * (X - np.mean(X, axis=0)).T @ (X - np.mean(X, axis=0))\n",
    "MSMT = M @ SS @ M.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8a446a3d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 0.01046467, 0.05123085, 0.22795013, 2.77191092])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lambdas, U = np.linalg.eigh(SS)\n",
    "lambdas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e5fce7cb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.        ,  0.        ,  0.        , -0.        ],\n",
       "       [ 0.        ,  0.09641117,  0.67287668,  0.65835816, -0.32327447],\n",
       "       [ 0.        , -0.10310285, -0.63395111,  0.74710188,  0.17121498],\n",
       "       [ 0.        , -0.40254057, -0.27341635, -0.08837738, -0.86913407],\n",
       "       [ 0.        ,  0.9044533 , -0.26568112, -0.02434645, -0.33284383]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e8b90fea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 0.01046467, 0.05123085, 0.22795013, 2.77191092])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lambdas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "144f244b",
   "metadata": {},
   "source": [
    "#### Plot data after PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e3a22672",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAniUlEQVR4nO3df5RkdXnn8ffTM4NjB50euhWEoatQMCsSacJAS7IcNZITdBEWAkbsQ/SgO0iSqZrGk0XSuzM9OenEJMrQIzmb9GrCSaYFxV0XYthlg8aFcwINA7byK0TEbhiiBkZmMBkhzNSzf1RVc7v61q3qqlt1b1V9XufcM31v3Xv7W8Nwn/v9Pt8f5u6IiIhU05d0AUREJN0UKEREJJIChYiIRFKgEBGRSAoUIiISaW3SBYjb0NCQZ7PZpIshItJRHnzwwefd/Q1hn3VdoMhms+zduzfpYoiIdBQzW6z2mZqeREQkkgKFiIhEUqAQEZFIChQiIhJJgaLDVc7Vpbm7RCRuChQdbHJykvHx8aXg4O6Mj48zOTmZbMFEpKsoUHQod+fAgQNMT08vBYvx8XGmp6c5cOCAahYiEpuuG0fRK8yMXbt2ATA9Pc309DQA+XyeXbt2YWZJFk9Euoh125vn5s2bvZcG3Lk7fX2vVgwLhYKChIismpk96O6bwz5T01MHKzc3BQVzFiIicVCg6FDBnEQ+n6dQKJDP55flLERE4qAcRYcyMwYGBpblJMo5i4GBATU/iUhslKPocO6+LChU7ouI1EM5ii5WGRQUJEQkbgoUIiISSYFCREQiKVCIiEgkBQoREYmkQCEiIpEUKEREJJIChYiIRFKgSAEtPiQiaaZAkTAtPiQiaadAkSAtPiQinUCTAiZIiw+JSCdIdFJAM/tz4ALgn939tJDPDZgG3g8cAj7q7g9F3bMTJwXU4kMikrQ0Twp4E3B+xOfvA04pbVuA/9aGMrVVrcWH1PwkIklLtOnJ3e82s2zEKRcBf+nFp+V9ZjZgZm9y9x+0p4StVbn40IYNG7j99tuXmqCuv/56rrnmGgYGBpTcFpHEJF2jqOUE4JnA/r7Ssa4QXHzo+uuv5+DBg8zPzzMyMsKGDRu45pprlNgWkcR1RTLbzLZQbJpieHg44dKszuTk5NJiQ8HE9vz8PKDEtogkL/EV7kpNT1+rksz+M+Cb7n5zaf8J4N1RTU+dmMwOUmJbRJKQ5mR2LbcDv25F7wQOdkt+IkytxLaISBISDRRmdjNwL/CzZrbPzD5mZp8ws0+UTrkDeAp4EvjvwG8kVNSWq0xsFwoF8vn8ssF4IiJJSLrX0+U1PnfgN9tUnEQFE9vlnEQ5ZzEwMKDmJ5EuMzs7y8TEBE8//TTDw8NMTU0xNjaWdLFCJZ6jiFs35CiCQaFyX0Q63+zsLFu2bOHQoUNLx/r7+5mZmUksWETlKBQoRETaLJvNsri4uOJ4JpNhYWGh/QWis5PZIiJd5+mnn17V8aQpUIiItFm18V5pHQemQCEi0mZTU1P09/cvO9bf38/U1FRCJYqmQCEi0mZjY2PMzMyQyWQwMzKZTKKJ7FqUzBYRESWzRUSkcQoUPayyNtlttUsRiYcCRY+anJxcsUDS+Pi41r0QkRUUKHqQu3PgwIFl80iV55nS2hciUqkr1qOQ1alc+6K8op7WvhCRMKpR9KhgsChTkBBJr9nZWbLZLH19fWSzWWZnZ9v2uxUoepTWvhDpHOVJBBcXF3F3FhcX2bJlS9uChQJFSrSzB5LWvhDpLBMTE8tmmgU4dOgQExMTbfn9ylGkwOTkJAcOHFjWFDQ+Ps7AwAA7duyIvTlIa1+IdJakJxHUyOyEBd/uR0dHGR0dBWD37t3kcjkANm7c2JJuq1r7QqQztGNa8qiR2apRJKz8Nu/u7N69m7m5OYClILF7927y+XxLHuKV91OQEEmnqamp0IWO2jWJoGoUKeHu9PWtTBnF1WVVtQeRztbqpVM111PKhfVAKosjSGgUtkjjkuyWGjQ2NsbCwgKFQoGFhYW2zjSrQJGwYI4il8stNTmVbdu2raleSBqFLdK4pLulpoWanlJgcnKSF154AViexJ6bm2Nubq7p5qdgcCjTKGyR2tK4tnWrRDU9KVCkhLuzc+fOqt1km20mqsyBFAoFBQmRGvr6+kJr3WZGoVBY9f1anWdoRlSgwN27ajvzzDO9kxUKhcj9Ru+Zz+cdWNry+Xws9xbpZplMZtn/N+Utk8ms+l579uzx/v7+Zffp7+/3PXv2xF/wBgB7vcpzVTmKlIm7y6prFLZIw+Jc2zrp0dXN0DiKLqdR2CKNKzcLxdFclPTo6mYoR9EjXOMoRBKV9sS4xlGIRmGLJCzOZqx2U6DoQJW1wG6rFYp0o7GxMWZmZshkMpgZmUyGmZmZ1PR6iqKmpw4TnGnWzJaS1XF0oRWR3qWmpy7hGmUtkoi0TONRTcvLV63fbKdunT6OohaNiRBpnz179vjg4OCKcRTVxj/s2bPHM5mMm5lnMpm2jJGIa3wGEeMoEn2oA+cDTwBPAp8K+fyjwHPAfGn7eK17dnugcC8Gi+A/CgUJkfiFPYCJGHSX1IC6uAYFRgWKxJqezGwN8CfA+4BTgcvN7NSQU7/k7iOl7fNtLWSCPCRhXd601rVI64UNkAuqHP+Q1IC6dozPSHLA3dnAk+7+FICZ3QJcBDyWYJlSISxhfc455wAwOjq6YuLA8mR/muRPJD61HrTDw8N1nd/qAXXDw8Oh4zMqy9eMJJPZJwDPBPb3lY5V+lUz+46ZfcXMTgy7kZltMbO9Zrb3ueeea0VZ28ZDEtbbtm1bmkl2bm5u2ep3o6Oj5HI5jbIWiVnUgzZs/EO181f7wF5tYrot4zOqtUm1egMuBT4f2L8CuLHinEHgNaWfrwK+Ueu+3ZCjCEtY53I5z+VyoUls5ShE4lctRzE4OFg1kd1sjqLRe8SRRCeNyWzgHODOwP51wHUR568BDta6bzcECvfwhLWS2CLttdoHcLMP7Dhnq12tqECRZI7iAeAUMzsJeBb4EPDh4Alm9iZ3/0Fp90Lg8fYWMRkekrDetm3bivPGx8eVlxBpobGxsbpGTleuM/FXf/VX3TVxYLUI0o4NeD/wj8D3gInSsd8FLiz9/AfAo8C3gb8D/l2te3Z6jSLY7DQ6Ouq5XM63bt269GYxMjLi27dvXzpHYyhEmtNsLSDObrFprVEkGihasXV6oHB337Fjx7KcRC6X8+OPP37pH2M+n/fDhw/7yMiIv+td70q6uCIdK46HfLcsbqRA0YHKOYnKpPbIyIgfOXIkthpFK1bUE+kUcTzkzSz0HmbWUJmSGN3tHh0oNClgyrkvX+s6KLgYUSM0waD0ujjWxE77OhP10qSAHar84K6mmSDhrgkGReIY+9DJ60zUrVpVo1O3bmp6CjYvHTlyxEdGRlaMrWi22UkTDEoviysnkFRzUZxQ01NnKjcNXX/99fzCL/wCc3NzjIyMcOGFF3LgwIGlkdn33ntvUzWLYNNWoVBQd1vpKZVdWxtdE7vTqempQ01OTq5oXjr33HNjyx94SNOWJhiUblc5RQbAwsIChUKBhYWFngwSNVWranTq1i1NT5UKhcKKKTzCmp7q7cVU2bQVti/SbZLsfpp2qOmpO3iNZqJavZjcfdn5O3bs4ODBg+r1JD2jW3ootULDTU9m9itm9jEzy1YcvzLG8kkdyg/xoGAzkXt0L6YdO3asOP/gwYNs2LBhKXiYGbt27VKQkK6V2iky0q5aVQP4feBu4AaKU2xsDXz2ULXrkt66semp3maiar2YwgboqZlJelG7psjoxF5QNDIyG3gYWFv6eQC4A9hV2v9WteuS3roxULi/Oq1HMCjkcjnfsWPHsvOqzTCrrrAi7clRdGoeJCpQRDU9rXX3w6VaxwHgA8DrzexW4KjG6i/SSh7RPFVuVgrSzLPSa8bGxpiZmSGTyWBmZDIZZmZmYu3plNSSqC1VLYIAXwPeFXL894BCteuS3rqxRlFP01Otc4LNT6hGIdIycc/91C402PT0WuC1VT47odp1SW/dGCjc62s62rFjR2jOImxacuUoRFojyanCmxEVKNQ9NmHur3ZZLf+3CO4Hm4bca4+iDrvGzDQBoEibzM7OsmXLlmXNT/39/bE3ccUtqntskivc9bzgw3vnzp288MILAEsP723btrFx40Z27NgBrFzlLmyFu8rAUd6vHEdRzlkoRyESr3Iw6KppQapVNTp165Smp2DzT3CRIkprTpx99tkO+NatW33r1q1L+6Ojo2o6ElmFTuyqmgSaWTPbzN4C7HP3l83s3cA7gL/0Yk8oaVCwF9L09PSyz+bn55d+vueee5ifn+e4444DYHR0FGDp2oGBAdUKRKqobAZaXFxky5YtAJ39ht9u1SJIeQPmKTZRnUxxfes/Bu6odV1SW6fUKMoqxz1EbZVzO6kmIRItjsRyIzWSTqzF0MxSqJRGYQO/TWl0Nhpw17Bgj6Swif6iNgUGkdVptqtqI4PnunHAXc1eT2Y2R3EajwngA+7+fTN7xN1PW13dpT3S3OupnLzesGHD0ipyn/vc5zjrrLN45ZVXlpqcRkZGljU/lTW79KlIr2l2EsBGru/UiQejej3VU6M4FdgNXF7aPwm4ttZ1SW1prVEEE9DlleqOO+64Zfunn366n3766aHnlI8HB9iJSLRm3+4bqZF044C7mgsXuftj7p5z95tL+9939z9cfbzqbeXkdT6fX6ot/PCHPwSKyet8Ps9DDz3ERRddRD6f58ILLySXy/Hss88yOjrK+vXryeVyDAwMAMWusdXGP3hFLbFyX6RXNDtlRyNrasexDnfqVIsgwJdLfz4MfKdyq3Zd0ltaaxRl1ZLXYUnqyjxGOZkd1TW22ujsyskDRaQ25ShK8aDqB/Cm0p+ZsK3adUlvaQ4UYdNwlLeo8RD1zvyqVetE4hfVg6naZ73Y6+nUkGPvrnVdUltaA0VYjqLyz1rBoloNpNrvqScIiaRZmh+4nVpzqCYqUNTMUQBfNrNrrei1ZvY54A/quE4CzIyBgYGl/EM+n+fBBx9cth82eK78HypqdbvK36PpxKUblAfLLS4u4u5Lg+VmZ2dXnJfNZunr6yObza74vFW6cjrxaqpFkPIG/AxwI3Av8AhwHdBX67qktrTWKMqCuYOwP4PKixUFcxS5XM5HR0fran5CNQrpYPUMlkvyrb5TezdVQ5M1ileAn1Kcdnw98H13L8QQo3pScFK+sD/L3ItrYO/evZu5uTlyuRwAu3fvZnR0dKkHVHCmWfdX18nO5XIUCgXy+fyydbRFOkU961u34q2+3hpKV/ZuqqKe2WMfAG4DzgKGgD81s19198taWrIeVzkX1NzcHPDqoLvyObB8FtqBgYGloLJz587QOaHcw6ciF0mT4eHh0IFrwQdxPcFkNVYzN9TU1FTodOJTU1MN/e5Uq1bVKG/A5pBjV9S6Lqkt7U1Pq1UriR3Ws6ncVBU2OE/dZ6VT1NOsFPciQbXuV5lcv/rqq1ObbF8tmun1tHQivBEYLm/1XlfjnucDTwBPAp8K+fw1wJdKn88B2Vr37KZA0Ui32HrPU/dZ6QS1ej3FnaOIyjt0Wy+nSk0FCuADwHeBfwW+DxSAR2tdV8d91wDfA94MHAV8m4quuMBvAH9a+vlDwJdq3bfRQBH2UG23ykF3lTWDqIe6us9Kr4qzC21UjaJTlzitV7OB4tvAIKUZY4H3AF+odV0d9z0HuDOwfx1wXcU5dwLnlH5eCzwPxYkMq22NBIo0NMeElWF0dHRpoaKocq324V9vUBHpNVG1hm7r5VSp2UCx118NGH3ln2tdV8d9LwU+H9i/Arix4pxHgE2B/e8BQyH32gLsBfYODw+v6i8nDc0xUWWotQbFasuvGoVItGo1FNUooh/odwFHA58Dbgamgb+vdV0d940tUAS3RmoUaXh4NlOGemtEaQiKIp1KOYroB/rPAH2lpp+PADlgsNZ1ddw3NU1P7q1tjqk3/9FMGer9HWloZhOpJa1Td6S1XHFoKlC0ais9+J+iuL5FOZn99opzfpPlyewv17pv2moUjbztt7pWk4bEvUg13f7mnlapDBTFcvF+iutwfw+YKB37XeDC0s/rgVspdo+9H3hzrXuuNlC0sjmm3nvXOu/IkSMr7ivSrbo9F5BWUYGi5lKonaaRpVCDI5vNDPfiVBgDAwNVFweqV/le09PTS8fyIUuaVivD/Pw8IyMjLSmbSBr19fUR9lwyMwoFzR7UKk0thVq5AScCv73a69q1pXEcxWrGOAQdOXIkttqOmpskrSrb/QcHB3u6RpFUHoRmm56AN1Ac/HYPxWaiz9RzXRJb2kZmN5t7iCN3oQS2pFVYPmLdunV+1FFH9WSOIsn8TEOBAngdxV5Od1Ickf1ZYF+189OypSlQxJX/aLY3lLrESlpVy0cMDg4uvVUPDg764OBgV/Y0qpRkfqbRQPFT4P8B58JSLuOpauenZUtToHBv/m0+jhpFGsaJSGdpV/NHrdHOvdYDKsnR340Gim3AfcDDwO8Ab1GgaEyj+YE4awOatkPq1c6Hc6036DjfsDthDETH1SiWTihO2vc7pYDxEnAt8NZa1yW1pTFQNCOO/IJqFLIa7XxY1QpKcb1hd0rNpFo52zGdeVOBYtnJcBowBTy5muvauXVboHBvrseSchSyWq1q/qj2Rh/1ph9X0OqksRlha160I8g12vR0MvCLIcf/PXByteuS3tIeKJLopqpeT7IarXio1nqjDz4cg8nrwcFBX7duXdMPyU6e+bVdQa7RQPE14OdCjv8c8NfVrkt6S3OgSPKBrXEUUq9WNNNEPezCfl9wO+qoo5ru9dTowzYNeY12BblGA8UDEZ89XO2zpLe0Bgo1AUknifsBGfWwq/YQj/PtuZHgl5a8RtprFN+N+Ew5igY0ssBQ1J8inSLqYVctiMT99nz11Vf7mjVrHPA1a9b41Vdf3XCZ26ldAavRQHEz8J9Cjn+cOpYkTWpLc6Bwr7+barmZavv27UsTAwb3w2afjdoXSVLUwy6tNYo05TXa0QTWaKA4Fvh74JsUR2V/luIAvHuB46pdl/SW5kBRb40ieN7IyEjon8HrlKyWThDV6ykqRxHH23MjtYO01CjapaFAsXRCcY3sraXtl2qdn/SWxkBRKBRClzeNylGEBZWw4KLch3SDYLOQmfnRRx8d69tzI7WDtOQo2qXRGsV6iqOzbwSuAtZWOzdNW9oCRfBtf8eOHZ7L5TyXy/mOHTtqvvlXNlNVa67SgDppRBp69JTL0YoHcvD7lYPQamsHafk7aodGA8WXgD2lIPG/gBuqnZumLU2BIuztPpfLLdsPe5hX1kAqt3KNpPKaenIfIu7peltu17iNVjRpdZNGA8XDgZ/XAg9VOzdNW5oChfvq3/aDtQ7Ajz322GV/lnMUo6OjS/cIrluhGoXUI03t761IGlf7fmvWrOmJ2kEjGg0UD0Xtp3VLW6BwX93CReUH/ujoqG/dunUpQJx++um+detWP/744x3ws88+2wuFgh85cmRFgls5CqklTT16WhG00vT9OkWjgeII8GJp+wlwOPDzi9WuS3pLW6BoZOxE5fmnn3561eanYO+o8tra6vUktaSpRtHukeASrqFA0albmgJFMzmK4D/uI0eOrAgSlYGnHCSC9xCpJk05inJ5opLGq00qp+37dQIFigTV6vVU/rksLN9QblqqFigUFKQRndKjp9GHfju+X6f8HdZDgSJh1cZRjI6OLtuvzDcE98tNS5VBQrkI6XZpbUbqtlpLVKDoQ1rOzDAzdu3aRT6fZ/fu3fT19TE3N7fsvGuuuYb5+XmGhob47Gc/S19fHx/4wAcYGhri9a9/PWa2dO7o6CiFQoF8Ps/09DTj4+PFyC/SZZ5++ulVHW+XiYkJDh06tOzYoUOHmJiYSKhEraNA0UblYBGUy+WWAsf09DQjIyM8//zzfPKTn8TdefHFF3n++ec544wzANi4cSO5XI577713WfAZGBhYFkhEusXw8PCqjrdLWgNYS1SranTqlsamp7KwHk2VTUn1jInQBIDSS9LaxJPWJrFGoRxF8qJ6QFUGhcpeTgoE0uvSmDROawBrVFSgUNNTm5gZAwMD5PP5Fc1P5XxDLpdjenqaM888c9nnyj9Ir5mdnSWbzdLX10c2mwVgYWGBQqHAwsICY2NjyRYQGBsbY2Zmhkwmg5mRyWSYmZlJRdniZt32ANq8ebPv3bs36WJU5e5LuYTJyUleeOEFbrjhBnbu3MmPf/xj7rnnHubn58nlcgDMzc0xNze3FGCUh5BuNzs7y5YtW5Ylivv7+7v2IZwWZvagu28O/UyBIlnlv//x8XGmp6c5++yzeec73wnA7t27lwLGxo0bmZycTKqYIm2TzWZZXFxccTyTybCwsND+AvWIqECReE4h7i3JHEUzSeaoqT6Uo5BO1GheQfM0JYOIHIVqFDGZnJzkwIEDS81D7s74+DgDAwN11wTcnb6+V9NGhUJBTU3SkZppPlKNIhlRNQols2Pg7hw4cGDZwLdyU9KBAwfqSkSXrwlSEls6VTOD0aampujv7192rL+/n6mpqVjLKKtQrarRyg04Bvhb4LulPzdWOe8IMF/abq/n3kk1PTWzypyWM5Vu02zzUSOTAKat+2ynIW1NT2b2R8CP3f3TZvYpioHi2pDz/sXdj17NvZNMZnsTTUdxNF2JpEU7m4/USyoeqUtmA08Abyr9/CbgiSrn/ctq792JNYrgPaL2RTpFOwejddsI6aSQtpHZwIHAzxbcrzjvMLAXuA/4jxH321I6b+/w8HD8f4M1qOlIZKW415ioRr2k4hEVKNbGUmcJYWZ3AceFfLQsm+W+9B86TMbdnzWzNwPfMLOH3f17lSe5+wwwA8WmpyaLvmqVo66Dk/9psj7pVWNjY1WbfiqbixYXF9myZcvSdasxPDwc2syV9KSB3SSpHMUTwLvd/Qdm9ibgm+7+szWuuQn4mrt/Jeq8pHMUwaBQuS8iRXHmMJSjiEcau8feDnyk9PNHgNsqTzCzjWb2mtLPQ8AvAo+1rYQNqAwKChIi4eKcoruX5lxKSlKB4tPAL5vZd4HzSvuY2WYz+3zpnLcBe83s28DfAZ9291QHCpFuVDlB3+zsbOhnQ0NDDA0NhZ5XKe41JsbGxlI3aWBXqZa86NQtrdOMi3SiqN5LYZ+Fnbfa+0oySNs4ilZKw6SArlyFdImoXAIQ+lnledVyDrOzs0xMTPD0008zPDzM1NSUagIJ0uyxbaSBc9JN+vr6CHtGlF98aj0/zIxCodCSskm80pjM7kruzc/5JJImUbmEevIJ6qLaHVo2jqIXBcdPTE9PMz09DaBFh6RjTU1NhXY9LU/QV/lZkCby6yLVkhedusWRzG52Ko1CoaA1r6VrRI2gDn42ODjog4ODmpivQ6Fkdv2azTGUzy/XJkA1ChFJP+Uo6uRN5hiC5+fzeQqFAvl8ftn9RGS5qHEakg7KUQQ0m2OoNedTJVe3Welxcc75JK2jpqcQ7s0tSVoZANydnTt3qttsD9OYgXBa9jQ91PS0CuUHeNBqm43Cgoq6zfau8lvz4uIi7r701qwmlnjnfJIWqpbl7tStmV5PrVxXIo6FjaQzaWGd6vR3kx5E9HpSjSKgWo4hn883va5EMF9Rpp5QvUFvzdVNTU3R39+/7JjGX6RQtQjSqVsaxlFUu6dqFL1Jb83R4lrpTpqDahSrE/e6Eu7qNtvL9NYcTVOEp5+6x7aBlkrtbeUHn3o9SadS99g2ctf04yKSTuoemxJaKlVEOpEChYiIRFKgEBGRSAoUIiISSYFCREQiKVCIiEgkjaMQkZ70yiuvsG/fPl566aWki9JW69evZ9OmTaxbt67uaxQoRKQn7du3j9e97nVks9me6aru7uzfv599+/Zx0kkn1X2dmp5EpCe99NJLDA4O9kyQgOLYrcHBwVXXohQoRKRn9VKQKGvkOytQiIhIJAUKkRSanZ0lm83S19fH0NAQQ0ND9PX1kc1mtTJeQirnxUvjPHkvv/wyv/Zrv8bJJ5/M6OhobMvJKlCIpEzl0qn79+9n//79WkY1QZOTk8uWBCgvHZC29e6/8IUvsHHjRp588knGx8e59tprY7mvAoVIykxMTHDo0KGqnx86dIiJiYk2lqi3uXtL1rzfvn07N9xww9L+xMQE09PTTZX1tttu4yMf+QgAl156KV//+tdjqfmoe6xIytSzRKqWUW2f4Pox09PTSw/z4Poyjbjyyiu55JJL2LZtG4VCgVtuuYX7779/xXnnnnsuP/nJT1Yc/8xnPsN555237Nizzz7LiSeeCMDatWvZsGED+/fvZ2hoqKEyliUSKMzsMmASeBtwtruHLiBhZucD08Aa4PPu/um2FVIkIcPDwywuLtY8R9qnHCyCb/zNrnmfzWYZHBzkW9/6Fj/60Y8444wzGBwcXHHePffc0/DviEtSTU+PAJcAd1c7wczWAH8CvA84FbjczE5tT/FEkhO2dGqQllFtv3JzU1Acyxh//OMf56abbuIv/uIvuPLKK0PPOffccxkZGVmx3XXXXSvOPeGEE3jmmWcAOHz4MAcPHgwNPqtWbTHtdmzAN4HNVT47B7gzsH8dcF2te5555pkNLy4ukhZ79uzxTCbjZuaDg4M+ODjoZuaZTMb37NmTdPG6wmOPPVbXeYVCwfP5vAOez+dD9xv18ssv+1vf+lY/6aST/PDhww3fp+zGG2/0q666yt3db775Zr/ssstCzwv77sBer/JcTXOO4gTgmcD+PmA07EQz2wJsAVXJpTuMjY1pTe2UaOWa90cddRTvec97GBgYYM2aNU2X9WMf+xhXXHEFJ598Mscccwy33HJL0/eEFuYozOwu4LiQjybc/bY4f5e7zwAzUFwzO857i4hMTk4uW+O+HCyaHdldKBS47777uPXWW+MoJuvXr4/tXkEtCxTufl7tsyI9C5wY2N9UOiYi0nZxr3n/2GOPccEFF3DxxRdzyimnNHWvVktz09MDwClmdhLFAPEh4MPJFqk1gm8qYfsi0n1OPfVUnnrqqaSLUZdEej2Z2cVmto9iwvpvzOzO0vHjzewOAHc/DPwWcCfwOPBld380ifK2UqeM+BSR3pVIoHD3r7r7Jnd/jbsf6+6/Ujr+T+7+/sB5d7j7W939Le7edf0BvUUjPkVE4pTmpqeu16oRnyIicdJcTwkLBosyBQkRSRMFioSVm5uC4hjxKSK95+677+bnf/7nWbt2LV/5yldiu68CRYKCOYl8Pk+hUCCfzy/LWYhIOgTXCEnruiDDw8PcdNNNfPjD8XYQVY4iQa0c8Ski8SmvEVKe/r28LgjQ8Aj67du3c8wxx7Bt2zagOM34G9/4RvL5fMPlzGazAPT1xVsHsG57a928ebPv3Rs6GW1qaRyFSPs9/vjjvO1tb6vr3Gw2GzqjbyaTaXgVuYWFBS655BIeeughCoUCp5xyCvfff/+KSfxWM8142Uc/+lEuuOACLr300tDPw767mT3o7pvDzleNIgXiHvEpIvGqtv5HM+uCaJpxEUmFTmhX7wTVJhttdhLSuKcZbxXVKES6VCva1XvV1NTUsr9LiGddkIsvvpjt27fzyiuv8MUvfjH0HNUoRKRlwtbe1nrbjRkbG2NmZoZMJoOZkclkmJmZaTrglqcZ/+AHPxjLNOMPPPAAmzZt4tZbb+Wqq67i7W9/e9P3BNUoRLpWK9rVe1kr1giJe5rxs846i3379sVyryDVKES6VKva1SUejz32GCeffDLvfe97Nc24iCSjVe3qEg9NMy4iiWtVu3o36bZxZPVo5DurRiHSxbT2dnXr169n//79DA4O9szYJXdn//79rF+/flXXKVCISE/atGkT+/bt47nnnku6KG21fv16Nm3atKprFChEpCetW7eOk046KelidATlKEREJJIChYiIRFKgEBGRSF03zbiZPQesnA+4NYaA59v0u9JI31/fX9+/e2Tc/Q1hH3RdoGgnM9tbbf72XqDvr++v798b319NTyIiEkmBQkREIilQNGcm6QIkTN+/t+n79wjlKEREJJJqFCIiEkmBQkREIilQNMnM/tjM/sHMvmNmXzWzgaTL1E5mdpmZPWpmBTPria6CAGZ2vpk9YWZPmtmnki5PO5nZn5vZP5vZI0mXJQlmdqKZ/Z2ZPVb6t59PukytpkDRvL8FTnP3dwD/CFyXcHna7RHgEuDupAvSLma2BvgT4H3AqcDlZnZqsqVqq5uA85MuRIIOA59091OBdwK/2e3//RUomuTu/9fdD5d27wNWN39vh3P3x939iaTL0WZnA0+6+1Pu/m/ALcBFCZepbdz9buDHSZcjKe7+A3d/qPTzT4DHgROSLVVrKVDE60rgfyddCGm5E4BnAvv76PIHhYQzsyxwBjCXcFFaSutR1MHM7gKOC/lowt1vK50zQbFKOtvOsrVDPd9fpNeY2dHA/wC2ufuLSZenlRQo6uDu50V9bmYfBS4A3utdODCl1vfvQc8CJwb2N5WOSY8ws3UUg8Ssu//PpMvTamp6apKZnQ/8Z+BCdz+UdHmkLR4ATjGzk8zsKOBDwO0Jl0naxIoLbH8BeNzdr0+6PO2gQNG8G4HXAX9rZvNm9qdJF6idzOxiM9sHnAP8jZndmXSZWq3UeeG3gDspJjK/7O6PJluq9jGzm4F7gZ81s31m9rGky9RmvwhcAfxS6f/5eTN7f9KFaiVN4SEiIpFUoxARkUgKFCIiEkmBQkREIilQiIhIJAUKERGJpEAhUmJmR0pdHR8xs1vNrL90/Dgzu8XMvmdmD5rZHWb21sB128zsJTPbEFM5jjezr6zi/Ckze8bM/iWO3y9SSYFC5FU/dfcRdz8N+DfgE6XBVV8Fvunub3H3MynOEHxs4LrLKQ7CuySOQrj7P7n7pau45K8pTlQo0hIKFCLh7gFOBt4DvOLuSwMp3f3b7n4PgJm9BTga+C8UA8YKZna0mX3dzB4ys4fN7KLS8bNK65isN7OfKa1tcJqZZctrPZjZ283s/lJN5ztmdkrl/d39Pnf/Qdx/ASJlmutJpIKZraW41sT/AU4DHow4/UMUpxm/h+JI5WPd/UcV57wEXOzuL5rZEHCfmd3u7g+Y2e3A7wGvBfa4+yOlGUnLPgFMu/tsabqQNXF8R5HVUI1C5FWvNbN5YC/wNMX5fGq5HLjF3QsUJ4m7LOQcA37fzL4D3EVxSvJy09XvAr8MbAb+KOTae4HfMbNrgYy7/7T+ryMSD9UoRF71U3cfCR4ws0eB0HyBmf0ccArFeb4AjgK+T3H+r6Ax4A3Ame7+ipktAOtLnw1SbLpaVzr2r8EL3f2LZjYH/AfgDjO7yt2/0egXFGmEahQi0b4BvMbMtpQPmNk7zOxcirWJSXfPlrbjgePNLFNxjw3AP5eCxHuA4Od/BvxXiuuY/GHlLzezNwNPuftu4DbgHXF+OZF6KFCIRCitL3IxcF6pe+yjwB8AP6SYn/hqxSVfLR0PmgU2m9nDwK8D/wBgZr9OMVH+ReDTwFlm9ksV134QeKTUJHYa8JeVZTSzPyrN4Ntfms11stHvKxJGs8eKiEgk1ShERCSSAoWIiERSoBARkUgKFCIiEkmBQkREIilQiIhIJAUKERGJ9P8BsoR7ws524EQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(pca.fit_transform(X)[:, 0][y == 0], pca.fit_transform(X)[:, 1][y == 0], \n",
    "            color='black', marker='x', label='y = 0')\n",
    "plt.scatter(pca.fit_transform(X)[:, 0][y == 1], pca.fit_transform(X)[:, 1][y == 1], \n",
    "            color='black', label='y = 1')\n",
    "plt.legend(loc='lower right')\n",
    "plt.xlabel('PCA axis 1')\n",
    "plt.ylabel('PCA axis 2')\n",
    "plt.savefig('PCA_iris.pdf', format='pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52cb4f9e",
   "metadata": {},
   "source": [
    "#### Get plotting limits (used later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b896b7ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_lim_lower = np.min(pca.fit_transform(X)[:, 0]) - .2\n",
    "x_lim_upper = np.max(pca.fit_transform(X)[:, 0]) + .2\n",
    "y_lim_lower = np.min(pca.fit_transform(X)[:, 1]) - .2\n",
    "y_lim_upper = np.max(pca.fit_transform(X)[:, 1]) + .2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "291b35fc",
   "metadata": {},
   "source": [
    "#### Prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "58fd7db1",
   "metadata": {},
   "outputs": [],
   "source": [
    "m0 = np.zeros(p)\n",
    "alpha = .01 # Isotropic Gaussian prior\n",
    "S0 = alpha**(-1) * np.eye(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7edeab0e",
   "metadata": {},
   "source": [
    "#### Newton-Raphson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "731a36aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(w):\n",
    "    f1 = -.5 * (w - m0).T @ inv(S0) @ (w - m0)\n",
    "    f2 = 0\n",
    "    for i in range(n):\n",
    "        eta = w.T @ X[i]\n",
    "        f2 += y[i]*norm().logcdf(eta) + (1 - y[i]) * norm().logsf(eta)\n",
    "    return f1 + f2\n",
    "\n",
    "\n",
    "def grad(w):\n",
    "    f1 = - inv(S0) @ (w - m0)\n",
    "    f2 = 0\n",
    "    for i in range(n):\n",
    "        eta = w.T @ X[i]\n",
    "        if y[i] == 0:\n",
    "            log_coeff = norm().logpdf(eta) - norm().logsf(eta)\n",
    "            f2 -= np.exp(log_coeff) * X[i]\n",
    "        if y[i] == 1:\n",
    "            log_coeff = norm().logpdf(eta) - norm().logcdf(eta)\n",
    "            f2 += np.exp(log_coeff) * X[i]\n",
    "    return f1 + f2\n",
    "\n",
    "\n",
    "def fisher(w):\n",
    "    f1 = inv(S0)\n",
    "    f2 = 0\n",
    "    for i in range(n):\n",
    "        eta = w.T @ X[i]\n",
    "        log_coeff = 2 * norm().logpdf(eta) - norm().logcdf(eta) - norm().logsf(eta)\n",
    "        f2 += np.exp(log_coeff) * np.outer(X[i], X[i])\n",
    "    return f1 + f2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1906b98",
   "metadata": {},
   "source": [
    "#### Compute $\\hat{w}_{\\mathrm{MAP}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a659f86d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w = [-0.33747436 -0.61395942 -1.58386446  2.76171258  1.53012308]\n"
     ]
    }
   ],
   "source": [
    "w = np.zeros(p)\n",
    "old_w = np.inf*np.ones(p)\n",
    "convergence_tol = 1e-6\n",
    "\n",
    "while np.sum(np.abs(w - old_w)) > convergence_tol:\n",
    "    old_w = w\n",
    "    w = w + inv(fisher(w)) @ grad(w)\n",
    "print(f'w = {w}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f961209a",
   "metadata": {},
   "source": [
    "#### Fisher information (approximate posterior covariance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "75dc6309",
   "metadata": {},
   "outputs": [],
   "source": [
    "Sn = inv(fisher(w))\n",
    "MSMT = M @ Sn @ M.T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98e3c9e2",
   "metadata": {},
   "source": [
    "#### Predictive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7e2e72cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAABxYElEQVR4nO1dd3hUVfp+zyR0pEiwAQFBdy2rYmUVVFAp+tvVRcoqoSi66FoIYddFDCUhBNRdCcGOICIEEEHWTlFXBVwbihXdRZSIawGUIkVKzu+PmTtz5uTUW2aG5L7Pc597z3e+UzKZOe/9vu8UQilFiBAhQoQIIUMk3R0IESJEiBCZjZAoQoQIESKEEiFRhAgRIkQIJUKiCBEiRIgQSoREESJEiBAhlMhOdwf8xmGH5dCWLduluxtpAyHp7oEYQfXLbb025Ux0ZTqqsqZlRHomMjbtNo9Sik8++RjZ2dk48YQTqjfKzpqUPQeRVsltdP3M96rvtZzHetds2LCFUtpSlFfjiKJly3YoKXnXcz1+D2ypGNBsy+j03Ob7OQDapG2e+bsqjxB1HqujS/N5/AUAkUj1u0iXzxeleX3VXaR3773lKCgYgeeffRaXXHRRdLCpqore2Wf+Lsvn5bI0Kwf0d1FZXq5KO8/8XZXH3k3zTJ/9SFvKSb9+G8WKoetJikxZXpLOfgTdtsnLYKb8H1SwJXMTK0NEZKLyIuKUyUR5Kr2ff96J0tKJuPjii3HJxReLB2N+4HVgMnCb5DsyPk8m5wdeVT9N6zK5i8iBbcfmWZTn9ofh1rriEBKFAn4OUl7qsi3r9TuUyVD1+VD4e7xYqiorSUUOonpMrilT7sGWLVsweeJEOSGwsCEBvozs2YZsZP0wkbN3/m/S6Yh0RfpuSELVDps2kenkCtQ415PfcD4/P1xRlLqvx7asqb5OT5Uvy7Ppq4mul89NV95r3Q68uAh1rjqR64uXs2mdRWFybd78A6ZMuQd9+/TB2WedlTzA8q4nQE0SIj2R68mkHhsy0pGQ86yS2dxlz24JwjYtk7mRcwgtihTD8P/ie9mgYPP9M+2/H39nJn1WJrEYVqZyLTkyG4tC53ISYfLkSdizZw8mlpTold0M4nxZ2bNp3aIysjZE+qzM5s5/Drr2ZJD1V5TWte9VLkBIFIYQ/a+81JUK+DUwZ9KgC5j/9tIFm4FZFyQX6arSNjEI2bVx41d46KEHce011+DXv/qV97d8wJ1V4Dyr7qK6ebmpni1JyAiBrVf27OVLnGKSAELXU9pAqTt3hW05t+0EhUzrz6EAr8Fr26u4eDwikQjGjx0bDEHwcpGeI9PdbQiCTYvkpm3ykJGXzbObtEruE0E4CC0KS4i+W17qSkU5E32djirf63fSr8/TBqlu062VoYpNyGQ27iZe/5NPPsacOXNw6y23oHWrVu6JoqpKPeXVph7bMiYkpiMt3V1WRqTHy0XPqjSvr9Mz1dXVzyC0KNIMSlNjWRxqsPn70vVZmA7EMpmJu4kv57e7iceYMYVo0qQJbv/b3/wdkE0HatUgrNLh8/k2RG3K7rrBX5fHw2TQV+WZlreV6/IY1EiLQkeemdaG23psypno6nTcfN9Mvvde2/azjJ8wDWKL8lPtbiIE+Pe/38Czzz6Dv912Gw4//PDkBmwJQlZGVZdf7ajadvJ4HVbOw4QkRJ+VTi7rl6oem/Zk+royAtRoi4LS4N80/WojU/qq00lFP720ler+6wLMunK64LXIDcXnm7qa1JYFxR133I6jjjoK+bfeGv2gvLh9dCutbeuzIRZH7kAkl8lUd53MzbNJ2kamkuvyFKjRRAFkzgDsdz2UUhBC4mWc9KECrySQSsIyhV/9MXUl6fJNSMTJe/HFF7Fy5Uo8cN99aNSwofzt1mYQF5VT1eXkyXTYfFndrFymz7ejuvOfgahOVb5KX9aOrF2TcgGQBFBDXU88PHw+GdWGg0WLijBnTgGcY2yrqigef7wAixYVacu6/V6awM13tzbA1soQPassDtN8GXlUVVWhsHA0OnTogOuvu0498APuCMKNdeDI2DyZHq+j0mPzdXcVgYj0ZG3wz7K0CZEEQRKa/FpBFMChQxb6/yfF7t3bsHRpeZws5swpwNKl5di1axv8OgPd7XfOS/M1nVBkbioVQfD6NoFrU70FC+bjww8/RElxMerUqWNuLYhgQwYifUcmu9sQBJs2bUv2N/E6sgFcpmNCErJ2ZfXr5Lo8k/wYarzriYXzeQTpsqDUe/2qOgghGDSoDACwdGk5li4tBwD06pWPQYPKjNxPfnwOtn+nH5+LaZ0iuap9J8/PPqpiF6JnG2vCj8C1c+3fvw/jx49Fx44d8cd+/dTxCFHMQrd7rNvLJJbB58vSrFz07EBHVF6fTdJeZSZ5Jvkcao1FwcLyM0pL/ao6WLJw4JCEn39bOj4nk9+RaV1BlPEKXfBa9OykTQLWsjKi8oQAM2c+gi+//BKTS0sRcfYpVw3MqnzZgO1Gz0s/ZPmOTAQVOfA6bp5N00GThKwNDWolUQCuPy+r+oOqw3E3sWBjFulGhnQjbVAN+LZl3cYeTKyQn3/+GRMnluDCCy9Ezx49zAd6FiqZavCXPZu060bO31UkIuufqD4/SYKHDXGo9HXlDFBricLBoUYWbEyiV698VFRUoVev/KSYhWm7JmTp5cXFL5nq9+cGXuvQWQMqmShf9KbP6wVxTZs2Fd9//z3uLC0FAaIfjG5FNX+Z5LE6umevl/N3qEjAVObAhDhU5GGTlsncyE3zHR0FalWMQgZKg4tb+FE3WwchBA0bNkuKSThuqIYNm8VjFEH+TbK+pRuivgTVPxurwTReweu47bfMumDTW7duwT/+8Xf84Yor8NtOndQDrUmezYCtelbVxes5MB3AVXfTPNNnN2kbmUmeH/kxhEQRQ6aTBYu+fYvArptwyCKodRR+9t9v4gwCqSAXWb7OmlDl6XRY2d133xl1PZWU2BGCKs+UJHREwN55menbuWrANSEJUTk3zyZ12sp0eSaDvyFBOAiJgoHz2QUxSHgdfPi+8aQgIgnTv8dET9V/WZ7p3+zls7EpGzTBqMAHsVWxCJU1YUsa/BWJAJs2fY37778PgwcNwsknnmh2frXKJSW6TM7IBuSWgnMXzX5SpUXPJnfTPNWzm7SNTCXX5dnoCBAShQBBDSjpHKgONRyqn5VNn00tDJNLVk6ULi4uAqUURePGRQWmg7/N5cDtgM7LdHo2A78NWZg++5F2K9fl2ehIUOuD2TJ4+EwDrddNeb/+FjffU68vUKb5XvVlcDvwiywHkzJuyFEXj2DTn322DrNnP4abb7oJbXNz7Qb+oAhFdDchGBUppZIkZP2QpfnyKj2VXFWXTXkDhEShQG0ji6D+3qDhN8kEDVVw20m7sSZMyGLcuDFo1KgRRo8aJR98+bTbwd9rnWwZ/lnWhiO3uetkps8maRuZSu7k+fXl1+iFRKGBD2QcSL3pIjEvLza2upk2wJtCFo/gB22+jK31YnNFIsDbb7+FJUuewl9HjkTLnJzEl1A0zZUdgN3EJ0TTZ22Jwu20WUBNTqJ8lUxFVOyzCWmYyFRyXR6bb/JjNhyIwhiFISjNPJ+5bZ/S9Tf41S5bj6zOVP+NXtryEpQW6YvqdZ4pjW4j3rJlS4wsKNAPsn5cfF22dTswSfNy3V0nM302SdvI3MhN893qIiQKK/g9CEVXUhNm8At+q3CTv0Gno8r3MoD7//lmDrmL4heytEldKsIQuZxWrFiBV199FdPKy9G4USN/F7nZkgSfJxrgRTqyOvhybHmbvFSRhBsi8IskLAnCQeh6sgT/MuMWixcXYe7cxEpqSinmzi3A4sVFaemPqF4v+X60EXR5r/BqTfBpG2vCxuVEaRXGjBmNdu3aYdj116tXSOvcUKxctwrbdiW2KXkBdnoqfTZP9czWJcsTEYJORyXX5Znk2+pJEBJFGuBsFb5sWXmcLJxtOXbv3gY3ezbZFAl6gPVSv+nLV6ZBZiXoLAm3cQldmr0WLXoS7733HiYUFaFe3bpmg6zNQM2n3ZZ3INLXyW3uLFREwLcrypOlTb/IOhJQIQUE4SB0PbmE89m7easkhGDgwDIAwLJl5Vi2rBwA0LNnPgYODG6FNQtKg3VBmerb1mFSt1918jAZ2EX5KrJg67WNS8jKs7L9+/dj3LgxOOWUUzDgqqu8k4Tqzdn2bZ29i+rh5WxaV58JWdg+6/JsZG7kpvmmOhb6oUXhEW7JmiULB9G0+xHOhxcH3+r00hfb32a6ISMPGUG4e7mwIwjn+bHHHsX69esxaeJEZEX9UNFMvwjDhCREcv6uK8vK2LKy+kR3FqI6ZfXK0rI62bSsbTckoSpnWofLOtNKFISQRwkhPxBCPpbkE0LINELIekLIh4SQM1LdRxO4GbycmAQLPmYRdF9Mv3d+tedG/1CCyN2k0g0qLuHEJvbs2Y2SkmJ07twZ/3fZZYl/uOoAIudZtzOs7JLVIWvbJC4BmOXr7rIybD7/7CbN163S4/smgg1BmP7ALH/86XY9PQbgPgCPS/IvBXB87OoE4MHYXYp0DUROuyaDhEMSy5aVo2fPfDRs2BRr1jwTd0Hl5U1BRcVINGrUDH36FLnqi5+uF7/r81q/l/4E/bfwULmMbOowtSRY2X33TcO3336LhfPnJ7YRT/fl9IO/iwZL1aAuy9fpq2QyHdO0Sm6ja5NvquNGl0FaLQpK6esAflSoXAHgcRrFmwCaEUKO1tfrVw+DgbNVeM+e+cjLm4Ldu7ejsnItcnM7okGDpqioGIlly7ydge33d0el4/Y34AdUv2tVX1L5HbENYutiEypXk3Pftu0n3H33Xfjd736HLp07qwdvv6bKAuq0SXm2DJ+WEQerJ5OZfFG8foH4fnttz4REUkASQPotCh1aAfiaSW+Kyb5llQghwwAMA4AWLXIBJD6TVL492rTbp08RnHUTbGC7snItgERgG4geb+rm77Ap57YNXXmR3ESmS6cbtkFtmyC2Ll/mamLv//jHXdi+fTtKJ0xIJgOV68fEdaS7TM66dks0TtrN3fRZl6fT9yLX5dnouNFVoEYEsyml0ymlZ1FKzzrssJZcXrr6pNdhz5MQBbb9mP3k59+fqXV5acfvfqiC1iauJ1tLQnT973/fYNq0cuQNGIBTTznFvYXgx2BvWoYvx6Z52JADX0b3rCMQE5JQ9dsLSejKu9U1QKYTxTcA2jDp1jGZFXz+zKzaNdOTB7Zt6/ICv19mvMLmN5FKiMjAhtN1QWteT9cmf584cQIOHjyICUVF7gb7oFxRpvmsXKTjyHR3EXHISEhHADryUclUclFfbMt70TVEphPFMwAGx2Y//RbAdkrpt7pCMqSDMPQvEcmB7TlzqtCzZ37SYjy2riD6wOu6hZvfiI2e6nfttk0/IRvUTYPYblxM/H39+v/g0Udn4sYbbsCxbdsmD/yqu+k52DaznXSEo8sH3JMTW1b07EBEGKq0jUwl1+Xxf7sONrouyqU1RkEImQ+gK4AcQsgmAOMB1AEASulDAF4AcBmA9QB2A7jWj3YpTU/sQuynTgS2HXeT44Ziz8DW1ZNKqPogy8uEfvsNlTVhE5OQuZpEz7K7c40fPxb169dH4ejR9oOsXxdgrsfqitKywdcmz+TZTVol1xEAg4pVq1D4xBOo3LIFuTk5KP3jH5HXubO8vE1bPpVJK1FQSq/W5FMANwfTdvSeysFLNliyge1on4gyRuF20DUtZ6Lnx8AvqiNIQnHqThVp2bYhIgxWrrvee28NnnxyIcYWFuLII44ADh6sPjjbrF9wSwxOWiXn81QDvKxe0V1VBwsbwrCRqeSCvIpVqzDskUewe98+AMDGLVsw7JFHAEqR16WLvB5dOz6Xy3TXU+Bw+1n73Z7JGdgm9bht36/6bcvqdC1+c1ayoCCKM5iUMQlei/JYncLC0WjRogX+OnKknCC8ym1JQzbY8/mqcmy+6M5C1B6fHxRJyNpUlCl84ok4STjYvW8fChcuVNfj5kvtthwyf3psSuB8dqmyLij1py2/6gmifVGeaX8z8fOxDVbzaZ3LySYuIcr7179exksvrcCUe+5Bk8MOqx5/YAd/0ztPGibkoZoeC4jz2bTo2ebu5lkm8+NNRDMwV27ZIpZv3Wpdl5v2TVHrLQoWqXzz9KutIL87QXweXvsrK5/G3xAAd5aEo+c2HpFog2LMmNFo06YN/nzDDWZWganF4MVFBeiJQJXnhSRUloOIEGT94fVMZCZ5sfzcnBxhVm6LFvr+aOr28wseEgUHnz9fbVvpqscvskgluZrC7z6pBn1RTEFXThbA5p9N3U7//OdTeOedd1A8fjzq16tnN/D7QSYicjCVB0ESLFSEYZIW1aHSVekL8kv790fDunWTshvWrYvS/v3dD0YB/ChD15MEzmcdtGvHr3Yota/DtIybuoOs321//AA/eJvoi60Ae/eTyPVUVXUAY8cW4qSTTsLggQPNpqvqVme7mdbquJRY1xJ/l5EDny8rL7q7eTZJy2Ru5Ip8J2BduHAhKrduRW6LFijt399u1pNJ2x4REoUGqRqU/GgnXWRhm2ciSzUZeGlPNPjr9FXBa1NyIQSYM2c2Pv/8cyxZvBhZhOgHfFW8wetMKOeDZO/OM6/HykX6qjsLWZ5bwlDJdVaCmzxEySKvS5f0Wg+aekLXkwECJmtf2zGtg99s0O3mg7bt+l2/agwIGjqXlOmgz9ZlYl041969ezBhQhE6deqEK37/e3srwA9y0FkBOoKQ6enqZGUsRHWLnlWwIQm+PdNybvVM2/W5ntCiMESq3nD9aEdXx+LFRdi9e1t8rUZVVXR1eMOGzdC3b5G3xn3qo6lOKqFyO9n0UxafENWvIooHH7wfmzZtwpzHHjPbRtx0Q0BbcnCr76QdmJCETFf1bJKWydzITfNNdfwo47GO0KKwgF8kbtJOUHU453UvXZrYImTuXLPzut2+PPn1UuW1/qChcjuZuJFUMQk+vWPHNtx55yT06tkTXS+4IDHY67bZkG3fYVJeFP+QkZCbC3BvUaieVWmZzI2czdP9UGwHEzdlfKwjtChcIBVvun60IaqD3SJk6dJyLF1aDgDo1SvY87r9/syC/h/YBKtFMpm7ic1ny5vGJJxrypR/4KeffsKkiRPNB2g3p8u5HeRFg5Iozctld9mzylJQ5dnKVHLTfFMdL/oB1RFaFC7hldxN2wiiDtW25l7btH1J8lLebRt+w4ZUVC4nVkdEHM7zd999i/LyMlz1xz/i9I4d5YTAWwwmJOLnBahlzrMJdGVlJKRKm8hk5U3Lmdbhh35QdSAkCs8ImjD8qL/670G/rXnQfZL1zTbftJzud24LlZUgkuuC1KpV1+xzJAJMnjwR+/btQ0lRkdrFZONW0um4IZWgLBb2WUQQNmmZjJWryMPkC2zzgzHQr1i9Gu3y8xEZOBDt8vNRsXq1uLyPA1NIFD7By//E7xlI4jYSdTsxiV698jF3bhV69cpPiln4TUw6edB6QcHE7STS11kMqmvDhi/wyCPT8afrr8dxHTqIB2c3A7vNoA6Y67G6zjMrN7nzMpMvmIwQVGVM5Lo+mNbhss6K1asxbOZMbNy6FRTAxq1bMWzmTFSsWhXoDyIkCh/h5v/01FNFqKgoQFVsoRKlFBUVBVi8eHy1uv34HjjbmrMxiYEDy9CrV37StubpJgtZGZMXOC/5gNqV5CYuIopJ8Hmm5FFcPA516tTB2MJC/cBvuo7Cb7eTA1laRByqu4w4ZF8ME0IwJRybPJO63dYZQ+HCheJNBJ980rw9FwiD2T7D+Z+bDCjODKRly8qxbt1rKClZg3nzRmLZsnLk5nYEAPTpU1ytfrdBXKdvptua69qy+Vt1/Qp6coBbyOII7LONy4nPU7mi+OvDD9di/vx5GD1qFI4+8sjoNuK2M51EchWJ8PmAeuM/FRHYEITsWWcl2Kb9yLPRcaPLQLhZoELuF0KiCAgmgx8hBAMGTMG6da+hsnIthgzJAgDk5nZEZeVanHjihUkDuk3dOphua+6VmIIgAFG9KlkqiEg3a0mkYxLEZq8xY+5A8+bN8be//lU8sNvs22RiUfCEIHs2sShUb/mmeSqrwYQQgrQgTOGSIJyyuS1aYKOAFJI2EQwAoespQJh8JyKRCEpK1iTJKivXomfPfOTlqQ8vMu9HsnJVlYcvqw8w9Qz4WX86oHNh2VwrV76GpUtfxOjbb0ezJk3UVoDpXk9uL8CcIGS6MouDl7F1yJAukuD/Th18+FKX9usn3kSwXz/3dTP1yxASRcDQf9co5s0bWU0+YMAUKUmY1g1EV2GzM5qcYPaiRUX6wobteH0h86MPqajDxOUkyvNqSQAUhYWjccwxx+CWm27SWwE8QaguP+ITzoerszD4f4SKEFSWhK1lIeuDijwEeRWrVqHd8OGIDBiAdsOHJ882UsGWUBRl8zp3xvShQ9G2RQsQAG1btMD0oUPtNxG07FvoekoBnP9DdddINHDtxCQqK9fG88aOPRMlJWsQiai5nFL5gMXGQIDoWom5c6Pt9eyZj6oqikjEzCejasck342+rY5tH3SQBaB5HbckIJoSy8qd2MTzzz+DN9/8N6Y/+CAa1KsHHDhgt4ra9q4iFUAeo+BlTtrNXfXsJq2Ta/IqVq3CsBkzko8snTEDAMSDtA/Wgwx5nTt7IwYXCC2KFIL/HxFC0KBB0zhJ9OyZj9mzD8bT8+aNhMlUWdlLgROk7tkzH8uWlWPQoEicJPxaYGcC0za8WC62cFuXKRHJSEYW5JZdVVUHMWbMHfjVr36Fa4cMMZ/JZKrjZlqscxe93cusA5M6WdhYFCYk4cKCYCGdbSQ6stTgyyVcC+HF8tDBY9010qLw+83ST/B9c2Y1nXjihfGYhDP7iZ2u6hYOWThWBQDXW3X4bVUcqjCxMFg9L9ZGRcVcfPrpp3jyiSeQnZUlXiDHE4EuRuGGJGSDPCC3KpxnFWwtCR10pGOTx+RLjyxl5YZ9ddZCxK2T2FoIUOreUpDBJ+KpkUQBJD6fTBy4+L716VMMdnZTJBJRBrJN643KxKuwWbKw+azckoVIbirTtZMOgpLHE/RkwZ+BzbuaHL19+/aiuHgczjrrLPTp3Tt5OixLAqars3VEIiMQG5eT6Nnkrnp2k5bJTPIE+bk5OdgoIIvcFi2sB2PVWghfiCIAq6TGu56CtOa8gu2X6XRVu/ppUkxizpyquBvKZsuOIGH5ew0Muo/bjctJRyYyPYcspk9/CJWVlbiztBSEH8BFz24IwNRy4ElAZFGwenwZ0Z2FijBs0qI6eLmLL53yyFJLBLYWIsDBrsZaFDzS8dZpgiD6lbAQoquw2ZiEsxmgyK3l5m3etl9+1h90P3iYupxU5UwD2z//vAOTJ5fikksuwcXduunJwGa2k98zndgP1k9LIs1WBAvnTb9w4UJUbtni7sjSWP2+r4VIwdtUrSEKgB1A09sPHkGSmOkqbNu+qPT8/nvY+lJN+KJgNJ/WEYFMrrrKy6dgy5YtmFRSYuY+kpGB6TRZr0RhY0mkmiQ8EASLvPPOQ9555xnpquov7dcPwx59NMn9ZL0WIgXkwKJWEYUD9jPOFNIwITF+lbZo1ba4Tju3lh+DsawOkTwoqyqo/y1PAqzMeeb1ZFaE6HCiLVt+QFnZPejbpw/OPvPMxHTYgweTB382Tak8X2V9uLU0nA9ZRQYm5KF6NknLZCq5ab6tnmG5uHXy5JOo3Lo1ap3066e3TlJMDixqJVGwSPUbqg6y/jz1VPT4UifI7azBaNiwGa68sihl/bDVSTX86hNPAra6NtYD73a6885S7NmzBxOLi8WWgm6NhJvtOmwIwvmg2We/rQhTQnBDEmkiCBZWayHSSBAOanww2wQZ8H9IQvWXqMTCuYqKArAL9XTHl4rqC7Kvbts1eXH0o14dbElGRRKs3NS6qKz8Cg899CCGXnstfn388eLZTPygr1pH4SdJ6MiD/9B1A7uOGGTlnLSs/oBIQnkOhEnbtnBRX8Ubb6BdQQEigwejXUEBKt54w5eu1HqLwoHz/8iUN2S2P4QQ5OWVAQCWLSuPr4nQ7QfF12f7t/nxdi6qI4i+eOmrlyA0m6fSE5GC6CopGY+srCyMKyyUT33lXUyqabF+zn5iCcDUgrCxKNykdXLTfI3OTbNm4aGXX4ajEV/7AFSLW1SsXm3vVrLppwQVb7yRFPvYuHUrhj36qLCPtkirRUEI6UUI+ZwQsp4Qcrsg/xpCyGZCyNrYdX3QffL7pcArEoSRIAsHtmst3PxdHn9fKa831RARBJtnc3388UeYO3cObr35ZrRu1Sp5IDcd3HUL7NxaFQ50JKHK1z2r0ry+iVxWj2UdFatXJ5GEA9HK7IrVqzHs0UeTDxZ69FH9vlA+DDyFTz4Z2FkVaSMKQkgWgPsBXArgJABXE0JOEqg+QSntGLtmpKp/6SIM0Wl3VVUUjruJheOGsqvfcxdTWm862jXlXpUl4uYaN64QTZo0we233aZeLGdqObi1GmQWhGiQF5GBTO7WqpDJdHJTgtDoFS5cWI0kHPBrH6wHax8HmiDPqkin6+kcAOsppRsAgBCyAMAVAD5NY5+qIfFGH3xbooB1cfG5AIAOHTphxYpp6NFjOADgiy/eirug3FgWNn+Prb7b8jo9v+qRQVRGNbDr8m1cTm++uRrPPfcsJk2ciMObNVOvujZZkS27Ox+QLXk45djyvEx0Vz3r8mQyN3K3elAPtPzaB6PBOqC3KlfrMwz7kk7XUysAXzPpTTEZjz6EkA8JIYsIIW1EFRFChhFC3iWEvLtz5+Yg+hr4G7MTsF6+PBGwnjt3BDZseAsbNryFL754C927R0li+fJp6NChE3r0GO56Pyjbv8ej9W48Bnj5nG3HIK9gycK5u7soCgtvx1FHHYXht9xiRwiyNRW6hXk2BOF8gCrikH3IKsLQEYisPi8k4eINXjbQEqDa2geZbnyrD4u2bQPTVmdVWPYl04PZzwKYTyn9hRByA4DZAC7ilSil0wFMB4Bjjz0rsCHd65u1CmwMYvnycixfHrUWuncfDkKi5LBhw1sAEkFsp1xNhuwzD/J/AejrVpEEn9ZZGEuXvoDVq1fhwfvuQ6MGDaLrJngikFkNOjLwMybhpE3ups+itAy2byJu2hCUES2QIwBuvPjiakFqXxbTwV1g2pFXC6Sz+i7fmEi69vshhJwLoIhS2jOWHg0AlNLJEv0sAD9SSpuq6j322LNocfG7fndX0J9g6qWUYsiQhKE3e3bUXcDLTM+R0MH27zDRN5lFpJKr0rpn2QAue1YN4qJ0VpbYjeTIs7Ki6ezsaDo7O6HjPGdnR/WyswGgCp06dcTu3bux7sMPUScrC9i/PzrAOwvtZHdnk0BnFpTsbrPADqjuojIlBxWh8HLTtFu5rY5G32Ymk1DXctZRu4ICoRupbYsW+KqszKou07+fDBmyhlJ6ligvnRbFOwCOJ4QcC+AbAFcBGMAqEEKOppR+G0teDmBdarsoh/PZ+0kYooD13LkjqrVRUVGAAQOmIBKJeG7f9s3cy5u8admgrQUbiOISjlyUL1pprboWLJiPjz76CPMrKlAnEkneIZYlAvZuGpvgLQ0bKwNQPztwY1GYpN3KTfMty5gukPODJACPgekAXv7TRhSU0gOEkFsALAOQBeBRSuknhJAJAN6llD4DYDgh5HIABwD8COCadPVXBr8GNYckli8vR/v2ndChQycAFCtWTAMAHHZYDrp1+zP27t2B5cvL8dlnr+H00y9H377FPrTtN+Fldn2AvdXD5ossFLasqbVy4MA+FBePxemnn47+ffrYxSFkMltCMHE56SwFL24nG5lKbqsj0Pe0/gGJqbF+rGMIMjDtBmldR0EpfYFS+itKaQdKaWlMNi5GEqCUjqaUnkwpPY1S2o1S+lk6+yuDH/8fZ6fXHj2Gx2c4AQTt25+DunUbYufOLVi79llcffU98RPw9uzZjqoqf74cNn+Dl7/X62elG3P8gAmJ8M82MQnnecaM6fjyyy8xeeJERAhRk4OKDPxYL8F+qKK0G+Lgn0VIJ0kwf5Pr9Q9MPX6uYwgyMO0GmR7MPmTghyvqyiuL4MSMCCHxgDaAODlce230X9ajRz4zjdZd2/ymguk6QzvVriZVbIOVi9xOtmQgurKygN27f8akSSXo2rUrelx0UdStxLudRK4nSt3HIFTkASTHJpy7zLUUhKsplQTBQTXIm56J7ec6BmVgOmBSECEkCp/hddBzBu68vLIkopgwYQ2uuSYrnnZzAh4L1SaDffoUua7XL6SaPAB3gX3TYDh/3XvvVPzwww94evFiEKC6NaFyManuQbicVOTAltM928pUctN8Q12rQV5Sj9/nTFTb1jwFloMM4aaAAcC7e6V6UHvcuDOT0qJV2abfI6+bDLLtuYHNOJEpUFkYsnyZlfHjj1twzz134w9XXIHfduokHvRl7iST4LWOEGTWgMqlJLrrrAoRRJaF7Auhqst20NToKtc/GLZp5S6yQRoJwkFIFAHB7f+WDWr36JGPxx47GHc75eZ2xGOPHUSPHsOTFubZwlmz4RyLOnhwJH5canR9hvmrte63HATc1CuyFkRuJ9lgz+vrrAeZyykSAe6+ezJ27dqF0uLihIvJ5GKJQjYDSkQOOnJREYqKaNhnlWViktbJTfNFugb6ykHesI68887D9KFD0bZFCxBEp7JOHzrU/WFHGUAQDkLXU8Cg1M6lkQhqRwftJUuKkZ1dF7m5HXH66ZfHYxLt23eSrsp2vluqdh2ycLYBAZLdWbb9lkFWj0n9rI5f/XELFamwzzoS2bSpEg8+eD+GDB6Mk044QTzgqywIXaDayQfUA7+TZu+sjihfZjmoLAof3/p9rYuD9DAhy0He9Sl4DjKEGHiERJEC2A5ybFB79+5t2LDhbXTvfiuuvLIIFRUFWLFiGrp3H47evcdz7ZifgCfbZNBr7KMmQmVR8GkRObBWxcSJRQCAorFj5bEHERHIAtOmbifAnBhY6GS2A5tM3y+SsOiPaDrsV1OnpmewzlCCcKAkCkJITwCtAbxMKf2KkQ+llD4acN9qFNxYFgCStvVYseJeAMkznhzIgtONGjWLE0/CWojuI7V8+bS4u8mJUThtBrnJYCaCHfhFMptZT7LYxGeffYo5c2ZjxPDhyG3dOrEC283llSRUrqCgLAk3JBGQFeHnmgdPyHCCcCCNURBCJgEoBHAKgJcJIbcy2bcE3bGaCDcuR5NzKEQbCjpxjl27tmHx4vHV4hlffPEWOnToFK/LiVmw7qxUf4fT8ZuRkYMNMbCXI3PiEexzUdEYNGrUCKNvu00cf2DvstiEijhkVohtXEJ052UqwhCRhow4VORhY0VYfnmCPLtBC/4zPwSgCmb/HsBFlNIRAM4EcCkhxBmxDvH3x/TC7iVJfw6FM9D36JGP5cvLMWRIJB4MHzBgCvbs2V5thpNDFHwd/Pnbflj8mfR7kMUaVLps2iYmwZLHu+++hX/+cwlu+8tfkHP44cluJhNrwe0iOj8sCZV7yq3LyI8vjIfBNsizGxxU2/119erM+jFYQEUU2ZTSAwBAKd2GKHE0IYQ8CaCuolwIn8DPgJo9uypOBjKyYJGXV4ZIJIIBA8QznAYOnJpkmaQ6NpFJvxkRKTh3GwuDJ4joRTFmzO044ogjUDB8uDpoLcpTWQkyPRsyAdSWBK8jSsssCVEdfpGEW1BqNh3WA5zdX5NWes+a5dsZ1qmGiii+IIRc6CQopQcppdcB+BzAiYH3rIbD5GWIEIKNG6PTYgcMmAJCCAYMmILc3I7YuHFtNfeTzPKIllO7r3R99YpMsExE8QY+3yQuoQpW89Nhs7KAV15ZjtdeexVj77gDjRs2tJ8Sq3M/2RKG82HaWB8qUrBxNclgah14cdkwZQNb8xCD1LW1aJEv9acaKqLoB+BtXkgpHQNAeIBQCHuofzsUbdtG11DMmzcSlFLMmzcSlZVr0bZtx7hFobM8qqqqfDlG1cvfkg6IAtQ2Zdk7X59J8DqqW4UxY0ajXbt2GDZ0qFmAWrbuwTTWYOqOYp9VFoWOFFRpnVyXx+p4tCJY+LrmQdBWKlxbqYR01hOldI8i75tgulN7kDwLCQCqT2WVHWbEz3ri116w5Ro2bIp580bGSWTgQHcznCg99Gc2iSBzH8nyZHKeMJz7U089iffffx9zZs9G3exsu6C029iEzsWkIwc/SUI3uJuShFsoynpe86Boy+/tPNKNtB1cFBRSdXCRFyxZEp3KOmBAYirrvHkj0LBhc/TuPT4uc+5A9YOLxAvtxOso+KmzQGJfJz54rYIJUah0bILGujyTu2pwj8Q+Tt4C4A8ocmYs8Xf28KFIBKhTJ5pm71VV+3H66SehQYMGeP/tt5HluJD27Yve9++PupUOHEjcq6rEBxOpAt7sHah+l7mdnGf2rpPxzyqZG7lbPb/K+dRWxRtvYNisWdVOupt+7bWpnX5rgUw9uKhWgp3KCgADBpRhwoRzsWHDW2jf/hzs3v0Trr56CubP/wsaNGiK3bu3Yf36N5PqkC2Mk6X5dRSOxRFE8DpTLQ+T2U42gWp+KqyIUGbNmon169fjuX/+E1mAehqsiBDY3WFVBKFzMankqrsDG6vCrdytnl/lfG4rvvvrokWJBX19+2YsSeigtSgIIR0AbIqdW90VwKkAHo/NhMo4HAoWRdSCKEjaHRZIbCV+2GE52LlzC9q0OQ1ff/0BAKB9+3MwfvybSbEIt4O9V7+9W71UWRQyKwKQH3Eq2xpcZEmwFkWdOmKL4sCB3TjhhA44rkMHvP7SSyCsJSGyKFhLQnbUqWoKrciisLEkdBaFKG0jM8mz0XFRruKNN3w5fS6lZJRCqCwKk00BFwM4SAg5DsB0RAPZ83zsX62DaBZS9+7DUVm5FgCwc+cWAMDXX3+AJk2OAgAcd9xvASC+XkK2z5MJDoWXNRFs/lyRro3rzMR1pZr1dP/90/Ddd99hckkJCD/IqwLYOiLQxSFEpMCnTWMSLESkYeuCSjNJCA8mspiuWrF6dXRdxJAhaDdy5CE71dUNTFxPVbFjS3sDuJdSei8h5P2gO1bTwG+hUVExwqjcjh3foUeP4cjLmxovH+7H5B0mbiYVEagsj+3bf8Tdd9+J3/3f/6HLb3+rXoHtxnKwuQD/LYlDyIpwoDyYSGdVUFot5uCsiwDUW35UvPFGjXA/mVgU+wkhVwMYAuC5mKxOcF2qeViypAjz5kWno0anqo7AihXT0L59J8yadRDduw+Pn40tghP0duAHSaTbOggCKitC5Jri9VSzmkzjFZEIcM89d2HHjh2YNGGC2RYcNlaDDUEAwVgSOh3TPBsdURmLcq6nq8baKFy0yHpdhEMuNWHRnQlRXAvgXACllNIvCSHHApgTbLdqDtjg9YQJ52LevIJ4cLp9+3Mwb14BvvgimnbcTIcdlpNUx7hxZ6LK8TuHqAZZTEMX++DdS7xcpCcjB+f69ttvcN990zBwwACcctJJdu4mNxYEIB/4TQdSGWFkKkm4gKuV2Ezf3BCNG3LJVGiJglL6KaV0OKV0fiz9JaX0ruC7VjPgxCO6dx+ODRvewooV0/Dll2/jkktuBSEkZkkQXHLJreja9U/Ize2InTu3IDe3Iy6/fCzatDktacGdavJB9RPv/P8RZoIlIprBZKIvK+PW5cRv/JeVBZSWFuPgwYMoHjMm4W5ypr+abPZna4GYup10JKN6FqVlMpM8tl+2cFnOaiW2oA03RFOTFt1JYxSEkIWU0v6EkI8AVPvPUEpPDbRnNQjR6ahTk9xLL72U2DLcCWwTQkAIwQknXIirr56Cp5+egBNOuBAnnHAhGjZsBgDS9Q+qM7Bt1kqkG7YzpkT5IveSzprwI2bx3/9+jlmzHsXNf/4zjs3NNd/51YYAAL07SRZr4BEESWRAPEKE+HRV1awnRf2lffsK10WU9u0rLVOTFt2pgtn5sfvvUtGRmgxnOqwIfPzhyiuL45bA7t3bsGLFNPToET2kiJ0aywfH2bUZzvkSIt1Mg9du2RKIKN+ULERbdLBEMWHCWNSvXx+Fo0aZkYLK5cQ+s7CJTZjcVc+itExmijSQhAPlSmxN/W7WRbghl0yFaguPb2OPjSiln7J5sfUUG4PrVs0Bu2aie/fhAJBkWVRUjEia0QQkgtXJ23dEy4jWT5hu9VETYON2Eul6iU3wzyxJvP/+u1i06EmMKyzEETk5ibUStkFsW8uCvTvP6SKJgPySopPofJs5ZNFnGdHIZjbVpEV3JgvuPkY0eH03gPqx+1mU0nOD7549MnHB3ZIlRdi16ycAiB9jCkQPD9qw4a24+0k8oFOj7TuAKCmZ6gLu3uZtBmgTuU3wWXdnB3p2iw7nzg/yIjeSaGGdbIuOrCygbt3o/Xe/644PPliLLz7+GE0aNUosptu3L+GCkm3ZwU+P5e9+LqyTyWQ6snybPJN8SRln/UO1rTC4DfxcLabzgdgOxa06ZPC64K4Toovs3gDwDoD/AejsX/dqPnr3LkJe3lQ0atQ89pY/FXl5UzFu3L+li+ecwLXprq82upkCvwwdFcHo2jFxO+lcTv/610t4+eWXUHj77WjSuLFZUNrWanCeVTEJkb5MxkJVh0pmkucWsTpNTqJztZjOpz7XpJlNKpgsuNsPYA+ABohaFF9SSsO5mpYghKB376Jq8QKRJeFYIIQAy5dHYxSURi0QNg7BbhqYiElEF+c5aVb3UIIutmCio3I52QSuVS4nQijGjh2N3Nxc3Hj99foZTDaxCZ4MbF1OPHRkIEOqLQkGJjOHrBfT+UhsNWZmk+YzMSGKdwA8DeBsADkAHiKE9KGU+nPCRy2DbiM/JzDtLMhzSMJxWR13XKe4BcLOdIpuMx7VXbKkmNlmPGGtyHaXNe+7ex0/eEoWc2CfdS4pE4tBdgCR6HrmmaewZs27mDV9OurXqZNwNYlWY5vuAsuSiY2Vweb5GZeQIUBLwoHJzCGrwdqwz6Yrqg/5mU2Gn4cJUVxHKXWc/t8CuIIQMshtv0Kowe4DtXx5OTZseAtAIjDt6IhmOs2dOyI2Syo/LnOIQDR9dt68zJk+a0swpu4kWRkRgcj2b5JZE1VVBzBuXCFOOukkDLr6avPgtcqC4Ad+wNzlFARJmFgnMvhAJKX9+gljFOz6B+PB2oIkdNt1OESycetWECBp/cAhM7PJ4v9jsuDuXQAghBxBCMklhOQCeM197xIghPQihHxOCFlPCLldkF+PEPJELP8tQkg7P9rNdKiOLuW3CndOshsyJBInCV6XJRUnbjFvXvTwot27t/kWx3AbxLapW1aXSZDb1L2kC3g7V0XFY/j8888xacIEZBGidjOJnk3iFDYuJ7+RapIQlDE5ic5oMZ1Ff3RxB3ZrDiBKEs7Xsm2LFodGINvy/2My6+n3AKYAOAbADwDaAlhHKT3ZZRederMA/AdAdwCbEHVxXc1OxSWE3ATgVErpjYSQqwD0ppT+UVWvzawnr64YPyHaNJCdRiub6mo604k9LtVBz55202d1am6JQvfmL7vrXEoiMnAOJNJtIy7bPtyZ5VRVtQcnn3w82rRujTdefRWEn+XEbyVeVaXfRpyf5WRCJoDYylC5o/hnUVomU8ltdbzoc9DOerJ5ex4ypPoKY0TJoGr2bLQbOVJowbRt0QJfTZli3/l0QPB5kGuu8TTraSKA3wL4D6X0WAAXA3hTXcQI5wBYTyndQCndB2ABgCs4nSsAzI49LwJwMfFpJGc36gMQf8tesqTIj+pd94UlifbtO1U7/5oldpuZTuxaCwd+koQNTNxG7LOMJFg9W0vCxGqQxSWys4GHHroP33zzDe4sKYmeNWESlxDFKUyJwIQETAhBJEs3SfiAvPPOw1dlZah6/HF8VVbmmiQA/XYdh3QA26UFarR7LKV0K4AIISRCKf0XACHrWKIVgK+Z9KaYTKhDKT0AYDuAav9FQsgwQsi7hJB3d+7crG2YdcU4A7SzKM5PV4wJ+L4A0dlNAHDccZ0AiM+gYC2EHj3ylYTC6rMwnT6bqiC2ysowqU8W7HbjelKts9ixYxvuumsyLu3VCxd26aKeDsu7nNy6mxy4HdRtvtMZ4HLyDS7qLu3bV+zKisUdXG0weIjDJJi9jRDSGMDrACoIIT8A2BVst+xAKZ2O6KFKOPbYs7TfDD5gzK5kli98CwbyvqjPoCCExGY6JdxHoplOQDKpOO6miopojEJUt/3fYC5XuZlUdaQ6NiGzLiIRoKzs7/jpp58S24jrSMKGGFjILApWZoI0vOEbQ9E3306ks4RuRXVN2prDFCYxikaIrqOIAMgD0BRARczKcN8wIecCKKKU9oylRwMApXQyo7MspvNvQkg2gO8AtKSKTtvGKK65JmFUPfaYeiWzG5jGQfi+6FZV29b/1FNF2LPHftNAP60Jk1iE7FkUj2CfTVxLTozCGfgJEa/AZmMTTkzCuW/d+i1OPPE4XPH732PerFnVYxKyY05FsQmTabImLinnWXRXPYvSMplKbppvUcZ0Rbbv/THEIXkgkW6sV8QotBYFpdSxHqqQiBf4gXcAHB873+IbAFcBGMDpPIPogUn/BtAXwCsqkrCBaKO+efMKfLUoliyJTkl16mSnpPbuXaTsS0VFQdJ0WFZXdYiRrO9RMkiUdSwQr1t8+OFy0tWtIhE3loTO5aS6Jk8uwb59+zBh7Fj1QjrZ4jpbt5NJ7EHmpuIRdGzCR3g6kS4FUG4wmInw+P8ziVEEgljM4RYAywCsA7CQUvoJIWQCIeTymNpMAC0IIesBjARQbQqty7bjMYkePfLx2GMJ/z4b4PbahkkchO8LG2soLj4Xc+eOSNKtqCjAU08VWfcnOiCakYpNnaZyW2tCdneeVVaHqC0VWZhOh/3yy/WYMeMRDPvTn3Dcscea7wqriz/IXEssUhGb8AKf2/EcMM5kd9shCJMYRWCglL4A4AVONo553gvA9xXgrH/fedt34gSifZfctmESB5HHGijWr48edOTI3G8dTgEkdE3K+umBc+NyktUjIgwb15OMBFiXlOwqKRmHunXrYszf/iY+iMgmXmG7jTifbwLbwTLV1oSi3kN+xbMPyCT3ljZGUa0AIW0AXEUp/XswXfKGTFtHoYuDOI+ivgCotvbBdutwx/1lE5vw2+XkhzUhsxr43WF5onBiEexdtjusat3Ep5+uxTnnnI47Ro1C6fjx5rEJNkah2x3WzboJk7sDnaXhhSg8xif4wPVlHTti9sqV7mMUh7hF4fuutAafh9d1FCCEtCSE3EQIWQngVQBH2vQxU+G3K4aHLA4iW+fAp72ufQCi7q9lyxJTZp3ZTl6mAWeqy4l9NolLiKwP2TV27Gg0b94ct40YoQ5C285yUrmdVP8fHSnI9DMQot1fZ69ciSHnn5+0InvI+eej8MknERk8GO0KCtS7wx7i8H1XWo9jm+oo1MMAXIlogPlXAJ4CcCyltLWnFmsJ+NjDgAFl8TRgNuDL1j6YkgVLNMuWlcenw8pWZHuxJEx0dQSjyte5nGRyE0IQydjzsFetehXLli3F3++6C82aNFFvIS5yLYlIQeZ20j3r4CMhVKxejcKFC1G5ZUvU9dG/P/I6+3/CgCxw/cLatfiqLPr95WdBOVuJAzi0gsqGyLRFfSqL4gcAQxFdmd2eUvoXAPsU+iEYyOIgsvMneNguqEtumx08zawSr+4mHSnILAOZFcKTgKycihRMLlVMIjsbyMqKbiPeqlUr3DxsmHj1tclCO1u3kgNZDIPN8ysWwaBi9WoMmzEDG7dsSZzxMHMmKlavtmvLAF63ErdFxRtvoN3IkYgMGYJ2I0dmpGUSyKI+D1aFiihGA6gH4AEAowkhHVy3UkvRu3dRtcB1Xl6Z0W6tsiC3jmh4scmK7EywJEyIQ1bO5jIlEOd67rmn8dZbb6JozBg0qFfPm8sJsHMnifJU8GKRcDqFCxeKB+aFC/V1WcJkULR+w5Z86dgN/eIEOGtWxpGFbHX4ZaedlhaSkxIFpXQqpfS3SOy/9E8AxxBCRhFCfpWKztUEJEjCucxH3CuvLEp6+9cRjYwkli2Lrsh+/PEq9OyZnxSzMPsb5HLblxSVvsxa4NtTWRs6cuBlKvKg9CDGjr0Dv/71r3HNwIF6d5OINADxXeaC4nVYpDDOULlli1ju0fVRsXo12o0YgcigQfE4g8nur369YR8qJ9LlnXcepl97bXKMpksXzF61yhvJubQqtMHs2KZ9kyilpyC6x1MTcFNaQ8jhDEyivZfMyvMuour/aNmA7VglbEwiL68MPXtGrZJIhGi/N7aWhMxVxNYlcyGJdPy2ENitOHSupwUL5mDdunUoLS5GNiFylxP7zMcoRJYGYE4CIsJIAXJzcsRyD66PitWrhUeWAsD0oUPRonHjuC4BkD93bjxwfVnHjvqtxHkIvryZ5vtXIe+88/DVlCmomj0bX02Zghc++MAfknNBFlKiIIQcRwhJilxRSj8G8CKAXtYt1UI4/4+nnipKeoP3snBO1oYMMqukTx9920FbEjqSELXn1u1kuwp73769mDBhPM4680xcecUV7l1OgF7uwG3MIQCU9u8vHpj793ddpy7OsIfJ27VvH7b+/LNyFpTVdh4xuLFMMiWm4SvJWf6AVQvupiIap+CxHUAZgN9btVTL4PwfRCfRuV84J25Dr5esGInoC/oZk5BZFSIdG4tE53qSzWxSkYVjbTz00IOorKzEow8/DMKSgq3LySY2kUFwZjf5OetJNdCJSIQFPwvKLWw39DM57S5V8H0RYtTVYaSqIoojKaUf8UJK6Ue15aQ5N6g+ACZmHfErtGXTXB3ykN1F7bjtn42OmzZVJGFrUdi6o9y6qHbt2oE77yxF90suwcVduyZcS7wrycuaCbfWRAqJJ69z5ygx+FSfaqAzeSv2wz2Ud955WP3f/2L6q6/iYFUVsiIRDOnSRTroq2IaqSaKQHatNSQLFVE0U+Q1sO1PTYdqEHXIgl1hLSMJ52zrhg2bYvfu7RgwYArmzRuJBg2aYs+e7WjUKHlFtenqci8EIcsztSRMyUL17MalRIj+RDtRbGLatHuwdetWTCopEcclZLu/6shDBK+xiQy0RmRQnX9d+OSTQhJh4erNOWHaA4haCLNXrcLBqioAwMGqKsxetQqdjz9eOPBnUkxDt/25axgMDiqieJcQ8idK6SPJdZLrAazx1rOaA5MB2HThHOumys3tiMrKtfjss9dQWbk2nu7ZM+GuckhFtz1HKkjCpozsLnr2concTDrrYsuWH1BWdg/69e2Ls04/Pbo1hyg4LSMEQD5jyYs1cQgRggxxd5bkjAmeRFhoA9eGsLUQMm3PqXTtWqsiihEAlhBC8pAghrMA1AXQO+B+HRLQDZb8fk3OYUTs/k2iQDOAeH5l5dr4nZ295JAKe/gQOxWWUmoUi9D9HV4sCV3dfBmZRcHKTEhBRhIm1513TsTevXsxcfx4d+dE+GFNpAMW/mov5ePuLBaUJlxC//oXDlZVgQBoVL8+du3d6+uhRbYWQm08pEgEKVFQSr8HcB4hpBuA38TEz1NKX0lJzzIUpgMf+7YfXTg3HJQCS5YUS0+ii9Zf3U3lQEYqou05vAasTa0IFSnYxiRMiEK1FkJHDKx7SeR2qqz8Eg8//BCuGzoUv+rQwWxXWJPtxAE5aYisiRQSSMWqVYmAdU5OYNt0KEEIKlavxuyVK+MuIQqgqqoKc2680RNB8JsNHt6oEbbuqn5Ap8xCCMzdc4hBtddTfQA3AjgOwEcAZsbOkKiVMCUIQDzTae7cEVixYhp69MiPy3iSYC0QEebOHREjgcSs5ry8sjhJOOkgSMJGz8SiEJECLxfpiQhDZ2mYEkpJyXhkZWVh3B13qIPXJlNiHZjK2DxV2kdUrFqFYTNmJGb0bNmCYTNmAEAyWZhaGx6skiAOKhLtD1U3Oxt1srKw/+DBuJ7OQjjkDikKACrX02wA+wGsBHApgBMRdUfVOtiQRFTffqbTU08VYdeun0AIsHz5tHhMom7dhti3bzeaNj0Ky5dPw2efvY4zzrgcV15ZhLlzR+CLL95KqqeiogADB7o/uc4PS4KXm8YmTEnCJBZhe33yyUeYN28ubvvLX9Dq6KOj24W73bNJ53LKEKi26Ui1VRFE0FhEPvsOHECLxo3RuF69Wm0h2EJFFCfFVmODEDITwNup6VJmwJYcqpc3n+lEKcWePduwYsU0dOjQCT16DMf69W8CAPbt2402bU5DdnY9bN/+HSor1+KEEy7A3LkjsHz5NABAjx7DMWjQVMydWxC3LkRkofubvJKEqWXg3HWuJ53loHMzOTq8u0nkdho37g40bdoUo0aOVG/6ZzoN1rl7DWKzMNXTvdnH8qXbdIjkXmMYGgQRNJaRzI8//4wt99/vut7aCBVR7HceKKUH/D6rIVPh15+pm+nkzFyKDnrJ8QbHSuje/VYQQuKE4MBJd+jQCR06dMKgQVNBCMHAgarYh7yvJi4kXdqWJFhdnQUhkumC17Jn0fXvf6/C888/h8mlpTi8WTM9Mci2EBfdRcgQSyM3JwcbBaQg277DCC4JRTV11i2U5ON8sTLkf5HpkJ5wRwg5CMCJ+hBE107sjj1TSmmTlPTQEjYn3LHwkwf5LcIbNmyK9957BpWVa2Pbjk/B/Pkjq01jpZRi8OBE/OHxx6OBPVbGYs6cqljfE50XraPw29XEp1Uk4TzLLAkgOlizeTwhODqimAN/cp3JCXbOPbqVOMVFF52PL7/cgPXr1qFhnTpRK4I9qY49sc5xSTmWxoED0cFGtrZCFuMA5Cu5+TT/7KTZu+pZkuZjFEDsFLXrrxe7nmwGVRcDcMXq1dKps0D1wLRuJhQfowAEp+T5RBSZdGypW5AhQ6Qn3KlmPWUF16XMQRCGUvJZFNEFc85aiAYNmmL+/JFYtqw8NhMqsep6zpz8pHrmzMlXxhrmzi2IWxFs24lnXT9N/x552g9Lgk2LLjafn/XkNTaxdOnzeOON1XjwvvvQsH79BBGYbOone2ZlutlOOllQIAR5XboAiG3T4Qxw6Zj1FEPS1Fnus3BzcFF8xpKKXHxwqWXSNh9BwfrM7EyHqUWRCk8aSwLz5hUkzU5y3EYDB04FABQV/RYbNryN9u3PQVHRm5gzJx8rVtwLIOqCikQIli1LxCQIAZYtm4aePfOrxSPcEoSNq4lNm5KFSUzCxMUki0mILArHgogeQpS416kDUHoQnTp1xN69e/Hp+++jTiSSIAr+zp+BLVqt7VgWshlSbmMc/DOvw8v4Z1FaJlPJdXledDVl2xUUCN1IbVu08Lz/k6xNU7QbOVLetylTvPYqZXBlUdREpDrM4gzekQipNo21Q4dOWL58GqIn303B1q2VAIADB/YluY+aNj0KgwZNxZIlE9CzZ3QtRqNGzXHlleMBkGrxCD9iESKZzIoQ6fhlSXghDpNYxfz58/Hxxx9jQUUF6mRluV9Q58DPly5RXaK3X9kbMS/34c3Z13os20nJVhou/7ZM2uYjKNQKokg1QfDtigLbQNQyYBfLOVNihwzJiucPHBhdN9GnTxEc688hBltLQtVHlcyEYEwtCVbfLUnodFRk4Vz79+9DcfFYnH766ejXp4944z/bLcQd6NxOtQVeSYUpn7KtNFz0OdO2+QgC4ihpDQDv3khXu7JT5viZTABQUpK8hdagQVOTFtdFZ0kRLq3/O2X5IkIwefvny9q6m0yIwGTAF62yFuXzVyQCzJz5ML766itMLilBhA1G84cQ6WZAAXL3EIugLA8vsDEx04VYX0xOwfO1TYvPQHZsaU3a5qNGWhTptiCSZeJT5iil1RbLjR17ZlLaCVabno9tmu/GijAhBudZRhg60lDliUhERzAywtm1aydKS0vQrVs39Lj4YrnFoDtfQkQWNvCbMNLkFvJNlwM/C2rI+efjhbVrjWc9eYZh32vDNh81kihSDd2AfeWVRRBNW/3ii7fQoUMntG9/Dj7/fGV8ZtSvf30+Nmx4W7p4zoQIgyAJEVmY1O3V5WRjdcisEPZ+331TsXnzZkwuKUkcSmQaeAbkd+dZ5nLKFEuChWww1A2SAZOFc2wqO5No9sqVrk618wQLsqhJxMCjxrqeUgEbC5V3GTVq1Bw9egxHhw6dsGLFvcjOrovc3I7Izq6LFSvuRYcOndCz53CrYLVOx5YkdC4rUytBlC+qR0YMfGxCJdORxY8/bsE99/wdvf/wB3Q65xwxMZgGtQG928ktObg1F/0uk8r6GOiOTU0pbH7oNRShRWEJ0ffF9PAgFnxwmp0R5Ux5dfJk7Zr0TSY3tSJkz6q7DRGorAURMYgW3KlIgp82e9ddk7Br1y6UFhfLd4T1MvNJBFF8Iijrgn0DTpUrSoGK1auT12n062e0TiMjZxJlwOeZLoQWhSFkLxVPPVWEioqC+KDvBK+feqpIW48TnOYXzTmuJjZYreuXiduHlZumbUlC1C8vl5vYhMyq+PrrSjz44P24ZsgQnPjrXycHq/0gBZVOpgwwNm/GXt5OEHMfzZyJjVu3giKxSK5i9WpttbIZQ0nbb6QDtdSySAtREEIOJ4SsIIT8N3ZvLtE7SAhZG7ueSXU/E/0Qy9nDgxyycGY47d69DfxiRpk1MnduQZJs7tyCmFXirl+yPBFxiAZ3Pk9mJcjusjptLQ23pMAThFNm4sQiEEIwfswY9y4m51kW5OZ1/IbtQOVxsPcyMEp3pzVwH6V0lpMtdG9vNRDpcj3dDuBlSumdhJDbY+lRAr09lNKOKe1ZDGa/L/XhQY5FIINDEs7U2YEDo6fULV1aDkLEO8Dq+mZiQYhkOiuCfVaRg+jZlBBExGC7wE7mevrss08xZ85sjBg+HLmtWydWWJscRiTbxgMI3lIgJNqGc1fp+NGOj2W9uI+kx6Y6bivnS8a0a7sXlGcI+lBTkS6iuAJA19jzbACvQkwUKYf9C1v1Vdfmhwclps5GF9apd4DV9S8okuDvzrNJvkzPjTVhSxJsHePGFaJx48a4Y9Qoe3eTAxMrorZCQBZeF6IJj02VtOtmLyjf4AdRZzjSFaM4klL6bez5OwBHSvTqE0LeJYS8SQj5Q5AdcmtNyrYTN91Dq0+fojhJRPsRJYs+fYqM+yfLExGCiauJfdbVqxv8VXlurQ0+X5f3zjtv4umn/4nb/vpXtGje3PuMJhs3VCYOIPZvQ66aKe3fX+w+6t/f/Q9OBELSP0vKr78lQxGYRUEIeQnAUYKsQjZBKaWEENmvqS2l9BtCSHsArxBCPqKUfiFoaxiAYQDQokWuRR+NVYXgV10PHFimPTxI1L7o3GyTPqqIQyfTEURy//R3N8Qgcy/J8hxXEiHJBxOJZjc5d0IoCgtvxxFHHIERt9zizd2ksyhSDfZNVvYsStvITPIkOnH3kWJ3WuHW4i52r82IWVLODyDd34sAEBhRUEovkeURQr4nhBxNKf2WEHI0gB8kdXwTu28ghLwK4HQA1YiCUjodwHQgunusrm/+vcgku450hwfZ9iFdJGFKFjJ9EXGI8nQEokuLXEz8tWLFMrz++mu4d+pUNG7UKHF+hKlFIbIOTAaCTLcqZDAhBIvyKveRaFFd3F1kSRYZtd9SDSSMtGwzTgj5O4CtTDD7cErp3zid5gB2U0p/IYTkAPg3gCsopZ+q6j722LPohAnibcb9tg4T3wfzdRRBEYQozybttyUBRAdw/q6zKGTxCMeCyM5O3FnLgT+oKDsbAKpw7rlnYvv27fjso49QNzu7+nbhttuH89uI82nWIgEScucZkBOSyOWluque3aR1cl2epV67/Hz59txTp5q1EwNPOoDgkKJ04BAiC9U24+mKUdwJoDsh5L8ALomlQQg5ixAyI6ZzIoB3CSEfAPgXgDt1JCEDO3j5Ab4+levItg9urQhbUpARgw05yOqUkYeNJWHijuLJhb8WLVqItWvXoqS4OEoSusV1gHpmkyywLcMhNEioULFqFdoNH47IgAFol5+vXAdRsXo12uXnIzJwoFbXT3dRXufOmD50KNq2aAGCKNmknSQA/wefNKFGHlw0YcK7gVkPQZVzSxA6mVsrQiSzIQPAzpIgJGE1sM9szIFN85YFf6d0P0455UQ0atQI77/9NiIOMTiuJ7eWhOxgItUxp35bFG6eRWmZjJGbHJfqrL7euGVL9JxkppqGdeti+nXXCV1JfloUUmTa+JZp/WGQiRZFYAiCwE3rY0k3Wkb/pdD1Nx0koXrr95pvYk3IZHxMQmVNPProDHzxxReYVFKCCCHV3UKmloQodsHn83A7GAT95mnypeHk0kVzCxcCiLl8ZszAxi1bAFT/xrO6PKSzovxcVJdpb/OZ1h9D1Dii8BM2pLN4cVFshXV0RbWzmG7x4iJXdcvyRXKZTJYvswpUEJUR1eGGLGSEILI6RGX4eMbu3bswceIEnH/++bisVy+z4LWpW0nlmlLBjzdJ20HGh0GpMkYA1eQxS0BEJDJdHnmdO2P6ddclu4uuuy5+lrdvyLTBOYi32YARbgrIwc3/j1KKPXuiW3kQgvg02aVLy9GrV35ScNt0QDaR21oVNhYFfzclAyDhcpIRgi7WIAtqE1J9WqzoeuCBafjuu++w6IknEtuIi6a/mp5kJ7MoWJjELPyEYyWZyk3zOd3cnJy4tcDC2XNJRiRJujk50jzprCibfprA+XJmkusnE/skQUgUMbgl+OhAl5gWu3RpOZYuLQcA9OqVn7SWQmdB2OT5QRI6spARB5uv0nNjUaimupoEs3/66Uf8/e934fe/+x06n3tuIqYgG/RFridVbCAV4AdJL4OmSVmFTmn//sIYRWn//gAgJRKRrjX8Joug6vSKQ4AwQtcT2P9T8j9KF+hPHjhVu8AGSxJ8/aI0LzexKERtyAjClhxUJCGzLGR5bPruu+/Ejh07MKmkxM5aAMxIwtTC8AOmby9u33IM6sg7/3xMv/56tM3JASEEbXNyooHsmHtIFGdwaorrdu7s7U3MbwRRpx/I1H6hFlsU/P9k8eIi7N69LT64OzGGhg2bJW2nIftfynaBHTSoDImfjroPujyVxaBLmxCDiWXAp1m5rctJRwKifZtUOv/73ybcf/+9GDRwIH5z0knVz71243IC1MQgi1fYkojpm64fb8SWdeR16SKNGzjywoULUbllC3JzcqqtvnbbrudyujqBzHuLZ3+oGdS3WkcUogGY3S4cQNJWHD17JmIMOpJwYhJsjAIABg2qvpWHrC5TuRuSsCULWTkRcYjK2loUOrlqbQUhQGnpBBw8eBDF48fbWxGy4LTXH2uqfuzsYMoPrF4HWk35vC5dzFdSS+qqdsARTzZBkEWQ9fqBDOpbrSEK9dt7wm3EbxfObtinKt+wYTP06pUfJ4WoJZHYyoMnG371tg1xmJKErk7d3XmWWRZ8GVuisFlMpwuE/+c/n2PWrEdx8003oV1urnzFtMzNpMqX6aQSfg8aovpUbejat+kfp+sccJS0lcfMmQAQkkWGWD41bsFd+/bJW3ioCIIHpRSDBiXCNnPmVGlJgm1HtpXHokVRt9bgwQm31uOPF6BRo2bo27fIN5IwIQw3loTIgrCxElRuI0LEG/2JFtyp7gMG9MOyZUvxxeef44icnMTCOXaBHXsXLawTrdrWLbBzLnahneBe8fLLKHzsMVRu3ozcli1ROngw8rp2db4o1e8qd5ap9WObdiu31eH0lAvvysvdt2GLTB8LA+5frVpwB4jffHUQxRhMtgtPHkAJl0cARN1aS5eW4/HHC+IksXRp9BQ8CBblifouepMXDe5sWqZjShI6IjANMKviDX5ca9a8g8WLF+EvI0dGSYIf7EWDu8rScO62P0xJfKLilVcwrLwcG3/4AZRSbPzhBwy7915U/Otf8rpsvrypgp99Yuqy3sojqM/GdtBINdwMbD6h1rieVGBPmhPFGETbhZv8r6I6BIMHR91Q7NTZSy/Nj1sYqnpt0jrLxOYuIw1RnaaWhsrqcHtFIsDYsXcgJycHI0eMMItNyFxOMteSyuWkIpNYXuGsWdj9yy9JWbt/+QWFjz+OvG7d5OVNQIj7N003ZXVlNPn7s7KwqVUr7K1XL/7FWbZsGQ46FhiDrEgE6w4/XN5WplsAGYr6e/eidWUl6hw8aFwmJApE3/wbNWqWtO5BtV24OUkk6h88uCxOEgCqkYSoTlvSELdt1meddcKnVWRiQxI2i/BE95dffgkvv/wSyu65B00aNzZbN6EKXrMyFWnwMgUqN28Wyw0WqwkhG4z9Gvi9kI+m/KZWrXBYq1Zod9hh8e//ES1bYuOWLahiykRIdCpui8aN3fejtsDif0UpxdadO7EJwLFffmlcrlYTBTsA9ulThOQV1MR4sZyoPhaUUsyZU5Ake/zxAgweLA6Uu7UidM86C0L27KSd6a+iu4gI3JCALI4hnukUPZQoNzcXfx42TB1DkMlk1oUOFj/O3JYtsfGH6keuqFYs+wJ+wLYhAJWuST0Snb316iWRBIA4GXzz00/Yd+AA6mZno1Xz5iFJmML5LA3+t4QQtDjsMGyuX9+qiRoZo9CBf0tOyHnLwT+SePHFclx6aT7mz6/CpZfm48UXEzELVT2mJCHS4Qd9UZ5IR5RvWo8Xi8Jkqiyrv2TJYqxZswYTiopQr25dM3dT9J+ijk/I9HgYDrql116LhvXqJcka1quH0sGDjcrHYfPPN6kjaIh/ZNV+Z0CULE5t0wZnHXssTm3TJiQJN5ANbNXU7L8DtcqisP18vBCEg0gkOnWWjUk4MYtGjZopyUg1LphYErxMRAisjsyS4GV+XSJrQUUarP7BgwcwblwhTj75ZAzMy7OLQThQPXsgBh55F10EUFp91lO3bsH72T28/Xu2KkKkBxYWhilqBVG4eYnyShJsXr9+CbdWdNAjGDJE7dYyrVv0rJKprAM+n03rLlNdXRBb5p7ir8cem4X//Oc/eHrJEmQRol9lDahl/LMMHsgiHrg2dW+ZIl2DthcSMoTN6ZHpwi+//ILBN92ENR9+iBbNm+OJGTOia3nSDR8Jo0a7ngwtsWr6ukFapSPLi0SSV3bL9oCSDdQ2z3w/Uzmo6y6droggeJ1fftmDCROKcO655+L3l11WPQYhOr3OdLsOQJ3n5PMIcqD2e2C0eTPxo22XdRTddRcKxoyJu2cppSgYMwZFd93lvU8+YmZFBZo3a4b177yDghtvxKji4nR3KRkmA5sGNZIo3HwmflgQst+fye9SRRC8vqlFobuz7XghFJvLD4vigQfuw//+9z/cOWlSdBctGwIQ3fk6RLCVBwWTAV305Up1fzyCUoptO3ag/OGH42RRMGYMyh9+GNt27KgW2zPFuMmTMfWhh+LpwtJSlD/8sKe+Pv3iixhy1VUAgL6XX46XV6503b/A4fL/VStcTyqYfm46kjCVm5KGTZ4tWagIg0/rLhtd02mxqrLbt2/DXXdNxqW9euGCLl3E02EdqALWsviECiZl3AwQhNjHCLxCVLfbfpj00/JvIYSgbOJEAED5ww/HB/P8G25A2cSJcOt+GpqXhyuHDMGIG29EVVUVFixZgreXL6+md/7vfoedP/9cTf6P4mJccuGFSbJvvv0WbVq1AgBkZ2ejaZMm2Prjj8hp0cJVHwOHi8+u1hKFVwtCle8HQfBpvywK1bPpwO9Mi5WleV2btRGy6bGO7J577sZPP/2EyRMnJruVZIcSmVyA3hrxG6qB2nQQVw2+fJ7pQB0kaVkOUA5ZsG/8XkgCANrl5qLF4Yfj/Q8/xPebN+P0U05BC8GivpXPPee6jUMG7Oeo+d/WKqKw+X5lOkmYyNwQBiuzsSJEMtn6Ch1JyK7vvvsW5eVTMeDqq3Haqae6G/wduQMdIcjcV+xzproZUokArB/H3cSiYMwYz2Rx/cCBeGzBAnz3/fcYOmCAUMfGomh19NH4+ptv0PqYY3DgwAFs37FDSD4ZDc3nWWuIwvR75SdByOSmpKDLcyNTWRK8zMvlZYsOWdlJk0qwf/9+TCgqMicJ59mBKN8WmRKvALwN0LZldfo+kgUbk3DcTU4a8GZZ9P6//8O4O+/E/gMHMG/6dKGOjUVxea9emL1gAc49+2wseuYZXHT++Z6ILBNRo4kiXRaETO7VimCfTSwLlaWgS7sZyFXrIthnncuJlTnpDRvWY8aMRzDsT39Ch/btE7u9mrqXTOMYMrdTuqwG0eArG5B5uS5tW78uzyTfEIQQNGvSJCkm4cQsmjVp4mkgrlu3Lrp16YJmTZsiKyvLc1+vy8vDoJtuwnFnn43DmzXDgkce8VxnpqHGEkXQFoQsz5Q4/CAJUXkZOYjqCsKSUJGKyvWkIiBCgPHjx6Ju3boYe8cdZrEFUfDaxuXElpWl/YDpoJ9qpLt9AEWjRoHfVser2wkAqqqq8OaaNXgyduaFV9SvXx9PPvqoL3VlKmrs9Fg/9PwgCdmArEqrnlUEILuz+qaXmzIqkrDR49Nr176PJ55YgIIRI3DUkUcGG5PgcajFH7x8qVPVB6uqiDJti08//xzHnX02Lj7/fBzfoYOnumoTaqxFoYLJd83UMlDJVVYDn5blqXREeioi4eUyPVldfl2mU2Kda+zYO3D44Yfjtr/8xS4mIctn05nidnLe4GP3ildfReGcOYlzqAcORN4FF1TXN63Xq44X/QzCSb/+NTasWZPubhxyqJEWhQz8ICnLlw3IpnKZTJavG7jZ56AsANMdX/kyulXWIj2T9hzd1157FcuWLcXoUaPQtEmT5MFctVOsDRGIINMTWSk+o+LVVzHs/vuxcfNmUEqxcfNmDCorA+ndG+2GDUPFa68lF/DjDd7kTcemTd2PLcQhhVphUXixIGR5ppaFW6tCRxYmd14mIiSdReGVkGxmP/G6QHQb8datW+PmP/9ZH4OQ3UUWgi2ByPICIIzCOXOqHXTktLJx82YMe/BBAEi2MGxgYxGodA9hyyKEHWq0RWHyUqPSUVkRJrpurAoTi8L0bjKw63Rt6vKbPJ599mm89dZbKBo3Dg3q17e3GGSuJhUyYODTHWi0+5dfUFhRIVfQfRFTidCqqBGokUQRJEGIfoM6mWzQFeWxcjZfd9no2gzgMheSrh7ZJn8mM58iEaCq6gDGjLkDJ5xwAoYMHmzvWpJZCiZWBK/PPwfsfjI50KgamdgOyLIvuKmuTbshWRzySAtREEL6EUI+IYRUEULOUuj1IoR8TghZTwi53Z+23X3vTX9XNi9zKmLQleXrURGOjrRsL1VZtzOd+Gvu3DlYt24dSktKkJ2VZW9J2FgTugC2LSF4JJDSQYOqHXTEQ0smNl9Er6jlRPD6G2/gjG7dkH3kkVj0zDPp7k4gSJdF8TGAKwG8LlMghGQBuB/ApQBOAnA1IeQkL43qBm1TK0JUl2pAlqVN6jcduPl2ZO3LSEV0T8cFAL/8shfFxeNx9tlno/cVV5i7mGQEYYp0zHYSIK9rV0y/+Wa0bdkSAMB/PRrWq4fSvLxgBmibNyU/6tWgYtEitOvYEZGWLdGuY0dULFrkrR8BILd1azx2330Y0KdPursSGNISzKaUrgMAzZzocwCsp5RuiOkuAHAFgE9t29N9R21+G6YyWb5JWyLLwKR+E8uBz/diPdhaCTrXlnM9/PCD+Prrr/HYo49W30bcSbNykcw0gM3DLTn4TCp5Xbsir2tXgFJUvPZa8lTZvDzkcfsNuQIh/vXbpC5LsqhYtAjDCgqwe88eAMDGTZswrKAAAJDXt6+rbo6bPBmHN2+OETfeCCC6zfgROTnIv+EGV/UBiB9SFInUSE8+gMye9dQKwNdMehOATiJFQsgwAMMAICcnl5GrG/BKECK5Kq17VslMLQD22SQtG9BVeiZbhIumyOrKRyLAjh3bMXlyKbpfcgku6tat+jbibk6y081mMrEkbGdEiWQuB+a8Cy+MEoMoNsLWqatflG8qc1u/BxROnBgnCQe79+xB4cSJrokiiG3GawMCIwpCyEsAjhJkFVJKn/azLUrpdADTAaB9+7OoyYuL31aEV8JwSxaidmSWBZv205Lw8yoruwdbt27F5NJSsxiEKjbBykwsCZMBPxXuKL8GXJ8H7sDqlKDym2+s5CYItxl3h8CIglJ6iccqvgHQhkm3jsmU8NOKkMndkgKbjv7eKIDEEanJ+9pU1xfdVTI3BJFOkvj+++9RVjYF/fv3x5lnnGFGEKLB360lwednKrwM1qZl3bZhUs6w7txWrbBx0yah3Av83ma8NiCTXU/vADieEHIsogRxFQDxf1UDFXkERRB8mn+eP78Iu3Ztw/XXlwEgAChmzixAo0bNkJdXZEwU/ODOypy0zpXEX6KDiUzL8fq6qbHsNXlyKfbu3YuS4uLqJCBzOQHm7ieVxcGm+WeRyycIIhENoLYDNq/v48Dtug8uUTpmTFKMAgAaNmiAUu6MClv4vc14bUC6psf2JoRsAnAugOcJIcti8mMIIS8AAKX0AIBbACwDsA7AQkrpJ/Zt2ecFTRKUUuzatQ3PPluOGTMKAFDMmFGAZ54px65d26A7b1dGEplmIdhcX365AQ8//BCuGzoUvzr+eDNrwY0lwerxMtFzEFB9KQ8F+NF/gzry+vbF9LIytG3dGoQQtG3dGtPLylzHJxw424z3v+IKX7YZf+e999D6lFPw5DPP4Ia//AUnd+7suc5MA8nYQ8BdokOHs+ikSe8C8IckbElDlle9nig5PPtseVxy+eX5GDasDISQauVMLQodOQAJi4G96ywEVk92z8qqLmctCjaflw8ZMgiLFy/CF//5D445+uhEENs51tSxKFi5KG1711ksOsvEJO3lLpPJdHR6orRK7tY9J8lbd9xxOLF9e7M6AkJVVRXOuOgiPDlzZq3dQXbdhg048bPPkmRk8OA1lFLhurYaOZ+LHRRFclWeqg6Vji6Pl0ciBH/6U1lS/SYk4fVyM2OJH9x1mwCK6lA9f/TRh5g3rwL5t94aJQnbAdvN5cDEAkkHTN5Q/Kzbja4uz6SdFFtX4Tbj7pDJMQpfYWpBmMpMfz8ysqCUxtxOCTzySAGGDSsDICcL51llWbBpma6obDouABgz5g40bdoUo267Tf/m7kA1wJu4oEzcUyZ6mYjoF0yelsnShRT2Jdxm3B1qpEXBwwtJ6CwNUVpXv0MSzzxTjssvz8dzz1Xhiivy8fTT5Zg+vQCI7RUqG8j5tk0IwpYcTHT8uFatWonnn38et48ahebNm9u5dXiZLPYgIw1bV0sq4PYN24838xS/3Yc4dFCjLQpbq9nGauDzZWXF1gVBo0bNkmISUUsCaNy4GUTuJ75OlVyXpyIRv0hC1J/qhEUxZsxoHH300bj15pvF1oOKLGwHdZ0VoiubCTB9+zaxKvxoU1ev2/6GyCjUSKKwsSBkch1pmFgRsgEbAAYOLAKlFJEIickJbrihjEmLy5kM9E7aZIW1KN/mDAlRWd3Rpk76hReex+rVq/HQAw+gYYMGyUFr0wCybUzCKa96VgWW3SKIgTqIwdVNWyFZ1HjUSKIQwZQkvKRVxCAuSzTpxF1Wn6mOjUXhpyUhu6qqDqKwcDSOP/54DL3mGjsXkwgmMQkdQdjCz0B4OgbJcGAOYYhaG6NgBzaZno4UTCwJXuZm0OXblLXt56DvRx2yugBg/vx5+PjjjzFxwgTUqVPH3MVkOovJZBB3M0geYgNrxeuvo90NNyDSty/a3XgjKl6XbtjsDbI3MdN8Wz2fUHTXXfjHfff5WmfXyy/Hu++/72udKnxVWYnfdOkSaBs12qJQDdg2Mlsrgn1WyVR3WZ4oX1be6+VmRbeJfP/+X1BUNA5nnHEG+vbpY7eGATB3MbG6sjQrZ+8iBOGC4tOm5URyTqdi5UoMe+ih+LGqG7dswbCHHgLAHKMqqtekLZs+svkmMP1MQihx4MABZGf7M8TXSKKQfR9NSMImbfLMy0zvfFnZmzn7zOd7eev3m3TY65FHpuOrr77Cww88gAg7KJhaECqo4hIqgjCV+wWTwdDjgFlYUVHt7O3d+/ahcN48/XnbbsnCAiPuuANrP/5YruCinY6nnIKppaVKndIpUzB7wQIc0bIl2hxzDM487TQAwBdffombR43C5i1b0LBhQzxSVoYTjj8e3//wA27861+xYeNGAMCDf/87zjvnHGUbc558EtcXFODAgQN4dNo0nHPGGfjxp58wdPhwbNi4EQ0bNMD0KVNw6skno+iuu9C4USP89ZZbAAC/6dIFz82bBwC49I9/RJdOnfDGO++g1dFH4+k5c9CgQQOsWbsWQ/PzAQA9unaNt/vY/Pl46rnn8POuXTh48CBee/ZZ689QhFrjehIN+LakYEMSsjd90zs/YMv6oBvkTeo3JQ63ZVj5zz/vxMSJJejWrRu6X3KJucvJgaklwerK8vjnmvAWy3xBZGdvV27dmpL2M6quGNasXYsFS5Zg7auv4oUFC/DO2rXxvGEjR+LeyZOx5pVX8I/iYtx0220AgOGjR+PC887DB6+9hvdeeQUnn3CCtp3du3dj7auv4oG778bQ4cMBAOPvugunn3IKPnz9dUwqLMTgm2/W1vPfDRtw83XX4ZPVq9GsaVMsjg381w4fjnsnT8YHr71Wrcx7H36IRbNm+UYSQA21KFiYftd0pCHLE1kCqnzTOy8zIQzTdlXPXi0JWVlHXl5ehs2bN+POSZNACKnudnJgQwgifT5P92yinykwfKvPzcnBxs2bq8tbtHBVn9+YOmmSXsnnfq188030vuwyNGzYEABwec+eAICff/4Zb7zzDvpdd11c95eYNfbKqlV4/IEHAABZWVlo2qSJtp2rr7wSAHDBeedhx86d2LZ9O1a99RYWz5oFALjogguw9ccfsWPnTmU9x+bmouMppwAAzjztNHz19dfYtn07tm3fjgvOOw8AMKh/f7z48svxMt27dsXhzZvrPwyL/3uNJQrRQG8i01kVps+6gd+527z969KOTLT7q0jO57NpNyShKheJAFu2bMY99/wDV/bujXPOPrv6oUQ2sQg3sQnRswq2g1SGkUppXh6GPfhgkvupYb16KBVtrS0aNGQDiW6A8ZN4UkRiVZSiWZMmWPvqq77Up5vRyCI7OxtVVVXx9F7m/1WPOTs9KxLBngMHtG03ipGgYUejd81nXCNdT36QhOwNXfYsGvRFdZvqyKwEG0tCRkQ6i8Jva8LJmzx5Enbt2oWJEyZ4sx5EsLUmMhmyL5ulTt4FF2D6n/+Mtjk5IISgbU4Opt94o90xqibtBw0f27rg3HPxzxdfxJ49e7Bz5048u2wZAKDJYYfh2LZt8eTTTwMAKKX4IBY/ufj88/FgzBI4ePAgtu/YoW3niX/+EwCw6s030bRJEzRt0gTn//a38TO/X121CjktWqDJYYehXZs2eO/DDwEA733wAb6MxUJkaNa0KZo1bYpVb74JAP6cI675jGskUfBwY0nI8nS/T7fWgohA/Big3ejbtm9Sd2XlRjz44AO49pprcOIJJ+iDzl6tCVWaldWk+IQAeRdcgK8efhhVixbhq4cf1gexbWD6o8kgnHHaafjjH/6A0y68EJdedRXOPv30eF7FQw9hZkUFTrvwQpzcuTOefvFFAED5pEn416pVOOX883HmxRfj088/BwBcdtVV+N+33wrbqV+/Pk7v1g03/vWvmFleDgAo+tvfsOaDD3DqBRfg9pISzI5Ny+3z+9/jx23bcHLnzrhv5kz8ymCzwlnTpuHmUaPQsWtXqHYAV/XRBjVym/G77xZvM25LGDbPtjLVm70qX+Q+4tO6S7RtOJ/WyUV30Tbjzv26667FggXz8d/PPkObVq2iA7Nq62/V9uGA/bbhbDm2PE84Tr6JO8xJs3LRswMZOal03TybpE11VHJdniB/XYcOyduM26KGjVfpwroNG3BijPAckEGDwm3GdSQhepu2eeZlpm/uMjLRWR1uCMLWCrCtn/8snevTTz/BnDmP45abb0abNrHTbW2tBBYqXZPYhJsBMdOhe4M3ecOX6WSodRAidahxRGFiNchkojwTS0JEDiZ9VNXBy2Vp07sNAfDtuSEetvzYsYVo3LgxRo8a5d3dJIKOXNwQgFfScDO4+jkgu63LtpwfBOVneyECQY2d9QS4syp0ZU2Iw/ausxR0aVuiMNks0C1J8OXffPPfePrppzGxpAQtDj9cP/ibkoLOmjCxTGTIBMuCkOD7YdOGSldXD5dPKYVqFpCnvoTQwk24ocZZFA50JCF7OxfVoyIJ04FaVEZVt0laRxYyclHVYaOjKgtEv5B33HE7jjzySOTfequ9xWCjw8t4mMQFDgWYfGkzGPV/+QVbd+50NViF8A5KKbbu3In6e/dalauRFoWJJSHLVxGGSseWLHiZDWGYEINJvgmpeLmWL1+G119/Hffdey8aN26sDhR7dTXJ9L3IagP8ejs3rKf1t99iE4DNzPqAEBr4/N2sv3cvWm/aZFWmRhIFCxPLws2zSmZqUYjqEQ3UsjRfzvRtPxUkQWkV7rjjdrRv3x5/uu4680A1D1trgs3jn01n+djqeIVukDUdzIOuxwupxMrWOXgQx6oGKTffDzfw8/+aiu9Iml9kaixR2BAEmzaReyULm2c2bWMpmFy2u8PaXE888QQ++OADVMyZg7r8NuJuLgcyGZ/Py1UwIRub+rzAdjBW6YvygiQLP4jK+WL7ETsJolzQdaWzDQVqJFHYkASfp9OREYAqTzbAs+VEA78oLarXlDx07flBPoQA+/fvw7hxY3Dqqafiqj/+0burCdDHIkzK82VM2pDVY4o0/8BrBWoTWQBp+T7VSKJgoSIFt4Qh0lWRhEhHRxiqtGxwl9WtKmNDEqaYOXMmNmzYgOeffRaRSCQRmzCBbfBaF5swtRZqCvgBy/QtPlUDZqYR56FGFqlsh22yps0+IIRsBrAx4GZyAIj3cK49CD+D8DMAws8AqDmfQVtKaUtRRo0jilSAEPKubKl7bUH4GYSfARB+BkDt+Axq7DqKECFChAjhD0KiCBEiRIgQSoRE4Q7T092BDED4GYSfARB+BkAt+AzCGEWIECFChFAitChChAgRIoQSIVGECBEiRAglQqJwCULI3wkhnxFCPiSELCGENEt3n1INQkg/QsgnhJAqQkiNnh7IghDSixDyOSFkPSHk9nT3Jx0ghDxKCPmBEPJxuvuSDhBC2hBC/kUI+TT2G8hPd5+CREgU7rECwG8opacC+A+A0WnuTzrwMYArAbye7o6kCoSQLAD3A7gUwEkAriaEnJTeXqUFjwHole5OpBEHAPyFUnoSgN8CuLkmfw9ConAJSulySumBWPJNAK3T2Z90gFK6jlL6uV6zRuEcAOsppRsopfsALABwRZr7lHJQSl8H8GO6+5EuUEq/pZS+F3veCWAdgFbp7VVwCInCHwwF8GK6OxEiJWgF4GsmvQk1eIAIoQchpB2A0wG8leauBIYavymgFxBCXgJwlCCrkFL6dEynEFEztCKVfUsVTD6DECFqKwghjQEsBjCCUroj3f0JCiFRKEApvUSVTwi5BsDvAFxMa+iCFN1nUAvxDYA2TLp1TBailoEQUgdRkqiglD6V7v4EidD15BKEkF4A/gbgckrp7nT3J0TK8A6A4wkhxxJC6gK4CsAzae5TiBSDEEIAzASwjlI6Jd39CRohUbjHfQAOA7CCELKWEPJQujuUahBCehNCNgE4F8DzhJBl6e5T0IhNYLgFwDJEA5gLKaWfpLdXqQchZD6AfwP4NSFkEyHkunT3KcXoDGAQgItiv/+1hJDL0t2poBBu4REiRIgQIZQILYoQIUKECKFESBQhQoQIEUKJkChChAgRIoQSIVGECBEiRAglQqIIESJEiBBKhEQRIkQMhJCDsWmOHxNCniSENIzJjyKELCCEfEEIWUMIeYEQ8ium3AhCyF5CSFOf+nEMIWSRhX4pIeRrQsjPfrQfIgSPkChChEhgD6W0I6X0NwD2AbgxtrBqCYBXKaUdKKVnIrpT8JFMuasRXYh3pR+doJT+j1La16LIs4huVhgiRCAIiSJECDFWAjgOQDcA+yml8QWVlNIPKKUrAYAQ0gFAYwBjECWMaiCENCaEvEwIeY8Q8hEh5IqY/OzYeSb1CSGNYuca/IYQ0s4554EQcjIh5O2YpfMhIeR4vn5K6ZuU0m/9/gBChHAQ7vUUIgQHQkg2oudNLAXwGwBrFOpXIbrV+EpEVykfSSn9ntPZC6A3pXQHISQHwJuEkGcope8QQp4BMBFAAwBzKaUfx3YjdXAjgHJKaUVsy5AsP/7GECFsEFoUIUIk0IAQshbAuwAqEd3LR4erASyglFYhukFcP4EOATCJEPIhgJcQ3ZbccV1NANAdwFkA7haU/TeAOwghowC0pZTuMf9zQoTwB6FFESJEAnsopR1ZASHkEwDCeAEh5BQAxyO63xcA1AXwJaL7gLHIA9ASwJmU0v2EkK8A1I/ltUDUdVUnJtvFFqSUziOEvAXg/wC8QAi5gVL6its/MEQINwgtihAh1HgFQD1CyDBHQAg5lRByPqLWRBGltF3sOgbAMYSQtlwdTQH8ECOJbgDY/IcBjEX0PJO7+MYJIe0BbKCUTgPwNIBT/fzjQoQwQUgUIUIoEDtnpDeAS2LTYz8BMBnAd4jGJ5ZwRZbE5CwqAJxFCPkIwGAAnwEAIWQwooHyeQDuBHA2IeQirmx/AB/HXGK/AfA430dCyN2xXXwbxnZyLXL794YIIUK4e2yIECFChFAitChChAgRIoQSIVGECBEiRAglQqIIESJEiBBKhEQRIkSIECGUCIkiRIgQIUIoERJFiBAhQoRQIiSKECFChAihxP8D/IkzG0jR0t4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create grid\n",
    "mplot = 500\n",
    "xspace = np.linspace(x_lim_lower, x_lim_upper, mplot)\n",
    "yspace = np.linspace(y_lim_lower, y_lim_upper, mplot)\n",
    "xmesh, ymesh = np.meshgrid(xspace, yspace)\n",
    "\n",
    "# mu and sigma2 (evaluated at grid)\n",
    "mu = (M @ w)[0] * xmesh + (M @ w)[1] * ymesh\n",
    "sigma2 = MSMT[0][0] * xmesh**2 + MSMT[0][1] * xmesh*ymesh \\\n",
    "            + MSMT[1][0] * xmesh*ymesh + MSMT[1][1]*ymesh**2\n",
    "\n",
    "# Predictive distribution (evaluated at grid)\n",
    "pred = norm().cdf(mu / np.sqrt(1 + sigma2))\n",
    "\n",
    "\n",
    "# Decision boundary\n",
    "vec = M @ w\n",
    "yplot = - (vec[0] / vec[1]) * xspace\n",
    "\n",
    "# Plot\n",
    "plt.pcolormesh(xspace, yspace, pred, alpha=.5, vmin=0, vmax=1, cmap='bwr')\n",
    "plt.scatter(pca.fit_transform(X)[:, 0][y == 0], pca.fit_transform(X)[:, 1][y == 0], \n",
    "            color='black', marker='x', label='y = 0')\n",
    "plt.scatter(pca.fit_transform(X)[:, 0][y == 1], pca.fit_transform(X)[:, 1][y == 1], \n",
    "            color='black', label='y = 1')\n",
    "plt.plot(xspace, yplot, color='black', label='dec. boundr.')\n",
    "plt.xlim(x_lim_lower, x_lim_upper)\n",
    "plt.ylim(y_lim_lower, y_lim_upper)\n",
    "\n",
    "plt.xlabel('PCA axis 1')\n",
    "plt.ylabel('PCA axis 2')\n",
    "plt.legend(loc='lower right')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}