notes-aaron-2017-03-10.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": "true"
   },
   "source": [
    "# Table of Contents\n",
    " <p><div class=\"lev1 toc-item\"><a href=\"#Introduction-&amp;-Set-up\" data-toc-modified-id=\"Introduction-&amp;-Set-up-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Introduction &amp; Set-up</a></div><div class=\"lev2 toc-item\"><a href=\"#Note-on-code-included\" data-toc-modified-id=\"Note-on-code-included-11\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Note on code included</a></div><div class=\"lev2 toc-item\"><a href=\"#Starting-the-tutorial\" data-toc-modified-id=\"Starting-the-tutorial-12\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Starting the tutorial</a></div><div class=\"lev2 toc-item\"><a href=\"#Re-cap-from-last-time\" data-toc-modified-id=\"Re-cap-from-last-time-13\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>Re-cap from last time</a></div><div class=\"lev2 toc-item\"><a href=\"#Model-assumptions\" data-toc-modified-id=\"Model-assumptions-14\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>Model assumptions</a></div><div class=\"lev2 toc-item\"><a href=\"#Create-fake-data\" data-toc-modified-id=\"Create-fake-data-15\"><span class=\"toc-item-num\">1.5&nbsp;&nbsp;</span>Create fake data</a></div><div class=\"lev2 toc-item\"><a href=\"#Writing-a-basic-Logistic-Regressor-class\" data-toc-modified-id=\"Writing-a-basic-Logistic-Regressor-class-16\"><span class=\"toc-item-num\">1.6&nbsp;&nbsp;</span>Writing a basic Logistic Regressor class</a></div><div class=\"lev3 toc-item\"><a href=\"#A-black-box\" data-toc-modified-id=\"A-black-box-161\"><span class=\"toc-item-num\">1.6.1&nbsp;&nbsp;</span>A black box</a></div><div class=\"lev3 toc-item\"><a href=\"#Logistic-Regressor-shell-class\" data-toc-modified-id=\"Logistic-Regressor-shell-class-162\"><span class=\"toc-item-num\">1.6.2&nbsp;&nbsp;</span>Logistic Regressor shell class</a></div><div class=\"lev2 toc-item\"><a href=\"#Test-run\" data-toc-modified-id=\"Test-run-17\"><span class=\"toc-item-num\">1.7&nbsp;&nbsp;</span>Test run</a></div><div class=\"lev2 toc-item\"><a href=\"#Visualizing-Progress\" data-toc-modified-id=\"Visualizing-Progress-18\"><span class=\"toc-item-num\">1.8&nbsp;&nbsp;</span>Visualizing Progress</a></div><div class=\"lev1 toc-item\"><a href=\"#Setting-alpha\" data-toc-modified-id=\"Setting-alpha-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Setting alpha</a></div><div class=\"lev2 toc-item\"><a href=\"#From-above/last-time\" data-toc-modified-id=\"From-above/last-time-21\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>From above/last time</a></div><div class=\"lev2 toc-item\"><a href=\"#Lipschitz-constant\" data-toc-modified-id=\"Lipschitz-constant-22\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Lipschitz constant</a></div><div class=\"lev2 toc-item\"><a href=\"#Adaptive-step-size\" data-toc-modified-id=\"Adaptive-step-size-23\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Adaptive step-size</a></div><div class=\"lev2 toc-item\"><a href=\"#Backtracking-line-search\" data-toc-modified-id=\"Backtracking-line-search-24\"><span class=\"toc-item-num\">2.4&nbsp;&nbsp;</span>Backtracking line-search</a></div><div class=\"lev2 toc-item\"><a href=\"#Comparison-of-the-three\" data-toc-modified-id=\"Comparison-of-the-three-25\"><span class=\"toc-item-num\">2.5&nbsp;&nbsp;</span>Comparison of the three</a></div><div class=\"lev1 toc-item\"><a href=\"#Coordinate-optimization\" data-toc-modified-id=\"Coordinate-optimization-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Coordinate optimization</a></div><div class=\"lev2 toc-item\"><a href=\"#Lipschitz-and-Uniform-Sampling\" data-toc-modified-id=\"Lipschitz-and-Uniform-Sampling-31\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Lipschitz and Uniform Sampling</a></div><div class=\"lev2 toc-item\"><a href=\"#So-which-is-better?\" data-toc-modified-id=\"So-which-is-better?-32\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>So which is better?</a></div><div class=\"lev1 toc-item\"><a href=\"#Exercises\" data-toc-modified-id=\"Exercises-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Exercises</a></div><div class=\"lev2 toc-item\"><a href=\"#Compare-run-time-on-larger-data-sets\" data-toc-modified-id=\"Compare-run-time-on-larger-data-sets-41\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span>Compare run time on larger data sets</a></div><div class=\"lev2 toc-item\"><a href=\"#Use-cython-or-numba-to-speed-up-the-optimizations\" data-toc-modified-id=\"Use-cython-or-numba-to-speed-up-the-optimizations-42\"><span class=\"toc-item-num\">4.2&nbsp;&nbsp;</span>Use <code>cython</code> or <code>numba</code> to speed up the optimizations</a></div><div class=\"lev3 toc-item\"><a href=\"#numba-implementation\" data-toc-modified-id=\"numba-implementation-421\"><span class=\"toc-item-num\">4.2.1&nbsp;&nbsp;</span><code>numba</code> implementation</a></div><div class=\"lev3 toc-item\"><a href=\"#cython-implementation\" data-toc-modified-id=\"cython-implementation-422\"><span class=\"toc-item-num\">4.2.2&nbsp;&nbsp;</span><code>cython</code> implementation</a></div><div class=\"lev2 toc-item\"><a href=\"#UMass-Amherst-smoking-data-set-for-logistic-regression\" data-toc-modified-id=\"UMass-Amherst-smoking-data-set-for-logistic-regression-43\"><span class=\"toc-item-num\">4.3&nbsp;&nbsp;</span>UMass Amherst smoking data set for logistic regression</a></div><div class=\"lev1 toc-item\"><a href=\"#Stochastic-Gradient-Descent\" data-toc-modified-id=\"Stochastic-Gradient-Descent-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Stochastic Gradient Descent</a></div><div class=\"lev2 toc-item\"><a href=\"#Mini-batch-SGD\" data-toc-modified-id=\"Mini-batch-SGD-51\"><span class=\"toc-item-num\">5.1&nbsp;&nbsp;</span>Mini-batch SGD</a></div><div class=\"lev2 toc-item\"><a href=\"#Stochastic-Average-Gradient\" data-toc-modified-id=\"Stochastic-Average-Gradient-52\"><span class=\"toc-item-num\">5.2&nbsp;&nbsp;</span>Stochastic Average Gradient</a></div><div class=\"lev3 toc-item\"><a href=\"#References\" data-toc-modified-id=\"References-521\"><span class=\"toc-item-num\">5.2.1&nbsp;&nbsp;</span>References</a></div><div class=\"lev2 toc-item\"><a href=\"#SVRG\" data-toc-modified-id=\"SVRG-53\"><span class=\"toc-item-num\">5.3&nbsp;&nbsp;</span>SVRG</a></div><div class=\"lev3 toc-item\"><a href=\"#References\" data-toc-modified-id=\"References-531\"><span class=\"toc-item-num\">5.3.1&nbsp;&nbsp;</span>References</a></div><div class=\"lev1 toc-item\"><a href=\"#Further-reading\" data-toc-modified-id=\"Further-reading-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Further reading</a></div><div class=\"lev2 toc-item\"><a href=\"#Regularization\" data-toc-modified-id=\"Regularization-61\"><span class=\"toc-item-num\">6.1&nbsp;&nbsp;</span>Regularization</a></div><div class=\"lev3 toc-item\"><a href=\"#$L^2$-regularization\" data-toc-modified-id=\"$L^2$-regularization-611\"><span class=\"toc-item-num\">6.1.1&nbsp;&nbsp;</span><span class=\"MathJax_Preview\" style=\"color: inherit;\"><span class=\"MJXp-math\" id=\"MJXp-Span-1812\"><span class=\"MJXp-msubsup\" id=\"MJXp-Span-1813\"><span class=\"MJXp-mi MJXp-italic\" id=\"MJXp-Span-1814\" style=\"margin-right: 0.05em;\">L</span><span class=\"MJXp-mn MJXp-script\" id=\"MJXp-Span-1815\" style=\"vertical-align: 0.5em;\">2</span></span></span></span><script type=\"math/tex\" id=\"MathJax-Element-129\">L^2</script> regularization</a></div><div class=\"lev3 toc-item\"><a href=\"#$L^1$-regularization\" data-toc-modified-id=\"$L^1$-regularization-612\"><span class=\"toc-item-num\">6.1.2&nbsp;&nbsp;</span><span class=\"MathJax_Preview\" style=\"color: inherit;\"><span class=\"MJXp-math\" id=\"MJXp-Span-1816\"><span class=\"MJXp-msubsup\" id=\"MJXp-Span-1817\"><span class=\"MJXp-mi MJXp-italic\" id=\"MJXp-Span-1818\" style=\"margin-right: 0.05em;\">L</span><span class=\"MJXp-mn MJXp-script\" id=\"MJXp-Span-1819\" style=\"vertical-align: 0.5em;\">1</span></span></span></span><script type=\"math/tex\" id=\"MathJax-Element-130\">L^1</script> regularization</a></div><div class=\"lev3 toc-item\"><a href=\"#Block-sparsity\" data-toc-modified-id=\"Block-sparsity-613\"><span class=\"toc-item-num\">6.1.3&nbsp;&nbsp;</span>Block sparsity</a></div><div class=\"lev2 toc-item\"><a href=\"#Constrained-Optimization\" data-toc-modified-id=\"Constrained-Optimization-62\"><span class=\"toc-item-num\">6.2&nbsp;&nbsp;</span>Constrained Optimization</a></div><div class=\"lev3 toc-item\"><a href=\"#Projected-Gradient\" data-toc-modified-id=\"Projected-Gradient-621\"><span class=\"toc-item-num\">6.2.1&nbsp;&nbsp;</span>Projected Gradient</a></div><div class=\"lev3 toc-item\"><a href=\"#Proximal-Gradient\" data-toc-modified-id=\"Proximal-Gradient-622\"><span class=\"toc-item-num\">6.2.2&nbsp;&nbsp;</span>Proximal Gradient</a></div><div class=\"lev3 toc-item\"><a href=\"#Newton's-method\" data-toc-modified-id=\"Newton's-method-623\"><span class=\"toc-item-num\">6.2.3&nbsp;&nbsp;</span>Newton's method</a></div><div class=\"lev3 toc-item\"><a href=\"#References:\" data-toc-modified-id=\"References:-624\"><span class=\"toc-item-num\">6.2.4&nbsp;&nbsp;</span>References:</a></div><div class=\"lev2 toc-item\"><a href=\"#Dual-methods\" data-toc-modified-id=\"Dual-methods-63\"><span class=\"toc-item-num\">6.3&nbsp;&nbsp;</span>Dual methods</a></div><div class=\"lev3 toc-item\"><a href=\"#The-Fenchel-Dual,-geometric-multipliers-and-the-KKT-conditions\" data-toc-modified-id=\"The-Fenchel-Dual,-geometric-multipliers-and-the-KKT-conditions-631\"><span class=\"toc-item-num\">6.3.1&nbsp;&nbsp;</span>The Fenchel Dual, geometric multipliers and the KKT conditions</a></div><div class=\"lev3 toc-item\"><a href=\"#Dual-coordinate-ascent\" data-toc-modified-id=\"Dual-coordinate-ascent-632\"><span class=\"toc-item-num\">6.3.2&nbsp;&nbsp;</span>Dual coordinate ascent</a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\reals}{\\mathbb{R}}\n",
    "\\newcommand{\\ip}[1]{\\langle #1 \\rangle}\n",
    "\\DeclareMathOperator*{\\argmin}{\\arg\\min}\n",
    "\\newcommand{\\E}{\\mathbb{E}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(101)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction & Set-up\n",
    "\n",
    "## Note on code included\n",
    "\n",
    "For a package-style Python implementation of the code used in this file, you can `git clone` the repo found [here](https://github.com/asberk/simplML). \n",
    "\n",
    "## Starting the tutorial\n",
    "\n",
    "To start the tutorial, go to your preferred directory and enter\n",
    "> `git clone https://github.com/asberk/simplML.git`  \n",
    "> `cd simplML/doc`  \n",
    "> `jupyter notebook`\n",
    "\n",
    "Next, open the `Gradient descent (skeleton).ipynb` file (or `(master)`) file, respectively.\n",
    "\n",
    "## Re-cap from last time\n",
    "\n",
    "For a real-valued continuously differentiable function $f \\in \\mathcal{C}^1(\\reals^n)$, the gradient of $f$ — denoted $\\nabla f$ — gives the direction of \"steepest ascent\". This is the direction in which $f$ increases most quickly. This statement can be verified by optimizing the directional derivative of $f$ to see which direction $u \\in \\reals^n$ gives the largest d.d. \n",
    "\n",
    "The directional derivative of $f$ at a point $w^0 \\in \\reals^n$ in the direction $u \\in \\reals^n$ is given by $\\nabla f(w^0) u \\equiv \\ip{\\nabla f(w^0), u}$ and denotes the rate of change of $f$ in the direction $u$ when $\\|u\\|_2 = 1$. In this case, it follows from the Cauchy-Schwarz inequality that\n",
    "$$\n",
    "\\nabla f(w^0) u \n",
    "\\leq \\|\\nabla f(w^0)\\|_2 \\|u\\|_2 \n",
    "= \\|\\nabla f(w^0)\\|_2 \n",
    "$$ \n",
    "The two sides achieve equality precisely when $u = \\nabla f(w^0) / \\|\\nabla f(w^0)\\|_2$.\n",
    "\n",
    "Correspondingly, the direction $v \\in \\reals^n$ of steepest **descent** of $f$ at $w^0$ is *opposite* the direction of the gradient, *i.e.,* $v = -\\nabla f(w^0)/\\|\\nabla f(w^0)\\|_2$. Hence, if we're looking to minimize a function $f$, then a seemingly reasonable approach is to *march* in the direction opposite the gradient until [hopefully] we find the smallest value of $f$. \n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/f/ff/Gradient_descent.svg\" height=\"300px\" width=\"300px\" />\n",
    "\n",
    "## Model assumptions \n",
    "\n",
    "This is a quick re-cap from last time and for consistency of notation.\n",
    "\n",
    "The **data** is comprised of the **features**/**covariates** $X \\in \\reals^{n\\times d}$ and the **labels**/**response** $y \\in \\reals^n$. Suppose for the time being that $f(w; X,y): \\reals^d \\to \\reals$ defines a convex objective (\"cost\") function parametrized by the data. We think of $w$ as a vector of weights or **parameters** for a linear combination of the $d$ features. \n",
    "\n",
    "For example, the objective function for $L^2$ regularized least squares regression is given by \n",
    "$$\n",
    "\\argmin_{w\\in \\reals^d} f(X,y; w,\\lambda) = \\argmin_{w\\in\\reals^d}\\big\\{\\frac{1}{2} \\|Xw - y\\|_2^2 + \\frac{\\lambda}{2}\\|w\\|_2^2\\big\\}\n",
    "$$\n",
    "\n",
    "Since the function $f$ is assumed to be convex, a global minimizer $w^*$ exists. To find $w^*$ we need simply start at some point $w^0$ and *march* for a sufficiently long time in the direction opposite the gradient (*cf.* picture above):\n",
    "\\begin{align*}\n",
    "w^1 &:= w^0 - \\alpha_0 \\nabla f(w^0)\\\\\n",
    "w^2 &:= w^1 - \\alpha_1 \\nabla f(w^1)\\\\\n",
    "&\\,\\,\\vdots\\\\\n",
    "w^{j+1} &:= w^j - \\alpha_j \\nabla f(w^j)\n",
    "\\end{align*}\n",
    "It may not always be sensible while marching to take a step of size equal to the magnitude of the gradient. In general, one instead takes walks in the direction of the gradient scaled by some a parameter $\\alpha_t, t\\geq 0$. But how does one choose $\\alpha_t$? Should they all be the same? Is there an optimal $\\alpha$? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create fake data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Throughout this example, we'll be performing logistic regression on a two-dimensional data set, which will allow us to visualize the results. We'll be examining how long our code takes to fit and predict, by looking at both run-time and number of iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 500\n",
    "d = 2\n",
    "X = np.random.randn(n,d)\n",
    "y = np.zeros((n,1))\n",
    "for j in range(n):\n",
    "    if X[j,1] > X[j,0]:\n",
    "        if np.random.rand(1) >= .9:\n",
    "            y[j] = 1\n",
    "        else:\n",
    "            y[j] = -1\n",
    "    else:\n",
    "        if np.random.rand(1) < .1:\n",
    "            y[j] = -1\n",
    "        else:\n",
    "            y[j] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYVMcaB+DfoVcRVFBRFAtYgr23oNHYe0tibDGaWGI0\n0WhivK7daOxGg4oVFRUVGwoWUFREEFAQkN57r8uyu9/9g7iKy8ICS3Xe5+G5es6cmW8xdxjmzHzD\nEREYhmGYuk2ppgNgGIZhKo915gzDMPUA68wZhmHqAdaZMwzD1AOsM2cYhqkHWGfOMAxTDyisM+c4\nTonjOC+O424oqk6GYRhGPoocmf8MwF+B9TEMwzByUkhnznFcCwBjABxXRH0MwzBM+ShqZL4XwGoA\nbDspwzBMDah0Z85x3FgAiUTkA4D774thGIapRlxlc7NwHLcNwLcAhAA0AegCuEpEcz4qx0btDMMw\nFUBEZQ6SKz0yJ6I/iMiEiNoA+ArAw4878g/K1vqvDRs21HgMLE4WI4uTxfnuS15snTnDMEw9oKLI\nyojoEYBHiqyTYRiGKRsbmX/E0tKypkOQC4tTcepCjACLU9HqSpzyqvQLULkb4jiqrrYYhmHqC47j\nQNXxApRhGIapeawzZxiGqQdYZ84wDFMPsM6cYRimHmCdOcMwTD3AOnOGYZh6gHXmDMMw9QDrzBmG\nYeoB1pkzDMPUA6wzZxiGqQdYZ84wDFMPsM6cYRimHmCdOcMwTD3AOnOGYZh6gHXmDMMw9QDrzBmG\nYeoB1pkzDMPUA5U+A5TjOHUAjwGo/VefHRFtrGy9DMMwjPwUcmwcx3FaRJTHcZwygKcAlhPRi4/K\nsGPjGIZhyqlaj40jorz//qiOotE567UZhmGqkUI6c47jlDiO8waQAOAeEXkool6mbkhIyEFQUGpN\nh8EwnzRFjczFRNQdQAsAfTmO66SIepm6oWfPo+jc+TBCQ9NqOhSG+WRV+gXoh4goi+M4ZwCjAPh/\nfJ/H40n+bGlpCUtLS0U2z9SQzz4zREhIGho21KjpUBimznNxcYGLi0u5n6v0C1CO4xoDKCSiTI7j\nNAE4AthBRA4flWMvQBmGYcpJ3hegihiZNwNwmuM4JRRN21z8uCNnGIZhqpZClibK1RAbmTMMw5Rb\ntS5NZBimfmIDsLqDdeYMw5TI1TUSWlrbsGqVU02HwsiBdeYMU8dERGQgI4Nf5e0kJ+eBzxciJiar\nyttiKo/NmTNMHRIYmAILiyNo394A/v5Lq6W91q0bQkNDoauYmXKoztUsDMNUE11dNRgYaKJVq4bV\n0l6HDo2rpR2m8tjInGGYWiclJQ9eXvEYMaINOK7MQWm9xlazMAxTZ82da4+RI21w7pxvTYdSZ7Bp\nFoZhah1Ly1YICUmDhYVhTYdSZ7BpFoYpg79/MqZNu4TZs7vg998H13Q4zCeGTbMwtdq6dQ8wf/51\nCASimg6lTH5+SQgISMGDB+E1HQrDyMRG5ky1E4sJGhpbUFgoxtu3y2Bm1qimQyoVEcHJKRTdujWF\nkZFOTYfDfGLkHZmzzpxRqPPnfRETk4XVqweUugrh3r1QJCfn4ZtvLKoxOoape1hnztQIVdXNEArF\neP36R1hYGNV0OAxT57FNQ0yNOHhwNKKjM9G5c9WsQtix4wmiozOxb98oqKoqV0kbDFMXsZE5U2cQ\nEdTVi+ba/fwWV9kPDIapTdjInKl3OI7DlSszEB+fwzpyhvkIG5kzDMPUYmydOcMwzCeEdeYMU4OW\nL7+D9u0PIjw8vaZDYeq4SnfmHMe14DjuIcdxbziO8+U4brkiAmOYT8HTp9EICUlDZGRmTYfC1HGV\nnjPnOK4pgKZE5MNxnA6AlwAmElHgR+XYnDnDfCQxMQfh4Rno169FsevBDg4IdXLC0M2boa6rW0PR\nMbVBtc2ZE1ECEfn89+ccAAEAjCtbL1N7CQQijB59DpMnX4RYXPkf0EKhGEuW3AaP5yJ1b+fOp7Cy\n8qx0G7WVkZGOVEcOAA/++APu+/cj+PbtGoiKqYsUupqF47jWAFwAfPZfx/7hPTYyryeSknJhbLwH\nysocUlN/g7a2WqXqCw5OhZnZISgrcygo+BPKykVjjPDwdLRpcwAAkJ+/7pM6uizs/n2E3b+PIevX\nQ01bu6bDYWpQta8z/2+KxQ7Azx935O/weDzJny0tLWFpaamo5plymjvXHunp+Th6dBz8/VMwdGhr\nuU90MTTUxrNn30FVVbnSHTkAtG/fCKdPT0KjRpqSjhwATE31sW3bMOjpaXxSHTkAtBk+HG2GD6/p\nMJga4OLiAhcXl3I/p5CROcdxKgBuAbhDRPtllGEj81pCIBChQYPtEAhEGD68De7dC8OZM5Mwe3bX\nmg6tThMIRHjwIAyff94aWlqqNR0OU09U9zrzEwD8ZXXkTO2ipqYMV9f5ePhwLoYPb4PWrRuyHZUK\nsH27K8aMOY+1a+/XdCjMJ0gRq1kGAngMwBcA/ff1BxHd/agcG5kz9dqdO8FYssQB27d/ga+++qym\nw2HqCZYCl2EYph5g2/mZWo2IcHHKFJz54gsI+fyaDqfOUcSSUKZ+YZ15HRIQkIzs7IKaDkMhRAIB\nQu7eRcSjR8hPL3kre35+IdzcolFXfqPz80sqc1t+QkIORo60wZ49bhVuJzQ0DY0a7cT48RcqXAdT\n/7DOvI54+DAcnTodxqRJF2s6FIVQUVfH9+7uWODmBt1mzUoss2yZAwYMOIF//639m4aiojLRvbsV\nevY8Wuqo2dMzDk5OoTh1yqfCbeXkCJCdXYD4+OwK18HUP6wzryOaNNFCw4YaaNtWv6ZDURgjCwsY\n9+4t876FhREMDDTRvn3tPvAZAPT1NWBhYYh+/VpASUn29ObYse1hYzMZ06Z1Qv/+1ggISJYq4+ER\nC3PzQ/jnnxcl1tG1a1OEhi6Hs/NchcXP1H3sBWg9c/y4F/75xwMnT05Et25NazocRoYxY87hzp0Q\nHDs2Ht9/36PYPSsrT/z4420MHNgSfL4Qu3aNwNChpjUUKVPT2ElDn5Dk5FyMGHEWPXo0Q0YGHz4+\nCXBzi2adeS129Oh4uLpGYvr0zlL3Fi7siZYt9WBvH4hjx7zg6BjKOnOmTGxkXsM8PeNgY/Mav/8+\nCEZGOhWqw8cnAd27W6FVKz28eLEQrq6RmDixA1RU2CxaVfPxSYC5eSNoaip+x2dmJh83bwZh0qQO\n0NGpfNoEpm5i68zriIkTbXHjxlts3/4F1q4dVOF6nj+PQbNmOmjVqqECo2NKc+GCL7755ipmzuwM\nW9tpVdrWjRtvYW7eCObmjau0Hab2YdMsdcS6dYPRsmUDzJ1bubwoJaVRZcrH0zMOkZEZmDq1k1zl\nW7bUQ4MG6jAzU/wL2sDAFDRtqoOGDTXg7ByOiRNtYW7eCIGByxTeFlM/sN/Da1ifPsY4dGgMmjWr\nPQcQ2NsHwtMzrsrqF4nEyMkRVFn9FTVmzDlMm3a51M+ekJADB4dgEBEGDmyJ+fO7ISdHoNC18O7u\nMejU6R+MHGkDAOjc2RBDhrTCjBnS8+tM9SvIzoaosLCmw5DCOnOmGG/veEyefBGjRtlUWRtTplxC\n48Y74eOTUK7n4uOzcePG23LvfoyIyICbW3SZ5RYu7IHx481gbi57pD1r1lWMHXsely69QUpKHvbv\nd8e+fc+Rnq64XayNGmmhcWMtyYjf0FAbjx7Nw6ZNQxXWBlMx6eHh2N2sGU4OHlzToUhh0yxMMe3b\nN8K4cWbo2FG+uVkezwWvXiXizJlJ0NVVl+sZgUAEkYggFIrlKs/PzISGnh7mzbsOJ6dQnD49CXPm\nyD8tZWl5CpGRmfD2/qHUFT5bt34hdS0iIgP6+hrQ09MAAIwc2RZJSbno2rUpmjTRhp3ddCgpcTAw\n0JQ7nrK0a2eApKTVCquPUSAigAgkEtV0JFLYC1CmUoyM/kZSUi7c3b9Hnz7ynRZYWChCZmYBGjfW\nKrOs94kTuLFgAYZt2wYvraE4edIH589PRadOTeSO8fvvb8DTMw737s1Gkybyn9rz5k0SunWzQqdO\nTfDq1Y9yP1dd/PySsGqVE37+uS9Gj25f0+F8MviZmVDV1ISyWvWsMGIvQJlqcfPm1wgLSy+1I582\n7RKePo2Gu/v3MDHRg6qqslwdOVCUw+Xd//78ez/8/HO/csd4/PiEcpXPyytE797HoKLCoUEDdRgb\nK+59hkAgwqBBJwAAT58WndZUUTduvIWjYygaNtRgnXk10tDTq+kQSsQ6c6ZC0tPzsWbNfYwZ077E\n3N3x8dk4e/Y15s7tirdvU5GUlIv09HyYmJTv/wi9fvwR5hMnyszfUhXy8wsRHp4OdXUVxMX9Ak1N\nVSQl5eLVqwQ8eBCOly/jcenSNOjrl39qJT+/EG/eJIPjAD5fKNWZv36diE6dmsi1R+Cnn/pAR0cN\nkyZ1KHccJQl78ACqWlpo2b+/QupjqhfrzJkKcXAIxrFjXvD0jAPHFaVknTy5o+T+li2PcfiwJ+Lj\ns/Ho0TwkJ+dWeI10dXbkQNELSFvbaTh37jXi4rLRtq0BvvnmCh48CIeBgSbS0vLx9m1qhZaD6ulp\nwNd3MQBIvWPYv/85VqxwxC+/9MPu3SPLrEtXVx3Ll/ctdwwlyYqJwdkRI6Csqoq1WVlQUZfv/QdT\ne7DOnKmQyZM7YsOGNPTta4yxY8+DCEhKWiWZk548uSMCA1Mwa1YXGBhoVugFIZ8vRP/+1lBTU8aT\nJ/MrNSVRXpcv++PSJX80baqDHTuGY9gwU8TGZmP79i+gpMRBLCb8/fcz/Pxz3zLjysjgo39/axgb\n6+L+/Tlo06bkZGktWjSAmppyhTd+kVgMTqliC9S0DQ3RcfJkqOvpsY68jmIvQJlKISL88osjRCLC\nvn2jJBkDe/c+Bi+veLx48T169mxeobrT0/PRosVeKCtziI//Fdra1bel/e3bFPy19SHOnPNHZ1M1\nvAr5o9j9Tp3+QUBACq5dm1nmNEdUVCbatTuAhg01kJi4ChxX5ruscotwcYHNyJHo+/PPGLFzp8Lr\nZ2oOewHKVKlkf3883rIF/VaswN69o6Tum5o2RFhYOho21KhwG/r6mvDzWwwlJU4hHblYTKWmp/2Q\nuXljzBquDZuzheBHx0jd//33Qbh/PxyWlq3LrMvERA9+fkugra1aJR05AOSlpkIkECA7NrZK6mfq\nACKq9BcAawCJAF6XUoaY+sPpt9+IB5D9d9+VWTYyMoN27HCllJRchceRHR9PDsuXU6ynZ6nlNmxw\nJhWVTXT7dpDcdYuEQnLe8TcF3XWsbJhERJQVF0dnv/ySnu/fr5D6Ppby9i0V8vnleuaNnR3ZjBpF\nKUHyf1+Y6vVf31lmP6yoHaAnAZT9xoapN/qvXIkh69djyJ9/lll2wwYXrF37oFJHpcnifeIEXhw4\nANctW0otl5CQA6FQjJSUPLnrVlJWhuWaX9F+5JcAivKljBt3Hg4OwTKfsZs5E8d69wY/M1PqXqy7\nO0KdnOBz8qTcMeza9RSHD3vIVbaRmVm557t9TpxAyN27CLlzR67yqUFBENfCDTMMFDMyL/rhgVZg\nI/NPnr19AK1Zc4/y8wsl1x49iqAJEy7Q69cJctXx6lUCff/9dQoKSimzbFZsLN1eupRi3N1LLScQ\nCCkwMFnm/SVLblGDBtvJyytOZpnt210J4NHUqRdllvmrUSOayPWgIQP+pbCwtGL3xCIR+Zw5Q0lv\n3pQa6zuRkRkE8AjgUU5OgVzPlFdqSAi5HzxIgry8Msu+PHaMeAA5LF9eJbEwJYOcI3M2Z/4JCgxM\nga2tH5Yv76vQbegA8MsvTggLS8fgwSYYO9YMADBkSCsMGdJK7jr273+OEyd8oKengb///rLUsrrN\nm2PMoUNl1qmqqgxz88b49VdHuLpGwczMAN980wVjxhRttgkLy0BWVgGSknJl1rF0aW+oqyuX+sLz\ne3d32M9wgNuzBDx+HAlT0w9WrnAc/JS6QknUCPLsXzUx0cPOncOhra1WZS9/Ddq2RZ9ly5Ds7w89\nExOo6cjOqa/ZqBE4JSVoGxpWSSxM5ShsNQvHca0A3CSiLjLu04YNGyR/t7S0hKWlpULaZspn2rRL\nuHIlAFu3DsMffyg2YdD164F49iwaGzcOhYZGxcYKb9+mwMrqJVas6FfuTUZladv2AMLC0gEA/fu3\nwLNnCwAUHZIcHZ2Jjh3lTxMgS0hIGh4/jsTs2V2KLVu8ezcEo0efg76+Bvz8lqB589qRKTPUyQk2\nI0ei7ciR+Pbu3VLLioVCKKmwMWBVcnFxgYuLi+TvGzdurN7DKeTpzBXVFlM5jx5F4OhRL2zePFTm\nmufqEhGRgTVr7uO777ph5Mh2cj0TGZkBVVXlCnWGwcGp8PdPxsuX8Rg9uh36929Z7jrkEefpCQ19\nfRi0bSu5lpSUCwuLI0hKysWmTZZYv/7zKmm7vBJ8fHDK0hIWs2Zh7D//lFlekJuL3KQk6JvW3FF2\nBdnZsJ0wAfpt22LC8eM1Fkd1kHdpoiLnzFsD8C3lflVNKTE15PnzaNLT206//eYks0xiYg69eBEj\n8/7u3c8I4NHYseek7t27F0r9+h2nhw/DJNeSknJIU3MLGRj8RQUFwsp9gBK4uITTpEm29OZNktS9\n1NQ8mjLlIh08+H5+vm/fY2RispfS0t7POacEBRGP4+ivRo2k6rCze0PKyhupX7/jCo+9srLi4kgk\nLPt7enrYMOJxHEU+eVINUZUsOSCAeADtaNiwxmKoLqjO1Swcx50H8AyAGcdxURzHzVdEvUz14vOF\n8PCQb51yRgYfU6deQmZmAYKC0mSWGzXKBn36HIera2SJ92fP7oKlS3tj584RUvdu3QrC8+cxuHUr\nSHJNU1MVLVvqoW1b/Uqdcfr8eQxWrnTErl1Pix1GceyYF+ztA3Hxop/UM+7uMbh6NQAHD74AUDQQ\niozMRGJiDnJz3x9WoG1oiOY9e6LNF9IpdTt3NoSqqjJ8fRPRrt0B5ObKf0hHVmwsDrZvD/u5c6Xu\niYVCpIeFSf4u5PPx8tgxeB45gr8MDOBz+nSpdQfa22NP8+a4vWRJmXHoNG0KNW3tGk041bhDB8x1\ndsb8J09qLIZaR54eXxFfYCPzcvv99/s0ceIFys6umpUMH/v+++sE8OjIEY8yy0ZEpJOKyiZq2HAH\nZWfLXtv84483ydR0n9TKjsdbt5Lzhg304483ZbaZnp5PJ054UUZGvtyfISengMaPP0+//iq9NtzK\nypNmz75KmZl8GjHijGSlSMeOhyRlwsLSaNu2x5SaKr26QyQS0+7dz+jQIXcSi8VEVPSbQmRkRomx\nhIam0dat0nUlJmaTgcEO0tbeWq6197GensQD6FCHDlL3bixcSDyAXp8/T0REz/fvJx5AB83NiQeQ\n0+rVpdYddPs2bVRSojs//yxXLGKRSO64mcoBW81SPoWFIgwceAKFhWK4uS2o8Ms7RbKyeom0tHwE\nBaWiR4+qTzbVoUNj6OmplzqPHhmZgWfPojFjRmf4+i6Gjo4adHRkr20+cmSc1DV+ZiYerlsHADD9\n0x6amiolvuhs2FAD8+d3L9dnCA1Nx82bQdDTU5daCbN9+xNERGRg1iwLrF8/BCYmekhMzMXYse/T\nx5qa6uP334teCufkCKCj834ViZISB2trb/j7J0NXVx1z5nQtNT/6unUPYWvrh8JCETZssJRcNzTU\nQUDAMgiFYkRHZ+HKlQAsWNAdysol/6ZBRHi0aRN0mjbFD97e0GkqfcCGVuPG4JSVodGwKK9LmxEj\n0HroUPRYuBB6JiYw7tOn1O9b+zFjsDYzs9TVLB+qaA4YpgrJ0+Mr4gu1fGSek1NAurrbSEtrK6Wn\nyz8SrEovXsTQlSv+VVK3tbUXNW36N925E1yu54YOPUUAj86c8alU+94nT5KnlVWl6pDF3j6gxHn6\nx48jaP/+5yQUlj2qXL/+IXEcjy5d8it2XVl5IwE8unjRT8aT7zk7h9PUqRcpIED2+vbPPjtMAI/s\n7QNklnk3P7xRSanUEbFQIJB5L/TePXqwbp1c68mJiAr5fLKdNImuf/+9XOWZqgM2Mi8fbW01+Pou\nhlhMlconokjNm+uic+eqWdP78mUcEhJy8Pp1IkaNkm8VCQDMmNEZBQUiuVeBvH6diL17n+O33wYU\nW/bXbd688oYst4kTS14HPnhwKwwe/H69O58vRGYmH0ZGRaNRkUiM589j0Lu3MfLyCkEE5OcLi9Vh\nazsNUVGZmD69k+QaEcHRMRQ9ejSDoeH7kbqlZesyc7csWdILjo6hGDBA9vezkbk5vtixAzpGRqWO\niJVVVWXec1y5Ekl+fjDq2hWdp08vNSYAyImPR+D161BWU8O4I0fYcsS6QJ4eXxFfqOUj89rG2zue\nVFQ20YAB1lVSf16egB4+DJNrlFoZixbdIIBHK1bcqdJ2KsLS8hSpqGyS7PrcsuURATz69VdHEonE\nFB2dKfWMUCii776zp59+cpBcs7F5RQCPRo2ykbvtwBs3KPju3cp/iA/42trSFg2NEn/jCbx+nW4v\nW0b8rCy56wtxcqLo588VGWKpUoODKT08vNraqytQzblZGAXT1FSBhoYK9PWr5rcETU1VDB1qKnOe\nVlGWLeuDVq2K5qYVQSQSY+fOp7h9O6jswv/Jjo9HarB0PhVtbVWoqSlDXb1o1NmpUxPo62vgs88M\noaTEoUWLBlLPJCbm4sQJHxw+7CFZiWJhYQgDA03k5xdKlZcVj+3EiTg3ZgzsZs5EsIOD3J9FlsOH\nPTB55VvE8hsgPTxc6r75hAkYc/Ag1HXlX5vfdsQItOirmMMvypKfloYjXbrgSJcuKMyTP38O8x7r\nzGspc/PGSE39DTdvfl3ToZSpsFA68VJmJh9RUZlQU1NGZGQmLl16g+zsArnrfPgwHG3bHoCNzWsA\nQG6uAESEJ0+isGbNfcyff12ueogIVt264XCnTsiIiCh27+bNr5Ga+hvi4rJx/34YJk/uiLS0NbAQ\nvMA/HTsi1kM6wVXz5rqwt58JB4dZki32TZvqIi0tH48fRyIgIBk9ex7FqlVOUs/affUVrPv3h4qm\nJrp/9x1MBg3Cm0uX8HzvXrm/L7LcuROCkHiC3vQVGFZG0jEACAtLx6FDL5CXJ98PoKqmqqWFJp06\nwcjCotoOSq5v2ERYLaamVn0n61TU6tVO2LfPHXZ203HzZhAGDmyJ+fO7o18/a4SGpknyeOflFSI+\nPkfqqDRZPDxiERaWDlfXSGhoKOPrr69i7tyuWLq0NxYu7IG+fWUfIP0hjuNgaGGBzKgoqH+0Lprj\nOOTmCjBypA0AICVlNfT1NRH+4AFSAgMR5+EB4969per8eE7e0FAbtrZToayshOTkPHh5xSM/v7DY\nahoiQqijIwqyspCfmooJx4+jIDsb7vv3w2z8eLk+S2n2bu4F7VuboHY5CIXWy8ocga9YcRc3bwZB\nKBRjxYryH5KtaCoaGljk6VniPRKL4X3iBJr17Ilm3cu3uulTwjpzplJSU/MhFIrx4kUsrK298eBB\nOObP7w5jY12kpORBW1sVHh4LkZqaDzOzRnLVaWPzGjduBOHIkbFQVuYwfbodACA6OhM9ex5Fz57N\ncfRoyR1gcHAqJk26iBkzOkmWA865f19mW/r6mpgzpys4ruh8Th+fBNxVGY9vDk9Az0VfScpFRmYg\nNTVf5hLRmTPfH2p999YMmHUqvnxQJBCg07Rp0G/TRrLFX11XV64UwvJo26Ut5q8cDRXNKXJNpcyf\n3w0CgQijR8v38jszOhr+dnbo/t131b5ZKPjOHdxcuBCNO3TA0oCAam27TpFnYl0RX2AvQOslgUBI\nYWFpJBSK6O+/nxbbev9uY015jRt3ngAeHTrkTnZ2b4jjil6gvn2bQk2b/k3ffnuFYmNLfpFnbx9A\nAI8GDTpRobYXLLgueQn6oWbN/iaO45W6zJCIKOLRI9qkokI3Fi4sdj3EyYl4AO1q2pQyoqIqFNuH\ncpOTyef0abmXGlbWlVmziAeQM49XLe19KDc5mS5MnEjPDxyo9rZrA8j5ApR15kyl7NjhSr17H5Xa\n4VkeUU+fkse//5LXiRPEz8ykqKgMOnXKW5J75eNc3tOnXyKAR6dPl7zW3ckphOLi5F+18aHg4FRa\nvdqJoqKK7+qcMeMydelypMSdodnx8RTrUbSDNfjuXeJxHNl99ZXkvv/Vq7TT0JBsRo+mTWpqtF1P\njwrzK7eXwX7ePOIB9Gjz5krVI68QJyc6N2YMJbx+XS3tMe+xzpypEqtXO9GSJbckSxr79TtOAI+u\nXi2+uSk+PpuWL3cgD49YIiI6etST/v77aYl17mzcmHhA0bbz334jIqLCQtlLJtu1208Aj+bMuaaI\nj1Rphy0siAdQ1LNnRESUGRNTbAOP84YNxAPo7ooVdNjCgo736ydXQqvS+F+9Skd796ZoN7dK1aNI\nmza50NChpygxMaemQ6lXWGfOKFxOTgEpKRXtgIyJKVqDHRGRTvb2AVJTKtu2PSaARxMnXqC8PIEk\nD0p4eDqtX/+Q9PS205MnkURE9GjzZjo5ZAgd7dWLwl1caO7ca6SpuYVeviz51B8Xl3CaPv0ShYSk\nKuRz/f33U7KyKv0M0dLcWLiQDrRrV2z6JPr5c0lWQVFhIYU7O5OwQHE5drJiY+mfTp2kpnNqkpnZ\nQQJ4dP9+aE2HUq+wzpwpRiwW0/btrmRj86pcz6Wm5tGFC76SY+Bu3nwrtcW9JLGxWbR06W1ydy/a\nVn/gwHPasuUREb2fJrlwwZeIikbxHxo92oY4jkfOzuGSa5cvv6GVK+9SXp7sLesfevMmiQSCske/\nERHpkh80ublFdVtZedK4ceel4pJXQXY2bVJVpY3KypSTJJ1KtzRCoYh8rt6W2jwTev9+0TTHq6J/\nvxh3d0kirdKEu7hQ3MuX5Yqhovz9k8jOTr4j8Rj5sc68Hjp79hWNGHGmQiPS168TCOCRisommjbt\nEg0caF1qNkY+v5CcncNpzpxrBPBo69bHlQm9mJSUXDp3rmju9V0+8927n0nu5+UJKDw8vdgzpqb7\nCODR7dty5CV3AAAgAElEQVRlnyJvbe1FAI8WL74lVzxLltwiI6NddPjwCyIiatfuAAE8Wr68YrtW\nxSIRXZ45k2wnTSKhQEDp4eFkO2kSBVyTPS0UG5tFGYkp9EurwcQDaGer9sXuX/32W+IB9PB//5Nc\ni37+nDJjZOeKTw8PJx7H0RZNTRIVvj+T1fvUKTrevz8l+pX9Q5mpefJ25mzTUB1y+vQr3LsXhocP\npXf4laVzZ0P8+edgHDw4Go6OIXBziyn1pPq1a+9j6NDTEAiE6N69KYYObV3xwD/y44+3MGvWVfzy\ni6MkJ7mq6vv/FDU1VdG6dVH2PyKCSCTGoUNjsH79EAwf3qbM+o2MtKGqqgRjY12EhaWXmDOciODr\nmwixmGBhYYTExFw4OYWCzxciMTEHAPD4cQQ6dz6MV68SyvX5OCUlTLO1xcxr16Csqoqg27cRaG8P\nzyNHAAC5SUkQFb7frPPENQItW+7FiP4H0CDSFSIowaDvkGJ1KqmpQb1BA7QfM0ZyrUXfvmhgLHu9\nvbaREdqNHIlO06YVy60SYGeHGDc3RLm6lutzlVdWTAwuz5iBgGvXqrQdpghbZ16HHDkyFs7O4Zg7\nt1u5n1VS4rB58zAARQmgsrMLJB1mSbp3bwYjI218840Fxo83l1kuMjID48ZdwNix7bFjx/BSY0j2\n90dOYiKEwqLjAx0dQ/HmzRI0baqDoKBUCIViqKgoITIyA3/+6YzJ3fJw+F9v+GU3hZ/fEsnhy6W2\nkZyL4cPboKDgT1y7Foi2bQ/g889bwcVlXrFy27c/wbp1D7Fhw+cYObItlJQ4pKbmAyja5s9xHAwM\ntPDwYThevoxH167SaWfl1X3+fIgKCtBu9GhEurritKUlOkyahBlXroCfmYmrUydCiSYgLywQAKBv\n0hILLhY/Ci09NBQFWVkoyMyUu11VTU3MunNH6vrYf/9FZ2dndJ45U+66MqOjcX3ePHSYMgV9li6V\n65mgW7fgf/ky8lJS0HHyZLnbYiqGdeZ1SLt2BmjXzqDS9XTo0LjMMnPmdMWcOV3LLBcUlAo/vySo\nqSljx47Sy576/HPkpaRgj7s3LCwMJRtWZs++BoFAhHbtDPDVV5/h4sU3uGbzDG1tdmMgOHhq/gmB\nQDplwDtiMWHXrqdQVlbCH388wMCBJti0yRLTpl2CsjIHY+PiOVZEIjGMjXWhoqIEPl8IXV11aGio\nwNBQGxoaKggJWY631+3hsW8vftj7G6bNe//D09bWD+vXO+Pw4TEYMaIt5KGqpYX+v/wCAIh2cwOn\nrAwVjaKcOyKBAI1zQ7Be+xDEOZn4bNYsjNqzR6qOmdeuIT0sDM179pSrzdI0MDZGl2+/LdczMW5u\nCH/4EILcXLk7865z5qAgOxvtRo2qSJhMObHOnKmUESPa4sGDOXLt7uwwZQpSAwPRzNwUm/oUdZBE\nBE1NFQgEIixYcB19+hhj4cIeyM3OR0uPIOjq6yJw7yo0bSp7V+PLl3FYu/YBNDVVoKKiBE1NFVy4\n4AdlZSXMnNkZAwe2xLp1D7B58zBMmXIRDx+Gw8NjIQYPNsFffz3FwIEtkZ6+BmpqykhMzMGRI55o\n42OHBI/n6DzlLZSU3p+l+/hxJEJC0uD2LFLuzvxDLfv3x5q0NKhqaQEAtJs0gbKGBoSZmfjRxxtG\nXbqA46TP7tXU14emAjryiuo4dSomnTlTrsRbqlpaGLh6dRVGxRQjz8S6Ir7AXoDWaWFhaTR/vj35\n+iYovO6kpByysDhMamqby9xhWZLCQhGtXu1ER496Um6ugKKjMyUrVNLT80hFZRMBPPL1TaQBA6xJ\nWXkjvXwZR2vX3iMTk73k55coqWvt2nsE8Gj+jJPkdeKE1A7LYFd3Wt2wLW1UVlZYCttz48bRVi0t\nsh4woNLrz5n6B9W5mgXAKACBAIIArJFRpuo/NVNlOnX6R+q8TEXKyxNQYmIOZWcX0JEjHhXewfnO\nP/+8oGPHipbkWVt70aRJF8jPL5FOnvQmf3/ZywXfvEmihb2X0DaDJhR0+7bUfbe9eyUbnN6dt/mx\n8qYxyE1Ops3q6rRVW5tuL11Kr2zkz4teUtvRbm5UyJd9Lus73t7xlJlZdrnq8vrcOYpxd6/pMGod\neTvzSq9m4ThOCcAhACMBdAbwNcdxJR/1wtQ5mZl8FBaKMGmSOXR11bBwYY9K1ykWCqWuaWqqwtBQ\nG3v2uGHx4ttYs0Z2cix5aGqqYOPGR7h3LxQpKXmwt3+LKVMuYv786zh69KXM5zp1aoJJfZUgSEtG\n4uvXUvd7/vADun33HSafPQuLr6XTE58eNgx/GxoiOz5e6p77wYO4Nns2CrKzi13XatwYP756hTGH\nD8Pjn39wf82aCnzi/9o4cADW/fvDceXKUss5OASje3crzJhxucJtKVK0mxuuzpoF24kTazqUOksR\nSxP7AAgmokgiKgRgC4D9i9QD/v7JaNp0N7780gaPH0chJ0eA3r3lSz0ri5e1NTarqcHz339LvD9u\nnBksLVtj5szOctVXWCjC+fO+iIsr3kG6u8ciJiYLXl7xGDfODF9+2RazZlmga1cjfPll6XPdiT4+\nUNfTg0UJLwmDHRzgc+IE3D54SZnw6hX+7doVnlZWiPfyQl5KCl6fPSv17LOdO/HaxgZxJaR6bWxu\nji7ffAPLTZsw7qPvjc2oUTjQti3yUlJkxsznC5GQkINGZmbQ0NdHk86lf/+MjXXRpIkWOnduUmq5\n6mL42WfoMHkyei1ZUtOh1FmKeAFqDCD6g7/HoKiDZ+qRdu0M4O+fjMaNtSpVT15yMkCE3OTkEu/3\n6NEMzs5z5a7v2DEvLF3qgPHjzXDjxtcIDEzBy5dx2LVrBKZM6Yhhw0yhoqIER8eijvl//7Mss86c\nxEQI8/MhKpA+TMNk0CCYjR+P9mPHSq7Fursj8fVrhNy5g54//IBnu3ahcceOUs9OvXABib6+aP35\n5yW2q6Sigs/Xr5e6nuTnh9zERPAzM6HVuOSVSCNGnIW7eww8PBZiTVqa5LpYJELiq1do2q1bsfND\nu3ZtiqSk2vNyUl1XFzOvXq3pMOq0al3NwuPxJH+2tLSEpaVldTbPlFOnTk0QH/8rtLVVoapavoMy\niAje3gno0sVIsjEIAAatXYuOU6bAoH3Za8blMWRIK/Tta4wpU4o6z4kjrBAUI8TRbV2w8Pey1zaL\nhUKE3rsHk0GDJHnAR+7ZA387O6g3kD42TsfICMO2bi3WqXZfsABaTZrAZOBAaBsaYtiWLSUermwy\naBBMBg0q92dc5OmJgqwsSR70kuir5WMQnkKcNQXA+zXxzv/7H55s24ZhW7di8B9/lLttpvq5uLjA\nxcWl/A/KM7Fe2heAfgDufvD3tSjhJSjYC9BPyp49Rdv0f/ml+IqPGzcCqU+fY/T0aeVzepdkfo8l\n1AazaemEnTKzNH7o6a5dxAPo+oIFkmunhg4lHkDuh6Rf9ib6+hKP42hf69aUHhFBOxs3JtvJkysd\nd9TTp5XKc+7w00/EA+j2smXFrr88fpy2aGqS74ULlQ2RqSGQ8wWoIkbmHgDacRzXCkA8gK8A1P6D\nK5kq1bp1Q2hqqqBt2+KbnK5ff4sXL2Lh5BSKAQNaVrodsZjA47nAyEgbS5f2wRHn7QhyeYouE18A\nN+5h2rROaNVK9k7X5r16oaGpabER87AtWxBw9SosvvlGqrxWkyYwaNsWTbt1Az8jA/lpaUgPC6vU\nZ4h98QInBg5EIzMzLHv7tkJ1dJk9GxkREVKbgXosWIAeCxZUKj6mjpCnxy/rC0VLE98CCAawVkaZ\nqv8RxtR6KSm5dOqUt9SBEyUpKBDSwYPu5OtbtA48Pj6bevU6Wmy0HxCQTACPOI4nybFORHTwoDtt\n21Z6crDQ+/cpNSSEMiIj6clff1FemvQBG4V8vszMh8HBqdTU6C+aP+dyseueVlbkefRoic/Ee3vT\ntTlziiW5yoqLo8MWFnTzhx9kxvrKxoae7NwptezR8+hRcli+vFzpdfNSU+nRli2UGhIi9zOVlR4e\nXuJST6ZsYFkTPw3x8dl0+fKbUg9zqK3y8wspIUE6zaxYLKanT6Po+PGXBPBowABrIiJ68iSSAB6Z\nmR0sVn7Xrqd05kzJpw69w8/MJJeNGynep6hc1NOnxAPoQLt2dHX27KKMhOvXSz13auhQ2qSqSlFu\nbmQ9cCBdnDpVcs/VtSievn2PERFRekQEuWzaJFmH/i6f+YduLFpEPIC2amvT8f796eH69bSvTRva\n07JliT9M3n0/NiorEw+QynS4Q1+/6GCMp2VPKb3jsnEj8QC6PHOm3M9U1pGuXYkHUIijY9mFmWLk\n7czZdv5aIiurACKRGPr6muV6btGim7h5MwirVw9Ar17NMWOGfEv6ZMnNFUBbW61SdXyMiHDihDc6\ndWqC/v3fT62MGXMOrq5RcHf/HkqvHRERnIAV5zSgr68BH59EzJjRCRMmmGPatI4QCsV48iQKf/89\nQvKyc9kyB7x5kwx7+5nQ09Mose2kpFzk5AiQfNsGLhs2IPrpU3zr6Aj9tm3RvHdvtBw4EB0mTgQ/\nIwOdZ8yQel5FQwNKKioQZGcjxs0NqtraEItEUFJWRrP8t1imfgIjh30PAHiwdi38bG2ha2yM7NhY\n+NnaQllVFaLCQpgMHAgAGLRmDZSUleFlbY30sDA83rwZ4DioqKujMC8Pmvr6UjFwHIfRBw8iOzYW\nTT5aJTPp9GmkBAaiRb9+cv97dJ45E4mvX6PXjz/K/UxlmU+YAGVVVTTp1Kna2vzkyNPjK+ILbGQu\nE59fSM2a/U36+jsoPb18Z0NaWXlSjx5WpKm5hQAevX5d8e32VlaeBPBo3z7FHkXm7BxOAI9attwj\nuVaYn0+r9VrTQq4lvXoZRTyOo8UwJCWljdSixR5q1Ogv2rPnfY5zB4cgAnhkarpPcs3QcBcBPPLx\niSdf30T66is7evEihm7efEvR0ZkUeOMGNVb/jToqz6JXzi/pyqxZFHyn5BzlZ0aMoF1GRpQVG1vs\nulgkIkFuLhERxXp4UHJAgOSe98mTxAPIdvJk2t28OdlOnkznxoyhYAcHcvz1V0r09aXNamolHlKR\nERVFucnJ9GjzZvI4ckSul58ioZDifXxILBZTTmJimeWZ+gFsZF53KClx0NIqSruqrCydZKk0ixb1\nxKJFPfHHHw8QFpZeqayKeXlFObZzcwvLKFk+PXo0w4wZndGv3/sNR4V5edATJEJfRwUdOzSC4PBh\n5CQmYtn0H9GkiTaaNNEuVoeFhRF69myGKZYNkJOYCB0jIzg7z0VcXDa6dm2KVaucYGvrh9TUPNy7\nF4beXQ0wT2iFwQVK6ApfRNg1xhQbG5kxZkZGIj8tDYKcnGLXOSUlSVKs5r16FbvXbd48NOvRAxlR\nUQi8dg389HTMdXYGALQbPRokFqPD5MkQFRRAo2Hxl7B6LYt+Qxny558yYxILhcXykD/44w8827kT\nbb78EmFOThh39Ch6Llwo83nm08IVdfzV0BDHUXW1VRcJhWKIxQQ1tfKt564ssbjo8Id368ijojJh\nYqJXLW2nhYZCSVkZDVu3lllm3jx7xMfnYMiQVvj7zytYzh1Cg6aG+DUurli5uLhsWFl54osvTLFy\npRNapz5Gl8hzMBk8GFnR0fhyzx50nDwZWTExuL1kCTpNn46us2dLns9PT0dBZmapsZQm7uVL6Ldp\nU+I0SWkKsrNx6vPPoW1oiG/v3pVcz4yOxr9du8KoSxfM+2/Nscfhw3BcuRIdJk/Gm4sXMcHaGt2/\n+65C8TJ1B8dxIKIyR3msM//E9ep1FBERGXjzZgmMjHSqtW2BQITNmx+hd29jTJggfQCGWExo0GA7\n8vIK4eIyD1b7nNHl+QY069gOcx48KLXuN5cuwd/ODqP27YNu8+aS66/OnoX9nDloOWAAvnv6VGGf\nRVhQAIelS9HIzAwDf/tN7ucyo6NxoE0bqOnoYHVKCpSUlZEcEIDUt29hN3MmDNq3xxI/P6nn8tPS\noGlQ+dz2TO0nb2fO5szrmHfL9SozN/6hdu0OkKbmFoqKylBIfR/KzRVQZKTseu/cCSaAR4aGO2WW\n8faOJxeXcKnrYlHFVu8ICwrI/dChUs+/TA0OphODB9PLY8fkrtfrxAniAbRFS4tSU/Po+PGXlJEh\n3/uPRD8/SgsLI6KipZBbtbVpk4oKxXl5ET+rePbIqGfPyHXHDrmyIjL1A9gZoPXTpUtv8NNPd7Bo\n0S2F1OfltQgRESvQsqXip1bGj78AU9P9cHOLLvE+/febmlBIWLnyLry9pTMNduvWFJ9/3rrYtZyE\nBOxu1gwnBw8ud0zKamros3QpDEtJRBXp6oooV1f42drKXW+chwcAoP3o0diwwRnff38TW7cWnbF5\n4/vvsVlVFZ5WViU+a9i5M/RNTSXxmQ4bhpYDBqBxhw6SFAPv3F68GA/WrkWgvb3csTGfBvYCtJb5\n7bd7CAxMwfnzU6GjI71EcNgwU0ycaC5ZnldZurrq0NVVV0hdALBz51PY2vrh998HoUkTLWhqqpT4\nOQBg4EATfPONBbKy+Ni3zx2RkZm4cGEqDh16gaFDTdGjRzNJ2YiIDOjpqUNfXxOF+fngZ2aWmGa2\nJPzMTKjp6EBJWb73EV1nz4aSigpaleOHxZD169HQ1BTd5s2DvncmHl5xQ/OEx+Bn9kHUkycQC4V4\ndfo0mvXoAW9ra/T9+Wck+fnBbOxYyQtWoOhX6pYDBkBdTw+qmpogomInDw1auxYhd+6g7YgRcsem\nKC/++Qdex45hio0NDD/7rNrbZ8ogz/BdEV9g0yxyadhwh2S5nTw8PGJpyJCTdOvWWyIqWuZ48qR3\nlUybyOPdqT7vliHKs5kpJiaTFi++Re7uMXT27CsCeNSnz/spjsDAZFJR2VRss1BmTAzlZ5T9GWM9\nPGizmhpdmDixAp+mYjJjYoo2DnEcbVJTI+tBg+jmDz9Qkr8/XZo+nXgAHevbl3gA3f/jj2LPpoWF\nSTYdXZ09m7bp6FDCq1fVFntpzo0ZQzyAvE6cqPa2fS9ckHkYSH0HtjSxbnJy+haxsdlynwbv4BCM\nx48jYWraEGPHmuHYMS/89NMdTJhgjuvXv6riaKWdOzcZ1tbeWLGiaBPLhxkTZTE2boDDh4tSyrZq\npYeZMztj4sT3L0R1ddVhYKCJVq3eTwU1MJYvr7rkP3aR7AOhFa2BsTEmnTmDsAcP8ObCBaioq+NI\n9BD8b9ZjXDr0Gxq0aAGjLl0gyMlB648yh+qbmmL4X39BXU8PAXZ2EOTmIv+DlLblYWXliRMnfHD6\n9CS5DvEuSUFWliR75ARra0S7ucF8woQK1VVReampuPLfQSBthg+HdpPakYO91pGnx1fEF9jIvEpk\nZxeQlZUnxccXbYt/8yaJBg8+QTY28o/m8vIEdOaMDyUl5VRVmNUqLSyMTg4ZQi8OHyaiolwkosJC\nIiJ6vHUrnRg8WGpzkKIlvH5dtG1fR4cE/ALS199BKiqbSn0h/LGCnBxKDQ4u8V5SUk6xI9/SQkMl\n54cmvHpFCa9f08iRZwngkbW1l9xtigoLJXleHm3eTDyAfM6ckfv5qiAWi+ne2rXk9Ntv5T6Srz4A\ny83CyGvbtscE8Gj27KvV2q5QKKLRo22oW7d/KTExh+Ljs2nixAt04oTszicnMVGyI7MkGRn55Gp1\nmngAnRg8uNi9eG9v2tuqFfGAEpM+OSxfTrubN6ckf/+Kf6j/CHJzyXbyZHJavZqIiMLD0yUJw4iI\n/C5epNfnztG1OXPo+f795ao7Pj6bdHS2kYnJXhKLxZKVNHd+/pny09Nps7o6bVZXp+A30XTG2p1i\nvErPW/OOUCCg/W3b0u7mzYmflUXXf11HrTGXplvKXm3EVD15O3M2zfIJEghEiIjIgJlZIwDA8OFt\nYG//ttjURnXw8orHnTshAAA3t2jw+UJcv/4WsbHZmD+/u1T5lMBA/Nu1KwwtLLCohKPX+PkCmJsd\ngKCQcMPqNDoPf5/WlohwdsQI5KWk4Ivt29HAxARWPXqg+4IF6LN0KYCi4+Ky4+KQFRMjlQPlY257\n9+L+mjWYdvEiOk6WPgRDVUur2Mk5rVsX7QAVFhQg8Pp1XJk5E+A4gAhBt2+j7/LlcnzH/qtbVQla\nWqrQ01MHx3HQaNgQnJISNA0MoKajg1ZDhgAA3uxaj/AzZ3BcLMa0S5fQefr00ismQmFeHoR8Pkgk\nQqvZyxCx2wq5b+QOjalBrDOvB7KyCjB27Hm0bauPU6cmlVl+0aKbOH36FWxtp2LmzM/Qu7cx3N2/\nV1g8CxZcR2JiLuzsZkBDQ/Z/Yj16NMOiRT0gFIoxbpwZxGLCP/+MwaBBJiWWV1ZXh7K6utRyvXce\n/r4WhUlCqOg1QZevZqJBA3UIcnNxtGdPFOblIS8lBdqGhui1eDECrl5Fgrc3gm/dknTmM+3tkREe\njmY9Sj+0Oi00FCmBgRAXFiI7NlbO7wqQERmJFwcPwm33bjTv0wcgQpyHB4z79i1WLvbFCzzZvh1D\n1q8vMZZGjbQQE7MSyspF7yM6Tp6Mdfn5UFYrWjU028kJAHBp6lSQWAwVDY1iG6dkUVZTw7LAQIhF\nImg0bAiT9HCs6eGH3lOqf+UMU36sM68H4uKy8fRpFPz9Sz5X82PNm+tCVVWp0ud5lkQkEuPCBT/w\n+ULEx2fD1FT29nZlZSVYWY3H77/fR8+eR3H37iwsWdJbZnl9U1OsTk4u8Ug2AGjQ1BDLuD8xxcoW\nDRoULbcU8vnIiokBp6QEFW1t5CYlwf/yZXSbPx8aenpo0b+/5HlNff0yt+Onh4fjn44dodW4MRb7\n+sq9RM/94EHcXb4cHadNQ4MWLTD4999h3LcvPA4fRvf584uV9Tl1CoH29tBr3VrmD5aPj/F715F/\naOqFC8hNSkKDFi3kihFAsaPyIh89gqaXHXIapADrFsldR0ZEBABUODUCU0HyzMUo4gtszrxKubpG\nUkBAstzlRaKqe5Hk6RlL9++Hyl1eTW0zATxavdqpXO1kRkfTGzu7YrtB373o/FBGVBRlxcW9f6F3\n+rRUmdyUFLL76ivyPnmy1DZzU1Jon6kpnf7iC7njFBYU0PMDB4gH0LPdu8ssnxUbS4+3bqXsePmW\np1YVoUBAL48do5SgILmfyc/IoK1aWrRNR0dq9ypTMWAvQJna5upVf/L2lu6g/ph/mAY1mkehnr7l\nqu/dWZ2vbGyk7hXm55Pzhg0Uev9+sev8zMwS6/K/coV4AB3+7LNyxfDOKxsb2qqtXeIPitNffEHr\nVTSpt8ka+maSdYXqrysK+Xyy6tGDjvbqVa7TjxjZ5O3M2XZ+Rm4BAcnw9U2s0LPnz/tiypRLGDz4\nJAIDU4rd65JxD8NTTyHuXvlSFHSYNAnNe/eGcW/pqZlgBwc82rgRd5YtK3b9w2mED5mNH4+Re/di\ngrV1uWJ4Jy0kBIW5uUgLCZG6p6ymhmzowCNKE3Y3IiAWKybhXNKbN7j67beI9/Yu13NioRBPd+5E\nWBnJyipCRV0di16+xEIPjxKnfpiqw7ImMnLJyytEkya7IBSKERv7i8z59vDwdMyffx0TJpghODgN\n48ebY+DAljA03AWBQAwAmDXLAjY2UyTPpAYFwd/ODr2XLJHK+11RBdnZuL9mDVoOHIjwhw/RtFs3\n9P3pJwBAeno+IiIy0L17s1LrEPL5SPDxkesUHxKLEe/lhabdu0ulDRCLREgJCsLBH3ei2+TRGPWN\nJYIdHNB5xoxiW/nL6+6KFXDfvx89Fi7E+KNH5X4uxNER50aNgl6rVljx3/w2U3tVS9ZEANMA+AEQ\nAehRRtmq/V2EqVJCoYi+/PIstWt3gDp0OEQeHiVvurGxKdqOb25+kAAe9e59lAQCIX3xxWnq3/84\njRlzjqysPKst7ojHj4kH0C5DQ8m1wYNPEMCje/dKn9e/vmAB8QDyOHKk0nGcHz+eeBxH4c7OdPXb\nb4kH0OOtW4uVEYtElBEZKbOOQj6fwp2dSVRYSJFPntBmNTWy6tGD0sPDyxVLQXY23Vi4sMz3A0zt\ngGqaZvEFMBnAo0rWw9RyyspKcHT8Ft26NUVgYAo8PEpekvfVV5/h7NnJMDdvDB0dNfz6a3+oqirj\n/v05uH/vWzy4H4LFi28jLi672HMpb9+WOEUR+fgxHq5fD0FuboXiNhk0CKMPHcLUCxck17p2NUKz\nZjowNi55ieM7TTp3hoa+PgzatZO7vQfr1uHBunVS1zX19aGkqoqoJ09gNmECWvTvj7ZfflmsjOMv\nv2Bfq1Z4LeNEpAe//47TQ4fiyY4dKMjKgkgggH6bNuVeNaKmo4PxR4+i27x55XquLI82b8bRnj0l\nq1mYaiZPj1/WFwBnsJH5JyElJZdu3AgkobD0BFo9elgRwCM3t2giIuJnZdGeFi1onbYhTRh5nPj8\n96tO8lJTaYuGBm3V1qaC7Oxi9Vj17Fn0kvPs2QrFe/q0D5mbH6THjyMq9Lwsb2/elNpqnx0fL0mS\nlZuSIvXM+QkTiAfQ8wMHSqzzwbp1tFFJifyvlrwT1+fMGdplZESB168TEVF6eHitymt+tHdv4gH0\n9uZNIir6dz07ciQ583g1HFndhupczcI6cyYjI59GjDhDq1Y5EhFRampese3reamptEVTs+gA5EmT\nij1byOfT8X796OSQIVJLC9/evEk3f/yR8tLSKhTX3LnXCODRrl1PSy0XfPcuOW/YILNzDL5zhy7P\nmEHp4eEUev8+8QD6p1OnYmXCnZ2JB9D2Bg2knk8LDaUtmpq0TVeX4r29ZcYhyMuT41PJJiospLe3\nbsmVUVLRMiIji6VJiHj0iHgA7WvdutpjqU/k7czL3DTEcdw9AEYfXgJAANYR0c3y/BbA4/Ekf7a0\ntITlRxnjPnVBQalo0kQL+vqaNR2KhKtrJA4d8sCmTZYwN5edeS8kJA337oXB1zcJu3Z9CQMDTRgY\nvARc9vwAABRiSURBVP8cmgYGmHHlCq589RUaf7RVXkVdHQvc3Eqs12zcOJiNGwcAeHvjBhxXrsSX\nu3ejw6Syd7oCwP79ozBlSkeMHl36VInD0qVIDw1F8169JO196MXBgwh2cIBx37747Ouv0XLAALQe\nNqxYmZYDBqD/qlVo2q2b1PMigQAkEkGvTRvcXrwYXefNQ68ffpAqp6pZuX/7F4cOwXHlSnT59ltM\nPnu2UnWVl56JCfRM3u/ebTVkCKZduoTG5tWbJqKuc3Fxgct/576Wizw9fllfYCPzSvPyiiMlpY3U\nrdu/NR1KMV9/bUcAj/73v4dllr158y29eqWY4+xKcm/NGuIB5LhqlcLr9rt4kW4vW0YFOSVnjkx6\n84Yebdkic526PLITEsjz33+JB9CZESMqXE9pIp88oQPt29PL48erpH6m+kHOkblCliZyHOcMYBUR\nvSylDCmirfoqIiIDffocw+eft8bly2UkRKoi/IwM2M+di+Z9+mDIfy/xQkLScP68L5Ys6V0l2//L\nQ8jnI+z+fZh+8UWlR7A1RSwSwf/yZbQcOBB6LVvK/ZxIICj31nymfpB3aWKlOnOO4yYBOAigMYAM\nAD5ENFpGWdaZVwE/vyQMHHgCU6Z0xMmTEytVV6SrK04NGYIGLVtiZVSUgiIs4rBsGTKjojDt4sU6\n2xGXJeDaNQAoMYuiLC9exGLkSBv88ENP7NgxXGa5i1OmINDeHrPv3UObL76odKxM3SFvZ16pRFtE\nZA+AnSxbg9LS8pGdXYDIyIxK19Vq8GBMOXcOjRQ8x0lE8Dl5EoV5eciKjkYjMzOF1l8b5KWk4NLU\nqQCA1cnJ0GrUSK7nYmOzkJHBx9u3qaWW09DTg7KqaqU2GdVHia9fIzs+Hu1GjqzpUGoc2wFaDbZu\nfQw9PQ0sW9anSuoPDk5Fs2a6Mg9OloebWzQ6dGhcJS9fX509i6Dbt2Hx9dfoMHEi4r29EXTrFvqv\nXAk1HR2Ft1cZia9fw2nVKnw2cya6L1gg93NEBMeVKwGOw8g9e8BxHAqys3F52jQYtG+PMYcOyXz2\n1asEtGtnAG3t0v/9hAUFUFFX3OHb9cGuJk2Ql5KCH1+9glGXLjUdTpWolh2g5fnCJ/oCNCwsjQAe\nATzKz5fO6Fcb2NsHEMCjESOq5niw3c2bEw+g0Hv3iOh9giyPI0fI98IF+svAoEKH9QqFItq3z02h\na8gdV62SrBUv787KjyX6+REPoB0NGyomOEaKw08/0Znhwyv1Yrq2AztpqHYwNdXHrl0j0LChRqkH\nNdSkNm300aJFA/TqVfYBBhUx7uhRJPj4SA4v7v/rr2jQogXMJ06Ex+HDyE9LQ/Kb8h9n4+QUihUr\nHNHWtAFCwlZKrvucPg3/y5cx7t9/S3xhmBoUBIN27cApSW+AHvjbbwi5cwcqGhrQalyxQ5DfMezc\nGbPu3oVOU/kO52bKb/SBA5V6/vn+/Yjz8MDYI0dkHnpSZ8jT4yviC5/oyFxR/v77KQ0YYE3h4ek1\nHYpCiQoLKdLVtcQ85GXJzOTT4s9XE49TosfbtkmuWw8YIHPX6Lu84vfWrq1QvB/mTq+NSosvPTyc\nAuztP8lDkWXZbWxMPIDCnZ1rOhSZwFLg1i92dgF49iwaL1/G1XQoCqWkogKTQYOgpFL+31oaNFDH\nglntARIjL+V9Wt3xx4//v707D46yzPMA/v1xJCAQBGJACDfIXcDIBLk2QXBgYZCjLBCsooB1XRxF\n2BV0HKzdTMadovCgHA7dVQnCykBFLFiGQ2CSDkpCBCUQIGIMJGQ4ApuDQ0gnnf7tH90wIN2dvt/O\nm++niqrut5/3fX6VkF8/7/M+ByZ/8AEGzJr1wDkPxcYCImgRF+dzfXuXLMF/Nm+O4q++8uk8VUXx\nV1+h5tYtn+t0p/rmTRTu2wd7be3dY9dKSvD2I49g01Out3lLmzkTW6dNww87fZrrZ2rPbNmCyR9+\niK6JiUaHEjhvMn4w/oEt84CcPVuuaWmn2KpyoaygQGttNq/L+3MXoKqaNmuWJgOav327T+flrF6t\nyYBuX7DAr3pd2fH885oM6KG337577Gp+vqY0baprBwxwec7BP/5RPxo+XCuKgveMgUIPbJmbS/fu\nbfDMM/0hUvdD7XC4cfEiPhg0CLt+8xuP5a7m52PjuHE4lZYWslja9ur1wBrinvhzFwAA0z79FIsK\nCtB3qm/j+WP79UOLuDiX0/z91XnECLTu2vW+PUJj+/bFK4WF+KesLJfnjHnjDTx/+DAe7to1aHFQ\n5ODQRPLL33Jy8MkTTyC2b1+8lJ/vttzh99/Hl0uWoM/UqXh2u/mmJNRWV0MaN/bpy4TIF2GZAeoL\nJnPzKcnKQusuXTxOMbdVVeH4xo3oOWFCxLcIbVVVaBwV5XKUiys3S0uxbsAAtO7cGf/i49ZtkeKb\nNWtQWVSE8StW+H3HQqHlbTJnNwu5VHH2LFLHjME3a9e6LdN55Mg61wpp0qwZHn/hhZAn8trqaty4\ndMnv88sLC/F2XJzbh4eu2GtqYKuqgvXGjboL+8hmtaLIYrnvAWco7Fu6FNnvvosrJ0+GtB4KPSZz\ncqkkOxvnv/4ap7ZuNToUjy7n5uL4xo3YOmMGVsXH+zzS5I7a6mrUWq0+JeaY+HgsKSrCwtxcnD1w\nAIX79vlVtysHfvvbu7sKhdKMzz7DhFWr0H7w4JDWQ6HH+ypyaeCzz0JEED9ihNGheJQ2cybKCwrQ\nNTERjZo29XsRr0f69cO/lpQgyouJI3abDT/u3YsuY8bgodhYVFVW4n+ca4MsLS31a7JRbXU1ygsL\n8YhzrfdHhw5Fi7i4kE9R7+9cT4bqP7bMyaVGjRtj0Jw5aNO9+33HLx49ioI9ewK+/vlDh3DhyJGA\nr9N3+nREx8Sg+7hxeL2iAh2HDfP7Wi3i4rz6MshZvRp/njIFe195BQAQHRODIfPnY/DcuWjWpo1f\ndf9l4UKs698fJz77DAAweO5cLC0tRZ8pUwAApXl52PXSSx7317xWUoLs995DVWXgi65R/cNkTj75\ndOxYbJ40CVc9jGCpy09XrmBDYiJSR4/2e6PmOzoMGQLr9esotljCtrRup4QEtO3dG12dyxNIo0Z4\n+uOPMTU11e9RLTHx8WgcFeV2MlP2O+/g6Lp1OPrhh26vkb58Ofa9+iqyV61y+fmFI0fYN25i7GYh\nnwyZPx+V587dtz2Yr5q1aYPekyahafPmASfggbNmoelDDwXUIvdVl1GjsOiHH4J6zbEpKUj6/e/d\nziMYuWwZomJiMOzFF91eY/Dcubh19Sr6zZjxwGc3Ll7EJyNGoEl0NF4rL+fqiybEoYn11MGDxRg4\nMO6+fTbp746lpuLEpk2Y8tFHaNuzp6GxfPvtRdy+bcPo0f5/AQbKVlWFzZMnI7p1a8zcti2gyWf7\nli5Fwe7deG73bjzcrVvwgiSXODTRxD7//DQSEzdgzpxtYamvJCsLtysqwlJXsGS/8w6KMjJQsGuX\n2zLXSkpwLj09pHFYrTaMHp2KxMQNuHDhekjr8qRJs2aY+9e/YtYXXwQ8i7goIwP/l5+P8sLCIEVH\nwcBkXg/16dMO3bs/jJEjvd9D0l/5X3yB9aNGYdvs2SGvK5ju3AV6Glny51//GhvHjcO5jIy7xzaO\nG4f/HjYs4L78O6Kjm2D69L6YOLGX4XuoBsvsnTsxLzOT29dFGPaZ10ODBrXH2bOLw1JX2169ENO5\nMzolhGaXpFCZun49ijIz0c/D0LueEyYAImjXuzcA4Mcvv3S01EVQVVmJqBYtghLL5s3mGv7XqmNH\ntOoYmrXvyX+Bbui8EsAUAFYAhQDmq6rLe0n2mVOkK7JYsHH8ePSZMgWznJszExktLGuziMh4AOmq\naheRFXAs1fiGm7JM5hTxrNevI6pVq4hZnZIoLA9AVfWAqtqdbw8D8LxQB1GEi46JYSKneimYD0AX\nAAh8aiBRhDh58gpat16B+fN3+H2NMzt3YmVsLI6tXx/EyIgeVOcDUBHZD6D9vYcAKIDlqrrTWWY5\ngBpV3ezpWsnJyXdfJyUlIck5g44oElVU3MaNG1YUF/s/Pf7q6dO4XVaG0ry8IEbmPev169i3bBm6\nJSVhUD0bkdRQWSwWWCwWn88LeNKQiMwD8M8AnlRVq4dy7DOngFz67jucy8hAwssvB20GY2ZKCtRu\nR5KzoVGSlYXqmzfR81e/AgAUFJTh0UdboWXLKL+ur3Y7SrKy0CkhAY2j/LtGIL7fsQNbp01Du8ce\nw8tnzoS9fgqct33mge7rORHAKQDtvCgb4E541NB9lJCgyYAe37QpKNe7VVamyYAmA3qztFRtVqv+\nISpKk0W0srg4KHUE27mMDH23Y0fNWbPGq/I2q1Uz33pLz6an33c8Z/VqXd2nj144ejQUYVIQwcs9\nQAMdZ74aQBSA/c6HRodV1fOmkER+Gr5kCc5s344ePmwg4Unztm0xNTUV9trauwtcDZ43D1Xl5WjZ\noUNQ6gi2K6dO4cbFi7h09KhX5RtHReEfli9/4Pi59HSUnTmDy8eOoePjjwc7TDIA12YhqkfUbkfx\nwYPo+MtfBjSp6VZZGS7k5KDXxIleb5NHxuAeoEREJsCFtoiIGhAmcyIiE2Ayp7DL+dOfsG32bFRd\nuxayOnI3bMAfoqKQu2FDyOogiiRM5hR22e+9h5NbtuBiEPYAdefa+fOw19Tg2vnzIauDKJLwASiF\nXUl2NkpPnMDjL7wQsnVQ1G5HaV4e2g8axNEaVK9xNAs1SLfLy3EqLQ0DZs5E8zZtjA6HKGAczUIN\nUmZKCnYtXIjMlBSjQwHg6O75Zu3aoO1cROQOdxoiU+k7fTouHzuGftOnGx0KAMfmx6fT0lDz008Y\n9dprRodDJsZkTqbSLTER8zIzjQ7jrkHPPYeqigr0njzZ6FDI5NhnTkQUwdhnTkTUgDCZN2B2mw0f\nDx+O/xo6FLaqKqPDIaIAsM+8AbNZrbianw+7zYaaW7fQpFkzo0MiIj+xz7yBqywuhtbWok2PHkaH\nQkQucNIQEZEJ8AEoEVEDwmRORGQCASVzEUkRkeMickxE9opIZG6cSERkcgH1mYtIS1W96Xy9CEB/\nVX3RTVn2mRMR+SgsfeZ3ErlTCwD2QK5HRET+CXicuYi8BWAugEoAYwOOiIiIfFZnN4uI7AfQ/t5D\nABTAclXdeU+51wE0V9VkN9dhNwsRkY+87Waps2Wuqk95WedmALsBJLsrkJz894+SkpKQlJTk5aWJ\niBoGi8UCi8Xi83mBPgDtpao/Ol8vAjBGVWe6KcuWOfntyLp1KM3Lw8RVq7jsADUoYZkBKiKfA3gM\njgefxQAWquolN2WZzMlvK9u1w+3yciw4dAidR440OhyisOF0fjKVgj17UHbmDIYvXhyyTaCJIhGT\nORGRCXBtFiKiBoTJnIjIBJjMiYhMgMmciMgEmMyJiEyAyZyIyASYzImITIDJnIjIBJjMiYhMgMmc\niMgEmMyJiEyAyZyIyASYzImITIDJnIjIBJjMiYhMgMmciMgEmMyJiEyAyZyIyASCksxF5FURsYtI\n22Bcj4iIfBNwMheReABPASgOPBzjWSwWo0PwCuMMnvoQI8A4g62+xOmtYLTMVwFYFoTrRIT68gtm\nnMFTH2IEGGew1Zc4vRVQMheRpwGUqGpekOIhIiI/NKmrgIjsB9D+3kMAFMCbAH4HRxfLvZ8REVGY\niar6d6LIQAAHANyCI4nHA7gAIEFVr7go719FREQNnKrW2VD2O5k/cCGRcwB+oaoVQbkgERF5LZjj\nzBXsZiEiMkTQWuZERGQcQ2aARvokIxFJEZHjInJMRPaKSAejY/o5EVkpIvkikisi20QkxuiYXBGR\nZ0TkpIjUisgvjI7n50Rkooh8LyI/iMjrRsfjioh8IiKlInLC6Fg8EZF4EUkXkVMikicirxgd08+J\nSLSI5Dj/tvNE5D+MjskTEWkkIt+JyP/WVTbsybyeTDJaqaqDVXUogF0AIvEXvg/AAFUdAqAAwBsG\nx+NOHoDpADKNDuTnRKQRgDUAJgAYAGC2iPQ1NiqXUuGIMdLZAPybqg4AMALAS5H281RVK4Cxzr/t\nIQD+UUQSDA7Lk8UATntT0IiWecRPMlLVm/e8bQHAblQs7qjqAVW9E9dhOEYTRRxVPaOqBYjM5ykJ\nAApUtVhVawBsATDV4JgeoKpfA4j4gQWqellVc52vbwLIB9DJ2KgepKq3nC+j4RieHZF9zc6G7yQA\nH3tTPqzJvD5NMhKRt0TkPIA5AP7d6HjqsADAHqODqIc6ASi55/3fEIHJpz4SkW5wtHxzjI3kQc6u\ni2MALgPYr6pHjI7JjTsNX6++bOqcNOSr+jLJyEOcy1V1p6q+CeBNZz/qIgDJkRajs8xyADWqujnc\n8d0Nyos4qeEQkZYAPgew+Gd3uRHBeUc71PmcabuI9FdVr7oywkVEJgMoVdVcEUmCF7ky6MlcVZ9y\nddw5yagbgOMicmeS0bci4nKSUai5i9OFzQB2w4BkXleMIjIPjtuwJ8MSkBs+/CwjzQUAXe55f2fi\nG/lJRJrAkcg3qeoOo+PxRFWvi0gGgInwsl86jEYBeFpEJgFoDqCViGxU1bnuTghbN4uqnlTVDqra\nQ1W7w3FLO9SIRF4XEel1z9tpcPT9RRQRmQjHLdjTzoc69UGk9ZsfAdBLRLqKSBSAZwHUOWrAIILI\n+/m5sh7AaVV93+hAXBGRWBFp7XzdHI6egu+NjepBqvo7Ve2iqj3g+H+Z7imRA8ZuThHJk4xWiMgJ\nEckFMB6OJ8qRZjWAlgD2O4curTM6IFdEZJqIlAB4AsBfRCRi+vZVtRbAy3CMDDoFYIuqRuIX92YA\nWQAeE5HzIjLf6JhcEZFRAJ4D8KRz6N93zkZHJHkUQIbzbzsHwJequtvgmIKCk4aIiEyA28YREZkA\nkzkRkQkwmRMRmQCTORGRCTCZExGZAJM5EZEJMJkTEZkAkzkRkQn8P8OYDKiXrHefAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111fe2898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[:,0], X[:,1], c=y, s=5, lw=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#{y == 1} = 251\n",
      "#{y == -1} = 249\n"
     ]
    }
   ],
   "source": [
    "print('#{y == 1} = %d\\n#{y == -1} = %d' % ((y == 1).sum(), (y == -1).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing a basic Logistic Regressor class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last time, we wrote a very simple gradient descent method. I've re-written it here in a slightly more general form so that it cooperates with what we'll be defining below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "code_folding": [],
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def vanillaGradientDescent(funObj, w, X, y, \n",
    "                           objectiveParms={}, \n",
    "                           gdParms={}):\n",
    "    # optimality tolerance\n",
    "    optTol = gdParms.get('optTol')\n",
    "    if optTol is None:\n",
    "        optTol = 1e-2\n",
    "    maxEvals = gdParms.get('maxEvals')\n",
    "    if maxEvals is None:\n",
    "        maxEvals = 500\n",
    "    verbose = gdParms.get('verbose')\n",
    "    if verbose is None:\n",
    "        verbose = False\n",
    "    # initial objective and gradient value\n",
    "    f,g = funObj(w, X, y, **objectiveParms)\n",
    "    # this took a single function evaluation\n",
    "    funEvals = 1\n",
    "    \n",
    "    # initialize the step size\n",
    "    alpha = gdParms.get('alpha')\n",
    "    if alpha is None:\n",
    "        alpha = .01\n",
    "    \n",
    "    # Gradient Descent loop\n",
    "    while True:\n",
    "        # Step in search direction\n",
    "        w = w - alpha*g\n",
    "        # Update objective and gradient values\n",
    "        f, g = funObj(w, X, y, **objectiveParms)\n",
    "        # This costs one more function evaluation\n",
    "        funEvals += 1\n",
    "        # Test terminate conditions — largest magnitude of derivative\n",
    "        optCond = np.linalg.norm(g, np.inf)\n",
    "        if verbose:\n",
    "            print('%6d %15.5e %15.5e %15.5e' % (funEvals, alpha, f, optCond))\n",
    "        if optCond < optTol:\n",
    "            if verbose:\n",
    "                print('Problem solved up to optimality tolerance')\n",
    "            break\n",
    "        if funEvals >= maxEvals:\n",
    "            if verbose:\n",
    "                print('At maximum number of function evaluations')\n",
    "            break\n",
    "    return (w, f, funEvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A black box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also be using a function [translated and modified from Mark Schmidt's Matlab version] that can be — for now — treated as a black box. This function is not available, to my knowledge, on his website; however, you can find `minFunc` [through his website](https://www.cs.ubc.ca/~schmidtm/Software/minFunc.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "code_folding": [
     1
    ],
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# checks if a number v is in the appropriate real range\n",
    "def isLegal(v):\n",
    "    # assumes v is a scalar or vector\n",
    "    return any(np.imag(x) for x in [v]) and any(np.isnan(x) for x in [v]) and any(np.isinf(x) for x in [v])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "code_folding": [],
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def findMin(funObj, w, X, y, \n",
    "            objectiveParms={}, \n",
    "            gdParms={}):\n",
    "    \"\"\"\n",
    "    findMin computes the optimal parameter vector w to \n",
    "    funObj given the data X, y and execution parameters \n",
    "    objectiveParms & gdParms using gradient descent with \n",
    "    backtracking via Armijo rule, setting alpha via cubic\n",
    "    Hermite interpolation (I think)\n",
    "    \"\"\"\n",
    "    optTol = gdParms.get('optTol')\n",
    "    if optTol is None:\n",
    "        optTol = 1e-2\n",
    "    maxEvals = gdParms.get('maxEvals')\n",
    "    if maxEvals is None:\n",
    "        maxEvals = 500\n",
    "    verbose = gdParms.get('verbose')\n",
    "    if verbose is None:\n",
    "        verbose = False\n",
    "    gamma = gdParms.get('gamma')\n",
    "    if gamma is None:\n",
    "        gamma = 1e-4\n",
    "    \n",
    "    f,g = funObj(w, X, y, **objectiveParms)\n",
    "    funEvals = 1\n",
    "    \n",
    "    alpha = 1\n",
    "    while True:\n",
    "        # Line-search to find an acceptable value of alpha\n",
    "        w_new = w - alpha*g\n",
    "        f_new, g_new = funObj(w_new, X, y, **objectiveParms)\n",
    "        funEvals += 1\n",
    "        \n",
    "        gg = g.T @ g\n",
    "        while f_new > f - gamma*alpha*gg:\n",
    "            if verbose:\n",
    "                print('Backtracking...')\n",
    "            alpha = alpha**2 * gg/(2*(f_new - f + alpha*gg))\n",
    "            w_new = w - alpha*g\n",
    "            f_new, g_new = funObj(w_new, X, y, **objectiveParms)\n",
    "            funEvals += 1\n",
    "        # Update step size for next iteration\n",
    "        dg = g_new - g\n",
    "        alpha = -alpha*(dg.T @ g)/(dg.T @ dg)\n",
    "        # Sanity check on step-size\n",
    "        if (not isLegal(alpha)) or (alpha < 1e-10) or (alpha > 1e10):\n",
    "            alpha = 1\n",
    "        # Update parameters/function/gradient\n",
    "        w = w_new\n",
    "        f = f_new\n",
    "        g = g_new\n",
    "        # Test terminate conditions\n",
    "        optCond = np.linalg.norm(g, np.inf)\n",
    "        if verbose:\n",
    "            print('%6d %15.5e %15.5e %15.5e' % (funEvals, alpha, f, optCond))\n",
    "        if optCond < optTol:\n",
    "            if verbose:\n",
    "                print('Problem solved up to optimality tolerance')\n",
    "            break\n",
    "        if funEvals >= maxEvals:\n",
    "            if verbose:\n",
    "                print('At maximum number of function evaluations')\n",
    "            break\n",
    "    return (w, f, funEvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regressor shell class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we construct a class called `LogisticRegressor` which has `fit` and `predict` methods. One major difference between this and `sklearn`'s `LogisticRegression` class object is that I more easily know how to alter these moving parts without breaking the whole thing. This will let us alter the gradient descent method and the objective function so that we can compare run time, iteration cost, *etc.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The logistic regression objective function is given by \n",
    "$$\n",
    "f(w; X,y) := \\sum_{i=1}^n \\log (1 + \\exp(-y^i w^T x^i)) + \\frac{\\lambda}{2} \\|w\\|_2^2\n",
    "$$\n",
    "and its gradient is given by \n",
    "$$\n",
    "\\nabla_w f(w; X,y) = g(w; X,y) = X^T r + \\lambda w, \\quad r_i := -y^i \\sigma(-y^i w^T x^i)\n",
    "$$\n",
    "where $\\sigma(x) := (1 + e^{-x})^{-1}$ is the *logistic function*.\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/8/88/Logistic-curve.svg/320px-Logistic-curve.svg.png\" width=\"300px\" height=\"200px\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "code_folding": [],
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class LogisticRegressor:\n",
    "    \"\"\"\n",
    "    LogisticRegressor takes data [X, y] and optimization parameters and computes a \n",
    "    logistic regression linear model to perform binary classification.\n",
    "    X: an n-by-d numpy.ndarray object whose rows are samples/observations of d features\n",
    "    y: an n-by-1 numpy.ndarray vector whose elements lie in the set {-1, +1}\n",
    "    objectiveParms: A dictionary of parameters to pass to the objective function \n",
    "                    (e.g., an L2-regularization parameter lam)\n",
    "    gdParms: A dictionary of parameters to pass to the gradient descent method\n",
    "             (e.g., the optimality tolerance, the gradient descent step size, etc.)\n",
    "    kwargs:\n",
    "        intercept: whether to fit a linear model a.T * X or affine model a.T * X + a0\n",
    "        maxIter: maximum number of iterations before loop breaks\n",
    "        verbose: whether to print progress during the optimization\n",
    "        objective: the objective function to use in the optimization\n",
    "        gd: the gradient descent method to use in the optimization\n",
    "    \"\"\"\n",
    "    def __init__(self,objectiveParms={},gdParms={},**kwargs):\n",
    "        # # # Optional Parameters # # #\n",
    "        # Whether to include the intercept term \n",
    "        # (default: yes)\n",
    "        self.intercept = kwargs.get('intercept')\n",
    "        # Whether to print out information along the way \n",
    "        # (default: yes)\n",
    "        self.verbose = gdParms.get('verbose')\n",
    "        if self.verbose is None:\n",
    "            self.verbose = False\n",
    "            gdParms['verbose'] = False\n",
    "        # Objective function for logistic regression \n",
    "        # (default: L2 regularized LR)\n",
    "        self.objective = kwargs.get('objective')\n",
    "        self.objectiveParms = objectiveParms\n",
    "        if self.objective is None:\n",
    "            self.objective = LogisticRegressor.simpleObjective\n",
    "            if (objectiveParms is None) or (objectiveParms.get('lam') is None):\n",
    "                self.objectiveParms = {'lam': 1}\n",
    "        # Gradient descent algorithm\n",
    "        # (default: vanillaGradientDescent)\n",
    "        self.gd = kwargs.get('gd')\n",
    "        self.gdParms = gdParms\n",
    "        if self.gd is None:\n",
    "            self.gd = vanillaGradientDescent\n",
    "            if (gdParms is None) or (gdParms.get('alpha') is None):\n",
    "                self.gdParms = gdParms = {'alpha': .01}\n",
    "        return\n",
    "    def fit(self, X, y, w=None):\n",
    "        \"\"\"\n",
    "        LogisticRegressor.fit fits the objective function `objective` to the \n",
    "        data `[X, y]` using `w` as the initial starting point for the parameters\n",
    "        and lam as the regularization parameter.\n",
    "        \"\"\"\n",
    "        # # # Necessary parameters # # #\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        self.n, self.d = self.X.shape\n",
    "        self.w = w\n",
    "        \n",
    "        if self.w is None:\n",
    "            self.w = np.zeros((self.d, 1))\n",
    "        # Whether to fit intercept\n",
    "        if self.intercept:\n",
    "            self.X = np.hstack((np.ones((self.n, 1)), self.X))\n",
    "            self.w = np.vstack((np.ones((1,1)), self.w))\n",
    "        # # # Gradient descent algorithm # # #\n",
    "        if self.verbose:\n",
    "            # start time\n",
    "            st = timeit.time.clock()\n",
    "        w, oM, fE = self.gd(self.objective, self.w, \n",
    "                            self.X, self.y, \n",
    "                            self.objectiveParms,\n",
    "                            self.gdParms)\n",
    "        if self.verbose:\n",
    "            et = timeit.time.clock()\n",
    "            print('total elapsed time: %15.5g s' % (et-st))\n",
    "        self.w = w\n",
    "        self.objMin = oM\n",
    "        self.funEvals = fE\n",
    "        return\n",
    "    def predict(self, Xhat):\n",
    "        \"\"\"\n",
    "        LogisticRegressor.predict computes the predicted values `yhat` from the\n",
    "        data `Xhat` and the parameters `w` computed from `fit`.\n",
    "        \"\"\"\n",
    "        t,d = Xhat.shape\n",
    "        if self.intercept:\n",
    "            np.hstack((np.ones((t,1)), Xhat))\n",
    "        return np.sign(np.dot(Xhat, self.w))\n",
    "    def simpleObjective(w,X,y,lam):\n",
    "        \"\"\"\n",
    "        LogisticRegressor.simpleObjective returns the objective value and its \n",
    "        gradient given the parameters `w`, data `[X, y]` and regularization\n",
    "        parameter `lam`. It does this without regard to number of computations\n",
    "        required for matrix multiplication, etc. (hence, 'simple')\n",
    "        \"\"\"\n",
    "        yXw = y*(X @ w)\n",
    "        nll = np.sum(np.log(1 + np.exp(-yXw))) + .5*lam*np.dot(w.T, w)\n",
    "        sigmoid = 1/(1+np.exp(-yXw))\n",
    "        g = -np.dot(X.T, y*(1-sigmoid)) + lam*w\n",
    "        return (nll, g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "LR_vanilla1 = LogisticRegressor(gdParms={'alpha': 1e-1})\n",
    "LR_vanilla2 = LogisticRegressor(gdParms={'alpha': 1e-2})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 40.4 ms per loop\n",
      "100 loops, best of 3: 2.39 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit LR_vanilla1.fit(X,y)\n",
    "%timeit LR_vanilla2.fit(X,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "500"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LR_vanilla1.funEvals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LR_vanilla2.funEvals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 2 ms per loop\n"
     ]
    }
   ],
   "source": [
    "# Fill in the missing argument(s) in the parentheses so that the gradient descent method used is findMin\n",
    "LR_findMin = LogisticRegressor(gd=findMin) \n",
    "%timeit LR_findMin.fit(X,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Number of function evaluations\n",
    "LR_findMin.funEvals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing Progress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "code_folding": [
     0
    ],
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def binaryClassifier2DPlot(model):\n",
    "    \"\"\"\n",
    "    binaryClassifier2DPlot makes a plot of the two dimensional data X\n",
    "    used by model, coloured by class membership y. The background of the \n",
    "    plot is coloured according to the decision function computed by \n",
    "    fitting the model to the data [X, y].\n",
    "    \"\"\"\n",
    "    increment = 500\n",
    "    X = model.X\n",
    "    if model.intercept:\n",
    "        X = X[:, 1:]\n",
    "    # Get the label values\n",
    "    yu = np.unique(y).tolist()\n",
    "    # This sets the axis limits\n",
    "    plt.scatter(X[y.ravel() == yu[0], 0], X[y.ravel() == yu[0], 1], s=20, lw=1, alpha=.5, c='g', marker='+');\n",
    "    plt.hold(True)\n",
    "    plt.scatter(X[y.ravel() == yu[1], 0], X[y.ravel() == yu[1], 1], s=10, lw=0, alpha=.5, c='b', marker='o');\n",
    "    # Fetch axis limits\n",
    "    xLim = plt.xlim()\n",
    "    yLim = plt.ylim()\n",
    "    # Domain on which to compute the decision function\n",
    "    domain1 = np.linspace(*xLim, increment+1)\n",
    "    domain2 = np.linspace(*yLim, increment+1)\n",
    "    d1, d2 = np.meshgrid(domain1, domain2)\n",
    "    d12 = np.array([d1.ravel(), d2.ravel()]).T\n",
    "    # Compute the decision function\n",
    "    vals = model.predict(d12)\n",
    "    zData = vals.reshape(d1.shape)\n",
    "    # Set the colour map\n",
    "    if all(zData.ravel() == yu[0]):\n",
    "        cm = [0, .4, 0]\n",
    "    elif all(zData.ravel() == yu[1]):\n",
    "        cm = [0, 0, .5]\n",
    "    else:\n",
    "        cm = np.array([[0, .4, 0], [0, 0, .5]])\n",
    "    # Filled contour plot of the decision function\n",
    "    plt.contourf(d1,d2,zData+np.random.rand(*zData.shape)/1000,colors=cm)\n",
    "    # Plot the data [X, y] over top (it has been covered by the filled contours)\n",
    "    plt.scatter(X[y.ravel() == yu[1], 0], X[y.ravel() == yu[1], 1], s=10, lw=0, c=[0,.5,1], marker='o');\n",
    "    plt.scatter(X[y.ravel() == yu[0], 0], X[y.ravel() == yu[0], 1], s=20, lw=1, c=[0,1,0], marker='+');\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAAJZCAYAAACwSNHoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvX+4HEd55/t9hSXHPgZHCRsQOLZvbhZjWeyG7LNceJyN\nxvy0QVLI1c9sll9hb8SPJBAlTkAOOSMnNl6cxeQJGEs3xEm4YSVL2jg6MiaYtUYQExIWcCL5CCe7\nWWzIGpPEsljLXiSjun/UlKamprq7uru6p7vn+3me80zPTHdVn2P3V9966623RCkFQgghhBASxpJp\n3wAhhBBCSJugeSKEEEIIyQHNEyGEEEJIDmieCCGEEEJyQPNECCGEEJIDmidCCCGEkBzQPJFciMhH\nReS64fFqEfm69d3/EJGXT+/uCCGEkOqheSK5UEq9XSl1g/1RjHZFZKeIfFVEvisibyzZ1kYRuU9E\nTorIvTHujxAyu8TUFBF5rYh8TkSOi8j/FJFdIjIX615JPdA8kaZwP4C3A/hShLb+CcAtAN4foS1C\nCImpKc8C8BsAVgC4HMBFAG6O0C6pEZqnDiIivyIie53PfltEPjQ8frOILIrIt0Xkv4nIz1rnrRaR\nr4vINhF5VET+XkTebH1/u4hcH3AP/1pEPj8cXf29iPyOiJyTdL5S6qNKqUMAvuNpS0TkPcN7/QcR\n2S0i35vS1r1KqX0AHsm6T0JIO6lS51yyNEVEXjqMTB0Xka+IyOqUtnYrpT6tlPrfSqkTAP5fAFfm\n++3JtKF56ia7AVxjQsEisgTARgB/NPz+UQCvVUo9C8BbANwiIj9iXf9cAM8E8DwA/x7AR0Tkwpz3\n8F0A7wbwfQBeBuDlAN5R7NfBLwBYB+DfDO/pOIBbC7ZFCOkGTdA5iMjzABwEcL1SajmAXwawX0S+\nP7CJ1QAeyNsvmS40Tx1EKfUwgC8D+MnhR68AcFIp9cXh93crpb42PP4cgE9DGxPDKQC/oZT6rlLq\nbgBPALgs5z18WSn1l0rzMIBd0CJRhK0ArlNKPaKUOg3gegAbhmJJCJlBmqBzQ/4dgLuUUn867Ou/\nAPivAF6bdaGIvArAGwC8r0C/ZIrwH5/u8p8A/NTw+KcAfMJ8ISLXiMifi8g/ichxANcAeLZ17T8p\npc5Y758EcEGezkXkn4vIgog8IiKPA7jB6SMPlwD4YxF5TEQeA7AI4DSA5wxX//2vYWj+PQXbJ4S0\nk6nq3JBLAGwy+jTs60oAK0Tkxyx9OmJfJCIvhY6SrVdK/fcC/ZIpkpiDQlrPXgC/JSLPhx6ZvRQA\nRGQZgH3Qo6U/UUqdEZE/BiCR+/8o9Khws1LqSRF5F4D1Bdt6GMDPKKX+3PPd24c/hJDZY9o6BwBf\nB/CHSqmtCd8/0/1ARF4M4E4Ab1ZKDSq4J1IxjDx1FKXUPwI4DOB2AH+nlHpw+NWy4c8/DgXlGgCv\nruAWngng20Pj9EJkGBwRWSoi3wMtbstE5FwRMUK3E8CNInLx8Nx/JiLrUtpaIiLnAlgK4BnDtjhQ\nIKRj1KVzGZry/wFYKyKvHp73PcOE9OcltLUKwN0Afl4p9cmi90SmC81Tt/kEdB6ASaCEUuoJ6ATs\nvcMpsC0A/iSjndBaTvZ5vwzgp0Xk29DmZ3fGtZ+GDpu/bHj+kxjlJ/z28B4/LSInAHwewEtS2noD\ngKcAfATAjw3b2hX4OxBC2kUdOpeoKUqpbwD4CQDbAfwDgIeg9S/p39dt0NOHHxtO6f0vd0qPNB9R\nKkqNQ0IIIYSQmYCRJ0IIIYSQHNA8EUIIIYTkgOaJEEIIISQHNE+EEEIIITmobfm2iDAznZAZRClV\nRW2dWqF+ETKbJOlXrZEnpVTQz/z8fPC5sX7YJ/tsU39t6bNLdO2/Tdv6Y5/ss+7+0uC0HSGEEEJI\nDmieCCGEEEJy0Ejz1Ov12Cf7bFWfs/A7TqvPNjIL/21m4Xdkn93qM2Z/tVUYFxFVV1+EkGYgIlAd\nSRinfhEyW6TpVyMjT4QQQgghTSWaeRruJv1lETkQq01CCKkD6hchJA8xI0/vArAYsT1CCKkL6hch\nJJgo5klELgLwWgC/G6M9QgipC+oXISQvsSJPtwC4FgAzKgkhbYP6RQjJRentWUTkdQAeVUrdLyI9\nAIkra/r9/tnjXq+XuGxQtrZ+cQ4py67+tO+AJKDUfOJ3g8EAg8GgvpspCfWLkA4R8O9GLP0qXapA\nRG4E8O8APA3gPADPBPCflVJvdM4LXupL8ZlxaJwaTZr4uDS9VAH1i5COEPjvRiz9Kj1tp5TarpS6\nWCn1QwC2ALjXFR5CCGki1C9COsAUBtys80QIIYSQdjKlmYrSOU82SqnDAA7HbJMQQuqA+kVIy5hi\nigcjT4QQQghpF1POjaV5IoQQQkh7aMCiIponQgghhLSDBhgngOaJEEIIIW2gIcYJoHkihBBCSNNp\nkHECaJ4IIYQQ0mQaZpwAmidCCCGENJUGGieA5okQQgghTaShxgmgeSJNo8EPCyGEkJpo+L8FNE+E\nEEIIaQ4NN04AzRMhhBBCmkILjBNA80QIIYQQkguaJ0IIIYRMn5ZEnQCaJ9IkWvTgEEIIiUjL9J/m\niRBCCCHTo2XGCaB5IoQQQsi0aKFxAmieCCGEEDINWmqcAJonQgghhNRNi40TQPNECCGEEJILmidC\nCCGE1EfLo04AcE7ZBkTkXACfBbBs2N4+pdSOsu2SGaMDDxNpJ9QwQmqkI1pf2jwppb4jIlcppZ4U\nkWcAuE9E7lZK/WWE+yOEkEqhhhFSEx0xTkCkaTul1JPDw3OhDZmK0S4hhNQBNYyQiumQcQIimScR\nWSIiXwHwTQD3KKW+GKNdEoE1074BQpoPNazBUMPaT8eMExBh2g4AlFJnALxYRJ4F4E4RWamUWnTP\n6/f7Z497vR56vV6M7kkaawEcnPZNkFlhMBhgMBhM+zZyE6Jh1K8pQQ1rNy0yTnn0S5SKG50WkfcB\nOKmU+qDzuQrtS7ZK1HuaSdZAi45hAc0WoBY9YLOOUvPB54oIlFKteqB9Gkb9mgJt0zAySQN1PZZ+\nlZ62E5Fni8iFw+PzALwKwFfLtktKchDA1uHxVlB0CEmAGtZQqGHtpoHGKSYxpu1WAPgDEVkCbcb2\nKKU+GaFdEoOFad8AIY2HGtZkqGHto+PGCYhTquAIgB+NcC+kCjhaIyQValjDoYa1ixkwTgArjBMS\nxvyhad8BIYQUo64VizNinACaJ0LC6B+e9h0QQkgx1mafUpoZMk4AzRNpAlU9dDGiRfOHANXXx6rP\nCBQhpD7KRozWANg5PN4Zob0kZsw4ATRPpMvEiBbtuAqQvj6Wvn5PCCF1UDZiVMeKxRk0TgDNE+ki\nVUSL+qvLt0EIISHEjhhxxWJ0aJ5I96giWlRFxIlTgIQQH7EjRlVEnNZgZqNOAM0T6TJNjxYxCZ0Q\nkkaTI0Z1JKE3GJonMj2qXj7b1PwkJqET0g2q1rAm1riypxRnWL+ibAxMSCFmdcPPHVfpH9UfTS8S\nQtrHLGrYur5+nXH9YuSJ1FdAze7PToac0ZFL46cVCWkL09awuvufBrv64zlOM65fNE+k/rlrNxmy\nqdNrVTOrvzchsZm2hs1a9AmYef2ieZoVfCOjaY+empwMSQhpFtSw6THDq+qSoHmaFXwjs2mPnmZx\ntEYIKQY1bDrQOHmheeo6ISOz0NHTLMzrE0KaRSwNo37lh8YpEZqnrhMyMgsdPVWRV8CHkxCSRiwN\nm/G6RLmhNqdC8zQrlJmbn3ZeASGEFNUw6ld+aJwyoXmaFcrMzU87r4AQQorqDvUrHzROQdA8kXBm\nZWUJIaR7UL+yoXEKhuaJhNO0EdusFtckhOSnafrVtOlDGqdclDZPInKRiNwrIg+IyBER+YUYN0Zm\ngLIPKzfWJRGghpGp0KQEdhqn3MSIPD0NYJtS6goALwPwThF5YYR2CfGTtLHusqXTuiPSbqhhpD6S\nEthlbjr3Q+NUiNLmSSn1TaXU/cPjJwAcA/D8su2SkhQNCbvXNS20DOhtAcyGlNLX7zesB7Zv16+E\n5IAa1kC6rF++BPblu4EVT+jXOqFxKkzUnCcRuRTAjwD4i5jtkgIUDQm71zUptOxiNqZcthRYtUof\nr1rFCBQpDDWsIcyCfpkEdpkDztusj8/bXF8EisapFNHMk4hcAGAfgHcNR28kNmuc16RzitQ0ca/b\nVrCdOjEbU546DRw9qo+PHtXvCckJNaxiqF/jmAR2dRJ4ao8+fmqPfl81NE6lOSdGIyJyDrTofFwp\n9SdJ5/X7/bPHvV4PvV4vRvezw1roB868+jg4/NmJUWg4hKTr8rYTi/lD+Xbt3rcfOHCAxmnKDAYD\nDAaDad9GbkI0jPpVklnSrzXIt7rv+Bbg8bfWY5xIInn0S5RSpTsUkT8E8I9KqW0p56jQvmSrlL6n\nTrEG/vDzApIf0LwPb9J1RdsJIW30o/qjvCbSKJSaDz5XRKCUavwDnaVh1K8SdFW/0piWaQthxqNO\nsfSrtHkSkSsBfBbAEQBq+LNdKfUp5zyKT1nMA9nkBzOUpAd4/tB4CYL+6nwRKFI5XTNPIRpG/YpA\nl/QrCdcophnEaTDjxglokHnKcRMUn7KYUdS0RlOxYeSplXTNPIVA/YpA1/QrjSYaRBonAPH0ixXG\n28RB5zUm06oxkoRZSUcI6QazpF9N2wqGxik6NE9kejVG0kiaqsu7JQu3cCGk2zRRv9IMYp6VfzFW\nCdI4VQLN07RoytLZrBojVd1n1gOdZHrybsnCLVwIic+s61cWaf3mqT1Vtk4VjVNl0DxNi6YUb8uq\nMTKt+3RNT9KWLPb3ec4nhBSH+pWOr9+sGlZrcpwbAo1TpdA81U2eB6gujm8BHrlAv9r3MY0ic0mm\nx96Sxby3MWbLd77ZwoUQUo4QXahbw5qkX2n92tuymPc2a63zfVu4kEbB1XbTwrcawyxzbdIqjapW\njWSNinyr7XxlDIDJKJV9Xd5imyQTrrbLOHeW9atIgcsqmda9JPXrK2UAz2cHrfOLGCdGnRLhars2\nkDbacVdj2A9Vk7YTmNaqEd9qO180yY1IAZMRK0JIfvLoF6D1q2lbokxLv5L69UWUDjrnuxGovNA4\n1QLNU5Wkzbe7VXCbWlitzvuw85LSTI/PWPVXc5qOkJjk0a+dzvdNmWqq+x5CTY9rrmwDVeZvR+NU\nGzRPVZB3vt0ejTTJONVN6Mo4nzEyn/VXMzmckDKU0a+taF6NozoJTVD3abxtoJgg3nhonqqgaLJf\n141TkiDEWhln8ptYnoCQ4pTRL3N9F0kzNDES1O38pryrBGmcaofmqUryjsCKboRZlrpyE4wg+B70\ntJV0ofQPszxBEsuWTvsOSNugfo2TZGiM6Sm7Os7kOuU1YbNgnBqoXzRPVVLHCCxrhBLy8IWMckIF\nKmnpsi0IrqnpHx7/LK/xsSNXBuY9jdiwHti+Xb8SEkqX9Cu0rRD9cs8panp87a9FvtynWTBODdUv\nmqe2EvqwpglLSBumYm+oQPnOc0dlO67Shsc2Pf3Do0TwvMbHXYXHffFGLFsKrFqlj1etauQIjswg\ndekXkE/DQvTLbG5cxvRktR8S9ZsF49Rg/aJ5aitZYeKQYpxZbSzfDWz55TCBChEyIwimXpOv9ECW\n8UmLSJlrXeM1y9N3p04DR4/q46NH9XtCpk0d+gWEa1ge/TKro4uaniQNta/1/T1sZsE4AY3WLxbJ\nbDtZRdSSirXZn/vakDm92aZhXsKKzWUVpVux2l/oMjTS5CueWcU1XWPZ0mjCwyKZGedSv8KpSr+A\nYhqWpV9JRS7zRJqKFO60r5kV42TTQP1i5KntZO3e7Y6EtmFyhOVrw90zKjR5NO28NfAbpxDSVuQl\nRZbKruJrUIi4NA0asRFylqr0CyimYVn65TNOoSRFt/Ks4luRQzOpX5VC89RlfPP3lyF8VcjxLcC/\nWqZfQ0dWaecdxCgC1F89KisQWlogqQhm0vVl9rdraJIiITODmRqzzUUe/QKAK61970I0LEu/3Hp8\na5GvrIDv3rOKkZpr8mgY9atyaJ66iJvI6CY37gTwYGhbFTn+3tfGV8hlRYaMQbKjVaGRJXNN6Eis\nwUmKhHQeV6t8q9lC9WstdAQqNi+w7gfD420B9wKM506FrtDLo2HUr1qgeeoivr2SgPFRzwcz2qhi\nV/I1GJmd3kOT03a+UZVrkNzzQyJLO67KNxJrcJIiIZ3HjrYY3NVs09Iv0+ZlGNfYrcPPQu7FEFob\nalc/n4ZRv2qBCeNdReaA153UorMVxffPK5LcmJaH4CaMf+1C4NITo/dmOs8lLenbVBZPYtlSLTqG\nG28ME5SISYpdggnjGedSv+KwdimwcFpr0AKaoV+ujpoImG2cku4t6V7S+jPJ4UU0jPrlpVEJ4yLy\nMRF5VET+OkZ7pCTLd+tVJvftHl8NYo9yQklKijS1U3ykzeG//9XAYPg/7w4FXPYrYdGjtLB1Vh5A\n3pGY6YPCMxNQvxrI8t3Al07pV2NGmqBfd82N69ctczoKFhJBWkjoN8s4Afk0jPpVC7Gm7W4H8JpI\nbZEyyBxw3mZ9fN5m/bDbGDEps4GlMWfLd49/HhIqP3Ua+MgGLUBGBOYPTU7huSYpNGydlPe0b78e\nre3bn3wtENYHcwi6BvWrSaRp2LT1S50EPnW51q+n9uj35jzXqLlG6b6Eft17APzlCEI0jPpVG1HM\nk1LqzwAcj9EWKYm7PNckS9pCEFr00ocrbLZApM3h2+3t2w+8ZtlIBEzBTIMRgD2Pjj4LTYJMW7kX\nEnHK6oOrWDoH9athZGlYaC5T0o4JZfXr+BZg92/pV2Bk5OzzjUHbuD67X5u1SK/jlBVxon7VBhPG\n20Ro0uNxa3muwX7AQ0LMSSO7JGGzWcBohV9Seybi5K6WswVg00fHQ9BHjwKr+/6wddmaTnYfALC4\nONkHV7EQUpw8SdtJGpZnA16fhsXSLxNx8hk52yi9cr9+b/pd3ff3a7dF/WoF0RLGReQSAAtKqX+R\n8D0TLstSNPnRl2iZlKQYmlhuBCHtXs31dnu+hHA3GXzPo9o4+a7JqhYeo5r4po3AypVaiNwQ+Yb1\nWnh8380IXUwYp37VQBH9AvJX9Q7RsKL65WvL93ttXK+Nk3td2t9gV5/6VQOx9OucaHcUQL/fP3vc\n6/XQ6/Xq7L692GJgVp7kKVrpe2iTjFPS+fb3QLLwuMJllhibVX+7PMndbr7T5ucAb7gB+M51wLk3\nWHlRwyk51U9elVd2U+BlS7XwAFpkDhwYH8Ht2z/5WWxavEpmMBhgMBhM+zYqgfpVkDL6BaRrktvP\nwYzzzTlF9cuHLyl9735g3xxw28nRamc7SpX0N6B+TZU8+hUz8nQp9MjtRQnfc+RWlrwjN9vsJB0n\ntW/Osc/N078RCru9dX1/WQH7M/t4z6PAsVv1KOmKn5vc2iVPxfA8hIzOqhKIFowMOxp5uhTUr2op\nWzbAp0lp7ZtptIMp5yTh06+kvpO0dflu4OpjOsH8yi3pEayY+9VRv1JpWqmCTwD4PIAXiMjDIvKW\nGO0Sh7x7KdkP68GEzwH/3L29dUCRgnML0DVQXAHzJXSPmaLh8bKl2jgB+mF8/6snt3bJIm1Of9nS\n5O/37Qduvjn54a8q6ZI5CVOB+lUTRfaC82lYiH6Z89dmnJN2r3buU9pWKm6kChjlPR3uj1YMpk3X\n+ch6/pM0zKzKO3DAfx31KwosklknacXQYraXlufk+9zM/7uRJ/fcpNC1ex/mvT16c9sz4emk1XH9\n1cADH9YP4eXv0AmQeSJPvhGQGW2Z74D8eQFFC26GMkMjtzZB/RoSU8PyFKP0lSjIq18h6QjuZ2m5\nT+79wPn+vt3aOF25QU/j+e5pXd9/fZIOlNUw6lezIk+dIEb5/izybCBZpr0XJHzuW6mSVPPEd27S\nyNG9j7WY3PcJmCyG6W6v4n5vRlCbPjrcRNjKB+gfBg7d7r+fufMnR0BmtLVp4+g7+3tD1uip6q0P\nQutREeLSNg1Lamsb/ManrH6ZtkLuw9Uwd3sY3/24/ZkVgyZx3NXQtfCvqkvSoBgaRv2KBs2TIbax\nsSky7ZXUTkh7WTuPmwfYhJZX94Gbtvjbsx923+jMvo9tGBccg30P/dUjwfCVF3A3/v3OdaPvzfWA\nNli9hyb72rAeuPZa4MRwyxcjFEZMVq7UUSyD+d7gExfXQFUtEC1NtiRTpi0alke/fInVRfULmIyQ\nu/dhf2ZjDyDTdNjubw10wrj9vWugfNHzdetGx7YGJWmYW5YgS8OoX1HgtF3o0vwYFF2qm3S9L0HS\n3XfJt4Gmm9R43mZde+SmLcXuL+m+3HwnM7dvluPaK+iA8H3tDt0+bpzMdW5I+uabgZNP6mM3nGzE\nZN26sPB4g8PQdcNpu4xz6562a6uGFdUvc66JPMXWL/uztKTz0L+72767v6ete2nTaj4Ne/3rk8sS\nUMO8cNouFnmKrpXlwexTvCSN1HyjKvt3SRIe+4G3i9HlTeg0uNcl1WE5dPt4pAnw72vnhrLd5buv\n+Vn/de6IyxgnYDKJ0ghSWnh7xhIgSUupS8PcyErea9OiNEC4fgEjDatCv9zPzN/Tjqy7q/Fs4+RG\n0uy2dvW1Xp17g37v7ueZNq3mW8hilyXwpRhQwyqD5slQ9MHLw2XZp3hJEsckkUz6XZJMmKl5kiW6\nSaF69zrfyr5dfR0tCtoE2Eki923dsmmjvyZKWkj6498YX2WSNf9fVX4ABYxUQdUaZip8F8GnYXn1\nC/BrWKh+met995b0mT3QTEqHcJPJ3elTN/JuNGzvO/z3mKRhG9YDB0+G61foOXmhfgHgtF09xAqr\nx1rp4o6aQvsLCdmblS/u72wzuAS46i2T03aDS8an48x5hrnzdT6TYXERuGNvxg0NWbZU50/tGP4/\naIfDs2qexKyJ0oEQOqftMs6lfqW3FUvD8upXyHUh+vU3S4D/eGbUR9J59rSjMU5FNcxM582L1rA8\n+hV6TgjUr7Mw8lQHscLqRQ2XPdraZh1nJX4aUXDD1UnX2Ctf7N95wTreoYCPDH8Rd7Xd4NLx9mwj\n5SaCAzpkHTIKshPP5wV46XsmEyzNeT5iRpwYQidtI+a0YN5r1zivoVoEjJuakOuS9GsrRtGw9+wG\n/tN3R6v7XJ2zMTMNdsSpqIa999NauwD9+t5Pj76z9alKDaN+jUHzVCd1TA26rMVk6DkLNzSetXpv\nDZJ3DTchbZkDBvN6ZYz74JnpN9tMGVQfuP6zo4f2wguBr35VH4eGoe0cgx0K+MJN/gc/qeZULJJC\n6DMuQqQlTEu/7Nci+rVteJ1basA+P02/TP8Hl+qK4e455jzbSBl2YnLD8yIa9v5Xj6LmO5R+76NK\nDUubApxBDeO0XVdJCycb0sLXvuW6SdeYULi98uX4lsl7uOPtumbTxiM65OvbqgUYfW6vsvvFd2vR\nOXECuOVDxcLQ9nYvdsg5dNUfECf8bbfRwjA4p+0yzqV+laesfgGjPeSyVvAl6VfSfexQ4xqXVmDT\n6JevqGVeLdmwHnjnPuAjGya1IlTDYuuXua8WaVgs/aJ56jpufpO9zNYXQg8RrbSK5Xd5divfCf0w\nuw93/3DyDuK2UBw4UK4qrv2wp4lH1o7msUWi6mq/FUHzlHEu9SseRq/c1zzVyV3y6pe5D18bSTlU\njx2pX7+AdA2rwuS0UMOY80SSSVsqm1RGANBhaHsPJl/Okn2tLxfCJzwLGJ+Ss42UKY5pMPs1rVql\n94VatUrXMrGxi8hl4e7jlPZgp+1onjXfn5RrkEbV1X4JaSO+fKQF5zXJOKVpmCGvfpnr3TaTcqj+\n6APT0S8gWcOq0C9zPzOqYTRPXcTdzNe8t1njHJtkSZMwaWpS2aM0Q1YVXxdb6Ew42VeywAjFunWj\nB3JxcVTLxBCarJg3wTFtr7wskSiaazBD2xkQEoRvM3IfbjK5rWEb108aLkNe/QImt3dJSqJ/7Mh4\nUnid+gUka1hV+gXMrIbRPHUJ1zD5EiTdBExzbJIlMTw2iZkhq2OK1oeyR0muUBw4oB/IO+8cPfT2\nliuhS3PNNSdOlB8V+UTCt8VMXmZotEZIInn0yz42C1JsDXvlfh2BArL1KU99KHcgaRuv379hPCn8\n5ptnQ7+AmdQw5jxNizz1TnxbBKS15+4u7pvf9/G11cDXekBvx+gzO/KUNMef9buYnCc3CdzFnZO3\n3x84MKqY63tQr/8s8Os/Pt6OO+pz5+PdtsokU2blS3UI5jxlnDsL+gWEa5itQSEaVlS//u5y4OFN\n+fXL7ieJtG1lgPFyBEazAOpXA2HOU9vJs4mnW27AxZf4aEia3zev9vGlh4Hdv+UPR5v2Q3clN/dl\nRpL9w+MjnEO3T55vj4zcSNTSYcjaFYf5Q8DmTcD77tWik7SBpjvac3MJ3Pd5ScuXIqSLhGqYL9pt\ns8b5PkS/fFr2Q8e0fr1tbvRZiH5l3Zs9fejLczLGCRjfBsrWL8BvbuYP6d0SqF+tg5GnuslTrdc3\n2nLPz2rPHVH5RoG+lSQYfm9qpPjad/v2LQN2r7fJGumYEdiJEzoMbo/ggFHVcJvfeDlw/D+nlzXw\nbSBsV/0tu2IkZkXyBsLIU8a5XdYvIFzDkqJF7vlJK9lMG0l6ZkeT3OuBbP3y3aP7fZp+LQBY10/4\nEtlRKGA0oDRQvyqHkae2kqdar1t0zVfcLUvEkvZZOghg7dJJ4bH7kLn0ApkHMZ6U6ROZy5A8qjFz\n7Enz7GYjzAsv1O9XrRqNrvY8OmmctiwCv/mK0fkXXugXAt8GwqErRrKSNsuOAAlpOqEalhTtts1P\nmv4AyfoFAOLRrwcx2jdvDbIL/Nr35rYP53qXtcjWL18Uavt2rV+ucdqhgP9wNfWrJTDyNC1i5jyZ\nYnB5jNXy3bparp0fYJ/nEyU3qmT6eBCTo7sXIHnEdu4N2viYyFNofSUfZq+nl74XuPp7xs/PqmeS\nN2cgq90W1jwpAiNPGefOgn4B8XKefLWbiuiX0aGtnuvT9MstpJmlXweXAmtOh+sXkKxhRr9W94GP\nbtL73FEC7xdwAAAgAElEQVS/KoVFMsmINBHzJTfKnF7Oa5jP+Hu/zVM4Likkb0TKLmon/ZHAHD0I\n7PsicPk79Hx+nsre69aNHv4f/EE9MlvdB27bDOy5Y/L8mA9+qLC0rNpuEWieMs6lfuUjhn65Azib\nwbzOhbI1LC0B3ZgpX3HOX70c+Phifv0C/Bq27hZdD8rdIJj6VRmctiMj3Dl6G1+CpDqptxcA9OtC\nwnk7FPCZ9dnGyZ26c+uzqD4wuEQ/wFf8V/3ZsVv1/kxmxLZDAQ982PfbaU6dHg+Dm9D24T7g+0ct\n9ogptBjcjNY8IaQwMfTrbxKuG8wDn7o83Ti5fZjvbB17ENq0/dBwb7u8+gX4Nezw8Po775w8NybU\nr+hEMU8icrWIfFVE/kZEfjVGm6Qg7ugraUR3fAvwyAX61eQImMKYW4fHj1wA7PU8RHYuw4NIXk1j\nzpE+MLh0lKM0L3rkduq0XmkyGI4E7EJwSXPzp07rH7MKBQjfmdyQdm5anZNQYQmp4UIaBTWsIZTR\nL19V8d2/Ndqnzm7Tp18udnt/A+C2oQGbF+Cl7xmVHvDpF5CuYS5LqV9to7R5EpElAD4M4DUArgDw\nUyLywrLtkpzcgMmdxLNQJ8d3Bv8gRqOwy5C8VYFhYXiNK1rmO2CULO7uCv6Gi/RDuHLlaPS1uKgf\nXJO0uGnjeH/2Q3vH3nEDFbrlQVZCZFal3bIjQiZkNg5qWAOIoV/ASHfWIn27FXOuT798OnYQOlcT\n0Pr1hZtGemT0yybpOTfX2JGgEyf0arkQTaB+NYYYkaeXAPhbpdRDSqnTAHYD+IkI7ZI8PNs63orJ\n5EgfZjuD5bvHP0/as8nGl6dgJ5mbLRZMDoCpljuYH1XLdSvo3rF3ss6JMVC+h9YOdYdEra7/bPJ2\nB7Eq7aZRZLsFUgfUsGnzbIybl2nql29bmF19PUAzUSYz9eWrAj53vv85dzXMt5rY7O3pg/rVKGKY\np+cD+Lr1/hvDz0gd+Jb7JgmGjb2dwXmbRyO40GXIJtS9zdO/79plS0dz/GYJrvkMGH3mRpBWrkwW\nI988ftrI6H33jp9vk7TfXkxmeBPNhkMNmxa+PezSTI8hln4llUuwr9/VH9+sHNA5S8Ckhm3aOF5z\nyZBkPNwyA+vWheuXmTYEqF9TgAnjbSerFpSLESU36dINbydV4nXF7jKMV/T1JZ+blSPug+d+Bkwu\n582qY5JWldwIiz0q2/siYMkr9LFPpKqutMuETEJGZNWCcomtXzZGv+zrjXGytcqkFwCTn7ubACdp\nnU/D3HpQafpltq9yNYz6VRvnRGjj7wFcbL2/aPjZBP1+/+xxr9dDr9eL0D0BELYzODDKBQB0IuWP\nvQFY8Iwg0grfuXWgbjsJHH82gH+cvM6tspvF0aP+Srz79o+/N9hLeo1AmevN5zuu0ith9r4I2HhE\nt7PdEim73SpGbC4dHrENBgMMBoNp30ZegjSM+lUhrmlJIqZ+ufWd7Kk88+ou3V+yRBukDesnDcSZ\nMyMD5eoX4Ncwo1/msyT9MvuC7lD6u2WO0TLtUr9KkUe/Std5EpFnQK9ZeAWARwD8JYCfUkodc85j\nnZRpklR0LmmjzCy2AbhlThunHUqvQHHrQbn1WGxuvFG/2rVHAL/oJJFUkyRpO4PVfR12v/HG8Xor\nHEXlomt1nkI0jPo1ZWLr1xroCuVfOjWqc+fq1+/fkL0NCjCpYW7NpiRC9ct89qenRtOG1LDCNKbO\nk1LquwB+DsCnATwAYLdrnEgD8IXHQxMrfcdmNd7fXT4Sn9tOjp+jTo5WxC0upk/bAdrc2NN2voRE\nOyk8KXnRzgUw7xcXtfCYvn3h59gJkDOeUNkWqGEtILZ+HYSuEn7lhtFntn7t6o+XRFlc9KcP+DTM\nlE5Jev7Nd6H6ZfjIhvG+XQ2jftUKK4w3hTzbtcTsJ2nkZm+6aYuVbyuDu4YRqLN1nYYjuOW7dTKn\nb2NfM7oyD+i6daNptX37deKlCX8bcTAjtcVFvdLu9a+fPMc+z3xuX5c0IoxdWbfDlXq7FnkKgfqV\nQdP0y5y7FuP65du6xWy3Ys77/RtGi0/cjclNbqWtXebY1jAzvedqjq0L5pxQ/UqLylO/gmlM5IlE\nwiztrxojPKavBee9zI3ExR7ZmWOzlx0wSuxUJ0ftnF0+vG+0GsZeiguMEh1/8d369ePf0KID6Nc7\nvjVKvFy1Sq+2c0sYbN8+Eic3emSP6OyVeknFNGMvweWSXjJr1K1fGPaXpF+2QbL1ayf0vnV2FGvh\n9Kidx46MasyZZ9gtJWCv6jUJ2+/99LiGvX24XdTKlcCWoQ66uuAml/vOcVca+6B+TQWap6rJEhR7\n9UdeASojVkZYDlrvl+8GtvxycqVwc3wZtNiY/s1ocN2F1vLh9cBTQ1Nj6qC4K+qMKB27Fdj8gD7e\n/ACw+JHx/q+9Vo/q7KKYBtcQFdlxPPYSXC7pJV2hqfpl+nP1a8UT+jhp+xWjX+ZzY7QOArjkleMD\nLaM3tn6ZKJNtZFat0kWANx7R73eo8cKZL3zhZMFfm7TyK9SvxsJpu6oJSWhMSoaM0XZWX2kbambt\nLg6MNs1838PAb1wMzD8BLJnTy4ePb9EhcDex0Q2HmxUqJqE7iRtvHE3VGZLCynPna+GxQ+wh2w7E\n3oyzg8LDabuMc6lf1emXr780svTNtLP5AWDPFaNdEOwpOl9ittEwc+6BA+OLUgDg298GnvWsUX92\nkrfRPrsvG18OFPUrCrH0i+apKvIKijs/H7NtH0a4zIq4eRkJx7yM2rRzn8xqFF9RTJu7LgQWngZ2\nXpu8m7dramwxMtx8M3DNNZPCZYRl6VJtkAxuW26uAokCzVPGudSvuG0nYTTJp18mtwnw65ebu2kz\nmAdes2z03qdfwHi0CBjladrn+vIybe3L0jCA+hUZmqe2kGd0FZp06Uvmzovd1/LdwNs+Ctz2duDx\ntwKvOzn6zjVKnxlOx9mCZJcreN/DwOnP6+k7X8JjmgjYozw70dI3AkpLCvflEtiiR0pB85RxLvUr\n7Lwy+mW3s3w3cPUx4OmjwL5PheuX/b3RL1NHCRivTxeaPG0MVJZ+Afk0jPoVjVj6FaNIJkkjreib\nKzahoy8z1x9aGNOH3dfxLcClAI4fHn3nC42bkd0jFwAYJomrk1qQAOAzm4F/uHxU28kUb7OLw6WF\ngs3y22c8Q4vHnkeBzc/x1zyxExrvvnsyV2HlyvHIE4WHkPxkaYytYXXql93f8S3A7uEgbi/C9WtX\nH1hxaGRU9r5jsnbTqlWj6t92xClJT+7YOzJQG9YDV/ycv3Clq2Gf+YxfwwDqV0NhwnjVpAlK6Ly9\nwU7OzJo6S2sjrU2zm7lbV8UYpKf2jwrJmd9t734tSHv3jG+bYD/0WfvOGTZtBC6/fHj80VFofO78\n0TluQuM114y+O3pUC9ji4iifiiFvQoqRZYjyaFgM/TLtuO9vOzlq1xi6NP3a1dfvd1ylB2g33ghs\n+oHkrVSMjmVp2LKl4yuF+4dHnyctbDlxAnj3u8cT08+c0cfUr8ZC81QHWYYldNWJKwh58gTc5btJ\nmORKt5TB3v1aeM5bP7mLOTBemff4Fn8Byqzlr0Z4VvdHhTe/cx3wqaf0irtffPfoOt9+UIB+bwtY\nUnmConDZLplFYmhYGf0CRgO7EMOWpl+PHZk8343sJBXRDV3Cb2uY6uvq4KYEgt3HzTePbyx8883j\nmhZbv8zvQUpD81QH7sPuE5E8y3aLhLvt5bu22Nn3Aut7u5TBGgx3MR+OtuxdzF3M564YhSx/Necc\n7o9KF+xQwBdu0scXXjga9dl7QpkRm3mtaqltSOSsKBQ00mSyNCwPRafrLoPfsLlTgLbWnS1lsHSk\nX2nGx40O2eTVsG2P68/s0gUrV2oDZfo5+eS4fp3O2ES4LFVp2AzqFxPGqyRrVYmdL2AnT4YmXha5\nB9992H264XSzKmUrRhXDTRkCF/v7//BV//2ELrkFxvdy8pUxsJcIu5/blYCTmD8UvpFm0qrBGLSs\nmi8TxjPO7Yp+AeEaVpV+ATri5JYbSNoVIUnrdgJ4z+7xhSwueRa2ZD37c+fraLmtW66G1alf5r6r\n0LAZ1S9GnqokK0xtojruaCpvLlToPQDJy4IPWt9vtY7taNWVW3Ruk2ucZG4YmTIFMjeP5yjZhDys\nJqL00U2jz3o7Jotk2kt53c9DBM7kI4RSxWiQ1XxJkwmZaqtSvwDgg849JEWuTAQqaQ+8m7YAS16R\nXBfOreydRMizf/LJ0X6aht6O8XPq1C/TXmwNm2H9onmqg7QwtS1OrllJm8rLW513Aen1VOxKu2us\nc13htHObgFFV3+/92ChR/OmH9ajLDQ3nfbDOnBnPHdhzhRY/VwBMfoJPGOYPTfZ//Wd1HgKgX+1z\nfPe7Yb2OggGTeRBlYTVf0gbSDEsd+gWMtoVK0rBtGNcv9/6kD/z6j09et2G91it7+iyGft2xV7dl\na9i8AJe/Qx8vLmbrFzDSp7z6Za4x+rVhvT+XqwwzrF+ctmsKSVN4SZStkeJrDxhVDLfbNvdmNvw1\nyJzezsWMrh65QH/23EdH59x8sx6FZYV23ZGWG2KeFy1+gH/T4KR2VF9f5yucufdFozZd3JorpgZM\nVfVWWlTNl9N2GedSv+rXL7tNt+1dfW0w3v9qf9kTOz3gQx/SK98Mofpl2grVsM2b9Kpitz2fDqi+\n3vrFp1/n3pCsG3ahTerXWTht1zXspPGkUZ6JDhVZqWeu931m5zklJWSe3fDXWmmnTo5C0U8NyxRc\neOvoezOCszfZ9IV2fUmM7oimv1ofb9k8vhTYtOVW+50/ND46e+c+fWxv8jmY1yv4XOxQ9NvvGB81\nvvfTk+fHoCXCQ4iXaekXoCNObgkE07YpSfDAh5MTpY2GHT0KPGklcJ86FaZfQH4NM+VYXA1zI05G\nw/a+SEewXP1K+p2MhrlRL+pXNBh5ahJZozH7+yIjt7Rr0iJPZgsXwyMX6Cq+biLpfftGK1pcTATH\nN9JKS2K0BcWM1mzs/aJ8o8KkyJONGV3a2CNNQIvXxiOtSIisGkaeMs6lfmV/H1u/7O/NqzFOSRqz\n51FdS85wx9v1JuU+kvQrrX37e1Mg2N2b05yfpWFu5MnFF1GyNe+D30v9GsIK413CTrI0+y65K1rc\n7/Ms981qH8PPXmYd25jCl2YlnToJ3DUHfPmJUTj4m88BnusYJ1dwfKHddetGx745c3vfKNc42Vsn\nmFe7EjAwGu25Vc7f+U4tQidOTBon+3xAC+Ng3t++TYtC14REown6BeicqG3D741xAkYRIGNOjAYc\nuxXYcavWsAs+oKNMNqdOAcuWZa9+C9Ew396d5nwgW8Nc/TJmK6lPQF9z993696J+RYeRpyaRZ+Q2\njfZNzpN5NXtKfepyvQJvuRN5MlsdJG3LkrVpsHueWep77Bjwx388XvW3yFLZufP9xsklpP2WLdct\nAyNPGedSv4p9X7Z9c05SHmPSZuSXv0NXGU/anNxoRFENc/UL0APLO+8sp2FuqkIS1K8xuDFwF8mq\nj1K2fop7vXnv1lF5EHp5sA+31tOEoRoaKPshTHswkzbH9F0/L7p4ptlw06bqUVNa+1XWgGogNE8Z\n51K/in2ft337vathg0uAq94y2YarL67JAUZTa7E0zNavjUeSoz9Vahj16yxMGO8iWcISIjxpCZhJ\nofQPQheQw/A1yTi5tZzsCJRJJj++QedEGYHJqgNiL531nbtsqU72tssV+Jbmpj3soUuM085La3+G\nl+sScpY69QvQ+mXO/yC0McHw1WecfPri7lm3bp0enNnL+cto2Nz54/q190XJSdtlNYz6VSs0T13D\nFhSXNdaPnUOwdqk2Q19bnb71ir3pr8l9cg3Vkh8Adg5zB+YPJT+Yvm0Qkjbk/Oim0Y7oOxTwW68L\n/GMgfDuCstsWxK6fQsgsEqpfO53zf/8GbVa+tjp5RVyaFvlMlRmk+a5z20/SsJNPltMvQOd7ZmkT\n9at2Sk3bicgGAH0AlwP410qpL6ecy7B3ldiGCPAnVSbtZL4A4L6MrVds3HpPZirv6YeBcy4ehavN\nSjdgcpSXp+YToEsKrLsFuGf7eBJn0vlG3ELC0TMWti5Ll6btQjWM+lUxefRrAZNVzPurdTmCotur\n2PlOPv2yrytS8ymvfpnP3dV5Pm2ifuWiKdN2RwD8JICc+1yQ6JitCQz2CM6t5eS79vhw65XH3zr+\nnS8K5VYZP75Fr7Y752L9/p37Jivgpo3yXNyR3bKlemXc4b4WHvta34jLfPb616f80k5/ZcPWM7Qt\nQceghjWBPPq11jlX+nqPt337dYK3r5SAje/5NitrAb9+metCtyOxNSyvfgGjz23jZCqS+/qiftVO\nqVIFSqkHAUBEONRqAmakZjby9X1uOP5sYPk/6mOz/Pe+j41Hn7I2ArY58y3gyg261tNHNgBXyeTI\nDfAvGwbSI1OnTo/qm9jLh4HJJb72ZytXjpdLSBMVeylwXtKSUEmjoYY1iDz6dcXlAI7pY9VPjjyF\nriIzuvTOCPrl69dcY0grUWAbNMPion+hjIH6VTus89Q1fEmZ9qqUtQDeNges+Af9fl6A9+wDHn8T\nsMLKXTrxC+O5TI+/dTLi5PLK/cC+C4B9w5wnU2MJGH8gk2qWmF3GfYbIFIZbtmx8+bBPyOzPfGKQ\nJA5FR2z2/S5ZklxMjxCSTlpSuZmue9scsGK4SbhZgXvnncD2lIFUVo0jQD+ve3do4wQU0y9f0viB\nA1oXbMz1Pv1yP3fvm/rVCDLNk4jcA+A59kcAFIDrlFJ5Sp2h3++fPe71euj1enkubydll+cW6ce3\nKuWg9bk6CTy1X5cUGMzr18ffNF4I88y3xt9n9W1GhbedBJ57SIfRd1ylP/ON/nyF41at0kXdsjht\niYRvxOV+ljYiLIstdCbKZX6XoiPBFjMYDDAYDKZ9G2PE0rCZ1C+gfg3z9WVHo3ZeC2waPmuDef16\n553ZA6k05g8B/eHMrYlkFdEve/rPsHTp+PTbsWN+I2aT9Dn1q1Ly6FeUOk8icgjALzFh3EMVG2CG\n9pOVhGlqMtnTcr7Nf7/3Y2HTd+62CEB6ETn7cyCs1pN9Xl6qTKzMk0zacrqUMG7I0rCZ1S+gWRpm\nmxrfBuFJi0fStj+x8SWJ59UvILnWk72xb159oH5FoVFFMofC88tKqS+lnDNb4hOyeqQsMuffY863\nSiVJ/Fyz5Pve3dcu6XwzcrTNE5D8QNqfZ4Wm8wpHUmi7DnHoeM5Ah81ToobNnH4B1WuY0Z7QlXa+\nyuEhz1oe7Zg/NDJnhiL65bu3ufPHt4BJuw/qV2U0wjyJyOsB/A6AZwN4HMD9SqlrEs6dPfEBwkdt\neUPjbjJ3Wj++trNMU1pfWbjmCUhfjhv6uSscPqHznRfaJwmiS+YpVMNmVr+AMA0rq19Z/axYPfms\n53mOy5qOGPrl3scVP0f9mgKNME95mFnxCRWVPKFxXzTodSfDxSuvGTJ9hpotn3nKQ5p42MLhWw3T\nlZonDRbILpmnUGZWv4AwDSurXyYC5evHpydF94Kr45kKNT/Ur6nQlDpPxJBUFTdLdOwaJjtT2rFR\nJ3VZAGBU6TvUOPm2WEmqKO72GYO589O/T9qixWAq//rqsJjv277VQNlqwYTkJU130rQlln4l9ZMU\nxXa3QAkhhhaU1S9Ab89C/Wo9jDzFoo4dw93z35YjGmTzAw/pgpZPPwyc/vNRFCqkHAGQPhJNizr9\n4rt1suSJE8AtH0o+zx65AcmjON/IzdDgkU8qLRh5MvKUcS71K+z8LP1K0hJ78YidgJ1ndVjSlH8a\n1K9sZki/GHkqS5GRl4/QBdN2f7edzN+fzI0qgZ9z8XgUymzum4W7NUJI5Gru/FGtpgsvTB/BmX2W\n3JpP7gjOrsPi0rAHNpgiI09WByZFaZt+uQUkjaasWpUv2mHKEtjtpkH9CmOG9IvmqSwHMRpxbUXx\nncNDp92K9Gfjbu7r1nBK2xjYJ7TLd/tNl/tAnHxSj9gA/WqKXCZhNgVOexDzjhyTaNrDm2eTzhkJ\nkZOKKKonrobF1q+kqJOrCW79prQtUwD/lH/SM2S302T9cu912syIfnHaLhZ5VpvEqJtStnCdnQCe\np5YToO9/hxoV2jQ8coEuYJeWMDl3/rjwhC41LjoSy6ou3ua6JlMKkXPaLuPcrusXUF7Dik79G+zn\nOE8dJ0A/83tfBGw8oqNDvmcoSRdc/XLvJetei9BVDWu5fjHyFIvQiFOMEHlof2nYuQbq5Ghj4Czj\nJHO6qi8wLLA5fFhN4mfWxpm28ISOOtyNgkNx23ffh27y2VS6kFxKmkGeiFMMDStjnIDx/9dPnQ6P\ndphnfjA/XnTXvGZt/usapxANs01eXrqsYS3XL5qnOik75ZaXpOm3pM9DksXVSeBTl+vjp/YAxzeM\nm67QByLPQ18ktOtbkeP21/KHF0C+EDkhZalSw1zjlKYJvu9Cnl/zzB/uj5559xmihtVHi/WLGwNP\ng1w7AhYkqZZTkRpPLse3jK/MM69G/EJ2+E7aFNPFt8lmHpE07Z980t9fmd3Im0Kb7520k9gaFror\nQdZ3Ifie+dC95WyoYXFo6X0z56mLJBWhy7PVShGKFMd0cxdib0mQlS8Qcg0ZgzlPGedSv/Lhizgl\n5cI0cSl80zSM+pVKLP1i5KkLuNW/zYo6E2GyI0S+z+sg6YE2n6WJS5mRlXtNVhttTb4kpK242pAW\n0QmN9tRxn/Y9Ac3QMOpXbTDnqe0klQpISgAPTQwPIaS+E5A93x+SO1BHvZA2J18S0kYeO+LXhrRc\nmFh5Mnme7zo0jPrVKmie2oxvqxWbpMhSjIjT8n1hRTVDRcVOeixCjHohbU++JKRN/P4N6dqQlW9U\nhk0bw/Uir4YtLua/P+pX66B5ajNuwcu6puGW7xvVd0orqgmEP9AHDujRJJBPRMzeUVn7SYXS4tUf\nhLSGXf3p/WO/aSOwcqU+DonQhNznsqVaMxYXddt5TJBvFR31q/EwYbwLuDlP7vu812edayedP7Vf\nlysA0hPG05IYzTy9ER5DVjJo0h5S9rErIkymLA0TxjPOpX6l40sQz7uoI+36rHPthPPFReCOveX6\nKapf9rX2Hn0A9atCWCSTjLCNT1IOVBJ5z3d5/E36NWulXcgy3pUrtQAB2SNR3/LfrP2kQrZhIIRU\nh08n3GrZeSLPZae77rwz/NykiFMR/XKvvfBC4OabqV8tguapS2TlQJU9H0ieKiz6ALsh8Tv2hoWd\nfaH0tP2kkvIWioovBYuQfGQNsPImPBdJkE7ThyIU1S/ftSefpH61CE7bTZuye9S55C2CWbRopj3V\nZ9ooszw2Kxw9f8i/kWbSdb7P3WW8RWvGzPhyYE7bZZxL/ZoktAZc3mer6LNo60OM57mofiVdS/2q\njFj6RfM0bWJsEgxMbvRbVc6T71o7B8p9gEMKyIWg+oD0i11rU3ZTzSYW6asZmqeMc6lf44QaJ/Ns\nVpnz5Ls27Xkum49loH41BhbJbDtrAKwdHu+E3u6gaATKjR7lNUJVrdJLSujOM9qZPwT0D+tj1Qf6\nq5NHcCEU2YbBvX5aRfoIaQqh+lVVxMmmqmfQvSezSi/PPVK/OgsjT9OmbOSp6i1XQvBN27kjHJsi\no51YI7dYzPCqF0aeMs6lfmnyRJymGQ3xGTf3nr76VeCFLyx+j9SvxtCI1XYi8gEROSYi94vIfhF5\nVpn2ZpKyG2xOq9aTzfEtk0mSbtG4svRXl28jJjMqPF2DGlaSJP3Ks8/ltIs7+mojufplG6ciRTCp\nX52jVORJRF4J4F6l1BkRuQmAUkq9N+FcjtyqpEzeUgySxNKMcGY8SbFLdCnyFKph1K8cFNkgHGhm\nNMTVrzx1oUgjaUTkSSn1GaXUmeHbLwC4qEx7pATTNE5pGDF0R3fzh8Lb4LJaUhHUsMgUNU5A84wT\nMKlfxjjl0S+AGtZBYtZ5+hkAd0dsj3QNWxxNEmUWMfZ8IiQMalgZyhinNlBEvwBqWEfJXG0nIvcA\neI79EQAF4Dql1MLwnOsAnFZKfSKtrX6/f/a41+uh1+vlv2PSbvKsPvFVEY89Om3iVEGLGQwGGAwG\n076NMWJpGPUrha4bJ0Pe1XNVaxj1Kyp59Kv0ajsReTOA/wfAy5VS30k5jzkDbSYrpypkXzv7QQ9d\nfVJlrhTzsArRpZwnIEzDqF8ptMU4FTUavlp1eVbPVaUz1K9CNKLOk4hcDeBaAD+eZpxIy8mqQp4m\nnr6NL/ftD199kreOiU2aWPpGhABHcTMGNawkbTFORY1GUq26fo6cpyo0LCmixUhUbZTNefodABcA\nuEdEviwit0a4J9Ikiux/Z3A3vgRG+zLlKRRXRAyy8gzc5dHr1jEvYTahhhWlLcapyB54vuvs4/e/\nOt89xNYwX3kH5lbVStnVdv9cKXWJUupHhz/viHVjpCGUqSNlP+AnTujXOuq4hIqlWUGTtpM56TTU\nsIK0xTgBxetIudfVWYsqRMPsFcxFDSIpDLdnIdkc3wI8/la/ccoSUTtkXVdIOc8WBOY7bllASBht\nMk6GolNn7nVVLFrxEaph5nNuu1I73J6FlKPJQlrnBqMzRNcSxkOgflk0+ZnvGtSw6DSiSCYhjSav\niFB0CEmHxqleqGGNheaJ5EsCT4Jz7IR0m6YapxjaQ/0iOWHO06yTVYYgBNYbIaTbNNU4xdAe6hcp\nACNPs0yZMgQGrvIgpNs01TjF0B7qFykIzdMsU6YMgeHUab3TOMBVHoR0jaYaJ6B4CQK3DeoXKQBX\n25HsrVfSeOyIHrEtLo52HI8NV5A0Cq62yzi3K/rVZONkU0YfzJRdlfoFUMMaRCO2ZyEdoahxkrlR\nyHvlymoEgvkIhNRPW4wTUFxz7Cm7qvQLoIZ1FE7bkeKok9VW3WU+AiH10ybjVIYY035ZUMM6CyNP\npHJsUZEAACAASURBVBxlNr3MglVzCamXWTFOhir1C6CGdRiaJ1KeKgWhanEjhGhmzTgZqtYWalgn\n4bRd24hR0LJtUHQIqZY6jdMsTl1RwzoHzVObWL4bWPGEfm0CszpSJaRL1Pkcb1gPbN+uXwlpMTRP\nbSFGQUtCCLGpO+LE5GnSEWie2kKMgpaEEGKoO3Jcx+o2QmqCCeNt4vgW4PG30jgRQsoxrSl3Jk+T\njsDIU9ugcSKElGHauYo0TqQD0DwRQsisMG3jREhHKGWeROR6EfkrEfmKiHxKRJ4b68YIIaRqZkrD\naJwIiUbZyNMHlFL/Uin1YgB3AQjfcY+0Gwox6QazoWF8XgmJSinzpJR6wno7B+BMudshreFn+9O+\nA0JKMxMaRuNESHRKr7YTkd8E8EYAjwO4qvQdEUJIjVDDCCF5yYw8icg9IvLX1s+R4etaAFBK/ZpS\n6mIAfwTg56u+YUIIycNMaxijToRUQmbkSSn1qsC2PgHgkwD6SSf0+6Over0eer1eYNOkcVCUiYfB\nYIDBYDDt2xgjloa1Tr/4jBKSizz6JUqpwh2JyA8rpf7b8PjnAfwbpdSmhHNVaF+yVQrfE6kAmZus\nL0VhnlmUCs+pFhEopRr7QIdqWOv0i8/nOMuWsr4UARBPv8rmPN0kIi+ATrJ8CMDbSrZHmsby3Xov\nvaf26ArnhHSL7mkYjdM4G9brvfSOHtUVzgmJQCnzpJTaEOtGSANxNyPm1jCkY3ROw2icxnE3I+bW\nMCQSrDBOkuFmxIS0BxqnSbgZMakIbgxM0uFmxIQ0HxqnZLgZMakARp5INjROhDQXGqdsaJxIZGie\nCCGkrdA4ETIVaJ5IfijYhEwfPoeETA2aJ0IIaRs0ToRMFZonQghpEzROhEwdmidCCGkLNE6ENAKa\nJ5IPijchhJAZh+aJEELaAAcuhDQGmidCCGk6NE6ENAqaJ0IIaTI0ToQ0DponQghpKjROhDQSmidC\nCGkiNE6ENBaaJ0IIaRo0ToQ0GponQgghhJAc0DwRQkiTYNSJkMZD80QIIU2BxomQVkDzRAghTYDG\niZDWQPNEwqG4E1INfLYIaRVRzJOI/JKInBGR74vRHiGE1MlUNYzGiZDWUdo8ichFAF4F4KHyt0MI\nIfUyVQ2jcSKklcSIPN0C4NoI7RBCyDSYjobROBHSWkqZJxFZB+DrSqkjke6HEEJqY2oaRuNESKs5\nJ+sEEbkHwHPsjwAoAL8GYDt0uNv+jnQRij1pKY3TMD5LhLSeTPOklHqV73MRWQXgUgB/JSIC4CIA\nXxKRlyilvuW7pt/vnz3u9Xro9Xr+PneqrNsidbNz2jdA2sBgMMBgMJj2bYwRS8Oi6RefJUIaSR79\nEqXiGBUR+R8AflQpdTzhexWrL0JIOxARKKVaEZFO0zDqFyGzR5p+xazzpMBpO0JIe6GGEUKCiBZ5\nyuyIIzdCZo42RZ7SoH4RMnvUFXkihBBCCOk8NE+EEEIIITmgeSKEEEIIyUEjzdM0ljqzT/bZpv5m\nqc82Mgv/bWbhd2Sf3eozZn80T+yTfbawv1nqs43Mwn+bWfgd2We3+uy8eSKEEEIIaSo0T4QQQggh\nOai1zlMtHRFCGkVX6jxN+x4IIfWTpF+1mSdCCCGEkC7AaTtCCCGEkBzQPBFCCCGE5IDmiRBCCCEk\nB403TyLySyJyRkS+r4a+rheRvxKRr4jIp0TkuTX0+QEROSYi94vIfhF5Vg19bhCRoyLyXRH50Qr7\nuVpEvioifyMiv1pVP06fHxORR0Xkr2vq7yIRuVdEHhCRIyLyCzX0ea6I/MXw/9MjIjJfdZ/DfpeI\nyJdF5EAd/XWFujSM+lVJX7VqGPWr8r6jaVijzZOIXATgVQAeqqnLDyil/qVS6sUA7gJQx3/UTwO4\nQin1IwD+FsB7a+jzCICfBHC4qg5EZAmADwN4DYArAPyUiLywqv4sbh/2WRdPA9imlLoCwMsAvLPq\n31Mp9R0AVw3/P/0RANeIyEuq7HPIuwAs1tBPZ6hZw6hfEZmShlG/qiWahjXaPAG4BcC1dXWmlHrC\nejsH4EwNfX5GKWX6+QKAi2ro80Gl1N8CqHIJ+UsA/K1S6iGl1GkAuwH8RIX9AQCUUn8G4HjV/Vj9\nfVMpdf/w+AkAxwA8v4Z+nxwengvgHACVLpsdmoDXAvjdKvvpILVpGPUrOrVrGPWrOmJrWGPNk4is\nA/B1pdSRmvv9TRF5GMC/BfDrdfYN4GcA3F1zn1XxfABft95/AzU8lNNERC6FHkn9RQ19LRGRrwD4\nJoB7lFJfrLhLYwJY2ySQaWgY9SsqM6VhHdcvILKGnROjkaKIyD0AnmN/BP2L/RqA7dDhbvu7Kvu8\nTim1oJT6NQC/Npzf/nkA/ar7HJ5zHYDTSqlPlO0vtE8SDxG5AMA+AO9yIgCVMBztv3iYY3KniKxU\nSlUypSYirwPwqFLqfhHpofoRf2uoW8OoX9SvKuiyfgHVaNhUzZNS6lW+z0VkFYBLAfyViAh0KPhL\nIvISpdS3qujTwycAfBIRxCerTxF5M3Q48eVl+wrtswb+HsDF1vuLhp91DhE5B1p4Pq6U+pM6+1ZK\nfVtEDgG4GtXlI10JYJ2IvBbAeQCeKSJ/qJR6Y0X9tYa6NYz6VSszoWEzoF9ABRrWyGk7pdRRpdRz\nlVI/pJT6P6DDpS8ua5yyEJEftt6+Hnr+t1JE5GroUOK6YSJd3VQVRfgigB8WkUtEZBmALQDqWqUl\nqDc68nsAFpVSv11HZyLybBG5cHh8HnR046tV9aeU2q6Uulgp9UPQ/x3vpXFKZxoaRv2KzrQ0jPoV\nmSo0rJHmyYNCPf8z3SQify0i9wN4JXRmftX8DoALANwzXEJ5a9UdisjrReTrAF4K4KCIRM9TUEp9\nF8DPQa/GeQDAbqVUHWL+CQCfB/ACEXlYRN5ScX9XAvhpAC8fLr398vAflCpZAeDQ8P/TvwDwp0qp\nT1bcJylHHRpG/YrINDSM+tUeuLcdIYQQQkgO2hJ5IoQQQghpBDRPhBBCCCE5oHkihBBCCMkBzRMh\nhBBCSA5ongghhBBCckDzRAghhBCSA5onQgghhJAc0DwRQgghhOSA5okQQgghJAc0T4QQQgghOaB5\nIoQQQgjJAc0TIYQQQkgOaJ4IIYQQQnJA80QIIYQQkgOaJ0IIIYSQHNA8EUIIIYTkgOaJEEIIISQH\nNE+EEEIIITmgeSKEEEIIyQHNEyGEEEJIDmieCCGEEEJyQPNECCGEEJIDmidCCCGEkBzQPJExROQF\nIvIVEfm2iDwtItcVbOdNIvK5Atd9UkTeUKRPQgghpA7OmfYNkMbxKwDuVUq9OEJbyhyIyBkA3wLw\nPKXUmeFn5wD4nwC+Xyn1DABQSr02Qr+EEEJIZTDyRFwuAfBARW0fB3CN9f4aAI9V1BchhBBSCTRP\n5Cwi8l8AXAXgw8Npuz8SkeuH360Wka+LyDYReVRE/l5E3mxd+30ickBETojIFwD8n54uPg7gTdb7\nNwL4A+ceDonIzwyP3yQinxORm0XkMRH57yJyddzfmhBCCMkHzRM5i1LqFQA+B+CdSqlnATjlnPJc\nAM8E8DwA/x7AR0TkwuF3twJ4EsBzALwVwM+4zQO4E8CPi8izROR7AfwYgD/JuK2XADgG4PsB3Azg\nYwV+NUIIISQaNE/EhyR8fgrAbyilvquUuhvAEwAuE5ElAP5vAO9TSv1vpdQDcCJKQ/43gAMAtgDY\nPDz+Tsa9PKSU+j2llBq2+VwR+YH8vxIhhBASB5onkod/MsneQ54EcAGAfwbgGQC+YX33kHOtMWQf\nh56uewOAPwzo85vmQCn11LCdC/LdNiGEEBIPmicSg38A8F0AP2h9drHvRKXU5wCsAPADSqn7arg3\nQgghJCo0T6Q0w2jUfgB9ETlPRFZiPDHcZQ2An7DeJ00TEkIIIY2D5om4qOxTvOf+PHQy+SMAfm/4\n4z1XKXVMKXUsoZ2s/vPcHyGEEBId0Xm4hBBCCCEkBEaeCCGEEEJyQPNECCGEEJIDmidCCCGEkBzU\ntjGwiDC5ipAZRCnF1ZSEkE5Ra+RJKRX0Mz8/H3xurB/2yT7b1F9b+iSEkC7CaTtCCCGEkBzQPBFC\nCCGE5KCR5qnX67FP9tmqPmfhd5xWn4QQ0jRqK5IpIoo5EITMFiICxYRxQkjHaGTkiRBCCCGkqUQz\nTyKyRES+LCIHYrVJCCGEENI0Ykae3gVgMWJ7hBBCCCGNI4p5EpGLALwWwO/GaI8QQgghpKnEijzd\nAuBaAMwIJ4QQQkinKb09i4i8DsCjSqn7RaQHoPTKGtnKxTmETIVd/cxTlJqv/j4IIaTBxNjb7koA\n60TktQDOA/BMEflDpdQb3RP7/f7Z416vx5oxhDSJAOOUxWAwwGAwKN0OIYQ0mah1nkRkNYBfUkqt\n83wXXOeJkSdCaiaHccoTeWKdJ0JIF2GdJ0JmnQgRJ0IImSViTNudRSl1GMDhmG0SQiqExokQQnLD\nyBMhswqNEyGEFILmiZBZhMaJEEIKQ/NEyKxB40QIIaWgeSJklqBxIoSQ0tA8ETIr0DgRQkgUaJ4I\nmQVonAghJBo0T4R0HRonQgiJCs0TIV2GxokQQqJD80RIV6FxIoSQSqB5IqSL0DgRQkhl0DwR0jVo\nnAghpFJongghhBBCckDzREiXYNSJEEIqh+aJkK5A40QIIbVA80RIF6BxIoSQ2qB5IqTt0DgRQkit\n0DwR0mZonAghpHZonghpKzROhBAyFWieCGkjNE6EEDI1aJ4IIYQQQnJwTtkGRORcAJ8FsGzY3j6l\n1I6y7RJCEmDUiRBCpkpp86SU+o6IXKWUelJEngHgPhG5Wyn1lxHujxBiQ+NECCFTJ8q0nVLqyeHh\nudCGTMVolxBiQeNECCGNIIp5EpElIvIVAN8EcI9S6osx2iURWDPtGyBRoHEihJDGECvydEYp9WIA\nFwH4v0RkZYx2SQTWTvsGSGlonAghpFGUznmyUUp9W0QOAbgawKL7fb/fP3vc6/XQ6/Vidk9s1mBk\nnHYCWABwcHq3QwrSMuM0GAwwGAymfRuEEFIpolS59CQReTaA00qpEyJyHoA/BXCTUuqTznkqtC/Z\nKqXuiVjsBLB12jdBCtFQ46TUfPC5IgKlFB9oQkiniBF5WgHgD0RkCfQ04B7XOJEpsjDtGyCFaKhx\nIoQQEqdUwREAPxrhXkgVcKqufdA4EUJIo2GFcUJCqGvVIo0TIYQ0HponQkKoY9UijRMhhLQCmifS\nXWJEi9ZAJ91j+FpVBIrGiRBCWgPNE+kuMaJFBzFarbgV1eSQ0TgRQkiroHki3aOKaBFXLRJCCBkS\ntUgmIY3g4PAnZo2rKiJOawCs61fQMCGEkCph5Il0l6ZHi7h1DiGEtBKaJ9Jdmlrjyp5WVH1g/tAU\nb4YQQkheOG1HpscaNNfgVImZqlN9QPpTvBFCCCFFoHki0zMxa6fU77RwV9X1V0/lNgghhJSD03ak\n/tybumonNQlfOYIdV9V+G4QQQspD8zQr+AzKtExMHbWTCCGEkIqgeZoVfNGlaZuYpq+GiwWLYBJC\nSKegeeo6IdGlUBMTOzI1CxEnGidCCOkcNE9dJyS6FGpiWJcoHzROhBDSSWieZoUyU2SzmOBdFhon\nQgjpLDRPs0KZKbJp50a1DRonQgjpNDRPJJxZSfAuA40TIYR0HponEk7TIk5Nmz6kcSKEkJmA5om0\nlyYlsNM4EULIzFDaPInIRSJyr4g8ICJHROQXYtwYIYkkJbDL3HTuh8aJEEJmihiRp6cBbFNKXQHg\nZQDeKSIvjNAuKUPRKS33uqZNjQH+BPblu4EVT+jXOqFxIoSQmaO0eVJKfVMpdf/w+AkAxwA8v2y7\npCRFp7Tc65o0NeZiEthlDjhvsz4+b3N9ESgaJ0IImUmi5jyJyKUAfgTAX8RslwwJiQIVrcnkXret\nYDt1YhLY1UngqT36+Kk9+n3V0DgRQsjMck6shkTkAgD7ALxrGIGaoN/vnz3u9Xro9Xqxup8N1kIb\nhjVIXvl2cPizE6OprRCSrsvbzrQ4vgV4/K00TlNmMBhgMBhM+zYIIaRSRClVvhGRc6D/6b1bKfXb\nCeeo0L5kq5S+p06xBpPTZ1mGJs1g5bmuaDtlmVa/Icy4eVJqPvhcEYFSig80IaRTxJq2+z0Ai0nG\niZTETpA2ZE2lFTUe7nXTMjBNzbWaceNECCEkTqmCKwH8NICXi8hXROTLInJ1+VsjEyyg+9ukNHkf\nPRonQgghiJDzpJS6D8AzItwLycKYpSq2SZG5evKFsiias1U1NE6EEEKGsMJ4G4kdcZpWjaQ0kgxi\n3khUjMgVjRMhhBALmqdp0ZTpqKwaSU25T0PeXKiyuVM0ToQQQhxonqZFUxKis2okTes+3X6zcqF8\n78vmTtE4EUII8UDzVDdNTIg+vgV45AL9apjWfSb1a684XMDk1KVrtnxbuBBCCCERiFLnKagj1nka\nx5cQbRuFpjCtxO2kfk3Nq63Oe8ODAD7onF/k78moUyKs80QImXUYeZoWvoTotWjOdJ6hipV9Rfu1\njZKJSrkRpsuca2icCCGERIaRpyoJjXr4Koj7pqa6TujfyxeV2oZx41T070fjlAkjT4SQWYeRpyoJ\njSK5FcRnNUcn9O/li0p9EOk5USHQOBFCCAmA5qkKiiZbL2B602R1kPR3yPv3SjJGa6D/fkWmPmmc\nCCGEBELzVAVFV3qZ6tp5iLEKrq6VdGmmpmzUyLTv5kSFMAvGadnSad8BIYR0BpqnKqkjihQjwTyk\njVAj4jsvK7Jk+i8aNbLbN4SY1l392TBOG9YD27frV0IIIaWheaqSKvOWQqe60kxPnumyUFPjO88X\niVvj6b9I1MjXfpenPvOybCmwapU+XrWKEShCCIkAzVNbCZ0aTDM9IW2sXRpu0rLOs03NWk//ZYta\nmvbda333MgsRJwA4dRo4elQfHz2q3xNCCCkFSxW0naTl/W75AzeXyL4uqY3lu/V+d0/tAW7aElYs\nM6uopu++EHAvCPw+5J5mxTjZLFsazTixVAEhZNZh5KntJBmJgxgZEzuSY6Iwa51zXdwNgw8GTvdk\nTZnZfflMVpYxypMT5YuG5TFOXZriYsSJEEKiQfPUZdZCb1fifhYyDeduGLwQ+I9vVtRop/M+tKp6\n2rRg0u/gTguu6wd0NIRJ1oQQQhKgeeoittG4DJPJ2Yas3KLjW4B/tWx8w+AyuJv72qYpJEk8KScq\ny3wtIH/EiUnWhBBCEqB56iL2lB0wMhdFVqStiTjdYxu4tZ57SDNydjkDX3tp5stEnEJNEJOsCSGE\npMCE8a4ic8DrTmrTYecWZSV0G7ISzstgDM97dgNv2ZK+J13IfaT9TibitGG9jiIdPQrs2x92nxGT\nrLsEE8YJIbMOI09dZPluYMUTwH27R5EaN0qzLaONrDIGMpd8bWq9KCv6c9MW4G+Xpvfj3sddnn6T\nomjGOOWdhjPf0zgRQgjxEMU8icjHRORREfnrGO2RErir5IzZcE3IZZ5rffiMiTFny3f7r0nLQVo4\nrSNOgH5dOD3ak87GNWcLKf26kSpgPMcpzzQcE8UJIYRkECvydDuA10Rqi5TBXSWnTupjYyoeRL7K\n5L6Ik23ObJMTuiLu+BZd+sAkopuCmQZjkjZaBuaulH5t1sKfHL5vP3DjjelTdqERKiaQE0LITBPF\nPCml/gzA8RhtkRRCtyw5vgV45ILxVXImGvRBlKtMnmTOgMnVdGntmYiTa7Zsc/bK/SOTlNYvMN6W\n6gPzhybvPWsaLiRCxcgUIYTMPNESxkXkEgALSql/kfA9E8bLEprsbZOUcF20MrlB5iYNjK8NtyRB\nSML3xvXaOPmu+aUlwH884+93V18bJ+n7vw8lKedp2VJtnAw33jiTeVFMGCeEzDrn1NlZv98/e9zr\n9dDr9ersvr3YZmQn8q18Ozj8cQ1KknFKOt8lxDgBo5IE7qo/GzdKtXc/sG8OuO3k5DUvSDBOhv7q\n9O9DWLfOvzLPRKbMd1UZpxav8hsMBhgMBtO+DUIIqRRGntpE3shTyP51Se2b8/O2YbeFjPay7nMN\ndK6TOqlXB9pJ7g9CT0EaYu1XFxJdqnI1XpGSCjXDyBMhZNaJWapAhj+kKkILWxqS9q9zc6d8uUfm\nfDeKFMoCxu/X156vXfv4Pmt1nZ2rtYAw4xRaksAmJO9p3bpq8p5Y2ZwQQlpBrFIFnwDweQAvEJGH\nReQtMdrtHKEJ30n4IjZJ/aStqHMNzF1zk0nkbhtp7fneuxsR++4prY+1S8dX161dOl6d3FeSwMaX\n2G2bkU0bixmgKg0OK5sTQkgriLXa7t8qpZ6nlDpXKXWxUur2GO3WSlljE0KeyE0V7fkMTFrtJNtQ\npa3Qc+9jrfVqjn1FN9P6WDg9Wl135Qb93t1yZkVCfpPP4NhmatNGYOXK8e/TrrWp2uCElFQghBAy\nVVhh3BDb2NiE7sEW0k5oez6z46vWfd5mYHVfV/v2tbeQcOy7j23Oe1jHps2FhPv19WFKLpiVd/ae\nfVkr6myDA4wbImOcAGBxcdwA+cyRa6CqNjiMOBFCSKPh3nZV7uHmUqTUQNr1Se2Zz9PKEZjPl+/W\nBuqpPdpAFbm/pPtyk8Z995f2t0/6fsVqYMdV/rIEdsL1gQMjI2J/DujjxUXgjr3+1W3msxYkcNcN\nE8YJIbMOI09Ze7jF5MGC1yVFmrIiQSHYBTXzJqQb3OsWrFf7d3YjUy9IaM+ubu7+tzE5Tsp6NQUx\n3Sk3GxMtOnBgdHzH3uSilybixARuQgghDjRPhqLGIQ+h+8m5JBm8rE1006qH21NnpmZTUePouw/z\nav/Olzn3l1TtfKIa+fDVGKf3v3oUcZK+jkIB2flIH//GyCiFmKOq8ptowgghpNVw2q4OYk0NhtZZ\nKls9PKmdkP5N5XG3Lxu3RlPSNJ593q6+Nh2mgOXiInD5O0bGySZpGu471wE7hv8PmvpNIdNyMYtW\ndmAakNN2hJBZh5GnOog1NZinqrhhjfO5PY3mRqBcjJFZ47xPwl655+5zZ453KOD23aN23WlG++9k\nola7+qPpNRMpWrkSOHar/z5cozN/SBsnAJgXbbrMOXbyt28/PF97ReE0ICGEdAKapzqpY2rQxTVI\nIVOHrqlZi+zVffaGvudtHm3oa7ZmkTlgMD/+vWsqTfvu6/WfncxhArSBCjEgO64aTfPtUMDed4x/\nb8xR/3B2W2VImgakiSKEkFZB81QnVSaju9gGaC3Gk7UNSVEwuySAi++aNdBTdaYu01N7JqfubjsJ\nPH1Ul0Wwv4fVl89MSR/49R8fNx2Li6Pj0KjQsqUj8+ZGfeYP+RPQk9opg1vmIClhnRBCSGOheeoq\nrgG6zHq/FdlRsKQpurSK5fbKPXMPthnaux/o7Rh9bxfQtFkA8NgRbXZcU7FkiV4ll6fOkslZ+sgG\n/d41XXZkyk5Ad4lldOyIE6fxCCGkdZwz7RsgFXIQuhyAvcrN/s5F5oA1p4A1TjTHJJabyJW99Yox\nQDvNeScn213wnAvr2E5clzlgw3Zg+yrg8CpgFYC77x7Pddq0URuoENwEbbv2k0s/oWI5MGl00toJ\nxUzjmftjcUxCCGkFjDx1ETfHyRdlcqNHJtn7S6eA9wwTuk2E6qBzvr15cEgifNZWLHbBzhVP6BV1\n9jTdNdeMtxea6+SL7KQZlKSIE5BdtiBtqi8NbsdCCCGtg+api5gkcd9GumudV2B8E15g/DgkTytr\nCjBkuxc74dxEdm68cTzqZAiN0sSu05RmdMokmzPiRAghrYLmqQ2E7IW3BpOGyc5xAvwr2dZAT9M9\ntUcnc6/u6+X89veGouUWkvKn7DZ2XjtpdNatA669FjhxYvR5knm5/rPj70MjU3nxRZxCk80JIYR0\nAhbJnBahBS+BybwgX1trMW6StlrfjZkUjPadM2UEbN42LHKZtF+d776T9tjzFeVc6znXVA4HRlNr\ny5bq5GzDhz4EPH5iso/5Q7re054rgI1HtLEyeU6Li+ObAJvCmIC/WGWZYpi+ffY6CotkEkJmHUae\npkVWwUkgeerNPcdOxF6D8ekw30o282rnIRluOzneRtp+ee5325zv3RV/az3n2sYJGJkXe8rtxAng\n3e/2r3LrH9bG6f9v735j5bjOOo7/nhA7ShwwQRFpVLcFVEjjJIL2RcQ/4XVoIW3jS1BjOwUJAW+i\nJoXSlAo1Dty9hYTSQFPUNCER/5EqJ7alYCdNqat4DQVRoG0ax3ZDEdCGCoqA9qI4VXxpDi/OHs/Z\nuTM7M7szc3d3vh/J2ru7M3PGadT7yznPPEeSDlwj3fXJ0eLyrLYGWbVQ0z5JFxeb89QcACw0wlPb\n8jb5zZIVPuLj87Zbyeswrui7x+Rrncb1fgqNLvOW69JF4FkNOPNC4hVaH5zSDh6S7rlH2rrVv7/6\namnLRf7neLks2HtS2vfDyTLf6mp2W4N0LVS4dngdF37yvgvF5vRtAoCFR3hqW9WtWuIANS685C3r\n5YWXS/ZLF9+R8+XwvMufT87Pe2LvBvktX9bVUWl9MEsLNULj6oTOvDA6A/Xud/tgcvK+ZJ+6Fed7\nQj36qA83IWxt3Zq/FBcXf5ctLC8KRvRtAoBOIDxtlCpbtYzr+B2uldX1O2+GKzzZdryfBBBFx4bz\n4sLxLLuGf+JNfuN7SXcLj/e3k5KGlEVPqmXNQIWQErqGP7DHh54qT9jF3xW1DCgTjOp+ug8AMJMo\nGF80t2s0zOQVc1+y3weorz8s/dzNo0tuobB7xfkAlT4/vVyYFgrSY89KuneLr6lacf6pvr/6K+lX\nn0yO6e8Y32sprxD89Gnp4UdGj52m+Lto/LjIPEsTY88QCsYBdF0t4cnMrpf0QfmZrD9wzv1WxjGE\npzakw9K4p/ps+GRdOC4ORP8l6c4t0pvP5J9ftCwXhKf04us/8jZp72U+kBy4xgeqolAiJcEkswgK\nYwAAGw1JREFUBBnJh6myHcenVUcwmvNwRXgC0HVTL9uZ2XmS7pP045KukvRWM3vNtNdFRfHGv/HT\nbONqqly0lUq6turS4ffjzj8y/JNellP0Pt70d/cJ//Pekz44haWwrA17s5bF4tBx+HDyediypaxx\ntUhFfZqmDT0UlAPA3Kuj5ulaSV9wzn3RObcmab+kn6jhuqjiHzUafj4w5thYeKJO8gEnq/g7S5jR\nSveQyjpGkv74riQohe1Vzq75IvDjfX9MqBPas9sHjDgQpUPH2bWkDYE0umXLuHBUFF6m6RRehIJy\nAFgIdYSnl0t6Lnr/b8PP0KZQvP1shXPCXnKX7E8++4DKPQ0YluDiGa9YfO5D/dGgtLqaNMIMReCS\nn03aszupZwozSnmh45ED6/s4jQtH7/3L/PDSRqdwCsoBYCHwtN28S7cEuELltnOJ95K7cO/oDJSU\n/3Rf+im+K+S7kkvJMl18bujllNdCIN1vKS4Ej9/nhY64j1PRzM6vPrl+vGBlZ9IhPDwF2AQ2AgaA\nuXd+Ddf4sqRXRu+3DT9bp9/vn/u51+up1+vVMHzHhaWzEGiynqyLhaU0d8Y/aXf9aeljV47WP4Xr\nlhlP8k/QffVSSf+VfV68RLd1azLzJPkQ8cQTvp/T5k0+1Fx9dXJsCEsHD/mZqazZmnRX8vBEXPh8\n+ViyHHfgGunXr5P0u36GKl2kHncKb8oCzzgNBgMNBoONvg0AaNTUT9uZ2TfJLxb9qKR/l/R3kt7q\nnDudOo6n7ZoUZpuKmm7m7VlX1e0abT2wbMm+eMFD/eJ95sL3cVgKIanMU2lZx6Q/C/vkLZu/13vu\n8c020/eCUnjaDkDXTb1s55z7hqS3S/q4pJOS9qeDE1qQLt5Oy9qHruw2MfE1givkg9I/X5k00wz7\n4kk+OMXLaFn7zMXfb93qez+F93nipbi8+qZw7fj9qVO+WP2ZZ0a7llN7BACoiCaZXVN15il+Yu5B\nJQ00gyOSHh/OQIXrxHvWhQLw0IspzAqF1/D96qr0gW/17QwOHspuSBnPYj3+eP7sUfrc+Ly4H1TW\nDFWdQWrO+znlYeYJQNfVUfOEOoxrZlmndCF4eB+PHy8BhqAUvz4rP/N0LnSdSa4TZpzCklwIRtu3\n+/dxmAlLddtvk3bf788PNUkvRb2fQk+neBYrXgKMZ4/SReNPPDF6Xhxo4mBTtnt4WXVfDwAwM3ja\nrmlllsMkH0rKHlv12rGsgBTGl/xTd6HtwYPRd6EFwi3ywelIavzHJG27LllK27N7dElO8u+3XLT+\n81Mf9o0zJT/zdPK+0XvetGl9X6dY3DAz/QRfmSW6uvsv0c8JABYa4alp4/aAk0ZrkaoGqKJrlzk3\nXQv1e2fWHxsCU9xAMwSscI3/OeFbAWTVOJ09619XV0fDzOqqfw3F5INlf/6VVyZjr6765bmb3jLa\n1ynICkRhI+HQvqCoPUDd/Zfo5wQAC42ap6ak93I7ovxluSrHTnL8uHMf2yR9+mzyJJqUFIDfIuku\n+a1a4rGKQttgWXpgjw87Wy7Krk1K1z7Fe9UFn/+89JrXZJ8r+RmpMy8k36evFT/BV2bpjJqnUqh5\nAtB1zDw15TGV69Qdjj1S8tig7LWL7uvImu/3NFiW/u9L/vNPvCW5n0uj41ec9Nf7R/ezS1tx0u/t\nTWqc8pbN0q9hdigce+qUtP/h/HOXlpIZKal4ubDM0lndQWcBgxMAgJmn5lUpBC97bHhCbpoi8/S5\nuzb5IGXDXk3pGapgsCzt/21/zOU7pN/8MenFfT407ehLD96cP1tUpXeTNFoEntW3KUj3bQrLgFVn\nnlAKM08Auo6n7ZpWtPwWf18UhOJAE9oGVJW1oa8k3bDmrxeaXMadxG8Zvq446euH/DHn2hGsSY+8\nzf/4wB7pkYdHnzSLZ4vKPoF2440+/Fx5q7T3svWhK91J/I1vTM4N1073eQIAoCbMPG2kSbt7T3pe\n1rnpGaZn5TcHjr9/bPj61welC9+SHX6KeialZ4vyunrHGwMvm3TBXX6JLit0hYCUdd2mWgUsaB1T\nFcw8Aeg6ap7akH6CLv2EW9WWA5PMOJXtKH7F8DVuZXCDfCPMC4f1RVk1RPESW/w+/r5My4Dt2/3y\nXyhaf3GfdNvBZNwtFyXHnl3L3ly4qVYBeR3NAQCdwsxTG/JmiuLP000q626YGW/km3cvabdE390i\n345g3GxOmdmeopmbuBv4w1clTwDu6EvH+/7ndC1T1tN1Ur0zT2Vnzqa5/pzMaDHzBKDrmHlqUtEM\nUzyDtCvn5zrvIf4sLesJugc1ev9XvT2/X1LZ2Z6igBCeunv0UV+cHvRWkp+Lmm6GruTjejtJ0vKx\n8feS1lTvJma0AGCuEJ6aVNSuICyJxQGl7HJe2aW++B7y7kPRZ3GIuiX62frSys785bas5bOs48oI\ny3GH3zm6hLds/n1orpnuIB5/Hq4RB6QwfnjtHx9/H+G4EG6k4kBWFd3IAWDusGzXhjLLcPESXpmC\n8KpF4+ntWMZ9nz728h0+OGUtLT38Fen0/aNPuYWn5bI29c1bRisqMF82H+CKGmSmr+P6/rx4aW/p\n3tGZrP7w7xeLlw9DHdbxfv3LdfFYc9JSgWU7AF3HzFMbytQvHcn5OXaDJi82z2pPEIu3W9kVHftQ\n3weLrKWlzZukPQ/4n8OsSQhO4bMtFxXPrmRdOz2T1d/hf77xRn9s3CAzHB+/Lh/zwUnyr6HofOtW\nH4JCLdXtX1sfnOL7DdvM9Faa22qlaPsYAMBMITzNiqwZn7QQaibtLp4VtNI1UXEwu3wYWLLCz/Ix\n/ySc5GeFrrzV/xyCk5TsS7e0lF8vNC5YxaFiZae0d8/o9ePj04FsZaefcZL864dvSu4pGCz7MBWe\n4Avi4Hblrb5wXZIOXFO9TqqsOSkWBwDQJHO2xDM+sazmmJO0K8i6ftwMU0oaYobgIa1vSnl2zXcW\nt0GyJ97dd/tZoVhcwH333b6IOx0SlpaSn7NmdsL7PbtHNwyOjVv2CjNWBw8l42/eJN12m5+BCpsV\np507/jLpwAkfnHafGD87NEdPzAEAJkd4mgVZ4SjdeTzu9l339aVkw9+H+tLlGbMr6fARAtVgOZml\niWeFpKReKHyfDhbxrJPkr58l9H8KsloShNd0QIuX5OKlvXs/6GecsoJTfHy4x8Fy9vWDOatbAgBM\njvA0C8qGo0lmm8ZdPxSy366kOabrS4NXSStaL6tz9+HN0tlhoXiYnZJGC8iXlnydUjpYZM1oSetn\ncMLs1I7+cAuYA6PHZF2jjHHBKX2P6o9v7jkuwAEAFgpP282SOppjjrtGuhHnLiVh6pL90vtuHr80\nldcoMg5UYfYoDkJFzSXjIJQOZ/H5yybtPenDU9a95QWWsstp01yjQzNPPG0HoOsoGJ8ldXQVL2qw\nGZ7Yi5fxdm2SLtwr/euO4gaX6cLv9KzLpihkLB/L35Yla+PerOLxs2vS9tuSXk8PXyXd9cnse8tS\ntgFl0XFlm3sueHACAEwZnszsJjN7xsy+YWavq+umMIEyLQx2abQlQfAPP+iDzZ8Mipe+0iEhDkfh\n6boQQEITyvQ5eUElL2jt+XY/4yT51gJrR4vD0OZN5RtQ1tWokqU6AOiEaWeeTkj6SUkFrZrRuHQn\n8XgWK2uLltjKziTgpIu2s4JEOiQcPCTdc0/ydN1tB0d7LIUZqHC9cUEl7z5eeilpLZA+N69v1NJS\n0ppgdTU/3JTZtLgIncEBoDOmKhh3zj0rSWZGTcMsiGeb4qfqstoRpIvE+zukk/eN1u1UqeM584Lv\niXT6ft9Taacl3b1jZYrEl5bW1z1dfbV0PHoyL5ybVSMVh7Ng69bxdUvx04RVZd0Ds1AAsLB42m6R\njHuqLvagpEfeJu2/f3Tbkzui0PHEE9WfINvzgHTBtuS40GNJGg0U6bYHcVg6fHj9uGn33OPDWt5T\nbnE4C9+VmVGadMYpvofzzlu/NQ0AYKEULtuZ2VEzezr6c2L4WlSajDbFASnd0iDdmfyCu/wMkeSX\nwvbsXr90FW+4m7fRbxBvhfLivqQLd+ixlLf9Svg8XXCednZtdPkttBjIW26La6yyCrnrXGKL7yH0\ntQp/F5byAGAh1dKqwMyOSXqXc+4zY45xy8vJI869Xk+9Xi/72EVqVVBH+4EysnpExU/VST5UfXq4\nCe6e3aONJ0MLgawNetPLaHmylunyWhWkP5fylwvTx546NdquoMoyWVMtBcI9dKBlwbhWBYPBQIPB\n4Nz7lZUVWhUAWDh1hqdfds59eswx3ezzNGlX8DJsi/TmM+sDUjqshXt4qD/6eQhQ437Rl+nTFCwf\nW7/JrpQfKNL9odKhLX6fF/aqqPJ3mcaC1zzR5wlA100VnszsRkkfknSppK9Jeso598acY7sVnrJm\nfcbNQFWdobpkv+/N9PWHpa/ePD6k3SBpqT/6WfgFX+YXfR2zKXnjVPk8Hfbywtq463ZgZqhphCcA\nXUeH8aaVnXmqMkNlW6TLn0/e//vFfgYqL3ylZ5wmCRBtzqaMu7/4PrKWCYvOT18DlRGeAHQdHcbr\nktWUUirej65Mc8s0d8bPOEn+1VUITll9lsoUNtcVNorGKrq/s2ujBeqhj9S489PmIThRbA4AM4uZ\np7pMW9tU9fwbJD2+xQenPOngFLzzl3zfo9VV6bnn8uuO6lZ2xis+Tso/Z9KZp1k34/fPzBOArmPm\naVqTzBxlKZqhStul8cEpz+ZNSZfurVtHZ2nK7AEnjc70hGuWGbfsFihxl/Fx58R9pLLOn8HgUWiS\nrWKYpQKAVhGephX3T7pF5Yq+swJW2WLxsmEtb9Yp3Rsp3cOpzC/s/vHkmLx96qb9hX52rXjblLxi\n8XD+PKq6VUzZjY8BALVh2a4uVZ6Wq6N9wbhr5AWnWFw0XbaX0/KxZLNfSTpwq3Tqw8n78Oh/mdYE\n4fMyxdt1FnhnXWsWC8jL/nNpo/VCCst2ALqO8NSmqu0Liq6VdW6Z4JSn7C/sF/dJK8P/LUNX7byG\nlulf6PEYVWp7Jg04RePNeH1RoQ24f8ITgK5jb7s25e09N+m10rKerKvSW6lMODm75vfFk7I3ws3b\n+Dc9Rt6+dFkmDQjpJpxZe+ZV3b9v1kyzoTEAYCKEp41QtTi8jLK9nOqYqdh7mbT57tHAFCvzC70o\nZAVVQlbReVnjlbmHWTev9w0Ac4plu0WQNeNUZj+5lmpkxkrXXtXZFbxsjdUs1jzNMJbtAHQdM0+L\nKG9Wp+xsTxOKlgrHBaRJl6ayziuzXEmYAgCMwczTvHuoX8++cZMoe50y26XUMSNWx99r3gvIW8DM\nE4Cuo8/TPHuoP77PT16QqCM47dldrr9Q2e1S4t5Gk6ij39EkDSoBAJ1DeJpXYcZpI37Z79nt2xOU\nGbds08fDh/2Mk1QtBIW97+r451C1QSUAoJOoeZpHoUA8r4ap6vJVleM3b0qCk+T7PBWdW1SzFJbK\nQs8oqdyTdfESW/qfw6RLljz6DwAoQHiaN+kn69K/7KvW7FQ9/uya31A4bCz8yIFy912mFcH27aNN\nN8cFmPRsU9gLr2qH83CtoqJyAACGWLabJ+P2q5OqL19NugltvLFwHXvYxUtljxwot6lv1hJbmHHK\n+jvlfc7ecACAighPGy1vY9+0MtuuVK3ZmaTGJ++caULUwUOjgSnrPpaPFZ837v6yPp+mVopicgDo\nLFoVbLQyW7VMstFvUzVPWeeMW/qrqy2C60vWn+z+xn0+SWuCjrczoFUBgK5j5mmj3CAfnDR8zZuB\nKhOc0ktPVcPKJOGmzFJh+r62XFR9nOVjPjhJ/jVrBmrc/RV9njV7NQ7tDACg8ygY3yhlNgkuO+M0\ni5vbpu/rFa9Iiszv/WD566zs9H+qzjxVUeWf10Z2aQcAzISpZp7M7P1mdtrMnjKzQ2b2LXXdWGfk\nbRJcJjhJG9+bqEyN0enTo0Xmk8xA9XdMf691qTpbBQBYKFPVPJnZ6yU96Zx7yczeJ8k5596Tcyw1\nT2WVDU6xjd6PrajG6J2/NNnME2YONU8Aum6qZTvn3Ceit38riee9pzVJcJI2fvmoqMbo3g/6Gacz\nL7R3TwAANKDOgvGfl/REjdfrnkmD07yIg1PZwm+JomwAwEwpDE9mdtTMno7+nBi+7oqO2SdpzTn3\nkUbvdpEtenBK6x8vd1zTTSwJZgCAigqX7Zxzbxj3vZn9rKQ3Sbqu6Fr9fv/cz71eT71er+iUbpiH\n4DRNTVU4d/Mm6T0fT4KT6/tC8JWd+ec1+SRhx/s1NWEwGGgwGGz0bQBAo6YtGL9e0u9I+hHn3H8X\nHEvBeJZ5CE7ThIxwbtgPL1yjbOuBacYeF/g2b/IzWsHdd/vXja4dmwMUjAPoumlrnj4k6WJJR83s\nM2Z2fw331B3zEJym3cIknBtaFYRrlG09MGlbgKLlvnSLhaUl9rgDAJQyVXhyzn23c+5VzrnXDf/c\nWteNLbx5CE7SdH2k4nNXV0evkbdUl3edKsoGvhDMDh+mazgAoDQ6jG+EeQlOwcFDk9cbxee21Yuq\nShfw8B1dwwEAJbEx8EaYt/A0r9rYILmDqHkC0HVsDNw2glN72tggGQDQOYSnNs1qcKqjxoc6IQBA\nR1Dz1JZZDU519DqiXxIAoEOYeWrDrAanadoQ1HkNAADmCOGpabManKTp2hDUeQ0AAOYIT9s1aZaD\nU2zap8zCbBPBqRN42g5A11Hz1JR5CU7SdKGnjXonWggAAGYIy3ZNmKfgNI026p2KtlkBAKBlhKe6\ndSU4Sc3XO1GMDgCYQSzb1alLwSmYZuuWIlW2WQEAoCUUjNelreDUxfqfLv6dZxgF4wC6jpmnOrQV\nnLrajJLgBACYIdQ8TavNGSfqfwAA2HCEp3lBM0oAAGYCNU+T2qjicOp/sMGoeQLQdcw8TWIjn6oj\nOAEAsKEIT1V1sR0BAAA4h/BUBcEJAIDOmyo8mdl7zexzZvZZM/uYmb2srhubOQQnAACg6Wee3u+c\n+17n3GslPS6pfCXpPCE4AQCAoanCk3Pu+ejtFkkvTXc7AAAAs23qDuNm9huSfkbS1yTtnPqOZg2z\nTgAAIFIYnszsqKTL4o8kOUn7nHNHnHN3SrrTzH5F0i9I6uddq99Pvur1eur1ehPddGsITkAlg8FA\ng8Fgo28DABpVW5NMM3uFpI86567J+X6+mmQSnBI05kSEJpkAum6qZTsze7Vz7p+Gb2+UdHr6W5oB\nBKdEVzcjBgAgx7RP273PzJ42s6ckvV7SO2q4p41FcEqwGTEAAOtMNfPknLuprhuZCQSnUWEz4jDz\nxNIdAABsDHwOwSkfNU+IUPMEoOvYnkUiOBUhOAEAcA7hieAEAAAq6HZ4IjgBAICKuhueCE4AAGAC\n3QxPBCcAADCh7oUnghMAAJhC98ITAADAFLoVnph1AgAAU+pOeCI4AQCAGnQjPBGcAABATRY/PBGc\nAABAjRY7PBGcAABAzRY3PBGcAABAAxY3PAEAADRgMcMTs04AAKAhixeeCE4AAKBBixWeCE4AAKBh\nixOeCE4AAKAFixGeCE4AAKAltYQnM3uXmb1kZt9Wx/UqITgBAIAWTR2ezGybpDdI+uL0t1MRwQkA\nALSsjpmneyW9u4brVENwAgAAG2Cq8GRmS5Kec86dqOl+yiE4AQCADXJ+0QFmdlTSZfFHkpykOyXd\nIb9kF3+Xq9/vn/u51+up1+uVv9OA4ATMrMFgoMFgsNG3AQCNMufcZCeaXS3pE5JekA9N2yR9WdK1\nzrn/zDjeTToWgPlkZnLOjf2PKgCYNxOHp3UXMvsXSa9zzn0153vCE9AxhCcAi6jOPk9OBct2AAAA\n8662mafCgZh5AjqHmScAi2gxOowDAAC0hPAEAABQAeEJAACggpkMTxvRJ4YxGXOexuvSmAAwawhP\njMmYczhel8YEgFkzk+EJAABgVhGeAAAAKmi1z1MrAwGYKfR5ArBoWgtPAAAAi4BlOwAAgAoITwAA\nABXMfHgys3eZ2Utm9m0tjPVeM/ucmX3WzD5mZi9rYcz3m9lpM3vKzA6Z2be0MOZNZvaMmX3DzF7X\n4DjXm9nnzewfzexXmhonNeYfmNlXzOzplsbbZmZPmtlJMzthZr/YwpgXmNmnhv+enjCz5abHHI57\nnpl9xswOtzEeAMyqmQ5PZrZN0hskfbGlId/vnPte59xrJT0uqY1fSh+XdJVz7vskfUHSe1oY84Sk\nn5R0vKkBzOw8SfdJ+nFJV0l6q5m9pqnxIn80HLMt/yfpdufcVZJ+QNJtTf89nXMvSto5/Pf0+yS9\n0cyubXLMoXdIOtXCOAAw02Y6PEm6V9K72xrMOfd89HaLpJdaGPMTzrkwzt9K2tbCmM86574gqcmn\noK6V9AXn3Bedc2uS9kv6iQbHkyQ55z4p6atNjxON9x/OuaeGPz8v6bSkl7cw7gvDHy+QdL6kRp/8\nGP6HzJsk/X6T4wDAPJjZ8GRmS5Kec86daHnc3zCzL0n6KUm/1ubYkn5e0hMtj9mUl0t6Lnr/b2oh\nVGwkM/sO+ZmgT7Uw1nlm9llJ/yHpqHPu7xseMvyHDI/nAui88zdycDM7Kumy+CP5/3O+U9Id8kt2\n8XdNjrnPOXfEOXenpDuHNTq/IKnf9JjDY/ZJWnPOfWTa8cqOifqY2cWSDkp6R2oGsxHD2crXDmvk\nHjWz7c65RpbUzOzNkr7inHvKzHpqdsYSAGbehoYn59wbsj43s6slfYekz5mZyS9lfdrMrnXO/WcT\nY2b4iKSPqobwVDSmmf2s/JLIddOOVXbMFnxZ0iuj99uGny0cMztfPjj9mXPuz9sc2zn3v2Z2TNL1\naq4e6YckLZnZmyRdKOmbzexPnXM/09B4ADDTZnLZzjn3jHPuZc6573LOfaf8ks9rpw1ORczs1dHb\nG+XrVxplZtfLL4csDQuB29bULMLfS3q1mb3KzDZLullSW09pmdqdHflDSaecc7/bxmBmdqmZbR3+\nfKH8DO3nmxrPOXeHc+6Vzrnvkv/f8UmCE4Aum8nwlMGpnV+G7zOzp83sKUmvl3+6qGkfknSxpKPD\nx8Dvb3pAM7vRzJ6T9P2SHjOz2uusnHPfkPR2+acJT0ra75xrI4x+RNLfSPoeM/uSmf1cw+P9kKSf\nlnTdsHXAZ4aBuEmXSzo2/Pf0U5L+wjn30YbHBAAMsT0LAABABfMy8wQAADATCE8AAAAVEJ4AAAAq\nIDwBAABUQHgCAACogPAEAABQAeEJAACgAsITAABABf8PXs/YclTSe8EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11233cc88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Fill in the list of models and list of names below (3 each)\n",
    "plt.figure(figsize=(10,10))\n",
    "for j, model, pTitle in zip(range(221,224), \n",
    "                            [LR_vanilla1, LR_vanilla2, LR_findMin],\n",
    "                            ['vanilla 1e-1', 'van 1e-2', 'findMin']):\n",
    "    plt.subplot(j)\n",
    "    binaryClassifier2DPlot(model)\n",
    "    plt.title(pTitle)\n",
    "    plt.axis('tight');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setting alpha\n",
    "\n",
    "## From above/last time\n",
    "\n",
    "Last time, we used the approach with $\\alpha = \\alpha_0 = \\alpha_1 = \\ldots$ and we saw that the results obtained were highly dependent on the [sequence of] learning rates.\n",
    "\n",
    "## Lipschitz constant\n",
    "\n",
    "If we can determine the Lipschitz constant of our objective function, then we get a very nice progress bound under certain assumptions on the objective function.\n",
    "\n",
    "Namely, assume that $f$ is strongly smooth. In particular, there exists $L$ such that for any $v, z$, \n",
    "$$\n",
    "v^T \\nabla^2 f(z) v \\leq L\\|v\\|_2^2\n",
    "$$\n",
    "For some $z$, it then follows by Taylor's theorem that \n",
    "$$\n",
    "f(y) = f(x^t) + \\nabla f(x^t) (y-x^t) + \\frac{1}{2}(y-x^t)^T \\nabla^2 f(z) (y-x^t) %\n",
    "\\leq f(x^t) + \\nabla f(x^t) (y-x^t) + \\frac{L}{2}\\|y-x^t\\|_2^2\n",
    "$$\n",
    "Setting $y = x^{t+1}$ to be the minimizer of this equation yields\n",
    "$$\n",
    "x^{t+1} = x^t - \\frac{1}{L} \\|\\nabla f(x^t)\\|_2^2\n",
    "$$\n",
    "In particular, gradient descent with $\\alpha_t := 1/L$ minimizes the upper bound (*i.e.*, maximizes the *worst-case* progress). Substituting this into the expression above, \n",
    "$$\n",
    "f(x^{t+1}) \\leq f(x^t) - \\frac{1}{2L} \\|\\nabla f(x^t)\\|_2^2\n",
    "$$\n",
    "*Caveat* *via* a direct quote from Mark Schmidt: \"In practice, you should **never use $\\alpha = 1/L$**.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def logisticLipschitz(X, lam=0):\n",
    "    return .25 * np.max(np.linalg.eigvals(X.T @ X)) + lam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lipschitz constant is 127.45, implying alpha is 7.846e-03\n"
     ]
    }
   ],
   "source": [
    "objectiveParms_Lip = {'lam': 1}\n",
    "L_logist = logisticLipschitz(X, **objectiveParms_Lip)\n",
    "\n",
    "print('Lipschitz constant is {:5.5g}, implying alpha is {:5.3e}'.format(L_logist, 1/L_logist))\n",
    "\n",
    "# construct a model LR_Lip that uses the 'correct' value of alpha on each iteration\n",
    "LR_Lip = LogisticRegressor(objectiveParms=objectiveParms_Lip, gdParms={'alpha': 1/L_logist})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 3.18 ms per loop\n"
     ]
    }
   ],
   "source": [
    "# fit the model and time how long it takes\n",
    "%timeit LR_Lip.fit(X,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "40"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the number of function evaluations\n",
    "LR_Lip.funEvals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive step-size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This approach often gives bigger steps and faster progress, but step size never increases.\n",
    "\n",
    "1. Start with a small guess for $L$  \n",
    "$ L \\leftarrow 1$\n",
    "2. Double $L$ if the *progress inequality* from above  \n",
    "$f(x^{t+1}) \\leq f(x^t) - \\frac{1}{2L}\\|\\nabla f(x^t)\\|_2^2$\n",
    "is not satisfied."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that this \"backtracking\" approach makes sense when it is cheap to evaluate our function and expensive to compute the matrix product $Xw$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Look at Armijo below before writing this section of code. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def adaptiveGradientDescent(funObj, w, X, y,\n",
    "                            objectiveParms={}, \n",
    "                            gdParms={}):\n",
    "    # optimality tolerance\n",
    "    optTol = gdParms.get('optTol')\n",
    "    if optTol is None:\n",
    "        optTol = 1e-2\n",
    "    maxEvals = gdParms.get('maxEvals')\n",
    "    if maxEvals is None:\n",
    "        maxEvals = 500\n",
    "    verbose = gdParms.get('verbose')\n",
    "    if verbose is None:\n",
    "        verbose = True\n",
    "    # Armijo Rule.\n",
    "    # Set the default value for gamma for the Armijo rule\n",
    "    gamma = gdParms.get('gamma')\n",
    "    if gamma is None:\n",
    "        gamma = 1\n",
    "    # initial objective and gradient value\n",
    "    f,g = funObj(w, X, y, **objectiveParms)\n",
    "    # this took a single function evaluation\n",
    "    funEvals = 1\n",
    "    \n",
    "    # initialize the step size\n",
    "    L0 = gdParms.get('L0')\n",
    "    if L0 is None:\n",
    "        L0 = 1\n",
    "    # what should alpha be set to?\n",
    "    alpha = 1/L0\n",
    "    \n",
    "    # Gradient Descent loop\n",
    "    while True:\n",
    "        # update w \n",
    "        w_new = w - alpha*g\n",
    "        # update f and gradient \n",
    "        f_new, g_new = funObj(w_new, X, y, **objectiveParms)\n",
    "        funEvals += 1\n",
    "        gg = g.T @ g\n",
    "        while f_new > f - gamma * alpha * gg:\n",
    "            if verbose:\n",
    "                print('Backtracking...')\n",
    "            alpha /= 2\n",
    "            w_new = w - alpha*g\n",
    "            f_new, g_new = funObj(w_new, X, y, **objectiveParms)\n",
    "            funEvals += 1\n",
    "        # Update parameters/function/gradient\n",
    "        w = w_new\n",
    "        f = f_new\n",
    "        g = g_new\n",
    "        # Test terminate conditions\n",
    "        optCond = np.linalg.norm(g, np.inf)\n",
    "        if verbose:\n",
    "            print('%6d %15.5e %15.5e %15.5e' % (funEvals, alpha, f, optCond))\n",
    "        alpha = 1/L0 \n",
    "        if optCond < optTol:\n",
    "            if verbose:\n",
    "                print('Problem solved up to optimality tolerance')\n",
    "            break\n",
    "        if funEvals >= maxEvals:\n",
    "            if verbose:\n",
    "                print('At maximum number of function evaluations')\n",
    "            break\n",
    "    return (w, f, funEvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "LR_adap = LogisticRegressor(objectiveParms={'lam':1},\n",
    "                            gdParms={'L0':1},\n",
    "                            gd=adaptiveGradientDescent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 32.7 ms per loop\n"
     ]
    }
   ],
   "source": [
    "# Timed fit of the model above to the data [X, y]\n",
    "%timeit LR_adap.fit(X,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "506"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the number of function evaluations required by LR_adap\n",
    "LR_adap.funEvals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backtracking line-search\n",
    "\n",
    "Checks against the Armijo condition\n",
    "$$\n",
    "f(x^{t+1}) \\leq f(x^t) - \\alpha\\gamma \\|\\nabla f(x^t)\\|_2^2, \\quad \\gamma \\in (0, 1/2]\n",
    "$$\n",
    "See Wikipedia for a more detailed discussion on the more technical [Wolfe conditions](https://en.wikipedia.org/wiki/Wolfe_conditions) (of which the Armijo rule is a special case)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this version was implemented above in `adaptiveGradientDescent` the only step required for this section is to choose `gamma` and fit the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a model named LR_armijo that uses the Armijo rule\n",
    "# ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Timed fit of LR_armijo to the data [X, y]\n",
    "# ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Print the number of function evaluations required by LR_armijo\n",
    "# ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** We now have the tools to describe how `findMin` works. `findMin` uses backtracking to choose `alpha` according to the Armijo rule: it starts with `alpha = 1` and iteratively uses *cubic Hermite interpolation* to select the next `alpha` if the previous one does not satisfy the Armijo rule. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of the three"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "for j, model, pTitle in zip(range(221,224), \n",
    "                            [LR_Lip, LR_adap, LR_armijo], \n",
    "                            ['Lipschitz', 'Adaptive', 'Armijo']):\n",
    "    plt.subplot(j)\n",
    "    binaryClassifier2DPlot(model)\n",
    "    plt.title(pTitle)\n",
    "    plt.axis('tight');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coordinate optimization\n",
    "\n",
    "In this case, we must assume that our function $f(w; X,y)$ is the sum of a smooth function and a separable function\n",
    "$$\n",
    "f(w; X,y) = g(w; X,y) + \\sum_{j=1}^d h_j(w_i; X,y)\n",
    "$$\n",
    "Then at each iteration $t$ select an index $j_t$ and update that index via \n",
    "$$\n",
    "x^{t+1}_{j_t} = x^t_{j_t} + \\gamma_t \\quad\\iff \\quad x^{t+1} = x^t + \\gamma_t e_{j_t}\n",
    "$$\n",
    "for some constant $\\gamma_t$ and where $e_{j_t}$ is the $j_t$-th unit vector.\n",
    "\n",
    "In the case of logistic regression, one of the key tools here is that this allows us to track updates via \n",
    "$$\n",
    "Xw^{t+1} = Xw^t + \\gamma_t X e_{j_t} = Xw^t + \\gamma_t x_{j_t}\n",
    "$$\n",
    "so that it is only necessary to keep track of the index $j_t$ that is selected and the update constant $\\gamma_t$ — no more than one expensive matrix multiplication is required [to compute $Xw^0$]. But how does one go about choosing $j_t$ on iteration $t$? We'll mention two ways. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lipschitz and Uniform Sampling\n",
    "\n",
    "In the first selection method that we'll discuss, $j_t$ is chosen uniformly at random on each iteration from $[d] := \\{1, 2, \\ldots, d\\}$ – *i.e.*, $j_t \\sim \\mathrm{Unif}([d])$ and $j_t$ chosen according to $P(j_t = k) = L_k/Z$ for $k \\in [d]$, where $L_k$ denotes the $k$th coordinate-wise Lipschitz constant and $Z$ a normalizing factor so that $Z^{-1}\\sum_k L_k = 1$ gives a pmf on $[d]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def weightedIndexSample(weights):\n",
    "    \"\"\"\n",
    "    weightedIndexSample returns a random integer in {1, 2, ..., weights.size}\n",
    "    according to the probability mass function induced by weights\n",
    "    \"\"\"\n",
    "    return np.random.choice(range(weights.size), p=weights/weights.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def uniformIndexSample(v):\n",
    "    \"\"\"\n",
    "    uniformIndexSample returns a random integer in {1, 2, ..., v.size}\n",
    "    uniformly at random. \n",
    "    \"\"\"\n",
    "    # Generate an index uniformly at random, given the vector v\n",
    "    # Hint: \"return np.random.choice(...)\"\n",
    "    # ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def lipschitzLogisticLipschitz(X, lam=0):\n",
    "    \"\"\"\n",
    "    lipschitzLogisticLipschitz returns a vector of the coordinate-wise\n",
    "    Lipschitz constants for the logistic function given data X and \n",
    "    regularizer lam. (i.e., each L_j, j in [d] satisfies the Lipschitz\n",
    "    continuity inequality for coordinate j)\n",
    "    \"\"\"\n",
    "    return .25 * (X**2).sum(axis=0) + lam\n",
    "def uniformLogisticLipschitz(X, lam=0):\n",
    "    \"\"\"\n",
    "    uniformLogisticLipschitz returns the uniform coordinate-wise Lipschitz \n",
    "    constant for the logistic function given data X and regularizer lam \n",
    "    (i.e., the smallest Lipschitz constant for which the Lipschitz continuity \n",
    "    inequality is satisfied for each coordinate)\n",
    "    \"\"\"\n",
    "    # Which of the coordinate Lipschitz constants works for all of the coordinates?\n",
    "    # Hint: \"return FUNCTION(lipschitzLogisticLipschitz(X, lam))\"\n",
    "    # ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def randomizedCoordinateDescent(funObj, w, X, y, \n",
    "                                objectiveParms={}, \n",
    "                                gdParms={}):\n",
    "    \"\"\"\n",
    "    randomizedCoordinateDescent performs randomized coordinate \n",
    "    descent on the convex function funObj to find the optimal \n",
    "    parameter vector w given labeled training data [X, y]. \n",
    "    \"\"\"\n",
    "    # progress tolerance\n",
    "    progTol = gdParms.get('progTol')\n",
    "    if progTol is None:\n",
    "        progTol = 1e-4\n",
    "    # maximum number of passes through the data\n",
    "    maxPasses = gdParms.get('maxPasses')\n",
    "    if maxPasses is None:\n",
    "        maxPasses = 500\n",
    "    # whether to print output along the way\n",
    "    verbose = gdParms.get('verbose')\n",
    "    if verbose is None:\n",
    "        verbose = False\n",
    "    # Which Lipschitz method to use?\n",
    "    # default: uniformLogisticLipschitz\n",
    "    lipschitzMethod = gdParms.get('lipschitzMethod')\n",
    "    if lipschitzMethod is None:\n",
    "        L = uniformLogisticLipschitz(X, **objectiveParms)\n",
    "    else:\n",
    "        L = lipschitzMethod(X, **objectiveParms)\n",
    "    # How to choose the coordinate/index jt on iteration t\n",
    "    # default: weighted index sampling\n",
    "    # default: if L is a uniform Lipschitz constant, then \n",
    "    #          uniform index sampling is performed.\n",
    "    coordSample = gdParms.get('coordSample')\n",
    "    if coordSample is None:\n",
    "        coordSample = weightedIndexSample\n",
    "    if isinstance(L, float):\n",
    "        unif = np.ones(w.shape).ravel()\n",
    "    funEvals = 0\n",
    "    # For fast mulitiplication\n",
    "    Xw = X @ w\n",
    "    # For tracking progress\n",
    "    w_old = w.copy()\n",
    "    \n",
    "    for t in range(maxPasses*d):\n",
    "        # Choose index and Lipschitz coefficient\n",
    "        if isinstance(L, float):\n",
    "            # Pick j...\n",
    "            # Pick Lt...\n",
    "        else:\n",
    "            # Pick j...\n",
    "            # Pick Lt...\n",
    "        # Compute objective value and partial derivative g_j\n",
    "        Xj = # ...\n",
    "        f, g_j = funObj(w, Xw, Xj, y, j, **objectiveParms)\n",
    "        funEvals += 1\n",
    "        # Variable update\n",
    "        # ... for product ... \n",
    "        # ... for w ... \n",
    "        # Check for lack of progress after each pass\n",
    "        if np.mod(t,d) == 0:\n",
    "            change = # ... \n",
    "            if verbose:\n",
    "                print('Passes = %d, function = %.4e, change = %15.4e' % (t/d, f, change))\n",
    "            if change < progTol:\n",
    "                if verbose:\n",
    "                    print('Parameters changed by less than progTol on pass %d.' % t)\n",
    "                break\n",
    "            w_old = w\n",
    "    return (w, f, funEvals)\n",
    "    \n",
    "    \n",
    "def coordinateObjective(w,Xw,Xj,y,j,lam):\n",
    "        \"\"\"\n",
    "        coordinateObjective computes the negative log-likelihood of the MAP estimate\n",
    "        and the jth gradient for the logistic function using the matrix-product Xw,\n",
    "        the labels y, an index j and regularizer lam\n",
    "        \"\"\"\n",
    "        yXw = y*Xw\n",
    "        invsigmoid = 1+np.exp(-yXw)\n",
    "        nll = np.sum(np.log(invsigmoid)) + .5*lam*np.linalg.norm(w,2)\n",
    "        g_j = -Xj.T @ (y*(1-1/invsigmoid)) + lam*w[j, 0]\n",
    "        return (nll, g_j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "coordModels = {}\n",
    "for name1, coordSamp in zip(['weighted', 'uniform'], \n",
    "                            [weightedIndexSample, uniformIndexSample]):\n",
    "    for name2, lipMethod in zip(['coordLip', 'unifLip'], \n",
    "                                [lipschitzLogisticLipschitz, uniformLogisticLipschitz]):\n",
    "        coordModels[(name1, name2)] = LogisticRegressor(objectiveParms={'lam':1},\n",
    "                                                        gdParms={'coordSample': coordSamp,\n",
    "                                                                 'lipschitzMethod': lipMethod},\n",
    "                                                        objective=coordinateObjective,\n",
    "                                                        gd=randomizedCoordinateDescent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for key, value in coordModels.items():\n",
    "    print(key)\n",
    "    %timeit value.fit(X,y)\n",
    "    print('Number of Passes: {}\\n'.format(value.funEvals))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "j = 221\n",
    "for pTitle, model in coordModels.items():\n",
    "    plt.subplot(j)\n",
    "    j += 1\n",
    "    binaryClassifier2DPlot(model)\n",
    "    plt.title(pTitle)\n",
    "    plt.axis('tight');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## So which is better?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** For a more concrete/realistic example of which is better, see the exercises below. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that where uniform sampling is used - $\\forall j \\in [d], \\,\\, p(j) = d^{-1}$ - we achieve the guaranteed progress bound\n",
    "$$\n",
    "  \\E (f(x^{t+1})) \\leq f(x^t) - \\frac{1}{2dL} \\|\\nabla f(x^t)\\|^2\n",
    "$$\n",
    "where $[d] := \\{1, \\ldots, d\\}$. Correspondingly, for Lipschitz sampling, the\n",
    "pmf is given by\n",
    "$$\n",
    "\\forall j \\in [d],\\quad p(j) = \\frac{L_j}{\\sum_{j \\in [d]} L_j} =: \\frac{L_j}{d \\overline{L}}\n",
    "$$\n",
    "where $d\\overline{L}$ is the normalizing constant for the pmf, with\n",
    "$\\overline{L}$ the arithmetic mean of the coordinate-wise Lipschitz constants\n",
    "$L_j$. Hence, the guaranteed progress bound for Lipschitz sampling is\n",
    "$$\n",
    "\\E f(x^{t+1}) \\leq \\E \\big( f(x^t) - \\frac{1}{2L} |\\nabla_{j_t} f(x^t)|^2\\big) %\n",
    "= \\sum_{j=1}^d \\frac{L_{j_t}}{d\\overline{L}} \\big( f(x^t) - \\frac{1}{2L} |\\nabla_j f(x^t)|^2 \\big)\n",
    "$$\n",
    "$$\n",
    "= f(x^t) - \\frac{1}{2d L} \\sum_{j\\in d} \\frac{L_{j_t}}{\\overline{L}} |\\nabla_{j_t} f(x^t)|^2\n",
    "$$\n",
    "In particular, by comparing the two equations above, uniform sampling is on average preferred to Lipschitz sampling if"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\ip{\\mathbf{1}, |\\nabla_{j_t} f(x^t)|^2} = \\|\\nabla f(x^t)\\|^2 > \\sum_{j_t\\in [d]} \\frac{L_{j_t}}{\\overline{L}}\n",
    "|\\nabla_{j_t} f(x^t)|^2 = \\ip{\\frac{L_{j_t}}{\\overline{L}}, |\\nabla_{j_t}\n",
    "  f(x^t)|^2}.\n",
    "$$\n",
    "Firstly, by the Cauchy-Schwarz inequality, this certainly occurs if\n",
    "$\\|L_{j_t}\\|_2 < \\overline{L}$. Moreover, as an example, perhaps\n",
    "Lipschitz sampling is sub-optimal on its second pass through the data, where\n",
    "the gradient on the highly-preferred indices is already quite small. This would\n",
    "cause the two vectors to be closer to orthogonal, thereby being more likely to\n",
    "satisfy the above inequality. Note: more generally, it is desirable to have a\n",
    "sampling scheme $p(j_t)$ which is well-aligned with the squared modulus of the $j_t$-gradient vector\n",
    "$\\big(|\\nabla_{j_t} f(x^t)|^2 \\big)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import requests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare run time on larger data sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a test data set that includes 500 observations of 100 features (*i.e.,* $X \\in \\reals^{n\\times d}$ with $n = 500$, $d = 100$; correspondingly, $y \\in \\{-1,1\\}^{500}$). The labels `y` are stored in the very last column. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bigData = pd.read_csv(filepath_or_buffer='./data/logisticData.csv',header=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('bigData.shape = {}'.format(bigData.shape))\n",
    "bigData.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Xbig = bigData.values[:,:-1]\n",
    "ybig = bigData.values[:, -1].reshape(-1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "coordModelsBig = {}\n",
    "for name1, coordSamp in zip(['weighted', 'uniform'], \n",
    "                            [weightedIndexSample, uniformIndexSample]):\n",
    "    for name2, lipMethod in zip(['coordLip', 'unifLip'], \n",
    "                                [lipschitzLogisticLipschitz, uniformLogisticLipschitz]):\n",
    "        coordModelsBig[(name1, name2)] = LogisticRegressor(objectiveParms={'lam':1},\n",
    "                                                           gdParms={'coordSample': coordSamp,\n",
    "                                                                    'lipschitzMethod': lipMethod},\n",
    "                                                           objective=coordinateObjective,\n",
    "                                                           gd=randomizedCoordinateDescent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for key, value in coordModelsBig.items():\n",
    "    print(key)\n",
    "    %timeit value.fit(Xbig,ybig)\n",
    "    print('Number of Passes: {}\\n'.format(value.funEvals))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for key,value in coordModelsBig.items():\n",
    "    print(key)\n",
    "    print(value.w[value.w != 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use `cython` or `numba` to speed up the optimizations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `numba` implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numba\n",
    "from numba import jit, float64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Insert code below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `cython` implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "# Insert code below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## UMass Amherst smoking data set for logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could try testing out these methods on real data. For this:\n",
    "* Need to think about missing values\n",
    "* May have to think about using a kernel (argument may not longer be best represented as a linear combination of the input data)\n",
    "\n",
    "For example, we include below a data set from the UMass Amherst website. Many more data sets can be found on this site – and for more than just logistic regression. The ones specifically for logistic regression are [here](https://www.umass.edu/statdata/statdata/stat-logistic.html). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "descr = requests.get('https://www.umass.edu/statdata/statdata/data/nhanes3.txt')\n",
    "print(descr.text)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "smoking = pd.read_csv(filepath_or_buffer='https://www.umass.edu/statdata/statdata/data/nhanes3.dat', sep=' *', header=None, na_values='.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('smoking.shape = {}'.format(smoking.shape))\n",
    "smoking.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "smoking.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** something's up with the last row there; should remove it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic Gradient Descent\n",
    "\n",
    "So far what we've talked about has had a big effect on speed-up in the case where $d$ is large. When instead $n$ is very large we should look to Stochastic gradient descent. Note that Randomized coordinate descent is an example of stochastic gradient descent. Here, instead, we talk about \n",
    "\n",
    "* mini-batch stochast gradient descent (mini-batch SGD)\n",
    "* stochastic average gradient descent (SAG)\n",
    "* stochastic variance-reduced gradient methods (SVRG)\n",
    "\n",
    "Generally, SGD operates as follows. We seek to minimize an objective function \n",
    "$$\n",
    "f(w; X,y) = \\frac{1}{n}\\sum_{i=1}^n f_i(w)\n",
    "$$\n",
    "where, typically, $f_i$ is associated to the $i$-th observation. In usual gradient descent, we march *via* the updates\n",
    "$$\n",
    "w \\leftarrow w - \\alpha \\nabla_w f(w; X,y)\n",
    "$$\n",
    "whereas in SGD we approximate the gradient $\\nabla_w f(w; X,y)$ by $\\nabla_w f_i(w)$ so that the SGD algorithm is given by \n",
    "1. Let $w^0$ be an initial parameter vector; $(\\alpha_t)$ a sequence of learning rates\n",
    "2. Repeat until stopping condition satisfied:\n",
    "  1. Randomly shuffle $[n]$; \n",
    "  2. For each $i$ (in the new order):\n",
    "    1. $w^{t+1} \\leftarrow w^t - \\alpha_t \\nabla_w f_i(w)$\n",
    "\n",
    "3. Return $w^T$\n",
    "\n",
    "## Mini-batch SGD\n",
    "\n",
    "In mini-batch SGD we take *batches* of observations (say of size $m_t$ on iteration $t$) and approximate $\\nabla_w f(w; X,y)$ *via*\n",
    "$$\n",
    "\\nabla_w f(w; X,y) \\approx \\frac{1}{m_t} \\sum_{i=1}^{m_t} \\nabla_w f_i(w)\n",
    "$$\n",
    "\n",
    "**Exercise:** Using the `randomizedCoordinateDescent` function as a starting point, write your own function that performs mini-batch SGD.\n",
    "\n",
    "## Stochastic Average Gradient\n",
    "\n",
    "$w$ is computed according to the updates\n",
    "$$\n",
    "w^{t+1} \\leftarrow w^t - \\frac{\\alpha_t}{n} \\sum_{i=1}^n g_i^t\n",
    "$$\n",
    "where $g_i^t$ is defined in the following way. On each iteration $t$, an index $i_t \\in [n]$ is chosen randomly and we define\n",
    "$$\n",
    "g_i^t := \\begin{cases}\n",
    "\\nabla_i f(w^t) & i = i_t\\\\\n",
    "g_i^{t-1} & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "where $\\nabla_i$ corresponds to the $i$-th element of the gradient (*i.e.*, $\\partial/\\partial w_i$)\n",
    "\n",
    "### References\n",
    "1. [Original paper](https://arxiv.org/pdf/1202.6258.pdf)\n",
    "2. [Minimizing Finite Sums with SAG](https://arxiv.org/pdf/1309.2388.pdf)\n",
    "2. [SAG software](https://www.cs.ubc.ca/~schmidtm/Software/SAG.html)\n",
    "3. [SAG slides](https://www.cs.ubc.ca/~schmidtm/Documents/2014_Google_SAG.pdf)\n",
    "\n",
    "## SVRG\n",
    "\n",
    "We refer the interested reader to reference 3. below instead of repeating the derivation and intuition here.\n",
    "\n",
    "### References\n",
    "1. [Accelerating SGD using predictive variance reduction](https://papers.nips.cc/paper/4937-accelerating-stochastic-gradient-descent-using-predictive-variance-reduction.pdf)\n",
    "2. [Practical SVRG](https://arxiv.org/pdf/1511.01942.pdf)\n",
    "3. [SVRG intro from Stanford](http://cs.stanford.edu/~ppasupat/a9online/1321.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further reading\n",
    "\n",
    "We've barely scratched the surface. Some things for which there was not enough time to write notes...\n",
    "\n",
    "## Regularization\n",
    "\n",
    "**Idea:** avoid *overfitting* the model to the data by constraining the size/magnitude/norm/behaviour of the parameters\n",
    "\n",
    "### $L^2$ regularization\n",
    "\n",
    "*cf.* Ridge regression, elastic net, *etc.*\n",
    "\n",
    "### $L^1$ regularization\n",
    "\n",
    "Robust linear regression is an excellent way of demonstrating the difference between optimizing the $L^1$ norm and optimizing the $L^2$ norm. \n",
    "\n",
    "Also *cf.* Compressed Sensing; the $\\ell^1$ norm is *sparsity promoting*. That is, under certain conditions, the solution to the constrained $\\ell^1$ minimization problem \n",
    "$$\n",
    "\\hat x := \\argmin \\|x\\|_1 \\quad \\text{s.t.}\\quad \\|Ax - y\\|_2^2 \\leq \\eta\n",
    "$$\n",
    "is unique and $s$-sparse if elements of $A$ come from a certain distribution and $m \\geq C s \\log (N/s)$ for a known [and small] constant $C$.\n",
    "\n",
    "### Block sparsity\n",
    "\n",
    "This constrains the structure of the parameters - *e.g.* if some parameters are zero then it forces other parameters to be zero as well. \n",
    "\n",
    "## Constrained Optimization\n",
    "### Projected Gradient\n",
    "\n",
    "### Proximal Gradient\n",
    "\n",
    "### Newton's method\n",
    "\n",
    "From a second-order Taylor expansion, one sees \n",
    "$$\n",
    "f(w^*) = f(w^j) + \\nabla f(w^j) p^j + \\frac{1}{2} (p^j)^T \\nabla^2 f(w^j) p_j + O (\\|p^j\\|^3)\n",
    "$$\n",
    "for an unknown vector $p_j$. Since $\\nabla^2 f(w^j)$ is PSD, seek a particular vector that \"cancels\" this second-order term with the gradient term. Gets superlinear convergence. \n",
    "\n",
    "### References: \n",
    "1. Boyd & Vandenberghe's [Convex Optimization](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf)\n",
    "2. Bertsekas's Convex Optimization Algorithms\n",
    "3. [Ascher & Greif](http://gw2jh3xr2c.search.serialssolutions.com/?sid=sersol&SS_jc=TC0001261310&title=A%20first%20course%20in%20numerical%20methods) (pg. 261–265)\n",
    "\n",
    "## Dual methods\n",
    "### The Fenchel Dual, geometric multipliers and the KKT conditions\n",
    "### Dual coordinate ascent"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "toc": {
   "nav_menu": {
    "height": "336px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": true,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}