PLSDA-cancer-classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PLSDA Cancer Classification on GBM Cell Lines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import data from csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import csv\n",
    "import time\n",
    "import scipy\n",
    "from scipy import stats\n",
    "from sklearn.cross_decomposition import PLSRegression\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PARSING DATA\n",
    "After performing preprocessing on our images, the input data for this program will be a '.csv' file. The following code parses the data in that file and stores it into two arrays. X will be the array that contains the pixel color values for all of the pixels in all 100 of our images. Y will be the array which lists the label (0) mes or (1) pro for each of the cell images. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "csv_filename = '/Users/chiara.mooney/Desktop/data.csv'\n",
    "# csv_filename = '/Users/Ryan/Dropbox/Sharing/project data/data.csv'\n",
    "# csv_filename = 'D:/Ryan/Dropbox/Sharing/project data/data.csv'\n",
    "\n",
    "data = []\n",
    "with open(csv_filename, 'r') as f:\n",
    "    reader = csv.reader(f)\n",
    "    for row in reader:\n",
    "        data.append(row[1:])\n",
    "        \n",
    "data = np.array(data)\n",
    "X = data[:, 0:-1]\n",
    "Y = data[:, -1]\n",
    "X = X.astype(int)\n",
    "Y = Y.astype(int)\n",
    "\n",
    "# normalize data\n",
    "x_zscore = scipy.stats.mstats.zscore(X, axis=0)\n",
    "Y_mean = np.mean(Y)\n",
    "Y_std = np.std(Y)\n",
    "y_zscore = scipy.stats.zscore(Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inputs matrix size: (100, 262143)\n",
      "Outputs matrix size: (100,)\n",
      "\n",
      "\n",
      "\n",
      "Completed at: 1552161200.9483578\n"
     ]
    }
   ],
   "source": [
    "# Print data shape\n",
    "print(\"Inputs matrix size: {}\".format(X.shape))\n",
    "print(\"Outputs matrix size: {}\".format(Y.shape))\n",
    "print(\"\\n\\n\")\n",
    "\n",
    "print(\"Completed at: {}\".format(time.time()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EXAMINING DATA\n",
    "Let's take a look at the data we are working with. Based on our earlier results we have chosen to use 2 PCs. Let's build our scores and loadings plot so we can better examine our data set in connection to these principle components. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PLSRegression(copy=True, max_iter=500, n_components=2, scale=True, tol=1e-06)"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plsda = PLSRegression(n_components=2)\n",
    "plsda.fit(x_zscore,y_zscore)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcgAAAG5CAYAAAAK36SOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XucHGWd7/Hvb0hIZlnBAIEIgSRAdoVAjO6AIEauEvClIqyuKCusN46uSjhm2QUviKyKGuMSWAVROaBHQXgpkqOsIKASFIRwzImEyxoyGRhwIBBALkPIkN/5o6qSTufpe1VXdffn/Xr1a7qrqruemknmO8+lnsfcXQAAYEt9eRcAAIAiIiABAAggIAEACCAgAQAIICABAAggIAEACCAggQ5jZmvM7Oi8ywF0OwISPcHM3mhmvzOzZ8xsnZn91swOzLtclZiZm9nzZvacmT1iZl83s20a/IzDzWy4xjHnmtmG+DzJY6+S/Zea2QNmttHM/qnsvaea2d1m9hczGzazr5rZuIYutHK51pjZaEmZbizZt7+Z3WBmT5iZl71vgpl918yGzOxZM/uDmR2XRpnQewhIdD0z217SzyRdJGlHSbtL+ryk9Smfp6EAq8Nr3P2vJR0l6b2SPpzy5yd+5O5/XfJYXbLv/0n6Z0n/N/C+v5J0hqSdJb0+Lue/pFiut5WU6ZiS7RskXS3pg4H3jJP0sKTDJO0g6bOSrjaz6SmWCz2CgEQv+BtJcvcr3f1ldx919xvdfUVygJl92Mzui2sd95rZ6+Lt+5rZr83saTNbaWZvL3nP5WZ2sZldb2bPSzoirsF8zcweMrPHzOwSM+uPj9/ZzH4Wf9Y6M1tqZjX/D7r7/ZKWStq/fF98vgvM7NH4cUG8bTtJ/yVpt5Ja2G6NfuPc/RvufrOkFwP7Lnb3pe7+krs/IukHkg4NfY6ZvSGu8e0Rv35N/H14dRNlesDdvytpZWDf8+5+rruvcfeN7v4zSYOS/q7R8wAEJHrBf0t62cyuMLPjzGxS6U4ze5ekcyWdIml7SW+X9KSZjZf0fyTdKGkXSZ+Q9AMz+9uSt79X0hclvULSbZK+oiiQ50jaR1Ft9Zz42AWShiVNlrSrpE9JqjnXo5ntJ2mupD8Edn9a0sHx+V4j6SBJn3H35yUdJ+nRklrYoxVO8bY4sFea2UdrlaeKNykQWpLk7r+T9C1JV8R/MHw/Luf9VT7vB2a21sxuNLPXNFMgM9tV0c8jWC6gKnfnwaPrH5L2lXS5ooAak7RE0q7xvhskzQ+8Z66kEUl9JduulHRu/PxySd8r2WeSnpe0d8m2QyQNxs/Pk3SdpH3qKK9L+oukpyQ9KOkLSTkkrZF0dPz8QUlvKXnfPElr4ueHSxqucZ79JO0maRtJb5D0Z0nvCRx3m6R/qvI574+/tztXOWa8pLsl/VHSLyRZlWMPldSvqBn37Pjn8MqyY/aJfoVVPd9Nkr6V978/Hp35oAaJnuDu97n7P7n7VEVNlbtJuiDevYeioCm3m6SH3X1jybYhRbXCxMMlzycr+oV+d9x8+LSiIJgc718oaZWkG81stZmdVaPYr3P3Se6+t7t/pqwcpWUcKitf3U2p7n6vuz/qUdPz7yQtlvTOet8vSWb2DklflnScuz9R5VwbFP1Rsb+kRe5esfbs7r/1qCn8BXc/X9LTiv5gqbdMfYpqqS9J+ni97wNKEZDoOR41612uzX16D0vaO3Doo5L2KOsn3FPSI6UfV/L8CUmjkma5+yvjxw4eDbSRuz/r7gvcfS9Jb5P0STM7qsXLeVTStLLyJU2pzSzV44pqwnUxs2MlfVvRgJo/1jh2d0mfk/S/JC0yswlZlMvMTNJ3FTVj/30czEDDCEh0PTN7tZktMLOp8es9JL1H0h3xId+R9C9m9ncW2cfMpkn6vaIm0381s/FmdriiYLsqdJ64hvdtSf9hZrvE59rdzObFz98af7Ypaj59OX604kpJnzGzyWa2s6L+zv8d73tM0k5mtkOlN5vZ8WY2Kb7ugySdrqgZONm/rZlNVBRO481sYvIHg5kdqWhgzt+7+53VChlf8+WKguuDippy/73CsXua2aHJuc3sTEUjZX+bfFZcpm3j1xPLwvZiRU3qb3P30WrlAqrKu42XB4+sH4qaRK9WVPN7Pv76LUnblxzzEUkPSHpO0j2SXhtvnyXpN5KekXSvpBNK3nO5pC+UnWuipC9JWq0oBO+TdHq8738q6j98XlF/3WerlNlVoa9SW/ZBTpR0oaLA+XP8fGLJsZdJelJRE+Vugc+6Mt7/nKT7k7KW7P91XJbSx+Hxvl8p6s99ruTxXxXKPF/SCknbxq93k7RW0tzAsbPiY5+Py3azpIGS/dMDZVoT75sWv36xrFwn5/3vkEfnPcydBZMBAChHEysAAAEEJAAAAQQkAAABBCQAAAGpzLxfVDvvvLNPnz4972IAAArk7rvvfsLdJ9c6rqsDcvr06Vq2bFnexQAAFIiZDdU+iiZWAACCCEgAAAIISAAAArq6DxIAUNmGDRs0PDysF1/caj3srjBx4kRNnTpV48ePb+r9BCQA9Kjh4WG94hWv0PTp0xXNJ9893F1PPvmkhoeHNWPGjKY+gyZWAOhRL774onbaaaeuC0dJMjPttNNOLdWOCUgA6GHdGI6JVq+NgAQAIICABAA0Zt1g3iVoCwISAFC/pYukC+dEX7scAQkAqM/SRdKtC6Pnty5MJSTXrFmjV7/61frQhz6k/fffXyeffLJuuukmHXrooZo5c6buvPNOPf/88/rABz6gAw88UK997Wt13XXXSZJWrlypgw46SHPmzNHs2bP1pz/9qeXylOI2D6DbrBuUdmxuWDtQURKOG0aj1xtGN4fl3AUtffSqVat0zTXX6NJLL9WBBx6oH/7wh7rtttu0ZMkSfelLX9J+++2nI488UpdddpmefvppHXTQQTr66KN1ySWXaP78+Tr55JP10ksv6eWXX27xIrdEQALdZOki6ebzpKPOafmXFrBJeTgmUgrJGTNm6IADDpAkzZo1S0cddZTMTAcccIDWrFmj4eFhLVmyRF/72tckRbenPPTQQzrkkEP0xS9+UcPDwzrxxBM1c+bMpssQQhMr0C0yaP4CtG4w+qOrPBwTG0aj/S0M3JkwYcKm5319fZte9/X1aWxsTO6uH//4x1q+fLmWL1+uhx56SPvuu6/e+973asmSJerv79e8efN0yy23NF2GEAIS6AaVmr8ISbRqxxlRi8T4/vD+8f3R/gyb9efNm6eLLrpI7i5J+sMf/iBJWr16tfbaay+dfvrpevvb364VK1akel4CEuh0tZq/sgzJHhnu3/PmLpDedObWITm+P9qecXP+Zz/7WW3YsEGzZ8/W/vvvr89+9rOSpB/96Efaf//9NWfOHN1///065ZRTUj2vJYncjQYGBpwFk9HV1g1GQ+5rOX15+n/h09/Z8e677z7tu+++9b+h9I+xNoVjq0LXaGZ3u/tArfdSgwQ6WV7NX/R39qakJil1RDi2ilGsQKdLfkmVN7Nm9Rd+hsP90QHmLpBmndgTtxJRgwS6QXkfUT3h2Ez/YZ79nSiOHghHiYAEukcjzV/NTBfWhuH+QJEQkEA3mbsgGpBTKxyb6T8swHB/oJ0ISKDbVAuoVu+XzHm4P9BOBCTQK9LqP2ymvxNdZWh0fd5FaAsCEugFafcf9thwf2y2eGhEr7/jPi0eGmnbOdOehLxeBCTQC7LoP6ynvxNdZfHQiC5Y85gk6YI1j6USkslyV6eeeqpmz56td77znXrhhRc0ffp0nXfeeXrjG9+oa665RsuXL9fBBx+s2bNn64QTTtBTTz3V8rlrISCBXpFF/yEDcnpGEo6jG6PZ10Y3emoh+cADD+i0007TihUrtP322+ub3/ymJGnixIm67bbbdNJJJ+mUU07RV77yFa1YsUIHHHCAPv/5z7d83loISKCX0H+IJpSHYyKtkNxjjz106KGHSpL+8R//Ubfddpsk6d3vfrck6ZlnntHTTz+tww47TJJ06qmn6tZbb23pnPUgIIFeQ/8hGjA0ul7nrx7ZKhwToxtd568eaWngjpkFX2+33XZNf2YaCEigF9F/iDpN65+gs/eaov4+C+7v7zOdvdcUTeufENxfj4ceeki33367JOnKK6/UG9/4xi3277DDDpo0aZKWLl0qSfr+97+/qTaZJQIS6FX0H6JO86dN0RnTd90qJPv7TGdM31Xzp01p6fP33XdfXXHFFZo9e7bWrVunj370o1sdc8UVV+jMM8/U7NmztXz5cp1zzjktnbMeTFYOVLJukBABYkkIJn2RaYWjJPX19emSSy7ZYtuaNWu2eD1nzhzdcccdLZ+rEdQggZBm5ioFulxSk5SUWjgWGTVIoFz5XKVSOn111EjRBeZPm6J37DKppT7HUtOnT9c999yTymeljRokUKrVuUqrfS41UhSQe3h0ajVphWPWmrm2UrkGpJldZmaPm9k9Jdt2NLNfmtmf4q+T4u1mZhea2SozW2Fmr8uv5OhKWa112OzqGd2GZbAKZ+LEiXryySdbDpIicnc9+eSTmjhxYtOfkXcT6+WS/lPS90q2nSXpZnf/spmdFb/+N0nHSZoZP14v6eL4K9C6ZK7SSpK5ShtdSb1SjVTqrVssli6Kvn9HndNb111wU6dO1fDwsNauXZt3UTIxceJETZ06ten35xqQ7n6rmU0v23y8pMPj51dI+rWigDxe0vc8+lPnDjN7pZm9yt3/3J7Soqslc5WGapDS5hlnWgnHRK+FZFZ9umjZ+PHjNWMG/eKVFLEPctck9OKvu8Tbd5f0cMlxw/G2LZjZaWa2zMyWdetfRchImnOVpr16RqfKqk8XaIMiBmQloWkctmo4d/dL3X3A3QcmT57chmKhq6Q1V2kWq2d0mqz6dIE2KWJAPmZmr5Kk+Ovj8fZhSXuUHDdV0qNtLht6QVpzlWaxekanoAaNLlDEgFwi6dT4+amSrivZfko8mvVgSc/Q/4jMpDVXaa+unkENGl0g79s8rpR0u6S/NbNhM/ugpC9LerOZ/UnSm+PXknS9pNWSVkn6tqR/zqHI6CVp/fIu6uoZWdfekuvuG7/l9r7xxfo+ABXkPYr1PRV2HRU41iV9LNsSARmZu6DxW0SylNx28YZPSMd8IeOTlQ8V6L577tCditjECnSnIoVjcrvF7y6Svn9CtufZOLbl9o1jDNJBRyAggV4SGln64C3phySDdNAFCEigV1S67UJKPyQZpIMuQEACvaBWjU6KQvKGz6R3vl6+zQVdgYAEesGOM6IBObXcflHrzZ6lK5f06m0u6Ap5T1YOoF2O+YL02MqophjSzHyz5arNu3rzeYQjOgoBCfSS910b9TWWh2QaNbtaK5cU6TYXoA40sQK95n3XSnsfufl1FuGYKJ13lXBEhyEggV70vmulQ+I+yVbDkVs60KVoYgV61bwvSAd+MJ2a3d5HZtu3CeSAgAR6WauhlUxZVz7faqk9D2FgDjoSTawAmrN0kfSbr0TPN26ofNyDt9C8io5EDRJA45JBOWPrqx9H8yo6GDVIAI2pNmVdKSYFQIejBgmgfsmI1VrGTSQc0fGoQQKoX61JyPviv7kP+1fCER2PGiSAxiTBV97MmjSpMmMOugQBCaAxyUod0uaQpL8RXYgmVgD1C63UIdUXjtzqgQ5DDRLoVusG023qvOEz0rLvRM8bnYQ8mVDgqHOoZaJjUIMEulFpTS8N3z8hWiuyfKWOeiYhL18Cq5kyUftEDghIoNukEUilQstjSVuGZK2yhIK1XmmHPVAnAhLoJmkEUqkbP1N5EvLk8yut1FHPEli1pB32QAMISKBbVAuk33y18XBZuki66zvVjxnfH/UrljezprEEVtphDzSIgAS6Qa1AGnuxsTUZa31eYuBD4UE3tSYUqBSsiTRqn0CLCEigG9QKJCma5WblT9L7vL2PjNaUrCS5DaT8M0L3TJYGNwswoyAISKBbVAqkxMaxxmpf1T5v7yOl913b+GeEwrF8EE6rtU8gJQQkkJYi1GiSQKq0gHGjTZShgHvDJ+oLx/LPkDZPRZeoNAinkdonkBECEkhDkW5FmHVi9QWMG22iLA+4Y6o0q1b7jNOXR8+T71OtQTj11D6BDDGTDtCq8lqQlO8v8aSJ8jdfCS9o3MwixvXOmFPNyp9s/v786nxJHjX7lkpCMjln8n28+TzCEW1n7p53GTIzMDDgy5Yty7sY6Gah0ZZFqekUqWz1LrJc6vTlmwM57Wnz0NPM7G53H6h1HE2s6Hx59f0V/VaEojRRNhqOoUE4hCNyQECis+XV91e0WxEqnafRFTfSVu/9lImi1L4BEZDoZHlOQ1akWxFu+Ez1PxKSATJ5hM6OM6JbQirpGx/dnym1PxyLMOoYhUZAojMVYRqyItyKkKyyIVW//ryaKJcukh66PbxvfL90xNnSEZ+KXrczHIs06hiFxShWdJ5afX9S+37RJudJytPucCydSDyP60+EBtHU6nvc85DN5Wx1hGwjijbqGIVFDRKdpWh9f1I+/XytLEGVtlBtrJ6+xwdv2fxzanc4MgE66kBAorMUqe+vVDv7+VpZgiptlfqBi/hzKvqoYxQOAYnOU4S+v5B2/LJfNyj97qLaxx3yiezL0+hMOIk8fk5FbHlA4RGQ6ExFucev3dJYZSMN9dbGivJzKmKNFoVHQKJz5X2PX17SWGWjFY3Wxtrxc6qn5jd3gXTgh7Kp0VLz7EoEJDpbnvf45SkUku0IR6m52liWP6d6b9lYuihqnt7zkHRrtNwy0rUISHS+Xm0WK62ZNboEVVrnbqQ2Vv5zSqPWVe9kEaXHPXR7FJJSOuGY12QVyBwBCXSypGbWzBJUaZy72f7FNGpd9d6yETruodujgUxphCO3jHQtAhLodHnWoJvpX0yj1lXvIKFqxy37TvNhxi0jPYHlrgC0rt7lqJYukn7zVWnsxc3bGu0HXDcY1T5rOfVn0hVvrX1c6bJaaZ6/0c9F27DcFYD2qTccf3X+luEoNV7rqneQ0Iy52dzawS0jPYOABJC95NaQjRvC+xu9Ub/eQUJZTVZQpEkQkBkmKweQvZU/iZa12jgW3p8ESyO1rnonis9qQvk8J6pHW1CDBJCtZEBLpXDsG998sNQ7SCiryQp6dbKKHsEgHQDZadeAlnoHCdV7XFbnRyEwSAdA/moNaOkbl86Alnrfn1WIEY5diT5IANkq76tLjJsgHfZvNEuisAhIANljQAs6EE2sANojNKCl1flYWUUDGSIgAbRP6aoerc7HyioayBgBmaGh0fV5FwEojqS2t+OM1udjZRUNtAEBmZHFQyN6/R33afHQSN5FAfJXWttrdRUMVtFAmzBIJwOLh0b09TWPSZIuiL/OnzYlzyIB+Smt7f3qfEm+9aQBSchJ0qwTK982UWsVDYmBP0gNNciULR4a0cLBEa3fGE3AMLrRdcGax6hJojeVB9rGDZVn1EnmY63Ur5jM51oejuXvZ+AOUkJApigJx7GyyYkISfSkSrW9SvrGRxMHSOEmU1bRQJsRkCkZGl2v81dvHY6J0Y2u81ePMHAHvaFWba9c33ht0fRaqV+RVTTQRgRkSn76+FNVO3THmXT2XlM0rX9C28oE5KbmFHMltcXycEzUG5KEIzJCQKZgU+2xyjFjLr1jl0ltKxOQu2q1vSPOlo74VPS6nn7JwaXhz5YIR2SGgEzBtP4JOnuvKervs+D+caL2iB5VrbaXTBpQq19x7yOlK94arkkmkw4AGSAgUzJ/2hSdMX3XrUJynKQz95rCbR7oXdVqezvOqF7T3PMQ6aHbo9eVBu4AGeE+yBQlIXjBmsc0utE1waRPziAcAc1dUP3+xtBk5kk4lk8IUHo8kCFqkClLapIS4QhsoVZtr7SmWR6OCWbNQRuZe4X7ErrAwMCAL1u2LJdzD42up88RaMbg0qjPsZbTl9PEiqaY2d3uPlDrOGqQGSEcgSbNmMuEACgEAhJA8TAhAAqAgCywIs66U8QyoUsxIQByRkA2qF0BUcTlsiqVidBEZpgQADkiIBuQRWiFwmXx0MimZbKKMsl5pTIVMcjRZZgQADnpuIA0s2PN7AEzW2VmZ7XrvFmEVihckvOMFmi5rEplevfyVfr6YFSuWmWklomWMCAHOeioiQLMbBtJ35D0ZknDku4ysyXufm+W560UEFLzCyGXB26i9DyJNM7XrPJrLy3Tb556rq4yLh4a0fmrR3Q2MwoB6CAdFZCSDpK0yt1XS5KZXSXpeEmZBWS1gGg2tEKB+/XBEa2vcktqslzWO3aZ1LZbSJJJ2OsV+p6E/hAgJAF0gk5rYt1d0sMlr4fjbZuY2WlmtszMlq1du7alkyUBUR6OiWbWeKwUuOs9+mtlXHi+c/X3WaoTntdT5lqTsIeUfk+K2FwMAPXqtIAM/abeImnc/VJ3H3D3gcmTJ7d0sloB0Who1QrcMUXLYk0oO19/n+mM6bumVvNqZGBNpUnYK0lWLvnp409VrXkTkgCKrtMCcljSHiWvp0p6NMsTVgqIZkKr3sD9ZMn5sgjHRgcbJd+DCXVk5JikA7ffLvWaNwC0W6cF5F2SZprZDDPbVtJJkpZkfdLykGwltOoJ3NIJz7MIx2aaPOdPm6L3T925rvPsPnHbVGveAJCHjhqk4+5jZvZxSTdI2kbSZe6+sh3nTkLq/NUjLYdW+bJYocCdP21KqgNy0hhsdO4+U/WHv7yg3z/zQsVjTNJPH39qq2tMNPLHBRO+A8hTRwWkJLn79ZKuz+PcaYZWPYGb5oCcaqNRQyNkQ+E0NLq+ajhKUYdwcm9kPX8IVMKtIcjNukHuu4SkzmtizV2aNZr506bo9wfvm3kANDrYqNIgnnpHta73zf2bzTQXF3EmIfSIpYukC+ew3iQkdWANstu0owlxaHT9pmBaODiisZJW1vJaXa37FpPnX109opernLO0VtpIzTuLSRmAuixdFC3GLG3+yvR2PY0aZJfbqjZYEo7jTMFwrGcQj9UY0VpeK20mHBPcGoLMJeG4YTR6vWE0ek1NsqdRg+xipbXBhYMjkke3YWziWx9bKZye2jCmc/eZWvfsOjuN26ahGl8z/aRAKsrDMZGEpERNskdRg+xS5YE3Vh6Oil5fsOYxnbtquOZ9i5c8/IQ+v+oRTeufoNfv8Fc1zz/80pg+v+qRusub9qQMQF3WDUo3n7d1OCY2jEb71w22t1woBAIyQ3ndCF+pNhiShN9H95hcc/DN5Y88oXcvX6UVz1b4ZVLm4ofXNvQ9SHNSBqAuO86Qjjpn86LM5cb3R/sZ1dqTCMiM5LVOYq3p7MolNbPP7bN7zSnlkhU86v3swyb9dV3HlUpzUgagLsmizOUhOb6fRZp7HAGZgTxvU2hkgvHy8Gl03tVafvvUc039kZDVTEJAReUhSThCDNJJXRFuU6g0i80406aBOpVqZvOnTdFTG8Z0ycNPtFyOpM/z601cf9ozCQE1JWF483mEIyQRkKnKYu3IZlWaxUYKz95TOnNOPVPKNWJ9k9dPOKLt5i6QZp1InyMk0cSamizWjmxVqKmydPaepCzl/aX1TCnXKO5lRMcgHBEjIFNS1NsUQtPZTeufsCkU3718VbC/tNGFkuvBMlcAOgkBmaKi3qZQHsqlg4hKR6WObnQtHIyCU1KqA3akqA+0/I8EwhJAURGQKSv6bQq17pFM5mldGM9qc9AO26V38rJT5nUrDADUg4DMQFFvU2hkAoExRROS//bp57bYPs6iNR/Lt9Uz2mtM2tTE2sqtMNQ6AbQDAZmRdi1lVa/fPfVsQxMISNLL0hYrfyh+vY1Fq1VLUQ35zBlTdGYdfZZJP+xPH3+qrknRQ0FIrRNdjSntCoWAzFBRblNYPDSiE5c/qMMm/XUqfYpjrk1LXZWOjj1j+q4Va5Klt5nUs2JHKAhZJxJdjbUoC4f7IAum9H7ENJSGyp3PPK+DdthOdz7zfEM1yXJJ2AVv5I8nI+hT9CidlODA7bfTicsfrPi5ySjXpze8rMsfiSYqSMqePGedSHQl1qIsJGqQBZJ282FoVp8kJEsHESVzppb/tTRO8ew7JUoHHZWGY3KupEl2o7Sps/KM6bvq6Q0v16zFJmW5/JEnthpZuzDQPMy9legKrEVZWARkQaTdfFhtVp8kJKUovH40Zx/9/uB9t+hH7O8znbnXFJ05Y0rNEbmVzjXm0oQ+01WPPqmLH14rSVsFdKK/zyrWbkNLdZVeD/dWomPVWouSkMwVTawFkPb8rfUsPvybp57TT+bsrTdMeoWkqL80OVdoKrrQtnrOtX6ja/DFDVucu7ypt7/PdOruOzU1/2sS2kXp7wXqlqxFWUmyFiVT3+WGGmTOas3f2kxNst5ZfZJwLBUafVttRG4jq4ckQrXYc/eZWvVzQreSFO0eU6AhrEVZeARkjrKcv7WVWX1CtbFqNbRmlskqrcXWWm4rdCsJ4YiuwFqUhUZA5ijr+VvTntWnWlA3E5If3WPyVrXYamUu6gQMQEtYi7KwCMicZT1/a1qhUs8I29Jz1brn0iS9cvw2wX3Vyly0CRiAVCQhKRGOBWLuzd8PV3QDAwO+bNmyvItRl9K+yCyaD5u9v3JodP0WM99UvQey7FyLh0a0cHBkq9l4ErWuM+17QoHCWzdIn2MbmNnd7j5Q8zgCsjgWD43o/NUjOnuvKYWoISXlGWdbTjmXvK5VzqHR9ZtWBqnm9wfvSxACaJt6A5Im1gIpUvPh4qGRTSt6hOZjlaIVP6o1uSZ9rNWmn6vUx8p9jQDyRkAWTBFqUouHRvT1wZGKN+cnxiR9fbB6SM6fFo0+nVBlRp7Q+ZmQHEDeCEhsIekLXV9ny/t6rz3zz/xpU/TJOmbkKT2/xITkKLBmVt1gpY6OQ0Bik1r3ZVZSz/2a9YymrTSjECGJQmlm1Y08VuogkFtGQGKTZmbFkaIZbuq5X7NaH2sWMwoBqStfdaOewGvmPa1i6axUEJDYQqX7MqtGZgN5WmlATlYzCgGpaWbVjXrek3ZNL49A7lIEJLZSHpITTKrW6DrmainAsp5RCGhZM6tu1POetGt6LJ2VKlbzQFDpyh6dfr0iAAAR4UlEQVSfnBE9DzWBSumsqJGcr/wczLmK3DWz6ka97xk3MXqdxiLJtQK51c/vQdQgUVFpn2HWU+Il50tz7lggFc2sulHrPUm/xNiL0ddWa3pJIJeHYyIJZAbuNISARFWltcJ2BBgTkqOQmll1o9J7ZAp2WrQSkiydlQmaWNGQaosqp3mOanO9ArlIQjBpxqxn1Y3y99g2kr9c+fhWFkkuP1eC1UGaRkD2oFYnAW9HgBGOKKQkZG4+r/7QKX2P9VUPyCTMmq3pNRPiqIgm1h6T1jRuBBh61twF0unLGwudWSdGXzduqH7cwIdaDzOWzkoNAdlDmMYNSEmjNbyag3Yk7X2kNO8LrZUr0UyIYysEZI9gGjcgZxUH7SgKx/ddW/39jY5AZUBOywjIHsA0bkBBlIfk+H7pDZ+oHY5MHZcLArLLMY0bUDDlfYTH1GhWZeq43DCKtcsl07hlOQsOgAbNXVDfrRyVpo5LPgOZogbZA9oxCw6ABjUajgnmV20bArJHMI0b0EGYOq4QCMgewjRuQIdg6rhCoA+yxzCNG9AhmDoud9QgexDhCHSI0G0hhGPbEJAAUGRMHZcbmlgBoOjqvS0EqaIGCQBFVTpKlXBsOwISAIqI6eVyR0ACQNEwvVwh1AxIM9vezPYObJ+dTZEAoIdVml6OkGy7qgFpZv8g6X5JPzazlWZ2YMnuy7MsGAD0HKaXK5RaNchPSfo7d58j6f2Svm9m8dLYsspvAwA0hOnlCqdWQG7j7n+WJHe/U9IRkj5tZqdLCq+fBABoHNPLFU6tgHy2tP8xDsvDJR0vaVaG5QKA3lM+c06CGXRyUWuigI+qrCnV3Z81s2Ml/UNmpQKAXlU+ByvhmJtaNcjnJe0a2H6wpDvSLw4AgOnliqFWQF4g6dnA9tF4HwAgC3MXSKcvJxxzVCsgp7v7ivKN7r5M0vRMSgQAiDAgJ1e1AnJilX0VhloBAND5agXkXWb24fKNZvZBSXdnUyQAQMO4PzJ1tUaxniHpWjM7WZsDcUDStpJOyLJgAIA6LV0UTSJw1Dn0WaaoakC6+2OS3mBmR0jaP978c3e/JfOSAQBqK5/YXCIkU1I1IM1soqSPSNpH0h8lfdfdx9pRMABADZUmNpcIyRTU6oO8QlGT6h8lHSfpa5mXCABQGxObZ65WH+R+7n6AJJnZdyXdmX2RAABVJRObV5JMbD7rRG4VaUGtGuSG5AlNqwBQEExs3ha1apCvMbO/xM9NUn/82iS5u2+faekAAGHlc7YmmLs1NbVGsW7TroIAABrExOaZqtXEmgkze5eZrTSzjWY2ULbvbDNbZWYPmNm8ku3HxttWmdlZ7S81ABQQE5tnplYTa1bukXSipG+VbjSz/SSdpGityd0k3WRmfxPv/oakN0saVjTDzxJ3v7d9RQaAgpq7gAE5GcglIN39Pkkys/Jdx0u6yt3XSxo0s1WSDor3rXL31fH7roqPJSABQCIcM5BLE2sVu0t6uOT1cLyt0vatmNlpZrbMzJatXbs2s4ICALpbZjVIM7tJ0pTArk+7+3WV3hbY5goHuYc+wN0vlXSpJA0MDASPAQCglswC0t2PbuJtw5L2KHk9VdKj8fNK2wEASF3RmliXSDrJzCaY2QxJMxXN3nOXpJlmNsPMtlU0kGdJjuUEAHS5XAbpmNkJki6SNFnSz81subvPc/eVZna1osE3Y5I+5u4vx+/5uKQbJG0j6TJ3X5lH2QEAvcHcu7ebbmBgwJctW5Z3MQAgmj+VkaaFYGZ3u/tAreOK1sQKAN1n6SLpwjmssNFhCEgAyFL5gsaEZMcgIAEgK5UWNCYkOwIBCQBZYEHjjkdAAkDakgWNy8MxkSxovG6wveVCQwhIAEgbCxp3BQISALKQLENVHpKs2dgxCEgAyEp5SBKOHYWABIAssaBxx8prwWQA6B0saNyRqEECQDsQjh2HgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACCAgAQAIICABAAggIAEACAgl4A0s4Vmdr+ZrTCza83slSX7zjazVWb2gJnNK9l+bLxtlZmdlUe5AQC9I68a5C8l7e/usyX9t6SzJcnM9pN0kqRZko6V9E0z28bMtpH0DUnHSdpP0nviYwEAvWLdYFtPl0tAuvuN7j4Wv7xD0tT4+fGSrnL39e4+KGmVpIPixyp3X+3uL0m6Kj4WANALli6SLpwTfW2TIvRBfkDSf8XPd5f0cMm+4Xhbpe1bMbPTzGyZmS1bu3ZtBsUFALTV0kXSrQuj57cubFtIjsvqg83sJklTArs+7e7Xxcd8WtKYpB8kbwsc7woHuYfO6+6XSrpUkgYGBoLHAAA6RBKOG0aj1xtGN4fl3AWZnjqzgHT3o6vtN7NTJb1V0lHungTZsKQ9Sg6bKunR+Hml7QCAblQejok2hWReo1iPlfRvkt7u7i+U7Foi6SQzm2BmMyTNlHSnpLskzTSzGWa2raKBPEvaXW4AQJusG5RuPm/rcExsGI32ZzhwJ68+yP+U9ApJvzSz5WZ2iSS5+0pJV0u6V9IvJH3M3V+OB/R8XNINku6TdHV8LACgG+04QzrqHGl8f3j/+P5o/44zMiuCbW7d7D4DAwO+bNmyvIsBAGhWqJl1fL/0pjObbl41s7vdfaDWcUUYxQoAQNjcBVEYJjXJFsOxEQQkAKDYkpCU2haOUoajWAEASM3cBdKsEzPtcyxHDRIA0BnaGI4SAQkAQBABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQBAAAEJAEAAAQkAQAABCQDIxrrBvEvQEgISAJC+pYukC+dEXzsUAQkASNfSRdKtC6Pnty7s2JAkIAEA6UnCccNo9HrDaMeGJAEJAEhHeTgmOjQkCUgAQOvWDUo3n7d1OCY2jEb7O2jgDgEJAGjdjjOko86RxveH94/vj/bvOKO95WoBAQkASMfcBdKbztw6JMf3R9vnLsinXE0iIAEA6SkPyQ4NR4mABACkLQlJqWPDUZLG5V0AAEAXmrtAmnViR/U5lqMGCQDIRgeHo0RAAgAQREACABBAQAIAEEBAAgAQQEACABBAQAIAEEBAAgAQQEACABBAQAIAEEBAAgAQQEACABBAQAIAEEBAAgAQQEACTRoaXZ93EQBkiIAEmrB4aESvv+M+LR4aybsoADJCQAINWjw0ogvWPCZJumDNY4Qk0KUISKABSTiObnRJ0uhGJySBLkVAAnUqD8cEIQl0JwISqMPQ6Hqdv3pkq3BMjG50nb96hIE7QBchIIE6TOufoLP3mqL+Pgvu7+8znb3XFE3rn9DmkgHICgEJ1Gn+tCk6Y/quW4Vkf5/pjOm7av60KTmVDEAWCEigAeUhSTgC3YuABBqUhKQkwhHoYuPyLgDQieZPm6J37DKJPkegi1GDBJpEOALdjYAEACCAgAQAICCXgDSzfzezFWa23MxuNLPd4u1mZhea2ap4/+tK3nOqmf0pfpyaR7kBAL0jrxrkQnef7e5zJP1M0jnx9uMkzYwfp0m6WJLMbEdJn5P0ekkHSfqcmU1qe6kBAD0jl4B097+UvNxOUjJ/1/GSvueROyS90sxeJWmepF+6+zp3f0rSLyUd29ZCAwB6Sm63eZjZFyWdIukZSUfEm3eX9HDJYcPxtkrbQ597mqLap/bcc890Cw0A6BmZ1SDN7CYzuyfwOF6S3P3T7r6HpB9I+njytsBHeZXtW290v9TdB9x9YPLkyWlcCgCgB2VWg3T3o+s89IeSfq6oj3FY0h4l+6ZKejTefnjZ9l+3XEgAACrIaxTrzJKXb5d0f/x8iaRT4tGsB0t6xt3/LOkGSceY2aR4cM4x8TYAADKRVx/kl83sbyVtlDQk6SPx9uslvUXSKkkvSHq/JLn7OjP7d0l3xced5+7r2ltkAEAvySUg3f3vK2x3SR+rsO8ySZdlWS4AABLMpAMAQAABCQBAAAEJAECARd1+3cnM1ioaBNRpdpb0RN6FaINeuM5euEaJ6+w23X6d09y95o3yXR2QncrMlrn7QN7lyFovXGcvXKPEdXabXrnOWmhiBQAggIAEACCAgCymS/MuQJv0wnX2wjVKXGe36ZXrrIo+SAAAAqhBAgAQQEACABBAQObIzN5lZivNbKOZDZTtO9vMVpnZA2Y2r2T7sfG2VWZ2VvtL3RozO9fMHjGz5fHjLSX7gtfcqTr9Z1WNma0xsz/GP8Nl8bYdzeyXZvan+OukvMvZKDO7zMweN7N7SrYFrytedejC+Oe7wsxel1/JG1PhOnvm/2a9CMh83SPpREm3lm40s/0knSRplqRjJX3TzLYxs20kfUPScZL2k/Se+NhO8x/uPid+XC9VvuY8C9mKLvpZVXNE/DNM/rg7S9LN7j5T0s3x605zuaJ/f6UqXddxkmbGj9MkXdymMqbhcm19nVIP/N9sBAGZI3e/z90fCOw6XtJV7r7e3QcVLf91UPxY5e6r3f0lSVfFx3aDStfcqbr5Z1XJ8ZKuiJ9fIekdOZalKe5+q6TypfQqXdfxkr7nkTskvdLMXtWekramwnVW0m3/N+tGQBbT7pIeLnk9HG+rtL3TfDxukrqspBmuW64t0W3XU84l3Whmd5vZafG2XeMFzhV/3SW30qWr0nV148+4F/5v1o2AzJiZ3WRm9wQe1WoTFtjmVbYXSo1rvljS3pLmSPqzpEXJ2wIfVbhra0C3XU+5Q939dYqaGT9mZm/Ku0A56Lafca/836xbLgsm9xJ3P7qJtw1L2qPk9VRJj8bPK20vjHqv2cy+Leln8ctq19yJuu16tuDuj8ZfHzezaxU1uT1mZq9y9z/HTY2P51rI9FS6rq76Gbv7Y8nzLv+/WTdqkMW0RNJJZjbBzGYoGgRwp6S7JM00sxlmtq2ijvMlOZazYWV9NCcoGqgkVb7mTtXxP6tKzGw7M3tF8lzSMYp+jksknRofdqqk6/IpYeoqXdcSSafEo1kPlvRM0hTbiXro/2bdqEHmyMxOkHSRpMmSfm5my919nruvNLOrJd0raUzSx9z95fg9H5d0g6RtJF3m7itzKn6zvmpmcxQ10ayR9D8kqdo1dyJ3H+uCn1Ulu0q61syk6HfID939F2Z2l6SrzeyDkh6S9K4cy9gUM7tS0uGSdjazYUmfk/Rlha/reklvUTRo5QVJ7297gZtU4ToP74X/m41gqjkAAAJoYgUAIICABAAggIAEACCAgAQAIICABAAggIAEuoCZvRyvwHCPmV1jZn8Vb59iZleZ2YNmdq+ZXW9mfxPv+4WZPW1mP6v+6UBvIiCB7jAar8Cwv6SXJH3EohsVr5X0a3ff2933k/QpRfcxStJCSe/Lp7hA8RGQQPdZKmkfSUdI2uDulyQ73H25uy+Nn98s6dl8iggUHwEJdBEzG6doAvE/Stpf0t35lgjoXAQk0B36zWy5pGWKpkP7bs7lAToec7EC3WHU3eeUbjCzlZLemVN5gI5HDRLoXrdImmBmH042mNmBZnZYjmUCOgYBCXQpj1YiOEHSm+PbPFZKOlfxWn5mtlTSNZKOMrNhM5uXW2GBAmI1DwAAAqhBAgAQQEACABBAQAIAEEBAAgAQQEACABBAQAIAEEBAAgAQ8P8Bmv5he4VW8IcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "scores_PC1 = plsda.x_scores_[:, 0]\n",
    "scores_PC2 = plsda.x_scores_[:, 1]\n",
    "\n",
    "plt.figure(figsize =(7,7))\n",
    "plt.title(\"Scores Plot 512 x 512\")\n",
    "plt.scatter(scores_PC1[0:50], scores_PC2[0:50], c='tab:orange',marker='D', label=\"mes\")\n",
    "plt.scatter(scores_PC1[50:], scores_PC2[50:], c='tab:cyan',marker='D', label=\"pro\")\n",
    "plt.xlabel(\"PC1\")\n",
    "plt.ylabel(\"PC2\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x_loadings_PC1 = plsda.x_loadings_[:, 0]\n",
    "# x_loadings_PC2 = plsda.x_loadings_[:, 1] \n",
    "\n",
    "# y_loadings_PC1 = plsda.y_loadings_[:, 0]\n",
    "# y_loadings_PC2 = plsda.y_loadings_[:, 1] \n",
    "\n",
    "# y_weights_PC1  = plsda.y_weights_[:,0]\n",
    "# y_weights_PC2  = plsda.y_weights_[:,1]\n",
    "\n",
    "# # we are looking to plot w*c\n",
    "# # first lets find w*, where w* = w(p'w)^-1\n",
    "\n",
    "# wStar =(np.matmul(plsda.x_weights_, np.linalg.inv(np.matmul(plsda.x_loadings_.T, plsda.x_weights_))))\n",
    "\n",
    "# x_wStar_PC1  = wStar[:,0]\n",
    "# x_wStar_PC2  = wStar[:,1]\n",
    "\n",
    "# plt.figure(figsize=(4, 4), dpi=100)\n",
    "\n",
    "# # need to scale \n",
    "# plt.scatter(0.15*y_loadings_PC1,0.15*y_loadings_PC2, c='b', label = \"Response\") \n",
    "# plt.scatter(np.multiply(x_wStar_PC1, y_weights_PC1), np.multiply(x_wStar_PC2, y_weights_PC2), c='y', label=\"Signals\")\n",
    "# plt.legend()\n",
    "# plt.title(\"PC1 vs PC2 -- Loadings Plot\")\n",
    "# plt.xlabel(\"PC1\")\n",
    "# plt.ylabel(\"PC2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PERFORMING PLSDA\n",
    "Now that we have processed our csv file and correctly stored our oberservations into an array, we will begin to construct our model using PLSDA. After doing some research online we discovered a method for implementing PLSDA via PLSR. We will build a model using PLSR, where our yobs is an array of 0s (mesanchymal) or 1s (proneural). PLSR will make predictions which will result in continous values from 0 to 1. However, we need our data to be discrete so that we can correctly classify each image into a cell type. Thus, following PLSR we will need to threshold our ypred."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculation of R2X and R2Y\n",
    "This will be used to determin how many primary components to use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "92.81525762243389 % of variation for y accounted for using 2 PCs.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xu8XPO9//HXWxKSkiYIPUg0qKpLE5wdtBRpEFTj0rprq269HKqc9iEaR1PqVPU4PfXjtLQ0tJT0lsatigin6pJoCBKXVIsdRJrIBQmJfH5/rDVjMmbPXnvPXjN79n4/H4/92LPun1mzZj5rfb9rfb+KCMzMzADWaXQAZmbWfTgpmJlZkZOCmZkVOSmYmVmRk4KZmRU5KZiZWZGTQpOTNEnSdxu0bUn6uaTXJD3ciBiaiaTbJX2h0XH0VpJOlPTnRsfR3TkpdDFJ/5C0QNL6JeNOkTS9gWHlZS9gf2BoROxWPjH9Er4j6XVJyyQ9KumQkunvl/Q/kl5I55mXDg+p55uoEPex6eeosvF9Jb1a+h46IiIOiohruybKrpN+TiHpqEbHUo2k6ZJOyWndw9N98HrJ32N5bKu7c1LIR1/gzEYH0VGS+nRwkQ8C/4iIN6rM80BEbAAMBq4GJkvaSNK6wN3AjsCBwPuBjwOLgPckmDr7PUm8+5SNPxAI4I8dWVl6RdWdv2tfABan/3u7wRGxQfo3sqMLS+qbR1D11J0P1Gb2A+AbkgaXTyg5I+lbMq54BpSetd0v6YeSlkh6TtLH0/Evpmeq5V/eIZLulLRc0r2SPliy7o+k0xZLerr0bDAtevqxpNskvQGMrhDv5pKmpsvPk3RqOv5k4GfAx9Kzqu9U2yERsQa4BhgAbA18HtgSODwi5kTEmoh4NSIujIjb0m2cI2l++r6eljSmQnx7SHqlNKFJOlzS7PT1bpJmplcqCyT9d7U401hXApPTGEt9Hrg+IlZL2lDSLZIWpsVnt0gaWhLDdEkXSbofeBPYuuxz3kbSNEmLJP1T0vWlx0t6pfINSbMlLZV0k6T+JdMPTa+8lkn6m6QD0/GDJF0t6eV03323WrJPj5V9gNOAsZI+UDa9re1spKTo8KX0/U8pWebU9FhZnB47m6fjsxz7f5b0X+k6/y7poHTaRcAngMvT4+3ydHy143vjdPvLlBRvbtPWfqhG0jqSzpP0fPr9u07SoLL3dLKkF4BpndlGtxIR/uvCP+AfwH7A74DvpuNOAaanr4eTnG32LVlmOnBK+vpEYDXwRaAP8F3gBeAKYD3gAGA5sEE6/6R0eO90+o+AP6fT1gdeTNfVF9gV+CewY8myS4E9SU4Q+ld4P/cC/wv0B3YGFgJjSmL9c5V9cWJJLIWrp+XAIOBG4Noqy26Xxr55yX7bpo15/wbsXzL8a2B8+voB4HPp6w2APTJ+jnsCy4AB6fAgYAWwczq8MfAZ4H3AwHSbU8o+0xdIroT6Av3KPucPkRS9rQdsAtwH/E/ZcfQwsDmwETAX+HI6bbf0c9s//dy2AD6STpsCXJl+9pum6/hSlff5H8DD6evHgbNLplXbzq3ATcCG6XvbJx3/yfQY2zV9b/8PuK8Dx/4q4FSSY/8rwEuAyufNeHzfSJLc1wd2AubTxvFaKbaSaScB80hOZjYg+W7/omy569LtDGj0b1Ctfw0PoKf98W5S2Cn9Qm1Cx5PCsyXTPprO/4GScYt498dpEnBjybQNgHeAYcDRwP+VxXcl8O2SZa+r8l6GpesaWDLue8CkkljbSwqrgSXpl/VBYL902p3AxVWW/RDwarov+7Wzz78LXJO+Hgi8AXwwHb4P+A4wpBOf5bPAcenrU4HHqsy7M/Ba2Wd6Qdk8xc+5wvKHAbPKjqMTSoYvAX5S8hn+sMI6PgC8VfrDBBwL3NPOe/x6+vrc0vdYZTubAWuADStMuxq4pOx4XEVy3A+n/WN/Xsm096Xz/0ul/Vft+CZJKqtIk1g67T/bOl5LYltS8veNdNrdwFdL5t0uXXffkuW27ujx1V3/XHyUk4h4ArgFGN+JxReUvF6Rrq983AYlwy+WbPd1kvLhzUnK/HdXUgy1RNIS4HjgXyotW8HmwOKIWF4y7nmSM8asHoyIwRExJCL2iIi70vGLSH5cKoqIecDXgYnAq5JuLBRDVHADcISk9YAjgL9GxPPptJOBDwNPSZqhjlUSX8e7RUifA4qVxJLeJ+nKtEhhGUnyGVxWVNPmvpW0afqe5qfL/xIor2B/peT1m7z7mQ8juToq90GSs/aXSz7vK0muGCrFsCewFckZNST78aOSdm5nO8NIjovXKkzbnOQYAYrH4yKyHzPF9xwRb6YvN2hj3mrH9yYkP9qln8HzFdZRbkh6vA6OiP9Kx631ntLXfUmScEG171FTcVLI17dJzjBLvxCFStn3lYwr/ZHujGGFF5I2IClueInkQL235CAvVKJ9pWTZas3kvgRsJGlgybgtSS7Da3UXSRn2+m3NEBE3RMReJF/+AL7fxnxzSL6oBwHHkfy4FaY9GxHHkvwwfh/4TbVtlrkOGCPpY8AepesF/p3kjHH3iHg/SfEdQOkdS9X27ffS6SPS5U8oW7aaF6lcPv4iyZVC6Q/b+yNixzbW84V0m49KegV4KB1fSITVtrORKtSZkRwzpXVa65MUtc2n9mO/fH9WO74XklylDiuZf8sObKvUWu8pXc9q1j556zHNTTsp5Cg9270J+FrJuIUkX5ATJPWRdBKdrAArcbCkvZTc0XMh8FBEvEhypfJhSZ+T1C/9GyVp+4zxvwj8BfiepP6SRpCceV9fY7wAvyD5Uv82rSxcJ60Y/JakgyVtJ+mT6dn/SpKro3eqrO8Gkv28N0n5PgCSTpC0SSQV3UvS0dXWU5RebfwZ+BVwZ0SUnrkPTGNaImkjkhOAjhgIvJ4uvwXwzQ4sezXwRUlj0v22haSPRMTLwJ+AS5Xc7ruOkgrt8ruoSCutjyKpYN655O8M4Pi0Mrjadm4H/ldJhXs/SYWkeEO6zM7pZ/efJMfjP7rg2F9AUq5f0ObxHRHvkJT9T0yv6nag83dX/Qo4S9JW6UnXfwI3RcTqTq6vW3NSyN8FJBVQpU4l+RFYRFIR+Zcat3EDyY/SYuBfSS6hSYt9DgCOITnbeYXkbHm9Dqz7WJJy05dIbtX8dkTcWWO8RMRbJPUFT5HULywjqRQdQnLGuh5wMUldxCskZ/rfqrLKXwH7AtMi4p8l4w8EnpT0Okkl/DGR3F1EehfLJ9oJ9VqSs8Trysb/D8mdVIW6kg7dpkpSz7ErSb3TrSQ/YJlExMMklas/TJe/l3fPZD8PrAvMAV4DfkPlYrrDSJLadRHxSuGPJBH0AQ5sZzufIylXf4qk7ufraWx3k1Re/xZ4meRH/5iS7dZy7P8I+Gx6Z9JlGY7v00mKnl4hqT/7eQe2VeoakpOY+4C/k5yknNHJdXV7hVp9MzMzXymYmdm7nBTMzKzIScHMzIqcFMzMrKjpGm8aMmRIDB8+vNFhmJk1lUceeeSfEbFJe/M1XVIYPnw4M2fObHQYZmZNRVKWJ7pdfGRmZu9yUjAzsyInBTMzK3JSMDOzIicFMzMrclIwM7Oi3JKCpGvS/kyfaGO6JF2mpC/X2ZJ2zSsW6wFmT4Yf7gQTByf/Z09unuUde3Mu38yx1yDP5xQmAZfz3iaHCw4Ctk3/dgd+nP63vMyeDHdfAEtbYdBQGHM+jDiq/eW6Yvlal735a7BqRTK89MVkGLKto5HLO/bmXL6ZY69RblcKEXEfSfv+bTmUpC33iIgHSboybLN7RqtR4SBb+iIQ7x5kWc8+alm+1m3ffcG7X46CVSuS8d19ecfenMs3c+w1amSdwhas3a9pK2304yrpNEkzJc1cuHBhXYLrlmq5nGzmL9jS1o6N707LO/bmXL6ZY69RI5NCpf5oK/b4ExFXRURLRLRsskm7TXf0TLWebTfzF2zQ0I6N707LO/bmXL6ZY69RI5NCK2t3qj2UpEs9q6TWs+1m/oKNOR/6DVh7XL8Byfjuvrxjb87lmzn2GjUyKUwFPp/ehbQHsDTtELxn62wRUK1n2838BRtxFHz6Mhg0DFDy/9OXZa9wa+Tyjr05l2/m2GuUWx/NkgodqQ8BFpB0LN8PICJ+IkkkdycdCLwJfDEi2m3+tKWlJZq2ldTyOwog+XHM8mH/cKe06KjMoGFwVsW7fitvvxnvPjKzmkl6JCJa2p0vr6SQl6ZOCrX8sNeSUMys18uaFPxEcz3VUgTUwMtJM+s9mq6TnaY2aGgbVwoZK1xHHOUkYGa58pVCPTXwjgIzsyycFOrJRUBm1s25+KjeXARkZt2YrxTMzKzISaGjGtScrZlZPbj4qCMa2JytmVk9+EqhIxrYnK2ZWT04KXREA5uzNTOrByeFjmhgc7ZmZvXgpNARfvjMzHo4J4WO8MNnZtbD+e6jjvLDZ2bWg/lKwczMipwUzMysyEnBzMyKnBTMzKzIScHMzIqcFMzMrMhJwczMipwUzMysyEnBzMyKnBTMzKzIScHMzIqcFMzMrMhJwczMipwUzMysyEnBzMyKnBTMzKzIScHMzIpyTQqSDpT0tKR5ksZXmL6lpHskzZI0W9LBecZjZmbV5ZYUJPUBrgAOAnYAjpW0Q9ls5wGTI2IX4Bjgf/OKp2j2ZPjhTjBxcPJ/9uTcN2lm1izyvFLYDZgXEc9FxNvAjcChZfME8P709SDgpRzjSRLAzV+DpS8mm176YjLsxGBmBrSTFCStI+njnVz3FsCLJcOt6bhSE4ETJLUCtwFntBHHaZJmSpq5cOHCToYD3H0BrFqx9rhVK5LxZmZWPSlExBrg0k6uW5VWWTZ8LDApIoYCBwO/kPSemCLiqohoiYiWTTbZpJPhAEtbOzbezKyXyVJ89CdJn5FU6Ue+mlZgWMnwUN5bPHQyMBkgIh4A+gNDOrid7AYN7dh4M7NeJktSOBv4NfC2pGWSlktalmG5GcC2kraStC5JRfLUsnleAMYASNqeJCnUUD7UjjHnQ78Ba4/rNyAZb2Zm9G1vhogY2JkVR8RqSacDdwB9gGsi4klJFwAzI2Iq8O/ATyWdRVK0dGJElBcxdZ0RRyX/774gKTIaNDRJCIXxZma9nLL8BksaB+ydDk6PiFtyjaqKlpaWmDlzZqM2b2bWlCQ9EhEt7c3XbvGRpIuBM4E56d+Z6TgzM+th2i0+IrkraOf0TiQkXQvMAt7zhLKZmTW3rA+vDS55PSiPQMzMrPGyXCl8D5gl6R6SZw/2Bs7NNSozM2uILHcf/UrSdGAUSVI4JyJeyTswMzOrvzaLj9LbSQs2ioipEfEHJwQzs56rWp3CSSWvf5F3IGZm1nhZK5o72sSFmZk1oWp1CoMlHU6SON4v6YjSiRHxu1wjMzOzuquWFO4FxqWv7wM+XTItACcFM7Meps2kEBFfrGcgZmbWeLn20WxmZs3FScHMzIqcFMzMrKjNOoXyu43K+e4jM7Oep9rdR4W7jTYFPg5MS4dHA9Px3UdmZj1Ou3cfSboF2CEiXk6HNwOuqE94ZmZWT1nqFIYXEkJqAfDhnOIxM7MGytJ09nRJdwC/Inlo7RjgnlyjMjOzhsjSdPbpaXMXhT6ar4qI3+cblpmZNUKWKwWAvwLLI+IuSe+TNDAilucZmJmZ1V+7dQqSTgV+A1yZjtoCmJJnUGZm1hhZKpr/DdgTWAYQEc+S3KZqZmY9TJak8FZEvF0YkNSXpMLZzMx6mCxJ4V5J3wIGSNof+DVwc75hmZlZI2RJCuOBhcDjwJeA2yJiQq5RmZlZQ2S5++iMiPgR8NPCCElnpuPMzKwHyXKl8IUK407s4jjMzKwbqNZK6rHAccBWkqaWTBoILMo7MDMzq79qxUd/AV4GhgCXloxfDszOMygzM2uMaq2kPg88D3ysfuGYmVkjZXmi+QhJz0paKmmZpOWSlmVZuaQDJT0taZ6k8W3Mc5SkOZKelHRDR9+AmZl1nSx3H10CfDoi5nZkxZL6kPS7sD/QCsyQNDUi5pTMsy1wLrBnRLwmyU9Km5k1UJa7jxZ0NCGkdgPmRcRz6RPRNwKHls1zKnBFRLwGEBGvdmI7ZmbWRbJcKcyUdBNJI3hvFUZm6KN5C+DFkuFWYPeyeT4MIOl+oA8wMSL+WL4iSacBpwFsueWWGUI2M7POyJIU3g+8CRxQMi5ov49mVRhX3mZSX2BbYF9gKPB/knaKiCVrLRRxFXAVQEtLi9tdMjPLSZZOdr7YyXW3AsNKhocCL1WY58GIWAX8XdLTJEliRie3aWZmNWg3KUj6ORVaRY2Ik9pZdAawraStgPkk3XgeVzbPFOBYYJKkISTFSc9liNvMzHKQpfjolpLX/YHDee8Z/3tExGpJpwN3kNQXXBMRT0q6AJgZEVPTaQdImgO8A3wzIvy0tJlZgyiiY0X0ktYB7oqIT+YTUnUtLS0xc+bMRmzazKxpSXokIlramy/LLanltgV8C5CZWQ+UpU5hOUmdgtL/rwDn5ByXmZk1QJa7jwbWIxAzM2u8LBXNSBoH7J0OTo+IW6rNb2ZmzSlLg3gXA2cCc9K/MyV9L+/AzMys/rJcKRwM7BwRawAkXQvMImnIzszMepCsdx8NLnk9KI9AzMys8bJcKXwPmCXpHpI7kPbGVwlmZj1SlruPfiVpOjCKJCmcExGv5B2YmZnVX5aK5sOBNyNiakT8AVgp6bD8QzMzs3rLUqfw7YhYWhhIm7X+dn4hmZlZo2RJCpXmyfR8g5mZNZcsSWGmpP+WtI2krSX9EHgk78DMzKz+siSFM4C3gZuAycAK4N/yDMrMzBojy91HbwDj6xCLmZk1WJa7j+6UNLhkeENJd+QblpmZNUKW4qMh6R1HAETEa8Cm+YVkZmaNkiUprJFU7FRH0gep0GezmZk1vyy3lk4A/izp3nR4b+BL+YVkZmaNkqWi+Y+SdgX2IGnm4qyI+GfukZmZWd1laiU1Iv6ZdqwzB/iypCfyDcvMzBohy91Hm0n6uqSHgSeBPsCxuUdmZmZ112ZSkHSqpGnAvcAQ4BTg5Yj4TkQ8Xq8AzcysfqrVKVwBPAAcFxEzAST5riMzsx6sWlLYHDgS+G9JHyBp4qJfXaIyM7OGaLP4KK1c/nFE7A2MAZYCr0qaK+k/6xahmZnVTda7j1oj4r8i4l+Bw4C38g3LzMwaocP9IkTE08B3cojFzMwaLNOVgpmZ9Q5OCmZmVpTl4TVJOkHS+enwlpJ2y7JySQdKelrSPElt9skg6bOSQlJL9tDNzKyrZblS+F/gY7z7FPNykmcYqpLUJ53vIGAH4FhJO1SYbyDwNeChjDGbmVlOsiSF3SPi34CVUOxPYd0My+0GzIuI5yLibeBG4NAK810IXFJYv5mZNU6WpLAqPesPAEmbAGsyLLcF8GLJcGs6rkjSLsCwtLE9MzNrsCxJ4TLg98Cmki4C/gxkeXhNFcYVm8mQtA7wQ+Df212RdJqkmZJmLly4MMOmzcysM7L0p3C9pEdInmoWcFhEzM2w7lZgWMnwUOClkuGBwE7AdEkA/wJMlTSu0NZSSQxXAVcBtLS0uP0lM7OctJsUJO0BPBkRV6TDAyXtHhHtVQzPALaVtBUwHzgGOK4wMSKWkrS+WtjOdOAb5QnBzMzqJ0vx0Y+B10uG30jHVRURq4HTgTuAucDkiHhS0gWSxnUmWDMzy1eWZi4UEcUim4hYIylT8xgRcRtwW9m489uYd98s6zQzs/xkuVJ4TtLXJPVL/84Enss7MDMzq78sSeHLwMdJ6gVagd2B0/IMyszMGiPL3UevklQSm5lZnUyZNZ8f3PE0Ly1ZweaDB/DNsdtx2C5btL9gjbLcfbQJcCowvHT+iDgpv7DMzHqvKbPmc+7vHmfFqncAmL9kBef+7nGA3BNDlgrjPwD/B9wFvJNrNGZmxg/ueLqYEApWrHqHH9zxdLdICu+LiHNyjcLMrIeppfjnpSUrOjS+K2WpaL5F0sG5R2Jm1kMUin/mL1lB8G7xz5RZ8zMtv/ngAR0a35WyJIUzSRLDCknLJC2XtCzvwMzMGmnKrPnsefE0thp/K3tePC3zDzpUL/7J4ptjt2NAvz5rjRvQrw/fHLtd5hg6K8vdRwNzj8LMrBuptaK31uKfwja65d1HAJI2BLYF+hfGRcR9eQVlZtZItVb0bj54APMrJICOFP8ctssWdUkC5bLcknoKSRHSUOBRYA/gAeCT+YZmZlabzlb21nqm/82x2611pQH1K/6pVdY6hVHA8xExGtgFcKcGZtat1VLZW2tF72G7bMH3jvgoWwwegIAtBg/ge0d8tCFn/h2VpfhoZUSslISk9SLiKUndP92ZWa9WSxFQV5zpN6r4p1ZZkkKrpMHAFOBOSa+xdmc5Zma5aNS9/o2s6G20LHcfHZ6+nCjpHmAQ8MdcozKzXq/WO4Bqrext1jP9WrVZpyDp/en/jQp/wOMkfTRvUKf4zKyXauZ7/ZtZtSuFG4BDgEeAIOmfufT/1rlHZ2ZNrZFNPfTmIqBatJkUIuIQSQL2iYgX6hiTmfUAjS7+KWzHSaBjqt6SmnbD+fs6xWJmPYiLf5pTlruPHpQ0KiJm5B6NmXUrLv7pfbIkhdHAlyQ9D7xBWqcQESNyjczMGsrFP71TlieaDwK2IWnW4tMklc+fzjMoM2s8F//0TlmeU3geQNKmlDSIZ2bNoVHt/7j4pzllaRBvHHApsDnwKvBBYC6wY76hmVmtaikCcvFP75Sl+OhCkpZRn4mIrYAxwP25RmVmXaKWIiAX//ROWSqaV0XEIknrSFonIu6R9P3cIzMzwO3/WH1lSQpLJG0A3AdcL+lVYHW+YZkZNP4OIBf/9D5Zio8OBVYAZ5E0hPc3fPeRWV34DiCrtzavFCRdDtwQEX8pGX1t/iGZ9Sx+AMyaSbXio2eBSyVtBtwE/CoiHq1PWGY9Q6OLfwrbcRKwrNosPoqIH0XEx4B9gMXAzyXNlXS+pA/XLUKzJubiH2s27dYpRMTzEfH9iNgFOA44nOQ5hXZJOlDS05LmSRpfYfrZkuZImi3pbkkf7PA7MKuDKbPms+fF09hq/K3sefG0TP38QtcU/zRrX7/WnLI8vNYPOBA4huQZhXuB72RYrg9wBbA/0ArMkDQ1IuaUzDYLaImINyV9BbgEOLrD78IsR34AzHqTaj2v7S/pGpIf9NOA24BtIuLoiJiSYd27AfMi4rmIeBu4keROpqKIuCci3kwHHwSGduZNmOXJD4BZb1LtSuFbJL2vfSMiFndi3VsAL5YMtwK7V5n/ZOD2ShMknUaSmNhyyy07EYpZ5/kBMOtNqvW8NrrGdavSaivOKJ0AtJBUaleK5SrgKoCWlpaK6zCrppbbQv0AmPUmWR5e66xWYFjJ8FDgpfKZJO0HTADGRcRbOcZjvVShTmD+khUE79YJZK0sdhGQ9SZ5JoUZwLaStpK0LklF9dTSGSTtAlxJkhBezTEW68VqvS3UdwBZb5Kl7aNOiYjVkk4H7gD6ANdExJOSLgBmRsRU4AfABsCvJQG8EBHj8orJmlcjnwoGFwFZ75FbUgCIiNtI7loqHXd+yev98ty+9Qzd4algs94iz+Ijsy7hp4LN6ifXKwWzAjcKZ9YcnBQsd92h+Md1AmbZuPjIcufiH7Pm4SsFy52Lf6zUqlWraG1tZeXKlY0OpUfq378/Q4cOpV+/fp1a3knBMutsvYCLf6xUa2srAwcOZPjw4aS3olsXiQgWLVpEa2srW221VafW4eIjy6SWp4Jd/GOlVq5cycYbb+yEkANJbLzxxjVdhTkpWCa11Av4iWAr54SQn1r3rYuPLJOuqBdwEjDr/nylYJm0Vf7vp4LNKps6dSoXX3xxo8PoMF8p9CK1PED2zbHbrfWsAbhewOqjluO2UVavXs24ceMYN675mnLzlUIvUWvz0a4XsEao9bhty4wZMxgxYgQrV67kjTfeYMcdd+SJJ54oTl+6dCnDhw9nzZo1ALz55psMGzaMVatW8dOf/pRRo0YxcuRIPvOZz/Dmm0nnkSeeeCJnn302o0eP5pxzzmHSpEmcfvrpANx8883svvvu7LLLLuy3334sWLAAgIkTJ3LSSSex7777svXWW3PZZZcVY7juuusYMWIEI0eO5HOf+xwACxcu5DOf+QyjRo1i1KhR3H///TXth0p8pdBLVKsozvrD7noBq7euOG4rGTVqFOPGjeO8885jxYoVnHDCCey0007F6YMGDWLkyJHce++9jB49mptvvpmxY8fSr18/jjjiCE499VQAzjvvPK6++mrOOOMMAJ555hnuuusu+vTpw6RJk4rr22uvvXjwwQeRxM9+9jMuueQSLr30UgCeeuop7rnnHpYvX852223HV77yFZ555hkuuugi7r//foYMGcLixUnnl2eeeSZnnXUWe+21Fy+88AJjx45l7ty5nd4PlTgp9BJd0Xy0Wb3ledyef/75jBo1iv79+691hl5w9NFHc9NNNzF69GhuvPFGvvrVrwLwxBNPcN5557FkyRJef/11xo4dW1zmyCOPpE+fPu9ZV2trK0cffTQvv/wyb7/99lrPEHzqU59ivfXWY7311mPTTTdlwYIFTJs2jc9+9rMMGTIEgI022giAu+66izlz5hSXXbZsGcuXL2fgwIE1748CFx81kSmz5rPnxdPYavyt7HnxtA5dQrui2JpRnsft4sWLef3111m+fDkrV65kwoQJ7Lzzzuy8884AjBs3jttvv53FixfzyCOP8MlPfhJIiokuv/xyHn/8cb797W+v9UzA+uuvX3FbZ5xxBqeffjqPP/44V1555VrLrLfeesXXffr0YfXq1URExVtL16xZwwMPPMCjjz7Ko48+yvz587s0IYCTQtNwl5LWG+V53J522mlceOGFHH/88ZxzzjlcdNFFxR9bgA022IDddtuNM888k0MOOaR4BbB8+XI222wzVq1axfXXX59pW0uXLmWLLZLirmuvvbbd+ceMGcPkyZNZtGgRQLH46IADDuDyyy8vzleItSs5KTQJdylpvVFex+11111H3759Oe644xg/fjwzZsxg2rRp75nv6KOP5pe//CVHH310cdyFF17I7rvvzv5kBiRdAAAMVUlEQVT7789HPvKRTNubOHEiRx55JJ/4xCeKRULV7LjjjkyYMIF99tmHkSNHcvbZZwNw2WWXMXPmTEaMGMEOO+zAT37yk4zvODtFRJevNE8tLS0xc+bMRodRd1uNv5VKn5SAv1/8qXqHY9Zpc+fOZfvtt290GD1apX0s6ZGIaGlvWVc011kjG5UzM2uPi4/qyI3KmVl356RQR25Uzsy6Oxcf1ZEblTOz7s5JoYNqaYfF9QJm1t25+KgD/KyAmfV0Tgod4GcFzCwrN53dC3RFOyyuFzDroNmT4e4LYGkrDBoKY86HEUc1Oqqq3HR2E3H7QWZNZPZkuPlrsPRFIJL/N38tGV+D//iP/+BHP/pRcXjChAlrNYrXm5vO7lVJwXUCZk3m7gtgVdmV+KoVyfganHzyycU2iNasWcONN97I8ccfX5xe2nQ28J6ms2fMmMFjjz3G9ttvz9VXX11crtB0dqFZ7IJC09mzZs3imGOO4ZJLLilOe+qpp7jjjjt4+OGH+c53vsOqVat48sknueiii5g2bRqPPfZYMYEVms6eMWMGv/3tbznllFNq2g+V9Krio1rbZi/M02y9QJk1raWtHRuf0fDhw9l4442ZNWsWCxYsYJdddmHjjTdea57e2nR2r0oKrhMwazKDhqZFRxXG1+iUU05h0qRJvPLKK5x00klMmDCBW2+9FUhaHx03bhznnntuxaazp0yZwsiRI5k0aRLTp08vrrNa09lnn30248aNY/r06UycOLE4rTNNZw8YkF+Rda8qPnKdgFmTGXM+9Cv7fvYbkIyv0eGHH84f//hHZsyYwdixY910dirXpCDpQElPS5onaXyF6etJuimd/pCk4XnG4zoBsyYz4ij49GUwaBig5P+nL+uSu4/WXXddRo8ezVFHHVWxyAfcdHbXrljqAzwD7A+0AjOAYyNiTsk8XwVGRMSXJR0DHB4RR1dcYarWprNreSLZzGrXXZrOXrNmDbvuuiu//vWv2XbbbRsdTpfqrk1n7wbMi4jn0oBuBA4F5pTMcygwMX39G+BySYocO3lwnYCZzZkzh0MOOYTDDz+8xyWEWuWZFLYASmuIWoHd25onIlZLWgpsDPyzdCZJpwGnAWy55ZZ5xWtmvcQOO+zAc8891+gwuqU86xTeW3XOezoPyzIPEXFVRLRERMsmm2zSJcGZWeM0W4+PzaTWfZtnUmgFhpUMDwVeamseSX2BQcDiHGMyswbr378/ixYtcmLIQUSwaNEi+vfv3+l15Fl8NAPYVtJWwHzgGOC4snmmAl8AHgA+C0zLsz7BzBpv6NChtLa2snDhwkaH0iP179+foUM7/xxHbkkhrSM4HbgD6ANcExFPSroAmBkRU4GrgV9ImkdyhXBMXvGYWffQr1+/tZ7ote4l1yeaI+I24LayceeXvF4JHJlnDGZmll2veqLZzMyqc1IwM7Oi3J5ozoukhcDz6eAQyp5p6EYcW+d059ige8fn2DqnO8cGXRffByOi3Xv6my4plJI0M8tj243g2DqnO8cG3Ts+x9Y53Tk2qH98Lj4yM7MiJwUzMytq9qRwVaMDqMKxdU53jg26d3yOrXO6c2xQ5/iauk7BzMy6VrNfKZiZWRdyUjAzs6JunxS6W5eeZdseJukeSXMlPSnpzArz7CtpqaRH07/aO5fNHt8/JD2ebvc93dUpcVm672ZL2rVOcW1Xsj8elbRM0tfL5qnrfpN0jaRXJT1RMm4jSXdKejb9v2Eby34hnedZSV+oU2w/kPRU+rn9XtLgNpategzkFNtESfNLPruD21i26nc7p9huKonrH5IqdnJch/1W8bejWxxzEdFt/0ga0vsbsDWwLvAYsEPZPF8FfpK+Pga4qY7xbQbsmr4eSNL9aHl8+wK3NGj//QMYUmX6wcDtJP1a7AE81KDP+BWSB2satt+AvYFdgSdKxl0CjE9fjwe+X2G5jYDn0v8bpq83rENsBwB909ffrxRblmMgp9gmAt/I8LlX/W7nEVvZ9EuB8xu03yr+dnSHY667XykUu/SMiLeBQpeepQ4Frk1f/wYYI6lS5z1dLiJejoi/pq+XA3NJepNrFocC10XiQWCwpM3qHMMY4G8R8Xy7c+YoIu7jvX15lB5b1wKHVVh0LHBnRCyOiNeAO4ED844tIv4UEavTwQdJ+iupuzb2WxZZvtu5xZb+RhwF/Kort5lVld+Ohh9z3T0pVOrSs/xHd60uPYFCl551lRZb7QI8VGHyxyQ9Jul2STvWMawA/iTpESVdmpbLsn/zdgxtfzEbtd8KPhARL0PyJQY2rTBPd9iHJ5Fc8VXS3jGQl9PToq1r2igCafR++wSwICKebWN63fZb2W9Hw4+57p4UuqxLzzxJ2gD4LfD1iFhWNvmvJEUjI4H/B0ypY2h7RsSuwEHAv0nau2x6Q/edpHWBccCvK0xu5H7riEbvwwnAauD6NmZp7xjIw4+BbYCdgZdJimnKNfp7eyzVrxLqst/a+e1oc7EK47ps33X3pNDtu/SU1I/kQ70+In5XPj0ilkXE6+nr24B+kobUI7aIeCn9/yrwe5JL9lJZ9m+eDgL+GhELyic0cr+VWFAoTkv/v1phnobtw7SC8RDg+EgLm8tlOAa6XEQsiIh3ImIN8NM2ttnI/dYXOAK4qa156rHf2vjtaPgx192TQrFLz/Ss8hiSLjxLFbr0hDp36ZmWS14NzI2I/25jnn8p1HFI2o1kny+qQ2zrSxpYeE1SMflE2WxTgc8rsQewtHDpWidtnq01ar+VKT22vgD8ocI8dwAHSNowLSY5IB2XK0kHAucA4yLizTbmyXIM5BFbab3U4W1sM8t3Oy/7AU9FRGulifXYb1V+Oxp/zOVVu95VfyR3yDxDcqfChHTcBSRfBoD+JMUP84CHga3rGNteJJdts4FH07+DgS8DX07nOR14kuTuigeBj9cptq3TbT6Wbr+w70pjE3BFum8fB1rquO/eR/IjP6hkXMP2G0lyehlYRXImdjJJ3dTdwLPp/43SeVuAn5Use1J6/M0Dvlin2OaRlCsXjrvCHXibA7dVOwbqENsv0uNpNsmP3GblsaXD7/lu5x1bOn5S4Tgrmbfe+62t346GH3Nu5sLMzIq6e/GRmZnVkZOCmZkVOSmYmVmRk4KZmRU5KZiZWZGTgjUtSSHp0pLhb0ia2EXrniTps12xrna2c2TaUuY9ZeOHS1qRttI5R9JPJK2TTvuwpNvS1kXnSpos6QN5x2q9g5OCNbO3gCMa8KRzVZL6dGD2k4GvRsToCtP+FhE7AyNIWtA8TFJ/4FbgxxHxoYjYnqRZiU1qjdsMnBSsua0m6b/2rPIJ5Wf6kl5P/+8r6d707PoZSRdLOl7Sw0raz9+mZDX7Sfq/dL5D0uX7KOnLYEba4NuXStZ7j6QbSB7cKo/n2HT9T0j6fjrufJKHmH4i6QdtvclIGnr8C/Ah4DjggYi4uWT6PRHxhKQd0/fxaBrbttl3pVmib6MDMKvRFcBsSZd0YJmRwPYkbWQ9R/Kk6G5KOjo5Ayh0+DMc2Iekcbd7JH0I+DxJcyCjJK0H3C/pT+n8uwE7RcTfSzcmaXOSPg/+FXiNpPXNwyLiAkmfJOl7oM2OXCS9j6SJ8fOB/YFH2pj1y8CPIuL6tOmIjlyxmAG+UrAmF0nLktcBX+vAYjMiac/+LZImFgo/6o+TJIKCyRGxJpLmlZ8DPkLSzsznlfTY9RBJswSFM/KHyxNCahQwPSIWpmf915N0ANOebdLt3A/cGhFtNY9d8ADwLUnnkLQwuyLDNszW4isF6wn+h6Sp7Z+XjFtNetKTNj62bsm0t0perykZXsPa34nyNmCCpL2oMyJirQbIJO0LvNFGfJ3t9KlQp1DqSZKrl/eIiBskPQR8CrhD0ikRMa2T27ZeylcK1vQiYjEwmaTStuAfJMU1kPRm1a8Tqz5S0jppPcPWwNMkrVF+JW32uHAn0PrtrOchYB9JQ9JK6GOBezsRD8ANwMclfaowQklfxx+VtDXwXERcRtIQ3YhObsN6MScF6ykuBUrvQvopyQ/xw8DutH0WX83TJD/et5O0qrkS+BkwB/irkg7hr6SdK+5ImiM/F7iHpOXNv0ZEpSaR25UWCR0CnKGk0/Y5wIkk7e4fDTyRFjl9hKRYzaxD3EqqmZkV+UrBzMyKnBTMzKzIScHMzIqcFMzMrMhJwczMipwUzMysyEnBzMyK/j8dIVovvdvv7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed at: 1552161369.2739499\n"
     ]
    }
   ],
   "source": [
    "import scipy as sp, numpy as np\n",
    "import scipy.io as sio\n",
    "from sklearn.preprocessing import scale, StandardScaler\n",
    "from sklearn.cross_decomposition import PLSRegression\n",
    "from sklearn.model_selection import LeaveOneGroupOut, LeaveOneOut\n",
    "from sklearn.metrics import r2_score\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import scipy\n",
    "from scipy import stats\n",
    "\n",
    "\n",
    "\n",
    "X_R2_data = []\n",
    "Y_R2_data = []\n",
    "\n",
    "def make_params_R2 (scores, loadings, observed):\n",
    "    loadings_transpose = np.ndarray.transpose(loadings)\n",
    "    x_model = np.squeeze(np.dot(scores, loadings_transpose))\n",
    "    # when is squeeze necessary\n",
    "    return (R2(x_model, observed))\n",
    "\n",
    "def R2 (x_model, observed):\n",
    "    sum_sq_err = np.sum(np.square(np.squeeze(np.subtract(x_model, observed))))\n",
    "    sq_obs = np.sum(np.square(observed))\n",
    "    return 1 - (sum_sq_err/sq_obs)\n",
    "\n",
    "max_pcs = 20\n",
    "reach_90_x = False\n",
    "reach_90_y = False\n",
    "\n",
    "for i in range(max_pcs):\n",
    "    pls= PLSRegression(n_components=(i+1))\n",
    "    pls.fit(x_zscore,y_zscore)\n",
    "    new_R2_X = make_params_R2(pls.x_scores_, pls.x_loadings_, x_zscore)\n",
    "    X_R2_data.append(new_R2_X)\n",
    "    if (new_R2_X > 0.9) and not reach_90_x:\n",
    "        print((100*new_R2_X), \"% of variation for x accounted for using\",(i+1), \"PCs.\")\n",
    "        reach_90_x = True\n",
    "    \n",
    "    Y_model = np.squeeze(pls.predict(x_zscore))\n",
    "    new_R2_Y = R2(Y_model, y_zscore)\n",
    "    Y_R2_data.append(new_R2_Y)\n",
    "    if (new_R2_Y > 0.9) and not reach_90_y:\n",
    "        print((100*new_R2_Y), \"% of variation for y accounted for using\", (i+1), \"PCs.\")\n",
    "        reach_90_y = True\n",
    "\n",
    "#plt.figure(figsize =(1300,80), dpi=100)\n",
    "plt.title(\"Number of PCs vs. Variance Accounted For\")\n",
    "plt.scatter(np.arange(1,max_pcs+1), X_R2_data, label=\"x-variance\")\n",
    "plt.scatter(np.arange(1,max_pcs+1), Y_R2_data, label=\"y-variance\")\n",
    "plt.xlabel(\"Number of PCs\")\n",
    "plt.ylabel(\"Variance Accounted For\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "print(\"Completed at: {}\".format(time.time()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Leave One Out Cross Validation\n",
    "Perform LOO cross validation on our model to evaluate the predicted error of our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed at: 1552162220.2221\n"
     ]
    }
   ],
   "source": [
    "# build model and run LOO cross validation\n",
    "num_of_comp = 2\n",
    "loo = LeaveOneOut()\n",
    "\n",
    "Y_raw_loo = []\n",
    "Y_test_loo = []\n",
    "pls_loo = PLSRegression(n_components=(num_of_comp))\n",
    "\n",
    "for train_index, test_index in loo.split(X, Y):\n",
    "    \n",
    "    pls_loo.fit_transform(X[train_index], Y[train_index])\n",
    "    pls_loo.transform(X[test_index])\n",
    "    new_y_pred = pls_loo.predict(X[test_index])\n",
    "    Y_raw_loo.append(np.squeeze(new_y_pred))\n",
    "    Y_test_loo.append((Y[test_index]))\n",
    "    \n",
    "print(\"Completed at: {}\".format(time.time()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LOO Cross Validation Results [512 CASE]:\n",
      "mes cells correctly predicted:  44\n",
      "mes cells incorrecly predicted:  7\n",
      "pro cells correctly predicted:  43\n",
      "pro cells incorrectly predicted:  6\n"
     ]
    }
   ],
   "source": [
    "# threshold our results to parse our data into two categories: (0) mes and (1) pro\n",
    "Y_predict_loo = []\n",
    "for i in range(100):\n",
    "    if Y_raw_loo[i] >= 0.45:\n",
    "        Y_predict_loo.append(np.squeeze(1))\n",
    "    else:\n",
    "        Y_predict_loo.append(np.squeeze(0))\n",
    "    \n",
    "np.squeeze(Y_predict_loo)\n",
    "np.squeeze(Y_test_loo)\n",
    "\n",
    "# count number of correct predictions made by our model\n",
    "mes_correct = 0\n",
    "mes_false_pos = 0\n",
    "pro_false_pos = 0\n",
    "pro_correct = 0\n",
    "for i in range(100):\n",
    "    if (Y_predict_loo[i] == Y_test_loo[i]):\n",
    "        if Y_predict_loo[i] == 0:\n",
    "            mes_correct += 1\n",
    "        else:\n",
    "            pro_correct += 1\n",
    "    elif (Y_predict_loo[i] != Y_test_loo[i]):\n",
    "        if Y_predict_loo[i] == 0:\n",
    "            mes_false_pos += 1\n",
    "        else:\n",
    "            pro_false_pos += 1\n",
    "print(\"LOO Cross Validation Results [512 CASE]:\")\n",
    "print(\"mes cells correctly predicted: \", mes_correct)\n",
    "print(\"mes cells incorrecly predicted: \", mes_false_pos)\n",
    "print(\"pro cells correctly predicted: \", pro_correct)\n",
    "print(\"pro cells incorrectly predicted: \", pro_false_pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### K-Fold Crossvalidation\n",
    "Run K-fold with 5-10 observations held out to compare with performance of LOOCV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_of_comp = 2\n",
    "logo = LeaveOneGroupOut()\n",
    "\n",
    "group_size = 8\n",
    "groups = []\n",
    "for i in range(X.shape[0]):\n",
    "    groups.append(i % group_size)\n",
    "groups = np.array(groups)\n",
    "\n",
    "X = np.array(X)\n",
    "Y = np.array(Y)\n",
    "X = X.astype(float)\n",
    "\n",
    "Y_raw_preds = []\n",
    "Y_test_logo = []\n",
    "pls_logo = PLSRegression(n_components=(num_of_comp))\n",
    "\n",
    "for train_index, test_index in logo.split(X, groups=groups):\n",
    "    X_train = X[train_index]\n",
    "    Y_train = Y[train_index]\n",
    "    \n",
    "    # Scale\n",
    "    scaler = StandardScaler()\n",
    "    X_train_zs = scaler.fit_transform(X_train)\n",
    "    X_test_zs = scaler.transform(X[test_index])\n",
    "    \n",
    "    # Fit model\n",
    "    pls_logo.fit(X_train_zs, Y_train)\n",
    "    Y_test_pred = pls_logo.predict(X_test_zs)\n",
    "    Y_raw_preds.append(np.squeeze(np.array(Y_test_pred)))\n",
    "    Y_test_logo.append(Y[test_index])\n",
    "    \n",
    "    # Log\n",
    "    # print(np.squeeze(np.array(Y_test_pred)))\n",
    "    # print(Y[test_index])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100,)\n",
      "(100,)\n",
      "Completed at: 1552162239.877781\n"
     ]
    }
   ],
   "source": [
    "# Manually flatten\n",
    "Y_test_logo_flat = []\n",
    "for i in range(len(Y_test_logo)):\n",
    "    for j in range(len(Y_test_logo[i])):\n",
    "        Y_test_logo_flat.append(Y_test_logo[i][j])\n",
    "        \n",
    "Y_raw_preds_flat = []\n",
    "for i in range(len(Y_raw_preds)):\n",
    "    for j in range(len(Y_raw_preds[i])):\n",
    "        Y_raw_preds_flat.append(Y_raw_preds[i][j])\n",
    "        \n",
    "Y_test_logo_flat = np.array(Y_test_logo_flat)\n",
    "Y_raw_preds_flat = np.array(Y_raw_preds_flat)\n",
    "\n",
    "print(Y_raw_preds_flat.shape)\n",
    "print(Y_test_logo_flat.shape)\n",
    "print(\"Completed at: {}\".format(time.time()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1]\n",
      "[0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1\n",
      " 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1\n",
      " 1 1 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1]\n",
      "8-Fold Crossvalidation Results [512 CASE]:\n",
      "mes cells correctly predicted:  44\n",
      "mes cells incorrecly predicted:  7\n",
      "pro cells correctly predicted:  43\n",
      "pro cells incorrectly predicted:  6\n",
      "Completed at: 1552162239.8969312\n"
     ]
    }
   ],
   "source": [
    "Y_class_preds = []\n",
    "for i in range(100):\n",
    "    if Y_raw_preds_flat[i] >= 0.45:\n",
    "        Y_class_preds.append(1)\n",
    "    else:\n",
    "        Y_class_preds.append(0)\n",
    "\n",
    "print(Y_class_preds)   \n",
    "print(Y_test_logo_flat)   \n",
    "\n",
    "mes_correct = 0\n",
    "mes_false_pos = 0\n",
    "pro_false_pos = 0\n",
    "pro_correct = 0\n",
    "for i in range(100):\n",
    "    if (Y_predict_loo[i] == Y_test_loo[i]):\n",
    "        if Y_predict_loo[i] == 0:\n",
    "            mes_correct += 1\n",
    "        else:\n",
    "            pro_correct += 1\n",
    "    elif (Y_predict_loo[i] != Y_test_loo[i]):\n",
    "        if Y_predict_loo[i] == 0:\n",
    "            mes_false_pos += 1\n",
    "        else:\n",
    "            pro_false_pos += 1\n",
    "\n",
    "\n",
    "\n",
    "print(\"{}-Fold Crossvalidation Results [512 CASE]:\".format(group_size))\n",
    "print(\"mes cells correctly predicted: \", mes_correct)\n",
    "print(\"mes cells incorrecly predicted: \", mes_false_pos)\n",
    "print(\"pro cells correctly predicted: \", pro_correct)\n",
    "print(\"pro cells incorrectly predicted: \", pro_false_pos)\n",
    "\n",
    "print(\"Completed at: {}\".format(time.time()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Random Model Control\n",
    "Build set of predictions made by the random choice of 0 or 1. Use the accuracy of this model to serve as a baseline for our PLSDA model. If our PLSDA model does better than random, we have a success."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Random Prediction Results: \n",
      "mes cells correctly predicted:  27\n",
      "mes cells incorrecly predicted:  25\n",
      "pro cells correctly predicted:  25\n",
      "pro cells incorrectly predicted:  23\n"
     ]
    }
   ],
   "source": [
    "# build data set of random values from [0,1]\n",
    "Y_random = []\n",
    "for i in range(100):\n",
    "    Y_random.append(np.random.randint(0,2))\n",
    "Y_random = np.array(Y_random)\n",
    "\n",
    "mes_correct = 0\n",
    "mes_false_pos = 0\n",
    "pro_false_pos = 0\n",
    "pro_correct = 0\n",
    "\n",
    "# evaluate success of random prediction\n",
    "for i in range(100):\n",
    "    if (Y_random[i] == Y[i]):\n",
    "        if Y_random[i] == 0:\n",
    "            mes_correct += 1\n",
    "        else:\n",
    "            pro_correct += 1\n",
    "    elif (Y_random[i] != Y[i]):\n",
    "        if Y_random[i] == 0:\n",
    "            mes_false_pos += 1\n",
    "        else:\n",
    "            pro_false_pos += 1\n",
    "\n",
    "print(\"Random Prediction Results: \")\n",
    "print(\"mes cells correctly predicted: \", mes_correct)\n",
    "print(\"mes cells incorrecly predicted: \", mes_false_pos)\n",
    "print(\"pro cells correctly predicted: \", pro_correct)\n",
    "print(\"pro cells incorrectly predicted: \", pro_false_pos)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}