unit05_optim/demo05_2_grad_descent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo:  Gradient Descent Optimization\n",
    "\n",
    "In the [breast cancer demo](../logistic/breast_cancer.ipynb), we used the `sklearn` built-in `LogisticRegression` class to find the weights for the logistic regression problem.   The `fit` routine in that class has an *optimizer* to select the weights to best match the data.  To understand how that optimizer works, in this demo, we will build a very simple gradient descent optimizer from scratch.  You will learn to:\n",
    "* Compute the gradients of a simple loss function and implement the gradient calculations in python\n",
    "* Implement a simple gradient descent optimizer\n",
    "* Visualize the effect of the learning rate in gradient descent\n",
    "* Implement an adaptive learning rate algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the Breast Cancer Data\n",
    "\n",
    "We first load the standard packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "%matplotlib inline\n",
    "from sklearn import preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next load the data from the [breast cancer demo](../logistic/breast_cancer.ipynb).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id</th>\n",
       "      <th>thick</th>\n",
       "      <th>size_unif</th>\n",
       "      <th>shape_unif</th>\n",
       "      <th>marg</th>\n",
       "      <th>cell_size</th>\n",
       "      <th>bare</th>\n",
       "      <th>chrom</th>\n",
       "      <th>normal</th>\n",
       "      <th>mit</th>\n",
       "      <th>class</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1000025</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1002945</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>7</td>\n",
       "      <td>10.0</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1015425</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1016277</td>\n",
       "      <td>6</td>\n",
       "      <td>8</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>4.0</td>\n",
       "      <td>3</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1017023</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        id  thick  size_unif  shape_unif  marg  cell_size  bare  chrom  \\\n",
       "0  1000025      5          1           1     1          2   1.0      3   \n",
       "1  1002945      5          4           4     5          7  10.0      3   \n",
       "2  1015425      3          1           1     1          2   2.0      3   \n",
       "3  1016277      6          8           8     1          3   4.0      3   \n",
       "4  1017023      4          1           1     3          2   1.0      3   \n",
       "\n",
       "   normal  mit  class  \n",
       "0       1    1      2  \n",
       "1       2    1      2  \n",
       "2       1    1      2  \n",
       "3       7    1      2  \n",
       "4       1    1      2  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "names = ['id','thick','size_unif','shape_unif','marg','cell_size','bare',\n",
    "         'chrom','normal','mit','class']\n",
    "df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/' +\n",
    "                 'breast-cancer-wisconsin/breast-cancer-wisconsin.data',\n",
    "                names=names,na_values='?',header=None)\n",
    "df = df.dropna()\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the breast cancer demo, we create a data matrix `X` of various features of the breast cancer sample.  The response vector `y` is a binary indicating if each sample is benign or malignant.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Get the predictors\n",
    "xnames = names[1:-1]\n",
    "X = np.array(df[xnames])\n",
    "Xs = preprocessing.scale(X)\n",
    "\n",
    "# Get the response.  Convert to a zero-one indicator \n",
    "yraw = np.array(df['class'])\n",
    "BEN_VAL = 2   # value in the 'class' label for benign samples\n",
    "MAL_VAL = 4   # value in the 'class' label for malignant samples\n",
    "y = (yraw == MAL_VAL).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to learn the classification rule to predict `y` from `X`.  To this end, we first split the data into training and test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ntr = 400        # number of training samples\n",
    "Xtr = X[:ntr,:]  # training\n",
    "ytr = y[:ntr]\n",
    "Xts = X[ntr:,:]  # test\n",
    "yts = y[ntr:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing the Gradient and Loss Function\n",
    "\n",
    "Recall that training a logistic function means finding a weight vector `w` for the classification rule:\n",
    "\n",
    "    P(y=1|x,w) = 1/(1+\\exp(-z)), z = w[0] + w[1]*x[0] + ... + w[k]x[k-1]\n",
    "    \n",
    "From the class notes in logistic regression, we saw that the weight vector can be found by minimizing  the negative log likelihood.  The negative log likelihood is called the *loss* function.  For the logistic regression problem, the loss function simplifies to:\n",
    "\n",
    "\n",
    "    f(w) = \\sum_i -y[i]z[i] + ln(1+exp(z[i])),  z[i] = w[0] + w[1]X[i,0] + ... + w[k]X[i,k-1].\n",
    "    \n",
    "To run a numerical optimizer in python, we need to implement a python function that returns `f(w)` and its gradient.  That is, we would like a function of the form:\n",
    "\n",
    "    def feval(w):\n",
    "        ...\n",
    "        return f, fgrad\n",
    "\n",
    "One issue is that the function `f(w)` also depends parameters `X` and `y`.  So, we need a way to pass those parameters to the function.  There are two approaches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 1.  Use a class method\n",
    "\n",
    "One approach to pass the parameters is to implement the function as a method in a class that contains the parameters.  In the code below, we create a class `LogisticFun` where the constructor takes the parameters `X` and `y`.  Then, the class method `feval` takes a variable `w` and computes the loss and its gradient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class LogisticFun(object):\n",
    "    def __init__(self,X,y):\n",
    "        \"\"\"\n",
    "        Class for computes the loss and gradient for a logistic regression problem.\n",
    "        \n",
    "        The constructor takes the data matrix `X` and response vector y for training.\n",
    "        \"\"\"\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        n = X.shape[0]\n",
    "        self.A = np.column_stack((np.ones(n,), X))\n",
    "\n",
    "    def feval(self,w):\n",
    "        \"\"\"\n",
    "        Compute the loss and gradient for a given weight vector\n",
    "        \"\"\"\n",
    "        # The loss is the binary cross entropy\n",
    "        z = self.A.dot(w)\n",
    "        py = 1/(1+np.exp(-z))\n",
    "        f = np.sum((1-self.y)*z - np.log(py))\n",
    "        \n",
    "        # Gradient\n",
    "        df_dz = py-self.y\n",
    "        fgrad = self.A.T.dot(df_dz)\n",
    "        return f, fgrad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have defined the class we can create an *instance* of the class with the training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "log_fun = LogisticFun(Xtr,ytr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we can call the method with only the value `w`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some random point\n",
    "p = X.shape[1]+1\n",
    "w0 = np.random.randn(p)\n",
    "\n",
    "# Call the function\n",
    "f, fgrad = log_fun.feval(w0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 2.  Using Lambda functions\n",
    "\n",
    "An alternative method is to use python `lambda` functions.  The code below creates a function `feval_param` that has all the parameters.  Then, we use the `lambda` construction to fix the parameters `X` and `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a function with all the parameters\n",
    "def feval_param(w,X,y):\n",
    "    \"\"\"\n",
    "    Compute the loss and gradient given w,X,y\n",
    "    \"\"\"\n",
    "    # Construct transform matrix\n",
    "    n = X.shape[0]\n",
    "    A = np.column_stack((np.ones(n,), X))\n",
    "\n",
    "    # The loss is the binary cross entropy\n",
    "    z = A.dot(w)\n",
    "    py = 1/(1+np.exp(-z))\n",
    "    f = np.sum((1-y)*z - np.log(py))\n",
    "\n",
    "    # Gradient\n",
    "    df_dz = py-y\n",
    "    fgrad = A.T.dot(df_dz)\n",
    "    return f, fgrad\n",
    "\n",
    "# Create a function with X,y fixed\n",
    "feval = lambda w: feval_param(w,Xtr,ytr)\n",
    "\n",
    "# You can now pass a parameter like w0\n",
    "f0, fgrad0 = feval(w0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing the gradient\n",
    "\n",
    "Either way you create the function, it is very important to test if the gradient is correct.  This is the number one reason people's code does not work with numerical optimizers.  The simplest method is to take two points `w0` and `w1` that are close to one another and then verify that \n",
    "\n",
    "    f(w1) - f(w0) \\approx fgrad(w0).dot(w1 - w0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Actual f1-f0    =  -4.6214e-03\n",
      "Predicted f1-f0 =  -4.6214e-03\n"
     ]
    }
   ],
   "source": [
    "# Take a random initial point\n",
    "p = X.shape[1]+1\n",
    "w0 = np.random.randn(p)\n",
    "\n",
    "# Perturb the point\n",
    "step = 1e-6\n",
    "w1 = w0 + step*np.random.randn(p)\n",
    "\n",
    "# Measure the function and gradient at w0 and w1\n",
    "f0, fgrad0 = log_fun.feval(w0)\n",
    "f1, fgrad1 = log_fun.feval(w1)\n",
    "\n",
    "# Predict the amount the function should have changed based on the gradient\n",
    "df_est = fgrad0.dot(w1-w0)\n",
    "\n",
    "# Print the two values to see if they are close\n",
    "print(\"Actual f1-f0    = %12.4e\" % (f1-f0))\n",
    "print(\"Predicted f1-f0 = %12.4e\" % df_est)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the two agree well.  You would have got the same answer had you used the `feval` constructed in method 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Simple Gradient Descent Optimizer\n",
    "\n",
    "Now, we build a simple gradient descent optimizer function with a fixed learning rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def grad_opt_simp(feval, winit, lr=1e-3,nit=1000):\n",
    "    \"\"\"\n",
    "    Simple gradient descent optimization\n",
    "    \n",
    "    feval:  A function that returns f, fgrad, the objective\n",
    "            function and its gradient\n",
    "    winit:  Initial estimate\n",
    "    lr:     learning rate\n",
    "    nit:    Number of iterations\n",
    "    \"\"\"\n",
    "    # Initialize\n",
    "    w0 = winit\n",
    "    \n",
    "    # Create history dictionary for tracking progress per iteration.\n",
    "    # This isn't necessary if you just want the final answer, but it \n",
    "    # is useful for debugging\n",
    "    hist = {'w': [], 'f': []}\n",
    "    \n",
    "    # Loop over iterations\n",
    "    for it in range(nit):\n",
    "\n",
    "        # Evaluate the function and gradient\n",
    "        f0, fgrad0 = feval(w0)\n",
    "\n",
    "        # Take a gradient step\n",
    "        w0 = w0 - lr*fgrad0\n",
    "        \n",
    "         # Save history\n",
    "        hist['f'].append(f0)\n",
    "        hist['w'].append(w0)\n",
    "\n",
    "    # Convert to numpy arrays\n",
    "    for elem in ('f', 'w'):\n",
    "        hist[elem] = np.array(hist[elem])\n",
    "    return w0, f0, hist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now run the gradient descent starting from a random initial condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEECAYAAAAoDUMLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHohJREFUeJzt3XmYVOWVx/HvQQREkU0FBAVcUFwiRkM0mthqALcAjsoD\njoFGCBmTKDHGCJoRo9HBmIm4wrg2mMUlBkUEQYTSMUOUZwJBIwIT0wRciFFAgrggZ/54b2OJvVR3\nV9WtW/f3eZ56uuqte7tO+2Kfvufc+15zd0RERFrEHYCIiJQGJQQREQGUEEREJKKEICIigBKCiIhE\nlBBERATIMSGYWbWZ/cnMlprZi9FYRzObb2YrzWyembXP2n6ima02sxVmNjBr/ItmttzMVpnZlPz/\nOCIi0lS5HiFsByrc/Wh37x+NTQAWuPshwEJgIoCZHQYMA/oCpwN3mplF+0wFxrh7H6CPmQ3K088h\nIiLNlGtCsFq2HQJMj55PB4ZGzwcDD7r7NnevBlYD/c2sK9DO3ZdE283I2kdERGKWa0Jw4GkzW2Jm\nY6OxLu6+HsDd3wL2ica7A2uz9n09GusOrMsaXxeNiYhICWiZ43YnuPubZrY3MN/MVhKSRLa8rYFh\nZlpPQ0SkCdzdGt6qdjkdIbj7m9HXt4HHgP7AejPrAhCVg/4ebf46sF/W7j2isbrG6/rMWB+TJk0q\nie/XmP1y2ba+bep6rzHj+f7vluT5K+bcldP8lcLcFXv+GvteXds3V4MJwczamtke0fPdgYHAS8As\noDLabBTwePR8FjDczFqZWW/gIOBFD2WlTWbWP2oyj8zap+RUVFSUxPdrzH65bFvfNnW915jx6urq\nBmMohlKYv2LOXX3vJ23+SmHuGrtfof7fq+u9fP83qmENZZXol/pMQkmoJfArd59sZp2Ahwl/9a8B\nhrn7xmificAY4GNgvLvPj8aPAaqANsAcdx9fx2d6PrKdFF9lZSVVVVVxhyFNpPlLNjPDm1EyajAh\nxEEJIbkymUzB/nqRwtP8JZsSgoiIAM1PCFq6QvIqk8nEHYI0g+Yv3ZQQREQEUMlIRKRsqGQkIiJ5\noYQgeaUadLJp/tJNCUFERAD1EEREyoZ6CCIikhdKCJJXqkEnm+Yv3ZQQREQEUA9BRKRsqIcgIiJ5\noYQgeaUadLJp/tJNCUFERAD1EEREyoZ6CCIikhdKCJJXqkEnm+Yv3ZQQREQEUA9BRKRsqIcgIiJ5\noYQgeaUadLJp/tJNCUFERAD1EEREyoZ6CCIikhdKCJJXqkEnm+Yv3ZQQREQEUA9BRKRsqIcgIiJ5\noYQgeaUadLJp/tJNCUFERAD1EEREyoZ6CCIikhdKCJJXqkEnm+Yv3ZQQREQEKOEewscfOy1bxh2J\niEhyFK2HYGYtzOyPZjYret3RzOab2Uozm2dm7bO2nWhmq81shZkNzBr/opktN7NVZjalvs/71rdg\n+/am/EgiItIUjSkZjQdeyXo9AVjg7ocAC4GJAGZ2GDAM6AucDtxpZjUZayowxt37AH3MbFBdH7Zq\nFVx2GZTgAYzUQzXoZNP8pVtOCcHMegBnAPdkDQ8BpkfPpwNDo+eDgQfdfZu7VwOrgf5m1hVo5+5L\nou1mZO3zObNnw8KFcN11uf4oIiLSHLlW6W8GLgfaZ411cff1AO7+lpntE413BxZnbfd6NLYNWJc1\nvi4ar1XHjjB/Ppx4InToAJdckmOkEquKioq4Q5Bm0PylW4MJwczOBNa7+zIzq6hn07wWdyorK+nV\nqxdnnglXX92BN9/sx3/8R/j4msPamn+8eq3Xeq3XaXydyWSoqqoCoFevXjRXg2cZmdkNwAWEv/B3\nA9oBM4FjgQp3Xx+Vgxa5e18zmwC4u98Y7f8UMAlYU7NNND4cOMndL6rlMz9zpfKKFXDKKTBtGgwZ\n0uyfWQook8ns+IcryaP5S7aCn2Xk7le6+/7ufgAwHFjo7t8EngAqo81GAY9Hz2cBw82slZn1Bg4C\nXnT3t4BNZtY/ajKPzNqnXn37hp7Ct74V+goiIpJ/jboOwcxOAi5z98Fm1gl4GNiP8Nf/MHffGG03\nERgDfAyMd/f50fgxQBXQBpjj7uPr+Jxa1zJ69lk477yQHPr3z/2HFBFJg+YeIZTshWl1xTV7Nowd\nCwsWwBFHFDkwEZESlrrF7c46C37xCzjtNHjttbijkZ3VNLwkmTR/6ZbIxSHOPx82bYIBA+D556Fb\nt7gjEhFJvsSVjLLdcAP85jeht9CpUxECExEpYanrIWRzhx/9CP77v0NPYY89ihCciEiJSl0PIZsZ\n/OxncOSRMHQofPBB3BGJatDJpvlLt0QnBAhJYdq0UDIaMQK2bYs7IhGRZEp0ySjbRx/B4MHQtSvc\ndx+0SHyqExFpnFSXjLK1agWPPgqrV8MPfqBls0VEGqtsEgLA7ruHC9cWLYJrr407mnRSDTrZNH/p\nlsjrEOqz87LZ42tdHENERHZWNj2Ena1ZA1/9arjBzqhReQpMRKSENbeHUHZHCDV69oR588Ky2e3b\nh9NSRUSkbmXVQ9hZzbLZ48bBM8/EHU06qAadbJq/dCvrhABwzDHwyCMwfDi88ELc0YiIlK6y7SHs\nbPZsGDMmHClo2WwRKUe6DiFHZ50FN9+sZbNFROqSmoQAYdnsq64Ky2a/8Ubc0ZQn1aCTTfOXbmV7\nllFdLroINmyAgQPDstmdO8cdkYhIaUhNDyFbzbLZzz0Xls1u165gHyUiUjSpvh9Cc7iH01Ffew2e\nfBLatCnox4mIFJyayk1Us2x2585aNjufVINONs1fuqU2IQDssgv88pewdSuMHQvbt8cdkYhIfFJb\nMsq2ZQsMGgTHHhtOTbUmH3CJiMRHJaM8qFk2O5PRstkikl5KCJEOHcJieL/6FdxyS9zRJJdq0Mmm\n+Uu31F2HUJ8uXeDpp8Oy2R07wsiRcUckIlI86iHUYsWKkBSeew4OOyy2MEREGkU9hALo2xeuvz7c\nWEeno4pIWigh1GHcOOjUCSZPjjuSZFENOtk0f+mmhFAHM7j3Xrj1Vli2LO5oREQKTz2EBsyYAT//\nOSxZAq1bxx2NiEjd1EMosG9+Ew44QNcniEj5U0JogBn813+F8pFuwdkw1aCTTfOXbkoIOejSBW67\nLVyX8P77cUcjIlIY6iE0wogR0LVrWO9IRKTUFLyHYGatzewFM1tqZi+Z2aRovKOZzTezlWY2z8za\nZ+0z0cxWm9kKMxuYNf5FM1tuZqvMbEpTg47L7bfDww+HO62JiJSbBhOCu38InOzuRwP9gNPNrD8w\nAVjg7ocAC4GJAGZ2GDAM6AucDtxptmP90KnAGHfvA/Qxs0H5/oEKqXPn0E8YPRo2b447mtKkGnSy\naf7SLacegrvXVM5bE9Y/cmAIMD0anw4MjZ4PBh50923uXg2sBvqbWVegnbsvibabkbVPYpx1FlRU\nwOWXxx2JiEh+5ZQQzKyFmS0F3gKejn6pd3H39QDu/hawT7R5d2Bt1u6vR2PdgXVZ4+uiscS5+WaY\nOzesjiqfVVFREXcI0gyav3TL9Qhhe1Qy6kH4a/9wwlHCZzbLd3Clqn37cBrq2LGwcWPc0YiI5Eej\nlr929/fMLAOcBqw3sy7uvj4qB/092ux1YL+s3XpEY3WN16qyspJevXoB0KFDB/r167fjr5eaOmec\nr1u2hCFDKhg/HkaPjj+eUnmdXYMuhXj0WvNXzq8zmQxVVVUAO35fNkeDp52a2V7Ax+6+ycx2A+YB\nk4GTgHfd/UYzuwLo6O4Toqbyr4AvE0pCTwMHu7ub2R+AS4AlwJPAre7+VC2fWZKnne5syxY46qiw\ntMXQxHVDCiOTyez4hyvJo/lLtuaedppLQjiS0DRuET0ecvfrzawT8DDhr/41wDB33xjtMxEYA3wM\njHf3+dH4MUAV0AaY4+7j6/jMRCQEgOefh/POg+XLYe+9445GRNKs4AkhDklKCBDOOKquDtcoWJOn\nQkSkebS4XQm47jr485/hoYfijiR+2TVoSR7NX7opIeRBmzZhmezx4+HNN+OORkSkaVQyyqOrr4Y/\n/hGeeEKlIxEpPpWMSsiPfwyvvw733x93JCIijaeEkEetWoXS0RVXwJo1cUcTD9Wgk03zl25KCHl2\n5JFw2WVw4YWwfXvc0YiI5E49hALYtg1OPDHcfvO73407GhFJC12HUKJWroQTToDFi+Hgg+OORkTS\nQE3lEnXIIfDv/w6VlfDJJ3FHUzyqQSeb5i/dlBAK6OKLYddddctNEUkGlYwK7K9/hf79IZOBww+P\nOxoRKWcqGZW43r3D0hZjx6ardCQiyaOEUATjxoXS0dSpcUdSeKpBJ5vmL92UEIqgRQu46y645hpY\nu7bBzUVEYqEeQhFdey0sWQKzZmmtIxHJP/UQEmTCBHjtNXjkkbgjERH5PCWEImrVCu65B77/fdiw\nIe5oCkM16GTT/KWbEkKRHX88/Mu/hLusiYiUEvUQYvDee+GahBkz4OST445GRMqFeggJtOeecMcd\n4XTUrVvjjkZEJFBCiMngwdCvX7horZyoBp1smr90U0KI0W23hSbz8uVxRyIioh5C7O65J1y0tngx\n7LJL3NGISJKph5BwY8ZA27bhaEFEJE5KCDEzC0cIP/1pedyHWTXoZNP8pZsSQgno0wcuvRQuughS\nUikTkRKkHkKJ+OgjOPZYmDgRRoyIOxoRSSLdU7mMvPACDB0KL78MnTvHHY2IJI2aymXky1+GYcPg\nssvijqTpVINONs1fuikhlJif/hQWLYIFC+KORETSRiWjEvTkk3DJJfDSS+GUVBGRXKiHUKZGjID9\n94cbb4w7EhFJCvUQytSUKVBVBUuXxh1J46gGnWyav3RTQihRXbrATTdBZWU4JVVEpNBUMiph7jBk\nCBx1VPmtiioi+aceQpl7882wTPbs2fClL8UdjYiUsoL3EMysh5ktNLM/m9lLZnZJNN7RzOab2Uoz\nm2dm7bP2mWhmq81shZkNzBr/opktN7NVZjalqUGnSbduoZ9QWQkffBB3NA1TDTrZNH/plksPYRvw\nA3c/HDge+K6ZHQpMABa4+yHAQmAigJkdBgwD+gKnA3eaWU3GmgqMcfc+QB8zG5TXn6ZMDR8Ohx4K\n11wTdyQiUs4aXTIys8eA26PHSe6+3sy6Ahl3P9TMJgDu7jdG288FrgHWAAvd/bBofHi0/0W1fIZK\nRjv5+9/hC1+Axx6D446LOxoRKUVFPe3UzHoB/YA/AF3cfT2Au78F7BNt1h1Ym7Xb69FYd2Bd1vi6\naExysM8+cPvtMGqU7sMsIoXRMtcNzWwP4LfAeHf/p5nt/Cd8Xv+kr6yspFevXgB06NCBfv36UVFR\nAXxa50zb63PPreDRR2HUqAzf+U788dT2OrsGXQrx6LXmr5xfZzIZqqqqAHb8vmyOnEpGZtYSmA3M\ndfdborEVQEVWyWiRu/etpWT0FDCJUDJa5O59o3GVjJrgnXfgyCPhoYfgq1+NO5rPy2QyO/7hSvJo\n/pKtWCWj+4BXapJBZBZQGT0fBTyeNT7czFqZWW/gIODFqKy0ycz6R03mkVn7SI46d4apU2H0aNiy\nJe5oPk+/TJJN85duDR4hmNkJwHPAS4SykANXAi8CDwP7Ef76H+buG6N9JgJjgI8JJab50fgxQBXQ\nBpjj7uPr+EwdITRg5Eho3173YhaRT+nCtJTasCGUjh54AE4+Oe5oPqWSQ7Jp/pJNi9ulVMeOcNdd\ncOGFsHlz3NGISDnQEULCjRkDu+4K06bFHYmIxE0lo5TbtClcsHb33TBwYMPbi0j5Usko5dq3h3vu\ngbFjYePGuKPRWjhJp/lLNyWEMjBgAJx9djjzaPv2uKMRkaRSyahMfPRRONvojDPgqqvijkZE4qAe\nguzwxhvhngn33QeDtI6sSOqohyA77Lsv/OY3oXRUXR1PDKpBJ5vmL92UEMrM174GEybAOedoVVQR\naRyVjMqQO4wYAW3bwr33gjX5AFJEkkQlI/kcs3Aq6osvhusTRERyoYRQpvbYA373O/jxj0NiKBbV\noJNN85duSghlrE+fcIRw7rnw9ttxRyMipU49hBS48kp44QWYNw9a5nyPPBFJGvUQpEHXXQctWuiC\nNRGpnxJCCuyyS7g+YeZMuPxy+OSTwn2WatDJpvlLNyWElNhrL1i8GJYuhTPPDDfYERHJph5Cymzb\nBj/8IcyZA7NmwaGHxh2RiOSLegjSKC1bwpQpMHFiuKr5ySfjjkhESoUSQkqNHg2PPw7jxsHkyeHq\n5nxQDTrZNH/ppoSQYscfH05HffRROP98eP/9uCMSkTiphyBs3Qr/9m+waFG4ZmH0aGjdOu6oRKSx\n1EOQZtttN5g+HR55JDSaDz4Ypk2DDz+MOzIRKSYlBNnhy18OZx81JzGoBp1smr90U0KQz9k5MRx4\nYLjK+dVX445MRApJPQRp0PLlMGMG/PrX0KMHfPObMHw47L133JGJSDbdU1mKZts2eOYZeOABmD0b\nTjgBvvENOOuskChEJF5KCBKLzZvDRW2zZ8PcubD//iE5dOuW4dvfrqCFipGJlMlkqKioiDsMaSKd\nZSSxaNculI1++UtYvx5uvTU0nydPhn33hTFj4LHHYMuWuCMVkVzpCEHy7i9/CUcOs2eHC9++8hU4\n/XQ444xw5pKIFIZKRlLS3nsPFiwIZaU5c8I1D2ecERJERUV4LSL5oYQgJaW+GrR7OGNp7tzwWLoU\nTjzx06OHAw8sbqzyeeohJJt6CJIYZnDUUTBhAjz7LPztb3DhhZ8mhj594Pvfh4ULwxlNIlJcOkKQ\nkrB9O/zpT+HMpcceg+rqcDrr0KEwcCC0bRt3hCKlTyUjKUtr14bEMHMm/O//wqmnwjnnwODB4Qwn\nEfk8lYykpORrLZz99oOLLw7lo9degyFDwn2he/SAYcPgd7+DDz7Iy0dJFq1llG4NJgQzu9fM1pvZ\n8qyxjmY238xWmtk8M2uf9d5EM1ttZivMbGDW+BfNbLmZrTKzKfn/UaRcde4Mo0aF01hfew2+/nW4\n/Xbo1i2MP/NMKDmJSPM0WDIysxOBfwIz3P0L0diNwDvu/jMzuwLo6O4TzOww4FfAl4AewALgYHd3\nM3sB+J67LzGzOcAt7j6vjs9UyUga9MYb8NBDUFUVrpwePRoqK8PRhUgaFbxk5O7PAxt2Gh4CTI+e\nTweGRs8HAw+6+zZ3rwZWA/3NrCvQzt2XRNvNyNpHpEn23RcuvRSWLYOHH4Y334R+/eC008JKrbqf\ng0jjNLWHsI+7rwdw97eAfaLx7sDarO1ej8a6A+uyxtdFY1Jm4qhBm8Gxx8Kdd4Zm9AUXwNSpod/w\ngx/A//1f0UNKLPUQ0q1lnr5P3us7lZWV9OrVC4AOHTrQr1+/HRfM1Pyj1Wu93vl127bQo0eGq6+G\n/far4O674ZhjMhx2GFxzTQUDBsBzz5VOvHqt1815nclkqKqqAtjx+7I5cjrt1Mx6Ak9k9RBWABXu\nvj4qBy1y975mNgFwd78x2u4pYBKwpmabaHw4cJK7X1TH56mHIHmzdWu4l0PNAnzf+15oRuv0VSk3\nxTrt1KJHjVlAZfR8FPB41vhwM2tlZr2Bg4AXo7LSJjPrb2YGjMzaR6SgdtstrL66bBncdRdkMtCz\nJ1x+eWhMi0iQy2mnvwb+B+hjZn8zs9HAZGCAma0ETo1e4+6vAA8DrwBzgO9k/an/XeBeYBWw2t2f\nyvcPI/GrOZwtRWbwta/Bb38blsv46CM44ggYNw5Wr447utJQyvMnhacrlSWvMglbHO0f/4DbbgsN\n6VNOCessHX103FHFJ2nzJ5+lpStE8mDzZrj7bvjFL8JRw5VXhqMJkSRRQhDJow8/DHeBu+EG6N0b\nrr023OBHJAm0lpGUlKTXoFu3Dg3oV1+FESPg/PPD/RqWLGl433KQ9PmT5lFCEKnFrruGxLBqVVhh\n9eyzw9dly+KOTKRwVDISycEHH4RTVidPDiWkn/wEDj887qhEPkslI5EiaNMGLrkkLINx3HHhjKQL\nL4R16xreVyQplBAkr8q9Bt22Lfzwh6GU1K1buCXoFVfAhp2Xf0yocp8/qZ8SgkgTtG8P118Py5fD\nu+/CIYfAz3+um/ZIsqmHIJIHr7wSrl1YuhSuuw7+9V9hl13ijkrSRtchiJSQ55+HH/0I/vlPuPHG\ncG8Ga/L/niKNo6aylJS016BPPBF+//twFtKll8KppybrGoa0z1/aKSGI5JlZuG7h5Zdh+HAYOhTO\nOy80okVKmUpGIgX2/vtwyy3wn/8J554LkyaFM5RE8k0lI5ES17YtTJwIK1fCHnuExfOuugo2bYo7\nMpHPUkKQvFINum6dO4dTU5cuDTfmOfjgsLpqKZ2qqvlLNyUEkSLbf3+4/35YuDDcve2QQ2DGDPjk\nk7gjk7RTD0EkZs8/H6523rQp9BfOOQda6E81aQJdhyBSBtxh7txwuuqWLXD11aEBrcQgjaGmspQU\n1aCbxgzOOAP+8Ae46abQWzjySHjooeKWkjR/6aaEIFJCzMINeRYvDqepTpkSEsP994e7uYkUkkpG\nIiXMHZ55JiSHP/0JLr4Yvv1t6NQp7sikFKlkJFLGzODrXw/9hXnzwtXOBx0U7s3wl7/EHZ2UGyUE\nySvVoAunpnT08svhYrfjjoNBg2DmTNi2LT+foflLNyUEkYTZd99wK8+1a2HkyNCA7tkznLK6Zk3c\n0UmSqYcgUgZefhmmTYMHHwxLY1xwQThttUOHuCOTYtJ1CCKyw4cfhn7DAw/AggUwcCCcf34oLbVt\nG3d0UmhqKktJUQ06Xq1bh+W2H30UqqthwAC4446wuurZZ4clMt59t+79NX/ppoQgUqY6doRx48KR\nwl//GhLCzJnQqxecckq4o9vSpbB9e9yRSqlQyUgkZbZsCdc2PP00zJ8PGzaEI4kBA8Id3w48ULf9\nTCr1EESkWaqrQ3J45pmw0N62bfCVr8AJJ4SvRx2l/kNSKCFISclkMlRUVMQdhjTRokUZDjiggt//\nPtwbevFiWLECeveGo4+Gfv3C44gjoGtXHUmUmuYmhJb5DEZEks0sXNPQs2c4Owngo4/glVdg2bLw\nmD07JImtW8NNfvr0CY+DDw73eujRA7p3hzZt4v1ZpPF0hCAiTbJhA6xeHR6rVoWv69aFC+beeAP2\n3DMkhx49oEsX2Gsv2Hvv8DX70akTtGsHLfXnabOpZCQiJWf7dnj77U8TxNtvwz/+UfvjnXdg8+Zw\nRLHnniE57LnnZx/t2oX7Ue+2W+hn5PI1+3nLlukobykhSElRDyHZ4po/d3j/fXjvvc8+Nm8OXzdt\nCmdHbd0atqvta33vuYdrNFq3Dokn+2tTxprzfQp506PE9RDM7DRgCuEaiHvd/cZixyCFs2zZMiWE\nBItr/sxg993Do1u3/H//bdvCVdwffggffPD557mObdjQ/O/TsmXdSSb7Udd4Xdvtvnvz/zsVNSGY\nWQvgduBU4A1giZk97u6vFjMOKZyNGzfGHYI0Q7nOX8uW4ZGPX5rN4Q4ff1x3wsh+1DVe89i8OZTi\nal7n5cjD3Yv2AI4D5ma9ngBcUct2HrdFixaVxPdrzH65bFvfNnW915jxSZMmNRhDMZTC/BVz7up7\nP2nzVwpz19j9CvX/Xl3v1bV99Luzyb+ji710RXdgbdbrddFYycn3mi5N/X6N2S+Xbevbpq73GjNe\nXV3dYAzFUArzV8y5q+/9pM1fKcxdY/cr1P97db1XqDWnitpUNrNzgEHuPi56fQHQ390v2Wk7dZRF\nRJrAE9RUfh3YP+t1j2jsM5rzA4mISNMUu2S0BDjIzHqaWStgODCryDGIiEgtinqE4O6fmNn3gPl8\netrpimLGICIitSvJC9NERKT4dIMcEREBlBBERCSSiPUFzawtcCfwIfCsu/865pCkEcysN3AVsKe7\nD4s7HsmdmQ0BzgTaAfe5+9MxhySNYGaHAuOBzsBCd59W7/ZJ6CFE1ytscPcnzexBdx8ed0zSeGb2\nsBJCMplZB+Amd/9W3LFI45mZAdPdfWR928VSMjKze81svZkt32n8NDN71cxWmdkVWW/14NMrnD8p\nWqBSqybMn5SIZszdj4E7ihOl1KUp82dm3wBmA3Ma+v5x9RDuBwZlD2QtfDcIOBwYER3uQEgGPWo2\nLVaQUqfGzt+OzYoTntSj0XNnZpOBOe6+rJiBSq0aPX/u/oS7nwlc0NA3jyUhuPvzwIadhvsDq919\njbt/DDwIDInemwmca2Z3AE8UL1KpTWPnz8w6mdlUoJ+OHOLVhLm7mLA68blmNq6owcrnNGH+TjKz\nW8xsGvBkQ9+/lJrKtS181x/A3d8HLowjKMlZffP3LnBRHEFJTuqbu9uA2+IISnJW3/w9Czyb6zfS\naaciIgKUVkLIaeE7KVmav+TS3CVb3uYvzoRgfLbJqIXvkkXzl1yau2Qr2PzFddrpr4H/AfqY2d/M\nbLS7fwJcTFj47s/Ag1r4rjRp/pJLc5dshZ6/RFyYJiIihVdKPQQREYmREoKIiABKCCIiElFCEBER\nQAlBREQiSggiIgIoIYiISEQJQUREAPh/A/wjTQDtO94AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x21d7c44d780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Initial condition\n",
    "winit = np.random.randn(p)\n",
    "\n",
    "# Parameters\n",
    "feval = log_fun.feval\n",
    "nit = 1000\n",
    "lr = 1e-4\n",
    "\n",
    "# Run the gradient descent\n",
    "w, f0, hist = grad_opt_simp(feval, winit, lr=lr, nit=nit)\n",
    "\n",
    "# Plot the training loss\n",
    "t = np.arange(nit)\n",
    "plt.semilogx(t, hist['f'])\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can measure the accuracy of the final estimate by creating a predict method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy = 0.968198\n"
     ]
    }
   ],
   "source": [
    "def predict(X,w):\n",
    "    z = X.dot(w[1:]) + w[0]\n",
    "    yhat = (z > 0)\n",
    "    return yhat\n",
    "\n",
    "yhat = predict(Xts,w)\n",
    "acc = np.mean(yhat == yts)\n",
    "print(\"Test accuracy = %f\" % acc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not quite as good as the test accuracy for the `sklearn` method, which was about 98.5%.  The reason is that the learning rate was somewhat slow and we didn't yet fully converge.\n",
    "\n",
    "To see the effect of the learning rate, the code below tries different learning rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lr=    1.00e-05  Test accuracy = 0.558304\n",
      "lr=    1.00e-04  Test accuracy = 0.968198\n",
      "lr=    1.00e-03  Test accuracy = 0.989399\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x21d7ca6e208>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEECAYAAAAoDUMLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4VMXXgN8JCSGBhBR6R3pHQJqiEYWAAqIiAlIFUbGA\nqB9BEcEGogKKP3oXREGkCNIh2JCiIFVAkV4EQgklCSHz/TEJSTBls9ns3c2e93nus/fOnXvnbCZ7\nz51zZs5RWmsEQRAEwctqAQRBEATXQBSCIAiCAIhCEARBEBIRhSAIgiAAohAEQRCEREQhCIIgCICN\nCkEpdVgp9YdSartSaktiWbBSarVSar9SapVSqmCK+oOVUgeVUvuUUi1TlNdTSu1USh1QSo11/NcR\nBEEQ7MXWEUICEKa1vlNr3TCxLAJYq7WuAqwHBgMopaoDHYFqQGtgvFJKJV4zAeitta4MVFZKhTvo\newiCIAjZxFaFoNKo+wgwK3F/FtA+cb8d8JXWOl5rfRg4CDRUShUDArTWWxPrzU5xjSAIgmAxtioE\nDaxRSm1VSvVJLCuqtT4DoLU+DRRJLC8JHEtx7YnEspLA8RTlxxPLBEEQBBfA28Z6d2utTymlCgOr\nlVL7MUoiJQ6LgaGUkngagiAIdqC1VpnXShubRgha61OJn2eBxUBD4IxSqihAojno38TqJ4DSKS4v\nlViWXnl6bVq6vf322y5xv6xcZ0vdjOqkdy4r5Y7+u7lz/zmz73JT/7lC3zm7/7J6Lr362SVThaCU\n8ldKFUjczw+0BHYBS4GeidV6AEsS95cCnZRSeZVS5YGKwBZtzEqXlFINE53M3VNc43KEhYW5xP2y\ncp0tdTOqk965rJQfPnw4UxmcgSv0nzP7LqPz7tZ/rtB3Wb0up3576Z1z9N8oCZWZVkl8qC/CmIS8\ngbla65FKqRBgPuat/wjQUWt9MfGawUBv4AbQX2u9OrG8PjATyAd8r7Xun06b2hHaTnA+PXv2ZObM\nmVaLIdiJ9J97o5RCZ8NklKlCsAKllC5fXlO2LJQpw38+y5QBPz+rpRTSIjIyMsfeXoScR/rPvcm1\nCuHAAc3Ro3DkCLc+k/aPH4eCBdNWFmXLmi0kBJTdfxZBEAT3I9cqhIzkSkiAM2dSK4vblcaNG8mj\nibSURsmS4G3rHCvBZuQNM+coV64cR44csVoMwQUoW7Zsmv6e7CoEt3wkenlB8eJma9w47TqXL/9X\nWXz/ffLxmTNQrFjayqJMGShXDvLnd+rXEoQMOXLkiENmkgjuj8oh84dbjhAcwY0bxvR0u9JI+Vm+\nPDRokLzVrSu+C8E6Et/+rBZDcAHS+1/wSJORM7hxA3bvhm3bkrd9+6BSJbjrrmQlUasW+PpaKqrg\nIYhCEJIQheACxMTArl2plcTBg1C9euqRRI0a4ONjtbTWID6EnEMUgpBETikEt/QhWEW+fGZ0cNdd\nyWXXrsEffxjl8OOPMGYMHD5sRg4plUS1apAnj2WiC4IgZIqMEHKA6GjYsSN5FLF1K0RFQcuW0Lo1\nhIdDkSKZ30cQUuLKI4Ty5cszbdo0mjdvbrUoHkFOjRAkY1oOEBAAzZrBK6/A3Llw4ABs3w733w+L\nFkHlymaUMXQobNoEN29aLbEg5A6GDh1K7dq18fHx4Z133sm0/qBBgyhUqBCFCxcmIiIi1bkjR47Q\nvHlz8ufPT/Xq1Vm3bl22ZMuorXLlyuHv709gYCCBgYG0atUqW23ZiygEJ1G6NDzzDHz7Lfz7L3z0\nEcTGwrPPQtGi0LkzzJ5tzrkzkZGRVosguBg3nfjGU6lSJT766CPatGmTad1JkyaxdOlSdu3axc6d\nO/nuu++YPHnyrfOdO3emfv36REVF8d5779GhQwfOnz9vl1yZtaWUYvny5Vy+fJnLly+zcuVKu9rJ\nLqIQLCBvXggLgw8/hJ07jXnpgQdgyRIzemjQAN56C375RUYPgvsxfPhwnnjiCbp160ZQUBCzZs3K\n/CIH0a1bN8LDwylQoECmdWfPns2rr75K8eLFKV68OK+99tqtOE4HDhxg+/btDBs2DF9fXx577DFq\n167NwoULb10/ffp0qlevTmhoKK1bt+bo0aN2tZWEK5gDRSG4AKVKQZ8+sHAhnD0Ln3xipr0+95xZ\nKDdokJkC6w7IDCMBYOnSpXTs2JGLFy/y1FNPMW/ePIKDgwkJCSE4ODjVfkhICMePH8/8pg5mz549\n1KlT59ZxnTp12LNnDwB79+7ljjvuIH+K1akpzy9ZsoSRI0eyePFizp49S7NmzejcubNdbSXx1FNP\nUbRoUVq1asXOnTsd8h2ziigEF8PHB+67D0aONKOH1avNyuzWraFePTOL6cwZq6UUXBGlHLM5giZN\nmtC2bVsAfH196dy5MxcuXCAqKooLFy6k2o+KiqJUqVKOaTgLXLlyhYIFC946DgwM5MqVK2meSzof\nHR0NGBPQ4MGDqVy5Ml5eXkRERLBjxw6OHTtGWmTUFsCXX37J4cOHOXLkCGFhYYSHh3P58mWHfVdb\nEYXg4lSvDiNGmKmsH31kprhWqQIPPQRffQXXr1stYWrEh2AdWjtmcwSlS5fOvJLFFChQINVD99Kl\nS7dMTbefSzofEBAAGIdz//79CQkJISQkhNDQUJRSnDhxghEjRhAQEEBgYCD9+vXLtC0wCtTX15d8\n+fIRERFBUFAQP/74Y4599/QQheAm5Mlj/AwzZ8KJE9ClC8yYASVKQO/esHGj437MgpBdbo+18+WX\nX956SKbcksqsMBnVqFGDP/7449bxjh07qFGjxq1zhw4d4urVq7fO//HHH7fOly5dmkmTJhEVFXVr\npHPlyhUaN27M4MGDiY6O5vLly4wfPz7TttLCqinGohDckPz5oWtXWLXK+BaqVoV+/aBmTZg2zayo\ntgrxIQhp0aVLl1sPyZRbUpmjTEbx8fHExMSQkJDAjRs3iI2NJSEhIc263bt3Z/To0Zw8eZITJ04w\nevRoevXqBZjZSnXr1mX48OHExsby7bffsnv3bh5//HEAnnvuOT744AP27t0LmDf+b775Jl25Mmrr\n2LFj/PLLL7fk/eijjzh//jx33323Q/4mWcKRuUsdtRmxhKyQkKD1mjVat26tddGiWg8bpvWZM1ZL\nJTgSV/5dlC9fXq9bt05rrfWwYcN0t27dLJGjZ8+eWimlvby8bm2zZs3SWmv9448/6oCAgFT1Bw0a\npENCQnRoaKiOiIhIde7IkSM6LCxM+/n56apVq+r169enOj9nzhxdq1YtXbBgQV2mTBndu3fvDGVL\nr609e/bo2rVr6wIFCuhChQrpBx98UP/+++8Z3iu9/4XEcrufvbJSOReybx+MHQvz50OHDjBggImv\n5AwkllHO4corlQXnIiuVBZupVg0mTTIrpMuUMb6HVq3MjCV5ngiCkB4yQvAAYmJg3jwzSyk4GD74\nwExtFdwLGSEISUj4ayHb3LxpFMPbb0OFCvD++6kjtwqujSgEIQkxGQnZJk8eMztp3z547DF49FHz\neduCyWwh6xAEwX0RheCB5M1rwmIcPAh33w3Nm0O3bvDPP1ZLJgiClYhC8GD8/ODVV41iqFDBmI+G\nD8/e6meZYSQI7osoBIHAQBg2DH7/3aQIrVEDli2zWipBEJyNKAThFmXKwDffwMSJZuTQti0cOpS1\ne4gPQRDcF1EIwn9o2dJEWr37bmNGGjbMJPMRhPQoX74869evt1oMIZuIQhDSxNcXIiJM6s/t203S\nnt9/z/w68SEIVuKuKTST2LhxI15eXgwdOjRbbdmLKAQhQ8qUgcWLTZKeVq3MGoa4OKulEtwJSaGZ\neVtgAvMNGDCAxo0b29WGIxCFIGSKUmb9wo4dsG0bNGxo8jKkhfgQBEmhmbW2kvjkk08IDw+natWq\ndn13RyAKQbCZEiXM7KMBA+DBB81KZ8n5LKSFpNC0vS0w5qkZM2YwdOhQS1eje1vWsuCWKAU9e5qA\ned26wYYNMHcuFC1qzosPwTrUcMfkv9RvZ/+BlFYKzYwemFZgTwrNkydPAqlTaAJERETw/vvvc+zY\nsTSzxWWWQrN///689957+Pv7O+4L2oEoBMEuSpeGtWvhnXdMrufZs42SEKzDEQ9yR+EpKTRfffVV\nwOSVSUqhOWfOHD744AOUUnTt2pXx48dn2NZ3331HdHQ0HTp0yNHvawuiEAS78fY2CuHee81o4Zln\n4N57I3nggTCrRRMsJq0Ums8+++x/ypMepHv37nVY1jRbSUpr2aBBAyD9FJpJZqM//viDrl27Akbh\nDRkyJM1RT1IaTVvbWr9+Pb/99hvFixcHjLLw9vZm165dLFq0KAe+eQZkJ7tOTm24cGYoIW1OndK6\neXOt69bdoP/912ppcieu/LsoV66cS2RMu3Hjhr5+/bru0qWLHjJkiI6JidE3b95Ms+7EiRN19erV\n9YkTJ/Tx48d19erV9eTJk2+db9KkiX799dd1TEyMXrhwoQ4ODtbnzp3TWmu9aNEiXbNmTb1nzx6t\ntdYXL17UCxYsSFeujNq6cuWKPnPmzK3tySef1AMHDtQXLlxI937p/S+QzYxplj/80xQK9I2bN9L9\nYwiuSXy81oMHa12unNY7dlgtTe7DlRWCpNC0P4VmWt/hrbfeyvBeOaUQbM6HoJTyArYBx7XW7ZRS\nwcDXQFngMNBRa30pse5g4GkgHuivtV6dWF4PmAnkA77XWg9Ipy3dc3FPprWbhpeSiVDuxtdfw4sv\nwoQJJoWn4BgkH4KQhCvkQ+gP7E1xHAGs1VpXAdYDgxMFqg50BKoBrYHxKtlwOAHorbWuDFRWSoWn\n19iB8wd4ddWr8gNwMyIjI3nySZOu87XX4K23ICHBaqkEQbAFmxSCUqoU8BAwNUXxI0DSipNZQPvE\n/XbAV1rreK31YeAg0FApVQwI0FpvTaw3O8U1/2FZ52WsP7yed39419bvIrgQd94JW7bAxo3QsWP2\nQmoLguAcbB0hjAFeB1K+rhfVWp8B0FqfBooklpcEjqWodyKxrCSQcvXJ8cSyNAn2C2Z119V8sfML\nPtv8mY1iClaTch1CkSKwZo1JyPPgg2Dnqn9BEJxEptNOlVIPA2e01juUUmEZVHWobadnz56UK1eO\nh2MeZuiIoZxqc4oRfUYAyeERkh4+cuy6x76+0KdPJFOnQtOmYaxYAUePuo587nQsCLcTGRl5KwRG\nuXLlsn2/TJ3KSqkPgK4YB7EfEAAsAhoAYVrrM4nmoA1a62pKqQiMp/vDxOtXAm8DR5LqJJZ3Au7T\nWj+fRps6pVz7zu6j+ezmTHx4Io9UfSTbX1rIOSIjI9NdrTx+PLz3HixdaqKnCllDnMpCEpY5lbXW\nb2ity2it7wA6Aeu11t2A74CeidV6AEsS95cCnZRSeZVS5YGKwJZEs9IlpVTDRCdz9xTXZEi1wtVY\n1nkZz3z3DOv/kZjr7kq/fib5zkMPgbz0CoLrkZ05nSOBFkqp/cADicdorfcC8zEzkr4H+qV43X8B\nmAYcAA5qrVfa2lj9EvVZ8MQCOn3TiS0ntmRDbCEnySyWUbt2Zlpqx46wfLlzZBIEwTZsXofgTG43\nGaVk2YFl9Fnah7Xd11KzSE0nSyY4ii1bjHL49FN48kmrpXEPxGQkJOEK6xBcgjaV2zA6fDSt5rTi\n0IUsJvwVchxbHaANG5oZSAMHwvTpOSuTkPNICs3cgdspBIAutbrwZrM3afFFC05Fn7JaHMFOatUy\nvoS334Zp06yWRsgNuGsKzebNm1OkSBGCgoK48847Wbp0abbaspvsxL3IqQ0bY7a8/8P7uub4mvr8\ntfM21RdckwMHtC5VSuupU62WxLWx9XdhBSmD291OfHy80+SYPXu2XrlypW7fvr0ePnx4hnUnTpyo\nq1atqk+ePKlPnjypq1evridNmnTrfJMmTfRrr712K7hdUFDQreB2WSWztnbu3Knj4uK01lpv3rxZ\nBwQE6NOnT6d7v/T+F8itwe1sISEhQb+26jXdaEojHR0bbdM1gmuyf78ohcxwF4UwbNgw3aFDB921\na1ddsGBBPW3aNKfL07Vr10wVQtOmTfWUKVNuHU+fPl03adJEa631/v37db58+fSVK1dunb/33ntT\nPcSnTZumq1WrpkNCQnSrVq30kSNH7GrrdjZv3qz9/Pz01q1b071fTikEtzQZJaGUYlSLUdQqUov2\nX7UnJj7GapE8HnsXUVWuDOvWGfPR3LmOlUlwPpJC0/a2kmjbti1+fn40btyY+++//1buBGfi1goB\njFKY2GYiIX4hdF7YmfiEeKtFEuykcmVYtQpefVWmpNqFUo7ZHEBaKTQvXLhAVFQUFy5cSLUfFRXl\n9OQ4YF8KzejoaCB1Ck0vLy8iIiLYsWMHx44dIy0yS6EJJnPalStXWLFiBS1atHDId8wqbq8QAPJ4\n5WHOY3O4fuM6fZb2IUFLeE2ryG5O5Ro1YMkS6NULfvzRMTJ5DMYGnP3NAXhKCs2QkBBCQkIIDQ29\nlUJzxIgRBAQEEBgYSL9+/TJtKyV58uQhPDycVatWsWzZMod/58zIFQoBIG+evCzsuJCDUQcZuGpg\nki9CcEMaNTJmo8cfh+3brZZGsIe0UmgmPSRTbkllVpiMktJaJpFeCs0k/vjjj1vnS5cuzaRJk4iK\niro10rly5cqt9JnR0dFcvnyZ8ePHZ9pWWsTHx/P333879PvaQq5RCAD58+ZnWedlbDi8gXc2Zj7l\nTHA8jgrE1qKFiX308MNw8KBDbilYSJcuXW49JFNuSWWOMhnFx8cTExNDQkICN27cIDY2loR0EnJ0\n796d0aNHc/LkSU6cOMHo0aPp1asXAJUqVaJu3boMHz6c2NhYvv32W3bv3s3jjz8OwHPPPccHH3zA\n3r0mRcylS5f45ptv0pUro7b279/PypUriYmJIT4+njlz5vDjjz9y3333OeRvkiWy45HOqY1szqY4\nHX1aV/ysoh67aWy27iNknQ0bNjj0flOmmJScx4459LZuSXZ/FzmJpNC0P4Xmvn37dKNGjXRgYKAO\nDg7WDRs21EuWLMnwXun9L+CsFJrOJKPQFbZy5OIRms1oxrv3v0uPuj0cJJlgBaNGwaxZ8MMPEBpq\ntTTWIaErhCRyKnRFrlUIkBw2e8LDE2hfNd3kbIIbMGiQcTKvXw/58lktjTWIQhCSkFhGdpAUNrvv\nd31Zdyh7y84F28ipZC4jRkCZMmb2keRoFoScIVcrBEgRNnthJzYf32y1OIKdeHnBjBlw+DAMG2a1\nNIKQO8nVJqOULDuwjN5Le7Ou+zoJm+3G/PsvNG5slEL37lZL41zEZCQkISajbNKmchvGhI+RsNlu\nTpEisGwZvPaacTILguA4PEYhQOqw2SejT1otTq7EGQnhq1c3C9c6dpQ1CoLgSLytFsDZPH/X81yI\nuUDLL1qysedGQv09eB6jG9OiBbzzjlm49uuvEBJitUQ5T9myZf+zAljwTMqWLZsj9/UYH0JKtNb8\n35r/44ejP7C221oCfANyrC0hZxk4EHbvhhUrIE8eq6URBGuRdQh2orWm73d9OXTxEMu7LCeft4dO\nbndz4uMhPBzuugtGjrRaGkGwFnEq20lS2OxQv1AJm+1AnOFDSIm3N3z9NXz1Fcyf79SmcyXO7j/B\ntfBYhQASNju3UKgQfPstvPAC7NpltTSC4L54rMkoJVfjrhI+J5wGJRowJnyMOO7clC++gOHDYetW\nCA62WhpBcD7iQ3AQF2MuEjYzjEerPsrbYW87tW3BcQwYAAcOwHffiZNZ8DzEh+AggvIFsarrKubu\nmsunv35qtThui9U26I8+gmvXzJRUIetY3X+CtXjcOoSMKFqgKGu6raHZjGYE+wXTvY6HxUbIBfj4\nGAdz/frQrBk8+KDVEgmC+yAmozTYd3YfzWY044deP1C9cHXL5BDsZ/166NoVfvsNihe3WhpBcA5i\nMsoBqhWuxvvN36fH4h4yHdVNad4cnnsOOnc2axUEQcgcUQjp0Ld+X0L8Qhj5k6x2ygquZIN+801j\nQho+3GpJ3AdX6j/B+YhCSAelFNPaTeOzzZ+x4/QOq8UR7CBPHpgzB6ZPh9WrrZZGEFwf8SFkwuw/\nZvPxLx+z9Zmt+Hr7Wi2OYAeRkcZ09Pvv4k8QcjfiQ8hhutXuxh3Bd/DORpnH6K6EhcEzz0DPnpJ+\nUxAyQhRCJiilmNRmEtO2T5MUnDbgqjbot96CS5dg3DirJXFtXLX/BOcgCsEGihYoyrjW4+i+uDvX\nblyzWhzBDnx8TFKd996DnTutlkYQXBPxIWSBzgs7Uyx/Mca0GmO1KIKdzJoFH38MW7aAn5/V0giC\nY8lxH4JSylcptVkptV0ptUsp9XZiebBSarVSar9SapVSqmCKawYrpQ4qpfYppVqmKK+nlNqplDqg\nlBprr9BW8Xnrz5m/dz4bD2+0WhTBTrp3Nyk4IyKslkQQXI9MFYLWOha4X2t9J1AXaK2UaghEAGu1\n1lWA9cBgAKVUdaAjUA1oDYxXyeFDJwC9tdaVgcpKqXBHf6GcJNQ/lEltJtFrSS+iY6OtFsclcXUb\ntFIwcSIsWmSyrAmpcfX+E3IWm3wIWuskw7kvJv6RBh4BZiWWzwLaJ+63A77SWsdrrQ8DB4GGSqli\nQIDWemtivdkprnEb2lRuQ1i5MF5f87rVogh2EhwMs2dD797w779WSyMIroNNCkEp5aWU2g6cBtYk\nPtSLaq3PAGitTwNFEquXBI6luPxEYllJ4HiK8uOJZW7HmPAxrPhrBav+WmW1KC5HWFiY1SLYRFgY\n9OhhlIILuqssw136T8gZbB0hJCSajEph3vZrYEYJqao5WjhXpWC+gkxrN40+3/XhYsxFq8UR7GT4\ncDh1CqZMsVoSQXANshT+Wmt9WSkVCbQCziilimqtzySag5IG3yeA0ikuK5VYll55mvTs2ZNy5coB\nEBQURN26dW+9vSTZOa089sabR6o8Qv+V/ekV1MtyeVzlOKUN2hXkyex41ixo2jSSwEDo1Ml6eaw+\ndrf+8/TjyMhIZs6cCXDreZkdMp12qpQqBNzQWl9SSvkBq4CRwH1AlNb6Q6XUICBYax2R6FSeCzTC\nmITWAJW01lop9SvwMrAVWA58prVemUabLjnt9Hauxl2lzsQ6fNzyY9pXdTt3SI4QGRl56x/XXRg1\nClauhLVrwcvDV+a4Y/8JyeR4Ck2lVC2M09grcftaa/2+UioEmI956z8CdNRaX0y8ZjDQG7gB9Nda\nr04srw/MBPIB32ut+6fTplsoBICfjv7EEwueYOdzOymcv7DV4gh2cPMm3HOPyZ/wwgtWSyMI9iM5\nlV2A11e/zuFLh5nfYT7JM2wFd2L/fqMUNm2CihWtlkYQ7EOC27kA7zZ/lz3/7uHrPV9bLYrlpLRB\nuxNVqpj8Cb16eXYAPHftP8ExiEJwAPm88zH70dn0X9mfU9GnrBZHsJOXXzYL1z791GpJBMEaxGTk\nQIZuGMrvp37nu87fienITfn7b2jUCH7+2YwaBMGdEJORCzHk3iGciD7BjB0zrBZFsJMKFcz6hB49\njLNZEDwJUQgOJG+evMxuP5tBawdx5OIRq8WxhNxgg37+ecif30RF9TRyQ/8J9iMKwcHUKlqLV5u8\nytNLnyZBe7B30o3x8jJ5mD/+GPbssVoaQXAe4kPIAeIT4rln+j10q92NFxrKxHZ3ZfJkmDYNfvkF\n8uSxWhpByBzxIbgg3l7ezGo/i7cj3+bg+YNWiyPYyTPPGNORzDoSPAVRCDlElUJVeOvet+i5pCc3\nEzzHO5mbbNBKmcB3H3wAf/1ltTTOITf1n5B1RCHkIC81egkfLx/G/CopN92VChXgjTfMaMGTF6wJ\nnoH4EHKYfy78Q8OpDYnsEUmNIjWsFkewg5s3oWlTkzuhb1+rpRGE9JFYRm7AxG0TmfXHLH7q9RN5\nvMQ76Y7s2WOS6mzfDqVKWS2NIKSNOJXdgL71++Lj5cOEbROsFiXHya026Bo14KWX4LnncneGtdza\nf4JtiEJwAl7Ki8ltJzMschjHLh3L/ALBJYmIgKNHYd48qyURhJxBTEZO5J2N77D15FaWdloqsY7c\nlK1boW1b2LkTihTJvL4gOBMxGbkREfdEcOjCIRbsXWC1KIKd3HUXdO9uIqMKQm5DFIITyZsnL1Pb\nTmXAygFcuH7BanFyBE+wQQ8fDr/9BkuWWC2J4/GE/hPSRxSCk2lSugmPVXuM19e8brUogp34+ZmQ\nFi+8ABcvWi2NIDgO8SFYwOXYy9QYX4PZ7Wdzf/n7rRZHsJN+/SAuDqZOtVoSQTCID8ENCfQN5H8P\n/Y++y/py/cZ1q8UR7GTkSFizBtats1oSQXAMohAsol2VdtQtVpd3f3jXalEciifZoAMDYeJEE9bi\n6lWrpXEMntR/wn8RhWAh41qPY+rvU9l5ZqfVogh20ro13HMPvPmm1ZIIQvYRH4LFTP19KpN/m8ym\n3pskrIWbcv481KoF33xjYh4JglWID8HN6X1nb/x9/Bm3ZZzVogh2EhoK48bB00/DdXEJCW6MKASL\nUUoxue1k3vvhvVyRh9lTbdCPPw61a5s1Cu6Mp/afYBCF4AJUDq3MK41f4fnlz+MpprLcyOefw8yZ\nJryFILgjrutDuHoV/P2tFsVpxN2Mo8HkBgy+ZzCda3W2WhzBTubNg/ffNyuZfX2tlkbwNHKvD2Hw\nYKslcCp58+RlStspDFw9kPPXzlstjmAnnTqZLGsffGC1JIKQdVxXIXzzDWzYYLUUTqVRqUZ0rN6R\nV1e/arUoduPpNmilYMIEs+3YYbU0WcfT+8/TcV2FMHmymbYRHW21JE7lvebvseHwBtYeWmu1KIKd\nlCgBo0aZf98bN6yWRhBsx3V9CFpDnz7g5WWUgwex/MByXl75Mrue34W/j+f4UXITWsNDD8miNcG5\n5O6cypcvm7l8EydCq1ZWi+VUOi/sTJnAMnzY4kOrRRHs5NgxqFcPIiNNCk5ByGlyr1MZTLCYadPM\nSOFC7swfkB5jw8cy84+ZbD+13WpRsoTYoJMpXdrMOHr6aYiPt1oa25D+82xcWyEAPPAAtG/vcSmq\nihYoykd+t6bHAAAgAElEQVQtPqLnkp7E3YyzWhzBTp55BgoUgLFjrZZEEDLHtU1GSVy9CnXqwEcf\nwaOPWieYk9Fa88hXj1CnaB3ebZ67oqJ6EocOQaNG8PPPULmy1dIIuZnc7UNIyc8/Q4cOJrt54cLW\nCGYBp6JPUXdSXZZ1XsZdJe+yWhzBTsaNg6+/hh9+MPMkBCEnyHEfglKqlFJqvVJqj1Jql1Lq5cTy\nYKXUaqXUfqXUKqVUwRTXDFZKHVRK7VNKtUxRXk8ptVMpdUAplbVB9N13Q9eu8PzzZgqHh1A8oDhj\nw8fSc0lPYuJjrBYnU8QGnTYvvGDWKHz+udWSZIz0n2djy7tKPDBQa10DaAK8oJSqCkQAa7XWVYD1\nwGAApVR1oCNQDWgNjFdKJWmsCUBvrXVloLJSKjxL0r77LuzdC199laXL3J1ONTtRtVBVhkUOs1oU\nwU68vGD6dPMvfPCg1dIIQtpk2WSklFoMfJ643ae1PqOUKgZEaq2rKqUiAK21/jCx/gpgGHAEWK+1\nrp5Y3inx+ufTaCP9fAjbtsHDD5tloMWLZ0l2d+bfq/9Se0JtFndaTONSja0WR7CTceNMvKMff4Q8\nkv5CcDBOnXaqlCoH1AV+BYpqrc8AaK1PA0USq5UEjqW47ERiWUngeIry44llWaNBA3j2WTN9w4NM\nR0XyF+Hzhz6nx+IekofZjXnhBciXD0aPtloSQfgv3rZWVEoVAL4B+mutryilbn8aO/Tp3LNnT8qV\nKwdAUFAQdevWJSwsDIDIZs3gyy8JmzkTevW6Zfe8dT6XHncI68DCfQvpMaYH/Rr2s1yetI5T2qBd\nQR5XO/bygr59I3n2WXjooTBq1HAt+aT/3Os4MjKSmTNnAtx6XmYLrXWmG0ZxrMQog6SyfZhRAkAx\nYF/ifgQwKEW9lUCjlHUSyzsBE9JpT2fKzp1aFyqk9eHDmdfNRZy7ek4X/7i4/uHwD1aLkiYbNmyw\nWgS3YMoUrevV0zouzmpJUiP9594kPjtteq6ntdlqMpoO7NVaf5qibCnQM3G/B7AkRXknpVRepVR5\noCKwRRuz0iWlVMNEJ3P3FNdknVq1YOBA6N0bEhLsvo27EeofyoSHJ9BrSS+uxl21Wpz/kPQWI2RM\n795QtCiMGGG1JKmR/vNsMnUqK6XuBn4AdmHMQhp4A9gCzAdKYxzGHbXWFxOvGQz0Bm5gRhWrE8vr\nAzOBfMD3Wuv+6bSpM5MLMPEA7r4bevSAfv0yr5+L6L6oOwV9CzLuIcnF7K6cOAF33gkrV5qYR4KQ\nXTxnYVp6/PmnCSn5669QsWLOCuZCXLh+gVoTavHFo19wf/n7rRbnFpGRkfKWmQXmzIEPPzST51wh\nw5r0n3uTu4Pb2ULVqjBkCPTqBTdvWi2N0wj2C2Zy28k8vfRpomM9K2dEbuKpp6BSJRg2zGpJBCE3\njBDA+BDuvx/atYNX3TfbmD30XtIbnzw+TGwz0WpRBDv5918T5X3xYmgsS0yEbCAmoySSIoj98ANU\nq5Yzgrkgl2IuUXtibaa0nULLCi0zv0BwSRYuhIgI2L7dREcVBHsQk1ESd9xh4gJ07+4+wecdQMF8\nBZnadip9lvbhYsxFq8VJNY9dsJ3HH4dmzayP8i7959nkHoUAZgVzcDCMHGm1JE6lRYUWPFr1Ubov\n6k6C9pwpuLmNzz6Dn34yUVEFwQpyj8koiWPHoH59WL0a6tZ1rGAuTNzNOO6fdT8PVXyIN++VJL7u\nym+/QevWsGULOGLhqeBZiMnodkqXNol0evSAOM/JNJY3T14WPLGA8dvGs+qvVVaLI9hJ/frw+usm\n0rsHWT4FFyH3KQQwfoRy5eCdd6yWxKmUCCjBvMfn0X1xdw5fPGyJDGKDzj6vvgp+fvD2285vW/rP\ns8mdCkEpmDQJpkwxY28P4t6y9xJxdwSPz39coqK6KV5eMHcufPEFLLE/uIsgZJnc50NIyfz55jXr\n99/NK5eHoLWm88LO+Pv4M63dNJLzEwnuxObN0LatyZ1QpYrV0gjugPgQMqJjR6hTx6xk9iCUUkxt\nN5UtJ7Yw5fcpVosj2EmjRvD++/DYYxAti9EFJ5C7FQLA//6XnKLKgyiQtwDfPvktQ9YPYcsJ55nN\nxAbtWJ55xsRv7NnTOUF9pf88m9yvEEJDYeJE84u6csVqaZxK5dDKTGk7hQ7zO3D26lmrxRHsZNw4\nOHcOXnvNakmE3E7u9iGkpGdP8PeH8eMde1834I11b7D5xGZWdV2Ft5fNSfKS+eknE1FWsIwLF0wX\n9O5t0oAIQlpILCNbuXjRRBCbNg1atHDsvV2cmwk3aTW3FfWK1ePDFh9m7eLTp6F4cbh+3SQDFizj\n6FFjPnrnHRPcVxBuR5zKthIUBFOnmlesS5eslsap5PHKw7zH57Hoz0W8vvp1biZkIUz4zp3m8/Rp\nm6qLDTrnKFMG1q6Ft94yM6pzAuk/z8ZzFAJAy5bw8MPwyitWS+J0CvkXYlPvTWw/vZ2Hv3yYC9cv\n2Hbhrl3m88SJnBNOsJkqVWDDBhPH8bPPrJZGyG14lkIAE9Zi40b47jurJXE6of6hrOy6kqqFqtJo\naiP+PPdn5hft2mUW+p08aVMbWc62tWwZLFqUtWs8nEqVzL/wxInQv79j80JJtjTPxvMUQoECMGOG\niYx6/rzV0jgdby9vxrYay+B7BnPvjHtZfmB5xhfs2gUNGtisELLMmjUgZoosU748/PIL7NkDbdrA\nWZlEJjgAz1MIAPfeC506wYsvWi2JZfS6sxdLOi2h77K+jPxpJGk68W/ehH37jKnNRoWQZRv0qVNm\nCo2QZYKCYMUKs/byzjuNKSm7iA/Bs/FMhQBmCej27bBggdWSWEaT0k3Y3GczC/ctpMu3Xbh241rq\nCn/9BcWKGcN1To0QTp40M8AEu/DxMek/pk83+Zlffx2uXrVaKsFd8VyF4OcHs2bBSy/BmTNWS2MZ\npQJL8UPPH8ibJy9VP6/KxG0TiY2PNSd37YJataBECZudylm2QcsIwSG0bAk7dpg/Z82asDwTS2B6\niA/Bs/GcdQjp8eabxhC7aJFxnnowm49vZvjG4ez+dzdvNHuDPkuO4a2VefV89FH40wYndFbQ2iwW\nrFgxeTaTkG3WrIF+/UxW2REjoF49qyUSnIWsQ8guQ4fCoUMwZ47VklhOo1KN+P6p71nwxAKW7l/K\nmiVjWON/mtgioTnjQ7h0CWJiZITgYFq0MO84jzxiHM6dOsHu3bZdKz4Ez0YUgq8vzJ5tspIcP261\nNC5BkmK4PzqUb7z3U2HmncTdiOHA4d8c29DJk1CypPgQcoC8ec0o4eBBk0m2RQuTmnP9ejMwE4S0\nEJNREu++a2L2rFzp8aYjwHgmCxeGS5fYGbWP4vXDeOwpb2IrlqNb7W50qtmJwvkLZ6+NdevgvfdM\nJNrr142HVMgRYmJM0p2PPzZ/5j59TJrOkBCrJRMciZiMHEVEBERF5VxMAHdj716oXBl8fKhdtDaF\nK9YmsuWXvHv/u2w+sZlK4yrx8JcPM3HbRI5ftnNkdfKkcVgHBckoIYfJl89EbdmzB8aMgV9/NT6G\nzp1NOAxHLm4T3Bc7Ql/mUnx8zKyj++4z4+vy5a2WyFqSZhglUaIEeU6fIfyBpwivGE50bDTLDy5n\n2YFlvLn+TcoULEPbym0pfq44z3Z4Fi9lw7vGqVNGIQQHGz9C4WyOOIRM8fKCBx4wW1SUGTUMGmQm\nkXXoAJUrR/Lii2F4yauiRyLdnpLq1c2vo1cv52QjcWV27TLRYZMoUSKVYznAN4BONTsx57E5nHnt\nDJ+1+ozY+FhG/jSSEp+UoPeS3iz+czFX4zKYFH/ypImkKiMESwgJMbOuf/vNWO2KF4exY6F0aRgw\nAH74AeLjrZZScCaiEG7nlVfMr+Cjjzzb+7Zz539GCOnNNPL28qZZ2WZ82OJDjnx6hJ+f/pnaRWvz\nv63/o/gnxWk1pxWf/vopB88fTH3h7SMEwTIqVTIzsA8dCmPtWtMl/fsbJdGzp5mVLQvecj/iVE6L\nv/822c0DA8201NatPc/RXKSIWcldsqQ5/vprWLgQ5s/P0m0ux15m7aG1rDi4gu//+h4/bz8eqvQQ\nrSu2pmWv98jz/gcmzeljj5n5kYJLceQILF0KS5bAli0m6ssjj5ifR7FiVksn3I44lXOCChXMxO1X\nXzUmpLvuMr8KF1SeOcKZM2aUVKJEclkGI4SU3D6PPdA3kMeqPcaUdlM4/spxFnZcSImAEoz8eSRH\n9m2iz7ah7L5xgrMnDqZ9Q8Gp3N5/Zcsas9LatSZBz1NPmclh1apBkyYmbMb+/dbIKjgeUQjp4eUF\nTzwBf/xhxtJvv22WfH77be73LyQ5lFOOimxUCBmhlKJOsTpE3BPBxh6RlL+ej4fv7cNRr2imrBtF\n5XGVGbByAOv/WU98ghivXY2gIDMr6auvzDvDO+/AsWPQvLlREG+8AVu3es57U25ETEa2orWJ3f/O\nO2ZS91tvweOPQ548VktmOwcOmKmkmTFmjDGbff55ctm1a8YLef26Y8xnFy+aFGCXL8PIkeioKHa8\n0pnlB5ez+M/FHL54mDaV29C+antaVmiJv49/9tsUcoSEBNi2zfgZFi2CK1egfXsT7eTee2V5iTMR\nk5GzUMoYTrdsgVGjYPRo8xb95ZfuMYn7t99M1NJNmzKve/uUUzAxh/z8HOf8PXXKeCwBgoNRFy9y\nZ/E7GXLvELb13cb2Z7dTv3h9Ptv8GcU/Kc5jXz/G3J1ziY6Ndkz7gsPw8oKGDU3cpD//NOalkiVh\n8GDjZ+je3SiKmBirJRUyQxRCVlHKOJk3bYJPP4Xx48101dmzXXuO3uDB0KgRfPJJ5nXTUghgfuWZ\nRD21ORZO0qI0SHPaaemCpXmp0Uus77GeQy8f4pEqjzBv9zxKjSlFxwUd+Xbft8TEyxPG0TgillHV\nqubfbcsWY3Ft1AjGjTPd3asXrF7t2j8VTyZThaCUmqaUOqOU2pmiLFgptVoptV8ptUopVTDFucFK\nqYNKqX1KqZYpyusppXYqpQ4opcY6/qs4GaXMArYffzS5DGfMML+E6dPhxg2rpUvNunUmgN/KlSb3\n4qFD6de9edOsUq5Z87/nHOBHuEXSlFPIdNppqH8oPer2YFmXZRx6+RAP3vEgn2/5nOKfFKfH4h6s\nO7SOBJ3L/TpuSqlS8MILJobS7t0mmc+QIebd4sUX4eefc79Lzp2wZYQwAwi/rSwCWKu1rgKsBwYD\nKKWqAx2BakBrYLxStwzOE4DeWuvKQGWl1O33dE+UgvvvN+mqpk+HefOMnX7yZIiLs1o64/uIiDAx\ng4KCTBCbsRno40OHzJTTwMD/nrNBIdgcTz9pURpkaWFaqH8ofev3ZX2P9ezpt4e6ResycPVAKn5W\nkXc3vsuxS8dsa19Ik5zMh1CihFnwtmWLUQTFikHfvmZSX5KDWrCWTBWC1von4PbXt0eAWYn7s4D2\nifvtgK+01vFa68PAQaChUqoYEKC13ppYb3aKa3IP995rgtHPnWuMphUrGpOSlcbThQvNW3/Hjub4\npZdMqO/03sjTMxeBZSOE9CgRUIJXmrzCjmd3MP+J+Zy6coq6k+rSak4rFuxZkJzoR3A5KlY0I4Xd\nu+Gbb8yspbp1oVUrk8QwVrrOEuz1IRTRWp8B0FqfBooklpcEUur5E4llJYGUEdCOJ5blTpo2Nclu\nv/nGfFasaMxKzjacxsebKbMjR3IrOE2JEsY5PmlS2tdkphAc6UOwY4SQFkopGpRowPiHx3PslWN0\nrd2VCdsmUGpMKQauGshfUX/ZfW9Pw9n5EJSC+vXN2sTjx6FbN5gwwZiaBgxwfE4mIWMcFdzO4XNE\ne/bsSbly5QAICgqibt26t4azSf+0bnH83XdETpoEEycSNm4cfPIJkfnyOaf9AwegVCkifXwgMjL5\nfLNmEBFB2MCBkDdv6ut37SKyevXU9ZPOlygBa9c6Rr59+wjr188c//EHXLhAmNagVLbu7+/jT6mo\nUgwtO5TSbUsz5fcp1H+jPtULVWdYz2G0qNCCHzb+kDN/bznO9vFTT0HJkpGcPAl79oQRFgalSkXy\n6KMQERFGnjyuJa/Vx5GRkcycORPg1vMyW2itM92AssDOFMf7gKKJ+8WAfYn7EcCgFPVWAo1S1kks\n7wRMyKA9netISNB66VKtK1fWOjxc6127cra9q1e1LlFC6y1b0j7/wANaz5793/LKlbXevTvta379\nVeu77nKMfBUqaL1/f/Jx/vxaX7rkmHvfxrW4a3rqb1N17Qm1dZVxVfS4zeP05ZjLOdKW4FhiYrT+\n4gutGzbUulw5rUeN0vrcOaulcl0Sn502PdfT2mw1GanELYmlQM/E/R7AkhTlnZRSeZVS5YGKwBZt\nzEqXlFINE53M3VNc4xkkrWPYvRseftjEH37uOWM8zQk++8yYru66K+3zr75qpqCmXAB47Zrx7KW3\neM1RPgStU69DAONHyKGIp34+fvSu15sdz+5gctvJRB6OpOzYsry++nVORjvIJyLkCL6+JpHP5s0m\nnNbu3cYC++yzJhuc4GAy0xjAl8BJIBY4CvQCgoG1wH5gNRCUov5g4C/MKKJlivL6wC6Mo/nTTNrM\nWTXqCkRFaT1woNahoVp/8IHW16457t7nz2tdqFDqN/DbSUjQunp1rdeuTS7bulXr2rXTvyYuTmsf\nH63j49OtsmHDhszlu3hR6wIFUpfVqqX1jh2ZX+sgDl84rF/+/mUdPDJYP7P0GX3g3AGnte3K2NR/\nFnPmjNZvvWX+xR9/XOvNm62WyHUgp0cIWusuWusSWmtfrXUZrfUMrfUFrfWDWusqWuuWWuuLKeqP\n0FpX1FpX01qvTlH+m9a6lta6kta6f3YVmdsTHGze0H/91az7r1rVTFl1RMiOkSNN9NCMwlQoBQMH\npl6olpFDGUwMguBg+Pff7MmXcoZREk7OiVA2qCyftv6UAy8doHiB4jSd3pQnv3mS7ae2O00GwT6K\nFDHTVP/5x0zs69jRzPxesULiKGWb7GiTnNrwhBHC7WzcqHX9+lo3aqT1zz/bf59jx7QOCdH6xInM\n616/rnXRolrv2WOOX3lF65EjM76mbl2tt22zXz6ttV63Tuv77ktd1q6d1osWZe++2eByzGX9yS+f\n6JKflNThX4TrjYc3WiaLkDXi4rSeO1frOnXMAHfBAq1v3rRaKmvAST4EIae5916zYueFF+DJJ832\nzz9Zv8+wYWa1z+1v4GmRLx/062fiMsF/s6SlRcmS2fcjuMAI4XYCfAMY2GQgf7/8N09Uf4JeS3rx\n4OwH+eXYL5bJJNiGjw906WLSd4wYYXJb1aljUnfIKuisIQrBlfDyMhOx9+83oSMaNDBJY2bONA/R\nzNi3z+RtGDTI9jaff94sXjtzJnOTEWTqWE6aEpchKdcgJOEiWdN8vX3pXa83f77wJ51rdqbLwi60\nntuarSe2Zn5xLsCm/nNRlIKHHjJW2I8+So4/+fXX7hF/0hUQheCK+Pub8Nr79pl4ScuXQ40a5rVn\n0CATJiOtsBhDhsDrr5u3bVspXNiMRoYPN8tDS2ayXtARM41ccIRwOz55fOhdrzcHXjpAu8rtePTr\nR2k3rx07Tu+wWjQhE5QyK543bTJKYcwYoxi++kpGDJkh+RDchfh4Y1JatcoEqfvzT7jvPvOfHx4O\n585Bhw4m54GfX9bu/eefZkTStKnJrJ4RkycbOaZOtf+7dOpk8jB27pxc9umnJgfDZ5/Zf98cJCY+\nhsm/TWbkTyNpWropw8OGU6NIDavFEmxAaxNhdehQ887zwQe5Nyuu5EPwFLy9zQN7+HAzKfvvv43h\ndMsWuOceCAszWd2yqgzAzHBq3TpzcxE4ZoSQlsnIxUYIt5PPOx8vN3qZv17+i8alGtN8dnOeXvI0\nxy8fz/xiwVKUMu9Mv/5qfiKvvWZcdj/9ZLVkrocoBHelUKFk/8LJkybw/NNP23+/adOMmSozMnEq\n22SDTstk5CI+hMzw9/HntaavceBFM121zsQ6DFoziAvXXV92W3BnH0JmKGWyuO3aBb17mwVvbdrA\nzp2ZX+spiELIDShl1hx4ZaM7ixQx8YgzI7sjBK3THyG4gUJIomC+grz/wPvsfG4nUdejqPJ5FT7+\n5WNJ2uMG5MkDPXuauRstW5rtqafMoNvTEYUgZI3ChY1pJ534xEkBuNIlOtooroCA1OU5GLoiJykZ\nWJIp7aYQ2TOSn47+RJXPqzD7j9ncTHDPaS2Z9l8uwtcXXn7ZhMCoWtVkduvXD06ftloy6xCFIGQN\nLy8oWtT+X01aowNwG5NRelQvXJ3FnRYz97G5TNw2kTsn3cmKgyuQyRGuT0CAsZb++adxwdWoYSbs\nXbpktWTORxSCkHUyMBtlaoNOy38ALu9UtpV7ytzDz0//zPCw4byy6hUemP2AW61hyM0+hMwoVMhE\nctm+3fx7V6pkjq3Mb+VsRCEIWSc7q5X//NNkP7md/PnN2gpXSDuaTZRSPFrtUXb3202nmp1o/3V7\nnljwBAfOH7BaNMEGypQx2XAjI81MpEqVzJwLZ+e3sgJRCELWyWCEkKENOi7OLCF95pn/nlMq14wS\nkvD28qZv/b4cfOkg9YrVo+m0pjy37DlORduw6twiPMmHkBnVq5tMuAsWwBdfmFnZ336buwPoiUIQ\nso4NqTTTZPJkqFLFLKhLCzf3I6SHv48/g5sNZv+L+ymQtwA1J9TkzXVvcinGA43UbkjjxiY4wJgx\nJspq0nFuRBSCkHXs8SFcuQLvv2+ij6VHLhsh3E6ofygft/yY7c9u5+SVk1QaV4nRm0a71FRVT/Yh\nZERSOIzffze5nvv0Mcfbc1m0dFEIQtaxZy3CmDHQvDnUrZt+nVw6QridMgXLMOORGazvsZ7Iw5Fu\nP1XVk/DyMhFX9u2Ddu1M8sPOneGvv6yWzDGIQhCyTlZ9CGfPmlhF77yT8X3tGSH8+KPbKpGaRWqy\ntPNS5j42l0m/TaLOxDos2LOABG1dBDbxIdhG3rxmzcLBgyYMWOPGJnCwLUGJXRlRCELWyeosoxEj\nTJiNChUyrpfRCEHr/07z0NosOV2wwHZZXJB7ytzDT71+YlSLUXy86WNqT6jN/D3zLVUMgm3kzw9v\nvmlWPefPb9YwDBhgn4vNFRCFIGSdoCCzUvnq1f+c+o8N+sgRmDXLtjhJGYWvGDXKJP5JycGDcOiQ\nCfbn5iileKjSQ/za+1c+avERozeNptaEWny9+2unmpLEh2AfoaHw8cewZ4+JQ1mrFjz3nH05rqxE\nFIKQdZSyzY+QkAAvvWSywBUtmvl90wtfoTVMmWJGAleuJJd//z3Uq5crFEISSilaV2rNpt6b+KTl\nJ4zdPJZaE2oxY/sMYuPTDhciuA7FixvFcOCAWejWoAH06GGW37gDohAE+0hHIaSyQQ8dah7wQ4bY\nds/0TEabNpk8ic2awZIlyeUrVpiEQf/8A5cvZ01+F0cpRauKrfjl6V/4rPVnzN87n/KflmfEjyOI\nuh6VY+2KD8ExFCoE771nAuZVqmTCbbdtC+vXu/Y6BlEIgn1kNkL48kuzLVxoPHC2kJ5T+YsvTGrR\np54y9wRjrvrlFzP3r25d2LYt69/BDVBK8eAdD7LiqRWs6rqKA1EHqPhZRV5e8TJ/R0l4TlcnKMi8\nDx0+bBTCiy+axIfTp7tmSAxRCIJ9lCyZpucsMjLSJO3p39+8zRcubPs90xohxMaabOldu5osaz//\nbGYtrV8Pd90FgYEmTGUuMhulR62itZjxyAx299uNv48/jac1JnxOOIv2LSI+wTFxFcSHkDP4+xsX\n2J49xqT0zTdQtiwMHuxaYbdFIQj2kd4I4exZeOwxE/zFlgxsKUlrhLBsmXmlKlMGChQwWdTnzzfm\notatTZ3bFcLp03DtWtbadiNKBJRg5IMjOfbKMbrX7s7oX0dTdmxZ3t7wNkcuHrFaPCEDlDL5F77/\nHjZuNNFcmjQxS3TmzbN+1CA5lQX7+PJLk86zTRvzoE7apk6FJ56AiIis3/Ovv0yuw5SvTO3bm5FB\nr17mePlys+L55EmzX6OGGY83aZKsoOrVM/cZOTLbX9Nd2P3vbiZum8hXu7+iZpGadK3dlQ7VOxCU\nL8hq0YRMiI01g+mpU83K506dzGK3xo2znvMquzmVRSEI9nHpEnz9tUl4c+VK8lamjFEG9mQwP3fO\nxDo6fz75uGJFOHrUmIYAbtwwUzn8/c2UVqWMl65YMdi61aQS7dfP1D1yJHtZ5NyQ2PhYVvy1gi92\nfsHaQ2tpWaElXWp2IbxiOP4+/laLJ2TCP//A3LlmtHDlilEOnToZN5ktPylRCIJLERkZaf9Mlfh4\nyJfPjKOvXDGzlM6eNb+QlLz0Ety8CePHJ5e1a2ccz6NGGcPsu+/C6NFw//12fxd358L1CyzYu4D5\ne+az9eRWmpdvzqNVH6VN5TaE+IWkeU22+k9wGFqb3M/z5sFXX5lJdu3ama1pU7PWIS1EIQguRbYf\nKAEB8PjjZgz94INmlXPFiqnrxMUZheDnl1z2/vvGUxcfb0YJY8YYD9706WZUMWmSWSmU3i8plxN1\nPYplB5ax6M9FrDu0jgYlGhBeIZyWFVpSp1gdvJQZSYlCcD20ht9+g+++M9uRI8Z91rat+YmEhibX\nFYUg5C6eftqEuHj66bRTbabH2rXQooVxOD/xhAkqU726mQn1+utmYduoUSauABjF4aHK4WrcVdb9\ns441f69h9aHVXLh+gRYVWtDijhbcU+YeKgRXQNlj8hOcwrFjZq7FsmUmlFelSvDAA0Y5hIeLQhAE\nY2IaMsSYiZL8BuHhxgS1f79Zy9C6tRk9XL5sfkERESbLOpgA9/fcY8bmHsbhi4dZ8/ca1v2zjp+O\n/kR8QjxNSzfl7tJ307R0U+oUqyP+BxclLs7M8l67FnbsgCVLRCEILoRLmRzmzDFhM7ZsMc7qt94y\nqwM4oPYAAAjySURBVJ737TOKYPRoM6JYtcpMDu/VCyZONCOJNWtg6VJjmvrgA/i//zPrJFatMooj\nf37jASxXznj74uJsX4DnwmzYsIE77ryDn4/9zM9Hf2bT8U3sO7eP8kHlubP4ndQtWpe6xepSs0hN\nihUoJiMJF0NMRoJL4VIKISEBzpxJNj1dv25iCPTvbxa6rV5tprVWq2b8D48+atZC/PuviWkcGwtR\nUcaxXaKEMdpGRJjPHj3MtR9/bBbfvfiiWS39ySfGFPV//2cisX75pfFfNG1qRit//23G9mvXmvUa\n69aZkcyOHXDnnWbqbIkSxqSVN6/ZtLZv1pYdpNV/cTfj2Ht2LztO72DH6R1sP72dfWf3cT3+OpVC\nKlE5tDKVQytTKaQSZQqWoVRgKUoGliSfdz6nyCwkIwpBELLDhg1m3ULBgnD8uBlFjBpljh95xIwA\nxo0zju7ffjMP8F69zBSQadPMSOPmTaNg5s0zHr7r143ZKl8+oxxOnjT7YMxVtWqZkUqrVmaF0hNP\nwOLF5nP9erOmIibG+EHatjUjnQEDTHTXhARj7tqzxyiOmjWNwsqb1xiT8+UzyqRgQeN09/bOMYVy\n4foFDkYd5OD5gxw4f4CDUQc5fvk4xy4f42T0SQJ9AykVWIpSgaUomr8ohfwLUdi/MIX8C6XaQvxC\nCPANwNvLM306jkQUgiA4g7g4s+YiNNSskzh61LzR//qrmcV0zz0mAVDnzmZUMW2aSQrUtq2ZPjtv\nnnko33EHzJhhZjwNHWpGK6NGwZNPmmiud91lTFre3mY7ejR5Bbefn1EUXl5GCXl7G5+Hl5dZmR0U\nZBRDXJxRCiVLQkiIURJ+fmYrUsSMaEJCTP3gYPMZFGSUSGCgWePh45MtJZKgEzh79ewtBXH26lnO\nXTtntuvnkvevneP8tfNEx0WTzzsfgb6BBOQNINA3MNUWkDeAAnkL4Ofjh7+PP37eiZ8ZHKfc9/by\n9gjzligEwaVwKZORK5Hy/zk2Fnx9zZt9aKh58y9TBnbvNuato0fNg79IEVi50sye+vprKF3aPNTn\nzTPhOs6dMz6NSpXMqOTvv43SCggwwf9iY8198uQxCsQGIoGwtE5Ur25mf5UubZRJqVJGiZQoYRSM\nv7/5TFJkWXz4aq25duMal2Mvp9qi46K5HHuZSzGXuHrjKtdvXOfajWtcj7/t88b1VPu3n9NofPP4\n4uvtSz7vfPjmSfz09k21n3TO19uXfHnSKEtjP9P73FYvaYpvTuB2CkEp1QoYi4mjNE1r/WEadUQh\nuCljx45lQNLUTiHn0do87PPkMUrg2jXzME5IMArjzBkzgrlxwyiNCxfMdvGiWW1++bJxoP/7L2B+\nmDnSe76+ZsRSuLBRKiVKmBwZRYuasiJFjHIsWNA47H19zXdy0Ft9fEI8sfGxxN6MJSY+5j/7MfEx\nxN6MTbWfab2b9t3H28s7XSWT8jOloknr/O318vvkp0vtLtlSCE412imlvIDPgQeAk8BWpdQSrbWb\npI8QMuNiVnMiC9lDqeT1FEnxpJIoWtTEesoCF4cNg7ffNgolLs4okZgYs3/tmjGHnT9v9s+eNQol\nKsqs97h2zYx2oqP/G/gwNtZkt3NUhrugIPNdQ0KMAilY0Iyykkxe/v7JZrJ8+fD288M7MJD8SY56\nHx+zeXtD3kQTmU8e8Es0w+XJk7x5eSVvSqXezyJaa24k3EhXYaT8TKlc0vqMjo3mbPzZW2WOGHk4\n24vTEDiotT4CoJT6CngEcDmF4GjTh733y8p1ttTNqE5657Ja7gq4Qv85s+8yOp/l/lPKPAiTHqiO\nRmuzJSSYkUtMjNmuXoUrV4xc5cqZkUxUlBnpnDtnFFBsrHmIx8WZh3J8vPlMSDAjnqgoM2KKjzef\nN28SefYsYcHBps24uOQR040bpk5SWdJ1ifuRsbGEJcXKSkhI3pKOk/5WXl7G1JY0oklni7x5kzAf\nH/IqRV6lCEh57sYNwnx9074uNpYwP7/UZTExhPn7p65boABTs9k1zlYIJYFjKY6PY5SEy+EKD5Ss\nXucKCuHw4cM2yZrTuEL/uaNCcEr/JT3AvLzMw/02pRO5eDFh7ds7rLnIYcMIGzbM8dclKQatTd0h\nQ5KVXRpb5IgRhA0alPa5UaMIe+21tM+NHk3YgAGpy8aOJezll1PXVcqst8kGTvUhKKUeB8K11n0T\nj7sCDbXWL99WTxwIgiAIduA2PgTgBFAmxXGpxLJUZOcLCYIgCPbh7GDxW4GKSqmySqm8QCdgqZNl\nEARBENLAqSMErfVNpdSLwGqSp53uc6YMgiAIQtq45MI0QRAEwfl4Vn5BQRAEIV1EIQiCIAiA82cZ\n2YVSyh8YD8QCG7XWX1oskpAFlFLlgTeBQK11R6vlEWxHKfUI8DAQAEzXWq+xWCQhCyilqgL9gVBg\nvdZ6Yob13cGHkLhe4YLWerlS6iutdSerZRKyjlJqvigE90QpFQR8pLV+xmpZhKyjTKjXWVrr7hnV\ns8RkpJSappQ6o5TaeVt5K6XUn0qpA0qpQSlOlSJ5hbNtYRuFHMOO/hNchGz03RDgf86RUkgPe/pP\nKdUWWAZ8n9n9rfIhzADCUxakCHwXDtQAOicOd8Aog1JJVZ0lpJAuWe2/W9WcI56QAVnuO6XUSOB7\nrfUOZwoqpEmW+09r/Z3W+mGga2Y3t0QhaK1/Ai7cVnwr8J3W+gaQFPgOYBHQQSn1P+A750kqpEVW\n+08pFaKUmgDUlZGDtdjRdy9hohN3UEr1daqwwn+wo//uU0p9qpSaCCzP7P6u5FRON/Cd1voa8LQV\nQgk2k1H/RQHPWyGUYBMZ9d04YJwVQgk2k1H/bQQ22nojmXYqCIIgAK6lEGwKfCe4LNJ/7ov0nXvj\nsP6zUiEoUjsZJfCdeyH9575I37k3OdZ/Vk07/RL4BaislDqqlOqltb4JvIQJfLcH+EoC37km0n/u\ni/Sde5PT/ecWC9MEQRCEnMeVfAiCIAiChYhCEARBEABRCIIgCEIiohAEQRAEQBSCIAiCkIgoBEEQ\nBAEQhSAIgiAkIgpBEARBAOD/AYqgqgLzbOKwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x21d7ca6e5c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Learning rate values to test\n",
    "lr_test = [1e-5,1e-4,1e-3]\n",
    "ntest = len(lr_test)\n",
    "\n",
    "# Strings for the legend\n",
    "leg_str = []\n",
    "\n",
    "for i, lr in enumerate(lr_test):\n",
    "    # Run the optimizer\n",
    "    w, f0, hist = grad_opt_simp(feval, winit, lr=lr, nit=nit)    \n",
    "    \n",
    "    # Plot the results\n",
    "    plt.semilogx(t, hist['f'])\n",
    "    leg_str.append(\"lr=%12.2e\" % lr)\n",
    "    \n",
    "    # Measure the test accuracy\n",
    "    yhat = predict(Xts,w)\n",
    "    acc = np.mean(yhat == yts)\n",
    "    print(\"lr=%12.2e  Test accuracy = %f\" % (lr, acc))\n",
    "    \n",
    "plt.grid()\n",
    "plt.legend(leg_str, loc='upper right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that increasing the learning rate, speeds the convergence time, but the optimization is beginning to go unstable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive Step Size\n",
    "\n",
    "The above example shows that gradient descent is sensitive to the step size.  We now consider a slight variant of gradient descent with an adaptive step-size using the Armijo rule in class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def grad_opt_adapt(feval, winit, nit=1000, lr_init=1e-3):\n",
    "    \"\"\"\n",
    "    Gradient descent optimization with adaptive step size\n",
    "    \n",
    "    feval:  A function that returns f, fgrad, the objective\n",
    "            function and its gradient\n",
    "    winit:  Initial estimate\n",
    "    nit:    Number of iterations\n",
    "    lr:     Initial learning rate\n",
    "    \"\"\"\n",
    "    \n",
    "    # Set initial point\n",
    "    w0 = winit\n",
    "    f0, fgrad0 = feval(w0)\n",
    "    lr = lr_init\n",
    "    \n",
    "    # Create history dictionary for tracking progress per iteration.\n",
    "    # This isn't necessary if you just want the final answer, but it \n",
    "    # is useful for debugging\n",
    "    hist = {'lr': [], 'w': [], 'f': []}\n",
    "\n",
    "    for it in range(nit):\n",
    "\n",
    "        # Take a gradient step\n",
    "        w1 = w0 - lr*fgrad0\n",
    "\n",
    "        # Evaluate the test point by computing the objective function, f1,\n",
    "        # at the test point and the predicted decrease, df_est\n",
    "        f1, fgrad1 = feval(w1)\n",
    "        df_est = fgrad0.dot(w1-w0)\n",
    "        \n",
    "        # Check if test point passes the Armijo condition\n",
    "        alpha = 0.5\n",
    "        if (f1-f0 < alpha*df_est) and (f1 < f0):\n",
    "            # If descent is sufficient, accept the point and increase the\n",
    "            # learning rate\n",
    "            lr = lr*2\n",
    "            f0 = f1\n",
    "            fgrad0 = fgrad1\n",
    "            w0 = w1\n",
    "        else:\n",
    "            # Otherwise, decrease the learning rate\n",
    "            lr = lr/2            \n",
    "            \n",
    "        # Save history\n",
    "        hist['f'].append(f0)\n",
    "        hist['lr'].append(lr)\n",
    "        hist['w'].append(w0)\n",
    "\n",
    "    # Convert to numpy arrays\n",
    "    for elem in ('f', 'lr', 'w'):\n",
    "        hist[elem] = np.array(hist[elem])\n",
    "    return w0, f0, hist\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now run the new optimizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm8XeP1/98r880skTkhQQgxpFS0NV0zVaIoMTZqqpq1\nv6Kt5tKq6duixqpZq6pVJZWSIAchJEgEMaQyhwxknm5yb9bvj3W2s+/JOfeeebrr/Xrt1z772Xs/\n+7nnSc5nr/WsZz2iqjiO4zhOqdGi2A1wHMdxnES4QDmO4zgliQuU4ziOU5K4QDmO4zgliQuU4ziO\nU5K4QDmO4zglSV4FSkT6i8jLIvKhiLwvIhdHy0eLyAIReTe6HRm652oRmSkiH4nI4aHyPUVkuoh8\nKiK35bPdjuM4TvGRfM6DEpHeQG9VnSYiHYF3gBHAycBqVf1D3PU7A48DewP9gReBwaqqIvIWcJGq\nThGRscDtqvpC3hrvOI7jFJW8WlCqukhVp0U/rwE+AvpFT0uCW0YAT6hqnarOAWYCw6NC10lVp0Sv\nexQ4Lp9tdxzHcYpLwcagRGQgMAx4K1p0kYhME5H7RaRLtKwfMD9028JoWT9gQah8ATGhcxzHcSqQ\nVoV4SNS990/gUlVdIyJ3A9dFXXe/BX4PnJOjZ3nuJsdxnBJBVRN5y1Ii7xaUiLTCxOkxVX0GQFWX\namzw68/A8OjnhcCA0O39o2XJyhOiqiWxjR49uqTqTPfeVK5v6prGzic7l6g8H99lc+3bVK/1vi2N\nOkupb9Pt82wphIvvQWCGqt4eFETHlAKOBz6Ifn4WGCkibURkELADMFlVFwErRWS4iAhwJvBMAdqe\nFdXV1SVVZ7r3pnJ9U9c0dj7ZuXx8b7mmnPs21Wu9b0ujzlLq20z6PBvyHcW3L/Aq8D6g0e0XwKnY\neNRmYA5wvqoujt5zNXA2sAlzCY6Llu8FPAy0A8aq6qVJnqn5/Juc4jBq1CgefvjhYjfDyQPet5WL\niKBZuPjyOgalqq8DLROcer6Re24AbkhQ/g6wW+5a55QTw4YNK3YTnDzhfeskI68WVDFwC8pxHKc0\nyNaC8lRHjuM4TkniAuWUBZFIpNhNcPKE962TDBcox3EcpyTxMSjHcRwnL/gYlOM4jlORuEA5ZYGP\nU1Qu3rdOMlygHMdxnJLEx6Acx3GcvOBjUI7jOE5FUugl3y+Jlm8lIuNE5BMReSG0HpQv+e4kxMcp\nKhfvWycZ+bag6oArVHUo8G3gQhEZAlwFvKiqOwEvA1cDiMguwEnAzsBRwN3R7OUA9wBnq+qOwI4i\nckSe2+44juMUkYKOQYnIv4E7o9uBqro4uvRGRFWHiMhVgKrqTdHr/wvUAHOBl1V1l2j5yOj9FyR4\nho9BOY7jlABlMwYVWvL9TaBXsLxGdK2nntHLfMl3x3EcByiQQMUv+Y6tCxUmpybP3Lm5rM0pBXyc\nonLxvnWSkdf1oCDxku/AYhHpFXLxLYmW52TJ9113HcW55w6kc2fo2rUrw4YN+3q1x+A/gx+X13FA\nqbTHj3N3PG3atJJqjx9nfnzbbbcxbdo0Bg4cSC7I+xiUiDwKfKmqV4TKbgKWqepNInIlsJWqXhUN\nkvgrsA/mwhsPDFZVFZE3gUuAKcBzwB9VdYuFD0VEb75Zuf9+iESgT5+8/nmO4zhOErIdgyrWku+T\ngScxq2gucJKqrojek5Ml33/7W/jb30ykevTI25/oOI7jJKGkBaoYhKP4rrkGnn0WXn4ZuncvcsOc\nrIhEIl+7EZzKwvu2cimbKL5icN11cMQRcPjhsGJFsVvjOI7jpENFW1AAqnD55TBpEowfD507F7Fx\njuM4zQh38cWRaKKuKlx4IUyfDs8/Dx07FqlxjuM4zQh38aWACNx5JwwZAsccA+vWFbtFTrrEh5s7\nlYP3rZOMZiFQAC1awJ/+BAMGwHHHwYYNxW6R4ziO0xjNwsUXpq4OTjsN1q6Ff/0L2rQpYOMcx3Ga\nET4GFUcqyWI3bYKTT7axqSefhNatC9Q4x3GcZoSPQWVA69bwxBMmVKefblaVU9r4OEXl4n3rJKNZ\nChSYa++f/7T5UaNGQX19sVvkOI7jhGmWLr4w69bB974HAwfC/fdbMIXjOI6TPe7iy5L27WHMGJg5\n0+ZKVZheO47jlC15FSgReUBEFovI9FDZaBFZICLvRrcjQ+euFpGZIvKRiBweKt9TRKaLyKcicluu\n29mhAzz3HEydCpdd5iJVivg4ReXifeskI98W1EPAEQnK/6Cqe0a35wFEZGfgJGBn4CjgbhEJTMN7\ngLNVdUdgRxFJVGdWdO5sWSZefx2uvNJFynEcp9jkVaBUdSKwPMGpRD7JEcATqlqnqnOAmcDw6IKG\nnVR1SvS6R4Hj8tHerl1h3Dh44QX49a/z8QQnUzzbdeXifesko1hjUBeJyDQRuV9EukTL+gHzQ9cs\njJb1AxaEyhdEy/JCt27w4os2ife3v83XUxzHcZymyPuS7wm4G7guukrub4HfA+fk8gGjRo36esnh\nTJd8f+mlag48EObPj3DKKaWzpHJzPQ7KSqU9fpy742nTpnHZZZeVTHv8uMyWfBeR7YEFqlorItXA\n7sCjwSq4Tdy7LTBGVXdv7JyIXAWoqt4UPfc8MBpbcXeCqu4cLR8JHKiqFyR5Xlph5o2xcCEccABc\nfLEFTzjFI+KL2lUs3reVS6HCzJ8C6kVkB+A+bKn2x1O8VwiNOUXHlAKOBz6Ifn4WGCkibURkELAD\nMFlVFwErRWR4NGjiTOCZFJ+dFf362Wq8t98O99xTiCc6yfAfsMrF+9ZJRqouvs2qWici3wfuUNU7\nRGRqUzeJyONANdBdROZhFtFBIjIM2AzMAc4HUNUZIvIkMAPYBPwkZApdCDwMtAPGBpF/hWDbbeGl\nl6C62rJPnH12oZ7sOI7TvEnVxfcWcBvwS+AYVZ0tIh+o6q75bmC65NLFF2bmTDjoILjhBjjjjJxX\n7zSBu4EqF+/byiVbF1+qFtRZwI+B66PiNAh4LNOHliODB1sI+qGHmiV18snFbpHjOE5lk3YuPhHZ\nChigqtObvLgI5MuCCpg+HQ4/3Makvv/9vD3GcRyn7CmIBSUiEeDY6PXvAEtE5HVVvSLTB5cru+8O\nY8fCUUfZsh3f+16xW+Q4jlOZpBrF10VVV2FRd4+q6j7AoflrVmmz556WYPZHPzK3n5N/wvOhnMrC\n+9ZJRqoC1UpE+mC58v6Tx/aUDcOHw9NP24KHL79c7NY4juNUHqlG8f0AuAZ4XVUvEJHtgFtU9YR8\nNzBd8j0GFU8kAj/4gaVG2n//gj3WcRyn5Ml2DKrZL1iYC158EU49FZ59Fr71rYI+2nEcp2QpSCYJ\nEekvIk+LyJLo9pSI9M/0oZXGoYfCww/DscfCO+8UuzWViY9TVC7et04yUh2DeghLRdQ3uo2JljlR\nvvtd+POf4eij4b33it0ax3Gc8ifVMahpqjqsqbJSoBguvjD/+Adccom5/YYOLVozHMdxik6hksV+\nJSKni0jL6HY68FUKjUu05PtWIjJORD4RkRdC60EVbcn3XPKDH8Dvf2+TeT/5pNitcRzHKV9StaC2\nBe4Avg0o8AZwsarOb+K+/YA12Nyp3aNlNwFfqerNInIlsJWqXiUiuwB/BfYG+gMvAoOj60a9BVyk\nqlNEZCxwu6q+kOSZRbWgAh5+2Cypbt1yW29Vlc3DGj4c9tkHhg2Ddu1y+4xSxPO1VS7et5VLQTJJ\nqOpcLJNE+MGXYQlkG7tvYlTcwowADox+fgSIAFdF639CVeuAOSISLPk+l8RLvicUqFJh1CjLNrF+\nfW7rXb0a3n4bJk+GRx4xK22XXUys9tnHhGvwYGhRrLWSHcdxckTGYeYiMk9Vt0nhugYLForIMlXt\nFjq/TFW7icgdwCRVfTxafj8wFluw8AZVPTxavh/wc1U9Nv5Z0fMlYUEVinXrYOpUeOst2yZPhhUr\nYO+9G4pWz57FbqnjOM2NQmUzT/jsLO4Nk3M1ycWS7+VyPHmyHV9xRez8smXQpk01b70FNTURPv4Y\nevSo5jvfgQEDInzzm3DiiaXRfj/2Yz+unOOiLPme8MbMLaiPgGpVXRxdXXeCqu5ciku+VwqbN9t6\nVq+9BuPHW4Rhnz4WyHH44basffv2xW5l40R8nKJi8b6tXPIaxSciq0VkVYJtNTYfKqU20tDaehYY\nFf38Q2LLt5fcku+VQosWsNNOcM458Pe/w5Il8NBD0L27LcDYqxcccgjcdJO5CzdvLnaLHcdx8pzq\nKLzkO7AYs4j+DfwDGIBZRyep6oro9VcDZ2NLvl+qquOi5XvRcMn3Sxt5pltQabJ6teUUHDfOtuXL\nLTvGz39uUYKO4ziZ4Ln44nCByp45c+Df/zbr6rTT4NproVOnYrfKcZxyo1ATdZ1mxMCBcNll8MEH\nZk3tsgs89RQUU/eDwVin8vC+dZLhAuUkpUcPG6v661/h17+2PIOzZhW7VY7jNBfcxeekxMaNcOut\ncMstcPnl8LOfQdu2xW6V4ziljLv4nILQpg1ceaVlsXjzTdhjD5gwoditchynknGBctJi4EBbmPHG\nG+GHP4SLLy5MWLqPU1Qu3rdOMlygnLQRgeOOsyCK996z+VX19cVuleM4lYaPQTlZsXYtfO97sM02\n8OCD0LJlsVvkOE6p4GNQTlHp0AGeew4WLjSXX11dsVvkOE6l4ALlZE379jBmjKVQOuOM/IiUj1NU\nLt63TjJcoJycUFUFzzxjS32ceips2lTsFjmOU+4UbQxKROYAK4HNwCZVHS4iWwF/B7YF5mB5+lZG\nr78a+BFQRyhPX4J6fQyqiGzYACeeaGHpTzxhe8dxmiflPAa1GVt24xuqOjxadhXwoqruBLwMXA0Q\nXQ7+JGBn4Cjg7mhmc6fEaNfO0iLV18NJJ9kEX8dxnEwopkBJguePwJaBJ7o/Lvr56+XgVXUOMBMY\njlOStG0L//gHtGoFXbtC587Zb+3bR7YoGzgQHn7Ylwcpd3wMyklGNivqZosC40WkHviTqt4P9FLV\nxQCqukhEgoXK+wGTQvcujJY5JUqbNiZSq1fnpr7XXoP9929Y9uGHltT23nvhjjtsmXvHcSqHYgrU\nvqr6hYj0AMaJyCdsufx7RoNJzWnJ91I+FoF3381NfUcfveX5b38bbrghwrhxcOyx1Rx1FBx+eITe\nvUvj7/fj1I8DSqU9fpzZccks+Z5LRGQ0sAY4hzSWg1fVtxLU5UESzZBVq+D66+H+++Gb37TsFsce\n6wltHaeYlGWQhIi0F5GO0c8dgMOB90lzOfiCNtopKvFv2vF07mxL1i9YYBOG77kHBgyA3/3OQ95L\nnab61mm+FCtIohcwUUSmAm8CY6Jh4zcBh0XdfYcANwKo6gzgSWAGMBb4iZtJTiKqqmwe1ssvw8SJ\n8OqrMHw4TJtW7JY5jpMuJeHiyyXu4nPCqMIjj8DPfw7nn29JblNhxx19mXvHyZZsXXwuUE6z4PPP\nbZHFTz5p+tr6els5eN99YcQIOOEEW13YcZz0cIGKwwWqMolEIl9HChWCVavg+efh3/+GsWPh4IPh\n7LPhu9+15Uac3FHovnUKR1kGSThOqdO5s2XCePxxmDsXjjoKrr4aDjsMPvus2K1znOaBW1COkyJ1\ndXDbbbaa8IknQseOia9r2xZ2392CMwYOdIvLab64iy8OFygn38yaBU8/bQEYiVi7FqZOhbfeMhH7\n7nehTx8Tql12sYwYXbsWts2OUwxcoOJwgapMynGcQhXeew9eeAGWL7fgi6lTYdKkhkl0q6osTdO+\n+9p2wAFW1lwox751UiNbgSpmqiPHqWhEYNgw28Js3mxiFbBqlVlbr79u2TBGjoTDD7dxsHhatoTt\ntzdLbJddYNttoYWPJDsViltQjlNifPEFjB8PtbVbntu4Ef73P5gxw7bly+GQQ2DXXWNjXSLQu3dM\nyAYM8HEwpzi4iy8OFyinObFsmYXBz5kTK9u8GRYutGjDDz6wVY7btIHWrWNb27Y2t6tvXxgyBNq3\nt+VRWrWyNb26dIlt7dvb1rGjbW3axOpz4XMawwUqDheoysTHKTJDFdavt3yEwVZXZysfL10K8+fb\n5OXa2ti59eth5crYtn49rFsHa9bYtnGjbXV1JlJVVSZqVVWxrX37hvtgC8Qt2KqqYP78CLvtVv11\nPcEWPm7b1p7VqlVs36qVlbdp40JZqjSrMSgRORK4DZu/9UCQ3dypfKZNm+YClQEiJhKJGDw4u7pV\nTag2bDARC2/r1iX+vGlTTOA2bLD1wt59dxrr1lWzYUOsrkSf6+tjIlpXZ59ra608ELFAuMIiFliM\nwTWJ9o2Vhbc2bWJiGWzxohtf1rJldt9zc6ZsBEpEWgB3YklkPwemiMgzqvpxcVvmFIIVK1YUuwlO\nHCKxH+4uXTKvp6ZmBTU1md9fX28CtnFjTMDCQrZxowlZba1dF7+PL1u92qzL8PlAVIN6whZpWHTD\nWyCgIg1FKyxuYZFLVB4IbFNl8dZlsn2651q2jG2JjoOyfFmwZSNQ2BLvM1V1LoCIPIEtEV+yApUP\nt1Q2daZ7byrXN3VNY+eTnSsHd145922q15ZL37ZsCR062JYLct3G+np48cUI3/pW9ddiFi9u8YIX\nFtpp0yLsuGN1wnPh47VrYdasCH36VDcQ6Xirs64OFi2K0KVLddJrVq2K0Lp1NfX1NNjWr4/QooWV\n19XFylu2hBYt7J6wiGVLOQlUP2B+6HgBJlolSzn/iKV6faF+xOaEowBKgHLu21ITqErv25YtYdKk\nCEcckVmd770X4bTTUru3piZCTU3T1zZ1XbLzicpVTaRGj45w9dUNRa1nz5SanZSyCZIQkROAI1T1\nvOjx6cBwVb0k7rry+IMcx3GaAc0lSGIhsE3ouH+0rAHZfBmO4zhO6VBOc9CnADuIyLYi0gYYiS0F\n7ziO41QgZWNBqWq9iFwEjCMWZv5RkZvlOI7j5ImyGYNyHMdxmhfl5OJzHMdxmhFl4+LLFBFpD9wN\n1AKvqOrjRW6SkyNEZBDwS6Czqp5U7PY4uUNERgBHA52AB1V1fJGb5OQAERkCXAp0B15W1Xsbvb7S\nXXzRcPTlqvqciDyhqiOL3SYnt4jIky5QlYmIdAVuUdVzi90WJ3eIiACPqOqZjV1Xdi4+EXlARBaL\nyPS48iNF5GMR+VRErgyd6k9sgm89TsmSQd86ZUIWffsr4K7CtNJJl0z6VUSOAf4DjG2q/rITKOAh\n4IhwQShP3xHAUOCUqCkJJk79g0sL1UgnI9Lt268vK0zznCxIu29F5EZgrKpOK2RDnbRIu19VdYyq\nHg2c3lTlZSdQqjoRWB5X/HWePlXdBAR5+gCeBk4UkbuAMYVrqZMu6fatiHQTkXuAYW5ZlTYZ9O3F\nWGLoE0XkvII21kmZDPr1QBG5XUTuBZ5rqv5KCZJImqdPVdcBPypGo5yc0FjfLgMuKEajnJzQWN/e\nAdxRjEY5WdNYv74CvJJqRWVnQTmO4zjNg0oRqJTy9Dllifdt5eJ9W5nkrF/LVaCEhgPjnqevcvC+\nrVy8byuTvPVr2QmUiDwOvAHsKCLzROQsVa0HLsby9H0IPOF5+soP79vKxfu2Msl3v1b8RF3HcRyn\nPCk7C8pxHMdpHrhAOY7jOCWJC5TjOI5TkrhAOY7jOCWJC5TjOI5TkrhAOY7jOCWJC5TjOI5TkrhA\nOY7jOCWJC5TjOI5TkrhAOY7jOCWJC5TjOI5TkrhAOY7jOCWJC5TjOI5TkrhAOY7jOCWJC5TjOI5T\nkpSNQInIIBG5X0SeLHZbHMdxnPxTNgKlqrNV9Zxit8NxHMcpDAUXKBF5QEQWi8j0uPIjReRjEflU\nRK4sdLscx3Gc0qIYFtRDwBHhAhFpAdwZLR8KnCIiQ6LnzhCRP4hIn+DyQjbWcRzHKQ4FFyhVnQgs\njyseDsxU1bmqugl4AhgRvf4xVb0CqBWRe4BhbmE5juNUPq2K3YAo/YD5oeMFmGh9jaouAy5oqiIR\n0dw2zXEcx8kUVc3Y61U2QRLpMHr0aCZMmICqFnUbPXp0SdWZ7r2pXN/UNY2dT3YuUXk+vsvm2rep\nXut9Wxp1llLfptrnEyZMYPTo0Vn/lpeKBbUQ2CZ03D9alhE1NTXZticnVFdXl1Sd6d6byvVNXdPY\n+WTnEpXPmTOnybYUknLu21Sv9b4tjTpLqW9T7fPq6mqqq6u59tprU2pPMkS18B4xERkIjFHV3aLH\nLYFPgEOAL4DJwCmq+lEGdevo0aO//oKcymDUqFE8/PDDxW6Gkwe8byuPSCRCJBLh2muvRbNw8RVc\noETkcaAa6A4sBkar6kMichRwG+Z2fEBVb8ywfi2G6Dr5JRKJ+AtHheJ9W7mISHkJVL5xC8pxHKe4\nlK0FlW/cgqpM/C27cvG+rVyytaAqMoqvpqaGSCRS7GY4juM0SyKRSE6C1dyCchzHcfKCW1AJcAvK\ncRyneLgFlQS3oCoTH6eoXLxvKxe3oBzHcZyKpCIFyl18lYe/YVcu3reVh7v4kuAuPsdxnNLAXXxO\ns8At4srF+9ZJhguU4ziOU5KkLVAi0l5ErhGRP0ePB4vI93LftMzxMajKw8cpKhfv28qjaGNQIvJ3\n4B3gTFXdVUTaA2+o6rCsW5MDfAzKcRynNCjGGNT2qnozsAlAVdcBGTfAcVLBLeLKxfvWSUYmArVR\nRKoABRCR7YHanLbKcRzHafZk4uI7HPglsAswDtgXOEtVJ+S+eenjLj7HcZzSoCjrQYlId+BbmGvv\nTVX9MtMG5BoXKMdxnNKg4GNQIvKSqn6lqs+p6n9U9UsReSnTBuQDj+KrPLw/Kxfv28qj4FF8ItIO\naA9MwJZsD1SxM/C8qg7JujU5wC2oysQTilYu3reVS8FcfCJyKXAZ0BdYSEygVgF/VtU7M21ELnGB\nchzHKR7PPAPHHgsiRRiDEpGLVfWOTB+Yb1ygHMdxiocIrF0L7dsXYQxKVe8QkV1F5CQROTPYMm2A\n46SCj1NULt63lYfkaGZsq/QfLKOxMahdgLHAUcBE4NHcNCnpc0cARwOdgAdVdXw+n+c4juMUl0xc\nfO8DewBTVXUPEekF/EVVD8tHAxM8vytwi6qem+S8u/gcx3GKgCq0aFFEFx+wXlU3A3Ui0hlYAgxI\n9WYReUBEFovI9LjyI0XkYxH5VESubKSKXwF3ZdBux3EcJ48EtkGubIRMBOrtqBXzZyxp7LvApDTu\nfwg4IlwgIi2AO6PlQ4FTRGRI9NwZIvIHEekrIjcCY1V1WgbtdsoYH6eoXLxvK4fNm22fK4FKawxK\nRAS4QVVXAPeKyPNAZ1Wd3sStX6OqE0Vk27ji4cBMVZ0bfc4TwAjgY1V9DHhMRC4GDgE6i8gOqnpf\nOm13HMdx8ksgTIFQZUtaAqWqKiJjgd2ix3Ny0wz6AfNDxwsw0Qo/+w4gpfD2UaNGMXDgQAC6du3K\nsGHDvp4IGLyt+bEf+3HpHAeUSnv8OLPj2267DZjG7343kHbtyJpMgiQeAe5U1SkZP9QsqDGqunv0\n+ATgCFU9L3p8OjBcVS/JoG4dPXo01dXVX39pjuM4Tv7ZsAGqqmDMmAhvvx3h2muvLXiQxD7AJBH5\nTESmi8j78QEPGbAQ2CZ03D9alhE1NTUuThVG/Ju2Uzl435YfjzwSZIpouFVV2fljjqlm3ryarJ+T\n9jwo4gIcMkRouMjhFGCHqGX1BTASOCXTygOBcpFyHMfJPQ8/3NQVER56KJL1c9IWqCCQIVNE5HFs\nom93EZkHjFbVh6JBEOMwq+4BVf0om+c4lYW/bFQu3rdOMjJaD6qU8Ym6juM4+eWggyA1z2zhJ+o6\nTsHxcYrKxfvWSUZFCpQvWOg4jlNMIkBN1rVkEma+Goi/aSXwNvBTVZ2VdauywF18juM4+aWUXXy3\nAf8Pm1zbH/gZ8DjwBPBgpg3JJW5BOY7jFJMIxbKg3lPVPeLKpqnqsETnCo1bUJVJxJcFr1i8b8uP\nQllQmcyDWiciJwH/jB6fCGyIfnZlcL5m7VqorYVu3YrdEqcc+OIL6NkTWra044ULoV+/xu/5/HPo\n0ye9BfLWr4dZs2CXXWDVKntex46xNvTqBUuWQI8etu/dO1b/559D374N61q/PvV/44sW2d/YIo+j\n/59/bvsNG6BzZ1i+HLp0gRUroK7O/pYuXeyaujpYvLhhkldV+z5at4ZNm2zy7fvvW1l9PXz5JcyY\nkb/2N0BV09qA7YAxwJfA0ujnHYAqYL9068v1Bujo0aN1woQJ6hSP2lrVAw5Q3XVX1fXri90apxw4\n8EDVSCR23LGj6urVjd+zzz6qkyen95zjj7ef4TFjVK+6SvWWW2Ln9t1X9c03VQ87THXCBNVvf1t1\nyhQ79+WXqr17N6zrwQdVzz039WdXV6tOnJhee9MlJjPF3CYojFaTmMx/zzOZqDsLOCbJ6YnpS2Tu\nqampKXYTmjWqcNFF9pbWqxdceSXcfnuxW+WUOsuW2ds+wMaNsGYNrF4ds24SsXKlWUHpMCsaxrV8\nuT2jbdvYuTVrzPIP9sEGsG5d7HPAunVmQaXK+vVm2VQ+1dHt2qxqyWTJ9x7AucDA8P2q+qOsWuJU\nDHfeCZMmwRtvmAth2DA48kg46qjM6/Rxisol6NvVq02QILZfs6bxe9etsy0TVE0samtjZbW1sW3j\nxtgGDT8HbNq0ZVlj1NXZ5qRGJmNQzwCvAS8C9bltjlPujB8P119vAtWpk5U9+iiccgpMnWoWleMk\nIixQgTDlU6DAhChs0YQFKixUEBMo1diY1MaNJlKpsmlTetc3dzIRqPaq2tiS7EXHk8UWh08/hdNP\nh3/8AwYNipUfeCD86Edw1lnw3HPpDWgHeF9WLkHfrlqVvkCtXZu5QIlsaUFt3BgTpfA+OKdqgQKt\nor+cbkElIxLdsiOTWJL/iMh3s35yHvHlNgrP8uVwzDHw29/CAQdseX70aBtjuCOlJSed5kZtrf3Y\nB4KUiotv82Yb08nGxRdYSeF2hC2neAsqvA8+p2MRNR+BqiYX86AyEahLMZFaLyKrRGS1iKQ5TOlU\nEnV1MHJ7FJoZAAAgAElEQVSkjTOde27ia1q3hr/+FX7zGwtZTRefeF25RCKRpGNPjQlU4JpLV6CC\nkOrNm62Oxlx88WNQ0FCQ0nXZuYsvPdIWKFXtpKotVLVKVTtHjzvno3FOefD//p+9jf7+941ft/32\n8H//Z+NR6UQ+OZVPEIkXL1DxUXNhwtF16RAIUmAhJbKgAmEKn09mQZWSi6/SchSkLFAiMiS63zPR\nlr8mOqXMAw/YuNLf/x7zyzfGmWfCbruZqKWDu2wrlyCCD9KL4guEKVOBCqynQIBUtxSmRBZUWJDS\ntYjyLVD1FRa2lo4FdUV0//sE2//luF1OGRCJwNVXw5gxsNVWqd0jAvfcA//5jwlbIfnnP2H27Ozr\nmTcPHnss+3qSUV9vlmaqvPIKvPVW4nOzZ9vf3dizbr11y/KxY1PPFvDAA5alIODhh+Grr2L133Zb\nw+vffnvLNDmZuPjSFag33oA334z9iAfitGED/Pe/MH26lQfitHatidayZbDttjBtmp3/y1+gutrG\nW5cuNYHafns47ji4+GI47TQ4+mi47jr49a/h4IPt/8mhh9rfc/XV0KaNXdenD3ToYNGtXbsmXkY9\nna1169S+i7Ihm1m+pbjhmSQKwsSJqltvrfryy5nd//LLqgMGqK5aldr12fbn+vWqW22l+pOfZFWN\nqqrecINq376qmzdnX1ciPv7YZuN/9VVq1591lurllyc+d++9lqEhGXPn2rM2bmxYfvzxqrfemtrz\nd9yxYQaI3XZTHT/ePi9erNqihWp9fex8TY3qhRfGjidMmKDPPafavr1lclBVvflma9cvf5n8uVOm\n2DVnnZVaO6+8UvVXv1I9/HC777rrVIcMUd1/f9XTT1f9v/+z8muusf0ll9j+0ENt37LllhkTunZV\nHTq02FkbSnGbkJNMEhllhBKR74jIqSJyZrDlVDWzxKP48stbb8H3v29BDwcdlFkdBx1kb5TXXJPb\ntiXjqadgm23MFRkec8iEF16wnG2ffJKbtsXzwQe2/+ij1K7/7LPkluHcuTB/fvJ7v/jC9osXNyxf\ntGjLsmR89ZVZGQHhjBArVlgwQtgSSpT9YfVqy3EXtqS6ds2tBRV26R16aOxzYEXFj4MF++DfSyL3\n2caNPp6amGqKEsUnIo9hLr39gL2j2zezbolTFrzzDhx7rLlxDj88u7puuQWeeAImT2762mxfOO67\nz9wtu+6anWtxzRpzUf3gB/DSS1k1KSkffmj7VF1sjQnUvHmwYEEsci2eILFoIFQBixdbotSm2LzZ\nxChw6YEdBy6/QKhWroydX7my4XEwBtW3b0PXXu/e+RGoDRssDVf4uLY2JlDxQtVYEMSmTY0HcjjZ\nkYkF9U1gX1X9iapeHN0uyXXDnNJj2jTzrd93H3w3BzPhune3yL/zzstv6O3HH9sk4mOOgR/+EB55\nJPO6XnkF9t7bRPrFF3PXxjAffAB77JGaBbV+vVk7s2cnjuCaO9d+YJcuTXx/IFCLFjUsX7QoNYEK\nLKTAgqqtNcEIW1DhPSS2oFatsszl4TGoPn0aF6i1ay1PXyYC1bXrliHlmQpUNpksnMbJRKA+AHrn\nuiFNISJDROQeEXlSRH5c6Oc3dz74wHLp3XUXjBiRu3pPPdUGiOMH0uPJZh7U/ffDqFE2gHziiSYy\nyX6wm+KFF8xyPPhgG+jPR9TUhx+ahZaKBTVnDmy3nQ2Qh91sAXPn2lIQydx8geUUtqCCBKmpuPgC\nyyl4diBMTQlU2IIK5kHFu/hSsaC23jp1gQgStQYCFbaeGhOoplzCbkHlj0wEamtghoi8ICLPBluu\nGxaPqn6sqhcAJwPfyffznBgffWQ/yrfeCieckNu6g6i+m26KZZnOJbW1lgvwnHPsuFMn+N734G9/\ny6y+cePgiCPs7b5vX3j33dy1FextfdYsiwhLxYL67DOLIBs0aEs336ZNZgnts09ygfr8cxOCsEAt\nXmxTBlKxoOIFKl6oUrWgVq+2dZI2bbIw7FQsqHXrbM2mbFx8wfhTY2NQTWUfT+Y+dbInE4GqAY4D\nfkfDUPOUEJEHRGSxiEyPKz9SRD4WkU9FJGGuPxE5BvgPMDaDdjsZ8OmncNhhJiAjR+bnGdttBz//\nOVxwQfKJhpmOQT39tLnLtt8+Vpapm2/uXPsB3iO6ZvTBB+d+HOrTTy2kecgQs/KaykXXmEAtXGjW\n6aBBjQvUnns2FKhFi+z5S5Y0PfGzKQsqXqjAhCB+DGrVKltcr2NHE4ZUxqDWrk3Pgop38a1daxZw\nvAWVbkZ1J3+kJVAi0hKoUdVX4rc0qnkIOCKu3hbAndHyocApoYnBZ4jIH0Skj6qOUdWjgdPTabeT\nGZ99ZtFO110HZ5yR32ddfrn9MGZq2STjvvu2TL908MFmJQTRcqkybpyJdbAa6iGH5F6gPvjAAjla\ntoQdd7Txs8aYNSu5QM2da2I3YEDjLr699mo4BrV4sb00tGjRtPvqq69MSMIWVJs2WwpTfJBEIguq\nU6eYQOXDxbdhQ8zN17VrrE1NBUm4QBWPtARKVeuBzSLSJdMHqupEYHlc8XBgpqrOVdVNwBPAiOj1\nj6nqFcCOInK7iNwLFHiKZ/Nj7lz7Af7lLy0Teb5p3Rr+/Gf46U8Tj6VkMgY1c6aN5xx3XMPyli0t\n6/qjj6ZXX+DeC6iutmVFcrkA3YcfwtCh9nmXXZoeh2rMgpo3r2mBSmZB9e5tLremxqG++goGD45Z\nUsuXw8CBDV173bpt6eILZwkPxqA6dzaRWrOmMC6+lSvteYEFFQhWvAXlQRDFIxMX3xrg/air7o/B\nlmU7+gHh/0ILomVfE7XULlXVH6vqPVk+z2mE6dNtiYyf/hTOP79wzx0+3IIDrr46N/Xdd5+lVmrT\nZstzP/yhzeNKNcihrs6spcMOi5V17Woi8uabjd979NEWnh9w2GGxUPJ4AgsKYOedUxOo7bZLbkFt\ns01ygdq40X6Ud999yzGoXr1MoIJxqLvvtkCWvfe2oJNg/tpXX5mlF7agttuuoYtv4EALLhk1KpYa\nqHt3+Na3rL5588xq6dTJtrAFtXBh8qwJgYtv+fLUsixMnmwW1MaNMYHq0MHGkGbNsja0axezoNLJ\nsefkh0zWg/pXdCtZRo0axcCBAwHo2rUrw4YN+3oMI3gT9+PEx9dfH+Hmm+FPf6pm5MjCP/+wwyKc\neSZcdFE1u+2WeX17713NQw/BHXdEiEQSX9+zp50fNqzp+lq1qmbgQPjkkwiffBI7P2hQhEceSX7/\nuHERXngBDjmkmr32ghdeiPDyy/D669UMHbrl9W+9FeH73weoZpdd4A9/SN7++nr47LMICxfCoEHV\nzJ7d8PzcudCpk52fP3/L+7/4Arp0ifC//8HixdWowiuvRHj3XTjyyGp69YKXXopQWwsffVTN+vXw\n9tsRunaFLl2svunTI1EXnx1PnRqhXTtYvtyOZ82K0KEDTJ1azaZNMHZshKoqu99SB0VYsiTm4qur\nizBxIqxZU03v3nbeqI7uY8fr1sHSpcnPJzpetChCq1ZQVVXNihXQqlXsvE0OjrBhA3ToUM3atSAS\niY7DpVa/H98GTAMGstNOOZjMnk0aikw3YFtgeuj4W8DzoeOrgCszrNtTHWXA5s2WXqZvX9U33yxu\nW+64w9LLZJNK6E9/Uh0xovFrbrxR9bzzUqvv6qsTp90ZM0b1kEOS3zdxoqV+OflkO540yY5//OMt\nr121SrWqSrWuzo5nzFDdYYfkdc+ebemiVFXXrlVt27ZhSqHDDlMdO1a1tla1detYvQGvv666zz72\nuWtX1S+/tM/HHaf6z3+qnnOOfY+q1v5997W277mnanW1lZ90kurf/mb1r1+vetFFljKoZUvrv+HD\nVa+4wu7baSfVzz5THThQ9RvfiKXFeeop1e23V/30U+uzf/3L0iNt2tR4Op1TT1X9y1/SS8Gz9daq\nXbqovveeaqtW1qbw+R12sH3v3rbv3Nn2rVoVK2VQuW4TdMiQ0UqhUx2JyGAR+aeIzBCRWcGWbjXR\nLWAKsIOIbCsibYCRQMah657qKD1qa22c6fHHzV21zz7Fbc/551v2g7GhWM10xqBU4Y9/hEuamD5+\n8smWAimVScLPPZd4cvJ++1nqp2TuoEmTLNAkyJbx9tsWITd16pbXfvihuQxbtrTjHXYw11yyMa6Z\nM+0agPbtzeUYdtUFQRJt2phLLX4y7sKFFioPNt4T3BsegwpcfEuWxMLeP/4YvvzSPn/1ldXdrZu5\n2pYvt/uqqswqWrHCXHwQm//UpYuNNwVMnhz52sXXsaM9q02bprPjr11rf3c6rFhhbry2bc1t265d\nw/NBu4J9x44N923bpve85ks1Q4bUZF1LJmNQDwH3AHXAQcCjwF9SvVlEHgfewIIe5onIWWrBFxcD\n44APgSdUNcVMZFtSU1PjC9ylyNKl9gO6ciVMnGjjFcWmdWvL5v3Tn2aWYWLCBNs3lSdw4EAbPxk/\nvvHrFiywLZFwd+1qIhEeYwozaZKNdy1bZt/1O+/A2Wfboo3x41/Tp9t4UEDr1jae8+mniev+3/8s\nQCEgPA6lGguSgMTjUJ9/btkboKFABWNQvXrFBGrp0tg407p1iQVq2TLbunWz7yUQrKYEau3amIuv\nUydrVyAIjbFsWfoCVVdnIhMIU7zgBO3q1KnhPl6onKaI8NFHNVnXkolAVanqS4CoRd3VAEenerOq\nnqqqfVW1rapuo6oPRcv/q6o7qepgVb0xg3Y5afLBB/aje8ABtiRDhw7FblGM737XBvj/9Cc7Tsci\n/uMfbdkDkaavPeWUpkPb//tfi94LLJt4DjjAslPEo2pLPOy7L3zzmzBlillQBx9sghDvn48XKLBA\niWQTdsMWFDQUqCVLrD+DPk0kUAsXxgSqd2+znFRtHwRJBFF88ZN2v/rKro0XqOXLbemVrbaK5eQL\nRHL9ehO2Ll1sC9h662o2bjSx6dTJnh8IQ2MsXZrZv9l27WIClcyCiheoYJ/oeYmCcJym59ClQiYC\nVRudtzRTRC4Ske8DJfVe4S6+phk71n4or7sOrr8+NrenVBCxPH2/+U0sIiwVZs82S/D0FGfK/eAH\ntp5VYxmpk7n3Ag44AF59dcvyuXNtP3CgRShGIhZ1t+uu8I1vbOnmSyRQjYWaz5yZ3IKaN88EPqB/\n/8QCFe/iW7PGvvuOHWMuviDXXljwN20yqycQqO7dY1nNu3Uzgfr8cxP13qHEaPPmbWlBLVxoAiBi\n+y++SM1SWbo0fQsKYi4+2NKCCoSoKUsqTCm92JUO1ey0U03WtWTys3Qp0B64BNgLmzT7w6xb4hQE\nVUtZdM458Mwzqf+QF4PddrM5TL/5TepjUHfdBWedlfqPRu/eZt0ky3BeW2suwyOPTF7HAQfA669v\n6bKbNAm+/W374d17b8sAP3SovXHHC5Squf12261hHY0JVLyLb7vtYgIVjD8FpOriC9x7IjGBWrbM\nfqC7dzehA/u8cKEJVYcOiS2o2bOJRvzFnhkIVLhs+vRIAyH4/PPULKhMXHwtWzZuQSWznBpz8WUi\nkk5qpC1QqjpFVdcAy1T1LFU9QVWbmAlSWHwMKjEbN1rm8Iceiv14ljrXXWcTahcsaPratWtNBC68\nML1nnHKKLfuRiNdeM5HYeuvk9/foYT/0wYqrAeHvePhwe+P/ZnRhmniBWrDAfix79GhYRzIXX319\nLFFsQNiCSkWgwi6+QKCCAAkwoVq82Nrdo4dtO+9sVse229rYWPfuJmbdupkFFS9QW23VcMwnkQX1\n5ZcNLZZULSjV9K2XQJySWVCZjEG5BZWICB9/XJN1LZlE8X1bRGYAH0eP9xCRu7NuSQ4pBRffrFm2\nrHkuswxkw5IllvB1yRJ72w//eJUyvXrBz34GTz1V3eS1f/kL7L9/bFA+VY4/3gIl4tPvgLlCU1la\n5MADt3TzvfFGTKD69TN3WrxABX76RO49gJ12MrdgXV3D8nnzTDCqqmJl6QiUakMXXzAGFVhQYKKz\ncqVZND17xgRq663t8yefmEAF186bZz/4bdqYMM2aZRYU2L57d7umc2cTqcBluHJldQMhWLw49WCE\ndK2XQKBErJ1NjUHFC1NTAlVxS65nTDWDB9dkXUsmLr7bsJx5XwGo6nvAAVm3JIcUw4JatswCDc4/\n31LPfOc7cMMN9qPxu9+lN46SS1ThscfMdbT//vCvf6XmPiklLrvMot8ay9hQXw9/+ANcemn69W+1\nlX03Y8Y0LFe1slQEKn4cat06c83ttVes7NprbckSMEFo29Z+sCG5QFVVmXXz2WcNy//3v4YBEmAi\ntGiRWcpNCdSqVTbuGPwgJ7KgWrY04fnoIxOk8883C/yqq0xs4gXqf/+zPTR08YFlB9lzz5gFVV1t\nrubzzjMLKmyx1NfHjh991KYDnHyyrec1YoT1RzBm2r69ZQz5xjdM/Pfe2/7uvfaytgwdatcMGGDJ\njuPde+3aQU2NjcdCTIA6dIiNiQXtCp+/4grbt2hh/XjMMfYdtW5t///Drtp+DXLiGIFrtqoquRhv\ntVXi8jBBXwUimStrrlUrexFp184iXeMZNGjLsrZt7fs3IvTvX5N9Q9KdOAW8Fd1PDZW9l81krFxu\n9iflnw0bVCdMUP3FL1T33lu1Y0fVI49U/f3vVadPj00ynT5d9cwzVbfaSvXyy1Xnzi1I81TVnnXU\nUap77KH69tuFe24+uOKKCXrYYcnPP/WUTQrNdHLvgw+qfv/7DctmzFDt3z+1OhcsUO3ePTZR9pVX\nrD2NcdRRqk8/bZ9HjlR95JHk1z3zTMOyu+5KPMl44EDVmTO37PNNm2wybW2tHX/4oU1SDVi+XLVT\nJ9VrrlEdPTpWvuuuquefb1uYSy5R3W+/2Hf297+rDh6suvvudvzHP9qE2FNOid3zwx/aJM7HH4+V\nvf22Teo8/ng7fvVVu+bCCxN/FwHHHacqkn5/b7+9TS5WVe3ZM/acVavsuffcYxOeb7hBtV071Wuv\ntfKf/tT2551nk4hVVe+7z6454AD7fu67z/7mZAwdGpvIumFDeu1Oh/AE49/9zvYPPBB79n335e/Z\nW7alwBN1gfki8h1ARaS1iPwMyHjOUrmgamHZt95qb3A9etgSEaq2FMWXX1o48hVX2NtT4L7YbTdb\n2uG99+xta9gwyw83fXrjz8uGzZstd9pee9lE0ilTGr7JlyNHHmlRa6+9tuU5VbjxRnuzTyW0PBEj\nRliuvXD27meftZVzU6mzXz+zFoKAhlTG+MLjUO+/n9iCAhKmjImP4AsI3HxBHr6AVq3MdResoPv5\n5zH3HphVE6xFFY6669XL/t3Hj40lsqBmz25oQa1c2dAKCAIjwgES6UTLhenQwSyjdPs7bEElmg/V\ntm3MTdmmTWJXX+DGa9PGPrdu3fBzKjQ1CTnXlFqUbqpk0uwfAxdiyVwXAsOAn+SyUdmSKxffF1+Y\ni+HMM+0H6NhjbRb9j35kA9STJ5v77qCDmp5hPmCATT6dNcsG3Y880tw9EybkZr5AwCef2HjIX/9q\nLqdf/KIy/OKHHVbNr39tSUrjv69IxFxW2az0262bzQn7739jZc88Y32eKuH5UOHxp2TsuacJVG2t\nufB23jnxdYkEKpGLD0ygpk+3OuMDO8JuvnCABNgPfZ8+FugRjEGBjT0lEqitt7bgibBA1dXFBCnY\nBy4+SCxQJgDVSaPnkhEIVLpUVTV08QX/b4P/I4E4BULVvr39uIddf8G8p0DEwuKU6v+1ZHPq8kWh\nBSoSiVBTU5N1PZlE8X2pqqepai9V7amqpwNnZt2SHJJpkMTatTYofvnlNldl6FD7kfrOd+zNfdYs\nmzh64omxN8V06drV3vRnz7Z6LrjAIryefHLLgfB02LTJrIj99jN//WuvJf/BK1fOOMPe/INMEQE3\n3WTWbLb/CU84wcYRwQbqZ8ywsZJUCQIlVNOzoD76yMYtkr3kpGtBvfKKjcPEWxfxAhW2oMAE6uOP\nG1pQPXuaJdSzZ8NrA/ELC1R4n6pANZVaKBnt2+dmkm7wWSQmSuEtkSUVCFRYlNIVqEJTaEGsrq4u\njkAl4Yoc1VNQ6uvNCrr+evsh6t0bbr7Z/vM9+KC9IT71FPz4xw1XZM0FbdtaypsZM+BXv4Lbb7cf\norvvTn/9malT7e1/wgRz5110Ufma9MmIRCwL9ejRDa2oqVPtDf+007J/xnHHwfPPW+Tlc89Z1GM6\nudeCQIlZs+xHrKm0UYMG2WTXl19O7t4Dy90XFqhEIebhOl99NXGU5oABsXD98ByogD59rO6wBRV8\nTuTiC+/jhSl+DzGxCoeYV1VBixaRjF186RIvUOH+TSRMweewcMZbUJm4+ApNpq7vYpOrn7Gy+fM/\n+wzuvdfelnv0MHfdl1/aG/gXX5i76Je/NKumEG8dLVqYa+r1182dOG6c/chcd11sEbhkbNhgLrwj\njrDoteefTz/EutwYOdIiIl94wY5vvtmi/HKRxLNXL1vOfdy49N17YN9969Y25hhM0G0MERuTfPTR\nLSfohund2/o6yIU3b14sIWs8gwaZxZNIoMLZJOJdfMFzYEsXX3gfEG9Bdepk4yrpWlAi9ncEwmSC\n1bSLr3377AUqPAYVHCeyoNq2bdyCysTF56RGrgQqh6Mo2RMegwqHf2+3nbnA3njDROGDDxoGPhQ7\nEeS++8K//20umnnzzIVz8cVbLkQHls5njz3M1TN9uiUkLde3pFQIXLYtW1pY8DXX2MvG+PEWqpwr\nTjzRwvInTEgtvDyMiFlRd96Z+iTob3zDAmgas6BEGrr5krn3IBb+m8yCasrF17FjQ9dZIExNWVDB\nZN14YYoXqHDodqyu6q+tquB8KhZUpi6+4GWmMQsqLE7xFlQ4SCJenDwvn1HwMSgRWS0iqxJsq4G+\nTVZQQKqraxg/vprhw+2t9sEHzU3y7LPm2ggCH+L/g5YKQ4bYqqUffmj/Ifbe27IdvPuuuYQuusjG\nmW64Af7xj4ZjBs2BE0+0iLPjj7cXj7DLKFuOP97cusE8mnQ54ACz8NIRKGhcoKChQCULkAD7t9Cu\nXcMIvoCwQCVz8cX/WwoEKhCigHiBAvu+gu8ssFTCAtW1q4lPvPs5WOo9IMhq3hi5cvElsqDCQRKB\nCCVy8cWPPQVilYxCvkDmMvAqEwo+BqWqnVS1c4Ktk6oWOGiycX7+c9vffLONI4UDH8rJyujTx0Ro\n1iybhHjssfbDs369WX7HH1/sFhaOcFRmixY26fXTT5te8yld+vWzpciPOy6z+6ur7Ydtzz1Tu36v\nvczqCHLcJSM8DtWYBSVi46WJJlIGAlVfbxlF+vRpeL5//y1f2vr2jU1ADVNVZcISjhTs0aPx47CA\nNWxzZIuxqrAbMBGdO2f2YhK2EOOtsOA4EL9g69DBBLOqytoWuFbbt7fPwT7YklGo7C09e8bGJ+Nd\nrtDwpaHkyWYSVSluFGiibjGorbVVR5sj8Sskb95sk2PzwdKlNrE1U2bPTu/6OXOavubJJ21yqqrq\n0Uer/vvfya+dOzfxBNb6epusO2uWTVKNZ9Mm1fnzU29ffPn8+Q1X7U00KT1R2ZNPTtjivqYm4G7Y\noLpoUePXJGLpUtU1a+zzokUNJ8wGz12wwP6O4O8JysN7Vfs+582zlYhXr44dJ2PNGru2sWtywfLl\nqitXNmx/8PfNm5fdStXpQpYTdUWLbQvmGBHRSvubHGf6dAsQmTHDrKmnnrJpEOmy7bY2H++GG8xl\n7Dj5RERQ1Yz9VhUWjOw4lcngwebqra21EPNMpz30729L1CfKD+c4pUZFCpQvt1F5NPf+DJLGvvqq\njTHEZ+FOlf79be5fKQUINfe+rUSKlkmiHCiF5TYcJ9fstJMt4ZIsQCIV+ve3zPBuQTn5pNQySThO\nXvEXjtwIVL9+lqmklCwo71snGS5QjlMm7LSTjUMlmwOVCkE4u1tQTjlQVgIlIu1FZIqIpDnH3yl3\nfJzCBAqyd/FBaQmU962TjLISKOBK4O/FboRTeKZNm1bsJhSdIUNsnwsLqpRcfN63TjIKLlAi8oCI\nLBaR6XHlR4rIxyLyqYhcmeC+Q4EZwFLKKDmtkxtWrFhR7CYUnb59bR2xbDLr9+ljdcSnLiom3rdO\nMophQT0EHBEuEJEWwJ3R8qHAKSIyJHruDBG5FTgF2Ac4FTinoC3OkHy4LrKpM917U7m+qWsaO5/s\nXDm4fIrRtyK2oGKiEPNU29O6NVx5ZSSllF/et6VRZzr3pnptpn2bSZ9nQ8EFSlUnAsvjiocDM1V1\nrqpuAp4ARkSvf0xVL1fVs1X1CuCvwJ8L2ugMKed/6KleX6gfsTlz5jTZlkJSzn2b7x+xxs553+b3\n3koTqKKkOhKRbYExqrp79PgE4AhVPS96fDowXFXTTgUqIp7nyHEcp0TIJtVRSWUhzwXZfBmO4zhO\n6VAqUXwLgfAKNv2jZY7jOE4zpVgCJTSMxJsC7CAi24pIG2Ak8GxRWuY4juOUBMUIM38ceAPYUUTm\nichZqloPXAyMAz4EnlDVjwrdNsdxHKd0qLj1oBzHcZzKoFTGoBzHcRynARUXxRePiLQH7gZqgVdU\n9fEiN8nJESIyCPgl0FlVTyp2e5zcISIjgKOBTsCDqjq+yE1yckA0AcOlQHfgZVW9t9HrK93FF51T\ntVxVnxORJ1R1ZLHb5OQWEXnSBaoyEZGuwC2qem6x2+LkDhER4BFVPbOx68rOxZdBLr/+wPzo5/qC\nNdRJm0zzNDqlTxZ9+yvgrsK00kmXTPpVRI4B/gOMbar+shMo0szlh4lT/+DSQjXSyYh0+/brywrT\nPCcL0u5bEbkRGKuqnu68dEm7X1V1jKoeDZzeVOVlJ1Dp5vIDngZOFJG7gDGFa6mTLun2rYh0E5F7\ngGFuWZU2GfTtxcAh2P/d8wraWCdlMujXA0XkdhG5F3iuqforJUiiHzE3HsAC7EtCVdcBPypGo5yc\n0B9oRtoAAAK9SURBVFjfLgMuKEajnJzQWN/eAdxRjEY5WdNYv74CvJJqRWVnQTmO4zjNg0oRKM/l\nV7l431Yu3reVSc76tVwFynP5VS7et5WL921lkrd+LTuB8lx+lYv3beXifVuZ5LtfK36iruM4jlOe\nlJ0F5TiO4zQPXKAcx3GcksQFynEcxylJXKAcx3GcksQFynEcxylJXKAcx3GcksQFynEcxylJXKAc\nJweIyOroflsROSXHdV8ddzwxl/U7TqniAuU4uSGY8T4IODWdG0WkZROX/KLBg1T3S6d+xylXXKAc\nJ7fcAOwnIu+KyKUi0kJEbhaRt0RkmoicC1+vi/OqiDyDpYNBRJ4WkSki8r6InBMtuwGoitb3WLRs\ndfAwEbklev17InJSqO4JIvIPEfkouM9xyo1KWQ/KcUqFq4CfquqxAFFBWqGq+0QTZ74uIuOi134D\nGKqq86LHZ6nqChFpB0wRkadU9WoRuVBV9ww9Q6N1nwDsrqq7iUjP6D3BWjvDgF2ARdFnfkdV38jn\nH+44ucYtKMfJL4cDZ4rIVOAtoBswOHpuckicAC4TkWnAm9gSBYNpnH2BvwGo6hIgAuwdqvsLtWSb\n04CB2f8pjlNY3IJynPwiwMWqOr5BociBwNq444OBfVS1VkQmAO1CdaT6rIDa0Od6/P+6U4a4BeU4\nuSEQh9VAp1D5C8BPRKQVgIgMFpH2Ce7vAiyPitMQ4FuhcxuD++Oe9RpwcnScqwewPzA5B3+L45QE\n/lblOLkhiOKbDmyOuvQeVtXbRWQg8K6ICLAEOC7B/c8DPxaRD4FPgEmhc/cB00XkHVU9I3iWqj4t\nIt8C3gM2A/9PVZeIyM5J2uY4ZYWvB+U4juOUJO7icxzHcUoSFyjHcRynJHGBchzHcUoSFyjHcRyn\nJHGBchzHcUoSFyjHcRynJHGBchzHcUqS/w84XBqgH347GwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x21d7ccd6d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "winit = np.random.randn(p)\n",
    "nit = 1000\n",
    "w, f0, hist = grad_opt_adapt(feval, winit, nit=nit)\n",
    "\n",
    "t = np.arange(nit)\n",
    "plt.subplot(2,1,1)\n",
    "plt.semilogx(t, hist['f'])\n",
    "plt.grid()\n",
    "plt.ylabel('Loss')\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.loglog(t, hist['lr'])\n",
    "plt.grid()\n",
    "plt.ylabel('Learning rate')\n",
    "plt.xlabel('Iteration')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we measure the accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy = 0.989399\n"
     ]
    }
   ],
   "source": [
    "yhat = predict(Xts,w)\n",
    "acc = np.mean(yhat == yts)\n",
    "print(\"Test accuracy = %f\" % acc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In-Class Exercise\n",
    "\n",
    "Try to a build a simple optimizer to minimize:\n",
    "\n",
    "    f(w) = a[0] + a[1]*w + a[2]*w^2 + ... + a[d]*w^d\n",
    "    \n",
    "for the coefficients `a = [0,0.5,-2,0,1]`.\n",
    "\n",
    "* Plot the function f(w)\n",
    "* Can you see where the minima is?\n",
    "* Write a function that outputs f(w) and its gradient.\n",
    "* Run the optimizer on the function to see if it finds the minima.\n",
    "* Print the funciton value and number of iterations.\n",
    "* Bonus:  Instead of writing the function for a specific coefficient vector `a`, create a class that works for an arbitrary vector `a`.\n",
    "\n",
    "You may wish to use the `poly.polyval(w,a)` method to evaluate the polynomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy.polynomial.polynomial as poly\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}