{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Name:  **Your name here**  \n",
    "UID:  **Your student ID num here**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Homework 6:  Gradient methods \n",
    "\n",
    "Put the file `grad2d.py` into your directory, along with `logreg.py`.  Then run the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from grad2d import grad2d, divergence2d\n",
    "from logreg import logistic_loss, logreg_objective, logreg_objective_grad, create_classification_problem\n",
    "import numpy as np\n",
    "from numpy import sqrt, sum, abs, max, maximum, logspace, exp, log, log10, zeros\n",
    "from numpy.linalg import norm\n",
    "from numpy.random import randn, rand, normal, randint\n",
    "import urllib, io\n",
    "import matplotlib.pyplot as plt\n",
    "np.random.seed(0)\n",
    "def good_job(path):\n",
    "    f = urllib.request.urlopen(path)\n",
    "    a = plt.imread(io.BytesIO(f.read()))\n",
    "    fig = plt.imshow(a)\n",
    "    fig.axes.get_xaxis().set_visible(False)\n",
    "    fig.axes.get_yaxis().set_visible(False)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 1:  Gradient descent\n",
    "Write a function that estimates the Lipschitz constant of a function $g$.  This can be done using the formula\n",
    "$$L \\approx \\frac{\\|g(x) - g(y)\\|}{\\|x-y\\|}.$$\n",
    "The inputs should be a function $g:\\mathbb{R}^n \\to \\mathbb{R}^m,$ and an initial vector $x$ in the domain of $g$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_lipschitz(g, x):\n",
    "    # Your work here\n",
    "    y = x+normal(size=x.shape)\n",
    "    L = norm(g(x)-g(y))/norm(x-y)\n",
    "    return L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now, run this unit test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Great - your Lipschitz estimator works!\n"
     ]
    }
   ],
   "source": [
    "g = lambda x: 10*x\n",
    "x = randn(3,4,5)\n",
    "L = estimate_lipschitz(g,x)\n",
    "assert abs(L-10)<1e-10, \"Your Lipschitz estimator is broken!\"\n",
    "print(\"Great - your Lipschitz estimator works!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write a routine that minimizes a function using gradient descent.\n",
    "The inputs $f$ and $grad$ are function handles.  The function $f: \\mathbb{R}^N\\to \\mathbb{R}$ is an arbitrary objective function, and  $grad: \\mathbb{R}^N \\to \\mathbb{R}^N$ is its gradient.  The method should minimize $f$ using gradient descent, and terminate when the gradient of $f$ is small.  I suggest stopping when\n",
    " $$\\|\\nabla f(x^k)\\|<\\|\\nabla f(x^0)\\|*tol$$\n",
    " where $x^0$ is an initial guess and $tol$ is a small tolerance parameter (a typical value would be $10^{-4}$).  \n",
    " \n",
    "  Use a backtracking line search to guarantee convergence.   The stepsize should be monotonically decreasing.  Each iteration should begin by trying the stepsize that was used on the previous iteration, and then backtrack until the Armijo condition holds:\n",
    "  $$f(x^{k+1}) \\le f(x^k) + \\alpha \\langle x^{k+1} - x^k, \\nabla f(x^k)\\rangle,$$\n",
    "  where $\\alpha \\in (0,1),$ and $\\alpha=0.1$ is suggested.\n",
    "\n",
    "  The function returns the solution vector $x_{sol}$, and also a vector $res$ containing the norm of the residual (i.e., the norm of the gradient) at each iteration.\n",
    "\n",
    "This initial stepsize should be $10L$, where $L$ is an estimate of the Lipschitz constant for the gradient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_descent(f, grad, x0, max_iters=10000, tol=1e-4):\n",
    "    # Your work here\n",
    "    t = 10/estimate_lipschitz(grad,x0)\n",
    "    res = zeros((max_iters,1),dtype=np.double)\n",
    "    x=x0\n",
    "    for iter in range(max_iters):\n",
    "        g = grad(x)\n",
    "        res[iter] = norm(g)\n",
    "        if res[iter]/res[0] < tol:\n",
    "            res = res[0:iter+1]\n",
    "            break\n",
    "        y = x - t*g\n",
    "        while f(y) > f(x) - 0.1*t*norm(g)**2:\n",
    "            t = t/2\n",
    "            y = x - t*g \n",
    "        x = y\n",
    "    return x, res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now run this unit test.  It will use your routine to fit a logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver terminated in 245 steps\n",
      "Terrific!  Your routine works!\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXzU9b3v8dcn+8oWwpKNgCCCCEmIqIiI1aq4EHDD7dTWtnQ/Pcfb06u391Y97Tm12p7WnvbYoqXaFndUQNwXREGRAAEEBCJrEiDsJITs3/vHTDDFZAhhwm9m8n4+Hnkk88vMbz75iXnnu/7MOYeIiEh7orwuQEREQpuCQkREAlJQiIhIQAoKEREJSEEhIiIBxXhdQFfo27evy83N9boMEZGwsnz58r3OufTjj0dkUOTm5lJcXOx1GSIiYcXMtrV1XF1PIiISUEQFhZlda2YzDx065HUpIiIRI6KCwjk33zk3o2fPnl6XIiISMSJyjEJEuoeGhgbKysqora31upSwkpCQQFZWFrGxsR16voJCRMJWWVkZqamp5ObmYmZelxMWnHPs27ePsrIyBg8e3KHXRFTXk4h0L7W1taSlpSkkToKZkZaWdlKtMAWFiIQ1hcTJO9lrpqBo5YUVZcxe2uY0YhGRbktB0cr8VRU8/fEOr8sQkW4sNzeXvXv3AjB+/PhOn+fxxx+noqIiKDUpKFpJjIumpr7R6zJEJMI0Nnbu98qSJUs6/Z7BDIqImvVkZtcC1w4dOrRTr0+MjeFofVNwixKRiPezn/2M2bNnk52dTd++fRk7diwvv/wy48ePZ/HixUyZMoUzzzyTn//859TX15OWlsbs2bPp378/+/bt45ZbbmHPnj2MGzeO1ncdTUlJobq6GoCHHnqIZ599lrq6OqZNm8b999/P1q1bmTx5MhMmTGDJkiVkZmYyd+5cFixYQHFxMbfddhuJiYl8+OGHJCYmdvrni6igcM7NB+YXFhZ+szOvT4qL5miDgkIkHN0/fy3rKg4H9ZwjM3pw77VnB3xOcXExc+bMYeXKlTQ2NlJQUMDYsWMBOHjwIO+99x4ABw4c4KOPPsLMeOyxx3jwwQf59a9/zf3338+ECRP46U9/yoIFC5g5c+YX3uONN95g06ZNfPzxxzjnmDJlCosWLSInJ4dNmzbx1FNP8eijj3LTTTcxZ84cbr/9dn7/+9/zq1/9isLCwlO+DhEVFKcqKS6aGrUoROQkfPDBBxQVFR37i/3aa6899r3p06cf+7qsrIzp06ezc+dO6uvrj61hWLRoES+88AIAV199Nb179/7Ce7zxxhu88cYb5OfnA1BdXc2mTZvIyclh8ODB5OXlATB27Fi2bt0a9J9RQdFKQmw0dY3NNDU7oqM05U4knJzoL/+u0rqr6HjJycnHvv7BD37AXXfdxZQpU1i4cCH33Xffse+daLqqc4577rmHb33rW/9wfOvWrcTHxx97HB0dzdGjR0/yJzgxDWa3khQXDUCtup9EpIMmTJjA/Pnzqa2tpbq6mgULFrT5vEOHDpGZmQnAE088cez4xIkTmT17NgCvvvoqBw4c+MJrr7jiCmbNmnVsvKK8vJzKysqAdaWmplJVVdWpn+l4alG00hIUNfVNJMfr0ojIiZ177rlMmTKFMWPGMGjQIAoLC2lrY9L77ruPG2+8kczMTM4//3y2bNkCwL333sstt9xCQUEBF198MTk5OV947eWXX8769eu54IILAN8g99///neio6PbreurX/0q3/72t4MymG2Bmk3hqrCw0HXmxkXPLy/jR8+tYtG/XUJOWlIXVCYiwbR+/XpGjBjhdRlUV1eTkpJCTU0NEydOZObMmRQUFHhdVkBtXTszW+6c+8Lot/5sbiUx1t+iaNBaChHpuBkzZrBu3Tpqa2u54447Qj4kTpaCopWWrietpRCRk/Hkk096XUKX0mB2K4kKCpGwE4nd513tZK+ZgqKVY11PCgqRsJCQkMC+ffsUFieh5X4UCQkJHX5NyHc9mVky8D9APbDQOTe7q97r2KwnTY8VCQtZWVmUlZWxZ88er0sJKy13uOsoT4LCzGYB1wCVzrlRrY5fCTwMRAOPOeceAK4DnnfOzTezZ4AuC4qWrqdatShEwkJsbGyH79ImnedV19PjwJWtD5hZNPAHYDIwErjFzEYCWUDL3t9d+hs8Kc6Xm9pBVkTkc54EhXNuEbD/uMPjgFLn3GbnXD3wNFAElOELC+jiej+fHqsWhYhIi1AazM7k85YD+AIiE3gBuN7MHgHmt/diM5thZsVmVtzZ/sqE2CjMNOtJRKS1UBrMbmtXLOecOwJ87UQvds7NBGaCb2V2pwowIzE2WkEhItJKKLUoyoDsVo+zgJO6PZOZXWtmMw8dOtTpIpLiotX1JCLSSigFxTJgmJkNNrM44GZg3smcwDk33zk3o60NuToqQS0KEZF/4ElQmNlTwIfAcDMrM7OvO+cage8DrwPrgWedc2tP8rzBaVFo1pOIyDGejFE4525p5/grwCuncN5TuhUqQHJ8jFZmi4i0EkpdTyEhJT6Gqlq1KEREWkRUUASj6yklPobqOgWFiEiLiAqKYAxmp8THUK0WhYjIMREVFMGQkhDDEbUoRESOUVAcJzU+hur6RpqbtW2xiAhEWFAEZYwiIQbntN+TiEiLiAqK4IxRxAJonEJExC+igiIYUhJ8S0uq6xo8rkREJDREVFAEo+spNd4XFFpLISLiE1FBEZSup2MtCgWFiAhEWFAEQ4q/RaExChERHwXFcVqCokotChERQEHxBakJalGIiLQWUUERjMHs5HiNUYiItBZRQRGMwezY6CgSYqMUFCIifhEVFMGSEh9LVa3WUYiIgIKiTT0SYzh8VC0KERFQULSpd1IcB2rqvS5DRCQkKCja0CsxloM16noSEYEIC4pgzHoC6JkUy6GjCgoREYiwoAjGrCdQ15OISGsRFRTB0isxlpr6JuoadU8KEREFRRt6JccBcEjjFCIiCoq29Er03bzooMYpREQUFG3pleQPCrUoREQUFG3pneTretKAtoiIgqJNPf1dTxqjEBFRULSppetJLQoRkQgLimAtuEuJjyEuJoq91XVBqkxEJHxFVFAEa8GdmTGkbzKlldVBqkxEJHxFVFAE0/ABqWzcraAQEVFQtGP4gFTKDx7VfSlEpNtTULRjeP9UADburvK4EhERbyko2nGmPyjW7VRQiEj3pqBoR1bvRAb0SGBJ6V6vSxER8ZSCoh1mxiVnpfP+pr00NDV7XY6IiGcUFAFMGt6P6rpGlm3Z73UpIiKeUVAEcNGwviTFRTNvVYXXpYiIeEZBEUBSXAxXjhrAgtU7qW3QTYxEpHtSUJzADWOzqKprVKtCRLqtkA8KMxtiZn82s+e9eP8LhqRx1oBUHl20meZm50UJIiKe6tKgMLNZZlZpZp8cd/xKM9tgZqVmdnegczjnNjvnvt6VdQZiZnz74jPYVFnNuxsqvSpDRMQzXd2ieBy4svUBM4sG/gBMBkYCt5jZSDM7x8xePu6jXxfX1yFXjx5IRs8E/rRos9eliIicdl0aFM65RcDxc0vHAaX+lkI98DRQ5Jxb45y75riPDv8Jb2YzzKzYzIr37NkTxJ8CYqOj+MZFQ/h4y37e2xjcc4uIhDovxigygR2tHpf5j7XJzNLM7I9Avpnd097znHMznXOFzrnC9PT04FXrd9v5OQzum8x989ZS16gZUCLSfXgRFNbGsXZHiZ1z+5xz33bOneGc+0XAEwfpxkVtiY+J5t5rR7Jl7xH+/MGWoJ9fRCRUeREUZUB2q8dZQFDmngbrxkXtmTS8H5eP7M/Db23i012Hu+Q9RERCjRdBsQwYZmaDzSwOuBmY50EdnfLzaaNITYjle7NXcKSu0etyRES6XFdPj30K+BAYbmZlZvZ151wj8H3gdWA98Kxzbm2Q3q/Lup5a9EtN4Hc357F57xHufmGN1laISMQz5yLvF11hYaErLi7u0vf4n4WlPPjaBr41cQj3XDWiS99LROR0MLPlzrnC44/HeFFMJPjOxWew61Atf1q0mZ5JsXx30lCvSxIR6RIRFRRmdi1w7dChXf9L28y499qzOVjTwIOvbeBIXSM/unw4Zm1N6hIRCV8hv9fTyejqWU/Hi44yfjM9j1vGZfOHdz/jvnlrNWYhIhEnoloUXoiOMv5z2jmkJsQyc9Fmdh6q5TfT80iO16UVkcgQUS2K0zHrqZ335Z7JZ/HTa0by1vrdXP/IEnbsrzmtNYiIdJWICorT3fXUmplx54TBPHHnOCoOHmXK7z/gnU93n/Y6RESCLaKCIhRcNCydud+fwICeidz5eDH/+cp66hubvS5LRKTTFBRdYHDfZF787nj+6fxBzFy0mRv/9KG6okQkbEVUUHg1RtGWhNhofjZ1FP9zWwGbK6uZ/PD7PLNsO5G4wFFEIlvAldlmVhDoxc65FUGvKAhOx8rsk7Fjfw0/em4VS7fs59Kz+vGL68+hX2qC12WJiPyD9lZmnygo3g1wTuec+1Iwigu2UAsKgOZmx1+WbOXB1z4lMS6an08dxTWjM7wuS0TkmE4FRbgKxaBoUVpZzf96toRVZYe44uz+/HvRKPr3UOtCRLx3ykFhZqPw3eP62G8159xfg1ZhEIVyUAA0NjXz6Ptb+O1bG4mLjuJ/Tz6LW8flEBWl7T9ExDvtBUWHBrPN7F7gv/0flwAPAlOCWmEQhNJgdiAx0VF8Z9IZvP4vEzknqyf/96VPmD7zQ0orq7wuTUTkCzo66+kG4FJgl3Pua8AYIL7LquokLxfcdUZu32Rmf+M8HrphNJv8M6N+8+ZG3ZNbREJKR4PiqHOuGWg0sx5AJTCk68rqPsyMGwuzeeuui7nqnIE8/PYmJv/2fd7buMfr0kREgI4HRbGZ9QIeBZYDK4CPu6yqbqhvSjwP35zPE3eOwwF3zPqYb/2tWAv1RMRzJz3rycxygR7OudVdUVAwhPpg9onUNTbx2Ptb+P07pTQ7x/cuGcqMiUNIiI32ujQRiWCnNOvJzCa2ddw5tygItQVduAdFi4qDR/mPV9azYPVOcvok8dNrRnLZyP5elyUiEepUg2J+q4cJwDhgeagtuGt1h7tvbtq0yetygmZx6V7unbeW0spqLhmezk+uHsHQfqlelyUiESaoC+7MLBt40Dl3SzCKC7ZIaVG01tDUzOOLt/K7tzdR09DEreNy+JfLhpGWEnKTz0QkTJ3SOoo2lAGjTq0kORmx0VF8c+IQFv7bJG4dl8OTH29n0kMLeWThZ9Q2aDqtiHSdjnY9/TfQ8sQoIA/Y6py7vQtr67RIbFEcr7Syil+88ilvf1pJZq9EfnzlcKaMycBMq7tFpHNOdYzijlYPG/GFxOIg1hdU3SEoWiwu3cvPF6xn/c7DjMnuxf+7egSFuX28LktEwpA2BYxgTc2OOSvK+NXrG6isquOyEf35tyuGM3yABrxFpOM6u834Gj7vcvoC59zo4JQXXN0tKFrU1Dfy5/e3MHPRZqrrG5mWl8m/fvlMsvskeV2aiISBzgbFIP+X3/N//pv/821AjXPu34NaZZB016BoceBIPX987zMeX7KVZue4dVwO3//SMNJTNUNKRNp3qmMUi51zF57omNcidR1FZ+06VMvv3tnEM8t2EBcdxZ0Tcpkx8Qx6JsZ6XZqIhKBTnR6bbGYTWp1sPJAcrOKCJdx2j+1qA3om8J/TzuGtuy7mspH9+cO7nzHxwXf5w7ulVNc1el2eiISJjrYoxgKzgJbfwAeBO3XP7PCytuIQv35jI+98WkmvpFi+edEQ7hifS0p8jNeliUgICMqsJ/8W4+acC+k7AykoAlu14yC/fWsj727YQ++kWL45cQhfuUCBIdLddXYw+3bn3N/N7K62vu+c+68g1hg0CoqOKdlxkIdbBcaMiWfwlQsGkazAEOmWOjtG0TIOkdrOh4SxvOxe/OVr43jxu+MZk92LX772KRN++Q6PLPyMIxrDEBE/LbiTY1ZsP8DDb23ivY176JUUy9fGD+aO8YPolRTndWkichqc0qwnM3vQzHqYWayZvW1me80sJPd5ks4ryOnNE3f6WhiFg/rwm7c2cuED7/CLV9ZTWVXrdXki4pGOznoqcc7lmdk0YCrwr8C7zrkxXV1gZ6hFERzrdx7mkYWf8fLqCmKio5hemM2MiUO00lskQp3qOoqWFVpXAU855/YHrTIJWSMG9uB3t+Tzzv+axPUFmTy9bDuTfrWQu54tobSyyuvyROQ06WiL4gF8LYmj+O5u1wt42Tl3XteW1zlqUXSNnYeO8tj7W3hy6XZqG5u4YuQAvjPpDMZk9/K6NBEJglNeR2FmvYHDzrkmM0sGUp1zu4JcZ1AoKLrWvuo6Hl+ylceXbKWqtpHzBvfhWxcPYdKZ/YiK0v0wRMLVqe71lATcBeQ452aY2TBguHPu5eCX+oX3ngpcDfQD/uCce+NEr1FQnB5VtQ08s2wHsz7YQsWhWob2S2HGRUMoys8gPiba6/JE5CSd6hjFX4B6YLz/cRnw8w686SwzqzSzT447fqWZbTCzUjO7O9A5nHMvOee+CXwVmN7BeuU0SE2I5RsXDeG9H1/Cb6fnERsdxY/nrGbCL337SR2qafC6RBEJgo62KIqdc4VmttI5l+8/tupEs57MbCJQDfzVOTfKfywa2Ah8GV/gLANuAaKBXxx3ijudc5X+1/0amN2R/aXUovCGc44PSvcyc9Fm3t+0l6S4aKafm83XJwwmq7dmSomEuvZaFB3dq6HezBLx38TIzM4A6k70IufcIjPLPe7wOKDUObfZf66ngSLn3C+Aa9oo3IAHgFcDhYSZzQBmAOTk5HTgR5JgMzMuGpbORcPSWVdxmMfe38zfPtzGXz/cxlXnDOSbFw1mdJYGvkXCzQm7nvy/qP8IvAZkm9ls4G3gx518z0xgR6vHZf5j7fkBcBlwg5l9u70nOedmOucKnXOF6enpnSxNgmVkRg/+a3oei358CXdemMu7n1Yy5feLueGRJSxYvZPGpmavSxSRDjphi8I558zsh8DlwPmAAT90zu3t5Hu2NS0m0O1Wfwf8rpPvJR7L6JXIT64eyQ8uHcZzxWU8sWQr33tyBRk9E/jK+FxuPjdbW4SIhLiOdj19BAxxzi0IwnuWAdmtHmcBFUE4b+s73AXjdBJEPRJi+fqEwXx1fC5vr9/NXxZv5YFXP+XhtzZxXUEmX7swl6H9tM+kSCjq6GD2OuBMYBtwBF+rwDnnRnfgtbn4Fue1DGbH4BvMvhQoxzeYfatzbm3nfoQv0mB2eFhXcZjHl2zhpZIK6hubmXhmOndemMvEYelajyHigVNdRzGorePOuW0neN1TwCSgL7AbuNc592czuwr4Lb6ZTrOcc/9xwiI6QPfMDk/7qut4cul2/vrRNvZU1TEkPZmvXTiY6wsySYrTvTFETpeg3OEuXKhFEZ7qG5t5Zc1OZi3ewuqyQ/RIiOGmwmxuP38QuX1D7hbtIhFHQSFhwznHiu0HmLV4K69/sovGZsek4el85YJBXHxmP6LVLSXSJbpFUKjrKfLsPlzLk0u389TH26msqiO7TyK3nzeImwqz6Z2s2VIiwdQtgqKFWhSRp6GpmdfX7uJvH25j6Zb9xMVEMWVMBl+5YJAW8YkEiYJCIsaGXVX89cOtvLiynJr6JsZk9+Ir5w/i6tEDSYjVZoQindUtgkJdT93L4doGXlhext8+2sZne47QJzmOmwqzue28HN2FT6QTukVQtFCLontxzvHhZ/t44sOtvLluNw649Kx+3Hpejga/RU7CqW4KKBKyzIzxQ/syfmhfKg4e5cml23l62Q7eWl9MZq9Epp+bzfRzs+nfI8HrUkXCkloUEpEampp5c91unly6nQ9K9xIdZVw2oh+3njeIi4b21cpvkTZ0i64njVFIW7buPcJTy7bzXHEZ+4/Uk90nkZvPzeGmwmzSU+O9Lk8kZHSLoGihFoW0pa6xidfX7ubJpdv4aPN+YqKMy8/uz23nDeKCIWlqZUi3p6AQaeWzPdU8tXQ7z68o42BNA7lpSdwyLocbxmaRlqJWhnRPCgqRNtQ2NPHaJ7uYvXQby7YeIC46iitGDeDWcTmcP6QPvvt2iXQPCgqRE9i4u4onl27nhRVlHK5tZEjfZG46N5vrCjLpl6oZUxL5ukVQaDBbguFofRML1uzk2WU7+HjrfqKjjEvP6sfN47KZOCydmOgT3kFYJCx1i6BooRaFBEtpZTXPFe9gzooy9lbXM6BHAjeMzeKmwmxy0rT6WyKLgkLkFDQ0NfP2+kqeWbad9zbuodnBhUPTmH5uDpeP7K89piQiKChEgqTi4FGeX17Gs8U7KDtwlF5JsUzNy+TmcdmcNaCH1+WJdJqCQiTImpsdSz7bx9PLtvPG2t3UNzUzJrsXN5+bzbVjMkiJ1w45El4UFCJd6MCRel5cWc4zy3awYXcVSXHRXH3OQG4el01BTm9Ns5Ww0C2CQrOexGvOOUp2HOTZ4h3MK6ngSH0TQ/ulcFNhFtPys7RliIS0bhEULdSikFBwpK6RBat38vSy7azYfpDoKOOS4encMDabL53Vj7gYTbOV0KKgEPFQaWU1zy8v44UVZVRW1dEnOY6peZncWJjFiIEaAJfQoKAQCQGNTc28v2kvzy3fwVvrKqlvamZUZg9uKMiiKC+T3slxXpco3ZiCQiTEHDhSz9yScp5bXsbaisPERUdx2ch+3Dg2m4uG9dUKcDntFBQiIWxdxWGeW76DuSUV7D9ST/8e8UzLz+LGwizOSE/xujzpJhQUImGgvrGZdz6t5PnlO3h3wx6amh0FOb24sTCba0YPJDUh1usSJYIpKETCTGVVLS+tLOe54jI2VVaTEBvF5FEDuXFsFufrRkvSBRQUImHKOceqskM8V7yDeasqqKptJLNXItePzeKGgixtTihB0y2CQgvuJNLVNjTx+tpdPL+8jA9K9+IcjMvtw/VjM5l8zkB6qGtKTkG3CIoWalFId1Bx8CgvrixnzooyNu85QnxMFFecPYDrCjKZMFSzpuTkKShEIlRL19Sc5WXMW1XBoaMN9EuNZ2p+JtcXZDF8QKrXJUqYUFCIdAN1jU28+2klzy8vZ+GGShqbHWdn9OD6giym5GXQN0V7TUn7FBQi3cy+6jrmrarghRXlrCk/REyUMWl4OtcVZHHpiH7Ex+hmS/KPFBQi3djG3VXMWVHGSyvL2X24jp6JsVwzeiDXj80iP7uXtkEXQEEhIkBTs2Nx6V7mrCjj9bW7qG1oZkjfZK4ryGRqfiZZvTXVtjtTUIjIP6iqbeDVNbuYs6KMpVv2A3DBkDSuK/BNtdUd+rofBYWItGvH/ppjU2237ashMTaaK0cN4PqCLC44I41orQLvFhQUInJCzjmWbzvAnBXlvLzatwp8YM8EpuZncl1+JsP6a6ptJFNQiMhJqW1o4q31u3lhRTnvbfRtUHh2Rg+m5WcyJS+DfqkJXpcoQRa2QWFmI4AfAn2Bt51zj5zoNQoKkeDaU1XH/FUVvFRSzuqyQ0QZTBiWzrT8DK44ewBJcRrPiASeBIWZzQKuASqdc6NaHb8SeBiIBh5zzj3QgXNFAY86575+oucqKES6TmllFS+trODFleWUHzxKUlw0V5w9gGn5mVw4tK/GM8KYV0ExEagG/toSFGYWDWwEvgyUAcuAW/CFxi+OO8WdzrlKM5sC3A383jn35IneV0Eh0vWamx3F2w7w4soyXl69k6raRtJT4ykak8HU/EzOzuih9RlhxrOuJzPLBV5uFRQXAPc5567wP74HwDl3fEi0da4FzrmrT/Q8BYXI6VXb4Ns65IWVvq1DGpocZ/ZPYWp+JlPzMsnoleh1idIB7QWFFx2LmcCOVo/LgPPae7KZTQKuA+KBVwI8bwYwAyAnJycYdYpIByXERjP5nIFMPmcgB47U8/Kanby0spwHX9vAQ69v4LzBfbguP4srzxmgrdDDkBctihuBK5xz3/A//idgnHPuB8F6T7UoRELDtn1H/OMZZWzdV0N8TBSXjezPtLxMLh6eTqy2Qg8podSiKAOyWz3OAiqCceJWNy4KxulE5BQNSkvmh5cN458vHUrJjoO8uLKc+asqWLB6J32S47hm9ECm5WeSp/2mQpoXLYoYfIPZlwLl+Aazb3XOrQ3We6pFIRK6GpqaeW/DHl4sKefNdbupb2xmcN9kpuZlMjU/g0FpyV6X2G15NevpKWASvjUQu4F7nXN/NrOrgN/im+k0yzn3H0F6P90KVSSMHK5t4LU1u3hhZRkfbfbtNzV2UG+m5mdyzTkD6Z0c53GF3UvYLrjrDLUoRMJP+cGjzC0p58UV5WyqrCY22rhkeD+m5WfyJd0/47RQUIhIWHDOsbbiMC+uLGduSQV7q+tITYjhqlEDKcrP4PzBaURpUV+X6BZBoa4nkcjS2NTMks/28eLKcl5fu4ua+iYG9EhgSl4GRXkZjByoRX3B1C2CooVaFCKR52h9E2+u383clb5NChubHcP6+Rb1TRmTQXYf3XTpVCkoRCRi7D9Sz4I1O5m7spzibQcA/yB4XgZXj86gjwbBO0VBISIRacf+GuatquCllb5B8JgoY+KZ6RTlZfDlkf21s+1J6BZBoTEKke7LOcf6nVXMLSln3qoKdh6qPbazbVFeBhOG9iVGK8ED6hZB0UItCpHurbnZsXTLfuaWlPPKmp0crm0kzb8SvCg/k3ytBG+TgkJEuqW6xibe/XQP81aV89b6SuobmxmUlkTRmAyK8jM5Iz3F6xJDhoJCRLq9w7UNvPbJLuaWlLPks304B+dk9qQoL4Nrx2TQv0f3vr1rtwgKjVGISEftPlzL/FUVzC2pYE257/auF5yRRlFeJleO6p7boXeLoGihFoWInIzSymrmlZTzUkkF2/fXEBcTxWUj+lGUl8mk4endZvsQBYWIyAk451i54yBzV5bz8uqd7DtST4+EGK4ePZApYzI5b3CfiN4+REEhInISGpqa+aB0L3NXlvPGut3U1DcxsGcCU8ZkUJSXyYiBqRE3c6pbBIXGKESkK9TUN/Lmut3MLalgkX/7kDP7p1CUl0lRXgZZvSNj+6S5CKIAAAjASURBVJBuERQt1KIQka6yr7qOV9bs5KWSCpb7tw8pHNSborwMrjpnIGkp8R5X2HkKChGRINuxv4a5Jb7t0DdVVhMdZUwY2peivAwuP3sAKfHhtX2IgkJEpIs45/h0VxXzVlUwr6SC8oNHiY+J4rKR/ZkyJiNsZk4pKEREToPmZseK7QeYW1LBgjU72e+fOTV51ECK8jI4b0ga0SE6c0pBISJymjU0NbO4dC/zSip4fe0ujtQ30S81nmtG+268NDqrZ0jNnOoWQaFZTyISqo7WN/HOp5XMLSln4YY91Dc1k5uWxJQxGUzJy2RoP+/3nOoWQdFCLQoRCWWHjjbw+ie7mLvq8z2nzs7oQVFeBteMziCjV6IndSkoRERCUOXhWl5evZO5qypYteMgAOMG9/FNtx01kN6n8W59CgoRkRC3de8R5q2qYG5JOZ/tOUJMlHHxmelMycvgshH9Se7i6bYKChGRMOGcY93Ow8wrqTh2t77E2Gi+PLI/RXkZXDQsnbiY4N+tT0EhIhKGmpsdxdsOMLeknAVrdnKwpoFeSbHHptuOyw3eRoUKChGRMFff2MwHpXuYV1JxbKPCAT0SuHbMQIryMjk7o8cpTbdVUIiIRJCa+kbeWl/JvJIK3ttYSUOTY0h6Mn+6fSzD+qd26pztBUV4bUQiIiIAJMXF+NZgjMngYE09r36yi9c+2dUlO9lGVItCC+5ERDqvvRZF8IfNPeScm++cm9GzZ0+vSxERiRgRFRQiIhJ8CgoREQlIQSEiIgEpKEREJCAFhYiIBKSgEBGRgBQUIiISUEQtuGthZnuAbZ18eV9gbxDLOV3Cse5wrBlU9+kWjnWHY80Ag5xz6ccfjMigOBVmVtzWysRQF451h2PNoLpPt3CsOxxrDkRdTyIiEpCCQkREAlJQfNFMrwvopHCsOxxrBtV9uoVj3eFYc7s0RiEiIgGpRSEiIgEpKEREJCAFRStmdqWZbTCzUjO72+t62mNmW81sjZmVmFmx/1gfM3vTzDb5P/cOgTpnmVmlmX3S6libdZrP7/zXfrWZFYRY3feZWbn/mpeY2VWtvnePv+4NZnaFRzVnm9m7ZrbezNaa2Q/9x0P6egeoO9Svd4KZfWxmq/x13+8/PtjMlvqv9zNmFuc/Hu9/XOr/fq4XdXeac04fvnGaaOAzYAgQB6wCRnpdVzu1bgX6HnfsQeBu/9d3A78MgTonAgXAJyeqE7gKeBUw4HxgaYjVfR/wozaeO9L/byUeGOz/NxTtQc0DgQL/16nARn9tIX29A9Qd6tfbgBT/17HAUv91fBa42X/8j8B3/F9/F/ij/+ubgWe8uN6d/VCL4nPjgFLn3GbnXD3wNFDkcU0nowh4wv/1E8BUD2sBwDm3CNh/3OH26iwC/up8PgJ6mdnA01PpP2qn7vYUAU875+qcc1uAUnz/lk4r59xO59wK/9dVwHogkxC/3gHqbk+oXG/nnKv2P4z1fzjgS8Dz/uPHX++W/w7PA5eamZ2mck+ZguJzmcCOVo/LCPwP1ksOeMPMlpvZDP+x/s65neD7nw/o51l1gbVXZzhc/+/7u2lmteraC7m6/d0a+fj+yg2b631c3RDi19vMos2sBKgE3sTXujnonGtso7Zjdfu/fwhIO70Vd56C4nNtpXuozh2+0DlXAEwGvmdmE70uKAhC/fo/ApwB5AE7gV/7j4dU3WaWAswB/sU5dzjQU9s4Fkp1h/z1ds41OefygCx8rZoRbT3N/zlk6u4MBcXnyoDsVo+zgAqPagnIOVfh/1wJvIjvH+nulq4D/+dK7yoMqL06Q/r6O+d2+38xNAOP8nl3R8jUbWax+H7ZznbOveA/HPLXu626w+F6t3DOHQQW4huj6GVmMf5vta7tWN3+7/ek492bnlNQfG4ZMMw/ayEO34DTPI9r+gIzSzaz1JavgcuBT/DVeof/aXcAc72p8ITaq3Me8BX/bJzzgUMtXSah4Lj++2n4rjn46r7ZP6tlMDAM+NiD+gz4M7DeOfdfrb4V0te7vbrD4Hqnm1kv/9eJwGX4xlfeBW7wP+34693y3+EG4B3nH9kOC16PpofSB76ZIBvx9TX+xOt62qlxCL5ZH6uAtS114uvvfBvY5P/cJwRqfQpft0EDvr+ovt5enfia5n/wX/s1QGGI1f03f12r8f1PP7DV83/ir3sDMNmjmifg68pYDZT4P64K9esdoO5Qv96jgZX++j4Bfuo/PgRfcJUCzwHx/uMJ/sel/u8P8erfd2c+tIWHiIgEpK4nEREJSEEhIiIBKShERCQgBYWIiASkoBARkYAUFCIBmNkS/+dcM7s1yOf+P229l0io0fRYkQ4ws0n4djO95iReE+2cawrw/WrnXEow6hPpSmpRiARgZi07hD4AXOS/N8K/+jeEe8jMlvk3rvuW//mT/PdXeBLfgjHM7CX/Bo5rWzZxNLMHgET/+Wa3fi//aumHzOwT8913ZHqrcy80s+fN7FMzmx1OO5BK+Io58VNEBN+9HI61KPy/8A855841s3hgsZm94X/uOGCU822DDXCnc26/f6uHZWY2xzl3t5l93/k2lTvedfg2wxsD9PW/ZpH/e/nA2fj2EFoMXAh8EPwfV+RzalGIdM7l+PZKKsG3LXYavn2HAD5uFRIA/2xmq4CP8G0MN4zAJgBPOd+meLuB94BzW527zPk2yysBcoPy04gEoBaFSOcY8APn3Ov/cNA3lnHkuMeXARc452rMbCG+fX9OdO721LX6ugn9PyyngVoUIh1The9WnS1eB77j3yIbMzvTv5vv8XoCB/whcRa+rahbNLS8/jiLgOn+cZB0fLdmPe07pIq00F8jIh2zGmj0dyE9DjyMr9tnhX9AeQ9t3372NeDbZrYa326nH7X63kxgtZmtcM7d1ur4i8AF+HYIdsCPnXO7/EEjctppeqyIiASkricREQlIQSEiIgEpKEREJCAFhYiIBKSgEBGRgBQUIiISkIJCREQC+v9JGJYrPA/G5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define a classification problem\n",
    "X, y = create_classification_problem(100, 10, cond_number=10)\n",
    "# Define the logistic loss function, and its gradient\n",
    "f = lambda w: logreg_objective(w,X,y)\n",
    "grad = lambda w: logreg_objective_grad(w,X,y)\n",
    "# Pick the initial guess\n",
    "w0 = zeros((10,1))\n",
    "\n",
    "# Now, solve the minimization problem\n",
    "w, res = grad_descent(f,grad,w0)\n",
    "\n",
    "# Check the solution\n",
    "assert norm(grad(w))/norm(grad(w0))<1e-4, \"ERROR:  your gradient descent routine did not the minimize the function\"\n",
    "print(\"Solver terminated in %d steps\"%len(res))\n",
    "print(\"Terrific!  Your routine works!\")\n",
    "\n",
    "_, res_g = grad_descent(f,grad,w0)\n",
    "n_g = list(range(len(res_g)))\n",
    "plt.semilogy(n_g,res_g)\n",
    "plt.legend(('gradient',))\n",
    "plt.xlabel('iteration')\n",
    "plt.ylabel('residual')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now, implement a gradient solver that begins each iterations using a Barzilai-Borwein stepsize \n",
    "  $$\\tau = \\frac{\\langle x^{k} - x^{k-1} ,x^{k} - x^{k-1}   \\rangle}{\\langle x^{k} - x^{k-1} ,\\nabla f(x^{k}) - \\nabla f(x^{k-1})   \\rangle}.$$\n",
    "  Your routine should perform only 1 gradient computation per iteration (i.e., your method should not require any \"extra\" work).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_descent_bb(f, grad, x0, max_iters=10000, tol=1e-4):\n",
    "    # Your work here!\n",
    "    t = 10/estimate_lipschitz(grad,x0)\n",
    "    res = zeros((max_iters,1),dtype=np.double)\n",
    "    x=x0\n",
    "    grad_x = grad(x)\n",
    "    for iter in range(max_iters):\n",
    "        res[iter] = norm(grad_x)\n",
    "        if res[iter]/res[0] < tol:\n",
    "            res = res[0:iter+1]\n",
    "            break\n",
    "        y = x - t*grad_x\n",
    "        while f(y) > f(x) - 0.1*t*norm(grad_x)**2:\n",
    "            t = t/2\n",
    "            y = x - t*grad_x \n",
    "        grad_y = grad(y)\n",
    "        delta_x = y-x\n",
    "        delta_g = grad_y-grad_x\n",
    "        t = norm(delta_x)**2/np.sum( delta_x.ravel()*delta_g.ravel() )\n",
    "        \n",
    "        x = y\n",
    "        grad_x = grad_y\n",
    "    return x, res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now run this unit test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver terminated in 42 steps\n",
      "Terrific!  Your BB routine works!\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3yV5fn48c+VvfdeJOwRwgrIEnCzcStVf462qK1+u1vtcHz1W63WtlZrLYqr4p5MRUVAFIWwlxiWEBL2Jjvcvz+e55AQsgjn5Ixc79crr+c59znnOVceQq7cW4wxKKWUUo3xc3cASimlPJsmCqWUUk3SRKGUUqpJmiiUUko1SROFUkqpJgW4OwBXSEhIMNnZ2e4OQymlvMry5cv3G2MS65f7ZKLIzs6moKDA3WEopZRXEZHvGyrXpiellFJN8qlEISITRGTqkSNH3B2KUkr5DJ9KFMaYmcaYKdHR0e4ORSmlfIZP9lEopdqHqqoqioqKKC8vd3coXiUkJISMjAwCAwNb9HpNFEopr1VUVERkZCTZ2dmIiLvD8QrGGA4cOEBRURE5OTkteo9PNT0ppdqX8vJy4uPjNUmcBREhPj7+rGphmiiUUl5Nk8TZO9t7pomijhUf/5evX33Q3WEopZRH0URRR/XGOXTc/JK7w1BKtWPZ2dns378fgKFDh7b6Oi+99BLFxcVOiUkTRR0nA8MJQUdPKKWcq7q6ulXv++qrr1r9mc5MFD416klEJgATOnfu3Kr3m8Bwwkw5GAPa7qmUaqGHHnqI6dOnk5mZSUJCAgMGDGDWrFkMHTqUL7/8kokTJ9K1a1cefvhhKisriY+PZ/r06SQnJ3PgwAEmT57Mvn37GDRoEHV3HY2IiOD48eMAPP7447z11ltUVFRwxRVX8OCDD7J9+3bGjBnD8OHD+eqrr0hPT+fDDz9k9uzZFBQUcMMNNxAaGsqSJUsIDQ1t9ffnU4nCGDMTmJmfn//jVl0gOJwAOUllRRlBIWHODU4p5VIPzlzPhuKjTr1mz7Qo7p/Qq8nXFBQU8O6777Jy5Uqqq6vp378/AwYMAODw4cMsXLgQgEOHDvH1118jIjz//PM89thjPPHEEzz44IMMHz6c++67j9mzZzN16tQzPmPevHkUFhaydOlSjDFMnDiRRYsWkZWVRWFhIa+//jrPPfcc1157Le+++y433ngjTz/9NH/961/Jz88/5/vgU4ninAVFAlB2/IgmCqVUiyxevJhJkyad+ot9woQJp5677rrrTp0XFRVx3XXXUVJSQmVl5ak5DIsWLeK9994DYNy4ccTGxp7xGfPmzWPevHn069cPgOPHj1NYWEhWVhY5OTn07dsXgAEDBrB9+3anf4+aKOrwCw4HoPzEUaITUt0cjVLqbDT3l7+r1G0qqi88PPzU+d13380vf/lLJk6cyIIFC3jggQdOPdfccFVjDPfeey+33377aeXbt28nODj41GN/f3/KysrO8jtonnZm1+EXYtUoKkqdW31VSvmu4cOHM3PmTMrLyzl+/DizZ89u8HVHjhwhPT0dgJdffvlU+YgRI5g+fToAc+fO5dChQ2e897LLLuOFF1441V+xa9cu9u7d22RckZGRHDt2rFXfU31ao6gjICQCgMoy59xcpZTvGzhwIBMnTqRPnz506NCB/Px8GlqY9IEHHuCaa64hPT2dwYMHs23bNgDuv/9+Jk+eTP/+/Rk5ciRZWVlnvPfSSy9l48aNDBkyBLA6uV999VX8/f0bjeuWW27hjjvucEpntjRVbfJW+fn5pjUbF61d8hG9P76ODRe/Qs/hk1wQmVLKmTZu3EiPHj3cHQbHjx8nIiKC0tJSRowYwdSpU+nfv7+7w2pSQ/dORJYbY87o/dYaRR1BoVbTU3W51iiUUi03ZcoUNmzYQHl5OTfffLPHJ4mzpYmijuAwO1GUHXdzJEopb/Laa6+5OwSX0s7sOoLDrXZFU6E1CqWUcvD4GoWIhAPPAJXAAmPMdFd9VohdozhZoTUKpZRycEuNQkReEJG9IrKuXvloEdkkIptF5B67+ErgHWPMj4GJrowrLMJKFFSecOXHKKWUV3FX09NLwOi6BSLiD/wLGAP0BCaLSE8gA9hpv6zGlUEFBQRwwgRrolBKqTrckiiMMYuAg/WKBwGbjTFbjTGVwBvAJKAIK1mAi+MVEcokBL8qTRRKqZbZvn07ubm5Z5TXXS7c23lSZ3Y6tTUHsBJEOvAecJWI/BuY2dibRWSKiBSISMG+fftaHUQZofhVlbb6/Uop5Ws8KVE0tNiJMcacMMbcaoy5s6mObGPMVGNMvjEmPzExsdVBlPuF4l+jiUIp1XLV1dXcfPPN5OXlcfXVV1Naav0Oefzxxxk0aBCDBg1i8+bNbo6y9Txp1FMRkFnncQZwVrtunOt+FAAVfqEEVmuiUMrrzL0Hdq917jVTesOYR5t92aZNm5g2bRrDhg3jtttu45lnngEgKiqKpUuX8sorr/Dzn/+cWbNmOTe+NuJJNYplQBcRyRGRIOB6YMbZXMAYM9MYM6WhdVZaqsovjMAa7aNQSrVcZmYmw4YNA+DGG29k8eLFAEyePPnUccmSJW6L71y5pUYhIq8Do4AEESkC7jfGTBORu4CPAX/gBWPM+rO87jnXKCoDwgmp2N3q9yul3KQFf/m7Sv1lwh2P65Y3t5S4J3PXqKfJxphUY0ygMSbDGDPNLp9jjOlqjOlkjPm/Vlz33GsUgZGEndQJd0qpltuxY8epGsPrr7/O8OHDAXjzzTdPHR0rv3ojT2p68gg1QZGEGe2jUEq1XI8ePXj55ZfJy8vj4MGD3HnnnQBUVFRw3nnn8eSTT/L3v//dzVG2nid1Zp8zZzQ9meBoQqiE6goICG7+DUqpdi07O5sNGzacUe7YkvT+++9v44icz6dqFM5oeiI4CoCq0sNOikoppbybTyUKZ/ALtZLMiSP1J44rpVT75FOJQkQmiMjUI0eOtPoa/mExAJQdPeCssJRSLuSLu3S62tneM59KFM5oegoMtxJF+fEzNzhXSnmWkJAQDhw4oMniLBhjOHDgACEhIS1+j091ZjtDUHgcAJUnNFEo5ekyMjIoKiriXNZ3a49CQkLIyMho/oU2n0oUzhj1FBoZC0DVidY3Xyml2kZgYCA5OTnuDsPnadNTPWFR8QCcLNNRT0opBT6WKJwhIjKaGiOcLNMahVJKgSaKM0SEBnKMMKRcE4VSSoEmijME+vtxnDD8Ko+6OxSllPIIPpUonDGPAuCERBCgiUIppQAfSxROWcIDOOIfS1ilb+x1q5RS58qnEoWzHA5MIqZqr7vDUEopj6CJogEnQlKIPnkYqsrdHYpSSrmdJooGVISlWidHd1nHY7vh5En3BaSUUm6kiaIB1ZFp1snRXbBrBTzRDVZNd29QSinlJrqER0PXibbWQKk8uJOgoq+swsPfn2N0SinlnXyqRuGsUU+BsVaiKN+/HTbOtAqDws8xOqWU8k4+lSicJSoyigMmEr/vF4Njhnal7qOtlGqfNFE0IC48iEKTQXjJ17WFVZoolFLtkyaKBsSFB7LmZEfE2COdAsOg8oR7g1JKKTfxqc5sZ4kJC2LtSXuN+5gsQLRGoZRqt7RG0YCY0EBWm07Wg6SeVke21iiUUu2UJooGBPj7cSQ4nX3BWZAzwmp60hqFUqqd8qlE4azVYwGSokK4N+0FGPJTCArTUU9KqXbLpxKFs+ZRAHRPjWJjyTHrQWA4VGnTk1KqffKpROFMPVIj2XW4jCNlVVqjUEq1a5ooGtEzNQqAjSVHtY9CKdWuaaJoxGmJQkc9KaXaMU0UjUiMDCY1OoRvth6snXBnjLvDUkqpNqeJohEiwoXdk1hUuI8q/1AwNVBT6e6wlFKqzWmiaMIlPZMpraxh2xG7JqHNT0qpdkgTRROGdIonOjSQghK7JqEd2kqpdkgTRROCA/y5ol86y3ZVWAU6RFYp1Q5pomjG9YMyOXYyyHqgk+6UUu2QxycKEekoItNE5B13fH73lCi6ZCQBUF56zB0hKKWUW7k0UYjICyKyV0TW1SsfLSKbRGSziNzT1DWMMVuNMT90ZZzNmTSwCwAL1mxzZxhKKeUWrq5RvASMrlsgIv7Av4AxQE9gsoj0FJHeIjKr3leSi+Nrke6drL0plqz7jvKqGjdHo5RSbculicIYswg4WK94ELDZrilUAm8Ak4wxa40x4+t97W3pZ4nIFBEpEJGCffv2OfG7ACKSAQir2M/T8ze37D2lB+Hrf+skPaWU13NHH0U6sLPO4yK7rEEiEi8izwL9ROTexl5njJlqjMk3xuQnJiY6L1qwlvAIimRIchXPLtzCpt0t6KuY9yf46B7YttC5sSilVBtzR6KQBsoa/bPbGHPAGHOHMaaTMeaRJi/sxP0ozhCZzHmJ1USFBnLPe2uorjnZ9Osds7iPFjs/FqWUakPuSBRFQGadxxmAU36bOnM/ijNEpBBcto8HJ/Zi5Y7DPDr326ZfH55gHU/sd34sSinVhtyRKJYBXUQkR0SCgOuBGc64sEtrFBFJcHw3E/qkccvQbJ5fvI0PVu5q/PXBkdaxVBOFUsq7uXp47OvAEqCbiBSJyA+NMdXAXcDHwEbgLWPMemd8nktrFJEpcGwPAH8Y14PzcuL47TtrWPRdIx3n1fZs7iNNJBOllPICrh71NNkYk2qMCTTGZBhjptnlc4wxXe1+h/9zZQxOE5FszcyuOEZg+SGm/qAPnZMimPLfApZsOXDm66vKrOORoraNUymlnMzjZ2afDdd2ZqdYxwWPwuMdiX5+EP+9OY+M2DBufnEpc9cUw+d/hn2brNc5FhDURKGU8nI+lShc2vQU18k6LptmHY/sJH7bLN66fQi5aVH87Y2ZsPAv8MGd1vPV5dbxWDHUVDs/HqWUaiM+lShcKn0ARGVAdRkMvRuSekLBNOLCg3jtx4O5Md3qq9h2uMaave1oejInoeKoGwNXSqlz41OJwqVNT35+0Psq67zLpdbX7rVwsoaQQH9uyrQmkX93NIBrnl1CWenx2vc6ahdKKeWFfCpRuLTpCWDo/8AlD0GHYRDfyZpUd8SaZO5XtAyAQUkn2X7gBBt27Kl9n6N2oZRSXsinEoXLhSfAsP8BP//aPosDm60NjfZZE/BizWFm3T2c6ICqU287dlybnpRS3ksTRWvFOxLFVti7weqLiM6CE/vpEB9Op5gAKv3DAPjVq0v4crNOvFNKeSefShQu7aOoLyIZgiLg4BarrwKg0wVWx3VVOVJdSlCEtYxHpH81Nzz/Db97Zw2HSytdH5tSSjmRTyUKl/dR1CUCcR3hgJ0ogqMhvb/13Il9UFUOYXEA/HlCZ24f2ZF3VhRx8d8WMnN1MUaXH1dKeQmfShRtLrkX7FoOJashJffUvhWc2Gt1YNuJItiUc++YHsy4axhpMaHc/fpKbntpGUWHSt0YvFJKtYwminPR6SIoOwi7CiBrMITbG/Kd2G/NzA61EoVj1FOvtGje/8kw/jS+J99sO8ilf1/Eswu3UFndzJLlSinlRj6VKNq0jwKg04Wc2l6j340QYW+YdHgHmJpTNQqqa4fH+vsJPxyew7xfjGBIx3genfsto/+xiAWbWryZn1JKtSmfShRt2kcBEB4POedDt3FWf0VkGvgFwt6N1vP1ahR1ZcSGMe2Wgbx4y0AMcMuLy/jRywXsOKDNUUopzxLQ1JMi0r+p540xK5wbjhe64V2rYxvAPwBiO1jDZaG2RtHEhLsLuicxtHM8LyzezlPzC7n47wu5fURH7hzVibCgJv95lFKqTTT3m+iJJp4zwIVOjMU7BQSd/jiuI+z42joPjrRqGM3MzA4O8OfOUZ24ol86j8zdyFPzN/N2QRG/urQrV/bPwN+vod1jlVKqbTSZKIwxF7RVID4jriMUzrPOA0MhMKzFS3ikRIfw5PX9uHFwBx6etYHfvLOGaYu38YdxPTi/S6ILg1ZKqca1uG1DRHKBnkCIo8wY84orgvJqjqU9AAJCITCkdm+KFhqYHcf7PxnGrLUlPPbRt9w0bSkjuiZy75ju9EiNcnLASinVtBZ1ZovI/cBT9tcFwGPARBfG1SptPuqpIXEda88DQ62vViwK6OcnTOyTxme/GskfxvZg1Y5DjP3nF/zm7dXsPqKr0Sql2k5LRz1dDVwE7DbG3Ar0AYJdFlUrtfmop4ZkDqo9D46wmp6qW796bHCAPz8e0ZFFv72AHw7L4cNVxYz66+c8Mncjh07ociBKKddraaIoM8acBKpFJArYC3Rs5j3tU0gU3FsE17wEqX0hIMQpy4zHhAXxx/E9+fSXIxndK4Wpi7Yy4rHP+cen33GsvKr5CyilVCu1NFEUiEgM8BywHFgBLHVZVN4uOBJ6XWEtR95UZ/aeDXDy7GZlZ8WH8Y/r+/HRz0YwtHM8//i0kBGPfc5/Fm6hrLLGCcErpdTpWpQojDE/McYcNsY8C1wC3Gw3QanmBIY23Jl9bA/8ewi8P6VVl+2WEsl/bspnxl3D6J0RwyNzv2Xk45/zypLtuiSIUsqpWtqZPcLxBWQBMfa5ak5gqLWSLFhrQJUetM5LD1jHtW9D5YlWXz4vI4ZXbhvEm1MGkx0fzn0frufCJxbwxtIdmjCUUk7R0uGxv6lzHgIMwmqC0gl3zalbo3j7FgiNgetehco6e2pvXwxdLzunjzmvYzxv3j6YRYX7+du8Tdzz3lqemr+ZO0Z14tr8DIID/M/p+kqp9qtFicIYM6HuYxHJxBoiq5pTd3jsoe+tWgWcnigqjlnHkzVQXQFBYa36KBFhZNdERnRJYOF3+/jnZ4X86YN1PD2/kDtGdmLyoCxCAjVhKKXOTmsXBSwCcp0ZiDN4xDyK+gLDrERwtARK98OxYqu8ok6icDQ9zbgb/px6zh8pIozqlsS7dw5l+o/OIzs+nAdnbmD4Xz7nuUVbKa2sPufPUEq1Hy2qUYjIU1hrO4GVXPoCq10VVGsZY2YCM/Pz83/s7lhO6XIJLJsGL0+wmqCqsBJD3X4JR9PUqunWsfIEBIWf80eLCMM6JzCscwJfbz3AU/ML+b85G/n3wi386PwcbhrcgciQwHP+HKWUb2tpH0VBnfNq4HVjzJcuiMf3dLoQhvwEFv+9tuxoyelNT3XPwdpK1QmJoq7BHeMZ3DGe5d8f4qn5hTz20Sb+vWALNw7uwK3DskmKDGn+IkqpdqmlfRQvuzoQnxadefrjo7vqJQq7RhEUYZUf3wex2S4JZUCHWF66dRBrig7z7MItPLtwC9MWb+Oq/hlMGdGRnATnJiillPdrbj+KtdQ2OZ3BGJPn9Ih80RmJotjqoxA/CIqsbXoKDLMSxYl9Lg8pLyOGZ24YwPb9J5j6xVbeWV7EG8t2MLpXCreP7ETfzBiXx6CU8g7N1SjG28ef2sf/2scbAN2KraWi009/fKzY7oeIsEY4OforAkOtYxskCofshHD+fEVvfnFxV176ahv/XfI9c9ftZnDHOG4f2YlRXRMR0f0wlGrPmtuP4nsAERlmjBlW56l7RORL4H9dGZzPiM6oPQ8IsWoU1eVWoggMrU0Ujn6JNkwUDomRwfzmsu7cOaozbyzdwfNfbOPWF5fRPSWSH5/fkQl90ggK8Kmdc5VSLdTS//nhIjLc8UBEhgLamN1SIdEQHGUlifjOcGRX7cimoPDapqeT9rBVNyQKh4jgAH50vrVa7V+v6UPNScOv3l7NsL/M56nPCjmoK9Yq1e60dNTTD4EXRMSxfvdh4DbXhOSjotKt+RQxWdbEu6g0axnygJDaGoWjU9uNicIhKMCPqwdkcFX/dBYV7ueFxdt44pPvePrzzVzRL53bhufQNTnS3WEqpdpAS0c9LQf62EuMizHGg2a0eYmk7tas7OhMa8mOkGir6ck/CMoPW6+pshOGByQKB8ds75FdEyncc4wXvtzOeyuKeGPZTs7vksBtw3MY2SURP93XWymf1dyopxuNMa+KyC/rlQNgjPmbC2PzLROetJboWDUdKo7C0SJI6gX+AVafBdTWKI57TqKoq0tyJI9c2ZvfXNaN15fu4OWvtnPri8volBjOrcNyuLJ/OmFBLd5dVynlJZrro3D0Q0Q28tUmRORyEXlORD4UkUvb6nOdKiQawuKspieAwzvsPooIq+npZA3UVFjPOVaW9VBx4UH89ILOLP7dhfzjur6EBQXwxw/WMeSR+TwydyM7D+qAOKV8SXOjnv5jHx9s7QeIyAtYw2z3GmNy65SPBp4E/IHnjTGPNhHHB8AHIhIL/BWY19p43K7unIrgCBB/q8nJ0U8hflaNwwsEBfhxeb90JvVNo+D7Q0z7YhvPLdrK1EVbuah7EjcNyeb8zgnaLKWUl2vpWk+PAQ8DZcBHWHtm/9wY82oL3v4S8DTwSp3r+QP/wtoEqQhYJiIzsJLGI/Xef5sxZq99/kf7fd7LUaMAqzYhYjU5OUY+RaZZzVLVFRDgcduSN0hEGJgdx8DsOIoPl/HaNzt4Y9kOPt24lJyEcG4c3IGrB2QQHarrSinljVo6PPZSY8xRrJpBEdCV0/eoaJQxZhFwsF7xIGCzMWarMaYSeAOYZIxZa4wZX+9rr1j+Asw1xqxo6HNEZIqIFIhIwb59ntnGD0BYPITGWuf+QRAYDtVltavJRqZYx3LvqFXUlxYTyq8v68aX91jNUrFhgTw0awOD//wZ9763lo0l3vl9KdWetTRROP4UHIu1IGD9X/xnKx3YWedxkV3WmLuBi4GrReSOhl5gjJlqjMk3xuQnJiaeY3guJALXvwYRKZDev3aSXam9T0WUvcx4uXcPLAsO8Ofyfum895NhzLp7OBP6pPLeiiLGPPkF1zz7FTNXF+sOfEp5iZYOUZkpIt9iNT39REQSgfJz+NyGGq2bWlPqn8A/m72oyARgQufOnc8htDbQYSj8epN1vux56+gYEhuZZh0rWpAo9qyHwzuh22jnx+hEuenRPHZ1H34/tgdvFezk1a93cPfrK0mMDOb6gZlcNzCTjNjWbdaklHK9FtUojDH3AEOAfGNMFdY6T5PO4XOLgLor5WUAxedwPcDaj8IYMyU6Orr5F3uKwHrLdpxqemokURzZVXv+1dMw6+eui83JYsKCmDKiEwt+PYoXbxlIr7Qonv58M+c/9jm3vriUTzbsobpGaxlKeZqWdmaHYS0MmAVMAdKAbsCsVn7uMqCLiOQAu4DrgR+08lrezbHtqWPuRJTdAtdQoti5DKZdDD9dBoldofKYV/Zl+PkJF3RP4oLuSew8WMpbBTt5c9lOfvxKASlRIVw7MJPrB2aSFhPq7lCVUrS8j+JFoBIYaj8uwhoF1SwReR1YAnQTkSIR+aExphq4C/gY2Ai8ZYxZf1aRN/xZnrcVanPiOlrHJfZgrlN9FA0kgMPfW0fHdqpVZdbQ2hrv3do0My6MX11qdX4/e+MAuqVE8tT8Qob/ZT4/fGkZn27YQ83JRlsllVJtoKV9FJ2MMdeJyGQAY0yZtHDtaWPM5EbK5wBzWvj5LeKRW6E2J6U3DLkLljxtPY6s05k993ew+TO46E/Qc5K1VhScuTZU5bHakVReKtDfj9G5KYzOTWHnwVLeWLaDtwqK+OyVAlKjQ7jO7stIjdZahlJtraU1ikoRCcXucBaRTkCFy6JqJa+sUQBcdH/teUSyNemu7BAUvAgHCmHLfOs5x0Q8R6JwzL1wJBAfkRkXxm8u685X91zIszf2p3NSBP/4tJBhj87nRy8XMP9brWUo1ZaarVHYNYdnsSbaZYrIdGAYcItrQzt7XlmjAAgIgl9+ay8WGGUtSb5n/ZlLejgSguPoo4nCwaplpDI6N5UdB2prGZ9u3ENKVAhXD8jgmvwMOsTrivdKuVKzicIYY0TkZ8ClwGCsoa0/M8bsd3Vw7UpUKuRdY52HREPRMutc/KDUnrZSXr9GUWYdfTRR1JUVH8ZvR3fn5xd35bONe3irYCfPLNjM059vZnDHOK4bmMnoXqmEBvm7O1SlfE5L+yi+BjoaY2a7Mphz5TXzKJoTEl3bcZ015Mwaxak+ihOnlxd+CondIKbeHt0+JCjAjzG9UxnTO5WSI2W8t2IXbxXs5Bdvrua+kPVM7JPGdQMz6Z0erVu4KuUkLU0UFwC3i8j3wAmsWoUxxuS5LLJW8Nqmp/pComuPCV3g20Lr8ak+Cnu5j1M1Crt8+lVWp/bvtrdZqO6UGh3KTy/ozJ0jO/HNtoO8XbCTd1cUMf2bHXRPieTa/Eyu6JdObHiQu0NVyqu1NFGMcWkU6nQ5I2D7FxCVYa0NVXoAjDk9UZw8aa0RBVaNwjFEtuyQe2J2Iz8/YUineIZ0iueBSb2YsaqYtwt28r+zNvDo3G+5pGcy1w7MZHjnBPx1JVulzlpLd7j73tWBqDpG/MZaByoiGbYtAlNjDZet20fhSBJgJQpHLaOdiwoJ5MbBHbhxcAe+3X2Ut5YV8f7KImavLSEt2tEBnklmnC4ZolRL+dR2ZD7TRyECnS+2zvfY8xBLD5zeR1FZZ3OgimO1/RXqlO4pUdw3oSe/G9ONTzfs5a2CnTz1+Wb+OX8zg3LiuLp/BmPzUokI9qn/Bko5nU/9D/GZPoq6wuKtY+nB05ueqjRRtFRwgD/j8lIZl5dK8eEy3l+5i3eWF/Hbd9dw/4z1jM5N4ar+GQzpFK9NU0o1wKcShU8Ki7OOdWsUFfUTxdHapidp6RzK9iktxuoA/8moTqzYcZh3VxQxc3Ux76/cRWp0CFf0S+eqARl0Soxwd6hKeQxNFJ7OUaM4VgLV9srulScar1H46wiflhARBnSIZUCHWO4b35NPN+7h3eVFPLtwC88s2ELfzBiuGpDBxLw0osN0Zz7VvvlUovCZPoq6HIlif2FtWf0+ivI6NQo//aV2tkIC/Rmfl8b4vDT2Hivnw5XFvLO8iD99sI6HZm7g4p5JXNU/g5FdEwnw1xqban98KlH4ZB9FcCQk9oCNM+3H0XYfhT3qKSTabnpy1Cg0UZyLpMgQfjyiIz86P4f1xUd5Z3kRM1YXM2ftbhIigri8r9U01SM1yt2hKtVmfCpR+Kxuo2Hx363z6AzYv8laXhwgOtPqv3DUKDRROIWIkK8HW8MAABeSSURBVJseTW56NL8f24MFm/by7ooiXl6ynecXb6NnapTVNNUnjcTIYHeHq5RLaT3aG3S15zuGxkL3sXCyGsoOW2UxHaxNjxw1Cm16crqgAD8u7ZXCf27K55vfX8wDE3ri7yc8NGsDgx/5jJtfWMoHK3dRWum9+4Io1RStUXiDzEFw9YvWjO1171pljq1TY7Ot2sXxvdZjf/0ndaW48CBuGZbDLcNyKNxzjPdX7uLDVcX8/M1VhAX5M7pXCpf3S2dop3jtz1A+Q3+reAMRyL3SOg+228Z3r7GOsdnW8dA262h0z+m20iU5kt+O7s6vL+3Gsu0H+WDVLmatKeG9lbtIjAxmYp80ruiXTq+0KF2gUHk1n0oUPjnqqb6ul0F4Um3ndkyWdTy41TpWV7onrnbMz084r2M853WM5/4JvViwaS/vr9zFK0u2M23xNjonRXBFv3Qm9U0jI1aXDlHeR4zxvZ3C8vPzTUFBgbvDcJ3CT2D61VbCuOEtmDoKgiJrt0T93XY4vAPeuhlueBvCE9wdcbt0uLSSOWt388HKXSzdbu0pMignjiv6pTM2N1XnZyiPIyLLjTH5Z5RrovBSlaXWENnqcvh7z9rywHD4QzF8OwfemAy3fgQdhrgvTgXAzoOlzFhdzHsritiy7wRB/n5c2D2Jy/ulc0H3RIIDdMMl5X6NJQqfanpqV4LCrK/6TU019mPHLO4aj9vavF3KjAs7tXTIul1HeX/lLmasLuaj9buJCglgXF4aV/ZPZ0BWLH663pTyMJoovF1AEITEQLk9XPZklb1XhZ0gqjVReBIRoXdGNL0zovn92O58ueUAH6zcxQcrd/H60h1kxIYysU8ak/qm0y0l0t3hKgVoovANEcm1iQKsWoVjvwpHzUJ5nAB/P0Z2TWRk10QevryaTzbs4f2Vu/jPoq08s2AL3ZIjmdg3jYl90nT/DOVWmih8Qf+bYN4fax/XVNapUegoKG8QHhzA5f3SubxfOgeOVzBnbQkfrirm8Y838fjHm+ifFcPEPmmMy9OZ4KrtaaLwBUPvhrzrYNVr8On9dqKwaxJao/A68RHB3DQkm5uGZFN0qJSZq0v4cNUuHpi5gf+dtYFhnROY2CeNy3JTiArRkVPK9XwqUbSLeRSNiUiC0BjrvLqitkahndleLSM2jDtHdeLOUZ34bs8xZqwqZsbqYn7zzhr+8ME6LuyWxKS+aVzQPYmQQB05pVzDpxKFT64eezb87SaJmoo6NQpNFL6ia3Ikv76sG7+6tCurdh5mxupiZq4u4aP1u4kIDuCyXilM7JvGMF0+RDmZTyWKdi/A3rSopkpHPfkwEaFfViz9smL547ieLNlygBmrdzF33W7eXVFEQkQQY3unMqlvGv2zYnX5EHXONFH4EkeNolprFO2Fv58wvEsCw7sk8NDluSzYtI8Zq4p5c9lOXlnyPekxoUzsm8akvml0T9E9NFTraKLwJY5tUGsqoUon3LU3wQH+XNYrhct6pXC8opp563czY3UxUxdt5d8LttA1OYJJfdOZkJdGVrwOt1Utp4nClwTUSRRao2jXIoIDuLJ/Blf2zzg13HbG6trhtn0yY5iQl8q4vFRSo0PdHa7ycJoofMlpTU+OPgo7YaycDruWw/i/uSc25Tb1h9vOWlPCrDXFPDx7Iw/P3sjA7Fgm9EljTG6qztFQDdJE4UsarFHYE+62fg7bFmmiaOcyYsO4Y2Qn7hjZiW37TzBrdTGz1pRw34freWDGeoZ0imd8Xhqje6UQGx7k7nCVh9BE4UscfRQN1SiqyrQZSp0mJyGcuy/qwt0XdeG7PceYtbqYmWtKuPe9tfzpg3UM75LAhLw0LumVrBP72jlNFL7k1DyKyjNXj60qrV1ZVql6uiZH8stLu/GLS7qyvvgoM9cUM2t1Cb96ezVB7/kxslsiE/qkcXGPJMKC9NdGe6P/4r7ktKanevMoKku1RqGaJSLkpkeTmx7NPaO7s3LnYWatLmH22mI+2bCHkEA/LuqRzIS8VEZ109ng7YXHJwoR6QH8DEgAPjPG/NvNIXmu0zqzy2rPwapRmBo4WQN++p9bNU9E6J8VS/+sWP44rgfLth9k5ppi5q7dzew1JUQEB3BJz2Qm9ElleOdEggJ0NrivcmmiEJEXgPHAXmNMbp3y0cCTgD/wvDHm0cauYYzZCNwhIn7Ac66M1+v5N1GjqKqTOIJ0DL06O3X3BX9gQi+WbD3ArNUlzF1XwvsrdxEdGsjoXimM75PKkI66hIivcXWN4iXgaeAVR4GI+AP/Ai4BioBlIjIDK2k8Uu/9txlj9orIROAe+1qqMY6mp70baxNDTb1EUVMBaKJQrRfg78f5XRI5v0siD12ey+LN+5i5uoTZa0t4s2An8eFBjOmdwvi8NAZlx+mOfT7ApYnCGLNIRLLrFQ8CNhtjtgKIyBvAJGPMI1i1j4auMwOYISKzgdcaeo2ITAGmAGRlZTklfq8TFAGdL4HlL9aW1W16qvtYKScICvDjwu7JXNg9mfKqGhZs2sfMNcW8s7yIV7/eQXJUMGN7pzKhTxr9MmN03Skv5Y4+inRgZ53HRcB5jb1YREYBVwLBwJzGXmeMmQpMBcjPzzfOCNTriMD1r8HDSYB9CxpqelLKBUIC/Rmdm8Lo3BROVFTz2bd7mbW6mOlf7+DFL7eTHhPK+LxUxvZOJS8jWpOGF3FHomjop6PRX+zGmAXAghZduD3vR+EQEAThCXBin/W4usLeQ9vR9KRDZJXrhQcHMLGPtY3r0fIqPlm/h5lripm2eBv/WbSVjNhQxuWlMr53GrnpUZo0PJw7EkURkFnncQZQ7IwLt/v9KBzCE+skivLTd7nTGoVqY1EhgVw1IIOrBmRwuLSSeRv2MHtNCdO+2MZ/Fm4lKy6Msb1TGZ+XSq80TRqeyB2JYhnQRURygF3A9cAP3BCH7wqLt47iZ68kW1b7nK4mq9woJiyIa/MzuTY/00oa6/cwa20Jz32xlWcXbqFDfBjjelvNU5o0PIerh8e+DowCEkSkCLjfGDNNRO4CPsYa6fSCMWa9kz5Pm54AwuKsY3CU1Ynt6MiG2rWflHKzmLAgrh2YybUDMzl0opJ5G3Yza00J/1m0lWcWbCE7Poxxdp9Gz1RNGu7k6lFPkxspn0MTHdPn8Hna9AQQlmAdRewaRZ1EoTUK5YFiw4O4bmAW1w3M4uCJSj5eb03qe3bhVv71+RZyEsIZ19taFr17SqQmjTbm8TOzz4bWKGzhdqKoOGYdyw7XPqd9FMrDxYUHMXlQFpMHZXHgeAUfr9/D7LXFPLNgM09/vpmOibVJo1uyJo224FOJQmsUtsRu1jE8CY4VQ9mh2uc0USgvEh8RzA/Oy+IH52Wx/3jFqZrGvz7fzFPzN9MpMZxxeWmMz0ula3Kku8P1WT6VKJSt15XWch5HdsFHvzs9UejwWOWlEiKCueG8DtxwXgf2Havgo/W7mb2mmKfmF/LPzwrpnBRxqqahScO5NFH4IhHoMQFW2CunaI1C+ZjEyGBuGtyBmwZ3YO+xcj5eZ3WE/3N+IU9+VkiXpAhrnkZeKp2TNGmcK59KFNpHUU9kmnXcu6G2TDuzlY9Jigw5tdXr3qPlfLTeShpPflbIPz4tpFtyJGPtmkbnpAh3h+uVxBjfW+0iPz/fFBQUuDsM9ys9CI/lQEwHOPy9VTb6LzD4DvfGpVQb2HO0nLlrS5izdjfLvj+IMdA9JdKap5GXSqdETRr1ichyY0x+/XKfqlGoesLiIL4LHCisLdMahWonkqNCuGVYDrcMy2H3kXLmrith9poSnvjkO5745DtNGmfBpxKFNj01IGPg6YlC+yhUO5QSHcKtw3K4dVgOJUfKrM2X1tYmDUfz1NjeKXTRjvAzaNOTr9swA966yToXPxj+S7joT+6NSSkPUXKkjI/W7WbO2hIKvj+EMdA5KeJU0mhv8zQaa3rSRNEeLH8ZdhXAmrdh0I/g0ofdHZFSHmfP0fJT8zSWbrf6NDomhjM211pGpEeq7ycNTRQKHs2CvOth7GPujkQpj7b3WDnz1u9hztoSvt56gJMGsuPDGNM7lXE+vGChdmYr8A/WzmylWiApMoQbB3fgxsEdOHC8gnkbrKQxddFW/r1gC5lxoadqGu1hEyafShTamd2MgGBdPVapsxQfEXxq7SnHKrdz1u4+tQlTekwoY3JTGJuXSt+MGJ/cI1ybntqTpwZAah+4+gV3R6KU1ztcWsknG/Ywd91uvijcR1WNITU6hNG5KYzrnUr/rFivSxra9KSspicdHquUU8SEBXFNfibX5GdypKyKzzbuYc7a3af2CE+OCmZMbipjclPIz47D38uSRl2aKNqTgCBNFEq5QHRoIFf2z+DK/hkcK69i/rd7mbO2hNeX7uClr7aTGBnM6F4pjOmdwqDsOAL8/dwd8lnRRNGeaGe2Ui4XGRLIpL7pTOqbzvGKaj7/di9z15Xw9vKd/Pfr74kPD+Ky3BTG5qYyuKN3JA2fShTamd2MgCDtzFaqDUUEBzChTxoT+qRRWlnNgk37mLO2hA9W7uK1b3YQGxbIZb1SGNM7laGd4gn00KShndntyfRr4VgJ3PGFuyNRql0rq6xh4Xf7mLuuhE837OFEZQ3RoYFc2jOZsb1TGdY5gaCAtk8a2pmtIL4TbFsENdXgr//0SrlLaJA/o3NTGJ2bQnlVDV8U7mfu2hI+Wrebt5cXERkSwCU9kxnXO5XhXRIIDvB3a7z626I9Se0L1WWw/ztI7unuaJRSQEigP5f0TOaSnslUVNfw5eb9zF6zm0827Oa9FbuIDA7goh5JjO2dyoiuiYQEtn3S0ETRnqT1tY4lqzRRKOWBggP8ubB7Mhd2T6ayujdfbdnPnLUlzNuwhw9WFRMe5M+FPZIZm5vCqG5JhAa1TdLQRNGexHeGwHAoXgV9f+DuaJRSTQgK8GNUtyRGdUvi/2pOsmTLAeauK+Hj9XuYubqY0EB/LuieyJjcVC7onkREsOt+nWuiaE/8/CGpO+zf5O5IlFJnIdDfjxFdExnRNZGHJp1k6baDzFlXwkfrrEl+QQF+jOyayNjeKVzcI5nIkECnfr4mivYmLB6O73V3FEqpVgrw92No5wSGdk7gwYm5LP/+EHPsjvBPNuzhtR+dx9DOCc79TKdeTXm+kBjYpzUKpXyBv58wKCeOQTlx3De+J6uKDpOXHu30z/HM2R2tJCITRGTqkSNH3B2K5wqNhfLD7o5CKeVkfn5C/6xYl8z09qlEYYyZaYyZEh3t/IzqM0JjoPwonDzp7kiUUl7CpxKFaoGQGMBAhda6lFIto4mivQmNsY5lh9wbh1LKa2iiaG9CY61jmfZTKKVaRhNFexNi1yi0Q1sp1UKaKNobbXpSSp0lTRTtjaNGoU1PSqkW0kTR3oRq05NS6uxoomhvAkMhIERrFEqpFvOKRCEi4SKyXETGuzsWnxASo30USqkWc2miEJEXRGSviKyrVz5aRDaJyGYRuacFl/od8JZromyHekyAlDx3R6GU8hKuXhTwJeBp4BVHgYj4A/8CLgGKgGUiMgPwBx6p9/7bgDxgAxDi4ljbj3F/dXcESikv4tJEYYxZJCLZ9YoHAZuNMVsBROQNYJIx5hHgjKYlEbkACAd6AmUiMscYowsVKaVUG3HHMuPpwM46j4uA8xp7sTHmDwAicguwv7EkISJTgCkAWVlZzopVKaXaPXd0ZksDZaa5NxljXjLGzGri+anGmHxjTH5iYuI5BaiUUqqWOxJFEZBZ53EGUOyMC+t+FEop5XzuSBTLgC4ikiMiQcD1wAxnXFj3o1BKKedz9fDY14ElQDcRKRKRHxpjqoG7gI+BjcBbxpj1Tvo8rVEopZSTiTHNdg94nfz8fFNQUODuMJRSyquIyHJjTH79cq+Yma2UUsp93DE81mVEZAIwATgqIoWtvEwCsN95UbUZjbtteWPc3hgzaNxtqUNDhT7Z9HQuRKSgoaqXp9O425Y3xu2NMYPG7Qm06UkppVSTNFEopZRqkiaKM011dwCtpHG3LW+M2xtjBo3b7bSPQimlVJO0RqGUUqpJmiiUUko1SRNFHa3Yec8tRGS7iKwVkVUiUmCXxYnIJyJSaB9jPSDOM3Y4bCxOsfzTvvdrRKS/h8X9gIjssu/5KhEZW+e5e+24N4nIZe6JGkQkU0Q+F5GNIrJeRH5ml3v0PW8ibo++5yISIiJLRWS1HfeDdnmOiHxj3+837TXtEJFg+/Fm+/lsd8TdKsYY/bL6afyBLUBHIAhYDfR0d1yNxLodSKhX9hhwj31+D/AXD4hzBNAfWNdcnMBYYC7WMvSDgW88LO4HgF838Nqe9s9KMJBj/wz5uynuVKC/fR4JfGfH59H3vIm4Pfqe2/ctwj4PBL6x7+NbwPV2+bPAnfb5T4Bn7fPrgTfdcb9b86U1ilqndt4zxlQCbwCT3BzT2ZgEvGyfvwxc7sZYAGuHQ+BgveLG4pwEvGIsXwMxIpLaNpGerpG4GzMJeMMYU2GM2QZsxvpZanPGmBJjzAr7/BjWopvpePg9byLuxnjEPbfv23H7YaD9ZYALgXfs8vr32/Hv8A5wkYg0tD+Px9FEUauhnfea+mF1JwPME5Hl9s5+AMnGmBKw/uMBSW6LrmmNxekN9/8uu4nmhTpNex4Zt92s0Q/rr1yvuef14gYPv+ci4i8iq4C9wCdYtZvDxlolu35sp+K2nz8CxLdtxK2jiaJWq3bec5Nhxpj+wBjgpyIywt0BOYGn3/9/A52AvkAJ8IRd7nFxi0gE8C7wc2PM0aZe2kCZ22JvIG6Pv+fGmBpjTF+sDdgGAT0aepl99Ji4z5Ymilou23nP2YwxxfZxL/A+1g/oHkezgX3c674Im9RYnB59/40xe+xfCieB56ht6vCouEUkEOuX7XRjzHt2scff84bi9pZ7DmCMOQwswOqjiBERx4KrdWM7Fbf9fDQtb+J0K00UtVy2854ziUi4iEQ6zoFLgXVYsd5sv+xm4EP3RNisxuKcAfw/eyTOYOCIo7nEE9Rru78C656DFff19oiWHKALsLSt4wNrFBMwDdhojPlbnac8+p43Fren33MRSRSRGPs8FLgYq3/lc+Bq+2X177fj3+FqYL6xe7Y9nrt70z3pC2sUyHdY7Yx/cHc8jcTYEWvEx2pgvSNOrLbOz4BC+xjnAbG+jtVkUIX119QPG4sTq1r+L/verwXyPSzu/9pxrcH6D59a5/V/sOPeBIxxY9zDsZoy1gCr7K+xnn7Pm4jbo+85kAestONbB9xnl3fESlybgbeBYLs8xH682X6+o7t+Vs72S5fwUEop1SRtelJKKdUkTRRKKaWapIlCKaVUkzRRKKWUapImCqWUUk3SRKFUE0TkK/uYLSI/cPK1f9/QZynlaXR4rFItICKjsFYyHX8W7/E3xtQ08fxxY0yEM+JTypW0RqFUE0TEsTroo8D59r4Iv7AXg3tcRJbZi9bdbr9+lL23wmtYk8UQkQ/sBRzXOxZxFJFHgVD7etPrfpY9U/pxEVkn1r4j19W59gIReUdEvhWR6d6y+qjybgHNv0QphbWPw6kahf0L/4gxZqCIBANfisg8+7WDgFxjLYENcJsx5qC9zMMyEXnXGHOPiNxlrAXl6rsSayG8PkCC/Z5F9nP9gF5Y6wd9CQwDFjv/21WqltYolGqdS7HWSVqFtSR2PNaaQwBL6yQJgP8RkdXA11iLwnWhacOB1421IN4eYCEwsM61i4y1UN4qINsp341STdAahVKtI8DdxpiPTyu0+jJO1Ht8MTDEGFMqIguw1vxp7tqNqahzXoP+H1ZtQGsUSrXMMaxtOh0+Bu60l8dGRLraq/nWFw0cspNEd6xlqB2qHO+vZxFwnd0Pkoi1NatbVqRVCvSvEaVaag1QbTchvQQ8idXss8LuUN5Hw9vPfgTcISJrsFY6/brOc1OBNSKywhhzQ53y94EhWCsEG+C3xpjddqJRqs3p8FillFJN0qYnpZRSTdJEoZRSqkmaKJRSSjVJE4VSSqkmaaJQSinVJE0USimlmqSJQimlVJP+P2Rl5Nj9jrrYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Minimize that logistic loss again\n",
    "f = lambda w: logreg_objective(w,X,y)\n",
    "grad = lambda w: logreg_objective_grad(w,X,y)\n",
    "\n",
    "w, res = grad_descent_bb(f,grad,w0)\n",
    "assert res[-1]/res[0]<1e-4, \"ERROR:  Your BB routine did not the minimize the function\"\n",
    "print(\"Solver terminated in %d steps\"%len(res))\n",
    "print(\"Terrific!  Your BB routine works!\")\n",
    "\n",
    "_, res_g = grad_descent(f,grad,w0)\n",
    "_, res_b = grad_descent_bb(f,grad,w0)\n",
    "n_g = list(range(len(res_g)))\n",
    "n_b = list(range(len(res_b)))\n",
    "\n",
    "plt.semilogy(n_g,res_g,n_b,res_b)\n",
    "plt.legend(('gradient','bb'))\n",
    "plt.xlabel('iteration')\n",
    "plt.ylabel('residual')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write a routine that uses Nesterov's accelerated gradient method\n",
    "\\begin{align}\n",
    "x^{k+1} &= y^k - \\tau \\nabla f(y^k)\\\\\n",
    "\\delta^{k+1} &= \\frac{1+\\sqrt{1+4(\\delta^k)^2}}{2}\\\\\n",
    "y^{k+1} &= x^{k+1}+\\frac{\\delta^k-1}{\\delta^{k+1}}(x^{k+1}-x^{k})\n",
    "\\end{align}\n",
    "The stepsize restriction for Nesterov's methods is $\\tau<1/L,$ however when $L$ is not known exactly you can use the line search condition\n",
    " $$f(x^{k+1}) \\le f(y^{k}) + \\alpha (x^{k+1}-y^k)^T\\nabla f(y^k), $$\n",
    " where $\\alpha \\in [1/2,1).$  I suggest choosing $\\alpha=1/2.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_descent_nesterov(f, grad, x0, max_iters=10000, tol=1e-4):\n",
    "    # Your work here!\n",
    "    t = 10/estimate_lipschitz(grad,x0)\n",
    "    res = zeros((max_iters,1),dtype=np.double)\n",
    "    y=x0\n",
    "    x=x0\n",
    "    delta = 1\n",
    "    for iter in range(max_iters):\n",
    "        g = grad(y)\n",
    "        res[iter] = norm(g)\n",
    "        if res[iter]/res[0] < tol:\n",
    "            res = res[0:iter+1]\n",
    "            break\n",
    "        xp1 = y - t*g\n",
    "        while f(xp1) > f(y) + 0.5*np.sum((xp1-y)*g):\n",
    "            t = t/2\n",
    "            xp1 = y - t*g\n",
    "        deltap1 = (1+np.sqrt(1+4*delta*delta))/2\n",
    "        yp1 = xp1 + (delta-1)/deltap1*(xp1-x)\n",
    "        \n",
    "        x = xp1\n",
    "        y = yp1\n",
    "        delta = deltap1\n",
    "    return x, res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now run this unit test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver terminated in 200 steps\n",
      "Terrific!  Your Nesterov routine works!\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUZdrA4d+bSSaTHtJISAiJ1BQSSgApBkGqNBsglpXVlXVt67rrZ1mxra6uva4uouCqi6CogICAKCpFegtNipSEhJAEQnp9vz9mUoAkBDKTmUye+7rmmjlnzpzzzAnkyduV1hohhBCiPi72DkAIIYRjk0QhhBCiQZIohBBCNEgShRBCiAZJohBCCNEgV3sHYAtBQUE6KirK3mEIIUSLsnnz5iytdfC5+50yUURFRbFp0yZ7hyGEEC2KUupIXful6kkIIUSDnCpRKKXGKaVm5Obm2jsUIYRwGk6VKLTWi7TW0/z8/OwdihBCOA2nbKMQQji/srIyUlNTKS4utncoLY7JZCIiIgI3N7dGHS+JQgjRIqWmpuLj40NUVBRKKXuH02JorcnOziY1NZXo6OhGfcapqp6EEK1HcXExgYGBkiQuklKKwMDAiyqJSaIQQrRYkiQuzcXeN0kUtcxYMJ1nPr7Z3mEIIYRDkURRy9YTq/ihdJu9wxBCtEJRUVFkZWUBMGDAgEs+z+zZszl+/Li1wgIkUZzFXRkplpKsEMJKysvLL+lza9euveRr2iJROFWvJ6XUOGBcp06dLunzRhcTRS4KrbXUfQohLugf//gHn376Ke3btycoKIjevXvzzTffMGDAANasWcP48ePp0qULzz77LKWlpQQGBvLpp5/Stm1bsrOzmTJlCidPnqRv377UXm3U29ub/Px8AF566SXmzZtHSUkJ1157LU8//TSHDx9m9OjRDBo0iLVr1xIeHs6CBQtYvHgxmzZt4uabb8bDw4N169bh4eHR5O/pVIlCa70IWJSUlHTnpXze5OJBhVIUluTjZfKxcnRCCFt5etEudh8/Y9Vzxrbz5clxcfW+v2nTJubPn8/WrVspLy+nV69e9O7dG4DTp0/z448/AnDq1Cl++eUXlFLMnDmTF198kVdeeYWnn36aQYMG8cQTT7B48WJmzJhx3jWWL1/O/v372bBhA1prxo8fz08//URkZCT79+9nzpw5vP/++0yaNIn58+dzyy238Pbbb/Pyyy+TlJRktXvhVImiqUwGE2jIzs2URCGEaNDq1auZMGFC9V/s48aNq35v8uTJ1a9TU1OZPHky6enplJaWVo9d+Omnn/jyyy8BGDNmDG3atDnvGsuXL2f58uX07NkTgPz8fPbv309kZCTR0dH06NEDgN69e3P48GGbfE+QRHEWk6sXlEFufha07WjvcIQQjdTQX/62Uruq6FxeXl7Vr++77z4efPBBxo8fz6pVq3jqqaeq37tQFbfWmkcffZQ//vGPZ+0/fPgw7u7u1dsGg4GioqKL/AaNJ43ZtXi4eQNwOj/bzpEIIRzdoEGDWLRoEcXFxeTn57N48eI6j8vNzSU8PByAjz76qHp/cnIyn376KQBLly7l1KlT53125MiRfPjhh9XtFWlpaWRmZjYYl4+PD3l5eZf0neojJYpaPIw+UAh5Ref/wIQQorY+ffowfvx4EhMT6dChA0lJSdQ1IelTTz3FxIkTCQ8P5/LLL+e3334D4Mknn2TKlCn06tWLwYMHExkZed5nR4wYwZ49e+jfvz9gbuT+5JNPMBgM9cY1depU7rrrLqs2ZquGik8tVVJSkr6UhYs+X/oqz2TO4vHIu5g85B4bRCaEsJY9e/YQExNj1xjy8/Px9vamsLCQ5ORkZsyYQa9evewaU2PVdf+UUpu11ue1gkuJohYvD/NfAwXF1u09IYRwTtOmTWP37t0UFxdz2223tZgkcbEkUdTi4xEAQGGpJAohxIX973//s3cIzUIas2vx8bQkirJ8O0cihBCOQxJFLX4+gQAUlxXYORIhhHAcDl/1pJTyAv4NlAKrtNaf2upafl6BKK0pLi+01SWEEKLFsUuJQin1oVIqUymVcs7+UUqpfUqpA0qpRyy7rwO+0FrfCYy3ZVye3n54aE1Rpe0GrgghREtjr6qn2cCo2juUUgbgHWA0EAtMUUrFAhHAMcthFbYMyt3dhGelpqRC1uAVQlzY4cOHiY+PP29/7SnDnYFdEoXW+icg55zdfYEDWutDWutS4DNgApCKOVlAA/EqpaYppTYppTadPHnykuJSLi6YtKZEl1zS54UQwhk5UmN2ODUlBzAniHDgS+B6pdS7wKL6Pqy1nqG1TtJaJwUHB19yEB6VimJddsmfF0K0LuXl5dx2220kJCRwww03UFhobuN86aWX6Nu3L3379uXAgQN2jrJpHKkxu67ZsbTWugD4faNO0MT1KAB8KhWnDFKiEKJFWfoIZOy07jlDu8PoFy542L59+/jggw8YOHAgt99+O//+978B8PX1ZcOGDfz3v//lgQce4JtvvrFufM3IkUoUqUD7WtsRwEUt06S1XqS1nlbXfCuN5Vvhwhl1aatSCSFan/bt2zNw4EAAbrnlFlavXg3AlClTqp/XrVtnt/iswZFKFBuBzkqpaCANuBG4qbmD8NVunHEpk1XuhGhJGvGXv62c+3uiarv2/pb+u8Re3WPnAOuArkqpVKXUHVrrcuBeYBmwB5intd51kecdp5SakZube8mx+WCiTEF+WT5smwOvxkGlTTtbCSFasKNHj1aXGObMmcOgQYMAmDt3bvVz1eyvLZW9ej1N0VqHaa3dtNYRWusPLPuXaK27aK07aq2fu4TzNrnqycfFvOBIdlE2rHgCzqRC6sZLPp8QwrnFxMTw0UcfkZCQQE5ODn/6058AKCkpoV+/frzxxhu89tprdo6yaRyp6qnJrNGY7WvwBTLIKsoiKvoKSJkPe7+ByMutFqcQwjlERUWxe/fu8/ZXLUv65JNPNnNEtuFIjdlNZo0ShZ/RH4C0M5ng4mbeeaRlN0QJIURTOFWisIYAD/MYjNTTaVBumcqjXLrLCiFaL0kU5wjwCsFFazJy06DMkigqSu0blBBC2JFTJQpr9Hry8gnEv7KS7IITkiiEEAInSxTWaKPw9A0gsKKCUyXZUG6ZHLBCpvQQQrReTpUorMHLL5DAigpyy05DWVWikBKFEKL1cqpEYY2qJw9vfwIrKsmryK9pzJZEIYSwkW3btrFkyRJ7h9Egp0oU1qh6Uh7+BFZUkK8L0WVS9SSEsK1LSRTl5c07H51TJQqrcPclqKKCclVBYdWSqFKiEELU4fDhw8TExHDnnXcSFxfHiBEjKCoq4uDBg4waNYrevXtzxRVXsHfvXgA+//xz4uPjSUxMJDk5mdLSUp544gnmzp1Ljx49mDt3LgUFBdx+++306dOHnj17smDBAgBmz57NxIkTGTduHCNGjEBrzUMPPUR8fDzdu3evnjJk8uTJZyWeqVOnMn/+/CZ9T6camW0VRi/8LVM7ZelSvAAqy0BraOETewnhrP614V/szdlr1XN2C+jGw30fvuBx+/fvZ86cObz//vtMmjSJ+fPnM2vWLN577z06d+7M+vXrufvuu/n+++955plnWLZsGeHh4Zw+fRqj0cgzzzzDpk2bePvttwF47LHHGDp0KB9++CGnT5+mb9++DBs2DIB169axY8cOAgICmD9/Ptu2bWP79u1kZWXRp08fkpOTufHGG5k7dy5XX301paWlrFy5knfffbdJ90ISxbmUwkP5ApBdWUYHFKDN1U+uRvvGJoRwONHR0fTo0QOA3r17c/jwYdauXcvEiROrjykpMQ/aHThwIFOnTmXSpElcd911dZ5v+fLlLFy4kJdffhmA4uJijh49CsDw4cMJCAgAYPXq1UyZMgWDwUDbtm0ZPHgwGzduZPTo0dx///2UlJTw7bffkpycjIeHR5O+o1MlCmvM9QRgcgsCTpFtcAGTLxTnmqufJFEI4ZAa85e/rbi7u1e/NhgMnDhxAn9/f7Zt23bese+99x7r169n8eLF9OjRo85jtNbMnz+frl27nrV//fr1eHl5nXVcXUwmE1deeSXLli1j7ty51etiNIVTtVFYozEbwNvTvER3tsEA7pZzSTuFEKIRfH19iY6O5vPPPwfMv9C3b98OwMGDB+nXrx/PPPMMQUFBHDt2DB8fH/Ly8qo/P3LkSN56663qRLB169Y6r5OcnMzcuXOpqKjg5MmT/PTTT/Tt2xeAG2+8kVmzZvHzzz8zcuTIJn8np0oU1uLnF43SmiyDwVyiAOn5JIRotE8//ZQPPviAxMRE4uLiqhukH3roIbp37058fDzJyckkJiYyZMgQdu/eXd2YPX36dMrKykhISCA+Pp7p06fXeY1rr72WhIQEEhMTGTp0KC+++CKhoaEAjBgxgp9++olhw4ZhNDa9JkTVV3xpyZKSkvSmTZsu+fPHvp/JLYdfZWhBIU96x8KRNfBACvi3v/CHhRDNYs+ePcTExNg7jBarrvunlNqstU4691gpUdQhsH1nAisqLFVPVSUKqXoSQrROkijq4Nm2S02ikKonIUQr51SJwhpTeADgE0pQRSXZBgNaShRCOCxnrDpvDhd735wqUVir1xNKEVhRQZbBhTPaZN4nJQohHIrJZCI7O1uSxUXSWpOdnY3JZGr0Z5xqHIU1+bv6UuIC+/Ir6QtSohDCwURERJCamsrJkyftHUqLYzKZiIiIaPTxkijqERzZH7LXsfVMuSQKIRyQm5sb0dHR9g6jVXCqqidrCu45FYBVeZY+yFL1JIRopSRR1CPQIxCAjJIC8w4pUQghWilJFPUI8ggCoMxgnsxLEoUQorWSRFEPf3d/XJQLymRJEFL1JIRopZwqUVhtHAVgcDEQYArAz6/SvENKFEKIVsqpEoXVxlFYBJoCMXqYq55y8vKtck4hhGhpnCpRWFuQRxDarQiAwydO2zkaIYSwD0kUDQj0CKSg4gwAx7IkUQghWidJFA0I9AjkZHE2lcDxrDMyVYAQolWSRNGAcK9wyivLOWkwUFRczLGcInuHJIQQzU4SRQPCfcIBSHVzw02Vs+5Qlp0jEkKI5ieJogER3uZJs1KNJvyMmjUHsu0ckRBCND9JFA1o590OheK4mxuRfq6sPZgl7RRCiFZHEkUDjAYjwZ7BpLq5EuHjSlZ+KXsz8qrfr6yUpCGEcH6SKC4gwjuCVFc3wj0rAPh+byYAezPOcNljS1h7QNothBDOTRLFBXTy78R+V4Wp5ASJEX4s330CgI/WHgFg85FT9gxPCCFszuEThVLqMqXUB0qpL+xx/ZjAGPKUJjXvOCPiQtl+7DSHTuaz+oB5VS1/Tzd7hCWEEM3GpolCKfWhUipTKZVyzv5RSql9SqkDSqlHGjqH1vqQ1voOW8bZkJjAGAD2lGYzqXcE7q4uTHhnTfWYioLSCnuFJoQQzcLWJYrZwKjaO5RSBuAdYDQQC0xRSsUqpborpb455xFi4/guqLN/Z1xxYbdBE+xayB+uiCavuJzoIC8ACiVRCCGcnE3XzNZa/6SUijpnd1/ggNb6EIBS6jNggtb6eWDspV5LKTUNmAYQGRl5qac5j9FgJMYrjE3FRXDmOA+NjGd8YjhtPN0Y/NIqikrLrXYtIYRwRPZoowgHjtXaTrXsq5NSKlAp9R7QUyn1aH3Haa1naK2TtNZJwcHB1osW6B/ckxR3I2dyDgLQNdSHEF8TXu4GqXoSQjg9eyQKVce+egckaK2ztdZ3aa07Wkod9Z/YigsX1TYgYjAVSrEhY/1Z+z2MBookUQghnJw9EkUq0L7WdgRw3BontvbCRVUSIpPxq6hgefYOyMuAmcPh+2fxMroSlrsVvnvaqtcTQghHYo9EsRHorJSKVkoZgRuBhdY4sa1KFG5unows0awqPEbhqucgdQPsXoCH0cD/HX8AVr9q1esJIYQjsXX32DnAOqCrUipVKXWH1rocuBdYBuwB5mmtd1njerYqUQCMNbShiEqW5e437/Bog5exVl+ACmnUFkI4J1v3eppSz/4lwBJbXtvaenhF0Kl0H/P0Sa4FKC3Ew8tQc0BFCRhsejuFEMIuHH5k9sWwVdUTgPINY2J+MSmqlF1GNygrwMtYK1GUl1j9mkII4QicKlHYsuoJnzDG5ZzAQ2vm+fqYSxS1q54kUQghnJRTJQqb8gnDR2uuzitgiZcnZ8qLzi5RVEiiEEI4J6dKFLasesInDIDJeXkUu7iwyKjxdKt1+6REIYRwUk6VKGxa9eTbDoCY0jISSsuZ5+OFl1tlzfvlxda/phBCOACnShQ25Vczy8hE7cUhoxunS7fXvF9eaoeghBDC9iRRNJbJv/rlKFM4vhUV7MhfUfO+lCiEEE7KqRKFTdsoVM0UVSbfcCbkF7CrcAtZBsstlMZsIYSTcqpEYdM2itp8wph0Jp8KXcGX3t4ApGfbIDkJIYQDcKpEYXOuJvOzTyhR5eX08+/GF77eVADvrtzNyTwpVQghnI8kiothspRUvM0L700M6kW6qysbTO6UFBdyy8z1kiyEEE6nwUShlOrV0KO5gnQYt8yH3lOhTTQAg91D8KisZLmXJ3f2D+doTiGTZ6zjWE6hfeMUQggrutAsdq808J4GhloxliZTSo0DxnXq1Mk2FwjtDuPegJxDAJiKcrmysIjvvDz5exsDH93elz98tJHxb6/m3zf3pn/HQNvEIYQQzajBEoXWekgDD4dKEtCMjdluXubnwmxGFhRy2mBgQ94h+rYpYJPfw8R4nOKWD9Yz8+dDVFbWu3ifEEK0CI1uo1BKxSulJimlflf1sGVgDs3oaX4uzGFQURGelZUsP3MAds7DmPsbs+N3MLRbCM8u3sOtH67n+Oki+8YrhBBN0KhEoZR6EnjL8hgCvAiMt2Fcjs3NkigKMnHXMKSwiO8Kj1Lm0QYAY3EOM27tzfPXdWfr0dOMfP0nPv7lCBVSuhBCtECNLVHcAFwFZGitfw8kAu42i8rRuRgguBsc+A6AkQWF5FaWsqEw1fx+YRZKKab0jWTJ/VcQ386P6V+ncN27a9m/c71MICiEaFEamyiKtNaVQLlSyhfIBC6zXVgtQOyE6pcDiorwUgaW51hWdC3Iqn4vKsiL/93Zj9cn9+B0zkmivhjNl7NfJStfkoUQomVobKLYpJTyB94HNgNbgA02i+oS2XQKj3PFXVf90l3DINdAfsz/jUqAwqyzDlVKcU3PcBZN64mbqmDf4WNc+dIq3vnhAMVlFbaPVQghmqBRiUJrfbfW+rTW+j1gOHCbpQrKoTRbryeAkG7QfZL5tcHIYIMf2RVF7DIaIT8T9PntEb6u5QD8cWAE/TsG8tKyfQx9eRVfbU2V3lFCCIfV2Mbs5KoHEAn4W163btfNgAf3QJsokpUnBhQ/eHmYZ5ItOnX+8Za2iQATvP+7JD6bdjkB3kb+Mnc7V7/5M8t3ZaDrSDBCCGFPFxpwV+WhWq9NQF/MVVAON5aiWSllXtDI1R2/igp6uPryo0cJ95/KhZI88Aw4+/iqqcgtCePyywJZeM8gFu04zuvf7Wfax5tJbO/P30Z0YVCnIFStGWuFEMJeGlv1NK7WYzgQD5ywbWgtiMEdyou50sWXX92NHHc1QGm++b3KCsjab35d1dupoqz6oy4uigk9wlnxl2RevD6BrLwSbv1gA5Nn/MLGwznN/EWEEOJ8lzopYCrmZCHAPKtseSnJ2txjeI2HCUryze0UM4fB20lw+ljNmhV1rF3hanBhUp/2fP+3wTwzIY7fsgqY+N46bv1gvSQMIYRdNarqSSn1Fua5ncCcXHoA2+v/RCtj9ISMFKJNPoSYyvnFZGJiaR5k7oHjW8zHFGbVKlHUv2yqu6uB3/WPYmLv9nz8y2H+8+MhJr63jssvC+D+oZ3p3zFQqqSEEM2q0d1jMbdJbAbWAQ9rrW+xWVQtzYD7oSgHdXIfl5cr1nuYqCw+AwWZNceUFtRqo7jw+toeRgPTkjvy88NDeHxMDIdOFnDTzPXc8N46Vu3LlEZvIUSzaVSJQmv9ka0DsQabzx5bn+grzDPLpm7kchcfFhoK2ZN7iDhq/TIvLWxUieJcnkZX/nDFZdxyeQc+33SMd1cdZOqsjSRE+HHf0M4MiwmREoYQwqYutB7FTqXUjvoezRVkYzXrOIpz+bUHoL9HOwB+Of0rFNZqWyjNrylRXESiqGJyM3Br/yhWPTSEF67rzunCMu787yZGv/EzC7alUV5R2eSvIIQQdblQ1dNYYBzwreVxs+WxBPjCtqG1MP7mRBHkHUan0lLW5f8Ghdk175cWXFKJ4lxGVxdu7BvJ938dzCsTEymrqOTPn21j8EurmLXmNwpLy5vyLYQQ4jwXWo/iiNb6CDBQa/1/WuudlscjwMjmCbGFsJQoqCijf3EpW4szKa7dRlFWq+rJCpMCuhpcuL53BCv+MpiZv0sizM/E04t2M+CF73l1+T6ZS0oIYTWNbcz2UkoNqtpQSg0AvGwTUgtVlSgKs7m83IVSKtmSfxR8w837z6p6Kqv7HJfAxUUxLLYtX/xpAF/c1Z8+UQG8+f0BBr7wPY9/vZMj2QVWu5YQonVq7MjsO4APlVJVlf+ngdttE1IL5RVkfi7MJglPXFH8UnyCAT5hkJdhrnpSBvMxdYyjsIakqACSogI4kJnP+z8dYt7GVP63/iij48OYlnwZie39bXJdIYRza+zI7M1a60QgAUjUWvfQWm+xbWgtTEgseIfC0Ol4unvTQ3mytvKMOYEYvc7uHtuENorG6BTizb9uSODnh4cwLbkjP/16kgnvrOGGd9eyZGe6NHwLIS5KgyUKpdQtWutPlFIPnrMfAK31qzaMrWUxesLf9plf//gvBuoS3nAp4KTJh+CqROFm6S7biHEU1tDW18Qjo7txz5COzNuUyuy1v3H3p1sI9/fgtgEdmNwnEj8Pt2aJRQjRcl2oRFHVDuFTz0PUxd2HQSXmdSbWqJI6ShTN29DsY3LjjkHRrPrbEP5za28i2njwzyV76f/8Sp5YkMKhk/nNGo8QomVpsEShtf6P5fnp5gnHSbh70/XgSoLah7PGpYxr3DzNicJgNL9vxcbsi2FwUYyMC2VkXCgpabnMWnOYzzYc47/rjnBVtxBuHxTNAJkiRAhxjsauR/GiUspXKeWmlFqplMpSSskUHvXxa4/CvETq2sKjVBi9zN1jK6zXPbap4sP9eGVSIqsfGcKfr+rMtmOnuXnmeka9/jP/W39UxmMIIao1tnvsCK31GcwD8FKBLpy9RoVNKaWuUUq9r5RaoJQa0VzXvWTJD4FHAFcExHGmNI+dbpZpx+uYZtzeQnxM/GV4F9Y8MpSXbkjAxUXx2Fc76ffPlTy9aBcHpVpKiFavsYmiqsXzamCO1rrR814rpT5USmUqpVLO2T9KKbVPKXVAKfVIQ+fQWn+ttb4TmApMbuy17cbDH/6yi/4TZuHm4sYyVWzXNorGMLkZmJjUniX3D+KLu/oztFsIn/xyhKte+ZFbZq5n2a4M6S0lRCvV2HEUi5RSe4Ei4G6lVDBQ3MjPzgbeBv5btUMpZQDewbz+diqwUSm1EDAAz5/z+du11lVDnB+3fM7xGT3xM3oyOGIwS479yF9LS3C1whQetqaUqh6P8fiYWOZuPMqn64/yx483087PxE39IpncJ5JgH3d7hyqEaCaqsdNVK6XaAGe01hVKKS/AR2ud0cjPRgHfaK3jLdv9gae01iMt248CaK3PTRJVn1fAC8AKrfV3F7peUlKS3rRpU2NCs7lVx1Zx3/f38a9ThVzt1QHSNpvfmJ4Nhsbmafsqr6hk5d5MPl53hNUHsnAzKEbHh/G7/h3o3aGNNH4L4SSUUpu11knn7m/swkWewD1AJDANaAd0Bb65xHjCgWO1tlOBfg0cfx8wDPBTSnXSWr9XR4zTLLERGRl5iWFZX3JEMp3c/JjhVcbwvIzqOjwqSltMonA1uFT3ljp4Mp9PfjnCF5tTWbj9ODFhvtxyeSQTeoTj7d4yvo8Q4uI0to1iFlAKDLBspwLPNuG6df0JWm/RRmv9pta6t9b6rrqShOWYGVrrJK11UnBwcBNCsy4X5cK9sbdx0OjGhy55NW84YDtFY3QM9ubJcXGsf+wqnr+uOwB//yqFvs99x6Nf7mBH6mlZVEkIJ9PYRNFRa/0iUAagtS6i7l/2jZUKtK+1HQEcb8L5APPCRUqpGbm5uU09lVVdlfgHRpfAO/5+rPT0MO90oJ5Pl8LT6MqUvpEsuX8QX909gLEJYXy99Tjj317DmDdX8/EvRzhT3LK/oxDCrLGJolQp5YHlr36lVEegKX8SbwQ6K6WilVJG4EZgYRPOB9h54aKGKMXTA/9B9woX/i8kmM3u7g4xlsIalFL0jGzDizcksv7vV/GPa+IBmP51Cv2eW8n/fbGdrUdPSSlDiBbsgo3ZlobkWzHPIBsLLAcGAlO11qsueAGl5gBXAkHACeBJrfUHSqmrgdcx93T6UGv93KV/jeprVS2Feuf+/fubejqrO1V8it99fS0ni07yQfKrxHV0ziU9tNbsTMtlzoajLNh2nMLSCrqF+nBTP3NbhswvJYRjqq8xu1G9npRSm4ERwOWYq5x+0VpnWT1KK3GkXk/nytgyi6lbXiTf058PR39ElzZd7B2STeWXlLNw23HmbDjKzrRcTG4ujE1ox5S+kfSK9JceU0I4kKYmineA2VrrjbYIztocOVGwdzHHvriVqZ3iKUXzYvKL9G/X395RNYudqbnM2XiUBVvTKCitoHOIN5OS2nNNz3AZlyGEA2hqotiNedqOI0AB5lKF1lonWDvQpnD0qicA9q+AT2/g6LXvcs+h/3H4zBFiA2Np79OeNu5tCPAIIKltEkltk5z2r+2CknIWbT/OvE3H2HL0NK4uiiHdQpiU1J4ruwbjZmhs05kQwpqamig61LXfsp62w3HoEkX6dvhPMri4UaTLmRsSwU+X9SWzMJOc4hzySs1daIe0H8Jzg57Dx+jcs7kfyMzj802pzN+SRlZ+CUHe7lzXK5yJvSPo3Na5v7sQjqZJiaKlcehEAbD0YVhfazjIUzXdeQvLCpm7b8YI09MAACAASURBVC5vbnmTrgFd+XDkh3i6edohyOZVVlHJj/tOMm/TMb7fm0l5paZnpD+TktozNiEMH5M0gAtha60iUbSIqieAzbNh0Z9rtp86f9zHD0d/4IFVDzC8w3BeSn7Jaauh6pKVX8LXW9OYt+kYv57Ix+TmwtXxYUxMak+/6ABcXFrPvRCiObWKRFHF4UsU2+bA13fVbNeRKAA+2PkBr295nb/2/itT46c2T2wORGvNjtRc5m06xsLtx8krLqd9gAcTe7fn+t4RhPt72DtEIZxKk+Z6ElZmaFw1yu3xt7Mrexevb3mdHiE96BHSw8aBORalFInt/Uls78/0sbEs25XBvE3HeHXFr7z23a8M6BjIdT0jGBUfipfMMyWEzUiJwh72LIK5tRYIrKdEAXCm9AyTFk2iUlfy+bjP8XN3sFHndnAsp5D5W1L5amsaR7IL8XAzMDo+lOt6RdC/YyAGqZoS4pK0iqqnFtNG8esy+N+kmu0GEgVASlYKty69leTwZF4f8nqraq9oiNaaLUdP8cXmNL7ZYa6aCvMzcU3PcK7vFU6nEOk1JcTFaBWJoorDlygO/gAfX1OzfYFEAfDx7o95ceOL3BxzMw/3eViSxTmKyypYuSeTL7eksurXk1RUahIi/LiuZzjjEtsR6C0D+oS4EGmjcCSu5/zSqii/4NoUt8TcQnpBOh/v/hiF4m9Jf8PgYrBhkC2Lyc3AmIQwxiSEkZVfwoJtx/lySypPLdrNs4v3MKRbCNf3CmdItxDcXeW+CXExJFHYg8F49nZ5ERgariZRSvFQ0kNorflkzyek5qXyQvILeLl52TDQlinI2507BkVzx6Bo9mac4astaXy1NY0Vu0/g5+HGuMQwrusVQc/2MteUEI0hVU/2kL4D/nNFzfZffwWfto3++Gd7P+OFDS8Q7RfNW0PfIsInwgZBOpeKSs3qA1l8uSWVZbsyKC6rJDrIi+t6hnNNz3DaBzj/oEYhLqRVtFG0mMbszL3w71orv96/DQKiL+oU646v468//hWDMvDala+RFHrez1bUI6+4jKUpGXy5JZVfDuUA0LtDG67p0Y4xCe0I8DJe4AxCOKdWkSiqOHyJIucQvNmzZvtP66Bt7EWf5siZI9y78l5S81J5ZuAzjOs4zopBtg6ppwpZuP04C7YeZ9+JPFxdFMldgpnQox0jYkPxMEp7hmg9pDHbkRjOacwuK7qk03Tw7cCnYz7lwVUP8viax/E1+jK4/WArBNh6RLTx5O4rO3H3lZ3Yk36Gr7elsXDbcb7fm4mX0cDIuFAm9AxnYMdAXGVWW9FKSYnCHvJPwsudarZvWwTRyZd8usKyQm5fdjsHTx9k1qhZxAfFWyHI1quyUrP+txwWbEtjyc50zhSXE+RtZGxCO67pGU5ihJ80ggunJFVPjqQ4F16IrNm+aR50adqyqNlF2dy85GZKK0qZM2YObb0a3zgu6ldSXsEPe0+yYFsaK/dmUlpubgQfn2hOGtFB0utMOA9JFI6krAieC63Znjgb4q5t8ml/PfUrty65lWi/aGaPmo3J1dTkc4oauUVlLEvJ4Ottaaw7lI3WkBjhxzU9wxmb0E5W6RMtXqtIFC2m11NlJTzTpmb7mnehx01WOfUPR3/gzz/8mVFRo/hX8r+kisRGMnKLWbg9ja+3Hmd3+hkMLoqBnYK4pkc7RsSF4i2TFIoWqFUkiioOX6IAeCYQKsvNr69+GfreabVTz9w5kze2vMF9Pe9jWsI0q51X1G3/iTy+3pbGgm3HST1VhMnNhau6tWVcYjuu7BqMyU16TomWQXo9ORqDsSZRXGKvp/rcEX8HB04f4K2tb9HRryNXdbjKqucXZ+vc1oeHRnbjbyO6svnIKRZuP86Sneks3pmOj7srI+JCGZcYxsBOQbIeuGiRJFHYi8EIZYXm11ZOFEopnh7wNEfPHOXR1Y/ysc/HdA3oatVriPMppUiKCiApKoAnxsay7lA2C7cd59tdGczfkkqAl5HR8aGMT2xHnyhZqU+0HFL1ZC8vdYaCTPPrgX+G4c9Y/RInC09y4+IbMSgDc8bMIdAj0OrXEBdWUl7Bj/tOsnD7cb7bc4LiskpCfU2MTQhjfI92dA+X7rbCMUgbhaN5LR5yj5lf97kTxrxsk8vsyt7F1KVTiQmMYeaImRjPnZBQNKuCknK+23OCRdvT+fHXTMoqNFGBnoxLbMf4xHZ0bitraAj7kUThaN7sBTkHQblAwo1w7bs2u9S3h7/loR8fYmTUSJ4f9DxujVyKVdhWbmEZy3ZlsHD7cdYezKJSQ7dQH8YltmNcQjsiA2WiQtG8pDHb0VT9Ze/RBsoKbHqpUVGjyMjP4JXNr3C6+DQvD34Zf5P/WcdorTlVcopTxafwMfoQ7BEs1SE25ufpxqQ+7ZnUpz2ZecUs3WlOGi8t28dLy/bRo70/4xPbMTYhjBBfGRMj7MepShQtZhwFwH+SIX07BHaCNtFwyxc2v+Sig4t4Yu0TeLt5c23na+no15GMggz25uxlS+YWcopzqo+N8I7gxm43clPMTbi5SAmkOaWeKuSbHeks3GYeo6EUXB4dyPge7RgVF0obmd1W2IhUPTmamcMgdSOEJ5lXvPv9kma57L6cfby59U3WpK2hQlcA5qTQM6QnsYGxBJgCOFVyipVHV7IxYyOxgbE8f8XzXOZ3WbPEJ852IDOfRduPs2j7cQ5lFVTPbjs2IYxhsW3xNUkSF9YjicLRzBoDR1ZDp2FQkAV//LFZL19YVkhWURaBHoH1rpK34sgKnv3lWSp1Jf8Z/h9iAy9+KnRhHVprdh0/w6Idx/lmezppp4swGlwY3NWcNK6KaSujwUWTSaJwNB9fC7/9DN3GQOZuuHej9c69bymUFkD3G5p8qmNnjnHH8jvIL81n1qhZMh7DAWit2XrsNN9sT2fJznQyzhTj7urCkK4hjEkI46qYEDyNkjTExZNE4Wj+N9mcKOKugUOr4MHd1jv3U36W51yrnC49P51bltyCq4src8bOIcAUYJXziqarrNRsPnqKxTvMI8FP5pVUTyEyNiGMK7uGyOJLotGk15OjMRjB1Qhunua//h1YmHcYbwx9g6nfTuWvq/7KzBEzMbjILx9H4OKi6BMVQJ+oAKaPjWXj4RwW70hnaYo5cXgaDQyLacuYhDAGd5F5p8SlkURhLwajeaU7o2fNVB4OLD4onumXT+fxNY/zYcqH3JlgvUkMhXUYXBSXXxbI5ZcF8uS4WDb8lsOiHel8m5LOwu3H8XZ3ZXhsW8Z0D+OKLkG4u0rSEI0jicJeAjvCmTRw84KKUqgoh8xdUF4C7fvaO7o6je84ntVpq/n3tn8zoN0A4oLi7B2SqIerwYUBnYIY0CmIZybEse5gNot3pPPtrgy+2pqGj8mVEbGhjE0MY2DHIIyuMlmhqJ+0Udjb2rdh+d/hkaM1q941tW3Bym0UteWW5HLDohtwN7gzb+w8PN1k9HBLUlpeyZqDWSzekc6yXRnkFZfj5+HGqLhQxiSE0b9joMxw24rV10Yh/yLszWj5RVu7naL4jHXOXVFmnfPU4ufuxz8H/ZOjZ47y4sYXrX5+YVtGS++olycmsunxYXxwWxJDu4WweGc6v/twA32f+45Hv9zJmgNZlFdU2jtc4SCk6sneqsYwnDlesy9jB0QNavq5SwvAw//Cx12kPqF9uD3+dj5I+YArwq+Q9S5aKHdXA1fFtOWqmLYUl1Xw468nWbwjnQXb0piz4ShB3kZGxYcyNsE8LbpBpkVvtRw+USilYoA/A0HASq217WbPs4eqEkX6tpp96dutkyjKCm2SKADu6XEP69LXMX3tdKL8oujo39Em1xHNw+RmYGRcKCPjQikqrWDVvky+2ZHOF5tT+eSXo4T4uDM6PpSru4eRJEmj1bFp1ZNS6kOlVKZSKuWc/aOUUvuUUgeUUo80dA6t9R6t9V3AJOC8urMWr6qOP317zb5MK42pKLVdbyo3gxuvDH4Fd4M7d313FxkFGTa7lmheHkYDo7uH8c7NvdgyfThv39STXpFt+GzjMSbP+IXLn1/J9K9TWHcwm4pK52vjFOezdYliNvA28N+qHUopA/AOMBxIBTYqpRYCBuD5cz5/u9Y6Uyk1HnjEci7nYrRUPR1dDy5uENTFPKWHNdh4VtoInwj+fdW/+f2y33PT4pu4r+d9dA/qTmZhJjuydrD95HYOnD5ASXkJPkYf+ob15daYW7nMX+aNaik8ja6MTWjH2IR2FJSU8/3eTJbsTOfzzcf4+JcjBHm7Myq+LVd3D6NvVACu0hDulGze60kpFQV8o7WOt2z3B57SWo+0bD8KoLU+N0nUda7FWusx9bw3DZgGEBkZ2fvIkSNWid/mMlLgvYHm111GmRugi0/Dnd9f+jmrej3dvgwiL296jBewL2cff1/9d/ad2le9T6Ho6N+RmIAYPN08OVl4knXp6yitKGVKtyk80PsB3A3uNo9N2EZhaTk/7D3Jkp3pfL83k6KyCgK9jIyMD2VM9zD6RUvSaIkcaWR2OHCs1nYq0K++g5VSVwLXAe5AvVOsaq1nADPA3D3WGoE2i+Bu4BkIhdnQfSLsXwHZlinSsw9C/gnoMABOHYEf/gnjXgc3j/rPVzvxN9OI764BXfl83OekZKWQmp9KG1Mb4gLj8DGevVpbTnEO72x9h0/2fMKOkzt4bchrhHiGNEuMwro8ja6MSQhjTEIYhaXlrNpnThpfb03jf+uPEuBlZGScuaTR/7JASRotnD0SRV2tYPX+YtdarwJW2SoYuzO4wh9WwsaZ0G0sHN9aU/X0wz/h8M/wt19h1fOw4zPoOBQSJ9d/vtpdYptxxLdSiu7B3eke3L3eYwJMAUzvP53+7frz2OrHuGXJLcwYPoMov6hmi1NYn6fRlau7h3F19zCKSiv48ddMFu/MYMG248zZcIw2nm6MiA3l6oQwBsg4jRbJHokiFWhfazsCOF7PsRel1sJF1jhd8wmIhpHPmV97BZt/wZcWwKnD5hJFWTG4+5rfz0tv+FwVJTWvbdiY3RTDOgyjnXc7/vTdn7jt29t4d9i7MoW5k/AwGhgVH8ao+LDqLrdLdqbzzY7jzN10DH9PN0bEmksaA2REeIthj5/SRqCzUipaKWUEbgQWWuPEWutFWutpfn5+1jidfXgFm58LTkJuqvn1mTSomoQv+0DDnz+rROG4kw3GBsby0aiPcDe4c/uy29mYYcVp1oVDqOpy+8aNPdk8fTgzbu3NkK4hLN2ZwdRZG+nz3Hf87fPt/LA3k9JyGdznyGzdPXYOsA7oqpRKVUrdobUuB+4FlgF7gHla611Wut44pdSM3FzrT13RbKoSRW4a5Fu6nJ5Jq6mOqqvrbGUlvN0Xdn5hnjeqioOWKKpE+UXx39H/JcQzhLtW3MWKIyvqPba8spyCsgLKKq0/2lzYnsnNwIi4UF6b3INN04cx83dJXNUthGUpGfx+9kZ6P7uCB+dtY+WeE5SUV9g7XHEOm1Y9aa2n1LN/CQ00TDfheouARUlJSS13alOvQPNz7QF4uWnmxm6AzL3mBmtVq6kna5/58c2DcNfPNftbwKy0oV6hfDTqI+5ZeQ8PrnqQKyOuZETUCHyMPqTlp7E3Zy87Tu7gUO4hAFxdXElqm8Tv437PgPABdo5eXAp3VwPDYtsyLLYtJeUVrDmQxeIdGazYncGXW9LwcXdlmKV66orOQTI1ugNw+JHZrU5ViSJtc82+M6lQaClRlBVASR64+9Qki6O/mJ+Du5xd9eTg61xUaWNqw+xRs5m9azaf7P6EVamrqt/zd/cnITiB4R2G4+3mTXZxNssOL+OP3/2RSV0m8Wi/R3F1kX/GLZW7q4Gh3doytFtbSsu7s+ZgFkt2pLN89wm+2pqGt7srw2JCuLp7GMmynobdONX/sBbbmF2bTzswesOeRZYdytxWUZANygV0JZz6Df6TDGNegT5/qEkUJv+zq55aQImiitFgZFrCNO6Iv4NDuYcoLi8mzDuMQFMgSp3dUe6+nvfx9ta3mbVrFnmlebyQ/AIuShpFW7qqCQuHdA3hufJK1h7MYslOc9L4ettxvIzmuamu7h7GlV0laTQnp0oUTlH1ZHA1j5vYv9w8vsI3HHJ+M1c9hcTCiRRI22I+dsn/mRNFxk7zdvHpFtHrqSEGFwOd23Ru8BijwciDSQ/i6+7LG1ve4DL/y7gr8a5milA0B6OrC1d2DeHKriE8V1HJuoPZLNlpnhp94fbjeBoNDO0Wwuj4MIZ0C5Y1wm1M7q4jCks0J4rLhoC7N2z9FCrLICTGnCiqEoO2NPrlnzA/F+acXfVUXtS8cTezO+Lv4MDpA7y3/T0GthvY4BgO0XK5GVxI7hJMcpdg/nFNPOsP5bB4ZzrLd2XwzY50TG4uDO4SzOj4MIbGhOBrcrN3yE7HqRKFU1Q9AfS4GX79FoY8BodXw+bZ5v0hMebnE7XmWDx1pKahu+iUeYW8Ki2wRHExlFL8vd/f2Zi+kec3PM+nV396XjWVcC5uBhcGdQ5iUOcgnr0mng2/5fBtSjpLUzJYtusERsv7o+JDGRHbFn9Po71DdgpOlSicouoJzAPw7lptfl27naHTMPj+2ZoSBcDBlYAGjzbmqqfaiaIFtVFcKh+jD/f2vJcn1j7BssPLGBU96qI+/3Pqz8z7dR7HzhyjnXc7ru9yPUPbD5WE0wIYXBT9OwbSv2MgT46LY+uxUyzdmcHSlAy+35vJY5b3R8eHMSKuLUHeMrfYpZIWQEcXHANGH0icYq6S8gw8OwEcXmN+DupqbuiuKl2Y/KDMuaueqozvOJ7ObTrz+pbXKa3dmN+A8spynl73NHevvJs92XuI8oviUO4hHvjhAR5f8zjlleU2jlpYk4uLoneHAB4fG8vqh4ew8N6B/OGKyziWU8hjX+2k73PfMfk/6/ho7WEycovtHW6L41QlCqdkcIWH9oOrybztFWIete3fwVx6OGJJFMFd4dgvNYP0TP6tJlEYXAz8pddfuHvl3Sw6uIjru1zf4PFaa5795Vnm75/PHfF3cE+Pe3AzuFFeWc6MHTN4d/u7uLm48dSAp5rnCwirUkqREOFPQoQ/D4/qyt6MPJamZLB0ZzpPLtzFkwt30SvSn6u7hzEqPpSINrLu+4U4VaJwmjaKc9WeLdY3DDJ3gXeIOXkctgywC+5qfs6zNGyb/FrMOAprGBQ+iNjAWD5I+YAJnSY0OLZi/v75zN8/nzu738n9ve6v3u/q4srdPe6mrLKMmTtn0i+sH6OjRzdH+MJGlFLEhPkSE+bLg8O7cCAzr7p66tnFe3h28R66h/sxunsoo+PDiA7ysnfIDsmpqp6cYq6nC+lxs/k5c695kaMqQZZEUdUDysO/VbRRVFFKMa37NI7lHWP54eX1HpdVlMUrm16hX1g/7u15b53H3NPjHhKCE3hhwwvklebZKmRhB51CfLjvqs4s+fMV/PjQlTwyuhsuLooXv93HkJdXMer1n3jju/38eiIPW6/V05I4VaJoFWKvgS6jzbPNtutZsz/Qsmrc6aPmZ1PrShQAQyKH0NGvI+/vfJ9KXfckc29seYPiimIe7/d4vYP0XF1ceazfY+QU5zArZZYtQxZ21CHQi7sGd2TBPQNZ88hQpo+Nxcfkyusrf2XEaz8x7NUfeXnZPlLSclt90pBE0dK4uMBNn0Hv28wN3L4R0CYa/CLBxRUy95iP82g9bRRVXJQLd3Q3j6348diP572/8+ROvj7wNbfG3nrBNTDiAuMY0WEEn+39jPzSfBtFLBxFuL8HdwyK5vO7BrD+0av4x4Q42vqa+PeqA4x9azWDX1rF80v2sPXoqVaZNGy+FGpzqtVGcef+/fvtHU7zqCg3D6xz94G3kmpWxxtwH6x9C6ZnmxvEW4nyynLGfjWWAFPAWeMqKnUlNy++mROFJ1h07SK83C5cF70rexc3fnMj/9fn/7g19lZbhy4cUE5BKSt2Z7BkZwZrD2ZRVqFp52diZLy5TaN3hzYYXJynK3V9S6E6VYmiVbRRnMvgak4SAIG1GvFNlnvg5KOzz+Xq4srt8bezM2sn6zPWV+9fcGABKdkp/KX3XxqVJMBcquge1J0v93/ZKv+KFBDgZWRyn0g+ur0vm/4+nFcmJhLbzo9P1x9l0n/WcfnzK3n8652sOZBFeYXzrqnhVImi1QvsWPPa5G9+dvLR2XWZ0GkCwR7BzNwxE4Dsomxe3/I6PYJ7MPaysRd1rms7X8uB0wfYlW2VJVNEC+bn6cb1vSOYeVsSW6YP580pPekT1Yb5m9O4eeZ6+jz3HQ9/sYNV+5xvIabWUyfRGgTVmkzP6G1+bmUN2gDuBndui7uNlze9zCe7P+H7Y9+TX5rPE/2fuOgR1yM6jOCfv/yT5YeXEx8Ub6OIRUvj7e7K+MR2jE9sV71O+NKUDBbvTGfupmP4mFwZHtOW0U6ypoYkCmcSOwFO7jOPr6gae9HKGrSrTO46mZ/TfuZfG/+Fq4srzwx45oKz0tbFz92Pfu36seLICv7S+y8ytYc4T+11wkvKK1i9P4ulKRms2H2CL7em4WU0MKRbSPX06C1xptuWF7Gon0cbGPW8+fWvlrEErTRRmFxNvDvsXbZnbifEM4RI38hLPtewyGE8ve5pDuUeoqN/xwt/QLRa7q7mNTOuimlLmWV69KUpGefNdDsqPpSrYtq2mJlunSpROO3I7EtRXaJofVVPVdxc3EgKPa8Dx0Ub2G4gAGvS1kiiEI1We3r02jPdfrvLPNOtm0ExsFMQo+JCGR7blkAHnrTQqbrHVklKStKbNm2ydxj2lbYZ3h8KN82DLiPtHU2LN/7r8bTzasd7w9+zdyiihaus1GxLPc23KRksTUnnWE4RLgr6RgcwOj6MkXGhhPqZ7BJbfd1jnapEIWpxs0x01opLFNZ0edjlfH3ga8oqy3BzaRnVBcIxubgoekW2oVdkGx4d3Y3d6Wf4NiWDb1Myqict7Bnpz+j4UEbFhREZaP9JCyVROKtW3phtbb3a9mLO3jnsy9knvZ+E1SiliGvnR1w7P/46oisHMvNZtstc0vjnkr38c8leYsN8GRUfyuj4UDq39bFLnJIonFXVoLJVL0C3MTUD8MQl6RXSC4AtJ7ZIohA20ynEm04hnbhnSCeO5RRakkYGr674lVdX/ErHYC9L0ggjrp1vs/XCkzYKZ6U1rJhunsbj5vnQeZi9I2rxRs0fRWxgLK9e+aq9QxGtzIkzxSy3JI31v+VQUamJaOPBqLhQRsWH0iuyDS5WmEpE2ihaG6Ug6Q5zoijItHc0TqF7UHd2nNxh7zBEK9TW18St/aO4tX8UOQWlfLf7BEtT0vlo3WFmrv6NYB93Rsa1ZXR8GP2iA3A1WHfSDUkUzsw7xPxctUaFaJK4wDi+Pfwtp4pP0cbUxt7hiFYqwMvIpD7tmdSnPWeKy/hhbybfpmQwf3Man/xylFlT+zCkW4hVrymJwpkZvcxTeeRLicIaYgNjAdidvZuB4QPtHI0Q4GtyY0KPcCb0CK+eSmRAp0CrX8epJgVUSo1TSs3Izc21dyiOwztEShRWEhMYAyATBAqHVDWViLur9eeVcqpE0SqnGb8Q77ZSorASH6MP9/S4h54hPS98sBBORKqenJ13iHl9bWEVdyXeZe8QhGh2TlWiEHXwkqonIUTTSKJwdt5tofg0lJfYOxIhRAslicLZhSVC/PVQXmzvSIQQLZS0UTi7LiPMDyGEuERSohBCCNEgSRRCCCEaJIlCCCFEgyRRCCGEaFCLSBRKKS+l1Gal1Fh7xyKEEK2NTROFUupDpVSmUirlnP2jlFL7lFIHlFKPNOJUDwPzbBOlEEKIhti6e+xs4G3gv1U7lFIG4B1gOJAKbFRKLQQMwPPnfP52IAHYDdhntXEhhGjlbJootNY/KaWiztndFzigtT4EoJT6DJigtX4eOK9qSSk1BPACYoEipdQSrXVlHcdNA6YBREZGWvNrCCFEq2aPAXfhwLFa26lAv/oO1lr/HUApNRXIqitJWI6bAcywHHtSKXXkEuMLArIu8bPNRWK0DonROiTGpnOU+DrUtdMeiaKuhV0vuHC31np2Yy+gtQ6+mIBqU0ptqmvNWEciMVqHxGgdEmPTOXp89uj1lAq0r7UdARy3QxxCCCEawR6JYiPQWSkVrZQyAjcCC+0QhxBCiEawdffYOcA6oKtSKlUpdYfWuhy4F1gG7AHmaa0daW3JGfYOoBEkRuuQGK1DYmw6h45PaX3B5gEhhBCtWIsYmS2EEMJ+JFEIIYRokCSKWi5hapFmoZQ6rJTaqZTappTaZNkXoJRaoZTab3lu08wxnTc9S30xKbM3Lfd1h1Kql53ie0oplWa5j9uUUlfXeu9RS3z7lFIjbR2f5ZrtlVI/KKX2KKV2KaX+bNnvSPexvhgd5l4qpUxKqQ1Kqe2WGJ+27I9WSq233Me5ls4zKKXcLdsHLO9H2THG2Uqp32rdxx6W/c3+s26Q1loe5nYaA3AQuAwwAtuBWHvHZYntMBB0zr4XgUcsrx8B/tXMMSUDvYCUC8UEXA0sxTyG5nJgvZ3iewr4Wx3Hxlp+3u5AtOXfgaEZYgwDelle+wC/WmJxpPtYX4wOcy8t98Pb8toNWG+5P/OAGy373wP+ZHl9N/Ce5fWNwNxmuI/1xTgbuKGO45v9Z93QQ0oUNaqnFtFalwKfARPsHFNDJgAfWV5/BFzTnBfXWv8E5DQypgnAf7XZL4C/UirMDvHVZwLwmda6RGv9G3AA878Hm9Jap2utt1he52HuBRiOY93H+mKsT7PfS8v9yLdsulkeGhgKfGHZf+59rLq/XwBXKaXqGgjcHDHWp9l/1g2RRFGjrqlFGvoP0Zw0sFyZp1qfZtnXVmudDub/zECI3aKrUV9MjnRv77UU5T+sVV1n9/gs1R89Mf+l6ZD38ZwYwYHu2wfjRgAABHlJREFUpVLKoJTaBmQCKzCXZE5rc3f8c+OojtHyfi4Q2Nwxaq2r7uNzlvv4mlLK/dwY64i/2UmiqHFJU4s0k4Fa617AaOAepVSyvQO6SI5yb98FOgI9gHTgFct+u8anlPIG5gMPaK3PNHRoHfuaJc46YnSoe6m1rtBa98A800NfIKaBOBwiRqVUPPAo0A3oAwRgXlLBbjHWRxJFDYedWkRrfdzynAl8hfk/womqoqjlOdN+EVarLyaHuLda6xOW/6yVwPvUVInYLT6llBvmX8Cfaq2/tOx2qPtYV4yOeC8tcZ0GVmGu1/dXSlXNZ1c7juoYLe/70fhqSmvGOMpStae11iXALBzkPp5LEkUNh5xaRJlX9/Opeg2MAFIwx3ab5bDbgAX2ifAs9cW0EPidpSfH5UBuVdVKczqnjvdazPexKr4bLb1hooHOwIZmiEcBHwB7tNav1nrLYe5jfTE60r1USgUrpfwtrz2AYZjbUn4AbrAcdu59rLq/NwDfa0sLcjPHuLfWHwQKcxtK7fto9/8z1ezZku5oD8w9DX7FXL/5d3vHY4npMsy9SLYDu6riwlynuhLYb3kOaOa45mCucij7//buJ8SmMIzj+PeJQikLLOwmRYqSQk1GzUKKrGym2LFAoWwkynpqVhY2VophgVgaC43JaDI1xjWKmqUFKZI/JabH4nmuOW5z3xnTnTuz+H3qNnPP3/ee6c5z33PO/b3Ep5/jzdpEdKOv5nF9BexcpPbdyP3XiDfihsryl7J9b4EDbTqGXcTphBowno+DS+w4NmvjkjmWxOBmL7ItE8DlnL6RKFKTwB1gRU5fmc8nc/7GRWzj4zyOE8BNpu+MavvfuvRQhIeIiBTp1JOIiBSpUIiISJEKhYiIFKlQiIhIkQqFiIgUqVCIFJjZs/zZYWZHWrztizPtS2Sp0e2xInNgZt1EWuqh/1hnmbtPFeZ/c/fVrWifyEJSj0KkwMzqiZ+9wN4cM+BcBrz1mdloBrqdyOW7LcZvuEV8UQoze5CBjq/roY5m1gusyu31V/eV38btM7MJi3FIeirbHjSzu2b2xsz6Fzr1VARg+eyLiAgxLsTfHkX+w//i7rsy8XPYzB7lsruBbR4x2wDH3P1TRjeMmtk9d79gZqc9QuIaHSbC9rYD63KdoZy3A9hK5P4MA3uAp61/uSLT1KMQmZ/9RBbPOBG7vZbINQJ4XikSAGfN7CUwQgS9baKsC7jtEbr3AXhCpIvWt/3OI4xvHOhoyasRKVCPQmR+DDjj7gP/TIxrGd8bnu8DOt39h5kNEllDs227mZ+V36fQe1jaQD0Kkbn5SgwFWjcAnMoIbsxsc6b7NloDfM4isYWIv677VV+/wRDQk9dB1hPDui54mq1IM/o0IjI3NeB3nkK6DlwhTvuM5QXlj8w8HO1D4KSZ1Yg01ZHKvGtAzczG3P1oZfp9oJNIDHbgvLu/z0Ij0na6PVZERIp06klERIpUKEREpEiFQkREilQoRESkSIVCRESKVChERKRIhUJERIr+AHyMZjDVoo60AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Minimize the logistic loss using Nesterov's method\n",
    "f = lambda w: logreg_objective(w,X,y)\n",
    "grad = lambda w: logreg_objective_grad(w,X,y)\n",
    "\n",
    "w, res = grad_descent_nesterov(f,grad,w0)\n",
    "assert res[-1]/res[0]<1e-4, \"ERROR:  Your Nesterov routine did not the minimize the function\"\n",
    "print(\"Solver terminated in %d steps\"%len(res))\n",
    "print(\"Terrific!  Your Nesterov routine works!\")\n",
    "\n",
    "_, res_g = grad_descent(f,grad,w0)\n",
    "_, res_b = grad_descent_bb(f,grad,w0)\n",
    "_, res_n = grad_descent_nesterov(f,grad,w0)\n",
    "n_g = list(range(len(res_g)))\n",
    "n_b = list(range(len(res_b)))\n",
    "n_n = list(range(len(res_n)))\n",
    "\n",
    "plt.semilogy(n_g,res_g,n_b,res_b,n_n,res_n)\n",
    "plt.legend(('gradient','bb','nesterov'))\n",
    "plt.xlabel('iteration')\n",
    "plt.ylabel('residual')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 2:  Image denoising\n",
    "Consider this noisy test image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here is a noisy image...\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD6CAYAAABnLjEDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deZBlV33fv+ftS7/3eu9ZejbNaCQNSGgZBBKLZAwYg42wA46JK4VdKuNKJVU4OLFxnKTKVakUrlQFKJPYVsBlJXFZYHCCLEMwJpIxAYQEGkmjGaTZ157p6f3t68kf06Pu7+93Znqkkd6MdH+fqqnpc98555577j3vvt/v/BbnvYdhGK9/Yld7AIZh9Adb7IYREWyxG0ZEsMVuGBHBFrthRARb7IYREa5osTvn3uece945d9A596lXalCGYbzyuJe7z+6ciwN4AcB7AJwE8ASAj3rv912sTSqR89lk6cWyT8Z1v+0ulX1CfB85p9p4WaWjr6mX5HbuMi7bx7lNrNW75OcA4LrcsbyediEZOI/oQ4wtXteDjbU6a/Ybb3K7WJPb9LIJPZaYuOYmj7+X1vdMzksvxTfE62lSr5l4o6frrHGPXC9QQT7PotzN63mC6EfNQeB5UvdenlZPE5y4xNC8yDrqGkNzuYpGbR7tVjVYS9/ty+dOAAe994cBwDn3EID7AFx0sWeTJdy19VdfLLc2FFWd1NkKldsjeSqHFlgnxzObnmmoOvX1WSrH2mtPYqvA/eZPcr+toZRqk1pqUzl5epHK0++c0OcZlF8q/Pnw/qZqkz08S+Uz716v6gwe4o4yh85RufoGPZZOhlfhwFG+H0s7CqpN4UiV+92cE33qye2m+NjgC3VVR31xioUba/CXFwC4Js8/2lxnYbe+5oT4omlneQ6ys/o8zRI/G7GO/Fz/aE5V+DyheUnW+BrjdfGCSVx6te/5+89d9LMr+Rm/EcCJVeWTy8cMw7gGuZI3e+grRv3ecc59HMDHASCT0G9ywzD6w5Us9pMANq0qTwI4LSt57x8A8AAADAxv8ou3jLz4WeFwRVbH1E+NUnlsT43Ki9szqk16kX/qnLtjQNWJiV93oz9aonJ1q25TH+EfPrkzXE7Pid/b0D+zfJZ/6heP6TaJnwi5OMnnkWIKALTXlahcOtpWdarrWD5tDvJP/fSC/nkqv67lz/bctB5/Y5xFpPxxvmexhh7b/C2DVG4XQzoHvq+JKvdz+h7uAwBKR/mapPwtRTEA6ObEMhBz0E3rH8DJKo8tWebzxjr6eqRY0izo5Zdp8LOQWuD57mX0s9AqrvQT1I9cGNPFP1qTJwBc75zb5pxLAfhlAA9fQX+GYbyKvOw3u/e+45z7FwC+CSAO4E+998+9YiMzDOMV5Up+xsN7/3UAX3+FxmIYxquIWdAZRkS4ojf7SyXW7pGCpzWSVXWGDrBCIl7hPeZUJa3ayP3KwsmuqiOVI0s7WfGUXtRtxvbwuRNzvJ987q2sTASA0SfmqVy5nhVpwf38Et8GuQedmddj62a4Tbyp6+TOij1bofAKKZ6kcqp4mK+5vJX30AFtCFLZyrYRyqYBQOkgK/GkIu18v9yuJ4ywhp/Xir9EVSjohF1GZZNW8IaMZlaTrAQMfsRwa+tYEZsIGEIl6nyPisd1v40RVuyFFH2S+vDKvITsUF7sa82eDMN4XWCL3TAigi12w4gIfZXZu+kYljavyEwhWak9wDJHdoCt7lp5/f0UD8iEqo6QV4tTLDM21mlZtLyZ5btBoT/Izmk5eWkXG3oM/pDtjKo3advseEt6SKgqiubQ2rcuWeHxTd/BOpKxp7WBCWJ88vo6YTBzSreRPgKL21jOzMxp2TQ+U+Y2t42rOhlhk14fE/3Oapl9YQffs/QSnzveWntypdNUuxDwahF2+tLYaHGr1g30hC4mFTRqEk45QgSXTkYA0MmtVJJOYauxN7thRARb7IYREWyxG0ZE6KvMDoD2J+tj+rumeJzlmMw5loVcT/uQJ8vCkSTg8yuPxRrcb3NQ+2ln5ngs5R2sPwgFwMjMsBzZ3MZ78cmKltPKW9h2YPjHc1T2KX2bYm1uM3W3tlmY/DbvkQ8dELJoXvfbTQudiZBFy1v0eYpHWP8xKq6xm9T3efaudVRO1rRcv7iN7/XwPj7P3Bu0nkXub2fPCjuNo+zTDwDl3eyZ3c7xeGvr9fNUPMrjrY/yWJOBfXYpT8e6+pqlc1hlGztoFV7g+AgA4Horz+6ldBL2ZjeMiGCL3TAigi12w4gIttgNIyL0VUEXb/RQOrCiZClv08oeGfCwJSKYpBa1gks6koScKmTgvs4gK3dyZ7WBRifPxhQy+GKITpEVZ7UJVtzkpwLBI8/xNdW2sfNMKPKqNPQY26MjyDRHhGGH0N3EulqZk57l8ckoLiGHocXtPJcyaGLh2WnVxvWGqSzvIQDkp/lc8nmJ60tGJyvuvXAMae7QRk2FZ3h8yW0jVE5V9NgGjrBRUG2SFWnNkn4G00s8L40RrWyOi/lOCMWl6+hnYXXwy545whiGYYvdMCKCLXbDiAh9ldl9zKGzypAjFCR/eB8bgrSLLNeUN+ngFUP7WZY+8W5tIJOf4u+1zDQbaFQ2a/2BNIJovIGNaorHtFOIjAwrDXM6gSwsrSLL36klllWlcQkAZIVzSW1CBznIzgj9hhDRpcMHALo/AJCssC6jm9HznxKRVuMNmdIm4EgiGDikjUVmbxuisnScGn5Gt5FUtvGzkKhrmXf+zWzgMyAi0CYCRkEySET2LCe58An9PFXHeR5yM1r/0RjmOvIe9ZI6om7hyMq5pcPXauzNbhgRwRa7YUQEW+yGERFssRtGROiv11sM6K2KaCqVcQBQn2BDkI7IqBmKtFrZwkYdWx6eVXVclRUo7Y1s1FHav6DaNDawcifWllk411Y8Sarr9JQXj7EhizTECUUqlcZF3ul+ZTbb+V1s+JENKIjaIiVRZZLHEoqEklkQHmAibdbitjHVJjctI91qperAKbaaSc/wPTy3W0TuBTD2A47um5viNpVNWnEmvRdbg6zsDBk1tYb4OZXzVnhBP08De4UVUForXpduYgVcXETaCXlNyuy2F8Pe7IYREWyxG0ZEsMVuGBGhrzK763hytOjk1852kV5guTLe0HJmUxjeqCwsANpZNtAY2iPk+s7aBg7xJn83yiguAJA6yVFmZu/mNMm9wCXLLCDxJstgscBXsszmEnQQGmIDmIyYy1B039QiG9H4mIjAUtXya3kTP0aZea4TYzEaAFDcJ+Zpt86ukxCyclxE55GGRQBQ38Kyv5S3m4N6Mof3soFVL833vTWol4mMfJQTDk7eaYOl+g18jd1ApFiZ7rp8HWfXSc3raybnmEvI7/ZmN4yIYIvdMCKCLXbDiAj93Wf3Hq69Sm70AecNsS/azXEduQcNAAXhkFLboLNxyAio3RLvt8r9fUDLzjJbaayt5fy5u1hGz84KOXlKy1xdYUugMsIE4hGk5llf4APyXzvJt7ejsrbqsVTE3I1/8xiVazdzJFYAKJxgfYG0JZDBLACgtpX3k0eenFF1mhvY8aibXtuuQe5DO5FtZzQQ5KM+wc+CF6tCBj4BgGbA8Wg1rYm8OlYbE5l3A5Fg6+t5LK28iAi8U9sJFI6vXJOPX/z9bW92w4gIttgNIyLYYjeMiLDmYnfO/alzbto5t3fVsWHn3LeccweW/x+6VB+GYVx9LkdB92cAPg/gv6869ikA3/bef9o596nl8u+s2VPMoZdZUWzI6K0A0BgRKX2EcULIqKMpUu+k53SkWBUBVTiXyCiegDbIqGwUkW6XtJImLR11hJFDapYVhQBQ38gOKuk54RizQStlYi2RznhEK4RaBf4uL73AEVGXdoQi+ggHmwGet5AjT1sokUb3spK1F1AaVTZKZx89luQSK9NqG3keXCA6rnROamwUhjjT+tmQ/UjHmMawvmYZEVimUo4FosAO7eOISq1hrRTOTHGdXoKVlCFl7ep15C+hw1zzze69/w6AOXH4PgAPLv/9IIAPrdWPYRhXl5crs09476cAYPn/8YtVdM593Dn3pHPuyXZbu7QahtEfXnUFnff+Ae/9bu/97mRS/9Q0DKM/vFyjmrPOufXe+ynn3HoAOuVHAB9zaBdW5NxYKyAnC68PWac+GsgcIrO5BDLCDD/Fkkj5RtYptgf09172DMvOQ88L45ENOtKqTD0snR1aUicBoCNSBHfTWkaXtEbWrpOQUV7F3Iai+8Y6XGfmLg48IQNtAEBTBHuQ6YtDKbSHn2LvmOr2oqrTybIAmpnmc4cNoXj+C0eFkVZGPz9y/nPivqfmtCHO2Tv5xTV4iK8xM6MjDzfHeLwhY53GOu534Bj/Gl7azvodAEis0mNJw6/VvNw3+8MAPrb898cAfO1l9mMYRp+4nK23vwDwfQA3OOdOOufuB/BpAO9xzh0A8J7lsmEY1zBr/oz33n/0Ih/99Cs8FsMwXkX6G7yi62nvdLX8foGBg5zlY+5Wlq3zU1p+ksEF8ifrus4EyzpSjiwc4D1oAJi5neXI0lE+d+Ggzuq6eKPYLxbianpB7/PK8TZHWBcgM6sCQEPIfyH5WzpayDaDB/Sef2eA70lxP8vWnUGtK8icE3oKsdedmtfyq3IUCewfJ+psszB9B7eZeFzfs9lb+D4P/YTntrpeP3PSTqCT4XkaeloHj+yleSy9FPcRyiS8OtsqADS26LGMPc33xAl7ilB22PyplWfK9Sx4hWFEHlvshhERbLEbRkSwxW4YEaG/kWoEIQOZzBn+/pGpcaXSBtBpkstbteFKssJKJJkWubZFW/cNHWTFUnU9K85C0UGlUix3mvsIRbepTfJ480dZ8bd0fcBhRaQVzgb6bY6KCCxCt9Ma0lF/lDKwzQqiuTfouR17nJV43RwryWJLWmHaE2NrDGoPDhlZRyrkeoH016XDrESNN3leiof1WCpbhIONMEypbdMGP+M/4nmSUY9DEWNibX42Cqd1ROCaeMbi4h6ll7QCrjG20qaXsEg1hhF5bLEbRkSwxW4YEaGvMns3HUN5lXwUysjaFMYWpYNsZBAyVpDOD+n5QKaTeZaxytv4POnFQEbTAZ4emcFUGlIAOjNLa5BlLpkJFgCaRZZXc0JPUdpzTrVZvJUdVKROAgC8MBxa2sLXM/QTbawzcwvPS2ae5dn8GT1P5Z2cgSd3Sug6btDZXmQ2mtURUi/QzfCE10UQj7h09AHQGOK5HJxmR5LOBu1IUjzEdWbfyHUKR7ScDxmcQtyz+Jw2uGqMTVBZOvYAQH2dzGLM/cqxAkCrtPKMhQJ6XMDe7IYREWyxG0ZEsMVuGBGhrzJ7rOORmVuRaaVMCQAtIb9Kh5X6SGhvVTg7BDLCLG5jp4PmMPe75as6I0lnmOXXWJPP3SrpsaisMUKG+oXP/Z1q8/dz11P5lyaeoPKRpo76Ve7yNWZi2sFmX5mz0/zC4BEq/82Zm1WbX57YR+VSnOf272ZvUm1uKpyhcqXDe8VLHe08c+/gT6j8/aUdqs7nNz5O5T+c30LliaR2UPn8v/7HVO4WeSwyiCWgM81kpS4pEDyyvEMECZ0XDivbtN1GQgQ5VZl/oPVPcRn8clw/26uDa5gjjGEYttgNIyrYYjeMiGCL3TAiQl8VdD7G6YmlAQqgU+7KrMIhBd3CDlYApctaoSINeHpTrKCbv10bfuRE9pBElceWPqcjvSzsYgOTuHB+OFDXyrae8FD55hwrzt5SOqzatEXqj+mWdpYZTLFy7YXaOiq/e5yVZAAQF1qjvzx9B5X/0fofqzYPHHw7ld+3aT+VpTIOAI41eb6zca04u/8493vfCJ/7C6ffqdrEhCNSu8hGTS6gbEvNshFQZyMrwSqBiK71UX5PxtvCmOf7J1Wb+bdvonKroBWX2Vl+xqSzVSjjS3XTSj+9/eYIYxiRxxa7YUQEW+yGERH6LLM7tFfJ7NkzWmbvJVkoqa1jmSs7o9s0RriNDBJwvl+Wiwf2snPJ0q1alm4PcL/tgjhPSwd/kI4L0sGmExC6dhXZKOWW3HEqP1G5TrV5X+kZKn+5caeqsy69ROXRBDtnHG2MqDbNHhsfSX3CptSsavPR656k8kyb9QfP1TaqNkNJdujYntFJhZY6m6l8uMX36NRSIKjEKQ5wsXjTIJVTOiCtClySnZG6Gm2w1ElzEI/saWHYdcsG1aYoIhhXNwd0AUInVTzCeqFeWj8/bpVzVSjL0oufXfQTwzBeV9hiN4yIYIvdMCJCfx1huh6pVXvgoeB4DZEBVDoYyD1RAOhmWPaRQSYAIH+E5aWlN7H8JwMuAEC7yPKr3GeXwQEBYPCQ7mc195b2q2Nn2yxXtj3flvEUy94AcKbDbQaTOsDCeuEocrzJMvrNeb0XvL/OsuYnt/4tlR+eu021+egIO6w0MjxvX5h6h2pzz/AB7vfsm1SdD0w8y3WmuM7t46dUmwObdlE50RCBRhf0fn6sKfa2c/wMhoJHysCinQFu08nqNr2UWG6BLDhpEQh16m5+tif/Rus25t68YrPQtX12wzBssRtGRLDFbhgRwRa7YUSE/maE8RxddXUmiwvIqKMdYdgSb2tDlsKzrLSob9fGIrVNbOghFSyxgIOEjGTbFMrD1oD+rmzn2IlCXk+I481hKpfibHDy9BI7UADAU54NTm4rHVd1pELug6WnqPzpE+9Xbd42cpDK35hnpdg7Sy+oNn8x+xYqj6XYeGfngFYqyUg7/2nbV1WdP5lhR5dfm/x/VP76nI60I+9ZaomVb6FnLj3LbcqbuU5mTkfUlYpjmTY8dN8b49xvo6Sfn1SVxzL5yFkqV3bpZzt/esXoJ2RQ9uJnF/3EMIzXFbbYDSMirLnYnXObnHOPOuf2O+eec859Yvn4sHPuW865A8v/D736wzUM4+VyOTJ7B8Bvee9/7JwrAPiRc+5bAH4VwLe99592zn0KwKcA/M6lOvIJh+bQyimlYwAALNzIsvXI96aofPoD2qmilOfLyJ4MeDs4loUW3sBBJlILeiriwqlABttILuk2s29kuWz8CXZkON3W34nvKe6l8sPzt1M5ZDDz3kE2OPlBRUdnvTF7mspH2xww4ufHn1ZtZJTaxQ47fPz1jDZ+ua14gsqH63yekaTOYrIhxZlff//kz6k67xrmoBffWbqBym8pcbRcADh3jI2CKjtKqo6kPiGz866dObg+xoZD8vlKV3Sb5AJngEmWtf5JOmzVr2N9TuacNgqqr9IFyKzHq1nzze69n/Le/3j57zKA/QA2ArgPwIPL1R4E8KG1+jIM4+rxkmR259xWALcBeBzAhPd+Cjj/hQBA+4gahnHNcNmL3Tk3AOCrAH7Te6+NtS/e7uPOuSedc0+2G/rnnGEY/eGyFrtzLonzC/3Pvfd/tXz4rHNu/fLn6wHozVQA3vsHvPe7vfe7kxmdJcMwjP6wpoLOOecAfBHAfu/9f1710cMAPgbg08v/f22tvuLlFoa+u6LMaW0dU3WSIkXO4u0cEXXdY3Oqzdl3sNIr1tERQGrjQok3wwqU5Kz+1TF/KytHhp7iHzSnfkGPf+xpVnB54dl3osF9AsD35rdT+Q5hILMYSJ/0P8/eReVPbvhbVee/nv0pKv/iKEdn3VufVG0OV1m5NtvgL+i3jmql2NNL3M+mHCvfmj39mGUcz9O7R7Q3oIx4c2iJx/ZLwz9Ubb68maPXDBxcpHI3r41qpHdj7vSlPRcBICuzaAu9WKypFXTVzTl1TJKoslJYpQX3ei6Le1YiHcVrOqrOi32veXbgbQD+KYBnnXN7lo/9G5xf5F92zt0P4DiAj1xGX4ZhXCXWXOze++8i6HkLAPjpV3Y4hmG8WpgFnWFEhP5Gl00nSE7vZgLpLYQdf7LGMkxlpzaSyJ0TTiyBtLXFQ2yYsrSd5eD5nTojTOEEy12tdSxDrvsBG0kAQG2C5T8vjC125djQBQDeUWTnkr9deCOVY4Hcvj8z8hzXkbmioR1SHprmCLQ/P6qNaqR8nY6zIdHpBkfIAYD1GZaL78xzBpvPH2PdAQCca7FeZaqm7+s9YxzN5pYhjkzzwJl7VJuFHTz+3Ame/05BGMMASJZZzpVOUY1RnSZZ0hziZzkt0z4DqKznOrlpfc+kIVf+eY7mO30vp+EGAJ9Y2fXuzV58Sdub3TAigi12w4gIttgNIyL0VWbvJRwFgOikAxE4hUiVqLG8mpnVjgDxGsuVnbyWy+Jllq8Lx/jcyZp2SpCBAKoiO43M1gEAWREIIX2W9+/jAdl6X52de04L+XXrgM7C8uj8jVSulfT+8Vyb98jzIlPqwcaEahNzfM2zTe6j09P37CObeL/7qfpWKn/u+i+pNl+Y4YizO8e0TZbUVTx+jvv9yKTOKHv2Wc6e0x5i3UwoY4oMaCGdUZoFfc1J8VwmhH1I5ox+NiZmuZ/yNm1k1smyXD97N9uZxJuB4BRrx0cBYG92w4gMttgNIyLYYjeMiGCL3TAiQn/TP3U80rMrBgzddVoplhNpnGNd1j7I6KEA0M3wZUglBwDESmwYUR/nc4eigSZFuicfT4hyQMGYEBFxdrGy7R1ZNjgBgE8v/gyV/93mv6byVxd2qzbvHBKGODO7VJ3tAzNUPlVng5hCXDt87MhwNNNjNXbcec/wPtXmsyfeQ+UPrdtD5f8xx047AHB9lhVycx2trCr3+B69ffwQlZNOp+9OLgkFbo8VctN36jTP0igrJVIwOW0fo1OAH2LDoup1+jytPD8v6SWtLEwJA5/MNF9jc0wb+JDO9xLKOnuzG0ZEsMVuGBHBFrthRIS+yuyu06MgEcmiPn1jlA1iMnMss7iAbN0W2Ti66YBcn2XDCSmjp+cCTv/iqzA7xc40vZTWDXQz3Gj4Bxwd90dNHR13vsWGH882ORiEjBILAGc7rAu4taTTL48nOdiGdHL57hwHzQCADVlusyXHwUJCRkH3b/wulY+12KloNMkOOYCWt+c7OrDDxjQHwXi2zHMXmpdujq+xLSIPj/1Ij6W+gedfpuauj+jnNCGMWzpFlqWTS4HosmU+VhvXxl+9JOspEjU+d2pBP6ftgHNPCHuzG0ZEsMVuGBHBFrthRIS+yuydfAKzd6xkoYwHMk42iyxvZ2a4jsrEASAjMmqm5vX+a3WSZaqWcG5I1LUsWp7kczWHeGylw1ou6wh9Qe0GDqe/v65l9ttLnFHlf5+9jcqf2aoznH55ieucaOhMM0dqLDu/f/gZKn+9fYtqc2eB97JlhphvBjKnHi5zZtGbBnmvfqmt94ZvL3JQzZtzWudQ67GeZSLNmX7+5ChneQWA4qIIKCLtNBL6/ZY9zbqYylbe88/O6udpaRM/G4lqQpR1m16adTwyMAsA5I6z41RjPesyWoP6+W/nV/q9oowwhmG8PrDFbhgRwRa7YUQEW+yGERH67giTWVhRamXP6FTECzewciTeYCVYL6EVFI0hcRmyDCC1yP0kGqzICDnYpMoi0udZLq9WjKz0KyKYVFjB9ejZnarNm0ePUfltIwep/Hsnf161edsQ16l39bzcmGdFmTSIubmglWJTIqV0TTijbM/JVCjARJoNcc422Qlke1632ZBkg5kTbZ0pp9Jlxd7NOVZkfuQGnRHm3+Z/ncr1dazkGziijWpmb+Hxpio8TyFDrsyieBYK/CyEFHTS4Eoq4wCgJ5y6ZJvWgH4/F06sKCVDDl0vfnbRTwzDeF1hi90wIoItdsOICH12hPFIz6wyeuhpowIp80JkVCkc0wEXnMgAE6treamxjo0T0vMsS1fX6eisA8dZp7Cwk/sIJCfF0EFh1CGMHP7JJi1nHqhzlNfH57dROZTh9KnKZioXEzo7zb4KZw95Y5Zl3rbXOoe0MKKZSHJQhj1VPi8AvK/0LJX/T48Nb861OJMOAKxP8lzWunr+D1U5S+5QgmXco22dxUc6hTRK/D7LFPV5Bg/zM1We5DrV9Vofkp0RgU2EzkdGrAWA1KIIRDGhnX9k8JN4k9dILKN1S6vH233ajGoMI/LYYjeMiGCL3TAiQn8zwiQdGuMr8kV6QX/XdFPCkUQEpXSBbcTis5wxpXq93rNtCtmtslFcekDUkZleZcCC/CGtP5CymszcubeqHWFkZteqkF8fmdYOKz8zxoEff1zWsnSjy9f4g8oOKh+tsQMLAAynOJPJewdZHp9pcvZVADjQ5KwlPc+TuS2r99m/M8/2Bu8VWWkBHRAzKSI/fq98vWoj1RDZOZZ5Q8FIM2f5mmNi/tMLgYysG1mOLxznQJe9UAAVuWde0GPJn2HdS0847rhAm8yqQBm2z24Yhi12w4gKttgNIyKsudidcxnn3A+dc087555zzv3+8vFtzrnHnXMHnHNfcs7pjA+GYVwzXI6CrgngXd77inMuCeC7zrlvAPgkgM947x9yzv0xgPsB/NGlOoq1PUVoPXe7VvYkhG+MTJucPaONRzojIptIQNkmnRtiXf6eCxnIjHzvDJXLN3PUGWnwAwCJGitzBg9x+e1FzuQCACdbrCibEWmS/+PW/6XaPLx0K5V/VkShAYDHFm+i8vdn2FhnJKMdMe4qsoPN8002zNmYXVBtpHLwt9d9k8q/8fyvqDaTA9zPD8vXqTpNoWDsiRt738hTqs0Lz2+i8vyb2WAppGyLVfiZyswJxWxZR3T1k8JhJaD4k8jAvKlKINJRlvutTnB59Ada2Tl754rxkTTKWc2ab3Z/nguuQsnlfx7AuwB8Zfn4gwA+tFZfhmFcPS5LZnfOxZ1zewBMA/gWgEMAFrz3F+z/TgLQe0rn237cOfekc+7Jdlu/SQzD6A+Xtdi9913v/a0AJgHcCeCmULWLtH3Ae7/be787mdTJ+wzD6A8vyajGe7/gnHsMwFsBDDrnEstv90kAOj2HbB93aA2vBCQYeU4bpcivn/oY6/1C0UHnt7MRxPBPdFCMxDwbTvSybBQhAxgAQGsDZz2V+gPX1Y48nZzQUwoR6pHZN6k2p6p8njcMchaZWk87Yjw6zUYp6zdqWVpGdf3FDZxdteH17T/XYaeV6RbPy615jgoLAEfi7LDy0OKbqSyDcwDAXJu/+E/WBlWdsRxebUAAABRaSURBVDQHmtiY4Ws83BQ6FAD1HewcIx2rOpmAbL2Br7E+yvNSOhd4nuoiaq3oNndC/4pd2sFzWzii63TzfK9TFX7eF27Tzj+ro99ekVGNc27MOTe4/HcWwLsB7AfwKIAPL1f7GICvrdWXYRhXj8t5s68H8KBzLo7zXw5f9t4/4pzbB+Ah59x/APAUgC++iuM0DOMKWXOxe++fAXBb4PhhnJffDcN4DWAWdIYREfru9VZdt3JKF7BkGXqOI5XGReSR+rg21Bs8xEYR8zdkVZ1EnZVVaZFSN6aD26A5cmllWygibbLK/Xay/H16z6A2qpnKs3LqUI0VXt8oa6+324c56ky5p1MsvX3wAJW/t8gpmmX0WQDoOWFsJDzY5jraEOqjpSep/Icz91L5aEV7131yExvenG7r9FUtofVqeL4fj83doNp0RSSXzDl+NtolrexsidTh0vCmO6CfOfmcVrax8s01tcGMTHe2cIOey9WRYs93pKooljavjF96ja7G3uyGERFssRtGRLDFbhgRoa8yu48BnVXi9MhebVRTvo7lmPoIfx+FDHGWtrK8mlnQxi7tLMsy6VmWjaQDCwCk9nHGlO4WdqqobNHRQYvPczRWv56v53hLy6/fmeYIMr+x+TtU/tb8G1Sb24tsqLK/ukHV2ZFjmbzWYdlTpmMGgMdm2Vjnc1u/QuWvlLVR0L88+mEq3zbI+oSO1++UL8++hcq/NvoPqs5zTb4mmZ3mjQVtx3WwdCOV2zl+Nkov6IwwvbRItywMsKbv1vcsNyOMX5ZY6dMLyPkyumyiquXrWJuf3aRw4AoEBIaPrRx0+jFe6fviHxmG8XrCFrthRARb7IYREfqbxbXtMTC1IlSUt+isGXK/O1XmvcnmsN4nlSJh7rR2XIDIGtMuskwV6rfxdg72cDn77PVJ3m9tD7CQtT6pHVbu2/A0lb+7xHLzjQPsGAMAGcfy9j8be0zVeUTsz39wjB1hznZKqs0dg+zo8t/m7qZySUYXATCRKVO50uH7ujU/p9qsS7Nu46H5t6g6co9/a4ajCJ9paeel0e9P81h2sbzdDmSESU+JPfOdvOffHghk+D0k5O+yjC6rl5bcz89O60Asa5Gab6ljPrZyTbGuRZc1jMhji90wIoItdsOICLbYDSMi9FVBJyke1sqe+gQbQcSkEUHAzl8qJVqDWglT2cCXmjvH1gfpBe0JI6OKdkXkz/md2nCiNs6Kvtw0G0X81a+/R7Wpj/N4C8+ykmnPrRxJFgASNe73zwMRSqQTjjS4SC1qo5r2gIyayn3E69pgKb+fjXeevofDEWbntKXHXqEwzT6q0z813snGRN9Zz2MbPKCfn+ZNfE9iLT5PbUIrYtNn+Ro7Ik1Tfkpfcy956fdkdYN+BmVkmvI2HaYtPc/PYW1cRM1Z0gq61MLKfXSmoDMMwxa7YUQEW+yGERH6KrO7LpBYlQUjXtFGBb1JkY1DyIiJqpb/pHGL62gZK1UR6ZZ/MkPl6o06ame8LmReIWcWTuqxSAebxevZWaaT11PezvH42xs5mEVjUH8np8U1l/Zqw5XGJBud1Mb43AvbdZAPqWMYemaeypXrtSFOawsbruTPsC6gk9PeG9lplj1bb9XRyVX65XmebxVcBNqgqpvha5ZZWQCgfAPPd/Ycj78xou9Zc5CPtQdY/g6dR0aOzU3p5781yNc0vIfnvzsQ0EdNrhzr7r34+9ve7IYREWyxG0ZEsMVuGBGhv8ErEg6t0sopZRBIQDv4S0eSdk7vk6ZE8Mh4WwtMMlNGY9swlfOH2TEDABoi8ETmFDtMNEa5DwBojrBMlRF7zK2ill8Hf8KOJF2RVSYo/wlxtb5Zy9LCjwSpKneUndU6h0SV579b4Otp5/T7IXtGBO8UWXvqw7pNTATSCGUfzcywXF/ZxGNJ1gL73ymeX9kmFNwhO8MyemWjkJuf1c9GV2QUKm9j/UdmVtttNIe4Tf6oDqQh9ULTd7FTTmZOX/NqvUqiFoicuoy92Q0jIthiN4yIYIvdMCKCLXbDiAh9VdDF6x0U968oE1xdGxXEBlkplmiwsqRV1Aq6RF1k8Ag4KeSPsxOCNMSpXqcVXLkjrJCbu4ONR1JlrSzJHeU29c0cuaa4L5BaeRcrYUqPn6JyYnJStZHZRaQyDoB2GpI+El47TTRGeb7jTb5GmS0FADpCaSqdgcb/gR17AKB+HSs32/nAe0cckpFhYzXtFIIWK9tyaVbYpU+wkQoA9IqsXCsJJVd9g3ZYkcSbMp23ntuB5/ncnSEdnbi6kZXWI89ypNvWkDYkaq5bWTO9Y4Hws8vYm90wIoItdsOICLbYDSMi9DeLazqO2rYV2biX1IJm9gxnfHEt4QgTyNwiHV+WdmgHj9Jh0UYYL8gyALTGWVbLnhMGC06Pv3wjy/7NggiMkNPZSivrWc7q3c0yen5Ky6arjZMAoD6qb2XxCDuFtEVG3JCDSrIsoqYK+XVpm57bwnGuI51nzv7UuGozcIrvY2mv1mVUdgo9inB8SS1qp5CuCDwhdTN+mzaEahV4HrppbjMYGtsOHpt07AlFK65dxw43MkgGABQPsIFVdSvrsPJH+HMAqG9YVecSWV/tzW4YEcEWu2FEhMte7M65uHPuKefcI8vlbc65x51zB5xzX3LO6T0BwzCuGV6KzP4JAPsBXIiI8AcAPuO9f8g598cA7gfwR5fqwLV7yJxZ2TesTer9SxlwMt5geVzKZAAwcJr3tocOaMElMcvy6+xulp3TS3rPPFlhWbQxyjJiek4HbMzKwA1pvp6QI4kMGJFoysyd+nryh/ia27cMqjpSpyDtEWoT+vYnatymskk4eASCR0p9x+IuHsvIXh0YsryF56W+W+syQo46q4k39OexVsBraBUhWTouglImK9xvL6PbSKTeKOTYIwNrVDbrffYzd/PcDR4Wz1jANqK+KihJSA92gct6szvnJgF8AMAXlssOwLsAXMjn+yCAD11OX4ZhXB0u92f8ZwH8NoALX18jABa89xdefScBbAw1dM593Dn3pHPuyXanFqpiGEYfWHOxO+d+DsC09/5Hqw8HqgYDVnvvH/De7/be704m9M8WwzD6w+XI7G8D8EHn3PsBZHBeZv8sgEHnXGL57T4J4PSrN0zDMK6UNRe79/53AfwuADjn7gXwr7z3v+Kc+0sAHwbwEICPAfjammdzDr3Uyiml8wAANAfZwEFGmEkvaKXY6XePUTlR1/0mRMQS2e/AYW2sUN7BTiwya4x0EgkhU1BLIyEAqI/wD6zCSf48GcgC0hlixVkoasvMLcLBQyh7QvMvM51kRIaS1KxWtsnIOvmTbBjVGtYbNSmhEE1WdYQVafQjs9HM79QGPvlp7kdmxUlW9DzJlN/SkSfk1JKd4mtsCwct2ScALOxkhXTpQFXVSS9Ixy9hPCVSgp/vZ0U8lgrt1VzJPvvvAPikc+4gzsvwX7yCvgzDeJV5Seay3vvHADy2/PdhAHe+8kMyDOPVwCzoDCMi9NURppuOURTO7LSWv9PzwsBBZFI9d7s2xMnOiKip53S/zSG+VBl4QkYlBYDCATZcqW1ieambDmRqmeGAHPkzIkPMVh1RNyfG3yxJ5xl9zTKLa8ipYuD0pY1oEo2AbkM4GkkjDRlVFQCaI3ws1paBHFQTwHGdkFOOvCYvjIQKp0JZaLmfzKzQU1R1m8pWnt/CSb5n8XogUqyIIlwV2VaH92t5XGYXXtqu72vhGOtEemJuEzV9z2LtVRMcMLp5sd5FPzEM43WFLXbDiAi22A0jIvQ34GSzh+LBFVmmuklb1KXFvm5tA8u4pSNa5pJyZuqszrQR63AQgMok72emZ/V4WyM8vtyjz1G5fu8u1aYjAkRUNojMnTIABoDsMRGwYDtnX5UyMACk51iurI9q+a9ZZBk3vSgy2Z7RAT/bA/xIlCe5HG9qmV3K5M5L24hAVt0Fth2Q9xnQ8rbMnJqs6n47aa7THRc2AKd0m+ILrJuRASPiKf1OTM3x+KXjS3tAz5O0l4i1XvrykzosAJi9ZUWX1D1iAScNI/LYYjeMiGCL3TAigi12w4gIfVXQwQE+uaJASJVDkUb4WGOQFSwTT06pNufuYVf62AbteBETp5JOIO2CblMVyrXuPUIhF/A5kEqkzDyfWBpfAEDuMNfJnWSDjFg9oJQR2WmKR7WyrSGNXYTzT3JGxxeobOSIMaNPi0w6cf1+aJX4PMklHm9tvY4CG+sIZ5NAdN/66KXrJBb0NceHZRtRIRARuJvne5+Z5n4T57STVGXXKJVbA8IQKmDklMyK7DQzevxVERlIRUHuBLIQTa88PyFl7oufXfQTwzBeV9hiN4yIYIvdMCJCX2V27xy6q7JqysACAFCeFNFM51lGmbtrg2ojo4MWD2knhOYIG22kv/EE97Fzu+63zsYVUp4KRZdNLV5cZgKAoee1nLbwJpkdlmX40Dxl59g4RxqynD/IRRn4o1PUhiwy+0xjlOvE2msbyEi5PiS/Slk6lIW2cITvY+IUWz5Vb9FhD2W/2bMcZKKyWQe8kHLu4jaWrcee0QYyMlKsm+B+O4Eowg0RmCVZ0QYwcRHcRDr/+ETgWVgVsdkF7s8F7M1uGBHBFrthRARb7IYREWyxG0ZE6K9RjSCUIqedF8dEmE5lJAFgeB97uUmlEqAVH+337qZybUArSxJVVpTJKKPlzdoQRyoU0/OsFJNpkwFg8FlOcSyVMt2Svp5uhsdbH9NjGTjE81LdxgrHUKT/zAwr22QKrOQCK7wAIHaGFWfd9WxwMrRfKzKrm9hLL7WovQHbKrLLZiqXDmpFbGyer7k7ytGFeoEnvniIjYvyp3j+Q89TfR0r5PIH+R6Wb9SpoXNL/Dx1Q950Yh7kfV3YqT1F82dW2oQUeBewN7thRARb7IYREWyxG0ZE6KvM3ks7LG1dkUFkhFcAGH9SOF4IuSYWiPRZ28jyk3RKAIAMi1SorBORVgPZUWT65fS0MKSQMjCARP3SWWJkxhUA6BZYNu3kWa4PRV7NnGPZWTrgAECswuNNVnie5m7Qcn6ifulHorVVX/NAk+epO8D9Jpa0nC8J6WJkVNeUiOYr5WYAiAnjKek4Mvi8lvMb49wmJ+TvziZ9HvnsVnayA5HSPQFo5/k+Fo7qeYk3Wa6Pi4wwrqf7zR5ffPFv6Ui2GnuzG0ZEsMVuGBHBFrthRIT+Rpet9zC8b0VmWtyh9wxlNsxWSQSDmF17b1JmlQGgAk0MLQl58OScarKwex2VZTTT/HEt/9U28jWlxb51KIuMzO7ZEdu6Qy8EAlOIvd/clJb/ekWWNeW+bv6Mlu+kQ0pNBNsonNQZZZ3QDbS3DnK5oOV8affQCugcZOCSoef5PD5gpxEXEVznb+ZIvd2s1n9kz3C/5V3smJSb1tccr/CxxZ28n5+ZDcjO4tZ3AhFoZXCNdl5k1Z3T/VauX5nv7imLLmsYkccWu2FEBFvshhERbLEbRkTor1FNMoba+hWlUSgtUFdENZFGKq1iIDrrlEhzm9ZKiqaIOhpvcL+L9+gIOO0cK0vGT/B5unmtYJHRa+rrWGGXrGqjoLRwkMiLVMshhyFJKALqzN0TVE4JBWM3pfvNi8iwKVFeCqSc9nFWZGam2bGklwoY/LR4HubeWFR18lNcJ9YUqcHGtIIXQ3xPisdYcVkf1YZEMoqRfDZCSMOhwf2cQmru5tKafQweCBj4CMVrUqTm9oHXc6sQW/X5xZ8Ve7MbRkSwxW4YEcEWu2FEBOdDUUlfrZM5dw7AMQCjAGb6duIr47U0VuC1Nd7X0liB18Z4t3jvx0If9HWxv3hS55703u9eu+bV57U0VuC1Nd7X0liB1954JfYz3jAigi12w4gIV2uxP3CVzvtyeC2NFXhtjfe1NFbgtTde4qrI7IZh9B/7GW8YEaGvi9059z7n3PPOuYPOuU/189yXg3PuT51z0865vauODTvnvuWcO7D8/9Cl+ugXzrlNzrlHnXP7nXPPOec+sXz8Wh1vxjn3Q+fc08vj/f3l49ucc48vj/dLzjltz3qVcM7FnXNPOeceWS5fs2O9HPq22J1zcQD/BcDPAtgF4KPOuV39Ov9l8mcA3ieOfQrAt7331wP49nL5WqAD4Le89zcBeCuAf748n9fqeJsA3uW9fxOAWwG8zzn3VgB/AOAzy+OdB3D/VRyj5BMA9q8qX8tjXZN+vtnvBHDQe3/Ye98C8BCA+/p4/jXx3n8HgAxZcx+AB5f/fhDAh/o6qIvgvZ/y3v94+e8yzj+UG3Htjtd77y+ka0ku//MA3gXgK8vHr5nxOucmAXwAwBeWyw7X6Fgvl34u9o0ATqwqn1w+dq0z4b2fAs4vMADjV3k8CufcVgC3AXgc1/B4l38W7wEwDeBbAA4BWPDeX3Bnu5aeic8C+G2sBDQbwbU71suin4s95HtnWwFXiHNuAMBXAfym935prfpXE+9913t/K4BJnP+ld1OoWn9HpXHO/RyAae/9j1YfDlS96mN9KfTTn/0kgE2rypMATvfx/C+Xs8659d77Kefcepx/K10TOOeSOL/Q/9x7/1fLh6/Z8V7Ae7/gnHsM53UNg865xPIb81p5Jt4G4IPOufcDyAAo4vyb/loc62XTzzf7EwCuX9ZopgD8MoCH+3j+l8vDAD62/PfHAHztKo7lRZZlyC8C2O+9/8+rPrpWxzvmnBtc/jsL4N04r2d4FMCHl6tdE+P13v+u937Se78V55/T/+u9/xVcg2N9SXjv+/YPwPsBvIDzstrv9fPclzm+vwAwBaCN879E7sd5We3bAA4s/z98tce5PNa34/zPyGcA7Fn+9/5reLy3AHhqebx7Afz75ePXAfghgIMA/hJA+mqPVYz7XgCPvBbGutY/s6AzjIhgFnSGERFssRtGRLDFbhgRwRa7YUQEW+yGERFssRtGRLDFbhgRwRa7YUSE/w+54nGSI92S0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Don't modify this block\n",
    "image = zeros((50,50))\n",
    "image[15:35,15:35]=1\n",
    "image = image+0.1*randn(50,50)\n",
    "print('Here is a noisy image...')\n",
    "plt.imshow(image)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Here's the TV loss function from a previous homework..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Don't modify this block\n",
    "mu = .1\n",
    "def h(z, eps=.01):\n",
    "    \"\"\"The hyperbolic approximation to L1\"\"\"\n",
    "    return sum(sqrt(z*z+eps*eps).ravel())\n",
    "def tv_denoise_objective(x,mu,b):\n",
    "    return mu*h(grad2d(x)) + 0.5*norm(x-b)**2\n",
    "def h_grad(z, eps=.01):\n",
    "    \"\"\"The gradient of h\"\"\"\n",
    "    return z/sqrt(z*z+eps*eps)\n",
    "def tv_denoise_grad(x,mu,b):\n",
    "    \"\"\"The gradient of the TV objective\"\"\"\n",
    "    return mu*divergence2d(h_grad(grad2d(x))) + x-b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use your BB solver to minimize the TV objective, and denoise the test image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver terminated in 111 steps\n",
      "Here's the denoised image...\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD6CAYAAABnLjEDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdcklEQVR4nO2dbahlZ3XH/+veuXduZiY2xsYwzYSaghT90GoZQor9IFFLGsWkYEGRMoVAvrQQUdGxhVJpP8Qv6ocWZTDiFMT4WhJSSwlpggiSOJpojYNODLQOGZwWTczbzNyX1Q93Jznn/6x7n3Wee+659/b5/2CYu/fZ+3nWfllnn7X2ejF3hxDi/z9zOy2AEGI2SNmF6AQpuxCdIGUXohOk7EJ0gpRdiE7YkrKb2U1m9hMze8LMjk9LKCHE9LHW9+xmNg/gpwDeAeAsgO8CeJ+7/3ijfRbnlvyyuUNN843MnNgmOqbMfhXW1ppmHpNiLvh+rYm2FozK182mcHytbJMs2xEDYnORbFOQtxii9T6dfJNRXlx7FpfWLoST75tsqDGuB/CEuz8JAGZ2N4BbAGyo7JfNHcIfHrplC1MCCC8WESlHZj/GxhXTX3yxuouvbv6FMHfwsnLl/PzmY164WK5cXR1fXlgot+EvJ1aeaN7aPpEisyz8hZZR/kCxnceNrmsNuu5z+/dXt+HrDg+uKW/DY8xtfk03HJfPw4TH/J1n/mXDz7byM/4aAD8fWT47rBNC7EK28mSPvq6LryEzux3A7QCwZAe3MJ0QYitsRdnPArh2ZPkIgKd4I3c/AeAEAPzGvqt87Kdj9DOGf1ryz5j5xh8j/PNoeWV8ueFnfotN6TwvAFR++ofwT+XAn2B0LlPyRj6FyjzFNeN5onkTP+2Nt5kfX7bFxWKfbcn1qJhZIdG9zbApAAC+Wq6rMXrvbnJat/Iz/rsAXm9m15nZIoD3Arh3C+MJIbaR5ie7u6+Y2V8B+HcA8wA+7+6PT00yIcRU2crPeLj7NwF8c0qyCCG2EUXQCdEJW3qyT4yBHGENjo8MkZOm8r4yfD++Nu5MY+cPO8AAwDHuYCmcTNOCxw0ca1NxVrFDLpindoyhHPwOPaLhfX31fEfOtuI9+zZds8Jx2eCMixh7p7+x7HqyC9EJUnYhOkHKLkQnzNZmhwH7Xpkysq8K+45tuzAQIRHAQOMUcdcJMvY3b1PY+aE/YfKgmsJfkAkK4nMXnIPi/NeCbCJINovmSfgcivOdOUY+Lw3XObzHavA9GI5B20R+otoxRp8nA8L0ZBeiE6TsQnSClF2ITpixzd5AYeNOKQe4peBCg/2aes/O4xZ+isQY+4JLeWmZ5qn7BqryRvYhn++EnZzyOdRs5ygpiq5rET+xSucEgC1SLQDKRbdgnlrdgpQfKaLmo4oSqcau2cb3vp7sQnSClF2ITpCyC9EJUnYhOmG2Djr3MaeRR06lzBg1ImdPzWfUUk0lCoZpKbbI02bOCx/jtAI0Woo6FpVqxmUJq/PQebJMNdaWoqFMFMhF8hmdAw/mZQdd4cRrqW4D5IpdMqOybHL59GQXohOk7EJ0gpRdiE6YcVCNjwcNtNiMmX2iYAveL1GUoSkJpDJGGLRSa1LQaquuVeSP7MHaXFGgS+E/2HyIUJTA52D7uIIu2dIrZYBM9bomrmmZjFVPXnJaDs9iS4JNsU+QVLQy4nPYxKelJ7sQnSBlF6ITpOxCdIKUXYhOmHFQDTk/OCsLQXACOyjC9sXkLHmh7Hpa7QgaODaqGWBNLYwSwRatnWoZdsAVrbUSzjaWP+Nk4qyxhbJNE5YvkSz1CrRF1Z8oQIZXsOMscf6nUZU37LzLBLLYwuYqGcqWrHSkJ7sQnSBlF6ITpOxCdMLsK9WM2mGRzVusI1u7pcIogqqv7C+IbMaW9sVFddZEhVSeu8Uez5AZl7cpTkFDm+3UPFNKykm0smZqNnroG6hVF4ooqvPUK+BEVXLKnXLnSU92ITpByi5EJ0jZheiEXVddNvMutWAt8f616MA6ebGBwu6PtuHCE4U9HhVyaEh0ydizLcUrCPZtFJVYM7IkEm4yBS5SVGz0sBNQJlmJ4PuHE2EyiVWpe5viGiy4TYu5N0BPdiE6QcouRCdI2YXohKqym9nnzey8mf1oZN2VZna/mZ0Z/n/19oophNgqGQfdFwD8I4B/Hll3HMAD7n6nmR0flj+amnHUKTGtNr1EGCTBDriW6p8ki2HyMSLZrCE+pupwBBKVSgP5a87OlMO0knQUkXDGNbXSStDU8qoFrm4TyVqreBO2bM4dc3Urd/8WgF/S6lsAnBz+Pgng1tRsQogdo9Vmv9rdzwHA8P9rN9rQzG43s1NmduqSX2icTgixVbbdQefuJ9z9qLsfXbSl7Z5OCLEBrUE1vzCzw+5+zswOAzg/LYEKW5TNzOZKq1MoYtCSmDENgmNusfOb4PMUtipuaF+cCISaGUXl4cR1p+CiootMIngnloWCdVgfMp1zNhq6cb97ARwb/j4G4J5mCYQQMyHz6u1LAL4D4HfN7KyZ3QbgTgDvMLMzAN4xLAshdjHVn/Hu/r4NPnrblGURQmwjO5sIkwngbymMEFGxl8L33/zuehrJJ0F3kWryT1hYo6EQJG+TKLhQFFNosK3DczuNQpxB8kyRoJKRv+jAQ+/Do241m4i5EZkYkmIu1pFMMYsNULisEJ0gZReiE6TsQnSClF2ITthZB10iEYDJOHtSFTl5jIzLhZ1IUbBIxQmWcXA1VeuJqAYgBU6xmrxRpZqVoMrM2JiZxKTENWPnVaYiLXeVCc5lKtmE8Bcp9Luli0zgeK06LrcQ2KUnuxCdIGUXohOk7EJ0wmxtdhtPQAmLTFQSAaYnS62wQ8MYEUUXkHoXHAu62xa7LFEGYWTXU6VbP0D7RP4RlmWFAm+iKrBclOHieAdTW5jSM2UfBzkFQUFG9wv7E4L7iZNWiqIkmW4vU+g8A9T9M1vRBz3ZhegEKbsQnSBlF6ITZmuz++Q2Ryphgu1t7soCBMkONC6/N42gBInQZi+6to7L+/eP/Guxyzz1lrl6ftxmnw/suFU6j5FVuVw51ctBbMFC2OfmFTKeDT4r84EZeiARO3DAxt/pv+Dj5+XJlcVin7/7gz8eX8HzRPcfrSvu0UhWjhHZyeIbSfRkF6ITpOxCdIKUXYhOkLIL0Qm7vlINd80I3SBFgkrkhCEX1gptk0mq4OXAV2iL5DQi+ReCsrDsoGOH3P7AEbgGdtCVx7zY5DMa3+kCncuoGPgFmvrgXF1+Zi24Zhd9PCBmnmT75eqh6riFIy2652pO46hzEVfE4Q0SgTizrqirJ7sQnSBlF6ITpOxCdMLO2uwJCrsmtIUSRSWIoupoVDyhVsE13GfzYIuDViaSXEn77LegQAQxT7Its08iQWTnL9O5WyL5D1jpqLi8YR62vyNnDNvxy4WfItjJN7+uoXXOdnym6nFRAZgDrqID4kqxgdOHrmu12uwE6MkuRCdI2YXoBCm7EJ2w+2x2Ll6RKRyQSa5hW6eYJ+j6UXTwIHtwX5mIUbPLDs+X+/B79bmG7+DUu2xKY4k8A7VEmIgoUWeMBn9CxCW6zpfPBclLtaIkmffsLe+/OR4kGKNIsEkUtjQuQEKFQQDkuulAT3YhukHKLkQnSNmF6AQpuxCdMHsH3RaCAjaEnXgtDpZILnZ8FEE1UVBEpSJOA+xYi9gXZeUQCxSsc9HLKrbsbOOKOBG1bTIOx2iLixgPQFok2ZYsqMJbtHWm7jrBNat14JlGVdjsNoxfusSDlBsldUpPdiE6QcouRCdUld3MrjWzB83stJk9bmZ3DOuvNLP7zezM8P+rt19cIUQrGZt9BcCH3P37ZnY5gO+Z2f0A/gLAA+5+p5kdB3AcwEcnmj0TDMP2VFA5trCpokAc7nTCdlnUUbaYKBGgsW/yrpts007LRp+vBNpEY6xQndqMLDUyNns0zxz7D+h8/9Z8EGDS0OV01kUkXqalI2ukM0n5q1fB3c+5+/eHv58FcBrANQBuAXBy2OwkgFtTMwohdoSJbHYzex2ANwN4GMDV7n4OWP9CAPDaaQsnhJgeaWU3s0MAvg7gA+7+6wn2u93MTpnZqUueaMQghNgWUspuZgtYV/Qvuvs3htW/MLPDw+eHAZyP9nX3E+5+1N2PLlpUrlAIMQuqDjpb917cBeC0u39y5KN7ARwDcOfw/z21sdwdPtJCNwpwKKBtogCHwtkWjNvkhKkFxETyB9lzo0SOKG7TVM0iQ1lddiGoILNaqdjDY6zLskrL42Ow0wwIqs4U8wTOt4a3vux2je4e2z+eVbj27HMTz5OhqHqcub8yAVbstEtUxx3Vqc0CgDLe+LcA+HMA/2lmjw3r/hrrSv4VM7sNwH8D+LPEWEKIHaKq7O7+bWxQrh3A26YrjhBiu1AEnRCdMNNEGDMbt9MDO6clCSFDkexAwTlhdVmem+wpv1C+XbAFqv9SdIQpLc2WoJq5ip28Po7T8uQBMpGNXoOP8YKXFXVXOXgncZ35zAVpMMCBy8aXp2Czp3xLGXucxwmDasgzscbVcjeXxS5tLIee7EJ0gpRdiE6QsgvRCbMtXmFkK2fsnNW6LZeqQFvsRParN7wnDexMXyFZ5up2MtvSXAwieu9eJKwEsrR0iWH/AdvWEWyj8xgHrKyoy7KtBd1tuZMM3wrRk8p/9QxtVH+elVWEK52AIookqeC+XSPfRaa67EJdRW3Um7HJkHqyC9EJUnYhOkHKLkQnSNmF6ITZOuh8PLjFEoEhqSCbFidMxulSJCVMHmDC80YJIOyg22/76PO6k7IlYCaVcEOOyzjhZly+VHwJz+3leWHnIFeXDavaVpy1YfDUFChagEfntuH+SY2bRE92ITpByi5EJ0jZheiE3deymW1pbpMcmaaZIIiWziy1fRJ2GdtcUReWcohx+aOkl4wdX+vuEtm8vE8tYAZoPLWJ4B0uisEztzSCDltzT8GOz/iWins3E1Qzxcq3erIL0QlSdiE6QcouRCfM1mafM9jiwubbsI2SqTJYKawYjpt4z27c3SVR/K+AdnkhSE5ZINt0FeOdOw/Z/uo0mW6rbH9H7+ZrnWY4ASc7d02W8J0/vXufJ3mfnVJhkyKeoqVYBSdSBQlcZeeihu6wkZ0/di+reIUQ3SNlF6ITpOxCdIKUXYhOmHkizJgzJHKsZZwjTBGIM6UAmqK67NbbFz8bVDC5nMRfoE2eQ9maeAF1ZxsHrnBwzlyi7XNtTABFtEiLwy4DJwhdHYlP909xVaNEmZZ7rtJuOQqGKRx00f3ETmC+LyNZk62f9WQXohOk7EJ0gpRdiE7YfYkwNSI7n22WqFIs2zoJs77oIlMEOATflWyH0Tb7E/MuU5LL/kDY5VQayBot1W27TKeZcp/JK9LWKupG26Q6v3JSSxHIEuzDdjzfK5E/p2YnB/sUiTAtvoLaXKouK4SQsgvRCVJ2ITph5h1hxm2ZwGYp7GL6fC36fiKbKyxeQe9fWwoWkP3nK2V3UmM7nhJ/Wt7UR91eFozt5HIbXpcpnLHfxuVtseE5qSXqTLNMx8RxAxFswy9nEqASFL4ZHjey89nebqmkMWP0ZBeiE6TsQnSClF2ITqgqu5ktmdkjZvYDM3vczD4+rL/OzB42szNm9mWzoFWnEGLXkHHQXQRwo7s/Z2YLAL5tZv8G4IMAPuXud5vZZwHcBuAzE80+lwkqYCdM1BEm4URaG/egFKNEY1wih1bCqcfjGrVwDt2LlcSRQ3NldR92VkUOLg68YYcWO/kiWrrRcPDLBR+vvNM6TwquHpu6N8b3yUhitXOXTE6pT8TVkYJxk87m6la+znPD4sLwzwHcCOBrw/qTAG5NzSiE2BFSXwlmNm9mjwE4D+B+AD8D8LS7v/Tu6SyAazbY93YzO2Vmpy6tXZiGzEKIBlLK7u6r7v4mAEcAXA/gDdFmG+x7wt2PuvvRxbmldkmFEFtioqAad3/azB4CcAOAK8xs3/B0PwLgqdQgowEXLUERYXVNOoyg60eR7FBUFE18702hKMPFyOVAy0sJM7MMoqmfS7bRMx1lM8knvA0H0XBnF6CUPwoKivYbmyd6vvB1jO4FptYFOCoyUekW69HnLd2GmS10r8l4468ysyuGvy8D8HYApwE8COA9w2bHANzTLIUQYtvJPNkPAzhpZvNY/3L4irvfZ2Y/BnC3mf0DgEcB3LWNcgohtkhV2d39hwDeHKx/Euv2uxBiD6AIOiE6YcZZbwbb98qUUVvbmuMjbJPMRE4Mrl5TVLdJON94n0wFk6k49crsurnEeWAHV8qpRy2XuLxK2LKZ5mF5I+fbkk1epYWDgqJKvUU2GmUmhvcck6giHDrgRml00BUZeNU9MKZTm+mHnuxCdIKUXYhOkLIL0Qk7W102srlq3V0iW4j3yQTr8DjBsGyXGduD80H7aZ7bOXgkCjAZZ4ls3Ch4JKo4OylxdZtxaRYwnswYVZ15Zm3zRJfIul2jtQcSNnwkL8NBKV4kkpT+j5pfZcwmfnkXrlqbKFXD22T8T4lxfdSvssmh6MkuRCdI2YXoBCm7EJ0w4y6uDl8eKQiRKV7Bdn2i4AK/Ww3HSVDY6Ow/CKvYbm6HrQbdahbpXfYFEvVgo3les3GjqrW8zzM+npYceUO4Zu1q5lU2zXMhiCWo8byXt+/a8y+MLVcTS4DyvfoC+WKia7rMiVUNSV2Jd/HF+/zgmlnl85eHnkA0IcQeRsouRCdI2YXoBCm7EJ0wcwcdlkcdMYFTZoFEYgdL5IDIOGHIyeJcwSTjYKEkkbml4LuSZaFxP3zDn9anIdlsf6JKd+t5qY3DlVdXAqcSO0Sp5VXoHOVqvxcuFpvYIh03B7csR8E8JC8fT5SMUkuOyTh3E0kutXbeKQKnno+OKwedEELKLkQnSNmF6IQZ2+zU5jiwkwsrM9PRg8nYWMsUCtJgc/nFwGZkednOjGRrOcZM+2i22TlIKLLp+RjZlo5k5eQf7qQT7VNLMorGZRs9OpctySYVP4tH9wav43s5USQjKoCRCgJionMXoCe7EJ0gZReiE6TsQnTC7N+zj9g2YfE/tj3JHolsmmKchqSElP3EdlpUJGONO5IkihqsVOy76DxliiLyfnxuE++/i2Pk44vGmUvsk6Hie0kVj8zA4/B1z1zDFls7A9vjwTGP+ms2Oyd6sgvRCVJ2ITpByi5EJ0jZheiEmTroHJTkETi42L3Abg9vdYRUnCypYAZ2loQOJJonEzBTOIQ4+SQRMBNRCy7KVOFtqdzL0zTKX6vgWu3KEsyTceqlzm1DEkuTQzHjHEyiJ7sQnSBlF6ITpOxCdMLsO8KM2LBmiU4bLQn/YVLFFAIwEgkSxh1kuYJuYMOXXUzomKOCBbXOOUCTfV12oW0I3qHyssU5CeZJ2bN0nsLkmcr9YlEwUu2eShU2SfgCMgkrlbm2EkikJ7sQnSBlF6IT0spuZvNm9qiZ3TcsX2dmD5vZGTP7spklCqUJIXaKSWz2OwCcBvCqYfkTAD7l7neb2WcB3AbgMxPNHtm821AkMaR4Zx7s02IfFUkhZGdO4/gAYJXeXQedRquEHW34XDYk9jCJgh0WdaVlu57vl0RHlYk/nxbRda4l3ABlN5qi23CQsDVyL292f6WO3MyOAHgngM8NywbgRgBfGzY5CeDWzFhCiJ0h+zX3aQAfwSt1el8D4Gn3lxt0nQVwTbSjmd1uZqfM7NQy9Q0TQsyOqrKb2bsAnHf3742uDjYNf/O6+wl3P+ruRxdsqVFMIcRWyRh6bwHwbjO7GcAS1m32TwO4wsz2DU/3IwCe2j4xhRBbpars7v4xAB8DADN7K4APu/v7zeyrAN4D4G4AxwDcUxvLQA6EyJlQcaCEDgjuIhPBjg12/mScPQlHIDuRjJM3gsCKqtOuod00ABg57YouOJm5KKhmatVhGo9pjExQDdMQvBPC8xTtvRs6xETrJg0k2qZKNR8F8EEzewLrNvxdWxhLCLHNTPS+xt0fAvDQ8PeTAK6fvkhCiO1AEXRCdMJsE2HMxm2bhF1T2LMt9nkA256hLcpdS2odVoDS5mI7OQhkcX6RwcknLR1jEHUwbaiaSjQFBc2SWqeWTCJM5rzwPrV5M2NE+2XGHd1mq0E1Qoi9j5RdiE6QsgvRCbMvXjHC1Io88rvgzPvkBrhwYvjOs9Y1hjuRAvDlxLjFPNytJrEP26JBIciCokhGw7v6APYnpDr9FPM0FKKIZKn5ZlrI2OMNnYtC1MVVCDGKlF2ITpCyC9EJUnYhOmEHqsvWKomw4ykRCEJOJIv8FZnEFx6WHR9JR8hm84bUOp1EFX24Ms3Fi+W4RbtfcmQuBw66IsBkudxmUlrOW0QRsJQIfplCB5jIUcj7lO29gzET3Wmqsmz66eboyS5EJ0jZhegEKbsQnbCjQTWhLce2KAelhMknbI/XiwKsPf/82HJhAwNlUkFDR4/CxopsyCkEgsQbNVTdrYwbzlvxS4R2KB1zOG7hy6gnCHGxkNq8g4CbbhOOWUmkCqvlckXdlpiaTLDORrs2TCeE2INI2YXoBCm7EJ0gZReiE2ZeqWasfW+mAktR2SVwRlArpDCDihwbcwcP1ueuBWRknCWZirSVeTIOrrBCCTuWKo6ocJ9EIAtnAxbzBG2fbXHy1oDc+jlyxBZVf4o2UwG1FkstrZajNtUtFO3DAllG5d1EpfRkF6ITpOxCdIKUXYhOmH1QTdQmeBS2/9jea+2OUgmQyVSXTQWp1LaJkloS2xQk9ikqsFDgUOjbqNn1kf1aqXjjwTUrjrmFlqq706oOw4E3fN7Cdti0zVpDpaAtoCe7EJ0gZReiE6TsQnTC7G32URt8LfiuCd7JVoecRmfRyJZj+7T2PjZDYCfX3jkX77GBVKeQqVRJzYxR6zSasTsTfopU8g/5ByzTQYhp6eLaMk9jp5+CsfOtjjBCdI+UXYhOkLIL0QlSdiE6YWeDahqccc3T1pwuicQMXx6vtJpxDKacZEWwBY0bjZFxSi4s1LepwRVYoqAaPneNgU8FRbukSgvqaF2RoJJo2dUCy7aYGHNaATOjx6xEGCGElF2ITpCyC9EJNpWAlOxkZv8D4L8A/CaA/53ZxFtjL8kK7C1595KswN6Q97fd/arog5kq+8uTmp1y96Mzn7iBvSQrsLfk3UuyAntPXkY/44XoBCm7EJ2wU8p+YofmbWEvyQrsLXn3kqzA3pN3jB2x2YUQs0c/44XohJkqu5ndZGY/MbMnzOz4LOfOYGafN7PzZvajkXVXmtn9ZnZm+P/VOynjS5jZtWb2oJmdNrPHzeyOYf1ulXfJzB4xsx8M8n58WH+dmT08yPtlM5u8oPw2YWbzZvaomd03LO9aWTPMTNltvbr9PwH4EwBvBPA+M3vjrOZP8gUAN9G64wAecPfXA3hgWN4NrAD4kLu/AcANAP5yOJ+7Vd6LAG50998H8CYAN5nZDQA+AeBTg7y/AnDbDsrI3AHg9Mjybpa1yiyf7NcDeMLdn3T3SwDuBnDLDOev4u7fAvBLWn0LgJPD3ycB3DpToTbA3c+5+/eHv5/F+k15DXavvO7uzw2LC8M/B3AjgK8N63eNvGZ2BMA7AXxuWDbsUlmzzFLZrwHw85Hls8O63c7V7n4OWFcwAK/dYXkKzOx1AN4M4GHsYnmHn8WPATgP4H4APwPwtLu/VHdrN90TnwbwEbzSb+w12L2yppilskfJd3oVsEXM7BCArwP4gLv/eqfl2Qx3X3X3NwE4gvVfem+INputVCVm9i4A5939e6Org013XNZJmGU++1kA144sHwHw1Aznb+UXZnbY3c+Z2WGsP5V2BWa2gHVF/6K7f2NYvWvlfQl3f9rMHsK6r+EKM9s3PDF3yz3xFgDvNrObASwBeBXWn/S7UdY0s3yyfxfA6weP5iKA9wK4d4bzt3IvgGPD38cA3LODsrzMYEPeBeC0u39y5KPdKu9VZnbF8PdlAN6OdT/DgwDeM2y2K+R194+5+xF3fx3W79P/cPf3YxfKOhHuPrN/AG4G8FOs22p/M8u5k/J9CcA5AMtY/yVyG9ZttQcAnBn+v3Kn5Rxk/SOs/4z8IYDHhn8372J5fw/Ao4O8PwLwt8P63wHwCIAnAHwVwP6dlpXkfiuA+/aCrLV/iqATohMUQSdEJ0jZhegEKbsQnSBlF6ITpOxCdIKUXYhOkLIL0QlSdiE64f8AF9CkTcUempMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Your work here!\n",
    "x0 = zeros((50,50))\n",
    "f = lambda x: tv_denoise_objective(x,mu,image)\n",
    "grad = lambda x: tv_denoise_grad(x,mu,image)\n",
    "x, res = grad_descent_bb(f,grad,x0)\n",
    "print(\"Solver terminated in %d steps\"%len(res))\n",
    "\n",
    "print(\"Here's the denoised image...\")\n",
    "plt.imshow(x)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}