diff --git a/20-InteriorPoint/Seminar20en.ipynb b/20-InteriorPoint/Seminar20en.ipynb new file mode 100644 index 0000000..279bc91 --- /dev/null +++ b/20-InteriorPoint/Seminar20en.ipynb @@ -0,0 +1,1043 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Seminar 20\n", + "# Interior point method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Reminder\n", + "\n", + "- Projected gradient descent\n", + "- Frank-Wolfe method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Convex optimization problem with equality constraints\n", + "\n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "&\\min f(x) \\\\ \n", + "\\text{s.t. } & Ax = b,\n", + "\\end{split}\n", + "\\end{equation*}\n", + "where $f$ is convex and twice differentiable, $A \\in \\mathbb{R}^{p \\times n}$ and $\\mathrm{rank} \\; A = p < n$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Dual problem\n", + "Relation between dual and conjugate functions \n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "g(\\mu) & = -b^{\\top}\\mu + \\inf_x(f(x) + \\mu^{\\top}Ax) \\\\\n", + "& = -b^{\\top}\\mu - \\sup_x((-A^{\\top}\\mu)^{\\top}x -f(x)) \\\\\n", + "& = -b^{\\top}\\mu - f^*(-A^{\\top}\\mu)\n", + "\\end{split}\n", + "\\end{equation*}\n", + "\n", + "Dual problem\n", + "\n", + "$$\n", + "\\max_\\mu -b^{\\top}\\mu - f^*(-A^{\\top}\\mu)\n", + "$$\n", + "\n", + "**Approach 1**: find conjugate function and \n", + "\n", + "solve unconstrained optimization problem" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "**Issues**\n", + "- it may be not so easy to find solution of the primal problem from the dual one\n", + "- conjugate function $f^*$ has to be twice differentiable for fast solving of dual problem, but this is not always hold" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Optimality conditions\n", + "\n", + "- $Ax^* = b$\n", + "- $f'(x^*) + A^{\\top}\\mu^* = 0$\n", + "\n", + "or\n", + "\n", + "$$ \\begin{bmatrix} f' & A^{\\top} \\\\ A & 0 \\end{bmatrix} \\begin{bmatrix} x^{\\\\*} \\\\ \\mu^{\\\\*} \\end{bmatrix} = \\begin{bmatrix} 0 \\\\ b \\end{bmatrix} $$\n", + "\n", + "**Approach 2**: solve generally non-linear system with Newton method.\n", + "\n", + "**Q**: in what case the system becomes linear?" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Newton method for convex optimization problem with equality constraints\n", + "\n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "& \\min_v f(x) + f'(x)^{\\top}v + \\frac{1}{2}v^{\\top}f''(x)v\\\\\n", + "\\text{s.t. } & A(x + v) = b\n", + "\\end{split}\n", + "\\end{equation*}\n", + "\n", + "From the optimality condition follows \n", + "\n", + "$$ \\begin{bmatrix} f''(x) & A^{\\top} \\\\ A & 0 \\end{bmatrix} \\begin{bmatrix} v \\\\ w \\end{bmatrix} = \\begin{bmatrix} -f'(x) \\\\ 0 \\end{bmatrix} $$\n", + "\n", + "**Newton direction $v$ is defined only for non-singular matrix!**\n", + "\n", + "**Q:** how direction $w$ can be interpreted?" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "**Exercise**. \n", + "\n", + "Estimate number of iterations required for convergence of \n", + "\n", + "Newton method for quadratic objective and linear equality constraints." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Linearization of optimality conditions\n", + "\n", + "- $A(x + v) = b \\rightarrow Av = 0$\n", + "- $f'(x + v) + A^{\\top}w \\approx f'(x) + f''(x)v + A^{\\top}w = 0$\n", + "\n", + "or\n", + "\n", + "- $f''(x)v + A^{\\top}w = -f'(x)$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Pseudocode\n", + "**Important note:** initial point has to lie inside the feasible set!\n", + "\n", + "```python\n", + "def NewtonEqualityFeasible(f, gradf, hessf, A, b, \n", + " \n", + " stop_crit, line_search, x0, \n", + " \n", + " tol):\n", + " \n", + " x = x0\n", + " \n", + " n = x.shape[0]\n", + " \n", + " while True:\n", + " \n", + " newton_matrix = [[hessf(x), A.T], [A, 0]]\n", + " \n", + " rhs = [-gradf(x), 0]\n", + " \n", + " w = solve_lin_sys(newton_matrix, rhs)\n", + " \n", + " h = w[:n]\n", + " \n", + " if stop_crit(x, h, gradf(x), **kwargs) < tol:\n", + " \n", + " break\n", + " \n", + " alpha = line_search(x, h, f, gradf(x), **kwargs)\n", + " \n", + " x = x + alpha * h\n", + " \n", + " return x\n", + "\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Stopping criterion\n", + "\n", + "Estimate the following difference\n", + "\n", + "$$\n", + "f(x) - \\inf_v(\\hat{f}(x + v) \\; | \\; A(x+v) = b),\n", + "$$\n", + "\n", + "where $\\hat{f}$ is quadratic approximation of function $f$.\n", + "\n", + "To do this multiply both side by $h^{\\top}$ from the left \n", + "\n", + "$$\n", + "\\langle h^{\\top} \\rvert \\cdot \\quad f''(x)h + A^{\\top}w = -f'(x)\n", + "$$\n", + "\n", + "and use constraint $Ah = 0$\n", + "\n", + "$$\n", + "h^{\\top}f''(x)h = -f'(x)^{\\top}h\n", + "$$\n", + "\n", + "Then \n", + "\n", + "$$\n", + "\\inf_v(\\hat{f}(x + v) \\; | \\; A(x+v) = b) = f(x) - \\frac{1}{2}h^{\\top}f''(x)h\n", + "$$\n", + "\n", + "**Summary:** value of $h^{\\top}f''(x)h$ is the most natural stopping criterion of Newton method." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Convergence theorem\n", + "\n", + "Convergence of the equality constrained Newton method is equivalent \n", + "\n", + "to convergence classical Newton method for unconstrained optimization problem.\n", + "\n", + "**Theorem**\n", + "Assume the following conditions hold\n", + "- level set $S = \\{ x \\; | \\; x \\in D(f), \\; f(x) \\leq f(x_0), \\; Ax = b \\}$ is closed and $x_0 \\in D(f), \\; Ax_0 = b$\n", + "- for any $x \\in S$ and $\\tilde{x} \\in S$ hessian $f''(x)$ is Lipschitz\n", + "- in the set $S$ $\\|f''(x)\\|_2 \\leq M $ and norm of the inverse matrix of the KKT system is bounded above\n", + "\n", + "Then, Newton method converges to the pair $(x^*, \\mu^*)$ \n", + "\n", + "- linearly when iterands far from the solution\n", + "- quadratically in sufficiently small neighbourhood of the solution" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Infeasible starting point\n", + "\n", + "- Newton method requires that starting point is feasible\n", + "- But what to do if this requirement is violated? Example of hard case is the problem in which domain of the objective function is not $\\mathbb{R}^n$.\n", + "- Assume starting point is infeasible, then KKT conditions can be written as\n", + "\n", + "$$\n", + "\\begin{bmatrix}\n", + "f''(x) & A^{\\top}\\\\\n", + "A & 0\n", + "\\end{bmatrix}\n", + "\\begin{bmatrix}\n", + "v\\\\\n", + "w\n", + "\\end{bmatrix}\n", + " = -\n", + "\\begin{bmatrix}\n", + "f'(x)\\\\\n", + "{\\color{red}{Ax - b}}\n", + "\\end{bmatrix}\n", + "$$\n", + "\n", + "- If $x$ is feasible, then the system is equivalent to the system for Newton method with feasible starting point" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Primal-dual interpretation\n", + "\n", + "- Method is called *primal-dual*, if every iteration updates both primal and dual variables\n", + "- In particular, re-write optimality condition in the following form\n", + "\n", + "$$\n", + "r(x^*, \\mu^*) = (r_p(x^*, \\mu^*), r_d(x^*, \\mu^*)) = 0,\n", + "$$\n", + "\n", + "where $r_p(x, \\mu) = Ax - b$ and $r_d(x, \\mu) = f'(x) + A^{\\top}\\mu$\n", + "- Solve system with Newton method:\n", + "\n", + "$$\n", + "r(y + z) \\approx r(y) + Dr(y)z = 0\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "- Primal-dual direction in Newton method is defined as the solution of the following system\n", + "\n", + "$$\n", + "Dr(y)z = -r(y)\n", + "$$\n", + "\n", + "or more detailed\n", + "\n", + "$$\n", + "\\begin{bmatrix}\n", + "f''(x) & A^{\\top}\\\\\n", + "A & 0\n", + "\\end{bmatrix}\n", + "\\begin{bmatrix}\n", + "z_p\\\\\n", + "z_d\n", + "\\end{bmatrix}\n", + " = -\n", + "\\begin{bmatrix}\n", + "r_p(x, \\mu)\\\\\n", + "r_d(x, \\mu)\n", + "\\end{bmatrix}\n", + "= - \n", + "\\begin{bmatrix}\n", + "f'(x) + A^{\\top}\\mu\\\\\n", + "Ax - b\n", + "\\end{bmatrix}\n", + "$$\n", + "\n", + "- Replace $z_d^+ = \\mu + z_d$ and obtain\n", + "\n", + "$$\n", + "\\begin{bmatrix}\n", + "f''(x) & A^{\\top}\\\\\n", + "A & 0\n", + "\\end{bmatrix}\n", + "\\begin{bmatrix}\n", + "z_p\\\\\n", + "z_d^+\n", + "\\end{bmatrix}\n", + "= - \n", + "\\begin{bmatrix}\n", + "f'(x)\\\\\n", + "Ax - b\n", + "\\end{bmatrix}\n", + "$$\n", + "\n", + "- This system is equivalent to the previous one in the following notation\n", + "\n", + "$$\n", + "v = z_p \\qquad w = z_d^+ = \\mu + z_d \n", + "$$\n", + "\n", + "- Newton method gives direction for update primal variable and updated value for dual variable" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Pseudocode\n", + "```python\n", + "def NewtonEqualityInfeasible(f, grad, hessf, A, b, \n", + " \n", + " stop_crit, line_search, x0, \n", + " \n", + " mu0, tol):\n", + " \n", + " x = x0\n", + " \n", + " mu = mu0\n", + " \n", + " n = x.shape[0]\n", + " \n", + " while True:\n", + " \n", + " z_p, z_d = ComputeNewtonStep(hessf(x), A, b)\n", + " \n", + " if stop_crit(x, z_p, z_d, grad(x), **kwargs) < tol:\n", + " \n", + " break\n", + " \n", + " alpha = line_search(x, z_p, z_d, \n", + " \n", + " f, grad(x), **kwargs)\n", + " \n", + " x = x + alpha * z_p\n", + " \n", + " mu = mu + alpha * z_d\n", + " \n", + " return x\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Stopping criterion and backtracking\n", + "\n", + "- Update $r_p$ after step $z_p$\n", + "\n", + "$$\n", + "A(x + \\alpha z_p) - b = [A(x + z_p) = b] = Ax + \\alpha(b - Ax) - b = (1 - \\alpha)(Ax - b)\n", + "$$\n", + "\n", + "- Total update after $k$ steps\n", + "\n", + "$$\n", + "r^{(k)} = \\prod_{i=0}^{k-1}(1 - \\alpha^{(i)})r^{(0)}\n", + "$$\n", + "\n", + "- Stopping criterion: $Ax = b$ and $\\|r(x, \\mu)\\|_2 \\leq \\varepsilon$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "- Backtracking: $c \\in (0, 1/2)$, $\\beta = (0, 1)$\n", + "```python\n", + "def linesearch(r, x, mu, z_p, z_d, c, beta):\n", + " alpha = 1\n", + " while norm(r(x + alpha * z_p, mu + alpha * z_d)) >= \n", + " (1 - c * alpha) * norm(r(x, mu)): \n", + " alpha *= beta \n", + " return alpha\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Convergence theorem\n", + "\n", + "The theorem result is similar to the case of feasible starting point\n", + "\n", + "**Theorem.** Assume that\n", + "- sublevel set $S = \\{(x, \\mu) \\; | \\; x \\in D(f), \\; \\| r(x, \\mu) \\|_2 \\leq \\| r(x_0, \\mu_0) \\|_2 \\}$ is closed\n", + "- in the set $S$ norm of the inverse KKT matrix is bounded\n", + "- hessian is Lipschitz in $S$.\n", + "\n", + "Then convergence is \n", + "- linear far from the solution and\n", + "- quadratic in the sufficiently small neighbourhood." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## General convex optimization problem\n", + "\n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "& \\min_{x \\in \\mathbb{R}^n} f_0(x)\\\\\n", + "\\text{s.t. } & f_i (x) \\leq 0 \\qquad i=1,\\ldots,m\\\\\n", + "& Ax = b,\n", + "\\end{split}\n", + "\\end{equation*}\n", + "where $f_i$ are convex and twice smoothly differentiable, $A \\in \\mathbb{R}^{p \\times n}$ and $\\mathrm{rank} \\; A = p < n$. \n", + "\n", + "Assume that the problem is strictly feasible, i.e. Slater conditions are satisfied." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Optimality conditions\n", + "\n", + "- Primal feasibility\n", + "\n", + "$$\n", + "Ax^* = b, \\; f_i(x^*) \\leq 0, \\; i = 1,\\ldots,m\n", + "$$\n", + "\n", + "- Dual feasibility\n", + "\n", + "$$\n", + "\\lambda^* \\geq 0\n", + "$$\n", + "\n", + "- Lagrangian stationarity\n", + "\n", + "$$\n", + "f'_0(x^*) + \\sum_{i=1}^m \\lambda^*_if'_i(x^*) + A^{\\top}\\mu^* = 0\n", + "$$\n", + "\n", + "- Complementary slackness condition\n", + "\n", + "$$\n", + "\\lambda^*_i f_i(x^*) = 0, \\qquad i = 1,\\ldots, m\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Idea\n", + "\n", + "- Reduce the problem with **inequality** constraints to sequence of **equality** constrained problems\n", + "- Use methods for solving equality constrained problems" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "\\begin{equation*}\n", + "\\begin{split}\n", + "& \\min f_0(x) + \\sum_{i=1}^m I_-(f_i(x))\\\\\n", + "\\text{s.t. } & Ax = b,\n", + "\\end{split}\n", + "\\end{equation*}\n", + "where $I_-$ is an indicator function\n", + "\n", + "$$\n", + "I_-(u) = \n", + "\\begin{cases}\n", + "0, & u \\leq 0\\\\\n", + "\\infty, & u > 0\n", + "\\end{cases}\n", + "$$\n", + "\n", + "**Issue.** Now objective function **is not differentiable**." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Logarithmic barier\n", + "\n", + "**Idea.** Approximate function $I_-(u)$ with function\n", + "\n", + "$$\n", + "\\hat{I}_-(u) = -t\\log(-u),\n", + "$$\n", + "\n", + "where $t > 0$ is fixed parameter.\n", + "\n", + "- Functions $I_-(u)$ and $\\hat{I}_-(u)$ are convex and non-decreasing\n", + "- But $\\hat{I}_-(u)$ is **differentiable** and approximates $I_-(u)$ while $t \\to 0$" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGLCAYAAAD05of+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xt03Hd95//nV3PTXaORLFmyJSvjxHES5yY7IRcgBGSu4bfQWixld9vSHqLdhXZ7fmyjDbvbFtjfuvbZc5YWWlYOtLTdsk0lynIpBaxCgUIAOwqO49iJ7fFdlixpLrrN/fv5/TEjWbJlSZZG0oz8epwjvqPPfC8f43Pwi8/n/f18LGMMIiIiIpJbRWvdAREREZH1SCFLREREZAUoZImIiIisAIUsERERkRWgkCUiIiKyAhSyRERERFaAQpaIiIjIClDIEhEREVkBClkiIiIiK8C51h2ora01LS0ta90NERERkQW9+OKLw8aYDYs5d81DVktLC4cPH17rboiIiIgsyLKsc4s9V9OFIiIiIitAIUtERERkBShkiYiIiKwAhSwRERGRFaCQJSIiIrICFLJEREREVoBCloiIiMgKUMgSERERWQFrvhjpzYjH4wSDQcbGxkin02vdHVkFDoeDiooKfD4fHo9nrbsjIiKyaAUTsuLxOOfPn6e6upqWlhZcLheWZa11t2QFGWNIJpOMjo5y/vx5mpubFbRERKRgFMx0YTAYpLq6mtraWtxutwLWLcCyLNxuN7W1tVRXVxMMBte6SyIiIotWMCFrbGyMysrKte6GrJHKykrGxsbWuhsiIpKn7LRNOmWvdTdmKZiQlU6ncblca90NWSMul0t1eCIickP/+BfH+fInf7bW3ZilYEIWoCnCW5j+7kVEpNAUVMgSERERKRQKWSIiIiIrQCFLREREZAUUzDpZkh/6+vo4fPgwfr+fQCCA3++nra1t0df39PTg9Xpv6hoREZFCpJBVgPr6+mhvb+f06dOr+txAIEBnZycHDx6cbmtvb8fv9+P3+xe8vre3l4985CN0d3evZDdFROQWZMxa9+B6mi4sQM8//zxer3fVn9vV1UVHR8esto6ODjo7O+e9LhAI0NHRQSAQwOfzrWQXRUTkFpZv76ErZBWg3t7eNZlu6+npobW1dVbbrl276Onpmfc6v99PV1cXTz/99Ep2T0REJK8oZBWgvr4+du/evarPDIfDc45ETY2oBQKBVe2PiIhIvlsXNVmf/MYxXu0fXetuzOvuxkp+/733LPn63t5euru7p8NMd3c33d3ddHR0XDe6tBKm9g280TTlVBG8iIiIZKyLkHUraGtro62tjf379xMMBunq6lrUde3t7YTD4UU/p6ura86wdDP3EBERkXUSspYzQlRoDh48eFP1WHqTT0REZG2oJqvA9Pb2rno91kwa0RIRkbyVZ68XKmQVkL6+PiDzRt9qm5pCnKrNmjIVulSPJSIiMtu6mC68VfT29tLa2jpdfB4OhxdcLytXNVlerxe/33/dvYLB4PR3IiIictWiQpZlWXsAP7A1e+wyxvRcc04rsAsIZM8JGGN6c9vdW9uhQ4dmjWIdOHCAZ555Zt5rclmT1dbWxuHDh2e9zdjX16ctckREROaw4HRhNmAFjDH7jTEdQDuwz7Ksp2ec4wf2GWMOGGN6jTEHgI5su+TQzp07gbVZkHTfvn3Xhbauri727ds3/Xs4HGbr1q0cOHBgznsEg0HVdYmIyC1hMSNZ/pmjVsaYsGVZ+4AuYOpf0o7s7zN1AfvIhDLJgWeffZa9e/fi8/nwer2rsj7WTF6vl3379tHZ2clDDz00vZfhtVOFc9Vt7d27l0AgQDgcnt7/cPfu3ezZs2c1/wgiIiKrZt6QZVmWF/iXlmUdMMbMHH7ozX7vN8YEgD1cH7IOAweRnGltbV3zJRlaW1vnDXder5dQKHRd28zRLhERkZzLwx2i550uzAYrf/ZnTtkg5gdmDV9MhTJNGYqIiMhqsKz8WsNhwZosY0y1MabvmuY2IJwdxfJlz7tRoY1CloiIiNxylrpOVgewN/t5/jUE5mBZ1tOWZR22LOvw0NDQErsgIiIikr9uOmRl3yoMGmP2L/Wh2bcQdxljdm3YsGGptxERERHJWzcVsrL1VR3GmOv2dcnWZomIiIgINz+StQ942zVtgezRN7NxRugKICIiIrKC8u/dwpsIWZZldQGd1xa4Z38PcH1tlo+rxfEiIiIit5RFhaxsHda+mYHJsqy2Gcsz9JLZUmem1my7iIiIyC1nsdvqAHgty2rN/rQB7TNCVyfXr+zekW0XERERueUsZsX3Gy0xPj2qld1qpzO73c4hMmtj7dNUoYiIiNyq5g1Z2XqrRS2fml2w9NpFS0VERERuSUtdjFREREQkf+Th64UKWSIiIrIu5NnWhfNPF4pcq6+vj8OHD+P3+wkEAvj9ftra2ua9JhwOc+DAAfbs2YPP5yMYDNLV1cXu3bsXvFZERKRQKWQVoL6+Ptrb2zl9+vSqPjcQCNDZ2cnBgwen29rb2/H7/fj9N94HPBgM0tnZSWdn5mVTr9fLc889p4AlIiLrmqYLC9Dzzz+P17v6uxh1dXXR0dExq62jo2M6PM3n4MGDhEIhTp8+TSgUYs+ePQteIyIiUsgUsgpQb2/vmowC9fT00NraOqtt165d9PT0LOp6r9c774iXiIjIeqKQVYD6+vrYvfu6PbpXVDgcJhAI4PPN2qJyekQtENCSaCIiIjOtj5qsf/hPMHB0rXsxv433wrv+cMmX9/b20t3dPR1muru76e7upqOj47rRpZUQDAYBbjhNOVUEfyOBQGB6xCsYDOLz+TRlKCIiOWPycAmH9RGybgFtbW20tbWxf//+6bfzFqO9vZ1wOLzwiVldXV1zhqWbuce1pka/Zoaq9vb269pERETWk/URspYxQlRoDh48eFP1WN3dN9oVafV4vV6efvrpWW0dHR10dHQoZImIyLqlmqwC09vbu+r1WDMtZ0Rrpql1tnJ1PxERkXyjkFVA+voyW0Pu2rVr1Z89NYU4VZs1ZSokzVePtX///uvapqYQVTAvIiLr1fqYLrxF9Pb20traOl18Hg6HF1wvK1c1WVPLL1x7r2AwOO/SDFMLmO7Zs2fWOVNhTUs6iIjIeqWQVUAOHTo0axTrwIEDPPPMM/Nek8uarLa2Ng4fPjzrbca+vr55a8T8fv+cwe3awCgiIrI8+fd6oaYLC8zOnTuBtVmQdN++fdeFtq6uLvbt2zf9ezgcZuvWrRw4cGC6zefzzZoWDIfDdHV18dxzz618p0VE5NaRZztEaySrgDz77LPs3bsXn8+H1+tdlfWxZvJ6vezbt4/Ozk4eeuih6anAa0eprq3b2rNnDz09PfT09DAyMkI4HKa7u1tThSIisq4pZBWQ1tbWNV+SobW1dd5w5/V6CYVC17VrqQYREbnVaLpQREREZAUoZImIiIisAIUsERERKXz593KhQpaIiIisD3n2cqFCloiIiBS+PBzIUsgSERGRdSAPU5ZCloiIiBQ8Y4ymC0VERERWRn6lLIUsERERWR/yK2MpZImIiEjhM6rJEhEREVkZqskSERERybU8HMpSyJIl6enpobe3d627ISIiAmRWcLDybChLIasA9fX1sXXr1jV7fm9vLx/5yEfW7PkiIiLXyb+BLIWsQvT888/j9XpX/bmBQICOjg4CgQA+n2/Vny8iIjKfPBvIUsgqRL29vbS1ta36c/1+P11dXTz99NOr/mwREZH55GFJlkJWIerr62P37t1r3Q0REZE8kn8pSyGrQPT29tLR0TEdrrq7u+no6KCvr2+NeyYiIpIHMpXva92LWZxr3YFc2PfzfZwInljrbsxru287nQ93Lvn6trY22tra2L9/P8FgkK6urkVd197eTjgcXvRzurq68Pv9S+2miIjImsjDjLU+Qtat5ODBgzdVj9Xd3b2CvREREZEbWRchazkjRIWmt7eXzs5b588rIiKyKHlY+b4uQtatYqr+ateuXWvcExERkfySHh3DjibXuhuzKGQVkN7eXlpbW6fXyAqHwwuul6WaLBERuRUkBwZITsTWuhuzKGQVkEOHDs0axTpw4ADPPPPMvNeoJktERG4FhjyrekdLOBScnTt3Amu3IOmUYDB4UyNkIiIiKy7P6rI0klVAnn32Wfbu3YvP58Pr9dLa2rqqzw+Hw+zdu5dAIEA4HKazs5ODBw+ye/du9uzZs6p9ERERyXcKWQWktbV1Taf/vF4v+/btW7Pni4iIzCffJgw1XSgiIiIFL1OTlV/ThQpZIiIiUvjyK18BClkiIiKybuRX0lLIEhEREVkBClkiIiKyLhRs4btlWXssy7puYSbLsryWZT1jWZY/+9lvWda+uc4VERERWQlmxn/mi0WFrGxgeu4GX/uAfcBpIAS8CBwyxvTmpIciIiIiC8q3cawF1smyLMsPdJIJTsF5Tt0NHAZ8xphA7ronIiIisjj5FrPmDVnZwNQBYFlW5wLnhgHtsyIiIiJrpACnC0VERETyWX7Fq4xcbavjtyxravM6HxA0xvTk6N4iIiIiC8uzpJWLkBUEmBmqLMvqtiyLGwUty7KeBp4GaG5uzkEXRERE5NZmYeVZylr2dKExJmyMOXBNcxeZNw5vdM0BY8wuY8yuDRs2LLcLIiIicovLr3iVsVI1WQEyU4jeFbq/rLGenh56exe3Skc4HGb//v0EAgHC4TCBQIDOzs5FXy8iIrI4+RW1lh2yLMt6Zo7mqeUe/Mu9v1yvr6+PrVu3rtnze3t7+chHPrLo84PBIJ2dnWzdupXq6mp27tzJQw89RFub1qsVEZFcsfJuCYdlhazsOlr7sseZfNmj1sxaAc8//zxe7+oPEgYCATo6OggEAvh8voUvmOHgwYOEQiFOnz5NKBRiz549C18kIiKySAbArKORrKl1tOZYgLQN6MuunSU51tvbuyajQH6/n66uLp5++uklXe/1evH7NbgpIiK5Z0xhF777gLmGT4IzR7KydVgdwOLnk+Sm9PX1sXv37rXuhoiISF7Jt5C10LY6XuBZMrVVXjJTg7uBg1PLMxhjerKbR+8BarLntWt7ndzq7e2lu7ubQCDzX2t3dzfd3d10dHTQ2tq6xr1bWCAQoKcns6JHMBjE5/NpylBERHLGYJFvhe8LbasTJrN34bzWeuHRgf/+34kfP7GWXViQ567tbPzEJ5Z8fVtbG21tbezfv59gMEhXV9eirmtvbyccXvysbVdXV86n9Kbqt2aGqvb29uvaRERElsrk4TpZuVrxXVbJwYMHb6oeq7u7ewV7szher/e6Oq6Ojg46OjoUskREJCcMBTZdWCiWM0JUaHp7e+nsXHBwMe/5/f7pdbPW4k1JERFZX/JxJEsbRBeQvr4+AHbt2rXGPbk5+/fvv65tagpxqsZMRERkOYwpsJosyS+9vb20trZOj/wsZhRorWuyplZ337Nnz6z7BoOZ9Wq1pIOIiOSKlWfrZClkFZBDhw7NGsU6cOAAzzwz14L7V611TdbU2lrXhqlrA6OIiMhyGCwsK79ClqYLC8zOnTuBtVuQdEowGJxzhCwcDrN161YOHLi6Z7jP55s1LRgOh+nq6uK5555blb6KiMj6l481WRrJKiDPPvsse/fuxefz4fV6V319rHA4zN69e6cL1js7Ozl48CC7d++e9Zbg1FTglD179tDT00NPTw8jIyOEw2G6u7s1VSgiIjljsPJuWx2FrALS2tq6ptN/Xq+Xffv2LXhOKBS6rl1LNYiIyErLt5EsTReKiIhIwcvH6UKFLBERESl4mZCVXxSyREREpODl496FClkiIiKyLmi6UERERCTHNF0oIiIisgJU+C4iIiKyAlSTJSIiIpJjxhjQtjoiIiIiuTW90HuerfiukCUiIiIFzU7bABRhr3FPZlPIEhERkYJmpzMjWEWqyRIRERHJnamQZWkkS0RERCR3ro5k5VfIcq51B6Sw9PT0EAgEOH36NIFAgI6ODvbs2bPW3RIRkVvYVE2WZRSyZJn6+vpob2/n9OnTq/rcnp4e/H7/dKgKh8Ps3LmTYDDI008/vap9ERERmZJOZUeyrPwKWZouLEDPP/88Xq931Z8bCARobW2d/t3r9dLZ2UlHR8eq90VERGTK9EiWCt9luXp7e2lra1vVZ4bDYZ5//nnC4fCs9ql+BAKBVe2PiIjIlOmarDybLlTIKkB9fX3s3r17VZ/p9XoJBAIKUyIiknem3y7Ms+nCdVGT9aO/fZ3hC+Nr3Y151TaV86YPbFvy9b29vXR3d0+HnO7ubrq7u+no6Jg1hbeSQqHQnP3yer34/f5V6YOIiMi18nUka12ErFtBW1sbbW1t7N+/n2AwSFdX16Kua29vv26Kbz5dXV03FZi6urp49tlnF32+iIhIruXriu/rImQtZ4So0Bw8ePCm6rG6u7tXrC8HDhzA5/PxzDPPrNgzREREFpKeGsnSBtGyHL29vatejzWXQCBAV1cXBw8eXOuuiIjILS5f18lSyCogfX19AOzatWuNewKdnZ384z/+41p3Q0RE5GpNVp6NZK2L6cJbRW9vL62trdNrZIXD4QXXy1qJmqyOjg727du3Jmt1iYiIXOvq3oUKWbJEhw4dmjWKdeDAgQXroXJdk3XgwAE6OztnBbHe3l78fr/eMBQRkTWRTqnwXXJg586dwNosSNrT0wNkRtCmpi6DwSDd3d2LfttRREQk19LJqZCVXuOezKaQVUCeffZZ9u7di8/nw+v1rtr6WJAJVu3t7XN+pxEsERFZS6lsyHKQWuOezKaQVUBaW1tXdEmG+Xi9XozJr7luERERuDqS5cizwne9XSgiIiIFLZXMTBM6HApZIiIiIjmTSmRHsooUskRERERyJpW0KTIprCLHWndlFoUsERERKWjpRBqHSUGRtdZdmUUhS0RERAqaRrJEREREVkAqaWdHsvIr1uRXbxagJQRuXfq7FxGRG0lPjWQ58ivW5Fdv5uFwOEgmk2vdDVkjyWQShyO/hoFFRCQ/pJLZmiwrv2JNfvVmHhUVFYyOjq51N2SNjI6OUlFRsdbdEBGRPJRK2NjpGMeCx9e6K7MUTMjy+XyEQiGGh4dJJBKaProFGGNIJBIMDw8TCoXw+Xxr3SUREclDyXgaKx1lJBFc667MUjDb6ng8HpqbmwkGg5w9e5Z0Or82gZSV4XA4qKiooLm5GY/Hs9bdERGRPJSMp3Gk45j8WsGhcEIWZIJWQ0MDDQ0Na90VERERyROJaIqydIy4I79SVsFMF4qIiIjMJRFL4UjHCncky7KsPUDYGNM7x3etwC4gAPiBwFzniYiIiOSSbRtSCRtHOo5t5VfKWlTIsiyrDXgOaJ/jOz+wzxize0Zbt2VZAWNMIGc9FREREblGMpYCwJmK5t1I1rzThZZl+S3L6iIzOnWjkv0OoOuati5g3/K7JyIiInJjiVjmRThHOo7JsyKoebtjjAkYYzqMMQfmOW0P0HdN2+Fsu4iIiMiKSUyPZMUweTZduKzMZ1mWlzlGuYwx4ez3/uXcX0RERGQ+yamRrEKbLlwEH1wNVXNQyBIREZEVk4hmRrJcyRjpdbaEgzcnvRARERFZgvhkJmS5E+PrLmQtiWVZT1uWddiyrMNDQ0Nr0QURERFZB6LjSQDcsTHs9RiysrVZi2aMOWCM2WWM2bVhw4ZcdEFERERuQbGJbMhKTq67kaypdbBm7dw7I3RpnSwRERFZMbHxJJ5SJ0XGJu1Y697MtqyQlS14D3B9bZaPzOrwClkiIiKyYmITSTylmXS13kayAHrJbKkzU2u2XURERGTFxMYTFJdkQlaqgEOWj7nfJuzk+u12OrLtIiIiIismOp6cDln5Vvg+796F2dqqZ8msd+UF9lmWtRs4aIzpgcyUoWVZnZZl7QMOZc/dp6lCERERWWmxiSQ+nwfIv+nCeUNWtuZqwREpY0wf12+tIyIiIrJijDHExpJ4PMVA/oWsPNtKUURERGRx4pMpUkmb0tJMuFLIEhEREcmBiXAcgNKSTLjKt5oshSwREREpSBORTMgqKTaARrJEREREcmJ6JMudBiDlVMgSERERWbaJcAKAYisGQMKdX7Emv3ojIiIiskgT4TjFZS6sZCZkpd35ta+OQpaIiIgUpPFwnDKvGxPLhKykRrJERERElm9sJEpFTQl2NApAUiNZIiIiIstjjCEyHKOqtmR6JCvtnneN9VWnkCUiIiIFJzqWJBVPU7mhBDuaCVkpTReKiIiILE9kKDNFWFlbjB2LYltgO/Mr1uRXb0REREQWYXRoEoCqDSWYaIyU24Fl5Vesya/eiIiIiCxCZDgGFlTWlGDHoqTcRRQpZImIiIgsT3hgggpfMQ5XEfb4BAmPAwut+C4iIiKyLMHLE9Q0lgGQHhslXuLEshSyRERERJYsnbYJDUziaywHwB4dI1aikSwRERGRZYkMRrHTBt81I1mqyRIRERFZhpH+cYDpkGWPjhEtza/V3kEhS0RERApMsH8Cy4LqjaUApMfGiBc7VJMlIiIishxXzo1R3VCG0+XAJBKYaJRoiUPThSIiIiJLZYzhytlR6lsqgcwoFkBche8iIiIiSzc6HCM2kaRuKmSFQgBESzVdKCIiIrJkV86OAkyPZKWGRwAYK3dqJEtERERkqQbPjuJwFeHblHmzMDU8DMB4hWqyRERERJas/2SY+pZKHI5MhEmPZELWWIVGskRERESWJDaRZOjCGJvurJ5uSw2PgMvFuMfgLHKuYe+up5AlIiIiBaH/ZBgMbL7TO92WGh7G6fORNCmFLBEREZGluPRaCKeriPqWqum21PAwzpoa0iatkCUiIiKyFBdfC7FxaxUO19X4krzcj7OxgZSdwmHl19Y6ClkiIiKS90aHowT7J9iyo2a6zRhD8lI/7k2bSdmaLhQRERG5aWeOZN4ivO3+2um2dDCIiUZxbdpE2k7jKnKtVffmpJAlIiIiee/My8NUN5RRtaF0ui156RIArs2bSRlNF4qIiIjclPhkkssnw7NGsWBGyNq0iaSd1HShiIiIyM04/dIQtm3w379hVnvi/AWA6elChSwRERGRm/D6zweoqiuhrqViVnv89CmcGzfiKC/LvF1YpOlCERERkUUZC8a49HqYO9+wEcuavW1O4tRpPLffDkDKpHBaGskSERERWZSThwbBwLaH62e1G9smHgjguf12knYS29h4HJ416uXcFLJEREQkLxnb8Oo/99OwtWrWW4WQKXo3sRie27cST8UBKHYWr0U3b0ghS0RERPLSheNBIkNRdrxl03XfxV97DQDP7bcTS8cAKHGWrGr/FqKQJSIiInnp6D9dpKTSzdYH6677LvryUXA68WzfTjQVBTSSJSIiIrKgyFCUs6+McM8bG3E4r48r0aMvU7xtG0XFxcRSmZGsYodCloiIiMi8Xjp4niKHxY43Xz9VaGyb2NFXKL7/PgCNZImIiIgsxngozvGf9HPXY42Uea9/YzBx5gz2+Dgl92ZC1tRIlmqyRERERObx0sFzGBta39485/eTP/85AKU7WwGmC981XSgiIiJyA+OhOK/+qJ87H66nsnbukamJn7yAq7ERV3MmhE0kJwAodZXOef5aUcgSERGRvPGzbwSwjeGhp26b83uTTjPxs59R+tij0yvAj8ZHAah0V65aPxdDIUtERETywtCFMU68cJn7nmy64ShW7Ngx7NFRyh55dLotkogAUOlRyBIRERGZxRjDT75yCk+pk53v3HLD88a+9z0oKqLs8cem20bjoxQ7irWtjoiIiMi1Th4a5OKJEA8/5ae4zHXD88a+e5DShx7CWV093TaaGM27qUJQyBIREZE1FhtP8s/dJ6lrqWTHE9evizUlfuoUiUCAirfvntUeiUeodFdAdg/DfJGTkGVZlteyrGcsy/JnP/sty9pnWVZbLu4vIiIi69eP/+4U8YkUT/7r7RQVWTc8b/Qfvg1ARdvskBWOh6kauwJ/8oYV7efNytVIlg/YB5wGQsCLwCFjTG+O7i8iIiLr0JkjQ5z4yWUe2N1M7ebyG55nbJvIV79K2WOP4qqfvZfhlckr1BW5V7qrNy2X04W7gWpgqzGm2hjTk8N7i4iIyDozEYnzvb86QW1TOQ/fYMmG6XNfeIFkfz9Vv/zLs9qNMVyZvEK95QLrxqNga8GZy5sZY8JAOJf3FBERkfXH2IZ//NKrpOJp3v6b9+BwzT/uE/nKV3BUVVHRNrsSKRKPkLAT1DldQH6FLBW+i4iIyKo79K2zXDge4o0fuIPqjWXznpu8dInR73yXqve/nyLP7GUaBicHAahjfY9k+S3L2pP97AOCmjIUERGRawV+McShb55h+yMbufuNjQueH/zLvwTLwver/+a676ZCVr3lJt9GsnIVsoIAM0OVZVndlmUxV9CyLOtp4GmA5ua5N38UERGR9SfYP0Hvn79K3ZYKnvhXd05vjXMj6XCYUHcPVe95N67G6wPZ2chZALYU5ddCpJCj6UJjTNgYc+Ca5i4ybxzOdf4BY8wuY8yuDRs25KILIiIikucmwnG++bkjOD0O3vVv78Xpcix4zcgX/wwzOYnvN35zzu/PjZ6jylNFNY68my5cyZqsAJkpRO8KPkNEREQKQHwyyTc+e4TYRJKnPnof5dXFC16THBwk+Jd/SeV730vxndvmPOfs6Fm2VG4BDPk2XZirxUifmaM5mD36c/EMERERKUypRJpvff4ooYEJ3tVxL3VbFrcFztBnPwu2zYb/8B9ueM6ZyBlaKlvAmPU3kmVZlh/Ylz3O5MseA8t9hoiIiBSmTMB6mf5TYdp+/W6a7vYtfBEQPXKEyFf+juoPfQj35rm32hmODjMUHWK7b3u2ZZ2FLGNMAOjIHmdqA/qya2eJiIjILSaZSPP3f/oyF06EeOu/uYs7Hqpf1HUmmeTy7/0+zro6an/rYzc879jwMQDuqbmHzHRhfsnZ24WWZfmngla2DqsD+EiO7i8iIiIFJBFL8a3Pv8yl18O87dfuYvsjDYu+duRLXyL+2mts/txncZTfeKudYyPHKLKKMiNZeThdmJOQZYzpsSxrT3adrBrAC7TPMbolIiKf7R5fAAAgAElEQVQi69xEJM7f/8nLDF8cp+3X7+bON2xc9LXRV44x9MefpeLtb79udfdr/eLKL9jq3UqpqzTbsg5DFjDnelgiIiJyawkNTPCNzx4hOpbg3f/uXlrurV30tfbEBP0f/zjOmhoaPvXJec+Np+P0XemjfVt7pmG9jmSJiIiIXHotxD8cOEpRkcX7P9666LcIIbPR8+Xf/wMSFy7Q/KU/x+GdfwWol668RDwd59HGR2e0KmSJiIjIOmKM4eXvXeTHXzmFt66E93z0fqo2lNzUPUYOPMfoN7/Jht/5HcoefnjB83986cc4LSc763dmO2HnW8ZSyBIREZGlSybSfP+vTnDy0CD+Bzbwtl+7C3fJzcWLsd5ehv7n/6Tyqaeo6Xh6wfONMXz37Hd5pPERylzZzaXtJBS5lvJHWDEKWSIiIrIkI/3jHPziMUb6J3jDv/Cz8x1bsIpubjhp4mc/59LH/yPF991Hw3/79IJ7GQK8PPwy/RP9/PsH/v3VRjsFDoUsERERKWDGGF75wSV+/JVTuIsdvPdj99N8T81N3yd65AgX/92/w9W0maau/0VR8cJb7QB88/Q3cRW5eGvzW682plNQlF+xJr96IyIiInltcjTB9//qOGePjtB8Tw1v+7W7KK103/R9okePcv7pDhy1tTR/8c9wVlcv6rrxxDhfP/113tHyDircFVe/sFPg9Nx0P1aSQpaIiIgsyBjD6z8b4EfdJ0nFbd74gTu478nNi5reu9bECy9w4aMfw+nz0fxnf4arvm7R13799NeZTE3yoe0fmv2FnYKispvuy0pSyBIREZF5jY5E+cFfv8b5V4Ns9Ffy5L++C1/j0gLN6He/S//H/yPulhaavvCFmwpYyXSSvzj2F9xXex/3brh39pd2UjVZIiIiUhjSaZuj37/Iz75xBoA3fuAO7n3LZopusrgdMiNhI11dDH3mjyh54AGa/tfnF1wL61pfOfkV+if6+a+P/tfrv7TTqskSERGR/Hf+1RH++W9PEhqYpPmeGp740DYqa25u7asp9uQk/f/5PzP2D9+m8r3vpeHTn1p0kfuU8cQ4B14+QGtdK483Pn79Cam4RrJEREQkf0WGJvnn7lOcfXmYyg0lvPvf30fLvTVLqr0CiJ86xaX/9+PET56k7j9+HN9v/uaS7vW5X3yO4egwf/TkH819fXIS3KrJEhERkTwzOZrg8D+c5dgPL+FwFvHo+7dy/1ubcLiKlnQ/Ywzh5/+Wwb17KSoro+lAF+VvetOS7nV06ChfPv5lPrj9g9fXYk1JjIO7fEn3XykKWSIiIrew+GSSlw6e58j3LpJO2tz1eAMPP3UbZVVLXw4hNTTEwKc+xdjBXsoee4zGfX+Ic8OGJd1rLDHGMz98hvqyen7rwd+a+yRjIDGhkSwRERFZe4lYild+cIm+75wjPpnijl11PPxeP9760iXf0xhD5O++yuD+/ZholLrf/V18H/51rKKlj4Z98oVPcnniMl9655dmr4s1UzqRWcJBIUtERETWSmw8yZHvX+Do9y8Sn0zRfE8Nj/wLPxuabxBgFilx9iyXP/lJJl/4KSW7dtLwqU/j8d+2rHt+/sjn+c7Z7/A7rb/DA3UP3PjE+FjmeKMQtkYUskRERG4BE+E4L/We59iP+knF0/gf2EDrO7dQ31K5rPumR0cZ/tPPE/zrv6bI7WbjH/w+3g98YMmjV1O+evKrfP7I53nf7e/jN3b8xvwnj1/JHMuXNiW5UhSyRERE1rGh82O8/L0LvH54EGPDHQ/V0fqOLdQ0Lq9I3KRShHt6GPqjPyYdDuPd88ts+O3fXnLt1UxfO/U1/uCFP+Cxxsf4vUd/b+G3EccHM8fy+mU/O5cUskRERNYZO21z5sgwR753gcunIjg9Du5+vJEH2pqp2rC0ta6mmHSa0W99i+HP/QmJc+co3bWL+k88S/Hdd+ek7z2v9/CpFz7FGxrewGee/AyuokWsfTU9kqWQJSIiIitgcjTBiRcuc/QHFxkPxqmoKebxPbdz12MNeEqXt1CnsW3GvvMdhj73JyROn8Zz551s/txnKX/b25a8htZMtrH5474/5ouvfJE3bnojn3nyM3gci3zDMXIhc1TIEhERkVwxtuHC8SCv/rifM0eGsdOGTdu8vOkD22i5r3ZJW+DMZCcSjH7jm4z8+Z+ROHUa9+1b2fSZz1Dx9t3LrruaMp4Y5xP//Am+f+H77Nm2h088/AlcN7N6+/BJqNwEHq2TJSIiIss0Foxx4oXLHP/xZcaCMYrLXNz75GbufqxxyZs3z5SORAj9zfME//dfkR4axrN9O43/439Q+a53YjkcOfgTZBwZOkLnDzu5PHGZ//Twf+JD2z908yNjw69D7R0561OuKGSJiIgUiEQsReClIV7/+QAXT4QwBpruqubRX9qK//4NS16dfabYa68Tfv5vCP/fr2EmJyl7/HF8f/iHlD32WE6mBack00m+8MoX6DrSRX1pPV9655d4sO7Bm79RKgFXjsOuD+esb7mikCUiIpLH0imb868Gef3nA5w9MkwqaVNZW8zOd7Vw12MNVNYur5AdwI7FGPvOdwj9zfNEX3oJy+2m8l3vwvcbH6b4zjtz8KeY7cXBF/n0C5/mdOQ0777t3fyXR/7LjRcaXcjlI5CKQvMjue1kDihkiYiI5BnbNlw+FebU4SucfHGQ+ESK4jIX2x9rYNvDG9nor1z2qJIxhvjx40S+9jXC//dr2JEI7i1bqOvspOp9/wJndXWO/jRXDUwM8NmXPsvXT3+dxrJG/uRtf8KbN795eTc9+8PMsfnR5XcwxxSyRERE8kA6ZXPptRCnXxrizJEhomNJnK4ibntgA9serqfpbh8Ox/KnA5OXLxP5xjcZ/cbXiZ88BS4XFW97G9Uf/JeUvuENOZ0SnBKJR/jiK1/ky8e/jG1sPrzjw/zb+/4tpa6lb+Ez7dWvw6adUF63/HvlmEKWiIjIGkkl0px/NUjgpSHOHh0mPpnC5XHQcm8N/gfraL7Hh7t4+f9Up0Ihxr/3PSJf/waTP/85GEPJAw+w8fd/j4p3vnNFRq0AQrEQ/+fE/+F/H//fjCfGecr/FB998KNsKt+UmwcMvQ6XfwG7P52b++WYQpaIiMgqGg/FOffKMOdeGeHCiRCpeBpPqZPb7qvF31pH013VOF3Lf3svNTTEWG8vo9/9LpM/PwTpNK7mZmo/+lGq/p/34m5uzsGfZm794/38xbG/4Kunvko0FeUtTW/hYw98jDt9Oa7v+umfgsMD938wt/fNEYUsERGRFWTbhitnRzn3yghnjw4zfGEcgPJqD3e+YSNbH9hA453enEwFJi5eYvwfexn97kGifX1gDO6WFmp+8zep2L2b4h33rMh0IGQWE/1p/0/pfr2b71/4PhYW7/G/hw/v+DBbvVtz/8DQWfjFlzMBKw+nCkEhS0REJOcmRxNcfC3I+VeCnDs2Qmw8iWXBxq1VPPI+Py331uJrLFt+8XoiwWTfS4z/8IeM//AHJE6dBsCzbRu1H/0oFW/fjeeOO1YsWAEMR4f5+umv0/1aNxfHL1LtqeZX7/5VPnTXh9hYtnFlHmoMfPtZKHLCE50r84wcUMgSERFZpmQ8Tf+pMBeOB7l4PMTIpcxolafMyZZ7athybw3Nd9dQXLa8rW0AkoNXmPjRDxn/wQ+Z+MlPsCcmwOWi7KFdVLe3U/7EE7hbWpb9nPmMJ8b53oXv8a3At/jp5Z+SNml21e/itx78Ldq2tOF2uFf0+bz45/DatzK1WFU5qu9aAQpZIiIiN8lO21w5P8bF40EuHA8xEIhgpw1FTouGrV4eeZ+fprt81DZVLHtbm3Q4zMShQ0y+8FMmfvYzEqczo1XOjRupfM97KH/izZQ98ghFZctf5X0+E8kJftL/E7595tv84OIPiKfjbCrfxG/s+A2e8j+F3+tf0edPO/ND+IdOuL0NHv3Y6jxziRSyREREFpBO2Vw5N0b/yRD9JyMMnA6TiKUBqG0q5/63NtF0l4+Nt1fhci+vaN2emGCyr4+Jn/6UyRd+Suz4cTAGq6SE0l278L7/fZS96U14tm1b0WlAgCuTV/inC//E9y98n59d/hlJO0m1p5r33/5+3uN/D/dvuH/F+zDLuZ/Alz8IPj/80nOQo70TV4pCloiIyDWSiTSDgQj9J8P0nwozGBgllbQBqG4o446H6tl0ZzWb76ympGJ5U2OpkREm+/qIvtjH5Et9xI69CqkUuFyU3n8/tR/7KGWPPELJvfdiuVd2Gi5pJzk6dJQXLr/Ajy/9mKPDRwHYXL6ZD27/IE82PcmDdQ/iLFqD+PCLL8M3/gN4t8Cvfg1Kfavfh5ukkCUiIre8ydEEA4EIg2cywerKuTHstMGyoLapgnvetInGO7w03F61rFBljCFx5gzRvj4mX+wj2tdH4tw5ACy3m+J776Xmwx+m9JE3UNraSlHJ8rfMWag/ZyJneOHyC7zQ/wKHBg4xmZqkyCpiR80OfvvB3+bJpifZ6t26uiNWM0XDmSL3I1+G256A9i8VRMAChSwREbnFpFM2wxfGGTgTYfDMKAOBCGMjMQCKHBZ1Wyp4oK2Jxjuq2bi1Ck/J0v+pTAWDxI4eJXr0lczx5ZdJh0IAOLxeSlpb8X6gnZIHWynecQ9FKzxSlbbTvB56nb4rfbx05SX6BvsYig4B0FzRzFP+p3i08VEe2vgQVZ6qFe3LgmwbjnZD7x/A+CC8+XczbxI6lv/ywGpRyBIRkXXLGMNEOM5AYJTBMxEGAqMMnR8jncpM/ZV5PWz0V3LvWzaz0V/FhqZynEusqUqPTxA7dozYK0eJvnyU2NGjJPv7M19aFu6tfsrf8hZKd7ZS0tqK+7bbVnx0aCI5wbHhY9Oh6sjQESaSEwA0lDWwa+MuHt74MI82Ppq7VdiXy7bh5Hfh+/8fDLwMDffDB/8aNrWudc9umkKWiIisC5lAlWDo/ChXzo8xdH6MK+fGiI4mAHA4i9jQXMGOt2xi421VbPRXUl5dvKRnpUIh4idOEDt+gvhrJ4i+coxEIJBZvwlwbd5M8f33Uf2v/hXF9+6g+O57cJSv7Nt/0VSU14KvcWzkGMeGj3Fs5BhnImcwGCws7qi+g6f8T9Fa18qDdQ/SUN6wov25aYkJePlvM6u4D78O3mb4pS/Ajl/O+wL3G1HIEhGRgrNQoLKsTIF6890+NjRXsPG2KmqbynE4b+4fa2PbJM+fJ3biNWInjhM/foLYiROkBgenz3HW1VF8991UvvtdlNx7L8U7duD0rWzN0HhinJPhk7wefJ1Xg6/yyvArnA6fJm0ybzzWFNewo3YH77ztneyo2cH9dfdT6a5c0T4tiW3DuX+GI38Dr34NEuOZkatf+gLc876Cmhqci0KWiIjkNds2jA5FGb44zvDFMYYvjjN0bozJGYHKu/FqoKprrqC2qQKX5+am/VLBIPGTp4ifOkn81CniJ14j9tprmMnJzAkOBx6/n9KHH6Z4+3aK79qOZ/v2FQ1UKTvF+bHzvB56ndeDr3MydJKT4ZNcGr80fU6Vp4odNTt4YvMT3FN7DztqdlBXWrd2heoLScbg7I8yi4m+9m0Y6wd3BdzzfnjgQ9D8aOYvdR1QyBIRkbwRj6YYuTjO8MVxRi5ljsH+cVKJTA2VVWThrS+l6S4fG7YsLVDNDFOJ06ezn09NF6QDFFVUUHznnXh/6Zemw5Tn9tsp8nhy/mcGSKaTXBi7wJnIGc6MnuFM5AwnQycJRALE03EAHJaDlsoW7qu9j1++45fZVr2NO6rvoKGsIX8DFWRGq64cg7M/ziwkGvgnSE6Aqwxufyvc/Wm4893gLl3rnuacQpaIiKy6maNTU2Fq5OI4Y8HY9DmeMie1m8u5+42N1G4up2ZTOb7GMpyuhQOVsW1SAwPEz5whcfbs1TB1+jTpYHD6vKLycjy3305F29twb92K5/Y78NxxO866lRkJCsfCnB09mwlTUz+jZ7g4dnF6qg+gvrQef5WfX9n+K9xRfQfbqrfhr/Kv/HY1uZCYgIGjcOnFzOKh534M0WyA9W7JbOh857uh5Y3gWlpN3FxeuRRhNJrksdtrc3bP5VLIEhGRFWOnbSJDUYKXJwhdniDYP0Hw8iThwcnpN/wsC7z1pWz0V3LPmxup2VRO7eYKyrzuBYNOKhTKhKiz50hMBaqzZ0mcO4eJx6fPmw5Tb3vrioYp29gMTQ5xYewCF8cvcmHsQubzWOZzOB6ePtdd5GZL1Ra2VW/jHS3v4Laq27it6jZaKlsoc61skXzOJCbgygm4/BJcegn6X4Kh42Ayf7dU3wbb3wMtb4Itj4O3acW68oUfBXjpQpgf/O6TK/aMm6WQJSIiy5ZO20QGs2FqYCpMTRAenMROm+nzKmqK8TVm6qeqG8qo2VSGr6Fs3mUT0mNjJC9cIHH+wtUQdfYsiTNnSEciV090OnE3NeFuaaHs8cdxt7RM/zjrNuQsTE0mJxmYHODS2KXrQtTF8YvT03sARVYRDWUNNFU00baljZbKlukw1VjWiKNoeVvwrJpkNPPG35UTmRB1JfsTPnf1nNIaaGzNhKpNrdD4IFRsXL0u2gbHMveJzDWFLBERWbTYRJLwYGYkKjQ4SWRwkuBA5mjb2TBlQWVNMb6GMlruraG6IROkqjeWzVk7ZdJpkpcukbhwkcSF8yQvXCR5MROqkhcuzA5SgLO+HndLCxXvfCfu2zIhytPSgmvTJizX8t5GM8YQjAW5PHGZyxOX6R/vZ2BigP7x/um2maNRACXOEjZXbGZL5RYe3/Q4TRVN0z8N5Q24igrkDbnEJITOQjAAoTOZ49RP5OLV0akiJ9TckQlSD/5r2LA980agt3lNC9bTaYMrz5Z6UMgSEZFZ0kmb8NAkkcEoocEJwleihAcmCV+ZJDaenD7PKrKo2lCCt76U2+6rxddQiq+xHO/G0lmbJBtjsEdHSQZeJ3rxIsnzF0hcvEByKlT1X4bk1fvidOLa1Ih7cxPFO+7B3dSMq2lzZpSquZmisqVNpRljiMQjDE4OMhQd4srkFQYnB6dD1MDEAJcnLs8aiYJMiGosa6ShvIEdtTtoLG9kY9lGNpVvoqmiiZrimvwuPJ8Si0DkUiYwjV7MHCOXIHIhE6TGLs8+v6Q6sxHz5ofh/l/JhKm6u8C3FZz5VxuWTNs4Hfn196CQJSJyC7LTNuOhOJGh6PTIVPhK5jg2EptaUxOA0ko33vpS/A9swFtfire+lOr6Uipqi3E4ijDpNKmhIZL9/SRfvkzk2/0k+y+RvHyZVH8/yUv92FPLIGQVVVXhbmrKrC/19nfgam7C3dSEa3MTro31WM6b++cplopxZfLK9M9QdIjBycHM58nM56HJIRJ24rpra0tqaShrYFv1Np7Y/AQN5Q00lDXQWN5IQ1kDle7K/A5R6SRMDGW2nhm/AmMDmc+jl2YEqYuQGJt9neWAykao2gz+JzOByndb5qf6toLZH3BKcDJBdWl+hT+FLBGRdSqZSDM6FCUyFGV0OHvM/j42Ers6vQc4PQ68dSXUt1Sy7Q0bqc6GKW9dKU6SJAcGMiGq/yip45eZvNRPpL8/0zY4CKnUrGc7qqpwbmrEtWULpY88iquxEVdDQ2ZEavNmHFUL74uXttOE4iFGoiOMxEYyx+zn4egwI9GR6RGp0cToddcXO4qpK62jrrSO+zfcP/15Q+kG6kvrM59LNuTnG3vpJEwGYXIEJocz4Wl88JoglW2bHAHM9fcorc0EqJqtcNubM5+rNkFVE1RuytRLFUpN2AJs2xAYmuBdO1avBmwxFLJERAqUMYbYRHJWeBodihLJBqrJyOxRG0+pk8raEjY0V7B1Zx2VtcVUlNiUMY57/AqpK+dIDQyS+vkgyYFBxgcHCQ0OYl9TE0VREc66OlyNjZQ8+CCVDQ24NjVmglQ2TN1oSi9tpwnHwwwHX7suOI1Es+Ep+zkUD2FP1QHN4C5yU1tSS01JDU0VTeys3zkdoOpK66grqaOurI4KV0V+jEDZNsTC2cA01082TE0MX/09Hpn7Xg43lNdDeR1Ut0DTw5nfK+qz7dnvyuvBuTJreuWjbx8bIBJN8ngeLd8AClkiInktHk0xNhJldDjG2Ej2JxhjdCQTqBKx9Kzzy6rcVG4oofmuaspLDWWOKKWpMCUTV3CEBkgODJI6MUhycIDU4BUmYjEmZt7AsnDU1OCqr8e1eTOlu3birKvHubE+G6I24aqvw3K5MMYQTUUJxUMMxcKE4iFCseOEzv6EcDxMMBYkHA8TioUIxUOEY2Eiicicwcnj8FBTXENNSQ2N5Y3cW3svNSU11BTXTAeqqe/LXeWrH55ScYiGM3VNsewxGs5+Dt/gu6nfR5lzpAnAWZwZcSr1Zd7Oq27JHEtrrraV1lwNTyXV62Y19FxIpGx6XrzIf/v7V7mroTLvRrIsY27wF3+zN7KsVmAXEAD8QMAY07vQdbt27TKHDx/OSR9ERApNfDLJ6MjVADU6Ep0OUmMjMeKTs6fhnO4iyiudlJfalDkTlDJOaSJI8fggntBFGBogNTREKhiE9OwAhsuFq64OZ309ro31V8NTfT2O+nrSNVVMVLkZtSeJxCOzA9JUaIqFZwWoa4vEp/tpOfEWe/F6vFQXV1Ptqaa6uBqvxzsrME2FqDJX2coEp1QC4mOZeqT41M84xEez7eMz2mf8JMZnB6VUbP7nOEugxAvFXiiuyn6uuvr7dGCaOtZmjutwlfOVZIzhylicX1wI86OTQ3z32CBXxuI8fJuPz/7Kg9RX5m5x0xuxLOtFY8yuxZybk5Esy7L8wD5jzO4Zbd2WZQWMMYFcPENEpNAY2zA5lmA8GGcsGGM8NHMkKvM5Eb0mRDmhrNimzBFjs2OcYncQz9gAnvAlXINncAQvc10UKSqiqMaH2bABR00trm1+qC4nXl3ORHUxo1VOQhVFBItTRJKjjCZGicQjROKvMJr4CZHRCJHhCCk7de2dp1W4KvAWZwJTXWkd26q34Sv2ZdpmBKiptiVN1aVTme1WEhOZ5QSS2WNiYsbncUhOXv1+OjSNzR2mbhACZ7PAXQ6eCvBkj+5yqNs+OyhNh6hrg1TVLTU1t5KMMYzHUwyOxhmIxOiPRBmIxLgciRIYmuC1wTHCk5k3UcvcDh67vZYPvaGZJ+7YQFGerZEFuZsu7AC6rmnrAvYB7Tl6hohI3piqhxoPxhkPxRgPZY5j07/HmAgnZi3ECeCw0pQVRSlJj9IQywQod/AixeODFMdGcCUnroYojxvj85LyVZDwlRG5rYmJqq2MVjiJlMFQaZrBkgQD7kkiyTEiiQHGEidnd3Qy+5N9O7/MVUaVu4oqTxWV7kq2erdS5ama1VblufrZV+zD6/HicrgyxdjJaOYnFc1s9JuczIzyxKMwfhGSr2dDUTYczfl5jrCUmID09W/+zctVOiMcZX+qmq6GpJntUz/ucvBUzg5UrjLIs/WVCo0xhkTaZjKeZiKRIppIM5FIMxlPZY6JFJOJNBPxzHE0miQ0mSQSTRCeTBKaTBCJJglPJknZ18+w1Za7afaV8q4dDdxZX87djVU80OTF7czvv7dchaw9XB+yDgMHc3R/EZFVFY+mGA/ODE8xxq6MMz48wXg4wcRYmrQ9+/85WyZNcXocTyxI6cQw1fEQxbEQnniI4njm6LCjJCuLiVa4mSh3MlpqEd4Gw2VprpQ6GPC4GC6zCZVB1JMGKwgEZz3HaTmpcFdQ5a6k0l2Bz1VFS9lGqpylVDk8VBV5qLRcVBU5qcRBFQ6qjKHCNrjSiRlBKQbRKCRDmc/JyUxwSmW/n/6c/c5cM/24GA5PZkrMVZY5ussyn8vrroYkd+kcn8uy55bOPk59dpXecsEobRuSaZuUbUinDUnbJpXOtKVtQ8q2SaZNps22r56fvvrdVFsiZRNPXXtMX/972iaetKeP8bRNPDm7PZZMM5lIk54jHN1IictBdamLqlI33hIXd26soKrEjbfUhbfERV2lh4aqEhqrSqir9FC8iP0q89GyQ5ZlWV4yNViz/lfAGBO2LAvLsvyaMhSRfJFO2kxE4owPTzB2OcTY4CgTw5OMh6OMjyWJRiGWdJE21/zPo7HxJCJ4YiGK42GqpoJTLAQmQsIVZtw9RrgkzVApjNZApNQiUgaj2WOkFCaKLUqKbCqK0lRYFuU4qLAcVBgHXjw0mTLKbUOlbVOeTlORTmV+UgkqUgnKk3FKknGs1LmlhR4ALHCVZH6cJZlNemd+LvFmCrLnO8dVeoNzyjDuUmxnKbarFNtyYAzYxmBnj8ae+j3TZmZ+R+Z1/KvXzDgnAXbcYJs0xoyRtg1pY7BtM+MzpGwb2xjSdiaYZD5fPV79zOzrs/dI29k2c/WYzvZ59vWZYyo981yu6c+Me07fO/PclG1IZUPQVGCaGaSmAlLaznyfoxLqG3IWWbidRXicRdmj45rfi6hyu/BUeKZ/92TPK/M4KHU7KXU7KHM7KXHPbit1O2f97nLcGgE5FyNZPsiEqht87ydTDC8iczAm8w+Kmfk7ZNtMto3p48w2M/N8O/PdVLs9fd/sccbnqe+m7mNPP3PGcfr+899r5rOYvv81f44Zn2/4rGvvZWywk5BKYNkJLDuJlU5mPqeTFNlx0qko8dQEicQoqdER0qFJUqMp7HELO+qEuAeTKsG2y7EpJ11URdpZft3fgWUn8cQjeBIRyuIRKlIRjAmTskLEi8KMe0KMlUS4WGqI+QyJEkOq2GBKbCg2FFs2FbadCUa2wWfbbLEzbRW2TbltUzlpKB/PfJ7aZCVpuUlZblKWi1SRO/u7i5SV+Zz5KSdpuRnGxWXLRcLpJuV0kbDcJMgckziJ4yaOhxhu4paHOG5ilocYLmIUEzcuoriJGg8J48BgTYcYEwc7NiPozBVwptpmBaCp86PYZnK6bb1xFFk4LIuiIuKeGvoAABbpSURBVLJHa0abhbPIosjKthVZFFlkjzPb/v/27qXHceSwA/i/ig9RpNSvmX3YCBykN7nlNHaAADkkgGc/QICJDd9ymskn8MKfwJjcHWAWyAeIx0cfAswCvgTIwbuLILcg2XEMxAnWj5lZT0t8s3KoIkVRVOvdkqj/DyBIFYuPEru6/01SlJhajyMlLCnhSAHbErAtPW1JCccyZWa6KpMStiVuL5N6XbYlzDaEWZ+ELQUcS5fVA5JrQtShfe9fF2wjZF2suoAQ4jGAxwDwjW98Ywu7cLt//JdfIEyyxh+y6T8SZWHbH7dmGep/3CZFU+ucrHKyTl1v+o9XVapml5/a7lRZc52TMjS3M6ctmCmrWjdVhqmyyTqbbZnZbm15VVsnZsqm/9jfut3Ge4qWspn2qsn8udutvSdt72l9HUtvt6rX/j7v+j/SJokCDjLYyOEgg4McNnLYIoNryvU8Pd8Wea1+bVq0rAM5XJHBbqnvigwOMlhIAZGhkCmUzJDLFJnMkYscqcyRyQKpLBArBZUpIAWQAEUsIGIbMj6HTM9hpeew8jM94BwC5yjsC2T2u0h6fwwlGpcUVAGneAsr+wp2/jvYeAmFNyisr1DYb1A4XyF3v0LeG2PkSbx1JeDasGDBLmxYhQOrsNEvhvCLe3in6CGFDjdp5CCJHMTCRaocJMLB/8JBAheJcHRAgi5P4SIVjl7OhKkUDoQUEACkEBACEACE+SMshajK9XRZbzItBaZei8br5vIDITCsz5dl/el1LN6eKZMr1m+uX07qC7TUkcvtgyXETKiZCka1QFTVbQSn+nLTgYjBg9a3l+dkKaU+BvAxoB/hsOvt/cPP/gu/Gy13Q2X9F52YKtMzmmViqsws0yyrrdPMNo85EdXjTsQt66y2V9VtXyfmbEfUtlPfz2p9LdtFrU6zfdX7Ud+mLMtEta8Lt9tShqm2Ndrbtt2Z4zRdhtb3dHqdM9uFgkQOW+WwkMMyY6lySJXpsnJQGeS8eiqHNOV6OV2mlzHzVQpLZbWy5nQKqTLIYlImy6EoX6eNsrSqr8tSiHnP6KlRABIBREIiFAKhFIiEQCgkIil0mSkfC4mx1PPGUuC1kBhJiVBKjKVEUgggFpCxgBULuJHCIAKCCBiECoMQCGIf/ewMXnYGR51B4gwQZ0jcMyTu0IzPkLpDvYMOkDtAeYFMFiFsdQNbjtFzfo3A+z+4A8C7cOHfDzD8+j2c/8H7GF78KSynD2k5EJY9+RluBpvGz+xM8BGTvkVEtIythSwhxMUtlwz36l9/8G0A03+EAUz9sqU5lAKKHCgyQJlxkc8pyyZDXpsu0kmdPJ2uV9VPa+tYpn5L2VbWlQItD0rcOWEBlqOf5ixtPS0dU6aHzLIQWQ5CaeuxsBBKywSgehACIgFEAEIUCAUQqgIRCkQq19MqQ6gyhCpHVKQIixRRkaKoTkkq9FIdivwYCGLAjxSCqJwGBpHCWWzh/VhiGAkEEeDFNqzUh5IDxO5ZFZTKoSobniE5t9H810eKAn23QL8vcBVY8M9cBFc+gvsDDN4Zwr/nwz9z4Z+5sI/0RlgiOh3bCFnl/VZXAKqQZW6Ir8/fm1s/4qmU/iM7Exby2TJVLFEnn57ftlxrnbKsFlaaZXNDzjJ1Mv3VDgvX3aizj8BRJ6QOHdIxYxNGyumq3AYsezItHf3cGhksWb82WI6pO2+7k6EQEhGAWBSITJCJlTLjDLHSwSZWGaIiQ6xyRCpDVKSICzPOE0RFjDALEWURoizS03mEMA0R5ro8LdLp90ZBn9Zpu/fZhKSrtIfL1MVl6uIstfC12MIwsTCIHQSxAz9S8KIcXpijN87gjFM44xjWKIYoChTCQuoM9JklZ4jEHSJ1hkjcAVLvHGn/Emlwht/YAX4l+8hbfqUIAXh9CX/o4OLCQ3DZh3/e02Hp3K1Ck3/eg+tZ/KeHiDpj45BlPkX4ErP3Zl0BeHMQnyz80Z8Do19PQoWqB5k9h4g2VbCw9RkOaU3+6AurFhis2+vYPUD67XVuXXdz+7WyW7dvtwaRVYPL9LDaJ1CyIkOcx4iySI/zCHEWV9NRFk2XNepFeYQ4H+npJJozP66mZ4LPCnpWDz2rB8/y4Nke+nYfnu3Bs3q4J4c4Sx2cKQtBIeArAT8VOhTFBXphBnecwhmnsMcxrFEEeRMCN2PgZgR1MzJP+y4fkqQpCKROgMQZIh3eQz58B2lwhXB4id/fHyKxB0hkHzE8xIWNNG8/WyQtAf/MRX/o4nzooD/U0/2Bng7OywDVgzdweF8LEZ2kbV0u/AT6K3U+r5U9MOX79ycP9QPvWsNBrWzqtd1eVg9A85arh46Z5eohR7Yvd8T/yReqQJIniPN4ZqynY8Tp7xtlpl4xKUvzdHZ+6/qSqeCTqflPrL6NgIBne5PgY6Y9y0PP7uHKuarCUNv8ntWDJ3vwlIV+DPTCXIehKIcTZbDHCewwgRUmkKMIYhSiGI1Q3Nwgv3mL4maE4uYrFDe/Qn5zA2SL2yF8Hzi7RH7xLtLB15C9f4msf4HUHeogZfWRoIe4cBCnFqIEiCPVetO9EIA3mISlq3pwakz7QxcOzzgRES20le8uNJcGnze+VucFgCeLzmTxuwu3pww4SZEgyROkeYq0SHUoKeaHlEUhprVszvpu+1qOZbnSRc/qwbWmx21lruXCtVwddqweevbkzFCzbG44Ei6sKIUaj3XwKcdVCLoxIehtbfoGxdu3yEfT00gXn9kSrgs5GEAOBrDMWA4GUMEQmX+JzDtD6gyROj4SOQlKSW4hTiXiWCEKC0Q3KbJ0zplYAXi+Y4JTIzANJmGpP3TRP3PQ83m2iYhoGXf+3YXmkuFHQoinAH4O/WyspwdxqXBHlFJVmEnypAoz5bgedKp6hXldm99WlhbpwuXatrPuWZy68oyOa7lwpdsacHzHvzX8NEPSKoHJkQ6kuP0Socrz6SBUTr8dNcp+g2I00kFoKjyNEY9GCE2ZCsPl3hzHmQpFVhDAef998zqANRhCDgYQQQD4Q6TuAKnVRyI8pKKHOLcRZxJxWCAapYhuUoQ3ZjxKkb3Ngbftm+75Ct5AwAtsDO45uD9w4A305Tlv4MALnMn0gKGJiOgQbO3ThUqpzzF9ufBg/Pg/foyb9GYmEM2c8WkJTfPC0jbO2JRsYcOxHDjSmQo3juVU06504Tv+5LUJJFP1pTO9TK3OvIDTDDm2sLd6GUgpBZWmKEYjqDBEcTNGEYYoRq8ngWg0QjIeIzLTbaGoHp6WDkUApO9DBsFk8H047703U9b2Wvg+CtfXQQku4gSIxxmicYZolCIepYjGmR6PUsSjDNGXKaJxiiIrzxCHZphwPcuEIX3D99XXA33GyYSlybRrApQNeSJPRyYi6pK9PCfrrv3o336EV5H+1h9LWDpMSPvWwNJ3+lMBZ6qOKVsmGLWVOdZkGUc6sOT+P4qulIKKIuTj35sQNIYKxzrohKEej804HOtLa9XrSXkxHkM1yvUN2MuZCTxBsDgUBT6kX5sOAlhBANHvQ0gJVSgkUYZolCEe60A0KkNROR6niH+bIvplWSdGPBqhuOXx1bYr4QX6rJEX2Lh4z4f3RzZ6gQlLgYNeYKPnO1NnnKwD/0JTIiLajpMIWT/9659WoeoQAs0mVJahiKLqjE4VZka1kNMSiqryUT0UhSYU6elVHkMuXBey34cIfMi+r4NPvw/n3fcg/T6Eb8r7fT3P12NRvg4CE4xmQ1GbIi8QhxnicYYkzDAeZeZ1inicIX5Tvh4hGX+lzzCNdYCKx+mtTXM8C56vA5EXOAguBvACE5Zq5WVo0sHK5nOaiIjoVicRsoblE6N3TOU5ijCCikIdYMIQKop02InMdBjpM0FhpMvCUJdFoS4Lw/bpSL9e5sbquirU1MeBD+vqarp8KgQFprwRjPxJmBL2aj86SimkcV6FpCocfZkhHsWIwxHicYrEXI6bqjPOkMa3nw2TUlRnjdy+Dc+3cXbfm4Qj34SjwIHn29V9S73AhsVLcUREtAMnEbIAQBWFDjQm6KhwXAWXqaBTD0UtQUeFY1M30meAorI8hEqW++qeOtHrQXqeDjKeB+H3IT0daqz79/R034Pw+jr4VNNeFXraQpDs9289M7T6+6eQJjnCMEP86whJmCMJTRgygSgJTUAaTwekONRlt116A/S9Sq6vg1Kvb+Psfh+98rVv66Fvwy1f9yfzbFfykQJERHRQTiJk/edf/hWyL79ceTnhOPqyl+eZAGSm+x6sy0sTZDwThOrTZWiqTZchqj7teRDW7i85qUIhiadDUVIb9Ot8Mh21zIsyLPr6O2mLKiD1zNmi83f96rXr2/DMmaZeYFflPd+B61m8uZuIiDrlJELW5fe+B5VnkF4tCPl9CK8xXQUhH9LrrXxJbBeKQk0CUVQGnzmBKZpTHueLA5IUcPs23L6lQ5A5k6TLzBkkb3p+fV7Pt2E5PJtERERU2n+KuAP3/+7JnW5PKYUsKZBEGdIon5wZinKkkQlDUfl6Mj81Z5vKcRLlC+9FAswZpCoE6SB0/k5/KghNB6PZoGQzIBEREW3VSYSsZSilkGcFkjBHGk8ukaW1QNR8nZqzRM0AlUbZUh/UE1Lo+5A8G44Ze4GD4b2+DkK9SSBqO3tUBiZ+yo2IiOjwnETI+uyf/xvhTapDUZTXAlEtPIX5whuzAQAC+oyRZ8ExY9ezMLjoVUGpHLt9G07PhCQzruZ5Fi+vERERddhJhKx//9n/II3yqWDkeLV7jnoWnL5dnVVq1nM9qwpMTo9fjEtERESLnUTI+tsf/gUEv8eNiIiI7tBJfGaeAYuIiIju2kmELCIiIqK7xpBFREREtAMMWUREREQ7wJBFREREtAMMWUREREQ7wJBFREREtAMMWUREREQ7wJBFREREtAMMWUREREQ7wJBFREREtAMMWUREREQ7wJBFREREtANCKbXfHRDiNwB+eQebug/gt3ewnUN0ym0HTrv9bPvpOuX2n3LbgdNu/120/Q+VUu8sU3HvIeuuCCE+VUp9a9/7sQ+n3HbgtNvPtp9m24HTbv8ptx047fYfWtt5uZCIiIhoBxiyiIiIiHbglELWx/vegT065bYDp91+tv10nXL7T7ntwGm3/6DafjL3ZBERERHdpVM6k0VERER0ZxiyiIiIiHbA3vcOEBEREdUJIR4A+BaAlwCuAbxUSn2yjeXWXfc6OhOyhBCPoN+sD8z4mVLqJ0ssd1AHZBPmPXiz5A/iUwAvAHyqlHozp84FgMcAfgLgFYArAE8AvDi09q/Y9qXadSzHHVit/bX6c/tLV4+9qd+JPr/OPnal36/Z9s70+w0CyFH0eyHENYCnSqkPa2XPhRAvlVIvN1lu3XWvTSl19AOARwAe1F5fAPgCwOMFy11D//DUy54DuF6lziEMAB4CeA3g4ZL1XwBQc4Yvam2vl78G8Gjfbd1C2xe261iO+5rtX9hfOn7sj77Pr7uPXej3G7S9E/1+g/YfTb8H8LTl2DwE8HzT5dZd97pDV+7JulZKfV6+UPo/tKcAni1Y7klLnWdm2VXq7I0Q4loI8Qy6c7xaYdGXAD4E8E3o/2rK4QmAv6nV+xDAJYAPlFKXaomzg3dlg7YDi9t10Mcd2Kj9y/aXLh77o+/zxrr7ePT9Hpsdn6Pv91h/H4+p3z8C8Hmj7FNTvuly6657PftI4ltOvBcAPgNw0SgvE/ncdA+d4q9b1qdWqXMog9nXZf+jbz3Lh9n/apZa376HFdu+sF3HdNxXaf+y/aXDx74TfX7dfexCv9+g7Z3o9+vs4zH1+7ItzX018+b+TV9muXXXvclw9GeylE7j12ZYmrn2PPNfsFlf+Z/ywjrr7/l+KaVmHtgmhHjcVn5qOn7c1+ovXdCVPr/JPh57v9/l8enysT+yfn8FTNrUYl4blllu3XWvrRM3viulLluKH0LfDDvvRral3+wFdbZ/o9weCCEeQp8ybbo2N0sC+j17pQ7r0sG6bmvXMj8bR3vcV+gvXTv2XenzW/v5PMJ+v2nbj73fr72PR9TvL3a43LrrXlsnQtYcTwD88Jb5B3lA9uiBUurvG2WvAEBNf/rkuRACB/ILd12L2nVKx73U7C9dPPZd6fPb3Mdj6/ebtL0L/X7b+3gK/X6vjv5yYRshxGPo9N385UEtzH8tM//9KKXetFxGOLSbQFfW1Xatq62/8D3qPvZ7AB1o17oOvd+bS6M7WW7dda/joM5krXjN+1XbKVOzjidKqW8uuc2LW069Ll1nU9to+wZ+AODbS9Z9CX06eWvvyZ7bXqraVRbcxXE329lb+1fsL5049ofS58121m7/FvbxaPv9FvfjKPv9pvu4736/xPYAfalyqs2N+ess92qJOlt1MCGrfEDYCov8HEDbmaqnWO4Xx8EckC22fZ1tX0BfMmgLrN+fdykB+tp/82Ow62z/ztu+RLvW7eTr7Mvejr3R2l86euwPps+bda7b/o1/Po+436/d9o70+23t4976/SJKqTdCiJeYvTR6hVvus152uXXWvZFtf1xxnwP0ac2lP4IJ/VHYB42yawCvV6lzKANW+Ch7bZlHbW3BnEdg1MpnPgJ7DG1ftl3HdNw3OPat/aWrx37Z43oMx37TfTzmfr9O27vU77dw7A++35t9fNwoe4TFDyNduNy661536Mw9Web68lNVS6JCiIcLTsl+Av3VBHUPTPkqdY7Zn6H9voyX0KeTm/MeAvhc3c1p461boV2dPu639ZeuHnujK31+03085n6/cts71u/X3scj6vcfYfrhuIC+Sf+j8oUQ4kII8YVp09LLLVlne+4qme449T6C/r6lB7XhIfT3MpV12r5C4AKzX0/wAtNfsbGwzqEMmPMVCG1tr8173mxf431tvhefofFf1CEMq7R9mXYd03Ffs/2L+ktXj30n+vwK7ehcv1+37V3p9xu2/2j6vdm/p2afvo/G2Wqzb6+bP9+Lllu2zrYGYTZ4tMy16NdzZr9USn1Qq/cLAB+p2qcnhP6ize9CX/O+hk7sbV8We2udfTHt+gH0fpWfFvoEuhP+pFZnpu1m3jMAUEo9mbP+8gtF70H/UE/9F7RPm7R9mXYd8nEH1mv/sv3FLNvVY3/Ufb60aB+72u+B9dvehX4PrN7+rvT7Y3T0IYuIiIjoEHXmniwiIiKiQ8KQRURERLQDDFlEREREO8CQRURERLQDDFlEREREO8CQRURERLQDDFlEREREO8CQRURERLQDDFlEREREO8CQRURERLQDDFlEREREO8CQRURERLQDDFlEdPSEEI+EEM+EEC+EEI9a5l8IIV4LIR7uY/+I6DQxZBHRURNCPABwpZR6AuC5GZq+A+ACwKu73DciOm0MWUR07J4opT4209+cU+dDAFBKfX43u0RExJBFREdMCHEB4Ita0XcAfNJS9SEABiwiulMMWUR0tJRSbwB8DFSXDS/QuFxogtgF2sMXEdHOMGQR0VEzQQsAvmvGP25UKW92f3E3e0REpDFkEVFXPALweS10lcr7sXgmi4juFEMWER09IcQ1gGsscT+WuXxIRLRzDFlE1AXXZvzFnHmfAlXAenxXO0VEp40hi4i6oHz+1ct6oRCiDFSfmfFD8AZ4IrojDFlEdPTM86/ewNx/BVSfNvzAvCxD2Id8VhYR3RWhlNr3PhARbczcl/UMk7NZnymlPjZfs/NdU/5PDFlEdFcYsoiIiIh2gJcLiYiIiHaAIYuIiIhoBxiyiIiIiHaAIYuIiIhoBxiyiIiIiHaAIYuIiIhoBxiyiIiIiHaAIYuIiIhoBxiyiIiIiHaAIYuIiIhoB/4fUEBsW2kz3TMAAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "plt.rc(\"text\", usetex=True)\n", + "import numpy as np\n", + "\n", + "x = np.linspace(-2, 0, 100000, endpoint=False)\n", + "plt.figure(figsize=(10, 6))\n", + "for t in [0.1, 0.5, 1, 1.5, 2]:\n", + " plt.plot(x, -t * np.log(-x), label=r\"$t = \" + str(t) + \"$\")\n", + "plt.legend(fontsize=20)\n", + "plt.xticks(fontsize=20)\n", + "plt.yticks(fontsize=20)\n", + "_ = plt.xlabel(\"$u$\", fontsize=26)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### \"Constrained\" problem\n", + "\n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "& \\min f_0(x) + \\sum_{i=1}^m -t \\log(-f_i(x))\\\\\n", + "\\text{s.t. } & Ax = b,\n", + "\\end{split}\n", + "\\end{equation*}\n", + "\n", + "- The problem is still **convex**\n", + "- Function \n", + "\n", + "$$\n", + "\\phi(x) = -\\sum\\limits_{i=1}^m \\log(-f_i(x))\n", + "$$ \n", + "\n", + "is called *logarithmic barier*. \n", + "\n", + "Its domain is a set of points such that the inequality constraints are strictly feasible.\n", + "\n", + "**Exercise.** Find gradiebnt and hessian of $\\phi(x)$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Central path\n", + "\n", + "For every $t > 0$ \"constrained\" problem has unique solution $x^*(t)$.\n", + "\n", + "**Definition.** Sequence $x^*(t)$ for $t > 0$ is formed *central path*." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Optimality conditions for \"constrained\" problem\n", + "\n", + "- Primal feasibility\n", + "\n", + "$$\n", + "Ax^*(t) = b, \\; f_i(x^*) < 0, \\; i = 1,\\ldots,m\n", + "$$\n", + "\n", + "- Lagrangian stationarity\n", + "\n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "& f'_0(x^*(t)) + \\phi'(x^*(t)) + A^{\\top}\\hat{\\mu} = \\\\\n", + "& = f'_0(x^*(t)) - t\\sum_{i=1}^m \\frac{f_i'(x^*(t))}{f_i(x^*(t))} + A^{\\top}\\hat{\\mu} = 0\n", + "\\end{split}\n", + "\\end{equation*}" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "- Introduce the following notation\n", + "\n", + "$$\n", + "\\lambda^*_i(t) = -\\frac{t}{f_i(x^*(t))} \\; i=1,\\ldots,m \\text{ и } \\mu^* = \\hat{\\mu}\n", + "$$\n", + "\n", + "- Then optimality condition can be re-written as\n", + "\n", + "$$\n", + "f'_0(x^*(t)) + \\sum_{i=1}^m \\lambda^*_i(t)f_i'(x^*(t)) + A^{\\top}\\mu^* = 0\n", + "$$\n", + "\n", + "- Then $x^*(t)$ is minimizer of the following Lagrangian\n", + "\n", + "$$\n", + "L = f_0(x) + \\sum_{i=1}^m \\lambda_if_i(x) + \\mu^{\\top}(Ax - b)\n", + "$$\n", + "\n", + "where $\\lambda = \\lambda^*(t)$ and $\\mu = \\mu^*$." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Duality gap\n", + "\n", + "- Dual function $g(\\lambda^*(t), \\mu^*)$ is finite and is represented in the following way\n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "g(\\lambda^*(t), \\mu^*) & = f_0(x^*(t)) + \\sum_{i=1}^m \\lambda^*_i(t)f_i(x^*(t)) + (\\mu^*)^{\\top}(Ax^*(t) - b)\\\\\n", + "& = f_0(x^*(t)) - mt\n", + "\\end{split}\n", + "\\end{equation*}\n", + "\n", + "- Duality gap\n", + "\n", + "$$\n", + "f_0(x^*(t)) - p^* \\leq mt\n", + "$$\n", + "\n", + "- While $t \\to 0$ duality gap is 0 and central path converges to the solution of the original problem" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## KKT interpretation\n", + "\n", + "Optimality conditions for \"constrained\" problem is equivalent to optimality conditions for original problem if\n", + "\n", + "$$\n", + "-\\lambda_i f_i(x) = 0 \\Rightarrow - \\lambda_i f_i(x) = t \\quad i = 1,\\ldots, m\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Physical interpretation\n", + "\n", + "- Assume that we do not have equality constraints \n", + "- Consider classical particle in the field of forces\n", + "- Every constraint $f_i(x) \\leq 0$ corresponds to the following force\n", + "\n", + "$$\n", + "F_i(x) = -\\nabla(-\\log(-f_i(x))) = \\frac{f'_i(x)}{f_i(x)}\n", + "$$\n", + "\n", + "- Objective function corresponds to some force, too\n", + "\n", + "$$\n", + "F_0(x) = -\\frac{f'_0(x)}{t}\n", + "$$\n", + "\n", + "- Every point from the central path $x^*(t)$ is a equilibrium state of particle where sum of forces is zero\n", + "- While decreasing $t$ forces $F_0(x)$ dominates forces $F_i(x)$ and particle aims to get state which corresponding to optimal value of objective\n", + "- As far as forces $F_i(x)$ go to infinity when particle is close to the bounds, particle will never leave feasible region" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Barier method\n", + "\n", + "- $x_0$ has to be feasible\n", + "- $t_0 > 0$ is initial value of parameter\n", + "- $\\alpha \\in (0, 1)$ is multiplier for decreasing $t_0$\n", + "\n", + "```python\n", + "def BarrierMethod(f, x0, t0, tol, alpha, **kwargs):\n", + " \n", + " x = x0\n", + " \n", + " t = t0\n", + " \n", + " while True:\n", + " \n", + " x = SolveBarrierProblem(f, t, x, **kwargs)\n", + " \n", + " if m * t < tol:\n", + " \n", + " break\n", + " \n", + " t *= alpha\n", + " \n", + " return x\n", + "\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Parameters selection\n", + "\n", + "- Multiplier $\\alpha$\n", + " - In the case of $\\alpha \\sim 1$, \"constrained\" problem is solved after **small** number of iterations, but central path consists of **large** number of points\n", + " - In the case $\\alpha \\sim 10^{-5}$ the situation is completely different: **large** number of iterations for solving \"constrained\" problems, but **small** number of points for central path\n", + "- Initialization $t_0$\n", + " - Alternatives are similar to parameter $\\alpha$\n", + " - Parameter $t_0$ affects the initial point in central path" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Almost convergence theorem\n", + "\n", + "- As it was shown above, while $t \\to 0$ barrier method converges to the solution of the original problem\n", + "- Convergence speed is directly affected by parameters $\\alpha$ and $t_0$, as was shown above\n", + "- The main difficulty is fast solving auxilliary problems with Newton methods" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Problem of finding feasible starting point\n", + "\n", + "- Barier method requires feasible starting point\n", + "- Method is splitted in two phases\n", + " - The first phase gives feasible starting point\n", + " - The second phase uses this starting point to run barier method" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### The first phase method\n", + "\n", + "Simple method to find strictly feasible point\n", + "\n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "& \\min s\\\\\n", + "\\text{s.t. } & f_i(x) \\leq s\\\\\n", + "& Ax = b\n", + "\\end{split}\n", + "\\end{equation*}\n", + "\n", + "- this problem always has strictly feasible starting point\n", + "- if $s^* < 0$, then $x^*$ is strictly feasible and can be used in barier method\n", + "- if $s^* > 0$, then original problem is infeasible and feasible set is empty" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Sum of inconsistencies\n", + "\n", + "\\begin{equation*}\n", + "\\begin{split}\n", + "& \\min s_1 + \\ldots + s_m\\\\\n", + "\\text{s.t. } & f_i(x) \\leq s_i\\\\\n", + "& Ax = b\\\\\n", + "& s \\geq 0\n", + "\\end{split}\n", + "\\end{equation*}\n", + "\n", + "- optimal objective value is 0 and it attains in the case of consistency of inequality constraints\n", + "- if the problem is infeasible, it is possible to identify what constraints are violated, i.e. $s_i > 0$ " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Second phase\n", + "\n", + "- After starting point $x_0$ is found, one can run standard Newton method for equality constrained problem" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Primal-dual method\n", + "\n", + "It is similar to barier method, but\n", + "- every iteration updates both primal and dual variables\n", + "- Newton direction is taken from the modified KKT system\n", + "- iterands in primal-dual method may be not feasible\n", + "- it works even if the problem is not strictly feasible" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Recap\n", + "\n", + "- Newton method for convex optimization problem with equality constraint\n", + "- The case of infeasible initial point\n", + "- Primal barier method\n", + "- Primal-dual method" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Python 3 (cvxpy)", + "language": "python", + "name": "cvxpy" + }, + "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.6.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}