polynomial_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By Qi Goh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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.model_selection import KFold\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_hot = pd.read_csv(\"one_hot_responses.csv\")\n",
    "one_hot = one_hot.drop(columns=[\"Unnamed: 0\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Very masculine, Somewhat masculine, Not very masculine, Not at all masculine\n",
    "#masculinity_scale = [col for col in one_hot.columns if col.startswith('q0001')]\n",
    "# Straight, Gay, Bisexual, Other\n",
    "#sexuality_scale = [col for col in one_hot.columns if col.startswith('q0026')]\n",
    "# print(masculinity_scale)\n",
    "#print(sexuality_scale)\n",
    "\n",
    "# top n pca\n",
    "n = 30\n",
    "predictors = [one_hot.columns[i] for i in range(len(one_hot.columns))]\n",
    "X = one_hot[predictors].values.T\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",
    "pcs = PCA(n_components=n).fit_transform(S.T)\n",
    "df_pcs = pd.DataFrame(pcs, columns=['PC'+str(i+1) for i in range(n)])\n",
    "\n",
    "# predict whether the person was masculine based on their survey responses\n",
    "df_pcs['y'] = one_hot['q0001'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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_mse</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.027619</td>\n",
       "      <td>0.027591</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>y ~ 1 + PC1 + PC2</td>\n",
       "      <td>0.027679</td>\n",
       "      <td>0.027568</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3</td>\n",
       "      <td>0.025894</td>\n",
       "      <td>0.025692</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.024912</td>\n",
       "      <td>0.024705</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.024190</td>\n",
       "      <td>0.023842</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.024338</td>\n",
       "      <td>0.023831</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.024420</td>\n",
       "      <td>0.023766</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.023849</td>\n",
       "      <td>0.023225</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.023147</td>\n",
       "      <td>0.022460</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.022716</td>\n",
       "      <td>0.022252</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.022764</td>\n",
       "      <td>0.022094</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.022957</td>\n",
       "      <td>0.022072</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.022738</td>\n",
       "      <td>0.022051</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.022730</td>\n",
       "      <td>0.022024</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.022404</td>\n",
       "      <td>0.021716</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.022885</td>\n",
       "      <td>0.021543</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.022468</td>\n",
       "      <td>0.021554</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.022365</td>\n",
       "      <td>0.021534</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.022386</td>\n",
       "      <td>0.021443</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.022196</td>\n",
       "      <td>0.021431</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.022401</td>\n",
       "      <td>0.021255</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.022048</td>\n",
       "      <td>0.020980</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.022185</td>\n",
       "      <td>0.020966</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.022020</td>\n",
       "      <td>0.020835</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.021864</td>\n",
       "      <td>0.020676</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.021536</td>\n",
       "      <td>0.020673</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.021429</td>\n",
       "      <td>0.020310</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.021425</td>\n",
       "      <td>0.020300</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.021431</td>\n",
       "      <td>0.020235</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.020982</td>\n",
       "      <td>0.019872</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              formula  mean_test_mse  \\\n",
       "0                                         y ~ 1 + PC1       0.027619   \n",
       "1                                   y ~ 1 + PC1 + PC2       0.027679   \n",
       "2                             y ~ 1 + PC1 + PC2 + PC3       0.025894   \n",
       "3                       y ~ 1 + PC1 + PC2 + PC3 + PC4       0.024912   \n",
       "4                 y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5       0.024190   \n",
       "5           y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6       0.024338   \n",
       "6     y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7       0.024420   \n",
       "7   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.023849   \n",
       "8   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.023147   \n",
       "9   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022716   \n",
       "10  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022764   \n",
       "11  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022957   \n",
       "12  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022738   \n",
       "13  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022730   \n",
       "14  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022404   \n",
       "15  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022885   \n",
       "16  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022468   \n",
       "17  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022365   \n",
       "18  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022386   \n",
       "19  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022196   \n",
       "20  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022401   \n",
       "21  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022048   \n",
       "22  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022185   \n",
       "23  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022020   \n",
       "24  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.021864   \n",
       "25  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.021536   \n",
       "26  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.021429   \n",
       "27  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.021425   \n",
       "28  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.021431   \n",
       "29  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.020982   \n",
       "\n",
       "    mean_train_mse  num_pc  \n",
       "0         0.027591       1  \n",
       "1         0.027568       2  \n",
       "2         0.025692       3  \n",
       "3         0.024705       4  \n",
       "4         0.023842       5  \n",
       "5         0.023831       6  \n",
       "6         0.023766       7  \n",
       "7         0.023225       8  \n",
       "8         0.022460       9  \n",
       "9         0.022252      10  \n",
       "10        0.022094      11  \n",
       "11        0.022072      12  \n",
       "12        0.022051      13  \n",
       "13        0.022024      14  \n",
       "14        0.021716      15  \n",
       "15        0.021543      16  \n",
       "16        0.021554      17  \n",
       "17        0.021534      18  \n",
       "18        0.021443      19  \n",
       "19        0.021431      20  \n",
       "20        0.021255      21  \n",
       "21        0.020980      22  \n",
       "22        0.020966      23  \n",
       "23        0.020835      24  \n",
       "24        0.020676      25  \n",
       "25        0.020673      26  \n",
       "26        0.020310      27  \n",
       "27        0.020300      28  \n",
       "28        0.020235      29  \n",
       "29        0.019872      30  "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 10 fold cross validation, linear regression\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",
    "    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.ols(formula=formula, data=df_train).fit()\n",
    "\n",
    "        # MSE of training set\n",
    "        train_pred = res.predict(df_train)\n",
    "        mse_train = np.mean((train_pred - df_train.y)**2)\n",
    "\n",
    "        # MSE of test set\n",
    "        test_pred = res.predict(df_test)\n",
    "        mse_test = np.mean((test_pred - df_test.y)**2)\n",
    "\n",
    "        # store results\n",
    "        train_mses.append(mse_train)\n",
    "        test_mses.append(mse_test)\n",
    "\n",
    "\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",
    "    \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",
    "    \n",
    "    \n",
    "df_res = pd.DataFrame(mean_mses)\n",
    "df_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl4FfXVwPHvyc1GFgiEsIaQSEBlDwQUFReUTUFUFIG6Votr9dVXW2mrVlpftbbVWlfcl+KCCtKKQhVEZQ8E2ZGELQGEQNgCZL3n/WMmGGNILiQ3N8v5PM88987Mb2bOIN7DzG8TVcUYY4w5WUGBDsAYY0z9ZonEGGNMtVgiMcYYUy2WSIwxxlSLJRJjjDHVYonEGGNMtVgiMcYYUy1+TSQiMkxENohIhog8UMH+MBF5392/WEQS3e2DRWSZiKxyPweVOWacu32liHwuIi39eQ/GGGMqJ/7qkCgiHuB7YDCQDSwFxqnq2jJlbgd6quqtIjIWuFxVrxaRFGCXqu4Qke7ALFVtLyLBwA6gq6ruEZG/AEdU9Y9+uQljjDFVCvbjufsDGaq6CUBE3gNGAWvLlBkF/NH9/iHwrIiIqqaXKbMGCBeRMMALCBApInuBpkBGVYG0bNlSExMTq3c3xhjTyCxbtmyPqsZVVc6fiaQ9kFVmPRs443hlVLVYRA4AscCeMmVGA+mqWgAgIrcBq4DDwEbgjqoCSUxMJC0t7SRvwxhjGicR2epLOX/WkUgF28q/R6u0jIh0A54AbnHXQ4DbgBSgHbASmFjhxUUmiEiaiKTl5OScePTGGGN84s9Ekg10KLMej1O/UWEZt/6jGZDrrscD04DrVDXTLd8bQFUz1anc+QA4q6KLq+pkVU1V1dS4uCqfzIwxxpwkfyaSpUBnEUkSkVBgLDCjXJkZwPXu9yuBOaqqIhIDfApMVNX5ZcpvB7qKSGlmGAys89sdGGOMqZLf6kjcOo87gVmAB3hNVdeIyCQgTVVnAK8Cb4tIBs6TyFj38DuBZOBBEXnQ3TbEbcX1CPC1iBQBW4Eb/HUPxpjaV1RURHZ2Nvn5+YEOpdEIDw8nPj6ekJCQkzreb81/65LU1FS1ynZj6ofNmzcTHR1NbGwsIhVVo5qapKrs3buXQ4cOkZSU9JN9IrJMVVOrOof1bDfG1Cn5+fmWRGqRiBAbG1utJ0BLJMaYOseSSO2q7p+3P/uRNGyqkLcL9mZCbiYU5UO/myHIcrMxpnGxRFIZVTiy98dksTcT9mY433M3Q2HeT8u37gaJZwcmVmNMjdi/fz9Tpkzh9ttvP+FjL774YqZMmUJMTMxxyzz00EOce+65XHTRRdUJ84RNnz6dLl260LVr1xo/tyWSShx55gwi9m04tq7igZgEJLYTJJwFsZ2gRSeIbAmTz4NtCy2RGFPP7d+/n+eff77CRFJSUoLH4znusTNnzqzy/JMmTapWfCdr+vTpjBgxwhJJbXvtyLnkFvVns7Zhi7YhS+MoPhpM6O4gosKDiQpzl/CjPBPakbitizj+XzFjTH3wwAMPkJmZSe/evRk8eDCXXHIJjzzyCG3btmXFihWsXbuWyy67jKysLPLz87n77ruZMGEC8ONwTHl5eQwfPpxzzjmHBQsW0L59ez755BOaNGnCDTfcwIgRI7jyyitJTEzk+uuv59///jdFRUVMnTqV0047jZycHMaPH8/evXvp168fn3/+OcuWLaNlyx8HOy8pKeGmm24iLS0NEeGXv/wl99xzD5mZmdxxxx3k5OQQERHByy+/TG5uLjNmzGDevHn8+c9/5qOPPqJTp0419mdmiaQSF97wEPsOF3KooJi8/GLyCpzlUH4xeQVFx7btPlTAnCOncNW2xXi8XqsnMaaGPPLvNazdcbBGz9m1XVMeHtntuPsff/xxVq9ezYoVKwD46quvWLJkCatXrz7WPPa1116jRYsWHD16lH79+jF69GhiY2N/cp6NGzfy7rvv8vLLLzNmzBg++ugjrrnmmp9dr2XLlixfvpznn3+ev/71r7zyyis88sgjDBo0iIkTJ/L5558zefLknx23YsUKtm/fzurVqwHnSQpgwoQJvPjii3Tu3JnFixdz++23M2fOHC699NJjCaymWSKpxOltm/pULq+gmD9OOpXxRXMhZz20rvlHR2NM4PTv3/8nfSyeeeYZpk2bBkBWVhYbN278WSJJSkqid+/eAPTt25ctW7ZUeO4rrrjiWJmPP/4YgG+//fbY+YcNG0bz5s1/dtwpp5zCpk2b+PWvf80ll1zCkCFDyMvLY8GCBVx11VXHyhUUFJzkXfvOEkkNiAoLZm+LPnAIyFpkicSYGlLZk0NtioyMPPb9q6++4osvvmDhwoVERERw/vnnV9gHIyws7Nh3j8fD0aNHKzx3aTmPx0NxcTHgdBKsSvPmzfnuu++YNWsWzz33HB988AFPP/00MTExx56maou9g6khbZNOZw/N0K0LAx2KMaYaoqOjOXTo0HH3HzhwgObNmxMREcH69etZtGhRjcdwzjnn8MEHHwAwe/Zs9u3b97Mye/bswev1Mnr0aP70pz+xfPlymjZtSlJSElOnTgWchPTdd9/5dF/VYYmkhqQkNGdpSReKt9T8XypjTO2JjY3l7LPPpnv37tx///0/2z9s2DCKi4vp2bMnDz74IGeeeWaNx/Dwww8ze/Zs+vTpw2effUbbtm2Jjo7+SZnt27dz/vnn07t3b2644QYee+wxAP71r3/x6quv0qtXL7p168Ynn3wCwNixY3nyySdJSUkhMzPzZ9esDhtrq4Zk5uQx5enf8mDIO3Dvemja1q/XM6ahWrduHaeffnqgwwiogoICPB4PwcHBLFy4kNtuu83vr6sq+nP3dawtqyOpIUmxkawPcetGshZBt8sDG5Axpt7atm0bY8aMwev1EhoayssvvxzokCpliaSGBAUJ4Qkp5G8LJXzbYkskxpiT1rlzZ9LT0wMdhs+sjqQG9UyII92bTMnWBYEOxRhjao0lkhrUp2MMS71dCNq1Ggryqj7AGGMaAEskNahXhxiWeU9FtAS220RaxpjGwa+JRESGicgGEckQkQcq2B8mIu+7+xeLSKK7fbCILBORVe7nIHd7tIisKLPsEZGn/XkPJ6JpeAgHWvbGi8C2xYEOxxhjaoXfEomIeIDngOFAV2CciJTv8n0TsE9Vk4GngCfc7XuAkaraA7geeBtAVQ+pau/SBWfO9o/9dQ8n47SO8WwkAd1mHRONqY9KR/89WU8//TRHjhypdhxfffUVCxbUj/pWfz6R9AcyVHWTqhYC7wGjypUZBbzpfv8QuFBERFXTVXWHu30NEC4iYWUPFJHOQCvgG7/dwUlISYhhSXFnNGsJlBQHOhxjzAmyRHLi/JlI2gNZZdaz3W0VllHVYuAAEFuuzGggXVXLjzw2Dnhfj9OjUkQmiEiaiKTl5OSc5C2cuD4JzVnqPZWgosOwe02tXdcYUzPKDiNf2rP9ySefpF+/fvTs2ZOHH34YgMOHD3PJJZfQq1cvunfvzvvvv88zzzzDjh07uOCCC7jgggsqPHfXrl3p2bMn9913HwA5OTmMHj2afv360a9fP+bPn8+WLVt48cUXeeqpp+jduzfffFOn/r38M/7sR1LRJMDlf/QrLSMi3XBedw2poNxY4NrjXVxVJwOTwenZXlWwNaVTXBTrQ903eNsWQ9tetXVpYxqezx6AH1bV7Dnb9IDhjx93d/lh5GfPns3GjRtZsmQJqsqll17K119/TU5ODu3atePTTz8FnDG4mjVrxt///nfmzp37k7lDAHJzc5k2bRrr169HRI4N+3733Xdzzz33cM4557Bt2zaGDh3KunXruPXWW4mKijqWcOoyfyaSbKBDmfV4YMdxymSLSDDQDMgFEJF4YBpwnar+ZGAYEekFBKvqMj/FftKCgoTW8cnkbG9JXNYiOGNCoEMyxlTD7NmzmT17NikpKQDk5eWxceNGBg4cyH333cdvf/tbRowYwcCBAys9T9OmTQkPD+fmm2/mkksuYcSIEQB88cUXrF279li5gwcP+m1wRX/xZyJZCnQWkSRgO84TxPhyZWbgVKYvBK4E5qiqikgM8CkwUVXnV3DuccC7fou8mlI6tmDR1mQu2brQ2lcbUx2VPDnUFlVl4sSJ3HLLLT/bt2zZMmbOnMnEiRMZMmQIDz300HHPExwczJIlS/jyyy957733ePbZZ5kzZw5er5eFCxfSpEkTf96GX/ntd86t87gTmAWsAz5Q1TUiMklELnWLvQrEikgGcC9Q2kT4TiAZeLBMU99WZU4/hjqcSPokxDj1JId2wP6sqg8wxtQZ5YdbHzp0KK+99hp5eU4n4+3bt7N792527NhBREQE11xzDffddx/Lly+v8PhSeXl5HDhwgIsvvpinn3762KuzIUOG8Oyzzx4rV7rdn8O+1zS/jrWlqjOBmeW2PVTmez5wVQXH/Rn4cyXnPaUGw6xxKR2a86T3VGdl2yKI6VD5AcaYOqPsMPLDhw/nySefZN26dQwYMACAqKgo3nnnHTIyMrj//vsJCgoiJCSEF154AXCmuh0+fDht27Zl7ty5x8576NAhRo0aRX5+PqrKU089BTizLd5xxx307NmT4uJizj33XF588UVGjhzJlVdeySeffMI///nPKl+dBZINI+8ng//6JTMOX0OT1PFwyd9q9drG1Gc2jHxgVGcYeXuF7ye9O7ZkhSZbx0RjTINnicRPUhKas6ioM+xaC/kHAh2OMcb4jSUSP+nTMYaleiqCQvbSQIdjTL3SGF651yXV/fO2ROInnVtFszHkNErwOBXuxhifhIeHs3fvXksmtURV2bt3L+Hh4Sd9Dpsh0U88QUKXDq3ZtDOJzpZIjPFZfHw82dnZ1ObQRo1deHg48fHxJ328JRI/SunQnPlbk0nOnoeUFIEnJNAhGVPnhYSEkJSUFOgwzAmwV1t+1KdjDEtKuiDFR+GHlYEOxxhj/MISiR+ldGhOWtmOicYY0wBZIvGj5pGhRLaMJye4rSUSY0yDZYnEz1ISYlha0hnNWgzWCsUY0wBZIvGzlITmzC9MRvJ2wb7NgQ7HGGNqnCUSPysdCRhwJroyxpgGxhKJn53aOprtIQkc9USDjbtljGmALJH4WbAniB7xzVkTdCpk2ROJMabhsURSC1ISmjPvaCfIWQ9HcgMdjjHG1Ci/JhIRGSYiG0QkQ0QeqGB/mIi87+5fLCKJ7vbBIrJMRFa5n4PKHBMqIpNF5HsRWS8io/15DzWhT0JzFpd0cVaylgQ2GGOMqWF+SyQi4gGeA4YDXYFxItK1XLGbgH2qmgw8BTzhbt8DjFTVHjhzur9d5pjfA7tVtYt73nn+uoeakpIQw3faiRIJtnoSY0yD488nkv5AhqpuUtVC4D1gVLkyo4A33e8fAheKiKhquqrucLevAcJFJMxd/yXwGICqelV1jx/voUa0jAqjdYsYtoZ2tnoSY0yD489E0h7IKrOe7W6rsIyqFgMHgNhyZUYD6apaICIx7rY/ichyEZkqIq1rPvSa1ychhgVFndHty6G4INDhGGNMjfFnIpEKtpXv2l1pGRHphvO66xZ3UzAQD8xX1T7AQuCvFV5cZIKIpIlIWl0YjjoloTnf5HdCSgpgx4pAh2OMMTXGn4kkG+hQZj0e2HG8MiISDDQDct31eGAacJ2qZrrl9wJH3O0AU4E+FV1cVSeraqqqpsbFxVX/bqqpT0Jz0rylFe427pYxpuHwZyJZCnQWkSQRCQXGAjPKlZmBU5kOcCUwR1XVfYX1KTBRVeeXFlZnyrR/A+e7my4E1vrvFmrOaW2jORzSnJzwjvD97ECHY4wxNcZvicSt87gTmAWsAz5Q1TUiMklELnWLvQrEikgGcC9Q2kT4TiAZeFBEVrhLK3ffb4E/ishK4Frgf/11DzUpxBNEz/YxzPQMgq3fwu51gQ7JGGNqhDSGeZFTU1M1LS0t0GHw2Mx1TJu/isVNfo30Hg8jngp0SMYYc1wiskxVU6sqZz3ba1FKQnN2l0SyN2kkfPceHN0f6JCMMabaLJHUoj4JTuvlb2Iuh6Ij8N27AY7IGGOqzxJJLWrVNJwOLZrwye44iO8PS14GrzfQYRljTLVUmkhExCMi79RWMI3B5b3bM+/7HPZ2ux5yM2HTnECHZIwx1VJpIlHVEiDObb5rasD4MzoSJMIre3tCZCvnqcQYY+qxYB/KbAHmi8gM4HDpRlX9u7+CasjaNAtnaLfWvLv8B/73rOsI/vZvkLsZWiQFOjRjjDkpvtSR7AD+45aNLrOYk3TtmYnsP1LEZ2HDIcgDaa8GOiRjjDlpVT6RqOojACIS7axqnt+jauDOPKUFXVpH8dKKo4w4fSSy/G04/3cQGhHo0Iwx5oRV+UQiIt1FJB1YDaxxJ5rq5v/QGi4R4doBiazefpDvO46D/P2w+sNAh2WMMSfFl1dbk4F7VbWjqnbEGZLEaoir6YqU9kSHBfNCZhy07g5LJkMjGGXAGNPw+JJIIlV1bumKqn4FRPotokYiMiyY0X3jmbl6F4d63gg/rLJJr4wx9ZIviWSTiDwoIonu8gdgs78DawyuObMjhSVe/nWkP4Q3g8UvBTokY4w5Yb4kkl8CccDH7tISuNGfQTUWya2iOCe5JW+m5eDtfQ2smwEHdwY6LGOMOSFV9mwHfqeqd6lqH3f5H1XdV0vxNXjXDejIzgP5fBtzKXhLYNkbgQ7JGGNOiC892/vWUiyN0oWnt6Z9TBNeXAV0HgzLXofiwkCHZYwxPvPl1Va6iMwQkWtF5IrSxe+RNRKeIOEXZyawIHMv2ztfA3m7nFdcxhhTT/iSSFrgzJU+CBjpLiN8ObmIDBORDSKSISIPVLA/TETed/cvFpFEd/tgt7/KKvdzUJljvnLPWX7mxHrr6tQOhHqCeGl7EjRPsvG3jDH1SqU92906kpWqesJT+bnHPgcMBrKBpSIyQ1XLzrF+E7BPVZNFZCzwBHA1sAcYqao7RKQ7znS97csc9wtVDfyUhzUkNiqMEb3a8lH6Dn530S8Jn/Mg7PwO2vYKdGjGGFMlX+pILq2sTCX6AxmquklVC4H3gFHlyowC3nS/fwhcKCKiqumqusPdvgYIF5Gwk4yjXrhuQCKHC0v4WM+HkAh7KjHG1Bu+vNpaICLPishAEelTuvhwXHsgq8x6Nj99qvhJGVUtBg4AseXKjAbSVbWgzLbX3ddaD4qI+BBLnde7Qwy94pvxalou2vNqWDUVjuQGOixjjKmSL4nkLKAbMAn4m7v81YfjKvqBLz8GSKVl3DG9ngBuKbP/F6raAxjoLtdWeHGRCSKSJiJpOTk5PoQbeNcOSCQz5zAr2lwFxfmQbnOKGWPqvioTiapeUMEyqKrjcJ5AOpRZj8cZkr7CMiISDDQDct31eGAacJ2qZpaJZ7v7eQiYgvMKraK4J6tqqqqmxsXF+RBu4I3o2ZbmESG8uD4MOp4DS192+pYYY0wd5svov61F5FUR+cxd7yoiN/lw7qVAZxFJcmdYHAuUb9c6A7je/X4lMEdVVURigE+Biao6v0wswSLS0v0egtN6bLUPsdQL4SEeru6XwH/X7mJvjxth/zZYMSXQYRljTKV8ebX1Bk6rqXbu+vfA/1R1kFvncad77DrgA1VdIyKTRKS0Av9VIFZEMoB7gdImwncCycCD5Zr5hgGzRGQlsALYTgMbifgXZyQA8NqebhDfH76cBPkHAxyVMcYcn2gVQ5eLyFJV7Sci6aqa4m5boaq9ayXCGpCamqppafWntfDNb6axfNs+Fl4XQ9jrF8HZ/wODHwl0WMaYRkZElqlqalXlfHkiOSwisbiV4CJyJk7rKuMn15/VkdzDhXy6ty30GgeLnofcTYEOyxhjKuRLIrkXpy6jk4jMB94Cfu3XqBq5szu15JSWkby1cCtc+DAEhcDsBwMdljHGVMiXVlvLgfNwmgHfAnRT1ZX+DqwxCwoSrh3QkRVZ+1m2LxwG3gPr/wObvw50aMYY8zO+PJGgqsWqukZVV6tqkb+DMjAmtQMxESE8PzcDBtwJzRLg84nWHNgYU+f4lEhM7YsMC+bGs5L4cv1u1uYUwZBJsGs1LH+z6oONMaYWWSKpw244K5GosGCe+yoDul4GCWfBnD/D0f2BDs0YY445biIpO65WRUttBtlYNYsI4ZozOzJz1U4y9xyG4Y874299/WSgQzPGmGMqeyIpHVfrOWAxMBmn899i4Bn/h2YAbjoniVBPEC98lekMK9/nWlj8IuzJCHRoxhgDVJJISsfVArYCfdxxq/oCKYD9itWSuOgwxvVPYHr6drL3HYFBD0JwE5j9+0CHZowxgG91JKep6qrSFVVdDdSbXu0NwYRzT0EEXpq3CaJawXn3w/efQ8YXNXMBVafD4+qP4b8Pw1uj4C+nwH/uqZnzG2MatEpnSHStE5FXgHdwerdfgzN2lqkl7WKacEVKPO+nZfHrQcm0OuNWSHsdPv8d3HYeeEJ8P5nXC7mZzgyMO9Kdz50rocAdrCAoBFp3hVZdIe01p4K/51X+uTFjTIPgSyK5EbgNuNtd/xp4wW8RmQrddn4npi7L4pVvN/O7i0+HoY/Ce+OdH/szbqn84JIiyJwLqz6ADZ9D4SFnuycMWneDHqOd+pe2vaHV6RAcBiXF8MYlzlNJfCq0SPL/TRpj6qUqB20EEJEmQIKqbvB/SDWvvg3aeDx3vZvOF+t2Mf+3g2geEeK8gtr5HdyVDhEtflpYFbKXwsoPYM3HcGQvNGkOp4+EDmc4SSPu1MqfZvZvgxfPgdhk+OWsE3vyMcbUezU2aKM75PsK4HN3vbeIlJ9XxNSCOy5I5khhCa8v2AIiMOwxKDgIXz32Y6Gc752+Js/0hlcHQ/rbkHQujH0X/vd7uPSfkHINtOledWKISYCRz8D2ZTD3Ub/emzGm/vLl1dbDOLMQfgWgqitEJNF/IZnjObVNNIO7tuaN+Zv51cAkolt3g743wtJXoUkLpwJ+5wqQIEg6D877LZw2AsKbnvxFu10Gm26Ab592ztnpghq7H2NMw+BLq61iVbVh4+uIOy9I5mB+Me8s2uZsuOD3EBYF8x531of+H9y7Dq6bDr3HVy+JlBr6GLTsAtNugbyc6p/PGNOg+JJIVovIeMAjIp1F5J/AAj/HZY6jV4cYBnZuyavfbiK/qAQiY+HmOXBnGtwyDwbcAdFtavaioRFw5WvO0CzTb3NafhljjMuXRPJroBtQAEzBmdSqyql2AURkmIhsEJEMEXmggv1hIvK+u39x6SszERksIstEZJX7OaiCY2eISIOZr/1E3HFBMnvyCnl/aZazoWUytOzs34u26e60FMv4r9Oz3hhjXJUmEhHxAI+o6u9VtZ+7/EFV86s6sXvsc8BwoCswTkS6lit2E7BPVZOBp4An3O17gJGq2gO4Hni73LmvAPKqvr2G6YykFqR2bM5L8zIpLK7Fp4N+N8OpF8N/H4IdK2rvusaYOq3SRKKqJUDfkzx3fyBDVTepaiHwHjCqXJlRQOm46B8CF4qIqGq6qu5wt68BwkUkDEBEonBmbfzzScZV74kIdwxKZseBfKanb6/NC8Oo5yAyDj66CQoabS43xpThy6utdPc10rUickXp4sNx7YGsMuvZ7rYKy6hqMc5rs9hyZUYD6apa4K7/CWcwySOVXVxEJohImoik5eQ0vAri87vE0b19U16Yl0mJt+q+QDUmogVcMRn2ZsJnv6296xpj6ixfEkkLYC8wCBjpLiN8OE4q2Fb+F6/SMiLSDed11y3uem8gWVWnVXVxVZ3sDjSZGhcX50O49YuIcMf5yWzec5hPV+2s3YsnDYRz74MV78CqD2v32saYOqfKfiSqeuNJnjsb6FBmPR7YcZwy2SISDDQDcgFEJB6YBlynqplu+QFAXxHZ4sbeSkS+UtXzTzLGem1otzYkt4ri+bkZjOjRlqCgivKyn5z3AGya5wyh0r6vDaFiTCNWZSIRkXCcSvFuQHjpdlX9ZRWHLgU6i0gSsB0YC4wvV2YGTmX6QuBKYI6qqojEAJ8CE1V1fplrvoA7zpfbwus/jTWJAAQFCbef34l7P/iOZ+dm0KV1FCGeIEKDgwj1BBHifoYFBx3bHhkaTLOIGhjqxBMMo1+BFwfCRze7Q6j40r/VGNPQ+PJ//tvAemAoMAn4BT6M/quqxSJyJzAL8ACvqeoaEZkEpKnqDOBV4G0RycB5EhnrHn4nkAw8KCIPutuGqOpu32+tcbi0Vzv+OSeDv//3e5+PaR4RQqe4KGdpFXnse3zzJgR7TmD25eYdYcTfnYr3tFerHjzSGNMgVTloo4ikq2qKiKxU1Z4iEgLMUtWf9e2oqxrKoI3Hk1dQzK6D+RQWe52lxEtRsZeCEu+xbUXu90P5xWzac5jMnDw25eSxJ6/w2HlCPUEktow4llgu79OeTnFRlV9cFd6+HLYvh18vg6iGVx9lTGPl66CNvjyRFLmf+0WkO/ADkFiN2EwNiwoLJqqqH/zj2H+kkMwcJ7Fk5uSRufswG344xOy1u/ggLYtP7xpIXHTY8U8gAsP/Ai8MgC8fgVHPnuRdGGPqK18SyWQRaQ48iFOnEQU85NeoTK2JiQilb8dQ+nZs/pPta3cc5PLn53P3e+m8fdMZeCqryI/rAmfcCgufg9Qbncp3Y0yjUeULcVV9RVX3qeo8VT1FVVupqo2R0cB1bdeUP13WnQWZe3n6Cx/qX877rTMN8Mz7bSwuYxoZX1ptVfj0oaqTaj4cU5eMSe1A2pZc/jkngz4dm3PBqa2OXzi8KQye5IwQvOJf0Ofa2gvUGBNQvjTROVxmKcEZOyvRjzGZOmTSqO6c3rYp97y/gux9lQ4mAD2vdmZf/OKPzkjBxphGwZdXW38rszwKnM/PhzoxDVR4iIcXftGHkhLljinpFBSXHL9wacX7kb3w1eO1F6QxJqBOoNPAMRHAKTUdiKm7EltG8uRVvfguaz+PflpFF6J2vaHvDbBkMuxaWyvxGWMCy5c521eJyEp3WQNsAP7h/9BMXTKsext+NTCJtxZu5ZMVVYw4fOFDTp3JZ79x+pkYYxo0X5r/lh2gsRjY5Y7UaxqZ3wxLwXAoAAAgAElEQVQ7jfRt+5n48Sq6tWtKcqvoigtGtIBBf4BP/xfWTIPuvgwWbYypr3x5tXWozHIUaCoiLUoXv0Zn6pQQTxDPju9DkxAPt76znMMFlfx7ou+N0KYHzP4DFB6uvSCNMbXOl0SyHMgBvgc2ut+XuUvDHXfEVKhNs3CeGZfCppw8fj9tFccdYifIA8OfhIPb4Zu/126Qxpha5Usi+Rxn2tuWqhqL86rrY1VNUlWrdG+Ezk5uyb2DuzB9xQ7+tXjb8Qt2HAA9xsCCZyB3U+0FaIypVb4kkn6qOrN0RVU/A87zX0imPrj9/GTOPzWOSf9ey8rsSvqMDJ4EnlD4/He1F5wxplb5kkj2iMgfRCRRRDqKyO9xZkw0jVhQkPDUmN7ERYdx2zvLyTtefUnTtnDeb+D7z+D72bUbpDGmVviSSMYBcTizFU53v4/zZ1CmfmgeGcoz43qzff9RJn9dyaurM26D2M7w+W+huKD2AjTG1ApferbnqurdqpoCpAIPqWqu/0Mz9UHfji24pGdbXvlmE7sP5VdcKDgUhj/u1JN8/WTtBmiM8TtfOiROEZGmIhIJrAE2iMj9vpxcRIaJyAYRyRCRByrYHyYi77v7F7vT5yIig0VkmdsZcpmIDCpzzOci8p2IrBGRF0XE4+vNGv+4f8ipFBZ7eebLjccvlHwR9BrvJJLFL9VecMYYv/Pl1VZXVT0IXAbMBBKAKod2dX/gn8MZ5LErME5EupYrdhOwT1WTgaeAJ9zte3BaivXAmdP97TLHjFHVXkB3nNdsV/lwD8aPEltGMv6MBN5dksWmnLzjF7z0GThthNPjfflbtRegMcavfEkkIe70upcBn6hqEeDLuBf9gQxV3aSqhcB7wKhyZUYBb7rfPwQuFBFR1XRV3eFuXwOEi0gYgJvUwOmVH+pjLMbP7rqwM+HBQfx19objF/KEwJWvOU8nM+6ClVNrL0BjjN/4kkheArYAkcDXItIROFjpEY72QFaZ9Wx+PmrwsTLusCsHgNhyZUYD6ap6rJZWRGYBu3F623/oQyzGz1pGhfGrc09h5qofSN+27/gFg8Pg6ncg8Rxn7pK1M2ovSGOMX/hS2f6MqrZX1YvV6ca8DbjAh3NXNDdr+aeHSsuISDec1123lItpKNAWCAMGUQERmSAiaSKSlpOT40O4prp+NfAUWkaF8dhn64/f4x0gpAmMe8+ZkvfDX/q/WXDOBph+BzzbDxY8C0XHaRRgjDkpJzyMvDp8GbQxG+hQZj0e2HG8MiISDDQDct31eJwmx9epamYFceTjzCFf/nVZ6f7JqpqqqqlxcXE+hGuqKzIsmLsv6sySzbnM3bC78sJhUfCLqdC6K7x/DWyaV/MBbVsM746D5/rD6o8gLBpm/x7+2QeWvQkldXTs0cIqJhAzpo45mflIfLUU6CwiSSISCozF+eEvawZOZTrAlcAcVVURiQE+BSaq6vzSwiISJSJt3e/BwMXAej/egzlBY/t1IKllJE98toESbxXVV01i4NrpENvJ+cHftqj6AXi9sOEzeHUovDbEOef5E+GeNfCrOXDdDIhuA/++C54/A1Z/XLfmmP/iEXgiEdbPrLKoMXWF3xKJ+9RyJzALWAd8oKprRGSSiFzqFnsViBWRDOBeoLSJ8J1AMvCgiKxwl1Y49TQzRGQl8B1OPcmL/roHc+JCPEHcP/RUNuw6xMfLs6s+IKKFk0yatoV/XQXbl5/chYsLIf1f8MIAeHcsHNzhzNZ4z2o4/wGIdKveTjkPbv4Sxk6BoBD48EaYfB5s/CLwc6ekvQ7f/h1CwuGD62DD54GNxxgfSaXvsksLiZyFM0/7sflLVLXetN9MTU3VtDQbqLi2qCqXPb+A3QfzmXvf+YSH+NDV58B2eH0Y5B+EGz6FNt19u1j+AVj+Nix63hlpuHV3OPtu6Ha500qsMt4SWDUV5j4K+7dBx7Phwoch4Qzfrl2TNn4BU8ZA8oVw+UvwzmjYtdppmNBlaO3HYwwgIstUNbXKclUlEhF5G+gErABKJ+xWVb2r2lHWEksktW/Rpr2MnbyIB4afxq3ndfLtoNzN8PrF4C2CG2ZCbDIc3g0HsuFAlvvpLvu3OZ9H3UEWEgfC2f/j/BBLRW04KlFcCMvecDpLHt4NXYY7P96eEOepxRPsfAYFu9uCf9wXFQctqjkI9g+r4bVh0CIRbvzMqcs5uh/evgx2rYGr/wVdhlTvGsachJpMJOtwOiXW2/4alkgC45dvLCVtSy5f/+YCYiJCfTtoz0Z4fTgU5IG32EkqZYVGQ0wHaBb/45J0PsT3rX7AhYdh0Qsw/xkoOOD7cefe79TDBJ3EIAsHd8IrFzqv1X71JTRt9+O+o/vgrctg91oY+y50vujEz29MNdRkIpkK3KWqO2squNpmiSQw1v9wkOH/+IZfDTyF3118uu8H7l7n/KA3ae4mCzdxxHSA8Gb+C7hU0VHnR7ykyElmJUVOQispcl6HHfteBKunwYp3oNOFMPoVp87HVwV5TtLM3eQ8ibTt+fMyR/fBW6Ng93oYN8XpzGlMLanJRDIX6A0sAY51ClTVS497UB1jiSRw7pv6HTO+28Hc+86nfUyTQIdT81Sd12Kf/cZpDXb1O9C2V9XHeUvgvfGwcTaMe7/yV1dHcp1kkrMBxr3rvL4zphb4mkh8abX1R5zhUf4P+FuZxZgq3Tu4CwB/n/19gCPxExFIvdF5oigphleHwIp3qz7u84nw/edw8ZNV139EtIDrPoG4Lk7yyZxTM7EbU0N86dk+r6KlNoIz9V+7mCbceFYiH6dns26nLyPr1FPxqXDL1xDfD6bfCp/+r1OJX5FFL8CSl2DAndDvZt/OH9HC6QMT29npc5M5t+ZiN6aafBlG/kwRWSoieSJSKCIlItKAfxFMTbv9/GSiw4L5y+cNvO9oVJzTJ+asX8PSV+CNS5z+LGWt/9R5GjltBAz+04mdv/TJJDbZ6Suz6asaC92Y6vDl1dazODMibgSaADe724zxSbOIEO64IJm5G3L4+vsGPu6ZJxiG/BmuesNpuvvSebDFHZxh+3L46GZo3weueBmCTqI/cGSs82TSohNMGQsLn4Nda+tW73zT6PhS2Z6mqqkislJVe7rbFqjqWbUSYQ2wyvbAyy8qYfBT89i+7yijerfnzkHJdIqLCnRY/rV7nTOOWO5mZ976tNec0Y9v/hKiWlXv3If3OB0Yty9z1pu0gMSzoeM5zsjKrbqeXKIypoyabLX1NXAR8ArwA7ATuMGdXKpesERSN+zNK2DyN5t4a8FWCopLGkdCyT8A02+H9f+BsGZw02xodVrNnX/fVtg6H7Z8C1u+cTpqgtN0uuPZTlJJPAdadQMtgYJDUJjnfBbkQeGhMt/znM/2faDToBPv2FmeKmyY6Yy43PEsuOD3ltzqmZpMJB2BXTiTSN2DM0Lv86qaUROB1gZLJHVLo0soXi98N8V5Smjfx7/X2r/NeZV2LLFsdbYHBTt9YnzVugecfZdvQ82U5/U6iXPeX2DXKoiIhSN74fSRcPlkCI04sfOZgKmxROKerAmQoKqVTH9Xd1kiqZv25BXw8tebeGthI0kogbA/y3li2b0OQiOd4VdCo5xh/MOinZECwqKd9dAo59XbmmlO7/49G5zOoAPugJRrnTKV8Xph3Scw70nYvcapxzn3fuhxFSyZDLN+B+16O3PRRLepnfs31VKTTyQjgb8CoaqaJCK9gUnWIdHUFEsodZDX63SWnP8P2LYAwmOcpspn3PLz+h1vCayd7iSQnHVOE+Vz74fuo53GB6U2fAYf3uS8dhv/HrTpUbv3ZE5YTSaSZTizEH6lqinutmMV7/WBJZL6YU9eAZO/3sRbC7eQX+Tl7ORYrurbgWHd2/g2grDxj6ylsOAfsO4/4AmF3uNgwK+hRZIzn8vXTzpPLy1PdRoVdLv8+OOO7fzOaW1WcBCufN0Go6zjajKRLFbVM0Qk3RKJqQ178gqYsngbU5dlkZV7lOjwYEb2aseY1A70im+GVLcS2JycPRmw8FlYMQVKCp3XU4d2QtzpTgLpeplvlekHd8CUq51h8oc97jzlnIiCQ04fmmYdnFdlxm9qMpG8CnyJM+nUaOAuIERVb62JQGuDJZL6yetVFm/OZWpaFjNX7yS/yEuX1lFc1bcDl6W0Jy46LNAhNk55u506j+3Loe/1cNrIE2+NVXjY6VOzYSb0nwBDH/vpa7Dy8g86Q8qsmQ4ZX0CJO+xf297Q9wbocaVT12NqVE0mkgjg98AQQHBmPPyTO2d6VccOA/4BeIBXVPXxcvvDgLeAvsBe4GpV3SIig4HHcVqKFQL3q+ocN5apOPOjlAD/VtUHqIIlkvrvUH4R/1m5k6lpWSzftp/gIOGC01oxJrUDfTs2JyosmNBga1par3hL4L8POU85yYPhytcgvOmP+/MPOPUqa6ZD5pfuU1A76DoKTrsEctY7s0ruXgMhkU4y6XsDtEupftNlA9Rwq62TDMADfA8MBrJx5nAfp6pry5S5HeipqreKyFjgclW9WkRSgF2qukNEugOzVLW9m0jOUNW57jzwXwL/p6qfVRaLJZKGJWP3IaamZfPR8u3syTs2IDWhwUFEhwUTFR5MZKjzWboeFeYske4SFeY59j0yNJjIMM+x/VFhwYQFB9krtNqS9rozNlncaXDFZPhhpZs85jhD9TeNd5JHt8ugfepPn35UnU6Zy1536muKjkCbnu5TylU/TUzmhFU7kYjIjMoOrKrVlogMAP6oqkPd9YnucY+VKTPLLbNQRIJxOjzGlZ1ES5z/m/cA7VS1oNw1/gGsVtWXK4vFEknDVFTi5duMPWzdc5i8gmLyCkrIKygiL7+YvIJiDrmfhwt+XC8o9m0okRCP0KxJCE3DQ2jaJMT53iSEZk2Cj21v1iSE+OYRnJ0ca0mnujLnwAfXO5XwAM0SoOulTsV9uz6+vTrLPwArP3CG9d+1GkIioPsV0OcGZ1BN+290wnxNJJW8lGQAkAW8CyzGea11Itq7x5fKBspPhn2sjKoWi8gBIBYncZQaDaRXkERigJE4r85MIxTiCeKCU1vBqb4fU1Ti5UhBCXmFPyaYIwUlxxLO4cIfk86Bo0UcPFrEgaNF7D9SyLbcIxxw10u8P/4D7OIebXh8dE+ahp9gxz3zo06DnKFjNsx0pk1u3+fEf/jDm0H/XznNlLcvd59SPoL0d5wWZb3HQ8+roWlb/9xDI1ZZImmD81pqHDAe+BR4V1XX+Hjuiv4WlH/8qbSMiHQDnsCpn6HM9mCcBPeMqm6q8OIiE4AJAAkJCT6GbBq6EE8QzSKCaBZx8j/6qsqRwhIOHC1ixnc7eHLWBlZt/4Znx/WhV4eYGoy2kYnr4izVJeJMvRzfF4b+n5NMvnsXvngYvnzEmWWy93joMhxCwqt/PeNzz/YwnITyJE5nxH/6cEy1Xm2JSDwwB7hRVeeXO/drQJ6q3uXLTdqrLeNPy7bu465309l9KJ8Hhp/OL89OtFddddGeDGeomhXvwqEdTifLHlc6SaXdSTwBNQI1UtnuJpBLcJJIIjADeE1Vt/sQQDBOZfuFwHacyvbxZZ9oROQOoEeZyvYrVHWM+9pqHk7S+qjcef8MnA5cpao+vfC2RGL8bf+RQu7/cCX/XbuLi05vzV+v6klMRGigwzIV8ZY4/VBWTHHGBCvOd/rC9B4Pfa51et4boGYq298EugOfAe+p6uqTCOJi4Gmc5r+vqeqjIjIJSFPVGSISDrwNpAC5wFhV3SQifwAm4syBUmoITnPgLGA9P84f/6yqvlJZHJZITG1QVV6fv4XHPltHXFQY/xyfQt+OLQIdlqnM0f3O2GIrpkD2EohJgKvfgbb1ZnBzv6qJROIFDrurZQsJoKpab9rVWSIxtWll9n7unJLO9v1HuW/Iqdxy7ikEBdlrkzovaylMvd4ZqXjkM9Dr6kBHFHC+JpLjtqlT1SBVjXaXpmWW6PqURIypbT3jY/jPXecwrFsbnvh8PTe8sfQn/V1MHdWhH0yYB/H9YNoEmPkbKCkKdFT1gnUFNsYPmoaH8Oz4FB69vDuLNu3l4n98w39W7qCoxKbErdOi4uDa6TDgTljyErx5KRzaFeio6jy/9WyvS+zVlgmktTsOctd76WTszqN10zDG9+/IuP4daNXUmp7Waas+hE/uhCYxMOYt6NA/0BHVuoAPkVKXWCIxgVbiVeas381bC7fwzcY9BAcJw7q34boBifRLbG7NheuqH1bD+7+AA9vh4r9A3xsbVTNhSyRlWCIxdcmmnDzeWeQMk38ov5jT2kRz7YCOXNa7PZFhlfURNgFxdB989CvI+C+kXAMX/63RdGS0RFKGJRJTFx0pLGbGih28tXAra3ceJDosmNF947k8pT2hwUEUFnspLPFSVOylwP0sLPFSVOJ19ykhQUJEWDCRoR4i3MEnI0KDiQj1EBkaTESYhxCPVYVWm7cEvnocvv6LM7rwmLchpkOgo/I7SyRlWCIxdZmqsnzbPt5auJWZq3ZSVFKz/0+GeoKIiQihZ3wzUhKak5IQQ6/4GHv6ORnrP4WPbwH1QtJASDoPks6FVl1PfE6WesASSRmWSEx9sSevgEWb9hIcJIR4gggNDjr2GVpuPSRIKPYqRwqLOVxQwuHCYo4WlnC4sIQjBcUcLizhaKHzuftgASuy9pGZ43QNCxLo0jqaPh2bk9Ihhj4dm5MUG2n9XXxROlPk5q8hN9PZFtHyp4mlxSkNoi7FEkkZlkiMcRw4UkR61j7St+0nPWs/6dv2cSi/GIBmTULo3SGGdjHhhHqCCAvxOJ/BQYSFBP10W0gQsZFh9E9qgacxJ58D2U5C2TQPNs9zph4GZxrgpHOdpWk7kCAQj/sZ5Dy9lN8m4rxC8xY7n1r63V0v3YdC4jm1MiOkJZIyLJEYUzGvV9m0J4/l25ykkr5tP7mHCykoduphCopL8FbyE9GhRROuPbMjY1I72NhiqrA3w0kom+bBlm+cinp/aNMTbviPM3S+H1kiKcMSiTEnr7jES0Gx9yfJpaDYy/e7DvHWwq0s2ZxLeEgQl/Vuz/VnJXJ6Wxv4AgCvF3avhfz97hOG13nKUHW+H9tWZntQsLt43MVdF8+P2/dmwvTbnH4t13wEIU38dguWSMqwRGKM/6zdcZC3F21hWvp28ou89E9swfVnJTKkW2trMeYvqz6Ej26GLsOcQSY9/mk4YYmkDEskxvjf/iOFTE3L5q1FW8jKPUqbpuH84owExvZPIC46LNDhNTxLX3Hmuu81HkY955dWY5ZIyrBEYkztKfEqc9fv5s2FTi/+UE8Q1w7oyH1DTqVJqCfQ4TUs8/4Ccx+FM++AoY/WeEuxmpiz3RhjTpgnSLioa2su6tqajN15vPLNJl79djNfrtvFX67sRf8km6Olxpx7vzPs/aLnIKIFnHtfQMKwF5jGGL9JbhXF46N78u6vzqRElasnL+SPM9ZwpLA40KE1DCIw9DHoeTXM+ROkvRaQMCyRGGP8bkCnWGb9z7lcPyCRNxZsYdjT37Bo095Ah9UwBAU5dSSdh8J/7nVmfKztEPx5chEZJiIbRCRDRB6oYH+YiLzv7l8sIonu9sEiskxEVrmfg8oc86iIZIlInj9jN8bUrIjQYP54aTfem3AmIjB28iIe+mQ1hwvs6aTaPCFw1RuQcKY7wOSXtXp5vyUSEfEAzwHDga7AOBHpWq7YTcA+VU0GngKecLfvAUaqag/gepx53Uv9G2h8EwMY00CceUosn909kBvPTuTtRVsZ+vTXLMjcE+iw6r/QCBj3HsSdBu9f40wdXEv8+UTSH8hQ1U2qWgi8B4wqV2YU8Kb7/UPgQhERVU1X1R3u9jVAuIiEAajqIlXd6ce4jTF+FhEazMMju/HBLQMI8QQx/uXF/GH6KvLs6aR6msQ4nRSjWsOUq2D3ulq5rD9bbbUHssqsZwNnHK+MqhaLyAEgFueJpNRoIF1VT2jSaxGZAEwASEhIOLHIjTG1ol9iC2beNZC/zd7Aq/M388mKHcRGhhIe4nHH+HI+j60HewgPcdZPb9uUi05vZUOzlBfdGq6bDq8Ohbcvh1/NhaZt/XpJfyaSiho0l++0UmkZEemG87pryIleXFUnA5PB6UdyoscbY2pHk1APfxjRleE92jA1LZsjhSUUFJeQX+QMx5JXUMyevEJnaBZ325FCZ/EECWee0oKh3dowpGsb2jRrHBNOVal5Ilw7DZa9AZFxfr+cPxNJNlB25pd4YMdxymSLSDDQDMgFEJF4YBpwnapm+jFOY0wd0LdjC/p29K2PiaqyMvsAs9b8wOdrfuChT9bw0Cdr6N0hhqHd2jC0W2tOiYvyc8R1XOuuzvTAtcBvPdvdxPA9cCGwHVgKjFfVNWXK3AH0UNVbRWQscIWqjhGRGGAeMElVPzrO+fNU1ae/Kdaz3ZiGLWP3IWat2cWsNT+wMvsAAJ1bRTG0WxuGdW9Dt3ZNkQYwP0htqxNDpIjIxcDTgAd4TVUfFZFJQJqqzhCRcJwWWSk4TyJjVXWTiPwBmAhsLHO6Iaq6W0T+AowH2uE84byiqn+sLA5LJMY0Htv3H2X2mh+YteYHlmzOxatwWptoxqR24LKU9rSItDoVX9WJRFJXWCIxpnHKPVzIp6t28mFaFt9lHyDEI1x0emvGpHZgYOeWBNvoxJWyRFKGJRJjzPofDjI1LZtp6dvJPVxI66ZhjO4Tz1WpHUhqGRno8OokSyRlWCIxxpQqLPYyZ/1upqZlMXfDbrwK/RKbc1VqB05rE02QCCIQJOIuIO5n6bZmESE0axIS6FvxO0skZVgiMcZUZPfBfD5O384HaVlsyjns83GhwUE8NaY3l/T0b/+MQLNEUoYlEmNMZUqbE+/JK8Cr4FVFVY999yruuuL1wrtLtpG2dR+/v/h0bh6Y1GBbhNl8JMYY4yMRoVeHGJ/LX9KzLfd+sIJHZ65j+/6jPDiiK56ghplMfGFNFowx5gSFh3h4dlwfbj4niTcWbOG2d5ZxtLAk0GEFjCUSY4w5CUFBwh9GdOXhkV3577pdjHt5EXvzTmhIwAbDEokxxlTDjWcn8cIv+rJu50GueGEBm/f4XmnfUFgiMcaYahrWvQ3vTjiTQ/nFjH5hAcu27gt0SLXKEokxxtSAPgnN+fi2s2gaHsz4lxfx+eofAh1SrbFEYowxNSSxZSQf3XYWXds15bZ/LeP1+ZsDHVKtsH4kxhhTw44WlnD3e+nMXruLllGhhHqCCAkOIjhICPEEERocRIgniBCPu+4JIiIsmKSWkXSKiyS5VRSntIyiSagnoPdh/UiMMSZAmoR6eOGavrw+fzOb9hymqNhLsVcpLPFSVOylqMRLUYmznldQTHGJcmB3Hp+u3IG3zL/t28c0oVOrKDrFRdIpLspZWkXSKrpuTeBlicQYY/zAEyTcPPCUEzomv6iErXuPkJmTR+buPOcz5zDvb8nlSJl+KqN6t+OJ0T0JDwnsE0spSyTGGFNHhId4OLVNNKe2if7JdlXlh4P5ZO4+zLcZe3hxXiY7D+Qz+dq+dWLOeqtsN8aYOk5EaNusCed0bskDw0/jmXEprNi2n9EvLCAr90igw/NvIhGRYSKyQUQyROSBCvaHicj77v7FIpLobh8sIstEZJX7OajMMX3d7Rki8ow01NHSjDHmOC7t1Y63b+pPzqECLn9+Aavc6YUDxW+JREQ8wHPAcKArME5EupYrdhOwT1WTgaeAJ9zte4CRqtoDuB5nOt5SLwATgM7uMsxf92CMMXXVGafE8vHtZxEWHMSYlxYyZ/2ugMXizyeS/kCGqm5S1ULgPWBUuTKjgDfd7x8CF4qIqGq6qu5wt68Bwt2nl7ZAU1VdqE675beAy/x4D8YYU2clt4pm2h1n0alVJDe/mcaUxdsCEoc/E0l7IKvMera7rcIyqloMHABiy5UZDaSraoFbPruKcxpjTKPRKjqc9ycM4Lwucfxu2iqenLWe2u4f6M9EUlHdRfm7q7SMiHTDed11ywmcs/TYCSKSJiJpOTk5PoRrjDH1U2RYMC9fl8q4/h14bm4m97y/gsJib61d35+JJBvoUGY9HthxvDIiEgw0A3Ld9XhgGnCdqmaWKR9fxTkBUNXJqpqqqqlxcXHVvBVjjKnbgj1B/N/lPbh/6KlMX7GD619bwoGjRbVybX8mkqVAZxFJEpFQYCwwo1yZGTiV6QBXAnNUVUUkBvgUmKiq80sLq+pO4JCInOm21roO+MSP92CMMfWGiHDHBck8dXUv0rbmctWLC9h9KN/v1/VbInHrPO4EZgHrgA9UdY2ITBKRS91irwKxIpIB3AuUNhG+E0gGHhSRFe7Syt13G/AKkAFkAp/56x6MMaY+ujwlnjdv7E9ibCQxTfzfYdEGbTTGGFMhXwdttJ7txpj/b+/uY66s6ziOvz/hAwgWgmlOMdRsyxiDNO3BDF06tRY4K3KytJppRWmuKatVJGtzjbKtP2CaLB9IxATEbAWNm6Q5HuTBmwezMKkMBjaTwJgofPvj97vlcHbO/eB1bs59nfN5bYxz/8718Ptev/s+33Nd1znfn1khTiRmZlaIE4mZmRXiRGJmZoU4kZiZWSFOJGZmVogTiZmZFeJEYmZmhbTFFxIlvQT8vaLpRNKcJ62k1WJqtXig9WJqtXig9WIqGs+7I6LHYoVtkUiqSXq6N9/WLJNWi6nV4oHWi6nV4oHWi+lIxeNLW2ZmVogTiZmZFdKuieTuZnegH7RaTK0WD7ReTK0WD7ReTEcknra8R2JmZo3TrmckZmbWIG2XSCRdLuk5SVslTet5jYFP0jZJG/MEYKWbeEXSHEm7JG2qaBshaamkv+b/T2hmH/uqTkzTJf2rYrK2K5vZx76QNEpSh6RnJW2WdHNuL+U4dRNPmcdosKTVkp7JMf0wt58haVUeo4fzjLWN3Xc7XdqSNNUp0e4AAAa2SURBVAj4C3Apaf73NcA1EbGlqR0rSNI24LyIKOXn3yVdBOwF7o+IMbntx8DLEXFnTvgnRMTtzexnX9SJaTqwNyJmNrNvb4WkU4BTImKdpOOBtcAk4HpKOE7dxPM5yjtGAoZGxF5JRwN/Am4mzT67ICLmSZoNPBMRsxq573Y7Izkf2BoRf4uI/cA8YGKT+9T2IuJJ4OWq5onAffnxfaQ/8tKoE1NpRcSOiFiXH+8hTZ99KiUdp27iKa1I9uYfj87/ArgE+HVu75cxardEcirwz4qfX6TkvzxZAEskrZX0lWZ3pkFOjogdkP7ogZOa3J9GmSqpM1/6KsVloGqSRgPjgVW0wDhVxQMlHiNJgyRtAHYBS4HngVci4o28SL+85rVbIlGNtla4tvfRiPgAcAXw9XxZxQaeWcBZwDhgB/CT5nan7yQNAx4FbomI/za7P0XViKfUYxQRByJiHHAa6QrM+2ot1uj9tlsieREYVfHzacD2JvWlYSJie/5/F7CQ9AtUdjvzdeyu69m7mtyfwiJiZ/5DPwjcQ8nGKV93fxSYGxELcnNpx6lWPGUfoy4R8QqwHPgQMFzSUfmpfnnNa7dEsgY4O3+K4Rjg88DiJvepEElD881CJA0FLgM2db9WKSwGrsuPrwMea2JfGqLrBTe7ihKNU76Rey/wbET8tOKpUo5TvXhKPkbvlDQ8Px4CfIJ076cD+ExerF/GqK0+tQWQP873M2AQMCciftTkLhUi6UzSWQjAUcCvyhaTpIeACaRKpTuBHwCLgPnA6cA/gM9GRGluXteJaQLpkkkA24Abu+4vDHSSLgRWABuBg7n5O6T7CqUbp27iuYbyjtFY0s30QaSThPkRcUd+jZgHjADWA1Mi4rWG7rvdEomZmTVWu13aMjOzBnMiMTOzQpxIzMysECcSMzMrxInEzMwKcSKxtiRpuaR+n8ta0jdzhdm5Ve3jGllZVtKBXK12k6RHJB2X298laZ6k5yVtkfRbSe9t1H576NMROcbWfE4kZn1U8S3h3vgacGVEXFvVPg5oZInyfRExLlca3g/clL90txBYHhFnRcQ5pO9KnNzA/Zo5kdjAJWl0fjd/T55fYUn+xu5h73YlnZhL6SPpekmLJD0u6QVJUyXdKmm9pJWSRlTsYoqkp/K7+PPz+kNzsb41eZ2JFdt9RNLjwJIafb01b2eTpFty22zgTGCxpG9VLHsMcAcwOZ9FTFaa12NRLha4Mn+5rGt+jAckLVOaT+KGXhy6FcB7gIuB1yNidtcTEbEhIlbUOM6V86Z8W6nkfdcZ1Zbcr3k9HKMh+eynU9LDwJBe9NVaQF/eWZk1w9mkOWNukDQfuBp4sId1xpCquQ4GtgK3R8R4SXcBXyBVNoA0d8NHcpHLOXm97wLLIuJLudzEakl/yMt/GBhb/c1tSecCXwQuIBUGXSXpjxFxk6TLgYsr54qJiP2Svk+aQ2Zq3sbPgfURMUnSJcD9pLMWgLGkmklDgfWSnuiqr1Ytny1dAfwux7O2h2PVk2nAGRHxWlf5DeofoxuB/0XE2JwI1xXct5WEz0hsoHshIjbkx2uB0b1YpyMi9kTES8Bu4PHcvrFq/YfgzblD3p5fFC8DpimV4l5OSkan5+WX1in/cSGwMCJezfNBLAA+1rvwDtvGA7k/y4CRkt6Rn3ssIvblZNRB7UKCQ3KfnyaVKrm3j/uvpxOYK2kK0FWKvN4xuoic5COiM69rbcBnJDbQVdYEOsChyyVvcOiN0OBu1jlY8fNBDv+dr64PFKQziqsj4rnKJyRdALxap4+1pifoq+6mOKjVz2r7cvnwQxuUNnOoWF93Ko8lHH48P0lKEJ8Gvifp/dQ/RvX6Zi3OZyRWVtuAc/Pj3rxY1jIZ3izgtzsidgO/B76Rb1QjaXwvtvMkMEnScUoVmK8i3afozh7g+KptXJv3OQH4d8V8HxOV5uMeSSr8uKYXfQJYBhxbeV9F0gclfbxquZ3ASZJGSjoW+FRe9m3AqIjoAG4DhgPDqH+MKmMYQ7okZ23AicTKaibwVUlPkSrsvhX/yevPBr6c22aQpijtzDegZ/S0kTxl6y+B1aRquL+IiPU9rNYBnNN1sx2YDpwnqRO4k0Ol2cnbfQJYCcyod3+kRr+ClNQuzR//3Zz3s71quddJN/9XAb8B/pyfGgQ8KGkjqWrsXXmei3rHaBYwLMdwW+63tQFX/zUbwPKnp/ZGxMxm98WsHp+RmJlZIT4jMTOzQnxGYmZmhTiRmJlZIU4kZmZWiBOJmZkV4kRiZmaFOJGYmVkh/weLzPRDblaw7QAAAABJRU5ErkJggg==\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:8: 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",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "formula           y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...\n",
       "mean_test_mse                                             0.0209818\n",
       "mean_train_mse                                            0.0198717\n",
       "num_pc                                                           30\n",
       "Name: 29, dtype: object"
      ]
     },
     "execution_count": 36,
     "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_xlabel('number of top PC used')\n",
    "ax.set_ylabel('Mean squared error')\n",
    "plt.show()\n",
    "best_model = df_res.loc[df_res.mean_test_mse.argmin()]\n",
    "best_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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_mse</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.028483</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.028483</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.026625</td>\n",
       "      <td>0.026625</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.026006</td>\n",
       "      <td>0.026161</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.027864</td>\n",
       "      <td>0.026625</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.027864</td>\n",
       "      <td>0.026316</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.027245</td>\n",
       "      <td>0.026625</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.026006</td>\n",
       "      <td>0.023375</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.028483</td>\n",
       "      <td>0.024458</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.025387</td>\n",
       "      <td>0.024458</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.026625</td>\n",
       "      <td>0.022755</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.027864</td>\n",
       "      <td>0.023065</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.026625</td>\n",
       "      <td>0.023065</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.026006</td>\n",
       "      <td>0.022601</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.024768</td>\n",
       "      <td>0.022446</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.022910</td>\n",
       "      <td>0.022755</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.026006</td>\n",
       "      <td>0.021517</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.023529</td>\n",
       "      <td>0.020898</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.024768</td>\n",
       "      <td>0.020898</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.026006</td>\n",
       "      <td>0.020898</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.027245</td>\n",
       "      <td>0.019659</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.024149</td>\n",
       "      <td>0.019505</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.026006</td>\n",
       "      <td>0.019969</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.027245</td>\n",
       "      <td>0.018885</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.024768</td>\n",
       "      <td>0.018576</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.027245</td>\n",
       "      <td>0.017957</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.027245</td>\n",
       "      <td>0.018731</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.026006</td>\n",
       "      <td>0.017957</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.029102</td>\n",
       "      <td>0.017183</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.024768</td>\n",
       "      <td>0.014861</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              formula  mean_test_mse  \\\n",
       "0                                         y ~ 1 + PC1       0.028483   \n",
       "1                                   y ~ 1 + PC1 + PC2       0.028483   \n",
       "2                             y ~ 1 + PC1 + PC2 + PC3       0.026625   \n",
       "3                       y ~ 1 + PC1 + PC2 + PC3 + PC4       0.026006   \n",
       "4                 y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5       0.027864   \n",
       "5           y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6       0.027864   \n",
       "6     y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7       0.027245   \n",
       "7   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.026006   \n",
       "8   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.028483   \n",
       "9   y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.025387   \n",
       "10  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.026625   \n",
       "11  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.027864   \n",
       "12  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.026625   \n",
       "13  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.026006   \n",
       "14  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.024768   \n",
       "15  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.022910   \n",
       "16  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.026006   \n",
       "17  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.023529   \n",
       "18  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.024768   \n",
       "19  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.026006   \n",
       "20  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.027245   \n",
       "21  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.024149   \n",
       "22  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.026006   \n",
       "23  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.027245   \n",
       "24  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.024768   \n",
       "25  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.027245   \n",
       "26  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.027245   \n",
       "27  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.026006   \n",
       "28  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.029102   \n",
       "29  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.024768   \n",
       "\n",
       "    mean_train_mse  num_pc  \n",
       "0         0.028483       1  \n",
       "1         0.028483       2  \n",
       "2         0.026625       3  \n",
       "3         0.026161       4  \n",
       "4         0.026625       5  \n",
       "5         0.026316       6  \n",
       "6         0.026625       7  \n",
       "7         0.023375       8  \n",
       "8         0.024458       9  \n",
       "9         0.024458      10  \n",
       "10        0.022755      11  \n",
       "11        0.023065      12  \n",
       "12        0.023065      13  \n",
       "13        0.022601      14  \n",
       "14        0.022446      15  \n",
       "15        0.022755      16  \n",
       "16        0.021517      17  \n",
       "17        0.020898      18  \n",
       "18        0.020898      19  \n",
       "19        0.020898      20  \n",
       "20        0.019659      21  \n",
       "21        0.019505      22  \n",
       "22        0.019969      23  \n",
       "23        0.018885      24  \n",
       "24        0.018576      25  \n",
       "25        0.017957      26  \n",
       "26        0.018731      27  \n",
       "27        0.017957      28  \n",
       "28        0.017183      29  \n",
       "29        0.014861      30  "
      ]
     },
     "execution_count": 33,
     "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=5, 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",
    "    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",
    "        # store results\n",
    "        train_mses.append(mse_train)\n",
    "        test_mses.append(mse_test)\n",
    "\n",
    "\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",
    "    \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",
    "    \n",
    "    \n",
    "df_res = pd.DataFrame(mean_mses)\n",
    "df_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXdYVFfawH8vXRCkCKKigF0EG9g1MRoTMcVsTExZN6aaTTbZzaZs4pfedtN7NdFNNs0Y00ysSYwxaizYK1ZQ7AVRUaSd748zo4gD3IEpIOf3PPcZ5t5zz30HmPve81ZRSmEwGAwGQ3Xx8bYABoPBYKjbGEViMBgMhhphFInBYDAYaoRRJAaDwWCoEUaRGAwGg6FGGEViMBgMhhphFInBYDAYaoRRJAaDwWCoEUaRGAwGg6FG+HlbAE/QuHFjlZCQ4G0xDAaDoU6xdOnSA0qp6KrG1QtFkpCQQEZGhrfFMBgMhjqFiGRbGWdMWwaDwWCoEUaRGAwGg6FGGEViMBgMhhphFInBYDAYaoRRJAaDwWCoEUaRGAwGg6FGGEViMBgMhhphFInBYDDUJTZMhcPbvS3FGRhFYjAYDHWFE4fhy1Ew71VvS3IGRpEYDAZDXWHHIlClsHedtyU5A6NIDAaDoa6QNU+/7lsPSnlXljIYRWIwGAx1hewF+vVkHhzZ6V1ZymAUicFgMNQFTh6DXcshvp9+X4vMW0aRGAwGQ10gZzGoEuhxi36/zygSg8FgMDhD1nwQX2h7EYQ1N4rEYDAYDE6SPR+adYXAUIhJqlWmrXrR2Kq6rHu2H02Kdlgae1xCaHDbDBo3i3ezVHWYb8ZA81Todbu3JTEY6hZFJ2Dn0tPfnZiOsO03KCkGX+/fxr0vQS0mr2k/jubvrXKcqBJ6HvqBP2a8TeObX/CAZHWQvBxY9aV+ijKKxGBwjpwMKCmE+P76fZNO+v2hLRDd3ruyYRRJpfRxQimsem4QrbZPprjoGfz8A9woVR0lc7p+3bsGTuRCgwjvymMw1CWyFwACLXvr9zFJ+nXv2lqhSIyPxEUUd7uJJhxk9ZyvvC1K7SRzOvgGAAq2L/S2NAZD3SJ7HsQmQ4Nw/b5xO+1437feu3LZcKsiEZGhIpIpIptF5CEHxwNF5Evb8UUikmDbP0RElorIatvroDLnXGfbv0pEZohIY3d+Bqt0HnQN+4jEd+l/vS1K7aPgCGybC91Hg2/g6excg8FQNcWFsGPJabMWgH8QRLWuNZFbblMkIuILvA2kA0nAdSKSVG7YLUCuUqoN8CrwvG3/AeAypVQKMBr4xDanH/A6cIFSqjOwCrjLXZ/BGfz8A9jSYgTJJzLYubV2PCXUGrbMhtIi6PQniEs7nZ1rMBiqZtdyKD4B8X3P3B+TpE1btQB3rkh6ApuVUluVUoXARGB4uTHDgY9tP08GBouIKKWWK6V22favBYJEJBAQ2xYiIgKEAbuoJbS6+A4UsP3nd7wtSu0iczo0iIQWvfSXYfdKOHnU21IZDHWD7Pn61ZEiyc2CwnyPi1QedyqS5kDZ2Nkc2z6HY5RSxUAeEFVuzAhguVLqpFKqCLgDWI1WIEnAeNeLXj2axLVmVUgf2u/6jsKTBd4Wp3ZQUgybZkK7i3WYYnw/nZ27Y5G3JTMY6gbZ8yG6A4SUs+I3SQIU7N/gFbHK4s6oLXGwr3y5ykrHiEgntLnrItt7f7Qi6QZsBd4ExgLPnHVxkTHAGICWLVs6L3018e15K5Fzbibj509Ju+RWz1x0/0Y4kGltbFA4JPQHcfSrdwM7Fuoorfbp+n2LnuDjp7N021zoGRkASkvgwEYdf28w1BVKimH7Iug88uxjpyK31un8LC/iTkWSA7Qo8z6Os81Q9jE5Nv9HI+AQgIjEAd8CNyilttjGdwWwvxeRScBZTnzbmHHAOIC0tDSP1VtOHnAFu34bS4OVH4EnFEnxSfhvOhw/YP2cztfCZa9rh527sUdrtbbFSwSEQLNunveTzHoUFr4N138F7S7y7LUNhuqyZxUUHj3brAUQkQB+DWpF5JY7FckSoK2IJAI7gWuB68uNmYJ2pv8BXAXMVkopEQkHpgJjlVLzy4zfCSSJSLRSaj8wBPD+b7EMPr6+ZCdeQ5+tb5C9YRnxHbq794Lrf9BKZPg70LSzhfE/wm/PwcFNcM1nENbUfbIppduCJp6nyzrYie8Lf7wDhcchINh917ezZzUseg8QmP4AJA4A/wbuv67BUFPsD1z2ir9l8fGFmA6wz/sOd7f5SGw+j7uAmeib/SSl1FoReUpELrcNGw9Eichm4F5Ory7uAtoAj4rICtsWY3PAPwnMFZFV6BXKv931GapLu4tvp1D5snv2u+6/2JLxEJEIXa6D2JSqtwvGwjWfwr4NMG4g5Cx1n2wHNkLuttNmLTvx/XUU184M913bTmkp/HivToAc+T/tnPz9Ffdf12BwBdnzIbJVxQ98MZ1qRc0tt+aRKKWmKaXaKaVaK6Wete17TCk1xfZzgVLqaqVUG6VUT6XUVtv+Z5RSIUqprmW2fbZj7ymlOiqlOiulLlNKHXTnZ6gOUU3iWB12Pkn7pnIi343RSfvWw/YFkHYT+Djxp+x4GdwyC/wCtFls5ZfukS9zmn5tV06RtOwF4qP9JO5mxae6/PZFT0PS5ZAyEua/Bgc2u//aBkNNKC3VKxJHqxE7MR0hfx/kO2HadgMms91NNOg7hjDyWT3rI/ddJOO/2v/Q9c/OnxubDLfNgbge8O0Y7UMoLXGtfJnToWlXaFQuWC+okV4dZbtZkeQfhJ8eg5Z99YoN4KJntF152n21qlWpwXAW+9ZBweHKFUmTpNNjvYhRJG6iY6+LyfZpQfja/7nnAoX5sHIiJA0/OyzQKiFRcMN30ONWWPAGfH4NFOS5Rr5j+2DHYmg/zPHx+H6Qs0QHC7iLnx/X+SqXvHw6Si20CQx+FLbOgbXfuu/aBkNNsftHEipbkXTSr142bxlF4ibEx4fdba+jXfFGNq90Q0mQNd/ovs1pN9dsHl9/faO95BXY+it8MNg1Zp+NMwF1tn/ETnw/KC6Anctqfi1HbF8Eyz+B3neefmqzk3YzNO0CM8bq8i0GQ20kex40agHhlaQvNIzRyb5mRXLu0nHo7ZxQARz87X3XT54xAaI7Qss+rpmvxy1ww/dw4hB8MAg2/1yz+TKnQ1icNmE5wi63O8xbJcUw9V59/fMfPPu4jy9c+ioc2wtz/uP66xsMNUWpqv0joFfaTTp5XZGYMvJupFFEYxZHXEjKwZkczTtEaKNI10y8aznsWgbpL7o2sTChP9z2K0y8Hj67Gi5/C7pVw/9SdELX1+o2qmL5QqJ0QlX2fOD+Gol9Fovf1+Xqr/kUAhs6HtM8Va9MFr2n/SdWQqcNzqMU/PKkDgFvPajq8a5m5zK9Mh32kn6AcAW7lsPsZ3WFhqoQH+j3D/35neHAJsjf7zh/pDwxSbDiM+2cdyboxoWYFYmbCT/vrwTLSdbN+NB1k2ZMAP9g6HKN6+a0ExEPN8/UT0LT/wV5O52fY+tvushchwr8I3bi+2oTVElx9WR1xJFd8Ou/dV/rDpdWPnbwo9osMPU+/SU0uJ41X8O8V3XekKcpLoRvb9ffF1eW5Fn4nq5gffJo1due1TD5Fjhx2Llr2FfqCf0rHwc6cqvwGORZ6+bqDowicTNtuw5gs29rYjI/Q7niZlWQB6snQ/IIHf3kDgIbwuVvQmkxzBzr/PmZ0yAg9Myy146I7wdF+bqIo6uY+X9a7vQXql6tNYjQUVw5i/VTq8G1FOTpvwfoHjSujgqsij/e1LlMyOlQ9Jpirx2XNBxu/bnq7fovdcLwr886d53s+dCwic4hqYomNoe7F81bRpG4GfHx4VDHUSSWZpGZ8UvNJ1w1CYqO19zJXhWRiXDe/bDue9jkhL+ktBQ2zoC2F+o8lcqw239d5SfZ/IuOxBpwn5bfCl2u1XL8/LgOFza4jl//raP3+tyly3zsWeW5a+dmw28v6lVp60GnO3TWlB2LzqwdVxXNuumoyCUfapOYFZTSOVbx/ayZrqM76FcvlpQ3isQDdLr4Zo6pBhydP65mEymll+nNukFzN5deAej7d4hqA9Pu134PK+xapp3YFYX9liW0iZ7fFYqkqEDLGdla26StIqKj1k4e1crE4Bp2rYDF4/RNtM/f9D5P1leb8ZD+2w59Tt/0D27WfoeakjlN5261GWz9nAsehuDGusKClVVZbhYc3WXNPwIQFAaNWnq15pZRJB4gJDSctdHpdD78K4cP7Kn+RDsW6eWru1cjdvwC9U02dxvMe83aOZnTdAtQq5V94/tB9h81N3vMfx0ObYVLXtJyO0NMR32zW/6JaQPsCkpLddRccBQMegTCmukyPp6oZACwYZr+Pxz4EIS3OL16qKl5Syk9R8KAM2vHVUWDcLj4Wf2QtfSjqsefyh+x4B+x0yTJmLbqAzEX/JVAKWLDjBqEAi8ZD4Fh2j/iKVoNhOSrtMP04JaqRmsTQnxfCLYYoRbfT+fD1GRZfnAL/P4ydLqy+pFB5/1Lhwv/eK9rnf/1kWUfw86lcNGzp3uMJ/TT5XzcHdRQmA/TH9Tmnt536n2N4iC2c83NWwc26YcVq2atsqRcrRXQL0/Csf2Vj82er5Ww3WRlhZiO2h9UXOi8bC7AKBIPkdipF+v9k2i+ZSKlJdV4+s4/COu+0zb9gBDXC1gZFz+rl/PT7q+8rMihbfqpyIpZy449a7e6Zg+lYNoDWr6La1C/M7AhpD+vK6kudkPeT33h2H74+Ql90yzbQyO+n/Yt7Hez+WXuS5C3XSfY+vqf3t9+mF7R16QmVebU03M5i92EWnhcl+2pjKx5Os/KmdD+mE46yOSgd2rIGUXiQfJTbqCF2sW6BT86f/KKz6Ck0KFZq7C4lKIS1z7pFRSVUGyfMzRWmyi2zNbKrCI2ztCv7Ydav1CjOJ25m13N7P/1U2DLLzDo4ZqXxO9wCbS9WDuJqxP2bNB+psJjZ5algTKBFW70k+zPhAVvQpfrzy4r0j4dVKmt4kI1yZyuKyKUrx1nlej20PduWPl5xWa+vBw4nO2cWQu8XnPLKBIPkjzkBnIJJWT2w/yR4UT59tJSWPpfXXywTIe/0lLFZ4uySXvmJ3r9+xf+M2092Qer379ZKcXy7bncN2klXZ6cxVM/lvmn7HGrNg/MGFtxv/XMaTrb3krIYlni++sbjLNFFE8ehekP6ez5Hrc5d64jRPSqpLRYJ2UeKd+HzVAp2Qv0A0/fu/VNsyzhLbXpMMsN5YJA/+9MvU/3txny1NnHm3aBsObV95Mc21957TirnPeAdoxPvRdKis4+fqr/iEVHu52otrrzqJcit4wi8SBBDULYd9G7NFaHaP/D5bwybjy78yxEQ237Tdtmy6xG1uzM48p3F/Dwt2vo2DSMHgkRfDhvG+e/OIcbJixm1to9p1cUVZB/spjPF23n0jfn8ad3FjBjzW4iQwKYk1nGluvrp8uKHN0DvzooK3IiVz9lVcd+HN8Xjh/UT5TOMOc5Hd1yyStaPlcQmQhXf6xNBOMGwo4lrpn3XKekSN/IG7XUN8vyiOhVQnUeGKyw+ivI+h0GPw4Nox1fv326XlUXFTg//6YqasdZJSAYhr2g+6wvdJCkmT0fAhtBk2Tn5vUL0MrES5FbRpF4mPZ9LyPwzjkQ3Ji7dz7ABy8/zIe/b638pp8xQTvfki7naEERT/6wlsvfmseOQ8d5ZWQXJo7pzft/SWP+g4O458K2ZO45wphPljLghV9545dN7Dvi+Iuzce9RHvt+Db3//Qv/9+1qSkoVT1+RzKKHL+SW/olsP3T8zHPj0iB1tC4rsmf1mZNt+lmXjKjOE1tCNfJJ9q6Fhe9C99G6D7wraT9UJ5P5N4CPhsGKz107/7nIwne1WSX9+Yp9ePF9de8MV9vxTxyGmQ/rsjepN1Y8rn26zsHaNtf5a2RO1yuaWBeU0mmfrr8nc56Dw+Wy0bPmQ8ve1Svn0iTJa90SjSLxAgEx7Yj8x1yKEwbymM94gmbez5/enMOy7blnDz6yGzZMRXX9Mz+uO8iFr/zGRwuyuK5nS2bfN5Aru8chNlt0bKMg7rmwHfMfHMR7o1JpE9OQV37aSN/nZnPnZ0tZsPkAJ4tL+H7FTka+9wcXvTqXiYt3cGFSEyb/tQ/T/zGAv/SOp2GgH2kJOuoqI7ucTIMf15E45cuKZE6DkBj9ZXaWiEQIbWpdkZzqehgOFz7h/PWsENNR1x1r2Ru+u0PfqEw0l2PycvRNsV165WVx7JUOXG3emv2Mzh6/5OXKb8AJAyCgofPmLXvtuPbprqttl/68XpnNeOj0vmP7dAvsysrGV0ZMRzi8vWLTsxsxRRu9RVAjGoz+CvXLU4ya/xpJeXu47d27uahHMg8ObU94sC0rfPknoEp4YFt3Js9eTqdmYbw3KpVuLSMqnNrP14ehybEMTY5l24F8Pl+UzVdLc5i2eg8Bvj4UlpTSMjKYsekduDqtBZEhZ2egd2oWRpC/D0uyDjEspYwTOzgShjwN39+puw92v0GHHG7+GTpdUb2icSK2fJL5+stV1Zd15eewYyEMf9t6mHF1CI6EUd9oJfLHW9pscNV4XVrFcJoZD2lHdvrzlY+Laq0fNrJtXT1dwa7lOmu85206Ubcy/AJ1IuHGGc4VONw2V69kamrWKkt4Szj/XzoceONMaHfx6Qepqir+VoS9N8m+DdCih2vktIhZkXgTH19kyJNw5Qd089nML6FPsHLpPAa//BtfL82h4ORJjs7/kN9LOzNjVzCPX5bE93/rV6kSKU9i4xAeviSJhWMH88rILlydFsfHN/dkzv0Duf381g6VCIC/rw9dW4SztPyKBHS13JZ9dBhj/kH9BTh5pGaOyPi+cHS39gVVxvFDuptji946Osfd+Pprm/Zlb+gbyocXuiZD+lxh4yxY/wOc/4Au+FkZp/wk813jJykt0SvTkGidPW6F9sP0/9lui+VKADZM1bXjEgZUT86K6HMXNG6vw9cLj2sF6x+iAwOqw6nILc+bt9yqSERkqIhkishmEXnIwfFAEfnSdnyRiCTY9g8RkaUistr2OqjMOQEiMk5ENorIBhHxYHaem+g8Erl5OuGBwo8NnmJkyHLu+2ol9/37JUIL97K++VX8ct/53NQvET/f6v3Jgvx9ubJ7HM/+KYXz20Xj41P1Ej0tPpK1u46Qf7KcScfHRzu4C47AL09o+7FfA0g8v1qyAafDHasKD/35CV0M8JKXPVsyO3U0jP5B2+M/GASbfvLctWsrRSdg+gPQuB30udvaOfH94MhOHeJaU5Z+pLPFL/736cTHqmh7ka68YDU50V47rs1g5ysmVIVfgP4/PpwN817R/pEWPc/Mf3GGRi21IvJCt0S3fRNFxBd4G0gHkoDrRKRcqzpuAXKVUm2AVwH72vgAcJlSKgUYDZQtzfowsE8p1c4272/u+gwepXkqjJmDT5MkHjzyLD+kzOO2oNmcbNCEMbfeSZOwII+LlJYQQUmpYuUOByWwmyRBnzth2f9g1URofYGOSKkujdvpekSV+Ul2LNFZ073v0D3nPU18HxgzRz95f3Y1zH+jfvd9//0VXRdq2EtVF+i0Yzfb1LRcyrF92iyUMABSrrJ+XnCkXk1bVSS7lluvHVcdEgdA52t0CaJ9a6vvHwH9YBXT0Su5JO58pOsJbFZKbVVKFQITgeHlxgwHPrb9PBkYLCKilFqulLIH8a8FgkTE/jhwM/AfAKVUqVKqBqmqtYzQWLhxKnS5jpRN79D1ZAaBPW+s/hNKDekeH4EILMlyYN4COP8hHclSkFdz+7GINm9VpEhKimHqPyG0ma6h5C3CW+h+LUnD4adHtSPe0+XR7az/wXuF+g5shvmvQcpIaOXESjS6g/Yx1TQx8afHtDnoklecd4C3T9eNz3ItrIrstePaDqmenFa46BndXwiq7x+xY6+55eEHHHcqkuZA2di2HNs+h2OUUsVAHhBVbswIYLlS6qSI2NevT4vIMhH5SkSaOLq4iIwRkQwRydi/v4raNrUJ/yC44l1dpyi6Q+XhjG4mLMif9k1Cycg+5HhAYEOdWxLZ2jVPbPH9dNRJ+ZBIgCUf6JDjof9xrmCeOwgIgas/0op05RfWCvG5mqITMPlm7S/yNErBtPvAL0jfBJ3Bx8cWWFGDyK3cbFg5EXr/FaLbOX++/aHHXomhMjKn6xWMO4M6GsboQIUmydWLeixLTJLOyTq2zzWyWcSdisTRY0J5NVnpGBHphDZ33W7b5QfEAfOVUt2BP4CXHF1cKTVOKZWmlEqLjnaQoFSbEYG+d8HfFunKqV4kLSGCZdm5Fee5tLsY/r4MQhrX/GIV1d06slu3Nm09WK8EagMiemWUeJ6tEJ9nv7jsXKpL5mz7DU4e8+y1134LW+fAoEd1KwBnie+nTWLVLUOz7GP9++95e9VjHRHVWju5qwoDzs3S5iZXRmtVRNfr4I75NffDxHinVIo7FUkO0KLM+zigfM2JU2NExA9oBByyvY8DvgVuUErZy84eBI7b9gN8BXigMUf9pUdCJPmFJWzY44HY9JhOuutj+afVWQ/rm+YwF/eorykiMMxiIT5XY/cxlBTqHAdPUXBEl8lp2gV63FK9OezlP6pj3iouhGWf6Jpo4S2qHl8R7dN1PktBXsVjMmecHltX8FK3RHcqkiVAWxFJFJEA4FpgSrkxU9DOdICrgNlKKWUzYU0FxiqlThnNlVIK+AEYaNs1GPBeEf56QGq8DjV2GAbsanx8dD2xso7YLb/qvt8D7tVPkrWN6Ha6kdbKL9xXR8oR2fP102dQuOu6/1lhzn+08/mSV6uXfQ26NlpgWPUammVO1dnxNe3J036YrqlWWfRd5lRtXq6N/3cVEdJYh0N7OHLLbYrE5vO4C5gJrAcmKaXWishTInK5bdh4IEpENgP3AnYv6l1AG+BREVlh22Jsxx4EnhCRVcBfgPvc9RkM0Dy8AU0bBbEkqwI/iatJ6AeHtuiaXsUnbV0PW0G/ezxz/eow4D6dYDb1PseF+FxNcaEuIJh4njYtbpzhmaz73at0eZy0myCuBrZ8H19dMaA6iiRjgv5dO9Oh0BFxafqGW5ESrkntOG8T4/kmV27NbFdKTQOmldv3WJmfC4CrHZz3DODQi6eUygbOc62khooQEdISIlmy7RBKqVPlWNzGKbPHfDi4VddlGvWNDkKorQQEQ/qL8MU18Mfb0N/NSm/3Cig+oX0NqgRWfQk5i52vGOsMpaVaUTaIhMEuMOPF94NNs7RvqWFM1eNBJ4Jum6uvX93VkB0fX62E1/2glX/5yMjNv1S/dpy3adJJB4A4k71fQ0xmu6FK0uIj2HOkgJ2HLfZtrwmxXXQ9pBVfwO8vQdIVNX/69ATth0L7S+C35x1HnbkSuwktvq8OQPDxr3kb2apY/olWVhc945oSMdXpT5LxX10qvdtfan590EriZJ5jGTKn6RVLTaOovEFMR13SJXebxy5pFImhStISPOgn8fWDFr1g80/6pjHUQcn62kr6c/p1hpvzXLIXaNt9SGMICtMmLnf6SfIP6oZVLfvqDp2uoFlXnTth1bxVdEL3Oul4mfUVTFW0GqhDmMv/7ooLdTXrdhfXfOXjDU7V3PJcjlGlikREfEXkU08JY6iddIgNo2Ggn2f9JAADx3o9/Nkp7IX4NvxYs058lVFSDNsXnmnGap+uTYDuqgH28+O6ouyl1Uj+qwhff10OxOqKZO13UHC45k72sgSEaGWSOe3MBL7tC/RKpS6atQBibL3ePegnqVSRKKVKgGhb1JWhnuLrI3RrGU5GRRnurqb7aF0/qddfPXM9V9L7b3q1MO1+HRbsavauhsKjZ2ZA2x3C7jBvbV+kzVp9/nZGd06XEN9f95U5buEBJWOCbtzk6sKJ7dN1rauyN90N0/RKpdUFrr2WpwgIgYgEj3ZLtGLaygLmi8ijInKvfXOzXIZaRlp8JJl7j5J3wgNRSSGN9Y3LVV0PPcmpQnzbdSE+V5PloNR4ozjdcGmDixVJSbFuCRsWB+f9y7Vzg21VpfQKqzL2rNb+mbSbXZ9H1K6cElZKm7pa1bB2nLeJ6VR7TFs2dgE/2saGltkM9YgeCREoBcsdNd8ynElCf+h8rS7E52pzU/YCHQ4d1vTM/e2HwY5FkO/C0nOL39c1qdKf0+VwXE3zVPANrNpPkjFBrxBc5Z8pS2gTaJ522k+ydy3kba+bYb9laZKkzZ3FJz1yuSoViVLqSaXUk8ArwMtl3hvqEV1bhuPrI54zb9V1LnpaO5On3ue6Anqlpdp+7yjMt306oFznmzmyC379ty673uFS18xZHv8gnc9RmSI5eRRWTYJOV7qv3lX7dF1y5uie0wql3VD3XMtTxHTU4csHNnrkclUqEhFJFpHlwBpgra0/SCf3i2aoTQQH+NGpWVjFBRwNZ9IwBgY/qmthrfnaNXPuX68T5ewta8vStIuuxOwqP8nM/9OZ3+kvuLcsTXw/2L1Sl15xxOqvoPCYa53s5bE71TfO0L+/5mnVqyFWm7BHbnkow92KaWsccK9SKl4pFY/OJP/AvWIZaiNp8ZGs2HGYwuIKCjgaziTtZt3+deb/VV7TySqn/CMOViQi+sl6y2woKqjZdTb/ogszDrgfIhNrNldVxPfVbXp3LD77mFKwZIIuqRKX5j4ZYjpCePzpRll13awFuqyLb4DHuiVaUSQhSqlf7W+UUnOAELdJZKi1pCVEUFBUytpdLrgp1gd8fHW/jGP7tJmopmTPh0YtKm5p2z5dJ6Jtm1v9axQV6IizqDbQ7+/Vn8cqLXrqfCFHZeVzMnSUmjuc7GUR0auSXbb2u3U17Lcsvv66WZyHHO5WFMlWW8RWgm17BPBcyqSh1pDmyQKO5wrNu+squYvHaRNOdVFKK5LKyqAkDNBVAWpi3pr/Ohzaaut66OLWso4ICNGrNkcdEzMm6M+TclYVJddjX4VEJLg+zNlbxCTVKtPWzUA08I1tawzc5E6hDLWTmLAgWkYGey4x8Vxh0KO6RtXsZ6s/x4FNkL+/8g56foG6nEzmdO2Yd5ZDW+H3l7Vju7VxDDLBAAAgAElEQVQHcyji+2mTUtm8m+OHYO030HmkZxqZxfeF0Kb6s9emVgU1oUkSHMmBEw5aZbuYKjPbgf9TSv1dKdXdtt2jlDKPpPWUtIQIlmbnoupzr3JnaRCuO11umqXzS6pDtoP8EUe0HwbH9sDu5c7NrxRMe0Db1S92gRnOGeL7acd+Thk/ycqJUFzgXid7WXz94a4lcMHDnrmeJ+j6Z/jHSl2y381YyWyvg1XLDO4iLT6SA8cKyTrohqztc5nUG/WT7rL/Ve/87PnQsEnVvTHaXqR7jDtbe2v9FNj8Mwx65OwcFXfTsheIz+lyKUpps1ZcT+1o9xSBoXUzCbYiGsZoU50HKgBbucJyEZkiIn8RkSvtm9slM9RKetgKOGZUw7xVWqrIP1lsaTvnIsPCW+ib/LL/Od+zRCntQ4jvW7XZJThS9xh3RpGcPArTH9I37R63OiebKwhqpK9t95Nk/Q4HN3luNWKoMVbUbyS6xe2gMvsU2l9iqGe0jm5Iowb+ZGTlcnWa9VanxwuLuerdP1i3u4J8gXKEBPgy85/nERdRh8tUlCftZvh8JGyYCp2usH5ebhYc3VW1WctO+3Tdnjg3u+IIr7LMeQ6O7oaR//PeE3l8f1jyoc7EzpigOz868zsyeJVK/2tsPpJVSqlXPSSPoZbj4yOkxUc4nZj4zNT1rN9zhL8PakPDoMpvVsWlipdnbeSj+Vk8cmlSTcStXbS5UIfvZkxw7iZpN/k4q0g2zoBet1c+du9aWPgupI6GFj2sy+Rq4vvCwrf1Smr9D9DzdvBv4D15DE5R6TdaKVVia4trFInhFKkJEfyyYR+H8guJDKm6MPQv6/fy+aLtjDmvFfde1N7SNdbvPsrEJTv4x4VtCQ3yr/qEuoCPr75hz34GDmyGxm2snZc9X0d9RXewNj6qNTRur1c+lSmS0lL48V4dDDD4cWtzuwt7WPOMh7TjPc0EhtYlrPhIFojIWyIyQES62zcrk4vIUBHJFJHNInJWtx8RCRSRL23HF4lIgm3/EFspltW210EOzp0iImusyGFwLT0SdM0jK/kkB4+d5MGvV9EhNpT7Lmpn+Rq39E/k2MliJmXkVFvOWkm3G3QC3tL/Wj/Hnj/ijNO0fbo+r7LQz5Wfw46FMOQp99WxskpwpM57OLpbN+pq3Na78hicwsp/Zl+gE/AU8LJte6mqk2xmsbeBdCAJuE5EytspbgFylVJt0Kue5237DwCXKaVSgNHAJ+XmvhI4ZkF2gxtIad6IAF+fKh3uSike+mY1R04U89q1XQn0s95trmuLcHokRPDf+dsoLjmHHO+hTXQRxBWf6a5/VZG3U/tIrJq17LQfpp/sN//s+PjxQzDrUe2Y73K9c3O7C/tnNE72OoeV6r8XONjOWiE4oCewWSm1VSlVCEwEhpcbMxz42PbzZGCwiIhSarlSapdt/1ogSEQCAUSkIXAv8IwFGQxuIMjfl5S4RlUmJn6VkcNP6/bywMXt6RDrfCz7Lf1bkZN7glnr9lZX1NpJ2s26+OK676see8o/UklGuyPi0nTP8Yqit35+Qtf/uuRlj4SHWqL7DXprf4m3JTE4iZXqv01EZLyITLe9TxKRWyzM3RzYUeZ9jm2fwzFKqWIgD4gqN2YEsFwpZS+s/zR6VWQSGbxIWkIEq3fmUVBU4vD49oPHefKHtfRpFcUt/atX+G9IUhNaRgYzft45VpEn8TxdyypjQtVjs+fphDJn8yl8fHXP8U0/nR1uvGMJLPsYet8BTWpRIe+mneHyN3VzMEOdwsqjyEfATMDePHsjcI+F8xwFvJdPh650jK1c/fPA7bb3XYE2Sqlvq7y4yBgRyRCRjP3791sQ1+AMafGRFJUoVuWcXcCxuKSUf05agY+P8NLILvj4VK/khK+PcFO/BJZm57LsXGqoJQKpN+lGVHuqcPNlL4CWvbVicJb2w3Tv8bJ90UuK4cd/QmgzGDjW+TkNBgdYUSSNlVKTgFI4tXJw/Bh6JjlA2USDOHS3RYdjRMQPaAQcsr2PA74FblBKbbGN7wOkikgWMA9oJyJzHF1cKTVOKZWmlEqLjo62IK7BGVJtBRwdhQG/P3crS7NzeXp4Ms3DaxbCeXVaC0KD/M69VUnX63V3wMpWJcf26cZEzvpH7LQaqDsLljVvLflAV9R1V9dDQ73EiiLJF5EobCsFEemNNkFVxRKgrYgkikgAcC0wpdyYKWhnOsBVwGyllBKRcGAqMFYpdaosqFLqXaVUM6VUAtAf2KiUGmhBFoOLiQwJoHV0yFkdE1fn5PHqTxu5tHNThndtVsHZ1mkY6Md1PVsyY80ecnLPIWtmcCQkXwmrvtSZ5Y5wNn+kPAEhWplkTtXZ8Ud268KRbS6EjpdXb06DwQFWFMm96Bt+axGZD/wPuLuqk2wrl7vQZrH1wCSl1FoRecqWmwIwHogSkc2269hDhO8C2gCPisgK2xbjzAczuJ8eCZEszc6ltFRbIwuKSrjny+U0bhjIM1ckIy6qojq6bwIAHy/Icsl8tYa0m3X3v9WTHR/Pnq/b9TbrWv1rtE/XhSL3rdMNtkoKYdiL506FW0OtwErU1jLgfHQY8O1AJ6XUKiuTK6WmKaXaKaVaK6Wete17TCk1xfZzgVLqaqVUG6VUT6XUVtv+Z5RSIUqprmW2feXmzlJKJTv3cQ2uJDU+grwTRWzeryOxn5u+gS3783np6i6EB7vOYdo8vAHDUpoycfEOjp0sdtm8XieuBzRJ1uYtR9WUsxfoxk++NUjIbGfrszHzYV2WfcB9ENmq+vMZDA6wFPenlCpWSq1VSq1RSjlZcc5wrmJPTMzIymXuxv18tCCLm/ol0L9tY5df65b+iRw9WcykJTuqHlxXENEZ3HtWwc5lZx47fkiXL3HUn90ZQpvoHuRbf9UKpN8/ajafweCAWhJAbqiLxEcF07hhID+v38v9X62kbUxDHhxqsYyHk3RtEU5afAQT5m+jpPQc6oWSMhL8QyBj/Jn7ty8ElPP5I47oYMvLGPYS+AfVfD6DoRxGkRiqjYgu4Dh7wz5yjxfy6jVdCfKvRpiqRW4dkKgTFNfucds1PE5QmO4CuOZrnaRoJ3u+jupq7oJ2QL3vgFt+1t0TDQY3UKEiKVtXy9HmSSENtZc0W3+Sfw5pR3LzRm691pCkWFpENjj3QoHTbtLdAFdOPL0va57OTnfFCsK/gXcr+xrOeSqr/vuy7TUISANWohMIOwOL0OG3hnrOyB4tCGvgz4jucW6/lq+PcHO/RJ78YR3Lt+fSrWWE26/pEZp20X6MjAnQ6686HHjPKhhwv7clMxgsUeGKxF5XC8gGutuS+1KBbsBmTwloqN2EBfkzMq0FvtXMXneWq9NaEBp4DiYopt2skw+z5+uMd1XqGv+IweABrPhIOiilVtvfKKXWADUIbDcYqk/DQD+u69WS6Wv2sPOwheq5dYVOf9ItZ5eM18rEx0+H/hoMdQArfTXXi8iHwKfo7PZR6ARDg8ErjO6bwPh52/h4QRb/N6xjleOPFBTx7bKdzN24nxJH+RoOGNA2mhv6xOPv66F4lIBg6PpnWPwBRCZCs246M91gqANYUSQ3AXcA9gD0ucC7bpPIYKiC5uENSE+O5YtF2/n74LY0DHT8b7xmZx6fLcrmu+W7OFFUQqvoEEIrGFuWE0UlPP3jOr7K2MGzf0omNd5DTZ9Sb4KF72gTVz8rdVENhtpBld8qpVSBiLwHTFNKZXpAJoOhSm4d0IofV+3mq4wd3NTvdJn6gqISfli5i88WbWfFjsME+ftweZdmjOodT+e4cEtzK6WYtW4vT0xZy4h3/+CatBY8lN6BCAtthWtEdDtIGABZv1e/vpbB4AWqVCS2ulgvAgFAoq2U+1NKKVP1zeA1yiYo3tAngeyD+Xy2aDuTl+aQd6KI1tEhPHZpEiO6x9Eo2LkSIyLCxZ1i6d+mMa//sonx87Yxa90exg7ryFXd46pdFt8SA+6FgsMQ38d91zAYXIyoKmzGIrIUGATMUUp1s+1bpZTq7AH5XEJaWprKyMjwthgGFzNjzW7++ukyOjYNY/3uI/j5CBcnxzKqVzy9W0W6rGjkhj1HeOTbNWRk59IjIYJnrkihfWyoS+YuT/bBfDbtPcbgjjEuk99gqC4islQplVbVOCs+kmKlVJ75pzbUNoYkxdKuSUOOnCji/ovaMbJHC2JCXV8CpENsGJNu78PkpTn8Z/p6Lnnjd24ZkMg/BrclOMDKV6hyiktK+WXDPj5btJ25G3UTtnf/3J30lKY1nttg8ARWViTjgV/QJd5HAH8H/JVSf3W/eK7BrEjOXYpLSvERca+5qQyH8gt5fvoGvszYQfPwBowd1oHz20UTGuR8hd69RwqYuHgHE5dsZ3deAbFhQbr3yto95OYX8vN951cYSGAweAKrKxIriiQYeBi4yLZrJvCMUqqgxlJ6CKNIDK4mI+sQD3+7hsy9uilVq+gQusSFk9K8EZ3jGtGpWSMaBJxdd6y0VLFgy0E+W5TNrHV7KSlVDGjbmFG94xncIQY/Xx+WZucy4t0F3DYgkYcvSfL0RzMYTuESRSIivsBzSqkHXCmcpzGKxOAOikpKWbDlIKtzDrMyJ4/VOXnsOaKfr3wE2jUJpXNcI1LiwunULIxl2bl8tmg72w7kExGsKwJc17MlCY3Pzhd56OtVfLU0h6l/70+H2DBPfzSDAXDtimS2UmqQyyTzAkaRGDzFviMFrMrJY1XOYVbtzGNVTh6H8gtPHU+Nj2BU75akJzettFJybn4hg16eQ+vohky6vY/HTHcGQ1lc6WxfLiJTgK+AfPtOpdQ3NZDPYDgniQkL4sKkIC5MagLonJSdh0+wdtcRWkYG07GptdVFREgAY9M78q+vVzF5WQ4j01q4U2yDoUZYUSSRwEF0CLAdBRhFYjBUgYgQFxFMXESw0+delRrHlxk7+M+09Qzp2MT9CZEGQzWx0rP9JgfbzVYmF5GhIpIpIptF5CEHxwNF5Evb8UUikmDbP0RElorIatvrINv+YBGZKiIbRGStiDzn3Mc1GOoOPj7CM1ckc6SgmBdmmqIShtpLlYpERIJE5G8i8o6ITLBvFs7zBd4G0oEk4DoRKR+CcguQq5RqA7wKPG/bfwC4TCmVAowGPilzzktKqQ7ocvb9RCS9KlkMhrpKx6Zh3NQ3gS8Wb2fZ9tyqTzAYvICV0qafALHAxcBvQBxw1MJ5PYHNSqmtSqlCYCIwvNyY4cDHtp8nA4NFRJRSy5VSu2z71wJBIhKolDqulPoVwDbnMps8BsM5yz1D2hEbFsQj366huKTU2+IYDGdhRZG0UUo9CuQrpT4GLgFSLJzXHNhR5n2ObZ/DMUqpYiAPiCo3ZgSwXCl1suxOEQkHLkMnSxoM5ywNA/147LIk1u0+wicLs70tjsFwFlYUSZHt9bCIJAONgAQL5zmKVywfa1zpGBHphDZ33X7GSSJ+wBfAG0qprQ4vLjJGRDJEJGP//v0WxDUYai/pybGc1y6al2dtZO+ROpMLbKgnWFEk40QkAngUmAKsA16wcF4OUDZmMQ7YVdEYm3JoBByyvY8DvgVuUEptKS8TsEkp9VpFF1dKjbO1B06Ljo62IK7BUHsREZ66vBOFJaU8M9X0lTPULqxEbX2olMpVSv2mlGqllIpRSr1nYe4lQFsRSRSRAOBatCIqyxS0Mx3gKmC2UkrZzFZTgbFKqfllTxCRZ9AKx3T+MdQrEhqHcOfA1vywchfzNh3wtjgGwymsZLY/5mi/UuqpKicXGQa8BvgCE5RSz4rIU0CGUmqKiAShnfnd0CuRa5VSW0XkEWAssKnMdBehe6LsADYAdp/JW0qpDyuTw2S2G84VCopKuPi1ufiKMP2eAQT6VZwdbzDUFFeWSLmvzNsg4FJgvdVcktqAUSSGc4nfNu5n9ITF3DekHXcPbuttcQznMC4rkaKUerncxC9xtonKYDB4iPPbRXNJSlPe+nUzw7s2p2WU81nzBoMrseJsL08w0MrVghgMBus8emkSfj7Cjf9dzJb9x7wtjqGeYyWzfbWIrLJta4FM4HX3i2YwGCoitlEQ/72pJ3knirji7fnMydznbZEM9RgrK5JL0Yl/l6Ed3s2UUm+5VSqDwVAlPRMj+f6ufsRFBHPzR0v4YO5WqvJ5GgzuwIoiOVpmOwGEiUikfXOrdAaDoVLiIoL5+o4+DE2O5dlp67nvq5UUFJV4WyxDPcNKGfll6KTBXHQmejiw3XZMYfwlBoNXCQ7w4+3ru/Pm7M288tNGtu7P5/2/pNIkLMjbohnqCVZWJDPQlXgbK6Wi0Kaub5RSiUopo0QMhlqAiPD3wW15b1QqG/ce5fK35rFyx2Fvi2WoJ1hRJD2UUtPsb5RS04Hz3SeSwWCoLkOTY/nmzr74+/pw9ft/8N3ynd4WyVAPsKJIDojIIyKSICLxIvIwumOiwWCohXSIDWPKXf3p1iKce75cwX+mr6ek1DjhDe7Dio/kOuBxdAFF0D1JrnObRAaDocZEhgTw6a29eOqHdbz/21YysnKJd3HiYlxEMH8f1AY/3+qkoxnOJaxkth8C/gGnuh6GKKWOuFswg8FQM/x9fXj6imQ6Ng3jg9+3urT8vFLwzeGdNGrgzy39E102r6FuUqUiEZHPgb8CJcBSoJGIvKKUetHdwhkMhppzfa+WXN+rpUvnVEpx00dLeGVWJpekNCW2kYkQq89YWZMm2VYgVwDTgJbAX9wqlcFgqNWICE9e3omiUsXTU9d5WxyDl7GiSPxFxB+tSL5XShVxdqdDg8FQz4iPCuFvA9swddVu5m40XUjrM1YUyftAFhACzBWReMD4SAwGA7ef34rExiE89v0ak1Ffj7HSIfENpVRzpdQwpQv5bAcucL9oBoOhthPk78uTl3ci6+Bxxs3d6m1xDF7C6bg9pSl2hzAGg6HucV67aC7prPujZB/M97Y4Bi9gAsANBkONeezSJAJ8fXjs+7WmAnE9xCgSg8FQY5qEBfHPIe34beN+Zq7d421xDB7GkiIRkb4icr2I3GDfLJ43VEQyRWSziDzk4HigiHxpO75IRBJs+4eIyFJbU62lIjKozDmptv2bReQNERFrH9VgMLiT0X3i6dg0jCd/WEf+SWP9rk9Y6ZD4CfAS0B/oYduqbAZvy4J/G0gHkoDrRCSp3LBbgFylVBvgVeB52/4D6IrDKcBo4JMy57wLjAHa2rahVcliMBjcj5+vD89ckczuvAJe/2WTt8UxeBArtbbS0EmJzho+ewKblVJbAURkIjAcKJu9NBx4wvbzZOAtERGl1PIyY9YCQSISCEQCYUqpP2xz/g+d3zLdSdkMBoMbSI2P4NoeLRg/bxsjusfRPjbU2yIZPIAV09YaILYaczcHdpR5n2Pb53CMLRIsD4gqN2YEsFwpddI2PqeKOQ0Ggxd5cGgHwoL8eOS71W5xvO86fIKHvl7FT+v2unxuQ/WwokgaA+tEZKaITLFvFs5z5Lso/19V6RgR6YQ2d93uxJz2c8eISIaIZOzfb7JuDQZPERESwNj0jizJymXy0pyqT7BIUUkp4+Zu4cJXfmPikh3888sV7Dh03GXzG6qPFUXyBNp89G/g5TJbVeSgW/TaiQN2VTRGRPyARsAh2/s4dOn6G5RSW8qMj6tiTgCUUuOUUmlKqbTo6GgL4hoMBldxVWocqfER/Gf6Bg4fL6zxfBlZh7jszXn8e9oG+rSKYuKY3ghw36SVptdKLcBKZvtvjjYLcy8B2opIoogEANcC5VcyU9DOdICrgNlKKSUi4cBUYKxSan4ZWXYDR0Wkty1a6wbgewuyGAwGD+LjIzxzRTJ5J4p4YWZmtefJzS/kwcmruOq9Pzhyooj3/5LKh6PT6N0qiieHd2Jx1iGTUV8LsBK11VtElojIMREpFJESEamy1pbN53EXMBNYD0xSSq0VkadE5HLbsPFAlIhsBu4F7CHCdwFtgEdFZIVti7EduwP4ENgMbME42g2GWknHpmHc1DeBLxZvZ1LGDnYdPmHZZ1JaqpiUsYNBL8/h62U53H5eK36693wu7hSLPeL/T92aMywllld+ymTNzjx3fhRDFUhVf1gRyUCvJr5CR3DdALRVSv2f+8VzDWlpaSojI8PbYhgM9Y5jJ4u59I3fyTqofRmNGwaQ0rwRnePC6RynX6NDA884J3PPUR75bjVLsnJJi4/gmT8l0yE2zOH8ufmFXPzaXBo18OeHu/sT5O/r9s9UnxCRpUqpKtM9rIT/opTaLCK+SqkS4L8isqDGEhoMhnOehoF+zPzneazffZRVOYdZlZPH6pw8ftu4Cbtro2mjIFKaN6JLi3By8wv5aEEWoUF+vDCiM1elxuHjU3HOcURIAC9e3YXRExbzwoxMHrusfKqawRNYUSTHbT6OFSLyArAbXVLeYDAYqiTQz5euLcLp2iL81L78k8Ws232ElTsOs3qnVi6zbOG816S14MH0DkSGBFia//x20YzuE8+E+dsY3DGGfm0au+VzGCrGimkrHtgLBAD/REdWvaOU2ux+8VyDMW0ZDLWfvBNFHDtZTPPwBk6fe6KwhEvf/J38kyXMvOc8GgX7u0HC+odV05aVqK1sdP5GU6XUk0qpe+uSEjEYDHWDRg38q6VEABoE+PLaNd04cOwkj3y/xsWSGarCStTWZcAKYIbtfVeLCYkGg8HgMVLiGnHPhW35YeUuvl+x09vi1CusJiT2BA4DKKVWAAnuE8lgMBiqx1/Pb01qfASPfLeGXYdPeFuceoMVRVKslDJB2gaDodbj5+vDKyO7UFqquP+rlZSarHePYKloo4hcD/iKSFsReRMw4b8Gg6FWEh8VwmOXJbFgy0EmzN/mbXHqBVYUyd1AJ+Ak8AVwBLjHnUIZDAZDTRiZ1oILOzbhhZmZZO456m1xznmsRG0dV0o9rJTqYSuC+LBSqsATwhkMBkN1EBGeG5FCWJAff/t8GXM37jdmLjdSYUJiVZFZSqnLKztuMBgM3qRxw0Beu6Ybf5+4nBsmLCY+Kpg/92rJVaktLCc7GqxRYUKiiOxHN536AlhEuV4gFisA1wpMQqLBUH85WVzCjDV7+GzhdhZnHSLAz4dLUpoyqndLureMOFUE0tUopfhyyQ46x4WT1MxxrbDajtWExMoUiS8wBLgO6Iwu6/6FUmqtKwX1BEaRGAwG0AUhP1uUzTfLdnLsZDEdYkMZ1TueK7o1p2GgpdKDllBK8ezU9Xw4bxuhQX58dmsvOseFV31iLaPGiqTcZIFohfIi8JRS6s2ai+g5jCIxGAxlyT9ZzPcrdvHpwmzW7T5CSIAvV3Rrzh0DWxMXEVyjuZVSvDAzk3fnbOHq1Dj+2HqQowXFfH5bLzo1a+SiT+AZXKJIbArkErQSSUA3opqglKpTaaNGkRgMBkcopVix4zCfLtzOj6t20SDAl3f+3J2+ratf+PHVnzby+i+buL5XS569Ipmc3BNc8/4fnCgqYeKYPrSPDXXhJ3AvrjBtfQwkoxtHTVRK1dkCNkaRGAyGqsg6kM+t/8tg24F8nrgsiVG94532n7w1exMvzdrI1alxPD+i86kS+FkH8rlm3B+UlComjulNm5i6oUxcUbTxL0A74B/AAhE5YtuOWumQaDAYDHWJhMYhfHtnXwa2i+bR79fy8HdrKCwutXz++79t4aVZG7myW3OeK6NE7HN/fltvQLjug0Vs3X/MDZ/Ae1SoSJRSPkqpUNsWVmYLVUrVzRAEg8FgqITQIH/G3ZDGnQNb8/mi7Yz6cBEHj52s8rwJ87bxn+kbuKxLM168ugu+DppxtY5uyBe39aK0VHH9B4vIPpjvjo/gFaxkthsMBkO9wddH+NfQDrx+bVdW5hzm8rfms25XxUaYT/7I4qkf15GeHMsrIx0rETttm4Ty6a29KCgu4foPFpGTe9wNn8DzuFWRiMhQEckUkc0i8pCD44Ei8qXt+CIRSbDtjxKRX0XkmIi8Ve6c60RktYisEpEZImLaoRkMBpczvGtzvvprH0pKFSPeXcD01bvPGjNx8XYe/X4tF3ZswuvXdsPft+pbasemYXx6Sy+OFhRx3QcL2Z1X96sUu02R2PJQ3gbSgSTgOhEp31D5FiBXKdUGeBV43ra/AHgUuL/cnH7A68AFSqnOwCrgLnd9BoPBUL/pHBfOlLv60aFpKHd8toxXf9p4qtTK5KU5jP12NQPbR/P2n7sR4Gf9dprcvBGf3NKLw/lFXP/BIvYeqdtVp9y5IukJbFZKbVVKFQITgeHlxgwHPrb9PBkYLCKilMpXSs1DK5SyiG0LER1OEQbsctsnMBgM9Z6YsCAmjunNValxvP7LJu74bCkTF2/nX5NX0r9NY94blUqgn6/T83ZpEc5HN/dg35ECrv9gIfuPVu2Lqa24U5E0R5dYsZNj2+dwjFKqGMgDoiqaUClVBNwBrEYrkCRgvKOxIjJGRDJEJGP//v3V/QwGg8FAoJ8vL17VmUcu6chP6/by0Der6ZkYybi/pBHk77wSsZMaH8mEG3uw63ABf/5woSXHfm3EnYrEkcepfNKKlTGnB4v4oxVJN6AZ2rQ11tFYpdQ4W7XitOjoaGsSGwwGQwWICLcOaMX/bu7FjX0TGD+6Bw0Cqq9E7PRqFcX40WlkHzzOqPGLOXy80AXSehZ3KpIcoEWZ93GcbYY6Ncbm/2gEHKpkzq4ASqktSmdSTgL6ukpgg8FgqIr+bRvzxOWdCHFhba6+bRrzwQ1pbNl/jL+MX0zeiSKXze0J3KlIlgBtRSRRRAKAa9ElVsoyBRht+/kqYLaqvPjXTiBJROxLjCHAehfKbDAYDF7hvHbRvD8qlQ17jjB6wmKOFtQdZeI2RWLzedwFzETf7CcppdaKyFMiYu9lMh6IEpHNwL3AqRBhEckCXgFuFJEcEUlSSu0CngTmisgq9Arl3+76DAaDweBJLugQw9vXd13kS+YAABjNSURBVGfNzjxu/O8S8k8We1skS1iq/lvXMbW2DAZDXWL66t3c9cVyUuMj+OimHgQHuM6M5gyuqLVlMBgMBi+QntKUV0Z2ISPrELf9L4OCohJvi1QpRpEYDAZDLWR41+a8eFUXFmw5yJhPltZqZWIUicFgMNRSRqTG8fyVnZm7cT9/+2yZU9WIPYlRJAaDwVCLGdmjBc/+KZlfNuzj7i+WUVRS+5SJUSQGg8FQy/lzr3ieuCyJmWv3cs+XKyiuZcrEO6EABoPBYHCKG/slUlyqeGbqejrGhnLXoLbeFukUZkViMBgMdYRbB7SiT6sovl2+09uinIFRJAaDwVCHGJYSy5b9+Wzae9TbopzCKBKDwWCoQ1zcKRYRmLZ6j7dFOYVRJAaDwVCHiAkLIi0+gulrzu7Y6C3qrbO9qKiInJwcCgrqdmeyukJQUBBxcXH4+/t7WxSDoc4zNLkpT/+4jm0H8klsHOJtceqvIsnJySE0NJSEhAR0s0WDu1BKcfDgQXJyckhMTPS2OAZDnWdocixP/7iO6Wt2c+fANt4Wp/6atgoKCoiKijJKxAOICFFRUWb1ZzC4iObhDejSIpwZa2qHn6TeKhLAKBEPYn7XBoNrSU+OZVVOHjsOHfe2KPVbkXiTw4cP884771Tr3GHDhnH48OFKxzz22GP8/PPP1Zq/Jnz33XesW7fO49c1GOob6cmxAMxc6/1ViVEkXqIyRVJSUnmVz2nTphEeHl7pmKeeeooLL7yw2vJVF6NIDAbPEB8VQlLTMKbXAvOWUSRe4qGHHmLLli107dqVBx54gDlz5nDBBRdw/fXXk5KSAsAVV1xBamoqnTp1Yty4cafOTUhI4MCBA2RlZdGxY0duu+02OnXqxEUXXcSJEycAuPHGG5k8efKp8Y8//jjdu3cnJSWFDRs2ALB//36GDBlC9+7duf3224mPj+fAgQNnyFlSUsKNN95IcnIyKSkpvPrqqwBs2bKFoUOHkpqayoABA9iwYQMLFixgypQpPPDAA3Tt2pUtW7a4/fdoMNRn0pNjWZqdy5487/of623UVlme/GEt63YdcemcSc3CePyyThUef+6551izZg0rVqwAYM6cOSxevJg1a9acimyaMGECkZGRnDhxgh49ejBixAiioqLOmGfTpk188cUXfPDBB4wcOZKvv/6aUaNGnXW9xo0bs2zZMt555x1eeuklPvzwQ5588kkGDRrE2LFjmTFjxhnKys6KFSvYuXMna9asAThlUhszZgzvvfcebdu2ZdGiRdx5553Mnj2byy+/nEsvvZSrrrqqer84g8FgmfSUprz800Zmrt3D6L4JXpPDrSsSERkqIpkisllEHnJwPFBEvrQdXyQiCbb9USLyq4gcE5G3yp0TICLjRGSjiGwQkRHu/AyepGfPnmeEx77xxht06dKF3r17s2PHDjZt2nTWOYmJiXTt2hWA1NRUsrKyHM595ZVXnjVm3rx5XHvttQAMHTqUiIiIs85r1aoVW7du5e6772bGjBmEhYVx7NgxFixYwNVXX03Xrl25/fbb2b279iRHGQz1hTYxDWkb09DryYluW5GIiC/wNjAEyAGWiMgUpVRZA/otQK5Sqo2IXAs8D1wDFACPAsm2rSwPA/uUUu1ExAeIrKmsla0cPElIyOnEojlz5vDzzz/zxx9/EBwczMCBAx2GzwYGBp762dfX95Rpq6Jxvr6+FBcXAzq/oyoiIiJYuXIlM2fO5O2332bSpEm89tprhIeHn1pNGQwG75GeHMtbv27mwLGTNG4YWPUJbsCdK5KewGal1FalVCEwERhebsxw4GPbz5OBwSIiSql8pdQ8tEIpz83AfwCUUqVKqQMOxtR6QkNDOXq04qJreXl5REREEBwczIYNG1i4cKHLZejfvz+TJk0CYNasWeTm5p415sCBA5SWljJixAiefvppli1bRlhYGImJiXz11VeAVkgrV6609LkMBoNrSU9pSqmCWWv3ek0GdyqS5sCOMu9zbPscjlFKFQN5QBQVICL2UKWnRWSZiHwlIk1cJ7LniIqKol+/fiQnJ/PAAw+cdXzo0KEUFxfTuXNnHn30UXr37u1yGR5//HFmzZrF/7d39tFdlNkd/1wIGN5FXiQSXKhFBWISkLjSAoIvsAEELVrZaqvruquuLognrlJbdNfTqotVjnJKzip0W6GCdovoAhoVECpqBMHwEhU5IvImiLsCYpTI7R/PE/LLMDO/H/7yAsn9nDPnNzO/79znmefOzJ15ZubOgAEDWLJkCVlZWbRr166GZseOHQwbNoz8/HxuuOEGHnzwQQDmzp3LrFmzyMvLo1+/fixcuBCACRMmMG3aNPr372832w2jHji3Wzt6dmrdsN1bqlonA3A18FTC9N8DTwQ0G4HshOktQKeE6RuAGQnTnQEFxvvpO4GnI8r/ObAaWH3mmWdqkE2bNh0zr6lRUVGhhw8fVlXVVatWaV5eXp2WZ21uGHXDg4vL9awpi/RPX31Tq3aB1ZrC8b4ur0i2Az0SprOBnVEaEckAOgBfxNjcBxwCFvjp54ABYUJV/Z2qDlTVgV26dDn+2jcBtm3bRkFBAXl5eUycOJEnn3yyoatkGMb3YNR53ag8oryyqWG6t+ry8d93gN4i0gvYAUwA/i6geQG4HngTuApY6qNgKKqqIvIiMAxYClwC2Ntv35PevXuzdu3ahq6GYRhpcl73DnQ/tRUvbdjN1QN7JF+glqmzQKKqlSJyO/Ay0ByYraobReQ3uMulF4BZwNMi8hHuSmRC1fIishVoD7QUkSuAEeqe+LrbLzMd2Av8pK7WwTAM42RARPhRTjeefvMTDlQcpl1m/X6uoU5fSFTVxcDiwLypCeMVuHspYcv2jJj/CTC09mppGIZx8jPqvG7M+r+PWfr+HsblB59rqlssRYphGEYjoH+Pjpze/hSWNMAneC2QGIZhNAKaNRNG9uvG8g/3cOjbyvotu15LM46SThp5gOnTp3PoUPrfIVi+fDmrVq1K245hGA1PYU4WFYePsPyDvfVargWSBsICiWEYtc0FvU6jU5uW9Z5a3gJJAxFMIw8wbdo0CgoKyM3N5b777gPgq6++YvTo0eTl5ZGTk8P8+fN5/PHH2blzJ8OHD2f48OGhtvv27Utubi5FRUWASxk/fvx4CgoKKCgo4I033mDr1q0UFxfz2GOPkZ+fz8qVK+uvAQzDqHWaNxNG9DudpeWfUXE4/rtGtYmlkQdYcg/sXl+7NrudB4UPRf4dTCNfUlLC5s2bKS0tRVUZO3YsK1asYO/evZxxxhksWrQIcDm4OnTowKOPPsqyZcvo3LlzDbtffPEFCxYs4P3330dEjqZ9nzRpEpMnT2bw4MFs27aNkSNHUl5ezi233ELbtm2PBhzDME5uCnOyeKb0U1Zu/pzL+tZPBikLJCcIJSUllJSU0L9/fwAOHjzI5s2bGTJkCEVFRdx9992MGTOGIUOGxNpp3749mZmZ3HTTTYwePZoxY8YA8Oqrr9b4cuH+/fstuaJhNEIGndWJDq1asGT9Lgsk9UrMlUN9oapMmTKFm2+++Zj/1qxZw+LFi5kyZQojRoxg6tSpIRYcGRkZlJaW8tprrzFv3jxmzJjB0qVLOXLkCG+++SatWrWqy9UwDKOBadG8GZf2OZ2STbv5tvIILTPq/g6G3SNpIILp1keOHMns2bM5ePAg4LLu7tmzh507d9K6dWuuu+46ioqKePfdd0OXr+LgwYN8+eWXjBo1iunTpx/tOhsxYgQzZlR/I6xqvqV9N4zGx6jzunGgopJVW+rnKxt2RdJAJKaRLywsZNq0aZSXlzNo0CAA2rZty5w5c/joo4+46667aNasGS1atGDmzJmA+9RtYWEhWVlZLFu27KjdAwcOMG7cOCoqKlDVo99Yf/zxx7ntttvIzc2lsrKSoUOHUlxczOWXX85VV13FwoULeeKJJ5J2nRmGceIzuHdn2p6SwZL1uxl2Ttc6L09iciQ2GgYOHKirV6+uMa+8vJw+ffo0UI2aJtbmhlF/TJq3lhUf7uWdey8lo/n363wSkTWqOjCZzq5IDMMwGiFX5HenuQgHKirp2KZlnZZlgcQwDKMRMvzcrgw/t+67tcButhuGYRhp0qQDSVO4P3SiYG1tGI2XJhtIMjMz2bdvnx3g6gFVZd++fWRmZjZ0VQzDqAOa7D2S7Oxstm/fzt699Zsls6mSmZlJdnZ2Q1fDMIw6oMkGkhYtWtCrV6+GroZhGMZJT5Pt2jIMwzBqBwskhmEYRlpYIDEMwzDSokmkSBGRvcAn33PxzkAqmc9S1ZlNs2k2zeaJaDOMH6hql6QqVbUhZgBW16bObJpNs2k2T0Sb6QzWtWUYhmGkhQUSwzAMIy0skCTnd7WsM5tm02yazRPR5vemSdxsNwzDMOoOuyIxDMMw0qM+7uifjAMwG9gDbEii6wEsA8qBjcCkGG0mUAq857W/TmK7ObAW+GMS3VZgPbCOmKc0gFOB/wHe9/UdFKE7x9uqGvYDd0RoJ/t12QA8A2RG6CZ5zcagrbC2Bk4DXgE2+9+OMdqrvd0jwMAY3TS/7mXAAuDUGO0DXrcOKAHOiNsmgCJAgc4xNu8HdiS066gom8AvgQ/8ev02xub8BHtb/W+YLh94q2obAS6IsZkHvOm3qReB9kRs5yF+yonQhfkoymbQT/0idGE+it0fE/yUG2EzzEeRNgN+mhlhM8xHUese9NOYCF2Yj0KPL0Av4G3vo/lAy1o/Xta2wcYyAEOBASQPJFnAAD/eDvgQ6BuhFaCtH2/hnXthjO07gf8mtUDSOYV1+k/gJj/eEn8gTbJMc2A37nny4H/dgY+BVn76WeCGEF0OLoi0xuV3exXoHdfWwG+Be/z4PcDDMdo+uOC3nOqDVJhuBJDhxx9OYrN9wvhEoDhqm/AHhZdx7yp1jrF5P1CUbDsDhvs2OsVPd01lmwT+DZgaYbMEKPTjo4DlMeW/A1zkx2/EHbBDt/MQP82I0IX5KMpm0E9RNsN8FLk/BvzUL8JmmI+i6hn0U05U2SE+irIZ9NOqCF2Yj0KPL7j9coKfXwzcmmy/P97BurYiUNUVwBcp6Hap6rt+/ADuzKF7hFZV9aCfbOGH0JtUIpINjAaeOv7ah9prjztozPJ1+VZV/5zCopcAW1Q16oXODKCViGTgAsXOEE0f4C1VPaSqlcDrwJVVf0a09Thc4MP/XhGlVdVyVf0gMC9MV+LLB3fWlx2j3Z8w2cbNitwmHgN+RYIvj2P7CdPdCjykqt94zZ5kNkVEgL8FnonQKe6sFaAD3k8R2nOAFX78FWB8zHYe9NNlYboIH4XaDPFTxwhdmI/i9sdEP312HPttlM2gnzbE2Qz4KMpm0E9bI3RhPoo6vlyM64mAhH2pNrFAUouISE+gP+5MIErTXETW4boTXlHVKO103EZ/JIWiFSgRkTUi8vMIzV8Ae4H/EJG1IvKUiLRJwfYEXJfVsYWq7gAeAbYBu4AvVbUkRLoBGCoinUSkNdXdBXGcrqq7fDm7gNr+ZuiNwJI4gYj8i4h8ClyLO4sM04wFdqjqeymWe7uIlInIbBHpGKE5GxgiIm+LyOsiUpCC3SG4g+PmiP/vAKb59XkEmBJjawMw1o9fTcBXge080k+p7A8paGv4KaiL81GiNs5PIWVH+iigjfRTxPqE+iigjfRTQBfqo+DxBdgC/DkhMG8nImCmgwWSWkJE2gJ/wPX/74/Sqep3qpqPOxu+QERyQmyNAfao6poUi/9rVR0AFAK3icjQEE0Grgtjpqr2B77CdUXErVNL3Mb6XMT/HXFnpL1w/dNtROS6oE5Vy3FdFK8AL+H6cCuDuvpCRO715c+N06nqvaraw+tuD7HTGriXiCATwkzgLFw/+C5cN0cYGUBHXLfEXcCz/mw2jh8TEfA9twKT/fpMxl+ZRnAjbjtag+tO+bbqj1S381R1cdqgn8J0UT5K1HoboX4KsRnpoxBtqJ9i1v0YH4VoQ/0Uogv1UfD4gusNCFL7j+rWdl9ZYxqAniS5R6LV/ZEvA3cep/37CPTH+vkP4s4ctuLuTxwC5qRo8/4Im91wl8lV00OARUlsjQNKYv6/GpiVMP0PwL+nUMd/BX4R19a4G5hZfjwL+CCZX0jof4/SAdfjblK2TtXXwA+q/kvUAefhzvy2+qESd3XWLQWbPcNs+umXgGEJ01uALjHrlAF8BmTHtOeXVD/uL8D+FNf9bKA0ajsP81OYLsZHodqgn+JshviohjbGT9lJbPaMshnjp6yI9QnzUZjNY/yUwrof9VFg/n24APc51fecBgEvJ9tHj3ewK5I08WeKs4ByVX00ibaLiJzqx1sBl+KeTqmBqk5R1WxV7YnrWlqqqsec6Xs7bUSkXdU47kblhhCbu4FPReQcP+sSYFOS1Ut2lrsNuFBEWvt2uATXhxtWz67+90zgb5LYBXgBdzDB/y5Mok+KiPwIuBsYq6qHkmh7J0yOJdxP61W1q6r29L7ajrsxujvCZlbC5JWE+MnzPK5fGxE5G/dgRFzivUuB91V1e4xmJ3CRH78Y9wRPKAm+agb8E1Acs52H+SnV/SHUZtBPMbpjfBSmjfIT7oQmaPMYH8Wse5ifHo5Y9xo+irEZ5qewdQ/zUdjxpRz31NdVftFa2ZeOobYjU2MZcAe6XcBh3Ib30wjdYNylYtVjiOuAURHaXNzjvGW4g8jUFOoxjJintnD3Pt6j+pG/e2O0+bhHCstwO0HHGG1rYB/QIUn9fo07yG4AnsY/wRKiW4kLXO8BlyRra6AT8BpuR3oNOC1Ge6Uf/wZ31vdyhO4j4NMEPxXH2PyDX6cy3OOV3ZNtEyQ8PRdh82nc45pluANwVoSuJTDHl/8ucHHcNgn8HrglSXsOBtb49n8bOD9GOwn3dNCHwEO4M+PQ7TzET4URujAfRdkM+un5CF2Yj5Luj95PoyNshvkoqp5BP02MKjvER1E2g376aYQuzEehxxfcMaLUt+tzROyj6Qz2ZrthGIaRFta1ZRiGYaSFBRLDMAwjLSyQGIZhGGlhgcQwDMNICwskhmEYRlpYIDGaJCKyXEQG1kM5E0WkXETmBubni8ioWiznOxFZJyIbROQ5/9Y9ItJNROaJyBYR2SQii/07D3VOfbWx0fBYIDGM40RcgspU+QXuXYJrA/Pzce8D1BZfq2q+qubg0mXc4l96W4DL9HuWqvYF/hE4vRbLNQwLJMaJi4j09GfzT4rIRhEp8W/s1jjbFZHOIrLVj98gIs+LyIsi8rGI3C4id4pLVPmWiJyWUMR1IrLKn8Vf4Jdv45P1veOXGZdg9zkReRGX6jtY1zu9nQ0icoefV4x7GewFEZmcoG0J/Aa4xl9FXCMip/l6l/l65nrt/SLytIgsFZHNIvKzFJpuJfCXuDTnh1W1uOoPVV2nqitD2nlDwnSRiNzvxyf6K5kyEZmXpI1a+aufMhGZD7RKoa5GI+B4zqwMoyHoDfxYVX8mIs8C43FvE8eRg8uSmol7m/duVe0vIo/h8oFN97o2qvpX4pJczvbL3YtLSXOjTzdRKiKvev0gIFdVa6RdF5HzgZ8AP8S9Yfy2iLyuqrf4dB/DVfVoihNV/VZEpuJyTt3ubTwBrFXVK0TkYuC/cFct4N5YvhCXKn2tiCxS1bB0/VVXS4W4PFA5uLek0+EeoJeqflOVfoPoNroZOKSquT4Qvptm2cZJgl2RGCc6H6vqOj++BpdILxnLVPWAqu7FJcF70c9fH1j+GTj6TY72/qA4ArhHXCru5bhgdKbXvxIMIp7BwAJV/Urd9yD+F5cU83gYjEvPgaouBTqJSAf/30JV/doHo2W4rK5BWvk6r8blQIvL7ns8lAFzxWV1rsrYHNVGQ/FBXlXL/LJGE8CuSIwTnW8Sxr+jurukkuoTocyYZY4kTB+h5jYfzA+kuCuK8Rr4CJOI/BCXej+MZCneUyHMhgZ+g/MT+Vpd+vBqgyIbqU7WF0diW0LN9hyNCxBjgX8WkX5Et1FU3YxGjl2RGCcrW4Hz/XgqB8swrgEQkcG4j3J9iUsm+Et/oxoR6Z+CnRXAFeKyILfBJShcmWSZA7jvSCTauNaXOQz4XKu/ZTFORDJFpBMuiec7KdQJYClwSuJ9FREpEJGLArrPgK7iPjx2Cu474VWZZXuo6jLcR9ZOBdoS3UaJ65CD65IzmgAWSIyTlUeAW0VkFdD5e9r4k1++GJdlFdy3r1sAZf4G9APJjKj7FOrvcRlW3waeUtW1SRZbBvStutmO+47MQBEpw2VzvT5BWwoswn129oGo+yMh9VJcULvMP/670ZezM6A7jLv5/zbwR6pT5jcH5ojIelxW2cfUfZ45qo1mAm39OvzK19toAlj2X8M4gfFPTx1U1Ucaui6GEYVdkRiGYRhpYVckhmEYRlrYFYlhGIaRFhZIDMMwjLSwQGIYhmGkhQUSwzAMIy0skBiGYRhpYYHEMAzDSIv/B74CZgoocBRfAAAAAElFTkSuQmCC\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:9: 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",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "formula           y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...\n",
       "mean_test_mse                                             0.0229102\n",
       "mean_train_mse                                            0.0227554\n",
       "num_pc                                                           16\n",
       "Name: 15, dtype: object"
      ]
     },
     "execution_count": 34,
     "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.show()\n",
    "best_model = df_res.loc[df_res.mean_test_mse.argmin()]\n",
    "best_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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>degree</th>\n",
       "      <th>formula</th>\n",
       "      <th>mean_test_mse</th>\n",
       "      <th>mean_train_mse</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.023533</td>\n",
       "      <td>0.024149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.053838</td>\n",
       "      <td>0.008531</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...</td>\n",
       "      <td>0.232179</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   degree                                            formula  mean_test_mse  \\\n",
       "0       1  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.023533   \n",
       "1       2  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.053838   \n",
       "2       3  y ~ 1 + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC...       0.232179   \n",
       "\n",
       "   mean_train_mse  \n",
       "0        0.024149  \n",
       "1        0.008531  \n",
       "2        0.000000  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# try doing polynomial fit with different degrees \n",
    "# and find the optimal degree (the model with the lowest testing MSE of 10-fold cross-validation)\n",
    "\n",
    "degrees = [1, 2, 3]\n",
    "kf = KFold(n_splits=10, shuffle=True)\n",
    "mean_mses = []\n",
    "\n",
    "# 10 PCs\n",
    "formula = 'y ~ 1'\n",
    "for num in range(1, 10+1):\n",
    "    formula += (' + PC%d' % num)\n",
    "            \n",
    "for degree in degrees:\n",
    "    # do a d=degree polynomial fit with 10-fold cross-validation and report training and testing MSE\n",
    "    train_mses = []\n",
    "    test_mses = []\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",
    "        # MSE of training set\n",
    "        train_X = PolynomialFeatures(degree=degree).fit_transform(df_train.drop(columns = 'y'))\n",
    "        train_model = sm.OLS(df_train.y, train_X).fit()\n",
    "        train_yhat = train_model.predict()\n",
    "        for k in range(0, train_yhat.shape[0]):\n",
    "            if train_yhat[k] > 0.8:\n",
    "                train_yhat[k] = 1\n",
    "            else: \n",
    "                train_yhat[k] = 0\n",
    "        mse_train = np.mean((train_yhat - df_train.y)**2)\n",
    "\n",
    "        # MSE of test set\n",
    "        test_X = PolynomialFeatures(degree=degree).fit_transform(df_test.drop(columns = 'y'))\n",
    "        test_yhat = train_model.predict(test_X)\n",
    "        for k in range(0, test_yhat.shape[0]):\n",
    "            if test_yhat[k] > 0.8:\n",
    "                test_yhat[k] = 1\n",
    "            else: \n",
    "                test_yhat[k] = 0\n",
    "        mse_test = np.mean((test_yhat - df_test.y)**2)\n",
    "\n",
    "        # store results\n",
    "        train_mses.append(mse_train)\n",
    "        test_mses.append(mse_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",
    "    \n",
    "    mean_mses.append({'degree': degree,\n",
    "                      'formula': formula,\n",
    "                    'mean_train_mse': mean_train_mse,\n",
    "                     'mean_test_mse': mean_test_mse})\n",
    "df_res = pd.DataFrame(mean_mses)\n",
    "df_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "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
}