Final_Project_Code_py.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final Project Code: Linear Models in Python\n",
    "# Jacob Tiede"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "from sklearn.preprocessing import Normalizer\n",
    "from sklearn.model_selection import KFold\n",
    "\n",
    "#-----------\n",
    "#Read in the data and preprocess\n",
    "#-----------\n",
    "data = pd.read_csv('./Data/train_features.csv')\n",
    "data = data.drop(['sig_id'], axis=1)\n",
    "\n",
    "#Some columns are coded with strings, need them to be numeric\n",
    "label_enc = LabelEncoder()\n",
    "data['cp_type'] = label_enc.fit_transform(data.cp_type.values)\n",
    "data['cp_time'] = label_enc.fit_transform(data.cp_time.values)\n",
    "data['cp_dose'] = label_enc.fit_transform(data.cp_dose.values)\n",
    "\n",
    "#convert to numpy:\n",
    "col_names_data = list(data.columns) #keep a list of column names for later\n",
    "data = data.to_numpy()\n",
    "\n",
    "#Column-wise normalization usually helps with these problems:\n",
    "transformer = Normalizer().fit(data)\n",
    "data = transformer.transform(data)\n",
    "#add the bias column: NOT REQUIRED, but it helps with some of the logic of the program (this could be done in pytorch automatically, by allowing the linear layer to have a bias)\n",
    "data = np.append(np.ones(data.shape[0]).reshape(-1,1), data, 1)\n",
    "\n",
    "#load in the response:\n",
    "response = pd.read_csv('./Data/train_targets_scored.csv')\n",
    "response = response.drop(['sig_id'], axis=1) #Drop ID since it cannot have a meaningful impact on explanation\n",
    "col_names_response = list(response.columns) #save a list of response column names just in case\n",
    "response = response.to_numpy()\n",
    "\n",
    "#Find k folds for k-fold cross validation:\n",
    "kf = KFold(n_splits=10, random_state = 0, shuffle = True)\n",
    "kf.get_n_splits(data)\n",
    "kcross = kf.split(data, response) #creates an iteratable that we will use later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<torch._C.Generator at 0x2144b068cf0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#create the neural net class NOTE: this assumes you have a working GPU and have installed pytorch with gpu libraries\n",
    "\"\"\"\n",
    "log_reg -- a class that stores the parameters of a multivariate multiple logistic regress\n",
    "Note: This class uses a sigmoid link function instead of a logit link function\n",
    "\n",
    "@params: --input_features: int, the number of features which the model will use as predictors\n",
    "         --output_features: int, the number of expected outputs, this is the number of logistic regression models to be trained simultaniously\n",
    "\"\"\"\n",
    "class log_reg(nn.Module):\n",
    "    def __init__(self, input_features, output_features):\n",
    "        super(log_reg, self).__init__()\n",
    "        self.fc_layer = nn.Linear(input_features, output_features, bias = False)\n",
    "    def forward(self, x):\n",
    "        x = F.sigmoid(self.fc_layer(x))\n",
    "        return x\n",
    "torch.manual_seed(35) #set seed for reproducability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\JAT\\anaconda3\\envs\\pytorch\\lib\\site-packages\\torch\\nn\\functional.py:1625: UserWarning: nn.functional.sigmoid is deprecated. Use torch.sigmoid instead.\n",
      "  warnings.warn(\"nn.functional.sigmoid is deprecated. Use torch.sigmoid instead.\")\n",
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "  1% (103 of 10020) |                    | Elapsed Time: 0:00:00 ETA:   0:00:09"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 1 were:\n",
      "accuracy: 0.9971635915063078\n",
      "loss: 0.015549170781863842\n",
      "The train statistics for fold 1 were:\n",
      "accuracy: 0.9970685187614802\n",
      "loss: 0.01800676988991531\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 2 were:\n",
      "accuracy: 0.9971545077321087\n",
      "loss: 0.015628117741738074\n",
      "The train statistics for fold 2 were:\n",
      "accuracy: 0.9971956179480451\n",
      "loss: 0.017190761436280365\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "  1% (109 of 10020) |                    | Elapsed Time: 0:00:00 ETA:   0:00:09"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 3 were:\n",
      "accuracy: 0.9971685875821173\n",
      "loss: 0.015561347609820838\n",
      "The train statistics for fold 3 were:\n",
      "accuracy: 0.9970644187877198\n",
      "loss: 0.017855714629026683\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "  1% (104 of 10020) |                    | Elapsed Time: 0:00:00 ETA:   0:00:09"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 4 were:\n",
      "accuracy: 0.9971617747514677\n",
      "loss: 0.01554152829583712\n",
      "The train statistics for fold 4 were:\n",
      "accuracy: 0.9971320683547625\n",
      "loss: 0.018008565681206214\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 5 were:\n",
      "accuracy: 0.9971590496192082\n",
      "loss: 0.015583967087333074\n",
      "The train statistics for fold 5 were:\n",
      "accuracy: 0.9971689681186041\n",
      "loss: 0.017476386541651713\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 6 were:\n",
      "accuracy: 0.9971556432038834\n",
      "loss: 0.015624830691690395\n",
      "The train statistics for fold 6 were:\n",
      "accuracy: 0.9972079178693257\n",
      "loss: 0.017235161107335542\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 7 were:\n",
      "accuracy: 0.9971592767135632\n",
      "loss: 0.015570653716804918\n",
      "The train statistics for fold 7 were:\n",
      "accuracy: 0.9971751180792442\n",
      "loss: 0.01787060801241849\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 8 were:\n",
      "accuracy: 0.997171766903087\n",
      "loss: 0.01553338622774728\n",
      "The train statistics for fold 8 were:\n",
      "accuracy: 0.9970377689582787\n",
      "loss: 0.01818494118649412\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 9 were:\n",
      "accuracy: 0.997160412185338\n",
      "loss: 0.015548926223247887\n",
      "The train statistics for fold 9 were:\n",
      "accuracy: 0.9970808186827606\n",
      "loss: 0.01803042455198797\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 10 were:\n",
      "accuracy: 0.9971658624498577\n",
      "loss: 0.015575567316152379\n",
      "The train statistics for fold 10 were:\n",
      "accuracy: 0.9970767187090002\n",
      "loss: 0.017709220533032675\n"
     ]
    }
   ],
   "source": [
    "import progressbar #This is the progressbar2 library if it must be installed\n",
    "\n",
    "fold = 1 #keeps track of which fold we are on\n",
    "\n",
    "#do training loop for each fold in the K-fold cross:\n",
    "for train_index, test_index in kcross:\n",
    "    X_train, X_test = data[train_index], data[test_index]\n",
    "    y_train, y_test = response[train_index], response[test_index]\n",
    "    \n",
    "    #Initialize necessary constants like batch size, etc as well as the model:\n",
    "    model = log_reg(data.shape[1], response.shape[1]).cuda() #again requires a working gpu\n",
    "    opt = torch.optim.Adam(model.parameters()) #Initializes and tells the optimizer what parameters to keep track of\n",
    "    batch_size = 64 #This could probably be increased, but 64 seems to work well enough\n",
    "    epochs = 30\n",
    "    loss = nn.BCELoss() #Binary cross entropy loss, the loss function defined in the write up\n",
    "    batchesPerEpoch = int(np.floor(X_train.shape[0]/batch_size))\n",
    "    \n",
    "    #-----------\n",
    "    #training loop\n",
    "    #-----------\n",
    "    with progressbar.ProgressBar(max_value=epochs*batchesPerEpoch) as bar:\n",
    "        for epoch in range(epochs):\n",
    "            for i in range(batchesPerEpoch):\n",
    "                opt.zero_grad()\n",
    "                batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "                batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "                #forward pass:\n",
    "                out = model(batch_x)\n",
    "                batch_loss = torch.mean(loss(out, batch_y.float()))\n",
    "                batch_loss.backward()\n",
    "                opt.step()\n",
    "                bar.update(i+batchesPerEpoch*epoch)\n",
    "    #-----------            \n",
    "    #Get train stats:\n",
    "    #-----------\n",
    "    accuracy = []\n",
    "    test_loss = []\n",
    "    model.eval()\n",
    "    with torch.no_grad():\n",
    "        for i in range(batchesPerEpoch):\n",
    "            batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            out = model(batch_x)\n",
    "            test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "            accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "    print(\"The train statistics for fold \" + str(fold) + \" were:\")\n",
    "    print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "    print(\"loss: \" + str(np.mean(test_loss)))    \n",
    "    \n",
    "    #-----------\n",
    "    #Get test stats:\n",
    "    #-----------\n",
    "    accuracy = []\n",
    "    test_loss = []\n",
    "    batches = int(np.floor(X_test.shape[0]/batch_size))\n",
    "    with torch.no_grad():\n",
    "        for i in range(batches):\n",
    "            batch_x = torch.from_numpy(X_test[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_test[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            out = model(batch_x)\n",
    "            test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "            accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "    print(\"The train statistics for fold \" + str(fold) + \" were:\")\n",
    "    print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "    print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "    fold += 1\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (50 of 50) |########################| Elapsed Time: 0:09:06 Time:  0:09:06\n"
     ]
    }
   ],
   "source": [
    "# To get p-values of the coefficients we need standard errors, to do this we will use bootstrapping:\n",
    "n = 50 #number of bootstrapped samples to consider\n",
    "m = int(data.shape[0]*99/100) #number of bootstrapped samples\n",
    "boot_coeff = torch.zeros(206,876,n) #3D array of coefficients dimension each slice along dim 2 is a mxn matrix of parameters\n",
    "\n",
    "#Also going to use this to test my chi-squared test for model sufficiency, for which we'll need to record some stuff:\n",
    "boot_pred = np.zeros((m,response.shape[1], n))\n",
    "boot_response = np.zeros((m,response.shape[1], n))\n",
    "\"\"\"\n",
    "generateBootSample -- generates a sample of rows from X and y with replacement\n",
    "\n",
    "@params -- X a numpy array to sample \n",
    "        -- y a second numpy array to generate the same sample from\n",
    "        -- size an int for the size of the bootstrapped sample\n",
    "\"\"\"\n",
    "def generateBootSample(X,y, size):\n",
    "    samples = torch.randint(0, int(X.shape[0]), (size, 1)).view(-1).numpy()\n",
    "    \n",
    "    return X[samples, :], y[samples,:]\n",
    "\n",
    "#-----------\n",
    "#train a model for each bootstrapped sample required, store the parameters in boot_coeff for each one\n",
    "#-----------\n",
    "with progressbar.ProgressBar(max_value=n) as bar:\n",
    "    for sample in range(n):\n",
    "        X_train, y_train = generateBootSample(data, response, m)\n",
    "        #Initialize necessary constants like batch size, etc as well as the model:\n",
    "        model = log_reg(data.shape[1], response.shape[1]).cuda()\n",
    "        opt = torch.optim.Adam(model.parameters())\n",
    "        batch_size = 64\n",
    "        epochs = 30\n",
    "        loss = nn.BCELoss()\n",
    "        batchesPerEpoch = int(np.floor(X_train.shape[0]/batch_size))\n",
    "\n",
    "        for epoch in range(epochs):\n",
    "            for i in range(batchesPerEpoch):\n",
    "                opt.zero_grad()\n",
    "                batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "                batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "                #forward pass:\n",
    "                out = model(batch_x)\n",
    "                batch_loss = torch.mean(loss(out, batch_y.float()))\n",
    "                batch_loss.backward()\n",
    "                opt.step()\n",
    "        \n",
    "        #-----------\n",
    "        #Get stuff required to evaluate the chi-squared test later\n",
    "        #-----------\n",
    "        boot_response[:,:,sample] = y_train\n",
    "        model.eval()\n",
    "        with torch.no_grad():\n",
    "            for i in range(batchesPerEpoch):\n",
    "                batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "                batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "                out = model(batch_x)\n",
    "                boot_pred[i*batch_size:(i+1)*batch_size,:, sample] = out.detach().cpu().numpy()\n",
    "                \n",
    "        bar.update(sample)\n",
    "        for param in model.parameters():\n",
    "            boot_coeff[:,:,sample] = param.detach().cpu()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:11 Time:  0:00:11\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for the full model:\n",
      "accuracy: 0.9971699501482473\n",
      "loss: 0.01553173159377393\n",
      "The train statistics for the full model were:\n",
      "accuracy: 0.9970726187352399\n",
      "loss: 0.017991432441851578\n"
     ]
    }
   ],
   "source": [
    "standard_errors = torch.std(boot_coeff, dim = 2) #estimate the standard error of the coefficients using the bootstrapped sample\n",
    "\n",
    "#-----------\n",
    "#train a final model:\n",
    "#-----------\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(data, response, test_size=0.1, random_state=0) #create a final split\n",
    "\n",
    "#Initialize necessary constants like batch size, etc as well as the model:\n",
    "model = log_reg(data.shape[1], response.shape[1]).cuda()\n",
    "opt = torch.optim.Adam(model.parameters())\n",
    "batch_size = 64\n",
    "epochs = 30\n",
    "loss = nn.BCELoss()\n",
    "batchesPerEpoch = int(np.floor(X_train.shape[0]/batch_size))\n",
    "\n",
    "#-----------\n",
    "#training loop\n",
    "#-----------\n",
    "with progressbar.ProgressBar(max_value=epochs*batchesPerEpoch) as bar:\n",
    "    for epoch in range(epochs):\n",
    "        for i in range(batchesPerEpoch):\n",
    "            opt.zero_grad()\n",
    "            batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            #forward pass:\n",
    "            out = model(batch_x)\n",
    "            batch_loss = torch.mean(loss(out, batch_y.float()))\n",
    "            batch_loss.backward()\n",
    "            opt.step()\n",
    "            bar.update(i+batchesPerEpoch*epoch)\n",
    "\n",
    "#-----------\n",
    "#Get train stats:\n",
    "#-----------\n",
    "accuracy = []\n",
    "test_loss = []\n",
    "model.eval()\n",
    "pred_probs = np.zeros(y_train.shape)\n",
    "with torch.no_grad():\n",
    "    for i in range(batchesPerEpoch):\n",
    "        batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "        batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "        out = model(batch_x)\n",
    "        pred_probs[i*batch_size:(i+1)*batch_size,:] = out.cpu().numpy()\n",
    "        test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "        accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "print(\"The train statistics for the full model:\")\n",
    "print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "\n",
    "#-----------\n",
    "#Get test stats:\n",
    "#-----------\n",
    "accuracy = []\n",
    "test_loss = []\n",
    "batches = int(np.floor(X_test.shape[0]/batch_size))\n",
    "pred_probs_test = np.zeros(y_test.shape)\n",
    "with torch.no_grad():\n",
    "    for i in range(batches):\n",
    "        batch_x = torch.from_numpy(X_test[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "        batch_y = torch.from_numpy(y_test[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "        out = model(batch_x)\n",
    "        pred_probs_test[i*batch_size:(i+1)*batch_size,:] = out.cpu().numpy()\n",
    "        test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "        accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "print(\"The train statistics for the full model were:\")\n",
    "print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "\n",
    "#-----------\n",
    "#get p-values for coefficients:\n",
    "#-----------\n",
    "#get z-scores of the paramters\n",
    "for param in model.parameters():\n",
    "    z_scores = torch.abs(torch.div(param.detach().cpu(), standard_errors)).numpy()\n",
    "\n",
    "#get the p-values\n",
    "from scipy.stats import norm\n",
    "p_vals = 1 - norm.cdf(z_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The lowest p-value:\n",
      "0.0\n",
      "The highest p-value:\n",
      "0.499997732682243\n",
      "The proportion of parameters that are significant (across all models):\n",
      "0.6904009841734273\n",
      "The number of predictors that are non-significant in all models:\n",
      "0\n",
      "The number of predictors that are non-significant in all but 10 or less models:\n",
      "1\n",
      "The number of predictors that are non-significant in all but 50 or less models:\n",
      "108\n"
     ]
    }
   ],
   "source": [
    "print(\"The lowest p-value:\")\n",
    "print(np.min(p_vals))\n",
    "print(\"The highest p-value:\")\n",
    "print(np.max(p_vals))\n",
    "print(\"The proportion of parameters that are significant (across all models):\")\n",
    "print(np.sum(p_vals < .05)/(p_vals.shape[0]*p_vals.shape[1]))\n",
    "\n",
    "count = 0\n",
    "for i in range(p_vals.shape[1]):\n",
    "    if np.sum(p_vals[:,i]<.05) == 0:\n",
    "        count += 1\n",
    "print(\"The number of predictors that are non-significant in all models:\")\n",
    "print(count)\n",
    "count = 0\n",
    "for i in range(p_vals.shape[1]):\n",
    "    if np.sum(p_vals[:,i]<.05) <= 10:\n",
    "        count += 1\n",
    "print(\"The number of predictors that are non-significant in all but 10 or less models:\")\n",
    "print(count)\n",
    "\n",
    "count = 0\n",
    "for i in range(p_vals.shape[1]):\n",
    "    if np.sum(p_vals[:,i]<.05) <= 50:\n",
    "        count += 1\n",
    "print(\"The number of predictors that are non-significant in all but 50 or less models:\")\n",
    "print(count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of insufficient models: 0\n"
     ]
    }
   ],
   "source": [
    "#Lastly let's code a chi squared test for this:\n",
    "#Inspiration for some of this code: \n",
    "#https://stackoverflow.com/questions/50975774/calculate-residual-deviance-from-scikit-learn-logistic-regression-model\n",
    "from sklearn.metrics import log_loss\n",
    "\"\"\"\n",
    "resid_dev -- given predicted probabilities and actual realizations return the residual deviance\n",
    "\n",
    "@params: --pred_probs - a numpy array of predicted probabiliities\n",
    "         --y - a numpy array of actual classifications\n",
    "\"\"\"\n",
    "def resid_dev(pred_probs, y):\n",
    "    out = np.zeros(y.shape[1])\n",
    "    for i in range(y.shape[1]):\n",
    "        out[i] = 2*log_loss(y[:,i], pred_probs[:,i], normalize=False, labels = [0,1])\n",
    "    return out\n",
    "\n",
    "deviances = resid_dev(pred_probs, y_train) #find the residual deviances of the model\n",
    "\n",
    "#this is the residual deviance so we can construct chi squared tests for all of them (to see if the current model is sufficient)\n",
    "dof = response.shape[0] - data.shape[1] #degrees of freedom\n",
    "\n",
    "#-----------\n",
    "#Do the chi-squared test for model sufficiency\n",
    "#-----------\n",
    "from scipy.stats import chi2\n",
    "test_pvals = 1-chi2.cdf(deviances, dof)\n",
    "print(\"Number of insufficient models: \" + str(np.sum(test_pvals<.05)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEICAYAAAC0+DhzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABCWklEQVR4nO3deXxU9bn48c8zk0DYwxIQCJggAQ2y7wqIC7K4oLZVca9ai5Xba229anuvYm/ttbfX6++6olartip69bZSi4obAiKSIIgCIgEChDVkYQtkmXl+f5wTOkwmySQkOTPJ83695pWZc77fc55z5uQ8c77fs4iqYowxxoTyeR2AMcaY2GPJwRhjTBWWHIwxxlRhycEYY0wVlhyMMcZUYcnBGGNMFZYcGpCIvCgiv3HfTxSRjU00XxWR/vWse5OILKth/LsicmP9o6tTLMeXQ0TaiMjfROSAiPxvI883zZ13wklMY66I/Lkh42pstX33XhGRX4rIH5pwfteKyKKmml+8aHHJQURyReSoiBwWkb0i8kcRad/Q81HVpao6MIp4Gv0fVESmisgSETkkIvki8qmIXBpNXVWdrqovudO5SESWiUixiOwRkedEpEMjhf19oAfQVVV/0EjziJr7PX0tIiXusj8tIslex9VYQhLm4ZD/lXdEZEpjz1tVf6uqt4bFUa/EXd3/l7sfuMCd3yuqemEU0zr+468laHHJwXWJqrYHRgCjgX8NL3AyvyJjiYh8H/hf4GUgFWeHez9wST0m1wn4DdALOMOd3u8bJtIqTgW+U9WKRpp+1ETk58DvgLtx1sE4nPg+EJFWTRiHF9tksvu/MhT4APiLiNzkQRzNWkzub1S1Rb2AXOCCkM+/B95x3ytwB7AJ2OoOuxhYAxQDy4EhIXWHA18Ch4DXgfnAb9xxk4G8kLJ9gP8D8oEC4AmcHewxIAAcBordsq2B/wK2A3uBeUCbkGndDewGdgE3u3H3j7Cs4k7j7hrWx03AMnd+RcBWYHrI+MXArdXUvQL4upZ1fR+w3p32H4Gk2pYDeBAoA8rd9XJLhGmPAT53v5fd7vpsFTJegdnud1kEPAmIO87vLu9+YIv7nSuQEGE+Hd0Yrgwb3h7YB9zsfp4LvOluB4fc7WJoSPl7gJ3uuI3A+e5wH3AvsNndLt4Aurjj0ty4bnG/xyXAe8CcsFi+Aq5w35+OsxMvdOdzZUi5rsAC4CCwEvh3YFk1313lvBPChv8CZ5v0uZ97AW/hbNdbgZ+GlJ3rLs/L7nKvA0ZFsU7mAn9232934zjsvs5xl21wyHS6A0eBlOq275r2A6FlcP5nHnW/2wPAWuBM4Dac7bHMjeNvbvkzcP5Hit3luzRsff/NXd9ZOD+sloWMj7S/+R9gh1tnFTAxbH3+L/Bnd519DQzA+R/b59a7sMH2lQ01oXh5hW0Ufdwv9N9DvqwPgC5AG5wji33AWJwdyo1u/dZAK2Ab8DMgEacZpJwIycGt+5W70bUDkoAJ1W28wP/D+SfuAnRwN7D/cMdNw/nnPNOd1qtUnxxOd8el17A+bnLj/pEb5+04O+vKHeliqk8O/w+YX8u6/sZdz12Az0LWT43LQcgOopppj8T5BZ+AsyPbANwZ9o/3DpAM9MXZeU1zx80Gvg2J6xOqTw7TgIpqxr0EvBYSb7m7HSTi7ES3uu8H4vzj9nLLpgGnue/vBFbgHIW1Bp4JmWaaG9fL7jpqA9wAfBYSQybOjqm1W2YH8EN3vYzASYCD3LLzcXbW7dz1vpO6J4d+7vAzcBLbKpwj0VbuuC3A1JB1cgyYgbNt/Qewwh1X0zo5/t1HigN4CvhdyOd/xt1ZV7N91yU5THWXKRknUZwB9HTHvYi7/bqfE4Ec4Jfu8p+Hs9MeGLK+5wNt3e9pB1WTw/H9jTvsOpykkgD8HNiD+4MqZH1Odce/jLON/cqN5Ue4SaZB9pUNNaF4ebkbxWGcf6ht7oZW+cUocF5I2adxE0fIsI04v14mEbITdcctJ3JyGI+zc4q0gzlh43U3yCOV/ygh9be6718AHg4ZN4Dqk8PZ7rikGtbHTUBOyOe2bp1T3M+LiZAcgCk4v8gH1LKuZ4d8ngFsjmY5qCU5RJjXncBfQj4rbgJ2P78B3Ou+/zgsrgupPjlcB+ypZp4PAx+ExLsiZJwP54hmIs7R0D7gAiAxbBobcH8xu5974iSZyqSnQL+Q8R3c7eNU9/NDwAvu+6uApWHTfwZ4AGfnXA6cHjLut9Q9OSS5w8/G+dG0PWz8fcAfQ9bJhyHjMoGj7vua1snx7z5SHO58d/CPo5dswo7swrbvCpz/99BXkMjJ4TzgO5wfHr6wab3IiclhIs7O2xcy7DU3/sr1PTBkXKQjh/MixR1Spgj3CNSd7gch4y7B2Zf5Q7YNxWkKPOl9ZUvtc7hMVZNV9VRV/YmqHg0ZtyPk/anAz90O2GIRKcb5tdnLfe1U91txbatmfn2AbRpd+3kKzg56Vcg833OH4843NMbq5glOMwU4O5ya7Kl8o6ol7ttqO+lFZBzOL/3vq+p3tUw7PNZe7vu6LEekGAa4HaR7ROQgzo6uW1ixPSHvS/jHMtVl3vuBbtW0Cfd0x1c6Pk1VDQJ5OL+Mc3CS11xgn4jMF5HK9XAqTjt+5Xe9AaeZsUc10z0E/B242h10NfBKyLTGhm2v1wKn4Gw/CXVY7ur0dv8WuvPrFTa/X4bFHv4dJIlIQi3rpEaq+gVOgjxHRE7HSTQLaqiywv1/P/7Caa6KNO2PcZoonwT2isizItKxmun2Ana433WlbTjrKNL6Dn0fcZiI/FxENrhn6RXj9HGFbtd7Q94fBfaraiDkM9Twv1sXLTU51CR0Z78DeChsw2qrqq/h/CrsLSISUr5vNdPcAfStZgejYZ/343zJg0Lm2UmdTkHc+faJYp7gHOXsAL5XQ5k6EZHhOP+IN6vqR1FUCY91l/u+LssRydM4TUMZqtoRZ6ckNVc5ri7z/hwoxelfOU5E2gHTgdB10CdkvA+nqWgXgKq+qqoTcHaoitPBDc73Mz1sG0tS1Z0h0w3fRl4DZonIeJympk9CpvVp2LTaq+rtOEeuFXVY7upcjvOLv3Lb2ho2vw6qOiOaCdWwTk4oVk31l3CO6q4H3lTVY3VdkBriekxVRwKDcI5o764mll1AH/e7rtQXp7mucn2nhowLXffHZ1f5RkQm4vTDXAl0dpPYAaLfrhuUJYeaPQfMFpGx4mjnns7ZAWenUQH8VEQSROQKnE7SSFbi7JAedqeRJCJnu+P2AqmVZ724v0KeAx4Vke4AItJbRKa65d8AbhKRTBFpi9NkEJF7VHMX8G8i8kMR6SgiPhGZICLP1nVliMiZOEcx/6Sqf4uy2h0ikioiXXB24K/XdTmq0QGn0+6w++vx9jrUfQPne0sVkc44HcIRqeoBnA7yx0VkmogkikgaTsdgHvCnkOIjReQK90fAnThJZYWIDBSR80SkNU6b8VGcowNwTjZ4SEROBRCRFBGZWUv8C3F2qL8GXg/55foOMEBErnfjTBSR0SJyhvvr8v+AuSLSVkQycfrQoiIiPURkDs73dJ87z5XAQRG5R5zrUvwicqaIjI5iejWtk1D5OE1A/cKG/wknUV2H0/beINz1NVZEEnGOTipPGAHnfzU0jsojmH9x1/VknKae+RHW9+k4/UU16YCzT8kHEkTkfpwTIjxhyaEGqpqN08nzBE7bXw5O+ySqWobza/Imd9xVOBtDpOkEcDaa/jiHs3lueXDav9cBe0SksoniHndeK9wmkw9xOvBQ1XdxOoI/dst8XMsyvOnO62acXzp7cdo+345uLZzg5ziHy8/LP85/X1dLnVeBRTgdlVvcedd5OSL4BXANTgfgc/wj6UTjOeB9nJMEvqSa762Sqv4nTmL7L5yE9AXOr+bzVbU0pOjbOOu6COcX7RWqWo7TWfwwzlHhHpyza37p1vkfnCOxRSJyCKdzemwt8ZS6MV+As34rhx/C6T+5Gue73oPza7y1W2QOTpPDHpz28z/WNB9XsYgcwTkzZgbwA1V9wZ1f5XY9DKdjdD/wB5ymkNrUtE5Cl7UEp1/lM7fpapw7PA/nu1NgaRTzi1ZHnO2jCKeJqADnewd4Hsh04/iruw+4FOcIcj9O/+UNqvqtW34OzrrYg5PMXsP5wVCd94F3cfo8tuEkpkhNUU2i8owUYxqciOTidGZ/6HUspvkRkReAXapa5TqlWCQiv8M50SPqIzYvxd6FF8YYUwu3ae8KnGuNYpLblNQK56hrNM71Krd6GlQdWLOSMSauiMi/41w/83tV3ep1PDXogNP8dwSnn+sR6tec6wlrVjLGGFOFHTkYY4ypoln0OXTr1k3T0tK8DsMYY+LKqlWr9qtqSqRxzSI5pKWlkZ2d7XUYxhgTV0Sk2qvkrVnJGGNMFZYcjDHGVGHJwRhjTBXNos/BmOagvLycvLw8jh1rsHvIGQNAUlISqampJCYmRl3HkoMxMSIvL48OHTqQlpbGiTf7Nab+VJWCggLy8vJIT0+Pup41KxkTI44dO0bXrl0tMZgGJSJ07dq1zkeklhyMiSGWGExjqM92ZcmhBVNVKgLB2gsaY1ocSw4tTDCoLPhqF1c+8zkD//U9+v/qXSb//hP+6/2NFByu6VbzpiXw+/0MGzaMM888kx/84AeUlJTUXqkaN910E2+++SYAt956K+vXr6+27OLFi1m+fHmd55GWlsb+/ftrLxhBbm4uZ5555vHPs2bNYsiQITz66KN1nlZ5eTkjR46sMvxXv/oVffr0oX37E5/cuX37ds4991yGDx/OkCFDWLhw4fFxL730EhkZGWRkZPDSSy816XKEsg7pFmTvwWP87PU1LN9cQHq3dtww/lTatk7gqx3FPLU4h1e+2MZ/XDGEaWee4nWoxiNt2rRhzZo1AFx77bXMmzePu+666/j4QCCA3++v83T/8Ic/1Dh+8eLFtG/fnrPOOqvO024Ie/bsYfny5WzbVp/HasOyZcsixn7JJZcwZ84cMjIyThj+m9/8hiuvvJLbb7+d9evXM2PGDHJzcyksLOTBBx8kOzsbEWHkyJFceumldO7cuUmWI5QdObQQOwpL+P685azZUczDVwzmo7vO4V8vzuSuKQN46eYxvHfnJPp2acvtr6zixc9i+S7IpqlMnDiRnJwcFi9ezLnnnss111zD4MGDCQQC3H333YwePZohQ4bwzDPPAE4z5Zw5c8jMzOSiiy5i3759x6c1efLk47e4ee+99xgxYgRDhw7l/PPPJzc3l3nz5vHoo48ybNgwli5dSn5+Pt/73vcYPXo0o0eP5rPPPgOgoKCACy+8kOHDh/PjH/+Y6u4q3b59e37+858zYsQIzj//fPLz8wFYtWoVQ4cOZfz48Tz55JPHy1944YXs27fv+PwrBQIB+vXrh6pSXFyMz+djyZIlJ6yfymWaPn16lTjGjRtHz549qwwXEQ4ePAjAgQMH6NWrFwDvv/8+U6ZMoUuXLnTu3JkpU6bw3nvvValf1+WoDztyaAGKjpRx3fNfcPBoBa/9aBxD+yRXKTOgRwde//F4fvraaub+bT0dkhL53sjUqhMzTeLBv61j/a6DDTrNzF4deeCSQVGVraio4N1332XatGkArFy5km+++Yb09HSeffZZOnXqRFZWFqWlpZx99tlceOGFrF69mo0bN/L111+zd+9eMjMzufnmm0+Ybn5+Pj/60Y9YsmQJ6enpFBYW0qVLF2bPnk379u35xS9+AcA111zDz372MyZMmMD27duZOnUqGzZs4MEHH2TChAncf//9/P3vf+fZZyM/Cv3IkSOMGDGCRx55hF//+tc8+OCDPPHEE/zwhz/k8ccf55xzzuHuu+8+Xn7BggVcfPHFx4+aKvn9fgYMGMD69evZunUrI0eOZOnSpYwdO5a8vDz69+8PwCeffMIDD0T/GPS5c+dy4YUX8vjjj3PkyBE+/NB5WOLOnTvp06fP8XKpqans3LmzSv26Lkd92JFDMxcMKv/02mp2Fx/jhZtGR0wMlZIS/TxxzQjO7t+Ve95ay5odxU0Wp4kNR48eZdiwYYwaNYq+fftyyy23ADBmzJjj58gvWrSIl19+mWHDhjF27FgKCgrYtGkTS5YsYdasWfj9fnr16sV5551XZforVqxg0qRJx6fVpUuXiHF8+OGHzJkzh2HDhnHppZdy8OBBDh06xJIlS7juuusAuOiii6ptbvH5fFx1lfOY9uuuu45ly5Zx4MABiouLOeeccwC4/vrro1onEydOZMmSJSxZsoT77ruPZcuWkZWVxejRowHYtWsXXbp0oW3btlFND+C1117jpptuIi8vj4ULF3L99dcTDAYjHgmFn2lU3+WoKztyaOZe/jyXZTn7efiKwYw8tfZ2y1YJPp66ZiQzHlvKP89fzTv/NIEOSdFfVWkaRrS/8BtaaJ9DqHbt2h1/r6o8/vjjTJ069YQyCxcurPWUSVWN6rTKYDDI559/Tps2baqMq89pmSIS9bzDTZw4kXnz5rFr1y5+/etf8/vf/57FixczadIkAN59990q66I2zz///PHmovHjx3Ps2DH2799PamoqixcvPl4uLy+PyZMnn1C3vstRV3bk0Izl7j/Cw+99y3mnd+eq0X1qr+Dq1DaR/7l6GDsKS3jo7xsaMUITj6ZOncrTTz9NeXk5AN999x1Hjhxh0qRJzJ8/n0AgwO7du/nkk0+q1B0/fjyffvopW7c6/VqFhYUAdOjQgUOHDh0vd+GFF/LEE08c/1yZsCZNmsQrr7wCODvloqKiiDEGg8HjZ0q9+uqrTJgwgeTkZDp16sSyZcsAjk+nNmPHjmX58uX4fD6SkpIYNmwYzzzzDBMnTgSq72+oSd++ffnoo48A2LBhA8eOHSMlJYWpU6eyaNEiioqKKCoqYtGiRVUST32Xo64sOTRjv/n7ehJ9Pn57+eA6/9IYldaFWyf2Y37WDr7cHvkf0LRMt956K5mZmYwYMYIzzzyTH//4x1RUVHD55ZeTkZHB4MGDuf322483e4RKSUnh2Wef5YorrmDo0KHHm34uueQS/vKXvxzvSH3sscfIzs5myJAhZGZmMm/ePAAeeOABlixZwogRI1i0aBF9+/aNGGO7du1Yt24dI0eO5OOPP+b+++8H4I9//CN33HEH48ePj3hUEknr1q3p06cP48aNA5wjiUOHDh3vnN+0aROnn356xLr/8i//QmpqKiUlJaSmpjJ37lwAHnnkEZ577jmGDh3KrFmzePHFFxERunTpwr/9278d74i///77Iza91Wc56qpZPEN61KhRag/7OdHyzfu55rkvuGfa6dw++bR6TeNwaQXnP7KYlA6tefuOCfh9dvVuY9qwYQNnnHGG12E0C+3bt+fw4cONPp9ly5bx5z//+XjyimWRti8RWaWqoyKVj+rIQUSmichGEckRkXsjjBcRecwdv1ZERrjD+4jIJyKyQUTWicg/h9TpIiIfiMgm92/nkHH3udPaKCJ1a8wzqCr/sfBbeie34Ydnp9V7Ou1bJ/DLGWfwzc6DLPiq6hkTxrR0EyZMiIvEUB+1JgcR8QNPAtOBTGCWiGSGFZsOZLiv24Cn3eEVwM9V9QxgHHBHSN17gY9UNQP4yP2MO/5qYBAwDXjKjcFEafHGfL7eeYA7L8ggKfHkVt0lQ3pxRs+OPPrBJsrtVhsmTjTFUUNzF82RwxggR1W3qGoZMB+YGVZmJvCyOlYAySLSU1V3q+qXAKp6CNgA9A6pU3lt+EvAZSHD56tqqapuBXLcGEyUnlqcQ+/kNlw2vHfthWvh8wl3Tx3A9sIS3sje0QDRmZo0h2ZeE3vqs11Fkxx6A6F7hTz+sYOPuoyIpAHDgS/cQT1UdTeA+7d7HeaHiNwmItkikl159aOBrNxCsnKL+NHEdBL9DXO+wbkDuzOsTzLzPt1MIGg7r8aSlJREQUGBJQjToCqf55CUlFSnetFc5xCpFzJ8662xjIi0B94C7lTV2i77jGZ+qOqzwLPgdEjXMs0W45lPN9OlXSuuGh35LI76EBFmn9OP2X/+kvfX7WHG4Kq3AzAnLzU1lby8POzHjmlolU+Cq4tokkMeEHqSfCqwK9oyIpKIkxheUdX/Cymzt7LpSUR6Avtqm5ap2Y7CEj76dh9zzu1Pm1YN200zJfMU0rq25ZklW5h+5in23IFGkJiYWKcndRnTmKJpd8gCMkQkXURa4XQWLwgrswC4wT1raRxwwN3pC/A8sEFV/ztCnRvd9zcCb4cMv1pEWotIOk4n98o6L1kLND9rOwJcPabhjhoq+X3CLRP78dWOYlZuLWzw6RtjYkutyUFVK4A5wPs4HcpvqOo6EZktIrPdYguBLTidx88BP3GHnw1cD5wnImvc1wx33MPAFBHZBExxP6Oq64A3gPXAe8Adqho4+UVt3soDQV7PyuO807vTO7lxLor5/ohUOrVJ5OUVJ387YGNMbIvq3kqquhAnAYQOmxfyXoE7ItRbRuQ+BFS1ADi/mnEPAQ9FE5txfLB+L/sPl3Lt2FMbbR5tWvn53ohU/rQil/2HS+nWvnWjzcsY4y27fUYzMT9rB706JTFpQEqjzueasX0oDyhvrspr1PkYY7xlyaEZyD9UyrJN+Vw2vHej3+Kif/cOjEnvwmsrtxO001qNabYsOTQD76zdRVDh8ga46C0a147ty7aCEpZvLmiS+Rljmp4lh2bgr6t3ktmzIxk9OjTJ/KYOOoXktonMz9reJPMzxjQ9Sw5xbuv+I3yVd4DLhvdqsnkmJfq5ZEgvPli/l0PHyptsvsaYpmPJIc79dfVORODSoU3TpFTpsuG9Ka0I8v66vU06X2NM07DkEOf+tnYX49K7ckqnut035WSN6JvMqV3b8pfVdtaSMc2RJYc4lrPvEFvyjzBj8ClNPm8R4bJhvVm+uYA9B441+fyNMY3LkkMcq2zSmZLZ9MkBnKYlVexBQMY0Q5Yc4th73+xheN/kJm9SqpTerR1D+yTz9hq7L6IxzY0lhzi1s/goX+88wNRB3hw1VLp4cE/W7TrI9oIST+MwxjQsSw5xatG6PQCeJ4dpZzrzf/eb3Z7GYYxpWJYc4tR73+xhYI8OpHdr52kcfbq0ZUhqJxZ+bcnBmObEkkMcOlBSTva2IqZk9vA6FABmDO7JV3kHyCuypiVjmgtLDnFoaU4+gaBy7umNewfWaE13m5be+2aPx5EYYxpKVMlBRKaJyEYRyRGReyOMFxF5zB2/VkRGhIx7QUT2icg3YXVeD3kAUK6IrHGHp4nI0ZBx8zAnWLwxn05tEhnWp7PXoQBwatd2DOrV0ZqWjGlGak0OIuIHngSmA5nALBHJDCs2HedxnhnAbcDTIeNeBKaFT1dVr1LVYao6DOcZ06HPl95cOU5VZ4fXbcmCQWXxxnwmDUhp9Ntz18WMwT35cnsxuw8c9ToUY0wDiObIYQyQo6pbVLUMmA/MDCszE3hZHSuAZBHpCaCqS4BqHzrsPmf6SuC1+ixAS7N+90H2Hy5lciM/1Keupg5y+j8+2rDP40iMMQ0hmuTQG9gR8jnPHVbXMtWZCOxV1U0hw9JFZLWIfCoiEyNVEpHbRCRbRLLz8/OjnFX8++RbZ+d7zsDYSg6npbTn1K5t+XCD3YjPmOYgmuQQqe0i/BFg0ZSpzixOPGrYDfRV1eHAXcCrItKxysRVn1XVUao6KiUltnaUjemTjfsYktop5p7fLCJccEYPlm8u4EhphdfhGGNOUjTJIQ/oE/I5FQi/X0I0ZaoQkQTgCuD1ymGqWqqqBe77VcBmYEAUcTZ7RUfKWLOjmMkDu3sdSkTnn9GdsoogSzft9zoUY8xJiiY5ZAEZIpIuIq2Aq4EFYWUWADe4Zy2NAw6oajSnrlwAfKuqx+/7LCIpbic4ItIPp5N7SxTTavaWbMonqDA5xpqUKo1O60LHpAQ+sqYlY+JeQm0FVLVCROYA7wN+4AVVXScis93x84CFwAwgBygBflhZX0ReAyYD3UQkD3hAVZ93R19N1Y7oScCvRaQCCACzVbXaDu2W5NPv8uncNpGhqclehxJRot/H5IHd+fjbfQSCGlNnUxlj6qbW5ACgqgtxEkDosHkh7xW4o5q6s2qY7k0Rhr2Fc2qrCaGqfJazn7P6d4vpne4FmT1Y8NUu1uwoZuSpsXEdhjGm7uwK6TixOf8Iew+WcvZp3bwOpUbnDEghwSd21pIxcc6SQ5xYvtnp5D27f1ePI6lZpzaJjE7rYv0OxsQ5Sw5xYtmm/fRObkPfLm29DqVWF2T24Lu9h+0ZD8bEMUsOcSAQVFZsKeDs/l1xLiiPbeed7pxqu/g7u1ramHhlySEOfLPzAAePVXB2/9jub6iU3q0dp3Zty6cbW86V68Y0N5Yc4sBnbn/DWTHeGR3qnAEpLN9cwLHygNehGGPqwZJDHFieU8DAHh1I6RBbt8yoyTkDUjhaHiA7t8jrUIwx9WDJIcYdKw+QlVvIWTF+llK48ad1pZXfx+KN1u9gTDyy5BDjvtxeRGlFMOavbwjXtlUCY9K78Ol31u9gTDyy5BDjlucU4PcJY/t18TqUOps8MIVN+w6zs9geAGRMvLHkEOM+31LA4N6d6JCU6HUodXaO+0AiO2vJmPhjySGGHSsPsDavmLHp8XfUANC/e3t6J7fhU7vewZi4Y8khhn21o5jygDIqLT6Tg4gwaUAKn+UUUFYR9DocY0wdWHKIYVm5zp3KR8Xx3U0nD0zhcGkFX263U1qNiSeWHGJYVm4RA3q0p3O7Vl6HUm9nndaVBJ/YWUvGxBlLDjEqEFS+3FYUt01KlTokJTIqrTOLrVPamLgSVXIQkWkislFEckTk3gjjRUQec8evFZERIeNeEJF9IvJNWJ25IrJTRNa4rxkh4+5zp7VRRKaezALGqw27D3KotIIxcZ4cACZmpLBh90HyD5V6HYoxJkq1Jgf3ec5PAtOBTGCWiGSGFZuO86znDOA24OmQcS8C06qZ/KOqOsx9LXTnl4nz+NBBbr2nKp8p3ZJku/0No+P0TKVQEzOcC/g+y9nvcSTGmGhFc+QwBshR1S2qWgbMB2aGlZkJvKyOFUCyiPQEUNUlQF2eAT0TmK+qpaq6Fee51GPqUL9ZyMotolenJHont/E6lJM2qFcnktsmsnSTJQdj4kU0yaE3sCPkc547rK5lIpnjNkO9ICKVp+RENS0RuU1EskUkOz+/ebVnqyorcwubxVEDgN8nnH1aN5bl5OM8btwYE+uiSQ6Rni4T/h8eTZlwTwOnAcOA3cAjdZmWqj6rqqNUdVRKSkots4ov2wtLyD9Uyuhm0N9QaUJGN/YeLCVn32GvQzHGRCGa5JAH9An5nArsqkeZE6jqXlUNqGoQeI5/NB3VeVrNzcqtTivcmGZy5AAwwX1QkTUtGRMfokkOWUCGiKSLSCuczuIFYWUWADe4Zy2NAw6o6u6aJlrZJ+G6HKg8m2kBcLWItBaRdJxO7pVRxNlsZOUW0qlNIv1T2nsdSoPp06Ut6d3ascw6pY2JCwm1FVDVChGZA7wP+IEXVHWdiMx2x88DFgIzcDqPS4AfVtYXkdeAyUA3EckDHlDV54H/FJFhOE1GucCP3emtE5E3gPVABXCHqraox4ll5xYxOq0zPl/sPy+6Lib078ZbX+ZRVhGkVYJdYmNMLKs1OQC4p5kuDBs2L+S9AndUU3dWNcOvr2F+DwEPRRNbc5N/qJQt+49w1eg+tReOMxMyuvGnFdtYvb2Isf3i6+FFxrQ09vMtxlRe3xDvV0ZHMv60rvh9Yv0OxsQBSw4xJiu3iKREH4N7d/I6lAbXMSmRoamdWGr9DsbEPEsOMSYrt5BhfZKbbZv8hIwUvs4r5kBJudehGGNq0Dz3QHHqcGkF63YdaFbXN4SbmNGNoMLyzXb0YEwss+QQQ1ZvLyKoNOvkMKxPMu1bJ1jTkjExzpJDDMnaWohPYEQcP9ynNol+H+P6dWGZdUobE9MsOcSQlbmFZPbqSPvWUZ1hHLcm9O/G9sIStheUeB2KMaYalhxiRFlFkDU7ipt1k1KlCRnOvbCW5jSvGyYa05xYcogR3+w6wLHyYLN4uE9tTktpR89OSda0ZEwMs+QQI7K2Nt+L38KJCBP6d+OznP0EgnYLb2NikSWHGJGVW0R6t3akdGjtdShNYkJGNw4eq2BtXrHXoRhjIrDkEAOCQSV7WyGj05rvWUrhKm/hbU1LxsQmSw4xICf/MMUl5S2iM7pS1/atGdSro13vYEyMsuQQA7Lcm+21pOQATtPS6u1FHC6t8DoUY0wYSw4xIGtrISkdWnNq17Zeh9KkJmWkUB5QvthS4HUoxpgwUSUHEZkmIhtFJEdE7o0wXkTkMXf8WhEZETLuBRHZJyLfhNX5vYh865b/i4gku8PTROSoiKxxX/No5rJyixiT1gWR5vVwn9qMPLUzSYk+u4W3MTGo1uQgIn7gSWA6kAnMEpHMsGLTcR7nmQHcBjwdMu5FYFqESX8AnKmqQ4DvgPtCxm1W1WHua3aUyxKXdhYfZWfxUUa1oM7oSkmJfsakd2XpJrsYzphYE82RwxggR1W3qGoZMB+YGVZmJvCyOlYAyZXPiFbVJUBh+ERVdZGqVjY2rwBS67sQ8Sy7hfY3VJqU0Y3N+UfYVXzU61CMMSGiSQ69gR0hn/PcYXUtU5ObgXdDPqeLyGoR+VREJkaqICK3iUi2iGTn58fvL8+VWwtp3zqBM3p29DoUT0zIsFNajYlF0SSHSA3h4Ze1RlMm8sRFfgVUAK+4g3YDfVV1OHAX8KqIVNlzquqzqjpKVUelpKREM6uYlJ1bxIhTO+P3taz+hkoDe3QgpUNrlljTkjExJZrkkAeEPu0+FdhVjzJViMiNwMXAtaqqAKpaqqoF7vtVwGZgQBRxxp3ikjI27j3EmBbY31BJRJiY4dxKI2i30jAmZkSTHLKADBFJF5FWwNXAgrAyC4Ab3LOWxgEHVHV3TRMVkWnAPcClqloSMjzF7QRHRPrhdHJviXqJ4kh2bhHQcvsbKk3M6EZRSTnrdh30OhRjjKvWBweoaoWIzAHeB/zAC6q6TkRmu+PnAQuBGUAOUAL8sLK+iLwGTAa6iUge8ICqPg88AbQGPnBP4Vzhnpk0Cfi1iFQAAWC2qlbp0G4OsrYVkugXhvZJ9joUT53t3kpjyaZ8Bqd28jiaRjDXo2Wae8Cb+ZpmIaqnyqjqQpwEEDpsXsh7Be6opu6saob3r2b4W8Bb0cQV77K2FjIkNZmkRL/XoXiqe4ckTj+lA8s27eeOcyNuFsaYJmZXSHvkWHmAr3ceaPFNSpUmDUghe1shJWV2Kw1jYoElB4+s2VFMeUBb1J1YazIxo5tzK42tzbIF0Zi4Y8nBI1lbCxGBUafakQM4nfKtEnws/c6udzAmFlhy8MjK3EIG9uhAp7aJXocSE5IS/YxN78Iye660MTHBkoMHKgJBvtxWZP0NYSb078Z3ew+z58Axr0MxpsWz5OCBb/cc4khZoEXebK8mEzOcK93tRnzGeM+SgwdWup2uY9LtyCHU6ad0oFv7Viyzp8MZ4zlLDh7I3lZIauc29OzUxutQYorPJ0zo341lm+xWGsZ4zZJDE1NVVm61/obqTMxIoeBIGRv22K00jPGSJYcmlltQwv7DpZYcqlF5C297Opwx3rLk0MSyjvc3WGd0JD06JjGwRwfrlDbGY5YcmtjK3EK6tGvFaSntvQ4lZk3I6EZWbhFHywJeh2JMi2XJoYmt3FrI6LTOuHeiNRFMzOhGWUWQlbl2Kw1jvGLJoQntOXCM7YUl1t9Qi7HpXWnl97H0O2taMsYrlhyaUOUv4bHpXT2OJLa1aeVndHpnPrXkYIxnLDk0oaythbRr5eeMnh28DiXmnTuwO5v2HSavqKT2wsaYBhdVchCRaSKyUURyROTeCONFRB5zx68VkREh414QkX0i8k1YnS4i8oGIbHL/dg4Zd587rY0iMvVkFjCWrNxayMi0LiT4LSfXZvLA7gAs3mhHD8Z4oda9lPs85yeB6UAmMEtEMsOKTcd51nMGcBvwdMi4F4FpESZ9L/CRqmYAH7mfcad9NTDIrfdU5TOl41lxSRkb9x5ijN1PKSqnpbQjtXMbSw7GeCSan7BjgBxV3aKqZcB8YGZYmZnAy+pYASSLSE8AVV0CRDrtZCbwkvv+JeCykOHzVbVUVbfiPJd6TB2WKSZl5RYBMMb6G6IiIpw7sDuf5eyntMJOaTWmqUXzDOnewI6Qz3nA2CjK9AZ21zDdHqq6G0BVd4tI95BprYgwrROIyG04Ryn07du39qXwWFZuIa38PoakevSw+XBePfQeon7w/bmnp/CnFdtYubXw+B1bT26+MbLujYkD0Rw5RDohP/yuaNGUiVZU01LVZ1V1lKqOSklpgB1HI/tiayHD+iSTlBj3LWRNZny/brRK8FnTkjEeiCY55AF9Qj6nArvqUSbc3sqmJ/fvvpOYVkw7UlrBup0HGG23zKiTNq38jOvXlU827qu9sDGmQUWTHLKADBFJF5FWOJ3FC8LKLABucM9aGgccqGwyqsEC4Eb3/Y3A2yHDrxaR1iKSjtPJvTKKOGPW6u3FVATV+hvq4dyBKWzJP8L2Ajul1ZimVGtyUNUKYA7wPrABeENV14nIbBGZ7RZbCGzB6Tx+DvhJZX0ReQ34HBgoInkicos76mFgiohsAqa4n1HVdcAbwHrgPeAOVY3rHsmVuYX4BEb0TfY6lLhz/JTW7+zowZimFE2HNKq6ECcBhA6bF/JegTuqqTurmuEFwPnVjHsIeCia2OLByq0FZPbqSIekRK9DiTvp3dqR1rUtn3y7jxvGp3kdjjEthl2N1cjKKoKs3l7MmDRrUqqvyQO7s3xzAcfK4/oA0pi4YsmhkX29s5jSiqA9v+EkTB6YQmlFkBVbCrwOxZgWw5JDI1u51bn4ze7EWn/j+nUlKdHHJ99av4MxTcWSQyP7fEsBGd3b07V9a69DiVtJiX4m9O/Ghxv24XRvGWMamyWHRlQeCJKdW8j406y/4WRdcEYPdhYfZcPuQ16HYkyLENXZSqZ+1uYVU1IWYHw/Sw4nqMdtLM7XjghP8eFTPyUz4S+NEJQxJpQdOTSi5TlOB+pYSw4nLUUOMlxy+CAw0utQjGkRLDk0os+3FHBGz450adfK61CahQv8q/ha+7FbrXPfmMZmyaGRlFYEWLWtyJqUGtCFvlUAfBgYUUtJY8zJsuTQSFZvd65vsM7ohnOa7CJddvNB0JqWjGlslhwayeebC/AJjEm3JpCGIgIX+FbxeXAQh7SN1+EY06xZcmgkn28pYFCvTnRqY/dTakhT/KsoJ4ElwSFeh2JMs2bJoREcLQuwensRZ1mTUoMbIZvozCE7a8mYRmbJoRGs2lZEeUAZZ8mhwSVIkPN8X/JxcDjlak/VM6axWHJoBJ9v2Y/fJ3Y/pUYyxb+Kg7QjKzjQ61CMabaiSg4iMk1ENopIjojcG2G8iMhj7vi1IjKitroi8rqIrHFfuSKyxh2eJiJHQ8bNC59frPt8cwFDUjvRvrVdgN4YzvGtpQ3HWBgc63UoxjRbtSYHEfEDTwLTgUxglohkhhWbjvM4zwzgNuDp2uqq6lWqOkxVhwFvAf8XMr3NleNUdTZx5OCxcr7KO2D9DY2ojZRxru8r3guMJqDidTjGNEvRHDmMAXJUdYuqlgHzgZlhZWYCL6tjBZAsIj2jqSsiAlwJvHaSyxITPt9cQCCoTMpI8TqUZm2GfwX7SSZLrWnJmMYQTXLoDewI+ZznDoumTDR1JwJ7VXVTyLB0EVktIp+KyMRIQYnIbSKSLSLZ+fn5USxG01i6KZ92rfwM72sP92lM5/rWkEQpCwPjvA7FmGYpmuQQ6bg9/Kb61ZWJpu4sTjxq2A30VdXhwF3AqyLSscpEVJ9V1VGqOiolJXZ+pS/dtJ/xp3WlVYL19TemdlLKZN9XvBsYTdCaloxpcNHswfKAPiGfU4FdUZapsa6IJABXAK9XDlPVUlUtcN+vAjYDA6KI03PbCo6wraCEidak1CRm+L8gn85ka1xsHsbElWiSQxaQISLpItIKuBpYEFZmAXCDe9bSOOCAqu6Oou4FwLeqmlc5QERS3I5sRKQfTif3lnouX5Naumk/ABMzunkcSctwnm81rSljYcDOWjKmodWaHFS1ApgDvA9sAN5Q1XUiMltEKs8kWoizA88BngN+UlPdkMlfTdWO6EnAWhH5CngTmK2qhfVcvia1dFM+vZPbkN6tndehtAjt5Rjn+L7i3cAYa1oypoFFdSK+qi7ESQChw+aFvFfgjmjrhoy7KcKwt3BObY0rFYEgy3MKuHhoT5wTsExTuMj/BYuCo1mlAxgtG70Ox5hmw3pNG8hXecUcKq2w/oYmdr7vS5IoZUFgvNehGNOsWHJoIEu+249PsIvfmlh7OcYFvi95JzDe7rVkTAOy5NBAlm7KZ0hqMslt7ZGgTe1y/zKK6GC38TamAVlyaABFR8pYs6OYSXaWkicm+dbSmUP8NXC216EY02xYcmgAn36XT1Dh/DN6eB1Ki5QoAS7yr+CD4EgOa5LX4RjTLFhyaAAffbuPbu1bM7h3J69DabEu9y/jGK15Pzja61CMaRYsOZyk8kCQTzfu49yBKfh8dgqrV0bIJvrIPmtaMqaBWHI4Sau2FXHwWAXnn9Hd61BaNBGY6fuMz4Jnsk/tCM6Yk2XJ4SR9/O0+Ev3CBLu+wXOX+T8jiI+37ejBmJNmyeEkfbRhL+P6dbWnvsWA/r5dDJdNvB6YjIbf+9cYUyeWHE5C7v4jbM4/wnmnW5NSrLja/wk5msqXmuF1KMbENUsOJ+Hjb/cBWHKIIRf5V9CWY7wemOx1KMbENUsOJ+HDDXvp3709p3a1u7DGivZyjEv8n/NOYLxd82DMSbDkUE+FR8r4YmshUwfZhW+x5kr/YkpI4h17hKgx9Wa9qPX04fq9BILK9DN7eh2KCTNCNtFf8ng9cC5XJyz2OpyWZ65HpxLPPeDNfJspO3Kop3e/2U1q5zYM6lXl8dbGYyJwtX8xqzWDDcE+tZY3xlQVVXIQkWkislFEckTk3gjjRUQec8evFZERtdUVkbkislNE1rivGSHj7nPLbxSRqSe7kA3t4LFyluXsZ9qgU+zBPjHq+/5PaU0ZLwcu9DoUY+JSrcnBfZ7zk8B0IBOYJSKZYcWm4zzrOQO4DXg6yrqPquow97XQrZOJ8/jQQcA04KnKZ0rHik++3Ud5QJk++BSvQzHVSJYjXOb/jL8GzuaA2gkDxtRVNEcOY4AcVd2iqmXAfGBmWJmZwMvqWAEki0jPKOuGmwnMV9VSVd2K81zqMXVYpkb37td76N6hNcP7dPY6FFODG/yLOEoS/xs4x+tQjIk70SSH3sCOkM957rBoytRWd47bDPWCiFTuaaOZHyJym4hki0h2fn5+FIvRMErKKlj83T6mDjrFbrQX4wb5tjFavuXlwBQCat+VMXURTXKI9F8VfnOC6srUVPdp4DRgGLAbeKQO80NVn1XVUao6KiWl6e5r9MH6vRwrDzJjsJ2lFA9uTHif7dqDxcFhXodiTFyJJjnkAaGnfKQCu6IsU21dVd2rqgFVDQLP8Y+mo2jm55m31+yiZ6ckxqZ38ToUE4Wpvmx6UMiLgZg7r8GYmBZNcsgCMkQkXURa4XQWLwgrswC4wT1raRxwQFV311TX7ZOodDnwTci0rhaR1iKSjtPJvbKey9egCo+UseS7fC4d2sualOJEogS4IeEDlgaHsD7Y1+twjIkbtSYHVa0A5gDvAxuAN1R1nYjMFpHZbrGFwBaczuPngJ/UVNet858i8rWIrAXOBX7m1lkHvAGsB94D7lDVQEMs7Mn6+9e7qQgqM4dV6QIxMew6/4e04yjPVFzidSjGxI2orpB2TzNdGDZsXsh7Be6Itq47/Poa5vcQ8FA0sTWlt1fvZECP9pzRs4PXoZg66CRHuMb/ES8EpvOL4Bv08TXdCQzGxCu7QjpKOwpLyN5WxMxhve3Ctzh0S8K7+FCeC1zkdSjGxAVLDlFa8JXTJ37p0F4eR2Lq4xQp4jL/Mt4InEOB2pGfMbWx5BCFYFB5I3sHY9K70KdLW6/DMfX0Y/87lJLIHyumeR2KMTHPkkMUVmwtYFtBCbPG2E3c4ll/3y6m+7J4MTCVIm3vdTjGxDRLDlF4PWsHHZIS7PbczcCdCW9yhCSeqbjY61CMiWmWHGpRXFLGu9/s4fLhvUlKjKn7/5l6GODbyaW+z3kpcCH5ardbN6Y6lhxq8dfVOymrCHLVaGtSai7+OeEtykjk6YpLvQ7FmJhlyaEGqsorX2xncO9ODOrl0dOtTIPr59vDFf6l/DlwAXnazetwjIlJlhxq8FlOAZv2HebGs9K8DsU0sJ8lvIkA/1l+ldehGBOTLDnU4I+fbaVru1ZcPMQ6opubXlLIbf6/syB4Nl8G+3sdjjExx5JDNXL3H+Hjjfu4dmxf64hupmYnLKA7Rfx7+fVolZvCG9OyWXKoxsufb8MvwrXjTvU6FNNI2kkpv0h4g9WawdvBs70Ox5iYYskhguKSMl7P2s7FQ3rSo2OS1+GYRvQ9/xKGyGZ+U36tPWvamBCWHCL442e5HCkLcPtka4tu7vyi/DbxDxTSkd9VWOe0MZUsOYQ5XFrBi8tzmZLZg4Gn2A3aWoIzfdu42f8urwYuIDs4wOtwjIkJUSUHEZkmIhtFJEdE7o0wXkTkMXf8WhEZUVtdEfm9iHzrlv+LiCS7w9NE5KiIrHFf88Ln15heWbGNA0fLmXOuHTW0JD9LeJPe5HNf+a0c00SvwzHGc7UmBxHxA08C04FMYJaIZIYVm47zOM8M4Dbg6SjqfgCcqapDgO+A+0Kmt1lVh7mv2TSRkrIKnlu6lYkZ3RjaJ7mpZmtiQDsp5beJz7NJU/mviiu9DscYz0Vz5DAGyFHVLapaBswHZoaVmQm8rI4VQLL7jOhq66rqIvcxogArgNQGWJ6T8vzSrew/XMqdF1jTQkt0jn8t1/sX8YfARSwPhP/+MaZliSY59AZ2hHzOc4dFUyaaugA3A++GfE4XkdUi8qmITIwUlIjcJiLZIpKdn3/yj30sOFzKM0u2MHVQD0ae2vmkp2fi0y8TXqWf7OLn5bfb2UumRYsmOUR6Jmb4JUPVlam1roj8CqgAXnEH7Qb6qupw4C7gVRGpcvtMVX1WVUep6qiUlJRaFqF2T3ySQ0lZBXdPPf2kp2XiVxsp49HEp8inEz8vn01Q7ZGwpmWKJjnkAaG3JE0FdkVZpsa6InIjcDFwrapzjaqqlqpqgft+FbAZaNR2ni35h/nzim1cOaoP/bvbQ2BauqG+Lfxrwp/5MDiSpwN251bTMiVEUSYLyBCRdGAncDVwTViZBcAcEZkPjAUOqOpuEcmvrq6ITAPuAc5R1ZLKCYlIClCoqgER6YfTyb3lZBayJqrK/W+vIylwhLvWXgxfH2isWZk4cqN/EauDGfxXxQ8YLFuY5P/a65Dqbq7dSdjUX61HDm6n8RzgfWAD8IaqrhOR2SJSeSbRQpwdeA7wHPCTmuq6dZ4AOgAfhJ2yOglYKyJfAW8Cs1W18OQXNbJ31u5mWc5+fpHwBt3FEoNxiMB/JP6BAZLHnPKfsikYqavMmOZLtBnccWzUqFGanZ1d53qHjpVz/iOf0r1ja97efwl+if91YRrWjmA3Li/7Na0p5y+tH6C7FHsdkqnOXPtxV1ciskpVR0Ua16KvkN607zBBVX5z2WBLDCaiPr79/LHV7ymiAzeX3c0hbeN1SMY0iRadHEb07cyye85jmF3wZmow2LeVJxMf41vtw01l/8JhtZsxmuavRScHwJ7VYKJyrn8Njyc+zhrtz01l91iCMM1ei08OxkRruj+LxxMfZ7X257qy+yhQuzGjab4sORhTBzP8K3kq8X/YoKfy/bK5bA929zokYxqFJQdj6miqP5tXWz1EkbbnirIHWRG0q+pN82PJwZh6GOnbxFut5tJRjnBt2a94rmKGPYfaNCuWHIypp9N8u3m71b8xxbeKhyqu45byX7BXk70Oy5gG0aIvgjuB3WrA1JMqvBiYyu8qrqYVFTyY+CKX+T5D7J59pjE1wEV/dhGcMY1IBH6Y8D4LW91Hf9nJz8rvYFb5r1gf7Ot1aMbUmyUHYxpIP98e/rfVg/x7wgtsDPbl4rLfcm/5reRpN69DM6bOLDkY04D8olyf8CGLW9/Fjf73eSswicml/82/lP+IzcGeXodnTNQsORjTCDrJER5I/BOftr6T6/wf8nbgbM4ve4Rry37JwsAYytWuzDexLZrnORhj6qmXFDI38WXuSPgrrwfO5bWK8/hJ8E46c4gL/dnM8H3BWb51JErA61CNOYElB2OaQIocZE7C29zuX8CnwaG8HTibdwLjeD1wLh05wlm+dZzt+4azfOvoJ7vtTCfjOUsOxjQhvyjn+ddwnn8NxzSRpcHBLAqOYnlgEO8FxwDQjWKG+LYwWLYyxLeF033b6UkhPrutvGlCUSUH95Ge/wP4gT+o6sNh48UdPwMoAW5S1S9rqisiXYDXgTQgF7hSVYvccfcBtwAB4Keq+v5JLaUxMShJypni/5Ip/i/RBNimPfgsOIhVwQF8o+ksDg4jGHC6BZMoJU320E/2kCZ76CX76SFFnCJF9JAiunLAnkliGlStyUFE/MCTwBQgD8gSkQWquj6k2HScZz1n4DxD+mlgbC117wU+UtWHReRe9/M9IpKJ86zpQUAv4EMRGaCq1ihrmi0RSJO9pPn2ci0fA1CirVmnaXwX7M1W7clW7cm32odFwZFUhP3r+gnQiSMky2E6coROcoRO7t8OHCVJymhDKUmU0Zpy2ojzvg1ltJZyEqgggSB+AiQcfwXxi/PeT/D4Xz9BJ2YUQU94L5HeWxNZXIrmyGEMkKOqWwBEZD4wEwhNDjOBl9W53HqFiCSLSE+co4Lq6s4EJrv1XwIWA/e4w+eraimwVURy3Bg+r/9iGhN/2kopo2Ujo30bTxheoT4K6Mge7cJe7Xz8VUQHDmg7DtCOYu3ANk7hQLAdh2hTJZl4QQhGTB6c8L4p4mi6I6xGzYv3vwfA9DN78siVQxt88tFsMb2BHSGf83CODmor07uWuj1UdTeAqu4Wkcp7H/cGVkSY1glE5DbgNvfjYRHZGF6miXUD9nscQ23iIUaIjzg9jrEY2B5NwXhYlxAfccZkjBuA/77qhEF1ifPU6kZEkxwiJb/w1FtdmWjq1md+qOqzwLO1TKvJiEh2dfcoiRXxECPER5zxECNYnA0pHmKEhoszmovg8oA+IZ9TgV1Rlqmp7l636Qn37746zM8YY0wjiiY5ZAEZIpIuIq1wOosXhJVZANwgjnHAAbfJqKa6C4Ab3fc3Am+HDL9aRFqLSDpOJ/fKei6fMcaYeqi1WUlVK0RkDvA+zumoL6jqOhGZ7Y6fByzEOY01B+dU1h/WVNed9MPAGyJyC07j6Q/cOutE5A2cTusK4I44OVMpZpq4ahAPMUJ8xBkPMYLF2ZDiIUZooDibxfMcjDHGNCy78Z4xxpgqLDkYY4ypwpJDFETkBRHZJyLfhAybKyI7RWSN+5oRMu4+EckRkY0iMtXjOF8PiTFXRNa4w9NE5GjIuHlNFGMfEflERDaIyDoR+Wd3eBcR+UBENrl/O4fUafL1WUOcvxeRb0VkrYj8RUSS3eFNvj5riDGmts0a4oyZbVNEkkRkpYh85cb4oDs81rbL6uJs+O1SVe1VywuYBIwAvgkZNhf4RYSymcBXQGsgHdgM+L2KM2z8I8D97vu06so1cow9gRHu+w7Ad+46+0/gXnf4vcDvvFyfNcR5IZDgDv9dSJxNvj5riDGmts3q4oylbRPn+qr27vtE4AtgXAxul9XF2eDbpR05REFVlwCFURY/fvsPVd2KcwbXmEYLLkRNcYqIAFcCrzVFLNVR1d3q3pRRVQ/hXODZG2e9veQWewm4zH3vyfqsLk5VXaSqFW6xFTjX4XiihnVZnZhal5XjY2HbVMdh92Oi+1Jib7uMGGdjbJeWHE7OHPcw7oWQw83qbiXitYnAXlXdFDIsXURWi8inIjKxqQMSkTRgOM6vnxNupwKE3k7F0/UZFmeom4F3Qz57tj4jxBiT22Y16zImtk0R8btNW/uAD1Q1JrfLauIM1SDbpSWH+nsaOA0YBuzGOSyG+t0ypCnM4sRfZruBvqo6HLgLeFVEOjZVMCLSHngLuFNVD9ZUNMKwJluf1cUpIr/CuQ7nFXeQZ+szQowxuW3W8J3HxLapqgFVHYbzq3uMiJxZQ3HP1mVNcTbkdmnJoZ5Uda/7JQWB5/jHIWXM3f5DRBKAK3CenwGAezhc4L5fhdNmOqCJ4knE2Um8oqr/5w6OudupVBMnInIjcDFwrboNu16tz0gxxuK2WcO6jKlt051nMc5doqcRg9tlpbA4G3y7tORQT5UbjOtyoPIMoVi8/ccFwLeqmlc5QERSxHneBiLSDyfOLY0diNu+/DywQVX/O2RUTN1Opbo4xXl41T3ApapaEjK8yddnDTHG1LZZw3cOMbJtuvNMdt+3qYyL2NsuI8bZKNtlQ/akN9cXziHvbqAc5xfDLcCfgK+BtTgbSs+Q8r/CydAbgelexukOfxGYHVb2e8A6nDMuvgQuaaIYJ+Acfq8F1rivGUBX4CNgk/u3i5frs4Y4c3DamiuHzfNqfdYQY0xtm9XFGUvbJjAEWO3G+A3/OHMq1rbL6uJs8O3Sbp9hjDGmCmtWMsYYU4UlB2OMMVVYcjDGGFOFJQdjjDFVWHIwxhhThSUHY4wxVVhyMMYYU8X/B6CH5shPx19RAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "#-----------\n",
    "#A way to check this test is to see if a chi squared distribution does look like the density given by the bootstrapped samples from above:\n",
    "#-----------\n",
    "devs = []\n",
    "for i in range(boot_pred.shape[2]):\n",
    "    devs.append(resid_dev(boot_pred[:,1,i].reshape(-1,1), boot_response[:,1,i].reshape(-1,1))[0])\n",
    "\n",
    "x = np.linspace(np.min(devs)-1, np.max(devs)+1, 1000)\n",
    "plt.plot(x, chi2.pdf(x, 180), label='Predicted pdf w/ 180 df')\n",
    "plt.hist(devs, density = True)\n",
    "plt.legend(loc = 'best')\n",
    "plt.title(\"Predicted Chi2 pdf and Observed Density Histogram\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEICAYAAAC0+DhzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABF3UlEQVR4nO3deXwU9f348dd7NwnhDkdAICBBA3LIfcohiHJ4QKUWxfsuKv221Vq1/X4ttrWX9WfrURG864F4VbRcgiKgcioi4b4JEAh3uHLsvn9/zASXzbUJSWaTvJ+PRx7ZnfnMzHtmd+a98/nMzEdUFWOMMSaUz+sAjDHGRB9LDsYYY/Kx5GCMMSYfSw7GGGPyseRgjDEmH0sOxhhj8rHkUIZE5FUR+aP7eqCIrK+g5aqInF/KaW8VkUVFjJ8pIreUProSxXJ6PUSkpoh8LCJHROTdcl5ua3fZMWcxj4ki8kZZxlXeivvsvSIivxGRFytweTeIyJyKWl5lUe2Sg4hsE5GTInJMRPaKyCsiUqesl6OqC1W1XQTxlPsOKiLDRWSBiGSKSIaIfCEioyKZVlVHqupr7nyuEJFFInJYRNJFZIqI1C2nsK8BmgKNVPUn5bSMiLmf0/cicsJd9+dFJMHruMpLSMI8FrKvfCIil5X3slX1T6p6Z1gcpUrche1f7nHgUnd5b6rqsAjmdfrHX3VQ7ZKD6ypVrQN0B3oB/xte4Gx+RUYTEbkGeBd4HUjCOeA+ClxVitnVB/4INAfau/N7omwizedcYIOq5pbT/CMmIg8AfwUexNkGfXHi+1RE4iowDi++kwnuvtIF+BT4UERu9SCOKi0qjzeqWq3+gG3ApSHvnwA+cV8rcB+wEdjqDrsSWAkcBr4COodM2w34BsgE3gGmAn90xw0G0kLKtgQ+ADKAA8CzOAfYU0AAOAYcdsvWAP4O7AD2ApOAmiHzehDYA+wGbnfjPr+AdRV3Hg8WsT1uBRa5yzsEbAVGhoyfD9xZyLRjgO+L2daPAGvceb8CxBe3HsBjQDaQ426XOwqYd2/ga/dz2eNuz7iQ8QqMdz/LQ8BzgLjj/O767ge2uJ+5AjEFLKeeG8PYsOF1gH3A7e77icB77vcg0/1edAkp/xCwyx23HhjqDvcBDwOb3e/FNKChO661G9cd7ue4AJgFTAiL5TtgjPv6ApyD+EF3OWNDyjUCpgNHgaXAH4BFhXx2ecuOCRv+K5zvpM993xx4H+d7vRX4n5CyE931ed1d71SgZwTbZCLwhvt6hxvHMffvYnfdLgyZTxPgJJBY2Pe7qONAaBmcfeYp97M9AqwCOgF343wfs904PnbLt8fZRw676zcqbHt/7G7vZTg/rBaFjC/oePNPYKc7zQpgYNj2fBd4w91m3wNtcfaxfe50w8rsWFlWM6osf2FfipbuB/qHkA/rU6AhUBPnzGIf0AfngHKLO30NIA7YDvwSiMWpBsmhgOTgTvud+6WrDcQDAwr78gL/wNmJGwJ13S/Yn91xI3B2zk7uvN6i8ORwgTsuuYjtcasb911unPfgHKzzDqTzKTw5/AOYWsy2Xu1u54bAlyHbp8j1IOQAUci8e+D8go/BOZCtBX4RtuN9AiQArXAOXiPcceOBdSFxfU7hyWEEkFvIuNeAt0PizXG/B7E4B9Gt7ut2ODtuc7dsa+A89/UvgMU4Z2E1gBdC5tnajet1dxvVBG4GvgyJoQPOgamGW2YncJu7XbrjJMCObtmpOAfr2u5230XJk0Mbd3h7nMS2AudMNM4dtwUYHrJNTgGX43y3/gwsdscVtU1Of/YFxQH8C/hryPuf4x6sC/l+lyQ5DHfXKQEnUbQHmrnjXsX9/rrvY4FNwG/c9b8E56DdLmR7TwVquZ/TTvInh9PHG3fYjThJJQZ4AEjH/UEVsj2Hu+Nfx/mO/daN5S7cJFMmx8qymlFl+XO/FMdwdqjt7hct74NR4JKQss/jJo6QYetxfr0MIuQg6o77ioKTQz+cg1NBB5gzvrzuF/J43o4SMv1W9/XLwF9CxrWl8OTQ3x0XX8T2uBXYFPK+ljvNOe77+RSQHIDLcH6Rty1mW48PeX85sDmS9aCY5FDAsn4BfBjyXnETsPt+GvCw+/qzsLiGUXhyuBFIL2SZfwE+DYl3ccg4H84ZzUCcs6F9wKVAbNg81uL+YnbfN8NJMnlJT4E2IePrut+Pc933jwMvu6+vBRaGzf8F4Hc4B+cc4IKQcX+i5Mkh3h3eH+dH046w8Y8Ar4Rsk7kh4zoAJ93XRW2T0599QXG4y93JD2cvywk7swv7fufi7O+hf0EKTg6XABtwfnj4wub1Kmcmh4E4B29fyLC33fjztne7kHEFnTlcUlDcIWUO4Z6BuvP9NGTcVTjHMn/Id0NxqgLP+lhZXdscfqSqCap6rqreq6onQ8btDHl9LvCA2wB7WEQO4/zabO7+7VL3U3FtL2R5LYHtGln9eSLOAXpFyDJnucNxlxsaY2HLBKeaApwDTlHS816o6gn3ZaGN9CLSF+eX/jWquqGYeYfH2tx9XZL1KCiGtm4DabqIHMU50DUOK5Ye8voEP6xTSZa9H2hcSJ1wM3d8ntPzVNUgkIbzy3gTTvKaCOwTkakikrcdzsWpx8/7rNfiVDM2LWS+mcB/gevcQdcBb4bMq0/Y9/UG4Byc709MCda7MC3c/wfd5TUPW95vwmIP/wziRSSmmG1SJFVdgpMgLxaRC3ASzfQiJlns7u+n/3Cqqwqa92c4VZTPAXtFZLKI1Ctkvs2Bne5nnWc7zjYqaHuHvi5wmIg8ICJr3av0DuO0cYV+r/eGvD4J7FfVQMh7KGLfLYnqmhyKEnqw3wk8HvbFqqWqb+P8KmwhIhJSvlUh89wJtCrkAKNh7/fjfMgdQ5ZZX51GQdzltoxgmeCc5ewEflxEmRIRkW44O+LtqjovgknCY93tvi7JehTkeZyqoRRVrYdzUJKiJzmtJMv+GsjCaV85TURqAyOB0G3QMmS8D6eqaDeAqr6lqgNwDqiK08ANzuczMuw7Fq+qu0LmG/4deRsYJyL9cKqaPg+Z1xdh86qjqvfgnLnmlmC9C3M1zi/+vO/W1rDl1VXVyyOZURHb5IxihUz+Gs5Z3U3Ae6p6qqQrUkRcT6tqD6Ajzhntg4XEshto6X7WeVrhVNflbe+kkHGh2/704vJeiMhAnHaYsUADN4kdIfLvdZmy5FC0KcB4Eekjjtru5Zx1cQ4aucD/iEiMiIzBaSQtyFKcA9Jf3HnEi0h/d9xeICnvqhf3V8gU4CkRaQIgIi1EZLhbfhpwq4h0EJFaOFUGBXLPau4H/k9EbhOReiLiE5EBIjK5pBtDRDrhnMX8TFU/jnCy+0QkSUQa4hzA3ynpehSiLk6j3TH31+M9JZh2Gs7nliQiDXAahAukqkdwGsifEZERIhIrIq1xGgbTgH+HFO8hImPcHwG/wEkqi0WknYhcIiI1cOqMT+KcHYBzscHjInIugIgkisjoYuKfgXNA/T3wTsgv10+AtiJykxtnrIj0EpH27q/LD4CJIlJLRDrgtKFFRESaisgEnM/pEXeZS4GjIvKQOPel+EWkk4j0imB+RW2TUBk4VUBtwob/GydR3YhT914m3O3VR0Ricc5O8i4YAWdfDY0j7wzm1+62HoxT1TO1gO19AU57UVHq4hxTMoAYEXkU54IIT1hyKIKqLsdp5HkWp+5vE079JKqajfNr8lZ33LU4X4aC5hPA+dKcj3M6m+aWB6f+OxVIF5G8KoqH3GUtdqtM5uI04KGqM3Eagj9zy3xWzDq85y7rdpxfOntx6j4/imwrnOEBnNPll+SH699Ti5nmLWAOTkPlFnfZJV6PAvwKuB6nAXAKPySdSEwBZuNcJPANhXxueVT1bziJ7e84CWkJzq/moaqaFVL0I5xtfQjnF+0YVc3BaSz+C85ZYTrO1TW/caf5J86Z2BwRycRpnO5TTDxZbsyX4mzfvOGZOO0n1+F81uk4v8ZruEUm4FQ5pOPUn79S1HJch0XkOM6VMZcDP1HVl93l5X2vu+I0jO4HXsSpCilOUdskdF1P4LSrfOlWXfV1h6fhfHYKLIxgeZGqh/P9OIRTRXQA53MHeAno4MbxH/cYMArnDHI/Tvvlzaq6zi0/AWdbpOMks7dxfjAUZjYwE6fNYztOYiqoKqpC5F2RYkyZE5FtOI3Zc72OxVQ9IvIysFtV892nFI1E5K84F3pEfMbmpei78cIYY4rhVu2NwbnXKCq5VUlxOGddvXDuV7nT06BKwKqVjDGVioj8Aef+mSdUdavX8RShLk7133Gcdq4nKV11riesWskYY0w+duZgjDEmnyrR5tC4cWNt3bq112EYY0ylsmLFiv2qmljQuCqRHFq3bs3y5cu9DsMYYyoVESn0LnmrVjLGGJOPJQdjjDH5WHIwxhiTT5VocyhITk4OaWlpnDpVZs/jMua0+Ph4kpKSiI2N9ToUY8pFlU0OaWlp1K1bl9atW3Pmg1ONOTuqyoEDB0hLSyM5OdnrcIwpF1W2WunUqVM0atTIEoMpcyJCo0aN7KzUVGlVNjkAlhhMubHvlqnqqnRyMKYqUVVyA8HiCxpTBiw5lCO/30/Xrl3p1KkTP/nJTzhx4kTxExXi1ltv5b333gPgzjvvZM2aNYWWnT9/Pl999VWJl9G6dWv2799ffMECbNu2jU6dOp1+P27cODp37sxTTz1V4nnl5OTQo0ePfMN/+9vf0rJlS+rUObMXxB07djBkyBC6detG586dmTFjxulxr732GikpKaSkpPDaa69V6HqcrVM5AT5auYuf/ns5vR+fS5vfzCDlf2fS78/zGP/vFXz4bRpZuQX1j2PM2auyDdLRoGbNmqxcuRKAG264gUmTJnH//fefHh8IBPD7/SWe74svvljk+Pnz51OnTh0uuuiiEs+7LKSnp/PVV1+xfXtpuiiGRYsWFRj7VVddxYQJE0hJSTlj+B//+EfGjh3LPffcw5o1a7j88svZtm0bBw8e5LHHHmP58uWICD169GDUqFE0aNCgQtajtIJB5b1v0nhi9noyMrNokVCT/uc3JqlBTQTYeegkX23ez6zUdP4ycx0/H9qW63q1xOezqi5TduzMoYIMHDiQTZs2MX/+fIYMGcL111/PhRdeSCAQ4MEHH6RXr1507tyZF154AXCqECZMmECHDh244oor2Ldv3+l5DR48+PTjQmbNmkX37t3p0qULQ4cOZdu2bUyaNImnnnqKrl27snDhQjIyMvjxj39Mr1696NWrF19++SUABw4cYNiwYXTr1o2f/vSnFPaE3jp16vDAAw/QvXt3hg4dSkZGBgArVqygS5cu9OvXj+eee+50+WHDhrFv377Ty88TCARo06YNqsrhw4fx+XwsWLDgjO2Tt04jR47MF0ffvn1p1qxZvuEiwtGjRwE4cuQIzZs7/dTPnj2byy67jIYNG9KgQQMuu+wyZs2alW/6kq5Hedp/LIubX17Kr99bRVKDmrxxRx8W/noIT13blQeGteP+Ye146tqufP3wUP59R29aNazFbz78nnFTFrPnyMniF2BMhKrFmcNjH6eyZvfRMp1nh+b1+N1VHSMqm5uby8yZMxkxYgQAS5cuZfXq1SQnJzN58mTq16/PsmXLyMrKon///gwbNoxvv/2W9evX8/3337N37146dOjA7bfffsZ8MzIyuOuuu1iwYAHJyckcPHiQhg0bMn78eOrUqcOvfvUrAK6//np++ctfMmDAAHbs2MHw4cNZu3Ytjz32GAMGDODRRx/lv//9L5MnF9yt9PHjx+nevTtPPvkkv//973nsscd49tlnue2223jmmWe4+OKLefDBB0+Xnz59OldeeeXps6Y8fr+ftm3bsmbNGrZu3UqPHj1YuHAhffr0IS0tjfPPPx+Azz//nN/9LvIupSdOnMiwYcN45plnOH78OHPnOh3P7dq1i5Ytf+jTPSkpiV27duWbvqTrUV427cvkxheXcvBENo9f3Ynre7cqtOHb5xMGpiQy4PzGvLs8jcc+TmXUs1/ywk096N4qsjMjY4piZw7l6OTJk3Tt2pWePXvSqlUr7rjjDgB69+59+vr4OXPm8Prrr9O1a1f69OnDgQMH2LhxIwsWLGDcuHH4/X6aN2/OJZdckm/+ixcvZtCgQafn1bBhwwLjmDt3LhMmTKBr166MGjWKo0ePkpmZyYIFC7jxxhsBuOKKKwqtbvH5fFx7rdPl9Y033siiRYs4cuQIhw8f5uKLLwbgpptuimibDBw4kAULFrBgwQIeeeQRFi1axLJly+jVy+mTfvfu3TRs2JBatWpFND+At99+m1tvvZW0tDRmzJjBTTfdRDAYLPBMKPxgW9r1KGurdx1h7AuLyQ0qH9xzETf0OTeiK6JEhLG9WvLhff2pGevnxheXsGTLgQqI2FR11eLMIdJf+GUttM0hVO3atU+/VlWeeeYZhg8ffkaZGTNmFHtwUNWIDiDBYJCvv/6amjVr5htXmksyRSTiZYcbOHAgkyZNYvfu3fz+97/niSeeYP78+QwaNAiAmTNn5tsWxXnppZdOVxf169ePU6dOsX//fpKSkpg/f/7pcmlpaQwePPiMaUu7HmVp2/7j3PTSEmrFxfDGnX1Ibly7+InCtG1al/fu6ce4yYu57dVlvHpbb3onF/xjwZhI2JmDx4YPH87zzz9PTk4OABs2bOD48eMMGjSIqVOnEggE2LNnD59//nm+afv168cXX3zB1q1OT4kHDx4EoG7dumRmZp4uN2zYMJ599tnT7/MS1qBBg3jzzTcB56B86NChAmMMBoOnr5R66623GDBgAAkJCdSvX59FixYBnJ5Pcfr06cNXX32Fz+cjPj6erl278sILLzBw4ECg8PaGorRq1Yp58+YBsHbtWk6dOkViYiLDhw9nzpw5HDp0iEOHDjFnzpx8iae061FWDh7P5tZXliIivFnKxJCnSd143r67L+fUj+eu15ezJeNYGUZqqhtLDh6788476dChA927d6dTp0789Kc/JTc3l6uvvpqUlBQuvPBC7rnnntPVHqESExOZPHkyY8aMoUuXLqerfq666io+/PDD0w2pTz/9NMuXL6dz58506NCBSZMmAfC73/2OBQsW0L17d+bMmUOrVq0KjLF27dqkpqbSo0cPPvvsMx599FEAXnnlFe677z769etX4FlJQWrUqEHLli3p27cv4JxJZGZmnm6c37hxIxdccEGB0/76178mKSmJEydOkJSUxMSJEwF48sknmTJlCl26dGHcuHG8+uqriAgNGzbk//7v/043xD/66KMFVr2VZj3KQiCoTHjrG3YfOcWUm3vQ+iwSQ54mdeN59dbe+H3CHa8t5/CJ7DKI1FRHVaIP6Z49e2p4Zz9r166lffv2HkVUtdSpU4djx8r/V+iiRYt44403TievaHe237FnP9vI3+ds4K8/vpBrexWcmEtr+baDjJuymEEpibx4S0/Pq85MdBKRFaras6BxEZ05iMgIEVkvIptE5OECxouIPO2OXyUi3d3hLUXkcxFZKyKpIvLzkGkaisinIrLR/d8gZNwj7rzWi0jJKqBNpTVgwIBKkxjO1ortB3lq7kZGdWnO2J4ti5+ghHq2bsgjI9szb90+XvlyW5nP31R9xSYHEfEDzwEjgQ7AOBHpEFZsJJDi/t0NPO8OzwUeUNX2QF/gvpBpHwbmqWoKMM99jzv+OqAjMAL4lxuD8UhFnDVUJ1m5AR58bxXn1Ivn8as7lduv+tv6t+bS9k3488y1ZX4pt6n6Ijlz6A1sUtUtqpoNTAVGh5UZDbyujsVAgog0U9U9qvoNgKpmAmuBFiHT5D3P4DXgRyHDp6pqlqpuBTa5MZRYVagyM9HpbL5bz32+mS0Zx/nTmAupG19+/UGICE9c04X6NWN56P1V9lwmUyKRJIcWwM6Q92n8cICPuIyItAa6AUvcQU1VdQ+A+79JCZaHiNwtIstFZHneHbuh4uPjOXDggCUIU+by+nOIj48v8bQb9mby/PxN/Khrcy5um1gO0Z2pQe04HhvVie93HeGlRVvLfXmm6ojkPoeCznnDj7hFlhGROsD7wC9Utbjz20iWh6pOBiaD0yAdPj4pKYm0tDQKShzGnK28nuBKQlX53Uep1K4Rw/9dGV4zW34uv/AchnVoyv/7dAPDOp5zVpfLmuojkuSQBoS2mCUBuyMtIyKxOInhTVX9IKTM3ryqJxFpBuwrbl4lERsba710magyb+0+vt5ygMdGdaRRnRoVtlwR4Q8/6sTQJ7/gj5+s4aVbe1XYsk3lFUm10jIgRUSSRSQOp7F4eliZ6cDN7lVLfYEj7kFfgJeAtar6/wqY5hb39S3ARyHDrxORGiKSjNPIvbTEa2ZMFMkJBPnTzLW0SazN9X3K9rLVSDStF8/PLjmfeev2MX/9vuInMNVesclBVXOBCcBsnAblaaqaKiLjRWS8W2wGsAWn8XgKcK87vD9wE3CJiKx0/y53x/0FuExENgKXue9R1VRgGrAGmAXcp6r20HpTqU1duoMtGcd5ZGR7Yv3e3Ht6W/9kkhvX5g+frCHHGqdNMarsTXDGRItTOQEG/u1z2jSuzdS7+3p6Q9q8tXu547Xl/O6qDtzW36pdq7uiboKrFg/eM8ZLby3ZQUZmFs+O61axiWFi/XyDLlG4yPcbnv34CGPn9KG2ZJXDco+U/TxNhbNnKxlTjk7lBJj0xWb6tmlInzaNvA4HEfhVzDQOUJ9XAiO8DsdEMUsOxpSjqUt3sC8zi58Pbet1KKd1923iUt9yXsi9ksNql7WagllyMKacnMoJ8PwXm+md3JB+53l/1hDqgZh3OUZNXsi90utQTJSy5GBMOflo5S72Hs3i50NTvA4ln/a+nVzl+5pXA8M5qHW9DsdEIUsOxpQDVeXFhVvp0KweF0XZWUOen8X8h5PE82quPfjY5GfJwZhy8MWGDDbuO8Zdg5Kjti+FFN8uhvuW8mpgOJlacZ0cmcrBkoMx5eDFhVtpWq8GV1zY3OtQinRvzHSOUps3A0O9DsVEGUsOxpSxtXuOsmjTfm65qDVxMdG9i3XxbWGgbxUv5l7OKS2/x4ebyie6v7nGVEIvL9pKzVg/N/Q+1+tQInKv/yP2k8C0wGCvQzFRxJKDMWXoyIkcPl61m6u7t6B+rcrxS7yvby3dZCMvBi4noNHZPmIqniUHY8rQ+9+kcSonyA0ePHm1tETgjpgZ7NCmzAt29zocEyUsORhTRlSVN5dsp1urBDo2z/9co2g2wreM5uznZXukhnFZcjCmjCzecpDNGce5oU/laGsIFSNBbomZw+JgR1KDlS9+U/YsORhTRt5csp36NWO5snMzr0Mplev8n1OTU/ZAPgNEmBxEZISIrBeRTSLycAHjRUSedsevEpHuIeNeFpF9IrI6bJp3QjoA2iYiK93hrUXkZMi4SWe5jsaUu4zMLGanpnNNjyTiY/1eh1Mq9eU4P/EvYHrgIvZp5aoWM2Wv2OQgIn7gOWAk0AEYJyLhvaOPxOnOMwW4G3g+ZNyrQL6fIqp6rap2VdWuOH1Mh/YvvTlvnKqOD5/WmGjz/jdp5ATUky5Ay9Kt/llkE8sbuZd6HYrxWCRnDr2BTaq6RVWzganA6LAyo4HX1bEYSBCRZgCqugA4WNjM3X6mxwJvl2YFjPGaqvLeijR6ntuA8xLreB3OWWnjS2eI71veDlxCjlbOMyBTNiJJDi2AnSHv09xhJS1TmIHAXlXdGDIsWUS+FZEvRGRgQROJyN0islxElmdkZES4KGPK3ndpR9i07xjX9EjyOpQycaN/Lhk04NNgD69DMR6KJDkUdFdMeMfTkZQpzDjOPGvYA7RS1W7A/cBbIlIv38xVJ6tqT1XtmZiYGOGijCl7763YSXysjysqaUN0uMG+lbQggzcDVrVUnUWSHNKAliHvk4DdpSiTj4jEAGOAd/KGqWqWqh5wX68ANgPR042WMSFO5QSYvnI3Izs1o2585bgjujh+UcbFfMaXwU5sCZ7jdTjGI5Ekh2VAiogki0gccB0wPazMdOBm96qlvsARVd0TwbwvBdapalreABFJdBvBEZE2OI3cWyKYlzEV7tM1ezl6KrfKVCnlGeufTwy5vGVPa622ik0OqpoLTABmA2uBaaqaKiLjRSTvSqIZOAfwTcAU4N686UXkbeBroJ2IpInIHSGzv478DdGDgFUi8h3wHjBeVQtt0DbGS++tSKN5/Xj6tYnODn1Kq4kcYbhvOe8FBtnTWqupmEgKqeoMnAQQOmxSyGsF7itk2nFFzPfWAoa9j3NpqzFRLf3IKRZuzOC+Iefj81W9B9bd4J/Lf4N9+W+wLz/2L/Q6HFPB7A5pY0rpPyt3EVT4cfeqVaWUp59vDW1kN2/mWtVSdWTJwZhSmr5yN11bJtC6cW2vQykXInC9fx7faFvWBVsWP4GpUiw5GFMKm/ZlsmbPUUZ1ie5uQM/WGP8iYsllWuBir0MxFcySgzGlMH3lbnxCpX3IXqQaSiaX+VbwYWAg2XbHdLViycGYElJVpn+3m37nNaJJvXivwyl3Y/2fc4i6zLU7pqsVSw7GlNCqtCNsO3CC0V0ifUJM5TbQ9z3NOGB9TFczlhyMKaHp3+0mzu9jeKfqcfewX5Rr/AtYEOzMHm3odTimglhyMKYEAkHlk1W7ubhdIvVrVp+bw37i/4IgPt4PFPgcTFMFWXIwpgSWbD3A3qNZjO5ata9SCtfKt49+vlSmBQYT1Kp3w5/Jz5KDMSXw8Xe7qR3nZ+gFTb0OpcKN9c9nhzZlSfACr0MxFcCSgzERys4NMuP7dC7r0JSacdXvss6RvqXU5bg1TFcTlhyMidDiLQc4cjKHKztXryqlPPGSwyj/18wI9uGo1vQ6HFPOLDkYE6FZqenUivMzIKWx16F45hr/F2QRx8xAH69DMeXMkoMxEQgElTmpexnSrgnxsdWvSilPV9lMG9ltVy1VA5YcjInAtzsOsf9YVrW5t6EwIvBj/0KWant2Bq173qrMkoMxEZidmk6c38eQdnZA/JF/EQAfBgd4HIkpTxElBxEZISLrRWSTiDxcwHgRkafd8atEpHvIuJdFZJ+IrA6bZqKI7BKRle7f5SHjHnHntV5Ehp/NChpztlSVWanpXHR+oyrTT/TZaCEH6OdL5YPAQFS9jsaUl2KTg9uf83PASKADME5EOoQVG4nT13MKcDfwfMi4V4ERhcz+KVXt6v7NcJfXAaf70I7udP/K61PaGC+s3ZPJzoMnGdGxelcphRrjW8g2PYdvNMXrUEw5ieTMoTewSVW3qGo2MBUYHVZmNPC6OhYDCSLSDEBVFwAl6QN6NDBVVbNUdStOv9S9SzC9MWVqdmo6PoFLO1S/G98KM9K/lHiy+MAapqusSJJDC2BnyPs0d1hJyxRkglsN9bKINCjJvETkbhFZLiLLMzIyIliUMaUzOzWdnq0b0rhODa9DiRp15BQjfMv4ONCPLI2oK3pTyUSSHAp6kEp4TWMkZcI9D5wHdAX2AE+WZF6qOllVe6pqz8REayQ05WPb/uOsS89kuFUp5TPGv5Cj1OazYDevQzHlIJLkkAaEdiCbBOwuRZkzqOpeVQ2oahCYwg9VRyWelzHlZXZqOgDDO1qVUrj+vtU04ZDd81BFRZIclgEpIpIsInE4jcXTw8pMB252r1rqCxxR1T1FzTSvTcJ1NZB3NdN04DoRqSEiyTiN3EsjiNOYMjc7NZ1OLeqR1KCW16FEHb8oV/sXMT/YlQNa1+twTBkrNjmoai4wAZgNrAWmqWqqiIwXkfFusRnAFpzG4ynAvXnTi8jbwNdAOxFJE5E73FF/E5HvRWQVMAT4pbu8VGAasAaYBdynqoGzX1VjSmbv0VN8s+MwwztYlVJhxvgXkksMHwf6eR2KKWMRtSS5l5nOCBs2KeS1AvcVMu24QobfVMTyHgcejyQ2Y8rLnDV7ARhRze+KLko7XxodZSsfBAZya8wcr8MxZcjukDamEHNS02nTuDbnN6njdShRbYx/Iav0PDYGq0ef2tWFJQdjCnDkRA5fbz7A8E7nIGI9nxVllP8r/ASsYbqKseRgTAHmrdtLblDtEtYIJMpRBvlW8VGgv3UhWoVYcjCmALNWp3NOvXg6t6jvdSiVwtX+ReyhEYuD7b0OxZQRSw7GhDmRncuCjRkM79gUn89+CUdimG85dTnBB0GrWqoqLDkYE2bBhgxO5QStSqkE4iWHkf4lzAz05mS2XXleFVhyMCbM7NS9JNSKpXdyQ69DqVSu9i3iODWZsybd61BMGbDkYEyI7Nwgc9fu5dL2TYnx2+5REn1862hBBh98s8vrUEwZsG+/MSEWbzlA5qlc67uhFHyi/Mj/JQs3ZrDv6CmvwzFnyZKDMSFmp6ZTK87PgJTGXodSKV3tX0RQYfp39qzMys6SgzGuYFCZs2Yvg9slEh9rnQ+Wxvm+3XRJqm9VS1WAJQdjXN/uPERGZpZdpXSWru7WgjV7jrIu/ajXoZizYMnBGNes1enE+oUhFzTxOpRK7aouzYnxCR/a2UOlZsnBGEBVmZ26l/7nN6ZefKzX4VRqjerUYHC7RP6zcheBYHEdQppoZcnBGGDtnkx2HDxhVUpl5OpuSew9msXXmw94HYoppYiSg4iMEJH1IrJJRB4uYLyIyNPu+FUi0j1k3Msisk9EVodN84SIrHPLfygiCe7w1iJyUkRWun+TMKaczU5NRwQu62DdgZaFoe2bUDc+hg++SfM6FFNKxSYHEfEDzwEjgQ7AOBHpEFZsJE53ninA3cDzIeNeBUYUMOtPgU6q2hnYADwSMm6zqnZ1/8YXMK0xZWp2ajq9zm1I4zo1vA6lSoiP9XPFhc2YlZrOiexcr8MxpRDJmUNvYJOqblHVbGAqMDqszGjgdXUsBhLy+ohW1QXAwfCZquoctwtSgMVAUmlXwpizsf3AcdalZzKso501lKUx3ZM4kR1gdqo9TqMyiiQ5tAB2hrxPc4eVtExRbgdmhrxPFpFvReQLESnwMY8icreILBeR5RkZGSVYlDFnyjt4WXtD2ep5bgOSGtS0ex4qqUiSQ0HPLA6/BCGSMgXPXOS3QC7wpjtoD9BKVbsB9wNviUi9fDNXnayqPVW1Z2JiYiSLMqZAs1an07F5PVo2rOV1KFWKzydc3a0FX27az157nEalE0lySANahrxPAsLvjY+kTD4icgtwJXCDqiqAqmap6gH39QpgM9A2gjiNKbF9R0/xzY7D9iylcnJ1txYEFT5aaWcPlU0kyWEZkCIiySISB1wHTA8rMx242b1qqS9wRFX3FDVTERkBPASMUtUTIcMT3UZwRKQNTiP3lojXyJgSmL1mLwDDO1lyKA9tEuvQtWWCVS1VQsUmB7fReAIwG1gLTFPVVBEZLyJ5VxLNwDmAbwKmAPfmTS8ibwNfA+1EJE1E7nBHPQvUBT4Nu2R1ELBKRL4D3gPGq2q+Bm1jysKc1HTaNK5NSpM6XodSZY3p3oJ16Zms2W2P06hMYiIppKozcBJA6LBJIa8VuK+QaccVMvz8Qoa/D7wfSVzGnI0jJ3L4evMB7hzYBhHrDrS8XNm5OX/4ZA0ffptGh+bhV8GbaBVRcjCmKpq3bi+5QWV4eV/COrF++c4/yjWsHcfgdk34z8rdPDTiAutEqZKwT8lUW7NT0zmnXjxdkhK8DqXKG9OtBRmZWXxpj9OoNCw5mGrpZHaALzZkMKxjU3w+q1Iqb5e0b0K9+Bg+tMdpVBqWHEy19MWGDE7lBO0S1gpSI8bPlV2aMzt1L8ez7HEalYElB1MtzU5NJ6FWLL2TG3odSrUxplsLTuYEmLXaHqdRGVhyMNVOdm6QeWv3MvSCptY4WoF6nNuAVg1r8cG3VrVUGdieYaqdxVsOcPRULiPsxrcKJeI8TuOrzQfYc+Sk1+GYYlhyMNXO7NR0asX5GZjS2OtQqp2ru7VAFT5aWezTdYzHLDmYaiUYVOas2cvgdonEx/q9Dqfaad24Nt1bJfDBN2m4j1MzUcqSg6lWvt15iIzMLHs8t4fGdE9iw95jpNrjNKKaJQdTrcxO3UusXxhyQROvQ6m2ruzcjDi/jw+/tYfxRTNLDqbaUFVmrt7DRec1pl58rNfhVFsJteIYckEiH63cTW4g6HU4phCWHEy1sWbPUXYePMlIu0rJc2O6J7H/WBYLN+33OhRTCEsOptqYvTodn8BlHayvaK8NadeEhFqxfGj9PEQtSw6m2pi5Op3eyQ1pVKeG16FUe3ExPq7s3Iw5a9LJPJXjdTimAPbIblMtbNp3jI37jnHDwVdh4hyvwzE4VUtvLN7BzNXpjO3ZsvgJTIWK6MxBREaIyHoR2SQiDxcwXkTkaXf8KhHpHjLuZRHZJyKrw6ZpKCKfishG93+DkHGPuPNaLyLDz2YFjQHnxjeA4f7lHkdi8nRrmUCbxrV5b7k9TiMaFZsc3P6cnwNGAh2AcSIS3p3TSJy+nlOAu4HnQ8a9CowoYNYPA/NUNQWY577Hnfd1QEd3un/l9SltTGnNWp1OV9lIM7EeZ6OFiDC2V0uWbjvI5oxjXodjwkRy5tAb2KSqW1Q1G5gKjA4rMxp4XR2LgQQRaQagqguAgvbI0cBr7uvXgB+FDJ+qqlmquhWnX+reJVgnY86QdugE3+86wkj/Mq9DMWHGdG+B3ydMW77T61BMmEiSQwsg9JNLc4eVtEy4pqq6B8D9n3dXUkTzEpG7RWS5iCzPyMgodiVM9ZX3iOgRvqUeR2LCNakbz9ALmvD+ijRy7J6HqBJJciiom6zwh6JEUiZSEc1LVSerak9V7ZmYmFjKRZnqYHZqOu2b1eNc3z6vQzEFuLZXS/Yfy2beWvt8okkkySENCL2UIAkIf6RiJGXC7c2renL/530zSjMvYwq0L/MUy7cfsh7fotjFbRNpWq+GVS1FmUiSwzIgRUSSRSQOp7F4eliZ6cDN7lVLfYEjeVVGRZgO3OK+vgX4KGT4dSJSQ0SScRq5rT7AlMqc1L2owsgLLTlEqxi/j2t6JDF//T7Sj5zyOhzjKjY5qGouMAGYDawFpqlqqoiMF5HxbrEZwBacxuMpwL1504vI28DXQDsRSRORO9xRfwEuE5GNwGXue1Q1FZgGrAFmAfepauCs19RUS7NT02nTuDYpTep4HYopwtieLQkqvLfCzh6iRUQ3wanqDJwEEDpsUshrBe4rZNpxhQw/AAwtZNzjwOORxGZMYQ6fyObrzQe4a1AbRApqyjLR4txGtenXphHvLN/JvYPPx+ezz8tr9vgMU2V9umYvuUG1B+1VEtf1bsnOgydZvOWA16EYLDmYKmzG93tokVCTC1vU9zoUE4HhHc+hXnwMU5dZ1VI0sORgqqQjJ3JYuHE/V3ZuZlVKlUR8rJ+ru7VgVmo6h09kex1OtWfJwVRJs1PTyQ0qV3Ru5nUopgTG9mpJdm6QD+xR3p6zp7Kaijex/Kt5Psl+iFZyDhdOObfg2ypNVOrYvD5dWybw5pLt3Na/tZ31ecjOHEyVc1Dr8mWwE1f4FmPHlsrnxr7nsjnjOF9bw7SnLDmYKmd2oCcB/FzpX+x1KKYUruzcjIRasbyxeLvXoVRrlhxMlfNJsB/JsocOYgeXyig+1s/Yni2ZnbqXvUftjmmvWHIwVcp+rcfXwQ5WpVTJXd+7FYGgMnWpXdbqFUsOpkqZFehFEJ9VKVVyrRvXZlDbRN5euoNce5S3Jyw5mCrlk2BfzpNdtBP7xVnZ3dinFelHTzHXHuXtCUsOpsrYp/VZEmxvVUpVxCUXNKF5/XhrmPaIJQdTZcwI9EGtSqnKiPH7uL5PKxZt2s8W62O6wllyMFXGfwL9uUC209Znd9dWFWN7tSTGJ7yxeIfXoVQ7lhxMlbA1eA4rNYWr/V96HYopQ03qxnP5hc14d/lOjmXleh1OtRJRchCRESKyXkQ2icjDBYwXEXnaHb9KRLoXN62IvCMiK92/bSKy0h3eWkROhoybFL48Y8J9FLwIIcgo/1deh2LK2O0DksnMyuVd60a0QhX7bCUR8QPP4fTWlgYsE5HpqrompNhInO48U4A+wPNAn6KmVdVrQ5bxJHAkZH6bVbXrWa2ZqTZU4T+BAfTzraGZHPQ6HFPGurZMoMe5DXjly23c3K81fusIqEJEcubQG9ikqltUNRuYCowOKzMaeF0di4EEEWkWybTiPFlrLPD2Wa6LqaZW6nls03P4kc+qlKqq2/sns+PgCeat3et1KNVGJMmhBRB6PpfmDoukTCTTDgT2qurGkGHJIvKtiHwhIgMLCkpE7haR5SKyPCMjI4LVMFXVR4H+xJHNCP9Sr0Mx5WR4x6a0SKjJS4u2eh1KtRFJcijoHE4jLBPJtOM486xhD9BKVbsB9wNviUi9fDNRnayqPVW1Z2JiYqHBm6otR/18HOjHZb5vqCcnvQ7HlJMYv49bLjqXJVsPsnrXkeInMGctkuSQBrQMeZ8E7I6wTJHTikgMMAZ4J2+Yqmap6gH39QpgM9A2gjhNNbQoeCEHqM+P/Iu8DsWUs2t7taJWnJ+Xv7Szh4oQSXJYBqSISLKIxAHXAdPDykwHbnavWuoLHFHVPRFMeymwTlXT8gaISKLbkI2ItMFp5N5SyvUzVdyHgf4kkMnFvu+8DsWUs/o1Y/lJjyQ+/m43++xpreWu2OSgqrnABGA2sBaYpqqpIjJeRMa7xWbgHMA3AVOAe4uaNmT215G/IXoQsEpEvgPeA8arql2CYvI5pvF8GuzBFf4lxEnA63BMBbitfzK5QeXVr7Z5HUqVF1E3oao6AycBhA6bFPJagfsinTZk3K0FDHsfeD+SuEz1NiPQh5PEM8a/0OtQTAVp3bg2l3dqxr+/3s74wedRLz7W65CqLLtD2lRa0wKDOU920V02Fl/YVBn3DD6PzKxc3rRHapQrSw6mUtocbMZybcdY/xf2BNZqplOL+gxMacxLi7ZyKseqE8uLJQdTKb0buBg/Aa62KqVq6Z7B57H/WBbvrUgrvrApFUsOptLJVR/vBwYyxPctTcSuea+O+rVpRNeWCUxesMV6iisnlhxMpfNFsAsZNGCs/wuvQzEeERHuGXweOw6e4L/f7/E6nCrJkoOpdKYFLqYxRxjiW+l1KMZDl7VvyvlN6vD8/M04F0yasmTJwVQq+7Ue84LdGeNfSKzd21Ct+XzCvYPPY116JrNT7YF8Zc2Sg6lUPggMJJcYfmJVSgYY1aU5bRrX5h9zNxAM2tlDWbLkYCqNoApvBobSS9aRYl2BGpwH8v3P0BTWpWcyKzXd63CqlIjukDYmGiwMXsh2PYf7Y9/1OhRTlIn1K3RxV6nwjPyNf8ytw4iO5+Cr6M6AKnh9f1hu+V6pZ2cOptJ4I3ApjTjCCN8yr0MxUcQvys9jPmDD3mN25VIZsuRgKoVd2oh5we5c659PDbGO5s2ZrvAtpm3TOvxj7gYC1vZQJiw5mEphau4lKDDOP8/rUEwU8ovy86Ft2ZxxnI9WWntUWbDkYKJetvp5OzCES3wraenb73U4JkqN7HQOnVrU48k5G+yZS2XAkoOJenOCPdlPAjf653odioliPp/wm5Ht2XX4JK9/vc3rcCo9Sw4m6r2cO5JWspdB1tubKcZF5zdmcLtEnv1sE4dPZHsdTqUWUXIQkREisl5ENonIwwWMFxF52h2/SkS6FzetiEwUkV0istL9uzxk3CNu+fUiMvxsV9JUXiuCKXyjbbnDPxO/WEOjKd7DIy8gMyuX5z7f5HUolVqxycHtz/k5YCTQARgnIh3Cio3E6es5BbgbeD7CaZ9S1a7u3wx3mg443Yd2BEYA/8rrU9pUPy/lXk49jnON3RFtInTBOfX4cfckXvtqOzsPnvA6nEorkjOH3sAmVd2iqtnAVGB0WJnRwOvqWAwkiEizCKcNNxqYqqpZqroVp1/q3iVYJ1NF7AwmMivYi+v986gtWV6HYyqR+y9riwg8MXu916FUWpEkhxbAzpD3ae6wSMoUN+0EtxrqZRFpUILlISJ3i8hyEVmekZERwWqYyuaVwHB8BLklZo7XoZhKpnlCTX46qA3Tv9vNki0HvA6nUookORR0L3p45W9hZYqa9nngPKArsAd4sgTLQ1Unq2pPVe2ZmJhYwCSmMjuqNXknMISrfF/TTA56HY6phO4ZfD4tEmryu+mp1iFQKUSSHNKAliHvk4DdEZYpdFpV3auqAVUNAlP4oeookuWZKu7twFCOU5M7YmZ6HYqppGrG+fnfK9qzLj2TN5fs8DqcSieS5LAMSBGRZBGJw2ksnh5WZjpws3vVUl/giKruKWpat00iz9XA6pB5XSciNUQkGaeRe2kp189UQqc0lhdzR9Lft5pOvm1eh2MqsRGdzmHA+Y15cs569h+zdquSKDY5qGouMAGYDawFpqlqqoiMF5HxbrEZwBacxuMpwL1FTetO8zcR+V5EVgFDgF+606QC04A1wCzgPlW12x2rkXcCQ8igAT+L+cDrUEwlJyJMHNWBE9kBnphljdMlEdEju93LTGeEDZsU8lqB+yKd1h1+UxHLexx4PJLYzFnw6lHDRchWP5Nyr6KXrKOPrPM6HFMFnN+kLrcPSGbygi1c0zOJXq0beh1SpWB3SJuo8kFgIHtoxISY/yAV/Fh+U3X94tIUkhrU5KH3V9lzlyJkycFEjVz18a/AaDrLZgb5VnkdjqlCasXF8OcxF7Il4zjPfLbR63AqBUsOJmp8EBjIDm3KfTEf2VmDKXMDUxK5pkcSL3yxhdTd5duLWlVgycFEhSyN4Z+5Y+gsmxnmW+51OKaK+t8r2pNQK46H3l9Fjt37UCRLDiYqvBm4lF0k8uuYd+yswZSbhFpx/GF0R1bvOsozn9mD+YpiycF47pjG81zuaPr7VjPAv7r4CYw5CyMvbMaY7i149rONrNh+yOtwopYlB+O5lwIjOUB9Hox5x+tQTDXx2KiONE+oyS/fWcmxLOuTvCCWHIynMrQeU3KvYIRvKV19m70Ox1QTdeNjeerarqQdOsFj01OLn6AasuRgPPVE7rVkEcevY6Z6HYqpZnq1bsi9g8/n3RVp/OfbXV6HE3UsORjPrAom827gYm73z6SNL93rcEw19ItLU+id3JBHPvieDXszvQ4nqlhyMJ5QhYk5t9CIo0yI+Y/X4ZhqKsbv49lx3ahdI4bx/15B5qkcr0OKGpYcjCf+E+zPN9qWh2KmUldOeh2Oqcaa1Ivnueu7sf3gCR56fxXOo+KMJQdT4Q5rbR7PuYEusokf+xd6HY4x9GnTiIdGtGPG9+k897nd/wARPpXVmLL0p9zrOURdXov9Kz6xX2kmOtw1sA3r9mTy9zkbaN24Nld2bu51SJ6y5BDKq0dYT6w+z3n5MtCRaYEh3OP/iI6+7V6HY6qSs9x/BfizxrBTfsP9b2XT/L0/0t1Xfc8iIqpWEpERIrJeRDaJyMMFjBcRedodv0pEuhc3rYg8ISLr3PIfikiCO7y1iJwUkZXu36Tw5ZnK6YTW4JHcO0mWPfzcOvIxUaiG5PJC3FM0k4Pcnf0A24NNvA7JM8UmBxHxA88BI4EOwDgR6RBWbCROd54pwN3A8xFM+ynQSVU7AxuAR0Lmt1lVu7p/4zFVwuO5N7BDm/KnmBeJF7sqxESnhpLJS7F/J4CPG3N+Q7o28DokT0Ry5tAb2KSqW1Q1G5gKjA4rMxp4XR2LgQS3j+hCp1XVOW43ogCLgaQyWB8TpeYGuvNm4FLu8n9CP/9ar8Mxpkjn+3bzWtxfOKR1uCn7EQ5qXa9DqnCRJIcWwM6Q92nusEjKRDItwO3AzJD3ySLyrYh8ISIDCwpKRO4WkeUisjwjIyOC1TBeydB6PJRzF+1lG7+KmeZ1OMZEpLNvKy/G/Z0d2oRbs3/NUa3pdUgVKpLkUNADlMMvMSmsTLHTishvgVzgTXfQHqCVqnYD7gfeEpF6+WaiOllVe6pqz8TExGJWwXgloML9OfdyjJo8HfscNcQecmYqj76+dfwr9p+s1XO5Ifu3HNI6XodUYSJJDmlAy5D3ScDuCMsUOa2I3AJcCdyg7p0nqpqlqgfc1yuAzUDbSFbGRJ8nc8eyMNiZx2JeJcVnz68xlc9Q/7e8EPv/WK9JjMv+XzI032/VKimS5LAMSBGRZBGJA64DpoeVmQ7c7F611Bc4oqp7ippWREYADwGjVPVE3oxEJNFtyEZE2uA0cm85q7U0npgd6Mm/AqO5zv8Z18XM9zocY0rtEv9KXol9gu3ahGuzH2WXNvI6pHJXbHJwG40nALOBtcA0VU0VkfEikncl0QycA/gmYApwb1HTutM8C9QFPg27ZHUQsEpEvgPeA8ar6sGzX1VTkdYFW/JAzng6y2YmxrzmdTjGnLX+/lRej/srGZrA1Vm/5/tgstchlSupCs8R6dmzpy5fXgb9Dle3m+DKaX3TtQFXZ/2eIMKHNR6luVhuN1XHhmALbsv+NQepy9Oxz3KZ/xtvAimD44aIrFDVngWNs2crmTKVqTW5LftBMqnJK3F/s8Rgqpy2vl18WONR2sou7s65n+dyRxPUqtfxuSUHU2ZOaA3uyP4VG7Ql/4r9Jx18O7wOyZhy0USOMDXuD1zhW8ITuddyV84DHNHaXodVpiw5mDJxUuO4I+dXLNd2/CP2OQb5v/c6JGPKVU3J5pnYZ3gs5lUWBDtzRfbjfBM83+uwyowlB3PWjmsN7sz5FUuC7Xkq9l9c5V/sdUjGVAgRuCVmDtPiHkNVuCZ7Ik/kjCVb/V6HdtbsqazRwKuG8DKwX+txe/aDpGpr/h47idH+r7wOyZgK1823mVk1HuYPuTfyXOBHzAt24++xL9DJt83r0ErNzhxMqW0NnsOPsyeyQZOYEvskY/yLvA7JGM/UlZP8LXYKL8b+nQNaj6uy/8j/5tzG4UraFmHJwZTK3EB3RmX/kaNam7fiHucS/0qvQzImKlzq/4a5NR7kFv8c3goM5ZKsJ3k7dwi5WrkOt5UrWuO5HPXzZM5PuDPnV7SWdD6u8dtq3SGKMQWpLyeYGPs6n8T9hjayh0dy72JY9t/4ONC30lz2asnBRGxjsAVjsh/jmcDV/MQ/n3fjHiNJ9nsdljFRq4NvB+/GPcbk2CeJJZef5fwPV2Q/zkeBi8iJ8kZrSw6mWKc0ludyR3NF9uOkaSLPxz7FE7GTrcMeYyIgAsP8K5gR9wj/iH2OLOL4ec4ELs56iim5l0fto8DtaiVTKFWYG+zOH3JvYoc2ZYRvKX+IfZlEOep1aMZUOn5RfuT/klG+r/g82JXJuVfweO6N/L/ca7jct4RrY+bTS9YjUVLrZMnB5KMKXwQ783TuGL7Rtpwvafw79k8M9K/2OjRjKj2fKEP93zLU/y3fB5N5K3AJHwf68X72xbSWdEb5vmKEfyntZYenicIevBeqEt9vUBZOaSz/Dfbl9dxhfKfn0YIMxsd8zHX+z4mVgNfhGVNlndAazAr24t3AxSwJtieIj1aylxG+ZQz2raSHb0P+jrLK+cF7duZQzanC95rM9MBFvBcYxGHq0kZ286eYF7nG/wVxlhSMKXe1JIsx/kWM8S9iv9ZjbqA7M4O9eSUwgsmBK4kni16+9fT3raaPby0dZDs1yjkmO3MIVU3OHE5pLCuCbVkQ7MyMYB92ahNiyGWYbwU3+j+ln29N1NR7GlOdZWpNlgTbsyjYiS+DndioSQDEkUPHVol0bZnAoLaJDGnXpFTztzOHam6fJvB9MJlVwTas0LYsC7YjizhiyKW/bzU/83/IMP9yEuS416EaY0LUlZNc6v+GS90+I/ZpAiuCKXwbTGGl72beXrqDIydzSp0cihJRcnC79Pwn4AdeVNW/hI0Xd/zlwAngVlX9pqhpRaQh8A7QGtgGjFXVQ+64R4A7gADwP6o6+6zWsopThSPUZp8mkKaJbNVz2KbnsFWbsTHYgr00BMBHkLaSxg3+uQzwraa3bx115JTH0RtjItVEDjPSv4yR/mUw/nlyAkGOZ+UWP2EpFJsc3P6cnwMuA9KAZSIyXVXXhBQbidPXcwrQB3ge6FPMtA8D81T1LyLysPv+IRHpgNPXdEegOTBXRNqqatRUfquCIih5/898TSHD816DEETIIYYsjSGHGLKJdf/HkKU/vD5FHJlai0xqcUxrkklNjrrvD2g99pHAPk0gm7gzYqzLcdpIOhf5Uunk20pn31Y6yDZqS1YFby1jTHmJ9ftIqBVXfMFSiOTMoTewSVW3AIjIVGA0EJocRgOvq9OAsVhEEkSkGc5ZQWHTjgYGu9O/BswHHnKHT1XVLGCriGxyY/i69KtZsO/TjjD2ha9R1Dng574KhRzUAdTjewZ9BKnDSepygrpygkaSSS/W08R3mEQ5TBM5THM5QLLsoSGZ1m5gjCm1SJJDC2BnyPs0nLOD4sq0KGbapqq6B0BV94hIXqVZC2Bx2DQtwoMSkbuBu923x0RkfQTrEqoxUJme/VDZ4gWLuSJUtnjBYi4bjxX76y+SmM8tbEQkyaGgCMIvcSqsTCTTlmZ5qOpkYHIx8yp8ISLLC2ulj0aVLV6wmCtCZYsXLOaKcrYxR1JPkga0DHmfBOyOsExR0+51q55w/+8rwfKMMcaUo0iSwzIgRUSSRSQOp7F4eliZ6cDN4ugLHHGrjIqadjpwi/v6FuCjkOHXiUgNEUnGaeReWsr1M8YYUwrFViupaq6ITABm41yO+rKqporIeHf8JGAGzmWsm3AuZb2tqGndWf8FmCYidwA7gJ+406SKyDScRutc4L5yulKp1FVSHqls8YLFXBEqW7xgMVeUs4q5StwhbYwxpmxZfw7GGGPyseRgjDEmnyqZHETkZRHZJyKrQ4Y1FJFPRWSj+79ByLhHRGSTiKwXkeFRFPMTIrJORFaJyIcikhDtMYeM+5WIqIg0DhnmacyFxSsiP3NjShWRv0VLvG4MBX0vuorIYhFZKSLLRaR3tMQsIi1F5HMRWetuz5+7w6N2/ysi5qjd/wqLOWT82e9/qlrl/oBBQHdgdciwvwEPu68fBv7qvu4AfAfUAJKBzYA/SmIeBsS4r/9aGWJ2h7fEuQhhO9A4WmIuZBsPAeYCNdz3TaIl3iJingOMdF9fDsyPlpiBZkB393VdYIMbV9Tuf0XEHLX7X2Exu+/LZP+rkmcOqroAOBg2eDTOYzpw//8oZPhUVc1S1a04V1z1poIVFLOqzlHVvKdqLca55wOiOGbXU8CvOfPmRc9jLiTee4C/qPO4FlQ1734bz+N14ykoZgXqua/r88N9QJ7HrKp71H3opqpmAmtxnnAQtftfYTFH8/5XxHaGMtr/qmRyKMQZj+sAQh/XUdCjP6LN7cBM93XUxiwio4Bdqvpd2KhojbktMFBElojIFyLSyx0erfEC/AJ4QkR2An8HHnGHR1XMItIa6AYsoZLsf2Exh4ra/S805rLc/6w/h9I94qNCichvce75eDNvUAHFPI9ZRGoBv8U5Hc83uoBhnseMsw80APoCvXDuvWlD9MYLztnOL1X1fREZC7wEXEoUxSwidYD3gV+o6lEp/CmQURtzyPCo3f9CY8aJscz2v+p05lApH9chIrcAVwI3qFt5SPTGfB5OfeZ3IrINJ65vROQcojfmNOADdSwFgjgPLIvWeMF5osAH7ut3+aF6ICpiFpFYnAPWm6qaF2dU73+FxBzV+18BMZft/leRjSgV3GDTmjMb8Z7gzAaxv7mvO3JmQ80WPGh4LCTmETh3iieGlYvamMPGbeOHBrGoiLmAbTwe+L37ui3OqbdES7yFxLwWGOy+HgqsiJZt7G6714F/hA2P2v2viJijdv8rLOawMme1/1X4F72CNtzbwB4gBydj3gE0AuYBG93/DUPK/xan9X497lUgURLzJvdgtdL9mxTtMRf25YyGmAvZxnHAG8Bq4BvgkmiJt4iYBwAr3J19CdAjWmJ2Y1NgVcj39vJo3v+KiDlq97/CYg4rc1b7nz0+wxhjTD7Vqc3BGGNMhCw5GGOMyceSgzHGmHwsORhjjMnHkoMxxph8LDkYY4zJx5KDMcaYfP4/XbkuapafzPIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "devs = []\n",
    "for i in range(boot_pred.shape[2]):\n",
    "    devs.append(resid_dev(boot_pred[:,0,i].reshape(-1,1), boot_response[:,0,i].reshape(-1,1))[0])\n",
    "\n",
    "x = np.linspace(np.min(devs)-1, np.max(devs)+1, 1000)\n",
    "plt.plot(x, chi2.pdf(x, 180), label='Predicted pdf w/ 180 df')\n",
    "plt.hist(devs, density = True)\n",
    "plt.legend(loc = 'best')\n",
    "plt.title(\"Predicted Chi2 pdf and Observed Density Histogram\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of insufficient models: 206\n"
     ]
    }
   ],
   "source": [
    "#A quick sanity check that giving random data will produce only insufficient models:\n",
    "deviances2 = resid_dev(np.random.rand(response.shape[0],response.shape[1]), response)\n",
    "\n",
    "#this is the residual deviance so we can construct chi squared tests for all of them (to see if the current model is sufficient)\n",
    "dof = response.shape[0] - data.shape[1]\n",
    "\n",
    "from scipy.stats import chi2\n",
    "test_pvals = 1-chi2.cdf(deviances2, dof)\n",
    "print(\"Number of insufficient models: \" + str(np.sum(test_pvals<.05)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The confusion matrix for this model on the train set is:\n",
      "[[21357.57281553   104.95145631]\n",
      " [  163.83980583   117.75242718]]\n",
      "The confusion matrix for this model on the test set is:\n",
      "[[2373.64563107  196.67475728]\n",
      " [ 203.44660194  197.95145631]]\n"
     ]
    }
   ],
   "source": [
    "#-----------\n",
    "#Look at an aggregated confusion matrix of the models (take average values for the matrices)\n",
    "#-----------\n",
    "from sklearn.metrics import confusion_matrix\n",
    "conf_matrix = np.zeros((2,2,pred_probs.shape[1]))\n",
    "for i in range(pred_probs.shape[1]):\n",
    "    conf_matrix[:,:,i] = confusion_matrix(y_train[:,i], np.round(pred_probs[:,i]))\n",
    "print(\"The confusion matrix for this model on the train set is:\")\n",
    "print(np.mean(conf_matrix, axis = 2))\n",
    "\n",
    "conf_matrix = np.zeros((2,2,pred_probs_test.shape[1]))\n",
    "for i in range(pred_probs_test.shape[1]):\n",
    "    conf_matrix[:,:,i] = confusion_matrix(y_test[:,i], np.round(pred_probs_test[:,i]))\n",
    "print(\"The confusion matrix for this model on the test set is:\")\n",
    "print(np.mean(conf_matrix, axis = 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimensionality Reduction Via VIF\n",
    "Results indicate that any split of the data is sufficient and will yeild simular training. Now look at VIF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (875 of 875) |######################| Elapsed Time: 0:12:30 Time:  0:12:30\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\"\"\"\n",
    "VIF -- find the variance inflation factor for each column of a given matrix\n",
    "\n",
    "@params: --X: a numpy array which we will find VIFs for\n",
    "\"\"\"\n",
    "def VIF(X):\n",
    "    vifs = []\n",
    "    with progressbar.ProgressBar(max_value=int(X.shape[1])) as bar:\n",
    "        for i in range(X.shape[1]):\n",
    "            #Create an array which excludes the ith column\n",
    "            relCols = [j for j in range(X.shape[1]) if j != i] \n",
    "            x = X[:, relCols]\n",
    "            \n",
    "            #fit a linear model onto the excluded column\n",
    "            reg = LinearRegression().fit(x, X[:,i])\n",
    "            vifs.append(1/(1-reg.score(x, X[:,i])))\n",
    "            bar.update(i)\n",
    "    return vifs\n",
    "\n",
    "\n",
    "vifs = VIF(X_train[:, 1:X_train.shape[0]]) #drop bias column for this calculation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The max VIF for the base model was 8.269421633323343\n",
      "The condition number of the base matrix is: 41.52752960126006\n"
     ]
    }
   ],
   "source": [
    "print(\"The max VIF for the base model was \" + str(max(vifs)))\n",
    "print(\"The condition number of the base matrix is: \" + str(np.linalg.cond(X_train[:,1:X_train.shape[1]])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 40.86190704387075\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (874 of 874) |######################| Elapsed Time: 0:12:12 Time:  0:12:12\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 40.38941847312909\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (873 of 873) |######################| Elapsed Time: 0:11:28 Time:  0:11:28\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 40.21960258021064\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (872 of 872) |######################| Elapsed Time: 0:11:00 Time:  0:11:00\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 40.14979175921387\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (871 of 871) |######################| Elapsed Time: 0:11:25 Time:  0:11:25\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 40.05845990783548\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (870 of 870) |######################| Elapsed Time: 0:10:45 Time:  0:10:45\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 40.05594799971745\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (869 of 869) |######################| Elapsed Time: 0:10:54 Time:  0:10:54\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 40.0444410025309\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (868 of 868) |######################| Elapsed Time: 0:10:48 Time:  0:10:48\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 40.00297380345969\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (867 of 867) |######################| Elapsed Time: 0:10:52 Time:  0:10:52\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 39.97294259689092\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (866 of 866) |######################| Elapsed Time: 0:10:48 Time:  0:10:48\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 39.908841039099755\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (865 of 865) |######################| Elapsed Time: 0:10:47 Time:  0:10:47\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 39.89389454520559\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (864 of 864) |######################| Elapsed Time: 0:10:51 Time:  0:10:51\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 39.67396218099479\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (863 of 863) |######################| Elapsed Time: 0:10:51 Time:  0:10:51\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 39.581316225658\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (862 of 862) |######################| Elapsed Time: 0:10:48 Time:  0:10:48\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 39.544635136354486\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (861 of 861) |######################| Elapsed Time: 0:10:39 Time:  0:10:39\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 39.466460826740864\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (860 of 860) |######################| Elapsed Time: 0:10:36 Time:  0:10:36\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 39.45079956018348\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (859 of 859) |######################| Elapsed Time: 0:10:30 Time:  0:10:30\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 37.808918661018815\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (858 of 858) |######################| Elapsed Time: 0:11:01 Time:  0:11:01\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 37.80349283672296\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (857 of 857) |######################| Elapsed Time: 0:10:40 Time:  0:10:40\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 37.68915466769941\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (856 of 856) |######################| Elapsed Time: 0:10:47 Time:  0:10:47\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 37.44954606196908\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (855 of 855) |######################| Elapsed Time: 0:10:37 Time:  0:10:37\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 37.39867452052772\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (854 of 854) |######################| Elapsed Time: 0:10:39 Time:  0:10:39\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 37.3960878026293\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (853 of 853) |######################| Elapsed Time: 0:10:37 Time:  0:10:37\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 37.32014799340859\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (852 of 852) |######################| Elapsed Time: 0:10:54 Time:  0:10:54\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 37.21154368549519\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (851 of 851) |######################| Elapsed Time: 0:10:28 Time:  0:10:28\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.943413039218896\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (850 of 850) |######################| Elapsed Time: 0:10:32 Time:  0:10:32\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.91771296668034\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (849 of 849) |######################| Elapsed Time: 0:10:44 Time:  0:10:44\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.87552574448633\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (848 of 848) |######################| Elapsed Time: 0:10:30 Time:  0:10:30\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.85580142937162\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (847 of 847) |######################| Elapsed Time: 0:10:21 Time:  0:10:21\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.841538516154365\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (846 of 846) |######################| Elapsed Time: 0:09:55 Time:  0:09:55\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.83177917131003\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (845 of 845) |######################| Elapsed Time: 0:09:56 Time:  0:09:56\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.82201148370746\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (844 of 844) |######################| Elapsed Time: 0:09:55 Time:  0:09:55\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.7820195375577\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (843 of 843) |######################| Elapsed Time: 0:09:53 Time:  0:09:53\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.7459424377119\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (842 of 842) |######################| Elapsed Time: 0:09:52 Time:  0:09:52\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.22708376583285\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (841 of 841) |######################| Elapsed Time: 0:09:48 Time:  0:09:48\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.214851163712474\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (840 of 840) |######################| Elapsed Time: 0:09:48 Time:  0:09:48\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 34.198268993021145\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (839 of 839) |######################| Elapsed Time: 0:09:45 Time:  0:09:45\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 31.604936092780413\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (838 of 838) |######################| Elapsed Time: 0:09:40 Time:  0:09:40\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 31.57925959616027\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (837 of 837) |######################| Elapsed Time: 0:09:38 Time:  0:09:38\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 31.46114708658209\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (836 of 836) |######################| Elapsed Time: 0:09:40 Time:  0:09:40\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 31.248830394099308\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (835 of 835) |######################| Elapsed Time: 0:09:32 Time:  0:09:32\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 31.104170839960517\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (834 of 834) |######################| Elapsed Time: 0:09:30 Time:  0:09:30\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 31.07533593715494\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (833 of 833) |######################| Elapsed Time: 0:09:37 Time:  0:09:37\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 31.038386272269086\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (832 of 832) |######################| Elapsed Time: 0:09:49 Time:  0:09:49\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.87911122139734\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (831 of 831) |######################| Elapsed Time: 0:09:36 Time:  0:09:36\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.830366147029878\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (830 of 830) |######################| Elapsed Time: 0:09:32 Time:  0:09:32\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.813922957338246\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (829 of 829) |######################| Elapsed Time: 0:09:28 Time:  0:09:28\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.80222812661647\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (828 of 828) |######################| Elapsed Time: 0:09:49 Time:  0:09:49\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.762265383874848\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (827 of 827) |######################| Elapsed Time: 0:09:49 Time:  0:09:49\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.756112399905653\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (826 of 826) |######################| Elapsed Time: 0:09:40 Time:  0:09:40\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.74732781732603\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (825 of 825) |######################| Elapsed Time: 0:09:28 Time:  0:09:28\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.73320643442421\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (824 of 824) |######################| Elapsed Time: 0:09:26 Time:  0:09:26\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.699152468562513\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (823 of 823) |######################| Elapsed Time: 0:09:21 Time:  0:09:21\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.652380734534955\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (822 of 822) |######################| Elapsed Time: 0:09:21 Time:  0:09:21\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.644152567214746\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (821 of 821) |######################| Elapsed Time: 0:09:18 Time:  0:09:18\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.64235511009651\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (820 of 820) |######################| Elapsed Time: 0:09:17 Time:  0:09:17\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.633506734998505\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (819 of 819) |######################| Elapsed Time: 0:09:17 Time:  0:09:17\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.603331641683273\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (818 of 818) |######################| Elapsed Time: 0:09:15 Time:  0:09:15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.59648496332336\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (817 of 817) |######################| Elapsed Time: 0:09:11 Time:  0:09:11\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.55164412814701\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (816 of 816) |######################| Elapsed Time: 0:09:12 Time:  0:09:12\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.54494680371747\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (815 of 815) |######################| Elapsed Time: 0:09:09 Time:  0:09:09\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.385216262276312\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (814 of 814) |######################| Elapsed Time: 0:09:06 Time:  0:09:06\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.380242428712414\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (813 of 813) |######################| Elapsed Time: 0:09:04 Time:  0:09:04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.364088327821364\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (812 of 812) |######################| Elapsed Time: 0:09:03 Time:  0:09:03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.349574091635525\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (811 of 811) |######################| Elapsed Time: 0:09:02 Time:  0:09:02\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.315323537934816\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (810 of 810) |######################| Elapsed Time: 0:09:00 Time:  0:09:00\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.314340834770896\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (809 of 809) |######################| Elapsed Time: 0:08:58 Time:  0:08:58\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.258727625329506\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (808 of 808) |######################| Elapsed Time: 0:09:01 Time:  0:09:01\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.253485303337282\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (807 of 807) |######################| Elapsed Time: 0:08:56 Time:  0:08:56\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 30.24661439522932\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (806 of 806) |######################| Elapsed Time: 0:08:56 Time:  0:08:56\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number of this iteration is: 29.351846752477485\n",
      "['c-83', 'c-13', 'g-50', 'g-100', 'g-75', 'g-441', 'g-651', 'g-369', 'g-64', 'g-534', 'g-300', 'c-18', 'g-392', 'g-635', 'g-629', 'g-253', 'c-52', 'g-502', 'g-37', 'c-73', 'g-195', 'cp_dose', 'g-121', 'g-744', 'c-26', 'g-80', 'g-410', 'g-498', 'g-139', 'g-206', 'g-727', 'g-102', 'g-615', 'c-42', 'g-385', 'g-770', 'c-94', 'g-664', 'g-38', 'c-55', 'c-10', 'g-351', 'g-314', 'c-6', 'g-327', 'g-178', 'g-140', 'g-761', 'g-566', 'g-553', 'g-287', 'g-406', 'g-374', 'g-235', 'g-175', 'g-76', 'g-439', 'g-588', 'g-63', 'g-172', 'c-63', 'g-404', 'g-248', 'g-509', 'g-196', 'g-454', 'g-503', 'g-522', 'g-21', 'c-38']\n"
     ]
    }
   ],
   "source": [
    "#eliminate the highest VIF and see what the condition number was, then repeat:\n",
    "\n",
    "n = 100 #max number of columns to remove\n",
    "col_to_elim = [] #a list of column names which are eliminated in this process\n",
    "temp_data = X_train[:, 1:X_train.shape[0]] #Eliminate the bias column while we do this calculation\n",
    "col_names_data_temp = col_names_data.copy()\n",
    "\n",
    "#-----------\n",
    "#Dimensionality reduction using the VIF\n",
    "#Stopping criterian: if the condition number is less than 30 or if we have eliminated 100 columns (this is due to time constraints)\n",
    "#-----------\n",
    "for i in range(n):\n",
    "    topInd = np.argsort(vifs)[-1] #find the index of the highest VIF\n",
    "    col_to_elim.append(col_names_data_temp[topInd]) #record the name of the col to be eliminated\n",
    "    \n",
    "    #-----------\n",
    "    #Eliminate highest VIF column from both the column names and the temp data\n",
    "    #-----------\n",
    "    col_names_data_temp.pop(topInd) \n",
    "    relCols = [j for j in range(temp_data.shape[1]) if j != topInd] \n",
    "    temp_data = temp_data[:, relCols]\n",
    "    \n",
    "    #find an state the condition number\n",
    "    cond = np.linalg.cond(temp_data)\n",
    "    print(\"The condition number of this iteration is: \" + str(cond))\n",
    "    \n",
    "    #stopping condition\n",
    "    if cond < 30:\n",
    "        break\n",
    "    #calculate new VIFs\n",
    "    vifs = VIF(temp_data)\n",
    "print(col_to_elim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of columns eliminated to acheive a condition number less than 30:\n",
      "70\n"
     ]
    }
   ],
   "source": [
    "col_to_elim = ['c-83', 'c-13', 'g-50', 'g-100', 'g-75', 'g-441', 'g-651', 'g-369', 'g-64', 'g-534', 'g-300', 'c-18', \n",
    "                'g-392', 'g-635', 'g-629', 'g-253', 'c-52', 'g-502', 'g-37', 'c-73', 'g-195', 'cp_dose', 'g-121', 'g-744', \n",
    "                'c-26', 'g-80', 'g-410', 'g-498', 'g-139', 'g-206', 'g-727', 'g-102', 'g-615', 'c-42', 'g-385', 'g-770', 'c-94', \n",
    "                'g-664', 'g-38', 'c-55', 'c-10', 'g-351', 'g-314', 'c-6', 'g-327', 'g-178', 'g-140', 'g-761', 'g-566', 'g-553', 'g-287', 'g-406', 'g-374', \n",
    "                'g-235', 'g-175', 'g-76', 'g-439', 'g-588', 'g-63', 'g-172', 'c-63', 'g-404', 'g-248', 'g-509', 'g-196', 'g-454', 'g-503', 'g-522', 'g-21', 'c-38']\n",
    "print(\"Number of columns eliminated to acheive a condition number less than 30:\")\n",
    "print(len(col_to_elim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\JAT\\anaconda3\\envs\\pytorch\\lib\\site-packages\\torch\\nn\\functional.py:1625: UserWarning: nn.functional.sigmoid is deprecated. Use torch.sigmoid instead.\n",
      "  warnings.warn(\"nn.functional.sigmoid is deprecated. Use torch.sigmoid instead.\")\n",
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:08 Time:  0:00:08\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 1 were:\n",
      "accuracy: 0.9971476949014593\n",
      "loss: 0.015840276557126803\n",
      "The train statistics for fold 1 were:\n",
      "accuracy: 0.9970562188401996\n",
      "loss: 0.018172875978052616\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:08 Time:  0:00:08\n",
      "  1% (118 of 10020) |                    | Elapsed Time: 0:00:00 ETA:   0:00:08"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 2 were:\n",
      "accuracy: 0.9971352047119352\n",
      "loss: 0.01592043881209757\n",
      "The train statistics for fold 2 were:\n",
      "accuracy: 0.9971853680136444\n",
      "loss: 0.017392672220798763\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:08 Time:  0:00:08\n",
      "  1% (113 of 10020) |                    | Elapsed Time: 0:00:00 ETA:   0:00:08"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 3 were:\n",
      "accuracy: 0.9971486032788792\n",
      "loss: 0.015854994797086467\n",
      "The train statistics for fold 3 were:\n",
      "accuracy: 0.9970582688270797\n",
      "loss: 0.018030608903515984\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:08 Time:  0:00:08\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 4 were:\n",
      "accuracy: 0.9971436072030696\n",
      "loss: 0.01583442405078136\n",
      "The train statistics for fold 4 were:\n",
      "accuracy: 0.997115668459722\n",
      "loss: 0.018159668187837343\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:08 Time:  0:00:08\n",
      "  1% (111 of 10020) |                    | Elapsed Time: 0:00:00 ETA:   0:00:08"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 5 were:\n",
      "accuracy: 0.9971411091651649\n",
      "loss: 0.01587850388842161\n",
      "The train statistics for fold 5 were:\n",
      "accuracy: 0.9971546182104436\n",
      "loss: 0.017614141408656095\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 6 were:\n",
      "accuracy: 0.9971352047119354\n",
      "loss: 0.01591739144655848\n",
      "The train statistics for fold 6 were:\n",
      "accuracy: 0.997193567961165\n",
      "loss: 0.0174160413695751\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 7 were:\n",
      "accuracy: 0.9971406549764549\n",
      "loss: 0.015856526872912745\n",
      "The train statistics for fold 7 were:\n",
      "accuracy: 0.9971628181579638\n",
      "loss: 0.018024778622831847\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:10 Time:  0:00:10\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 8 were:\n",
      "accuracy: 0.9971554161095285\n",
      "loss: 0.015824282645041536\n",
      "The train statistics for fold 8 were:\n",
      "accuracy: 0.997023419050118\n",
      "loss: 0.01835083991691873\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:10 Time:  0:00:10\n",
      "N/A% (0 of 10020) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 9 were:\n",
      "accuracy: 0.9971420175425847\n",
      "loss: 0.01584605443776546\n",
      "The train statistics for fold 9 were:\n",
      "accuracy: 0.9970767187090002\n",
      "loss: 0.018226616915214707\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:09 Time:  0:00:09\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 10 were:\n",
      "accuracy: 0.9971495116562991\n",
      "loss: 0.015865012983830567\n",
      "The train statistics for fold 10 were:\n",
      "accuracy: 0.9970644187877198\n",
      "loss: 0.01789466938557657\n"
     ]
    }
   ],
   "source": [
    "#since I don't want to have to rerun the code above I will have to create another data set:\n",
    "ind_to_keep_bool = np.ones(data.shape[1], dtype = bool)\n",
    "for i in range(data.shape[1]):\n",
    "    if i != 0:\n",
    "        if col_names_data[i-1] in col_to_elim:\n",
    "            ind_to_keep_bool[i] = 0\n",
    "\n",
    "red_data = data[:, ind_to_keep_bool]\n",
    "\n",
    "#Find k folds for k-fold cross validation:\n",
    "kf = KFold(n_splits=10, random_state = 0, shuffle = True)\n",
    "kf.get_n_splits(red_data)\n",
    "kcross = kf.split(red_data, response) #creates an iteratable that we will use later\n",
    "\n",
    "fold = 1 #keeps track of which fold we are on\n",
    "\n",
    "#do training loop for each fold in the K-fold cross:\n",
    "for train_index, test_index in kcross:\n",
    "    X_train, X_test = red_data[train_index], red_data[test_index]\n",
    "    y_train, y_test = response[train_index], response[test_index]\n",
    "    \n",
    "    #Initialize necessary constants like batch size, etc as well as the model:\n",
    "    model = log_reg(red_data.shape[1], response.shape[1]).cuda() #again requires a working gpu\n",
    "    opt = torch.optim.Adam(model.parameters()) #Initializes and tells the optimizer what parameters to keep track of\n",
    "    batch_size = 64 #This could probably be increased, but 64 seems to work well enough\n",
    "    epochs = 30\n",
    "    loss = nn.BCELoss() #Binary cross entropy loss, the loss function defined in the write up\n",
    "    batchesPerEpoch = int(np.floor(X_train.shape[0]/batch_size))\n",
    "    \n",
    "    #-----------\n",
    "    #training loop\n",
    "    #-----------\n",
    "    with progressbar.ProgressBar(max_value=epochs*batchesPerEpoch) as bar:\n",
    "        for epoch in range(epochs):\n",
    "            for i in range(batchesPerEpoch):\n",
    "                opt.zero_grad()\n",
    "                batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "                batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "                #forward pass:\n",
    "                out = model(batch_x)\n",
    "                batch_loss = torch.mean(loss(out, batch_y.float()))\n",
    "                batch_loss.backward()\n",
    "                opt.step()\n",
    "                bar.update(i+batchesPerEpoch*epoch)\n",
    "    #-----------            \n",
    "    #Get train stats:\n",
    "    #-----------\n",
    "    accuracy = []\n",
    "    test_loss = []\n",
    "    model.eval()\n",
    "    with torch.no_grad():\n",
    "        for i in range(batchesPerEpoch):\n",
    "            batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            out = model(batch_x)\n",
    "            test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "            accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "    print(\"The train statistics for fold \" + str(fold) + \" were:\")\n",
    "    print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "    print(\"loss: \" + str(np.mean(test_loss)))    \n",
    "    \n",
    "    #-----------\n",
    "    #Get test stats:\n",
    "    #-----------\n",
    "    accuracy = []\n",
    "    test_loss = []\n",
    "    batches = int(np.floor(X_test.shape[0]/batch_size))\n",
    "    with torch.no_grad():\n",
    "        for i in range(batches):\n",
    "            batch_x = torch.from_numpy(X_test[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_test[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            out = model(batch_x)\n",
    "            test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "            accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "    print(\"The train statistics for fold \" + str(fold) + \" were:\")\n",
    "    print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "    print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "    fold += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (50 of 50) |########################| Elapsed Time: 0:08:09 Time:  0:08:09\n",
      "100% (10020 of 10020) |##################| Elapsed Time: 0:00:10 Time:  0:00:10\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for the full model:\n",
      "accuracy: 0.9971540535433987\n",
      "loss: 0.015822401142004363\n",
      "The train statistics for the full model were:\n",
      "accuracy: 0.9970582688270795\n",
      "loss: 0.01816255523747689\n",
      "The train statistics for the full model were:\n",
      "accuracy: 0.9970582688270795\n",
      "loss: 0.01816255523747689\n",
      "The lowest p-value:\n",
      "0.0\n",
      "The highest p-value:\n",
      "0.4999999211531815\n",
      "The proportion of parameters that are significant (across all models):\n",
      "0.7214158375292106\n",
      "The number of predictors that are non-significant in all models:\n",
      "0\n",
      "The number of predictors that are non-significant in all but 10 or less models:\n",
      "0\n",
      "The number of predictors that are non-significant in all but 50 or less models:\n",
      "74\n",
      "Number of insufficient models: 0\n",
      "The confusion matrix for this model on the train set is:\n",
      "[[21357.5631068    104.96116505]\n",
      " [  164.16990291   117.4223301 ]]\n",
      "The confusion matrix for this model on the test set is:\n",
      "[[2373.66019417  196.66019417]\n",
      " [ 203.49514563  197.90291262]]\n"
     ]
    }
   ],
   "source": [
    "# To get p-values of the coefficients we need standard errors, to do this we will use bootstrapping:\n",
    "n = 50 #number of bootstrapped samples to consider\n",
    "boot_coeff = torch.zeros(206,806,n) #3D array of coefficients dimension each slice along dim 2 is a mxn matrix of parameters\n",
    "\"\"\"\n",
    "generateBootSample -- generates a sample of rows from X and y with replacement\n",
    "\n",
    "@params -- X a numpy array to sample \n",
    "        -- y a second numpy array to generate the same sample from\n",
    "        -- size an int for the size of the bootstrapped sample\n",
    "\"\"\"\n",
    "def generateBootSample(X,y, size):\n",
    "    samples = torch.randint(0, int(X.shape[0]), (size, 1)).view(-1).numpy()\n",
    "    \n",
    "    return X[samples, :], y[samples,:]\n",
    "\n",
    "#-----------\n",
    "#train a model for each bootstrapped sample required, store the parameters in boot_coeff for each one\n",
    "#-----------\n",
    "with progressbar.ProgressBar(max_value=n) as bar:\n",
    "    for sample in range(n):\n",
    "        X_train, y_train = generateBootSample(red_data, response, int(red_data.shape[0]*9/10))\n",
    "        #Initialize necessary constants like batch size, etc as well as the model:\n",
    "        model = log_reg(red_data.shape[1], response.shape[1]).cuda()\n",
    "        opt = torch.optim.Adam(model.parameters())\n",
    "        batch_size = 64\n",
    "        epochs = 30\n",
    "        loss = nn.BCELoss()\n",
    "        batchesPerEpoch = int(np.floor(X_train.shape[0]/batch_size))\n",
    "\n",
    "        for epoch in range(epochs):\n",
    "            for i in range(batchesPerEpoch):\n",
    "                opt.zero_grad()\n",
    "                batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "                batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "                #forward pass:\n",
    "                out = model(batch_x)\n",
    "                batch_loss = torch.mean(loss(out, batch_y.float()))\n",
    "                batch_loss.backward()\n",
    "                opt.step()\n",
    "                \n",
    "        bar.update(sample)\n",
    "        for param in model.parameters():\n",
    "            boot_coeff[:,:,sample] = param.detach().cpu()\n",
    "standard_errors = torch.std(boot_coeff, dim = 2) #estimate the standard error of the coefficients using the bootstrapped sample\n",
    "\n",
    "#-----------\n",
    "#train a final model:\n",
    "#-----------\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(red_data, response, test_size=0.1, random_state=0) #create a final split\n",
    "\n",
    "#Initialize necessary constants like batch size, etc as well as the model:\n",
    "model = log_reg(red_data.shape[1], response.shape[1]).cuda()\n",
    "opt = torch.optim.Adam(model.parameters())\n",
    "batch_size = 64\n",
    "epochs = 30\n",
    "loss = nn.BCELoss()\n",
    "batchesPerEpoch = int(np.floor(X_train.shape[0]/batch_size))\n",
    "\n",
    "#-----------\n",
    "#training loop\n",
    "#-----------\n",
    "with progressbar.ProgressBar(max_value=epochs*batchesPerEpoch) as bar:\n",
    "    for epoch in range(epochs):\n",
    "        for i in range(batchesPerEpoch):\n",
    "            opt.zero_grad()\n",
    "            batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            #forward pass:\n",
    "            out = model(batch_x)\n",
    "            batch_loss = torch.mean(loss(out, batch_y.float()))\n",
    "            batch_loss.backward()\n",
    "            opt.step()\n",
    "            bar.update(i+batchesPerEpoch*epoch)\n",
    "\n",
    "#-----------\n",
    "#Get train stats:\n",
    "#-----------\n",
    "accuracy = []\n",
    "test_loss = []\n",
    "model.eval()\n",
    "pred_probs = np.zeros(y_train.shape)\n",
    "with torch.no_grad():\n",
    "    for i in range(batchesPerEpoch):\n",
    "        batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "        batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "        out = model(batch_x)\n",
    "        pred_probs[i*batch_size:(i+1)*batch_size,:] = out.cpu().numpy()\n",
    "        test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "        accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "print(\"The train statistics for the full model:\")\n",
    "print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "\n",
    "#-----------\n",
    "#Get test stats:\n",
    "#-----------\n",
    "accuracy = []\n",
    "test_loss = []\n",
    "batches = int(np.floor(X_test.shape[0]/batch_size))\n",
    "pred_probs_test = np.zeros(y_test.shape)\n",
    "with torch.no_grad():\n",
    "    for i in range(batches):\n",
    "        batch_x = torch.from_numpy(X_test[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "        batch_y = torch.from_numpy(y_test[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "        out = model(batch_x)\n",
    "        pred_probs_test[i*batch_size:(i+1)*batch_size,:] = out.cpu().numpy()\n",
    "        test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "        accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "print(\"The train statistics for the full model were:\")\n",
    "print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "\n",
    "#-----------\n",
    "#Get test stats:\n",
    "#-----------\n",
    "accuracy = []\n",
    "test_loss = []\n",
    "batches = int(np.floor(X_test.shape[0]/batch_size))\n",
    "with torch.no_grad():\n",
    "    for i in range(batches):\n",
    "        batch_x = torch.from_numpy(X_test[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "        batch_y = torch.from_numpy(y_test[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "        out = model(batch_x)\n",
    "        test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "        accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "print(\"The train statistics for the full model were:\")\n",
    "print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "\n",
    "#-----------\n",
    "#get p-values for coefficients:\n",
    "#-----------\n",
    "#get z-scores of the paramters\n",
    "for param in model.parameters():\n",
    "    z_scores = torch.abs(torch.div(param.detach().cpu(), standard_errors)).numpy()\n",
    "\n",
    "#get the p-values\n",
    "from scipy.stats import norm\n",
    "p_vals = 1 - norm.cdf(z_scores)\n",
    "\n",
    "print(\"The lowest p-value:\")\n",
    "print(np.min(p_vals))\n",
    "print(\"The highest p-value:\")\n",
    "print(np.max(p_vals))\n",
    "print(\"The proportion of parameters that are significant (across all models):\")\n",
    "print(np.sum(p_vals < .05)/(p_vals.shape[0]*p_vals.shape[1]))\n",
    "\n",
    "count = 0\n",
    "for i in range(p_vals.shape[1]):\n",
    "    if np.sum(p_vals[:,i]<.05) == 0:\n",
    "        count += 1\n",
    "print(\"The number of predictors that are non-significant in all models:\")\n",
    "print(count)\n",
    "count = 0\n",
    "for i in range(p_vals.shape[1]):\n",
    "    if np.sum(p_vals[:,i]<.05) <= 10:\n",
    "        count += 1\n",
    "print(\"The number of predictors that are non-significant in all but 10 or less models:\")\n",
    "print(count)\n",
    "\n",
    "count = 0\n",
    "for i in range(p_vals.shape[1]):\n",
    "    if np.sum(p_vals[:,i]<.05) <= 50:\n",
    "        count += 1\n",
    "print(\"The number of predictors that are non-significant in all but 50 or less models:\")\n",
    "print(count)\n",
    "\n",
    "deviances = resid_dev(pred_probs, y_train) #find the residual deviances of the model\n",
    "\n",
    "#this is the residual deviance so we can construct chi squared tests for all of them (to see if the current model is sufficient)\n",
    "dof = response.shape[0] - red_data.shape[1] #degrees of freedom\n",
    "\n",
    "#-----------\n",
    "#Do the chi-squared test for model sufficiency\n",
    "#-----------\n",
    "from scipy.stats import chi2\n",
    "test_pvals = 1-chi2.cdf(deviances, dof)\n",
    "print(\"Number of insufficient models: \" + str(np.sum(test_pvals<.05)))\n",
    "\n",
    "conf_matrix = np.zeros((2,2,pred_probs.shape[1]))\n",
    "for i in range(pred_probs.shape[1]):\n",
    "    conf_matrix[:,:,i] = confusion_matrix(y_train[:,i], np.round(pred_probs[:,i]))\n",
    "print(\"The confusion matrix for this model on the train set is:\")\n",
    "print(np.mean(conf_matrix, axis = 2))\n",
    "\n",
    "conf_matrix = np.zeros((2,2,pred_probs_test.shape[1]))\n",
    "for i in range(pred_probs_test.shape[1]):\n",
    "    conf_matrix[:,:,i] = confusion_matrix(y_test[:,i], np.round(pred_probs_test[:,i]))\n",
    "print(\"The confusion matrix for this model on the test set is:\")\n",
    "print(np.mean(conf_matrix, axis = 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Faster Dim Reducation using Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\JAT\\anaconda3\\envs\\pytorch\\lib\\site-packages\\torch\\nn\\functional.py:1625: UserWarning: nn.functional.sigmoid is deprecated. Use torch.sigmoid instead.\n",
      "  warnings.warn(\"nn.functional.sigmoid is deprecated. Use torch.sigmoid instead.\")\n",
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:52 Time:  0:00:52\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 1 where:\n",
      "accuracy: 0.9965692855793269\n",
      "loss: 0.23948609190965128\n",
      "The train statistics for fold 1 where:\n",
      "accuracy: 0.9965375721595382\n",
      "loss: 0.23953108449239988\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:52 Time:  0:00:52\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 2 where:\n",
      "accuracy: 0.9965672417301321\n",
      "loss: 0.23949746631398172\n",
      "The train statistics for fold 2 where:\n",
      "accuracy: 0.9965560220414589\n",
      "loss: 0.2395093598075815\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:53 Time:  0:00:53\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 3 where:\n",
      "accuracy: 0.996565197880937\n",
      "loss: 0.2394983736340871\n",
      "The train statistics for fold 3 where:\n",
      "accuracy: 0.9965744719233798\n",
      "loss: 0.23948241729994077\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:52 Time:  0:00:52\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 4 where:\n",
      "accuracy: 0.9965654249752922\n",
      "loss: 0.2394793821487598\n",
      "The train statistics for fold 4 where:\n",
      "accuracy: 0.9965683219627393\n",
      "loss: 0.23947950150515582\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:52 Time:  0:00:52\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 5 where:\n",
      "accuracy: 0.9965574766728679\n",
      "loss: 0.23952607315279054\n",
      "The train statistics for fold 5 where:\n",
      "accuracy: 0.9966441714773024\n",
      "loss: 0.23941944136812882\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:52 Time:  0:00:52\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 6 where:\n",
      "accuracy: 0.9965636082204523\n",
      "loss: 0.2394842043697477\n",
      "The train statistics for fold 6 where:\n",
      "accuracy: 0.9965929218053005\n",
      "loss: 0.23944164570924398\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:52 Time:  0:00:52\n",
      "  0% (63 of 33400) |                     | Elapsed Time: 0:00:00 ETA:   0:00:53"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 7 where:\n",
      "accuracy: 0.9965617914656125\n",
      "loss: 0.23948985382825314\n",
      "The train statistics for fold 7 where:\n",
      "accuracy: 0.9966011217528208\n",
      "loss: 0.23944759932724205\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:53 Time:  0:00:53\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 8 where:\n",
      "accuracy: 0.9965699668623917\n",
      "loss: 0.2394895251104218\n",
      "The train statistics for fold 8 where:\n",
      "accuracy: 0.9965293722120177\n",
      "loss: 0.2395499376831828\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:54 Time:  0:00:54\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 9 where:\n",
      "accuracy: 0.9965690584849718\n",
      "loss: 0.23947146343078443\n",
      "The train statistics for fold 9 where:\n",
      "accuracy: 0.9965396221464181\n",
      "loss: 0.23951299528817874\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:54 Time:  0:00:54\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 10 where:\n",
      "accuracy: 0.9965729190890064\n",
      "loss: 0.23946867196145885\n",
      "The train statistics for fold 10 where:\n",
      "accuracy: 0.9965109223300967\n",
      "loss: 0.23953692373391744\n"
     ]
    }
   ],
   "source": [
    "#Find k folds for k-fold cross validation:\n",
    "kf = KFold(n_splits=10, random_state = 0, shuffle = True)\n",
    "kf.get_n_splits(data)\n",
    "kcross = kf.split(data, response)\n",
    "import progressbar\n",
    "fold = 1\n",
    "\n",
    "#-----------\n",
    "#do training loop for each fold in the K-fold cross:\n",
    "#-----------\n",
    "for train_index, test_index in kcross:\n",
    "    X_train, X_test = data[train_index], data[test_index]\n",
    "    y_train, y_test = response[train_index], response[test_index]\n",
    "    \n",
    "    #Initialize necessary constants like batch size, etc as well as the model:\n",
    "    model = log_reg(data.shape[1], response.shape[1]).cuda()\n",
    "    opt = torch.optim.Adam(model.parameters())\n",
    "    batch_size = 64\n",
    "    epochs = 100\n",
    "    loss = nn.BCELoss()\n",
    "    l1 = torch.nn.L1Loss(size_average=False)\n",
    "    batchesPerEpoch = int(np.floor(X_train.shape[0]/batch_size))\n",
    "    with progressbar.ProgressBar(max_value=epochs*batchesPerEpoch) as bar:\n",
    "        for epoch in range(epochs):\n",
    "            for i in range(batchesPerEpoch):\n",
    "                opt.zero_grad()\n",
    "                batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "                batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "                #forward pass:\n",
    "                out = model(batch_x)\n",
    "                \n",
    "                #-----------\n",
    "                #create an l1 penalty for the loss\n",
    "                #-----------\n",
    "                l1_loss = 0\n",
    "                for param in model.parameters():\n",
    "                    l1_target = torch.zeros(param.shape).cuda()\n",
    "                    l1_loss += l1(param, l1_target)\n",
    "                scale_l1 = .001 #this is the lambda in the l1 penalty\n",
    "                batch_loss = torch.mean(loss(out, batch_y.float())) + scale_l1*l1_loss\n",
    "                \n",
    "                batch_loss.backward()\n",
    "                opt.step()\n",
    "                bar.update(i+batchesPerEpoch*epoch)\n",
    "                \n",
    "    #-----------\n",
    "    #Get train stats:\n",
    "    #-----------\n",
    "    accuracy = []\n",
    "    test_loss = []\n",
    "    model.eval()\n",
    "    with torch.no_grad():\n",
    "        for i in range(batchesPerEpoch):\n",
    "            batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            out = model(batch_x)\n",
    "            test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "            accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "    print(\"The train statistics for fold \" + str(fold) + \" where:\")\n",
    "    print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "    print(\"loss: \" + str(np.mean(test_loss)))    \n",
    "    \n",
    "    #-----------\n",
    "    #Get test stats:\n",
    "    #-----------\n",
    "    accuracy = []\n",
    "    test_loss = []\n",
    "    batches = int(np.floor(X_test.shape[0]/batch_size))\n",
    "    with torch.no_grad():\n",
    "        for i in range(batches):\n",
    "            batch_x = torch.from_numpy(X_test[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_test[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            out = model(batch_x)\n",
    "            test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "            accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "    print(\"The train statistics for fold \" + str(fold) + \" where:\")\n",
    "    print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "    print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "    fold += 1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage of zeros in this parameter set:\n",
      "0.006766192312807554\n"
     ]
    }
   ],
   "source": [
    "for param in model.parameters():\n",
    "    temp = param.detach().cpu().numpy()\n",
    "    print(\"Percentage of zeros in this parameter set:\")\n",
    "    print(np.sum(np.isclose(temp, 0, atol = 1e-06))/(temp.shape[0]*temp.shape[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:59 Time:  0:00:59\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 1 where:\n",
      "accuracy: 0.996963067190861\n",
      "loss: 0.018101531683194068\n",
      "The train statistics for fold 1 where:\n",
      "accuracy: 0.9969270696667539\n",
      "loss: 0.018469580449163914\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:58 Time:  0:00:58\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 2 where:\n",
      "accuracy: 0.9969508040956921\n",
      "loss: 0.01817018351333584\n",
      "The train statistics for fold 2 where:\n",
      "accuracy: 0.9970213690632381\n",
      "loss: 0.018039979317502397\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:01:02 Time:  0:01:02\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 3 where:\n",
      "accuracy: 0.9969596607755363\n",
      "loss: 0.018146387987329576\n",
      "The train statistics for fold 3 where:\n",
      "accuracy: 0.9969434695617947\n",
      "loss: 0.01806398800800781\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:59 Time:  0:00:59\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 4 where:\n",
      "accuracy: 0.996959433681181\n",
      "loss: 0.01812299818335893\n",
      "The train statistics for fold 4 where:\n",
      "accuracy: 0.9969475695355549\n",
      "loss: 0.01822873021199091\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:01:02 Time:  0:01:02\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 5 where:\n",
      "accuracy: 0.9969521666618221\n",
      "loss: 0.018125124808267026\n",
      "The train statistics for fold 5 where:\n",
      "accuracy: 0.9970049691681971\n",
      "loss: 0.018184709941615928\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:01:01 Time:  0:01:01\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 6 where:\n",
      "accuracy: 0.996953075039242\n",
      "loss: 0.018190557683656315\n",
      "The train statistics for fold 6 where:\n",
      "accuracy: 0.9970070191550773\n",
      "loss: 0.01782207007243021\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:01:03 Time:  0:01:03\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 7 where:\n",
      "accuracy: 0.9969542105110167\n",
      "loss: 0.01813397277425417\n",
      "The train statistics for fold 7 where:\n",
      "accuracy: 0.9969885692731566\n",
      "loss: 0.01819022186100483\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:59 Time:  0:00:59\n",
      "  0% (56 of 33400) |                     | Elapsed Time: 0:00:00 ETA:   0:00:59"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 8 where:\n",
      "accuracy: 0.9969648839457009\n",
      "loss: 0.01808035687381725\n",
      "The train statistics for fold 8 where:\n",
      "accuracy: 0.9968901699029128\n",
      "loss: 0.018620314957523667\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:58 Time:  0:00:58\n",
      "N/A% (0 of 33400) |                      | Elapsed Time: 0:00:00 ETA:  --:--:--"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 9 where:\n",
      "accuracy: 0.9969653381344107\n",
      "loss: 0.0181316642709701\n",
      "The train statistics for fold 9 where:\n",
      "accuracy: 0.9968819699553922\n",
      "loss: 0.018370094306364253\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% (33400 of 33400) |##################| Elapsed Time: 0:00:58 Time:  0:00:58\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The train statistics for fold 10 where:\n",
      "accuracy: 0.996961023341666\n",
      "loss: 0.018150247307215445\n",
      "The train statistics for fold 10 where:\n",
      "accuracy: 0.9969106697717136\n",
      "loss: 0.01811806668804304\n"
     ]
    }
   ],
   "source": [
    "#Find k folds for k-fold cross validation:\n",
    "kf = KFold(n_splits=10, random_state = 0, shuffle = True)\n",
    "kf.get_n_splits(data)\n",
    "kcross = kf.split(data, response)\n",
    "import progressbar\n",
    "fold = 1\n",
    "\n",
    "#-----------\n",
    "#do training loop for each fold in the K-fold cross:\n",
    "#-----------\n",
    "for train_index, test_index in kcross:\n",
    "    X_train, X_test = data[train_index], data[test_index]\n",
    "    y_train, y_test = response[train_index], response[test_index]\n",
    "    \n",
    "    #Initialize necessary constants like batch size, etc as well as the model:\n",
    "    model = log_reg(data.shape[1], response.shape[1]).cuda()\n",
    "    opt = torch.optim.Adam(model.parameters())\n",
    "    batch_size = 64\n",
    "    epochs = 100\n",
    "    loss = nn.BCELoss()\n",
    "    l1 = torch.nn.L1Loss(size_average=False)\n",
    "    batchesPerEpoch = int(np.floor(X_train.shape[0]/batch_size))\n",
    "    with progressbar.ProgressBar(max_value=epochs*batchesPerEpoch) as bar:\n",
    "        for epoch in range(epochs):\n",
    "            for i in range(batchesPerEpoch):\n",
    "                opt.zero_grad()\n",
    "                batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "                batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "                #forward pass:\n",
    "                out = model(batch_x)\n",
    "                \n",
    "                #-----------\n",
    "                #create an l1 penalty for the loss\n",
    "                #-----------\n",
    "                l1_loss = 0\n",
    "                for param in model.parameters():\n",
    "                    l1_target = torch.zeros(param.shape).cuda()\n",
    "                    l1_loss += l1(param, l1_target)\n",
    "                scale_l1 = .000001 #this is the lambda in the l1 penalty\n",
    "                batch_loss = torch.mean(loss(out, batch_y.float())) + scale_l1*l1_loss\n",
    "                \n",
    "                batch_loss.backward()\n",
    "                opt.step()\n",
    "                bar.update(i+batchesPerEpoch*epoch)\n",
    "                \n",
    "    #-----------\n",
    "    #Get train stats:\n",
    "    #-----------\n",
    "    accuracy = []\n",
    "    test_loss = []\n",
    "    model.eval()\n",
    "    with torch.no_grad():\n",
    "        for i in range(batchesPerEpoch):\n",
    "            batch_x = torch.from_numpy(X_train[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_train[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            out = model(batch_x)\n",
    "            test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "            accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "    print(\"The train statistics for fold \" + str(fold) + \" where:\")\n",
    "    print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "    print(\"loss: \" + str(np.mean(test_loss)))    \n",
    "    \n",
    "    #-----------\n",
    "    #Get test stats:\n",
    "    #-----------\n",
    "    accuracy = []\n",
    "    test_loss = []\n",
    "    batches = int(np.floor(X_test.shape[0]/batch_size))\n",
    "    with torch.no_grad():\n",
    "        for i in range(batches):\n",
    "            batch_x = torch.from_numpy(X_test[i*batch_size:(i+1)*batch_size,:]).float().cuda()\n",
    "            batch_y = torch.from_numpy(y_test[i*batch_size:(i+1)*batch_size,:]).cuda()\n",
    "            out = model(batch_x)\n",
    "            test_loss.append(torch.mean(loss(out, batch_y.float())).item()) \n",
    "            accuracy.append(torch.sum(torch.round(out) == batch_y).item()/(out.shape[0]*out.shape[1]))\n",
    "    print(\"The train statistics for fold \" + str(fold) + \" where:\")\n",
    "    print(\"accuracy: \" + str(np.mean(accuracy)))\n",
    "    print(\"loss: \" + str(np.mean(test_loss)))  \n",
    "    fold += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of zeros in the parameter set:\n",
      "545\n",
      "Proportion of zeros in this parameter set:\n",
      "0.0030201267899100056\n"
     ]
    }
   ],
   "source": [
    "for param in model.parameters():\n",
    "    temp = param.detach().cpu().numpy()\n",
    "    print(\"Number of zeros in the parameter set:\")\n",
    "    print(np.sum(np.isclose(temp, 0, atol = 1e-06)))\n",
    "    print(\"Proportion of zeros in this parameter set:\")\n",
    "    print(np.sum(np.isclose(temp, 0, atol = 1e-06))/(temp.shape[0]*temp.shape[1]))"
   ]
  },
  {
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}