notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook was created by Michael Scherbela for the purpose of teaching the first day of the ESI-DCAFM-TACO-VDSP Summer School 2022 on \"Machine Learning for Materials Hard and Soft\". After permission, I have used and adjusted this notebook as a material for the day two of the same event (made available in repository https://github.com/paloha/MLSummerSchoolVienna2022). The aim was to smoothly follow up on the work from the previous day. Not to confuse the students with yet another dataset and model. Changelog: Pavol Harar added saving of the trained models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_digits\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.metrics import accuracy_score\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "For a first example, we will use a small dataset of 8x8 pixel images of handwritten digits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%% Load and inspect the data\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X.shape = (1797, 8, 8)\n",
      "y.shape = (1797,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8cAAAJOCAYAAACax3X+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df4zn+V0X8NeLbFTgYL7bYgkQuLkGTMDKfc9GaZF4c0SCYM2OgS1Sa3YakKqJ3jSN6cmvnZUiVyS5OaUWq3Zn/UHxSryZIEJptLNItFrIzUqtVeFulh+nSHs3Cy0UKbz9Y6aw3W7vduE185mb1+ORTOby+X7vOa/d+b4/n89zP98fOcYIAAAA6OyTph4AAAAApqYcAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwfosz8k5n59sz81cz8YGa+MzP/1NRzwRQy88sz8y2Z+b7M/FBm/mJmbmXmi6eeDabgGAEfKzPvyMz1zHwyMz+cmTuZ+RenngumkplLmTk+wddLpp7vJFKOD0lm/omI+PGI+OSI+MsHX38oIv5dZr50ytlgIn8tIhYj4uGI+OqIuD8iXhAR78rML59wLjhyjhFwU/86Is5FxIWI+KqIeHdEvDUzXzHpVDC9b46Il97w9Z5JJzqhcowx9QwnUmb+aETMI+KFY4xfO9j2aRHxeET8zzGGqwO0kpkvGGP83xu23RERPxMR7xlj/JlpJoOj5xgBHyszvzoifjgiXjHGeOt1238sIv5oRHzeGOO3ppoPppCZSxHxzog4O8b4wYnHacGV48PzpyJi+6MnPRERY4xfjf0rBV+amZ812WQwgRuL8cG2D0bEeyPic49+IpiUYwR8rL8QER+MiLfdsP1iRHx2RHzJkU8EtKMcH54/EBG/cZPtH932x45wFjiWMnMhIv54RPy3qWeBI+YYAR/rRRHx38cYH7lh+3+97nbo6o2Z+ZHM/JWD96r4sqkHOqmU48Pz3oh4SWb+zt9xZp6K3/2Xz+dPMhUcL2+MiE+NiO+cehA4Yo4R8LGeHxFP3WT7U9fdDt1ci/33anl1RNwX++/X8rkRsZ2ZXznlYCeVcnx4/kFE/JGI+N7M/JzM/NyI+L6IuPPg9t+ebDI4BjLzOyLiL0XEa8YYPzX1PHDEHCPg4z3TG+F4kxzaGWM8NsZYHWNsjjH+wxjjYkR8aUT874j47onHO5GU40MyxnhLRDwQ++9A+gsR8XMR8UUR8T0Hd/nFiUaDyWXm+Yj41oj4ljHG9049Dxw1xwj4OB+Im18dft7B95tdVYZ2xhh7EfFvIuKLM/OTp57npFGOD9EY4w0R8Rmx/9qxxTHGl0bE6Yj4UES4UkZLB8V4LSLWxhh/d+JxYDKOEfAxfjoivvDg5QXX++jr731sDfyuPPjuGRXFlONDNsb4jTHGe8YYVzPz8yLi6yLiH48xfn3q2eCoZea3xX4xfv0Y48LE48DkHCPgdzwaEXdExNfcsP1cRDwZEf/5yCeCYygzT0fEyyJiZ4zx4annOWlu/Nc5imTmi2J/B/+Tsf/uo3fH/lPo/ldEfNuEo8EkMvO1EfF3IuJHI+KHM/Ml198+xnjXJIPBBBwj4GONMX4kM98REW/KzE+PiJ+JiK+PiD8bEa/0Gcd0lJnfH/svu/nJiHh/RHxBRLw2Ij4zIlamm+zkyjFcjT8MmflHIuIfx/5HD9wR+w/sH4iIB8cYH5pyNphCZm5HxL2f6PYxRn6i2+CkcYyAj5eZd8T+pxe8PPZfa/y+iPiuMcYPTDoYTCQzH4j9ZxTdFfvHiqci4idif128e8rZTirlGAAAgPa85hgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaO/Ubdz3WH8g8s7OTmne8vJyad7Vq1dL8+6///7SvPX19dK8Q5BTD3ATx3pNbG9vl+atrKyU5lWvic3NzdK8M2fOlOYdAmviNm1sbJTmra2tleatrq6W5s1ms9K86n3AIbAmJvbyl7+8NO/pp58uzXvkkUdK806fPl2adwisidtUvZ/b3d0tzbt8+XJpXrXz58+X5lUfZ+MW1oQrxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAeznGuNX73vIdpzCfz4913tLSUmneq171qtK8hx56qDRvdXW1NC8isjqwQOma2N7eroyL++67rzTvuDt//nxp3traWmneITjxa2Jvb68yrny/Xv0YWV9fL827cuVKad7TTz9dmjebzUrzosGaqFb9O33e855Xmve6172uNK/agw8+OPUIz+bEr4lu505nzpwpzavuEwsLC6V51f0pbmFNuHIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALR3aqofvLW1VZp35cqV0ryNjY3SvNlsVpp35513luZVz8ftq37MLSwslOatrKyU5j388MOleUtLS6V5TO/y5culeYuLi6V51WtiZ2enNK/6uLi7u1uaN5/PS/M6eNvb3laad/r06dK8s2fPluY9+OCDpXkvfvGLS/Oq5+P2VZ873XvvvaV56+vrpXn2m4fPlWMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoL0cY9zqfW/5jlNYWVkpzbt06VJp3nH3zne+szRvaWmpNC8isjqwQOma2Nvbq4wrz1teXi7Nq55vd3e3NO85wJq4TYuLi6V5s9nsWOdV/3nn83lp3traWmleNFgTb3jDGyrj4umnny7Ne/zxx0vzHnnkkdK86r+/F77whaV5Z8+eLc2LBmui+nyzej+8ublZmsfv27OuCVeOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2Tk09QJWNjY1jnVdtaWmpNO/atWuledy+2WxWmre3t1ead+XKldK8ixcvluZx8hz3NXHcra2tTT0CxU6fPl2a91M/9VOleY888khpXrWzZ8+W5r3hDW8ozauer4OVlZXSvFe96lWleYuLi6V51X9ex4mP58oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANDeqakH4Pdmb2+vNG9hYaE0j+ltbm5OPcIzWllZKc3b2toqzXviiSdK81ZXV0vz4NnMZrPSvOrjDrfv7NmzpXkPPvhgaV43jz/++NQjtFd9LnH69OnSvEcffbQ078KFC6V5y8vLpXnz+bw0bwquHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7Z2a6gfv7OyU5q2trZXmra+vl+ZduXKlNG93d7c0b2lpqTSP6e3t7U09wjOaz+eledVr7OLFi6V53L7qx3D1Y257e7s0bzableZtbm4e6zxu3+nTp0vzHnjggdK8r/iKryjNe93rXlea9+Y3v7k0j5NnYWGhNG9xcbE0j8PnyjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0F6OMW71vrd8xylsb2+X5q2trZXmXb58uTTv3LlzpXkbGxuleYcgpx7gJo71mtjd3S3Nm8/npXnXrl0rzbt48WJp3srKSmneIbAmblP1cWJ1dbU0r3rNXrp0qTTvzJkzpXmHwJqY2Ktf/erSvDe/+c2lebPZrDSver6zZ8+W5oU1cduqz3WuXLlSmle9H97c3CzNew541jXhyjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0F6OMaaeAQAAACblyjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwfksxcyszxCb5eMvV8MIXMvCczNzPzycz8tcx8X2Z+e2Z+ytSzwVHLzC/PzLccrIMPZeYvZuZWZr546tlgCpm58QznTs6faCcz55n5w5n5c5n565n5VGb+p8x85dSznVSnph6ggW+OiHfesO09UwwCU8rML4qI/xgR/yMiViPi/RHxpyPi2yPixRFxZrrpYBJ/LSKeHxEPR8R7I+IPR8RrI+JdmfmVY4x/P+VwMIHviIjvu8n2H4qI34iIdx/tODC5WUT8fES8NSJ+MSI+NSL+UkT888xcHGO8fsrhTqIcY0w9w4mUmUuxX4rPjjF+cOJxYHKZ+fqI+JaI+Pwxxs9et/0fRcQ3RcTzxhhPTzUfHLXMfMEY4//esO2OiPiZiHjPGOPPTDMZHB+ZeW9EbEfE68cY3zbxOHAsZOa7IuKzxxifN/UsJ42nVQNH5TcPvl+7YfteRPx2RPy/ox0HpnVjMT7Y9sHYv4r8uUc/ERxL3xARIyLeMvUgcIy8PyI+MvUQJ5FyfPjemJkfycxfycy3Z+aXTT0QTORS7BfhN2XmCzPz0zLzZRHx6oh44xjjQ9OOB9PLzIWI+OMR8d+mngWmdrAevjYi/t0Y44mp54GpZOYnZeapzPzDmfnXI+IrI+INU891EnnN8eG5FvuvI9uOiA9ExOdHxN+KiO3M/HNjjLdPOBscuTHGbma+NCIejYifve6mvx/7r0EGIt4Y+68p+86pB4Fj4Osj4pMj4p9OPQhM7B/G/sWEiP1n2v3NMcY/mnCeE8trjo9QZs4i4qcj4qkxxt1TzwNHKTMXI+IdEfFLEfFQRPxyRHxJRHxrRPzgGOMbJhsOjoHM/I7YXw9/Y4zxvVPPA1PLzHdHxF0R8TljjN+Yeh6YSmZ+XkS84ODrz8f+e7W8bozxPZMOdgIpx0csM98UEX81Ij5ljPHrU88DRyUzfyAi7ouIF17/FOrMfFXsv5ZsaYxxear5YEqZeT4i1iLiW8YYf3ficWBymfnFEXElIh4eY3h2EVznoE98Y+y/KdcvTz3PSeI1x0cvD777Vwm6mUfEe2/y2uKPfjTHi454HjgWrivGa4ox/I6PPpvon0w6BRxP/yX2Xx77wqkHOWmU4yOUmacj4mURsTPG+PDU88ARezIi/ujBR9Vc76UH33/hiOeByWXmt8V+MX79GOPCxOPAsZCZfzAiXhkR/2WM8Z6p54Fj6L7Y/6SPx6ce5KTxhlyHJDO/PyJ+LiJ+Mvbfbv0LIuK1EfGZEbEy3WQwmfWI2IyId2TmQ7G/Ll4SEX879j+65kcmnA2OXGa+NiL+TkT8aET8cGa+5PrbxxjvmmQwmN5yRDwvXDWmucx8c0T8SuxfKf6liPiMiDgbEV8XEX/PU6rrec3xIcnMB2L/gXtXRNwREU9FxE9ExHeNMd79TP8vnFSZeV9EPBARXxwRCxHx8xHxQ7G/Lj4w5Wxw1DJzOyLu/US3jzHyE90GJ1lm/lhEfGlEfNYY41enngemcvC+LK+KiC+MiFlEfDD2X4v/T8YY/2LK2U4q5RgAAID2vOYYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGjv1G3c91h/IPLy8nJp3tbWVmletTNnzpTmbW5uluYdgpx6gJs41mui+jH80EMPleZdvny5NK96vtXV1dK8Q2BN3Ka1tbXSvAsXLpTmnT9/vjSv+s/7HGBNTGx7e7s0r3o/fOXKldK8hYWF0rydnZ3SvMXFRWviNlXvNzc2NkrzqlXPt7S0VJp3CJ51TbhyDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0d2rqAarMZrOpR3hGDz30UGne8vJyaR7T29raKs27//77S/M2NjZK8zY3N0vz1tbWSvPuuuuu0rwzZ86U5nH7dnd3px7hGV24cKE0r/o4MZ/PS/OYXvWaqH7MLS4ulubdfffdpXlXrlwpzdve3i7NW1lZKc07jqrPnap/B9XnTjs7O6V51edO1X9/U3DlGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABo79RUP3hvb68079KlS6V5999/f2ne6upqaR4nz8LCQmnebDYrzdvd3S3N29zcLM1bX18vzbv77rtL85je2tpaad7Kysqxzqs+znLyVO+Hq9dY9bnTzs5Oad4999xTmre4uFia10H1udPGxkZpXvV+uHq+CxculOadBK4cAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADt5RjjVu97y3ecwnw+L827cuVKad758+dL89bW1krzngNy6gFuonRN7O3tVcaVr4mrV6+W5j3xxBOleYuLi6V5zwEnfk0cd9X74Y2NjdK83d3d0rznAGtiYtvb21OP8Iyq53sOrFlrYmKZtb+ChYWF0rzqx9xsNivNOwTP+gtx5RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaO/U1ANU2dnZKc3b3t4uzVteXi7NW1paOtZ53L7d3d3SvL29vdK8e++9tzRvfX39WOdx8mxtbZXmXbhwoTTv/PnzpXnwbBYXF0vzrl69Wpp33F28eHHqETjmHnvssdK8jY2N0rz5fF6aV93HZrNZad6tcOUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGgvxxi3et9bviMfb3FxsTRvZWWlNG9tba007xDk1APcROma2NraqoyLhx56qDRve3u7NK/6MTyfz0vzVldXS/MOwYlfE7u7u5Vx5Y+R2WxWmlf9mKuer/rvrzovGqyJasvLy6V51ec6e3t7pXnV+5Tq4+IhsCZ4Ruvr66V51Wv2EPrJs64JV44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPZOTT1AlZ2dndK8zc3N0ryrV6+W5i0tLZXmMb0777yzNK96TWxvb5fmra2tlebN5/PSvNXV1dI8bt+VK1dK865du1aaV+01r3lNad7CwkJpXvU+gOlVn+tUq96vb2xslOZx8lSfO+3t7ZXmVateEysrK6V5U3DlGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABo79TUA1RZXl4uzbt69Wpp3v3331+at7S0VJrH9ObzeWnepUuXSvNWVlZK86qtra1NPQLFzpw5U5p37ty50rydnZ3SvOrjWPWaXVxcLM3j5Nne3i7Nqz4uVudx8lSfS2xtbZXmLSwslOZVnytWH7en4MoxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANBejjGmngEAAAAm5coxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccH5HM/MbMHJn5walngSlk5qdl5ndn5o9l5i8frIe1qeeCqWTmHZm5nplPZuaHM3MnM//i1HPBVDLznszcPFgTv5aZ78vMb8/MT5l6NphCZi4dnC/d7OslU893Ep2aeoAOMvNzIuJ7IuLJiFiYeByYyvMj4psi4kpEbEbEN047DkzuX0fEn4iIByLif0bEKyLirZn5SWOM7590MjhimflFEfEfI+J/RMRqRLw/Iv50RHx7RLw4Is5MNx1M7psj4p03bHvPFIOcdMrx0fi+iPjxiHgqIr524llgKlcj4vQYY2TmZ4RyTGOZ+dUR8RUR8YoxxlsPNr8zM++MiL+Xmf9qjPFb000IR+4VEfGHIuJrxhg/e7Dt32fmZ0XEN2Xm6THG09ONB5P6X2OMd009RAeeVn3IMvOVEXFvRPz1qWeBKY0DU88Bx8RfiIgPRsTbbth+MSI+OyK+5Mgngmn95sH3azds34uI346I/3e04wAdKceHKDNfEBHrEfHAGOMXpp4HgGPjRRHx38cYH7lh+3+97nbo5FLsF+E3ZeYLD96n4mUR8eqIeOMY40PTjgeTemNmfiQzfyUz356ZXzb1QCeVp1Ufrn8Y+6+dedPUgwBwrDw/Ih6/yfanrrsd2hhj7GbmSyPi0Yj42etu+vux/xpk6OhaRDwcEdsR8YGI+PyI+FsRsZ2Zf26M8fYJZzuRlONDkplfExF/PiLu8VRSAG7imY4Njhu0kpmLEfFDEfFLsf/+LL8c+y8v+NaIuCMivmGq2WAqY4zHIuKx6zb9h8x8NCJ+OiK+OyKU42LK8SHIzDsi4o0R8Q8i4snMnB3c9AcObp9FxG96ihBAWx+Im18dft7B96duchucZA9GxKdHxPy686Mfz8z3R8RbMvOfjTEuTzceHA9jjL3M/DcR8Vcz85PHGL8+9UwnidccH47PiIjPjIjXRsTT1319fUR86sF//8vJpgNgaj8dEV+YmTf+I/UfO/juIzroZh4R773JhYN3H3z3Onz4XXnw3bOMirlyfDj+T0Tcd5PtD8T+O1d/Vex/fh8APT0aEX8lIr4mIv7VddvPRcSTEfGfpxgKJvRkRLwoM+8YY3zwuu0vPfjujU0hIjLzdES8LCJ2xhgfnnqek0Y5PgQHD9TtG7dn5kpE/NYY4+Nugw4y86ti/9kTn3aw6Ysy86Of/f1vxxi/Ns1kcLTGGD+Sme+I/Xfm/fSI+JnYf3bRn42IV/qMYxpaj4jNiHhHZj4U+xcRXhIRfzsi3hsRPzLhbDCJzPz+iPi5iPjJ2F8TXxD7z0z9zIhYmW6ykyu9V9TRycyNiPjaMcYdU88CU8jM3Yi48xPcfNcYY/fopoFpHbw/xXdGxMtj/7XG74uI7xpj/MCkg8FEMvO+2H+W3RdHxEJE/Hzsv0nXd40xPjDlbDCFzHwgIr4uIu6K/TemeyoifiL218S7n+n/5fdGOQYAAKA9b8gFAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALR36jbu2+oDkdfW1krzNjY2SvN2d3dL854DcuoBbsKa+H3Y3NwszateY/P5vDTvEFgTt2llZaU0b3t7uzSv+jFXvWatid+T0jWxt7dXGReLi4ulebPZrDSveo1V/3mfA078mtja2qqMi3PnzpXmPQf2m6WqjztLS0uleXELa8KVYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgvVNTD1Bld3e3NO/ChQuleefOnSvNg2eztrZWmle9Jqptbm6W5s3n89I8bt/29nZp3qVLl0rz7r777tK85eXlY523s7NTmjebzUrzOqjez1X/DqrPxaod98fw4uJiaV4Hj89Qc/gAAAZsSURBVD76aGlet/1S9blO9RpbWloqzbsVrhwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO3lGONW73vLd5zC0tJSad7ly5dL855++unSvN3d3dK8q1evluadOXOmNC8isjqwQOmaOO6P4ePuiSeeKM1bXFwszTsEJ35NbG9vV8bFfffdV5p33B9zy8vLpXnV+6jV1dXSvGiwJqqP/XfddVdp3mOPPVaaN5vNSvOq18TGxkZp3nw+L82LBmuier9U/Zjb3NwszeP37VnXhCvHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7p6b6wVtbW6V5ly9fLs27ePFiad7e3l5p3j333FOad//995fmnTlzpjSvg/l8fqzzHn744dK8c+fOleYtLi6W5jG97e3tqUd4Rsf9MTebzUrz7rrrrtI8bl/1Y676WL28vFyaV33utLKyUppXfZzl9r3mNa8pzas+N6l+zG1sbJTm8fFcOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2js11Q++ePHiVD/6lqysrJTmra2tleZVW1xcnHqE9tbX10vzdnd3S/M2NjZK8zzm4HBtb2+X5lUfF5ne0tJSad7W1lZp3pkzZ0rzqo+zTK/6MVK937znnntK81ZXV0vz5vN5ad5J4MoxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANDeqal+8Gw2m+pH35LMnHqEZ3TmzJnSvJWVldI8pre7u1uad+3atdK8474PgKO2t7dXmnf16tXSvPl8XprH9NbW1krz7r777tK8nZ2d0jx4NtX7uYWFhdK86uMEH8+VYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgvVNT/eD19fXSvPl8Xpq3u7tbmvfwww+X5q2urpbmzWaz0jymd9x/pzs7O1OPwDG3tLRUmnfhwoXSvL29vdK8lZWV0ryFhYXSvOO+T+mg+tykWvW53X333VeaV33cqT735PZV74c3NjZK865du1aax+Fz5RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaE85BgAAoD3lGAAAgPaUYwAAANpTjgEAAGhPOQYAAKA95RgAAID2lGMAAADaU44BAABoTzkGAACgPeUYAACA9pRjAAAA2lOOAQAAaC/HGLd631u+40mwsrJSmre3t1eat7m5WZr3HJBTD3ATrdbE+vp6ad7GxkZpXvWaWFxcLM07BNbEbVpeXi7N29raKs279957S/Oq1+x8Pi/NOwQnfk1Un0tU7+euXbtWmle9JqqPE7PZrDTvEJz4NbG9vV0ZV36cqN5vVv95G3rWNeHKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7SnHAAAAtKccAwAA0J5yDAAAQHvKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQXo4xpp4BAAAAJuXKMQAAAO0pxwAAALSnHAMAANCecgwAAEB7yjEAAADtKccAAAC0pxwDAADQnnIMAABAe8oxAAAA7f1/xRomnFvCwrcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1000x600 with 15 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = load_digits()\n",
    "X = data.images\n",
    "y = data.target\n",
    "print(f\"X.shape = {X.shape}\")\n",
    "print(f\"y.shape = {y.shape}\")\n",
    "\n",
    "n_rows, n_cols = 3, 5\n",
    "fig, axes = plt.subplots(n_rows, n_cols, figsize=(2*n_cols,2*n_rows), dpi=100)\n",
    "for row in range(n_rows):\n",
    "    for col in range(n_cols):\n",
    "        ind = np.random.randint(len(data.images))\n",
    "        axes[row][col].imshow(data.images[ind, ...], cmap='binary')\n",
    "        axes[row][col].axis(\"off\")\n",
    "        axes[row][col].set_title(str(data.target[ind]))\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "To prepare the dataset, we flatten each 8x8 image into a 64-element vector.\n",
    "To avoid overfitting on the training set, we split the dataset in 2:\n",
    " - A training set, on which we will train our model\n",
    " - A test set, which we will use to test accuracy of our model on unseen data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%% Prepare the data\n"
    }
   },
   "outputs": [],
   "source": [
    "X = X.reshape([-1, 64])\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=1400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "All models in sckit-learn share a common interface:\n",
    "\n",
    "- A **fit()** function to train the model\n",
    "- A **predict()** function to evaluate the model on (new) data\n",
    "\n",
    "This makes it easy to quickly test different methods on the same dataset.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "As a first model we will try a Naive Bayesian Classifier:\n",
    "It uses Bayes' Theorem and makes the following assumptions:\n",
    " - The input features are independent of each other\n",
    " - The individual likelihoods are Gaussian\n",
    "\n",
    "$$p(y|\\boldsymbol{x}) = \\frac{p(y) p(\\boldsymbol{x} | y)}{p(\\boldsymbol{x})} \\qquad $$\n",
    "\n",
    "$$p(\\boldsymbol{x} | y) \\stackrel{\\text{indep.}}{=} \\prod_i p(x_i | y) \\stackrel{\\text{Gauss.}}{=} \\prod_i e^{-\\frac{1}{2} \\left(\\frac{x_i-\\mu_{y,i}}{\\sigma_{y,i}} \\right)^2 }$$\n",
    "\n",
    "$$p(y|\\boldsymbol{x}) \\propto p_y\\, e^{-\\frac{1}{2}\\sum_i \\left(\\frac{x_i-\\mu_{y,i}}{\\sigma_{y,i}} \\right)^2 }$$\n",
    "\n",
    "Effectively this means for each class we get 2 matrices:\n",
    "\n",
    "- $\\mu_i$ containing the average pixel value of all training samples in that class\n",
    "- $\\sigma_i$ containing the standard devaiation of pixel values within that class\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Naive Bayes classifier on raw data:\n",
      "Training accuracy: 87%\n",
      "Test accuracy    : 86%\n"
     ]
    }
   ],
   "source": [
    "print(\"Naive Bayes classifier on raw data:\")\n",
    "model = GaussianNB()\n",
    "model.fit(X_train, y_train)\n",
    "\n",
    "def evaluate_model(trained_model):\n",
    "    y_pred_train = trained_model.predict(X_train)\n",
    "    y_pred_test = trained_model.predict(X_test)\n",
    "    acc_train = accuracy_score(y_train, y_pred_train)\n",
    "    acc_test = accuracy_score(y_test, y_pred_test)\n",
    "    print(f\"Training accuracy: {acc_train:.0%}\")\n",
    "    print(f\"Test accuracy    : {acc_test:.0%}\")\n",
    "    return y_pred_train, y_pred_test\n",
    "\n",
    "_,_ = evaluate_model(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Assuming that the individual pixel values are independent is apprently too naive.\n",
    "Let's instead calculate the PCA first and then apply naive Bayes on these transformed features.\n",
    "\n",
    "Since all sklearn models share the same interface, it is easy to chain them together into Pipelines, applying one data processing step after the other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PCA + Naive Bayes classifier:\n",
      "Training accuracy: 92%\n",
      "Test accuracy    : 91%\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print(\"PCA + Naive Bayes classifier:\")\n",
    "model = Pipeline([('pca', PCA(n_components=10)),\n",
    "                  ('naive_bayes', GaussianNB())\n",
    "                  ])\n",
    "model.fit(X_train, y_train)\n",
    "_,_ = evaluate_model(model)\n",
    "\n",
    "# Save\n",
    "pickle.dump(model, open('pretrained_nb.model', 'wb'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Lastly we can try a full linear model (logistic regression):\n",
    "\n",
    "$$ p(y | \\boldsymbol{x}) = \\frac{ \\exp(\\boldsymbol{w_y}^T \\boldsymbol{x})}{\\sum_{y'}\\exp(\\boldsymbol{w_{y'}}^T \\boldsymbol{x})}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Linear classifier on raw input data:\n",
      "Training accuracy: 98%\n",
      "Test accuracy    : 98%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9IAAAGMCAYAAAA2iJpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debhlV10n/O8vhLEDKYF2BFLg1PaLpNJOqA1VYZJJq1rpdgDNVelXWltSyqjNUAEVbWlyy7YdXkQqiq39aNsVZRCVNzeINm8zVQBFUaGCYU6gIgEN03r/2PumTq6Vyl1V6+ZW3fp8nuc895699/2dde456+z9PXvvtau1FgAAAGB9ztrsBgAAAMDpRJAGAACADoI0AAAAdBCkAQAAoIMgDQAAAB0EaQAAAOggSAMAAEAHQRoAAAA6CNIAAADQQZAGYENU1eOqqlXVtx9j3lXzvG86xry/rao33zat3Jrm/+2+hftL87TtnXUevVhnzbzDVXXgJJoJAKctQRqAjbKSpCW5cHFiVd09yVcm+fgx5t0ryf2SXHHbNPGM8YokX5/k/Z1/9+gkz72Fef8myfNPplEAcLo6e7MbAMDW1Fq7tqrenmTXmlk7k3w6yUuyJkgv3L/FIF1VleROrbV/GNTUU0JV3S7J2a21G0fXbq19OMmHB9d8y8h6AHA6sUcagI10RZIvr6ovWJi2K8kbkrwyyVdV1V3XzPtMkj9ZnTAfkvzzVfWkqnpHkhuTXDTPu3tV/UJVvbeqPllV76qqn6yqOy42YqHGd1fVO6rqE/Ph5Y9d2+Cq2l1Vb62qG+d6F1fVvqpqt/Zkq2qlqt5eVQ+qqtdX1T/MbXv+HJRXl9s+t+npVfWsqnr3/LwunOffrapeWFXvnp/Xe6tquar+2ZrHu1tVvbiqrquqG6rqD6rqy47RrmMe2l1Vj6yq11TV9fP/5B1V9WPzvANJfmjh/9cWaxzr0O6quk9VvayqPjT//95RVU+pqrMWlll97k+tqh+dn+MNVfW/q+qBa+rdr6p+q6reN9f74NzeHbf2WgDARrJHGoCNdEWSJ2cKyL85T7swycuT/GmmQ78flClUr857c2vt+jV19szLPS/JB5J8qKruNNf/4kyHH791XubHkuxI8pg1NR6T5GuSPCfJDUmenuR/VdWXt9belUzBMsnvJnltkm/PtJ58apLP63jOn5/kt5L89PxYj0nyrCSfk+Q/rln2yUneOT/G3yf566q6S5Irk9wryU/Nz+v/mp/7V1bVw1prbd4zfzDJN8zz3pDkG5O8aj2NrKrvT/Li+bGelORDSb4syf3nRZ6f5J8leVymw8JXHfPw8Kr650n+LMkdkjw7yeEkj03ywkyv0Q+u+ZMfSvKXSfYuPN4rq+q+C6//K5PcLtNr9Z4k95yf77b1PEcA2DCttS13y7Rhtp7brs66d03yPzJt9Hws04bYnyf58SR3OYn2XrOmXTckeX2Sx59EzbOS/ECSN2faOLs20/mKj9zs18fttr9tVJ+Ya1eSi5P8VaY9au/OtBF99km0d3ifWFP/rBwNccub/fps5Vum8PiZJL88379Hks8m+ab5/v+X5Gfn3+89vyY/s6ZGS3Ikyeesmf4D87x/u2b60+fpD19T4wNJ7rpwf119Isk582doW8fz/T/z3/3t3ObrkrwuyR/O/4f7zMttn5f7myS3X1PjmfOyX70w7Vh94h3z709e8/c/Pk/ftzBtaZ62fb5/zlzjA0neNPfdluRea2r9/C0970xB+cDC/RfMNb52zXK/ML/mX7bmub81ye0Wlvuaefp3LLxXWpKLN/t9fCbcevpEZ93TadvpZbfwnN++2a+P221/28A+8cRbqffUE2yvPHEb37bqHumvX3P/2Zn2cjxkzfS/6Kx7h0xvzP+SKSy0TOf6XZLkwUke2d3So16b5Bnz7/dO8rQkL6uqu7TWXnwC9X4i016Z/zbXvXOmPR+vqqrdrbXfO4m2cvrZqD6RTHvcnpNpI/qPk3xdpr1jX5B/ugeqx+g+sejJmTbm2WCttY9W1VU5ep70zkwh8U/n+1fm6PvweOdH/7+ttY+umfaQTAOW/c6a6QeS/EyShyb5o4XpV7TWPjb/vtonXj634U7z499prveief5ftNZuqKrfzxRGb83dM4XGX0nyxkx7tL8ryRPm+Q/OtLG+6vdaa59aU+OxSd6e5FBVLa6nX5dpj/NvJPm9JPvn6WtP0/rvSX7yVtr5DZn2Nn8iydWZAs7OW/mbW/OQTP+v/7Nm+oEk/2Ge/86F6a9orX1m4f5b55/nzT8/kukLiafNh8VfkeSq1tpnT7KdHJttp8kNSR6+ZtonTrAWp7eN6hOXZ/qMX+sn5/oHO+stkiduQ1sySLfWXr94v6o+nOSza6efQN3rknzHmsl/VFV3TvKjVXWf1tp7TrD8Rxfa9/qqek2mjZsfzXToXa/vS7LSWrvpMMKq+qNMh+5dlGkjjDPERvWJ+VDOH0vyS621Z82TV+bzU/dV1XJr7Z23XOG4RveJ1TbfL9PK6gmZDuFl412R6TPyCzNtJLyptXbDPO/KJE+pqnPneZ/OFBjXOtbhxPdI8oE2f22+qrX2oar69Dx/0XULy7w+Sarqhkwh8mOZNvArU1hb2zc+eKvP8uhyrbX2goVpr5rPK/7XSf75muWP9bw+L8mXJFkbsL9o/vnJ1tpvV9W3ZHof/0CS5YXlPrCOdq6249taa39SVc/MyQfpe2TaS73W+xbmL7pu8U5r7cbpaPXceb7fquqhmb6oe3qmIPaRqvqNJP9p4UsRBrDtdJPPnOxzZmvYwD7xTwZ/rGmskAcmubK19jcnUV6euA0ZbGyM1c7w6VEFW2sfSfLXOfrNfK9PJVl7juE/Zjp87x9Pommw6NFJ7pjkpWumvzRTINkz6oEG9InV0Z5/JdOlgH5/UNO4dat7mHfNtysX5q2G5gfP896wELIXtWNMuy7J582v602q6nMzfVF8bWc7Vw+HO9b50J+/zhqfuoW/f9f8c+3n77Ge17VJ3pbpUOfV24cy/d++JkcvObUaUNf2ifW0dXW9da91LLte12U6EmWtL5x/9r4eaa1d3Vr7/tba5yf58iSXZjrS5WdPuJWcKk7FbSfYLN+Z5C6ZtlGGkSc2liCdpKpeN++9WO/yVVVnz6OlPirJjyT59dba+27tbzse4w6ZDsn48MK0J84jnT7hlv/yJvuTPGYeqXXbvCdoOVMn/a+j2snW1NEn7p/5/LHFia21v8t0fuj9j/VHJ9imk+0TybTnbkeSHx7VLtbltZkO535cpkGzVlZntGlQqUOZvtnenr7rR78m07m+a7+w+Z6F+b3emGTP/H5LklTVOUke31HjrvPe4kWPmH++cu3Cx/DyTINzXddae2Nr7Y2ZNmaOzPcPz8utjmx+08ZMVT0x0+Gzt+bPMm0cPWntFxFr3DjXvfM6ar4myb+sqn+1Zvr3ZPqc6Hlt/4nW2jtbaz+R6UuGtY/BbWyLbjslyTnz6PCfqaq/q6r/WlWfM6qNbF29fWKN78+03fQ/BzZJnthgW/LQ7hPw6fR9I/r4JL++cP9XMo14ejJq4Vy4e2Xa43DPTOedrvpspo3RWz0/rLX2wqr6RJJfztG9hdcleaxDlliH9faJeyT5RGvtWN9KfiT/9FDOHkP7RFXdO8l/TrK3tfbBNeeesoFaa39fVW/OFHg/m6PnR6+6MkdHbu4JW7+WaeTny+ZDp9+W6fDpH0/yytbaH59Ac5+T6YiFV1fV/kwjRj8t0/tsve+Z65L8YlXdJ9M5wc/ItJf4ta21q9fx98tJvi3Ja6vq0kznDt8xyXlV9duZ9sq+L9Oo4klyt5ouWfXGJN+8ngbO530/JdP6649z9ND176yq8xYO43vb/PMZVfWqTP+Ht7bWPnmMspdmCs2vqKrnZDqc8DGZ9iD/Yu9pHlX1gEyDnf12pj0qn8x0buIDMo2IzubacttOSd6SqR+9PdNRVRdm+my6sKq+trXmXGmOp7dPJEmq6v5JvjbJf7uF7anOcvLEbWazRzu7LW6ZBjq5YWC9uyf56kwr9GdlGsXud5PUCdZbO8peyzSAznLWjObaUfOJmfZS/Oe5nY/KNGrmDVkYydbtzLyN6hNJfjXJx25h3t8mefkJ1t2IPvGqTANWrd4/O0btvs1umQb/apkO3V47b/c878YcYxTfed7P30Lduyf5xUzB8lOZztH9qSR3XE+NefkDi30iU+B/69yeqzMF4f1JPrKO57mSaSN8Z6bLUX1yfuw/z8JI9jk6cvUxR2fNNBDY8zNdHurGTBs8x+oTvzD3w4/O9/8w0yHQLccZtXth+qPmNq+O2v1XSZ6+MP8Omc6r+9BCG7Yv/u/W1LtPpsHQrp2f+19murTXWet57ovtTvK5mTbc3pFp3fWxJFdlCja3O97r4Dakz97UJwbVO+W3nW7hcb59rv3Dm/2auG3ubXSfWKh76fwe23GSdeSJ2/hW8z9pS6uqA0ke11o7Z4PqPz7TKKzf0lrrPu+yqq7JdAmUp+bom/5d7djf+K+n3j2T/F2my83sXTPvdUk+r7X2pSdSm61hVJ+oqp9N8pQkd26t3bhm3keT/H5r7XuO+cfHrzu6T3xHpr2X3zDXTaY9jdcm+aVMA6bd0Fobdq4ep5fj9Ymqun2mw8/f21p7xNr5a5ZdSXLP1tr9q+rRmYLCq+faawcP62nf0D5xjPrPzLTH4t6ttWtG1OT0dqZtOx3ncW6XKTT8bmut5xQPtpiN6BPzodfvTXJ1a+2rT7KWPHEbc2jjGKuX+viyk6hxpE3nwY3wLzJdwuUNx5j3xiQXV9Wd2skfPgJvy3T42/0zXYs2SVJV90qyLce+vMN6jewT909y+xy7Tzxpvn1zpnNTOcNV1UsyXTbr/ZkOyX5Skq/IdL309dZYDdGvyXSd6xMO0QtG9gnYbKfattOtcdk1NsLuTIdeP3tQPXniNiRIj7F67dOTGa5+pNWBOx6Y6RC7JDeNWPx1Sa49k9/0DPWqTIeELmUhSOfoYaQncy3EkV6S6TzQRbebp/3PTOdhvm3tH3HGumuSF2a6RNSnkrw5yaPb+s+5PidTiF5J8q2j95DBFnGqbTvdkm/PFCacD8pG+P5M1yn/zc1uyDHIE7fCqN2ZDsWrqlt9I1TVD1XVr1fVd1fVhVX12PnQ1v2ZRk99+cKyXzKPiDdsGPuq+r6q+nRVfdfxlmutvSvTdd1+sKpeVFWPqKpvzrRh98BM1+KEW7TePtGmayG+INN77XlVtbOqnp7pm9VfaguDC21yn3h3a21l8Zajl2C6Zp523XFKcAZprf271tq9Wmt3bK2d01p7cJJnrqdPJHlupstfvTfTgFgXVNUDF253XV1wM/vEvOznVtXjqupxOTrC/qPnaQ8a1Sa2pq227VRVX1xVf1JV/7GqHllVj6qqn8n0RexbM41FALdovX1iYfl7J3l4kt9p0xUsjrWMPHEKs0d6cnbW97+4KtNJ9j+d6TCMT2cakfX5SS5trX1mYdnV8yfeP7CdZ2Xai7aeL0C+I9Mlfp6Q5N9nGvTlr5J8Z2vttwa2ia1pvX0iSZ6X6TI6P5hpUKb3J/nJTAM+LdrsPgEnY7194uGZ9l7dL8cegfxBOXrt7M3uEw/INCL2ol+ef74mycMGtoutZ6ttO12fadyMp2b6MuysTAPqLSd5QWvtHwa2ia2pZ9spSb430/vseCF5s9cT8sRxnBGDjW2GqnpypjBx39batZvdHths+gTcnD4BN6dPwM3pE6c2e3E2zoWZvmn1poeJPgE3p0/AzekTcHP6xCnMHmkAAADoYI80AAAAdBCkAQAAoIMgDQAAAB0EaQAAAOggSAMAAECHnouGn9LDex86dGhovX379g2tN7p927dvH1pveXl5aL0dO3YMrZekRhcc4JTuE6MtLS0Nrbdt27ah9Ua/h08D+sQWs2fPnqH1VlZWhtY7ePDg0Hq7du0aWi/6RLfR75HR205HjhwZWm/08x29HtsAW75PjN6+Hr2tc9VVVw2td/755w+tN/r/dxoY2ifskQYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0qNbaepdd94Lrcfjw4ZHlsmfPnqH1duzYMbTe0tLS0HoHDhwYWm/067GysjK0XpIaXXCAoX1itEOHDg2tt2vXrqH1tm/fPrTeJZdcMrTe7t27h9bbAPrEJhvdxy644IKh9UbbuXPn0HrWE5tv9Gtw4YUXDq032sGDB4fWs544IUP7xN69e0eWG94nRueT0ds6o/vE6PXEtm3bhtbL4D5hjzQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHc7e7AaMsn379qH19u3bN7Te6PYdPHhwaL0dO3YMrcfmW1paGlpv9Ht4ZWVlaL1du3YNrbd79+6h9dh6Rn8On3feeUPrXX311UPrsfWMfg+PdtFFFw2t53N961leXh5ab/T2/4EDB4bWG23Pnj1D643+TDnV+6w90gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0OHuzHnj79u1D6x08eHBovUOHDg2tt7S0NLTeaMvLy5vdBAa76qqrhta79NJLh9bbtm3b0Hqjny/cmgsuuGBovcOHDw+td9lllw2td6qvx84E+/btG1pv//79Q+uNduTIkaH1VlZWhtbbtWvX0HpsPVdfffXQeueee+7QeqPzzuh8d6qzRxoAAAA6CNIAAADQQZAGAACADoI0AAAAdBCkAQAAoIMgDQAAAB0EaQAAAOggSAMAAEAHQRoAAAA6CNIAAADQQZAGAACADoI0AAAAdBCkAQAAoIMgDQAAAB0EaQAAAOggSAMAAEAHQRoAAAA6CNIAAADQQZAGAACADtVaW++y615wMywvLw+t9yM/8iND64327ne/e2i97du3D623AWqzG3AMp3SfqBr7L9u9e/fQevv27Rtab9euXUPrHTlyZGi9DbDl+8To1+DAgQND662srAytd/nllw+td/755w+td+jQoaH1NsCW7xOjjX4Pj653ySWXDK032qWXXjq03t69e4fWiz7R7fDhw0Pr3fe+9x1a7zR4z53qhvYJe6QBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6HD2ZjdglF27dg2t11obWm/79u1D6x05cmRoPbae8847b2i90X3sggsuGFpv586dQ+ux+fbt2ze03v79+4fWG210n11ZWRlaj61n9Of66HpLS0tD641u3/Ly8tB6e/fuHVqPfgcPHtzsJhzXjh07NrsJLLBHGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOZ292A0bZsWPH0HpHjhwZWm/btm2ndD22nj179gytd+DAgaH1zj333KH1lpaWhtZj841+TUd/bo7uE1dfffXQevv27Rtab3l5eWg9tp6VlZWh9Uavx66//vqh9Xbv3j20Hptv9OfmxRdfPLTerl27htbj5NgjDQAAAB0EaQAAAOggSAMAAEAHQRoAAAA6CNIAAADQQZAGAACADoI0AAAAdBCkAQAAoIMgDQAAAB0EaQAAAOggSAMAAEAHQRoAAAA6CNIAAADQQZAGAACADoI0AAAAdBCkAQAAoIMgDQAAAB0EaQAAAOggSAMAAECHaq2td9l1L7gehw4dGlkuKysrQ+stLy8Prbd///6h9Xbv3j203mmgNrsBxzC0T5zqduzYcUrXO3DgwNB6pwF9YpMdOXJkaL2lpaWh9S6//PKh9Xbu3Dm03uj1dvSJbqPfI3v27Bla77zzzhtab3Qf27dv39B6G2DL94nDhw+PLDd822T059zo9p2BhvYJe6QBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6FCttc1uAwAAAJw27JEGAACADoI0AAAAdBCkAQAAoIMgDQAAAB0EaQAAAOggSAMAAEAHQRoAAAA6CNIAAADQQZAGAACADoI0AAAAdBCkAQAAoIMgDQAAAB22ZJCuqrbO264TqP0FVfULVfWuqvqHqjpcVS+uqnufRHuvWdOuG6rq9VX1+JOo+eCqeklVvamqbpzr3utE63F628g+MdffXlUHqur98/vtmqr6nROsdfYx2nV9VV1RVY88kZpz3d1V9etV9faq+nRVffpEa3H62+D1xBdW1a9V1Yer6h+r6qqqWjqJtm5Un3jdLTznl59oTU5f1hNJVd2xqp5RVX9eVR+vqg9U1Suq6oEnWhPYus7e7AZskK9fc//ZSS5M8pA10/+ip2hV3TnJ65Kck+S5Sf4yyb9IckmSR1TVv2ytffyEWpy8Nskz5t/vneRpSV5WVXdprb34BOo9PNPzfUuSjyXZeYLtYmvYkD6RJFX1gCRXJPnrJD+a5L1JvjDJw/qbeTP/I8lypi/8vjhTm19RVY9urb36BOp9a5KvzdQnPpnkASfZPk5vG7We+Jwkf5rpffvUJB9M8l1JXlpVd2ut/dyJNTfJ+D6RTP32e9ZM++iJN5HTmPVE8tIk/y7JTyVZSXKPJD+W5LVV9fWttTedZHuBLaRaa5vdhg1XVQeSPK61ds5J1nlkklclWWqtXbYw/buT/FqSb2mt/f4J1L0myRtba3sWpt09ydVJrmmtfcUJ1DyrtfbZ+fdnJnlBknu31q7prcXWM7BPVJK3JvnHJN/YWvvkgLadneRTSfa31vYuTP+yJH+V5A9aa486gbqLfeKXkjyxtbZVv0yk08A+8awkz0tyQWvtqoXpr8n0Rc4XttY+1llzo/rE65Kc01rb0fu3bH1n2nqiqu6SacfDr7XWvndh+r2TvCfJi1prTznZtgNbx5Y8tHsDfWr+ef2a6Ufmn/846oFaax/J9M3teSf4958d1RY4jockuX+S5REbR8fTWntnpj1l+gSnsm9M8r7FED17eaajmb5p1AOdbJ+A28jpsp74TJKWf7qNd/08fdg2HrA1CNK56Tyx9Zwv+dpMh4U+r6q+qqrOqaqvSvITSd6Q6bClUW26Q6ZDvD+8MO2J83lATxj1OHAsHX3iwfPPj1fVH8zng95QVb8/7xkY2aZ7JNmWm/eJn5j7xL8e+ViwVkefuEOSG48xfXXasFMKBvWJL6+qj9Q0bsDfVNXzq+pOo9rI1rXV1hOttRuT/FKS76uqb6mqu1XVfZP8P0k+kuRXRrYVOP0J0pNPz7fjaq19KtP5Qu9J8sZMhwC9MdMH9iNaayczeFHNg2ecXVXbk7wkyT2T/MbCMp/N9I2pPWtstHX1iSRfNP+8LNOpCI9J8h+SnJ/kdVX1eSfRhtU+cfuq+ook/z1J5eZ94jM5uhcBNtJ6+8RfJDmvqr5ozfTVjfh7nEQbRveJK5PszTR+wO4kr850Pugr5sNx4Xi24nrih5P8XJKDmfZEvyvJ1yV5SGvt3SfRTmALco50X507Jvm9JF+e5PlJ3plpcIv/lOnw7oe21v7+BOpek6MrmlWfSPLiJE+bA/zJtNs50tzMwD7xq0m+N8krWmuPXZj+VZm+ZLqktbavs+bquW9rHUnywtbaT554i296DOdIczMD+8T9k7wpyeszhYUPJnl8kp/NtLf651trP9xZc8P7xMJjPSPJT+cEx/xg6zgT1xNVtS/TIIEvyDS47LZM4foBSR5+jFM2gDOYjcg+/z7JIzINInNonvYnVfWnmUL1Dyc50Q2aKzN9eLckH0/yro0+lwgGuG7+ebPRUVtrb6qqDyf5VydR+zeTvChTn/j7TH3iMydRDzZca+3tVfVtSX4xyZ/Pk9+T6fP95zKNVnyibos+8bJMQfqBSQRpRjgt1hNV9ZWZrsjyI6215YXpr8p0lZYXZroiCkASQbrXjiSfWgjRSZLW2l9X1fWZBtM4UUdaa288qdbBbe+tx5lXObnTED6kT3A6aq29vKpekeRLM51C9c5Me6WTaayNE3Vb9InVQ7qdQsQop8t64vz55xsWJ7bWPllVb810iDfATZwj3ed9SW4/H450k/m8nHOTOGyaM80rk/xDkptdZqSqvibTOf6v34xGwWZrk3e21v4yye2TPDnTZQ7/bJObdmtWrymt7zLK6bKeeN/884GLE+fB9y6IbTxgDUE6SVWtVNV6Lmvwq5kOHfpfVfUDVbWrqp6YaSVxQ6aRHVdrfsk8SuSwUR6r6vvmkVW/ax3Lfm5VPa6qHpeje8ofPU970Kg2sTWtt0+01q5LckmSR1XVS6rqm6pqKcnvJjmcaQTU1ZoPm/vEswa285K5T3zjOpbdvtAn7jdNmu6v/XIM1lpvn6iqs6rq56rqW+d1xPdn2sN1Xo6G1NVlN61PVNWFVfWqqvq/53Z88zx2wPOS/GGm9Rrcoi24nrgyyZuTPL+qnltVD51P03h1kvsk+S+j2gRsDQ7tngkXnQwAAAikSURBVJyddfwvWmuH529Qn51pZNPPT/KBJP8702AZf72w+OrgHO8f2M6zktwu6/sC5AFJfnvNtF+ef74mycMGtoutZ119Iklaaz9TVR/NtMftCZm+bHpVkme21o4sLLoRfeJ28209Iww/LNMAfotW+8hLkjxxYLvYetbdJ5JsT/Lzmfa2XZvkD5Lsa629Z81ym9kn3pfpvNLnzO1smQ5Bf1aSF7UzYSRSTtaWWk+01j5TVQ9N8vQk3zH/vCHJ25M8srX26uP9PXDmOSNG7d4MVfXkTAOP3be1du1mtwc2W1W9KMm3JflSA+mBPgFr6RPA6cSh3RvnwiSXCtFwkwszHblh4wgm+gTcnD4BnDbskQYAAIAO9kgDAABAB0EaAAAAOgjSAAAA0EGQBgAAgA4915E+o0Ylu/zyy4fWu+iii4bW2759+9B6Bw8eHFpvdPuyvusE39bOqD6xZ8+eofVG97HnPve5Q+vt27dvaL0NoE90Wl5eHlpv9Hvk+uuvH1pv586dQ+utrKwMrbcB9IlNdqqvJ0b3idHbTtu2bRtaL6dmn4Atwx5pAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADqcvdkNOFVdeumlQ+vt27dvaL2VlZVTut7S0tLQevQ7cODA0HqXX3750Hrnn3/+0Hp79uwZWo/Nt2vXrqH1rrzyyqH1du/ePbTe6D525MiRofXYeka/50bXu+iii4bWu+yyy4bWO3To0NB6oz/zgI1ljzQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHc7e7Aacqnbs2DG03sGDB4fWO3To0NB6Bw4cGFqPzXf48OHNbsJx7dmzZ2i90X2WzTf6c260yy+/fGi98847b2i9Sy65ZGg9Nt/oz/WLLrrolK43ettk+/btQ+tZ78CZzR5pAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADqcvdkNGOXIkSND623btm1ovSuvvHJovSuuuGJovdHPl8138ODBzW7CcS0tLW12EzjF7d27d2i9Sy65ZGi9nTt3Dq23srIytB5bz1VXXTW03vXXXz+03uhtidHbdvv27RtaDziz2SMNAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQIezN+uBDx06NLTerl27htbbu3fv0Hrnn3/+0Hqjny9bz+HDhze7Ccd11VVXDa23ffv2ofXYfPv27Rtab2lp6ZSuN3q9s7y8PLQem+8tb3nLZjfhuPbv37/ZTTgufQIYyR5pAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADpUa229y657wfVYXl4eWS4HDx4cWm9lZWVovaWlpaH1duzYMbTe3r17h9bbALXZDTiGoX3i8OHDI8vlvve979B6p7rRnwG7d+8eWm8DbPk+caob/bk5er1z6NChofVOA1u+T5zq64mLLrpoaL3LLrtsaL2XvvSlQ+uN3rbbAKdin4Atwx5pAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADqcvVkPvHfv3qH1Dh48OLTenj17htbbtm3b0HpHjhwZWo/NN/o9ctFFFw2td9lllw2tN9rFF188tN7u3buH1mPzHT58eGi9/fv3D6333Oc+d2g9tp7t27cPrTf6c+5UX0/YdgJGskcaAAAAOgjSAAAA0EGQBgAAgA6CNAAAAHQQpAEAAKCDIA0AAAAdBGkAAADoIEgDAABAB0EaAAAAOgjSAAAA0EGQBgAAgA6CNAAAAHQQpAEAAKCDIA0AAAAdBGkAAADoIEgDAABAB0EaAAAAOgjSAAAA0EGQBgAAgA7VWlvvsutecCs4cODA0HrLy8tD6+3YsWNovdHPdwPUZjfgGM6oPrFr166h9a688sqh9c4999yh9Y4cOTK03gbQJzotLS0NrXf48OGh9Q4ePDi03rZt24bWOw3oE51Gf86N3tYZ3ScOHTo0tN5p4FTsE7Bl2CMNAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQAdBGgAAADoI0gAAANBBkAYAAIAOgjQAAAB0EKQBAACggyANAAAAHQRpAAAA6CBIAwAAQIdqrW12GwAAAOC0YY80AAAAdBCkAQAAoIMgDQAAAB0EaQAAAOggSAMAAEAHQRoAAAA6CNIAAADQQZAGAACADoI0AAAAdPj/AfPvztCmAC3JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1000x400 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"Linear classifier on raw input data:\")\n",
    "model = LogisticRegression(C=0.005, max_iter=1000)\n",
    "model.fit(X_train, y_train)\n",
    "_, y_pred_test = evaluate_model(model)\n",
    "ind_wrong = (y_pred_test != y_test)\n",
    "\n",
    "n_rows = 2\n",
    "n_cols = 5\n",
    "fig, axes = plt.subplots(n_rows, n_cols, figsize=(2*n_cols,2*n_rows), dpi=100)\n",
    "for i, (img, y_true, y_pred) in enumerate(zip(X_test[ind_wrong],\n",
    "                                              y_test[ind_wrong],\n",
    "                                              y_pred_test[ind_wrong])):\n",
    "    if i >= (n_rows * n_cols):\n",
    "        break\n",
    "    ax = axes[i//n_cols][i%n_cols]\n",
    "    ax.imshow(img.reshape([8,8]), cmap='binary')\n",
    "    ax.set_title(f\"T: {y_true}, P: {y_pred}\")\n",
    "for ax in axes.flatten():\n",
    "    ax.axis(\"off\")\n",
    "fig.suptitle(\"Wrong predictions\")\n",
    "fig.tight_layout()\n",
    "\n",
    "# Save the model\n",
    "pickle.dump(model, open('pretrained_lr.model', 'wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}