preliminaries/dcp_autodiff/dcp_autodiff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# DCP & Autograd \n",
    "# или \n",
    "# как правильно ставить задачи и автоматически вычислять градиенты"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Задача оптимизации\n",
    "\n",
    "\\begin{align*}\n",
    "& \\min_{x \\in \\mathbb{R}^n} \\; f_0(x)\\\\\n",
    "\\text{s.t. } & f_i(x) \\leq 0, \\; i = 1,\\ldots, m \\\\\n",
    "& h_j(x) = 0, \\; j = 1,\\ldots,p \n",
    "\\end{align*}\n",
    "\n",
    "- Возможность эффективного решения сильно зависит от свойств $f_0, f_i, h_j$\n",
    "- Если $f_0, f_i, h_j$ аффинны, то это задача линейного программирования (LP), которая может быть решена крайне быстро\n",
    "- Простые задачи с нелинейными $f_i, h_j$ могут быть очень сложными для решения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Задача выпуклой оптимизации\n",
    "\n",
    "\\begin{align*}\n",
    "& \\min_{x \\in \\mathbb{R}^n} f_0(x)\\\\\n",
    "\\text{s.t. } & f_i(x) \\leq 0, \\; i = 1,\\ldots,m \\\\\n",
    "& Ax = b \\\\\n",
    "\\end{align*}\n",
    "\n",
    "- $f_0, f_i$ выпуклые функции: для всех $x, y$ и $\\alpha \\in [0, 1]$\n",
    "\n",
    "$$\n",
    "f(\\alpha x + (1 - \\alpha)y) \\leq \\alpha f(x) + (1 - \\alpha)f(y)\n",
    "$$\n",
    "\n",
    "- Ограничения типа равенств аффинны"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Свойства задач выпуклой оптимизации\n",
    "\n",
    "- Подмножество задач оптимизации: LP - частный случай\n",
    "- Могут выглядеть очень сложно, однако решаются также эффективно как и задача LP\n",
    "- Встречаются гораздо чаще, чем можно было бы подумать\n",
    "- Очень много приложений, часть обсудили на прошлом семинаре"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Общие подходы к использованию выпуклости\n",
    "\n",
    "- Надеяться/предполагать/делать вид, что $f_i$ выпуклы\n",
    "    - Просто для пользователя\n",
    "    - Теряется часть преимуществ выпуклых задач\n",
    "- Проверка выпуклости задачи перед решением\n",
    "    - в общем случае может быть затруднительна  \n",
    "- Построение выпуклой задачи из элементарных блоков\n",
    "    - пользователь следует фиксированному набору правил при определении $f_i$\n",
    "    - выпуклость проверяется автоматически"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Как проверить выпуклость?\n",
    "\n",
    "- Определение, критерии первого или второго порядка, например $\\nabla^2 f(x) \\succeq 0$\n",
    "- Исчисление выпуклых функций: построение $f$ специальным образом\n",
    "    - Дан набор простых функций, выпуклость которых известна\n",
    "    - Даны сочетания и преобразования, не меняющие выпуклость"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Примеры простых выпуклых функций\n",
    "\n",
    "- При $x > 0$: $x^p$ для $p < 0, \\; p \\geq 1$ и $x ^{-p}$ для $p \\in [0, 1]$\n",
    "- $e^x, \\; -\\log x, \\; x\\log x$\n",
    "- $\\langle a, x \\rangle + b$\n",
    "- $\\|x\\|$ - любая норма\n",
    "- $\\max \\{x_1, \\ldots, x_n \\}$ и $\\log(e^{x_1} + \\ldots + e^{x_n})$\n",
    "- $\\log\\det X^{-1}$ для $X \\in \\mathbb{S}^n_{+}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Правила исчисления выпуклых функций\n",
    "\n",
    "- Умножение на неотрицательную константу: $f$ выпукла и $\\alpha > 0$, тогда $\\alpha f$ выпукла\n",
    "- Сложение: $f, g$ выпуклы, тогда $f + g$ выпукла\n",
    "- Композиция с аффинной функцией: $f$ выпукла, тогда $f(Ax + b)$ также выпукла\n",
    "- Взятие максимума: $f_1,\\ldots, f_m$ выпуклы, тогда $\\max\\limits_{i=1,\\ldots,m} \\{f_i(x)\\}$ выпукла\n",
    "- Композиция: если $h$ выпукла и возрастает, $f$ выпукла, тогда $g(x) = h(f(x))$ выпукла\n",
    "- И многие другие..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Примеры\n",
    "\n",
    "- $f(x) = \\max\\limits_{i=1,\\ldots,m} (\\langle a_i, x \\rangle + b_i)$\n",
    "- $\\ell_1$ регуляризация задачи наименьших квадратов\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{1}{2} \\|Ax - b\\|_2^2 + \\lambda \\|x\\|_1, \\quad \\lambda > 0\n",
    "$$\n",
    "\n",
    "- Максимальное собственное значение $A \\in \\mathbb{S}^n$: \n",
    "\n",
    "$$\n",
    "\\lambda_{\\max}(A) = \\sup\\limits_{\\|x\\|_2 = 1} (x^{\\top}Ax )\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Как решать задачу выпуклой оптимизации?\n",
    "\n",
    "- Использовать \"стандартный\" солвер (для LP, QP, SDP...), примеры таких солверов были на прошлом семинаре (например [sdpt](https://github.com/SQLP/SDPT3) или [MOSEK](https://www.mosek.com/))\n",
    "    - лёгкий путь\n",
    "    - задача **должна быть** в стандартной форме для выбранного солвера\n",
    "    - сложность разработки компенсируется количеством пользователей\n",
    "- Придумать и/или реализовать метод самостоятельно\n",
    "    - Трудоёмко\n",
    "    - Может быть эффективнее для конкретной задаче\n",
    "- Преобразовать задачу к стандартному виду и использовать стандартный солвер\n",
    "    - Расширяет множество задач, подходящих для решения\n",
    "    - Преобразование может быть громоздким"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Методы внутренней точки (IPM) для выпуклых задач\n",
    "- [Interior-Point Polynomial Algorithms in Convex Programming, Y. Nesterov, A. Nemirovskii, 1994](https://epubs.siam.org/doi/abs/10.1137/1.9781611970791)\n",
    "- Обзор про IPM см. [тут](https://people.kth.se/~andersf/doc/sirev41494.pdf)\n",
    "- Применим для **гладких** $f_i$ и задач в конической форме (SOCP, SDP)\n",
    "\n",
    "\\begin{align*}\n",
    "& \\min f_0(x)\\\\\n",
    "\\text{s.t. } & Ax = b\\\\\n",
    "& x \\in \\mathcal{K},\n",
    "\\end{align*}\n",
    "\n",
    "где $\\mathcal{K}$ - некоторый конус (обычно $\\mathbb{R}^n_{+}$, $\\mathbb{S}^n_{+}$ или конус второго порядка) или декартово произведение конусов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### А если IPM нельзя применить к задаче?\n",
    "\n",
    "\n",
    "- Пример: $\\ell_1$ регуляризация задачи наименьших квадратов\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{1}{2} \\|Ax - b\\|_2^2 + \\lambda \\|x\\|_1, \\quad \\lambda > 0\n",
    "$$\n",
    "\n",
    "- Задача выпукла, но $f$ негладкая!\n",
    "- **Основная идея**: изменить задачу так, чтобы IPM можно было применять\n",
    "- Даже если в новой задаче будет больше переменных и ограничений, она может быть эффективна решена с помощью IPM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Пример\n",
    "\n",
    "- Исходная задача: $n$ переменных, нет ограничений \n",
    "\n",
    "$$\n",
    "\\min_x \\frac{1}{2} \\|Ax - b\\|_2^2 + \\lambda \\|x\\|_1, \\quad \\lambda > 0\n",
    "$$\n",
    "\n",
    "- Введём новую переменную $t \\in \\mathbb{R}^n$ и новые ограничения $|x_i| \\leq t_i$:\n",
    "\n",
    "\\begin{align*}\n",
    "& \\min_{(x, t)} \\frac{1}{2} \\|Ax - b\\|_2^2 + \\lambda \\mathbf{1}^{\\top}t\\\\\n",
    "\\text{s.t. } & -t \\leq x \\leq t\n",
    "\\end{align*}\n",
    "\n",
    "- В новой задаче $2n$ переменных и $2n$ ограничений, но она **гладкая**! \n",
    "- **Важно**: задачи эквивалентны! Решив одну, получим решение другой и наоборот"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Преобразование задачи и эффективность решения\n",
    "\n",
    "- Дана выпуклая задача $P_0$\n",
    "- Выполняются последовательные эквивалентные преобразования\n",
    "\n",
    "$$ \n",
    "P_0 \\to P_1 \\to \\ldots \\to P_K,\n",
    "$$\n",
    "\n",
    "где $P_K$ - задача, которую можно решать IPM\n",
    "- Эффективное решение $P_K$\n",
    "- Обратное преобразование решения $P_K$ в решение $P_0$\n",
    "- $P_K$ может иметь больше ограничений и/или переменных, но наличие определённой структуры и высокая эффективность IPM компенсируют это"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Примеры преобразований задач\n",
    "\n",
    "- Правила преобразования выпуклых функций порождают преобразования задач\n",
    "- $\\max \\{ f_1(x), f_2(x) \\}$   \n",
    "    - Вводим новую переменную $t = \\max \\{ f_1(x), f_2(x) \\}$\n",
    "    - Добавляем ограничения $f_1(x) \\leq t, \\; f_2(x) \\leq t$\n",
    "- $h(f(x))$    \n",
    "    - Вводим новую переменную $t = f(x)$\n",
    "    - Добавляем ограничение $f(x) \\leq t$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### От доказательства выпуклости к применимости IPM\n",
    "\\begin{align*}\n",
    "& \\min_{x \\in \\mathbb{R}^n} f_0(x)\\\\\n",
    "\\text{s.t. } & f_i(x) \\leq 0, \\; i = 1,\\ldots,m \\\\\n",
    "& Ax = b \\\\\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "- Построение $f_i$ из элементарных функций и правил преобразований даёт доказательство выпуклости\n",
    "- Аналогичный разбор даёт преобразование задачи к форме, состоящей из элементарных функций и аффинных равенств \n",
    "- Если элементарные функции подходят для IPM, преобразование автоматически даёт форму задачи, которая может быть решена IPM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Disciplined convex programming (DCP)\n",
    "\n",
    "- Задаются искомые переменные и фиксированные параметры\n",
    "- Целевая функция и ограничения строятся из элементарных функций с помощью правил композиций и сочетаний\n",
    "- Задача выпукла по построению\n",
    "- Автоматически разбирается на элементы\n",
    "- Приводится к форме для запуска IPM\n",
    "- Решается некоторым стандартным пакетом для IPM\n",
    "- Восстанавливается решение исходной задачи\n",
    "\n",
    "\n",
    "Визуализацию и результат разбора можно посмотреть на [сайте](http://dcp.stanford.edu/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"./dcp_parse.png\" width=800>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### История\n",
    "\n",
    "- Системы [AMPL](https://ampl.com/), [GAMS](https://en.wikipedia.org/wiki/General_Algebraic_Modeling_System) - 1970-ые\n",
    "- Пакеты для задач SDP/LMI: <tt>sdpsol</tt> (Wu, Boyd),\n",
    "<tt>lmilab</tt> (Gahinet, Nemirovsky), <tt>lmitool</tt> (El Ghaoui) - 1990-ые \n",
    "- <tt>yalmip </tt> (L&ouml;fberg 2000–)\n",
    "- automated convexity checking (Crusius PhD thesis 2002)\n",
    "- disciplined convex programming (DCP) (Grant, Boyd, Ye 2004)\n",
    "- <tt>cvx</tt> (Grant, Boyd, Ye 2005) для MATLAB\n",
    "- <tt>cvxopt</tt> (Dahl, Vandenberghe 2005)\n",
    "- <tt>cvxpy</tt> (Diamond, Boyd 2016) для Python\n",
    "- disciplined geometric programming (DGP) (Agarwal et al 2018) - [paper](https://web.stanford.edu/~boyd/papers/pdf/dgp.pdf)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Главное по DCP\n",
    "Pro:\n",
    "\n",
    "- Проверка выпуклости и генерация преобразования задачи для IPM\n",
    "- Построение задачи: элементарные выпуклые функции + правила композиций и преобразований \n",
    "- Очень похоже на математическую нотацию\n",
    "\n",
    "Contra:\n",
    "\n",
    "- Не про \"plug & play\" или \"try my code\"\n",
    "\n",
    "<img src=dcp_mem.jpg >\n",
    "\n",
    "- Нельзя записать произвольную задачу и надеяться, что она будет выпукла"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Солверы для решения задач оптимизации oбщего вида\n",
    "\n",
    "- [ipopt](https://projects.coin-or.org/Ipopt)\n",
    "- [Pyomo](http://www.pyomo.org/)\n",
    "- [Gurobi](http://www.gurobi.com/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['ECOS', 'ECOS_BB', 'CVXOPT', 'GLPK', 'GLPK_MI', 'SCS', 'MOSEK', 'OSQP']\n"
     ]
    }
   ],
   "source": [
    "import cvxpy as cvx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "print(cvx.installed_solvers())\n",
    "\n",
    "USE_COLAB = False\n",
    "if USE_COLAB == False:\n",
    "    plt.rc(\"text\", usetex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## DCP на примере CVXPy\n",
    "\n",
    "Модельная задача - [compressed sensing](https://www.raeng.org.uk/publications/other/candes-presentation-frontiers-of-engineering)\n",
    "\n",
    "\\begin{align*}\n",
    "& \\min_x \\|x\\|_1\\\\\n",
    "\\text{s.t. } & Ax = b,\n",
    "\\end{align*}\n",
    "\n",
    "где $A \\in \\mathbb{R}^{m \\times n}$ и $n \\gg m$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Генерация данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Num of nnz in x = 42\n"
     ]
    }
   ],
   "source": [
    "n = 1000\n",
    "m = 10\n",
    "x_true = np.random.randn(n)\n",
    "x_true[np.abs(x_true) > 0.05] = 0\n",
    "print(\"Num of nnz in x = {}\".format(np.sum(x_true != 0)))\n",
    "A = np.random.randn(m, n)\n",
    "b = A.dot(x_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### CVXPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------------------------\n",
      "           OSQP v0.5.0  -  Operator Splitting QP Solver\n",
      "              (c) Bartolomeo Stellato,  Goran Banjac\n",
      "        University of Oxford  -  Stanford University 2018\n",
      "-----------------------------------------------------------------\n",
      "problem:  variables n = 2000, constraints m = 2010\n",
      "          nnz(P) + nnz(A) = 14000\n",
      "settings: linear system solver = qdldl,\n",
      "          eps_abs = 1.0e-03, eps_rel = 1.0e-03,\n",
      "          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
      "          rho = 1.00e-01 (adaptive),\n",
      "          sigma = 1.00e-06, alpha = 1.60, max_iter = 3000\n",
      "          check_termination: on (interval 25),\n",
      "          scaling: on, scaled_termination: off\n",
      "          warm start: on, polish: on\n",
      "\n",
      "iter   objective    pri res    dua res    rho        time\n",
      "   1  -8.0000e+03   8.00e+00   6.68e+01   1.00e-01   4.68e-03s\n",
      " 200   2.9266e-01   2.64e-03   2.38e-02   3.27e+00   7.63e-03s\n",
      " 400   2.6957e-01   3.28e-04   1.08e-02   3.27e+00   2.68e-02s\n",
      " 600   2.7603e-01   1.36e-04   8.75e-03   3.27e+00   4.36e-02s\n",
      " 750   2.7066e-01   3.04e-04   1.77e-03   3.27e+00   5.68e-02s\n",
      "\n",
      "status:               solved\n",
      "solution polish:      unsuccessful\n",
      "number of iterations: 750\n",
      "optimal objective:    0.2707\n",
      "run time:             5.91e-02s\n",
      "optimal rho estimate: 2.75e+00\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.2706647563479417"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = cvx.Variable(n)\n",
    "objective = cvx.norm1(x)\n",
    "constr = [A*x == b]\n",
    "problem = cvx.Problem(cvx.Minimize(objective), constr) \n",
    "problem.solve(verbose=True, max_iter=3000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Проверим решение"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.366577785575698e-07\n",
      "0.02806119193399998\n",
      "0.2170992263162404\n",
      "Num nnz = 200\n"
     ]
    }
   ],
   "source": [
    "tol= 1e-4\n",
    "print(np.linalg.norm(A.dot(x.value) - b))\n",
    "print(np.linalg.norm(A[:, np.abs(x.value) > tol].dot(x.value[np.abs(x.value) > tol]) - b))\n",
    "print(np.linalg.norm(x_true - x.value))\n",
    "print(\"Num nnz = {}\".format(np.sum(np.abs(x.value) > tol)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-0.1 , -0.08, -0.06, -0.04, -0.02,  0.  ,  0.02,  0.04,  0.06]),\n",
       " <a list of 9 Text yticklabel objects>)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoUAAAHVCAYAAACHT4G+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XmQXdd9H/jv7y29Lw8NNADuYLdESaZkiw1Q1mJbdrEhyU4caxKAdMbO2PGE6NgzrkQZBzAnSQ0TuwpuTE25JlUjp0EvWcqJyG5nEluKY6M5tuUtEYH2JlGUKbQMieICgI2H3rvfcuaPe+/r97rfct9dz7n3+6mCmnrLvedu5/7e7yxXlFIgIiIionTLxF0AIiIiIoofg0IiIiIiYlBIRERERAwKiYiIiAgMComIiIgIDAqJiIiICAwKiYiIiAgMComIiIgIDAqJiIiICEAu7gKY6MiRI+rEiRNxF4OIiIioo2vXrt1WSo13+hyDQg9OnDiBq1evxl0MIiIioo5E5Iabz7H5mIiIiIgYFBIRERERg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiSrDb6zs48dOfw39cei3uohARaY9BIREl1vKtDQDAf/jC12MuCRGR/hgUEhERERGDQiIiIiJiUEhEREREYFBIRERERGBQSERERERgUEhEREREYFBIRERERGBQSERERERgUEhERERBef1Pgdf/JO5SkEe5uAtARERECXH5o9bfZ+/GWw7yhJlCIiIiImJQSEREREQMConCVa0Av/EPgduvxl0SIiKithgUEoXp5svAtV8B5n807pIQERG1xaCQiIiIiBgUEhERERGDQiIiIiICg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCovT64n8E7n4z7lIQEZEmGBQSpVGlDCz8XeBXvjfukhARkSYYFBKlkrL+rDJTSEREllxcKxaRKQAzAK7bLx0GcFEpVYx6WSIyAeCCUmqm23UTERERJUEsQaGITAOYVUqdrHutAOBFEXmim8AwoGXNA1h2vwVEhlMq7hIQEZFm4mo+ngfwdP0LdvA2Z78X2bJE5ByAQpfrJDIcg0IiImoUeVAoImcAQCm11OTtRQDTdnNu6MuyM4oAs4SUNswUEhHRPnFkCp9CiyBMKeW8fiaiZZ1TSl12uS6iBGFQSEQpt/o68I0vxF0KrcQRFE4BWGnzfhHA42Evy84yLrpcD1GyMFNIRGn3L6eAXzoddym0EkdQOAYrWGtlBe77+Hlalt1sPNGi2bkpETknIldF5OqtW7fcfo1IUwwKiSjlyltxl0A7cQSFbgK+sZCXdU4pdcnlOgAASqnLSqlTSqlT4+Pj3XyVSD9OppAZQyIisuk6eXWQo4EblmVPYeM6Q0iUTAwGiYioka5BYdcTWHexrCmlFPsSUroxQ0hERPvoGhS2GzzieVkich4ARxsTpSRTqBj8EhG5FkdQuIz2zcOdBo94WpY9X+Gyl8foESVOSoKldGwlEVEw4njM3SKA6TbvFwBcCWFZUwCeEpGnmnxuGkBRRJwnoFzsZmQykXnscEkk3mIQEZE24ggKrwA41+yNuqePvBD0spRSCwAWWnz2OoAlpdRZl+slMltaMoXp2EwiokBE3nxsB2dFexTwfjMAFvc38YpIwZ4nsKGp2MuyiCKla1Sia7kCptiATETkWlwDTc4CmK1/wc7snYEVzO33HIA5+6/fZbUS5DQ4RDZdgxJdyxWwlGwmEVEQ4mg+hlJqUUSeFpE5ANftlycBnGyR2XseVr+/5wNYFoDaU02eg9XXcALAhIhcgzUYhc3IFAxdM3K6litg6dhKIqJgxBIUAoA9kMNVJq9df8Bul1X3nSKsLCNRiHQNS3QtV7BSEvsSEQVC13kKiZJB16hE13IFjH0KiYjcY1BIFCpdgxJdy0VERHFhUEgUJl1jr7RkCtOxmUREgWBQSBQqXaMSXcsVrHRsJVEEvvAccOOP4i4FhSy2gSZEqaBrqkrXcgWMzz4mCsh/+Snr77N34y0HhYqZQqJQ6RqU6FquYKVjK4mIgsGgkChMumaqdC1X0FKymUREQWBQSJRK6YiWOCUNEXWUlh/JLjAoJAqVU9lIrKU4gJUgEZGF9WENg0KiMNUqG90qHd3KEw7W9UTUGSsKB4NColBpWtmkJFpKyWYSkR+sKGoYFBKFSdvKRtdyBSsdW0lE/rCmcDAoJAqVppWNtsFqsDhPIRF1xHqihkEhUZhY2cSKe5+IOmNN4WBQSBQqTSublASrzmaKbqO/iUgfKakP3WBQSBQmbSsbXcsVNGX/b1q2l4i6x/rBwaCQKFSaVjbaBqtERBFjfVjDoJAoTNpWNrqWK1ja7n4i0ggrCgeDQqJQaVrZpCRaSsdWEpEvKakP3WBQSBQmbSsbXcsVLG13PxHpQ1XjLoE2GBRS9G59BdjdiLsUEdE0KklJtGTaAJMTP/05/OxnX467GEQpY1Y9ESYGhRStagX4fz4APP/DcZckGtrWNdoWLFAmTknzi3/wtbiLQJQuKfmR7AaDQoqWk6b/2ufjLUdkNK1sUlIJpmMricgf1hQOBoUUrZQEIzXabq+u5QqHac3IRBQhbevp6DEopIil7eLTdHvTUglWdtGL3bhLQWS2tNQXxKCQIpa2UV7aVqa6litY3/n7P4yv9P1o3MUgMpu29VhAkr59XWBQSNHixaeHlByHQ3f+Iu4iECVA0uuLpG+fewwKKVppyxTWKhvdRr+yEiQil5L+IzLp29cFBoUUsZRdfLXKRrPtTlklaMKUNCplx4RMkvRzM+nb5x6DQopWajOFutG1XOlV5SEhXSW93uYPshoGhRQtnS++T38I+IOfD3aZum6vpsUKiwlT0jBTSNpK/LmZ9O1zj0EhRUvnX5w3XwYWnw14obpWNh7L9fKvA8u/G2hJyMJMIekr4Sdn4oNe93JxF4Ao0XStbLyW64W/Y/199m5wZSEAZmQzKaV0rceConOyImLMFFK0Unfx6VqZ6lqu9Er6fZcMlvh6mxefg0EhRSttdz5dt1fXcqUYDwnpK+EnJy++GgaFFK3E/+LcT9fKxi5XSipDE6akqabkWJCBvJybRp3PJpU1XAwKKWIpu/h0rRh1LRcA/PGnga//97hLETmNjwilXsKDQpPKGjIONKFoMVOoCV3LBeC3nrH+BjiYxYRBHMwUkrY81dsmnc8mlTVczBRStEy48b31MrB+M5hl6bq9upYrxXhISFtJbz42qawhY1BI0TIhU/gLHwL+5VTcpQgZK0HdcPJqorjw2nMwKKSIGXLx7a4Fsxxdb/S6livFeEhIW55OToNOaF58NQwKKVomZAoDpWtlo2u50ot9CklbXuptns9Gim2giYhMAZgBcN1+6TCAi0qpYtjLEpFzACYBTNj/lgFcUEotd7tu6pKuFUVY5Urb9mrKhClp0nVEyCzMFKZFLEGhiEwDmFVKnax7rQDgRRF5opvAsNtlicgsgCtKqcv7XrsuImeVUgvet4w60jVTGFq5dK1sdC1XejFTSNpK+kAT1oc1cTUfzwN4uv4FO3ibs98LZVl2APmSUmpx3+cvAFgAMC8iE12un5JA12A1LEZV2P6ZMCWNCUWktGKmMC0iDwpF5AwAKKWWmry9CGDabWDmYVntMoEX7b9n3KybPNI1+Aq9+Vi35ktWgrqp8pCQrhLfp9CksoYrjkzhU7D68B1Q16fPbWDW7bKeFJHrzT5ft5zTLtdNXuhaUaSt+VjX4xAWA7bXiGwmpVPiRx9rmqyIQRxB4RSAlTbvFwE8HtKyrrpYZtcDXagLul58YZWrVpnqVkHqVh5bSMGb6Lq9dZgpJH0lvE+hSWUNWRxB4RjaB14rAAphLEspdVopNdnis6fsvy81e1NEzonIVRG5euvWLZfFo4N0vfjCKpem26tpsdIcFHLyatJW0jOFRpU1XHEEhW4CvrEYljUDoKiUutTsTaXUZaXUKaXUqfHxcZeLpANSmynUTbrKZcKErNqeKkRJ71NoUllDFts8hR24zRQGsix7VPIZsD9h+HS9+NinUA8hHQeBpj9G6uh6SIgSP/rYqLKGS9cf0EH262u7LHtOw3kAM/unqqEQaJspTNnk1bpWgmkLzutwoAlpS9t6LCBJ374u6BoUths8EvSy5mE9/eRyh89RIDS9+NIWjOhaCaY4U8iBJqSvhA800bWejkEcQeEy2jfpdho8EtiyRGQO1tNNmvYjpBAwU6gJTcsVVlCo6ebW40AT0lbSB5rw2quJIyhchPW84VYKAK6EvSwROQ/gOgPCiGl78aVt9LGu5WKmkEg7fMxdasQRFF5Bi0Cu7ukjL4S5rLonoRwICO1gkcKia0XB0cd6CGl/ZXTd3gYmlJHSKeHnprb1dPQiDwrtx8wV7RG/+80AWLSfXVwjIgV7nsCGpmKPy5oGMNYmQ3jY7baQF5pefOxTqIfQMoX6Y6aQtMVMYWrENSXNWQCzAE46L9iZvVbTwjxX995Zr8sSkSn7s883yQgehvWElLnuN4dcY59CTdjlEs3CpdCCQl2Pwx5tTxUiT9elQSe0QUUNWyxBoVJqUUSetgd6OM8ingRwcn9mz/Y8gGn7r59lvQirn+FUm+Jd6GJTqFu63vl0DVbDou1xCKtc+h/fqq7HhIijj1MjtsmrlVJLsJp43Xx2AcCC32UppQ65LiCFQ9fgi30K9RDSccgYEBRqe6oQcfRxaug6TyElliYX3/XfATZu173A0cdaCC04D2exQWKmkPSV8EyhrsmKGDAopGjpcPFVq8C/+yTwb/7G3mvMFLr4ShTbEtboYw3OOyJTJb1PoVFlDReDQoqWDkGSqlh/b71S91rKggZdRxOm7TjUYaYwAK9eAV75XNylSB5d64ugmFTWkMXWp5BSSoebfrVs/a0feRtapeAsV7NRvp4yhREcu9D6FOpf6fO+FIBfPWP9ffZuvOVInIT3KTSqrOFippAipsHF5wSFiCAo1PVO7+mXP4PCMDFTSNpK+rmZ9O3rAoNCipa2mcKUTV7ticnNxxqcdx0k6UyhhPFyXRoVaJlU1nAxKKRo6VBRVJygsP70DztTqMF210tZplB0OO86UAaUkdIq4c3HvPZqGBRStHS4+Jo2H6csU+jpl38UQWE4+4tPNCHywVNMaNIJbVJZw8WgkCKmwcUXZfOxthWjrpnCsIJC/ZuPU//s4z/9D8DK1+IuBTXFTGFaMCikaOnUpzCKgSa6Voy6TjER2rOP9Zfq5mOlgP/094FffCLuklAz7FOYGgwKKVo6VBRVe57C+j6FzBS6+IrBfQqZKdSbc9w33463HNQcH3OXGgwKKVo6ZQqjCAp1rRhTNtDEhEpf6XquRKGWvSc9adqyEBiTyhouBoUUMQ0uvlpQWP9iyuYp1OE4NJO24LyOtqdKFJzsPemJmcLUYFBI0dIpU5jq0ceaZgpD2l+iw3nXQarvS4pBodaS3qfQpLKGjEEhRUuHi6/WpzDFTzRJWZ9CbYPzOql+ogmbjzWX9HMz6dvnHoNCipYOGZumfQo5+jiU73S9jrAGmmh6HOroX8IQVTWoF6g1XeuLoJhU1pAxKKSIaXDxRdl8rMHmNpeyTKEBlX6qM4VsPtZcwvsUGlXWcDEopGhplSlkn8LuvmPu5NUmBIW6niqRYPOx3nSot8NkQv0QEQaFFC0dLr5mmUKOPu7uO6EFb2mep1DXcyUCHH2st6Q3H6f6F1kjBoUULR0qiignr9a1svGbKTSsmVd0OO86MKCI4Ul48/FuWf8fJe0lvPk41RdfIwaFFDENLj4++xi++xSGtr8MW26AdD1TIpHgTOGNtzfwyD/9TfzatdfiLop3zBSmBoNCipYON+dIn32sKWd7u9nuhqDQrOZjEyp9Nh8n01feXAMA/OYX34y5JD5w8urUYFBI0dLh4otjoElQ2/3pDwP/7m8Gs6xu1W+DYRm9jAF9CnW4NGKT8OZj8zFTmBa5uAtAKaNFptC5AUXZfBxQpXPzS9Y/v/w2BxkWFJpQ5yujbqIBC3r0sVKNP/pilIij6jdTWNq2jnHvUGBFClSar719mCmkiGlw8TWbvDq0cgWTKXzxy2/hxE9/LoDyOPxOSWNa87EGP0Y60ODKiE/Qzcc6/PhMFJ8/Ij/9QeDifcEVJ3CpvvoaMCikaOlQWdeCwrrXQp8fz9/yP/vnb/gvSz2/v/xN6/tnQCYg1X0Kg24+1qiPoh75Sp88Xe915/OdrwVWlFCk+drbh0EhRUuHiy/KJ5oElCkMvmkxBaOP6/aZCfMU6nBpxCbwTKE+QWEiDmtD1xGXW6T7Cd1QvpjKuvE2cPOVeNbdAoNCipZWmcIIRh8H3acwKGmYp7A+KNT9BoWUZwqDCArr959GmcJk8DJxvebncxR9pDv5hQ8Bn/72eNbdAoNCipgGFUUck1drd8P3O9DEgD6FUfSBpGAEMdCkPhDUKFOYjOZjD1k17eq8/SKozzpZfyue9bbBoJCipVOm0MTRx0HxnSk0ICis2+cZ3fZ/E6nOFAYRxNUvQ6NMYSKOaxStBFFLwnEJAYNCipYOF2Kz5mPNRx8HT9fm45AyhQbcyLQ7RaIURBDXkCnUZ2eWKvqUxbsENh/rkCnUEINCipYOFx8zhf4zhSaMEq4rr+i2/5uo6l/E8ATSp1DP5uNKEg5sEpuPdRho4tAos82gkCKmQUVR61Nozujj4Gk6JU1Io4+1OO86SPXk1UEEcVU9m49LFf2z1J0FlCms6rQvNMoUVnbjXX8dBoUULR2a8aol62+ngRNBVBQBZgoDfVRbGp5owuZjcwSSKaw/3voEheVEZAo9XEvNNtupe3Wg049GBoWUWjrc+crb1t+G5qYmFV1oWStv8gjyUWAp6FNYt40mNB8rA8oYmqBHH2uUKUxGUOglgGryOY2CH70yhfoEywwKKVo6ZGxKdlDYqWN6EGUNMFOYQ4A3Ot+ZQgNGH9f3KYy70nchCbGDZ0GPPtYpU5iI5uM6fiav1ij4YaawOQaFFDEN7nzlLetvp4ETgQQoquGPH4EGhdpmCgMMPJVhmUL9ixieIPqa1Z+TGvVdK9ujj8XkCQsDyxRqFBRqlSlkUEhppVOmsFOQo1GmUAHII8ApN3R9okmQFbVhk1cnYj47rxI8eXXJDlDNPrwerstmn2OfwubKDAoprXSoGcvNmo/D6lMY3OjjXH2fQt83UZ+jj0ObkibAQK7hMXca/BjpQIMrIz6BNx/rc7wrSZin0NMPQvYpdE2j/cKgkKKlQ2VdcpqPI+hTGBClgLzUBYJ+KxFPMWHETzTxvf81ygS4kOopaYKevFqjgSalJHQWDWqewkqQg+V80ilTyKCQ0kuDCtJ1plCfKWkqVdXYfOy7b46mzcdBrsOwgSaN3Sn1L2+ggmg+1nVKGnugSVmjfo7d8zFPYX09q1Hwo9WUVRr1tWRQSNGK++ID6jKFEfQpDKj5uFSpNg408XsT1bVPYaBBoUaZABfq+xQmIbnUlSDOJ00zhc4TTcomNyP7ma2gPuDRqU+hVs3HO/Guv04urhWLyBSAGQDX7ZcOA7iolCqGvawg101d0qFeLDcZaBLW6OOAMoWlSlWDTGHEU9IEminU4MdIBwczhSYPV+1S4I+50+d4O88+NjpT6Oe6rM8Oeq23PvNDwLu+F3jsh719vxmdfjRqlEGNJSgUkWkAs0qpk3WvFQC8KCJPdBOcdbusINdNHuhQWTuZwiiajwPKFJarat9AE59BobZPNAln9HEm7krfhVRnChM9ebXdfGxyptBX83F9X2gP9Va1ArzyWetfkEFhvdgzhfpkUONqPp4H8HT9C3YwNme/F+ayglw3dU2DijHKJ5oElCncLe9rPg4sU9hFuUxrPjZsoEm91D3dJMGTVzuZQqMHnPgZaNKQKfSQEdtc6f47bsSdKdS0r2XkQaGInAEApdRSk7cXAUyLyEQYywpy3eSRc6OP85dZ00xh2JNX+88U5iWEPoVdlSviKWmCbD42IMiqzxTGnbiIXODPPtagRcJWsTOFFZObj/1kChv6FHqot9bf7P47rsTcp7Bc148w5ZnCpwAsN3tDKeW8fiakZQW57lSoVlWwIyHrMmcll49/8lIGpRSqrX6ZN+tT2CEorC+DUl2UJ8A+hXFmCitVhWrFfZ+t7dLeZ3fLVby9frAj9f59WKkqqH039t1yteXnO9rX3N3qfNguVdoue6dcwdZupfX51HTVqjbAwP13mv93s2UHrdLFNaaUws3Vbdef39wtd77WA5mSpv6JJvpkCp1m4/rm43bn0lur2/j1P3v9wOvdnk9B2S5V2gbcdzZ2sb7TJNjzmSmsnV/rb7n+jlvbpQqq1ca6pVSp4tf/7HW8tbod+Pqaqt8XZX0GmkjUUx+IyHUAy0qp0y3evwNgUSl1NuhlBbXuU6dOqatXr3YqnmfqtWt46/d/Bd+8s4WtUgXlShWZjKBcURjpz2F1u4zXi1uYODKEwd4svnpzHfcW+tGfz+Lm2jZEBEO9OeyWq9itVLFTquDtjV3kcnlMjA/h/rEBfPWtVVQrZQz25PDm6jbWt8sY6sthuC+PnXIFI315bJUqWN0uoVSu4qHDg+jLZ3FrbQeFgTz68hkoBby1CZR2tjA+3IueXAarWyWUKwrb5Qq2q1nc3djGaF8WAz1Z5LMZTO58GSd2/xIA8O/Vx1FVCqP9eeSzGQz35ZDPZlDcLGGkP4eVDeui2SpV0JfL4vhIH0rVKkrlKoZ681jfKWFj19q2/nwW5WoVgz05DPXmcGdzF1sqj3eO5fH1lU2sbZfQm8vgkWPDmHprobavP4OP4/5DAzi29VW8c/svGo7D7/ZNQ/UMYmywB1+7vQEBUBjowZ3NXfTkMjhxeBBfX9lELmsdm7tbJfTlMnjo8CDyWcFrd7bwreU/xz27N7ArPfhvo9+Ht9d3MTbYg91KFb05ax/msmLdNBWQzQh68jnc2qjg2GAGPbkMdspVfPmNVdwntzGd/RMAwNeOfQyHj96L14tbWN0uI5/NoFSpQimFclWhL5fFdrmCt1Z30J/P4p7RPtzdKmGzDBwa7MHht5fwzspXAQCfL3wSw315lCtV3N6soAqBVEsY7Mnh/kP9WNsp48tvrOKhzG18uHoNAPDno9+Db+wMolK1KtN7RvuxsVtGfz6L3lwGX3lzDWNDPRjIZ3FjZRMA8MixYZSrChkBipslFO39VVEKQ705rG6V8KGht/DQ+p9a+3/g47ixWsXYYA+UsoKLwd4cjo/0obhVwsZOGSP9efTmMtjcrWC0P4+NnTLWtst4YGwAI7KFyTc+CwD4w8qjeAUPIavKuOfQIFa3dlHoy6K4WcLmbhlHhnox2m+d80O9OVSVwspmCSXksb6xUbtB3VvoRy4jGBvswVapglJFoVKt4vZm1cpGVivoy2dxe30HuUwGk0cHkRFBLiN44+42NnbKOHFkEJu7Zdxa20FvzrouMrk8bq1uoVI3j9u9hX4UN3fRm8tic7eMod48cllBcXMXAz1WGXfKVQz15tCXzyKXEWzsljFeGEGltI2313exXapgfaeM4b48Bnqy6Mtn0J/PYrtURRkZVCpVlMpl3N0qQQQ4MtSL3XIVo/15a3BTNoPN3TJ6shmMj/QBuV4sv/E2VjZ20ZvLoC+fxUh/Hj3ZDAZ7c1jfKaNaVVjfKdf24zfubKInm8GhgR4UBvJ47c4WBnpyGBvMY2Ongp5cBu8s/yWOrX0JAPDi0Pfj0GAvdspVrG2XMdyXQ6lilSmbEShl3dRvr1vXYX8+i9fubGKifx2nNv/AOqflcXytVMCRoV4M9+WglPWYuaOFYaxtbuPNu5vIZzNY3S6h0J9H1V5mfz6LR44PY2VjFzulKipKYXOnjPWdMu4t9OPO5i7uK/QDANa2y+jJZfDwkUHceHsTm7tlq+9vRrC5W8FupYr7Cv345p0t3Nm06rKxwR6s2dfrO48OYadcQVUBKxu7GOrNob8ni1feXIVSwD2j/ejPZzDQm8PGThlfu72BYyN9GO3Po1xVGOnLobibwer6OipVhUMDPVBQ2NqtYrgvh7XtslWf5ARbu1VUlUJPLoPenJUL2tyt4M7mLo4N9yGfyyArgtsbO8iKIJsRjPbncXerhNtrO3is7w08WrLqx9/rP42dTJ/VpSWbwdsbuxjsyeKRY8P4+somPrH5GwCALw59BDsD92Cguob33P4tAMBbR78Dr+wcwdp2CfcW+nFzbQdHhnrQk83g1ZvrKAz0YGygB2+tbeOttV1INofvKtzGieJ/BwB84cjfxO2NXdzdLOHwUC8KA3nslCpY3bbOt56+fvSoEtZ2SshmrPtRTy6DXEZwaKAHuazg9voOvrGyhV61gydzvwcA+FL+UfzZ7v21Hy+HB3vQk8uiqhTubpWxVgKODmasayiXRTYjtX3cm8vg9eI2hnqzmDw6hK/d3sDmrnVeb+yUce/hUaxubqNUKuPYSC9eL25b53O+go9u/FcAwJ0n/k8c+s5zrmMEL0TkmlLqVMfPxRAUtg28OgVufpblZ90icg7AOQB48MEHT964caNT8Tz7/H96Du/9k2f31g1x3ccoI3ZDmV0JZuwHblaqVRyS9YbPFjHUtEN7NiPW62qvb5OIHMgMDGMLeamgggzuqoG974sgKwpDylrfFnpRkl5UlUJVKYzJOm6pEWRRrZU5I4KKnYFrtq569fsjI1JrdstlrMquXK1iENvolTKqSlDEYMP3K8hgsXIS09lryKJaW9+YrOOuGsCobOKOGnK1x52y1JdjvzFZRxEjqNq/sPcfz4wIRKxER5/soh/Wr8aiGkR13wjUMVnHlurBFnqal0cEAqspUiD2sbT2e0aAAuxjonrQL7sNx6EPJQxI83U751EBa9hCL7ZVT227Vd15Ui8rAshehsMpm7Obslnr/+/PIo3JOu5gBAJlPy9WkLGL4nw2l8lAxGpW33u8tHXuZOrWmRWFUWxgR+USBeQEAAAgAElEQVTRK43Z1U3Vix3pqZ1r9ccBAEaxjowo7KgcNtC3t4/t4yf2PsmpMobF6pKwrvqxi2zTY5Oxb7bONgis9fWqHfSLFTTcUUP2+di4T0UEWdk7lvXqj/mIbCGHCsoqg43MEIC9jKVzXSulMCC76IO1zhU11HDt7C+ziKBSraKADWREYVdlsZ0drGUX6899pyzOdwCgJ5vBbt0xzmYEgFjZ97ptHJN1rGAUoqrW/nX2s8v6T0QwgG3kUUZRDR54P4sqRsX6gbKperEN68douU2W1DknAOnY/CsiDXVmuzph/zZ1U8c7RrGBrCiUkMMG+qHQmEXOyl6nCee9rF3P1m9fVR2si/ZnJQWCQ7KGuzJS2w/OOePU2wAwgk3kxHp/RQ3Vvp9HBcOy1fBaOyKCQ1gDAGyrPPqkdKBO3l92Z90llcWGDLTe9/a1VKkqVAEclrXaNdBshPiYfd/cUXlsoBdi14v1xyubkaaZ3BwqGKnVDX3Y3Te217lP/eGjz+IjZz/Vecf44DYojGP0ccHFZ8ZCWpbndSulLgO4DFiZQhfL8ezR6f8Jv33vx/GxR49jpC+HXDZjncBKYcNO05cqCq8XtzDSn8c9o33oy2dRqlStGzGsSiCb2bup3767Bvz8/Q3rGfln30DZ/pXqLDubEUyM7124VlpdIZ8VXL+1gXsLfShVrEpg89+fxehrv4vsve9H+W//Jl59ax2T40M4PtoHrL0J/F/vAgD0ft9F9H/gfwZgNcUtvb6Kdx8ftrIdVYWMXc7dchWv3dnEw0cG8Y2VLdxc28bx0T7cf2gA1arCbsX6tduby+Jzf/EG7h3tw7c9UMDvvHITH33XOHpz1s14u1RB6f/9CfS+/BlkRu7B6Kdexu31HWzuVjA+3AsAeOi1IspHBrG6W8Gkvb3fWLGyB7m+HAaygnwmg1K1CoGV6Xn97haODPWiL29lbla3yhgf7sWmnSHbKVexU65ipC+H2+u7GO7LoT+fBTKCAoC7WyUM9FgZHaWA3375TTz24CEcG+mr7Wtc+9fAZ/8hAGD97/0xMDSOI0NWmVe3SviDt9Zx/6F+PP1vr+LVm+voyWbwn//Xj+D4SB8ODe4FitYNF3ZmxTqGPZVN4OJ9AIDME/8Ef3LiR/Ct9xdQrlZRrQKVL/wisHjeOjee+QqWVwWr2yUcGezFA2P9tSCoH0CfUnjtzhbuP9SPb6xs4e2NHbz3vlFkRfDqzXUcGshjzN6u7VIFGRHks9Z2V5SVTXGWd3NtG2MDPdjYqeB3vnITDx4ewGMPFGrv17v6VyvIZARTDx6qbacVKAG7lSoGe7IQEWyXKri1toOR/jzQn8fGZ34Cva/8KlTvCGRnFQCQ+9iz6P/wj2O7VEV/z14gt12qWNfOz78bWH8L8m0/iOtTP4NH7x3Bxk4FhYE8dstVDPZa52/lxn8D/s0nAAADf+czGJj4KL5Z3EJfPlvLWueygqPD1v7Y2q2g184A9/dkUfr/LgKf/znrmP/kK3jgyDAA68a8vlPGaH++YR9s7JRRVQpvre5guC+HYyN9tebq3Gd+EHj1t1A59l6M/sTvA7Ca9+rPDWtH/jLw2U9BDRyB/C9fxmB/HhkB3ri7jXtG+/Dm6jYODfSgL2/tl9/+0pv4zv/6MfSv3UD2vZ/EyNlfhlIK3yxu4b5CP1a3yrizuYsTR/aCsd/+0pv4xp0t/NhHTmBlZQu31nfw2AOF2vXu7IvX725BAJR7czg60odq1foxILIXQK9tl7G+XcZAbxaVqsJgbw6bu1ZdeHiwt6GuA4DRqsIffvU2PjhxGDfXtjHQk8PdWzcw+q+te+Lm9/wMCt/5NHLZDHbLVYgAb97dRjYjuPH2Jr71/lGUKlUUBvb22+LLb+Hv/VurheiBsX585tyHcHezhN/5yk18cOIwTj60d04653ipYv0tVxV2yhVsl6q4ubaNdx8fwdp2CQLBm6vbeJednezLZ3BzdQeHh3pQrQJ/dP023nF0CCeODOKN4jbuKfRhfbuMz796C9//R2eAmy8j/+5PoPCDv4pqVeHW+g6++M27+Mg7jtSOHQAUN3exuVvBfYV+lCpVbJUqyIpgsDeH14tbGO3P442725gct47f6nYZGbHqj75cFmvbZZR6sxjNZnBrbQd3NnfxyLHh2vJXt0v4g1dv4+Nf+DHgG3+Eux/437DxwZ/C0ZFe5DMZfH1lE9u9Ofzs517Gj374BN7/QAG31nawZmf2j4304p7RfuyUK7V6XP3sMUh5G2vf9mP46rf/E7z0VyuYHB/Cdz0yjtXtEoZ7c9iwryUAyPzKx4HXvoD8O78HhR/+Nbzy5ioqVYXN3QqOj/RhsDeHl/5qBd/1znH092Rrx+mLb66hVKnisQcPoVJVVrm2S8hnMxgf7kX15yeR2S6i9yN/H3c/+E/R15NFRgQ7pQqG7Ra10f48tnYr+I0/fx3jw72YeuAQ8jnB9T/9fbzvv3zSug4/8WnceeCJ2n77/F/ewrdPHMbbO2V8aGTvR2fcYpunsAM3wVtYywpy3Z4cHurFD37gwYbXshlBFtJQSTkBjiOfbd1F9MhQ/4HXRKxf8cdHW5+QIlbzAwC84+i+X3n99volg6PDfbWbnv1i7b/qbwK9uWzthr7/vZ5cphaQPnh4AA8eHmj4XF9mr5L7G992b+2/P/bo8YZi9eWzQG+utpHZjNQCL8eHJ48c2NYHxgYOvNZbt877D+29P9CTw0CPtY7hPuumnctmMGjvkmb7tP7mLgJ84r33NLwvIkD9+saGgMG9dfblszhqb8ch+zyYeqiA99wzcmBd9fu1dgyr9cchh8fs45B11jlQt79zvXjH0ebZSGeZzv7af6zedXy44bP1NycRILMv++mcN6MDGXzysftarhMATp1o/M2WyUhteT25vfO/L59tOJ5jdmAt2b1j0NPbB4g0BIQN5RXrb08+h8ft9TrH3LnWMhlBJrdXjWayGSAjDevu72m89pz1OX/z2b31PzC2F1Q5TXj7DfY2nneAdTxyWQHEKldvXZkOBIQAMPqA9b19799rN43eM9pY5o89ehx4MWeXK1Nbp3NNjA7kMTqQP/gd2/5zxNHfk639KHPUn7vOfh4b7MHYvu0Y6m19+8pmBN/1yDiAvet2rLK3b48M9wL2sp3zxjlmzj7Yb/pbjuGhwwO48fYmHnvgEO4r9OO+Qj++5d7G66/xnLT/ZgQ9uQyG+/bqbacud/abs30njuxt1/e+b6+OcPbfocEe/MD77wP+2DlP987FYyN9B+o6Z10Fe/fns5mGe4WzvfX1+/7zrv7Yjg/3Hrj3jPTl8X3vuwe4KrXvj9ZdA86Phf/7Bx+rvXZ0pA9H95XTCQgB1DJy48N9GL9vFO+9b7RhfcC+cyDj1PnWMt59/GC9+PG6c9I5TvXLzWYEx0f7GuvvzF59cLRu3zrrds6f/p4snjz1QMP63le37AfGBvDA/XvhhXN9NLvG46TrE02CnCuw22Ulc55COXioJePz8DvLbLJsNGR54piEV/b9NYTUBSjtjo84f7rYvvrj1CQLh2zdTTerV0UVmPp9kOmwjeLyHKo/TtK86bjDipqs06N212S90fvbv9902YZeUzXB1Ul+D1MgnEJ0OtZxCGIHdXu+uT33u+UEmxkPObTY74Pd0/BsAgAEOTFRt8sKaVKkmIVRi9UqpWbLDvBG50XbsmmsIXBrfXnWqsuuNq9DBVUfCJq23zpqcoPJts6ENny2075oCOQ9BIUNh8VvUOgyUHCCwkf/h24W3rgO04j/OkmrLdeyjguyLF1un/PjzMs12LYY9vK8BIVx3wc9iKP5eBntm2jHADSbRzCIZQW5brOEEhQyUxi4+gqtTdbJ0/2gU8CZ6z34WlI022GdsqFuMw8uj1nH9QTBbZl7h4FPvQwM7W/Aa7dsQ6+pmgAzhf4KEoywMmOBiCNT2NicHhinHMwUhmYRQLsJogsAroS0rCDXTW2zErpkCqNftS8NzZttgkJ7w7prPu5wTDpmzkzmIVMoB/6jxefcZXddrMi/bgKF0fu67CagY2aqC0FkCnXadh2DwkD3T7eZwmzj36B56m5lXqYwjrPpCloEZnVPE3khpGUFuW5Cm6Aw9l9IhmY14swUJrUfIdAiU9gpKHR503XbD7RT2YIQZpMiM4V739biBq/x8QikT2GXQa80DhALHDOF4VBKLQAoish0k7dnYM0j2DDYQ0QKInJORBqafrtdlpd1UxuuL1b2KXTNZf+0vc3rKips/3bqMoWdKnmX55Df5uO4MoXdL9z+Y9g15Uhcn0INM4VBcpupd4SeKfTbpzCwkoQqrrPpLIDZ+hfsTN0ZWMHZfs8BmLP/+l1Wt5+nVtinMHgNAUa7gSZO83EXGm6KzTKFCQ4KQ80Uumvyd/V9v2plZqbwIPYpDF2czce1/aFRUBj7fbB7scxTqJRaFJGnRWQOwHX75UkAJ1tk6p4HMG3/9bUsD+umVtpWSnH3KQzx5him+gqtTdm9NR936lOY4OZjT30KXQZBfvsUBtp8HEWmMIRFRyGATKFWvzW1bg2JY6BJyEGypwykeX0KY5u8Wim1BJeZObvZd6HN+66X5eXz1EK7SinuX0g6z+HVTpcVj+c92zRTyNHHjd9x+cMiLc3HxmcK6yVgG7T84RtmhrqDWvNxSHV+WBlIzRh2xySt6Jwp1OonfRdc3sydvoTeO7wzU6jPPIUhBIWhnPc6Z6ZcCLBPYVej/sOiY/OxI9Dd4zZTyIEmQdDwbCJj6NynUOumlTZcBhW+Q960TUnTbIqijk808TJPYcxT0rSbESCgRZtyczvIf52kx6hjW6Iyt21oMyVNOiavZlBIPrSrlOK+GAytMF3+yvUd86ZtoElNN83HUfUpDHKgSZjdJgz9oeUIok9hMF8PiM5dZEL4odPxY8wUBkHHs4lM0TZT6HdCX59Skin0XtGkrPk4snkKNWk+Zp/CJvyXW6st17H5ONRR750+F9Jj7hwpGWii0dlExmmXlYj9AjD0Buaygs/U+hQGuJ5EZwpD7FNY37E9yQNNEpUp9Ld/tNgDWg40sYUxIXsnmZCDZN8/+DQ8Tk0wKCTv2mbjYv6FZGqmsNvmY8/rSWufwm6ajw/8R4vP6ZgpDDNjY9g1VRPAQBOdNl3LGRZCHODU8WMhPfvYwT6FRB20uwFp8wvJjAuxxvV0CmFkChPcfNw0UxjGQJOYp60INXtk6A8tR4B1kha7QMfm45o4MoUaDjTR5j7ono5nE5lC5ylpUpMpDLBPoWn7yotu+hS6npKmvvnYwz40ZaCJ6ZlCv8cJfq63EIQ6/ZBHYf4Y6fixsAeacJ5Cog7c9ink6GPXup2SxnNMaNh+8atZQNPpl7/bAEvL5mP2KTwowEyhFvWKjs3Htjj6FIY90MRvsGnIdRPbE00oAdr+Uo07U+jcHKNftS+RTUnT4ouHHgbe/0MeF6ozP5nRTgNN+EQTIwQxebVOm67lQJM4M4VOpcjmYz8YFJJ3JkxebciFWOM6UygNf7vX4nv/4E89Lk9z9VH0+LuBW6908d1OmcK695kp1FjS+hRqeDzinJLGEdZj7lIy0IRBIXmnc59CnZtW2ukyU+g9JjRsv/hWt8N+7LeAtTe7/2rL9+vPdU0mrw71GbRm3NwOCHDyai1oOfrYEcIPHdefZ6bQDwaF5J3reQo50MQ1l79yQ5mSJsnqz4f+gvXP9Xe7uCklufmYmcK9b+uwC7QcfRxj87GDk1f7otPZRKZxOyUNH3PnnttMob1dGa/7VqsbSRT8nA9dfCfu5uNQn31s6DXlCKRPoUb7QMfRx444BpqE3TrEyauJOnDdfMFMoWtuKx7fm2fYfvHLz/nQVaYw5KCz46KYKWzNvBt0ezo3Hwepy2PFyat9SfrZRKFyeZNgptA9lxWa761L/I1kPx97LOzz17gnmhgqiEyhv68HS8vmY0ccmUIV3DqbSUmfQh3PJjKF6+YLZgpdcz3QRBr+dr8ew/aLX362N+ybbqADTZgpbMnAG3RbOk5JE+fk1WHzlCk0D4NC8s7tDSjWTKFhup282ut6tMwuhMnP+RD2uWRaptDQa6uez3kKtdoDOl7LsfQpDImyM5Ce9rM0/U+daXg2kTFM6FNoGrfNx/6jQq9fNJOX86F2MzCw+ThMpl5bDfwFhVrRKihMcLcFDjQh6sBtViKWSsuMC/CArjOFbD52J6KBJp4YklVJVKbQ59d12AVaH48QRtTHxdd8kBxoQmniuvk4/KIcXKehp7bLPoXOVDR89rFLngYFR3TTNaapzfA+hQ28DjTRaNujymR3I9YBThrtBwczhZQubm+acTQfR7/KQHQ7JY3nFZm6g7zSeaBJGM3H7FPYls9jqlVwqNUP4DgHmoQ8+tgTZgopTVw3H3OgiWtdTl7tPVOYskvf1+hjAweahHKDTFCm0OdAEy3oPMOCMdlvF971fdbfnqHuv8tMIaWK68CCA01cc9unsJYpZJ/C0HGgib3sBGUKk/CYO4dOP/CSOCXNJ34O+EevdPdozBpmCilNXHf1iDNTaMaFWNPt5NXMFLqUsilpQqFxZqpbPiev1oqW13KCMoXZHDByj7fvMlNIqWJEplDHfiZtuKwA/bccmVFBBSYtk1eHeVyZKQzo20HTqzTBS/r26YVBIXnn3MxUh8CLfQoDJ34zoVpmF8Jk76dO52rTr7L52F64/ScB11aSngSk47VszDkdNjYfU5oYkSlMJt+ZwoTvnwN0zhSG0nzM0cft+e1TqNE+0CpoinNKGg2x+ZjSxWVkwkxh4Pw/bivZ++cgjfsUhpJV4ejjthLVp1DHUgVZJh23zwMtj9NBDArJO7fNx8wUhoBT0nQlNZnCKM77JFxbSdgGm07XcqyTV2uImUJKFbeVETOFgeOUNN3SeJ7CIG/qfPaxOz7nKdRqF2hVGFugZdJw+1xjn0JKE9c1ZByZwmSf2s4ezTBT6I6vCtmg5uNImFbeZpKwDY4kbUsTxl0fdZgppHRxeZLHcVGbXJG4sBePe44KAyuLGVLSfBwq51m78ZYiED77FPIxd50wU2gxL1OYi7sAZDCtH1RuxgXY1N/6JeDII20/4vumZEgFFRidH3Nn3CPBknDueG0+1nDbdQoK4+xTWHjQ+jt8PPgyeGVgppBBIYUvjkpLx8rbrfed6fgRTknTLY0zhaYdC9PK24zPbdBqF+hVmBAW6XKZH5gBxiaBd54OvgwpotFPDEosDjQJnO9mLJ2yC1EI+xzsHfWxT007V00rbzNJmpJGw2s5joEmmQzwyMf0DZK1KldrzBSSf5ySJnJOMxYfc+eWj9Sqm5vuP37V29NSAPPOVdPK24zf3hfBlMIf53zTKSjklDSNDCw7g0IKn3nXhfZ8T15tYGXlS9h9CnO93pcf5AXiNTDtShLOHX9T0uhFy0IFyOTtMy9TqNFPDDKWllPSmHEBeuU0G3Py6giwT2Ej08rbjM9jqsUu0HLSxBCYvH0GDjThnYH865SdYJ/CwHFKmiiZNHk1Rx+74nlKGo22PZKssEecvNrGTCFRE8wUBk32/e1+ASm79NPymLsoJOLa8jv6OAn7IAzsU9iAmUJKpU4XLTOFgRO/UaHJFa0nKZmnMBKmlbcJr/vcb1/eIGl93gQ596bJYQozhURNxJEpTPapXRt9zClp3GGmMDiG3NzaS8A26Nh8HMq5kYBjZZCU3RkoFnzMXeBqiUJOSeOSn+01qE9hJBJw7vh8zF0SdkGojHtKT0gMbD6OfEoaEZkCMAPguv3SYQAXlVLFsJclIucATAKYsP8tA7iglFrudt3UDTYfB45T0nQnLY+5i4Jp5W0qAdug5XFgprCRec3HkQaFIjINYFYpdbLutQKAF0XkiW4Cw26XJSKzAK4opS7ve+26iJxVSi143zJqi5nCwPmfkibZ++cge3u9NLmx+Xgf08rbhNdMYe3HWAL2QaiYKQRgZKYw6naLeQBP179gB29z9nuhLMsOIF9SSi3u+/wFAAsA5kVkosv1k2vMFAbN/80p2fvnAF83FmYKG5hW3qaSsA0aYp/CfczLFEYWFIrIGQBQSi01eXsRwLTbwMzDstplAi/af8+4WTd5wExh4Hz3KTSuH5tfGjcfG3fTM628TficpzDh1YsPnJKmATOFbT0Fqw/fAXV9+twGZt0u60kRud7s83XLOd1uhSJyTkSuisjVW7duuSwmWZgpDBofc9clrfsUGhagJ+Lc8dd8TB1w8mqLMFPYzhSAlTbvFwE8HtKyrrpYZtv+jEqpy0qpU0qpU+Pj4y4WRzXMFAZuL2PBKWnc0XhKmmSfqnryWT/wkEUoMXW5GdsR5Z1hDO0DrxUAhTCWpZQ6rZSabPHZU/bfl1yum7oVSwBixgXoF6ekcUnnPoWmHYtE3KSZKQwXB5ocYMh2RHm3dhPwjcWwrBkARaXUJZefp64xUxgWTl7tls6ZQtPOVdPK24TfTGECdkEoONDEeJHPU9iB20xhIMuyRyWfQYf+hOQTH3MXGk5J45LOfQpNO1fTdu7U4VQ0LnHy6ibM2A7d0gVdT2DtdVn2nIbzAGb2T1VDCZCYiqQ9763Hul36YXNG5njYY6FnCk07Fgm4tnz3KUzAPggFM4UtGXJPalsbici8iCiP/655KE+7wSNBL2se1tNPLnf4HPnFgSah8f7003Tsnxqd+xSadq6aVt6m2KcwXMwUHmTGdrRtPlZKnQ1wXcto36Q7BqDZvIOBL0tE5mA93YT9CCPB5uOg+a4nE1PRhsh5+gmbj/cxrbxN+MzO8vJpgX0KWzPkpImy3WIR1vOGWykAuBL2skTkPIDrDAgjxEyhfoxrsvSLA00CY1p5m+GUNOEKtE9h2uqqeEW5t6+gRSBX9/SRF8JcVt2TUA4EhHawSKFgppBi5qkvYe25MYEW5eB6TLvpJeHaSsI26IhPNGnJkO2IrDayHzNXtEf87jcDYNF+dnGNiBTsJ4k0NBV7XNY0gLE2GcLDbreFusRMoX6MC0RiFPq+MuxcTcK15fUxd7VHCSVgH4SKTzQ5yIztiHpKmrMAZgGcdF6wM3utpoV5ru69/f0bXS9LRKbszz7fJCN4GNYTUua63xxyh5lC7aTtpuZrSprgitF8+QGu4IFvt/6++68Ht8wDknDuJGEbNBRGvZKUusqQ7Yg0KFRKLYrI0/ZAD+dZxJMATu7P7NmeBzBt//WzrBdh9TOcalO8C11sCnWDmUL9pC5TqPP5EGDZjj0KPHs3uOU1k4Rry2umcN9f2o8DTUwX+eTVSqklWE28bj67AGDB77KUUodcF5BCwEyhflK2f3QOZHQuW1OmlbeZJGxDShh3fbRixnakLV1AcWCmUD/MFOrDtGORhGuLj7kzSEJ2tiEnjWG1EZkpjqCQp3ZbhlRQgXG2V3mf7pscSTh3krANOgvwOktMXWXGdvDOmSbvC3IucgDD91h/Cw+1/xyffayhtO0fjbfXtJueaeVtho+5CwcHmrRmyHYwKEyTT/4rvPw9vxTc8t7z/cAP/RrwwZ/o8EE2H2snbZlUrc8HncvWjGnlbSYB23DsUevv6P3xlqOZQDPyCThWAEzZjsgHmlCMsjlUeoaCW54I8M5mU0Xu/1wcAYgZF2BstA6SwqDx9pp2LEwrbzNJ6FP4oZ8ETnwHcN/Jzp+NDDOFLRmyHSlLF1AszR4caBI43z/EmSnUR9qOhRY0Ph/cymQ0CwjrMVN4kBnbwdoobVLTv8+MC9Av71uZjv2zx8P2OpF36GNTDDsWOgfYbvHZx+EIZZpC7u0oMShMGUlL1o71SHtpq2h1fjyZjmVqy7TyUuTYp/AgQ65zBoVpk5qmKjMuwNgYUkEFx8P21gLJYEvSZEVhryBYSTh3ktCnUEvsU9iaGduRlgiB4pSW7GQMOOteAhj3Qy0J15bfoDAJ+yBMQc5TaNr10YIh50xC9ja5FU9llpZ1RseQ+kUfOu8wncvWjGnlbSYpgYZuQjk3EnC+ATBlO3hlpE1asnZJuHFRgHQ+H3QuWzOmlbcJ1g/mSMqxMmQ7GBSmTGoyhcwEUD2dK2Sdy9aMaeVtKgnboDEONGnCjO3gnTN10pK1M+MCpKjofD7oXLZmTCtvExxoEhIONGnJkO1gUJg2aQnQDLkAKSI6nw6mZbUTcW0lYRt0xkzhAYZcN4bVRuRXap5okpSKhALi4Xz4rn9s/T38jmCLsp8hN4s9ppW3Cd+TVydgH4QhjHPZuOvDbAwK0yaW6ystgagBPvFz4Qc5OvJyPrz7rwHP3gX6RoMvj8kScW0lYRs0xj6FxmJQmDISR1MVM4X6+OCPAz95Le5SxEDj88G4IMu08jbBPoUhYabQdAwKUyclARorEqqn8/nAPoUxSMI26CzIyat5rKJkWG1EfsXz7OM4TjNWJFTPPh8CbdYKCs/VyPnuU0hNcfJq4zEoTJu0NB/z1yXV0/l80LlsicV9Hqogf3zx+ogUg8KUSc3k1VpmhCg+Ot9YdC5bQrFPYUiYKTQdg0IKH2tQipvO56DOZUss7vNQcEoa4zEopAjwog4aE6Hd0vgcNG2gSRJwnkKDcF9HibVRykgmJf37+OuS6jnng5bnhY5lSjg2H5uDOztSDApTJyUDTYhMweuDkibQgSYMU6LEvZ0y8Qw0oaDxMHZL5x2mc9mIusGBJqZjUJgysTQfE8VN5yha57IRecIpaUzFoDBlFH91USppfN6zecw4bHFpgZNXG4+1Ucpw1BylktY3cZ3LRuQBJ682FoPClGHzMaWTxuc9b3rG4RFrhZlC0zEoTBlmCimVdA68dC4bkSfMFJqKQWHa8AKjVDLhvDehjASwGm2JfQqNl4u7ABS1lFxgfOQH1dP9tJ9+FnjHdNylIAoG+xQai0FhyrBPIaWT5uf9d3wq7hJQFzQ/m2Lk7Bn+KDcVm49TJyXVGVnHUxQAABZzSURBVH9dUj2eD0REHTEoTBkONKF04nlPweE8hZRUDArThs3HlEbOTZx9TYnCw2DZeAwKU4aZQkonnvcUHMY+HfDHl7EYFKYMB5pQKvEuTgHi2dQK94zpGBSmjPCQExFRqALIFOYH/S+DusYpaVKHv+QoxZgxpCDwPGouyP3yI78OvPyfgYGx4JZJHUUeFIrIFIAZANftlw4DuKiUKka9LBGZAHBBKTXT7bqNxeZjIiIKUxB9Cg9PAt/5j/wvh7oSaVAoItMAZpVSJ+teKwB4UUSe6CYwDGhZ8wCW3W8BERGlHX9aU1JF3cFsHsDT9S/Ywduc/V5kyxKRcwAKXa7TeJxfi4iIwsXRx6aKLCgUkTMAoJRaavL2IoBpuzk39GXZGUUghVlCkZQMNHG2M5OPtxxElDj8bd0Cd4zxoowQnkKLIEwp5bx+JqJlnVNKXXa5LgBWZlFErorI1Vu3bnXzVa2k5pq9dwr4yD8A/tYvxl0SIqJ04TyFxooyKJwCsNLm/SKAx8Nelp1lXHS5nhql1GWl1Cml1Knx8fFuv66RlESFmQxw+l8Ao/fFXRIiShg+BKAV7hfTRRkUjsEK1lpZgfs+fp6WZTcbT7Rodk6HKEcff/f/Dhx7b3TrIyIiIs+iDArdBHxuJyTyuqxzSqlLLteRSJFOXv3dF4Af/8Po1pcik+NDAICHj3CCVyLXJBvMYpgQ64DNx6bSbfLqIEcDNyzLnsImvRlCG0cfJ8MPvP9ePHxkEN96/2jcRSEyx6e+CKy/FXcpkov3F+PpNhS16wmsu1jWlFKq676EScOgMBlEBN/2QIHHk6gbI/cC9z7mezG86jrgQBNjtQ0KRWReRJTHf9c8lKfd4BHPyxKR8wC6Gm2cWKzNiIgoFLzBmK5t87FS6myA61pG++bhMbhv3nW9LHu+wmUvj9FLokj7FBIRJRAT9J0wU2iqKPsULgKYbvN+AcCVEJY1BeApEXmqyeemARRFxHkCysXEj0xmbUZE5AunpGmB9xfjRRkUXgFwrtkbdU8feSHoZSmlFgAstPjsdQBLAWdEtcY+aJRq7Ovk3dik9bf/ULzlIP3xOjNWZG2JdnBWtEcB7zcDYHF/E6+IFOwniTQ0FXtZFtkYFBKRF6f/OfBDC8ADH4i7JPFjNdoCd4zpou5gdhbAbP0LdmbvDKxgbr/nAMzZf/0uq5Ugp8HRHjOFlGo8/73L9QLvPB13KcgIzBSaKtJ5CpVSiyLytIjMAbhuvzwJ4GSLzN7zsPr9PR/AsgDUnmryHKy+hhMAJuyR0stpaEZmXxgiIn9Yi7bAH13Gi3zyansgh6tMXrv+gN0uq+47RVhZxlSSKB9zR0RERMbg/CSpw6CQiMgPdsPpgK3HxmJQmDIiPORERBQGBsumY4SQNrxmiYh8YTXaCVOFpmJQmDJs9iAiolDw/mI8BoUpw6CQUmnwqPX3kY/HWw5KBFajHXDyamNFPvqY4sU+hZRKw8eAn3oVGDgcd0mIEozRsukYFKYNf+JSWg0djbsElBCsRjthptBUTBulDCevJiKiUBw6Yf0tPBhrMcg7ZgpThpNXExH5wx/XLbz/fwRG7wMe/mjcJSGPGBSmjDA5TETkC5uPWxABJr477lKQD4wQUsbJFO4o/h4gIiKiPYwMUuhi6W/jd6rvx2/HXRAiIiLSBoPClBEB5irfH3cxiIiISDMMComIiLrAhwCQaz+5BJR34i6FawwKU4aj5oiIiCJyeDLuEnSFA01Shj9wiYj8YTVKScWgMGVYmREREVEzDApThn1hiIj8YTVKScWgMGVYlxEREVEzDAqJiIi6wAF7lFQMClOGzR5ERETUDIPClHH6FP71b70n5pIQEZmJP64pqThPYQot/bPTGO7joSciIqI9jAxSaGywJ+4iEBEZi4lCSio2HxMREXWBzceUVAwKiYiIiIhBIRERUXeYKqRkYlBIRERERAwKiYiIusE+hZRUDAqJiIiIiEEhERFRN5gopKTiPIURUkphbW0Nq6ur2NzcRKVSibtIlHDZbBYDAwMYGRnB8PBw7Yk2RERE+zEojIhSCjdv3sTGxgbGxsZw/PhxZLNZ3qQpNEopVCoVrK+v4/bt29ja2sLRo0d5zhH5xGuIkopBYUTW1tawsbGBhx56CNlsNu7iUAqICHK5HAqFAoaHh3Hjxg2sra1hZGQk7qIREZGG2KcwIqurqxgbG2NASLHIZrMYGxvD6upq3EUhMh7zhJRUDAojsrm5iaGhobiLQSk2NDSEzc3NuItBRESaYlAYkUqlwiwhxSqbzXJwE1EA2KWQkopBYYTYOZnixPOPKBi8lCipGBQSEREREYNCIiKibgiHmlBCMSgkIiIiIgaFREREXWGikBIq8smrRWQKwAyA6/ZLhwFcVEoVo1iWiBQAzAKYALAMYMz+zlK36yciIiJKikiDQhGZBjCrlDpZ91oBwIsi8kQ3gaGXZdlB5DyAGaXUYt3rc7CCS0qRpaUlTE1NxV0MIiIiLUTdfDwP4On6F+zgbc5+L7RlicgEgBcBnN0XEJ4BcM4OMilFzp49G3cRiIiItBFZUGgHX2jRTLsIYNoO3MJa1hyAF5p8ZxnAkv2XEu7SpUtYWrJOgYkJ6xRZXl7GpUuX4iwWERFR7KLMFD6FFoGXUsp5/UwYy7KbjadhBYb7P7+klDpZ972mROSciFwVkau3bt1yWUzSzfT0NC5cuIBLly5hamoKly9fxszMDJuRiYgo9aIMCqcArLR5vwjg8ZCWNQO0zCy6opS6rJQ6pZQ6NT4+7nUxFLOpqSlcuXIFhUIBly5dQrFYxJUrVzA9zd4DRESUblEGhWOwgrVWVgAUQlrWNOzMop3xO2//nbWziJQSS0tLOH36NIrFIs6fPw8AOH36dK1JmYiIKK2iHH3sJuAbC2lZEwCWROQcrH6FtYBSRK6JyEWl1ILLdZPBFhcXMTs7i6mpKZw+fbqWJVxcXGQTMhG5o+IuAFE4Ip+nsAO3mUIvyyoAWG4yVc0FAPMisuhlrkQyi5MdBICrV68CsJqUGRASEVHa6RYUBhmU7V/WRP1UNA6l1KI9v+E5ALEOQf3nv/ElvPz6apxFCNy33DuC/+P7Hw1t+cvLy5idnUWhUMDy8jKeeeaZ2gCS69ev1z43Ozt74LsvvvhiaOUiogTjE00oodr2KRSReRFRHv9d81CedoNH/CyriM5Tzrgd5EKaWFpawtzcHObm5jA7O4uZmRk88cQTWFhYwNjYWC0QrJ+Gph6zg0RERHvaZgqVUkHO7ruM9s3DY7DmCwxjWZ0CwiKsfoexCjOjlkQXL17E/PzePOUTExMoFou4cOFCLUu4vLzM5mEiIiIXomw+XoQ1CriVAoArIS3rKoAnO3yek1cbZHFxEU899VTDa0428MKFC7XX6oNGIiIiai3KoPAKrH57B9Q9feSFkJY1D+tRdoX9g0ns/oQA8JLLdZMGms0r+NJL1iF88sl28T8RERE1E9k8hfaUL8UWzxieAXBg9K+IFOz5BBuairtdlj3AZBnNA8knARSVUnzOmeEWFhYwNTWFQiHIQexERETpEOXk1QBwFkDDMFA7s3cG9lNH9nkO1qPpngtgWWcBPFP/TGT7v2ft98hgxWIRy8vLfDIJERGRR5FOSWNP//K0iMwBcOYLmQRwssUcgc/D6jv4vN9lKaWWROQJALMisoK9ya07PveY9PfCC1ZvgccfbxxEvrhozULEYJGIiKi9yOcptJ8/3CyT1+yzCwBaPmmkm2XVfZ5ZwQQ4efIkJiYmagNJnOln9o8ynp+fx9zcXOTlIyIiMk3UzcdEvi0tLWFpaQljY1ayd2FhoRYULi/vJX0vXLiAmRnXvxmIiIhSTbcnmhB1NDU1hfPnz+Pw4cO4cOECJicnce7cOVy7dg1zc3O4cuVKbb7CiYnYp58kIiIyAoNCMlKzx9ZNTU2xqZiIiMgjNh8TEREREYNCIiKirqi4C0AUDgaFRERERMSgkIiIqCsSdwGIwsGgkIiIiIgYFBIRERERg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiKi7qi4C0AUDgaFRERERMSgkIiIqCsSdwGIwsGgkIiIiIiQi7sARN0qFot4+OGHPX33ySefxNzcXMAlIiIiMh+DQjLSxMREw/8fGxsDACwuLtbedz6zsrJS+9zk5GREJSQiIjILg0IyTqFQwLVr15q+Nzk5ieXlZczOzuLMmTMRl4yIiMhc7FNIRERERAwKiYiIiIhBIRERERGBQSERERERgQNN9PKbPw28+RdxlyJYx98HfO/PxV2KBsViESdPnsTU1BTm5+exsLCAixcvolgs4vr16ygWizh06BAmJiZw/fr1hu9evnwZMzMzOH/+PGZnZw8se2FhAXNzc7h69SrGxsYwPT2N2dlZFAqFrsroDJYpFApYXl7GM888g6mpKVy+fLmhTM3KQERE5AUzhZRKy8vLWF5exszMDM6ePYulpSVMTU35WqazrKtXr2J6ehqFQgGXL1/GyZMnu1rO0tIS5ubmMDc3h9nZWczMzOCJJ57AwsICxsbGaoHgpUuXsLS05KvMREREDmYKdaJZRi3plpaWsLS0hNnZWZw/f97XshYWFnD58mVMT0/jypUrDa+fPXsWly5dcr2OixcvYn5+vvb/JyYmUCwWceHChVqWcHl5GVNTU74DWSIiIgczhZRq58+f9x0QAsCFCxcAoCGYA4AzZ85gYmLC9VNUFhcX8dRTTzW85mQDnXU462k1VyMREZEXzBRSqj3zzDOBLGd5eRkAcOjQIV/LmZ6ePvDaSy+9BMB6RB8REVFYGBRSahUKha4HgDTjBIRTU1OYmZnxvbz9FhYWMDU1FUhZiYiIWmFQSKnlPC/Zr/rnMJ87dy6QZTqKxSKWl5cDaeImIiJqh30Kibqwf4oax9TUFJaWlmpZw6C88MILAIDHH3+84fXFxUUsLi4Gui4iIko3BoVE+zjNtCsrKw2vF4tFXL58uel3nGlizp49i2Kx2PDe8vJywyCRdk6ePImzZ88eWO7+Ucbz8/NN+x8SERF5xaCQqAlnGpjTp0/jwoULmJmZwcMPP9yyydmZpHppaQmHDh3C2bNnMTMzg5MnT2JychILCwsd1+lMkeOsY2FhoRYU1mcgnfIQEREFiUEhURNzc3OYmJjA4uIiFhYWsLKygvn5eczNzWFqagqTk5MHvnP+/Hlcu3YNZ86cweLiYi2rODs727LZud7U1BTOnz+PyclJXLhwASsrKzhz5gyuXbuG+fn5WjA4MzPD+QmJiChwopSKuwzGOXXqlLp69WpX3/nyl7+M97znPSGViMgdnodE3v3IL38Bv/eXt/Arf/dxfM+7jsZdHCLXROSaUupUp89FPvpYRKYAzABwUieHAVxUShVbfyuYZYnINACnw5bTDnhRKcVnhREREVGqRRoU2kHZrFLqZN1rBQAvisgT3QSG3S5LRM4BgFJqZt/nnxOR55VSnTt9ERERESVU1H0K5wE8Xf+CHbzN2e+FsiwRmQBwUil1ef/nlVJnAcx2uW4iIiKiRIksKBSRMwDQoql2EcC0HbyFsaxp7DUxN7NkN0UTERERpVKUmcKnADSd2Vcp5bx+JsRlPdVmeROtlucQkXMiclVErt66dctlMYmIiIjMEGVQOAVgpc37RQCPt3nfz7JeADAlIlfsfoQ1dt/Eq536MyqlLiulTimlTo2Pj7ssJhEREZEZogwKx2AFa62sACi0ed/zsuyA7yysZuQ7zqATZzRy/eATIiIiojSKMih0E/A1f1xEAMuyRxefhBVMzonIdQATDAiJiIiI9HuiidtMoZ9lXbb/TcAKDt32YyQiohT7mR94L773vcfx4cnDcReFKBS6BYVdT2Dtdlkich7AtFLqgp0dPA1rcMm8iHBKGiIiauvBwwP4hR8+id5cNu6iEIWibVAoIvMiojz+u+ahPO0Gj3help0NnFRKXXJeU0otKqUmYWUNzzv9DMPERwpSnHj+ERFRO22faGJP7ByUZbRvHh4D4PZxc90uaxZWf8IDlFIzIgJYj8u73OwzQchms6hUKsjlIn+yIBEAoFKpIJtlhoOIiJqLsvl4EVY/vlYKAK4EvSxnCpoOU85c6LA83wYGBrC+vh7mKojaWl9fx8DAQNzFICIiTUUZFF5Bi8Cr7ukjLwS9LDsYdDOq+arLdXsyMjKClZUVVCqVMFdD1FSlUsHKygpGRkbiLgoREWkqsqDQnhKmaM8NuN8MgMX92TwRKdhPEmloKvawrMUOo4zPwXpmcmiGh4cxODiIGzduoFgsolwus48XhUophXK5jGKxiBs3bmBwcBDDw8NxF4uIiDQVdQe3s9jXv8/O7J2BNRp4v+fq3tvfv7GbZT0N4EURcQJK1H3nHPYNQgmDiODo0aNYW1vD6uoqbt68yawhhS6bzWJgYABHjhzB8PAw7P6zREREB0QaFCqlFkXkaRGZA3DdfnkSwMkWff6eh/UUkuf9LMv+/yftrOM8Gkc5zyulQhtgUk9EMDIywiY8IiIi0o6wCbN7p06dUlevhtoFkYiIiCgQInJNKXWq0+d0m7yaiIiIiGLAoJCIiIiIGBQSEREREYNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiIiIgKDQiIiIiICg0IiIiIiAoNCIiIiIgIfc+eJiNwCcCPk1RwBcDvkdVD3eFz0xOOiHx4TPfG46CeKY/KQUmq804cYFGpKRK66eU4hRYvHRU88LvrhMdETj4t+dDombD4mIiIiIgaFRERERMSgUGeX4y4ANcXjoiceF/3wmOiJx0U/2hwT9ikkIiIiImYKiYiIiIhBIRERERGBQSERERERgUGhlkTkjIhMi8j5uMuSNiJSsPf/GRGZrXv9wDHhcYoHj4teROT/b++ObhOHgjAK/yOlAIvt4KYDlO3AJUQpAUpZbTqAFtgOQgdZ0kFcAqKDyQPjYCFr1y94rsT5pEi2FSnWPQImxiTL/jEzOEaXRIO1Xo0co8mMYn3fro5NapHRh6GwMv0Tq7vvJZ3MrE0+pXvzImnh7n8kycxWY03olCPWucQ2XeqwjsdLMbNCl1yxtl2sddcP7RJNMsT6fpvaIqsPQ2F9fkrqYruTtEw8l7vj7lt37/88QJG013gTOs3MzIouay7RJV1ciTqYWXH3V3fvRJdsfyXtzGwpqbj7h2hSk6ktUvowFNanudr/kXIWdy4GkGO8yI01odP8SvTo0SXfY3wdzWxjZo3oksrdT5I2knY6t5FoUpOpLVL6MBTW5yRpkX0S0LO7r2N7rAmdZmRm7fXbMKJLLT5jEDlIWokuqeJtx727Pw72aVKPqS1S+jzM/QPxX++6/IZQJL3943txA2b27O6vsd1qvEkzcgy3c4wWjc73ri1Flxq86/LC1ej8Qja8uk6X+ZX+nmhJv3S+T5rHSj2mtkjpw5XCygxu2G5j//rqCG4o1v23mR3M7CCNN6HTvNz9I9Z4oXiipEu+WO9msN5buqTbxgfkWkkvNMkVV2qfBh8cmdQiqw//5g4AAABcKQQAAABDIQAAAMRQCAAAADEUAgAAQAyFAAAAEEMhAAAAxFAIAAAASV+ycBzz8940QgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 8))\n",
    "plt.plot(x.value, label=r\"$x^*$\")\n",
    "plt.plot(x_true, label=r\"True $x$\")\n",
    "plt.legend(fontsize=26)\n",
    "plt.yticks(fontsize=26)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-200.,    0.,  200.,  400.,  600.,  800., 1000., 1200.]),\n",
       " <a list of 8 Text xticklabel objects>)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAHJCAYAAABHfXcUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3WlwXOd95/vfv7HvjY07KRKAJGu3IciKt9ixwMSphJNJTEqZyc1NzWQEJrlVdzK5KTJ+kVszU5NxwGTuTE0lUyZUWWqSjCNRTm6i8STXpLzEVmyLixZbsiyS4Cau2IkdBPDcF+c01AAbQAM43ecc4Pup6gJx+vQ5Tz8tsX98nuf8jznnBAAAgLVLhN0AAACA9YJgBQAAEBCCFQAAQEAIVgAAAAEhWAEAAASEYAUAABAQghUAAEBACFYAAAABIVgBAAAEpDDsBqxHDQ0Nbvfu3WE3AwAABOT06dO9zrnG5fYjWOXA7t27derUqbCbAQAAAmJml7LZj6lAAACAgBCsFmFmTWbWYWbHw24LAACIB6YCMzCzVkltkpKS6kJuDgAAiAmCVQbOuTOSzpjZ/rDbAgAA4iOSwcrMkpI6JR13zr24xH4d8kaVBv2f3UvtDwAAkEuRClZmdlTe1NtJSR2SFl3fZGadkvqcc0fSt5lZU/o2AACAfIlUsHLOHUz92Q9OGZlZk6QO51ztgtcfNrMBM+tyzg3msKkAAAB3iVSwWoGDkk4s8ly3vNGuI5Lkr5N6JotjPs80IgAAWIu4Bqt2LR2s9soPVn5YIjABAICci2sdq1ZJfYs81y2vVAIAAEBexTVYLScZ0HGyrmHlFxM9ZWanenp6Ajo9AACIk9gFK78Ug+SVWMjVOZr8xfMHJbX6Vxt2LPUa51yXc67NOdfW2LjsPRoBAMA6FNc1VjnlnOuWdDjsdgAAgHiJ3YhVWhmFoKb7AmNm+8ysa2hoKOymAACAEMQuWGUhVYk975xzLznnOmpqanJy/FvDE3r35nBOjg0AANYursHqhKTmRZ6rk3Qqj23Jm//68ln9fNd3wm4GAABYRFyD1RlJTYs81yTpWB7bMifXU4FlRQUam5rOybEBAMDaxTVYHVWGWlX+FYOtkl7Ie4uU+6nAsuJCTdyZ1eysy8nxAQDA2kQ9WGWsI+VftdeV4X6CnZIOr9f7BJYVFUiSJqZnQm4JAADIJFLlFszskKQn5I06SVKnmR2QV019XmDyb7jc4b9mUN6aq+Nh3u/PzPZJ2tfS0pKT45cXe8FqfGpG5cWR+ugAAIAiFqycc0dWuH9XrtqyGs65lyS91NbW9mwujp8asRq/w4gVAABRFPWpQKQpSxuxAgAA0UOwihFGrAAAiDaCVYByXW4htcZqjBErAAAiiWAVoFyXWygtZsQKAIAoI1jFSDlrrAAAiDSCVYxU+CUWLvSOhtwSAACQCcEqRnbUlunh7dX669euht0UAACQAcEqQLlevG5menxXrW7ensjJ8QEAwNoQrAKU68XrktRQWaLhiWlNsIAdAIDIIVjFTENViSSpb3Qq5JYAAICFCFYx01jpBave4cmQWwIAABYiWMVMfWWxJKlvlGAFAEDUEKwClOvF65JUVVokSRqemM7ZOQAAwOoQrAKUj8XrVaVeLavRSRavAwAQNQSrmKko8YLVyOSdkFsCAAAWIljFTHlRgcykEaYCAQCIHIJVzCQSpsriQg1PEqwAAIgaglUMVZYWMmIFAEAEEaxiqKKkUKNTBCsAAKKGYBWgfJRbkKTKkkLKLQAAEEEEqwDlo9yCJJUWJTR5Zzan5wAAACtHsIqhwkRC07MEKwAAooZgFUOFBaaZWRd2MwAAwAIEqxgqTJimCVYAAEQOwSqGChKMWAEAEEUEqxjy1lgRrAAAiBqCVQwxYgUAQDQRrGLIW2PFVYEAAEQNwSqGChKmmRlGrAAAiBqCVYDyVXm9sICrAgEAiCKCVYDyVXmdNVYAAEQTwSqGuCoQAIBoIljFECNWAABEE8Eqhrw1VlwVCABA1BCsYqiQESsAACKJYBVDBayxAgAgkghWMVSYMDknRq0AAIiYwrAbEFVmdsj/Y7MkOecOhticeQoSJkmanp1VQaIg5NYAAIAUglUGZtbpnDuc9vtRMzvunNsbZrtSCv1gxYgVAADRwlTgAmaWlNTq/0w5KqndzJpCatY8749YEawAAIiSSI5Y+aGmU9Jx59yLS+zXISkpadD/2b3U/ivQJqlJ0hn/927/ZzLz7vk1N2LF/QIBAIiUSAUrMzsqqU7SSUkdko4vsW+npD7n3JH0bWbWlL5tpZxzg5JqF2xu9392KwIKCryBRkasAACIlkgFq/QF4n5wysifkutwzs0LQM65w2Y2YGZdfkAKyuckHQz4mKvGGisAAKIpUsFqBQ5KOrHIc93yRruOSJKZ7Zf0TBbHfD7TNKIf8I4657pW2dbApV8VCAAAoiOuwapdSwervfKDlR+WVrXuyg9l56MUqiRGrAAAiKq4XhXYKqlvkee65S0+XxMza5ekVKgys2RUrgosZI0VAACRFNdgtZw1Xb1nZq3ywtsZM2vyf/+cpP4gGrdWjFgBABBNsZsKTKsvlZOF5P7xX5YXzuYtoE8vGprhdR3y1nZp165duWjanLk1VpRbAAAgUtbriNWqOecGnXO1zjlb+FjmdV3OuTbnXFtjY2NO21hU4DVlaobF6wAAREnsglVayYNIFOtMZ2b7zKxraGgop+dpqCyRJN26PZHT8wAAgJWJXbDKQqoSe945515yznXU1NTk9DzbkmWSpOtDBCsAAKIkrsHqhKTmRZ6rk3Qqj23Ju/qKYhUXJnRtcDzspgAAgDRxDVZn5N3LL5MmScfy2JY5+ZoKNDNtrSllxAoAgIiJa7A6qgy1qvwr+lolvZD3Fil/U4GSVFSQoPI6AAARE/VgVZdpo3OuW1JXhvsJdko6HJV7+uVSwiRyFQAA0RKpOlZmdkjSE/JGnSSp08wOyKumPi8w+Tdc7vBfMyhvzdXxTPf7yxcz2ydpX0tLS+7PJZMTdawAAIiSSAUr59yRFe4fqXv4OedekvRSW1vbs7k+l5nkyFUAAERK1KcCsQgzE3e0AQAgWghWMeXd1YZkBQBAlBCsApSvcgveucSIFQAAEUOwClA+yy0kzORYZAUAQKQQrGLKxIgVAABRQ7CKKTNjhRUAABFDsIopr9wC0QoAgCghWAUor4vXRR0rAACihmAVoHwvXp8lWQEAECkEq5ii8joAANFDsIopb/E6yQoAgCghWMUU5RYAAIgeglWA8rl4PWHGHW0AAIgYglWA8rl43bulDckKAIAoIVjFVIICoQAARA7BKqYYsQIAIHoIVjFlZpRbAAAgYghWMeVVXidZAQAQJQSrmOKiQAAAoodgFaB8l1tgjRUAANFCsApQXsstiFvaAAAQNQSrmGLxOgAA0UOwiinKLQAAED0Eq5hKWNgtAAAACxGsYsrE4nUAAKKGYBVTiQSL1wEAiBqCVUwxYgUAQPQQrGKKAqEAAEQPwSqmKLcAAED0EKwClM/K69wrEACA6CFYBSifldcTJs2SqwAAiBSCVUyZmRyrrAAAiBSCVUyZUW4BAICoIVjFlInF6wAARA3BKqYSxuJ1AACihmAVU8bidQAAIodgFVMJFq8DABA5BKuYYsQKAIDoKQy7AVFkZklJT/u/NktqkvSsc24wvFbNR+V1AACihxGrzDolnXLOdTnnDkvql3Qs5DbNQ+V1AACih2CVWZOk9rTfz0tqC6ktGXETZgAAoieSU4H+VFynpOPOuReX2K9DUlLSoP+ze6n9s+Wc27tg0xOSTqz1uEFKmDFiBQBAxEQqWJnZUUl1kk5K6pB0fIl9OyX1OeeOpG8zs6b0bQG0ab+80HYgqGMGwcTidQAAoiZSwco5dzD1Zz84ZWRmTZI6nHO1C15/2MwGzKwriIXmaSNiZ9Z6rKAZI1YAAEROpILVChzU4lNz3fJGu45IcyNOz2RxzOcXTiM657r8Y3RIumBme6JyZSD3CgQAIHriGqzatXSw2is/WPlhKet1V/76rs9J+nxaiDohb+SqfSXHyiWvQCgAAIiSuF4V2Cqpb5HnurW2K/iaJB2St9YrJen/jMRolZRaY0W0AgAgSuIarJaTXH6XzJxzZyQdcc51p21+RtIZ59yiVwaaWYeZnTKzUz09Pas9fdYSCQqEAgAQNbGbCvSn6qTcjh593swOpf2elPTUUi/w12N1SVJbW1vOIw8jVgAARE/sglU++GurVlyywcz2SdrX0tISfKPuPhdrrAAAiJjYTQWmLShf9XRfrjjnXnLOddTU1OT8XN5VgUQrAACiJHbBKgupSuzrmnevwLBbAQAA0sU1WJ2Q1LzIc3WSTuWxLXPMbJ+ZdQ0NDeX8XJRbAAAgeuIarM7IK4uQSZOkY3lsy5x8TwWyeB0AgGiJ6+L1o5JOL9zoXzHYKumFvLcoz7xb2oTdCgAAkC7qI1Z1mTb6Naa6MtxPsFPS4bBuO5PPqUDzf7KAHQCA6IhUsDKzQ2Z2zMzO+5s6zey4mR1Nq18lybvhsqTz/ms6/JB13Dm34jIJQcnnVGDCzD9nzk8FAACyFKmpwJWGotRNkjciP1dp1jkl5savAABAmCI1YoXsJfwsxYAVAADRQbAKUF7XWPlDVlwZCABAdBCsApTvcgveOXN+KgAAkCWCVUyZWLwOAEDUEKxiam7EilVWAABEBsEqphJMBQIAEDkEqwDlt0Aoi9cBAIgaglWAQlm8nvMzAQCAbBGsYipVbsHNhtwQAAAwh2AVUwkWrwMAEDkEq5hK3cRmllwFAEBkEKwClM/F64lEqo4VyQoAgKggWAUor4vX/Z+MWAEAEB0Eq5iaW7zOGisAACKDYBVT3CsQAIDoIVjFFPcKBAAgeghWMUW5BQAAoodgFVOpqcA3ruT+CkQAAJAdglWA8llu4aFt3pWHL5y6kvNzAQCA7BCsApTPcgsPb6/RLzy5S9/t7tP0DPe1AQAgCghWMfahXbUanZrRlYHxsJsCAABEsIq15sYKSdK5WyMhtwQAAEgEq1hraqyUJH2nuy/klgAAAIlgFWs1ZUXa99g2/dG3Luj42zfDbg4AABsewSrmfm//oypMmF451xt2UwAA2PAIVjFXWlSgx3Ym9eZ7g2E3BQCADY9gtQ586r5Gnbk8qFcv9IfdFAAANjSC1TrwLz++R4UJ09d/eCvspgAAsKEVht0ArF1FSaEe25nUn7xyUcMT09r/+A49tjMZdrMAANhwGLEKUD5vabPQ7/zsw3qyqU5//t1L+tn/9or+6FsXNDvLDZoBAMgnc44v36C1tbW5U6dOhXLu4Yk7+rW/OKNvnu1Vc2OFfu1TLWp/YLNqyotCaQ8AAOuBmZ12zrUtux/BKnhhBitJmrgzo//55nX95+Pv6urguAoTptZ7avVLH9mtzzy8RQUJC61tAADEEcEqRGEHq5TZWac3rw7pK2/d0N99/4Yu9I5qV1252u6p1fbaMt27uUqPbq/RjtoyFRYwKwwAwGIIViGKSrBKNzPr9JW3buh/vHpZ3T2juj40rtQSrNKihBqrSpQsK9aehgrdt7lSj+xI6tHtNaqtKA634QAAREC2wYqrAjeIgoTpJx/Zqp98ZKsk6c7MrN65Pqy3rw/p7M0R9Y9OqW90SqcvDehv37g297qGymI9uiOp9gc2q213re7dVCkzphIBAMiEYLVBFRUk9MiOGj2yo+au50Ymp/XqhT5194zqhzeG9Q9ne/TVd7waWVtrSvWj9zbqk/c36mMtDaopY1E8AAApTAXmQBSnAtdidtbpysCYvn2+T994t0ffOter4YlpFSRMH9qZ1I/e16gP7kyqtrxYzZsqVF5MXgcArC+ssQqQmR13zu3Ndv/1FqwWmp6Z1etXBvWNd3v0jXd79L2rQ0r9Z2Qm7awt1/1bqvTAlirtqCuf+72O9VoAgJgiWAXEzPZLOuacy3ph0XoPVgv1jUzqfM+o+kcn9e7NEf3w5rDeuX5bF3pH5xbIm0mbqkpUX1GiPY0VenBrtR7cVq1Ht9eovrIk3DcAAMAyWLweADNLSqoLux1RV19ZMheOPvPw+9snp2d0c2hSl/pH9drlQb03MKae4Um9+d6gvvzmdUlSYcK0u6FCu+rK9dHmen1oV1IPbatRaVFBGG8FAIA1iWSw8gNNp6TjzrkXl9ivQ1JS0qD/s3up/VfhaUkvSDoa4DE3jJLCAu2qL9eu+nJ94t7Gec8Njd/RD67f1lffuaUr/WNzf5a8sHX/lio9sr1Gj+1MqnVXrbYlS1VVykJ5AEC0ZQxWZvaCpD05Ome3c+6ZRc57VN4I0UlJHZKOL3YQM+uU1OecO5K+zcya0retlpm1Szqx1uMgs5qyIv1IU71+pKl+btut2xN6470hvXFlUG+8N6j/9b3r+suTV+aeryguULK8WFtrSvXw9ho9tM2bTtzTwIJ5AEA0ZFxjZWY1yt0UWL9zbtm7FJuZk3Qg0wiUmTVJOu2cq83w3ICkPc65wbU00sz2O+de9EfPBlhjlX/OOX3v6pAu94/p2uC4bgxNanB8Slf6x/TWtdsam5qZ23drTamaGiu0p6FCTQ2V2t1QrpqyIu2sLdem6tIQ3wUAYD1Y0xorP/gsG35CdFCLjyZ1yxvtOiLNLT7POEK2wPOpEGdmHc65riAaitUzMz26I6lHdyTvem5m1uli36jeuT6sC70j6u4ZVXfvqP729Wu6PTGddgzpgzuT2lRVorqKYt27qUrbkmXa01ChHbVlqihhpAsAEJy4fqssNU3XLWmv/GDlh6Ws112ZWaskhpsiriBham6sVHNj5bztzjn1jkzpvYExDY3f0bfP9+n1K4O62DumkxcH9MVXr8zbv66iWM2NFdqe9ELWlupSNTVWamuyVM2NlRRABQCsSFyDVauk5xd5LjVitVp1klr9NVaS1CxJZnZIwS+OR8DMTI1VJWqs8q5S/NT9m+Y9f2t4QrduT+rsrWHdGJrU5f5Rnbs1ojOXBzUyOa3+0am5fYsLEtpRV6aGyhI1NVSouqxID22r1q66cm1LlmkzU4wAgAXWFKz8tVhPOef+aol9npK3yPz1tZxrhe6eO8qSc+6E0kbD/BGsjuUWxPtXKHZI0q5du1Z7euTYpqpSbaryFr9nMjo5rYt9o7o+OKFvd/fpxu0J3Rya0Ikf3NTtiWlNTc/O7fvg1mptrSlVdVmRqkoLVVdRrG01ZdpUXaLy4kKVFiVUVlSg0qIClRUXqK68WIkE91kEgPUsiBGrn5f0V5JkZs9KOizptKSjzrmvOudeNrM9ZvZzSwWwbPmLySWvxEJOpa/P8q9YPOYHr7v4a7K6JG/xeq7bhtyoKCnUQ9tq9NC2GrU/uHneczOzTj+4flu3hid05tKgvn9tSNeHJvTDm8MamZzW0PgdLVVvt7y4QC2bKtVQWaLWXcm5+l01ZUXalixTUUEix+8OAJBrawpW/iL3pxds3itvDdQRM9sjrw7UcUlNazlXGFaxPmufpH0tLS25axRCU5Awf6SrRp/+wOa7np+antXN2xO6NTypiTszGp+a0cS093NsakYXekd1vmdEV/rH5mp2pZQUJrSzrlzJsiJVlhaqorhQDZXFemBrtZo3Veq+TVWqKWe9FwBEXdBrrLrllTp4TtJzfrDaL29U6y+DOIFzbtDMpDVM9+WKc+4lSS+1tbU9G3ZbkH/FfjjaWVe+7L7DE3f03sC4LveP6fb4Hf3wxrCuDo5raPyO+kendLl/TDeHJjSaVlKiqrRQlSWFun9L1Vwtrx+9t1E7asvk/z8BAAhZoMEqbdrv0/404AVJvxfkObKQqsQORFZVaZEe2FqkB7ZWL7rPzKzTtcFxnb01rLM3R3R9aEIDY1M6d2tEb743NHeFY1GBqba8WLXlxaopL1J1aZFqyop0/5ZKfWCLN+K1PVmWr7cGABvaqoKVmX1wscXofpi6sKZWLe+E/Kv1MqhTSOUSmApEkAoSNjcCtnDq0Tmnc7dG9J0L/XpvYEyDo3c0MDalofE7ujY4ru9fHdKXzrw3t//m6hLVlBVpe7JMjVUlSpYXK1lepNryYjU3VqplU6WSZUUsrgeANVpxsDKz35V0yMy+4pz7zILnAlmgnoUz8kouZNKkkO7tx1Qg8sXMdO/mKt27uWrRfXpHJtXdM6pTl/p1oWdUQ+N3dGVgXD+4PqyBsSlNpl3hKHn3aKyvLNY99RXaU1+h+spi7a6vUENVsWrKilVbXqT6yhJqewHAElYzYtUn6YCkJzI895ofvF7NccA6Ku/Kw3n8KwZb5S2YBza0hsoSNVSW6MN7Mt+dauLOjHpHJvXuzWFd7B1T78ikeoYndaF3VF/94S0NjE5pevbuyxwT5pWtaN5UoR3Jcm2vLdPWmlI1VpVoU5X3s66iWAWMfgHYgFYTrOqdc78n6UsLn/CnAX/LzD6bWme1xvZl/EZwznWbWZeZdTrnDqc91Snp8FrvEwhsBKVFBdpRW64dtZkX20/PzOra4IT6x6Y0MDalwbEp9QxPamj8jq4OjOtS/5hefueWekcm73ptQcJUX1E8V6x1k/+zsbJEm6tLtaO2XJuqvQBGmQkA60nGmzAv+QKzD8krsfAF59ylJfb7vHPucys89iF5I2Gt8qb0BuWtl+pWhsDkF+VMLVZvlnQyzMroaWusnj179mxYzQDyauLOjG7dntSt4Qn1DE+qZ2RSt25Pvv9nf3vvyJRmMoyA1ZYXqaGyRFtqSvXknjrtrCvXlupSba0p0+aaEpUUFoTwrgBgvmxvwrziYOUf/FlJX5BXn+q4pJcXLmY3s990zv3+ig++DrS1tblTp7jdIJBudtapf2xKN4Ym9N7AuHpHJuemH3tHJnWxd0w/vDl81+saKkt0/5ZKPbClWu0PbtYDW6pVXVZIiQkAeZWzYOWvoWqXN4qUGllKHeSMvBGmJnn31fvVFR18nSBYAaszPHFHN29P6PqQ97gxNKHL/WN69+aw3rkxPHdLofQSE8nyItVVFCtZXqwdtWXaUl3q/+5t31JTyqgXgDXLNlitqtzCwgP7NyzeK+mzkg5K6lzpNOB6QLkFYG2qSotUVVqklk13X+04NHZHJy/260LvqPpGvTVf/aNTGhy7o3O3RjQwNqXekam7XldUYGrZVKV7N1XOBa5kWZFq08JYU0MFI2AAArGqESvn3G8t8Xyr3l+DdXFtzYsnRqyAcIxOTqt3ZFL9o96C+76RKXX3juqta7d1sXdUA2NTGp6Yvut1JYUJba4u1ZbqUrVsrlTbPbWqKvVGvHbWerW/CF7AxpbLqcDPSnLLlVNYzeL19YJgBUTX9MyshsbvaGDsjgbHpnS+Z0Tnbo3o5u1J3bw9odevDN5V46ukMKHttWXaWVuuHbVlam6s1GM7k7qnvlwNlSUhvRMA+ZSzqUDn3Jf8cgqfl3R04aiUmf0rSf3y6l0BQKQUFiRUX1miej8Qte2eX9VlbGpaN29PamRiWj0j3kL79wbGdaV/TFcGxvTGe4MaHLszt/+9myq1u6FCm6q8UhKP7KjRrrpyVZUWqr6ihHpewAaz2jVWXzKzE/Km/J5b8PTTkp6SdGSNbQOAvCsvLtSehtRfjTUZ97nSP6Zzt0b0zo1hnbzYryv9Yzp1sV8DaYFL8q5ofGBrlbZUl+rRnUk1NVSosapE9f76LkIXsP6sqtzCsgc1e8o593LgB4446lgBG9v41IxeuzygW8OTGhib0qlLA7rqj3gtLKRqJiXLvBpeH9qV1H2bq3Tf5io9sr1GtRXFIb0DAIvJaR0rLI01VgDSOed0pX9c14bG1TPsLa7vG51S/+ik3hsY15vvDal/1LuisaQwoY801+veTZXeTbhry9XcWKktNaUqLqRKPRCWnJZbAABkz8y0q75cu+oz3z5IkvpGJvXOjWF98dXL6u4Z1TfP9s6rVG8mbasp032bK+dGtx6/p1a76sqVYEoRiAyCFQBEQH1liT7WUqKPtTRIkianZ3R7fFrdPSO61D+mqwPjutg3qh/eGNYr5/o0NeNduViQMNVVFKuuvFi1FUWqr/DuwVhbUayGymLtrC3XlppSJcuLtLmqlBAG5FjGYGVmL0jak6NzdjvnnsnRsQFgXSgpLFBjVYEaq0r0ZFP9vOemZ2bV3Tuqkxf7dXVgfG5qcXBsSj+4cXuucOpCVSWF2lFXru1J716MW5Ol2lZTpq01pdpVX66tNWX5envAurXYiNWzkuoWeW6t+nN03NBReR1APhQWJOamAxczPTOr/rEpXe4b8+7HODqld28M69rguK4OTujkxQENjc8PX9tqSvXgtmrVlhfrI8312lRVqoe3VytZzmJ6IFssXs8BFq8DiIOxqWldG5zQ9aFxvXtzRKcv9evszRHduD0xr0J9VUmhkhVF2lpdpp115dpdX67GqhI1VJaoeVOlNleXqKyogOr0WNe4KjBEBCsAcTZxZ0aX+sbUOzKpM5cG1D82pYHRKV0b9G6KfeP2xF2vSZYX6ZHtNXpsR1J1FcVq2VSp+7dUqba8mKsZsS5wVSAAYFVKiwp0/5Yq3a+qucX06SbuzKh/dEo3b0/o3K0R9Y5M6VLfqF67PKg//Po5pf973Uxqf2CzPtpcr+bGSrVsqtTWmlJGt7BuEawAACtSWlSgbckybUuW6UO7auc9NzPrNDg2pTOXB3Xz9oR+cP22Xnrjmo6/fXNun/LiAjU1VqilsVLNjZX6keZ6Pbi1WhUlfCUh/pgKzAGmAgHgfc459Y5M6dytkbmbXp/vGdH5WyO6NvT+tOL2ZJkev6dWH29p0P1bvMX5ZcUFIbYceB9TgQCASDAzNVaVqLGqRB9pnl86Ymjsjr57oU9n/Xsv/uP5Pv3tG9ckeTW6Giu9122pKdXHmuv19BM7VV7MVxeiixGrHGDECgBWZ3bWqbvXK4T6g+u3dfP2hHpGJnW5b0zdvaMykz7W3KAfaarTTz26TfdQeR55wlWBIeAmzACQG845fbu7T1//YY9eeuOabtyekHNScUFirtDptmSZdtSW6Wc+uE1NjZVhNxnrDMEqRIxYAUBuXekf09d/eEtXByfVc3tTAAAgAElEQVR0bXD8/Ye/ZqtlU6Ue2V6jqtJCba4u1Qe2VClZXqwP7UwywoVVYY0VAGDd2llXrl/8yO67tl/sHdXfv3VD3zrbq1OX+jU8MT3v9j4tmyr1qfsadU9DhR7YUqXdDRWqryim/AMCw4hVDjBiBQDR0T/q1dl669ptffnN6zp9aWDuJtaS1FhVoo8116ttd50+eV+jdtaVh9haRBVTgSEiWAFAdM3OOt0cntBbV2/rysCYTl8a0He6+9U7MilJenRHjX7ioS16eHuNPtHSwNQhJBGsQkWwAoB4cc7pQu+ovvL2Tf3d967rjfeGJEmlRQk9uLVav/TR3fqJh7aotIi6WhsVwSpEBCsAiLf+0Sl961yvXr88qG+8e0vne0ZVXlygT97XqMfvqdWDW6u1uaZUdeXFqikrYlRrAyBYhYhgBQDrx+ys0yvne/X337+h42/f1K3hyXnPlxQmtKuuXA9uq9YTu+v05J46tWyqZEH8OkOwChHBCgDWJ+ec+ken9Na12+ofnVLf6JRuDI3rQu+Y3nhvUD1+6KotL9KH99TpE/c26p99eJcKGNGKPYJViAhWALDxOOd0qW9Mr17s18kL/fr6uz3qGZ5UwqRtyTJ95qEtemBrtXY3lKt1Vy0jWjFDsAoBldcBACnOOX35e9f1/au39b2rg/pOd79mZr3v3A/uTOp//8g9erKpXtuTZSG3FNkgWIWIESsAwEJT07O6MjCmb5/v0x9+7ZyuD02oMGH6px/artryIn20uUE/9oFNYTcTiyBYhYhgBQBYyuT0jF4516v/cuKsbt6e0MDYHU1Nz+qnHt2qn39ipz68p04lhZR2iBJuaQMAQESVFBbo0x/YrE9/YLMkL2j9p6+8qy++ellffvO67t1UqX+z9z7dv6VKzdxQOlYYscoBRqwAAKsxcWdG/+9rV/U7X/6BhienJUm/8OQuPd22U9uSZWqsKgm5hRsXU4EhIlgBANZi4s6MzveM6M+/c0lffPWKJMlM+s0fv1+//PE9VIAPAcEqRAQrAEBQLvSO6uzNYf3Zdy7pm2d7VVyY0L/82B791CNb1VBVrPqKEhUXJsJu5rpHsAoRwQoAEDTnnL59vk9/8LVz+sfzfXPbiwsT+mhzvfY+uFlPfWCzttSUhtjK9YtgtQZm1i7poKTjkrol7ZV00jn3YjavJ1gBAHJldtbpzatDunV7Qr0jU/re1SF990KfuntGJUk/17pdv7H3Pu2oLQ+5pesLVwWuTVJSq6T98oJVZ7ahCgCAXEokTB/cmZy3zTmn164M6ovfvay/fu2q/ub1a7pvc5X+xcd2a3/rDm4SnUeMWGVgZvslnXDODa7m9YxYAQDCcm1wXH958oq+/OY1ne8Z1WM7k/qx+xvVuqtWj99Tq4oSxlRWI9ZTgWaWlNQp6fhSI0Vm1iFvdGnQ/9kdxMgSwQoAEHfOOb1w6or++FsX9e6tYTknlRcX6NBP3K9f/Mhubgy9QrGcCjSzo5LqJJ2U1CFvjdNi+3ZK6nPOHUnfZmZN6dvW4Gkz6/f//IRz7nAAxwQAIC/MTM88sUvPPLFLtyfu6PXLg/qjb13Qv33pbX3j3R793oHH1FBJXaygRXLESpLMzEk6kGkEysyaJJ12ztVmeG5A0p7Vjjb5x2iVNOic6/Z/75C01zl3IJvXM2IFAIiimVmnP/3Hi/oPX35bZUUFOvD4Dv0fP9aiTdVcSbicWI5YrcBBSScWea5b3mjXEWluWu+ZLI75fCrEOefOLHjuBUlHzSy5lsAGAECYChKmX/74Hj25p05d/9Ctv/juZf3Fdy9rW7JMj2yv0b7HtuozD28Nu5mxFtcRq9Py1kDdNT1nZsckJZ1ze9dw7v0Lz+u35/EMoesujFgBAOLgYu+oXjh1RZf7x/TdC/3qGZ7UE7trtauuQk+37dCTTfVhNzEy1vuIVauk5xd5LjVitSr+wvljZtacNhWYuq61e7XHBQAganY3VOjQZz4gSZqemdV//eo5ffNsj/6/t27oS2fe02/95Af0K59sDrmV8bJea+Anl98lM3+q72AqVPk6JL241DSgmXWY2SkzO9XT07Pa0wMAEIrCgoR+Y+99+utf+5j+8XOf1mM7avS7f/eOfuP51xXV2a0oit2IVdroUS7XOr1gZofSfq9fbuG6c65LUpfkTQXmsG0AAORUdWmRvvSrH9Wv/PkZ/dVrV1VSVKBf+WSTdtWVy4wyDUuJXbDKB39kasUlG8xsn6R9LS0twTcKAIA8KixI6OgvPq7/84uv6YuvXtYXX72spz6wSb9/4DHVVhSH3bzIit1UYNp03Kqn+3LFOfeSc66jpqYm7KYAALBmBQnTH/5Cq17+vz6pX/rIPXr5nVv6xJGv6fSl/uVfvEHFLlhlIVWJHQAABKC5sVL/7mce1t/960+otCihP/za+bCbFFlxDVYnJC12mUKdpFBqHZjZPjPrGhoaCuP0AADk1ANbq/VzrTv01Xdu6dP/6ev6s29fDLtJkRPXYHVGUtMizzVJOpbHtsxhKhAAsN79evu9+s0fv08TUzP67b95S//m+df199+/rumZ2bCbFglxLRCa8ZY2/hWDA5Jqw6yQToFQAMB6Nzg2pV9//nWdvjig4clpNVSW6NBn7teBx3esyysHsy0QGvURq7pMG/0aU13+jZjTdUo6zG1nAADIrWR5sf70X3xYp397r/6fpx9TbXmRDr34pn7pT07q779/Q5PTM2E3MRSRGrHya0c9Ia+yepO8Rein5FU8vysw+TdHTi1Wb5Z0MtMIV76klVt49uzZs2E1AwCAvJuddfrCP5zXH3/rgnpHpvTzT+zU73720bCbFZhsR6wiFazWC6YCAQAb1fjUjP7Vfz+pkxcG9Na//wkVFUR9ciw762UqEAAAxEhZcYF++tFtmpqZ1U/8l3/Q29duh92kvCJYBYhyCwAASJ9t3aEj+x/VrduT+qd/+IpOXdw4BUWZCswBpgIBAJCuDY7rnz/3HV0fmtDHWxr0H372YW2tKQu7WavCVCAAAAjVtmSZ/sezP6LPPr5D/3C2Rwe+8G39x//1g3U9gkWwAgAAObMtWab/+LOP6M9++UltrSnVn75yUU8f/bZePP2eJu6sv5IMTAXmAFOBAABkNjg2pf1f+LbO3RrRh/fU6YWDHwm7SVlhKjAELF4HAGBpyfJifeXXf1S/sfc+vXqhX3/8rQuamV0/gzwEqwBxr0AAAJaXSJh++eN79OSeOv37//m2Dn/pzbCbFBiCFQAAyLuKkkL991/+sD7aXK+/ff2ahsbuhN2kQBCsAABAKEoKC/Svn7pXUzOzenWdXClIsAIAAKF5bGdSJYUJvXKuN+ymBIJgBQAAQlNaVKBP3teoY6eu6HzPSNjNWTOCVYC4KhAAgJX7tR9rUWFBQr/yZ6djX9uKYBUgrgoEAGDlPrgzqT/45x/S2VsjOvqN7rCbsyYEKwAAELpP3Nuox3bU6D+feFcnY7yQnWAFAAAi4Xd+9hFJ0tfeuRVyS1aPYAUAACLh4e01emhbtf7b18/r5R/cDLs5q0KwAgAAkdH52UfV3Fih//tv3tLY1HTYzVkxghUAAIiMh7fX6N/9k4d1fWhc/+y572pgdCrsJq0IwSpAlFsAAGDtPn5vg77wvz2ut68N6Q++di7s5qwIwSpAlFsAACAYP/7QFn28pUFfefuGnHNhNydrBCsAABBJex/coiv943r3ZnwqshOsAABAJLU/sEkJk/76tathNyVrBCsAABBJm6pL9ZMPb9Vz3+zWH3/rQtjNyQrBCgAARNa//ScPqaK4QH/wtXOanY3+WiuCFQAAiKzGqhL99k8/qP7RKb1zYzjs5iyLYAUAACLtU/d7a62Onb4S+SsECVYAACDSGqtK9NOPbtOfvHJRz32zO+zmLIlgBQAAIu/3DzymzdUl+sa7PWE3ZUkEqwBReR0AgNwoLkzoY80N6u4ZDbspSyJYBYjK6wAA5E7zpkpdH5rQ7Yk7YTdlUQQrAAAQCx/alZQkvXK2N+SWLI5gBQAAYqHtnjrVlhfpV//ijL76zs2wm5MRwQoAAMRCcWFCLxz8iCTpzfeiuZ6ZYAUAAGLj3s1Vqiot1OBYNNdZEawAAECs1FUUq390KuxmZESwAgAAsZIsL9bAWDSDVWHYDYgyM+uUdN7/td8592KY7QEAAFJteZH6RghWsWFmSUkvS3rKOTdoZq2STkuycFsGAADqyot19uZI2M3IiKnAzD4n6Xnn3KAkOefOSNobbpMAAIAkbaou1a3hCc3ORu+GzJEcsfJHjDolHV9q+s3MOiQlJQ36P7sDmq47JOlxM2uS1OScO+GcOxHAcQEAwBptT5bqzoxT78ikNlWXht2ceSIVrMzsqKQ6SScldUg6vsS+nZL6nHNH0reZWVP6tlW0ocn/Y5OkbkndfruOEa4AAAjf9toySdJ7g+MEq6U45w6m/uwHp4z88NPhnKtd8PrDZjZgZl2pabxVSAWrQX8KUGZ2WNIFSbWLvgoAAOTF9mS5JOlK/5had0XrqzlSwWoFDkpabPSoW95o1xFJMrP9kp7J4pjP+9OI/f7vp1JP+AvYk2bWzqgVAADh2tNQoeLChN66dls/88HtYTdnnrgGq3YtHaz2yg9Wflhaybqrbv9nnby1WymDen80CwAAhKS4MKGHtlXruxf65ZyTWXQu2o/rVYGtkvoWea5bUttqD+xPIXbr7hCVVNooFgAACM++R7fpjSuD+t7VaN0zMK7BajnJNb7+sNLKK/jTiSdSa64yMbMOMztlZqd6enrWeHoAALCUj7bUS5LeGxgPuSXzxW4q0C/FIM2fpguUc+5FM6szs0P+pnrn3JJ1rJxzXZK6JKmtrS16hTUAAFhHkmXFkhS5W9vELljlix+UVsTM9kna19LSkoMWAQCAlGR5kSRpcOxOyC2ZL3ZTgWllFNY63Rc459xLzrmOmpqasJsCAMC6VlpUoNKihAYjNmIVu2CVhVQldgAAsI7VlhczYhWQE5KaF3muTiFdvWdm+8ysa2goWlcoAACwHtWUFWmAYBWIM1q8plSTpGN5bMscpgIBAMif2vJiDY0zFRiEo8pQq8q/YrBV0gt5bxEAAMirZDkjVitVl2mjc65bUleG+wl2Sjq8hvsEAgCAmEiWF0VujVWkyi34daOekDfqJEmdZnZAXiX0eYHJv+Fyh/+aQXlrro77t7AJBeUWAADIn2R5sQbHpiJ1W5tIBSvn3JEV7r/iWlO55Jx7SdJLbW1tz4bdFgAA1rtkWZGmZ51Gp2ZUWRKNSBP1qUAAAICMasv96uuj0VnATrACAACxFMXq6wSrAFHHCgCA/KmriN79AglWAaKOFQAA+VNLsAIAAAhGnb/Gqp81VgAAAGtTXVakhLF4fd1ijRUAAPlTkDAly4vVR7Ban1hjBQBAftVXFKtneDLsZswhWAEAgNjaUlOqm7cnwm7GHIIVAACIra01pbo+RLACAABYsy01ZeoZmdSdmdmwmyKJYAUAAGJsU1WJnItOyQWCVYC4KhAAgPwqLy6QJE3cmQm5JR6CVYC4KhAAgPwqK/KC1TjBCgAAYG1KU8FqimAFAACwJqWMWAEAAASjzF9jNXmHqwIBAADWpLTIizKMWAEAAKxRGWus1i/KLQAAkF+pYDUxTbBadyi3AABAfpUWM2IFAAAQiNJCCoQCAAAEoqjAVJAwTXBVIAAAwNqYmcqKCrgqEAAAIAilRQmCFQAAQBBKiwo0weJ1AACAtSsrKqDcAgAAQBBKiwootwAAABAEFq+vU1ReBwAg/0qLCyi3sB5ReR0AgPwrLUxQIBQAACAIZcVMBQIAAASirKiAESsAAIAgcFUgAABAQEqLWLwOAAAQiNKihKZmZjUz68JuCsEKAADEW3GhF2fuzIQ/alUYdgOiyMyOSTou6ZSkwfTnnHPdoTQKAABkVFzwfrAqLSoItS0Eq8xaJe3PsP1FSQfy3BYAALCEIj9YTU2HP2LFVGBmLzrnLP0h6aBzjlAFAEDEFM2NWIW/xiqSI1ZmlpTUKem4c+7FJfbrkJSUN12XlNS91P4rcHTBedrlTQsCAICIKSowSayxuouZHZVUJ+mkpA5565wW27dTUp9z7kj6NjNrSt+2GunrqPyQ1+ScO7GWYwIAgNxILV6fIljN55w7mPqzH5wyMrMmSR3OudoFrz9sZgNm1uWcG1zk5SvVmd4uAAAQLemL18MWqWC1AgclLTaC1C1vtOuIJJnZfknPZHHM5xdOI5pZ61oaCQAAcm9ujdU0a6xWq11LB6u98oOVH5ZWu+7qoKTzq3wtAADIg6K5qcDwb2sT16sCWyX1LfJct6S2gM7TrgV1rAAAQLSkFq9PRWDEKq7BajnJgI7TJC+oLcvMOszslJmd6unpCej0AABgOVFaYxW7YOVfpSflZySpW1J/Njs657qcc23OubbGxsYcNwsAAKRwS5uYcM41h90GAACwtCJGrFYvrYxCUNN9gTGzfWbWNTQ0FHZTAADYMOZuaROByuuxC1ZZSFVizzvn3EvOuY6ampowTg8AwIY0t8aKewWu2glJi03T1YnbzwAAsGEUFfpXBTIVuGpn5F2xl0mTpGN5bMscpgIBAMg/1lit3VFlqFXlXzHYKumFvLdITAUCABCGuXsFMhW4rLpMG/2bJHdluJ9gp6TDAd4nEAAARNz7dazCX7xuzoXfiBQzOyTpCXmjTk3yFqGfkldP6q7AZGYden+xerOkkwvv95dPZrZP0r6WlpZnz549G1YzAADYUGZnnb5zoU/31Fdoe7IsJ+cws9POuWXv7BKpYLVetLW1uVOnWD8PAMB6kW2wivpUIAAAQGwQrAAAAAJCsAoQ5RYAANjYCFYBotwCAAAbG8EKAAAgIAQrAACAgBCsAsQaKwAANjaCVYBYYwUAwMZGsAIAAAgIwQoAACAgBCsAAICAEKwAAAACQrAKEFcFAgCwsRGsAsRVgQAAbGwEKwAAgICYcy7sNqw7ZtYj6VKODt8gqTdHx8bd6O/8oa/zh77OL/o7f3LZ1/c45xqX24lgFTNmdso51xZ2OzYK+jt/6Ov8oa/zi/7Onyj0NVOBAAAAASFYAQAABIRgFT9dYTdgg6G/84e+zh/6Or/o7/wJva9ZYwUAABAQRqwAAAACQrACAAAICMEKAAAgIAQrAACAgBSG3QBkx8w6JCUlDfo/u51zL4bbqvgws0OS6iW1+puOLtZ/K+lrPpfsmNlx59zeRZ6jv9fI75fmtE19zrkji+xHX69S2t8jSX/T+Uz97O9LX2fBzJKSOiUdX+o956o/c9H3XBUYA2bWqQV/UWbahsz8vvq8c27Q/71V0suSTjjnDmTYN6u+5nPJjpntl3TMOWcZnqO/18jMjkk6uaBfjsn7oupK20Zfr4GZHZXU6ZzrTtvWLunwwn800NfL8/uzTtJJecHqwBLhJyf9mbO+d87xiPBDUpOkgUWeG5CUDLuNUX5I6pDUlGF7uyQnaf9q+prPJev+T0o65v1Vc9dz9Pfa+7dT3uhrpj7hv+3g+rk9vT8zfAbt9PWa+tct0b856c9c9j1rrKLvoKQTizzXLS84YHHNLu1fmCnOuRPyhn4Ppm1eSV/zuWSnQ9LzizxHf6+BmTVJOiTvi30e51ytm/+vf/p6bfbK+/sik5N6f4mBRF8HLVf9mbO+J1hFX7u8DzmTbnn/w2Nxh/xpkUxOSUq/WedK+prPZRn+lGsqwGZCf6/NQXnrQRbrl3T09dr0KUOA9e2VdCbtd/o6WLnqz5z1PcEq+lrl/U+dSbfmBwPcbbF/kUje/H5/2u8r6Ws+l+W1O+fOLPE8/b027fK/0M0saWb7/VGsTOjrtemS1GRm5/1/MEia+8dDagQ8hb4OVq76M2d9T7CKv+Tyu2xczrm9bsEC9TSpEZVsraSvN/Tn4l9ps9Z7dtHfS2uV1O9fHNAu77/lpJkdWyJgLYa+XoLzLnx5St4/xk6bWaff703OuYNLv/ou9HWwctWfq+57yi1EmH8ZqrT4VApWyf/il/zh/ZX0NZ/L0vwv9X7/y2ixfejvYNRJGkwbMTljZoflffnvcc4N0tfBcM6dMbM9ki7IW9vWLWnhVcX0dYBy1Z+57ntGrLDhpNVNOZDl+hSszH63QWrwhCXti6F1wTSU/P+muyV9Lu8NW8f8ab/PSdojL1ClRq/2h9owRA7BKsLS/sXPcHCwjsmrPTP35b+SvuZzWZz/JbNsqKK/1yatTxZbw9Ytaf+CfenrVUqFKufcYefcoP93xx5506/HUuGKvg5Wrvoz131PsIq3VLVYZMkv/nbMpRVOzNJK+npDfi7+KEpdQKOA9PfyBjX/4ot0/fLq9GSDvl7ec5KeTd/gB6y98tYSPpflcejrYOWqP9fU96yxir4Tmn+rinR18koGIAv+7SjOLxGqVtLXfC5365D0hJk9vmB7kzRXaVnygu0J0d9r1S3vvS8m/YuBvl6l1LTrYmsGnXMHzexpM0v6+9DXwcpVf+as7wlW0XdG84vPpWuSdHSR55DGH6ofXBiqzKw9feGvsu9rPpcF3OL3TNsvr/TCwqun6O+1eV7zC9yma9L8Lwb6Ore6F0zP0tfByVV/5q7vwy5lz2PphxYpuy9vqNJpg9zyYI19uOTtKFbT13wuK+r//VrBLW3o76z7NfXeM92yySmL26zQ11n39flM/ZzWL8fo6zX174pvabPW/sxl37PGKuKct16ly18blK5T3gJs5uCX4C86Pej/eX/ao8OfGpxbvLiSvuZzWZE6ad6VbJLo77Xy3/NhLfiXtT/l2uXSrhakr9dsr7xF6vPWrfm/z1t/RV+vWsZp7Vz1Zy773vyEhojz6y6lFtQ1y7ubPZe0L8PMBrT0lR+H3d13PM+6r/lcFudPAT4jb8QwKW/ovVvSs+l/adHfa+P3Sfq6ttNukXWE9PXqpZVpkd5fv9a38O+PtP3p6yX4/7B9Qt50XJO8935K3t8RdwWbXPVnLvqeYAUAABAQpgIBAAACQrACAAAICMEKAAAgIAQrAACAgBCsAAAAAkKwAgAACAjBCgAAICAEKwAAgIAQrAAAAAJCsAIAAAgIwQrYYMysycyOmdl5M3P+z2Nm1h5229bCfy/Hw25HStTakw0zS/r/LQz47T+W4/O1ZrgJLhBrBCtgA/FvOHpe0n55Nzvtknfz0f2Sjvs3Tt4QzGx/PsJDzJyWd9PsFyQdkfffSODMrDWbbUAcFYbdAAD5YWZNko7K+7J8PP3u8WaWlNSpHH2RIvrM7JCkJknNzrlc/3fwjD9SdVhSn//nVjM7kP7fJRBHBCtg4ziY+rnwy8v//eDdLwmOPy3W7pyzXJ4nW865FyVFoi0RsVeS8hCq5Jw77If55yS1SjrgnDuc6/MC+cBUILBxpKZaToXaCkRVXb5O5I+ePifppKQTkj63kaahsb4RrICNIzUS8XSorQC80dHDkl6UdN45d0DSXn8UC4g1ghWwcaSuvjpqZoey/RIzsw4zO+4v9D6d6Sou/+qu8/6+STPr9K8sazezo2bm5C2KTl0tl3okFxznkH+OAf+cd12puODKtfP+8ZtW2hl+m136+1nwPtr9tqSunDy0yHFW3J6l3qffd87MjmZ4XefCNmfxPpf8/PznnfwRzbTPpiPbc/iv2++3rzN1jtRVf/77nTuvc+6wP+WYlNTsb7trihqIJeccDx48NshDUockl/Y47m9LLrL/sbT9OtN+P53+Gnlfyk7e4vjz/p/P+9tb/XOktnekHgvOddp//ph/rtT+7Wn7JBccP71NTtLxFfRFqs2dGbadljTgv5+jacffv+AYK25Plu/z+MLzpbdtBe9x2c9P3oL1TJ9P0wrOs39BW0/7x+hMe95Jas3w2qzPw4NHHB6hN4AHDx75ffhf0Kkv7vTHwtCQCmGHFmxvT4WoBcdMHef8Il+gpyW5Rdp0KPWlvmD7+fQg4YcDJ+nYIm0KKljNCwFpwWDheVfUnhW8z6S8YOf84JP6fUCLhOAM7y/rz2+5z2eZ8yTT+9DflgpwTWn9lPVnw4NHnB+hN4AHDx7hPfwv2fQRmfRRkwF5618yvW7hF2fGQLLgNUsFq4wjMWnhIJnWpruO4X+5Bzpilc3xV9qebN/nwrak9Xf7Ct5f1p/fcp/PMue5a3TLD4oEKR4b8sEaK2ADc86dcM4dlHTA33RUmrtqKynviq1MUhXFFxZ1POOcO7OSNqStR2q1+euvUlOLktTmr8dKSlrR8Vdp3pWT7v21P3NXzq20Pdm+z7RznpG3yLtV3ojZEefcYp9HpnOt5vNbMedcl7u7RENT2jmADYU6VgDknHvRzLolNfmBIfWFe36Rl6S+SBcu0l5NKYfUMV6U9Pwi+5zS+6EmKkVMV9qebN/nHOdcl7/oO6n3w1c2Vvv5rVnaQvysQiCw3hCsAKSc0ftftKkv3uZF9l24X8pqrupKhYl+5xXtzMhsrpZnVC7J7/d/ZtuerN5nOv92O0l5/XpM0uNZnmu1n18QUoVG8zGyCEQOU4HABpHFvdjaJQ065wbTvhQXuzFzaupwzV+e/jTboJapr5W2X9tS++XLStuT7ftM8csd7Jc3HfiUvCnErEat8vn5ZbB/4XHNu/E3BUCxIRCsgI3jmF9v6a7pH/8LOynp82mbD8ubGlxY92i/vC/sFzOsrVlKt//6TNNPhyUlF6ndlF5zq8vfb2EdpxXVXArQStuT1ftMu6/jGX8N0xl5N0XuWEFACfrzy8ivkdXh/7lV3mjYwinhg9mO0gGxF/bqeR48eOTnofn1lc7LW1x8TO9f1p+p5lKqLMNp+ZfMp70+0xVsnUuc/1DasY76523P0L7UudJrYqWuPkwvQ3A8rU2Lvocl2rPUVYFHM+y/2NWCK2pPlu9z3u9pr824fYn3mNXn5++74qsC9X7phv3+7wg8nTYAAADtSURBVJ1aUFpBC2pc8eCx3h+MWAEbhPNuG/K4vFGW1BRWu7zRhYPOub0ZXrNX3siH5AWjOnlXpzW7lVfJ7pI3RdTqn7dLaSMbfvsOylu7dMjf54ykx50/suKfc4+8hdHt8kaGzvjbzijPC9tX057l3qe/rqpJ3mey8PWpzyirK+4C/vwyHf+EvJG01MjYUefdTLk7rdr6oGO0ChuIOefCbgMAAMC6wIgVAABAQAhWAAAAASFYAQAABIRgBQAAEBCCFQAAQEAIVgAAAAEhWAEAAASEYAUAABAQghUAAEBACFYAAAABIVgBAAAE5P8HPwQxegHBOroAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(9, 7))\n",
    "plt.semilogy(np.sort(np.abs(x.value))[::-1])\n",
    "plt.ylabel(\"$|x_i|$\", fontsize=24)\n",
    "plt.yticks(fontsize=24)\n",
    "plt.xlabel(r\"Sorted index of $x^*$\", fontsize=24)\n",
    "plt.xticks(fontsize=24)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Проверим DCP правила"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "ename": "DCPError",
     "evalue": "Problem does not follow DCP rules.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mDCPError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-9-53bf8440ab91>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMinimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpower\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mproblem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mProblem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpower\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mproblem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/anaconda3/envs/cvxpy/lib/python3.6/site-packages/cvxpy-1.0.14-py3.6-macosx-10.7-x86_64.egg/cvxpy/problems/problem.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    266\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    267\u001b[0m             \u001b[0msolve_func\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mProblem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_solve\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 268\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0msolve_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    269\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    270\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mclassmethod\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/cvxpy/lib/python3.6/site-packages/cvxpy-1.0.14-py3.6-macosx-10.7-x86_64.egg/cvxpy/problems/problem.py\u001b[0m in \u001b[0;36m_solve\u001b[0;34m(self, solver, ignore_dcp, warm_start, verbose, parallel, gp, **kwargs)\u001b[0m\n\u001b[1;32m    376\u001b[0m                                                               gp=gp)\n\u001b[1;32m    377\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 378\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    379\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cached_chain_key\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mchain_key\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    380\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/cvxpy/lib/python3.6/site-packages/cvxpy-1.0.14-py3.6-macosx-10.7-x86_64.egg/cvxpy/problems/problem.py\u001b[0m in \u001b[0;36m_solve\u001b[0;34m(self, solver, ignore_dcp, warm_start, verbose, parallel, gp, **kwargs)\u001b[0m\n\u001b[1;32m    374\u001b[0m                 self._solving_chain = construct_solving_chain(self,\n\u001b[1;32m    375\u001b[0m                                                               \u001b[0msolver\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msolver\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 376\u001b[0;31m                                                               gp=gp)\n\u001b[0m\u001b[1;32m    377\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    378\u001b[0m                 \u001b[0;32mraise\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/cvxpy/lib/python3.6/site-packages/cvxpy-1.0.14-py3.6-macosx-10.7-x86_64.egg/cvxpy/reductions/solvers/solving_chain.py\u001b[0m in \u001b[0;36mconstruct_solving_chain\u001b[0;34m(problem, solver, gp)\u001b[0m\n\u001b[1;32m     77\u001b[0m             append = (\" However, the problem does follow DGP rules. \"\n\u001b[1;32m     78\u001b[0m                       \"Consider calling this function with `gp=True`.\")\n\u001b[0;32m---> 79\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mDCPError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Problem does not follow DCP rules.\"\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mappend\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     80\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mgp\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mproblem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_dgp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     81\u001b[0m         \u001b[0mappend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mDCPError\u001b[0m: Problem does not follow DCP rules."
     ]
    }
   ],
   "source": [
    "# Non-affine equality constraint\n",
    "y = cvx.Variable(1)\n",
    "obj = cvx.Minimize(cvx.power(y, 3))\n",
    "problem = cvx.Problem(obj, [cvx.power(y - 3, 2) == 0])\n",
    "problem.solve(verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "ename": "DCPError",
     "evalue": "Problem does not follow DCP rules.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mDCPError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-fb37b750a34a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMinimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mproblem\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mProblem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mproblem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/anaconda3/envs/cvxpy/lib/python3.6/site-packages/cvxpy-1.0.14-py3.6-macosx-10.7-x86_64.egg/cvxpy/problems/problem.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    266\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    267\u001b[0m             \u001b[0msolve_func\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mProblem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_solve\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 268\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0msolve_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    269\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    270\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mclassmethod\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/cvxpy/lib/python3.6/site-packages/cvxpy-1.0.14-py3.6-macosx-10.7-x86_64.egg/cvxpy/problems/problem.py\u001b[0m in \u001b[0;36m_solve\u001b[0;34m(self, solver, ignore_dcp, warm_start, verbose, parallel, gp, **kwargs)\u001b[0m\n\u001b[1;32m    376\u001b[0m                                                               gp=gp)\n\u001b[1;32m    377\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 378\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    379\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cached_chain_key\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mchain_key\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    380\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/cvxpy/lib/python3.6/site-packages/cvxpy-1.0.14-py3.6-macosx-10.7-x86_64.egg/cvxpy/problems/problem.py\u001b[0m in \u001b[0;36m_solve\u001b[0;34m(self, solver, ignore_dcp, warm_start, verbose, parallel, gp, **kwargs)\u001b[0m\n\u001b[1;32m    374\u001b[0m                 self._solving_chain = construct_solving_chain(self,\n\u001b[1;32m    375\u001b[0m                                                               \u001b[0msolver\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msolver\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 376\u001b[0;31m                                                               gp=gp)\n\u001b[0m\u001b[1;32m    377\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    378\u001b[0m                 \u001b[0;32mraise\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/cvxpy/lib/python3.6/site-packages/cvxpy-1.0.14-py3.6-macosx-10.7-x86_64.egg/cvxpy/reductions/solvers/solving_chain.py\u001b[0m in \u001b[0;36mconstruct_solving_chain\u001b[0;34m(problem, solver, gp)\u001b[0m\n\u001b[1;32m     77\u001b[0m             append = (\" However, the problem does follow DGP rules. \"\n\u001b[1;32m     78\u001b[0m                       \"Consider calling this function with `gp=True`.\")\n\u001b[0;32m---> 79\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mDCPError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Problem does not follow DCP rules.\"\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mappend\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     80\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mgp\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mproblem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_dgp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     81\u001b[0m         \u001b[0mappend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mDCPError\u001b[0m: Problem does not follow DCP rules."
     ]
    }
   ],
   "source": [
    "# Non-convex objective function \n",
    "y = cvx.Variable(2)\n",
    "obj = cvx.Minimize(y[0]**2 - y[1]**2)\n",
    "problem = cvx.Problem(obj)\n",
    "problem.solve(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Autodiff на примере PyTorch'a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Вычисление производных (градиентов)\n",
    "\n",
    "- [Теорема Баура-Штрассена](https://core.ac.uk/download/pdf/82480031.pdf)\n",
    "- [Лекция с пояснениями](https://courses.cs.washington.edu/courses/cse446/18wi/slides/auto-diff.pdf)\n",
    "- Граф вычислений\n",
    "    - Статический\n",
    "    - Динамический"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Существующие реализации\n",
    "\n",
    "- Статическая (Tensorflow, Theano): сначала строится граф вычислений, после чего выполняется вычисление \n",
    "    - Pros: \n",
    "        - Процесс вычислений в статическом графе можно оптимизировать\n",
    "        - Возможна улучшенная поддержка на уровне железа\n",
    "    - Cons:\n",
    "        - Сложно отлаживать и модифицировать функцию\n",
    "        - Изучение языка построения графа\n",
    "- Динамический (PyTorch, Autograd): граф строится непосредственно в процессе вычислений\n",
    "    - Плюсы:\n",
    "        - Можно писать на Вашем любимом языке, без дополнительных усилий\n",
    "        - Проще отлаживать, поскольку есть доступ к промежуточным вычислениям\n",
    "    - Минусы:\n",
    "        - сложнее оптимизировать процедуру вычисления, чем в статическом графе"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"./plot_autodiff.png\" width=700>\n",
    "\n",
    "Картинка [отсюда](https://ai.googleblog.com/2017/11/tangent-source-to-source-debuggable.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Пример: квадратичная функция $f(x) = \\frac{1}{2} x^{\\top}Ax - b^{\\top}x$\n",
    "\n",
    "- Градиент по $x$: $f'(x) = \\frac{1}{2}(A + A^{\\top})x - b$\n",
    "- Градиент по $A$: $f'(A)$ - ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(-3.7108, grad_fn=<SubBackward0>)\n",
      "-3.7108304500579834\n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "\n",
    "n = 5\n",
    "A = torch.randn((n, n), requires_grad=True)\n",
    "b = torch.randn((n,))\n",
    "x = torch.randn((n,), requires_grad=True)\n",
    "\n",
    "f = 0.5 * x @ A @ x - b @ x\n",
    "f.backward()\n",
    "print(f)\n",
    "print(f.item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([-0.4549,  2.2341, -5.7215,  1.8369, -2.1107])\n",
      "tensor([-0.4549,  2.2341, -5.7215,  1.8369, -2.1107])\n"
     ]
    }
   ],
   "source": [
    "manual_grad_x = 0.5 * (A + A.t()) @ x - b\n",
    "\n",
    "print(manual_grad_x.data)\n",
    "print(x.grad.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[ 1.2039, -0.3909,  1.1023,  1.5957,  0.8493],\n",
      "        [-0.3909,  0.1269, -0.3579, -0.5181, -0.2758],\n",
      "        [ 1.1023, -0.3579,  1.0092,  1.4609,  0.7776],\n",
      "        [ 1.5957, -0.5181,  1.4609,  2.1149,  1.1257],\n",
      "        [ 0.8493, -0.2758,  0.7776,  1.1257,  0.5992]])\n",
      "tensor([[ 1.2039, -0.3909,  1.1023,  1.5957,  0.8493],\n",
      "        [-0.3909,  0.1269, -0.3579, -0.5181, -0.2758],\n",
      "        [ 1.1023, -0.3579,  1.0092,  1.4609,  0.7776],\n",
      "        [ 1.5957, -0.5181,  1.4609,  2.1149,  1.1257],\n",
      "        [ 0.8493, -0.2758,  0.7776,  1.1257,  0.5992]])\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "manual_grad_A = 0.5 * torch.ger(x, x)\n",
    "\n",
    "print(manual_grad_A.data)\n",
    "print(A.grad.data)\n",
    "print(torch.norm(manual_grad_A.data - A.grad.data).item())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Peзюме\n",
    "\n",
    "- Для удобства решения задач выпуклой оптимизации нужно их правильно записывать\n",
    "- Специально написанные солверы могут работать лучше стандартных пакетов, нацеленных на решение произвольной задачи\n",
    "- Градиенты можно вычислять автоматически, используя информацию, сохранённую в процессе вычисления самой функции"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python (pytorch)",
   "language": "python",
   "name": "pytorch"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}