\n",
"Theorem: If $f : \\R^d \\rightarrow \\R$ is convex, differentiable and L-Lipschitz on all of $\\R^d$ and $||x_0 - x^*|| \\leq R$ where $x^*$ is the global minimum, then there is a step size $\\eta > 0$ such that the iterates of gradient descent satisfy\n",
"$$\n",
"f\\left(\\frac{1}{t}\\sum_{i=0}^{t}x_i\\right) - f(x^*) \\leq \\frac{R L}{\\sqrt{t}}\n",
"$$\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Problem: Tuning the Step Size\n",
"\n",
"Last time, we were able to show that for any step size $\\eta > 0$,\n",
"\n",
"$$f\\left(\\frac{1}{t}\\sum_{i=0}^{t}x_i\\right) - f(x^*) \\leq \\frac{R^2}{2\\eta t} + \\frac{\\eta L^2}{2}$$\n",
"\n",
"
\n",
"Problem: Find the choice of step size $\\eta > 0$ that gives the best possible bound on the convergence rate. What is the corresponding bound?\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Solution: Tuning the Step Size\n",
"\n",
"Define $g(\\eta) = \\frac{R^2}{2\\eta t} + \\frac{\\eta L^2}{2}$ and differentiate:\n",
"\n",
"$$\n",
"\\frac{\\partial g}{\\partial \\eta}\n",
"= \\frac{-R^2}{2\\eta^2 t} + \\frac{L^2}{2}\n",
"$$\n",
"\n",
"Setting the derivative to zero and solving for $\\eta$ gives $\\eta = \\frac{R}{L\\sqrt{t}}$. Plugging this into $g$ proves the theorem!"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Choosing the Step Size\n",
"\n",
"In general, there are several strategies for choosing the step size:\n",
"\n",
"* **Constant:** Choose the same step size $\\eta_t = \\eta$ for all times.\n",
"* **Theoretical:** Choose a step size that guarantees a nice convergence rate.\n",
"* **Shrinking:** Choose some sequence $\\eta_t \\rightarrow 0$ to guarantee convergence.\n",
"* **Line Search:** Solve a nested minimization to pick the best step size in the current direction: $$\\eta_t = \\arg\\min_{\\eta \\in (0,1]} f(x_{t-1} - \\eta \\nabla f(x_{t-1}))$$\n",
"* **Adaptive Step Size:** Use more information about the function/problem to adjust the step size as needed. For example, *Newton's method* adjusts the step size based on the curvature of $f$ at the current iterate."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Remark: Projected Gradient Descent\n",
"\n",
"Suppose we want to minimize $f : \\Omega \\rightarrow \\R$ within some convex set $\\Omega \\subset \\R^d$ instead of the entire space $\\R^d$. We can still apply gradient descent as long as we **project** back onto the convex domain $\\Omega$ after each iteration:\n",
"\n",
"
\n",
"Here, $\\Pi_\\Omega(x) = \\min_{z \\in \\Omega} \\norm{z-x}_2^2$ projects onto the convex set $\\Omega$.\n",
"\n",
"
\n",
"Problem: Argue that Projected Gradient Descent achieves the same convergence rate by making a minor adjusment to the proof.\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Assumptions about the Objective\n",
"\n",
"When minimizing an objective $f : \\Omega \\rightarrow \\R$ using gradient descent, we can obtain stronger convergence guarantees by making assumptions about how \"nice\" of a function $f$ is:\n",
"\n",
"* **Convex:** $f(\\theta x + (1-\\theta) y) \\leq \\theta f(x) + (1-\\theta) f(y)$ for all $x,y \\in \\Omega$ and $\\theta \\in [0,1]$\n",
"* **$L$-Lipschitz:** $\\norm{\\nabla f(x)} \\leq L$\n",
"* **$\\alpha$-Strongly Convex:** $\\nabla^2 f(x) \\succcurlyeq \\alpha I$\n",
"* **$\\beta$-Smooth:** $\\nabla^2 f(x) \\preccurlyeq \\beta I$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Strong Convexity & Smoothness"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Problem: Matrix Inequalities\n",
"\n",
"Remember the notation $C \\succcurlyeq D$ means $C - D$ is positive-semidefinite. Below, let $A \\in \\R^{n \\times n}$ be symmetric.\n",
"\n",
"
\n",
"Part A: For some $\\alpha \\in \\R$, when is $\\alpha I \\succcurlyeq 0$? When is $\\alpha I \\preccurlyeq 0$? \n",
"\n",
"Part B: If $A$ has eigenpairs $A v_k = \\lambda_k v_k$ for $1 \\leq k \\leq n$, what are the eigenvalues and eigenvectors of $A - \\alpha I$ when $\\alpha \\in \\R$? \n",
"\n",
"Part B: Let $A \\in \\R^{n \\times n}$ be symmetric. For fixed $\\alpha > 0$, which conditions guarantee that $A \\succcurlyeq \\alpha I$? What about $A \\preccurlyeq \\alpha I$?\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Solution: Matrix Inequalities\n",
"\n",
"**Part A:** The eigenvalues of $\\alpha I$ are all equal to $\\alpha$, so this matrix is psd when $\\alpha \\geq 0$ and nsd when $\\alpha \\leq 0$.\n",
"\n",
"**Part B:** The matrix $A - \\alpha I$ has eigenvectors $v_k$ and eigenvalues $\\mu_k = \\lambda_k - \\alpha$, since\n",
"$$\n",
"(A - \\alpha I) v_k = A v_k - \\alpha v_k = (\\lambda_k - \\alpha) v_k\n",
"$$\n",
"\n",
"**Part C:** By definition, $A \\succcurlyeq \\alpha I$ iff $A - \\alpha I$ is positive semidefinite. By the previous part, this requires $\\lambda_k \\geq \\alpha$ for all $k$. Similarly, $A \\preccurlyeq \\alpha I$ whenever $\\lambda_k \\leq \\alpha$ for all $k$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Def: Strongly Convex\n",
"\n",
"Remember that a differentiable function $f : \\Omega \\rightarrow \\R$ is convex iff it lies above its linear approximation at every point, that is,\n",
" $$\n",
" f(y) \\geq f(x) + \\inner{ \\nabla f(x), y-x } \\quad \\forall\\, x,y \\in \\Omega\n",
" $$\n",
" \n",
"A function is **$\\alpha$-strongly convex** if it also lies above a *quadratic* approximation at every point, where $\\alpha$ controls the steepness:\n",
" $$\n",
" f(y) \\geq f(x) + \\inner{ \\nabla f(x), y-x } + \\tfrac{\\alpha}{2} \\norm{y - x}^2\n",
" $$\n",
" \n",
"
\n",
"Problem: Prove that every strongly convex function is convex.\n",
"
\n",
"Problem: Show that if $f$ is twice-differentiable and $\\alpha$-strongly convex, then $\\nabla^2 f \\succcurlyeq \\alpha I$.\n",
"
\n",
"\n",
"> *Hint:* First, show that $f$ is $\\alpha$-strongly convex if and only if the function $g(x) = f(x) - \\frac{\\alpha}{2}\\norm{x}_2^2$ is convex. Then, use the second-order condition for convexity on $g$.\n",
"\n",
" \n",
"**Reminder**: A function is **$\\alpha$-strongly convex** if it also lies above a *quadratic* approximation at every point, where $\\alpha$ controls the steepness:\n",
" $$\n",
" f(y) \\geq f(x) + \\inner{ \\nabla f(x), y-x } + \\tfrac{\\alpha}{2} \\norm{y - x}^2\n",
" $$\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Def: $\\beta$-Smoothness\n",
"\n",
"A function $f : \\Omega \\rightarrow \\R$ is $\\beta$-smooth if it lies *below* a quadratic at each point:\n",
"\n",
"$$\n",
"f(y) \\leq f(x) + \\inner{ \\nabla f(x), y-x } + \\tfrac{\\beta}{2} \\norm{y-x}_2^2\n",
"$$\n",
"\n",
"Similarly to strong-convexity, a twice-differentiable function is $\\beta$-smooth if and only if $\\nabla^2 f(x) \\preccurlyeq \\beta I$ for all $x \\in \\Omega$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Problem: Quadratic Forms\n",
"\n",
"Let $A \\in \\R^{d \\times d}$ be symmetric and positive-definite and define $f(x) = \\frac{1}{2} x^T A x$. \n",
"\n",
"* Is $f$ smooth? If so, with what smoothness constant?\n",
"* Is $f$ strongly convex? If so, with what constant?\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Gradient Descent for Smooth and Strongly-Convex Functions\n",
"\n",
"Let $f(x)$ be a $\\beta$-smooth and $\\lambda$-strongly convex function with minimizer $x^*$. Assume our initial point is $x_1$, and repeatedly perform gradient descent: $x_{t+1} = x_t - \\eta \\nabla f(x_t)$. Then it can be shown that\n",
"$$\n",
" f(x_T) - f(x^*) \\leq \\frac{\\beta}{2} \\exp(-4T/(\\kappa + 1)) \\|x_1 - x^*\\|^2.\n",
"$$\n",
"where $\\kappa = \\frac{\\beta}{\\lambda}$ is the *condition number*, and the learning rate needs to be set as $\\eta = \\frac{2}{\\lambda + \\beta}$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Frank-Wolfe / Conditional Gradient\n",
"\n",
"Let $f(x)$ be a $\\beta$-smooth function, and assume we want to solve $\\min_{x \\in \\Omega} f(x)$ for some constrained set $\\Omega$ with diameter $D$. The *Frank-Wolfe* algorithm does the following:\n",
"\n",
"+ Initialize $x_1$ to some arbitrary point in $\\Omega$\n",
"+ For $t=1,2,\\ldots, T$:\n",
" + compute $v_t + \\arg\\min_{v \\in \\Omega} v^\\top\\nabla f(x_t)$\n",
" + update $x_{t+1} = x_t + \\eta_t (v_t - x_t)$\n",
" \n",
"**Theorem**: after $T$ rounds of Frank-Wolfe, we have that\n",
"$$ f(x_t) - f(x^*) \\leq \\frac{2 \\beta D^2}{t + 2}$$\n",
"for step size $\\eta_t = 2/(t+2)$.\n",
"\n",
"[more](http://fa.bianp.net/blog/2018/notes-on-the-frank-wolfe-algorithm-part-i/) [information](https://ee227c.github.io/notes/ee227c-lecture05.pdf)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Implementation\n",
"\n",
"> Adapted from [Moritz Hardt's lecture notebook](https://ee227c.github.io/code/lecture4.html)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
"## Projected Gradient Descent\n",
"\n",
"We start with a basic implementation of projected gradient descent. Note that this implementation keeps around all points computed along the way. This is clearly not what you would do on large instances. We do this for illustrative purposes to be able to easily inspect the computed sequence of points."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def gradient_descent(init, steps, grad, proj=lambda x: x):\n",
" \"\"\"Projected gradient descent.\n",
" \n",
" Inputs:\n",
" initial: starting point\n",
" steps: list of scalar step sizes\n",
" grad: function mapping points to gradients\n",
" proj (optional): function mapping points to points\n",
" \n",
" Returns:\n",
" List of all points computed by projected gradient descent.\n",
" \"\"\"\n",
" xs = [init]\n",
" for step in steps:\n",
" x_step = xs[-1]\n",
" x_update = None\n",
" xs.append(x_update)\n",
" return xs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Warm-Up: Optimizing a Quadratic\n",
"\n",
"As a toy example, let's optimize $f(x) = \\frac{1}{2} \\norm{x}^2$, which has gradient $\\nabla f(x) = x$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def quadratic(x):\n",
" return 0.5*x.dot(x)\n",
"\n",
"\n",
"# What is the gradient of this function?\n",
"def quadratic_gradient(x):\n",
" return None\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note the function is 1-smooth and 1-strongly convex. Our theorems would then suggest that we use a constant step size of 1. If you think about it, for this step size the algorithm will actually find the optimal solution in just one step."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"x0 = np.random.normal(0, 1, (1000))\n",
"_, x1 = gradient_descent(x0, [1.0], quadratic_gradient)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Indeed, it does:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1.all() == 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's see what happens if we don't have the right learning rate."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAH/CAYAAACW6Z2MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4m+XZBfBzS94zw87eey87HowyQykQNnH2tg1llRYopUBL6QcUChRoKbazE5I4ZTZQ9izgETs7IXtvO068p/R8f1h+Jad2EtuyH43zuy5dlm69kk5acA7veCRKKRARERG5KpPuAERERETnw7JCRERELo1lhYiIiFwaywoRERG5NJYVIiIicmksK0REROTSWFaIiIjIpbGsEBERkUtjWSEiIiKX5qM7gLcTkYkAJoaGhiYOGjRIdxwiIqI2k5ubm6+UirzQdsLl9l1DdHS0ysnJ0R2DiIiozYhIrlIq+kLb8TAQERERuTSWFSIiInJpLCtERETk0lhWiIiIyKWxrGgmIhNFJLWwsFB3FCIiIpfEsqKZUmqtUiopPDxcdxQiIiKXxLJCRERELo1lhYiIiFwaywoRERG5NJYVzXiCLRER0fmxrGjGE2yJiIjOj2WFiIiIXBrLChEREbk0lhUiIiJyaSwrRERE5NJYVjxQcUU1nvvPT6iotuiOQkRE1GIsK5o5+9Ll0soazF68Dinf7cP8pTkor2JhISIi98ayopmzL12+d+V65B48AwD4fk8+5izJRmlljVPem4iISAeWFQ+T/LP+CPIzG48z9xVg1qJsFFdUa0xFRETUfCwrHia+f0csmxuDEH8fY5Zz8AymL8xGYRkLCxERuR+WFQ8U3acDls+LQWiAvbBsOnwW0xZm4kxplcZkRERETcey4qHG9mqPVYlxaBfka8y2Hi3ClLRM5JdUakxGRETUNCwrHmxE93CsSoxDh2A/Y7bjRDGmpGbiVFGFxmREREQXj2XFww3tGobVSXGICPE3ZrtPlWByaiZOFLKwEBGR62NZ0czZ66w0ZFDnUKQnx6FzmL2w7MsvxaSUDBw5U9Zqn0tEROQMLCuaOXudlcb0jwzBmuR4dG8XaMwOFZQhISUTh06zsBARketiWfEivTsGIz05Dj072AvL0bPlSEjNwP78Uo3JiIiIGsey4mV6tA/CmuR49I0INmbHCyuQkJKBPaeKNSYjIiJqGMuKF+oaHoj0pDj0j7QXllPFlZicmomdJ1hYiIjItbCseKlOYQFYnRSPwZ1DjVl+SRUmp2Zg69HWO9mXiIioqVhWvFhkqD9WJcVhWNcwY3amrBpT0zKx6fBZjcmIiIjsWFa8XIdgP6xKjMPoHvarkYoqajB9QRZyDxZoTEZERFSLZYUQHuSL5fNjEdW7vTErrqzBzIXZyNp3WmMyIiIilhWyCQvwxdK5MYjp28GYlVZZMHvxOvywJ19jMiIi8nYsK2QI8ffBkjnjcemAjsasvNqCuUvW4dtdeRqTERGRN2NZ0awtlttviiA/HyycNR5XDIo0ZpU1ViQuzcEX209qTEZERN6KZUWztlpuvykCfM1InRmFa4d2MmZVFivuXpGLT7Ye15iMiIi8EcsKNcjfx4w3pkXhFyO6GLMaq8K9Kzdg7aZjGpMREZG3YVmhRvn5mPD6lLGYOLqbMbNYFR5cvQHvrj+iMRkREXkTlhU6Lx+zCX9LGIPbx3Y3ZlYF/OZfm7Bm3WGNyYiIyFuwrNAFmU2CF+8ajYTonsZMKeDRdzZjReZBjcmIiMgbsKzQRTGbBM/dPhLT43rVmz/x/lYs/mG/plREROQNWFbooplMgmduGYG5l/atN3967XakfLtXUyoiIvJ0LCvUJCKCJ28airuv6F9v/tzHO/D6l7s1pSIiIk/GskJNJiL47fWD8cDVA+rNX/p8F17+bCeUUpqSERGRJ2JZoWYREfz6usH4zYRB9eavfbUHf/mEhYWIiJyHZYVa5P5rBuJ3vxhSb/bmt3vxzIc/sbAQEZFTsKxQiyVf0R9P3TSs3mzRD/vx1AfbYLWysBARUcuwrJBTzL2sL565dUS92fLMg3j8vS0sLERE1CIsK+Q0M+J64y93jISIfbZ63WE88vZmWFhYiIiomVhWyKkSxvfCS3eNhsmhsLyz/gh+lb4RNRarvmBEROS2WFY0E5GJIpJaWFioO4rT3D6uB/42eSzMDo1l7aZjuH/VBlTVsLAQEVHTsKxoppRaq5RKCg8P1x3FqW4e3Q1/nzIWPg6F5eOtJ/DLt3JRWWPRmIyIiNwNywq1ml+M7Io3p0fBz2z/x+yLn04haVkuKqpZWIiI6OKwrFCrunZYZ6TOjIKfj/0ftW935WHe0nUor2JhISKiC2NZoVZ35eBOWDx7PAJ87f+4/bDnNGYvzkZpZY3GZERE5A5YVqhNXDogAkvmxCDIz2zMsvYXYOaibBRVVGtMRkREro5lhdpMXL+OWD4vBqH+PsYs9+AZzFiQhcIyFhYiImoYywq1qajeHbBifizCAuyFZdORQkxdkImC0iqNyYiIyFWxrFCbG92zHVYmxqFdkK8x23asCFPTMpFfUqkxGRERuSKWFdJiRPdwrE6KQ8dgP2O240QxJqdm4lRRhcZkRETkalhWSJshXcKwOikOkaH+xmzPqRIkpGbieGG5xmRERORKWFZIq4GdQ7EmOR5dwwOM2f78UiSkZOLImTKNyYiIyFWwrJB2fSOCkZ4Uj+7tAo3ZoYIyJKRk4uDpUo3JiIjIFbCsOIGIjBeR10Rki4iUiMgxEflQRKJ1Z3MXvToGIT05Dr06BBmzo2fLkZCSib15JRqTERGRbiwrzvFbAJMAfAPgIQCvABgMIEtEbtCYy630aF9bWPpFBBuzE0UVSEjJxO6TxRqTERGRTqKU0p3B7YnIJQBylFJVDrP2ALYBOKGUGneh94iOjlY5OTmtmNJ9nCqqwLQFWdh9yr5HpWOwH1bMj8XQrmEakxERkTOJSK5S6oJHIbhnxQmUUj86FhXb7AyArwEM05PKfXUKC8CqpDgM6RJqzE6XVmFKWia2Hi3UmIyIiHRgWWld3QDk6w7hjiJC/LEqMQ4jutv3pJwtq8bUtExsPHxWYzIiImprHldWRCRURG4WkWdE5GMRyRcRZbsNucj36CIir4rIXhGpEJGTIrJWRK5pQo5LAVwBIL25fxZv1z7YD2/Nj8Ponu2MWVFFDaYvyELOgQKNyYiIqC15XFkBcA2ADwA8AeB6AB2b8mIRGQVgK4AHAPQDUAkgAsBNAD4Xkccu4j06AVgJ4CCAp5vy+VRfeKAvVsyLQVTv9saspLIGMxdlI3PfaY3JiIiorXhiWQGAUwD+g9qikHSxLxKRQAD/Rm3B2QBghFIqHEB7AC8BEADPich153mPUNtnhwKYqJQqau4fgmqFBvhi2dwYxPbtYMzKqiyYvTgb3+/mUTYiIk/niWVlrVKqs1LqRqXUHwF83oTXJgPoDaAEtUVjGwAopYqUUg8DeN+23XMNvdhWdtYCGArgJqXU1mb+Gegcwf4+WDInBpcNiDBmFdVWzF26Dl/vPKUxGRERtTaPKytKKUsLXj7N9nOlUupoA8+/aPs57tzzX0TEF8DbAC4BcIdS6scW5KAGBPqZsWBWNK4cHGnMqmqsSF6Wi8+3n9SYjIiIWpPHlZXmsh2+ibI9/LSRzTIB1F07e7XDa00A3kLtOTIzlVKftFZObxfga0bKjChcO7SzMauyWHHPilx8vOW4xmRERNRaWFbshqL2nBSgdjG3/6GUsgLYaXvouH7KXwHcBeBLAD4iMt3x1lqBvZW/jxlvTBuHX4zoYsxqrAr3rdqADzY2tEOMiIjcmY/uAC6kq8P9Y+fZru45x+3rVqidYLuda0ULclED/HxMeH3KWPx6zSb8e1Pt/yUWq8JD6RtRbVG4M6qH5oREROQs3LNiF+xwv/w825XZfobUDZRSVyqlpLFbY28kIkkikiMiOXl5eS2M7318zCa8kjAGd4yzFxOrAh55exNWZx/SmIyIiJyJZcWu0VLRWpRSqUqpaKVUdGRk5IVfQP/DbBK8eOcoTInpacyUAh57dwuWZRzQlouIiJyHZcWuxOF+4Hm2C2pge9LIZBL8360jMTO+d735Ux9sw4L/7tOUioiInIVlxc7xPJVu59mu7jmnXHoiIhNFJLWwkF/Q1xImk+Dpm4dj3mV9683//NFP+Oc3ezWlIiIiZ2BZsdsBQNnuD29oA9slyoNtD7c740OVUmuVUknh4eHOeDuvJiJ44sahuOfK/vXmf/lkB177cremVERE1FIsKzZKqWIAObaHDV3RAwCxAOpaxZetHoqaTETw6M8H48FrBtabv/z5Lvz1051QSjXySiIiclUsK/WttP2cJiJdG3j+YdvPXKXUzgaeJxcgInhowiA8fN2gevO/f70Hz328g4WFiMjNeGRZEZGIuhtqv4SwTjvH52yHdRyloPabkkMBfCgiw2zvFyoiLwC43bbd4639Z6CWu+/qgXj8hnrfioDU7/bh6bXbWViIiNyIpy4K19iiJRnnPO4L4EDdA6VUuYjcgtpDPOMAbBORItSuqWJC7TktjyulPnNWUBGZCGDigAEDnPWW5CDpZ/3hazbh6bX2U4yW/HgA1RYrnrllBEymNr9inYiImsgj96y0hFJqE4ARAF4DsA+AP4DTAD4CMEEp9byTP48n2LayOZf2xf/dNqLe7K2sQ3js3c2wWLmHhYjI1XnknpXzrRp7ka8/AeBB2408wLTY3vA1mfDbdzej7gjQmpwjqLYovHjnKPiY2duJiFwVf0OT15g0videnjQajkd+3ttwFL9K34hqi1VfMCIiOi+WFfIqt43tgVcnj4XZobF8uPk47lu5HlU1LCxERK6IZUUzrmDb9iaO7oZ/TB0HX7O9sHy67STuWZGLyhqLxmRERNQQlhXNeIKtHteP6II3p0fBz+FclS93nELislxUVLOwEBG5EpYV8lrXDO2MtFnR8Pex/2vw3a48zF2yDmVVNRqTERGRI5YV8mpXDIrE4tnjEehrNmY/7j2N2YvWoaSShYWIyBWwrJDXu2RABJbOjUGwn72wZB8owMyFWSiqqNaYjIiIAJYVIgBATN8OWDYvFqH+9qWH1h86ixkLslBYxsJCRKQTy4pmvBrIdUT1bo8V82MRFmAvLJuOFGJKWiYKSqs0JiMi8m4sK5rxaiDXMrpnO6xMjEP7IF9jtv14EaakZiKvuFJjMiIi78WyQnSOEd3DsTopHhEhfsZs58liTE7NwKmiCo3JiIi8E8sKUQMGdwnF6qQ4dAr1N2Z780qRkJqJ44XlGpMREXkflhWiRgzoFIr05Hh0DQ8wZvvzSzEpJQOHC8o0JiMi8i4sK0Tn0TciGGuS49G9XaAxO1xQjsmpmTh4ulRjMiIi78GyQnQBPTsEYc3d8ejdMciYHT1bjkkpGdibV6IxGRGRd2BZ0YyXLruH7u0CkZ4Uj34RwcbsZFElElIysetkscZkRESej2VFM1667D66hAdgdXIcBnYKMWb5JZWYnJqJ7ceKNCYjIvJsLCtETdApNACrk+IwpEuoMSsorcLUBZnYepR7x4iIWgPLClETdQzxx6rEOIzoHmbMzpZVY0paJjYcOqMxGRGRZ2JZIWqG9sF+eGt+HMb0bGfMiitqMGNhNtYdKNCYjIjI87CsEDVTeKAvls+Lwfg+7Y1ZSWUNZi3KRsbe0xqTERF5FpYVohYIDfDFkjkxiOvXwZiVVVkwZ0k2/rs7T2MyIiLPwbJC1ELB/j5YPDsGlw+MMGYV1VbMW5qDr3ec0piMiMgzsKxoxnVWPEOgnxlpM6Nx1eBIY1ZVY0XS8hx8tu2ExmRERO6PZUUzrrPiOQJ8zXhzRhQmDOtszKotCr98az0+2nxcYzIiIvfGskLkRP4+ZrwxbRxuHNnVmNVYFe5ftR4fbDyqMRkRkftiWSFyMl+zCa9OHoNbx3QzZlYF/Cp9I97OPaIxGRGRe2JZIWoFPmYTXpo0BndG9TBmSgGPvL0Jq7IPaUxGROR+WFaIWonZJHjhjlGYEtPLmCkF/O7dLViWcUBbLiIid8OyQtSKTCbBs7eNwKz43vXmT32wDQv+u09TKiIi98KyQtTKRAR/vHk4Ei/vW2/+549+whvf7NGUiojIfbCsELUBEcHjNwzFvVf1rzd/4ZOdePWL3VBKaUpGROT6WFaI2oiI4OHrBuOhawfVm7/yxS789bOdLCxERI1gWSFqQyKCB68diEevH1xv/o+v9+K5j3ewsBARNYBlRTMut++dfnnlADxx49B6s9Tv9uHptdtZWIiIzsGyohmX2/de8y/vh6dvHl5vtuTHA/j9+1thtbKwEBHVYVkh0mjWJX3w7G0jIWKfrcw6hN++sxkWFhYiIgAsK0TaTY3thRfuGFWvsPwr9wh+s2YjaixWfcGIiFwEywqRC7gruif+ljAGJofC8v7GY3hw9UZUs7AQkZdjWSFyEbeM6Y7Xp4yDj0Nj+WjLcdz71npU1lg0JiMi0otlhciF3DiqK96YNg6+Znth+Wz7Sdy9PBcV1SwsROSdWFaIXMx1w7sgdUY0/Hzs/3p+vTMPictyUF7FwkJE3odlhcgFXTWkExbOioa/Q2H57+58zFmSjdLKGo3JiIjaHssKkYu6fGAkFs8Zj0BfszHL3FeAWYuyUVxRrTEZEVHbYlkhcmGX9I/AsnkxCPH3MWY5B89gxsJsFJazsBCRd2BZIXJx4/t0wPJ5MQgNsBeWjYfPYvqCLJwtq9KYjIiobbCsELmBsb3aY+X8OIQH+hqzLUcLMSUtCwWlLCxE5NlYVojcxMge4ViVGIcOwX7G7KfjRZicmoG84kqNyYiIWhfLCpEbGdYtDKuT4hAR4m/Mdp0sweTUDJwsqtCYjIio9bCsaCYiE0UktbCwUHcUchODOociPTkOncPshWVvXikSUjJw7Gy5xmRERK2DZUUzpdRapVRSeHi47ijkRvpHhiA9KR7dwgOM2YHTZZiUkoHDBWUakxEROR/LCpGb6hMRjPTkePRoH2jMjpwpR0JKBg7kl2pMRkTkXCwrRG6sZ4cgrEmOR5+OQcbsWGEFElIzsOdUicZkRETOw7JC5Oa6tQtEenI8+kUGG7OTRZWYnJqJnSeKNSYjInIOlhUiD9A5LADpSfEY1DnEmOWXVGJKWia2HePJ20Tk3lhWiDxEZKg/ViXGYWjXMGNWUFqFqWlZ2HzkrMZkREQtw7JC5EE6hvhjVWIsRna3X11WWF6NaWlZWH/ojMZkRETNx7JC5GHaBflhxfxYjO3VzpgVV9ZgxoIsrDtQoDEZEVHzsKwQeaDwQF8snxeL8X3aG7PSKgtmLszGj3vzNSYjImo6lhUiDxXi74Olc2MQ36+jMSuvtmDO4nX4bleexmRERE3DskLkwYL8fLBo9nhcPjDCmFXWWDF/aQ6+2nFSYzIioovHskLk4QL9zEibGY2rh3QyZlUWK5KX5+KTrSc0JiMiujgsK04gIiEi8rSI/EdE8kREichjunMR1QnwNePN6VH4+fDOxqzaonDvyvX4cPMxjcmIiC6MZcU5IgA8BWAkgA2asxA1yM/HhL9PHYcbR3U1ZharwgOrNuC9DUc0JiMiOj8f3QE8xHEA3ZVSx0SkD4D9euMQNczXbMKrCWPgZzbhvQ1HAQBWBfx6zSZUWxQmRffUnJCI6H9xz4oTKKUqlVLcl05uwcdswl/vGo1J0T2MmVLAo29vxltZBzUmIyJqGMsKkRcymwTP3z4K02J71Zv//r2tWPIDdwwSkWvxuLIiIqEicrOIPCMiH4tIvu2EVyUiQy7yPbqIyKsisldEKkTkpIisFZFrWjs/UVsxmQR/vnUEZl/Sp978j2u3I/W7vXpCERE1wBPPWbkGwHvNfbGIjALwFYC6lbSKUHsC7U0AbhSRx5VSz7c4JZELEBH8YeIw+PmYkPrdPmP+7H921F4tdNUAjemIiGp53J4Vm1MA/gPgaQBJF/siEQkE8G/UFpUNAEYopcIBtAfwEgAB8JyIXOf0xESaiAh+94shuO+cYvLipzvxyue7oJTSlIyIqJYn7llZq5R6v+6B7eqci5UMoDeAEgATlVJHAUApVQTgYRHpD+BWAM8B+MxZgYl0ExE8/PPB8PMx4eXPdxnzV7/cjSqLFY/+fDBERGNCIvJmHrdnRSllacHLp9l+rqwrKud40fZz3MWe/0LkTh64ZiAevX5wvdk/v9mLP3/0E/ewEJE2HldWmktEQgFE2R5+2shmmQAKbfevbvVQRBr88soBeOLGofVmC7/fjz/8exusVhYWImp7nngYqLmGovacFADY1tAGSimriOwEEANgmONzInIfgHa2GwBcJSJ1//u+rpQqBJGbmH95P/j5mPDUB/Z/FZZlHERVjRXP3jYSJhMPCRFR22FZsevqcP98C7zVPdf1nPnDqD3fpc51thsArIB9jwyRW5gZ3we+ZhMef28L6o4ArV53GNUWhRfuHAUzCwsRtRGWFbtgh/vl59muzPYzxHGolOrT1A8UkSTYrlbq1avXBbYmantTYnrB12zCo29vQt0RoHfWH0GN1YqX7hoNHzOPJBNR6+NvGrs2/89EpVSqUipaKRUdGRnZ1h9PdFHujOqBVxLG1NuT8sHGY3hg9QZUW6wakxGRt2BZsStxuB94nu2CGtieyKPdMqY7Xp8yFj4OheU/W07gnhXrUVnTkgvwiIgujGXFzvE8lW7n2a7uueOtmIXI5dwwsivemDYOvmZ7Yfnip5NIXp6LimoWFiJqPSwrdjsA1F2XObyhDUTEBKBuEYrtzvhQEZkoIqmFhTz/llzfdcO7IHVmNPx87L86vtmZh/lLc1BexcJCRK2DZcVGKVUMIMf2cEIjm8UCCLfd/9JJn7tWKZUUHh5+4Y2JXMBVgzth0azxCPC1//r4fk8+5izJRmlljcZkROSpWFbqW2n7OU1Ezr00Gai9PBkAcpVSO9soE5HLuWxgBBbPjkGQn9mYZe4rwKxF2SiuqNaYjIg8kUeWFRGJqLuh9ksI67RzfM52WMdRCoCDAEIBfCgiw2zvFyoiLwC43bbd4639ZyBydfH9O2LZ3BiE+NtXQMg5eAbTF2ajsIyFhYicxyPLCoA8h9t6h3nGOc/VW9xEKVUO4BYApwGMA7BNRAoBnAXwCGrPafmdUsppX2LIc1bInUX36YDl82IQGmAvLJsOn8W0hZk4U1qlMRkReRJPLSvNppTaBGAEgNcA7APgj9ry8hGACUqp5538eTxnhdza2F7tsSoxDu2CfI3Z1qNFmJKWifySSo3JiMhTCL9J1TVER0ernJycC29I5KJ+Ol6EaQuyUOCwR2VgpxC8NT8WncICNCYjIlclIrlKqegLbcc9K0TkFEO7hmF1UhwiQvyN2e5TJZicmokThRUakxGRu2NZISKnGdQ5FOnJcegcZi8s+/JLkZCagaNnz/eVW0REjWNZ0Ywn2JKn6R8ZgjXJ8ejezv6tFQdPl2HSmxk4dLrsPK8kImoYy4pmPMGWPFHvjsFIT45Dzw72wnL0bDkSUjOwP79UYzIickcsK0TUKnq0D8Ka5Hj0jQg2ZscLK5CQkoE9p4o1JiMid8OyQkStpmt4INKT4tA/0l5YThVXYnJqJnaeYGEhoovDskJErapTWABWJ8VjcOdQY5ZfUoXJqRnYepTnahHRhbGsEFGriwz1x6qkOAzrGmbMzpRVY2paJjYdPqsxGRG5A5YVzXg1EHmLDsF+WJUYh9E97CeTF1XUYPqCLOQePKMxGRG5OpYVzXg1EHmT8CBfLJ8fi6je9u8XLa6swcyFWcjad1pjMiJyZSwrRNSmwgJ8sXRuDGL7djBmpVUWzF68Dj/sydeYjIhcFcsKEbW5EH8fLJ4zHpcO6GjMyqstmLtkHb7dlacxGRG5IpYVItIiyM8HC2eNxxWDIo1ZZY0ViUtz8OVPJzUmIyJXw7JCRNoE+JqROjMK1w7tZMyqLFbcvSIXn2w9oTEZEbkSlhUi0srfx4w3pkXh+uFdjFm1ReHeleuxdtMxjcmIyFWwrGjGS5eJAD8fE16fOhYTR3czZharwoOrN+Dd9Uc0JiMiV9DssiIiD9hu3S68NTWGly4T1fI1m/C3hDG4fWx3Y2ZVwG/+tQnp6w5pTEZEuvm04LWvALAAeNNJWYjIy5lNghfvGg1fswnpOYcBAEoBv31nC6osCjPiemtOSEQ6tOQwUD6AYqVUlbPCEBGZTYLnbh+J6XG96s2ffH8rFn2/X1MqItKpJWVlPYBwEYm84JZERE1gMgmeuWUE5lzap978Tx9uR8q3e/WEIiJtWlJWXrO9/kknZSEiMogInrppGJKv6Fdv/tzHO/D6l7s1pSIiHZpdVpRSHwN4GMDdIrJcREY7LxYRUW1heez6IXjg6gH15i99vgsvf7YTSilNyYioLTX7BFsR2We7WwNgKoCpIlIO4DRqT7xtiFJK9W/uZxKR9xER/Pq6wfA1m/DS57uM+Wtf7UGlxYrHrh8CEdGYkIhaW0uuBurTwCzIdmsM/zPoHCIyEcDEAQMGXHBbIm92/zUD4etjwvMf7zBmKd/uQ3WNwpM3DWVhIfJgLSkrVzkthRdTSq0FsDY6OjpRdxYiV3f3Ff3hazbhmQ+3G7NFP+xHtcWKp28eDpOJhYXIEzW7rCilvnVmECKiizHvsr7w8zHhyfe3GrPlmQdRbbHi2dtGsrAQeSAut09EbmdGXG/85Y6RcDzys3rdYTz89iZYrDzaTORpWnIYqB6pPWA8GEDduit5AHYqnq5PRK0gYXwv+JpNePhfm1DXT95dfxTVFoWXJ9WugktEnqHF/zaLyAARWQKgEMA2AN/YbtsAFIrIYhHh2aNE5HS3j+uBv00eC7PDoZ+1m47h/pUbUFVj1ZiMiJypRWVFRG4GsAHADAAhAOScWwiAmQA2iMhNLYtKRPS/bh7dDX+fMhY+DoXlk20n8Mu3clFZ09gqCkTkTlryrcv9AawGEAxgH4BkAAMBBAIIsN2/G8Be2zZrbK8hInKqX4zsijenR8HP4dDPFz+dQtKyXFRUs7AQubuW7Fl5FLWl5GsAo5RSaUqpvUq843ZOAAAgAElEQVSpSqVUle1+KoDRAL4F4A/gkZZHJiL6X9cO64zUmVHw87H/Wvt2Vx7mLV2HsqoajcmIqKVaUlYmoHaRt2SlVHljG9meS0btYaHrWvB5RETndeXgTlg8ezwCfO2/2n7YcxqzF69DSSULC5G7aklZ6QqgUCm150IbKqV2AThrew0RUau5dEAElsyJQZCf2Zhl7y/AzIVZKKqo1piMiJqrJWWlDECQiPheaEMR8UPteSuN7oHxViIyUURSCwsLdUch8hhx/Tpi+bwYhPrbV2dYf+gsZizIQmEZCwuRu2lJWdkCwBfArIvYdpZt280t+DyPpJRaq5RKCg8P1x2FyKNE9e6A5fNjERZgLyybjhRi6oJMFJRWaUxGRE3VkrKyHLXnobwmIvOlgW8RE5EAEXkAwGuoPb9laQs+j4ioScb0bIeViXFoF2TfAbztWBGmpmUiv6RSYzIiagpp7gKztnLyCewn2p4A8F8AR1F75U9vALEAOqK21HwG4Bdc0bZh0dHRKicnR3cMIo+040QRpqVl4bTDHpUBnUKwcn4sOoUFaExG5N1EJFcpFX2h7Zq9Z8VWOm4FkIrastIVwCQAvwJwD4AbAUTYnnsTwG0sKkSkw5AuYVidFIfIUH9jtudUCRJSM3G8kKfSEbm6Zu9ZqfcmIr0A3A5gHOp/N9B6AO8qpQ61+EM8HPesELW+fXklmJqWhRNFFcasV4cgrEyMRY/2QRqTEXmni92z4pSyQi3HskLUNg6dLsOUtEwcPWvfo9K9XSBWJsaid8dgjcmIvE+rHwYSkfUikisi/Zr7HkREba1XxyCkJ8ehVwf7npSjZ8uRkJKJfXklGpMRUWNacjXQMAADlVL7nBWGiKgt9GhfW1j6Rdj3pJwoqkBCaiZ2nyzWmIyIGtKSsnIUtVf5EBG5na7hgVidFIeBnUKMWV5xJSanZuKn40UakxHRuVpSVj5F7Qq2sc4KQ0TUljqFBWBVUhyGdAk1ZqdLqzAlLRNbj3JVaSJX0ZKy8mcApwG8KSIRTspDRNSmIkL8sSoxDiO6hxmzs2XVmJqWiY2Hz2pMRkR1WrIo3M8ADAbwEoBqAMsAZKD2kmVLY69TSn3XrA/0cLwaiEivwvJqzFyUjU0OBSXE3wdL5oxHdJ8OGpMRea5Wv3RZRKyoXfANqD135WLeSCmlfC68mfdhWSHSr7iiGrMXr0PuwTPGLMjPjEWzxyOuX0eNyYg8U6tfugzgkMPt4DmPG7sdbsHnERG1qtAAXyybG4PYvvY9KWVVFsxenI3vd+drTEbk3bgonGYiMhHAxAEDBiTu3r1bdxwiAlBeZUHishx8v8deUPx8TEiZEYWrBnfSmIzIs7TFnhVyAqXUWqVUUnh4uO4oRGQT6GfGglnRuHJwpDGrqrEieVkuPt9+UmMyIu/UkhVsz4jIaa5gS0SeKMDXjJQZUZgwrLMxq7JYcc+KXHy85bjGZETepyV7VvwAmLmCLRF5Kn8fM96YNg43jOxizGqsCvet2oAPNh7VmIzIu7T0BFs/ZwUhInJFvmYTXps8FjeP7mbMLFaFh9I34p3cIxqTEXmPlpSVfwPwF5EJzgpDROSKfMwmvJIwBneM62HMrAp4+O1NSF93SGMyIu/QkrLyLIADANJEZKhz4hARuSazSfDinaMwJaanMVMK+O07W7A844C2XETeoCULtN0C4J8AngKwQUQ+xsWtYLusBZ9JRKSNyST4v1tHwtdswrKMg8b8yQ+2ocqiMO+yvhrTEXmulpSVJahdtbbum5dvtt0uhGWFiNyWySR4+ubh8DObsOD7/cb8mQ+3o6rGinuu7K8xHZFnaklZ+Q4Xt8Q+EZFHERH8/sah8PUx4Z/f7DXmf/lkB6otVjxwzUCN6Yg8T7PLilLqSifmICJyKyKCR38+GH5mE1790r769Muf70K1xYpfTxgEETnPOxDRxeIKtkREzSQieGjCIDzy88H15q9/tQfPf7wD/DoTIue46LIiIg+IyLxGngsRkbALvP4VEVnY1IBERK7u3qsG4Pc31L8oMuW7fXh67XYWFiInaMqelb8B+FMjz+0GUHCB108GMLsJn0dE5DYSf9YPf5w4rN5syY8H8MT7W2G1srAQtURTDwOd7wCsVx+cFRFfEfmTiBwSkQoR2SwiU3XnIqK2M/vSvvi/20bUm72VdQiPvbsZFhYWombjOSvOkwrg9wDeB3A/gMMA3hKRmVpTEVGbmhbbGy/cOQqO59auyTmCR/61CTUWq75gRG6MZcUJRGQcag9x/VEp9YBSKg3ATQC+AfCiiPhrjEdEbWxSdE+8PGk0TA6F5d0NR/Gr9I2oZmEhajKWFeeYBMAK4B91A1V7Vt3fAXQCcKWeWESky21je+C1KWNhdmgsH24+jvtWrkdVDQsLUVOwrDjHOAB7lVLnnmSc5fA8EXmZm0Z1wz+mjoOv2V5YPt12EvesyEVlTaPfSkJE5/C4siIioSJys4g8IyIfi0i+iCjbbchFvkcXEXlVRPbaTpY9KSJrReSaRl7SDcDxBubHHJ4nIi90/YgueHN6FPzM9l+3X+44hcRluaioZmEhuhgeV1YAXAPgAwBPALgeQMemvFhERgHYCuABAP0AVAKIQO05KJ+LyGMNvCzQtl09SikrgGrb80Tkpa4Z2hlps6Lh72P/lfvdrjzMXbIOZVU1GpMRuYemLrffQUS+amgOAI08V2+bNnIKQA6AdQCOovZKnQsSkUAA/0ZtwdkAYIZSapttwbunAPwGwHMisl4p9ZnDS8sB/M9JtCJiAuBre56IvNgVgyKxePZ4zFuag3LbHpUf957G7EXrsGjOeIT4t+Sr2og8m1zs6ooi4owzwpRSyuyE92mUiJiVUhaHx30A1H016lCl1I7zvPZXAF4BUAJgiFLq6DnPvwfgVgDrlVJRDvPPAfRWSg06Z/seqL2E+XGl1HPnyx0dHa1ycnIu/AckIreWvb8AcxZno7TKfghoXK92WDI3BmEBvhqTEbU9EclVSkVfaLumVPmlLcjTZhyLSjNMs/1ceW5RsXkRtWVlnIgMcSg+uQCuFpEO55xkG2v7ub4FmYjIg8T07YBl82Ixe1E2iitrDwGtP3QWMxZkYdncWIQHsbAQneui96y4q4vdsyIioQAKUbsS7x1KqXcb2MaE2q8VCAdwr1LqDds8CrWHnZ5USv3ZNhMAXwEYDqCnUup/zmlxxD0rRN5l0+GzmLEwC0UV9nNWhnUNw4r5segQ7KcxGVHbudg9K554gm1zDYX9KwO2NbSB7YTZnbaHwxzmuQCWA3jadhXRfABrUbu+ym8vVFSIyPuM7tkOq5Li0N5hT8r240WYkpqJvGL+yiByxLJi19Xh/rFGt7I/1/Wc+XwAzwK4HbWLw/VG7Qm6ixt7IxFJEpEcEcnJy8trRmQicmfDu4VjdVI8IkLse1J2nizG5NQMnCyq0JiMyLWwrNgFO9w/39U7ZbafIY5DpVSVUupJpVRPpZS/UmqkUmrF+T5QKZWqlIpWSkVHRkY2MzYRubPBXUKxOikenULtFxTuzStFQkoGjp3lhYREAMuKI6/+1mgi0mdApxCkJ8eja3iAMTtwugwJqRk4XFB2nlcSeQeWFbsSh/vnW8QtqIHtiYhapG9EMNYkx6N7O/uvn8MF5ZicmomDp0s1JiPSj2XFzvE8lfMtj1/3XEPL6zeZiEwUkdTCwkJnvB0RubGeHYKw5u549O4YZMyOni3HpJQM7M3jfx+R92JZsdsBoO467uENbWC7dHmw7eF2Z3yoUmqtUiopPDzcGW9HRG6ue7tApCfFo1+k/TS6k0WVSEjJxK6TxRqTEenDsmKjlCpG7VopADChkc1iUbvGCgB82eqhiMgrdQkPwOqkOAzsZD+PP7+kEpNTM7H9WJHGZER6sKzUt9L2c5qInHtpMgA8bPuZq5Ta2cDzRERO0Sm0trAM6RJqzApKqzAlLRNbjvCwMXkXjywrIhJRdwPQ3uGpdo7P2Q7rOEoBcBBAKIAPRWSY7f1CReQF1K6hAgCPt/afgYioY4g/ViXGYUT3MGNWWF6NqQsyseHQGY3JiNqWR5YVAHkON8fv5ck457leji9SSpUDuAXAaQDjAGwTkUIAZwE8gtpzWn53zjcutwhPsCWi82kf7Ie35sdhTM92xqy4ogYzFmZj3YGC87ySyHN4allpNqXUJgAjALwGYB8Af9SWl48ATFBKPe/kz+MJtkR0XuGBvlg+Lwbj+9h3FJdU1mDWomxk7D2tMRlR2/D4LzJ0F/wiQyK6kNLKGsxfmoOMffaCEuBrQtrMaFw+kKtgk/vhFxkSEXmYYH8fLJo9HpcPjDBmFdVWzFuag693nNKYjKh1sawQEbmRQD8z0mZG46rB9j0pVTVWJC3PwWfbTmhMRtR6WFY04wm2RNRUAb5mvDkjChOGdTZm1RaFX761Hh9tdsri2kQuhWVFM55gS0TN4e9jxhvTxuHGkfYloWqsCvevWo8PNh7VmIzI+VhWiIjclK/ZhFcnj8GtY+xfZ2ZVwK/SN+JfOYc1JiNyLpYVIiI35mM24aVJY3BnVA9jphTwyNubsTLrkMZkRM7DskJE5ObMJsELd4zClJh661zi8fe2YOmPB/SEInIilhXNeIItETmDySR49rYRmBXfu978D//ehgX/3acpFZFzsKxoxhNsichZRAR/vHk4Ei/vW2/+549+wj++3qMpFVHLsawQEXkQEcHjNwzFvVf1rzd/8dOd+NsXu8BVy8kdsawQEXkYEcHD1w3GQ9cOqjf/2xe78eKnO1lYyO2wrBAReSARwYPXDsQjPx9cb/7GN3vx7H9+YmEht8KyQkTkwe69agCeuHFovVnaf/fj6bXbWVjIbbCsEBF5uPmX98PTNw+vN1vy4wE8/t5WWK0sLOT6WFY046XLRNQWZl3SB8/eNhIi9tmq7EN49J3NsLCwkItjWdGMly4TUVuZGtsLL9wxql5heTv3CH6zZiNqLFZ9wYgugGWFiMiL3BXdE69MGgOTQ2F5f+MxPLh6I6pZWMhFsawQEXmZW8d2x+tTxsHHobF8tOU47n1rPapqWFjI9bCsEBF5oRtHdcUb08bB12wvLJ9tP4m7V+SiotqiMRnR/2JZISLyUtcN74KUGVHw87H/VfDVjlNIXJaD8ioWFnIdLCtERF7s6iGdsWBmNPwdCst/d+dj7pJ1KKuq0ZiMyI5lhYjIy/1sUCQWzxmPQF+zMcvYdxqzFmWjpJKFhfRjWSEiIlzSPwJL58Yg2M9eWNYdOIMZC7NQWF6tMRkRy4p2XBSOiFxFTN8OWD4/FqEBPsZsw6GzmL4gC2fLqjQmI2/HsqIZF4UjIlcyrld7rJwfh/BAX2O25WghpqRloaCUhYX0YFkhIqJ6RvYIx6rEOHQI9jNmPx0vwuTUDOQVV2pMRt6KZYWIiP7HsG5hWJUYh4gQf2O262QJJqdm4GRRhcZk5I1YVoiIqEGDu4RidVIcOoXaC8vevFIkpGTg2NlyjcnI27CsEBFRowZ0CsGa5Hh0Cw8wZgdOl2FSSgYOF5RpTEbehGWFiIjOq09EMNKT49GjfaAxO3KmHAkpGTiQX6oxGXkLlhUiIrqgnh2CsCY5Hr07BhmzY4UVmJSSgT2nSjQmI2/AskJERBelW7tArEmOR7/IYGN2qrgSk1MzsfNEscZk5OlYVoiI6KJ1DgtAelI8BnUOMWb5JZWYkpaJbce4uCW1DpYVzbiCLRG5m8hQf6xKjMPQrmHGrKC0ClPTsrD5yFmNychTsaxoxhVsicgddQzxx6rEWIzqYf/dVVhejWlpWVh/6IzGZOSJWFaIiKhZ2gX5YcX8WIzt1c6YFVfWYMaCLGTvL9CYjDwNywoRETVbWIAvls+Lxfg+7Y1ZaZUFsxZl48e9+RqTkSdhWSEiohYJ8ffB0rkxiO/X0ZiVV1swZ/E6fLcrT2My8hQsK0RE1GJBfj5YNHs8Lh8YYcwqa6yYvzQHX+04qTEZeQKWFSIicopAPzPSZkbj6iGdjFmVxYrk5bn4ZOsJjcnI3bGsEBGR0wT4mvHm9Cj8fHhnY1ZtUbh35Xp8uPmYxmTkzlhWiIjIqfx8TPj71HG4cVRXY2axKjywagPe23BEYzJyVywrRETkdL5mE15NGIPbxnY3ZlYF/HrNJqzJOawxGbkjlhUiImoVPmYT/nrXaNwV1cOYKQU8+vZmvJV1UGMycjcsK0RE1GrMJsFf7hiFabG96s1//95WLPlhv6ZU5G5YVoiIqFWZTII/3zoCsy/pU2/+x7XbkfrdXj2hyK2wrBARUasTEfxh4jAk/6xfvfmz/9mBv3+1W1MqchcsK0RE1CZEBI/9Ygjuu2pAvflfP9uFlz/fBaWUpmTk6lhWiIiozYgIHv75YPx6wqB689e+3I0XPt3JwkINYlnRTEQmikhqYWGh7ihERG3mgWsG4rfXD6k3++c3e/Hnj35iYaH/wbKimVJqrVIqKTw8XHcUIqI2dc+V/fHkTcPqzRZ+vx9/+Pc2WK0sLGTHskJERNrMu6wvnrlleL3ZsoyDePy9LSwsZGBZISIirWbE98Hzt4+EiH22et1hPPL2ZlhYWAgsK0RE5AImx/TCi3eOhsmhsLyz/ggeSt+IGotVXzByCSwrRETkEu6M6oFXEsbA7NBY/r3pGO5ftQHVLCxejWWFiIhcxi1juuP1KWPh41BYPt56AvesWI/KGovGZKQTywoREbmUG0Z2xRvTxsHXbC8sX/x0EsnLc1FRzcLijVhWiIjI5Vw3vAtSZ0bDz8f+19Q3O/Mwf2kOyqtYWLwNywoREbmkqwZ3wqJZ4xHga/+r6vs9+ZizJBullTUak1FbY1khIiKXddnACCyZE4MgP7Mxy9xXgFmLslFcUa0xGbUllhUiInJpcf06YtncGIT4+xiznINnMGNhNgrLWVi8AcsKERG5vOg+HbB8XgxCA+yFZePhs5i2IBNnSqs0JqO2wLLiBCISIiJPi8h/RCRPRJSIPKY7FxGRJxnbqz1WJcahXZCvMdt6tAhT0jJxuqRSYzJqbSwrzhEB4CkAIwFs0JyFiMhjjegejlWJcegY7GfMdpwoxuTUTJwqrtCYjFoTy4pzHAfQXSnVE0CS7jBERJ5saNcwrE6KQ2SovzHbfaoEk1MycaKQhcUTsaw4gVKqUil1THcOIiJvMbBzKNKT4tAlLMCY7csvRUJqBo6eLdeYjFoDywoREbmlfpEhSE+OQ/d2gcbs4OkyJKRk4HBBmcZk5GwuW1ZEJFREbhaRZ0TkYxHJt524qkRkyEW+RxcReVVE9opIhYicFJG1InJNa+cnIqLW17tjMNKT49Czg72wHDlTjkkpGdifX6oxGTmTy5YVANcA+ADAEwCuB9CxKS8WkVEAtgJ4AEA/AJWoPRH2JgCf82odIiLP0KN9ENYkx6NvRLAxO15YgYSUDOw5VawxGTmLK5cVADgF4D8AnkYTTlwVkUAA/0ZtwdkAYIRSKhxAewAvARAAz4nIdee87kqHvTcXul3rpD8jERG1UNfwQKQnxaF/pL2wnCquxOTUTOw8wcLi7nwuvIk2a5VS79c9EJE+TXhtMoDeAEoATFRKHQUApVQRgIdFpD+AWwE8B+Azh9ftAnDPRX7GT03IQ0REraxTWABWJ8Vj+oIs7DxZW1DyS6owOTUDy+fFYkT3cM0JqblctqwopVrytZrTbD9X1hWVc7yI2rIyTkSGKKV22D7zGIA3W/C5RESkUWSoP1YlxWH6gixsP14EADhTVo2paZlYPi8Wo3u205yQmsPVDwM1mYiEAoiyPfy0kc0yARTa7l/d6qGIiKjNdAj2w6rEOIzuYd+TUlRRg+kLspB78IzGZNRcHldWAAxF7TkpALCtoQ2UUlYAO20PhznjQ0XkPhF5AsB9ttFVIvKE7cZ9j0REbSg8yBfL58diXC/7npTiyhrMXJiFrH2nNSaj5vDEstLV4f75Fmqre67rebZpiocBPAPgN7bH19keP4PaE3uJiKgNhQX4Ytm8WMT07WDMSqssmL14HX7Yk68xGTWVJ5aVYIf751vGsG7FoBBnfKhSqo9SShq5HWjoNSKSJCI5IpKTl5fnjBhEROQgxN8HS+aMxyX97atflFdbMHfJOny7i7933YUnlhW58CauQSmVqpSKVkpFR0ZG6o5DROSRgvx8sGj2ePxskP33bGWNFYlLc/DF9pMak9HF8sSyUuJwP7DRrYCgBrYnIiIPFOBrRuqMKFwzpJMxq7JYcfeKXHyy9bjGZHQxPLGsOJ6n0u0829U9x39KiYi8QICvGf+cHoXrh3cxZjVWhXtXbsDaTfwuWlfmiWVlBwBluz+8oQ1ExARgsO3h9rYI1RgRmSgiqYWFhRfemIiIWsTPx4TXp47FxNH2/5a1WBUeXL0B764/ojEZnY/HlRWlVDGAHNvDCY1sFgug7nLiL1s91HkopdYqpZLCw3l1MxFRW/A1m/C3hDG4fWx3Y2ZVwG/+tQlr1h3WmIwa43FlxWal7ec0EWno0uSHbT9zlVI7G3ieiIg8mNkkePGu0UiI7mnMlAIefWczlmce1JiMGuLSZUVEIupuqL9WSTvH52yHdRylADgIIBTAhyIyzPZ+oSLyAoDbbds93tp/BiIick1mk+C520dielyvevMn39+KRd/v15SKGuKy3w1k09hF8BnnPO4L4EDdA6VUuYjcgtpDPOMAbBORItSuqWJC7TktjyulPoNmIjIRwMQBAwbojkJE5HVMJsEzt4yAr9mExT8cMOZ/+nA7qi1WJF/RX184Mrj0npWWUEptAjACwGsA9gHwB3AawEcAJiilntcYz8BzVoiI9BIRPHXTMCRf0a/e/LmPd+D1L3drSkWOXHrPilKqRQu8KaVOAHjQdiMiImqQiOCx64fA32zCa1/tMeYvfb4L1RYrHpowCCJus+aox/HYPStERERNISL49XWD8ZsJg+rNX/tqD57/ZAeUUo28klobywoREZGD+68ZiN/9Yki9Wcq3+/CnD7ezsGjCsqIZF4UjInI9yVf0x1M3Das3W/zDATz5wVZYrSwsbY1lRTOeYEtE5JrmXtYXz9w6ot5sReYh/O7dLSwsbYxlhYiIqBEz4nrjL3eMhOO5tek5h/Hw25tgYWFpMywrRERE55Ewvhdeums0TA6F5d31R/Gr9I2otlj1BfMiLCtEREQXcPu4Hvjb5LEwOzSWtZuO4f6VG1BVw8LS2lhWiIiILsLNo7vh71PGwsehsHyy7QR++VYuKmssGpN5PpYVzXg1EBGR+/jFyK54c3oU/Mz2vz6/+OkUEpfloqKahaW1sKxoxquBiIjcy7XDOiN1ZhT8fex/hX63Kw9zl6xDWVWNxmSei2WFiIioia4c3AmLZo9HgK/9r9Ef957G7MXrUFLJwuJsLCtERETNcOmACCydE4NgP7Mxy95fgJkLs1BUUa0xmedhWSEiImqm2H4dsWxeDEL97d8LvP7QWcxYkIXCMhYWZ2FZISIiaoGo3h2wYn4swgLshWXTkUJMXZCJgtIqjck8B8sKERFRC43u2Q4rE+PQLsjXmG07VoSpaZnIL6nUmMwzsKxoxkuXiYg8w4ju4VidFIeOwX7GbMeJYkxOzcSpogqNydwfy4pmvHSZiMhzDOkShtVJcYgM9Tdme06VICE1E8cLyzUmc28sK0RERE40sHMo0pPi0CUswJjtzy/FpJQMHC4o05jMfbGsEBEROVm/yBCsSY5H93aBxuxwQTkmp2bi4OlSjcncE8sKERFRK+jVMQjpyXHo1SHImB09W46ElEzsyyvRmMz9sKwQERG1kh7tg7AmOR79IoKN2YmiCiSkZmL3yWKNydwLywoREVEr6hIegNVJcRjYKcSY5RVXYnJqJn46XqQxmftgWSEiImplncICsCopDkO6hBqz06VVmJKWia1HuXTFhbCsaMZ1VoiIvENEiD9WJcZhRPcwY3a2rBpT0zKx8fBZjclcH8uKZlxnhYjIe7QP9sNb8+Mwumc7Y1ZUUYPpC7KQe7BAYzLXxrJCRETUhsIDfbFiXgyie7c3ZiWVNZixMBuZ+05rTOa6WFaIiIjaWGiAL5bOjUFs3w7GrKzKgtmLs/HDnnyNyVwTywoREZEGwf4+WDInBpcNiDBmFdVWzF2yDt/sPKUxmethWSEiItIk0M+MBbOiceXgSGNWWWNF0rJcfLH9pMZkroVlhYiISKMAXzNSZkRhwrDOxqzKYsXdK3Lx8ZbjGpO5DpYVIiIizfx9zHhj2jjcMLKLMauxKty3agM+2HhUYzLXwLJCRETkAnzNJrw2eSxuHt3NmFmsCg+lb8TbuUc0JtOPZYWIiMhF+JhNeCVhDO4Y18OYWRXwyNubsDr7kMZkerGsEBERuRCzSfDinaMweXxPY6YU8Ni7W7A844C2XDqxrGjG5faJiOhcJpPg2dtGYkZc73rzJz/YhoXf79eUSh+WFc243D4RETXEZBL86ZbhmHdZ33rzZz7cjn9+s1dTKj1YVoiIiFyUiOCJG4fi7iv615v/5ZMdeO3L3ZpStT2WFSIiIhcmIvjt9YPxwDUD681f/nwX/vrpTiilNCVrOywrRERELk5E8OsJg/DwdYPqzf/+9R48//EOjy8sLCtERERu4r6rB+LxG4bUm6V8tw9/+nC7RxcWlhUiIiI3kvSz/vjDxGH1Zot/OIAn3t8Kq9UzCwvLChERkZuZc2lf/PnWEfVmb2UdwmPvbobFAwsLywoREZEbmh7XGy/cMQoi9tmanCN4+F+bUGOx6gvWClhWiIiI3NSk8T3x8qTRMDkUlvc2HMWv0jei2oMKC8sKERGRG7ttbA+8OnkszA6N5cPNx3HfyvWoqvGMwsKyQkRE5OYmju6Gf0wdC1+zvbB8uu0k7lmRi4pqi8ZkzsGyQkRE5AGuH9EVb06Pgp/Z/lf7lztOIWm5+xcWlhFMovgAABS6SURBVBUiIiIPcc3QzkibFQ1/H/tf79/tysPcJetQVlWjMVnLsKwQERF5kCsGRWLx7PEI9DUbsx/3nsbsRetQUumehYVlRTMRmSgiqYWFhbqjEBGRh7hkQASWzo1BsJ+9sGQfKMDMhVkoqqjWmKx5WFY0U0qtVUolhYeH645CREQeJKZvByybF4tQfx9jtv7QWUxfkIXCMvcqLCwrREREHiqqd3usmB+LsAB7Ydl8pBBT0jJRUFqlMVnTsKwQERF5sNE922FVUhzaB/kas+3HizAlNRN5xZUak108lhUiIiIPN7xbOFYnxSMixM+Y7TxZjMmpGThZVKEx2cVhWSEiIvIC/9/evYfJUZV5HP/+ciNXkhgCAQQSQQIhAkKQgLIGEHTRIOqjqCgGgcDuougu3tBHV10FdRVRcSGoi8hFWB9EQFS84GWXRAgCKsjFYAAJBJJA7jfIu3+c05mi0z3TMz3T3dPz+zxPPdVVdU736Tc1lberTtWZOmkM3597GDuO2W7rukVPr+XES+az5Nn1TWxZ15ysmJmZDRB77Tiaa844jJ3HDt+6bvHydZw4bz6PrVjXxJZ1zsmKmZnZADJlh1Fce8Zh7DpuxNZ1j61Yz4mXzGfxsrVNbFl1TlbMzMwGmN1eNJJrzzyMPSaM3LpuycoNnDhvPoueXtPEllXmZMXMzGwA2nXcCK6ZexgvmThq67qlqzZy4iULeHDp6ia2bFtOVszMzAaoSWOH8/25M9l7p9Fb1y1bs5G3z1vAfUtWNbFlL+RkxczMbADbccxwrj59JvvuvP3WdSvWbuIdly7gT39vjaFgnKyYmZkNcBNGb8fVpx/Ky3btGPpl5frNvPNbC7jr0Wea2LLEyUovkHSIpK9J+pOkNZKWSLpJ0oxmt83MzKwW40YO44rTDuXA3cZtXbd6w3O8+9u3c8fiFU1smZOV3vIR4G3Ar4EPAhcAU4HfSzquie0yMzOr2dgRQ/neqa/gkMnjt65bs/E53vOd25m/aHnT2qWIaNqHtwtJhwMLI2JTYd144F7gyYg4qKv3mDFjRixcuLAPW2lmZlabtRuf47TvLmT+wx0JyvChg7j05Bkc8dKJvfY5ku6MiC6vQvjMSi+IiNuKiUpe9wxwKzCtOa0yMzPrmVHbDeE7cw7hiJfusHXdhs1bOPW7C7n1/qca3h4nK31rF2BZsxthZmbWXSOGDebSk2dw5NSOMylbtgTPbWn8FZmWTVYkjZF0vKTPSvqJpGWSIk/71PgekyRdKGmRpA2Slkq6UdLRDWj/K4FXA9f09WeZmZn1heFDB3Pxuw/mmGk7MXiQ+Po7Xs4x03ZqeDtats+KpBOAH1bZvG9E3N9F/f2BXwET8qpVwGhSghbAuRFxfi81t/yzdwTuALYAB0REl0/WcZ8VMzNrVZuf38IfHnmGQ18yoevC3dAufVaeAm4GPg3MrbWSpBHADaRE5S5gekSMBcYDXwYEnCfp2LJ6swpnb7qaXlPls8fkNo8BZteSqJiZmbWyoYMH9Xqi0h1DmvbJXbsxIq4vLUia3I26ZwB7AGtICcPjADlxOEfSnsAJwHnALYV6DwL/VONn/KV8RU6SbgT2BY6JiD93o81mZmZWQcsmKxHxfB3VT8rzq0qJSpkvkZKVgyTtU7qkFBFLgIt78oGShgI/AA4Hjo+I23ryPmZmZvZCrX4ZqNvyZZiD8+LPqhRbAJQGPDiqFz5zEHAl8Drg5Ij4ab3vaWZmZknLnlmpw76kPimQHsq2jYjYIukB4BX0znNQ/hN4K/BzYIikd5V93hW98BlmZmYDUjsmKzsXXi/ppFxp286dlKlV6Qm1x+SpnJMVMzOzHmrHZGVU4fX6Tsqty/PR9X5gRMzqST1Jc8l3Oe2+++71NsPMzKwttV2fFTouAbW8iJgXETMiYsbEib031oKZmVk7acdkZU3h9YhOyo2sUN7MzMxaTDsmK8V+Krt0Uq607Yk+bIuZmZnVqR2TlftJj9MH2K9SgXyr8dS8eF8jGlWNpNmS5q1cubLrwmZmZgNQ2yUrEbEaKA2yU+nOHIBDgbH59S/7vFGdiIgbI2Lu2LFjuy5sZmY2ALVdspJdlecnSap0a/I5eX5nRDzQoDaZmZlZD7R0siJph9JEGoSwZFxxW76sU3QJ8AhpMMGbJE3L7zdG0heBN+dy5/b1dzAzM7P6tPpzVp6usn5+2fIUYHFpISLWS3oj6RLPQcC9klaRnqkyiNSn5dyIuIUmkzQbmL3XXns1uylmZmYtqaXPrNQjIu4BpgNfAx4GtgOWAz8mjYh8fhObt5X7rJiZmXWupc+sRERdD3iLiCeBs/NkZmZm/ZAioutS1uckPU3qZ9NbdgCW9eL7DVSOY/0cw/o5hvVzDOvXFzHcIyK6fIS7k5U2JWlhRMxodjv6O8exfo5h/RzD+jmG9WtmDNu2z4qZmZm1BycrZmZm1tKcrLSvec1uQJtwHOvnGNbPMayfY1i/psXQfVbMzMyspfnMipmZmbU0JytmZmbW0pysmJmZWUtzstJGJE2SdKGkRZI2SFoq6UZJRze7ba0gD2R5vKTPSvqJpGWSIk/71FBfkuZKmi/pWUmrJd0l6UOShjXiOzSbpN0lfSDvV49K2pjjcI+k86uMcl6sP0zShyXdLWlNjuP8HNe6nljdn0iakffDn0r6q6SVOZaPS/qRpBO6qO84lpE0WtJjhb/pOZ2UdfwASXMK8ao2remkfuOOiRHhqQ0mYH/SkwUjTyuB5/PrLcBHm93GZk/ACYX4lE/7dFF3KGlcqVL5jcC6wvLtwOhmf8c+jt9ueV+Ksv3sucLyCuDIKvW3BxYWyq7NcSwt3wgMafb3bFAsLy6L42pgfdm6HwBDHceaY/rVsvjN8X7YZczm5O+8CXiyyrSoSt2GHhN9ZqUNSBoB3ABMAO4CpkfEWGA88GVAwHmSjm1eK1vGU8DNwKeBud2o9x/AccAG0h/4SGAUMJv0H/QhwCW92dAWNDjPfwy8FXhR3s9GkmLzN9I+d72kSRXqXwocTIrXbNIo6CNJ8dwAvIH07zIQzAc+SIrHmIgYExEjgN2BL+UybwE+WqGu41hG0kHAWcDvayju+G3rtoiYVGXas0qdxh4Tm53Zeap/Aj5Ax6+zXSts/2Hefmez29rkOA0uW55MDWdWgEn5DzKA91fY/kY6zmDt3+zv2YfxGwsc0Mn2feg4O/Cpsm0vL8T6+Ap1z87b1gE7Nvu7NnsCvpfjsahsveO47XceBNxBOsNXjM+cCmUdvxd+3zn5+/66m/Uafkz0mZX2cFKeXxURj1fYXvqldlAtfTPaVUQ838OqbwG2I13y2OahSBHxI+BB0hmsd/a4gS0uIlZGxD2dbL8fWJAXDy7bXIrLAxFxQ4Xq80jxHQG8ud62toE78nyXsvWO47beB8wA/isi7uqirOPXOxp+THSy0s9JGkPHfww/q1JsAWmnAjiqzxvVfo7M899GxIYqZW7J84Ee3+V5PrhsfSmGt1BBRKwHfpcXB3oMAQ7P87+VrXccCyTtCnwWWAp8ooYqjl/vaPgx0clK/7cvKXsFuLdSgYjYAjyQF6c1olFtphSzivHN7svzfQfS3QRFkoYAr8yLfy6sF+kSEdQWwwG5j+a7WfaXdBFwYl79jcJ2x3FbXwfGAOdExMrOCjp+ndpP0r2S1uc7ev4s6QJJU6qUb/gxcUi9b2BNV7xVdEkn5UrbOr211CoqxayW+I7O0+o+bVFr+hfStewtwOWF9duTOt6B99EXkPRi4LEKmzYAn4+IbxbWOY4FkmYDbyL1t7iihiqOX3U7kG7QeIYUp/3ydIak0yLiqrLyDT8m+sxK/zeq8Hp9J+XW5fnoPmxLuyrFuJb4wgCMsaT9gc/nxW9ERPEXl/fR6p4nXcJYSrp9FFJH0fMonFXJHMdM0ihSfDaTkuRaOH7bWgJ8CpgODI+ICaTv/XrSmZERwOWS/qGsXsOPiU5W+r8BecmhSTzqZwX5QXDXk25dvBP4SHmRwmvHsCAinoh8iyjpP4appLNSnwbulrRfobjj2OEzpNu8L4iI+7oqnDl+ZSLiloj4TETcGxGb8rqNEXEzqd/UX0n9z86v9hYNaqqTlTZQfLrgiE7KjaxQ3mqzNs9HdlKmuG3AxFjSi0gd6aYADwGvr9DhrhiPWmI4YOJXFBFbIuLBiDgV+ArpP+MrJJWO044jIOlA0i3Gj5GSllo5ft2Q+wCVzpbOlDSxsLnhx0QnK/1f8Zph+W2OVNj2RB+2pV2VYlxLfNcwQA5yksaS7kCbDjwKvCYillYouoqOg5v30dp8Pc8PJD0bBBzHkgtJv/Y/Tuo3O7o4Fcptl9eV/tN0/Lqv9JA9kZ5LVdLwY6KTlf7vfjpOxe1XqUD+ZTY1L9Z6ytQ6lGJWMb5ZqXf8XyI/Famd5T4DN5Oeb/EkKVF5tFLZHI+/5MVaYuh9FIrPS9oTHMeCPfL8clKnzfKp5OK8fB84fj1U7dJZw4+JTlb6uYhYTRrnAuCYKsUOJT19FOCXfd6o9nNrnh8haXiVMqXYt3188/AON5KuaS8nJSoPdVGtFMOK+2iO6xF5se1jWIPiLaPFX6WOY30cv+55ReH1I4XXDT8mOllpD6Xbyk6qMurtOXl+Z0Q8UGG7de460iBd44DTyjfmWyinkn55XN3YpjVWHkn1OtJDoZ4Fji2786eaUlz2kfSGCttPJyXU60nDQ7QtSYNreO7Eh/L8OdI4QiUDPo4RMTkiVG0qFD0lr5tcWDfg41fS1T4oaXs6xqa6PSKeLmxu/DGxmeMSeOqdidSxdnHeMe4EpuX1Y4Av0jEWxrHNbmuzJ9LzBEpTcZyQmWXbBpXV+wIdY4a8mzzOEGkgr9Jo11c2+/v1cewGk0YCDtL1/5ndrH9NrrsMOK7wnifTMVrr55r9PRsQx8n57/S9wIsL6weR+qhcWdgvv+I4dju+XY267PjF1v1wAXAqsHth/TDgdcCfciyeB46qUL+hx0TlN7d+TtIBpNNtE/KqVaR72weRdppzI6La7WcDhqRad/gpEbG4UG8o6fbc4/KqjaQ/4lLnvTuAoyNdlmtL+VkLv8mLG+gYwqGSxyLikLL62wO/omN4iHWk/yS2y8s3AW+KiOd6rdEtSNJkXvgY/Q2kSz1j6IgFwGXA6eXxcBw7V/gbPyUiLquw3fGj6n64lvRQuKF53TrgzIj4XoX6DT0m+jJQm4g0wNx04GvAw6Q/vOXAj4FjnKjUJyI2k4Y+P5P0a2QjKQm8m/RckVe1c6KSFY8Xw4GdOpkmlleOiFWkfi4fBe4hxW8jKZ5nkEbBbev/ILIlpMfpzyPtPytJp9M3kzoufpu0P51SKR6OY30cv62WAu8HriUNx7KOdAlsHakf5BdIZ+m3SVSg8cdEn1kxMzOzluYzK2ZmZtbSnKyYmZlZS3OyYmZmZi3NyYqZmZm1NCcrZmZm1tKcrJiZmVlLc7JiZmZmLc3JipmZmbU0JytmZmbW0oY0uwFmZl2RNAR4F/B24ADSGFhrgSdJw0v8FvhVRNxRqHMgcAKwuNIYMWbWf/hx+2bW0iRNBG4GZhRWbyCNRbI9UBrqfmVEjCvUmwP8N/CbiJjVkMaaWZ/wZSAza3VXkBKV1cCHgZ0jYkROTMYCxwDfBJ5tXhPNrC/5MpCZtSxJ+wDH5sX3RsQPitvzqK6/AH4h6ZxGt8/MGsNnVsyslb2s8PqmzgpGxPrSa0lBugQE8GpJUTbNKq8v6VWSvi/p75I2Slou6ReS3iFJFcrPyu+1OC/PlnSrpGckrZE0X9I7e/CdzayMz6yYWX+xK7CoxrJLgRGkPi2bgRVl2zcVFyR9gXSJqWQ1MA44Ok/HSzopIrZU+jBJZwNfBQJYmT97JjBT0mER8b4a221mFfjMipm1sjsLry/KnW27FBGTgLPz4m0RMalsuq1UNicaHwaeBv4ZGB8R2wOjgLcBT5DuQvpIlY+bCHwJuJzUn2Y8sAPw5bz9LJ9hMauP7wYys5Ym6bvAyXlxE/A7YAFwBykRebpKvTl0cTeQpHHAY8Bw4JURcXuFMjOB20gdeCdFxKa8fhZway72c+C1UXZAlXQZ8B7gr8De5dvNrDY+s2Jmre504CukRGUY6bLMx4Hrgack3S7ppEr9SmrwFmA08L+VEhWAiFhAepbLeODgKu9zXpVE5HN5vhfp+TBm1gNOVsyspUXEpoj4N2A34EzgauAhUv8QgENItzdfI6m7x7TD8/xQSU9Wm4Ddc7ndKrzHZuD/qrT9IdJlJICDutk2M8vcwdbM+oWIeAq4JE9I2gmYDXySlES8lZQ0XNiNt905z0fkqSsjK6xbVro0VMXj+XNq6m9jZtvymRUz65ciYmlEfIt0xmJpXv3ebr5N6Rh4QUSohumyHjS1J5enzKzAyYqZ9WsRsQz4UV7cu5vVS0nOtDqasIOkYZ1sL529qdgR2My65mTFzNrB2jwvXo4pPROlszMb8/P81ZIm9PCzhwKHVdogaS9gl7z4hx6+v9mA52TFzFqWpCmS9uyizEjS6MoAdxc2rcrzcVT3P6REZzjpWSmdfc74TjZ/rMrdSB/L84eAezp7fzOrzsmKmbWy/YAHJF0n6W2SSpdUkDRK0mzSc1em5NXFzrX35vk0SYdWevOIWE5HQnGKpGslTS98xvD8GP6LqHLHD7AOOAr4tqQdc71x+am4pT40/+5nrJj1nB8KZ2YtS9JrgZ+WrV5PutwztrDueeCTEfH5svq/Af4hL64gPUYf4O35+Smlcp8APkPHJaN1wMb8GaUfdYsjYkqhzizSQ+EeIT1q/wLS7dTPltW7KCLOqvlLm9k2nKyYWUuTtDfpFuVXAdNJYwQNIyUeDwO/Bb4VEfdWqDuBlIT8Y6EewJER8euysi8DzgKOBF4MDCZ1iv0jcANwXb59ulR+FjlZiYjJ+SzPvwIvJ/Vj+SPwjYi4su4gmA1wTlbMzHqgPFlpbmvM2pv7rJiZmVlLc7JiZmZmLc3JipmZmbU0JytmZmbW0tzB1szMzFqaz6yYmZlZS3OyYmZmZi3NyYqZmZm1NCcrZmZm1tKcrJiZmVlL+3+IqeDNhUc/XQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"xs = gradient_descent(x0, [0.1]*50, quadratic_gradient)\n",
"error_plot([quadratic(x) for x in xs])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Constrained Optimization\n",
"\n",
"Let's say we want to optimize the function inside some affine subspace. Recall that affine subspaces are convex sets. Below we pick a random low dimensional affine subspace b+U and define the corresponding linear projection operator."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# U is an orthonormal basis of a random 100-dimensional subspace.\n",
"U = np.linalg.qr(np.random.normal(0, 1, (1000, 100)))[0]\n",
"b = np.random.normal(0, 1, 1000)\n",
"\n",
"def proj(x):\n",
" \"\"\"Projection of x onto an affine subspace\"\"\"\n",
" return None\n",
"# What is this???"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x0 = np.random.normal(0, 1, (1000))\n",
"xs = gradient_descent(x0, [0.1]*50, quadratic_gradient, proj)\n",
"# the optimal solution is the projection of the origin\n",
"x_opt = proj(0)\n",
"error_plot([quadratic(x) for x in xs])\n",
"plt.plot(range(len(xs)), [quadratic(x_opt)]*len(xs),\n",
" label='$\\\\frac{1}{2}|\\!|x_{\\mathrm{opt}}|\\!|^2$')\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The orangle line shows the optimal error, which the algorithm reaches quickly. The iterates also converge to the optimal solution in domain as the following plot shows."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"error_plot([np.linalg.norm(x_opt-x)**2 for x in xs])"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python [default]",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
![](\"images/l10-projected-gd.png\")
\n",
"