{ "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": [ "