{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "# Interior Point Methods" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "This tour explores the use of [interior point methods](https://en.wikipedia.org/wiki/Interior-point_method) for constraint minimization under positivity constraints.\n", "\n", "The definite reference for this Numerical Tour is the book [\"Convex Optimization\"](https://web.stanford.edu/~boyd/cvxbook/) of Boyd and Vandenberghe, which is a must read.\n", "\n", "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n", "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n", "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n", "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n", "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n", "$\\newcommand{\\choice}[1]{ \\left\\{ \\begin{array}{l} #1 \\end{array} \\right. }$\n", "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n", "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n", "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n", "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n", "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n", "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n", "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n", "$\\newcommand{\\CC}{\\mathbb{C}}$\n", "$\\newcommand{\\RR}{\\mathbb{R}}$\n", "$\\newcommand{\\EE}{\\mathbb{E}}$\n", "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n", "$\\newcommand{\\Ww}{\\mathcal{W}}$\n", "$\\newcommand{\\Vv}{\\mathcal{V}}$\n", "$\\newcommand{\\Nn}{\\mathcal{N}}$\n", "$\\newcommand{\\NN}{\\mathcal{N}}$\n", "$\\newcommand{\\Hh}{\\mathcal{H}}$\n", "$\\newcommand{\\Bb}{\\mathcal{B}}$\n", "$\\newcommand{\\Ee}{\\mathcal{E}}$\n", "$\\newcommand{\\Cc}{\\mathcal{C}}$\n", "$\\newcommand{\\Gg}{\\mathcal{G}}$\n", "$\\newcommand{\\Ss}{\\mathcal{S}}$\n", "$\\newcommand{\\Pp}{\\mathcal{P}}$\n", "$\\newcommand{\\Ff}{\\mathcal{F}}$\n", "$\\newcommand{\\Xx}{\\mathcal{X}}$\n", "$\\newcommand{\\Mm}{\\mathcal{M}}$\n", "$\\newcommand{\\Ii}{\\mathcal{I}}$\n", "$\\newcommand{\\Dd}{\\mathcal{D}}$\n", "$\\newcommand{\\Ll}{\\mathcal{L}}$\n", "$\\newcommand{\\Tt}{\\mathcal{T}}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\al}{\\alpha}$\n", "$\\newcommand{\\la}{\\lambda}$\n", "$\\newcommand{\\ga}{\\gamma}$\n", "$\\newcommand{\\Ga}{\\Gamma}$\n", "$\\newcommand{\\La}{\\Lambda}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\Si}{\\Sigma}$\n", "$\\newcommand{\\be}{\\beta}$\n", "$\\newcommand{\\de}{\\delta}$\n", "$\\newcommand{\\De}{\\Delta}$\n", "$\\newcommand{\\phi}{\\varphi}$\n", "$\\newcommand{\\th}{\\theta}$\n", "$\\newcommand{\\om}{\\omega}$\n", "$\\newcommand{\\Om}{\\Omega}$\n", "$\\DeclareMathOperator{\\eqdef}{\\overset{\\tiny def}{=}}$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "import math as m\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Useful helpers." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "def dotp(x,y):\n", " return np.sum( x.flatten()*y.flatten() )\n", "np.random.seed(123) # to ensure reproductibitily " ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Problem formulation\n", "\n", "The goal is to solve problem of the form:\n", "$$\n", " (\\mathcal{S}_\\infty) \\qquad \\umin{x\\in \\RR^d, A x \\leq b} f(x)\n", "$$\n", "for $A \\in \\RR^{m \\times d}$.\n", "\n", "This can be generalized for instance by replacing $A x$ by a matrix and $\\leq$ by PSD matrix inequalities." ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Lasso Primal Problem\n", "\n", "The [Lasso problem](https://en.wikipedia.org/wiki/Lasso_(statistics)) (applied to the regression problem $Bw \\approx y$ for the design matrix $B \\in \\RR^{n \\times p}$)\n", "$$\n", "\t\t(\\mathcal{P}_\\la) \\qquad \n", " \\umin{w \\in \\RR^p} \\frac{1}{2}\\norm{Bw-y}^2 + \\la \\norm{w}_1\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Set the parameters $(n,p)$ of the Lasso problem and generate a random Gaussian matrix (so that the problem is a [compressed sensing](https://en.wikipedia.org/wiki/Compressed_sensing) problem, with $n
0] = np.log(w[v>0])\n", " return np.sum(w)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "The gradient and hessian of $f_t$ read\n", "$$\n", "\t\\nabla f_t(x) = \\nabla f(x) + \\frac{1}{t} A^\\top \\frac{1}{y-Ax}\n", "\t\\qandq\n", "\t\\partial^2 f_t(x) = \\partial^2 f(x) + \\frac{1}{t} A^\\top \\text{diag}\\pa{\\frac{1}{(y-Ax)^2}} A.\n", "$$" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "def ft(x,t):\n", " if t<0:\n", " return f(x)\n", " else:\n", " return f(x) - 1/t * Log(b - A@x)\n", "def nablaft(x,t):\n", " return nablaf(x) + 1/t * A.T @ (1/(b - A@x))\n", "def hessianft(x,t):\n", " return hessianf(x) + 1/t * A.T @ np.diag( 1/(b-(A@x))[:,0]**2 ) @ A" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Should be 0: 6.905096885932568e-07\n", "Should be 0: 2.9706266473930562e-09\n" ] } ], "source": [ "t = 10\n", "x = .5 + np.random.rand(d,1)\n", "u = .5*np.random.randn(d,1)\n", "m2 = (ft(x+tau*u,t)-ft(x,t))/tau\n", "m1 = dotp( u, nablaft(x,t) );\n", "print('Should be 0: ' + str(abs(m1-m2)/abs(m1)) )\n", "m2 = (nablaft(x+tau*u,t)-nablaft(x,t))/tau\n", "m1 = hessianft(x,t)@u;\n", "print('Should be 0: ' + str(np.linalg.norm(m1-m2)/np.linalg.norm(m1)) )" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Initialize the algorithm using a feasible point, i.e. here $x>0$." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "x = .01*np.ones((d,1))" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "For a fixed $t$, one can solve $(\\mathcal{S}_t)$ using [Newton method](https://en.wikipedia.org/wiki/Newton%27s_method) with some [line search](https://en.wikipedia.org/wiki/Line_search) procedure to select the step size $0<\\tau_\\ell\\leq 1$\n", "$$\n", "\t(N) \\qquad x_{k+1} \\eqdef x_k + \\tau_\\ell d_k \n", " \\quad\\text{where}\\quad\n", " d_k \\eqdef -[ \\partial^2 f_t(x_k)]^{-1} \\nabla f_t(x_k)\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "dk = - np.linalg.solve(hessianft(x,t),nablaft(x,t))" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "The use of a backtracking is extremely important in our case to ensure that the iterate stays within the constraints." ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "The simplest linesearch rule is [Armijo](https://en.wikipedia.org/wiki/Backtracking_line_search), which imposes that, for some $0 < \\al < 1/2$ one has\n", "$$\n", " \\text{(AC)}\\qquad\n", " \\phi_k(\\tau_k) \\eqdef f(x_k + \\tau_k d_k )\n", " <\n", " \\psi_k(\\tau_k) \\eqdef f(x_k) + \\alpha \\tau \\dotp{d_k}{\\nabla f(x_k)}. \n", "$$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4lGX2xvHvAQIhCAgYxQQwKLqgCIhZAV0r2LCABUEXKxZWlBXXsrv6s67Yy4quLDYUFQsidsVGsRuadCwghKAirgUVRHl+f5zJBlkgA8zMO+X+XFcuQzKTOa+BkzdPuR8LISAiIpmvRtQFiIhIYqihi4hkCTV0EZEsoYYuIpIl1NBFRLKEGrqISJZQQxcRyRJq6CIiWUINXUQkS9RK5YtttdVWoaSkJJUvKSKS8SZNmvRVCKGwuseltKGXlJRQVlaWypcUEcl4ZvZZPI/TkIuISJZQQxcRyRJq6CIiWSKlY+giIsm0atUqysvLWbFiRdSlbJL8/HyaNWtGXl7eJj1fDV1EskZ5eTn169enpKQEM4u6nI0SQmDZsmWUl5fTsmXLTfoaGnIRkayxYsUKmjRpknHNHMDMaNKkyWb9dqGGLiJZJRObeaXNrT0jGvqUCx5m4qn3wQ8/RF2KiEjaSvuGHgL8MmIkew/vxy9Ni2HgQJg5M+qyRETSTto3dDPYfuazHLP1RJ4LhxP+/W9o2xb23hsefhgydDZbRCTR0r6hAzTZyrhgzB/otfIhTjtoMeGGG+Hzz6FvX2jWDC68ED76KOoyRUR+Y/To0QwcODBlr5cRDR2gSxe49loY/txW3JF/AcydC6+8AvvvD7fdBjvtBN26wahRsGpV1OWKiDBlyhQ6duyYstfLqHXof/kLTJjg/+3SpQal3bp5E1+yBO67D4YNg169oGlT6NcPzjgDttsu6rJFJMfMmzePAQMG8O6779KkSRO++eYbzjvvvKS/roUQkv4ilUpLS8Pmpi1+/TXsthvUrAmTJ8OWW67xyV9/hZdegn//G55/3mdUDz0U+veH7t39SSKStWbPnk2bNm0AOO88mDo1sV+/QwcfENiQlStX0qlTJ0aMGEGPHj14++232XnnnamoqCA/P7/a11jzGiqZ2aQQQml1z82YIZdKjRvDY4/BokVw2mnes/+rZk047DB45hmYPx8uvRSmTIEjj4SWLeHqq6GiIrLaRST7vfLKK7Rv356ioiIaNGhA06ZNyc/P59dff+W0005L6mtn1JBLpc6d4frrfehlyBBfyfg/WrSAq66C//s/ePZZv2u/7DK48kpv8P37+3BNjYz7mSYicajuTjpZpk6dyq677sq0adNo164dX375JfXr18fMaNiwIW+88QYvvfQSV155ZVx37BsjY7vZoEFwxBFwwQXw/vsbeGBeHhx9NLz8sq+E+ctfYOJEOPhg2HFHuOEGWLo0ZXWLSHZr0KABc+bM4cMPP6Rdu3ZcfvnlDBgwgMmTJzN16lTmzp3L9ddfn/BmDhnc0M1g+HAoKoLeveE//4njSa1a+a19eTk88gg0bw4XXwzFxXD88TB+/FpjOCIiG6dv37589NFHXHXVVdx11100btyYc889lw8++IBOnTpRr169pL12XA3dzAaZ2Uwzm2FmI80s38y6mtlkM5tqZm+aWaukVbkelePp5eXrGE/fkDp1vIGPGwezZsHZZ/tk6n77wc47wz//GedPCBGR32rcuDHjx4+nefPmvP/++1xzzTWYGbNmzWLw4MHMmTOHiRMnJufFQwgbfAOKgflA3difHwdOAeYBbWIfOxsYXt3X2n333UMy3HJLCBDCrbduxhf54YcQ7r8/hE6d/Ivl54dw8skhvPNOCKtXJ6hSEUmmWbNmRV1CCCGEFStWhJYtW27Sc9d1DUBZqKa/hhDiHnKpBdQ1s1pAAVABBKBB7PMNYx+LxHnnQY8evmH03Xc38YsUFMApp/gXmDLF33/ySd/RtNtuMHQofP99AqsWkWxVp04dPv3005S/brUNPYSwGLgJWAgsAb4NIYwFTgdeMLNy4ETgumQWuiFmcP/9PiTeq1cC5jg7dIC77vIljkOH+gv86U+w7bZw1lne8EVE0ky1Dd3MGgE9gJZAEVDPzPoCg4DuIYRmwP3ALet5/plmVmZmZUuTuJqkUSPf9b90KZxwgu8x2mz163sDnzwZ3nsPjjsORoyAjh2hUyf/KfLjjwl4IRGRzRfPkEs3YH4IYWkIYRUwGtgLaB9CeC/2mMeAPdf15BDCsBBCaQihtLCwMCFFr0/HjnDnnfDqq3DFFQn8wmawxx4eL7B4sU+afv+9z8QWFSnSV0TSQjwNfSHQ2cwKzI/T6ArMAhqa2U6xxxwIzE5SjRulXz/vs//4h+/+T7hGjaoa+IQJvjO1MtJ3n318OeTKlUl4YRGRDYtnDP09YBQwGZgee84w4AzgSTObho+hX5jEOjfKHXf4PGbfvp4AkBRmVZns5eW+QamiAv74R0X6ikgk4lrlEkK4PITQOoTQNoRwYghhZQjhqRDCriGE9iGE/UIIqZ/SXY+6dX08HeCYY1JwBkZhoTfwefM80nfffeHWWz3S98ADfbWMIn1FJMkydqdodbbf3ucvp0yBc89N0YvWqFGVyb5woYeBzZ0Lxx7r2TKXXgqffZaiYkQk12RtQwc4/HC45BK45x6fz0ypoiJv4PPnw3PPQWkpDB7sqY+HH+4fS8hSHBERl9UNHTxcsWtX390fyfLxykjfZ5/15n7JJb4M8ogjFOkrkqU6dOjAF198waWXXsoDDzzAuHHj6NOnT9JfNyPjczdGzZowcqQvaTzmGJg0yReqRGK77byBX3aZN/ihQ6sifXv08Ejfrl0V6SuSCBGdcPHLL7/w9ddfs8022zBt2jSOPfZYJkyYQPv27RNbyzrkROcoLIQnnvDFKCedBKtXR1xQZaTv2LG+Eub8830J5EEH+USqIn1FMtacOXP+e+LQrFmz2Hnnnf+bjZ7sAy6qDXtJ5FuywrniNWSI52794x+RlrFuK1aE8MgjIeyzjxdZu3YIxx8fwvjxCgcTiVM6hHONGjUqnHvuuWHZsmWhXbt2IYQQ2rVrFxYsWBDOO++88Prrr4eLLroo/PTTT+t8firCubLCgAG+TPz//s/TctNKZaTv+PG+aal/f3jhBV8CucsuivQVyRC1a9dmzpw5lJWV0b59ex566CFKSkpYtGiRDrhIJDMYNgzatfPe+cknUVe0HpWZ7BUVnhfToIGPBxYVwamneq6MDuIQSUuHHHIIrVu35o9//CPjxo2jrKyMBx98MH0OuMgmBQXw1FM+79izJyxfHnVFG7BmpO/kyXDyyb7GvXNnn+VVpK9I2snLy+P222/n8MMP59577+W2226jYcOGKTngIucaOvhqwUcf9cOK+vXLkJvdykz2ykjfEDzSt6jIh2cSPZsvIpul8kzRSnfffTc1atTgmmuuYe+9907Ka+ZkQwffkX/ttfD443DTTVFXsxEqI32nTPE792OPhQce8IbfubMifUXSxKRJk9hmm21S+po529DB41d69YK//tUjWDKKWVUme0WFj7l/911VpO+f/+y/gohIzsjphm7mkQBt2kCfPklMZky2NSN9x4+H7t39xKVddvFVMiNHKtJXJAfkdEMH2GILGDPGNxsdfXSGj1aYVWWyl5fD9df7gRwnnOCRvhddBB9/HHWVIkkVMmJSbN02t/acb+gArVp5rPm0aXDmmRkySVqdrbf2Bj5vnu9I3WcfuOUW2HFH35GqSF/JQvn5+Sxbtiwjm3oIgWXLlm3W+nRL5YWXlpaGsrKylL3exrrmGg9IvPVWX/addSoqfIxp2DBYtAiaNoXTT4czzvB4X5EMt2rVKsrLy1mR9EMQkiM/P59mzZqRl5f3m4+b2aQQQml1z1dDX8Pq1R7g9eyzPkm6//5RV5Qkv/4KL77oyx9feMGHarp399Uzhx7qiWYikjbibegacllDjRq+AnDHHeG44/yMiqxUs2ZVJvv8+fD3v0NZmUf6br+9H8i6ZEnUVYrIRlJDX0uDBj5J+vPPnmj7ww9RV5RklZG+Cxf6LtSddvKwmxYtfI37K6+kQTyliMQjroZuZoPMbKaZzTCzkWaWb2YTzWxq7K3CzMYku9hU+d3vfKHItGkenZKB8ysbLy/Px5teecUjfQcNgnHjqiJ9b7xRkb4iaa7ahm5mxcBAoDSE0BaoCfQJIewdQugQQugAvAOMTm6pqXXYYXDddZ6jPnhw1NWkWKtWnsleXu7Lf4qKfMVMs2YeVzlxYo78lBPJLPEOudQC6ppZLaAA+O+ZaWZWHzgAyJo79EoXXuj969JL4emno64mAvn5voZ9wgSYMcMzY55/3pdA7rIL3H67In1F0ki1DT2EsBi4CVgILAG+DSGMXeMhRwGvhRC+S06J0TGDu++G3/8e+vaF6dOjrihClZnslUsf69f3eIHiYo8bUKSvSOTiGXJpBPQAWgJFQD0z67vGQ44HRm7g+WeaWZmZlS3NwDHYunU9brd+fZ8k/eqrqCuKWEFBVSb7pElw4omecFYZ6fvvfyvSVyQi8Qy5dAPmhxCWhhBW4WPlewKYWRNgD+D59T05hDAshFAaQigtLCxMRM0pV1zsTb2iwpczaoNlTGUDr6jw7JgQfFimqMijfRXpK5JS8TT0hUBnMyswMwO6ArNjn+sFPBdCyMxtWRuhUyffYPnGG74ARNbQoIE38ilT4J13fLnj8OFVkb7Dh2d4SI5IZohnDP09YBQwGZgee86w2Kf7sIHhlmxz0klwwQVw551+YyprMavKZK+ogNtug2+/9SGa4mLPU5g9u/qvIyKbRFv/N9Kvv/omy1dfhdde8wUfsgEh+CqZoUOrAsH22cfv6I8+2g/HFpEN0tb/JKlZ0+PFd9jB9+EsWBB1RWnOrCqTvTLSt7y8KtL34ovT+LRukcyihr4JttwSnnnGbzZ79Ejzg6bTSWWk70cfwcsvw957w803+0amgw6C0aM14yyyGdTQN9FOO8Fjj/l+mxNPVNzJRqlRo6qBL1wIV13lY+vHHOPZMpddlsXJaCLJo4a+GQ4+2M+MGDMGLrkk6moyVFGRh4HNn++/9uy2m6c9tmzp6Y/PP+8TFyJSLTX0zTRwoM/vXXedr86TTVSrVlUD//RT+Nvf4IMPfAZ6++399BFF+opskBr6ZjLzSJOuXf34ugkToq4oC5SU+F36okWejrbTTh6oUxnp++qrGuMSWQc19ATIy/O+s/32cNRROoc5YfLyqjLZ583zdezjxsGBB3rG8U03KYtBZA1q6AnSqJEfAAQ+cvDNN9HWk3V23NEz2cvL4aGHYNttPQ6zuFiRviIxaugJ1KqVL9z45BNlviRNfr438MpI37POqor0bdsWhgzRT1PJWWroCbbvvh4L8MorPmGqm8YkqsxkX7wY7r0X6tXz/+lFRR7p+/77+gZITlFDT4JTT/X9M0OH+g2jJFm9elUNfM1I306dYPfdPVVNkb6SA9TQk+Taa6FnT09mfPHFqKvJIWtG+v7rX76G/ayzqiJ9p02LukKRpFFDT5IaNXzurn176N3bh3slhRo0qMpkf+cd34U6fDh06ABdusADD8BPP0VdpUhCqaEnUb16vvlxiy18f8yXX0ZdUQ6qjPQdPtzH2m+91c9BPeUUv2tXpK9kETX0JGvWzJv6l196kJduCiPUuHFVAx83Dg45xIdldt4Z9tsPHn0UVq6MukqRTaaGngKlpTBihB/DedJJ2uQYubUjfa+7zsPAjj8emjeHv/7V4wdEMowaeoocc4zvixk1yvuFpImtt/ZM9o8/hpdegj/8wXeg7rCDp6899ZQ2FEjGUENPofPPh7PP9sZ+111RVyO/UaOGN/DRo+Gzz+DKK2HWLD9VqaTEI30XLYq6SpENiquhm9kgM5tpZjPMbKSZ5Zu7xszmmdlsMxuY7GIznRn8859w2GFwzjnwwgtRVyTrVFzsDbwy0rdDBw8LKymBI4/0b5wifSUNVXumqJkVA28CO4cQfjKzx4EXAAP2B04JIaw2s61DCBtcx5ENZ4omwvLlvlN93jyPINltt6grkmotWAB33+07Ur/4wg/iOOMM6NcPmjaNujrJcok+U7QWUNfMagEFQAXwJ+CqEMJqgOqauVTZYgsP8mrc2O/W9Zt8Bigp8Uz2hQs9WrNVK4/0bd4cevXyE8M12y0Rq7ahhxAWAzcBC4ElwLchhLHADkBvMyszsxfNbMfklppdioo8U+qHH6B7d/j226grkrjUrl2VyT53Lvz5z/D669CtG7RurUhfiVS1Dd3MGgE9gJZAEVDPzPoCdYAVsV8D7gbuW8/zz4w1/bKlS5cmrvIssOuu8OSTMGeO3+RpMUWG2Wknb+CLF/u24G22qYr07dsX3nxT4WCSUvEMuXQD5ocQloYQVgGjgT2BcuDJ2GOeAtqt68khhGEhhNIQQmlhYWEias4q3bpVpTP+6U/695+RKiN9J06E6dP96Kpnn4W99/af2nfcoUhfSYl4GvpCoLOZFZiZAV2B2cAY4IDYY/YF5iWnxOx32mk+HHvvvR7qJRmsMpO9osK/oQUFcO65PsbWr5+fk6qf2pIk1a5yATCzK4HewC/AFOB0oC7wMNACWA70DyFsMMpOq1zWLwT/Lf2RR+Dhh+GEE6KuSBJm0iT/NeyRR3zSpGNHT4A84QSfIRepRryrXOJq6Imihr5hK1f63pZ33oGxY313umSR777zn9Z33eVDM/Xr+0/x/v2h3TpHLEWAxC9blBSoU8d3mu+wgwd5TZ8edUWSUJWRvtOmwdtv+4ni993nGct77qlIX9lsauhpplEjjxSpVw8OPVRr1LOSWVUme0WFR/p+/bVH+hYX+6koc+ZEXaVkIDX0NNSihZ9y9P333tT/85+oK5KkWTPS94034KCD4M47oU0b2H9/eOwx+PnnqKuUDKGGnqbatYMxYzweoGdPWLEi6ookqcyqMtkXLfLlTp99Bn36eKi+In0lDmroaWz//eHBB2HCBD/3WDvLc8Q223gDr4z03WuvqkjfQw7xn/S//BJ1lZKG1NDTXJ8+cMstnqM+aJCWMOeUykjfp57ycLArrvDDaY86ysPBLr9ckyzyG2roGWDQIM9Sv/12v1GTHNSsmTfwBQvg6ad9ZczVV3toWI8eivQVQA09Y9x4o9+tX3SRL2WWHFWrVlUm+6ef+tDMu+96bOcOO8DgwfD551FXKRFRQ88QNWr4wfX77w+nnuphf5LjKiN9Fy2Cxx/3hn7JJR7pe9xxngKpMbqcooaeQSo3HrVu7SejTZ0adUWSFmrXrspknzMHBg7097t2hd/9Dm6+GZYti7pKSQE19AzTsKGvUW/UyNeoL1gQdUWSViob+OLFMGKEH4J9wQW+YenEE+Gtt3TXnsXU0DNQcbGvZqvMflHMvPyP/PyqTPYPP4TTT/fzUf/wh6pIX52qknXU0DNUmzYeub1woZ949P33UVckaauygVdUwD33QN26VZG+p58OCszLGmroGWyvvXx9+pQpvjR55cqoK5K0Vq9eVSb7Bx94fO/IkfD738Puu/sh2MuXR12lbAY19Ax32GEe2Pfaa/4btpYiS1xKS72BV1R4dsyqVX7SUlERDBjgwzSScdTQs8BJJ/k82KhRcM45mvOSjdCwIZx9tkf6vvWW/6p3771Vkb4PPqhI3wyihp4lzj/f95gMHeobCkU2illVJvvixZ43sWwZnHyyz8Kffz7MnRt1lVINNfQsMniwz3FdfbXHBIhskiZNqjLZX38dDjzQJ1Vbt1akb5pTQ88iZn662VFHwZ//7EdYimwys6oGXhnpu2CBZ1A0bw5/+5sifdNMXA3dzAaZ2Uwzm2FmI80s38yGm9l8M5sae+uQ7GKlerVqeSPfbz//bfmll6KuSLJCZaTvJ5/4zrYuXeCGG6BVK0X6ppFqG7qZFQMDgdIQQlugJtAn9ukLQwgdYm/aiJ4m8vM9kG/XXeGYY/zQaZGEqFGjqoF/9plP2FRG+paUeMRveXnUVeaseIdcagF1zawWUABUJK8kSYQGDfxGqqjIlzbOnBl1RZJ11oz0HTPG7yCuusqz2nv08L+AWkebUtU29BDCYuAmYCGwBPg2hDA29ulrzOxDM7vVzOqs6/lmdqaZlZlZ2VLtUU+pbbaBsWP9jv3gg/2GSiThatWqauCffAIXX+yRvt27+5DMtdfCF19EXWVOiGfIpRHQA2gJFAH1zKwv8DegNfB7oDFw8bqeH0IYFkIoDSGUFhYWJqxwiU/LlvDyy/DDD9Ctm6KyJclatvTlVosW+WRqy5bw97/73bwifZMuniGXbsD8EMLSEMIqYDSwZwhhSXArgfuBPZJZqGy6XXf18xCWLPFD5b/+OuqKJOvVrl3VwCsjfV991SN9W7euWucuCRVPQ18IdDazAjMzoCsw28y2BYh9rCcwI3llyubq0sUnSufOVZiXpNiakb4PPghbbQV/+YtvWDrpJEX6JlA8Y+jvAaOAycD02HOGAQ+b2fTYx7YC/pHEOiUBunb1g23KynzIc8WKqCuSnFK3blUme2Wk75gxHunbrp1nyijSd7NYSOFPxtLS0lCmqM7IPfyw/7s67DAYPRry8qKuSHLW8uXw6KOeWTFpEhQUwPHHQ//+HiAmAJjZpBBCtf9DtFM0B/3xj/Cvf8Fzz/lvvFpZJpHZYouqTPa1I31LSz2/XZG+cVNDz1H9+8P11/vNUf/+GsKUNLBmpO8dd3jA/xlnVEX6Tp8edYVpTw09h110kR8Sf889fuykmrqkhYYNqzLZ33oLevb0SN927fxUlxEjFOm7HmroOe7qq/00sltu8fdF0kZlpO+DD/oKmZtvhq++8nHCZs0U6bsOaug5zgxuu82DvC6/3N8XSTtNmngDr4z07dYNhgzxNe0HHODLtxTpq4Yunrd0zz0e5DVokB9pJ5KW1o70HTwY5s+H3r090vfvf/c/5yg1dAE8juPhhz3z5Ywz/N+LSFpr2tQz2T/+2LdCd+niM/077ACHHuo76XIs0lcNXf6rTh1fl77XXr608amnoq5IJA41a3oDr4z0vewyn1Dt2TPnIn3V0OU3Cgrg+ed9GXDv3v6+SMZo1swb+Gef/W+kb8+efuLL6tVRV5k0aujyP+rX9yTUdu18XP2VV6KuSGQjrRnp+/HHvkb37bf9Tj6LI33V0GWdttzSs9R/9zv/dzF+fNQViWyi7bf3Bl5e7pNDJSU+edq8uf8a+sYbWbMJQw1d1qtxY088bdnSc1/efjvqikQ2w5qRvrNnwznn+K+fBxzgyx9vvTXjs6XV0GWDCgu9qRcV+W+rH3wQdUUiCVCZyb54MTzwgEf6nn++/0U/6SS/e8nAu3Y1dKnWttv6TU2TJn5AxlQdBy7Zom7dqkz2adOgXz+fTN1rL2jf3lPsMijSVw1d4tKsmTf1Bg18k94MHWci2aYyk72iAoYN81zpAQP8rv2MMzzeN82poUvcSkrgtdd8vXrXrr4LWyTrbLFFVQP/4APPZ3/kEU+D/P3vPSjshx+irnKd1NBlo7Rq5U3dzOeSPv446opEkqgyk70y0vennzy/vajIJ1XTLNJXDV02WuvWPlH688/e1BcsiLoikSSrjPSdPh3efBOOPNIbfbt2foTeiBFpcaZjXA3dzAaZ2Uwzm2FmI80sf43PDTEzHSmSY9q29aa+fLlnJampS04wq8pkr4z0/fJLn1gtLvbDr+fNi6y8ahu6mRUDA4HSEEJboCbQJ/a5UmDLpFYoaatDB1/G+803auqSgyojfefO9XHIrl3h9tt9N17XrvDEEymP9I13yKUWUNfMagEFQIWZ1QRuBC5KVnGS/nbf3e/U1dQlZ1VOKD3+uEf6XnMNfPKJb2JKcaRvtQ09hLAYuAlYCCwBvg0hjAXOAZ4JISxJbomS7tTURWKaNvUG/sknHunbuXNVpG9ZWdJfPp4hl0ZAD6AlUATUM7OTgF7AkDief6aZlZlZ2dKlSze3XklTauoia6iM9H36af/HcP310LFj0l/WQjXbW82sF3BICKFf7M8nAVcCdYHKad0WwKchhFYb+lqlpaWhLAU/pSQ6kyb5xqMtt4Rx4zy1VEQ2j5lNCiGUVve4eMbQFwKdzazAzAzoCtwSQmgaQigJIZQAP1bXzCU3rHmnvt9+HkstIqkRzxj6e8AoYDIwPfacYUmuSzLY7rtXrX5RUxdJnbhWuYQQLg8htA4htA0hnBhCWLnW57dITnmSqUpL1dRFUk07RSVp1NRFUksNXZJq7aau1S8iyaOGLkm3ZlPfZx8Feokkixq6pERpqR/d+OOP3tQVvSuSeGrokjIdOvja9NWrYd99dUiGSKKpoUtKtW3rTb1WLR9TnzIl6opEsocauqRc69YwfjwUFHimkQ6eFkkMNXSJRKtWMGECNGrkSaNvvRV1RSKZTw1dIlNS4k29aVM4+GAfihGRTaeGLpFq1syHX7bbDrp39+WNIrJp1NAlcttu60sad9wRjjgCnn8+6opEMpMauqSFrbeG11+HXXaBo46Cp56KuiKRzKOGLmmjSRM/mrFjR+jVCx56KOqKRDKLGrqklS239HH0ffaBE0+Ef/0r6opEMocauqSd+vX9OMYjjoABA+C666KuSCQzqKFLWsrPhyefhOOPh7/9zd+qOS1RJOfViroAkfXJy4MRI/yO/brr4LvvYMgQqKHbEJF1UkOXtFazJgwdCg0bwo03elO//37PghGR39I/C0l7ZnD99d7UL70Uvv8eHn3Uh2VEpEpcv7ya2SAzm2lmM8xspJnlm9m9ZjbNzD40s1FmpnNFJWnM4JJL4Pbb4emn4fDDYfnyqKsSSS/VNnQzKwYGAqUhhLZATaAPMCiE0D6E0A5YCJyT1EpFgHPPheHDfWfpQQf5KUgi4uKdXqoF1DWzWkABUBFC+A7AzAyoC2gNgqTEySfDE09AWZlnqn/+edQViaSHaht6CGExcBN+F74E+DaEMBbAzO4HPgdaA0PW9XwzO9PMysysbOnSpQkrXHLb0UfDc8/BRx/BXnvBJ59EXZFI9OIZcmkE9ABaAkVAPTPrCxBCODX2sdlA73US4E9VAAAK+0lEQVQ9P4QwLIRQGkIoLSwsTFjhIgcd5Pkv33zjTV2nH0mui2fIpRswP4SwNISwChgN7Fn5yRDCr8BjwDHJKVFk/Tp1gjffhNq1/ZxSZapLLounoS8EOptZQWy8vCsw28xawX/H0I8AdI67RKJNGz/xqFkzPyhj9OioKxKJRjxj6O8Bo4DJwPTYc4YBD5jZ9NjHtgWuSmKdIhvUvDlMnFiV1DhsWNQViaReXBuLQgiXA5ev9eG9El+OyKZr0gRefdUb+llnwRdf+EYks6grE0kNpWJIVqlXzzcenXgiXHaZr1tfvTrqqkRSQ1v/Jevk5fnmo623hptvhq++ggcegDp1oq5MJLnU0CUr1agBN90E22wDF10Ey5b5ZGn9+lFXJpI8GnKRrHbhhZ7O+MYbfgpSRUXUFYkkjxq6ZL1TToFnn/VdpV26wMyZUVckkhxq6JITDj0UJkyAn3/2XaXagCTZSA1dckbHjvDOO1BU5BuQHn006opEEksNXXJKSYnvKu3c2c8rveEGnVUq2UMNXXJOo0Ywdiz07g0XXwznnAO//hp1VSKbT8sWJSfVqQOPPAItWvhZpYsX+58LCqKuTGTT6Q5dclaNGj7kcscd8MwzcMABoMh+yWRq6JLzBgzwTUfTpvnY+uzZUVcksmnU0EWAnj1h/Hj44Qdfq/7KK1FXJLLx1NBFYvbYA95/38fVDz0U7ror6opENo4ausgaWrTwZY2HHAJnnw3nnacVMJI51NBF1lK/vkfwnnce/POf0KMHfP991FWJVE8NXWQdataEW2/1YZeXXvK4gM8+i7oqkQ1TQxfZgP794cUXYeFCP5D6vfeirkhk/eJq6GY2yMxmmtkMMxtpZvlm9rCZzY197D4zy0t2sSJROPBAz4ApKIB994WHHoq6IpF1q7ahm1kxMBAoDSG0BWoCfYCHgdbArkBd4PQk1ikSqTZt/O68c2c/3u6CC+CXX6KuSuS34h1yqQXUNbNaQAFQEUJ4IcQA7wPNklWkSDooLPT16QMG+NF23bvD119HXZVIlWobeghhMXATsBBYAnwbQhhb+fnYUMuJwEvJKlIkXeTleVTA3Xd7pvoee+jADEkf8Qy5NAJ6AC2BIqCemfVd4yH/AiaEECau5/lnmlmZmZUtVVCGZInTT6/aWdq5M4wZE3VFIvENuXQD5ocQloYQVgGjgT0BzOxyoBA4f31PDiEMCyGUhhBKCwsLE1GzSFro0gXKynx8/aij4MorYfXqqKuSXBZPQ18IdDazAjMzoCsw28xOBw4Gjg8h6K+x5KTiYj/a7uST4Yor4Oij4Ztvoq5KclU8Y+jvAaOAycD02HOGAUOBbYB3zGyqmV2WzEJF0lV+Ptx/P9x2Gzz/PJSWwtSpUVcluchCCs/fKi0tDWVlZSl7PZFUe/ttOO44+OoruPNO6Ncv6ookG5jZpBBCaXWP005RkQTac0+YPBn23tsnTk87DX78MeqqJFeooYsk2NZbe/7LZZfB8OE+efrRR1FXJblADV0kCWrW9FUvL7wA5eU+rv7441FXJdlODV0kiQ45BKZM8aWNvXv7mPry5VFXJdlKDV0kyVq0gIkT4ZJLfDVMx46+fl0k0dTQRVIgLw/+8Q944w346ScfV7/hBm1EksRSQxdJoX33hQ8/9EOpL74YunWD+fOjrkqyhRq6SIo1auQTpPfc40Mvu+7qgV+bcrf+0Ue+mWnlysTXKZlHDV0kAmY+QTpjBvzhD3DuubDffhu/vPHJJ+Hww2HVqqSUKRlGDV0kQi1a+BF3998P06dDu3a+3DHezUgLFsBWW8EWWyS1TMkQaugiETODU07xXPUjj/SQrzZt4IknoLpkjvnzYbvtUlGlZAI1dJE0UVQEjz3mB2c0auSZMF26+Bj5uhr7ypV+1unuu6e8VElTaugiaWbffWHSJD8V6fPPfYx8t91gyBD/M3iDv+oq+P576NUr2nolfaihi6ShmjU93Oujj+C++3xYZuBA2HZbz2DfemsYPNgPrO7aNepqJV3UiroAEVm/vDw49VR/mz4dXn7Zx9pr1IADDoA+fbzZi4AaukjG2HVXfxNZHw25iIhkCTV0EZEsEVdDN7NBZjbTzGaY2Ugzyzezc8zsYzMLZrZVsgsVEZENq7ahm1kxMBAoDSG0BWoCfYC3gG7AZ0mtUERE4hLvpGgtoK6ZrQIKgIoQwhQA0xS7iEhaqPYOPYSwGLgJWAgsAb4NIYxNdmEiIrJx4hlyaQT0AFoCRUA9M+sb7wuY2ZlmVmZmZUuXLt30SkVEZIPimRTtBswPISwNIawCRgN7xvsCIYRhIYTSEEJpYWHhptYpIiLViGcMfSHQ2cwKgJ+ArsAmnYg4adKkr8xsUydRtwK+2sTnZipdc27QNWe/zb3euDI1LVSXzwmY2ZVAb+AXYApwOnAWcBHQFPgSeCGEcPqmVhtHDWUhhNJkff10pGvODbrm7Jeq641rlUsI4XLg8rU+fHvsTURE0oB2ioqIZIlMaujDoi4gArrm3KBrzn4pud64xtBFRCT9ZdIduoiIbEDaNXQzO8TM5saCv/66js/XMbPHYp9/z8xKUl9lYsVxzeeb2Swz+9DMXjOzjD8WuLprXuNxx8YC4DJ6RUQ812tmx8W+zzPN7JFU15hocfy9bmFmb5jZlNjf7e5R1JlIZnafmX1pZjPW83kzs9tj/08+NLOOCS0ghJA2b3jw1yfA9kBtYBqw81qPORsYGnu/D/BY1HWn4Jr3Bwpi7/8pF6459rj6wATgXTwcLvLak/g93hFfEtwo9ueto647Bdc8DPhT7P2dgQVR152A694H6AjMWM/nuwMvAgZ0Bt5L5Oun2x36HsDHIYRPQwg/A4/isQNr6gE8EHt/FNDVMjshrNprDiG8EUL4MfbHd4FmKa4x0eL5PgNcDdwArEhlcUkQz/WeAdwZQvgPQAjhyxTXmGjxXHMAGsTebwhUpLC+pAghTAC+3sBDegAPBvcusKWZbZuo10+3hl4MLFrjz+Wxj63zMSGEX4BvgSYpqS454rnmNfXDf8Jnsmqv2cx2A5qHEJ5LZWFJEs/3eCdgJzN7y8zeNbNDUlZdcsRzzVcAfc2sHHgBODc1pUVqY/+9b5R0O1N0XXfaay/DiecxmSTu64mFopUC+ya1ouTb4DWbWQ3gVuCUVBWUZPF8j2vhwy774b+BTTSztiGEb5JcW7LEc83HA8NDCDebWRdgROyaVye/vMgktX+l2x16OdB8jT83439/DfvvY8ysFv6r2oZ+xUl38VwzZtYNuAQ4MoSwMkW1JUt111wfaAuMM7MF+FjjMxk8MRrv3+unQwirQgjzgbl4g89U8VxzP+BxgBDCO0A+nnmSzeL6976p0q2hfwDsaGYtzaw2Pun5zFqPeQY4Ofb+scDrITbbkKGqvebY8MO/8Wae6WOrUM01hxC+DSFsFUIoCSGU4PMGR4YQNikULg3E8/d6DD75TexIx52AT1NaZWLFc80L8bA/zKwN3tCzPWP7GeCk2GqXzvj5EksS9tWjnhVezyzwPHyG/JLYx67C/0GDf9OfAD4G3ge2j7rmFFzzq8AXwNTY2zNR15zsa17rsePI4FUucX6PDbgFmAVMB/pEXXMKrnln/CjLabG/1wdFXXMCrnkkfhDQKvxuvB/QH+i/xvf5ztj/k+mJ/nutnaIiIlki3YZcRERkE6mhi4hkCTV0EZEsoYYuIpIl1NBFRLKEGrqISJZQQxcRyRJq6CIiWeL/ARDtmbQO13YYAAAAAElFTkSuQmCC\n", "text/plain": [ "