NB7_CVII-logreg_mnist.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook 7: Logistic Regression and SoftMax for MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Goal\n",
    "\n",
    "The goal of this notebook is to familiarize the reader with SoftMax regression (a generalization of logistic regression to more than two categories), categorical predictions, and the MNIST handwritten dataset. The reader will understand how to use the Scikit Logistic regression package and visualize learned weights.\n",
    "\n",
    "## Overview\n",
    "### The MNIST dataset:\n",
    "The MNIST classification problem is one of the classical ML problems for learning classification on high-dimensional data with a fairly sizable number of examples (60000). Yann LeCun and collaborators collected and processed $70000$ handwritten digits (60000 are used for training and 10000 for testing) to produce what became known as one of the most widely used datasets in ML: the [MNIST](http://yann.lecun.com/exdb/mnist/) dataset. Each handwritten digit comes in a grayscale square image in the shape of a $28\\times 28$ pixel grid. Every pixel takes a value in the range $[0,255]$, representing $256$ nuances of the gray color. The problem of image classification finds applications in a wide range of fields and is important for numerous industry applications of ML. \n",
    "### SoftMax regression:\n",
    "We will use SoftMax regression, which can be thought of as a statistical model which assigns a probability that a given input image corresponds to any of the 10 handwritten digits. The model is a generalization of the logistic regression and reads:\n",
    "\\begin{align}\n",
    "p(y=i|\\boldsymbol{x};W) = \\frac{e^{\\boldsymbol{w}_i^T \\boldsymbol{x}}}{\\sum_{j=0}^9 e^{\\boldsymbol{w}_j^T}},\n",
    "\\end{align}\n",
    "Where $p(y=i|\\boldsymbol{x};W)$ is the probability that input $\\boldsymbol{x}$ is the $i$-th digit, $i\\in[0,9]$.\n",
    "The model also has 10 weight vectors $\\boldsymbol{w}_i$ which we will train below. Finally, one can use this information for prediction by taking the value of $y$ for which this probability is maximized:\n",
    "\\begin{align}\n",
    "y_{pred}=\\arg\\max_i p(y=i|\\boldsymbol{x})\n",
    "\\end{align}\n",
    "\n",
    "## Numerical Experiments\n",
    "\n",
    "The reader is invited to check out the code below to build up their intuition about SoftMax regression. The following notebook is a slight modification of [this Scikit tutorial](http://scikit-learn.org/dev/auto_examples/linear_model/plot_sparse_logistic_regression_mnist.html) by Arthur Mensch on studying the MNIST problem using Logistic Regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Automatically created module for IPython interactive environment\n",
      "Example run in 47.833 s\n",
      "Sparsity with L2 penalty: 9.18%\n",
      "Test score with L2 penalty: 0.8948\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.datasets import fetch_openml # MNIST data\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.utils import check_random_state\n",
    "\n",
    "print(__doc__)\n",
    "\n",
    "# Turn down for faster convergence\n",
    "t0 = time.time()\n",
    "train_size = 50000\n",
    "test_size = 10000\n",
    "\n",
    "### load MNIST data from https://www.openml.org/d/554\n",
    "X, y = fetch_openml('mnist_784', version=1, return_X_y=True)\n",
    "\n",
    "# shuffle data\n",
    "random_state = check_random_state(0)\n",
    "permutation = random_state.permutation(X.shape[0])\n",
    "X = X[permutation]\n",
    "y = y[permutation]\n",
    "X = X.reshape((X.shape[0], -1))\n",
    "\n",
    "# pick training and test data sets \n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y,train_size=train_size,test_size=test_size)\n",
    "\n",
    "# scale data to have zero mean and unit variance [required by regressor]\n",
    "scaler = StandardScaler()\n",
    "X_train = scaler.fit_transform(X_train)\n",
    "X_test = scaler.transform(X_test)\n",
    "\n",
    "# apply logistic regressor with 'sag' solver, C is the inverse regularization strength\n",
    "clf = LogisticRegression(C=1e5,\n",
    "                         multi_class='multinomial',\n",
    "                         penalty='l2', solver='sag', tol=0.1)\n",
    "# fit data\n",
    "clf.fit(X_train, y_train)\n",
    "# percentage of nonzero weights\n",
    "sparsity = np.mean(clf.coef_ == 0) * 100\n",
    "# compute accuracy\n",
    "score = clf.score(X_test, y_test)\n",
    "\n",
    "#display run time\n",
    "run_time = time.time() - t0\n",
    "print('Example run in %.3f s' % run_time)\n",
    "\n",
    "print(\"Sparsity with L2 penalty: %.2f%%\" % sparsity)\n",
    "print(\"Test score with L2 penalty: %.4f\" % score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAFCCAYAAADyjdmjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztvXvQZVd53vm8rTuSkETrfmshQAIhR4JAikliwJhkZpIyDg6eYE9MSNmeco08jMkYYzvOlDPjlE1mPL7E9jgFvjCuYLDBIXHiCa4xyBCKGMdYxkg4ktClW61LSy2pWxd0X/PHuejZj77z9D6nz/5Od+/nV9VV+3xrn73XWu9a66x+b6taawghhBBCGCs7Nl2BEEIIIYRNks1QCCGEEEZNNkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXZDIW1U1V3VNVbB3r2jVX1Zvp8ZVX9aVU9UlXv0fKh3nsksUzdhpTNWNAxN8Dzf72qfoI+LyPftYzT7R4nVfXVqrpmu94XgnL8pisQwjK01l4tf/ohANe31l4z/fzzh/uOqroDwPe01v4/894jhnXWbau2b5IjrT5TdMwNyjLy5XuP0L7bktbaqzZdhzBuohkKRzu7ANy46UqEI4+qGuo/eyuNuQHrE0I4TLIZCitTVZdU1e9U1f1Vtb+qfmGLe364qr42NSncVFVvl/L3V9Xeafl/qapvPsTf5+r7qvo0gG8C8AtV9WhVXSHlC+u3qF5V9RsALgXwu9Nn/tAW731VVV1fVQ9PzRJvkzbdUVU/WFVfrqoDVfWxqjp5i775h1X1u/T51qr6Lfq8p6qunV5fWFWfmLbldjbPSN1eSyac356++yfQ5dqt6rZV2xfJYQsZf1z+9nNV9fP02dV/SzktqE+fvn9/VX0ZwGO6AZk+55Lp9fdVVauq86af31dVH9L2yfe3GnML63So+kzveU1VfWnaxx8DcLKU95bv7N5F41iee8j5O71vrXN4i+f/91X1qUV9HsK20FrLv/xb+h+A4wD8GYCfAXAqJgv4X5+W3QHgrdPrbwdwISYb778H4DEAF0zLrgSwB8CF08+XAXjZor/rs6efr8fEFAAud/XrUa/OO+S5JwC4FcCPAjgRwFsAPALgSrn3i9PnvwTAVwF83xZ9eDmAh6d1uADAnQD2UtlD07IdAP4EwP86feflAG4D8F9L3U6cPuN/ntbz2wA8BeAn+tZNZLdQDtKOXQAeB/BiGhv3AHjD9PPC+veQE9enb9/fAOASAKdsUdc9AK4CUAD+HMAtAF45/XwrgGt6jP3rMR1zh6pTj/rMZPbe6bPeAeDpLWS2jHzfqtfLzN+h5/AWdflJAD+96TUt/8b9L5qhsCp/BZMF8n2ttcdaa0+01v6j3tRa++3W2t2ttedaax/D5Mfnr0yLnwVwEoCrquqE1todrbWvmb+vrX6HqJfjDQBOA/BTrbWnWmufBvDvAHyH3Pfz0+c/COB3AVyrD2qt3YbJD+e1AN4E4FMA9lbVK6efP9daew7A6wGc01r736bvvA3ABwG8c4u6HT9999Ottd/BZOOjHLJuU3rJobV2J4AvAfg70z+9BcDjrbX/NP3s6t9rHFH7+vb9ntba17d4xsPTZ/xNTDYxXwFwJoD/BpON6J8BQFV9T1VdvaAey9bJ1ecNmGxsfnYqs48D+GPzrj7y7UPvft+GOXw1JhvTEDZGNkNhVS4BcGdr7Rl3U1W9q6pumJoQHsZk4TsbAFprtwL4AQA/DmBfVX20qi5c9Pd11s/V6xBcCGDPdJMy404AF8l999L145j8YG7FHwJ4M4A3Tq+vx2Qj9KbpZ2CieblwVtdpfX8UwHlb1G1va63R3/Zs8c5edVtSDh/B8xuA75x+nuHq32scTenb91u1ecZDmLT3BwD8HICDAM4C8D+CnO9bax9qrX1lTXVy9dlKZncuca97tqN3v2/DHL4ak01pCBsjm6GwKnsAXLqVD8SMqtqFiQbg+wHsbK2dicmiV7N7Wmsfaa39dUx+MBuAD7i/r6N+PerV9DvE3QAuqSqeO5cC2Ltk/WbMNkPfOL3+Q7xwM7QHwO2ttTPp3+mttb8lz7oHwEVVVfS3S5asT6ftS8jhtwG8uaouBvB2dDdDrv6HGkdcn7597+T3MIDXYWLmuR6TzdA1mPwgf3J2U1Vdb57B9KmTq89WMrt0iXudfN17Dzl/geHncFWdNm3DTa4eIQxNNkNhVb6IyeL8U1V1alWdXFV/Te45FZNF8H5g4jCMyY8Opp+vrKq3VNVJAJ4A8HUAzy76+xrrZ+sF4D5M/Fq24o8w8Zn4oao6oSY5Xb4FwEeXrN+MP8TEIfeU1tpdAD6HiclmJ4A/pbYcnDqknlJVx1XV1VX1ennWFzDpp++vquOr6lvRz/THzNu+jBxaa/djotX6NUw2Pl+lYlf/Q40jlsU6+v4hTPxzZlqgg5j44PzL1tqz03afjon5sg+HW6cvAHgGwHumMvs2LJbZsvJ147jP/AWGn8OvxmS8PG7aEcLgZDMUVmL6w/EtAF4OYDeAuzBxruR7bgLw05gs4vcB+AYAn6dbTgLwUwAewMR0cy4m5pNFf19L/XrU6ycB/NjULPCD8tynALwNwH87rd8vAXhXa+0vlqkfPe9mAI9isglCa+0gJs7Fn5/9OFNbrgVw+/S9HwJwxhZ1+zYA342JBuTvY+K/8uQSVZq3HZP+WkYOH8HE0Ze1Qrb+PcYR1+c9OPy+fwgTv5tZHQ9i4jP0Qbqnt9nmcMcDyezd07r9PQC/c4h7+8rXjeNDzt/pfUPP4ZjIwhFBdc3PIYRjiar6IwC/3Fr7tU3X5Wihqv4HAI+21j5yyJs3zNEu36r6OQAPttb+6abrEsZNNEMhHENU1Zuq6vypGeUfAPhLAP7Dput1lPENOEKjm44l+U79hf42gM9sui4hJCNqCMcWVwL4LUwipr4G4B2ttXs2W6Wjjm8AsJLZcxs4JuQ79a36VwA+jKmJOIRNEjNZCCFMqapPArihtfbjm65LCGH7yGYohBBCCKMmPkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXZDIUQQghh1GQzFEIIIYRRk81QCCGEEEZNNkMhhBBCGDXHL3Pzaaed1nbu3DlUXcIh2L9/Px599NFax7Miy82yTlkCkeemydw8dogsjy127979QGvtnEPdt9RmaOfOnXj/+9+/eq22ieOP37pZzzzzzDbXZL184AMfWNuzIsvNsk5ZApHnpsncfJ7I8nkiy81z3XXX3dnnvpjJQgghhDBqltIMrUrVYo1ja63XfVp23HHHza+fffbZThl/5u/t2LF478f1OBzW9Zxjjb6yfe655zpli/5n8tRTT3U+s8xPPfXUThnLXceKYwyydHJx97Fc9H+VJ5544vxa+5A/8zNUzjwODh482Ck75ZRT5tennXZap4zXBYWf6WR7tMq9ryyXeYZbZ3WuzlhmnV21zn1ldCzKsu8YVznw9/qOf33GySefPL/W8cDfWzQ2tsLVmRlaSxXNUAghhBBGTTZDIYQQQhg12QyFEEIIYdSszWeIbZyPP/54p+xFL3rRwu+xjVDtmGzjVDspfz7hhBMW1sX5JC26D+jaQ519dRnbaF9/lWWeOTTr8EMAfB86OSzy83C+KOpPxD4s6t/i3u1s4EeSjA4Ft9H5hDzxxBOdMrbROz+FAwcOLPyesihq5emnn15Yr5NOOqlTxvNd28OydjJyvmlaF563WpftZtX52Nd3UuF79XuL5o6Olb4+o4rz8VrH3Ny0P5Gbl4xbE/W3z/1m8thVWfLnvr51Ok9cv7vfU+fL5PyJ3O/BKkQzFEIIIYRRk81QCCGEEEbNIKH1rKoGvIqdcWYsVa/zZ1XPLVL5Ka5erH7U+/qq31WN6FShfVXC242+25k1HSwv7Qd+pjN5sqpVTTF8H4deA91wUGcmU7ny+9SE5EyeXDdV7aoJb9NwXbmfgG4btW8ee+yx+fWjjz7aKeO+0u99/etfn1/z+FGZuXHG95555pmdMmdecyZRN9+ffPLJLZ+x1XM2CY/thx9+uFPGstU2cPu0PTw+WHZAdyzzu104ts7vRfXXe/WZLFv9HtfZuT8oR5IsmWXSHfBvr85nN8YXyc+ZILWMf++0Xlzm0qK4cH0dO+t2Y4hmKIQQQgijJpuhEEIIIYyabIZCCCGEMGoG8RlyvjjOL8iF+y0TUsq2bb52oYDu3ep7wPZVZ4NWu+WqRx9sN33DzR06BlyILqO+OfyZfRu0b9lWrjJhnxZ9/iOPPDK/dm1VvwT2VdEjIbgu+r6Z3XtIGS/zbCdr7jfue6Dbb+w/BHTnnLafy9h3wIXb6vzjvle5sD+R+kxwig/1a+Qxo+1xvgnbPVfdmsUpDjS9Cfe1+tuxjJYp4/dzX2sqFe4jlSV/T9vDstRnsvx0Xed6qj8b94OuE5tcd7UNjNaLx4COfy7T9rmUFFzGfasycX6Uzl+W+12/x2XOz1b9Cl0ajVVSYEQzFEIIIYRRk81QCCGEEEbNIGYydzKxM6Gpyo/vXeY0W/7M6mINDV0U5gt0VXCqnuPPWuZC5F0mbq6zqjCXOWl9HfTN7Kn1ZPm50E0NkeQyl0Kh77vVbOWyU3OZmoL6Zl/V77Gcl8nGOxR9Tbl6H88dbaMzNziTDL+vb2oLxYV7s1xU5e7kyeuLmtBYrb/p0HqupzNBqpmMcePBrZGuLoxb/1Wup59++vzapWRx4eX6TG67mkqdvLZ7nWW0fTzGXUZ37TNeB3XO8tx78Ytf3CnjfuLrZdLKcL+redKV8Tx1KUx27tzZKePfXjWLrSLLaIZCCCGEMGqyGQohhBDCqMlmKIQQQgijZhCfIeffozZHtvcuk46d7Yx9fY3cibhqq3T+Ei7kn+usbeXnqP8QP1NtxNvtl8BtWOYEchcOzXboU089deH71M+D+0z9sxguO+OMMzpli8K5ga6c1T7uUtuzX5ILJV4U9jukTF2Yuo5lFzrtQl5d6gluv/qcLDpGwR0doLBc1D+M6+WOBHC+FtpWlyphlbT/hwO3Sccyy8/NI3cch44Pfo7KZFFoszu2w/mxsP8Q0F0ntJ/5mTr+3Mnt/HkZn7Khcb6Ezk9tUeoOwKeE0TWY+/4lL3nJ/Fp9i7j/9MgXfqb+Xt93333zax0D/FnnLMvdpXZwxzP1JZqhEEIIIYyabIZCCCGEMGpWNpM51Z3LMu1OLFfVGqvkXOZXVa0tOm3bnUquKmdWP6o61Z3Q7E4Udqq7VbNtrwN9H7fXhcm6cFdVT7t0BKza1bJF2YxdSLWqYVmFqvVy6Q6caY/NM2pycSkbnAn5cJnJQ+Xi6sNmQp0ffK8zReg4Zzlp3yzKVOxOrVcVOL9bzWsuHJvnuMqBzQQu9YN+z5mRD4dZ3TVEnuv50EMPdcoOHjw4v9Ys2iwj7TNnDua261zhMcF97cKqdW7u379/fu1M2FovnY8M36um/E2uszp2XDbxRSlFgO4Y1/HnzGQsP5Ulr2d8raZLrpf2H8vr/vvv75Tx/HamZe0HXoNcBup1mDijGQohhBDCqMlmKIQQQgijJpuhEEIIIYyapXyGnnvuubnN14U5uyMIXDi280NSm6ALB+UTtdl+rfZ3treqHZPrpbZstse78EXF1Zlts+sIEzwcnP2171EdWrbIj0txY4B9X/bu3du5j5+padvdac0chr/MCehsV9fxznV2KSGGwh0to+POzT/uD53vfK/Kk/0RtN8WHb/g5rc75VpxobjcHvVz4nvVz8kdbTH03NT1hceP881Uvyd+joZL81jWsrPOOmt+feaZZ3bKuJ/YX+n222/v3HfXXXfNr93RLfwMvdeldtA55o6TWCaEfV3M6reqv6DKkj+7Y49c23UO8e8kjxWdvywHDa1nPyH1Z+Pv6TO5Xjr3GPUT4zVJx0eO4wghhBBCWJJshkIIIYQwapYyk+3YseMFats+OPUjq7NUVe5OFGeTCZvFAODAgQPza6fiZnW+hhCyal6fwSo5LeO2qprShXizOnq7mNV1GXMXt8+ZwlwGXBfu6syMrMq98847O2U8Hnbt2tUpu/jii+fXagZgtaxTR6uKm/vBhcVqe2Z9NqR5Rd/JstC6LsoiDHTnhPYbq6hd3zjVNr9P5zer3N380zHIc9plnFcTGpe5+bDptBfOxMyyVLMmmyZ0rbn00kvn1xdeeGGnjM3Iat5gme3Zs2d+fcstt3TuY3OKmmdYXjpW3NxkXDi2W6O2K+P07D3adp6nrn0uM75zK1DYdMW/kUC377mP9N1cpubQG2+8cX6tqR1cWD+3R81kvN/Q8cFy1z3AKpnhoxkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmpWP43D2eLUz9v2e2vnYN0CPX2D/kQcffLBTxv4+bIvl03gB4Oyzz55fq62SwzzvvffeThn7MKhvDNsu3XEcaiN2oaJD27a1ni5l/ar+X+xfoH3tQiv5OWxzv/XWWzv3cb+fe+65nTL2d1F/CX6+C/t1Zeqf4Y5omI3x7Ty13vmxcVi1jjueH9qn7D+iMuOxrWWLjvFQHwP2p1B/Ig21Z3g+qu8Uv0/nH3/PnZC+3afUu/nn/KV0THJYPPsIAcAll1wyv2aZ63N0DvCazD586jO0e/fu+bWOv3POOWd+zWMR6MpPy7g96svEMnIh185HdQi0Li6tBZfpWHVHIrlUBe7oHfYPZD8x9RVkvyP2EQK6cndHLjkfSx0f7mgVxh110pdohkIIIYQwarIZCiGEEMKoWdlMpjjTCofHqQmBVa2q2uJnuizTqn5k1R6r9y+//PLOfawSVlUkqxhV5ccZVTULJ6siNZSS+8WZVlbJnrkKM7XwMmZNxpkBtcyphHl8aIqDRdmTb7755s597nRtVqmrnFleGsLN5hgnE30fy1bbOtQp54zLpu6y0qoJkcOsNfswq6813UbfDMDcpzqPeO5o/7os3iwz7WseSy4TvprseC1YlCphKJbJys/jTmV53nnnza81QzuHLDvz4QMPPNAp+9KXvjS//sIXvjC//upXv9q5j2XpzHBqCuN6nn/++Z0yHo8qZ5afjhV36sDQLONywGg9uQ26ZnHIvJ4cz3NM5yynH3G/u/xMTW/CJjQ3712qAC1zJk+ei2o6XyUFRjRDIYQQQhg12QyFEEIIYdRkMxRCCCGEUbO20HoXGs72O7Vxss+Qs/Npum1+jvrfXHDBBfPrq666an7NdnOg64ugtnK2o6vvD79b/Qv4s7aV7aEa9siftzu0Xul7srNLda+4FO8sB/Ub4L5ge7j2H/uDaRg4y1LtzixblSXbodVuz/3gnql9Muuz7ZbpDNf3HOYMdH02NBU+zxf1P2BbvjtOgvtA68VyV/8X59/Dfa/14nVC/ZD4mZrGg8edrhOzZw4lT+dfp3OFfS+0jOXnjgnS08Z53f2Lv/iLTtnnP//5Lcu0/3gcqb8S+xDxug14PyceH33XHaA7/nSdHcpXczY2XAoThceZ+kSxTDStDKeB2bdvX6eM10z9LVwUtq6y5HdzehugO/ecr5v6IfWdO7rOrvvYnGiGQgghhDBqshkKIYQQwqhZ2UymKj+XHdKFArrQelaDqVmE369hvy972cvm1xwerKp+VpU7lalmOOXvqcqZ6+zMhfo9Vh2qen+Ik7JbayudoM7tcyGfLlxfTQ2s3te2sxmEVfhqBuCT6jUM12XRdaG2LozZ9Vkf88+Qp5+7LO86JvmzzuG+JlE1I7M6W1X8XMbXagbnZ2of8vjRscTjQk1hnHJDn8l1cZmWNVx5Zm4YKkzbmQLUPcClr+D6qbzYDKimiPvuu29+rSHzd9999/yaxxHPRaBrClOZcLZjd2K59oNLUeHmFtdzO9JcAM/XR39/HGyiV1MVy0RPR2DTGMsH6MpW108eE2x603fv379/fq2/5fwMDXXnZ6p5jftF1yf+7dW53ve3ti/RDIUQQghh1GQzFEIIIYRRs7SZbKZ+cmYxVXeyilZVa2weUtUaf0/VdZypWKNgFh3+p2plboM7JFbVqeoNz7BqV1V3/Ex3UOuiCKR1UlUrqRIXRQP1ed8MFzmkbeesqXytplGWOZtGga4aVscYjwGnanV95cyFm2DVuqr5hE1XOv8WHbgKdPtY+5vfwfNI38311P50h0LzuqDznT+7LNb6TDblaPTOTI0/lNlT5cXj1R1EqzLh9ur6xfdqGZs0dL7zwdccCaZZ5Ln/nLlVI9m47S5yzh3u7A4F1b4dOvO/W190nnC/OzPZnj17OmVsNtNnsslJxz+Pa36muojs3bt3y3oA3d92dXfguvB4ALpy0Pc5dws3PtxhzouIZiiEEEIIoyaboRBCCCGMmmyGQgghhDBq1paBmu13astm+6c7RVhto+xHoPZdtklqxmEOoWV7ssv4qTZGDuXncELghVlvGfZPUV8HtkmrnZvr5sL1twN+X99s1EC33s5Wr+1hPwW1299+++3za7ZJq22ZfRbUn4jf7cK03Unmrh/U/8tlj52VbWdoPbdLQ9idbxyHQasPAOOy1KqPzaIwbuffo3JxGcv5e1pn58vFdVHfFa6z+jbNxu665+hsfLiTzrXfuQ4u3YGupTzONfSdMxXr3OQUFjzmnW+mjr8HHnhgfs39DHTXWa0Xz3GVs1vnuf+2KwP1DB1/XBf11eL+1PHI9+ozeQzoGsl9pmsW//5x+gj9Xbzjjjvm1/fcc0+nzKUwcX5cLoUJt8+F1mvqnVWIZiiEEEIIoyaboRBCCCGMmpXNZKq+dZlsWWWmarBFobb6Wc1K/FnVgYxTjXO9VEXLpjE98I7Vt87sonXmPlNVMvef9tHQh0Auc+Ad19Op8LVfeEy4gwd1DHCYJ6ts1Txy0UUXbfkurZcrcziToEsVoMzaPqTp06W9cHNMzcHcZhfyqmZjVllrxmY2izBq2uTwW2cicW1Vk7wL5ecyne+33Xbb/FrDxocye87GhzuIllNNcF2AF/aLM9/z91TO7HKgh6Xy+unWBTZvayg4h2drZuK+h7Hq+HAZ7d3cHIrZXHFt0HBzlpGudXxwsZo8WV7uxAWVM48BloOGz7Pbgr6bxwenXdDPXEegKyP9zXSpT/g3VNc1lzpjEdEMhRBCCGHUZDMUQgghhFGTzVAIIYQQRs3SPkOLTjp3NmNXxmiZC7njZ6p9kH1LXLg336d+Amzb1vBgrovzZXIhzn3TjAPP+zOtO/RzVj99rrPTMu6UaYVt0iov9itRnxK2WbMdXU+mZzu0C2HVMv7sxpjixpU7tX6rOq2L2Xu1Hfwu9ZVh/wA9KofHofNx0xBvfo7OK55LPH70+fwMDall/wCXKsCFamvIMI8DDWXm+aHynLVn3UewzGSo7+O+Vn8sd8I3f9YynjvqS+L8MRetz9oX7igQ/qzj1qV9YPnp/OM+c0dD6fvcmnw4zOruft/Uv4fbq78HLD/1v+Fnqi8V+x6530L2HdTfBvbXUz+nV7ziFfNrPRKJv+eOT3HrokvJoqgPYh+iGQohhBDCqMlmKIQQQgijZuVT693fXSZfVd25E5NZdedUZBpiyqpDl3mVQ4A1myaroJ2a0oUcr2om0+/NVH7rVONW1dyspSYDlqXLVuuyBissBzVtsOlGzWQsZw7PPPvsszv3sVpZQ3QZNelw21Ul7E6td+pcbo/KcvY+7dd1MFN967NZBa/qcZaL9o0r47Gs/cbyVZPoogy5OnZYpa/rCT9DTXssT2c+17WGTWhaxuH0u3bt6pTNVP4uxH8VZs/T8cNmRp1HzlS3yHUA6K6XanbhtU6/x3OA66njz5n2eA3WMcbv1jHm1nUeS268D51xGpj00WyMuHBvlR2PJ01pwKH1l1xySaeM+0x/a3kdUHP53XffvWVd9HQHXne1Xpdffvn8+pxzzumUsWlM5eX6he/V8eFSL6ySkTqaoRBCCCGMmmyGQgghhDBqshkKIYQQwqhZObRebfXupG62LWv4Ndsu1cbZN7RSU4Yz7JegdlK2uWs4LT9fQwG5rS6cXG2j7sRkZjtOqW+tzX0jhggpdeG1esQA24LVX4Pl506mZ9mqLFleakt2vgd9T1NW+DkunHbdzN7rQv11jnFf6Thw84/L1ObP9nsN414kQ/UxYP8GfQa3T+c0z1tdo7ie7rgKPf7AyWw2Ptd9zMOsv1VeW90zw51mzp+1ru4YDw6fdmH23Lfqs7d37975tc4/vlfnH6+7Ojb5s9ZrmTWZGWJuttbmY9TNSx2r7Kem7XPHavBn7QeWrfrf8DrFZW79Ur9N9jdzv5nOR8itszqm1+13Gc1QCCGEEEZNNkMhhBBCGDUrx4NqKJtTZzGqdmO1nqq9+LOGkWo4PcOqV6dmZhW7qu5YpajZQVn9qGo97hdtq8vsvGom5HWgz3dt4LppnzmTE/e1hvKzel9NIqyKdWYrfreaybjfnYpW+8Flq+V7nfpWnzlE5mlF38l970y3qtrum7lb1wKec5qlltMjcGiu3sd10efz/FazC39WWbv0Hy5M12XPHRqXJkHbwPV2pktnXtAyzvTu1lKe05oeg81kuo6zTNTkw+ZKLXNmTbdeujm93Ti3CZeVmeWu45HTCqiZzGWi53WX56U+n+e9psPg57u0BTrG+JlaL36OW4PX4Y4QzVAIIYQQRk02QyGEEEIYNdkMhRBCCGHUrOwzpGnb3Snyzn7o7Phsk3QnY6uNmt/PIYoansl2RfVVcbZ5bqvzGdJnuuMd2I6qPjXumIt1sIx9lz/rqeAuJb87hZnlwqnmga78+Bn6budf5kL59OdzAAAgAElEQVR0nUwY50+kcD0XHU0zhO/Q7Jnq3+N8rVzfOL8gd4wOj1ddJzi9P/uDsZy1nuofyKkZtD0sF+17/qxl7hR09knRuTjEsSrA87JUHwqut/qE8Lqnaw8/R/uT2+D8s9yazzLRuemORuC+1Tpzmfptso+Z+jI5H0HnSzK0D5E+n+eb+z1wfjQqL77XrVnOn4595HS8a1oUhtcLTRXgfjPdETrcPrc+KfEZCiGEEEJYkmyGQgghhDBqVjaTOfW0yySq5i4XQuhOSHcZcFkty+o6d4qvC0N05i5V3a0a1snf0/cNzTKqY1Zbqvp70SnW+lnbxyaSCy+8sFPGKlseO6oi5XB6Ndvw+NMQXdfXzpTlskwvMh8ALzQHDYHWx6maeU6406RVHc/3OtMmh9IDXTMZy1bl4FIguDQNzrzAstc+4n7QOczrhI4fNQesGw1t5vdr3/JYc+HYKkvuQz69HOiavHSNZJk5MwiX6fjn8ajrs8u0zDLR9/H817byvcuYXYbAmSddSDm3z50eoHPIvY9xJ8xzXdw4cvsDrbPO4UXv0++pqftwiWYohBBCCKMmm6EQQgghjJpshkIIIYQwatZmdHO2erY1uxTeagNk26WWsQ1Z/TLY9utOmOe6aL34XvVncCn5nZ+FC4lk26jax7fblu3stCxbF8Lt/Gi0jPva2fHdkSXcn84GrX3r/KWcnwXXU8t4vAztU8LM2qL9604l53G/TPoAfofOKw7B1pBoDpfmeeTS6asvCYdV61rDz9Rx4NJesN+TO5ZET2Qf6qicReOSxx37XAFd2Wrb+XvaBnfsCofh69hZ5NOj8uK+1nWWy7TN/D6XXsQdDaV+My7Fx9A4P0rFhcGrTyTTN9zcrc98n84vfreubTzmnC+f83PS9YLHgOu/daRFiGYohBBCCKMmm6EQQgghjJq1mclcWLU7ldyFNrM6TVVrzuTEakWui6q0nRmO3+dOdV+XqnW7VbYOZ7pkNaaqWl0IO6vOVY3O4bZqcllkXuPTtIGu+UVVu86050JYXTZXRscmj/ftCKVXtB3cb9r33C41D6ksFn1P73NZrbluLgSa5eTMnnpiOdfFZUF3Y3f//v2dMmcudSe5DwGvlxoG70yLnJZC5dzX7OKyA7PJROXFZk2XFdyF7uu7eW46E4mbt+s46fxwcKfW89hV+bjwdr5X+4zf4X7v+Nr9LuncdnOPP2tbnZnWjQ8e09qeZKAOIYQQQliSbIZCCCGEMGqyGQohhBDCqFlvPuspy4Qhs61PbYkuHFt9Upi+J4K7U9b5s/Mzcc9092lbt9te3RdtO4e4Ohu/ysfZfhk9roXfz32m8nL+Bcwyx8g4WfY94X7Tfgn6TvU/YDm5VAnq28dt1L5hGbo+ZZzPkH6Hy5zvg+L8X/h7O3fu7JRx6g43p7ebK6+8svOZfZ203116E+fP4cKXuS/YD0n9lXiMubQFKhP2PXJHL7iQfJfC5Eii72+KousLf3Y+gO5YIicT9oF0MtH0Dfwc59P5yCOPdMqcv6A7NmQVohkKIYQQwqjJZiiEEEIIo2ZtZjKnymMVoKr1nEmDcepAx6qnyDPLZA51OLMct+dICrN3GahdCgXNCs7tc7JzZiXXR67PVlVB962X4sx5m8CZ+zgkWueKS23h5OlChhfJ3pkv3fxzZpBlMtaqeZZhE4Izw20H3J9q1uTM34qbc65NTpaLzMFurjh56Vxx88ilTHG/N0fSOusyQi/KvK8sM/d4fe5reu6bqRro/pY71xJ3MoOa452LjXNxWIVohkIIIYQwarIZCiGEEMKoyWYohBBCCKNmkNB6hW2LfX2EgOV8S1bB2WyZvn5Hy7zvSA3xVJaR1zIhoH3L+vaZk6Wzoy96xqHK3PEzR7Jsl/Gjcb5Wq/qgLJLTqmPHPb9vWL97xpGMW3uW8XHrmwKk7/vWJUvnn9L3OZv2C+rLMn3GfpvaPi47cOBAp6yvb1jf9BR9/Tu1ns7fc5k1eN1EMxRCCCGEUZPNUAghhBBGTS2jeqqq+wHcOVx1wiHY1Vo7Zx0Piiw3ztpkCUSeRwCZm8cOkeWxRS95LrUZCiGEEEI41oiZLIQQQgijJpuhEEIIIYyao2YzVFXnV9VHq+prVXVTVf1eVV1RVZdV1VcGeudJVfWxqrq1qv6oqi4b4j1jY0OyfGNVfamqnqmqdwzxjrGyIXn+o+m7vlxVf1BVu4Z4z9jYkCy/r6r+vKpuqKr/WFVXDfGeMbIJedK731FVrapeN+R71sVRsRmqSdKCfw3g+tbay1prVwH4UQDnDfzq7wbwUGvt5QB+BsAHBn7fMc8GZbkbwLsBfGTg94yKDcrzTwG8rrX2lwB8HMA/H/h9xzwblOVHWmvf0Fq7FhM5/l8Dv28UbFCeqKrTAbwHwB8N/a51cVRshgB8E4CnW2u/PPtDa+2G1trn+KbpbvdzUw3Al6rqr07/fkFVfXb6P4+vVNU3VtVxVfXr089/XlXv3eK93wrgw9PrjwP45lpHtsdxsxFZttbuaK19GcDRkYXt6GFT8vxMa212yux/AnDxgG0cC5uS5UH6eCqARPWsh039bgLA/47JxvaJoRq3brYlA/UauBrAn/S4bx+Av9Fae6KqXgHgNwG8DsB3AvhUa+2fVdVxAF4E4FoAF7XWrgaAqjpzi+ddBGAPALTWnqmqAwB2AnjgcBs0YjYlyzAMR4I8vxvA/7tqA8Kcjcmyqq4D8I8AnAjgLYfdkgBsSJ5V9RoAl7TW/l1V/eCa2jI4R8tmqC8nAPiFqroWwLMArpj+/Y8B/GpVnQDgk621G6rqNgCXV9W/APDvAfz+Fs/bSguU/7VsD+uWZdgsg8izqv4+Jgv3mwatfWDWLsvW2i8C+MWq+k4APwbgHwzdiDBnbfKsqh2YuJS8e7sqvy6OFjPZjQD+co/73gvgPgDXYLJAnggArbXPAngjgL0AfqOq3tVae2h63/UArgPwoS2edxeASwCgqo4HcAaABw+nIWFjsgzDsDF5VtVbAfxjAG9rrT15eM0IODLm5kcB/J1VKh9ewCbkeTomGqnrq+oOAG8A8G/rKHCiPlo2Q58GcFJVfe/sD1X1+qrS/w2eAeCe1tpzAL4LwHHTe3cB2Nda+yCAXwHw2qo6G8CO1tonAPwTAK/d4r3/Fs//D+UdAD7dkqXycNmULMMwbESeU1X8v8RkI7RvgHaNkU3J8hX08W8DuGWNbRoz2y7P1tqB1trZrbXLWmuXYeLP97bW2n8eponr46gwk7XWWlW9HcDPVtUPY+KUdQeAH5BbfwnAJ6rq2wF8BsBj07+/GcD7quppAI8CeBcm/kC/NlXrAcCPbPHqX8FkR3wrJhqhd66tUSNlU7KsqtdjEllxFoBvqap/2lp79TrbNkY2ODf/DwCnAfjtmsQ07G6tvW1d7RojG5Tl90+1fE8DeAgxka2FDcrzqCTHcYQQQghh1BwtZrIQQgghhEHIZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmmyGQgghhDBqshkKIYQQwqjJZiiEEEIIoyaboRBCCCGMmuOXufm0005rO3fuHKou4RDs378fjz76aK3jWUeLLJ9++ukt/37CCSdsc03WyzplCRw98jxWGePcPFaJLI8tdu/e/UBr7ZxD3bfUZmjnzp14//vfv3qttonjj9+6Wc8888w212S9fOADH1jbs44WWd53331b/v28887b5pqsl3XKEjh65Jm5eWgiy80SWT7P0S5LALjuuuvu7HPfUpuhVanqt8l29+3YsWPh59Zar+c7bYI+g+vS9/kA8Oyzz/a6b5lnHkm4enPbdRJx3x88eLBTtm/fvvn1GWecsfD5PGGfeOKJhfe9+MUv7nx++OGH59fPPfdcp+ypp56aX+v4O/XUU+fXxx13XKeM79Wxeayw6rzl/tD+1s9bfUdZZm66sqN1zvXlSFlnTzzxxM5n/t661lkeR+57R6vM+8pymWfwGta3XxZtkoAXzmU37x196zK0LI/NVTyEEEIIoSfZDIUQQghh1GQzFEIIIYRRszafIbZPOnunsxmrXwajtmy2ZWoZP6dvvdTGyf4vrsy1R/1mXFvdM7cbrov6WT3yyCPza/apAbrte/zxxztl3Cb1q+J+uvfeeztlX//617d8hsqcy1TO7BeksGy1Paeccsr8Wn2ZzjrrrC3vA4CTTjppfn0kOCCuw/9Axyv3m/Mz6etPteqYX/V7OqePFr+TvuuZwnJYxmeI57+OgUXPdP512pcsB+dnovOI79Uyt3YvesYmWHVeOn9FbpPKgT/r9xaNq2Xq6HzD+LObe+63tu+7VyWaoRBCCCGMmmyGQgghhDBqBgmtV5UVm1bOPPPMTplTt7MpzIUJaignP4fvUxXck08+ubDOrJ7TMG7+7EK19ZlshnHmAw1n3KQ6V01MJ5988vxaVZiPPfbY/FqTJXJfaJkL32Q1Pfe7kwmPN6Ar5wMHDnTKuA1s+gKACy+8cH7Npi8AePDBB7d8BgBccMEF82tV4Z922mk4kuD50VclDfSfm+sw0TmcerxvWP+hyvq+b7txJpJlTJcsL52L/FlN5vzZuS30NZXq+OO+5jkMdNcahb/n3BGUTa6zOk/6jjP9npODM3nyZ2dOc64sXGdd91wZ/x4sSrK7FU6Wq8zTaIZCCCGEMGqyGQohhBDCqMlmKIQQQgijZmWfIecnoLAvhrNVqr2afYHUnut8hti/g+2kandmm6P6knBouIaJM2qb5FBw7SN+v5bxERLaR7P2DO1/MYPl4FIHKOxPpG1w4ZPsR6O+OVwXZ9NnGanPEH9m+QDdMaDjiOul/hLcntNPP71TxnZvPRpkNs6W8c85XFzfO5zvAH/WebvIl8SV9fVFALo+By6sWsv4s3umymY7ZaXofOOx5Y5B0f50Mln0DKA7XnTusN8cr226zvIzdY696EUvml/r3Ofvqe+iS6XB7VMfFH6OPnP2nKF8h5ZZv134vCvj/tW+5jHgfjO5zPmQOT8n7VueXzo+nEzcuGJ0jq7yWxnNUAghhBBGTTZDIYQQQhg1awutd5kvF5mtgK7azWUmdiHsqj5jda47SddlzORnqpmMw7hVdexUeax+VPUm95H2g6o0txNnHlIVqgvTZnW1M7moyWmROVRNIE4m/NmZVVyKBh073B4d725Mz8btpkK0XZizC0V3IbZOjc+mU6BremQTicta7OqlZhD+rHORPy/zPZdBebtNaDwf3Vrq5qbLAu36TF0J2EzGZRr27twDeN07++yzO2W8Fug44u9pGaPycS4Ps0zy2+WO4HApDZwZi8u0X/Q5i963yGSmz9AyZ47n32ude25e8tr96KOPdspYli5cvy/RDIUQQghh1GQzFEIIIYRRk81QCCGEEEbNyj5DLjTV2Sqd3VL9Btjup7ZeZ4NcFPar9lW2jS6TKqCvf4o7ksL5HvBxDvy+TaSL13ouStsOdPtd2879q3LgU99daD3bqBeFxQIv9Dvi77lwUPUvYFu2ypnbp/Zq5+8y85npe0TBOuA26ntdKnzuD+03bpfKk+e/PnORD5j6H7hT0F0YvJtjjDuqQNvq1pqhfYa0De70eXcUA5epXxCvZ+7Yi4MHD3bKOGWF8wtyY92tsyxnXf/5feqj5MYtP/OMM87Ysi5D+Qw5WTofUj7GCej62vE14H35eM7q+vnQQw/Nr/fu3Tu/vu+++zr3sRz03ZxC59xzz+2UnXPOOfNrXu8B7wPV18dRx8Aq62s0QyGEEEIYNdkMhRBCCGHULG0mm6kQl8l+quo0xpkb2DShKlpW7arKb1GotmYD5jI1Z7iMrc60wnXRk9Xvvffe+TWHpQLdfjj//POxnTjTkQub1n5xIc9sOlGTCH9WMyq/z5l0+Plnnnlmp4zlrGY4Z4plM4DWy6V24DGhY2A7wnadOVXNfS5br8vCzm12qRJUxc9z8CUvecn8Wk0WLvuwGwcu5QF/1nHgTHY8H5zpe53MxokzI7l0BC4zvo5XDll2GYB1vrP82PShZmpXZ66nZo53J51zPXVsclv1fVxPnQsuRP9w6CNLHXM8b3RucD3V5MS4tAJ33313p+yWW26ZX990003z6y9+8Yud++6///75taZCeMMb3jC/vvrqqztlPDbZnAb4TP8ulN/9Rus47kM0QyGEEEIYNdkMhRBCCGHUZDMUQgghhFGzcmi9C6te5pgGRu3CbON0fgNqA2cbJNtU+e9A1y7rQhtdqgD9nrNzc/ifhixyXa688spOGdvm18ki/y/G+VLpGHBHsjh7OaM+Nu4kZIZlpD5q7EOkdmd+pvpZ9E3foG1z/myz8ThEaP1s3u3fv3/hPe50Z/W94M/qs+fS3XM/qm8czyuemzrOuA/1Xby+qKz5s/qO8fzX93E/OL8156czBM4304WNu9PanZx1/jFuXWc/kJe//OWd+9iHyB2bwD6VALBnz54t6w9011KXdsWFs+u6PvOdWne6hD5+tsscueF82Lh/ta937949v1afoQceeGB+zXPjla98Zec+/qw+Q5wSRn9rua3uSBSdz268rzvVTDRDIYQQQhg12QyFEEIIYdSs7dR6py52akuXSZRDPlV968J++R2slnch/hoCzCp2rTOHvqsJ65577plfqymM23rHHXd0ym677bb59Zve9KZO2XafdL6q6p/Vy87kpCphp9pdZHJSFfeFF144v1ZZssrZmTVV7crt0XqxLLWt7sTpoTIWt9ZeYC6e4UxaListo33qVN0O7is+6VxPROcQXp37bHbRNBRsrlGTqDNtuqzWLqP4rG+HDrHf6rOOSa6b9hmHkavsHn744fm1rsE8z3TOsSmE+1rXWb5Px6jLZM7t4zoCXfOr9gPLT8ctm2Z1vM/uXbfpczY2nCwV7heVl5vPLGf9/eE5pXLm38lrrrlmfr1z587OfSxLrb87fV7TJjAuJYvOYcalkkgG6hBCCCGEJclmKIQQQgijJpuhEEIIIYyatfkMudOi2X6t9ni2XapdkX2GtIw/Oxs421DVjsi+B5oinE/ZVVs5t+G8887rlF166aXzaw5XBLpp1bXOX/7yl+fXau+f+Z1s12nK/B53HIc7isT5XahdmG3G6lPA/cT3adg0+waoPwi/z40/DR/nk5x5LALdtmofqZ8C0yedweGi7+d2OZ8l/Z6TC88Jl35Bj8BheJ7qsQLsw6D+RA6WhY5r1x5+v/PF0bbO/Fy2a24uqhfQXev0qAn+rGsp94U+k+efjh3+HvsJaWoHboP2O9flrrvu6pSx3LU9fVMAuHGl6/rQvpnLyNKlDnBy5rnu0mGoL85ll102v37Vq141v9b5y2NAf9/4feyfBHTlrkfhcBtcqgytMz/T/d70JZqhEEIIIYyabIZCCCGEMGrWZiZjVB3Iamen7tRQR/6sql1Ww6lKjlXlrErTMHhWG2oINKtQNZsm48JI1fTGZgh9Jt+r35v137rVt4tCPvk9LhxUVd78WU0NfVXELjsu94uGfLK89BnuVG5WR+s44iytaibjOusYcKp4Hf/roqrm5kFtowsNdxl5WYYu+7CqvTkzrY5lHgfuJHWWBWciBrrmNa0Xy8LJRb/HanZVufPnZcKj14HOGzc3uZ7OFK1rnesXHgM6B9gswiZlHePscqAmZX6mM+uoTNhkou1h0472H48dXYNnbhNDy3SGCynnflFzFK9ZLh2NwiZDde/gzNL8O6nmSDZlqlmT56nWmeWn85Lnvq7r3EdqJnOuOasQzVAIIYQQRk02QyGEEEIYNdkMhRBCCGHUrOwz5GzZamd0Kf/ZJu2OSlCbMdskNe04+zOw/4I7NdgdH6H+Svw9rRc/R+3vbIu9+uqrO2V8r57ePPPlGOKk863gfnJ9prJ0oY7Ovsv2ctfX7Aug9n62navdnMej83PSE9ZZDu7IF60Lp1Dga+D5uTDk0SrqQ+HS93Pf6Hjlzzr22P/g4osv7pRdcskl82tt/6KjSnTN2Lt37/zahc3q+GR/BPVl0jDrRWh/uXE9q/e6T9Ce4Xy83DvVL4P9Fd06or5+HN6uPj233377lnXR/nNpSlgmzv9L02UwnCIFeOH4Z3jN0rDx7V5n3REwLkSe1zddl5xvGPcTnzAPdP1x+PdU/YJ4XurvFN/rjnVR3x/nr8d9pPPZoeO4D9EMhRBCCGHUZDMUQgghhFGzspnMqW9dlltVd7IqVM0ubH7QUD1Wj6vajVXzfKr1ueee27mPTzpXFZwLx+bPau5g1aSqW13oOatstW+HPhF7mVBSF77rTnl3qlBWqWpfs2mRxwOH8gL+5G3XPpaJmtdYBe1C6xV3+vts/A8Rvjt7ppqcePzomHQmBb5X5weHwPI8ArpzTkPrea7yWuBSKqgZkk0r2o9cT20by8ylk3CmfDXlbLdpxcmSzSC6JrKZTE1VfK/OTQ6XvuOOOzplizK0q9mKTd9sQgW6Id6aVZ7l5+aUyoTRdYjrqd8b0nS9FSw/t/YoPHZdigiXvkV/h9nkxb+1nF4E6JrQ1EzGZSpLXsfZbKr10rHD49bN2XWYqaMZCiGEEMKoyWYohBBCCKNmKTNZa22ujlIVrVNB872quuMyVfmxylvVhmyi2bVrV6eMVXJXXXXV/PqKK65YeJ/WmdWWmpmYzTCqVmaznzuwUyME+Dmq3l/FM34ZFh0+eagyZ+rRfmF5uYMHndmDTVp6GCSr7FXdzXJ2ZiE1L7kMuPwcfSabBVQlPHvmkCp5p/5XEwn3qarOuV06ltlMplE5PMfVvMYmGTd+eI6pXPoeFuzMuLp+8RzTevEccBGvQ+CidrUuLFuVCZstVJYsI11ruO/d3OF6aXQvv0/rxZGIfMg10D9yVcdH3+hRbY9z71gH+nyeb9oGnsM6Z7mv3W+Fmkq5D3X95DKXJd8daMzjSE3nl19++fxas1/z+NA6c/vcXE8G6hBCCCGEwySboRBCCCGMmmyGQgghhDBq1nZqPWcq1fA4tgNqmTtlmrNHa1i8+vEw7KfB/kQa7sc2W818zJ/vueeeThn7vKjd3vk68L1a5vwZtusU5a3qojZcrpv6mLgssWxfVps095mOAYbHGF8DXf8QDed2aQu4Li48U+32XE/1g3Ah6UOF1rM/n/oMaX8wXA9tI8vTnQCvtnyeH9pOtvNzf+scu/nmm+fXnDZBn+/8o9RniO9VHwNuq84/l6l+u0Lqt0L7ltun8uJ5rH3mfFe4r10mcB4PevI4p1rga6CbCVmzlXPfqu8Pr8/qT8T11DnNPi8qy6GyiM/6Sd/H/ae+tOyzqmsw97X+bvF6pvJyfcbP5PXioosu6tzHWaY17J7XwZe+9KWdMvYZ0vWS/YKWkQn3n/72rCLLaIZCCCGEMGqyGQohhBDCqFnKTFZVC00hrHp16mlVX7GKTFWFbCZTdR2rAJ26mNV/qtJmFbCGe3PooYbBc/ucWs+plfV7qgrdTpzJRstY/moCctmpuX9V5c2qXZe126Uf4PB5PYSQ66lZU7leTqWu457DQXXcsiw1s+xQtNbmJg4XOq1tdNlsuY91HLB8NeM3j21tP5ex2v7WW2/t3MfZjtV049YM/qwHs7qs9Q53YPTQofXa726uuuz3bKZ2a53OjxtvvHF+feedd3bKeAxwWDynMwG6a7fKxGWmX/QuoDuOXToALeM+cgcYr5OZLHTu8Wc1UXNf6Drr0thwm7TtLn3EovVT5wmbyVSWbg3mcHpdu3kO6W+Dm1/rnnvRDIUQQghh1GQzFEIIIYRRk81QCCGEEEbN0qH1i47jYJ8KtUeyHdCFqatfBtsk1TfAhckyXKa+B+zPoKegc+i+C4N371YbcV+brTJkODbgw4PdKecaIskycadMq62efRbUb4DDqrmP1C+GUy+ojf3BBx+cX/OJzEA3RF9lyc9R+zi/36UD0DEws8GvW5Y7duyY19GFf7ujctQ/gNvs/BbUzs+ozxD7h7FfkPqj8H3a97wW6LrAY1J9E9w457mp45PLtvt0el0TWX7aPuebyTLSVAW7d++eX992222dMvbl0lPr+f0cOv2yl72scx+HaqssuV7qe8Zt1XQc/FnHGK897hgdZSjZzsaPS9ng/PW0zixbnbPse6S/W/w9lQOnneExp7+L3Nf6m8V+QppeoW/friNEfqvn9CGaoRBCCCGMmmyGQgghhDBqljaTzdR5qr5itZSq6Vldp6rQffv2za9VFcoqTs1ArWq+Rd9zmY9ZBajqQHdy76LnK85MpmpRlzV1Vrb2UMKp6lKfy3Vzql3NGMsqWjU1uBBzF97OZhAXzs7jSk1hPMbURMBt1THFn9UUxm11stQw5kVhtofLc889N3+XmhBZvlrXRVmE9bOq413GZjfn2EzOpjGVC8tXT7l2J7BriPIidN46s64zoQ1lWlE5bfU+lx1eTZfcPs0cvGfPnvn13r17O2XsLqDz4+Uvf/n8+tWvfvX8mk1mQHesuPGnc4XniJbxuNLfFG6rO9V9qIzTzLPPPjtPIaHrpcugzPV2aSD0RAf+nVQTIY8P7ReWLfe1/i7y9/REh4svvnh+7bLyu/QULh2AK9NnJgN1CCGEEMKSZDMUQgghhFGTzTVuQAAAAA2hSURBVFAIIYQQRs3SPkOLwgTZRqdlbBNUHwL2JVG7MNsr9TgODt3re3q6+sbwuzUMkXGhtlrGdmC123O9tI4uFf1QKf8XpUnoG/btwkHVZsv9pCfOs5+ClrHPAvuDaH/xfepDwJ+1zi5M24WWOzt+X/v4OtmxY8fcr8n586lsWWY6Jl3or0ttwe/XU7XZn4vXAvVh4KN4zjnnnE6ZC9VeVA/9rHOTP6tvArdvu0Lr+4RjuzBkHXc8B3QN5vY6vyD1A+EQ+te+9rULn8Frq/pmuaNbuJ7uyA3nH7XMsStD+BAdd9xxc7+eZcYOy1bnBq9Fui7xvVrGc1bXAe5rXoN1jLFPko6H888/f+G7WSaun3Xcspz1N9r5f62yzkYzFEIIIYRRk81QCCGEEEbN0mayPqg60IUouwykbCbTzNWsKld13aLQW1U3slpW381lqlZmk48+k9WBqtrluqgq2YXdz1SAQ5pYGHdKN6sjVTXJqlENI+V3uOy4fCoy0DWzcEiwhnw6UxCHn/K4AbomGJUJt8eF6GofuZO4V8mMuizOtOJOQV8mqzyPEZdKQ00fLGvuU10XOLutjiU3p13fuzr3/d6mcSeWcz3VfM9rmM4/Xi81NQHPAU1xwKHUPK90HLmM/bzu6pxmtwknE+dioLh1aGj0fc4dwY1VnpcuA7XK0rkx8Pjg/lTXAZ6nOi95/XSmcx2b7tR6/j3V32F+psss3pdohkIIIYQwarIZCiGEEMKoyWYohBBCCKNmaZ+hRXZOlyaebYJ6mi3fqyn52V9Ew+o4/E9txmyTdMdqsD+Dvpvt1Wr/5D5Q+6ezY7IN1/mgLDqpet0nnc9Q+zHLUv2euL0aJsu+TtoG9tt56UtfuvB9muKdfYjYv0D9Eti2racpsz+DPp9t4CoTlybeHT/DctKyI8n/RNG6cZudf4r2Dc8dnR/cx7wWuNPnNUzX9S+z6Fibrerl2ETaiz5z3b1bfS9YRu5oFfXdcqlCOE2Cm5sutJ6/547j0P5wKRV4nXW/Rc6vZZ3M6u7kpWXcT4t8SLf6nvMbc79bPD7YJ0/HA89L52uqz+fPOo54funvPI8JHTu8LujvzTLze0Y0QyGEEEIYNdkMhRBCCGHUrBxa78L9nLpRT9llNZyq5Nh0pao1FyJ58ODBLa9dSLeq+vn5rq2qKmS1rIY28metC6sAtyP8mnHyciYQhftJ1dNsHtExcNVVV82vOeMt0JWf6yN3MjZ/dv2uqlUucypuVeFzXZYJ+10XTv3vTpjX73F/9M3yDnTnlb6PZc/9pjLjuqh63K01LuTapb1w4d99Q6CHwMlS10SXYdiFY7vT7h988MH5tZ5oz24Mzj2AzRn6bq6zjhX+PVAzKq9Dup64MOtNzMdFOLMSo3PDZdh2c53HgL6P5e7kxWU6L1kmau5yKXS4zpqKg+vpTIIJrQ8hhBBCOEyyGQohhBDCqMlmKIQQQgijZmmfoZmdztnR1e7Mdj+12bK9V0OinS2bP+sz1R9A6z5DfZQYto2qzZbf504B1+9xvdSmeiSFXLvQTbUTMy58l0NhNSyWP6v9n4/L6HucgrOHq88TP2eZ0NpFY0xxKRS2CxcG79rvfLS4XW7+9bXlqzxZTjrH+Pnql7foPq2nC7t3c9GdFL9OFq2z7lgD/qzpCLgv3HEZrn3q63HLLbfMr/fv3z+/5tBs/czzGfBpL7jt2g/uqJy+R6ssc4r84dBHls5PTX8rXFnfIyr0ffxbxd/T8cBrAvtzAt0xoGOFv+d++5xPqvM3W0famWiGQgghhDBqshkKIYQQwqhZW2i9U2myKlvVlqwWc6dt6/tYlefUnfw9VR2z2k1Vuy4LMz9T1Yh9T1PWZ5599tlbPmOre4fGydKFgzpTFX9WNSlnCVfV7qLQWxdmqaY8/qxlXC9nOlFZOhW3M6NuAm6j9i/PPz01nPvKqfgdi7Kpb1XG8Nx0J2dzRlx9/jIh8i70nNu63XPRhdY7meh6xmgb1AWBYTPWrl27OmV8ajmPI32+SzfCz1C3BTafu2z+6o7A/eJMgtt9ar0zv2obeM5quDmnnNH1zJl7nXmZfxv5Wucer9333HNPp0zXD4ZlouOD3+f2ADrXebysY15GMxRCCCGEUZPNUAghhBBGTTZDIYQQQhg1K/sMKS48ju3C6rfjYDuq2pOdb84ivwR9N9sj1U+AbZDOHql2Z7apsi+MPkf9Zo7UNPFKXzu787NSfyLuC/WDcP5Li57v/Cyc74+2re/p0K5PjoSUCTy2tf18pILz7+FjGYBu6KyOc36OOzLCzSuXboFDsHltAfzxPu64GMeROjfd+qK+JOwjor4X3BcujYhLs8F+LRpWzc/Q5/fF+X/pOHLpP44kWTKufTxHga7c9beW+0LLWEbuxHn3m8mozxA/U33WWO7aVvaXveCCCzplPPe1zs4XchWiGQohhBDCqMlmKIQQQgijZhAzmYbjMWpSYPWZqulZRbeM+WSRaWKZdABcLxf2qGpXZ9rrexr9dod8LoMzD7kTk5162pm/VjGTubK+9wGrm7+4nkeCLLkOanJidbaq49XcwfA4V1W6yyLbN/UEq/RdmgYtY/OaZrRn84nOae4HNaEdOHBgfq0mpu2G+8+F1qupyKWe4HVJ+8WlMOHvcb9r2LZLTeBk2deErWXcL26MbXpuskz27dvXKXNZwRk1R3HbVQ48H7Rf2PTGMtLxwHLWjOEuIzTPKR1HbOrWuedOmuC1bB2yjGYohBBCCKMmm6EQQgghjJpshkIIIYQwatbmM8S4NOBqS2Q7sZaxbXEZf5FFxyG4tO19/Y60zN3nwhePhJDrPizjY8M4O77SV7buqIW+z1OZuPD5RfU4FEeybDW1P6N9wzZ59U1gHwOVNfseqSzYB4VTKrjxor4kXOaOb9E6sw+f+mGwb4LWZdN+QotwfjRunVXcOtt3DXb0nQ8ufYXz/VnXvN0k6t/GuDVLx8Cq8uJ+cr4/7ugrRtM+sF+hrkE8F7U9vAY5n7J1EM1QCCGEEEZNNkMhhBBCGDW1jBqxqu4HcOdw1QmHYFdr7Zx1PCiy3DhrkyUQeR4BZG4eO0SWxxa95LnUZiiEEEII4VgjZrIQQgghjJpshkIIIYQwao6azVBVnV9VH62qr1XVTVX1e1V1RVVdVlVfGeid766q+6vqhum/7xniPWNjE7Kcvve/m77vxqr6yFDvGRsbmps/Q/Py5qp6eIj3jI0NyfLSqvpMVf1pVX25qv7WEO8ZIxuS566q+oOpLK+vqouHeM+6GSTP0LqpSbKDfw3gw621d07/di2A8wDsGfj1H2utff/A7xgNm5JlVb0CwI8A+GuttYeq6tyh3jUmNiXP1tp7qQ7/E4DXDPWusbDBdfbHAPxWa+3/rqqrAPwegMsGfN8o2KA8/08A/09r7cNV9RYAPwnguwZ831o4WjRD3wTg6dbaL8/+0Fq7obX2Ob5putv9XFV9afrvr07/fkFVfXb6v8ivVNU3VtVxVfXr089/XlXvRdgONiXL7wXwi621h6bv3LfFPWF5joS5+R0AfnPtLRsfm5JlAzDLbHkGgLsHat/Y2JQ8rwLwB9PrzwD41oHat1aOCs0QgKsB/EmP+/YB+ButtSemmoDfBPA6AN8J4FOttX9WVccBeBGAawFc1Fq7GgCq6swFz/y7VfVGADcDeG9rbWhN1LHOpmR5xbTs8wCOA/DjrbX/cNitCZucm6iqXQBeCuDTh9eMgM3J8scB/P5Uw3cqgLcedksCsDl5/hmAvwvg5wC8HcDpVbWztbb/sFs0IEfLZqgvJwD4hakq8FlMfwAB/DGAX62qEwB8srV2Q1XdBuDyqvoXAP49gN/f4nm/C+A3W2tPVtX3AfgwgLcM3ooArF+WxwN4BYA3A7gYwOeq6urWWnxNtod1y3PGOwF8vLW2+LyGsG7WLcvvAPDrrbWfrqr/CsBvTOfmes9bCItYtzx/cPq8dwP4LIC9AJ7Z4r4jiqPFTHYjgL/c4773ArgPwDWY7GxPBIDW2mcBvBETofxGVb1rai65BsD1AK4D8CF9WGttf2ttdmjRB3vWIXg2IksAdwH4N621p1trtwP4L5hsjsLhsSl5zngnYiJbF5uS5XcD+K3pM74A4GQAZx9OQwKAzf1u3t1a+7bW2msA/OPp3w4cdmsG5mjZDH0awElV9b2zP1TV66vqTXLfGQDumf6P4rswMYfMVOn7WmsfBPArAF5bVWcD2NFa+wSAfwLgtfrSqrqAPr4NwFfX2KaxshFZAvgkJjZ0TO+/AsBta23ZONmUPFFVVwI4C8AX1tymsbIpWe4G8M3TZ7wKk83Q/Wtt2TjZ1O/m2VU121v8CIBfXXO7BuGoMJO11lpVvR3Az1bVDwN4AsAdAH5Abv0lAJ+oqm/HxHHrsenf3wzgfVX1NIBHAbwLwEUAfk2Eprynqt6GiYrvQQDvXlebxsoGZfkpAH+zqm7CRBX8viPdhn00sEF5AhPzykdb0uivhQ3K8n8B8MGpM24D8O7I9PDZoDzfDOAnq6phYia7bl1tGpIcxxFCCCGEUXO0mMlCCCGEEAYhm6EQQgghjJpshkIIIYQwarIZCiGEEMKoyWYohBBCCKMmm6EQQgghjJpshkIIIYQwarIZCiGEEMKo+f8B2DGD4F+RWYgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 10 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plot weights vs the pixel position\n",
    "coef = clf.coef_.copy()\n",
    "plt.figure(figsize=(10, 5))\n",
    "scale = np.abs(coef).max()\n",
    "for i in range(10):\n",
    "    l2_plot = plt.subplot(2, 5, i + 1)\n",
    "    l2_plot.imshow(coef[i].reshape(28, 28), interpolation='nearest',\n",
    "                   cmap=plt.cm.Greys, vmin=-scale, vmax=scale)\n",
    "    l2_plot.set_xticks(())\n",
    "    l2_plot.set_yticks(())\n",
    "    l2_plot.set_xlabel('Class %i' % i)\n",
    "plt.suptitle('classification weights vector $w_j$ for digit class $j$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Exercises: ###  \n",
    "<ul>\n",
    "<li> The snippet above invoked the \"sag\" solver which only deals with $L2$ regularization. Try another solver, e.g. 'liblinear', in scikit that supports $L1$ regularization. You can find more details on <a href=\"http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html\"> http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html.</a>\n",
    "<li> Visualize the misclassified samples. Analyze on how logistic regression performs on MNIST, and note this since we will study the MNIST problem later on using deep neural nets.\n",
    "\n",
    "<li> Try to vary the size of the training set and see how that affects the accuracy.\n",
    "<li>  [Theory] In this notebook, we considered a $K=10$ multi-class logistic regression (MLR) also known as the SoftMax regression. Recall that the our model reads:\n",
    "$$\n",
    "p(y=i|\\boldsymbol{x};W) = \\frac{e^{\\boldsymbol{w}_i^T \\boldsymbol{x}}}{\\sum_{j=0}^9 e^{\\boldsymbol{w}_j^T \\boldsymbol{x}}}\n",
    "$$\n",
    "Use this to construct the likelihood function $P(y|\\boldsymbol{x})$.\n",
    "\n",
    "<li> [Theory] Now consider we're given a set of i.i.d. data $\\mathcal{D}=\\{(\\boldsymbol{x}_1,y_1),\\cdots, (\\boldsymbol{x}_n,y_n)\\}$. As we argued in the main text, the cross-entropy error used in logistic regression is motivated by maximizing the likelihood (or equivalently, minimizing the negative log-likelihood), find the corresponding error function $E(\\boldsymbol{w}_0,\\cdots,\\boldsymbol{w}_9)$ in MLR using the hypothesis above. \n",
    "\n",
    "<li> [Theory] Using the error function derived above, calculate its gradient $\\frac{\\partial E}{\\partial \\boldsymbol{w}_k}$. Can you find an elegant way of writing the gradient?\n",
    "</ul>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python mlreview",
   "language": "python",
   "name": "mlreview"
  },
  "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": 1
}