pca_logistic_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By Shirley Yu and Angela Zhang"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "import statsmodels.formula.api as smf\n",
    "from sklearn.cluster import KMeans\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_hot = pd.read_csv(\"one_hot_responses.csv\")\n",
    "one_hot = one_hot.drop(columns=[\"Unnamed: 0\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA to Reduce Dimensionality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "# top n pca\n",
    "n = 50\n",
    "predictors = [one_hot.columns[i] for i in range(len(one_hot.columns))]\n",
    "X = one_hot[predictors].values.T\n",
    "\n",
    "#don't need to standardize since they're dummy variables? idk I asked on piazza\n",
    "#but when standardizing, if we try to have more than 30 PCs it causes a LinAlgError: Singular matrix error\n",
    "#n_i, n_j = X.shape\n",
    "#stds = np.std(X, axis=1)\n",
    "#stds = np.repeat(np.reshape(stds, (-1,1)), n_j, axis=1)\n",
    "#S = X/stds\n",
    "\n",
    "pca = PCA(n_components=n)\n",
    "pcs = pca.fit_transform(X.T)\n",
    "df_pcs = pd.DataFrame(pcs, columns=['PC'+str(i+1) for i in range(n)])\n",
    "\n",
    "# predict whether the person was gay based on their survey responses\n",
    "df_pcs['y'] = one_hot['q0001'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression for Predicting Masculinity Using Various Numbers of Top PCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accuracy(predictions, true):\n",
    "    correct = (predictions == true)\n",
    "    acc = sum(correct) / len(correct)\n",
    "    return acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "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>formula</th>\n",
       "      <th>mean_test_acc</th>\n",
       "      <th>mean_test_mse</th>\n",
       "      <th>mean_train_acc</th>\n",
       "      <th>mean_train_mse</th>\n",
       "      <th>num_pc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>y ~ 1 + PC1</td>\n",
       "      <td>0.971532</td>\n",
       "      <td>0.028468</td>\n",
       "      <td>0.971517</td>\n",
       "      <td>0.028483</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>y ~ 1 + PC1 + PC2</td>\n",
       "      <td>0.971517</td>\n",
       "      <td>0.028483</td>\n",
       "      <td>0.971517</td>\n",
       "      <td>0.028483</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3</td>\n",
       "      <td>0.971540</td>\n",
       "      <td>0.028460</td>\n",
       "      <td>0.971517</td>\n",
       "      <td>0.028483</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4</td>\n",
       "      <td>0.971509</td>\n",
       "      <td>0.028491</td>\n",
       "      <td>0.971517</td>\n",
       "      <td>0.028483</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5</td>\n",
       "      <td>0.971513</td>\n",
       "      <td>0.028487</td>\n",
       "      <td>0.971517</td>\n",
       "      <td>0.028483</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6</td>\n",
       "      <td>0.971501</td>\n",
       "      <td>0.028499</td>\n",
       "      <td>0.971517</td>\n",
       "      <td>0.028483</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7</td>\n",
       "      <td>0.970915</td>\n",
       "      <td>0.029085</td>\n",
       "      <td>0.971655</td>\n",
       "      <td>0.028345</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.970903</td>\n",
       "      <td>0.029097</td>\n",
       "      <td>0.971586</td>\n",
       "      <td>0.028414</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.970873</td>\n",
       "      <td>0.029127</td>\n",
       "      <td>0.971517</td>\n",
       "      <td>0.028483</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.972126</td>\n",
       "      <td>0.027874</td>\n",
       "      <td>0.972824</td>\n",
       "      <td>0.027176</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.970286</td>\n",
       "      <td>0.029714</td>\n",
       "      <td>0.971999</td>\n",
       "      <td>0.028001</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.970899</td>\n",
       "      <td>0.029101</td>\n",
       "      <td>0.971861</td>\n",
       "      <td>0.028139</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.970911</td>\n",
       "      <td>0.029089</td>\n",
       "      <td>0.972549</td>\n",
       "      <td>0.027451</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.969665</td>\n",
       "      <td>0.030335</td>\n",
       "      <td>0.972067</td>\n",
       "      <td>0.027933</td>\n",
       "      <td>14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.968407</td>\n",
       "      <td>0.031593</td>\n",
       "      <td>0.972205</td>\n",
       "      <td>0.027795</td>\n",
       "      <td>15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.969044</td>\n",
       "      <td>0.030956</td>\n",
       "      <td>0.972343</td>\n",
       "      <td>0.027657</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.970888</td>\n",
       "      <td>0.029112</td>\n",
       "      <td>0.971998</td>\n",
       "      <td>0.028002</td>\n",
       "      <td>17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.969669</td>\n",
       "      <td>0.030331</td>\n",
       "      <td>0.971724</td>\n",
       "      <td>0.028276</td>\n",
       "      <td>18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.969665</td>\n",
       "      <td>0.030335</td>\n",
       "      <td>0.971655</td>\n",
       "      <td>0.028345</td>\n",
       "      <td>19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.967196</td>\n",
       "      <td>0.032804</td>\n",
       "      <td>0.971311</td>\n",
       "      <td>0.028689</td>\n",
       "      <td>20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.969032</td>\n",
       "      <td>0.030968</td>\n",
       "      <td>0.971586</td>\n",
       "      <td>0.028414</td>\n",
       "      <td>21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.966536</td>\n",
       "      <td>0.033464</td>\n",
       "      <td>0.972205</td>\n",
       "      <td>0.027795</td>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.965332</td>\n",
       "      <td>0.034668</td>\n",
       "      <td>0.971586</td>\n",
       "      <td>0.028414</td>\n",
       "      <td>23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.968400</td>\n",
       "      <td>0.031600</td>\n",
       "      <td>0.971311</td>\n",
       "      <td>0.028689</td>\n",
       "      <td>24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.968438</td>\n",
       "      <td>0.031562</td>\n",
       "      <td>0.970898</td>\n",
       "      <td>0.029102</td>\n",
       "      <td>25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.967184</td>\n",
       "      <td>0.032816</td>\n",
       "      <td>0.970485</td>\n",
       "      <td>0.029515</td>\n",
       "      <td>26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.967180</td>\n",
       "      <td>0.032820</td>\n",
       "      <td>0.970898</td>\n",
       "      <td>0.029102</td>\n",
       "      <td>27</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.969051</td>\n",
       "      <td>0.030949</td>\n",
       "      <td>0.971586</td>\n",
       "      <td>0.028414</td>\n",
       "      <td>28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.967157</td>\n",
       "      <td>0.032843</td>\n",
       "      <td>0.971448</td>\n",
       "      <td>0.028552</td>\n",
       "      <td>29</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.968415</td>\n",
       "      <td>0.031585</td>\n",
       "      <td>0.971930</td>\n",
       "      <td>0.028070</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.967786</td>\n",
       "      <td>0.032214</td>\n",
       "      <td>0.972274</td>\n",
       "      <td>0.027726</td>\n",
       "      <td>31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.968423</td>\n",
       "      <td>0.031577</td>\n",
       "      <td>0.972411</td>\n",
       "      <td>0.027589</td>\n",
       "      <td>32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.967157</td>\n",
       "      <td>0.032843</td>\n",
       "      <td>0.972618</td>\n",
       "      <td>0.027382</td>\n",
       "      <td>33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.966559</td>\n",
       "      <td>0.033441</td>\n",
       "      <td>0.972480</td>\n",
       "      <td>0.027520</td>\n",
       "      <td>34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.969661</td>\n",
       "      <td>0.030339</td>\n",
       "      <td>0.972480</td>\n",
       "      <td>0.027520</td>\n",
       "      <td>35</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.969017</td>\n",
       "      <td>0.030983</td>\n",
       "      <td>0.973030</td>\n",
       "      <td>0.026970</td>\n",
       "      <td>36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.965942</td>\n",
       "      <td>0.034058</td>\n",
       "      <td>0.973099</td>\n",
       "      <td>0.026901</td>\n",
       "      <td>37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.967794</td>\n",
       "      <td>0.032206</td>\n",
       "      <td>0.973237</td>\n",
       "      <td>0.026763</td>\n",
       "      <td>38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.965965</td>\n",
       "      <td>0.034035</td>\n",
       "      <td>0.973788</td>\n",
       "      <td>0.026212</td>\n",
       "      <td>39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.966582</td>\n",
       "      <td>0.033418</td>\n",
       "      <td>0.973581</td>\n",
       "      <td>0.026419</td>\n",
       "      <td>40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.967196</td>\n",
       "      <td>0.032804</td>\n",
       "      <td>0.974200</td>\n",
       "      <td>0.025800</td>\n",
       "      <td>41</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.966544</td>\n",
       "      <td>0.033456</td>\n",
       "      <td>0.973374</td>\n",
       "      <td>0.026626</td>\n",
       "      <td>42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.963465</td>\n",
       "      <td>0.036535</td>\n",
       "      <td>0.973856</td>\n",
       "      <td>0.026144</td>\n",
       "      <td>43</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.964707</td>\n",
       "      <td>0.035293</td>\n",
       "      <td>0.973375</td>\n",
       "      <td>0.026625</td>\n",
       "      <td>44</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.965313</td>\n",
       "      <td>0.034687</td>\n",
       "      <td>0.973994</td>\n",
       "      <td>0.026006</td>\n",
       "      <td>45</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.964704</td>\n",
       "      <td>0.035296</td>\n",
       "      <td>0.973650</td>\n",
       "      <td>0.026350</td>\n",
       "      <td>46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.965317</td>\n",
       "      <td>0.034683</td>\n",
       "      <td>0.973443</td>\n",
       "      <td>0.026557</td>\n",
       "      <td>47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.965317</td>\n",
       "      <td>0.034683</td>\n",
       "      <td>0.973512</td>\n",
       "      <td>0.026488</td>\n",
       "      <td>48</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>48</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.964083</td>\n",
       "      <td>0.035917</td>\n",
       "      <td>0.973719</td>\n",
       "      <td>0.026281</td>\n",
       "      <td>49</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.965332</td>\n",
       "      <td>0.034668</td>\n",
       "      <td>0.973168</td>\n",
       "      <td>0.026832</td>\n",
       "      <td>50</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              formula  mean_test_acc  \\\n",
       "0                                         y ~ 1 + PC1       0.971532   \n",
       "1                                   y ~ 1 + PC1 + PC2       0.971517   \n",
       "2                             y ~ 1 + PC1 + PC2 + PC3       0.971540   \n",
       "3                       y ~ 1 + PC1 + PC2 + PC3 + PC4       0.971509   \n",
       "4                 y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5       0.971513   \n",
       "5           y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6       0.971501   \n",
       "6     y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7       0.970915   \n",
       "7   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.970903   \n",
       "8   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.970873   \n",
       "9   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.972126   \n",
       "10  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.970286   \n",
       "11  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.970899   \n",
       "12  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.970911   \n",
       "13  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.969665   \n",
       "14  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.968407   \n",
       "15  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.969044   \n",
       "16  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.970888   \n",
       "17  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.969669   \n",
       "18  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.969665   \n",
       "19  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.967196   \n",
       "20  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.969032   \n",
       "21  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.966536   \n",
       "22  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.965332   \n",
       "23  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.968400   \n",
       "24  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.968438   \n",
       "25  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.967184   \n",
       "26  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.967180   \n",
       "27  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.969051   \n",
       "28  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.967157   \n",
       "29  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.968415   \n",
       "30  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.967786   \n",
       "31  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.968423   \n",
       "32  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.967157   \n",
       "33  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.966559   \n",
       "34  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.969661   \n",
       "35  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.969017   \n",
       "36  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.965942   \n",
       "37  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.967794   \n",
       "38  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.965965   \n",
       "39  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.966582   \n",
       "40  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.967196   \n",
       "41  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.966544   \n",
       "42  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.963465   \n",
       "43  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.964707   \n",
       "44  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.965313   \n",
       "45  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.964704   \n",
       "46  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.965317   \n",
       "47  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.965317   \n",
       "48  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.964083   \n",
       "49  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.965332   \n",
       "\n",
       "    mean_test_mse  mean_train_acc  mean_train_mse  num_pc  \n",
       "0        0.028468        0.971517        0.028483       1  \n",
       "1        0.028483        0.971517        0.028483       2  \n",
       "2        0.028460        0.971517        0.028483       3  \n",
       "3        0.028491        0.971517        0.028483       4  \n",
       "4        0.028487        0.971517        0.028483       5  \n",
       "5        0.028499        0.971517        0.028483       6  \n",
       "6        0.029085        0.971655        0.028345       7  \n",
       "7        0.029097        0.971586        0.028414       8  \n",
       "8        0.029127        0.971517        0.028483       9  \n",
       "9        0.027874        0.972824        0.027176      10  \n",
       "10       0.029714        0.971999        0.028001      11  \n",
       "11       0.029101        0.971861        0.028139      12  \n",
       "12       0.029089        0.972549        0.027451      13  \n",
       "13       0.030335        0.972067        0.027933      14  \n",
       "14       0.031593        0.972205        0.027795      15  \n",
       "15       0.030956        0.972343        0.027657      16  \n",
       "16       0.029112        0.971998        0.028002      17  \n",
       "17       0.030331        0.971724        0.028276      18  \n",
       "18       0.030335        0.971655        0.028345      19  \n",
       "19       0.032804        0.971311        0.028689      20  \n",
       "20       0.030968        0.971586        0.028414      21  \n",
       "21       0.033464        0.972205        0.027795      22  \n",
       "22       0.034668        0.971586        0.028414      23  \n",
       "23       0.031600        0.971311        0.028689      24  \n",
       "24       0.031562        0.970898        0.029102      25  \n",
       "25       0.032816        0.970485        0.029515      26  \n",
       "26       0.032820        0.970898        0.029102      27  \n",
       "27       0.030949        0.971586        0.028414      28  \n",
       "28       0.032843        0.971448        0.028552      29  \n",
       "29       0.031585        0.971930        0.028070      30  \n",
       "30       0.032214        0.972274        0.027726      31  \n",
       "31       0.031577        0.972411        0.027589      32  \n",
       "32       0.032843        0.972618        0.027382      33  \n",
       "33       0.033441        0.972480        0.027520      34  \n",
       "34       0.030339        0.972480        0.027520      35  \n",
       "35       0.030983        0.973030        0.026970      36  \n",
       "36       0.034058        0.973099        0.026901      37  \n",
       "37       0.032206        0.973237        0.026763      38  \n",
       "38       0.034035        0.973788        0.026212      39  \n",
       "39       0.033418        0.973581        0.026419      40  \n",
       "40       0.032804        0.974200        0.025800      41  \n",
       "41       0.033456        0.973374        0.026626      42  \n",
       "42       0.036535        0.973856        0.026144      43  \n",
       "43       0.035293        0.973375        0.026625      44  \n",
       "44       0.034687        0.973994        0.026006      45  \n",
       "45       0.035296        0.973650        0.026350      46  \n",
       "46       0.034683        0.973443        0.026557      47  \n",
       "47       0.034683        0.973512        0.026488      48  \n",
       "48       0.035917        0.973719        0.026281      49  \n",
       "49       0.034668        0.973168        0.026832      50  "
      ]
     },
     "execution_count": 126,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 10 fold cross validation, logistic regression\n",
    "df_pcs = df_pcs.fillna(0)\n",
    "mean_mses = []\n",
    "kf = KFold(n_splits=10, shuffle=True)\n",
    "pc_numbers = [i+1 for i in range(n)]\n",
    "# linear models using top k=1,2,3,...,20\n",
    "for num_pc in pc_numbers:\n",
    "    # do 10-fold cross validation\n",
    "    train_mses = []\n",
    "    test_mses = []\n",
    "    train_accs = []\n",
    "    test_accs = []\n",
    "    for train_idx, test_idx in kf.split(df_pcs):\n",
    "        df_train = df_pcs.iloc[train_idx, :]\n",
    "        df_test = df_pcs.iloc[test_idx, :]\n",
    "        \n",
    "        # generate formula\n",
    "        formula = 'y ~ 1'\n",
    "        for num in range(1, num_pc+1):\n",
    "            formula += (' + PC%d' % num)\n",
    "            \n",
    "        # fit model with df_train\n",
    "        res = smf.logit(formula=formula, data=df_train).fit(disp=0)\n",
    "\n",
    "        # MSE of training set\n",
    "        train_pred = 1*(res.predict(df_train) > 0.5)\n",
    "        mse_train = np.mean((train_pred - df_train.y)**2)\n",
    "\n",
    "        # MSE of test set\n",
    "        test_pred = 1*(res.predict(df_test) > 0.5)\n",
    "        mse_test = np.mean((test_pred - df_test.y)**2)\n",
    "        \n",
    "        acc_train = accuracy(train_pred, df_train.y)\n",
    "        acc_test = accuracy(test_pred, df_test.y)\n",
    "\n",
    "        # store results\n",
    "        train_mses.append(mse_train)\n",
    "        test_mses.append(mse_test)\n",
    "        train_accs.append(acc_train)\n",
    "        test_accs.append(acc_test)\n",
    "\n",
    "    # report average training MSE and testing MSE\n",
    "    mean_train_mse = np.mean(train_mses)\n",
    "    mean_test_mse = np.mean(test_mses)\n",
    "    mean_train_acc = np.mean(train_accs)\n",
    "    mean_test_acc = np.mean(test_accs)\n",
    "    \n",
    "    mean_mses.append({'num_pc': num_pc,\n",
    "                      'formula': formula,\n",
    "                    'mean_train_mse': mean_train_mse,\n",
    "                     'mean_test_mse': mean_test_mse,\n",
    "                     'mean_train_acc': mean_train_acc,\n",
    "                     'mean_test_acc': mean_test_acc})\n",
    "    \n",
    "    \n",
    "df_res = pd.DataFrame(mean_mses)\n",
    "df_res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXd8VeX5wL9PBoQACRAgbMJSmbJFBWWJuLVu6/rVukdtq1XrqlqrVVutVetu1bqtg7YOynSAbJCwhIQVAiEhIQkJ2c/vj/fc5HJzb+65SW7m+/187ufee8573vOcXDjPeZ8pqorFYrFYLLUlorEFsFgsFkvzxioSi8VisdQJq0gsFovFUiesIrFYLBZLnbCKxGKxWCx1wioSi8VisdQJq0gsDYaIJIrI1yKSLyJ/amx5QkVEkkRERSSqsWXxRkTaici/RSRXRD5sbHmaGyJyg4jMa2w5mjNWkbRSRGSHiJSISFef7Wudm2WS872PiPxLRLKcG9V6Ebna2ee5sR7yeV0c4LTXAVlAnKr+uh6u4Wrn/Hf6bE8Tkal1nb8ZcQGQCCSo6oXeO0TkRa/fpURESr2+f1GfQjg35DJn7jwRWS0is732dxaRv4rIbmfMNhF5SkS6+JnrGBEp87P9PRG5rz7lttQdq0haN9uBSz1fRGQk0M5nzFvAbqA/kABcCWT4jOmkqh28Xu8HOF9/YKPWIgu2hlVANnCXiMSFOmdjUs+rmv7Aj6pa7carqjd4fhfgD8D7Xr/TafUog4dFzrk6A+8CH4pIBxFpBywCBgEzgThgMlAIjA2DHJYGxCqS1s1bGMXg4SrgTZ8xE4B/qGqBqpap6hpVDflJVkT+4cz/G+dpdKaItBWRZ0Qk3Xk9IyJtnfFTnZXFXSKyD/h7gKk3AUuBXwY6r4j83uv7VBFJ8/q+Q0TuFJEfRKRARF5zTHBfOCa4eSLS2Wfanzny7hWRX3vNFSEid4tIiogcEJEPPE/bXqu3a0RkF7BARGJE5J/O2IMiskJEEgNcx1ARWeSM2yAiZzvbHwIeAC52/q7X1PAz+EVEzheRjc7c80RkiNe+fSLyGxHZLCLZIvKy5zeqCVUtB14HOgBJwDWYB5HzVXWLqlao6j5VfUBVa2VWclYt3zor5UwRedNr3wgRWSAiOSKySUTO9drXXUQ+d1ZNSzGK2FIHrCJp3XwPxDk3qUjgYuCffsY8LyKXiEi/2p5IVa8G3gaecJ6G5wH3ApOA0cCxwETA22zRA+iC+Y9+XQ3T3w/80p+JxCXnA6cARwFnAV8AvwW6Yv6P3OYzfhowBJgF3C0iM53ttwHnAicDvYAc4HmfY08GhgKnYhRrPNAXc5O9ATjsK5yIRAP/BuYC3YFbgbdF5GhVfZAjVxqvhXLhIjIC+AdwkzP3YmCOz4rpUmA6cDQwBriTIDjHXwPkYla+M4H/qmq166sDjwGfAp2AfsBLzrnjgP8Br2F+wyuB10VksHPcy5iVbCJwI/CzepSpVWIVicWzKjkF2Azs8dl/IfAN5ma9XYwPZYLPmCznadbzGury3D8FHlbV/aqaCTwEXOG1vwJ4UFWLa7oBqepazE32Lpfn9eWvqpqhqnsw17rMWXkVA59gbp7ePOSs0NZjVkoe8+D1wL2qmuYc+zvgAp+b8u+cYw8DpRgFMlhVy1V1larm+ZFvEubJ/nFVLVHVBcB/vM5bFy4FPlHVRapaglFKXYHxXmP+oqrpzm/0WJDzniwiB4F9wDnAuapagLnOvfUgrzelmNVOD1U9rKrfOdvPA5JV9W3n77oCo4jPF5EY4GzgPueYtZgHHEsdsIrE8hZwGXA11c1aqGqOqt6tqsMxT3BrgU9FRLyGdVXVTl6vTS7P3QvY6fV9p7PNQ6aqFrmc6wHgRhHp4XK8N94+n8N+vnfwGb/b67O3zP2BTzwKFWN2K8f83fwd+xbwFfCeYyp7wll9+NIL2K2qFT7n7V3zZbniiN/AMUnt8Zk70PX6Y7Hzb6Crqp6oqouc7QeAniHIVQZEiIjvPSoao0DAmDNjgTWOafJyZ3t/4CTvhxvMqrMnZpUrfq7JUgesImnlqOpOjOnhdODjIGOzgKcwN5LampG8SedI+3Q/Z1vlKd1OpKqbMfL/1mdXAeZm46E2isaXvl6fvWXeDZzmo1RjnJVOpaheMpeq6kOqOgw4ATiTI31WHtKBvj431X5UXz3WhiN+A8fE2dtn7kDXGwrzgNOdFYEbPDd6X//FAJwbv6ruUdWfYRTEbRjzVT/n2Lk+v0MHVb0ds1JSP9dkqQNWkVjA2LKnOyaIIxCRPzqOyygR6YixKW9T1QP1cN53gftEpJuYMOQHqO6jCYWHgP/D2Mw9rMXcwLo4q5Xb6zC/h/tFJFZEhjvn80SpvQg8KiL9AZzrOifQJCIyTURGOjfvPMyTdrmfocswCvE3IhItJrT5LOC9eriW94HzROQkZzV0N2b1sNJrzG0i0tP5je6m6npD4TWMX+JDETlKDN1E5EERmeE72DENfgY8JiZsOFpM2Hl/jP8DEblYRHo5UYAHnUPLMH6TMc7+aBFpIyKTROQoZ4X7b+AhMfk3ozAmVksdsIrEgqqmqOrKALtjMX6Cg0Aq5j/y2T5jDsqReSS/cnnq32NuWD8A64HVzrZaoarbMeai9l6b3wLWATswfpTa3AR9WQxsA+YDT6nqXGf7X4A5wFwRyccEKhxXwzw9gI8wSmSTM281Rer4Ls4GTsPk4bwAXOmswuqEqv6AeZB4CcgEZgDn+IQSvwcsBLZifqcnanGew8BUzGpiAZCPibZrj/nd/XEtUARswJgbf4ZZ8XkeYo4HVonIIeBD4DrHl5ODCWb4P4xfJh3z78pjNrweY27McK47UESgxSViG1tZLJZAiAm9vkBVv21sWSxNF7sisVgsFkudsIrEYrFYLHXCmrYsFovFUifsisRisVgsdaJJlcMOF127dtWkpKTGFsNisViaFatWrcpS1W7BxrUKRZKUlMTKlYGiWy0Wi8XiDxFxlfVvTVsWi8ViqRNWkVgsFoulTlhFYrFYLJY60Sp8JP4oLS0lLS2NoiK3xWUtdSEmJoY+ffoQHe2vuK3FYmnOtFpFkpaWRseOHUlKSuLIiuiW+kZVOXDgAGlpaQwYMKCxxbFYLPVMqzVtFRUVkZCQYJVIAyAiJCQk2NWfxdJCabWKBLBKpAGxf2uLpeXSqhWJxWKxhMSWLyBnR2NL0eSwiqSROHjwIC+88EKtjj399NM5ePBgjWMeeOAB5s2bV6v568Knn37Kxo0bG/y8FkvYqaiAD66E+Y80tiRNDqtIGomaFEl5ub8meVV8/vnndOrUqcYxDz/8MDNnzqy1fLXFKhJLi6UgE8pLIHWRUSqWSqwiaSTuvvtuUlJSGD16NHfeeSeLFi1i2rRpXHbZZYwcORKAc889l3HjxjF8+HBefvnlymOTkpLIyspix44dDB06lGuvvZbhw4cza9YsDh8+DMDVV1/NRx99VDn+wQcfZOzYsYwcOZLNm01jvczMTE455RTGjh3L9ddfT//+/cnKyjpCzvLycq6++mpGjBjByJEjefrppwFISUlh9uzZjBs3jilTprB582aWLFnCnDlzuPPOOxk9ejQpKSlh/ztaLA1GntPGvjALMpIbV5YmRqsN//XmoX9vYGN6Xr3OOaxXHA+eNTzg/scff5zk5GTWrl0LwKJFi1i+fDnJycmVIbKvv/46Xbp04fDhw0yYMIHzzz+fhISEI+bZunUr7777Lq+88goXXXQR//rXv7j88surna9r166sXr2aF154gaeeeopXX32Vhx56iOnTp3PPPffw5ZdfHqGsPKxdu5Y9e/aQnGz+43hMatdddx0vvvgiQ4YMYdmyZdx0000sWLCAs88+mzPPPJMLLrigdn84i6Wpkpde9TllAfQc1XiyNDGsImlCTJw48Yg8i2effZZPPvkEgN27d7N169ZqimTAgAGMHj0agHHjxrFjxw6/c//kJz+pHPPxxx8D8O2331bOP3v2bDp37lztuIEDB5Kamsqtt97KGWecwaxZszh06BBLlizhwgsvrBxXXFxcy6u2WJoJHkXSoQekLoTJtzeuPE0Iq0igxpVDQ9K+ffvKz4sWLWLevHksXbqU2NhYpk6d6jcPo23btpWfIyMjK01bgcZFRkZSVlYGmETBYHTu3Jl169bx1Vdf8fzzz/PBBx/wzDPP0KlTp8rVlMXSKsjbAxHRMOJ8WPEqlB6G6HaNLVWTwPpIGomOHTuSn58fcH9ubi6dO3cmNjaWzZs38/3339e7DJMnT+aDDz4AYO7cueTk5FQbk5WVRUVFBeeffz6PPPIIq1evJi4ujgEDBvDhhx8CRiGtW7fO1XVZLM2WvHSI6wmDp0N5Mexc0tgSNRnCqkhEZLaIbBGRbSJyt5/9bUXkfWf/MhFJcrZPFJG1zmudiJzndUwnEflIRDaLyCYROT6c1xAuEhISOPHEExkxYgR33nlntf2zZ8+mrKyMUaNGcf/99zNp0qR6l+HBBx9k7ty5jB07li+++IKePXvSsWPHI8bs2bOHqVOnMnr0aK6++moee+wxAN5++21ee+01jj32WIYPH85nn30GwCWXXMKTTz7JmDFjrLPd0rLIS4e43tDvBIhsY/wkFiCMPdtFJBL4ETgFSANWAJeq6kavMTcBo1T1BhG5BDhPVS8WkVigRFXLRKQnsA7o5Xx/A/hGVV8VkTZArKrWmFQxfvx49W1stWnTJoYOHVqPV9z8KC4uJjIykqioKJYuXcqNN94YVnOV/ZtbmjV/GQ29x8IFr8MbZ0PhAbjxu8aWKqyIyCpVHR9sXDh9JBOBbaqa6gj0HnAO4J1kcA7wO+fzR8BzIiKqWug1JgZQZ4444CTgagBVLQFKwncJLZtdu3Zx0UUXUVFRQZs2bXjllVcaWySLpWmialYkQ8803wdNg3m/g/wM6JjYqKI1BcJp2uoN7Pb6nuZs8ztGVcuAXCABQESOE5ENwHrgBmf/QCAT+LuIrBGRV0WkPX4QketEZKWIrMzMzKzP62oxDBkyhDVr1rBu3TpWrFjBhAkTGlski6VpUpht/CJxzi1s0HTznrqo0USqpLwUShu3IGo4FYm/Kn2+drSAY1R1maoOByYA94hIDGYFNRb4m6qOAQqAar4X5/iXVXW8qo7v1i1o73qLxWIJTF6aefcoksSRENu1afhJPrsF3jy7UUUIpyJJA/p6fe8DpAcaIyJRQDyQ7T1AVTdhFMYIZ3yaqi5zdn+EUSwWi8USPjw5JB5FEhEBA6eafJIw+ZldUV4GWz6H3cvg0P5GEyOcimQFMEREBjhO8UuAOT5j5gBXOZ8vABaoqjrHRAGISH/gaGCHqu4DdovI0c4xMzjS52KxWCz1j6c8Slyvqm2DpsGhDNjfiLeg9NVQ7FTlaEQzW9gUiePTuAX4CtgEfKCqG0TkYRHxrMNeAxJEZBvwK6rMVJOBdSKyFvgEuElVPUWgbgXeFpEfgNHAH8J1DRaLxQKYFYlEQofuVdsGTjPvjWneSlkACLSNh5SFjSZGWPNIVPVzVT1KVQep6qPOtgdUdY7zuUhVL1TVwao60RPhpapvqepwVR2tqmNV9VOvOdc6vo9RqnquqlbPomsG1KWMPMAzzzxDYWFh8IFBWLRoEUuW2MQqi6VG8tKhY0+IiKzaFt8buh7dqDdwUhaYkOTB0xvVzGYz2xsJq0gslmZE3p4jzVoeBk2Hnd81TtTU4YOQttLIMHAa5O+FzM0NLwdWkTQavmXkAZ588kkmTJjAqFGjePDBBwEoKCjgjDPO4Nhjj2XEiBG8//77PPvss6SnpzNt2jSmTZvmd+5hw4YxatQo7rjjDsCUjD///POZMGECEyZM4LvvvmPHjh28+OKLPP3004wePZpvvvmm4f4AFktzIi89gCKZBmVFsLv+SxgFZcc3oOVGkQzymNkaZ3VkizYCfHE37Ftfv3P2GAmnPR5wt28Z+blz57J161aWL1+OqnL22Wfz9ddfk5mZSa9evfjvf/8LmBpc8fHx/PnPf2bhwoV07dr1iHmzs7P55JNP2Lx5MyJSWfb9F7/4Bb/85S+ZPHkyu3bt4tRTT2XTpk3ccMMNdOjQoVLhWCwWHzzJiENOrb6v/4mmkGPKAhPF1ZCkLIA2HaDPBIiMhi6DjHnr+JsaVg6sImkyzJ07l7lz5zJmzBgADh06xNatW5kyZQp33HEHd911F2eeeSZTpkypcZ64uDhiYmL4+c9/zhlnnMGZZ5pM3Hnz5h3RuTAvL88WV7RY3FB0EEoL/a9I2naAvseZlcApDSxXygIYcJJRImBWJWvfhbISiGrToKJYRQI1rhwaClXlnnvu4frrr6+2b9WqVXz++efcc889zJo1iwceeCDgPFFRUSxfvpz58+fz3nvv8dxzz7FgwQIqKipYunQp7drZstcWS0hU5pD4USRgbuALHoFDmdChgZKfs1MhZwccf0vVtoHTTHn7tOWQNLlh5HCwPpJGwrfc+qmnnsrrr7/OoUOHAFN1d//+/aSnpxMbG8vll1/OHXfcwerVq/0e7+HQoUPk5uZy+umn88wzz1SazmbNmsVzzz1XOc6z3ZZ9t1iC4JuM6IvHP7F9ccPIA1Uhx55SLQADppgQ5Ubwk1hF0kj4lpGfNWsWl112GccffzwjR47kggsuID8/n/Xr1zNx4kRGjx7No48+yn333QeYVrennXZaNWd7fn4+Z555JqNGjeLkk0+u7LH+7LPPsnLlSkaNGsWwYcN48cUXATjrrLP45JNPrLPdYgmEv2REb3qOhnadYdt893PmpsGXv4Wi3NrJlLIQOvWDLgOrtsXEQ5/xxk/SwIStjHxTwpaRbxrYv7mlWbLwD7D4Cbg/s8of4ctnt8D6D+GWldCpr/8x3rx/BWyaAyfcBrMeCU2e8lJ4YiCM+Amc9RcfWR+Dr5+AO1Mgtkto8/rBbRl5uyKxWCyWmsjbAx0SAysRgKlOUY75Dwefb9f3RonEdoVlL8HBXaHJs2eVKYvibdbyMGgaaAVs/zq0OeuIVSQWi8VSE4FySLyJ7wPH3wzrP4A9qwOPU4W595ks+Z99CSIwP8QVScoCkAgTseVL73HQpmODm7datSJpDWa9poL9W1uaLW4UCcCJt0P7bkZRBPr3vvFTSFsB0+6FrkNg0o1G+aSH0Jk0ZYFRGO06V98XGW2c7g3scG+1iiQmJoYDBw7YG1wDoKocOHCAmJiYxhbFYgkdT6/2YMTEwdR7TMmULZ9X319WbLoqdh8Ooy8z2yb/EmIT4H/3u6uTdTjHmLb8mbU8DJwGB3eaEOEGotXmkfTp04e0tDRs98SGISYmhj59+jS2GBZLaBTlGX+EmxUJwNirjN/jfw/AkFlH+lVWvGpyPy7/uKr4Y0w8nHwXfPEb2Po/OGpWzfNv/9r4QGpSJN7lUryjusJIq1Uk0dHRDBgwoLHFsFgsTZn8veY93uVDUGSUicJ65yJY+Toc5yQYH84xkV+DpsPgGUceM+7/YNmLRvkMmm7mCETKAmgbZ0xbgUgYDHF9jJ9kwjXu5K4jrda0ZbFYLEEJlkPijyGzjCN80eOmQi/A10+ZnJFT/DjWo9rAjAchcxOseyfwvKqwzacsij9EYNBUs3qpKHcvdx2wisRisVgCEaw8ij9EYNajZhXyzZ8gezssfxnG/BR6jPB/zLBzoM9EWPAolBT4H5OdCrm7qkxXNTFwmlFc6Wvcy10HrCKxWCyWQHgUSceeoR3Xc5RxqC97Ef59G0REmUitQIjArN/DoX2w9Hn/Y/yVRQnEwGmANFj3xlbrI7FY6oXSIoi20Wgtlrw9JqQ3qm3ox06/D5I/Niamk34TfFXT7zgYehZ89xcT2is+z/nr3oXOSe4c6O0TjDJLWQgn/yZ02UPEKhKLpbZkbYW/nQBXzoH+xze2NJZwkBugM6Ib4nrBjPthzT/hxNvcHTPzIdj+DXweoD+Qd7XfYAycBkufg+J8aNvR/XG1wCoSi6W2bPkcyktMUzSrSFomeenQuX/tjz/+Zph0kzFduSFhEPxqIxQf8r+/Q3f35z72EtNgLyL8t3mrSCyW2rJtnnn3RPZYWh55e+r+kOBWiXho09686kr3oebVAFhnu8VSG4oPmeJ7UOWQtbQsSgpMd8TamrZaEVaRWCy1Yce3xqwVFWMVSUslz0lGdFMepZVjFYnFUhu2zYPoWBg805q2Wiq1SUZspVhFYrHUhm3zIGmKCcXMS3dXcM/SvAjWYtdSiVUkFkuoHEiBnO1mNRLXG8qLoTC7saWy1DeeFUmoyYitEKtILJZQ8WQLD54B8c7Tal5a48ljCQ956SYxsE1sY0vS5LGKxGIJlW3zqjKMPfZz63BvebjtQ2IJryIRkdkiskVEtonI3X72txWR9539y0Qkydk+UUTWOq91InKez3GRIrJGRP4TTvktlmqUFZvM48EzTX6A50ZjHe4tj7w6ZLW3MsKmSEQkEngeOA0YBlwqIsN8hl0D5KjqYOBp4I/O9mRgvKqOBmYDL4mId/LkL4BN4ZLdYgnIru+htAAGOT0l2nczmcN2RVK/HNoPr0w3lXPdsnt5aC1rg+G2xa4lrCuSicA2VU1V1RLgPeAcnzHnAG84nz8CZoiIqGqhqpY522OAypAYEekDnAG8GkbZLRb/bJsHEU5fbDCd7jr2tIqkvtmz2rSUTZnv/phPb4R3L4GSwrqfv7QICrOsacsl4VQkvYHdXt/TnG1+xziKIxdIABCR40RkA7AeuMFLsTwD/AaoqOnkInKdiKwUkZW2na6l3tg2H/pNOrIIXlwva9qqbzx/z33J7saXFJhouvy98P0LdT9/vk1GDIVwKhJ/BWZ8g+0DjlHVZao6HJgA3CMiMSJyJrBfVVcFO7mqvqyq41V1fLdu3UKV3WKpTt5e2L/B+Ee8ietlVyT1jefvmeFSkWRsBBTad4dvn4GCrPo5vzVtuSKciiQN6Ov1vQ/g+7+tcozjA4kHjgjIV9VNQAEwAjgROFtEdmBMZdNF5J/hEN5iqYbHzOLbczuutyk3bpMS649KRbIRKmo0Phj2/WDez/sblBaa/uj1cX67InFFOBXJCmCIiAwQkTbAJcAcnzFzgKuczxcAC1RVnWOiAESkP3A0sENV71HVPqqa5My3QFUvD+M1WCxVbJsHHRIh0addalxvKDtsWqta6gePaau0wCR/BiMjGWLiTRDE2Cth5WvG1FXX88fZZEQ31KhInDDbWj3xOz6NW4CvMBFWH6jqBhF5WETOdoa9BiSIyDbgV4AnRHgysE5E1gKfADepah3XqhZLHagoN93mPGG/3thckvonLx069TOf3Zi39iUbBS8CU++ByLYw/+G6nb9tfNgbQrUUalQkqloOdHNWFCGjqp+r6lGqOkhVH3W2PaCqc5zPRap6oaoOVtWJqprqbH9LVYer6mhVHauqn/qZe5GqnlkbuSyWkNmz2pQU9zVrgVcuiVUk9YKqWREMnmnazWZsqHl8RYUZ02Ok+d4xEU64FTZ+CmkrayeDzSEJCTemrR3AdyJyv4j8yvMKs1wWS9Ni2zxATPtSXypXJDZyq14oOmj8HF0GQcLg4JFbOduNCczb5HjCLcbxPvf+2vmubA5JSLhRJOnAf5yxHb1eFkvrYds86D0OYrtU39ch0Tw52xVJ/eAdMZU4HDLW1zx+n7O/h5ciadsRpt4Nu5bAli9qJ4NVJK4J2mpXVR8CEJGO5qsGaCZssbRQCrMhfTWc9Bv/+yOjbFJifeIdMZU4AjZ8AkW5xpnuj4xkkEjo5tNWduyV8P3fYN7vYMgs8zu5obwUDmXYiK0QCLoiEZERIrIGU7Zkg4isEpHh4RfNYmkibF8MWuHfP+IhrpetAAzm6X/vD3Wbw2MijO9d5ffI2Bh4/L5k6DoEomOO3B4ZDTN/B1lbYG0IMUO5uwG1K5IQcGPaehn4lar2V9X+wK+BV8IrlsXShMjcAgj0GBV4jE1KhKxt8P7lsPiPwcfWRF66MRV2SDSmLag5cisjuUrh+HLMGdB3Eiz8AxS7NKaseM2cv/+JocndinGjSNqr6kLPF1VdBLQPm0QWS1MjZ4dRFL5PvN7YpET43wNQUVb3oIO8PUaJREabv2tMp8CKpDDbrCB8c3s8iMCsR4ypaulzLs6dDitehVGXQNfBtb+GVoYbRZLqRGwlOa/7gBBKcloszZzs7dB5QM1j4nqZyKHivIaRqamx/WvY8l+Ibm9KydQFb0e3iFltBIrc8oQG9wigSAD6ToShZ8N3z0J+Rs3n/vopowxPDuAPs/jFjSL5GdAN+Nh5dQX+L5xCWSxNipwdppFVTbTmpMSKcvjqtxDfDyZea57+y0trP59vxFTicNi/0ZzHF89KJTGAacvDzN+ZlsiLHgs8JmcHrH7TOOm7BHlwsBxB0Mx24LeqepuTGDhWVW9XVVsLwtI6KCmEQ/ugS1LN4+L6mPfWmEuy9h0TgjvzQdM1EoX8fbWfz7czYeIIk1eSs6P62H3JJl+kY2LNcyYMgvHXGEWRucX/mMVPGN/ISXfWWvTWipvM9nENJIvF0vTw3LzcmLbA+ElaE8WHYMEj0GcCjDi/7ln+RXnGPOi9IvGYrfb5ySfJWF+zWcubk++CNu3hfw9W35e1Fda9CxN+bqO1aoEb09YaEZkjIleIyE88r7BLZrE0Bdwqko49AGl9pq3vnjGmrFMfc1oPO0UOa7sy89cHpNsxTqkUHz9JeSns3xTY0e5L+wSY/Ev48QvY8e2R+xb+AaLamf2WkHGjSLoAB4DpwFnOy9a4srQOPJVng9nMI6NNpFFrMm3lpsGSv5qVSN8JZltdfUW5aUfOAxDdDhKGVK+5lbUVyksCh/76Y9KNRknNva+qPP2+ZNjwMUy6ATrY3kW1ocZUT8dH8oOqPt1A8lgsTYvs7dA2Dtp1Dj62teWSzH/YhDvP/F3VtphOEB1btbIIlUANpXqMgN0rjtxW6Wh3uSIBo5Sm32fa8m74GEZeYFYjbeNNoUdLrXDjIzm7pjEWS4vGE7HlWzreH61JkaStgh/eh+Nvrir3Do55qw6thz1/v44+fUASh0PuLlMqxcOtOJM6AAAgAElEQVS+9aZcfNchoZ1j1MUmymv+Q7BziQlbPuFWdw8LFr+4MW0tEZHnRGSKiIz1vMIumcXSFMjZ7j4UNK533RRJWbHpPV6UZ5pkFWbDoUwozg9tni1fQuaPtZcjGKom3Ld9N5jipxB4XRRq3h4zb1TbI7d7wnu9zVv71kP3Y4xZMRQiImHWw3BwF7xzCcQmGLOWpda4qWJ2gvPu3SVGMT4Ti6X58fWTJmR03FU1j6soNzebY85wN298byjONTf+UBsiLX3e3Jz9ERENt62BTn397/emohw+vBoGTYNL3w1NBrdkboHd38PsP/q/zrje1Z3ZbglUdbcycisZ+ju3pIxkGHJq7c4zaLrpppgyH2Y9ahtY1RE31X/9NGCwWJoxy1+FDi4USV66ceYGi9jy4B362u3o0GRK/pfpvTHmCvPELJHmvSALvn4C9qx0p0iyU03b3+3fQHmZ+4q3oZC6yLwfc7r//R17Gh9JRbm5hlDIS/ef/NmxpzE9efwi+RlQkOk+9NcfZ/4ZVv0DJlxT+zksgAtFIiKJwB+AXqp6mogMA45X1dfCLp3FUt94EgwPZ5vw0ZrMIp6IrWBZ7R68G1yFokgO50D6GpMIN/n2I/eVFsE3fzImneHnBZ/LY/opyYc9q6Dfce7lcEvqQpN46O0b8SaulykzUpDphEWHQN6eqhWHNyLGqe5RJJ4eJaE42n3pnHRkoICl1rjxkfwD03fds978Ebg94GiLpSnjyQspLzHho27GuvaR1DL0dfs3pky9v+6L0THGmRys3ayHjA0m5wIxN/z6przUmK0GTg08prZJiSUFpjtioITAHiNNOfmK8qraW3VZkVjqDTeKpKuqfgBUAKhqGeCn6I3F0gzI8ao3WlNpcjChvxFRVeVPguGJNAr1Bpq6ENp0hD7j/e9PHB5cVg/7NxoTWa/RVSao+mTPKig5FESR1FKh5vlJRvQmcbgx22VvN472+L420qqJ4EaRFIhIAsbBjohMAnJrPsRiaaJkO4pEIv2X3PAmZ4e5Wbn1M0S1NRFHoYa+piyEpMmBzWyJw43Tv8jFf7uMZOg+zNzo01aEHvEVjNRFgMCAkwKPqe2KxPN3C7Qi8ZixMtab66yLWctSr7hRJL8C5gCDROQ74E3AZu5Ymic5203LVjdP+aGE/noINQQ4Z4c5z6AaYlo8N8z9m2qeqzjfzJc4wpjJKspqHz0ViJSF0GtMzSuB2ASIbBO6Qg2UjOih2zHmASBtpTFLWrNWkyGoIlHV1cDJmDDg64HhqlrHXpoWSyPh6S1SU4+LI8YmhTZ/qIokxfFjDJwaeIynS2CwFdT+zVXj+x5nakfVp3mrKM+scgZOrXlcRIRxstf3isTjL9rwCWh5aKVRLGHFzYoEVS1T1Q2qmqyqdWg0YLE0Mp5VRuIIKNgPh/b7H3c4xzh+3Yb+eojrVVUvyg2pC6FjL+h6VA1zeroEBnG4V5YMGWZuuv2Pr1JU9cHOJeYGPnBq8LG1Sc7MS4d2XUwZk0AkjqhSONa01WRwpUgslhZBeZnxNXQeUHNpcgg9YstDXC+jgEoKgo+tKIfUxcasVVMJlsrQ1yCKZP9G47SPd8JyB06DrC31V9o+dZFZ5fR1EVJcmzIpeXsCO9o9eH63Nh1CV/KWsGEViaX1kJdm/AaeFQkE9pNUlo9PCu0clY5mF0UL9641Ssdf2K8vlV0CKwKPydgA3Yca0xJU+V22Lw4+vxtSF5pVTk296z3E9TJJiaH0sM/bE7wXiOd36z6s6jotjU7AX8K7rpa/V0MKabHUC56Irc4DILaLuekH8pNkh5iM6ME7KTEYHv/FwKnBxyYON2G3B3f6369qFInHnwLQfTjEdg1u3lr/Ebx+Ws2rqLy9kLnZnaxg/rZlRcZE6JZA5VG88SgS62hvUtQU1/gn5z0GGA+sAwQYBSwDJodXNIulnvHtLeKdKe1vbGzX0GswhZJDkbLQFCN00wOjcgW1wb+5LS/drG68FUlEhLnxpy4yisaf+aw4H764Cwqz4Pu/wUl3+D+/Z1UzcGpwWeFIhRrbJfj40iIoPBDctNWxh6kAcIxtidSUCLgiUdVpTp2tncBYVR2vquOAMcA2N5OLyGwR2SIi20Tkbj/724rI+87+ZSKS5GyfKCJrndc6ETnP2d5XRBaKyCYR2SAivwj9ki2tluztpux4R+cm12MEZP1oqu76krMjdP8IeJm2gqxISgph9zIYNNXdvN2PASSw4vP4T7wVCRjzVsF+Yxbzx5LnjBJJHAnf/cVUHPZH6iIT1pvoMlKqY4hJiflBQn89iJh+Ir1Gu5vX0iC4MTIeo6qVHklVTQaC/opOU6zngdOAYcClTp0ub64BclR1MPA08EdnezIwXlVHA7OBl0QkCigDfq2qQ4FJwM1+5rRY/JOzHTr3r7KtJ44wPpPMzdXHZu+onTM3OsbccIPdQHcuMWVaBk51N2+b9qa+VSBFst9RJN19/jt45vdn3jq0H5Y+B8POgZ+8bFYn3/rpYadqjh9wsnu/RCgmPgieQ2Jp0rj5V7FJRF4VkakicrKIvAIEyYwCYCKwTVVTVbUEeA84x2fMOcAbzuePgBkiIqpa6JRiAWNaUwBV3evktaCq+Y4cQdbCFouDr3Lw5CH4+knKSoxjPlT/iAc3EUupC03SXj8/BQoD0aOGyK2MDaaUS7tOR26P72Pa1PrLJ/n6SSg9DNMfMCHDoy6G5S9XV4KZW0yhy4FT3cvaIdHU/HK7IvGMi3dZjsbSpHCjSP4P2AD8AlOscaOzLRi9gd1e39OoftOvHOMojlwgAUBEjhORDcB64AYvxYKzPwljZlvm7+Qicp2IrBSRlZmZmS7EtbRoVJ0VSVLVti4DTTir71N+7m5TRLE2pi1wciiCKJKUhdBvErSJdT9v4ghjnis+VH1fxsbqZi0Pg6bBzu+ONOFlb4eVf4exV0LXwWbbtHtMSPLiPx55fChBAR4io6BDD3fRa1D19/LtjGhpFrjJbC8CXgTuVtXzVPVpZ1sw/AXG+8YCBhyjqstUdTgwAbhHRCpjDkWkA/Av4HZVzQsg98uOX2d8t24unJmWlk1Blol68lYOEZHmSdw3l8Q7uqs2BOsQmJ9hTFFuwn69SRwOaHVTXFmJyRcJpEgGToXSQti9vGrbwkdNQcqT76ra1jkJxv8MVr8FWV5u0NRFRul27h+avKHkkuSlm9I1bTuEdg5LkyCoIhGRs4G1wJfO99EiMsfF3GmAdyeePoDv/67KMY4PJB44wtunqpuAAmCEMy4ao0TeVtWPXchhsXj1FvFRDp7ILe98h1D7kPgS18tEIJUGeN7yPOHXVF/LHx5F4buCyvrR+HoCKZKkyaZGlee8e9fB+g9h0o0Q57MCOOkOiIoxigbclY0PRCgtd/PSg0dsWZosbkxbD2L8HQcBVHUtkOTiuBXAEBEZICJtgEswxR+9mQN42tRdACxQVXWOiQIQkf7A0cAOERHgNWCTqv7ZhQwWiyHbJ/TXQ4+RJtfB+4aXs8PcTENtyuTBU3Y+P8BNNHWRKQXS49jQ5o3vZzLXff0knoisQIokJt6UqPf0J5n3kCm6eKKfoMcO3eH4m2DDx5C+1ikbn98AisRFMqKlyeJGkZSpashl4x2fxi2YplibgA9UdYOIPOyscsAohQQR2YapMuwJEZ4MrBORtcAnwE2qmgWcCFwBTPcKDw7Q79Ni8SJnOyDQycc84y/DPWeHWY3UVLakJmrKJVE1N/SBIURAeYiIcExxPiuSjGTT1z1hcOBjB04zXRg3zjF9yqf8urpj3sMJtxpFM//hqrLxSVNCkxXM36Ek3xR7DIabZERLk8VNo4VkEbkMiBSRIcBtwBI3k6vq58DnPtse8PpcBFzo57i3gLf8bP8W/34Vi6VmsrebG5VveQ/vyrpHnVo1ti51nDwmGn81rjK3mNIhA6fWbu7E4bD+X0cmGGZsNCXWa2obPHAqLH4cPr3RrJgmXBt4bEw8TP4V/O9+83fpNcZdUqEv3n1JYuICjysrMaHI1rTVbHHzSHQrMBwoBt7BRFbZVruWurFtXvAy7vVJTgDlEBNnVimeFYlq7ZMRPXj8Dmvegm+fMSVIdn0PB3eb64bQHe0eEkdAce6RFYZ9S6P4o894YxYrOQTTfhu8XtbEa01SYcH+2is9t7kk+XsBtSuSZkyNKxInqfAhVb0TuLdhRLK0eMpL4YOrjA/ipmXuOxDWheztcNQs//u8e5MUZEJpQe0d7WCSB485E7Z/DTu+qb6/NhFQHrxLpXTqazLR89ONyasmIqPN9Wf9CMdeEvw80e1MOPCcW2HwzNrJ6lEM+UFCgG0yYrOnxv/BqlouIuMaShhLKyFtpXkyPrAN1r4N464KfkxdKD5knqwDmasSR8CWz03ZkrqG/nq45G3zXpRnnshz95gkx9w9poJubek+1LxnJMPRs4M72r0590XTTyQi0t25xlwBvcbWvkCi2x72lQ2trGmrueLmUXCNE+77ISYMFwAbemupNakLTdZz4nBY9DiMuqjmZkZ1JVhvkR4jTALi/k21Lx8fiJg48/IogPqYr1P/qsgtz3t3F4okqk1o5xKpW5XdqLam8GUw05ZdkTR73PhIugAHgOnAWc7Llt601J6UheZJd/bjxiyz/JXwni9QDomHSnPR+qrortqanhoC76rFGRtMKHFtQ5XDjZsQ4Lx006iqbQ0OeUuTJuiKRFXdlEOxWNxRlGtyE6b8yiTKDZoB3/7ZmLdi4sNzzkA5JB469TeO6H3JpnBhXG/zNN1USRwOP35h6mR5HO21DVUON3G9TcmZmvDkkDTVa7AExU1me4yI3CwiL4jI655XQwhnaYHs+PbIvt8zHjAJgUv+Gr5z5mw3Pc/bdfa/P8Ixs2UkV+WQNGUShzumuI3GHOfGP9JYuCmTYnNImj1uTFtvAT2AU4HFmFIn+eEUytKCSVkI0e2hz0TzvddoGP4TWPq8ySUIB9nbg4fzeirrZqdCl6TwyFFfeKoWb/7cRJg1dUVyOMesngKRl15VDcDSLHGjSAar6v1Agaq+AZwBuOxuY7H4kLoIkk480vE7/T5TmfbrJ8NzzkA5JN4kjoDiPCe6Kyk8ctQXnZMgOhbWf2C+u3G0NxbeSYn+KC8zJertiqRZ40aRlDrvB0VkBKawYlLYJLK0XHLT4MDW6gluCYNg7BWmrLknaqq+KC81iYBBVyRez0Z1Df0NNxGRJgrs4C5AnO6JTZS4ICHAhzKMmc4qkmaNG0Xysoh0Bu7HFFncCDwRVqksLRNPlz5/Wd0n32VukAsfq99z5u42PplgyqH7UCqr79Qlq72h8JizugwwCZBNlWArksrQX5tD0pxx04/kVVXNUdXFqjpQVbur6osNIZylhZG6yHTO85dTEdcLjrsefng/cBfA2hAsYstDm/ZmZQRNf0UCVSHLTdk/Al5JiQEc7pXJiHZF0pwJGv4rIg/4266qD9e/OJYWS0WFUSSDZwQO8zzxdlj5D1j4h6rM8LoSLIfEmx4jTYmUQNFdTQmPAkmsQ8JgQ9C2gwnrDrQi2b3cJKfaFrvNGjeZ7QVen2MwyYhuerZbLFVkJENhVs3FCmO7wMgL4IcPjOIJtcy6P7K3Q2Rbdy1cp91nyoI0h3yGXmNMDs4xZzS2JMGJ6+2/3lZeOqx8DUZdErikvaVZ4CYh8U/e30XkKao3qLJYasZt3++eo8zN5eDO+vFVePJC3CilroOr+pc3ddq0hyuaSZWiQLkkXz9pesRPvav6PkuzojaPfLHAwPoWxNLCSV1oemb4tnb1JdGJnvJtJ1tb3OSQWMJLx57VTVvZqbD6TRh3ddMPt7YExU1m+3oR+cF5bQC2AH8Jv2iWFkNpEexc4q4HR/ehxmZeH71KPL1FmoPzvCUT19skm5aVVG1b9Ljp6njSHY0nl6XecOMj8S7QWAZkOG10LRZ37F4GZUUwyIUiaRMLXQbVz4rk0H6T+W1XJI1LXC9ATeJhp36mo+MPH5ie8U212KQlJNwoEt9yKHHi5YxU1ex6lcjS8khdCBFR0P8Ed+N7jDD9xetKZUl4q0gaFe9ckk79YOGj0LajUSSWFoEbRbIa6AvkYDK2OgG7nH2K9ZdYgpGy0NTWatvR3fjE4bDhE9MUqqZe38HIcZlDYgkv3i1301bB5v+YCLna9IG3NEncONu/BM5S1a6qmoAxdX2sqgNU1SoRS80UZsPede7MWh48DndP97/aku30FunUr27zWOpGZZmUvbDgYdPsatINjSuTpV5xo0gmqOrnni+q+gVwcvhEsrQoti8GNHjYrzeernz71tft3DnbTaJbU+4t0hqI6VRVZDJ1EUz5tfvVqaVZ4Ma0lSUi9wH/xJiyLsd0TLRYgpOyENrGm46IbonrbW4+dXW4Z2+3oaVNARFj3tq7zvy243/W2BJZ6hk3iuRS4EHgE+f7YmebxRKc1EUwYApEuvmn5iBiypW4CQFWhV3fmxLwFeWmQKNWmM8HtsLQs2otuqUeiesFB7aZ4pzRMY0tjaWecZPZng38AkBEIoH2qpoXbsEsLYDsVJOhfsKtoR+bOAJWv2EUQkRk4HHbF8Ob59Qwj22d0yToOdr4y0Zf1tiSWMKAm6KN7wA3AOXAKiBeRP6sqmHqQmRpMXjKxg+aHvqxicOhtNCE8Hqq8vrjx7mmltZV/4bIaKN0JNK8R7aBLjYepElwysMw48HQVqaWZoObX3WYquaJyE+Bz4G7MArFKhJLzaQsgPh+tbuZezvca1IkKfOh//HQ77jayWhpGESsEmnBuInaihaRaOBc4DNVLcU43YMiIrNFZIuIbBORu/3sbysi7zv7l4lIkrN9ooisdV7rROQ8t3NamgjlZbD9Gxg0tXbVdLsNNSuLmhzuuXsgc7OpgmuxWBoNN4rkJWAH0B74WkT6A0F9JI4/5XngNGAYcKmIDPMZdg2Qo6qDgaeBPzrbk4HxqjoamA28JCJRLue0NAXS10Bxbu3MWmAcsl2H1OxwT1lg3gdbRWKxNCZuOiQ+q6q9VfV0VVVMVrub7LKJwDZVTVXVEuA9wNcreg7whvP5I2CGiIiqFnrV84qhagXkZk5LUyBlASAwoA4pR4kjal6RpMyHDj2gu32WsFgak5DLyKvBTdHG3sBur+9pzja/Y5w5c4EEABE5zqk2vB64wdnvZk5LUyB1IfQaXbcyGD1GmJ7rh3Oq76soN878QdObRyMqi6UFUw8t6ALi73+3r28l4BhVXaaqw4EJwD0iEuNyTjOxyHUislJEVmZmZoYgtqXOFOVB2gp3ZeNrwtNGNsNPqZT0tVB00Jq1LJYmQDgVSRqm2KOHPoBv4+bKMSISBcQDR1QTVtVNmHa/I1zO6TnuZVUdr6rju3XrVofLsITMjm+hoqz2/hEPlYrEj3krZT4gdVdWFoulzriKxxORE4Ak7/Gq+maQw1YAQ0RkALAHuATwzUaaA1wFLAUuABaoqjrH7FbVMse5fzTG4X/QxZyWxiZ1oamt1Hdi3ebp2ANiE/zX3No235jO2ifU7RwWi6XOuElIfAsYBKzFJCWCMSfVqEgcJXAL8BUQCbyuqhtE5GFgparOAV4D3hKRbZiVyCXO4ZOBu0WkFKgAblLVLEeeanOGcsGWBiBlIfQ/se7FEkX8O9yLco3pbPIv6za/xWKpF9ysSMZjkhJd5Y5441QN/txn2wNen4uAC/0c9xbwlts5LfXMoschfy+cVYuOygd3mxpX9VWYr8dIWPGqyUvxJLSlLjY1tepqOrNYLPWCGx9JMmD7YbYmfvgAVv0Dtn8d+rGpnrIo9eS7SBxh2vRmp1RtS1kAbTrW3XRmsVjqBTeKpCuwUUS+EpE5nle4BbM0EiUFptgiwNz7oaIitONTFkLHntDtmPqRJ3G4eff4SVSNo33ASaa2lsViaXTcmLZ+F24hLE2IjI2AwrBzYONnkPwvGFXN+uifigpTNv6o2fWX29HtaNPvPWMDjLwADqTAwV0trt93SVkFOYUlJMbZEuuW5oebMvKLG0IQSxMhw3nyn/mQWZnMf9j09HDTQ2LfOjicXX9mLTAO+65HVzncU+ab9xZUX6u0vIKrXl/O6l05fHjD8Yzq06mxRbJYQiKoaUtEJonIChE5JCIlIlIuIrYfSUtlX7LpaNg5CU55BHJ3wYpX3B3rKRs/cGr9ytRjRFXNrW3zTTXhLgPq9xyNyO//s5GlqQdo1yaS695cxf68osYWyWIJCTc+kucwHRG3Au2AnzvbLC2RjGTjlxAxK4vBM+HrJ01TomCkLDDO8Q7d61emxBGQnw75+2DHNy1qNfLu8l28sXQn104ZwDs/n0Tu4VJu+OcqisvKgx9ssTQRXGW2q+o2IFJVy1X178DUsEplaRwqKowvwtMLBExDouJ8+OZPNR9bUgi7l9WvWcuDR54Vr5lmVy0k7HfFjmwe+CyZk47qxt2nDWVYrzj+fNGxrN51kHs/SaYWEfcWS6PgRpEUikgbYK2IPCEiv8SUlLe0NA7ugJJDVaVJwKxORl8Gy1823QoDsXMJlJeEp2SJR54VrxjH+4Ap9X+OBmbPwcPc8NYq+nSO5a+XjCEywgQnnDayJ7fNGMJHq9J4/bsdjSukxeISN4rkCmfcLZiaV32B88MplKWR8PghvFckANPuNU2m5j8S+NiUBablbf8T6l+uDt2hfXdTBbjvJGjbsf7P0YAUlpRx7RsrKSmr4JUrxxMfe2QY8+0zhnDq8EQe/e9GvtlqC45amj5u+pHsxFTd7amqD6nqrxxTl6U5UF7qvwy7PzKSQSKq9/eI6wXH3wzJH8Ge1f6PTV1oWt5Gt6ubvIHwKLfBzdusparc+eEPbNqXx7OXjmFw9w7VxkRECH++aDRHJXbklnfWsD2roBEktVjc46bW1lnAU0AbYICIjAYeVtWzwy2cxQ/F+bB1rvFJlBWZV6nzXlIABfvh0H4oyDTvhx0n+SXvwDFn1Dz3vmRIGOxfGZz4C5Pt/umNJhkwKsa8omPMamX/Rhh1cb1fbiWJI8yqp5k72p9fuI3/rt/LPacdw7RjAgcltG8bxStXjufs577l2jdX8u9bJtOuTWQDSmqxuMdtQuJEYBGAqq719Fa3NDCq8P7lJunPF4kwFXfbdzOmoK5DTOHEDt3h+xdg83+DK5KM9dB7vP99MXFwxp/gq3tNCRWPEvM+/1Gza31pQRlzhckp6TEqfOcIM3M37OOpuT9y7uheXHfSwKDj+3aJ5dlLx3DFa8t5e9lOfj4l+DEWS2PgRpGUqWqu2C50jc+WL4wSmX4/jLqoalUQFWPKhQT6jfZvNMepBh5TlGsyxsddHfj8w881Lw+qUFYMZYfNqiQmrpYX5oJuR8H0+8I3f5jZsi+fX76/llF94nn8/FG4/f80ZUg3ThiUwEtfp3L5pP7ERNtViaXp4apoo4hcBkSKyBAR+SuwJMxyWXwpK4a595os7xN/AZ36mdVGTBxEtam5JMnAaZC3B7K2Bh7j6UKYOCLwGF9EjGmrXefwKpFmTk5BCde+uZLYtlG8fMX4kJXBbTOGkJlfzHvLd4VJQoulbrhRJLcCw4Fi4F0gD7g9nEJZ/LDsJVOyZPYfQi9WOHCqefdU5vWHpwRJKIrEEpTS8gpufmc1+3KLeOmKcfSID72W1qSBCUxM6sKLi1NtoqKlSeImaqtQVe9V1QlO69p7nT4ilobi0H5Y/AQMOdVkmodKlwGm5ElKDYpk33qzsojrVWsxLdV59L+bWJJygD/8ZCRj+3Wu9Ty3zRjCvrwiPlyZVo/SWSz1Q0AfSbBS8TZqqwFZ8IjxQ5z6aO3nGDgN1n9kwoH9rWgyks1qxPrC6o33V+ziH0t2cM3kAVwwrk+d5jpxcAJj+nXib4tSuGh8X9pEuSpKYbE0CDX9azwe6AN8gwn//ZPPy9IQ7F0Hq9+CidebSKzaMmgalORD2srq+yrKjY+kx8jaz9+MqKjQsJcfWbkjm/s+TWbKkK7cc1rde7OICLfNGMKeg4f5ZI1dlViaFjUpkh7Ab4ERwF+AU4AsVV1sS8s3EKrw5T0Q2wVO/k3d5hpwkgnR9ecnyU41K55W4B85cKiYM//6LRe/9D3788NjoV2WeoBr31xJn86xPHfpWKIi62f1MPWobozqE8/zC1MoKw+x4ZjFEkYC/gt3CjR+qapXAZOAbcAiEbm1waRr7Wz8DHZ+Z8Je29WxR0W7ztBrjH8/iaf7oG9plBZGTkEJP311GalZh1i/J5ez//odP6QdrNdzvL1sJz99dRmd27fh71dPqFb+pC6ICLdOH8Ku7EI+W5teb/NaLHWlxkclEWkrIj8B/gncDDwLfNwQgrV6Sovgf/ebVcLYq+pnzoHTYM8qkzPiTUayKYZYX+1xmyC5haVc/toytmcV8OqVE/jXjScQGSFc+OJSPl2zp87zl5ZXcN+n67n3k2QmD+nKpzefSFLX+q9tOnNod4b2jOP5hdsor6hunsstLOXtZTvZvM+2DLI0HDU529/AmLW+AB5S1eQGk6qJ8EPaQQ4XlxFRXkREaQGRZYeJKCsgsqwQVClrE09Z23jK28SjkW2OOFbKi4k+nEV0URbRh7OIKnZKlUgkKhEgEajzWSqcc5QVEVFeTER5EbHZG0k4uAuu+jdE1D4JraC4jJ0HCjmmR0ciBk2Db56CHd8emeW+Lxm6HmUyx1sguYdLueL1ZWzNOMQrV41n8pCuAMy55URuens1t7+/lk178/jN7GMqq/CGwoFDxdz09mqWbc/m+pMH8ptTazePG0SE26YP5sa3V/OfH9I5Z3RvVJV1abn88/ud/HtdOsVlFfSKj+GL208ivp3ta28JPzVltl+BqfZ7FHCbVyauAKqqLT4Drd3r0xhRvp0ICe6YLdC2HKQDxRpNguQRL4V1One5Cm9XzGTxN7H8tDSTKYO7EhHk5lRUWvPUhfYAACAASURBVM7GvXn8sPsgP+zJZX1aLtsyD6EKT198LOeNnGDKqKQsPFKRZCSbciotkPyiUq56fTmb9ubx0hXjOPmobpX7Ejq05Z8/P46H/r2Bl75OZUtGPn+5ZExIN99Ne/P4+RsryTpUzDMXj+bcMb3DcRlHcOrwHhyV2IHnFmyjsKScf36/kw3pecS2ieT8cX2YmNSFX3+4jgc/S+aZS8aEXR6LJaAiUdVWH1/Ybvxl7Dl8kPKoWCqiYimPiqU8uj0VUaaoYVRxLlEleUSVHCSqJI/o4oO0KS+iMKYLB2O6URrTldJ2Xc17TBdAQMsRrUC0ovKzRkZTEdmWisgY8x4VQ4VEk7Yti1UrdjN3Ywb9E2K5bGI/Lhzfly7t21BSVsGPGfmsSzvI+rRcfkjL5ceMfMocc0fXDm05tk88Z47qxZtLd7B4SybnjeljFIa3w70w22S9t0D/SEFxGf/39xUk78nlhZ+OZfoxidXGREdG8PtzRzK0ZxwPfraBS1/+njm3nOjKQZ6aeYjz/7aEjjFRfHD98Rzbt2F6rUdECLdMH8Jt767hno/Xc0yPjjxy7gjOHd2LjjFGCe48UMjT835kxtBEzjrW5gZZwoubWlutlj6n39mo5z9xcFdunzmEL5P38fb3u3jsi838ae6PDO7egW2ZhygpM5E78e2iGdUnnuuOHsioPp04tm88PeJiKus5pWQeYknKAVQVGTQNvvofHNwNnfq2uIz2g4UlrN9jFOvn6/eyeV8+z106hlnDe9R43E+P6098u2hueWcNH6xM47Lj+gU91x+/3IwAn958Ij3jw1Q+PwBnjOxJbmEJw3rFMbZf52q1u26eNohFP+7n3k/WMz6pc4PLZ2ldWEXSxGkbFck5o3tzzuje/JiRz9vf7yQls4CrhyQxqk88o3p3om+XdjUWATxhUAJz1qWTklnAYE8Hw9SFMPZKr2ZWzTOHpKi0nHeW7WL1rhzW78ll54Eqk2JSQix/vXQMp43s6WquM0b25O/9d/D0vB85d0wvYtsE/u+xamc2X23I4NenHNUoN+nICOGK45MC7o+KjODpi0Zz+rPfcMeH63jrZ8cFNY1aLLXFKpJmxFGJHXnonNBXDicMMs7lpSlZDJ40FDr0MH6SsVeaHu3tu5sCkM2M9IOHue6tlSTvyaN3p3aM6hPPxRP6Mqp3J0b2jg859FZE+O3px3D+35by6jfbuW2G/wRQVeUPn2+me8e2XDNlQH1cSlhI6tqeB84cxt0fr+f177bbMvSWsGEVSSugb5d29O7UjiUpB8xT7MCpsO1/UFFhepAkDm9kCUNn5Y5sbvjnKopKK3j1yvHMHFbd/1EbxvXvwqnDE3lpcQqXHdePrh2qR7LN3ZjBqp05PPaTkTWuWpoCF0/oy7xN+3niqy1MGdKNo3s07zbFlqZJWB3qIjJbRLaIyDYRudvP/rYi8r6zf5mnYZaInCIiq0RkvfM+3euYS53tP4jIlyLSNZzX0BIQEU4YlMDS1ANUVKgpl1J4APaugf2bgzraVZUl27K4+e3VXPLy0kavQPve8l1c+sr3dGgbxac3n1BvSsTDb2YfQ1FZBc/Or152v7S8gj9+sZlB3dpzYR3rZzUEIsLj548kLiaKX7y3ptF/O0vLJGyKREQigeeB04BhwKUi4tMMnGuAHFUdDDwN/NHZngWcpaojgauAt5w5ozDlWqap6ijgB+CWcF1DS+KEwQkcLCxl0768qrLyy1+F8mJI9O8fOVhYwqvfpDLjT4u57NVlLP4xk+9Ts/loVePUeiotr+CBz5K5++P1TBqYwGc3T2Zw9/p/wh7UrQOXTOjLO8t2VeuX/v6K3aRmFXD3aUPrrfRJuOnaoS1PXDCKzfvyefyLzY0tjqUFEs7/CROBbaqaqqolwHvAOT5jzgHecD5/BMwQEVHVNarqqQGxAYgRkbaYHBYB2ovxLscBtlaEC44f6PGTHICOPaDbUFj/odnpsyLZnlXArz9Yx3F/mM/v/7uJzu3b8OeLjmXlfTMZ3bcTLyxMobSBaz3lFJRw5WvLeXPpTq6dMqDey4/48ouZQ2gTFcGTX1XdeAuKy3hm3lYmJHVm5tDm5VOafkwiV5+QxN+/2+F3pWWx1IVwKpLewG6v72nONr9jVLUMyAUSfMacD6xR1WJVLQVuBNZjFMgw4DV/JxeR60RkpYiszMzMrOu1NHt6xMcwsFt7lqQcMBsGTYOKUohsY7LaHVSVG95axRfJe7lgXB8+v20K/7rxBH4ytg8x0ZHcNmOwqUC72l1ZkcKSsjpX2j1YaGpkrdqVw58uPJZ7zxgW9tVA944xXHfSQD5fv481u3IAeOWbVLIOFfP/7d15XFVl/sDxzxcQUFSQxQ0EFdxxS0SjzGWytEzbs2xvnPxNTrZnv3Is/bVM095MM9Nii1qalWubmTqWuWG4r7igiIlL4pKCwPf3xz0oIggGF+Te7/v18sU95z7nnOc5XPne5zznfJ8nrmhT5qlyzyd/HdCWay+I5JXvNvHW/NSqro7xIO7831jc/7Sif1HOWkZE2uG63HWvs1wDVyDpDDTGdWnrieIOrqpvOxNxJURERBRXxOskxYaxZOt+V2+i4DbgiFanzU+ybPuvbNxzmNFXteXZa9rTtvHpCQx6t6pPfGRd/jk/tdQMtOt3H6LL2Dlc9uoCPli4jUPHT5xznbOOneC295aSmnmEd25P4LpKHJcY2qM54bUDeP6rDWQePs7bC7bSP75huSaoqko+PsLfr+/IoE6NefGbjbyzYGtVV8l4CHcGknSgSaHlKM68DHWyjDP+EQwccJajgKnA7aq6xSnfCUBVt6jra+6nQJK7GuBpkmLDOZqTx6r0LGh6kas30rDDaWUmLE6jTqBfiU9DF2SgTdv/GzNWlnxV8fiJPB6YtIKgAD9q+fvy9Mx1dHv2ex7/bBWr07NK3K6wQ8dPcPu4pWz45cz0JpUhKMCPBy5twdLtB7j7g2Xk5Obz6OWtKrUOFc3XR3j5ho5c2aERz361nnE/bqvqKhkP4M57F5cBLUSkGbALGAzcUqTMDFyD6YuA64G5qqoiEgJ8CTyhqgsLld8FtBWRCFXdi2uOlPVubINH6d7cddVw0ZZ9dIlpAbdMhtDYk+/vO5LN12t2M6RbzFlva+3bpgGtG9bhH/NSGdQpstgEhS99u5GNew7z/l1d6d2qPqvTs5i4JI3pKzKYnLyTDlHBDOkWzcCOkdT0PzMp5ZHsXO4ct5S1u7L4161d6N26asYkburahHE/bmPNrkPc1j2G5hG1q6QeFcnP14fXbupEXp4yZtY6avie/eFGY0rjth6JM+YxHPgW1x/7T1V1rYiMEZGCaXrfA8JEJBV4CCi4RXg4EAeMEpEVzr/6zgD8M8ACEVmFq4fynLva4GlCg/xp06huoXGSPlAv5uT7U5LTOZGnDCklPYiPj6tXsnXvUb5cvfuM9xem7uPdH7dxW/cYerdyBYD2UcG8cF0Hljz5B8YMasfxE3k8/vlqEp+bw9Mz1rJ5z+GT2/+Wk8vd7y9jZXoWb97cmb4VfHvvuajh68PYq+Pp2rReiQ8oVkc1fH144+bOXNqmPqOmr+XDn7a7bg035ncQd085ej5ISEjQ5ORippj1QmNnrWP84jRWjb6MwBqnegL5+UrPl+bROLgmk++9sNT95Ocrl7+2ABH4ZsQlJ9NvZP12gstfW0CtAF++/EuPYnsb4BrUT077lQmL0/h69S/k5OWT2CyUId2imbR0J0u27ef1wZ0t4aCbZefmMWz8cuZt3HtGYlBjRGS5qiaUVq563AhvKkxSbBg5ufn87NyJVGDB5r3sPHCMW7vHlLDl6VwZaOPYtOcI36795eT6p6avOZlSvaQgAq6xlq5NQ3l9cGcWPdGHkf1bszvrGCMmrWDxtv28cmMnCyKVIMDPl7dvT+CNmzvToG4gz3+9ge7Pfc8Dk1JI3n7A7XPbG89wfud3MBUusVkovj7Coi37T+bgApiweAfhtf25vJQsuYUN6NCY1+ds5o25qfSLb8iMlRnMXJnBI5e1pENU2VOqh9UOYFjPWP7Uozk/pO7Dz0e4KM4SFlSWGr4+DOzYmIEdG7Npz2E+XrKDz5enM21FBq0b1uGJK9pU+o0OpnqxQOJl6gTWoH1kMD9t2c/DzrqMg8eYu2EPw3rG4u9X9k6qr4/w595xPDJlJR8tSuOl2RvpElOPYT1jS9+4GD4+Yn+wqljLBnV4emA7HuvXipkrM/jPf7dyx7ilXNm+EaMGtKVhcGCJ2x7JzmXWygwyso4X+35CTD0usd+vR7JA4oWSYsN4e8FWjmTnUjvAj0lLd6DAzYmlz8FR1KBOjXnj+82MnrGWIH9fXr2xU7VJHWJKVsvfj5u6RnN150jeWbCVN+emMn9jJg/2bcmdSU1P+x2vyzjEhCVpTE/ZxdGcknN5+fkIX97fwxJHeiALJF4oKTact+ZvYdn2A1wcF86kZTvp1TKCJqG1znlfNXx9GN4njsc+W8Xoq9oRHXbu+zDnrwA/X4b3acHAjpGMnrGG//tyPZ8tT2f0Ve3IOHiMCUvSSNlxkAA/H67q2Jgh3aLp1CTkjCf/DxzNoc/L8xk1bQ2T7+1eLTMDmJJZIPFCXWLq4e/rw6It+zmek0fm4WyeL+Mge3FuTGjChc3DflcgMtVDdFgtxt3ZlW/X7uGZmWu5+Z3FADSPCGLUgLZcf0HUWXOfhQb5M7Jfa0Z+sZrPlqdzQ0KTEsua6scCiReq6e9L5+gQftqyj7UZWUSG1KRXq/I98GdBxPOJCP3iG9KjRThfpOwiLqI23ZuHlrl3cWNCEz5N3snzX2/g0jYNqGe3GHsMu5jtpZJiw1mbcYiFqfu5ObFJsU+nG1OcoAA/busew4WxYed0icrHR3j2mvZkHTvBi99udGMNTWWzQOKlkuLCUHUNgN7Y1S4zmMrRplFd7kpqyidLd5zxLJOpviyQeKmOUSHUCfCjX3xD6tcp+ZZOYyraA31b0rBuIE9OXVNqBmlTPVgg8VL+fj5MvS+JZ68pfnZEY9yldoAfo69qy/rdh/hwUVpVV8dUAAskXiyufh2Ca7pvlkFjStIvviG9WkXwyuyN/FLCA4ym+rBAYoypdCLCmIHx5OYrT01bzZa9Ryz7cDVmt/8aY6pEdFgtHuzbkhe+3sCc9ZnUDvAjPrIuHaNCaB8VTMeoELutvJqwQGKMqTLDesbSp3V9Vuw8yOr0LFalH+T9hdvJcQbhu8TUY0i3aK5o3+i0aQ9M6VSVvHytlJRFNh+JMea8kpObz8ZfDrNo6z4mLd3J1n1HCalVgxu6RHFLtxiahQdVdRWrhRkrM3hrXiof3p1Ig7q/787Mss5HYj0SY8x5xd/Ph/ZRwbSPCmZoj+Ys2rKfiUt28P7C7bzzwzYujgvn6YFtiatvyR9Lcuj4CcbOWkej4EDCawe4/XgWSIwx5y0RISkunKS4cDIPHefT5J28++M2RkxawfT7LrJM0yV4ZfYm9h3JZtwdXSsla4X9Fowx1UL9uoEM79OC565pz9qMQ4xf7L5nUM71kr+qnjezSa7ZlcVHi7ZzW/cY2kcFV8oxLZAYY6qV/vEN6dkygpdnb2LPoYp/BiVt/1EuemEuM1dmlKl8Xr4y9KNkhn60vMqDSV6+8uTU1YQGBfDwZa0q7bgWSIwx1YqI8MzAduTk5TN21roK3/+ExWlkZB3nf79YTfqvv5Va/t0ftjJnfSZz1u9h7obMch07ZcevDP0omUenrGT8ou2s3HmQ7NySJwsr6uOlO1iZnsWoAW0q9WFjGyMxxlQ7TcODuK9XHK/O2cRNXffSo0XFTOF7/EQeU5ank9g0lHW7D/Hwpyv5eGj3EscZ1mZk8dLsjVzWtgGbM4/wwtcb6Nky4pzHbrJ+O8Hfvt3AJ0t3EFrLHwWmLE8HoIav0KphHTpEhXB1p0gSm4UWu4+9h7N58ZsNXBQXxsCOjc/p+OVlPRJjTLV0b8/mNAsP4q/T13L8RMnf2o9m55J5uGyXwL5avZuDv51gxKUtGH1VW5ZsO8C7P2wttuzxE3k8OHkFIbX8eeG6Djx2eSs2Zx7h85/Ty9wGVeXz5en0eXk+k5bu4O6LmvHfx3qz/KlL+fHx3vxryAXcc3FzgmvWYOaKDG78zyIemryCfUeyz9jXc1+tJ/tEPmMGxVf6DJTWIzHGVEuBNXwZM6gdt723lP/8dysjLm1x2vuqysxVuxk7ax05ufnMf6RXqZNpTVicRvPwIJJiwwD4fn0mL83eSI8WEbRtXPe0si9+s5FNe47wwV1dCQ3yp198QzpHh/DKd5sY2DGSmv5nf4By857DPDVtDUu2HaBzdAgf3ZNIu8anBsej6tUiql4t+rdvBMCxnDz+MW8zby/Yypz1e3i0X2tuSYzG10f4acs+pqbs4i994oiNqF3mc1hRrEdijKm2erSIYECHRvxzfirb9x09uX7r3iPc9t5S7v8khfDaARw+foJ/zEs9677WZRzi5x0HuaVbNCKCiPDcte0JqeXPA5NTTuv1/Lh5H+MWbuOOC2NOzi4qIjzRvw17DmUzbuG2sx5r+opdXPHGD2z45TDPX9uez4clnRZEilPT35dHL2/N1yMuoV3jYEZNW8O1by3k5x2/MmraGqJDa3Ff77jSTplbWCAxxlRrowa0xd/Xh7/OcF3iemX2Rvq99gMrdx5kzKB2zPrLxdzQpQnjF6Wx80DJg+cTl6QR4OfD9V2iTq4LDfLn79d3YNOeI7z4jWtWx4O/5fDIlJXERgQxsn+b0/aR2CyUS9s04N/zt3DgaE6xx5m1KoMHJ6/gguh6zH24JzcnRuNzDs96xNWvzcdDu/H64E7sOnica9/6iS17j/LMoHZVlkbGAokxplprUDeQhy9ryYJNe7n4b/N4Y24qV7RvyPeP9OT2C5vi6yM82LclPj7w0uzip/g9kp3LtJRdDOjQmJBap1/+6tWqPndcGMO4hdv4cfM+npy2hn1Hsnl9cOdiL1893q8VR3NyeXPu5jPe+2bNbkZMWkGXmHqMu7MrYb/zqXMRYVCnSL5/uCd/vLgZ/9Mrlt5Oz6gqWCAxxlR7t3WP4YLoEIJr+vHxH7vx2uDOp8382TA4kHsubsb0FRms2ZV1xvbTUnZxNCePId2ji93/yP5tiI0I4t7xyXy5ajcP9m1JfGTxl6JaNKjDTV2bMGFxGjv2n+oBfbduD8M/TqFjVDDv35VIUED5h6iDa9bgqQFtebxf63LvqzzcGkhEpJ+IbBSRVBEZWcz7ASIy2Xl/iYg0ddb3FZHlIrLa+dmn0Db+IvK2iGwSkQ0icp0722CMOf/5+fowZVgScx7qSVJceLFl7u0ZS71aNXj+6/WnPTioqkxcsoO2jerSuUlIsdvW9Pfl9cGdyc7NJyGmHsN6xp61Pg9c2hJfH+HvTg9o3oZM/jxxOe0ig/nw7kRqV0AQOZ+4rTUi4gv8E+gLpAPLRGSGqhZ+guge4FdVjRORwcDfgJuAfcBVqpohIvHAt0Cks82TQKaqthQRH6D4m6qNMV6ltJxSdQNrcP8fWvDMzHUs2LyPni1dz578vOMg63cf4tlrzn7bbHxkMF+N6EGj4MBSj9WgbiBDezTnzbmptGtcl1e+20TrhnX56O5E6gR63qyk7uyRJAKpqrpVVXOAScCgImUGAR86rz8D/iAioqopqlqQn2AtECgiBRcT7waeB1DVfFXd58Y2GGM8yJBuMUSH1uL5r9aT58zIOHFJGrUD/BjUKbKUraFlgzplDgR/uqQ5oUH+vPD1BuIiajP+nkSPndranYEkEthZaDmdU72KM8qoai6QBYQVKXMdkKKq2SJS0O8cKyI/i8gUEWlQ3MFF5E8ikiwiyXv37i1vW4wxHsDfz4dHL2/Fhl8OMy1lF78ezWHWqt1c3blxhV9uqhNYgzGD2tG7VQQT/tjtjEF8T+LOQFJc369oRrOzlhGRdrgud93rrPIDooCFqnoBsAh4qbiDq+rbqpqgqgkRERWTPsEYU/1d2b4RHaKCeXn2RiYuSSMnN59bu8e45VgDOjTm/bsSCS3lQcjqzp2BJB1oUmg5CiiaTvNkGRHxA4KBA85yFDAVuF1Vtzjl9wO/OesBpgAXuKPyxhjP5OMjjOzfmoys47zy3SYSYurRumHd0jc0JXJnIFkGtBCRZiLiDwwGZhQpMwO4w3l9PTBXVdW5hPUl8ISqLiworK5bLWYCvZxVfwAqPv2nMcajJcWG07tVBPlKibf8mrJz211bqporIsNx3XHlC4xT1bUiMgZIVtUZwHvAeBFJxdUTGexsPhyIA0aJyChn3WWqmgk87mzzGrAXuMtdbTDGeK5nBsbTfNF2rnByWZnfT6p6IpbKkJCQoMnJyVVdDWOMqVZEZLmqJpRWzp5sN8YYUy4WSIwxxpSLBRJjjDHlYoHEGGNMuVggMcYYUy4WSIwxxpSLBRJjjDHlYoHEGGNMuXjFA4kishdIK6VYOK55ULyNtdu7WLu9S3nbHaOqpWa99YpAUhYiklyWJzg9jbXbu1i7vUtltdsubRljjCkXCyTGGGPKxQLJKW9XdQWqiLXbu1i7vUultNvGSIwxxpSL9UiMMcaUiwUSY4wx5eL1gURE+onIRhFJFZGRVV0fdxKRcSKSKSJrCq0LFZHvRGSz87NeVdaxoolIExGZJyLrRWStiIxw1nt0uwFEJFBElorISqftzzjrm4nIEqftk52psD2KiPiKSIqIzHKWPb7NACKyXURWi8gKEUl21rn9s+7VgUREfIF/Av2BtsDNItK2amvlVh8A/YqsGwl8r6otgO+dZU+SCzysqm2A7sB9zu/Y09sNkA30UdWOQCegn4h0B/4GvOq0/Vfgniqso7uMANYXWvaGNhforaqdCj0/4vbPulcHEiARSFXVraqaA0wCBlVxndxGVRcAB4qsHgR86Lz+ELi6UivlZqq6W1V/dl4fxvXHJRIPbzeAuhxxFms4/xToA3zmrPe4totIFHAl8K6zLHh4m0vh9s+6tweSSGBnoeV0Z503aaCqu8H1RxeoX8X1cRsRaQp0BpbgJe12LvGsADKB74AtwEFVzXWKeOJn/jXgMSDfWQ7D89tcQIHZIrJcRP7krHP7Z92vondYzUgx6+x+aA8kIrWBz4EHVPWQ60uq51PVPKCTiIQAU4E2xRWr3Fq5j4gMADJVdbmI9CpYXUxRj2lzERepaoaI1Ae+E5ENlXFQb++RpANNCi1HARlVVJeqskdEGgE4PzOruD4VTkRq4AoiE1X1C2e1x7e7MFU9CMzHNU4UIiIFXyI97TN/ETBQRLbjulTdB1cPxZPbfJKqZjg/M3F9cUikEj7r3h5IlgEtnDs6/IHBwIwqrlNlmwHc4by+A5hehXWpcM718feA9ar6SqG3PLrdACIS4fREEJGawKW4xojmAdc7xTyq7ar6hKpGqWpTXP+f56rqEDy4zQVEJEhE6hS8Bi4D1lAJn3Wvf7JdRK7A9Y3FFxinqs9WcZXcRkQ+AXrhSi29BxgNTAM+BaKBHcANqlp0QL7aEpGLgR+A1Zy6Zv6/uMZJPLbdACLSAdfgqi+uL42fquoYEWmO69t6KJAC3Kqq2VVXU/dwLm09oqoDvKHNThunOot+wMeq+qyIhOHmz7rXBxJjjDHl4+2XtowxxpSTBRJjjDHlYoHEGGNMuVggMcYYUy4WSIwxxpSLBRLjlURkvogklF6y3Me538k8PLHI+k7OrecVdZw8J+PrGhGZIiK1nPUNRWSSiGwRkXUi8pWItKyo45ZSp0o5x6bqWSAx5hwVekK6LP4MXOE8FFdYJ6DCAglwzMn4Gg/kAMOchzGnAvNVNVZV2+J6hqZBBR7XGAsk5vwlIk2db/PvOPNpzHae0D7t266IhDspMRCRO0VkmojMFJFtIjJcRB5y5qZYLCKhhQ5xq4j85HyLT3S2DxLXvC3LnG0GFdrvFBGZCcwupq4POftZIyIPOOv+DTQHZojIg4XK+gNjgJucXsRNzpwR00RklVPPDk7Zp0VkvIjMdeaTGFqGU/cDEAf0Bk6o6r8L3lDVFar6QzHnufAcNY+IyNPO6/udnswqEZlUyjmq6fR+VonIZKBmGepqPIC3J200578WwM2qOlREPgWuAyaUsk08riy/gUAq8LiqdhaRV4HbcWUyAAhS1SQRuQQY52z3JK60Gnc76UWWisgcp/yFQIeiTwWLSBfgLqAbrgSBS0Tkv6o6TET64ZofYl9BeVXNEZG/AgmqOtzZx5tAiqpeLSJ9gI9w9VoAOuDKkRUEpIjIlwU5lYpyekv9gW+c9iwv5VyVZiTQTFWzC9KtUPI5uhf4TVU7OIHw53Ie21QT1iMx57ttqrrCeb0caFqGbeap6mFV3QtkATOd9auLbP8JnJynpa7zR/EyYKS4Uq/PxxWMop3y35WQWuJiYKqqHnXm//gC6FG25p22j/FOfeYCYSIS7Lw3XVWPOcFoHq5EfEXVdOqcjCsNxnvnePySrAImisituCYJg5LP0SU4QV5VVznbGi9gPRJzviucDymPU5dLcjn1RSjwLNvkF1rO5/TPfNH8QIqrR3Gdqm4s/IaIdAOOllDHishJf7ZU58XVs6hjqtqp8AoRWcupRIVnU/hcwunn80pcAWIgMEpE2lHyOSqpbsbDWY/EVFfbgS7O67L8sSzOTXAysWOWqmYB3wJ/cQaqEZHOZdjPAuBqEaklrqyr1+Aapzibw0CdIvsY4hyzF7BPVQ857w0S1/zrYbiSbi4rQ50A5gIBhcdVRKSriPQsUm4PUF9EwkQkABjglPUBmqjqPFwTRYUAtSn5HBVuQzyuS3LGC1ggMdXVS8D/iMhPuLIZ/x6/Otv/m1NzeI/FNSXtKmcAemxpO3Gm8v0AWIorq/C7qppSymbzgLYFg+3A00CCiKwCXuBU2m+c/X4JyUzGYAAAAIVJREFULAbGljQ+Uky9FFdQ6+vc/rvWOU5GkXIncA3+LwFmAQWTIfkCE0RkNa6Mua8685qUdI7+BdR22vCYU2/jBSz7rzHnMefuqSOq+lJV18WYkliPxBhjTLlYj8QYY0y5WI/EGGNMuVggMcYYUy4WSIwxxpSLBRJjjDHlYoHEGGNMufw/lU8wGCYVl7AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: 'argmin' is deprecated, use 'idxmin' instead. The behavior of 'argmin'\n",
      "will be corrected to return the positional minimum in the future.\n",
      "Use 'series.values.argmin' to get the position of the minimum now.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "formula           y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...\n",
       "mean_test_acc                                              0.972126\n",
       "mean_test_mse                                             0.0278736\n",
       "mean_train_acc                                             0.972824\n",
       "mean_train_mse                                            0.0271759\n",
       "num_pc                                                           10\n",
       "Name: 9, dtype: object"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(df_res.num_pc, df_res.mean_train_mse, label='training set')\n",
    "ax.plot(df_res.num_pc, df_res.mean_test_mse, label='test set')\n",
    "ax.legend()\n",
    "#ax.set_xticks(df_res.num_pc.values)\n",
    "ax.set_xlabel('number of top PC used')\n",
    "ax.set_ylabel('Mean squared error')\n",
    "plt.title(\"MSE for Numbers of Top PC Used\")\n",
    "plt.show()\n",
    "best_model = df_res.loc[df_res.mean_test_mse.argmin()]\n",
    "best_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXd4VFXawH9vCgmBEEpCDZCE3ruA0oKKICoqqOhasKGurq5+2Na16+qqa6+o2MWOi64gIF2aQAApAUKAEEJJIY305Hx/nDvJZDKTmZRJApzf88yTmXvPPffcO5P7nrec9xWlFAaDwWAwVBef+h6AwWAwGE5tjCAxGAwGQ40wgsRgMBgMNcIIEoPBYDDUCCNIDAaDwVAjjCAxGAwGQ40wgsRwSiAi54jIXhHJFpFL63s8VUVEZojI6voehyMi0kNEYkQkS0Turu/xnGqIyFci8s/6Hkd9YwTJaYCILBeREyISUN9j8SJPAW8qpZoqpX6saWci8rGIKBE5y25bVxE50xZWPQAsV0oFK6Vet98hIjsswZ0tIsUikmf3+R+1OQjrgZxv9Z0mIgtFpJvd/t4i8oOIpIpIuohsEZG7RUSc9HW7iCxxsv2oiIyqzXEbNEaQnOKISAQwGlDAJXV8br86PF1nYEd1DqxknGnAM9UeUT1Ry/fd5X1VSvWxBHdTYBVwl+2zUupftTgGG09b5+oEZAHvA4hIT2AtsBvoo5RqDlyD/t2fzpOnUwYjSE59rgfWAR8DN9jvEJHGIvIfETkoIhkislpEGlv7RonIGmt2d0hEZljbl4vILXZ9lDPJWLP4O0VkL7DX2vaa1UemiGwSkdF27X1F5B8iss8yn2wSkY4i8paI/MdhvD+JyN8dL1BE9gFRwE/WjDVARNqLyHxr9honIrfatX9CRL4Tkc9FJBOY4eLefQL0F5GxznaKyAEROc+h38+t9xHWvbjRuvYT1kx4mIhss+7rmxW7lDes7yJWRM612xEiIh+KyBEROSwiz4iIr9138LuIvCIiacATlva0wuorRUS+dnGNiMgllnaRbn2/vaztS4Fo4E3rvnZ31YeLfn1F5EkRSRCRYyIyR0SCrX09RaTIuidHRCRJRP7mSb9KqWzgK6CvtekZYIlS6mGl1FGrzU6l1BVKqbyqjNlu7FOs7yDL+v7uttt3md13uEpEetvtO0tEtlrHfQ40qs75TzeMIDn1uR74wnpdICJt7Pa9BAwBzgZaos0YJSLSCVgAvAGEAQOBLVU456XAcMD2D/aH1UdL4EvgWxEJtPbdB1wNXAg0A24CctAP8atFxAdAREKBc4G5jidTSnUBEoCLrdlwvtUuEWgPTAP+Zf9gBqYA3wHNrXvjjBzgX8CzVbh2R4YD3YCrgFeBR4DzgD7AlQ5CajgQD4QCjwM/iEhLa98nQBHQFRgETABucXJsa2u8TwOLgBZAOPq7rIAlHOYCf0d/17+gBXIjpdR4ymsae6p47bcBV6I1g27W2F622+8LjERPAiYDT4oHpiURaYb+zcRYm85Df5e1yRzgeqVUMPq3u8o69wjgbeBGoBXwGfCjiPhZv+kfgffQv/UF1LEVoKFiBMkpjPVP2Rn4Rim1CdiHVvmxHtA3AfcopQ4rpYqVUmush/Bf0DO8uUqpQqVUqlKqKoLkOaVUmlIqF0Ap9bnVR5FS6j9oc0MPq+0twD+VUruVZqvVdgOQgRYeANPRtvpjHlx3R2AU8KBSKs8a+wfAdXbN1iqlflRKldjG6YL3gE4iMqkK12/P09YYFgEngblKqeNKqcPoh9Mgu7bHgVete/412lQz2RL+k4C/K6VOKqWOA6+g74mNJKXUG9Y9zgUK0d99e+v8rhz5VwH/U0otVkoVoicXjdGTi5ryF+BFpdRBpVQmWoj+RaSc3+JxpVSuUioG+BwtIFzxiIiko++LL3CLpZWFAEdqYbz2FAF9RCTY+j3ahNZtaF/cJut/Zjb69zwEGAPkK6Xetr7DL4BttTyuUxIjSE5tbgAWKaVSrM9fUmbeCgUC0cLFkY4utnvKIfsPIvJ/IrLLMrOko//xQz041yfAtdb7a9GzP09oD6QppbLsth0EOrgaoysswfq09arguPUAe8GX6+RzU7vPh1X5LKkH0dfSGfAHjljmlHS0gGtt19bxeh6wxrvBMlvd5GJ87a3zAKCUKrH66uCifVUo17f1vjF6tm7jkMP+9pX096xSqrlSqp1S6jJLQBWjJxztqjCuIvT9dMQfLYBBa9VTgQQRWSoiQ63tnYF/2L4H67sIQ9+v9mgt2J6DGIwgOVUR7eu4EhgrOhrlKHAvMEBEBgApQB7Qxcnhh1xsBz2rDrL73NZJm9KHoWh/yIPWWFpYjtAMyh7KlZ3rc2CKNd5eaLOBJyQBLW32eItOwGFnY/SAj9DC7zKH7Z7ci6rQwWG23gl9LYeAfCDUepA2V0o1U0r1sWtb7nqUUkeVUrcqpdqjZ9Fvi0hXJ+dMQj8cAe2kQQv3w07aVpVyfVvXk4sOYrDR0WF/UjXOswT90PeUBIdx2cxlLbEe/EqptUqpi4A2aBOhzaR6CHjM7ntorpQKUkr9gNaKwh3O1anKV3MaYgTJqculQDHaTzHQevVCm1Out2aec4CXRTumfUVkpOgQ4S+A80TkSsv220pEBlr9bgEuF5Eg68F0s5txBKNngMmAn4g8hvaF2PgAeFpEuommv4i0AlBKJaL9K58B37sxQZWilDoErAGeE5FAEelvjdOVL8Rdf0XAE2iBaM8WYLqI+Fsz1mnV6d+O1sDdVn9XoL+vX5RSR9APs/+ISDMR8RGRLuIiCABARK4QEdtD7QRa0BQ7afoN2nx2roj4A/+HFlprangtoB++s0SkkyXUnwG+dNC6Hhcd9DEAbXp0GRRQCY8C54vIszYfoOj1L1/b+eLsWY3+Ld4nOjAjGHgBWKGUOioiTURkuiVcCtERYrZ7Nxv4m4gMtX6vTUUHKwQBK4FAK4DAT0SuBvpX43pOO4wgOXW5AfhIKZVgzU6PWhEtb6Lt1H7ALOBP9MM6Dfg34KOUSkA7v//P2r4FGGD1+wpQgDbRfIL7h/OvaKfjHvRsL4/y5oyX0Q+zRUAm8CHa/GHjE6Afnpu1bFwNRKBnuPPQtvjFVezDnrlUtMM/itamTgBPok2HNWE92imdgnaYT1NKpVr7rkdHAO20zvcdlZtzhgHrRSQbmI/2he13bKSU2o02G75hnfdidNBCQQ2vBeAd4Ae0UNqH/i3dZ7e/GH3N+4GFwFNKqZVVPYlSahfap9Mb2GWZm75CT5rynbTPQf++L0T/Pvaigy6usWt2E/r3moG+9zdYx/4O3I02Laajf9fX6F0qF621/hX9HU0Gfqrq9ZyOiDKFrQz1iIiMQZu4IiwtynAaIHrtx3alVF2uNTLUE0YjMdQblqnlHuADI0QMhlMXI0gM9YLoRXHpaPPNq/U8HIPBUAOMactgMBgMNcJoJAaDwWCoEWeEIyw0NFRFRETU9zAMBoPhlGLTpk0pSqkwd+3OCEESERHBxo0b63sYBoPBcEohIh6t3DemLYPBYDDUCCNIDAaDwVAjjCAxGAwGQ404I3wkzigsLCQxMZG8vGrVxTFUkcDAQMLDw/H3d5aU1WAwnMqcsYIkMTGR4OBgIiIikIplnw21iFKK1NRUEhMTiYyMrO/hGAyGWuaMNW3l5eXRqlUrI0TqABGhVatWRvszGE5TzlhBAhghUoeYe20wnL6c0YLEYDAYTnUSUnNYvNNthWqvYgRJPZGens7bb79drWMvvPBC0tPTK23z2GOPsWTJkmr1XxN+/PFHdu7cWefnNRjOVF5ctJvbP99Edn5RvY3BCJJ6ojJBUlzsrNBdGb/88gvNmzevtM1TTz3FeeedV+3xVRcjSAyGuqOkRPF7XArFJYpNB0/U2ziMIKknHnroIfbt28fAgQO5//77Wb58OdHR0VxzzTX069cPgEsvvZQhQ4bQp08fZs+eXXpsREQEKSkpHDhwgF69enHrrbfSp08fJkyYQG6urlY7Y8YMvvvuu9L2jz/+OIMHD6Zfv37ExsYCkJyczPnnn8/gwYO57bbb6Ny5MykpKeXGWVxczIwZM+jbty/9+vXjlVdeAWDfvn1MnDiRIUOGMHr0aGJjY1mzZg3z58/n/vvvZ+DAgezbt8/r99FgOJPZeSSTtJO62OW6+FQ3rb3HGRv+a8+TP+1gZ1JmrfbZu30zHr+4j8v9zz//PNu3b2fLli0ALF++nA0bNrB9+/bSENk5c+bQsmVLcnNzGTZsGFOnTqVVq1bl+tm7dy9z587l/fff58orr+T777/n2muvrXC+0NBQNm/ezNtvv81LL73EBx98wJNPPsn48eN5+OGHWbhwYTlhZWPLli0cPnyY7du3A5Sa1GbOnMm7775Lt27dWL9+PX/9619ZunQpl1xyCRdddBHTptW0vLnBYHDH6jg98YsMbWIEiUFz1llnlVtn8frrrzNv3jwADh06xN69eysIksjISAYOHAjAkCFDOHDggNO+L7/88tI2P/zwAwCrV68u7X/ixIm0aNGiwnFRUVHEx8fzt7/9jcmTJzNhwgSys7NZs2YNV1xxRWm7/PwKpbMNBoOXWbU3mR5tgjmvd2veXRHPyfwimgTU/WPdCBKoVHOoS5o0aVL6fvny5SxZsoS1a9cSFBTEuHHjnK7DCAgIKH3v6+tbatpy1c7X15eiIu2U86SoWYsWLdi6dSu//vorb731Ft988w2vvvoqzZs3L9WmDAZD3ZNXWMwfB05w3YjOjIhqxVvL9rHp4AnGdHeb9b3WMT6SeiI4OJisrCyX+zMyMmjRogVBQUHExsaybt26Wh/DqFGj+OabbwBYtGgRJ05UdNalpKRQUlLC1KlTefrpp9m8eTPNmjUjMjKSb7/9FtACaevWrR5dl8FgqB027E+joKiEUd1CGdK5BX4+Um/mLSNI6olWrVpxzjnn0LdvX+6///4K+ydOnEhRURH9+/fn0UcfZcSIEbU+hscff5xFixYxePBgFixYQLt27QgODi7X5vDhw4wbN46BAwcyY8YMnnvuOQC++OILPvzwQwYMGECfPn3473//C8D06dN58cUXGTRokHG2GwxeZHVcCo18fRge2ZKgRn70Dw+pN0FyRtRsHzp0qHIsbLVr1y569epVTyNqGOTn5+Pr64ufnx9r167ljjvu8Kq5ytxzg6H2mPTaKkIa+/HVzJEAvLAwltkr49n2xASCGtWO10JENimlhrpr51WNREQmishuEYkTkYec7O8sIr+JyDYRWS4i4db2aBHZYvfKE5FLrX0iIs+KyB4R2SUid3vzGk5nEhISGDZsGAMGDODuu+/m/fffr+8hGQwGD0jOymfXkUxGdyvzh4yIakVRPa0n8ZqzXUR8gbeA84FE4A8Rma+Usl+t9hLwqVLqExEZDzwHXKeUWgYMtPppCcQBi6xjZgAdgZ5KqRIRae2tazjd6datGzExMfU9DIPBUEXW7NNhv6O6hpZuG9K5Bb6Wn8RewNQF3tRIzgLilFLxSqkC4CtgikOb3sBv1vtlTvYDTAMWKKVyrM93AE8ppUoAlFLHa33kBoPB0IBZtTeFkMb+9O0QUrqtSYDNT5JW5+PxpiDpAByy+5xobbNnKzDVen8ZECwirRzaTAfm2n3uAlwlIhtFZIGIdHN2chGZabXZmJycXO2LMBgMhoaEUorVe1M4p2srfH3KZ9UeEdWKrYfSySmo27xb3hQkzvKGO3r2ZwFjRSQGGAscBkrvgIi0A/oBv9odEwDkWQ6g94E5zk6ulJqtlBqqlBoaFlb3cdUGg8HgDfYlZ3M0M49RXSs+12x+ks0HK0/qWtt4U5Akon0ZNsKBJPsGSqkkpdTlSqlBwCPWtgy7JlcC85RShQ79fm+9nwf0r+2BGwwGgzNyC4rrfLbvyKq92j8yultohX32fpK6xJuC5A+gm4hEikgjtIlqvn0DEQkVEdsYHqaidnE15c1aAD8C4633Y4E9tTrqOqImaeQBXn31VXJyctw3dMPy5ctZs2ZNjfsxGM4Eps9ey8CnFnPjRxuYuyGB41l1X/Vz9d4UIloF0bFlUIV9TQP86Neh7teTeE2QKKWKgLvQZqldwDdKqR0i8pSIXGI1GwfsFpE9QBvgWdvxIhKB1mhWOHT9PDBVRP5ER3nd4q1r8CZGkBgMpxaH0nLYmpjBwPDm7D2ezcM//Mnwf/3GZW//ztvL40hKd56eqDYpLC5hXXwqo5xoIzZGRLVia2I6uQWVl6OoTby6jkQp9YtSqrtSqotS6llr22NKqfnW+++UUt2sNrcopfLtjj2glOpgi86y256ulJqslOqnlBqplNrqzWvwFo5p5AFefPFFhg0bRv/+/Xn88ccBOHnyJJMnT2bAgAH07duXr7/+mtdff52kpCSio6OJjo522nfv3r3p378/s2bNAnTK+KlTpzJs2DCGDRvG77//zoEDB3j33Xd55ZVXGDhwIKtWraq7G2AwnGIs260DRP89rT+rHohmwT2jufe87hQVK15YuJtp76yhuMS7C7xjEtI5WVDs1D9iY3hUSwqLFZsT6m49iUnaCLDgITj6Z+322bYfTHre5W7HNPKLFi1i7969bNiwAaUUl1xyCStXriQ5OZn27dvzv//9D9A5uEJCQnj55ZdZtmwZoaHlZyZpaWnMmzeP2NhYRKQ07fs999zDvffey6hRo0hISOCCCy5g165d3H777TRt2rRU4BgMBucsjT1OZGgTIkN1ctVe7ZrRq10z7j63Gz/GHObvX2/hjwNpjIhyDDytPVbtTcZHYGQX1+cYaucnOaera82lNjGCpIGwaNEiFi1axKBBgwDIzs5m7969jB49mlmzZvHggw9y0UUXMXr06Er7adasGYGBgdxyyy1MnjyZiy66CIAlS5aUq1yYmZlpkisaDB6SU1DEmn2pXDu8s9P95/duQ4CfDwu3H/WyIElhQMfmhDT2d9kmOFCvL6lLP4kRJFCp5lBXKKV4+OGHue222yrs27RpE7/88gsPP/wwEyZM4LHHHnPZj5+fHxs2bOC3337jq6++4s0332Tp0qWUlJSwdu1aGjdu7M3LMBgaHBm5hfyZmFGpX8Eda+JSKSgq4dxezhNpNAnwY0z3MH7dcZTHLuqNj4+z1Q+esfFAGn6+PvTvEFKun4ycQrYlpnNXdFe3fYyIasmc1fvJLSimcSPfao/FU0z233rCMd36BRdcwJw5c8jOzgZ01t3jx4+TlJREUFAQ1157LbNmzWLz5s1Oj7eRnZ1NRkYGF154Ia+++mqp6WzChAm8+eabpe1s203ad8Ppzvsr47n2w/V8s/GQ+8YuWLr7OE0D/BgW0dJlm4l92nIkI4+tidVbw6GU4rUle5n27loufet3Rj7/G4/M+5Plu4+TX1TM2vgUShSM8iD9yYjIVhQWK2LqyE9iNJJ6wj6N/KRJk3jxxRfZtWsXI0fqTJ5Nmzbl888/Jy4ujvvvvx8fHx/8/f155513AF3qdtKkSbRr145ly5aV9puVlcWUKVPIy8tDKVVaY/3111/nzjvvpH///hQVFTFmzBjeffddLr74YqZNm8Z///tf3njjDbemM4PhVGOtZeL557ztdGvdlEGdKlYCrQylFMtijzOqayiN/FzPvc/r1QY/H2Hh9qNVPkdJieKpn3fy8ZoDTB0czjldW7F45zHmxRzmi/UJNA3wo3mQP00a+TKoU3O3/Q2NaIGP6DruZ9eBn8SkkTfUGeaeG+qanIIi+j+xiGlDwvl9XwoFRSX8dNcoWjcL9LiPnUmZXPj6Kl6Y2p8rh3WstO31czZwIOUkK+4fh4hn5q3C4hIe+G4b82IOc/OoSB65sFepSSuvsJi1+1JZtPMYS2OPMa57a/49zbM12FPeXE2Any/f3D7So/bOaBBp5A0Gg6E+2XTwBEUliol92zL7uqFk5hZx2+ebyC/yfI2FLex3XE/3JqVJfduSkJbDziOZHvWdV1jMHZ9vYl7MYWZN6M4/J/cq5xcJ9Pclumdrnru8H+v/cZ7HQgRgeFQrthxKJ6/Q++tJjCAxGAynLevj0/D1EYZGtKRXu2a8fOUAYhLSeezHHXhqjflt1zH6h4fQOti9FjOhdxt8BBZuP+q2bWZeIdfP2cBvscd5+tK+3DW+m8dajCeMiGpJoL8P+1NO1lqfrjijBcmZYNZrKJh7bagP1sWn0q9DCE0DtDt4Ur92/G18V77eeIjP1h10e3zayQJiDqUT3cOzsketmgZwVmRLFrgRJFl5hVzz/jo2HzzBq1cN5LoRzsOKa8KYbmHEPDaBXu2a1XrfjpyxgiQwMJDU1FTzgKsDlFKkpqYSGOi5XdpgqCk5BUVsTUxneFT5SKt7z+vOeb1a89RPO92utVix5zhK4TLs1xmT+rYj7ng2ccddR0M+vyCWHUmZzL5+CFMGOlbXqB38fH0qpJn3Fmds1FZ4eDiJiYmYWiV1Q2BgIOHh4fU9DMMZxOaD6RQWqwoLBH18hFeuGsilb/3OX7/YzPy7ziG8RcUEiAC/7TpOaNMA+rYPcbrfGRf0acvj83ewcPtR7hofXGH/+vhUvlifwM2jIhnfs03VLqqBcsYKEn9/fyIjI+t7GAaDwUus35+q/SOdK4biBgf68/71Q5ny1u/M/HQT399xdoWFe0XFJazck8wFfdpWaYFh25BABndqzoLtR7lrfPm6e3mFxTz0w590bNmY/5vQvXoX1gA5Y01bBoPh9GZdfCp9O4QQHOg8nUhUWFNenz6IXUczuf+7rRXM3JsOniAzr6hKZi0bE/u2ZUdSJgmp5TN0v/bbXvannOS5y/oT1Oj0mccbQWIwGE47cguK2XIonRGRrleiA0T3bM39F/Tg521HeHdFfLl9S2OP4+8rHq0kd2RS33YALNxxpHTb9sMZzF4Zz5VDw2uUrqUhcvqIRIPBcMqxIymD+77eSnZ+xaqDAX4+PHtZv0oz3bpic8IJp/4RZ9wxtgs7kzJ54ddYerYLLo3QWhp7nLMiW5ZGfFWFji2D6NO+GQu2H2XmmC6liw5bNmnEIxf2rnJ/DR2vaiQiMlFEdotInIg85GR/ZxH5TUS2ichyEQm3tkeLyBa7V56IXOpw7Bsiku3N8RsMBu+hlOLJn3ZyPCuPkV1aVXjlFRbz6H+3U1Rc4r4zB9bHp+IjOlWIO0SEF6b1p1fbZtw9N4b45GwOpeWw93h2jZzhk/q2JSYhnSMZuby/Kp6dRzJ5ekpfQoJcZ+49VfGaRiIivsBbwPnoOut/iMh8pdROu2YvAZ8qpT4RkfHoiofXKaWWAQOtfloCccAiu76HAu4TzhgMhgbL8t3JbNifxlNT+nD9yIgK+xduP8rtn2/i+82JXDWsU5X6XhefRr9K/COOBDXyY/b1Q7jkzd+59dONXD5YRxiO71l1/4iNiX3b8dKiPby7fB9z/zjEpL5tmdi3bbX7a8h4UyM5C4hTSsUrpQqAr4ApDm16A79Z75c52Q8wDViglMqBUgH1IvCAV0ZtMBi8TnGJ4t8LY+ncKojpLoTEBX3aMKhTc15ZvLdKZWNt/pHhVawLEt4iiLeuGcyB1BxeWrS7XBGr6tC1dVO6tW7KJ2sPEujnw5NT+lS7r4aONwVJB8A+b3Oitc2ercBU6/1lQLCIOH7704G5dp/vAuYrpY5gMBhOSX6MOUzs0SxmTejhMqOuiPDQxJ4czczj4zUHPO47JuEEBcUljIiq3NHujJFdWvHYRb1RCo9Xs1fGJEsDefSi3h6lWDlV8aaz3VngteMy8lnAmyIyA1gJHAZKvW4i0g7oB/xqfW4PXAGMc3tykZnATIBOnaqmFhsMBu+RV1jMy4v30K9DCJP7tau07fCoVpzbszVvL4/j6rM60jyokdv+1+1Ps/wjVRckANeP7Ezr4ACGuYn48oSbR0fRrU0wF/Wv/DpPdbypkSQC9jmXw4Ek+wZKqSSl1OVKqUHAI9a2DLsmVwLzlFKF1udBQFcgTkQOAEEiEufs5Eqp2UqpoUqpoWFhVQ/fMxgM3uGztQc5nJ7Lw5N6erTQ74GJPcnOL+Lt5fs86n9dfCp92ofQzEP/iCMiwqR+7QhtGlCt4+0JaezPxQPa12oyxoaINwXJH0A3EYkUkUZoE9V8+wYiEioitjE8DMxx6ONq7MxaSqn/KaXaKqUilFIRQI5Syn3dSYPB0CDIyC3kzWVxjOke5nHBpR5tg5k6OJyP1xzgcHpupW3zCovZkpBeLbOWofp4TZAopYrQ/oxfgV3AN0qpHSLylIhcYjUbB+wWkT1AG+BZ2/EiEoHWaFZ4a4wGg6FueXfFPjJyC3lwYo8qHXfv+TqdyMuL9lTaLiYh3fKPVH3tiaH6eHVBolLqF+AXh22P2b3/DvjOxbEHqOicd2zTtOajNBgMdcHRjDzmrN7PpQPb06cKSRABOjRvzIyzI3h/VTy3jomkZ1vnqdHXla4fMRpJXWJSpBgMhjrh1SV7KFGK/5tQNW3Exl/HdaFpgB8vLNztss26+FR6t29GSOPTb9FfQ8akSDEYDE5RSrF2XypDI1q6DNF1ZH/KSTYfPFFhe05BEd9sPMQNZ0fQsaXzlO3uaB7UiL+O68q/F8ay4M8jFbLy5hUWE3Moneu9UCTKUDlGkBgMBqes2JPMjI/+YNaE7hXSoTujuERx88d/EO+itGto0wDuiq5ZbMyN50Tw5YaD3PHFZto0C+C8Xm04v3cbRnbR9ckLiox/pD4wgsRgMDjl3RU63PbjNQe4ZXQUgf6+lbZfvPMo8Skn+ddl/RjlJCKrVdNGNKlGAkR7Av19+flvo/lt1zEW7zzGvJjDfLE+gaYBfrRpFoAItbL+w1A1jCAxGAwV2HIonXXxaUzo3YZFO4/x/eZE/jLctclIKcU7K+Lp3CqIq4Z19GqJ15DG/lw+OJzLB4eTV1jMmn0pLN55jMU7j3NOl1DjH6kHjLPdYDhD2JeczcLtnmUWem/FPoID/Xj5qoEMCA/h/ZXxFJc4JqYoY/3+NLYeSufW0VF1VicctIYyvmcbnru8P388ci6f3XxWnZ3bUIYRJAbDGUBuQTE3fvQHt3++mTX7Uiptuz+U0+ySAAAgAElEQVTlJAt3HOW6EZ1pGuDH7WO7cCA1h193HHV5zHsr9tGqSSOmDQmv7aF7jIic9ivIGypGkBgMZwAvL95NQloOYcEBPPzDn5Vm031/VTz+vj7MOCcCgAl92hLRKoj3VuyrUI4WIPZoJst2JzPj7Ai3fhTD6YkRJAbDac7WQ+l8uHo/1wzvxOvTB3EwNYdXlzhfIX48K4/vNiUydXB4abZaXx/h1jFRbE3MYG18aoVjZq+IJ6iRL9eNNGG3ZypGkJzJbP4Mdv3sefv0BFj4Dygq8N6YDLVKQVEJD36/jdbBgTw0qScju7Ti6rM68f6qeLYlpldo/8maAxQWl3Dr6Mhy26cODie0aSPec6hrfjg9l/lbk5g+rJNHmXkNpydGkJypFBfBr/+A31/1/Jjt38O6t+Dg794bl6FWeXfFPmKPZvHMpX1Ls+E+fGFPwoIDeOC7bRQUlZWxzc4v4rO1B7mgd1uiwspnHwr09+XGcyJZsSeZnUmZpds/XLUfgJsdBI/hzMIIkjOVpM2QnwnJe8CJ3dspyZY5JH6ZR80Likp4+Idt/LA5sZqDNNSEvceyeGPpXi4e0J7zepfVHm8W6M8zl/Yj9mgW760oS83+1YYEMvOKuG1slNP+rh3emSaNfJm9Uh+TnlPAV38kcMmA9nRo3ti7F2No0BhBcqayzxIG+RmQ5ToapxzJseWPdcMTP+1g7oZD3PfN1tLFbYa6obhE8eD322ga4MfjF/eusP/83m24qH873lgaR9zxLAqKSvhw9X6GR7ZkUKcWTvsMCfLn6rM68dO2IySeyOHzdQfJKShmpgvBYzhzMILkTCV+GfhaNu2UsiR4mXmFPDF/B3fPjSm/bkApSNmjjzm6DU5WHkL6+bqDfLk+gZljorh4QHueXxDL8wtinUb9NBS+WH+Q6bPXsvFAWn0PpcZ8uvYAmxPSefziPi4LND1xSR+CAnx54Ltt/LjlMEcy8rh9XJdK+71pVCQCvLUsjo9+P0B0jzCXmXgNZw5mZfuZSH4WJP4B/a+CLV9A8m5U5Fh+2naEp3/eSXJWPgD9w0O4ZbQ128w8DAXZMPAv+pj45dBvmtPuN+xP44n5OxjXI4wHJ/YEIKSxn1WLooBnLu1Xp4vW3KGU4o2lcby8eA8Bfj5Me3ct04d15MGJPWnR5NRzIB9Ky+GFhbsZ1yOMKQPbu2wX2jSAxy7qzX3fbGXnkUx6tAlmXPfKq4m2b96YSwa2Z+6GQwDcNrZywWM4MzAayZnIgdVQUgQDpkNACJmHdnD9nA3cPTeGdiGB/HTXKMb3bM1/Fu0hITVHH5NsaS39r4TA5i79JEnpufz1i010ahnEa9MH4esj+PoIT0/py9/Gd2XuhkPcPTeG/KLy6xgSUnP4YFU8V723loveWMXKPcnevAOllJQonvp5Jy8v3sPlgzuw4ZHzuG1MFN9uSuTcl1fw3abEBq1FOZKVV8jfv96Cj8Czl/Vzu0DvskEdGNs9jLzCEm4bG+XRgr7bxmjhMaBjc4abvFYGvKyRiMhE4DXAF/hAKfW8w/7O6PK6YUAacK1SKlFEooFX7Jr2BKYrpX4UkS+AoUAhsAG4za6mu8ET9i0D/yDy2w0ltVFnDv25gS1M5qkpffjL8M74+gjPXNqXCa+s5OF52/j85uFIiuVob90HIsfAvuXa3CXl03jP/Gwj+YUlzJ45tFzOIxHh/yb0IKSxP8/8bxeZeYX8/bxuLN+dzOKdx4g9mgVAjzbB5BcVc/2cDVw8oD2PTu5F62aBXrkNRcUlPPD9Nn7YfJgbz4ng0cm98fERHr6wF5cN7sAj87Yz69utfLvxEM9e1peurYO9Mg4bSqkarcxOzc5nxkd/sOtIJq9NH+SRA1xE+M+VA1iw/SiXDHCtvdjTo20wL0zrT9/2IWYluQEA8dZsS0R8gT3A+UAiuob71UqpnXZtvgV+Vkp9IiLjgRuVUtc59NMSiAPClVI5InIhsMDa/SWwUin1TmVjGTp0qNq4cWNtXdqpz5vDyG8azmUZ93FD8ktMCthK7j2xtHF4YH++7iD//HE7L0ztz5VHX4Kd/4UH9sOmj+Dne+GujRCq04srpbj36y38d2sSH1w/lHN7tXF2ZgC+25TIg99vo7hE4SMwLKIl5/fW6cA7t2pCXmEx762I563lcQT4+nD/xB6lAs5GclZ+aQbY9fvTeGFafy7s187jW5BXWMxdX8awZNcx7ju/O38b37XCQ7GkRPHNxkM8tyCW7PwiRkS15PxebTivdxvCW1SsqZFbUMyqvVowLtt9nJFdQnntqoHlama44uPf9/POin28cfVgzqrGLD8pPZdrP1zP4RO5vHPtYMb3dH3/DQZPEZFNSqmhbtt5UZCMBJ5QSl1gfX4YQCn1nF2bHcAFlhYiQIZSqplDPzOBsUqpvzg5x71AqFLqkcrGYgSJHRmH4ZXevNXoRt7Mnch/B22m+7Z/awERVP4BVlKimP7+OmKPZLKxw8s08gFu/hXS9sPrA2HSizB8JgCzV+7jX7/Eely7YnPCCRJScxjbPcylH2J/ykke/XE7q+NS6B8ewr3ndyf2SBaLdx4l5lA6SlFu1r101lgC/Nyn6MjKK+TWTzeyLj6Np6b04fqREZW2T8nOZ87q/SzaeYy449kA9G7XjPN7t2FsjzDijmWzaOcxVsclk1dYQnCgH73bNWP9/jTuPa8795xX+f1YtTeZG+ZswNdH8BGpsiDYl5zNdR+sJyu/iDkzhjHMlJk11BKeChJv+kg6AIfsPidSsQb7VmCq9f4yIFhEHKvSTAfmOnYuIv7AdcBCZycXkZkislFENiYn1429/VTg6BatzC0t6MMXtw6ne1/rN5JcsXypj4/w/OX9yCsqoeDITgizSqS2jIQWERC/jOz8Ip76aSfPLYjlwn5tudPDwkWDO7Xg0kEdKnVmR4Y24bObz+K16QNJSs/jxo/+4N8LYykoLuHv53bnl7tHs/rBaP49tT+H03P5fF2C2/MqpfjrF5vZeOAEr00f6FaIgHZKPzCxJ0vuG8uyWeP4x4U9CWrky+tL93L522t44Ptt7EzK4MqhHfns5rPY9M/z+WrmCC4b1IFXluxhUSXJDg+mnuSuL2Po1jqYpf83jh5tg5n56SZ+jDnsdlwA2w9ncMW7aykoLuGrmSOMEDHUC970kTjT5x3Vn1nAmyIyA1gJHAaKSjsQaQf0A3510tfbaLPWKmcnV0rNBmaD1kiqOvjTkZiEExxd9gPDaM6/bruKHu2awYnuemfKbug8ssIxUWFNeXhMGE3XZLKrqB29rO0qKprird9wwUu/kZRdxDVndeKRyb1q3WYuIkwZ2IFxPVrze1wKAzo2r2D7H9UtlNHdQnlz6V6uGBpeuoLbGd9uSmTV3hSemtKHKQMd5zXuiQxtwswxXZg5pgsp2fms3ZdKZGgT+rRvVuHan7u8H/uSs7n36y38eOc5dGtT3sdyMr+ImZ9uAmD29UPo2DKIL28dwa2fbOTvX28hI7eQG86OcDmWNXEpzPxsEyGN/fn8luFEhjap8vUYDLWBNzWSRKCj3edwIMm+gVIqSSl1uVJqEPCItS3DrsmVwDxHZ7qIPI520N/njYGfjqzam8y1H6xlBNtp3ONcLUQAQjqBX2OnGomN67rmAvDWdj8ycgo5lJbDWwc74ld0kpGBB/jhjrN59rJ+BDXy3rwkpLE/F/Zr59KB/ODEnpzIKSy3UtuR45l5PPPzTs6KaMm1lRRp8pTQpgFcPKA9fTs4dzoH+vvy3nVDaNzIj1s/3UhGTtnPuKRE8X/fbGXv8SzeumYwnVtpIdA0wI+PbhzG+b3b8Pj8Hby2ZG9p1FhJiWLroXRe+nU3F7yykms+WE+7kEC+v+NsI0QM9Yo3NZI/gG4iEonWNKYD19g3EJFQIE0pVQI8jI7gsudqa7v9MbcAFwDnWsedUWTk6gd5r3bNPFqLUVKi+GlbErO+3cr5LY7TIjsDep1f1sDHRzvMKxEkfqk6YmtLXhuun7Oe3ceyaC6d+auPDy8MSsXHxUrouqRvhxAuGdCeD1fv5/qRERUCBwAen7+DvKISnp/azyMHeG3QLqQx7147mKvfX8fdX8UwZ8YwfH2EN5fFsXDHUf45uRejupUvSxvo78s7fxnMQz/8yStL9pCUnoufr7Bk1zGOZeaXBij8c3Ivrhja0VQENNQ7XhMkSqkiEbkLbZbyBeYopXaIyFPARqXUfGAc8JyIKLRp607b8SISgdZoVjh0/S5wEFhrzQJ/UEo95a3raAgcTs9l8Y6jLN51jPXxaRSVKFo1acT4nq05v3cbRncLo3GjMiezrfzooh3HWLLrOCnZ+Qzt3IL/9IzVdzNqXPkThPWEg2tcDyBlDzRqypQRw3hreTwT+7TlsYt74/PtINi/HPhH7V90NZg1oQcLth/h1SV7ee7yfuX2LfjzCAu2H+WBiT0qJCT0NkMjWvLkJX35x7w/eeHXWIZ2bqnXrQzqwM2jnCc79PP14YWp/Qlp7M+Hq/cT1MiXsd3DOL93G6J7tD4lF0oaTl+8FrXVkKhu1NbOdQvJTfPM6VnbKAXHMvPYlphB4gltWmrTLID+4c1pGxLIrqRMdhzJILegBH8/oWebZnRu35r5WT1ZEZdGTkExTQP8GNsjjAm923BBn7YEzp2q82rdua78yVa+CEufgYcTIcDJWolPp0BeBiW3LONA6smyB/FvT8PqV+DB/RAY4uU74hlPzN/BZ+sO8uvfx9C1tR5nRk4h5768gjbNAvjxznPw962fdbiPzPuTL9YnEODnQ4+2wXxz20i3haCUUiSk5dCmWaApGmWoczyN2jIpUiqhcMV/GJK7oV7HMBnANvnMQ6+oAUbYGtj2JetXus9FtBr8EOf3bsuIqJZl4bCFeZCwFobeVPEkoVY0Vsoe6DCk4v7kPRA5Bh8fKT+b7xINq17SK+V7Tq7+RdYifxvfle82JfLSr7t59zp9Lc/8bycncgr4+MZh9SZEAB6/uA9xx7PZn3KS964b4pFgEJFS/4nB0FAxgqQS2lz9NgdPZrpv6CVCghrR3EP7t0KR8/t7XL31I4i4CLqXN+2QsBaK8iAquuLBYTofFslOBEleJmQllYX+2hN+Fvg30SvlG4ggadU0gJljonh58R42HTxBbkEx325K5K/jutC3Q/1qTY38fPjiluHkF5XQJMD86xlOH8yvuRLadnK/sK6hIECTS16EjDj46R4I7Q7hdkIhfhn4+EPEORUPbhkJPn7lsgCXYkuN4kyQ+DXS/XlYn6SuuHlUJJ+uPciz/9vJ8ax8okKbcPe5DeO79PP1wa8etSKDwRuYX/TphK8/XPEJBLeBr6+FrGNl+/Ytg47DoZETM4mvP7Tq6jxyy7bNprU4EhUNqXG6DG8DoUmAH/ec143NCekknsjl+an9jX/BYPAiRpCcbjRpBdO/hLx0+OY6KMqH7GRdQ6TLONfHhXZ3IUhidQ2S5i7WXXSxTGUeFruqK6YP68iwiBbcGd2lWrmrDAaD5xhBcjrSth9c+jYcWg+/zIL9VgR11HjXx4T1hBP7tVPenpQ90Kob+Lqwgob1hOB2Dc685e/rw7e3n839F7jQpAwGQ61hfCSnK30ug6N/wqr/aG0hsDm0H+i6fVgPUCWQtg/a9CnbnhwL7Qe5Pk5Er0vZ8yuUlOgFjgaD4YzC/NefzkT/E7pdABmHdA0Rn0r8BDZnuq0uO0BhLpw46No/YiMqGnLT4OjWmo+5tlBKr3E58Ht9j6Thk3UMljwJxUXu2xoMTjCC5HTGxwemvg/dJ8KQGypv26orIDoE2EZqHKC0/6Qyup6rI8K2fVvTEdcea96AJU/Axg/reyQNn13zYfXLcGRLfY/EcIpiBMnpTmAIXPM1dD2v8nb+jaFF5/IhwKURW05Cf+1pEgq9LoKtX2otpr6JWwJLHtfvM+onM8EpRUai/ltJvjWDoTKMIDGUEdaz/MMkeTeIj6WtuGHoTZB7QldRrE9S98F3N0FYL+h1cdlD0uCaUkESW3k7g8EFRpAYygjtrs1ZNlt5ciy0iAS/APfHRozWAmejYwLnOiQ/C766Rgu/q7/UqV+yjhjbvzsyLa0tZU/l7QwGF7gVJCJyl4jUf55wg/cJ6wnFBXDigP6csse9o92GCAy5UYccH9vhtSG6pKQE5t0OKXvhio91BceQcFDFkO26QqEBY9oy1BhPNJK2wB8i8o2ITJTaLoFnaDjYfCEpu6G4UJuJwtw42u0ZeA34BsDGj7wzvspY+QLE/gwTnilLkx8Srv8aP4lrSoohM0kvOj1xoGH4uAynHG4FiVLqn0A34ENgBrBXRP4lIl28PDZDXWOLzkqOhbT9UFLouUYCENQS+lwK276GgpPeGaMzdv0My5+DAdfAiDvKtpcKkkN1N5ZTjayjWmvrNAJQVqSewVA1PPKRKF205Kj1KgJaAN+JyAteHJuhrglsBsHtdQiwLXrLXeivI0NvgvxM2P599caw+VOYHa01Ik/IToZ5t0H7wXDRK9rEZqOZVZPdONxdY7s3Xc7Vfz0xbykFH02uH83T0CDxxEdyt4hsAl4Afgf6KaXuAIYAU90cO1FEdotInIg85GR/ZxH5TUS2ichyEQm3tkeLyBa7V56IXGrtixSR9SKyV0S+FhFTKq42CeuhNRJbBE9VBUnH4TpiqroPmdj/QdJm2L3As/Yxn0FBNlz2Lvg7lNcNbAYBIWXOZENFMi1BEjVOByl4IkiyjsDB1bB3sTdHZjiF8EQjCQUuV0pdoJT6VilVCGDVS7/I1UEi4gu8BUwCegNXi0hvh2YvAZ8qpfoDTwHPWX0vU0oNVEoNBMYDOcAi65h/A68opboBJ4CbPbtUg0eE9dAO6+TdENIRAqpYllZEayVJmyEppmrHKlV2zCYPBFFJCWz6WEeMuVrrEhJuNJLKsN2bllE6Qs+TEGDbd2TChQ0WngiSX4A02wcRCRaR4QBKqV2VHHcWEKeUildKFQBfAVMc2vQGfrPeL3OyH2AasEAplWM5+scD31n7PgEu9eAaDJ4S2h0KT0L88qprIzb6Xwl+jauulWQdgexjENIJ9i2FtPjK28cvhfSDMPRG121COhhBUhkZh7XWFtjMmkR4EAJsEyQn9uvs0oYzHk8EyTtAtt3nk9Y2d3QA7L2cidY2e7ZSZh67DAgWkVYObaYDc633rYB0pZRtYYCzPgEQkZkislFENiYnJ3swXANQ5lw/mVw1R7s9jZtDv6nw53e6wqKn2B5QE54G8dXaRmVs/AiCQqHnxa7bGI2kcjIStbAFLUhS49z7p2zfkyrRkX2GMx5PBIlYznag1KTlSdZgZ2HCyuHzLGCsiMQAY4HDaGe+7kCkHdAP+LUKfdrGOVspNVQpNTQsLMyD4RqA8iaiqoT+OjL0Jq3Z/PmN58ckxWg7fbcJ0GMSxHzhesabmaT9KIOu1ZUaXRESrhNKFuRUbfxnChmHyqLbQntASZGO2HOFUpC0BdoN0J+NecuAZ4Ik3nK4+1uvewA3NgdAawsd7T6HA0n2DZRSSUqpy5VSg4BHrG0Zdk2uBObZ/DJACtBcRGyCrEKfhhrSJBSCLKWwuhoJ6Ciqtv1h48f64eMJSVu0o75RkDZX5aTArp+ct938mQ5bdZeMspn1kDQOd+dkHi4TJPbriFyRkai/l35XAGJWwxsAzwTJ7cDZaG0hERgOzPTguD+AblaUVSO0iWq+fQMRCRUR2xgeBhzza1xNmVnLFoa8DO03AbgBqOfkTqchodYDpbo+Eihzuh/7ExI3um9vc7Tbap9EjddVGZ2Zt4qLYPMn0GW8dhJXRm2sJSkpgcRNngtE0KWHT6ZW/5x1QUEO5KSWhUnbryNyhc2s1Wmkzh5QHxrJiQM6r5uhweDJgsTjSqnpSqnWSqk2SqlrlFLHPTiuCLgLbZbaBXyjlNohIk+JyCVWs3HAbhHZA7QBnrUdLyIRaI1mhUPXDwL3iUgc2mdi8oTXNuFD9AM6qIYlavtNA/8msO0r921tM11b8S0fHxgyAw6sKp/aHiBusZ5JD73Jfb8htbCWZOlT8MH4skqT7lAKPr4Ifrq7+uesCzItZT7EMhwENNUanOP9ticpBnz8dPGzsB6Vt/UWn10Gvz5S9+c1uMStr0NEAtEhtn2A0kB9pZTb/2Kl1C/oqC/7bY/Zvf+Osggsx2MP4MSRrpSKR0eEGbzF+Edh9Kya9xMQDBGjPKvnbpvp2ldjHHQtLPuX1kom/qts+8Y50LStrrPijuD2gFQ/Tcr273WBLNDp6aPGuT8mZa+OJstJ1Y5rX//qndvb2LS0ELt/M9s6Ilcc2QKte+myA6HddXRdcZHrUsy1TUmJ1vbEE2OKoa7w5Nv4DJ1v6wK0dhAOZHlzUIZ6xi9AR17VBl2idfne9ITK2x3ZUjbTtdG0ta5zsuWLshxQ6Ql6Idzg6z17QPs1guC21dNIjmyDH++EjiO0KWffcs+Os9WvL8iGxD+qft66wnZPbOY/KFtHVFJSsb3N/NjO0hptST7TD3p/rDby0nVAQGoc5GW4b2+oEzwRJF2VUo8CJ5VSnwCT0ZFUBoN7oqL1X3daSVJM2UzXnqE36YeHrc7Jpk+0/2Xw9Z6PoVmHshXcnnIyRaekD2oJV32mI8mO/QnZbq26+lqD2+lZsyfaWH2ReRgQS2uzCOsBRbmQ4UTwpx/Uvgmb1lhanrkOswbb3/8j2+ruvIZK8USQ2CKm0kWkLxACRHhtRIbTi7Ae+kEVX8kD1dHRbo99nZPiQp0SpdsEaN6xYltXVHUtSXEhfHODXktz1edaM+piCcR4N36S4kLt1+kxSUeuVXbd9U3GIWjapnz4tC3Qwpnvw9H86IlzvrY5aSdIqpo5weA1PBEks616JP9ER13tRKcpMRjcI6IfwvHLdcpyZzjOdB2Pt9U5WfmSXvnuiZPdHpsg8TTq6td/6FxSF78OHQbrbW0HQOOW2idQGYl/aJNWVLS+7sObIDe9auOtKzIOlzdrQeUhwElbwMe/zPxoS/JZlyHANo1EfI0gaUBUKkis0NxMpdQJpdRKpVSUFb31Xh2Nz3A6EBWtBcWRrc73O3O022Orc7LieR1h5K7+vCMh4VCUBzlp7ttu/gw2zIaRd8GAq8q2+/hA1FitYVQmkPYt0yatyDH6ulUJ7F9ZtfHWFfar2m0EtYQmYc61jKQYLUTsK2aGda9jjcTKUtFppBEkDYhKBYm1iv2uOhqL4XQlapz+68rMY5vptnbM6Wlhq3MCMPgG8PGt2vk9XUuSFAP/u08LgPOerLg/KlrnA6vMJxC/TJu0GjeH8GE6/LkhmreUsgSJExNhqJOwXtuKdlt4to2wnrptVdbY1ITs41ob6RKtc315Yz2JUrD61dpJ/5J5BFa+6Dx44TTCE9PWYhGZJSIdRaSl7eX1kRlOH5qGQZt+rh3Pzma6jpx9t56FulvJ7gzbgjt3q9u3zNWRY9PmOA9nLfWTuLiO3HRtyrK182vkefhzXZN7QjvVmzlJVRfWQwtLe+GQFg/5GRW1RluSz7rKZ3byuNaYbCZHV1puTcg+DksehzVv1Lyv7d/B0mcgdW/N+2rAeCJIbgLuBFYCm6yXB0uVDQY7uozTfg7HnFelM10XZi0bbfvCTQu147uq2Gbd7h52ttBWVwsxm3eCll1cC4YDq7QpyxapBnr1/Yn9ejV2Q6J0DUl4xX1hPbTQyD5Wtu3IFv3X8XuypdGpLK1KbZKdrCcmthBkb5i3Tli5xmpDk7T95hra91/LeLKyPdLJy01eCoPBgahovebg4Jry213NdGuTJqHax1KZaau4CI7+6X4cXaLhwGooKqi4b98ybcoKH1a+vW1fQ8K2QNOVIIHyvo+kGH0Pw3q5aFtHguTkcWjSWgv7FhFeEiQHyv5WlsDSE4wg0YjI9c5edTE4w2lE57P1g8hxlufO0V4biFh1SSoxbaXs1qYeRx+AI1HR2pTjbKFh/DJtyioXTtvdffhzfeBsMaINZyHASVu0VuiYablJqI5mqytBkp1cppW2H6THVdvYC4+afm+2+1xTgdTA8cS0NczuNRp4AriksgMMhgr4N4ZOIyrOzI9s0QKmdS/nx9UW7taSeCrQIkdrZ6/jA+bEQa1ddYkuv700/HmF6/Dn+iAzUd/3oNCK+4Lb6mJXNo2kpMRKHe9CyIb1rJsQYKXKfCSgv6v0g55F41WFEwd0zrFm4TXXJI1GolFK/c3udSswCDB10g1Vp0s0HN8BWXa2d9tM19v5qJqFV+5sT4qBRsHaB1IZgSHQYUjFB4xNsERFVzwmKlqvzj/ihdlzdclIhGbtdVizIyI6rNcmHNLioSDLtZC1hQB7O3IrL0ObR20aibf8JCf2Q8tI7dfbv7L6E4DCXJ2I1NbnaUx1Mp/lAN1qeyCGMwDbQzZ+uf5rm+l606xlIyRch+66qv6XFKPNWs4erI50idY16e1DT21pUZzVjo8aV9amoZCR6NysZSO0R5m5yp22FtZT34uTKbU7Rkdsa0ia2ATJgPLjqy3S9mv/i20CUF3zmS27clCo1kjqKkS6HvDER/KTiMy3Xj8DuzE1QAzVoW1/XTTLNntP21f5TLc2CemgI6qyjlTcV1wIR7e794/YKF1ouEp/LinWKeajovVs3hFb+LNNgDoj5gt4++y6s6VnHHa+hsRGWA9tRspJ0w9qv0DXhc5qkirleCy8NdyzHGa2Nk0t01bj5lqDrE1BUnBSX3eLCLv1T26yGbjCFtwRMUoviM06WgsDbJh4opG8BPzHej0HjFFKPeTVURlOT3x8IHKsnpnb8mtB3Wkk4NzhfnwXFOe79gE4Ej5Um8FsAvHIVj0jd/SP2NNlHCSs0w8qRxLWwU/3aLPfV3+B/GzPxlFdiosgK6niqnZ7SlOl7NHfU9v+rlPF1yQEOH65FkCeJGC05dlqYo1ES7EAACAASURBVBcC3n5Q7a4lsfkyWkbqQIK2/T3P+uyI7bcWMap836chngiSBGC9UmqFUup3INUqOmUwVJ0u0ZB9VD88kraAX+OyKCFvUtlakqoKNF//8gsNS/0j41wfExUNJYUVw58zDsPX1+kklNM+guRd8OMd3jWDZB/VGlVlpi2bIDm+E45uq1xba9ZeC9bqRG7ZhI8nFSyzLdNWUwdBknGobF9NsT3sW0Tov12i9fqn6gh322+t8zlW36evn8QTQfItYL++v9ja5hYRmSgiu0UkTkQqaDEi0llEfhORbSKyXETC7fZ1EpFFIrJLRHbahJeInCsim0Vki4isFpGunozF0ECwTyufFANt+9VNUaTS1e0uBElAiPuyvfbYUnScOKCvpU3fyhdL2sKf7f0khbnw9V+gMAemz4W+l8P5T8Gu+TpBpbewPeCaVSJIQjppIR/7P52EsjIhKwKh3aonSGzHuMs6AFojER9tHrVhE3C1FchgMy22iNR/XU0APCEzUWdXbtVVj/sM10j8lFKlq6+s926jtkTEF3gLmAT0Bq4WEcdkSi8Bnyql+gNPoU1nNj4FXlRK9UJXRLQZUd8B/qKUGgh8ic5KbDhVaN5R/2Pt+02bJOrCrAW6jGxgc9caSfuBzv0brrAJxNhf9Iw1alzl7W3hzzbtRSltzkqKgctnQ2vLPDTyLuh/FSx7BnYv8Hw8VaGyNSQ2fHy0cLBlO3b3PYX1rJkg8STFSvZxLUTsc6217Q9I7flJThzQk4rGLfTnTiO1f6g660kyEvUExq+RFtqn8VoSTwRJsl2NdURkCuBJeMZZQJxSKt4SPl8BUxza9AZ+s94vs+23BI6fUmoxgFIqWylly62hgGbW+xAgyYOxGBoSUdEQ95te2FdXggScryUpyodjO6o+jtBu+iGx+hUdklqZf8RGl2htKso6Cmvfgm1fw7h/QM/JZW1E4OLXtL/m+1u9s9CvVJBU4iMBbd5SJeAfVOZQd9m2uzaZVaVq4cnUsvBYTwTJyeTy/hHQqexDu9XewsQT+6FlRNmkwj/Qqo5ZTUFiE9YtOp/xGsntwD9EJEFEEoAHgds8OK4DYG/4TKRiDfatwFTr/WVAsIi0ArqjC2n9ICIxIvKipeEA3AL8IiKJwHXA885OLiIzRWSjiGxMTq4l+6mhdugSjZ4PUA+CxMGEcnynNl14GrFlQ0QLxJPHwbcRdDrb/TE2Lea3p2Hxo9DrYhhzf8V2/o1h+hf6ITb36tqvZ5KRqNfDBARX3s7mu2o3wH3GZZvD3VlBLFfY/CONW3iukdgituxpP6h2NRKbf8RGl2jtu8p0EvHnCqXKR8a1jDyzfSRKqX1KqRFo7aGPUupspVScB307sxM4ehBnAWNFJAYYCxwGigA/9Cr6WegV9VHADOuYe4ELlVLhwEfAyy7GPVspNVQpNTQszMmPz1B/RIzSq8P9m+jZZF3RrENFp25NIsdsWkinEdAoyH17W/jzls/1g/fSd12vWwkJhys/0zXqv7sJ4pZU7VVZOG2mm9BfGzaHuyfRbNUJAba1jYrWY3KXat2WZ8uRdgN1FFpNw2tLinWGApt/xEbp+qcqaCW5J7TGbdP6WkRojSo/q2ZjbKC49XKKyL+AF5RS6dbnFsD/KaXc+SYSAftfazgOZiil1P+3d+5hVtxlnv98ae7QQAMdbs1dIMEEISF4jaEzq5Ook6isGzPjLa66M07Wcdw8Thyfzbhx8zgz0dFnHWccM4MazU6MMYkxoiZDQDPjGkMCkgAhkoBCOkAjBNLNteHdP36/6i4O51KnT58+fXk/z8NzTlX9qupX1Yd66723AO+Mxx0LrDKzw1Hb2Ghmz8dt9wOvkfQA8Cozeywe4jvAj0tdg9PHGDk+CpMh5fcWqYTxTSHB7ERb8JlAECSjGmDC7PKPN29lsJ8vvCrb+CFDQpvg7T8KGkcyh0LMfi285W/hwT8PPqVyaLoUPvRv+bcd3p2/fHwu05YEgT87g7bVMCcEE5QTAtz6bHiZmPUa2HJvMHMVClgwO7vOVprkJaBlEyy6Mvv5cznSErTTXI1kyoUhqfC5daHJWqZj5RTFTITTod+ESg4DjCzhMleZ2V8mC2Z2SNJbKO3kfhxYIGkuQdN4N3DWX0HSZOBgbKD1KWB1at8GSY1m1gpcQShdfwgYL2mhmT0LvAnYluEanL7Gtd/q/XMmb+FHXuh6204y68txtCeMmQwf2xgic7Ly1i/Amz6b30STj+UfhNlvCAIwK9t/BP/+dyGYIcn+TnN4DzStKH2chjnw8c3ZhM6QuvIjtw5sD/ukG48VEiQn20JRzTF57tvUi8JLyYsVCpLE9DQxRyMZMiS8NDy/Pgi0LL+V3Mi4RDgd2jloBUmdpBFmdgJA0iigSAeigJl1SLoB+AlQB6w2sy2SbgE2mNkDwErgc5KM0O/kT+O+pyXdCKyVJEIPlNvjMT8MfE/SGYJgKbOBt9MnGDm+98+ZmBkO7w6C5NTx4CN53ce6f8xx08sbP3xM+FcOjSUc3blMXgC/+AfY8HX4gy+dve1kezC7lHK0JxSL7DrnvAtDY6+stG6HOZednSw645L8Yzuz2vMImhFjgz+nUj9Jbg5JmvnNoUnVvi3ZBEFuZFwinAaowz2LIPk24YH+9bh8PfDNLAc3szXAmpx1N6e+3wPcU2Dfh4EledbfB9yX5fyOcxa52e37tsCZjt51+PcGoxrgwlXw1HfhzZ8926ne2Yckg4+kXBrPhy33heZlpXxGx49EzXBhtsZjuXW2cpm+LJj/smoM+Ti4M3TIzJdfk/aTZBUkdcO7NKhRDeHlaYCGAGdxtv8t8L+BCwgO9x8D3TAoO06NqZ8GqOuB1fJk+BxoggSCSexkWxAmaY5kyCHpLo0LAcvWVvZAHNN4fnjIDh1VPCkxt85WLtOXho6O+WqpZeXQrtAFM1+C7PgZQePKGgacr7pyw9wBq5Fkrf67l5Ddvgr4Pdwv4fRH6oYFYdIpSDYFJ2o1Hqq1ZsYloVDkhtVnl1vptN1nNG2VQ2cIcAY/SeKUn7woNh5rKl4mJV+drTRph3t3ObQzv1krYV5zyHA/dbz0sfJFxjXMGbAhwAUFiaSFkm6WtA34e0JOiMys2cz+vtdm6Dg9yfimrrfyFytwtPd1JFh+fWgf/MKTXesP7wFUvm8nCxPnhyivLIKk9Zlg+kke3ONnFDdtJbW0xuRpxAUhskp1lflJDu06N/Q3zfzm4PDf/VjhMQlJVnuaiXNDOHdfanDWQxTTSJ4haB9/YGZvMLMvE+psOU7/JXlgnTwaqv6Wm4jYn7joXSG8dsPqrnWHXwgdEKvRSGzo8PCwzJJL0vpsKJWTmJHyJYumaY/lUQrNe/jo0GWzu4Lk2EshCKGYRjLnDcGHUiqf5MzpEEqcq+k2zAk+uSzJl/2MYoJkFcGktU7S7ZJ+j/xJho7Tf0geWHufAjs9MP0jCSPHwZJ3wdPf68qOP7y7uqa8qReFyK1S1YsPbD+77Mr4mcHH0XEy//i2AsmIZ517Cex7urz5JqTLxxdiRD1MvziU/S/Gy3vDb+scQTJwI7cKChIzu8/MrgXOB9YTMsqnSPpHSW/upfk5Ts8yrin0HtnxcFgeyIIE4JLrgzlm83fC8pEXqitI5q0MDu9i5q1Tx8PDNN0oa9wMwEKGej7aW0vn3jQuCucup95XQuK7KKaRQFf/k2LmqUJFMdO5JAOMLFFb7WZ2p5m9jZCdvgnwxlZO/yT5z73twZBIWD+ttvOpNtOXhrfoDV+P9Z/y2O57kizlRH63IxSDTOfHdIZmFzD7ZNFIkiTTcup9JRwsQ5CcOtoVdZaPQpFx45uCaWwwaST5MLODZvZPZnZFtSbkOFUl+c/dum3gOtpzWX59uN5nHgwtX6uRQ5LQMDv0dSkWJpv4UNINzYp1sISokZQQJImprDudGg/tChF8pQpZdkaHFfHFFIqMG1IXwosHYC5JWYLEcfo96bfErK11+zsXroIR4+CRW8Ny1qz27jKvGXb9e2F/x4FnQ0mTSamedONSVQdyOdkecmLylUdJk9T76k7v+EM7i/tHEiYvCAEMRQXJC6Gnychx524boLkkLkicwcXoSaHQIgx8/0jC8DGhWVZrTP+qdt7M/OZQ+XbP4/m3tz4THvrDRqbmODr8bfKZtoqVR0nTWe+rG6atfOXjC51j2pLSGkkhYT1Ac0lckDiDC6nr7Xcgh/7msvz6ru/FWuz2BHMuCxpHIT9J67NnO9o75zUjf3Z7qfIoaRoXlW/a6jgZHv7FckjSTF8Wov5Od+TfXiwybuLcEAxw7FB5c+zjuCBxBh/jm4KTvX5qrWfSe0x5Jcx8ddDGCiX19RSjJoTM+nx+ktMdwdmer+Pi+JklNJIMFZMnLwql2k8dyz7fw7uD8z+LRgJBkHQcKyywikXGJcKqN/wk2x6Ev5lbnS6bObggcQYfzZ8O7WwHG2/9Alz95d4JMJjXHGqZ5b55H9oZen40Ljp3n/Ez8jvbS5VHSZPU+yoWVZVLofLxhSjmcD95FI7+rnBkXGcI8K7s8+subXvh2MHgH6syLkicwcesV8PC36/1LHqfqRfBkv/SO+eaf0V4y9/56Nnrk7fjvIKkCU4cPjcPpLM8SgaNpJx6XwnFysfnY+J8GF6fX5AciXkwhSLjejOXpFRZmR7EBYnjOD1P0/LwsM31k3SG/uYzbRUIAW7fDyMnhBIspUjqfZXjJzm4M5j8xmY0dQ4ZEhqG5RMkSdRZIWf7iLFBIPaGRtK+H0ZNrE45nBxckDiO0/PUDQu1qXL9JAeeDc7+fPkaSRBArsO9bX/piK2EocNDHks5IcCHdoVWy0PKeBxOXwp7n4bTp85en9tiNx8Nc3vHR1LOfauQqgoSSVdK2i5ph6RzsuElzZa0VtJmSeslNaW2zZL0kKRtkrZKmhPXS9Ktkp6N2ypob+c4TtWY3xxMOOm379ZnCnd8TLfcTdPems0/ktC4qLwQ4EO7svtHEqYvC6V29ud01EiqK9cXqa7cMCcEBFSb9tZs5sAeoGqCRFId8BXgKkJDrOskLc4Z9nngDjNbAtwCfC617Q7gNjO7AFgBRI8bHwBmAufHbXdV6xocx6mApFxKopWcOROc4JPz+EcgRNGp7lzTVtv+7D3uIQiSg8+dqy3kw6x0+fh8FHK4H94dSu8UM8NNnBvKqBRK2OwpBohGsgLYYWbPm9lJwgP/mpwxi4G18fu6ZHsUOENju13MrM3MjsZxfwLcYmZn4rb9OI7T95i8IEQvJX6Sw7tDnap8jnYIyX7jpp8bAlyuRjJ5USjXfvD50mPbD4Ss+ayO9oSJ80L2+jmCJENRzIY5IRChWCOvfBw92FXFOQvl3rcKqKYgmUFohpWwJ65L8ytCuXqAdwD1kiYBC4GXJN0raaOk26KGAzAfuFbSBkk/krQg38klfSSO2dDa2tpjF+U4TkakoJU8/9NQLfdANDcVEiQQG4+lNJJTx+HEkTI1kmg6yxK5laV8fD4kmJ7H4V4sqz2hu7kkd78P7v9otrEnjwYBWc59q4BqCpJ8weq5TQpuBC6XtBG4HHgB6ACGApfF7ZcC8wgmLYARwHEzWw7cDqwmD2b2NTNbbmbLGxt752Y6jpPD/GY4/lJogdsZ+psnqz1h3Iyz39TLySFJmFyOIMlY9Tcf05fBvi3QcSIsm+VvsZtLd0KAT5+C3b8MyZxZ6M59q4BqCpI9BF9GQhNwVrMBM2sxs3ea2TLg03Hd4bjvxmgW6wDuBy5OHfd78ft9wJLqXYLjOBUx9/Lw+fwjwdE+ejKMnlh4fNJ47MyZsJzkQpRj6x8+BsbPyhYCnGgkE2ZnP37C9GUhuXL/1rB87FAw3ZUybdVPDeHG5YQA798WnPsv7802vjv3rQKqKUgeBxZImitpOPBu4IH0AEmTJSVz+BRd2sXjQIOkRJW4Aoh/Le6PyxC0mG5UaHMcp1cY2xgSIZ9bH0xbxcxaEB7CZ0511dfq7pt146JsIcAHd4YIq3QByazkOtwTTapUvxcpRm7tyn6uFzeFzxOHQzXkUnTet35u2oqaxA3AT4BtwN1mtkXSLZKujsNWAtslPQtMAW6N+54mmLXWSnqKYCa7Pe7z18CquP5zwIeqdQ2O4/QA85ph92Owb2s2QQJdDvdy6mylaVwEB3Z0aTaF6E7ob8KE2TCqISVIMuSQJJSbS5L2xWTRSrJWTO4hhlbz4Ga2BliTs+7m1Pd7gHsK7PswecxWZvYS8NaenanjOFVjfjP8/P/AyVOFQ38TkofwkT3AJd3XSCYvDIUVD/+2uP/j0M5QzqU7SKGnTacgKdAZMR8Nc2Dnz4JfJUvts5aNwRzWcTwIkknzi49vL6OsTA/gme2O41SXWa/t6gFTSiPpbHCVaCStoehguaanLDW3Th0LPd7LzSFJM31Z8F+cOh6EX92I4AcqxcS5oWdLe4aI0o4TIYt+3sqw/PKLpfdp2w8jx8PQEaXH9gAuSBzHqS7DRgVhAqUFyaiG0IEwESTt+7v3Vp0lBDjJLu9OxFbC9GUhZ2Xflq7Q3yylVsqpArx/a/AbLbwyLGcxbbVn6HHfg1TVtOU4jgPAxe8FLPSBKYYUI7dSGkl37PyjGkKGeTFB8twj4XPqReUfP6HT4f5kmHMpR3tCojG1bIKZK4qPTUxn81YGzS6TRtIarr+XcI3EcZzqc+EqeN/3s/kDxs+oXCOB4CcpFAJsBhtWQ9OlcF6RvJZSjG8KpqyWTTGrvUQOSULD7KCVFOoimaZlUxCMDXNC6HDbvtL7tJdZVqZCXJA4jtO3SGe3V1IvqnFR0EgsNw8a+M1/wO9+DZdcf+62cpBCJeA9j8PLLaWz2tPMaw79WkrVBGvZGDQfKWh0maK2eq88CrggcRynrzGuKbx1n2gLWfHdfSA2nh/Kq+R78G5YHZzRr3xHZXOF8JA/sD3Uz8oSsZUwvxlOvgx7NhQec+p48JEkJrT6qaVNW6eOh3wT10gcxxm0JA/jF38VPrv7QExKpeSat9paYesD8Ko/hOGju3fsNMlDHrp6qmRh7htBQ4qbt/ZtCc78aUvDchaNpDP01zUSx3EGK4kgaXkyfFaikcC5DvdNd4YoqOUVmrUS0oKkHI1kVEPYN7f5V5rkHiTnGDslFGM88XLhfdp7NxkRXJA4jtPX6BQkMVqpuw/EsecF81VakJw5A098HWa/vnQoclbqp3VFSJXjI4HgJ3nhiXP71Ce8uCk485N7kkS9FdNK2lwjcRxnsDMudhd8IdFIumnakoJWkhYkO9eH3I1Kney555m2NAitfC2EizG/Gex0cLrno2VTl6Mdgo8EivtJ2rtZVqYCXJA4jtO3GDYqvIUnZdYrMdHkhgBvWA2jJ8Hiqwvv0x0u/wt4yxfK369pRUjAzOcnOXk0ZM1PX9q1LpNG0rsl5MEFieM4fZHElDNsTCgL310aFwXn89GDcORFeGYNLP3Dni8d0nQJLHlX+fsNHQ5zXp/fT7Lv6aCtpH0wmTSSbpaVqQAXJI7j9D0SQVKpwzjtcN/47fBg7kmzVk8wrzn0mH/pt2evT3xEaUEyoj4I11IaSS8Va0xwQeI4Tt+jpwRJEgK8fws88Y3QaKtU5dzeZn5z+MzVSlo2BSd+uqyMFHNJSgiSXozYAhckjuP0RRJBUumb9fiZMGw0PPZPoTrv8g9WPreepvH8ICxy/SQtG4MTP7esTKlckkrKynSTqgoSSVdK2i5ph6Sb8myfLWmtpM2S1ktqSm2bJekhSdskbZU0J2ffL0tqq+b8HcepET2lkQwZApMXhO6MY86D8/tgKyMpFGR8/qddjbhOtIUggbRZK6FUdvtA0kgk1QFfAa4CFgPXSVqcM+zzwB1mtgS4hdDxMOEO4DYzuwBYAexPHXs5MKFac3ccp8YkGeI9EXmUNNO6+L1QN6zy41WDec1w7CDsjdn8e58KJVcKCpK9+WuIdZysrKxMN6mmRrIC2GFmz5vZSeAu4JqcMYuBtfH7umR7FDhDY5dEzKzNzI7GbXXAbcAnqzh3x3FqycS5MGRoZb1CEqa9Khzr4vdXfqxqMW9l+Ez8JJ2O9qXnjq2fGro/5ktiTMqj9GIOCVRXkMwAdqeW98R1aX4FrIrf3wHUS5oELAReknSvpI2SbosCBEIf+AfMLENRfsdx+iVjJsNHH4OLuhFSm8ulHwrHaphd+bGqRf0UOO+VXX6SFzcFX0gS7nvW2CK5JN1tTVwh1RQk+RoP5OpiNwKXS9oIXA68AHQQGm5dFrdfCswDPiBpOvAu4MslTy59RNIGSRtaWzO0s3Qcp28x+RVQ1wO994aNDMfq68xvht/+IiQiJqXj81EslyQpjzJQfCQEDSTd5aUJaEkPMLMWM3unmS0DPh3XHY77boxmsQ7gfuBiYBnwCmCHpF3AaEk78p3czL5mZsvNbHljY++qeY7jOGUzrxlOn4RfPwQHfl1EkGTRSAaOaetxYIGkuZKGA+8GHkgPkDRZUjKHTwGrU/s2SEruxhXAVjP7oZlNNbM5ZjYHOGpm/eBVw3EcpwSzXwd1w+E/vgRYYUGSFIjMq5H0fuVfqKIgiZrEDcBPgG3A3Wa2RdItkpJCNyuB7ZKeBaYAt8Z9TxPMWmslPUUwk91erbk6juPUnOGjYearuxzt0/I42gFGjA0lUPK13G1vrbysTDfoAQNkYcxsDbAmZ93Nqe/3APcU2PdhYEmJ44/tgWk6juP0DeY3w65HQ/hzscirQrkkbb3bqz3BM9sdx3H6CvNiuZR8Yb9pCpVJad/f6xFb4ILEcRyn7zBtKcx6LVxQosx9/bTCUVu97B+BKpu2HMdxnDIYMgQ++OPS49LZ7elaXO37YdZrqje/ArhG4jiO098YOzWECh871LXudEfou1IDjcQFieM4Tn8jX1Li0QOA9XoOCbggcRzH6X90JiWmBEmNckjABYnjOE7/o1MjSUVu1ajOFrggcRzH6X/kM23VqM4WuCBxHMfpfwwbBSMnwMup7PYa1dkCFySO4zj9k9xckrb9MHQkjKjv9am4IHEcx+mP5Ga3t7cG/0huj/dewAWJ4zhOf6R+2tmCpEZ1tsAFieM4Tv+kfiq07YUzZ8JyopHUABckjuM4/ZH6qXCmA47+Liy7RuI4juOURToE+MzpkNnuGonjOI6TmXTL3aMHwc7UJIcEqixIJF0pabukHZJuyrN9tqS1kjZLWi+pKbVtlqSHJG2TtFXSnLj+znjMpyWtljSsmtfgOI7TJ0lrJDXMIYEqChJJdcBXgKuAxcB1khbnDPs8cIeZLQFuAT6X2nYHcJuZXQCsAOKd4k7gfOAiYBTwoWpdg+M4Tp+ls3f73lSdrSk1mUo1NZIVwA4ze97MTgJ3AdfkjFkMrI3f1yXbo8AZGtvtYmZtZnY0fl9jEeCXQBOO4ziDjaEjYPSkELnVXrvyKFBdQTID2J1a3hPXpfkVsCp+fwdQL2kSsBB4SdK9kjZKui1qOJ1Ek9Z7gbxdYCR9RNIGSRtaW1t74HIcx3H6GEkuSdsANW0B+dIrLWf5RuBySRuBy4EXgA5C58bL4vZLgXnAB3L2/QfgZ2b2aL6Tm9nXzGy5mS1vbKzNzXUcx6kq9VO7fCR1w2Hk+JpMo5qCZA8wM7XcBLSkB5hZi5m908yWAZ+O6w7HfTdGs1gHcD9wcbKfpL8CGoFPVHH+juM4fZukTEpb7cqjQHUFyePAAklzJQ0H3g08kB4gabKkZA6fAlan9m2QlKgSVwBb4z4fAn4fuM7MzlRx/o7jOH2b+mnQti9oJTVKRoQqCpKoSdwA/ATYBtxtZlsk3SLp6jhsJbBd0rPAFODWuO9pgllrraSnCGay2+M+X41j/5+kTZJurtY1OI7j9GnGTgn5I/u31iwZEYIvomqY2RpgTc66m1Pf7wHuKbDvw8CSPOurOmfHcZx+Q5KU2LZvYGokjuM4TpVJBAnUVCNxQeI4jtNfSbLboWY5JOCCxHEcp/8y9jw6My1qlEMCLkgcx3H6L3XDugSIaySO4zhOt0jMW+4jcRzHcbpF4nB3jcRxHMfpFvVTYchQGDmhZlPwnAzHcZz+zCUfgPMWw5Da6QUuSBzHcfozMy4O/2qIm7Ycx3GcinBB4jiO41SECxLHcRynIlyQOI7jOBXhgsRxHMepCBckjuM4TkW4IHEcx3EqwgWJ4ziOUxEys1rPoepIagV+U2LYZOBAL0ynr+HXPbjw6x5cVHrds82sZH36QSFIsiBpg5ktr/U8ehu/7sGFX/fgoreu201bjuM4TkW4IHEcx3EqwgVJF1+r9QRqhF/34MKve3DRK9ftPhLHcRynIlwjcRzHcSrCBYnjOI5TEYNekEi6UtJ2STsk3VTr+VQTSasl7Zf0dGrdREkPS/p1/Gyo5Rx7GkkzJa2TtE3SFkl/FtcP6OsGkDRS0i8l/Spe+/+K6+dKeixe+3ckDa/1XHsaSXWSNkp6MC4P+GsGkLRL0lOSNknaENdV/bc+qAWJpDrgK8BVwGLgOkmLazurqvIN4MqcdTcBa81sAbA2Lg8kOoD/YWYXAK8B/jT+jQf6dQOcAK4ws1cBS4ErJb0G+Bvgi/HaDwH/tYZzrBZ/BmxLLQ+Ga05oNrOlqfyRqv/WB7UgAVYAO8zseTM7CdwFXFPjOVUNM/sZcDBn9TXAN+P3bwJv79VJVRkze9HMnozfXyY8XGYwwK8bwAJtcXFY/GfAFcA9cf2Au3ZJTcBbgX+Oy2KAX3MJqv5bH+yCZAawO7W8J64bTEwxsxchPHSB82o8n6ohaQ6wDHiMQXLd0cSzCdgPPAw8B7xkZh1xyED8zX8J+CRwJi5PYuBfc4IBD0l6QtJH4rqq/9aH9vQB+xnKs87joQcgksYC3wM+bmZHwkvqwMfMTgNLJU0A7gMuyDesd2dVPSS9DdhvZk9IWpmszjN0ia2JLAAABYBJREFUwFxzDq83sxZJ5wEPS3qmN0462DWSPcDM1HIT0FKjudSKfZKmAcTP/TWeT48jaRhBiNxpZvfG1QP+utOY2UvAeoKfaIKk5CVyoP3mXw9cLWkXwVR9BUFDGcjX3ImZtcTP/YQXhxX0wm99sAuSx4EFMaJjOPBu4IEaz6m3eQB4f/z+fuD7NZxLjxPt4/8CbDOzv0ttGtDXDSCpMWoiSBoF/CeCj2gd8J/jsAF17Wb2KTNrMrM5hP/Pj5jZHzGArzlB0hhJ9cl34M3A0/TCb33QZ7ZLegvhjaUOWG1mt9Z4SlVD0r8CKwmlpfcBfwXcD9wNzAJ+C7zLzHId8v0WSW8AHgWeostm/pcEP8mAvW4ASUsIztU6wkvj3WZ2i6R5hLf1icBG4D1mdqJ2M60O0bR1o5m9bTBcc7zG++LiUOD/mtmtkiZR5d/6oBckjuM4TmUMdtOW4ziOUyEuSBzHcZyKcEHiOI7jVIQLEsdxHKciXJA4juM4FeGCxBmUSFovaXnpkRWf52Ox8vCdOeuXxtDznjrP6Vjx9WlJ35U0Oq6fKukuSc9J2ippjaSFPXXeEnPqlXvs1B4XJI5TJqkM6Sx8FHhLTIpLsxToMUECHIsVXy8ETgJ/HJMx7wPWm9l8M1tMyKGZ0oPndRwXJE7fRdKc+DZ/e+yn8VDM0D7rbVfS5FgSA0kfkHS/pB9I2inpBkmfiL0pfiFpYuoU75H08/gWvyLuP0ahb8vjcZ9rUsf9rqQfAA/lmesn4nGelvTxuO6rwDzgAUl/nho7HLgFuDZqEdfGnhH3S9oc57kkjv2MpG9JeiT2k/hwhlv3KPAKoBk4ZWZfTTaY2SYzezTPfU73qLlR0mfi949FTWazpLtK3KNRUfvZLOk7wKgMc3UGAIO9aKPT91kAXGdmH5Z0N7AK+HaJfS4kVPkdCewA/sLMlkn6IvA+QiUDgDFm9jpJbwRWx/0+TSir8cFYXuSXkv4tjn8tsCQ3K1jSJcD1wKsJBQIfk/RTM/tjSVcS+kMcSMab2UlJNwPLzeyGeIwvAxvN7O2SrgDuIGgtAEsINbLGABsl/TCpqZRL1JauAn4cr+eJEveqFDcBc83sRFJuhcL36L8BR81sSRSET1Z4bqef4BqJ09fZaWab4vcngDkZ9llnZi+bWStwGPhBXP9Uzv7/Cp19WsbFh+KbgZsUSq+vJwijWXH8wwVKS7wBuM/M2mP/j3uBy7Jd3lnH+FaczyPAJEnj47bvm9mxKIzWEQrx5TIqznkDoQzGv5R5/kJsBu6U9B5CkzAofI/eSBTyZrY57usMAlwjcfo66XpIp+kyl3TQ9SI0ssg+Z1LLZzj7N59bH8gIGsUqM9ue3iDp1UB7gTn2RE36YqXO880zl2NmtjS9QtIWugoVFiN9L+Hs+/lWgoC4Gvifkl5J4XtUaG7OAMc1Eqe/sgu4JH7P8rDMx7XQWdjxsJkdBn4C/PfoqEbSsgzH+RnwdkmjFaquvoPgpyjGy0B9zjH+KJ5zJXDAzI7Ebdco9F+fRCi6+XiGOQE8AoxI+1UkXSrp8pxx+4DzJE2SNAJ4Wxw7BJhpZusIjaImAGMpfI/S13AhwSTnDAJckDj9lc8DfyLp54Rqxt3hUNz/q3T18P4soSXt5uiA/mypg8RWvt8AfkmoKvzPZraxxG7rgMWJsx34DLBc0mbgr+kq+0087g+BXwCfLeQfyTMvIwi1N8Xw3y3xPC05404RnP+PAQ8CSTOkOuDbkp4iVMz9YuxrUuge/SMwNl7DJ+O8nUGAV/91nD5MjJ5qM7PP13oujlMI10gcx3GcinCNxHEcx6kI10gcx3GcinBB4jiO41SECxLHcRynIlyQOI7jOBXhgsRxHMepiP8P4+WK4ozLGSEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: 'argmax' is deprecated, use 'idxmax' instead. The behavior of 'argmax'\n",
      "will be corrected to return the positional maximum in the future.\n",
      "Use 'series.values.argmax' to get the position of the maximum now.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "formula           y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...\n",
       "mean_test_acc                                              0.972755\n",
       "mean_test_mse                                             0.0272446\n",
       "mean_train_acc                                             0.972136\n",
       "mean_train_mse                                            0.0278638\n",
       "num_pc                                                           10\n",
       "Name: 9, dtype: object"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(df_res.num_pc, df_res.mean_train_acc, label='training set')\n",
    "ax.plot(df_res.num_pc, df_res.mean_test_acc, label='test set')\n",
    "ax.legend()\n",
    "#ax.set_xticks(df_res.num_pc.values)\n",
    "ax.set_xlabel('number of top PC used')\n",
    "ax.set_ylabel('Accuracy')\n",
    "plt.title(\"Accuracy for Numbers of Top PC Used\")\n",
    "plt.show()\n",
    "best_model = df_res.loc[df_res.mean_test_acc.argmax()]\n",
    "best_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explained Variance and Most Imporant Feature for each PC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.13063776 0.05241703 0.0446671  0.02854428 0.02477324 0.02168846\n",
      " 0.01967302 0.01873948 0.01759851 0.01754495 0.01530759 0.01482918\n",
      " 0.01430431 0.01343961 0.01316621 0.01263012 0.0120022  0.01189397\n",
      " 0.01155233 0.01141751 0.01114379 0.01103642 0.01072142 0.01055201\n",
      " 0.01033964 0.00989237 0.0098532  0.0096393  0.00932993 0.00904221\n",
      " 0.00880386 0.00869686 0.00850667 0.00846557 0.00831656 0.00816874\n",
      " 0.00805011 0.00779247 0.00777166 0.00756868 0.00740771 0.00727339\n",
      " 0.007166   0.00700152 0.00689503 0.00666616 0.00648185 0.0063594\n",
      " 0.00619028 0.00599489]\n"
     ]
    }
   ],
   "source": [
    "print(pca.explained_variance_ratio_) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_importance = pd.DataFrame(pca.components_,columns=one_hot.columns,index = ['PC'+str(i+1) for i in range(n)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PC1                            q0009_Not employed-retired\n",
       "PC2                                q0025_0003_No children\n",
       "PC3                                q0025_0003_No children\n",
       "PC4                               q0025_0003_Not selected\n",
       "PC5                               q0004_0002_Not selected\n",
       "PC6                               q0011_0004_Not selected\n",
       "PC7     q0020_0001_Read their physical body language t...\n",
       "PC8                               q0020_0005_Not selected\n",
       "PC9                                             q0005_Yes\n",
       "PC10                              q0020_0004_Not selected\n",
       "PC11                                             q0005_No\n",
       "PC12              q0020_0004_Every situation is different\n",
       "PC13                                    q0007_0004_Rarely\n",
       "PC14                              q0020_0004_Not selected\n",
       "PC15                                             q0017_No\n",
       "PC16                              q0020_0004_Not selected\n",
       "PC17                    q0004_0006_Other (please specify)\n",
       "PC18                              q0020_0005_Not selected\n",
       "PC19                                   q0004_0005_Friends\n",
       "PC20                                    q0007_0005_Rarely\n",
       "PC21                     q0007_0010_Never, but open to it\n",
       "PC22                                            q0005_Yes\n",
       "PC23                               q0008_0002_Your weight\n",
       "PC24                              q0020_0002_Not selected\n",
       "PC25                                    q0007_0004_Rarely\n",
       "PC26                 q0007_0005_Never, and not open to it\n",
       "PC27                             q0002_Somewhat important\n",
       "PC28                                     q0007_0008_Often\n",
       "PC29                                 q0007_0004_Sometimes\n",
       "PC30                               q0029_College graduate\n",
       "PC31                                 q0007_0003_Sometimes\n",
       "PC32                                    q0007_0011_Rarely\n",
       "PC33                              q0019_0003_Not selected\n",
       "PC34    q0020_0002_Ask for a verbal confirmation of co...\n",
       "PC35                                     q0007_0009_Often\n",
       "PC36                              q0019_0005_Not selected\n",
       "PC37                               q0029_College graduate\n",
       "PC38                                       age3_65 and up\n",
       "PC39                 q0007_0007_Never, and not open to it\n",
       "PC40                                 q0007_0003_Sometimes\n",
       "PC41    q0020_0003_Make a physical move to see how the...\n",
       "PC42                                 q0007_0009_Sometimes\n",
       "PC43       q0008_0007_Sexual performance or amount of sex\n",
       "PC44                                        q0024_Married\n",
       "PC45                               q0004_0004_Pop culture\n",
       "PC46    q0019_0004_Societal expectations, even though ...\n",
       "PC47                    q0004_0006_Other (please specify)\n",
       "PC48                                         q0018_Always\n",
       "PC49                                    q0007_0001_Rarely\n",
       "PC50                                    q0007_0009_Rarely\n",
       "dtype: object"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "feature_importance.idxmax(axis=1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  },
  "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": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}