{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Using SGD on MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### ... about machine learning (a reminder from lesson 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "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",
    "1. 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",
    "2. 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": {
    "hidden": true
   },
   "source": [
    "### About The Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "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": {
    "hidden": true
   },
   "source": [
    "<img src=\"images/mnist.png\" alt=\"\" style=\"width: 60%\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "A matrix can represent an image, by creating a grid where each entry corresponds to a different pixel.\n",
    "\n",
    "<img src=\"images/digit.gif\" alt=\"digit\" style=\"width: 55%\"/>\n",
    "  (Source: [Adam Geitgey\n",
    "](https://medium.com/@ageitgey/machine-learning-is-fun-part-3-deep-learning-and-convolutional-neural-networks-f40359318721))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports and data"
   ]
  },
  {
   "cell_type": "markdown",
   "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](https://kb.iu.edu/d/abbe) (often called *symlinks*) to make it possible to import these from your current directory.  For instance, I ran:\n",
    "\n",
    "    ln -s ../../fastai\n",
    "    \n",
    "in the terminal, within the directory I'm working in, `home/fastai/courses/ml1`."
   ]
  },
  {
   "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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/paperspace/anaconda3/envs/fastai/lib/python3.6/site-packages/sklearn/ensemble/weight_boosting.py:29: DeprecationWarning: numpy.core.umath_tests is an internal NumPy module and should not be imported. It will be removed in a future NumPy release.\n",
      "  from numpy.core.umath_tests import inner1d\n"
     ]
    }
   ],
   "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": "markdown",
   "metadata": {},
   "source": [
    "Let's download, unzip, and format the data."
   ]
  },
  {
   "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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mnist.pkl.gz\r\n"
     ]
    }
   ],
   "source": [
    "!ls {path}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": 7,
   "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": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(numpy.ndarray, (50000, 784), numpy.ndarray, (50000,))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(x), x.shape, type(y), y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalize"
   ]
  },
  {
   "cell_type": "markdown",
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.13044983, 0.3072898, -3.1638146e-07, 0.99999934)"
      ]
     },
     "execution_count": 9,
     "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": "markdown",
   "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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.005850922, 0.99243325)"
      ]
     },
     "execution_count": 10,
     "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": "markdown",
   "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": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show(img, title=None):\n",
    "    plt.imshow(img, cmap=\"gray\")\n",
    "    if title is not None: \n",
    "        plt.title(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000, 784)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_valid.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000, 28, 28)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_imgs = np.reshape(x_valid, (-1,28,28)); x_imgs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAEICAYAAACQ6CLfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADl1JREFUeJzt3X3IXHV6xvHr8mUxMaLRVE1iNLtPhL5htAapKEVd3NqtEFdwMWBJoyVrWaGrVSpBUBTBlu5qK1SJGMzirltN3FVWRcXa+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": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(x_imgs[0], y_valid[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000,)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_valid.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's the digit 3!  And that's stored in the y value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_valid[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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,  2.34669,  2.80432,  2.32126, -0.05587, -0.42452],\n",
       "       [-0.42452, -0.42452, -0.42452,  0.78312,  2.43567,  2.80432,  2.68991,  0.40176, -0.42452, -0.42452],\n",
       "       [-0.42452, -0.27197,  1.20261,  2.77889,  2.80432,  2.5755 ,  0.08396, -0.42452, -0.42452, -0.42452],\n",
       "       [-0.42452,  1.76194,  2.80432,  2.80432,  1.73651,  0.31278, -0.42452, -0.42452, -0.42452, -0.42452],\n",
       "       [-0.42452,  2.20685,  2.80432,  2.80432,  0.40176, -0.42452, -0.42452, -0.42452, -0.42452, -0.42452],\n",
       "       [-0.42452,  1.31702,  2.80432,  2.80432,  2.76618,  1.43143, -0.09401, -0.42452, -0.42452, -0.42452],\n",
       "       [-0.42452, -0.31011,  1.77465,  2.42296,  2.80432,  2.80432,  2.49923,  0.47803, -0.42452, -0.42452],\n",
       "       [-0.42452, -0.42452, -0.32282, -0.27197,  2.80432,  2.80432,  2.80432,  2.70262,  0.89752, -0.42452],\n",
       "       [-0.42452, -0.42452, -0.42452, -0.42452,  0.16023,  1.97804,  2.80432,  2.80432,  2.42296, -0.42452],\n",
       "       [-0.42452, -0.42452, -0.42452, -0.42452, -0.42452, -0.20841,  1.80007,  2.80432,  2.80432, -0.10672]],\n",
       "      dtype=float32)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_imgs[0,10:20,10:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAACt9JREFUeJzt3V+IVnUex/HPp5nKUTcSdm/USc3+rURiDVENBKUXtUVC7IVBkd1IsJVJGLUEXXTRTZRdlDFYEjTUxeRFRKRCdbEE0qRB2ViIteOkkhJbEsooffdingVrdZ6jc36deb69XxA40+nbF5m35zzPHM84IgQgp/OaXgBAOQQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGLdJYba5va4Qnp6eorMXbx4cZG5F154YZG5x48fr33m119/XftMSTpx4kSRuRHhdscUCRzlXHHFFUXmDg0NFZl72WWXFZlbIsYVK1bUPlOS9u/fX2RuFVyiA4kROJAYgQOJETiQGIEDiRE4kFilwG3fZvsr23ttP1F6KQD1aBu47S5JL0m6XdISSffYXlJ6MQBTV+UMfr2kvRGxLyLGJb0laWXZtQDUoUrg8ySdeivOWOtzv2J7je1h28N1LQdgaqrcqnq6+13/717ziBiQNCBxLzowXVQ5g49J6j3l4/mSDpRZB0CdqgT+iaTLbS+yfYGkVZLeKbsWgDq0vUSPiJO2H5K0VVKXpNciYnfxzQBMWaW/LhoR70l6r/AuAGrGnWxAYgQOJEbgQGIEDiRG4EBiPHSxkPvvv7/I3GeeeabI3N7e3vYHnYNjx44Vmfvggw/WPrPJhyOWwhkcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiso56qOnv27CJz169fX/vMp556qvaZknTeeWX+TP7hhx+KzO3v7y8yd8+ePUXmZsMZHEiMwIHECBxIjMCBxAgcSIzAgcTaBm671/aHtkds77a99vdYDMDUVfk++ElJj0XETtt/kvSp7e0R8WXh3QBMUdszeEQcjIidrV8flTQiaV7pxQBM3Vm9Bre9UNIySTtKLAOgXpVvVbU9W9Lbkh6NiJ9O8+/XSFpT424ApqhS4LbP10TcgxGx5XTHRMSApIHW8VHbhgDOWZV30S3pVUkjEfF8+ZUA1KXKa/B+SfdJutX2Z61//lZ4LwA1aHuJHhH/kuTfYRcANeNONiAxAgcSI3AgMQIHEiNwILGOeuji66+/XmTu3XffXWRuCUNDQ0XmvvDCC0Xm8nDEZnEGBxIjcCAxAgcSI3AgMQIHEiNwIDECBxIjcCAxAgcSI3AgMQIHEiNwIDECBxIjcCAxAgcSI3AgMQIHEiNwIDECBxIjcCAxAgcS66inqi5evLjpFRr38ssvF5n78ccfF5mLZnEGBxIjcCAxAgcSI3AgMQIHEiNwIDECBxKrHLjtLtu7bL9bciEA9TmbM/haSSOlFgFQv0qB254v6Q5Jm8quA6BOVc/gGyQ9LumXMx1ge43tYdvDtWwGYMraBm77TknfR8Snkx0XEQMR0RcRfbVtB2BKqpzB+yXdZftbSW9JutX2G0W3AlCLtoFHxJMRMT8iFkpaJemDiLi3+GYApozvgwOJndXfB4+IjyR9VGQTALXjDA4kRuBAYgQOJEbgQGIEDiTWUU9V3bZtW5G5S5cuLTK3hO3btxeZu3HjxiJzn3322SJzDxw4UGRuNpzBgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHEHBH1D7XrHyqpp6enxFgNDg7WPvO6666rfaYkXXLJJUXmlnLo0KEic1evXl37zK1bt9Y+s6SIcLtjOIMDiRE4kBiBA4kROJAYgQOJETiQWKXAbV9se8j2Htsjtm8svRiAqav600VflPR+RPzd9gWSZhbcCUBN2gZu+yJJN0taLUkRMS5pvOxaAOpQ5RL9UkmHJW22vcv2JtuzCu8FoAZVAu+WdK2kjRGxTNLPkp747UG219getj1c844AzlGVwMckjUXEjtbHQ5oI/lciYiAi+iKir84FAZy7toFHxCFJ+21f2frUcklfFt0KQC2qvov+sKTB1jvo+yQ9UG4lAHWpFHhEfCaJS2+gw3AnG5AYgQOJETiQGIEDiRE4kBiBA4l11FNVS5kxY0btM7u7q95icHaOHj1aZG6nOX78eO0z161bV/tMSXrllVeKzOWpqsAfHIEDiRE4kBiBA4kROJAYgQOJETiQGIEDiRE4kBiBA4kROJAYgQOJETiQGIEDiRE4kBiBA4kROJAYgQOJETiQGIEDifHQxQ5zzTXXFJm7YcOGInNvueWWInNLGB0dLTJ3wYIFReby0EXgD47AgcQIHEiMwIHECBxIjMCBxAgcSKxS4LbX2d5t+wvbb9qu/6f1Aahd28Btz5P0iKS+iLhaUpekVaUXAzB1VS/RuyX12O6WNFPSgXIrAahL28Aj4jtJz0kalXRQ0o8Rse23x9leY3vY9nD9awI4F1Uu0edIWilpkaS5kmbZvve3x0XEQET0RURf/WsCOBdVLtFXSPomIg5HxAlJWyTdVHYtAHWoEviopBtsz7RtScsljZRdC0AdqrwG3yFpSNJOSZ+3/puBwnsBqEF3lYMi4mlJTxfeBUDNuJMNSIzAgcQIHEiMwIHECBxIjKeqQpI0Z86cInM3b95cZO7KlSuLzC1h7ty5tc88cuSIxsfHeaoq8EdG4EBiBA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4kVuqpqocl/bvCoX+WdKT2BcrppH07aVeps/adDrsuiIi/tDuoSOBV2R6OiL7GFjhLnbRvJ+0qdda+nbQrl+hAYgQOJNZ04AMN///PVift20m7Sp21b8fs2uhrcABlNX0GB1BQY4Hbvs32V7b32n6iqT3asd1r+0PbI7Z3217b9E5V2O6yvcv2u03vMhnbF9sesr2n9Xt8Y9M7Tcb2utbXwRe237Q9o+mdJtNI4La7JL0k6XZJSyTdY3tJE7tUcFLSYxHxV0k3SPrHNN71VGsljTS9RAUvSno/Iq6StFTTeGfb8yQ9IqkvIq6W1CVpVbNbTa6pM/j1kvZGxL6IGJf0lqRp+fNgI+JgROxs/fqoJr4A5zW71eRsz5d0h6RNTe8yGdsXSbpZ0quSFBHjEfGfZrdqq1tSj+1uSTMlHWh4n0k1Ffg8SftP+XhM0zwaSbK9UNIySTua3aStDZIel/RL04u0camkw5I2t15ObLI9q+mlziQivpP0nKRRSQcl/RgR25rdanJNBX66H1w+rd/Otz1b0tuSHo2In5re50xs3ynp+4j4tOldKuiWdK2kjRGxTNLPkqbz+zFzNHGluUjSXEmzbN/b7FaTayrwMUm9p3w8X9P4Usf2+ZqIezAitjS9Txv9ku6y/a0mXvrcavuNZlc6ozFJYxHxvyuiIU0EP12tkPRNRByOiBOStki6qeGdJtVU4J9Iutz2ItsXaOKNinca2mVStq2J14gjEfF80/u0ExFPRsT8iFioid/XDyJiWp5lIuKQpP22r2x9armkLxtcqZ1RSTfYntn6uliuafymoDRxifS7i4iTth+StFUT70S+FhG7m9ilgn5J90n63PZnrc/9MyLea3CnTB6WNNj6g36fpAca3ueMImKH7SFJOzXx3ZVdmuZ3tXEnG5AYd7IBiRE4kBiBA4kROJAYgQOJETiQGIEDiRE4kNh/Ac2gXO9TZq9jAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(x_imgs[0,10:20,10:20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAF0CAYAAADGqzQSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xm8ndPd///3R0iVIEITESJISkOJiqFEb0OJqEcirbZS3MQQlFu5dYiguPGLTiJFozElao4kRFEltDHcaISm0ZjahpAJEUlMmT7fP/bO706z1sm5ztnj2tfr+XjkkXPe2Xtfa5/zPjsrO9e6lrm7AAAAgDxar9YDAAAAAGqFyTAAAAByi8kwAAAAcovJMAAAAHKLyTAAAAByi8kwAAAAcovJcI2ZWV8ze9zM5pnZZ2b2tpndY2Y9az02YF3MbH8z+6OZLTCzxWY2zcxOqvW4gCzM7Agzm2JmS4v9nWpmB9d6XMC6mNlBZvaUmX1iZgvN7Hdm1qnW40odk+Ha6yDpBUlnSTpM0vmSdpH0rJltV8uBAU0xs90kPSZpA0mnSvqWpL9IusnMzqjl2IDmmNlpku5X4bV3oKRvSxonaaNajgtYFzM7QNIfJS1S4TX3B5K+JmmymX2ulmNLnbHpRv0xs50kvSLph+7+q1qPB1ibmf1/kn4oqYO7L10jf1aSu/tXazY4YB3MrJukmZLOd/erazsaIDsze0xSN0k7u/uKYraXpOclnenuv6nh8JLGO8P16f3i78trOgqgaW1V6Ocna+WLxOsK6ttJklZJur7WAwFaaF9Jj66eCEuSu/9FhTnDwJqNqgHwl1adMLM2ZtbWzHpI+q2keZLuqvGwgKaMKf7+azPb2szam9mpkg6RNKJ2wwKa1UeF/3k7xsz+YWYrzOwNMzuz1gMDmrFS0rJI/pmkXas8lobCaRJ1wsymStqz+Okbkvq7+8waDglYp+J/z02U1KUYLZd0hrvfVLtRAetmZq9I2lqFCcQwSf9Q4Zzh0yWd4+4jazg8oElm9rwKp6Hts0a2naR/SVru7pw33EpMhuuEmX1J0qaSdlDhXMxOkvq4+6xajguIKf4PxmQVzr28RoXTJQZIOkPSie5+ew2HBzTJzF6T1EPSt9x9whr5w5L2kNTZ+YsRdcjMjpV0m6QrJP1ahQX4oyXtp8Jk+PM1HF7SmAzXITNrL2mWpLvc/fQaDwcImNk4SV9RYSHH8jXy2yX1ldTR3VfVanxAU8zsf1U493JTd1+yRn6upKskdXH3ObUaH7AuZnaZCm+YbSjJJd0taWNJu7r7DrUcW8o4Z7gOufsiFU6V6F7rsQBN+LKkv645ES56XtIWkjpWf0hAJi83kVvxd/4Rh7rl7hdJ2lLSbir8L8YgFf6n46maDixxTIbrUPEC2jurcC4bUI/mSeplZm3XyveR9KmkhdUfEpDJxOLvfdfK+0p6293nVXk8QIu4+0fu/jd3n29mh6swX+DqKCVYv9YDyDszmyhpmqTpkhZL+qKkcyWtkMQ1hlGvrlVhk4IHzOw3Kpwz3F/SIEkj3D224hmoBw9JekLSb81sS0n/lHS0CpseDa7lwIB1MbM9JPVTYc4gFa6M8iNJP3f3Z2o2sAbAOcM1ZmY/kfQdSTuqcO3W2ZL+JGk4i+dQz8ysn6SfqLBj4oYq/E/GaEm/dfeVtRwbsC5mtqmk4SpMgjdX4VJrV7r7HTUdGLAOZraLCpde3VXS51RcwOzut9R0YA2AyTAAAAByi3OGAQAAkFtMhgEAAJBbTIYBAACQW0yGAQAAkFslXVqteH27kZLaSLrR3a9s5vas1kM5vOfuXyjlAeguasHdrflbrVtLuktvUSa85iJVmbrb6neGzayNpOtUuOZdT0mDzKxnax8PaIE3S7kz3UWq6C5qhNdcpCpTd0s5TWJvSW+4+z+LF9i/S9KAEh4PqBa6i1TRXaSI3qKulTIZ7qLCBhGrvV3M/o2ZDTGzqWY2tYRjAeVEd5GqZrtLb1GHeM1FXSvlnOHYuW/BOT7uPlqFXak4Bwj1gu4iVc12l96iDvGai7pWyjvDb0vado3Pt5E0p7ThAFVBd5EquosU0VvUtVImw3+R1MPMtjeztpKOkTSpPMMCKoruIlV0Fymit6hrrT5Nwt1XmNlZkh5R4VIpN7v7y2UbGVAhdBeportIEb1FvTP36p2WwzlAKJMX3L13NQ9Id1EO5bjOcEvQW5QJr7lIVabusgMdAAAAcqukHegAAEhZu3btguzkk08OsgED4pfF7d+/f5AtXbq09IEBqBreGQYAAEBuMRkGAABAbjEZBgAAQG4xGQYAAEBuMRkGAABAbnE1CQBAbp1wwglBNmLEiMz332WXXYLsueeeK2lMAKqLd4YBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQWyygK8Huu+8eZOeee270tjvuuGOQbbTRRkE2bNiwINtss82C7OGHH44eZ8mSJdEcAPLuxBNPDLKrr746yJYvXx5kv/zlL6OPOW3atJLHBaC2eGcYAAAAucVkGAAAALnFZBgAAAC5xWQYAAAAuWXu3vo7m82StETSSkkr3L13M7dv/cFqrF27dkH21ltvBVn79u2rMRy988470Ty2gO/ee++t9HCq7YXmutacPHU3JtbTgQMHRm+7xx57BFmfPn2CLPYzsnDhwiDbaqutoseZN29ekI0ZMybIbrjhhiBbuXJl9DHrjbtbqY/Rku42Wm9bon///kE2ceLEIPv444+D7Kc//WmQtWRXugbEay5Slam75biaxEHu/l4ZHgeoNrqLVNFdpIjeoi5xmgQAAAByq9TJsEv6o5m9YGZDYjcwsyFmNtXMppZ4LKCc6C5Stc7u0lvUKV5zUbdKPU1if3efY2YdJT1qZq+4+5Q1b+DuoyWNljgHCHWF7iJV6+wuvUWd4jUXdaukybC7zyn+vsDMJkraW9KUdd8rTWbhupeXX345yN5///3o/V988cUgiy1M2m677YJs2223DbIOHTpEj/Pzn/88yJ588skgmz9/fvT+eZGn7m6zzTZBdt999wVZrI9NWbx4cZDFOr7BBhsE2aJFi6KP2bVr1yC77rrrguyDDz4IsilTwm/d3Llzo8dJXZ66m0Xbtm2j+Xe/+90gi72OP/fcc0GW88VyFUFvUc9afZqEmW1sZpus/ljSYZJmlGtgQKXQXaSK7iJF9Bb1rpR3hjtJmlj8l/b6ku5w9z+UZVRAZdFdpIruIkX0FnWt1ZNhd/+npN3LOBagKuguUkV3kSJ6i3rHpdUAAACQW+XYdCMXlixZEmQHHHBAVY695ZZbBtmPfvSj6G1j+eGHHx5kY8eOLX1gSMKkSZOCbPfdwzdpHnvssej9zzvvvCB7773wuvmxHeRa4gtf+EKQPfzww0F24403BtnQoUODLLb4Do2nqdfCQYMGBdltt90WZCeddFLZxwSsrXPnzkH2/e9/P3rbWL58+fIgi+2Ce8UVVwTZtGnToseZPXt2NM8j3hkGAABAbjEZBgAAQG4xGQYAAEBuMRkGAABAbjEZBgAAQG5xNYkExFbuP/3009HbxlZWx7bZ5WoSjSm2YrlXr15Bds899wTZscceG33MlStXlj6wDN59990ge/XVV4Ms1uemfh7QWHr37h1kF110UfS2r732WpANHjw4yKrVb+THDjvsEGSjRo0KskMPPbSk43Ts2DHIJk6cGGSfffZZ9P77779/kDV15YlGxzvDAAAAyC0mwwAAAMgtJsMAAADILSbDAAAAyC0W0CVg8803D7Jhw4Zlvv/WW29dzuGgjsUWy5lZkM2ZMyfIar2QaN999w2y2Ja6TzzxRJDFnvdLL71UnoGhJtZbL3yvJrbtdtu2baP3f+CBB4Ks1h1H4+nSpUuQzZgxI8jWXz+cbo0YMSL6mNdcc02m4+y8885B9otf/CLI2rdvHz1ObCF17HU4toi/0fDOMAAAAHKLyTAAAAByi8kwAAAAcovJMAAAAHKr2QV0ZnazpCMlLXD3XYtZB0l3S+omaZak77j7B5UbZn7svvvuQTZu3Lgg6969e/T+sV2XzjvvvNIHlqA8dvfhhx8OMncPsu9973tBdvXVV0cf86233ip9YGvYZJNNovkNN9wQZPfdd1+QxXbKi+34lLI8dndtnTp1CrJvfvObme//5ptvlnM4yCCPvf3xj38cZG3atAmyU045JchuvfXWzMeZNWtWkMV23txwww2DrKnX9tjr5pQpU4Istqhu8eLF0cdMVZZ3hsdIOnytbKikye7eQ9Lk4udAvRkjuos0jRHdRXrGiN4iQc1Oht19iqSFa8UDJI0tfjxW0lFlHhdQMrqLVNFdpIjeIlWtvc5wJ3efK0nuPtfMOjZ1QzMbImlIK48DlBvdRaoydZfeos7wmou6V/FNN9x9tKTRkmRm4cmLQJ2iu0gRvUWq6C5qpbWT4flm1rn4r7zOkhaUc1B5ccIJJwTZ//zP/wTZtttuG2SffPJJ9DHPOOOMIJs9e3YrRtewctfdWKcuuuiiIPvDH/4QvX/fvn2DrJROjR8/Ppp/8YtfDLLYDnSx7r/88sutHk9CctXdww9f+9TTuMceeyyajxo1qpzDQes1RG833XTTaH7MMccEWWxnuZYslivFddddF2Rnn3129LY9evQIstiudhdffHGQNdrC/NZeWm2SpNUzuRMk3V+e4QAVR3eRKrqLFNFb1L1mJ8Nmdqek/5W0k5m9bWYnS7pS0qFm9rqkQ4ufA3WF7iJVdBcpordIVbOnSbh7+P+UBYeUeSxAWdFdpIruIkX0FqliBzoAAADkFpNhAAAA5JbFtmqt2MFycKmUdu3aRfMf/vCHQXbhhRcG2Xrrhf8+Wbhw7WuYS3369Ike55VXXmluiI3gBXfvXc0Dptzd2PacY8eODbKjjz46ev833ngjyA488MAgmzt3bpD95je/CbIhQ+KXEf3Rj34UZLFV2Slzd6vm8VLp7frrh2fszZw5M8i22267INt+++2jj/nOO++UPjCslvvX3L333juaP/vss0F26KGHBtnkyZPLPqasBg4cGM0nTJgQZLE54aJFi4IsdiWK999/vxWjq7hM3eWdYQAAAOQWk2EAAADkFpNhAAAA5BaTYQAAAORWa7djRhPGjBkTzb/5zW9muv+9994bZFdffXWQ5WShHMrg008/DbJTTjklyDp27Bi9/3/8x38E2Z///OcgGzduXJAdd9xxQdbUdsyNtlgO2cUWb+64445BFttuvtYL5WLbRvfv3z/IYtud//GPfwyy2M8ram+PPfbIfNsXX3yxgiNpuYceeiiaxxZHx37uYp386KOPSh9YHeGdYQAAAOQWk2EAAADkFpNhAAAA5BaTYQAAAOQWC+jKLHbyeUuMGjUqyJ555pmSHhNY25IlS4JswIAB0dtecsklQXbOOecE2dChQzMd+5prrsl0O+RH165dM92ubdu2FR5J00488cRoHttlMbbr4+mnnx5ksZ297rvvvuhxTjrppGZGiEp66qmnovmqVauC7NFHHw2yI488Mshiu3ZWwk477RTNYz3t27dvkG200UZB1mgLPXlnGAAAALnFZBgAAAC5xWQYAAAAucVkGAAAALnV7AI6M7tZ0pGSFrj7rsXsEkmnSnq3eLNh7h7f4iRnYjsKSdLuu+/e6vvHFtVdeeWV0fvPmTMn03HygO62zOLFi6P5T3/60yA79NBDg6xnz56ZjvP1r389mje1QCWP8tbd7t27Z7pdtXbebN++fZBdddVV0dvGFiGtWLEiyGKLqvr06RNksV0bpTQW0DVyb19++eVo/vvf/z7IYouRZ86cGWSxXQml+C6djz/+eJB16dIlyGKL5WK72EpS586dgyzW3fvvvz96/0aS5Z3hMZLC/SalEe7eq/gruWIjF8aI7iJNY0R3kZ4xordIULOTYXefImlhFcYClBXdRaroLlJEb5GqUs4ZPsvMppvZzWa2eVM3MrMhZjbVzKaWcCygnOguUtVsd+kt6hCvuahrrZ0Mj5K0o6RekuZK+lVTN3T30e7e2917t/JYQDnRXaQqU3fpLeoMr7moe63agc7d56/+2MxukBSeQZ5Tsd26pPgCkT333DPIYjsxnXXWWUF29NFHR48zePDgIHvkkUeit80juttyBxxwQJD16NGj1Y/3k5/8JJq/+eabQXbLLbe0+jiNppG7G1sING/evCCLLUKrhNhuc7FFdZJ02223BdnIkSOD7K233gqy2AKqL3/5yxlGmI5G7q0kDRo0KMiGDx8eZGeffXaQfec734k+ZixfuDA8+6RDhw5ZhtgilXjMFLTqnWEzW3MJ4kBJM8ozHKCy6C5SRXeRInqLFGS5tNqdkg6UtKWZvS3pYkkHmlkvSS5plqTTKjhGoFXoLlJFd5EieotUNTsZdvfw/wCkmyowFqCs6C5SRXeRInqLVLEDHQAAAHKLyTAAAAByq1VXk0DTPvnkk2h+7LHHBtn664df/qa2xF3bVlttFc0nTpwYZP/93/8dZNdff32m4wAHHXRQkLl7kA0cODDIYiugY9uXSvFtx997770ge+CBB6L3R7r22WefIFu2bFkNRtJyc+bMCbJtttkmyEaPHh1kX/nKV4KMq/+kJfZ3/jnnnBNk99xzT5DF5gVN6dSpU6bbLV++PMhiP1+StP322wfZxx9/nHlMjYR3hgEAAJBbTIYBAACQW0yGAQAAkFtMhgEAAJBbLKCrkk8//TTT7Xr16hVkI0aMCLLYoiZJ2nDDDYNs6NChQcYCOqxtt912i+Y/+MEPgiy22G3SpEmZjnPqqadG85tuCi9HescddwTZLrvsEmSxrW6RjgkTJgTZkUceWZVjm1mmrCk//vGPM90utuj02muvDbJhw4ZlPjbS8cwzz2TKKuF3v/tdNO/WrVuQrVy5ssKjqU+8MwwAAIDcYjIMAACA3GIyDAAAgNxiMgwAAIDcYgFdRhtttFGQVWKnlunTpwfZ0UcfHWQ333xz9P4DBgwIsq5duwZZ586dg2zu3LlZhogGtckmm0Tz2E6J9957b6uPM27cuGi+3XbbBdnPfvazINtzzz2DjAV0jad9+/ZBFlsIdNttt0XvH+vtMcccE2QdOnQIsn79+mUZoiTpo48+CrKnnnoqyH7+858H2RNPPJH5OEA17LjjjrUeQk3wzjAAAAByi8kwAAAAcovJMAAAAHKLyTAAAAByq9kFdGa2raRbJW0laZWk0e4+0sw6SLpbUjdJsyR9x90/qNxQqyd2AnlsQcSDDz4YZDNmzIg+Zmxx2sknnxxkG2ywQZB16dIlyLp37x49Tsw//vGPTONpNHnsbiliux9K0rx584Is9vNQqthuXLHd6s4888wgmzhxYtnHUyt57O2LL74YZKecckqQHXvssZmyUi1evDjImlr4efnllwfZm2++WfYxpSCP3U3B0qVLaz2EupflneEVks5z9y9J2lfSmWbWU9JQSZPdvYekycXPgXpCd5EieotU0V0kqdnJsLvPdfdpxY+XSJopqYukAZLGFm82VtJRlRok0Bp0Fymit0gV3UWqWnSdYTPrJmkPSc9J6uTuc6XCD4CZdWziPkMkDSltmEBp6C5SRG+RKrqLlGSeDJtZO0njJZ3j7ovNLNP93H20pNHFx/DWDBIoBd1FiugtUkV3kZpMk2Ez20CFYt/u7hOK8Xwz61z8V15nSQsqNchq+/a3vx1kW221VZCddNJJZT927EXDPftrQuxE+dNPP72kMaUsb90tRWynQkl6/vnnq3L8ZcuWBdkHH4RrbA444IAgi+0itnDhwvIMrAby1ts77rgjyGI7b77++utB1qZNm+hjNpWv7fbbbw+yWbNmBVlsITJCeetuCqZMmRLNTzvttCDr2DH6pn3Da/acYSvMzm6SNNPdr1rjjyZJOqH48QmS7i//8IDWo7tIEb1FquguUpXlneH9JR0v6W9m9lIxGybpSkn3mNnJkt6SFL6dCtQW3UWK6C1SRXeRpGYnw+7+lKSmTvg5pLzDAcqH7iJF9BaportIFTvQAQAAILeYDAMAACC3WnSd4bzYYostaj2EfzN+/Pggu+yyy6K3XbAgXKQb204XWFtTVy3p06dPkB1zzDFB9vjjjwdZu3btgqxt27bR4+y8885BttdeewXZddddF2QpXzkC0ocffhhkhxzC/6oD5bDeevH3PWNXr4rNIfKAd4YBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQWyygixg2bFiQPfbYY0F23HHHBdnWW28dfczYApGYa665JsiefPLJIFuxYkWmxwOymjlzZjSPbXUc2z73/fffD7KWLKCLLeZ4+umng+ySSy6J3h8AEFq1alU0b2rRdB7xzjAAAAByi8kwAAAAcovJMAAAAHKLyTAAAAByiwV0EcuXLw+yRx55JFMGpOoPf/hDNL/22muDLLYrXa9evUo6/gUXXBBkN998c5Cx2xwAVMZhhx0WZKNGjarBSKqLd4YBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQW80uoDOzbSXdKmkrSaskjXb3kWZ2iaRTJb1bvOkwd3+oUgMFWorutsz8+fOj+Q9+8IMqjyTf6C1SRXfr09KlSzPfdv3183ldhSzPeoWk89x9mpltIukFM3u0+Gcj3P2XlRseUBK6ixTRW6SK7iJJzU6G3X2upLnFj5eY2UxJXSo9MKBUdBcpordIFd1Fqlp0zrCZdZO0h6TnitFZZjbdzG42s82buM8QM5tqZlNLGilQArqLFNFbpIruIiWZJ8Nm1k7SeEnnuPtiSaMk7Siplwr/EvxV7H7uPtrde7t77zKMF2gxuosU0Vukiu4iNZkmw2a2gQrFvt3dJ0iSu89395XuvkrSDZL2rtwwgdahu0gRvUWq6C5SlOVqEibpJkkz3f2qNfLOxfODJGmgpBmVGSLQOnQXKaK3SBXdrU9PPvlk5tsefPDBFRxJ/cpyNYn9JR0v6W9m9lIxGyZpkJn1kuSSZkk6rSIjBFqP7iJF9BaportIkrl79Q5mVr2DoZG9UO1zyuguysHdrZrHo7coE15zE9a+fftovnDhwiD75JNPgmzjjTcu+5iqKFN32YEOAAAAucVkGAAAALmVz333AAAAcmDRokXRfL31eD90Nb4SAAAAyC0mwwAAAMgtJsMAAADILSbDAAAAyK1qL6B7T9KbxY+3LH7eKHg+1bNdDY65urv1/HVpDZ5P9dSyt1J9f21ag+dTPXS3vHg+1ZOpu1XddOPfDmw2tdoX8a4knk8+NNrXheeTH432teH55EejfW14PvWH0yQAAACQW0yGAQAAkFu1nAyPruGxK4Hnkw+N9nXh+eRHo31teD750WhfG55PnanZOcMAAABArXGaBAAAAHKLyTAAAAByq+qTYTM73MxeNbM3zGxotY9fKjO72cwWmNmMNbIOZvaomb1e/H3zWo6xJcxsWzN7wsxmmtnLZvaDYp7sc6oUultf6G42qfdWaqzu0tvsUu9uI/VWauzuVnUybGZtJF0nqZ+knpIGmVnPao6hDMZIOnytbKikye7eQ9Lk4uepWCHpPHf/kqR9JZ1Z/J6k/JzKju7WJbrbjAbprdRY3aW3GTRId8eocXorNXB3q/3O8N6S3nD3f7r7Mkl3SRpQ5TGUxN2nSFq4VjxA0tjix2MlHVXVQZXA3ee6+7Tix0skzZTURQk/pwqhu3WG7maSfG+lxuouvc0s+e42Um+lxu5utSfDXSTNXuPzt4tZ6jq5+1ypUBZJHWs8nlYxs26S9pD0nBrkOZUR3a1jdLdJjdpbqQG+z/R2nRq1uw3xfW607lZ7MmyRjGu71QEzaydpvKRz3H1xrcdTh+hunaK760Rv6xS9bRbdrVON2N1qT4bflrTtGp9vI2lOlcdQCfPNrLMkFX9fUOPxtIiZbaBCsW939wnFOOnnVAF0tw7R3WY1am+lhL/P9DaTRu1u0t/nRu1utSfDf5HUw8y2N7O2ko6RNKnKY6iESZJOKH58gqT7aziWFjEzk3STpJnuftUaf5Tsc6oQultn6G4mjdpbKdHvM73NrFG7m+z3uaG76+5V/SXpCEmvSfqHpAuqffwyjP9OSXMlLVfhX64nS9pChRWUrxd/79DKr8sUSUslLZY0VdLBVXg+fVT4r6fpkl4q/jqiHM+p0X7R3WYf/w/FLl1epedDd7N9nZLubfE5lK27kg4s9mbtX4uq9FzobfavVdLdLfdrrqS+kh6XNE/SZ8XHvEdSzyo9n4btLtsx1wEzO03StcVfD6nwjn0vSS+7++9rOTYgCzMbJOkqSVtJusLdL6zxkIAoMztQ0hOSzlbh3cfVVrj71JoMCsig+Dr7FRUWrb0rqasKlzHbVtKX3f3NGg4vaevXegB5V1yRebWkH7n71Wv80SM1GRDQQmbWXtIISedKuqPGwwGymunuz9Z6EEBW7n6nCu82///M7HlJr0g6WtKvajGuRsB2zLV3kqRVkq6v9UCAVvq5Cv+LcWeztwQAlNP7xd+X13QUiWMyXHt9VPhX3TFm9g8zW1HcevLMWg8MaI6Z9ZH0n5K+X+uxAC10u5mtNLP3zewOM+ta6wEBWZhZGzNra2Y9JP1WhXOI76rxsJLGaRK1t3Xx1y8kDVNhocC3JV1rZuu7+8haDg5oSvESO7+V9Et3f7XW4wEy+lCF/07+swqLlfdQ4bX3f81sD3dP7rJQyJ3nJO1Z/PgNFRbb09sSsICuxszsNUk9JH3L/++afTKzh1V4ke7sfJNQh8zsQhVO89nF3T8pZi4W0CExZvYVSc9LupLuot6Z2ZckbSppB0k/lNRJUh93n1XLcaWM0yRqb/X5Po+ulf9RhYJ3ru5wgOYV/0v5AkkXSfqcmbUvLqTTGp+3qd0IgezcfZoKl/Daq9ZjAZrj7jPd/bniOo1DJLVT4aoSaCUmw7X3chP56q0oV1VrIEAL7CBpQ0m3SfpgjV9S4Z2KDyR9uTZDA1rFxHa/SIy7L1LhVInutR5LypgM197E4u9918r7Snrb3edVeTxAFi9JOijySypMkA9S4QUaqHtm1lvSF1U4FxNIhpl1krSzCuuN0EosoKu9h1S4APxvzWxLSf9U4Xp1Fe1TAAAaI0lEQVSBh0kaXMuBAU0pvhvxp7Xzwm6detPdgz8D6oGZ3S7pX5KmSVqkwtqM8yW9I+maGg4NWCczm6hCb6ersPjziypc332FuMZwSZgM15i7u5kdJWm4pEslba7CpdaOdXc2MACA8pohaZCk/5K0kQqXpZog6WJ3f6+WAwOa8ayk70g6T1JbSbNVeFNiOIvnSsPVJAAAAJBbnDMMAACA3GIyDAAAgNxiMgwAAIDcYjIMAACA3CrpahJmdrikkZLaSLrR3a9s5vas1kM5vOfuXyjlAeguasHdrflbrVtLuktvUSa85iJVmbrb6neGi1utXiepn6SekgaZWc/WPh7QAm+Wcme6i1TRXdQIr7lIVabulnKaxN6S3nD3f7r7Mkl3SRpQwuMB1UJ3kSq6ixTRW9S1UibDXVS44PNqbxezf2NmQ8xsqplNLeFYQDnRXaSq2e7SW9QhXnNR10o5Zzh27ltwjo+7j5Y0WuIcINQNuotUNdtdeos6xGsu6lop7wy/LWnbNT7fRtKc0oYDVAXdRaroLlJEb1HXSpkM/0VSDzPb3szaSjpG0qTyDAuoKLqLVNFdpIjeoq61+jQJd19hZmdJekSFS6Xc7O4vl21kQIXQXaSK7iJF9Bb1ztyrd1oO5wChTF5w997VPCDdRTmU4zrDLUFvUSa85iJVmbrLDnQAAADILSbDAAAAyC0mwwAAAMgtJsMAAADILSbDAAAAyC0mwwAAAMgtJsMAAADIrVZvuoGW6devX5Cde+65QXbooYcGWexa0K+//nr0OPfcc0+QjRo1KsjmzGEnTAAAAN4ZBgAAQG4xGQYAAEBuMRkGAABAbjEZBgAAQG5ZbHFWxQ5mVr2D1cgZZ5wRzUeMGBFkbdu2rfRwJElPPPFEkB133HFBNnfu3GoMpxxecPfe1TxgHrqLynN3q+bx6C3KhNfcFhg7dmyQHX/88UH24IMPRu8/fvz4IHvmmWeCbPbs2ZnGs2zZsmi+cuXKTPdPXKbu8s4wAAAAcovJMAAAAHKLyTAAAAByi8kwAAAAcqukHejMbJakJZJWSlpR7RPsa+0b3/hGkP3yl7+M3ja2WO7FF18MsqFDhwbZyy+/nHlMJ598cpBdeumlQXb++ecH2dlnn535OKnLe3c33njjIBs2bFj0thdeeGGQxRbeXnbZZUG2++67B1n//v2zDBFNyHt3kaY89faVV14JslWrVgVZbA6xrry1brnllmh+2mmnBdmKFSvKeuxUlGM75oPc/b0yPA5QbXQXqaK7SBG9RV3iNAkAAADkVqmTYZf0RzN7wcyGxG5gZkPMbKqZTS3xWEA50V2kap3dpbeoU7zmom6VeprE/u4+x8w6SnrUzF5x9ylr3sDdR0saLaV9EW00HLqLVK2zu/QWdYrXXNStkibD7j6n+PsCM5soaW9JU9Z9rzQdeeSRQXbnnXcG2ec///no/e+7774gi+1WN3/+/FaM7v9cfvnlQRbbWe6www4r6Tipy1N3Y7bYYosga2oBXWw3pGnTpmU6zte+9rUg69SpU/S2pXY/L/Le3UYS+1no3r179LYbbrhhkA0aNCjIbr/99iBrageyp59+urkhlk2eejt8+PAg+9vf/hZkffv2zfyYe+21V5B17do1yGJzkMGDB0cf87bbbguy2I61edDq0yTMbGMz22T1x5IOkzSjXAMDKoXuIlV0Fymit6h3pbwz3EnSRDNb/Th3uPsfyjIqoLLoLlJFd5Eieou61urJsLv/U1J4EVGgztFdpIruIkX0FvWOS6sBAAAgt8qx6UbDWX/98MsS28UttovX9OnTo48Z2+nl3XffbcXo1i22M9iNN94YZBMnTiz7sZGObt26lf0xly9fHmSbbbZZkPXs2TN6fxbQoVHsuuuuQfbd7343yE466aQg69y5c/QxY6/tMU0tlopp06ZN5tuiNL///e8zZaXq169fkD344IPR2x5xxBFBxgI6AAAAIGeYDAMAACC3mAwDAAAgt5gMAwAAILeYDAMAACC3uJpExKmnnhpke+yxR5B99tlnQXbiiSdGH7MSV44oxfvvv1/rIaCGvvrVr5b9Me+///4gi12FpXfv3tH753UVM9LQq1evaH7uuecG2de//vUg22qrrco+ppglS5YE2eOPP16VY6O6OnToEGQXX3xxkK1YsSJ6/6auMpFHvDMMAACA3GIyDAAAgNxiMgwAAIDcYjIMAACA3GIBXcR//dd/Zbrd6aefHmQvvfRSuYcDlCS25eq3vvWtIFu1alX0/k0tvgBaKrbVvSRtuOGGQbZ06dJKD0dSfEHnLbfcEmQ77rhj9P6f+9znyj6mmL///e9BduGFFwZZbHH0U089VZExoTSbbLJJNO/Tp0+QtW3bNsguuOCCIIv1+dZbb40e509/+lMzI8wP3hkGAABAbjEZBgAAQG4xGQYAAEBuMRkGAABAbjW7gM7MbpZ0pKQF7r5rMesg6W5J3STNkvQdd/+gcsOsT2+//Xath4B1oLsFnTp1CrK99toryP71r39F7z99+vRMx1m+fHmQrVy5Msi6d++e6fHyrFG7G9sdS5KOOuqoIBs/fnyQXXLJJZmPtdtuuwXZT37ykyCLLSbdYIMNgszMosdx98xjyiL2vCXpP//zP4Psk08+KeuxS9WovW2pdu3aBdnw4cODLNY9qbTdCp977rkgu/LKK1v9eHmR5Z3hMZIOXysbKmmyu/eQNLn4OVBvxojuIk1jRHeRnjGit0hQs5Nhd58iaeFa8QBJY4sfj5UU/rMeqDG6i1TRXaSI3iJVrb3OcCd3nytJ7j7XzDo2dUMzGyJpSCuPA5Qb3UWqMnWX3qLO8JqLulfxTTfcfbSk0ZJkZuU9uQqoILqLFNFbpIruolZaOxmeb2adi//K6yxpQTkHVS2xBRaS1KNHjyBbsmRJkL366qtlHxMqriG6Wwmvv/56Sfd/4403gmz27NlB1qtXr5KOk2NJdXfTTTcNsuOPPz56265duwbZLrvsEmSxhUk77bRT9DG/8Y1vNDfEFmlqAV1MbBe43/3ud0E2YcKEIGvA3eKS6m057L///kF25plnVuXYsZ+RpnYXxf9p7aXVJkk6ofjxCZLuL89wgIqju0gV3UWK6C3qXrOTYTO7U9L/StrJzN42s5MlXSnpUDN7XdKhxc+BukJ3kSq6ixTRW6Sq2dMk3H1QE390SJnHApQV3UWq6C5SRG+RKnagAwAAQG4xGQYAAEBuVfzSavVs/fXjT79NmzZB9vHHHwcZ2zEjBQcffHCm240YMaKk48R+nmI/S507d47eP3b1gcWLF5c0JtROhw4dgmzjjTeO3jbrlsbnnntukFVim+S//OUvQXb33XdHb/vQQw8F2dKlS4PsnXfeafV4kJY+ffqUdP8FC8ILbowaNSrI1lsvfD/zoosuCrLYVtCSdMoppwTZBx809E7ZTeKdYQAAAOQWk2EAAADkFpNhAAAA5BaTYQAAAORWrhfQ1doWW2wRZEceeWSQnXfeeZkfc9asWUHWrVu3IJs3b16Q3XvvvUF2yy23RI+zfPnyzGNCbe23335BNn/+/CB78sknSzpObJHpgw8+GGSnn3569P6bbbZZkLGALl2x16J33303etvYYrtqueyyy4Ls17/+dZAtXLiwGsNBA7j00kuD7IUXXgiyjz76KHr/P//5z0G2bNmyIIstHh03blyQTZ48OXqcG2+8MchOPvnkIFu0aFH0/o2Ed4YBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQWyygyyi2wKN3795BNnXq1Oj9u3fvHmSPPfZYkHXt2jXIPvnkkyD761//Gj1ObNFKLBs8eHCQff3rXw+yvn37Ro/zrW99K5qjtmI7fB1xxBFBFluM0dRijlLkYeEFsmtqIc9OO+3U6secMmVKNB8/fnyQ3XHHHUEW23Fr1apVrR4PsGLFiiC77777yn6c2C6LM2bMCLJTTz01ev+JEycG2RNPPBFk1157bStGlxbeGQYAAEBuMRkGAABAbjEZBgAAQG4xGQYAAEBuNbuAzsxulnSkpAXuvmsxu0TSqZJWbyc0zN0fqtQgK6WpHYU+/PDDIIvtjhXLdthhh+hjPv7440G2zTbbBFlsgcmZZ54ZZK+99lr0OFlNmjQpyGIn0++8884lHaeWGrm7Tdloo42CbLvttguy2bNnV2M40Z+lpsR+nqo1znrTqN09//zzo3ls583YYuKYAw88sJQhoYwatbepi/19L0l33XVXkMV+Ru++++4ga2o3yVRleWd4jKTDI/kId+9V/EWxUY/GiO4iTWNEd5GeMaK3SFCzk2F3nyKJTdmRHLqLVNFdpIjeIlWlnDN8lplNN7ObzWzzpm5kZkPMbKqZxS/AC1Qf3UWqmu0uvUUd4jUXda21k+FRknaU1EvSXEm/auqG7j7a3Xu7e7hDBVB9dBepytRdeos6w2su6l6rdqBz9/mrPzazGyT9vmwjqqLYzmySNHfu3CCLLe753ve+F2Q9e/aMPmZssVxsB7qBAwcGWSV2Bosd+8Ybbwyyww47rOzHrqVG6W6p2rZtG2R77rln9LaffvppkMUWn37+858PstgOSU0ZNWpUkB188MFBtnz58syP2UgaobtLly6N5rGFPMcdd1yQdenSJcjmzZsXfcxx48YF2cUXXxxkTS2kRnk0Qm8b1ciRI4Ns0KBBQTZkyJAgu+KKKyoyplpp1TvDZtZ5jU8HSgr3/wPqEN1FquguUkRvkYIsl1a7U9KBkrY0s7clXSzpQDPrJcklzZJ0WgXHCLQK3UWq6C5SRG+RqmYnw+4evmcu3VSBsQBlRXeRKrqLFNFbpIod6AAAAJBbTIYBAACQW626mkSju//++4Msti3x4MGDMz9m7EoN55xzTpB9/PHHmR+z3Dp16hRknTt3jtwyvlXqW2+9VfYxoTJi3+upU+OX9lyxYkWQxa4KELtCRWx76Kb06dMnyL7xjW8E2X333Zf5MZGG2Baw06dPD7Lrr78+yDp27Bh9zO9///tB1qtXryDr379/kH3wwQfRxwQaSexn7MUXXwyy7t27V2M4NcU7wwAAAMgtJsMAAADILSbDAAAAyC0mwwAAAMgtFtBF/OxnPwuy2BaFsUVkTfn73/8eZLVcLBcT2/r2c5/7XPS2sa13UXuxrWWHDx8eZLEFS01Zf/3wZaJ9+/YtG1gGsQV8DzzwQNmPgzTceeedQfb0008HWWxxsiQdcsghQbbffvsF2ZQpU4Ls29/+dpC98sor0eMgH2JbgUvxRZ1HH310kH322WdlH1OpPv300yCbOHFikJ12WrhPSuzvgEWLFpVnYDXAO8MAAADILSbDAAAAyC0mwwAAAMgtJsMAAADILRbQRcROAj/zzDOD7K677gqyjTfeOPqYl112WZDts88+QXbFFVcE2YwZM6KPWYp+/foF2RZbbBFkr732WvT+r776atnHhNKtXLkyyC666KIgu+aaa4Ksqe7GuhJbEBrLvvzlLwfZI488Ej3OsmXLgiz2fJBfsV0uf/rTn0Zvu9tuuwXZF77whSDr2bNnkI0ZMybIzjrrrOhxmtq5EY1lgw02iOaxXTJHjx4dZD/+8Y+DbP78+aUPrApiiwdZQAcAAAA0CCbDAAAAyC0mwwAAAMgtJsMAAADIrWYX0JnZtpJulbSVpFWSRrv7SDPrIOluSd0kzZL0HXf/oHJDra0HH3wwyGK7zMR2r5Piizm++93vBln//v2D7JRTTgmy2EISSZo1a1aQ9enTJ8hGjhwZZLHFSs8//3z0OCmguwWx7+vcuXMz3z+22C6rrbfeutX3zSt62zLPPvtsND/99NODbPz48Zkec6+99gqy2GJQiQV0a2rk7sYW+ErSRx99FGTHH398kO27775BFuuoJD355JNBtmLFiuaG2GIDBw4MssGDBwfZO++8E2QffJDUt69ZWd4ZXiHpPHf/kqR9JZ1pZj0lDZU02d17SJpc/ByoJ3QXKaK3SBXdRZKanQy7+1x3n1b8eImkmZK6SBogaWzxZmMlHVWpQQKtQXeRInqLVNFdpKpF1xk2s26S9pD0nKRO7j5XKvwAmFnHJu4zRNKQ0oYJlIbuIkX0Fqmiu0hJ5smwmbWTNF7SOe6+2Mwy3c/dR0saXXwMb80ggVLQXaSI3iJVdBepyTQZNrMNVCj27e4+oRjPN7POxX/ldZa0oFKDrFexnbSmTZsWve1JJ50UZLEdaTbffPMgu/3221sxunWLnYwf2/3u0ksvLfuxq4nu1tabb74ZZO+//370tt27dw+yzTbbLMg+/PDD0gdW5+htdmeccUY0v+6668p6nAMOOCCa33LLLWU9Tuoatbtz5syJ5rEFZ/fcc0+Q9ejRI8gmT54cfczYznTu4b8N7r///iAbMGBA9DFjOnToEGRt27YNsssvvzzIGu11uNlzhq3wT7qbJM1096vW+KNJkk4ofnyCpPC7AtQQ3UWK6C1SRXeRqizvDO8v6XhJfzOzl4rZMElXSrrHzE6W9Jakb1dmiECr0V2kiN4iVXQXSWp2MuzuT0lq6oSfQ8o7HKB86C5SRG+RKrqLVLEDHQAAAHKLyTAAAAByy2IrFCt2MC6V8m9iKzljK6Nj2z7vvvvumY8ze/bsILv++uuDbPjw4Zkfs8ZecPfe1Twg3S2fZ555Jpp/9atfDbLYds4t2Uq63rh7tmtMlUnKve3bt2+QnX/++UH2ta99LXr/cv/ddtZZZ0XzUaNGlfU4dYrX3Ca0adMmyA4//PAgGzo03HRv//33L+nYsUvWldr7G2+8McguuOCCIHv33XdLOk4VZeou7wwDAAAgt5gMAwAAILeYDAMAACC3mAwDAAAgt1hAhxSxmCNh5557bjS/6qqrguyoo44KstgWpKnI+wK6fv36RfMhQ4YEWWwRUmyr2NgiIin7QqLLLrssyKZNmxZkkyZNyvR4DYrX3BKtt1743uPee+8dvW1s0fx+++0XZPvuu2+QLVu2LMjGjRsXPc7IkSODLNb9VatWRe+fCBbQAQAAAOvCZBgAAAC5xWQYAAAAucVkGAAAALnFAjqkiMUcCdtnn32i+bPPPhtkf/rTn4LsoIMOKveQqiZPC+hOOeWUIGtql8vYbpwxixYtCrKnnnoqetu//vWvQTZhwoQgmz59epAlvmCoEnjNRapYQAcAAACsC5NhAAAA5BaTYQAAAOQWk2EAAADk1vrN3cDMtpV0q6StJK2SNNrdR5rZJZJOlfRu8abD3P2hSg0UaCm6W5+ee+65aN7UTmJ50yi9/fvf/x5ko0ePjt72wQcfzPSYCxYsCLI33nijZQNDxTRKd5E/zU6GJa2QdJ67TzOzTSS9YGaPFv9shLv/snLDA0pCd5EieotU0V0kqdnJsLvPlTS3+PESM5spqUulBwaUiu4iRfQWqaK7SFWLzhk2s26S9pC0+v85zzKz6WZ2s5lt3sR9hpjZVDObWtJIgRLQXaSI3iJVdBcpyTwZNrN2ksZLOsfdF0saJWlHSb1U+Jfgr2L3c/fR7t672hfsBlaju0gRvUWq6C5Sk2kybGYbqFDs2919giS5+3x3X+nuqyTdIGnvyg0TaB26ixTRW6SK7iJFWa4mYZJukjTT3a9aI+9cPD9IkgZKmlGZIQKtQ3eRokbp7TPPPJMpQ+NolO4if7JcTWJ/ScdL+puZvVTMhkkaZGa9JLmkWZJOq8gIgdaju0gRvUWq6C6SZO5evYOZVe9gaGQvVPucMrqLcnD3ql5Mmd6iTHjNRaoydZcd6AAAAJBbTIYBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQW0yGAQAAkFtMhgEAAJBbWTbdKKf3JL1Z/HjL4ueNgudTPdvV4Jiru1vPX5fW4PlUTy17K9X316Y1eD7VQ3fLi+dTPZm6W9VNN/7twGZTq30R70ri+eRDo31deD750WhfG55PfjTa14bnU384TQIAAAC5xWQYAAAAuVXLyfDoGh67Eng++dBoXxeeT3402teG55Mfjfa14fnUmZqdMwwAAADUGqdJAAAAILeYDAMAACC3qj4ZNrPDzexVM3vDzIZW+/ilMrObzWyBmc1YI+tgZo+a2evF3zev5Rhbwsy2NbMnzGymmb1sZj8o5sk+p0qhu/WF7maTem+lxuouvc0u9e42Um+lxu5uVSfDZtZG0nWS+knqKWmQmfWs5hjKYIykw9fKhkqa7O49JE0ufp6KFZLOc/cvSdpX0pnF70nKz6ns6G5dorvNaJDeSo3VXXqbQYN0d4wap7dSA3e32u8M7y3pDXf/p7svk3SXpAFVHkNJ3H2KpIVrxQMkjS1+PFbSUVUdVAncfa67Tyt+vETSTEldlPBzqhC6W2fobibJ91ZqrO7S28yS724j9VZq7O5WezLcRdLsNT5/u5ilrpO7z5UKZZHUscbjaRUz6yZpD0nPqUGeUxnR3TpGd5vUqL2VGuD7TG/XqVG72xDf50brbrUnwxbJuLZbHTCzdpLGSzrH3RfXejx1iO7WKbq7TvS2TtHbZtHdOtWI3a32ZPhtSduu8fk2kuZUeQyVMN/MOktS8fcFNR5Pi5jZBioU+3Z3n1CMk35OFUB36xDdbVaj9lZK+PtMbzNp1O4m/X1u1O5WezL8F0k9zGx7M2sr6RhJk6o8hkqYJOmE4scnSLq/hmNpETMzSTdJmunuV63xR8k+pwqhu3WG7mbSqL2VEv0+09vMGrW7yX6fG7q77l7VX5KOkPSapH9IuqDaxy/D+O+UNFfSchX+5XqypC1UWEH5evH3DrUeZwueTx8V/utpuqSXir+OSPk5VfBrRXfr6Bfdzfx1Srq3xefQMN2lty36WiXd3UbqbfH5NGx32Y4ZAAAAucUOdAAAAMgtJsMAAADILSbDAAAAyC0mwwAAAMgtJsMAAADILSbDAAAAyC0mwwAAAMit/wdPV5KE15UD1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plots(x_imgs[:8], titles=y_valid[:8])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Neural Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**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! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Functions, parameters, and training"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "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": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<img src=\"images/sgd2.gif\" alt=\"\" style=\"width: 70%\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "In the above gif from fast.ai's deep learning course, [intro to SGD notebook](https://github.com/fastai/courses/blob/master/deeplearning1/nbs/sgd-intro.ipynb)), 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": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### PyTorch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We will be using the open source [deep learning library, fastai](https://github.com/fastai/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](http://pytorch.org/), 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](http://pytorch.org/tutorials/beginner/blitz/tensor_tutorial.html):\n",
    "\n",
    "<img src=\"images/what_is_pytorch.png\" alt=\"pytorch\" style=\"width: 80%\"/>\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](https://www.youtube.com/watch?v=Z15cBAuY7Sc) 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](http://pytorch.org/tutorials/beginner/deep_learning_60min_blitz.html) or this [tutorial to learn by examples](http://pytorch.org/tutorials/beginner/pytorch_with_examples.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### About GPUs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "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",
    "- **Use crestle, through your browser**: [Crestle](https://www.crestle.com/) 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",
    "- **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](http://course.fast.ai/lessons/aws.html).]  AWS charges 90 cents per hour for this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Neural Net for Logistic Regression in PyTorch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "hidden": true
   },
   "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": "markdown",
   "metadata": {
    "hidden": true
   },
   "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": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "net = nn.Sequential(\n",
    "    nn.Linear(28*28, 100),\n",
    "    nn.LogSoftmax()\n",
    ").cuda()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "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": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "md = ImageClassifierData.from_arrays(path, (x,y), (x_valid, y_valid))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "loss=nn.NLLLoss()\n",
    "metrics=[accuracy]\n",
    "opt=optim.SGD(net.parameters(), lr=1e-2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Loss functions and metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "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:\n",
    "`-(y * log(p) + (1-y) * log (1-p))`\n",
    "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": 25,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def binary_loss(y, p): #also called Negative Log Loss\n",
    "    return np.mean(-(y * np.log(p) + (1-y)*np.log(1-p)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.164252033486018"
      ]
     },
     "execution_count": 26,
     "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": {
    "hidden": true
   },
   "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",
    "\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": {
    "hidden": true
   },
   "source": [
    "### Fitting the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "*Fitting* is the process by which the neural net learns the best parameters for the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2710165d0f92423ea74793e5f37da5cf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Epoch', max=1, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch      trn_loss   val_loss   accuracy                     \n",
      "    0      0.392233   0.334639   0.9091    \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[array([0.33464]), 0.9091]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit(net, md, n_epochs=1, crit=loss, opt=opt, metrics=metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "preds = predict(net, md.val_dl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000, 100)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**Question**: Why does our output have length 10 (for each image)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 8, 6, ..., 5, 6, 8])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds.argmax(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "preds = preds.argmax(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "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": 32,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9091"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(preds == y_valid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's see how some of our predictions look!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAF0CAYAAADGqzQSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xm81VW9//H3R4ZMMRENRBwwIQ1NIXFIseuQIuYDpKwk9SoOqOE1vTYgaupVf9okkhqGE5g4IaCYmilaOFw1RCMIpwpFOYCKCDgxfX5/nO3jEmsdzvecPa79fT0fDx6c82bv/V37nPfZLDbf9V3m7gIAAADyaKNqDwAAAACoFibDAAAAyC0mwwAAAMgtJsMAAADILSbDAAAAyC0mwwAAAMgtJsNVZmb9zewxM1toZp+Y2ZtmdreZ9ar22IANMbP9zeyPZrbYzJaZ2UwzO6na4wKyMLMjzGy6ma0o9HeGmR1c7XEBG2JmB5nZk2b2kZktMbPfmVmXao8rdUyGq6+TpOclnSnpMEnnSdpV0jNmtkM1BwY0xcx2l/SopHaSTpX0LUl/kXSTmZ1RzbEBzTGz0yTdp8bX3sGSvi1poqRNqjkuYEPM7ABJf5S0VI2vuT+Q9DVJ08zsM9UcW+qMTTdqj5ntLOklST90919VezzA+szs/0n6oaRO7r5infwZSe7uX63a4IANMLPukuZKOs/dr67uaIDszOxRSd0l7eLuqwvZXpKekzTc3X9TxeEljXeGa9O7hd9XVXUUQNPaq7GfH62XLxWvK6htJ0laK+n6ag8EaKF9JT3y6URYktz9L2qcMwyu2qjqAH9p1Qgza2Nm7c2sp6TfSloo6c4qDwtoyrjC7782s23MrKOZnSrpEEmjqjcsoFn91Pg/b8eY2T/MbLWZvWZmw6s9MKAZayStjOSfSNqtwmOpK5wmUSPMbIakPQufviZpoLvPreKQgA0q/PfcFEndCtEqSWe4+03VGxWwYWb2kqRt1DiBGCnpH2o8Z/h0SWe7++gqDg9okpk9p8bT0PZZJ9tB0r8krXJ3zhtuJSbDNcLMviTpc5K+oMZzMbtI6ufu86o5LiCm8D8Y09R47uU1ajxdYpCkMySd6O4Tqjg8oElm9oqknpK+5e6T18kfktRHUlfnL0bUIDM7VtJtki6X9Gs1LsAfK2k/NU6GP1vF4SWNyXANMrOOkuZJutPdT6/ycICAmU2U9BU1LuRYtU4+QVJ/SZ3dfW21xgc0xcz+V43nXn7O3Zevk58j6SpJ3dx9QbXGB2yImV2qxjfMNpbkku6StKmk3dz9C9UcW8o4Z7gGuftSNZ4q0aPaYwGa8GVJf113IlzwnKQtJXWu/JCATOY0kVvhd/4Rh5rl7hdK2krS7mr8X4whavyfjierOrDEMRmuQYULaO+ixnPZgFq0UFJvM2u/Xr6PpI8lLan8kIBMphR+779e3l/Sm+6+sMLjAVrE3T9w97+5+yIzO1yN8wWujlKEttUeQN6Z2RRJMyXNkrRM0hclnSNptSSuMYxada0aNym438x+o8ZzhgdKGiJplLvHVjwDteBBSY9L+q2ZbSXpn5KOVuOmR0OrOTBgQ8ysj6QBapwzSI1XRvmRpJ+7+9NVG1gd4JzhKjOzn0j6jqSd1Hjt1vmS/iTpChbPoZaZ2QBJP1Hjjokbq/F/MsZK+q27r6nm2IANMbPPSbpCjZPgLdR4qbUr3f32qg4M2AAz21WNl17dTdJnVFjA7O63VHVgdYDJMAAAAHKLc4YBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQW0VdWq1wfbvRktpIutHdr2zm9qzWQym84+6fL+YB6C6qwd2t+VttWEu6S29RIrzmIlWZutvqd4bNrI2k69R4zbtekoaYWa/WPh7QAq8Xc2e6i1TRXVQJr7lIVabuFnOaxN6SXnP3fxYusH+npEFFPB5QKXQXqaK7SBG9RU0rZjLcTY0bRHzqzUL2b8xsmJnNMLMZRRwLKCW6i1Q12116ixrEay5qWjHnDMfOfQvO8XH3sWrclYpzgFAr6C5S1Wx36S1qEK+5qGnFvDP8pqTt1vl8W0kLihsOUBF0F6miu0gRvUVNK2Yy/BdJPc1sRzNrL+kYSVNLMyygrOguUkV3kSJ6i5rW6tMk3H21mZ0p6WE1XirlZnefU7KRAWVCd5EquosU0VvUOnOv3Gk5nAOEEnne3ftW8oB0F6VQiusMtwS9RYnwmotUZeouO9ABAAAgt4ragQ4AgJR16NAhyE4++eQgGzQoflncgQMHBtmKFSuKHxiAiuGdYQAAAOQWk2EAAADkFpNhAAAA5BaTYQAAAOQWk2EAAADkFleTAADk1gknnBBko0aNynz/XXfdNcieffbZosYEoLJ4ZxgAAAC5xWQYAAAAucVkGAAAALnFZBgAAAC5xQK6Iuyxxx5Bds4550Rvu9NOOwXZJptsEmQjR44Mss033zzIHnrooehxli9fHs0BIO9OPPHEILv66quDbNWqVUH2y1/+MvqYM2fOLHpcAKqLd4YBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQW+burb+z2TxJyyWtkbTa3fs2c/vWH6zKOnToEGRvvPFGkHXs2LESw9Fbb70VzWML+O65555yD6fSnm+ua83JU3djYj0dPHhw9LZ9+vQJsn79+gVZ7GdkyZIlQbb11ltHj7Nw4cIgGzduXJDdcMMNQbZmzZroY9Yad7diH6Ml3a233rbEwIEDg2zKlClB9uGHHwbZT3/60yBrya50dYjXXKQqU3dLcTWJg9z9nRI8DlBpdBeportIEb1FTeI0CQAAAORWsZNhl/RHM3vezIbFbmBmw8xshpnNKPJYQCnRXaRqg92lt6hRvOaiZhV7msT+7r7AzDpLesTMXnL36evewN3HShorcQ4QagrdRao22F16ixrFay5qVlGTYXdfUPh9sZlNkbS3pOkbvleazMJ1L3PmzAmyd999N3r/F154IchiC5N22GGHINtuu+2CrFOnTtHj/PznPw+yJ554IsgWLVoUvX9e5Km72267bZDde++9QRbrY1OWLVsWZLGOt2vXLsiWLl0afcztt98+yK677roge++994Js+vTwW9fQ0BA9Tury1N0s2rdvH82/+93vBlnsdfzZZ58NspwvlisLeota1urTJMxsUzPb7NOPJR0maXapBgaUC91FquguUkRvUeuKeWe4i6QphX9pt5V0u7v/oSSjAsqL7iJVdBcporeoaa2eDLv7PyXtUcKxABVBd5EquosU0VvUOi6tBgAAgNwqxaYbubB8+fIgO+CAAypy7K222irIfvSjH0VvG8sPP/zwIBs/fnzxA0MSpk6dGmR77BG+SfPoo49G73/uuecG2TvvhNfNj+0g1xKf//zng+yhhx4KshtvvDHIRowYEWSxxXeoP029Fg4ZMiTIbrvttiA76aSTSj4mYH1du3YNsu9///vR28byVatWBVlsF9zLL788yGbOnBk9zvz586N5HvHOMAAAAHKLyTAAAAByi8kwAAAAcovJMAAAAHKLyTAAAAByi6tJJCC2cv+pp56K3ja2sjq2zS5Xk6hPsRXLvXv3DrK77747yI499tjoY65Zs6b4gWXw9ttvB9nLL78cZLE+N/XzgPrSt2/fILvwwgujt33llVeCbOjQoUFWqX4jP77whS8E2ZgxY4Ls0EMPLeo4nTt3DrIpU6YE2SeffBK9//777x9kTV15ot7xzjAAAAByi8kwAAAAcovJMAAAAHKLyTAAAAByiwV0Cdhiiy2CbOTIkZnvv80225RyOKhhscVyZhZkCxYsCLJqLyTad999gyy2pe7jjz8eZLHn/eKLL5ZmYKiKjTYK36uJbbvdvn376P3vv//+IKt2x1F/unXrFmSzZ88OsrZtw+nWqFGjoo95zTXXZDrOLrvsEmS/+MUvgqxjx47R48QWUsdeh2OL+OsN7wwDAAAgt5gMAwAAILeYDAMAACC3mAwDAAAgt5pdQGdmN0s6UtJid9+tkHWSdJek7pLmSfqOu79XvmHmxx577BFkEydODLIePXpE7x/bdencc88tfmAJymN3H3rooSBz9yD73ve+F2RXX3119DHfeOON4ge2js022yya33DDDUF27733Bllsp7zYjk8py2N319elS5cg++Y3v5n5/q+//noph4MM8tjbH//4x0HWpk2bIDvllFOC7NZbb818nHnz5gVZbOfNjTfeOMiaem2PvW5Onz49yGKL6pYtWxZ9zFRleWd4nKTD18tGSJrm7j0lTSt8DtSacaK7SNM40V2kZ5zoLRLU7GTY3adLWrJePEjS+MLH4yUdVeJxAUWju0gV3UWK6C1S1drrDHdx9wZJcvcGM+vc1A3NbJikYa08DlBqdBepytRdeosaw2sual7ZN91w97GSxkqSmYUnLwI1iu4iRfQWqaK7qJbWToYXmVnXwr/yukpaXMpB5cUJJ5wQZP/zP/8TZNttt12QffTRR9HHPOOMM4Js/vz5rRhd3cpdd2OduvDCC4PsD3/4Q/T+/fv3D7JiOjVp0qRo/sUvfjHIYjvQxbo/Z86cVo8nIbnq7uGHr3/qadyjjz4azceMGVPK4aD16qK3n/vc56L5McccE2SxneVasliuGNddd12QnXXWWdHb9uzZM8hiu9pddNFFQVZvC/Nbe2m1qZI+ncmdIOm+0gwHKDu6i1TRXaSI3qLmNTsZNrM7JP2vpJ3N7E0zO1nSlZIONbNXJR1a+ByoKXQXqaK7SBG9RaqaPU3C3cP/p2x0SInHApQU3UWq6C5SRG+RKnagAwAAQG4xGQYAAEBuWWyr1rIdLAeXSunQoUM0/+EPfxhkF1xwQZBttFH475MlS9a/hrnUr1+/6HFeeuml5oZYD553976VPGDK3Y1tzzl+/PggO/roo6P3f+2114LswAMPDLKGhoYg+81vfhNkw4bFLyP6ox/9KMhiq7JT5u5WyeOl0tu2bcMz9ubOnRtkO+ywQ5DtuOOO0cd86623ih8YPpX719y99947mj/zzDNBduihhwbZtGnTSj6mrAYPHhzNJ0+eHGSxOeHSpUuDLHYlinfffbcVoyu7TN3lnWEAAADkFpNhAAAA5BaTYQAAAOQWk2EAAADkVmu3Y0YTxo0bF82/+c1vZrr/PffcE2RXX311kOVkoRxK4OOPPw6yU045Jcg6d+4cvf9//Md/BNmf//znIJs4cWKQHXfccUHW1HbM9bZYDtnFFm/utNNOQRbbbr7aC+Vi20YPHDgwyGLbnf/xj38MstjPK6qvT58+mW/7wgsvlHEkLffggw9G89ji6NjPXayTH3zwQfEDqyG8MwwAAIDcYjIMAACA3GIyDAAAgNxiMgwAAIDcYgFdicVOPm+JMWPGBNnTTz9d1GMC61u+fHmQDRo0KHrbiy++OMjOPvvsIBsxYkSmY19zzTWZbof82H777TPdrn379mUeSdNOPPHEaB7bZTG26+Ppp58eZLGdve69997ocU466aRmRohyevLJJ6P52rVrg+yRRx4JsiOPPDLIYrt2lsPOO+8czWM97d+/f5BtsskmQVZvCz15ZxgAAAC5xWQYAAAAucVkGAAAALnFZBgAAAC51ewCOjO7WdKRkha7+26F7GJJp0p6u3Czke4e3+IkZ2I7CknSHnvs0er7xxbVXXnlldH7L1iwINNx8oDutsyyZcui+U9/+tMgO/TQQ4OsV69emY7z9a9/PZo3tUAlj/LW3R49emS6XaV23uzYsWOQXXXVVdHbxhYhrV69Oshii6r69esXZLFdG6U0FtDVc2/nzJkTzX//+98HWWwx8ty5c4MstiuhFN+l87HHHguybt26BVlssVxsF1tJ6tq1a5DFunvfffdF719PsrwzPE5SuN+kNMrdexd+JVds5MI40V2kaZzoLtIzTvQWCWp2Muzu0yUtqcBYgJKiu0gV3UWK6C1SVcw5w2ea2Swzu9nMtmjqRmY2zMxmmNmMIo4FlBLdRaqa7S69RQ3iNRc1rbWT4TGSdpLUW1KDpF81dUN3H+vufd29byuPBZQS3UWqMnWX3qLG8JqLmteqHejcfdGnH5vZDZLCM8hzKrZblxRfILLnnnsGWWwnpjPPPDPIjj766Ohxhg4dGmQPP/xw9LZ5RHdb7oADDgiynj17tvrxfvKTn0Tz119/PchuueWWVh+n3tRzd2MLgRYuXBhksUVo5RDbbS62qE6SbrvttiAbPXp0kL3xxhtBFltA9eUvfznDCNNRz72VpCFDhgTZFVdcEWRnnXVWkH3nO9+JPmYsX7IkPPukU6dOWYbYIuV4zBS06p1hM1t3CeJgSbNLMxygvOguUkV3kSJ6ixRkubTaHZIOlLSVmb0p6SJJB5pZb0kuaZ6k08o4RqBV6C5SRXeRInqLVDU7GXb38P8ApJvKMBagpOguUkV3kSJ6i1SxAx0AAAByi8kwAAAAcqtVV5NA0z766KNofuyxxwZZ27bhl7+pLXHXt/XWW0fzKVOmBNl///d/B9n111+f6TjAQQcdFGTuHmSDBw8OstgK6Nj2pVJ82/F33nknyO6///7o/ZGuffbZJ8hWrlxZhZG03IIFC4Js2223DbKxY8cG2Ve+8pUg4+o/aYn9nX/22WcH2d133x1ksXlBU7p06ZLpdqtWrQqy2M+XJO24445B9uGHH2YeUz3hnWEAAADkFpNhAAAA5BaTYQAAAOQWk2EAAADkFgvoKuTjjz/OdLvevXsH2ahRo4IstqhJkjbeeOMgGzFiRJCxgA7r23333aP5D37wgyCLLXabOnVqpuOceuqp0fymm8LLkd5+++1BtuuuuwZZbKtbpGPy5MlBduSRR1bk2GaWKWvKj3/840y3iy06vfbaa4Ns5MiRmY+NdDz99NOZsnL43e9+F827d+8eZGvWrCnzaGoT7wwDAAAgt5gMAwAAILeYDAMAACC3mAwDAAAgt1hAl9Emm2wSZOXYqWXWrFlBdvTRRwfZzTffHL3/oEGDgmz77bcPsq5duwZZQ0NDliGiTm222WbRPLZT4j333NPq40ycODGa77DDDkH2s5/9LMj23HPPIGMBXf3p2LFjkMUWAt12223R+8d6e8wxxwRZp06dgmzAgAFZhihJ+uCDD4LsySefDLKf//znQfb4449nPg5QCTvttFO1h1AVvDMMAACA3GIyDAAAgNxiMgwAAIDcYjIMAACA3Gp2AZ2ZbSfpVklbS1oraay7jzazTpLuktRd0jxJ33H398o31MqJnUAeWxDxwAMPBNns2bOjjxlbnHbyyScHWbt27YKsW7duQdajR4/ocWL+8Y9/ZBpPvcljd4sR2/1QkhYuXBhksZ+HYsV244rtVjd8+PAgmzJlSsnHUy157O0LL7wQZKecckqQHXvssZmyYi1btizImlr4edlllwXZ66+/XvIxpSCP3U3BihUrqj2EmpflneHVks519y9J2lfScDPrJWmEpGnu3lPStMLnQC2hu0gRvUWq6C6S1Oxk2N0b3H1m4ePlkuZK6iZpkKTxhZuNl3RUuQYJtAbdRYroLVJFd5GqFl1n2My6S+oj6VlJXdy9QWr8ATCzzk3cZ5ikYcUNEygO3UWK6C1SRXeRksyTYTPrIGmSpLPdfZmZZbqfu4+VNLbwGN6aQQLFoLtIEb1FquguUpNpMmxm7dRY7AnuPrkQLzKzroV/5XWVtLhcg6y0b3/720G29dZbB9lJJ51U8mPHXjTcs78mxE6UP/3004saU8ry1t1ixHYqlKTnnnuuIsdfuXJlkL33XrjG5oADDgiy2C5iS5YsKc3AqiBvvb399tuDLLbz5quvvhpkbdq0iT5mU/n6JkyYEGTz5s0LsthCZITy1t0UTJ8+PZqfdtppQda5c/RN+7rX7DnD1jg7u0nSXHe/ap0/mirphMLHJ0i6r/TDA1qP7iJF9BaportIVZZ3hveXdLykv5nZi4VspKQrJd1tZidLekNS+HYqUF10Fymit0gV3UWSmp0Mu/uTkpo64eeQ0g4HKB26ixTRW6SK7iJV7EAHAACA3GIyDAAAgNxq0XWG82LLLbes9hD+zaRJk4Ls0ksvjd528eJwkW5sO11gfU1dtaRfv35BdswxxwTZY489FmQdOnQIsvbt20ePs8suuwTZXnvtFWTXXXddkKV85QhI77//fpAdcgj/qw6UwkYbxd/3jF29KjaHyAPeGQYAAEBuMRkGAABAbjEZBgAAQG4xGQYAAEBusYAuYuTIkUH26KOPBtlxxx0XZNtss030MWMLRGKuueaaIHviiSeCbPXq1ZkeD8hq7ty50Ty21XFs+9x33303yFqygC62mOOpp54Ksosvvjh6fwBAaO3atdG8qUXTecQ7wwAAAMgtJsMAAADILSbDAAAAyC0mwwAAAMgtFtBFrFq1KsgefvjhTBmQqj/84Q/R/Nprrw2y2K50vXv3Lur4559/fpDdfPPNQcZucwBQHocddliQjRkzpgojqSzeGQYAAEBuMRkGAABAbjEZBgAAQG4xGQYAAEBuNbuAzsy2k3SrpK0lrZU01t1Hm9nFkk6V9HbhpiPd/cFyDRRoKbrbMosWLYrmP/jBDyo8knyjt0gV3a1NK1asyHzbtm3zeV2FLM96taRz3X2mmW0m6Xkze6TwZ6Pc/ZflGx5QFLqLFNFbpIruIknNTobdvUFSQ+Hj5WY2V1K3cg8MKBbdRYroLVJFd5GqFp0zbGbdJfWR9GwhOtPMZpnZzWa2RRP3GWZmM8xsRlEjBYpAd5EieotU0V2kJPNk2Mw6SJok6Wx3XyZpjKSdJPVW478EfxW7n7uPdfe+7t63BOMFWozuIkX0Fqmiu0hNpsmwmbVTY7EnuPtkSXL3Re6+xt3XSrpB0t7lGybQOnQXKaK3SBXdRYqyXE3CJN0kaa67X7VO3rVwfpAkDZY0uzxDBFqH7iJF9Baporu16Yknnsh824MPPriMI6ldWa4msb+k4yX9zcxeLGQjJQ0xs96SXNI8SaeVZYRA69FdpIjeIlV0F0kyd6/cwcwqdzDUs+crfU4Z3UUpuLtV8nj0FiXCa27COnbsGM2XLFkSZB999FGQbbrppiUfUwVl6i470AEAACC3mAwDAAAgt/K57x4AAEAOLF26NJpvtBHvh36KrwQAAAByi8kwAAAAcovJMAAAAHKLyTAAAAByq9IL6N6R9Hrh460Kn9cLnk/l7FCFY37a3Vr+urQGz6dyqtlbqba/Nq3B86kcultaPJ/KydTdim668W8HNptR6Yt4lxPPJx/q7evC88mPevva8Hzyo96+Njyf2sNpEgAAAMgtJsMAAADIrWpOhsdW8djlwPPJh3r7uvB88qPevjY8n/yot68Nz6fGVO2cYQAAAKDaOE0CAAAAucVkGAAAALlV8cmwmR1uZi+b2WtmNqLSxy+Wmd1sZovNbPY6WScze8TMXi38vkU1x9gSZradmT1uZnPNbI6Z/aCQJ/ucyoXu1ha6m03qvZXqq7v0NrvUu1tPvZXqu7sVnQybWRtJ10kaIKmXpCFm1quSYyiBcZIOXy8bIWmau/eUNK3weSpWSzrX3b8kaV9Jwwvfk5SfU8nR3ZpEd5tRJ72V6qu79DaDOunuONVPb6U67m6l3xneW9Jr7v5Pd18p6U5Jgyo8hqK4+3RJS9aLB0kaX/h4vKSjKjqoIrh7g7vPLHy8XNJcSd2U8HMqE7pbY+huJsn3Vqqv7tLbzJLvbj31Vqrv7lZ6MtxN0vx1Pn+zkKWui7s3SI1lkdS5yuNpFTPrLqmPpGdVJ8+phOhuDaO7TarX3kp18H2mtxtUr92ti+9zvXW30pNhi2Rc260GmFkHSZMkne3uy6o9nhpEd2sU3d0geluj6G2z6G6NqsfuVnoy/Kak7db5fFtJCyo8hnJYZGZdJanw++Iqj6dFzKydGos9wd0nF+Kkn1MZ0N0aRHebVa+9lRL+PtPbTOq1u0l/n+u1u5WeDP9FUk8z29HM2ks6RtLUCo+hHKZKOqHw8QmS7qviWFrEzEzSTZLmuvtV6/xRss+pTOhujaG7mdRrb6VEv8/0NrN67W6y3+e67q67V/SXpCMkvSLpH5LOr/TxSzD+OyQ1SFqlxn+5nixpSzWuoHy18HunVn5dpktaIWmZpBmSDq7A8+mnxv96miXpxcKvI0rxnOrtF91t9vH/UOjSZRV6PnQ329cp6d4WnkPJuivpwEJv1v+1tELPhd5m/1ol3d1Sv+bS3fL9YjvmGmBmp0m6tvDrQTW+Y99b0hx3/301xwZkYWZDJF0laWtJl7v7BVUeEhBlZgdKelzSWWp89/FTq919RlUGBWRAd8unbbUHkHeFFZlXS/qRu1+9zh89XJUBAS1kZh0ljZJ0jqTbqzwcIKu57v5MtQcBtALdLTG2Y66+kyStlXR9tQcCtNLP1fi/GHdUeyAAALQUk+Hq6yfpJUnHmNk/zGx1YevJ4dUeGNAcM+sn6T8lfb/aYwFaaIKZrTGzd83sdjPbvtoDAjKiuyXGaRLVt03h1y8kjVTjQoFvS7rWzNq6++hqDg5oSuESO7+V9Et3f7na4wEyel/SryT9WY2Llfuo8bX3f82sj7snd1ko5AbdLRMW0FWZmb0iqaekb/n/XbNPZvaQGove1fkmoQaZ2QVqPM1nV3f/qJC5WECHxJjZVyQ9J+lKuouU0N3S4DSJ6nu38Psj6+V/lNRFUtfKDgdoXuG/5c6XdKGkz5hZx8JCOq3zeZvqjRDIzt1nqvESXntVeyxAS9Dd0mAyXH1zmsg/3YpybaUGArTAFyRtLOk2Se+t80uSflj4+MvVGRrQKia2+0Wa6G6RmAxX35TC7/3Xy/tLetPdF1Z4PEAWL0o6KPJLapwgHyTpteoMDWgZM+sr6YuSnq32WICWoLulwQK66ntQjRfR/q2ZbSXpn5KOlnSYpKHVHBjQFHdfKulP6+eNu3XqdXcP/gyoBWY2QdK/JM2UtFSNazPOk/SWpGuqODRgg+hu+TAZrjJ3dzM7StIVki6RtIUaL7V2rLuzgQEAlNZsSUMk/ZekTSQtlDRZ0kXu/k41BwY0g+6WCVeTAAAAQG5xzjAAAAByi8kwAAAAcovJMAAAAHLKY/CvAAAZeElEQVSLyTAAAAByq6irSZjZ4ZJGS2oj6UZ3v7KZ27NaD6Xwjrt/vpgHoLuoBne35m+1YS3pLr1FifCai1Rl6m6r3xkubLV6naQBknpJGmJmvVr7eEALvF7MnekuUkV3USW85iJVmbpbzGkSe0t6zd3/6e4rJd0paVARjwdUCt1FquguUkRvUdOKmQx3kzR/nc/fLGT/xsyGmdkMM5tRxLGAUqK7SFWz3aW3qEG85qKmFXPOcOzct+AcH3cfK2msxDlAqBl0F6lqtrv0FjWI11zUtGLeGX5T0nbrfL6tpAXFDQeoCLqLVNFdpIjeoqYVMxn+i6SeZrajmbWXdIykqaUZFlBWdBeportIEb1FTWv1aRLuvtrMzpT0sBovlXKzu88p2ciAMqG7SBXdRYroLWqduVfutBzOAUKJPO/ufSt5QLqLUijFdYZbgt6iRHjNRaoydZcd6AAAAJBbTIYBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQW0yGAQAAkFtMhgEAAJBbrd50Ay0zYMCAIDvnnHOC7NBDDw2y2LWgX3311ehx7r777iAbM2ZMkC1YwE6YAAAAvDMMAACA3GIyDAAAgNxiMgwAAIDcYjIMAACA3LLY4qyyHcyscgerkjPOOCOajxo1Ksjat29f7uFIkh5//PEgO+6444KsoaGhEsMphefdvW8lD5iH7qL83N0qeTx6ixLhNbcFxo8fH2THH398kD3wwAPR+0+aNCnInn766SCbP39+pvGsXLkymq9ZsybT/ROXqbu8MwwAAIDcYjIMAACA3GIyDAAAgNxiMgwAAIDcKmoHOjObJ2m5pDWSVlf6BPtq+8Y3vhFkv/zlL6O3jS2We+GFF4JsxIgRQTZnzpzMYzr55JOD7JJLLgmy8847L8jOOuuszMdJXd67u+mmmwbZyJEjo7e94IILgiy28PbSSy8Nsj322CPIBg4cmGWIaELeu4s05am3L730UpCtXbs2yGJziA3lrXXLLbdE89NOOy3IVq9eXdJjp6IU2zEf5O7vlOBxgEqju0gV3UWK6C1qEqdJAAAAILeKnQy7pD+a2fNmNix2AzMbZmYzzGxGkccCSonuIlUb7C69RY3iNRc1q9jTJPZ39wVm1lnSI2b2krtPX/cG7j5W0lgp7Ytoo+7QXaRqg92lt6hRvOaiZhU1GXb3BYXfF5vZFEl7S5q+4Xul6cgjjwyyO+64I8g++9nPRu9/7733Bllst7pFixa1YnT/57LLLguy2M5yhx12WFHHSV2euhuz5ZZbBllTC+hiuyHNnDkz03G+9rWvBVmXLl2ity22+3mR9+7Wk9jPQo8ePaK33XjjjYNsyJAhQTZhwoQga2oHsqeeeqq5IZZMnnp7xRVXBNnf/va3IOvfv3/mx9xrr72CbPvttw+y2Bxk6NCh0ce87bbbgiy2Y20etPo0CTPb1Mw2+/RjSYdJml2qgQHlQneRKrqLFNFb1Lpi3hnuImmKmX36OLe7+x9KMiqgvOguUkV3kSJ6i5rW6smwu/9TUngRUaDG0V2kiu4iRfQWtY5LqwEAACC3SrHpRt1p2zb8ssR2cYvt4jVr1qzoY8Z2enn77bdbMboNi+0MduONNwbZlClTSn5spKN79+4lf8xVq1YF2eabbx5kvXr1it6fBXSoF7vttluQffe73w2yk046Kci6du0afczYa3tMU4ulYtq0aZP5tijO73//+0xZsQYMGBBkDzzwQPS2RxxxRJCxgA4AAADIGSbDAAAAyC0mwwAAAMgtJsMAAADILSbDAAAAyC2uJhFx6qmnBlmfPn2C7JNPPgmyE088MfqY5bhyRDHefffdag8BVfTVr3615I953333BVnsKix9+/aN3j+vq5iRht69e0fzc845J8i+/vWvB9nWW29d8jHFLF++PMgee+yxihwbldWpU6cgu+iii4Js9erV0fs3dZWJPOKdYQAAAOQWk2EAAADkFpNhAAAA5BaTYQAAAOQWC+gi/uu//ivT7U4//fQge/HFF0s9HKAosS1Xv/WtbwXZ2rVro/dvavEF0FKxre4laeONNw6yFStWlHs4kuILOm+55ZYg22mnnaL3/8xnPlPyMcX8/e9/D7ILLrggyGKLo5988smyjAnF2WyzzaJ5v379gqx9+/ZBdv755wdZrM+33npr9Dh/+tOfmhlhfvDOMAAAAHKLyTAAAAByi8kwAAAAcovJMAAAAHKr2QV0ZnazpCMlLXb33QpZJ0l3SeouaZ6k77j7e+UbZm168803qz0EbADdbdSlS5cg22uvvYLsX//6V/T+s2bNynScVatWBdmaNWuCrEePHpkeL8/qtbux3bEk6aijjgqySZMmBdnFF1+c+Vi77757kP3kJz8Jsthi0nbt2gWZmUWP4+6Zx5RF7HlL0n/+538G2UcffVTSYxerXnvbUh06dAiyK664Ishi3ZOK263w2WefDbIrr7yy1Y+XF1neGR4n6fD1shGSprl7T0nTCp8DtWac6C7SNE50F+kZJ3qLBDU7GXb36ZKWrBcPkjS+8PF4SeE/64Eqo7tIFd1FiugtUtXa6wx3cfcGSXL3BjPr3NQNzWyYpGGtPA5QanQXqcrUXXqLGsNrLmpe2TfdcPexksZKkpmV9uQqoIzoLlJEb5Equotqae1keJGZdS38K6+rpMWlHFSlxBZYSFLPnj2DbPny5UH28ssvl3xMKLu66G45vPrqq0Xd/7XXXguy+fPnB1nv3r2LOk6OJdXdz33uc0F2/PHHR2+7/fbbB9muu+4aZLGFSTvvvHP0Mb/xjW80N8QWaWoBXUxsF7jf/e53QTZ58uQgq8Pd4pLqbSnsv//+QTZ8+PCKHDv2M9LU7qL4P629tNpUSScUPj5B0n2lGQ5QdnQXqaK7SBG9Rc1rdjJsZndI+l9JO5vZm2Z2sqQrJR1qZq9KOrTwOVBT6C5SRXeRInqLVDV7moS7D2nijw4p8ViAkqK7SBXdRYroLVLFDnQAAADILSbDAAAAyK2yX1qtlrVtG3/6bdq0CbIPP/wwyNiOGSk4+OCDM91u1KhRRR0n9vMU+1nq2rVr9P6xqw8sW7asqDGhejp16hRkm266afS2Wbc0Puecc4KsHNsk/+Uvfwmyu+66K3rbBx98MMhWrFgRZG+99Varx4O09OvXr6j7L14cXnBjzJgxQbbRRuH7mRdeeGGQxbaClqRTTjklyN57r653ym4S7wwDAAAgt5gMAwAAILeYDAMAACC3mAwDAAAgt3K9gK7attxyyyA78sgjg+zcc8/N/Jjz5s0Lsu7duwfZwoULg+yee+4JsltuuSV6nFWrVmUeE6prv/32C7JFixYF2RNPPFHUcWKLTB944IEgO/3006P333zzzYOMBXTpir0Wvf3229HbxhbbVcqll14aZL/+9a+DbMmSJZUYDurAJZdcEmTPP/98kH3wwQfR+//5z38OspUrVwZZbPHoxIkTg2zatGnR49x4441BdvLJJwfZ0qVLo/evJ7wzDAAAgNxiMgwAAIDcYjIMAACA3GIyDAAAgNxiAV1GsQUeffv2DbIZM2ZE79+jR48ge/TRR4Ns++23D7KPPvooyP76179GjxNbtBLLhg4dGmRf//rXg6x///7R43zrW9+K5qiu2A5fRxxxRJDFFmM0tZijGHlYeIHsmlrIs/POO7f6MadPnx7NJ02aFGS33357kMV23Fq7dm2rxwOsXr06yO69996SHye2y+Ls2bOD7NRTT43ef8qUKUH2+OOPB9m1117bitGlhXeGAQAAkFtMhgEAAJBbTIYBAACQW0yGAQAAkFvNLqAzs5slHSlpsbvvVsgulnSqpE+3Exrp7g+Wa5Dl0tSOQu+//36QxXbHimVf+MIXoo/52GOPBdm2224bZLEFJsOHDw+yV155JXqcrKZOnRpksZPpd9lll6KOU0313N2mbLLJJkG2ww47BNn8+fMrMZzoz1JTYj9PlRpnranX7p533nnRPLbzZmwxccyBBx5YzJBQQvXa29TF/r6XpDvvvDPIYj+jd911V5A1tZtkqrK8MzxO0uGRfJS79y78otioReNEd5GmcaK7SM840VskqNnJsLtPl8Sm7EgO3UWq6C5SRG+RqmLOGT7TzGaZ2c1mtkVTNzKzYWY2w8ziF+AFKo/uIlXNdpfeogbxmoua1trJ8BhJO0nqLalB0q+auqG7j3X3vu4e7lABVB7dRaoydZfeosbwmoua16od6Nx90acfm9kNkn5fshFVUGxnNklqaGgIstjinu9973tB1qtXr+hjxhbLxXagGzx4cJCVY2ew2LFvvPHGIDvssMNKfuxqqpfuFqt9+/ZBtueee0Zv+/HHHwdZbPHpZz/72SCL7ZDUlDFjxgTZwQcfHGSrVq3K/Jj1pB66u2LFimgeW8hz3HHHBVm3bt2CbOHChdHHnDhxYpBddNFFQdbUQmqURj30tl6NHj06yIYMGRJkw4YNC7LLL7+8LGOqlla9M2xmXdf5dLCkcP8/oAbRXaSK7iJF9BYpyHJptTskHShpKzN7U9JFkg40s96SXNI8SaeVcYxAq9BdpIruIkX0FqlqdjLs7uF75tJNZRgLUFJ0F6miu0gRvUWq2IEOAAAAucVkGAAAALnVqqtJ1Lv77rsvyGLbEg8dOjTzY8au1HD22WcH2Ycffpj5MUutS5cuQda1a9fILeNbpb7xxhslHxPKI/a9njEjfmnP1atXB1nsqgCxK1TEtoduSr9+/YLsG9/4RpDde++9mR8TaYhtATtr1qwgu/7664Osc+fO0cf8/ve/H2S9e/cOsoEDBwbZe++9F31MoJ7EfsZeeOGFIOvRo0clhlNVvDMMAACA3GIyDAAAgNxiMgwAAIDcYjIMAACA3GIBXcTPfvazIIttURhbRNaUv//970FWzcVyMbGtbz/zmc9EbxvbehfVF9ta9oorrgiy2IKlprRtG75MdOzYsWUDyyC2gO/+++8v+XGQhjvuuCPInnrqqSCLLU6WpEMOOSTI9ttvvyCbPn16kH37298Ospdeeil6HORDbCtwKb6o8+ijjw6yTz75pORjKtbHH38cZFOmTAmy004L90mJ/R2wdOnS0gysCnhnGAAAALnFZBgAAAC5xWQYAAAAucVkGAAAALnFArqI2Engw4cPD7I777wzyDbddNPoY1566aVBts8++wTZ5ZdfHmSzZ8+OPmYxBgwYEGRbbrllkL3yyivR+7/88sslHxOKt2bNmiC78MILg+yaa64Jsqa6G+tKbEFoLPvyl78cZA8//HD0OCtXrgyy2PNBfsV2ufzpT38ave3uu+8eZJ///OeDrFevXkE2bty4IDvzzDOjx2lq50bUl3bt2kXz2C6ZY8eODbIf//jHQbZo0aLiB1YBscWDLKADAAAA6gSTYQAAAOQWk2EAAADkFpNhAAAA5FazC+jMbDtJt0raWtJaSWPdfbSZdZJ0l6TukuZJ+o67v1e+oVbXAw88EGSxXWZiu9dJ8cUc3/3ud4Ns4MCBQXbKKacEWWwhiSTNmzcvyPr16xdko0ePDrLYYqXnnnsuepwU0N1Gse9rQ0ND5vvHFttltc0227T6vnlFb1vmmWeeieann356kE2aNCnTY+61115BFlsMKrGAbl313N3YAl9J+uCDD4Ls+OOPD7J99903yGIdlaQnnngiyFavXt3cEFts8ODBQTZ06NAge+utt4LsvfeS+vY1K8s7w6slnevuX5K0r6ThZtZL0ghJ09y9p6Rphc+BWkJ3kSJ6i1TRXSSp2cmwuze4+8zCx8slzZXUTdIgSeMLNxsv6ahyDRJoDbqLFNFbpIruIlUtus6wmXWX1EfSs5K6uHuD1PgDYGadm7jPMEnDihsmUBy6ixTRW6SK7iIlmSfDZtZB0iRJZ7v7MjPLdD93HytpbOExvDWDBIpBd5EieotU0V2kJtNk2MzaqbHYE9x9ciFeZGZdC//K6yppcbkGWatiO2nNnDkzetuTTjopyGI70myxxRZBNmHChFaMbsNiJ+PHdr+75JJLSn7sSqK71fX6668H2bvvvhu9bY8ePYJs8803D7L333+/+IHVOHqb3RlnnBHNr7vuupIe54ADDojmt9xyS0mPk7p67e6CBQuieWzB2d133x1kPXv2DLJp06ZFHzO2M517+G+D++67L8gGDRoUfcyYTp06BVn79u2D7LLLLguyensdbvacYWv8J91Nkua6+1Xr/NFUSScUPj5BUvhdAaqI7iJF9BaportIVZZ3hveXdLykv5nZi4VspKQrJd1tZidLekPSt8szRKDV6C5SRG+RKrqLJDU7GXb3JyU1dcLPIaUdDlA6dBcpordIFd1FqtiBDgAAALnFZBgAAAC5ZbEVimU7GJdK+TexlZyxldGxbZ/32GOPzMeZP39+kF1//fVBdsUVV2R+zCp73t37VvKAdLd0nn766Wj+1a9+Nchi2zm3ZCvpWuPu2a4xVSIp97Z///5Bdt555wXZ1772tej9S/1325lnnhnNx4wZU9Lj1Chec5vQpk2bIDv88MODbMSIcNO9/fffv6hjxy5ZV2zvb7zxxiA7//zzg+ztt98u6jgVlKm7vDMMAACA3GIyDAAAgNxiMgwAAIDcYjIMAACA3GIBHVLEYo6EnXPOOdH8qquuCrKjjjoqyGJbkKYi7wvoBgwYEM2HDRsWZLFFSLGtYmOLiKTsC4kuvfTSIJs5c2aQTZ06NdPj1Slec4u00Ubhe49777139LaxRfP77bdfkO27775BtnLlyiCbOHFi9DijR48Oslj3165dG71/IlhABwAAAGwIk2EAAADkFpNhAAAA5BaTYQAAAOQWC+iQIhZzJGyfffaJ5s8880yQ/elPfwqygw46qNRDqpg8LaA75ZRTgqypXS5ju3HGLF26NMiefPLJ6G3/+te/BtnkyZODbNasWUGW+IKhcuA1F6liAR0AAACwIUyGAQAAkFtMhgEAAJBbTIYBAACQW22bu4GZbSfpVklbS1oraay7jzaziyWdKuntwk1HuvuD5Roo0FJ0tzY9++yz0bypncTypl56+/e//z3Ixo4dG73tAw88kOkxFy9eHGSvvfZaywaGsqmX7iJ/mp0MS1ot6Vx3n2lmm0l63sweKfzZKHf/ZfmGBxSF7iJF9BaportIUrOTYXdvkNRQ+Hi5mc2V1K3cAwOKRXeRInqLVNFdpKpF5wybWXdJfSR9+v+cZ5rZLDO72cy2aOI+w8xshpnNKGqkQBHoLlJEb5EquouUZJ4Mm1kHSZMkne3uyySNkbSTpN5q/Jfgr2L3c/ex7t630hfsBj5Fd5EieotU0V2kJtNk2MzaqbHYE9x9siS5+yJ3X+PuayXdIGnv8g0TaB26ixTRW6SK7iJFWa4mYZJukjTX3a9aJ+9aOD9IkgZLml2eIQKtQ3eRonrp7dNPP50pQ/2ol+4if7JcTWJ/ScdL+puZvVjIRkoaYma9JbmkeZJOK8sIgdaju0gRvUWq6C6SZO5euYOZVe5gqGfPV/qcMrqLUnD3il5Mmd6iRHjNRaoydZcd6AAAAJBbTIYBAACQW0yGAQAAkFtMhgEAAJBbTIYBAACQW0yGAQAAkFtMhgEAAJBbWTbdKKV3JL1e+Hirwuf1gudTOTtU4ZifdreWvy6twfOpnGr2Vqrtr01r8Hwqh+6WFs+ncjJ1t6Kbbvzbgc1mVPoi3uXE88mHevu68Hzyo96+Njyf/Ki3rw3Pp/ZwmgQAAAByi8kwAAAAcquak+GxVTx2OfB88qHevi48n/yot68Nzyc/6u1rw/OpMVU7ZxgAAACoNk6TAAAAQG4xGQYAAEBuVXwybGaHm9nLZvaamY2o9PGLZWY3m9liM5u9TtbJzB4xs1cLv29RzTG2hJltZ2aPm9lcM5tjZj8o5Mk+p3Khu7WF7maTem+l+uouvc0u9e7WU2+l+u5uRSfDZtZG0nWSBkjqJWmImfWq5BhKYJykw9fLRkia5u49JU0rfJ6K1ZLOdfcvSdpX0vDC9yTl51RydLcm0d1m1ElvpfrqLr3NoE66O07101upjrtb6XeG95b0mrv/091XSrpT0qAKj6Eo7j5d0pL14kGSxhc+Hi/pqIoOqgju3uDuMwsfL5c0V1I3JfycyoTu1hi6m0nyvZXqq7v0NrPku1tPvZXqu7uVngx3kzR/nc/fLGSp6+LuDVJjWSR1rvJ4WsXMukvqI+lZ1clzKiG6W8PobpPqtbdSHXyf6e0G1Wt36+L7XG/drfRk2CIZ13arAWbWQdIkSWe7+7Jqj6cG0d0aRXc3iN7WKHrbLLpbo+qxu5WeDL8pabt1Pt9W0oIKj6EcFplZV0kq/L64yuNpETNrp8ZiT3D3yYU46edUBnS3BtHdZtVrb6WEv8/0NpN67W7S3+d67W6lJ8N/kdTTzHY0s/aSjpE0tcJjKIepkk4ofHyCpPuqOJYWMTOTdJOkue5+1Tp/lOxzKhO6W2Pobib12lsp0e8zvc2sXrub7Pe5rrvr7hX9JekISa9I+oek8yt9/BKM/w5JDZJWqfFfridL2lKNKyhfLfzeqdrjbMHz6afG/3qaJenFwq8jUn5OZfxa0d0a+kV3M3+dku5t4TnUTXfpbYu+Vkl3t556W3g+ddtdtmMGAABAbrEDHQAAAHKLyTAAAAByi8kwAAAAcovJMAAAAHKLyTAAAAByi8kwAAAAcovJMAAAAHLr/wPTnYimLPzMGAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plots(x_imgs[:8], titles=preds[:8])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining Logistic Regression Ourselves"
   ]
  },
  {
   "cell_type": "markdown",
   "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": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_weights(*dims): \n",
    "    return nn.Parameter(torch.randn(dims)/dims[0])\n",
    "def softmax(x): \n",
    "    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": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = LogReg().cuda()\n",
    "opt=optim.Adam(m.parameters())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6f444429795b4e1489de25441ca39ea4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Epoch', max=1, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch      trn_loss   val_loss   accuracy                     \n",
      "    0      0.315913   0.282962   0.9215    \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[array([0.28296]), 0.9215]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit(m, md, n_epochs=1, crit=loss, opt=opt, metrics=metrics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression from Sci-kit Learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 23s, sys: 0 ns, total: 1min 23s\n",
      "Wall time: 1min 23s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n",
       "          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,\n",
       "          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,\n",
       "          verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "lr = LogisticRegression()\n",
    "%time lr.fit(x[:10000], y[:10000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 68 ms, sys: 0 ns, total: 68 ms\n",
      "Wall time: 22.3 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(10000,)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time preds = lr.predict(x_valid)\n",
    "preds.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3, 8, 6, 9, 6])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.899"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(preds == y_valid).mean()"
   ]
  }
 ],
 "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.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}