Spring2022/intro_gd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "196b8a50-3d29-45c3-82b9-f4a09b49491d"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Введение в численные методы оптимизации (Ю. Е. Нестеров Введение в выпуклую оптимизацию, гл. 1 $\\S$ 1.1)\n",
    "\n",
    " 1. Обзор материала весеннего семестра\n",
    " 2. Постановка задачи\n",
    " 3. Общая схема решения\n",
    " 4. Сравнение методов оптимизации\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "2a573842-172b-4931-b0dd-9c9d3c47a450"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Обзор материала весеннего семестра\n",
    "\n",
    "Также на [странице курса](https://github.com/amkatrutsa/MIPT-Opt/blob/master/Spring2021/README.md).\n",
    "\n",
    "1. Методы решения задач **безусловной** оптимизации\n",
    "    - Градиентный спуск и способы его ускорения\n",
    "    - Метод Ньютона\n",
    "    - Квазиньютоновские методы\n",
    "    - Метод сопряжённых градиентов  \n",
    "    - Решение задачи наименьших квадратов\n",
    "    - Безградиентные методы\n",
    "    - Стохастические методы\n",
    "2. Методы решения задач **условной** оптимизации\n",
    "    - Методы проекции градиента и условного градиента\n",
    "    - Проксимальные методы\n",
    "    - Методы штрафных и барьерных функций\n",
    "    - Метод модифицированой функции Лагранжа"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Постановка задачи\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "& \\min_{x \\in S} f_0(x)\\\\\n",
    "\\text{s.t. } & f_j(x) = 0, \\; j = 1,\\ldots,m\\\\\n",
    "& g_k(x) \\leq 0, \\; k = 1,\\ldots,p\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "где $S \\subseteq \\mathbb{R}^n$, $f_j: S \\rightarrow \\mathbb{R}, \\; j = 0,\\ldots,m$, $g_k: S \\rightarrow \\mathbb{R}, \\; k=1,\\ldots,p$\n",
    "\n",
    "Все функции как минимум непрерывны. \n",
    "\n",
    "Важный факт</span>: задачи **нелинейной** оптимизации \n",
    "\n",
    "в их самой общей форме являются **численно неразрешимыми**!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Аналитические результаты\n",
    "- Необходимое условие первого порядка: \n",
    "\n",
    "если $x^*$ точка локального минимума дифференцируемой функции $f(x)$, тогда \n",
    "\n",
    "$$\n",
    "f'(x^*) = 0\n",
    "$$\n",
    "\n",
    "- Необходимое условие второго порядка \n",
    "\n",
    "если $x^*$ точка локального минимума дважды дифференцируемой функции $f(x)$, тогда \n",
    "\n",
    "$$\n",
    "f'(x^*) = 0 \\quad \\text{и} \\quad f''(x^*) \\succeq 0\n",
    "$$\n",
    "\n",
    "- Достаточное условие: \n",
    "\n",
    "пусть $f(x)$ дважды дифференцируемая функция, и пусть точка $x^*$ удовлетворяет условиям\n",
    "\n",
    "$$\n",
    "f'(x^*) = 0 \\quad f''(x^*) \\succ 0,\n",
    "$$\n",
    "\n",
    "тогда $x^*$ является точкой строго локального минимума функции $f(x)$.\n",
    "\n",
    "**Замечание**: убедитесь, что Вы понимаете, как доказывать эти\n",
    "\n",
    "результаты!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Особенности численного решения\n",
    "\n",
    "1. Точно решить задачу принципиально невозможно из-за погрешности машинной арифметики\n",
    "2. Необходимо задать критерий обнаружения решения\n",
    "3. Необходимо определить, какую информацию о задаче использовать"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Общая итеративная схема\n",
    "\n",
    "Дано: начальное приближение $x$, требуемая точность $\\varepsilon$.\n",
    "\n",
    "```python\n",
    "def GeneralScheme(x, epsilon):\n",
    "    \n",
    "    while StopCriterion(x) > epsilon:\n",
    "        \n",
    "        OracleResponse = RequestOracle(x)\n",
    "        \n",
    "        UpdateInformation(I, x, OracleResponse)\n",
    "        \n",
    "        x = NextPoint(I, x)\n",
    "        \n",
    "    return x\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Вопросы\n",
    "1. Какие критерии остановки могут быть?\n",
    "2. Что такое оракул и зачем он нужен?\n",
    "3. Что такое информационная модель?\n",
    "4. Как вычисляется новая точка?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Критерии остановки\n",
    "1. Сходимость по аргументу: \n",
    "$$\n",
    "\\| x_k - x^*  \\|_2 < \\varepsilon\n",
    "$$ \n",
    "2. Сходимость по функции: \n",
    "$$\n",
    "\\| f_k - f^* \\|_2 < \\varepsilon\n",
    "$$ \n",
    "3. Выполнение необходимого условия \n",
    "$$\n",
    "\\| f'(x_k) \\|_2 < \\varepsilon\n",
    "$$\n",
    "\n",
    "Но ведь $x^*$ неизвестна!\n",
    "\n",
    "Тогда\n",
    "\n",
    "\\begin{align*}\n",
    "& \\|x_{k+1} - x_k \\| = \\|x_{k+1} - x_k + x^* - x^* \\| \\leq \\\\\n",
    "& \\|x_{k+1} - x^* \\| + \\| x_k - x^* \\| \\leq 2\\varepsilon\n",
    "\\end{align*}\n",
    "\n",
    "Аналогично для сходимости по функции, \n",
    "\n",
    "однако иногда можно оценить $f^*$! \n",
    "\n",
    "**Замечание**: лучше использовать относительные изменения \n",
    "\n",
    "этих величин! \n",
    "\n",
    "Например $\\dfrac{\\|x_{k+1} - x_k \\|_2}{\\| x_k \\|_2}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Что такое оракул?\n",
    "**Определение**: оракулом называют некоторое абстрактное \n",
    "\n",
    "устройство, которое отвечает на последовательные вопросы \n",
    "\n",
    "метода\n",
    "\n",
    "Аналогия из ООП: \n",
    "\n",
    "- оракул - это виртуальный метод базового класса\n",
    "- каждая задача - производный класс\n",
    "- оракул определяется для каждой задачи отдельно согласно общему определению в базовом классе\n",
    "\n",
    "**Концепция чёрного ящика**\n",
    "1. Единственной информацией, получаемой в ходе работы итеративного метода, являются ответы оракула\n",
    "2. Ответы оракула являются *локальными*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Информация о задаче\n",
    "1. Каждый ответ оракула даёт **локальную** информацию о поведении функции в точке\n",
    "2. Агрегируя все полученные ответы оракула, обновляем информацию о **глобальном** виде целевой функции:\n",
    "    - кривизна\n",
    "    - направление убывания\n",
    "    - etc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Вычисление следующей точки\n",
    "\n",
    "$$\n",
    "x_{k+1} = x_{k} + \\alpha_k h_k\n",
    "$$\n",
    "\n",
    "- **Линейный поиск**: фиксируется направление $h_k$ и производится поиск по этому направлению \"оптимального\" значения  $\\alpha_k$\n",
    "- **Метод доверительных областей**: фиксируется допустимый размер *области* по некоторой норме $\\| \\cdot  \\| \\leq \\alpha$ и *модель* целевой функции, которая хорошо её аппроксимирует в выбранной области. \n",
    "    \n",
    "    Далее производится поиск направления $h_k$, минимизирующего модель целевой функции и не выводящего точку $x_k + h_k$ за пределы доверительной области"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Вопросы\n",
    "1. Как выбрать $\\alpha_k$?\n",
    "2. Как выбрать $h_k$?\n",
    "3. Как выбрать модель?\n",
    "4. Как выбрать область?\n",
    "5. Как выбрать размер области? \n",
    "\n",
    "<span style=\"color:red\">\n",
    "    В курсе рассматривается только линейный поиск!</span> \n",
    "    \n",
    "Однако несколько раз копцепция метода доверительных областей \n",
    "\n",
    "будет использована."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Как сравнивать методы оптимизации?\n",
    "Для заданного класса задач сравнивают следующие величины:\n",
    "1. Сложность\n",
    "    - аналитическая: число обращений к оракулу для решения задачи с точностью $\\varepsilon$\n",
    "    - арифметическая: общее число всех вычислений, необходимых для решения задачи с точностью $\\varepsilon$\n",
    "2. Скорость сходимости\n",
    "3. Эксперименты"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Скорости сходимости \n",
    "\n",
    "_1._ Сублинейная\n",
    "\n",
    "$$\n",
    "\\| x_{k+1} - x^* \\|_2 \\leq C k^{\\alpha},\n",
    "$$\n",
    "\n",
    "где $\\alpha < 0$ и $ 0 < C < \\infty$\n",
    "\n",
    "_2._ Линейная (геометрическая прогрессия)\n",
    "\n",
    "$$\n",
    "\\| x_{k+1} - x^* \\|_2 \\leq Cq^k, \n",
    "$$\n",
    "\n",
    "где $q \\in (0, 1)$ и $ 0 < C < \\infty$\n",
    "\n",
    "_3._ Сверхлинейная \n",
    "\n",
    "$$\n",
    "\\| x_{k+1} - x^* \\|_2 \\leq Cq^{k^p}, \n",
    "$$\n",
    "\n",
    "где $q \\in (0, 1)$, $ 0 < C < \\infty$ и $p > 1$\n",
    "\n",
    "_4._ Квадратичная\n",
    "\n",
    "$$\n",
    "\\| x_{k+1} - x^* \\|_2 \\leq C\\| x_k - x^* \\|^2_2, \\qquad \\text{или} \\qquad \\| x_{k+1} - x^* \\|_2 \\leq C q^{2^k}\n",
    "$$\n",
    "\n",
    "где $q \\in (0, 1)$ и $ 0 < C < \\infty$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxEAAAIHCAYAAADgn9FUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3Xt8W9WZN/rfI8mW4rudu8MlcVpCeoa2JIa20M4MkEChJQGa0Os50xkIzBlo55wBEiBMgQ6BJEDfM6W0kwB9+55pmQlJAwklQElgLqU3nFAonzYB4tC0OBc7tiVfZUn7ef/YW7buluwtbVn6feezR9LaS9qPTUv189prLVFVEBERERERZcvldAFERERERDS1MEQQEREREVFOGCKIiIiIiCgnDBFERERERJQThggiIiIiIsoJQwQREREREeWEIYKIiIiIiHLCEEFERERERDlhiCAiIiIiopx4nC6AxjdjxgydP3++02UQERERUQnbv39/l6rOzKYvQ8QUMH/+fLS1tTldBhERERGVMBH5Q7Z9eTsTERERERHlhCGCiIiIiIhywhBBREREREQ5YYggIiIiIqKcMEQQEREREVFOGCKIiIiIiCgnDBFERERERJQThggiIiIiIsoJQwQREREREeWEIYKIiIiIiHLCEEFERERERDlhiCAiIiIiopwwRBARERERUU4YIoiIiIiIKCcMEURERERElBOP0wVQcfn94bew5cd3Ql0uGC4PYD2qyw243BCXByJuQFzmc5cbEA9cLjck+ujyQFweuNzma5erAq7oa+twuyrMR3cFKlwuuF0ueNwueFxueFwuVLjNtgqXG263oMLlts6bbR6XwON2o8Idfe2C2+2CQOAS81FEkl/HtAEYt68LrtFzKftEX4s4/E+OiIiIqHAYIijOqd//Ejd+/5205w2xDlfCY7q2rPpK3Dm12sMCDI33+a7498Sfl+zfk+r9CeciLiDsFoRdQNg9dkRGXwsibjHb3C7AJRC4gJjAMfp/4oIAQDSoAFZbYtCxwkyKAOMSGX2PK7GPjPV1j75ODj4uuADBWFiKqTNdWBq9vvXebPumCmPAWJBzSUItoz9n/Gdn8/lJITFDMMwYJBN+1uh7swqhWX5+9LPiassUWNME2LjfY7rXGf65ERER5YIhguIYIyMZz7vUPGDYeVW188MK8LnZiQgQiYaNhOCR3GaGk4z93TFBJkWfSMr3CIKudNcce0/EDSi/SJY1W4JjhrCWKuCkDXbjvTdDAE4Kjhl+jsSglW1ATRVqswnA6QJqpgCcLrBmG1BTvSfV5433zz5TAE78Q0AuATjx95D2Z0/1z5wBmMhRDBEUJxIKOV1CyXAr4A4DleFseucSePITjsIJYSNzABKMVADDlcBwzGOwUsZeVwLBCmA42hbTPlLB0FJsFIqIRqIviKYMM0xYI7pJAdQFsYKHOfIr8f1H+4wXxuJHi12jIXhsNDg+HFnXjQm15rnYAGeOQ4tr4gF4UsEsQwDOFPwyfV7KgJyHUdhMfwyw7Wcd/Wec+g8b2fbN9NlTHUMExTHCmUciqHR5DPNAVjky3bfM7L99RoNFMC6ISNzr4GjwSA4iwQqJCS/WYwWAEvgXMxFlT6HQaABGhCGYphCx/r8r7eMnZ6/AI5++y5nyxsEQQXFamlsQ8LogqhADo48u/kuZbOYLmQcGY1snF04MIOUIyXBlTDipSAghsWGlwgwy0fbouREPGE6IiMhmav3/SNxjrA5/f0ErygVDBMVZcNHngDc+l/KcGgZgGOZjJAKNGIARgUYiZrv1iEhkrM/oo1p9E94TDgPhEDQ8Yj2Gxh4j1vNIOOZcGBoJAZGw1ScMjUSg4RA0HIYRCcMIh0efazgCtR6NiPkcEbNWNSLmcyNat1WbYQCGmm2jh471iRgQQ822iJo/WwRQwzxgCNTgF04nuBATTuKkCiHZJ2OFYqQCCFWYt3GZhyBYIdYtXIJAlQv+Khf81W70VrnRU+VGT5UHfdNcCIsLBlyIQGBAEIFYf4DSmAOAjD0XaNxrwEh+T9zrxPek6hP/WoR/HSAiKmaC4v0+wRBRQCKyFsAOVW13upaJEJcLcLmK+D/ORcAwgMgINDQEjAxBhwegQwPQ4KD5PDgMDQ5Ah4ehI0PQ4SFocBgIDUODQejIMHRkZOwxFIKGRqAj1mMobLaFrYBkPSIcgYatcBRRaMQwHw2xwk2qx7HnYOjJSCDwhgBvXDiJ/YKf6c0Kj8+AxxeB22eMPlefC5FpHoSnVWLY58XwtGkY9FZhUKowiGkYgA8DUoVB+DCAaeah5vN++NCPKvO1+jCkFVARGIZC1azKUPO5oWpWmvQ61c8yfvCAFVQk9n2jYUQT2mLaJSboxL3HGP09xf9Oo9dJ/IxUoSv5tWTRx3wdXSUiMVile0/yz5v8nhS/ExveM/b7SP69Jr8n/vea/PtIfE/yz5v2PYm/j3HfM14Atv7zlPDPggGYCEU9d4IhokBEZBWATQB2OF0L5ZHLBbh8kAofUNXobOCyAg0iQSBsHZERIDwc89xs1/AwEByCDg/GhJsh6EjQfAwOQ0NBM+iERqDDQzCCQejgEIzhIIzhEfMIhs1jJAIjpNCwC0ZYkg6NlNE+lyoID7kRHnJn6DSCCoyg3tWL6d7YwBGxQocBz7QI3F4Dnmlmu6tCx+6wEjdQWQNU1gLeWsBbY772Rtuir63zlTXQympoZS20sgaG9dyoqAEqq2GIazRwGNb3t2j4iIYTjXltRF9rzHmM9dG4ttg+0c8aex4fcqKfHX89JH5WqmtpfJjC6OfHXDeutpjzmeqG+UIBM7QBmT8L0euP/3maVJv5GjGfP/r7SPzcdO9P/F2n+dzYn2Psd5z8zyf2n40m1mU9R2ItcXUg/rNTXD/6gyX9TAnXMH9/xlg9UKgaMT+bed6AWtcwxn526zD7WzWIMVqLOayso/+H6GeI9ahWMNLodcfORUOTasx7Yf4AY300pg4kB+BxwyySA3DKPrGvU10rRZ/x/ogAWCEu2/eM1SoZ+yR+psa/J2XwTnw9znvS/LEi+jzl7368P0RAE34f4/1zQI4hmCGibIlIC4B1APY7XQuVGSvQoMI3bldJeLRFJASEBoHQkPk4MvZcgwPQ/l4Y/QHzGOiD0T9gPg4OwhgcgjE0BB0atkJKCEbQOkYiMEYMGCMGNIzSCidGNoHDJC6FOyZkmM/98Ph6kkJIXOCIvh9j/7yTrlZRNRZCvLUJQSQ5kCSFFV8d4KsHvPXmfw6JKKP4gJk6RI31TQ5hsf2R8J7koDZ2IrEtMRBnfc2Yc0n9Ez4/8XxijUj4zPjrxP8sVq/RnyX2sxN/nujnIPH8RGpI+nkTfr7Y96WpwVAD0WA7GnRVYSBihXLz3P8xtwHFiiEiz6xbl24EABHZ4nA5RIXjrgDc9eaXyQTRL7CT/nqZEFR0qA9Gnx9Gfy+0z28+H+izjn4YA4MwBgdgDA3DGBqGDg+nHEXRkAFjRGGEASNUnHNc1BCEBz0ID47fNzFwpBrpcMcGjtCg+XsdODnJKgXwWoEim2NaQ/zrylqGECoL5vKgo6+cLIUoawwRRDR1JQQVgflX9fH/jp+DSAg6FIDR1wMj0G0Gk4AVUgb6YPQHEPH3Itzdg0i3H2F/H8K9Awj3DSPSNwJjxNadGSck18AxGjKmReDxGmPPfQY8XuvWKp8Blyd5hCPhykDQbx7+iVQuY6MavnrA15AifKRqi4aQGoYQIqI8YYggIsrEXQGpmQ53zXS45+b+dmNoCOFTpxDu7ETk1CmEu8zn4c7jiHSeRLirE+FTPYj09MIYCtpff47UEIQGPQhlEzjcmnJUo6I6DG9DGN66EFyT+l8ZBYb95jER4ko9EjKtIXP4iA0hRTypkYjISSUXIqwVkNpVNeMEZhFpAHCH9fIUgOkADqvq1jyXSERlxDVtGipPOw2Vp502bl9jcNAKHF0In+oyQ4f1PNzVhUjXKfN8Vxd0aKgA1WemEUFowIPQQJoOoqisDcPXEIa3IQRvQwi+hhA804zCfDdXAxjuNY+JELc1EpIucKRpj96WVVHFEEJEJatkQoQ1gXkTgFWw5iBk6NsAc6LzalU9ENO+VkS2q+rqvBZLRJSCq6oKlVVVqDz99HH7GgMDo4EiesSGjEhX11jgGB4uQPUpqGAkUIGRQAVwdNpos7vSiAsV3oYwvPUhuGy9D80GGgGGesxjIlweK1Q0AtOagKqmmMfGhNcxj1kshkBE5LQpHyJE5AYAqwEcALANZogYz3aY+zUciG1U1c0iclhEbuCIBBEVM1d1NSqrq1F5xhkZ+6kqjIFBRLo6rVBxCuGumFuruqxRD+u5BvN/S1VkxIXBk14MnvSONQrgnV4BbxPMcFE7AG9NHzy+Ao1a5IMRBgZPmUcuKqqsUJEqfKQJI1wJi4gKbMqHCOvL/lZgdDQiI6vPMpjLrqaywzq3NeY9DTCDR7a2M4QQUTEQEbhrquGuqUbl/PkZ+5qBYyB+/kZCyDCDSCciXaegIyP2FapAsCuEYFe0oRpANdwNDfB+4Ez45jfDe8ZM+Jrr4Z3pg4T7x+ZLJB295mM2EzuKUXR1rMCfsn+PuMZGPJJGOTjqQUT2m/IhYgJWAUDiKESM1wCsFZGW6M7SqtoLYHmB6iMicoQZOGrgrqkBFizI2FdVYfT3I9zZhciprrG5HCdPIvjuuxg+dBDhjmOTrinS24vBtl4Mtr0x1ujxwNvSAu/Zi+BbtBjeRWfDd/YieGbMiH9zeAQIBuKDRbpjKMX5sPPzTrKmBkc9iKigyjFELAeQaZZdu/W4DDGjEURENEZE4K6thbu2FmhJHTgifj+GDx1C8OAhDL9tPgbfeWfyt0yFwwi+/TaCb7+NAJ4dbXbPmAHfokXwLloE39mL4F10NrwtCyDVMzJ8WKbrBIHhQPIIR9rXiSHEobkouZjQqIfbnDyeyzyPqibA4x3/s4loyijHENEEoDvD+WjAWJin6xfv1oNERDZy19ej+vzzUX3++aNtGg5j5A9/QPDQIQwfPIThQwcRPHgI4RMnJn29SFcXBrq6MPDqq2ONFRXwLlxohouzz7bCxSJ4mprG/0CPF6iZaR4TERo2g8ZQDzDYDQx1p3jsSW7XyMSuVygaGRv1yGXgo6IaqJoOVM8Aqmdax4yYx5j2qhmApzJvPwIRTV45hoiCfomPWUp2idX0mIi0AdiS4Zaq6ITxGwDgjHEmThIRTRXi8cC7cCG8Cxei7oorRtvDPT0IHnobwUMHMXzwEIKHDiH47ruTn3cRCiF48CCCBw8Cu3aNNntmzhwLFWeZIxeVCxZAPDb+z2KFD6iYA9TOyf49quYtWJlCRqoQMtJvX935EhoA/AOA/2h2/X31KcJGzOuqmLZpjbzFiqjARFWdrsE21qTpwwBuTDexWUR6AHSrasqRhpjP2KqqGZeKLZTW1lZta2tzugwiooLSUAgj771nhYqDGD70NoIHDyLc2ZmX60llJbwf+IAZLhadNTrXwt0wBQaQw8HcRzyGeop/1CNb4rZGORIDR+Jz6zU3EiRKSUT2q2prNn3LdSQi0+1MRERUBKSiAt4PfhDeD34QuPKzo+3h7m4ED46FiuFDhxA8fBgIhSZ1PR0ZwfDvfofh3/0OsXtke+bMgXfRWfBZocJ79tmoPPNMiLuINrbweM0Rj1xGPQzDHPXIasSjeyykpN1d0EEaAQZOmkc2PL7xA0dVzK1WnM9BlKQcQ8QEty4lIqJi4GlqgueCC1B9wQWjbToyguCR98Zuh7LCReRUjqsVpRA+fhzh48cx8J//NdomXq8ZcM5eBN+is82QcfbZcNfVTfp6BeNyWROkG8zZgtkKB7MLG7HtQz3mClLFIjwM+P9oHtnw1ifM35iZInzMiLm1qogCJlGelGOI6EZ2/7pk2CAimiKkshK+RWfBt+gs1K8Yaw93dcXcDmWtENXeDoTDk7qeBoMYfustDL/1VvyoRfPcuFDhXbTIHLUopfv1PV6gbq55ZMswgKDfDBUDnTHHqYTXXebj4CkARXS7ddBvHt2Hx+8rrvFvraqybq2qmQ14a/JfP1EelGOIaAeQaVO66LnXClALERHlkWfGDNR8cgZqPnnhaJsxMoKRw4dHJ3BHV4iK9PRM+nrhjmPo7ziG/ldeGath1izUXflZ1K9cCd9ZZ036GlOSK7oZXiMwPYvFD41IisARDRhdY8+j7cFA/n+GbKkxVls2KqqB2tlAzZyxx5pZ5q1pNbOtxznmMrmcx0FFpBxDxAEAy0SkwdpELlE0ROwtYE1ERFQgrspK+BYvhm/x4tE2VUX4ZCeCbx/C8EEzVAwfOoiRI+8BkclNPg6fPInuJ76P7ie+D+/ixahfuQL1n/kMPDMnuHRsOXC5c1teNzRshYvOhIARfR177iQQsXG39ckKDQDd7eaRiavCDBVJASP6OGssgLgrClM7lbVyXJ0p2me5qiYFBRHZAmBZutWbnMDVmYiInGEEgwi+++5oqAhak7kjfv/4b87E7Ub1hRegfsVK1F5yMVzTptlTMI1PFQj2pQgcCeFj8NTYYzHN5xiXmLdTpQoYiSMelVVOF0tFJpfVmcouRFj9tsNc5jVpCVdrCdg1qrojf5XmhiGCiKh4qCrCJ06MjliYoxeHMPLee+a9/zlyVVej9rLLUL9yJarOay2t+ROlwIiYE8NT3VoV99wa8QhOMmAWUmVtcrConW2NeMSMdExr5K1UZaKcQ8QSAPsBrFPVzRn6NVj9Vsdu+CYimwC0qOrqvBebA4YIIqLiZwwNmaMW0d24D/4eQ795I6dJ3J7muai/cgXqV66AtyXT9D0qWuFgQrBIFz46gf6TgDG5pYkLwu21QkWKgBF9rJltThp3l+Od8qWjrEKEiCwDsA7mXIbYf+O2W8f2VKMSMTtJA8ApAAsB7M80guEUhggioqkp3N2NwJ7n4d+9G8NvvpnTe33nnIP6lStR95kr4GlszFOF5CjDMEc5+k8A/ceBvoTH/pNA33Hz/FTYlVxc5spTSaMaKW6lqvA5XS2lUFYhohwwRBARTX3B9nb4d++Gf/duhDuOZf9Gjwc1n/oU6leuRM1FfwmXlxuflaVgvxkmoqEi9nls2+Dk90YpCF998qhG3Txr+eB5QF2zGTY4slFQDBElhiGCiKh0qGFgsK0N/l270PfCizAGst8B2lVXh7pPfxr1V63EtHPPhfA+dUoUHjFXoBod1bBGNBJHOgZOAsbk9kvJO3EB1bPMQFHXnBwy6pqB2maOatiIIaLEMEQQEZUmY2gIfS+/DP+uXRh49ec5LSdbcfrpqF9hzp+oPOOMPFZJJckwzN3E+46nuJUqYaQjNOh0tZlNa4oJFgkho24eUDsX8E2h3eQdxBBRYhgiiIhKX7izE/7nnoN/924Ef/f7nN477dxzUb9yBeo+/Wm4GxryVCGVpeiSuOluoYod6Ria/IaNeVNZGxMsmpNDRt08bugHhoiSwxBBRFReht9+G4Hdu+Hf/SzCJ09m/T6pqEDNX/4l6q9aiZpPfQpSWZnHKokShINWsIhODD+RPFk8eitVMe694fZmDhl1zeaeGy6305XmDUNEiWGIICIqTxqJYPBXv4J/1y4EXtoLHcz+thJ3QwPqrrgC9VethO+cczh/goqHETEngMeOZvQdB/o6gEAHEHgfCBwzdyEvNuI2J4LHzdNojg8atXMBz9QM8AwRJYYhgoiIjIEB9O3dC/+u3Rj4xS/M20yyVDl/PuqvWon6K69Exbx5eaySyEahYaDvmBks+o5Z4aIj5vGYObpRjKMa1TPTh4zoBPHKaqerTMIQUWIYIoiIKFboxAkEfvIT+J/ZheA77+T03qrzzkP9yhWovewyuGtr81QhUYFEwuZoRqAjYSTDChmB980AEhlxutJkvvrMIaOuGfA1FHSeBkNEiWGIICKiVFQVwYMH4X9mF/zPPYdIV/a3f4jXi9pLLkb9ypWovvBCiIfr8VOJMgzz9qlooBgNGbHH+8W5ClXLXwL/166CXY4hosQwRBAR0Xg0HMbAz38O/67d6Nu7FxoMZv1e9/TpqP/sZ1C3YgV8H/oQ509Q+VEFhv0pQoY1PyP6fLi3sHX92Spg1RMFuxxDRIlhiCAiolxE+vvR9+KL8O/ajcFf/zqn93o/+AHUrVhhzp+YMydPFRJNUSODKeZnHIu/jWog+xXVxnXB14FL/8m+zxsHQ0SJYYggIqKJCr3/PvzP/gT+XbswcuRI9m8UQdXHP4b6lStRt3w5XNXFNwmUqCiFR8wJ36nmZ0RHOPqOAZrF5pKf3gh8/P/Of80WhogSwxBBRESTpaoYfust+J/ZhcBzzyHSm/1tGTJtGmqXLzPnT3z84xB36a6TT1QQRgQY6EwzP6NjbP7GNVuBD60sWFkMESWGIYKIiOykIyPo/9nP4N+1G/0vvwwNhbJ+r2fWLNRd+VnUr1gJ36Kz8lglUZlTNZevLeDmdgwRJYYhgoiI8iXi9yPwwovw79qFoQMHcnqvd/Fi1K9YgfrPfgaemTPzVCERFQpDRIlhiCAiokIYOXoU/t3Pwr97N0JHj2b/RpcL1RdeiPqVK1F7ycVwTZuWvyKJKG8YIkoMQwQRERWSqmLo9d/Av3sXAnuehxEIZP1eV3U1ai+7DPUrVqDq/PMgLlceKyUiOzFElBiGCCIicooxMoL+V/4D/t270f+f/wmEw1m/19M8F/VXrkD9yhXwtrTksUoisgNDRIlhiCAiomIQ7ulBYM8e+HftxvCbb+b0Xt8555jzJ66+Cu6amjxVSESTwRBRYhgiiIio2ATbj5i3O+1+FqGOjqzf554xA80bN6LmkxfmsToimgiGiBLDEEFERMVKDQODbW3w796NvudfgDEwkNX7mq77G8z6+7+HVFbmuUIiyhZDRIlhiCAioqnAGB5G/8svo3fXLgz87FUgknlHXt+HP4x5Dz+EytNPL1CFRJQJQ0SJYYggIqKpJtzVhcBzz5nzJ373u7T9XNXVmHPvvaj/7GcKWB0RpZJLiOC6a0RERGQ7z4wZaPqrv8KCnT/Ggt27UH3BJ1L2MwYG0HHrrehYvx7G4GCBqySiiWKIICIiorzynXUWTn/8ccy85R8AjydlH/+Pd+LI51Zh+Pe/L3B1RDQRDBFERESUd+JyYcaaNZj/w39Fxbx5KfuMHDmC9679PLp/+CPwdmui4sYQQURERAUz7aMfxYKnd6L28k+nPK+hEE7cdx/+dPPXEO7pKXB1RJQthog8E5FNIrJERBqcroWIiKgYuOvqMO9b38Lc+/4J4vOl7NO/bx+OXH0NBl97rcDVEVE2GCLybxmA/QB6REQTju1OF0dEROQEEUHDqlVYsGM7vGedlbJP+Phx/OGvvorOR74DDYcLXCERZcIQkX9tAJYDWJpwbAawxsG6iIiIHOf9wAcw/6ltaPzSF1N3MAx0Pfoo/vDVryJ07FhhiyOitBgi8m+/qu5V1QPRA0AvgNdUtdfp4oiIiJzm8vkw5xvfwLxHvg1XfX3KPkNt+9F+1dXo27evwNURUSoMEXmmqltTNK9S1R0FL4aIiKiI1S1fjpand2La0qUpzxt+P/500804/s1/ghEMFrg6IorFEFFgIrIKQKpgQUREVPYqmptx5v/6AWb83d8BrtRfU3qefBLvXft5BA8fLnB1RBRVciFCRNZaX9TH69dgrZy0yXrPJhG5oQAltvA2JiIiovTE48HMr38NZ/zgf8Ize3bKPsFDh3Bk1Wr07tjBPSWIHFAyIUJEWqzVjjYBaBqnbwPMFZO2qeo6Vd2squsANORzxSQRWQtgb74+n4iIqJRUn38+FjzzNGouuijleR0awrG7/hEdt9yCSF9fgasjKm9TPkSIyA0i8hKAGwFsy/Jt2wHssCY5j1LVzQCW5HFE4sbEaxIREVF6nsZGnPbdRzF7/XpIRUXKPoE9z+PI1ddg6I03ClwdUfmSUhoCFJEWAIdhfllPOe8gps/SVF/oRWQTzInPC2PaGmAGj2xtT7y+iCwBsE9VG3P4HABAa2urtrW15fo2IiKikjL8+9/j/X+4BSNHjqTu4PFg5t9/HdOvuw6SZj4FEaUnIvtVtTWbvp58F1OEVgFAhhGB1wCsFZEWVW23+vbC3OthMpYBaJ/kZxAREZUt3+LFWLBjO45vuB/+nTuTO4TD6Hz4Wxj8xS/RvGkjPDNnFr5IojJRjjF9Ocx9GtKJftFfZvN1zwPQbfNnEhERlRVXdTWa79+A5ocegqu6OmWfgZ//HO1XXY3+//5ZgasjKh/lGCKakPnLfDRgLMzQZyIakDm8xLHmerSJSFtnZ6fNpRAREU1t9Z/9DBY8vRO+P/uzlOcjp07hj2vW4MSDD0JHRgpcHVHpK8cQ0eDQdXuRw+1MqrpVVVtVtXUmh2OJiIiSVJ5xBuY/+SM0Xfc3aft0P/F9vPflr2Dk6NECVkZU+soxRGRc/jWGrWFDVVdby8gSERGRTaSyErNvuw2nP/YY3NOnp+wz/Nvf4sjV18D/k+cKXB1R6SrHEOHUSAQRERHlSc2nPomWZ55G9QUXpDxvDAyg49Zb0XHnehiDgwWujqj0lGOI4G7RREREJcgzcyZOf/wxzLr1FsCTegFK/86dOPK5VRj+/e8LXB1RaSnHENGN7G5pYtggIiKaYsTlwvTrr8f8H/0QFfPmpewzcuQI3rv28+j+1x+ilPbLIiqkcgwR7ch8S1OL9fhaAWohIiKiPJj2kY9gwTNPo+6Ky1Oe11AIJzZswJ9uuhnhnp4CV0c09ZVjiDgAjO5CnUo0ROwtTDlERESUD+7aWjQ//DDm3vdPEJ8vZZ/+l1/GkauuxsCvf13g6oimtnIMEVusx3Rbei8F0G7tUk1ERERTmIigYdUqLPjxDngXLUrZJ3ziBI5+9a/R+e1HoOFwgSskmppSzzpKQ0SusevCqppiv/r8U9V2EdkBYDVSjzZcC2BNYasiIiKifPIuXIj5T23DyU2b0fPkk8kdDANd3/0uBn79K8x78EFUzJ1b+CKJppCcQgSAHQDsmoHktulzYjUkPKazBsB+EVmiqgeijSKyCcBeVd2Rh9qIiIjSDl0fAAAgAElEQVTIQS6vF3O+8Y+ovuAT6Fh/Fwy/P6nPUNt+tF91NZo33IfaZcscqJJoapBcViUQkZeQOUS0wvwCfyDFuQaY8w16Yd4ulO52opyIyDIA66zPbok51W4d21V1a4r3NQC4w3p5CsBCAPtT9XVaa2urtrW1OV0GERFRyQgdO4b3b7sNQ2370/Zp/NIXMWvtWrjSzKcgKjUisj/b7+g5hYhxLloPoAfAqnS3Kll/6b9GVT9oy0XLBEMEERGR/TQcRtf3/gVd3/seYBgp+3jPOgvz/se34F24sMDVERVeLiHCzonVmwAczjTXQVXXAZguIvfbeF0iIiKinInHg5lfuxln/OB/wjN7dso+wbffxpHPrULP9u3cU4Iohp0hohXm7UPjaYM5qZmIiIjIcdXnn48FzzyNmosvTnleh4dx/B+/gY5bbkGkr6/A1REVJztDxBKkXzY1VuLcBSIiIiJHeRobcdqj38Hsu+6CVFSk7BPY8zyOXHU1hn7zmwJXR1R87AwRrwNoyHSrkoisgRkgUk28JiIiInKMiKDpK1/G/Ke2oXLBgpR9Qu+/j/e+8n+i67HHoGnmURCVAztDxAMABMA6Efl3EblYROqs42IR+R6Af4G5utOWjJ9ERERE5BDf4sVY8OMdqP9cmu2xwmF0Pvwt/PH66xHu7CxscURFwrbVmQBARNYC2IjUy8CK9bhJVe9IcZ7S4OpMREREzvD/5Dkcv/tuGAMDKc+7m5rQvGkjaj71qQJXRmQ/p1ZngqpuBvABADsBHIEZHMR6vgPAQgYIIiIimirqP/sZLHh6J3znnJPyfKS7G39ccwNObH4QOjJS4OqInGPrSATlB0ciiIiInKUjIzj5z/+M7ie+n7aP78/+DPMefgiVZ55ZwMqI7OPYSAQRERFRKZLKSsy+7Tac/thjcE+fnrLP8Ftv4cg1n4P/2Z8UuDqiwrN9JEJE5sNcgakhU79Mm9JRPI5EEBERFY9wZyc61t2OgZ//PG2f+quvxpy71sNVXV3AyogmJ5eRCLsnVr8Gc7+Icamq27YLlziGCCIiouKihoHu738fJ/+/fwbC4ZR9KufPx7z/8S34Fi8ucHVEE+PI7Uwi8i8AlmJsIvXr4xxEREREU5K4XJh+/fWY/6MfouK001L2GXnvPbx37efR/f//KzgHlUqNbSMRIvIugAUAlqvqy7Z8KAHgSAQREVExi/T14fjd9yCwZ0/aPjUXXYS592+Ap7GxgJUR5capidUtAA4wQBAREVE5cdfWovnhhzB3w32QadNS9ul/5RUcWXkVBn716wJXR5QfdoaIAwC6bfw8IiIioilBRNDwuc9hwY7t8C5alLJP+ORJHP3qV9H57W9D08yjIJoq7AwRGwEsExEujkxERERlybtwIeY/tQ2NX/5y6g6q6Pru9/CHv/oqQh0dhS2OyEZ2hoifAngMwAERuUVEPioi89MdNl6XiIiIqGi4vF7M+ce7cNqj34G7vj5ln6H9+9F+9TUZl4klKmZ2TqyO5NBdVdVjy4XLACdWExERTU2hY8fw/m23Yahtf8rzMm0aFuz8MbwLFhS4MqJkTk2slhwO7pRNREREJa9i7lyc+YMfYMbNNwOu5K8/OjSEY7ffwTkSNOXY9mVeVV25HHZdl4iIiKiYiceDmTffhDP/1w/gmT076fzQG2/g1BPfd6Ayoonjl3kiIiKiAqg67zwseOZpVJxxRtK5zu98B8MHDzpQFdHEMEQQERERFYinsRHND9wPiMSfCIXQse52GCMjzhRGlCPbJjeLyDW59FfVnXZdm4iIiGiqqFq6FE1/89foTriFKXjoELoe+Q5m3fIPDlVGlD07V0jaASCXpZ7cNl6biIiIaMqY+fWvY+C//hvBd96Jaz/1xBOouegiVC0516HKiLJj5+1Mr49zRFdmOgJgn43XJSIiIppSXF4vmjdtBDwJf881DHTccTuMwUFnCiPKkp2rMy1V1dYMhwvAgwAaAayy67pEREREU5HvQx/CzJv+Lqk99IejOPnQww5URJS9gk6sVtV1AHoArCnkdQtBRNaKSEuac6tEZJPVZ4uIMEQRERERpq9ZA9+HP5zU3vPkk+h/9VUHKiLKjhOrMx0A8AUHrps3VijYlOFcr6quU9XNqnojgMdEZFlBiyQiIqKiIx4Pmjc+APF6k84du3M9IoGAA1URjc+pJV5T/sV+qhGRFhHZAqApQ7fPA+hOaGsDcGPeCiMiIqIpw9vSglm33JLUHj5xAic2bHCgIqLx2RYiRGR+FsfFAJYBaLDruk5S1XZVvVFVt47TNdWoQ2KwICIiojLV+JUvo+pjH0tq9+/ajcBLLzlQEVFmdi7xejjLfgKgbP7boKqrUzS3Is3tT0RERFR+xOVC8/0b0L5iJYyBgbhzx+++B1VLlsAzfbpD1REls/N2Jsni8MPcT+JvbbzulCIimwBsVdW9TtdCRERExaNi3jzMvvPOpPZIdzeO3X03VHPZjosov2wbibCWcM0rEVkLoF1Vd4zTrwHAHdbLUwCmAzicxW1HeWNNsF4Oa5K1U3UQERFR8aq/5mr07d2L/ldeiWvv37sP/l270HDVVQ5VRhTPqYnVObEmMG+HeQtQpknM0QCxH8C2mBWR1gFosD7DEaq6w1qZaYuI7E+3HCwRERGVLxHB3G/eC3dD8vTRE/dtQKijw4GqiJLlNUREJ1RP4v03iMhLMFcy2pbl27YD2KGqB2IbVXUzgCUicsNE67GDqrYD2AuzTiIiIqI4npkzMefee5Pajf5+dKxfDzUMB6oiimfnxGoAgIjUAXgMMbtSiwhgTqZerap92X6WdfvRVuszxv3LvdVnGYB0twvtsM5tjXlPA3L7Qr/dhtuiXgOwVkQaVLV3kp9FREREJabuskvRd+WVCDz7bFz74C9+iZ4n/w1NX/myQ5URmWwNESJyCYCnADRaTb0YW871UgA9IrJcVV9J9X4brAKAxFGIGNEv7y3WiACsL/HL81GMiCwBsA/A0uj1LNHg0BTznIiIiGjUnLvWY/DXv0b4xIm49pMPPYTqCy+Ad8EChyojsnefiHqYf9FvhPnX/kZVbbImXC+EOTrhArDdGq3Ih+XI/KU8+kW+kLtFtyF5T4glMCdYt6foT0RERAR3fT3m3ndfUrsOD+PY7XdAw2EHqiIy2Tkn4g4A9QBWqeqDquqPnlDVI9ak4tth/vV9o43XjdWEzJu4RQPGwjxdP24WlDUikmpPjBsBrMn0QdZ8kDYRaevs7LSxRCIiIpoqaj71STR84fNJ7UNvvIFTjz/hQEVEJjtDxDKYy6/uTNfBmtzcC+A8G68bq6A7YYtIg4hssiZ/A8BjIrLFuo0JwOjPfK3Vb621QtSN4y1Tq6pbVbVVVVtnzpyZx5+CiIiIitns225DxRlnJLV3Pvoohg8edKAiInvnRCxBdjtRt1t982G8kYgoW8KGNZ9i3D0fnNyfgoiIiKY2V3U1mjc+gD98+StA7IZzoRA61q7D/B3b4aqsdK5AKkt2jkQcAJDN3gctVt98KOhIBBEREVEhVC1ZgunX/U1Se/Dtt9H1yHccqIjKnZ0hYh+AFhG5Ol0HEbkN5hf9NhuvG4srHREREVFJmvH1r8P7wQ8mtZ964gkMHnjdgYqonNkZIu4HEACwQ0RuiV2BSUTqROR7MCdU9yCLW4AmqBvj7GhtYdggIiKiKcVVWYnmTRsBT8Ld6IaBjttvhzE46ExhVJZsCxHWakyXwAwSm2HuCXFKRCIwg8ONAATmhnMBu66boB2Zb2mK3m71Wp6uT0RERJQ3vg99CDNvvimpPXT0KE4+9JADFVG5snMkIrqk6XwAO2EGhkbrUQDsBbBQVV+285oJDgCju1CnEg0Re/NYAxEREVHeTL/+evg+8uGk9p4n/w39r77qQEVUjmwNEYA5IqGqq2M2mVuoqi5VvVRVj9h9vQRbrMfWNOeXwlyGlrczERER0ZQkHg+aH9gI8fmSzh27cz0igXzd8EE0xvYQEcvaZC7fwSH2eu0AdgBYnabLtcjffAwiIiKigvC2LMCsf/iHpPbwiRM4sWGDAxVRuclriLBZQ8JjOmsALIvd8A0ARGQTgL3jbfJGRERENBU0fuXLqPrYx5La/bt2I/DTnzpQEZUT0dhNS+z6UHNlplaMzUFoU9XfTOBzlsEcOWhB/B4U7daxPdVGbtaciDusl6dg3la1f6pu+tba2qptbflaFZeIiIimqtD776N9xUoYAwNx7e7GRrQ8uxueGTMcqoymIhHZr6rppgXE97UzRIjIfACbAKxKcVoBrFPVh227YJlgiCAiIqJ0enc+jWN33pnUXnPJJTjtO49ARByoiqaiXEKEbbczicgCAPthBggB4Ie5Ad0+67kLwGYRecGuaxIRERGVu/qrr0LNxRcntffv2wf/M7scqIjKgZ1zIjbBXNL1CIDlqtpkrch0qao2wby96T0Ay0XkFhuvS0RERFS2RARzv3kv3I2NSedObNiAUEeHA1VRqbMzRCyDecvSKlXdl3jS2kNiKcxRis/beF0iIiKisuaZMQNz7rknqd3o70fH+vVQwyh8UVTS7AwR3QAOZJpAbe3P8BjMMEFERERENqm77FLUXXllUvvgL36Jnh896UBFVMrsDBEHsuz3LsyVlYiIiIjIRnPuWg/P7NlJ7ScffhjBIwXbuovKgJ0h4ikAS0Skdpx+yzG2szQRERER2cRdX4+5KTab0+FhdNx+OzQcdqAqKkU5hQgRmZ/uANAGc1L1PhH5aJrjewAuAfCS7T8JEREREaHmkxei4YtfSGoffuNNnHr8CQcqolKU0z4RIhKx45oADFX12PBZZYH7RBAREVEujMFBtF91NUJHj8afqKjAgqe2wbd4sTOFUVHL5z4RAZh7Pkzm6LU+h4iIiIjywFVVheaNDwCJG82FQuhYuw7GyIgzhVHJyClEqGqjtf/DpI98/UBEREREBFQtWYLp11+X1B585x10PfKIAxVRKbFzYjURERERFZEZX/savGedldR+6onvY/DA6w5URKWCIYKIiIioRLkqK9G8aSNQURF/wjDQcfvtMAYHnSmMpjyGCCIiIqIS5lu8GDNv+ruk9tDRozj50EMOVESlgCGCiIiIqMRNv/56+D7y4aT2nif/Df0/e9WBimiqY4ggIiIiKnHi8aD5gY0Qny/p3LH16xHx+x2oiqYyhggiIiKiMuBtWYBZt9yS1B4+cQLHU+xyTZQJQwQRERFRmWj88pdQ9fGPJ7UHdj+LwE9/6kBFNFUxRBARERGVCXG50Hz/BrhqapLOHb/7HoS7uhyoiqYihggiIiKiMlLR3IzZd96Z1B7p6cGxu++BqjpQFU01toUIEblVRK636/OIiIiIKD/qr74KNRdfnNTev28f/E8/40BFNNXYORJxJ4B/sfHziIiIiCgPRARzv3kv3I2NSedO3H8/Qh0dDlRFU4mdIWI7ABGRi2z8TCIiIiLKA8+MGZhz7z1J7UZ/PzruXA81jMIXRVOGnSFiEwABsEVEam38XCIiIiLKg7pLL0XdiiuT2gd/+Uv0/OhJByqiqcK2EKGq7QBWA/gAgPdE5BYR+aiIzE912HVdIiIiIpq4OXfdBc/s2UntJx9+GMH2Iw5URFOBnROrIwCeAqAAGgFsBrAfwOEUx7t2XZeIiIiIJs5dV4e5KTab0+FhdNx+OzQcdqAqKnZ23s4UAOC3jt6Y56mOgI3XLQoislZEWtKcWyYiW6w+m0SkodD1EREREaVT88kL0filLya1D7/5Jk49/rgDFVGx89j1QaqaPL2/TIjIKphzQnakOLcMwDpVXR7Tth/A0sJVSERERJTZrFtvRf+rryL0h6Nx7Z2Pfhc1f/EX8C1e7FBlVIy42dwkiEiLiGwB0JSh2xaYASPWXhG5IX+VEREREeXGVVWF5gc2Aq6Er4ehEDrWroMxMuJMYVSUGCImQVXbVfVGVd2a6rx121ILgLaEU68BuDHf9RERERHlomrJuZh+3d8ktQffeQddjzziQEVUrGwPEdaKTN8TkddE5OKEcxtF5Lt2X7OItQCAqvYmtLdHzxEREREVkxlf+xq8Z52V1H7q8ScweOCAAxVRMbI1RIjINpgrMt0IYAmAxAnEAuDGMtqQrh0YHZGI1YTk3w0RERGR41yVlWjetBGoqIg/oYqO2++AMTDgTGFUVOxc4nUNzH0i9gFYDjMwxFHVdVb7hG7lsVY3WpVFvwZrFaRNMSsiFXwOgjUCcQBAa8KpJdE6C10TERER0Xh8ixdj5k03JbWHjh7FiYcecqAiKjZ2jkSshrlHxCpYf4FP4wByXJnImsC8HeYE5UyTmKNfzPcD2Kaq61R1sxVeGqzPKLTViAlN1jKwvUDK25yIiIiIisL066/DtI98JKm999/+Hf0/e9WBiqiY2BkiWgEcUNXx9oDoRpbzAUTkBhF5CeaX8G1Z1rEdwA5VjbtpT1U3A1hS6BEJayfvNSKyyhpFaYD5O2CAICIioqIlHg/mbnwA4vMlnTu2fj0ifr8DVVGxsG2fCJhfjLPRCnM0YlzWqkdbgdG/4Gdk9VkGYF2aLjusc1tj3tMAM3hka3u61ZjSsUYcRveQEJFWZB6tISIiInKcd8ECzLrlFpxI2NE6fOIEjt+3AfMe3OxQZeQ0O0PE6wCuEZEz03UQkY0A6gHstfG6sVYBQOIoRIzXAKwVkRZrhCD6BX95mv75shTZj6wQEREROabxy19C3759GPzlL+PaA88+i9ply1B32aUOVUZOsvN2puik6ZcAnGu1KTC67OuLAG6DeRvPAzZeN9ZyZL5NKPrX/2V5un4SEXkpxWTwa63bq4iIiIiKmrhcaL5/A1w1NUnnjt99N8JdXQ5URU6zLURYf9m/FsAHYN4epAAeF5EIzInOywH4ASzLYt7ERDUh821V0YCxME/XT7XaUnSFJgCAtcP1mvE+yJoP0iYibZ2dnTaWSERERJSbiuZmzF6/Pqk90tuLY9+4G6rqQFXkJFv3iVDVHTC/yO+EGRgaYY5O+GHOQ1igqq/bec0EBV0yNWYp2ZespsdEZIuILInptg7AquiSszDnVOxI/rR4qrpVVVtVtXXmzJn5KJ+IiIgoa/VXrUTNJZcktfe//DL8Tz/jQEXkJDvnRAAYnWOwGgBEpN5qK9T0/fFGIqJsCRvWz5puEne0TzsA3rpEREREU5qIYO4370X7gQOI9PTEnTtx//2o/tj5qJg3z6HqqNBsHYlIpKr+AgYIgLtAExEREeWNZ/p0zLn3nqR2o78fHXeuhxpG4YsiR+Q1RDiAey8QERER5VHdpZeifuWKpPbBX/0KPT/8kQMVkRPyEiJEpE5ELhaR663jo/m4TgrdGGdHawvDBhEREdEEzV6/Hp45c5LaTz78MILtRxyoiArN1hAhIvNFZBuAHphLvW6xjv0iEhGRW+y8XgrtyHxLU3TDutfyXAcRERFRyXLX1WHuhvuS2jUYRMftt0PDYQeqokKyLUSIyAKYS7muwtiKTPusw2+1bRaRF+y6ZgoHrFrSBYloiMjXZndEREREZaHmwgvR+KUvJrUPv/kmTj3+uAMVUSHZORKxCWNLuq5S1SZVvdQ6mgD8rXVueR5HJLZYj61pzi8F0G6tqkREREREkzDr1ltRceYZSe2d33kUw7/7nQMVUaHYGSJWwdxgbrmq7kw8qapbYW5GJwC+YON1Y6/RDmAHrCVmU7gW4yzJSkRERETZcVVVofmBjYAr4StlOIyOdbfDGBlxpjDKOztDRDvMv/LvS9fB2mStHcCSdH0yaEh4TGcNgGUJG77B2uhtbzYbvRERERFRdqqWnIvp112X1B585x10ffvbDlREhWBniDiA7FY9arf6jktElonISyJyGOZ8CwDYJCKHrfYbEt9j3aq0FMDnrV2i14rIFgCHVTXdCAURERERTdCMr90M76JFSe2nnvg+Bg9k9bWPphhRVXs+SGQZgBcBNKhqX4Z+3QDWqipn3GSptbVV29ranC6DiIiIKK3hgwdxZPW1QCgU115xxhloeXonXNXVDlVG2RKR/aqabm5xHDtHIt4F8DqA7SLy0TTHiwBOAWhLdd7GWoiIiIiogHxnn42ZN92U1B46ehQnHnzQgYoonzw2ftbhmOfLx+m7P0Wbwt56iIiIiKiApl9/HfpfeQVDb7wR197779tQe8ky1Hzqkw5VRnazcyQiAHM/iIkeARtrISIiIqICE48Hczc+APH5ks4dW78eEb/fgaooH2wLEaraaO0NMeHDrlqIiIiIyBneBQsw69Zbk9rDJ0/i+H0bHKiI8sHOkQgiIiIiIjR+6Yuo+sTHk9oDzz6LwAsvOlAR2S0vIUJE6kTkYhG53jo4aZqIiIioTIjLheb774erpibp3PF77kG4s9OBqshOtoYIEZkvItsA9AB4CcAW69gvIhERucXO6xERERFRcaqYOxez169Pao/09uLYN+6GXdsMkDNsCxEisgDmqkurAAjMydL7rMNvtW0WkRfsuiYRERERFa/6q1aiZtklSe39r7wC/86nHaiI7GLnSMQmAI0ww8Iqa7L0pdbRBOBvrXPLOSJBREREVPpEBHPvvRfupuT1c07cfz9C77/vQFVkBztDxCqYez0sV9WdiSdVdSuAa2EGiS/YeF0iIiIiKlKe6dMx5957ktqNgQF03LkeahiFL4omzc4Q0Q6gXVX3peugqjusfktsvC4RERERFbG65ctRv3JFUvvgr36Fnh/+yIGKaLLsDBEHAPRm0a/d6ktEREREZWL2+vXwzJmT1H7y4YcRbG93oCKaDDtDxFYAS0Skdpx+rTBXbCIiIiKiMuGuq8PcDfcltWswiI7b74CGww5URRNlZ4h4F8DrALaLyEfTHC8COAWgLdV5G2shIiIioiJTc+GFaPzSl5Lah998E6cee8yBimiixK41ekUkMsmPUFX12FJMiWltbdW2tjanyyAiIiKaNGNwEO1XX43QH47Gn/B4sOCpbfB96EPOFEYQkf2q2ppNXztHIgIw94OY6BGwsRYiIiIiKkKuqio0b9wIuBK+hobDOLH5QWeKopzZ9pd/VW2067OIiIiIqHRVnXsupl93XdItTIO//CVCx4+jIsUEbCoudo5EEBERERFlZcbXbkbFGWcktQdeeMGBaihXDBFEREREVHCuykrUX3llUntgz/MOVEO5su12JhG5Jpf+qXa1JiIiIqLyUXfF5eh69NG4tuE338TIH/+IytNPd6gqyoadqyHtAJDLUk9uG69NRERERFOMd+FCeBctQvDQobj2wPMvYMYNaxyqirJhZ4h4HZlDxBLrsd06iIiIiKjM1V1+OToTQ8SePQwRRc62ORGqulRVWzMcLgAPAmgEsMqu6xIRERHR1FV3xeVJbcGDBxFs59+ci1lBJ1ar6joAPQAYLYmIiIgIlWecAd855yS1c4J1cXNidaYDAL7gwHWJiIiIqAjVXZ48GhHYswequUy3pUJyaonXFoeuS0RERERFpu7yTye1jbS3I/j22w5UQ9mwc4nX+Vl0awGwDEC9Xdd1moisgvlzLbQe16nqgZjzDQBuALBDVdtFZAmAzwPYoqq82Y+IiIjKXsXcuZi2ZAmGDhyIaw88twe+RYscqooysXN1psNZ9hMAL9l4XcdYAaJdVXdYr1sAHBaRpTFBogXAJgCbRAQwV6a6kQGCiIiIaEzdFVckh4jnn8fM//f/gfUdioqInbczSRaHH+Z+En9r43Wd1BI76mAFg80AHkvotxTmqlQLVXWhqu4tYI1ERERERa/usksBV/xX09Af/4jht95yqCLKxM4lXl1ZHE2qeq2qHrHruk6xblNabj3Geg1je2KMUtVejj4QERERpeaZORNV55+f1B54bo8D1dB4nJpYXSpawUniRERERLZItWdE4IUXoIbhQDWUCUPEBFkjC42xtzNZzoO5jG2sFhFZZR03iMiyApVJRERENGXULl8OeOKn7IaPH8fQ6687VBGlY+fE6rwTkbWImcicoV8DgDusl6cATAdwWFW35rlEwFyJKXYzvV4AiK1ZRPaLyJoUAYSIiIiobHkaG1F9wScw8F//Hdce2PM8qpYudagqSmVKjESISIuIbIe5ylHTOH0bAOwHsE1V16nqZmun7AbrM/JZ5yYAW2MDg6qmCj1bYP4sRERERBSj7oorktoCL7wADYcdqIbSKeoQYd368xKAGwFsy/Jt22HuyRD3V35V3QxgiYjcYHOZAABr/4clVmAZTzvM+RREREREFKP2kksgFRVxbZFTpzD42msOVUSpFPXtTNbtR1uB0T0YMrL6LAOQ7ov8Duvc1pj3NMAMHtnannhbVPT2KVVdnqKmtVaAidUNIHFVJyIiIqKy566tRfVf/Dn69+6Law/seR7Vn/iEQ1VRoqIOEROwCgAyzDV4DcBaEWmJLreqqr0Akr7852gTYuZBiEiDqvZaoWaTiOxIWN61CdZcCSIiIiKKV3f55Ukhou+nP8Wcb/xj0igFOaOob2eagOXI/OU8+kXettWRrMne66wwEnUtMLr5XKrdqZcDeMquGoiIiIhKSe1FF0GmTYtri/j9GPjFLxyqiBKVWohognmrUDrRL/oL7biYiKyCGUxaRGSJdSyDuUN1bL+WhOerkP6Wq2i/G0SkTUTaOjs77SiXiIiIaEpwVVWh5i//IqmdG88Vj1K7nalg8wzGmUsRuzrTVisQNMBcarYBwNKEkYsksfNBWltb1Z6qiYiIiKaGuiuuQN/zL8S19e3bByMYhMvrdagqirI9RIjIR2GuptQK8zafl2PObQRQp6p/Z/d1LeONRERNOmxYIUCy7FuI/SmIiIiISkbNn/85XNXVMAYGRtuM/n4M/Pd/o3YZ9+11mq23M4nINph7NNwIYAmSv6wLgBtF5CI7rxuDKx4RERERlQCX14vaZZcktQf28NPdYnAAACAASURBVJamYmBbiBCRNQBWA9gHc+Jw0l/prT0UBGbIyAeueERERERUIlJtPNf3yn/AGBx0oBqKZedIxGoACnPScOJqRLEOIGHisY26Mc6O1haGDSIiIqIiV/2JT8BVXx/XpkND6P+P/3CmIBplZ4hoBXBAVQPj9OsGMO7GcRPUjsy3NEWvyy0PiYiIiIqcVFai7tLk7bwCzz/vQDUUy84Qkc2EZsAKGzZeN9YBYHTlpFSiIWJvnq5PRERERDZKdUtT/3/+FyJ9fQ5UQ1F2hojXASwRkTPTdbBWZ6pH/r7Eb7EeW9OcXwqgfbzlVYmIiIioOFSddx7c06fHtenICPr27UvzDioEO0NEdNL0SwDOtdoUMJd9FZEXAdwGcz7CAzZed5S1M/QOmPMzUrkW42zyRkRERETFQzwe1F12WVI7b2lylm0hwvoCfy2AD8DchE0BPC4iEZjLvi4H4AewLIt5E6k0JDymswbAMhFZEtsoIpsA7FXVHanfRkRERETFqO6Ky5PaBl79OcI9PQ5UQ4DN+0RYX9CbAOyEGRgaYY5O+GHuvrxAVV/P9vNEZJmIvCQih2EGEQDYJCKHrfYbUtTQC/O2pc+LyCYRWSsiWwAcVtV0IxREREREVKSmLVkCz+zZ8Y3hMPr2cpqrU2zfsdr6Er8aAESk3mrzT/Cz9mIC8yesGnjbEhEREVEJEJcLdZdfju4f/CCuPbBnDxpX82/ETrBzs7mPisj82DZV9ccGCBFZICIPiEidXdclIiIiotKX6pamwV/9GuGuLgeqITtvZ9oPYO04fRpgjhCssfG6RERERFTifOecg4rTTotvNAwEXnzRmYLKnJ0hQqwjLWs+RDuAL9h4XSIiIiIqcSKCusuTRyMCe7hKkxNsnVidg3ztWE1EREREJaruM8kbzw3t34/Q8eMOVFPeJjWxWkS2JTQtS9EWdx7mLU3c7I2IiIiIcuJdtAiVLS0YaW+Paw88/wKm//VXnSmqTE12dabE6fALrWM8XDmJiIiIiHISvaWp69FH49oDe/YwRBSYXSFCADwFc7fqTJu5dcPc8G1CS75SblQVfX19CAQCGBwcRCQScbokIiKyuN1uVFVVoa6uDrW1tRDJOK2QiCx1VySHiOHf/hYjf/wjKk8/3aGqys+kQoSq/jj63PqXX7uqPjbZomjyVBUnT57EwMAAmpqaMGfOHLjdbv6PFBFREVBVRCIR9Pf3o6urC0NDQ5g1axb/HU2UBe/ChfAuWoTgoUNx7YE9z2PGjUn7EFOe2DmxejOA7TZ+Hk1CX18fBgYGcOaZZ6KhoQEej4f/40REVCREBB6PBw0NDTjzzDMxMDCAvr4+p8simjLqrkieYB14nqs0FZJtIUJVb1fVl+36PJqcQCCApqYmuN1up0shIqIM3G43mpqaEAgEnC6FaMpItfFc8OBBBA8fdqCa8uTUEq+UZ4ODg6ipqXG6DCIiykJNTQ0GBwedLoNoyqg8/XT4zjknqZ17RhTOZCdWJxGROgCtGGcvCFV93O5r05hIJMJRCCKiKcLtdnPxC6Ic1V1xBYZ/+9u4tsDzz2PGzTfxFu4CsDVEiMjnYK7SlA2GiDzjf4GIiKYG/vuaKHd1n74MJzdtimsbaW9H8NAh+M4+26GqyodttzOJyLkwJ1YLgNcB7Et4Hn0NZF4GloiIiIgoo4q5czFt6dKkdt7SVBh2zonYBEABLFHVVlW91Gq/X1UvtV5fCjNIbLHxukRERERUhlJNsA7s2QNVdaCa8mJniGgFsENVf5PQPjo3QlX3AjgAgIv4EhEREdGk1F12GeCK/zob+tOfMPzWWw5VVD7yvTpTO4CFKdqW5/m6RERERFTiPDNmoOpj5ye1B57b40A15cXOENEN4H+3d78xbqT3neC/jzQjeTRST3V7xk5w9q5VdLyePW/skN1Z/4FjOyqOLc29uOwVpRebwwF3aVaA+wPcJcfazgLB3uFutcUDFotLDhtSCwT3JgeJ3Cz2xUpjk9qMZwEb625WADvribNmyXY2f+yxyJqWNP+l516wHqpIVrFJdpHFZn8/QKFHD4tVD4s13c+vnuf3PNmhstsALg+VZQFoCZ6XaGHK5TJyuRyEEFhfX0c+n4dt2/B9P5Hj+74PIQQymeHYO9n3EB2G67rwPC/tahARAQDWLkYMabp1C/LRoxRqc3wkGUT8SwB6MMWrUgOwLoT4x0KITwkh/jl6w5v414eOFM/zkMlkYNs2XNeFrvdG6TWbTZTLZayvry9Vo8p1XQghOOMLDfA8D5ZloVwuo1wuw7KsmQLgSqWCTCaDTCaDfD6PfD6PXC6HTCaDep3zZhDRYp3L54EnBiccfe/HP8abf/zHKdXoeEgyiPg99JKm+/kOQQ7EbQA2gBYAK3jJGXk30RIrFArwPA/FYhHdbhftdhvdbhfdbhdOML3cxsZGyrUkiue6LgqFAiqVCkqlEkqlEmzbRi6XmymQMAwDnU4HzWYTnudhc3MTjUYDpmnOofZERPGeWF/H05/77Eg5hzTNV2LrREgp7yA6KDHRWxPiF9DrgahIKf9lUuclmrd6vQ7XdWEYBiqVwYnFNE1DqVSCaZrQtOUZpafrOhqNRtrVoCWiAogwXddhWRa2t7dRq9UmPpamaby/iGiprF28iAdff2WgbP+rX8UHf2sH4onE11YmzGHF6mFSytcBFOZ9HqJ52d3dBdBrhMVRw5uWhaZpMAwj7WrQklC9BVH3hGma/byeZQqEiYimce7CBfz1qVOQ77zTL3t49y7e2N3F05/5TIo1W11JLjb3m0KIX0vqeETLIqmkaaK01Gq12EBXld+4cWORVSIiStTJc+fw9C99fqR8/yaHNM1LkjkRv4VeXgTRSsnnezMSO47DgIKOpL29vbG9ZZqmodVqLbBGRETJe+bSpZGye19rDPROUHKSDCJqAIQQ4ksJHpModaZpIpvNwvM85HI51Ov1iYIJNfXq+vr6yGvVahVCCFiWFfHOx1Ti6/r6OnK5HJrN5kR1jpv2VZXncjkAgGVZWF9f709X67pu7DGbzSby+TzW19eRyWQiZ/ZRybuZTKZ/ftu2x9ah2Wz2Z/aZh2q1ikKhANu2Ydv2gdd8FXmeN3ao0sbGxkyzi1WrVZTLZVSrVViWtVQzlBHR8XP2i1+EeOqpgbKHr7+OB9/8Zko1Wm1J5kQ4ALYBVIQQOSnlvQSPvbSEECaALQB30VtYryGljJzjUAhRQm9V79T/0n7kH/ybtKswdz/4Jy8mdqxarYZ8Pg/P8/q5EYZhoFAooFhMfgF2z/P6wcfm5iY8z4Prusjn82g0GofOd/B9H7lcrj+rzt7eXr8x3+12Rxqc5XK5HwwYhgHP81CtVnHjxg3cuXMHmqb1jwkA2WwWpmnCdV2Uy2W4rhuZiGvbNsrlcv89SfJ9HxcuXMDm5uZA0rAKzMJP3l3XhaZpS5fbkpRJgt5OpzPV8arV6sC9r77/Wq2W+HdJRDSJE2fO4NyXvoj9m7cGyvdv3sLZL3whpVqtrsR6IoKGcQHARwH8QAjxG8HaEB+J2pI6b5qCAMKXUtpSyrKU0gJwTQgx0sIL9uXUtkeUrutot9solUoDa0RYloVMJjP2Cf6sdnZ20O120Wg00O12+1NnJvEk3fM8bGxsDBw/bmy867qwbRuGYfT3b7fbqNVq8H0fV69eBdAbEuM4DqSUaLVaqNVqaLfbMAwDzWZzZP0AFWCYpolut5v4cJpcLgdd10dmJHIcB67rolqt9ssqlcpEAYRlWf31N6bd1tfX53KfJGWaoXq2bY8Ez5qmwbKssRMQEBHN21rUkKZmE4/efjuF2qy2JBOrHwK4AUACWAdQRm9tiHbE9v2kzpuyK+it1B22h8frYUAIoQshKgC4iMAKcBwH7XYb7Xa73/D0PA8XLlxINF9CTR0bVqvVoGlav1fisIYb1yo4abfbA+WqB8JxHHQ6HXieB8/zkM1moev6QHAwXOfwcdUsV2GGYfQ/V5Js24bnebh27Vrk61HBxSQqlQqklDNt3W53ZZ7QxwVcqpeKC84RUVqe/vznceLs2YGyRw8e4P4rr8S8g2aVZE7EPoDXg80P/XfUtp/gedMWNa6kH1hIKT0ppSWlrEbsR0eUrusoFov9J+2+74+M/Z+Hy5cvA+glyh7WcENQNeSHgyGVh6HyFsKbCijGUceNCnxmacgfRPVwGIYRG5yo4A9Afzz/cZdEIKfuqaiAkYhoEU6cPo1zFy6MlN+7dStibzqMJBebG80eXXFSyqh++01w2NKxUqlUkMlkcOPGjbk0isNUzsFwb8G0Jm0wqoZ2NpvFzs7ORO+p1+uoVCrY29s7sHdmHjkI6jsYFxioHA7f99FoNOaS13JYhUJh5if62Wx2YHiY+rzjJLHi+riAkYhoUdYuXcTr//pfD5Td+6OX8eiNN3DizJmUarV6uIRfgoQQDoCqlHKyKXRSlGTS8XGn6/pAo3RVF+xSORlxwsna2WwWxWIRmUwmtpdmXknMqudkkuTzq1evwnGWM+ZXOSezGL4HNzc3xyZOxy1EF6VQKPSDr0nPT0S0SE9/5jM4+cwzePj66/0y+eabuP/yy5E5EzSbIxVEBLMbeXGzH4X20wCox6Z3AbwfQHteQ4qCpOk8giTreZyDlpsKHubdeFJPl+c1Heow1dCf5Mny9vY2PM9DrVYbCDgmnZY2KWo600m+i0wmM1UwY1nWQEL2NDRNw+3bt6fKi0jqfspmswf2aqj1UA6iZrIaZ2tra+K6ERElTZw6hXMvvAA/NDMfALx+8yaDiAQlmRMxN0Fycg29YUJj+9yDAKIF4Hpo1iQbgBYcI3FSynowM1NFCNESQqzmPJHHlJqidNzrQO9p76RmnYlomqfsSVHnUp8zjqrbQT0W86br+sSBwbTDmI5qYrWanjiKurcnvadM04y9f1Wgssj7k4goytqliyNlD77+Ch7eOxYrECzEUgcRQoiiEKKB3mxH1yd8Ww29tRgGWn1SyjKArBBiboOfg2lum0EdaEVUKhXkcjlYljUQTPi+P7B+wnA+RFyislpjYRzf90ee4KsZh0zTXOh6Bupz2bY9Uqd6vd5/gq3G1A8/8Z5luJDv+/A8b6bhPGqGoDiu6x67MfuGYYzMpKVcv34dpmmO9C54nte/58Ly+Xxsr0alUukvzkhElKYzv/iLOPnsswNl8t13ce/27ZRqtHqWOoiQUlallPmgJ+HAv/pBD4CB+ICjDmBguJEQQhNCNKbYDgpCdtELVjgoeEU4jgNN01CtVpHL5Qbm/Q8HEFENe1WmVnjO5/PIZDITNbLy+Twsy+ovjlYul/trMSySruv9xdry+TxyuVx/5epCodBvZKpE5kKhMLDPLMOZtre3+8nq01LXZ7ihq3Izms1mPzBSdV/EzFppq9VqsG17IDBT07FGTYWrFgIcvjaGYWB3d3fk+tq2jU6nEzutLhHRIomTJ7H2wgsj5fs3b6ZQm9V0pHIiJmACwHAvRMgugJIQQlerRkspffTyGaYihMgCuA0gN7QCtfoLvRH6bzrCTNOEaZqoVquo1Wr9aU11XUc2m4XjOLE9A7Vard/Qbjab/VlzwitfR9E0Da1WC7Zt9xvSpmni2rVrqSStmqaJdrsN27bhum7/81uW1R8SVCqVcPfuXVSrVbiu20+uzufzE4+3T4Kmabhz5w62t7cHkn81TcPOzk7/+pVKJViWBV3Xlza5OknZbBaNRgO2bfdzatrtNlqtVuQ9deXKFTSbTVy5cmXkNcdx+lPjdjod+L4/MiMUEVHa1l68hO4f/MFA2YNvfBPvdbt4Yv3YTSqaOCGlTLsOEwl6GdoAYtdcCIY+bcZNNxs0/FvjjjFFfbLo5WgUgkBElZcA7ETVQQghAWSGgo4DbW5uymnXBXj11Vfx/PPPT/UeomXheR4ymQxarRaHxtCxwd/bRMmSjx7h+798Ae/99V8PlP/M//6/YT1Yd4kGCSFaUsqJkjyXejjTDDYwuoJ0mGrsH3pqm6C3I2qOQwvA9pi3cpgT0QHUDEAMIIiIaFbixAmsfeUrI+X7N7nwXBJWLYhYaAM9SNa+LIRwhBClYPYnKzwFbZBz4QS9JABwTQhRCXoyYgVJ5XtCiL3XXnttjp+CaPns7u5yhh8iIjq0tRdHp3R941vfwntsWx3aquVEHNQToSQWbBw0LCoY6jR11mZw3CrQG840W+2Ijqatra3IsfhERETTeN8nPoEnP/xhvPvnf/648NEj7H/1a9j41b+fXsVWQGI9EUKI3xRC/FpSx5sRhwoRrQBOE0pEREkQQmDt4uiaEfu3OKTpsJIczvRbAH4vwePNgrMhEREREVFf1JCmN1stvPtXf5VCbVZHkkFEDYAQQnwpwWNOq4MDVrQOMNggIiIiOgZOf+xjOBUxFfv+rZdSqM3qSDKIcAAIABUhxLkEjzsND+OHNKk7aHcBdSEiIiKilAkhsHZptDeCQ5oOJ7EgIlj7oADgowB+IIT4DSHEp4QQH4nakjrvEBfozYgU87oKIqZfQpeIiIiIjqS1S6N5EW995zt450c/SqE2qyHJxOqHAG4AkADWAZTRW9itHbF9P6nzDqkEP+MWycgB8MKLwxERERHRajut6zj98Y+PlHNI0+ySHM60D+D1YPND/x217Sd43r6gN6SOXo9IlMuYYbpVIiIiIjraIoc03byZQk1WQ5LDmdallBuTbjOcQhv6GWcbgDG8mJsQwgHQDC8ER0RERETHw9rF0dWr3/7e9/B2u51CbY6+ua5Yfdj8ByGEIYRoCCHa6A2NAgBHCNEOyovD7wmGKuUAXAmtJF0B0JZSxvVQEBEREdEKO/XhD+N9P//zI+X7N5lgPYvEV6wWQqwBuAbADJUBQAPAZSnlxEOZpJRNzJAEPesq0URERES0utYuXsRb3/72QNn+zZt49n/471V7lSaUaE+EEOICgDvo5SQI9PIfRLC9AKCT8joSRERERHRMRQ1peufOHbz9ve+lUJujLcnZmZ5Bb8G5dfR6AVSOxAkAGfR6J04AqAW9FUREREREC/Pkz/wMntrMjZTv/xsmWE8ryZ6IHQDPADCllP+XlPJ19YKU8o6U0gLwD9BbUfqfJHheIiIiIqKJrF0cXTNi/9YtSClTqM3RlWQQYaC3BsMfxu0gpSyjN/3rVoLnJSIiIiKayNqXvwycGGwCv/uf/hPe+s53UqrR0ZRkEJEF4E2wnxfsS0RERES0UE88+yzO/N1fHCnnkKbpJBlEuAD0CfbTg32JiIiIiBYucuG5l16CfPQohdocTUkGEbcB6EKIX4nbQQjxv6K3WNxeguclIiIiIprYWj4PPDG40sF7P/4x3nT5nHtSSQYR/xjAPoC6EOI3wjMwCSHWhBD/HL2E6i64hgMdUeVyGblcDkIIrK+vI5/Pw7Zt+L6fdtUOzfd9CCGQyWTSrgotGdd14XmTjFYlIjoaTmoanv7cZ0fKufDc5BILIoLZmC6gF0iUAXSFEHeFEA/RCxws9NaLKEyz4BzRMvA8D5lMBrZtw3Vd6Hpv5F6z2US5XMb6+vrKNrJc14UQgovwHEGe58GyLJTLZZTLZViWNVPAW6lUkMlkkMlkkM/nkc/nkcvlkMlkUK/X51Dzo6Fer8OyrP71KBQKK/t7gGgVPRM1pOmrX4V8770UanP0JLrYnJTSBfARAH+IXsCwjseLzTUBZKSU/zbJcxItgmocFItFdLtdtNttdLtddLtdOI4DANjY2Ei5lkSPua6LQqGASqWCUqmEUqkE27aRy+VmCiQMw0Cn00Gz2YTnedjc3ESj0YBpmnOo/fIrl8vQNA2VSgWNRgOtVgsbGxvHPrAiOkrOXrgAcerUQNnDu3fxxre+lVKNjpYnDt5lOkGPRAEAhBDng7I7SZ+HaFHq9Tpc14VhGKhUKgOvaZqGUqkE0zShaVpKNZwvXdfRaDTSrgZNSQUQYbquw7IsbG9vo1arTXwsTdN4D4So3kjDMAbKK5UKPM/D9vY2DMNY2d8JRKvi5NmzOPuFX8K9RnOgfP/WLTz92dGhTjQoyRWrf1MI8WvhsmCROQYQdKTt7u4C6DXK4qjhTatI0zQYhjHSYKLlpXoLor4z0zRRr9dXIo8nLZVKJbYHRuVIVavVBdeKiGYRufDc1xqQ77yTQm2OliSHM/0WgN9L8HhES4GNLTpqarVabGCrym/cuLHIKq2UZrMZOwHB5uYmgMcPH4houZ394hchnnpqoOzR66/jwTe/mVKNjo4kg4gaACGE+FKCxyRKXT6fBwA4jsOAgo6Evb29sb1jmqah1WotsEarRdM0eJ439vcBf1cQHQ0nzpzBuS+NNl33b3LhuYMkmRPhANgGUBFC5KSU9xI8NiXtHz2Tdg3m7x+9nshhTNNENpuF67rI5XJwHGei8c6+72N9fR2apqHb7Q68Vq1WYVkWisXiwLh19R5d19Fut1Eul3H9+nV4ngdd17GzsxM5jKLZbMJxHOzt7WFjYwOGYcBxnIE6qmNns1m0Wi00m83+0It2u33g51B1ijueZVn9p9ubm5twHAfZbPTi9JPUF+iNPb969Wp/ilFd12GaZj+Z/TCfa1q+78O2H89OncvlUCwWUa1WsbGxsVQJxupaxdnY2JhpFqFqtQrf9/tBiG3bhx7Kp74rz/NQq9WgaRrK5fLAdzec25G2VqvVvw7D1HWNu/eJaPmsXbo4EjTca97Go7ffxonTp1Oq1fJLcopXD72E6o8C+EGwVsSnhBAfidqSOi/RIqjhIZ7noVAo9NeImNe4Z8/zsL6+jqtXrwLoNWDVbDvN5mACWLlcRj6fR7PZ7A+lqFarOH/+fOzTUNu2kc/n4bruoZM/fd9HLpfDjRs3+udvNpuxswBNWl913Hq9Dk3T+o109f55f64wz/P6U3hWKhVUKhVsbGz0p01dtnyRSZ6CdzqdqY5XrVZRLBZRKpVQLBbhOE7/Ws+qXq8jn8/3E5Jt24ZlWTBNs3+d1f9zyybu/rp+/ToAwLKsRVaHiA7h6c9/HifOnh0oe/TgAe6/8kpKNToakkysfgjgBgCJ3tSuZQAtAO2I7ftJnZdoEdRT+FKpNLBGhGVZyGQyh2pIxdnZ2UG320Wr1UK32+0/2Qw3TlzXhW3bMAwD3W4XjUYD7XYbtVoNvu/3g5Aw13VRLpdhmmb/+IfheR42Njb65+92u7Hj7qepr6ZpcBwHUkq0Wi3UajW0220YhoFmszkyjWbSn0tRwcxwsGCaJmzbhqZpEwcslmX119yYdltfX0/0PptmuI1t2ygWiwNlmqbBsqxDNfAbjUb/mnY6HVSr1ZHejUKhMHEi+DJc32q1CtM0V3qyBaJVc+L0aZy7cGGknEOaxksyJ2IfwOvB5of+O2rjYnN0JDmOg3a7jXa7jUql0u+duHDhQqJjoNXUseF/3759G0Cv0a4aO2p4jeM46HQ68DwPnuchm81C1/XY+eoNw+gPHUnC8HATFegMDyeatr7hazB87KjE1aQ/FwBsb2/H1kWdc1KVSgVSypm2cCC5aHENYsMw4HnezOsiDA+3KxaLI+dS66/s7e0deLy0r69lWdB1farpc4loOay9OLrw3P2Xv45Hb7yRQm2OhiSHM61LKTcm3ZI6L1EadF1HsVjsPxkfHi8/D5qm9Z8GqwaVGtqkVg8Ob6qBHiXpMebDDT/VOBwOrGatb9Sxo54aJ/25XNdFvV7H5cuXI18DEDu0atklEWip733WmYh2dnYAPL4voq6lOvayJyo3m03cuHGjH+wT0dHy9Gc+g5PPDOaLyjffxL0/+qOUarT8El9sjo6IhJKOqddwzWQyuHHjxtwTQHO5HIDeE/5wAqdqjE0qyaEWkzZGZ61vvV5HpVLB3t7egQ3JpIeQqO8zasiOavguWz4E0PtODrpWSaywPi6gm+b9aiG7qGupjr3MC7f5vg/LstBqtZa6nkQUTzz5JM698AL8oZ7E/Zu38MyLL6ZUq+WWWBAhhPhNAL6U8l8kdUyio0DX9X6jLW7GlnmbZmagZRirPUl9VS6CGu5ULBaRyWRie33m8blUoKASwMMajQY0TVuK6zlsc3NzbOJ03EJ0UQqFAnzfH7ti9WHv+Waz2f//aJjqdVvGYE25cOHC2LU5iOhoWHvx0kgQ8eCVV/Dw3j2cPHcupVotLy42R5QAFTzMO4BQycKZTKbfYJlHUvc8TFvf7e3t/rSfrVYLjuOgWCwuNC/A87zY7zU8u9SkFpX4m81mD+yJmHQYluu6B87ktLW1NdGxxp0jKkhoNpv9XIlJpJFYXSgUcO3atZH7ctmHXxHRqDNbWzj57LMDZfLdd3GvyWGKUbjYHNEByuXy2MZFuVwGEP20Os4sMwepaTaBx09l1U9Vh2U3TX1VL0Ca6y/E9TTMmg+xqMTffD4fm1+i6j7pk33TNGPvV5VQfZhegnH5ECpBObwuyDiLTqxWU9JGvW9e0z8T0fyIkyex9uUvj5Tv3+IsTVGSDCIcAAK9xebY50Mro1Kp9Kf4DAcTvu+jXC73h9YM50PEJRd7nndgAyMcMKh/XwimnwtPH6nOadv2yPoRag7+ZTJNfdWY/eGZfyZtUIapxcymfTocFxiq73xZh9gYhhE7O9f169dhmuZI74pap2E4+Mjn87GzL1Uqlf5ijJMcK4oaJjXc2+G6LqrVKiqVylLmGdTrdWxtbS3tPUBEs4mapenBN76J94YWjSVg5qc2URuA/wrAIwB3AfwGgE8B+EjUluR5V33L5XJyWt/97nenfg9Fq9VqUtM0id4aKJFbpVKJfK+u6xKA1HVdFotFaRiGBCCz2awEIIvF4sD+3W63f0xN02Q2m+2/R5W12+2R+qnX1f6qvrqujxw7XBZ17uHXDyrXNG3kWJVKJfLzTVNfx3Fi9wMgDcOY6HNJKaVpmmO/pzitVksCkN1ut19WKpX63+syHZgpoAAAIABJREFUa7VaUtf1gbq32+2RMkVdI9M0R14rlUqyVquNlGWz2bHHUt/RONlsVmqaNnCvdLtdqev61N/XoqhrWywWI7dsNjtyvSbB39tE6Xv08KH8sy98UX73b318YOtcv5521RYCwJ6csH2aZGL1QxWX4PFic3EkODMUHRGmacI0TVSrVdRqtf5UpLquI5vNwnGc2ITKWq2GQqEAz/PQbDaRzWbRarUOXIVX0zTcuXOn/8Re0zQYhoFr166NPJU1TRPtdhu2bcN13X7dLMuaeCz5Ik1a31KphLt376JarcJ13X5ydT6fX1gPSzabRa1Ww/b2dv87tiwL1Wp16Z9AZ7NZNBoN2LaNTCYDoDerV9wMQleuXEGz2cSVK1dGXnMcB9VqFZZlodPpwPf9/r0cRa1IPsmq2K7r9hfuU2uAdDodVCqVpb3Gal2YcT2KTLImOprEiRNYu3gRnd///YHy/Zu3sB4x3fdxJnpBRwIHEqKLXnAwEbkia0UIIUwAW+j1vmQANKSU9dDrDoDrADwp5UyZdpubm3KShZbCXn31VTz//POznI5S5Ps+1tfXoWkauuw6TYznechkMmi1WodOzFbHchwndgE66imXy2Ovkeu6yOVyqFQqSxnwLhp/bxMthze/8x38oDAUMJw4gZ/7+st44rnn0qnUggghWlLKiZI8E+sNkFKuJ3WsoyIIIHwppR0q6wohfCmlGvBtACgFrw0foi6ljH8cTUSJcF0XmqYlMrOTyuVIa/XoVTJuCl0iorS87xOfwJMf/jDe/fM/f1z46BH2v/o1bPzq30+vYksmycTq4+gKgOH++j0A1tC/8wByQ1sZwPYC6kh07O3u7iY2NGbcwmg0HXUtGZAR0TIRQmDt0miC9f5NztIUxiDi8KJaEuHAoiWlbEopXbUB8AHszjq8iYims7W1NfWq3lF83+/nttB4tm2PHaKkriWDMSJaRmuXLo6Uvem6ePev/iqF2iynqYIIIcTaYU8ohPh7werWR56UsiClHE4g30RvzQy1T1TmnRnOmyCi+YqahnQaavXs8+fPw/d9uK6LTCbDtQBieJ6Hra2t2KlZbdvG+npvBGyz2UQmk4mdRpaIKA2nP/YxnAompQjbv/VSCrVZTlMlVgshHqE39dMvxrz+KQC6lPIPxxzj+wDOSylPTlvZZRckUSOcIxGxjwmgOU0vBBOriYhWH39vEy2X1373/8FPf/d3B8re94lP4Hy9FvOOo2+axOpZhjONZAeH/BZCT+GTJoQoBY3wg/bThBBOsJWCn3Ob+kMIYQohKsD4ACKgcxgTERER0XKLGtL01p/8Cd750Y9SqM3yORI5EUIIXQhRQ29V7LFTwwohNAAtANellLaUshw07LXgGImTUtallBZ6q3W3hBCRE4QLIUoAmlGvEREREdHyOK3rOB3RO7h/81YKtVk+Sx1ECCGKQogGerMdXZ/wbTX0pk51w4VB7kJ2nj0SUkoPvSAhLlixhutFRERERMtp7eJob8T+LQYRwJKvGh0kJVeBXm/EQfsH+xgA4oYU1YPXqqH3aJhuCFYtJlla2QVQEkJo4WFLQogsDuhFISIiIqLlsXbpIl77p/90oOzt730Pb3//+zj90Y+mVKvlsNRBxAxMABjztF818PWg1wBBQz8/7YmCoOA2gJw6VkAFDhuh/wZ6wU14PyIiIiJaYqc+9CG87+d/Hm99+9sD5fs3b+G5/+l/TKlWy2GphzPNII/Bhvsw1YhPamLyPYwuNpdFbxXr4YBhK2JfIiIiIlpiUQnW+7duYZoZTlfRqgURGxjfUFcBxujEv1MKejsaES9ZiF6JWsP4AGdAkA+yJ4TYe+2112asJREREREdxtrFi4AYnJz0nTt38Paf/mlKNVoOqxZERK9sNCdBsvbl0FSyNfSSp6NWTfIxxXAmKWVVSrkppdx87rnnkqoyEREREU3hyQ9+EE/lRhcsPe6zNM2SE5EVQuzGvKYDwJjX1T7z6v85qCdCSSzYOCDJOrxfIalzEhEREdHirF26hDf3WgNl+zdv4rn/5X+GEOOWUFtdswQRAkDugH0Oen1eNDDvgIiIiIgStPbCC/jx//F/Ao8e9cve/Yu/wFvf/jae+uQnU6xZeqYNIiZ66p4irgRNRERERIl64tln8fSn/y4efOObA+X7N28xiJiElPLX51WRhHQw2VoMDDaIiIiIaGJrly6NBhEvvYQP2CWIE6uWZnywVfvEHsbnO6gF68blbBARERERDThnGMATg8/f3/vxj/GmG7c82WpbtSDCBfqrUEdRQURzMdUhIiIiolVwUtNw9nOfGynfv3kzhdqkb9WCiErwczPm9RwAL1ilmmgi+XweQgiUy+W0q0JHhOu68DwuUE9EtGoiF5576auQ772XQm3StVJBRLBKdB1A3HSqlwHYi6sRHUeu60IIcWynfDvKPM+DZVkol8sol8uwLAu+P/0zh0qlgkwmg0wmg3w+j3w+j1wuh0wmg3o9ahkZIiI6Cs5euABx6tRA2cNOB29861sp1Sg9s0zxmhZt6GecbQAtIUQ2WFUaACCEcAA0YxaCI6JjznVdbG9vo9V6PA+453nI5XJotVrQtOmWlzEMA3t7e/A8D7quwzAM2LYNXdcPfjMRES2lk2fP4uwXvoB7jcZA+es3b+Lpz342pVqlY6mDCCGEgV7PgY7H+QyOEMJCL4m6NrzYm5TSF0LkAOwIIa4AuAsgA6AlpWQvBM2drutoDP1yoeVXKBRQqVQGynRdh2VZ2N7eRq1Wm/hYmqbxHiAiWlFrly6OBBH3Gk3I3/7tkV6KVbbUQYSUsokZkqCDnAcGDJQKTdNgGEba1aApNJtNeJ4X+b2ZpgnbtuH7/tS9EUREtHrOfuELEGfOQL7xRr/s0euv4/43voFzX/xiehVbsJXKiSAimkWtVosdZqTKb9y4scgqERHRkjpx5kxksHDv1q3FVyZFDCKI6Njb29sbm6ugadpArgQRER1vay9eGim717yNR2+9lUJt0rHUw5lofv7O//t30q7C3H3nv/lOKuf1fR/r6+vQdR3tdnukPJvNotVqwbKs/tPtzc1NOI6DbDYbecxmswnHcbC3t4eNjQ0YhgHHcUaG17iui6tXr/anGNV1HaZpwnGc2Ho0m83+cJ1wfZO4Drb9eFRhLpdDsVhEtVrFxsYGTNNM7FyHpa5VnI2NjZmmbK1Wq/1hUK1Wi4nVREQr4unPfx4nzp7Fo/v3+2WPHjzA/VdewdoLL6RYs8VhTwTRgvm+j1wuhxs3bmBzs7ekSbPZRC6Xi5xOtFwuI5/Po9ls9vevVqs4f/78wP7quPV6HZqm9Rvp6v1RbNtGPp+H67qJjvdXsxqpZOVKpYKNjY3+tKnLljMyyTSunU5nquNVq1UUi0WUSiUUi0U4jtO/1kREdLSdOHWqt4L1kP1jNKSJQQTRgnmeh42NDXS7XTQaDXS73dhx967rwrZtGIbR37/dbqNWq8H3fVy9erW/r6ZpcBwHUkq0Wi3UajW0220YhoFmszmyPoHruiiXyzBNE91uN7HhOiqYGQ4WVIKypmkTByyWZfXX3Jh2W19fT7TBPs16EbZto1gsDpRpmgbLslAoxC1jQ0RER0nUkKb7f/QyHj14kEJtFo9BBFEKhqcStSwLAEaGE6nhQI7joNPpwPM8eJ6HbDYLXddHAoNSqTRyLnXs3d3dkdcMw0CtVku0F2J7ezu2Luqck6pUKpBSzrR1u93Y4WHzFjdkyTAMeJ7HBeeIiFbA05/+NE4O/f2Ub72Fey+/nE6FFoxBBFEKhhuZqhE//LS72ezNcKxWOw5vKqA4iDp21FP54WDmsFzXRb1ex+XLlyNfAxA7tGrZJRFoqe89KqAjIqKjRTz5JM5F5D/s3zweQ5qYWH1MpZV0TJM3RlWAkM1msbOzM/Hx6/U6KpUK9vb2DhyCk3SSrwpKoobsqIBo2fIhgN53ctC12tjYSOQ8QHRAR0RER8/apYvwh4YiP3jlFTzc38fJtbWUarUYDCKIjoBJZjJSuQhquFOxWEQmkxmZJUmZxyxBKlBQCeBhjUYDmqYt5exEm5ubYxOn4xaii1IoFOD7/tgVq7loHRHRajiztYWTzz6Lhz/9ab9Mvvsu7t3+t9B+5b9MsWbzx+FMREtKNbYnfWq9vb0Nz/NQq9XQarXgOA6KxeJC8wI8z4tNnA7PLjWpRSVWZ7PZA3siJh2G5brugTM5bW1tTXQsIiJabuLkSax95Ssj5fs3b6ZQm8ViEEG0xNTT73K5fOC+qhcgzfUX4noaZs2HWFRidT6fj80vUXWftCfCNM3Yma5UQvUyDukiIqLZrF26OFL24JvfxHvdbgq1WRwGEURLTOUY2LbdDxKUer0+0ChXY/aHZ/4JLzQ3Kd/34XneVNOaAtHDmIDHs0wta+PZMIzI2a4A4Pr16zBNc6R3xfM82LY9Enzk8/nY2ZcqlQpM00xt1igiIkreU5/6FJ742Z8dLHzvPdz7Wvyw1lXAIIJoQpVKBblcLnabtsE9CV3XUavVAPQap7lcDvl8Huvr6ygUCgMNWDWVa6FQGNhvOPiYxPb2NjKZzMi6FQdxHAeu6w5ci3BDe5kbz7Varb9yt6KmY7127drI/rZto1wuj+SbGIaB3d3dkUDCtm10Op3IYxER0dElTpw4lkOamFhNNKFJplOdB9M00W63Yds2XNeF53nQdR2WZQ0saFYqlXD37l1Uq1W4rttPrs7n8wubVjWbzaJWq2F7e7s/rMmyLFSr1aXthVCy2SwajQZs20YmkwHQW7ej1WpF5nhcuXIFzWYTV65cGXnNcRxUq1VYloVOpwPf95HNZhNb0I+IiJbL2qVL6Pz+7w+UvfGtb+Hdn/wET37gAynVar6ElDLtOtABNjc35d7e3lTvefXVV/H888/PqUa06jzPQyaTQavVOnTvgTqW4zixC9AREX9vEx1lUkq0v/wVvPujHw2Uf/Af/kNs/Ne/mlKtpieEaEkpJ5oFhcOZiGiE67rQNC2R4UdqONUyD2UiIiI6DCEE1i6OJljv31rdhecYRBDRiN3d3cSGH6n1EpZ9OBMREdFhrF26NFL2puvi3b/8yxRqM38MIohoxNbW1lSrZMfxfR/NZpO9EEREtPJOf+zncCrIqQvbv/VSCrWZPwYRRDTisNOQqtWzz58/D9/34bouMpkMqtVqgrUkIiJaHkKIyDUjVnVIE2dnIqLEaZrGmYiIiOjYWbt4CT/9nd8dKHvrT/4E7/zwhzj1N/9mSrWaD/ZEEBEREREl4LR+HqcjZllbxSFNDCKIiIiIiBISOaRpBReeYxBBRERERJSQqKle3/6zP8Pb3/9+CrWZHwYRREREREQJOfWhD+F9n/z5kfL9m6uVYM0ggoiIiIgoQc9ErBmxf/MmpJQp1GY+GEQQERERESXo3Fe+AggxUPbOD36At//0T1OqUfIYRCRICKELIUoR5aYQwhFClIQQFSGEmUb9iIiIiGj+nvzgB3EmlxspX6UEawYRyaoAeH+4IAgYfCmlLaUsSyktANeEEEYqNSQiIiKiuTsXOUvTrZUZ0sQgIiFCiLjlfa8A6AyV7QGw5lsjIiIiIkrL2pe/DJwYbGq/+xd/gbe+/e2UapQsBhHJ0QF4Ma9F9ToMBxZEREREtCKeeP/78fSnPz1SvipDmhhEJCAYstSMek1KWZBSloeKNwHU5l4xIiIiIkpN5MJzt16CfPQohdoki0FEQqSU/iT7CSEcAFUpZWTQQURERESr4ZxhAE8+OVD23k9+gjdbrZRqlJwn0q7ANIKZjzwpZf2A/TQAO8E/76KX7NyWUlbnUCfzoPqo/QDkESRZJ10PIiIiIlouJzUNZz/7Wdz/+tcHyvdv3cKZra2UapWMI9ETEUydWgPgANg4YF8NQAvA9dCMSDYALThGkvXSAEzUAyGlrAczM1WEEC0hhJ5kXYiIiIho+ay9GLHw3EtfhXzvvRRqk5ylDiKEEEUhRAO9mYyuT/i2GoC6lNINFwZ5CVkhRDHBKl6edliSlNJDL3+CORFHWLlcRi6XgxAC6+vryOVyqFYT7+hKVLVahRACts2OsOPMdV14XtwcEERElLSzv/zLEKdPD5Q97HTw4N//+5RqlIylHs4UDD+qAr3eiIP2D/YxAMS1kurBa9XQezRM16CvSSmrwZSue1O8L2wXQEkIoU2aS0HLwXVdFAqFfiMsm832yy3LQqVSQa1Wg66zoymO67rIBQvwrMpc2YvgeR4cx0EmkwEAtNttOI4DTdOmOk6lUkG1WoWu6/37tNPpwPd9OI4D0+RamERESTp59izO/tIv4V6jMVC+f+sWzn7ucynV6vCWOoiYgQkAw70QIarxrgc9AiohOj/DuTYBZIQQV0JlBgA/SJ5uoDeN620AOXW+gAocNjDhcChKX7jxWyqVsLOz02/A+b6P7e1t1Ot15HI53LlzZ+rGHVEc13Wxvb2NVigRz/M85HI5tFqtqe81wzCwt7cHz/Og6zoMw4Bt2wx+iYjmZO3FSyNBxL2vNSB/+7chTp1KqVaHs2pBRB7jG+WqIW8g1Bsxi6gk7aAnxFOJ06HeiuE1IbLoJVhzTMERUigUAACO46BUKg28pmkaarUaLMtCtVpFoVBAY+iXBfXous5rM6VCoYBKpTJQpus6LMvC9vY2arXJO1M1TeP1JyJasLNf+ALEmTOQb7zRL3u0v4/73/gGzn3xi+lV7BCWOidiBhsYv4ibCjAyC6iL6hGJ+mttAdge994gH2RPCLH32muvzaV+NLlyuQzP85DNZkcCiLBKpQJN09BsNlGvHzhp18qybRtCiMhroGkaDMOAYUStwUjDms0mPM+LvF6maaJer8P32aFJRLTMTjz1FM596Usj5Ud54blVCyJSGT8ihMgKISro9XAUhRAVlcMRJHRfFkI4QohSMEOUddC0sFLKqpRyU0q5+dxzz83/Q9BYV69eBdDrhTjIzs7OwHuIDmNcjo0qv3HjxiKrREREM4haeO5+8zYevfVWCrU5vFULIsZO/xqSaLAhpXSllJaUcj3YrPBQpSAgUNPNFrjQ3NHiui583+8/QT+I6qlQ7yM6jL29vbG5CpqmDeRKEBHRcnr685/HiXPnBsoevfEG7r/ySko1OpxVy4nQMH44EwVe/fjzaVdh7p7/01cTOU6z2Yv5Njc3J36PpmnwfR/NZrM/243v+1hfX4emaeh2uwP7V6tVWJaFYrE4MvbddV1cvXq1PzWnruswTXNsr4ht26jX6+h0Otjc3IRlWbH7qnpls1m0Wi00m03Ytg3f99Fut6eqg8oJUVQeCdAbelOr1frn03W9f/wwz/Ng2zaazeZA8OY4TqKJv77vD0x3m8vlUCwWUa1WsbGxsTSzFKnrHWdjY2OmKVur1Wr/+rZaLSZWExHN2YlTp3DOMPD6v/pXA+X7N29h7YUXUqrV7FYtiOBjX0qcauhO08Da3NxEs9lEp3O4mNb3/f6MUNlsFqZpwnVdlMtluK4bmSBbKBT6uQhqCtpCoTDRDD62baNcLg+8d5o65PO9ic7C4/jVdVPHGKfZbPaPoWkastksPM9DvV6H67qRQccsPM9DPp9HpVLp9y7V63WUy2XYtj0S5KVpkt6sae4z3/dRrVZRLBYHynK5HGq1Wv97JyKi5K1dujgSRNx/+WU8evAAJ55+OqVazWbVgogOJhvSxGCDJqYaaGp+/kmohvNhhzNpmhY5G1Q+n+8nb4efmIcTuhuNRr+BrBqJ4+rjui5c14Vpmrh27Vo/6JimDqZpwjTNfo+EZVkTP9H3fb8fQITrruqmeoQOS12LnZ2dgXOYpgkhBDRNmyjgGu51mYamabh9+3ZiDfZp7rOoHgdN02BZFgqFQmKBGhERjXr605/GSU3Dw9DvbfnWW7j3Ry/jmf/ixRRrNr1Vy4nwMD7fQf3l3F1AXWhFbGz04tJpGleH7YEIi5oNSg1P2t0dvJXV8KLwE3bgcSPxIIZhoFarjTSip6nDrNTQouG6AzhwVqxpbG/3JkaLO96ks0ZVKhVIKWfaut1uak/843rUDMPo9/oQEdF8iCefxLmIoUtHcZamVQsiXKC/CnUU9deTic00MdUDMc24c7XvvMaYq0a+6w6uq6ie1l++fHmm4w7nY8xSh1mpGYbCw2yS5rou6vV65PVRn0P1hhwlSSxsqO7VpIJCIiKKtnbp0kjZg3/37/Bwfz+F2sxu1YYzVQCU0FtNOipQyKG3GNyxH86UVNLxcaCeTO/t7U38HtUgTWothHq9jkqlgr29vdihK6p80uE4UcYFPZPU4TB83597Yq8KksIJ34oKwJZt/QqVpD+O6i077HmA5IJCIiKKdmZrEyefexYPX/tpv0y++y7uNW9D+3u/kmLNprNSPRHBtKp1AKMthJ7LAOyY14giZbPZgdmWwmzbhm3bA70Uapy8ruuHfkLs+z4ymQwKhQI6nU5/9qZJ1quYVlwDfhF1UI3kJIeBRRk301aj0YCmaUs3Q9Hm5ubY63LQ7E1hhULhwJ6WJHo1iIgonjh5Emtf/spI+f6tWynUZnZHKYjQhn7G2QZgCCEGBhwLIRwAzYMWeSOKoobYhKcEBYD3v//9KJfLA0+2VeM6iUb29vY2PM9DrVZDq9WC4zgoFouR4+lV4y/pXoJp6jAr1Xsy73U1PM+L7alpNptTTeNrWRaEEDNt6+vrEz/xz2azB16XSYdgua57YKC2tbU10bGIiGh2kUOavvENvLdEswMeZKmDCCGEIYRoCCHaANRqSo4Qoh2UjwyeDoYq5QBcCa0SXQHQllLG9VAQjeU4DjRNg+u6AzPylEol6Lren/K0UCgMrKMwqbjFwtST80mPpYbiRM0adPfu3Ynrc5g6zGpc3ZMS19MwSz7EohKr8/l8bD7OtMPmTNOMvddUQvWyDeciIlpFT33qk3jiZ392sPDhQ9z76tfSqdAMljqIkFI2pZR5KWVGSilCWyYoj2xtSCn9YIVotUq0Fbcv0aRu374NoPcEOtwjUavVADxe4E3TtMj1G+J6CjzPi204q7HuwzPmxPVyqHrZtj3wpFsFObOYtg7A488avg4Hzfpz7do1AL3rOzxsrNlsjqwz4fs+PM+bqvcirqdBXbdlbECrtTairt/169dhmuZIz4pasG84+Mjn87HfQ6VSgWmaXCeCiGgBxIkTWLt4caT8SA1pmvVJGrfFbblcTk7ru9/97tTvoYO1Wi2paZoEIAHIbDYrs9ls/98AZKvVin2/rusSgNR1XRaLRWkYRv84AGSxWBzY33GcgXMZhjFwfsMwRs5hmubIe8LnKJVKA/t3u91+naLMUodGozHwHrX/Qeer1Wr99+m6PnJtoz5npVKJvd7DWq2WBCC73W6/rFQq9b+XZdVqtaSu6wP1brfbI2WKujamaY68ViqVZK1WGynLZrORx6LF4e9touPlje/8ifzu3/r44Pbx5+U7P/5xanUCsCcnbJ8udU8E0bLJZrPodrtwHAfZbBau6/ZXZlZPg8fN4lSr1aDrOjzPQ7PZhKZpaLVa2NnZidy/VCqhVCr1h1KpxOaono7wOdQwK/UkutVq4cqVKzN95lnqYBhGP4/E930YhhE7jCbMNE202+3+0CnP8/prRCSxinQ2m0WtVsP29nY/Kd6yLHQ6naXshVCy2SwajUZ/RfFyuQzHcdBqtSLzO65cuQJN0yK/c8dx0Ol0+ovLqSFcccciIqL5eN9//rfx5N/4G4OFUh6ZIU2iF3TQMtvc3JTTTC8KAK+++iqef/75OdWIojSbzX6DrN1uL90sP6vG8zxkMhm0Wq1DDcFRx4lalZtokfh7m+j4+ck/+2e4+3uDazQ99Qu/gI/8f3+QSn2EEC0p5USzjLAngighhmH0n6BHrUNAyXJdF5qmHXoMv8q/YC4AEREt2trF0Vma3vzjP8a7f/mXKdRmOgwiiBJ07dq1/rCf4elgKVm7u7uJDEFSw7KWeTgTERGtptMf+zmc+mhmpHz/1ksp1GY6q7ZiNVGqNE1LZOw+HWxra2vmPA9FLSDIXggiIkqDEAJrly7hp//37wyU79+8iff/d/9tSrWaDHsiiOhIOsx0pL7vI5fL4fz58/B9H67rIpPJzHWNCiIioihRU72+9R/+A9754Q9TqM3k2BNBRMeOmhWLiIgobafPn8fpv/083v7uqwPl+7du4dlf//WUanUw9kQQEREREaUocuG5m8u98ByDCCIiIiKiFEXN0vT2n/0Z3v6P/zGF2kyGQQQRERERUYpOfeg/w1Of/ORI+f6t5e2NYBBBRERERJSytUvRQ5qWdWFoBhFERERERCk795WvAEIMFp48iYdLOnU8Z2daYVJKiOGbkYiIls6yPmkkosV58oMfxJlcDu/+5CdYu3gRa5cu4fTHfm5p23IMIlbUyZMn8fDhQzzxBL9iIqJl9/DhQ5w8eTLtahBRyj70u7+DE888s7SBQxiHM62oM2fO4P79+2lXg4iIJnD//n2cOXMm7WoQUcpOatqRCCAABhEra21tDZ1OBw8fPky7KkRENMbDhw/R6XSwtraWdlWIiCbGsS4r6ty5c3jzzTfxwx/+EBsbGzh79ixOnjx5ZKJbIqJVJqXEw4cPcf/+fXQ6HTz99NM4d+5c2tUiIpoYg4gVJYTABz7wAdy7dw/7+/v4yU9+wl4JIqIlcvLkSZw5cwbPPvsszp07x4c8RHSkMIhYYUIIrK2tsYuciIiIiBLFnAgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpoKgwgiIiIiIpqKkFKmXQc6gBDiNQA/TLsehGcB/DTtStBS4r1B4/D+oDi8NyhOWvfG35RSPjfJjgwiiCYkhNj22kFPAAAOBElEQVSTUm6mXQ9aPrw3aBzeHxSH9wbFOQr3BoczERERERHRVBhEEBERERHRVBhEEE2umnYFaGnx3qBxeH9QHN4bFGfp7w3mRBARERER0VTYE0FERERERFNhEEFERERERFNhEEFERERERFNhEEFERERERFNhEEFERERERFNhEEFERERERFN5Iu0KEC0rIUQWgAVgA0AWgA/gupSynGrFaGkJIRwADSllM+26ULqEEEUABfR+b3SCnxUppZdqxShVQggNgBP8cyP46QG4KqX006kVLZoQogTAk1LWD9hPA7AT/PMugPcDaEspl2INCQYRRBGCBoAmpbRCZVkANSGEBSDHX/gUJoTQAZQA7KZdF0pP8Ef/NoCmlDIfKq8AaADIpFU3SlfwN2QHwHb474cQwgDQEkLw78qKC/5OOABM9B5SjttXA9ACUJBSuqHykhCiJqUszLWyE+BwJqIhwf/k2nCPQ/A/cR6Ajl4jgSjMOXgXOgZa6PU42BGvbUSU0fFxTUpZGA4Ugp5LG0AtnWrRvAkhikKIBnqBw/UJ31YDUA8HEAAQtE2ywcPOVDGIIBplIWa5+WAoQh29/4GzC60VLa3gl/mkfxhoRQXD2RA11EBKaUkp1xdfK1oGwd+L2KFswbCWzcXViBZJSlmVUuaDhwvuQfsHDzMNxP9dqaMXeKaKQQTRKAPAneB/4ii7of3omAu6nIExDQRafcF9UAJQSbsutJRUbh3RJEygPwIiyi4AfUw7ZSEYRBCN6gDQgi0Kx6xS2M6yJLlRqi4HP5lUT1H20Gv0RQaZQggTMT3gdCzlMb6toR5apfowk4nVREOklHkhhDYmwS0X/DywS5JWW5AQyWFMBPT+6ENK6Qb3hZrRTQvKOavbMSal9IUQNgBHCLGJXrKsB/SHrlxZhkRZWhob6D3QjKPaJ6lO1MAggijCATNkXAbgcxpPApBl45ACWaA/9n1gYgYhhCmEaAG4wNl3ji8pZVkI4aM35K0dBBUegA0GEDQkbiTEUuFwJqIpBF3OGoDttOtC6Qrm+ebwA1LUzEvG8Nzvwb87AK4tvFa0VIKhjypgcNC7J8Y9cabjadKZ3FINNhhEEE0o6HKuAagetEAMrTY10wqfKlOI+mMeN8yxAcBMOxGS0hU8fNCllAK92XU09NYfYkI+hbEngmjFqABi7AIxdCxcYSBJQ1RAuRfzugouOKvbMSWEqAFw1VC34GcGvXtDrSNABByRCVwYRBBNIHhKtMcAgoI1IfjUkIZ1gAPzqYDHEzPQMRL0QHjDuXRSSk9KmQNQBmAswwJitBQ6mGxIU6rBBoMIogMEv/zBAIKCtQA0NasKUciks7UdiSeMlLgdjHn4ECxCVkcwyxcdex7GD2lSwyJ3x+wzdwwiiMYIngplogKI0CJjdHxcBpAXQjSGNzxOmt0JlfMeOT52gX6+TBQtvB8dO5M8fLiOIzIWnubOBca2M1QQkeoskQwiiGIEc72PCyDY7XzMSCmrUsp81IZekiQAXA2V86nz8aFm6tqMeX0r+Mn1ZY4nb4Kk+g30EvCJVK9V3O+THJZgcg8GEUQRgqeJ+aCLOcom2BggokDwx7wKIG7YowmgzKFwx1YFvSldxymgN6SJjrng90Qdj6cDHnYZjx9cpYaLzRENCU3l2oyZdm8DvRlWmCBJYRyGQDaAlhCiNLTYXA29p4ap/9GndAQLzdWCe2E7/AQ5+JvjAHAYZB4L2tDPONvo/T7JSin7Dy2FEA6A5jLMECiklGnXgWipBCvLxo1r7gvm+aZjLviFbmDwnnHB2byOpWCo4w4ej1nWADS4sjkB/WGyw78XPPSGQXL444oKvncbvd8L4WFtXrDVgoUIh9+nfp8AwF30pgRuRe2bBgYRREREREQ0FeZEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEEBERERHRVBhEENFKEkJUhBAytGlj9jWCfbqLrOM0hBDdoI562nWZNyFEVgjRCH3mdtp1onQJIbTQ/cB7gmgJMIggouPCSbsCdDAhhAGgBcAA4AFoAtiY8L1aqJGZnWM1U3VcPmeYlNKXUq4DKARFbpr1ISIGEUS0+lRjoziuN4KWRg2ADyAjpcxJKfNB45EIAFRPXCPVWhARgwgiWnl76D3NBoCdNCtC4wVP1TUATSmlN+37pZQ+gHywTf3+o+K4fM4Y+eDnjVRrQUQMIojoWFBDmYqp1oIOop4yd2Y9gJSyGWx+uFwIYQbDf2qHquECTFLXuM95DBgAvGP4uYmWDoMIIlp5Usomek9sNSFEKe36UKyJch/oeArlfzTH7khEC8EggoiOCzv4ySFNREeTEfxkPgTREmAQQUTHgpSyjse9Eeak7wvNgjMytaoQohW8Zsa8RwuGpqj9ukKImkrwHprKtBv894FTuAbHrAkh2qFjjn2fEEIfek9DCBE5vEtNpRl6XyUom6oXJ/j8Tujzt4J/a0P7mcH5KkFRMXQNpz3nwJS+wfkkegnbAGCKwal/s0PvT/Q6BeVO6HjqmCOzKk1T13FTF0963aOONXS/toUQlZj3FIP91LSr6j6c+P+tGah8iH5PROizOsH155BFokWRUnLjxo3bym3oNUglgEqorBiUtYf2NYLybsRxZLDpEa+1gtfMmPeoOjQAtEPlLQBm6L9bode6ALSIc3VV3ce8b6SOwXvNoXOH31cbc67w+ySA4hTXPztU5/A16AIwQvvqwbVqhPavBJsx6TmHrr0W+m7jjl0JX7N5XKdQWRu94KAWKhu+b6ap68DnnOW6H3C/hj97a2h/J3Q89ZnaUfsm/P/0wP+76CXhV4KfWlCf9rzOz40bt8Et9Qpw48aN2zw2RAQRQblqYIUbsfMKIgYa9qHzyIg66KG6jTTWQ6+1wg3HoNGoGnBRDV0t1KDUh8pVYzUbcy71vmkb8lroGMPXRjW4R4IlPA7yRj7HFOeOa1yPPfa8rhOAUkRdsnH1nPQ6RL3/ENd9kvtVj/jcw8fRpr1XpvheVV0qof9fwg8I1DUtzeP83LhxG904nImIjpurwc9FLD53VYamKpW9BG+1bkUz+Ld6zcPjaSszY465LUMz00gpXQBW8E8zYsjKteCnDfSH1+joJTGra2Ahmg8gF67nhK6h16C0ZW8YWV/wbzt4/VrEe9Myl+skpSzLoZmEgu9M3Qebh614yGGve9T9qv6tA73hQ8Ex/IjP5c9wr0xKDWVSQ/4sKWX/+5BSulJKIaUsz+n8RDSEQQQRHTfV4Gc2alx6wqIaVGr60qjkUNUom2pRvKDhpt473ChVyahqyEl4a8S8R7k63FCckBoXX415vTq03zJY9HVS98GBOTBTOOx1j7pfB+7J4HP66OUWtYIcikUs4qi+ny30elnscTsT0fw9kXYFiIgWSUrpCyHK6A0zcfD4Cec8jGtYJr1I2B56Da1+ozT81BjA9pj3xtVz6qfKoQTvkSfVSvAdqIaoLmdYWC5J875OQggDQAFD30+SErrukwZCFwDcRm8IUS04v49eD0hcAHNYWfTqlw3OBwDVGYNcIkoAgwgiOo6uohdEGJPMhrQCtOHhLRM6TAPtoKfTi3h6Pa1Er1MQnNTQCx589IartdELIHcQNIgTNvfrLqV0hRDn0cvdyKP3+TQAFSFELjzMKAlBEAb0ggY7dF2d4HzumLcT0ZxwOBMRHTvB00vVWFxEbsQiqKE2e6ogNPQECxi6pc7Zf7o9ZkpRPWr/tMzxOl1Dr4FdllKuSymtIEeijkOsyh1l0dc9yH8oSynzUkqBXk8L0JueN+nAvJ8Poc6Nx/kpA+u+LGhoFRGBQQQRHV9qTLWJ2Z/OLqRhfpCg4avGrA8/lVXJ2osMltTwnrhz2kP7LYN5XCd1f1xP8JjjpHbdg8BIBeZJ/39hBOcI13skCAvWiDgOPYtES4FBBBEdS8GTWNUoGbeKtRqqMtA4EUKk1YMxMFQkNLQDeNxIDLPR+wxGVJ3VQl1zqmNxeLG44N9qQbACFkc1OtXQGLXYn/pe53Gd1NP+4e/MCNdjhrrGWch1DxKqjaEyLVTfgV6R0ExXs5xLQy8oGQ58NobPBSDPoU1Ei8MggoiOM9XoHvfkVDVeasGKvDUhRBu9BlkaQ3GKwYrHpWA14S56AU4zanrLYOjHBfQayKXQismtYHXkEhLOTwgCtHxwTic4Z0sI0UXvKbmPXoNvkUmx6nvUgpWY2+ituWEFdZ7HdQqvwN0O7p0WomfmmriucRZx3YNeryx6U62q1bBbCNafAFAfasjv4PGiebNQgcnANQs+qw/gbqheB11XIkoQgwgiOraG5uuPY6PXqFNPWtXqyueRThCRQ6+xZKEXyLjozYoTO8tU8DnPAyijV2eVCFtHr1GZaCJscM5mcM5qcM5s8LMK4Pwc1xOIq48aR++j9xTbR++7vRraJ9HrFAzxsYJj6cHxVEM/NoF7krqOee9cr3twjTJ4fI304BzqPpxH75Kq/7ACgHwQTBtznBmKiCIIKWXadSAiIqIVFfSEVLm2A9FqYU8EERERzUWQC6FhccnlRLQgDCKIiIhoXrLoLYDHhGeiFcMggoiIiObFw/hVwInoiGJOBBERERERTYU9EURERERENBUGEURERERENBUGEURERERENBUGEURERERENBUGEURERERENBUGEURERERENBUGEURERERENBUGEURERERENJX/HxfWFgD8q6nUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "USE_COLAB = False\n",
    "if not USE_COLAB:\n",
    "    plt.rc(\"text\", usetex=True)\n",
    "\n",
    "import numpy as np\n",
    "C = 10\n",
    "alpha = -0.5\n",
    "q = 0.9\n",
    "num_iter = 10\n",
    "sublinear = np.array([C * k**alpha for k in range(1, num_iter + 1)])\n",
    "linear = np.array([C * q**k for k in range(1, num_iter + 1)])\n",
    "superlinear = np.array([C * q**(k**2) for k in range(1, num_iter + 1)])\n",
    "quadratic = np.array([C * q**(2**k) for k in range(1, num_iter + 1)])\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.semilogy(np.arange(1, num_iter+1), sublinear, \n",
    "             label=r\"Sublinear, $\\alpha = -0.5$\", linewidth=5)\n",
    "plt.semilogy(np.arange(1, num_iter+1), superlinear, linewidth=5,\n",
    "             label=r\"Superlinear, $q = 0.5, p=2$\")\n",
    "plt.semilogy(np.arange(1, num_iter+1), linear, \n",
    "             label=r\"Linear, $q = 0.5$\", linewidth=5)\n",
    "plt.semilogy(np.arange(1, num_iter+1), quadratic, \n",
    "             label=r\"Quadratic, $q = 0.5$\", linewidth=5)\n",
    "plt.xlabel(\"Number of iterations, $k$\", fontsize=28)\n",
    "plt.ylabel(\"Error rate upper bound\", fontsize=28)\n",
    "plt.legend(loc=\"best\", fontsize=26)\n",
    "plt.xticks(fontsize = 28)\n",
    "_ = plt.yticks(fontsize = 28)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Значение теорем сходимости (Б.Т. Поляк Введение в оптимизацию, гл. 1, $\\S$ 6)\n",
    "1. Что дают теоремы сходимости\n",
    "     - класс задач, для которых можно рассчитывать на применимость метода (важно не завышать условия!)\n",
    "         - выпуклость\n",
    "         - гладкость\n",
    "     - качественное поведение метода\n",
    "         - существенно ли начальное приближение\n",
    "         - по какому функционалу есть сходимость\n",
    "     - оценку скорости сходимости\n",
    "         - теоретическая оценка поведения метода без проведения экспериментов\n",
    "         - определение факторов, которые влияют на сходимость (обусловленность, размерность, etc)\n",
    "         - иногда заранее можно выбрать число итераций для достижения заданной точности \n",
    "\n",
    "2. Что **НЕ** дают теоремы сходимости\n",
    "    - сходимость метода **ничего не говорит** о целесообразности его применения\n",
    "    - оценки сходимости зависят от неизвестных констант - неконструктивный характер\n",
    "    - учёт ошибок округления и точности решения вспомогательных задач\n",
    "    \n",
    "**Мораль**: нужно проявлять разумную осторожность \n",
    "\n",
    "и здравый смысл!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Классификация задач\n",
    "1. Безусловная оптимизация\n",
    "    - целевая функция липшицева\n",
    "    - градиент целевой функции липшицев\n",
    "2. Условная оптимизация\n",
    "    - многогранник\n",
    "    - множество простой структуры\n",
    "    - общего вида"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Классификация методов\n",
    "\n",
    "### Какой размер истории нужно хранить для обновления?\n",
    "\n",
    "1. Одношаговые методы \n",
    "\n",
    "$$\n",
    "x_{k+1} = \\Phi(x_k)\n",
    "$$\n",
    "\n",
    "2. Многошаговые методы\n",
    "\n",
    "$$\n",
    "x_{k+1} = \\Phi(x_k, x_{k-1}, ...)\n",
    "$$\n",
    "\n",
    "### Какой порядок поизводных нужно вычислить?\n",
    "1. Методы нулевого порядка: оракул возвращает только значение функции $f(x)$\n",
    "\n",
    "2. Методы первого порядка: оракул возвращает значение функции $f(x)$ и её градиент $f'(x)$\n",
    "\n",
    "3. Методы второго порядка: оракул возвращает значение функции $f(x)$, её градиент $f'(x)$ и гессиан $f''(x)$.\n",
    "\n",
    "**Q**: существуют ли методы более высокого порядка?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**А**: [Implementable tensor methods in unconstrained convex optimization](https://link.springer.com/content/pdf/10.1007/s10107-019-01449-1.pdf) by Y. Nesterov, 2019 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Методы спуска. Градиентный спуск"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Что такое методы спуска?\n",
    "\n",
    "Последовательность $x_k$ генерируется по правилу\n",
    "\n",
    "$$\n",
    "x_{k+1} = x_k + \\alpha_k h_k\n",
    "$$\n",
    "\n",
    "так что\n",
    "\n",
    "$$\n",
    "f(x_{k+1}) < f(x_k)\n",
    "$$\n",
    "\n",
    "Направление $h_k$ называется *направлением убывания*.\n",
    "\n",
    "**Замечание**: существуют методы, которые не требуют монотонного убывания функции от итерации к итерации."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "```python\n",
    "def DescentMethod(f, x0, epsilon, **kwargs):\n",
    "    \n",
    "    x = x0\n",
    "    \n",
    "    while StopCriterion(x, f, **kwargs) > epsilon:\n",
    "        \n",
    "        h = ComputeDescentDirection(x, f, **kwargs)\n",
    "        \n",
    "        alpha = SelectStepSize(x, h, f, **kwargs)\n",
    "        \n",
    "        x = x + alpha * h\n",
    "        \n",
    "    return x\n",
    "\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Способ 1: направление убывания\n",
    "Рассмотрим линейную аппроксимацию дифференцируемой функции $f$ вдоль некоторого направления убывания $h, \\|h\\|_2 = 1$:\n",
    "\n",
    "$$\n",
    "f(x + \\alpha h) = f(x) + \\alpha \\langle f'(x), h \\rangle + o(\\alpha)\n",
    "$$\n",
    "\n",
    "Из условия убывания\n",
    "\n",
    "$$\n",
    "f(x) + \\alpha \\langle f'(x), h \\rangle + o(\\alpha) < f(x)\n",
    "$$\n",
    "\n",
    "и переходя к пределу при $\\alpha \\rightarrow 0$:\n",
    "\n",
    "$$\n",
    "\\langle f'(x), h \\rangle \\leq 0\n",
    "$$\n",
    "\n",
    "Также из неравенства Коши-Буняковского-Шварца\n",
    "\n",
    "$$\n",
    "\\langle f'(x), h \\rangle \\geq -\\| f'(x) \\|_2 \\| h \\|_2 = -\\| f'(x) \\|_2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Таким образом, направление антиградиента \n",
    "\n",
    "$$\n",
    "h = -\\dfrac{f'(x)}{\\|f'(x)\\|_2}\n",
    "$$\n",
    "\n",
    "даёт направление **наискорейшего локального** убывания функции$~f$.\n",
    "\n",
    "В итоге метод имеет вид\n",
    "\n",
    "$$\n",
    "x_{k+1} = x_k - \\alpha f'(x_k)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Способ 2: схема Эйлера решения ОДУ\n",
    "\n",
    "Рассмотрим обыкновенное диференциальное уравнение вида:\n",
    "\n",
    "$$\n",
    "\\frac{dx}{dt} = -f'(x(t))\n",
    "$$\n",
    "\n",
    "и дискретизуем его на равномерной сетке с шагом $\\alpha$:\n",
    "\n",
    "$$\n",
    "\\frac{x_{k+1} - x_k}{\\alpha} = -f'(x_k),\n",
    "$$\n",
    "\n",
    "где $x_k \\equiv x(t_k)$ и $\\alpha = t_{k+1} - t_k$ - шаг сетки.\n",
    "\n",
    "Отсюда получаем выражение для $x_{k+1}$\n",
    "\n",
    "$$\n",
    "x_{k+1} = x_k - \\alpha f'(x_k),\n",
    "$$\n",
    "\n",
    "которое в точности совпадает с выражением для градиентного спуска.\n",
    "\n",
    "Такая схема называется явной или прямой схемой Эйлера.\n",
    "\n",
    "**Q:** какая схема называется неявной или обратной?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Способ 3: минимизация квадратичной оценки сверху \n",
    "#### (А. В. Гасников \"Метод универсального градиентного спуска\" https://arxiv.org/abs/1711.00394)\n",
    "\n",
    "Глобальная оценка сверху на функцию $f$ в точке $x_k$:\n",
    "\n",
    "$$\n",
    "f(y) \\leq f(x_k) + \\langle f'(x_k), y - x_k \\rangle + \\frac{L}{2} \\|y - x_k \\|_2^2 = g(y), \n",
    "$$\n",
    "\n",
    "где $\\lambda_{\\max}(f''(x)) \\leq L$ для всех допустимых $x$.\n",
    "\n",
    "Справа &mdash; квадратичная форма, точка минимума которой имеет аналитическое выражение:\n",
    "\n",
    "\\begin{align*}\n",
    "& g'(y^*) = 0 \\\\\n",
    "& f'(x_k) + L (y^* - x_k) = 0 \\\\\n",
    "& y^* = x_k - \\frac{1}{L}f'(x_k) = x_{k+1}\n",
    "\\end{align*}\n",
    "\n",
    "Этот способ позволяет оценить значение шага как $\\frac{1}{L}$. Однако часто константа $L$ неизвестна."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Итого: метод градиентного спуска &mdash; дёшево и сердито\n",
    "\n",
    "```python\n",
    "def GradientDescentMethod(f, x0, epsilon, **kwargs):\n",
    "    \n",
    "    x = x0\n",
    "    \n",
    "    while StopCriterion(x, f, **kwargs) > epsilon:\n",
    "        \n",
    "        h = ComputeGradient(x, f, **kwargs)\n",
    "        \n",
    "        alpha = SelectStepSize(x, h, f, **kwargs)\n",
    "        \n",
    "        x = x - alpha * h\n",
    "        \n",
    "    return x\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Как выбрать шаг $\\alpha_k$? (J. Nocedal, S. Wright Numerical Optimization, $\\S$ 3.1.)\n",
    "\n",
    "Список подходов:\n",
    "- Постоянный шаг \n",
    "\n",
    "$$\n",
    "\\alpha_k = \\overline{\\alpha}\n",
    "$$\n",
    "\n",
    "- Априорно заданная последовательность, например\n",
    "\n",
    "$$\n",
    "\\alpha_k = \\dfrac{\\overline{\\alpha}}{\\sqrt{k+1}}\n",
    "$$\n",
    "\n",
    "- Наискорейший спуск\n",
    "\n",
    "$$\n",
    "\\alpha_k = \\arg\\min_{\\alpha \\geq 0} f(x_k - \\alpha f'(x_k))\n",
    "$$\n",
    "\n",
    "- Требование **достаточного** убывания, требование **существенного** убывания и условие кривизны: для некоторых $\\beta_1, \\beta_2$, таких что $0 < \\beta_1 < \\beta_2 < 1$ найти $x_{k+1}$ такую что\n",
    "\n",
    "    - Достаточное убывание: $f(x_{k+1}) \\leq f(x_k) + \\beta_1 \\alpha_k \\langle f'(x_k), h_k \\rangle$ или\n",
    "    $ f(x_k) - f(x_{k+1}) \\geq \\beta_1 \\alpha_k \\langle f'(x_k), h_k \\rangle\n",
    "    $\n",
    "    - Существенное убывание: $f(x_{k+1}) \\geq f(x_k) + \\beta_2 \\alpha_k \\langle f'(x_k), h_k \\rangle$ или\n",
    "    $\n",
    "    f(x_k) - f(x_{k+1}) \\leq \\beta_2 \\alpha_k \\langle f'(x_k), h_k \\rangle\n",
    "    $\n",
    "    - Условие кривизны: $\\langle f'(x_{k+1}), h_k \\rangle \\geq \\beta_2 \\langle f'(x_k), h_k \\rangle$\n",
    "\n",
    "Обычно коэффициенты выбирают так: $\\beta_1 \\in (0, 0.3)$, а $\\beta_2 \\in (0.9, 1)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Анализ и мотивация подходов к выбору шага $\\alpha_k$\n",
    "- Постоянный шаг: самое простое и неэффективное решение\n",
    "- Априорно заданная последовательность: немногим лучше постоянного шага\n",
    "- Наискорейший спуск: самое лучшее решение, но применимо только если вспомогательная задача решается аналитически или ооооооочень быстро. <br></br>\n",
    "То есть почти всегда неприменимо :)\n",
    "- Требование достаточного убывания, требование существенного убывания и условие кривизны:\n",
    "    - требование достаточного убывания гарантирует, что функция в точке $x_{k+1}$ не превосходит линейной аппроксимации с коэффициентом наклона $\\beta_1$\n",
    "    - требование существенного убывания гарантирует, что функция в точке $x_{k+1}$ убывает не меньше, чем линейная аппроксимация c коэффициентом наклона $\\beta_2$\n",
    "    - условие кривизны гарантирует, что угол наклона касательной в точке $x_{k+1}$ не меньше, чем угол наклона касательной в точке $x_k$, <br></br>\n",
    "умноженный на $\\beta_2$ \n",
    "\n",
    "Требование существенного убывания и условие кривизны обеспечивают убывание функции по выбранному направлению $h_k$. Обычно выбирают одно из них."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Альтернативные названия\n",
    "- Требование достаточного убывания $\\equiv$ правило Армихо\n",
    "- Требование достаточного убывания + условие кривизны $\\equiv$ правило Вольфа\n",
    "- Требование достаточного убывания + требование существенного убывания $\\equiv$ правило Гольдштейна"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Зачем нужно условие существенного убывания?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc(\"text\", usetex=True)\n",
    "import ipywidgets as ipywidg\n",
    "import numpy as np\n",
    "import liboptpy.unconstr_solvers as methods\n",
    "import liboptpy.step_size as ss\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "/* Put everything inside the global mpl namespace */\n",
       "window.mpl = {};\n",
       "\n",
       "\n",
       "mpl.get_websocket_type = function() {\n",
       "    if (typeof(WebSocket) !== 'undefined') {\n",
       "        return WebSocket;\n",
       "    } else if (typeof(MozWebSocket) !== 'undefined') {\n",
       "        return MozWebSocket;\n",
       "    } else {\n",
       "        alert('Your browser does not have WebSocket support.' +\n",
       "              'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
       "              'Firefox 4 and 5 are also supported but you ' +\n",
       "              'have to enable WebSockets in about:config.');\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
       "    this.id = figure_id;\n",
       "\n",
       "    this.ws = websocket;\n",
       "\n",
       "    this.supports_binary = (this.ws.binaryType != undefined);\n",
       "\n",
       "    if (!this.supports_binary) {\n",
       "        var warnings = document.getElementById(\"mpl-warnings\");\n",
       "        if (warnings) {\n",
       "            warnings.style.display = 'block';\n",
       "            warnings.textContent = (\n",
       "                \"This browser does not support binary websocket messages. \" +\n",
       "                    \"Performance may be slow.\");\n",
       "        }\n",
       "    }\n",
       "\n",
       "    this.imageObj = new Image();\n",
       "\n",
       "    this.context = undefined;\n",
       "    this.message = undefined;\n",
       "    this.canvas = undefined;\n",
       "    this.rubberband_canvas = undefined;\n",
       "    this.rubberband_context = undefined;\n",
       "    this.format_dropdown = undefined;\n",
       "\n",
       "    this.image_mode = 'full';\n",
       "\n",
       "    this.root = $('<div/>');\n",
       "    this._root_extra_style(this.root)\n",
       "    this.root.attr('style', 'display: inline-block');\n",
       "\n",
       "    $(parent_element).append(this.root);\n",
       "\n",
       "    this._init_header(this);\n",
       "    this._init_canvas(this);\n",
       "    this._init_toolbar(this);\n",
       "\n",
       "    var fig = this;\n",
       "\n",
       "    this.waiting = false;\n",
       "\n",
       "    this.ws.onopen =  function () {\n",
       "            fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
       "            fig.send_message(\"send_image_mode\", {});\n",
       "            if (mpl.ratio != 1) {\n",
       "                fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n",
       "            }\n",
       "            fig.send_message(\"refresh\", {});\n",
       "        }\n",
       "\n",
       "    this.imageObj.onload = function() {\n",
       "            if (fig.image_mode == 'full') {\n",
       "                // Full images could contain transparency (where diff images\n",
       "                // almost always do), so we need to clear the canvas so that\n",
       "                // there is no ghosting.\n",
       "                fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
       "            }\n",
       "            fig.context.drawImage(fig.imageObj, 0, 0);\n",
       "        };\n",
       "\n",
       "    this.imageObj.onunload = function() {\n",
       "        fig.ws.close();\n",
       "    }\n",
       "\n",
       "    this.ws.onmessage = this._make_on_message_function(this);\n",
       "\n",
       "    this.ondownload = ondownload;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_header = function() {\n",
       "    var titlebar = $(\n",
       "        '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
       "        'ui-helper-clearfix\"/>');\n",
       "    var titletext = $(\n",
       "        '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n",
       "        'text-align: center; padding: 3px;\"/>');\n",
       "    titlebar.append(titletext)\n",
       "    this.root.append(titlebar);\n",
       "    this.header = titletext[0];\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_canvas = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var canvas_div = $('<div/>');\n",
       "\n",
       "    canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n",
       "\n",
       "    function canvas_keyboard_event(event) {\n",
       "        return fig.key_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    canvas_div.keydown('key_press', canvas_keyboard_event);\n",
       "    canvas_div.keyup('key_release', canvas_keyboard_event);\n",
       "    this.canvas_div = canvas_div\n",
       "    this._canvas_extra_style(canvas_div)\n",
       "    this.root.append(canvas_div);\n",
       "\n",
       "    var canvas = $('<canvas/>');\n",
       "    canvas.addClass('mpl-canvas');\n",
       "    canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n",
       "\n",
       "    this.canvas = canvas[0];\n",
       "    this.context = canvas[0].getContext(\"2d\");\n",
       "\n",
       "    var backingStore = this.context.backingStorePixelRatio ||\n",
       "\tthis.context.webkitBackingStorePixelRatio ||\n",
       "\tthis.context.mozBackingStorePixelRatio ||\n",
       "\tthis.context.msBackingStorePixelRatio ||\n",
       "\tthis.context.oBackingStorePixelRatio ||\n",
       "\tthis.context.backingStorePixelRatio || 1;\n",
       "\n",
       "    mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n",
       "\n",
       "    var rubberband = $('<canvas/>');\n",
       "    rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n",
       "\n",
       "    var pass_mouse_events = true;\n",
       "\n",
       "    canvas_div.resizable({\n",
       "        start: function(event, ui) {\n",
       "            pass_mouse_events = false;\n",
       "        },\n",
       "        resize: function(event, ui) {\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "        stop: function(event, ui) {\n",
       "            pass_mouse_events = true;\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "    });\n",
       "\n",
       "    function mouse_event_fn(event) {\n",
       "        if (pass_mouse_events)\n",
       "            return fig.mouse_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    rubberband.mousedown('button_press', mouse_event_fn);\n",
       "    rubberband.mouseup('button_release', mouse_event_fn);\n",
       "    // Throttle sequential mouse events to 1 every 20ms.\n",
       "    rubberband.mousemove('motion_notify', mouse_event_fn);\n",
       "\n",
       "    rubberband.mouseenter('figure_enter', mouse_event_fn);\n",
       "    rubberband.mouseleave('figure_leave', mouse_event_fn);\n",
       "\n",
       "    canvas_div.on(\"wheel\", function (event) {\n",
       "        event = event.originalEvent;\n",
       "        event['data'] = 'scroll'\n",
       "        if (event.deltaY < 0) {\n",
       "            event.step = 1;\n",
       "        } else {\n",
       "            event.step = -1;\n",
       "        }\n",
       "        mouse_event_fn(event);\n",
       "    });\n",
       "\n",
       "    canvas_div.append(canvas);\n",
       "    canvas_div.append(rubberband);\n",
       "\n",
       "    this.rubberband = rubberband;\n",
       "    this.rubberband_canvas = rubberband[0];\n",
       "    this.rubberband_context = rubberband[0].getContext(\"2d\");\n",
       "    this.rubberband_context.strokeStyle = \"#000000\";\n",
       "\n",
       "    this._resize_canvas = function(width, height) {\n",
       "        // Keep the size of the canvas, canvas container, and rubber band\n",
       "        // canvas in synch.\n",
       "        canvas_div.css('width', width)\n",
       "        canvas_div.css('height', height)\n",
       "\n",
       "        canvas.attr('width', width * mpl.ratio);\n",
       "        canvas.attr('height', height * mpl.ratio);\n",
       "        canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n",
       "\n",
       "        rubberband.attr('width', width);\n",
       "        rubberband.attr('height', height);\n",
       "    }\n",
       "\n",
       "    // Set the figure to an initial 600x600px, this will subsequently be updated\n",
       "    // upon first draw.\n",
       "    this._resize_canvas(600, 600);\n",
       "\n",
       "    // Disable right mouse context menu.\n",
       "    $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n",
       "        return false;\n",
       "    });\n",
       "\n",
       "    function set_focus () {\n",
       "        canvas.focus();\n",
       "        canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    window.setTimeout(set_focus, 100);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items) {\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) {\n",
       "            // put a spacer in here.\n",
       "            continue;\n",
       "        }\n",
       "        var button = $('<button/>');\n",
       "        button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
       "                        'ui-button-icon-only');\n",
       "        button.attr('role', 'button');\n",
       "        button.attr('aria-disabled', 'false');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "\n",
       "        var icon_img = $('<span/>');\n",
       "        icon_img.addClass('ui-button-icon-primary ui-icon');\n",
       "        icon_img.addClass(image);\n",
       "        icon_img.addClass('ui-corner-all');\n",
       "\n",
       "        var tooltip_span = $('<span/>');\n",
       "        tooltip_span.addClass('ui-button-text');\n",
       "        tooltip_span.html(tooltip);\n",
       "\n",
       "        button.append(icon_img);\n",
       "        button.append(tooltip_span);\n",
       "\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    var fmt_picker_span = $('<span/>');\n",
       "\n",
       "    var fmt_picker = $('<select/>');\n",
       "    fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
       "    fmt_picker_span.append(fmt_picker);\n",
       "    nav_element.append(fmt_picker_span);\n",
       "    this.format_dropdown = fmt_picker[0];\n",
       "\n",
       "    for (var ind in mpl.extensions) {\n",
       "        var fmt = mpl.extensions[ind];\n",
       "        var option = $(\n",
       "            '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n",
       "        fmt_picker.append(option)\n",
       "    }\n",
       "\n",
       "    // Add hover states to the ui-buttons\n",
       "    $( \".ui-button\" ).hover(\n",
       "        function() { $(this).addClass(\"ui-state-hover\");},\n",
       "        function() { $(this).removeClass(\"ui-state-hover\");}\n",
       "    );\n",
       "\n",
       "    var status_bar = $('<span class=\"mpl-message\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
       "    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
       "    // which will in turn request a refresh of the image.\n",
       "    this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_message = function(type, properties) {\n",
       "    properties['type'] = type;\n",
       "    properties['figure_id'] = this.id;\n",
       "    this.ws.send(JSON.stringify(properties));\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_draw_message = function() {\n",
       "    if (!this.waiting) {\n",
       "        this.waiting = true;\n",
       "        this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
       "    }\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    var format_dropdown = fig.format_dropdown;\n",
       "    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
       "    fig.ondownload(fig, format);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
       "    var size = msg['size'];\n",
       "    if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
       "        fig._resize_canvas(size[0], size[1]);\n",
       "        fig.send_message(\"refresh\", {});\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
       "    var x0 = msg['x0'] / mpl.ratio;\n",
       "    var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n",
       "    var x1 = msg['x1'] / mpl.ratio;\n",
       "    var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n",
       "    x0 = Math.floor(x0) + 0.5;\n",
       "    y0 = Math.floor(y0) + 0.5;\n",
       "    x1 = Math.floor(x1) + 0.5;\n",
       "    y1 = Math.floor(y1) + 0.5;\n",
       "    var min_x = Math.min(x0, x1);\n",
       "    var min_y = Math.min(y0, y1);\n",
       "    var width = Math.abs(x1 - x0);\n",
       "    var height = Math.abs(y1 - y0);\n",
       "\n",
       "    fig.rubberband_context.clearRect(\n",
       "        0, 0, fig.canvas.width, fig.canvas.height);\n",
       "\n",
       "    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
       "    // Updates the figure title.\n",
       "    fig.header.textContent = msg['label'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
       "    var cursor = msg['cursor'];\n",
       "    switch(cursor)\n",
       "    {\n",
       "    case 0:\n",
       "        cursor = 'pointer';\n",
       "        break;\n",
       "    case 1:\n",
       "        cursor = 'default';\n",
       "        break;\n",
       "    case 2:\n",
       "        cursor = 'crosshair';\n",
       "        break;\n",
       "    case 3:\n",
       "        cursor = 'move';\n",
       "        break;\n",
       "    }\n",
       "    fig.rubberband_canvas.style.cursor = cursor;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_message = function(fig, msg) {\n",
       "    fig.message.textContent = msg['message'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
       "    // Request the server to send over a new figure.\n",
       "    fig.send_draw_message();\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
       "    fig.image_mode = msg['mode'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Called whenever the canvas gets updated.\n",
       "    this.send_message(\"ack\", {});\n",
       "}\n",
       "\n",
       "// A function to construct a web socket function for onmessage handling.\n",
       "// Called in the figure constructor.\n",
       "mpl.figure.prototype._make_on_message_function = function(fig) {\n",
       "    return function socket_on_message(evt) {\n",
       "        if (evt.data instanceof Blob) {\n",
       "            /* FIXME: We get \"Resource interpreted as Image but\n",
       "             * transferred with MIME type text/plain:\" errors on\n",
       "             * Chrome.  But how to set the MIME type?  It doesn't seem\n",
       "             * to be part of the websocket stream */\n",
       "            evt.data.type = \"image/png\";\n",
       "\n",
       "            /* Free the memory for the previous frames */\n",
       "            if (fig.imageObj.src) {\n",
       "                (window.URL || window.webkitURL).revokeObjectURL(\n",
       "                    fig.imageObj.src);\n",
       "            }\n",
       "\n",
       "            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
       "                evt.data);\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "        else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
       "            fig.imageObj.src = evt.data;\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        var msg = JSON.parse(evt.data);\n",
       "        var msg_type = msg['type'];\n",
       "\n",
       "        // Call the  \"handle_{type}\" callback, which takes\n",
       "        // the figure and JSON message as its only arguments.\n",
       "        try {\n",
       "            var callback = fig[\"handle_\" + msg_type];\n",
       "        } catch (e) {\n",
       "            console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        if (callback) {\n",
       "            try {\n",
       "                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
       "                callback(fig, msg);\n",
       "            } catch (e) {\n",
       "                console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
       "            }\n",
       "        }\n",
       "    };\n",
       "}\n",
       "\n",
       "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
       "mpl.findpos = function(e) {\n",
       "    //this section is from http://www.quirksmode.org/js/events_properties.html\n",
       "    var targ;\n",
       "    if (!e)\n",
       "        e = window.event;\n",
       "    if (e.target)\n",
       "        targ = e.target;\n",
       "    else if (e.srcElement)\n",
       "        targ = e.srcElement;\n",
       "    if (targ.nodeType == 3) // defeat Safari bug\n",
       "        targ = targ.parentNode;\n",
       "\n",
       "    // jQuery normalizes the pageX and pageY\n",
       "    // pageX,Y are the mouse positions relative to the document\n",
       "    // offset() returns the position of the element relative to the document\n",
       "    var x = e.pageX - $(targ).offset().left;\n",
       "    var y = e.pageY - $(targ).offset().top;\n",
       "\n",
       "    return {\"x\": x, \"y\": y};\n",
       "};\n",
       "\n",
       "/*\n",
       " * return a copy of an object with only non-object keys\n",
       " * we need this to avoid circular references\n",
       " * http://stackoverflow.com/a/24161582/3208463\n",
       " */\n",
       "function simpleKeys (original) {\n",
       "  return Object.keys(original).reduce(function (obj, key) {\n",
       "    if (typeof original[key] !== 'object')\n",
       "        obj[key] = original[key]\n",
       "    return obj;\n",
       "  }, {});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.mouse_event = function(event, name) {\n",
       "    var canvas_pos = mpl.findpos(event)\n",
       "\n",
       "    if (name === 'button_press')\n",
       "    {\n",
       "        this.canvas.focus();\n",
       "        this.canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    var x = canvas_pos.x * mpl.ratio;\n",
       "    var y = canvas_pos.y * mpl.ratio;\n",
       "\n",
       "    this.send_message(name, {x: x, y: y, button: event.button,\n",
       "                             step: event.step,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "\n",
       "    /* This prevents the web browser from automatically changing to\n",
       "     * the text insertion cursor when the button is pressed.  We want\n",
       "     * to control all of the cursor setting manually through the\n",
       "     * 'cursor' event from matplotlib */\n",
       "    event.preventDefault();\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    // Handle any extra behaviour associated with a key event\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.key_event = function(event, name) {\n",
       "\n",
       "    // Prevent repeat events\n",
       "    if (name == 'key_press')\n",
       "    {\n",
       "        if (event.which === this._key)\n",
       "            return;\n",
       "        else\n",
       "            this._key = event.which;\n",
       "    }\n",
       "    if (name == 'key_release')\n",
       "        this._key = null;\n",
       "\n",
       "    var value = '';\n",
       "    if (event.ctrlKey && event.which != 17)\n",
       "        value += \"ctrl+\";\n",
       "    if (event.altKey && event.which != 18)\n",
       "        value += \"alt+\";\n",
       "    if (event.shiftKey && event.which != 16)\n",
       "        value += \"shift+\";\n",
       "\n",
       "    value += 'k';\n",
       "    value += event.which.toString();\n",
       "\n",
       "    this._key_event_extra(event, name);\n",
       "\n",
       "    this.send_message(name, {key: value,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
       "    if (name == 'download') {\n",
       "        this.handle_save(this, null);\n",
       "    } else {\n",
       "        this.send_message(\"toolbar_button\", {name: name});\n",
       "    }\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
       "    this.message.textContent = tooltip;\n",
       "};\n",
       "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
       "\n",
       "mpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n",
       "\n",
       "mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
       "    // Create a \"websocket\"-like object which calls the given IPython comm\n",
       "    // object with the appropriate methods. Currently this is a non binary\n",
       "    // socket, so there is still some room for performance tuning.\n",
       "    var ws = {};\n",
       "\n",
       "    ws.close = function() {\n",
       "        comm.close()\n",
       "    };\n",
       "    ws.send = function(m) {\n",
       "        //console.log('sending', m);\n",
       "        comm.send(m);\n",
       "    };\n",
       "    // Register the callback with on_msg.\n",
       "    comm.on_msg(function(msg) {\n",
       "        //console.log('receiving', msg['content']['data'], msg);\n",
       "        // Pass the mpl event to the overridden (by mpl) onmessage function.\n",
       "        ws.onmessage(msg['content']['data'])\n",
       "    });\n",
       "    return ws;\n",
       "}\n",
       "\n",
       "mpl.mpl_figure_comm = function(comm, msg) {\n",
       "    // This is the function which gets called when the mpl process\n",
       "    // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
       "\n",
       "    var id = msg.content.data.id;\n",
       "    // Get hold of the div created by the display call when the Comm\n",
       "    // socket was opened in Python.\n",
       "    var element = $(\"#\" + id);\n",
       "    var ws_proxy = comm_websocket_adapter(comm)\n",
       "\n",
       "    function ondownload(figure, format) {\n",
       "        window.open(figure.imageObj.src);\n",
       "    }\n",
       "\n",
       "    var fig = new mpl.figure(id, ws_proxy,\n",
       "                           ondownload,\n",
       "                           element.get(0));\n",
       "\n",
       "    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
       "    // web socket which is closed, not our websocket->open comm proxy.\n",
       "    ws_proxy.onopen();\n",
       "\n",
       "    fig.parent_element = element.get(0);\n",
       "    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n",
       "    if (!fig.cell_info) {\n",
       "        console.error(\"Failed to find cell for figure\", id, fig);\n",
       "        return;\n",
       "    }\n",
       "\n",
       "    var output_index = fig.cell_info[2]\n",
       "    var cell = fig.cell_info[0];\n",
       "\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.handle_close = function(fig, msg) {\n",
       "    var width = fig.canvas.width/mpl.ratio\n",
       "    fig.root.unbind('remove')\n",
       "\n",
       "    // Update the output cell to use the data from the current canvas.\n",
       "    fig.push_to_output();\n",
       "    var dataURL = fig.canvas.toDataURL();\n",
       "    // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
       "    // the notebook keyboard shortcuts fail.\n",
       "    IPython.keyboard_manager.enable()\n",
       "    $(fig.parent_element).html('<img src=\"' + dataURL + '\" width=\"' + width + '\">');\n",
       "    fig.close_ws(fig, msg);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.close_ws = function(fig, msg){\n",
       "    fig.send_message('closing', msg);\n",
       "    // fig.ws.close()\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
       "    // Turn the data on the canvas into data in the output cell.\n",
       "    var width = this.canvas.width/mpl.ratio\n",
       "    var dataURL = this.canvas.toDataURL();\n",
       "    this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Tell IPython that the notebook contents must change.\n",
       "    IPython.notebook.set_dirty(true);\n",
       "    this.send_message(\"ack\", {});\n",
       "    var fig = this;\n",
       "    // Wait a second, then push the new image to the DOM so\n",
       "    // that it is saved nicely (might be nice to debounce this).\n",
       "    setTimeout(function () { fig.push_to_output() }, 1000);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items){\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) { continue; };\n",
       "\n",
       "        var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></button>');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    // Add the status bar.\n",
       "    var status_bar = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "\n",
       "    // Add the close button to the window.\n",
       "    var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n",
       "    var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></button>');\n",
       "    button.click(function (evt) { fig.handle_close(fig, {}); } );\n",
       "    button.mouseover('Stop Interaction', toolbar_mouse_event);\n",
       "    buttongrp.append(button);\n",
       "    var titlebar = this.root.find($('.ui-dialog-titlebar'));\n",
       "    titlebar.prepend(buttongrp);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(el){\n",
       "    var fig = this\n",
       "    el.on(\"remove\", function(){\n",
       "\tfig.close_ws(fig, {});\n",
       "    });\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(el){\n",
       "    // this is important to make the div 'focusable\n",
       "    el.attr('tabindex', 0)\n",
       "    // reach out to IPython and tell the keyboard manager to turn it's self\n",
       "    // off when our div gets focus\n",
       "\n",
       "    // location in version 3\n",
       "    if (IPython.notebook.keyboard_manager) {\n",
       "        IPython.notebook.keyboard_manager.register_events(el);\n",
       "    }\n",
       "    else {\n",
       "        // location in version 2\n",
       "        IPython.keyboard_manager.register_events(el);\n",
       "    }\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    var manager = IPython.notebook.keyboard_manager;\n",
       "    if (!manager)\n",
       "        manager = IPython.keyboard_manager;\n",
       "\n",
       "    // Check for shift+enter\n",
       "    if (event.shiftKey && event.which == 13) {\n",
       "        this.canvas_div.blur();\n",
       "        event.shiftKey = false;\n",
       "        // Send a \"J\" for go to next cell\n",
       "        event.which = 74;\n",
       "        event.keyCode = 74;\n",
       "        manager.command_mode();\n",
       "        manager.handle_keydown(event);\n",
       "    }\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    fig.ondownload(fig, null);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.find_output_cell = function(html_output) {\n",
       "    // Return the cell and output element which can be found *uniquely* in the notebook.\n",
       "    // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n",
       "    // IPython event is triggered only after the cells have been serialised, which for\n",
       "    // our purposes (turning an active figure into a static one), is too late.\n",
       "    var cells = IPython.notebook.get_cells();\n",
       "    var ncells = cells.length;\n",
       "    for (var i=0; i<ncells; i++) {\n",
       "        var cell = cells[i];\n",
       "        if (cell.cell_type === 'code'){\n",
       "            for (var j=0; j<cell.output_area.outputs.length; j++) {\n",
       "                var data = cell.output_area.outputs[j];\n",
       "                if (data.data) {\n",
       "                    // IPython >= 3 moved mimebundle to data attribute of output\n",
       "                    data = data.data;\n",
       "                }\n",
       "                if (data['text/html'] == html_output) {\n",
       "                    return [cell, data, j];\n",
       "                }\n",
       "            }\n",
       "        }\n",
       "    }\n",
       "}\n",
       "\n",
       "// Register the function which deals with the matplotlib target/channel.\n",
       "// The kernel may be null if the page has been refreshed.\n",
       "if (IPython.notebook.kernel != null) {\n",
       "    IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAYAAACPNyggAAAAAXNSR0IArs4c6QAAIABJREFUeF7tnQt0XFd57z+932/J8kO2YymO7TiJ41ecx4UFxA6PllKKbaANlMeNVQptaQCZdAEtKTTXprRAcYsd0sdtu9rabpvAhULkFCjk4UiW83JsJ0hOItmyrLdkSdb7rv8en/E5M2c0M9KcmbPP+Z+1tDQanTl779+3tf/69v72t9NmZ2dnhRcJkAAJkAAJkEBSCaRRgJPKm4WRAAmQAAmQgCJAAWZHIAESIAESIIEUEKAApwA6iyQBEiABEiABCjD7AAmQAAmQAAmkgAAFOAXQWSQJkAAJkAAJUIDZB0iABEiABEggBQQowCmAziJJgARIgARIgALMPkACJEACJEACKSBAAU4BdBZJAiRAAiRAAhRg9gESIAESIAESSAEBCnAKoLNIEiABEiABEqAAsw+QAAmQAAmQQAoIUIBTAJ1FkgAJkAAJkAAFmH2ABEiABEiABFJAgAKcAugskgRIgARIgAQowOwDJEACJEACJJACAhTgFEBnkSRAAiRAAiRAAWYfIAESIAESIIEUEKAApwA6iyQBEiABEiABCjD7AAmQAAmQAAmkgAAFOAXQWSQJkAAJkAAJUIDZB0iABEiABEggBQQowCmAziJJgARIgARIgALMPkACJEACJEACKSBAAU4BdBZJAiRAAiRAAhRg9gESIAESIAESSAEBCnAKoLNIEiABEiABEqAAsw+QAAmQAAmQQAoIUIBTAJ1FkgAJkAAJkAAFmH2ABEiABEiABFJAgAKcAugskgRIgARIgAQowOwDJEACJEACJJACAhTgFEBnkSRAAiRAAiRAAWYfIAESIAESIIEUEKAApwA6iyQBEiABEiABCjD7AAmQAAmQAAmkgAAFOAXQWSQJkAAJkAAJUIDZB0iABEiABEggBQQowCmAziJJgARIgARIgALMPkACJEACJEACKSBAAU4BdBZJAiRAAiRAAhRg9gESIAESIAESSAEB3wlwcXGx1NTUpAB1coqclVkZnRyVkYkR9XVl6sq1gtNE8jLzpDC7UAqyC9RrXiRAAiSgK4GOjg4ZGhrStfriWgE+dOiQgnrixAk5ePCgen306FEpLS2VlpYWaWhoiPjeXNa48cYb5eWXX9bWYKEVn52dldM9p+Xx1sfV189e/5kSYONaXb5a7qm7R3295bq3SHFOsWfazoaQAAn4m4Du47krBfjYsWNSW1urvvbv36962Pbt26W5uVn27Nkj+H1bW5ts2bIl7D383usC3DPaI8fajgVF9/zw+WCTS3NLZXvtdrmn9h7ZUbdDriu9zt9/oWw9CZCAZwlQgB0wLbzfgYEB5eXC621sbJS6ujolyDt37lTiu3fvXtm6dWvYe0eOHPGcAI9PjcvTHU8HBbels0Uw1YwrIy1D7lh+hxJceLlblm6RjPQMB6zCR5IACZCAuwhQgB22R319vezYsUOJ8K5du5QnDHHGawhy6Hu4T3cPGNPKZ3rOBAS37XH56Ws/tUwrX19+fVBwMa1cklvisBX4eBIgATMB/I3ySh6BtLQ028IowA7aAJ7uvn371BowhNgQW7yPn80CbLwXKsAHDhwQfBlXf3+/dHZ2Oljr+T0a08pPtD0RFN2OoY7gg0pySuTu2ruD08q1ZbXzK4SfIgESWBCByclJaW9vl/Hx8QU9hx+Oj0BOTo4sX75csrKyLB+kAMfHMa67Mc0MAcaFtWAvTUFPTE/I0+1Xp5XbHpcTF05YppW31WwLerlbl22VzPTMuNjxZhIggcQTwD/6RUVFUlFRIZG8ssSX6u8nYraht7dXhoeHlQaYLwqwQ30Da7+YbkbUM4KuysvL1XesC5uDsELfc2sQFjrRK72vBD3cn5z7iYxMjgTp1ZXVBaOV33rdWzmt7FC/4mNJYL4E1NLQmTOyZs0aSU9Pn+9j+Ll5EJiZmZGzZ8/K2rVrLf/4UIDnATPaRyCqmG6G6Pb19ckDDzyghBfBWfgPyLwNye69uZ6fTIP1jfVZppXfGHwjWDVsB7p71d1KdHfU7pC68rpoWPh7EiCBFBIwBDhUBFJYJd8UHYl9MsdzJ2C7chuSEw01numkwTCt/EzHM8Fo5eYLzcFp5fS0dNm2bFvQy71t2W2cVnbS0Hw2CSSYAAU4wUDjeBwFOA5Ybr41kQKMTvFq36tBwf3Jaz+RyxOXg83HHty3171die7bVr1NsEeXFwmQgJ4E3CjA2BHy0EMPqTVpIzkRZgixdBe6XmpQx2ewlr1p0yZHDIFnowzkbUCuhkSUQwF2xFTJf+hCBXh4fFh+3PrjoOi+Pvh6sBFF2UUqWhlTyhBdrOsyUCP5NmaJJOAEATcKMAJV3//+96vmGkKH+BnkS5jriuWe+TI04nfwedTPyGQ43+fhcxTghdBz0WcXKsAvXHxBNhzcoFqUJmmyomSF+qrMq5TszGwZHB+UgSsD8q13fEsQvcyLBEjAGwTcKMBGjgSDMGJidu/erTzguS54qIcPH1aZBZ26UBd6wHPT5RpwnL0Pf4Qffeyj8mLXi9JyscX20x/e8GH57ru/K1kZ1j1rcRbF20mABFxEIEyAP/5xkVOnEl/D9etFHnkk6nMRrIp8CNiqaXi8+NnsccIbRSArAlqNaWfsLsGFQNdomQOjVSLS81E3XEZZ0Z4T7ff0gKMR0uT3C/WAzc1sbG2U+x+/X1669FJY6/Oz8gWBVnfW3Cl3rbhLbq+5XcrzyjWhxGqSAAmEEnCbAMOLDZ3iNYsqBBfCi7VYrAtDpPGz4R3bCTA+A1G1u/B587pypOejPPxTgHs3b96cEC+bAuyRv8dECjCQTM1MySMtj8gXf/JF6R7tlpyMHBVwhdzNmIo2X+sq18mdy++Uu5bfpb7fUHED14g90q/YDO8TcNsUtJEPwTyNHOoBwyrGOnFoMFQiPOC5np/IHkEBTiTNFD4r0QJsNGXwyqA89IuH5BvPfEPO/cE5qS6sVvmcn3zjSXmq4yl5qv0plYjDfFXkVSghNr62Lt0qeVk8ozeF3YNFk0BEAm4TYGQHhKiap3ntPGC8hxS98ILNIpwoDzjS8xPZlSjAiaSZwmc5JcBGk871n1N7f+3yNXePdCvPGGKMr6YLTXJl6kqQBtJNblqySU1bK095xV2ytGhpCmmxaBIgAYOA2wTYzts1v2ek721qalIH2uAyi7Xd5+OxdrTnx/OsaPdSgKMR0uT3TgtwPBiQuONk50klxk+2P6m+Ll6+aHnEypKVlmnrm6tvZgKPeCDzXhJIEAEdBNhNUdAJwq4eQwFOJM0UPstNAhyKAZ0M+4rVtDW85I6n5IWuF2RmdiZ4a0FWgeCgBmMdGcFdTPCRwg7Fon1DwC0CjLVfTClj/69dkotY9vjGco+bDEsBdpM1FlAXNwuwXbOQ+OP4+ePBaWtMYQ+NDwVvxV7k9YvWW6atmQBkAR2EHyWBCATcIsDRDJTqTFjR6jef31OA50PNhZ/RTYBDEU7PTMvL3S8HPWR4y639rZbbqvKrLNPWm5dultzMXBdag1UiAX0I6CLA+hCNvaYU4NhZufpO3QXYDm7X5S4V3GVEXOMQCKwvG1dWepZAhI1pawR4LS5c7Go7sXIk4DYCFODUWYQCnDr2CS3ZiwIcCmh8alxOdJ4ITlsjuOvSyCXLbatKV6koayPi+qZFN0lGekZCWfNhJOAlAhTg1FmTApw69gkt2Q8CHAoMnbetvy0oyAjuQipNbJcyLhwkgYAuI1EIAr1wZjEvEiCBAAEKcOp6AgU4dewTWrIfBdgOIBKHGMFd8JBxjrH5KEUEd2HLk5FKE8IMr5mnOyW0O/JhGhGgAKfOWBTg1LFPaMkUYHucCO5CTmuIsZEo5NzAOcvNWDdWWbuuJgpB0pCczJyE2ocPIwG3EnCLAOOfYKSfxCELyAeNvMuJOHM32dyRSxppMmM5EIICnGzrOFQeBTh2sJ3DnZZp6xMXTsjkzGTwAch7vWXpFks6zUUFi2IvgHeSgEYE7EQg7ctpCW/B7B9fWxqye3hZWZn09/erX0GA7777biVi5oMSEl6pqw9c6P5h8+cpwCI8jtCpnurB5yJtJiKsDQ8Z3nLPaI+lpdeXX2/ZAnVj1Y2SnpbuQRpskt8IuFGAYQMk5oAAJ+Lg+2g2DT1/ONr9ob83f54CTAGOt//wfhMBDEi/7PtlMJUmhPlUt/V81JKcEhXcZWyBwhGNRTlF5EgC2hFwqwDDC8axf62treqMX/wMoWtoaFDiDGHGMYR4D0cKIgc0xA8XfgfPOfRz+D3eQ+5oeK3IvIV7keoS7xnHDeIZxr1btmwJlofX9913nyoDz8f9mG42fx6/QxnG51BGJC+eU9Da/bnYV5hT0M4asn+sXwV0Gak0j3ccl5HJkWCh8IY3VG+wTFsj3zWDu5y1C5++cAJuF2AIHMQXwguRNUTPEFIQgAAa9xhEzO8Zn8P6MgQbog4BhsjiuRD6EydOWGDid/gcBNT82rgJn3niiSfUOcTmz+NeuzLsLEUBXnj/dcUTKMDJNQPOS0Y+a2PaGt+R79p84cQnI7gLe5NvXXyrZGdkJ7eiLI0EohBwqwCbvVwIJ7xdfH/ggQeU6EGYIXYPP/yw8mLNgowmQzwjfQ5T20buaQh6JAE2B1PV1dUp4cZpSb29vUEPGt5tqAAbnzOXQQH28J8iBTj1xu0Y6pCn2wPHMmId+eTFkwKhNi6kzcTZyMa09R3L75DK/MrUV5w18DUBtwowRM0QShgIEdKhl7HeanicEFOzBxz6OfP6bCwCbPaAIfA4KAKCiylvlGlMdVOArZZxdRCW+QBpdAj8Z2WsEeA/O0yFGFMmhw8flt27d6v/+Oa6KMDuG0NHJ0dVcJeRShPC3DfWZ6noDRU3WFJprq1cy+Au95nS0zVyiwCHbkOCF2mc82us72IchMeL6WVMDRuebqQ14NDPwSO2806NtWJzmcZUMupgBIShPMPTNgTcvNaMz+NnesDoVS68YDR0Ckxl4MLPRieDwfGFnxGSjwv/0dn95xfaNAqwC40dUiV0yVd6X7HsST7dc9pyV1lumcAzNvYkI7irILvA/Y1jDbUl4JZ9wG4DGE8083zrzjXg+ZJbwOcihbzjvzpDbOPdl0YBXoBBUvjR3tHesOCusamxYI0y0jLU2rGRShPfl5csT2GNWbTXCFCA7S1KAZ5/T3f1FLSdAMPYiOIzMr9goR+vMV2NKL1oFwU4GiE9fj85PSnPdz0fXEfG9PX54fOWytcU11gEGdHXWRlZejSQtXQdAQpw6kxCDzgF7O0EGGsG5gACo1p24fV2VaYAp8CQSSqyfbDdsif5uYvPyfTsdLD0/Kx8wVS1MW2NKezyvPIk1Y7F6E6AApw6C1KAU8DeToDNUXSYfkawAKajI2WDOXDggODLuJDCrbOzMwWtYZHJJjAyMSLPnn82uCcZwV0DVwYs1VhXuS64JxlR1wj24p7kZFtKj/IMEVizZo2kpzO7WzKtNjMzI2fPnpW1a9da/j51d6i0moI28p4aG8Ex7YxIOkT8wQPGFS0QS3eDJbPTe62smdkZOdNzxpJKE8Fe5qsir0IFdxlboJDrGp4zLxIAASyBFRUVSUVFBf9RS1KXwD8+2E88PDwclilL9/HctQIM7xbTzYiEhqhCZCHASG9mPj0DwouN501NTbZT06F9RHeDJanP+6YY5LLGnmTjFKimC02CnNfGlZmeKRsXb7RsgVpWvMw3fNhQK4HJyUlpb2+X8fFxokkigZycHFm+fLlkZVljOHQfz10rwE7ZVneDOcWFzw0QmJieEKwdG3uS8b3zsnXJAqkzVeau5XcqYca5yRBqXv4h4NLdm541QKRlId3HcwqwZ7ssG5YIAhhokTrTnEoT0deYzjaugqwC2VazTQV3IZUmDp8ozZ07IUwi6sZnkIDfCVCANesBuhtMM9yerO7w+HAwuAtT1093PC1D40PBtqZJmuAYRmMdGZ4yjmlkcJcnuwMblUICuo/n9IBT2HlYtDcIwBt+uftlyxYoHNNovqryqywnQCG4CzmveZEACcyfAAV4/uxS8kndDZYSaCw0bgKXRi5Zpq2R63p8+lrgTlZ6lmxeujk4bQ0veXHh4rjL4QdIwM8EdB/P6QH7ufey7UkjMD41Li2dLcE9yQju6hrpspS/qnSVJXPXTYtukoz0jKTVkQWRgG4EKMCaWUx3g2mGm9WNQADBXecGzln2JL/Y9aLMyrWzUYqyi1Rwl7GWvG3ZNinJLSFTEiCBqwR0H8/pAbMrk4BLCCCQ63jH8eCe5Gc6npHhieFg7RDchS1PRipNRFzDa2Zwl0sMyGoknQAFOOnIF1ag7gZbWOv5aZ0ITM9My0uXXrKk0mzrb7M0obqg2jJtvWnJJsnJzNGpmawrCcybgO7jOT3geZueHySB5BPoHO5U256wLxlboE5cOCGTM5PBimRnZAsirM1boBYVLEp+RVkiCSSBAAU4CZATWYTuBkskCz5LfwJImwkRNlJpQpi7R7stDcMeZJW5qyaQvWv9ovWSnsbDBPS3Plug+3hOD5h9mAQ8RADBXa39rYFUmu1PyVMdT8mpS6cswV0lOSUqW5eRShNHNBblFHmIApviFwIUYM0srbvBNMPN6rqAAI5gRECXkU4Tr0cmR4I1gzd8S/Utlmlr5LtmcJcLjMcqzElA9/GcHjA7OAn4jMDUzJRgy5OxjozvyHdtvpYWLbVMW29cslGwvsyLBNxEgALsJmvEUBfdDRZDE3kLCcRN4PzQeRXcZZwChaQhEGrjQtrMrUu3WtJpVuZXxl0OP0ACiSSg+3hODziRvYHPIgGPEBibHBOcjWw+Bap3rNfSuhsqbrBsgVpbuZbBXR6xvy7NoADrYqmr9dTdYJrhZnU9QgDBXa/0vmLJ3HW657SldWW5ZXLH8juC0dYI7irILvAIATbDjQR0H8/pAbuxV7FOJKABgb6xvmBwF7ZBPXv+WRmdHA3WPCMtQ25dfGtw2hp7k5eXLNegZayiLgQowLpYih6wZpZidXUjMDk9KS90vWDZk9w+1G5pRk1xTTC4C6k0N1RvkKyMLN2ayvq6hAAF2CWGiLUauhss1nbyPhJwA4H2wfZr68gdT8nJzpMyPTsdrFpeZp5gqtrI3IUp7PK8cjdUnXXQgIDu4zmnoDXoZKwiCXiFwMjESDC4C9PWT7c/Lf1X+i3NQzCXOZXmmoo13JPslQ6Q4HZQgBMM1OnH6W4wp/nw+SSQTAIzszNytuesZU/y2d6zlirAIzZSaWLaGrmu87Pyk1lNluVSArqP5/SAXdqxWC0S8CuBntEe5RkbqTQR3IWc18aVmZ4pGxdvtGyBWla8zK+4fN1uCrBm5tfdYJrhZnVJYMEEJqYn5LmLz1m2QF0YvmB57oqSFZZpa6TWhFDz8jYB3cdzV3vALS0tsmnTpmAPqq+vl3379snhw4dl9+7dUlpaKkePHlXfcW9DQ0PU3qa7waI2kDeQgMcJYE/yG4NvWAT5+a7nBdPZxlWQVSDbarYF9yQjuKs0t9TjZPzXPN3Hc9cK8LFjxwSC29raGuxVZWVl6jVEeM+ePUp0m5ub1Wvc39bWpl7PdeluMP/9ibHFJBCdwOWJy2ofspFKE1PYg+ODlg+ur1pvmbbGMY08cCI6Wzffoft47loBhtF37NghjY2NQfvD2925c2fw5/3790ttba16D+K7d+9eOXLkCAXYzX8xrBsJJIEAvOGXu1+2pNJ8te9VS8nIZW0cyYjvm5dslrysvCTUjkUkigAFOFEkbZ4TKsAQXExJG9PN8JB37dol27dvl4GBAfXaLNh2VdPdYA7i5qNJwNMELo1cCgZ3YQtU84VmGZ8eD7Y5Kz1LNi/dHJy2higvKVriaSa6N0738VwrD9joLIcOHVKCi+lpQ4DhAUOQKcC6/0mx/iSQHALjU+Ny8uLJ4BYoTF93jXRZCl9VusqSSvOmRTdJRnpGcirIUqISoABHRTT/G8weMKaf+/r6guu9mGquq6uLOgV94MABwZdx9ff3S2dn5/wrxU+SAAl4kgCCu14beM2SShOpNWdlNtjewuxCub3m9mDE9bZl26Qkt8STPHRoFAXYQSuZBRjTzljvRcQzPGBcW7ZsUcFXiH5mEJaDhuCjScCnBIbGh+R4x/HgnmQEdw1PDAdppEmawCs2Z+6qLatlcFeS+gsF2CHQ8HgRVIVpZUQ2G8JbXl4uTU1NKhIaF8QYwsxtSA4Zgo8lARIIEpiemZZT3acsmbva+tsshKoLqoPT1kZwV05mDik6QIAC7ABUJx+pu8GcZMNnkwAJxE/g4uWLKrgLgV3I3nWi84QgeYhxZWdkq/SZd9bcKUileUfNHVJdWB1/QfxEGAHdx3NXB2E50d90N5gTTPhMEiCBxBFA2swTF04Ep60R3NU92m0poK6sTokxRBle8o1VNzK4ax4m0H08pwDPw+j8CAmQAAnESgDBXa39rZbMXacunbIEdxXnFCvPWB06sfxOQXBXUU5RrEX49j4KsGam191gmuFmdUmABGwIDFwZCAZ3Yer6mY5nZGRyJHhnelq6IJ+14SHDW15ZspLBXSEsdR/P6QFzeCABEiCBFBOYmpmSly69FEylibVkbIkyX0sKl1gyd21cslGwvuzniwKsmfV1N5hmuFldEiCBeRLAiU/qSMarXy2dLTI5Mxl8Wm5mrmxdujU4bY0p7KqCqnmWpufHdB/P6QHr2e9YaxIgAZ8RGJscU+kzIchGxHXvWK+Fwg0VNwQE+WrE9drKtYLpbK9eFGDNLKu7wTTDzeqSAAk4RADBXThgQgnyG0/KUx1PqQMozBeOYDSCu5As5LZlt0lBdoFDNUr+Y3Ufz+kBJ7/PsEQSIAEScIRA/1i/PN3xdHDa+vj54zI6ORosKyMtQzYs3mDJ3LWiZIUjdUnGQynAyaCcwDJ0N1gCUfBRJEACHieA4K7nLz5v2ZPcPtRuaXVNcU1w2hrT17cuvlWyMrK0IKP7eE4PWItuxkqSAAmQQGIItA+2Ky/ZmLY+2XlSpmengw/Py8xTU9XGWcl3LL9DyvPKE1N4gp9CAU4wUKcfp7vBnObD55MACfiLwMjEiAruMgK7sKbcf6XfAgHBXEZgF4R5TcUaV+xJ1n08pwfsr781tpYESIAE5iQwMzsjr/S+EvCQsQWq4yk503PG8hl4xEa0Nb5vXbZV8rPyk06WApx05AsrUHeDLaz1/DQJkAAJxE+gZ7RHZesytkA9e/5ZQc5r48pMz5SNizcG9yQj4npZ8bL4C4rzE7qP5/SA4zQ4bycBEiABvxPAaU9GcBemrvGFxCHmC9HVxjoyviO1JoQ6kRcFOJE0k/As3Q2WBEQsggRIgATiIoA9yYiuhofc2NYoTeeb1J5kc3BXQVaBCu6CdwxBfvPKNy94T7Lu4zk94Li6GW8mARIgARKYiwCirDce3CjYk1yeXy65GbkqhSYOoBifHg9+9NQnTsmNi25cEEwK8ILwJf/Duhss+cRYIgmQAAlEJ4DgrXP95+TkxZPycMvD8njr42EfSpM02VazTVaVrpJ/+o1/WnCaTN3Hc3rA0fsV7yABEiABEjARwBowppixhxiC+9zF59TX8MRwRE47b9wpD77lQVlXtS5hLCnACUOZnAfpbrDkUGIpJEACJBAgMDQ+FBRYQ2xPXTplOZkJ9xVlF6ksWvi6vvx6+YMf/YH6/LtWv0u+8tavCI5PTPSl+3hODzjRPYLPIwESIAENCSCQqvNyp/Jq4c0aYtva3xrWmqVFS5XQYuuR8X1V2arglPJjZx6Tbx7/pnzlbV9RAVdOXRRgp8g69FzdDeYQFj6WBEjARwSmZ6bll32/DIqsIbaXRi5ZKGDNFkccmsUWr6sLq+ekhYArnMTk9KX7eE4P2OkewueTAAmQQAoJIGHGS5desni2L3S9ICOTI5Za5WTkyM3VN8ut1beq6WJ4t/i5MLswhbWfu2gKsGtNY18x3Q2mGW5WlwRIIIkEsPXHPH0Mz/Z092nLflxUB96pefoYgov8zrqcgmQg1X08pwecxD8OFkUCJEACiSBgJL5QYmuKRH598PWwxy8vXh70aI2pZGSpSktLS0RVUvoMCrCD+FtaWmTTpk3BEg4dOqRenzhxQg4ePKhe19fXy759++Tw4cOye/duKS2de91Bd4M5iJuPJgEScCEBnOl7tuds2Hpt31ifpbbpaemCU4vg2RreLQS3Ir/Cha1KTJV0H89d6wEfO3ZMiWtrayACDz/X1taqr/3796v3GhoapKysTL2GCO/ZsyeqVXU3WNQG8gYSIAFtCYxOjgrWZ82RyC9eetFy8AEahzN7NyzeEFyvhdDevOhmycvK07bt86m47uO5awUYxtixY4c0NjYqu8D7HRgYUKJ79OhR9T68YLzeuXNnzLbT3WAxN5Q3kgAJuJpA90h32HotjgFERinzVZlfGbZeu7p8tWSkZ7i6fcmonO7juTYCbDYmPGOIM4QX3jCmqTFdDXGOdulusGjt4+9JgATcRQDrtecGzgXEFntsuwLfzw+fD6soUjQaEcjGei323HphvdYJq+g+nmsnwG1tbWq62VgDNoxq9pDnMrTuBnOiE/OZJEACiSEwOT2pUjSaI5HxenB80FIAjuW7sepGy3otppSTsXc2MS11x1N0H8+1E+C9e/cqAcaF6ee+vj619os14iNHjoQJ84EDBwRfxtXf3y+dnZ3u6D2sBQmQgLYEhseH5fmu54OeLbb8nOo+JciTbL6wj1alaLy6vxav11etl5zMHG3b7paKU4AdtIR5DdgQ3O3bt6tIZwhueXm5CsrCz0aEdLRALN0N5iBuPpoESCACgYuXL1q2+0BskUkq9FpcuNiyXguxrSuvW/CpPzSMPQHdx3PXesDwbuHtYr0Xotrc3Cy7du1Soguv94EHHlBrvhBevNfU1BT0jOfqrLobjH+IJEACzhFAABSENXQfZh3BAAAgAElEQVS9tmukK6xQBEJhvdbs2UKAeSWPgO7juWsF2CkT6m4wp7jwuSTgNwLjU+MqRaN5vRZTypcnLltQZGdky02LbrKs195SfYsU5RT5DZnr2qv7eE4Bdl2XYoVIgAQSTQCHAxhn1mL6GFHIp3tOC5JcmK+SnJIwr3Zd5TrtUjQmmp9bn0cBdqtlItRLd4NphpvVJYGkEsCWH2zvMSeygOC+NvBaWD1qimvC1muvK72OW36SarGFFab7eE4PeGH256dJgARSRABH6iFxheHRYn8tvNye0R5LjZCiEQcNhK7XIsEFL70JUIA1s5/uBtMMN6tLAgkhgBSNxpF6xtm1SNk4NjVmeX5uZq5gfdZ80g+O1MvPyk9IPfgQdxHQfTynB+yu/sTakIDvCfSO9oYdPHCm50xYisbyvHJLYBQ8XBwejyQXvPxBgAKsmZ11N5hmuFldEohIAOu1OD7PvF6LKeT2ofawz2Bt1uzVYn8t1nCZotHfHUz38ZwesL/7L1tPAkkhgBSN8GKN6WPjO6KTzVdGWkYgRaNpf+2G6g1Slhc49YwXCZgJUIA16w+6G0wz3KyuDwlgH23okXpYvx2fHrfQKMgqUEfqmT3b9YvWC9ZxeZFALAR0H8/pAcdiZd5DAiRgS6DrclcwkYXh1b7a+6rMyqzl/kUFi8LWa68vv54pGtmvFkSAArwgfMn/sO4GSz4xlkgCogKg2vrbLAcPYL2283L4wSYQVuMoPcO7XVK0hBhJIOEEdB/P6QEnvEvwgSSgNwGc5nPq0inLeu3zF5+X4YlhS8Oy0rOCKRqV4C7ZqLYAFecU6w2AtdeGAAVYG1MFKqq7wTTDzeq6nMDglUF1pJ75oHicZzs5M2mpOUTV8GqN7+uq1gnyJPMigVQR0H08pwecqp7DckkgiQSw5efC8IWw9VpMK4dey4qWWcV2yUbBNiBklOJFAm4iQAF2kzViqIvuBouhibzF5wSQovHVvlfD1mu7R7stZNIkTdZUrrGs1yIqGQFTvEhABwK6j+f0gHXoZawjCUQgMDY5FkjRePFk0LvFFiCkbjRf2Npz86Kbr4ntko3q54LsArIlAW0JUIA1M53uBtMMN6ubQAJ9Y31hB8UjucX07LSllLLcMhUQZd5fC0+XKRoTaAw+yhUEdB/P6QG7ohuxEiRwjQDWa98YfMNyUDw8XLwXeq0sWRm2Xru8eDlTNLJD+YIABVgzM+tuMM1ws7pRCOBAeHix2FNrRCLjNbxd84UUjYg6Nkci4zUOJOBFAn4loPt4Tg/Yrz2X7U46gZGJEZWiUYntxZPq68WuF8NSNOLoPOQ/DiazWLJR1letl7ysvKTXmQWSgJsJUIDdbB2buuluMM1w+7a63SPdlsAoeLc4PD40RWNVflXYei0ySWWkZ/iWHRtOArES0H08pwccq6V5HwnYEMB67bmBc2r62ByJjD23oVddWV3Yeu2SwiVcr2XPIoF5EqAAzxNcqj6mu8FSxY3liiBF4+nu04Hp486T8lzXc2o6eWh8yIIHKRpxqo95vRZTyiW5JcRIAiSQQAK6j+f0gBPYGfgo7xCAqCL/sbFei++nuk8pETZfRdlFSmjNYovzbHMyc7wDgy0hAZcSoAA7aJiWlhbZtGlTsISjR49KaWmp4P2Ghgb1vt17c1VJd4M5iNu3j+4c7gw7KP6Xfb8M44HpYvNB8dhnu6psFVM0+rbnsOGpJqD7eO5aD/jYsWNSX18vra2tysYQ3ebmZtmzZ4/gd21tbbJly5aw9/B7CnCq/yzcWT6O1IOwqunjq5HI+N410mWpMFI0rq5YbUlkAQ+3urDanQ1jrUjApwQowA4afseOHdLY2KhK2L9/v9TW1srOnTuV+O7du1e2bt0a9t6RI0cowA7aRJdHX5m6Ynuk3sjkiKUJORk5cnP1zXJrdeA4PQgtjtQrzC7UpamsJwn4lgAF2EHTmwUY3vCuXbtk+/btMjAwoF5DkEPfMwQ7UrV0N5iDuLV9dP9Yv/JozV7t6Z7TgiQX5qs0tzTMq11buVayMrK0bTsrTgJ+JqD7eO7aKWh0qkgCDA8YgmwWYOM9CrB3/xyx5adjqCNsvfa1gdfCGo10jKHrtStKVnDLj3e7B1vmQwIUYAeNnogp6AMHDgi+jKu/v186OzsdrDUfnQgC8F6RuCJ0vbZ3rNfyeJxRCy/WfPAAppEr8isSUQ0+gwRIwMUEKMAOGscswAjCQvAVop/NQVih7zEIy0GDOPRoHJ2HlIzmRBb4eWxqzFJiXmaeWp+F2BrrtThSjykaHTIMH0sCLidAAXbIQNhehEArTDVDVLH96NChQ2ra2bwNye69uaqku8Ecwp20x/aM9oQdFH+296wgQtl8VeRVBFM0Gt7tDRU3MEVj0izFgkjA/QR0H89dvQbshPl1N5gTTJx4JtZrsTZrPngAr7GGG3qtKl1lWa/FFPKyomVcr3XCMHwmCXiIgO7jOQXYQ50xVU2ZnJ4URB2HrtcOjg9aqoQD4ZElyrxeu2HxBkF0Mi8SIAESiJcABTheYim+X3eDpRifXJ64rFI0mtdrX7r0UliKRuyjRf5j83otjtRjisZUW5Dlk4B3COg+ntMD9k5fTHhLLl6+aDkoHh4uMkmFHqm3uHDxtbNrFweSWdSV1zFFY8ItwgeSAAmYCVCANesPuhssGbg//tjH5Ye//KFAgEOvgqwCJazbarbJZ+74jBJbCDAvEiABEkg2Ad3Hc3rAye4xGpT3tn94m/z0tZ+GebpG1cvzyuX7H/y+YAtQUU6RBi1iFUmABLxIgAKsmVV1N1gycHePdKtcyI+eeVS++vOvqmP4Il2IVl5XtU7WVqwNfK9cK+sq1ymvOC0tLRnVZRkkQAI+JaD7eE4P2KcdN9ZmY3/u985+Twlx84Vm9bE7l98p1QXVcqbnjLza92pYzmXcU5JTosTYEGT1vWqd1JbVCqKheZEACZDAQglQgBdKMMmf191gScYVLA77ehvbGpUQYx34h7/1Q/U7bEFq629TYoytSMb3092nZXhiOKy6WelZ6qg/izBXrpM1lWt4AlGqjMtySUBTArqP5/SANe14qaz2s+efla1Lt845xQzB7rzcGRDk7mvCjJ/PD5+3rX5NcY2avg71muFtczo7lRZn2STgTgIUYHfaJWKtdDeYZrhtqzs0PiRne85aPGYIM7Y4hR4hiAcgUUfYdHblOllVtorT2V7oEGwDCcyTgO7jOT3geRqeH0s8AUxnt/a3Bj3mM73XvGe76ezsjGxZXR4ynV21TtZUrJGC7ILEV5BPJAEScBUBCrCrzBG9MrobLHoLvXcHprMvDF8IW2eG14z37S6c/WueyjZeLypYxOls73URtsinBHQfz+kB+7TjeqXZg1cGBacpha4zYzp7enY6rJlluWW20dk4ECIjPcMrWNgOEvAFAQqwZmbW3WCa4U5ZdSemJ6S1rzVsnRleM/JZh16YzsZxh6HR2XiP09kpMyMLJoE5Ceg+ntMDZgf3FQFMZyMK2/CYzdunELVtd60sWRkWmQ2hrsqv4nS2r3oPG+s2AhRgt1kkSn10N5hmuLWq7sCVAdvobHjSdtPZSMlpt858Xel1nM7WyvKsrK4EdB/P6QHr2vNY76QRGJ8at0RnGwlH4D2PTI6E1SMnIyd8OrtqnXovPys/afVmQSTgdQIUYM0srLvBNMPt6eoiTef5ofPX1pmRcOTq1qmuka6wtqdJmqwsDZnOvpp4pKqgytOs2DgScIKA7uM5PWAnegWf6XsC/WP9attUaIpOpO2EcIdeFXkVtuvMWH9mdLbvuxMBRCBAAdasa+huMM1ws7ohBDCdjQMsQsUZP49Ojobxys3MjRidnZeVR74k4GsCuo/n9IB93X3ZeLcQgFfcMdQRtp8Z682XRi7ZTmcj2Cs0bzZ+rsyvdEuzWA8ScJQABdhRvIl/uO4GSzwRPtHtBPrG+q55zKZ15nMD52ynsyHAofuZ8TPWn9PT0t3eXNaPBGImoPt4Tg84ZlPzRhJwF4ErU1fk1d7AdLb5KEgcdDE2NWY7nY082TiXeW1F4HxmCDOiszHVzYsEdCNAAdbMYrobTDPcrG4KCGA6+43BN2y95u7RbtvpbJwsZec1V+RXpKAFLJIEYiOg+3iujQfc1tYmdXV1UltbqyxTWloqJ06ckPr6etm3b58cPnxYdu/erd6f69LdYLF1S95FAvYEekd7baOzz/Wfk1mZDfsQsn3ZrTPjsAtOZ7OXpZqA7uO5NgJ87Ngx2b59u7I3xBhf+LmsrEy9BxHes2dP1P6gu8GiNpA3kMA8CGA6+5XeVwLT2aZ1Zhx0gd+FXnmZebKmco2su7qP2fi+umI1p7PnwZ8fmR8B3cdzbQTYbJ5Dhw4Fxfbo0aOyc+fOmK2nu8FibihvJIEEEMB09usDr9seBdkz2hNWArxinCwVus4MLxqpO3mRQCIJ6D6eayfA8HwHBgZk06ZNyo779+9Xr1taWqShoSGqbXU3WNQG8gYSSBIBCHBwP7PJa35t4DXb6WycxWy3zry8ZDmns5NkM68Vo/t4rp0A7927V003h17wiiHMoSJ84MABwZdx9ff3S2en/ak3XuucbA8JpILA2OTYtensntNB7xnR2ePT42FVQn5su+js1eWrJSczJxVNYJmaEKAAJ9lQmzdvVsFXuDD93NfXp6ajsUZ85MgROXjw4Jw10t1gScbN4kggYQSmZ6bl9cGr09nwmK9un8IWKux1Dr0wnV1bVhu2zozp7dLcuYMtE1ZpPsjVBHQfz7XygOHh3n333UEBxrQzoqIR+QwPGFe0QCzdDebqvwZWjgTmSaB7pNt2nRnrz3bR2dUF1bbrzMuLl/OM5nnaQMeP6T6eayfA9913n/J0jQvCW15eLk1NTbZT06GdSneD6fhHwjqTwHwJID+2XXQ23rObzi7IKogYnZ2dkT3favBzLiWg+3iulQAnog/obrBEMOAzSEB3ApjORrBXaBYwbKHqv9If1ryMtIzAdHZIFjBsnyrJLdEdh2/rr/t4TgH2bddlw0nAewRmZ2cF2b6C+5lNaTqx/mx3LS5cbLvOvKxoGaezXd5FKMAuNxCnoDUzEKtLAg4RGJkYUdPZ5rzZEGq8NzE9EVZqYXah2jYVunWqrrxOOJ3tkJHifCwFOE5gqb5dd4Olmh/LJwGvEcB0Nk6WCvWaIdQDVwZsp7OvL7/eIsyGUHM6O7m9Q/fxnFPQye0vLI0ESEATApjOxlnMoevM+BmHXdhdSwqX2K4zLy1ayulsB+xOAXYAqpOP1N1gTrLhs0mABGIjcHnicmA627Sf2ZjOnpyZDHtIUXZR+HR21TqpK6uTrIys2ArlXWEEdB/P6QGzU5MACZBAgghMzUwJTpYy1pmDyUa6T8vg+GBYKZnpmRI6nY3IbBx0UZxTnKBaefcxFGDNbKu7wTTDzeqSAAmICKazu0a6bKOz24fabRkhCtscAKZeV60TTHOnpaWRq4joPp7TA2Y3JgESIIEUEsB0NvJkm6OzX+p6SV7pe8W2VvCMjaCvB9/yoKwsXZnC2qe2aApwavnHXbruBou7wfwACZCANgQwhd0x1CFt/W3yV8/+lTx65tE56/7Gp98QnCbl10v38ZwesF97LttNAiSQdAKYisbBE9j2BJHFF9aM2wYCrxFdDRGOdOVm5spv3/Lbcu+Ge2VyelLect1bfD0dTQFOehdeWIG6G2xhreenSYAEnCYwPjWuTn0yBFaJrElwh8aHbKuQJmlSU1yjUmbiC8c0HmgKHKWKs5Q/tfVT8omtn5DK/Eqnm6DN83Ufz+kBa9PVWFESIAE3EDACqiIJ7Pmh87YnOKHuJTklQYFdVboq+BqCu6JkheX8428d/5Z8p/k7cv8d98u9t9wr8H55WQlQgDXrEbobTDPcrC4JaEkAaSvNXqvZi8WU8djUmG27sK1oZcnKiCJbllcWM49Xe18VpL3Euci87AnoPp7TA2bPJgES8B0BpJ9EsJNlLdY0TYwMWJGuqvwqi+dq9mSXFS8TiDCv5BCgACeHc8JK0d1gCQPBB5GAxwn0j/VH9GJfH3hd7DJWAQmmes2ianldtkpwSAMvdxDQfTynB+yOfsRakAAJxEkAJxghajjSWqzdQQpGEUhyYQQ74btZZKsLqzntG6ctUnU7BThV5OdZru4Gm2ez+TES0I6AcRhCpLVYTCHPzM7Ytgu5lyMJLBJXMKBJu+5gW2Hdx3N6wN7oh2wFCWhJYHRyVO2DjbQWi9/bXRlpGSpqOJLIlueV+3p/rJadYR6VpgDPA1oqP6K7wVLJjmWTQLwE4KFiW04kL/bi5YsRH4n9rpHWYpH9icFO8VrDe/frPp7TA/Zen2SLSCCpBAavDEYU2NcGXhOs1dpdORk5cl3pdbZe7KqyVTwNKKlW1LMwCrBmdtPdYJrhZnU9QAApDxHsFMmLRWrFSBdO7ok0TbykaAmDnTzQP1LZBN3Hc3rAqew9LJsEXEAAwU49oz2BvMTmHMVXX0N8IwU7FWQVWATWHFEM7zYvK88FLWQVvEqAAqyZZXU3mGa4WV2XEBibHBNMB0fyYnEknt2FLEwIdoq0Fot1Wp5N6xIj+7Aauo/n9IB92GnZZO8RgIfaOdwZUWAvDF+I2GhEDIfmJTZ+hvhmZWR5Dxhb5AkCFOAkmrG+vl727dsnhw8flt27d0tpaakcPXpUfW9paZGGhoaotdHdYFEbyBs8S2B4fDhsmtiYNsZWnvHpcdu2Z6VnzRnsVJpb6llmbJi3Ceg+nmvlAZeVBRKZQ4T37NmjRLe5uVm9PnbsmLS1tanXc126G8zbf07+bh3OgW0fbI+4Fot12kjX4sLFEb3YpUVLJSM9w99w2XpPEtB9PNdKgOHt7ty5M9iR9u/fL7W1teo9iO/evXvlyJEjFGBP/qnp3yjjMPbQYCfjZ+Qnnp6dtm0ozoYNTZloTBMj2Kkgu0B/QGwBCcRJgAIcJ7CF3A7B3bRpU3C6GVPSu3btku3bt8vAwIB63djYSAFeCGR+dkEEcBj7XMFOcx3GjuQSkdZicSA7g50WZBp+2IMEKMApMOqhQ4eU4La2tgYFGB4wBDlUgA8cOCD4Mq7+/n7p7OxMQa1ZpBcIwItF9qZIXmysh7Ebe2PNwU45mTleQMQ2kEDSCFCAk4Qa0899fX3B9V5MNdfV1XEKOkn8/VQMtuQgqMluX2y0w9gxHRzJi43nMHY/8WZbSWC+BCjA8yUX5+cQcIX1XkQ8wwPGtWXLFhV8hehnBmHFCdTHtxuHsUfyYuc6jB1TwZEEtqa4hsFOPu5XbHryCVCAk8gcwlteXi5NTU0qEhoX3oMwcxtSEg2hQVE4jH2uYKe5DmMPDXYyfkZ+Yh7GroHxWUXfEKAAa2Zq3Q2mGW7HqosE/4gaNkS2ta9VXu5+Wc72npXOy50S6Ri7NEmTZcXLInqx2M7DYCfHzMYHk0BCCeg+nmu1DSkRltPdYIlgoPMzvvSTL8k/PPcP0j7ULrMya9uUvMw8WVO5xlZkeRi7ztZn3UnASkD38ZwCzB6tFYH679fLwy0PRxRfHNT+5pVvFmzpqSmqEazLqtfFgdcVeRX0cLWyOCtLApEJUIA16x26G0wz3AmvLiKUczNzpetyl3yn+TtyqOWQzBU0FVoBnEFriLHxfXnxNYHGe1UFVTwmL+GW4wNJIPEEdB/P6QEnvk/wiUkkgMQXR14+It9+9tty/PxxuWnRTfLo+x+VjqEOyxemrI33uka65qwhcidjnTgo1EVWLxrvVxdUM+I5iXZmUSRgR4ACrFm/0N1gmuFOanWbzjfJt5u+LV940xdkdcXqiGUjgAunA4WKNH42hBonC0VaY8aDMdWNHMtmbzrUk8aB85npmUllwMJIwE8EdB/P6QH7qbeyrTETmJyeVBmvwkR6uEMdmID3IeKRcjejIJyli6jquTxpiHh2RnbM9eKNJEAC1whQgDXrDbobTDPcnq4uEnpgOnsuTxqpKSPtOTbgYDp7Lk8a0+FY9+ZFAiRgJaD7eE4PmD2aBBwkMDM7I90j3baetJryvupNRzrL16haZX7lnJ70sqJlPBHJQTvy0e4kQAF2p10i1kp3g2mGm9WNgQAOeOgd67X1pA3vGmvTkZKLGEWU5ZbN6UnDyy7KKYqhRryFBPQgoPt4Tg9Yj37GWvqcAER64MqAvUgPByK+4U0PTwzPSao4p3hOTxoiXZJTwr3SPu9vujSfAqyLpa7WU3eDaYab1U0yAZw3bLcmbfakIeRzXQVZBUGRNic0Ma9Tl+eVU6STbFsWF05A9/GcHjB7NQn4jACSmSA4zC7C23ivZ7RnTioICouUyMR4H+vWiATnRQJOEaAAO0XWoefqbjCHsPCxJGAhMDY5JueHbUTalOAkWkITbK9CcFgwHejV1KBmT7q6sJoizb43bwK6j+f0gOdten6QBPxNAFnIbBOaXF2ThjcdLaEJEpUYCU1CE5kYQo291Exo4u++Fqn1FGDN+oXuBtMMN6vrcwJIaILjIedal44locmSwiVzetIQ8ayMLJ/T9l/zdR/P6QH7r8+yxSTgKgJIaBIp65gh3NESmuCcZ0xnq+nukMM1DE8a0+E5mTmuajsrszACFOCF8Uv6p3U3WNKBsUAScAEBI6GJ+VANO686WkKTqvyqa8dT2qxJI+tYfla+C1rMKsRCQPfxnB5wLFbmPSRAAq4ngL3SiN6eK7obe6XHpsbmbAu2WGHNGXuma8tq5ZbqW+TGyhst50oXZhe6nocfKkgB1szKuhtMM9ysLgm4ioCR0GQuTxq/w1atuS6I89679sofvemPXNU+v1VG9/GcHrDfeizbSwI+JYC1ZqT87LrcpdacsY1KfcfPI4HvxnvI3x3pOEp4x/hquLNBvvzWL/uUpjuaTQF2hx1iroXuBou5obyRBHxAAGvDfWN914T0qrCaRdUQ20sjlwT3R7vyMvPUMZL4nHm6+p66e6R+c728+4Z3q4hrPIuJRqLRdPb3uo/n9ICd7R98OgmQQJwEME0MUTV7qMHXZq/18kWBqM51JrNRdE5GjhJVREqr7wXXvoe+h/XdwfFBWfr1perwio/d+jG5b/N9aj2Yl7sIUIDdZY+otdHdYFEbyBtIwIUEjLVXy7SveRp4pCs4NQxRjXaGMpqITFtBIS2stohqqNhizTYtLS1mMs90PKMOt3jP2veocni5k4Du47lWHvChQ4dULzhx4oQcPHhQva6vr5d9+/bJ4cOHZffu3VJaWjpnT9HdYO78M2Ct/EgAoorDH8LWU41pYJO3CuGdmJ6Iiglrq2ZRXVxg9VrN3mppbmlcohq1cN6gHQHdx3NtBPjYsWNSW1urvvbv3686SkNDg5SVlanXEOE9e/ZE7UC6GyxqA3kDCSyAAEQVRxqaA5XCgpZMwhpt3y2qkpGWoaZ+Y/FWcaZxPJ7qAprKj3qAgO7juTYCDO93YGBAie7Ro0elsbFRecF4vXPnzpi7ku4Gi7mhvJEETASwrcaI+I02DRxtnywei+CjRQWLYvJWsa+WwUrsjk4Q0H0810aAzcbDtPOOHTuU8MIb3rRpk7S0tChxjnbpbrBo7ePv50ngxRdFnnlG5L775vmA5H9sZGJEBSrZbasJDVoanRyNWkGkc6wqqIoapARvtiKvQjLSM6I+U9sbJidFvvpVkfvvFyku1rYZXq+47uO5dgLc1tamppuNNWCjg5k9ZHOnO3DggODLuPr7+6Wzs9Pr/ZLti5fArl0iR4+KfO1rIp/9bLyfTtj9OAbQzkO121YTLVkEKgVRxbm8lilgIwLYFBGM3+M+njokIuPjIh/4gMijj4r89V+LfOITCbMvH5RYAhTgxPKM+rS9e/cqAcaF6ee+vj619os14iNHjoQJc+gDdTdYVEC8YX4ERkZE3vtekcZGkS9/WeSLXxSJI2p2rkKvTF0JJnmYaz0VwougplgueKCxbKmBR0tRjYXo1XtGR0Xe9z6RH/1I5EtfEvmTP0lYP4ijFrw1RgK6j+daecAQ3O3bt6tIZwhueXm5CsrCz0aEdLRALN0NFmO/5G3zIXDlisju3SLf/77I5z6HyL6Igy+mdH/++s/l+PnjcuviWwWZkyzRwKapYewpjeVCAJKdqIZuqcHaK4/ei4VonPcMDYm85z0iP/2pyEMPiXz+83E+gLcnm4Du47k2AgzB3bVrlxJdeL0PPPCAWvOF8OK9pqamoGc8VyfQ3WDJ7uC+Kw9rfx/6kMi//Zv0fPDX5bv3bZaL4z3y2sBr0trfqg6Yh6BOzUzFhKYkpyQsAYRlm83VpBAQVR6VFxNSZ266dEnkne8UaWkR+eY3RX7/950ph09NKAHdx3NtBDhRVtPdYIniwOfMQWB6WuRTnxL5znfkP9eKfPB9IuMRznqHcMJzxZQwpntx5mxNSY2sKl0l15ddLzjeriK/gskc3NzhXn9dZMcOkbY2kb//e5F773VzbVk3EwHdx3MKMLszCdgRmJ2V7s9+Uqr+4m/k9ZuXy+/Xr5Qf9j4jU7Oxeb6hjyzKLlJBThBj9T0v5Hvo+/kVkpuZS9s4TeDkSZF3vUtkYCAQhPcrv+J0iXx+AglQgBMIMxmP0t1gyWDEMkwEEAX7e78ncsMN0vcf/yz/d/B/5JGTj8hLl15SwU0/+8jPxDhlB2fR9o72qjNpceqO5ftor8pvHOmEHTvmBVkFUQXbIub5FTxMPp7O++MfiyCHQE6OyPe+J3LnnfF8mve6gIDu4zk9YBd0IlbB5QQee0zkgx8M7Ad97DGZve02ab7QrIT49prb5SO3fiSmBkCoB64MhAmzWahDRRtCHssJPkYFcJJPPF42BBxC77vsUw8/LPK7vyuyYoXIf/2X+geLl34EKMCa2Ux3g2mG2zvVPX5c5Nd+TQSRsn/3d4F9oiKC1I1OihfEd/DK4JyiHfS2Td53rABifRwAAA5hSURBVEFiaAMOGwj1pCvzKu2nzK9OlWNK3cl2O9ZxsL6PCPe//EuRbdsCnu+iRY4Vxwc7S0D38ZwesLP9g0/3EgEE67z73SLImoU9on/8xyLp6a5roXFIgp0wB6fJx3oCom4S7VgOSzAam5WepTzt4Fo2hDnPZo3btLZdkluS2pSUg4Miv/VbIj/4QWBG42//ViSX6+yu68BxVIgCHAcsN9yqu8HcwNDXdRgeDgzi2CsMj/gf/9ETqQoh2sisFU20Q9e4kWQk1guHMiAvtDkYLZpol+WVRRTtuGYfTp8W+fVfF3nlFZEHHxT5wheYYCNWw7n4Pt3Hc3rALu5crJpLCczMBDIk/emfiqxZI/Kf/ymybp1LK+tstZCQxOxJh65h94yFB6bFkpfaqDUOccA2L7sIcvyzgGQot1TfIrctu03evPLNUldWp0TekqcaKSU//OHAI/EPE5Jt8PIEAQqwZmbU3WCa4fZ2dSG8GNhnZ0VwVvVv/qa325ug1iHfNcQzLGJ8jghyHJEYz4XAsiW5lbL/vzPkvT9oCwRZQYh9+o9SPOx0ulf38ZwesE69jXV1H4GzZ0VwkAPWhX/nd0T+4i9E8vLcV08X1QgBYggsQ0Q4vvqv9Adfq5/Hrv48fu01tnCp98cHJJa16ppBkX85KvK/2iVgn+9+1xNLBS4yoyuqQgF2hRlir4TuBou9pbwzaQSQwB+ZsxAdfdNNIv/6ryLr1yet+GQXhMjs4fFhq2iaRNQsoBZBvSq48XqzRvuQmKQ0t1RNSeM7vs4NnJMzPWeCCDYv2SwP9t0q79j/75J2eUTG/s9XJP8PP8f13mR3kiSVp/t4Tg84SR2FxfiAwD//c+DoOuST/vM/D+wzTdCJSomkh+Clsamxa55mFC80VESRCzuevclG3ZG4xBBOs4iaXyuBzbsmsIbgIoI6NDMY2rHp0Ca5NHJJ7r35XvnI9Ttl3UMPi2CPL6ac/+VfRDZtSiQ6PstlBCjALjNItOrobrBo7ePvU0wA+YQRJf3MM4H8wo88IrJ8eUyVwtQq9uTGcuHeoKc511TueGCaN1REJ2cmYynGcg/OFoYQhnqhtj/biGh+Vn5C9w6j/c+ef1a2126XjCefEvnIRwL5nD/+8cCBCgUFcbeRH9CLgO7jOT1gvfoba6sDgakpkf37A5HSWA9G0oePftTiDSMrVvtguzxx7gn5+Rs/l6bzTbK6YrX86g2/GnFq1yyi8GDncxVmF9p6obGIaHFOcWr38do1GOc4Y082GCOhBoLhsD2Mly8IUIA1M7PuBtMMt7+r+/zzgSjpF16Q9s2rZc+7RZrz+2VofCimQCI7eMY6qN1UbqiIhk7l4mhET50j3NgoUl8vcu5cINAKebsrK/3d53zWet3Hc3rAPuuwbG6SCWA9+Gtfk8k/+ZJMzU7LV98k8ud3Rj7eMJLoIvUjPFCIbFFO4LXxnvl16O+Mn3EPvF/sq9X+6uwMpJPEmvvSpSIHDgSSbPDyHQEKsGYm191gmuFmda8SGHnppOT83qcl86f/I6MrlsrffHit/En+syr7FC6IIwKJhiaGlIeMKGP1feLq9/FhFTi10AvlGMI9XyHH5wuyC5Iv5vhn5tvfDqQAvXw5EPD2Z38mUlKyUCz8vKYEdB/P6QFr2vFYbQ0JIGHH4cMi998vcuGCTO+4Wx7/nXvkwJX/kR+3/lie/d/PysYlGyM2bHJ6UgmynTibxToo4BP2Qo7fj0+PLwggArIg5sr7jsMjt7s/anAWuOHQhIaGQCrJ228PeL2McF6QDb3wYQqwZlbU3WCa4WZ17QggnzSCtLBVaWJCBWj1fvaTcrE8W9YvSs7+YURRxyzkZk/c5JFDyPE1n4hqMxZMi8OrNk+XG176xnNj8rn/1y/Zv3haZNmygMd7772uPASDnT35BHQfz+kBJ7/PsEQSCBBobw8cCoD8xFlZgYCiz38+sK6p0TU+NR42XW7rkUcRcnj3xjGKGy+IPPgTkV99VWSmqFDSP9cg8pnPiOTna0SGVXWaAAXYacIJfr7uBkswDj7ODQROnQpspfmP/xDJzhb52McC062rVrmhdkmrAxJrTPzsv5WXm9P43zKTmyPnP/xeqX7w65Jdrdc/JUmD5vOCdB/P6QH7vAOz+S4i8NxzgSnWo0cDU6y/8Rsin/60yB13uDKjVsLIYd/0v/+7yDe+EUhgUlgYCLD6wz8UWbIkYcXwQd4jQAHWzKa6G0wz3KzufAicOSPy9a8HpqbHx0U2bw5MT+MQeYiTV66OjkCmMByUgNfV1YH0nZ/8pEhFhVdayXY4SED38ZwesIOdg48mgQUR6O4WOXgwkN0J68VFRYGEEx/6kMib36xnIBIOrkBEM/65+NGPRHC28tatAdH9wAdEcnIWhIwf9hcBCrBm9tbdYJrhZnUTQWB6OiBWEOIf/lAEU7Y1NSLve19gmvquu0QyMhJRkjPPwJ5d1B/nJ3//+yKIAi8tDQjunj0iGyNvvXKmQnyqVwjoPp57wgM+evSolJaWSktLizQgeGWOS3eDeeUPh+2YJ4GeHpEjRwIn/fziFyLYI4vpWhz88I53iNx9d0CcU3nBq335ZRGkioTw/uxngal0RDC/852Bwyre9S56u6m0kUfK1n08116AIbrNzc2yZ88eOXbsmLS1tanXkS7dDeaRvxs2IxEEurpEHntM5Ac/EHniCREcTIDruutE3vSmwNTuli0iGzY4u30H/xScOBH4On488I9BX1+gLgiievvbRd7zHpF77nG2HolgymdoRUD38Vx7Ad6/f7/U1tbKzp07lfju3btXjsBDiHDpbjCt/jpY2eQRQEKPJ58MeJs//3kgmhjrrbhwJvHKlSI33ihy/fUiK1YEfkbQEw4vKC8PHN0HDxXR18aFqW6IOqaQe3tFILQXLoi88YbI66+LnD0rcvq0yKVL1z5TVRUQf3y99a0it9zi7Qju5FmYJdkQ0H08116A6+vrZdeuXbJ9+3YZGBhQrxsx9UUB5h+snwlAPBFN3dysTmNSQokviCemrSNdWEuGYOMerD1HurBfua5OZN26wBfSQiJaG+KOz/MigSQQoAAnAfJcRZgFGB4wfjYL8IEDBwRfxtXf3y+dOE2FFwn4kQAONMCWHwgxoqzxhelieMv4gidtXLm5Aa8Y3jE8ZXi38JoN79nsLfuRJduccgIU4BSbgFPQKTYAiycBEiCBFBGgAKcIvFEsgrAQfIXoZwZhpdgYLJ4ESIAEkkiAApxE2JGKOnTokArE4jYkFxiDVSABEiCBJBGgACcJdKKK0d1gieLA55AACZCA7gR0H8+1j4KOtwPpbrB428v7SYAESMCrBHQfzynAXu2ZbBcJkAAJeJwABVgzA+tuMM1ws7okQAIk4BgB3cdzesCOdQ0+mARIgARIwEkCFGAn6TrwbN0N5gASPpIESIAEtCSg+3hOD1jLbsdKkwAJkAAJUIA16wO6G0wz3KwuCZAACThGQPfx3HcecHFxsdRcPS8VeaHLysoc6xy6PJgcApYiB3II/Ztln3B3n+jo6JChoSFdhtqwevpOgM0EdP/vKVG9jhwCJMmBHEL/ptgn2CcSNc7aPYcC/PLLTvLV4tkcZDjI8B9T+z9V/m3wb8PJQZwCTAGm53f1L4yDLQdbesD8R8RJwQ19tq8FGOcEf/KTn0wmb1eWRQ4Bs5ADOYT+gbJPsE84OWj7WoCdBMtnkwAJkAAJkMBcBCjA7B8kQAIkQAIkkAICFOAUQGeRJEACJOB2AjhffdOmTcFqHj16VEpLS2M6d93tbXNL/XwvwKGdrL6+Xvbt2yeHDx+W3bt3qw7nh4t/bH6wcvQ2cpANMPLrOGD0kGPHjikGra2t6i2MD83NzbJnzx7B79ra2tRrXgsj4GsBDu1kQGkk5oAI+6WD8Y8t/I/IjwMwB9lr/cCP40DoX8GOHTuksbFRvb1//36pra2VnTt3KvHdu3evHDlyZGHqw0+LrwUY9jd3MvwMDwCdzG8X/9isFvfjAMxB9lof8Os4YP4rMI8J+Id0165dsn37dhkYGFCvDXH221iZyPZSgE3/5Rn/6WHdA95AQ0NDIlm7+ln8Y7Oax48DMAfZa30A/4z4cRyIRYDhAaOvUIAXPqRTgEME2EB66NAh9Z+eX0Q4kgD79Y/NjwOwWYD9avfQIdVv40AkAebsyMLF1u4JnhZg/PHYXea1XbPwwOvp6+sLBhpgjePgwYPOkE/iU+PlwD+2a8bx0wBMuwfs7tVxIN4hxzw2YkYQsSJwSBiEFS/JyPd7WoBjwRTayRBogMhnQ7T8Eojltz+2uf4p8esAzEE2MGKAg1/HAWPMxN8AAq0wK4Ix0BgTwcVvy3Ox6Mh87/G1AEfqZOXl5dLU1KS2I/nh4h+b1cp+HoDxjwkHWVH/gPttHPDDWOe2NvpagN1mDNbHPQQ4ALvHFqwJCXiVAAXYq5Zlu0iABEiABFxNgALsavOwciRAAiRAAl4lQAH2qmXZLhIgARIgAVcToAC72jysHAmQAAmQgFcJUIC9alm2iwRIgARIwNUEKMCuNg8rRwIkQAIk4FUCFGCvWpbtIgESIAEScDUBCrCrzcPKkQAJkAAJeJUABdirlmW7SIAESIAEXE2AAuxq87ByJEACJEACXiVAAfaqZdkuEiABEiABVxOgALvaPKwcCZAACZCAVwlQgL1qWbaLBEiABEjA1QQowK42DytHAiRAAiTgVQIUYK9alu0iARIgARJwNQEKsKvNw8qRAAmQAAl4lQAF2KuWZbtIgARIgARcTYAC7GrzsHIkQAIkQAJeJUAB9qpl2S4SIAESIAFXE6AAu9o8rBwJkAAJkIBXCVCAvWpZtosESIAESMDVBCjArjYPK0cCJEACJOBVAhRgr1qW7SIBEiABEnA1gf8PDNgIqHQxuGEAAAAASUVORK5CYII=\" width=\"640\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ca0ee1e13d704d76840a135522769a25",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.5, description='Initial point', max=4.0, min=-4.0), FloatSlider(valu…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = lambda x: np.power(x, 2)\n",
    "gradf = lambda x: 2 * x\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "\n",
    "\n",
    "def update(x0, step):\n",
    "    gd = methods.fo.GradientDescent(f, gradf, ss.ConstantStepSize(step))\n",
    "    _ = gd.solve(np.array([x0]), max_iter=10)\n",
    "    x_hist = gd.get_convergence()\n",
    "    x = np.linspace(-5, 5)\n",
    "    ax.clear()\n",
    "    ax.plot(x, f(x), color=\"r\", label=\"$f(x) = x^2$\")\n",
    "    y_hist = np.array([f(x) for x in x_hist])\n",
    "    x_hist = np.array(x_hist)\n",
    "    plt.quiver(x_hist[:-1], y_hist[:-1], x_hist[1:]-x_hist[:-1], y_hist[1:]-y_hist[:-1], \n",
    "               scale_units='xy', angles='xy', scale=1, width=0.005, color=\"green\", label=\"Descent path\")\n",
    "    ax.legend()\n",
    "    fig.canvas.draw()\n",
    "\n",
    "step_slider = ipywidg.FloatSlider(value=0.8, min=0, max=1.2, step=0.1, description=\"Step\")\n",
    "x0_slider = ipywidg.FloatSlider(value=1.5, min=-4, max=4, step=0.1, description=\"Initial point\")\n",
    "_ = ipywidg.interact(update, x0=x0_slider, step=step_slider)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def plot_alpha(f, grad, x, h, alphas, beta1, beta2):\n",
    "    df = np.zeros_like(alphas)\n",
    "    for i, alpha in enumerate(alphas):\n",
    "        df[i] = f(x + alpha * h)\n",
    "    upper_bound = f(x) + beta1 * alphas * grad(x) * h\n",
    "    lower_bound = f(x) + beta2 * alphas * grad(x) * h\n",
    "    plt.plot(alphas, df, label=r\"$f(x + \\alpha h)$\")\n",
    "    plt.plot(alphas, upper_bound, label=\"Upper bound\")\n",
    "    plt.plot(alphas, lower_bound, label=\"Lower bound\")\n",
    "    plt.xlabel(r\"$\\alpha$\", fontsize=18)\n",
    "    plt.legend(loc=\"best\", fontsize=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: x**2\n",
    "grad = lambda x: 2 * x\n",
    "beta1 = 0.1\n",
    "beta2 = 0.9\n",
    "x0 = 0.5\n",
    "plot_alpha(f, grad, x0, -grad(x0), np.linspace(1e-3, 1.01, 10), beta1, beta2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## $f(x) = x\\log x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x_range = np.linspace(1e-10, 4)\n",
    "plt.plot(x_range, x_range * np.log(x_range))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0 = 1\n",
    "f = lambda x: x * np.log(x)\n",
    "grad = lambda x: np.log(x) + 1\n",
    "beta1 = 0.3\n",
    "beta2 = 0.7\n",
    "plot_alpha(f, grad, x0, -grad(x0), np.linspace(1e-3, 0.9, 10), beta1, beta2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Backtracking \n",
    "\n",
    "```python\n",
    "def SelectStepSize(x, f, h, rho, alpha0, beta1, beta2):\n",
    "    \n",
    "    # 0 < rho < 1\n",
    "    \n",
    "    # alpha0 - initial guess of step size\n",
    "    \n",
    "    # beta1 and beta2 - constants from conditions\n",
    "    \n",
    "    alpha = alpha0\n",
    "    \n",
    "    # Check violating sufficient decrease and curvature conditions\n",
    "    \n",
    "    while (f(x - alpha * h) >= f(x) + beta1 * alpha grad_f(x_k).dot(h)) and \n",
    "    \n",
    "          (grad_f(x - alpha * h).dot(h) <= beta2 * grad_f(x_k).dot(h)):\n",
    "            \n",
    "        alpha *= rho\n",
    "        \n",
    "    return alpha\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Теоремы сходимости (Б.Т. Поляк Введение в оптимизацию, гл. 1, $\\S$ 4;  гл. 3, $\\S$ 1; Ю.Е. Нестеров Введение в выпуклую оптимизацию, $\\S$ 2.2)\n",
    "От общего к частному:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Теорема 1.** \n",
    "Пусть \n",
    "\n",
    "- $f(x)$ дифференцируема на $\\mathbb{R}^n$, \n",
    "- градиент $f(x)$ удовлетворяет условию Липшица с константой $L$\n",
    "- $f(x)$ ограничена снизу\n",
    "- $\\alpha = const$ и $0 < \\alpha < \\frac{2}{L}$\n",
    "\n",
    "Тогда для градиентного метода выполнено:\n",
    "\n",
    "$$\n",
    "\\lim\\limits_{k \\to \\infty} f'(x_k) = 0,\n",
    "$$\n",
    "\n",
    "а функция монотонно убывает $f(x_{k+1}) < f(x_k)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Теорема 2.** Пусть\n",
    "- $f(x)$ дифференцируема на $\\mathbb{R}^n$\n",
    "- $f(x)$ выпукла \n",
    "- $f'(x)$ удовлетворяет условию Липшица с константой $L$\n",
    "- $\\alpha = \\dfrac{1}{L}$\n",
    "\n",
    "Тогда \n",
    "\n",
    "$$\n",
    "f(x_k) - f^* \\leq \\dfrac{2L \\| x_0 - x^*\\|^2_2}{k+4}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Теорема 3.**\n",
    "Пусть \n",
    "- $f(x)$ дважды дифференцируема и $\\mu\\mathbf{I} \\preceq f''(x) \\preceq L\\mathbf{I}$ для всех $x$\n",
    "- $\\alpha = const$ и $0 < \\alpha < \\frac{2}{L}$\n",
    "\n",
    "Тогда \n",
    "\n",
    "$$\n",
    "\\| x_k - x^*\\|_2 \\leq \\|x_0 - x^*\\|_2 q^k, \\qquad q = \\max(|1 - \\alpha l|, |1 - \\alpha L|) < 1\n",
    "$$\n",
    "\n",
    "и минимальное $q^* = \\dfrac{L - \\mu}{L + \\mu}$ при $\\alpha^* = \\dfrac{2}{L + \\mu}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### От чего зависит $q^*$ и как это использовать?\n",
    "Из Теоремы 3 имеем \n",
    "\n",
    "$$\n",
    "q^* = \\dfrac{L - \\mu}{L + \\mu} = \\dfrac{L/\\mu - 1}{L/\\mu + 1} = \\dfrac{M - 1}{M + 1},\n",
    "$$\n",
    "\n",
    "где $M$ - оценка числа обусловленности $f''(x)$.\n",
    "\n",
    "**Вопрос**: что такое число обусловленности матрицы?\n",
    "\n",
    "- При $M \\gg 1$, $q^* \\to 1 \\Rightarrow$ оооочень **медленная** сходимости градиентного метода. Например при $M = 100$: $q^* \\approx 0.98 $\n",
    "- При $M \\simeq 1$, $q^* \\to 0 \\Rightarrow$ **ускорение** сходимости градиентного метода. Например при $M = 4$: $q^* = 0.6 $\n",
    "\n",
    "**Вопрос**: какая геометрия у этого требования?\n",
    "\n",
    "**Мораль**: необходимо сделать оценку $M$ как можно ближе к 1!\n",
    "\n",
    "О том, как это сделать, Вам будет предложено подумать в домашнем задании :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Вычислительный аспект и эксперименты\n",
    "1. Для каждого шага метода нужно хранить только текущую точку и вектор градиента: $O(n)$ памяти\n",
    "2. Поиск $\\alpha_k$:\n",
    "\n",
    "    - дан априори\n",
    "    - ищется из аналитического решения задачи наискорейшего спуска\n",
    "    - заканчивается за конечное число шагов\n",
    "3. Для каждого шага метода нужно вычислять линейную комбинацию векторов: $O(n)$ вычислений + высокопроизводительные реализации"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Pеализация градиентного спуска"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def GradientDescent(f, gradf, x0, epsilon, num_iter, line_search, \n",
    "                    disp=False, callback=None, **kwargs):\n",
    "    x = x0.copy()\n",
    "    iteration = 0\n",
    "    opt_arg = {\"f\": f, \"grad_f\": gradf}\n",
    "    for key in kwargs:\n",
    "        opt_arg[key] = kwargs[key]\n",
    "    while True:\n",
    "        gradient = gradf(x)\n",
    "        alpha = line_search(x, -gradient, **opt_arg)\n",
    "        x = x - alpha * gradient\n",
    "        if callback is not None:\n",
    "            callback(x)\n",
    "        iteration += 1\n",
    "        if disp:\n",
    "            print(\"Current function val =\", f(x))\n",
    "            print(\"Current gradient norm = \", np.linalg.norm(gradf(x)))\n",
    "        if np.linalg.norm(gradf(x)) < epsilon:\n",
    "            break\n",
    "        if iteration >= num_iter:\n",
    "            break\n",
    "    res = {\"x\": x, \"num_iter\": iteration, \"tol\": np.linalg.norm(gradf(x))}\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Выбор шага\n",
    "\n",
    "Реализации различных способов выбора шага приведены [тут](https://github.com/amkatrutsa/liboptpy/blob/master/step_size.py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Зависимость от обусловленности матрицы $f''(x)$\n",
    "Рассмотрим задачу \n",
    "$$\n",
    "\\min f(x),\n",
    "$$ \n",
    "где\n",
    "$$ f(x) = x^{\\top}Ax, \\; A = \\begin{bmatrix} 1 & 0\\\\ 0 & \\gamma \\end{bmatrix} $$\n",
    "\n",
    "$$\n",
    "f'(x) = 2Ax\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def my_f(x, A):\n",
    "    return 0.5 * x.dot(A.dot(x))\n",
    "\n",
    "def my_gradf(x, A):\n",
    "    return A.dot(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.rc(\"text\", usetex=True)\n",
    "\n",
    "gammas = [0.1, 0.5, 1, 2, 3, 4, 5, 10, 20, 50, 100, 1000, 5000, 10000]\n",
    "# gammas = [1]\n",
    "num_iter_converg = []\n",
    "for g in gammas:\n",
    "    A = np.array([[1, 0], \n",
    "                  [0, g]], dtype=np.float64)\n",
    "    f = lambda x: my_f(x, A)\n",
    "    gradf = lambda x: my_gradf(x, A)\n",
    "#     x0 = np.random.rand(A.shape[0])\n",
    "#     x0 = np.sort(x0)\n",
    "#     x0 = x0[::-1]\n",
    "    x0 = np.array([g, 1], dtype=np.float64)\n",
    "#     print x0[1] / x0[0]\n",
    "    gd = methods.fo.GradientDescent(f, gradf, ss.ExactLineSearch4Quad(A))\n",
    "    x = gd.solve(x0, tol=1e-7, max_iter=100)\n",
    "    num_iter_converg.append(len(gd.get_convergence()))\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.loglog(gammas, num_iter_converg)\n",
    "plt.xticks(fontsize = 20)\n",
    "plt.yticks(fontsize = 20)\n",
    "plt.xlabel(r\"$\\gamma$\", fontsize=20)\n",
    "plt.ylabel(r\"Number of iterations with $\\varepsilon = 10^{-7}$\", fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- При неудачном начальном приближении сходимость для плохо обусловенной задачи очень медленная\n",
    "- При случайном начальном приближении сходимость может быть гораздо быстрее теоретических оценок"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Эксперимент на многомерной задаче\n",
    "Пусть $A \\in \\mathbb{R}^{m \\times n}$. Рассмотрим систему линейных неравенств: $Ax \\leq 1$ при условии $|x_i| \\leq 1$ для всех $i$.\n",
    "\n",
    "**Определение.** Аналитическим центром системы неравенств $Ax \\leq 1$ при условии $|x_i| \\leq 1$ является решение задачи\n",
    "$$\n",
    "f(x) = - \\sum_{i=1}^m \\log(1 - a_i^{\\top}x) - \\sum_{i = 1}^n \\log (1 - x^2_i) \\to \\min_x\n",
    "$$\n",
    "$$\n",
    "f'(x) - ?\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Точное решение с помощью CVXPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "n = 1000\n",
    "m = 2000\n",
    "A = np.random.rand(n, m)\n",
    "\n",
    "x = cvx.Variable(n)\n",
    "\n",
    "obj = cvx.Minimize(cvx.sum(-cvx.log(1 - A.T * x)) - \n",
    "                   cvx.sum(cvx.log(1 - cvx.square(x))))\n",
    "prob = cvx.Problem(obj)\n",
    "prob.solve(solver=\"SCS\", verbose=True)\n",
    "x = x.value\n",
    "print(\"Optimal value =\", prob.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Решение с помощью градиентного спуска"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import cvxpy as cvx\n",
    "print(cvx.installed_solvers())\n",
    "\n",
    "# !pip install jax \n",
    "# !pip install jaxlib\n",
    "\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "\n",
    "# from jax.config import config\n",
    "# config.update(\"jax_enable_x64\", True)\n",
    "\n",
    "A = jnp.array(A)\n",
    "print(A.dtype)\n",
    "x0 = jnp.zeros(n)\n",
    "f = lambda x: -jnp.sum(jnp.log(1 - A.T@x)) - jnp.sum(jnp.log(1 - x*x))\n",
    "grad_f = lambda x: jnp.sum(A @ (jnp.diagflat(1 / (1 - A.T @ x))), \\\n",
    "                          axis=1) + 2 * x / (1 - jnp.power(x, 2))\n",
    "grad_f_jax = jax.grad(f)\n",
    "print(jnp.linalg.norm(grad_f(x0) - grad_f_jax(x0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Подробнее про jax, его возможности и особенности можно посмотреть например [тут](https://github.com/amkatrutsa/MIPT-Opt/blob/master/Fall2020/03-MatrixCalculus/jax_autodiff_tutorial.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "gd = methods.fo.GradientDescent(f, grad_f_jax, ss.Backtracking(\"Armijo\", rho=0.5, beta=0.1, init_alpha=1.))\n",
    "x = gd.solve(x0, tol=1e-5, max_iter=100, disp=True)\n",
    "\n",
    "x_conv = gd.get_convergence()\n",
    "grad_conv = [jnp.linalg.norm(grad_f_jax(x)) for x in x_conv]\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.semilogy(grad_conv, label=r\"$\\| f'(x_k) \\|_2$\")\n",
    "plt.semilogy([np.linalg.norm(x - np.array(x_k)) for x_k in x_conv], label=r\"$\\|x_k - x^*\\|_2$\")\n",
    "plt.semilogy([np.linalg.norm(prob.value - f(np.array(x_k))) for x_k in x_conv], label=r\"$\\|f(x_k) - f^*\\|_2$\")\n",
    "plt.semilogy([np.linalg.norm(np.array(x_conv[i]) - np.array(x_conv[i+1])) for i in range(len(x_conv) - 1)], label=r\"$\\|x_k - x_{k+1}\\|_2$\")\n",
    "plt.semilogy([np.linalg.norm(f(np.array(x_conv[i])) - f(np.array(x_conv[i+1]))) for i in range(len(x_conv) - 1)], label=r\"$\\|f(x_k) - f(x_{k+1})\\|_2$\")\n",
    "plt.xlabel(r\"Number of iteration, $k$\", fontsize=20)\n",
    "plt.ylabel(r\"Convergence rate\", fontsize=20)\n",
    "plt.xticks(fontsize = 20)\n",
    "plt.yticks(fontsize = 20)\n",
    "plt.legend(loc=\"best\", fontsize=20)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Pro & Contra\n",
    "\n",
    "Pro\n",
    "- легко реализовать\n",
    "- сходимость как минимум к стационарной точке\n",
    "- параметры при выборе шага влияют на сходимость не столь сильно\n",
    "- имеет многочисленные вариации\n",
    "\n",
    "Contra\n",
    "- линейная сходимость для сильно выпуклых функций\n",
    "- очень сильно зависит от числа обусловленности $f''(x)$, выбор начального приближения может помочь\n",
    "- не является оптимальным для выпуклых функций с липшицевым градиентом и сильновыпуклых функций (см. [ускорение Нестерова](https://blogs.princeton.edu/imabandit/2013/04/01/acceleratedgradientdescent/))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Резюме\n",
    "1. Методы спуска\n",
    "2. Направление убывания\n",
    "3. Метод градиентного спуска\n",
    "4. Правила выбора шага\n",
    "5. Теоремы сходимости\n",
    "6. Эксперименты"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8"
  },
  "nbpresent": {
   "slides": {
    "b738dfa0-30f4-4350-8338-c31d236608ec": {
     "id": "b738dfa0-30f4-4350-8338-c31d236608ec",
     "prev": null,
     "regions": {
      "96278b08-d857-478b-803c-5921bd08bbcd": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "196b8a50-3d29-45c3-82b9-f4a09b49491d",
        "part": "whole"
       },
       "id": "96278b08-d857-478b-803c-5921bd08bbcd"
      },
      "9b98a7cc-d9d4-4a02-ac7e-30168dcde7af": {
       "attrs": {
        "height": 0.4,
        "width": 0.8,
        "x": 0.1,
        "y": 0.5
       },
       "content": {
        "cell": "2a573842-172b-4931-b0dd-9c9d3c47a450",
        "part": "whole"
       },
       "id": "9b98a7cc-d9d4-4a02-ac7e-30168dcde7af"
      }
     }
    }
   },
   "themes": {}
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}