ps5/pset5_lda.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Discriminant Analysis (LDA) [50 pts]\n",
    "In this part of the exercise, you will re-visit the problem of predicting whether a student gets admitted into a university. However, in this part, you will build a linear discriminant analysis (LDA) classifier for this problem.\n",
    "\n",
    "LDA is a generative model for classification that assumes the class covariances are equal. Given a training dataset of positive and negative features (x, y) with y $\\in$ {0, 1} , LDA models the data x as generated from class-conditional Gaussians:\n",
    "\n",
    "$P(x, y) = P(x|y)P(y)$ where $P(y = 1) = \\pi$ and $P(x|y) = N(x;\\mu^y,\\Sigma)$\n",
    "\n",
    "where means $\\mu^y$ are class-dependent but the covariance matrix $\\Sigma$ is class-independent (the same for all classes).\n",
    "\n",
    "A novel feature $x$ is classified as a positive if $P(y = 1|x) > P(y = 0|x)$, whichis equivalent to $a(x)\\gt0$, where the linear classifier $a(x) = w^Tx+w_0$ has weights given by $w = \\Sigma^{-1}(\\mu^1-\\mu^0)$.\n",
    "\n",
    "In practice, and in this assignment, we use $a(x)\\gt$ some threshold, or equivalently, $w^Tx>T$ for some constant $T$.\n",
    "\n",
    "As we saw in lecture, LDA and logistic regression can be expressed in the same form\n",
    "\n",
    "$P(y=1|x) = \\frac{1}{1+e^{-\\theta^Tx}}.$\n",
    "\n",
    "However, they generally produce different solutions for the parameter theta."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "In this assignment, you can assume the prior probabilities for the two classes are the same (although the number of the positive and negative samples in the training data is not the same), and that the threshold $T$ is zero. As a bonus, you are encouraged to explore how the different prior probabilities shift the decision boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAF3CAYAAADgjOwXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X18nGWZ9//vkRZSastDQebmqbTRWhBKC80qVcCE8qS4gApa7Q8oghEXEFxcge3vZ6e+ti6u3LcLv91Fe4O061YKdkFwV5SnxIoE3VQqFEpuJLalawkPttBQGmlz3H/MNSFJk5lJMtdcD/N5v155ZebKNXMdOebpmPM8r/M0dxcAAADioSbqAAAAAPAOijMAAIAYoTgDAACIEYozAACAGKE4AwAAiBGKMwAAgBihOAMAAIgRijMAAIAYoTgDAACIEYozAACAGBkbdQCjcdBBB/mUKVNCPcabb76pd73rXaEeI+nIUWHkpzhyVBw5Ko4cFUZ+igs7R2vWrHnV3d9dbL9EF2dTpkxRW1tbqMdoaWlRQ0NDqMdIOnJUGPkpjhwVR46KI0eFkZ/iws6RmW0sZT+6NQEAAGKE4gwAACBGKM4AAABiJNFjzgAAqHZvv/22Nm/erJ07dxbcb7/99tP69esrFFUylStH48aN0+GHH6699tprRLcPrTgzs+9L+rikl9392GDbJEl3SZoiaYOkT7v7VjMzSTdL+pikHZIWuPtvw4oNAIC02Lx5syZOnKgpU6Yo93E6uO3bt2vixIkVjCx5ypEjd9drr72mzZs3a+rUqSO6jzC7NZdJOmvAtuslPeLu0yQ9ElyXpI9Kmhb8NEm6NcS4AABIjZ07d+rAAw8sWJihcsxMBx54YNGWzEJCK87cfbWkPw3YfK6k5cHl5ZLO67P9Xz3nCUn7m9khYcUGAECaUJjFy2gfj0qfEJBx9y2SFPw+ONh+mKQX++y3OdgGAAAS4N5775WZ6bnnnhv07wsWLNCqVatKvr8//vGPOv/88yVJa9eu1U9/+tPev7W0tOjxxx8fdoxTpkzRq6++OuzbVVpcTggYrMT0QXc0a1Ku61OZTEYtLS0hhiV1dXWFfoykI0eFkZ/iyFFx5Ki4as3Rfvvtp+3btxfdb/fu3f32++Y3v6m//du/LVscP/jBDzRnzhwtX7580Pt9++239dZbb5UUqyRNnDhRd9xxh7Zv364nnnhCv/3tb3XyySdLkn7+859rwoQJmjFjxrBidHd1dXWptrZ20L8PzNFo7Ny5c+TPR3cP7Ue5gf/r+lxvl3RIcPkQSe3B5e9J+uxg+xX6mT17toetubk59GO4u7/00r/5448f6c3N5o8/fqS/9NK/VeS45VCpHCUV+SmOHBVHjoqr1hw9++yzJe33xhtv9LueKwHKY/v27X7ooYd6e3u7T58+3d3de3p6/IorrvCjjz7aP/axj/lHP/pR/9GPfuTu7kceeaTfcMMNfuKJJ/rs2bN9zZo1fsYZZ3hdXZ3feuut7u7+hz/8wY855hjv7u72I444wg866CCfOXOm33jjjZ7JZPzQQw/1mTNn+urVq/3ll1/2T37yk15fX+/19fX+2GOPubv7q6++6qeffrrPmjXLm5qafPLkyf7KK6+UnKPRGOxxkdTmJdRPlW45u1/SxZJuDH7f12f7lWa2UtIHJb3uQfdnNejsXKH29ib19OyQJHV3b1R7e5MkKZOZH2VoAAAU9eMf/1hnnXWW3ve+92nSpEn67W9/qw0bNqi9vV1PP/20Ojs79f73v1+f//zne29zxBFHqLW1VV/5yle0YMEC/epXv9LOnTt1zDHH6PLLL+/db++999Y3vvENtbW16Z/+6Z8kSW+99ZYmTJigr371q5Kkz33uc/rKV76ik046SZs2bdKZZ56p9evXa/HixTrppJP09a9/Xf/5n/+ppUuXVjYxIxTmVBp3SmqQdJCZbZa0SLmi7G4zu1TSJkkXBLv/VLlpNH6v3FQal4QVVxx1dCzsLczyenp2qKNjIcUZAKBsstmsFi9e3Hs9P3B90aJFymazI77fO++8U9dcc40kad68ebrzzjv19ttv67Of/azGjBmjQw89VKeeemq/25xzzjmSpBkzZqirq0sTJ07UxIkTNW7cOG3btm1Yx3/44Yf17LPP9l5/4403tH37dq1evVr33HOPJOnss8/WAQccMOL/sZJCK87c/bND/GnuIPu6pCvCiiXuurs3DWs7AAAjkc1me4swM8sPJRqV1157TY8++qjWrVsnM9Pu3btlZvrEJz5R8KzF/LivmpqafmPAampqtGvXrmHF0NPTo9bWVu2zzz57/C2JZ7KyfFMM1NZOHtZ2AADiYtWqVbrooou0ceNGbdiwQS+++KKmTp2qSZMmaeXKldq9e7e2bNmi5ubmER9j4sSJ/QbqD7x+xhln9HZ5SrmzOyXplFNO0YoVKyRJDzzwgLZu3TriGCqJ4iwG6uqWqKZmfL9tNTXjVVe3JKKIAERpNN1LQKkWLVpUlvu588479YlPfKLftk996lN66aWXNG3aNM2YMUNf+tKX9JGPfGTEx2hsbNSzzz6rWbNm6a677tJf/uVf6t5779WsWbP0y1/+Urfccova2tp03HHH6f3vf7+++93vSsr9j6tXr9YJJ5ygBx98UJMnJ6PRw8rRpBmV+vp6b2trC/UYLS0tamhoCPUYUu6kgI6Oheru3qTa2smqq1uSmPFmfXPUt8kcOZV6DiUZOepvsO4mclRcteZo/fr1Ovroo4vux/JNxZUzR4M9Lma2xt3ri902LvOcVb1MZn5iirFCFi9eTHEGAMAo0K0JADGQzWZlZr2Dl/OX+bIDVB+KM4zasmXL+FABRimbzfadwLv3Mq8joPpQnGHUFixYwIcKIsPzLP14jFFtKM4AxF6hD+e+E2qmRbnOokuLND7GQCEUZygrPlQQhmr7cA6rpYgWKCAZKM5QVrz5oxIYPD8ySSpyeYyTxcx07bXX9l6/6aabij5WP/7xj/stuTQSU6ZM0auvvlry/vfff79uvPHGQY+/bNkybdkyvGW9N2zYoGOPPXZYtykFxRmAWCr04czg+fTjMQ5PZ+cKtbZOUUtLjVpbp6izc8Wo77O2tlb33HPPsAqlchRnw3XOOefo+uuvH/T4IynOwkJxVgFhvBCAtOPDuTxogUJfnZ0r1N7epO7ujZJc3d0b1d7eNOrPpbFjx6qpqUnf+c539vjbxo0bNXfuXB133HGaO3euNm3apMcff1z333+//uZv/kazZs3SCy+80O82P/nJT/TBD35Qxx9/vE477TR1dnZKyq3jecYZZ+j444/XF7/4xd73hw0bNuioo47SZZddpmOPPVbz58/Xww8/rA9/+MOaNm2afvOb30jKFWBXXnnlHsf/1re+pba2Nl122WWaNWuW3nrrLa1Zs0Yf+chHNHv2bJ155pm9hduaNWs0c+ZMzZkzR//8z/88qrwNheIsZGG9EADkMM6xsDQUuTzG5dPRsVA9PTv6bevp2aGOjoWjvu8rrrhCK1as0Ouvv95v+5VXXqmLLrpITz31lObPn68vf/nL+tCHPqRzzjlH3/72t7V27Vq95z3v6Xebk046SU888YSefPJJzZs3T//wD/8gKdc1f9JJJ+nJJ5/UOeeco02bNvXe5ve//72uvvpqPfXUU3ruuef0wx/+UI899phuuukmffOb3+x3/wOPf91116m+vl633Xab1q5dq7Fjx+qqq67SqlWrtGbNGn3+85/XwoW5HF1yySW65ZZb1NraOuqcDYUVAkJW6IWQhhUBgEoo9OGcpCIDI8NjXD7d3ZuGtX049t13X1100UW65ZZbtM8++/Rub21t1T333CNJuvDCC/W1r32t6H1t3rxZn/nMZ7Rlyxb9+c9/1tSpUyVJq1ev7r2vs88+WwcccEDvbaZOnaoZM2ZIko455hjNnTtXZqYZM2Zow4YNw/pf2tvbtW7dOp1++umSpN27d+uQQw7R66+/rm3btvWuE3rhhRfqgQceGNZ9l4KWs5CF+UIAqgUfzuVBCxRqawdf+Huo7cN1zTXX6Pbbb9ebb7455D75LvZCrrrqKl155ZV6+umn9b3vfU87d+4sevva2treyzU1Nb3Xa2pqtGvXrlL/BUm5FuZjjjlGa9eu1dq1a/X000/rwQcflLuXFP9oUZyFLOwXAgCUiiIXdXVLVFMzvt+2mprxqqtbUpb7nzRpkj796U/r9ttv7932oQ99SCtXrpQkrVixQieddJIkaeLEidq+ffug9/P666/rsMMOkyQtX768d/spp5yiFStyw4IeeOABbd26dcSxDjz+xIkT1dXVJUmaPn26Xnnlld6uy7ffflvPPPOM9t9/f+2333567LHHev+fMFCchSzsFwIAAKXKZOZr+vSlqq09UpKptvZITZ++tKzDbK699tp+Z23ecsstuuOOO3TcccfpBz/4gW6++WZJ0rx58/Ttb39bxx9//B4nBGSzWV1wwQU6+eSTddBBB/VuX7RokVavXq0TTjhBDz74oCZPHnlDx8DjL1iwQNdcc41mzZql3bt3a9WqVbruuus0c+ZMzZo1S48//rgk6Y477tAVV1yhOXPm9Ou+LSfLDxJNovr6em9rawv1GC0tLWpoaBjVfXR2rlBHx0J1d29Sbe1k1dUtSdV4s3LkKM3IT3HkqDhyVFy15mj9+vU6+uiji+63fft2TZw4sQIRJVc5czTY42Jma9y9vthtOSGgAjKZ+akqxgAAQHjo1gQAAIgRijMAAIAYoTgDACDhkjx+PI1G+3hQnAFAQjE1BiRp3Lhxeu211yjQYsLd9dprr2ncuHEjvg9OCACAhFq8eDEFGnT44Ydr8+bNeuWVVwrut3PnzlEVDNWgXDkaN26cDj/88BHfnuIMAIAE22uvvXqXNyqkpaVFxx9/fAUiSq645IhuTQBIkGw2KzPrXUImf5kWNCA9KM4AIEGy2azcvXd8Uf4yxdk7yAWSjuIMAJAqixcvjjoEYFQozgAgoRYtWhR1CABCQHEGAAlF9907stmsGhsbGYuHVOBsTQBA4mWzWTU0NKihoUFmxpxfSDRazgAAAGKE4gwAACBGKM4AAKPG2C6gfCIpzszsajNbZ2bPmNk1wbZJZvaQmT0f/D4gitgAAMMX9fQVnBCANKl4cWZmx0r6gqQPSJop6eNmNk3S9ZIecfdpkh4JrgMAUFQ2m1VzczOT8yIVomg5O1rSE+6+w913SfqFpE9IOlfS8mCf5ZLOiyA2AECJWEoKCEcUxdk6SaeY2YFmNl7SxyQdISnj7lskKfh9cASxAQCKyBdfcV1Kisl5kXQWxVwwZnappCskdUl6VtJbki5x9/377LPV3fcYd2ZmTZKaJCmTycxeuXJlqLF2dXVpwoQJoR4j6chRYVHkZ9myZVqwYEFFjzkaPIeKi1OOGhsb1dzcXHRbpcUpR3FEfooLO0eNjY1r3L2+2H6RFGf9AjD7pqTNkq6W1ODuW8zsEEkt7j690G3r6+u9ra0t1PhaWlrU0NAQ6jGSjhwVFkV+kjYJJ8+h4uKUo8GeX9lsNvIWszjlKI7IT3Fh58jMSirOojpb8+Dg92RJn5R0p6T7JV0c7HKxpPuiiA0AsKdi48uiLsyANIlqnrN/N7NnJf1E0hXuvlXSjZJON7PnJZ0eXAdQIgZnI0xxHV+G8PDYRieS4szdT3b397v7THd/JNj2mrvPdfdpwe8/RREbICXzTYkPTwDlFPXcddWMFQKAQfCmBAyNsyGBcFGcASnEhyfCRGtsejE8Ih4ozoBAmt6UkhgzgOgxPCIeKM6AAG9KQPnwuok/HqP4ojgDAJQd4zbjr9hjxPCI6FCcIXJx/PbGmxKAahfH9+ZqQXGGyMXxGzZvSsDwpWncZlql/TFKy/9BcQYAKAvGbcZf2h+jOH7ZHwmKM0Qi7d/eMHw89gCQQ3GGSKT92xuGLy3feJHDuM34S8tjlMYv+xRnAIBRGexDMMkfjNUiLY9RGr/sU5whcmn59obhy3/jbWxslJSOb7zDkZb/k1ZPoLwozhC5tHxAYfjy33ibm5slpeMb73BQ1ADllZYv+xRnAIBhS+M4HyRfWp5/FGcAYiEt33iLSXJR0zfGNI7zAeKC4gxALFTLh3qSixq6YYHKoDgDAIxKtbR6ApVCcQYAEUlCUVNKN2wSWv2AJBkbdQAAUK2SUNRks9neOM2stzsWQHhoOQMAAIgRijMAQEmS0A0LpAHFGQCgJEnohgXSgOIMAAAgRijOAAAAYoTiDAAAIEYozgAAAGKE4gwAACBGKM4AAABihOIMAAAgRijOAAAAYoTiDAAAIEYozgBUNWa9BxA3FGcAqtrixYujDgEA+qE4AwAAiJFIijMz+4qZPWNm68zsTjMbZ2ZTzezXZva8md1lZntHERuA9MtmszIzmZkk9V6mixNAHFS8ODOzwyR9WVK9ux8raYykeZK+Jek77j5N0lZJl1Y6NgDVIZvNyt3l7pLUe5niDEAcRNWtOVbSPmY2VtJ4SVsknSppVfD35ZLOiyg2AACAyFS8OHP3/5Z0k6RNyhVlr0taI2mbu+8Kdtss6bBKxwag+ixatCjqEACgH8s361fsgGYHSPp3SZ+RtE3Sj4Lri9z9vcE+R0j6qbvPGOT2TZKaJCmTycxeuXJlqPF2dXVpwoQJoR4j6chRYeSnOHJUHDkqjhwVRn6KCztHjY2Na9y9vth+Y0OLYGinSfqDu78iSWZ2j6QPSdrfzMYGrWeHS/rjYDd296WSlkpSfX29NzQ0hBpsS0uLwj5G0pGjwshPceSoOHJUHDkqjPwUF5ccRTHmbJOkE81svOVOlZor6VlJzZLOD/a5WNJ9EcSGCunsXKHW1ilqaalRa+sUdXauiDokAABiIYoxZ79WbuD/byU9HcSwVNJ1kv7azH4v6UBJt1c6NlRGZ+cKtbc3qbt7oyRXd/dGtbc3UaABAKCIztZ090XufpS7H+vuF7p7t7t3uPsH3P297n6Bu3dHERvC19GxUD09O/pt6+nZoY6OhRFFBKQTU4MAycQKAai47u5Nw9oOYGRYmgpIJoqzmEvj2Kza2snD2g6EhZYlAHFEcRZjaR2bVVe3RDU14/ttq6kZr7q6JRFFhGqVxpYllqYCko/iLMbSOjYrk5mv6dOXqrb2SEmm2tojNX36UmUy86MODUg8lqYCko/iLMbSPDYrk5mvOXM2qKGhR3PmbKAwQ8XQsgQg7ijOYoyxWUD5VVPLEktTIQ7S+NoKG8VZjDE2C8Bo8KGIOEjj2M6wUZzFGGOzgHDRsgQgjijOYo6xWUB4aFlCXzwfyoexnaNDcQYAgOh+K6dqGtsZBoozlCyNE+KiMnhDBoDSUZyhJGmdEDcK1Vio0CKBuKL7LXyM7Rw+ijOUJK0T4kaBQgWID7rfwkcuh4/iLGbi2nWY5glxEQ5aJBAGnj+oBhRnMRLnrkMmxB2daixUaJFAGMJseab7DXFBcRYjce46ZELc0aFQAeKP1yPiguIsRuLcdciEuBgNWiQwGtXY8ozqNjbqAPCO2trJQZfmntvjIJOZTzFWBtVYqPAhitHIZrO9zyEz622BBtKKlrMYoeuwOlCoAOXH6wppQnEWI3QdAkBhQ7U8M0UN0oRuzZih6xAAhkYLGaoBLWcAgEQaeKJAY2MjJwogFSjOAMQKH6wo1cApapqbm5miBqlAcZYCcV1VABgJxg4BqHYUZwkX51UFgHKjRQRDqcYpapBeFGcJF+dVBYBSLVu2rKRJRmlVw1Ao3JEmFGcJF+dVBYBSLViwgOWtAEQqTu83FGcJx4LkSDuW7gFQCXFqmac4SzhWFUDaDBw7xKLxAKoNxVnCsaoA0oaiC0ClxHWuvKIrBJjZ+yTdKinj7sea2XGSznH3vws9OpSEVQVQLTgjD0A5ZbPZ3kLMzNTc3KyGhoZIY5JKazn735JukPS2JLn7U5LmhRkUkiM3Zcc85lhDRUT9bRYAKqGUtTXHu/tv8k1+gV0hxYMEyc+xJuWm8sjPsSaJljwAQKLEqWW+lJazV83sPZJckszsfElbQo0KicAcawCAtIhTy3wpxdkVkr4n6Sgz+29J10i6PNSokAjMsQZgpOL0QQjETcHizMxqJNW7+2mS3i3pKHc/yd03jvSAZjbdzNb2+XnDzK4xs0lm9pCZPR/8PmCkx0BlMMcagJGK05xSQNwULM7cvUfSlcHlN919+2gP6O7t7j7L3WdJmq3cgKV7JV0v6RF3nybpkeA6Yow51hAntMQASItSujUfMrOvmtkRQevWJDObVKbjz5X0QtASd66k5cH25ZLOK9MxEJL8HGtSRsyxhqjREhOechW+rPYQLvKYHqUUZ59XbtzZaklrgp+2Mh1/nqQ7g8sZd98iScHvg8t0DIQoV4itVENDj+bM2UBhVma82SIOylX4stpDuPiCkh6Wf5FU/MBme0v6o6Rj3L3TzLa5+/59/r7V3fcYd2ZmTZKaJCmTycxeuXJlqHF2dXVpwoQJoR4j6chRYaPJT2Njo5qbm8scUfyMNEfLli3T8uXL99h+8cUXa8GCBWWILD6ifJ2F8TwM4z6r/b2oWE6rPT+lCDtHjY2Na9y9vuiO+W8uQ/1I2kvSlyWtCn6ulLRXsduVcL/nSnqwz/V2SYcElw+R1F7sPmbPnu1ha25uDv0YSUeOChtNfnIv0fQrx3Mo7bmq9Ots0aJFrtwUSv1+Fi1aVLb7L2XbcFTje9FwHqdqzM9whZ0jSW1eQo1USrfmrcoN3P+X4Gd2sG20Pqt3ujQl6X5JFweXL5Z0XxmOASQO43IQB2F3QQ52P3TLDR9dxe8Y7v8c5xyVUpz9hbtf7O6PBj+XSPqL0RzUzMZLOl3SPX023yjpdDN7PvjbjaM5BpBUvNmOTJxm9wZQecMt7uP8ZaCU4mx3sEKAJMnM6iTtHs1B3X2Hux/o7q/32faau89192nB7z+N5hgAqgvFa3jCLHxpKS4fvqCkRynF2d9IajazFjP7haRHJV0bblgAJN5sEQ9hFkq0FJdPNeZsuMV9Ur4MFF343N0fMbNpkqZLMknPuXt36JEBiN0bBgDESTab7X2fNLPeIr9c+0elaMuZmV0haR93f8rdfydpvJn9VfihAQCqCS3FQE4p3ZpfcPdt+SvuvlXSF8ILCQBQjWgpxmgMt7iP85eBUoqzGst3zkoyszGS9g4vJAAAgOFJ01QaRcecSfq5pLvN7LvKTW53uaSfhRoVAABAlSqlOLtOueWSvqTcCQEPSrotzKAAAACqVdFuTXfvcffvuvv5yo01a3X3Uc1zBpRTZ+cKtbZOUUtLjVpbp6izc0XUIQGx7jIBEG+lnK3ZYmb7mtkkSWsl3WFm/yv80IDiOjtXqL29Sd3dGyW5urs3qr29iQINkYvz7OMA4q2UEwL2c/c3JH1S0h3uPlvSaeGGBZSmo2Ohenp29NvW07NDHR0LI4oIAIDRKaU4G2tmh0j6tKT/CDkeYFi6uzcNazsQpqTMPo7qwvMveUopzr6h3Bmbv3f3/wrW1nw+3LCA0tTWTh7WdiBMLEWEOIpTFzuvhdKUckLAj9z9OHf/q+B6h7t/KvzQgOLq6paopmZ8v201NeNVV7ckoogAAEOJU6EYZ6W0nAGxlcnM1/TpS1Vbe6QkU23tkZo+fakymflRh4YqF+fZx5F+g3WxNzY20nKVEBRnGFSSpqfIZOZrzpwNamjo0Zw5GyjMhoE36vCQW0RpsC725ubmSJ6XjMUcPooz7IHpKaoHXQwAwsZYzOErWJyZ2VFmNtfMJgzYfla4YSFKTE8BAOlBF3vyDFmcmdmXJd0n6SpJ68zs3D5//mbYgSE6TE+RbnQxANUlTq9tCsXSFFpb8wuSZrt7l5lNkbTKzKa4+83KrbGJlKqtnRx0ae65HcmXzWZ736zNrLerAQDCFqdCMc4KdWuOcfcuSXL3DZIaJH00WLqJ4izFmJ4CAIDoFCrOXjKzWfkrQaH2cUkHSZoRdmCIDtNTVA+6GAAgfgp1a14kaVffDe6+S9JFZva9UKNC5DKZ+RRjVYAuBgCInyFbztx9s7u/NMTffhVeSACkZM01BwAon0ItZwAikp9rLj+lSX6uOUm0aAJAyjEJLRBDzDUHANWr5JYzM9u37/7u/qdQIgLAXHMAUMWKFmdm9kVJ35D0lqT8hEguqS7EuICqxlxzAFC9SunW/KqkY9x9irtPDX4ozIAQMdccAFSvUoqzFyTtKLoXgLJhrjkAqF6ljDm7QdLjZvZrSd35je7+5dCiAsBccwBQpUopzr4n6VFJT0vqCTccAACA6lZKt+Yud/9rd7/D3Zfnf0KPDAAShNUWgGSJ82u2lOKs2cyazOwQM5uU/wk9MgCsEpAgixcvjjoEAMMQ59dsKd2anwt+39BnG1NpACFjlQAAqE5FW876TJ8xtVxTaZjZ/ma2ysyeM7P1ZjYnaJF7yMyeD34fMJpjpAGtJtWNVQLiL5vNysxkZpLUeznO3SVANUvKa7akFQLM7FhJ75c0Lr/N3f91FMe9WdLP3P18M9tb0nhJfyvpEXe/0cyul3S9pOtGcYxEo9UErBIQf9lstvdN3czk7oVvACBSSXnNFm05M7NFkv7/4KdR0j9IOmekBwyWgTpF0u2S5O5/dvdtks6VlD/RYLmk80Z6jDSg1SSZytnaOdRqAKwSAADpVsoJAedLmivpJXe/RNJMSbWjOGadpFck3WFmT5rZbWb2LkkZd98iScHvg0dxjMSj1SR58q2duWWXvE9r58Mjur84rxIQty6AOFi0aFHUIQAYhji/Zq1Yk56Z/cbdP2Bma5RrOdsuaZ27HzOiA5rVS3pC0ofd/ddmdrOkNyRd5e7799lvq7vvMe7MzJokNUlSJpOZvXLlypGEUbKuri5NmDAh1GMMbp6kzkG2ZySF+z8PV3Q5ipvBH7OennerpubuEd7nw5Juk/Syct9XLpN02kgDLJvGxkY1NzeX7f54DhVHjoojR4WRn+LCzlFjY+Mad68vtl8pxdm/KDcebJ6kayV1SVobtKINm5n9D0lPuPuU4PrJyo0ve6+kBnffYmaHSGpx9+mF7qu+vt7b2tq+R097AAAXm0lEQVRGEkbJWlpa1NDQEOoxBjNwzJmUazWJ4xI+UeUoblpaapQ7kXkgU0NDuuZvLvdYDZ5DxZGj4shRYeSnuLBzZGYlFWelnK35V+6+zd2/K+l0SRePtDAL7u8lSS+aWb7wmivpWUn3S7o42HaxpPtGeow0YG3F5Bl6LFg6euiTcpYTACRd0bM1zexSd88P3t9gZmPMbJG7j2b2tqskrQjO1OyQdIlyheLdZnappE2SLhjF/acCaysmS13dkkFbO3t6LoswqvJJyllOAJB0pUylMdfMPiXpUkkHSrpD0i9Gc1B3XytpsGa9uaO5XyBK+UK6o2Ohurs3qbZ2surqlmj9+sMijgwAkCRFizN3/5yZfUa5hc93SPqsu/8q9MiABBqstXP9+pZogglRnM9yAoCkK2Wes2mSrpb075I2SLrQzMYXvBGAVGOcGQCEp5R5zn4i6f9z9y9K+oik5yX9V6hRAQAAVKlSxpx9wN3fkCTPjQD+n2Z2f7hhAQAAVKchW87M7GuS5O5vmNnAMydHPJUGAAAAhlaoW3Nen8s3DPjbWSHEAgBAKjFOE8NRqDizIS4Pdh0AAAxh8eLRTA2KalOoOPMhLg92HQAAAGVQqDibaWZvmNl2SccFl/PXZ1QoPmBYOjtXqLV1ilpaatTaOkWdnSuiDglAShXrqmTJM4zUkMWZu49x933dfaK7jw0u56/vVckggVLkF4vv7t4oydXdvVHt7U0UaABCUayrMpvNyt17lzrLX6Y4QzGlzHMGJEJHx8J+61pKUk/PDnV0LIwoIgAAho/iDKnR3b1pWNsBYLhG2lXJkmcYDoozpEZt7eRhbQeA4RppVyVdmRgOijOkRl3dEtXU9F/2taZmvOrqlkQUEQAAw0dxhtTIZOZr+vSlqq09UpKptvZITZ++VJnM/KhDA5BCdFUiLKWsrQkkRiYzn2IMQEXQVYmw0HIGAAAQIxRnAAAAMUJxBgAAECMUZwAAADFCcQYAABAjFGcAAAAxQnEGAAAQIxRnAAAAMUJxBgAxwaSmACSKMwCIjcWLF0cdAoAYoDgDUq6zc4VaW6eopaVGra1T1Nm5IuqQAAAFUJwBKdbZuULt7U3q7t4oydXdvVHt7U2pKdDSUHhms1mZmcxMknov08UJVC+KM2CU4lwgdHQsVE/Pjn7benp2qKNjYUQRlU9aCs9sNit3l7tLUu9lijOgelGcAaMQ9wKhu3vTsLYnSZoLTwDVjeIMGIW4Fwi1tZOHtT1J0lh4Llq0KOoQAMQAxRkwCnEvEOrqlqimZny/bTU141VXtySiiMonjYUnXZkAJIozYFTiXiBkMvM1ffpS1dYeKclUW3ukpk9fqkxmftShjVqaC08A1W1s1AEASVZXt0Tt7U39ujbjViBkMvNTUYwNlP+fOjoWqrt7k2prJ6uubkkq/1cA1YXiDBgFCoRopbXwBFDdIinOzGyDpO2Sdkva5e71ZjZJ0l2SpkjaIOnT7r41iviA4aBAqKzOzhUUwwBSLcoxZ43uPsvd64Pr10t6xN2nSXokuA4AveI+dQkAlEOcTgg4V9Ly4PJySedFGEvvxKLSqbGbWBTJE+eJapMk7lOXAEA5RFWcuaQHzWyNmTUF2zLuvkWSgt8HRxQb385RZg/zfCqTuE9dAgDlYPklQyp6ULND3f2PZnawpIckXSXpfnffv88+W939gEFu2ySpSZIymczslStXhhDhPEmdg2zPSArjeMnW1dWlCRMmRB1GbPX0fFo1Na8M8heeT3mlP4eq97XJ66w4clQY+Sku7Bw1Njau6TOca0iRFGf9AjDLSuqS9AVJDe6+xcwOkdTi7tML3ba+vt7b2trKHlNLS41yjXt7RKuGhp6yHy/pWlpa1NDQEHUYscXzqbhSn0P5Vu2BU5ekZe62QnidFUeOCiM/xYWdIzMrqTireLemmb3LzCbmL0s6Q9I6SfdLujjY7WJJ91U6try4TyyKpBm8h37s2EkVjiP50jypLgDkRTGVRkbSvWaWP/4P3f1nZvZfku42s0slbZJ0QQSxSUrGxKJIkstk9m25/7nf1l273lBn5woKi2Fi6hIAaVfxljN373D3mcHPMe6+JNj+mrvPdfdpwe8/VTq2PL6do7xOU03NxEG2v81ZhkDCsP4pKiFOU2nESiYzX3PmbJD0qObM2UBhhlHZvXvw7xqcZQgky+LFi6MOAVWA4gyoAMYxAgBKRXEGVEBd3RLV1Izvt41xjKhWSesazGazMjMFY6V7Lyft/0ByUJwBFcA4RuAdSesazGazcnflp57KX6Y4Q1giWfgcqEacZQgAKAUtZwnDGo3pxWOLNEtL1+CiRYuiDgFVgJazBBk4O3p+jUZJtMgkHI9tunV2rlBHx0J1d29Sbe1k1dUtqbrHNZvN9hZiZqaoV6cZqaQVk0gmWs4SpKNjYb+JcSWpp2cHc2WlAI9teuUL77QtfE+RAoSH4ixBhpoTi7myko/HNr3SWniPZlA/XYNAYRRnCcJcWenFY5teFN57otUNKIziLEGYKytd+p4AsHt3l6S9+v2dxzYd0lR4p2VQP6LHc6YwirMEYa6s9Bg4DmnXrtdkZhoz5kDx2KZLmr5UMd8XyiVpc91VGmdrJgxzZaXDYOOQ3P+ssWMn6OSTX40oKoQh/3qt9rM1Ea1sNquGhoaow0CJaDkDIsA4pOqSyczXnDkb1NDQozlzNqSiMGNQf7LEoaWKbvHSUZwBEUjTOCRUp7h+oMY1LtAtPhwUZ0AE0jQOCYiTOLQQxcXAlqrGxkZaqhKC4gyIACd3AAjbwJaq5ubm2LRU0S1eGMVZDLHGYnVI4zgkIAr5FqLGxkZJjGVKAh6bwijOinq4ooVSWpd6AYCw5FuImpubJTGWaTC0VCULxVkBuYLopooWSmld6gUAEB0K1WShOCsgVxB199sWdqHEFAsAMHJRthAlrQBKWrzVhOKsgHIWSqWOI2OKBQAYuSgLjqSdKZq0eKsJxVkB5SqUhjOOjCkWAACobhRnBeQKotp+20ZSKA1nHBlTLADJxxnXyTKaxytps94nLd5qxdqaBWQy87V+/XrV1v7bqNbEG273KOtnohp1dq5IxfqT+Zby/BeyfEu5pET+P2k32scrm832FjZm1junWFwlLd5qRXFW1GmaM+fvRnUPtbWTgy7NPbcDkKSH1d7+nVQUNIVaypP2v1QDHi/EEd2aFcA4MqCY21IzhQxnXCdLOR+vpM0llrR4qwnFWQUwjgwo5uVBtyaxoOGM62Qp5+OVtHFbSYu3mlCcVQhL9QCFHDzo1iQWNGluKU/jiQ5pfryQXBRnAGLgstR8QKa1pTytS8ul9fFCsnFCAIAYOE3Tpx+duLM1hzrDNI1nXKd54HwaHy8kG8UZgFhI2gdktU2ZwYkOQOXQrQkAIzCcyaXTIO0nOqRxPB2Si+IMAEag2lqS0jxwPq3j6ZBckRVnZjbGzJ40s/8Irk81s1+b2fNmdpeZ7R1VbABQTNpbkgZK88D5amsFRfxF2XJ2taT1fa5/S9J33H2apK2SLo0kKgAoQZpbkoaS1imBqq0VFPEXSXFmZodLOlvSbcF1k3SqpFXBLsslnRdFbABQijS3JFWbamsFRfxFdbbmP0r6mqSJwfUDJW1z913B9c2SDosiMAAoVdLOMMXg6uqW9DvzVkp/KyjizSq9Ir2ZfVzSx9z9r8ysQdJXJV0iqdXd3xvsc4Skn7r7jEFu3ySpSZIymczslStXhhpvV1eXJkyYEOoxko4cFZa8/DysXKP2y8rN3H+ZpNNCPWLyclR55Ki40eWo8s/7SuM5VFzYOWpsbFzj7vXF9ouiOPt7SRdK2iVpnKR9Jd0r6UxJ/8Pdd5nZHElZdz+z0H3V19d7W1tbqPG2tLSooaEh1GMkHTkqLEn5GTh3l5RrQQi7uy5JOYoKOSqOHBVGfooLO0dmVlJxVvExZ+5+g7sf7u5TJM2T9Ki7z5fULOn8YLeLJd1X6diAasdZawCGg/nhwhGnec6uk/TXZvZ75cag3R5xPEDV4aw1AKVifrjwRFqcuXuLu388uNzh7h9w9/e6+wXu3h1lbEA14qw1AKWipT08cWo5AxCxapy7C8DIhNHSTjdpDsUZgF7M3QWgVOVuaaeb9B1RzXMGIKaYuwtAKco9P1yhbtJqe0+i5QwAAAxbuVvaOSHpHbScAQCAESlnS3tt7eSgS3PP7dWGljMAABA5Tkh6B8UZAACIHCckvYNuTQAAEAuckJRDyxkAAECMUJwBAADECMUZAABAjFCcAVWKZVIAIJ44IQCoQvllUvKzceeXSZHEYFwAiBgtZ0AVKrRMCgAgWhRnQBVimRQAiC+KM6AKDbUcSjUukwIAcUNxBlQhlkkBgPiiOAOqEMukAEB8cbYmUKVYJgUA4omWMwAAgBihOAMAAIgRijMAAIAYoTgDAACIEYozAACAGKE4AwAAiBGKMwAAgBihOAMAAIgRijMAAIAYoTgDAACIEYozAACAGKE4AwAAiBGKMwAAgBihOAMAAIgRijMAQMk6O1eotXWKWlpq1No6RZ2dK6IOCUidsVEHAABIhs7OFWpvb1JPzw5JUnf3RrW3N0mSMpn5UYYGpErFW87MbJyZ/cbMfmdmz5jZ4mD7VDP7tZk9b2Z3mdnelY4NADC0jo6FvYVZXk/PDnV0LIwoIiCdoujW7JZ0qrvPlDRL0llmdqKkb0n6jrtPk7RV0qURxAYAGEJ396ZhbQcwMhUvzjynK7i6V/Djkk6VtCrYvlzSeZWODQAwtNraycPaDmBkIjkhwMzGmNlaSS9LekjSC5K2ufuuYJfNkg6LIjYAwODq6paopmZ8v201NeNVV7ckooiAdDJ3j+7gZvtLulfS1yXd4e7vDbYfIemn7j5jkNs0SWqSpEwmM3vlypWhxtjV1aUJEyaEeoykI0eFkZ/iyFFx8cnRw5JuU+679cGSLpN0WqQR5cUnR/FEfooLO0eNjY1r3L2+2H6Rnq3p7tvMrEXSiZL2N7OxQevZ4ZL+OMRtlkpaKkn19fXe0NAQaowtLS0K+xhJR44KIz/FkaPi4pOjBkl/F3UQg4pPjuKJ/BQXlxxFcbbmu4MWM5nZPsp95VovqVnS+cFuF0u6r9KxAQAARC2KlrNDJC03szHKFYd3u/t/mNmzklaa2d9JelLS7RHEBgAAEKmKF2fu/pSk4wfZ3iHpA5WOBwAAIE5YvgkAACBGKM4AAABihOIMAAAgRijOAAAAYoTiDAAAIEYozgAAAGKE4gwAACBGIl1bc7TM7BVJG0M+zEGSXg35GElHjgojP8WRo+LIUXHkqDDyU1zYOTrS3d9dbKdEF2eVYGZtpSxSWs3IUWHkpzhyVBw5Ko4cFUZ+iotLjujWBAAAiBGKMwAAgBihOCtuadQBJAA5Koz8FEeOiiNHxZGjwshPcbHIEWPOAAAAYoSWMwAAgBihOAuY2Tgz+42Z/c7MnjGzxcH2qWb2azN73szuMrO9o441amY2xsyeNLP/CK6Toz7MbIOZPW1ma82sLdg2ycweCnL0kJkdEHWcUTKz/c1slZk9Z2brzWwOOcoxs+nBcyf/84aZXUN++jOzrwTv1evM7M7gPZz3oj7M7OogP8+Y2TXBtqp+HpnZ983sZTNb12fboDmxnFvM7Pdm9pSZnVCpOCnO3tEt6VR3nylplqSzzOxESd+S9B13nyZpq6RLI4wxLq6WtL7PdXK0p0Z3n9XnlOzrJT0S5OiR4Ho1u1nSz9z9KEkzlXs+kSNJ7t4ePHdmSZotaYeke0V+epnZYZK+LKne3Y+VNEbSPPFe1MvMjpX0BUkfUO419nEzmyaeR8sknTVg21A5+aikacFPk6RbKxQjxVme53QFV/cKflzSqZJWBduXSzovgvBiw8wOl3S2pNuC6yZyVIpzlcuNVOU5MrN9JZ0i6XZJcvc/u/s2kaPBzJX0grtvFPkZaKykfcxsrKTxkraI96K+jpb0hLvvcPddkn4h6ROq8ueRu6+W9KcBm4fKybmS/jWoD56QtL+ZHVKJOCnO+gi669ZKelnSQ5JekLQteGJL0mZJh0UVX0z8o6SvSeoJrh8ocjSQS3rQzNaYWVOwLePuWyQp+H1wZNFFr07SK5LuCLrHbzOzd4kcDWaepDuDy+Qn4O7/LekmSZuUK8pel7RGvBf1tU7SKWZ2oJmNl/QxSUeI59FghsrJYZJe7LNfxZ5TFGd9uPvuoCvhcOWago8ebLfKRhUfZvZxSS+7+5q+mwfZtWpzFPiwu5+gXJP4FWZ2StQBxcxYSSdIutXdj5f0pqqva6WoYLzUOZJ+FHUscROMCTpX0lRJh0p6l3Kvt4Gq9r3I3dcr1837kKSfSfqdpF0Fb4SBIvt8ozgbRNDF0iLpROWaMccGfzpc0h+jiisGPizpHDPbIGmlcl0I/yhy1I+7/zH4/bJyY4U+IKkz3xwe/H45uggjt1nSZnf/dXB9lXLFGjnq76OSfuvuncF18vOO0yT9wd1fcfe3Jd0j6UPivagfd7/d3U9w91OU68p7XjyPBjNUTjYr19qYV7HnFMVZwMzebWb7B5f3Ue7Fv15Ss6Tzg90ulnRfNBFGz91vcPfD3X2Kct0tj7r7fJGjXmb2LjObmL8s6QzluhfuVy43UpXnyN1fkvSimU0PNs2V9KzI0UCf1TtdmhL56WuTpBPNbHww7jX/HOK9qA8zOzj4PVnSJ5V7PvE82tNQOblf0kXBWZsnSno93/0ZNiahDZjZccoNBByjXNF6t7t/w8zqlGslmiTpSUn/j7t3RxdpPJhZg6SvuvvHydE7glzcG1wdK+mH7r7EzA6UdLekycp9sFzg7gMHpVYNM5ul3Ekle0vqkHSJgtedyJGCMUIvSqpz99eDbTyH+rDcdEefUa6r7klJlyk3Hoj3ooCZ/VK5ccFvS/prd3+k2p9HZnanpAZJB0nqlLRI0o81SE6Cwv+flDu7c4ekS9y9rSJxUpwBAADEB92aAAAAMUJxBgAAECMUZwAAADFCcQYAABAjFGcAAAAxQnEGIFbMbLeZre3zU7HVA8zs+2b2spmtq9QxAWAgptIAECtm1uXuEyI69imSupRb7PjYCh1zjLvvrsSxACQDLWcAYs/M9jOz9vyqAmZ2p5l9Ibh8q5m1mdkzwcSk+dtsMLNvmllr8PcTzOznZvaCmV0+2HHcfbVyy9wUiuUCM1tnZr8zs9XBtjFmdpOZPW1mT5nZVcH2ucHi7k8HrXK1fWL7upk9JukCM3uPmf3MzNaY2S/N7Khy5A1AMo0tvgsAVNQ+Zra2z/W/d/e7zOxKScvM7GZJB7j7/w7+vjCYzXuMpEfM7Dh3fyr424vuPsfMviNpmXLrw46T9Iyk744wvq9LOtPd/zu/5JukJuUW4T7e3XeZ2SQzGxccc667/x8z+1dJX1JuPVpJ2unuJ0mSmT0i6XJ3f97MPijpX5RbuxZAFaI4AxA3b7n7rIEb3f0hM7tA0j9LmtnnT582sybl3s8OkfR+Sfni7P7g99OSJrj7dknbzWynme3v7ttGEN+vlCsS71ZuwW0ptxbvd919VxDrn8xspnKLc/+fYJ/lkq7QO8XZXZJkZhOUW7T7R7nVYiRJtSOIC0BKUJwBSAQzq5F0tKS3lFs7cbOZTZX0VUl/4e5bzWyZci1jefl1FXv6XM5fH9H7n7tfHrRunS1pbbBOqEkaOIDX9rhxf28Gv2skbRusIAVQnRhzBiApviJpvaTPSvq+me0laV/lipzXzSwj6aNhB2Fm73H3X7v71yW9KukISQ9KutzMxgb7TJL0nKQpZvbe4KYXSvrFwPtz9zck/SFoFZTlzBy4H4DqQXEGIG72GTCVxo1m9j5Jl0m61t1/KWm1pP/X3X8n6UnlxpB9X7kuxxEzszsltUqabmabzezSQXb7djDAf10Qx+8k3SZpk6SnzOx3kj7n7jslXaJcd+XTyrXWDTXObb6kS4PbPiPp3NH8HwCSjak0AAAAYoSWMwAAgBihOAMAAIgRijMAAIAYoTgDAACIEYozAACAGKE4AwAAiBGKMwAAgBihOAMAAIiR/wsBMIbVPnHv9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "datafile = 'data/ex2data1.txt'\n",
    "#!head $datafile\n",
    "cols = np.loadtxt(datafile,delimiter=',',usecols=(0,1,2),unpack=True) #Read in comma separated data\n",
    "##Form the usual \"X\" matrix and \"y\" vector\n",
    "X = np.transpose(np.array(cols[:-1]))\n",
    "y = np.transpose(np.array(cols[-1:]))\n",
    "m = y.size # number of training examples\n",
    "##Insert the usual column of 1's into the \"X\" matrix\n",
    "X = np.insert(X,0,1,axis=1)\n",
    "\n",
    "#Divide the sample into two: ones with positive classification, one with null classification\n",
    "pos = np.array([X[i] for i in range(X.shape[0]) if y[i] == 1])\n",
    "neg = np.array([X[i] for i in range(X.shape[0]) if y[i] == 0])\n",
    "\n",
    "def plotData():\n",
    "    plt.figure(figsize=(10,6))\n",
    "    plt.plot(pos[:,1],pos[:,2],'k+',label='Admitted')\n",
    "    plt.plot(neg[:,1],neg[:,2],'yo',label='Not admitted')\n",
    "    plt.xlabel('Exam 1 score')\n",
    "    plt.ylabel('Exam 2 score')\n",
    "    plt.legend()\n",
    "    plt.grid(True)\n",
    "    \n",
    "plotData()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement the LDA classifier by completing the code here. As an implementation detail, you should first center the positive and negative data separately, so that each has a mean equal to 0, before computing the covariance, as this tends to give a more accurate estimate. \n",
    "\n",
    "You should center the whole training data set before applying the classifier. Namely, subtract the middle value of the two classes’ means ($\\frac{1}{2}$(pos mean+neg mean)), which is on the separating plane when their prior probabilities are the same and becomes the ‘center’ of the data. [5 pts]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.         74.7189227  73.95640208] [ 1.        52.0323011 54.6203921]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAF3CAYAAAD3rnzeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt83HWd7/H3Jy1Mqa1A0Z0FSpsGa6RAaSGo0aqpRYrLblkXcasRqOhGlCLHxRWx55Cpe+rq4lmPHF2xu0BxT5bCdr3UCyu3ZisQxFQqBUrWGtvSYw1QKbSWRtt8zh/zm5C2yVySmfld5vV8PPLIzG9un3wnM/OZz/dm7i4AAABEW13YAQAAAKAwkjYAAIAYIGkDAACIAZI2AACAGCBpAwAAiAGSNgAAgBggaQMAAIgBkjYAAIAYIGkDAACIAZI2AACAGBgfdgDl8JrXvMbr6+vDDmPMfve73+lVr3pV2GFEFu1TGG1UGG1UGG2UH+1TGG2U34YNG55399eWertEJG319fXq7u4OO4wx6+zsVEtLS9hhRBbtUxhtVBhtVBhtlB/tUxhtlJ+ZbRvN7egeBQAAiAGSNgAAgBggaQMAAIiBRIxpAwCg1v3hD3/Qjh07tH///rBD0bHHHqvNmzeHHUboJkyYoKlTp+qoo44qy/2RtAEAkAA7duzQ5MmTVV9fLzMLNZY9e/Zo8uTJocYQNnfXrl27tGPHDs2YMaMs90n3KAAACbB//36dcMIJoSdsyDIznXDCCWWtfJK0AQCQECRs0VLu54OkDQAAlNX3vvc9mZmefvrpYS9fsmSJ1qxZU/T9/frXv9Z73/teSdLGjRv1wx/+cPCyzs5OPfzwwyXHWF9fr+eff77k24WJpA0AgBqWyWTKfp9r1qzRvHnztHr16rLc30knnTSY5JUraYsjkrYI6uvrUFdXvTo769TVVa++vo6wQwIAJNTy5cvLen979+7VI488oltuuWUwaXN3LV26VLNmzdKFF16oZ599dvD69fX1+uxnP6vm5mY1NTXpZz/7mRYuXKhTTz1VN998syRp69atOuOMM/T73/9eN9xwg+68807NmTNHX/ziF3XzzTfry1/+subMmaMf//jHeu6553TxxRfr3HPP1bnnnquHHnpIkrRr1y6df/75mjt3rj760Y/K3cv6d1cDs0cjpq+vQz09bRoY2CdJ6u/fpp6eNklSOt0aZmgAABT0ne98R+edd55e//rXa8qUKfrZz36mrVu3qqenR5s2bVJfX59mzZqlK664YvA2p5xyirq6uvTJT35SS5Ys0UMPPaT9+/fr9NNP15VXXjl4vaOPPlqf+9zn1N3dra9+9auSpJdfflmTJk3Spz71KUnSBz7wAX3yk5/UvHnztH37di1cuFCbN2/W8uXLNW/ePN1www36wQ9+oJUrV1a3YcqASlvE9PYuG0zYcgYG9qm3d1lIEQEAkiaTycjMBgfK506Xo6v0jjvu0MUXXyxJWrx4se644w6tX79e73//+zVu3DiddNJJeuc733nIbRYtWiRJOvPMM/WmN71JkydP1mtf+1pNmDBBu3fvLunx77vvPi1dulRz5szRokWL9NJLL2nPnj1av369PvjBD0qSLrzwQh1//PFj/lurjUpbxPT3by/pOAAApcpkMoMJmpmVratw165deuCBB7Rp0yZdc801OnjwoMxM73nPe/LOpEylUpKkurq6wdO58wcOHCgphoGBAXV1demYY4454rK4z66l0hYxqdS0ko4DABAVa9as0WWXXaYnn3xSW7du1TPPPKMZM2ZoypQpWr16tQ4ePKidO3dq3bp1o36MyZMna8+ePSOeP//88we7TqXsxAVJevvb366OjuwY8bvvvlsvvPDCqGMIC0lbxDQ0rFBd3cRDjtXVTVRDw4qQIgIQlkrM6gMO197eXrb7uuOOO/Se97znkGMXX3yxfvOb32jmzJk688wz9bGPfUzveMc7Rv0Y8+fP11NPPaU5c+bozjvv1J/92Z/p29/+9uBEhJtuuknd3d2aPXu2Zs2aNTiZob29XevXr9fZZ5+te+65R9Omxa8YYnGcPXG4pqYm7+7uDjuMMevs7FRLS4v6+jrU27tM/f3blUpNU0PDithNQhhaei+XXPtgZLRRYXFqo3J2W5UiTm0Uhqi2z+bNm3XaaaeFHYYktrEaarjnxcw2uHtTqffFmLYISqdbY5ekHW758uVUCQAAKCO6RwEgQio5qw9AvJG0oWz4sAHGLpPJyN0Hu0Vzp3kdASBpQ9nwYYMw8X+WfDzHqHUkbQBiI9+Hdrm34omCcs7qS4IkPsdAKUjaUBF82KASau1Du1KVJSpWQDyRtKEi+FBANTCOcnTilPzyHMeLmenaa68dPP+lL32p4HP1ne98R0899dSYHre+vl7PP/980ddfu3atvvCFLwz7+KtWrdKvf/3rkh4/t6F9pZG0AYi0fB/ajKNMPp7jyunr61BXV706O+vU1VWvvr6OMd9nKpXSt771Le3atavo25QjaSvVokWL9JnPfGbYxx9N0lYtJG1VVIkXCJB0fGiXBxUrDNXX16Genjb192+T5Orv36aenrYxfy6NHz9ebW1t+trXvnbEZdu2bdOCBQs0e/ZsLViwQNu3b9fDDz+stWvX6m/+5m80Z84c/fKXvzzkNt/73vf0pje9SXPnztV5552nvr4+Sdk9Ts8//3zNnTtXH/3oRwffH7Zu3ao3vOEN+shHPqIzzjhDra2tuu+++/TWt75VM2fO1KOPPiopm5gtXbr0iMf/4he/qO7ubrW2tmrOnDl6+eWXtWHDBr3jHe/QOeeco4ULF2rnzp2SpA0bNuiss85Sc3PzsH9vJZC0VUmlXiAAshhHmV8Skl+e4/Lp7V2mgYF9hxwbGNin3t5lY77vq666SnfddZdefPHFQ44vXbpUl112mR5//HG1trbqE5/4hN7ylrdo0aJFuvHGG7Vx40adeuqph9xm3rx5euSRR/TYY49p8eLF+vu//3tJ2S7+efPm6bHHHtOiRYu0ffv2wdts2bJF11xzjR5//HE9/fTT+td//Vc9+OCD+tKXvqTPf/7zh9z/4Y9/3XXXqampSR0dHdq4caPGjx+vq6++WmvWrNGGDRt0xRVXaNmybBt96EMf0k033aSurq4xt1mx2BGhSvK9QOK++wFQLfk+tOOUfGB0eI7Lp79/e0nHS/HqV79aixcv1k033aRjjjlm8HhXV5e+9a1vSZIuvfRSffrTny54Xzt27NBf/uVfaufOnfr973+vGTNmSJLWr18/eF8XXnihjj/++MHbzJgxQ2eeeaYk6fTTT9eCBQtkZjrzzDO1devWkv6Wnp4ePfHEE3rXu94lSTp48KBOPPFEvfjii9q9e/fgHqqXXnqp7r777pLuezSotFVJJV8gQK3gQ7s8qFghlRp+s/SRjpfq4x//uG655Rb97ne/G/E6ua76fK6++motXbpUmzZt0je+8Q3t37+/4O1TqdTg6bq6usHzdXV1OnDgQLF/gqRsRfr000/Xxo0btXHjRm3atEn33HOP3L2o+MsttKTNzCaY2aNm9nMze9LMlgfHZ5jZT8zsF2Z2p5kdHVaM5VTpFwgAFIvkFw0NK1RXN/GQY3V1E9XQsKIs9z9lyhS9733v0y233DJ47C1veYtWr14tSero6NC8efMkSZMnT9aePXuGvZ8XX3xRJ598siTp9ttvHzz+9re/XR0d2eFFd999t1544YVRx3r44w8939jYqOeee26wC/QPf/iDnnzySR133HE69thj9eCDDw7+PdUQZqWtX9I73f0sSXMkXWBmb5b0RUlfdveZkl6Q9OEQYyybSr9AAAAoVjrdqsbGlUqlpksypVLT1di4sqzDda699tpDluG46aabdNttt2n27Nn6l3/5F33lK1+RJC1evFg33nij5s6de8REhEwmo0suuURve9vb9JrXvGbweHt7u9avX6+zzz5b99xzj6ZNG30B5PDHX7Jkia688krNmTNHBw8e1Jo1a3TdddfprLPO0pw5c/Twww9Lkm677TZdddVVam5uPqQbuJIsNyg1TGY2UdKDkj4m6QeS/tjdD5hZs6SMuy/Md/umpibv7u6uQqRj09fXod7eZerv365UapoaGlYc8gLp7OxUS0tLeAFGHO1TGG1UGG1UGG2UX1TbZ/PmzTrttNPCDkOStGfPHk2ePDnsMCJhuOfFzDa4e1Op9xXqRAQzGydpg6TXSfqapF9K2u3uuU7nHZJODim8skunW5l0AAAARiUqlbbjJH1b0g2SbnP31wXHT5H0Q3c/c5jbtElqk6R0On1Orp88zvbu3atJkyaFHUZk0T6F0UaF0UaF0Ub5RbV9jj32WL3uda8LOwxJ2VmW48aNCzuMSNiyZcsRy5/Mnz8/fpW2HHffbWadkt4s6TgzGx9U26ZKGnZZYndfKWmllO0ejWKpulRRLblHBe1TGG1UGG1UGG2UX1TbZ/PmzZHpkqR79BUTJkzQ3Llzy3JfYc4efW1QYZOZHSPpPEmbJa2T9N7gapdL+m44EQIAEC9R6D3DK8r9fIQ5e/RESevM7HFJP5V0r7t/X9J1kv7azLZIOkHSLXnuAwBqHkt4QMpWdHbt2kXiFhHurl27dmnChAllu8/Qukfd/XFJR9QL3b1X0hurHxEAxNPy5ctJ3KCpU6dqx44deu6558IORfv37y9rshJXEyZM0NSpU8t2f5EY0wYAAMbmqKOOGtzmKWydnZ1lG8eFV7CNFQDEUCaTkZkNbqWTO03FDUgukjYAiKFMJiN3Hxy/lDtN0vYK2gJJQ9IGAEik5cuXhx0CUFYkbQAQc+3t7WGHAKAKSNoAIOboBnwFY/2QZMweBQAkRiaTGUzQzIw1y5AoVNoAAABigKQNAAAgBkjaAABlw9gxoHJI2gAAZRP2MhtMRECSMREBAJAYTERAklFpAwCMCdUtoDpI2gAAo5JLyqK6pRaLDiNpSNqABAv7QxPJFvb4tUL4/0fSkLQBCRb1D1UkD9UtoHJI2gAARSs0fo3qFlA5JG1AwjAoHJUU1fFrqBye2+ggaQPyiOObFR+qAMqJYRbRQdIG5MGbFTAyxq8B1UXSBiQYH6qoJKq3ycUwi2giaQMOk6Q3qzjGDCB8DLOIJpI24DC8WQHlw+sm+niO4oOkDQBQMYwLjb5CzxHDLKKDpA2REcVve7xZAah1UXxvrlUkbYiMKH4j580KKF2SxoUmVdKfo6T8HYcjaQMAlBXjQqMv6c9RFIsA5UDShlAl/dseRofnHwCORNKGUCX92x5GJ6nfkmsR40KjLynPUS0UAUjaAABlMdyHY5I+MJMqKc9RLRQBSNoQGUn5tofRqYVvySNJyt9IhRSoLJI2REZSPrgwOrXwLXkkJDtAeSW1CEDSBgAYtVqukCK6kvr/R9IGIHKS+i15qDgnO0NjrOUKKVBtJG0AIqcWPvDjnOzQnQuEI7SkzcxOMbN1ZrbZzJ40s2uC41PM7F4z+0Xw+/iwYgQAFK8WKqRAmMKstB2QdK27nybpzZKuMrNZkj4j6X53nynp/uA8ACRWHJKdYrpz41AlBOIstKTN3Xe6+8+C03skbZZ0sqSLJN0eXO12SX8eToQAUB1xSHbi3J0LJEUkxrSZWb2kuZJ+Iint7julbGIn6Y/CiwwAACAaLPetKbQAzCZJ+k9JK9z9W2a2292PG3L5C+5+xLg2M2uT1CZJ6XT6nNWrV1ct5krZu3evJk2aFHYYkUX7FEYbFUYbFVaojVatWqUlS5ZUL6CI4X+oMNoov/nz529w96ZSbxdq0mZmR0n6vqQfufs/BMd6JLW4+04zO1FSp7s35rufpqYm7+7urnzAFdbZ2amWlpaww4gs2qcw2qgw2qgw2ig/2qcw2ig/MxtV0hbm7FGTdIukzbmELbBW0uXB6cslfbfasQEAAETN+BAf+62SLpW0ycw2Bsc+K+kLku4ysw9L2i7pkpDiAwAAiIzQkjZ3f1CSjXDxgmrGAgAAEHWRmD0KAACA/EjaAAAAYoCkDQAAIAZI2gAAAGKApA0AACAGSNoAAABigKQNABSPTdsB1DaSNgCQtHz58rBDAIC8SNoAAABigKQNQM3KZDIyM2W3QtbgabpKAURRmHuPAkCoMpnMYIJmZnL3cAMCgDyotAEAAMQASRsASGpvbw87BADIi6QNAMSSHwCij6QNAAAgBkjaEKq+vg51ddWrs7NOXV316uvrCDskAAAiidmjCE1fX4d6eto0MLBPktTfv009PW2SpHS6NczQAACIHCptCE1v77LBhC1nYGCfenuXhRQRkEyM1wOSgaQNoenv317ScQCjwxZdQDKQtMVEEsd+pVLTSjoOVAqVKABxQNIWA7mxX/392yT54NivuCduDQ0rVFc38ZBjdXUT1dCwIqSIUKuSWIliiy4geUjaYiCpY7/S6VY1Nq5UKjVdkimVmq7GxpVMQgDKIJPJyN0Ht+bKnSZpA+KLpC0Gkjz2K51uVXPzVrW0DKi5eSsJG6qGShSAuCFpiwHGfgHlV0uVKLboQhQk8bVVbSRtMcDYLwBjwYcloiCJY0erjaQtBhj7BVQWlSgAcUDSFhOM/QIqh0oUhuL/oXwYO1peJG0AAAxBN1751NLY0WogaUPJkrjQL6qDN2oAGD2SNpQkqQv9hqEWExgqGIgquvEqj7GjY0fShpIkdaHfMJDAANFBN17l0ZZjR9IWYVHshkzyQr+oDCoYqAT+f1CLSNoiKqrdkCz0Oza1mMBQwUAlVLJSTTceooqkLaKi2g3JQr9jQwIDRB+vR0RVqEmbmd1qZs+a2RNDjk0xs3vN7BfB7+PDjDEsUe2GZKFfjAUVDIxFLVaqgaHCrrStknTBYcc+I+l+d58p6f7gfM2JcjckC/2WRy0mMHy4YiyoVKPWhZq0uft6Sb897PBFkm4PTt8u6c+rGlRE0A2ZfHzQAOXH6wpJFnalbThpd98pScHvPwo5nlDQDQkAIxupUs1SOkgyy5WZQwvArF7S9939jOD8bnc/bsjlL7j7EePazKxNUpskpdPpc1avXl2dgCto7969mjRpUthhRBbtUxhtVBhtVFic22j+/Plat25dRR8jzu1TLbRRfvPnz9/g7k0l3zA3JiCsH0n1kp4Ycr5H0onB6RMl9RS6j3POOceTYN26dWGHEGm0T2G0UWG0UWFxa6P29naXdMRPe3t7RR4vbu0TBtooP0ndPoqcKYrdo2slXR6cvlzSd0OMBUCIGJ+EYjBBAbUi7CU/7pDUJanRzHaY2YclfUHSu8zsF5LeFZxHEaK4gwIwFoxPAoBXhD179P3ufqK7H+XuU939Fnff5e4L3H1m8Pvw2aUYRlR3UAAqgQoKRlKLS+mgdkSxexSjENUdFIBSFbOAKhU4jISEHklG0pYQUd1BASgV45MAhCnK7zUkbQkR5R0UgHJgCyMA1RDlSj5JW0KwgwKSaOj4JCpwAGodSVtCsIMCkoiEDEA1xKWSP77QFczs9ZK+ruz2UmeY2WxJi9z9f1Y8OpQknW4lSUNNYIYggHLKZDKDCZqZDVb0o6aYSts/Sbpe0h8kyd0fl7S4kkEhnlgnDtUStW+/AFANBSttkia6+6O5kmHgQIXiQUzl1onLLTuSWydOEtU/AEBsRLmSX0yl7XkzO1XZvdxkZu+VtLOiUSF2WCcOAJAEUa7kF5O0XSXpG5LeYGb/T9J/k3RlRaNC7LBOHIDRiPIHJBA1eZM2M6uT1OTu50l6raQ3uPs8d99WlegQG6wTB2A0orwmFhA1eZM2dx+QtDQ4/Tt331OVqBA7rBOHqKByAyCpiukevdfMPmVmp5jZlNxPxSNDrLBOHKKCyk3llCshjsuaWHFFOyZXMUnbFcqOa1svaUPw013JoBBP6XSrmpu3qqVlQM3NW0nYyog3YURBuRJidreoLL64JFfBpM3dZwzz01CN4ABk8SacH5UbALWgYNJmZkeZ2SfMbE3ws9TMjqpGcABQDCo3lVPphHikNbF47krDF5faUEz36NclnSPpH4Ofc4JjACqIN2FEQaUT4pHuh+pyafji8opS/+Y4tVExSdu57n65uz8Q/HxI0rmVDgyodbwJj06UVzMHUHmlJvxx+oJQTNJ2MNgRQZJkZg2SDlYuJAAYPZLayql0Qkx1uTz44pJcxSRtfyNpnZl1mtl/SnpA0rWVDQvAULwJIwoqnTxRXS6PWmyvUhP+uH5BKLhhvLvfb2YzJTVKMklPu3t/xSMDMCjqbyQAEKZMJjP4Pmlmg4l/ua4fFcXMHr1K0jHu/ri7/1zSRDP7eOVDAwDUKqrLwJGK6R79K3ffnTvj7i9I+qvKhQQAqHVUlzFapSb8cfqCUEzSVme5Tl9JZjZO0tGVCwkAAGB0krzkR8ExbZJ+JOkuM7tZkku6UtJ/VDQqAAAAHKKYpO06SW2SPqbsRIR7JP1zJYMCAADAoYrZe3TA3W929/cqO5aty91Zpw2R09fXoa6uenV21qmrq159fR1hhwTEqusFQLQVM3u008xebWZTJG2UdJuZ/UPlQwOK19fXoZ6eNvX3b5Pk6u/fpp6eNhI3hC5Oq60DiLZiJiIc6+4vSfoLSbe5+zmSzqtsWEBpenuXaWBg3yHHBgb2qbd3WUgRAQBQXsUkbePN7ERJ75P0/QrHA4xKf//2ko4DlRTX1daRbPz/xV8xSdvnlJ1BusXdfxrsPfqLyoYFlCaVmlbScaCS2I4JURSlrnpeC6NTzESEf3P32e7+8eB8r7tfXPnQgOI1NKxQXd3EQ47V1U1UQ8OKkCICAIwkSglknBRTaQMiL51uVWPjSqVS0yWZUqnpamxcqXS6NezQUOPitNo6koeu+mQhaUNecVpGI51uVXPzVrW0DKi5eSsJWwl4A68c2hZhilJXPQnk2JG0YUQso1E76KoAUGlRSiDjKm/SZmZvMLMFZjbpsOMXVDas7GOYWY+ZbTGzz1T68XAkltEAgOSgqz7+RkzazOwTkr4r6WpJT5jZRUMu/nwlgwo2pf+apHdLmiXp/WY2q5KPiSOxjEay0VUB1JYovbZJIEcn396jfyXpHHffa2b1ktaYWb27f0XZPUgr6Y3KLjHSK0lmtlrSRZKeqvDjYohUalrQNXrkccRfJpMZfBM3s8EuCwCotCglkHGSr3t0nLvvlSR33yqpRdK7gy2sKp20nSzpmSHndwTHUEUsowEAQHTYSN+uzewBSX/t7huHHBsv6VZJre4+rmJBmV0iaaG7fyQ4f6mkN7r71UOu0yapTZLS6fQ5q1evrlQ4VbN3715NmjSp8BWr6j5J/yzpWUl/JOkjCmsXs2i2T7SMto1WrVqlJUuWlD+gCOL/qDDaKD/apzDaKL/58+dvcPemUm+XL2mbKumAu/9mmMve6u4PlR5mkUGZNUvKuPvC4Pz1kuTufzfc9Zuamry7u7tS4VRNZ2enWlpawg4jsmifwmijwmijwmij/Gifwmij/MxsVEnbiN2j7r5juIQtuKxiCVvgp5JmmtkMMzta0mJJayv8mEAkxWmtPABA5eSbiBAadz9gZkuV3fN0nKRb3f3JkMMCqi63Vl5u6ZXcWnmSWDwYAGpMZBfXdfcfuvvr3f1Ud2fkO2oSa+UBAHKKrrSZ2auHXt/df1uRiAAMYq08AEBOwaTNzD4q6XOSXpaUm7XgkhoqGBcAsVYeAOAVxXSPfkrS6e5e7+4zgh8SNqAKWCsPAJBTTNL2S0n7Cl4LQNml061qbFypVGq6JFMqNV2NjSuZhAAANaiYMW3XS3rYzH4iqT930N0/UbGoAAxKp1tJ0gAARSVt35D0gKRNkgYqGw4AAACGU0z36AF3/2t3v83db8/9VDwyAIghNsIG4iVOr9likrZ1ZtZmZiea2ZTcT8UjA3AIdkaIh+XLl4cdAoASxOk1W0z36AeC39cPOcaSH0AVsTMCAKBgpW3IMh8zWPIjHFRYwM4I0ZbJZGRmMjNJGjwdp24XoJbE9TVb1I4IZnaGpFmSJuSOufs3KxUUXkGFBRI7I0RdJpMZfLM3M7l7/hsACFVcX7MFK21m1i7p/wQ/8yX9vaRFFY4LASos8VXOCulIOyCwMwIA1I5iJiK8V9ICSb9x9w9JOktSqqJRYRAVlnjKVUizW1D5YIV0tIlbVHdGiHpXQhja29vDDgFACeL0mi0maXvZ3QckHQg2jX9WTEKoGios8VTuCmlUd0aI06yraiGRBeIlTq/ZYsa0dZvZcZL+SdIGSXslPVrRqDCooWHFIWPapGhUWJBfJSqk7IwAALWtmNmjH3f33e5+s6R3Sbo86CZFFUS1woL8klwhjeusKwCIu4KVNjP7sLvfIknuvtXMxplZu7vTL1IlVFjiJ8kV0rjOugKAuCtmTNsCM/thsCPCGZIekTS5wnEBsUaFFABQbgUrbe7+ATP7S2U3jN8n6f3u/lDFIwNirhYqpHGadQUAcVfMOm0zJV0j6d8lbZV0qZlNzHsjADWBcWwAUD3FdI9+T9L/cPePSnqHpF9I+mlFowIAAMAhilny443u/pIkeXbE8f8ys7WVDQsAAABDjVhpM7NPS5K7v2Rmlxx2MUt+AAAAVFG+7tHFQ05ff9hlF1QgFgAAEo1xoBiLfEmbjXB6uPMAAKAAtn7DWORL2nyE08OdBwAAQAXlS9rOMrOXzGyPpNnB6dz5M6sUHzAqfX0d6uqqV2dnnbq66tXX1xF2SAASqlCXJ1u/oVxGTNrcfZy7v9rdJ7v7+OB07vxR1QwSKEVfX4d6etrU379Nkqu/f5t6etpI3ABURKEuz0wmI3cf3PItd5qkDaUqZp02IFZ6e5cdsuenJA0M7FNv77KQIgIAYOxI2pA4/f3bSzoOAKUabZcnW79hLEjakDip1LSSjgNAqUbb5UmXKMaCpA2J09CwQnV1h26PW1c3UQ0NK0KKCACAsSNpQ+Kk061qbFypVGq6JFMqNV2NjSuVTreGHRqABKLLE9VSzN6jQOyk060kaQCqgi5PVEsolTYzu8TMnjSzATNrOuyy681si5n1mNnCMOIDAACImrAqbU9I+gtJ3xh60MxmKbvn6emSTpJ0n5m93t0PVj9EAAAtonIQAAAUcUlEQVSA6Ail0ubum929Z5iLLpK02t373f1XkrZIemN1owMAAIieqE1EOFnSM0PO7wiOAQAA1LSKdY+a2X2S/niYi5a5+3dHutkwx4bdnN7M2iS1SVI6nVZnZ+dowoyUvXv3JuLvqBTapzDaqDDaqDDaKD/apzDaqDIqlrS5+3mjuNkOSacMOT9V0q9HuP+VklZKUlNTk7e0tIzi4aKls7NTSfg7KoX2KYw2Kow2Kow2yo/2KYw2qoyodY+ulbTYzFJmNkPSTEmPhhwTAABA6MJa8uM9ZrZDUrOkH5jZjyTJ3Z+UdJekpyT9h6SrmDkKAAAQ3uzRb7v7VHdPuXva3RcOuWyFu5/q7o3ufncY8QFAmFisFcBwotY9CgA1b/ny5WGHACCCSNqAGtHX16Gurnp1dtapq6tefX0dYYcEACgBSRtQA/r6OtTT06b+/m2SXP3929TT05aYxC0JCWkmk5GZySy78lHuNF2lAHJI2oAyiXLi0Nu7TAMD+w45NjCwT729y0KKqHySkpBmMhm5u9yzS1PmTpO0AcghaQPKIOqJQ3//9pKOx0mSE1IAGIqkDSiDqCcOqdS0ko7HSRIT0vb29rBDABBBJG1AGUQ9cWhoWKG6uomHHKurm6iGhhUhRVQ+SUxI6RIFMBySNqAMop44pNOtamxcqVRquiRTKjVdjY0rlU63hh3amCU5IQWAoSq29yhQSxoaVqinp+2QLtKoJQ7pdGsikrTD5f6m3t5l6u/frlRqmhoaViTybwVQ20jagDIgcQhXUhNSABiKpA0oExKH6uvr6yBRBlAzSNoAxFJumZVcl3RumRVJJG4AEomJCEWI8qKpiBf+l8on6susAEC5UWkrgG/zKBf+l8or6susAEC5UWkrgG/zKBf+l8or6susAEC5kbQVwLd5lAv/S+XF+mwAag1JWwF8m0e5jPQ/M378lCpHkgxJXjAYAIZD0lYA3+ZRLg0NK2R29BHHDxx4iQkJo5ROt6q5eataWgbU3LyVhA1AopG0FcC3eZRLOt2qurrJw1zyB8a1ATHD/rAIA7NHi8CiqSiXgwd/O+xxxrUB8bJ8+XISN1QdlTagihgjCQAYLZI2oIoYIwm8Im6VqkwmIzPT/PnzJUlmJjOL3d+B+CJpA6qIMZLAK5YvXx52CCXJZDJyd61bt06S5O5yd5I2VA1j2oAqY4wkAGA0qLTFGPtYJhPPK5Is18VoZpLi28XY3t4edgioQVTaYop9LJOJ5zXZ+vo61Nu7TP3925VKTVNDw4qae14zmcxggmZmcvdwAxqluCWZSAYqbTHFPpbJxPOaXLmEvL9/myQfTMjjXkkleQGqh6QtptjHMpl4XpMrqQn5WCYT0MUIlIakLaZY7yuZeF6Ti4T8SFTpgNKQtMUU630lx9CJBwcP7pV01CGX87wmQ5IS8qRMJkA08H9TPJK2mGK9r2Q4fJzTgQO7ZGYaN+4E8bwmS5K+aOXWK8tNImC9MoxF3NbrCxOzR2OM9b7ib7hxTu6/1/jxk/S2tz0fUlSohNxrtdZnjyJcQ2fvIn6otAEhYpxTbUmnW9XcvFUtLQNqbt6aiISNyQTxEpWqFl3so0PSBoQoSeOcUJui+iEb1biQRRf76JC0ASFK0jgnIEqiUlGKAqpayRFK0mZmN5rZ02b2uJl928yOG3LZ9Wa2xcx6zGxhGPEB1cKEEgCVFvWqFl3sxQur0navpDPcfbak/5J0vSSZ2SxJiyWdLukCSf9oZuNCijF07EFZG5I4zgkIAxWleOL5KV4oSZu73+PuB4Kzj0iaGpy+SNJqd+93919J2iLpjWHEeLhqJ1BJ3fIGACol6hWlKKCqFW9RGNN2haS7g9MnS3pmyGU7gmOhCiOBSuqWNwCA8JDAxpvlvpGU/Y7N7pP0x8NctMzdvxtcZ5mkJkl/4e5uZl+T1OXu/ze4/BZJP3T3fx/m/tsktUlSOp0+Z/Xq1RX5O7IWS+ob5nhaUvked+/evZo0aVJw7p2ShntuTNIDZXvMODm0fTAc2qgw2qiwJLTRqlWrtGTJkorcd6H2qeRjV0Il4k3C/1AlzZ8/f4O7N5V6u4olbQUf2OxySVdKWuDu+4Jj10uSu/9dcP5HkjLu3pXvvpqamry7u7tisXZ21mmkBKqlZaCk++rr6xhxcc3Ozk61tLRIkrq66oPK3qFSqelqbt5a2h+QEEPbB8OjjQqjjQqjjfIr1D5mprA+W0ejEvHyP5SfmY0qaQtr9ugFkq6TtCiXsAXWSlpsZikzmyFppqRHw4hxqHKtpVVKNytLQQAAgKHCGtP2VUmTJd1rZhvN7GZJcvcnJd0l6SlJ/yHpKnc/GFKMg8qVQJUyTo2lIID4YwZ4vIzl+YrbzNW4xYusUPYedffX5blshaRIlZPKtWdgqVsWsbcoalW+YQRxkaus576o5SrrkmL3t9SCsT5fQ/f0jEP3aNziRRYbxhepHAlUKjVthHFqbFkE5CQl2clXWY/T31EreL4QB1FY8qNmME4NKCwpy92UWllHuMr5fMVtLbS4xVvLSNqqiHFqQGFJSXbKNYEJ1VHO5ytu48LiFm8tI2mrMrYsAvJLSrKT5Mp6EidYJPn5QnKQtAGIlKR8eCa1sp7ULfaS+nwhWZiIACBSyjVbu5pGmu2axBngSR6wn8TnC8lC0gYgcuL04ZmU2a7FSsqYQyCO6B4FgDFIymzXYiVlzOFIkjheD8lB0gYAY1BrlaekjDkcTlLH6yE5SNoAYAySXnk6XJIH7Nda1RTxw5g2ABiDhoYVh4xpk5JTeRpJnMYclqLWqqaIHyptADAGSa481Zpaq5oifqi0AcAYJbXyVGtqsWqKeKHSBmBYzKJDraFqiqij0gbgCLW29hiQQ9UUUUalDcARmEUHoBSHV+al+8IOKZFI2gAcgVl0AIo13Pp20pcYUlEBJG0AjsAsOgDFGq4yL/VTma8AkjYAR0jyqvcAyqsSlXkmQg2PpA3AEZhFB6BY5a7Ms53YyJg9CmBYzKIDUIzh1reTUqOuzOebCFXr70kkbQAAYNRyiVRv7zL1929XKjVN/f0fHHWCxUSokZG0AQCAMTm8Mt/Z2Tnq+8omfduGPV7rGNMGAAAig4lQIyNpAwAAkcFEqJHRPQoAACKFiVDDo9IGAAAQAyRtAAAAMUDSBgAAEAMkbUCNY7sYAIgHJiIANSy3XUxu9fHcdjGSGAQMABFDpQ2oYfm2iwEARAtJG1DD2C4GAOKDpA2oYSNtC8N2MQAQPaEkbWb2t2b2uJltNLN7zOyk4LiZ2U1mtiW4/Oww4gNqBdvFAEB8hFVpu9HdZ7v7HEnfl3RDcPzdkmYGP22Svh5SfEBNYLsYAIiPUGaPuvtLQ86+SpIHpy+S9E13d0mPmNlxZnaiu++sepBAjWC7GACIh9CW/DCzFZIuk/SipPnB4ZMlPTPkajuCYyRtAACgplm2qFWBOza7T9IfD3PRMnf/7pDrXS9pgru3m9kPJP2duz8YXHa/pE+7+4Zh7r9N2S5UpdPpc1avXl2JP6Oq9u7dq0mTJoUdRmTRPoXRRoXRRoXRRvnRPoXRRvnNnz9/g7s3lXq7iiVtRQdgNl3SD9z9DDP7hqROd78juKxHUkuh7tGmpibv7u6uQrSV1dnZqZaWlrDDiCzapzDaqDDaqDDaKD/apzDaKD8zG1XSFtbs0ZlDzi6S9HRweq2ky4JZpG+W9CLj2QAAAMIb0/YFM2uUNCBpm6Qrg+M/lPQnkrZI2ifpQ+GEBwAAEC1hzR69eITjLumqKocDAAAQeeyIAAAAEAMkbQAAADFA0gYAABADJG0AAAAxQNIGAAAQAyRtAAAAMUDSBgAAEAMkbQAAADFA0gYAKFlfX4e6uurV2Vmnrq569fV1hB0SkHhhbWMFAIipvr4O9fS0aWBgnySpv3+benraJEnpdGuYoQGJRqUNAFCS3t5lgwlbzsDAPvX2LgspIqA2kLQBAErS37+9pOMAyoOkDQBQklRqWknHAZQHSRsAoCQNDStUVzfxkGN1dRPV0LAipIiA2kDSBgAoSTrdqsbGlUqlpksypVLT1di4kkkIQIUxexQAULJ0upUkDagyKm0AAAAxQNIGAAAQAyRtAAAAMUDSBgAAEAMkbQAAADFA0gYAABADJG0AAAAxQNIGAAAQAyRtAAAAMUDSBgAAEAPm7mHHMGZm9pykbWHHUQavkfR82EFEGO1TGG1UGG1UGG2UH+1TGG2U33RJy9x9ZSk3SkTSlhRm1u3uTWHHEVW0T2G0UWG0UWG0UX60T2G0UWGjaSO6RwEAAGKApA0AACAGSNqipaS+7RpE+xRGGxVGGxVGG+VH+xRGGxVWchsxpg0AACAGqLQBAADEAElbyMzsb83scTPbaGb3mNlJwXEzs5vMbEtw+dlhxxoWM7vRzJ4O2uHbZnbckMuuD9qox8wWhhlnmMzsEjN70swGzKzpsMtoI0lmdkHQBlvM7DNhxxMFZnarmT1rZk8MOTbFzO41s18Ev48PM8awmdkpZrbOzDYHr7FrguO0kyQzm2Bmj5rZz4P2WR4cn2FmPwna504zOzrsWMNmZuPM7DEz+35wvuQ2ImkL343uPtvd50j6vqQbguPvljQz+GmT9PWQ4ouCeyWd4e6zJf2XpOslycxmSVos6XRJF0j6RzMbF1qU4XpC0l9IWj/0IG2UFfzNX1P2dTVL0vuDtql1q5T9vxjqM5Lud/eZku4PzteyA5KudffTJL1Z0lXB/w7tlNUv6Z3ufpakOZIuMLM3S/qipC8H7fOCpA+HGGNUXCNp85DzJbcRSVvI3P2lIWdfJSk3yPAiSd/0rEckHWdmJ1Y9wAhw93vc/UBw9hFJU4PTF0la7e797v4rSVskvTGMGMPm7pvdvWeYi2ijrDdK2uLuve7+e0mrlW2bmubu6yX99rDDF0m6PTh9u6Q/r2pQEePuO939Z8HpPcp+6J4s2kmSFHxG7Q3OHhX8uKR3SloTHK/Z9skxs6mSLpT0z8F50yjaiKQtAsxshZk9I6lVr1TaTpb0zJCr7QiO1borJN0dnKaNCqONsmiH4qXdfaeUTVgk/VHI8USGmdVLmivpJ6KdBgXdfhslPatsz8gvJe0e8mWb15v0vyV9WtJAcP4EjaKNSNqqwMzuM7Mnhvm5SJLcfZm7nyKpQ9LS3M2GuavETvUt1EbBdZYp21XRkTs0zF3VdBsNd7NhjiW2jfKgHTAmZjZJ0r9L+m+H9ZDUPHc/GAzxmapsVfu04a5W3aiiw8z+VNKz7r5h6OFhrlqwjcaXLSqMyN3PK/Kq/yrpB5Lalc26Txly2VRJvy5zaJFRqI3M7HJJfyppgb+yTg1tVFhNtVEetEPx+szsRHffGQzJeDbsgMJmZkcpm7B1uPu3gsO002HcfbeZdSo79u84MxsfVJJq/fX2VkmLzOxPJE2Q9GplK28ltxGVtpCZ2cwhZxdJejo4vVbSZcEs0jdLejFXiq81ZnaBpOskLXL3fUMuWitpsZmlzGyGspM2Hg0jxgijjbJ+KmlmMFvraGUnZ6wNOaaoWivp8uD05ZK+G2IsoQvGHt0iabO7/8OQi2gnSWb22tyMfjM7RtJ5yo77WyfpvcHVarZ9JMndr3f3qe5er+x7zwPu3qpRtBGL64bMzP5dUqOy/dzbJF3p7v8veKP4qrIzu/ZJ+pC7d4cXaXjMbIuklKRdwaFH3P3K4LJlyo5zO6Bst8Xdw99LspnZeyT9H0mvlbRb0kZ3XxhcRhtJCr7l/m9J4yTd6u4rQg4pdGZ2h6QWSa+R1Kdslf87ku6SNE3SdkmXuPvhkxVqhpnNk/RjSZv0ynikzyo7rq3m28nMZis7iH6csoWgu9z9c2bWoOyEnymSHpP0QXfvDy/SaDCzFkmfcvc/HU0bkbQBAADEAN2jAAAAMUDSBgAAEAMkbQAAADFA0gYAABADJG0AAAAxQNIGIJLM7KCZbRzyU7UNuc3sVjN71syeqNZjAkAhLPkBIJLMbK+7Twrpsd8uaa+kb7r7GVV6zHHufrAajwUgnqi0AYgNMzvWzHrMrDE4f4eZ/VVw+utm1m1mT5rZ8iG32WpmnzezruDys83sR2b2SzO7crjHcff1kvIulGpmlwR7v/7czNYHx8aZ2ZfMbJOZPW5mVwfHF5jZY8HxW80sNSS2G8zsQUmXmNmpZvYfZrbBzH5sZm8oR7sBSAb2HgUQVceY2cYh5//O3e80s6WSVpnZVyQd7+7/FFy+zN1/a2bjJN1vZrPd/fHgsmfcvdnMvixplbJ7AU6Q9KSkm0cZ3w2SFgY7mBwXHGuTNEPSXHc/YGZTzGxC8JgL3P2/zOybkj6m7O4MkrTf3edJkpndr+yuKL8wszdJ+kdJ7xxlfAAShqQNQFS97O5zDj/o7vea2SWSvibprCEXvc/M2pR9XztR0ixJuaQtt8/oJkmT3H2PpD1mtt/MjnP33aOI7yFlk8e7JOU2ET9P0s3BBtAKksizJP3K3f8ruM7tkq7SK0nbnZJkZpMkvUXSv2V3sZOU3b4NACSRtAGIGTOrk3SapJeV3bNvh5nNkPQpSee6+wtmtkrZSlpObj+/gSGnc+dH9T7o7lcG1bALJW00szmSTNLhA4XtiBsf6nfB7zpJu4dLVAFAYkwbgPj5pKTNkt4v6VYzO0rSq5VNfl40s7Skd1c6CDM71d1/4u43SHpe0imS7pF0pZmND64zRdLTkurN7HXBTS+V9J+H35+7vyTpV0EVUZZ11uHXA1C7SNoARNUxhy358QUze72kj0i61t1/LGm9pP/u7j+X9JiyY9RuVbbrctTM7A5JXZIazWyHmX14mKvdGEwseCKI4+eS/lnSdkmPm9nPJX3A3fdL+pCy3Z6blK3ujTSOrlXSh4PbPinporH8HQCShSU/AAAAYoBKGwAAQAyQtAEAAMQASRsAAEAMkLQBAADEAEkbAABADJC0AQAAxABJGwAAQAyQtAEAAMTA/wdhjRqEFelaxgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# IMPLEMENT THIS\n",
    "pos_mean = np.mean(pos,axis=0)\n",
    "neg_mean = np.mean(neg,axis=0)\n",
    "print(pos_mean,neg_mean)\n",
    "\n",
    "pos_data = pos - 0.5*(pos_mean+neg_mean)\n",
    "neg_data = neg - 0.5*(pos_mean+neg_mean)\n",
    "\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.plot(pos_data[:,1],pos_data[:,2],'k+',label='Admitted')\n",
    "plt.plot(neg_data[:,1],neg_data[:,2],'yo',label='Not admitted')\n",
    "plt.xlabel('Exam 1 score')\n",
    "plt.ylabel('Exam 2 score')\n",
    "plt.legend()\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement the LDA algorithm here (Compute the covariance on all data): [10 pts each for getting cov_all, w and y_lda]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cov_all shape: (2, 2)\n",
      "w: [0.0785182  0.07571361]\n",
      "y_lda shape: (100,)\n"
     ]
    }
   ],
   "source": [
    "# IMPLEMENT THIS\n",
    "X_data = np.concatenate((neg_data,pos_data),axis=0)\n",
    "label = np.sort(y,axis=None)\n",
    "\n",
    "cov_all = np.cov(pos_data.T[1:]) + np.cov(neg_data.T[1:]) # SHAPE: (2,2)\n",
    "print('cov_all shape:',cov_all.shape)\n",
    "\n",
    "w = np.linalg.inv(cov_all)@((pos_mean-neg_mean)[1:]) # w=cov_all^(-1)(pos_mean-neg_mean)\n",
    "print('w:',w)\n",
    "\n",
    "y_lda =  np.dot(X_data[:,1:],w) # SHAPE: (100,)\n",
    "print('y_lda shape:',y_lda.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Completing the code to compute the training set accuracy. You should get a training accuracy around 89%. [5 pts]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 0. 1. 1. 1. 1.\n",
      " 1. 0. 1. 0. 1. 1. 1. 0. 1. 1. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n",
      " 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n",
      " 1. 1. 0. 1.]\n",
      "0.89\n"
     ]
    }
   ],
   "source": [
    "# IMPLEMENT THIS\n",
    "y_lda[y_lda<0]=0\n",
    "y_lda[y_lda>0]=1\n",
    "print(y_lda)\n",
    "count = 0\n",
    "for i in range(m):\n",
    "    if label[i] == y_lda[i]:\n",
    "        count+=1\n",
    "\n",
    "accuracy = count/m\n",
    "print(accuracy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Written Problem [10 pts]\n",
    "\n",
    "Show that the log-odds decision function a(x) for LDA\n",
    "\n",
    "$a = \\ln \\frac{p(x|C_l)p(C_l)}{p(x|C_k)p(C_k)}$\n",
    "\n",
    "is linear in x, that is, we can express $a(x)=\\theta^Tx$ for some $\\theta$. Show all your steps. Hint: This is a binary problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ANSWER:\n",
    "\n",
    "We notice that this is a binary problem. So that there are two classes: class l and k. \n",
    "\n",
    "By definition, the class-conditional density is multi-Gaussian:\n",
    "$$ p(x|C_k) = \\frac{1}{(2\\pi)^{p/2}|\\Sigma_k|^{1/2}}e^{-\\frac{1}{2}(x-\\mu_k)^{T}\\Sigma_k^{-1}(x-\\mu_k)}$$\n",
    "\n",
    "And we assume $p(C_l)=p(C_k)=\\frac{1}{2}$, $\\Sigma=\\Sigma_l=\\Sigma_k$\n",
    "\n",
    "\\begin{align*}\n",
    "a &=\\log \\frac{p(x|C_l)p(C_l)}{p(x|C_k)p(C_k)} \\\\\n",
    " &= \\log p(x|C_l)-\\log p(x|C_k)+\\log \\frac{p(C_l)}{p(C_k)} \\\\\n",
    " &= \\log p(x|C_l)-\\log p(x|C_k) \\\\\n",
    " &= -\\log (2\\pi)^{p/2}|\\Sigma_l|^{1/2}-\\frac{1}{2}(x-\\mu_l)^{T}\\Sigma_l^{-1}(x-\\mu_l)+ \\log (2\\pi)^{p/2}|\\Sigma_k|^{1/2}+\\frac{1}{2}(x-\\mu_k)^{T}\\Sigma_k^{-1}(x-\\mu_k)\\\\\n",
    " &= -\\log (2\\pi)^{p/2}|\\Sigma|^{1/2}-\\frac{1}{2}(x-\\mu_l)^{T}\\Sigma^{-1}(x-\\mu_l)+ \\log (2\\pi)^{p/2}|\\Sigma|^{1/2}+\\frac{1}{2}(x-\\mu_k)^{T}\\Sigma^{-1}(x-\\mu_k)\\\\\n",
    " &= -\\frac{1}{2}(x-\\mu_l)^{T}\\Sigma^{-1}(x-\\mu_l)+\\frac{1}{2}(x-\\mu_k)^{T}\\Sigma^{-1}(x-\\mu_k)\\\\\n",
    " &= -\\frac{1}{2}(x^T\\Sigma^{-1}x-2x^T\\Sigma^{-1}\\mu_l+\\mu_l^T\\Sigma^{-1}\\mu_l)+\\frac{1}{2}(x^T\\Sigma^{-1}x-2x^T\\Sigma^{-1}\\mu_k+\\mu_k^T\\Sigma^{-1}\\mu_k)\\\\\n",
    " &= -\\frac{1}{2}(x^T\\Sigma^{-1}x-2x^T\\Sigma^{-1}\\mu_l+\\mu_l^T\\Sigma^{-1}\\mu_l-x^T\\Sigma^{-1}x+2x^T\\Sigma^{-1}\\mu_k-\\mu_k^T\\Sigma^{-1}\\mu_k)\\\\\n",
    " &= -\\frac{1}{2}(-2x^T\\Sigma^{-1}\\mu_l+2x^T\\Sigma^{-1}\\mu_k+\\mu_l^T\\Sigma^{-1}\\mu_l-\\mu_k^T\\Sigma^{-1}\\mu_k) \\\\\n",
    " &= -\\frac{1}{2}(-2x^T\\Sigma^{-1}(\\mu_l-\\mu_k)+\\mu_l^T\\Sigma^{-1}\\mu_l-\\mu_k^T\\Sigma^{-1}\\mu_k) \\\\\n",
    " &= x^T\\Sigma^{-1}(\\mu_l-\\mu_k) + const \n",
    "\\end{align*}\n",
    "\n",
    "which is linear in x. So that we can express $a(x)=\\theta^T x + \\theta_0$, $\\theta = \\Sigma^{-1}(\\mu_l-\\mu_k)$, $\\theta_0$ is the const.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CNN on MNIST using TensorFlow&trade; [50 pts]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note 1**: The following has been verified to work with the current latest version of TensorFlow (1.11)\n",
    "\n",
    "\\* Adapted from official TensorFlow&trade; tour guide.\n",
    "\n",
    "TensorFlow is a powerful library for doing large-scale numerical computation. One of the tasks at which it excels is implementing and training deep neural networks. In this assignment you will learn the basic building blocks of a TensorFlow model while constructing a deep convolutional MNIST classifier.\n",
    "\n",
    "What you are expected to implement in this tutorial:\n",
    "\n",
    "* Create a softmax regression function that is a model for recognizing MNIST digits, based on looking at every pixel in the image\n",
    "\n",
    "* Use Tensorflow to train the model to recognize digits by having it \"look\" at thousands of examples\n",
    "\n",
    "* Check the model's accuracy with MNIST test data\n",
    "\n",
    "* Build, train, and test a multilayer convolutional neural network to improve the results\n",
    "\n",
    "Here is a diagram, created with TensorBoard, of the model we will build:\n",
    "\n",
    "![tensorflow graph](data/graph.png)\n",
    "\n",
    "## Implement Utilities\n",
    "\n",
    "### Weight Initialization\n",
    "\n",
    "To create this model, we're going to need to create a lot of weights and biases. One should generally initialize weights with a small amount of noise for symmetry breaking, and to prevent 0 gradients. Since we're using ReLU neurons, it is also good practice to initialize them with a slightly positive initial bias to avoid \"dead neurons\". Instead of doing this repeatedly while we build the model, let's create two handy functions to do it for us.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\h5py\\__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "import tempfile\n",
    "\n",
    "import tensorflow as tf\n",
    "from tensorflow.examples.tutorials.mnist import input_data\n",
    "\n",
    "tf.logging.set_verbosity(tf.logging.ERROR)\n",
    "\n",
    "def weight_variable(shape):\n",
    "    \"\"\"weight_variable generates a weight variable of a given shape.\"\"\"\n",
    "    initial = tf.truncated_normal(shape, stddev=0.1)\n",
    "    return tf.Variable(initial)\n",
    "\n",
    "\n",
    "def bias_variable(shape):\n",
    "    \"\"\"bias_variable generates a bias variable of a given shape.\"\"\"\n",
    "    initial = tf.constant(0.1, shape=shape)\n",
    "    return tf.Variable(initial)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution and Pooling [5 pts]\n",
    "\n",
    "Our convolutions uses a stride of one and are zero padded so that the output is the same size as the input. Our pooling is plain old max pooling over 2x2 blocks.\n",
    "\n",
    "NOTE: FOR ALL THE FOLLOWING CODES, DO NOT IMPLEMENT YOUR OWN VERSION. USE THE BUILT-IN METHODS FROM TENSORFLOW.\n",
    "\n",
    "Take a look at [TensorFlow API Docs](https://www.tensorflow.org/api_docs/python/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IMPLEMENT THIS\n",
    "def conv2d(x, W):\n",
    "    conv2d = tf.layers.conv2d(inputs=x,\n",
    "                             filters=W,\n",
    "                             kernel_size=[5,5],\n",
    "                             padding=\"same\",\n",
    "                             )\n",
    "    return conv2d\n",
    "\n",
    "\n",
    "def max_pool_2x2(x):\n",
    "    max_pool = tf.layers.max_pooling2d(inputs=x,\n",
    "                                      pool_size=[2,2],\n",
    "                                      strides=2)\n",
    "    return max_pool"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the CNN\n",
    "\n",
    "### First Convolutional Layer[10 pts]\n",
    "\n",
    "We can now implement our first layer. It will consist of convolution, followed by max pooling. The convolution will compute 32 features for each 5x5 patch. Its weight tensor will have a shape of [5, 5, 1, 32]. The first two dimensions are the patch size, the next is the number of input channels, and the last is the number of output channels. We will also have a bias vector with a component for each output channel.\n",
    "\n",
    "To apply the layer, we first reshape x to a 4d tensor, with the second and third dimensions corresponding to image width and height, and the final dimension corresponding to the number of color channels.\n",
    "\n",
    "We then convolve x_image with the weight tensor, add the bias, apply the ReLU function, and finally max pool. The max_pool_2x2 method will reduce the image size to 14x14.\n",
    "\n",
    "### Second Convolutional Layer[5 pts]\n",
    "\n",
    "In order to build a deep network, we stack several layers of this type. The second layer will have 64 features for each 5x5 patch.\n",
    "\n",
    "### Fully Connected Layer[10 pts]\n",
    "\n",
    "Now that the image size has been reduced to 7x7, we add a fully-connected layer with 1024 neurons to allow processing on the entire image. We reshape the tensor from the pooling layer into a batch of vectors, multiply by a weight matrix, add a bias, and apply a ReLU.\n",
    "\n",
    "### SoftmaxLayer[5 pts]\n",
    "\n",
    "Finally, we add a layer of softmax regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deepnn(x):\n",
    "    \"\"\"\n",
    "    deepnn builds the graph for a deep net for classifying digits.\n",
    "    Args:\n",
    "      x: an input tensor with the dimensions (N_examples, 784), where 784 is the\n",
    "      number of pixels in a standard MNIST image.\n",
    "    Returns:\n",
    "      A tuple (y, keep_prob). y is a tensor of shape (N_examples, 10), with values\n",
    "      equal to the logits of classifying the digit into one of 10 classes (the\n",
    "      digits 0-9). keep_prob is a scalar placeholder for the probability of\n",
    "      dropout.\n",
    "    \"\"\"\n",
    "    # Reshape to use within a convolutional neural net.\n",
    "    # Last dimension is for \"features\" - there is only one here, since images are\n",
    "    # grayscale -- it would be 3 for an RGB image, 4 for RGBA, etc.\n",
    "    with tf.name_scope('reshape'):\n",
    "        x_image = tf.reshape(x,[-1,28,28,1])\n",
    "\n",
    "    # First convolutional layer - maps one grayscale image to 32 feature maps.\n",
    "    with tf.name_scope('conv1'):\n",
    "        h_conv1 = tf.nn.relu(conv2d(x_image,32))\n",
    "\n",
    "    # Pooling layer - downsamples by 2X.\n",
    "    with tf.name_scope('pool1'):\n",
    "        h_pool1 = max_pool_2x2(h_conv1)\n",
    "\n",
    "    # Second convolutional layer -- maps 32 feature maps to 64.\n",
    "    with tf.name_scope('conv2'):\n",
    "        h_conv2 = tf.nn.relu(conv2d(h_pool1,64))\n",
    "\n",
    "    # Second pooling layer.\n",
    "    with tf.name_scope('pool2'):\n",
    "        h_pool2 = max_pool_2x2(h_conv2)\n",
    "\n",
    "    # Fully connected layer 1 -- after 2 round of downsampling, our 28x28 image\n",
    "    # is down to 7x7x64 feature maps -- maps this to 1024 features.\n",
    "    with tf.name_scope('fc1'):\n",
    "        h_pool2_flat = tf.contrib.layers.flatten(h_pool2,scope='pool2flat')\n",
    "        h_fc1 = tf.layers.dense(inputs=h_pool2_flat,units=1024,activation=tf.nn.relu)\n",
    "\n",
    "    # Map the 1024 features to 10 classes, one for each digit\n",
    "    with tf.name_scope('fc2'):\n",
    "        y_conv = tf.layers.dense(inputs=h_fc1,units=10)\n",
    "    return y_conv\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete the Graph[10 pts]\n",
    "We start building the computation graph by creating nodes for the input images and target output classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracting /tmp/tensorflow/mnist/input_data\\train-images-idx3-ubyte.gz\n",
      "Extracting /tmp/tensorflow/mnist/input_data\\train-labels-idx1-ubyte.gz\n",
      "Extracting /tmp/tensorflow/mnist/input_data\\t10k-images-idx3-ubyte.gz\n",
      "Extracting /tmp/tensorflow/mnist/input_data\\t10k-labels-idx1-ubyte.gz\n"
     ]
    }
   ],
   "source": [
    "# Import data\n",
    "mnist = input_data.read_data_sets('/tmp/tensorflow/mnist/input_data', one_hot=True)\n",
    "\n",
    "x = tf.placeholder(tf.float32, [None, 784])\n",
    "y_ = tf.placeholder(tf.int32, [None, 10])\n",
    "\n",
    "# Build the graph for the deep net\n",
    "y_conv = deepnn(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can specify a loss function just as easily. Loss indicates how bad the model's prediction was on a single example; we try to minimize that while training across all the examples. Here, our loss function is the cross-entropy between the target and the softmax activation function applied to the model's prediction. As in the beginners tutorial, we use the stable formulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tf.name_scope('loss'):\n",
    "    y_dict = dict(labels=y_,logits=y_conv)\n",
    "    losses = tf.nn.softmax_cross_entropy_with_logits(**y_dict)\n",
    "    cross_entropy = tf.reduce_mean(losses)\n",
    "\n",
    "with tf.name_scope('adam_optimizer'):\n",
    "    trainer = tf.train.AdamOptimizer(learning_rate=0.001)\n",
    "    train_step = trainer.minimize(cross_entropy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll figure out where we predicted the correct label. tf.argmax is an extremely useful function which gives you the index of the highest entry in a tensor along some axis. For example, tf.argmax(y,1) is the label our model thinks is most likely for each input, while tf.argmax(y\\_,1) is the true label. We can use tf.equal to check if our prediction matches the truth.\n",
    "\n",
    "That gives us a list of booleans. To determine what fraction are correct, we cast to floating point numbers and then take the mean. For example, [True, False, True, True] would become [1,0,1,1] which would become 0.75."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tf.name_scope('accuracy'):\n",
    "    y_pred = tf.argmax(tf.nn.softmax(y_conv),axis=1)\n",
    "    y_true = tf.argmax(y_,axis=1)\n",
    "    correct_prediction = tf.equal(tf.cast(y_pred,tf.int64),y_true)    \n",
    "    accuracy = tf.reduce_mean(tf.cast(correct_prediction,tf.float32))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving graph to: ./graph\n"
     ]
    }
   ],
   "source": [
    "# For saving the graph, DO NOT CHANGE.\n",
    "graph_location = r'./graph'\n",
    "print('Saving graph to: %s' % graph_location)\n",
    "train_writer = tf.summary.FileWriter(graph_location)\n",
    "train_writer.add_graph(tf.get_default_graph())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train and Evaluate the Model[5 pts]\n",
    "\n",
    "We will use a more sophisticated ADAM optimizer instead of a Gradient Descent Optimizer.\n",
    "\n",
    "We will add logging to every 100th iteration in the training process.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feel free to run this code. Be aware that it does 20,000 training iterations and may take a while (possibly up to half an hour), depending on your processor.\n",
    "\n",
    "The final test set accuracy after running this code should be approximately 99.2%.\n",
    "\n",
    "We have learned how to quickly and easily build, train, and evaluate a fairly sophisticated deep learning model using TensorFlow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "step 0, training accuracy 0.08\n",
      "step 100, training accuracy 0.98\n",
      "step 200, training accuracy 0.92\n",
      "step 300, training accuracy 1\n",
      "step 400, training accuracy 0.98\n",
      "step 500, training accuracy 0.98\n",
      "step 600, training accuracy 1\n",
      "step 700, training accuracy 0.98\n",
      "step 800, training accuracy 1\n",
      "step 900, training accuracy 1\n",
      "step 1000, training accuracy 1\n",
      "step 1100, training accuracy 0.98\n",
      "step 1200, training accuracy 1\n",
      "step 1300, training accuracy 0.96\n",
      "step 1400, training accuracy 0.98\n",
      "step 1500, training accuracy 1\n",
      "step 1600, training accuracy 1\n",
      "step 1700, training accuracy 1\n",
      "step 1800, training accuracy 1\n",
      "step 1900, training accuracy 0.98\n",
      "test accuracy 0.9905\n"
     ]
    }
   ],
   "source": [
    "with tf.Session() as sess:\n",
    "    sess.run(tf.global_variables_initializer())\n",
    "    for i in range(2000):\n",
    "        batch = mnist.train.next_batch(50)\n",
    "        if i % 100 == 0:\n",
    "            train_accuracy = accuracy.eval(feed_dict={x: batch[0], y_: batch[1]})\n",
    "            print('step %d, training accuracy %g' % (i, train_accuracy))\n",
    "        train_step.run(feed_dict={x: batch[0], y_: batch[1]})\n",
    "\n",
    "    print('test accuracy %g' % accuracy.eval(feed_dict={x: mnist.test.images, y_: mnist.test.labels}))\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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}