0_Principal_Component_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2fab4b4d",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#PCA\" data-toc-modified-id=\"PCA-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>PCA</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a555dffb",
   "metadata": {},
   "source": [
    "## PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "58e49aad",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-06-21T07:54:55.421450Z",
     "start_time": "2023-06-21T07:54:54.288311Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Eigenvalue</th>\n",
       "      <th>PCA</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.219011</td>\n",
       "      <td>PCA1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.235273</td>\n",
       "      <td>PCA2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.443738</td>\n",
       "      <td>PCA3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.790786</td>\n",
       "      <td>PCA4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.385636</td>\n",
       "      <td>PCA5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.091869</td>\n",
       "      <td>PCA6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.985803</td>\n",
       "      <td>PCA7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.813795</td>\n",
       "      <td>PCA8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.626491</td>\n",
       "      <td>PCA9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.577521</td>\n",
       "      <td>PCA10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.567163</td>\n",
       "      <td>PCA11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0.479458</td>\n",
       "      <td>PCA12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>0.362640</td>\n",
       "      <td>PCA13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.302675</td>\n",
       "      <td>PCA14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>0.297949</td>\n",
       "      <td>PCA15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>0.254269</td>\n",
       "      <td>PCA16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>0.156621</td>\n",
       "      <td>PCA17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.149583</td>\n",
       "      <td>PCA18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>0.101770</td>\n",
       "      <td>PCA19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>0.089703</td>\n",
       "      <td>PCA20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>0.069329</td>\n",
       "      <td>PCA21</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Eigenvalue    PCA\n",
       "0     5.219011   PCA1\n",
       "1     3.235273   PCA2\n",
       "2     2.443738   PCA3\n",
       "3     1.790786   PCA4\n",
       "4     1.385636   PCA5\n",
       "5     1.091869   PCA6\n",
       "6     0.985803   PCA7\n",
       "7     0.813795   PCA8\n",
       "8     0.626491   PCA9\n",
       "9     0.577521  PCA10\n",
       "10    0.567163  PCA11\n",
       "11    0.479458  PCA12\n",
       "12    0.362640  PCA13\n",
       "13    0.302675  PCA14\n",
       "14    0.297949  PCA15\n",
       "15    0.254269  PCA16\n",
       "16    0.156621  PCA17\n",
       "17    0.149583  PCA18\n",
       "18    0.101770  PCA19\n",
       "19    0.089703  PCA20\n",
       "20    0.069329  PCA21"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAEdCAYAAAAhEEITAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhHklEQVR4nO3de5gkdX3v8fcXFmTdBbyAo6zKavAScAVhMBqN7uBdTJQnmJgYPCSEjdeAIUaUxMvxiJwoSoyRJEc8XiKul8iJggpeGLzFwOxyWRGviAgoilFkEIVdvuePqpFmnJntmflV1dTs+/U8/Wx3dfXn+/v1Tvd8p6q6OjITSZIklbFT1wOQJElaTmyuJEmSCrK5kiRJKsjmSpIkqSCbK0mSpIJsriRJkgpa0fUABu211165du3arocBwM0338yqVavM77CGc+g+v40azqH7/DZqOIfu89uo0ff8+dq0adMNmbn3r92RmUvmcsghh+RScf7555vfcQ3n0H1+GzWcQ/f5bdRwDt3nt1Gj7/nzBUzkDP2MuwUlSZIKsrmSJEkqyOZKkiSpIJsrSZKkgmyuJEmSCrK5kiRJKsjmSpIkqSCbK0mSpIKW1Bna27D2xHOGWu+EdVs5ejvrXnXK4SWGJEmSlhG3XEmSJBVkcyVJklSQzZUkSVJBNleSJEkFNXpAe0RcBdwEbAO2ZuZok/UkSZK61sanBccy84YW6kiSJHXO3YKSJEkFNd1cJXBeRGyKiA0N15IkSepcZGZz4RH7ZOZ1EXEv4FPASzLzc9PW2QBsABgZGTlk48aNjY0HYMu1Nw613shKuP6WuddZt2bPBY9jcnKS1atXL/jxXee3UcM5dJ/fRg3n0H1+GzWcQ/f5bdToe/58jY2NbZrpePJGm6s7FYp4DTCZmW+abZ3R0dGcmJhodBzzOUP7qVvmPiRtMWdoHx8fZ/369Qt+fNf5bdRwDt3nt1HDOXSf30YN59B9fhs1+p4/XxExY3PV2G7BiFgVEbtPXQeeDHylqXqSJElLQZOfFhwBzoqIqTpnZuYnG6wnSZLUucaaq8y8EjiwqXxJkqSlyFMxSJIkFWRzJUmSVJDNlSRJUkE2V5IkSQXZXEmSJBVkcyVJklSQzZUkSVJBNleSJEkF2VxJkiQVZHMlSZJUkM2VJElSQTZXkiRJBdlcSZIkFWRzJUmSVJDNlSRJUkE2V5IkSQXZXEmSJBVkcyVJklSQzZUkSVJBNleSJEkF2VxJkiQVZHMlSZJUkM2VJElSQTZXkiRJBdlcSZIkFWRzJUmSVJDNlSRJUkE2V5IkSQXZXEmSJBVkcyVJklSQzZUkSVJBjTdXEbFzRFwcEWc3XUuSJKlrbWy5Og64ooU6kiRJnWu0uYqI+wKHA+9oso4kSdJS0fSWq9OAvwFub7iOJEnSkhCZ2UxwxDOAp2fmCyNiPfDXmfmMGdbbAGwAGBkZOWTjxo2NjGfKlmtvHGq9kZVw/S1zr7NuzZ4LHsfk5CSrV69e8OO7zm+jhnPoPr+NGs6h+/w2ajiH7vPbqNH3/PkaGxvblJmj05c32Vy9ATgK2ArsBuwBfCQz/2S2x4yOjubExEQj45my9sRzhlrvhHVbOXXLijnXueqUwxc8jvHxcdavX7/gx3ed30YN59B9fhs1nEP3+W3UcA7d57dRo+/58xURMzZXje0WzMxXZOZ9M3Mt8Bzgs3M1VpIkScuB57mSJEkqaO79XoVk5jgw3kYtSZKkLrnlSpIkqSCbK0mSpIJsriRJkgqyuZIkSSrI5kqSJKkgmytJkqSCbK4kSZIKsrmSJEkqyOZKkiSpIJsrSZKkgmyuJEmSCrK5kiRJKsjmSpIkqSCbK0mSpIJsriRJkgqyuZIkSSrI5kqSJKkgmytJkqSCbK4kSZIKsrmSJEkqaKjmKiJGIuKMiPhEfXv/iDim2aFJkiT1z7Bbrt4FnAvsU9/+BnB8A+ORJEnqtWGbq70y84PA7QCZuRXY1tioJEmSemrY5urmiLgnkAAR8SjgxsZGJUmS1FMrhlzvr4CPAr8REV8E9gaObGxUkiRJPTVUc5WZmyPi8cBDgAC+npm3NTqynlp74jlDrXfCuq0cPcS6V51y+GKHJEmSWjRUcxURz5u26OCIIDPf08CYJEmSemvY3YKHDlzfDXgCsBmwuZIkSRow7G7Blwzejog9gfc2MiJJkqQeW+gZ2n8OPKjkQCRJkpaDYY+5+hj1aRioGrL9gQ82NShJkqS+GvaYqzcNXN8KfDczr5nrARGxG/A54C51nQ9n5qsXNEpJkqSeGPaYqwsWkP1L4LDMnIyIXYAvRMQnMvPLC8iSJEnqhTmbq4i4iTt2B97pLiAzc4/ZHpuZCUzWN3epLzNlSZIkLRtR9UANhUfsDGwC9gP+KTNfPsM6G4ANACMjI4ds3LixsfEAbLl2uG/tGVkJ198y9zrr1uzZaP5sNYYxOTnJ6tWrF/TYpVLDOXSf30YN59B9fhs1nEP3+W3U6Hv+fI2NjW3KzNHpy+fVXEXEvajOcwVAZl495OPuBpwFvCQzvzLbeqOjozkxMTH0eBZiPmdQP3XL3HtNZzp7esn82WoMY3x8nPXr1y/osUulhnPoPr+NGs6h+/w2ajiH7vPbqNH3/PmKiBmbq6FOxRARvxcR3wS+A1wAXAV8YtjimflTYBx46rCPkSRJ6qNhz3P1OuBRwDcy8wFUZ2j/4lwPiIi96y1WRMRK4InA1xY+VEmSpKVv2Obqtsz8MbBTROyUmecDB23nMfcBzo+Iy4CLgE9l5tkLH6okSdLSN+x5rn4aEaupzlv1voj4IdX5rmaVmZcBj1jk+CRJknpl2C1Xz6T6ypuXAp8Evg38blODkiRJ6qtht1xtAD5Un5X93Q2OR5IkqdeG3XK1B3BuRHw+Il4UESNNDkqSJKmvhmquMvO1mXkA8CJgH+CCiPh0oyOTJEnqoWG3XE35IfAD4MfAvcoPR5Ikqd+GPYnoCyJiHPgMsBdwbGY+vMmBSZIk9dGwB7TvCxyfmZc0OBZJkqTeG6q5yswTI2LniNhn8DHDfregJEnSjmKo5ioiXgy8BrgeuL1enIC7BiVJkgYMu1vweOAh9VfgSJIkaRbDflrwe8CNTQ5EkiRpORh2y9WVwHhEnAP8cmphZr65kVFJkiT11LDN1dX1Zdf6IkmSpBkM+2nB1wJExKrMvLnZIUmSJPXXsCcRfXREfBW4or59YES8vdGRSZIk9dCwB7SfBjyF6mtvyMxLgcc1NCZJkqTeGvaYKzLzexExuGhb+eFoGGtPPGe765ywbitHD7HeVaccXmJIkiSpNmxz9b2I+G0gI2JX4C+pdxFKkiTpDsPuFnw+8CJgDXANcFB9W5IkSQOG/bTgDcBzGx6LJElS7w373YJvnWHxjcBEZv5H2SFJkiT117C7BXej2hX4zfrycOAewDERcVojI5MkSeqhYQ9o3w84LDO3AkTE6cB5wJOALQ2NTZIkqXeG3XK1Blg1cHsVsE9mbmPguwYlSZJ2dMNuufp74JKIGAeC6gSiJ0fEKuDTDY1NkiSpd4b9tOAZEfFx4JFUzdUrM/O6+u6XNTU4SZKkvplzt2BEPLT+92DgPsD3gKuBe9fLJEmSNGB7W65OAI4FTp3hvgQOKz4iSZKkHpuzucrMY+t/x9oZjiRJUr9tb7fg3wxcf/a0+05ualCSJEl9tb1TMTxn4Porpt331LkeGBH3i4jzI+KKiLg8Io5b0AglSZJ6ZHvHXMUs12e6Pd1W4ITM3BwRuwObIuJTmfnV+Q5SkiSpL7a35SpnuT7T7Tvfmfn9zNxcX78JuILqZKSSJEnL1va2XB0YET+j2kq1sr5OfXu3YYtExFrgEcB/LWSQkiRJfRGZc26AWnyBiNXABcDrM/MjM9y/AdgAMDIycsjGjRsbHc+Wa28car2RlXD9LXOvs27Nno3mL6bGYvKHNTk5yerVqxf8+K7z26jR9/w2ajiH7vPbqOEcus9vo0bf8+drbGxsU2aOTl/eaHMVEbsAZwPnZuabt7f+6OhoTkxMNDYegLUnnjPUeies28qpW+besHfVKYc3mr+YGovJH9b4+Djr169f8OO7zm+jRt/z26jhHLrPb6OGc+g+v40afc+fr4iYsbka9oubF1IwgDOAK4ZprCRJkpaDxpor4DHAUcBhEXFJfXl6g/UkSZI6N9QXNy9EZn6B7Z+uQZIkaVlpcsuVJEnSDsfmSpIkqSCbK0mSpIIaO+ZK/TWf00kcPcS6izndgyRJfeOWK0mSpIJsriRJkgqyuZIkSSrI5kqSJKkgmytJkqSCbK4kSZIKsrmSJEkqyOZKkiSpIJsrSZKkgmyuJEmSCrK5kiRJKsjmSpIkqSCbK0mSpIJsriRJkgqyuZIkSSpoRdcD0I5n7YnnDLXeCeu2cvQQ6151yuGLHZIkScW45UqSJKkgmytJkqSCbK4kSZIKsrmSJEkqyOZKkiSpIJsrSZKkgmyuJEmSCrK5kiRJKsjmSpIkqSCbK0mSpIJsriRJkgpq7LsFI+KdwDOAH2bmw5qqI82k5PcX+t2FkqT5aHLL1buApzaYL0mStOQ01lxl5ueA/24qX5IkaSnymCtJkqSCIjObC49YC5w91zFXEbEB2AAwMjJyyMaNGxsbD8CWa28car2RlXD9LXOvs27Nno3mL6ZG0/mLqbFUnqNha3T1/zyMyclJVq9evaDHLpUazqH7/DZqOIfu89uo0ff8+RobG9uUmaPTl3feXA0aHR3NiYmJxsYD8zvQ+dQtcx/vP9OBziXzF1Oj6fzF1Fgqz9GwNbr6fx7G+Pg469evX9Bjl0oN59B9fhs1nEP3+W3U6Hv+fEXEjM2VuwUlSZIKaqy5ioj3A/8JPCQiromIY5qqJUmStFQ0dp6rzPyjprIlSZKWKncLSpIkFWRzJUmSVFBjuwUlLc6wnwrd3tf3gF/hI0ltcsuVJElSQTZXkiRJBdlcSZIkFWRzJUmSVJDNlSRJUkE2V5IkSQV5KgZpBzWfL5/2dA+SNDy3XEmSJBXklitJjXDLmKQdlc2VpN4q2cDZvEkqxd2CkiRJBdlcSZIkFeRuQUmahceNSVoIt1xJkiQVZHMlSZJUkM2VJElSQTZXkiRJBXlAuyR1aJiD5j1gXuoXt1xJkiQV5JYrSVrGmj6LvaerkH6dW64kSZIKcsuVJGlJ8zsk1TduuZIkSSrI5kqSJKkgmytJkqSCPOZKkrRDa+MTj57PbMfilitJkqSC3HIlSVLP+YnKpaXR5ioingr8A7Az8I7MPKXJepIkqTxPFjs/jTVXEbEz8E/Ak4BrgIsi4qOZ+dWmakqSpH5aTselNXnM1SOBb2XmlZl5K7AReGaD9SRJkjoXmdlMcMSRwFMz88/r20cBv5WZL5623gZgQ33zIcDXGxnQ/O0F3GB+pzWcQ/f5bdRwDt3nt1HDOXSf30aNvufP176Zuff0hU0ecxUzLPu1Ti4z/xX41wbHsSARMZGZo+Z3V8M5dJ/fRg3n0H1+GzWcQ/f5bdToe34pTe4WvAa438Dt+wLXNVhPkiSpc002VxcBD4qIB0TErsBzgI82WE+SJKlzje0WzMytEfFi4FyqUzG8MzMvb6peA5reVdn3/DZqOIfu89uo4Ry6z2+jhnPoPr+NGn3PL6KxA9olSZJ2RH79jSRJUkE2V5IkSQXZXEmSJBVkcyVJklSQzdV2RMRDC2btMsOyvQpl7xQRO9XXd42IgyPiHiWy56j5wgazV9dzuFvBzF0jIgZuj0XECRHxtEL5Dy+Rs50a9596TiJibUQcGREPa6DOaEQcERG/W/I1UGc/JSJOj4iPRsR/1NefWrLGHLVfVSjnKRFxTESsnbb8zwpkR0T8QUQ8u77+hIh4a0S8cOo13oSI+GzBrL2m3f6Teg4bBl+Di6xxxNR7XETsHRHviYgtEfGBiLhvgfw3R8RjFj/SWfPvERGviog/r/+fT4qIsyPijRFx94J1xiLibfVr7d8j4pSI2K9gfiev51Kv5ab4acHtiIirM/P+i8wYA94L3AW4GNiQmVfV923OzIMXmf8s4F+A24HnA68EbgYeDLwgMz+2mPy6xl9NXwS8AjgZIDPfvMj8t2fmC+vrjwXOBL4N7Af8RWZ+fDH5de6lwPrM/ElEvAw4Avg48HhgIjNfscj8bcB3gPcD7y/9JeURcSLwF8AvgTcBfw18EXgUcMZi/w/qGo8HTgV+ChxS598duA04KjO/t8j806h+Lt9DdaJhqE4w/Dzgm5l53GLyh6hf4vV8MvBYYDPwu8BpmfmP9X0lXs9vB+4F7Ar8jOp942PA04HrSzxHEXHZ9EVU/y9fB8jMRf2hMPg8RMTfAr9D9Zp+BnBNZr50Mfl17lczc//6+geALwMfAp4IPDczn7TI/B8B3wX2Bj5A9Zq+eHGjvlP+x4EtwB7Ab9bXPwg8CTgwMxf9XbwRcQowAnwGeBbV+9M3gBcCJ2fmhxaZfxodvZ5LvJYblZk7/AV46yyXfwR+ViD/IuCA+vqRwDeBR9W3Ly6QfzFwb+ABVG/GD6mX70vVNJR4jm6ieoN5FfDq+vKTqesF8jcPXD8fOLi+/sCCc/jKwPUJYGV9fQVwWaH/h4cBrwe+BVwKnAisLTT+y4GVwD3r/4+96+WrBudWYA5TuQ8AzqqvPwk4r0D+N2ZZHlRvxiXm8LNZLjcBWwvkbwFW1NfvRtWgv2Xq+SuRX/+7C/BjYNe84+d0S6Hn6KPAvwEPrd8n1gLfq6/vW+LnaOD6ZmDVwJxKzeHrA9c3TbvvklJzAB4E/F39+vta/Z734AL5l9T/BnBt6fEP/iwN/Px8sb5+9xLvGU2/npt+LTd5cbdg5U+BrwCbpl0mgFsL5O+a9QlUM/PDVH9BvDsijmCG71tciMz8QWZ+B7g6M6f++vwu5Xb9HkB1MthVwBsz87XATzLztfX1kvbIzM0AmXllXbeEnw3sQrsB2K2+voIyz1Nm5lcy86TM3A84lmoLxOcj4ksF8rdl5i1UW5VuofrFS2beXCB7ys6Z+aP6+tVUv2zJzE8Bawrk/yIiHjnD8kOBXxTIh+r5eVBm7jHtsjvw/QL5KzJzK0Bm/pRq69UeEfEhqq1NizWVfRtwUWbeWt/eCmwrkE9m/h7w71QnZDwwqy3pt2Xmd+v3jcVaGRGPiIhDqH6mbq7r3kahOQDjEfE/I2Jlff1Z8Ks9BTcWyE+AzPxmZr4uMw8A/oDqfWPRW9KBnerdf/cDVk/tYo6Ie1Lm5wjg9rjj8JB9qN9LM/MnzPz9v/PV9Ov5pzT7Wm5Mk1/c3CcXUXXxv/YLMCJeUyD/toi4d2b+ACAzL4+IJwBnA79RIJ+I2Ckzbwf+bGDZzhR6kWbm1cCREfFM4FMR8ZYSuQMeWu+qCGBtRNw9q913O1H9tVvC84H31bsHfwhMRMQFwMOBNxTIv9ObVWZeCFwYEScAjyuQvzkizqRqcD9D1aB/EngCUGoX5EREnFHnPxMYB4iIu1KmyT0aOD0idueO3Qj3o/pr9OgC+VDtotgXuH6G+84skP/tiHh8Zl4AkJnbgGMi4n8Bv18g/wcRsTozJzPzV8euRMS9KfPHHgCZeVZEnAe8LiL+nHK/0KH6xTe1m/q/I+I+mfn9unHYWqjGi4GTqHdlAi+NiJupdqEeVSD/15qPzLwMuIzqkIjFegPVljCo3rffEREJ7A+U+oP1ZODiiPg61VbKF0B1jBrVlvXFOppmX89Nv5Yb4zFXVAcWAr/IzJ83lP9E4EeZeem05XcDXpSZr19k/qFUm39/MW35WuCxmflvi8mfod5dqV78v5WZJZoGImLfaYu+n5m31gfGPi4zP1Kozs7Ak6mOE1hB9YZwbr0FYrHZf5yZjb3gI2IF8Gyqv6g/DDwS+GOqLUxvK/HzG9WHLo6leoO/lOprq7bVWwfuVWirxlSjsIbqF9g1U3949EH9XFBvRZx+35rMvLahuquodq/9sIHsA4FHZ+Y/l86eVmdn4C6l32sjYk+qLYo/Lpi5OjMnS+XNUmNnqt/DW+vX90FUuwiLbZWpf789EPhWife5WWr09vXcFJurOUTE/YDnZOYbze+mhnPoPr/FGg/NzK9tf82lW6Pv+W3UcA7d55euERG71Lt8B5ftlZk39CG/CR5zNU1E7BURL4iIz1HtEhkxv90azqH7/LZqTHNew/lt1Oh7fhs1nEP3+UVqRHWah2uA6yLivLjzqUmWfH6TPOYKqPcXH0G1i+XBwFnAAzNz0edKWQ75bdRwDt3nt1EjIt46211Un7xb8jX6nt9GDefQfX5LNf4eeEp9HPGRVMfjHpWZX6bMAfNN5zfG3YJARNwCXAj8LfCFzMyIuDIzH2h+OzWcQ/f5bdSIiJuAE6jO1TXdqZm56JPqNl2j7/lt1HAO3ee3USMiLs3MAwduHwB8hOoUNH+Xiz/nW6P5jcolcD6Iri/AS4H/ojodwyupPsF3pfnOYUfKb2kOnwV+e5b7vtOHGn3Pdw47Rn5Lc5gA7j1t2X2BS4Cblnp+k5fOB7CULlSfqDiJ6iSBvwBeToGTxS2XfOewY+Q3WQO4B3DXkmNtu0bf853DjpHf0hyeSHWetOnL7wactNTzG33uux7AUr0A66jOEfJt853DjprfYo37AS9rKr+NGn3Pdw47Rv5ymEMbz9FiL35aEIiI/WLaF3Rm5hbgk8BTdvT8Nmo4h+7z26oxUKv3n3jse34bNZxD9/lt1Oh7fmk2V5XTqL6raLqfAyXORN73/DZqNJ3fRo2+5zdeIyJ2j4jnRXVm+Qupvpj7gZn5G5n514vNb6NG3/PbqOEcus9vo0bf8xvV9aazpXBhji+wpMCXjPY93znsGPktzeEW4ALgd7jj08qlD8pvtEbf853DjpG/HObQxnPU1MUtV5Xd5rhvpfmt1HAO3ee3UeOVdY3TgVdERJHv1my5Rt/z26jhHLrPb6NG3/MbY3NVuSgijp2+MCKOATaZ30oN59B9fuM1MvMtmflbwO9RnQTw/wH7RMTLI+LBi81vo0bf89uo4Ry6z2+jRt/zm+RJRIGIGKE6E/Wt3PELZJTqW+KPyEV+CWXf89uo4Ry6z2+rxgw11wF/BPxhZjbyl2nTNfqe30YN59B9fhs1+p5fis3VgIgYAx5W37w8Mz9rfrs1nEP3+U3WiIj9gJHM/OK05Y8DrsvMby31Gn3Pb6OGc+g+v40afc9vkrsFgYjYLSKOB36f6i/200v+wup7fhs1nEP3+S3VOI2ef+JxGeS3UaPp/DZq9D2/jRp9z2+MzVXl3VS7PrYATwPeZH7rNZxD9/lt1FibmZdNX5iZE8DantToe34bNZxD9/lt1Oh7fmNWdD2AJWL/zFwHEBFnUJ1Pw/x2aziH7vPbqLEcPvHY9/w2ajiH7vPbqNH3/Ma45apy29SVzNxqfic1nEP3+W3U6P0nHpdBfhs1nEP3+W3U6Ht+YzygHYiIbcDNUzepOuKf19czM/fYkfPbqOEcus9vo8Zy+MRj3/PbqOEcus9vo0bf85tkcyWpdX3+xONyyW+jhnPoPr+NGn3Pb4LHXElqTUTsBjyf6jvCtgBnlN792HSNvue3UcM5dJ/fRo2+5zfJLVeSWhMRH6A6ruvzVJ9GvCozj+9Tjb7nt1HDOXSf30aNvuc3yeZKUmsiYsvApxFXABdm5sF9qtH3/DZqOIfu89uo0ff8JvlpQUltWg6feOx7fhs1nEP3+W3U6Ht+Y9xyJak1y+QTj73Ob6OGc+g+v40afc9vks2VJElSQe4WlCRJKsjmSpIkqSCbK2kZi4htEXFJRHwlIj4UEXedZb0vLTB/NCLeuojxTS70sX0SEcfP8dy/IyL2n2feDvG8SX3lMVfSMhYRk5m5ur7+PmBTZr554P6dM3PbUhjfchYRVwGjmXlDobwd4nmT+sotV9KO4/PAfhGxPiLOj4gzqc56/KstIfV94xHx4Yj4WkS8LyKivu/QiPhSRFwaERdGxO71+mfX978mIt4bEZ+NiG9G/YWrEbE6Ij4TEZsjYktEPHN7A42I50XEZXWt99bL9q1zLqv/vX+9/F0RcXo9pysj4vER8c6IuCIi3jWQORkRp9bj+ExE7F0vPygivlznnhURd6+Xj0fE/67n+o2I+J16+c4R8caIuKh+zF/M9dxFxF8C+wDnR8T5M8x1PCJGB8b4+nreX47qu9WIiAdExH/WNV837fEvGxjLa+tlR0TEp+v696nHf++hfkokLZrNlbQDiOoEfE+jbqaARwInZeZMu6MeARwP7A88EHhMROwKfAA4LjMPBJ4I3DLDYx8OHA48GnhVROwD/ILqS1YPBsaAU6catlnGegBwEnBYXeu4+q63Ae/JzIcD7wMGd0feHTgMeCnwMeAtwAHAuog4qF5nFbC5HscFwKvr5e8BXl7nbhlYDrAiMx9ZPx9Ty48BbszMQ4FDgWMj4gGzPXeZ+VbgOmAsM8dmm/fAGL9cz/tzwLH18n8ATq9r/urLaiPiycCDqP4/DwIOiYjHZeZZ9XovAv4P8Oql/CW30nJjcyUtbysj4hJgArgaOKNefmFmfmeWx1yYmddk5u3AJcBa4CHA9zPzIoDM/NksJ/X7j8y8pd79dT7VL/0ATo6Iy4BPA2uAkTnGfBjw4aldaJn53/XyRwNn1tffCzx24DEfy+oYhy3A9Zm5pR7/5fX4AW6nahAB/g14bETsCdwtMy+ol78beNxA7kfqfzcN5DwZeF79vP4XcE+qBgdmfu7m41bg7BlqPgZ4f339vQPrP7m+XAxsBh46MJaXAK8AfpmZ70dSa/ziZml5uyUzDxpcUG80unnGtSu/HLi+jep9IoBhDtCcvk4CzwX2Bg7JzNvq4492myNjIbWmxnw7dx7/7cz+PjdMjamsqedhanwvycxzB1eMiPXM/NzNx215x4Gw0x8/03gDeENm/ssM962hmv9IROxUN3ySWuCWK0nD+BqwT0QcClAfbzVT4/DMiNgtIu4JrAcuAvYEflg3VmPAvtup9RngD+oMIuIe9fIvAc+prz8X+MI857ATcGR9/Y+BL2TmjcBPpo6nAo6i2mU4l3OBF0TELvX4HhwRq7bzmJuA3ec53kFf5M5zHxzLn0XE1IcW1kTEver/m/9LNc8rgL9aRG1J8+SWK0nblZm3RsQfAv8YESupjrd64gyrXgicA9wfeF1mXhfVpxQ/FhETVLvKvradWpdHxOuBC6L6+ouLgaOBvwTeGREvA34E/Ok8p3EzcEBEbAJuBP6wXv4/gH+O6lQJVw6R+w6q3XWb62PHfgQ8azuP+VfgExHx/SGOu5rJccCZEXEc8O9TCzPzvIj4TeA/6y2Sk8CfAM8HPp+Zn693X14UEedk5hULqC1pnjwVg6QiIuI1wGRmvqnrscwkPH2BpJa4W1CSJKkgt1xJkiQV5JYrSZKkgmyuJEmSCrK5kiRJKsjmSpIkqSCbK0mSpIJsriRJkgr6/wXT45hTHWquAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "\n",
    "# Load dataset\n",
    "indata = \"/home/compass/CPB_vector.shp\"\n",
    "outpath = '/home/compass/out'\n",
    "\n",
    "# Read data to GeoDataFrame\n",
    "df = gpd.read_file(indata)\n",
    "\n",
    "# Define and select features\n",
    "col_selec = ['kNDVI',\n",
    "             'Tmax', 'Tmin', 'Prcp',\n",
    "             'Forest', 'Shrub', 'Crop', 'Wetland',\n",
    "             'Shannon',\n",
    "             'Elevation', 'Aspect', 'Slope', 'TWI',\n",
    "             'Clay_mean', 'Sand_mean', 'BD_mean', 'PH_mean', 'SOC_mean', 'TN_mean',\n",
    "             'Inundation',\n",
    "             'Distance'\n",
    "            ]\n",
    "\n",
    "df = df[col_selec]\n",
    "\n",
    "# Normalization\n",
    "sc = StandardScaler()\n",
    "df_scaled = sc.fit_transform(df)\n",
    "\n",
    "# Convert array to DataFrame\n",
    "df_scaled = pd.DataFrame(df_scaled, columns=col_selec)\n",
    "\n",
    "# PCA\n",
    "pca = PCA(n_components=df_scaled.shape[1])\n",
    "\n",
    "# Perform PCA and transform the data\n",
    "df_pca = pca.fit_transform(df_scaled)\n",
    "\n",
    "# Eigenvalue\n",
    "eigenvalue_df = pd.DataFrame(pca.explained_variance_, columns=['Eigenvalue'])\n",
    "eigenvalue_df[\"PCA\"] = [f'PCA{i}' for i in np.arange(1, len(eigenvalue_df) + 1)]\n",
    "display(eigenvalue_df)\n",
    "\n",
    "# Plot the Eigenvalues\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "plt.bar(x=eigenvalue_df[\"PCA\"],\n",
    "        height=eigenvalue_df['Eigenvalue'],\n",
    "        align='center'\n",
    "        )\n",
    "plt.xlabel('Principal component index')\n",
    "plt.xticks(rotation=90)\n",
    "plt.ylabel('Eigenvalue')\n",
    "\n",
    "# Add grid\n",
    "plt.grid()\n",
    "\n",
    "# Save the plot\n",
    "plt.savefig('/home/compass/out/Eigenvalue.png',\n",
    "            format='png', dpi=100, bbox_inches='tight', pad_inches=1)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b151ec7f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-06-21T08:01:56.468982Z",
     "start_time": "2023-06-21T08:01:56.257175Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.24851153, 0.40256418, 0.51892666, 0.6041978 , 0.67017707,\n",
       "       0.72216815, 0.76910872, 0.80785888, 0.83769022, 0.86518982,\n",
       "       0.89219618, 0.91502635, 0.93229404, 0.94670638, 0.96089369,\n",
       "       0.97300111, 0.98045885, 0.98758148, 0.99242744, 0.99669879,\n",
       "       1.        ])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAEhCAYAAAA3RXe6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABKCklEQVR4nO3deXhV5dX38e8ygoiRQdE4gICKImACBBFaRKJ1VpxQpA5VCwgVK9rBqlW0DrWvw2PRKsUJqxa0ok/RB4eioK2CMsgsKApoqiIgERAnwnr/2DvxEDKcJGfvc074fa7rXDl7WmvdO8nJnT3c29wdEREREYnXDukuQERERGR7pE6YiIiISBqoEyYiIiKSBuqEiYiIiKSBOmEiIiIiaaBOmIiIiEgaRNYJM7OHzexzM1tYxXIzs9FmtszM5ptZ96hqEREREck0UR4JGwccX83yE4AO4WsocH+EtYiIiIhklMg6Ye7+OvBFNaucCvzNAzOAFma2d1T1iIiIiGSSdF4Tti/wccJ0cThPREREpMHbMY25rZJ5lT5DycyGEpyyZOeddy5s06ZNlHUlbcuWLeywQ7T92KhzZHv8OHKoDemPH0eOqOM3/fhj3J2v99svshzZvo/iyKE2pDb+l99/yfrN62ufxKm8F1CJrzd/DcDOO+6cdPjG1pi8Jnm1rysC77333hp336OyZenshBUDib2p1sAnla3o7mOBsQA9evTwWbNmRV9dEqZNm0a/fv2yOke2x48jh9qQ/vhx5Ii8Df36UVJSQou5cyNLkfX7KIYc21sbHpj9AOMXjq9V/JKSElq0aJHUuh+s/ACAvm37RpYDYFCXQQwpHJL0+nF8n5NlZiurWpbOTtgkYISZTQAOB75090/TWI+ISHR+/3tWzptHi3TXIRmlzp2klS2SWvf1la8Dte8kJatv27617iBBZnWS0imyTpiZjQf6Aa3MrBgYBTQCcPcxwGTgRGAZsAm4KKpaRETS7ic/Yd2O6fy/V2qrLh0kyP5OkjpI8YnsE8HdB9Ww3IFLo8ovIpJR5s4ld9ky0B+3rDF+4XjmrZpHQV5BZDnUSdq+6d8yEZE4jBzJgSUlMHhwuitpMKI+lVfWAXv1Z6/WKoc6SZIsPbZIRESyUtmRqqgU5BUwqEu1J3VE6kVHwkREJBKZeKRKR6kkk+hImIiIREJHqkSqpyNhIiISGR2pEqmaOmEiInG49VY+nDOH7umuIxT1qUIg8jsLRbKdTkeKiMThRz9ifZcu6a6iXNSnCkGnC0VqoiNhIiJxePNNmi1cmFHjhOlUoUh66UiYiEgcrrmG/R98MN1ViEgG0ZEwEZEMFNfwDiKSPjoSJiIStSeegBkzYMMGaNcumK6BhncQafh0JExEJEpPPAFDh8K33wbTK1cG0wDnnlvtprpmS6Rhy75O2NKl217YevbZ8ItfwKZNcOKJ225z4YXBa80aGDBg2+XDh8PAgfDxx3D++dsu/9Wv4JRTgtyXXFI+u2tJCbRoAb//PfzkJzB3Lowcue32t94KP/oRvPkmXHPNtsvvvhu6doUpU+Dmm7da1LWkBJ58Eg4+GJ57Du68c9vtH3sM2rQJ1rv//m2XP/00tGoF48YFr4rx33wTmjaF++6Dp57advtp04Kvd9wBzz+/9bKdd4YXXgje33QTvPLKVos7b9kCr78eTFx9NUyfvvX2rVvD448H70eODPZhooMOgrFjg/dDh8J77229vGtXOO204P1550Fx8dbLe/eGP/4xeH/mmbB27dbLjz4arrsueH/CCfD111svP/lk6NEjeF/ZH7cU/ezt9Pnnlcev4mevXJI/e80WLoQbbth2eTU/ewD89a9J/ezt8eqrlcev5mcPgMmTk/7Z6/r448HvW5kafvbYfXeYODF4n8TPXtdp07aOX83P3icbPmX67pv4y7kdgvB/XcweX3y7VfjFBzbjwbMOgOlvckOrb+n+KeT8dzlz8wA2Mefun/H45ocA+OOd89jpuy1bbd/ngG/5z9m9gokkf/bKP5Mgks+9reJH9LlXniPJn73afu6Vx6/Fz16mfe4d2KrVDz8TEXzutTnkkB/iR/W5l5dX67+55ZL52YM6/c0FIvvZq0ykpyPN7HgzW2pmy8zsd5Usb2lmz5rZfDN728wy5/5tEZEqfP7V56zetDq5lb8JOmf/bQYf797oh/mbS6vdbJ9d99HpQpGGzt0jeQE5wAfA/kBjYB7QqcI6twOjwvcdgVdqiltYWOiZYurUqVmfI9vjx5FDbUh//Dhy1CZ+0bgiLxpXlNzKbdu6gzv41DvuKH/vbdvWpcxqZdI+ytQcakP648eRI442JAuY5VX0aaI8EtYTWObuH7r7d8AE4NQK63QCXgk7g0uAdmaWF2FNIiLxuuWW4LQX0LLstFLTpsF8EdmuRdkJ2xf4OGG6OJyXaB5wBoCZ9QTaAq0jrElEJF7nnhtc37PTTrSdMgXatg2ma7goX0QaPguOlEUQ2Ows4Dh3HxxOnw/0dPfLEtZpBvwZ6AYsIDglOdjd51WINRQYCpCXl1c4YcKESGqurY0bN5Kbm5vVObI9fhw51Ib0x69tjuc/fZ5XP0/+rkKA0tJScnJyklr3g40fcEDuAdxVcFfS8buOHElpaSkL7rmnVnXVxvb2fc7E+HHkyPb4ceSIow3JKioqmu3uPSpbFuXdkcVAm4Tp1sAniSu4+3rgIgAzM2B5+KLCemOBsQA9evTwTLkFO47bwaPOke3x48ihNqQ/fm1z/OHRP7Dy25W1Goy0pKSEFol3R1ajsEUhg7oMol9hcvUA0KIFJSUlWf19yLTvcybGjyNHtsePI0e2DNcSZSdsJtDBzNoD/wXOAX6auIKZtQA2hdeMDQZeDztmIiL1ojG2RCTTRXZNmLtvBkYALwHvAk+5+yIzG2Zmw8LVDgEWmdkS4ATg8qjqERHJBC+++CIHH3wwBx54ILfddts2y9etW8fpp59Ofn4+PXv2ZOHCheXLLr74Yvbcc0+6dNl6NJ8bbriBfffdl8GDB9O1a1cmT54ceTtEpP4iHazV3ScDkyvMG5PwfjrQIcoaREQywl//yuLp07n00kv517/+RevWrTnssMPo378/nTp1Kl/t1ltvpWvXrjz77LMsWbKESy+9lFfCwUAvvPBCRowYwQUXXLBN+CuuuIIePXroaJ5IFtGzI0VE4nDwwbzz1VcceOCB7L///jRu3JhzzjmHf/7zn1uttnjxYo4++mgAOnbsyIoVK1i1ahUAffv2Zbfddou9dBGJRvY9tkhEst4Dsx9g/MLxtdqmpKSEFitbJLXuvFXzanVRfiyee45v3niDNm1+uF+pdevWvPXWW1utVlBQwDPPPEOfPn14++23WblyJcXFxeTlVT+E4r333ktOTg5HHnkkd955Jy1btoykGSKSOjoSJiKxG79wPPNWzat5xToqyCvIvEf+3Hknu7/55jazgxvDf/C73/2OdevW0bVrV+655x66devGjjtW///y8OHD+eCDD3jggQfYe++9+dWvfpXS0kUkGjoSJiJpsT3evbhPo0Z8/PEPY1gXFxezzz77bLVOs2bNeOSRR4DgsXLt27enffv21cYtO0q2ww47MGTIEE4++eQUVy4iUdCRMBGRmHRv2pT333+f5cuX89133zFhwgT69++/1TolJSV89913ADz44IP07duXZs2aVRv3008/LX//7LPPbnP3pIhkJh0JExGJyY5m3HvvvRx33HGUlpZy8cUX07lzZ8aMCW4aHzZsGO+++y4XXHABOTk5dOrUiYceeqh8+0GDBjFt2jTWrFlD69atufHGG/n5z3/Ob3/7W+bOncumTZvo3Lkzf/3rX9PVRBGpBXXCRERidOKJJ3LiiSduNW/YsGHl73v37s37779f6bbjx1d+M8Njjz0GNIxTtiLbE3XCRGQrUd+5CBl692LUHnuMd6dPp3e66xCRjKFrwkRkK1HfuQgZevdi1Nq04ds990x3FSKSQXQkTES2sT3euRi5J59kj0WLQPtJREI6EiYiEof772ffSZPSXYWIZBB1wkRERETSQJ0wERERkTSItBNmZseb2VIzW2Zmv6tkeXMze87M5pnZIjO7KMp6RERERDJFZBfmm1kO8BfgGKAYmGlmk9x9ccJqlwKL3f0UM9sDWGpmT7j7d1HVJZLttsuHX4uINEBRHgnrCSxz9w/DTtUE4NQK6ziwqwVPsM0FvgA2R1iTSNbbLh9+3RA8/TSLbrwx3VWISAYxd48msNkA4Hh3HxxOnw8c7u4jEtbZFZgEdAR2BQa6+/9VEmsoMBQgLy+vcMKECZHUXFsbN24kNzc3q3Nke/w4cmRaG66cdyUAdxXcFUn8usi0fZSJ8ePIke3x48ihNqQ/fhw54mhDsoqKima7e49KF7p7JC/gLODBhOnzgXsqrDMA+B/AgAOB5UCz6uIWFhZ6ppg6dWrW58j2+HHkyLQ2FI0r8qJxRZHFr4tM20eZGN8fecTfveqqSFNk/T6KIYfakP74ceSIow3JAmZ5FX2aKE9HFgNtEqZbA59UWOci4JmwzmVhJ6xjhDWJiKTHuHHs9eKL6a5CRDJIlJ2wmUAHM2tvZo2BcwhOPSb6CDgawMzygIOBDyOsSURERCQjRHZ3pLtvNrMRwEtADvCwuy8ys2Hh8jHATcA4M1tAcEryKndfE1VNIiIiIpki0mdHuvtkYHKFeWMS3n8CHBtlDSIiIiKZSA/wFkmhuozhBRrHS0Rke1TjNWFm1trMnjWz1Wa2yswmmlnrOIoTyTZRj+EFGscra02ezPzbbkt3FSKSQZI5EvYI8HeCIScAzgvnHRNVUSLZrCCvgFd/9mqttpk2bRr9+vWLpiDJDE2bsqVJk3RXISIZJJm7I/dw90fcfXP4GgfsEXFdIiINy333sc///m+6qxCRDJLMkbA1ZnYeUHahyyBgbXQliYg0QE89xZ4lJemuQkQySDJHwi4GzgY+Az4lGOX+4iiLEhEREWnoajwS5u4fAf1jqEVERERku1FlJ8zMfuvu/8/M7gG2ecq3u/8y0spEREREGrDqjoS9G36dFUchIiIiItuTKjth7v5c+HaTu/8jcZmZnVXJJiIiUpVp05g7bRr90l2HiGSMZO6OvBr4RxLzRDJeXUa012j2IiIShequCTsBOBHY18xGJyxqBmyOujCRKJSNaB9VR0mj2UuV7riDNh98ABqUV0RC1R0J+4TgerD+wOyE+RuAK5IJbmbHA38GcoAH3f22Cst/A5ybUMshBIPDfpFU9SJ1UNsR7TWavaTE88+zu8YJE5EE1V0TNg+YZ2Z/d/fvaxvYzHKAvxA83qgYmGlmk9x9cUKO24Hbw/VPAa5QB0xERES2B8kM1trOzJ42s8Vm9mHZK4ntegLL3P1Dd/8OmACcWs36g/hhVH4RERGRBi2ZTtgjwP0E14EVAX8DHktiu32BjxOmi8N52zCzpsDxwMQk4oqIiIhkPXPfZhzWrVcwm+3uhWa2wN0PDef9292PqGG7s4Dj3H1wOH0+0NPdL6tk3YHAee5+ShWxhgJDAfLy8gonTJiQRNOit3HjRnJzc7M6R7bHr22OK+ddCcBdBXdFEr+usv37oH1Us0OvuorSzZtZfOedkeXI9n0URw61If3x48gRRxuSVVRUNNvde1S60N2rfQFvEBwxewYYAZwOLE1iu97ASwnTVwNXV7Hus8BPa4rp7hQWFnqmmDp1atbnyPb4tc1RNK7Ii8YVRRa/rrL9+6B9lBk5sj1+HDnUhvTHjyNHHG1IFjDLq+jTJHM6ciTQFPglUAicB/wsie1mAh3MrL2ZNQbOASZVXMnMmgNHAv9MIqaIiIhIg1DtYK3hHY5nu/tvgI3ARckGdvfNZjYCeIlgiIqH3X2RmQ0Ll48JVz0deNndv6pLA0REssJNN9F2+XKNEyYi5arthLl7qZkVmpmFh9Rqxd0nA5MrzBtTYXocMK62sUVEssorr9BS44SJSIJkHlv0DvBPM/sHUH60yt2fiawqERERkQYumU7YbsBa4KiEeU5wob5ISunZjiIisr2osRPm7klfByZSX3q2o4iIbC+SORImEis921EapN135/stW9JdhYhkEHXCRETiMHEii6ZNo1+66xCRjJHMOGEiIiIikmI1dsLMLM/MHjKzF8LpTmb28+hLExFpQK6+mvYPPJDuKkQkgyRzJGwcwYCr+4TT7xGMoi8iIsmaPp3mixaluwoRySDJdMJauftTwBYIRsIHSiOtSkRERKSBS6YT9pWZ7U4wNhhm1gv4MtKqRERERBq4ZO6OvJLgwdsHmNkbwB7AgEirEhEREWngkhmsdY6ZHQkcDBiw1N2/j7wyEZGGpHVrvm3UKN1ViEgGSebuyEuBXHdf5O4LgVwz+0X0pYmINCCPP867116b7ipEJIMkc03YEHcvKZtw93XAkGSCm9nxZrbUzJaZ2e+qWKefmc01s0Vm9lpSVYuIiIhkuWSuCdvBzMzdyy7MzwEa17RRuN5fgGOAYmCmmU1y98UJ67QA7gOOd/ePzGzPOrRBRCTzjRzJgcXFoEdsiUgomU7YS8BTZjaG4A7JYcCLSWzXE1jm7h8CmNkE4FRgccI6PwWecfePANz981rULjF7YPYDjF84vlbblJSU0GJli6TXj/Lh3SJpNXcuuSUl6a5CRDJIMqcjrwJeBYYDlwKvAL9NYrt9gY8TpovDeYkOAlqa2TQzm21mFyQRV9Jk/MLxzFs1L9IcBXkFDOoyKNIcIiIimcDCs4ypD2x2FnCcuw8Op88Herr7ZQnr3Av0AI4GdgamAye5+3sVYg0FhgLk5eUVTpgwIZKaa2vjxo3k5uZmdY7axL9y3pUA3FVwVyTx6yqT9lGm5sj2+HHkiDp+15EjKS0tZcE990SWI9v3URw51Ib0x48jRxxtSFZRUdFsd+9R6UJ3r/YF/Bj4F8Hjij4ElgMfJrFdb+ClhOmrgasrrPM74IaE6YeAs6qLW1hY6Jli6tSpWZ+jNvGLxhV50biiyOLXVSbto0zNke3x48gReRuOPNLXFRREmiLr91EMOdSG9MePI0ccbUgWMMur6NMkc03YQ8AVwGxq97iimUAHM2sP/Bc4h+AasET/BO41sx0JLvY/HPifWuQQEckOBx3Epk8+oUW66xCRjJFMJ+xLd3+htoHdfbOZjSC4sD8HeNjdF5nZsHD5GHd/18xeBOYTPJvyQQ/GIhMRaVjGjuW9adPYJ911iEjGSKYTNtXMbgeeAb4tm+nuc2ra0N0nA5MrzBtTYfp24PakqhURERFpIJLphB0efk28qMyBo1JfjohIAzV0KAd98sl2MU7Y999/T3FxMd98802tt23evDnvvvtuBFXFEz+OHNkeP44ccbShoiZNmtC6dWsa1eLxZMk8O7KoXlWJiAi89x5Nt5NxwoqLi9l1111p164dZlarbTds2MCuu+4aUWXRx48jR7bHjyNHHG1I5O6sXbuW4uJi2rdvn/R2yRwJw8xOAjoDTRIS/qHWVYqISIP3zTff1KkDJpKtzIzdd9+d1atX12q7ZB7gPQYYCFwGGHAW0LYuRYqIyPZBHTDZ3tTlZz6ZEfN/5O4XAOvc/UaC8b/a1DqTiIhITD777DPOOeccDjjgADp16sSJJ57Ie++9V/OG9dCvXz9mzZpV7Tp33303mzZtKp8+8cQTKcmg09T9+vVjzpzq77sbPHgwixcvrnadZLVr1441a9akJFaiESNGpKzGKCXTCfs6/LrJzPYBvgeSP+EpIiLQtSsbDzww3VVsF9yd008/nX79+vHBBx+wePFibr31VlatWpXu0rbphE2ePJkWLVqkr6A6ePDBB+nUqVO6y6hSaWkp9957b0bXWCaZTtjzZtaCYBiJOcAKIDOeGyQiki3uvptlI0aku4rtwtSpU2nUqBHDhg0rn9e1a1eOOOII/v3vf3PyySeXzx8xYgTjxo0DgqMy11xzDb1796ZHjx7MmTOH4447jgMOOIAxY4LRlaZNm1bl9omGDx9Ojx496Ny5M6NGjQJg9OjRfPLJJxQVFVFUVFSec82aNVx11VXcd9995dvfcMMN3HnnnQDcfvvtHHbYYeTn55fHqujll1+md+/edO/enbPOOouNGzeycuVKOnTowJo1a9iyZQtHHHEEL7/8MitWrKBjx4787Gc/Iz8/nwEDBmzVMayuDbD1Eb/c3FyuvfZaCgoK6NWrV3lHd/Xq1Zx55pkcdthhHHbYYbzxxhsArF27lmOPPZY+ffpwySWXlD0tZyv3338/v/3tD4+oHjduHJddFjzx8LTTTqOwsJDOnTszduzY8nVyc3O5/vrrOfzww5k+fTonnnhieY1VtaNdu3aMGjWK7t27c+ihh7JkyRIgeOTRRRddxKGHHkp+fj4TJ06sch/XV42dMHe/yd1L3H0iwbVgHd39unpnFhERicDChQspLCys07Zt2rRh+vTpHHHEEVx44YU8/fTTzJgxg+uvv75WcW655RZmzZrF/Pnzee2115g/fz6//OUv2WeffZg6dSpTp07dav1zzjmHJ598snz6qaee4qyzzuLll1/m/fff5+2332bu3LnMnj27vENTZs2aNdx8881MmTKFOXPm0KNHD+666y7atm3LVVddxbBhw7jzzjvp1KkTxx57LABLly5l6NChzJ8/n2bNmm3VAayuDRV99dVX9OrVi3nz5tG3b18eeOABAC6//HKuuOIKZs6cycSJExk8eDAAN954I3369OE///kP/fv356OPPtom5oABA3jmmWfKp5988kkGDhwIwMMPP8zs2bOZNWsWo0ePZu3ateV1dOnShbfeeos+ffok3Y5WrVoxZ84chg8fzh133AHATTfdRPPmzVmwYAHz58/nqKOOqnIf11eVd0ea2VHu/qqZnVHJMtz9mcq2ExGRSpx3HoesWrVdjBOW6MqXrmTuZ3OTXr+0tJScnJxq1+m6V1fuOq7+fwAr079/fwAOPfRQNm7cyK677squu+5KkyZNanXt1lNPPcXYsWPZvHkzn376KYsXLyY/P7/K9bt168bnn3/OJ598wurVq2nZsiX77bcfo0eP5uWXX6Zbt25AcJTmgw8+2GrbGTNmsHjxYn784x8D8N1339G7d28guH7rH//4B2PGjGHu3Lnl27Rp06Z8/fPOO4/Ro0fz61//utZtaNy4cfmRwcLCQv71r38BMGXKlK2uyVq/fj0bNmzg9ddfL+9gnXTSSbRs2XKbfbHHHnuw//77M2PGDDp06MDSpUvLax09ejTPPvssAB9//DHvv/8+u+++Ozk5OZx55pmV7tvq2nHGGWeU115W15QpU5gw4YcTfi1btuT555+vch/XR3VDVBwJvAqcUskyJxhBX0REklFczE4ZdAF2Q9a5c2eefvrpSpftuOOObNmypXy64oCyO+20EwA77LBD+fuy6c2bN9e4PcCKFSu44447mDlzJi1btuTCCy9MauDaAQMG8PTTT5ffVADB9W1XX301l1xySfl6GzZs2Go7d+eYY45h/Pjx28TctGkTxcXFAOWdStj2Tr6K08uXL0+qDY0aNSrfNicnh82bNwOwZcsWpk+fzs4777zNNsncRThw4ECeeuopOnbsyOmnn46ZMW3aNKZMmcL06dNp2rQp/fr1K6+pSZMmlXbea2pH2fc4sXZ336bG6vZxfVTZCXP3UWa2A/CCuz+V0qwiIrJdqO0Rq1QMsnnUUUdxzTXX8MADDzBkyBAAZs6cyaZNm2jTpg2LFy/m22+/5ZtvvuGVV17Z5vRVddq2bVvj9hs2bGCXXXahefPmrFq1ihdeeIF+4RHQXXfdlQ0bNtCqVattYp9zzjkMGTKENWvW8NprrwFw3HHHcd1113HuueeSm5vLf//7X7799tut9lGvXr249NJLWbZsGQceeGB5x+uggw7iqquu4txzz6Vt27YMGTKE559/HoCPPvqI6dOn07t3b8aPH79NG9avX19lG5Jx7LHHcu+99/Kb3/wGgLlz59K1a1f69u3LE088weWXX84LL7zAunXrKt3+jDPO4JZbbqFt27b86U9/AuDLL7+kZcuWNG3alCVLljBjxowa66hLO8pqv/vuuwFYt25dtfu4PqodrNXdt4QP4VYnLAs8MPsBxi+sXS+9pKSEFitbJLXuvFXzKMgrqENlIiLxMTOeffZZRo4cyW233UaTJk1o164dd999N3vttRdnn302+fn5dOjQofw0X7LatGlT4/aHHnoo3bp1o3Pnzuy///7lp7AAhg4dygknnMDee++9zXVhnTt3ZsOGDey7777svffeQNAhePfdd8tPfeXm5pbfJFBmjz32YNy4cQwaNIhvvw0e8XzzzTfz6aefMnPmTN544w1ycnKYOHEijzzyCEVFRRxyyCE8+uijXHLJJXTo0IHhw4dvFbOgoKDKNiRj9OjRXHrppeTn57N582b69u3LmDFjGDVqFIMGDeLpp5+mqKiI/fbbr9LtW7ZsSadOnVi8eDE9e/YE4Pjjj2fMmDHk5+dz8MEH06tXrxrrqEs7fv/733PppZfSpUsXcnJyGDVqFGeccUal+7i+nTCr7M6ErVYwu45gmIonga/K5rv7F/XKXEc9evTwmsZhicu0adNq9Z9B1DmOevSoWneUSkpKanV79KAugxhSOCTp9TNtH2Vi/DhyZHv8OHJE3oZ+/YLft4TrclItU/bRu+++yyGHHFKnHA3hcTaZ3oYVK1Zw8skns3DhwkjiJyPT91FdVfazb2az3b1HZesn89iii8OvlybMc2D/mjY0s+OBPwM5wIPufluF5f2AfwLLw1nP6HFI9VOQV8CrP3s16fXj+OMpIkDv3nz50Ue0SHcdIpIxknmAd50GZjWzHOAvwDFAMTDTzCa5e8UhbP/t7idvE0BEpCH54x9ZPm2anvkmadeuXbsqj4JJvJJ9gHcXoBNbP8D7bzVs1hNY5u4fhjEmAKcCmf8cAREREZGI1dgJM7NRQD+CTthk4ATgP0BNnbB9gY8TpouBwytZr7eZzQM+AX7t7otqLltEJMuceSadV6+G119PdyUikiGSuTB/AVAAvOPuBWaWR3B9V2XjhyVudxZwnLsPDqfPB3q6+2UJ6zQDtrj7RjM7Efizu3eoJNZQYChAXl5eYeIgaum0ceNGcnNzMybHlfOuBOCuguRvCY+6DZm2jzIxfhw5sj1+HDmijt915EhKS0tZcM89keXIlH3UvHlzDqzjczKTGay1PqKOH0eObI8fR4442lCZZcuW8eWXX241r6ioqMoL83H3al/A2+HX2UAzwIBFSWzXG3gpYfpq4OoatlkBtKpuncLCQs8UU6dOzagcReOKvGhcUWTx6yLT9lEmxo8jR7bHjyNH5G048khfV1AQaYpM2UeLFy+uc47169fXedtMiB9HjmyPH0eOONpQmcp+9oFZXkWfJpkHeM8KH+D9QNgRmwO8ncR2M4EOZtbezBoD5wCTElcws70sHJbWzHoSPMtybRKxRUREqlTbI4KJD+aeNGkSt912W7XrX3/99UyZMqXSOGeddVatcicqe6B3VBIfvl2VESNGbPXIofqIqj2DBw9OWY3plMzdkb8I344xsxeBZu6+7VM8t91uczjQ60sEQ1Q87O6LzGxYuHwMMAAYbmabCcYiOyfsNYqISAMxdvbYpNf95ptvaNKkSbXrDC0cWt+SqtW/f//yZ0hW5Q9/aLijKd17771pGWMrWaWlpTz44IPpLiMlajwSZmb/NLOfmtku7r4imQ5YGXef7O4HufsB7n5LOG9M2AHD3e91987uXuDuvdz9zbo3RUQkgx19NOu6d093FdudsrEQBwwYQMeOHfn5z39edvkLL774Ih07dqRPnz7lD28GGDduHCNGjODLL7+kXbt25c+KLHvs0ffff8+FF15Y/nzKquLccMMN3HHHHeXTXbp0YcWKFQCcdtppFBYW0rlzZ8aOrbmT+vLLL9O7d2+6d+/OBRdcwMaNG1m5ciUdOnRgzZo1bNmyhSOOOIKXX36ZFStW0LFjR372s5+Rn5/PgAED2LRp0zYxhw8fTo8ePejcuTOjRo0qn3/iiSeWHy3Lzc3l2muvpaCggF69erFq1SoAVq9ezZlnnslhhx3GYYcdxhtvvAHA2rVrOfbYY+nWrRuXXHIJlR1Xuf/++7nuuuu22t+XXXZZtfslNzeX66+/nsMPP5zp06dvdUSvqna0a9eOUaNG0b17dw499FCWLFkCBNc2XnTRRRx66KHk5+czceLEbfbxWWedxcaNG2v8vtRXMqcj7wL6AIvN7B9mNsDMqv83RUREtnbdday84IJ0V7Fdeuedd7j77rtZvHgxK1as4I033uCbb75hyJAhPPfcc/z73//ms88+22a75s2bU1BQUP4cx+eee47jjjuORo0ala+TTJzKPPzww8yePZtZs2YxevRo1q6t+kqcNWvWcPPNNzNlyhTmzJlDt27duOuuu2jbti1XXXUVw4YN484776RTp04ce+yxACxdupShQ4cyf/58mjVrxn333bdN3FtuuYVZs2Yxf/58XnvtNebP3/YYy1dffUWvXr2YN28effv25YEHHgDg8ssv54orrmDmzJlMnDiRwYMHA3DjjTfSp08f3nnnHfr3789HH320TcwBAwYwadIPVyc9+eSTDBw4sNr98tVXX9GlSxfeeuutbZ5zWV07WrVqxZw5cxg+fHh5h/imm26iefPmLFiwgPnz53PUUUdts4979OjBXXfV7rmndVFjJ8zdXwtPSe4PjAXOBj6PujAREZFU6NmzJ61bt2aHHXYgPz+fFStWsGTJEtq3b0+HDh0wM84777xKtx04cCBPPvkkABMmTCjvLJRJNk5Fo0ePLj+69PHHH/P+++9Xue6MGTNYvHgxP/7xj+natSt///vfWblyJRBcG7VhwwbGjBmz1VG3Nm3alD8n8bzzzuM///nPNnGfeuopunfvTrdu3Vi0aFGl11g1bty4/Fq5wsLC8iN5U6ZMYcSIEXTt2pX+/fuzfv16NmzYwOuvv16+D0466SRatmy5Tcw99tiDdu3aMWPGDNauXcvSpUvLa61qv+Tk5HDmmWdWun+qa8cZZ5xRae2XXvrDQ4Batmy5zT5+9NFHy/dxlJIdrHVn4BRgINAdeDTKokREGpwTTuDQL76At95KdyXbnZ122qn8/Q477MDmzZuB4EHfNenfvz9XX301X3zxBbNnz+aoo47aZp2q4uy4447lpzIhOGoGwSnSKVOmMH36dJo2bUq/fv3Kl1XG3TnmmGMYP348sPVzETdt2kRxcTEQnGYrm1+xporTy5cv54477mDmzJm0bNmSCy+8sNIaGjVqVL5tTk5O+b7bsmUL06dPZ+edd056fyQ688wzeeqpp+jYsSOnn346ZlbtfmnSpEmlQ07U1I6y731i7e6+TY0V93Fckrkm7EngXeAogscQHeAJY32JiEgSvv6anG+/TXcVEurYsSPLly/ngw8+AKjyj29ubi49e/bk8ssv5+STT96mI1BdnHbt2jFnzhwA5syZw/LlwWOSv/zyS1q2bEnTpk1ZsmQJM2bMqLbWXr168cYbb7Bs2TIg6Hi99957AFx11VWce+65/OEPf2DIkCHl23z00UdMnz69vKaKp/DWr1/PLrvsQvPmzVm1ahUvvPBCtTVUdOyxx3LvvfeWT88NH0zft29fnnjiCQBeeOEF1q1bV+n2p5xyCv/7v//L+PHjy48u1na/1LUdFWtft25dtfs4SslcE/YIQcdrmLu/6u5batxCREQkgzVp0oSxY8dy0kkn0adPH9q2rfqpngMHDuTxxx/f5lRkTXHOPPNMvvjiC7p27cr999/PQQcdBMDxxx/P5s2byc/P57rrrqNXr17V1rrHHnswbtw4Bg0aRH5+PkcffTRLlizhtddeY+bMmeUdscaNG/PII48AcMghh/Doo4+Sn5/PF198wfDhw7eKWVBQQLdu3ejcuTMXX3xx+enAZI0ePZpZs2aRn59Pp06dGDNmDACjRo3i9ddfp3v37rz88svst99+lW7fsmVLOnXqxMqVK+nZs2ed9ktd2/H73/+edevW0aVLFwoKCpg6deo2+7hXr17lF/JHqqoBxDL1pcFaq6bBWrMzfhw5sj1+HDk0WGvq4muw1vS2Yfny5d65c+fI4qdCuvdRVKIYrFVEREREUiypC/NFRKSeTj6ZtR98QIt01yENXrt27Vi4cGG6y5AkVNkJM7NqRxV09zmpL0dEpIH69a/5eNo0Dkh3HSKSMao7EnZn+LUJ0AOYR/Dw7nzgLYIBXEVERLbhlQwDINKQeR2euljlNWHuXuTuRcBKoLu793D3QqAbsKzOVYqIbI/69aPryJHpriIWTZo0Ye3atXX6oySSjdydtWvX1vjc04qSuSaso7svSEi00My61rI+ERHZTrRu3Zri4mJWr15d622TeYB3fUQdP44c2R4/jhxxtKGiJk2a0Lp161ptk0wn7F0zexB4HHDgPILBW0VERLbRqFEj2rdvX6dtp02bRrdu3VJcUXzx48iR7fHjyBFHG1IhmSEqLgIWAZcDI4HF4bwamdnxZrbUzJaZ2e+qWe8wMys1swHJxBURERHJdjUeCXP3b8xsDDDZ3ZcmG9jMcggec3QMUAzMNLNJ7r64kvX+BLxUq8pFREREsliNnTAz6w/cDjQG2ofXg/3B3fvXsGlPYJm7fxjGmQCcSnAkLdFlwETgsNqVnn0emP0A4xfW7uGgJSUltFjZIql1562aR0FeQR0qE5HInX02n7/3nsYJE5FyVtPdK2Y2m+Dh3dPcvVs4b76759ew3QDgeHcfHE6fDxzu7iMS1tkX+HsY/yHgeXd/upJYQ4GhAHl5eYUTJkxIvoUR2rhxI7m5uUmvf+W8K/lg4wcckJv8SEGlpaWVPjm+KkfteRQn731y0uvXtg21FXX8OHKoDemPH0cOtSH98ePIoTakP34cOeJoQ7KKiopmu3uPShdW9TyjshfwVvj1nYR585PY7izgwYTp84F7KqzzD6BX+H4cMKCmuNn87Mjt8dmOeqZgZuTI9vhx5Ii8DV995a+98EKkKbJ+H8WQQ21If/w4csTRhmRRzbMjk7k7cqGZ/RTIMbMOwC+BN5PYrhhokzDdGvikwjo9gAnhgH6tgBPNbLO7/28S8UVEsseJJ5JfUgLHH5/uSkQkQyRzd+RlQGfgW2A8sJ7gLsmazAQ6mFl7M2sMnANMSlzB3du7ezt3bwc8DfxCHTARERHZHiRzd+Qm4NrwlTR332xmIwjueswBHnb3RWY2LFw+pg71ioiIiDQIydwdeRDwa6Bd4vruflRN27r7ZGByhXmVdr7c/cKa4omIiIg0FMlcE/YPYAzwIFAabTkiIiIi24dkOmGb3f3+yCsREWnILryQz5Ys0ThhIlIumQvznzOzX5jZ3ma2W9kr8spERBqSCy/kM90ZKSIJkjkS9rPw628S5jmwf+rLERFpoNasodGXX6a7ChHJIMncHdk+jkJERBq0AQPoXFICp56a7kpEJENU2Qkzs6Pc/VUzO6Oy5e7+THRliYiIiDRs1R0JOxJ4FTilkmUOqBMmIiIiUkdVdsLcfVT49aL4yhERERHZPiRzYT5mdhLBo4ualM1z9z9EVZSIiIhIQ5fMiPljgKZAEcGArQOAtyOuS0SkYRk+nP8uWqRxwkSkXDLjhP3I3S8A1rn7jUBvoE20ZYmINDADB7L6qBqf9iYi25FkOmFfh183mdk+wPeAhq0QEamNjz9mp88/T3cVIpJBkumEPW9mLYDbgTnACmBCMsHN7HgzW2pmy8zsd5UsP9XM5pvZXDObZWZ9alG7iEj2OP98Drn11nRXISIZJJnBWm8K3040s+eBJu5e47DPZpYD/AU4BigGZprZJHdfnLDaK8Akd3czyweeAjrWthEiIiIi2aa6wVorHaQ1XJbMYK09gWXu/mG4zQTgVKC8E+buGxPW34Vg/DERERGRBq+6I2GVDdJaJpnBWvcFPk6YLgYOr7iSmZ0O/BHYEziphpgiIiIiDYK5R3PwyczOAo5z98Hh9PlAT3e/rIr1+wLXu/tPKlk2FBgKkJeXVzhhQlKXpEVu48aN5ObmJr3+lfOuBOCugrsiy1Fb2R4/jhxqQ/rjx5Ej6vhdR46ktLSUBffcE1mObN9HceRQG9IfP44ccbQhWUVFRbPdvUelC9292hewOzCa4KL82cCfgd2T2K438FLC9NXA1TVssxxoVd06hYWFnimmTp1aq/WLxhV50biiSHPUVrbHjyOH2pD++HHkiLwNkyb5/FtuiTRF1u+jGHKoDemPH0eOONqQLGCWV9GnSebuyAnAauBMgoFaVwNPJrHdTKCDmbU3s8bAOcCkxBXM7EAzs/B9d6AxsDaJ2CIi2eWUU1j7ox+luwoRySDJPLZoN//hDkmAm83stJo2cvfNZjYCeAnIAR5290VmNixcPoagY3eBmX1PMB7ZwLDXKCLSsCxdys4ffZTuKkQkgyTTCZtqZucQDB8BwdGw/0smuLtPBiZXmDcm4f2fgD8lV6qISBa75BIOLimBCy5IdyUikiGSOR15CfB34NvwNQG40sw2mNn6KIsTERERaaiSGax11zgKEREREdme1HgkzMx+XmE6x8xGRVeSiIiISMOXzOnIo81sspntbWaHAjMAHR0TERERqYdkTkf+1MwGAguATcAgd38j8spERBqS3/+elfPm0SLddYhIxqixE2ZmHYDLgYnAIcD5ZvaOu2+KurhMduVLVzJtyTRarGyR9DbzVs2jIK8guqJEJHP95Ces2zGZG9JFZHuRzCfCc8Cl7v5KOLDqlQQDsXaOtLI0Gzt7bLXL56+az3el3/HJhk9qjLXPrvsAUJBXwKAug1JSn4hkmblzyV22DPr1S3clIpIhkumE9XT39QDhQKp3mtmkGrZp8M7ufDa7rdqNL/K+qHHdoYVDY6hIRDLayJEcWFICgwenuxIRyRBVXphvZr8FcPf14cO4E10UaVUiIiIiDVx1d0eek/D+6grLjo+gFhEREZHtRnWdMKvifWXTIiIiIlIL1XXCvIr3lU2LiIiISC1Ud2F+QfhsSAN2TnhOpAFNIq9MRKQhufVWPpwzh+7prkNEMkaVR8LcPcfdm7n7ru6+Y/i+bLpRMsHN7HgzW2pmy8zsd5UsP9fM5oevN81Mg2iJSMP0ox+xvkuXdFchIhkkmccW1YmZ5QB/AU4AOgGDzKxThdWWA0e6ez5wE1D94FwiItnqzTdptnBhuqsQkQwSWScM6Aksc/cP3f07YAJwauIK7v6mu68LJ2cArSOsR0Qkfa65hv0ffDDdVYhIBomyE7Yv8HHCdHE4ryo/B16IsB4RERGRjGHBIPgRBA4GeD3O3QeH0+cTjL5/WSXrFgH3AX3cfW0ly4cCQwHy8vIKJ0yYEEnNidZsWlPjOjmbcyjdsbTG9Vo1bVXnOjZu3Ehubm6dt2/o8ePIoTakP34cOaKO33XkSEpLS1lwzz2R5cj2fRRHDrUh/fHjyBFHG5JVVFQ02917VLYsyqfJFgNtEqZbA9s8aNHM8oEHgRMq64ABuPtYwuvFevTo4f1iePZaTc+OBJJ+bNGAwgF1rmPatGlE2d5sjx9HDrUh/fHjyBF5G1q0oKSkJKvboO9zZuTI9vhx5IijDakQ5enImUAHM2tvZo0JRuDf6pmTZrYf8Axwvru/F2EtIiIiIhklsiNh7r7ZzEYALwE5wMPuvsjMhoXLxwDXA7sD95kZwOaqDtmJiGS1u+9m2axZ6ANORMpEeToSd58MTK4wb0zC+8HA4ChrEBHJCF27srGkJN1ViEgGifJ0pIiIlJkyhZazZ6e7ChHJIOqEiYjE4eabafvYY+muQkQyiDphIiIiImmgTpiIiIhIGqgTJiIiIpIGkd4dKVVLZjBYgN027VbjukMLh6aiJBEREYmROmEiInH4619Z+tZbHJ7uOkQkY6gTJiISh4MP5utPP013FSKSQXRNmIhIHJ57jt3ffDPdVYhIBlEnTEQkDnfeSZunnkp3FSKSQdQJExEREUkDdcJERERE0kCdMBEREZE0iPTuSDM7HvgzkAM86O63VVjeEXgE6A5c6+53RFnP9iSV45CBxiITERFJtcg6YWaWA/wFOAYoBmaa2SR3X5yw2hfAL4HToqpDRCQjPPYY706fTu901yEiGSPK05E9gWXu/qG7fwdMAE5NXMHdP3f3mcD3EdYhIpJ+bdrw7Z57prsKEckgUXbC9gU+TpguDueJiGx/nnySPV59Nd1ViEgGMXePJrDZWcBx7j44nD4f6Onul1Wy7g3AxqquCTOzocBQgLy8vMIJEyZEUnOiNZvW1LhOzuYcSncsrXG9Vk1b1Sl+sjmijl9VjmRs3LiR3NzcOm2bKTnUhvTHjyNH1PG7jhxJaWkpC+65J7Ic2b6P4sihNqQ/fhw54mhDsoqKima7e4/KlkV5YX4x0CZhujXwSV0CuftYYCxAjx49vF+/fvUuribJXKy+26rd+CLvixrXG1A4oE7xk80RdfyqciRj2rRpRP39ijqH2pD++HHkiLwNLVpQUlKS1W3Q9zkzcmR7/DhyxNGGVIjydORMoIOZtTezxsA5wKQI84mIiIhkjciOhLn7ZjMbAbxEMETFw+6+yMyGhcvHmNlewCygGbDFzEYCndx9fVR1SeokdbRQQ2CIiIhUKtJxwtx9MjC5wrwxCe8/IzhNKSIiIrJdibQTJiIioaefZtEbb/DjdNchIhlDjy0SEYlDq1Z837x5uqsQkQyiI2GSsfToJWlQxo1jryVLIAvu2BKReOhImIhIHMaNY68XX0x3FSKSQdQJExEREUkDnY6U7ZqG2RARkXRRJ0wkQqm8rk2dPBGRhkWnI0VERETSQEfCRLKY7iDNIpMnM//11+mb7jpEJGPoSJiISByaNmVLkybprkJEMoiOhIlItaK+eSHqo3kZc7TwvvvY5733NE6YiJRTJ0xEpJ6S6byd/NCfafn9jroBQ0TK6XSkiIiISBpEeiTMzI4H/gzkAA+6+20Vllu4/ERgE3Chu8+JsiYRkWwTxynVbD/tLJKNIuuEmVkO8BfgGKAYmGlmk9x9ccJqJwAdwtfhwP3hVxERkaSpkyfZKMojYT2BZe7+IYCZTQBOBRI7YacCf3N3B2aYWQsz29vdP42wLhERkVqLevBlHfGsOX5DY0H/J4LAZgOA4919cDh9PnC4u49IWOd54DZ3/084/QpwlbvPqhBrKFD23TgYWBpJ0bXXCliT5TmyPX4cOdSG9MePI4fakP74ceRQG9IfP44ccbQhWW3dfY/KFkR5JMwqmVexx5fMOrj7WCC57nWMzGyWu/fI5hzZHj+OHGpD+uPHkUNtSH/8OHKoDemPH0eOONqQClHeHVkMtEmYbg18Uod1RERERBqcKDthM4EOZtbezBoD5wCTKqwzCbjAAr2AL3U9mIiIiGwPIjsd6e6bzWwE8BLBEBUPu/siMxsWLh8DTCYYnmIZwRAVF0VVT0TiOEUadY5sjx9HDrUh/fHjyKE2pD9+HDnUhvTHjyNHxl3CVJnILswXERERkappxHwRERGRNFAnTERERCQN1AkTERERSQN1wkRERETSQJ2wFDCzjimO16iSea1SFHsHM9shfN/YzLqb2W6piF1Fvl9EFTuMnxu2oUWK4jUOHyxfNl1kZr8ysxNSET+MmZ+qWNXk2K9sn5hZOzMbYGZdUpyjh5mdbmanpPp3IIx/nJndb2aTzOyf4fvjU52nkrzXpzDWcWb2czNrV2H+xSmIbWZ2tpmdFb4/2sxGm9kvyn7HU83MXk1xvFYVps8L2zA08fewHvFPL/t8M7M9zOxvZrbAzJ40s9b1jR/GvcvMfpyKWFXE383MrjezweH3+Voze97MbjezlinKUWRm94a/ZxPN7DYzOzAVsRNyZP3vcxR0d2QKmNlH7r5fCuIUAY8BOwHvAEPdfUW4bI67d69n/NOAvwJbgGHANcBXwEHAcHd/rp7xr6w4C7gauBXA3e+qT/wwx33u/ovwfR/g78AHwIHAJe4+uZ7x5wH93H2dmf0GOJ1gKJUjgVnufnW9GhDkKAWWA+OB8RUeal9vZvY74BLgW+AO4NfAG0Av4KH6fh/M7EjgTqAEKAxjtwS+B85394/rEz/McTfBz+XfCAZ1hmAw5wuA99398vrmqCZ3qn6fbwX6AHOAU4C73f2ecFkqfp/vA/YEGgPrCT43niMY9mdVffeRmc2vOIvge7IUwN3r/c9E4n4ws98DRxD8Tp8MFLv7FfWMv9jdO4XvnwRmAP8AfgKc6+7H1Cd+GHc1sBLYA3iS4Hf6nfrGTYg/GVgANAMOCd8/BRwDFLj7qfWMfxuQB7wCnEbw2fQe8AvgVnf/R33ihznuJst/nyPj7nol8QJGV/G6B1ifohwzgc7h+wHA+0CvcPqdFMR/B9gLaE/woX1wOL8tQQejvvE3EHwIXQ+MCl/ryt6naB/NSXg/Fegevt8/RW1YmPB+FrBz+H5HYH6K2vAO0AW4hWCMvHnA74B2KYq/CNgZ2D38nuwRzt8lsX31rL8sZnvg2fD9McDLKWrDe1XMN4IP7frGX1/FawOwOUVtWADsGL5vQdCZ/5+yfZiK+OHXRsBaoLH/8LO6IAXxJwGPAx3Dz4h2wMfh+7Yp2kfvJLyfA+yS0KZUtGFpwvvZFZbNTWUbgA7AdeHv35Lwc++gFMSfG3414L+pbkPifg5/dt4I37dMxedFGCvrf5+jeul0ZPIuAhYCsyu8ZgHfpShHY3dfBODuTxP8V/KomZ1OJc/UrAt3/8zdlwMfuXvZf7QrSc2p6c4EA/PuAtzu7jcC69z9xvB9qjVz9zkA7v5hmLu+1iectlsDNAnf70jqTt+7uy9092vd/UBgCMERjX+b2ZspiF/q7l8THKn6muAPNO7+VQpiA+S4++rw/UcEf5Rx938B+6Yoxzdm1rOS+YcB36QgfgnQwd2bVXjtCqTqqR07uvtmAHcvITga1szM/kFw9Kq+ymJ/D8x09+/C6c1AaX2Du3t/YCLBoJcFHhyV/97dV4afGamws5l1M7NCgp+rr8Lc35OCNgDTzOwPZrZz+P40KD/r8GUK4kP42ezu77v7Te7eGTib4LOjXkfmQzuEpx3bALllp7bNbHdS83O0xX64JGUfws9Rd19H5c93rouG8PsciSgf4N3QzCT4r2CbP5JmdkOKcnxvZnu5+2cAHjxh4GjgeeCAVCQwsx3cfQtwccK8HFLwy+zuHwEDzOxU4F9m9j/1jVmJjuFpEgPamVlLD04d7kDw33N9DQOeCE9Lfg7MMrPXgHzgjymIDxU+2Nz9beBtM/sV0DcF8eeY2d8JOsOvEHTkXwSOBlJx6nOWmT0Uxj4VmAZgZk1JTUcY4ELgfjPblR9OX7Qh+O/2whTE/xtB53FVJcv+noL4AB+Y2ZHu/hqAu5cCPzezm4EzUxD/MzPLdfeN7l5+bY2Z7UWK/jF092fN7GXgJjMbTGr+6Cf6FCg7Pf6Fme3t7p+GHYzNKYg/AriW8BQqcIWZfUVw2vb8FMSHSjoq7j4fmE9wOUZ9/ZHgyBoEn9sPmpkDnYBU/HN7K/COmS0lOOo5HIJr6AiO0qfChWT/73MkdE1YksL/FL5x900R5vgJsNrd51WY3wK41N1vqWf8wwgOPX9TYX47oI+7P16f+BViNiX4gDjc3VPRsSiL27bCrE/d/bvwAt++7v5MCnLkAMcSXMOwI8GHxkvh0Yx6M7OfuntkHwxmtiNwFsF/6E8DPYGfEhy1ure+P8MW3DgyhOCPwDyCR5KVhkcb9kzhUZKyDsW+BH/oisv+QckG4f4gPCpZcdm+7v7fiPLuQnBa7/MUxy0AenvwyLlIhb+DO6Xy89bMmhMcnVybqphh3Fx335jKmJXkyCH4e705/P3uSnBqMiVHecK/b/sDy1L1OVdFnqz9fY6KOmH1ZGZtgHPc/fZszZHt8ePIoTakP35Cno7uvqTmNTMzfhw5sj1+HDnUhvjjm1mj8FRz4rxW7r4mG+JHQdeE1YGZtTKz4Wb2OsGpmLxsy5Ht8ePIoTakP34VXs7y+HHkyPb4ceRQG2KKb8EQGMXAJ2b2sm09ZEu9c0QdP0q6JixJ4bns0wlO6xwEPAvs7+4pGWsmjhzZHj+OHGpD+uOHOUZXtYjgTsOMjh9HjmyPH0cOtSH98UP/DzguvM55AME1w+e7+wxSc/F/1PEjo9ORSTKzr4G3gd8D/3F3N7MP3X3/bMmR7fHjyKE2pD9+mGMD8CuCsc4qutPd6zV4cdTx48iR7fHjyKE2pD9+mGOeuxckTHcGniEYmuc6r/+YeZHGj5RnwDgZ2fACrgDeIhim4hqCuxU/zKYc2R5fbdg+4oc5XgV+VMWy5ZkevyG0QfsoM3Jke/wwzixgrwrzWgNzgQ2ZHj/KV9oLyLYXwR0k1xIMxPgNcBUpGJAvzhzZHl9taPjxgd2Apqnc33HGbwht0D7KjBzZHj/M8ROCseYqzm8BXJvp8aN86XRkPZjZocAgYKC7p2Qcr7hzZHv8OHKoDemPn5An6+/wzPY2aB9lRo5sjx9Hjrju2K4P3R2ZJDM70Co8pNXdFwAvAsdlQ45sjx9HDrUh/fEryZf1d3hmexu0jzIjR7bHjyNHHG1IJXXCknc3wXOoKtoEpGpk+KhzZHv8OHJEHT+OHNkeHzPb1cwusGCk/7cJHtC+v7sf4O6/zvT4ceTI9vhx5FAb0h8/jhxxtCEy6T4fmi0vqnmQKSl40GwcObI9vtqwfcQP43wNvAYcwQ93cafy5oJI4zeENmgfZUaObI/fUNoQ1UtHwpLXpJplO2dJjmyPH0cOtSH98SG467IJcD9wtZml+jqzqOPHkSPb48eRQ21If/w4csTRhkioE5a8mWY2pOJMM/s5MDtLcmR7/DhyqA3pj4+7/4+7Hw70Jxhs8X+BfczsKjM7KNPjx5Ej2+PHkUNtSH/8OHLE0Yao6O7IJJlZHsHI4N/xwx+aHkBj4HRPwYNIo86R7fHjyKE2pD9+NXmz/g7PbG+D9lFm5Mj2+HHkiKMNqaBOWC2ZWRHQJZxc5O6vZluObI8fRw61Ib3xzexAIM/d36gwvy/wibsvy+T4ceTI9vhx5FAb0h8/jhxxtCEqOh2ZJDNrYmYjgTMJjgDcH8EftUhzZHv8OHKoDemPH7qbLL/DM4Yc2R4/jhxRx48jR7bHjyNH1PEjo05Y8h4lOOWyADgBuCMLc2R7/DhyqA3pjw/Qzt3nV5zp7rOAdlkQP44c2R4/jhxqQ/rjx5EjjjZEYsd0F5BFOrn7oQBm9hDBWCTZliPb48eRQ21If3xoGHd4ZnsbtI8yI0e2x48jRxxtiISOhCXv+7I37r45S3Nke/w4cqgN6Y8PDeAOzxhyZHv8OHKoDemPH0eOONoQCV2YnyQzKwW+Kpsk6F1vCt+7uzfL9BzZHj+OHGpD+uOHObL+Ds9sb4P2UWbkyPb4ceRI1x3bqaBOmIhkrGy+wzOuHNkeP44cakP648eRI442pJquCRORjGNmTYBhBM+AWwA8lMpTn1HHjyNHtsePI4fakP74ceSIow1R0ZEwEck4ZvYkwbVn/ya4A3OFu4/Mlvhx5Mj2+HHkUBvSHz+OHHG0ISrqhIlIxjGzBQl3YO4IvO3u3bMlfhw5sj1+HDnUhvTHjyNHHG2Iiu6OFJFM1BDu8Mz2NmgfZUaObI8fR4442hAJHQkTkYzTQO7wzOo2aB9lRo5sjx9HjjjaEBV1wkRERETSQKcjRURERNJAnTARERGRNFAnTEQws1Izm2tmC83sH2bWtIr13qxj/B5mNroe9W2s67bZxMxGVrPvHzSzTrWMt13sN5FspWvCRAQz2+juueH7J4DZ7n5XwvIcdy/NhPoaMjNbAfRw9zUpirdd7DeRbKUjYSJS0b+BA82sn5lNNbO/E4xCXX5kJVw2zcyeNrMlZvaEmVm47DAze9PM5pnZ22a2a7j+8+HyG8zsMTN71czet/DBu2aWa2avmNkcM1tgZqfWVKiZXWBm88Ncj4Xz2oZx5odf9wvnjzOz+8M2fWhmR5rZw2b2rpmNS4i50czuDOt4xcz2COd3NbMZYdxnzaxlOH+amf0pbOt7ZnZEOD/HzG43s5nhNpdUt+/M7JfAPsBUM5taSVunmVmPhBpvCds9w4Jn52Fm7c1sepjzpgrb/yahlhvDeaeb2ZQw/95h/Xsl9VMiIvWmTpiIlLNgoMMTCDtdQE/gWnev7DRYN2Ak0AnYH/ixmTUGngQud/cC4CfA15Vsmw+cBPQGrjezfYBvCB622x0oAu4s69hVUWtn4FrgqDDX5eGie4G/uXs+8ASQeBq0JXAUcAXwHPA/QGfgUDPrGq6zCzAnrOM1YFQ4/2/AVWHcBQnzAXZ0957h/iib/3PgS3c/DDgMGGJm7avad+4+GvgEKHL3oqranVDjjLDdrwNDwvl/Bu4Pc5Y/tNjMjgU6EHw/uwKFZtbX3Z8N17sUeAAYlckPOxZpaNQJExGAnc1sLjAL+Ah4KJz/trsvr2Kbt9292N23AHOBdsDBwKfuPhPA3ddXMXjiP9396/C021SCzoEBt5rZfGAKsC+QV03NRwFPl526c/cvwvm9gb+H7x8D+iRs85wH12AsAFa5+4Kw/kVh/QBbCDqSAI8DfcysOdDC3V8L5z8K9E2I+0z4dXZCnGOBC8L9+hawO0FHCCrfd7XxHfB8JTl/DIwP3z+WsP6x4esdYA7QMaGWy4CrgW/dfTwiEhs9wFtEAL52966JM8KDUF9Vunbg24T3pQSfJwYkc6FpxXUcOBfYAyh09+/D66OaVBOjLrnKat7C1vVvoerPw2RylMUq2w9l9V3m7i8lrmhm/ah839XG9/7DBb0Vt6+sXgP+6O5/rWTZvgTtzzOzHcKOoYjEQEfCRCSVlgD7mNlhAOH1YJV1ME41syZmtjvQD5gJNAc+DztgRUDbGnK9ApwdxsDMdgvnvwmcE74/F/hPLduwAzAgfP9T4D/u/iWwrux6L+B8glOV1XkJGG5mjcL6DjKzXWrYZgOway3rTfQGW7c9sZaLzazs5ot9zWzP8HvzCEE73wWurEduEaklHQkTkZRx9+/MbCBwj5ntTHA92E8qWfVt4P+A/YCb3P0TC+7KfM7MZhGcoltSQ65FZnYL8JoFjy15B7gQ+CXwsJn9BlgNXFTLZnwFdDaz2cCXwMBw/s+AMRYMIfFhEnEfJDhNOCe8tm01cFoN24wFXjCzT5O4LqwylwN/N7PLgYllM939ZTM7BJgeHuHcCJwHDAP+7e7/Dk+bzjSz/3P3d+uQW0RqSUNUiEiszOwGYKO735HuWipjGtZBRGKi05EiIiIiaaAjYSIiIiJpoCNhIiIiImmgTpiIiIhIGqgTJiIiIpIG6oSJiIiIpIE6YSIiIiJpoE6YiIiISBr8f9G+zdUfhyIqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Determine explained variance using explained_variance_ratio_ attribute\n",
    "exp_var_pca = pca.explained_variance_ratio_\n",
    "\n",
    "# Calculate cumulative sum of eigenvalues; This will be used to create step plot\n",
    "# for visualizing the variance explained by each principal component.\n",
    "cum_sum_eigenvalues = np.cumsum(exp_var_pca)\n",
    "display(cum_sum_eigenvalues)\n",
    "\n",
    "# Create the visualization plot\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "\n",
    "# Plot individual explained variance\n",
    "plt.bar([f'PCA{i}' for i in np.arange(1, len(eigenvalue_df) + 1)],\n",
    "        exp_var_pca,\n",
    "        align='center',\n",
    "        color='green',\n",
    "        alpha=0.4,\n",
    "        label='Individual explained variance'\n",
    "        )\n",
    "\n",
    "# Plot cumulative explained variance as a step plot\n",
    "plt.step(range(0, len(cum_sum_eigenvalues)),\n",
    "         cum_sum_eigenvalues,\n",
    "         where='mid',\n",
    "         color='green',\n",
    "         alpha=0.99,\n",
    "         label='Cumulative explained variance'\n",
    "         )\n",
    "\n",
    "plt.xlabel('Principal component index')\n",
    "plt.ylabel('Explained variance ratio')\n",
    "\n",
    "# Set x-axis labels\n",
    "plt.xticks([f'PCA{i}' for i in np.arange(1, len(eigenvalue_df) + 1)], rotation=90)\n",
    "\n",
    "# Set y-axis limits and tick positions\n",
    "plt.yticks(np.arange(0, 1.1, step=0.1))\n",
    "plt.ylim([0, 1])\n",
    "\n",
    "# Add gridlines\n",
    "plt.grid()\n",
    "\n",
    "# Add lines for PCA12 and 0.915\n",
    "plt.axvline(x='PCA12', color='r', linestyle='--')  # Vertical line\n",
    "plt.axhline(y=0.915, color='r', linestyle='--')  # Horizontal line\n",
    "\n",
    "# Add legend\n",
    "plt.legend(loc='center right')\n",
    "\n",
    "# Add a label to the plot\n",
    "plt.text('PCA12', 0.85, \"0.915\")\n",
    "plt.scatter('PCA12', 0.915, color='#ff0000', zorder=1)\n",
    "\n",
    "# Save the plot (optional, commented out)\n",
    "plt.savefig('/home/compass/out/Explained_variance.png',\n",
    "            format='png', dpi=100, bbox_inches='tight', pad_inches=1)\n",
    "# display\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "175ff9ab",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.15"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "294.901px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}