LogisticRegression_day44_pine.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages for uploading data\n",
    "import re\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# update plot parameters\n",
    "params = {'legend.fontsize': 16,\n",
    "          'figure.figsize': (10, 10),\n",
    "         'axes.labelsize': 20,\n",
    "         'axes.titlesize':24,\n",
    "         'xtick.labelsize':22,\n",
    "         'ytick.labelsize':22}\n",
    "plt.rcParams.update(params)\n",
    "plt.style.use('grayscale')\n",
    "\n",
    "from scipy.stats import ttest_ind as tt\n",
    "from scipy.stats import linregress, norm\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# import k-fold function\n",
    "from sklearn.model_selection import KFold\n",
    "\n",
    "# import RF classifier \n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import roc_curve, auc, confusion_matrix\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "np.random.seed(123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "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>Sample ID</th>\n",
       "      <th>Total Biomass</th>\n",
       "      <th>Fungi Richness</th>\n",
       "      <th>Bact Richness</th>\n",
       "      <th>Bact Diversity</th>\n",
       "      <th>DOC</th>\n",
       "      <th>DOC_cont</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>100_A</td>\n",
       "      <td>232.645161</td>\n",
       "      <td>29</td>\n",
       "      <td>99</td>\n",
       "      <td>2.952329</td>\n",
       "      <td>High</td>\n",
       "      <td>9.375</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>100_B</td>\n",
       "      <td>273.066667</td>\n",
       "      <td>24</td>\n",
       "      <td>70</td>\n",
       "      <td>2.575168</td>\n",
       "      <td>High</td>\n",
       "      <td>13.710</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>100_C</td>\n",
       "      <td>421.933333</td>\n",
       "      <td>12</td>\n",
       "      <td>63</td>\n",
       "      <td>2.539701</td>\n",
       "      <td>High</td>\n",
       "      <td>10.835</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>101_B</td>\n",
       "      <td>398.206897</td>\n",
       "      <td>25</td>\n",
       "      <td>70</td>\n",
       "      <td>3.061506</td>\n",
       "      <td>High</td>\n",
       "      <td>13.960</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>101_C</td>\n",
       "      <td>459.586207</td>\n",
       "      <td>21</td>\n",
       "      <td>95</td>\n",
       "      <td>3.438616</td>\n",
       "      <td>High</td>\n",
       "      <td>9.335</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Sample ID  Total Biomass  Fungi Richness  Bact Richness  Bact Diversity  \\\n",
       "0     100_A     232.645161              29             99        2.952329   \n",
       "1     100_B     273.066667              24             70        2.575168   \n",
       "2     100_C     421.933333              12             63        2.539701   \n",
       "3     101_B     398.206897              25             70        3.061506   \n",
       "4     101_C     459.586207              21             95        3.438616   \n",
       "\n",
       "    DOC  DOC_cont  \n",
       "0  High     9.375  \n",
       "1  High    13.710  \n",
       "2  High    10.835  \n",
       "3  High    13.960  \n",
       "4  High     9.335  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table = pd.read_csv('data/day_44_pine_features.csv')\n",
    "\n",
    "table.dropna(inplace=True)\n",
    "\n",
    "sample_names = table['Sample ID'].values  \n",
    "table.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up X to be the matrix of features \n",
    "feature_names = ['Total Biomass', 'Bact Richness', 'Fungi Richness']\n",
    "X = np.array(table[feature_names].values, np.double)\n",
    "data_df = table[['Sample ID']+feature_names]\n",
    "\n",
    "# set of DOC binary vector \n",
    "DOC = table['DOC'].values\n",
    "DOC_cont = table['DOC_cont'].values\n",
    "DOC_ctg = np.zeros(len(DOC))\n",
    "\n",
    "for i, doc in enumerate(DOC):\n",
    "    if doc == 'High':\n",
    "        DOC_ctg[i] = 1\n",
    "    else:\n",
    "        DOC_ctg[i] = 0\n",
    "        \n",
    "sample2target = {s:t for (s,t) in zip(sample_names, DOC_ctg)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eg. 100A --> 100 \n",
    "def sample_parser(sample_id):\n",
    "    return re.split(\"[ABCDEF]\", sample_id)[0]\n",
    "\n",
    "# function to create a dictionary mapping sample to a set of replicates\n",
    "def get_sample_dict(data_df, sample_parser):\n",
    "    # get the set of samples\n",
    "    sample_ids = data_df.iloc[:, 0].values\n",
    "    sample_ids_norep = [sample_parser(s) for s in sample_ids]\n",
    "    # determine unique set of samples\n",
    "    unique_samples = np.unique(sample_ids_norep)\n",
    "    # create sample_dict to map sample id without rep to set of replicate samples \n",
    "    sample_dict = {}\n",
    "    for s in unique_samples:\n",
    "        # get all replicates of sample s\n",
    "        sample_inds = np.in1d(sample_ids_norep, s)\n",
    "        replicates = data_df.iloc[sample_inds, :]\n",
    "        # store replicates in dict \n",
    "        sample_dict[s] = replicates\n",
    "    # create a dataframe that has the sample names\n",
    "    sample_df = pd.DataFrame()\n",
    "    sample_df[\"Samples\"] = unique_samples\n",
    "    \n",
    "    return sample_dict, sample_df\n",
    "\n",
    "# function to return dataframes \n",
    "def get_Xtrain_Xtest(sample_dict, train_samples, test_samples):\n",
    "    \n",
    "    # create train dataframe\n",
    "    for i, train_sample in enumerate(train_samples):\n",
    "        # pull the dataframe for the sample\n",
    "        if i==0:\n",
    "            Xtrain_df = sample_dict[train_sample]\n",
    "        else:\n",
    "            Xtrain_df = pd.concat((Xtrain_df, sample_dict[train_sample]))\n",
    "    \n",
    "    # create test dataframe \n",
    "    for i, test_sample in enumerate(test_samples):\n",
    "        # pull the dataframe for the sample\n",
    "        if i==0:\n",
    "            Xtest_df = sample_dict[test_sample]\n",
    "        else:\n",
    "            Xtest_df = pd.concat((Xtest_df, sample_dict[test_sample]))\n",
    "    \n",
    "    return Xtrain_df, Xtest_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# keep track of predictions\n",
    "oof = []\n",
    "ytrue = []\n",
    "\n",
    "# create a dictionary mapping sample to a set of replicates\n",
    "# k-folds will divide samples, keeping replicates contained within train and test folds\n",
    "sample_dict, sample_df = get_sample_dict(data_df, sample_parser)\n",
    "kf = KFold(n_splits=10)\n",
    "kf.get_n_splits(sample_df)\n",
    "for j,(train_index, test_index) in enumerate(kf.split(sample_df)):\n",
    "\n",
    "    train_samples = sample_df.iloc[train_index].values.flatten()\n",
    "    test_samples  = sample_df.iloc[test_index].values.flatten()\n",
    "    Xtrain_df, Xtest_df = get_Xtrain_Xtest(sample_dict, train_samples, test_samples) \n",
    "    Xtrain, Xtest = np.array(Xtrain_df.values[:,1:], np.float), np.array(Xtest_df.values[:,1:], np.float)\n",
    "    train_samples = Xtrain_df.iloc[:, 0].values\n",
    "    test_samples = Xtest_df.iloc[:, 0].values\n",
    "    \n",
    "    ytrain = np.array([sample2target[s] for s in train_samples], np.float)\n",
    "    ytest  = np.array([sample2target[s] for s in test_samples], np.float)   \n",
    "    \n",
    "    Xscaler = StandardScaler().fit(Xtrain)\n",
    "    Xtrain = Xscaler.transform(Xtrain)\n",
    "    Xtest  = Xscaler.transform(Xtest)\n",
    "    \n",
    "    LR_model = LogisticRegression()\n",
    "    LR_model.fit(X=Xtrain, y=ytrain)\n",
    "    preds = LR_model.predict_proba(Xtest)\n",
    "\n",
    "    oof.append(preds)\n",
    "    ytrue.append(ytest.ravel())\n",
    "    \n",
    "ytrue = np.concatenate(ytrue)\n",
    "ypred = np.concatenate(oof)\n",
    "ypredproba = np.array([y[1] for y in ypred])\n",
    "ypredmax   = np.round(ypredproba)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAGfCAYAAAByaPfAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3XecU1X6x/HPIzYEAcWCIrZ1QQF1FLAX7KKudQEFewEV15+6Cij2AlhQQZEVy1pZXbGiggoKothQwF5XdNcKisDQYZ7fH+fGCTHJJDOZuZmZ7/v1mlcm955775PkJnlyzrnnmLsjIiIiIoWxStwBiIiIiNQlSq5ERERECkjJlYiIiEgBKbkSERERKSAlVyIiIiIFpORKREREpICUXNUiZtbZzNzMZsYdi/yRmb0WvT7Hxx2LFB8zuzY6P+6OO5aaYmb7R4/5y7hjkXJxvC5m9lB0zEtr6phxUnIVAzO7LzrJUv/mmdl0M7vRzDaJO85CyPA4l5vZbDN71cwuMLO14o5TKmZmjczs/8xsvJl9b2ZLzGyWmb1tZtfUlXM2X2a2ipldGf01iTue6mJm65rZRdHr/52ZLTaz+Wb2mZk9aGZHmNmqccdZzMxsx+g8OTHuWAolOi+uNLPL446lmCi5itcy4Kfo72egMbA9cCHwgZntkVJ+IfAZ8FVNBlkg8yh/rPOB5sCewBDgXTPbIMbYCuUbwuszN+5ACs3MDgS+AG4F9gNaAAuAZkAn4FLgczM7L7Yg47MKcEX0ly25mkU4P36oiaAKycx6A18DNxBe/42BJYTH3ho4HniK8LnVLq44a4EdCedJ3MnVAgr3XbIu4TFVlFx9Hx1zdgGOWfSUXMVriru3iP42JCRXJwK/Eb60HjOzhonC7v62u2/t7vvFFG9V/F/SY12H8Ia8HHBga+D2WKMrAHfvGb0+Y+KOpZDM7AjgWWAj4HOgK9DY3dcF1gD2Al4CGgK3mNnVccVazNx9aHR+XBZ3LPkws6uAfxASxzeBI4Em7t7U3RsBGwInATMI7+Ud4opVcuPub0Tn4kE1eMy+0TH/UVPHjJOSqyLi7gvd/UHg3GhRC8IHWZ3j7nPc/Rog0f/kKDNbO86Y5I/MbFPgfmA1YArQwd1Hu/tCAHcvc/fJ7n4g4QsY4FIzq7EPbak+ZnYY5TUSdwK7u/vT7j4/Ucbdf3b3BwhJ1QWEGnmR+s3d9VfDf8B9hBqbiRnWrwGsiMrclLS8c7RsZpptJkbrTibUIFxJqIJdRGhyfAT4cwVxrQ8MAj4ASglVxx8C1wHrVvKxeiKuDOv/mlRmxyz72QAYnBLbB8C1wDoVxNAOGElo1loIzAHeB4ZmOmZljge8Fj2O45OWXRkte7OCGE+Iyv0ANEizfi/gUeB/hOaYXwi1Rd0z7G//aH9fRvcPBcZF50IZcE6Or9+d0X7mA5tUUHa16DlyYFq254fQLDyU0NS0OHpc/wBaVHCMtQlNkFMJza+LCLVpQ4GWFb0uwDrAjUnvjdlJ5ZoDpwBPAp9Gj7kU+Ai4KV1swENJ52+6v7uTyl6buixavmpS+U2AzYF7gO+i1/o/UcxrZ3leViUkNu9Hj2sWMAbYNXX/ebx3LXrsDryT7rzMtF2W83Av4Pno/F0ITAfOTt4mZV9bARcBLyedK78BbwDnA2tm2O706Ljjk95fr0bHdeCwaHkDQjPnbcB7hG4LSwhNWE8AnXN4vJsCt0TP1fzo72PCD8e907zGmf72qKn3fWq5lH2sQngfvBIdb1l0Pn0UnZcHpnlvZfq7NM175dIsz2UX4HHKz/0fCbWlA8jw/i7Wv9gDqI9/VJBcRWV+isqMTFrWmYqTq3OjDwmPPogWJp3ovwB/ynC8PZI+eDw6sZO3/RZoU4nHWlFy1TWpTKcMZfYiJESZYptJhsSR8AG8Iqlsacq24wt1PNInV39O2m7LLM/T81GZW1OWG+GLPfkDax7hgzJx/yFglZTtfv/wBPpF/5dFj2s5OSRXwJpJj3t4jq/3iUlx7ZTh+TmfkDB4tP/SpG1+zHSeEZLkb5PKLkvZ9hdglyyvy98JX9BOSEDmsXJydWvK8/xb9Fwlx9Y+Zd+3R8sTZX6O7if+hiSVzSW5Ohz4Nfp/bvQYE+veBFZN8/hWB15IeV4S5+9SVv4Bk09ytVfSdkfn+95Pcx6eRngvlkXPbfJzfVOG7acnlSlj5felA28RmqhTt/s9uQLuiP5fET23ZZQnVyUp+1ucck450DfL4+sWnUuJsgsJyVXifiLJaRCdD3Mp/0z5MeVvp6T9Vuv7nuzJ1aMpx/0tOo8S919LKvs0IfFKfo8k/52XVDZjckWoUPhXynHnsPL5nzEpK8a/2AOoj39UXHPVMOlNdEPS8s5UnFzNIXyBHBS9oVchdBz/b7T+32m23YzyD627gDbRdkb4QhsbrfuIHH+9Ju078cY4OcP6u5LKrJdm/ZZJH0j/ICQridjaU/6l8n5qbMCxSft+BNg6Wm6EDrknADcW8Hh/SK6i5VOj5QMyPAfNkz5EUhOSv0fLfwJ6A82SzpHuhJouBy5K2S7x4bmQ8KUyDNggadsKfwUCeyc9fwfk+Ho3TTp3+6asSzw/vxE+eA8hqrEA9qE88ZlBShJBqHFKJFajCRd+NEh6zR6O1n1P6A+U7rjzCYnxgURfSsBWKc/1dYQv3EbRsgZAR+DFaB/T0zzmnGqGyC25mkOomWgXrVuTkCgsidb3SrPf6yhPqv4GNIyWbwE8x8oJST7J1RVJ+22Yz/s+zXlYSviCvjXpPGwGDKc8Adg6zfb3En4w/glYPek5OYJQE+3A0DTbJZKr+YTzfwDQNOkcXT/6fxvCZ8OhhNrqxPm4IaE5dHm0fYc0x9iD8vft+OQy0b6OSvNar1SjluV5q9b3PRmSK2DfaPny6FxaO1qe+Mw8Bbg+ZZutEttU8JiyJVd3Jh338qSYVyVcMNEXOK0y52Bcf7EHUB//qDi5OofyD8NjkpZ3puLkaiFJXxhJ64+h/JfZ6inrEif9Hz6kovWrU/4L8q95Pta0yRXhy/JSyr+Ix2bY/pFo/ZAM69egvCnqyJTl30fLH8gj3kodL1qXKbm6IFr+UYZ9npnhg25dQnPkEtJ8uEdldo+2nU1SQpL04ZnX40/Z91lJ+9ggj+1mRtv8M8PzUwbsmma7tpT/Qj42Zd1gyhOrPzQhET78E4nveRmOuwTYppLPxZqEpkIn9DtKXlfI5GpG6vszKjMiWv9iyvKmlNcu/qGGhfDe/SCX+LK8F9KetznuI/k8HJHhdUs0PV6S5763InwZlwJrpKw7Pem4V1ch/quifdyVZt270bqXSVOjmGF/FSZXNfG+J3NydUm0fEyer0Olkytgu6SYT63sa1Vsf+rQXiQs2NzMLiRc7gzh0v58rzwb7e7pBoZ7hnDyrkF4MySO25DQNAdwc7oduvtSwpcawAF5xpMw1Mx+jP7mEKrnryF8uM4kfJGvxMwaE5JCCH0a0sW2hNBGnxrbAYSr25YTfvVUqIrHy+YRQkLR1sy2S7P+uOh2VMryrsBawKvu/m6GeF4nnCfNyXyV1o05xplq3aT/f81ju8Sl1s0zrH/F3d9IXejuHxP6O0Foykp2UnR7i0efyCnbOqFZATK/Ls+6+ycZo87C3RcTaicgfLFVlyHR+y3VU9Ft+5TlBxNqJBaS5orbaF9pz+UcJF6/fF77bAanLohet2eiu6mPLavoc+5ToBGhJjOd5YTasspKfP6u9JqbWXvCsAoQao+WV+EYqeJ838+Lbjc0s5rKDxLDUnzo7vfW0DGrnQZ8i9feZvaHL4rID4SakXQftNm8k26huy8zs58J1d3rJK3qSPh1C/CWmWXab2JIiFZ5xpPQhPRjAI0Furl7aZp1nSj/VT81z9h2iW7fc/cfc4yxKsfLyN2/N7NJhKavHoQmRQCigTf3jO6mJle7JW7NLNtjSCRBrfjj619KuCihMjI+AVXcbmKWdZMI/VgSX1yY2RaEK2cBnsjynkmcx5lelz8kdKnMrC3Qh9DfaHPCF3fq49m4ov1UQdr3L6GDL6z83oXyL9b3PLqCM43JVY6q6n52928yrMv02ACIrjw9lfD+bEH5+y9ZptfkM3fPmhxGgxifSWhq3CaKI/W7MXX/ic+XWZkSoCqI833/EqGpsxPwspndBbzs7tU5NlviuXy+Go9R45RcxWsZ5b8KnVAV/B/CCX63u8+pxD7nZ1m3OLpdLWnZRkn/b5jD/is7mvop7n4fgJk1J3x5DSFcHdIPSDf2TyI2q0RsifLf5hFjVY5XkVGE5OpYM7s4qfbl2Oh409z90wzxrJXjsdKVmZ2upidHvyT9vy6hs3YuKqrx+C7D8uR16yctSz5HcxlsNtNzNSvbRmbWk9Bkn/hcTHS8TvzAaUxIthrlEENlZXr/pnvvAqwX3Wb78vu+krEkXv91s5bKTb6fSwCY2R2sXKud+MxMDPewbrRdptekote8JSHZ3ypp8QJCP7UyQp+79dLsP/H5kClhrIrY3vfu/pmZ9SH01do7+sPMviZcdTjS3afnu98KVOazuuipWTBeyYOIbuTuW7n7ge5+YyUTq8pInANz3N1y+Otc1QO6+y/u/iShY/FCwrhIh2SJbVaOse2ftG1lal2qcryKjCZ8SW9G+S9TyNwkmBzPjTnG81CafazII8ZUyU1omZpdVmJmTQmXpkO4HD1f6V635M+pRjk8D1ul2QdkeS7MbEPCcB2rEl6LDoTL/NdNvEcJl+tnijEu1RlL4vVvnTyYcU0xs78QEqtEJ+c/EfpWNU96TRK1Rpmeh4rO/2GExOorQgf0ddy9sbtvEO0/MUtG6v6r83mP9X3v7ncRLoY4n9Bk+0t0/yzgPTPrV9l9Z1BM76eCUXIlP0W365hZi6wlCyzqMzEkuntrmnnJErGtb2brk59EdfpmeWxTleNl5e6/EZpAIUqozKw1ofnLCf2yMsXTtpCx5OFNwmXmkPtgtkdQ/mE5KUOZbM1qiV/tyTUOPyX9X13PxaGEGoAPCBckvOfuqYNh5lKbWdMSz9NGWcpkW5fNK9HtqoQa5pqW6At6p7tf4+7/SVMbU+nXxMzWBA6L7h7r7k9F79Nc9l+Zz5dcxf2+x91/dPdb3f0IQi3yLoREy4DroubzQqnO5zI2Sq5kKuGXIcDRMRz/FkI1/J8JA6Ame5vyX2BH5bnfN6PbHfJIGqtyvFwkaqe6RYlkotbqVXf/X5ryiT5C+5hZ2v4o1SnqxP1gdPfEiiZmNrPVCAM+Asxw97cyFN07y24S695LWvYl5Z3kq+scTTy2GemaU6LOvftk2LYsuWihA6vAtOh2xywToO+ZYXlFJlNee9XfzBrkspFl6ayYp8RrMi3dSjPbklCjUlkbUN5PL1NTV6ba6cTny/pm1jGPYybOlWzPUazv+1QevEW42OcHQlNp8ry3uTymbBLPZbrWi1pLyVU952Eai8TVb5dGzSNpmdmq0RV1hTz+HMqnwFnpA9zd51J+ldTl2WqT0sT2IuGDYFXKr76sKJaqHC8XYwh9T9YnjAqdrUkQwmB+Cwk1Ktdn23E1fggPIlxB1Bh4pILHPIzyK74uzlJuPzPbOXWhmbWhPKl9LLE8Snbuj+7+LSqXVnTVbdMsx84kMdn2thnWn0no4P4H7l5G+IEAYeymmjSOULu4FumvuF2N0LyTt+h57x/d7QTcnu0Ksui5P49wQUIhVPSaDKri/ucl/f+HKxWjHxN90m3o7h9S/gPghjS17hUdM9t5Etv73sxWz7TOwxWRiR/iayStSjymVSo5hdkD0W07MzutEtsXJSVXAuED9FdC88EUMzvKzH5/85jZVtGH5ieEqwsL7RbCm/ZPlCccCX0JnUtbAm+Y2REpsf3ZzC4gXJJdklju4SrLRC3KCWb2r8SXcvQlsJGZ9Taz1Mu0K3W8XLj7IsqTt4GEwVqXUT7MRWr5WYSxwADOiB5Du6RY1jSzPc1sBGFqj4Jz95mEgQOXEy5Hn2pmf030wTGzVcxsDzN7gZCAAAx097FpdxjMA540s4MTtRxmtjeh2TQxLtPjKdsMJAzZ0Rh41cxOSE70zGxTM+tFqOX4SyUe6kuE5tntzezWRIJmZk3NrD8hcfwly/YfRbcn5lrDUwjRD4Kh0d1BZtYnau7CzDYnTOFS2St8cfdnCM89hNf3teg9kfzcb2BmJxKe+1tI0zG9kl6Kbs82s5OjRBEz28zMHiI0G1a6b2rUBJi4yu4+M9s+2v8qZnYA5WMHZpKY/WEfYKyZdUisMLP1zew4M3sgZZvEebJtphqvmN/3N5jZv6PX+PcLGcyshZkNJ5xLZZQPS4K7z6b8YpdT8j2gu79P+Q/sf5jZZYkfttGP2NZmdrWZnVHJxxQPL4LBtghfMv9HGGTsU8oHlsxrwMo0++1BqNqeS7g0dSrhl8gqhYi7CnHdR5ZBRLNs15mKBxE9Ocv2M6MyndOs60S4UisxmNsyQlPM4qRlTjRXVh4xe0VxReUSI2x/nPr6ADtTPiJxtth2T7Pfvqw8ZcR8Kp7+plLHI8MgoillDk7ZxzM5PIdXpDyGBYRkOHlany9Stsk4vUUlz9kurDzNSxnlV20lli0C/p5lH+mmv1lA7tPf/JnygTw9evyzU15PB3rm+7pE5Yal7Cf5OX6OUFPipAwCGm17RtJ2CwlXkc0EBieVyWluwQyxZRyokVCLMD5pH0spH5V9CaEpNbFu/Uq+/mcTkuLk5+e36PVLXjaDpIFaczkPyTCwZvS43k7a93JWHm3+4kyvbaZ9pjn2bqw8fU1p0v3ZhL6GaZ/3aPuelI+en3jt/zD9Tco2ryetnx2dJzOBjjX1vs9UjjBWWvLrOTfN655usNrrktYnZkKYSdI0W2QfoX1Nwo/MbJ8vtWr6m2KpuTqLMNBbT0KiVeU2+yjLfphQ0zKZ8CuoNeHkGV2Tvy5rA3d/B9iaMCzCFMIbpBnhg2YqoXq6k7tn6qRcVYmmu20oH8gzEdtbhPPiYkJ/hNKU2AYRRjJ+PXWn7n4DodP4fYQvvNWi7d4nnHN/T7NNpY+Xg5dYeUiDTE2CyfFcRRjP6G5C/yMjXBr+A6G250xWvgKx4DzURG1FSIxeJjyGxoQP3ncJtRut3X1Ixp2Um0VI5odF/69GSOxHAiXu/lmGGL4g1BaeQ/gxMYcwQvlywpf6MMIQH/9Kt30Oj/FcwnM5nfCF2YDQ9PM3Qm1YxiuwPFxh1ZtQE7KCcMXkZpQPlVBtPAxs24VQU/sR4YtpOaED8p6U124kviwrc4w7CP2b+hM6uv9IGG/KCRNnP0joHL6DV3Kg1jTHXEKYjuUGwtRIicf1AnCIu1e1WRB3n0J47zxNOJ9WJTy2fxDOtaxjRbn7w4SO53cQpuOxKM6PCVN7pavJOSLa/9eEicg3i/7WTNl3HO/7mwgVHc8QXlcjJLnfEi662SP6TE11BeHz8gPC+ybxmHJqJnf3xe7+V0K3gGcJnwuNCcnnG4SR42vVAKOJeZTiDcLsdELiM5XwQX0PoWNrV3dP22RSwf6OIWTBPwJ7RR/KicutXyF8gZ/n7kMz70VECsnMXiM0LZ7g6S8fl2pgYRDOccBXnnmYChEpoKIYRNTd706+X4CLTRKdafslEqvoOD+Z2VmEX739zew2D51RRUTqnKhP24XR3ZeylRWRwimWZsGCia7w6EDod/BY6vqoWes7wjQKu6SuFxGpTcxs9agT8kHJV0pamP/ucUL/mqWUD4IqItWsKGquCiwx19ZHHq7OSucdwtVgOxD6F4mI1FarEK6c6wpgZvMIn+2Jca9WAGd7mBhbRGpAnau5onxQuWxzPiXmMKrKAHQiIsVgKeEq6GcIV2E2iP5mEsYQ6uTu98QWnUg9VBdrrhLjryzIUqY0uq3MgGciUgnuvkfFpSRfUb/RO6I/ESkCdTG5SvSGr/RlkNFghL0AGjVq1GHrrbcuRFwiIiL12meffcaiRYto2LDG5wJn8eLFLF++nFVXXZXly5fPdveCziGbrC4mV/Oj22zTdCTWzU+30t1HEsbcoWPHjj516tTCRSciIrEZOXIko0ZVOLycVJMGDRqwxx57MHHixBo53oIFC3B3GjduzDvvvMMvv/zCwQcfjJll6zpUZXWxz9XM6HazLGUS00HMzFJGRETqmFGjRjF9eqZ5mqW6lZSU0KNHjxo51vjx42nfvj39+vUDoFOnThx88ME1cuy6WHOVmEG9nZk1zHDFYKeUsiIiUk+UlJTUWM2J1Lw5c+Zw4YUXcu+999K6dWuOPfbYGo+hziVX7v5fM3uPMOVJV8pn3AZ+nyB2E8Lo7W/UfIQiIsWrrjebTZ8+nZKSvOZcl1rk1VdfpXv37syaNYv+/ftzxRVXsOaaa1a8YYHV2mZBMxtkZp+aWbr5pRLLrjezrZK22YDyK2oGa3R2EZGV1fVms5pslpKa17JlS7bcckvefvttBg0aFEtiBUVSc2VmO7LyZcRto9uBZpaYugF3Tx5RfSPC5Lobpe7P3Ueb2QjChNAfmNl4wuza+wFNgKcIEziLiNQ61Vm7lKjZUbOZ1AbuzoMPPsiECRO47777+NOf/sTrr78ed1jFkVwREp6d0yz/c2V36O5nRxPF9iFMAt0A+JQws/YI1VqJSG2VqF2qjuYt1exIbfHNN9/Qu3dvXnjhBXbbbTfmz59PkyZN4g4LKJLkyt0nUj4+Va7bnAycXEGZUUDd7TwgIvWWapekviorK+OOO+6gf//+ANx2222cffbZrLJK8fR0KorkSkREcm/uU6dsqc/mzZvHtddey+67786dd97J5ptvHndIf1A8aZ6ISD2Xa2dyNd1JfbNs2TLuuusuli9fTrNmzXj77bcZN25cUSZWoJorEZGiouY+kZVNmzaN0047jWnTptG8eXOOPvpoNt1007jDyko1VyIiIlJ0Fi9ezIABA+jUqRPff/89jz/+OEcffXTcYeVENVciIiJSdHr27MkTTzzBySefzM0338w666wTd0g5U3IlIkJxjEyujupS382fPx8zo3Hjxlx88cX07t2bAw88MO6w8qZmQRERimNkcnVUl/rshRdeoH379vTt2xeAjh071srEClRzJSL1VGpNlUYmF4nHr7/+yvnnn88DDzzA1ltvzfHHHx93SFWmmisRqZdSa6pUayRS8yZNmsQ222zDqFGjGDBgANOmTWO33XaLO6wqU82ViNRbqqkSiVerVq1o06YNw4YNq1P9DZVciUitUOgO5+o8LlLz3J3777+fCRMm8MADD7Dlllvy6quvxh1WwalZUERqhUJ3OFczoEjNmjlzJgcddBCnnHIK33zzDfPnz487pGqjmiuRWqYYhgyIgzqci9ROK1asYPjw4VxyySWYGcOHD+fMM88sqomWC63uPjKROqoYhgyIg2qaRGqn+fPnM2jQIPbcc08++ugjzj777DqdWIFqrkRqJdXgiEgxW7ZsGffccw+nnXba7xMtb7LJJphZ3KHVCCVXIkUu03hMIiLF6N133+XUU0/l/fffZ4MNNuDoo4+mVatWcYdVo+p2vZxIHaDxmESkNli0aBH9+/dn5513ZtasWTz55JO1ZqLlQlPNlUhMcu2Yro7cIlIb9OjRg6eeeorTTz+dG2+8kWbNmsUdUmxUcyUSk1w7pqumSkSK1bx5834fUmHAgAGMHz+eu+66q14nVqCaK5FYqUZKRGqrsWPH0rt3bw499FBGjBhBx44d4w6paKjmSkRERHI2e/ZsTjjhBA455BDWXnttTjrppLhDKjqquRIREZGcTJw4kW7dujFnzhwuu+wyBgwYwBprrBF3WEVHyZWIiIjkZLPNNqNt27YMGzaM7bbbLu5wipaSK5FKquo0NBqvSkSKnbtzzz338PLLL/Pwww+zxRZbqJ9oDtTnSqSSqjoNja4CFJFi9p///If999+fM844g++//75OT7RcaKq5EslDcm2Vxp8SkbpoxYoVDBs2jAEDBrDqqqty5513cvrpp9f5+QALSc+USB6Sa6tU8yQiddH8+fO58cYb2Xffffn444/p1auXEqs8qeZKJE+qrRKRumbp0qXcfffdnHHGGTRr1ox33nmHjTfeuN5MtFxoSkVFRETqsXfeeYeOHTvSp08fnn32WQBatmypxKoKlFyJiIjUQwsXLuSiiy5il1124ddff+WZZ57hqKOOijusOkHNgiIiIvVQz549eeqpp+jVqxc33HADTZs2jTukOkPJlYiISD0xd+5czIwmTZpw6aWXcu6557LPPvvEHVado2ZBERGRemDMmDG0bduWvn37AtChQwclVtVEyZWIiEgdNmvWLHr06MHhhx/Ouuuuy6mnnhp3SHWemgVFRETqqJdffpnu3bszd+5crrrqKvr378/qq68ed1h1npIrERGROmrLLbdk++23Z+jQobRr1y7ucOoNNQuKiIjUEWVlZYwcOZLjjjsOd2fzzTdn/PjxSqxqmJIrERGROuDLL79kv/32o3fv3vz000+UlpbGHVK9peRKRESkFlu+fDlDhgxhu+2247333uOuu+5iwoQJrL322nGHVm+pz5WIiEgtVlpaypAhQzjggAO44447aNmyZdwh1XtKrkSyGDlyJKNGjfr9/vTp0ykpKYkxIhERWLJkCSNHjuTMM8+kWbNmTJ06lY022kjzARYJNQuKZDFq1CimT5/++/2SkhJ69OgRY0QiUt+9+eab7Ljjjpx77rk899xzAGy88cZKrIqIaq5EKlBSUsLEiRPjDkNE6rkFCxZw2WWXceutt9KyZUuee+45DjnkkLjDkjSUXImkSG4KVDOgiBSLHj168Mwzz3D22WczaNAgmjRpEndIkoGaBUVSJDcFqhlQROL022+/MW/ePACuuOIKJk2axPDhw5VYFTnVXImkoaZAEYnb008/zVlnncVf/vIX7rzzTnbccce4Q5IcqeZKRESkiPz8888ce+yxHHnkkay//vqcccYZcYckeVLNlYiISJGYMGEC3bp1o7S0lGuuuYZ+/fqx2mqrxR2W5Kmoaq7MrIeZTTazuWZWamZTzayPmeUdp5mtY2YDzewDM1tgZkvM7Bsze9DM1ENZVjJy5EgxZBNOAAAgAElEQVQ6d+5M586dVxp6QUSkJm211VZ06NCBadOmcemllyqxqqWKJrkys+HAw0BHYDLwEtAauB0YbWYN8tjXpsB04GKgBfAKMAZYBhwPvGNmxxT0AUitpk7sIhKHsrIyRowYQbdu3XB3NttsM1588UXatm0bd2hSBUXRLBglOmcDPwJ7ufsX0fINCYnRUcA5wNAcdzkY2BR4Hujq7guj/a0CXA5cAdxpZs+4+7JCPhapHTKNvK5O7CJSUz7//HNOP/10Jk+ezP77709paanmA6wjiqXm6uLotl8isQJw95+As6K7/fNoHtwnur0mkVhF+ysDrgEWAc2BP1cpaqm1NPK6iMRl+fLl3HDDDWy//fZ88MEH3Hvvvbz44otKrOqQ2GuuzGwToAOwFHgsdb27TzKz74CWwC7AlBx2uyTHw8/ONU6pe1RTJSJxWLBgAUOHDqVLly4MHz6cjTbaKO6QpMCKoeZqh+j2I3dflKHMOyllKzIuur3UzNZKLLQw8dLlQEPgGXf/Od9gpfZSp3URicuSJUsYOnQoS5cupWnTprz77rs8/vjjSqzqqGJIrraIbr/JUubblLIVuRR4GzgU+MbMxpjZaOBzQhPkQ4SO7VKPqNO6iMThjTfeoKSkhPPOO4/nn38egBYtWmii5Tos9mZBoHF0uyBLmdLoNqcGaXefbWb7AsOBk4DDklZ/Bkxy9/n5BirFI7VDei7UaV1EalJpaSmXXnopw4YNo1WrVowdO5aDDz447rCkBhRDzVUidfeC7dBsa2AacBBwArAR0AzYj5DE3WVm92bZvlc0xtbUWbNmFSosKaDUDum5UG2ViNSknj17MnToUPr06cOHH36oxKoeKYaaq0QNUuMsZRLrKqxtMrNVgceBrYDd3f2NpNUvm9kBwMfAKWb2oLu/kroPdx8JjATo2LFjwZI+KSzVQolIsZkzZw6rrLIKTZs25corr+Siiy5ijz32iDssqWHFkFzNjG43y1KmVUrZbHYG2gL/SUmsAHD3X81sLHAysD9hHC0pQtma/hJNfCIixeLJJ5/k7LPP5vDDD+fOO+9khx1yvQZL6ppiaBacFt22M7OGGcp0SimbzabR7dwsZX6LbtfNYX8Sk2xNf2riE5Fi8eOPP9K1a1eOPvpoWrRoQe/eveMOSWIWe82Vu//XzN4DdgS6Ag8krzezvYFNCKO3/6EmKo3vo9utzayZu/+Wpswu0e3XlYtaaoqa/kSkmI0fP55u3bqxcOFCBg4cyIUXXqj5AKUoaq4ABkW315vZVomFZrYBcEd0d3A0wnpi3Tlm9qmZrZSMERKw7wljWd1jZk2StlnFzC4lJFfLCX2zREREKqV169bsvPPOTJ8+nYsvvliJlQBFUHMF4O6jzWwEYaqbD8xsPGGS5f2AJsBThAmck60HtCHUaCXva6mZnQw8DRwN7G1m7xCmvCkhjJVVBpzn7l9V24MSEZE6JzHR8iuvvMJjjz3GpptuytixY+MOS4pMsdRc4e5nAz2B94C9CcMofEmYsPkYd1+Rx75eArYH/gH8AnQmjHW1KvAI4SrC4YWMX0RE6rbPPvuMvfbai3POOYfS0lJKS0sr3kjqpaKouUpw91FATiNDuvuVwJVZ1n9B+aTPIiIilbJs2TJuuukmrrrqKtZaay3uu+8+TjzxRI2wLhnlnVxFfaKOBbYBGrn7kdHyTYDtgNfcfV5BoxQREYnJwoULuf322/nLX/7CbbfdRosWLeIOSYpcXsmVmfUFrk3aLnmAzYbAGEIz3oiCRCciIhKDxYsXM2LECPr06UPTpk1577332HDDDeMOS2qJnPtcmdlRwGBgCrAHMCR5fdQMNw04opABioiI1KTXXnuNkpISLrjgAsaNGwegxErykk/N1fmEEdIPdvfF0TQyqT4C9ipEYFJ/ZBqJXaOwi0hNmj9/PhdffDHDhw9n880354UXXuDAAw+MOyyphfK5WrAEGOvui7OU+R5Qei95yTQSu0ZhF5Ga1KNHD+644w7+7//+jw8++ECJlVRaPjVXDYClFZRZL4cyIn+gkdhFJA6//vorDRo0oGnTplxzzTVccskl7LrrrnGHJbVcPjVXX1E+bcwfWLgmdTfgk6oGJXXbyJEj6dy58+9/meYPFBGpTqNHj2abbbbhoosuAsKPPCVWUgj5JFejgZ3M7MwM688DtgYerXJUUqelNgOq+U9EatIPP/zAMcccQ9euXWnVqhV9+vSJOySpY/JpFhwCdAeGm1lXYDUAM7sS2JMwCvp0yucCFPldcqf1REd1NQOKSE176aWX6NatG4sXL+b666/nggsuYNVVi2o8bakDcq65cvcFhGlpniQkUnsABlwO7EOY/+8Ad1efK/mD5Noq1VSJSE1zD8MytmnTht12240ZM2bQt29fJVZSLfI6q9x9NvBXM2tJ6H/VHJgLvOnu31RDfFKHqLZKRGraihUrGD58OJMmTWL06NFsuummPPfcc3GHJXVcpVJ2d/8OeLzAsUgtl2m8KtCYVSJS8z755BNOO+003njjDbp06cKCBQto3Lhx3GFJPZDPCO3zzKxfBWUuMrO5VQ9LaqNM41WBmgJFpOYsW7aMa6+9lpKSEj7//HMefPBBnnvuOSVWUmPyqblqDKxRQZnVo3JST6ijuogUm0WLFnHnnXdy1FFHMWzYMDbYYIO4Q5J6Jp+hGHLRFFhS4H1KEVNHdREpBosWLeKmm25i6dKlNGnShPfee49HHnlEiZXEImvNlZntmLJo4zTLIIzevilwHPBFgWKTWkK1VSISp1dffZXTTz+dL774gtatW3P44Yez/vrrxx2W1GMVNQtOBTz634Ezor9MDLiiAHGJiIhkNW/ePPr378+IESPYYostGD9+PPvtt1/cYYlUmFzdTEiqDLgAeAOYkqbcCuAX4GV3f7egEYqIiKTRs2dPnnvuOc4//3yuueYaGjVqFHdIIkAFyZW7X5j438xOAp5095uqPSoREZE0fvnlFxo0aECzZs247rrrGDBgALvsknHaW5FY5DNC+/pKrEREJA7uzqOPPrrSRMvbbbedEispSoW+WlBERKSgvv/+e4488kiOPfZYNttsM/72t7/FHZJIVnmP0G5mXYCDgJakH/fK3f2IqgYm8co22noyjbwuItXphRdeoHv37ixZsoQbb7yR8847T/MBStHL+Qw1s1UJkzYfQujgnujonuBJy6WWS4xfVVHipLGtRKQ6uDtmRtu2bdlzzz255ZZb2GqrreIOSyQn+aT/FwKHArcAQ4GZwEDgPqAzcDUwEehdwPgkRhq/SkRq2ooVKxg2bBiTJk3iySefpFWrVowZMybusETykk+fq+OA99397+7+bbRsqbt/6e53A/sARwAnFTpIERGp+z788EN23313LrjgAlasWMGCBQviDkmkUvJJrrYCJifdd2C13++4fwY8S/ZBRkVERFaydOlSrrrqKnbccUe++uorRo0axTPPPKOJlqXWyie5WgGUJt0vBZqnlPmakISJiIjkZPHixdx999107dqVjz/+mOOOOw4zq3hDkSKVT3L1HbBJ0v0vgdQBRtoDv1U1KBERqdsWLlzIDTfc8PtEy9OmTePhhx/WnIBSJ+STXE0Bdk66/wywvZkNNbO9zewK4GBWbjqUWmTkyJF07tyZzp07M3369LjDEZE6auLEiWy33Xb069ePF154AYD11lsv5qhECief5OoR4Ccz2zy6fzPwEfA34GXChM3fAf0LGJ/UoMTwC6AhFkSk8ObOnUvv3r3ZZ599AHj55Zf5y1/+EnNUIoWX81AM7v4S8FLS/flm1gk4ltDPaiYw2t3nFjpIqTkafkFEqkvPnj0ZO3YsF154IVdddRVrrbVW3CGJVIsqDXPr7kuA+wsUi4iI1DGzZs1itdVWo1mzZgwcOJDLL7+cnXbaKe6wRKpVQecWtEDjXImI1HPuzqhRo/4w0bISK6kPCpZcmdkxwIfAvYXap4iI1D7/+9//OPzww+nZsydbbbUV5513XtwhidSoCpsFzawxYWDQTsAywtWA97n78mh9Z+AmYAfC3IIvVlewIiJS3MaNG0e3bt1Yvnw5N998M+eeey4NGjSIOyyRGpU1uTKzZsAbQGvKJ2k+HjgKONTMbgHOjdZNBC5z99erLVoRESlKiYmW27Vrx7777svNN9/MlltuGXdYIrGoqOaqH9AG+JwwFIMR5hg82MyeAI4EPgDOd/eXqzNQEREpPsuXL+fWW29l8uTJPPXUU7Rq1Yqnnnoq7rBEYlVRn6vDCGNX7eDuV7r7FUAH4EfCJM1PAR2UWImI1D8ffPABu+222+8d1jXRskhQUc3VFsAD7r4osSAa3+oZoBfQP9H3SorXyJEjGTVqVIXlpk+fTklJSQ1EJCK12ZIlSxg4cCADBw5knXXW4dFHH6Vr166aD1AkUlHN1VqEWqpUiWVfFjYcqQ7JI69no1HZRSQXS5Ys4d577+XYY4/l448/plu3bkqsRJJUdRDRskIFItVLI6+LSFUsWLCAYcOGccEFF9CkSROmT59O8+bN4w5LpCjlkly1NbOjU5cBmNlRlF9F+Dt3f6IAsUkauTbxJVNzn4hUxYQJEzjjjDP4+uuv2XbbbTnssMOUWIlkkUty1TX6S2XA6AzbaFCTapJo4ssnWVJzn4hUxm+//cZFF13E3XffzZ///GcmTZrEXnvtFXdYIkWvouTqCcBrIhDJnZr4RKQmHH/88YwbN45+/fpxxRVX0LBhw7hDEqkVsiZX7v7XmgpERETi9/PPP7PaaquxzjrrMGjQIK666io6dOgQd1gitUpBJ24WEZHayd156KGHVppoedttt1ViJVIJSq5EROq5b7/9lkMPPZQTTjiBNm3a8Pe//z3ukERqtSoNxSAiIrXb2LFj6datG2VlZQwdOpQ+ffpoomWRKlJyJSJSDyUmWt5222054IADGDJkCFtssUXcYYnUCUWVXJlZD+AsYDvCcA6fAv8ERlRmwFIzawCcDvQE2gGNgFnAdGCku48pUOgFlW0sK41ZJSJVsXz5coYMGcLkyZMZM2YMm2yyCU88oaEJRQqpaPpcmdlw4GGgIzAZeAloDdwOjI4SpXz2ty4wBfgHIVl7kzDR9LfAfoSJp4tStulqNGaViFTWjBkz2Hnnnenfvz9rrLGGJloWqSZFUXNlZscAZxPmLNzL3b+Ilm8IvAIcBZwDDM1xf6sAY4CdgLuB89x9QdL6xsDmBXwIOcl3AmWNZSUihbB48WKuvfZarr/+epo3b87o0aM55phj4g5LpM4qlpqri6PbfonECsDdfyI0EwL0j5KmXJwB7AZMAnolJ1bRfkvd/cMqxpw3TaAsInFYunQpDzzwAD179uTjjz9WYiVSzfKuuTKzfQh9mLYBGrl7SbS8NbA/8HiUFOW6v02ADsBS4LHU9e4+ycy+A1oCuxCa+ipyTnR7vbsX1QjzqpESkZpQWlrK0KFDufDCC2nSpAkzZsxgnXXWiTsskXohr+TKzO4AehPmFVzOynMILgSGAWsBN+Wx2x2i24/cfVGGMu8QkqsdqCC5MrMWQHtgGfCKmW0LHANsBPwCvOLuL+URX5UkNwWqM7qI1ISXXnqJXr16MXPmTEpKSjj00EOVWInUoJybBc3sdOBM4FFgE2Bg8np3/x+h0/ihecaQuPb3myxlvk0pm8120e1M4DJgBnAF0IvQ/PiimU0ys/XyjLNSkpsC1dwnItVpzpw5nHrqqRx44IGsscYaTJ48mUMPzfcjWUSqKp+aq97AR8Dx7l5mZuma2z4nNA3mo3F0m+2yldLodu0c9rdudLsFcAnwADAI+I5wJeJwYC/g38C+6XZgZr0IyRibbrppDodcWbraKjUFikh1O+GEExg3bhwXX3wxl19+OWuuuWbcIYnUS/l0aG8LjK9gvKkfgQ3yjMGi20L1jUo8plWBCe5+krt/6u7z3f0V4EBgEbCPme2dbgfuPtLdO7p7x/XXXz/vAFRbJSI15ccff2TOnDkADB48mLfffpuBAwcqsRKJUT41VyuA1SoosxHZa6DSmR/dNs5SJrFufpYyqfsDGJm60t3/Z2bPAX8ljHc1KZcg86XaKhGpTu7O/fffzwUXXMDRRx/N3XffTfv27eMOS0TIr+bqU0JzWlpmtjrQmdDHKR8zo9vNspRplVI2l/0BfJ2hTGJ5ixz2JyJSVGbOnMnBBx/MKaecQrt27bjoooviDklEkuSTXD0MtDez6zKsHwRsSujjlI9p0W07M2uYoUynlLLZfEp57VnzDGUSndlLM6wXESlKzz//PO3bt2fKlCncfvvtTJo0iTZt2sQdlogkySe5uoMwLU1/M/uCMLwBZnZfdP984CV3vy+fANz9v8B7wOpA19T1Ub+oTQj9ud7IYX/LgGeju/ul2d9qlNfATc0nVhGRuCSG7Nt+++3p0qULH374IX369GGVVYplLGgRScj5XRklLQcBtxCa09oTOqOfCGwcLT+8knEMim6vN7OtEgvNbANCUgcwOLkzvZmdY2afmlm6mrJBQBnQx8z2S9qmAXA98CfC1YNPVjJeEZEasWzZMgYOHMhhhx2Gu9OyZUsee+wxNtssW08KEYlTXj953H2xu19IGO6gE3AwsCvQ3N3/7u5LKxOEu48GRhCStg/MbIyZPQF8QbhK8SnCBM7J1gPaEJoiU/c3AzgPWJMwrtWbZjaaMFTE+cBcoGuWQUtFRGI3bdo0dtppJwYMGMBaa63FwoUL4w5JRHJQqfpkd1/m7u+6+4vu/pa7L65qIO5+NmFanfeAvQm1ZF8SprI5xt1X5Lm/2wjjWI0FtiLUqq1KuIKwxN0rbGIUEYnD4sWLueSSS+jUqRM//vgjjz/+OI899hiNGjWKOzQRyUHOQzGY2QTgPsLcgdXy88ndRwGjcix7JXBlBWUmAhOrGJaISI1atmwZDz/8MCeeeCJDhgzR1DUitUw+NVf7EJKrH83sn5kG4BQRkfzNnz+fq6++msWLF7P22mszY8YM7r33XiVWIrVQPsnVn4HrgNnAScDLZva1mV2V3AldRETyM27cONq3b8+VV17JhAkTAGjWrFnMUYlIZeVzteBX7n65u29J6Mv0AGEcqcuAz8zsVTM7zcyaVFOsIiJ1yi+//MJJJ51Ely5daNSoEa+//romWhapAyrboX2iu58CbEioxXoF2I3QWfz7woUnIlJ3nXjiiYwaNYrLLruMadOmseuuu8YdkogUQJVGn3P3Re7+IHAYcDGwHMg0yrqISL33ww8/8OuvvwJwww03MHXqVK6++mrWWGONmCMTkUKpUnJlZrub2UjgB2AwYWLntwoRmIhIXeLu/POf/6Rt27a/zwXYrl07tt9++5gjE5FCy3kohgQz24wwKvuJwJaEUdq/IwwCer+7f1bQCEVEarmvv/6aXr16MX78ePbaay/69esXd0giUo3yGefqZEL/qj0JNV6LgEeA+wlzCnp1BCgiUps9++yzdO/enQYNGjBixAh69eql+QBF6rh8aq7ujW6nEBKqR919XuFDEhGp/dwdM2OHHXbgsMMO46abbqJVq1ZxhyUiNSCf5GogcJ+7f1ldwYiI1HZLly7l+uuv5/XXX2fs2LG0bNmSRx99NO6wRKQG5TPO1aVKrEREMps6dSqdOnXi8ssvZ5111tFEyyL1lBr+RUSqaNGiRfTt25edd96Z2bNn8/TTT/Ovf/1LEy2L1FMZmwXN7H3AgcPd/Zvofi7c3XVtsYjUG8uXL+ff//43p512GjfccIOmrhGp57L1udqYkFw1SLkvIlLvzZs3j5tvvpn+/fuz9tpr8/7779OkiWb/EpEsyZW7r5ftvohIffX888/Tu3dvvv/+e3beeWe6dOmixEpEfqc+VyIiOZo9ezbHH388hx56KE2aNGHKlCl06dIl7rBEpMjknFyZ2TNmdmwFZbqZ2TNVD0tEpPicdNJJPProo1xxxRW899577LzzznGHJCJFKJ9xrg4DplZQ5s/AoZUPR0SkuHz33XesueaaNG/enBtvvJHBgwez7bbbxh2WiBSxQjcLrgksL/A+RURqnLtz11130bZtW/r27QtA27ZtlViJSIXyTa4yXi1oZs2BA4HvqxSRiEjMvvrqK/bbbz969epFhw4duOSSS+IOSURqkazNgmaWOnfgADO7KE3RBoRaK4AhhQhMRCQOY8aMoXv37qy22mqMHDmS008/HTOLOywRqUUq6nP1OeW1VTsCv5C+ZmpFtG4CcFvBohMRqSFlZWWsssoqdOjQgSOPPJIbbriBTTbZJO6wRKQWyppcuXvHxP9mVgbc6e5XV3tUIiI1ZOnSpQwaNIjXX3+dcePGsfHGGzNq1Ki4wxKRWiyfqwW3BX6urkBqs5EjR670YTx9+nRKSkpijEhEcvHOO+9w6qmn8uGHH9KjRw8WLVqk+QBFpMpy7tDu7h+5+6zqDKa2GjVqFNOnT//9fklJCT169IgxIhHJZtGiRVx44YXssssuzJkzhzFjxvDwww8rsRKRgsg2cfMF0b/3uvtvSfcr5O43VzmyWqakpISJEyfGHYaI5GDFihU88cQTnHHGGVx//fU0bdo07pBEpA7J1ix4E6Ez+7PAb0n3K7psxoF6l1yJSHGbO3cuN910EwMGDKBx48bMmDGDtddeO+6wRKQOypZc/SW6/W/KfRGRWmXMmDGceeaZ/Pjjj+y222506dJFiZWIVJuMyZW7P5ftvohIsZs1axbnnnsujzzyCNtuuy1PP/00HTt2rHhDEZEqKPT0NyIiReOkk07i8ccf5+qrr2bq1KlKrESkRuQ8FIOZtQT+BEx194XRslWAvwNHAguAm9z9xeoIVEQkF//73/9o2LAhzZs3Z8iQIZSVldGuXbu4wxKReiSfmqurgKdZeWLmfsD1wK7A/sCzZrZj4cIrXiNHjqRz58507tx5pWEYRCQeZWVl3HnnnbRt25aLLgqzdG2zzTZKrESkxuWTXO0GTHD3pQAWJts6F/gKaAvsCywBch6yoTZLHttK41qJxOuLL75g33335cwzz2SnnXbi0ksvjTskEanH8hmhvQUwNun+dsCGwGB3/xT41MyeJtRi1Qsa20okfs888wzdu3dnjTXW4O677+bUU0/VRMsiEqt8aq7WAJYl3d+dMKbVhKRl3wAbFSAuEZGsysrKAOjYsSPHHHMMH3/8MaeddpoSKxGJXT7J1f8I8wsmdAF+dfcPk5atB5QWIjARkXSWLFnC5ZdfzkEHHURZWRkbb7wxDz30EBtvvHHcoYmIAPklV+OAA83sSjPrDxxMGL09WRvg20IFJyKS7I033mCHHXbgmmuuYaONNmLx4sVxhyQi8gf5JFeDgR+By4GBwK/AlYmV0VANewCTCxifiAgLFy7kvPPOY/fdd6e0tJTnn3+eBx54gLXWWivu0ERE/iDnDu3u/oOZtQUOjRa95O6/JBVZD7iCMFyDiEjBlJWVMWbMGM466ywGDx6sqWtEpKjlc7Ug7j4feCTDuhnAjEIEJSLy22+/ceONN3LZZZf9PtFy48aN4w5LRKRClZr+xszWNbN9zOwoM9vXzNYtdGAiUn89/fTTtG3bluuvv55JkyYBKLESkVojr+TKzFqY2WjgJ2A8MBp4CfjZzEabWZ0ehkGjsotUr59++onu3btz5JFHssEGG/DWW29x0EEHxR2WiEheck6uzGw94HXgaGAW8CRwR3T7U7R8SlSuTtKo7CLV65RTTuGpp57iuuuu45133qFDhw5xhyQikrd8+lwNALYArgWuc/cliRVmtjpwCeFKwgHA+YUMsphoVHaRwvr2229Za621WG+99bj55ptxd7bZZpu4wxIRqbR8mgUPB15x98uTEysAd1/q7lcCrwBHFDA+EamjysrKGD58OO3ataNv374AbL311kqsRKTWyye5agm8WUGZNwENkywiWX322WfsvffenHPOOey6665cfvnlcYckIlIw+TQLzgc2qaBMy6iciEhaTz/9NN27d6dhw4b885//5KSTTtJ8gCJSp+RTczUF6GpmO6RbaWbbAV2jciIiK1mxYgUAO+20E927d+eTTz7h5JNPVmIlInVOvtPfrAq8YWYjzKybme1pZl3NbDihSXC1qFylmFkPM5tsZnPNrNTMpppZHzOr1HhcKfvuZWYe/d1e1f2JSG4WL17MpZdeyoEHHkhZWRkbbbQR999/Py1atIg7NBGRapFz0uLubwDHA0uB3sC/gImEEdvPApYDJ0bl8hYlaA8DHQnzE74EtAZuB0abWYPK7Dfa92bATYBXdh8ikr8pU6awww47cN1119GqVStNtCwi9UK+09/828xeIjT/7Qg0BeYC04B/u/ucygRhZscAZxMmht7L3b+Ilm9IuALxKOAcYGgl9m3APYRE8gHgpMrEKCK5W7BgARdffDG33347rVq1Yty4cRoMVETqjbySK4AogRpZ4Dgujm77JRKr6Fg/mdlZhBqy/mZ2m7uX5bnvM4H9gHOB5oUIVkQq9txzz9GnTx8GDhyoiZZFpF7JqVkwmkNwkJkNNLOCjmNlZpsAHQjNjY+lrnf3ScB3QAtglzz3vQVwA2FkefWzEqlGv/76K/3792fRokU0atSI999/n9tuu02JlYjUO1lrrqKR18cCnVOWvwJ0cfdlBYghcfXhR+6+KEOZdwjDPOxAjlcjRs2B9xIe42nu7roqSaR6PP744/Tp04fZs2ezzz77cNBBB9GoUaO4wxIRiUVFNVd/A/YBfgMeInQ4/y1adm6BYtgiuv0mS5lvU8rm4hxCUnilu39WibhEpAI//PADxxxzDH/961/ZeOONmTp1qvpWiUi9V1Fy1R2YB5S4+0nufiKhI3tptK4QGke3C7KUKY1uc2pfMLM/AYOAdwlXCeYlGrZhqplNnTVrVr6bi9Qbp556Ks899xyDBw/m7bffpqSkJLBH00sAACAASURBVO6QRERiV1GH9jbAaHf/b2KBu39jZk8ARxcohkRbXUGGSUhqDlwdONXdV+S7D3cfSdRpv2PHjhq+QSTJzJkzady4Meuttx633norAG3atIk5KhGR4lFRzVVj4L9pln9LeY1TVSWmy8m2v8S6XKbWORfYCxjk7u9XJTARKVdWVsZtt91G+/btf59ouU2bNkqsRERSVFRzZUC6oQ/yHQ4hm5nR7WZZyrRKKZvNUdHtAWa2d8q6zRNlzKw9UOruh+WwT5F67ZNPPuH0009nypQpHHzwwVxxxRVxhyQiUrRyGedqYzPbMXUZQDTP4B8uwXP39/KIYVp0287MGma4YrBTStlc7Jpl3cbR39w89idSLz355JMce+yxNG7cmAceeIDjjz9e8wGKiGSRS3J1RvSXyoCpaZZ7jvsNhd3/a2bvETrKdyWMol5+kFD7tAlh9PYKp9Zx986Z1pnZlcAVwHB3PyfXGEXqoxUrVtCgQQN23XVXevToweDBg9lwww3jDktEpOhVlAS9R83MxzeIMIDo9WY2xd2/BDCzDYA7ojKDk0dnN7NzCMMtvB1dxSgiBbBo0SKuuuoq3nrrLSZMmECLFi345z//GXdYIiK1Rtbkyt071kQQ7j7azEYQJoD+wMzGA8sI09Y0AZ7ijyOsr0e4mvHHmohRpD6YPHkyp59+Op9//jmnnnoqixcvZq211oo7LBGRWiWn6W9qgrufDfQk1JbtDRwEfEmonTqmMkMqiEhuSktL6dOnD3vttRdLly7lpZde4p577lFiJSJSCeauYZyy6dixo0+dGrqWde7cGYCJEyfGF5BINViwYAElJSUcdthhXHvttZq6RkTqNDN7tzpb54qm5kpEatYvv/xC3759f59oecaMGdxyyy1KrEREqkjJVQU+++wzOnfuTOfOnZk+fXrc4YhUmbvz2GOP0bZtW2655RYmT54MoCZAEZECUXJVgUWLyofdKikpoUePHjFGI1I1P/zwA0cffTTdunWjVatWvPvuuxx44IFxhyUiUqeoz1UF1l57bZ8/P5dZd0SKX5cuXZg4cSJXX301559/PquumvOQdCIidUZ197lSclUBJVdS23399dc0btyY9ddfn88//xyA1q1bxxyViEh81KFdRCplxYoV3HrrrbRv355+/foBIalSYiUiUr3ybhMws62AY4FtgEbufmS0fBNgO+A1d59X0ChFJC8ff/wxp512Gm+++SaHHHIIV111VdwhiYjUG3klV2bWF7g2abvkNsWGwBjCoJ8jChKdiOTtiSee4LjjjmPttdfmoYceokePHppoWUSkBuXcLGhmRwGDgSnAHsCQ5PXu/gUwDTiikAGKSG6WL18OwG677cYJJ5zAxx9/TM+ePZVYiYjUsHz6XJ0PzAQOdvcpQGmaMh8R5vsTkRqyaNEi+vbty/77709ZWRktWrTg7rvvZoMNNog7NBGReimf5KoEGOvui7OU+R7YsGohiUiuJk2axHbbbceNN95I69atWbJkSdwhiYjUe/kkVw2ApRWUWS+HMiJSRaWlpZx11ll07tyZsrIyJkyYwMiRI2nYsGHcoYmI1Hv5JFdfAbtkWmmhY8duwCdVDUpEsjMzJkyYwAUXXMAHH3zAvvvuG3dIIiISySe5Gg3sZGZnZlh/HrA18GiVoxKRP5g9ezYXXnjhShMtDxkyRHMCiogUmXySqyHAp8BwM5sA7AdgZldG928CpgN3FDxKkXrM3XnkkUfYZpttGDZsGK+99hqAmgBFRIpUzsmVuy8A9gaeBDoThmMw4HJgH+Ap4AB3V58rkQL57rvvOOKIIzjuuOPYYostePfddznggAPiDktERLLIaxBRd58N/NXMWv5/e3ceX0V1/nH88wCBBNl3cUcpylJcEKWKIIsLVatFalUEIQgKKGrr0gK/2iJFa/VXLAgVZRF3Ra3WH1WpQbCg7CibFZVFRdQiIhBIgOf3x0xCCDfJTXJv7k3yfb9e8xqYOXPuM/fc5D6ZOXMOQf+rhsD3wHvuvjEO8YlUaoMGDeKdd97hwQcfZMSIEVStWjXRIYmISBE0cXMRNHGzlLVPPvmE2rVr06RJEz7++GOqVKnCiSeemOiwREQqDE3cLFJJ7N+/n4ceeoh27dpx9913A9CyZUslViIi5UzUtwXN7OEoi7q7jyhhPCKV0qpVq0hPT2fRokVceumljBkzJtEhiYhICRWnz9XwIvY7QQd3B5RciUTpxRdf5JprrqFu3bo888wzXHXVVZoPUESkHCtOctWugO31gDOBu4EM4N7SBiVSGezbt49q1arRuXNnBgwYwNixY2nUqFGiwxIRkVKKWYd2M2sBrARudPenYlJpElCHdom13bt3M3r0aJYsWUJGRgZVqqjro4hIWSo3Hdrd/VPg78CvYlWnSEWTkZFBu3bteOihh2jdurUmWhYRqYBi/SfzFoIpcEQkjx9++IHBgwfTrVs3qlSpwty5c5k0aZJGWRcRqYBillyFEzefB+yMVZ0iFUXVqlWZN28ed9xxBytXrqRLly6JDklEROKkOEMxnF5IHccA6UAHYEYM4hIp977++mvGjRvH2LFjqVmzJitWrCA1NTXRYYmISJwV52nBJQTDLBTEwjJ3lCoikXLO3Xn66acZMWIEO3bsoFevXvTs2VOJlYhIJVGc5OohIidXB4DvgEVAhms+HanENm/ezE033cTrr7/O2WefzWOPPUabNm0SHZaIiJShqJMrd/91PAMRqQgGDx7MvHnz+Mtf/sLw4cM10bKISCUU9ThX4fQ3a919UnxDSi4a50qKsn79eurUqUOTJk1Yv349VapUoUWLFokOS0RECpBM41wNAY6LVyAi5c2+fft44IEHaNeuHXfddRcAJ510khIrEZFKrjh9rjYBDeMViEh5snLlStLT01m6dCk/+9nPGDt2bKJDEhGRJFGcK1fPAReaWe14BSNSHrzwwgt06NCBzZs38/zzz/Pyyy/TvHnzRIclIiJJojjJ1b3Af4C3zKyrmR0Rp5hEklJ2djYA5513HoMGDWLNmjX06dOHYPxcERGRQHE6tO8gSMZqcnBIht0cPjyDu3vdmEWYYOrQLrt27WLkyJEsXbqUuXPn6glAEZFyLt4d2ovT5+o/FD6IqEiFM2fOHG644QY2bNjAsGHDyM7OVnIlIiKFKs44V3HL8ESSzY4dO7jtttuYOnUqP/rRj5g3bx6dO3dOdFgiIlIOFNrnysz6mdmPyyoYkWRRrVo1FixYwN13383KlSuVWImISNSK6tA+Hbi8DOIQSbitW7cyYsQIdu/enTvR8rhx4zQnoIiIFEtxnhYUqZDcnZkzZ9K6dWsmT57MwoULAahRo0aCIxMRkfJIyZVUaps2baJXr17069ePVq1asWLFCrp3757osEREpBxTciWV2uDBg5k/fz4PP/ww8+fP55RTTkl0SCIiUs5F87RgPTM7tjiVuvumEsYjEncfffQR9erVo2nTpkycOJGqVaty/PHHJzosERGpIAodRNTMDlD8sa3c3YszflZS0yCiFUd2djYPPvgg99xzD1dffTXTpk1LdEgiIpIAyTCI6A5ge7wCECkLK1asYODAgSxfvpzevXszbty4RIckIiIVVDTJ1f+6+x/iHolInLzwwgtcffXVNGrUiBdffJHevXsnOiQREanA1KFdKqyciZa7dOnCjTfeyJo1a5RYiYhI3CVVcmVm15jZfDP73sx2mtkSMxtmZlHHaWZVzOwnZnZvWNfnZpZlZlvN7P/MTIOiVnA7d+7klltu4fzzz2f//v00adKECRMm0KBBg0SHJiIilUDSJFdmNhF4CugAzAfeAn4ETABeNLNoZ8ttAfwbGAm0BlYDLwEbgYuBl81smplZbM9AksGbb75J27ZtmTBhAqeffnru1SsREZGykhTJlZn1BoYCXwE/dvdL3P0KoCWwFrgCGB5ldQ68TZBINXH3C939l+7eEegK7AKuDxepIHbs2MGAAQO48MILSU1NzR27SlPXiIhIWSt0KIYyC8JsCXAG0N/dn8i3rwswlyDxOsrdD5TytUYBY4C33b3Iobg1FEP5kJmZyRlnnMEVV1zB6NGjlVSJiEiB4j0UQ8KvXJnZ0QSJVRbwQv797v4O8AXQDDg7Bi+5PFwfHYO6JIG2bNnC8OHD2bVrF2lpaSxfvpyxY8cqsRIRkYRKeHIFnBauV7t7ZgFlFucrWxotw/WWGNQlCeDuTJ8+ndatW/PYY4/x3nvvAZpoWUREkkMyJFcnhOuNhZTJmU7nhELKFMnMagK3hP+dVZq6JDE2bNjAhRdeyIABA2jbti0rV67URMsiIpJUkiG5qhWudxVSZme4rl3K13qEIEFbAzxayrokAW688UYWLlzIxIkTeeedd2jVqlWiQxIRETlEMswBmDMkQlx71pvZaKA/8D3wC3ffW0jZwcBg0K2mZLBu3Trq1atHs2bNeOSRR6hWrRrHHlusucRFRETKTDJcucp5FK9WIWVy9pXosT0zux34A8EVsIvdfXVh5d39UXfv4O4dUlJSSvKSEgPZ2dn88Y9/pH379vzmN78BoEWLFkqsREQkqSXDlasN4fq4Qsock69s1MzsZuBBIBO4xN0XFrcOKXvLli0jPT2dFStW0KdPH+67775EhyQiIhKVZLhylTM0QhszSyugzJn5ykbFzIYBDwN7gMvCYR0kyT333HN07NiRr776ipdeeonnn3+epk2bJjosERGRqCQ8uXL3zcAyoDrQJ//+cBDRowkGEY36qpOZ3Ugwdc5e4HJ3nxOTgCVusrKyAOjWrRtDhw5lzZo1XHHFFQmOSkREpHgSnlyFxoXr+83spJyNZtaE4Ak/gPvyjs5uZsPNbJ2ZHTKie7jvhvC4vcDP3f2N+IUupbVjxw6GDRtG165d2b9/P40bN+bhhx+mfv36iQ5NRESk2JKhzxXu/qKZTQJuAj40szlANtAdqAO8QnAVKq9GQCuCK1q5zOxU4G8ETyF+BvzCzH4R4WW/dfdfx/REpNhmz57NkCFD+PzzzxkxYgTZ2dlUrRrtHN0iIiLJJymSKwB3H2pm7wLDgC5AVWAdMBWYVIw5BetxcHiHk8Mlko2AkqsE+f7777n55puZOXMmrVu3ZsGCBZx9dixmNxIREUmspJi4OZlp4ub4yMzM5Mwzz6R379789re/1XhiIiJSZuI9cXPSXLmSiu/LL7/k3nvv5YEHHuCII45g2bJlVK9ePdFhiYiIxFSydGiXCszdefzxx2ndujXTpk3j/fffB1BiJSIiFZKSK4mrTz/9lJ49ezJo0CDat2/PBx98QLdu3RIdloiISNwouZK4Gjp0KIsWLWLSpElkZGTQsmXLRIckIiISV+rQXgR1aC++1atX07BhQ5o1a8Znn31GtWrVOOaYY4o+UEREpAzEu0O7rlxJzGRlZTFmzBhOO+00fvvb3wJwwgknKLESEZFKRU8LSkwsXryY9PR0PvzwQ66++mruv//+RIckIiKSELpyJaX27LPPcvbZZ7Nt2zZeffVVnn76aRo3bpzosERERBJCyZWU2N69ewHo3r07t9xyC6tXr+bSSy9NcFQiIiKJpQ7tRVCH9sPt2LGDO++8k5UrV/Luu+9qLkARESlX1KFdksrrr79OmzZtmDJlCj/5yU/Izs5OdEgiIiJJRcmVROX777+nb9++XHLJJdStW5cFCxbw4IMPkpqamujQREREkoqSK4lK9erVWblyJffccw/Lli3jrLPOSnRIIiIiSUnJlRToiy++YMiQIezcuZO0tDSWLVvG7373O80JKCIiUgglV3KYAwcO8Oijj9K6dWtmzpzJ4sWLAUhJSUlwZCIiIslPyZUcYv369XTv3p0hQ4Zwxhln8OGHH3L++ecnOiwREZFyQyO0yyFuvvlmli1bxpQpU0hPT8fMEh2SiIhIuaJxropQGca5WrVqFQ0bNuTII49kw4YNpKSkcNRRRyU6LBERkbjQOFcSN1lZWdxzzz2cfvrpuRMtH3/88UqsRERESkG3BSupRYsWMXDgQFavXs21117LAw88kOiQREREKgQlV5XQM888Q9++fWnevDn/+Mc/+OlPf5rokESkktqxYwdff/21ZnuQUqtWrRqpqak0btw44QNcK7mqRPbs2UNqaio9evTg1ltv5Xe/+x116tRJdFgiUknt2LGDrVu3ctRRR5GWlqYHaKTE3J19+/axc+dONm3aRNOmTalbt27C4lGH9iJUhA7t27dv54477uCDDz5gwYIFmmhZRJLC+vXrad68OTVr1kx0KFKB7N69m6+++ooWLVoUWEYd2qVUXn31Vdq0acPUqVPp2rUr+/btS3RIIiIAZGdnk5aWlugwpIJJS0tj7969CY1ByVUFtX37dn75y1/ys5/9jIYNG/L+++9z//33U6NGjUSHJiKSS7cCJdaS4TOl5KqCSk1NZe3atYwZM4YlS5bQoUPcrn6KiIhIHkquKpDNmzdzww03sHPnTlJTU1m6dCmjRo3SRMsiIiJlSMlVBXDgwAEmT55MmzZtePrpp1myZAkQPJYqIiJlb9CgQZgZt99+e8T9119/PUcffXTEfXPnzsXMmDNnziHbs7OzeeSRRzjnnHOoV68eNWrU4IQTTmDgwIEsW7Ys5udQmClTpnDyySdTo0YNWrVqxeTJk6M+dsaMGZxxxhnUqVOHxo0b07NnT+bPn39YuYyMDM4991zS0tJo0KAB1113HVu3bo3lacSNkqty7uOPP6Zbt27cdNNNnHXWWaxatYquXbsmOiwRkUorMzOTF154AYCnnnoqJg8S7dq1i+7du/OrX/2Kjh078tRTT/Hmm28yatQoPvvsM7p3717q14jWlClTGDJkCL179+af//wnffr0YejQoUyaNKnIYx999FGuv/56OnbsyKxZs3jsscfIysqiZ8+eLF++PLfc/PnzueCCC6hXrx6zZs1i/PjxzJs3j+7duye8s3pU3F1LIUutWrU8mV100UVet25df/zxx/3AgQOJDkdEJGpr1qxJdAhx8dRTTzngvXr1csBfe+21w8r079/fjzrqqIjHZ2RkOOBvvfVW7rb09HSvXr26L1iwIOIxL730UmyCL0J2drY3btzY+/Xrd8j2AQMGeMOGDT0rK6vQ4zt16uSdOnU6ZNuOHTs8JSXF77777txt3bt39xNPPNGzs7Nzty1atMgBnzhxYpFxFvXZApZ4HHMHXbkqh1auXMmWLVsAmDx5MmvWrGHgwIFJ8YSEiEhlN2PGDOrXr8/06dNJS0vjiSeeKFV9W7ZsYfr06dxwww106tQpYpkrrriiVK8RrYULF/LNN9/Qt2/fQ7Zfd911/Pe//+Xdd98t9PisrKzDBq+uWbMmKSkpHDhwIHfbe++9R8+ePQ/p3nLmmWfSsGFDXn755RicSXwpuSpH9u7dy+jRo+nQoQMjR44E4LjjjqN58+YJjkxERAC+/PJL5syZw1VXXUXjxo25/PLLefXVV/nuu+9KXGdGRgb79+/nsssuK3EdHo5gXtSSN8GJZPXq1QC0bdv2kO1t2rQBYM2aNYUeP3ToUObMmcPjjz/O9u3b+eKLLxg+fDgpKSmkp6fnlqtatWrEh7Fq1KjBqlWrojrnRFKP53Ji4cKFpKens3btWvr166eJlkWkQrr11ltZsWJFQmM49dRT+ctf/lKiY2fOnMmBAwfo168fAP379+eZZ57hueee48YbbyxRnZs3bwaCP6ZLasaMGQwYMKDIcv3792f69OkF7t+2bRsA9evXP2R7gwYNDtlfkIEDBwJBkjVo0CAAmjVrxltvvcWPfvSj3HKtWrXivffeO+TYjRs3smXLFlJSUoo8j0RTclUOPP300/Tt25djjjmG2bNnc9FFFyU6JBERieCJJ56gZcuWubfvevToQfPmzXniiSdKnFzFwqWXXsrixYuLLNeoUaNC93s4ZV5Ju6H8/e9/Z9iwYQwZMoTLLruMzMxMxo8fT69evcjIyMi9IjZixAj69u3LqFGjuOWWW9i2bRuDBw+mSpUqVKmS/DfdlFwlsczMTNLS0rjgggv49a9/zejRo6ldu3aiwxIRiZuSXjFKBosXL2bNmjXcddddbN++PXf7z3/+cyZMmMB//vOf3Ksz1apVY//+/RHrydme09/omGOOAYIrN61atSpRbA0aNIhqIuOiEpe8V6iOPPLI3O05V6xy9kfi7gwePJgrr7yS8ePH526/4IILOPnkkxk9enRuf6prr72WdevW8ec//5mxY8diZlx11VX06tWrXNwWTP70rxL67rvvSE9Pz50LsFGjRvzpT39SYiUiksRmzJgBwP3330/9+vVzlwkTJgAc0rG9SZMmfPvtt2RlZR1Wz5dffglA06ZNAejatStVq1bltddeK1VsKSkpRS45t+0KktO3KqfvVY6cvlatW7cu8NitW7fy9ddfc+aZZx6yvXr16rRv3561a9cesn3MmDF8++23fPDBB2zZsoVnnnmGjz/+mHPPPTfq804UXblKMi+//DJDhw7lm2++4Y477mD//v0aDFREJMllZWXx7LPPctZZZ3Hfffcdtv+2225j5syZjBkzBjPj/PPPZ9y4cbz66qtceeWVh5SdNWsWRx55ZO5VqubNm3P99dfz6KOPcs0110R8YvCVV17h8ssvLzC+WN0W7NSpE40aNeKpp56iR48euduffPJJGjRowDnnnFPgsfXr16dGjRosWrTokO1ZWVmsWLGCFi1aHHbMEUccQbt27QD45z//ybp163j88ceLPI+Ei+c4DxVhKatxrrZt2+Z9+vRxwE899VRfunRpmbyuiEiiVKRxrmbNmuWAT58+PeL+SZMmOeBvv/22u7sfOHDAe/bs6UcccYSPGTPG33zzTZ81a1bu98C0adMOOf6HH37wzp07e1pamt92223++uuv+zvvvOPTpk3zHj16eL169eJ9ioeci5n5yJEjPSMjw0ePHu1m5hMmTDik3MCBA71q1aqHbBs+fLgDfvPNN/sbb7zhr7zyivfo0cMBf+WVV3LLLVu2zMeOHeuzZ8/22bNn+8iRI7169ep+5513RhVjose5SnjykuxLWSVXmZmZ3r59ex87dmyRg7CJiFQEFSm5uuyyy7x27dq+a9euiPu3b9/uaWlp3r9//9xtmZmZPnLkSG/ZsqVXr17da9Wq5eeee+4hSUZeWVlZPmHCBO/UqZPXrl3bU1JS/Pjjj/f09HRfuXJlPE6rQJMnT86N+6STToo4sGf//v09uIZzUHZ2tv/1r3/19u3be61atbxRo0bepUsXf+ONNw4pt2rVKj/nnHO8bt26npqa6qeddppPnTo16vgSnVxZ8BpSkNq1a/sPP/wQl7o3btzI73//e8aPH0/t2rXZt2+fbgGKSKWxdu1aTjnllESHIRVQUZ8tM1vq7h3i9frq0J4ABw4cYOLEibRt25bnn38+d8JNJVYiIiLln5KrMvbRRx/RpUsXhg8fzjnnnMPq1avp0qVLosMSERGRGNGlkjJ26623snr1aqZPn06/fv00H6CIiEgFo+SqDKxYsYImTZrQvHlzJk+eTI0aNWjWrFmiwxIREZE40G3BONqzZw8jR46kQ4cOjBo1CgjmhlJiJSIiUnHpylWc/Pvf/yY9PZ2PPvqIAQMG8OCDDyY6JBGRpOPu6h4hMZUMoyDoylUcPPnkk3Tu3Jk9e/bwxhtvMHXq1MNmEBcRqexSUlLIzMxMdBhSwWRmZlKjRo2ExqDkKoZ2794NwMUXX8ydd97JqlWruOCCCxIclYhIcmrSpAlffPEFu3fvToqrDVJ+uTvZ2dls27aNzz//nIYNGyY0nqQaRNTMrgFuAn4MVAXWAdOASe5+oAT1XQTcDnQAUoFPgWeAP7v73mjqiGYQ0W3btnHbbbexZs0aFi5cqPGqRESitGPHDr7++muys7MTHYqUc9WqVSM1NZXGjRuTmppaaNl4DyKaNFmAmU0EhgJ7gH8B2UB3YALQ3cz6uPv+YtR3J3A/sB+YC3wHdAHuBS4xs+7uvru0cb/44osMGzaMbdu2cdddd3HgQLFzQBGRSqtOnTrUqVMn0WGIxFRSJFdm1psgsfoKOM/dPw63NwUygCuA4cD4KOvrANwH7Aa6ufv74fZawOvAecBY4LaSxrxt2zZuuOEGXnrpJU4//XTeeOMNTj311JJWJyIiIhVEsvS5+k24visnsQJw960EtwkB7jazaOO9GzDg/pzEKqxvJzAAOAAMNbN6JQ24Zs2afPbZZ9x33328//77SqxEREQESILkysyOBs4AsoAX8u9393eAL4BmwNlR1FcduDj871MR6vsUWAhUB3oVJ9YNGzZw/fXX88MPP5CamsrixYu566671MdKREREciU8uQJOC9er3b2gZ3IX5ytbmFZATWCbu39S2vpq1qzJ/v37efjhh2nbti2zZs1i+fLlAFStWjWKcERERKQySYZLLieE642FlNmUr2w09W0qpEzU9TVu3JjzzjuPBQsWcNFFF/G3v/2NY489NoowREREpDJKhuSqVrjeVUiZneG6dlnXt3nzZrZs2cITTzxB3759NZKwiIiIFCoZkqucbCVWA26Vuj4zGwwMDv+7F1jVr18/+vXrV9rYJDYaAd8mOgjJpfZIPmqT5KL2SD6t4ll5MiRXOSN01iqkTM6+wkfzjFF97v4o8CiAmS2J50BjUnxqk+Si9kg+apPkovZIPma2JJ71J0OH9g3h+rhCyhyTr2w09RXWMao49YmIiIhELRmSq+Xhuo2ZpRVQ5sx8ZQuzDsgEGpjZiQWU6ViM+kRERESilvDkyt03A8sIxp3qk3+/mXUBjiYYvX1hFPVlAbPD/14bob4WQCeCcbVejyLER6MoI2VLbZJc1B7JR22SXNQeySeubZIUEzeb2ZUEA4h+BXR29/Xh9iYE09+0Bm519/F5jhlOMCXO4GpIMwAAEZ9JREFUInfvl6++M4H3Ca5gne/ui8LttYB/EMwx+Bd3L/H0NyIiIiKRJPzKFYC7vwhMIhiF/UMze83MXgI+JkisXiGYwDmvRgS9/Q/rW+XuiwmmwKkJLDCzN83seeATgsTqfWBknE5HREREKrGkSK4A3H0owW28ZQQJ0IXAeoKrU73dfX8x6/sTwTQ4GQR9ti4FsgkGKz0F+NrMlpjZsGLMWXgIM7soTNy2mdluM1tlZiPNrEZJ6quMzOwaM5tvZt+b2c6StImZVTGzn5jZvWFdn5tZlpltNbP/M7PL43kOFUks2qOQugebmYdL/j+WpACxbhMzq2pmQ8xsnpn918z2mNnm8I/aS2Mdf0UTy/Yws/pm9kcz+9DMdpnZXjPbaGYzzUwT1hbBzFqZ2Qgze9LM1pnZgfD3y5WlrLf0bezulWIBJhKMfZVJcGvwZWBHuO0loGox67szPHYfMIfgtubX4baFQM1En3OyL7FqE+Ck8BgH/gu8ATwLLMqzfRrhbXAt8W2PAuo+LqzrQFjfhESfb3lY4vB7qwHBlXsHthP0O30W+DewG3gs0eeczEss24PgrsvG8NhvwvpeJLio4AQXA3on+pyTeQH+kud3fN7lykS3ccLfnDJqgN7hG7MFaJlne1NgTbhvRDHq6xB+SewCzsqzvRbwTljf/yb6vJN5iWWbACcC/wIuyv/BJ7gKujOsb0CizztZl1j/jOSr2wj+ANkJTFdylZg2IbhT8e/wuCnAEfn21wLaJvq8k3WJQ3s8HR7zOnn+GA/b6Z5w37dASqLPPVkXYBDwJ+AX4ffA3NIkVzH9Xkr0m1NGDbAkfFP6RdjXJc+bWSXK+l4Mj/mfCPtaAPsJRnavl+hzT9Yl1m1SxGuNCuv7V6LPO1mXeLYHcFN4/M15vjSUXJVxmwBDwmPmoqu4ydAeW8Jjzo6wryrBlUQHWif63MvLEoPkKmZtnDR9ruLFzI4GziAYeuGF/Pvd/R3gC4LO9GdHUV91gr5cAE9FqO9TgtuC1YFeJQ68Aot1m0QhZzyzo2NQV4UTz/YwsxMI/rL8N4c/lCIFiFObDA/X93v4bSHRiVN77I2ynKbNKQOxbuMKn1wBp4Xr1e6eWUCZxfnKFqYVwVOI29z9kxjUVxnFuk2K0jJcb4lBXRVRXNrDglnOpxJMs5WuL/RiiWmbmFkzoC1BP54MM2tnZveY2d/CDtU9Sx9yhRaPn5F/hutRZlYzZ2P4c/M/QBrwqrt/XdxgpURi2sbJMLdgvJ0QrjcWUmZTvrLR1LepkDLFqa8yinWbFCj8pXVL+N9ZpamrAotXewwHugJ3u/tHJYirMot1m/w4XG8ARgO/4eAk9wC/MbN5BB2odaXkcPH4GRlF8CX9U2Cjmb1HcDWrPcEDIE8CQ4sfqpRQTNu4Mly5ypmkeVchZXaG69oJqK8yKsv38BGCH4Q1aJTkgsS8PcKpp8YBS4E/lzy0SivWbdIgXJ8A/BaYSTAkTR2gG7AWOA94vtiRVg4x/xkJk9huwAyCcRsvIehQfRLwKfCOu/9QomilJGLaxpUhucr56yxWtyRiXV9lVCbvoZmNBvoD3wO/cPdo+zhUNjFtjzy3A6sDA72YY9QJEPufkZzf9dUIHuzo7+7r3P0Hd88ALiCc0SKcckwOFfPfWWZ2MkF/0AuB64AjgXpAd4Iv+ClmNjVWrydFimkbV4bkKifzr1VImZx90fyVEOv6KqO4v4dmdjvwB4K/NC5299UlqaeSiHV73EJwFWScu39QmsAqsXj93oIIV3Dd/XMOzrXaPYr6KpuYtoeZVSPopnAS8HN3f9Ldv3L37939baAnsBUYYGbnlyJuiV5M27gy9LnaEK6PK6TMMfnKRlPfYdPulLC+ymhDuI5VmxzCzG4GHiT4S/wSdy9ywu9KbkO4jlV7XBGue0a4CnJ8ThkzawvsdPdLoqizstkQrmP9ewvgswLK5GxvFkV9lc2GcB2r9jiLYGq3TyP9fnL3bWY2G7ge6EEw04jE14ZwHZM2rgzJVc5j+G3MLK2ApwDOzFe2MOsIvrQbmNmJBTwx2LEY9VVGsW6TXGY2DHgY2ANcFj4+K4WLV3t0KmRf83D5vhj1VSbx+L21CzgCaFhAmUbhemcB+yuzWLdHzh/nhX3+t4frBoWUkdiJaRtX+NuC7r6ZYL7C6kCf/PvDv6yPBr4iGJ+qqPqygNnhf6+NUF8Lgi+VLA5eZpc8Yt0meY67kWAspb3A5e4+JyYBV3Bx+Bnp6u4WaQF+HxabGG6rF7szqTji0CbZBFN5QITbfmaWQnArF4KBFCWPOPzO+jJcn2xmBf0M5IylVNCVRomhmLdxokdULYsFuJKDI6uelGd7E2A1EYa0J3iMfB3wRIT6zuTg9Dcd82yvxcERYjX9Tdm2yQ1hm+wBeiX6/MrbEuv2KOR17kEjtCekTQge8d9PMPJ39zzbqwIPhfV9DqQl+tyTcYllexB8gX8RHjMLqJNnXxUOziqRDZyY6HMvLwtRjNBO8BTzOoI+oaVu4wJfJ9FvRhm+6Y9wcDLG1wgmYPw+3PYyh89Jl/MlMLeA+vJO3PwmwSPMW8Nt76GJm8usTYBTOTgh8FqC+esiLX9O9Dkn8xLrn5ECXiPnGCVXCWgTgimIDhAkWe8RTOX1CQcncu6U6HNO5iWW7UHQaT1niptvCe6IvEQwDIOHbTQs0eeczAtwevg5zllyJlj+T97t+Y6ZHpaZHos2LmipDH2uAHD3oWb2LjCMYI6gqgTZ61RgkrsfKGZ9fzKzD4BfEVzJSiX4oXiY4Etcj/0XIYZtUo+Dj9GeHC6RbAR+XfKIK7ZY/4xI6cXh99ZfzexDgp+Dswm+nLYQPEE4zt03xDD8CieW7eHub5lZe+B2gvGuuob1fQU8C4x39/diewYVTh2ChwPyaxlhW1Ri1cYWZmoiIiIiEgMVvkO7iIiISFlSciUiIiISQ0quRERERGJIyZWIiIhIDCm5EhEREYkhJVciIiIiMaTkSkRERCSGlFyJVEJm1sPM3MxGJTqWZGFm1cL3pFhzUprZk+FxR8crNhEpX5RciSSp8Au7sOX6RMcYC2Z2b4Rz221mH5nZBDM7KkniOzeRcUTLzAZFeD/3mtlGM5tpZu1i9DpK0EUKUGmmvxEpx35fwPYVZRpF/GUA88J/NwYuJJiC4hdm1jHeU7O4+z4zO4VgQvbiuAO4l2DakmSyHHg1/Hdd4FygL3ClmZ2vqVVE4kfJlUiSc/d7Eh1DGXnb3e/N+Y+ZpRBMit4VGAncEO8A3H1dCY7ZQjA/X7JZlv+zY2aPAenAGIKJg0UkDnRbUKQCMLNWZna/mS0xs2/C20AbzOxvxbmtZmYnmtljZvaJmWWa2X/N7EMzm2Rm9SOUv9bM5prZdjPbY2ZrzOy3Zla9tOfk7tnAlPC/HfO9bvMwpo3huX5tZrPM7LQIMdYws1vNbLmZfWdmu8L35hUz65an3GF9rszsc4LEDmB+ntts+/KUOaTPlZl1Dv//fEHnZmYfh+9vvXzbLzaz2eH7vjdshz+ZWZ3o37lCPR6uz4wQU9SfITN7Engr/O+YfLcgz81XNm6fEZFkpStXIhVDH2Awwa21fwPZQDuCqz2XmFmH8ApLgcIv0MVALeD/gBeBNOAEoB8wHvguT/kZ4fZNYdnvgZ8AY4FuZnahu+8v5XlZuM6dYd7MTgTeBZoBc4CngWMJ3oOfmtkV7j47Tx0zw30fADOAPcBRQGfgAuDtQl7/IeDysOw0gnMFOFDQAe4+38w+AS4zs/ru/l3e/Wb2E+Ak4Dl3355n+x+A0cB/gdeAb4D2BLcdLzazn7j7D4XEGo2c9zM7wr7ifIZeIngPruPQ27lw8D0qq8+ISPJxdy1atCThQpBQOHBPhOX6fGWPBmpEqONigi/Bv+bb3iOse1SebbeF24ZFqKcWkJrn/4PCss/n3R7uG1NQPQWc5735Ywm3pwBzw31/y7P9X+G2u/KV7wzsJ0hKaobbGoTn/x5QJcJrN8zz72phvXMKiO/cAuJ/Mtx/dJ5to8NtN0Yo/7dw38V5tvUMt80H6uYrn/NePxDl+5lT/rEI+6aF+16OsK/Un6EC4ij1Z0SLlvK26MqVSPL7XYRt7wDTc/7j7p9HOtDdZ5vZOoLO4dHKjFDPznybRgBZwCB335Nv3++B4cC1wMRivG43M8v5ndQIuAg4kSBZGgdgZscD3YDPgAfzxTg/vBX3S4KrTU8TfIEbsNfdD7va5O7/LUZ8xfEEwfvQH5ics9HMUoFfEPTRejNP+VvC9SB3/z5fjI+Z2QiC9/OOYsRwupndE/67LkHyeQbweaR6YvwZgvh8RkTKBSVXIknO3a2oMmZmBLdo+gM/BuoDVfMU2R3FS/2d4IrCZDPrBbxBcHtorbvnvS1XG2gLbAVuD176MHuAU6J4zbzODxcIvpQ3AY8Af3T3L8LtOX2q5rn7Pg73NkFydRrwtLt/Z2azCW6rLSe4nTUfeN/dD0siY8XdN5rZXOB8M/uRu/8n3PUzoB4wxQ+9HdYJ2AtcXcD7WQ040szq5k++CnEaB9+vHBuAzpESqRh9hnLqitdnRKRcUHIlUjE8THAl4Evgn8AXBF9eAAOB5kVV4O6fmtlZBFfKLgR6h7s2mdkD7j4h/H+DcN2UyFfVckRKfgoz2vM8LViAuuG6oP5jOdvzdhS/ErgbuBr4Q7gt08xeAH7t7t8UM85oTSdIFvtzsFN8/3A9I1/ZBgRX2Ap7PyG4PRttcvW4uw8Kk6amBH2n/gC8ambnREguS/0ZyiNenxGRckHJlUg5Z2ZHEowHtZKgX9DOfPuvi7Yud19NMK5UNYLO1BcANwN/NbMf3H0GB7/cF7t7xwKqipec125WwP4j85XD3XcD/wP8j5kdC5wHDCDoaH0sB6+Wxdosglte15nZaKAJwfu5NHyf89oBZLl7k1gHEV51/Irgqb4GwK0Et+XuzCkTy89QKJGfEZGE01AMIuXfiQRXPd6I8KV4HHB8cSt0933uvtTdxxH0i4GgHxMePOH2EdAu/1ACZWB5uO5sZlUj7M9JlJZFOtjdN7n7kwRJzmdAVzOrG6lsHjm37yK9XoHcfRfBE3LHhHH1DevIf9UKgg73jc2sVXFeowR+R/A04i1hopmjJJ+hAt+XBH9GRBJOyZVI+bchXB+ScIT9Xh4lyp9zM+toZpGunDQN13n73DwEpAKPR0pOzKxBpDGnSsuDUdozCJKBm/O95jnAVQTJw9/DbU3N7LAxnYAjwiWbom9N5XR6P7bQUpFND9f9wiWboKN9fg+F68fCq0iHMLNa4S3bUnH3HcADQA0OvV23IVwX5zNU1PuSkM+ISDLQbUGRcs7dPzezFwn6Fi0LB8GsS3B1ZifwIdA6iqr6AYPN7B1gPbCdYDymSwn63ozP85qPmtkZBOMidTGzNwk6oDcAWhA8mTaFoA9PrA0hGOfqf83sYmApB8e52kcwTEXOFDbHAIvMbA3B1azPCd6bSwhu0z2Up2xB3iZ46vB+M2tP8L4ccPc/RhHrPIIrZFcTDC3xcqQnFN39TQvm6BsDfBx2wv+MoI/V8UAXgqTykihesygTgNuB/mZ2v7v/p4SfoTUEfdyuNbP9wGaC92mGu29O8GdEJLESPRaEFi1aIi+E41xFWfYIguEK1hMkQpsIvkTrEyQi+/KVjzTOVSeCYQM+ALYRDMmwHpgKtC7gdS8DXicYLiGboG/P+wRJQqsoY484zlURxxwdxrqJ4MnCbwmeBOyQr1x9gis0GQQdtPcSJAQZBFe5LE/ZiONchfv6E/RHygzL7Muz77BxrvIde09OWwI/K+K8ziO4lfhleF7fENwKfRA4Pcr3psBxrvKUyRnT7LmSfobCY84K38sdec7x3HxlSv0Z0aKlvC3mnvuEtYiIiIiUkvpciYiIiMSQkisRERGRGFJyJSIiIhJDSq5EREREYkjJlYiIiEgMKbkSERERiSElVyIiIiIxpORKREREJIaUXImIiIjEkJIrERERkRj6f+M0ewIWTMvPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fpr, tpr, threshold = roc_curve(np.array(ytrue, np.int), ypredproba)\n",
    "roc_auc = auc(fpr, tpr)\n",
    "\n",
    "plt.figure(figsize=(9, 6))\n",
    "plt.title('Pine Receiver Operating Characteristic')\n",
    "plt.plot(fpr, tpr, label = 'AUC = %0.2f' % roc_auc)\n",
    "plt.legend(loc = 'lower right')\n",
    "plt.plot([0, 1], [0, 1],'k--')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.savefig(\"results/logistic_regression/Logistic_Regression_ROC_Curve_Pine.png\", dpi=150)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import stats library\n",
    "from statsmodels.discrete.discrete_model import Logit\n",
    "import statsmodels.api as sm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.368989\n",
      "         Iterations 7\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>y</td>        <th>  No. Observations:  </th>  <td>   281</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>               <td>Logit</td>      <th>  Df Residuals:      </th>  <td>   277</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>MLE</td>       <th>  Df Model:          </th>  <td>     3</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Tue, 15 Jun 2021</td> <th>  Pseudo R-squ.:     </th>  <td>0.4651</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>09:19:24</td>     <th>  Log-Likelihood:    </th> <td> -103.69</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>           <td>True</td>       <th>  LL-Null:           </th> <td> -193.83</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th> </th>                      <td> </td>        <th>  LLR p-value:       </th> <td>7.627e-39</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   -0.4182</td> <td>    0.178</td> <td>   -2.350</td> <td> 0.019</td> <td>   -0.767</td> <td>   -0.069</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>   -1.4436</td> <td>    0.255</td> <td>   -5.667</td> <td> 0.000</td> <td>   -1.943</td> <td>   -0.944</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>    <td>   -2.5400</td> <td>    0.317</td> <td>   -8.010</td> <td> 0.000</td> <td>   -3.162</td> <td>   -1.918</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x3</th>    <td>   -0.1508</td> <td>    0.179</td> <td>   -0.842</td> <td> 0.400</td> <td>   -0.502</td> <td>    0.200</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   No. Observations:                  281\n",
       "Model:                          Logit   Df Residuals:                      277\n",
       "Method:                           MLE   Df Model:                            3\n",
       "Date:                Tue, 15 Jun 2021   Pseudo R-squ.:                  0.4651\n",
       "Time:                        09:19:24   Log-Likelihood:                -103.69\n",
       "converged:                       True   LL-Null:                       -193.83\n",
       "                                        LLR p-value:                 7.627e-39\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         -0.4182      0.178     -2.350      0.019      -0.767      -0.069\n",
       "x1            -1.4436      0.255     -5.667      0.000      -1.943      -0.944\n",
       "x2            -2.5400      0.317     -8.010      0.000      -3.162      -1.918\n",
       "x3            -0.1508      0.179     -0.842      0.400      -0.502       0.200\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# standardize X and add a constant for LR model \n",
    "X_std = StandardScaler().fit(X).transform(X)\n",
    "X_L = sm.add_constant(X_std)\n",
    "\n",
    "# fit LR model to data \n",
    "L = Logit(DOC_ctg, X_L)\n",
    "results = L.fit()\n",
    "\n",
    "coef = results.params\n",
    "std_err = np.diagonal(results.cov_params())**.5\n",
    "wald = (coef / std_err)**2\n",
    "pvalues = results.pvalues\n",
    "\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p value on out-of-fold predictions: 1.0347659251064246e-10\n"
     ]
    }
   ],
   "source": [
    "import scipy\n",
    "\n",
    "def Z_proportions(p1, p2, n1, n2):\n",
    "    P = (p1*n1 + p2*n2) / (n1 + n2)\n",
    "    Z = (p1 - p2) / np.sqrt(P*(1-P) * (1/n1 + 1/n2))\n",
    "    return Z\n",
    "\n",
    "baseline = np.max([len(ytrue)-np.sum(ytrue), np.sum(ytrue)]) / len(ytrue)\n",
    "acc = np.sum(ypredmax==ytrue) / len(ytrue)\n",
    "z_score = Z_proportions(acc, baseline, len(ypredmax), len(ytrue))\n",
    "model_p_value = scipy.stats.norm.sf(z_score)\n",
    "print(\"p value on out-of-fold predictions: {}\".format(model_p_value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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>Factor</th>\n",
       "      <th>Coefficient</th>\n",
       "      <th>Wald stat.</th>\n",
       "      <th>P values</th>\n",
       "      <th>Model P value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Total Biomass</td>\n",
       "      <td>-1.443560</td>\n",
       "      <td>32.117473</td>\n",
       "      <td>1.451262e-08</td>\n",
       "      <td>1.03477e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Bact Richness</td>\n",
       "      <td>-2.539996</td>\n",
       "      <td>64.153475</td>\n",
       "      <td>1.150948e-15</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Fungi Richness</td>\n",
       "      <td>-0.150824</td>\n",
       "      <td>0.708784</td>\n",
       "      <td>3.998477e-01</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Factor  Coefficient  Wald stat.      P values Model P value\n",
       "0   Total Biomass    -1.443560   32.117473  1.451262e-08   1.03477e-10\n",
       "1   Bact Richness    -2.539996   64.153475  1.150948e-15              \n",
       "2  Fungi Richness    -0.150824    0.708784  3.998477e-01              "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fs_df = pd.DataFrame()\n",
    "fs_df['Factor'] = feature_names\n",
    "fs_df['Coefficient'] = coef[1:]\n",
    "fs_df['Wald stat.'] = wald[1:]\n",
    "fs_df['P values'] = pvalues[1:]\n",
    "fs_df[\"Model P value\"] = [model_p_value] + [\"\"]*(len(feature_names)-1)\n",
    "\n",
    "fs_df.to_csv('results/logistic_regression/FeatureImportances_pine.csv', index=False)\n",
    "fs_df.head(10)"
   ]
  }
 ],
 "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}