SGD-MNIST1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using SGD on MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ... about machine learning (a reminder from lesson 1)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "The good news is that modern machine learning can be distilled down to a couple of key techniques that are of very wide applicability. Recent studies have shown that the vast majority of datasets can be best modeled with just two methods:\n",
    "\n",
    "Ensembles of decision trees (i.e. Random Forests and Gradient Boosting Machines), mainly for structured data (such as you might find in a database table at most companies). We looked at random forests in depth as we analyzed the Blue Book for Bulldozers dataset.\n",
    "\n",
    "Multi-layered neural networks learnt with SGD (i.e. shallow and/or deep learning), mainly for unstructured data (such as audio, vision, and natural language)\n",
    "\n",
    "In this lesson, we will start on the 2nd approach (a neural network with SGD) by analyzing the MNIST dataset. You may be surprised to learn that logistic regression is actually an example of a simple neural net!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## About The Data"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "In this lesson, we will be working with MNIST, a classic data set of hand-written digits. Solutions to this problem are used by banks to automatically recognize the amounts on checks, and by the postal service to automatically recognize zip codes on mail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports and data"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "We will be using the fastai library, which is still in pre-alpha. If you are accessing this course notebook, you probably already have it downloaded, as it is in the same Github repo as the course materials.\n",
    "\n",
    "We use symbolic links (often called symlinks) to make it possible to import these from your current directory. For instance, I ran:\n",
    "\n",
    "ln -s ../../fastai"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fastai.imports import *\n",
    "from fastai.torch_imports import *\n",
    "from fastai.io import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = 'data/mnist/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.makedirs(path, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "URL='http://deeplearning.net/data/mnist/'\n",
    "FILENAME='mnist.pkl.gz'\n",
    "\n",
    "def load_mnist(filename):\n",
    "    return pickle.load(gzip.open(filename, 'rb'), encoding='latin-1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_data(URL+FILENAME, path+FILENAME)\n",
    "((x, y), (x_valid, y_valid), _) = load_mnist(path+FILENAME)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(numpy.ndarray, (50000, 784), numpy.ndarray, (50000,))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(x), x.shape, type(y), y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalize"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Many machine learning algorithms behave better when the data is normalized, that is when the mean is 0 and the standard deviation is 1. We will subtract off the mean and standard deviation from our training set in order to normalize the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.13044983, 0.3072898, -3.1638146e-07, 0.99999934)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean = x.mean()\n",
    "std = x.std()\n",
    "\n",
    "x=(x-mean)/std\n",
    "mean, std, x.mean(), x.std()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Note that for consistency (with the parameters we learn when training), we subtract the mean and standard deviation of our training set from our validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.005850922, 0.99243325)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_valid = (x_valid-mean)/std\n",
    "x_valid.mean(), x_valid.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Look at the data"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "In any sort of data science work, it's important to look at your data, to make sure you understand the format, how it's stored, what type of values it holds, etc. To make it easier to work with, let's reshape it into 2d images from the flattened 1d format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show(img, title=None):\n",
    "    plt.imshow(img, cmap=\"gray\")\n",
    "    if title is not None: plt.title(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plots(ims, figsize=(12,6), rows=2, titles=None):\n",
    "    f = plt.figure(figsize=figsize)\n",
    "    cols = len(ims)//rows\n",
    "    for i in range(len(ims)):\n",
    "        sp = f.add_subplot(rows, cols, i+1)\n",
    "        sp.axis('Off')\n",
    "        if titles is not None: sp.set_title(titles[i], fontsize=16)\n",
    "        plt.imshow(ims[i], cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000, 784)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_valid.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000, 28, 28)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_imgs = np.reshape(x_valid, (-1,28,28)); x_imgs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAEICAYAAACQ6CLfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADl1JREFUeJzt3X3IXHV6xvHr8mUxMaLRVE1iNLtPhL5htAapKEVd3NqtEFdwMWBJoyVrWaGrVSpBUBTBlu5qK1SJGMzirltN3FVWRcXa+gZifGGNGzer4saYJ09Qi4mou43e/eM5WR7jzG8mM2fmTHJ/P/AwM+eeM+dmyJVzzvzOzM8RIQD57Nd0AwCaQfiBpAg/kBThB5Ii/EBShB9IivADSRF+tGT7btvjtrfb3mj775ruCfUyF/mgFdt/IumNiPit7T+U9N+S/joiXmy2M9SFPT9aiojXIuK3ux5Wf2MNtoSaEX60Zfs/bH8s6XVJ45Iebrgl1IjDfhTZ3l/SqZLOkPTPEfF/zXaEurDnR1FEfBYRz0g6RtLfN90P6kP40a0DxDn/PoXw40tsH2n7QtszbO9v+y8lLZH0X033hvpwzo8vsf0HktZIWqjJHcRvJP17RNzRaGOoFeEHkuKwH0iK8ANJEX4gKcIPJHXAMDdmm08XgQGLCHfzvL72/LbPsf0r22/Yvrqf1wIwXD0P9VXXfG+UdLakzZJekLQkIn5ZWIc9PzBgw9jzn6LJ73u/FRG/k/QTSYv7eD0AQ9RP+OdKemfK483Vsi+wvdz2Otvr+tgWgJr184Ffq0OLLx3WR8RKSSslDvuBUdLPnn+zpHlTHh8jaUt/7QAYln7C/4Kk421/1fZXJF0o6cF62gIwaD0f9kfETtuXSXpU0v6SVkXEa7V1BmCghvqtPs75gcEbykU+APZehB9IivADSRF+ICnCDyRF+IGkCD+QFOEHkiL8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSIrwA0kRfiApwg8kNdQputGbhQsXFuuXX35529rY2Fhx3enTpxfrK1asKNYPPfTQYv2RRx5pW9uxY0dxXQwWe34gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSIpZekfAjBkzivVNmzYV64cddlid7dTq3XffbVsrXZ8gSWvWrKm7nRS6naW3r4t8bL8taYekzyTtjIhF/bwegOGp4wq/MyPivRpeB8AQcc4PJNVv+EPSY7ZftL281RNsL7e9zva6PrcFoEb9HvafFhFbbB8p6XHbr0fEU1OfEBErJa2U+MAPGCV97fkjYkt1u03STyWdUkdTAAav5/DbPtj2IbvuS/qGpPV1NQZgsHoe57f9NU3u7aXJ04cfR8SNHdbhsL+FQw45pFh/+OGHi/X333+/be3ll18urnvSSScV68cdd1yxPm/evGJ92rRpbWsTExPFdU899dRivdP6WQ18nD8i3pJU/pUJACOLoT4gKcIPJEX4gaQIP5AU4QeS4iu96MusWbOK9auuuqqnmiQtW7asWF+9enWxnlW3Q33s+YGkCD+QFOEHkiL8QFKEH0iK8ANJEX4gKaboRl/ee6/8263PPvts21qncf5OXzdmnL8/7PmBpAg/kBThB5Ii/EBShB9IivADSRF+ICnG+dGXmTNnFusrVqzo+bXnzJnT87rojD0/kBThB5Ii/EBShB9IivADSRF+ICnCDyTF7/ajaOHC8kTM9913X7G+YMGCtrWNGzcW1z377LOL9XfeeadYz6q23+23vcr2Ntvrpyw73Pbjtn9d3Zav9AAwcro57L9L0jm7Lbta0hMRcbykJ6rHAPYiHcMfEU9J+mC3xYsl7foNpdWSzqu5LwAD1uu1/UdFxLgkRcS47SPbPdH2cknLe9wOgAEZ+Bd7ImKlpJUSH/gBo6TXob4J27MlqbrdVl9LAIah1/A/KGlpdX+ppAfqaQfAsHQc57d9j6QzJM2SNCHpWkk/k3SvpGMlbZJ0QUTs/qFgq9fisH/ELF26tFi//vrri/V58+YV65988knb2rnnnltc98knnyzW0Vq34/wdz/kjYkmb0tf3qCMAI4XLe4GkCD+QFOEHkiL8QFKEH0iKn+7eB8yYMaNt7corryyue8011xTr++1X3j988EF5hPf0009vW3v99deL62Kw2PMDSRF+ICnCDyRF+IGkCD+QFOEHkiL8QFKM8+8D7rrrrra1888/v6/XXrNmTbF+yy23FOuM5Y8u9vxAUoQfSIrwA0kRfiApwg8kRfiBpAg/kBTj/PuAsbGxgb32bbfdVqw/99xzA9s2Bos9P5AU4QeSIvxAUoQfSIrwA0kRfiApwg8kxTj/PuCxxx5rW1u4cOHAXlvqfB3ATTfd1La2ZcuWnnpCPTru+W2vsr3N9vopy66z/a7tV6q/bw62TQB16+aw/y5J57RYfnNEnFj9PVxvWwAGrWP4I+IpSeU5mQDsdfr5wO8y27+oTgtmtnuS7eW219le18e2ANSs1/DfJmlM0omSxiV9v90TI2JlRCyKiEU9bgvAAPQU/oiYiIjPIuJzSXdIOqXetgAMWk/htz17ysNvSVrf7rkARpMjovwE+x5JZ0iaJWlC0rXV4xMlhaS3JX0nIsY7bswubww9mTZtWtva3XffXVz35JNPLtaPPfbYnnraZevWrW1ry5YtK6776KOP9rXtrCLC3Tyv40U+EbGkxeI797gjACOFy3uBpAg/kBThB5Ii/EBShB9IquNQX60bY6hv6A466KBi/YADygM+27dvr7OdL/j000+L9SuuuKJYv/322+tsZ5/R7VAfe34gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSIpxfhSdcMIJxfrNN99crJ955pk9b3vTpk3F+vz583t+7X0Z4/wAigg/kBThB5Ii/EBShB9IivADSRF+ICnG+UfA9OnTi/WPP/54SJ3suZkz287UJklatWpV29rixYv72vbcuXOL9fHxjr8mv09inB9AEeEHkiL8QFKEH0iK8ANJEX4gKcIPJNVxll7b8yT9UNLRkj6XtDIi/s324ZL+U9J8TU7T/e2I+N/Btbr3GhsbK9afeeaZYv2hhx4q1tevX9+21mms+5JLLinWDzzwwGK901j7ggULivWSN998s1jPOo5fl272/Dsl/WNE/JGkP5f0Xdt/LOlqSU9ExPGSnqgeA9hLdAx/RIxHxEvV/R2SNkiaK2mxpNXV01ZLOm9QTQKo3x6d89ueL+kkSc9LOioixqXJ/yAkHVl3cwAGp+M5/y62Z0haK+l7EbHd7uryYdleLml5b+0BGJSu9vy2D9Rk8H8UEfdXiydsz67qsyVta7VuRKyMiEURsaiOhgHUo2P4PbmLv1PShoj4wZTSg5KWVveXSnqg/vYADEo3h/2nSfobSa/afqVatkLSTZLutX2JpE2SLhhMi3u/Cy4ovzVHH310sX7xxRfX2c4e6XR6189Xwj/66KNi/dJLL+35tdFZx/BHxDOS2v0L+Hq97QAYFq7wA5Ii/EBShB9IivADSRF+ICnCDyTV9eW96N0RRxzRdAsDs3bt2mL9hhtuaFvbtq3lRaG/t3Xr1p56QnfY8wNJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUkzRPQSdfv76rLPOKtYvuuiiYn3OnDltax9++GFx3U5uvfXWYv3pp58u1nfu3NnX9rHnmKIbQBHhB5Ii/EBShB9IivADSRF+ICnCDyTFOD+wj2GcH0AR4QeSIvxAUoQfSIrwA0kRfiApwg8k1TH8tufZftL2Btuv2f6Havl1tt+1/Ur1983BtwugLh0v8rE9W9LsiHjJ9iGSXpR0nqRvS/ooIv61641xkQ8wcN1e5NNxxp6IGJc0Xt3fYXuDpLn9tQegaXt0zm97vqSTJD1fLbrM9i9sr7I9s806y22vs72ur04B1Krra/ttz5D0P5JujIj7bR8l6T1JIekGTZ4aXNzhNTjsBwas28P+rsJv+0BJP5f0aET8oEV9vqSfR8Sfdngdwg8MWG1f7LFtSXdK2jA1+NUHgbt8S9L6PW0SQHO6+bT/dElPS3pV0ufV4hWSlkg6UZOH/W9L+k714WDptdjzAwNW62F/XQg/MHh8nx9AEeEHkiL8QFKEH0iK8ANJEX4gKcIPJEX4gaQIP5AU4QeSIvxAUoQfSIrwA0kRfiCpjj/gWbP3JP1myuNZ1bJRNKq9jWpfEr31qs7ejuv2iUP9Pv+XNm6vi4hFjTVQMKq9jWpfEr31qqneOOwHkiL8QFJNh39lw9svGdXeRrUvid561UhvjZ7zA2hO03t+AA0h/EBSjYTf9jm2f2X7DdtXN9FDO7bftv1qNe14o/MLVnMgbrO9fsqyw20/bvvX1W3LORIb6m0kpm0vTCvf6Hs3atPdD/2c3/b+kjZKOlvSZkkvSFoSEb8caiNt2H5b0qKIaPyCENt/IekjST/cNRWa7X+R9EFE3FT9xzkzIv5pRHq7Tns4bfuAems3rfzfqsH3rs7p7uvQxJ7/FElvRMRbEfE7ST+RtLiBPkZeRDwl6YPdFi+WtLq6v1qT/3iGrk1vIyEixiPiper+Dkm7ppVv9L0r9NWIJsI/V9I7Ux5vVoNvQAsh6THbL9pe3nQzLRy1a1q06vbIhvvZXcdp24dpt2nlR+a962W6+7o1Ef5WUwmN0njjaRHxZ5L+StJ3q8NbdOc2SWOanMNxXNL3m2ymmlZ+raTvRcT2JnuZqkVfjbxvTYR/s6R5Ux4fI2lLA320FBFbqtttkn6qydOUUTKxa4bk6nZbw/38XkRMRMRnEfG5pDvU4HtXTSu/VtKPIuL+anHj712rvpp635oI/wuSjrf9VdtfkXShpAcb6ONLbB9cfRAj2wdL+oZGb+rxByUtre4vlfRAg718wahM295uWnk1/N6N2nT3jVzhVw1l3CJpf0mrIuLGoTfRgu2vaXJvL01+3fnHTfZm+x5JZ2jyK58Tkq6V9DNJ90o6VtImSRdExNA/eGvT2xnaw2nbB9Rbu2nln1eD712d093X0g+X9wI5cYUfkBThB5Ii/EBShB9IivADSRF+ICnCDyT1/x2VQ9c6BSuMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(x_imgs[0], y_valid[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000,)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_valid.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_valid[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can look at part of an image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.42452, -0.42452, -0.42452, -0.42452,  0.17294],\n",
       "       [-0.42452, -0.42452, -0.42452,  0.78312,  2.43567],\n",
       "       [-0.42452, -0.27197,  1.20261,  2.77889,  2.80432],\n",
       "       [-0.42452,  1.76194,  2.80432,  2.80432,  1.73651],\n",
       "       [-0.42452,  2.20685,  2.80432,  2.80432,  0.40176]], dtype=float32)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_imgs[0,10:15,10:15]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAACPxJREFUeJzt3UGIVYUex/Hfz3lKgg9aPBfhyDOi4knwFEQCd9LCMgpcKRQtgtm8wCCMAjfRPtrUYqhIKIqoFhE9QiiJoFepWeSbCokeSYE+MqpNMvV/i7kL6TneM95z5tzz4/uBgbl6uP6K+c65987ljKtKADKt6XsAgO4QOBCMwIFgBA4EI3AgGIEDwQgcCEbgQDACB4L9qYs7tc3b4zA469ev73tCYxcvXtTi4qLHHddJ4MAQ3XTTTX1PaOyrr75qdBwP0YFgBA4EI3AgGIEDwQgcCEbgQDACB4IROBCMwIFgBA4EI3AgGIEDwQgcCEbgQDACB4IROBCsUeC299j+0vYZ2490PQpAO8YGbntG0lOSbpe0VdIB21u7HgZgck3O4Dslnamqr6vqoqSXJd3d7SwAbWgS+CZJ315y++zozwBMuSYXXbzclRv/76qptuckzU28CEBrmgR+VtLmS27PSvrujwdV1bykeYnLJgPToslD9I8l3Wj7etvrJO2X9Ea3swC0YewZvKoWbT8g6W1JM5Keq6rTnS8DMLFGv/igqt6S9FbHWwC0jHeyAcEIHAhG4EAwAgeCETgQjMCBYAQOBCNwIBiBA8EIHAhG4EAwAgeCETgQjMCBYAQOBCNwIBiBA8EaXdEFuFr33Xdf3xMae/zxx/ue0NjevXsbHccZHAhG4EAwAgeCETgQjMCBYAQOBCNwIBiBA8EIHAhG4EAwAgeCETgQjMCBYAQOBCNwIBiBA8EIHAg2NnDbz9k+Z/vz1RgEoD1NzuDPS9rT8Q4AHRgbeFW9J+mHVdgCoGU8BweCtXZVVdtzkubauj8Ak2st8KqalzQvSbarrfsFcPV4iA4Ea/JjspckfSDpZttnbd/f/SwAbRj7EL2qDqzGEADt4yE6EIzAgWAEDgQjcCAYgQPBCBwIRuBAMAIHghE4EIzAgWAEDgQjcCAYgQPBCBwIRuBAMAIHgrV2Tbah2rBhQ98TVuTQoUN9T1iRw4cP9z2hsTVrhnO+W7duXaPjhvNfBGDFCBwIRuBAMAIHghE4EIzAgWAEDgQjcCAYgQPBCBwIRuBAMAIHghE4EIzAgWAEDgQjcCAYgQPBxgZue7Ptd20v2D5t++BqDAMwuSaXbFqU9FBVnbT9Z0knbB+tqn93vA3AhMaewavq+6o6Ofr8Z0kLkjZ1PQzA5Fb0HNz2FknbJX3YxRgA7Wp8VVXbGyS9JunBqvrpMn8/J2muxW0AJtQocNtrtRT3i1X1+uWOqap5SfOj46u1hQCuWpNX0S3pWUkLVfVE95MAtKXJc/Bdku6VtNv2qdHHHR3vAtCCsQ/Rq+p9SV6FLQBaxjvZgGAEDgQjcCAYgQPBCBwIRuBAMAIHghE4EIzAgWAEDgQjcCAYgQPBCBwIRuBAMAIHghE4EIzAgWCNr6qa6siRI31PWJF9+/b1PSHWq6++2veExi5cuNDoOM7gQDACB4IROBCMwIFgBA4EI3AgGIEDwQgcCEbgQDACB4IROBCMwIFgBA4EI3AgGIEDwQgcCEbgQLCxgdu+xvZHtj+1fdr2Y6sxDMDkmlyy6VdJu6vqF9trJb1v+59V9a+OtwGY0NjAq6ok/TK6uXb0UV2OAtCORs/Bbc/YPiXpnKSjVfVht7MAtKFR4FX1W1VtkzQraaftW/54jO0528dtH297JICrs6JX0avqR0nHJO25zN/NV9WOqtrR0jYAE2ryKvpG29eOPl8v6TZJX3Q9DMDkmryKfp2kI7ZntPQN4ZWqerPbWQDa0ORV9M8kbV+FLQBaxjvZgGAEDgQjcCAYgQPBCBwIRuBAMAIHghE4EIzAgWAEDgQjcCAYgQPBCBwIRuBAMAIHghE4EKzJFV2i3XDDDX1PwJR4+umn+57Q2Pnz5xsdxxkcCEbgQDACB4IROBCMwIFgBA4EI3AgGIEDwQgcCEbgQDACB4IROBCMwIFgBA4EI3AgGIEDwQgcCNY4cNsztj+x/WaXgwC0ZyVn8IOSFroaAqB9jQK3PStpr6Rnup0DoE1Nz+BPSnpY0u8dbgHQsrGB275T0rmqOjHmuDnbx20fb20dgIk0OYPvknSX7W8kvSxpt+0X/nhQVc1X1Y6q2tHyRgBXaWzgVfVoVc1W1RZJ+yW9U1X3dL4MwMT4OTgQbEW/2aSqjkk61skSAK3jDA4EI3AgGIEDwQgcCEbgQDACB4IROBCMwIFgBA4EI3AgGIEDwQgcCEbgQDACB4IROBCMwIFgBA4Ec1W1f6f2eUn/aflu/yLpvy3fZ5eGtHdIW6Vh7e1q61+rauO4gzoJvAu2jw/piq1D2jukrdKw9va9lYfoQDACB4INKfD5vges0JD2DmmrNKy9vW4dzHNwACs3pDM4gBUaROC299j+0vYZ24/0vedKbD9n+5ztz/veMo7tzbbftb1g+7Ttg31vWo7ta2x/ZPvT0dbH+t7UhO0Z25/YfrOPf3/qA7c9I+kpSbdL2irpgO2t/a66oucl7el7REOLkh6qqr9JulXSP6b4/+2vknZX1d8lbZO0x/atPW9q4qCkhb7+8akPXNJOSWeq6uuquqil33B6d8+bllVV70n6oe8dTVTV91V1cvT5z1r6QtzU76rLqyW/jG6uHX1M9QtItmcl7ZX0TF8bhhD4JknfXnL7rKb0i3DIbG+RtF3Sh/0uWd7o4e4pSeckHa2qqd068qSkhyX93teAIQTuy/zZVH/nHhrbGyS9JunBqvqp7z3LqarfqmqbpFlJO23f0vem5di+U9K5qjrR544hBH5W0uZLbs9K+q6nLXFsr9VS3C9W1et972miqn7U0m+5nebXOnZJusv2N1p6Wrnb9gurPWIIgX8s6Ubb19teJ2m/pDd63hTBtiU9K2mhqp7oe8+V2N5o+9rR5+sl3Sbpi35XLa+qHq2q2araoqWv2Xeq6p7V3jH1gVfVoqQHJL2tpReBXqmq0/2uWp7tlyR9IOlm22dt39/3pivYJeleLZ1dTo0+7uh71DKuk/Su7c+09E3/aFX18qOnIeGdbECwqT+DA7h6BA4EI3AgGIEDwQgcCEbgQDACB4IROBDsf/DAx4Xphc/GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(x_imgs[0,10:15,10:15])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAF0CAYAAADGqzQSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu4lnPa//Hzq1raktCWSkU9QiV7a2ZsSqWO0siMyNCOEkY2ox0ifjFjVERmJrIrtLfJblI26UGm6FdKZSZFO0raqrVa398f93qOp1/nebeu1X2v+1r3+r5fx9ERH9fmxLWudXat+7y+znsvAAAAQIgOi7sAAAAAIC40wwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAgWzXDMnHPtnHNznHMbnHN7nHPfOecmO+dOjrs24GCcc+c75951zm1yzm1zzi10zvWKuy4gCufcpc65D51zOwqv38+dcxfFXRdwMM65C51z85xzu51zW5xzLzjnasVdV7ajGY5fDRH5l4jcJCKXiMhgEWkuIp845xrEWRiQjHPuNBGZLSIVRKSviFwuIgtE5GnnXP84awOK4py7QURelcS9t6uIXCEiU0Skcpx1AQfjnPuViLwrIlslcc/9o4j8WkTec84dHmdt2c6x6Ebp45xrKiLLReQO7/1f464HOJBz7v+IyB0iUsN7v2O//BMR8d77c2MrDjgI51xDEVkmIoO996PjrQaIzjk3W0Qaikgz731+YXamiHwmIgO890/GWF5W48lw6bS58Pe8WKsAksuRxPW5+4B8q3BfQenWS0QKROSpuAsBiukcEfnn/zTCIiLe+wWS6Bm6xlZVGcA3rVLCOVfOOZfjnDtRRP4mIhtE5OWYywKSebbw98ecc3Wdc9Wdc31F5GIRGRVfWUCRciXxk7crnXPfOOfynXOrnHMD4i4MKMI+Edlr5HtE5JQM11Km8DGJUsI597mItC7821Ui0tl7vyzGkoCDKvzx3AwRqVcY5YlIf+/90/FVBRycc265iNSVRAMxRES+kcRnhvuJyK3e+zExlgck5Zz7TBIfQzt7v6yBiPxHRPK893xu+BDRDJcSzrn/EpEjRKSRJD6LWUtEcr33q+OsC7AU/gTjPUl89vJxSXxcoouI9BeR67z3E2MsD0jKObdCRE4Ukcu999P3y98SkVYiUsfzjRGlkHPuahF5UUQeFJHHJDGA/3cROU8SzXClGMvLajTDpZBzrrqIrBaRl733/WIuB1Ccc1NE5HRJDHLk7ZdPFJF2IlLTe18QV31AMs65/5bEZy+P8N5v3y8fKCKPikg97/26uOoDDsY5N0ISD8wqiogXkVdEpIqInOK9bxRnbdmMzwyXQt77rZL4qESTuGsBkjhVRL7cvxEu9JmIHC0iNTNfEhDJ0iS5K/ydP8Sh1PLe3y0ix4jIaZL4KUZ3SfykY16shWU5muFSqPAF2s0k8Vk2oDTaICItnXM5B+Rni8gvIrIl8yUBkcwo/L3dAXk7EfnOe78hw/UAxeK93+m9/7/e+43OufaS6Bd4O0oKysddQOicczNEZKGILBaRbSJykogMFJF8EeEdwyitxkpikYLXnXNPSuIzw51FpLuIjPLeWxPPQGnwpojMFZG/OeeOEZF/i0g3SSx61DPOwoCDcc61EpEOkugZRBJvRrlTRP7svZ8fW2FlAJ8Zjplz7i4R+Z2INJbEu1vXisj7IjKS4TmUZs65DiJylyRWTKwoiZ9k/F1E/ua93xdnbcDBOOeOEJGRkmiCj5LEq9Ye8t5PirUw4CCcc80l8erVU0TkcCkcYPbeT4i1sDKAZhgAAADB4jPDAAAACBbNMAAAAIJFMwwAAIBg0QwDAAAgWBl9tZpzjmk9pMx774reKr24dpEOmb52uW6RDtxzka2iXrs8GQYAAECwaIYBAAAQLJphAAAABItmGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABItmGAAAAMGiGQYAAECwMroCHQAApUnVqlVV1rt3b5V16dLF3L9z584q27FjR+qFAcgYngwDAAAgWDTDAAAACBbNMAAAAIJFMwwAAIBg0QwDAAAgWLxNAgAQrGuvvVZlo0aNirx/8+bNVfbpp5+mVBOAzOLJMAAAAIJFMwwAAIBg0QwDAAAgWDTDAAAACBYDdClo0aKFygYOHGhu27hxY5VVrlxZZUOGDFHZkUceqbK33nrLPM/27dvNHABCd91116ls9OjRKsvLy1PZI488Yh5z4cKFKdcFIF48GQYAAECwaIYBAAAQLJphAAAABItmGAAAAMFy3vvMncy5zJ0szapWraqyNWvWqKx69eqZKEe+//57M7cG+KZOnVrS5WSU995l+pzZfO1arOu0a9eu5ratWrVSWW5ursqsr5EtW7aorHbt2uZ5NmzYoLJnn31WZf/4xz9Utm/fPvOYpU2mr92ydt0WR+fOnVU2Y8YMle3atUtl99xzj8qKsypdWcM9F9kq6rXLk2EAAAAEi2YYAAAAwaIZBgAAQLBohgEAABAsBugiqlatmsrefPNNlW3evNncf9GiRSqzBpMaNGigsuOPP15llSpVMs+zceNGlZ177rmRtssWDHMUz3HHHaeymTNnqsy6HpPZtm2byqxrvEKFCiqzvpZERGrWrKmyWrVqqeyqq65S2Ycffqiy9evXm+eJEwN06ZeTk2PmEyZMUFn37t1VNmfOHJW1adMm9cLKEO65yFYM0AEAAABFoBkGAABAsGiGAQAAECyaYQAAAASLAboscMwxx6jszjvvNLe18p49e6rsueeeS72wmDDMUTwLFy5UWYsWLVQ2e/Zsc//bb79dZT/++KPKrBXkiuPYY49V2VtvvaWypk2bqmzQoEEqe+KJJ1KqpyQwQJd+Q4cONfMRI0ao7MUXX1RZr169VJafn596YWUI99zU1alTR2U33nijua2V5+XlqcxaBffBBx9UmfU9QERk7dq1Zl6WMEAHAAAAFIFmGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABIu3SWSpzp07m7m1zO5jjz2msltvvTXtNWUKk83JWRPL33//vcomT56ssquvvto85r59+1Iv7BBNnDhRZVdeeaXKWrdurbIvvviiRGpKBW+TSM0ZZ5yhsnnz5pnbrl69WmXNmzdXWZzXd7bgnls8jRo1Utm4ceNU1rZt20yUI3v27DHz888/X2XJ3jyRrXibBAAAAFAEmmEAAAAEi2YYAAAAwaIZBgAAQLDKx10AinbUUUepbMiQIZH3r1u3bjrLQSnWsmVLlTmn5wfWrVunsrgHic455xyVde/eXWVz585VmfXvXRoH6BDdYYfpZzXWsts5OTnm/q+//rrK4r7GUfbUq1dPZUuWLFFZ+fK63Ro1apR5zMcffzzSeZo1a6ayv/zlLyqrXr26eR5rkNq6D//444/m/mUJT4YBAAAQLJphAAAABItmGAAAAMGiGQYAAECwWIGulGnRooXKpkyZorImTZqY+69YsUJl1io3a9euPYTqSgdWQyqegoIClW3atEllZ511lrn/mjVr0lpPtWrVzHz+/PkqW7lypcqslfKsFZ+WLl16CNWVLFagiy7qaorJ3HLLLSobO3ZsSjWFintucmPGjFFZv379VNa3b1+VPf/882mvZ8CAASobPXq0uW25cuVUtnz5cpVZQ3Xbtm07hOoyjxXoAAAAgCLQDAMAACBYNMMAAAAIFs0wAAAAgsUAXYyuvfZald1///0qO/7441W2e/du85idOnVSmbViVzZjmKN4hg8frrK7775bZV9//bW5f7t27VSWygDmu+++a+a/+c1vVNa6dWuVWas7ZQsG6KLr2bOnyp5++mmVzZ4929y/Q4cOKmMFukPDPVfkiCOOMHNryHfChAkqs1ZPzJRk9/YTTzwx0v7WSnm33357SjVlCgN0AAAAQBFohgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLB4m0SaVa1a1czvuOMOlQ0bNkxlhx2m/3yyZcsWleXm5prnsZZSLGuYbC6eihUrquy5555TWbdu3cz9V61apbILLrhAZevXr1fZk08+qbLrr7/ePM+dd96pMmuKOZvxNglb+fLlVbZs2TKVNWjQQGUnnHCCecziLN2Mg+Oem3y5+k8++URlbdu2Vdl7772X9pqi6tq1q5lPnz5dZVZPuHXrVpVZb6LYvHnzIVRXsnibBAAAAFAEmmEAAAAEi2YYAAAAwaIZBgAAQLD01AJS8uyzz5r5b3/720j7T506VWWjR49WWQiDckiPX375RWV9+vRRWc2aNc39rWWSP/jgA5VNmTJFZT169FDZtGnTzPOUtWE5RGcNbzZu3Fhl/fv3V1ncg3Lt27dXWefOnVX29ttvq8xamtz6ekX8WrVqFXnbRYsWlWAlxffmm2+auTUcbX3dWdfkzp07Uy+sFOHJMAAAAIJFMwwAAIBg0QwDAAAgWDTDAAAACBYDdGlmffi8OMaNG6ey+fPnp3RM4EDbt29XWZcuXcxthw8frrJbb71VZYMGDYp07scffzzSdghH/fr1I22Xk5NTwpUkd91115m5tcqitepjv379VGat7DVz5kzzPL169SqiQpSkefPmmXlBQYHK/vnPf6qsU6dOKrNW7SwJTZs2NXPrOm3Xrp3KKleurLKyNujJk2EAAAAEi2YYAAAAwaIZBgAAQLBohgEAABAsBujSzFpRSESkRYsWh7y/NVT30EMPmfuvW7cu0nmAA23bts3M77nnHpW1bdtWZSeffHKk87Rp08bMkw2ooOxr0qRJpO0ytfJm9erVVfboo4+a21pDSPn5+Sqzhqpyc3NVZq3aKMIAXdyWLl1q5m+88YbKrGHkZcuWqcxalVDEXqVzzpw5KqtXr57KrGE5axVbEZE6deqozLp2X331VXP/soQnwwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFjOe5+5kzmXuZPFpFKlSmb+4osvqqx169Yqi7oS04YNG8y8Z8+eKnvnnXciHTNbeO9dps8ZwrWbTIcOHVQ2Y8YMlVWoUCHS8fbu3WvmN954o8omTJgQ6ZjZItPXbrZct7NmzVJZq1atVFa3bt1MlGOusJhsgM66t48ZM0Zla9asUZk1QHXqqaea54lz9T3uuclZ3/NHjhypsltuuSWl82zZskVlNWrUSOmYliuuuEJl1kBftoh67fJkGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABItmGAAAAMFiOeY02717t5lfffXVKitfXv/nT7Yk7oFq165t5taU/2233aayp556KtJ5gAsvvFBl1ltounbtqjJrAtpavlTEXnb8xx9/VNnrr79u7o/sdfbZZ6ss2VtHSpt169ap7LjjjlPZ3//+d5WdfvrpKitrb/8p66zv+dbbSCZPnqwyqy9IplatWpG2y8vLU5n19SUicsIJJ6hs165dkWsqS3gyDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgsVyzKXMaaedprJRo0apzBpqSsZaBrRhw4bFqqs0YWnQkmFdeyIiCxYsUJk17GYNjVis5T5FRJ5++mmVOaf/Vzdv3lxl1jVeGrEcs80aLuvUqZPKSmI5Zusas67lv/71rymdx/pe++STT6psyJAh5v7bt29P6fyp4J6b3V544QUztwb42rdvr7J333037TVlCssxAwAAAEWgGQYAAECwaIYBAAAQLJphAAAABIsV6CKqXLmyykpipZbFixerrFu3bip75plnzP27dOmisvr166usTp06Klu/fn2UElFGVatWzcytlRKnTp16yOeZMmWKmTdo0EBlDz/8sMpat26tsmwZoEN01atXV5k1CPTiiy+a+1vX7ZVXXqmyGjVqqKxDhw5RShQRkZ07d6ps3rx5Kvvzn/+ssrlz50Y+D5AJjRs3jruEWPBkGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABIsBOoP1AXJrIGLWrFkqW7JkiXlMazitd+/eKqtQoYLK6tWrp7ImTZqY57F88803kepB2Fq2bGnmGzZsUJn19ZCqsWPHqqxv374qGzBggMpmzJiR9nqQOYsWLVJZnz59VGatmGVlqdq2bZvKkg1+PvDAAyr79ttv014TcKh27NgRdwmlHk+GAQAAECyaYQAAAASLZhgAAADBohkGAABAsBigM1xxxRUqq127tsp69eqV9nM751TmvY+8v/VB+X79+qVUE8JgrVQoIvLZZ59l5Px79+5V2U8//aSyX/3qVyqzVhHbsmVLegpDiZs0aZLKrJU3V65cqbJy5cqZx0yWH2jixIkqW716tcqsQWQgG3z44YdmfsMNN6isZs2aJV1OqcSTYQAAAASLZhgAAADBohkGAABAsGiGAQAAECyaYQAAAASLt0kYjj766LhL+P9MmzZNZSNGjDC33bRpk8qs5XSBAyV7a0lubq7KrrzySpXNmTNHZVWrVlVZTk6OeZ5mzZqp7Mwzz1TZE088oTLeHJHdfv75Z5VdfPHFMVQClD2HHWY/97TeXmX1ECHgyTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAgWA3SGIUOGqGz27Nkq69Gjh8rq1q1rHtMaELE8/vjjKvvoo49Ulp+fH+l4QFTLli0zc2upY2v53M2bN6usOAN01jDHxx9/rLLhw4eb+wMAtIKCAjNPNjQdIp4MAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYDFAZ8jLy1PZO++8EykDstXbb79t5mPHjlWZtSpdy5YtUzr/0KFDVfbMM8+ojNXmAKBkXHLJJSobN25cDJVkFk+GAQAAECyaYQAAAASLZhgAAADBohkGAABAsBigAyAiIhs3bjTzP/7xjxmuBACQLjt27Ii8bfnyYbaFPBkGAABAsGiGAQAAECyaYQAAAASLZhgAAADBohkGAABAsJz3PnMncy5zJ0OZ5b13mT4n1y7SIdPXLtct0oF7bnarXr26mVtL2+/evVtlVapUSXtNmRL12uXJMAAAAIJFMwwAAIBg0QwDAAAgWDTDAAAACBYDdMg6DHMgWzFAh2zEPRfZigE6AAAAoAg0wwAAAAgWzTAAAACCRTMMAACAYGV0gA4AAAAoTXgyDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMOlhHPuUufch865Hc65bc65z51zF8VdFxCVc+5t55x3zj0Qdy1AMs65Cwqv0wN/bY27NuBgnHPtnHNznHMbnHN7nHPfOecmO+dOjru2bFc+7gIg4py7QUTGFv4aIYk/pLQUkcpx1gVE5ZzrLiIt4q4DKIZbRGTBfn+fH1chQEQ1RORfIvKkiPwgIvVFZJCIfOKcO9V7/22cxWUzmuGYOecaishoEbnTez96v3/0TiwFAcXknKsuIqNEZKCITIq5HCCqZd77T+IuAojKe/+SiLy0f+ac+0xElotINxH5axx1lQV8TCJ+vUSkQESeirsQ4BD9WUSWFt6oAQCZs7nw97xYq8hyNMPxy5XEn+qudM5945zLd86tcs4NiLswoCjOuVwR+YOI3Bh3LUAxTXTO7XPObXbOTXLO1Y+7ICAK51w551yOc+5EEfmbiGwQkZdjLiur8TGJ+NUt/PUXERkiIt+IyBUiMtY5V957PybO4oBknHMVJHEjfsR7/3Xc9QAR/SyJHyd/ICLbRKSVJO69/+2ca+W93xRncUAEn4pI68K/XiUiF3HdpsZ57+OuIWjOuRUicqKIXO69n75f/pYkbtJ1PP+TUAo554ZJ4mM+zb33uwszLyIPeu+HxVocUAzOudNF5DMReYhrF6Wdc+6/ROQIEWkkIneISC0RyfXer46zrmzGxyTi9z+f9/nnAfm7krjA62S2HKBohT9SHioid4vI4c656oWDdLLf35eLr0IgOu/9QhFZISJnxl0LUBTv/TLv/aeFcxoXi0hVSbxVAoeIZjh+S5PkrvD3gkwVAhRDIxGpKCIvishP+/0SSTyp+ElETo2nNOCQOBHhp3DIKt77rZL4qESTuGvJZjTD8ZtR+Hu7A/J2IvKd935DhusBovhCRC40fokkGuQLJXGDBko959wZInKSJD6LCWQN51wtEWkmiXkjHCIG6OL3pojMFZG/OeeOEZF/S+J9gZeISM84CwOSKXwa8f6BuXNORORb7736Z0Bp4JybKCL/EZGFIrJVErMZg0XkexF5PMbSgINyzs2QxHW7WBLDnydJ4v3u+cI7hlNCMxwz7713zl0mIiNF5D4ROUoSr1q72nvPAgYAkF5LRKS7iNwsiVU+N4jIdBG513v/Y5yFAUX4RER+JyK3i0iOiKyVxEOJkQzPpYa3SQAAACBYfGYYAAAAwaIZBgAAQLBohgEAABAsmmEAAAAEK6NvkyhcqhVIiffeFb1VenHtIh0yfe1y3SIduOciW0W9dnkyDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIVvm4CwhFhw4dVDZw4ECVtW3bVmXee5WtXLnSPM/kyZNVNm7cOJWtW7fO3B8AACAkPBkGAABAsGiGAQAAECyaYQAAAASLZhgAAADBctZwVomdzLnMnSwm/fv3N/NRo0apLCcnp6TLERGRuXPnqqxHjx4qW79+fSbKSZn33mX6nCFcuyh5mb52uW6RDtxzi+e5555T2TXXXKOyWbNmmftPmzZNZfPnz1fZ2rVrI9Wzd+9eM9+3b1+k/bNZ1GuXJ8MAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYrECXgo4dO6rskUceMbe1huUWLVqkskGDBqls6dKlkWvq3bu3yu677z6VDR48WGW33HJL5PMgu1WpUkVlQ4YMMbcdNmyYyqzB2xEjRqisRYsWKuvcuXOUEgEgKy1fvlxlBQUFKrN6iIPlh2rChAlmfsMNN6gsPz8/refOFjwZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwWIFuog6deqkspdeekll1mCSiMjMmTNVZq1Wt3HjxkOo7n85pxdbsYbqLrnkEpX97ne/S+ncmcJqSKmrX7++yr799ltz29atW6ts4cKFKrMG6G6++WaVNW3a1DxPqtd+NmAFOuyvVq1aKmvSpIm5bcWKFVXWvXt3lU2cOFFlyVYg+/jjj4sqUUS456aD1UO0a9cu8v5nnnmmyqz7eKVKlVR25JFHmse8+OKLVWatWJvNWIEOAAAAKALNMAAAAIJFMwwAAIBg0QwDAAAgWKxAZyhfXv9nsVZxs4blFi9ebB7TWunlhx9+OITqDs4aiBw/frzKZsyYkfZzI3s0bNgw7cfMy8tTmTW4cfLJJ5v7hzBAhzCccsopKvv973+vsl69eqmsTp065jGjDrv37Nkz0nYiIuXKlYu8LVLzxhtvRMpS1aFDB5XNmjXL3PbSSy9VWVkboIuKJ8MAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFm+TMPTt21dlrVq1UtmePXtUdt1115nHLIk3R6Ri8+bNcZeAGJ177rlpP+arr76qMustLGeccYa5f6hTzMgOLVu2NPOBAweqrE2bNiqrXbt22muybN++XWVz5szJyLmRWTVq1FDZvffeq7L8/Hxz/2RvmQgRT4YBAAAQLJphAAAABItmGAAAAMGiGQYAAECwGKAz3HzzzZG269evn8q++OKLdJcDpMRacvXyyy9XWUFBgbl/suELoLispe5FRCpWrKiyHTt2lHQ5ImIPdE6YMEFljRs3Nvc//PDD016T5auvvlLZsGHDVGYNR8+bN69EakJqqlWrZua5ubkqy8nJUdnQoUNVZl3Pzz//vHme999/v4gKw8GTYQAAAASLZhgAAADBohkGAABAsGiGAQAAECwG6FLw3XffxV0CUKRatWqp7Mwzz1TZf/7zH3P/xYsXRzpPXl6eyvbt26eyJk2aRDoeyh5rdSwRkcsuu0xl06ZNU9nw4cMjn+u0005T2V133aUya5i0QoUKKnPOmefx3keuKQrr31tE5A9/+IPKdu/endZzIz2qVq2qspEjR6rMuvZEUlut8NNPP1XZQw89dMjHCwVPhgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLCCHqCzBixERE488USVbd++XWVff/112msC4rJy5cqU9l+1apXK1q5dq7KWLVumdB5khyOOOEJl11xzjblt/fr1Vda8eXOVWYNJTZs2NY/ZsWPHokoslmQDdBZrFbgXXnhBZdOnT1cZq8Vlv/PPP19lAwYMyMi5ra+RZKuL4n/xZBgAAADBohkGAABAsGiGAQAAECyaYQAAAASLZhgAAADBCvptEuXL2//65cqVU9muXbtUxnLMyAYXXXRRpO1GjRqV0nmsryfra6lOnTrm/tbbB7Zt25ZSTYhPjRo1VFalShVz26hLGg8cOFBlJbFM8oIFC1T2yiuvmNu++eabKtuxY4fKvv/++0OuB9klNzc3pf03bdqksnHjxqnssMP088y7775bZdZS0CIiffr0UdlPP/0UpcQyhyfDAAAACBbNMAAAAIJFMwwAAIBg0QwDAAAgWEEP0MXt6KOPVlmnTp1Udvvtt0c+5urVq1XWsGFDlW3YsEFlU6dOVdmECRPM8+Tl5UWuCfE677zzVLZx40aVffTRRymdxxoynTVrlsr69etn7n/kkUeqjAG67GXdi3744QdzW2vYLlNGjBihsscee0xlW7ZsyUQ5KAPuu+8+lf3rX/9S2c6dO839P/jgA5Xt3btXZdbw6JQpU1T23nvvmecZP368ynr37q2yrVu3mvuXJTwZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwWKALiJrwOOMM85Q2eeff27u36RJE5XNnj1bZfXr11fZ7t27Vfbll1+a57GGVqysZ8+eKmvTpo3K2rVrZ57n8ssvN3PEy1rh69JLL1WZNYyRbJgjFSEMXiAisjLlAAAIkklEQVS6ZIM8TZs2PeRjfvjhh2Y+bdo0lU2aNEll1opbBQUFh1wPkJ+fr7KZM2em/TzWKotLlixRWd++fc39Z8yYobK5c+eqbOzYsYdQXXbhyTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAhW0AN0yVYU+vnnn1VmrY5lZY0aNTKPOWfOHJUdd9xxKrMGTAYMGKCyFStWmOeJ6rXXXlOZ9WH6Zs2apXQeZFblypVV1qBBA5WtXbs2E+WYX0vJWF9PmaoTmTF48GAzt1betIaJLRdccEEqJQFlnvX9XkTk5ZdfVpn1NfrKK6+oLNlqktmKJ8MAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYQQ/QWSuziYisX79eZdZwz1VXXaWyk08+2TymNSxnrUDXtWtXlZXEymDWucePH6+ySy65JO3nRvxycnJU1rp1a3PbX375RWXW8GmlSpVUZq2QlMy4ceNUdtFFF6ksLy8v8jFRuuzYscPMrUGeHj16qKxevXoq27Bhg3nMKVOmqOzee+9VWbJBaqCsGzNmjMq6d++usuuvv15lDz74YInUFBeeDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYrjjT3imfzLnMnSwFI0eOVNldd92V0jGtNzXceuutKtu1a1dK50nFpEmTVNa+fXtz25YtW6pszZo1aa/J4r13GTnRfrLl2j322GNVtmnTppSOmZ+frzLrrQDWGyqs5aGLw3q7ysyZM1M6Zpwyfe1my3Vrsaban3rqKZVVq1bN3N/63jZ//nyVde7cWWU//fRTlBKDwT23bKpYsaLKPv74Y5UtXrxYZT179iyRmtIt6rXLk2EAAAAEi2YYAAAAwaIZBgAAQLBohgEAABAsBugM1atXV9kXX3yhsvr160c+5m233aay0aNHF6+wEmYta5psOOX0009X2ddff532miwMcyRXrlw5lY0YMUJlgwcPzkQ5xfL555+r7JxzzlHZvn37MlFOiWCALjXWPdcaThYRufjiiyMd86uvvlLZFVdcobLly5dHOl5ZxD3XXgpcxB7q7Natm8r27NmT9ppKwrBhw1R2ww03qOzUU09V2datW0ukplQwQAcAAAAUgWYYAAAAwaIZBgAAQLBohgEAABAsBugi6tixo8pefvlllVWpUsXcf+fOnSp74403VPbggw+qbMmSJVFKLJYOHTqo7LXXXlPZihUrzP2bN2+e9pqiYpijeKyhupo1a6os2bVrXSvW0JGVWUMW77zzjnkea3Ww888/39w2WzFAl37WkKWIvVKhtUKjZcGCBSq76aabzG2twc+yhnuuSMOGDc383//+t8peeOEFlf3pT39S2caNG1OuK92sAbr7779fZY0aNVLZ6tWrS6KklDBABwAAABSBZhgAAADBohkGAABAsGiGAQAAECwG6FLQrl07lT388MPmtqeddlqkY+7evVtlffr0UdmaNWvM/a0PsOfm5qpszJgxKrNW3nvppZfM8/Ts2dPMM4FhjuzRunVrlSUbOGKALv1Cvm4vu+wylU2bNu2Qj2fdh0VEJkyYcMjHzBbcc0Xq1q1r5tbKq9Yw8sqVK1XWr18/85gfffSRyvLz84sqsdi6du2qskceeURlOTk5KjvllFNU9vPPP6ensDRigA4AAAAoAs0wAAAAgkUzDAAAgGDRDAMAACBYDNClWbIVjnr16qUya0Wao446Ku01WawP41ur3913332ZKKdYGObIHsccc4zKli9fbm67b98+lZ100kkqK41DGlExQJd+/fv3N/Mnnngired59tlnzdy6t5c13HOT69atm8omT56c0jGtlemsXu3VV19VWZcuXSKfp0aNGiqzhuUeeOABld1zzz2RzxMnBugAAACAItAMAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYPE2iRhZk5zWZLQ1rdqiRYvI51m7dq3KnnrqKZWNHDky8jHjxGRzdrOWXRYROffcc1VmLYG6fv36tNeUKbxNIjprufvBgwer7Ne//rW5f7q/t910001mPm7cuLSepzTinptcuXLlVNa+fXuVDRo0SGWpLjfvnP7fkup1P378eJUNHTpUZT/88ENK58kU3iYBAAAAFIFmGAAAAMGiGQYAAECwaIYBAAAQLAbokHUY5shuAwcONPNHH31UZZdddpnKrCVIs0XoA3QdOnQw8+uvv15l1hCStVSsNUQkEn2QaMSIESpbuHChyl577bVIxyuLuOem7rDD9LPHs846y9zWGpo/77zzVHbOOeeobO/evSqbMmWKeZ4xY8aozLr2CwoKzP2zAQN0AAAAQBFohgEAABAsmmEAAAAEi2YYAAAAwWKADlmHYY7sdvbZZ5v5J598orL3339fZRdeeGG6S8qYkAbo+vTpo7Jkq1xaq3Fatm7dqrJ58+aZ23755Zcqmz59usoWL16ssmweGCoJ3HORrRigAwAAAIpAMwwAAIBg0QwDAAAgWDTDAAAACBYDdMg6DHMgW4U0QGetmNWxY0dz21mzZkU65qZNm1S2atWq4hWGYuOei2zFAB0AAABQBJphAAAABItmGAAAAMGiGQYAAECwaIYBAAAQLN4mgazDZDOyVUhvk0DZwT0X2Yq3SQAAAABFoBkGAABAsGiGAQAAECyaYQAAAASLZhgAAADBohkGAABAsGiGAQAAECyaYQAAAASLZhgAAADByugKdAAAAEBpwpNhAAAABItmGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABItmGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABItmGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABItmGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABItmGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABItmGAAAAMGiGQYAAECw/h+g8MKF0kXt3QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plots(x_imgs[:8], titles=y_valid[:8])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural Networks"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "We will take a deep look logistic regression and how we can program it ourselves. We are going to treat it as a specific example of a shallow neural net.\n",
    "\n",
    "What is a neural network?\n",
    "\n",
    "A neural network is an infinitely flexible function, consisting of layers. A layer is a linear function such as matrix multiplication followed by a non-linear function (the activation).\n",
    "\n",
    "One of the tricky parts of neural networks is just keeping track of all the vocabulary!\n",
    "\n",
    "Functions, parameters, and training\n",
    "A function takes inputs and returns outputs. For instance, $f(x) = 3x + 5$ is an example of a function. If we input $2$, the output is $3\\times 2 + 5 = 11$, or if we input $-1$, the output is $3\\times -1 + 5 = 2$\n",
    "\n",
    "Functions have parameters. The above function $f$ is $ax + b$, with parameters a and b set to $a=3$ and $b=5$.\n",
    "\n",
    "Machine learning is often about learning the best values for those parameters. For instance, suppose we have the data points on the chart below. What values should we choose for $a$ and $b$?"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "In the above gif from fast.ai's deep learning course, \n",
    "\n",
    "intro to SGD notebook - https://github.com/fastai/courses/blob/master/deeplearning1/nbs/sgd-intro.ipynb )\n",
    "\n",
    ", an algorithm called stochastic gradient descent is being used to learn the best parameters to fit the line to the data (note: in the gif, the algorithm is stopping before the absolute best parameters are found). This process is called training or fitting.\n",
    "\n",
    "Most datasets will not be well-represented by a line. We could use a more complicated function, such as $g(x) = ax^2 + bx + c + \\sin d$. Now we have 4 parameters to learn: $a$, $b$, $c$, and $d$. This function is more flexible than $f(x) = ax + b$ and will be able to accurately model more datasets.\n",
    "\n",
    "Neural networks take this to an extreme, and are infinitely flexible. They often have thousands, or even hundreds of thousands of parameters. However the core idea is the same as above. The neural network is a function, and we will learn the best parameters for modeling our data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyTorch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be using the open source deep learning library, fastai, which provides high level abstractions and best practices on top of PyTorch. This is the highest level, simplest way to get started with deep learning. Please note that fastai requires Python 3 to function. It is currently in pre-alpha, so items may move around and more documentation will be added in the future.\n",
    "\n",
    "The fastai deep learning library uses PyTorch, a Python framework for dynamic neural networks with GPU acceleration, which was released by Facebook's AI team.\n",
    "\n",
    "PyTorch has two overlapping, yet distinct, purposes. As described in the PyTorch documentation:\n",
    "\n",
    "It’s a Python based scientific computing package targeted at two sets of audiences:\n",
    "\n",
    "A replacement for NumPy to use the power of GPUs\n",
    "B a deep learning research platform that provides maximum flexibility and speed\n",
    "\n",
    "The neural network functionality of PyTorch is built on top of the Numpy-like functionality for fast matrix computations on a GPU. Although the neural network purpose receives way more attention, both are very useful. We'll implement a neural net from scratch today using PyTorch.\n",
    "\n",
    "Further learning: If you are curious to learn what dynamic neural networks are, you may want to watch this talk by Soumith Chintala, Facebook AI researcher and core PyTorch contributor.\n",
    "\n",
    "If you want to learn more PyTorch, you can try this introductory tutorial or this tutorial to learn by examples.\n",
    "\n",
    "https://pytorch.org/tutorials/beginner/deep_learning_60min_blitz.html\n",
    "\n",
    "https://pytorch.org/tutorials/beginner/pytorch_with_examples.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### About GPUs"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Graphical processing units (GPUs) allow for matrix computations to be done with much greater speed, as long as you have a library such as PyTorch that takes advantage of them. Advances in GPU technology in the last 10-20 years have been a key part of why neural networks are proving so much more powerful now than they did a few decades ago.\n",
    "\n",
    "You may own a computer that has a GPU which can be used. For the many people that either don't have a GPU (or have a GPU which can't be easily accessed by Python), there are a few differnt options:\n",
    "\n",
    "Don't use a GPU: For the sake of this tutorial, you don't have to use a GPU, although some computations will be slower.\n",
    "\n",
    "Use crestle, through your browser: Crestle is a service that gives you an already set up cloud service with all the popular scientific and deep learning frameworks already pre-installed and configured to run on a GPU in the cloud. It is easily accessed through your browser. New users get 10 hours and 1 GB of storage for free. After this, GPU usage is 34 cents per hour. I recommend this option to those who are new to AWS or new to using the console.\n",
    "\n",
    "Set up an AWS instance through your console: You can create an AWS instance with a GPU by following the steps in this fast.ai setup lesson.] AWS charges 90 cents per hour for this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural Net for Logistic Regression in PyTorch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fastai.metrics import *\n",
    "from fastai.model import *\n",
    "from fastai.dataset import *\n",
    "\n",
    "import torch.nn as nn"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "We will begin with the highest level abstraction: using a neural net defined by PyTorch's Sequential class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "net = nn.Sequential(\n",
    "    nn.Linear(28*28, 100),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(100, 100),\n",
    "    nn.ReLU(),\n",
    "    nn.Linear(100, 10),\n",
    "    nn.LogSoftmax()\n",
    ")#.cuda()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Each input is a vector of size 28*28 pixels and our output is of size 10 (since there are 10 digits: 0, 1, ..., 9).\n",
    "\n",
    "We use the output of the final layer to generate our predictions. Often for classification problems (like MNIST digit classification), the final layer has the same number of outputs as there are classes. In that case, this is 10: one for each digit from 0 to 9. These can be converted to comparative probabilities. For instance, it may be determined that a particular hand-written image is 80% likely to be a 4, 18% likely to be a 9, and 2% likely to be a 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "md = ImageClassifierData.from_arrays(path, (x,y), (x_valid, y_valid))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "#?? ImageClassifierData"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss=nn.NLLLoss()\n",
    "metrics=[accuracy]\n",
    "# opt=optim.SGD(net.parameters(), 1e-1, momentum=0.9)\n",
    "opt=optim.SGD(net.parameters(), 1e-1, momentum=0.9, weight_decay=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loss functions and metrics"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "In machine learning the loss function or cost function is representing the price paid for inaccuracy of predictions.\n",
    "\n",
    "The loss associated with one example in binary classification is given by: -(y * log(p) + (1-y) * log (1-p)) where y is the true label of x and p is the probability predicted by our model that the label is 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binary_loss(y, p):\n",
    "    return np.mean(-(y * np.log(p) + (1-y)*np.log(1-p)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.164252033486018"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acts = np.array([1, 0, 0, 1])\n",
    "preds = np.array([0.9, 0.1, 0.2, 0.8])\n",
    "binary_loss(acts, preds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in our toy example above our accuracy is 100% and our loss is 0.16. Compare that to a loss of 0.03 that we are getting while predicting cats and dogs. Exercise: play with preds to get a lower loss for this example.\n",
    "\n",
    "Example: Here is an example on how to compute the loss for one example of binary classification problem. Suppose for an image x with label 1 and your model gives it a prediction of 0.9. For this case the loss should be small because our model is predicting a label  $1$ with high probability.\n",
    "\n",
    "loss = -log(0.9) = 0.10\n",
    "\n",
    "Now suppose x has label 0 but our model is predicting 0.9. In this case our loss is should be much larger.\n",
    "\n",
    "loss = -log(1-0.9) = 2.30\n",
    "\n",
    "Exercise: look at the other cases and convince yourself that this make sense.\n",
    "Exercise: how would you rewrite binary_loss using if instead of * and +?\n",
    "Why not just maximize accuracy? The binary classification loss is an easier function to optimize.\n",
    "\n",
    "For multi-class classification, we use negative log liklihood (also known as categorical cross entropy) which is exactly the same thing, but summed up over all classes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting the model"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Fitting is the process by which the neural net learns the best parameters for the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4a69baa08d7f48098e8e32d449aa36d0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Epoch', max=5), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch      trn_loss   val_loss   accuracy                     \n",
      "    0      0.286935   0.257693   0.9278    \n",
      "    1      0.23487    0.181463   0.9514                       \n",
      "    2      0.183344   0.169012   0.9534                       \n",
      "    3      0.196265   0.171689   0.9488                       \n",
      "    4      0.169468   0.158852   0.9542                       \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[array([0.15885]), 0.9542]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit(net, md, n_epochs=5, crit=loss, opt=opt, metrics=metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_lrs(opt, 1e-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9a0d21cce6504b0ab0b50ab3059a3ae1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Epoch', max=3), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch      trn_loss   val_loss   accuracy                      \n",
      "    0      0.06382    0.086563   0.9749    \n",
      "    1      0.058824   0.081585   0.9756                        \n",
      "    2      0.048132   0.076393   0.9763                        \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[array([0.07639]), 0.9763]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit(net, md, n_epochs=3, crit=loss, opt=opt, metrics=metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1abd2e818b0b4774a532fc3bf4869300",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Epoch', max=5), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch      trn_loss   val_loss   accuracy                      \n",
      "    0      0.050011   0.075903   0.976     \n",
      "    1      0.033426   0.074766   0.9769                        \n",
      "    2      0.034524   0.072397   0.9783                        \n",
      "    3      0.030426   0.068434   0.9786                        \n",
      "    4      0.030118   0.070237   0.9789                        \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[array([0.07024]), 0.9789]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit(net, md, n_epochs=5, crit=loss, opt=opt, metrics=metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_lrs(opt, 1e-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2286df8b93264e2cb18d506dcc98a909",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Epoch', max=3), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch      trn_loss   val_loss   accuracy                      \n",
      "    0      0.031436   0.068336   0.9791    \n",
      "    1      0.026236   0.068711   0.9791                        \n",
      "    2      0.022527   0.070696   0.9787                        \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[array([0.0707]), 0.9787]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit(net, md, n_epochs=3, crit=loss, opt=opt, metrics=metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([78400, 100, 10000, 100, 1000, 10], 89610)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = [o.numel() for o in net.parameters()]\n",
    "t, sum(t)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "GPUs are great at handling lots of data at once (otherwise don't get performance benefit). We break the data up into batches, and that specifies how many samples from our dataset we want to send to the GPU at a time. The fastai library defaults to a batch size of 64. On each iteration of the training loop, the error on 1 batch of data will be calculated, and the optimizer will update the parameters based on that.\n",
    "\n",
    "An epoch is completed once each data sample has been used once in the training loop.\n",
    "\n",
    "Now that we have the parameters for our model, we can make predictions on our validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "preds = predict(net, md.val_dl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000, 10)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds.shape"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Question: Why does our output have length 10 (for each image)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 8, 6, 9, 6])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds.argmax(axis=1)[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-17.86899, -14.92701, -13.58845, ..., -14.86605, -11.01833, -12.49672],\n",
       "       [-13.82064, -15.14436, -11.94373, ..., -15.64362,  -0.00008, -14.28446],\n",
       "       [ -8.54279, -11.40613, -10.33679, ..., -12.00947, -11.27964, -11.14462],\n",
       "       ...,\n",
       "       [-16.57417, -15.06741, -19.47988, ..., -16.7648 , -14.50415, -11.59032],\n",
       "       [ -7.01595,  -8.94965,  -5.93617, ...,  -8.96127,  -7.05934,  -9.48682],\n",
       "       [ -7.43876,  -8.94534,  -7.01448, ...,  -8.94267,  -0.00304,  -7.07597]], dtype=float32)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "preds = preds.argmax(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check how accurate this approach is on our validation set. You may want to compare this against other implementations of logistic regression, such as the one in sklearn. In our testing, this simple pytorch version is faster and more accurate for this problem!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9787"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(preds == y_valid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how some of our predictions look!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAF0CAYAAADGqzQSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmUVeWV9/G9BSqMimhAwAABDLSogGicKokTAsoCaTEt0XQEQTEYI1FbQJyCvpjECCiKSWvQKEZBBAc0GgQHpFUMKu2AggmCMimIDCJUUc/7x6Xf9nXvS53i3rqnbj3fz1ou4m+dc55tPHXZHu4+j4YQBAAAAIjRPmkXAAAAAKSFZhgAAADRohkGAABAtGiGAQAAEC2aYQAAAESLZhgAAADRohlOmar2UtV5qrpWVXeo6seqOl1VD027NmBPVPUEVX1WVder6mZVXayqQ9KuC0hCVU9X1RdVdevu+/d1VT057bqAPVHVk1R1gapuV9WNqnq/qrZIu65iRzOcvmYi8ncRuUREThOR0SLSRUReUdW2aRYGZKOqR4jIXBGpJyLDROQsEVkkIveo6sVp1gZURlUvEpHHJPPZO0BEzhaRGSLSMM26gD1R1R+IyLMiskkyn7m/FJEfishzqvqtNGsrdsqmGzWPqnYSkaUickUI4fdp1wN8k6r+HxG5QkSahRC2fi1/RURCCOG41IoD9kBV24nIeyIyOoQwMd1qgORUda6ItBORziGE8t3Z0SLymoiMCCHcmWJ5RY0nwzXTht2/lqVaBZBdiWTuz+3fyDcJnyuo2YaISIWI3JV2IUAVHSsif/ufRlhEJISwSDI9w4DUqqoF+E2rhlDVOqpaoqqHiMgfRGStiDyUcllANvfu/vU2VW2lqk1VdZiInCIiE9IrC6hUqWT+5O0cVf1QVctVdbmqjki7MKASu0Rkp5PvEJHDClxLrcLXJGoIVX1dRHrs/tvlItIvhPBeiiUBe7T7j+dmiUjr3VGZiFwcQrgnvaqAPVPVpSLSSjINxBgR+VAy3xkeLiKXhRAmpVgekJWqviaZr6Ed87WsrYj8U0TKQgh8b3gv0QzXEKr6LyKyr4i0l8x3MVuISGkIYUWadQGe3X+C8Zxkvnt5u2S+LtFfRC4WkfNDCNNSLA/ISlU/EJFDROSsEMKjX8ufFpHuItIy8BsjaiBVPVdEHhCRm0TkNskM4P9RRI6XTDPcIMXyihrNcA2kqk1FZIWIPBRCGJ5yOYChqjNE5EjJDHKUfS2fJiK9RKR5CKEirfqAbFT1vyTz3ct9QwhbvpaPFJFbRaR1CGF1WvUBe6Kq4yTzwKy+iAQReVhEGonIYSGE9mnWVsz4znANFELYJJmvSnRMuxYgi8NF5K2vN8K7vSYiB4hI88KXBCTyTpZcd//Kf8ShxgohXCMiB4rIEZL5U4xBkvmTjgWpFlbkaIZroN0v0O4sme+yATXRWhHppqol38iPEZGvRGRj4UsCEpm1+9de38h7icjHIYS1Ba4HqJIQwrYQwn+HENapam/J9Au8HSUHddMuIHaqOktEFovIEhHZLCLfE5GRIlIuIrxjGDXVZMlsUvCEqt4pme8M9xORQSIyIYTgTTwDNcFTIjJfRP6gqgeKyD9EZKBkNj0anGZhwJ6oancR6SOZnkEk82aUK0XktyGEhakVVgvwneGUqepVIvJjEekgmXe3rhKR50VkPMNzqMlUtY+IXCWZHRPrS+ZPMv4oIn8IIexKszZgT1R1XxEZL5kmeH/JvGrt5hDCg6kWBuyBqnaRzKtXDxORb8nuAeYQwtRUC6sFaIYBAAAQLb4zDAAAgGjRDAMAACBaNMMAAACIFs0wAAAAolXQV6upKtN6yFkIQSs/Kr+4d5EPhb53uW+RD3zmolglvXd5MgwAAIBo0QwDAAAgWjTDAAAAiBbNMAAAAKJFMwwAAIBo0QwDAAAgWjTDAAAAiBbNMAAAAKJFMwwAAIBoFXQHOgAAapLGjRub7IILLjBZ//793fP79etnsq1bt+ZeGICC4ckwAAAAokUzDAAAgGjRDAMAACBaNMMAAACIFs0wAAAAosXbJAAA0frZz35msgkTJiQ+v0uXLiZ79dVXc6oJQGHxZBgAAADRohkGAABAtGiGAQAAEC2aYQAAAESLAbocdO3a1WQjR450j+3QoYPJGjZsaLIxY8aYbL/99jPZ008/7a6zZcsWNweA2J1//vkmmzhxosnKyspMdsstt7jXXLx4cc51AUgXT4YBAAAQLZphAAAARItmGAAAANGiGQYAAEC0NIRQuMVUC7dYnjVu3NhkK1euNFnTpk0LUY588sknbu4N8D3yyCPVXU5BhRC00GsW873r8e7TAQMGuMd2797dZKWlpSbzfkY2btxosoMOOshdZ+3atSa79957Tfaf//mfJtu1a5d7zZqm0Pdubbtvq6Jfv34mmzVrlsm+/PJLk1177bUmq8qudLUNn7koVknvXZ4MAwAAIFo0wwAAAIgWzTAAAACiRTMMAACAaDFAl1CTJk1M9tRTT5lsw4YN7vlvvPGGybzBpLZt25rsO9/5jskaNGjgrrNu3TqTHXfccYmOKxYMc1TNwQcfbLLZs2ebzLsfs9m8ebPJvHu8Xr16JvN+lkREmjdvbrIWLVqY7Cc/+YnJXnzxRZOtWbPGXSdNDNDlX0lJiZtPnTrVZIMGDTLZvHnzTHbqqafmXlgtwmcuihUDdAAAAEAlaIYBAAAQLZphAAAARItmGAAAANFigK4IHHjggSa78sor3WO9fPDgwSa77777ci8sJQxzVM3ixYtN1rVrV5PNnTvXPf/yyy832WeffWYybwe5qvj2t79tsqefftpknTp1MtmoUaNMdscdd+RUT3VggC7/rr76ajcfN26cyR544AGTDRkyxGTl5eW5F1aL8Jmbu5YtW5rs5z//uXusl5eVlZnM2wX3pptuMpn3e4CIyKpVq9y8NmGADgAAAKgEzTAAAACiRTMMAACAaNEMAwAAIFo0wwAAAIgWb5MoUv369XNzb5vd2267zWSXXXZZ3msqFCabs/Mmlj/55BOTTZ8+3WTnnnuue81du3blXthemjZtmsnOOecck/Xo0cNkb775ZrXUlAveJpGbo446ymQLFixwj12xYoXJunTpYrI07+9iwWdu1bRv395kU6ZMMVnPnj0LUY7s2LHDzU844QSTZXvzRLHibRIAAABAJWiGAQAAEC2aYQAAAESLZhgAAADRqpt2Aajc/vvvb7IxY8YkPr9Vq1b5LAc1WLdu3UymaucHVq9ebbK0B4mOPfZYkw0aNMhk8+fPN5n3z10TB+iQ3D772Gc13rbbJSUl7vlPPPGEydK+x1H7tG7d2mRvv/22yerWte3WhAkT3Gvefvvtidbp3LmzyX73u9+ZrGnTpu463iC19zn82WefuefXJjwZBgAAQLRohgEAABAtmmEAAABEi2YYAAAA0WIHuhqma9euJpsxY4bJOnbs6J7/wQcfmMzb5WbVqlV7UV3NwG5IVVNRUWGy9evXm+z73/++e/7KlSvzWk+TJk3cfOHChSZbtmyZybyd8rwdn9555529qK56sQNdckl3U8zm0ksvNdnkyZNzqilWfOZmN2nSJJMNHz7cZMOGDTPZn//857zXM2LECJNNnDjRPbZOnTomW7p0qcm8obrNmzfvRXWFxw50AAAAQCVohgEAABAtmmEAAABEi2YYAAAA0WKALkU/+9nPTPbrX//aZN/5zndMtn37dveaffv2NZm3Y1cxY5ijaq6//nqTXXPNNSZ7//333fN79eplslwGMJ999lk3/9GPfmSyHj16mMzb3alYMECX3ODBg012zz33mGzu3Lnu+X369DEZO9DtHT5zRfbdd18394Z8p06dajJv98RCyfbZfsghhyQ639sp7/LLL8+ppkJhgA4AAACoBM0wAAAAokUzDAAAgGjRDAMAACBaNMMAAACIFm+TyLPGjRu7+RVXXGGysWPHmmyffex/n2zcuNFkpaWl7jreVoq1DZPNVVO/fn2T3XfffSYbOHCge/7y5ctNduKJJ5pszZo1JrvzzjtNduGFF7rrXHnllSbzppiLGW+T8NWtW9dk7733nsnatm1rsu9+97vuNauydTP2jM/c7NvVv/LKKybr2bOnyZ577rm815TUgAED3PzRRx81mdcTbtq0yWTemyg2bNiwF9VVL94mAQAAAFSCZhgAAADRohkGAABAtGiGAQAAEC07tYCc3HvvvW7+r//6r4nOf+SRR0w2ceJEk8UwKIf8+Oqrr0w2dOhQkzVv3tw939sm+YUXXjDZjBkzTHbeeeeZbObMme46tW1YDsl5w5sdOnQw2cUXX2yytAflevfubbJ+/fqZ7K9//avJvK3JvZ9XpK979+6Jj33jjTeqsZKqe+qpp9zcG472fu68e3Lbtm25F1aD8GQYAAAA0aIZBgAAQLRohgEAABAtmmEAAABEiwG6PPO+fF4VU6ZMMdnChQtzuibwTVu2bDFZ//793WOvv/56k1122WUmGzVqVKK1b7/99kTHIR5t2rRJdFxJSUk1V5Ld+eef7+beLovero/Dhw83mbez1+zZs911hgwZUkmFqE4LFixw84qKCpP97W9/M1nfvn1N5u3aWR06derk5t592qtXL5M1bNjQZLVt0JMnwwAAAIgWzTAAAACiRTMMAACAaNEMAwAAIFoM0OWZt6OQiEjXrl33+nxvqO7mm292z1+9enWidYBv2rx5s5tfe+21JuvZs6fJDj300ETrnHrqqW6ebUAFtV/Hjh0THVeonTebNm1qsltvvdU91htCKi8vN5k3VFVaWmoyb9dGEQbo0vbOO++4+ZNPPmkybxj5vffeM5m3K6GIv0vnvHnzTNa6dWuTecNy3i62IiItW7Y0mXfvPvbYY+75tQlPhgEAABAtmmEAAABEi2YYAAAA0aIZBgAAQLQ0hFC4xVQLt1hKGjRo4OYPPPCAyXr06GGypDsxrV271s0HDx5ssmeeeSbRNYtFCEELvWYM9242ffr0MdmsWbNMVq9evUTX27lzp5v//Oc/N9nUqVMTXbNYFPreLZb7ds6cOSbr3r27yVq1alWIctwdFrMN0Hmf7ZMmTTLZypUrTeYNUB1++OHuOmnuvsdnbnbe7/njx4832aWXXprTOhs3bjRZs2bNcrqm5+yzzzaZN9BXLJLeuzwZBgAAQLRohgEAABAtmmEAAABEi2YYAAAA0aIZBgAAQLTYjjnPtm/f7ubnnnuuyerWtf/3Z9sS95sOOuggN/em/H/1q1+Z7K677kq0DnDSSSeZzHsLzYABA0zmTUB725eK+NuOf/bZZyZ74okn3PNRvI455hiTZXvrSE2zevVqkx188MEm++Mf/2iyI4880mS17e0/tZ33e773NpLp06ebzOsLsmnRokWi48rKykzm/XyJiHz3u9812Zdffpm4ptqEJ8MAAACIFs0wAAAAokUzDAAAgGjRDAMAACBabMdcwxxxxBEmmzBhgsm8oaZsvG1A27VrV6W6ahK2Bq0e3r0nIrJo0SKTecNu3tCIx9vuU0TknnvuMZmq/VfdpUsXk3n3eE3Edsw+b7isb9++JquO7Zi9e8y7l3//+9/ntI73e+2dd95psjFjxrjnb9myJaf1c8FnbnG7//773dwb4Ovdu7fJnn322bzXVChsxwwAAABUgmYYAAAA0aIZBgAAQLRohgEAABAtdqBLqGHDhiarjp1alixZYrKBAwea7E9/+pN7fv/+/U3Wpk0bk7Vs2dJka9asSVIiaqkmTZq4ubdT4iOPPLLX68yYMcPN27Zta7Lf/OY3JuvRo4fJimWADsk1bdrUZN4g0AMPPOCe792355xzjsmaNWtmsj59+iQpUUREtm3bZrIFCxaY7Le//a3J5s+fn3gdoBA6dOiQdgmp4MkwAAAAokUzDAAAgGjRDAMAACBaNMMAAACIFgN0Du8L5N5AxJw5c0z29ttvu9f0htMuuOACk9WrV89krVu3NlnHjh3ddTwffvhhonoQt27durn52rVrTeb9PORq8uTJJhs2bJjJRowYYbJZs2blvR4UzhtvvGGyoUOHmszbMcvLcrV582aTZRv8vPHGG0320Ucf5b0mYG9t3bo17RJqPJ4MAwAAIFo0wwAAAIgWzTAAAACiRTMMAACAaDFA5zj77LNNdtBBB5lsyJAheV9bVU0WQkh8vvdF+eHDh+dUE+Lg7VQoIvLaa68VZP2dO3ea7PPPPzfZD37wA5N5u4ht3LgxP4Wh2j344IMm83beXLZsmcnq1KnjXjNb/k3Tpk0z2YoVK0zmDSIDxeDFF19084suushkzZs3r+5yaiSeDAMAACBaNMMAAACIFs0wAAAAokUzDAAAgGjRDAMAACBavE3CccABB6Rdwv9n5syZJhs3bpx77Pr1603mbacLfFO2t5aUlpaa7JxzzjHZvHnzTNa4cWOTlZSUuOt07tzZZEcffbTJ7rjjDpPx5oji9sUXX5jslFNOSaESoPbZZx//uaf39iqvh4gBT4YBAAAQLZphAAAARItmGAAAANGiGQYAAEC0GKBzjBkzxmRz58412XnnnWeyVq1audf0BkQ8t99+u8leeuklk5WXlye6HpDUe++95+beVsfe9rkbNmwwWVUG6Lxhjpdfftlk119/vXs+AMCqqKhw82xD0zHiyTAAAACiRTMMAACAaNEMAwAAIFo0wwAAAIgWA3SOsrIykz3zzDOJMqBY/fWvf3XzyZMnm8zbla5bt245rX/11Veb7E9/+pPJ2G0OAKrHaaedZrIpU6akUElh8WQYAAAA0aIZBgAAQLRohgEAABAtmmEAAABEiwE6ACIism7dOjf/5S9/WeBKAAD5snXr1sTH1q0bZ1vIk2EAAABEi2YYAAAA0aIZBgAAQLRohgEAABAtmmEAAABES0MIhVtMtXCLodYKIWih1+TeRT4U+t7lvkU+8Jlb3Jo2berm3tb227dvN1mjRo3yXlOhJL13eTIMAACAaNEMAwAAIFo0wwAAAIgWzTAAAACixQAdig7DHChWDNChGPGZi2LFAB0AAABQCZphAAAARItmGAAAANGiGQYAAEC0CjpABwAAANQkPBkGAABAtGiGAQAAEC2aYQAAAESLZhgAAADRohkGAABAtGiGAQAAEC2aYQAAAESLZhgAAADRohkGAABAtGiGAQAAEC2aYQAAAESLZhgAAADRohkGAABAtGiGAQAAEC2a4RpCVU9X1RdVdauqblbV11X15LTrArJR1ZNUdYGqblfVjap6v6q2SLsuYE9U9URVDc5fm9KuDdgTVe2lqvNUda2q7lDVj1V1uqoemnZtxa5u2gVARFUvEpHJu/8aJ5n/SOkmIg3TrAvIRlV/ICLPisgzInKWiBwgIjeKyHOq2iOEsCPN+oAELhWRRV/7+/K0CgESaiYifxeRO0XkUxFpIyKjROQVVT08hPBRmsUVMw0hpF1D1FS1nYi8JyKjQwgT060GSEZV54pIOxHpHEIo350dLSKviciIEMKdKZYHZKWqJ4rIfBHpGUKYm3I5QE5UtZOILBWRK0IIv0+7nmLF1yTSN0REKkTkrrQLAargWBH52/80wiIiIYRFIrJBRAakVhUAxGXD7l/LUq2iyNEMp69UMv9Vd46qfqiq5aq6XFVHpF0YsAe7RGSnk+8QkcMKXAuwN6ap6i5V3aCqD6pqm7QLApJQ1TqqWqKqh4jIH0RkrYg8lHJZRY3vDKev1e6/ficiY0TkQxE5W0Qmq2rdEMKkNIsDsnhfMk+H/x9VbSsiLYUnFKjZvhCR34vICyKyWUS6S+az979UtXsIYX2axQEJvCoiPXb/7+UicjL3bW74znDKVPUDETlERM4KITz6tfxpyXxItwz8S0INo6rnisgDInKTiNwmmcGOP4rI8SJSFkJokGJ5QJWo6pGS+b77zSGEsWnXA+yJqv6LiOwrIu1F5AoRaSEipSGEFWnWVcxohlOmqv8lmSds+4YQtnwtHykit4pI6xDC6rTqA7JR1XGS+SCuLyJBRB4WkUYiclgIoX2atQFVparvisiqEEKvtGsBklLVpiKyQkQeCiEMT7mcosV3htP3TpZcd/9aUahCgKoIIVwjIgeKyBGS+ROMQZL5U44FqRYG7B2VzH/UAUUjhLBJMl+V6Jh2LcWMZjh9s3b/+s2nEb1E5OMQwtoC1wMkFkLYFkL47xDCOlXtLSKdhTejoMio6lEi8j3JfBcTKBq7NzrqLJl5I+wlBujS95Rk3nn5B1U9UET+ISIDReQ0ERmcZmFANqraXUT6iMji3VGpiFwpIr8NISxMrTCgEqo6TUT+KZl7d5NkZjNGi8gnInJ7iqUBe6SqsyRz3y6RzPDn90RkpGQ2jOEdwzngO8M1gKruKyLjJdME7y+ZV63dHEJ4MNXCgCxUtYtkXulzmIh8SzIbx9weQpiaamFAJVR1tIgMEpG2ktnlc62IPC0i14UQ1qRZG7AnqnqViPxYRDqISImIrBKR50VkPMNzuaEZBgAAQLT4zjAAAACiRTMMAACAaNEMAwAAIFo0wwAAAIhWQV+tpqpM6yFnIQSt/Kj84t5FPhT63uW+RT7wmYtilfTe5ckwAAAAokUzDAAAgGjRDAMAACBaNMMAAACIFs0wAAAAokUzDAAAgGjRDAMAACBaNMMAAACIFs0wAAAAokUzDAAAgGjRDAMAACBaNMMAAACIFs0wAAAAokUzDAAAgGjRDAMAACBaddMuIBZ9+vQx2ciRI03Ws2dPk4UQTLZs2TJ3nenTp5tsypQpJlu9erV7PgAAQEx4MgwAAIBo0QwDAAAgWjTDAAAAiBbNMAAAAKKl3nBWtS2mWrjFUnLxxRe7+YQJE0xWUlJS3eWIiMj8+fNNdt5555lszZo1hSgnZyEELfSaMdy7qH6Fvne5b5EPfOZWzX333Weyn/70pyabM2eOe/7MmTNNtnDhQpOtWrUqUT07d+508127diU6v5glvXd5MgwAAIBo0QwDAAAgWjTDAAAAiBbNMAAAAKLFDnQ5OOOMM0x2yy23uMd6w3JvvPGGyUaNGmWyd955J3FNF1xwgcluuOEGk40ePdpkl156aeJ1UNwaNWpksjFjxrjHjh071mTe4O24ceNM1rVrV5P169cvSYkAUJSWLl1qsoqKCpN5PcSe8r01depUN7/oootMVl5ente1iwVPhgEAABAtmmEAAABEi2YYAAAA0aIZBgAAQLTYgS6hvn37muwvf/mLybzBJBGR2bNnm8zbrW7dunV7Ud3/UrWbrXhDdaeddprJfvzjH+e0dqGwG1Lu2rRpY7KPPvrIPbZHjx4mW7x4scm8Abpf/OIXJuvUqZO7Tq73fjFgBzp8XYsWLUzWsWNH99j69eubbNCgQSabNm2aybLtQPbyyy9XVqKI8JmbD14P0atXr8TnH3300SbzPscbNGhgsv3228+95imnnGIyb8faYsYOdAAAAEAlaIYBAAAQLZphAAAARItmGAAAANFiBzpH3br2/xZvFzdvWG7JkiXuNb2dXj799NO9qG7PvIHIu+++22SzZs3K+9ooHu3atcv7NcvKykzmDW4ceuih7vkxDNAhDocddpjJ/u3f/s1kQ4YMMVnLli3dayYddh88eHCi40RE6tSpk/hY5ObJJ59MlOWqT58+JpszZ4577Omnn26y2jZAlxRPhgEAABAtmmEAAABEi2YYAAAA0aIZBgAAQLRohgEAABAt3ibhGDZsmMm6d+9ush07dpjs/PPPd69ZHW+OyMWGDRvSLgEpOu644/J+zccee8xk3ltYjjrqKPf8WKeYURy6devm5iNHjjTZqaeearKDDjoo7zV5tmzZYrJ58+YVZG0UVrNmzUx23XXXmay8vNw9P9tbJmLEk2EAAABEi2YYAAAA0aIZBgAAQLRohgEAABAtBugcv/jFLxIdN3z4cJO9+eab+S4HyIm35epZZ51lsoqKCvf8bMMXQFV5W92LiNSvX99kW7dure5yRMQf6Jw6darJOnTo4J7/rW99K+81ed59912TjR071mTecPSCBQuqpSbkpkmTJm5eWlpqspKSEpNdffXVJvPu5z//+c/uOs8//3wlFcaDJ8MAAACIFs0wAAAAokUzDAAAgGjRDAMAACBaDNDl4OOPP067BKBSLVq0MNnRRx9tsn/+85/u+UuWLEm0TllZmcl27dplso4dOya6Hmofb3csEZEzzzzTZDNnzjTZ9ddfn3itI444wmRXXXWVybxh0nr16plMVd11QgiJa0rC++cWEfn3f/93k23fvj2vayM/GjdubLLx48ebzLv3RHLbrfDVV1812c0337zX14sFT4YBAAAQLZphAAAARItmGAAAANGiGQYAAEC0oh6g8wYsREQOOeQQk23ZssVk77//ft5rAtKybNmynM5fvny5yVatWmWybt265bQOisO+++5rsp/+9KfusW3atDFZly5dTOYNJnXq1Mm95hlnnFFZiVWSbYDO4+0Cd//995vs0UcfNRm7xRW/E044wWQjRowoyNrez0i23UXxv3gyDAAAgGjRDAMAACBaNMMAAACIFs0wAAAAokUzDAAAgGhF/TaJunX9f/w6deqY7MsvvzQZ2zGjGJx88smJjpswYUJO63g/T97PUsuWLd3zvbcPbN68OaeakJ5mzZqZrFGjRu6xSbc0HjlypMmqY5vkRYsWmezhhx92j33qqadMtnXrVpN98skne10PiktpaWlO569fv95kU6ZMMdk++9jnmddcc43JvK2gRUSGDh1qss8//zxJibUOT4YBAAAQLZphAAAARItmGAAAANGiGQYAAEC0oh6gS9sBBxxgsr59+5rs8ssvT3zNFStWmKxdu3YmW7t2rckeeeQRk02dOtVdp6ysLHFNSNfxxx9vsnXr1pnspZdeymkdb8h0zpw5Jhs+fLh7/n777WcyBuiKl/dZ9Omnn7rHesN2hTJu3DiT3XbbbSbbuHFjIcpBLXDDDTeY7O9//7vJtm3b5p7/wgsvmGznzp0m84ZHZ8yYYbLnnnvOXefuu+822QUXXGCyTZs2uefXJjwZBgAAQLRohgEAABAtmmEAAABEi2aO05neAAAI90lEQVQYAAAA0WKALiFvwOOoo44y2euvv+6e37FjR5PNnTvXZG3atDHZ9u3bTfbWW2+563hDK142ePBgk5166qkm69Wrl7vOWWed5eZIl7fD1+mnn24ybxgj2zBHLmIYvEBy2QZ5OnXqtNfXfPHFF9185syZJnvwwQdN5u24VVFRsdf1AOXl5SabPXt23tfxdll8++23TTZs2DD3/FmzZpls/vz5Jps8efJeVFdceDIMAACAaNEMAwAAIFo0wwAAAIgWzTAAAACiFfUAXbYdhb744guTebtjeVn79u3da86bN89kBx98sMm8AZMRI0aY7IMPPnDXSerxxx83mfdl+s6dO+e0DgqrYcOGJmvbtq3JVq1aVYhy3J+lbLyfp0LVicIYPXq0m3s7b3rDxJ4TTzwxl5KAWs/7/V5E5KGHHjKZ9zP68MMPmyzbbpLFiifDAAAAiBbNMAAAAKJFMwwAAIBo0QwDAAAgWlEP0Hk7s4mIrFmzxmTecM9PfvITkx166KHuNb1hOW8HugEDBpisOnYG89a+++67TXbaaaflfW2kr6SkxGQ9evRwj/3qq69M5g2fNmjQwGTeDknZTJkyxWQnn3yyycrKyhJfEzXL1q1b3dwb5DnvvPNM1rp1a5OtXbvWveaMGTNMdt1115ks2yA1UNtNmjTJZIMGDTLZhRdeaLKbbrqpWmpKC0+GAQAAEC2aYQAAAESLZhgAAADRohkGAABAtGiGAQAAEC2tyrR3zoupFm6xHIwfP95kV111VU7X9N7UcNlll5nsyy+/zGmdXDz44IMm6927t3tst27dTLZy5cq81+QJIWhBFvqaYrl3v/3tb5ts/fr1OV2zvLzcZN5bAbw3VHjbQ1eF93aV2bNn53TNNBX63i2W+9bjTbXfddddJmvSpIl7vvd728KFC03Wr18/k33++edJSowGn7m1U/369U328ssvm2zJkiUmGzx4cLXUlG9J712eDAMAACBaNMMAAACIFs0wAAAAokUzDAAAgGgxQOdo2rSpyd58802TtWnTJvE1f/WrX5ls4sSJVSusmnnbmmYbTjnyyCNN9v777+e9Jg/DHNnVqVPHZOPGjTPZ6NGjC1FOlbz++usmO/bYY022a9euQpRTLRigy433mesNJ4uInHLKKYmu+e6775rs7LPPNtnSpUsTXa824jPX3wpcxB/qHDhwoMl27NiR95qqw9ixY0120UUXmezwww832aZNm6qlplwwQAcAAABUgmYYAAAA0aIZBgAAQLRohgEAABAtBugSOuOMM0z20EMPmaxRo0bu+du2bTPZk08+abKbbrrJZG+//XaSEqukT58+Jnv88cdN9sEHH7jnd+nSJe81JcUwR9V4Q3XNmzc3WbZ717tXvKEjL/OGLJ555hl3HW93sBNOOME9tlgxQJd/3pCliL9TobdDo2fRokUmu+SSS9xjvcHP2obPXJF27dq5+T/+8Q+T3X///Sb7j//4D5OtW7cu57ryzRug+/Wvf22y9u3bm2zFihXVUVJOGKADAAAAKkEzDAAAgGjRDAMAACBaNMMAAACIFgN0OejVq5fJfvOb37jHHnHEEYmuuX37dpMNHTrUZCtXrnTP977AXlpaarJJkyaZzNt57y9/+Yu7zuDBg928EBjmKB49evQwWbaBIwbo8i/m+/bMM8802cyZM/f6et7nsIjI1KlT9/qaxYLPXJFWrVq5ubfzqjeMvGzZMpMNHz7cveZLL71ksvLy8spKrLIBAwaY7JZbbjFZSUmJyQ477DCTffHFF/kpLI8YoAMAAAAqQTMMAACAaNEMAwAAIFo0wwAAAIgWA3R5lm2HoyFDhpjM25Fm//33z3tNHu/L+N7udzfccEMhyqkShjmKx4EHHmiypUuXusfu2rXLZN/73vdMVhOHNJJigC7/Lr74Yje/44478rrOvffe6+beZ3ttw2dudgMHDjTZ9OnTc7qmtzOd16s99thjJuvfv3/idZo1a2Yyb1juxhtvNNm1116beJ00MUAHAAAAVIJmGAAAANGiGQYAAEC0aIYBAAAQLZphAAAARIu3SaTIm+T0JqO9adWuXbsmXmfVqlUmu+uuu0w2fvz4xNdME5PNxc3bdllE5LjjjjOZtwXqmjVr8l5TofA2ieS87e5Hjx5tsh/+8Ifu+fn+ve2SSy5x8ylTpuR1nZqIz9zs6tSpY7LevXubbNSoUSbLdbt5VfuvJdf7/u677zbZ1VdfbbJPP/00p3UKhbdJAAAAAJWgGQYAAEC0aIYBAAAQLZphAAAARIsBOhQdhjmK28iRI9381ltvNdmZZ55pMm8L0mIR+wBdnz593PzCCy80mTeE5G0V6w0RiSQfJBo3bpzJFi9ebLLHH3880fVqIz5zc7fPPvbZ4/e//333WG9o/vjjjzfZsccea7KdO3eabMaMGe46kyZNMpl371dUVLjnFwMG6AAAAIBK0AwDAAAgWjTDAAAAiBbNMAAAAKLFAB2KDsMcxe2YY45x81deecVkzz//vMlOOumkfJdUMDEN0A0dOtRk2Xa59Hbj9GzatMlkCxYscI996623TPboo4+abMmSJSYr5oGh6sBnLooVA3QAAABAJWiGAQAAEC2aYQAAAESLZhgAAADRYoAORYdhDhSrmAbovB2zzjjjDPfYOXPmJLrm+vXrTbZ8+fKqFYYq4zMXxYoBOgAAAKASNMMAAACIFs0wAAAAokUzDAAAgGjRDAMAACBavE0CRYfJZhSrmN4mgdqDz1wUK94mAQAAAFSCZhgAAADRohkGAABAtGiGAQAAEC2aYQAAAESLZhgAAADRohkGAABAtGiGAQAAEC2aYQAAAESroDvQAQAAADUJT4YBAAAQLZphAAAARItmGAAAANGiGQYAAEC0aIYBAAAQLZphAAAARItmGAAAANGiGQYAAEC0aIYBAAAQLZphAAAARItmGAAAANGiGQYAAEC0aIYBAAAQLZphAAAARItmGAAAANGiGQYAAEC0aIYBAAAQLZphAAAARItmGAAAANGiGQYAAEC0aIYBAAAQLZphAAAARItmGAAAANH6v6Gn6gIIP1QkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plots(x_imgs[:8], titles=preds[:8])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining Logistic Regression Ourselves"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Above, we used pytorch's nn.Linear to create a linear layer. This is defined by a matrix multiplication and then an addition (these are also called affine transformations). Let's try defining this ourselves.\n",
    "\n",
    "Just as Numpy has np.matmul for matrix multiplication (in Python 3, this is equivalent to the @ operator), PyTorch has torch.matmul.\n",
    "\n",
    "Our PyTorch class needs two things: constructor (says what the parameters are) and a forward method (how to calculate a prediction using those parameters) The method forward describes how the neural net converts inputs to outputs.\n",
    "\n",
    "In PyTorch, the optimizer knows to try to optimize any attribute of type Parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_weights(*dims): return nn.Parameter(torch.randn(dims)/dims[0])\n",
    "def softmax(x): return torch.exp(x)/(torch.exp(x).sum(dim=1)[:,None])\n",
    "\n",
    "class LogReg(nn.Module):\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "        self.l1_w = get_weights(28*28, 10)  # Layer 1 weights\n",
    "        self.l1_b = get_weights(10)         # Layer 1 bias\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = x.view(x.size(0), -1)\n",
    "        x = (x @ self.l1_w) + self.l1_b  # Linear Layer\n",
    "        x = torch.log(softmax(x)) # Non-linear (LogSoftmax) Layer\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create our neural net and the optimizer. (We will use the same loss and metrics from above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "net2 = LogReg()#.cuda()\n",
    "opt=optim.Adam(net2.parameters())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9451d15473aa4aa1a1462961629b7f38",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Epoch', max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch      trn_loss   val_loss   accuracy                     \n",
      "    0      0.300772   0.28667    0.92      \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[array([0.28667]), 0.92]"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit(net2, md, n_epochs=1, crit=loss, opt=opt, metrics=metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "dl = iter(md.trn_dl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmb,ymb = next(dl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       "-0.4245 -0.4245 -0.4245  ...  -0.4245 -0.4245 -0.4245\n",
       "-0.4245 -0.4245 -0.4245  ...  -0.4245 -0.4245 -0.4245\n",
       "-0.4245 -0.4245 -0.4245  ...  -0.4245 -0.4245 -0.4245\n",
       "          ...             ⋱             ...          \n",
       "-0.4245 -0.4245 -0.4245  ...  -0.4245 -0.4245 -0.4245\n",
       "-0.4245 -0.4245 -0.4245  ...  -0.4245 -0.4245 -0.4245\n",
       "-0.4245 -0.4245 -0.4245  ...  -0.4245 -0.4245 -0.4245\n",
       "[torch.FloatTensor of size 64x784]"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vxmb = Variable(xmb)#.cuda())\n",
    "vxmb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       "\n",
       "Columns 0 to 5 \n",
       " 3.6239e-07  1.1758e-11  1.6963e-06  1.0115e-03  1.2447e-07  1.6883e-08\n",
       " 9.9958e-01  2.4000e-13  5.2309e-06  1.3328e-06  6.2637e-09  3.8851e-04\n",
       " 2.6365e-07  9.3130e-01  3.2980e-03  1.0834e-02  1.5949e-04  2.1343e-04\n",
       "\n",
       "Columns 6 to 9 \n",
       " 4.0005e-10  9.9896e-01  1.8901e-06  2.7488e-05\n",
       " 7.3608e-06  5.5254e-09  1.2999e-05  2.2519e-07\n",
       " 8.6490e-04  3.7555e-04  5.1849e-02  1.1058e-03\n",
       "[torch.FloatTensor of size 3x10]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds = net2(vxmb).exp(); preds[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " 7\n",
       " 0\n",
       " 1\n",
       " 0\n",
       " 5\n",
       " 1\n",
       " 2\n",
       " 5\n",
       " 7\n",
       " 6\n",
       " 5\n",
       " 4\n",
       " 7\n",
       " 0\n",
       " 0\n",
       " 1\n",
       " 7\n",
       " 3\n",
       " 9\n",
       " 4\n",
       " 3\n",
       " 4\n",
       " 1\n",
       " 8\n",
       " 3\n",
       " 9\n",
       " 4\n",
       " 9\n",
       " 0\n",
       " 9\n",
       " 1\n",
       " 4\n",
       " 9\n",
       " 8\n",
       " 8\n",
       " 6\n",
       " 8\n",
       " 1\n",
       " 9\n",
       " 0\n",
       " 5\n",
       " 1\n",
       " 0\n",
       " 2\n",
       " 6\n",
       " 9\n",
       " 7\n",
       " 5\n",
       " 1\n",
       " 5\n",
       " 5\n",
       " 7\n",
       " 3\n",
       " 4\n",
       " 1\n",
       " 1\n",
       " 3\n",
       " 1\n",
       " 5\n",
       " 1\n",
       " 3\n",
       " 6\n",
       " 1\n",
       " 3\n",
       "[torch.LongTensor of size 64]"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds = preds.data.max(1)[1]; preds"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Let us take a look at some predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAF0CAYAAADGqzQSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmUFeXV7/H9CLSMiqiMCggYvKIC4mwncQChhQUSMRHFaDMoiBpxeGVSUeSiiRFQFJOgOIHK7IBTEByQV8WAckEQMEFQJgWRUemmn/tH910vl70PXd3n9Kk+/Xw/a7mIv1TVs02qD9vi7Hqc914AAACAEB0WdwEAAABAXGiGAQAAECyaYQAAAASLZhgAAADBohkGAABAsGiGAQAAECya4Zg55zo65+Y55zY5535xzn3rnJvqnDs57tqAQ3HOne+ce8c5t8U5t8M5t9g51zvuuoAonHOXOuc+cM7tKrp/P3POXRR3XcChOOcudM4tcM7tdc5tc84975yrF3ddmY5mOH51RORfInKTiFwiIkNEpJWIfOycaxJnYUAizrnTRGSuiFQRkX4icrmILBKRp5xzA+KsDSiOc+4GEXlFCj97u4vIFSIyTUSqx1kXcCjOuV+LyDsisl0KP3P/JCK/EZF3nXOHx1lbpnNsulH+OOdaishKEbnDe//XuOsBDuac+98icoeI1PHe7zog/1hEvPf+3NiKAw7BOddURFaIyBDv/dh4qwGic87NFZGmInKS9z6/KDtTRD4VkYHe+ydiLC+j8WS4fNpa9GterFUAiWVJ4f2596B8u/C5gvKtt4gUiMiTcRcClNA5IvLP/9cIi4h47xdJYc/QPbaqKgB+0yonnHOVnHNZzrkTReRvIrJJRF6KuSwgkWeKfn3UOdfQOVfbOddPRC4WkTHxlQUUK1sK/+TtSufc1865fOfcGufcwLgLA4qxX0T2GfkvInJKmmupUPiaRDnhnPtMRNoV/e0aEenqvV8RY0nAIRX98dwsEWlUFOWJyADv/VPxVQUcmnNupYg0lMIGYqiIfC2F3xnuLyK3eu/HxVgekJBz7lMp/Bra2QdkTUTkPyKS573ne8OlRDNcTjjn/peIHCEizaTwu5j1RCTbe782zroAS9GfYLwrhd+9fEwKvy7RTUQGiMh13vvJMZYHJOScWyUiJ4rI5d77mQfkb4pIWxFp4PmNEeWQc+5qEXlBREaJyKNSOID/dxE5Twqb4WoxlpfRaIbLIedcbRFZKyIvee/7x1wOoDjnponI6VI4yJF3QD5ZRDqKSF3vfUFc9QGJOOf+Wwq/e3mE937nAfkgEXlERBp57zfEVR9wKM65kVL4wKyqiHgReVlEaojIKd77ZnHWlsn4znA55L3fLoVflWgRdy1AAqeKyBcHNsJFPhWRo0WkbvpLAiJZniB3Rb/yL3Eot7z3d4vIMSJymhT+KUZPKfyTjgWxFpbhaIbLoaIXaJ8khd9lA8qjTSLSxjmXdVB+toj8LCLb0l8SEMmsol87HpR3FJFvvfeb0lwPUCLe+93e+//jvd/snOskhf0Cb0dJQuW4Cwidc26WiCwWkaUiskNEfiUig0QkX0R4xzDKq/FSuEnBa865J6TwO8NdRaSniIzx3lsTz0B58IaIzBeRvznnjhGRf4tIDync9Cg3zsKAQ3HOtRWRHCnsGUQK34xyp4j82Xu/MLbCKgC+Mxwz59xdIvJ7EWkuhe9uXS8i74nIaIbnUJ4553JE5C4p3DGxqhT+ScbfReRv3vv9cdYGHIpz7ggRGS2FTfBRUviqtQe991NiLQw4BOdcKyl89eopInK4FA0we+8nxVpYBUAzDAAAgGDxnWEAAAAEi2YYAAAAwaIZBgAAQLBohgEAABCstL5azTnHtB6S5r13xR+VWty7SIV037vct0gFPnORqaLeuzwZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLDSugMdAADlSc2aNVXWp08flXXr1s08v2vXrirbtWtX8oUBSBueDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYvE0CABCsa6+9VmVjxoyJfH6rVq1U9sknnyRVE4D04skwAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFgN0SWjdurXKBg0aZB7bvHlzlVWvXl1lQ4cOVdmRRx6psjfffNNcZ+fOnWYOAKG77rrrVDZ27FiV5eXlqezhhx82r7l48eKk6wIQL54MAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYDnvffoWcy59i6VYzZo1VbZu3TqV1a5dOx3lyHfffWfm1gDf9OnTy7qctPLeu3Svmcn3rsW6T7t3724e27ZtW5VlZ2erzPoZ2bZtm8rq169vrrNp0yaVPfPMMyr7xz/+obL9+/eb1yxv0n3vVrT7tiS6du2qslmzZqlsz549KrvnnntUVpJd6SoaPnORqaLeuzwZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwWKALqJatWqp7I033lDZ1q1bzfOXLFmiMmswqUmTJio7/vjjVVatWjVznc2bN6vs3HPPjXRcpmCYo2SOO+44lc2ePVtl1v2YyI4dO1Rm3eNVqlRRmfWzJCJSt25dldWrV09lV111lco++OADlW3cuNFcJ04M0KVeVlaWmU+aNEllPXv2VNm8efNU1r59++QLq0D4zEWmYoAOAAAAKAbNMAAAAIJFMwwAAIBg0QwDAAAgWAzQZYBjjjlGZXfeead5rJXn5uaq7Nlnn02+sJgwzFEyixcvVlnr1q1VNnfuXPP822+/XWU//PCDyqwd5Eri2GOPVdmbb76pspYtW6ps8ODBKnv88ceTqqcsMECXesOGDTPzkSNHquyFF15QWe/evVWWn5+ffGEVCJ+5yWvQoIHKbrzxRvNYK8/Ly1OZtQvuqFGjVGb9HiAisn79ejOvSBigAwAAAIpBMwwAAIBg0QwDAAAgWDTDAAAACBbNMAAAAILF2yQyVNeuXc3c2mb30UcfVdmtt96a8prShcnmxKyJ5e+++05lU6dOVdnVV19tXnP//v3JF1ZKkydPVtmVV16psnbt2qns888/L5OaksHbJJJzxhlnqGzBggXmsWvXrlVZq1atVBbn/Z0p+MwtmWbNmqlswoQJKuvQoUM6ypFffvnFzM8//3yVJXrzRKbibRIAAABAMWiGAQAAECyaYQAAAASLZhgAAADBqhx3ASjeUUcdpbKhQ4dGPr9hw4apLAflWJs2bVTmnJ4f2LBhg8riHiQ655xzVNazZ0+VzZ8/X2XWP3d5HKBDdIcdpp/VWNtuZ2Vlmee/9tprKov7HkfF06hRI5UtW7ZMZZUr63ZrzJgx5jUfe+yxSOucdNJJKvvLX/6istq1a5vrWIPU1ufwDz/8YJ5fkfBkGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABIsd6MqZ1q1bq2zatGkqa9GihXn+qlWrVGbtcrN+/fpSVFc+sBtSyRQUFKhsy5YtKjvrrLPM89etW5fSemrVqmXmCxcuVNnq1atVZu2UZ+34tHz58lJUV7bYgS66qLspJnLLLbeobPz48UnVFCo+cxMbN26cyvr376+yfv36qey5555LeT0DBw5U2dixY81jK1WqpLKVK1eqzBqq27FjRymqSz92oAMAAACKQTMMAACAYNEMAwAAIFg0wwAAAAgWA3Qxuvbaa1V2//33q+z4449X2d69e81rdunSRWXWjl2ZjGGOkhkxYoTK7r77bpV99dVX5vkdO3ZUWTIDmO+8846Z//a3v1VZu3btVGbt7pQpGKCLLjc3V2VPPfWUyubOnWuen5OTozJ2oCsdPnNFjjjiCDO3hnwnTZqkMmv3xHRJ9Nl+4oknRjrf2inv9ttvT6qmdGGADgAAACgGzTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAgWb5NIsZo1a5r5HXfcobLhw4er7LDD9L+fbNu2TWXZ2dnmOtZWihUNk80lU7VqVZU9++yzKuvRo4d5/po1a1R2wQUXqGzjxo0qe+KJJ1R2/fXXm+vceeedKrOmmDMZb5OwVa5cWWUrVqxQWZMmTVR2wgknmNcsydbNODQ+cxNvV//xxx+rrEOHDip79913U15TVN27dzfzmTNnqszqCbdv364y600UW7duLUV1ZYu3SQAAAADFoBkGAABAsGiGAQAAECyaYQAAAARLTy0gKc8884yZ/+53v4t0/vTp01U2duxYlYUwKIfU+Pnnn1XWt29fldWtW9c839om+f3331fZtGnTVNarVy+VzZgxw1ynog3LITpreLN58+YqGzBggMriHpTr1KmTyrp27aqyt956S2XW1uTWzyvi17Zt28jHLlmypAwrKbk33njDzK3haOvnzrond+/enXxh5QhPhgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLAYoEsx68vnJTFhwgSVLVy4MKlrAgfbuXOnyrp162YeO2LECJXdeuutKhs8eHCktR977LFIxyEcjRs3jnRcVlZWGVeS2HXXXWfm1i6L1q6P/fv3V5m1s9fs2bPNdXr37l1MhShLCxYsMPOCggKV/fOf/1RZly5dVGbt2lkWWrZsaebWfdqxY0eVVa9eXWUVbdCTJ8MAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYDNClmLWjkIhI69atS32+NVT34IMPmudv2LAh0jrAwXbs2GHm99xzj8o6dOigspNPPjnSOu3btzfzRAMqqPhatGgR6bh07bxZu3ZtlT3yyCPmsdYQUn5+vsqsoars7GyVWbs2ijBAF7fly5eb+euvv64yaxh5xYoVKrN2JRSxd+mcN2+eyho1aqQya1jO2sVWRKRBgwYqs+7dV155xTy/IuHJMAAAAIJFMwwAAIBg0QwDAAAgWDTDAAAACJbz3qdvMefSt1hMqlWrZuYvvPCCytq1a6eyqDsxbdq0ycxzc3NV9vbbb0e6Zqbw3rt0rxnCvZtITk6OymbNmqWyKlWqRLrevn37zPzGG29U2aRJkyJdM1Ok+97NlPt2zpw5Kmvbtq3KGjZsmI5yzB0WEw3QWZ/t48aNU9m6detUZg1QnXrqqeY6ce6+x2duYtbv+aNHj1bZLbfcktQ627ZtU1mdOnWSuqbliiuuUJk10Jcpot67PBkGAABAsGiGAQAAECyaYQAAAASLZhgAAADBohkGAABAsNiOOcX27t1r5ldffbXKKlfW//Mn2hL3YPXr1zdza8r/tttuU9mTTz4ZaR3gwgsvVJn1Fpru3burzJqAtrYvFbG3Hf/hhx9U9tprr5nnI3OdffbZKkv01pHyZsOGDSo77rjjVPb3v/9dZaeffrrKKtrbfyo66/d8620kU6dOVZnVFyRSr169SMfl5eWpzPr5EhE54YQTVLZnz57INVUkPBkGAABAsGiGAQAAECyaYQAAAASLZhgAAADBYjvmcua0005T2ZgxY1RmDTUlYm0D2rRp0xLVVZ6wNWjZsO49EZFFixapzBp2s4ZGLNZ2nyIiTz31lMqc0/9Xt2rVSmXWPV4esR2zzRou69Kli8rKYjtm6x6z7uW//vWvSa1j/V77xBNPqGzo0KHm+Tt37kxq/WTwmZvZnn/+eTO3Bvg6deqksnfeeSflNaUL2zEDAAAAxaAZBgAAQLBohgEAABAsmmEAAAAEix3oIqpevbrKymKnlqVLl6qsR48eKnv66afN87t166ayxo0bq6xBgwYq27hxY5QSUUHVqlXLzK2dEqdPn17qdaZNm2bmTZo0UdlDDz2ksnbt2qksUwboEF3t2rVVZg0CvfDCC+b51n175ZVXqqxOnToqy8nJiVKiiIjs3r1bZQsWLFDZn//8Z5XNnz8/8jpAOjRv3jzuEmLBk2EAAAAEi2YYAAAAwaIZBgAAQLBohgEAABAsBugM1hfIrYGIOXPmqGzZsmXmNa3htD59+qisSpUqKmvUqJHKWrRoYa5j+frrryPVg7C1adPGzDdt2qQy6+chWePHj1dZv379VDZw4ECVzZo1K+X1IH2WLFmisr59+6rM2jHLypK1Y8cOlSUa/HzggQdU9s0336S8JqC0du3aFXcJ5R5PhgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLAYoDNcccUVKqtfv77KevfunfK1nXMq895HPt/6onz//v2TqglhsHYqFBH59NNP07L+vn37VPbjjz+q7Ne//rXKrF3Etm3blprCUOamTJmiMmvnzdWrV6usUqVK5jUT5QebPHmyytauXasyaxAZyAQffPCBmd9www0qq1u3blmXUy7xZBgAAADBohkGAABAsGiGAQAAECyaYQAAAASLZhgAAADB4m0ShqOPPjruEv4/M2bMUNnIkSPNY7ds2aIyaztd4GCJ3lqSnZ2tsiuvvFJl8+bNU1nNmjVVlpWVZa5z0kknqezMM89U2eOPP64y3hyR2X766SeVXXzxxTFUAlQ8hx1mP/e03l5l9RAh4MkwAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFgN0hqFDh6ps7ty5KuvVq5fKGjZsaF7TGhCxPPbYYyr78MMPVZafnx/pekBUK1asMHNrq2Nr+9ytW7eqrCQDdNYwx0cffaSyESNGmOcDALSCggIzTzQ0HSKeDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGAxQGfIy8tT2dtvvx0pAzLVW2+9Zebjx49XmbUrXZs2bZJaf9iwYSp7+umnVcZucwBQNi655BKVTZgwIYZK0osnwwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFgM0AEQEZHNmzeb+Z/+9Kc0VwIASJVdu3ZFPrZy5TDbQp4MAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFjOe5++xZxL32KosLz3Lt1rcu8iFdJ973LfIhX4zM1stWvXNnNra/u9e/eqrEaNGimvKV2i3rs8GQYAAECwaIYBAAAQLJphAAAABItmGAAAAMFigA4Zh2EOZCoG6JCJ+MxFpmKADgAAACgGzTAAAACCRTMMAACAYNEMAwAAIFhpHaADAAAAyhOeDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFs1wOeGcu9Q594Fzbpdzbodz7jPn3EVx1wVE5Zx7yznnnXMPxF0LkIhz7oKi+/Tgv7bHXRtwKNy7Zady3AVAxDl3g4iML/prpBT+S0obEakeZ11AVM65niLSOu46gBK4RUQWHfD3+XEVApQQ926K0QzHzDnXVETGisid3vuxB/xXb8dSEFBCzrnaIjJGRAaJyJSYywGiWuG9/zjuIoBS4N5NMb4mEb/eIlIgIk/GXQhQSn8WkeXe+xfjLgQAgJKiGY5ftoisFJErnXNfO+fynXNrnHMD4y4MKI5zLltE/igiN8ZdC1BCk51z+51zW51zU5xzjeMuCIiIezfF+JpE/BoW/fUXERkqIl+LyBUiMt45V9l7Py7O4oBEnHNVRORvIvKw9/6ruOsBIvpJRP4qIu+LyA4RaSuFn73/7Zxr673fEmdxwCFw75YR572Pu4agOedWiciJInK5937mAfmbUnijN/D8n4RyyDk3XAq/5tPKe7+3KPMiMsp7PzzW4oAScM6dLiKfisiD3LvIJNy7qcHXJOK3tejXfx6UvyMi9USkQXrLAYpX9Mdyw0TkbhE53DlXu2iQTg74+0rxVQhE571fLCKrROTMuGsBSoJ7NzVohuO3PEHuin4tSFchQAk0E5GqIvKCiPx4wF8iIncU/edT4ykNKBUnIvwpHDIR926SaIbjN6vo144H5R1F5Fvv/aY01wNE8bmIXGj8JVLYIF8oImviKQ0oGefcGSLyKxH5JO5agJLg3k0NBuji94aIzBeRvznnjhGRf4tIDxG5RERy4ywMSMR7v11E3js4d86JiHzjvVf/HVAeOOcmi8h/RGSxiGyXwtmMISLynYg8FmNpwCFx75YdmuGYee+9c+4yERktIveJyFFS+Kq1q733bGAAAKm1TER6isjNUrjL5yYRmSki93rvf4izMKAY3LtlhLdJAAAAIFh8ZxgAAADBohkGAABAsGiGAQAAECyaYQAAAAQrrW+TKNqqFUiK994Vf1Rqce8iFdJ973LfIhX4zEWminrv8mQYAAAAwaIZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwaIZBgAAQLBohgEAABCsynEXEIqcnByVDRo0SGUdOnRQmfdeZatXrzbXmTp1qsomTJigsg0bNpjnAwAAhIQnwwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFjOGs4qs8WcS99iMRkwYICZjxkzRmVZWVllXY6IiMyfP19lvXr1UtnGjRvTUU7SvPcu3WuGcO+i7KX73uW+RSrwmVsyzz77rMquueYalc2ZM8c8f8aMGSpbuHChytavXx+pnn379pn5/v37I52fyaLeuzwZBgAAQLBohgEAABAsmmEAAAAEi2YYAAAAwWIHuiR07txZZQ8//LB5rDUst2TJEpUNHjxYZcuXL49cU58+fVR23333qWzIkCEqu+WWWyKvg8xWo0YNlQ0dOtQ8dvjw4SqzBm9HjhypstatW6usa9euUUoEgIy0cuVKlRUUFKjM6iEOlZfWpEmTzPyGG25QWX5+fkrXzhQ8GQYAAECwaIYBAAAQLJphAAAABItmGAAAAMFiB7qIunTporIXX3xRZdZgkojI7NmzVWbtVrd58+ZSVPc/nNObrVhDdZdcconKfv/73ye1drqwG1LyGjdurLJvvvnGPLZdu3YqW7x4scqsAbqbb75ZZS1btjTXSfbezwTsQIcD1atXT2UtWrQwj61atarKevbsqbLJkyerLNEOZB999FFxJYoIn7mpYPUQHTt2jHz+mWeeqTLrc7xatWoqO/LII81rXnzxxSqzdqzNZOxABwAAABSDZhgAAADBohkGAABAsGiGAQAAECx2oDNUrqz/Z7F2cbOG5ZYuXWpe09rp5fvvvy9FdYdmDUROnDhRZbNmzUr52sgcTZs2Tfk18/LyVGYNbpx88snm+SEM0CEMp5xyisr+8Ic/qKx3794qa9CggXnNqMPuubm5kY4TEalUqVLkY5Gc119/PVKWrJycHJXNmTPHPPbSSy9VWUUboIuKJ8MAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFm+TMPTr109lbdu2Vdkvv/yisuuuu868Zlm8OSIZW7dujbsExOjcc89N+TVfeeUVlVlvYTnjjDPM80OdYkZmaNOmjZkPGjRIZe3bt1dZ/fr1U16TZefOnSqbN29eWtZGetWpU0dl9957r8ry8/PN8xO9ZSJEPBkGAABAsGiGAQAAECyaYQAAAASLZhgAAADBYoDOcPPNN0c6rn///ir7/PPPU10OkBRry9XLL79cZQUFBeb5iYYvgJKytroXEalatarKdu3aVdbliIg90Dlp0iSVNW/e3Dz/8MMPT3lNli+//FJlw4cPV5k1HL1gwYIyqQnJqVWrlplnZ2erLCsrS2XDhg1TmXU/P/fcc+Y67733XjEVhoMnwwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFgM0CXh22+/jbsEoFj16tVT2Zlnnqmy//znP+b5S5cujbROXl6eyvbv36+yFi1aRLoeKh5rdywRkcsuu0xlM2bMUNmIESMir3Xaaaep7K677lKZNUxapUoVlTnnzHW895FrisL65xYR+eMf/6iyvXv3pnRtpEbNmjVVNnr0aJVZ955IcrsVfvLJJyp78MEHS329UPBkGAAAAMGiGQYAAECwaIYBAAAQLJphAAAABCvoATprwEJE5MQTT1TZzp07VfbVV1+lvCYgLqtXr07q/DVr1qhs/fr1KmvTpk1S6yAzHHHEESq75pprzGMbN26sslatWqnMGkxq2bKlec3OnTsXV2KJJBqgs1i7wD3//PMqmzlzpsrYLS7znX/++SobOHBgWta2fkYS7S6K/8GTYQAAAASLZhgAAADBohkGAABAsGiGAQAAECyaYQAAAAQr6LdJVK5s/+NXqlRJZXv27FEZ2zEjE1x00UWRjhszZkxS61g/T9bPUoMGDczzrbcP7NixI6maEJ86deqorEaNGuaxUbc0HjRokMrKYpvkRYsWqezll182j33jjTdUtmvXLpV99913pa4HmSU7Ozup87ds2aKyCRMmqOyww/TzzLvvvltl1lbQIiJ9+/ZV2Y8//hilxAqHJ8MAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYQQ/Qxe3oo49WWZcuXVR2++23R77m2rVrVda0aVOVbdq0SWXTp09X2aRJk8x18vLyIteEeJ133nkq27x5s8o+/PDDpNaxhkznzJmjsv79+5vnH3nkkSpjgC5zWZ9F33//vXmsNWyXLiNHjlTZo48+qrJt27aloxxUAPfdd5/K/vWvf6ls9+7d5vnvv/++yvbt26cya3h02rRpKnv33XfNdSZOnKiyPn36qGz79u3m+RUJT4YBAAAQLJphAAAABItmGAAAAMGiGQYAAECwGKCLyBrwOOOMM1T22Wefmee3aNFCZXPnzlVZ48aNVbZ3716VffHFF+Y61tCKleXm5qqsffv2KuvYsaO5zuWXX27miJe1w9ell16qMmsYI9EwRzJCGLxAdIkGeVq2bFnqa37wwQdmPmPGDJVNmTJFZdaOWwUFBaWuB8jPz1fZ7NmzU76OtcvismXLVNavXz/z/FmzZqls/vz5Khs/fnwpqsssPBkGAABAsGiGAQAAECyaYQAAAASLZhgAAADBCnqALtGOQj/99JPKrN2xrKxZs2bmNefNm6ey4447TmXWgMnAgQNVtmrVKnOdqF599VWVWV+mP+mkk5JaB+lVvXp1lTVp0kRl69dfjbV6AAAH4UlEQVSvT0c55s9SItbPU7rqRHoMGTLEzK2dN61hYssFF1yQTElAhWf9fi8i8tJLL6nM+hl9+eWXVZZoN8lMxZNhAAAABItmGAAAAMGiGQYAAECwaIYBAAAQrKAH6Kyd2URENm7cqDJruOeqq65S2cknn2xe0xqWs3ag6969u8rKYmcwa+2JEyeq7JJLLkn52ohfVlaWytq1a2ce+/PPP6vMGj6tVq2ayqwdkhKZMGGCyi666CKV5eXlRb4mypddu3aZuTXI06tXL5U1atRIZZs2bTKvOW3aNJXde++9Kks0SA1UdOPGjVNZz549VXb99derbNSoUWVSU1x4MgwAAIBg0QwDAAAgWDTDAAAACBbNMAAAAIJFMwwAAIBguZJMeye9mHPpWywJo0ePVtldd92V1DWtNzXceuutKtuzZ09S6yRjypQpKuvUqZN5bJs2bVS2bt26lNdk8d67tCx0gEy5d4899liVbdmyJalr5ufnq8x6K4D1hgpre+iSsN6uMnv27KSuGad037uZct9arKn2J598UmW1atUyz7d+b1u4cKHKunbtqrIff/wxSonB4DO3YqpatarKPvroI5UtXbpUZbm5uWVSU6pFvXd5MgwAAIBg0QwDAAAgWDTDAAAACBbNMAAAAILFAJ2hdu3aKvv8889V1rhx48jXvO2221Q2duzYkhVWxqxtTRMNp5x++ukq++qrr1Jek4VhjsQqVaqkspEjR6psyJAh6SinRD777DOVnXPOOSrbv39/OsopEwzQJcf6zLWGk0VELr744kjX/PLLL1V2xRVXqGzlypWRrlcR8ZlrbwUuYg919ujRQ2W//PJLymsqC8OHD1fZDTfcoLJTTz1VZdu3by+TmpLBAB0AAABQDJphAAAABItmGAAAAMGiGQYAAECwGKCLqHPnzip76aWXVFajRg3z/N27d6vs9ddfV9moUaNUtmzZsigllkhOTo7KXn31VZWtWrXKPL9Vq1YprykqhjlKxhqqq1u3rsoS3bvWvWINHVmZNWTx9ttvm+tYu4Odf/755rGZigG61LOGLEXsnQqtHRotixYtUtlNN91kHmsNflY0fOaKNG3a1Mz//e9/q+z5559X2X/913+pbPPmzUnXlWrWAN3999+vsmbNmqls7dq1ZVFSUhigAwAAAIpBMwwAAIBg0QwDAAAgWDTDAAAACBYDdEno2LGjyh566CHz2NNOOy3SNffu3auyvn37qmzdunXm+dYX2LOzs1U2btw4lVk777344ovmOrm5uWaeDgxzZI527dqpLNHAEQN0qRfyfXvZZZepbMaMGaW+nvU5LCIyadKkUl8zU/CZK9KwYUMzt3ZetYaRV69erbL+/fub1/zwww9Vlp+fX1yJJda9e3eVPfzwwyrLyspS2SmnnKKyn376KTWFpRADdAAAAEAxaIYBAAAQLJphAAAABItmGAAAAMFigC7FEu1w1Lt3b5VZO9IcddRRKa/JYn0Z39r97r777ktHOSXCMEfmOOaYY1S2cuVK89j9+/er7Fe/+pXKyuOQRlQM0KXegAEDzPzxxx9P6TrPPPOMmVuf7RUNn7mJ9ejRQ2VTp05N6prWznRWr/bKK6+orFu3bpHXqVOnjsqsYbkHHnhAZffcc0/kdeLEAB0AAABQDJphAAAABItmGAAAAMGiGQYAAECwaIYBAAAQLN4mESNrktOajLamVVu3bh15nfXr16vsySefVNno0aMjXzNOTDZnNmvbZRGRc889V2XWFqgbN25MeU3pwtskorO2ux8yZIjKfvOb35jnp/r3tptuusnMJ0yYkNJ1yiM+cxOrVKmSyjp16qSywYMHqyzZ7ead0/+3JHvfT5w4UWXDhg1T2ffff5/UOunC2yQAAACAYtAMAwAAIFg0wwAAAAgWzTAAAACCxQAdMg7DHJlt0KBBZv7II4+o7LLLLlOZtQVppgh9gC4nJ8fMr7/+epVZQ0jWVrHWEJFI9EGikSNHqmzx4sUqe/XVVyNdryLiMzd5hx2mnz2eddZZ5rHW0Px5552nsnPOOUdl+/btU9m0adPMdcaNG6cy694vKCgwz88EDNABAAAAxaAZBgAAQLBohgEAABAsmmEAAAAEiwE6ZByGOTLb2WefbeYff/yxyt577z2VXXjhhakuKW1CGqDr27evyhLtcmntxmnZvn27yhYsWGAe+8UXX6hs5syZKlu6dKnKMnlgqCzwmYtMxQAdAAAAUAyaYQAAAASLZhgAAADBohkGAABAsBigQ8ZhmAOZKqQBOmvHrM6dO5vHzpkzJ9I1t2zZorI1a9aUrDCUGJ+5yFQM0AEAAADFoBkGAABAsGiGAQAAECyaYQAAAASLZhgAAADB4m0SyDhMNiNThfQ2CVQcfOYiU/E2CQAAAKAYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYNMMAAAAIFs0wAAAAgkUzDAAAgGDRDAMAACBYad2BDgAAAChPeDIMAACAYNEMAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAgWzTAAAACCRTMMAACAYNEMAwAAIFg0wwAAAAjW/wWeKbMQcyi29AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "preds = predict(net2, md.val_dl).argmax(1)\n",
    "plots(x_imgs[:8], titles=preds[:8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.92"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(preds == y_valid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aside about Broadcasting and Matrix Multiplication"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Now let's dig in to what we were doing with torch.matmul: matrix multiplication. First, let's start with a simpler building block: broadcasting.\n",
    "\n",
    "Element-wise operations\n",
    "Broadcasting and element-wise operations are supported in the same way by both numpy and pytorch.\n",
    "\n",
    "Operators (+,-,*,/,>,<,==) are usually element-wise.\n",
    "\n",
    "Examples of element-wise operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([10,  6, -4]), array([2, 8, 7]))"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([10, 6, -4])\n",
    "b = np.array([2, 8, 7])\n",
    "a,b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([12, 14,  3])"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6666666666666666"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(a < b).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Broadcasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The term broadcasting describes how arrays with different shapes are treated during arithmetic operations. The term broadcasting was first used by Numpy, although is now used in other libraries such as Tensorflow and Matlab; the rules can vary by library.\n",
    "\n",
    "From the Numpy Documentation:\n",
    "\n",
    "The term broadcasting describes how numpy treats arrays with \n",
    "different shapes during arithmetic operations. Subject to certain \n",
    "constraints, the smaller array is “broadcast” across the larger \n",
    "array so that they have compatible shapes. Broadcasting provides a \n",
    "means of vectorizing array operations so that looping occurs in C\n",
    "instead of Python. It does this without making needless copies of \n",
    "data and usually leads to efficient algorithm implementations.\n",
    "\n",
    "In addition to the efficiency of broadcasting, it allows developers to write less code, which typically leads to fewer errors.\n",
    "\n",
    "This section was adapted from Chapter 4 of the fast.ai Computational Linear Algebra course.\n",
    "\n",
    "http://nbviewer.jupyter.org/github/fastai/numerical-linear-algebra/blob/master/nbs/4.%20Compressed%20Sensing%20of%20CT%20Scans%20with%20Robust%20Regression.ipynb#4.-Compressed-Sensing-of-CT-Scans-with-Robust-Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10,  6, -4])"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True, False])"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a >0 "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "How are we able to do a > 0? 0 is being broadcast to have the same dimensions as a.\n",
    "\n",
    "Remember above when we normalized our dataset by subtracting the mean (a scalar) from the entire data set (a matrix) and dividing by the standard deviation (another scalar)? We were using broadcasting!\n",
    "\n",
    "Other examples of broadcasting with a scalar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11,  7, -3])"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 3],\n",
       "       [4, 5, 6],\n",
       "       [7, 8, 9]])"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.array([[1, 2, 3], [4,5,6], [7,8,9]]); m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2,  4,  6],\n",
       "       [ 8, 10, 12],\n",
       "       [14, 16, 18]])"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2*m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Broadcasting a vector to a matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We can also broadcast a vector to a matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 20, 30])"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = np.array([10,20,30]); c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 22, 33],\n",
       "       [14, 25, 36],\n",
       "       [17, 28, 39]])"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m + c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 22, 33],\n",
       "       [14, 25, 36],\n",
       "       [17, 28, 39]])"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c + m"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Although numpy does this automatically, you can also use the broadcast_to method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 10, 10],\n",
       "       [20, 20, 20],\n",
       "       [30, 30, 30]])"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(c[:,None], m.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 20, 30],\n",
       "       [10, 20, 30],\n",
       "       [10, 20, 30]])"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(np.expand_dims(c,0), (3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 3)"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.expand_dims(c,0).shape"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "The numpy expand_dims method lets us convert the 1-dimensional array c into a 2-dimensional array (although one of those dimensions has value 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 3)"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.expand_dims(c,0).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 22, 33],\n",
       "       [14, 25, 36],\n",
       "       [17, 28, 39]])"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m + np.expand_dims(c,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10],\n",
       "       [20],\n",
       "       [30]])"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.expand_dims(c,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 1)"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c[:, None].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 12, 13],\n",
       "       [24, 25, 26],\n",
       "       [37, 38, 39]])"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m + np.expand_dims(c,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 10, 10],\n",
       "       [20, 20, 20],\n",
       "       [30, 30, 30]])"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "np.broadcast_to(np.expand_dims(c,1), (3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Broadcasting Rules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 20, 30]])"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c[None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10],\n",
       "       [20],\n",
       "       [30]])"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c[:, None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False,  True,  True],\n",
       "       [False, False,  True],\n",
       "       [False, False, False]])"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c[None] > c[:,None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[0],\n",
       "        [1],\n",
       "        [2],\n",
       "        [3],\n",
       "        [4]]), array([[0, 1, 2, 3, 4]]))"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xg,yg = np.ogrid[0:5, 0:5]; xg,yg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4],\n",
       "       [1, 2, 3, 4, 5],\n",
       "       [2, 3, 4, 5, 6],\n",
       "       [3, 4, 5, 6, 7],\n",
       "       [4, 5, 6, 7, 8]])"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xg+yg"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "When operating on two arrays, Numpy/PyTorch compares their shapes element-wise. It starts with the trailing dimensions, and works its way forward. Two dimensions are compatible when\n",
    "\n",
    "they are equal, or\n",
    "one of them is 1\n",
    "Arrays do not need to have the same number of dimensions. For example, if you have a 256*256*3 array of RGB values, and you want to scale each color in the image by a different value, you can multiply the image by a one-dimensional array with 3 values. Lining up the sizes of the trailing axes of these arrays according to the broadcast rules, shows that they are compatible:\n",
    "\n",
    "Image  (3d array): 256 x 256 x 3\n",
    "Scale  (1d array):             3\n",
    "Result (3d array): 256 x 256 x 3"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "The numpy documentation includes several examples of what dimensions can and can not be broadcast together.\n",
    "https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html#general-broadcasting-rules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix Multiplication"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "We are going to use broadcasting to define matrix multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[1, 2, 3],\n",
       "        [4, 5, 6],\n",
       "        [7, 8, 9]]), array([10, 20, 30]))"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m, c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([140, 320, 500])"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m @ c  # np.matmul(m, c)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "We get the same answer using Torch matmul"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " 140\n",
       " 320\n",
       " 500\n",
       "[torch.LongTensor of size 3]"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T(m) @ T(c)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "The following is NOT matrix multiplication. What is it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[1, 2, 3],\n",
       "        [4, 5, 6],\n",
       "        [7, 8, 9]]), array([10, 20, 30]))"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m,c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 10,  40,  90],\n",
       "       [ 40, 100, 180],\n",
       "       [ 70, 160, 270]])"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m * c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([140, 320, 500])"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(m * c).sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 20, 30])"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 20, 30],\n",
       "       [10, 20, 30],\n",
       "       [10, 20, 30]])"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(c, (3,3))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "From a machine learning perspective, matrix multiplication is a way of creating features by saying how much we want to weight each input column. Different features are different weighted averages of the input columns.\n",
    "\n",
    "The website matrixmultiplication.xyz provides a nice visualization of matrix multiplcation\n",
    "http://matrixmultiplication.xyz/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 40],\n",
       "       [20,  0],\n",
       "       [30, -5]])"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = np.array([[10,40],[20,0],[30,-5]]); n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 3],\n",
       "       [4, 5, 6],\n",
       "       [7, 8, 9]])"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[140,  25],\n",
       "       [320, 130],\n",
       "       [500, 235]])"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m @ n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([140, 320, 500])"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(m * n[:,0]).sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 25, 130, 235])"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(m * n[:,1]).sum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing Our Own Training Loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a reminder, this is what we did above to write our own logistic regression class (as a pytorch neural net):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "33f81b7d7b0c41ed8ca630c11376f1ad",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Epoch', max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch      trn_loss   val_loss   accuracy                     \n",
      "    0      0.323784   0.291313   0.9157    \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[array([0.29131]), 0.9157]"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Our code from above\n",
    "class LogReg(nn.Module):\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "        self.l1_w = get_weights(28*28, 10)  # Layer 1 weights\n",
    "        self.l1_b = get_weights(10)         # Layer 1 bias\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = x.view(x.size(0), -1)\n",
    "        x = x @ self.l1_w + self.l1_b \n",
    "        return torch.log(softmax(x))\n",
    "\n",
    "net2 = LogReg()#.cuda()\n",
    "opt=optim.Adam(net2.parameters())\n",
    "\n",
    "fit(net2, md, n_epochs=1, crit=loss, opt=opt, metrics=metrics)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Above, we are using the fastai method fit to train our model. Now we will try writing the training loop ourselves.\n",
    "\n",
    "Review question: What does it mean to train a model?\n",
    "\n",
    "We will use the LogReg class we created, as well as the same loss function, learning rate, and optimizer as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "net2 = LogReg()#.cuda()\n",
    "loss=nn.NLLLoss()\n",
    "learning_rate = 1e-3\n",
    "optimizer=optim.Adam(net2.parameters(), lr=learning_rate)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "md is the ImageClassifierData object we created above. We want an iterable version of our training data (question: what does it mean for something to be iterable?):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "dl = iter(md.trn_dl) # Data loader"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "First, we will do a forward pass, which means computing the predicted y by passing x to the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "xt, yt = next(dl)\n",
    "y_pred = net2(Variable(xt))#.cuda())"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "We can check the loss:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variable containing:\n",
      " 2.3020\n",
      "[torch.FloatTensor of size 1]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "l = loss(y_pred, Variable(yt))#.cuda())\n",
    "print(l)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "We may also be interested in the accuracy. We don't expect our first predictions to be very good, because the weights of our network were initialized to random values. Our goal is to see the loss decrease (and the accuracy increase) as we train the network:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.15625"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(to_np(y_pred).argmax(axis=1) == to_np(yt))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Now we will use the optimizer to calculate which direction to step in. That is, how should we update our weights to try to decrease the loss?\n",
    "\n",
    "Pytorch has an automatic differentiation package (autograd) that takes derivatives for us, so we don't have to calculate the derivative ourselves! We just call .backward() on our loss to calculate the direction of steepest descent (the direction to lower the loss the most)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Before the backward pass, use the optimizer object to zero all of the\n",
    "# gradients for the variables it will update (which are the learnable weights\n",
    "# of the model)\n",
    "optimizer.zero_grad()\n",
    "\n",
    "# Backward pass: compute gradient of the loss with respect to model parameters\n",
    "l.backward()\n",
    "\n",
    "# Calling the step function on an Optimizer makes an update to its parameters\n",
    "optimizer.step()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Now, let's make another set of predictions and check if our loss is lower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "xt, yt = next(dl)\n",
    "y_pred = net2(Variable(xt))#.cuda())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variable containing:\n",
      " 2.2148\n",
      "[torch.FloatTensor of size 1]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "l = loss(y_pred, Variable(yt))#.cuda())\n",
    "print(l)\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Note that we are using stochastic gradient descent, so the loss is not guaranteed to be strictly better each time. The stochasticity comes from the fact that we are using mini-batches; we are just using 64 images to calculate our prediction and update the weights, not the whole dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.15625"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(to_np(y_pred).argmax(axis=1) == to_np(yt))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "\n",
    "If we run several iterations in a loop, we should see the loss decrease and the accuracy increase with time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loss:  2.2799649238586426 \t accuracy:  0.125\n",
      "loss:  1.3721935749053955 \t accuracy:  0.671875\n",
      "loss:  0.9640874266624451 \t accuracy:  0.78125\n",
      "loss:  0.7121612429618835 \t accuracy:  0.8125\n",
      "loss:  0.7122271060943604 \t accuracy:  0.78125\n",
      "loss:  0.6321592330932617 \t accuracy:  0.875\n",
      "loss:  0.6986966133117676 \t accuracy:  0.8125\n",
      "loss:  0.5703806281089783 \t accuracy:  0.828125\n",
      "loss:  0.4912716746330261 \t accuracy:  0.859375\n",
      "loss:  0.590144157409668 \t accuracy:  0.8125\n"
     ]
    }
   ],
   "source": [
    "for t in range(100):\n",
    "    xt, yt = next(dl)\n",
    "    y_pred = net2(Variable(xt))#.cuda())\n",
    "    l = loss(y_pred, Variable(yt))#.cuda())\n",
    "    \n",
    "    if t % 10 == 0:\n",
    "        accuracy = np.mean(to_np(y_pred).argmax(axis=1) == to_np(yt))\n",
    "        print(\"loss: \", l.data[0], \"\\t accuracy: \", accuracy)\n",
    "\n",
    "    optimizer.zero_grad()\n",
    "    l.backward()\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Put it all together in a training loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "def score(x, y):\n",
    "    y_pred = to_np(net2(V(x)))\n",
    "    return np.sum(y_pred.argmax(axis=1) == to_np(y))/len(y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "object of type 'generator' has no len()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-125-4f30b5fca26a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      7\u001b[0m     \u001b[0mlosses\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m     \u001b[0mdl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0miter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrn_dl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m     \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m         \u001b[0;31m# Forward pass: compute predicted y and loss by passing x to the model.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m         \u001b[0mxt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0myt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: object of type 'generator' has no len()"
     ]
    }
   ],
   "source": [
    "net2 = LogReg()#.cuda()\n",
    "loss=nn.NLLLoss()\n",
    "learning_rate = 1e-2\n",
    "optimizer=optim.SGD(net2.parameters(), lr=learning_rate)\n",
    "\n",
    "for epoch in range(1):\n",
    "    losses=[]\n",
    "    dl = iter(md.trn_dl)\n",
    "    for t in range(len(dl)):\n",
    "        # Forward pass: compute predicted y and loss by passing x to the model.\n",
    "        xt, yt = next(dl)\n",
    "        y_pred = net2(V(xt))\n",
    "        l = loss(y_pred, V(yt))\n",
    "        losses.append(l)\n",
    "\n",
    "        # Before the backward pass, use the optimizer object to zero all of the\n",
    "        # gradients for the variables it will update (which are the learnable weights of the model)\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "        # Backward pass: compute gradient of the loss with respect to model parameters\n",
    "        l.backward()\n",
    "\n",
    "        # Calling the step function on an Optimizer makes an update to its parameters\n",
    "        optimizer.step()\n",
    "    \n",
    "    val_dl = iter(md.val_dl)\n",
    "    val_scores = [score(*next(val_dl)) for i in range(len(val_dl))]\n",
    "    print(np.mean(val_scores))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stochastic Gradient Descent"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Nearly all of deep learning is powered by one very important algorithm: stochastic gradient descent (SGD). SGD can be seeing as an approximation of gradient descent (GD). In GD you have to run through all the samples in your training set to do a single itaration. In SGD you use only a subset of training samples to do the update for a parameter in a particular iteration. The subset used in each iteration is called a batch or minibatch.\n",
    "\n",
    "Now, instead of using the optimizer, we will do the optimization ourselves!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "object of type 'generator' has no len()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-128-8050ee202e7e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      7\u001b[0m     \u001b[0mlosses\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m     \u001b[0mdl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0miter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrn_dl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m     \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m         \u001b[0mxt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0myt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m         \u001b[0my_pred\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnet2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mV\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: object of type 'generator' has no len()"
     ]
    }
   ],
   "source": [
    "net2 = LogReg()#.cuda()\n",
    "loss_fn=nn.NLLLoss()\n",
    "lr = 1e-2\n",
    "w,b = net2.l1_w,net2.l1_b\n",
    "\n",
    "for epoch in range(1):\n",
    "    losses=[]\n",
    "    dl = iter(md.trn_dl)\n",
    "    for t in range(len(dl)):\n",
    "        xt, yt = next(dl)\n",
    "        y_pred = net2(V(xt))\n",
    "        l = loss(y_pred, Variable(yt))#.cuda())\n",
    "        losses.append(loss)\n",
    "\n",
    "        # Backward pass: compute gradient of the loss with respect to model parameters\n",
    "        l.backward()\n",
    "        w.data -= w.grad.data * lr\n",
    "        b.data -= b.grad.data * lr\n",
    "        \n",
    "        w.grad.data.zero_()\n",
    "        b.grad.data.zero_()   \n",
    "\n",
    "    val_dl = iter(md.val_dl)\n",
    "    val_scores = [score(*next(val_dl)) for i in range(len(val_dl))]\n",
    "    print(np.mean(val_scores))"
   ]
  },
  {
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}