{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2021-12-07 23:22:03.818622: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:03.818648: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.\n"
     ]
    }
   ],
   "source": [
    "import demjson\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "import json\n",
    "import random\n",
    "import scipy\n",
    "import sklearn\n",
    "import string\n",
    "import tensorflow as tf\n",
    "from collections import defaultdict # Dictionaries that take a default for missing entries\n",
    "from sklearn import linear_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data is available at http://cseweb.ucsd.edu/~jmcauley/pml/data/. Download and save to your own directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dataDir = \"/home/jmcauley/pml_data/\"\n",
    "dataDir = \"/home/julian/backup/jmcauley/pml_data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sigmoid function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    return 1.0 / (1.0 + math.exp(-x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4lPW5//H3nRAW2Q0gq6IsVXCXYhUtFaQup8qpioC1iBcIcgpCcTmAPRhwg1K0UrSIgpTKYhCPAhXxp4AoogIqICCLLCbsspSwxEBynz8y8gthgADJPLN8Xtc1V57nme/M85lMMvd8v89m7o6IiEhhSUEHEBGR6KQCISIiYalAiIhIWCoQIiISlgqEiIiEpQIhIiJhqUBITDGz35nZ+9G2XjOba2Zdj3OfmdlrZrbbzL4ouZRh1z3TzO6L5DolfpiOg5BoY2bXAX8GmgK5wEqgj7svDDTYCZjZXOB1d381zH3XA5OAn7n7/hLMkAY0dPd7S2odklhKBR1ApCAzqwTMAHoA6UBp4HrgxyBznaHzgA0lWRxESoKGmCTaNAZw90nunuvuB939fXdfCmBmnc3sk58am9mvzWyVmf3bzF4ys49+GuoJtZ1vZs+b2R4zW2dm14aWZ5jZ9oLDL2ZW2czGm9kOM9toZn8ys6TjrLeNmX0bWu9IwMK9GDPrArwKXGNm+8xsUOHnCrVzM2sYmh5nZi+a2b/MLMvMPjezBgXaNjWz/2dmu8xsm5kNMLObgQFA+9B6loTaHhn6MrOk0GvaGHrt482scui++qEM95nZ92b2g5k9ftrvosQFFQiJNquBXDP7h5ndYmZVj9fQzKoBbwL9gVRgFXBtoWZXA0tD908EJgM/BxoC9wIjzaxCqO3fgMrABUBLoBNw/3HWOxX4E1AN+A5oES6ju48BHgQWuHsFd3/iZL+AkI7AIKAqsBZ4OrTuisAHwHtA7dDr+NDd3wOeAd4IreeyMM/ZOXS7IfQaKwAjC7W5DvgZ0BoYaGYXFTGvxCEVCIkq7r6X/A8pB14BdpjZNDM7J0zzW4Hl7v6Wux8GRgBbC7VZ7+6vuXsu8AZQDxjs7j+6+/tADtDQzJKB9kB/d89y9w3AcOD3x1nvCnd/090PAX8Ns94z9Za7fxF6XROAy0PLfwNsdffh7p4dyvp5EZ/zd8Bz7r7O3feRX1g7mFnBoeZBoV7bEmAJEK7QSIJQgZCo4+4r3b2zu9cFLib/m/JfwzStDWQUeJwDmYXabCswfTDUrvCyCuT3BEoDGwvctxGoU8T1ZoRpdyYKFpwDoYyQX+C+O83nrM2xr68UULD4Hm+9koBUICSqufu3wDjyC0VhW4C6P82YmRWcP0U/AIfI36D8k3OBTcdZb71C660Xpt3x7AfOKvD4mqfw2AygwXHuO9kuiZs59vUd5ugiKnKECoREFTO70MweNrO6ofl65I/Hfxam+b+AS8zsP0PDJH8ATuXD9ojQEFQ68LSZVTSz84C+wOvHWW9TM7sjtN6HTnG9S0KPv9zMygJpp/DYGUBNM+tjZmVCWa8O3bcNqP/ThvUwJgF/NLPzQ9tdftpmcfgU1i8JRAVCok0W+RuWPzez/eQXhm+Ahws3dPcfgHbkHzOxE2gCLOL0d4ntRf63+3XAJ+Rv1B57gvUOCa23ETC/qCtx99XAYPI3Nq8Jrauoj80C2gC3kT8ctIb8jc4AU0I/d5rZl2EePhb4JzAPWA9kk/+aRcLSgXISN0LfnDOB37n7nKDziMQ69SAkppnZTWZWxczKkH8cgBF+OEpETpEKhMS6a8jfq+cH8odd/tPdDwYbSSQ+aIhJRETCUg9CRETCiumT9VWrVs3r168fdAyRY2zevJnatWsHHUMkrMWLF//g7tVP1i6mC0T9+vVZtGhR0DFEjpGWlkZaWlrQMUTCMrONJ2+lISYRETkOFQgREQlLBUJERMJSgRARkbBUIEREJCwVCBERCUsFQkREwlKBEBGRsFQgREQkLBUIEREJSwVCRETCUoEQEZGwVCBERCSsiBQIMxtrZtvN7Jvj3G9mNsLM1prZUjO7MhK5RETk+CLVgxgH3HyC+28BGoVu3YC/RyCTiIicQESuB+Hu88ys/gmatAXGe/71Tz8LXYS+lrtviUQ+EZFY5u64O7m5ueTl5R35WXC64LKiipYLBtUBMgrMZ4aWHVMgzKwb+b0Mzj333IiEExH5SW5uLvv27SMrK+uo24EDB8jOziY7O5uDBw+edDonJ4dDhw5x6NAhDh8+fMLp4y3Lzc0lNzeX/O/WxS9aCoSFWRb2Fbv7aGA0QLNmzUrmtyIicc3d2bt3Lzt27GDnzp1H3X744Yej5nfv3s3evXvJyspi3759HDhw4JTXV7p0acqWLUvZsmUpV64cZcqUoXTp0qSkpJCSkkKpUqVISUmhfPnyR6YLLg83/dMtKSmJpKQkkpOTw/4Mt6x79+5Fyh0tBSITqFdgvi6wOaAsIhLj9u/fz3fffcd3331HZmYmmzZtOua2f//+sI9NSkri7LPPJjU1ldTUVOrUqcOFF15IxYoVT3grX778UUXgp+kyZcqQnJwc4d/AicVagZgG9DSzycDVwL+1/UFETuTQoUOsXr2a5cuXs3bt2qNuW7Yc/fGRkpJC7dq1qVOnDpdddhm33nortWvXpkaNGlSrVu1IMUhNTaVKlSokJekIAIhQgTCzScCvgGpmlgk8AaQAuPso4F3gVmAtcAC4PxK5RCQ27N+/n8WLF7Nw4UKWLFnC0qVLWblyJTk5OUfa1KpVi4YNG3LzzTfTsGFDGjZsSIMGDahXrx7VqlXTh/5piNReTB1Pcr8Df4hEFhGJfhkZGcyePZsFCxbw+eefs2zZMnJzcwGoXbs2l156KTfddBOXXHIJF198MQ0bNqRChQoBp44/0TLEJCIJbNeuXXz44Yd8+OGHzJ49mzVr1gBQqVIlrr76avr378/VV19N8+bNqVGjRsBpE4cKhIgEYuPGjbzzzju8/fbbzJs3j9zcXCpWrEjLli3p0aMHrVq14pJLLtHQUIBUIEQkYjZv3szrr7/O5MmT+eqrrwBo2rQp/fr14z/+4z/4+c9/TqlS+liKFnonRKREZWdnM23aNMaNG8esWbPIy8vjF7/4BcOGDaNt27Y0atQo6IhyHCoQIlIiMjIyGDFiBK+++ip79uyhbt269O/fn/vuu09FIUaoQIhIsVq8eDHPPfcc6enp5OXlceedd9KtWzduuOGGqDtgTE5MBUJEisUnn3zC//zP/zB37lwqVqxIr169eOihh6hfv37Q0eQ0qUCIyBlZunQpjz/+ODNmzKBmzZr85S9/oWvXrlSuXDnoaHKGVCBE5LSsX7+egQMHMmHCBCpXrsyzzz7LQw89xFlnnRV0NCkmKhAickpycnIYNmwYTz75JElJSfz3f/83jz32GFWrVg06mhQzFQgRKbL58+fTvXt3li9fTrt27Xj++eepU6dO0LGkhOgQRRE5qT179tC9e3euu+46srKymD59Ounp6SoOcU49CBE5ofnz59OxY0c2bdpE3759GTRokE6MlyDUgxCRsPLy8njmmWdo2bIlKSkpLFiwgOHDh6s4JBD1IETkGFu3buX3v/89H3zwAe3bt+fll1/WbqsJSAVCRI4yd+5c2rdvT1ZWFq+88gpdunTBLNxl4yXeaYhJRI4YNWoUbdq04eyzz2bhwoV07dpVxSGBqUCICLm5ufTq1YsePXrw61//ms8++4ymTZsGHUsCpgIhkuAOHjzIXXfdxciRI+nbty/Tpk3T9gYBtA1CJKHt2rWL2267jQULFjBixAh69eoVdCSJIioQIglq+/bttGnThm+//ZYpU6Zw5513Bh1JoowKhEgC2rRpEzfeeCMbN25kxowZtGnTJuhIEoVUIEQSzObNm7nhhhvYunUrs2bN4vrrrw86kkQpFQiRBLJt2zZat27Nli1beP/997nmmmuCjiRRTAVCJEHs2rWL1q1b8/333/Pee++pOMhJqUCIJID9+/fzm9/8hjVr1jBz5kwNK0mRqECIxLlDhw7Rrl07Pv/8c6ZMmUKrVq2CjiQxQgVCJI65Ow8++CAzZ85k9OjR3HHHHUFHkhgSsSOpzexmM1tlZmvNrF+Y+881szlm9pWZLTWzWyOVTSReDRkyhLFjxzJw4EAeeOCBoONIjIlIgTCzZOBF4BagCdDRzJoUavYnIN3drwA6AC9FIptIvEpPT2fAgAHcc889pKWlBR1HYlCkehDNgbXuvs7dc4DJQNtCbRyoFJquDGyOUDaRuLNkyRI6d+5MixYtGDNmjM7IKqclUgWiDpBRYD4ztKygNOBeM8sE3gXCnhTGzLqZ2SIzW7Rjx46SyCoS03bv3s0dd9xB1apVmTp1KmXLlg06ksSoSBWIcF9fvNB8R2Ccu9cFbgX+aWbH5HP30e7ezN2bVa9evQSiisSuvLw8OnXqREZGBlOmTOGcc84JOpLEsEjtxZQJ1CswX5djh5C6ADcDuPsCMysLVAO2RyShSBx46qmnmDFjBiNHjuTaa68NOo7EuEj1IBYCjczsfDMrTf5G6GmF2nwPtAYws4uAsoDGkESKaObMmaSlpfH73/+e//qv/wo6jsSBiBQIdz8M9ARmASvJ31tpuZkNNrPbQ80eBh4wsyXAJKCzuxcehhKRMNatW8c999zDpZdeyqhRo7RRWopFxA6Uc/d3yd/4XHDZwALTK4AWkcojEi8OHTpEx44dAXjrrbc466yzAk4k8UJHUovEuCeffJIvvviCKVOmcMEFFwQdR+KIrkktEsM+/fRTnn76ae677z7uuuuuoONInFGBEIlRe/fu5d577+W8885jxIgRQceROKQhJpEY1bt3bzZu3Mi8efOoVKnSyR8gcorUgxCJQW+++Sbjxo1jwIABtGihfTukZKhAiMSYTZs20a1bN5o1a8bAgQNP/gCR06QCIRJD3J0ePXqQnZ3NhAkTSElJCTqSxDFtgxCJIW+++SbTp09n2LBhNG7cOOg4EufUgxCJEbt27aJnz55cddVV9OnTJ+g4kgDUgxCJEY888gg7d+5k1qxZlCqlf10peepBiMSADz74gNdee41HH32Uyy+/POg4kiBUIESi3IEDB+jWrRuNGjXSXksSUeqnikS5J554gvXr1zN37lzKlSsXdBxJIOpBiESxr7/+mueee45u3brRsmXLoONIglGBEIlS7s4f/vAHUlNTGTJkSNBxJAFpiEkkSv3zn//k008/ZcyYMVStWjXoOJKA1IMQiUL//ve/eeyxx7j66qvp3Llz0HEkQakHIRKF0tLS2L59OzNmzCApSd/jJBj6yxOJMt988w1/+9vfjpyQTyQoKhAiUcTd6dWrF5UrV+bpp58OOo4kOA0xiUSRN954g7lz5zJq1ChSU1ODjiMJTj0IkSixf/9+HnnkEa688kq6du0adBwR9SBEosXw4cPZtGkTb7zxBsnJyUHHEVEPQiQabN68maFDh9KuXTtdQlSihgqESBT405/+xOHDh3XEtEQVFQiRgH311VeMGzeO3r17c8EFFwQdR+QIFQiRALk7Dz/8MKmpqQwYMCDoOCJHiViBMLObzWyVma01s37HaXO3ma0ws+VmNjFS2USCMn36dObMmUNaWhpVqlQJOo7IUSKyF5OZJQMvAm2ATGChmU1z9xUF2jQC+gMt3H23mdWIRDaRoOTk5PDII49w4YUX0q1bt6DjiBwjUru5NgfWuvs6ADObDLQFVhRo8wDworvvBnD37RHKJhKIUaNGsWbNGmbMmEFKSkrQcUSOEakhpjpARoH5zNCyghoDjc1svpl9ZmY3h3siM+tmZovMbNGOHTtKKK5Iydq7dy+DBw+mdevW3HrrrUHHEQkrUgXCwizzQvOlgEbAr4COwKtmdsygrLuPdvdm7t6sevXqxR5UJBL+8pe/sHPnToYMGYJZuH8PkeBFqkBkAvUKzNcFNodp8467H3L39cAq8guGSFzZtm0bzz33HO3atdPZWiWqRapALAQamdn5ZlYa6ABMK9TmbeAGADOrRv6Q07oI5ROJmKeeeors7GyeeuqpoKOInFBECoS7HwZ6ArOAlUC6uy83s8Fmdnuo2Sxgp5mtAOYAj7r7zkjkE4mUdevW8fLLL9O1a1caN24cdByRE4rYyfrc/V3g3ULLBhaYdqBv6CYSlwYOHEipUqUYOHDgyRuLBExHUotEyJIlS5g4cSK9e/emdu3aQccROSkVCJEI6d+/P5UrV+axxx4LOopIkeh6ECIR8NFHHzFz5kyGDh1K1apVg44jUiTqQYiUMHenX79+1K5dm169egUdR6TITrkHYWblgWx3zy2BPCJxZ9q0aXz22WeMHj2acuXKBR1HpMhO2oMwsyQzu8fM/mVm24FvgS2hM64OC51kT0TCyM3NZcCAATRu3Jj7778/6Dgip6QoPYg5wAfkn2n1G3fPAzCzs8k/sG2Imf2vu79ecjFFYtP48eNZsWIFU6ZMoVQpbfKT2FKUv9gb3f1Q4YXuvguYCkw1M52KUqSQ7OxsnnjiCZo1a8add94ZdByRU3bSIaafioOZ/dWOc1axcAVEJNG99NJLZGRk6IR8ErNOZS+mfcC00EZqzOzXZja/ZGKJxLbs7GyeeeYZ2rRpQ+vWrYOOI3Jaijwo6u5/MrN7gLlm9iOwHwh76VCRRPfpp5+yc+dOnn322aCjiJy2IhcIM2tN/lXf9gO1gC7uvqqkgonEqm3btrFgwQLuvvturrrqqqDjiJy25LS0tCI1HDRo0BhgoLsPHjRo0EfAxEGDBq1OS0tbX5IBT2T48OFpO3fuJDU1lfT0dGbOnEnjxo0ZPnw4O3fuJCMjg9dff51atWrxj3/8g7lz51KvXj3++te/kpWVxerVq5k0aRL169fnpZde4osvviA1NZWRI0fy448/smTJEtLT048857JlyzjrrLP4+9//DuR/S5w6deqR+1etWoWZMXr0aEqXLs2HH37I22+/feT+DRs2cODAAcaMGUPFihWZPn06M2bMOHL/5s2b2blzJ+PGjdNriuHX9OKLL7J69WpuuOEGGjRoEBevKR7fpwR/TVvS0tJGn+wz1vJPonrqzKwWMNXdrz2tJygGzZo180WLFgW1epFjfPfdd1x44YVcdtll6G9TopWZLXb3k16tqigHyh1vz6UtQOsTtRFJNAMHDiQlJYWWLVsGHUXkjBVlL6Y5ZtbLzM4tuDB0ZbhrzOwfwH0lkk4khnz99ddMnDiRPn36ULFixaDjiJyxohSIm4FcYJKZbTGzFWa2HlgDdASed/dxJZhRJCb079+fqlWr6nTeEjdOuheTu2cDL5lZdeBZIBU46O57SjqcSKyYM2cO7733Hn/+85+pUqVK0HFEisWpnBxmIHAWcDbwpZlNUpEQ+f+n865bty49e/YMOo5IsTnV60FkA7OAesACM7u8+COJxJa33nqLL774gkGDBul03hJXTqUH8a27PxGaftPMxgGjgFbFnkokRhw+fJgBAwbQpEkTOnXqFHQckWJ1KgXiBzO7yt0XA7j76tB2CZGENXbsWFavXs0777yj03lL3DmVv+iHgMlmthhYBlwKBHYUtUjQDhw4QFpaGi1atOC2224LOo5IsTuVk/UtCW1zuBG4mPwLCU0qqWAi0e6FF15gy5YtpKen63TeEpdOqU/s7j8C/wrdRBLWzp07GTp0KLfddhvXXXdd0HFESsSp7sUkIsCzzz5LVlYWzzzzTNBRREqMCoTIKfr+++8ZOXIknTp14uKLLw46jkiJiViBMLObzWyVma01s+NeaMjM7jIzN7OTnmlQJAhPPJG/t/egQYMCTiJSsiJSIMwsGXgRuAVoAnQ0syZh2lUkf2+pzyORS+RUffPNN4wfP56ePXty7rnnnvwBIjEsUj2I5sBad1/n7jnAZKBtmHZPAn8m/4htkagzYMAAKlasSP/+/YOOIlLiIlUg6gAZBeYzQ8uOMLMrgHruPuNET2Rm3cxskZkt2rFjR/EnFTmO2bNnM336dPr160dqamrQcURKXKQKRLidxI9cys7MkoDngYdP9kTuPtrdm7l7s+rVdSC3REZubi4PP/ww5513Hn369Ak6jkhEROrcAJnkn+DvJ3WBzQXmK5J/8N3c0AFHNYFpZna7u+u6jRK48ePH8/XXXzNp0iTKli0bdByRiIhUD2Ih0MjMzg9dia4DMO2nO9393+5ezd3ru3t94DNAxUGiwr59+3j88cf5xS9+Qfv27YOOIxIxEelBuPthM+tJ/qnCk4Gx7r7czAYDi9x92omfQSQ4w4YNY8uWLUydOlWn1JCEErHTT7r7u8C7hZYNPE7bX0Uik8jJZGZmMmzYMNq3b88111wTdByRiNKR1CIn8Pjjj5OXl8eQIUOCjiIScSoQIsexePFixo8fT58+fahfv37QcUQiTgVCJAx3p0+fPlSrVk0HxUnC0iWwRMKYMGECn3zyCa+++iqVK1cOOo5IINSDEClk7969PProozRv3pz7778/6DgigVEPQqSQQYMGsW3bNqZNm0ZSkr5DSeLSX79IAcuXL+eFF16ga9eu/PznPw86jkigVCBEQtydXr16UalSJV0pTgQNMYkcMWXKFObMmcNLL71EtWrVgo4jEjj1IETIP99S3759ueKKK+jWrVvQcUSignoQIkBaWhqbNm1iypQpJCcnBx1HJCqoByEJb+HChTz//PN0795d51sSKUAFQhLaoUOH6Nq1KzVr1mTo0KFBxxGJKhpikoQ2bNgwli5dyttvv60jpkUKUQ9CEtaqVasYPHgwd911F23btg06jkjUUYGQhJSXl8cDDzxAuXLl+Nvf/hZ0HJGopCEmSUijR4/m448/ZsyYMdSsWTPoOCJRST0ISTgZGRk89thjtGrVSifjEzkBFQhJKHl5eXTq1Im8vDxGjx6ta0yLnICGmCShPPfcc8ydO5cxY8bQoEGDoOOIRDX1ICRhfP311wwYMIDf/va3GloSKQIVCEkIBw8e5He/+x3VqlXT0JJIEWmISRJCv379WLFiBe+9957O1CpSROpBSNybNWsWI0aM4KGHHuKmm24KOo5IzFCBkLi2detWOnfuTNOmTRkyZEjQcURiioaYJG4dOnSIu+++m7179/L+++9Trly5oCOJxBQVCIlbjz76KB9//DETJ07kkksuCTqOSMyJ2BCTmd1sZqvMbK2Z9Qtzf18zW2FmS83sQzM7L1LZJP5MmjSJF154gd69e9OxY8eg44jEpIgUCDNLBl4EbgGaAB3NrEmhZl8Bzdz9UuBN4M+RyCbxZ9myZXTt2pXrrruOYcOGBR1HJGZFqgfRHFjr7uvcPQeYDBx1fmV3n+PuB0KznwF1I5RN4siePXu44447qFSpEunp6aSkpAQdSSRmRWobRB0go8B8JnD1Cdp3AWaGu8PMugHdAM4999ziyidxIDc3l3vvvZcNGzYwd+5catWqFXQkkZgWqR5EuMNWPWxDs3uBZkDYsQF3H+3uzdy9WfXq1YsxosS6vn378q9//YsRI0bQokWLoOOIxLxI9SAygXoF5usCmws3MrMbgceBlu7+Y4SySRwYOXIkI0aM4I9//CM9evQIOo5IXIhUD2Ih0MjMzjez0kAHYFrBBmZ2BfAycLu7b49QLokDb7/9Nr1796Zt27baKC1SjCJSINz9MNATmAWsBNLdfbmZDTaz20PNhgEVgClm9rWZTTvO04kc8dFHH9GhQweaN2/OhAkTSE5ODjqSSNyI2IFy7v4u8G6hZQMLTN8YqSwSH7766ituv/12LrjgAmbMmEH58uWDjiQSV3QuJolJ33zzDW3atKFy5crMmjWL1NTUoCOJxB0VCIk5K1eupHXr1pQpU4bZs2dTr169kz9IRE6ZCoTElOXLl9OqVSuSkpKYPXs2DRs2DDqSSNxSgZCY8eWXX9KyZUvMjNmzZ/Ozn/0s6EgicU0FQmLCJ598QqtWrahQoQIff/wxF110UdCRROKeCoREvalTp3LjjTdSs2ZN5s2bR4MGDYKOJJIQVCAkark7L7zwAu3atePKK69k/vz5Ov+WSASpQEhUysnJoXv37vTp04e2bdvy4YcfaldWkQhTgZCos3XrVtq0acMrr7zCgAEDmDp1qi4XKhIAXXJUosq8efPo0KEDe/bsYcKECdxzzz1BRxJJWOpBSFQ4fPgwTz311JE9lT7//HMVB5GAqQchgfv222/p1KkTCxcupGPHjowaNYpKlSoFHUsk4akHIYHJy8vj+eef54orrmDdunWkp6czceJEFQeRKKEehARi5cqVPPjgg8ybN4/bbruN0aNHU7NmzaBjiUgB6kFIRO3evZs+ffpwySWXsGTJEl577TXeeecdFQeRKKQCIRGRm5vLyy+/TOPGjRkxYgRdunRhzZo1dO7cGbNwlywXkaCpQEiJysvLY+rUqVxxxRU8+OCDNGnShC+//JKXX36Z6tWrBx1PRE5ABUJKRG5uLpMmTeLSSy/lrrvuIjs7mzfeeIO5c+dy+eWXBx1PRIpABUKK1YEDBxg7dixNmjThnnvuwd2ZOHEiK1eu5O6779ZwkkgM0V5MUixWrVrFqFGjGDduHHv27OGyyy7jzTff5Le//S1JSfoeIhKLVCDktGVlZTF9+nTGjBnD7NmzSUlJ4Y477qBHjx788pe/VG9BJMapQMgpOXjwIO+++y6TJ09mxowZZGdnc9555/H000/TpUsXzjnnnKAjikgxUYGQk9q0aRMzZ85k5syZvP/+++zbt48aNWrQpUsXOnTowLXXXqthJJE4pAIhx9izZw+ffvopH330ETNnzmTZsmUA1KlTh44dO3L33Xfzq1/9ilKl9OcjEs/0H57g8vLy+O6771i0aBHz58/n448/ZtmyZbg7pUqV4rrrrmPo0KHccsstXHzxxdquIJJAVCASSFZWFitXrmTFihUsWbKEL7/8kq+++oqsrCwAypcvzzXXXENaWhrXX389zZs3p3z58gGnFpGgqEDEEXdn9+7drF+/nvXr17NhwwbWr1/P2rVrWbFiBZmZmUfalitXjssvv5xOnTpx5ZVXcuWVV3LxxRdr2EhEjtCnQQxwdw7a04LEAAAFQ0lEQVQcOMD27dvZunUrW7ZsOeZnRkYGGzZsYO/evUc9tkqVKjRo0IAbbriBiy66iCZNmnDRRRdxwQUXqBiIyAlF7BPCzG4GXgCSgVfdfUih+8sA44GrgJ1Ae3ffEKl8JcHdycnJYf/+/ezbt4/9+/cfuRWcz8rKYvfu3ezateuonwWnc3Jyjnn+pKQkatSoQc2aNalXrx6//OUvOf/884/c6tevT5UqVQJ45SISDyJSIMwsGXgRaANkAgvNbJq7ryjQrAuw290bmlkHYCjQ/kTPu2/fPmbPnk1ubm6x3w4fPszhw4f58ccfycnJOenPgtPZ2dlHikBubm6Rf0+VKlWiatWqnH322VStWpWmTZtStWrVI7caNWpQq1YtatasSa1atahevTrJycmn/H6IiBRFpHoQzYG17r4OwMwmA22BggWiLZAWmn4TGGlm5u5+vCddtWoVrVu3LpHASUlJpKSkUKZMGcqUKUPp0qXD/ixTpgyVKlU6ar5MmTKUL1/+qFuFChWOu6xChQpUqVJFQz4iElUi9YlUB8goMJ8JXH28Nu5+2Mz+DaQCPxRsZGbdgG4ANWrUID09neTk5GK9JSUlaXdOEUl4kSoQ4T5tC/cMitIGdx8NjAZo1qyZt2zZ8szTiYjIMSJ1foRMoF6B+brA5uO1MbNSQGVgV0TSiYjIMSJVIBYCjczsfDMrDXQAphVqMw24LzR9FzD7RNsfRESkZEVkiCm0TaEnMIv83VzHuvtyMxsMLHL3acAY4J9mtpb8nkOHSGQTEZHwIrbbjLu/C7xbaNnAAtPZQLtI5RERkRPTOZpFRCQsFQgREQlLBUJERMJSgRARkbBUIEREJCwVCBERCUsFQkREwlKBEBGRsFQgREQkLBUIEREJSwVCRETCUoEQEZGwLJbPqG1mO4CNQec4A9UodMU8iRt6b+NbrL+/57l79ZM1iukCEevMbJG7Nws6hxQ/vbfxLVHeXw0xiYhIWCoQIiISlgpEsEYHHUBKjN7b+JYQ76+2QYiISFjqQYiISFgqECIiEpYKRMDMLM3MNpnZ16HbrUFnkjNjZjeb2SozW2tm/YLOI8XLzDaY2bLQ/+uioPOUJG2DCJiZpQH73P0vQWeRM2dmycBqoA2QCSwEOrr7ikCDSbExsw1AM3eP5QPlikQ9CJHi1RxY6+7r3D0HmAy0DTiTyGlRgYgOPc1sqZmNNbOqQYeRM1IHyCgwnxlaJvHDgffNbLGZdQs6TElSgYgAM/vAzL4Jc2sL/B1oAFwObAGGBxpWzpSFWaZx3PjSwt2vBG4B/mBmvww6UEkpFXSARODuNxalnZm9Aswo4ThSsjKBegXm6wKbA8oiJcDdN4d+bjez/yV/WHFesKlKhnoQATOzWgVmfwt8E1QWKRYLgUZmdr6ZlQY6ANMCziTFxMzKm1nFn6aBXxPH/7PqQQTvz2Z2OfnDEBuA7sHGkTPh7ofNrCcwC0gGxrr78oBjSfE5B/hfM4P8z8+J7v5esJFKjnZzFRGRsDTEJCIiYalAiIhIWCoQIiISlgqEiIiEpQIhIiJhqUCIiEhYKhAiIhKWCoRIMTKzOWbWJjT9lJmNCDqTyOnSkdQixesJYLCZ1QCuAG4POI/IadOR1CLFzMw+AioAv3L3rKDziJwuDTGJFCMzuwSoBfyo4iCxTgVCpJiEzsw7gfwryO03s5sCjiRyRlQgRIqBmZ0FvAU87O4rgSeBtEBDiZwhbYMQEZGw1IMQEZGwVCBERCQsFQgREQlLBUJERMJSgRARkbBUIEREJCwVCBERCev/AFTFKsA7kbA6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = numpy.arange(-7,7.1,0.1)\n",
    "Y = [sigmoid(x) for x in X]\n",
    "plt.plot(X, Y, color='black')\n",
    "plt.plot([0,0],[-2,2], color = 'black', linewidth=0.5)\n",
    "plt.plot([-7,7],[0.5,0.5], color = 'k', linewidth=0.5, linestyle='--')\n",
    "plt.xlim(-7, 7)\n",
    "plt.ylim(-0.1, 1.1)\n",
    "plt.xticks([-5,0,5])\n",
    "plt.xlabel(\"$x$\")\n",
    "plt.ylabel(r\"$\\sigma(x)$\")\n",
    "plt.title(\"Sigmoid function\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementing a simple classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read a small dataset of beer reviews"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = dataDir + \"beer_50000.json\"\n",
    "f = open(path)\n",
    "\n",
    "data = []\n",
    "\n",
    "for l in f:\n",
    "    if 'user/gender' in l: # Discard users who didn't specify gender\n",
    "        d = eval(l) # demjson is more secure but much slower!\n",
    "        data.append(d)\n",
    "    \n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'review/appearance': 4.0,\n",
       " 'beer/style': 'American Double / Imperial IPA',\n",
       " 'review/palate': 4.0,\n",
       " 'review/taste': 4.5,\n",
       " 'beer/name': 'Cauldron DIPA',\n",
       " 'review/timeUnix': 1293735206,\n",
       " 'user/gender': 'Male',\n",
       " 'user/birthdayRaw': 'Jun 16, 1901',\n",
       " 'beer/ABV': 7.7,\n",
       " 'beer/beerId': '64883',\n",
       " 'user/birthdayUnix': -2163081600,\n",
       " 'beer/brewerId': '1075',\n",
       " 'review/timeStruct': {'isdst': 0,\n",
       "  'mday': 30,\n",
       "  'hour': 18,\n",
       "  'min': 53,\n",
       "  'sec': 26,\n",
       "  'mon': 12,\n",
       "  'year': 2010,\n",
       "  'yday': 364,\n",
       "  'wday': 3},\n",
       " 'user/ageInSeconds': 3581417047,\n",
       " 'review/overall': 4.0,\n",
       " 'review/text': \"According to the website, the style for the Caldera Cauldron changes every year. The current release is a DIPA, which frankly is the only cauldron I'm familiar with (it was an IPA/DIPA the last time I ordered a cauldron at the horsebrass several years back). In any event... at the Horse Brass yesterday.\\t\\tThe beer pours an orange copper color with good head retention and lacing. The nose is all hoppy IPA goodness, showcasing a huge aroma of dry citrus, pine and sandlewood. The flavor profile replicates the nose pretty closely in this West Coast all the way DIPA. This DIPA is not for the faint of heart and is a bit much even for a hophead like myslf. The finish is quite dry and hoppy, and there's barely enough sweet malt to balance and hold up the avalanche of hoppy bitterness in this beer. Mouthfeel is actually fairly light, with a long, persistentely bitter finish. Drinkability is good, with the alcohol barely noticeable in this well crafted beer. Still, this beer is so hugely hoppy/bitter, it's really hard for me to imagine ordering more than a single glass. Regardless, this is a very impressive beer from the folks at Caldera.\",\n",
       " 'user/profileName': 'johnmichaelsen',\n",
       " 'review/aroma': 4.5}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predict the user's gender from the length of their review"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [[1, len(d['review/text'])] for d in data]\n",
    "y = [d['user/gender'] == 'Female' for d in data]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegression()"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = sklearn.linear_model.LogisticRegression(fit_intercept=False)\n",
    "mod.fit(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the accuracy of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9849041807577317"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predictions = mod.predict(X) # Binary vector of predictions\n",
    "correct = predictions == y # Binary vector indicating which predictions were correct\n",
    "sum(correct) / len(correct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accuracy seems surprisingly high! Check against the number of positive labels..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9849041807577317"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 - (sum(y) / len(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accuracy is identical to the proportion of \"males\" in the data. Confirm that the model is never predicting positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing a balanced classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the class_weight='balanced' option to implement the balanced classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = sklearn.linear_model.LogisticRegression(class_weight='balanced')\n",
    "mod.fit(X,y)\n",
    "predictions = mod.predict(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple classification diagnostics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the accuracy of the balanced model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4225849139832378"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "correct = predictions == y\n",
    "sum(correct) / len(correct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### True positives, False positives (etc.), and balanced error rate (BER)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "TP = sum([(p and l) for (p,l) in zip(predictions, y)])\n",
    "FP = sum([(p and not l) for (p,l) in zip(predictions, y)])\n",
    "TN = sum([(not p and not l) for (p,l) in zip(predictions, y)])\n",
    "FN = sum([(not p and l) for (p,l) in zip(predictions, y)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TP = 199\n",
      "FP = 11672\n",
      "TN = 8423\n",
      "FN = 109\n"
     ]
    }
   ],
   "source": [
    "print(\"TP = \" + str(TP))\n",
    "print(\"FP = \" + str(FP))\n",
    "print(\"TN = \" + str(TN))\n",
    "print(\"FN = \" + str(FN))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can rewrite the accuracy in terms of these metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4225849139832378"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(TP + TN) / (TP + FP + TN + FN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### True positive and true negative rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPR = TP / (TP + FN)\n",
    "TNR = TN / (TN + FP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.6461038961038961, 0.4191589947748196)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TPR,TNR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Balanced error rate (BER)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4673685545606422"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BER = 1 - 1/2 * (TPR + TNR)\n",
    "BER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precision, recall, and F1 scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "precision = TP / (TP + FP)\n",
    "recall = TP / (TP + FN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0167635414034201, 0.6461038961038961)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "precision, recall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F1 score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.032679201904918305"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F1 = 2 * (precision*recall) / (precision + recall)\n",
    "F1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Significance testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = dataDir + \"beer_500.json\"\n",
    "f = open(path)\n",
    "\n",
    "data = []\n",
    "\n",
    "for l in f:\n",
    "    d = eval(l)\n",
    "    data.append(d)\n",
    "    \n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Randomly sort the data (so that train and test are iid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(0)\n",
    "random.shuffle(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predict overall rating from ABV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = [[1] for d in data] # Model *without* the feature\n",
    "X2 = [[1, d['beer/ABV']] for d in data] # Model *with* the feature\n",
    "y = [d['review/overall'] for d in data]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the two models (with and without the feature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "model1 = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model1.fit(X1[:250], y[:250]) # Train on first half\n",
    "residuals1 = model1.predict(X1[250:]) - y[250:] # Test on second half"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "model2 = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model2.fit(X2[:250], y[:250])\n",
    "residuals2 = model2.predict(X2[250:]) - y[250:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Residual sum of squares for both models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "rss1 = sum([r**2 for r in residuals1])\n",
    "rss2 = sum([r**2 for r in residuals2])\n",
    "k1,k2 = 1,2 # Number of parameters of each model\n",
    "n = len(residuals1) # Number of samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F statistic (results may vary for different random splits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = ((rss1 - rss2) / (k2 - k1)) / (rss2 / (n-k2))\n",
    "1 - scipy.stats.f.cdf(F,k2-k1,n-k2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression in tensorflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Small dataset of fantasy reviews"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = dataDir + \"fantasy_100.json\"\n",
    "f = open(path)\n",
    "\n",
    "data = []\n",
    "\n",
    "for l in f:\n",
    "    d = json.loads(l)\n",
    "    data.append(d)\n",
    "    \n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predict rating from review length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratings = [d['rating'] for d in data]\n",
    "lengths = [len(d['review_text']) for d in data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = numpy.matrix([[1,l] for l in lengths])\n",
    "y = numpy.matrix(ratings).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First check the coefficients if we fit the model using sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/dist-packages/sklearn/utils/validation.py:590: FutureWarning: np.matrix usage is deprecated in 1.0 and will raise a TypeError in 1.2. Please convert to a numpy array with np.asarray. For more information see: https://numpy.org/doc/stable/reference/generated/numpy.matrix.html\n",
      "  FutureWarning,\n",
      "/usr/local/lib/python3.7/dist-packages/sklearn/utils/validation.py:590: FutureWarning: np.matrix usage is deprecated in 1.0 and will raise a TypeError in 1.2. Please convert to a numpy array with np.asarray. For more information see: https://numpy.org/doc/stable/reference/generated/numpy.matrix.html\n",
      "  FutureWarning,\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[3.98394783e+00, 1.19363599e-04]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)\n",
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert features and labels to tensorflow structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2021-12-07 23:22:08.609713: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:08.609805: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcublas.so.11'; dlerror: libcublas.so.11: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:08.609867: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcublasLt.so.11'; dlerror: libcublasLt.so.11: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:08.609927: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcufft.so.10'; dlerror: libcufft.so.10: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:08.609988: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcurand.so.10'; dlerror: libcurand.so.10: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:08.610046: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcusolver.so.11'; dlerror: libcusolver.so.11: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:08.610105: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcusparse.so.11'; dlerror: libcusparse.so.11: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:08.610163: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudnn.so.8'; dlerror: libcudnn.so.8: cannot open shared object file: No such file or directory\n",
      "2021-12-07 23:22:08.610176: W tensorflow/core/common_runtime/gpu/gpu_device.cc:1850] Cannot dlopen some GPU libraries. Please make sure the missing libraries mentioned above are installed properly if you would like to use GPU. Follow the guide at https://www.tensorflow.org/install/gpu for how to download and setup the required libraries for your platform.\n",
      "Skipping registering GPU devices...\n"
     ]
    }
   ],
   "source": [
    "X = tf.constant(X, dtype=tf.float32)\n",
    "y = tf.constant(y, dtype=tf.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build tensorflow regression class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "class regressionModel(tf.keras.Model):\n",
    "    def __init__(self, M, lamb):\n",
    "        super(regressionModel, self).__init__()\n",
    "        # Initialize weights to zero\n",
    "        self.theta = tf.Variable(tf.constant([0.0]*M, shape=[M,1], dtype=tf.float32))\n",
    "        self.lamb = lamb\n",
    "\n",
    "    # Prediction (for a matrix of instances)\n",
    "    def predict(self, X):\n",
    "        return tf.matmul(X, self.theta)\n",
    "\n",
    "    # Mean Squared Error\n",
    "    def MSE(self, X, y):\n",
    "        return tf.reduce_mean((tf.matmul(X, self.theta) - y)**2)\n",
    "\n",
    "    # Regularizer\n",
    "    def reg(self):\n",
    "        return self.lamb * tf.reduce_sum(self.theta**2)\n",
    "    \n",
    "    # L1 regularizer\n",
    "    def reg1(self):\n",
    "        return self.lamb * tf.reduce_sum(tf.abs(self.theta))\n",
    "\n",
    "    # Loss\n",
    "    def call(self, X, y):\n",
    "        return self.MSE(X, y) + self.reg()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the model (lambda = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimizer = tf.keras.optimizers.Adam(0.1)\n",
    "model = regressionModel(len(X[0]), 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train for 1000 iterations of gradient descent (could implement more careful stopping criteria)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "for iteration in range(1000):\n",
    "    with tf.GradientTape() as tape:\n",
    "        loss = model(X, y)\n",
    "    gradients = tape.gradient(loss, model.trainable_variables)\n",
    "    optimizer.apply_gradients(zip(gradients, model.trainable_variables))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confirm that we get a similar result to what we got using sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Variable 'Variable:0' shape=(2, 1) dtype=float32, numpy=\n",
       "array([[3.9834061e+00],\n",
       "       [1.1961416e-04]], dtype=float32)>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a few predictions using the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(10, 1), dtype=float32, numpy=\n",
       "array([[4.232921 ],\n",
       "       [4.165339 ],\n",
       "       [4.1651   ],\n",
       "       [4.197635 ],\n",
       "       [4.194166 ],\n",
       "       [4.0396247],\n",
       "       [4.0818486],\n",
       "       [4.047041 ],\n",
       "       [4.0570884],\n",
       "       [4.0489545]], dtype=float32)>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.predict(X)[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification in Tensorflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predict whether rating is above 4 from length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = numpy.matrix([[1,l*0.0001] for l in lengths]) # Rescale the lengths for conditioning\n",
    "y_class = numpy.matrix([[r > 4 for r in ratings]]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert to tensorflow structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = tf.constant(X, dtype=tf.float32)\n",
    "y_class = tf.constant(y_class, dtype=tf.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tensorflow classification class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "class classificationModel(tf.keras.Model):\n",
    "    def __init__(self, M, lamb):\n",
    "        super(classificationModel, self).__init__()\n",
    "        self.theta = tf.Variable(tf.constant([0.0]*M, shape=[M,1], dtype=tf.float32))\n",
    "        self.lamb = lamb\n",
    "\n",
    "    # Probability (for a matrix of instances)\n",
    "    def predict(self, X):\n",
    "        return tf.math.sigmoid(tf.matmul(X, self.theta))\n",
    "\n",
    "    # Objective\n",
    "    def obj(self, X, y):\n",
    "        pred = self.predict(X)\n",
    "        pos = y*tf.math.log(pred)\n",
    "        neg = (1.0 - y)*tf.math.log(1.0 - pred)\n",
    "        return -tf.reduce_mean(pos + neg)\n",
    "    \n",
    "    # Same objective, using tensorflow short-hand\n",
    "    def obj_short(self, X, y):\n",
    "        pred = self.predict(X)\n",
    "        bce = tf.keras.losses.BinaryCrossentropy()\n",
    "        return tf.reduce_mean(bce(y, pred))\n",
    "\n",
    "    # Regularizer\n",
    "    def reg(self):\n",
    "        return self.lamb * tf.reduce_sum(self.theta**2)\n",
    "    \n",
    "    # Loss\n",
    "    def call(self, X, y):\n",
    "        return self.obj(X, y) + self.reg()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the model (lambda = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimizer = tf.keras.optimizers.Adam(0.1)\n",
    "model = classificationModel(len(X[0]), 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run for 1000 iterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "for iteration in range(1000):\n",
    "    with tf.GradientTape() as tape:\n",
    "        loss = model(X, y_class)\n",
    "    gradients = tape.gradient(loss, model.trainable_variables)\n",
    "    optimizer.apply_gradients(zip(gradients, model.trainable_variables))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Variable 'Variable:0' shape=(2, 1) dtype=float32, numpy=\n",
       "array([[-0.08591048],\n",
       "       [ 0.30873907]], dtype=float32)>"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model predictions (as probabilities via the sigmoid function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(10, 1), dtype=float32, numpy=\n",
       "array([[0.49462333],\n",
       "       [0.4902634 ],\n",
       "       [0.490248  ],\n",
       "       [0.49234676],\n",
       "       [0.49212298],\n",
       "       [0.48215765],\n",
       "       [0.4848793 ],\n",
       "       [0.4826356 ],\n",
       "       [0.4832832 ],\n",
       "       [0.48275894]], dtype=float32)>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.predict(X)[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regularization pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parseData(fname):\n",
    "    for l in open(fname):\n",
    "        yield eval(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just read the first 5000 reviews (deliberately making a model that will overfit if not carefully regularized)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = list(parseData(dataDir + \"beer_50000.json\"))[:5000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a simple bag-of-words model (see Chapter 8 for more details)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "wordCount = defaultdict(int)\n",
    "punctuation = set(string.punctuation)\n",
    "for d in data: # Strictly, should just use the *training* data to extract word counts\n",
    "    r = ''.join([c for c in d['review/text'].lower() if not c in punctuation])\n",
    "    for w in r.split():\n",
    "        wordCount[w] += 1\n",
    "\n",
    "counts = [(wordCount[w], w) for w in wordCount]\n",
    "counts.sort()\n",
    "counts.reverse()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1000 most popular words"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "words = [x[1] for x in counts[:1000]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "wordId = dict(zip(words, range(len(words))))\n",
    "wordSet = set(words)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bag-of-words features for 1000 most popular words"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feature(datum):\n",
    "    feat = [0]*len(words)\n",
    "    r = ''.join([c for c in datum['review/text'].lower() if not c in punctuation])\n",
    "    for w in r.split():\n",
    "        if w in words:\n",
    "            feat[wordId[w]] += 1\n",
    "    feat.append(1) # offset\n",
    "    return feat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.shuffle(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [feature(d) for d in data]\n",
    "y = [d['review/overall'] for d in data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ntrain,Nvalid,Ntest = 4000,500,500\n",
    "Xtrain,Xvalid,Xtest = X[:Ntrain],X[Ntrain:Ntrain+Nvalid],X[Ntrain+Nvalid:]\n",
    "ytrain,yvalid,ytest = y[:Ntrain],y[Ntrain:Ntrain+Nvalid],y[Ntrain+Nvalid:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unregularized model (train on training set, test on test set)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(Xtrain, ytrain)\n",
    "predictions = model.predict(Xtest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5621377319894529"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum((ytest - predictions)**2)/len(ytest) # Mean squared error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regularized model (\"ridge regression\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = linear_model.Ridge(1.0, fit_intercept=False) # MSE + 1.0 l2\n",
    "model.fit(Xtrain, ytrain)\n",
    "predictions = model.predict(Xtest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5553767774281313"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum((ytest - predictions)**2)/len(ytest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Complete regularization pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Track the model which works best on the validation set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestModel = None\n",
    "bestVal = None\n",
    "bestLamb = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train models for different values of lambda (or C). Keep track of the best model on the validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "l = 0.01, validation MSE = 0.4933728029983656\n",
      "l = 0.1, validation MSE = 0.49292702478338096\n",
      "l = 1, validation MSE = 0.48865630138060057\n",
      "l = 10, validation MSE = 0.4585804821929264\n",
      "l = 100, validation MSE = 0.39865226713206803\n",
      "l = 1000, validation MSE = 0.413825760476879\n",
      "l = 10000, validation MSE = 0.5041631912430448\n"
     ]
    }
   ],
   "source": [
    "ls = [0.01, 0.1, 1, 10, 100, 1000, 10000]\n",
    "errorTrain = []\n",
    "errorValid = []\n",
    "\n",
    "for l in ls:\n",
    "    model = sklearn.linear_model.Ridge(l)\n",
    "    model.fit(Xtrain, ytrain)\n",
    "    predictTrain = model.predict(Xtrain)\n",
    "    MSEtrain = sum((ytrain - predictTrain)**2)/len(ytrain)\n",
    "    errorTrain.append(MSEtrain)\n",
    "    predictValid = model.predict(Xvalid)\n",
    "    MSEvalid = sum((yvalid - predictValid)**2)/len(yvalid)\n",
    "    errorValid.append(MSEvalid)\n",
    "    print(\"l = \" + str(l) + \", validation MSE = \" + str(MSEvalid))\n",
    "    if bestVal == None or MSEvalid < bestVal:\n",
    "        bestVal = MSEvalid\n",
    "        bestModel = model\n",
    "        bestLamb = l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the best model from the validation set, compute the error on the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4380927014652703"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predictTest = bestModel.predict(Xtest)\n",
    "MSEtest = sum((ytest - predictTest)**2)/len(ytest)\n",
    "MSEtest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the train/validation/test error associated with this pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEcCAYAAAA2g5hwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VNX28PHvIpCE0CEgNQQVAYEICBqqoNJEqYIgSPEiiu3+rnpFlKKo71XsXrkoeqWIiIhSLCCCIIKgFAEpQuiEICShhFBCynr/mEnuECAzKZOZJOvzPPMw55y9z1mbSWblnH3O3qKqGGOMMVkp5usAjDHG+D9LFsYYY9yyZGGMMcYtSxbGGGPcsmRhjDHGLUsWxhhj3LJkYQoFEQkXERWR4s7lRSIyxJOyOTjWsyLyUW7izcEx24rIzjza134Rud35Pt/bYgomSxbGL4jI9yIy4TLre4jIX9n9YlfVrqo6PQ/iai8i0Zn2/f9UdXhu932ZYw0VkVQRSRSRBBHZJCJ3Oo/5s6rWy+tjeqstpvCxZGH8xTTgPhGRTOvvAz5V1ZT8D8kn1qhqaaA88F9gjohU9HFMxliyMH5jPlARaJu+QkQqAHcCM5zL3UTkd+df3YdE5Pkr7UxEVojIcOf7ABF5XUTiRGQv0C1T2WEiskNETovIXhF50Lm+FLAIqO78az9RRKqLyPMiMtOlfncR2SYiJ53HbeCybb+IPCUiW0TklIh8LiLB7v4zVDUN+BgoCVyd+QzHud/RIrJdRE6IyFTX/YrInc4zk5Mi8ouIRFzh/ymjLS6X54aIyEHn/9dzLmWLicgzIrJHROJFxBJZEWLJwvgFVT0HzAEGu6zuB/ypqpudy2ec28vj+MIfKSI9Pdj9AziSTlOgOXB3pu3HnNvLAsOAt0SkmaqeAboCMapa2vmKca0oItcBnwH/B1QGvgO+FpHATO3oAtQBIoCh7gJ2XnYbDiQCUVcoNhDoDFwDXAeMcdZthiPRPAhUAj4AFopIkLvjOrUB6gG3AeNckt/jQE/gFqA6cAKY5OE+TQFnycL4k+lAXxEp6Vwe7FwHgKquUNU/VDVNVbfg+JK+xYP99gPeVtVDqnoc+JfrRlX9VlX3qMNPwBJcznDcuAf4VlV/UNVk4HUcZwOtXMq8q6oxzmN/DTTJYn+RInIS+AsYAPRS1VNXKPueS5tedpYHR3L8QFV/VdVUZ99NEhDpYZteUNVzziS9GbjBuf5B4DlVjVbVJOB54O6c3ihgChb7kI3fUNVVIhIL9BCR34AWQO/07SJyM/AK0AgIBIKALzzYdXXgkMvyAdeNItIVGI/jr/NiQAjwh4dhV3fdn6qmicghoIZLmb9c3p911rmStaraxsNjZ25T+n5rA0NE5DGX7YFujusqc7ylXfY7T0TSXLanAlcBhz3ctymg7MzC+JsZOM4o7gOWqOpRl22zgIVALVUtB7wPZO4Qv5wjQC2X5bD0N85LM1/iOCO4SlXL47iUlL5fd8Myx+D4Ek3fnziPlR9fnpnblH6J7BDwsqqWd3mFqOpnuTzeIaBrpv0Gq6oliiLAkoXxNzOA23FcSsl862sZ4LiqnheRm4B7PdznHOBxEanp7DR/xmVb+hlKLJDiPMvo5LL9KFBJRMplse9uInKbiJQAnsRxyecXD2PLjUecbaoIPAt87lz/IfCQiNwsDqWcNweUyeXx3gdeFpHaACJSWUR65HKfpoCwZGH8iqrux/FFWwrHWYSrh4EJInIaGIfji9oTHwLf47j+vhH4yuV4p3F03M7B0WF7r+txVfVPHH0je513Fl10KUdVdwKDgH8DccBdwF2qesHD2HJjFo7+lb3O10vOmNbjSLbv4WjTbjzoVPfAOzj+b5Y4P4O1wM15sF9TAIhNfmRMwSMi+4HhqrrU17GYosHOLIwxxrhlycIYY4xbdhnKGGOMW3ZmYYwxxi1LFsYYY9wqNE9wh4aGanh4uK/DMMaYAmXDhg1xqlrZXblCkyzCw8NZv369r8MwxpgCRUQOuC9ll6GMMcZ4wJKFMcYYtyxZGGOMccuShTHGGLcsWRhjjHHLkoUxxhi3LFkYY0wBdvjwYY4ePeq+YC55NVmISBcR2Skiu0XkmctsHyoisSKyyfka7rJtiIhEOV9DvBmnMcYUROfPn+eLL77gyy+/xNvj/HntoTwRCQAmAR2BaGCdiCxU1e2Zin6uqo9mqlsRx5zIzXFMa7nBWfeEt+I1xpiCZvHixSQkJDBs2DAcM/p6jzfPLG4CdqvqXuesYbMBT6dg7Az8oKrHnQniB6CLl+I0xpgCZ9u2bWzevJm2bdtSq1Yt9xVyyZvJogaOCd7TRTvXZdZHRLaIyFwRSW+xp3WNMabISUhI4JtvvqFGjRq0a9cuX47pzWRxuXOizBfVvgbCVTUCWApMz0ZdRGSEiKwXkfWxsbG5CtYYYwoCVWX+/PmkpqbSq1cvAgIC8uW43kwW0YDruVFNIMa1gKrGq2qSc/FD4EZP6zrrT1HV5qravHJlt4MmGmNMgbd27Vr27dtH586dqVSpUr4d15vJYh1QV0TqiEgg0B9Y6FpARKq5LHYHdjjffw90EpEKIlIB6ORcZ4wxRdbRo0dZtmwZ9erVo1mzZvl6bK/dDaWqKSLyKI4v+QDgY1XdJiITgPWquhB4XES6AynAcWCos+5xEXkRR8IBmKCqx70VqzHG+LuUlBS++uorgoODueuuu7x+91NmhWYO7ubNm2tO5rM4f/488+bNQ0Q8egEerfO0bm7KXa5ssWLFCA4OpmTJkhn/BgUF5fsPljEmb33//fesXbuWe++9l7p16+bZfkVkg6o2d1eu0Ex+lFOqSkJCAqp62Vd6GXev7JTLbyJCcHDwJUkk87+XWxcYGGiJxhgf27t3L2vXrqVFixZ5miiyo8gni5IlS/Lggw/m6zFzmmg8KZuWlsb58+c5d+4c586dy3if+d9Tp05llMkqiYlIlgkmq0RTokQJSzTG5NK5c+eYP38+oaGhdOzY0WdxFPlk4Quul5B8TVW5cOFClgkm/X368okTJzLWZZVoihUrluNEU7x4cb/5PzLGV1SVb775hjNnzjBgwABKlCjhs1gsWRRxIkJQUBBBQUGUK1cuW3VVlaSkJI8TzZkzZ4iPj89Yl5WAgAAqV65My5YtadSoEcWK2ZiXpujZsmUL27dv57bbbqNatWruK3iRJQuTY659IeXLl89WXVW96GzlcokmKiqKefPmsXz5clq1akXTpk0pXtx+ZE3RcPLkSb777jvCwsJo1aqVr8OxZGF8I70vpGTJklSoUOGyZW6//XZ27tzJqlWr+O6771i5ciWRkZE0b96coKCgfI7YmPyTlpaWcZdmr169/OLM2pKF8VsiQv369alXrx779+9n1apVLF26lFWrVtGiRQtuvvlmSpUq5eswjclzq1ev5uDBg/Tq1SvbZ+3eYsnC+D0RoU6dOtSpU4fDhw+zevVqfv75Z9asWUOzZs1o1apVtvtbjPFXMTExrFixgoYNG9K4cWNfh5PBkoUpUGrUqEG/fv2IjY1l9erVrF+/nvXr1xMREUHr1q0JDQ31dYjG5FhycjJfffUVpUuXplu3bh7dEfjVV1/RtWtXSpYs6dXYfH8hzJgcqFy5Mj179uSxxx6jefPmbN26lUmTJjFnzhxiYi4Zc9KYAmHJkiXEx8fTs2dPj778f/nlF/r06cO7777r9diK/HAfpnA4c+YMv/76K7/99htJSUlcffXVtGnThvDwcHtewxQIUVFRzJo1i5YtW9KpUye35dPS0oiMjOTw4cPs3LmT0qVL5+i4NtyHKVJKlSrFrbfeSuvWrVm3bh1r165lxowZ1KhRgzZt2lCvXj1LGsZvnTlzhgULFlClShVuvfVWj+rMmjWLdevWMX369BwniuywMwtTKCUnJ7Np0yZ++eUXTp48SeXKlWnTpo094Gf8jqry+eefs3v3bh544AGuuuoqt3XS0tKoX78+5cqV49dff83Vz7SnZxaWLEyhlpaWxtatW1m1ahWxsbGUL1+eVq1a0aRJE58OnWBMug0bNvDNN9/QuXNnIiMjPa536NAhEhISaNiwYa6Ob8nCGBeqyq5du1i1ahXR0dGUKlWKyMhIWrRoYQ/4GZ+Jj4/ngw8+oFatWgwaNMijS6Xnzp0jODg4zy6rWp+FMS5EhHr16nHddddx4MABVq1axbJlyzIe8IuMjLQH/Ey+Sk1NZd68eQQEBNCjRw+Pv/xHjBjByZMnWbhwYb72w1myMEWKiBAeHk54eDgxMTGsWrWKVatWsXbtWnvAz+SrlStXcvjwYfr27UvZsmU9qvPrr78yc+ZMnn32WZspL6fsMpTJqbi4OFavXs2WLVsAaNy4Ma1bt6Zy5co+jswUVocOHWLq1KlERETQs2dPj+qoKq1atWL//v3s2rWLMmXK5EksdhnKGA+FhobSo0cP2rdvzy+//MLGjRvZvHkz9evXp02bNtSoUcPXIZpCJCkpiXnz5lGuXDm6du3qcb3Zs2ezdu1aPv744zxLFNlhZxbGZJL+gN+6des4f/68PeBn8tSCBQvYvHkzQ4cOJSwszKM6qkrLli25cOEC69evz9Pbv+1uKGNyKSkpifXr17N27VoSExPtAT+Tazt27GDOnDm0bdvW44fv0iUmJvLXX39x7bXX5mlMliyMySMpKSls2rSJ1atXZzzg17p1axo1akRAQICvwzMFxOnTp5k8eTIVKlTg/vvv9/hn59SpU4SEhHjtuSBPk4U9ymqMG8WLF6d58+Y89thj9O7dGxFh/vz5/Pvf/+a3334jOTnZ1yEaP6eqLFiwgOTkZHr16pWtPzL+/ve/07x5c1JSUrwYoXvWwW2Mh4oVK0bjxo1p1KhRxgN+ixYtYuXKldx88820aNGC4OBgX4dp/NBvv/3Gnj176NatW7aG0V+/fj3Tp09n1KhRPp9S2C5DGZNDqprxgN+ePXsICgqyB/zMJY4dO8aHH35InTp1GDBggMf9XapK27ZtiYqKIioqyuNnMbLLbp01xstcH/A7cuTIRQ/4NW3alFatWvnNlJjGN1JSUvjqq68IDAyke/fu2boxYu7cuaxevZoPP/zQa4kiO+zMwpg8ZA/4GVc//PADv/zyC/3796devXrZqturVy/27dvHhg0bvHojhd0NZYwPnTp1ijVr1rBx40aSk5OJjIykU6dOdsttEbJv3z5mzJjBjTfeyJ133pnt+qmpqRw9epTq1at7Ibr/sctQxvhQuXLl6NKlC+3atWP58uWsXbuWc+fO0b17d5tPowg4d+4c8+fPp2LFih7NeucqLi6OYsWKUbFiRa8niuywn1pjvCgkJIQ77riDDh06sHnzZubOnevzWyCN93333XecPn2a3r17ExgYmK26Tz/9NA0bNuTs2bNeii5nLFkY42UiQrt27ejcuTM7duxg9uzZ9mxGIfbHH3+wdetW2rdvn+1xxTZu3Mi0adMYNGgQISEhXoowZ7yaLESki4jsFJHdIvJMFuXuFhEVkebO5XAROScim5yv970ZpzH5ITIyku7du7N3715mzpzJ+fPnfR2SyWMnT57k22+/pVatWrRp0yZbdVWVf/zjH1SqVIkxY8Z4KcKc81qyEJEAYBLQFbgeGCAi11+mXBngceDXTJv2qGoT5+shb8VpTH5q2rQpffr0ITo6mhkzZnDmzBlfh2TySFpaGvPnz0dV6dWrV7b7pr766itWrlzJiy++6JdzqnjzzOImYLeq7lXVC8BsoMdlyr0ITATszyxTJDRs2JD+/fsTGxvLtGnTSEhI8HVIJg+sWbOGAwcO0LVrVypUqJDt+j/99BONGjVi+PDhXogu97yZLGoAh1yWo53rMohIU6CWqn5zmfp1ROR3EflJRNp6MU5j8l3dunUZOHAgCQkJTJ06lRMnTvg6JJMLR44c4ccff6RBgwbccMMNOdrHu+++y+rVq30+rMeVeDNZXO6G8oyHOkSkGPAW8ORlyh0BwlS1KfAEMEtELnmEUURGiMh6EVkfGxubR2Ebkz/Cw8MZMmQISUlJfPzxx9jPcMGUnJzMV199RUhICHfeeWe2n6U5duwYUVFRAH7xpPaVeDNZRAO1XJZrAjEuy2WARsAKEdkPRAILRaS5qiapajyAqm4A9gDXZT6Aqk5R1eaq2tyekDUFUfXq1Rk6dCgAU6dOJSYmJusKxu8sXbqUuLg4evbsmaM7mMaMGUOTJk04fvy4F6LLO95MFuuAuiJSR0QCgf7AwvSNqnpKVUNVNVxVw4G1QHdVXS8ilZ0d5IjI1UBdYK8XYzXGZ6pUqcKwYcMICgpi+vTpHDhwwNchGQ/t3r2b3377jZtvvplrrrkm2/U3bdrERx99xIgRI6hYsaIXIsw7XksWqpoCPAp8D+wA5qjqNhGZICLd3VRvB2wRkc3AXOAhVfXvtGtMLlSsWJFhw4ZRtmxZZs6cye7du30dknHj7NmzLFiwgMqVK3Pbbbdlu76q8sQTT1ChQgXGjRvnhQjzlo0NZYwfOXPmDDNnzuTYsWP06dOH66+/5G5z4wdUlTlz5hAVFcXw4cOpWrVqtvexYMECevbsyXvvvccjjzzihSg9YzPlGVMAlSpViiFDhlCjRg3mzp3Lpk2bfB2SuYxNmzbx559/cuutt+YoUQBERUVxww038OCDD+ZxdN5hycIYPxMcHMygQYOoU6cOCxYs4NdfMz+vmv8mTpzI8uXLL1q3fPlyJk6c6KOIfOf48eMsWrSI8PBwWrZsmeP9PPXUU6xbt85vb5XNzJKFMX4oMDCQAQMGUL9+fRYvXszKlSvx5SXjFi1a0K9fv4yEsXz5cvr160eLFi18FpMvpKWlMW/ePIoVK0bPnj1zNOR8bGxsxv9jiRIl8jpEr7FkYYyfKl68OH379iUiIoLly5ezdOlSnyWMDh06MGfOHPr168e4cePo168fc+bMoUOHDj6Jx1d+/vlnoqOjufPOO3M8JMe4cePo1KkT0dHReRyddxWM8x9jiqj0v2ADAwP55ZdfSEpKolu3bj6ZRKlDhw6MHDmSF198kbFjxxa5RBEdHc1PP/1E48aNadSoUY728ccffzBlyhQeffRRatasmccRepedWRjj50SEO+64g9atW7NhwwbmzZtHampqvsexfPlyJk+ezNixY5k8efIlfRiF2YULF5g3bx5ly5bljjvuyNE+0m+VLVeuHOPHj8/jCL3PziyMKQBEhNtvv53g4GCWLVvGhQsXuPvuu/OtczS9jyL90lOHDh2K1KWo77//nuPHjzNkyBCCg4NztI9vv/2WpUuX8s477/j9A3iXY2cWxhQgbdq04Y477mDnzp3MmjWLCxcu5Mtx161bd1FiSO/DWLduXb4c35f+/PNPNm7cSOvWrQkPD8/xfhISEmjVqhUjR47Mu+DykT2UZ0wBtHnzZhYsWECNGjW49957KVmypK9DKpQSExOZPHkyZcuWZfjw4QQEBORqf6rqk/6mrNhDecYUYjfccAN9+/blyJEjTJ8+3SZR8gJVZcGCBVy4cIHevXvnOFHExcUxbdo00tLS/C5RZIclC2MKqAYNGjBgwADi4+OZOnUqp06d8nVIhcr69evZvXs3HTt2JDejWj///PP87W9/yxiGvKCyZGFMAXbNNddw3333kZiYyNSpU4mPj/d1SIVCbGwsS5Ys4dprr83Vg4fbt2/n/fff56GHHqJevXp5GGH+s2RhTAEXFhbGkCFDSE5OZurUqRw9etTXIRVoqampzJs3jxIlStC9e/dcXTp68sknKV26NC+88EIeRugbliyMKQSqVavG0KFDKVasGNOmTePw4cO+DqnAWrFiBUeOHKF79+6UKVMmx/tZtGgRixcvZvz48YSGhuZhhL5hycKYQqJy5coMGzaMkiVLMmPGDPbv3+/rkAqcAwcOsGrVKpo2bUr9+vVzta/g4GC6devm0+HH85IlC2MKkQoVKjBs2DDKlSvHzJkz2bVrl69DKjDOnz/PvHnzqFChAl26dMn1/jp06MA333xDYGBgHkTne5YsjClkypQpw9ChQ6lSpQqff/45W7du9XVIBcKiRYtISEigd+/eufqCP378OBMmTCAxMTEPo/M9SxbGFEIhISEMGTKEWrVq8eWXX7Jx40Zfh+TXtm7dypYtW2jXrl2uB/ibMGECL7zwAnv37s2j6PyDJQtjCqmgoCAGDhzItddey9dff82aNWt8HZJfSkhI4Ntvv6VGjRq0a9cuV/v6888/mTRpEg888AARERF5FKF/sGRhTCFWokQJ+vfvz/XXX8+SJUtYsWKFTydR8jeqyvz580lNTaV3794UK5a7r8SnnnqKkJAQJkyYkEcR+g8bddaYQi4gIIA+ffoQGBjITz/9RFJSEp06dSrQQ0/klbVr17Jv3z7uuuuuXI8Eu2TJEr799ltee+01qlSpkkcR+g9LFsYUAcWKFaN79+4EBgaydu1akpKSuPPOO3P9l3RBdvToUZYtW0b9+vVp2rRprvdXvXp1Bg8ezGOPPZYH0fkfSxbGFBEiQpcuXQgODmblypVcuHCBXr165Xok1YLowoULfPXVV5QsWZK77rorT86yGjVqxPTp0/MgOv9kycKYIkRE6NChA0FBQfzwww9cuHCBvn37UqJECV+H5lXnz5/n0KFDHDx4kIMHD3L48GFSU1MZOHAgISEhudr3yZMneeaZZxg7diw1atTIo4j9j9tkISLBwJ1AW6A6cA7YCnyrqtu8G54xxhtatWpFUFAQ33zzDZ9++ikDBgwgKCjI12HlmcTERA4ePMiBAwc4ePAgR48ezZhLonr16tx0001cd911uZrMKN2LL77IlClTeOihhwp1sshy8iMReR64C1gBbACOAcHAdUAH5/snVXWLtwN1xyY/Mib7/vjjD+bNm0e1atXy5K9sX1BVTpw4cVFyOH78OADFixenZs2a1K5dm7CwMGrWrJmnT1Tv2rWLhg0bMnToUD788MM8229+8nTyI3fJopuqfpvF9ipAmKr6/FvakoUxObNz506++OILKlWqxKBBg3I1eF5+UFWOHTuWkRgOHjzI6dOnAcd4TGFhYYSFhVG7dm2qVavm1T6ZHj16sHz5cqKiorjqqqu8dhxv8jRZZHkZyk2iKK6qx3CcbRhjCqh69eoxcOBAPvvsM6ZOncrgwYMpX768r8PKkJqaSkxMTMaZw6FDhzh//jzgGNok/ayhdu3aVK5cOd9uCV66dCkLFy7klVdeKbCJIjvcnVmsUtU2zvefqOp9Lts2qmqzfIjRI3ZmYUzuREdH8+mnnxIYGMh9993ns2G1L1y4cFFndHR0NCkpKQBUqlQpIzGEhYVRvnx5nz0v8tdff/Hmm28yYcIEgoODfRJDXsiry1C/q2pT5/uLkoPrNn9gycKY3Pvrr7+YOXMmqsp9991H1apVvX7Ms2fPXtTfcOTIkYzO6KpVq2ZcVgoLC6N06dJej6eoyZPLUEBW4wLYmAHGFDJVq1Zl2LBhzJgxg2nTpjFw4EBq1aqVp8c4deoUBw4cyEgOcXFxgONJ8xo1atCmTRvCwsKoVauWX96hderUKQYPHsxLL71E48aNfR1OvnGXLMqLSC8cY0iVF5HezvUClHO3cxHpArwDBAAfqeorVyh3N/AF0CK9s1xERgN/A1KBx1X1ew/aY4zJpUqVKnH//fczY8YMPvnkE/r378/VV1+do32pKnFxcRd1Rp86dQpwDHRYq1YtIiIiqF27NtWrV6d4cf9/9Ovll1/m66+/Zty4cb4OJV+5uww1NavKqjosi7oBwC6gIxANrAMGqOr2TOXKAN8CgcCjqrpeRK4HPgNuwvFsx1LgOlVNvdLx7DKUMXkrMTGRTz75hPj4eO6++26PZo5LS0vjyJEjGYnhwIEDnDt3DoBSpUpd1BldpUqVAjfcyO7du7n++usZOHAgU6dm+fVYYOTV3VBXTAYeuAnYrap7nQHNBnoA2zOVexGYCDzlsq4HMFtVk4B9IrLbuT8bY9mYfFK6dGmGDh3Kp59+ypw5c+jZs+clw24nJycTHR2dkRwOHTpEcnIy4Ji1r169ehn9DRUrVizwgxc+/fTTBAYG8vLLL/s6lHyXZbIQkbuALap6wLk8DugDHAD+rqr7sqheAzjkshwN3Jxp/02BWqr6jYg8lanu2kx1C++jkcb4qZIlS3Lfffcxe/Zs5s2bx7lz56hQoULGZaWYmBjS0tIAuOqqq2jSpEnGmYO/P6+RXStXrmTevHm89NJLVK9e3dfh5Dt3FwhfBiIBROROYBAwAGgKvA90zqLu5f6EyLjmJSLFgLeAodmt67KPEcAIgLCwsCxCMcbkVFBQEPfeey9z585l8eLFgGMU2+rVqxMZGUnt2rWpVasWJUuW9HGk3tWiRQveeOMNRo4c6etQfMLt3VCqetb5vjfwX1XdAGwQkYfd1I0GXG+jqAnEuCyXARoBK5ynplWBhSLS3YO66cFNAaaAo8/CTTzGmBwqUaIE/fr1Y/v27ZQpU4YaNWoU+sEHMytZsiRPPPGEr8PwGXe9SyIipZ1nAbcBy1y2uXsKZR1QV0TqiEgg0B9YmL5RVU+paqiqhqtqOI7LTt2dd0MtBPqLSJCI1AHqAr9lq2XGmDwVEBBA48aNCQ8PL1KJIiEhgZYtW/Ljjz/6OhSfcndm8TawCUgAdrjc1toUOJJVRVVNEZFHge9x3Dr7sapuE5EJwHpVXZhF3W0iMgdHZ3gK8EhWd0IZY4y3/Otf/2Lt2rWFrg8mu7K8dRZARGoAVYDNqprmXFcNKKGqB70fomfs1lljTF7bt28fDRo0oF+/fsyYMcPX4XhFntw6KyKuYz81ucxtb36TLIwxJq+NGjWKgIAA/t//+3++DsXn3F2GWg9sA2Kdy67ZQoFbvRGUMcb42vr16/niiy944YUXqFmzpq/D8Tl3yeJJHM9VnANmA/NUNdHrURljjI/deOONfP7559x5552+DsUvZHk3lKq+5Ryi/FEct7IuE5E5ItIkX6IzxhgfSE1NRUTo169fgZw90Bs8GpjF+aT2AmAJjmE3rvNmUMYY4yuJiYk0atSI2bNn+zoUv+Kug/tqHM+Y9gQ6AAAgAElEQVRH9MAxdMds4GVVPZ8PsRljTL579dVX+fPPPwkPD/d1KH7FXZ/FbmALjrOKBCAMeDj9rihVfdOr0RljTD46cOAAr7/+Ovfeey+RkZG+DsevuEsWE/jfmEw2RZUxplB75plnEBFeeeWyU+8Uae6GKH8+n+Iwxhif2rVrF7Nnz2bcuHF5PjtgYeCuz2IMMElVT1xh+61AiKp+443gjDEmv1x33XX88ssvl8zZYRzcXYb6A/hGRM4DG3E8nBeMY2C/JjhmsLNHG40xBdrZs2cJCQmhZcuWvg7Fb7l7zmKBqrYGHsLxJHcAjo7umcBNqvoPVY3Nah/GGOPPzpw5Q4MGDXjzTbtfJysezY6uqlFAlJdjMcaYfPfaa69x8OBBbr75ZveFi7CCNVu6McbkoUOHDjFx4kTuueceWrdu7etw/JolC2NMkTV69GjS0tJ49dVXfR2K33ObLEQkQET+kR/BGGNMfomJieHLL7/kySefpHbt2r4Ox++57bNQ1VQR6QG8lQ/xGGNMvqhevTrbt28nNDTU16EUCB51cAOrReQ94HPgTPpKVd3olaiMMcaLoqKiuPbaa6lTp46vQykwPO2zaAU0xDH8xxvO1+veCsoYY7xlx44dNGnShH/961++DqVA8fTW2Q7eDsQYY7zt7Nmz9OvXj1KlSjF06FBfh1OgeJQsRKQcMB5o51z1EzBBVU95KzBjjMlrjz32GNu2bWPx4sVUr17d1+EUKJ5ehvoYOA30c74SgKneCsoYY/LajBkz+Pjjj3nuuefo1KmTr8MpcDzt4L5GVfu4LL8gIpu8EZAxxnhD5cqV6d27N+PHj/d1KAWSp2cW50SkTfqCiLQGznknJGOMyTuqjil5unbtypdffknx4p7+jWxcefq/9hAww9l3AXACGOKdkIwxJu+MGDGCunXr8vTTT/s6lALNkye4iwH1VPUGIAKIUNWmqrrF69EZY0wuTJs2jY8++ojExERfh1LguU0WqpoGPOp8n6CqCV6Pyhhjcmnbtm08/PDDdOjQwfop8oCnfRY/iMhTIlJLRCqmv7wamTHG5FBiYiJ9+/albNmyzJo1i4CAAF+HVOB52mdxv/PfR1zWKXB13oZjjDG599NPP7Fnzx4WLVpE1apVfR1OoeA2WTj7LAap6up8iMcYY3KtW7du7Nmzh5o1a/o6lELD0z4LGwfKGOP3/vjjDxYtWgRgiSKPedpnsURE+oiIZGfnItJFRHaKyG4ReeYy2x8SkT9EZJOIrBKR653rw0XknHP9JhF5PzvHNcYUPadPn6Zv374MHz6cs2fP+jqcQsfTPosngFJAqoicAwRQVS17pQoiEgBMAjoC0cA6EVmoqttdis1S1fed5bsDbwJdnNv2qGqTbLXGGFMkqSoPPfQQUVFRLFu2jJCQEF+HVOh4OupsmRzs+yZgt6ruBRCR2UAPICNZZLoNtxSOTnNjjMmWjz76iFmzZvHiiy/Svn17X4dTKHl0GUocBonIWOdyLRG5yU21GsAhl+Vo57rM+35ERPYAE4HHXTbVEZHfReQnEWnrSZzGmKJn3759PPbYY3Ts2JHRo0f7OpxCy9M+i/8ALYF7ncuJOC4xZeVy/RuXnDmo6iRVvQYYBYxxrj4ChKlqUxyXwGaJyCWXvERkhIisF5H1sbGxnrXEGFOohIeH89577zFz5kx7nsKLPE0WN6vqI8B5AFU9AQS6qRMN1HJZrgnEZFF+NtDTuf8kVY13vt8A7AGuy1xBVaeoanNVbV65cmUPm2KMKQxUlZiYGESE4cOHU6VKFV+HVKh5miySnR3WCiAilYE0N3XWAXVFpI6IBAL9gYWuBUSkrstiNyAqff/O4yEiVwN1gb0exmqMKQKmTJlCvXr12LZtm69DKRI8vRvqXWAeUEVEXgbu5n+XjC5LVVNE5FHgeyAA+FhVt4nIBGC9qi4EHhWR24FkLh7Jth0wQURSgFTgIVU9ns22GWMKqU2bNvH3v/+dW265hQYNGvg6nCJB0sd6d1tQpD5wG46+iGWqusObgWVX8+bNdf369b4OwxjjZQkJCdx4442cPXuWTZs2YZegc0dENqhqc3flPJ4FRFX/BP7MVVTGGJMLqsqIESPYt28fy5cvt0SRjzztszDGGJ9LTU2lXLlyvPTSS7Rta3fU5yebX9AYU2AUL16cDz74AE8vn5u8Y2cWxhi/d+rUKbp27crmzZsByOYwdSYPWLIwxvg1VWX48OH88MMPNj2qD9llKGOMX5s0aRJz587l1VdfpXXr1r4Op8iyMwtjjN9av349Tz75JN26deOpp57ydThFmiULY4zfeuutt7jqqquYPn06xYrZ15Uv2WUoY4zfmjZtGgcPHqRSpUq+DqXIs1RtjPE7X3/9NXFxcZQoUYJrrrnG1+EYLFkYY/zMunXr6NOnD888c8lMzMaHLFkYY/zGiRMn6NevH9WqVWPixIm+Dse4sD4LY4xfUFXuv/9+oqOj+fnnn6lYsaKvQzIuLFkYY/zCBx98wPz583njjTeIjIz0dTgmE0sWxhi/0KdPH+Lj4/nHP/7h61DMZVifhTHGp06fPk1ycjKVK1fmueees3Gf/JQlC2OMz6gqgwYNonPnzqSluZup2fiSJQtjjM+89dZbLFy4kB49etgT2n7OPh1jjE+sXbuWUaNG0atXLx5//HFfh2PcsGRhjMl3x48f55577qFWrVp8/PHH1k9RAFiyMMbku2PHjlGmTBnmzJlD+fLlfR2O8YDdOmuMyXf169dny5Yt1k9RgNgnZYzJN7/88guPPPII58+ft0RRwNinZYzJF/Hx8dxzzz0sXryYpKQkX4djsskuQxljvC4tLY3Bgwdz7Ngx1qxZQ7ly5XwdkskmSxbGGK977bXX+O6775g0aRLNmjXzdTgmBwp1skhOTiY6Oprz58/7OhSTA8HBwdSsWZMSJUr4OhSTCydPnuSVV16hb9++jBw50tfhmBwq1MkiOjqaMmXKEB4ebvdxFzCqSnx8PNHR0dSpU8fX4ZhcKF++PL/++itVq1a138MCrFB3cJ8/f55KlSrZD2gBJCJUqlTJzgoLsLS0NL7++mtUleuuu46yZcv6OiSTC4U6WQCWKAow++wKtldffZXu3buzZMkSX4di8kChTxbGmPz3888/M2bMGO655x46derk63BMHvBqshCRLiKyU0R2i8gls6+LyEMi8oeIbBKRVSJyvcu20c56O0Wkszfj9JaTJ0/yn//8J0d177jjDk6ePJllmXHjxrF06dIc7d8Ybzl27Bj9+/fnmmuuYcqUKXaGWEh4LVmISAAwCegKXA8McE0GTrNUtbGqNgEmAm86614P9AcaAl2A/zj3V6BklSxSU1OzrPvdd9+5HTNnwoQJ3H777TmOL7syx5ySkuJRPU/LmYJPVRkyZAjx8fHMmTPH+ikKEW+eWdwE7FbVvap6AZgN9HAtoKoJLoulAHW+7wHMVtUkVd0H7Hbur0B55pln2LNnD02aNOGf//wnK1asoEOHDtx77700btwYgJ49e3LjjTfSsGFDpkyZklE3PDycuLg49u/fT4MGDXjggQdo2LAhnTp14ty5cwAMHTqUuXPnZpQfP348zZo1o3Hjxvz5558AxMbG0rFjR5o1a8aDDz5I7dq1iYuLuyTWJUuW0LJlS5o1a0bfvn1JTEzM2O+ECRNo06YNX3zxBe3bt+fZZ5/llltu4Z133uHAgQPcdtttREREcNttt3Hw4MGM2J544gk6dOjAqFGjvPefbPyKiPDYY48xZcoUmjRp4utwTB7y5q2zNYBDLsvRwM2ZC4nII8ATQCBwq0vdtZnq1rhM3RHACICwsDC3AbVv3/6Sdf369ePhhx/m7Nmz3HHHHZdsHzp0KEOHDiUuLo677777om0rVqzI8nivvPIKW7duZdOmTRnlf/vtN7Zu3ZpxO+jHH39MxYoVOXfuHC1atKBPnz5UqlTpov1ERUXx2Wef8eGHH9KvXz++/PJLBg0adMnxQkND2bhxI//5z394/fXX+eijj3jhhRe49dZbGT16NIsXL74oIaWLi4vjpZdeYunSpZQqVYpXX32VN998k3HjxgGO5x1WrVoFwPvvv8/Jkyf56aefALjrrrsYPHgwQ4YM4eOPP+bxxx9n/vz5AOzatYulS5cSEFDgTgpNDpw9e5aQkJDL/h6Zgs+bZxaXu1Cpl6xQnaSq1wCjgDHZrDtFVZuravPKlSvnKtj8ctNNN1303MC7777LDTfcQGRkJIcOHSIqKuqSOnXq1Mn4K+3GG29k//79l9137969LymzatUq+vfvD0CXLl2oUKHCJfXWrl3L9u3bad26NU2aNGH69OkcOHAgY/s999xzUXnX5TVr1nDvvfcCcN9992UkFYC+fftaoigijh49Sv369fnvf//r61CMl3jzzCIaqOWyXBOIyaL8bGByDut6JKszgZCQkCy3h4aGuj2T8ESpUqUuimfp0qWsWbOGkJAQ2rdvf9nnCoKCgjLeBwQEZFyGulK5gICAjH4C1Uty7CVUlY4dO/LZZ5+5jflyy65cOzOzKmcKj9TUVAYNGkRsbCwtWrTwdTjGS7x5ZrEOqCsidUQkEEeH9ULXAiJS12WxG5D+Z/VCoL+IBIlIHaAu8JsXY/WKMmXKcPr06StuP3XqFBUqVCAkJIQ///yTtWvXXrFsTrVp04Y5c+YAjn6JEydOXFImMjKS1atXs3v3bsBxOWHXrl0e7b9Vq1bMnj0bgE8//ZQ2bdrkUeSmoHj55ZdZunQp7733HhEREb4Ox3iJ15KFqqYAjwLfAzuAOaq6TUQmiEh3Z7FHRWSbiGzC0W8xxFl3GzAH2A4sBh5R1axvH/JDlSpVonXr1jRq1Ih//vOfl2zv0qULKSkpREREMHbsWCIjI/M8hvHjx7NkyRKaNWvGokWLqFatGmXKlLmoTOXKlZk2bRoDBgwgIiKCyMjIjA5yd959912mTp1KREQEn3zyCe+8806et8H4rx9//JHnn3+eQYMGcf/99/s6HONF4sllioKgefPmun79+ovW7dixgwYNGvgoIv+QlJREQEAAxYsXZ82aNYwcOTKjw70gsM/Qv02aNIn333+fNWvWULp0aV+HY3JARDaoanN35Qr1QIIGDh48SL9+/UhLSyMwMJAPP/zQ1yGZQuSRRx7hgQceIDAw0NehGC+zZFHI1a1bl99//93XYZhC5s0336RBgwZ07drVEkURYWNDGWOyZdmyZTz11FMZNzaYosGShTHGY3/99RcDBw6kfv36OR73zBRMdhnKGOORuLg4BgwYwOnTp1m2bJk9R1PEWLIwxmRJVRERlixZwooVK5gxYwYNGzb0dVgmn9llKD+TfvthTEzMJWNRpWvfvj2ZbxPO7O233+bs2bMZy54MeW5MOlXl+++/p0uXLrz22msA3H333WzdupX77rvPx9EZX7Bk4aeqV6+eMaJsTmROFp4MeZ5XMg9J7ukQ5e6GbTfed+7cOaZMmULDhg3p0qULmzdvzniIMzAw0M4oijBLFl40atSoizoBn3/+ed544w0SExO57bbbMoYTX7BgwSV19+/fT6NGjQDHL3D//v2JiIjgnnvuuWhsqJEjR9K8eXMaNmzI+PHjAcdT1TExMXTo0IEOHToA/xvyHBy3PTZq1IhGjRrx9ttvZxzvSkOhu4qNjaVPnz60aNGCFi1asHr16oy2jRgxgk6dOjF48GCmTZtG3759ueuuu+jUqROqyj//+U8aNWpE48aN+fzzzwEuO2y78Z0hQ4bw4IMPEhwczIwZMzhw4AAjR470dVjGDxSZPovFixfz119/5ek+q1atSpcuXa64vX///vzf//0fDz/8MABz5sxh8eLFBAcHM2/ePMqWLUtcXByRkZF07979ijOKTZ48mZCQELZs2cKWLVto1qxZxraXX36ZihUrkpqaym233caWLVt4/PHHefPNN1m+fDmhoaEX7WvDhg1MnTqVX3/9FVXl5ptv5pZbbqFChQoeDYX+97//nX/84x+0adOGgwcP0rlzZ3bs2JGx71WrVlGyZEmmTZvGmjVr2LJlCxUrVuTLL79k06ZNbN68mbi4OFq0aEG7du0ALhm23eSfDRs28Pbbb/Pyyy8TFhbG008/zaOPPkrbtm1thjtzkSKTLHyhadOmHDt2jJiYGGJjY6lQoQJhYWEkJyfz7LPPsnLlSooVK8bhw4c5evQoVatWvex+Vq5cyeOPPw5ARETERYO1zZkzhylTppCSksKRI0fYvn17loO5rVq1il69emXcydK7d29+/vlnunfv7tFQ6EuXLmX79u0ZywkJCRmDJXbv3p2SJUtmbOvYsSMVK1bMOO6AAQMICAjgqquu4pZbbmHdunWULVv2kmHbjXelpqayYMEC3n77bX7++WdKly5N3759CQsLo3lzt6M+mCKqyCSLrM4AvOnuu+9m7ty5/PXXXxnzSnz66afExsayYcMGSpQoQXh4+GWHJnd1ub/y9u3bx+uvv866deuoUKECQ4cOdbufrMYC82Qo9LS0NNasWXNRUkiX1VDmWR3XbsHMPxcuXCAiIoKdO3dSu3Zt3njjDf72t79Rrlw5X4dm/Jz1WXhZ//79mT17NnPnzs24u+nUqVNUqVKFEiVKsHz58osmGrqcdu3a8emnnwKwdetWtmzZAjj+qi9VqhTlypXj6NGjLFq0KKPOlYZHb9euHfPnz+fs2bOcOXOGefPm0bZtW4/b06lTJ957772MZU8HJWzXrh2ff/45qampxMbGsnLlSm66qcDNlFsg7du3j0mTJgGOTurBgwczd+5cdu/ezRNPPGGJwnikyJxZ+ErDhg05ffo0NWrUoFq1agAMHDiQu+66i+bNm9OkSRPq16+f5T5GjhzJsGHDiIiIoEmTJhlfsjfccANNmzalYcOGXH311bRu3TqjzogRI+jatSvVqlVj+fLlGeubNWvG0KFDM/YxfPhwmjZtesXZ9zJ79913eeSRR4iIiCAlJYV27drx/vvvu63Xq1cv1qxZww033ICIMHHiRKpWrerxUOgme1SVVatW8fbbbzN//nyKFStGjx49qFmzJs8++6yvwzMFkA1RbvyafYbZt337dgYPHsyGDRuoWLEiDz74IA8//DA1a9b0dWjGD9kQ5cYUIfHx8Rw+fJiIiAiqV6+OiDB58mQGDx5MSEiIr8MzhYAlC2MKsB07dvD2228zY8YM6tevz8aNGylfvjzr1q3zdWimkLEObmMKoNWrV9OlSxeuv/56pk+fzqBBg/jkk0/s2QjjNXZmYUwBce7cOVSVkJAQoqKi2Lx5My+++CIPPvgglStX9nV4ppCzMwtj/NyRI0cYM2YMtWrVyrjzbODAgezfv58xY8ZYojD5ws4sjPFTGzdu5K233uLzzz8nJSWFHj160KpVKwBKlCjh4+hMUWNnFl508uTJXM0mlnnkWFP4ud7K/s9//pP58+czcuRIoqKimDdvHpGRkT6MzhRlliycJk6ceNHDawDLly9n4sSJOd6nr5NFTocK97ScyTsJCQm88847XH/99cTExADwwQcfEB0dzTvvvMM111zj4whNUWeXoZxatGhBv379mDNnDh06dGD58uUZyzn1zDPPsGfPHpo0aULHjh157bXXeO2115gzZw5JSUn06tWLF154gTNnztCvXz+io6NJTU1l7NixHD16NGOY8dDQ0EsS2YYNG3jiiSdITEwkNDSUadOmUa1aNdq3b0+rVq1YvXo13bt3548//qBixYr8/vvvNGvWjOeee47777+fvXv3EhISwpQpU4iIiOD5558nJiaG/fv3ExoayqxZs3L7X2o8sG/fPv7973/z3//+l4SEBFq3bk18fDzVq1fn2muv9XV4xvyPqhaK14033qiZbd++/ZJ1Wfnxxx81NDRUx44dq6Ghofrjjz9mq35m+/bt04YNG2Ysf//99/rAAw9oWlqapqamardu3fSnn37SuXPn6vDhwzPKnTx5UlVVa9eurbGxsZfs98KFC9qyZUs9duyYqqrOnj1bhw0bpqqqt9xyi44cOTKj7JAhQ7Rbt26akpKiqqqPPvqoPv/886qqumzZMr3hhhtUVXX8+PHarFkzPXv2bK7anNey+xkWJHFxcRoYGKjFixfXAQMG6K+//urrkEwRBKxXD75j7czCRYcOHRg5ciQvvvgiY8eOzZg4KK8sWbKEJUuW0LRpUwASExOJioqibdu2PPXUU4waNYo777zT7cB+O3fuZOvWrXTs2BFwDDmdPu4UwD333HNR+b59+xIQEAA4hgr/8ssvAbj11luJj4/n1KlTwKVDjJu8lZyczBdffMHGjRt5/fXXqVSpEv/9739p3769DcVh/J4lCxfLly9n8uTJjB07lsmTJ18001xeUFVGjx7Ngw8+eMm2DRs28N133zF69Gg6derEuHHjstxPw4YNWbNmzWW3Z3eo8PQHuWyo8LyR/n8sIpw/f54jR44we/Zs3nvvPWJiYqhfvz4vvPACpUqVumRyKWP8lSULJ9c+ivQk4bqcE5mHCe/cuTNjx45l4MCBlC5dmsOHD1OiRAlSUlKoWLEigwYNonTp0kybNu2i+plnu6tXrx6xsbGsWbOGli1bkpyczK5duzyaHzl9uPOxY8eyYsUKQkNDKVu2bI7a529SUlJISkoiKSmJwMBASpcuTXJyMtu2beP8+fMZ286fP0+DBg2oW7cuJ06cYNasWRdtS0pKomfPnkRGRrJ3717GjBlzSf3x48fTsWNH1qxZw4ABAzK2p//79ddf061bNxYvXkyvXr0Ax2RQH330EZ07d6ZYMbu3xBQsliyAM2fO8P333/POO+8QHh7Ovn37CA8PZ+bMmaxbt47mzZtnzF/tqkaNGgQGBnLq1CmOHz9+yfZatWrRunVrGjRoQNu2bRk9ejSdO3fmxhtvBKBChQp8+umn/P7774wZM4ZixYpRvHhxXnzxRfbt28cDDzxA165dCQ0N5ZNPPrlo3++++y6jRo3i1KlTJCUlMWTIEEJCQjh//jzR0dFUrlyZ2rVrA467svbu3QvA0KFDGTVqFA0aNKBcuXJMnz6dI0eOcOLECZKSkjLKBQYGZlwaOXz4MElJSRcdPzg4mOrVqwNw6NAhkpOTL9pesmTJjEtjBw8evOQOq1KlSnHVVVcBjvm/U1NTL9pepkwZqlSpAjjm0Mj8ZXzffffx3HPPkZiYyFVXXUVSUtJF+xg7diwTJkwgPj4+47Kfq1dffZWnn36a+Ph4Hn300Yu2BQYGcs011xAZGcm5c+dYt24dwcHBBAUFERQURHBwcMaXfYUKFWjfvv1F24KCgjLuXmratCkffPABrVu39iiZG+OvbIhyHJMRHTx48JL111xzDSEhIRw/fpzDhw9fsr1u3boEBwcTFxfHkSNHLtler149AgMDOXbsGEePHr1ke4MGDShevDh//fUXsbGxl2xv2LBhxrSrmZORiNCoUSPA8WV98uTJi7YHBARw/fXXA3DgwAESEhIu2h4YGEi9evUAxx05iYmJF20PDg6mbt26AOzZs+eSW3hDQkIyvhCjoqIumaGvTJkyhIeHA44+lszJpGzZsoSFhQGOzylzsihfvjw1a9Zkx44dPPTQQxQrVuyiL+Nu3boxePBgUlJSGD16dMaXdXqZm266icjISJKSkvjuu+8y6qX/W6tWLapUqUJKSgonTpzIqBsYGGjjK5kixdMhyi1ZGL9mn6Ex3uVpsvDqhVMR6SIiO0Vkt4g8c5ntT4jIdhHZIiLLRKS2y7ZUEdnkfC30ZpzGGGOy5rU+CxEJACYBHYFoYJ2ILFTV7S7Ffgeaq+pZERkJTATS7/s8p6pNchuHqtplhQKqsJz1GlMYePPM4iZgt6ruVdULwGygh2sBVV2uqukXw9cCeXqzeXBwMPHx8falUwCpKvHx8QQHB/s6FGMM3r0bqgZwyGU5Grg5i/J/Axa5LAeLyHogBXhFVednriAiI4ARQEZnqauaNWsSHR192c5j4/+Cg4PtYTVj/IQ3k8Xlrv1c9k98ERkENAducVkdpqoxInI18KOI/KGqey7ameoUYAo4Orgz77dEiRLUqVMnp/EbY4xx8uZlqGiglstyTSAmcyERuR14Duiuqhk386tqjPPfvcAK4NKb5Y0xxuQLbyaLdUBdEakjIoFAf+Ciu5pEpCnwAY5EccxlfQURCXK+DwVaA64d48YYY/KR1y5DqWqKiDwKfA8EAB+r6jYRmYBjlMOFwGtAaeAL5x1LB1W1O9AA+EBE0nAktFcy3UVljDEmHxWah/JE5BQQ5bKqHHDqCsvp713XhQKXjunhmczHyk6Zy63PKnbX5cu1KTftyCpOT8pkty3u3vvqM7nStoLYltz8fLm+L4i/K978TLKK05My/tSWuqpazm0pT8YxLwgvYIqny+nvM63zaEx3T46dnTKXW+9pW67Qphy3I7/b4u69rz6TwtSW3Px8ZfGzViDa4s3PpDC1xZN2qGqhmlb162wsf32FMnl17OyUudx6T9tyuTblVn62xZP3OZWbdlxpW0FsS25+vlzf28+XZ/F4Wsaf2uLRPgrNZajcEpH16sH4KP6usLQDrC3+qrC0pbC0A/KnLYXpzCK3pvg6gDxSWNoB1hZ/VVjaUljaAfnQFjuzMMYY45adWRhjjHHLkoUxxhi3LFkYY4xxy5KFGyLSU0Q+FJEFItLJ1/HkhohcLSL/FZG5vo4lJ0SklIhMd34eA30dT24U9M8iXSH7/WggIu+LyFzn/DoFmvP3ZYOI3JkX+yvUyUJEPhaRYyKyNdP6LGfwc6Wq81X1AWAo/5uYKd/lUVv2qurfvBtp9mSzXb2Buc7Po3u+B+tGdtrij59Fumy2wy9+P64km23ZoaoPAf1wjILtV3LwHTAKmJNnAeTmCUZ/fwHtgGbAVpd1AcAe4GogENgMXA80Bjt5PHYAAALBSURBVL7J9KriUu8NoFkhactcX382OWzXaKCJs8wsX8eem7b442eRy3b49Pcjr9qC44+QX4B7fR17btoC3I5j8NahwJ15cXxvzmfhc6q6UkTCM63OmMEPQERmAz1U9V/AJadr4hjh8BVgkapu9G7EV5YXbfFH2WkXjmHvawKb8MOz4my2xW8HxsxOO0RkB37w+3El2f1M1DHA6UIR+RaYlZ+xupPNtpQGSuFIHOdE5DtVTcvN8f3uFy4fXG4GvxpZlH8MR5a+W0Qe8mZgOZCttohIJRF5H2gqIqO9HVwuXKldXwF9RGQyeTdkg7ddti0F6LNId6XPxJ9/P67kSp9JexF5V0Q+AL7zTWjZdtm2qOpzqvp/OBLeh7lNFODdmfL8lccz+AGo6rvAu94LJ1ey25Z4oCD8Ql+2Xap6BhiW38Hk0pXaUlA+i3RXaoc//35cyZXasgLHRGsFSZbfAao6La8OVBTPLDyawa+AKExtcVWY2lVY2lJY2gHWlhwpisnC7Qx+BUhhaourwtSuwtKWwtIOsLbkjK97+L1898BnwBEgGUcG/ptz/R3ALhx3ETzn6ziLWlsKa7sKS1sKSzusLXn7soEEjTHGuFUUL0MZY4zJJksWxhhj3LJkYYwxxi1LFsYYY9yyZGGMMcYtSxbGGGPcsmRhjDHGLUsWxhhj3LJkYYwXiUhjETlQGGZeM0WbJQtjvEhV/8AxXs9gX8diTG5YsjDG+44BDX0dhDG5YcnCGO97BQgSkdq+DsSYnLJkYYwXiUgXHNNbfoudXZgCzJKFMV4iIsHAROBh4A+gkW8jMibnLFkY4z1jgP/f3h0TAQgEQRDcV4MS/KEBZTg5EvKNCL6qOzwDUxvdPTNPxILNiQX8YK11JDmTXN9JLNia50cAVJYFAJVYAFCJBQCVWABQiQUAlVgAUIkFAJVYAFC943iTANx9CncAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xticks([])\n",
    "plt.xlabel(r\"$\\lambda$\")\n",
    "plt.ylabel(r\"error (MSE)\")\n",
    "plt.title(r\"Validation Pipeline\")\n",
    "plt.xscale('log')\n",
    "plt.plot(ls, errorTrain, color='k', linestyle='--', label='training error')\n",
    "plt.plot(ls, errorValid, color='grey',zorder=4,label=\"validation error\")\n",
    "plt.plot([bestLamb], [MSEtest], linestyle='', marker='x', color='k', label=\"test error\")\n",
    "plt.legend(loc='lower left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Precision, recall, and ROC curves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same data as pipeline above, slightly bigger dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = list(parseData(dataDir + \"beer_50000.json\"))[:10000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple bag-of-words model (as in pipeline above, and in Chapter 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "wordCount = defaultdict(int)\n",
    "punctuation = set(string.punctuation)\n",
    "for d in data:\n",
    "    r = ''.join([c for c in d['review/text'].lower() if not c in punctuation])\n",
    "    for w in r.split():\n",
    "        wordCount[w] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts = [(wordCount[w], w) for w in wordCount]\n",
    "counts.sort()\n",
    "counts.reverse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "words = [x[1] for x in counts[:1000]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "wordId = dict(zip(words, range(len(words))))\n",
    "wordSet = set(words)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feature(datum):\n",
    "    feat = [0]*len(words)\n",
    "    r = ''.join([c for c in datum['review/text'].lower() if not c in punctuation])\n",
    "    ws = r.split()\n",
    "    for w in ws:\n",
    "        if w in words:\n",
    "            feat[wordId[w]] += 1\n",
    "    feat.append(1) # offset\n",
    "    return feat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predict whether the ABV is above 6.7 (roughly, above average) from the review text"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.shuffle(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [feature(d) for d in data]\n",
    "y = [d['beer/ABV'] > 6.7 for d in data]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train on first 9000 reviews, test on last 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = sklearn.linear_model.LogisticRegression(max_iter=5000)\n",
    "mod.fit(X[:9000],y[:9000])\n",
    "predictions = mod.predict(X[9000:]) # Binary vector of predictions\n",
    "correct = predictions == y[9000:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.832"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(correct) / len(correct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute precision and recall, we want the output probabilities (or scores) rather than the predicted labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "probs = mod.predict_proba(X[9000:]) # could also use mod.decision_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build a simple data structure that contains the score, the predicted label, and the actual label (on the test set)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "probY = list(zip([p[1] for p in probs], [p[1] > 0.5 for p in probs], y[9000:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.9999998762354996, True, True),\n",
       " (0.012720481291453863, False, False),\n",
       " (0.997069530448619, True, True),\n",
       " (0.00011239818179373418, False, False),\n",
       " (0.9506213023816044, True, True),\n",
       " (0.8511369245118185, True, False),\n",
       " (0.9508087266170231, True, True),\n",
       " (0.0019793193061717392, False, True),\n",
       " (0.903593342768345, True, True),\n",
       " (0.08480983386485495, False, True)]"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probY[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sort this so that the most confident predictions come first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "probY.sort(reverse=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1.0, True, True),\n",
       " (0.9999999999999885, True, True),\n",
       " (0.9999999999999716, True, True),\n",
       " (0.9999999999999378, True, True),\n",
       " (0.9999999999998981, True, True),\n",
       " (0.9999999999980203, True, True),\n",
       " (0.999999999996775, True, True),\n",
       " (0.9999999999825748, True, True),\n",
       " (0.9999999997872742, True, True),\n",
       " (0.9999999997680737, True, True)]"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probY[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Receiver operator characteristic (ROC) curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "xROC = []\n",
    "yROC = []\n",
    "\n",
    "for thresh in numpy.arange(0,1.01,0.01):\n",
    "    preds = [x[0] > thresh for x in probY]\n",
    "    labs = [x[2] for x in probY]\n",
    "    if len(labs) == 0:\n",
    "        continue\n",
    "    \n",
    "    TP = sum([(a and b) for (a,b) in zip(preds,labs)])\n",
    "    FP = sum([(a and not b) for (a,b) in zip(preds,labs)])\n",
    "    TN = sum([(not a and not b) for (a,b) in zip(preds,labs)])\n",
    "    FN = sum([(not a and b) for (a,b) in zip(preds,labs)])\n",
    "       \n",
    "    TPR = TP / (TP + FN) # True positive rate\n",
    "    FPR = FP / (TN + FP) # False positive rate\n",
    "    xROC.append(FPR)\n",
    "    yROC.append(TPR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW5+PHPk7CHsER2kE1UXACBECCoiIQktt661FasLe2trXbR/tRqtRd729L6ur1etdXWXqQtYr1FrK1Vaq1gRatmJSwhARITCEtYQoCwZSPL8/tjJukYskxIzpxZnvfrNS/mnPnOnOeE5Dzz/Z5znq+oKsYYYwxAlNsBGGOMCR6WFIwxxjSzpGCMMaaZJQVjjDHNLCkYY4xpZknBGGNMM0sKxhhjmllSMGFFRPaISLWInBGRwyKySkT6t2iTKCIbROS0iJwUkb+KyOUt2gwQkV+IyD7vZxV7l4e0sV0Rke+ISL6IVIpIqYi8KiJTnNxfY7qbJQUTjv5NVfsDVwHTge83vSAic4H1wBvAKGACkAukichEb5tewLvAFUAqMABIBI4BCW1s8xng/wHfAeKAS4DXgU93NngR6dHZ9xjTXcTuaDbhRET2AF9T1X94l58ArlDVT3uXPwTyVPVbLd73d6BcVZeIyNeAx4GLVPWMH9u8GCgA5qpqdhtt3gf+T1V/613+ijfOq73LCtwL3A/0ANYBZ1T1IZ/PeAP4p6o+LSKjgF8C1wJngJ+r6rN+/IiMaZf1FEzYEpExwA1AsXe5H55v/K+20vyPwCLv8yTgbX8SgtdCoLSthNAJNwOzgcuB1cDtIiIAIjIYSAbWiEgU8Fc8PZzR3u3fLyIpXdy+MZYUTFh6XUROA/uBI8APvevj8PzOH2rlPYeApvMFF7TRpi2dbd+W/1LV46paDXwIKHCN97XbgAxVPQjMAoaq6jJVPauqu4HfAIu7IQYT4SwpmHB0s6rGAtcBk/nXwb4CaARGtvKekcBR7/NjbbRpS2fbt2V/0xP1jOuuAe7wrvoC8Afv83HAKBE50fQA/gMY3g0xmAhnScGELVX9J7AKeNK7XAlkAJ9rpfnn8ZxcBvgHkCIiMX5u6l1gjIjEt9OmEujnszyitZBbLL8M3CYi4/AMK/3Zu34/UKKqg3wesar6KT/jNaZNlhRMuPsFsEhErvIuPwp82Xv5aKyIDBaRnwJzgR9727yE58D7ZxGZLCJRInKBiPyHiJxz4FXVIuDXwMsicp2I9BKRPiKyWEQe9TbbCtwqIv1EZBJwV0eBq+oWoBz4LbBOVU94X8oGTonIIyLSV0SiReRKEZl1Pj8gY3xZUjBhTVXLgd8DP/AufwSkALfiOQ+wF89lq1d7D+6oai2ek80FwDvAKTwH4iFAVhub+g7wK+A54ASwC7gFzwlhgJ8DZ4Ey4EX+NRTUkZe9saz22acG4N/wXHJbgmfY67fAQD8/05g22SWpxhhjmllPwRhjTDNLCsYYY5pZUjDGGNPMkoIxxphmIVd4a8iQITp+/Hi3wzDGmJCyadOmo6o6tKN2IZcUxo8fT05OjtthGGNMSBGRvf60s+EjY4wxzSwpGGOMaWZJwRhjTDNLCsYYY5pZUjDGGNPMsaQgIitF5IiI5LfxuojIs94J0beJyAynYjHGGOMfJ3sKq/BMet6WG4CLvY+7gf91MBZjjDF+cOw+BVX9QETGt9PkJuD33hmmMkVkkIiMVNXumNbQGGOCWkNDA1VVVVRVVVFZWdnhv5WVldx4443MmuXstBlu3rw2Gp/pB4FS77pzkoKI3I2nN8HYsWMDEpwxJrLV1dWdc1D29wDe2gG95bra2lq/Y7n00ksZNWoUI0aMCOukIK2sa3VyB1VdAawAiI+PtwkgjIlwqkpNTU23HJzb+re+vr5TMUVFRRETE0O/fv0+8W9MTAyDBg1q9bXW2vqui4qKorCwkIkTJ3LFFVcQFeX8tUFuJoVS4EKf5THAQZdiMca0Q1WprKykoqKC48ePU1FR0fyora2lvr7esUddXd0nhlmaDtydnSCsZ8+ebR6Ihw0b5tdBu62Dd0xMDL169UKkte+65/fzzsvL4+jRo6SmptK7d+9u+Vx/uJkU1gL3isgaPJOSn7TzCcY4S1U5efIkZWVllJeXNx/gWx7ofZebnnf2mzNAjx49uvTo168fPXr04MILL2z3gNzRwbtfv3707NnTgZ9o9zt58iQZGRlMnjyZqVOnBnz7jiUFEXkZuA4YIiKlwA+BngCquhx4C/gUUAxUAf/uVCzGBEpNTU2nxoq7S11dHUeOHKGsrKzdx5EjR9qMT0QYOHAgcXFxDB48mMGDBzN27Njm54MHD/7Ea3FxcQwaNIi+ffu2ekCPiorqtm/OkUBV2bx5M5WVlSxcuNC1JObk1Ud3dPC6At92avvGOKm2tpbCwkLy8/PZvn07+fn55OfnU1JS0ulhDadER0czbNgwhg8fzvDhw7n88subnw8fPpxhw4Y1H+Tj4uIYMGAA0dHRbocdkY4dO0Z2djbTpk1j1KhRrsYScqWzjXFSQ0MDJSUllJWVcfjw4Ta/ce/du5eGhgbAM0Ry6aWXEh8fz5IlSxgwYEDA426ZAJoO+IE4MWnOX2NjI9nZ2TQ2NpKcnBwUSdmSgol4u3fv5p133uGdd95hw4YNVFRUnNNmyJAhzQfc2bNns3jxYqZMmcKVV17JJZdcQq9evVyI3ISyw4cPs3nzZmbNmsXQoR3OfRMwlhRMxFBVTp8+zaFDh8jOzua9997j/fffp6SkBIAxY8Zw8803c/XVVzN69OjmJDB06FB69LA/FdM96uvryczMpHfv3qSmpgZdb85+003YOXr0KHl5eWzbto1t27axY8cODh06RFlZGTU1Nc3t4uLimD9/Pg888ACLFi3i0ksvtROjxlH79+8nPz+fOXPmMHjwYLfDaZUlBRMWcnNzefXVV3n99dfZvn178/qhQ4cyZcoUrr322k+cZJ0yZQpTp04Num9pJjydPXuW9PR0Bg0aRGpqalB/+bCkYELasWPHePjhh3nhhReIjo7m2muvZcmSJUyfPp2pU6cyfPhwt0M0EW7Xrl0UFxeTmJhIbGys2+F0yJKCCRkNDQ0UFhayfft2SkpK2LNnD6+++ioVFRU88sgjPPTQQwwZMsTtMI0BoLq6mrS0NEaNGkVycnJQ9w58WVIwQa2wsJAXXniBrKwscnJyOHPmTPNrgwcPZubMmTz11FOu3PlpTFt27tzJgQMHSExMpF+/fm6H0ymWFExQys/P5/HHH+eVV16hR48eTJ8+nS9/+cvMmjWLadOmMWHCBAYOHOh2mMZ8wpkzZ0hLS2PixIkkJSW5Hc55saRggsrWrVv5yU9+wmuvvUb//v155JFHePDBB4PqOm5jWlJVtm3bxvHjx7nuuusCWsCuu1lSMK5QVQ4ePEhBQUHzIzc3lw8//JCBAwfygx/8gPvvv5+4uDi3QzWmXSdOnCAzM5PLLruMadOmuR1Ol1lSMAH37rvv8q1vfYuPP/64eV1sbCyTJ09m2bJl3HfffQwaNMjFCI3pmKqyadMmqqurSUpKCpsbHMNjL0xQUVX27t3Lxo0bKS0t/cRrOTk5rF69mkmTJvHss89yxRVXMHnyZEaOHBkyV2cYc/ToUbKzs5k+fTojR450O5xuZUnBdJv33nuPp59+mqysLMrLy1tt07NnT5YuXcrSpUvp27dvgCM0pmsaGhrIzs4GICUlJSgK2HU3SwqmXUePHmX//v2fqBDqWyqiSUFBAS+//DJjxozhxhtvJCEhgYSEBCZOnPiJHkDv3r3p06dPIHfBmG5x6NAhtmzZQkJCQljfD2NJwZyjqqqK1157jRdeeIENGzb49Z7evXvz/e9/n8ceeyzkrss2pj319fVkZGTQt29fbrjhhrAf5rSkYD4hKyuL1NRUTpw4wYQJE/jhD3/YXC6i6RETE9Pqe62OkAk3+/btY/v27cydOzdiLn6wpGAoKSnh9ddf5/XXX+ejjz5i7NixvP7661xzzTV2oDcRqba2lvT0dOLi4oK+gF13s6QQoerr6/n1r3/NypUryc3NBWDKlCk89thjfO1rX+PCCy90OUJj3FFcXMyuXbuYN28e/fv3dzucgLOkEKHuu+8+li9fzpw5c3jyySe5+eabueiii9wOyxjXNBWwGz16NCkpKW6H4xpLChFo5cqVLF++nIceeognnngiorrGxrRmx44dHDx4kHnz5kX8pdKWFMJcXV0dv/jFL3j++eepq6sD4ODBgyxatIif/exnlhBMRDt9+jTp6elMmjQpZAvYdTdLCmGqoaGBDz74gAceeIDc3FyuvfZaJk6cCMCgQYN47LHHwvLGG2P8oark5uZy4sQJFixYQK9evdwOKWhYUggTu3btYsGCBRw/fpyePXty9uxZqqqqGDVqFH/+85+55ZZbrFdgDFBRUUFmZiZXXHEFV111ldvhBB1LCiFu3759/N///R9Lly4FICkpicsvv5yoqChmz57NjTfeGJFXUBjTUmNjI5s2baKmpoZFixaFTQG77mY/lRD30EMP8eqrr3LJJZfwwAMP8I1vfMPtkIwJOuXl5WzcuJEZM2YwYsQIt8MJapYUQlBDQwMZGRm88cYb/OlPf+LTn/40b775ptthGRN0GhoayMrKIioqitTUVLsZ0w+WFEJMWloa99xzD9u3bwc8lRrXrFnjclTGBJ+DBw+Sm5tLQkICF1xwgdvhhAxLCkFOVXnppZfYt28fRUVF/P73v2fcuHG89NJLLFq0iGHDhtkJZGN81NXVkZGRQUxMTMSVqOgOlhSCVGNjIz/4wQ/44IMP+OijjwDPXAQPPvggP/7xj+3ksTGt2Lt3Lzt37mTu3LkMHDjQ7XBCkiWFIHTvvffy0ksvcerUKaZOncrDDz/M448/TlRUlN1bYEwramtrSUtLY8iQIaSkpFjvoAscTQoikgo8A0QDv1XVn7V4fSzwIjDI2+ZRVX3LyZiCXWNjIy+99BKTJk3i61//Ovfcc4/9ghvTjqKiIkpKSkhMTLQedDdwLCmISDTwHLAIKAU2ishaVd3h0+wx4I+q+r8icjnwFjDeqZiCVWNjI7m5udTX17N3715OnTrF/fffz5e+9CW3QzMmaFVVVZGWlsaFF15IcnKy2+GEDSd7CglAsaruBhCRNcBNgG9SUGCA9/lA4KCD8QQlVeXOO+885wqiq6++2qWIjAluqsqOHTs4fPgw11xzjU3v2s2cTAqjgf0+y6XA7BZtfgSsF5H7gBig1YpUInI3cDfA2LFjuz1QtxQUFPD000+zZs0aHn74YebPnw/A0KFDmTBhgsvRGRN8mgrYXXzxxSxcuNDtcMKSk0mhtYFwbbF8B7BKVZ8SkbnASyJypao2fuJNqiuAFQDx8fEtPyMknThxgoSEBKqqqrjnnnv47//+bzt3YEwbVJWtW7dy8uRJK2DnMCeTQingO33XGM4dHroLSAVQ1QwR6QMMAY44GFdQeOGFFzh9+jRZWVkkJCS4HY4xQaupgN2UKVOYPn262+GEPSeTwkbgYhGZABwAFgNfaNFmH7AQWCUilwF9gHIHYwoK77//PkuXLmXhwoWWEIxpQ2NjIzk5OZw9e9YK2AWQYz9lVa0XkXuBdXguN12pqttFZBmQo6prge8CvxGRB/AMLX1FVcNieKgtf/vb37jtttuYOHEiq1evdjscY4LSkSNHyMnJYebMmQwfPtztcCKKhNoxOD4+XnNyctwO47z87W9/4+abb2batGm8/fbbDBkyxO2QjAkqDQ0NZGZm0qNHD2bNmmUF7LqRiGxS1fiO2ll/LIBefPFFYmJi2LBhAwMGDOj4DcZEkAMHDrBt2zZmz55NXFyc2+FELEsKAVRfX8/o0aMtIRjjo66ujvT0dGJjY62AXRCwpBAgxcXFvPnmm3z2s591OxRjgsaePXsoKCggMTHRviwFCUsKAVBWVsY999xDXV0dDz74oNvhGOO6mpoa0tPTGTp0qBWwCzKWFBy0Z88ennzySX73u99RX1/PY489Rnx8h+d5jAlrH3/8MXv37iUxMZGYmBi3wzEtWFLoZmfPnmXDhg289tprvPjii6gqS5Ys4Xvf+x6XXHKJ2+EZ45rKykrS09MZO3YsixYtcjsc0wZLCt1ow4YNLF68mPLycmJiYrjzzjtZtmwZY8aMcTs0Y1yjqmzfvp0jR45YAbsQYBcBd5NVq1aRkpLC0KFD+etf/8rRo0dZuXKlJQQT0U6dOsW6deuIiYnh+uuvt4QQAqyn0E2++93vMmvWLN5++227isJEPFVly5YtnD59moULF9KzZ0+3QzJ+8ispiEgvYKyqFjscT0jatGkTJ06cICkpyRKCiXjHjx8nKyuLqVOnMmPGDLfDMZ3UYVIQkU8DTwO9gAkichXwQ1W9xengQkF5eTm33norY8aM4b777nM7HGNc09jYyMaNG6mvryc5OdnmEw9R/vQUluGZHOc9AFXdKiKTHI0qRNTX17N48WLKyspIS0tj6NChbodkjCvKysrYtGkT8fHxDBs2zO1wTBf4kxTqVPVEi5tLQquKnkNeeeUVNmzYwIoVK5g5c6bb4RgTcPX19WRmZtKrVy9SU1OtgF0Y8Ccp7BSRzwNR3rkR/h+Q6WxYoaGiogKAW26xkTQTeUpLS8nLy7MCdmHGn7R+LzATaAReA2rwJAZjTAQ6e/Ys77//PuXl5aSmplpCCDP+9BRSVPUR4JGmFSJyK54EEdFCbS4KY7pq9+7dFBUVMXfuXLvSLkz501N4rJV1S7s7kFBTX1/PqlWrGD58OAMHDnQ7HGMcVVNTw7vvvkt1dTXJycmWEMJYmz0FEUkBUoHRIvK0z0sD8AwlRbTnnnuOzZs388c//tFuzDFhraCggP379zNv3jz69evndjjGYe0NHx0B8vGcQ9jus/408KiTQYWC5cuXc80113Dbbbe5HYoxjqisrCQtLY3x48dbAbsI0mZSUNUtwBYR+YOq1gQwpqBXUFBAQUEB9913n9WBN2FHVcnLy+Po0aPMnz+f3r17ux2SCSB/TjSPFpHHgcuB5mpWqhqxdaBLSkoAmD59usuRGNO9Tp48SUZGBpMnT2bq1Kluh2Nc4E9SWAX8FHgSuAH4dyL8nML69euJiopi9OjRbodiTLdQVTZv3kxlZaUVsItw/lx91E9V1wGo6i5VfQxY4GxYwWv16tU888wzfPWrX2Xs2LFuh2NMlx07doy3336bkSNHcu2111pCiHD+9BRqxTNwvktEvgEcACKyuMnatWtZsmQJ8+fP59lnn3U7HGO6pLGxkezsbBobG62AnWnmT1J4AOgPfAd4HBgIfNXJoIKRqvLggw9y5ZVXsnbtWvr27et2SMact8OHD7N582ZmzZplhRzNJ3SYFFQ1y/v0NPAlABGJuOnE8vLy2LVrF88//zyxsbFuh2PMeamvrycjI4M+ffpYATvTqnaTgojMAkYDH6nqURG5Ak+5i+uBiEkMtbW1PPqo59aMlJQUl6Mx5vzs37+f/Px85syZw+DBg90OxwSpNr8miMh/AX8A7gTeFpGleOZUyAUi5nLUM2fOkJSUxN///neWLVvGuHHj3A7JmE5pKmB37NgxUlNTLSGYdrXXU7gJmKaq1SISBxz0LhcGJrTg8MQTT/DRRx+xevVq7rjjDrfDMaZTdu3aRXFxMYmJiTbsafzSXlKoUdVqAFU9LiIFkZYQfvnLX/LTn/6U22+/3RKCCSnV1dWkpaUxatQokpOT7c5747f2ksJEEWkqjy3AeJ9lVPXWjj5cRFKBZ4Bo4Leq+rNW2nwe+BGe2dxyVfUL/ofvjMbGRr73ve/x1FNPcdNNN7Fy5Uq3QzLGbzt37uTAgQMkJiZaATvTae0lhc+2WP5VZz5YRKKB54BFQCmwUUTWquoOnzYXA98H5qlqhYgExf0Pd911F6tWreLb3/42zzzzjF2/bULCmTNnSEtL46KLLiIpKcntcEyIaq8g3rtd/OwEoFhVdwOIyBo85yl2+LT5OvCcqlZ4t3mki9vssj179rBq1Sruv/9+nn76aet2m6Cnqmzbto3jx49z3XXXWQE70yVOXqQ8Gtjvs1zqXefrEuASEUkTkUzvcNM5RORuEckRkZzy8nKHwvV44403APj2t79tCcEEvRMnTrBu3ToGDx7MggULLCGYLvPnjubz1doRteX8lT2Ai4Hr8Nz38KGIXKmqJz7xJtUVwAqA+Ph4R+fA3LNnD7GxsUyaNMnJzRjTJarKpk2bqK6uJikpiR49nPxTNpHE798kEemtqrWd+OxS4EKf5TF4Lmtt2SZTVeuAEhEpxJMkNnZiO93OeggmmB09epTs7GymT5/OyJEj3Q7HhJkOh49EJEFE8oAi7/I0EfmlH5+9EbhYRCaISC9gMbC2RZvX8VZcFZEheIaTdnci/m7X0NDg5uaNaVNDQwMZGRkUFRWRkpJiCcE4wp+ewrPAjXgO4Khqroh0WDpbVetF5F5gHZ5LUleq6nYRWQbkqOpa72vJIrIDaAAeVtVj57kv3aKoqIgJEya4GYIx5zh06BBbtmwhISGBIUOGuB2OCWP+JIUoVd3bYkjFr6/TqvoW8FaLdf/p81yBB72PoLB161arb2SCRlMBu759+3LDDTfY0KZxnD9JYb+IJADqvffgPuBjZ8Nyx+HDhzl8+DBXXXWV26EYw759+9ixYwdz5sxh0KBBbodjIoQ/SeGbeIaQxgJlwD+868JObm4ugCUF46ra2lrS09OJi4sjJSXFegcmoPxJCvWqutjxSILA1q1bAZg2bZrLkZhIVVxczO7du0lMTKR///5uh2MikD83r20UkbdE5MsiEtZlFrdu3cq4ceOstLAJuKqqKt555x3q6upITk62hGBc48/MaxeJSCKeS0p/LCJbgTWqusbx6AKsuLiYSy+91O0wTARRVXbu3MnBgwe5+uqrbZpX4zq/ylyoarqqfgeYAZzCM/lO2FFVevbs6XYYJkKcPn2a9evX07t3b5KSkiwhmKDQYU9BRPrjKWS3GLgMeANIdDguV3iukDXGWapKbm4uJ06cYMGCBfTq1cvtkIxp5s+J5nzgr8ATqvqhw/G4pqKiovnyP2OcUlFRQWZmJldeeaVd5WaCkj9JYaKqNjoeicv+8Ic/UFNTw1133eV2KCYMNTY2smnTJmpra1m0aJEVsDNBq83fTBF5SlW/C/xZRM4ZV/Fn5rVQoaqsWLGCmTNnMmPGDLfDMWGmvLycjRs3MmPGDEaMGOF2OMa0q72vK694/+3UjGuhKCsri7y8PJYvX+52KCaMNDQ0kJWVRVRUFKmpqURFOTl9iTHdo72Z17K9Ty9T1U8kBm+hu67OzBY0fvOb3xATE8MXvuD69NAmTBw8eJDc3FwSEhK44IIL3A7HGL/589Xlq62sC5uB95MnT7JmzRruuOMOYmPD+t48EwB1dXV88MEHHDp0iNTUVEsIJuS0d07hdjyXoU4Qkdd8XooFTrT+rtCzevVqqqqquPvuu90OxYS4vXv3snPnTubOncvAgQPdDseY89LeOYVs4BieGdOe81l/GtjiZFCB9MILL3DVVVcRHx/vdigmRNXW1pKWlsaQIUOsgJ0Jee2dUygBSvBURQ1LZ8+eZfPmzTzyyCP2h2zOS1FRESUlJcybN4+YmBi3wzGmy9o8pyAi//T+WyEix30eFSJyPHAhOqe4uJiGhgYuu+wyt0MxIaapgF1DQwPJycmWEEzYaG/4qGnKzbCd+2/nzp0AlhSM31SV7du3U1ZWxjXXXEOfPn3cDsmYbtVmT8HnLuYLgWhVbQDmAvcAYfG1qKCgAMAqoxq/NBWw69evHwsXLrSEYMKSP/favw7MEpGLgN8DfwNWAzc6GVggFBQUMGbMGKtdb9qlqmzdupVTp05ZATsT9vxJCo2qWicitwK/UNVnRSQsrj4qLCxk8uTJbodhglhTAbspU6Ywffp0t8MxxnF+TccpIp8DvgTc7F0X8pMOqCoFBQUsWbLE7VBMEGpsbCQnJ4ezZ89aATsTUfz5Tf8q8C08pbN3i8gE4GVnw3Le4cOHOX36tPUUzDmOHDlCTk4OM2fOZPjw4W6HY0xA+TMdZ76IfAeYJCKTgWJVfdz50JxVWFgI2Elm8y8NDQ1kZmbSo0cPK2BnIpY/M69dA7wEHAAEGCEiX1LVNKeDc5JdeWR8HThwgG3btjF79mzi4uLcDscY1/gzfPRz4FOqugNARC7DkyRCui5EYWEh/fr1Y8yYMW6HYlxUV1dHeno6sbGxpKam2p3tJuL5kxR6NSUEAFXdKSIhf01eQUEBl1xyiQ0RRLCSkhIKCwtJTExkwIABbodjTFDwJylsFpHn8fQOAO4kDAriFRYWMnv2bLfDMC6oqakhPT2dYcOGWQE7Y1rwJyl8A/gO8D085xQ+AH7pZFBOq6mpYc+ePXY5agQqLCxk3759JCYmWr0iY1rRblIQkSnARcBfVPWJwITkvKKiIlTVLkeNIJWVlaSnpzNu3DgWLVrkdjjGBK32Jtn5DzwzrG3GU+ZimaquDFhkDvr4448Bu/IoEqgq+fn5lJeXWwE7Y/zQ3lnWO4Gpqvo5YBbwzc5+uIikikihiBSLyKPttLtNRFREAnJF0/HjnsrfQ4cODcTmjEtOnTrFunXr6N+/P9dff70lBGP80N7wUa2qVgKoarmIdOoyHRGJxjNj2yKgFNgoImt9r2TytovFc84iq1ORd0F1dTUA/fr1C9QmTQCpKlu2bOH06dMsXLiQnj1DviqLMQHTXlKY6DM3swAX+c7VrKq3dvDZCXjuft4NICJrgJuAHS3a/QR4AnioM4F3RVVVFWBJIRwdO3aM7Oxspk6dyowZM9wOx5iQ015S+GyL5V918rNHA/t9lkuBT1wDKiLTgQtV9U0RaTMpiMjdwN0AY8eO7WQY56qqqkJE6N27d5c/ywSHxsZGNm7c2DwTWnR0tNshGROS2puj+d0ufnZrF39r84ue4aifA1/p6INUdQWwAiA+Pl47aN6h6upq+vbta9enh4mysjI2bdpEfHw8w4YNczscY0Kak/WAS/HM2tZkDHDQZzkWuBJ433uFZMAPAAATB0lEQVRwHgGsFZHPqGqOg3FRVVVF3759ndyECYD6+noyMzPp1auXFbAzpps4mRQ2Ahd7S20fABYDX2h6UVVP4jP/s4i8DzzkdEIAT0/BzieEttLSUvLy8pgzZw6DBw92OxxjwobfSUFEeqtqrb/tVbVeRO4F1gHRwEpV3S4iy4AcVV3b+XC7R1VVlSWFEHX27FnS09MZOHCgFbAzxgH+lM5OAH4HDATGisg04Guqel9H71XVt4C3Wqz7zzbaXudPwN3Bho9C0+7duykqKmLu3LlWwM4Yh/gzCPsscCNwDEBVc4EFTgblNBs+Ci01NTW8++67VFdXk5ycbAnBGAf5M3wUpap7W3TTGxyKJyBs+Ch0FBQUUFpaSmJiov2fGRMA/vQU9nuHkFREokXkfuBjh+NylA0fBb8zZ86wfv16oqOjSUpKsoRgTID401P4Jp4hpLFAGfAPzqMOUjCx4aPgpark5eVx7Ngx5s+fbzcYGhNgHSYFVT2C53LSsGHDR8Hp5MmTZGRkcNlllzF16lS3wzEmIvlz9dFv8LkTuYmq3u1IRAFgw0fBRVXZvHkzlZWVVsDOGJf5M3z0D5/nfYBb+GRNo5Bjw0fBo6mA3bRp0xg1apTb4RgT8fwZPnrFd1lEXgLecSwih6mq9RSCQGNjI9nZ2TQ2NloBO2OCyPmUuZgAjOvuQAKlpqYGsLLZbjp8+DBbtmwhPj7eJjoyJsj4c06hgn+dU4gCjgNtzqIW7GyCHffU19eTkZFBnz59SElJsQJ2xgShdpOCeO5Ym4anoB1Ao6p2uXS1m5om2LHho8Dav38/+fn5VsDOmCDXblJQVRWRv6jqzEAF5DTrKQRWUwG7wYMHWwE7Y0KAP+cUskVkhqpudjyaALCpOANn165dFBcXk5iYSGxsrNvhGGP80GZSEJEeqloPXA18XUR2AZV4ZlRTVQ3JCXBt+Mh51dXVpKWlMWrUKFJSUtwOxxjTCe31FLKBGcDNAYolIGz4yFk7d+7kwIEDzJs3zxKvMSGovaQgAKq6K0CxBIQNHznjzJkzpKenM3HiRJKSktwOxxhzntpLCkNF5MG2XlTVpx2Ix3E2fNS9VJVt27Zx/PhxK2BnTBhoLylEA/3x9hjChQ0fdZ8TJ06QkZHBFVdcwbRp09wOxxjTDdpLCodUdVnAIgkQ6yl0naqyadMmqqurWbRoET16nM+N8caYYNThOYVwY+cUuubo0aNkZ2czffp0Ro4c6XY4xphu1l5SWBiwKAKoafjIegqd09DQQHZ2NgApKSlWwM6YMNVmUlDV44EMJFCqqqro0aOH1ezvhEOHDrFlyxYSEhIYMmSI2+EYYxwUcYPBNuua/+rq6sjMzKRv377ccMMNVqLCmAgQcUnBJtjxz759+9ixYwdz5sxh0KBBbodjjAmQiEsKNsFO+2pra0lPTycuLo6UlBTrHRgTYSIuKVhPoW3FxcXs3r2bxMRE+vfv73Y4xhgXRNwsJ3ZO4VxVVVX84x//oK6ujuTkZEsIxkSwiOsp2PDRv6gqO3fu5NChQ1bAzhgDRGBPwYaPPE6fPs369evp3bs3CxcutIRgjAEitKcwevRot8NwjaqSm5vLiRMnWLBgAb169XI7JGNMEHG0pyAiqSJSKCLFIvJoK68/KCI7RGSbiLwrIuOcjAci+5xCRUUFb7/9NhdccAHXXXedJQRjzDkc6ymISDTwHLAIKAU2ishaVd3h02wLEK+qVSLyTeAJ4HanYoLIHD5qbGxk06ZN1NbWWgE7Y0y7nDw6JADFqrobQETWADcBzUlBVd/zaZ8JfNHBeIDIO9FcXl5OTk4O06dPZ8SIEW6HY4wJck4mhdHAfp/lUmB2O+3vAv7e2gsicjdwN8DYsWO7FFSkDB81NDSQlZVFdHQ0KSkpREVF3DUFxpjz4GRSaO1WWG21ocgXgXhgfmuvq+oKYAVAfHx8q5/hj/r6eurq6sI+KRw8eJDc3FwSEhK44IIL3A7HGBNCnEwKpcCFPstjgIMtG4lIErAUmK+qtQ7GE/Zls+vq6sjIyCAmJobU1FQrUWGM6TQnk8JG4GIRmQAcABYDX/BtICLTgeeBVFU94mAsQHhPsLNnzx4KCgqYO3cuAwcOdDscY0yIciwpqGq9iNwLrMMz3/NKVd0uIsuAHFVdC/wPnnmgX/V+q92nqp9xKqZwnJ+5traWtLQ0hg4dagXsjDFd5ui1iar6FvBWi3X/6fM8ycnttxRu8zMXFRVRUlLCvHnziImJcTscY0wYiKhLUsKlp1BVVcU777xDQ0MDycnJlhCMMd0mou5iCvWegqqyY8cOysrKuOaaa+jTp4/bIRljwkxE9RRC+UTzqVOnWL9+Pf369eP666+3hGCMcURE9RRCcfhIVdm6dSunTp2yAnbGGMdFVFIIteGj48ePk5WVxZQpU5g+fbrb4RhjIkBEJoVg7yk0NjaSk5NDXV2dFbAzxgRURB1tQmH46MiRI2zatIkZM2YwfPhwt8MxxkSYiEoKwTx81NDQQGZmJj169LACdsYY10RkUgi2K3cOHDjAtm3bmD17NnFxcW6HY4yJYBGVFKqrq+nbt2/QfAuvq6sjPT2d2NhYK2BnjAkKEZUUgmmCnZKSEgoLC0lMTGTAgAFuh2OMMUAEJgW3TzLX1NSQnp7OsGHDrICdMSboRFRSaBo+ckthYSH79u0jMTHR6hUZY4JScAyuB4hbPYXKykrWr1+PiLBo0SJLCMaYoBVxPYVAJgVVJT8/n/Lycq699tqgu+rJGGNaiqikEMgTzSdPniQjI4PJkyczZcqUgGzTGGO6KuKSwuDBgx3dhqqyZcsWzpw5w8KFC+nZs6ej2zPGmO4UUUnB6eGjY8eOkZ2dzdSpU5kxY4Zj2zHGGKdEVFJwaviosbGRjRs3Ns+EFh0d3e3bMMaYQIi4pNDdPYWysjI2b97MzJkzGTZsWLd+tjHGBFpEJYXuHD6qr68nMzOTXr16WQE7Y0zYiJikoKrdNnxUWlpKXl4ec+bMcfzEtTHGBFLEJIXa2lpUtUs9hbNnz5Kens7AgQOtgJ0xJixFTFJommDnfHsKu3fvpqioiMTERGJjY7szNGOMCRoRkxTOdyrOmpoa0tLSGDFiBMnJydY7MMaENUsK7SgoKKC0tJTExETXq6saY0wgRMwlM50ZPjpz5gzr168nOjqapKQkSwjGmIhhPQUfqkpeXh7Hjh1j/vz59O7dO1DhGWNMUIiYpNDUU2grKTQVsLvsssuYOnVqIEMzxpigETFJoamn0HL4SFXZvHkzlZWVVsDOGBPxIi4p+PYUjh49ysaNG5k2bRqjRo1yKzRjjAkajp5oFpFUESkUkWIRebSV13uLyCve17NEZLxTsfgOHzU2NpKZmUlRURHJycmWEIwxxsuxpCAi0cBzwA3A5cAdInJ5i2Z3ARWqOgn4OfDfTsXT1FOorKzk7bff5qKLLmLu3LlW0dQYY3w42VNIAIpVdbeqngXWADe1aHMT8KL3+Z+AheLQ3WGVlZVcffXVVFRUkJqaytChQ53YjDHGhDQnk8JoYL/Pcql3XattVLUeOAlc0PKDRORuEckRkZzy8vLzCmbSpEmMGjWKuXPnWkVTY4xpg5NHx9a+8et5tEFVV6hqvKrGn+83/JtuuolXXnnFri4yxph2OJkUSoELfZbHAAfbaiMiPYCBwHEHYzLGGNMOJ5PCRuBiEZkgIr2AxcDaFm3WAl/2Pr8N2KCq5/QUjDHGBIZj9ymoar2I3AusA6KBlaq6XUSWATmquhb4HfCSiBTj6SEsdioeY4wxHXP05jVVfQt4q8W6//R5XgN8zskYjDHG+M8uwzHGGNPMkoIxxphmlhSMMcY0s6RgjDGmmYTaFaAiUg7sPc+3DwGOdmM4ocD2OTLYPkeGruzzOFXt8O7fkEsKXSEiOaoa73YcgWT7HBlsnyNDIPbZho+MMcY0s6RgjDGmWaQlhRVuB+AC2+fIYPscGRzf54g6p2CMMaZ9kdZTMMYY0w5LCsYYY5qFZVIQkVQRKRSRYhF5tJXXe4vIK97Xs0RkfOCj7F5+7PODIrJDRLaJyLsiMs6NOLtTR/vs0+42EVERCfnLF/3ZZxH5vPf/eruIrA50jN3Nj9/tsSLynohs8f5+f8qNOLuLiKwUkSMikt/G6yIiz3p/HttEZEa3BqCqYfXAU6Z7FzAR6AXkApe3aPMtYLn3+WLgFbfjDsA+LwD6eZ9/MxL22dsuFvgAyATi3Y47AP/PFwNbgMHe5WFuxx2AfV4BfNP7/HJgj9txd3GfrwVmAPltvP4p4O94Zq6cA2R15/bDsaeQABSr6m5VPQusAW5q0eYm4EXv8z8BC0WktalBQ0WH+6yq76lqlXcxE89MeKHMn/9ngJ8ATwA1gQzOIf7s89eB51S1AkBVjwQ4xu7mzz4rMMD7fCDnzvAYUlT1A9qfgfIm4PfqkQkMEpGR3bX9cEwKo4H9Psul3nWttlHVeuAkcEFAonOGP/vs6y483zRCWYf7LCLTgQtV9c1ABuYgf/6fLwEuEZE0EckUkdSARecMf/b5R8AXRaQUz/wt9wUmNNd09u+9UxydZMclrX3jb3ndrT9tQonf+yMiXwTigfmORuS8dvdZRKKAnwNfCVRAAeDP/3MPPENI1+HpDX4oIleq6gmHY3OKP/t8B7BKVZ8Skbl4ZnO8UlUbnQ/PFY4ev8Kxp1AKXOizPIZzu5PNbUSkB54uZ3vdtWDnzz4jIknAUuAzqloboNic0tE+xwJXAu+LyB48Y69rQ/xks7+/22+oap2qlgCFeJJEqPJnn+8C/gigqhlAHzyF48KVX3/v5ysck8JG4GIRmSAivfCcSF7bos1a4Mve57cBG9R7BidEdbjP3qGU5/EkhFAfZ4YO9llVT6rqEFUdr6rj8ZxH+Yyq5rgTbrfw53f7dTwXFSAiQ/AMJ+0OaJTdy5993gcsBBCRy/AkhfKARhlYa4El3quQ5gAnVfVQd3142A0fqWq9iNwLrMNz5cJKVd0uIsuAHFVdC/wOTxezGE8PYbF7EXedn/v8P0B/4FXvOfV9qvoZ14LuIj/3Oaz4uc/rgGQR2QE0AA+r6jH3ou4aP/f5u8BvROQBPMMoXwnlL3ki8jKe4b8h3vMkPwR6AqjqcjznTT4FFANVwL936/ZD+GdnjDGmm4Xj8JExxpjzZEnBGGNMM0sKxhhjmllSMMYY08ySgjHGmGaWFEzQEZEGEdnq8xjfTtvxbVWT7OQ23/dW4sz1loi49Dw+4xsissT7/CsiMsrntd+KyOXdHOdGEbnKj/fcLyL9urptExksKZhgVK2qV/k89gRou3eq6jQ8xRL/p7NvVtXlqvp77+JXgFE+r31NVXd0S5T/ivPX+Bfn/YAlBeMXSwomJHh7BB+KyGbvI7GVNleISLa3d7FNRC72rv+iz/rnRSS6g819AEzyvneht05/nrfOfW/v+p/Jv+aneNK77kci8pCI3IanvtQfvNvs6/2GHy8i3xSRJ3xi/oqI/PI848zApxCaiPyviOSIZx6FH3vXfQdPcnpPRN7zrksWkQzvz/FVEenfwXZMBLGkYIJRX5+ho7941x0BFqnqDOB24NlW3vcN4BlVvQrPQbnUW/bgdmCed30DcGcH2/83IE9E+gCrgNtVdQqeCgDfFJE44BbgClWdCvzU982q+icgB883+qtUtdrn5T8Bt/os3w68cp5xpuIpa9FkqarGA1OB+SIyVVWfxVMXZ4GqLvCWvngMSPL+LHOABzvYjokgYVfmwoSFau+B0VdP4FfeMfQGPDV9WsoAlorIGOA1VS0SkYXATGCjt7xHXzwJpjV/EJFqYA+e8suXAiWq+rH39ReBbwO/wjM/w29F5G+A36W5VbVcRHZ7a9YUebeR5v3czsQZg6fsg++sW58Xkbvx/F2PxDPhzLYW753jXZ/m3U4vPD83YwBLCiZ0PACUAdPw9HDPmTRHVVeLSBbwaWCdiHwNT5nhF1X1+35s407fgnki0uocG956PAl4irAtBu4Fru/EvrwCfB4oAP6iqiqeI7TfceKZgexnwHPArSIyAXgImKWqFSKyCk9huJYEeEdV7+hEvCaC2PCRCRUDgUPeGvlfwvMt+RNEZCKw2ztkshbPMMq7wG0iMszbJk78n5+6ABgvIpO8y18C/ukdgx+oqm/hOYnb2hVAp/GU727Na8DNeOYBeMW7rlNxqmodnmGgOd6hpwFAJXBSRIYDN7QRSyYwr2mfRKSfiLTW6zIRypKCCRW/Br4sIpl4ho4qW2lzO5AvIluByXimLNyB5+C5XkS2Ae/gGVrpkKrW4KlA+aqI5AGNwHI8B9g3vZ/3Tzy9mJZWAcubTjS3+NwKYAcwTlWzves6Haf3XMVTwEOqmotnbubtwEo8Q1JNVgB/F5H3VLUcz5VRL3u3k4nnZ2UMYFVSjTHG+LCegjHGmGaWFIwxxjSzpGCMMaaZJQVjjDHNLCkYY4xpZknBGGNMM0sKxhhjmv1/GLX8Zxg3tCwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(xROC,yROC,color='k')\n",
    "plt.plot([0,1],[0,1], lw=0.5, color='grey')\n",
    "plt.xlabel(\"False Positive Rate\")\n",
    "plt.ylabel(\"True Positive Rate\")\n",
    "plt.title(\"ROC Curve\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precision recall curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "xPR = []\n",
    "yPR = []\n",
    "\n",
    "for i in range(1,len(probY)+1):\n",
    "    preds = [x[1] for x in probY[:i]]\n",
    "    labs = [x[2] for x in probY[:i]]\n",
    "    prec = sum(labs) / len(labs)\n",
    "    rec = sum(labs) / sum(y[9000:])\n",
    "    xPR.append(rec)\n",
    "    yPR.append(prec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8VdWZ//HPE5IAgZAACRDuICg3FTX1WhG1VaCOOkpVfjpKW8VxhupYa1unv3qrtmOn1qplWql4vyClv0HkZlu81RFagngLiiIYCETuBMJFSPL8/tg7Zw4hJEfIuSXf9+t1XtmXdfZ+VgLnOWutvdc2d0dERAQgI9kBiIhI6lBSEBGRCCUFERGJUFIQEZEIJQUREYlQUhARkQglBRERiVBSkLRiZqVmNrqJMn3NrMrM2iQorLgys9FmVh61/pmZfS2ZMUnLpaQgzSL8oNoTfhhvMLPHzaxjc5/H3Ye7+2tNlFnj7h3dvaY5zmlmPes+lOvV83MzeyIe9TwSZnaymc0zs+1mttXM/m5m30p2XJIelBSkOf2Du3cETgS+Avzf+gUskG7/7sYBC6LW6+o5EjgBuC0pUTXAzE4DXgFeBwYBXYEbgLGHebwW0dqS2KXbf05JA+6+DpgPjAAws9fM7F4z+x9gNzDQzPLMbJqZVZjZOjO7J/oDyMyuM7MPzWynmS03sxPD7ZGuk/AbcYmZ7QhbJ78Kt/c3MzezzHC9p5nNDr81rzSz66LOc6eZzTCzp8JzlZpZcb0qjQPmNVDPz4GXCZJD3fHamtkvzWxNGNPvzKx91P6LzOydMOZPzWxMuP1bUfVdZWbXH+av/z+BJ939Pnff7IGl7n5ZeJ6JZvZm9BvC39WgcPkJM/tt2NLYBdwWtoii/zb/aGbvhcsZZvajsC5bwt9ll8OMXVKAkoI0OzPrQ/BBuixq8z8Bk4BcoAx4Eqgm+DZ7AnAecG34/m8CdwJXA52AC4EtDZzqQeBBd+8EHAXMOERIzwPlQE9gPPAzMzs3av+FwHQgH5gN/CaqLlnAKODPDdSzN8E38JVRm+8DjiZIFIOAXsDtYfmTgaeAW8NzjQI+C9+3EbggrO+3gAfqEmGszCwHOA2Y+WXe14D/A9xL8Lf6JbALOKfe/ufC5RuBi4GzCH6/24ApR3h+SSZ310uvI34RfLhVAdsJPvT/C2gf7nsNuDuqbHfgi7r94bYJwKvh8svATY2c52vh8hvAXUBBvTL9AQcygT5ADZAbtf/nwBPh8p3AX6L2DQP2RK2fCyxsoJ47w3MsBPLDfUbwAXpUVPnTgNXh8iPAAzH+PmfV/Q6A0UB5Q7+Deu/pFcY0pJHjTgTerLfNgUHh8hPAU/X23wM8Fi7nhnXsF65/CJwbVbYI2A9kJvvfpF6H91JLQZrTxe6e7+793P1f3H1P1L61Ucv9gCygIhwM3U7wgdkt3N8H+DSG832H4Fv5R2a2xMwuaKBMT2Cru++M2lZG8AFa5/Oo5d1Au7quJxruOrrY3XMJPqyHAAXh9kIgB1gaVa8F4fZG62VmY81scdjFtT08b0FDZRuxDagl+GA+EmvrrT8HXGJmbYFLgLfdvSzc1w/476j6fkiQhLsfYQySJEoKkijRc7SvJWgpFIRJJN/dO7n78Kj9RzV5QPdP3H0CQTK5D5hpZh3qFVsPdDGz3KhtfYF1McY9Dph7iPO/TvDN+pfhps3AHmB4VL3yPBiUhkPUK/yw/WN4nO7unk+QiCzGGOvi2Q0sAi5tpNgugsRVd+4eDR2q3nGXEyTSsRzYdQRBncZG1Tff3dt5MK4kaUhJQRLO3SuAPwH3m1mncLDyKDM7KyzyKPB9MzspvFppkJn1q38cM7vKzArdvZag2wqCb6nR51oLvAX83MzamdlxBC2MZ5uK08wGAG3d/aNGiv0a+LqZjQzj+D3BeEC38Bi9zOz8sOw04Ftmdm5Y515mNgTIBtoCm4BqMxtLMMZyOH4ATDSzW82saxjD8WY2Pdz/LjDczEaaWTuC7rNYPEcwfjAK+EPU9t8B99b9fcys0MwuOszYJQUoKUiyXE3wYbicoNtjJmG3h7v/gWCg8zmCvvtZQENXtIwBSs2simDQ+Qp339tAuQkE4wzrgf8G7nD3gwaOG/ANGrjqKJq7byIYPP5JuOmHBAPPi81sB/AX4Jiw7N8JB5GBSoLLRvuFXVs3EgyUbyP4Nj47hvgaiuctgkHhc4BVZrYVmFpXD3f/GLg7jOsT4M1DHKq+5wm6y15x981R2x8MY/2Tme0EFgOnHE7skhrMXU9eE2mImc0DfuPujSYGkZZELQWRQ3sNeDXZQYgkkloKIiISoZaCiIhEZDZdJLUUFBR4//79kx2GiEhaWbp06WZ3L2yqXNySgpk9RnDb/kZ3H9HAfiO4cmEcwQ1DE9397aaO279/f0pKSpo7XBGRFs3MypouFd/uoycILhk8lLHA4PA1CfhtHGMREZEYxK2l4O5vmFn/RopcRDDHihNc051vZkXhjU1xMW7cONaurX8Hf2KZGT/5yU/45je/mdQ4REQakswxhV4cOMdKebjtoKRgZpMIWhP07dv3sE84YMAA2rdv33TBOJozZw4LFy5UUhCRlJTMpNDQvC4NXh/r7lMJ7sqkuLj4sK+hnTIl+TP69ujR0FQzIiKpIZlJoZxg1sg6vQmmIZDD4O7s3LmTbdu2RV5FRUUMGTIk2aGJSBpJZlKYDUwOJ+o6BaiM53hCunF3KisrqaioYMOGDQ2+Nm7ceEASqKk58JHEhYWFbNy4MUk1EJF0FM9LUusm0Cqw4KHndxDMoY+7/45ggq5xBJOH7SaYKKxV+eKLL/jkk0/46KOPWL16NWVlZXz22WeUlZVRVlbGzp07D3pPmzZtKCwspEePHhQWFjJw4EA6d+58wKtLly4899xzzJv35absqUtEmzdvjry2bNnC5s2b2bZtG9u3b4/8rK6uZurUqfTrd9Dkpc2mtraWqqoqdu7cyY4dO6iqqmLo0KF07Nix6TeLyGGJ59VHE5rY78C/xuv8qeyll17ilVde4dNPP6W2tjayPS8vj/79+zNw4EDOPvts+vXrR1FREd27d4+8unbtSkZG01cSv/XWW0DwQb9t2zbWr18feVVUVER+btq0KZIAtm7dSnV1dYPHy8jIID8/n/z8fLKyslixYgUlJSWNJoWamhq2bdvGli1b2L59e+RVWVl5yPXo5V27dh10zGuuuYYnnniiyfqLyOFJuzua092oUaMoLS1l6NChXH755QwdOpShQ4cycOBA8vLymvVcu3fvpn379nzxxRcH7cvPz6eoqIhu3boxdOhQunbtSkFBwQGvum1du3alU6dOBPcbwvvvv89xxx3HwoUL+eyzz9i8efMBySU6yTQ2t1ZmZiZ5eXnk5+dHfvbo0YO8vDzy8vLo1KkTnTp1Ijc3l9zcXG677TaWLFnCHXfcQWVlJZdeeilnnnnmEf2O3D1SLxFJwwnxiouLXXc0N62kpIQpU6bQrVs3evbsSVFRET179owsH8mluWvWrDmghZCVlRVJJIWFhQ0ml7pWRnQCyMnJ+VIfyGPHjmXBggWR9UsvvZSZM2dSU1NDZWUl27ZtY+vWrQf8bGhbXYuksrKSHTt28JOf/IQ77rijwXO6O3v37qW6uprc3NwGy4ikAzNb6u7FTZZTUpDDsXr1atydgoICcnNzE/Jtu7q6mqqqKnJzcxk5ciQrV66kXbt2VFZWNtoiycnJOWDMpS4x5eXl8eyzz9K9e3fOPvvsSMKIHjvZvn07+/bto02bNixfvpyjjz467vUUiYdYk4K6j+SwDBgwIOHnzMzMJD8/H4CbbrqJ1157jS5dukQG16MH2qN/tm3b9pDHXLt2LfPnz2fbtm0HtGb69+8fWd64cSOPP/4448ePZ+jQoUyfPl1dTtJiqaUgrVrdv//GPuQrKiq4/PLLKSsrY82aNZGxGpF0opaCSAxi+cZfVFTEG2+8wX333cePfvQjVq5cyc6dO9m0aVNkgL1uuaqqinvvvZehQ4cmIHqR5qekIBKjzMzgv8txxx130L66cYt169ZxzjnnKClI2lJSEInR1VdfTbt27cjLy4tcZVX3Mycnh82bN1NYWMjKlSt58cUX2bBhA5s2beLiiy9m+PDhyQ5fJCYaUxBpJjt27CA/P/+gK6FGjx7Nww8/zIgRBz1rSiRhYh1T0DOaRZpJp06dePXVV5k7dy4lJSWsXbuWYcOG8dprr3H55Zfz0UcfJTtEkSappSASRxs2bGDixIksWLCAnJwcXn/9dT7//PMDph1Zv3491dXVvPDCC5G72nft2nXAdCQbN27kggsuaPBS4P3791NVVUXnzp0TXT1JI7p5TSRFrF+/nttuu42nnnrqgO1mRrdu3WjXrh1lZWWMHDmS3bt3U1FR0eBkiAA//OEPqaio4PPPP4/83LRpEwAPPfQQN9xwQ2RAXCSakoJICtm0aRMvvvgihYWFkSlHunfvTlZWFqtWrWL8+PF06NCBoqKiyP665aKiIsaMGcO6devIzs6mR48e9OjRg6KiInr06IG7M3XqVADOO+88brnlFkaPHk12dnaSay2pRElBpAXZvXs3e/fupXPnzg3eWzFjxgwuv/zyyPoLL7zAZZddlsgQJcVpoFmkBcnJyaFLly6HvNnusssu4+233+bRRx8F4I033khkeNKCqKUg0oLU3SsBsG7dOnr27JnkiCRVqKUg0goVFBRw4403AtCrVy9uvPFGKir0lFuJnVoKIi1MbW0tRUVFBzyfu3///mzdupV58+ZxxhlnJDE6SRa1FERaqYyMDMrLy6msrOTMM8/k9NNPZ8SIEezYsYNPP/002eFJitMFzSItUFZWFllZWZEB51WrVjFnzhxKS0t55ZVXOOecc5IcoaQqdR+JtALr1q2jd+/ekfWLL74YCO6GfuaZZ2jfvn2jDyOS9KfnKYhIRK9evZg3bx5vvfUW99xzD7NmzYrs69y5M4MHD+bjjz9OYoSSKjSmINJKjB07lp/+9KdUVVVRW1vLhg0b+O53v8tXvvIVPvnkEzZv3hwpW1tby9q1ayNTaEjroaQg0sp06NAhMu/SQw89xEUXXQRAYWEh3/jGNxgyZAjt27enb9++dOvWjdtvv50rr7ySUaNGMX/+/AaPuX///oOmDJf0pDEFkVZu3759dO/encrKSo4//niOOuooBg4cyIwZMygrKyMjI4O+ffvy2WefAfDb3/6Wbdu28emnn7Jq1SpWrVpFWVkZgwYN4sMPP9SEfClKcx+JSMzqPgeip9GoqamhvLycoqIisrOzGT9+PH/84x8j+7t3787AgQMZOHAgzz77LABdunRhzJgxPPHEE2RlZSW2EtKolEgKZjYGeBBoAzzq7v9Rb38/4DGgENgKXOXu5Y0dU0lBJDn279/PX//6VwoLCxkwYAAdO3aM7FuxYgVDhgyJrOfk5PCzn/0Md+emm2465JxNkjhJTwpm1gb4GPg6UA4sASa4+/KoMn8A5rj7k2Z2DvAtd/+nxo6rpCCSuqLnXqrzxhtvUFNTQ3Fx8QGJRBIrFe5oPhlY6e6r3H0fMB24qF6ZYcDCcPnVBvaLSBopKCigoqKCRYsW8etf/xqAUaNGcfbZZ/Pggw8mOTqJRTxHhHoBa6PWy4FT6pV5F7iUoIvpH4FcM+vq7lviGJeIxFHdQ4AGDx7Mrl276NGjB5MmTWLXrl3JDk1iEM+WQkOdiPX7qr4PnGVmy4CzgHVA9UEHMptkZiVmVqLrpkXSQ9euXfn3f/93vv3tb5ORoavf00U8/1LlQJ+o9d7A+ugC7r7e3S9x9xOAH4fbKusfyN2nunuxuxfX768UkfQwd+5cvvnNb1Je3ui1JJJk8UwKS4DBZjbAzLKBK4DZ0QXMrMDM6mK4jeBKJBFpYY466ijee+89Zs6cSZ8+fXSndAqLW1Jw92pgMvAy8CEww91LzexuM7swLDYaWGFmHwPdgXvjFY+IJE9paSnV1dXk5+cDwQ1wkprieuuhu88D5tXbdnvU8kxgZjxjEJHkqxtT2LhxI9nZ2SxZsiTJEcmhaPRHRBKm7i7nOXPmcMoppzB//nxeeOEFPvnkE9ydNWvW8PLLL7Ns2bIkR9p6aZISEUmoRx55hOuvv56///3vjBs3LrK9Y8eOVFVVAcFU3xqQTg7NfSQiCVdbW8u0adPo1KkTc+bMoaKiguHDhzN06FDmzJnD3LlzmTZtGtdccw0VFRX07NlTl7UeoaRPcxEvSgoiLdvy5csZPnw4nTt3Zvfu3XzxxRcAnHTSSezdu5e//e1vdOjQIclRpp9UmOZCRORLGzZsGFdeeSWnnnoqkydPBqBTp04sXbqU0tJS7r33Xt5++21qa2uTHGnLpJaCiKQ0d8fM2L59O507d45sf/311xk1alQSI0svaimISItQN+12fn4+Tz75JLfccgvAAYPU0nzUUhCRtFJdXU1WVhb5+fls27Yt2eGkDbUURKRFyszM5NprryUnJ4cdO3awYsUKPR+6GSkpiEhaWr9+PXl5eQwZMoTMzExOP/109u3bl+yw0p6SgoiknSuuuIKrrrqKm2++GQieF71o0SJOPfVU1qxZk+To0puSgoiknXPPPZenn36aX/3qV7g7s2bNYsiQISxbtowFCxZQXX3QY1kkRkoKIpL2Tj75ZJ5++mkArr/+eh555JEkR5S+lBREpEU44YQTIs+F1lVJh0+XpIpIi1F3uSrAvn37IsuiS1JFpBVq06YNHTt2BCA7O5vp06cnOaL0o6QgIi2GmfH666/To0cPACZMmMBLL72U5KjSi5KCiLQoJ554IhUVFTz++OMAvPfee0mOKL0oKYhIi3TllVcCUFZWxmOPPcbkyZO54YYbNLtqEzTQLCItUnV1NW3bto0kATPD3dm4cSOFhYVJji7xNNAsIq1aZmYmM2bM4Pnnn2fFihU8/PDDAHTr1o0f//jHSY4udamlICKtwubNm+nfvz979+6lT58+rF69OtkhJZRaCiIiUQoKCqiqqoqMNUjDlBREpNVJtx6SRFJSEJFWxd0pKyuje/fu7NmzJ9nhpBwlBRFpVSZPngzAxo0bycnJ4brrrktyRKlFSUFEWpWTTz6Z7du3M2TIEPLy8li4cGGyQ0opcU0KZjbGzFaY2Uoz+1ED+/ua2atmtszM3jMzPYlbROIuLy+PDz/8kAsvvDDZoaScuCUFM2sDTAHGAsOACWY2rF6x/wvMcPcTgCuA/4pXPCIiDVm9ejXHHHMMixcv5qmnnmLLli3JDimp4tlSOBlY6e6r3H0fMB24qF4ZBzqFy3nA+jjGIyJygLqWwscff8xpp53GNddcw/Dhw1m6dGmrvUIpnkmhF7A2ar083BbtTuAqMysH5gHfbehAZjbJzErMrGTTpk3xiFVEWqHx48ezdu1abr/9dp555hkANmzYQHFxMUuWLElydMkRz6RgDWyrn3onAE+4e29gHPC0mR0Uk7tPdfdidy9ujXOWiEj89O7dm7vuuosrr7ySdevWcdNNNwHw1ltvJTmy5IhnUigH+kSt9+bg7qHvADMA3H0R0A4oiGNMIiKH1LNnT7797W8DcPPNNzN//ny2b9+e5KgSK55JYQkw2MwGmFk2wUDy7Hpl1gDnApjZUIKkoP4hEUmaY489lnPPPReAcePGceuttyY5osSKW1Jw92pgMvAy8CHBVUalZna3mdVdB3YLcJ2ZvQs8D0z01jq6IyIpwcz485//zGOPPUa3bt148803ufrqqyktLU12aAmhWVJFRA7h9NNPZ9GiRZH16upq2rRpk8SIDp9mSRUROUILFixgw4YNkfVdu3YlMZrEUFIQETmETp060a1bN+6///5kh5IwSgoiIhKhpCAiIhFKCiIiEqGkICIiEUoKIiISoaQgIhKj888/n/bt21NeXp7sUOImM9aCZtYL6Bf9Hnd/Ix5BiYikkv79+5OVlcXixYsB6NOnDzU1NWRktLzv1TElBTO7D7gcWA7UhJsdUFIQkRbvkksuYe/evQCRO5p37dpFbm5uMsOKi1jT3MXAMe4+zt3/IXzpOXYi0mpkZGSQkZHBL3/5y2SHElexJoVVQFY8AxERkeSLdUxhN/COmS0Evqjb6O43xiUqERFJiliTwmwOfhaCiIi0MDF1H7n7kwTPO1gavp4Lt4mItEr5+fmceOKJzJ8/P9mhNKuYkoKZjQY+AaYA/wV8bGaj4hiXiEhK+upXvwpAbW0ty5YtY9y4cUyZMiXJUTWfWAea7wfOc/ez3H0UcD7wQPzCEhFJTaeccgo1NTXs37+fH/zgBwC89tpr1NTUNPHO9BBrUshy9xV1K+7+MboaSURaqYyMDDIzM7nvvvsYNGgQM2fOJDMzk6uuuirtk0OsSaHEzKaZ2ejw9XuCsQURkVbte9/7Hl26dAHg2WefJTMzk08//TTJUR2+WJPCDUApcCNwE8Gdzf8cr6BERNLFDTfcwJYtW3j77bcZNSoYan3ooYeSHNXhM3dPdgxfSnFxsZeUlCQ7DBGRg3z88cccc8wxAPzpT3/i61//epIj+l9mttTdi5sq12hLwcxmhD/fN7P36r+aK1gRkZbg6KOP5rnnngPg1ltvTXI0h6epm9duCn9eEO9ARERagksuuQSAd999l71799KuXbskR/TlNNpScPeKcHEzsNbdy4C2wPHA+jjHJiKSdtq2bcudd94JwH333Zd2z16IdaD5DaBd+EyFhcC3gCfiFZSISDo744wzALjzzjvp06cPBQUFpMv4baxJwdx9N3AJ8LC7/yMwLH5hiYikr6997WssWLCAiRMnArBlyxby8vJYu3ZtcgOLQcxJwcxOA64E5obbYn5qm4hIa3P++efz+OOP89Zbb3H88cezc+dO7rrrrmSH1aRYk8K/AbcB/+3upWY2EHi1qTeZ2RgzW2FmK83sRw3sf8DM3glfH5vZ9i8XvohIajvttNOYOXMmANOmTWPWrFlJjqhxsc6S+rq7X+ju94Xrq5p6loKZtSGYQG8sQVfTBDM7oMvJ3W9295HuPhJ4GPh/h1MJEZFUNmjQIObODTpZnnwytSeYbrQLyMx+7e7/ZmYvETyT+QBNPJLzZGClu68KjzUduIjgbuiGTADuiClqEZE0c9ZZZwEwa9Ys9u/fT1ZWak4f19S4wNPhz8N5KGkvIHpUpRw4paGCZtYPGAC8coj9k4BJAH379j2MUEREkqtDhw585zvfYdq0abz55pucffbZyQ6pQU3dp1A36V0J8NewG+l14E1gSRPHtoYOeYiyVwAz3b3B6QXdfaq7F7t7cWFhYROnFRFJTSeeeCIA55xzDh999FGSo2lYrAPNC4GcqPX2wF+aeE850CdqvTeHvuHtCoInu4mItFjXXnst48ePB+C6667jgw8+SHJEB4s1KbRz96q6lXA5p5HyELQkBpvZADPLJvjgP+g5z2Z2DNAZWBRjLCIiaSk7O5u7776b7Oxs3nzzzchVSakk1qSwy8xOrFsxs5OAPY29wd2rgcnAy8CHwIzwcta7zSx6gHoCMN3T5XY/EZEjMHToUD7//HMA7rrrrpR7KE9MU2eb2VeA6fxv908RcHnUmEPCaOpsEWkJRo4cybvvvssvfvGLhMyoGuvU2THdlezuS8xsCHAMwQDyR+6+/whjFBFptZ555hmOPfZYSktLkx3KAWLqPjKzHOCHwE3u/j7Q38w0nbaIyGEaMWIE+fn5PP/889TW1iY7nIhYxxQeB/YBp4Xr5cA9cYlIRKSV2L17N/v27WPZsmXJDiUi1qRwlLv/AtgP4O57aPg+BBERidHjjz8OwL59+5Icyf+KNSnsM7P2hDefmdlRwBdxi0pEpBXo2rVrskM4SKxJ4Q5gAdDHzJ4luJntB3GLSkSkFTn99NO54ILUGKZt8pJUMzOCu5F3A6cSdBstdvfN8Q/vYLokVURaisrKSgYOHMjWrVuBoBspXhPlxXpJapMthfCmslnuvsXd57r7nGQlBBGRliQvL48tW7Zwzz2pc91OrN1Hi8Mb2EREpAWL9ZGaZwP/bGafAbsIupDc3Y+LV2AiIpJ4sSaFsXGNQkREUkJTT15rB/wzMAh4H5gWTnQnIiItUFNjCk8CxQQJYSxwf9wjEhGRpGmq+2iYux8LYGbTgL/HPyQREUmWploKkZlQ1W0kItLyNdVSON7MdoTLBrQP1+uuPuoU1+hERCShGk0K7t4mUYGIiEjyxXrzmoiItAJKCiIiEqGkICIiEUoKIiISoaQgIiIRSgoiIinimmuuYc+ePUmNQUlBRCTJzj//fACef/55TjnllKTGoqQgIpJkxcXFLFq0CEAtBRERgVNPPZWLL76YnJycpMYR16RgZmPMbIWZrTSzHx2izGVmttzMSs3suXjGIyIijYv1ITtfmpm1AaYAXwfKgSVmNtvdl0eVGQzcBpzh7tvMrFu84hERkabFs6VwMrDS3Ve5+z5gOnBRvTLXAVPcfRuAu2+MYzwiItKEeCaFXsDaqPXycFu0o4Gjzex/zGyxmY2JYzwiItKEeCYFa2Cb11vPBAYDo4EJwKNmln/QgcwmmVmJmZVs2rSp2QMVEUkFNTU1vPfee2zfvj1pMcQzKZQDfaLWewPrGyjzorvvd/fVwAqCJHEAd5/q7sXuXlxYWBi3gEVEkik7OxuAzp0788EHHyQlhngmhSXAYDMbYGbZwBXA7HplZgFnA5hZAUF30qo4xiQikrIee+yxyPIDDzyQlBjilhTCx3dOBl4GPgRmuHupmd1tZheGxV4GtpjZcuBV4FZ33xKvmEREUlmnTp2ora0lLy8vaTHE7ZJUAHefB8yrt+32qGUHvhe+RERaPTMjNzc3aefXHc0iIhKhpCAikmJ27NjBq6++StCZklhKCiIiKWbnzp2sXr2a0tLShJ9bSUFEJMXUXXm0e/fuhJ9bSUFEJMUMHnzQ7VoJo6QgIpKiamtrE35OJQURkRSTkRF8NJ922mkJ70JSUhARSTFnnnlmZPmkk05K6LmVFEREUkyHDh1Ys2YNAB999BE1NTUJO7eSgohICurTpw933313ws+rpCAiIhFKCiIiEqGkICIiEUoKIiIk2fLbAAAImklEQVQpbu7cuQk7l5KCiEiKqrscddasWQk7p5KCiEiKGjduHL169YrczJYISgoiIhKhpCAiksK2bt3KtGnT2L59e0LOp6QgIpLC+vbtC8Ann3ySkPMpKYiIpLD7778fgJ///OcJOZ+SgohIChsxYgQA77//fkLOp6QgIpLC+vXrx4UXXkjbtm0Tcj4lBRGRFLdnzx5KS0uZM2dO3M+lpCAikuImTZoEQHl5edzPpaQgIpLivvrVrybsXEoKIiISoaQgIiIRmfE8uJmNAR4E2gCPuvt/1Ns/EfhPYF246Tfu/mi84nnnnXcSdlegiEhzqaqqYvTo0ezZsyfu54pbUjCzNsAU4OtAObDEzGa7+/J6RV9w98nxiiPayJEjE3EaEZFmVVtby1lnnZWQy1Lj2X10MrDS3Ve5+z5gOnBRHM8nItIiZWRkkJubS3Z2dvzPFcdj9wLWRq2Xh9vqu9TM3jOzmWbWp6EDmdkkMysxs5JNmzbFI1YRESG+ScEa2Ob11l8C+rv7ccBfgCcbOpC7T3X3YncvLiwsbOYwRUSkTjyTQjkQ/c2/N7A+uoC7b3H3L8LV3wMnxTEeERFpQjyTwhJgsJkNMLNs4ApgdnQBMyuKWr0Q+DCO8YiISBPidvWRu1eb2WTgZYJLUh9z91IzuxsocffZwI1mdiFQDWwFJsYrHhERaZq51+/mT23FxcVeUlKS7DBERNKKmS119+KmyumOZhERiVBSEBGRCCUFERGJUFIQEZEIJQUREYlQUhARkQglBRERiVBSEBGRCCUFERGJUFIQEZEIJQUREYlQUhARkQglBRERiVBSEBGRCCUFERGJUFIQEZEIJQUREYlQUhARkQglBRERiVBSEBGRCCUFERGJUFIQEZEIJQUREYlQUhARkQglBRERiVBSEBGRCCUFERGJiGtSMLMxZrbCzFaa2Y8aKTfezNzMiuMZj4iINC5uScHM2gBTgLHAMGCCmQ1roFwucCPwt3jFIiIisYlnS+FkYKW7r3L3fcB04KIGyv0U+AWwN46xiIhIDDLjeOxewNqo9XLglOgCZnYC0Mfd55jZ9w91IDObBEwKV6vMbMVhxlQAbD7M96Yr1bl1UJ1bhyOpc79YCsUzKVgD2zyy0ywDeACY2NSB3H0qMPWIAzIrcfdWNW6hOrcOqnPrkIg6x7P7qBzoE7XeG1gftZ4LjABeM7PPgFOB2RpsFhFJnngmhSXAYDMbYGbZwBXA7Lqd7l7p7gXu3t/d+wOLgQvdvSSOMYmISCPilhTcvRqYDLwMfAjMcPdSM7vbzC6M13mbcMRdUGlIdW4dVOfWIe51NndvupSIiLQKuqNZREQilBRERCSiRSaFpqbXMLO2ZvZCuP9vZtY/8VE2rxjq/D0zW25m75nZQjOL6ZrlVNYap1GJpc5mdln4ty41s+cSHWNzi+Hfdl8ze9XMloX/vsclI87mYmaPmdlGM/vgEPvNzB4Kfx/vmdmJzRqAu7eoF9AG+BQYCGQD7wLD6pX5F+B34fIVwAvJjjsBdT4byAmXb2gNdQ7L5QJvEFzdVpzsuBPwdx4MLAM6h+vdkh13Auo8FbghXB4GfJbsuI+wzqOAE4EPDrF/HDCf4F6wU4G/Nef5W2JLIZbpNS4CngyXZwLnmllDN9uliybr7O6vuvvucHUxwX0j6aw1TqMSS52vA6a4+zYAd9+Y4BibWyx1dqBTuJzHgfdDpR13fwPY2kiRi4CnPLAYyDezouY6f0tMCg1Nr9HrUGU8uHS2EuiakOjiI5Y6R/sOwTeNdNZknaOnUUlkYHEUy9/5aOBoM/sfM1tsZmMSFl18xFLnO4GrzKwcmAd8NzGhJc2X/f/+pcRzmotkaXR6jS9RJp3EXB8zuwooBs6Ka0Tx12zTqKSRWP7OmQRdSKMJWoN/NbMR7r49zrHFSyx1ngA84e73m9lpwNNhnWvjH15SxPXzqyW2FJqaXuOAMmaWSdDkbKy5lupiqTNm9jXgxwR3jn+RoNjipTVOoxLrv+0X3X2/u68GVhAkiXQVS52/A8wAcPdFQDuCieNaqpj+vx+ulpgUGp1eIzQbuCZcHg+84uEITppqss5hV8ojBAkh3fuZoXVOoxLLv+1ZBBcVYGYFBN1JqxIaZfOKpc5rgHMBzGwoQVLYlNAoE2s2cHV4FdKpQKW7VzTXwVtc95G7V5tZ3fQabYDHPJxeAyhx99nANIIm5kqCFsIVyYv4yMVY5/8EOgJ/CMfU17h7sqYbOWIx1rlFibHOLwPnmdlyoAa41d23JC/qIxNjnW8Bfm9mNxN0o0xM5y95ZvY8QfdfQThOcgeQBeDuvyMYNxkHrAR2A99q1vOn8e9ORESaWUvsPhIRkcOkpCAiIhFKCiIiEqGkICIiEUoKIiISoaQgUo+Z1ZjZO2b2gZm9ZGb5zXz8iWb2m3D5TjP7fnMeX+RIKCmIHGyPu4909xEE97H8a7IDEkkUJQWRxi0iarIxM7vVzJaE89jfFbX96nDbu2b2dLjtH8LndSwzs7+YWfckxC/ypbS4O5pFmouZtSGYPmFauH4ewTxCJxNMSjbbzEYBWwjmlDrD3TebWZfwEG8Cp7q7m9m1wA8I7r4VSVlKCiIHa29m7wD9gaXAn8Pt54WvZeF6R4IkcTww0903A7h73eSKvYEXwrnus4HVCYle5Aio+0jkYHvcfSTQj+DDvG5MwYCfh+MNI919kLtPC7c3NF/Mw8Bv3P1Y4HqCidpEUpqSgsghuHslcCPwfTPLIpiU7dtm1hHAzHqZWTdgIXCZmXUNt9d1H+UB68LlaxBJA+o+EmmEuy8zs3eBK9z96XBq5kXhTLNVwFXhrJ33Aq+bWQ1B99JEgieC/cHM1hFM3T0gGXUQ+TI0S6qIiESo+0hERCKUFEREJEJJQUREIpQUREQkQklBREQilBRERCRCSUFERCL+PxmMKCKswrkNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(xPR,yPR,color='k')\n",
    "plt.xlabel(\"Recall\")\n",
    "plt.ylabel(\"Precision\")\n",
    "plt.ylim(0.4,1.01)\n",
    "plt.plot([0,1],[sum(y[9000:]) / 1000, sum(y[9000:]) / 1000], lw=0.5, color = 'grey')\n",
    "plt.title(\"Precision/Recall Curve\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = dataDir + \"beer_50000.json\"\n",
    "f = open(path)\n",
    "data = []\n",
    "for l in f:\n",
    "    data.append(eval(l))\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count occurrences of each style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "categoryCounts = defaultdict(int)\n",
    "for d in data:\n",
    "    categoryCounts[d['beer/style']] += 1\n",
    "\n",
    "categories = [c for c in categoryCounts if categoryCounts[c] > 1000]\n",
    "\n",
    "catID = dict(zip(list(categories),range(len(categories))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build one-hot encoding using common styles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feat(d):\n",
    "    feat = [0] * len(catID)\n",
    "    if d['beer/style'] in catID:\n",
    "        feat[catID[d['beer/style']]] = 1\n",
    "    return feat + [1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [feat(d) for d in data]\n",
    "y = [d['beer/ABV'] > 5 for d in data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegression()"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = sklearn.linear_model.LogisticRegression()\n",
    "mod.fit(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute and report metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "def metrics(y, ypred):\n",
    "    TP = sum([(a and b) for (a,b) in zip(y, ypred)])\n",
    "    TN = sum([(not a and not b) for (a,b) in zip(y, ypred)])\n",
    "    FP = sum([(not a and b) for (a,b) in zip(y, ypred)])\n",
    "    FN = sum([(a and not b) for (a,b) in zip(y, ypred)])\n",
    "\n",
    "    TPR = TP / (TP + FN)\n",
    "    TNR = TN / (TN + FP)\n",
    "\n",
    "    BER = 1 - 0.5*(TPR + TNR)\n",
    "    \n",
    "    print(\"TPR = \" + str(TPR))\n",
    "    print(\"TNR = \" + str(TNR))\n",
    "    print(\"BER = \" + str(BER))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TPR = 0.9943454210202828\n",
      "TNR = 0.27828418230563\n",
      "BER = 0.3636851983370436\n"
     ]
    }
   ],
   "source": [
    "ypred = mod.predict(X)\n",
    "metrics(y, ypred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Balance the classifier using the 'balanced' option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegression(class_weight='balanced')"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = sklearn.linear_model.LogisticRegression(class_weight='balanced')\n",
    "mod.fit(X,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TPR = 0.6779348494161033\n",
      "TNR = 0.9751206434316354\n",
      "BER = 0.1734722535761306\n"
     ]
    }
   ],
   "source": [
    "ypred = mod.predict(X)\n",
    "metrics(y, ypred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Precision/recall curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "probs = mod.predict_proba(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "probY = list(zip([p[1] for p in probs], [p[1] > 0.5 for p in probs], y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "probY.sort(reverse=True) # Sort data by confidence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "xPR = []\n",
    "yPR = []\n",
    "\n",
    "for i in range(1,len(probY)+1,100):\n",
    "    preds = [x[1] for x in probY[:i]]\n",
    "    labs = [x[2] for x in probY[:i]]\n",
    "    prec = sum(labs) / len(labs)\n",
    "    rec = sum(labs) / sum(y)\n",
    "    xPR.append(rec)\n",
    "    yPR.append(prec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG9NJREFUeJzt3XuUlfV97/H3hxm5CiqCxnA1CrFEjZcp6spaXmqigCnk4rFwjifiSaRNY4ymkmrsSawnxpVWTGNDavFoRVeNQc+KJZWGVuK9ahjqJQFDBKIyXiJeIJHbAPM9fzzPPG6GPbMfhv3sPTN8XmvtNc/lt5/n+5uB/dnPXRGBmZkZQL96F2BmZj2HQ8HMzDIOBTMzyzgUzMws41AwM7OMQ8HMzDIOBTMzyzgUrFeRtFLSmRXajJX0nqSGGpVVKElnSmopGX9J0sfrWZP1XQ4Fq4r0g2pr+mH8W0n/JOnAaq8nIj4SEQ9XaPNKRBwYEbuqsU5JH2z/UO7Qzzck3VFEP/eFpMmSlkjaKOkdST+XdHG967LewaFg1fTHEXEgcBLwh8BfdWygRG/7dzcN+GnJeHs/TwBOBK6uS1VlSDoN+BnwCHA0cCjwRWBqN5fXJ7a2LL/e9p/TeoGIeBX4N+BYAEkPS7pe0hPAFuBDkg6SdJuk1yW9KulbpR9Aki6R9IKk30taJemkdHq26yT9Rtws6Xfp1slN6fTxkkJSYzr+QUmL02/NayRdUrKeayUtknRnuq6Vkpo6dGkasKRMP98AlpKEQ/vyBki6UdIraU23SBpUMn+GpGfTmtdKmpJOv7ikv+sk/Wk3f/1/CyyMiO9ExFuRWBERF6TrmS3p8dI3pL+ro9PhOyT9Q7qlsRm4Ot0iKv3bfFrS8+lwP0lXpX15O/1dDu9m7dYDOBSs6iSNIfkgfaZk8v8E5gBDgZeBhcBOkm+zJwLnAF9I3//fgGuBzwHDgOnA22VW9T3gexExDDgKWNRJST8EWoAPAucD35Z0dsn86cA9wMHAYuD7JX05ADgd+I8y/RxN8g18Tcnk7wATSYLiaGAU8I20/WTgTmBuuq7TgZfS970JfDLt78XAd9uDMC9Jg4HTgPv25n1l/HfgepK/1Y3AZuCPOsy/Ox2+DPgUcAbJ7/ddYP4+rt/qKSL88mufXyQfbu8BG0k+9H8ADErnPQxcV9L2cGB7+/x02izgoXR4KfCVLtbz8XT4UeCvgREd2owHAmgExgC7gKEl828A7kiHrwUeLJk3CdhaMn42sKxMP3+frmMZcHA6TyQfoEeVtD8N+E06/I/Ad3P+Pu9v/x0AZwIt5X4HHd4zKq3pmC6WOxt4vMO0AI5Oh+8A7uww/1vA7enw0LSP49LxF4CzS9oeAewAGuv9b9Kv7r28pWDV9KmIODgixkXEn0fE1pJ560uGxwEHAK+nB0M3knxgHpbOHwOszbG+z5N8K/+VpOWSPlmmzQeBdyLi9yXTXib5AG33RsnwFmBg+64nyu86+lREDCX5sD4GGJFOHwkMBlaU9Oun6fQu+yVpqqSn0l1cG9P1jijXtgvvAm0kH8z7Yn2H8buBz0gaAHwG+K+IeDmdNw74cUl/XyAJ4cP3sQarE4eC1UrpPdrXk2wpjEhD5OCIGBYRHymZf1TFBUa8GBGzSMLkO8B9koZ0aPYaMFzS0JJpY4FXc9Y9DXigk/U/QvLN+sZ00lvAVuAjJf06KJKD0tBJv9IP2/+XLufwiDiYJIiUs8b2erYATwKf7aLZZpLgal/3B8otqsNyV5EE6VR233UESZ+mlvT34IgYGMlxJeuFHApWcxHxOvDvwDxJw9KDlUdJOiNt8n+BKyWdnJ6tdLSkcR2XI+lCSSMjoo1ktxUk31JL17Ue+E/gBkkDJR1PsoXxz5XqlHQkMCAiftVFs78DPiHphLSOW0mOBxyWLmOUpHPTtrcBF0s6O+3zKEnHAP2BAcAGYKekqSTHWLrja8BsSXMlHZrW8FFJ96TznwM+IukESQNJdp/lcTfJ8YPTgXtLpt8CXN/+95E0UtKMbtZuPYBDwerlcyQfhqtIdnvcR7rbIyLuJTnQeTfJvvv7gXJntEwBVkp6j+Sg88yI2Fam3SyS4wyvAT8GvhkRexw4LuM8ypx1VCoiNpAcPP7f6aS/JDnw/JSk3wEPAh9O2/6c9CAysInktNFx6a6ty0gOlL9L8m18cY76ytXznyQHhf8IWCfpHWBBez8i4tfAdWldLwKPd7Kojn5IsrvsZxHxVsn076W1/ruk3wNPAad0p3brGRThJ6+ZlSNpCfD9iOgyGMz6Em8pmHXuYeChehdhVkveUjAzs4y3FMzMLNNYuUnPMmLEiBg/fny9yzAz61VWrFjxVkSMrNSusFCQdDvJZftvRsSxZeaL5MyFaSQXDM2OiP+qtNzx48fT3Nxc7XLNzPo0SS9XblXs7qM7SE4Z7MxUYEL6mgP8Q4G1mJlZDoVtKUTEo5LGd9FkBsk9VoLknO6DJR2RXthUdQsXLuTGG2+s3LCKko2hvreuWq+vlus64IADOOigg2hs3PO/RmcnZZSbvjdtAc466yy+/vWv70WlZsWo5zGFUex+j5WWdNoeoSBpDsnWBGPHju3Wyg455BAmTpzYrfd2Ry3P6qr1GWR9uW/bt29n06ZNtLW1lZ3fWUCVm5637dq1a3nwwQd57LHHGDp0KG1tbUQEbW1tXQ4fdthhzJgxg+HDh2fTyr127dpVdvrw4cOZNm0aDQ1+ZIK9r9BTUtMthX/t5JjCA8ANEfF4Or4M+FpErOhqmU1NTeFjCtaXtLa2csUVV7Bs2TL69euHJPr161dxeO3atbz55pv7tO5TTjmFiy66iF27dnX6mjRpEueff36u5WV32uznExt7GkkrIqLjs0L2UM8thRaSu0a2G01yGwKz/Ur//v2ZP3/vH0Gwfft2VqxYkX0I781LEjfddBPz58/n6aefrriuQw89lCFDhrBz50527dqV/Swd3rlzJ21tbQwZMoRFixYxbdq07vw6rM7quaVwHnApydlHpwA3R8TkSsv0loJZdUQEb7zxBpJoaGgo+2pra2PevHm89NJL7Nixg8bGRhobG7P57cOl07797W9z3nnncf/999e7i1ai7lsKktpvoDVCyUPPv0lyD30i4haSG3RNI7l52BaSG4WZWY1I4ogjKj964Zprrtmr5T7wwAPs2LGju2VZnRV59tGsCvMD+FJR6zez+ujXr1/NTxCw6ul1VzSbWc8mqdOzt+x9bW1tbN++fbdXa2vrHtN27NjBqaeeytChQysvtAocCmZWVT19SyEiaG1tZevWrWzdupVt27aV/VluWscP7Eqvch/y7a+dO3fmrvmrX/0q8+bNK/C38j6HgplVVXe2FCKCrVu3smXLFjZv3szmzZuz4Y4/24e786G+bds2tm3btk+h1dDQwIABA7p89e/fn6FDh1Zs1962q/nnnnsua9fmeWR5dTgUzKyq+vXrx4svvsjcuXP36kN+b/Xv35+BAwcyaNAgBg0alA23/xw2bFjZ6R1/5mkzcODA7FXri/0mTZpU0/U5FMysqo455hieeOIJ5s+fz5AhQxg8ePBuP0eMGMHYsWN3m1auXVfTBg8e7CuxC+JQMLOquvXWW1mwYIGvau6lHApmVlWSan6DRqseR7mZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZptBQkDRF0mpJayRdVWb+OEnLJD0v6WFJo4usx8zMulZYKEhqAOYDU4FJwCxJkzo0uxG4MyKOB64DbiiqHjMzq6zILYXJwJqIWBcRrcA9wIwObSYBy9Lhh8rMNzOzGioyFEYB60vGW9JppZ4DPpsOfxoYKunQAmsyM7MuFBkKKjMtOoxfCZwh6RngDOBVYOceC5LmSGqW1Lxhw4bqV2pmZkCxodACjCkZHw28VtogIl6LiM9ExInANem0TR0XFBELIqIpIppGjhxZYMlmZvu3IkNhOTBB0pGS+gMzgcWlDSSNkNRew9XA7QXWY2ZmFRQWChGxE7gUWAq8ACyKiJWSrpM0PW12JrBa0q+Bw4Hri6rHzMwqayxy4RGxBFjSYdo3SobvA+4rsgYzM8vPVzSbmVnGoWBmZhmHgpmZZRwKZmaWcSiYmVnGoWBmZhmHgpmZZRwKZmaWcSiYmfVwmzZtorW1tSbrciiYmfVghxxyCA8//DAf+MAHeOeddwpfn0PBzKwHu/nmm5k8eTLvvvsutXh0gEPBzKwHO+6447j88ssBiOj4SJrqcyiYmfVwUrlnlhXDoWBm1kt4S8HMzLItBYeCmZl595GZmb2vsTF5Htpvf/vbwtflUDAz6+FOPPFEAFavXl34uhwKZmY93KBBg2q2LoeCmZllHApmZpZxKJiZWcahYGZmGYeCmZllHApmZpZprHcBtfTss8+ycePGepdhZrZXtm/fzty5cxk+fHjh69qvQuGEE06odwlmZt1y7rnn1mQ9he4+kjRF0mpJayRdVWb+WEkPSXpG0vOSphVZj5mZda2wUJDUAMwHpgKTgFmSJnVo9lfAoog4EZgJ/KCoeszMrLIitxQmA2siYl1EtAL3ADM6tAlgWDp8EPBagfWYmVkFRYbCKGB9yXhLOq3UtcCFklqAJcCXyy1I0hxJzZKaa/GMUjOz/VWRoVDuBuAdnxAxC7gjIkYD04C7JO1RU0QsiIimiGgaOXJkAaWamRkUGwotwJiS8dHsuXvo88AigIh4EhgIjCiwJjMz60KRobAcmCDpSEn9SQ4kL+7Q5hXgbABJf0ASCt4/ZGZWJ4WFQkTsBC4FlgIvkJxltFLSdZKmp83+ArhE0nPAD4HZUYuHkJqZWVmFXrwWEUtIDiCXTvtGyfAq4GNF1mBmZvn53kdmZpZxKJiZWcahYGZmGYeCmZllHApmZpZxKJiZWcahYGZmGYeCmZllcl+8JmkUMK70PRHxaBFFmZlZfeQKBUnfAf4EWAXsSicH4FAwM+tD8m4pfAr4cERsL7IYMzOrr7zHFNYBBxRZiJmZ1V/eLYUtwLOSlgHZ1kJEXFZIVWZmVhd5Q2Exez4LwczM+phcoRARC9MH5UxMJ62OiB3FlWVmZvWQ9+yjM4GFwEskz14eI+kin5JqZta35N19NA84JyJWA0iaSPKktJOLKszMzGov79lHB7QHAkBE/BqfjWRm1ufk3VJolnQbcFc6/j+AFcWUZGZm9ZI3FL4IfAm4jOSYwqPAD4oqyszM6iPv2UfbgZvSl5mZ9VFdhoKkRRFxgaRfkNzraDcRcXxhlZmZWc1V2lL4Svrzk0UXYmZm9dfl2UcR8Xo6+BawPiJeBgYAHwVeK7g2MzOrsbynpD4KDEyfqbAMuBi4o6iizMysPvKGgiJiC/AZ4O8j4tPApOLKMjOzesgdCpJOI7k+4YF0Wu6ntpmZWe+QNxQuB64GfhwRKyV9CHio0pskTZG0WtIaSVeVmf9dSc+mr19L2rh35ZuZWTXlvU7hEeCRkvF1JBeydUpSAzAf+ATQAiyXtDgiVpUs54qS9l8GTtyr6s3MrKoqXafwdxFxuaSfUP46heldvH0ysCYNECTdA8wgec5zObOAb+aq2szMClFpS6H9Xkc3dmPZo4D1JeMtwCnlGkoaBxwJ/KyT+XOAOQBjx47tRilmZpZHl6EQEe03vWsGtkZEG2S7hgZUWLbKLbKTtjOB+yJiVyd1LAAWADQ1NXW2DDMz20d5DzQvAwaXjA8CHqzwnhZgTMn4aDq/4G0myfMZzMysjvKGwsCIeK99JB0e3EV7gOXABElHpo/ynEmZ5zxL+jBwCPBkzlrMzKwgeUNhs6ST2kcknQxs7eoNEbETuBRYCrwALEpPZ71OUukB6lnAPRHh3UJmZnWW9wK0y4F7JbXv/jkC+JNKb4qIJcCSDtO+0WH82pw1mJlZwfJep7Bc0jHAh0kOIP8qInYUWpmZmdVcrt1HkgYDfwl8JSJ+AYyX5Ntpm5n1MXmPKfwT0Aqclo63AN8qpCIzM6ubvKFwVET8DbADICK2Uv46BDMz68XyhkKrpEGkF59JOgrYXlhVZmZWF3nPPvom8FNgjKR/Bj4GzC6qKDMzq4+KoSBJwK9IHrBzKsluo69ExFsF12ZmZjVWMRQiIiTdHxEn8/4DdszMrA/Ke0zhKUl/WGglZmZWd3mPKZwF/Jmkl4DNJLuQIiKOL6owMzOrvbyhMLXQKszMrEeo9OS1gcCfAUcDvwBuS290Z2ZmfVClYwoLgSaSQJgKzCu8IjMzq5tKu48mRcRxAJJuA35efElmZlYvlbYUsjuhereRmVnfV2lL4aOSfpcOCxiUjreffTSs0OrMzKymugyFiGioVSFmZlZ/eS9eMzOz/YBDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs0yhoSBpiqTVktZIuqqTNhdIWiVppaS7i6zHzMy6lvchO3tNUgMwH/gE0AIsl7Q4IlaVtJkAXA18LCLelXRYUfWYmVllRW4pTAbWRMS6iGgF7gFmdGhzCTA/It4FiIg3C6zHzMwqKDIURgHrS8Zb0mmlJgITJT0h6SlJUwqsx8zMKihs9xHJ7bU7ijLrnwCcCYwGHpN0bERs3G1B0hxgDsDYsWOrX6mZmQHFbim0AGNKxkcDr5Vp8y8RsSMifgOsJgmJ3UTEgohoioimkSNHFlawmdn+rshQWA5MkHSkpP7ATGBxhzb3A2cBSBpBsjtpXYE1mZlZFwoLhfTxnZcCS4EXgEURsVLSdZKmp82WAm9LWgU8BMyNiLeLqsnMzLqmiI67+Xu2pqamaG5urncZZma9iqQVEdFUqZ2vaDYzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLFBoKkqZIWi1pjaSrysyfLWmDpGfT1xeKrMfMzLrWWNSCJTUA84FPAC3AckmLI2JVh6Y/iohLi6rDzMzyK3JLYTKwJiLWRUQrcA8wo8D1mZnZPioyFEYB60vGW9JpHX1W0vOS7pM0ptyCJM2R1CypecOGDUXUamZmFBsKKjMtOoz/BBgfEccDDwILyy0oIhZERFNENI0cObLKZZqZWbsiQ6EFKP3mPxp4rbRBRLwdEdvT0VuBkwusx8zMKigyFJYDEyQdKak/MBNYXNpA0hElo9OBFwqsx8zMKijs7KOI2CnpUmAp0ADcHhErJV0HNEfEYuAySdOBncA7wOyi6jEzs8oU0XE3f8/W1NQUzc3N9S7DzKxXkbQiIpoqtfMVzWZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlnEomJlZxqFgZmaZQkNB0hRJqyWtkXRVF+3OlxSSmoqsx8zMulZYKEhqAOYDU4FJwCxJk8q0GwpcBjxdVC1mZpZPkVsKk4E1EbEuIlqBe4AZZdr9H+BvgG0F1mJmZjk0FrjsUcD6kvEW4JTSBpJOBMZExL9KurKzBUmaA8xJR9+TtLqbNY0A3urme3sr93n/4D7vH/alz+PyNCoyFFRmWmQzpX7Ad4HZlRYUEQuABftckNQcEfvVcQv3ef/gPu8fatHnIncftQBjSsZHA6+VjA8FjgUelvQScCqw2Aebzczqp8hQWA5MkHSkpP7ATGBx+8yI2BQRIyJifESMB54CpkdEc4E1mZlZFwoLhYjYCVwKLAVeABZFxEpJ10maXtR6K9jnXVC9kPu8f3Cf9w+F91kRUbmVmZntF3xFs5mZZRwKZmaW6ZOhUOn2GpIGSPpROv9pSeNrX2V15ejzVyWtkvS8pGWScp2z3JPtj7dRydNnSRekf+uVku6udY3VluPf9lhJD0l6Jv33Pa0edVaLpNslvSnpl53Ml6Sb09/H85JOqmoBEdGnXkADsBb4ENAfeA6Y1KHNnwO3pMMzgR/Vu+4a9PksYHA6/MX9oc9pu6HAoyRntzXVu+4a/J0nAM8Ah6Tjh9W77hr0eQHwxXR4EvBSvevexz6fDpwE/LKT+dOAfyO5FuxU4Olqrr8vbinkub3GDGBhOnwfcLakchfb9RYV+xwRD0XElnT0KZLrRnqz/fE2Knn6fAkwPyLeBYiIN2tcY7Xl6XMAw9Lhg9j9eqheJyIeBd7poskM4M5IPAUcLOmIaq2/L4ZCudtrjOqsTSSnzm4CDq1JdcXI0+dSnyf5ptGbVexz6W1UallYgfL8nScCEyU9IekpSVNqVl0x8vT5WuBCSS3AEuDLtSmtbvb2//teKfI2F/XS5e019qJNb5K7P5IuBJqAMwqtqHhVu41KL5Ln79xIsgvpTJKtwcckHRsRGwuurSh5+jwLuCMi5kk6Dbgr7XNb8eXVRaGfX31xS6HS7TV2ayOpkWSTs6vNtZ4uT5+R9HHgGpIrx7fXqLai7I+3Ucn7b/tfImJHRPwGWE0SEr1Vnj5/HlgEEBFPAgNJbhzXV+X6/95dfTEUury9RmoxcFE6fD7ws0iP4PRSFfuc7kr5R5JA6O37mWH/vI1Knn/b95OcVICkESS7k9bVtMrqytPnV4CzAST9AUkobKhplbW1GPhcehbSqcCmiHi9Wgvvc7uPImKnpPbbazQAt0d6ew2gOSIWA7eRbGKuIdlCmFm/ivddzj7/LXAgcG96TP2ViKjX7Ub2Wc4+9yk5+7wUOEfSKmAXMDci3q5f1fsmZ5//ArhV0hUku1Fm9+YveZJ+SLL7b0R6nOSbwAEAEXELyXGTacAaYAtwcVXX34t/d2ZmVmV9cfeRmZl1k0PBzMwyDgUzM8s4FMzMLONQMDOzjEPBrANJuyQ9K+mXkn4i6eAqL3+2pO+nw9dKurKayzfbFw4Fsz1tjYgTIuJYkutYvlTvgsxqxaFg1rUnKbnZmKS5kpan97H/65Lpn0unPSfprnTaH6fP63hG0oOSDq9D/WZ7pc9d0WxWLZIaSG6fcFs6fg7JfYQmk9yUbLGk04G3Se4p9bGIeEvS8HQRjwOnRkRI+gLwNZKrb816LIeC2Z4GSXoWGA+sAP4jnX5O+nomHT+QJCQ+CtwXEW8BRET7zRVHAz9K73XfH/hNTao32wfefWS2p60RcQIwjuTDvP2YgoAb0uMNJ0TE0RFxWzq93P1i/h74fkQcB/wpyY3azHo0h4JZJyJiE3AZcKWkA0huyva/JB0IIGmUpMOAZcAFkg5Np7fvPjoIeDUdvgizXsC7j8y6EBHPSHoOmBkRd6W3Zn4yvdPse8CF6V07rwcekbSLZPfSbJIngt0r6VWSW3cfWY8+mO0N3yXVzMwy3n1kZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWUcCmZmlvn/LcYMCX1zMPYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(xPR,yPR,color='k')\n",
    "plt.xlabel(\"Recall\")\n",
    "plt.ylabel(\"Precision\")\n",
    "plt.ylim(0.4,1.01)\n",
    "plt.plot([0,1],[sum(y) / len(y), sum(y) / len(y)], lw=0.5, color = 'grey')\n",
    "plt.title(\"Precision/Recall Curve\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataTrain = data[:25000]\n",
    "dataValid = data[25000:37500]\n",
    "dataTest = data[37500:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pipeline(reg):\n",
    "    mod = linear_model.LogisticRegression(C=reg, class_weight='balanced')\n",
    "    \n",
    "    X = [feat(d) for d in dataTrain]\n",
    "    y = [d['beer/ABV'] > 5 for d in dataTrain]\n",
    "\n",
    "    Xvalid = [feat(d) for d in dataValid]\n",
    "    yvalid = [d['beer/ABV'] > 5 for d in dataValid]\n",
    "    Xtest = [feat(d) for d in dataTest]\n",
    "    ytest = [d['beer/ABV'] > 5 for d in dataTest]\n",
    "    \n",
    "    mod.fit(X,y)\n",
    "    ypredValid = mod.predict(Xvalid)\n",
    "    ypredTest = mod.predict(Xtest)\n",
    "    \n",
    "    # validation\n",
    "    \n",
    "    TP = sum([(a and b) for (a,b) in zip(yvalid, ypredValid)])\n",
    "    TN = sum([(not a and not b) for (a,b) in zip(yvalid, ypredValid)])\n",
    "    FP = sum([(not a and b) for (a,b) in zip(yvalid, ypredValid)])\n",
    "    FN = sum([(a and not b) for (a,b) in zip(yvalid, ypredValid)])\n",
    "    \n",
    "    TPR = TP / (TP + FN)\n",
    "    TNR = TN / (TN + FP)\n",
    "    \n",
    "    BER = 1 - 0.5*(TPR + TNR)\n",
    "    \n",
    "    print(\"C = \" + str(reg) + \"; validation BER = \" + str(BER))\n",
    "    \n",
    "    # test\n",
    "\n",
    "    TP = sum([(a and b) for (a,b) in zip(ytest, ypredTest)])\n",
    "    TN = sum([(not a and not b) for (a,b) in zip(ytest, ypredTest)])\n",
    "    FP = sum([(not a and b) for (a,b) in zip(ytest, ypredTest)])\n",
    "    FN = sum([(a and not b) for (a,b) in zip(ytest, ypredTest)])\n",
    "    \n",
    "    TPR = TP / (TP + FN)\n",
    "    TNR = TN / (TN + FP)\n",
    "    \n",
    "    BER = 1 - 0.5*(TPR + TNR)\n",
    "    \n",
    "    print(\"C = \" + str(reg) + \"; test BER = \" + str(BER))\n",
    "\n",
    "    return mod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C = 1e-06; validation BER = 0.16646546520651895\n",
      "C = 1e-06; test BER = 0.2640150292967679\n",
      "C = 1e-05; validation BER = 0.28744502888678103\n",
      "C = 1e-05; test BER = 0.39631747406856366\n",
      "C = 0.0001; validation BER = 0.28744502888678103\n",
      "C = 0.0001; test BER = 0.39631747406856366\n",
      "C = 0.001; validation BER = 0.28744502888678103\n",
      "C = 0.001; test BER = 0.39631747406856366\n"
     ]
    }
   ],
   "source": [
    "for c in [0.000001, 0.00001, 0.0001, 0.001]:\n",
    "    pipeline(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the classification problem using a regular linear regressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_reg = [2.0*a - 1 for a in y] # Map data to {-1.0,1.0}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0]"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_reg[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(fit_intercept=False)"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y_reg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "yreg_pred = model.predict(X)\n",
    "yreg_pred = [a > 0 for a in yreg_pred] # Map the outputs back to binary predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TPR = 0.9943454210202828\n",
      "TNR = 0.27828418230563\n",
      "BER = 0.3636851983370436\n"
     ]
    }
   ],
   "source": [
    "metrics(y, yreg_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}