18-LinProgPrimalInterior/Seminar18.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Семинар 18\n",
    "\n",
    "# Линейное программирование. Прямой метод внутренней точки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## На прошлом семинаре...\n",
    "\n",
    "- Постановки задачи линейного программирования\n",
    "- Примеры приложений\n",
    "- Симплекс-метод и его недостатки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Обзор существующих библиотек для решения задачи линейного программирования\n",
    "\n",
    "- Решение задачи линейного прграммирования стало полноценной технологией, поэтому обычно нужно вызвать правильную функцию на правильном языке и получить ответ. \n",
    "- Обзор можно найти [тут](https://en.wikipedia.org/wiki/Linear_programming#Solvers_and_scripting_.28programming.29_languages)\n",
    "- Сравнение можно найти [тут](http://prod.sandia.gov/techlib/access-control.cgi/2013/138847.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Немного истории\n",
    "\n",
    "- В 1979 г. Л. Хачиян [предложил](http://www.mathnet.ru/links/18409a6ec8cd239985969a9e98c2f6ae/zvmmf5239.pdf) метод эллипсоидов и показал, что он решает любую задачу линейного программирования за **полиномиальное время**. Сообщение об этом появилось на [первой полосе New York Times 7 ноября 1979](http://www.nytimes.com/1979/11/07/archives/a-soviet-discovery-rocks-world-of-mathematics-russians-surprise.html?_r=0). Популярное изложение успехов советской школы математического программирования в статье [Коммерсанта](https://www.kommersant.ru/doc/2679807).\n",
    "- Однако сложность метода эллипсоидов $O(n^6 L^2)$, и он проигрывал симплекс-методу при решении реальных задач\n",
    "- В 1984 г. Н. Кармаркар [предложил](http://retis.sssup.it/~bini/teaching/optim2010/karmarkar.pdf) другой полиномиальный метод решения задачи линейного программирования, который был значительно быстрее метода эллипсоидов, а именно $O(n^{3.5} L^2)$.  \n",
    "- Одна из последних [работ](https://arxiv.org/pdf/1312.6677.pdf) по этой теме предлагает алгоритм решения задачи линейного программирования за $\\tilde{O}\\left(\\sqrt{rank(A)}L\\right)$. \n",
    "$f(n)\\in\\tilde{O}(g(n)) \\leftrightarrow \\exists k:f(n) \\in O(g(n) \\log^kg(n))$.\n",
    "- И метод эллипсоидов, и метод Кармаркара относятся к прямо-двойственным методам или методам внутренней точки, которые будут освещены далее..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Альтернатива симплекс-методу"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Симплекс-метод основан на подходе поиска решения среди вершин\n",
    "- Однако можно искать решение по-другому: двигаться по точкам внутри области к одной из вершин - решению задачи\n",
    "- Поэтому такие методы называют *методами внутренней точки* "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## \"Дёшево, но много\" vs. \"дорого, но мало\"\n",
    "\n",
    "- Один из основных водоразделов в методах оптимизации\n",
    "- Линейное программирование решается за полиномиальное время, эксплуатируя стратегию \"дорого, но мало\"\n",
    "- Далее будут показаны задачи и методы, которые работают по стратегии \"дёшево, но много\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Двойственность в задаче линейного программирования\n",
    "\n",
    "**Теорема.**\n",
    "\n",
    "- Если прямая (двойственная) задачи имеет конечное решение, то конечное решение имеет и двойственная (прямая).\n",
    "- Если прямая (двойственная) задача неограничена, то допустимое множество двойственной (прямой) задачи пусто."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Идея методов внутренней точки\n",
    "\n",
    "Исходная задача\n",
    "\\begin{align*}\n",
    "&\\min_x c^{\\top}x \\\\\n",
    "\\text{s.t. } & Ax = b\\\\\n",
    "& x_i \\geq 0, \\; i = 1,\\dots, n\n",
    "\\end{align*}\n",
    "\n",
    "Аппроксимированная задача\n",
    "\\begin{align*}\n",
    "&\\min_x c^{\\top}x {\\color{red}{- \\mu \\sum\\limits_{i=1}^n \\ln x_i}} \\\\\n",
    "\\text{s.t. } & Ax = b\\\\\n",
    "\\end{align*}\n",
    "для некоторого $\\mu > 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Барьерная функция\n",
    "\n",
    "**Определение.** Функция $B(x, \\mu) = -\\mu\\ln x$ называется *барьерной* для задачи с ограничением $x \\geq 0$.\n",
    "\n",
    "Более подробно о таких функциях будет рассказано в контексте нелинейной условной оптимизации..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Что произошло?\n",
    "\n",
    "- Сделали из линейной задачу нелинейную\n",
    "- Перенесли ограничение типа неравенства в целевую функцию\n",
    "- Ввели дополнительный параметр $\\mu$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Почему это хорошо?\n",
    "\n",
    "Переход к задаче с ограничениями типа равенств $\\to$ упрощение условий оптимальности, в частности\n",
    "\n",
    "- Исключено требование дополняющей нежёсткости\n",
    "- Исключено условие неотрицательности множителя Лагранжа для ограничения типа неравенства"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Условия оптимальности\n",
    "\n",
    "- Лагранжиан: $L = c^{\\top}x - \\mu\\sum\\limits_{i=1}^n \\ln x_i + \\lambda^{\\top}(Ax - b)$\n",
    "- Стационарная точка $L$: \n",
    "$$\n",
    "c - \\mu X^{-1}e + A^{\\top}\\lambda = 0,\n",
    "$$\n",
    "где $X = \\mathrm{diag}(x_1, \\dots, x_n)$ и $e = [1, \\dots, 1]$\n",
    "- Ограничение типа равенства: $Ax = b$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Пусть $s = \\mu X^{-1}e$, тогда условия оптимальности можно переписать так:\n",
    "- $A^{\\top}\\lambda + c - s = 0 $\n",
    "- $Xs = {\\color{red}{\\mu e}}$\n",
    "- $Ax = b$\n",
    "\n",
    "Также $x > 0 \\Rightarrow s > 0$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сравнение с условиями оптимальности для исходной задачи\n",
    "\n",
    "- Лагранжиан: $L = c^{\\top}x + \\lambda^{\\top}(Ax - b) - s^{\\top}x$\n",
    "- Условие стационарности: $c + A^{\\top}\\lambda - s = 0$\n",
    "- Допустимость прямой задачи: $Ax = b, \\; x \\geq 0$\n",
    "- Допустимость двойственной: $s \\geq 0$\n",
    "- Условие дополняющей нежёсткости: $s_ix_i = 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### После упрощения\n",
    "\n",
    "- $A^{\\top}\\lambda + c - s = 0$\n",
    "- $Ax = b$\n",
    "- $Xs = {\\color{red}{0}}$\n",
    "- $x \\geq 0, \\; s \\geq 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Вывод\n",
    "\n",
    "- Введение барьерной функции c множителем $\\mu$ эквивалентно релаксации условий дополняющей нежёсткости на параметр $\\mu$\n",
    "- При $\\mu \\to 0$ решения задач совпадают!\n",
    "- Идея: итеративно решать задачи с барьерной функцией, уменьшая $\\mu$. Последовательность решений сойдётся к вершине симплекса по траектории из точек, лежащих внутри симплекса."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Общая схема\n",
    "```python\n",
    "def GeneralInteriorPointLP(c, A, b, x0, mu0, rho, tol):\n",
    "    x = x0\n",
    "    mu = mu0\n",
    "    e = np.ones(c.shape[0])\n",
    "    while True:\n",
    "        primal_var, dual_var = StepInsideFeasibleSet(c, A, b, x, mu)\n",
    "        mu *= rho\n",
    "        if converge(primal_var, dual_var, c, A, b, tol) and mu < tol:\n",
    "            break\n",
    "    return x\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Как решать задачу с барьерной функцией?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Прямой метод - следующий слайд\n",
    "- Прямо-двойственный метод - через несколько недель"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Прямой метод\n",
    "\n",
    "Вспомним исходную задачу:\n",
    "\\begin{align*}\n",
    "&\\min_x c^{\\top}x - \\mu \\sum\\limits_{i=1}^n \\ln x_i \\\\\n",
    "\\text{s.t. } & Ax = b\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Идея: приблизим целевую функцию до второго порядка, как в методе Ньютона."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Реализация\n",
    "На $(k+1)$-ой итерации необходимо решить следующую задачу: \n",
    "\\begin{align*}\n",
    "&\\min_p \\frac{1}{2}p^{\\top}Hp + g^{\\top}p\\\\\n",
    "\\text{s.t. } & A(x_k + p) = b,\\\\\n",
    "\\end{align*}\n",
    "где $H = \\mu X^{-2}$ - гессиан, и $g = c - \\mu X^{-1}e$ - градиент."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Снова KKT\n",
    "\n",
    "Выпишем условия ККТ для этой задачи\n",
    "- $Hp + g + A^{\\top}\\lambda = 0$\n",
    "- $Ap = 0$\n",
    "\n",
    "или\n",
    "$$\\begin{bmatrix}\n",
    "H & A^{\\top}\\\\\n",
    "A & 0\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "p\\\\\n",
    "\\lambda\n",
    "\\end{bmatrix}\n",
    "= \n",
    "\\begin{bmatrix}\n",
    "-g\\\\\n",
    "0\n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Из первой строки:\n",
    "$$\n",
    "-\\mu X^{-2}p + A^{\\top}\\lambda = c - \\mu X^{-1}e\n",
    "$$\n",
    "$$\n",
    "-\\mu Ap + AX^{2}A^{\\top}\\lambda = AX^2c - \\mu AXe\n",
    "$$\n",
    "$$\n",
    "AX^{2}A^{\\top}\\lambda = AX^2c - \\mu AXe\n",
    "$$\n",
    "Так как $X \\in \\mathbb{S}^n_{++}$ и $A$ полного ранга, то уравнение имеет единственное решение $\\lambda^*$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Найдём направление $p$\n",
    "\n",
    "$$\n",
    "-\\mu p + X^2A^{\\top}\\lambda^* = X^2c - \\mu Xe = X^2c - \\mu x\n",
    "$$\n",
    "$$\n",
    "p = x + \\frac{1}{\\mu}X^2(A^{\\top}\\lambda^* - c)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Способы решения системы из ККТ\n",
    "\n",
    "1. Прямой способ: формирование матрицы $(n + m) \\times (n + m)$ и явное решение системы - $\\frac{1}{3}(n + m)^3$\n",
    "2. Последовательное исключение переменных:\n",
    "    - $Hp + A^{\\top}\\lambda = -g$, $p = -H^{-1}(g + A^{\\top}\\lambda)$\n",
    "    - $Ap = -AH^{-1}(g + A^{\\top}\\lambda) = -AH^{-1}A^{\\top}\\lambda - AH^{-1}g = 0$\n",
    "    \n",
    "    Здесь матрица $-AH^{-1}A^{\\top}$ есть *дополнение по Шуру* матрицы $H$.\n",
    "3. Алгоритм вычисления решения при последовательном исключении переменных\n",
    "    - Вычислить $H^{-1}g$ и $H^{-1}A^{\\top}$ - $f_H + (m+1)s_H$ операций\n",
    "    - Вычислить дополнение по Шуру $-AH^{-1}A^{\\top}$ - $\\mathcal{O}(m^2n)$\n",
    "    - Найти $\\lambda$ - $\\frac{1}{3}m^3$ операций\n",
    "    - Найти $p$ - $s_H + \\mathcal{O}(mn)$ операций\n",
    "4. Итого: $f_H + ms_H + \\frac{m^3}{3} + \\mathcal{O}(m^2n)$ уже гораздо быстрее прямого способа"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Используем структуру матрицы $H$\n",
    "\n",
    "- В нашем случае $H = \\mu X^{-2}$ - диагональная матрица!\n",
    "- $f_H$ - $n$ операций\n",
    "- $s_H$ - $n$ операций\n",
    "- Итоговая сложность $\\frac{m^3}{3} + \\mathcal{O}(m^2n)$ операций, где $m \\ll n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Поиск шага $\\alpha$\n",
    "\n",
    "- Обычный линейный поиск с условиями достаточного убывания\n",
    "- Условие $A(x_k + \\alpha p) = b$ выполняется автоматически"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Псевдокод прямого барьерного метода\n",
    "\n",
    "```python\n",
    "def PrimalBarrierLP(c, A, b, x0, mu0, rho, tol):\n",
    "    x = x0\n",
    "    mu = mu0\n",
    "    e = np.ones(x.shape[0])\n",
    "    while True:\n",
    "        p, lam = ComputeNewtonDirection(c, x, A, mu)\n",
    "        alpha = line_search(p, mu, c, x)\n",
    "        x = x + alpha * p\n",
    "        mu = rho * mu\n",
    "        if mu < tol and np.linalg.norm(x.dot(c - A.T.dot(lam)) - mu * e) < tol:\n",
    "            break\n",
    "    return x\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сравнение барьерного метода и прямого метода внутренней точки\n",
    "\n",
    "- Пример Klee-Minty c прошлого семинара\n",
    "\\begin{align*}\n",
    "& \\max_{x \\in \\mathbb{R}^n} 2^{n-1}x_1 + 2^{n-2}x_2 + \\dots + 2x_{n-1} + x_n\\\\\n",
    "\\text{s.t. } & x_1 \\leq 5\\\\\n",
    "& 4x_1 + x_2 \\leq 25\\\\\n",
    "& 8x_1 + 4x_2 + x_3 \\leq 125\\\\\n",
    "& \\ldots\\\\\n",
    "& 2^n x_1 + 2^{n-1}x_2 + 2^{n-2}x_3 + \\ldots + x_n \\leq 5^n\\\\\n",
    "& x \\geq 0\n",
    "\\end{align*}\n",
    "- Какая сложность работы симплекс метода?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Сведение к стандартной форме\n",
    "$$\n",
    "\\begin{align*}\n",
    "& \\min_{x, \\; z} -c^{\\top}x \\\\\n",
    "\\text{s.t. } & Ax + z = b\\\\\n",
    "& x \\geq 0, \\quad z \\geq 0\n",
    "\\end{align*}\n",
    "$$\n",
    "- Сравним скорось работы прямого барьерного метода и симплекс-метода"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.optimize as scopt\n",
    "import scipy.linalg as sclin\n",
    "plt.rc(\"text\", usetex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def NewtonLinConstraintsFeasible(f, gradf, hessf, A, x0, line_search, linsys_solver, args=(), \n",
    "                                 disp=False, disp_conv=False, callback=None, tol=1e-6, max_iter=100, **kwargs):\n",
    "    x = x0.copy()\n",
    "    n = x0.shape[0]\n",
    "    iteration = 0\n",
    "    lam = np.random.randn(A.shape[0])\n",
    "    while True:\n",
    "        gradient, hess = gradf(x, *args), hessf(x, *args)\n",
    "        h = linsys_solver(hess, A, gradient)\n",
    "        descent_dir = h[:n]\n",
    "        decrement = descent_dir.dot(hessf(x, *args).dot(descent_dir))\n",
    "        if decrement < tol:\n",
    "            if disp_conv:\n",
    "                print(\"Tolerance achieved! Decrement = {}\".format(decrement))\n",
    "            break\n",
    "        alpha = line_search(x, descent_dir, f, gradf, args, **kwargs)\n",
    "        if alpha < 1e-16:\n",
    "            if disp_conv:\n",
    "                print(\"Step is too small!\")\n",
    "        x = x + alpha * descent_dir\n",
    "        if callback is not None:\n",
    "            callback((descent_dir, x))\n",
    "        iteration += 1\n",
    "        if disp:\n",
    "            print(\"Current function val = {}\".format(f(x, *args)))\n",
    "            print(\"Newton decrement = {}\".format(decrement))\n",
    "        if iteration >= max_iter:\n",
    "            if disp_conv:\n",
    "                print(\"Maxiter exceeds!\")\n",
    "            break\n",
    "    res = {\"x\": x, \"num_iter\": iteration, \"tol\": decrement}\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def simple_solver(hess, A, gradient):\n",
    "    n = hess.shape[0]\n",
    "    n_lin_row, n_lin_col = A.shape\n",
    "    modified_hess = np.zeros((n + n_lin_row, n + n_lin_row))\n",
    "    modified_hess[:n, :n] = hess\n",
    "    modified_hess[n:n + n_lin_row, :n_lin_col] = A\n",
    "    modified_hess[:n_lin_col, n:n + n_lin_row] = A.T\n",
    "    rhs = np.zeros(n + n_lin_row)\n",
    "    rhs[:n] = -gradient\n",
    "    h = np.linalg.solve(modified_hess, rhs)\n",
    "    return h\n",
    "\n",
    "def elimination_solver(hess, A, gradient):\n",
    "    inv_hess_diag = np.divide(1.0, np.diag(hess))\n",
    "    inv_hess_grad = np.multiply(-inv_hess_diag, gradient)\n",
    "    rhs = A.dot(inv_hess_grad)\n",
    "    L_inv_hess = np.sqrt(inv_hess_diag)\n",
    "    AL_inv_hess = A * L_inv_hess\n",
    "    # print(AL_inv_hess.shape)\n",
    "    S = AL_inv_hess.dot(AL_inv_hess.T)\n",
    "    cho_S = sclin.cho_factor(S)\n",
    "    w = sclin.cho_solve(cho_S, rhs)\n",
    "#     w = np.linalg.solve(S, rhs)\n",
    "    v = np.subtract(inv_hess_grad, np.multiply(inv_hess_diag, A.T.dot(w)))\n",
    "#     h = np.zeros(hess.shape[1] + A.shape[0])\n",
    "#     h[:hess.shape[1]] = v\n",
    "#     h[hess.shape[1]:hess.shape[1] + A.shape[0]] = w\n",
    "    return v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def backtracking(x, descent_dir, f, grad_f, args, **kwargs):\n",
    "    beta1 = kwargs[\"beta1\"]\n",
    "    rho = kwargs[\"rho\"]\n",
    "    alpha = 1\n",
    "    while f(x + alpha * descent_dir, *args) >= f(x, *args) + beta1 * alpha * grad_f(x, *args).dot(descent_dir) \\\n",
    "            or np.isnan(f(x + alpha * descent_dir, *args)):\n",
    "        alpha *= rho\n",
    "        if alpha < 1e-16:\n",
    "            break\n",
    "    return alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def generate_KleeMinty_test_problem(n):\n",
    "    c = np.array([2**i for i in range(n)])\n",
    "    c = -c[::-1]\n",
    "    bounds = [(0, None) for i in range(n)]\n",
    "    b = np.array([5**(i+1) for i in range(n)])\n",
    "    a = np.array([1] + [2**(i+1) for i in range(1, n)])\n",
    "    A = np.zeros((n, n))\n",
    "    for i in range(n):\n",
    "        A[i:, i] = a[:n-i]\n",
    "    return c, A, b, bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n = 7\n",
    "c, A, b, _ = generate_KleeMinty_test_problem(n)\n",
    "eps = 1e-10\n",
    "def f(x, c, mu):\n",
    "    n = c.shape[0]\n",
    "    return c.dot(x[:n]) - mu * np.sum(np.log(eps + x))\n",
    "\n",
    "def gradf(x, c, mu):\n",
    "    grad = np.zeros(len(x))\n",
    "    n = c.shape[0]\n",
    "    grad[:n] = c - mu / (eps + x[:n])\n",
    "    grad[n:] = -mu / (eps + x[n:])\n",
    "    return grad\n",
    "\n",
    "def hessf(x, c, mu):\n",
    "    return mu * np.diag(1. / (eps + x)**2)\n",
    "\n",
    "A_lin = np.zeros((n, n + A.shape[0]))\n",
    "A_lin[:n, :n] = A\n",
    "A_lin[:n, n:n + A.shape[0]] = np.eye(A.shape[0])\n",
    "mu = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Проверим верно ли вычисляется градиент"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.207673676669207e-07"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scopt.check_grad(f, gradf, np.random.rand(n), c, mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Выбор начального приближения допустимого по ограничениям и области определения целевой функции"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.689582056016661e-13\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "x0 = np.zeros(2*n)\n",
    "x0[:n] = np.random.rand(n)\n",
    "x0[n:2*n] = b - A.dot(x0[:n])\n",
    "print(np.linalg.norm(A_lin.dot(x0) - b))\n",
    "print(np.sum(x0 <= 1e-6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Проверим сходимость"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Decrement value = 4.862770418574955e-09\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAFgCAYAAACmKdhBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl4XOWVJ/7vqVVbqbTa2izbsvCOjeUyGAgQQCRkcSckBtNZJr2ATE//kl4mwU3Smd+keya0Sff8ep7uNG2RTBbSnQCGJO0EmtgkBBJwQJbB+4LlXbKtxdrXqjq/P+4tuSyXbC2lureqvp/n0aO6S1UdIPHrc9/znldUFURERERERDR9DqsDICIiIiIiShVMsIiIiIiIiOKECRYREREREVGcMMEiIiIiIiKKEyZYREREREREccIEi4iIiIiIKE6YYBEREREREcUJEywiIiIiIqI4YYJFREREREQUJy6rA7CLoqIinTdvntVhEBGlpV27drWparHVcSQDjldERNaY6FjFBMs0b948NDQ0WB0GEVFaEpGTVseQLDheERFZY6JjFUsEiYiIiIiI4oQJFhERERERUZwwwSIiIkoCIrJOROq7urqsDoWIiK6CCRYREVESUNVtqlrn9/utDoWIiK6CCRYREZFFRGS7iKiIXBSRR62Oh4iIpo9dBImIiKyzBcD9qtppdSBERBQfnMEiIiKyEJMrIqLUwhksIiIi61SJSC2APAAFqlpvdUBERDQ9TLCIiIimSURqADymqvfHuFYHoMM8rFLVJ6Iu10dmsMz1WDtUtWnmIyYiopnCBIuIiGiKzMRqg3lYFeN6HQCo6lbzuEpEtqjqRvN8dHlgE4BaAJzFIiJKYkyw4uB4Wx++89vj8Dgd8Lod8Did8Lgc8LgccDkEIoCIQAA4xDgGAAGMaxDjYPScca/LKchwO5HlMX4y3MZP9OcILn1exOjnm58jcun+giwPXE4uvSMiigdVbQTQaCZatTFu2aiqq6PubzJLAmH+vj+SbCXCs2+fxoGWbmR6nMh0O5HhdiDTbYxZDhHjx2GMGWaMkeHpsjEnMrbAPM50O5HtdSHH60K214kcrwtup+PSWDXmvZHPNj4Bo+MkEVEqYIIVB609Q/iPd5sxHAxjOBhGMKxWhzQur8uBRSU+LCnJxZJSH5aW+VHs80YNlmbCF3U8VnTCluV2wZ/lTkzwRERJRETyANTEuNRpJldNAJ6LOh8AsGkmY9p18iJe3NeCwZEQRkL2G6siCViG24klpbm4vtyPFRXGz/yiHDgdTMKIyP5E1X5/wFohEAhoQ0NDXD4rFFYMB8MYCoYQViCsClVAVRHJvRTmORjnASDynyLyeyQcxsBwCAMjodHfgyMhAJc+M6yX3m987qUXl38HEAqHcbK9HwfPdeNgSw86+oan/c8qAnzk+lL82d3X4brZvml/HhGlJxHZpaoBq+OYKnMG66no2Srz3Cuqmj/m3u0AtqvqEyKyHkABjCYXjaq641rfFa/xKhgKYzBojDNDwZA5ThnjS3h0vNJLYxNgjik6em/k9eBICL1DQfQNhdA3FETvUBAjoTAUl4+B0Z8VMXodxheEFegdCmJ/cxf2ne3GgDnueV0OZLido++LrgYZS0TgMGfYHOYM2oLiHPzx++bjzkWz4GCiRkRTMNGxijNYM8DpEKP8wuO89s0WUVWc7x7CwZZudA4MX5HcxUq7VdUcAI1BNazAibY+/GDnSfx8bws+uqIMf3Z3NapnMdEiIoKROHXEON8JoBC4tDbrWsy1XHUAUFlZGZfgXE4HcpwO5Hjt+1eBUFhxrLUXe8504fC57tFZt9EHkzHeE534RR5ChsKK377Xhj/+XgOqZ+Xg4dvm42M3lF+WsBERxYt9/1SlGSUiKPFnoMSfMe3P2njHAnzr9SZ8940T+NmeZnx0RRk+vLwEq+flY5Zv+p9PRJTuVLVeRFoArPN4PKuv+YYU4XQIFs72YWEcKiRGQmH8fE8L6l9rwqbn9+IbLx/BA4EKrJ6bj+vL/ZiVy/GKiOKDCRZNW0G2B4/euxgP3VaFp15vwvffOIFt7zYDACoLsrB6bj5Wz83Hx24ogy+D67WIKK0UxDiXB6B9sh+kqtsAbAsEAg9PO6o05HY68PFV5fjYDWV441g7trzWhCd/fWy0cmOWz4sVFX6srSrEH946n+u9iGjKmGBR3BRke7Dp3sX4i9qF2NfchcaTF9Fw4iJeP9qGH+8+i5/taca/PbSWgxYRpYsGGMnUWAUAGif7YSKyDsC66urq6caV1kQEt1YX4dbqIvQNBXGgpRt7znRh39ku7DnTiR0HL+B89yC+8pGlVodKREmKCRbFncflQE1lPmoq8/HQbUb9+7MNp7Hp+b34PzuO4C8/sMjqEImIZpyqdopIk4jkjdnvKm8izSxifB5nsOIs2+vCmnkFWDPv0kTj//vTfXjq9eOonpWDDWvis96NiNILN0SiGSci2LCmEutXV+CffvUeXj/aanVIRETxFqsUEAA2A3gscmB2Fpx0cmW+d52I1Hd1dU3l7TRBX/3oUty+sBhf+fE+vHls0pWcRERMsChx/vZjy3HdrBz8+Y/ewfnuQavDISKaNhGpEpHNMBKpGhHZYnb8A2A0pwBwTERqzZbstVPdWFhVt6lqnd/vj0/wFJPL6cA/f2oV5hVl40/+bRdOtPVZHRIRJRkmWJQwmR4nvvmpGvQPh/D5H+5GMBS2OiQiomlR1SZV3aSqq1VVVHWjmVRF31OvqjtUdauqPmFVrDRxuRlufPtzAQiAP/re2+jqH7E6JCJKIkywKKGum+3D//z4crx1vAP/uOOo1eEQESUNlggm1tzCbGz5bACnO/rxp//eiJauAatDIqIkwQSLEu6TqyvwQKAC33z1Pfx49xmEw7G2iiQiomgsEUy8G+cX4Ov3XY/fvNeGmx//Je79x9ew+T8P4a3jHazCIKJxSWQ39HQXCAS0oaHB6jDSxsBwCOv/9Q3sb+7G/KJsfHbtXKwPVCCX+2QRpSUR2aWqAavjSAYcrxLvvQs9eOXgBfzq8AU0nLiIYFjhy3BhQXEOSv0ZKPVnoizP+J2b6UKG2wmvywGvy/hdnp8Jt5PPtImS3UTHKiZYJg5YiTccDOOlfS343hsn0HiqE9keJz65ugKfqKnA4hIfMtxOq0MkogRhgnVtUftgPXz0KEusrdI9OILfHm3D6++14VR7P5q7BtDSOYiBkdC471lamotnH7kZOV7ujkOUzJhgTRITLGvtOdOJ775xAj97twXDoTBcDkH1rBwsLcvF0tJcLJiVg0y3ExluJzLcxlPBgiwP/Fmc8SJKBUywJo7jlf2oKroGRtDcOYi+4SCGRsIYCoYwFAyjuXMAj790CHcumoX6z66GwyFWh0tEUzTRsYqPUsgWVlTk4X8/cAO+/OEl+F1TBw60dGF/czd+c7QNLzSejfkel0PwDw+sxMduKE9wtERERJeICPKyPMjL8sS87nII/se2A/j7XxzGo/cuTnB0RJRoKZ1giUgVgPVsi5s8inK8+MiKUnxkRenoudaeIZzq6MfQiPE0cHAkhMFgCD966zT+/Jl3MBQM44HAHAujJiKaeVElglaHQpP0uVvm4fD5XvzLq8ewcLYPH1/FB4NEqSylEywAmwB0Wh0ETU+xz4tin/eK8/cuK0Xd0w14dOsejITC+PRNcy2IjogoMVR1G4BtgUDgYatjockREXzt95ahqbUXjz6/B3MLs7CqMt/qsIhohqRsSxsRqQWwy+o4aOZkepx46r8EcNfiWfjKj/fhO789bnVIREREMXlcDjz5mdWYnetF3dO7uK8WUQpL2QTL1GF1ADSzMtxO/OtnVuODy2bja9sOYMuvj1kdEhERUUwF2R58+3Nr0D8UxB99twFtvUNWh0REM8DWCZaI1IjIc+NcqxOR9ebPo2Ou1arqjsRESVbzuBz450/VYN3KMjz+0iF89Sf7MBQcv10uERGRVRbO9uHJz6zG8bZerH/yDZxq77c6JCKKM1smWGZitRnABgBVMa7XAYCqblXVrQC2isgW81oVgKZExkvWczsd+McNN6Du9io8vfMk1j/5JgctIkopIrJOROq7urqsDoWm6faFxfi3h9aic2AEn3jyDew7y/+mRKnElgmWqjaq6iYAz4xzy0ZVrY+6vwlArXlYA6DWTMLuAVBjrseiFOd0CL784SWo/+xqnGzvw0f+6XW8vP+c1WEREcWFqm5T1Tq/3291KBQHq+fmY+sjN8PjFDxYvxNvvNdmdUhEFCe2TLCuRkTyYCRRY3WapYFbVbXeTMB2AWhkuWB6+cCyEvz8C7dhflE2Nj69C3/7swPoGRyxOiwiIqLLVM/y4fn/egvK8jLwue+8hRcazyAYClsdFhFNUzK2aa9C7NbrHTASrx2AUWYI4H4ABVyTlX7mFGThuUduxuMvHsK3f3Mc3/7NcSwozsbKijysqPBjxZw8VBZkIcfrQobbaXW4RESUpkr9mXhu4y146Ptv4y+ffRdf/ck+1MzNx5p5BbhxfgFumJPHcYooySRjglWA2N0BOwEURg5UtRFGieC4zDLCOgCorKyMY4hkB16XE//j95bhoytK8eaxdrx7phOvv9eGF3afvew+j9MBX4YLvgwj2XI7HXA6BG6nwOVwwOUUiAgcAjhE4BCBCCCA+ds8Nl+Pkst+XZPLIfC4HPC4HPC6nPC4HFhdmY/apbPj8u+DiIjsyZ/lxg8eugnbD5zHW8c78NbxDvzv7UcAAA4BcjPd8Js/uRnGb6/LAbfTGKPcTgfcToHT4YBDjJL5yHjliIxPYoxGl41b5vfL6Hh15YglMQYxr8uBDLdz9CfL40RNZT4yPUwEiYDkTLDixiwjrAeAQCCgFodDMyQwrwCBeQUAAFXFue5BvHu6Cxd6BtEzGET34Ah6BoPoHQxiYCSEYCiMYFgRDCmC4TAGgwpV470hVYTDQFjV/DxAYV6P+k6NXJ9EnMGQYjgYxnAojKGREIaCYfyrKl76s9uwuCQ3Tv82iIjIjrwuJz66ogwfXVEGAOjsH0bDiYvYc6YTHf3D6B4IomtgBN2DI2juGsBwMIxgSDESCmMkZIwdkfEppMa4lEifXTsXf/vx5Yn9UiKbStYEqyDGuTwA7YkOhJKLiKDUn4lSf6bVoVzTxb5h3P7Er/D3Lx/Btz4XsDocIrKYiKwDsK66utrqUCgB8rI8qF06e8pVDKqKsBoJV/TDQOOacRx5DcR+IKgxsrSwAsPBMAZHQuZPGP/y6nt4ofEMNn1oMXK8yfpXS6L4Scb/FzTASKbGKgDQmOBYiGZMfrYHdbdX4R+2H8Gukxexem6+1SERkYVUdRuAbYFA4GGrYyH7ExE4BXBOuFB96jbesQAv7TuHH+8+i8+unTvj30dkd0nXRVBVOwE0md0Eo+WxkQWlmj9633wU5XjwjZcPxXySSESpQUSqRORRq+MgmoqVFX4sL8/FD948ybGKCPZPsGKVAgLAZgCPRQ7MjoFMrijlZHtd+NM7q7GzqQOvH+UeKUQpbBOiGjURJRMRwWfXzsXh8z1oOHnR6nCILGfLBMt8krcZRiJVIyJbzI5/AEabUxwTkVoRWQ+gVlU3WhUv0Uz61E2VKM/LxDdePswng0QpSERqYezbSJS01q0sgy/DhaffPGl1KESWs+UaLFVtgvE072r31CcoHCJLeV1O/HntdfjS1j14ad85fPj6UqtDIqL46wCwwOogiKYqy+PCJ2sq8G+/O4m23qUoyvFaHRKRZWw5g0VEl/tETQWqZ+Xg739xGMFQ2OpwiGgMEakRkefGuVYnIuvNn0fHXKvl+mFKFZ9ZOxcjIcUzb5+2OhQiSzHBIkoCTofgix9YiKbWPrzQePbabyCihDATq80ANgCoinG9DgBUdauqbgWwVUS2mNeqADQlMl6imVQ9Kwc3VxXi3393CqEwS9opfTHBIkoSH1xWgpUVfvzD9sN49fAFrscisgFVbVTVTQCeGeeWjdEl7WYJfK15WAOg1kzC7oGx5rg2xmcQJY3P3jwXZzsH8OrhC1aHQmQZJlhESUJE8LWPLYcq8AffeRsf+P9eww/fOoXBkZDVoRFRDOZ2IjUxLnWapYFbVbXeTMB2AWhkuSAlu3uWzsYsnxdP72SzC0pftmxyQUSx3TAnD7/ZdBd+tqcZ33r9OB57YS++8fJhPBCYg7K8DLgcDrgcApdT4HQIHCIQAQSR34BctuekcSDj7EN5te0pRQQCwOEAMt0uZHudyPa6kO1xISfDhRwv/3ihtFcFoDPG+Q4YidcOYHSrkfsBFHBNFiU7t9OBB2+sxD/98ihOtfejsjDL6pCIEo5/AyJKMh6XA5+oqcB9q8qxs6kD3/5NE7a8dgx2qxgs9nmxuMSHpaW5WFKai8WlPuRleuAQAAI4xEwAgdEkEGK8znI74XJygp2SXgGMZGqsTkTteaWqjTBKBMdllhHWAUBlZWUcQySKv9+/cQ6++av3sPnlQ/i7T1wPX4bb6pCIEooJFlGSEhHcvKAQNy8oRN9QEAMjIYTCipFQ2PytUFUoAFVAoZclYZHXxh1XulbCFrkeUkX/cBD9QyH0DQfRNxRC9+AIjp7vxaFz3fjOb09geJKdD31eF+5eMgsfur4UdywsRobbOan3E6UaVa0XkRYA6zwez2qr4yG6mlJ/JjbeXoV/efUYdh5rx1/csxAPrpnDB2eUNphgEaWAbK8L2TYtyRsJhdHU2odD57rRNxSCQhE2s75QODoBxGjjjiPne/CLA+fxk3eake1x4q4ls/GR60tx95JZcHOApuRSEONcHoD2yX6Qqm4DsC0QCDw87aiIZtij9y7GvctL8D9/dhB//ZN9+N4bJ/DlDy/B+xcVQ8arSydKEfb8GxkRpQy304FFJT4sKvFN6n3/KxTGzqZ2vLj3HF7efw7b3m1GSW4GPn1TJR68sRLFPm5iSbbXACOZGqsAQONkP0xE1gFYV11dPd24iBJiRUUentm4Fi/vP4+/e+kg/vC7b2NuYRayPC54nAKX0wG3U+ByOIxScRE4RtcLy+g6YCMfG3/NcKx0LXKfRL3P43LAl+FCboYbuZlu+DJcmO3LQM3cfBRke+L8T0/pTNjq2RAIBLShocHqMIgohmAojF8dbsX33zyB14+2we0UfPj6Unz6prm4YU4ePC7OaiU7EdmlqgGr45gqs1HFU6q6esz5YwBWq2pn9DlVXTDV7+J4RcloOBjGv//uJH53vAMjIUUwHMZIKIyRkFHaHqlkgBpVDpHyddXokvYrXevvsdGXB4Mh9AwG0T0wguCYfboWFGdjzbwCrJlXgBUV/tHSdIfDSPourRk2mkY5ROAUQW6mizNyaWSiYxVnsIjI9lxOB+5ZOhv3LJ2NptZePL3zJLY2nMFP32mG2ym4bpYPS8tysawsF4tKfMjxukabaDgdAqOq8MoBcLwnoaMNOMQYXLM9TuRl8ekmXVWsUkAA2AzgMQCbgNFEbEpdAjmDRcnM43LgD26djz+4db7VoUBVMTASQvdAEKc6+tFwsgMNJy7ixb0t+NHbpyf1WcvKcvG5W+bh91aWcb0wjeIMlolPBImSS99QEL86fAH7m7txoLkb+5u70dY7NGPfV+bPwPJyP64v92N5uR8r5+SxpCSOknUGS0SqAGyEsXlwDYB6ALuiNxc2OwA2wSgXrFLVJ6bznRyviGZGOKw4eqEXB1u6EQwrwmo0iworzNe4rHlU71AQP33nLI6c70V+lhsb1lTiM2srUZHP1vSpaqJjFRMsEwcsouR3oWcQR871YnAkhJA5MIbCRqfDscb7sy/ScTEcvjSgdvQPY39zN/af7UJTWx8AIMPtwNZHbsHycv+M/jOli2RNsKzA8YrIPlQVbza14/tvnMQvDpwDAHx0RRm++IFF3AMsBbFEkIjSzixfBmb5Mmb0O3oGR7DvbDf+9N8b8TfbDuCZjWtZf08JwRJBIvsREdyyoAi3LCjC2c4BfP/NE/jeGyfw0r4WfPqmufj8XdUozGFTpnTDleFERJPgy3Dj5gWF+G8fWIi3TnTgxb3nrA6J0oSqblPVOr+fs6ZEdlSel4nHPrQEv/7SnVi/ugLff/ME7vjGq/jnXx5F/3DQ6vAogZhgERFNwYNrKrG4xIevv3gQgyMhq8MhIiKbmJ2bgcc/sQK/+IvbccuCQvz9L47ghr/Zjge2vIl/+MVhvH60FX1DTLhSGddgmVjTTkST9eaxdvz+Uzvxl/csxBfuvs7qcJIa12BdW1SJ4MNHjx61OhwimqDGUxfx0t4WvHW8A/uauxEKK1wOwdzCLLidDrPb7aXOt5fv/2W+HqcTbuR9kdfzCrOxtqoQa6sK2P12BnANFhHRDLt5QSE+fH0Jnnz1GO4PVKDUn2l1SJTCVHUbgG2BQOBhq2MhoomrqcxHTWU+AKPz4K6TF/HW8XY0tfYhZHYrDIZ19DVw+f5dseZCFEYTppFQeLTLYTCk2NnUju++cQIiwJKSXNy8oBAPrpmD62b7EvGPSiYmWERE0/DYh5Zgx8EL2PzSIfzjg6usDoeIiGwsx+vCHQuLccfC4hn5/OFgGO+e6cSbx9rx5rF2PL3zJP5z3zn8+kvvh8vJlUGJwn/TRETTMKcgC3W3VeEn7zRj18kOq8OhFCYi60Skvqury+pQiMimPC4H1swrwBfuvg4/rFuLf/79VTjbOYD/3M+GTInEBIuIaJr+5P0LMDvXi69tO4BwmOtaaWawiyARTdbdS2ZjXmEWnnr9+Lj7P1L8McEiIpqmbK8Lm+5djD1nuvAXz76DXnaHIiIiG3A6BH/8vvl493Qndp28aHU4aYMJFhFRHNy3qhz/7Z6F2PZuM9b902+w7yzLuIiIyHqfXF0Bf6Yb33r9uNWhpA0mWEREcSAi+Pzd1+GHD69F/3AQn/iXN/D9N0+wJIOIiCyV5XHh0zdV4uUD53Cyvc/qcNICEywioji6qaoQL37hNtxaXYj//tP9eOQHu9A9OGJ1WJQC2OSCiKbqc7fMg8sh+M5vT1gdSlpggkVEFGeFOV58+3Nr8JUPL8ErBy/gsRf2Wh0SpQA2uSCiqZqdm4F1K8vwbMNpdPXzod9MY4JFRDQDHA7Bw7dX4Qt3X4ef72nBrw5dsDokIiJKYw+9rwr9wyH8+1unrA4l5THBIiKaQY/csQDVs3Lw1z/Zh/5hdhckIiJrLC3Lxa3VhfjuG8cxHAxbHU5KY4JFRDSDPC4Hvn7f9TjbOYB/3HHU6nCIiCiNPXRbFc53D+HFvS1Wh5LSmGAREc2wG+cX4ME1c/Dt3xzH/mY2KCAiImvccV0xqmfl4Ju/eo8NmGYQEywiogT4qw8tRn6WG19+YS9CYbZup8ljF0Eimi6HQ/CVjyzBifY+fOqpnejoG7Y6pJTEBIuIKAHysjz46keX4t0zXfjBzpNWh0NJiF0EiSge7lw0C/WfDeDI+V48WP8mLvQMWh1SymGCRUSUIL+3sgy3XVeEb7x8GOe6OKARICLrRaRWRB4VkTqr4yGi9HDn4ln47h+swZmLA9iwZSeaOwesDimlMMEiIkoQEcH/+vj1GAmF8aWt7yIYYhendCYieQA2q+oOVX0CwCYRqbI6LiJKD7dUF+HpP74Rbb1DuP9f38TJ9j6rQ0oZTLCIiBKosjALf/OxZXj9aBu+/uIhq8MhC6lqJ4DVY053WBELEaWn1XML8MOH16J/OIgNW3ay8UWcMMEiIkqwDWsq8Ye3zsP//e1xPPM2N3xMZ2aSBXPmqjFyTESUKMvL/fi/f7AG57oH8f03TlgdTkpggkVEZIGvfHgJbruuCH/9k314+wQnLZKdiNSIyHPjXKsz11qtF5FHY1yvBbAZwKaZjpOIKJZVlfm4a/EsfPs3x9E3FLQ6nKTHBIuIyAIupwP//KkazMnPwiNP78KZi/1Wh0RTYCZWmwFsAHDF+qlI4wpV3aqqWwFsFZEt0feYa7DuB7Cda7CIyCr/z13VuNg/gn/7HTvdThcTLCIii/gz3XjqcwEMh8J46HsNfGqYhFS1UVU3AXhmnFs2qmp91P1NAGqB0eRsfdS9jQA2zliwRERXUVOZj/dVF6H+teMYHAlZHU5SS8kEyyzDqBORzWx7S0R2tqA4B9/8VA2OnO/BV3+6z+pwKI7MLoE1MS51mmWBVQAKos7nAXg7EbEREcXy+buq0dY7hB+9xfXB05FyCZZZXtGpqvXmU8WNIhJrgCMisoXbFxbjwRsr8dLecxhh6/ZUUgUgVtOKDgA1ZslgR+SBIIwmF1sTGiERUZSbqgpx47wCbHmtCUNBzmJNVcolWDCeFkaXWDTALMcgIrKr26qLMDASwp4zXVaHQvFTgNht1zsBFAKja7PqVXWT+VAwJjMJaxCRhtbW1hkKl4gI+Pzd1WjpGsQLjWetDiVppVyCZT79ezjqVABGXTsRkW3dON+oFNvZ1G5xJGRHZhIWUNVAcXGx1eEQUQp7X3URVs7Jw7+8+h6rKqbI1gnWVNvejtlXpENVdyQgXCKiKSvM8WLRbB8TrNRTEONcHoBJ/4cWkXUiUt/VxVlOIpo5IoLP31mN0x0D+I93mq0OJynZMsGKR9tb0yZVvWdmoyUiio+1VQVoOHGRTwxTRwOMZGqsAkyhskJVt6lqnd/vn3ZgRERXc/eSWVhSmotv/uo9hMJqdThJx5YJ1nTa3kaYs1qbzNexBjgiIltZW1XIdVgpxKymaIoxBuVNpbKCM1hElCgigi/cVY2mtj78kB0FJ82WCdbVTKDtLczf9ZFSQbDJBRElAa7DSmqxSgEBYDOAxyIHZlfbKZWtcwaLiBLp3uUluLW6EJtfOoRzXYNWh5NUki7BwjXa3pqD13YAF0VERUQBNCUyQCKiqeA6rOQjIlVmSftmGGPQluj9F81qi2MiUmtuKlyrqlPaTJgzWESUSCKCr993PUbCYXz1p/ugylLBiXJZHcAUXLXtrao2ApCJfJA5CNYBQGVlZdwCJCKaqrVVBXi24QxGQmG4ncn4DCy9mCXq47ZXN++pv9r1SXzXNgDbAoHAw9e8mYgoDuYWZuMv71mIr794CP+57xw+dH2p1SElhbRJTw3oAAAgAElEQVQevdn2lojshuuwaDycwSIiK/zRrfOxvDwX//0/9qOrf8TqcJJCsiZYcWt7S0RkJ1yHRePhGiwisoLL6cDffWIFOvqG8fhLB60OJykkY4IV17a3RER2wnVYRERkN8vL/Xjotvn40dun8caxNqvDsb2kS7Di3faWiMhuuB8WxcISQSKy0p/fvRBzC7Pw5Rf2YnAkZHU4tmb3BGvG294SEdkN12FRLCwRJCIrZXqc+Pp91+NEez9+vPus1eHYmi0TrES2vSUishuuwyIiIju6ZUEh8rLcePd0rB2TKMKWbdoT2faWiMhuCnO8WDg7Bzub2vGnd1ZbHQ7ZhIisA7Cuupr/myAia4gIlpf5sa+ZFRZXY8sZLCKidLe2qpDrsOgyLBEkIjtYXu7H4XM9GA5yfBoPEywiIhviOiwiIrKj5eW5GAkpjpzvsToU22KCRURkQ1yHRUREdrS8zJhF33eWDwDHwwSLiMiGiqLWYREREdnF3MIs+DJcXId1FUywiIhs6uaqQrx9ogNdAyNWh0I2wH2wiMgORATLynKx92y31aHYFhMsIiKbuj8wB4MjYfzwrVNWh0I2wCYXRGQXy8v8ONjSzUZM42CCRURkU8vL/bhlQSG++9sT7NZERES2sbzcj+FgGMdae60OxZaYYBER2djDt1XhXPcgfr632epQkoqI3CAiXxSRu6LO+UXkkyJyg5WxERElu+Xlxkz6Xna6jYkJFhGRjd2xsBjVs3Lw1GvHoapWh5MURORhAN8CUA3gCRF5GQBUtUtVnwewy8r4iIiS3fyibGR5nNjfzHVYsTDBIiKyMYdD8ND75uNASzfeOMaOghO0QFUDqvqIqgYA/ImIPCkiueZ1sTI4IqJk53QIlpbmslX7OJhgERHZ3MdXlaMox4OnXm+yOpRkcSz6QFWbVPVPANSJiB9AUk4FsosgEdnJ8nI/9jd3IxROyj9SZxQTLCIim8twO/Ffbp6HVw+34uj5HqvDSQYdInK3iDwTNWsFVf17AAEk6QwWuwgSkZ0sL/djYCSE421sdDEWEywioiTwmbVzkeF24FuvH7c6FNsz11k1AXhWVbvHXHsFwGpLAotBRNaLSJ2IbBaROqvjISKaqOXlxvOrfdwP6wpMsIiIkkBBtgefrKnAj3efRWvPkNXh2J6qHjcTrVjXdic6nlhEpApAp6rWq+omABtFpMbquIiIJqK6OAdelwN7uQ7rCkywiIiSxB+/bz5GwmE8/eYJq0Oh+KgBsDHquAFArUWxEBFNisvpwBI2uojJZXUAREQ0MVXFObh78Ww8vfMk/uud1chwO60OyTIi8iyAqS5GEhiNLp5T1W/FL6rJUdWtIrIj6lQAwHNWxUNENFnLy3Px093NCIcVDkdSLm+dEUywiIiSyINr5mDHwfN453Qn1lYVWh2OZVT1AatjiGaW9j2mqvfHuFYHoMM8rFLVJyLXVLXTvKcKQIeq7hj7fiIiu7q+3I8f7DyFkx39mF+UbXU4tjGlEsHorkxERJQ4q+fmAwAaT120OBICjMRKRDYD2ACgKsb1OsCYrVLVrQC2isiWGB+1SVXvmdloiYjia1mZUUjAMsHLTXUN1hMi8onIgdkO94Y4xUREROPIz/agqigbjSc7rQ6FAKhqo9mg4plxbtmoqvVR9zdhzDorEXkUwCbzdd5MxUpEFG8LZ/vgcTqwr5kJVrSplgg+Z7a6RVRiFQDwTlyiIiKica2qzMerhy9AVSGSnjXvybAGy0yWYnUF7BSRWlXdISK1AOojpYIwkq+tMxUTEVE8eVwOLCrxYT9btV9mwgmWiHwRwA5VfUdVXxGRu2DsM/JlAG8BeGWGYiQioig1c/PwfOMZnOrox9zC9Kx5t9sarHFUAYg11dgBoEZEOgBsBxCdKNtmjy4ioolYXp6LF/eeS+uHfmNNZgbrAwCqRSQAoB3ADhg14x+ckciIiCimmspL67DSNcFKEgW41NwiWieAQlVthDGbdk3mWq46AKisrIxbgERE07WszI8fvnUaZy4OYE5BltXh2MJk1mDdr6qPqGoAwAMwZq+Oi8jbIvKyOcNFREQzbOFsH3K8Lq7DSiPmZsQBVQ0UFxdbHQ4R0ajry41q7f1chzXqigRLRL4kIg+NPa+qXdGvVfV5M+FaAyPhYokgEVECOB2ClXP87CQIQETmmb9zbdrhtiDGuTwYlSCTIiLrRKS+q4t/iSEi+1hU4oPLIdh9mg/9IsabwRrtYjSR7oBmwrU7blEREdFV1VTm49C5HvQPB60OxWo15prgx6wOJIYGRI2nUQoANE72w1R1m6rW+f1T7e1BRBR/GW4n1lYV4sW9LVBVq8OxhVgJViGM9VURGxIUCxERTVBNZT5CYcW7p9N3NkNEPgljDdNG89QDdtoyxOwM2BSj9XreVDYU5gwWEdnVfavKcbpjAA0nWVkBxE6wtsDYCDEkIi/DeDp4p01LL4iI0tKqSuPv7GleJtho/hzHpeZLxy2KJVYpIABsRtTsmojU4PKHmBPGGSwisqt7l5cg0+3EC41nrQ7FFq5IsFT1uKpWA6gG8DyABQCegrFvR7vZ0OJJEXnITk8KiYjSSV6WB1XF2didxgmWqh4HsB7GJr9F5rmETu+ISJWIbIaRSNWIyBaz418kxnoAx0SkVkTWA6hV1Y3jfd41voszWERkS9leFz64bDZ+vqcZgyMhq8Ox3Lht2s2Bq15EqlT1rwBARFbB2FB4NYzGFk+IiAJ4FsAWVeVGw0RECVJTmY9fHkrvDYdV9RvmS0vWAatqE4BN17inPk7ftQ3AtkAg8HA8Po+IKJ7uq6nAT95pxq8OXcCHri+1OhxLTaRN++ORF6q6W1WfMrsHfkBVC2BspLgDwCOxug8SEdHMqKnMR0ffME6291sdim2kcjk7Z7CIyM5uXVCIYp8XL+xmmeA1E6xrlVtEWrbD3I2eiIgSo2Yu12FFE5FnMWYNlojcnSpJF9dgEZGduZwOfPyGMrx6+AI6+oatDsdSk9loeFwi4gfwHIz1WkRElADXzTI3HGaCFbEdRhn7KFV9BcA9XDNMRDTz7ltVgZGQ4md7mq0OxVJxSbDMWSyHqtpxHxIiopTkdAhumJOHxpPc3NHUAeCK6R2zyqIq8eHEF0sEicjulpblYnGJL+27CcYlwSIiImusqszDoXPd6BtK+w2HAaAJwC9F5KiIPC4id0WVByZ9gsUSQSJKBvetKsc7pzvR1NprdSiWYYJFRJTEairzEVbg3TOcxQJQB+B+AE8AyAdQD3OLEQD8F0RElAAfu6EcIsBP0rjZBRMsIqIkFtlwePcp5g8AGlX1lahut9UwNgD+OwANFsdGRJQWSvwZuHVBEX78zlmoqtXhWIIJFhFREotsONx4ko0uAHSIyLzoE6raae6VlfQlglyDRUTJ4r5V5TjdMYCGNB2bmGARESW5msp87D7dmbZPCiPMZhYLROSu6PMi0gHgHmuiih+uwSKiZHHv8hJkup34cZqWCU4qwUqVvUSIiFJJZMPhE9xwGGaJ4C/HnL4fxnosIiJKgGyvCzcvKETDiQ6rQ7GEa5L3PwfggzMRSDyJSC2AvMixqm61MBwiohk1uuHwyYuYX5RtcTT2Y+6FRURECbS0NBe/PtKKwZEQMtxOq8NJqMmWCMqMRBF/G1V1q5lYbRCRvGu+g4goSUU2HN59Oj1r3dMF12ARUTJZXOpDKKx470L6tWufbIJl+wJ/c/Yqej7ybQAPWBQOEdGMi2w4zE6CqY1rsIgomSwpNVYWHWzptjiSxEvFJhdVuHy/k04ACyyKhYgoIYwNh3vQP8wNh4mIyHrzCrOR4XbgYEuP1aEknK0TLBGpEZHnxrlWJyLrzZ9Hoy7lAWgfcztLBIkopdVU5iMUVrx7muVjRERkPadDsGi2D4fOcQbLFszEajOADYixd4mI1AFG8wpzndVWEdliXu4EUDjmLaybIaKUdsMcc8NhrsMCwK63RER2sKQ0FwdbutNuGxFbJliq2qiqmwA8M84tG1W1Pur+JgC15mETLp+xyoOxDouIKGXlZ3tQVZSNxpN8nmSKWf1ARESJs6Q0Fxf7R3C+e8jqUBLKlgnW1ZgdAWtiXOoUkVpV3QGgIOr8AgA7EhIcEZGFVlXm453TF9PuSeE4kqXr7YSxiyARJZvFJT4AwME0KxNMugQLVzaxiOjApcTrcXNtVi2A7arKR7pElPJWVeahrXcYpzsGrA7FDlIuy2QXQSJKNovTtJPgZDcatsMTwQJc3oY9YnTtlao2Ami81geZa7nqAKCysjKOIRIRJV5NZT4AYx1WZWGWxdHQRIlIFYzS901Wx0JEFE/+TDfK8zLTrpPgZGewUqqmXVXrVTWgqoHi4mKrwyEimpaFs3OQ5XGi8SQbXSQLEVkPYBPY7ZaIUlSk0UU6mVSCpapPzVQgk1QQ41ys9uxERGnD5XRgZUUeGrnhcNIwO+Gm1MNLIqJoS0p9aGrtxeBIyOpQEiYZ12A1IPaTvgJMoCyQiCiV1czNw8GWbgwMp89ARkRE9rWkNBdhBY6e77U6lIRJugTLbFjRZHYTjJZndhAkIkpbq+bkIxhW7D3LTnOJZO7fGHMmSkTqzMZL60Xk0UTHRkRkpSVp2OjC7glWrFJAANgM4LHIgYjUgK3YiYiwqtLccPgU12ElgplYbQawAUaX27HX6wCjFNAsB9wqIlsSHCYRkWUqC7KQ6XamVav2yXYRTIhIRyUYmwfXmIPRrsjmwqpabz4RrIVRLlilqhuti5iIyB4Kc7yYW5iFRiZYCel6G+laaz7oq41xy0ZVXR11f5M5dhERpQWnQ7CoxJdWM1gxEywRmaeqJ0QkFwBUNaH/RlS1CUZXpavdU5+gcIiIkkpNZT5+814bVBUidthdwxKWN44wS9lrYlzqFJFalrUTUbpYUpqLF/e2pM24NF6JYI2I3IWoMjwiIkoONZV5aO0ZwtnO9N1w2CZdb6tg7NE4VgfMxMuczdoIoDZSTkhElGqWlPrQNTCCc92DVoeSEFckWCLySRilFZGSuwdE5IaERkVERFO2ytxwmO3aLVcAI5kaqxNAIQCo6g5VvV9VF1ytMsMsi28QkYbW1tYZCpeIaGakW6OLWDNYjebPcRj7Su0wXxMRURJYXOJDhtvBRhcpRFXrVTWgqoHi4mKrwyEimpTFJT4AwMGWHosjSYwrEixVPQ5gPYBnABSZ59jvl4goSbicDqzghsN2Easbbh6MB5iTIiLrRKS+q4tDMhElF1+GG3MKMtN6Bguq+g1V3a2qf6WqJxIcExERTVNNZT4ONHdhcIQbDluoAUYyNVYBjEqRSVHVbapa5/f7px0YEVGiLS7JTe8Ei4iIkltNZR5GQor9zek92xHphmsFVe0E0GR2E4yWN5UOgpzBIqJktqQ0F8fb+tLiwV+sJhdfEpGHrAiGiIjiI9LoYjfLBJ8QkU9EDkTk7hlq3BSrFBAANiOqI6+5X9aU2rNzBouIktnSUh/CChw5n/rrsMabwRp92sYOgkREyafY50V5XibePZP2sx3PqeoLwGXjWSBeHy4iVSKyGUYiVSMiW6LbrZudAY+JSK2IrAdQq6obx/u8a3wXZ7CIKGktLkmfToKxNhouBPCjqOMNAN5JTDhERBQvKyr82HMm/WawROSLAHao6juq+oq5r2MTgC8DeAvAK/H6LlVtArDpGveM2359kt+1DcC2QCDwcDw+j4gokSoLspDtcaZFJ8FYCdYWANtFZD7MMgYRuRPALlVN/ZSTiChFrKjIw0v7zqGzfxh5WR6rw0mkDwCoFpEALm03sklVP2htWNMjIusArKuurrY6FCKiSXM4BItKfGkxgxWzTbuqVgOoBvA8gAUAngLQKSLtIvKyiDwpIg+xfJCIyL5WVhhrdfakX5ng/ar6iKoGADwAY/bquIi8bY5hX7Q4vinhGiwiSnZLSnNxoKUbqmp1KDNq3C6CZqJVD2CrqlarqgNALYCtAATGoPVLM+l6kskWEZG9LCuPJFjpVSYYvXejqnap6vNmwrUGxtgVtxJBIiKauOXlfvQMBnG6Y8DqUGZUrBLBsR6PvFDV3QB2w5jRAgCIiB9G4vWIiDSo6rfiHiUREU2aP9ONqqLstGp0ISK5VytnN5Ov3QkMKW5YIkhEyW55mfHgb19zFyoLsyyOZuZccx+s6CeB411X1ecBbI9bVEREFBfp1OhCRJ4FcHzMubut3AsrnlgiSETJbmFJDlwOwb6zqf3gLy4bDZuzWM/BWK9FREQ2saIiD+e7h3C+e9DqUBJhO8a0YFfVVwDcwzJ2IiLreV1OLJztw/7m1G50EZcEy5zFcqjqY9e+m4iIEmXlHGO2493TaTGL1QHgiukds8qiKvHhEBHRWMvKcrHvbFdKN7qIS4JFRET2tLTUD6dD0qWTYBOM5ktHReRxEbkrqjww6RMsbjRMRKlgebkf7X3DON89ZHUoM4YJFhFRCsv0GOUYe1K83t1UB+B+AE8AyAdQD3OLEQBJP4XHNVhElAqWlxvPvVJ5HRYTLCKiFLfSbHSRyuUYpkZVfUVVnzLbslcDKADwdwAaLI6NiIhg7IUlYnQSTFVMsIiIUtyKijx09o+k/L4jADpEZF70CVXtVNVvIAVKBImIUkGWx4UFxTnYdzZ1G10wwSIiSnErKsxGFynert1sZrFARO6KPi8iHQDusSaq+OEaLCJKFcvLcrGfM1hERJSsFpX44HU50mI/LLNE8JdjTt8PYz1WUuMaLCJKFcvL/WjpGkR7b2o2umCCRUSU4txOB5aW5eLd9OgkeAUz6dptdRxERGRYVmY8KErV/bCYYBERpYGVFXnYd7YLoXDKN7ogIiKbW1pmdhJM0TJBl9UBEBHRzLu+3I/vvnECx1p7sXC2z+pwpk1EnkWMTYUn+nYACuA5Vf1W/KIiIqKJ8Ge6UVmQhf0p2uiCCRYRURpYOcdsdHG6MyUSLFV9wOoYiIho6paX56bsDBZLBImI0kBVUQ5yvC7sSdN1WEREZC/Lyvw42d6ProERq0OJOyZYRERpwOEQLC/PTYtOgqmKbdqJKJUsLzcqKw6kYKMLlggSEaWJlRV5+M5vT2A4GIbHldzP19JxDZaqbgOwLRAIPGx1LERE07XMbHSxv7kLNy8otDia+GKCRUSUJlZU5GE4FMahc91YUZFndTjTkiprsESkFsDofwxV3WphOERECVOU40WpPyMlW7Un9yNMIiKasBUVZqMLrsOyk42qutVMrDaISHJnvkREk7CszI99Z1NvTGKCRUSUJiryM1GU48HOY+1Wh0IYnb3qiDr1NoCUmJkjIpqIZWW5ONbai/7hoNWhxBUTLCKiNCEiWLeyDL84cA6tPUNWh0NAFYDoriOdABZYFAsRUcItL/cjrMDBlh6rQ4krJlhERGnk0zfNxUhI8WzDaatDSSkiUiMiz41zrU5E1ps/j0ZdygMwdjqRJYJElDaWl19qdJFKmGAREaWR6lk5uLmqED986xRCYbU6nKRnJlabAWyAMSM19nodYDSvMNdZbRWRLeblTgBjW2exjz4RpY2S3AwUZntSbh0WEywiojTzmbVzcebiAF470mp1KElPVRtVdROAZ8a5ZaOq1kfd3wSg1jxswuUzVnkw1mEREaUFEcGycn/KdRJkgkVElGY+sGw2in1e/GDnSatDSWlmR8CaGJc6RaRWVXcAKIg6vwDAjoQER0RkE8vKcnHkfA+Gg2GrQ4kbJlhERGnG7XRgQ2AOfnn4As5c7Lc6nFQ2tolFRAcuJV6Pm2uzagFsV1WWCBJRWpmTn4WRkKK9L3WaLzHBIiJKQ79/UyUEwA/fOmV1KKmsAJe3YY8YXXtllhhuVdUdV9tk2GyU0SAiDa2tLO0kotRRlOMBALT1DFscSfykZIJlPg2sE5HNkQXGRER0SXleJu5aPAvPvH06pcoyUpW5jutrABo9Ho/V4RARxU2RzwsAaOMMln2JSBWATlWtNxcebxSRWDXwRERp7dNr56Ktdxgv7z9ndSiprCDGuVjt2a9JVbepap3f759+VERENlGUbSZYKbQ/Y8olWDDq2jdGHTfgUscmIiIy3XFdMeYUZLLZxcxpQOx9rQoANE72w0RknYjUd3WlVjtjIkpvRT6zRLCXJYK2ZdawPxx1KoApDGRERKnO4RB86sa5+N3xDrx3ocfqcFKO2bCiyewmGC3P7CA42c/jDBYRpZwsjwtZHifaezmDNS1T3PF+wiJdmMxywY6pDGREROnggUAFPE4HfrCTzS6mKVYpIABsBvBY5MAsWZ/SmMQZLCJKVYU5HrSlUILlSuSXmQPLBvPwqjvem8dVIrJFVTeax4/iyl3vI7bHSKQ2qeo9cQmeiCgFFeZ4cfvCIrx2lJ3ppsJ8kLcRRil6jYhsAbArsrmwqtabDw5rYZQLVkXGtMlS1W0AtgUCgYeveTMRURIpyvGmVIlgQhMsVW0E0GgmWrHWRW1U1dVR9zeZg1Lk+ImJfpeZjG0yX+dxbxEiotiWlvnxy0MXMDgSQobbaXU4SUVVm2CONVe5pz5B4RARJaWiHC9Od6TOvoy2WYN1rR3vJ/lZtQDqo5IqNrkgIhrH4hIfwgocPd9rdSh0FSwRJKJUVZTjSakZLNskWJjYjvfXZM6ObQdwUURURBRAU3xCJCJKPYtLfACAQ+e6LY6EroZNLogoVRXleNHRN4RQWK0OJS7slGBdc8f7iVDVRlWVMT8xuwiadfENItLQ2sr1B0SUnuYWZiPD7cDhc+wkaGecwSKiVFWU40VYgYv9qTGLZacEK+HMzYgDqhooLi62OhwiIks4HYLrZvlw+DwTLDvjDBYRparCHGMvrPYUKRO0W4IVtx3viYho4haV+HCwhQkWERElXlGOFwBSplW7nRKsuO54T0REE7e4xIe23qGU2uiRiIiSAxOsGRLvHe+JiGjiFpfkAgDXYdkY12ARUaoqHk2wWCI4HTO+4z0REU3cotFOgkyw7IprsIgoVeVmuuB2SsrMYCV0o+FE7nhPREQTV+zzojDbw1btRESUcCKCwmwv2nqYYE0ad7wnIrKvRSU+lggSEZElinwetPexRJCIiFLI4pJcHDnfi3CKbPSYargGi4hSWWG2N2VKBJlgERERAKOT4MBICKc6+q0OhWLgGiwiSmVFOalTIsgEi4iIAEQ3uuA6LCIiSqwinwdtfcNQTf4qCiZYREQEAFg42wcRdhIkIqLEK8r2YjgYRs9Q0OpQpo0JFhERAQAyPU7MLchiowsiIkq4Ip8HAFKiTJAJFhERjVpckssZLCIiSrgic7PhVOgkyASLiIhGLSrx4UR7HwaGQ1aHQkREaaQw20iwOINFREQpZXGJD6rA0QucxbIbtmknolQ2WiKYAq3amWAREdGoS50EmWDZDdu0E1EqK8jyQARo7WWJIBERpZC5hdnIcDtwqIUJVqKISJWIbLY6DiIiK7mcDuRnedDOGSwiIkolTodg4WwfDp/nXliJICLrAWwCkGd1LEREVivK8YxbIjg4EsKf/Wg3mlp7ExzV5DHBIiKiyyya7WOr9gRR1a0AnrM6DiIiOyjK8aJtnBLBd0534qfvNOOZt08nOKrJY4JFRESXWVTiQ1vvMFpToJMTERElj8Ic77glggdbjMqKXx9pTWRIU8IEi4iILrOkNBcAOIsFQERqRCTmDJOI1InIevPn0UTHRkSUaowSwdgzWJEE69C5HrR0DSQyrElzWR0AERHZy6VOgt1433VFFkdjDRGpAbDBPKyKcb0OGC3xizSq2KKqG83jRwEUjvPx21V1R/yjJiJKbkU5XvQOBTE4EkKG23nZtYMtPSjzZ6C5axCvHWnFhjWVFkV5bUywiIjoMkU5XhTleNJ6BktVGwE0molWbYxbNqrq6qj7m0SkNur4iQSESUSUUopyLu2FVZGfNXo+GArj8PkefO7mudj2bgt+bfMEiyWCRER0hUUlPhw+n74J1tWISB6AmhiXOqOTLCIimpyiHC8AXFEm2NTWh+FgGEtKc/H+RcV4/WgbgqGwFSFOCBMsIiK6wnWzfHjvQi9U1epQ7KgKQGeM8x2InXiNy0zINgKojZQdEhGlq9EEa0yTpcj6q6VlubhjYTF6BoPYfTrWH8P2wBJBIiK6QnleJvqHQ+geDMKf6bY6HLspgJFMjdWJ8dddxWSuxbrmeiwz+aoDgMpK+5bFEBFNR6FZItjed3mCdaClGx6nAwuKc1Dqz4TTIfj14VasmVdgRZjXxBksIiK6Qok/AwBwrmvQ4kgIAFS1HsDXADR6PB6rwyEimhHjlQgebOlB9awcuJ0O+DPdWF2Zj1ePXLAixAlhgkVERFcoNRMsu7fCtVCsx6Z5ANpn6gtVdZuq1vn9/pn6CiIiS2W4nfB5XVfsw3iwpXt0CxEAuGNRMfad7bbtfo1MsIiI6AqcwbqqBhjJ1FgFABpn6ktFZJ2I1Hd1dc3UVxARWa7I50Vb1GbDbb1DaO0ZwpJS3+i5OxYWAwBes+mmw0ywiIjoCrN8GRABWphgXUFVOwE0md0Eo+XN5P5WnMEionRQmO1Be1SJ4GiDi6gZrKWluSjK8eLXTLCIiChZeFwOFOV4cb477ROs8VZQbwbwWOTA3C9rRjcP5gwWEaWDopzLZ7AONBsJVnSJoMMhuH1hEV4/2opQ2H7dbplgERFRTKX+jLSdwRKRKhHZDCORqhGRLdFt1M2mE8dEpFZE1gOoVdWNMxkTZ7CIKB0U+TyXJVgHW7pR6s9AfvblDX7uWFiMi/0j2HPGfu3a2aadiIhimp2bgVPt/VaHYQlVbQKw6Rr31CcoHCKitFGY7cXF/hEEQ2G4nA4cbOm5bPYq4vbriiEC/PpIK1ZV5lsQ6fg4g0VERDEZM1jsImgXLBEkonRQ5DNatXf0DWMoGMKx1t7LGlxE5Gd7sLIiz5brsJhgERFRTKPuNGcAABLBSURBVCX+DHQPBtE3FLQ6FAJLBIkoPRSbmw239g7h6PleBMMacwYLMMoE3zndiYt9wzGvW4UJFhERxRTZC+scG10QEVGCFJqbDbf3Do92EBwvwXr/omKoAq+/15aw+CaCCRYREcVUkpsJgHth2QVLBIkoHRSZCVZb7xAOtvQgw+3AvMLsmPeuqMhDtseJXSc6EhniNTHBIiKimCIzWOnaSdBuWCJIROmgyCwRbOsdwoGWLiwqyYXTITHvdToE84uzcdxmDZmYYBERUUwlkRJBNrogIqIEyfG64HE50NY7jIMtPZdtMBzLvMJsnGjrS1B0E8MEi4iIYspwO5Gf5eYMlk2wRJCI0oGIoDjHi31nu9A1MIKlMToIRqsqysaZi/0YDoYTFOG1McEiIqJxlfgzuQbLJlgiSETpoijHg4YTFwGM3+AiYl5RNsIKnOqwT5kgEywiIhqXsRcWEywiIkqcwhwvhkPGjNTiCSRYAGxVJsgEi4iIxlXiz2CbdiIiSqhIo4vKgizkeF1XvbcqkmC1M8EiIqIkUJqbgY6+YQyOhKwOhYiI0kSkVfuSa6y/AoC8LA/ystxo4gwWERElg0gnwfOcxbIcm1wQUbooHE2wrl4eGGG3ToJMsIiIaFylfm42bBdsckFE6SJSInitFu0RVUVMsBJGRKpE5FGr4yAiSlaje2FxBouIiBLk5gWF+NgNZbh5QeGE7p9XlI3mrkHblLOndIIFYBOAif2XISKiK0QSLHYSJCKiRJnly8D/eXAVfBnuCd0/z2aNLlI2wRKRWgC7rI6DiCiZ5Xhd8GW4WCJIRES2Nb/QXq3aUzbBMnVYHQARUbIz9sIasDoMIiKimOYVZQEAjrfZY7NhSxIsEakRkefGuVYnIuvNnymtnxKRWlXdMb0oiYgIAGbnZnAGywbYRZCIKDZfhhtFOV4cb+u1OhQAwNV37oozEakBsME8rIpxvQ4AVHWreVwlIltUdaN5/CjGX1O1XVV3iEgVgKa4B09ElKZK/Rk4fK7H6jDSnqpuA/D/t3c/y21c2R3Hf0eSFf2piBAtJWI50VjQNptQeINAy+yozBOYfAOx/AQuzhuQ8wamdl5SrnmAkNxnQWgmlapJjS0KjhVVZI11sri3yVaz8YdgE9238f1UdZEAmt33oFE4vLfvn+96vd5XdZcFAJrm0b1b+mND7mDNtYLl7oeSDmNFq1+yy4a7P8ntP4hjqbLHv5viNKuSls1Mkp5I6nJHCwBm92Dppn54+14ffv2oz662vWf5fJnZmqRlSY8lHbn7Ts1FAoAkPbp3W3/4jx/qLoakOVewxjGzjkLlqGh4ngpSdvcrHnNd0pDKFQDMbmXphtylv/z8Xl90btZdnNaIPS6GuV4bB2a2HxsjAQDn8OW92/ph/7/08/99mHr2wcvSpKbIrqRhyfPHKq94jRXvkj2T1M/fBQMAnM/JWlhMdFG1VUkbucf7Ku/dAQCYIJtJ8E+v6+8m2Jg7WApdJMpm/RtqhrWsYgvg03H7xDtc65L08OHD854CABbCCmthXQp3f2Fm+R4WPUmlE0ABAMZ7dD9UsF79+L/6py+Wai1Lk+5gzZ2777h7z9179+/fr7s4ANBIK3dCt8BFnEnwsme9dfdhPFZX0jFd2gFgNr9ZPq1g1a1Jd7CkcBerqCPp9bwLAgAI7ty8ppufXV2oO1jzmPW28Nymu4/tdQEAGO3m9ataWbrRiMWGm1TB2leoTBUtS2LALwDUxMy0srRYa2HNadZbSSeVsc34eye7qwUAOJ9H927r1ev6K1iN6SIYE8ogziaY16HLBADU68HSDf2ZSS4kTZ719pzH6kvayVWqmOQCAGb05b3bjegiWFcFq6wroCRtSfo6exBbDqlcAUDNHizYHawJKpn1Nua4PUlvzMzNzCUNqikiACyeR5/f1vDdBw3f/VJrOebaRTAO4t1QaKFbNbNtSQfZworuvhMHDfcVugt2s/7sAID6rCzd0F9+fq9fP7quXrG6i1O3Sma9jd0Qp3ozmfUWACZ7dO90oot/fni9tnLMtYLl7gPFfuZj9mEVewBomAdLN/XXj67Xb9/r7+7cqLs4Cyfmxh1J6vV6XnNxAKCRvvykgnW3tnI0ZgwWAKC5Vu6wFlbB3Ge9NbN/NbOdn3766bJOAQBJe7h8S1dMtc8kSAULADDRAxYbzmPWWwBooOvXrugf7t7Sq9fvai0HFSwAwEQrsYL138wkWNust+7+nbuvLy0tXdYpACB5YSbBt7WWgQoWAGCi5dvXdf3qFf35fxbuDhaz3gJAQh59fkt//PGd3OsbrkoFCwAwkZkt1FTtZtY1sy2FitSqmW3HmfwknUw6cWRmfTNbk9S/7FlvGYMFAJN9ee+23r7/q358W99U7XOdRRAAkK6w2PBiVLCaOOutu38n6bter/fVPM8LACnJT9V+/2//ppYycAcLADCVB3cW5w5WE3EHCwAmyypYdc4kSAULADCVldhFsM5+7YuMSS4AYLIvOjd17Yrp1WsqWACAhvv7Ozf0y68fNXz3oe6iAABQ6trVK/rH5Vv6zxqnaqeCBQCYyvVrIWV8+Pix5pIsJroIAsB0bnx2Vb/8Wl+uooIFAEAC6CIIAGmgggUAAAAAFaGCBQBAAugiCABpoIIFAEAC6CIIAGmgggUAAAAAFaGCBQAAAAAVoYIFAAAAABWhggUAQAKY5AIA0kAFCwCABDDJBQCkgQoWAAAAAFSEChYAAAAAVIQKFgAAAIBWca/v3OZ1nr1BzOwHSX+6wCHuSfqxouI0AfE0V5tikYinyeYZy2/c/f6czpU08tUZbYqnTbFIxNNkbYpFml88U+UqKlgVMbN9d+/VXY6qEE9ztSkWiXiarE2x4FTbrmub4mlTLBLxNFmbYpGaFw9dBAEAAACgIlSwAAAAAKAiVLCqs1N3ASpGPM3Vplgk4mmyNsWCU227rm2Kp02xSMTTZG2KRWpYPIzBAgAAAICKcAcLAAAAACpCBQsAzsHMumbWr7scAACMQ76qz7W6C5A6M1uXdBwfdt39d3WW57zMbFXS1+7+rOS15GKLZZakJ/HnprsPC68nEVMsayc+fCxpy90HhdeTiKXIzHaLn7mE4lmV9Hsz60gaStpX+JwdZjskFMsJM3uuEM+xJLn7i9xrycWDs1K+juSqZsfU1nyVeK6SyFe1xUMF6wKyL8jswsaWgm1336i3ZJPFZPXb+LBb8npysZnZurvv5B9LOlD4sk8qJjN7nv9SMLM1SXtKMJai+NlbKzyXVDzuftfMOvl/iDKpxSKFfyIUku4gPnYzu+vuwxTjwVmpXkdyVfNjamu+akOukshXtcXj7mwzbpIOSp47qrtc54xhdUQcScWm0HK2XvL8G0n91GKSdCRpLfe4K8kldVKLpaSc/fDVk+bnLX9dRryeTCyxbOuSnhee66YaD9vI65z0dSRXNTqmVuar1HNVLBv5qqZ4GIM1o3i7dbXkpWHq/V0Tja0raTuWPW8gqZtgTE89d8tbIb6hhxaa1GI5YWZr7v6y8Fyy8RQlGsuWpPxnTX7aMphiPCho63VMNK625Sqphfmq7blKSjaeZPIVFazZdRX6fxYdq/wCpyS52Dz0J37iZ2+BdxUTlxKKyXN916NNSVk/8KRiycTuFoclLyUXj5n1c9vz3D9LScUSy92Jv6+lHg9Gaut1TC6utuUqqX35qk25SiJfqaZ4GIM1u2WdDqLLG0r6fM5lqVqSsXlu0KZ00g984O4vY+tFcjHFGJ4qDBjOWtOSvD4Kt/FflDyfWjyH0ietZgNJuwrXKbVYsoTU8dM+6/uSvlcYfJ9aPCjX1uuYZFxtzFVSq/JVW3KVRL6SaoqHO1hopdii8bWkf6m7LBfh7i88DM5cNbPtusszq9jdoixhJcfdB/kW2/h7N7Z6pmZZoUUwH89QCq2edRUKWBRtyVVSO/JVm3KVRL6qExWsi1kuea4j6fW8C3IJUo9tS9KzQjeMZGPyMEPTv+W+RJKJxcyyri/jJBPPCENJvfh7SrEMpNMklZPvUpFSPBitrdcx9bhalaukdPPVguQqiXw1F3QRnN2+Ttd8yFtWed/dlCQdm4X1ET5Zg0MJxRRblr5397uFlwYKt/W/USKxRH1JnWILk52uY/GtEoknJuAjd7fCS8dxS+ZzJoXWTLNiKCeyNVOSiQcjtfU6Jh1X6rlKal2+ak2ukshXObXEQwVrRnF2nEHJ2gKd4swzqUk5trgGwgv/dIHDfuzbnkpMy5J2Sp7PviyTuj6eW+8lY2Zb/um6KanEcyypbD2NnqSN1K5NdGhm3cI/eV1J+4nGg4K2XseU42pJrpJalK9alqsk8lWmlnjoIngxWwp9pyWdtOQ09UM5StntVCnB2GKr035+ys5CS1QSMZV9EeT6S38bfyYRyzkkEU9J14TsH6Vvc1/4ScSSsxk3SSflHeQG4qcWD8qlfh3JVQ20gPkqmVjIV/XGYx4W4cKM4od1oHBbsptv6WiyeOt4Q+GW+KpCC9SBn11dPonYslvhI16+mxsImURMceDzeu6pxyp0JUkllrz4T8SGpDWFtSy2swSdUjy5LiMd6WTMQf71ZGKRTmb/6saHn7v7ZuH1pOJBuRSvI7kqiZhal6/akqsk8lVd8VDBAgAAAICK0EUQAAAAACpCBQsAAAAAKkIFCwAAAAAqQgULAAAAACpCBQsAAAAAKkIFC41hZutm5nHrFF7rxOf7o/6+wnLsmtnuZZ/nPMysa2Z78T04qLs8qTKz7brLACB95KvRyFfVIF+ljQoWmmioT9fUgLQn6VBhfZGvRu1kZgdmtje3Uk2hCWUys35M9L06ywGgdchXZ5GvLlYG8lULXKu7AECJHYWVuBu92N28xJXIu8XF9EZoYotX7WVy95exNfBx3WUB0Crkqxzy1cWRr9qBO1hoom1JnbhaN05XLJ/I3XfcfSd7HFvC5rKa+KhzFctUo6eS/r3uQgBoFfLVp8hX1SBfJY4KFhrH3QeSXiq0CgJV6St8rgCgEuQrXBLyVeKoYKGpNiWtxu4GpeLg3r3Cc55vScwGAJvZczM7MrM38fdObhDuUcl5Oma2lfubM90G4utv4rZVOOd21o/azJ5PCrZQvt1s0HQ87m4utqMJxzkZ8Bzfm73c3xbfm3OVP/fYs/cx9zcjz1U2CHtUvPn9i9ds0ns44X1ZlXTs7sP4eC0OUq+9OwiA5JGvRL4iXyGPChYayd0PFQbJVtEqmH1JP1XoL78l6SD+fKIwSLk4C1Nf0lH8m28krRe+oHfjcZ/Fbc3M1gt/v6vQfWRsd4P4pbkRtyeSjiW9MrNO7Mf+LO56192n7pPt7k+zv3V3i9uLC5S/q/CPxOP4N1tZUhp3rvPEm9vtzDUzs6m7npQ4aQ00s/Vc2S5yTAAgX5GvyFc4y93Z2BqxKczE5LnHa5JcUiduLqmfe31X0l7hGC5prbDPm9zj7DjbY85bdtxdSUfx9248xmqhrAe5fV1hoO+kmFfL9lVIqNv592HK93BX0m7xPSzsU0n5pzlXcb8p4x11zfrjyjOhrHvxOq/Negw2Nja2bCNfffI8+ersNSNfLfjGLIJoLHd/YWZDhVbBby5wqP3cMYdmJoUvyMzxFMfY02krVba2yUE8VplDD33zJ+lJGpbs+zJ3vqrNVP7YYrcu6bcKSa+j8/cRnzbesmvW0ez6kpZ1sc8RAJQiX5GvJPIVTtFFEE33jcKX5HLFx50mSY0z9NOuBSdb7vVpklWdZin/gUIXiK/c/a6k0i4VTWNhsc+Buz+RNIx95C+S/ACgDPnqcpCvkBwqWGi6HYWWoI1JO17yl9AThT72Umi96owb0HwO+yOOVdUMQsOS585d/rhvV9KGh/EG056r6LLjLfM0O7a7Z+foSifJDACqQL66GPIV+ao1qGCh0TzMorOj0CpYdCypZ2bdOKD09xWeup99mcXBtOuKt+tjd4EdSbtxpqJunOXn3Ku/xy//F/FYa2a2Ggf0rioM0L2oQYxhLZa1P0v5c0lqK9tfZ7uEnDnXiONcON5YhmlnauorzhgVdSQNYtJsesstgESQry6MfEW+ag0qWEjBlsr7M2dTlh7F36vsr3woacPM3ih8kW54bpYhd9+I59yO5/9aM64A7+7PFL7EtyR9rxDr45isLyQmp0OFgbhbuednKf+m4uBi5VrZJp2rpExVxNsfd46CYt/7TYV/QLpTjjsAgGmRr2ZEvpJEvmoNc5/LotkAUCkzexP71gMA0Fjkq8XDHSwAyYndJb6tuxwAAIxDvlpM3MECkJw4hqEzZgAzAAC1I18tJipYAAAAAFARuggCAAAAQEWoYAEAAABARahgAQAAAEBFqGABAAAAQEWoYAEAAABARahgAQAAAEBFqGABAAAAQEX+H/P3QA+qxThgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1518d88438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist_conv = []\n",
    "def cl(x):\n",
    "    hist_conv.append(x)\n",
    "res = NewtonLinConstraintsFeasible(f, gradf, hessf, A_lin, x0, backtracking, elimination_solver, (c, mu), callback=cl,\n",
    "                           max_iter=2000, beta1=0.1, rho=0.7)\n",
    "print(\"Decrement value = {}\".format(res[\"tol\"]))\n",
    "fstar = f(res[\"x\"], c, mu)\n",
    "hist_conv_f = [np.abs(fstar - f(descdir_x[1], c, mu)) for descdir_x in hist_conv]\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.subplot(1,2,1)\n",
    "plt.semilogy(hist_conv_f)\n",
    "plt.xlabel(\"Number of iteration, $k$\", fontsize=18)\n",
    "plt.ylabel(\"$f^* - f_k$\", fontsize=18)\n",
    "plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "\n",
    "hist_conv_x = [np.linalg.norm(res[\"x\"] - x[1]) for x in hist_conv]\n",
    "plt.subplot(1,2,2)\n",
    "plt.semilogy(hist_conv_x)\n",
    "plt.xlabel(\"Number of iteration, $k$\", fontsize=18)\n",
    "plt.ylabel(\"$\\| x_k - x^*\\|_2$\", fontsize=18)\n",
    "plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Реализация барьерного метода"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, linesearch, linsys_solver, \n",
    "                           tol=1e-8, max_iter=500, disp_conv=False, **kwargs):\n",
    "    x = x0.copy()\n",
    "    n = x0.shape[0]\n",
    "    mu = mu0\n",
    "    while True:\n",
    "        res = NewtonLinConstraintsFeasible(f, gradf, hessf, A, x, linesearch, linsys_solver, (c, mu), \n",
    "                                           disp_conv=disp_conv, max_iter=max_iter, beta1=0.01, rho=0.5)\n",
    "        x = res[\"x\"].copy()\n",
    "        if n * mu < tol:\n",
    "            break\n",
    "        mu *= rho_mu\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n",
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.23 s ± 90.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22.1 ms ± 2.02 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "[9.21322571e-09 1.85264514e-08 3.71311496e-08 7.00335912e-08\n",
      " 9.71878844e-08 8.95342930e-08 7.81315161e+04]\n"
     ]
    }
   ],
   "source": [
    "mu0 = 5\n",
    "rho_mu = 0.5\n",
    "x = BarrierPrimalLinConstr(f, gradf, hessf, A_lin, c, x0, mu0, rho_mu, backtracking, elimination_solver, max_iter=100)\n",
    "%timeit BarrierPrimalLinConstr(f, gradf, hessf, A_lin, c, x0, mu0, rho_mu, backtracking, elimination_solver, max_iter=100)\n",
    "%timeit BarrierPrimalLinConstr(f, gradf, hessf, A_lin, c, x0, mu0, rho_mu, backtracking, simple_solver, max_iter=100)\n",
    "print(x[:n])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Сравнение времени работы"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current dimension = 3\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current dimension = 4\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current dimension = 5\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current dimension = 6\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current dimension = 7\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current dimension = 8\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current dimension = 9\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "mu0 = 2\n",
    "rho_mu = 0.5\n",
    "n_list = range(3, 10)\n",
    "n_iters = np.zeros(len(n_list))\n",
    "times_simplex = np.zeros(len(n_list))\n",
    "times_barrier_simple = np.zeros(len(n_list))\n",
    "for i, n in enumerate(n_list):\n",
    "    print(\"Current dimension = {}\".format(n))\n",
    "    c, A, b, bounds = generate_KleeMinty_test_problem(n)\n",
    "    res = scopt.linprog(c, A, b, bounds=bounds, options={\"maxiter\": 2**max(n_list) + 1})\n",
    "    time = %timeit -o -q scopt.linprog(c, A, b, bounds=bounds, options={\"maxiter\": 2**max(n_list) + 1})\n",
    "    times_simplex[i] = time.best\n",
    "    A_lin = np.zeros((n, n + A.shape[0]))\n",
    "    A_lin[:n, :n] = A\n",
    "    A_lin[:n, n:n + A.shape[0]] = np.eye(A.shape[0])\n",
    "    x0 = np.zeros(2*n)\n",
    "    x0[:n] = np.random.rand(n)\n",
    "    x0[n:2*n] = b - A.dot(x0[:n])\n",
    "    time = %timeit -o -q BarrierPrimalLinConstr(f, gradf, hessf, A_lin, c, x0, mu0, rho_mu, backtracking, simple_solver)\n",
    "    times_barrier_simple[i] = time.best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAFMCAYAAABF6v+HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XdcXGXa//HPnYQ0UoB0TR3SqwKJbdVVYW0bdU2Iur2ZuD7Ps8XfGoxlXXvAdV23uEJ21+27CVjWHiH2siqg6ZVJj6kwCSn0+/fHGRIgQBgywxmY7/v1mheZM2fOuWYCc665y3Ubay0iIiIiTenkdgAiIiIS3pQsiIiISLOULIiIiEizlCyIiIhIs5QsiIiISLOULIiIiEizlCyIiIhIs5QsiIiISLOULIiIiEizurgdQLjo37+/HTlypNthiIiItJmCgoL91toBp9pPyYLfyJEjyc/PdzsMERGRNmOM2dqS/dQNISIiIs1SsiAiIiLNUrIgIiIizVKyICIiIs1SsiAiIiLNUrIgIiIizVKyICIiIs1SnYUAlZeXU1xcTGlpKdXV1W6HIxGic+fO9O7dm7i4OLp16+Z2OCISYZQsBKC8vJxt27YRGxvLyJEjiYqKwhjjdljSwVlrqays5NChQ2zbto3hw4crYRCRNqVuiAAUFxcTGxtL//796dq1qxIFaRPGGLp27Ur//v2JjY2luLjY7ZBEpI1Za3l15eeUV7nToq1kIQClpaX06dPH7TAkgvXp04fS0lK3wxCRNnTgcDk/+HshP/hHIYs/2e5KDOqGCEB1dTVRUVFuhyERLCoqSmNlRCJI3po93PHsSg4dq2TBleP52jkjXIlDyUKA1PUgbtLvn0hkOFxexQMvrmFx/nYmDOnD378/g/GD3WvZVrIgIiISRj7eXMxtSz5jl+8Yt34xnh8lj6Fbl86uxqQxC9KkvLw8EhMTMcYQGxtLSkoKeXl5bR5HamoqqampbX5eEZG2VF5VzSOvrOWGrA/pZAxL5p3H/CvGu54ogJIFaUJOTg4pKSkkJydTUFBAdnY2MTExpKenH98nMTGRlJQUF6MUEekYVu86yDW/eZ/Md7zcNGM4r/7oQpJGxrkd1nHqhpBGpaWlMX/+/HrJQXJyMj6f7/j9efPmuRGaiEiHUV1jeertIn6Vt4GYnl15+tvTuWT8QLfDOknEJwvGmJnAzNGjR7sdSljxer3Ex8eftD0mJub4v+fOnduWIYmIdChb9h/h/2Uvp2BrCVdPGcKD100mNrqr22E1KuK7Iay1L1pr5/bt29ftUMKKx+MhMzPT7TBERDocay3/+GgrV/36XTbuKeWJG8/it189O2wTBVCyIE3IzMyksLCQ2NhY5s2b1+jAxoYDD2vvZ2RkEB8fT2xsLBkZGfh8PlJSUjDGEB8fT2FhYb3npKSkHH+OMaZF3RtpaWnExsYSGxtLWlra8W3GGLxe7/H9MjIyiI2NrbdNRMQtew+V8Z0/f8Jdz60iYXgsS39yEdeedWbYT4uO+G6IYLnvxdWs2XXI7TDqmXhGH+6dOalVz01OTqaoqIi0tDSWLFlCVlYWHo+H3NxcPB5Pk8/Lyclh+vTp5ObmkpmZSVpaGpmZmWRmZpKens7NN99MamoqRUVFx5+Tl5dHSkoK2dnZLF68+HjiMH/+/EbPkZqaSmFhIdnZ2YAzdiI+Pp709HQKCwtJTU2loKCAwsJC0tLSThmziEhbeGnFLu5+fhVlldXcd80kvnHuCDp1Cu8koZZaFqRJHo+H7OxsSkpKyM3Npbi4+JSzH2JiYpg/fz4ej4cFCxYATuKRnJxMQkIC8+bNO+lbfnJyMvPnzychIYH09HRmz57dZBeI1+slJyeH7Ozs48dNT08/vn92djY+n4958+Zx2WWXMX/+fJKTk4PwboiItM7Bo5X86N+f8r///JQR/aJ5+YcX8q3zR7abRAHUshA0rf0G314kJyezbNkyEhMTycvLa/ICnJSUdPzftYMhExMTj2+Lizv1VKCUlBRycnIafay2O6TuMeuKiYkhNzeX+Pj444mEiIhb3t24j9uzV7D/cDm3pYzl1i/G06Vz+/uermRBWiwhIQGg3vTJlmhJghCImJgYSkpKmny8tuXC6/Xi8/nqzeAQEWkLxyqqWfjqWv7y4VbiB0ST9c3zmTq0/X4Wtb/0RtpEY9/sa7/V1yYNoZKdnd3kOWprPdQdJFmX1+slNTWV3NxcEhISVPlRRNrcp9tKuPrX7/KXD7fy3QtG8fIPL2zXiQIoWZBG1F5wU1JSyMrKIicnh4yMDFJSUpg7d27QBwvm5eWRlZWF1+s9PvOidrxDQx6Ph7lz55KamkpeXt7xMQy1YylSU1OZM2cOycnJZGdn4/V6j8+WEBEJpcrqGn75+npmP/UhZZXV/PPmc/jZzIl0j3K/XPPpUjeEnMTj8VBUVER6ejrp6el4vV48Hg/p6elNzlA4HQkJCRQUFJCWlkZcXByZmZnMnj27yf0zMzPJyMg4PlgyISGBBQsWMG/ePHw+X73BkdnZ2SQmJjJ9+vRmjykicjo27inlJ0s+Y9XOQ1yfcCY/v2YSfbpHuR1W0BhrrdsxhIWkpCSbn5/f7D5r165lwoQJbRRRZEhNTcXn85Gbm+t2KO2Gfg9FwkdNjeXpD7aQ/to6enXrwsNfmcwVk4e4HVaLGWMKrLVJp9pPLQsiIiKtsNN3jJ8uWc6H3gNcNn4gj8yawsDe3d0OKySULIiIiATAWsszhTu574XV1FhL+qwpzEkaFvZVGE+HkgUREZEWOnC4nDufW8nS1XuYMTKOx+ZMY1hcT7fDCjklC+Kq2pLNIiLhLnfNHhY8u4JDx6q486rxfO8LHjq3oyqMp0PJgoiISDNKyyp54KU1LMnfwYQhffj796cxfnAft8NqU0oWREREmvCR9wD/L3s5u3zHuPWL8fwoeQzdurT/ugmBUrIgIiLSQFllNb/M3cCid70Mj+tJ9i3nkTgiuKXr2xMlCyIiInWs3nWQ2xYvZ/2eUr56znDuumoC0d0i+3IZ2a9eRETEr6q6hsx3vPwqbwMxPbvy9Hemc8m4gW6HFRaULIiISMTbsv8Ity35jMJtPq6eMoQHr5tMbHRXt8MKG0oWREQkYllr+cdH23jo5bVEdTY8ceNZXDPtjA5dYKk1lCyIiEhE2nOojPk5K3h7wz4uHNOfjNlTGdK3h9thhSUtUS0nycrKwhhT7xYfHx92Sz2npqaSmpp62vuISOR5acUuLv/VO3y0+QD3XzuJv3xnhhKFZihZkCaVlJRgrcVaS2ZmJjk5OcTHx7fouYmJiaSkpIQ4wvAUya9dJNwdPFrJD//1Kf/7z08Z0S+al394Id88bySdIqQSY2upG0JaJDk5mYKCAmJjY8nIyGD+/PnN7j9v3rw2iiz8RPJrFwln727cx+3ZK9h/uJzbUsZy6xfj6dJZ35lbQsmCtFhMTAzJycnk5uaeMlmYO3duG0UVfiL5tYuEo6MVVSx8dR1//XArowf2YtE3k5gytK/bYbUrShYkIDExMfh8PrfDEBFpkU+3lXDbkuVs3n+E714wivlXjKN7VOSVaz5dShaC5dU7YPdKt6Oob/AUuHJhUA+Zl5fHggULjt9PTU0lLi6O1NRU0tLSuOGGG5g/f/7xQYW1q0rW3p8+fTqZmZkUFxezYMEC5s6dS2pqKnl5eXg8HrKzs0lISDh+rrS0NAoLC4mJiWHBggWnbNFojM/nIy0tjZycHIqLi5kzZw6ZmZn1XlNz52nuNbbktQOkpaWRlZUFOC0P6enpzR5bRE5PZXUNv162kd+9uYnBfbrzz5vP4fz4/m6H1W6ps0ZaxOv1kpKSQlxc3EkXs7y8PFJTU5k3b16zTfA5OTkA5ObmMnfuXNLS0khMTCQtLY2CggJiYmLqzVzwer2kp6dTVFREdnb28Qt+oPLy8oiPjyc3N5cFCxaQlZVFRkZGQOdp6jW25LWnpqaSk5NDdnY22dnZ5OTkHE8cWnoMEWm5jXtK+cqT7/ObNzbxlbOH8tpPLlKicLpqR7t3tBvgAdJbun9iYqI9lTVr1pxyn44gMzPTAjYmJsbGxMRYwAJ29uzZtqSkpN6+s2fPtoAtKio6afvs2bPr3Y+JiTl+v6SkxAJ27ty5J523KY0ds+79pp6TnJx80jaPxxPQeZp6jad67UVFRRawBQUFxx/Pzs62CQkJzR6jOZHyeygSqOrqGrvonSI75q5X7Nn3v25fXfm52yGFPSDftuAa2SG7IYwxswHNXTtNy5Ytw+PxAM5YhaYkJCQc3685SUlJx/9de7zExMTj2+Li6q/o5vP5yMrKYvHixXi9Xnw+H8nJyQG9hsakpKTUazloyXmaeo2neu15eXlA/dcZ6DFE5NR2lBzlp9nL+a+3mOQJA3nk+qkM6N3N7bA6jA6ZLFhrc4wxPkDVeE6Dx+NpNkmou19rNUwQ6kpMTMTj8bBo0SISEhJITU0NyeDKlpynqdfYktceExNDSUlJk48rURBpPWstzxTu5L4XVlNjLRmzppKaNFTlmoOsQyYL0v4VFhbi9XrJzc0N+sW0oKDg+CDKUJ4HnPoUPp+PwsLC4+cUkeA4cLicO59bydLVe5gxMo7H5kxjWFxPt8PqkFwZ4GiMSTDGZDfx2FxjzGz/TcPCI1TthTUtLQ2v10tOTk6rBjeC0xVQ2x2QlZVFVlbW8RkdwTxPYzweT70ZH7XnUIVHkdOTu2YPl//qHd5ct487rxrPv+aeq0QhhNq0ZcEYkwDc4L970tc4Y8xccLoR/Pc9xphMa+08//35QL8mDp9rrc0LftTilvT0dNLS0sjLy2POnDmtHq+QkJBAZmbm8WmKmZmZzJ49O+jnaUpmZiYZGRnMmzcPr9dLQkJCvemnItJypWWVPPDSGpbk72DCkD784/tnMW5wb7fD6vCMMxiyjU/qJA2LrLWJDbYXNLKtyFrbsgUJ6j8vGUitTTROJSkpyebn5ze7z9q1a5kwYUKgoYgElX4PJVJ95D3A/8tezi7fMW65OJ4fJ4+laxdVADgd/utu0qn2C5sxC8aYGKCxTl2fMSZZrQYiIpGprLKax15fzx/e28zwuJ5k33IeiSOaHhwtwRc2yQJOt0RjQ92LcZKIFicL/laFeUCCMWautTbrVM8REZHws2rnQW5b8hkb9hzma+cM586rJhDdLZwuXZEhKO+4MeZ2YKG19nQKbsfhJAYN+Wh6nEKj/K0Qp0wu/GMk5gIMHz48kFOIiEgIlVVW89TbRfzuzU3E9uzK09+ZziXjBrodVsQKVnqWQ+OtAmHN3+KQBc6YBZfDERGJeNZaXli+i/RX17HrYBkzp53B/ddMIja6q9uhRbSgJAvW2s3AoiAcqrFOqBjgQBCOLSIiYaxgawkPvLSGz7b7mHRGH355w1mc6wmoYVlCJJw6fvJxEoOG4oDCNo6lSdZaVQYT17gxe0kk1HaUHCX9tfW8uHwXA3p349HZU5mVMJROnfRZGy4CShaMMd8HEq21P2iwPR/4vbX2j60NxFrrM8Z4jTEx1tq6XRox4TITonPnzlRWVtK1q5rDxB2VlZV07nw6Q4NEwseR8ip+/1YRi971AvB/l47mlovjNYAxDAX6P5LmvzW00L+9pclCU3Ne0oEFtefw12MIi0QBoHfv3hw6dIj+/bXUqbjj0KFD9O6tAjTSvlXXWJ4p2MGjr69nX2k51551BvOvGM+ZMT3cDk2aEGiyEA94G9leROM1EuoxxnhwpjQm40xrzAQKaqc2Wmuz/OWek3G6JDwtLarUWsaYmcDM0aNHn3LfuLg4tm3bBkCfPn2IiopSl4SEnLWWyspKDh06RElJiWbuSLv2YdEBHnhpDWs+P8TZw2PI/EYiCcNj3Q5LTiGgCo7GmE3A7dba5xpsnwWkW2tPfcUNUy2p4AhQXl5OcXExpaWlVFdXt0FkIk4XWO/evYmLi6NbNy27K+3Plv1HePiVtby+Zg9n9O3OHVdNYObUIfrC5bJQVXDMAv5gjPFZa9/0n+hS//ZHAg+z/enWrRtDhgxhyJAhbociIhL2Dh6r5LdvbOTPH2whqnMnfvqlsXz/Qg/dozT2pj0JKFmw1mYYY/oBy4wxdZskHrXW/iK4oYmISHtVVV3Dvz7exuN5Gyk5WkFq4lB++qVxDOzT3e3QpBUCHnJqrU0zxjyMM+4AIM9aezC4YYmISHv11vq9PPTyWjbuPcy5njjuvnoik8/s63ZYchoCThb80ydvAc6y1nbxb1sKPNVwLIOIiESOTXtLefDltby1fh8j+vUk8xuJfGniII1L6AACrbOwEJgF3AEsqfPQIuBOoN0lC4HMhhARkZMVH6ngV3kb+MdH2+jZtTN3XTWBb54/gm5dNC6howh0NkQxMNta+4Yxprp24ShjzChg02kuJOWqls6GEBERR0VVDX/9cAu/XraRIxXVfHXGcH6cPIZ+vTRjp70I1WwIgNrsom67kgfY3IpjiYhIO2OtJXfNHh5+ZS1bDhzlorEDuPvqCYwdpIJhHVWgycIiIMsYk4o/aTDGnA085b+JiEgHtnrXQR58aS0feg8wemAvLR0dIQKdOpnmr8JYCGCMOQDEAlmaOiki0nHtLS3jsaUbWFKwnZgeUdx/7SRumjGcqM6d3A5N2kBrpk6m+hOGy/ybcq21W4IalYiIhIWyymr++N5mnnxzExXVNXz3glH88NIx9O0Z5XZo0oYCnQ3xU8BrrX0W8BpjlgBPGWOKgDnW2s9CEaSIiLQtay0vrficha+uY6fvGCkTB3HnVRMY1T/a7dDEBYG2LCwAUgGMMbfjTKOcA9yIM55helCjExGRNvfZdh8PvLSGgq0lTBjSh0dnT+X80VptN5IFmizEcmLVyRQgx1r7jDFmM/BJUCNrI6qzICLi2OU7RsZr63j+s13079WN9FlTmJ04jM6dVFQp0gWaLPiAvsaYvjjlnmf7t8f6H2t3rLUvAi8mJSXd7HYsIiJuOFJeRebbRWS966XGwq1fjOfWS0bTq1trZtdLRxTob8JC4FOghBNjF8DpmlgWzMBERCS0amosz366k0eXrmPPoXK+PHUIaVeMZ1hcT7dDkzDTmlUnC4G+QF6dh3LxT6cUEZHw9/HmYh54aQ0rdx5k2rAYnvxaAokj4twOS8JUa6ZO5jWy7ZnghCMiIqG07cBRHnl1La+u2s2Qvt351Q1ncc20M+ikcQnSDHVIiYhEgENllfzuzU08/d4WOncy/CR5LHMv8tCja7td0kfakJIFEZEOrKq6hsX52/nl6xs4cKSCWQlDuf3ycQzu293t0KQdUbIgItJBvbtxHw++tJb1e0qZMTKOP39nIlOG9nU7LGmHlCyIiHQwRfsO8/DLa1m2bi/D4nrw+68lcMXkwRijcQnSOhGfLKgok4h0FL6jFfwqbyN//+9Wukd15o4rx/Pt80fSPUrjEuT0RHyyoKJMItLeVVbX8LcPt/LEso2UllVy44zh3JYylv69urkdmnQQEZ8siIi0V9Za3li3l4deXot3/xEuHNOfu66ewPjBfdwOTTqYoCQL/kWlFuKUfF4CZGoFShGR0Fm3+xAPvrSW9zbtxzMgmj99O4lLxg3UuAQJiWC1LOQAPmvtImPM2cDBIB1XRETq2H+4nMde38DiT7bRu3sU986cyNfPHUFU505uhyYdWFCSBWvtZpwlqrHWfhqMY4qIyAllldU8/f4WfvfmJsoqq/nW+SP50WVjiOnZ1e3QJAIEnCwYY74P3AKcZa3t4t+2FHjKWvtckOMTEYlo1lpeXbWbR15dy/biY1w2fiB3Xj2B+AG93A5NIkhAyYIxZiEwC7gDyK7z0CLgTkDJgohIkKzY4ePBl9by8ZZixg/uzd+/dw5fGNPf7bAkAgXasjAXmG2tfcMYY+tsLwASgheWiEjk2n2wjIyl63i2cCf9orvy8FemcMP0YXTWYk/iktaMWahNEur+1nqAzacfjohI5DpWUU3WO16eeruI6hrLLRfHc+sl8fTpHuV2aBLhAk0WFgFZxphU/EmDf/bDU/5bu6MKjiISDpat3cO9L6xmR8kxrpoymDuumMDwfj3dDksEAGOtPfVedZ9gTDbOuAWAEiAGyLLW/iDIsbWppKQkm5+f73YYIhJhdpQc5b4X15C7Zg9jBvbi/msnc158P7fDCg/Wgq2BTipXHSrGmAJrbdKp9gu4G8Jam2qM8QCX+Tfl+adOiohIC1VU1fDH9zbz62UbAbjjyvF894JRdO2ieglUlsGnf4P3fw2HdkL0AOg1EHoNcm69/T/rbus1CLpphkiotKrOgrXWC3iDHIuISET4sOgA9/xnFZv2HubySYP42cxJnBnTw+2w3FdxBPKfhg9+DYf3wNAZMDUVjuyD0j3Otj2r4cheqKk6+flR0U4C0XtwnURiIPQaXD+5iB4AnbXaQSACnTp5Nk45Z09jj1tr1VYkItKEfaXlPPzKWp77dCfD4nrwp28ncen4QW6H5b6yg/BxFnz4JBwrhlEXwaw/wMgLobHy1TU1zn6H/QlEbSJxeO+JbXvXgvct59gnMRDdv0EyUSe56F0nuejWp/EYIkygqVU2zvoPtwDFwQ9HRKTjqa6x/OOjrTy6dD1lldX836WjufWLo+nRNcK/Xx05AB/9Hj7KgvKDMOZyuOinMGxG88/r1Mm52Ef3h0GTmt+3suzkROJwg+Ri3wbnZ03lyc/v0qNBK8WgxlsuogdAl45bTTPQZMEDJGiRKBGRllm+3cfdz69i5c6DXDC6H/dfO1nVF0v3wIe/gU/+BJVHYMI1TpIwZFrwzxXVHWJHOLfmWAvHSvwJxO4TiURpnX8fKIKt7zv7NaZHXBNjKhokFz1i211rRaDJQh4QG4pAREQ6koNHK3n09XX846NtDOjVjd/cdDZfnjoksleF9G13xiMU/MX5Fj95Nlx4Gwyc4HZkzsW7Z5xzGzi++X2ryp1xFE11gRzeA9s+dB6rLj/5+Z27Nj5As2HLRfRAJ9kJA4EmC6lAvjEmFyikQVeEtfbZYAUmItIeWWt5tnAnD7+ylpKjFXz7/JHcljKW3pFcWKnYC+89Dp/9y7l/1k1wwY+hX7y7cbVWl27Qd6hza461zpiJ5rpASrbC9o/h6P7Gj9G9b/2WienfhxHnBf81nUKgycLNQLz/1pAFIrwDTkQi2YY9pdz9/Co+3lzM2cNj+Ov3ZjDpjL5uh+Wevevg3cdgVQ50ioKk78D5P4SYYW5H1jaMgR4xzm3A2Ob3ra480VpxvAukQXKxMx8mz2r+OCESaLJwJ5AOLLTWNjbEVEQk4hwpr+LXyzbyx/c206t7FxZeP4U5ScPoFKlrOXy+HN75Bax9EaJ6wnn/A+f9r9O8Lo3rHAV9znBuYSjQZKEYyFSiICLidDksXb2H+19cza6DZcxJGsodV04gLrrjjopv1vaPnSRh41JnyuFFP4VzfgDRqkjZ3gWaLGQCs4FfhCAWV2htCBFpjW0HjnLvC6t4c/0+xg/uza9vOpukkXFuh9X2rIUt78I7j8Lmd5wZAZfeAzNudvrbpUMIaG0IY8wSnHUhfMBJCylYay8PXmhtS2tDiEhLlFdVk/m2l9+9uYkunQw/SRnLt88fSZfOEVam2VrYlOckCds/cgbfnf9DSPy2yi63IyFbGwJ4phXPERFp997duI+f/Wc1m/cf4eqpQ7jn6okM7hseU9vaTE0NrH/ZSRI+Xw59h8FVv4CzvxE20/wk+AJKFqy1c0IViIhIuNpzqIwHXlrDSys+Z2S/nvz1uzO4aOwAt8NqW9VVsPo5Z3bDvrUQ54FrfgtTb+jQlQvFoZU0RESaUFVdw18+3MrjuRuoqK7hJ8ljmXexh+5RETRLvKoCViyG937p1EsYMAGu/wNM+ooWY4ogzf5PG2OWAtnW2j/47y9ubn9r7Q1BjE1ExDUFW4u567lVrNtdysVjB3D/tZMY0S/a7bDazvFlop+Ag9udUsw3/B3GXe2szSAR5VRpYcP5LhE6aVhEIkXJkQrSX1vHvz/ZzpC+3Xnq6wlcPmlw5JRpbmyZ6C8/DqOT2916BhI8zSYLDUdIasyCiHRUNTWW7ILtLHx1HaVlVcy7yMMPLxtDdLcIaWoPdJloiSgB/RUYY65vbP0HY8ylANbaN4IVmIhIW1mz6xB3P7+Swm0+po+M5cHrpjBucG+3w2obrV0mWiJKoClzNo2v/2CAhcD0045IRKSNlJZV8njuRv7y4Rb69ojiF6nTmJVwZmR0ObTlMtHS7gWaLDT1FxQDJJxmLCIibcJay8srP+eBl9awt7Scr84Yzu2XjyOmZwRMAQznZaIlbLUoWTDG1OCsKmmNMdWN7QIUBDMwEZFQ8O47zL0vrObdjfuZfGYfMr+RxFnDYtwOK/Q62jLR0qZa2rKQgpMQvA6kNnjMB3ittZuDGZiISDCVVVbz5JubeOptL926dOK+aybx9XNH0LmjrwwZ6ctES1C0KFmw1i4DMMbkWWtV7llE2pU31+3lZy+sYnvxMa476wzuvHoCA3t38NLEWiZagijQcs9fClUgIiLBtst3jPteXM3S1XuIHxDNP28+h/Pj+7sdVmhpmWgJgQiZQCwikaSyuoY/vbeZJ5ZtpMZabr98HDdf6KFrlw5aedBa2PKef5not7VMtARdxCcLxpiZwMzRo0e7HYqIBMFH3gPc859VbNhzmOQJg7h35kSGxfV0O6zQaGyZ6C89pGWiJeiMtdbtGMJCUlKSzc/PdzsMEWml/YfLefiVtTxbuJMzY3rw82smkTJxkNthhUZjy0Rf8CMtEy0BM8YUNKzW3JiIb1kQkfatusbyz4+38ehr6zhWWc3/XBLP/14yhh5dO+DKkFomWlyiZEFE2q2VOw5y9/MrWb7jIOd5+vHAdZMZPbADNr9rmWhxmX7LRKTdOXisksdeX8/f/ruVftHdeOLGs7hm2hkdr0yzlomWMBFwsmCMGQl4cEo819PYIlMiIsFireX5z3by0MtrKT5SwbfOG8ltXxpLn+5RbocWXFomWsJMoKtO3o6zYFRjv62WxheZEhE5bZv2lnL386v4r7eYacNi+PN3ZjD5zA42LfDIASh4Gj6WHnnjAAAgAElEQVT8nZaJlrASaMvCAuBR4BFr7cEQxCMiUs/Riip+88YmFr3jJbpbFx76ymRumj6cTh2lTHNlGWx4DZb/GzblQk2VlomWsBNoslAMPKVEQUTawuurd3Pfi2vY6TvG7MSh3HHlePr36uZ2WKfPWtj2X1jxb2d2Q9lB6DUYzv0BTPsqDJrodoQi9QSaLGQCs4FfhCAWEREAthcf5b4XV5O3di/jBvVmybzzmDEqzu2wTl+xF5YvdpKEki3Omg0TZjpTHz1fhE7qyZXwFGiyMB2YZYxZAJxUwchae3lQohKRiFReVc0f3t3Mb97YSCdjuOuqCXz7gpFEdW7HI/+PFjutBysWO1UWMc5YhIvvgAlfhm693Y5Q5JRaM3VSq06KSNC9v2k/9/xnFd59R7hqymDu+fJEhvTt4XZYrVNV4Yw/WP4v2LAUqitgwHhI/jlMmQN9z3Q7QpGABLrq5JxQBSIikWnvoTIefHktLyzfxfC4njz9nelcMm6g22EFzlrYWegkCKuecWYzRA+ApO/BtBudGgma0SDtVGvqLPQB0oFk/6ZcIM1aWxrMwESkY6uqruFv/93KL1/fQHlVDT+6bAw/+GI83aPaWb+9b5vTxbD833BgE3TuBuOvdhKE+EuhcwerASERKdA6C32BzYCXE90RyYDPGOOx1m4Ncnwi0gF9uq2Eu59fxepdh7hwTH/uv3Yyo/pHux1Wy5UdhDX/cQYrbn3P2TbiAmcxp4nXallo6XACbVnIApZYa2+pu9EYk40zU+KKYAUmIh2P72gF6a+t59+fbGNg72787qsJXDVlcPso01xdBUVvON0M61+BqjLoNxouuRumzoHYEW5HKBIygSYLKUBCI9sfBvJOPxwR6YhqaizPFO7gkVfXcfBYJd+7YBQ/ThlLr25hvjyNtc4S0CsWw8psOLIPesQ6S0FPuxHOTNQ4BIkIrSnKNArY0mC7x/+YiEg963Yf4p7nV/HJlhISR8TywLWTmXhGH7fDat7BnbByidPNsG8tdO4KYy+HaTfB6BQtBy0RpzXdEFnGmOTa8QnGmLNw1ovICXZwItJ+HSmv4ollG/nje5vp070LGbOmMjtxaPiWaS4/DGtfdLoZNr8DWGcBp6t/6SwF3bMDFIUSaaVAp05mGGPigc3GmBL/5hggx1q7IOjRiUi7Y63ltVVOmebdh8q4acYw5l8+ntjoMPw2XlMNm992ZjKsfREqj0LMCLh4vlNVsV+82xGKhIWAOwyttfOMMZlAon9TvrX20+CGJSLt0Zb9R7j3hdW8vWEfE4b04XdfSyBxRKzbYZ1sz2onQViZDaWfQ7e+ziDFqTfC8HM1DkGkgVaNLrLWFgKFQY5FRNqpsspqnnq7iCffKqJr50787MsT+eZ5I+gSTmWaS/c4ycGKf8PuldCpizP+4IqFMPYKiOrudoQiYavZZMEYsxTIttb+wX9/cXP7W2tvCGJsItIOvLNhHz/7zyq2HDjKl6cO4Z4vT2RQnzC58FYcdaY5Lv+XM+3R1sAZCXBlBkyeBdH93Y5QpF04VctCvwb3O1zbnDFmJjBz9OjRboci0q7sPljGAy+t4eWVnzOqfzR/+94MLhwzwO2woKbGKZS0fLFTOKmiFPoMhS/8xOlmGDDW7QhF2h1jrXU7hrCQlJRk8/NPWkhTRBqoqq7hzx9s4fHcDVTVWP73ktHMvdhDty4ul2net8HpYlixBA5uh669nWqK0250qit2CqMuEZEwYYwpsNYmnWq/QMs9X2+tfbaR7ZcCWGvfCOR4ItK+5G8p5u7nV7FudymXjBvAfddMZni/nu4FdGS/s2jT8n/DrkIwnZz1GJJ/DuOugq4uxibSgQQ6wDEbaOzrg8GptTD9tCMSkbBTfKSCha+uZUn+Dob07c5TX0/k8kmD3CnTXFkGG15zEoRNuVBTBYOnwJcegimzoffgto9JpIMLNFlo6pMhhsbLQItIO1ZTY1mSv52Fr63jcFkV8y728MNLxxDd1mWarYVt/3W6GVY9B+UHoddgOPdWp5th0KS2jUckwrToL94YUwNYwBpjqhvbBSgIZmAi4q41uw5x9/MrKdzmY8bIOB64bjLjBvdu2yAOFDnrMqxYDCVbIKonTJjpJAijLoZO7Ww5a5F2qqVfD1JwEoLXgdQGj/kAr7V2czADExF3lJZV8njuRv78wWZie3blsdRpXJ9wZtt1ORwthtXPOd0MOz4GDHguhovvcBKFbr3aJg4ROa5FyYK1dhmAMSbPWvtMaEMSETdYa3lpxec88NIa9h0u56szhnP75eOI6dkGZZqrKpzxB8v/BRuWQnUFDJgAyffBlFToe2boYxCRJgW6NsSXQhWIiLjHu+8w976wmnc37mfymX3I+mYSZw2LCf2JD30O7//Kme54rBiiB8D07zvdDIOnquyySJgIeJSSMeZ64AacQY31WGsvD0ZQItI2yiqrefLNTTz1tpduXTpx/7WT+No5I+gc6pUhj5XAe7+CjzKd2QwTZjrLP8dfCp3bePCkiJxSoHUWbgfScZajTuHEstSzgdzghiYiofTmur387IVVbC8+xnVnncGdV09gYO8Ql2muPOYkCO89DmUHnS6GS+6EuFGhPa+InJZAU/g7gGRr7RvGmGJr7RwAY8xs4LKgRyciQbfLd4z7XlzN0tV7iB8QzT9vPofz40O8RkJ1FXz2D3hrIZTuchZwSr7XqY8gImEv0GQhFijy/9trjJlmrV2OM20yE/hBMIMTkeCprK7hT+9t5ollG6mxltsvH8fNF3ro2iWEZZCthbUvwhsPwP4NcGYSzFoEI78QunOKSNAFmiwUAvHAVmAJkGGMmYfTNdEGo6FEpDU+8h7gnv+sYsOewyRPGMS9MycyLC7EpZA3vwt5P4ed+dB/LNzwdxj/ZQ1aFGmHAk0W0gAP8AZOS8I8nJYG439MRMLI/sPlPPLKOp4p3MGZMT1Y9M0kUiYOCu1JP18By+6DTXnQ+wy45jcw7asauCjSjgU6dXJZnX8fBOKNMWdbaz8NemQi0mrVNZZ/fbyNjNfWcayymlu/GM//Xjqanl1DeMEu3gxvPgQrs6F7DKTcDzPmQlSP0J1TRNrEaa86aa39VKtOioSPlTsOcvd/VrF8u4/zPP144LpJjB4YwjLNh/fCO49C/tPQqQt84SdwwY+hh3omRToKrTop0kEcPFbJL19fz9/+u5W46G786oazuPasM0JXprnsEHz4W/jgt1BVBgnfcEoy9xkSmvOJiGu06qRIO2et5T+f7eLBl9dSfKScb5w7gtu+NI6+PaJCc8Kqcsj/k9OacPQATLwOLr0b+o8JzflExHVadVKkHdu0t5R7nl/Nh94DTBval6e/PZ0pQ/uG5mQ11c54hDcfAt82GHURJP8czkwMzflEJGxo1UmRduhYRTW/eWMji9710iOqMw9eN5mbZgwPTZlma2Hj65B3H+xdDUOmwcwnwHOJpkGKRAitOinSzuSt2cO9L6xmp+8YsxKGsuCq8fTv1S00J9v2kVMrYdsHEDsKZv8JJn4FOoWwkJOIhB2tOinSTmwvPsp9L64hb+0exg7qxeK553KOp19oTrZ3LSy7H9a/AtED4erHIOFb0DlE4yBEJKwFPHWyuccbTqsUkdNXUVXDone9/OaNjRgMC64cz3e/MIqoziH4du/bDm89Asv/BV17OQMXz70VukYH/1wi0m4EOhsip4nt1v+zsWmVItJKHxTt557nV1G07whXTBrMPTMncmZMCIocHS2Gdx+DjxcB1kkQvnAbRIeo5UJE2pVAuyFO+ipjjIkBsoCHgxWUSKTbW1rGwy+v5fnPdjEsrgdPf3s6l4wfGPwTVRyB/z4J7/8aKg7DtJvgiwsgZljwzyUi7dZp13611vqMMWnAU8Dlpx+SSOSqrrH8/b9b+cXS9ZRX1fDDS0dz6yWj6R4V5Ea76koo/Au8nQGH98C4q+Cyn8HACcE9j4h0CMEsFJ8cxGOJRJzPtvu4+/mVrNp5iC+M7s/9107CM6BXcE9SUwNrnoM3HoRiLww/D+b8FYafG9zziEiHEugAx6VNPJSEU29BRAJ08GglGUvX8c+PtzGgVzd+c9PZfHnqkOCXaS56w5kG+flyGDgRbloMYy9XrQQROaVAWxYONrE9G2fJahFpIWstzxbu5OFX1lJytILvnD+Kn6SMoXf3IE9P3FnoJAmb34a+w+G6p2DqHOik8cgi0jKBDnCcE6pAgskYMxuIA+KBImttlsshidSzYU8pdz+3io+3FHP28Bj++r0ZTDojyGWa92+CNx6ANc9Dz35w+SMw/XvQJUQFnESkwwrh4vbuMMZ4AJ+1Nsd/v8AYk2+tLXQ5NBFKjlTw2zc38ZcPttCrexcWXj+FOUnD6BTMMs2HPoe306Hwr9ClO1w0H87/P+jeJ3jnEJGIEnCyYIz5KTAP8Pg3FQIPW2ufC2ZgpyEBuAHI89/Pxxl8qWRBXHOsopo/vb+Zp94q4khFFamJw0i7cjxx0V2DeBIfvP8E/Pf3UFPltCJcdDv0CsGUSxGJKIEOcHwKuBl4Bsjwb54D5Bhj5lpr/xjk+AJmrc0xxuTV2ZSEM6ZCpM1VVdeQXbCDx3M3sLe0nOQJg5h/xTjGDuodvJNUHnOKKb37GJT5YEoqXHIXxI0K3jlEJKIF2rIwF0i31i6os22RMSYTuANoUbJgjEkAFlhrG65giTFmLlDsv+ux1mY03OdUrLU+/7E8QLG1Nu8UTxEJKmstS1fvJmPperz7jpA4IpbffS2B6SPjgneS6iqnLPNbj8ChnTA6GS67F4ZMDd45REQIPFnwAosb2b6Qk5euPok/SbjBf9fTyONzwWkd8N/3GGMyrbXz/PfnA03Vn81tJClIs9amnCoukWD6yHuAR15dx2fbfYwe2IusbySSMnFQ8KZCWgvrXoJlD8D+9XBmInwlE0ZdGJzji4g0EGiykInT//9Zg+2X0XgSUY9/kGGhP2lorIjTPGttYp39vcaY5Dr3W9zK4E8s0vz/jqltbRAJlXW7D5Hx2nreWLeXwX26kz5rCrMShtIlmAs+bXnPmQa54xPoNwbm/A0mzFStBBEJqUCThRTgMmNMw2/ryThJQL2iTdbaFpd/9q8xkdDIQz5jTHIgXQn+BCOrToKQTNOLYImclh0lR3k8dyPPfrqD3t26kHbFeL59/kh6dA1iHYPdKyHvPtiUC73PgJm/hrO+Bp073IQmEQlDgX7S+HAGNzbU2LZAeWi8CmQxThLRomTB32qR6/937ebEJp8g0kolRyr43Zub+OuHW8HA3As9/OCL8cT0DOIMh5It8MZDsDLbmfqYfB+cMw+iQrDypIhIE8KpKFMcJwY21uWj6XEKJ/F3dbSoTdY/RmIuwPDhw1t6ColwDadBzkoYyk9SxnJGMJeOPrwP3nkU8v/kVFq84EfwhR9Dj9jgnUNEpIUiug3TX9kxCyApKcm6HI6EudppkL/K28CeQ+UkTxjI7ZePZ9zgIE6DLC+FD34LH/7WmRJ59tfhi3dAnzOCdw4RkQAFWmfhbGAJjcxkALDWnm4nbWPzymKAA6d5XJFWc6ZB7iFj6brj0yB/+9UgT4OsKof8p53WhKP7YcI1zpLR/ccE7xwiIq0UaMtCNk63wC003mVwOvJxEoOG4lD1RXHJR94DLHxtHZ9uC9E0yJoaZzzCmw+CbxuMvNAZlzBUw2xEJHwEmix4gARrbcOpk6fNWuszxngbmeYYo6JK0tZCPg3SWtiYC8vugz2rYPAU+PozEH+ZpkGKSNgJNFnIA4Ixwqqp9tt0YAEn6iO0eBaESDDs9B3jl69v4NlPd9DrdKdBWgsVh501G8p8dX6WwPJ/w9b3IXYkzPojTLoeOgWxHoOISBAFmiykAvnGmFycroF6XRHW2mebe7K//PI8nLoHCf4y0QW1S0hba7OMMXP9dRJicMo9zwswxoAYY2YCM0ePHh3K00iYKzlSwZNvbeIvH24F6kyD7BEF5YegpOEF3wdlBxvZVudn2UFnQafGRA+Aq34BCd+CLkGcaikiEgLG2pZPAvCvONlUFUUbhAGOrklKSrL5+fluhyGhUFPjXPAbuaBXHi5hxaYtFG3fQc/qw4ztU8XIXpV0rTh04oJva5o+tukMPWKge0xgP3sNVpIgIq4zxhRYa5NOtV+gLQt34nQVLLTWHmxVZCKtUVPdxDf5U3y7P+ZzEoUmLvhRwFTbmTGde9O9bxxde8VBjwHQfUzLLvpde2mMgYh0eIEmC8VAphIFaZXaC/6xklNf4Ov9POhc8GmmFaxTVP2LePQAZ9phg4u77d6Xj3dbFn1SzOriTowaegY/vupsZnhaXPdLRCTitGYhqdnAL0IQi3QUVeWwdw18vgJ2r3B+7l/vJArN6dy1/sW992AYML5l3/Cjep7yG37daZDxA0Zy3zfGB3capIhIBxVosjAdmGWMWYBTF6GeQBaOkg6i7CDsXnUiKdi9AvatOzGwr2svZ1rg5NnOt/1mL/ihWe9g/e5SMl5bx7J1exnUp1toVoMUEenAWlPuORiLRkl7VLrbnxAsP5EYlGw58Xj0QBgyFcZ8yfk5eCrEjnJtSmBQp0GKiESwcFpIyhWaOtmImhoo2QyfL3eWRq5tNTiy98Q+sSOdZODsr8PgaU5y0HuwayHX1XAa5M0Xerg12KtBiohEkICmTgIYY/rgzIhI9m/KBdKstaVBjq1NRezUyaoKp9ugbjfC7lVQ4f/v7NTFGTcweOqJ1oLBk6F7X3fjbsSximqe/mAzv3+riCPlzmqQP04Zy5nBXA1SRKQDCcnUSWNMX2Az4OVEd0Qy4DPGeKy1WwOOVNpOeemJ8QW1ycHetVBT6TweFe0kAtNudMYZDJkKAyZAVHd34z6Fquoacgp28HgoV4MUEYlggY5ZyAKWWGtvqbvRGJONM1PiimAFJqfp8L76Yws+XwHFXo5PP+zZz2klOO9Wf6vBNIjzQKf2059fuxrko0vXUbTvCAnDY/jNTQnMGBXE1SBFRCTgZCEFSGhk+8NoDQd3WOsMMqzXjbASSj8/sU/McCchmHajvxthCvQ5o10XE/p4czGPvLrWPw0ymsxvJPIlTYMUEQmJ1hRlGgVsabDdQ/CXrJaGqith/4b6rQW7V0K5v36B6QwDxsGoi+qMMZgCPYKx9ld40DRIEZG215puiCxjTHLt+ARjzFnAQiAn2MFFtIojsGe1f0ZCnfEF1eXO4116wKBJMGXWicRg4MSQ1Spw207fMR7P3cAzhZoGKSLS1gKdOplhjIkHNhtjSvybY4Aca+2CoEfXBsJi6uTR4vpJwe4VcGDTifUMesQ6CcE5c/3dCFOh32jo3JoyGe2LpkGKiLgv4KmTAMaYBCDRfzffWvtpUKNyQZtMnbQWDm5v0I2wAg7tPLFPn6EnpijW/uw7tF2PL2iNutMgD/unQf5E0yBFRIIqVKtOAmCtLQQKW/PciFFdBQc21hl06E8OynzO46YT9BsDI86vnxj0jOyR/JoGKSISfppMFowx1wPTW9K9YIz5BPi+tXZ5MINrt6oqIMNzorBR527O+IKJ1/qTgmnO/a493Y0zjFhreX3NHjJe0zRIEZFw01zLwp3AJy08zhtABqCFpAC6dIULb3OmJw6eCv3HRsT4gtb6eHMxC19dS6GmQYqIhKXmrmAeYH4Lj/M68P3TD6cDufA2tyMIe+t3l/Lo0nXkrXWmQS68fgqzEzUNUkQk3DSXLMTQ8toJxf79RU6pdhrks4U7iNY0SBGRsNdcsuAFLgM+a8FxktGARzkF39EKnnyriD9/sAWA72sapIhIu9BcspAF3GmMyWlugShjzChgAU7JZ5GTlFVW8/T7W3jyrU2aBiki0g41mSz4CzDdABQaY+Zba//YcB//jIlFwCfW2l+EMM6QCYuiTB3YR94DLHhuJd59RzQNUkSknWq2KJN/Seo/ALOAEpyuCS/O+IQk/89ngJuttQdDHm0ItUlRpghy8FglC19dx78+3sawuB48dN0ULho7wO2wRESkjqAUZfInAKnGGA+QhrOIVCJOwrAIWNwRqjdK8FhreW3Vbu59YTX7D5cz9yIPP04eQ8+umjoqItJetegT3FrrBeaFOBZp53YfLOOe/6wid80eJp3Rhz9+azpThvZ1OywRETlN+ronp62mxvKPj7eR8eo6KmtqWHDleL73hVGqlyAi0kEoWZDTsnFPKQueXUn+1hK+MLo/D31lMiP6RbsdloiIBJGSBWmV8qpqfv9WEb97cxPR3brwWOo0rk84UyWaRUQ6ICULErD8LcXc8exKNu09zLVnncE9X55I/17d3A5LRERCRMmCtFhpWSXpr63j7//dxpkxPXj6O9O5ZNxAt8MSEZEQU7IgLfL66t387D+r2VtaxncvGMX/+9JYorvp10dEJBJE/Ke9Kjg2b++hMu59YTWvrtrN+MG9yfxGItOGac0wEZFI0mwFx0iiCo711dRY/v3Jdh55dS3lVTX8OHkMN1/oIUrTIUVEOoygVHCUyFS07zALnl3Jx5uLOdcTxyPXT2VUf02HFBGJVEoW5LiKqhqy3ini129sonuXTqTPmsKcpGGaDikiEuGULAgAn24r4Y5nVrJ+TylXTx3CvTMnMrB3d7fDEhGRMKBkIcIdLq/iF0vX85cPtzC4T3f+8M0kkicOcjssEREJI0oWItiytXu45/lVfH6ojG+eO4LbrxhPL02HFBGRBnRliED7Ssu578XVvLTic8YO6kXOV88ncUSs22GJiEiYUrIQQay1ZOfv4KFX1nKsoprbUsZyy8XxdO2i6ZAiItI0JQsRYsv+I9z53Eo+KDrAjJFxPHz9FEYP7OV2WCIi0g4oWejgKqtrWPSulyfyNtK1cyce+spkbpo+nE6dNB1SRERaRslCB7Z8u487nl3J2s8PccWkwdx37SQG9dF0SBERCUzEJwsdcW2IoxVVPPb6Bp5+fzP9e3Xjqa8ncsXkwW6HJSIi7VTEJwvW2heBF5OSkm52O5ZgeGv9Xu56bhU7fcf42jnDSbtyPH26R7kdloiItGMRnyx0FAcOl/PAS2t4/rNdxA+IJvuW85g+Ms7tsEREpANQstDOWWt5tnAnD768hsPlVfzwsjH8zyXxdOvS2e3QRESkg1Cy0I5tO3CUu55fybsb95MwPIaFs6YydlBvt8MSEZEORslCO1RVXcOf3t/ML3M30KVTJx64dhJfO2eEpkOKiEhIKFloZ1btPMgdz65g1c5DJE8YxAPXTWJI3x5uhyUiIh2YkoV24lhFNb/K28Af3ttMXHRXnvxaAldOHowxak0QEZHQUrLQDry3cT93PreSbcVHuXH6MBZcOYG+PTUdUkRE2oaShTBWcqSCB19eyzOFOxjVP5p/3Xwu58X3czssERGJMEoWwpC1lheW7+L+F9dw8Fgl/3NJPP936Ri6R2k6pIiItD0lC2FmR8lR7npuFW9v2Me0YTH8/fopTBjSx+2wREQkgilZCBPVNZY/f7CFx15fD8C9MyfyzfNG0lnTIUVExGVKFsLAml2HWPDsCpbvOMgl4wbwwHWTGRrb0+2wREREACULriqrrOaJZRtZ9I6Xvj2ieOLGs7hm2hmaDikiImFFyYJLPijaz53PrmTLgaPMThzKXVdNIDa6q9thiYiInETJQhs7eLSSh15Zw5L8HYzo15N/fP8cLhjd3+2wREREmhTxyYIxZiYwc/To0SE9j7WWl1d+zs9fWEPJ0QpuuTieH102hh5dNR1SRETCW8QnC9baF4EXk5KSbg7VOXb5jnHP86tYtm4vU87sy5+/M53JZ/YN1elERESCKuKThVCqrrH87cMtPLp0PTUW7rpqAt+5YCRdOndyOzQREZEWU7IQIut3l3LHsyv4dJuPC8f05+GvTGFYnKZDiohI+6NkIQR8Ryv4ypPv061LJx6/YRrXnXWmpkOKiEi7pWQhBGJ6duWXc85i+shY+vXq5nY4IiIip0XJQohcMXmw2yGIiIgEhUbaiYiISLOULIiIiEizlCyIiIhIs5QsiIiISLOULIiIiEizlCyIiIhIs5QsiIiISLOULIiIiEizlCyIiIhIs5QsiIiISLOMtdbtGMKCMWYfsDXIh+0P7A/yMdsrvRf16f2oT+/HCXov6tP7UV+w348R1toBp9pJyUIIGWPyrbVJbscRDvRe1Kf3oz69HyfovahP70d9br0f6oYQERGRZilZEBERkWYpWQitLLcDCCN6L+rT+1Gf3o8T9F7Up/ejPlfeD41ZEBERkWapZUFERESapWRBpI0YYzzGmGS34xARCVQXtwPoiIwxc4EY/914IN1a63UxpLBgjMm21qa6HYeLEoBFxpgYwAfkA2nW2kJ3w3KXMWY+zvtRDGCtzXE3orZnjCnA+f0A572oVWytjXchJNf5P0fB+SztBzxirfU185QOq841xQfEWGsz2joGJQtBZoyZX/c/0hgzG8jFSRoiljEmAZjtdhxus9bGGmNiIvVDryFjTDZOwuT137fGmNgIfH/ygFT8CZOfx3+LOP4EMqvu74H/dyXivmwYY9KBotrrijFmtjEm3Vqb1pZxqBsi+Ob5E4RahYDH/20yksW5HUC4iMALYaP835Y+adDqFh9p74//s2GxtdZrrfXV3oCkSGxl8ZveyO+BN9I+R/2vd7619vgMCP/vxNymnxUaShaCL6XBH7gH8EXaB2BdxpjZ1to8t+OQsJMO1LsYRmJ3nT85qNcV5f/CscSlkMKBx98aWVcktsh5qN8tVau4rcc/qRsiyBr5sEsjApvOavn/4CO6T76uBn/gCTRoao0U/m9MMf5/z8b5QIzY96MRcRH+PtwMFBhj0qy1Gf6/m0y3gwojPk6Mi2sTShZCxP8BmIIzuDGSv1V7IrgptaFCOJFQGmO8QDbO70mkqf3GFFP7+2GMyQeWAYluBuY2f5KVkRYAAAWMSURBVPdMJLcqYK0tNMbE4yQMC4DUSBwI7H8faGScU5uPZ1E3RIhYa3OstfOABGNMRGbE/u4HJQp+/j5pb937NN7cGgnicL4Z1X0/fHBS60skSozwVgWMMR6cAdGjcCoW5taZHRFpbgYW1N7x/33kt3UQShZCzD+CdU6kfQD6/9gjrv+5FXxAJK6o54VGB3sWc2IKYcSp0yUT6dKstRn+8RxpOK1N6ZH2OQrHBzQuNsYk10kUPLRx9666IYLI/w1xmbU2tsFDXpym5kjqjkgGYhr+cdfOqa87ujcS+JOnImutafBQMfWny0UEa63XmIZvxXGRfLG8AfjE7SDc5P/MyK27zd8cn0rkfY4Czuuve98YE0cbfxlTshBccTS+yIcHKGrjWFzVWDLgnxvc5sVEwkQxMK+R7UlNbI8EhcYYT4NBwR5caGINIwnAYreDCFNe4IDbQbQ1f2tTXoNuury2njmkboggamwgY53+6IgesBTpGuuDrh3IFonTBf3S/Dfg+N+KNxIHstXR1FS5iOH/HL2hkYdmE5krUC6gfp2aedT5u2krWnUyyPxTwuoOxIn4cs/+THgezh97DpAZqTNE6pQ2joHjY1oilv9bU+2o7n5tXZUu3BhjinBqtUTs5wUc/xxdgNOSUPv3khOJ74v/89PDibLXmW68D0oWREREpFnqhhAREZFmKVkQERGRZilZEBERkWYpWRAREZFmKVkQERGRZilZEOkgjDFzjTG2kVuuv4KkiEirKFkQ6XhirbXGX1o6EWeeekEjpbcLjDG5jR7BZeEcm0gkUrlnkQ7MXw0x1b/yabYxZlSdapLhvBpqOMcmEnFUlEmkg/CXj87EaVlorLx0CZAV6VUSRSRw6oYQiRxLcEpui4gERMmCSOQo4MQ6DBhjso0x2Q3vG2PmG2OKjDEl/n/H+AdJWv/2hLoHNcak+/ctMcakn+p4DZ4737/d1n28YWx19q09TrZ//QACOV8g/ANG5/t/evw/5/rPoQGjElGULIhEDi/AKS50tS0PKTgr/KXjJBnpnBgsWS/B8D8n1X+b7e8OafJ4tef3D7hMBy7DWXAttTbGhvxjLub5b4k4S35vbpgwNHe+QPif4wXy/DEmW2uz/Euvf4ILq/6JuEnJgkjkqL1oFjezj89am+Ff1e4R/7Y8a22ef7BkZu1x/BfU2UCq//E8nIvovFMcrzaOGJzkw2ut9fqPkdMwIH9Lxlyc1Rjz/PvOw7mYpzfYvbnzBcLjfz1J/vjqLo0cT4QvIy2RR7MhRCJHPM7FtLkLXX7tP6y1PmMMOC0LteomGrVTMQv8+7X0eLWtAXn+45UYY/KA3CaW7E7yx92w1SGPk8dgNHe+FquzhHoKsLjBw7VLrotEDLUsiESOuTiDHAN1qpYI0/DWkoNaa33/v507umkYhqIwfLxBYAS6QaRu0BFSsQEjIHUERkhXaDYANkBswAywweXBN20piZugAGr8f1Jf2kSOn3Lqe20za8sPH4olg92Z2/7aSjGUSNqvptwchQkgC4QFIAPeeFho2lr7k6TitOFxLDNrzGzf89BxyUvPOF9e5FPzYFB4+aVVSWr8d3aWIBuEBWDGQgil/1u/V6z5T1Zr97LAVvGwp5XvGKiGnrzo19Z+bynpVh0Njv6ybnyc6mhOpUaGH3/GoTskusLIUtKjN1ZejxkbuGSEBWB+2q2IpthvUEha/MbSuTca1v55k7TR8NMXXxVfuDsdtnWue8ZZKwaGB0nPOsxpbPhpd2AMsdD3udT+/d1J0yMwa5zgCCArIYR3M7v67+cALgkrCwCy4eWOnzR5AlljZQFANnqaFgGcQVgAAABJlCEAAEASYQEAACQRFgAAQBJhAQAAJBEWAABAEmEBAAAkERYAAEDSJ7ZrsIfSbcZ8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1147f1048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 5))\n",
    "plt.semilogy(n_list, times_simplex, label=\"Simplex\")\n",
    "plt.semilogy(n_list, times_barrier_simple, label=\"Primal barrier\")\n",
    "plt.legend(fontsize=18)\n",
    "plt.xlabel(\"Dimension, $n$\", fontsize=18)\n",
    "plt.ylabel(\"Computation time, sec.\", fontsize=18)\n",
    "plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Комментарии\n",
    "\n",
    "- Было показано, что прямой метод эквивалентен методу Кармаркара\n",
    "- Использует информацию только о прямой задаче\n",
    "- Начальное приближение должно лежать в допустимом множестве - отдельная задача"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Сравнение способов решения линейной системы в методе Ньютона\n",
    "\n",
    "- Явно сформированная матрица\n",
    "- Использование блочной структуры и последовательное исключение переменных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def generate_linprog_problem(n, m=10):\n",
    "    x = np.random.rand(n)\n",
    "    A = np.random.randn(m, n)\n",
    "    b = A.dot(x)\n",
    "    c = np.random.randn(n)\n",
    "    c[c < 0] = 1\n",
    "    return c, A, b, x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "n = 20\n",
    "c, A, b, x0 = generate_linprog_problem(n)\n",
    "res = scopt.linprog(c, A_eq=A, b_eq=b, bounds=[(0, None) for i in range(n)])\n",
    "m = A.shape[0]\n",
    "# x0 = np.random.rand(n)\n",
    "# resid = b - A[:, m:].dot(x0[m:])\n",
    "# x0[:m] = np.linalg.solve(A[:m, :m], resid)\n",
    "print(np.linalg.norm(A.dot(x0) - b))\n",
    "mu0 = 1\n",
    "rho_mu = 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.454027024366833e-08\n",
      "8.297346244593269e-08\n",
      "3.4614053214801457\n",
      "3.461405402686456\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "x_bar_simple = BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, backtracking, \n",
    "                                      simple_solver, max_iter=100, tol=1e-8)\n",
    "x_bar_elim = BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, backtracking, \n",
    "                                    elimination_solver, tol=1e-8)\n",
    "print(np.linalg.norm(res[\"x\"] - x_bar_simple))\n",
    "print(np.linalg.norm(res[\"x\"] - x_bar_elim))\n",
    "print(c.dot(res[\"x\"]))\n",
    "print(c.dot(x_bar_elim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current dimension = 20\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 ms ± 569 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14.6 ms ± 325 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "Current dimension = 120\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "31.4 ms ± 392 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25.4 ms ± 2.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Current dimension = 220\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65.7 ms ± 4.73 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32.7 ms ± 2.05 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Current dimension = 320\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "147 ms ± 11.7 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65.5 ms ± 4.7 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Current dimension = 420\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "184 ms ± 16.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.18 s ± 80.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "Current dimension = 520\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "290 ms ± 8.25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "87.8 ms ± 7.53 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Current dimension = 620\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "719 ms ± 20.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/cvxpy/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "275 ms ± 55 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "n_list = [10*i for i in range(2, 63, 10)]\n",
    "times_simple = np.zeros(len(n_list))\n",
    "times_elim = np.zeros(len(n_list))\n",
    "mu0 = 5\n",
    "rho_mu = 0.5\n",
    "for i, n in enumerate(n_list):\n",
    "    print(\"Current dimension = {}\".format(n))\n",
    "    c, A, b, x0 = generate_linprog_problem(n)\n",
    "    time = %timeit -o BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, backtracking, simple_solver)\n",
    "    times_simple[i] = time.average\n",
    "    time = %timeit -o BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, backtracking, elimination_solver)\n",
    "    times_elim[i] = time.average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEWCAYAAACwtjr+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XlcVXX++PHXRxFx5YKa+3bRyiUrwExbrISWqaapQJsWtamgZubbLL9Jsplp9gyafSuwTbMmBJv2DXQqSzOBXEpL5QIuuePFBZHt8/vjnIsXvCwX7r3nXng/Hw8exb2Xc95H8b7vZznvt9JaI4QQQvhTN6sDEEII0flJshFCCOF3kmyEEEL4nSQbIYQQfifJRgghhN9JshFCCOF3kmyEEEL4nSQbIYQQfifJRgghhN+FWR1AsBg4cKAeM2aM1WEIIURIKSwsPKS1HtTa6yTZmMaMGUNBQYHVYQghREhRSpW15XUyjSaEEMLvJNkIIYTwO0k2Qggh/E6SjRBCCL+TZCOEEMLvJNkIIYTwu0659VkplQDYXN9rrXMtDEcIIbq8zjqySdVa55pJZo5SytbqTwgh/O6jbQepqqmzOgxhgU6XbMxRTbnbQ+uB2RaFI4QwrS0+zNxnP+Mfq7ZbHYqwQKdLNoAdcLp97wRiLIpFCGFaurbU/G8Zx6pqLI1FBF5QJxulVKxSKqeZ51KUUknm1wK3p2zA4SYvl2k0ISz0jfMk72/Zz+VnD+JYVS0vrttpdUgiwIIy2ZhJJh2YgzFSafp8ChgL/+a6TK5SKtN82gkMaPIjToQQlnlxXRn1WvOH70zmsvEDeebjElm76WKCMtlorYu01mlAdjMvSdVaZ7m93gEkmN86aDySsWGs2wghLFBVU8d/PtvFrHMHMzK6Nw/MjOHgsVOsKNptdWgigIIy2bTE3FkW6+Epp1IqQWudD0S7PR4D5AckOCHEGd7evJfyE9XMnzEGgOkxAzh/pI3MDx3U1tVbG5wImJBLNpy5AcClnNNJaJG5lpMA5GmtZRpNCIssWVOKfVAfLhlnzG4rpXhgZgw7yyt5+4t9FkcnAiUUk000jbc2uzSs1ZjTcLla6/yWbug0NxkUKKUKDh486Kdwhei6NuxysnF3BfOmj0Ep1fD41RMHEzOoD09+UIzW2sIIRaCEYrLxGa11ltY6XmsdP2hQq43mhBBeWrqmlL49w7g1bkSjx7t1U9w/M4ate4/ywTb5oNcVhGqyifbwmKctz0IIixw6foo3N+3l1tjh9O15ZmWsmy4YzrDICJ78oNiC6ESghWKyKcDzfTPRQFGAYxFCNOPlz3ZSXVfPXdPHeHw+PKwb915m57OScgrLPM2Mi84k5JKNudjv8FDvzGbuRBNCWKy2rp5ln+7k0nEDGXdW32Zfd9tFI4nq3UNGN11AsCcbT9NlAOnAQtc3SqlY2rm9WSl1o1Iqq6Kioj0/LoTwIG/LfvYdrWKeud25Ob3Dw5g/Yyz5Ww/w9b5jgQlOWCIok41Sym5WEEgHYpVSma6qAWAs7APFSqkEpVQSkKC1Tm3PubTWb2itUyIjI30TvBCCJWtLGW7rxVXnntXqa+fNGE3v8O489aGMbjqzoOxnY1YESGvlNVktPS+EsMZX+47yqaOch687l+7dVKuvt/UO5/aLRvHcmlJ+mng2I6N7ByBKEWhBObIRQoSupWvL6BnWjTnxI9v8M/deZqebgsWrHX6MTFhJko0QwmcqTtbw36I93HTBMKL6hLf554ZERnDLhSPIXr+LQ8dP+TFCYZUun2xkg4AQvpNbuJuTNXXMbWa7c0tSZ9qprqvnuU9KfB+YsFyXTzayQUAI36iv17ywtpS40VFMHu79vyf7oL5cN3mINFfrpLp8shFC+MaH2w9SeriSudNHt/sYD8wcJ83VrFBX6/dTSLIRQvjE0jWlDOrXk+smD233Mc4bESnN1QKp9hS88zDkzAM/F0T1SbJRSj2klJLfDCG6qNJDJ/hg20Fuv2gU4WEde1uR5moBcrgYnkmEdU9C5Eio9+9buK9GNrnA/T46VkDJBgEhOu6FT8vorhS3TxvV4WNJc7UA2JwLmTPhSBnc9hJc9zh09+9tlz5JNlrrEq31Yl8cK9Bkg4AQHVNZXcvygl1cO3kIg/tHdPh40lzNj2pOwhs/ghX3wFkT4P6P4dzrA3JqWbMRQnTIq59/w7Gq2oa2z74gzdX84ODXsPgqKHweLv0J3P022Np+421HeZVslFL3KqWe9PB4gVLqHt+FJYQIBVprlq4tZeLQ/sSNjvLZcaW5mo9teAmyroDjB+COFZDwa+jeI6AheDuySQPyPDz+OCG6ZiOEaL91JeV8te8Y82aMbtT22RekuZoPnDoOr6TCqw/A8Dhj2mx8giWheJtsYgBPxYuKgdiOhyOECCVL15Zi692Dmy4Y7vNjS3O1Dtq32RjNbF4OVzwCc1+D/u3flt5R3iYbBzDWw+N2ICRrTMhuNCHaZ2/FSd77cj9z4kcS0aO7X84hzdXaQWtY/wwsngWnjsHc1+GKNOjmn7+jtvI22WQBTyulrnQ9oJS6ynz8KV8GFiiyG02I9nlp3U7qtebOi9tfMaA10lzNS1UVkDMf3vopjL3MmDYbe5nVUQFeJhutdQbwNLBSKVVn3siZByzWWv/RHwEKIYLPqdo6/vPZTmade5bf+89Ic7U22lMIT10GW9+AhN/A7TnQd5DVUTXweuuz1joNiAJmm1/RWuuHfR2YECJ4vb15L4eOV7fa9tkXXM3VXt/4DbvKK/1+vpCjNaz9FzxzDeh6uPsduPTH0C247mzxOhql1L3ASiBba71Ca12hlHpPKXWz78MTQgSjJWvKsA/qwyUxAwNyPmmu1ozKcvjPd+G9R2D81ZD6EYyaZnVUHnl7n83jGNufFwHu+xwXA4/4MC4hRJDauMvJhl1O5l48mm5taPvsC9JczYOytfDUpVC8Eq7LgNtehN7RVkfVLG9HNilAqtZ6RZPHC5Gtz0J0CUvWltInvDu3xo0I6HlTpLmaob4eVv8Jnr8euofDPe/DtFTw8X1OvtaeST1X7Qj3KwvZrc9CiLY7fPwUb27cy61xI+gXEdg70GOkuZpRAWDZLbDytzDxJmPabNiFVkfVJt4mm8VAllLqAsyko5S6EGPbc0hufZb7bIRou5fX76K6rr5DDdI6oks3V3N8YEyb7VwLN/4Nkp6FiP5WR9Vm3m59TgM2AEWAUkodxphCWxmqW5/lPhsh2qa2rp4XPy3jknEDGHdWP0ti6JLN1epqYdUfYOl3ICIS7lsFcfODftqsqfZsfU4GxgGpwMOAXWstddGE6OTyt+7nm4oq5k4fY2kcXaq52tFvYOm34aMMuOB2SPkABk+yOqp28XY32s+UUrdorR1m/5pEoFgptc2cWhNCdFJL1pQx3NaLhAmDLY2jyzRX2/a+MW32zQa4ORO+828I72N1VO3m7chmIeAEoxU0cCvGjZ0bMdZzhBCd0Lb9x1jrOMydF4+me4C2Ozen0zdXq6uB938BLyVDv6GQ+iGcf5vVUXWYt31Aozhd9TkRyNVar1BKlQDrfRqZECJoLFlTSnhYN+ZMDVyzrZa4N1e7ccpQn7c3sMyRMsj9HuwpgPh74Jo/QI9eVkflE96ObJxApFIqEkgAss3Ho8znhBCdTMXJGl4p2sO3zx9GdJ9wq8MBOmlztS2vQ+ZlcGgbJD8PN/y50yQa8D7ZPA58jjG6cWitXzEfT8YoYSOE6GRWFO7mZE2dT9s++0Knaa5WUwVvPwTL74Jou3HvzKTOV/3Lq2k0rXWGUqoIiATy3Z7Kw9gOLYToROrrNS98WkbsKBuThwfX7QGu5mq/fXMLhWXlxI0O3lItzTpcbLQE2LcJLv6B0a45LDhGj77Wnq3P+a4CnG6PrdBah2QFAbmpU4jmfbT9ICWHTgSkunN7hHRztU05kHk5VOyC774M1z7WaRMNtK9cTaciN3UK0byla8sY2Lcn1022rp1wS0KyuVp1Jbz2Q3jlXhg82Whwds51Vkfld10+2QghPCs7fIL/fX2A26eNIjwseN8qQqq52oGtsPgq+HwZXPb/YP5bEBnYgqZWCd7fICGEpZZ9WkZ3pbhj2iirQ2lRSDRX0xqKXoCsK6HyENy5AmY9Ct29vfskdEmyEUKcobK6luz1u7hm8hAG94+wOpxW3XPZ2OBtrnbqGLySAq//EEZONabNxs2yOqqAk2QjhDjDaxu+4WhVbdBtd27O0Mhewdlcbe8myJwJX+TClT+Hu16FfkOsjsoSkmyEEI1orVmyppQJQ/sTPzrK6nDaLKiaq2kNny2GpxOgphLmvQEzF0C37lZHZhlJNkKIRtaXHuGrfceYN310SJWBCZrmaiedxg2ab/8Mxl5uTJuNudS6eIKET5KNUuohpVSdUuqwUupJqQAtROhasqaUyF49uOmC4VaH4jXLm6vtLjBKznz9DiT+Dm5fDn0GWhNLkPHVyCYXuF9rPQDIAuQOSSFC0L6KKt79ch9zpo6kV3joTflY1lytvh7W/AOevcboYXz3u3DJg9BNJo9cfPInobUuMfvboLX+PFSrCQjR1b20rox6rblzmjVtn30h4M3VThyG/9xmtAU4+1q4/yNj15loxOtko5S6VylVoJSqdXvsPaVU56scJ0QXcqq2jpc+28lV55zFqAG9rQ6n3QLaXK1sjdHgzPE/uO4JmLMMeoXOpopA8rZT5+NAGrCoyc8uBh7xYVwBI7XRhDC8s3kfh45XMzdEtjs3JyDN1epq4aMn4PnroUcE3JMH01IghDZUBJq3I5sUIFVrvQJjZtKlEIj1WVQBJLXRhDAsWVuKfWAfLhsX+gva7s3VtNat/0Bb7d8C7/0c/jwBVv0eJt0CKR/CMNkT1Zr21Epw/c25p3A7IOs0QoSoTbudfL7Tya9unEg3i9s++4KrudpDuZv4YNtBrjznrPYfrLIcNufChhdh7wboFmaszVx4F5x9TciPZmrr6qmsqaN/RA+/nsfbZLMYyFJKJWMmHaXUhcBT5pcQIgQtXVtG7/Du3BrXeYpC3nTBcP6St40nPyj2PtnU1cCOfNjwkrGNub4GhkyBa9PhvKROs535VG0dP/rPBvYerSL3/un06O6/3XPeNk9LU0rZMRulKaUOAzYgS2v9Rz/EJ4Tws8PHT/H6xm+YHT/C759uA6ldzdX2f2kkmE3ZcOIg9B4IF6XABd+FIef5P+gAOnGqltQXCvl4xyF+deNEvyYaaMc0mtY62Uw4rkpy+bLVWYjQlV2wi+raeuZNH2N1KD5320Uj+ceq7Tz5QTFPz2sm2Zw4bNQu2/Ai7N0I3XrAOdfCBXfAuATo3nkSsIuzspq7n1/Ppt0V/Cn5/ICMaNtV31pr7QCCsLyqEMIbtXX1vPjpTmbEDGD84H5Wh+NzruZqf8nfxtf7jnHOEPMa62pge56RYLa9Z0yTDT0frsuAyUnQZ4C1gfvRgaNV3PXMZ5QcOsG/74jlmkmBKQzqVbIx12eWY2wIOIPWOvRuORaiC1v51QH2OE/yyxsmWh2K38ybMZrMj4p56sNi/jIz7PQ0WeUh6DMIpqXC+d+FIZOtDtXvdpVXcucz6zh47BTP3T2VSwK489DbkU0O4ATuB8p9H44QIpCWrCllWGQECRM6sFsryNn0Uf4yei3Dv3wVtpaa02TXmdNkszrlNJkn2/Yf465n1lFVU8+L907jwlGBvfnU22RjB2K11hv8EYwQInC27z/GmuLDLLj2HML8vDgccHU1sP19YxSz7V2uqa9lM2N5c8RPuOH2B6F3GzYLdCIbdzmZ99xnhHfvxvLU6aenEwPI22STD0gtBiE6gaVrywgP68ac+JFWh+I7+zab02TLzWmys2Da/XDB7SxbXcerG/ZwcX0fOsfG5bZZU3yI+5YUEN03nBfvudiyUkTeJptkoEAplYex/bnRVJrW+hVfBSaE8J+jVTWsKNrNjVOGMaBvT6vD6ZjjB2FzjpFk9m+G7uHGNNn5tzeaJkuZeZzlhbt47pMSHrrmXIuDDoy8Lfv5wUtFjBnQmxfumWZpi29vk819QIz51ZQGZIOAECFgReFuKqvrQqbt8xlqq09Pk21/D+prYdiF8K0/wuRbPU6TuTdXu39mDP060T1Fnvz38938LGcTk4dH8vz8qUT1Cbc0Hm+TzSNAOvC41loqVwoRgurrNS+sLePCUTbOGxFiNQH3bjISzOblUHnYmCa7+AFjFDO49R11D8wcx9ub9/Hiup3cP9PTZ+bOYcmaUn71+pfMiBlA1tx4+vZs110uPuVtBOVApiQaIULXxzsO4Th0gr/OCZHikccPGsllw0uw/wtzmuxbxm6ymKuge9vfxtybq82fMYaIHp1rMkZrzT9W7eDPedu4euJg/v7dC4PmGr1NNplAEiClaYQIUUvXljKwbzjXnReYm/napbbamB7b8JIxXVZfC8NiW5wma6sHZsZw+9PrWFG0mztCuElcU1prfv/WVp75uIRbYoeTceuUoNpl6G2ymQrcqpRaCBQ0fVJrfY1PohJC+MXOw5Ws/OoA/3flOHqGBccn3gZaG+ViNrxkLPifLIe+g2H6D4xpsrN8s6jv3lxtTvzIoHpDbq/aunoWvrKZnMLdzJ8xhkdvCL7q3e2ZyFvh8ygspJS6Ebhx3LhxVocihN8tW1dGN6W4PZg+0R8/YGxV3vASHPjSmCY793pjmsx+pVfTZG3haq52/7JC3v5iH98+f5hPjx9orsrN7365jx/NGs+PE8ajgrDtgbdVn2f7KxCraK3fAN6Ij4+/z+pYhPCnk9V1ZK/fxbWThjAk0rotsADUnoJt78KG/xjTZLoOhsfB9X8yGpL5+aZL9+ZqN04ZGpRvzm1x4lQt9y8rZPX2Qzx6w0S+d+lYq0NqlvVbFIQQAfH6xj1UnKxh7nSLRjW11cY02ebl5jTZEeg7BGb80KfTZG3h0+ZqFqmorGH+85+xcZeTJ5KmkBzkN+e2mGyUUu8BOVrrp83vs1t6vdZ6jg9jE0L4iNaa59eUce6Qflw0NgClWo7tN26w3P8l7PvC+O+hr42F/u493abJrvD5NFlbdai5msUOHKti7jOf4Th4gn/fEce1k4N4s4eptb/lpnW2Q3OsKUQXV1B2hK17j7LolvN8O2VUW20kkX1fGNuS95uJ5cTB06/pPxwGTzJaKA+eZNzV38v6qlftaq4WBNwrNz87fyqXjg+N4jstJhutdXyT7zvdmo0QXcGSNaX0jwjjpgs6sBje0mgFjBHLWRNg/DVGuf7Bk2Dw5KAuetmm5mpBZPv+Y9xpVm5edu80YgNcubkjvO1nc4un+mdKqasAtNarfBWYEMI39h+t4t0v9jF/xhh6h7fhn3x7RitDzoPoGMumxNqr2eZqQWjjLifzn/uMsO7dyE69mHOH9Lc6JK+0p5+Np835Cngc4z4cIUQQeXHdTuq05i5PGwM64WjFW42aqwVpVYW1xYe5d8l6ovuGs+yeaYwe0MfqkLzmbbJpbrLXBsR2MBYhhI9V19bz0rqdJJwdxegaB2zo/KMVb9l6h3P7RaN4bk0pP008m5HR1pTgb07+lv18/6UiRkf3Ztm91lZu7og2/RYppeoxqjprpVSdp5cAhb4MTAjRTm6jlf1bP+OF6o2cs3MvPNV1RiveuueysSxZW8ri1Q5+e1PwtIduqNw8rD/P332R5ZWbO6KtH1kSMRLK+xg9bdw5AYfWusSXgQkhWtGGtZVe3Qayp8do1MXJMKTrjFa8NTSyF7dcOILs9bt4cNZ4BgZBj5+la0t59LUvmW4fwOJ5wVG5uSPaFL3WeiWAUipfa92pytUIERLasbaytX4U1y3+kkdvmMjFQXxnebBImWkPiuZqWmv+uWoHf8rbRuLEwfwjiCo3d4S35Wqu9lcgQggPdhfCGw8aoxaXNq6tPJuzkd7h3bk1bkSAgw5NwdBcTWvNH97aytMfl3DLhcPJSAquys0dEdrjMiE6q5qT8L/HYO0/od9QuPoPMHRKm9dWjpyo5rWN35AcN4LIXp27I6UvWdlcra5es/CVTSwvCN7KzR0hyUaIYLPzU3jtB3B4B8TNh8TfQoR3HTWzC3ZRXVvPvFBt+2wRq5qrnaqt48cvb+CdL/bx4Kzx/CRIKzd3ROcYnwnRGVSfgHfS4Nlroa4a5r4GN/7N60RTZ7Z9nm4fwNmDg/cmxWD1wMwYDh47xYqi3QE5X2V1LfcuKeCdL/bxyxsm8tPEsztdogFJNkIEh5KP4N/TYd1TcFEKPLDWKFLZDiu37meP8yTzZgRRz5oQ4t5crbau3q/nqqis4c6n1/HJjkNkJE3hnk68kUOSjRBWqjoKb/4EltwI3brD3e/AtzKgZ992H3Lp2jKGRUaQMGGwDwPtOlzN1XaWV/L2F/v8dp4Dx6qYk7WWL/Yc5d93xDI7yFsEdJQkGyGssiPfGM0UPg/Tfwj3fwKjZ3TskAeO8fGOQ9xx8ehOs4vJCu7N1bTWPj/+rvJKZj+1lp3llTw7fyrXTh7q83MEG683CCilxgB2jBI1jXgq0mkVpZQdSNVap1kdixCNnDwC7/0CNiyDgefA996Hkb4pK7h0bRnh3btx29TO/SnZ3/zZXG3HgWPc+fRnVFbXhlzl5o7wturzQxgFNz2tXmk8F+kMOKVUEkbVAyGCy9fvwBs/Nu7yv+z/weULoIdval0dq6phReFubjh/KAOC4A74UOeP5mqbdjuZ96xRuXn5/dNDrnJzR3g7zl4IPAFEaa27NfkKikQDoLXOxahQLURwqCyHFffCf26DPgPhvpUw61GfJRqAV4r2cKK6jnnTx/jsmF2Zq7naZyXlFJaVd/h4a4sPc/vidfTpGUZOatdKNOB9sikHntJaV/gjGCE6pS9fhX9dBF/+F65YCPf9D4Zd6NNTaK1ZsraUC0baOH/kGTPcop1uu2gkUb178OQHxR06Tv6W/cx77jOGRkaQe/8MxgwMvRYBHeVtsskEkjp6UqVUrFLK48hDKZWilEoyvxZ09FxCWOb4AVg+F3LmGSVmUj6EKx6GMN9X7v14xyEcB0/IdmcfczVXy996gK/3HWvXMV79fA+pywqZMKQfy1OnMyQyNFsEdJS3GwSmArcqpRYCBU2f1Fpf09IPK6VigTnmt3YPz6eYx8k1v7crpTK11qnm9wuAAc0cPk9rnd/WCxHCb7SGzbnwzgKoPm5Ml834kV8rLS9ZU8aAPuF867zOv6sp0DrSXO2FtaU8+vqXXDy2c1Ru7oj2XHm7qz5rrYuAIjPpJHh4SarWOs7t9Q6lVILb9xntPbcQAXF0r3HfzLZ3YMRUuOlfMOgcv55yV3klK7/azw+uGEfPsKBZOu002tNcTWvNv/63gz++v42ECYP55+2do3JzR3hb9Xm2vwJRSjXX7dOplEqQUYsIalrDhhfh3Ueg7pRROPPiB4wbNf1s2adldFOKOy4e5fdzdVXeNFfTWvPY21tZvLqEm83KzT3knifvb+pUSvVXSj2plNpufv1bKeWLAkx2jEZsTZXjZctpczSUCiS4puaE8BvnLlh2q1E8c/AkeGANzPhhQBLNyeo6Xl6/i2smDWZoZC+/n6+rcm+uduj4qWZfV1eveXjFZhavLmHe9NH8Kfl8STQmb++ziQRKAAenp9MSMEYfdq11WQdiicZILE05aX6dxiNzFCQjIeFf9fVQ9Dy8/yjoerjuCZh6L3QL3JvLGxu/oeJkDXNlu7PftdZc7VRtHT/J3sDbm/fx4FXj+EknLajZXt6u2WQBy7XW97s/aO4sywSu9VVggWCOelIARo2SKQjhhfISeP3/oHQ1jJ0J3/47RI0JaAhaa55fU8o5g/sxbWzrPW5Ex7TUXK2yupbUFwpZvf0Qv7h+Avdedsb+py7P249giRgVBJp6DGOnWkd5+hdjAw774Nhn0Fpnaa3jtdbxgwYN8scpRGdTXw+fPgVPzoBvNhgtAOa+FvBEA1BYdoQte48yb8YY+QQdIA/MHMexqlpeXLez4bGKyhrueuYzo3LzrVMk0TTD25FNOTAWKG3yuB3PU2DeKMBDvTWMBFTUwWML0XGHdsDrP4Sda2FcItz4V4i0ruXykrVl9IsI4zsXDrMshq6maXO1Y1W13PXMOooPHudft8dynWw9b5a3I5ssIEsp1XDnmFLqAozRTm5HAtFaOwGHuSvNnU12oglL1dfBJ3+Hpy6BA1vgO0/CHTmWJpoDR6t4Z/NeZsePpHd41713wwqu5mr/XLWD5KfWUHbYqNwsiaZl3m59zlBKxQAlSqkj5sM2IFdrvdCLQzU3wZyOUX8tDRpuAvVrolFK3QjcOG7cOH+eRoSqA1/Ba9+HPYVwzvVww5+h3xCro+Klz3ZSpzV3XSwVAwLN1Vztn//bQf+IMJbdO4240V2jcnNHqPb0ajCTgOvmywKt9edt/Dk75pZkjO3MWUCh1jrL7TUpGLvdbIA9UDdyxsfH64KCM4oiiK6qrgY++Rt8mA7hfeFbT8DkWyEI1kaqa+u5JH0Vk4f157m7L7I6nC5pbfFh0t/9ikW3nMeEoV2roGZTSqlCrXV8a69r1/jbVQmgHT/nwBy1tPCarJaeF8Lv9m2GV78P+zbBpJuNLc19g2cDybtf7uPgsVPMnTHG6lC6rOkxA3j1B5dYHUZIaTHZKKXeA3K01k+b32e39Hqt9ZyWnhciqNVWw+o/wuo/Qa9omP0CTPy21VGdYemaUsYM6M3M8cGTAIVoTWsjm6Y3U1o/hyCEP+wpMioAHNgCU26DaxdB7+C7d+WLPRUUlB3hlzdMpFs3+ecoQkeLyabpPJw/a6NZRTYIdHE1VfDBIljzd+g7GL6bDecE773JS9eW0qtHd5LirNsJJ0R7eLX1WSl1SzOPX6WUuso3IQWW1voNrXVKZGSk1aGIQNu5Dp66FD75K1xwB3z/06BONEdOVPPahm+4OXY4kb16tP4DQgQRbzcI5ACeqgsqjHttfFFFQAj/qq6EVb+HT/9t3Ctz138hJvg/Ky0v2MWp2nrmTpftziL0eJtsmpskbq49gBDBpfRjeO2HcKTEKJqZ8Gvo6Yui5f5VV6954dMyLrZHd7ne9aJX8yK/AAAgAElEQVRzaFOyUUrVAxrQSqk6Ty8BCn0ZmBA+deoY5P8a1j9t1DGb9yaMvczqqNps1VcH2H3kJD//1gSrQxGiXdo6sknESCjvA8lNnnMCDq11iS8DE8JnilfB6z+Cil1w8ffhql9AeB+ro/LK0rWlDI2MIHHiYKtDEaJd2pRstNYrAZRS+VrrdreFDkayG60Tq6qA938BRUthwHj43nswaprVUXltx4HjrN5+iJ9dfTZh0ohLhCivfnO11lf7KxCryG60Tuqrt+FfF8Pny+CSH8P9H4dkotFa88zHDsK7d+O2i6TnkghdXperMbc/z8FDOwCt9TW+CEqIdisvgXcfhm3vwlkT4bZlMDyu9Z8LQl/sqeC3b27hs5JyvnvRKAb27Wl1SEK0m7dtoR/CqMyci7GO42orkATk+TY0IbxQU2UUzvz4z9AtDK7+PUy7H7qH3v0oB45W8cR7X5NbtJvo3uH84ebJzIkfaXVYQnSItyObh4EErfUqpVS5q6KAUioJmOXz6IRoi+358PbPjO3Mk26Gax6D/qHXUKyqpo6nVzv49wfF1NZpUi6z84OrxtE/IvQSphBNeZtsooBi8/8dSqnztdYbMbY9ZwIP+DI4IVrk3AXvLYStb8CAcXDXqxBzpdVReU1rzRub9pL+zlfscZ7k2klDWPitcxk9ILR2zAnREm+TTREQA5QBy4EMpVQqxtSap5bOQU92o4Wg2mpY+0/46AnQGmY9CtN/CGGht6axYZeT3725hcKyI0wa1p8/zT6fi+1N698KEfq8ap6mlJoFjNVaP62UisRIPmMw7sFJ01o/4ZcoA0Cap4UIxwfw1s/g8HY49wajOrMt9HZp7a04Sca7X/Pfz/cwqF9PHrr6HG6NG0F3qeQsQoxfmqe57rcx/78CiFFKXdjWTp1CtNvRb+C9n8OXr0DUWLg9B84OvZ34ldW1ZH7oIPOjYuo1/ODKGB64Yhx9e7arj6EQIcPb3Wi3aK1fcX9Ma/25q+Kz1nqVL4MTgroaWJdptAGoq4ErFhr3zfSIsDoyr9TXa17dsIeMd79m39EqbpgylLRrz2VkdG+rQxMiIKTqswhepZ8Yu8wObIHxV8N16RBttzoqrxWUlvO7N7ewcXcF54+I5J+3X0j8mOBrzCaEP0nVZxF8jh+A938Jm16GyFFw20twzrdAhdZ6xq7ySh5/9yve2rSXIf0j+Muc87np/OHSYVN0SVL1WQSPulooeNboNVNTCZf9P7jsZxAeWlNNx0/V8uQHO1i8uoRuCn40azypM+30Dpd1GdF1dfmqz7L1OUjsWg9v/QT2bQb7FfCtP8LA8VZH5ZW6es2Kwt088f7XHDx2ipsvHM6Ca89haGQvq0MTwnLebn1+vzMW4wTZ+myZE4ch/1fw+QvQbxhc+xhM/E7ITZmtLT7M797cwpa9R4kdZePRGydxwciQvPVMCK/4a+tzp0w0wgL19VC0BFb+xmhsNuP/YGZaSHTNdFd2+ASPvb2V977cz3BbL/7+3Qu5ccpQVIglSyH8zeutzy0933RbtBAeffM5vPlT+KYIRl8K1/8RzgqtDpRHq2r416odPPdJKWHdFT+7+mzuvcxORA9PmzWFEN6uWOY287hrLk7+pYnmnTwCK39nbALoMwhuWQznJYfUlFltXT3ZBbv48/vbKK+sJil2BA9dcw5n9Q+t+36ECDRvp9HOaLamlLIBWcBjvgpKdDL19bDxJch71Eg401LhykcgIrQa1q3efpDfv7mVr/cf46Kx0Sy5YSKTh4fWNQhhlQ7vxdRaO5VSacBTgDRPE43t22zUMtv1KYy4CK7/EwydYnVUXik+eJzH3trKyq8OMDK6F0/eEcu1k4fIuowQXvDlxv8EHx5LhLqqCvjfIvgsC3rZ4KZ/wfm3QzevOpFbqqKyhr+t3M7StaVE9OjOw9edy/wZY2RdRoh28HaDwHvNPBWPcb+N6Oq0hs058P4vjEoA8XfDVb+E3qFTnqWmrp6X1u3kL/nbOHqyhjlTR/HTxLMZ1C/0WhgIESy8HdlUNPN4DkbzNNGVHfjKqGVWuhqGXQjf/Q8Mj7M6Kq/87+sD/P7NLRQfPMGMmAH84vqJTBzW3+qwhAh53m4QmO2vQKwiFQR84NRx+DAdPv03hPeFG/4CsfOgW+hMN23bf4zfv7WVj7YdZOzAPiyeG0/ChLNkXUYIH/GqgkBnJhUE2kFr2PIavPcIHN0DF94JCb+BPgOtjqzNyk9U85e8bbz02U56h3fnR7PGM3f6GMLDQmdtSQgr+aWCgHngnwGpgKvWexHwmNb6v94eS4SwQzvgnYegeBUMOQ+Sn4eRF1kdVZtV19azdG0pf1u5ncrqOu6YNoofJ5xNdJ/wFn+uoqKCQ4cOUV1dHZhAhbBAeHg4AwcOJDLSd1v7vd0g8BRwH7ACyDAfng3kKqVStNbP+CwyEZyqK2H1n2DN3yEsAq7LgPh7oHtoVDTWWpO3ZT+L3vmKkkMnuPzsQfzy+gmMH9x6mZyqqir279/PiBEj6NWrl0yxiU5Ja83JkyfZvXs3PXv2JCLCNzcse/sOkQKka60Xuj22WCmVCTwMSLLprLSGr9+Gdx6Gip0w5TZI/C30G2x1ZG22de9RfvfmFtYUHyZmUB+eu3sqV55zVpt//uDBgwwaNIjevUOr5YEQ3lBK0bt3bwYOHMjBgwcZOXKkT47rbbJxANkeHn+cM1sPiM6ivATeSYPt78GgCTD/bRhzidVRtdnBY6f4c97XZK/fRf9ePfjNtydx+7RR9Oju3bpMVVUVQ4YM8VOUQgSXfv36cfjwYZ8dz9tkk4lx8+aGJo/PwnMSEqGspgo++Rt8/GfoFgZX/x6m3Q/de1gdWZucqq3juU9K+eeqHVTV1DF/xlgenDUOW++W12WaU1tbS1hYaEwXCtFRYWFh1NbW+u54Xr4+EZillEps8ngCUNT0pk+ttZSvCVXb8+Dth+BICUy6Ba75A/QfZnVUrao4WcPm3RVs3O3k5fU72VV+koQJZ7HwWxOIGdS3w8eXdRrRVfj6d93bZOPE2BzQlKfHRChy7oJ3H4av3oQB4+GuVyHmSquj8qiyupYvvznKxl1ONu2uYPOeCkoOnWh4fvLw/jx2z3lcNn6QhVEKIUBu6hTutrwG/73f2Aww61GY/kMIC44SLdW19Xy17ygbd1ewyUwu2w8co968TWxoZATnDY8kKW4EU0ZEct7wyHZPl3VlWVlZpKamNvv8kSNHsNlsJCcbS7Q5OTmBCq3L6Wx/xjIBLQxFL8AbD8LweEh6BmyjLAulrl6z48BxNu52smm3kVi+2nuM6rp6AKJ692DKCBvXTBrMlBE2poyIlH4yPlZYWIjNdmZba0+PAcTFxREdHU1eXl67z+mLYwTjuYTB2/tsLgSWc/qGzka01qFTn0SctuYfRuHMmFkw5wUI7xOwU2utKTtcaSaWCjbtdvLFnqOcrKkDoG/PMM4bHsndl4xpSCwjouQeF3+z2+3NJhZPWhoNBfIYwXguYfB2ZJODsW5zP1Du+3ACr0vXRtMaVv3OuElz0s1wcxaE+W/qSWvN3oqqhqTi+u/RKmPHS8+wbkwa1p85U0cyZUQkU0bYsA/sQ7dukliCXUpKSlAco6n8/HwSExNpWpbLH+cSLfM22diBWK11063PIUtr/QbwRnx8/H1WxxJQ9XVGheaCZyFuPlz/Z58Xzjx8/BSb9lSwaZeRVDburuDQ8VMAhHVTnDOkH9dPGWYmlkjOHtzP63tfhBChwdt/2flAlD8CEQFUWw2v3Gckmkt+DDf8tcOJ5mhVDWt2HOKpD4v5/ouFXPL4KuJ+n8/dz63nryu3UVZeyeVnD+TXN07kle/P4IvfXMNbD17GolvO47sXjWLSsEhJNCEsOTm5YUHb/fuMjAxiYmKIiooiIyMDp9NJYmIiSiliYmIoKiry+hju8vPziYuLQyl1xvOJiYkkJhp3aSilUEqRm5vr8VxAo/MkJyfjdJ5u0dWWWJqTkZFBVFSUxxhbOqe71NRUYmJiGj3mdDpRSpGfn9/wWFpaGlFRUURFRZGWltYo/tTU1IY/r7bG7kvejmySgQKlVB5GAc5GU2la61d8FZjwk+pKWD4XduRBwq/h0p94fYiqmjq+/KaCjbuM7cYbdztxHDy95XhkdC8uGGVj3ozRnDfcxuTh/ekXERo3grbHb974ki3fHLU6jEYmDuvPr26c1O6fj4o68zNlSkoKmZltb1uVm5vL1KlTycvLIzMzk7S0NDIzM8nMzCQ9PZ377ruP5ORkiouLvTpGUlISdruxbOxwOEhPT8dut+NwOEhMTMRut5OUlEReXh65ubkkJyefMY3WlOuNODMzE7vdTnp6OmPHjqWkpKRh7aq1WDzJz88nLS2tYcOFw+FoSChtOad7fFlZWTgcjobzLV++HJvNRkKC0SQ5OTmZoqKiht1rrgTlmjLMz89n+fLlpKenM3t24DcWe5ts7gNizK+mNCAbBILZSSf85zbY+Snc+Ddj+qwVNXX1fL3vGBt3O82bJSvYtv8Ydeae47P69WTKCBs3XzCc88x1ltYqJ4vg52k3WnS0d91WbTYbCxYsAGDhwoVkZGSQkJDQ8OaYmpra6kK9p2O4v+G6r724kkx2djZJSUltjrOoqIisrCyKi4sbjpuZmUlBQUFDgmxLLJ44nU5sNlvDhgvXa9t6TpfY2FhsNltDona93pU0HA4Hubm5FBYWEhsbC0B6ejqLFi1q+DNyOByNzhdo3iabR4B04HGtdXNdO0UwOn4Alt1idNNMehYm33LGS7TWFB88wYZdTjabayxb9h6lutbYchzZqwdTRkQy69yYhgX8IZGy5bgjI4hg5e1uNE/i40+3OHEdKy7udOfWtiQvT8dwn2pyOp1kZWWRnZ3dMGpwJbO2KigoaJQIXBISEhqm3doSiycJCQlER0cTFRVFQkICiYmJLFiwoM3ndJeSkkJubi7p6ek4nU6KiopYvHgxQMNUmvufb1OxsbGWJRrwPtmUA5mSaEKMcycs/Q4c/Qa++zKMb/yPsa5e896X+8j6yMGGXcY/nt7h3Zk8PJJ500c3bDkeFd1bthyLDvF2dNSauLg47HY7ixcvJjY2tsV1DyvYbDaKi4vJzc0lOzubtLQ01q9f37CW5I3U1FQyMjIoKipqSFauUYzrXEeOHGn2561MNNC+QpxJwB/9EIvwh4NfG4mm5gTMfRVGXdzwVGV1LTkFu3nm4xJ2llcyekBvfnXjRC4ZN5CYQX3pLluORRArKirC4XCQl5fX7BtpW0Zn8fHxDSMF9zfv/Px8r0dJzUlKSiIpKalhDWnhwoVen9Nut2O328nOzqaoqKjRFGJCQoLH4wUTb5PNVOBWpdRC4IweylJ4M8jsKYJltxoVm+e/DUMmA0bJ/aVrS3nh0zKclTVcMNLGwuvO5epJQyTBCMCYWvI0CvHF9JqvuN5U09LSSE9Pp6ioiNzc3EZv1q4klJub2xB30zfz2NhYkpKSSE5ObthssGjRIoqKili5cmWHYszNzSUvL4/k5GSio6PJzs7Gbre3+5ypqakNycZ9Xcdut5OSkkJycnLDhgPXa4KlSkJ7ytVI0c1QULLa2AzQO9oopjkghh0HjvPMxw5WFO2hpq6ehAmDSb3cTtzoKJkeE400N82Tk5Pj1eK7v6Wnp5OWlkZ+fj6zZ88+I5G43tiTk5OJjY1tWFxvKicnh7S0NNLS0igvLyc+Pp7i4uIOJ9bY2Fiys7MbpvdiY2Mbdou155wpKSmkpaU1jHLcZWZmkpGRQWpqKg6Hg9jYWBYuXNjMkQJPtbYlsKuIj4/XBQVnDNZC01dvQ858iB6LvvMVPjscweLVDvK3HiA8rBtJcSO459KxPim535Vs3bqVCRMmWB2GEAHTlt95pVSh1jq+xRfRjpGNUqo/xo4010eIPCBNa33M22MJP9j4Mrz6ffTQ88m78F/8a1kpG3c5ierdgwdnjWfu9NEM7BsclZyFEF2Ht4U4I4ESjPbQrum0BMCplLJrrct8HJ/wxrpMeGcB30RfxPwjP2bbihLGDOjN774zmaTYEfQKl9ughBDW8HZkkwUs11rf7/6gUioHY6fatb4KTHhBa46//wf6rn2CVUzlgW++z8RRNp663k7iRFn0F0JYrz1toT3tq3sMo26aCLAd+yvYt/ynXHo4l5y6y1k1/he8NHM8caN9ez+DEEJ0RHtu6hwLlDZ53E4naTkQCrTWrCsp55kPt3Ot4/fc2n01nwycTdzsP5F8Vn+rwxNCiDO0ZxotSymV4FqfUUpdADwOeK6xEORCqZ9NbV097365j8UfOfhq90EyI/7JFd0LqLzkYS5JeBhk+7IQIkh5lWy01hlKqRigRCnlqotgA3K11sGzodsLodDP5sSpWpYX7OKZj0vYfeQkE6IVHw77F0PKC+C6J+g9TRpBCSGCm9dbn7XWqUqpTMBV8a1Aa/25b8MSAAeOVrFkbSnLPt1Jxcka4kZH8dvEIVxZ8H3U3k1wy2KYEvhS4UII4a32VBBAa12E0c9G+MH2/cdYvNrBq59/Q019PddMHMJ9l48lznYSXrgZnGVw24twznVWhyqEEG3SbLJRSt0CTG3L9JhSaj1wr9Z6oy+D60q01nzqKGfxagervjpAz7BuzJ46gnsvtTNmYB84XAzP3mT0pLlzBYy51OqQhRCizVoa2TwCrG/jcVYBGYAU4vRSbV0973yxj8WrHWzaXcGAPuH8JOFs7rx4FANcd/rv3WT0otH1MP9NGHaBtUELIYSXWmr6bgdy2nic94FWa+OI006cquXZj0uY+cQH/N9/PudYVS1/uHkynzx8FT9KGH860ZSthedvgO494e53JdGIgHD1qldKERUVRWJiYqNe94GSnJxMcnJywM9r9bk7o5aSjY223ztTbr5etOLA0Soy3v2K6YtW8ts3tzDMFkHWXXGs/OlM7pg2mogebiVltucZazR9B8H33oVBZ1sXuOgycnNzSUxMJCEhgcLCQnJycrDZbI0qJsfFxbWrAZjoulqaRnMAs4ANbThOArJhoEXb9h9j8UcOXttgLPpfO2kI911uJ3ZUlOcf+GIFvJICZ02EO18xEo4QAZCWlsaCBQsaJRdXcy6X1NRUK0ITIaylZJMFPKKUym2pwKZSaiywEKNkjXCjtWat4zCLP3Lwv68PEtGjG7ddNJLvXTLWWPRvTsGz8OZPYdR0uP1liIgMXNCiy3M4HMTExJzxuHufFfcukUK0RbPJxryBcw5QpJRaoLV+pulrzB1ri4H1WmtpFW2qravnrc17WbzawRd7jjKgTzg/TTybOy8eTXSf8OZ/UGv4+M+w8rcw/hqYvQR69Apc4EJgNBzLzMyUhCJ8qqU1G4CrMHaaLVZKHVZKrVdKZSul3lNKHcbYQLASkDsLgeOnannGXPT/0csbqDxVx2M3n8cnD1/Fg7PGt55o8h41Es15ycZ9NJJohAUyMzMpKioiKiqK1NRUjxsDmi6eu77PyMggJiaGqKgoMjIycDqdJCYmopQiJiaGoqKiRj+TmJjY8DNKqTZNz6WlpREVFUVUVBRpaWmtvj4jI4OoqKiGzQ4ZGRmNnnPF6+qm6UlqauoZoz2n04lSqtGfT0uxJScnN/x5xsXFNYqjK2jxpk6tdQWQrJSyA2kYRTjjMNZzFgPZUj0A9h+t4vk1pbz4aRlHq2q5aEw0v/72JGadexbd2lLev74O3vgRfP4CTL0XrnsCurX2OUAEjXcehn2brY6isSHnwXWPt+tHExISKC4uJi0tjeXLl5OVlYXdbicvL++MVsTucnNzmTp1Knl5eWRmZpKWlkZmZiaZmZmkp6dz3333kZycTHFxccPP5Ofnk5iYSE5ODtnZ2Q1v/gsWLPB4juTkZIqKihpaK7uSQHOjsPz8fNLS0igsLMRms+FwOBoSiuuNPzMzE7vdTnp6OmPHjqWkpOSM1sypqalkZWXhcDga/gyWL1+OzWZraEXdltjy8/NZvnw56enpzJ7dtT6jt6mCgNbaAciKoAf7Kqq4PON/1NbXc+3kIdx3mZ0Lm1v096T2FLxyH2x5DS5/CK78uRTUFJaz2+0Nb5r5+fkNoxD3RNGUzWZrSBILFy4kIyODhISEhjfj1NTUM0YuCQkJDT8TGxuLw+EgMzPTY7JxOBzk5uZSWFhIbKzR6SQ9PZ1FixY1m2ycTic2mw273d7wX4CioiKysrIoLi5ueCwzM5OCgoKGJOkuNjYWm83WkDhdr3cljLbG5nA4Gp2zK2lXuRpx2pDICB751rlcee5ZjB7QwqK/J9UnIPtOKF4FV/8BZvzQP0EK/2rnCCJUJCQksHLlSuLi4sjPz29IHk3Fx5++1c41MoiLi2t4LDq69R5LiYmJ5OZ6LiDvmq5yP2ZbYo+OjiYqKoqEhAQSExNZsGABBQUFjZKP++ubO39KSgq5ubmkp6fjdDopKipi8eLFXsUWGxvbJRMNtL5mI9pg/iVjvU80leWw9Dvg+AC+/U9JNCKouT6tN7em0Zy2JBhv2Gw2tNZnfLX0+uLi4oZ7hdLS0tp9o2ZqaioOh4OioqKGKTTXn0tbY+uqiQYk2Vjj2D54/nrYuwGSl0DsXVZHJEQDT5/sXZ/c3d9c/SEnJ6fZc7ju9XHfZNBWSUlJ5OTkkJOTQ25uLvHx8R6P1dLIzW63Y7fbyc7OJicnp9H0WEdi6yok2QRaeQk8ew0cKYM7cmDit62OSIgGDoejYX0mKyuL3NxcMjIySExMJCUlxeefzPPz8xsW3l0L9gsXeq79a7fbSUlJITk5mfz8/IZ1kpYqGeTm5jYct6ioiOzsbOx2O7GxsSQlJZGcnExubi5FRUUNC/zuN7M25TpWfn5+o/Wn9sTW1UiyCaT9W+DZa43KzfNeB/sVVkckRCN2u71hATs9PZ3k5OSGRfGmi+a+EBsbS2FhYcN6UGZmJklJSc2+PjMzs2GjQUxMDIsWLWpxu3RsbCzl5eUkJycTFxeHw+Fo2PiQk5NDUlISaWlpzJo1C6fTSXFx8Rk70dylpKRQVFTUMMrpSGxdjWppvrMriY+P1wUFBf47we4CWHYrhEXAXf+FwRP9dy7hF1u3bmXChAlWh9FpuO5rycvLszoU0Yy2/M4rpQq11q0WYu6Uu9GUUklANBADFGutsywNqPh/8PIdRn2zua9B1BhLwxFCiEDrdMnGvAHVqbXONb8vVEoVmN1FA2/L67DiHhgwHu56BfoNsSQMIYSwUmdcs4ml8Q2oBRhVqQOv6AXImQdDL4C735JEI4TosjrdyEZrnauUci/mFE/bm8D5zpp/wvs/h5irYM4yCPfyPhwhOjnXQr3oGixJNkqpWGCh1vqMu6uUUimcbtpm11p7Xa1Oa+00j2UHyrXWgWsxqDWs+j2s/iNM/A7ckgVhPQN2eiGECEYBTTZmkpljfnvGhn0z0eC23mJXSmVqrVPN7xcAA5o5fJ6HpJKmtQ7cRvf6enjnIVj/NMTOhRv+Ct26t/5zQgjRyQU02ZiL9EVm0vG0jpKqtY5ze71DKZXg9n2bRzlmYkoz/9/mGu34TV0NvPoAbM6BGQ9C4m+loGYnpLVGyd+r6AJ8fVtM0GwQUErZMBb3m3K6J5w2HisByHJLMP7dIFBdaWxt3pwDs34FV/9OEk0nFBYWRm1trdVhCBEQtbW1hIX5bjwSNMkGY1rN0+ijHM9JyCNz1JQHHFFKaaWUxui/4x9VFcbNmtvfhxv+Apf91G+nEtaKiIjg+PHjVochREAcO3aMiIgInx0vmHajRXN6Y4A7J82v05zBnKpr07DCXCNKARg1alRbT9HkhPVQexKSnoHJt7bvGCIkDBo0iJ07d9KzZ0969eol02miU9Jac/LkSQ4dOtT+90UPginZBJxZWSALjHI17TpIryi4d6VsBOgCIiIiGDx4MPv27ePUqVNWhyOE3/Ts2ZPBgwd32pENGKObpmzA4UAH4hVJNF1GZGQkkZGRVochRMgJpjWbAozE0lQ0IE0ihBAihAVNsjF3jjnMXWnubAG9KVMIIYTPWZVsmusVmw40dE4yd5b5NdEopW5USmVVVFT48zRCCNGlBTTZmBUB0jGSSqxSKtNVNQAaFuyLlVIJZpuABFf1AH/RWr+htU6ReXghhPCfQFcQcGDe1d/Ca6ztPSOEEMLngmbNRgghROfV5ZONrNkIIYT/KV8XWwtVSqmDQJmHpwYChwIcjj90lusAuZZg1FmuA+RavDVaaz2otRdJsmmF2VI63uo4OqqzXAfItQSjznIdINfiL11+Gk0IIYT/SbIRQgjhd5JsWtdZtmJ3lusAuZZg1FmuA+Ra/ELWbIQQQvidjGyEEEL4nSQbIUSnZ5bK8m97eNGiYOtnExTMem2urqF2rXWGlfG0xCxWulBrnezhuRavI9iu061OXpz53zSzGrj780F/PWYcrurlMUC6WarJ/fmgv46mlFI5TX/PQuhaYoHFZlV5J0ZLkzSzsy8QUtfiimcBxrWUA2itc92eC75r0VrLl9sXRpvoFLfv7UCm1XF5iDMWo6BpOlDo7XUE23W6x+IWX3GoXQ+woMn3SaF4Hc38vulQ/R0Dksz/2pr7/QuVazHPn4ORJFzfa9e1Beu1WPoLHIxfzbxxFwc6Di/ijW0m5havI5iuE2MUkOLh8SMYlb9D5nqAYtcbm/m9vckbQUhch4cYEjwkm5C5Fve/k2aeD6VrSeHMDzXuiScor0XWbNyYQ+xYD085Q2m+t7XrCMLrtAOZHhrnOQB7iF1PonabzsC4NqfW2hli19FAKZWkmzQwDNVr8SQEryUdcP8dQ5vTtMF8LZJsGrNjzIE2VY7nv6Bg1dp1BNV1amPePE67rc+Y7JgJhxC5Hu22NmNKA1zrHCFzHS7mmqCntuyheC0Jbl8L3D7chMy1mDHbzP9PCqVrkQ0CjUVzesesJ/oAAAWwSURBVNHMnRMYEOBYOqK16wi669RuC7Vg/EMCHFrrfPMTV0hdjxl/IsbmANeoIOT+XjCmZ3I9PB5q11IEjUYADox1j0RC61pcycLm+ntRShUAKzE21gTttcjIRgQd81PaQmCW1bG0l9Y6VxtdZmOVUplWx9Me5vSZp0QTcrTWDvdRp/n/dnPkFkqiMUY27tfiBGPkZlVQbSHJ5kzRHh6zAYcDHUgHtXYdwXyd6UByk2m1kLwebWwpne32RhAS16GUck1htiQkrqUFTsBVETlUrsUBpxOMG/dpsKC8FplGa6yA0/dHuIvG87x1sGrtOoL2Os17B9KbrH2ExPWYn5JXaq2jmjzlwJiuWUQIXIcpAbA1/bTsdm/HckLkWszEWay1Vk2eKje/QuL3C4wRmVJNL6OB6/6h4LwWK7buBfMXxtZVW9PHrI6rhXib2/rc4nUE43VibOm0N3nMtfU56K8H4w063cPjRzC3dofCdbRwfToUf8doeWu9PZSuxTxvoYd/J8VAbDBfi0yjnSkdY70AaPi0mt/8yy3naUgMrV9HUF2n+Qm6QLtt4WzyqTror0c32R7sFgcYIwEIgevwQkhciz5zysl1B/1yfXoEHRLXYkozv9xjcejTm2yC8lqk6rMH5i+iA+MTkeVlKTwxpwZSMT5Nx2KUEi/UWme5vabF6wiW63RNczTzdJQ+vQAa9Ndjbm5IcXuouXI1QX0dTeJJwPhdS8K4vyPTlVhD6VrcpgBt0LCe5v58KF1LEsbONIABWuu0Js8H3bVIshFCCOF3Mo0mhBDC7yTZCCGE8DtJNkIIIfxOko0QQgi/k2QjhBDC7yTZiC5DKZWilNIevvLM7ddCCD+RZCO6oiittdJG+ZI4jHsvCj2UZilUSuVZEmErgjk2ITyR2miiSzPvuk42KzPnKKXGut1xHszVmoM5NiHOIDd1ii7DvGs6E7eqBE2ePwJkNb0bWwjRcTKNJsRpyzFKsgghfEySjRCnFXK63hRKqRylVE7T7802vMVKqSOulrzmJgNtPt6oIZdSKt187RGlVHprx2vyswvMx7X7801jc3ut6zg5bq2C23w+b5gbLhaY/7Wb/00xzyEbLkQjkmyEOM1VcbqlN0rXyCcRo/hpOkaSSuf0ZoNGCcr8mWTzK8mczmv2eK7zmxsW0jE6lsaYP++xoZm55pRqfsVh9GkpaZpwWjqfN9yaq+WbMSZorbPMQrDrcatKLARIshHCnetN11OPdhen1jrDrOK8yHwsX2udb242yHQdx3xDTsLoOppvVkpOw0gILR3PFYcNI3k5tNHWOF97aNNsjqRSgETzNQ5ttKR2YCSC1uJvzyjEbl5PvBlflttzMWbcQjSQ3WhCnBaD8Wbc0htlget/tNZOs2tiodvz7onKtZW6sIXuip6O5xqN5JvHO6KUygfymikFH2/G3XTUk8+Za1Atna/N3Hr3JALZTZ52tSQQooGMbIQ4LYXTTc680dpISDX9astBtdZOrbVr+syJMeWV08qPBVoCbo23zNGcXXtoJCe6Nkk2QmAs4mN8wvflWkM+YGu6YcBbWutcrXXDmo+HlxQ0c55GicDXzMRic+sQCacbrLkafAkBSLIRXZxSKtYcLSzAWPPw2VqDOa2VhXGzaIK5YyuprXf+m6/NNH82FpiDhw0C5pt9rnmeJLdrisXL5GnG2NYdap6S2VQgz9yY0FzLctEFSbIRXZFrK7HGWG+xATH+mPoxF+ozza9ijN7vbb37vwjjDTuH09uyk5s5TzJGwkkHVnL6mrxNnq4dcG0Rw5nXkmk+ntJk04Do4qSCgBCiEaXUEa11lNVxiM5FRjZCiAbmdF17NkkI0SIZ2QghGjSz6C9Eh0myEUII4XcyjSaEEMLvJNkIIYTwO0k2Qggh/E6SjRBCCL+TZCOEEMLvJNkIIYTwO0k2Qggh/O7/A1jI7KX1G0ttAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1519130c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(n_list, times_elim, label=\"Elimination solver\")\n",
    "plt.semilogy(n_list, times_simple, label=\"Simple solver\")\n",
    "# plt.semilogy(n_list, 10**(-6)*np.array(n_list)**3 / 3)\n",
    "plt.legend(fontsize=18)\n",
    "plt.xlabel(\"Dimension, $n$\", fontsize=18)\n",
    "plt.ylabel(\"Computation time, sec.\", fontsize=18)\n",
    "plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Резюме \n",
    "\n",
    "1. История развития теории линейного программирования и текущее состояние\n",
    "2. Концепция методов внутренней точки\n",
    "3. Прямой метод решения задачи линейного программирования"
   ]
  }
 ],
 "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
}