NB_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"
     ]
    }
   ],
   "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 = 500 # 50000\n",
    "test_size = 10000 # 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Example run in 72.524 s\n",
      "Sparsity with L2 penalty: 26.79%\n",
      "Test score with L2 penalty: 0.8270\n"
     ]
    }
   ],
   "source": [
    "# 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) # 'sag'\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": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind = 0\n",
    "A1int = []\n",
    "A2pos = []\n",
    "A3pred  = []\n",
    "for i in y_test:\n",
    "    if i != clf.predict([ X_test[ind] ]):\n",
    "        #print(i, clf.predict([ X_test[ind] ])  )\n",
    "        A1int.append(i)\n",
    "        A2pos.append(ind)             \n",
    "        A3pred.append(clf.predict([ X_test[ind] ])[0])\n",
    "    ind += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '8, not 5')"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAFVCAYAAACTlUbFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAQaklEQVR4nO3da4ymZ1kH8Ouew85ut9vd0qZZlkOJPYC1tdQ2IaABDVQTY/iA8YiJfFJiNCEaooZQiWjUtMbDBxNDgtWkAQySAGKixgAGoomEcpTSSg/0QJFtu9vd7mnmfW8/zDQMlz1c9+y8OzPd3y9p6OHa+3n2nXf/88zs/rla7z0A+K65rb4BgO1GMAIkghEgEYwAiWAESAQjQCIYAZKFjf7A1lqLiEMRcWzzbgdg5vZFxCP9Of4Q94aDMVZD8aGz+PEAW+WlEfHws/3HswnGYxERP3zTO2NhYeksjgE4N1ZWTsdnP3drxPN8pXs2wbh6wMJSLCzsPttjALYNv/kCkAhGgEQwAiSCESARjACJYARIBCNAIhgBkrP+A94wE9MX6C6iubbVd0CBJ0aARDACJIIRIBGMAIlgBEgEI0AiGAESwQiQCEaARDACJCqBO92MqnPt2TdLbvj6beRWq9cfnW3FSl51LiL6SMtvYLhX70HNcNN5YgRIBCNAIhgBEsEIkAhGgEQwAiSCESARjACJYARINF/O1gyaJ+XWyei5kxm0SQbObNPp5l8/ImLg2Fk8CrS5+qHlNktEtIXauUPNmxEDjZqRn9csrr/pl96yKwNsU4IRIBGMAIlgBEgEI0AiGAESwQiQCEaARDACJIIRIFEJfAZtMtAxK1YC20h1cKAON5OaX0S0leJrsDKpX39GlcB+wVL93GJ1rZ04XT7yiZsuK8+e/oUnyrO/cuVnSnNvP/Bw+czv+4dfLc9e/XcnyrMjm876/PZ/Htv+dwhwjglGgEQwAiSCESARjACJYARIBCNAIhgBEsEIkAhGgGR7VgKr9bmBLWKzqPlFDFT9RrbpTeo1u6FK4EB9r51Zrg0urwxcf2B2cbE8+ujNB8uzy3tr75m5gVs9dUn9Y7D84IHy7K1Hf7w096GDj5fP/J03fbw8+/7/fHN59sBdx8qzVX1uftPPrPLECJAIRoBEMAIkghEgEYwAiWAESAQjQCIYARLBCJAIRoBke1YCq1W/kc17A4Y2+hUreUM1v+qGvoihSl4bqe+dqm3J66dOlY/sxTMjIk796HXl2elivRparfr1gTba4vH69S98sP4ssutYbfvhAz9xafnMj8zfUJ498TNHy7MX3Vrf1Dh3pvZrYTa/ums8MQIkghEgEYwAiWAESAQjQCIYARLBCJAIRoBEMAIk27P5MgNDS6NGFldNay2VndRmiYjoJ06W5qYnB5ovN76qPHv8UP2t2QZegmpHZf5M/T3woq/VXquIiIUv31uenV7xstLc/75mX/nMU4fqS8ZuftnXy7N37qk3ahaX6y2wreKJESARjACJYARIBCNAIhgBEsEIkAhGgEQwAiSCESARjADJuasEzmBxVZsM1Oz6QM1vYLZcHxy417YyUJkaqASOLK6qVv3mX3KwfOZDr7mwPDuiDXy45k/Xhi/9wrH6oV+4qzzal+pLox75sf21Mxfr75fHnrqgPPuR+24sz77yiRPl2T6//Z/Htv8dApxjghEgEYwAiWAESAQjQCIYARLBCJAIRoBEMAIkghEgOXeVwLnqfrbBqt8WK9cHi9sEI2KovtjPLNdnR7YPXnNFae7wDx4onzmtt+GGLByvzx782H2lucnhx8tnzu2/qDx75Oary7OT3dXB8pFx6u5azTAi4lV3PFGe7fPz9ZtYrGfBVvHECJAIRoBEMAIkghEgEYwAiWAESAQjQCIYARLBCJAIRoDk3FUCZ2BkO9xIzW4mRq4/siVwWp+de9HF5dlvve5Aae70gXq9q4+82wZerguPbH7d8vQbry8f+eTL6z+x6a766zUpVij7nvrP/8AXBp6F2kB1b+DYPnLuFvHECJAIRoBEMAIkghEgEYwAiWAESAQjQCIYARLBCJAIRoDk3FUCpwOb76qVoYFOYBuqIQ300WZRNZyMVALr159eUt9mV636TRfLR0ZfmE0t8+gV9c/vT76itv2wDXwI2kAjceT1mj9dmzv4yfqGvovvPFye7YsDm/8Gfn1VN2v22LrqoCdGgEQwAiSCESARjACJYARIBCNAIhgBEsEIkAhGgGRHL8PaSdqkXo/o04EqxYAnr9pXv4fqp8yBckKbDrQjluvnzq2M3EPxzIHrD7V/Bh5FLv/HI6W5uSdP1K+/NHCzM2IZFsAOJBgBEsEIkAhGgEQwAiSCESARjACJYARIBCNAIhgBkp1dCZwbqBYNLI0aWnBVnR0585ID5dEHf/LS+rkDL1e1Ztdn016M+TP12aHFVQOzVQsDlcSXf+CB+nC1Orc4o1/GI9W9WSyF20KeGAESwQiQCEaARDACJIIRIBGMAIlgBEgEI0AiGAESwQiQnLtK4Kzqe1Wz2kxWvNfp/r3lIx9488UD16+PjmzTq26zq27di4iYP12fnRuoBA6pvg0G3i4ru+uzj73hZeXZSz79YP3gHaQV64N95IOwyTwxAiSCESARjACJYARIBCNAIhgBEsEIkAhGgEQwAiSCESDZnlsCi/XBkQ11LWazJbDvXizNPfBT9ZrfZKl+/flT9drUdOCjvXCqNrfvwfoHYeFkfXbxeL2/uHTPt8uzp155sDR3+Nql8pkjjxeTpR20eW+rrz9SDR6pHFeO29TTAF4ABCNAIhgBEsEIkAhGgEQwAiSCESARjACJYARIBCNAsj0rgUXVbWMRMVRvapP67KOvr1X9Tl02Up2r15tG6oN7Hy6PxqVfPlmam//SN8pn9jP11X99uV4JXJlOyrPTaw/VBmf0yLD/3oFVidOBzusszGqzZtUm1/yGLr1lVwbYpgQjQCIYARLBCJAIRoBEMAIkghEgEYwAiWAESHZ082VoWc5Am2XEk1cWz12oX38ysIdpYbneDjj474/VD37sSG1uaVf5yDZf/zzcW70hMrnxuvLs8RfX3vIji8NGLB4pbhmLiJirvV59pCGy1W2WiPqv2/nZ3sZz8cQIkAhGgEQwAiSCESARjACJYARIBCNAIhgBEsEIkAhGgGR7VgKLlaE2m5bfkF5sWPX5gWVc1UMjYulIfbY9frQ8W62O9fmB3takvtxpcsMry7PfueGC8ux0sTg48N5aqO0Ni4iI9tC368O7dxcP3QaVwO1QNdxEnhgBEsEIkAhGgEQwAiSCESARjACJYARIBCNAIhgBEsEIkGzPSuAWW37RQMVsT7HmVm/DxdxKvV514cP1g3sfqCXuqm3/W7780vKZj1+zpzw72V1/DQYalOWPw0jD7cD/LNeHRyqU1a2KxW2CERF94CfWB84dMrLVcIt4YgRIBCNAIhgBEsEIkAhGgEQwAiSCESARjACJYARIBCNAohL4DCZ7Bmpbu2ods/kj9Zd67kz98keuqn9um+y6ojx78tJibWtkQd1KfXZuYHakvjctfmj3PTgpn7nnzgfKs22h/j7oi8XZanUwImJ+4MUaOLYP1PxGaolbxRMjQCIYARLBCJAIRoBEMAIkghEgEYwAiWAESAQjQKL58gyWHjlenr3kPy4uzT1+XX0R1UDvZshTh+qNg7nifqfqXEREG1gIFvWXa2h236O1RstFn763fGZbGPiILdWWjEVERLH5MrK0amjB1UhDZQcsuBrhiREgEYwAiWAESAQjQCIYARLBCJAIRoBEMAIkghEgEYwAyY6uBPaRRUwjS4Cm9dn999Y2Vx29ql4Fm+yud9zmlkeWEJVH6zW7gTreyPUXT9YPvuxfv1m/h+O1umfbs6d8Zuxeql+/uuAqIvp8sWq4MPB8M6MFV0N2QH3QEyNAIhgBEsEIkAhGgEQwAiSCESARjACJYARIBCNAIhgBkh1dCRypFo0snRupWO06/FRp7so7anMREcev3F+e7a229S4i4qmD9W12e4vb9OYm9Vf2wnuOlmfbY0fKs31lpX7u3r21wYFtfkM1v4HZ6vuwL85or+TIr6+RjYI7gCdGgEQwAiSCESARjACJYARIBCNAIhgBEsEIkAhGgEQwAiTbsxJYrCL1qNeQWkzLs30Wny/ql48L7x6ozi3X63AXfbFeHxw5dyYGNu+1trs8W67kLdRrdn2gQlre/BcRUdxsObT9cVY1vx2w+W+EJ0aARDACJIIRIBGMAIlgBEgEI0AiGAESwQiQCEaARDACJNuzEjgDfX7gc0Ab2ClYrEK1gW16feTT1UAdLaYDvcRJfUveLIzV0QYOLp7b5wYOHbn+wLnlqt95vM1vVjwxAiSCESARjACJYARIBCNAIhgBEsEIkAhGgEQwAiTnTfNlyEiToLyQq946aQPbjfrAbqVoAwueqqOz+tQ6q4ZGtfliwdR5zRMjQCIYARLBCJAIRoBEMAIkghEgEYwAiWAESAQjQCIYARKVwHNkZBmXhUURrQ8sJJuBmX0M1Px2BE+MAIlgBEgEI0AiGAESwQiQCEaARDACJIIRIBGMAMlZN19WVk5vxn2w3nRrWx/bgeYLs1DNq9Y3+AZsrb0kIh7a0A8G2Fov7b0//Gz/8WyCsUXEoYg4tsEbA9gK+yLikf4c4bfhYAR4ofKbLwCJYARIBCNAIhgBEsEIkAhGhrTWFlprf9Bau6+1drK1dm9r7ZbW2jl/L7XW7m+tvaMw96nWWk9/ffAc3CI7lJ0vjPrtiHh7RPxyRHw1Im6KiL+JiKMR8RdbeF/P530Rccu6fz65VTfC9ueJkVGvjYiP9t4/0Xu/v/f+4Yj4l1gNyJLW2ivWntre0lr7ZGvtRGvti62116a5n26tfbW1dnrt6fC31v23T0XE5RHxZ08/BT7PZU/03h9d99fR8s+Y845gZNRnIuKNrbWrIyJaa9dHxI9ExD9t4Kw/jIjbIuLVEXF3RHygtbawdu6NEfH3EfHBiLguIt4TEe9trb1t7ce+JVYrqbdExIvX/noub22tHV4L2ttaa/s2cL+cJ3wpzag/iYj9EXFXa20SEfMR8a7e+wc2cNZtvfdPRES01n4vVr80vzIi7oqI34yIf+u9v3dt9u7W2jUR8c6IuL33/vja9Y/13h99nuvcERH3RcSjEXFtRPxRRFwfETdv4J45DwhGRv1cRPxSRPxirAbZqyPiz1trj/Te/3bwrC+t+/tvrf3vZbEajN8fER9N85+NiHe01uZ775PqRXrv71v3j19prd0TEZ9rrf1Q7/3zg/fMeUAwMurWiPjj3vvTv6v75dba5RHxuxExGozL6/7+6e8RPv3tnbbu38W6f7cZPr927avW/h6+h+8xMuqCiJimfzeJzX8v/Xesfu9yvddFxN3rnhbPxOqX8qN+ICIW47tPqfA9PDEy6uMR8a7W2jdj9UvpG2L1+4Hv3+Tr/GlE/Fdr7d0R8aFY/d3wX4+IX1s3c39EvH7tzySe7r0fzoe01q6IiLfG6m8OHY6Ia9bOvjNWvzSH/8cTI6N+IyI+HBF/FRFfi9XfVf7riHj30wOttfe01u4/m4usfe/vZyPi5yPiKxHx+xFxS+/99nVjt0TEKyLiGxHxnWc56kxEvDEi/jkivh4Rfxmrf7zoTSPfp+T84v+PkU3XWrs9IqL3/ratvRPYGF9KMwtviIjXb/VNwEZ5YgRIfI8RIBGMAIlgBEgEI0AiGAESwQiQCEaARDACJIIRIPk/O4Nk9NkJ+JkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "check = 11\n",
    "\n",
    "plt.figure(dpi=100)\n",
    "misplot = plt.subplot()\n",
    "plt.imshow(X_test[ A2pos[check] ].reshape(28, 28), interpolation='nearest') \n",
    "misplot.set_xticks(())\n",
    "misplot.set_yticks(())\n",
    "misplot.set_xlabel(A1int[check] + \", not \" + A3pred[check])\n",
    "#print( A1int[check] + \", not \" + A3pred[check]  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAFCCAYAAAAe+Ly1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABNsElEQVR4nO2debAkV3Xmv6N9X2hJvaj3TS2phbqFZAkZAxayh7HD2GwDBhszgT1BjFg9YOFtAs/YAZoZj42NNXaADQwx7NiMsRmjwJIMGAkJoZbU2npf1Uit1tralzt/VL3ky0+vTtd7r+plvc7vF9ERWS+zMm/ec++t2+c7594opcAYY4wxpi0c0nQBjDHGGGOmE09+jDHGGNMqPPkxxhhjTKvw5McYY4wxrcKTH2OMMca0Ck9+jDHGGNMqPPkxxhhjTKvw5McMnIjYFhGXDunet0fEK+nzGRFxc0Q8GhHv0fPDeu4oMZGyDdM2bUHb3BDu/+mI+EP6PBH7DqSdTnc7iYg7I+Lc6XqeMYc1XQBjJkIp5Wz5028BuLaUsrb7+c+m+oyI2Abg10sp30qeOzIMsmzjvXuTjFp5umibGyoTsS9fO6J1Ny6llDObLoNpF/b8mJnOIgC3N10IM3pExLD+czepNjfE8hhjJognP2bSRMSCiPjbiNgbEfsi4uPjXPOhiNjclQjuiIjXyvnLI2J39/zdEfGqA/y9csdHxNUAfhrAxyNif0SslPM9y9erXBHxWQALAXy9e8/fGue5Z0bEtRHxUFdmeI2807aI+EBE3BoRD0fEFyPiqHHq5t9HxNfp86aI+BJ93hkRa7rH8yLiq9132cpyi5TtPJJkvtx99h+izprxyjbeu/eywzg2/or87WMR8Wf0OSv/uHbqUZ5+6v7yiLgVwGM64ejeZ0H3+J0RUSJidvfzByPik/p+8v3x2lzPMh2oPN1r1kbED7t1/EUAR8n5vu07dm2vdiz3PWD/7V430D48zv3fGhHf7FXnxgyFUor/+d+E/wE4FMAtAP4EwLHoDNgv657bBuDS7vEbAcxDZ6L9JgCPAZjbPXcGgJ0A5nU/LwawrNff9d7dz9ei49oHn8/K10e5as+Q+x4OYBOA3wFwBIBLADwK4Ay59obu/V8E4E4A7xynDpcCeKhbhrkAtgPYTece7J47BMBNAP5z95lLAWwB8G+kbEd07/HebjlfB+BpAH/Yb9nEdj3tIO+xCMDjAE6gtrEHwEXdzz3L34eduDz91v06AAsAHD1OWXcCOAtAALgNwEYAq7qfNwE4t4+2fy26be5AZeqjPGM2e3/3Xm8A8Mw4NpuIfS/V44n032H34XHK8hEAf9z0mOZ/7fpnz4+ZLD+BzoD4wVLKY6WUJ0sp39WLSilfLqXcU0p5vpTyRXR+bH6ie/o5AEcCOCsiDi+lbCulbE7+PrDyHaBcGRcBOA7AR0spT5dSrgbwDwB+Wa77s+79HwDwdQBr9EallC3o/FCuAfAKAN8EsDsiVnU/f6eU8jyACwCcWkr5L91nbgHwCQBvHqdsh3Wf/Uwp5W/RmegoByxbl77sUErZDuCHAH6p+6dLADxeSrm++zkrf1/tiN6v37rfWUp5Ypx7PNS9x8+iM2lZD+AkAK9GZ+J5CwBExK9HxOoe5ZhombLyXITOROZPuzb7CoAbk2f1Y99+6Lvep6EPr0ZnImrMtOHJj5ksCwBsL6U8m10UEW+LiHVdSeAhdAa6UwCglLIJwPsAfBjAfRHxhYiY1+vvgyxfVq4DMA/Azu6kZIztAE6X635Ex4+j8wM5Hv8C4JUAXt49vhadic8rup+Bjmdl3lhZu+X9HQCzxynb7lJKob/tHOeZfZVtgnb4HH78g/+W7ucxsvL31Y669Fv3473zGA+i877vA/AxAI8AOBnAfwQFy5dSPllKWT+gMmXlGc9m2ydwbXbvjL7rfRr68Gp0JqHGTBue/JjJshPAwvFiGMaIiEXo/A//XQBmlVJOQmeQi7FrSimfK6W8DJ0fyALgiuzvgyhfH+Uq+h3iHgALIoL7zkIAuydYvjHGJj8/1T3+F7xw8rMTwNZSykn07/hSys/JvfYAOD0igv62YILlqb37BOzwZQCvjIj5AF6L+uQnK/+B2hGXp9+6z+z3EIDz0ZFtrkVn8nMuOj/AXxu7KCKuTe7B9FOmrDzj2WzhBK7N7Js994D9Fxh+H46I47rvcEdWDmMGjSc/ZrLcgM5g/NGIODYijoqIn5RrjkVn0NsLdAJ80fmRQffzGRFxSUQcCeBJAE8AeK7X3wdYvrRcAO5FJy5lPL6PTszDb0XE4dFZU+UXAHxhguUb41/QCaA9upSyC8B30JFgZgG4md7lkW4A6dERcWhErI6IC+Re16FTT++KiMMi4hfRn5THVO8+ETuUUvai47X6FDoTnTvpdFb+A7UjtsUg6v5BdOJrxrw8j6ATQ/NXpZTnuu99PDpyZD9MtUzXAXgWwHu6NnsdettsovbN2nE//RcYfh8+G5328njyHsYMHE9+zKTo/lD8AoDlAHYA2IVOMCRfcweAP0Zn0L4XwDkA/pUuORLARwHcj44Ucxo6ckivvw+kfH2U6yMAfq/r5v+A3PdpAK8B8G+75bsSwNtKKXdNpHx0vw0A9qMz6UEp5RF0goH/dezHmN5lDYCt3ed+EsCJ45TtdQDegY6H41fQiT95agJFqt4dnfqaiB0+h05gLnt90vL30Y64PO/B1Ov+QXTiZsbK+Ag6MT+foGv6lmGm2h7IZm/vlu1NAP72ANf2a9+sHR+w/3avG3YftuRlGiHq8rEx5mAiIr4P4C9LKZ9quiwzhYj4DwD2l1I+d8CLG2am2zciPgbggVLKHzRdFtMu7Pkx5iAiIl4REXO6ssivAXgxgH9qulwzjHMwotlHB5N9u/E+Pw/gmqbLYtqHVxw15uDiDABfQiejaTOAN5RS9jRbpBnHOQAmJWNOAweFfbuxUf8HwGfQlXyNmU4sexljTJeI+BqAdaWUDzdcFGPMEPHkxxhjjDGtwjE/xhhjjGkVnvwYY4wxplV48mOMMcaYVuHJjzHGGGNahSc/xhhjjGkVnvwYY4wxplV48mOMMcaYVuHJjzHGGGNahSc/xhhjjGkVnvwYY4wxplV48mOMMcaYVuHJjzHGGGNahSc/xhhjjGkVnvwYY4wxplV48mOMMcaYVuHJjzHGGGNahSc/xhhjjGkVnvwYY4wxplV48mOMMcaYVuHJjzHGGGNahSc/xhhjjGkVnvwYY4wxplV48mOMMcaYVuHJjzHGGGNahSc/xhhjjGkVnvwYY4wxplV48mOMMcaYVuHJjzHGGGNahSc/xhhjjGkVnvwYY4wxplV48mOMMcaYVuHJjzHGGGNahSc/xhhjjGkVh03k4uOOO67MmjVrWGUxB2Dfvn3Yv39/DOJetmWzDNKWgO3ZNO6bBw+25cHFjh077i+lnKp/n9DkZ9asWbj88ssHVyozIa644oqB3cu2bJZB2hLo2PNDH/rQQO9p+uejH/3owO7lvtksHmcPLi677LLt4/3dspcxxhhjWsWEPD8zncMPP7znuWeffbbv+5RSBlEcMwUienuls3PK888/P4jimCmS9c2J9LeJ9GMzHA47rPfPynPPPdf3fTzONk/WLyfCM888M5D7DBJ7fowxxhjTKjz5McYYY0yr8OTHGGOMMa1ixsb8cFzHkUceWTvHmjNfp1o0a8pPPPFE7Rxr0xpDcsQRR1THqmU+/fTTByx729E4G67fQw89tHaOP/N1Gg/A9tK4gkMO+fEcX+/PmrZ+z/Ej/aH9iu2k/erOO++sjk855ZRxv6P3nD17du3cySefXB0fddRRtXNPPvlkdfzoo4/Wzj311FPjv4CpyOLlsnE2g/uqjpf8PH0291v93ijGkIwaE7Elj4O9xlxFf+uyeC5+nn6vqd9Me36MMcYY0yo8+THGGGNMq2hc9mK32tFHH107x+43lTlYejrmmGNq5x544IFxnzVnzpzaZ3arKvw8dqUDwGOPPVYdqzRy3HHH9bwnu2rVbcvPm6kpnlxutg+Qy4hsZ3XH9gvbQdvDSSedVB3ff//9tXPsup8/f37t3EMPPVQd79+/v+ezVUpjHn/88Z7nRh1+Z65DoF43as9du3ZVx9qnebVb7kcqobBctm3btto5lrq0vcydO7c6/tGPflQ7x+1CU3gXL15cHWd9c6bKLTwuffe7362dO+GEE6pjrc8TTzyxOt6+vb5W3L333lsd89h67LHH1q5jm2jfZFtmMko2nqjcwjbStnkw2DKT/7MQAratysXc//odg7XeuSzctwHgvvvuq45XrFhRO8dhEHfccUft3JYtW6rjtWvX1s5xH56MdGbPjzHGGGNahSc/xhhjjGkVnvwYY4wxplVMe8yPxtmwbscprEBds1T9nuNwdu7c2fOefJ3G57COvHHjxto51j1Vw+a4BdWiWXvU1NoXvehF45YLyOOBRhWN02ANWO3MMSOa6n7aaadVxxxHAAD33HNPdcyasmrWbAeOYQDq+jPHo2iZVTfm9qj69o4dO6rjDRs21M69+tWvxkxkz549tc8cA6XxA3v37q2Ot27dWjvH8SWabs4xQBybp22ev6dticuyaNGi2rndu3dXx9o3Tz31xxs7a+wJ9019HscrzJS+qe/AsRQ69nB8FfdFAFi/fn3Pe3L9cnvQ2KC77767Ol64cGHtHMdaqb2yZUp4TNbxhPvxwbAUSbakAI+PQH2cXbBgQe0c/6Zpf+Zn8LH2X439YzieUsvMY7K2D46L5D4K1McS/f3mNsfLZvSLPT/GGGOMaRWe/BhjjDGmVQxM9mI3l7qnWIZS1ym7wTVtj1NtNWWR3Z6PPPJI7Ry709ml+/DDD9euYzegpjKzBKCuPi6nutsyaYZd7ZoCzK7offv2oUn4HbJVVxV+J313vo/KUtw+1IWt9+l1HZdLv8OyVyZtaYonP0NlL24T55xzzrhlHBXYxa3vcdddd1XH6qpmCUn7JvcrrTfuOyoZs4SULU8wb9686lhtzbKzyjR8H5Wd+d1VDuGyaPvh9qkrVk83WZozj4MqIfE4pX2AbaKyM4cUcL0D9ZW3ua51bOP28eCDD9bO8fe0H3Gavco0/O76PiydqL3Y7vqu0w2PWdrGuc50DGY7axvne+pyDnxP7eu9lhzIbKn9nu+xefPm2jle2V1/T/k+KrNlSxro2DJR7PkxxhhjTKvw5McYY4wxrcKTH2OMMca0ioHF/LBmrholn1Odjq/lJfOBum6tOiFrg1n6KWuWql+y/qzbELC+qDo1v4PGKmhcEcMattaDpp82CdtEtWj+rDo8a756jm2kMWGsy6t+32updbV5pv9yO9LlFLIl2bkN6PtwarRq6/y8UbBrtrXI6tWrq2PV6TX1neH+qPfkdqD9ivvg8ccfXx3rtjBsd74OqMeM6fc4HiLbZTrbSkfpd/fy6YDbq8az8DIOGouYjVlsP63P8847rzrmGB+gPtZpjBHDNtH+wO+jW8+w/bi/AXWbaAwhv4OOXzr2NAmXTfsQj28ag8b1wr8pQL1f6m8R169+j+s361/cdrRcHDer6ez8rjrGsy31nvzbq8vdcBs//fTTMVHs+THGGGNMq/DkxxhjjDGtYlr8uSwLqETAbk91X/I5dVdmshe7YDkVUGU1Pqdub3bPano3uw/V3cbudF6BEgA2bdpUHS9btgwzEa4nlZpYUlKXPNskkxy0PnvJiLycAVB346prne+hbYVdruq652tVQuG2mu04Pepw2dWe/B4qPbEspdIC21frgt377KLXcYH7pvZ9do1rmbMUaP6sKbW883i2W/UooZIwv5/2Ab5WxzOWD7R/sF30eSyZcj/SpTpY9lJb8j01lZ7RcYDLpW2Tz6mUrZ9HlWy8zMa6bKkHHme1rtkO/DwNN8nCC/j3VcvF/UvbEbc57Wv8Dro6eL+7z/fCnh9jjDHGtApPfowxxhjTKjz5McYYY0yrmHTMj6YMs4auaXusX+o51hdVi840eo4pUV2S9eEshZXPqW7MuqrG7nA5VVdlLVX1c/6sqYC8g/Z0o+/Odav24s9ZjIzWS6YVc7yApuGyjTiVkrdeAOopmaobc7q1tods+QFuf1oujivIdrdvgve+9721z1deeWV1rDp5tvUF903efRuo90eN1+H7aBvh7Si4LWn98tYDutQEx9xpHADbQre94c861vD3tP1oH5hONCWZ++Y111xTO8ftV9+Bt3zQlGGuQ42T4nFpzpw5tXM6boyhMT8c57Nhw4baOU6xVpvwthuaAs3vp+2PP2s9NBm/pXE2bEsdl/hztpSG/sZwG9fvcf1qTA5fy+1dbaLxfQzHU2mMFvcvTZ/n39dsmxt9V42Fmij2/BhjjDGmVXjyY4wxxphWMWnZS13R2SqomYuS3WqaBskyh6aps5ShO7Kz5MLuPXWdsnSn7jYuC7tf9f66MzCnQKubLkvd5+9N967uKimyWzVzG6sLVF3mDNeZ3jNLTe1lZ21/69atq451dVu2u7rq2e2e7WCfpWDqKsbaVqebj33sY7XPWd+84447quN77rmndo6lLr0H952JLEPBn7n9qHud0XbFz9NVqXstcwHUbagucy6XlrnJndxVPrjqqquqYw094NWYlyxZUjvH/YWlBKAul2nf4WtVwuT+yNKF9k0+p/I+98fMztr/uE+rffhctnLydKPjukqaTCbPZf2ZJWcdh7ieVPZaunRpdcySm45tWUgEy9H6m7JixYrqWH+vefxXSSzbaYHbS7bCeC/s+THGGGNMq/DkxxhjjDGtwpMfY4wxxrSKad+uWLVA1t41doA1RY0F4ZQ+Td/tle6ncQXZztz8vAULFtTO7dy5szrWNHguS7bDu27j0HR6NMO6eJbmq6nE/L0stkvvmaVGc3vpFWMA1NN3zz777No5XoZf4x0y/ZxjvTSugL+XvU+TadJjcHm0D3B8xKpVq2rneAkJjS/h+2hsAev2+j2OreO4Bo0t4JgEbUvcznQncI7zUdvyM3Ss0e1Sej2vabg+NW6JP+t4yW1AlwbhcVBj4theGuPHdueYO40d41gXrUt+Hi+DANTjILXd8jihbYdjVjTua5TGWW7jOu5x29Uys221DfBvk7YB7uvaZ7n9s000TpbtoDFat912W3WscT1ZPCj3RW1jXEf6e8O/tY75McYYY4w5AJ78GGOMMaZVTFr26je1XT+rfMCygLrKOLVN3WHZ7uK93ILLly+vXbdy5crqWF2nmcTC5dJz7IbX9+F31TTHUdptmOta653dnuqKZpd2Jq9kz1OJg13y7NJV1/rdd99dHavsxXXda1Va4IUSCrtcszadLWnQhOylq8Vy/aoMsGjRoupYbcR9R/sHrxisdZq5zXuNG+qG57rX+7Ospjbj+2Tf03GIpRluc0Cehj9stL5e+tKX9ry235V5VfZi+UDbDksue/bsqZ27/vrrq2Puf9o3GbUXp9brWMP9Vs+xHKeyENdDkys6K2rLTKbMxilu41m6uY6lbJdsRwEez1RW43Jqn+VlMxRuY7pSOPdLfR9+B07HB164BM1EsefHGGOMMa3Ckx9jjDHGtApPfowxxhjTKoaS6p7pxnqO9UXV4fl7GqvA2qDGFcyePbs6PuOMM6rjZcuW1a5j7Vs1ctYzNaZo165d1bGmyHLsgGrR/K66S7Yu3d0kXLfZlg9KFgeTpS5z/WosFH+Pj7W+eMdkTZnluBbV3TV9kuFy6vuwbdXOTcaIjEcWg5TFzvG1GifCMReckgzU006zmDFOf81SkHXMyOKosvgVHhc0FoS/p2XOdrIeJbhetJ3zGKlLg3Bch46lHNN4ww031M7x1iLcH7WvsM152Qmg3t91qyBujzpe8vvpeMKxlmq7UUp1z7ZN4barsS38fvq7yEuD6FIrbIezzjqr5/O4Hek2GPzsG2+8sXaO24DGLHEb0y1UeJkSjZPlcUdtx/eZTMzsaI3SxhhjjDFDxpMfY4wxxrSKoche6pZmF7O6rrK0PXbBaqptljbOLtIzzzyzOs7SBzOpQl1q7Irbtm1b7Ry7IbO0vSx1f5Tc7Nk76Cq5/FntzLZUNzWfU3c9f+53xWy2OVBPXVaZi6WXbBVShdu0yiujlF6rqD1ZolI5hM9p/2Bba9/MVqdltzxLxtrmuZwsaer39H3Yva5yHKfn64qwLNWM0rITGZkkrX2T04lVyuC613R2lqK0f6hMOsbcuXNrn7k/6K7dLEXq6v+ZLMpjRra6r57TemmSbOV1bo933nln7Vy2NAiPkSqrcaq49g2Wsbk/q4S5ZcuW6jhbGkTT2TN5k++TjTPaBnRcmCj2/BhjjDGmVXjyY4wxxphWMWnZS92H7OpWGYBdWSpdsCtLJQ/ONFD3G7v71P3Km/9lcguXS93u/H7qBmfXupaLXXiaocCuRZUYuM5GaSNFdUOyC1vlnSzbJFt1NVuZlut+9+7d1THXJVBfvVs3omVbat1mbtVMfuRyartSWWG60ffolcmhnzXTQtsow25svWfmlmdXPGd58OrdCmfrAXW3vGaOcHvRlZo1Y41hOVRlvCaz97JxNst6y76n/S/bKJZtpBmW3AfZzrp6MI/B3IeB+qrD2qe5zCpXZ/1vsuPQsJmILbledMziutBQDu57+n7cv3fs2FE7x2MGS5EKS5H6W84hHzp2sqSp0irXS5aFme0qMRns+THGGGNMq/DkxxhjjDGtwpMfY4wxxrSKgYlorNtpqiNrfKoTcjyCxgoxGkfBKZOazthrRdt+V5sF6lq3lpnvo3ERnE6r6X7ZatZ8Hy1Lk6nvqrPyO2SptnqO20eW1qnaPqdFsi6tNudVnXUn4mxVV44X0NgBRvV6bqvZarpaD7oT8nTA76/xcdyvtE1y2TXugOtbU045ZoZXPAeADRs2VMcc06exJozeg+N8Vq5cWTvHNlS7cLnU1hwrpLFPbHuNI9IVdocNl+W6666rnXvVq15VHWer7Gt/yFL7Od5K7czvzuOXprPzddo3ub9rHArHlek4y/FBOp6wbXUs5RgZbe8TWepiELAtdTzjcULr7JZbbqmOdadzXgZC42c4Zkvfle3Ax3odrzbNv5FAPX1e+wX/hmnbzH7fOH43W5lcxzVNpx8Pe36MMcYY0yo8+THGGGNMq5iWjU2ZbJVQdoMDdReeyiHs+laXF7u72YWt7l5266vbncui5eL7qLTFEpm6K7lcKiOwC3aUVnjO0r9VVmDXbSaJZfdUO2/atKk6ZperpiNze1B3NrvPtW3yO2QSn7rWMxmB3bNNyFwZN910U+0zu6pVWs6WZmApil3TQL0PqGuc+xLLLdpe+J5avyxrcPsA8qUt2PYqM3CZtQ3yqrYXXnghRoVs4+QsPV9tyZ91POP6VBml17IXulJzNgZzmTWVnr+nSxqw/bKlVbKlM6Zb5srQ9sj9UscztpH2G96wVMcetrMu+8BjH/8W6phw/fXXV8cqD7NEqvZiKTRbrVs31uWxW+3FttVxph/s+THGGGNMq/DkxxhjjDGtwpMfY4wxxrSKocT8cGwGkKcasz6sKZf8WeMxOJVOtUe+J+uCWfqdLrvOn1Vz5Xvq7riswWosEuvbqqWO8k7gvVCtnWMsNKaJr9X2wPEj2nZ4mQSOAeNUa6CehqvxQBx7pSnr3CZUU+Zzah9uA3pulHcF17bMsQZab2wnTTfneKAsNXzXrl21c7xDNd/jnHPOqV3HMQm6dAbbXlNaue/rmMFxB1oP/A7aDnSpi1Fh7dq1tc/8TtrHuK61vfI5jVNkO+j32M78PY214rH0xhtvrJ3jGDAdnzlGReOIODYk2xZD28BMGWd5DNEtLLiv6HjGv0fab7LlAbTvj6H1zjE4fAzUx2DdYmjFihXVMW/dAdR/h7Xv8ecsXmwy2PNjjDHGmFbhyY8xxhhjWsXAZK9sZ252T6lrM0s/ZXdftvqupkiy3MT31NQ/dvWpe56v1XJladvs+tPvMeqqzXZonm76dQ2rndmW2W7wmkrMdtCURbYtu/XnzZvX8/7q6s7kAHYbqw34nlnqsNbXKC1VANTfS/sR15W2SZYc9f05FV3Pbdy4sTrWZSJYrliyZMm45QDqbUnLzM9T+ZjlK/0e21rd6+ze37p1a+3cqC1XMIa2ZR5bszFYbcLyhaYac1tWSYyfke0mznWt8grXuz6b76P2ymzJds9CLpqG+6X+hvG7864BQN0OOtbwPbU/8++p2pJlbT5W2YulVg09YFmZZS6g/luo5WKJL1t9XFdXz8bkfrDnxxhjjDGtwpMfY4wxxrQKT36MMcYY0yqmZXsL1vE0NZw1+0xv1vRhjhHg2ASgrn3zOY3p4FTKLBYp29pAl4pnvVnP8TM0RVE/jwpZ/Ipq+6wj67tzLIbq2xwDpFo0p2vy0u26pQLfc/bs2T3LrCnOXO+aOsntRd+V20B2bqq69KBRu7Cmn8V6ab1x7IbGcbC+z0ve67Uca7J8+fLaddxe9P533313zzJn8Xj8Dpq6z2XRWAZNKR8VdJzlusjiGzUWitu5bnugY3IvsnRlXipAY0jYJvo+/OwshlDfJ4s/GtVxVvsXbw+hsS76G9oLjsEB6rE2GivH8PiZpdLrVigca5vF6WXp7FrmLH5Lf0cmymiNzMYYY4wxQ8aTH2OMMca0iqHIXgq7IdWtxW5QXcmSv6eryPLOs7zyL1CXStgVp/fPVgJlF5u6Y9llqKs4s6yg9+T0bJXxRtUdm6VSquuUZaNM+lS3J0sO2j7Y9c11q+Vi9za3Db2HliuD31Wfx5+zHexHDW2TLElo3bB8pe2V5eTMva7netl6/vz5tevYZrrqd7YacyYts6Sq6cNsMx0nsmU2mkRtyf1P7cXjrLZXtoMuNcHymdYLtxfufyop7ty5szrOVvdViY1XCdZ35XE3G7uzcWKU4X6i9uL+pZJY1la5fWg9sNTF4/rpp59eu45trhI6k+3Ani0zofIfl4WXkQGmvozB6I7SxhhjjDFDwJMfY4wxxrQKT36MMcYY0yqGEvOT7RydnVMNj3VCTW1m7THbLoH1YE0RZO1R0+a4nBrXw7EEnJII1FNFNW1b332mo1o0x4+orsu21W0/uH6zVFv+nsYm8LM1doB1cLUla9+63AGj57htasyBvvsooen8HMehKeUcd6BxTFynGivEMSQaj8fpsRzno3bhPp1tSaKp09yPtZ1xrITGK2zfvr3nPTXWYFTI2l2vXbqBF46DXGe6gz1/1u/x+MbtSuMze43HQD0OUvs+209jWXTHcobbR1YPo4S+T7atCI9ZGtfGdabjc3ZPPsd9T8eybGsNLov+1mUxkty/dbzgLXAGjT0/xhhjjGkVnvwYY4wxplUMZVf3bOVfdYexS1ldeCw1qOubU/zUldor3U9dalxmlc7Y5a/SVraSMbv8tR4yF+wopWDyO01k1WOuQ3WJ8n1UQsmWHOAUTJZoNJWX24pKVOy6zVYqzsqVpVWOku3GI9shPZN3uT6WLl1aO8cSktYNtx9dHZmlDb6HtjO2r44LnP6qLnTum/o9Lqcuh8DyjtpTU7ybZNOmTdWxLiPA5VQ7c5/QVHSWkHncA+rtY9euXbVz3Md5VWwtF9tk4cKF6IXKi9z3VX7hcmm/5XMq76ik2STcL3Xc47rV8YzfSZcbyXa752dom+Y2z8/WHdgZvT/3WZVPs37Jz9ZxgNujyqlZmEI/2PNjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFQOL+WHdVTVzTr/THbezeAzWBjNtU1Nae6V8aqomxw3pMuG8rLfq56ydapl1N2WGtdSpLs09TFR/Zvgd9u7d2/O6bNkCvT/HX2isEGvFbAe1OS+br/bismi8Q79bUWg6d6atjzIau8PxHvoeHLuj2j/HDGgcB8duZLbgY+3fnAaf9RXV/TluSZc84L6p78MptVpHo2RfrieNj2A05oLHQV02hO2lMVQ8lmu/5c88fmof5rrWc1wWbUdc79u2baud4/FZ2wCXa1S3DQLymFMelzQOi99Xxy+uX41v4thVjZ/pFUOl9uL2p9vOcP/Sfs9tVdsf/y6vXLmyZ7kG/Ztpz48xxhhjWoUnP8YYY4xpFQOTvTL3IksU2UqxKmWwO0zTILM0N34Gu83U1c2fNZ09e59speZeOx3r89QtPUpkKfnsUub0VqBe79mu7ppSzed0t21uE7t3766OVb7iFGrdIZzLnKXgq1ufXdGaGpql/I8aWVvmdqhuZXa3a//jvqn1xq7xTK7mOtRns+Sh/ZbblsrMvBO4nuO2pG2cV9gd5VWBV69e3fMcp8GrrLFjx47qWHdW51AEXVGX+1WWPs/HugwFtzG9x+LFi6tjlX42b95cHascx21OZclsd/tRIpMtuc3r7w2PN/p7ynWhq+Bz3Wt/5muzJUVY6tJlXrJwkGwVf36G/q5nKf9TxZ4fY4wxxrQKT36MMcYY0yo8+THGGGNMqxjKru5Ktqsvf1YNMYsh4RRo1ehZE9VYHibTVbPtGLgs2c62WuaDAX4/jXVhXVfjTDJbsgas7ePee++tjlmXVrvyPVSnznYU1s8Mx7WMcsrsVFi2bFl1rHbhVHF9f94JXOub+3HWBzjOR/tRRrZ0Bt9Ht1ng/p7FFsxUli9fXh1rnBRvM6L9luOydDkQjhPRpQO2b99eHXNcj/ZhjsGbM2dO7RyPs9pWOH4rW2riYLCdkm27w++rNmF7aZ/le+pvLd+H7aexQdz3dKkTjjNbtWpV7RyPz9mO79M5ztrzY4wxxphW4cmPMcYYY1rFtMhejLq1MjcXu/vU9cduO3X9sTud3dvqHmV3vZ7LXKnsNh71NOdhopJRliLM9tKUVraXpifzPbMlE3qtHAzUXetqV36HTAJrA1pv3D9UdmBbaDo713Fmp0xCZVtk0rKS3bNNzJs3r+c5tQlLG9k4q/2j18rQukI2SzGZ5KH24rbSZltmv1uaLs9SLy/fANSXrtAlB9jOXO+6SjS3D5XEWILWMYGXYRgVW9rzY4wxxphW4cmPMcYYY1qFJz/GGGOMaRXTHvOToTqh6pIM656qRff6nmqnk90ldlQ0y5lEFkPFu7pntuR7ZEvDK1ksT9vjfPol235CybZ/6dU3J2LPDPfNA5PFN+o5jpnUvsLxJXys9s/aQ4ZteWC0H2bbeXCsZbZdCLcBvf9kbTKKtrTnxxhjjDGtwpMfY4wxxrSKmIg7KiL2Ath+wAvNsFhUSjn1wJcdGNuycQZmS8D2HAHcNw8ebMuDi3HtOaHJjzHGGGPMTMeylzHGGGNahSc/xhhjjGkVM2byExFzIuILEbE5Iu6IiG9ExMqIWBwR64f0zCMj4osRsSkivh8Ri4fxnLbRkC1fHhE/jIhnI+INw3hGW2nInr/ZfdatEfHPEbFoGM9pGw3Z8p0RcVtErIuI70bEWcN4Thtpwp707DdERImI84f5nMkyIyY/0VmU4O8AXFtKWVZKOQvA7wCYnX9zyrwDwIOllOUA/gTAFUN+3kFPg7bcAeDtAD435Oe0igbteTOA80spLwbwFQD/bcjPO+hp0JafK6WcU0pZg44d/+eQn9cKGrQnIuJ4AO8B8P1hP2uyzIjJD4CfBvBMKeUvx/5QSllXSvkOX9SdzX6n+z/8H0bExd2/z42Ib3f/Z7E+In4qIg6NiE93P98WEe8f57m/COAz3eOvAHhVeFW8qdKILUsp20optwLovdqimQxN2fOaUsrYrrfXA5g/xHdsC03Zknc6PhaAs3AGQ1O/mwDwX9GZyD7Z43zjjNQKzwmrAdzUx3X3AfiZUsqTEbECwOcBnA/gLQC+WUr5o4g4FMAxANYAOL2UshoAIuKkce53OoCdAFBKeTYiHgYwC8D941xr+qMpW5rhMAr2fAeA/ze54huiMVtGxGUAfhPAEQAumeJ7mA6N2DMi1gJYUEr5h4j4wEDeZAjMlMlPvxwO4OMRsQbAcwBWdv9+I4C/iYjDAXytlLIuIrYAWBoRfw7gHwFcNc79xvPy+H8l08OgbWmaZSj2jIhfQWegfsUwC29qDNyWpZS/APAXEfEWAL8H4NeG/A7mxwzMnhFxCDohIm+fprJPmpkie90O4CV9XPd+APcCOBedAfEIACilfBvAywHsBvDZiHhbKeXB7nXXArgMwCfHud8uAAsAICIOA3AigAem8iKmMVua4dCYPSPiUgC/C+A1pZSnxrvGTIhR6JtfAPBLkyi7eSFN2PN4dDxO10bENgAXAfj7GMGg55ky+bkawJER8Rtjf4iICyJC/7d3IoA9pZTnAfwqgEO71y4CcF8p5RMA/hrAeRFxCoBDSilfBfD7AM4b57l/jx//D+QNAK4uXhVyqjRlSzMcGrFn17X+V+hMfO4bwnu1kaZsuYI+/jyAjQN8pzYz7fYspTxcSjmllLK4lLIYnXi815RSfjCcV5w8M0L2KqWUiHgtgD+NiA+hE0S1DcD75NIrAXw1It4I4BoAY1vcvhLAByPiGQD7AbwNnXieT3XddADw2+M8+q/RmfFuQsfj8+ZBvVNbacqWEXEBOpkPJwP4hYj4g1LK2QN8tVbSYN/87wCOA/Dl6OQg7CilvGZAr9VKGrTlu7pevGcAPAhLXgOhQXvOCLy9hTHGGGNaxUyRvYwxxhhjBoInP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVePJjjDHGmFbhyY8xxhhjWoUnP8YYY4xpFZ78GGOMMaZVHDaRi4877rgya9asYZXFHIB9+/Zh//79MYh72ZbNMkhbArZn07hvHjzYlgcXO3bsuL+Ucqr+fUKTn1mzZuHyyy8fXKnMhLjiiisGdi/bslkGaUvA9mwa982DB9vy4OKyyy7bPt7fJzT5mekcfvjhPc8999xzfd/n+eefH0RxzBTYvXt3z3PHHHNM3/c5+eSTB1EcM0UOOaS3Al9K6fs+E7nWDAePswcPmS0n0teeffbZQRRnoDjmxxhjjDGtwpMfY4wxxrQKT36MMcYY0ypmTMxPRPT8fNRRR9XOHXroodVxFkvAPPzww7XPW7du7XnuzDPPrI6PP/742rmnn366r+e1GdXyWQ/esGFD7dyJJ55YHe/atas61npmm8+bN6927tRTXxDoX3H11VdXx/Pnz6+dW7lyZc/vmR+Taf9HHnlk7TP3x+x73EaeeeaZntdp/+YYBY0zcAzJgRnEOKv3YDs/9dRTPc/p9/j+2t+zNmE6ZLbUWB7+zPWe2VJtwP1L+yWPA/q9pn4z7fkxxhhjTKvw5McYY4wxraJx2Yvdascee2ztXOamZjfd0UcfXTt32GE/fi2WrDI3nboBlyxZUh0/8cQTtXN8n+OOO65nGTWtk7+XuW1nqkuXban1wq5Urc+9e/dWx/v376+d4/vw8bZt22rXPfroo+PeDwB27NhRHb/sZS+rnWOJ7Prrr6+d27RpU3V8/vnn18696EUvqo5ZjgPq7Wrp0qWYqbCLm+2n55Qjjjhi3OPs/v3eL7uHftY+nbnlue2qNJM9b6bA78vj44HgutdxtpfslY3bKp0xOu5xXWsbYBvp8/hzZq+DYZzVNp7JiHytytHcJvhY67ZfO6uUxddq++PPmcyWyWOTSaW358cYY4wxrcKTH2OMMca0Ck9+jDHGGNMqpj3mRzVK1l337NlTO8ea8imnnFI7N3fu3OpY9ct9+/ZVx5y+rM/mmBxNZ9+5c2d1/Mgjj/R8h+3b69uGsH65aNGi2jnWM2fPnl07d88992CmoXET/O6q7fO7Z2nql156ae0c2+WBBx6ojjX+gNuRxo6tXbu2OubUeQB48sknq+Nzzjmndo5jh771rW/VzvEWGrz0AVCPB5rJsH2177A9VcPna7N0W7Z7FpuhsXN8bZbOrt/L0rH5XNa2ZkqcSJaWrtu/8Ls//vjjtXMPPvhgdXzvvffWzvHYyv0qizXR2J0s7ovvo3bmGD8tMz9P41J0LJ8JaN9je6kt+12aQO3Av6FZv3/ssceq4/vvv792jjdw1RhBHvO1D/Fnjbdju+v7TDVF3p4fY4wxxrQKT36MMcYY0yoGJntlKXC8c7a6yji1md1mQN0Vp2463tVbv8euM3YLqkuN3W3qQuu14iVQT9VWNypfqzuPb9mypTpesGBB7RxLfpdccgmahOtJJUW2g7q3s3PsmtY2wBKgPu+0006rjtkmKi2xnU866aTaOf6epqWz+1zlMv7esmXLaufYRZ+191EgS3/tN61V4e+pzVg2ylZg5ntkqfRZSq32aZa6WNLUcxNZ7ZnfoWnZi+2nMhFLT1rvXE/33Xdf7RzLCVqfvWRKAPjRj35UHfOYqM+++eabq2Pt+9xXdVkIHk9U2uK2o+M/14uGM2j7b5LsN1PHIobfXa/j5UDUDmw/tSX3y6y98/fmzJlTO5fJcZlUzf1Ul0Hhc1pHXObJSGD2/BhjjDGmVXjyY4wxxphW4cmPMcYYY1rFwGJ+WI9TvY+1OdX0OIVd9UvW8TjNGajri6phs47MOj+nbQJ1HVJToDluQ7dcYB1U9Usu58aNG2vnuJzr16/vWeam4RgO1cizrT1Yo9e0RK5DtTN/T3Vqbi/cVvTZ/JnTMYF6HIPGCrF+ru32hBNO6FkuTvvVXeO5vT/00ENomixOhO2pcQBswyzVXfsH1zfXIQAcf/zx1TG3A617Hk9Uz+fnabn4HTQeL0u3zZbY1zGrSdh+2VITGofFW7xom+QYrYlsM8LPv/3226tjHauz2IwsZZ3HZI0VYhupvbj9afvgcaLprUq4jes7cD3pWMr9RvsX20T7czbWcZtYt25ddcxxXUDdttrG+DdMYyR52RdtR/xZbZLFRek4P1Hs+THGGGNMq/DkxxhjjDGtYigrPKtLmV166qZj95t+j11s2UqT6jJktyC78zT1nMui6cnZCrZ8f3Wts9tOXX/8rurGVUluVOG6UPclv5+6JFl+VNkrSx9m1yq3AZUJOSVSV+hl16nKm4sXL66OtR3xPTVlNpMKmk6HzlDXMcu2mhrOK6Wrrdm+6v7mc+p65/rnJQ5UjuT2olIJn9P+xxKBvg+3Hx1P+BmZ1DvK8HjJq5MDwKZNm6rjbEkHrTNOI9dxkKVflq9UUuG+qv2Pz+n9+X1UquP7qCzJZdG2c+GFF1bH3L6bJls2ROuM7Zct2aAhGTwO3nTTTbVzLFvybgPaHrJ+z0sVaLn4nmpnXs4kS8HPwk8mgz0/xhhjjGkVnvwYY4wxplV48mOMMcaYVjHpmB+NueA0RY0PYA1W0/Y4lkK1aNYoVdfNtp/gHdk5vVzjNrKl/TmuQDVX5owzzqh9Zq2d31s/q4at104nGu/EZdG6zWKhsu1CGNWKOWZG78nP47go1bP52do2Fy5c2PPZ3FazHZ/1edlu1FPdbXiqaOr9N77xjer4TW96U+3cD37wg+pY2yD31SxlX89lsRRcV9q2GLaLxglyuTQ9mr+nu7rzLuS6ND/bN9vWY7rRtszt/I477qidY7trvfD7ZVuJcPwFUO8v2VYHHO+hsY5cn9o2uO1wXAhQHwu0jfH7aCwIx6UsX768dk7bxHSi23BwWbKYH/1d5PdTm3Cav45D/L3vfe97tXO9tqbQsYTvr7+Z3B60D/H76RYqfB8dgzm2UutoqksVjE4vN8YYY4yZBjz5McYYY0yrmLTspW5IlbMYdrdlLmV1u/O1KoewO0xdqfy9bGd4Pqfp0ez603O8Smi2Mq2+K7sQ1e3++te/Hk2hafaaBsxk9stWC2bUHcttR+UlbmfcjrZt21a7jlPW1SYsW6ot+Xl6jt3pWq5sh/CmVwTWNOfzzz+/OmZJGKi7krUfsbtdXc7cDrRueq2wrt/L0su57nVcYHc3rxit57SPseyt32MbqsygbvrpRMdZHkO0nNzudOkAfodsaQaV+NlGKj9yf5w7d251vGLFitp1LMHxStNAPSxBx3heXTjb0VvTsflddXkTTp9fs2YNphPtX9kyAjyG7dq1q3aO7aXfY3vpzgc8LuhYx3bnvnH66afXruPfvs2bN9fO8Tiu8wFum7oCOJMtg6JjRJbm3w/2/BhjjDGmVXjyY4wxxphW4cmPMcYYY1rFUHZ1V02Z4z80BoC1x2yXVo3X4e0hdKuIXrsbawomp75mu0+rfpntrMzfy3aH5qXh9dqmdwLP0rgZTbNk+6kWzRq9xnBw3au+zfFIrFnzDu9APc5E03VVR2ay9HyOHdAyZ0urN71btMJtUpd74DapadVsQ40L436rSyVwDJDGj3FZsqUm2NbZMhQaj8P9Xfsf92ONy2J76nIITaZHKxzfksXmZTEdWcyFxstt3769OtbxmeOKuO9oW+HP2hf5s47jbAe1Ccc+ad/kdqVtQMelJuG2q+/H/VRjd3r9vgH1utAlALjvaTwctw8+p/bivqC/YfxsjQdie2Vjp7ZbvjaLu5wM9vwYY4wxplV48mOMMcaYVjFpH6C6D7NVGlkaUvmKv6eyA6dPqhzC0pC6fzm1kt1mmrbH7np13bMUo64/ljXUfc6uOJU/+LNKaZnkN2xUumObqOTA12Yrxaprk92zWi+8squu8souX/7e6tWra9dxW9H34brVtFh2DWsqKqeFZynbek9uE+qWng44RRiou/61rPwenMYK1PuRpk7rtQy3bZVReCd3frbKV9wf1e2frXDLbSST2HX8ytpnkzKmlpPblq6+y/KSjqWcVq0SP8sVKmXwPVWO7yWn3nDDDbXruN7VXtw2tZ5Z2l60aFHtHLdjrQd+H13pWuXy6UTHfJYf9TeGU/Q3btxYO7dkyZLqWO3M9tPxmcM81A4se3F9ar9cv359dayyGi9joL+n3AZ0DsDPyFYR17JY9jLGGGOMmQCe/BhjjDGmVXjyY4wxxphWMZS8P02l5HgB1elY/9PYHb6P6ousbaouzpovxxxoXA3HEanWyBqsxr1kOwrztfqurFnq91SbHlWy7S04bkLTTzkOTO/Bum5WD6x1a/wWxwDcdttttXPcrlRb5/RajVXjOAaN7eKYF31Xfp8mYn4UjqXQPsb1pnXPMQqaBs/1qLEFHNczb9682jnuj9kO1FlK95YtW6pj7bfc57Sd8fvomMGfdZwYtaULxtA0bq5rjS/J0qo5rke3P+GxLlvaguOBtI9xvIf2FbYXx+0B9R3ZdYkGbi/Z0gS61QbXy1S3R5gqXE8aG8e/bxqnx7+LOmZxX9E4Kf4d1thKbgPcVnSrHC6X1h/bJFseJhsTtV9m2+hkKfP9YM+PMcYYY1qFJz/GGGOMaRVDkb3UTcyuOHXv6WeGU1xVSmPXHLvBgbobnlf7VdfpggULqmNN8WQ3rqb1sls/22lW3efZCsijikoHLFNmO/eqrMfn9N2z3eC5fSxdurQ61rTbW2+9tTrOpDN1ufI7aLtlV7CuPsvlVEmMd5gfNbLlCXRHb5YrtJ+yC11X1ObP2a7r3K9UIub2o22Jy6JyHLdPXS02kzwyaVAluVGFbavjEteZLlvAErK23euuu646Vkns5ptvro5ZktB2tHXr1upYJbELLrigOlaJlCUyTY/mfsxp/HpOxxqViZqE7aXl4vFG2zGPfTrWcX/Qdrtp06bqWPsbj8F8fx0vOfRA+yVfq78N/FnbAI8Damfui7pMAafWTwZ7fowxxhjTKjz5McYYY0yrmLT2opHX7F7UKGw+p24tdtPpqpCZVMLSk8L3ZJeounTZPavuc5aoVBJjOY6lGP2eRumzVKduaXaBat0Om2wV5ywzS9+By61SEH/WlTrZdauuVHaJsh22bdtWu47dxOpyZVetunu5XOqu53acSS8ql3Eb0NVnpwPdtJBXfFaJil3JWm/cV9XW3Fe17/Bn/V6vbEFt81yn6tpXmYPhtqX2zMrBfVVtrfUynWTjbCZhKtyPdHxme2nGFdeTroB+4403Vsc8hmTjntqE3y/LNlSphEMW9Hn8ORvbpjvbS8vJ76eSIkvuWmdcn/oObFvN1GI5SyUklsi4PWjIB99fn82ypY4zWSZir82Ogfr4oX12qra058cYY4wxrcKTH2OMMca0Ck9+jDHGGNMqBpZvzdrtRRddVDvHaXzZ6p8aC/LiF7+4OtYYH059v/jii2vnWLfm5+n9ucybN2+uncvSrzllntM4gXpqXrYytKYvsraaadjTAZdbNViOmdGU4GzlTn4HjQPhWAKNt2D9mWNXVDdetWpVdaxxC6w3q525HenSBFkaPMcKaazaVVddVR2/+93vrp3TNjgdcAyQ2ohtrXXKmnq28qrWDdtQ40v4+dx+NN6Dd7XWlN09e/b0PMf31LizbIfobKf4bCdwba/DhscQbXdch2pntpHGTmRLQ7DdNfaEz3Ha8RlnnFG7jlPruZ/q52z1f22bbC8dZ7kNqJ15nNVxXWOahg2PiXfddVftHMen6rjEKzfrbwO/n8b8sN1POOGE2jm2Edc77wQP1JeV0Xg0jt3R7/E5/e1jdCzhMuvq4PzbfuWVV9bOvfWtb+35jOreB7zCGGOMMeYgwpMfY4wxxrSKaVnhmV1X6vJit7GeY7fWmWee2fN52aak7MLLVrVU1y+vXqqr27KbVdPgOcVPpQJ20euKnuxSnm6ZayKwLTM3tbrW+drMha3wtSzfaDo3SwCaZsn20g39+P7qVmV7qYSSrTz6xje+sTpuQubKyNKh1WbcflXi5PfK7KdLW3A/4/7Hq88C9T6g9+e2pH2T3esq/WRSNrcZrQduB6O02vNElqjIJDG+z/bt22vnWOZQyZj7C49fOo6zTVSS5v6n5crkU27HahN+H23v3K5GaZx9yUteUvucrUSdbbqd1Sdfq+2fZbBeYSNA3Sb67GzVbX6e3pPbkY7B3L+1DfDnfmQuxZ4fY4wxxrQKT36MMcYY0yo8+THGGGNMqxhKzA8vzQ3UdUiNneBlvXX7CdYeNTWP06U1XofjDFgXzLYhyHYT19Rs3VGc4XLqPbMl4PX9RgXVdbkuVGvnc9kWC1p/nJKpMQG9luXXXda5jXH6OlBvHxrTwBq2xqdkSxPwZ9XWeZds1bCbJrOnxnDxe2XLNmi6N8cHaZ/j+tiwYUN1rGm5HDOg5eK619gC3k5Ev5dte8DXZu86SmgsFMdoZcv969jDKd633XZb7RzHYmlbXrhwYXW8fPnynvfn+CMtM28VpOM/t50s7iuzs/bbXturNI2Wk2Oo9DeTx91sKRLtG9wv9TeNf++47WgsHseA6W8Wp+CrTbL4QX4fjQfl++jvKY9P2XINvRjNlmCMMcYYMyQ8+THGGGNMqxiY7MVuNHXF8Wd1o/EKyeoOY7edurz4WnWxsYuUjzXtmNMJWZYB6qmvu3btqp3LdilmF6W6MvncdO8oPBH63UlX5QA+p7sBc7rr7bffXjvH9lP3JUtYmaubV5jVcvFqwSqvsJ01BZNdsPo8bpu8gjTwQvmsabLVfflcJompm5z7n9YbX6v3ZFtw6rSmrGt/Z7htaf/jule5oN8+N6oyF5D3zQzuf5pGzfKjprPz87Rds8xx7rnn9nz2tm3bep5j2YslS0VtyRJ4luo+U8ZZbe9c12w7oP7uuio1S1H6Pb6nSpNcvzzW6VIuPO5pn2UZVOVufr9Mcs7atI5BOs+YKPb8GGOMMaZVePJjjDHGmFbhyY8xxhhjWsVQUt0zNKaDY22yuB5Ns+R4DI0V4mewtqm6Md8zW/Y/S5FXXZV1SX0e7wav+rku3T3T0TgQ3jVY09k55kDTxntto6B1y9tbaCo9p2LrTti8NYWmhvKzOZVXy6kxB6O0bL6ienu25D23e+1j3M61f3DMjNqJY664zWt/YJtp+ivbUNsZ2yyzg5Y5WxKDYxI07mBU0TrjWArd4oX7lcZQcTwbx+cA9eUsON5Dx3huV7oMBfcj7bcc06F9jD9rGjy3aa2HUY7nYric+g48Tmm6Ode1/jZxnWk8Ksdbcd/TGBwexzXWkWOANLaS30Fjdbjvaxvgcur4NNVlROz5McYYY0yr8OTHGGOMMa1iKLKXpgWzC49XvwXqaZCaEskp5iqVsFyh92RXGbvt1B3L7jd1u/Mq1bpiNd9H3bHsilP3Hn9vVFcaVdTNz+XW9zvppJOqY5UO2H579uypnVu1alXPe7It2R2rLk92yes9WBrRHd/Zda+yF8sm+rwszTLb5bxpsuUJ9Bx/1v7BcoW2Ee7/6sZmyYXd5LNnz65dx3KZLpvA0kwmH2crrGuZ+035HyW0nPy+LLED+fILXL9Lly6tnctkB5ZcuN61vvj+Wb1rv+H+rmMw31PlzV4r/I/3/FEhk6P1HNeZhgnwWKf9hvu6/v7w7xgvhaDyGKPPzpbU4HrX7/F4ocvR8G+KjsFZqEo/zIxfYGOMMcaYAeHJjzHGGGNahSc/xhhjjGkVQ4n5UZ2VdTzVnlnT0y0mWBtUDZt1fz3HWqDqngxrm3oP/qy7VvOzVW/md9U0RNY2NU5pVJdhV41cY5wYfj9NmeUYrZUrV9bOsZ05lRKo1y/HiGiaJde7xibwtbwkP1DXyNUGrDFruimj6dyjjOrkWawL9+MsDiZLKdfUZo6V4vRaTcvl+tZ+xDbT+MKsH3G70PimyW4b0ST6rtn2MjzW6ftxXIfGeHC96DmO0bn77rurYx7Tgbpt9dz69eurYx1bst3guc/p8hUc1zdT+qb+ZmaxeGx33aqExyztU/xZ2wfXPY+lOl7w75b+lnPf1vdhG+n4zNdm6eu6ncZUmRm93BhjjDFmQHjyY4wxxphWMTDZi6ULdauyK0tdV3ytpiGzi03d2/w9XXmUU+T5OnUTs9tOV4plF6G+D+90q3IIp3Sre2/ZsmXVsbqQp7pa5XSRrZCapfKz3VX2yqQndrtzKrrai+UVdceyBKfuWHYbq4uc7ZctW5CtjDwKcBvNVhLP3kP7H79/tgK69k1uB1kaLT9PV5fmsUYlaW6fE5G2+J6jLJWwjVSS5s8qh3CfUMmDban2YjkkCyG46667el7HUvbOnTtr57Zv314dL168uHaOfw+0/3E59XvZruQ6bjRJtvxAJkdz29VwDe2nve6p/YbbPI+l+vvGttQ+xPK3LivD9tJla7Id33mlaF1lX58xUez5McYYY0yr8OTHGGOMMa3Ckx9jjDHGtIqBxfz0mx6qejrH1mjcC2uIGlfAy+Zr7ABfy/dYsGBB7TrWprOduVVH5RgS1S+5LKpfcrzDqKa2A/XYAU0xZe1d9fPMXqwdq72WLFlSHasOznXPz9Yl0jmVUmMhdMl0hvVm1bf5ffQeXK4sjmYU0Fgehu2kdc/f03aQLWXP/fj++++vneM2w/Wm9cvxOtpXuO71HLcLjfnJdmfP0vpHCS639rF+0/W173C7Vzuz/XTs5rgitrluRfGDH/ygOuYYHy2Lbi/DdsjSo3fv3l07N3/+/Or4lltuqZ27+OKLMSpkW21k/YvrRdPZ2Q4a98V1pvEy/Ay+TpeqYNvqciYbNmyojtUm3MY0rozHdd0NnseFFStWYJDY82OMMcaYVuHJjzHGGGNaxcBkryztjKWhLNWdVwEG6tLWjh07auey1VrZtc4uUHXTcQqtSm4szeiqwOyuzOQQlYX4nKYojhKZ25/tpVIh20GlIJaQ1F7s9sxS3fl56oLvdR1Qr3eVMPlalREyuYgloqmmXA6brHz8HmoXtmGWUqt1yv1WZS92m7NUme3Orn2MV4LX1aXZZmq/rO1mK+qOEpnEyjK+poaznJW9X7bkRyZN8vN0lXYeZ1Va5vtnKes6LrCkk0myoyRzKdnK6Pw7qfXObSALS1B4HFA5a86cOdXxPffcUx1re+B73HrrrbVzvMq3/p7y+Kzj0erVq3uWa926ddXx2WefjUFiz48xxhhjWoUnP8YYY4xpFZ78GGOMMaZVDGVXd2Xfvn3VsWqUnOqoeiXrjatWrep5TvVgXrqb9UWNE+HYBF1Cn9MuszTLbLl0jSsY5RTafmH9WVM12X6qU3NMhWrd/L1syf5syXe+rt97KPq97HkHC9xGNaX2pptuqo41PZrrQzV8jgXQ1GZOV2UNn7ekAYCNGzdWx7wUAgDMmzevOtb+x7ES2Q722TYAMxWOl9DYSh77NL6K4z3UzrwVjdZZr7gUjUvietexlNsfx5oAwObNm6tjbWM8xmus0KhtLzMZ+H21PvmzbhfCdtY6499e7escm8d9Vpcf4Jg9jvEB6nE+vN0JUI/t0nvyszWOaJjjrj0/xhhjjGkVnvwYY4wxplVMi+zFqOTBn1V6YvesymWcZqluTj7H7j1Ns+R0dpWoOLU3u//B4GKdLPru7I6diCTGrk39Hp/rdazfy1ylWmZ25bfZlsAL623t2rXVsbrXr7/++upYd/Fmt7mmrvZ63ve+973auU2bNvX83uzZs6tjlbLZ1Z8tVXCwky1voPJVFhqgq7EzLHOwLTWVPktz5mtPOOGE2jlOgc76po7d+nmmk8leLFkC9b6oKylz3atczHZg6VPDTViyUrlx/fr11fGaNWtq53i5mGx8vuiii3qeGzT2/BhjjDGmVXjyY4wxxphW4cmPMcYYY1rFtMf8ZGRLt6uOq8tnM72W+B7UsvVtjw3phyy2RsnszrYcRqqybdkfCxYsqH3WuC3mZS97Wc/rtm7dWh1zbNCrX/3q2nUc/6dxglnf1y1yzAvJYnKULHaIbcT9yOPs9KH9i+PcdMd3jvPRVHf+/NKXvrQ61t9d3nVdY2gvvPDCnuXMYv+awp4fY4wxxrQKT36MMcYY0ypiIq7FiNgLYPsBLzTDYlEp5dRB3Mi2bJyB2RKwPUcA982DB9vy4GJce05o8mOMMcYYM9Ox7GWMMcaYVuHJjzHGGGNaxYyZ/ETEnIj4QkRsjog7IuIbEbEyIhZHxPoD32FSz3x7ROyNiHXdf78+jOe0jSZs2X3uv+s+7/aI+NywntM2Guqbf0L9ckNEPDSM57SNhmy5MCKuiYibI+LWiPi5YTynjTRkz0UR8c9dW14bEfOH8ZypMlLr/PQiOgtR/B2Az5RS3tz92xoAswHsTL46CL5YSnnXkJ/RGpqyZUSsAPDbAH6ylPJgRJw2rGe1iabsWUp5P5Xh3QDWJpebPmhwnP09AF8qpfyviDgLwDcALB7i81pBg/b8HwD+dynlMxFxCYCPAPjVIT5vUswUz89PA3imlPKXY38opawrpXyHL+rOZr8TET/s/ru4+/e5EfHt7v8S10fET0XEoRHx6e7n2yLi/TDTQVO2/A0Af1FKebD7zPvGucZMnFHom78M4PMDf7P20ZQtC4CxXU1PBHDPkN6vbTRlz7MA/HP3+BoAvzik95sSM8LzA2A1gJv6uO4+AD9TSnmy+z/9zwM4H8BbAHyzlPJHEXEogGMArAFweillNQBExEk97vn6iHg5gA0A3l9KGban6WCnKVuu7J77VwCHAvhwKeWfpvguptm+iYhYBGAJgKun8hIGQHO2/DCAq7oevGMBXDrF9zAdmrLnLQBeD+BjAF4L4PiImFVK2TfF9xkoM2Xy0y+HA/h417X3HLo/eABuBPA3EXE4gK+VUtZFxBYASyPizwH8I4Crxrnf1wF8vpTyVES8E8BnAFwy7JcwAAZvy8MArADwSgDzAXwnIlaXUh4a6luYMQZtzzHeDOArpZTB7Hdi+mHQtvxlAJ8upfxxRLwUwGe7ffP5ob+JAQZvzw907/d2AN8GsBvAYPY8GSAzRfa6HcBL+rju/QDuBXAuOjPXIwCglPJtAC9HxwifjYi3deWPcwFcC+AyAJ/Um5VS9pVSxjZP+USfZTA5jdgSwC4A/7eU8kwpZSuAu9GZDJmp0ZQ9x3gzLHkNiqZs+Q4AX+re4zoARwHwJm1Tp6nfzXtKKa8rpawF8Lvdvz085bcZMDNl8nM1gCMj4jfG/hARF0TEK+S6EwHs6f6P4VfRkTfGXOP3lVI+AeCvAZwXEacAOKSU8lUAvw/gPH1oRMylj68BcOcA36mtNGJLAF9DRwNH9/qVALYM8sVaSlP2REScAeBkANcN+J3aSlO23AHgVd17nInO5GfvONeZidHU7+YpETE2t/htAH8z4PcaCDNC9iqllIh4LYA/jYgPAXgSwDYA75NLrwTw1Yh4IzqBVo91//5KAB+MiGcA7AfwNgCnA/iUGEl5T0S8Bh2X3QMA3j6gV2otDdrymwB+NiLuQMe1+8FR06BnIg3aE+jIJV8oXqZ+IDRoy/8E4BPd4NkC4O226dRp0J6vBPCRiCjoyF6XDeiVBoq3tzDGGGNMq5gpspcxxhhjzEDw5McYY4wxrcKTH2OMMca0Ck9+jDHGGNMqPPkxxhhjTKvw5McYY4wxrcKTH2OMMca0Ck9+jDHGGNMq/j/jTmsIdTE5PAAAAABJRU5ErkJggg==\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",
    "<br><span style=\"color:blue\">\n",
    "    <b>sag</b>\n",
    "<code>Sparsity with L2 penalty: 9.57%\n",
    "Test score with L2 penalty: 0.8966</code>\n",
    "<br> <br>  \n",
    "<b>saga</b>\n",
    "    <code>Sparsity with L1 penalty: 9.57%\n",
    "Test score with L1 penalty: 0.8904</code>    \n",
    "</span>\n",
    "    \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",
    "![NB7_1.png](https://github.com/DBustamante98/ML/blob/main/NB7_1.png?raw=true)        \n",
    "\n",
    "<br><span style=\"color:blue\">\n",
    "What the logistic regression does for each image is accept 784 inputs and multiply them by weights to generate the prediction. Due to the direct mapping between input and output (i.e., no hidden layer), the value of each weight corresponds to how much each of the 784 inputs is taken into account when calculating the probability of each class. By taking the weights of each class and transforming them into a 28×28 matrix (the resolution of the image), we can tell which pixels are the most important for the calculation of each class. How does a person draw a 0? He draws a circular shape that is empty in the middle. That is exactly what the weights for class 0 represent. To recognize zeros, no fancy filters are needed. The same for 1: in most cases, it corresponds to a vertical line in the middle of the image. The rest of the numbers are a bit more difficult, which is what prevents logistic regression from achieving greater than 90% accuracy. \n",
    " </span>\n",
    "    \n",
    "<li> Try to vary the size of the training set and see how that affects the accuracy.\n",
    "<br><span style=\"color:blue\">    \n",
    "<b>train_size = 50000</b>\n",
    "<code>Sparsity with L2 penalty: 9.57%\n",
    "Test score with L2 penalty: 0.8966</code>\n",
    "<br> <br>  \n",
    "<b>train_size = 5000</b>\n",
    "    <code>Sparsity with L2 penalty: 16.07%\n",
    "Test score with L2 penalty: 0.8894</code>  \n",
    "<br> <br>  \n",
    "<b>train_size = 500</b>\n",
    "    <code>Sparsity with L2 penalty: 26.66%\n",
    "Test score with L2 penalty: 0.8293</code> \n",
    "    </span>    \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",
    "<br><span style=\"color:blue\">\n",
    "$$\n",
    "\\textbf{W} = \\{\\textbf{w}_k\\}^{K−1}_{k=0}\n",
    "$$    \n",
    "$$\n",
    "   P(\\mathcal{D}|\\textbf{W}) = \\prod_{i=1}^n \\prod_{m=0}^{K-1} \\left[P(y_{im} = 1|\\textbf{x}_i, \\textbf{W})\\right]^{y_{im}} \n",
    "$$\n",
    "</span>   \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",
    "<br><span style=\"color:blue\">\n",
    "\n",
    " $$ \n",
    "   E(\\textbf{w}) = −\\sum^n_{i=1} \\sum^{K−1}_{m=0} y_{im} \\log P(y_{im} = 1|\\textbf{x}_i, \\textbf{W})\n",
    "$$\n",
    "</span>\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",
    " <br><span style=\"color:blue\">\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "  \\dfrac{ \\partial P(y_{im;m \\neq k} = 1|\\textbf{x}_i, \\textbf{W}) }{\\partial \\textbf{w}_k} &=- \\textbf{x}_i P(y_{im} = 1|\\textbf{x}_i, \\textbf{W}) P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W}) \\\\ \n",
    "  \\dfrac{ \\partial P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W}) }{\\partial \\textbf{w}_k} &=\\textbf{x}_i P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W}) \\left( 1- P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W}) \\right) \\\\\n",
    " \\end{align}\n",
    "$$       \n",
    "Then\n",
    " $$ \n",
    "\\begin{align}\n",
    "  \\dfrac{ \\partial E(\\textbf{w}) }{\\partial \\textbf{w}_k} &= \\sum^n_{i=1} \\textbf{x}_i \\left[ -y_{ik}\\left( 1-P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W})\\right) + P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W}) \\sum^{K−1}_{m=0;m \\neq k} y_{im}   \\right] \\\\\n",
    " &= \\sum^n_{i=1} \\textbf{x}_i \\left[ -y_{ik}\\left( 1-P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W})\\right) + P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W}) (1- y_{ik})   \\right] \\\\\n",
    "  &= \\sum^n_{i=1} \\textbf{x}_i \\left[ P(y_{ik} = 1|\\textbf{x}_i, \\textbf{W}) -y_{ik}  \\right] \\\\\n",
    " \\end{align}\n",
    "$$   \n",
    "    \n",
    "</span>   \n",
    "</ul>\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}