Optimizations methods HW Sergienko

170325_opt_hw2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Домашнее задание 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Задача 1.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Численно решить задачу\n",
    "$$\\min_{Bx\\le d} x^TAx - q^Tx,$$\n",
    "где $x\\in\\mathbb R^{10}$, матрица $A$ является трехдиагональной вида $[0, \\ldots, 0, -1, 2, -1, 0, \\ldots, 0]$ (двойки на главной диагонали матрицы), $q=(1, \\ldots, 1)^T$, матрица $B$ является нижней треугольной матрицей, все элементы на главной диагонали и под ней равны 1, вектор $d = (1, 2, 3, \\ldots, 10)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выпишем модифицированную функцию Лагранжа:\n",
    "$$M(x, y, K) = f(x) + \\frac{1}{2K}\\sum_{j=1}^{10}(y_j+ Kg_j(x))^2_+-\\frac{1}{2K}\\lVert y\\rVert^2,$$\n",
    "где $f(x) = x^TAx - q^Tx$, $g(x) = Bx - d$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Применим итеративный метод поиска решения:\n",
    "$$x^{k+1} = x^k - \\gamma M_x'(x^k, y^k, K), \\quad y^{k+1} = (y^k + \\gamma M'_y(x^k, y^k, K))_+$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так как при $x\\neq 0$:\n",
    "$$x^TAx = x_1^2 + (x_2 - x_1)^2 + \\cdots + (x_{10} - x_9)^2 + x_{10}^2 > 0,$$\n",
    "то матрица $A$ положительно определена и функция $f$ сильно выпукла с коэффициентом, равным наименьшему собственному значению $A$, которое больше нуля."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Согласно теореме 2 параграфа 3 главы 9 учебника Б. Т. Поляка \"Введение в оптимизацию\" (остальные свойства тривиальны для линейных и квадратичных функций), этот метод сходится со скоростью геометрической прогрессии при $\\gamma < \\gamma_0$ для некоторого $\\gamma_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Найдем производные $M'_x$ и $M'_y$:\n",
    "$$\\frac{\\partial}{\\partial x_i} M(x, y, K) = \\frac{\\partial}{\\partial x_i} f(x) + \\sum_{j=1}^{10} \\frac{\\partial}{\\partial x_i} g_j(x)\\cdot (y_j + K g_j(x))_+= $$\n",
    "$$ = 4x_i - 2x_{i-1} - 2x_{i+1} - 1 + \\sum_{j=i}^{10} (y_j + K (x_1 + \\cdots + x_j - j))_+$$\n",
    "$$\\frac{\\partial}{\\partial y_i} M(x, y, K) = \\frac{(y_i + K(x_1 + \\cdots + x_i - i))_+ - y_i}{K} = \\max\\left(-\\frac{y_i}{K}, x_1 + \\cdots + x_i - i\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def step(x, y):\n",
    "    mx = np.array([\n",
    "        4 * x[i] - 2 * (i > 0 and x[i-1]) - 2 * (i+1 < n and x[i+1]) - 1 + \\\n",
    "        sum([\n",
    "            max(y[j] + K * (sum(x[:j]) - j), 0)\n",
    "            for j in range(i, n)\n",
    "        ])\n",
    "        for i in range(n)\n",
    "    ])\n",
    "    my = np.array([\n",
    "        max(-y[i] / K, sum(x[:i])) - i\n",
    "        for i in range(n)\n",
    "    ])\n",
    "    return [\n",
    "        x - gamma * mx,\n",
    "        np.max([y + gamma * my, [0] * len(y)], axis=0)\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 389,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_problem(n):\n",
    "    A = np.diag([-1] * (n - 1), 1) + np.diag([1] * n)\n",
    "    A = A + A.T\n",
    "    B = np.tri(n)\n",
    "    q = -np.array([1] * n)\n",
    "    d = np.linspace(1, n, n)\n",
    "    return A, q, B, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f(x) = -8.74771372994\n",
      "Bx = [ 0.38419134  1.08353801  2.02890779  3.15113947  4.38110829  5.64978043\n",
      "  6.88825091  8.02776293  8.99970111  9.73553077]\n",
      "x = [ 0.38419134  0.69934666  0.94536979  1.12223167  1.22996882  1.26867214\n",
      "  1.23847048  1.13951202  0.97193818  0.73582966]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXl8nFd97/8+zzqrZkYaSbYlK/KexIkTOw6QBBogLAll\nvaWUlKW0QG75FWih8GtvL9ul/H4t3La0tJT1RbnQsiQ0pSGErZAQIAmJ42xe4sS7Zdnat9Fsz3Lu\nH8/MaGY0smRbsmTpvP06r7N9z3m+M7Y/z5nznOccIaVEoVAoFMsLbbEdUCgUCsX8o8RdoVAoliFK\n3BUKhWIZosRdoVAoliFK3BUKhWIZosRdoVAoliFK3BUKhWIZosRdoVAoliFK3BUKhWIZYizWhdPp\ntOzu7l6syysUCsVFyaOPPjoopWydzW7RxL27u5tdu3Yt1uUVCoXiokQIcWwudmpaRqFQKJYhStwV\nCoViGaLEXaFQKJYhStwVCoViGaLEXaFQKJYhStwVCoViGaLEXaFQKJYhi7bOXaFYKUjfx/d9ZCn4\nldibViZ9Hykb20+lS+0kSOmDBIkE30dSVSYlMjAq1Qex9CUga+ur4mnpmjbl/qvtAKaO66w5urMq\nPe1ATzlDm2lmjfubq52sv3KNWX1/1T7N0Gge6Lh0K93bts9rn/UocVdc1Pi+h1so4DoObrGAW3Tw\nnCKuU8QtFvGKxak6p1RXrM17joPnefiei+8Gsee6Qd7zgrTrzmjjeR6+O2UrSwJcFmXFMkGIeevq\n2lf/lhJ3xfJD+j6FbJZCNkM+k6GQnSQ/GaSL2UmK+TxOIY+TzwXpqtjJ52vSrlM8L19000Q3TDTD\nQNd1NN1AM3R03UAzDDR9Km2YBloohFay043pNkGdjtA0NE1DaBpCVKUr5XU2demZ6jQx1VZoAoEA\nAUJoIASiFMploiRIQtMqsYCSrVZlN1UmKv0BiKrrVPdfFZdsapnK12hinUCKKjvmbDdD39OuK2aq\nmrG/hu0uUpS4K84LKSVOIU92bIzs2EgQj4/WxLnxUfKZQMAL2QyFbPaMP68BDMvGDIWwQiHMUDhI\nhyNEk81BWThSqgthWDaGZQXBDIJumRimjWGZ6KWyso1umoGNaS6b/8gKRT1K3BUzIqUkNzHOxNAg\nE0ODZIYGmRgaYGJ4iImhATJDQ2RGhnGLhYbt7UiUSCJBOJ4g1txMy9ouQtEYdjQaxJEodixWSYdi\nQWyGQmiafoE/rUKxvFDivsKRvs/E8BCjp08x2tcbxKdPMdoXxE4hX2Ov6Tqx5hbiLWnaN2xiQ3ML\n0USSSCJJpClBJJEk3JQg0pTAsKxF+lQXF9IvPez0yrFfl59eX0n7MnjW55cebvpQeuoZPFStqafy\nEBS/lK6yD56dTpVV2pbrqvKB46WoLl9fX7luuayBvTxj+5LvDcqq7eVM7Wcqq25T84GoM6rPN35A\nK+dgUya6o43YDR0sJErcVwhSSiZHhhk8cawShk4cY7DnOG5hauStGwaJtlUkV61m7eVXkmhfRbwl\nTbw5TawlTTSRrMzfNryOL3GKHk7Rp5Ar4Lk+vifxXYnv+3iuxC+Ved5UXTktS6IiS8I0FQfCVRNX\n2fol4fGrRC4Qg9J//IqDVXFJ1IQvEZ5EK4msVpUXpRDYgSjbS0pxfb7apiovQZSuKSSIkloF6aWL\nBBDlePoctqzL1yWQDafia+fZG/dRVSBqfWlU17C/mabc5vKF1/s5SzvZwJVyybTvAHBGCsTm4Mb5\noMR9mZIZHuL0oWdL4Rn6Dh8kn5mo1EcSSdJrL2Hbi19Oak0n8ZZ2rGgaw0rg5n0KOZdC1qGY8xg6\n5XDqkEshdwqn0INb9HCLPk4pdoteJe05C786xABMAYYAU4hSDKYmamJDgEFQrwvQCYJRlT7XOXev\nNMCtD16DMp/S4LkUy/q8bGxTYy+n11UG0iWfZF3eZ+rGVj9grrnHle1lbXn96HPaU5IzD05nHb3O\n3t8ZGpzltWe71oVmx6YUnQt8jVnFXQjxFeCVQL+U8ooG9W8C/qyUzQDvklI+Ma9eKs6I9H0Gjh+l\nZ/8eevbt4dSzT5MZGQaC1RHNa7pYs3kHofhqNLMVobXgFCyy40VOPFPkwKNFPGcMGGvYv9AEdtjA\nCuuYIQPT0jAsHTtqYlgapqVjWDpGqdywNHQjCJoupmJdQzNKsS4QAnTHQxQ9RNGHoofIe8i8i8y7\nkPfwcy6yOhS82b8QASJkoNk6WshA2DrC0hBmEGuWjjA1hKWXgoZWqivnhaUH9qaGMDSEIUqxBrpQ\nD2KXOdPWv8/zzeRC/POZy8j9q8A/AV+bof4IcKOUckQIcQvwReC58+OeohFSSoZ7ezj6+G5O7HuK\nk0/vrYzKQ/EWIk3dtMSfi+umyeeSTE6aTB4pNRYQjueJxH0iTSbJ9gSRuEW4ySLSZBGKmFgRAzts\nYEcMrLCBaetzFjMpJTLv4Y0V8CaKeBkHf6KIN17EyxSD9ISDN1HEz7nMNM4Xto4WNdEiBnrcQmuP\noEdMRLhKtEO1sRbSEbYRiLMSX8V5cKZllA2yDUsWm1nFXUp5vxCi+wz1D1RlH4IF/7WxInGLRU7s\ne4rDux/h0KOPMDHYB4BpN4PejRnpQDM7QWvCw6C5LUqiNUy8JVQKYeLNIWIpG904910npCfxxgt4\no0FwRwt4o/mqdKHh6FqYGlrcQo9bmK1h7PWJQLRjJnrURIuYaFEjiMNGMEJWKBTnzHzPub8d+ME8\n97licYtFjjzxKHvvu4+jjz+C5xZBGGhGF0bkJgx7Pa1da2i7JE5LR4zU6ijNq6OE4+e3flt6Pu5I\nAXcohzuYwxvKV9LuSCGYnK1CixroyRBGS5jQhiR60kZP2IF4x030uBVMjajRtEJxwZg3cRdCvIhA\n3J9/BpvbgNsAurq65uvSywrp+xx98nF2/+DHnNizC8/Ngwijm5cSa91C15Xb6NjcSnt3Ey2dUQzz\n3NeDS8fD6c/h9Gdx+yZxTmdxBrJ4I3mq50uEpWOkQ5hrYoSvbMVoDgUCXgqapdakKxRLjXkRdyHE\nNuDLwC1SyqGZ7KSUXySYk2fnzp2L/Lx6aTE+NMhD//49DjzwM4q5ERA2ur2Jji072XL9TrouS9O8\nJnpOo18pJd5IAac3Q7E3g3M6i9ufxR3KTT3o0QRGaxhrTQxjWytGSxgjHYzGtZh6k1OhuNg4b3EX\nQnQBdwJvkVI+c/4urRyklBzc9QS/+vbtDJ14CpDoVhdrt93Ctpt+g/VXr8IKnd1fkZQSbzhP8WQG\n52RJzE9m8LNuYKCB0RLGXBUhfFUrZnsEsz2CkQ4jdDXPrVAsF+ayFPKbwAuBtBCiB/goYAJIKT8P\nfARoAf65NLpzpZQ7F8rh5YDve+y6+2fsuvs/yI0dBxEiueb5bH/5zVz5oisw7blPc/hFj+KJCYrH\nxykeC+KKkOsCc1WU8BVpzDUxrI4Y5qoowlQirlAsd+ayWubWWerfAbxj3jxaxviex6/u+AG777kD\ntzCE0JOs2/5b3PiW19LSkZpTH95EkcLhMYrHxikcG8c5lanMjxttEUKXt2B1xbE64pjtEbXqRKFY\noag3VC8Avu/z8H/+jIf/85s4uT50K822l76TF/zuzYQi9pnb5lwKh8coHBolf2gUty8LBEsLrbVx\n4jeuxbqkCbsrjhYxL8THUSgUFwFK3BeYg4/u58ef/yy58aNoRoptL307L3zrKzGtxkIsPUnx+Dj5\nA8PkD47inMwEe5CYGlZ3E5HtbYQ2JDHXRNUcuUKhmBEl7gvE+OAod336C/Qd/CVC2Fz6/Dfyknf8\nNnZ4+kjdm3TIPzNC/ulh8s+MIHMuaAKrK078xV2ENiSxuuIgPfxsFj+boXh8EHw/OFqtHHs+SB+h\n66AbCNNA6DrCMMAwEOVgmohQ6IwbgCkUiosbJe4LwK/u+BG/vvMrSH+S9CXX8er3/SGp1S01Nu5o\nntxTg+T2DFE8Ng6AMH00ewxpncQfPUj2/tOM3zmIPzaGn80iHWde/SyLvBYKlWIbYYcQIRvNDiHC\noSAu58s2oXDFVguHgri6PBRCs21EOBzEoRDCttVyymVAeTdO3/enduasbLkra2wapWerX4p91dvM\nVHY2de3t7XR0qC1/LxpG+4b597/6NKOnHsMMreKlt/0Fl91wVaXeHc0zuauX7K6TeKNBmZ/pxel5\nFPf0k/ijxwGJlkhgpNMY6TThK69ETyTQolG0aAQtEgRh21A6eg1ND6ZotOCoNXwf6bhIzwXXRbou\n0vWQrgOehywW8fMFZCGPn8sHcb6AzOfxC3lkvoCXmUAODpbKSnX5PDKfb/zh50BF9KtuKCJko4XC\nVTeQunz1DSYU3EiEaSIsE2FatWnLRJgmmmVBKRZmUMYCnrrkeR6e5+G6bk2YrczzgnNWy6E6fy51\n1WLbSIDnYjNbO8X8cMMNNyhxv1jYdff93P+Nf0Z6Wbq23cJr/vSdWCEL6XiM/WQP2YdP4eejAHij\nx3FP7gKtn9CWDiKv2ITV/RKsS7qxutaihcOL/GlmRkqJLBZrxD64MeQa3DCq60s3jtzUDaRSX8jj\nD2ZwppUXYB5+rfhC4JgmbjiMGwnj2SG8kI1r2XiWhWeZuIaJZxq4hoGr63iajqNreJqGq2m4QuBq\nGg4CT4AnBB7BFr/zLXmaEOiahiYEWk2soWmlugZB14IN00QpnjrDVQTnupbOPdWq7Ur9iyrbahut\nqmxa++ozVWFaulHZudpeyL7qy+rLq+O5ltXX2faZF1LMB0rczxPXdbnzrz/Piad+iGGnufk9/5Mt\n111Jbu8xTn33CdzRMEIP4U9m8XNPEN4cJfmq7YSvfBV6U9Niu3/WCCGCXw22jZ5ILPj1pOtSzEyS\nHRslOz5OdmKCbDZLLpsjn89RKBSC4DgUikUKrkvRdSl4HkXPo+j7uGcx4tSkxPD9qeB56I6D4XrY\nnovuuuieh+56aK6L7rlojovmuWhFB9330Dw/sCmnG5Rp/lQQUk7FS210XDoEu/yrUFSXVQXRqKxs\nW/5FWT64mzrbqgO4a0PJtr69EFBvX/GLqcOvp9WJqgM/BFLMZFdKMr39VL5asBvZTX1/Df152csI\nv+618/AXNDNK3M+D/mOnueMvP0F+4ijNndfy2x9+P96v9tHzwdtBX430osj8UUJXJkm+9gWYbf9t\nsV1edKSU5HI5MpkMk5OTZDIZMpkM2Ww2EO1cjlwuV5N2Zhm927aNbduEQiHsWIy4bZMu50t19cE0\nTUzTxLKsmljXz2+fHOn7SLc0HeZ5dWkPvFLacUH6SM8LNmLzveDBuJTB1Fl1Wd0D80ob2aCsqk0Q\nV5WVbSunesjSPuRBHFy/trxyLJ8k8KNsW99Ho36r+/AblFXbl/pu1IeU/sz91pxUUq4L8pXDPs5k\nR9V0U82JJnV9VY72mqqTdXYSeeY+qs7j8yfGz/Wf2JxR4n6OjPYN829/8X58N8tVL/t9drauYuRj\nP0KE2pFeDCPeQ+rW6wltummxXb0geJ5HJpNhfHyc8fFxJiYmKsJdLeSTk5P4/vRd3IUQhMNhIpEI\n4XCYRCLB6tWrCYfDNeXV+VAohGmaaEto1Y/QNIRlgTo/VrHIKHE/Rx7/8S/x3XFufMHbad9TJGP7\nSGFjd47Q8ns3oceji+3ivOH7PplMhtHR0Yp4l8PY2Bjj4+NkMplpD9w0TSMajRKLxYjFYrS3t1fS\n5fJyHAqFlpRIKxQXO0rcz5HeA8+yo+VltJ9oBpHFbB8i/c6Xo8eW7sPQM5HL5RgdHWVkZISRkZFp\nac+rPYDDNE0SiQRNTU1s2LCBpqamaSEcDqvljwrFIqHE/RyxB/JsSj0HofXQ9oFbMFuTi+3SrLiu\ny/DwMAMDAwwODjIwMMDQ0BAjIyPk65Y4hkIhUqkU7e3tbNmyhVQqRTKZrAi6rdatKxRLGiXu54Dv\neqS0GJ7vsvZ/vR7NXlp7uuTzeQYHBysCXk4PDw/XTJ0kEgnS6TQdHR2kUqlKSCaThJfwckyFQjE7\nStzPgeN7D9Nst5N1RhZN2KWUZDKZGgEvxxMTExU7TdNoaWmhra2NrVu3kk6nK8FSD/0UimWLEvdz\n4NCup9hgtVHQ+hb8Wr7vMzo6Ok3ABwcHa6ZSLMsinU6zfv36ini3traSSqXOe3mfQqG4+FDifg6M\n7zuGYbYS3jB/8+yFQoGhoSGGhoZqplSGhoZqHmZGo1FaW1u54ooraG1trQh5U1OTmgNXKBQVlLif\nA8akA0mIX7txzm2KxWLNuu/x8XEGBwcrYl49lQKQSqVIp9Ns3LixMgpPp9NqLlyhUMwJJe5nieu6\nNGlRHK+AfnkXR/qOMDA6wERmguxkFj/v4+Zd3Jxb8wZmsVic1pdt26TTadatW0c6naalpYV0Ok1z\nczOmubQe0ioUiosLJe5nyfEnD9Fst9Pv9fHFT30M252+AVBRK1LQC/iWjxkyia6O0p5spyvdxbq2\ndaSaUsTjcaLRqJpKUSgUC4IS97Pk8MNPsNlq5wntYWzXpvmqZtLNaSKxCHbYxjVdsl6WvmwfJzMn\n6cn0cHj0MJnxDIyDccRgU2oTO9p3cHXb1exo20FbpG2xP5ZCoVhmKHE/SzLP9KBZq+kNjZOzi7z3\nde+dtY0vfU5OnGT/8H72D+/nyYEnufPZO/m3/f8GQEesgx1tgdhf034N6xPr1YheoVCcF0rczxIz\n6yMtyQSScHpuDzc1obG2aS1rm9bysu6XAeD4DgeGD7C7bzeP9T/GA70P8L3D3wMgaScDoW+7hu3t\n27m8+XJMXc3BKxSKuTOruAshvgK8EuiXUl7RoF4A/wC8AsgCb5NS7p5vR5cCbtEhYcQ57Q+iS4Pu\nS7rPuS9TM7kifQVXpK/grVvfipSS4xPH2d23m939u9ndt5v7TtwHQEgPsa11G9vbtgfTOa1XEzEj\n8/OhFArFsmQuI/evAv8EfG2G+luATaXwXOBzpXjZcfSJZ2m22nlG9gJw7ZZr561vIQSXNF3CJU2X\n8LpNrwNgMDdYI/ZfeupL+E/66EJnS/MWdrTtYEf7Dralt9EWaVNTOQqFosKs4i6lvF8I0X0Gk9cA\nX5PBpiUPCSGSQojVUspT8+TjkuHIr5/gMquTU/peinqRyzovW9DrpcNpXtb9sspUTqaY4cmBJ3m0\n/1Ee63+MO565g3/d/68V260tW9naspXLWy5na3or6XB6Qf1TKBRLl/mYc+8ATlTle0ply07ccwdP\nQ6iTfjOHSIkLvv94zIpxfcf1XN9xPQCO57BveB97Bvewb2gfewf3cn/P/ZUTX9oibVzecjmbU5vZ\nlNrEpuQmupq6MDU1f69QLHfmQ9wbzQU0PAhSCHEbcBtAV1fXPFz6wmIVIBsqIDFo72hebHcwdZOr\nWq/iqtarKmVZJ8v+4f3sHdzL3qG97Bvax/099+PL4PQjUzNZl1jHptQmNiY3sjG5ka6mLtbG1qqH\ntgrFMmI+xL0HWFuV7wR6GxlKKb8IfBFg586dS+wk4DNTLBRJGE0cl6cB2LZx2yJ71JiIGeGa9mu4\npv2aSlnBK3Bk7AjPjjzLs6PPcnDkII/2Pcr3D3+/YqMJjTXRNVzSdAldTV1BHA/i1dHVSvgViouM\n+RD3u4B3CyG+RfAgdWxZzrc/doBmu51nxbN4wuN5W5632C7NGVu3ubT5Ui5tvrSmfKI4weGxwxwf\nP86x8WOV8PjA40w6kxU7gSAdTrM6tpo10TWsjq1mdbQ2HTNj6oGuQrGEmMtSyG8CLwTSQoge4KOA\nCSCl/DxwD8EyyIMESyF/f6GcXUyOPfgElxuX0GeO4cQcQlZosV06b+JWfNq0DgR7xQ/lhyqif3ry\nNL2TvZyaPMW+oX389PhPcXynpk1ID5EOp2mNtJIOp4N0uCodaaUl1EIylFRz/orzJli/4VdikA3K\nQMpyXSlGQtkOCbJsf5ZtZ+yvukxW2ZT6LNWFw51EIusW9Duay2qZW2epl8AfzZtHS5TC0QHccCcZ\nTZJoTyy2OwuKEKIiyjvad0yr96XPcH6Y3kwg+KcnTzOQHWAgN8BgbpBDo4d46NRDTBQnGvQOMTNG\n0k4GIZScSttJUqEUCTtB0k4St+LEzTgxK0bMii2pm4KUHr5fwPPy+H4B3y8gpYsvXaRfDNK+i5QO\nvnSQvlsqC+qqbYM4KAv+83tI/EB4KgLkTaXxSyLiBSJTbVdp59eK07T21W2q4opAVQtmVbtpZbKq\nrE4My2WV/qYL5vS2M4l2bZuLnUu6/jsbN/6/C3oN9YbqHAk5Ov2RUYQQbF63ebHdWVQ0oVXEf1vr\nzM8e8m6ewdxgTRgtjE6F/Cgj+RGOjB1htDBaMxXUiJAeCoTeLAUrRtyKV9IxI0rYsIjqBmFdI6QJ\nbAGWBpbwMfEx8NDx0KUDsoDnF/D9PL4XCLTnl8W6GJSXhHtKxINyKd35/lrr0BBCK8WiFOtVaQ0Q\nCKEh0EEIBBoIbSquah+UVbVHBHmAqr4qfQojaIMI+i7XoYEQpbSoa0ulvtauqqx0vdq2DeymlYmS\nL3Vty/5VPnvV5zlj28DX8ueasqv1dapeq71W1fc1vW15erLWvrrOtlfP7z+XBihxnwOFbIGE0cQJ\n2Q/A9Zdfv8geXRyEjBCd8U46451ntJNS4nlZcsVhRrOnGM2fIpPvJ1sYpuCMUHDGcN0JPDeD702C\nP4yQp9BkEd1xMFwPI+djCYleN+3vlEKj24YnwZUCF4GHhoeOj4aPjhQGEgMpDBAmaBGESCA0C2FY\naJqNroXQ9BC6ZmPoYQzdxtBsdM1C120MUYo1G0O3MbUQhh4K7HQbUw+VQjiINRtNM5gSPIXi3FHi\nPgcOPbqXZnsVj2iPk7fzrE4t/F33YiQQ6QyOM4rjjJTiURx3tK5sBNcZK9WN4XkZpPQa9mmXAmgY\nRgzDimMYSQw9jm7E0PVIJWhaGCksPKHjYeCi4Uidgg8FCQVfkvMlOd8n63nk3CI5L0/OzZFzchT9\nIgWvQNGbiqvTBT9byZeXls43hjAwtCCYmllJl4MudDSh1cRCiGnlmlaVRkPXSum6tuWyaeWl9uWR\npiaCXwAaGgjQSr8ItMooOEiL0o2pki+lq8sq/dXnG9ieKV99/UodIhgkV90gy+U1+bp0+XPWpKvr\nxFQ/0+yEaFwnZrg2gia7iYS9sNO7StznQM+DT3GZvo4hI4ednr5/+3LG83IUi4OlMDSVdoaqyoYC\nwXbHzjhdYRhxTCOFaSYxzSThSDemmcDQ4xhGEHQjVknXlOtLa+9713crQl9zI/ALuL6L4zm40sX1\np4LjO2dOz2BfbetLH096NXE5eNKj6BeDcj/IS2RNvr69lLJhf54fpGX5T/lBoWJe+IMr/oD3XfO+\nBb2GEvc54JwYYTTWjq8J1q498xTDxYCUPo4zTKHQVxuKAzjFKtF2hvC8xvPghhHHstKYZgvR6AZM\nM1UKiRoBD8qSGEaiNOWwPCiPpFfaBm5SyinRLwl+fV5KiV96KHomm1lty3bVNlVtkDS8AVXfiBql\np9VJavoAprep7qdBm4Z1Da5XttmQ3LDgf1fL53/bAhLyDPrECAA7Nk9fPbKUcN3JklifplDsnybg\nxZKIS+nUtRSYZjOW1YJlpWlKXB2kzTSWlS6Vt1QEXddX1i8YRUB5Ckix9FHiPgv5yRwJM8FeBnA0\nh6u6r5q90QIhpcR1x8jnT5LPnyRXivO5HvL5XnL5k7ju6LR2uh7Dtldh220kU8+tpG27Hdtqx7bb\nsaw02hJaaqhQKM4PJe6zcGjXPpqtVQxox/ATProWjFome37G/iffw+rWV7Jm519VloKdD1JKis5Q\nRbzLoh0IeZD2vExNG12PEAp1EAp10JS4OkjbqwLhttuxrDYMI3revikUiosLJe6z0PvQHjboXWR0\nh+ZVU5uF9T79t4yF8oxNfIeee++lY/OfkE6/kFBozYx9+X6RQqGvJNi95Au9FfEOQi++n69pYxhN\nhEIdhMNdpFLXEQ51BgIe7iAc6sAwkkvqQaNCoVgaKHGfBefkKP2xBAjBpeum9mYZdJ6l2Y+w2uvm\nqHiSA898mAPPgGEkCdnt6KXRsu87uG6w7M91x6f1b5othEJriEW3kG55UUm4AwEPxDt+wT6rQqFY\nPihxn4WQa9DHCBLJdZdeB0Du9INkbY8O87msev6XaL/3E2Qe/TSjG7YyueUGCt4IvpcDwNBjRCPr\nMcwEppkiZK8mFFpDKLQG216Nrl/8e9QoFIqlhxL3M1DMF0iYSQ6I0+RDedJNwclGgwe/CkB6w1uD\nV4pf/GHiyUuI3/0+OO3B734bkmvP0LNCoVAsLBf2KKGLjCO795O02xjSM4RbwpXyofFfEy5oRNb8\nxpTxjrfCm74DYyfgyzdB72OL4LFCoVAEKHE/Az2/fIqc7uFqks7Sy0teboARa5y00eAlhA0vgrf/\nGHQb/uUV8PT3p9soFArFBUCJ+xkonhylXxsDpl5eGj7wL/iaIL36VY0btV0G7/xpEH/rTfDgZ0t7\nOysUCsWFQ4n7GQi5BqfFaM3LS0N9P0T3JMnNvzdzw1gb/N7dcNkr4Ud/Ad97Lzj5me0VCoVinlHi\nPgNu0SFhJDgthisvL0nfZ1Aeo9lpRrNiZ+7AisBvfw2e/37Y/TX4ystg5OgF8V2hUCiUuM/A0cef\nIWo3M64VKi8vZY7fQ8GCdOr5c+tE0+AlH4U3fhOGj8IXboRnfrxwTisUCkUJJe4zcOLnjzNi5EHA\nlu4tAAwd/SYALZvfcXadXfoK+O/3QWItfOO34Yd/oaZpFArFgqLEfQaKJ0cYEONIJNdfFpy8NJR9\ngljewG654uw7bF4P7/gJXPtOeOiz8MUboffxefZaoVAoApS4z4DtmJzWRijYwctLbqaXMStLi73l\n3Ds1w/CbfwNvvhPyY8F6+Pv+GtzC/DmuUCgUKHFviO96NJlN9ItxQung5aWRZ/8PUhO0rHn1+V9g\n403wrgdg6+vgvr+Cf74ODt17/v0qFApFCSXuDTj+5EE0O4ajeZWXl4b6f4LuShIb3jg/F4k0w299\nORjFI+FlUyIzAAAgAElEQVTrr4Xv/AGMnpif/hUKxYpmTuIuhLhZCHFACHFQCPHnDeq7hBD3CiEe\nE0I8KYR4xfy7euE4fu+jDBtZIHh5Sfo+Q/I4zd4clkCeLRtvgnc9CDf+Oey/G/7xGvjxhyA7PL/X\nUSgUK4pZxV0IoQOfBW4BLgduFUJcXmf2IeB2KeV24I3AP8+3oxeSQs8ofWIMH5eruq8i23sfeUvS\n3PTchbmgGYIX/Q94z6Nw5evhgX+Cz1wNv/hbyE/fJlihUChmYy4j9+cAB6WUh6WUReBbwGvqbCTQ\nVEongN75c/HCYzsmfdooTtJF13SGjnwDgJYNb17YCyfXwmv/Gd71K1j7XPjpx+HTVwRxZmBhr61Q\nKJYVcxH3DqB6IrinVFbNx4A3CyF6gHuA9zTqSAhxmxBilxBi18DA0hQr3/cJGzFGRZbmVS0ADE/s\nIlLQCK+67sI40b4V3nQH3HYfrL8RfvF38PdXwt3vg769F8YHhUJxUTMXcW90hlv9Tli3Al+VUnYC\nrwC+LhocKiql/KKUcqeUcmdra+vZe3sB6NlzCCdkBS8vrduClx9hxBynRW+wC+RCs2Y7/M7X4Y8e\nhit/Cx7/BnzuevjKzfDkHWoJpUKhmJG5iHsPUH3yRCfTp13eDtwOIKV8EAgB6flw8EJz/L92M6RP\nggxeXhp99uv4uqBl1c2L51TrZnjNZ+H9++Fln4CJ03DnO+Bvt8D3/gSOPQC+v3j+KRSKJcdcxP0R\nYJMQYp0QwiJ4YHpXnc1x4CYAIcRlBOK+NOddZqFwYpg+bQzfcEg3pRk6/QO02XaBvFBEmuH698B7\ndgdLKDe+BJ78NvzLLfAP2+AnH4WeR5XQKxSK2Y/Zk1K6Qoh3Az8CdOArUsq9QoiPA7uklHcBfwp8\nSQjxPoIpm7dJeXFuYm46BgPaEKLdBGDIPURSNqGHUovsWRWaFiyh3HgTFDJw4B548nZ44B/hV38P\nsXbYfDNseUUwZ2+GZ+9ToVAsK+Z0hqqU8h6CB6XVZR+pSu8Dbphf1y48vu+j22EKwqWjs4Nc38Ol\ng7CvWWzXZsaOwbY3BCE7DAf/KxD7PXfC7v8TnArV9VxY9xuw7sZgHl83F9trhUKxwKgDsqs4tf8Y\nBSv4SrZt3Mbwoa8D0NI9T2+lLjSR5imhd4tw9Bdw6Gdw5Ofws08AnwArDl3PC5Zadu6Ejh0QSiy2\n5wqFYp5R4l7F8Z88wpCeQUjJjvU7ePrguwlpgkjHTYvt2tljWFNTNwCTQ4HYH/k5HP0VHPxJyVBA\n66WB0K/ZDquuDI4ItOOL5rpCoTh/lLhXkTs+Qn9sEmk5GPgMG8O0yy6Etgy24Im2wNbXBgEgNwq9\nu6FnF/Q8Ak/fDY99fco+tS5Yb7/qSmi7HNKbgm2LDXtx/FcoFGeFEvcqDMdgSExgrY0ycfhOPEPQ\nnHjhYru1MISTsOHFQYDgEO+xE8FLUqf3QF8pPP19pl5rEJDsgpaNVWEDJC+BRGewjYJixSGlRBL8\nK/ElSGQpLoVSvV9K+3W2smTrl9JT9aV+ZVDn19mVbWr6q7pWpU7Wtanuq8b3avsz1c3cV3mdmi+r\n/Z3+Wa5NRHl+amF/HStxr8K3TXwh6e5ax/DJu0BKUhvesthuXRhESbiTXbDllqny4iQMHIChQzB0\nsBSehRO/hmKmto9oa3DaVHJtKe4KRD++CmKrgnrDmhd3pZQUpaToSwq+pOj7FKUk7/sU/XJ5KS2r\nbHxJQUoc38eV4EqJ40tcGQSnHPsST1LJu1XljpR4ZVufaW3LduVbYiWW5fz0hWSi9K6gEFNvDda/\nPShEnW2dXW1bUblWvSgCtULF7EJbEc0G9Yqz591dbUrcLxSnDxwnZwXTL8/Z8hyGH/v/iWNhJRfh\nzdSlhBUNHrp27KgtlxIyfYHYj55AjvaQHz9JbryP3HAvuRN7yUlBXrPIaSFyuh3EoWZykdYgtpPk\nrQQ5K07OiJDTQoG9MMkJg5zUyNUI9JRwFxdgpa0pBIYQmBoY5XSlTKBX5Q0R2IR0gSG0KVstiHUB\nupiS52lCXCXdZbEvjxLLXy/V+Sqb6vqpPqrrpmw1EVxJMCX+GgIhgpdcBAKt5ErZVivZaiUvK+1q\n6kXlJRlNTPVTvlbtdUXpWlO2jfucuiYC9Co/q/urtKn7DNXX0ur81xr0Vf4cNX2VfS1/puq+Gvhc\n7Ytg5r60mrqp73whUeJe4tgPH2ZIz6D7kkuamviFlaWLrYvt1oIhpSTvS8Zdj3HXY9LzyXgeWc8n\n4/lMeh4Z159Wnplmm2LSa2KSy/DiwFkORmy/SChfIOznCXsZwn6ekF8g7BVo8wvYAmxNYGsalm5g\n6waWYQRpw8QyLCzDwjZtLDOEZYWwrQiWFcG2I4GNJrBEqQ9NYJcEOBDhKeFWKJYTStxL5I6PMBAd\nR0Zcxg99A6kJUm0vX2y3ZiXr+Qw5LsOOy1AxiMdKgl0OY67HhOtPK5/r6DesaUR1jZhRinWdZtOg\nK6RXlQfpiK4R1jTCukZYE4R1jVAlr5XygoimEdI1dID8aLBGPzsMuWHIOpDLTeVzI1V1I1CYgMI4\n07c4aoAZASsWvMhlRYPYjFTF5XQpb0Xq6sPBuwKGVRfboFtBbNhBmW4FL5gpFEsAJe5lHMGYlqVp\nfQvDfT9GE5LkxjctmjtjjsvJgkNvweFUoUhv3uFUwaGv6DBUdCuCnvNnFriwJmgydJoMnYShkzJ1\nusNWTVk5HdXLAq5XBLws1gs+qg2ngtByFlNgUgbPAwoTVWGsLl+6CRQmwMmDMwlOLgiZ/qm0k50K\n54tm1gp/JbaDl8d0M7DR9FLamArV+Uq63rZBW6EFZUIDodfltRnqG9lW181UL6byiFJ5acKkEteX\nMTe7Sp/VZeoX1bmixL2EawdvbW5Zt4XhgS+QIHpBthwYc1yemMjxVCbHwWyeQ9kCh7IFhhy3xk4D\n2m2Tdsuk1TK5NBai2TRoKQfLoNk0aDZ1EoZBk6FhLedRpBDB27l2DFg9P31KCW4eimWxLwm/Vwx2\n4PQKwcthNXHhzHXltm4BfAc8B3yvdJ3MVL66znfAd6fnfXf2z7AsqRf8RjeB+jLmaFeeyae2bf11\na+rFmetntGWq/uo3wfP+cB6/o+kocQcGDp4gawXLAK5Z3cLeCZcN2sJsOTDquNw/kuHe4XEeHp3k\nUG5q295Wy2BD2ObmdBPrIyHWhizW2CZrbJM2y8S4EE9hVjJCTE3R0LLY3kxHyimR96oEX/pB8L1S\n2ivZVufL9Weq82foq5yXtXlkqcyfSlfKZIOyRnZzaevPoYw52tWXVb7cqbLyd13d95zr69Mz2F6A\nlwSVuAPH7nmYAW0C05dw6m4Amte+bt76z3k+Pxwc4/bTw9w/MoEnIWHoPC8Z5Q2rmrm6KcK2eJiU\nqf46FGdAiKmpHbUZnGIWlJoAk8eGGYiNI+KS4aH7MTRJvLv+JMGzZ9z1+JeeQb7Q08+w49Fhm/zR\n2jZemk6wPR5RI3GFQrFgKHEHHA/ywmHV5lUMc5KU14LQz/1lGykl/943wscO9jLouLy4Oc671rZx\nQyqGph4QKRSKC4ASd8Ap7QS5tV2nMAnN9nPOua9Rx+W9+4/z46FxdjRF+Pq29WxvisyXqwqFQjEn\nVry4jxzpJWN5CAmdxd0cBprX/c459XU4W+DWJw7RW3D4y40dvL0zrUbqCoViUVjx4n7k+79mQJvA\n9mFi4lFCmiC86vln3c+zk3le//hBHCn5j+0b2ZmILoC3CoVCMTeW8ULouTF+ZIghMYHeDCPGMM10\nnPUWv/0Fh9954hA+cKcSdoVCsQRY8SP3vOfhCZ+ujXlcQ5CKv+Cs2hd9n9/fc4QRx+OuHRu5NKqW\nqCkUisVnxY/c83bwFay3jwCQWn92R+r9/bE+Hh3P8veXreXKuHpwqlAolgYreuQ+3jvIpOlh+gLD\nfYaIp2G3XDHn9nsmsnzmWB+vb0/xmraF36pAoVAo5sqKHrkfvutBBrRxbOkzZoyR0tbOua2Uko8c\n7CVhGPzlpo4F9FKhUCjOnjmJuxDiZiHEASHEQSHEn89g8wYhxD4hxF4hxDfm182FYejQKUZFlsTq\nQTxDkGq5Yc5tP/PkCR4YzfDmlqTaNkChUCw5ZlUlIYQOfBZ4KdADPCKEuEtKua/KZhPwP4AbpJQj\nQoi2hXJ4Ppl0XRDQ2dkDQHLD3Lb4dTyfv+vpRwDHH+uHyzoX0EuFQqE4e+Yycn8OcFBKeVhKWQS+\nBdRvvPJO4LNSyhEAKWX//Lq5MOTt4AWjZv0YkYKOnbp0Tu0+/MgRChGD1X0FfvjUKQ4PZGZvpFAo\nFBeQuYh7B3CiKt9TKqtmM7BZCPErIcRDQoibG3UkhLhNCLFLCLFrYGDg3DyeJ7KDo2RMl4gnyFnj\npPSuObWbKDr868go4ZzHD15/Dbah8Y8/O7jA3ioUCsXZMRdxb/T+fP3xPwawCXghcCvwZSFEcloj\nKb8opdwppdzZ2tp6tr7OKwfveoBBLUMiOoinC1Itc3sr9YO/Powb0vnjzjbam0K87fp1fPfxk+w/\nNb7AHisUCsXcmYu49wDVy0g6gd4GNv8ppXSklEeAAwRiv2Q5/fQJsqJActVxAFLrb521zWje4a7s\nJImsx3uvCH68vOvGDcRtg0/98OkF9VehUCjOhrmI+yPAJiHEOiGEBbwRuKvO5rvAiwCEEGmCaZrD\n8+nofDPpOgC0NB8jmtexUltmbfOBhw/hWxofXLcKrbRFQSJi8q4XbuTeAwP8/JnFnWpSKBSKMrOK\nu5TSBd4N/AjYD9wupdwrhPi4EOLVJbMfAUNCiH3AvcAHpZRDC+X0fJC3BAIfPTxOyrhkVvvRvMM9\nuSzJrMc7Lq995PD7N3SzPh3lw9/dQ97xFsplhUKhmDNzWucupbxHSrlZSrlBSvn/lco+IqW8q5SW\nUsr3Sykvl1JeKaX81kI6fb7kRjNkTJfWyAi+LkilZ99Ppjxq/8C6VdPqQqbOJ153BceHs3zmp88u\nhMsKhUJxVqzIN1QPfe9XDGkZkolTACTXn3l9+5lG7WWu35Dmt6/p5PM/P8Suo8Pz7rNCoVCcDStS\n3E/uPUJeODS1HCeWN7CSG85oXx61/2n39FF7NR951eV0piL88bceZyzrzKfLCoVCcVasSHGfcByE\n8Agnhkga3We0rR61v3PrmfeQiYdMPnPrdvrG8/zpHU/g+/UrRhUKheLCsCLFPW9CPDYMOqRaf+OM\ntnMdtZe5em2SD/3mZfzX/j4+qZZHKhSKRWLF7XiVz+aYND1a48GyxeS6N8xoO150+EFukoTLrKP2\nan7v+m4ODUzyhfsP09US4U3PnX01jkKhUMwnK07cj3zvAYa0DKsSp4kUNKzkzO9afWTXUTxL592r\nW87qGkIIPvqqy+kZyfKh7+7BNnRef43aXEyhUFw4Vty0zJEnnqEoikSSfSS1mUfjRc/nzrFxwjmP\nP7ri7PdrN3SNz735Gm7YkOaD33mCO3f3nI/bCoVCcVasOHHPFIpEo6Nopksy9dwZ7T752DGKIZ23\ntqYqb6OeLSFT50tv3cn1G1p4/+1P8Ln7DiGlesiqUCgWnhUn7nlTkGgqzbd3/beGNr7v89X+Ycy8\nx//c0X1e1wtbOl9527W86qo1fPKHT/MX/7GHouufV58KhUIxGytqzr1YLJIxXdqaBrCLEGq7tqHd\nl/b1MhnWeYMVwdLP//5nGzr/8DtX05kK87n7DrGvd4x/+t0drG1WB2orFIqFYUWN3A/94CGGtQni\nidMkZRtihumWfzjWj170+MS16+bt2pom+LObL+Xzb97B4cFJfvMzv+CuJ3rVNI1CoVgQVpa4P7IH\nIzyOEcqRbNre0ObfD/YzHNF4sR2hyTLn3Yebr1jN99/zAta1xnjvNx/jnV/bxamx3LxfR6FQrGxW\nlLhP5Is0JYITAJOdr2po88lnTiJcn0/unL9Rez1dLRHufNf1fOg3L+OXBwd56d/dz2fvPah2lFQo\nFPPGihL3ggHJpgEMVxLtfOm0+icGJjhuw1XSYE0stKC+6JrgHS9Yz4//5Eau29DC//7RAV70N/dx\nx64TOJ564KpQKM6PFSPubtFhwnRJNvWTdJMIffqz5A8/fhQEfOzKC/dGaVdLhC+9dSffuu15tMZt\nPvidJ3nR39zH1x88qkbyCoXinFkx4n7ovx4hYw9iR8dIRrdOqx/OOTwiHVbn4HmrExfcv+etb+G7\n/88NfPmtO2mN23z4P/fy/E/+jE/98GlODGcvuD8KheLiZsUshXz6oceItfUBkFz98mn1H919BGlq\nvHdt24V2rYKmCV5yeTs3XdbGQ4eH+fIvDvP5nx/icz8/xAs2tfKGnZ28+NI2ItaK+WtTKBTnyIpR\niUy2QKKpHzyId7+2ps73fe4amyAC/N6Wue3+uJAIIbhuQwvXbWihdzTHtx45wbcfOc67v/EYYVPn\npsvaeOW2NbxwSyshU19sdxUKxRJkxYh73oR0YoCkE0GzYjV1/7jnJIWwztsisXPeamChWJMM8/6X\nbuaPb9rEr48McfeTp/jhntPc/eQpbEPjeetbeOGWVm7c3Mq6dBQhxGK7rFAolgArQtx93ydrZ4nG\nhkmKq6bVf+nEAJoBH7qh+8I7N0d0TXD9hjTXb0jz8Vdv5cHDQ/zs6X5+fmCA//W9fQB0psI8Z10z\n13Y3c213ig2tMSX2CsUKZUWI+zM/exgvcQIhJMnWG2vq7j85wmBE4wW+Sewimcs2dI0XbGrlBZta\n4VVwfCjLz5/p55cHB/n5gQHu3H0SgFTEZEdXiq0dCbauaWLrmiY6kmEl+ArFCuDiULPz5Olf7qKp\nawApIbGudrOwv97bA7rkY1dfvAdqdLVEeMt13bzlum6klBwdyvLI0WF2HR3mseOj3Hugn/KJf8mI\nydY1TWxqi7OhLcaG1igbW2O0xm0l+grFMmJO4i6EuBn4B0AHviyl/OsZ7F4P3AFcK6XcNW9enicT\nkwXiTQOYOQsjNnVoxmjB4TEc1hQEW1tiZ+jh4kEIwbp0lHXpKG/YuRaAXNHj6dPj7OkdZ1/vGHt7\nx7lj1wkmi1Pr6OO2wfrWKJe0ROlMhelMRehMhelIhelIhtWDW4XiImNWcRdC6MBngZcCPcAjQoi7\npJT76uziwHuBXy+Eo+dD3pC0xQdple015X/9+DGkqfHONelF8uzCELZ0tnel2N6VqpRJKekbL3Bo\nIBOE/gyHBiZ5/MQo9zx1CrfucO/WuM2aRIjWeIi2JpvWmE1bk01bPERb3Ka1FMx52EVToVCcP3MZ\nuT8HOCilPAwghPgW8BpgX53dXwKfAj4wrx7OA8XEAIZZJBnZWVN+5+AYlgbvvGzNInm2eAghWJUI\nsSoR4oaNtTc3z5f0jec5OZqjZyRLz3COnpEcp8bz9Ixkeez4CEOTxYb9xm2DZNQkFbFIRiySYZNU\nxCQZsUhFTFJRi0TYJB4yiNkmsZBBzA6CrqlpIYVivpiLuHcAJ6ryPUDNEUZCiO3AWinl3UKIJSXu\nB+57GC0RHHGX6HxFpfyeo4OMR3Ru0WwMNdqsQdcEa5Jh1iTDXNvd3NDG8XwGMwX6xwsMTBTonwji\nkWyR0WyRkazDaLbI0cFJRrJFJvLurNeNWHpF6KtFP2YbhCydsBmEkKkRMnXC1WVWua66TMM2dCxd\nwzI0dfNQrCjmIu6N/kdUfrMLITTg08DbZu1IiNuA2wC6urrm5uF5svfnD9O0aQC/aBJZ88JK+d88\nfRIMyYe3X7wPUhcTU9dYnQizOhGek73r+YzlHEayDmO5IpmCRybvkik4TORdMgW3lHeZKKUnCy7H\nJ7NM5F0Krkeu6JFzPPxz3AJfE2AZWkXsy7FZF9vldE25QNcEhqaVYoGhC3RNw9BEpawS60H5me2C\nvkxdoGkCXQg0IRACNCHQtFIsgl9aNfVaUF5jX7IRVe3q68t9KZY/cxH3HmBtVb4T6K3Kx4ErgPtK\n/2hWAXcJIV5d/1BVSvlF4IsAO3fuvCCnVGQyeVLxQcL5aOVwjv5sgf26R3dRY31SnYZ0ITB0jZaY\nTUvMPq9+pJQ4niTneOSdKcHPO1Vx0SdXyhddvxIcz6folfJeVZlbW5YpuJVyx5OVOt+XOJ6P50tc\nX1bii40psQ+EXpTKBKIUT5VTnRdTI73qdtS0m95P+Zrle8qs16nrh+p+qvqq7oeadrX25euXr11d\nUV9fbTN17cblU9/nTNeauf7mK1bz+mumFncsBHMR90eATUKIdcBJ4I3A75YrpZRjQGXSVghxH/CB\npbJaphgqEImO0aFNHan3ySeOIw2NP1zEfWQU54YQAssQWIZGIjz/h6mcLVJKfAmuXyX6nsQp5z1Z\ndzOovTlU3yyklPg++KU+y30H+VIo1Us5ZRfkg/6m8o3qp9K+rL2OLH0WKSmlQRLky5+zXA5TdVNl\npbycXieRld/6sq6v6jzVbev6of46M/RDnY81+ZnKq9JlRxrZlL+H2nxtPTPW1/Y3kXdYaGYVdyml\nK4R4N/AjgqWQX5FS7hVCfBzYJaW8a6GdPC+agx8ZibYXVYruHh7H0uCtmxd/HxnFxU0wXQK6ppaK\nKpYWc1rnLqW8B7inruwjM9i+8Pzdmh8O/vpx9ORJpBQ0db8OgAd7RxmL6LwIa8ntI6NQKBTzxbJW\ntz0/eYB40wDueBNGJJiC+dS+HpCSP7tiYee7FAqFYjFZ1uI+PpYhHh8k7gcv77iezyNugZac5OrW\npkX2TqFQKBaOZb23jEyMYBgu6zqCl5e+tL8X19Z5XVQJu0KhWN4s65E7LcHD1NQlrwbgq8cHEK7P\nB7atPVMrhUKhuOhZtuLes+cAZqIXrxAi1H4dJzN5jlmwxdNJhhZ/CZ1CoVAsJMtW3B/93r00JQZw\nR1sQmsbfPXUCdME71rXP3lihUCgucpatuE9MDBAOT5AwAjH/4fAERt7jdzepF5cUCsXyZ9mKO839\nAFy29eUcGJlkKCTYYai17QqFYmWwbJVOS50uvbz0aj69pwc0wR9uWr3YbikUCsUFYVmKe98zRzGb\n+nAmkhiRNn42nsHOebyie3kfyqFQKBRllqW4P3LXj4g3DeKNNbOrf5zxiM51odBiu6VQKBQXjGUp\n7uMTJzDNAgl7FX+/Lzio4z2XrrzTlhQKxcplWYq7SA4AcOXO1/CrbI5ozuOGNalZWikUCsXyYVmK\nu9bch+/p7LGvJxfWeUFEHcihUChWFstub5nBnpNYTacpjqf5+qFRAP5wi5qSUSgUK4tlN3J/+Pa7\nicaGcceaeXAySzjn8bzVicV2S6FQKC4oy07cxycPoeseRXsjExGda221SkahUKw8lp24lx+m/rL1\nFgDevlEdpadQKFYey27OXU/14RZD/MRLY3keL12rVskoFIqVx7IauU/0D2E19TEy1sVwSGOb2ktG\noVCsUJaV8j14+7eJREd5zHsOaIK3dKsdIBUKxcpkWYn76MQzCAG7w1eiFz1+a33rYrukUCgUi8Kc\nxF0IcbMQ4oAQ4qAQ4s8b1L9fCLFPCPGkEOKnQohL5t/VOfiZHMBDY2+sm/W+jqEvq3uXQqFQzJlZ\n1U8IoQOfBW4BLgduFUJcXmf2GLBTSrkN+A7wqfl2dC4YqX72FbbjGwYvb1Nr2xUKxcplLkPb5wAH\npZSHpZRF4FvAa6oNpJT3SimzpexDQOf8ujk7k2Pj2E19POZcB77kD7aovdsVCsXKZS7i3gGcqMr3\nlMpm4u3AD87HqXPhwW99FTs0yVPmVpJ5nzUx9fKSQqFYucxlnbtoUCYbGgrxZmAncOMM9bcBtwF0\ndXXN0cW5MTr+ND4t9IbaeJmw57VvhUKhuNiYy8i9B1hble8EeuuNhBAvAf4n8GopZaFRR1LKL0op\nd0opd7a2zu9KFpEY5HG5A4A3r1dLIBUKxcpmLuL+CLBJCLFOCGEBbwTuqjYQQmwHvkAg7P3z7+bs\n6IkBdvvPxSh4vKRTvZWqUChWNrOKu5TSBd4N/AjYD9wupdwrhPi4EOLVJbP/DcSAO4QQjwsh7pqh\nuwUhNzmJ3jTMfu1yNmGot1IVCsWKZ057y0gp7wHuqSv7SFX6JfPs11nx629/lRPdnRSEzU2tTYvp\nikKhUCwJlsUQd3B4D3u4EqTkzZvULpAKhUKxLMRdxIfZyzbiuQLdTeHFdkehUCgWnWUh7m5ynINs\nZrMVXWxXFAqFYklw0Yt7MZvnWDKFJwxubk8utjsKhUKxJLjoxf3hO7/OAXMzuu9x66b2xXZHoVAo\nlgQXvbj3nX6MvVzJ2kwf6bC12O4oFArFkuCiF/fJpizHWEeXHl9sVxQKhWLJcNGL++HWJqTQuGXD\nhsV2RaFQKJYMF7W4F4tFjv7f9u4+Ro66juP4+0NbLmmB9vog1sr1wYJJDQmtxeBD0QQttCoVn1I0\noVESQiKJxJgAISHE/9CoiZFIMDQCAamCjdVAoFEjJkq1rS1t5aHX59KzrS3Qkmrh6Nc/5ncwvd5e\nd6+7MzvbzyvZ7NzvZnc+95vZ783+5vZ+4y5kdPRznWddMjN7R6WL+7qVj9I7ejY9x/qY0DWm7Dhm\nZm2j0sV919717GAW048cLjuKmVlbqXRx3zm5i7c1hg9NmVF2FDOztlLp4r6jO/vXvl+7/PKSk5iZ\ntZfKFvf+/n62j53G1Lf2M2vC2LLjmJm1lcoW93W/+xXbzpnNrDf6yo5iZtZ2Klvcn3tlC8c0jpmv\nvVZ2FDOztlPZ4r5jUvaJ1E9eMrfkJGZm7aeyxX3P+G7Oi6N89mMLyo5iZtZ2Klvcd3VNpef4Ps+X\namY2hEpWxr/8/nH2aRo9Rw+VHcXMrC1Vsrj/dc+/OKFR9Bx5o+woZmZtqZLFfXd39nftC2ZfVnIS\nM7P2VFdxl3SNpJck9Uq6fYjvd0lakb6/RtKMZgfN233BRLpPHOYzn1rYys2YmVXWaYu7pFHAvcAi\nYJUkYx8AAAZWSURBVA5wvaQ5g1a7EXg1ImYDPwbuaXbQvN1d76PnzX2t3ISZWaXVc+b+EaA3IrZH\nxJvAY8CSQessAR5My48DV0lS82K+69nVq9h/znvpOeKLqWZmtdRT3KcBe3Jf701tQ64TEf3A68Ck\nZgQc7NntmwDoee1YK57ezKwj1FPchzoDjxGsg6SbJK2VtPbgwYP15DvFmDiHS49v4crZl47o8WZm\nZ4PRdayzF7go9/X7gcED3gPr7JU0GhgPnDKDRkTcD9wPMH/+/FOKfz1uu/kObhvJA83MziL1nLn/\nA7hY0kxJ5wJLgVWD1lkFLEvLXwb+GBEjKt5mZnbmTnvmHhH9km4BngZGAcsjYouk7wFrI2IV8ADw\nsKResjP2pa0MbWZmw6tnWIaIeBJ4clDbXbnl/wFfaW40MzMbqUp+QtXMzIbn4m5m1oFc3M3MOpCL\nu5lZB3JxNzPrQCrrz9ElHQR2jfDhk4H/NDFOs7RrLmjfbM7VGOdqTCfmmh4RU063UmnF/UxIWhsR\n88vOMVi75oL2zeZcjXGuxpzNuTwsY2bWgVzczcw6UFWL+/1lB6ihXXNB+2ZzrsY4V2PO2lyVHHM3\nM7PhVfXM3czMhlG54n66ybpbvO2LJP1J0guStkj6dmq/W9Irkjak2+LcY+5IWV+SdHULs+2UtClt\nf21qmyhptaSt6b47tUvST1Ku5yXNa1GmD+b6ZIOkI5JuLaO/JC2XdEDS5lxbw/0jaVlaf6ukZUNt\nqwm5fiDpxbTtlZImpPYZkv6b67f7co/5cNr/vSn7GU1zWSNXw/ut2a/XGrlW5DLtlLQhtRfZX7Vq\nQ3nHWERU5kb2L4e3AbOAc4GNwJwCtz8VmJeWzwdeJps0/G7gu0OsPydl7AJmpuyjWpRtJzB5UNv3\ngdvT8u3APWl5MfAU2QxaVwBrCtp3/waml9FfwJXAPGDzSPsHmAhsT/fdabm7BbkWAqPT8j25XDPy\n6w16nr8DH02ZnwIWtSBXQ/utFa/XoXIN+v4PgbtK6K9ataG0Y6xqZ+71TNbdMhHRFxHr0/JR4AVO\nnU82bwnwWEQcj4gdQC/Zz1CU/MTlDwJfyLU/FJnngAmSprY4y1XAtogY7oNrLeuviHiWU2cHa7R/\nrgZWR8ThiHgVWA1c0+xcEfFMZHMRAzxHNvtZTSnbBRHxt8gqxEO5n6VpuYZRa781/fU6XK509v1V\n4JfDPUeL+qtWbSjtGKtaca9nsu5CSJoBzAXWpKZb0tur5QNvvSg2bwDPSFon6abUdmFE9EF28AHv\nKSHXgKWc/KIru7+g8f4po9++SXaGN2CmpH9K+rOkBaltWspSRK5G9lvR/bUA2B8RW3NthffXoNpQ\n2jFWteJe10TcLQ8hnQc8AdwaEUeAnwEfAC4D+sjeGkKxeT8eEfOARcC3JF05zLqF9qOy6RmvBX6d\nmtqhv4ZTK0fR/XYn0A88kpr6gJ6ImAt8B3hU0gUF5mp0vxW9P6/n5BOIwvtriNpQc9UaGZqWrWrF\nvZ7JultK0hiynfdIRPwGICL2R8TbEXEC+DnvDiUUljci9qX7A8DKlGH/wHBLuj9QdK5kEbA+Ivan\njKX3V9Jo/xSWL11I+xzw9TR0QBr2OJSW15GNZ1+ScuWHblqSawT7rcj+Gg18EViRy1tofw1VGyjx\nGKtaca9nsu6WSWN6DwAvRMSPcu358errgIEr+auApZK6JM0ELia7kNPsXOMknT+wTHZBbjMnT1y+\nDPhtLtcN6Yr9FcDrA28dW+SkM6qy+yun0f55GlgoqTsNSSxMbU0l6RrgNuDaiDiWa58iaVRankXW\nP9tTtqOSrkjH6A25n6WZuRrdb0W+Xj8NvBgR7wy3FNlftWoDZR5jZ3KFuIwb2VXml8l+C99Z8LY/\nQfYW6XlgQ7otBh4GNqX2VcDU3GPuTFlf4gyvyA+TaxbZXyJsBLYM9AswCfgDsDXdT0ztAu5NuTYB\n81vYZ2OBQ8D4XFvh/UX2y6UPeIvs7OjGkfQP2Rh4b7p9o0W5esnGXQeOsfvSul9K+3cjsB74fO55\n5pMV223AT0kfUGxyrob3W7Nfr0PlSu2/AG4etG6R/VWrNpR2jPkTqmZmHahqwzJmZlYHF3czsw7k\n4m5m1oFc3M3MOpCLu5lZB3JxNzPrQC7uZmYdyMXdzKwD/R+JtAev0iccowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fd8bd3b3470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 10\n",
    "A, q, B, d = generate_problem(n)\n",
    "q = -q\n",
    "\n",
    "x = np.array([0] * n)\n",
    "y = np.array([0] * n)\n",
    "K = 20\n",
    "gamma = 0.005\n",
    "n_iter = 2000\n",
    "x_hist = [[] for _ in range(n)]\n",
    "for iter_n in range(n_iter):\n",
    "    for i in range(n):\n",
    "        x_hist[i].append(x[i])\n",
    "    x, y = step(x, y)\n",
    "args = []\n",
    "for a in x_hist:\n",
    "    args.append(range(n_iter))\n",
    "    args.append(a)\n",
    "plt.plot(*args)\n",
    "print('f(x) = ' + str(x.dot(A).dot(x) - q.dot(x)))\n",
    "print('Bx = ' + str(B.dot(x)))\n",
    "print('x = ' + str(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Результат, выданный Mathematica:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 315,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-8.7472161488419999"
      ]
     },
     "execution_count": 315,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1 = np.array([0.37777, 0.688865, 0.933308, 1.11108, 1.2222, 1.26667, 1.2445, 1.1556 ,1. ,0.75002])\n",
    "x1.dot(A).dot(x1) - q.dot(x1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Эксперимент 2.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Реализовать метод барьеров и метод одновременного решения прямой и двойственной задачи для решения задачи квадратичного программирования (без ограничений типа равенства, для простоты):\n",
    "$$\\min_{Ax\\preceq b} \\ \\frac{1}{2} x^TPx + q^Tx,$$\n",
    "где $A\\in\\mathbb R^{m\\times n}$. Можете предполагать, что для начальной точки все неравенства строгие. Проверьте код на несольких примерах. Изобразите графики основных показателей работы алгоритмов."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В методе барьеров без ограничений типа равенство шаг ищется по формуле\n",
    "$$(t\\nabla^2 f_0(x) + \\nabla^2\\phi(x)) \\Delta x = - t\\nabla f_0(x) - \\nabla\\phi(x),$$\n",
    "где (предполагаем $P$ симметричной)\n",
    "$$f_0(x) = \\frac{1}{2} x^TPx + q^Tx,$$\n",
    "$$\\nabla f_0(x) = P x + q,$$\n",
    "$$\\nabla^2 f_0(x) = P$$\n",
    "$$f_i(x) = A_i x - b_i$$\n",
    "$$\\nabla f_i(x) = A_i^T$$\n",
    "$$\\nabla^2 f_i(x) = 0$$\n",
    "$$\\phi(x) = -\\sum_{i=1}^n \\log(-f_i(x)) = -\\sum_{i=1}^n \\log(-A_ix+b_i)$$\n",
    "$$\\nabla\\phi(x) = \\sum_{i=1}^n\\frac{1}{-f_i(x)}\\nabla f_i(x) = \\sum_{i=1}^n \\frac{A_i^T}{-A_ix + b_i}$$\n",
    "$$\\nabla^2\\phi(x) = \\sum_{i=1}^n\\frac{1}{f_i(x)^2}\\nabla f_i(x) \\nabla f_i(x)^T + \\sum_{i=1}^n\\frac{1}{-f_i(x)}\\nabla^2f_i(x) =\n",
    "\\sum_{i=1}^n \\frac{A_i^T A_i}{(A_ix - b_i)^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 590,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def barrier_step(x, t, eps):\n",
    "    while True:\n",
    "        phi2 = sum([np.outer(A[i], A[i])/(A[i].dot(x) - b[i])**2 for i in range(n)])\n",
    "        lhs = t * P + phi2\n",
    "        phi1 = sum([-A[i].T / (A[i].dot(x) - b[i]) for i in range(n)])\n",
    "        df = P.dot(x) + q\n",
    "        rhs = - t * df - phi1\n",
    "        l = rhs.dot(np.linalg.inv(lhs)).dot(rhs)\n",
    "        f = x.dot(P).dot(x) / 2 + q.dot(x)\n",
    "        gaps.append(n / t)\n",
    "        if l < 2 * eps:\n",
    "            break\n",
    "        x += np.linalg.inv(lhs).dot(rhs)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Протестируем на матрицах разного размера из первого задания."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 603,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n = 2\n",
      "[ 0.49999947  0.49999962] -0.5\n",
      "n = 3\n",
      "[ 0.74999832  0.99999817  0.74999893] -1.25\n",
      "n = 10\n",
      "[ 0.37777712  0.68888792  0.93333234  1.11111032  1.22222183  1.26666681\n",
      "  1.24444518  1.15555685  1.00000152  0.75000067] -8.74722204135\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd8beea60f0>"
      ]
     },
     "execution_count": 603,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVOfZ+PHvMww7CLIJsggoiDvuK6AY10TNapKmjU2T\nmrZJmjTNYpq3afN737fN1iRN2jfGJE2s2YxNUs2iJoiK+76LuC8gIu4r+/P7Y7ABB2QQzpyBuT/X\ndS6Gc84z5/Yw3hyeVWmtEUII0fpZzA5ACCGEc0jCF0IINyEJXwgh3IQkfCGEcBOS8IUQwk1IwhdC\nCDchCV8IIdyEJHwhhHATkvCFEMJNWM0OoKawsDAdHx9vdhhCCNGibNiw4YTWOryh81wq4cfHx7N+\n/XqzwxBCiBZFKXXIkfOkSkcIIdyEJHwhhHATkvCFEMJNGFqHr5SKBf4JRAJVwAyt9V+NvKYQQlxL\neXk5+fn5lJSUmB1Ko/n4+BATE4Onp+d1lTe60bYC+K3WeqNSKhDYoJT6Xmu90+DrCiFEnfLz8wkM\nDCQ+Ph6llNnhOExrzcmTJ8nPzychIeG63sPQKh2tdaHWemP16/NALhBt5DWFEOJaSkpKCA0NbVHJ\nHkApRWhoaJP+MnFaHb5SKh7oDaxx1jWFEKIuLS3ZX9HUuJ3SD18pFQB8DjymtT531bGpwFSAuLi4\n63r/ootFzNk955rnDIoaRL/Iftf1/kII0RoYnvCVUp7Ykv1HWusvrj6utZ4BzADo16/fdS2wW3y5\nmBlbZ9R7XKP59sC3fHPLNy32N7sQonU4cuQI9957L8eOHcNisTB16lQeffRRp1zb6F46CngPyNVa\nv2rUdbqHdWfrlK31Hp+9azb/s+Z/OHjuIAlB19fYIYQQzcFqtfKXv/yFPn36cP78efr27cuoUaPo\n2rWr4dc2ug5/KPATIFMptbl6G2/wNe2kx6QDkJOf4+xLCyFELVFRUfTp0weAwMBAunTpQkFBgVOu\nbegTvtZ6OWB6HUpUQBRJbZPIyc9hSrcpZocjhHARz3+1g51HzzV8YiN0bd+GP0zo5tC5Bw8eZNOm\nTQwcOLBZY6iPS02eZqT06HRm7pjJzB0zsajaf9j4e/ozqeMkPCweJkUnhHA3Fy5c4LbbbuP111+n\nTZs2Trmm2yT80fGj+WDHB7yy/pU6j4f4hDA8drhzgxJCmMrRJ/HmVl5ezm233cY999zDrbfe6rTr\nuk3C7xralVU/WkV5VXmt/RVVFYz7fBw5+TmS8IUQhtNac//999OlSxcef/xxp17brSZP87X60sar\nTa0txCeEwe0Hk5Ofg9bX1StUCCEctmLFCmbNmkV2djapqamkpqby7bffOuXabvOEfy0ZMRksOryI\n3ad30zmks9nhCCFasWHDhpn2cCkJH0iLSQPggx0fMChqkN3xfpH9iA6QKYCEEC2bJHwgzDeMPhF9\n+Hr/13y9/2u74wOjBvLu6HdNiEwIIZqPJPxqM0bPoPhSsd3+f+78J3Py5nCh7AIBXgEmRCaEEM3D\nrRptr8Xbw5uYwBi7bUz8GCp0BSuPrjQ7RCGEaBJJ+A3oFd6LNl5tZFoGIUSL1zqqdM4fg62zf/je\n6gt9p4DVu8lvbbVYGRo9lGUFy1hyZInd8WDvYFIjUpt8HSGEMFrrSPjnCuD752rv8w6A1B81y9uP\njBvJ/APzeST7kTqPfzHxC5LaJjXLtYQQrVtJSQnp6emUlpZSUVHB7bffzvPPP++Ua7eOhB+VCr87\nanutNbzZF3YvbLaEP7rDaL6Y+AVlVWW19l8ou8AD3z3A0vylkvCFEA7x9vYmOzubgIAAysvLGTZs\nGOPGjWPQIPsu4c2tdSR8iwd4+f/wfdIo2DkXKsvB4/pWd69JKVVvQu8S0oWc/Bwe6PFAk68jhGj9\nlFIEBNh6/JWXl1NeXu60hZlaR8K/WvJY2DQLDq+ChHRDL5URm8GMrTM4U3KGYJ9gQ68lhGhm86fB\nsW3N+56RPWDcC9c8pbKykr59+7J3714eeughmR65MU5eKCUrt+g/3/urZG708EKtng6nD9p2xgyA\niJRmv3Z6dDrTt0zns92f0T+yf61jCkVy22T8PP2a/bpCiJbLw8ODzZs3c+bMGW655Ra2b99O9+7d\nDb+uM9a0HQv8FfAA3tVaX/tX33XIP32Zpz+v/Vu6b/wwovK+gbxvbDsiusGvmr8vfbewbkT4RvDm\npjfrPH5b0m38ccgfm/26Qohm0MCTuNGCg4MZPnw4CxYsaPkJXynlAfwdGAXkA+uUUvO01jub8zop\nUYGsnJYJgAYmvLmcl4J+x2t3VyfhTR/Ckj/DmcMQHNecl8aiLMwcN5PD5w/bHfvnjn+y+MhintPP\n2S26IoRwT8XFxXh6ehIcHMzly5fJysri6aefdsq1jX7CHwDs1VrvB1BKfQpMApo14XtbPWgf7Puf\n74d3Did713EqAvpi9bBA99tsCX/3Qhjw8+a8NMB/RuVe7eTlk6w4uoKdJ3fSPcz4395CCNdXWFjI\nlClTqKyspKqqismTJ3PTTTc55dpGJ/xo4EiN7/MBw1snMlMi+GJjAZuOnKF/fAiEdoK2CbDnO0MS\nfn2GRQ9DocjJz5GEL4QAoGfPnmzatMmUaxud8Ovqa1RrImil1FRgKkBcXPNUt6QlhWO1KGatOsTJ\nC6UADO0wksDtH8KOf4Oy2LpxdswEA7tDtfVpS8/wnmQfzmZM/Bi740HeQYT5hhl2fSGEqMnohJ8P\nxNb4PgY4WvMErfUMYAZAv379mmVVgCBfT4Z0CmPelqPM22K73F0RSbxQUQJzpvxw4o+/gE4jm+OS\n9RoeO5y/bvwrN8+92e6Yj4cP39/+vXTnFEI4hdEJfx2QpJRKAAqAu4DmGf7agLfu6cPhU5cA+HJT\nATNy4Le/XE24dyXoSnhvDOxeYHjCv6fLPcS3iadCV9TaX3SxiFfWv8Lyo8u5KdE59XdCCPdmaMLX\nWlcopR4GFmLrlvkPrfUOI695hb+3lS5RbQCwKMWMnP1kFflz94DqaqPEDFvCH/eSodU6vlZfbuhw\ng93+Kl3FP7b/g5z8HEn4QginMLyvoNb6W611sta6o9b6f42+Xl2S2wUQHexL9q7jNXaOsXXTLM4z\nIyQsysKw6GGsKFhBRVVFwwWEEKKJWsVI24YopchMieDzjfms3HsCpRTBAYPoArDhfUipfsJu1w38\nQpwWV3pMOvP2zWNFwQq6hXWr8xx/T398rb51HhNCiMZwi4QPMLpbO2atPsSP3l3zn33bo7sTsGY6\nrJlu29ExE37ypdNiGtJ+CFaLlYezH673nLbebfnu9u/wsfo4LS4hhLF+9rOf8fXXXxMREcH27dsB\nOHXqFHfeeScHDx4kPj6ezz77jLZt2zbrdZXWzdIxpln069dPr1+/3pD31lqzJf8sl8sqqdKa+z5Y\nxy/6BvJ4avW/f/PHsO0zeOoA+LQxJIa6rClcw6Fzh+o8ln8+n/d3vM//jfw/0mLSnBaTEK1Zbm4u\nXbp0MTWGnJwcAgICuPfee/+T8J966ilCQkKYNm0aL7zwAqdPn+bFF1+0K1tX/EqpDVrrfg1d122e\n8JVSpMb+0P1xUGIoX++7xOO3DK8+wQJbPob9i6HrJKfFNTBqIAOj6h6LVlpZyqd5n5KTnyMJX4hW\nJD09nYMHD9baN3fuXJYsWQLAlClTGD58eJ0JvyncJuFfLbNzOH/8aicHTlwkIcwfYgeCT5Bt+gUn\nJvxr8fbwZmDUQHLyc/id/p3T5swWwl28uPZFdp3a1azvmRKSwtMDGj83TlFREVFRUQBERUVx/Pjx\nBko0ntvO6JWZ0g6Af204ws6j59hZdImKxJG26RcKt9rmyL4ytbKJ0mPSOXrxKHmn8yirLKtzq9JV\nZocphGgB3PYJPy7Uj87tAvn74n38ffE+AH4X24WpF7+At2tUn/xypa33jknSo20LuNzx1R31nhPf\nJp65N8+VGTmFaKTreRI3Srt27SgsLCQqKorCwkIiIiKa/Rpum/AB3p3Sjx1HzwHw1ZajvJoL994z\nGx/KoaIUPn8Adn1rasJv59+OlzNeJv98fp3HD549yNx9c9l+Yjs9w3s6OTohRHOZOHEiM2fOZNq0\nacycOZNJk5q/atmtE35siB+xIbbVqAJ9rHyzrZBlujejutqqe1j9f7bRuBlPmhgljI0fW++xs6Vn\n+Wr/VyzNXyoJX4gW4u6772bJkiWcOHGCmJgYnn/+eaZNm8bkyZN57733iIuLY86cOc1+XbdO+DX1\njw8hwNtK9q6iHxJ+0hjbPPoXiiEg3NwA6xHkHURqeCrL8pfxSO9HzA5HCOGATz75pM79ixYtMvS6\nUulbzctqIS0pjOxdxyk4c5mjZy5zPi4T0LD9czibb9uqKs0O1U56TDq5p3I5fqn5W/WFEK2HPOHX\nkJkSwfztxxj6QjYAHhbIaxuJdcHTsKC6caf3T2DS30yM0l56TDqvb3ydcZ+Pw8PiUec5PcJ68N6Y\n95wcmRDClUjCr2FSajSeHhZKKyoprajiubk7+DLlVe6IPmk7Ydu/IPcruOl18HCdW9cpuBPTBkzj\n2MVjdR4/cPYAS/OXcvDsQeKD4p0bnBAuSGvdIse1NHVmBNfJWi7Ay2rh5t7R//l+9rojzCmwcseE\nG207vNvYFlDJXwcdBpsUpT2lFPd0uafe4wUXCliav5Sc/BxJ+MLt+fj4cPLkSUJDQ1tU0tdac/Lk\nSXx8rn9eLUn415CZEsHfF+/lzKUygv28oOMIsFhtPXdcKOE3JDogmk7BncjJz+HebveaHY4QpoqJ\niSE/P5/i4mKzQ2k0Hx8fYmJirru8JPxrGJESwZvZe8nKPc74HpFg8cc3bjBq90LIeMp2kocXeHia\nG6gD0mLSmLVjFhfKLhDgFWB2OEKYxtPTk4SEBLPDMIVhs2UqpV4GJgBlwD7gPq31mWuVMXK2zOtR\nWaUZ8L9ZnLxY9p99b3dcxZiCN384ySsAHtkAgZEmROi49cfWc9/C+wjxCcFqqfv3/KCoQfzvMFPW\nqBFCNIGjs2UamfBHA9nVyxy+CKC1vuY4ZldL+ABr9p9k8xHb76nvdhZxvPgES0flY6kqh9LzkPOy\nrRG3330mR3ptlVWV/H3z3zlZcrLO4wfOHmDz8c1kT84mzDfMydEJIZrC9OmRtdbf1fh2NXC7Udcy\n0sDEUAYmhgIQGeTDo5+eZkv03fSOawtaw5bZtgnXXDzhe1g8+HWfX9d7PPdkLpO/nszyguXc3Olm\nJ0YmhHAWZw28+hkw30nXMkxGcjgWxQ9r4yplWxt3/xIoLzE1tqZKCUkhwjeCnPwcs0MRQhikSQlf\nKZWllNpexzapxjnPAhXAR/W8x1Sl1Hql1HpXbzUP9vOib4e29ouhl1+Cg8vNC6wZKKVIi0lj1dFV\nlFeWmx2OEMIAhi5xqJSaAvwCGKm1vtTQ+a5Yh3+1t5bs48UFu/D0UPxqeCd+MzwOXkqAihJQ1aNc\nRzwDab81N9DrkH04m0cXP0p8m3g8r+p55GXx4k9pfyIxKNGk6IQQ9XG0Dt+wKh2l1FjgaWCiI8m+\npbh7QCyPjkyiU0Qgn60/grZ6w60zYOhjMOQRCO0Em+r8Y8blDY0eyq1Jt9IxuCNxgXG1tl2ndvH1\nvq/NDlEI0QRG9tLZC3gDV7qFrNZa/+JaZVrCE/4Vs9cd5unPtzH/0TS6RNVY9HztO/DtE/DwBgjr\nZF6AzeynC37KhbIL/Gviv8wORQhxFdOf8LXWnbTWsVrr1Ortmsm+pRnR2bYaTa36fICk0bavexY6\nOSJjZcRkkHc6r975eoQQrk+mR75OEW186BEdZJ/w23aA8C62xdBbkfQY21KL0otHiJZLplZoghEp\nEbyxaA/9/ieLXw7vyP3DqodrJ4+GFW/Ay0m27wPbwc8Wgpe/ecE2UWJQItEB0UzfMp1vD3xrd3xE\n7AimdJtiQmRCCEfJE34T/HhgHD8dEk+gj5UPVx/64cCAqTDg55ByI3QYAse22frqt2BKKX7Z65fE\nB8VjUZZaW+GFQqZvmU55lXTnFMKVGdots7FaUqNtTTNXHuQP83aw+InhJIRd9RRfUQYvJUL3W2Di\nm3W/QQu36PAiHlv8GO+Nfo8BUQPMDkcIt2N6o607yUyppwEXwOoFnTJhz/e2qRhaocFRg/G0eEr9\nvhAuThJ+M4gN8SMpIoDsXUV1n5A0Bs4XwrGtzg3MSfw8/ejXrh85BZLwhXBl0mjbTDJTInhn2X5G\n/mUJD2Z0ZHK/2B8OJo2yff1oMvi0ARQMfxq632ZKrEbIiM3ghbUv8MiiR7Co2s8Rvp6+PDPgGYK8\ng0yKTggB8oTfbO4Z2IGJvdpzuayS91ccrH0wIAJG/sG2Sla7blByBtbMMCVOo4zuMJrU8FSOXjxK\n/oX8/2yHzx/mm/3f8N2h7xp+EyGEoaTRtpm9vXQff56/i1XPZBIV5Fv3Sdn/C8tegSf3gV+IcwN0\nMq01474YR1JwEm+ObJ2N1kKYTRptTXLNBtwrkseCroK9WU6KyjxKKdKi01hzbA2llaVmhyOEW5OE\n38w6RQQQG+LL4msl/Pa9wT+81Y3GrU96TDqXKy6z7tg6s0MRwq1Jo20zU0qR2TmCT9Ye4fa3VvJA\nWgJju0fVPsligU6jYMeX8F713DtBsbZZNy0ezg/aYP0j++Pj4cMbG99g4UH7X3KjOoz6z9QNQgjj\nSMI3wD2DOnDo1CW25Z/l7Zz99gkfYOBUuHgcqirg8hnY/i8Y+CDEtr6BSz5WH+7sfCcLDy1kdeHq\nWsfOlp5la/FWSfhCOIE02hro9azd/HXRHtY/ewOhAd71n3j5NLzUEYb9Bkb+3nkBuoAPd37Ii+te\n5NtbvyU2MLbhAkIIO9Jo6wIyUyLQGpbkNbB0o29biB3oNnX6NcksnEI4jyR8A3VvH0R4oDfZeddo\nwL0ieQwUbYOzBcYH5kLi2sQR3yZeEr4QTmB4Hb5S6gngZSBca33C6Ou5EotFMaJzOF9tKeT+D9bx\n06HxpCWF131y8hjI+gPMvgf8I0BZbFU8cQOdG7QJ0mPS+WTXJ7y87mUUqtaxNt5t+Fn3n2G1SHOT\nEE1l6P8ipVQsMAo4bOR1XNmPBnZgd9EF1h44xaWyyvoTfngK9LwTivPgQhGc2AMeVrdI+OMTxzNv\n3zzm7J5Ta7/WmpLKErqFdmNo9FCTohOi9TD6sek14ClgrsHXcVmpscH8+6GhvDB/F+8u28+5knLa\n+Hjan6iUrVvmFV89Btvm2KZXtno5L2ATdAvtxrK7ltntL6koIe3TNHLycyThC9EMDKvDV0pNBAq0\n1luMukZLkpkSQUWVZvkeB2u1ksdA2QU4tMLYwFyYj9WHgVEDWZq/FFfqTSZES9WkhK+UylJKba9j\nmwQ8CzznwHtMVUqtV0qtLy5uoDdLC9YnLpggX08W5TrQgAuQkAFWH9jj3pOOpcekU3ChgANnD5gd\nihAtXpOqdLTWN9S1XynVA0gAtiilAGKAjUqpAVrrY1e9xwxgBtj64TclHldm9bCQkRzOdzuP8fhn\nmh8P6kCfuLb1F/Dyg/g02PoZXDpl2xfaCTKedE7ALuJKt803Nr1Bt9BudsdHxo0kMTjR2WEJ0SIZ\nUoevtd4GRFz5Xil1EOjnbr10rnZX/1i2FZzlm62FnL5Yxvv3NTCqtt/PYOHv4PAqKL8EWz+FHrdD\nSIJzAnYBkf6RDIgcwKLDi1h0eJHd8Q1FG5g+aroJkQnR8khfNyca0imMxU8M54/zdvDJ2sNcLqvE\n1+sac+ekjLdtACf3wZt9bFU8Ax90TsAu4t3R71JRVWG3/9UNrzI7bzaXyi/h5+lnQmRCtCxOGXil\ntY5396f7mjJTIiitqGLV/kbcktCOtiqd3QuMC8xFKaXw9PC02zJiMyivKrebn0cIUTcZaWuCgYkh\n+Hl5ON6Ae0XyWDi4HEovGBNYC9M3oi/+nv4ySlcIB0mVjgm8rR4M6xTGgu3H8LJauKt/HJ0jAxsu\nmDQaVv0N5v7K1oun//3GB+vCPD08GdJ+CEvzlzJ712y742F+YYyMG2lCZEK4Jkn4JrmjXyzrDp7i\nn6sOcexsCW/9uG/DheIGQ2RP2PM97JwLncdDmzqmXnYjo+NH8/2h7/mfNf9T5/G5N88lMUh68QgB\nkvBNM6prOzY9N5pnvtjKV1sKKauowsvaQA2b1Qt+sQyKdsBbQ2wNuH2nOCdgFzU2fiyDowZTXlVe\na//Jyye5/avbWZa/TBK+ENWkDt9kIzpHcKG0gvUHTzleKKIrtIlx+0FZVwR5BxHmG1Zr6xzSmaS2\nSSzNX2p2eEK4DEn4JhvaKQwvD8u1Fz2/mlK2qRf2LYYKWRi8PhkxGWwq2sS5snNmhyKES5AqHZP5\ne1sZ1DGUb7YV0sbXk5tTo4kLdaBPefIYWP8ezH/K1oDb/Vbjg21h0mPSeXfbu7y37T26hHSpfVDB\n4KjBBHkHmROcECaQhO8Cbu0dzeOfbebV73dz8MRFXr0zteFCCem2ap0NH8DGf0LicPALMTjSlqVn\nWE8ifCP4x/Z/1Hl8cvJkfj/YvZaUFO5N1rR1EVprfjN7Mzl7TrDu2RvwsChHCkHBBnh3JNz6DvSc\nbHygLczZ0rOcvHzSbv8r618h73QeWbdnUT3fkxAtlqxp28Iopcjs0o5TF8vYkn/G0ULQvg/4h7vl\neriOCPIOIjE40W4b1WEUxy8dJ+90ntkhCuE0kvBdSEZSOB4WRXZjRuBaLLYBWXuzoNJ+vhlRt7SY\nNEAWTxfuRerwXUiQnyd949ry9dajhAd6M65HJBGBPg0XTBoNmz+C7P+GxAzomGl8sC1cmG8Y3UO7\n8/2h7+ke2t3ueDv/dnQM7mhCZEIYR+rwXcys1Yf4/b+3A/CjgXH86ZYeDRcqOQev94CSM+DhBU/t\nB28Hpmpwc+9sfYc3Nr1R5zEvixdL7lxCoJfcR+H6HK3Dlyd8F/OTQR2Y2Ks9T8zZwuJdx9FaN9yo\n6NMGfpsHB5bCx5Nt/fO7TnROwC3YT7v9lIFRA6nSVbX2Hzh7gOdWPsfKoysZEz/GpOiEaH5Sh++C\ngnw9GdW1HYVnS8gtPO9YIU8fW1WOdxDskQZcR3h6eNIzvCepEam1tgkdJ9DGq43U74tWRxK+ixre\nORyA7F1Fjhfy8IROmbbJ1aqqGj5f1MlqsTIsehjLC5bbPf0L0ZIZmvCVUo8opfKUUjuUUi8Zea3W\nJiLQh54xQXy1pZC5mws4e6m84UIASWPgQhGsetPWR19cl/SYdE6VnGL+gfnsOLGj1pZ7MrfOFbiE\ncHWG1eErpUYAk4CeWutSpVREQ2VEbeO6R/Higl08+ulm7h+WwO9v6tpwoaRRYPWF758DT39bA66n\nAz19RC3DoodhtViZtmxanccf7/s493W/z8lRCdE0hvXSUUp9BszQWmc5WkZ66dRWVaXJP32ZaV9s\npfBsCYufGO5YwYsnbTNp/vsXcM/nkHSDoXG2Vnmn8jh28Zjd/tc3vo6/pz8fjv/QhKiEsOcKI22T\ngTSl1Bql1FKlVH8Dr9UqWSyKuFA/xnaP5MCJi+wvdnBpQ/9Q6Haz7UnfDdfAbS6dQzqTEZtht43q\nMIqtxVs5XXLa7BCFaJQmJXylVJZSansd2yRs1UVtgUHAk8Bnqo7+hUqpqUqp9Uqp9cXFxU0Jp9Ua\n0dlWG9aoKZQ9fW0Tqu1ZaJtzRzSb9Jh0NJrlBcvNDkWIRmlSHb7Wut66AqXUL4EvtK3OaK1SqgoI\nA2plda31DGAG2Kp0mhJPaxUb4kdyuwC+2lpIp4gABiaE4uvl0XDB5NGwe75tFG7sQAhLMj5YN9A1\ntCuhPqF8f+h7UsPtZzYN8Q3B39PfhMiEuDYjB179G8gEliilkgEv4ISB12vVxnSL5M3svfz0/XU8\nktmJ347u3HChpDFgscLch8C3LTyxx9Z1UzSJRVlIj0nny71fsvjIYrvjUf5RLLhtARYlvZ6FazGy\n0dYL+AeQCpQBT2its69VRhpt61dWUUXesfM8N287l8sqWfBYumMFTx2AvPmw8Bn46TcQP8zYQN3E\nqZJTrChYYbc/91Qus3bO4uPxH9Mj3IFpMYRoBqZPraC1LgN+bNT7uxsvq4UeMUGM6x7Jn77dRcGZ\ny0QH+zZcMCQB+vzE1k1z90JJ+M0kxCeECR0n2O1Pi07jo9yPyCnIkYQvXI78zdnCZKbYGnAXN6YB\n1zsQ4ofKnPlOEOwTTK/wXiw9IounC9cjCb+F6RgeQFyIH19vPcrmI2cor3Rw6H/SGDiRZ6veOVdo\nbJBuLj0mndxTuRw4e4CzpWdrbRfLL5odnnBjMj1yC/TfX+/kveUHAHhqbGd+NbxTw4VO7Yc3+gAa\nAtvDb3bYFk8RzW736d3cNu+2eo//dcRfyYyTNQtE83G0Dl8Sfgt0sbSCjYdP87/f5OLj6cG/Hxrq\nWMGjmyHvW1j6IjyQDTF9jQ3UjS04sIATl+07pc3YOoNB7QfxUrpMLSWaj+mNtsI4/t5W0pLCGd/j\nDK9l7ebkhVJCA7wbLtg+FYLjIOdl24AsSfiGGZswts79u07tIvtINhVVFVgt8t9POJf8Td+CZaZE\noDUsyWvECGW/EIjpLw24JkmPSed82Xm2FG8xOxThhiTht2Dd2rchItCbb7cVcvDERaqqHKyeSx4D\nhZvhyDrb8ojCaYa0H4JVWVmaL714hPNJHX4L98wX2/hk7WEAnp/YjSlD4hsuVLQT3hpsex3aCR5e\nDw0toyiazQMLH2DNsTV1HhsaPZTpN0x3ckSipZM6fDfx1JjODO4Yymvf72b+9kLHEn67rvCTf9u6\naK59G47n2vYJp3iy/5MsOrzIbv+OkzvIyc+h+FIx4X7hJkQmWjtJ+C1cW38vJvZqT27hOd7J2c+5\nknLa+DgwX07HERCeYkv4uxdIwneiziGd6RxiPxdS3qk8cvJzWF6wnFuSbjEhMtHaSR1+KzEyJYKK\nKs2y3Y2Yn65NFET2tC2WIkyX3DaZdn7tpH5fGEYSfivRO64twX6eLNxxjDOXyhwvmDwWjqyBk/ug\nohHlRLO0ytPlAAAa/0lEQVRTSpEek86qo6soq5SfhWh+UqXTSnhYFCM6R/DlpgLmbTnKy7f35I5+\nsQ0XTB4LOS/Bm30gKhUelKdLM6XHpDNn9xwyZmfYTa/s7eHNWze8VWd1kBCOkITfijw9NoU+Hdoy\nfck+vtlW6FjCj+4Dt79vq9bZ8oltCoaQROODFXUaGj2UB3s+yPmy87X2azRz8ubw7YFvJeGL6yYJ\nvxWJDPLhJ4M6sL/4Ah+vOczlssqGV8ZSCrrfCu172xL+7u9g0C+cE7Cw42nx5OHeD9d5bN+ZfeTk\n5/Cbvr9xclSitZA6/FYoMyWC0ooqVu5rRANuSAKEJcui5y4sPSadvWf2cvTCUbNDES2UYQlfKZWq\nlFqtlNpcvUj5AKOuJWobkBCCn5cHWbnHqXR09C1A0mg4tAIun5GFz11QeoxtlbOc/ByTIxEtlZFL\nHH4HvKa1nq+UGg88pbUefq0yMtK2+Tw4az0LdxRhUfB/9/RhbPeohgsdWAYzb7K9TsiAKfOMDVI0\nitaaG7+8kVMlp2jr3dbu+LiEcfy6z69NiEyYzdGRtkZW6WigTfXrIED+DnWiaeO68MToZEL8vfj3\nJgdvfYehMP4V6DweDiyFc/IjcyVKKZ7o9wQjYkeQGpFaa/P19OXjXR9TXlludpjChRnZaPsYsFAp\n9Qq2XyxDDLyWuEpCmD8PZyZRcKaEr7YcpayiCi9rA7/fLRYY8HNb4s/71tZzp+9PnRKvcExmXGad\ni6csObKER7IfYcPxDQyKGmRCZKIlaNITvlIqSym1vY5tEvBL4Dda61jgN8B79bzH1Oo6/vXFxY2Y\n5lc4JDMlggulFaw7eMrxQhFdIChOplBuQQZEDsDL4iVr6YpralLC11rfoLXuXsc2F5gCfFF96hyg\nzkZbrfUMrXU/rXW/8HCZMKq5De0UipfVQnZjFj1XCpJHw/4lUF5iWGyi+fh5+jEgagDLCpaZHYpw\nYUZW6RwFMoAlQCawx8BriXr4eVkZnBjK+ysO8K8N+fztR71JS3LgF2vyWFj3LryUYHt9x/vGByua\nJD0mnT+t+ROTv5psN0rXz9OPl9NfJtQ31KTohCswstH258BflFJbgD8BUw28lriGp8Z25v5hCWit\n+WJjgWOFEkfAiP+C6L6w899w8aSxQYomGxs/llEdRhHmG0aIT8h/tiDvINYdW0fWoSyzQxQmM+wJ\nX2u9HJBFU11At/ZBdGsfxMkLZSzOs/XN97A0sOCJhxUynoROmfBOJuzNgl53OidgcV3a+rTl1eGv\n2u3XWjP+i/HkFORwZ4r8DN2ZjLR1IyNSIjh9qZzNR047XiiqN/hH2BY9Fy2SUoqM2AzWFq6lpELa\nZNyZJHw3kp4cjodFNa4B12KxjcDdmwWVFcYFJwyVHp1OSWUJa4+tNTsUYSJZ09bNTH57FVuOnCEy\nyIdXJ/eib4eQhgvtnAuf3QuBUZByE9z4ivGBimZVVlnGsE+HEe4bTkxgjN3xiR0ncmPijSZEJpqD\nK4y0FS7oidGdmdCrPcXnS/nXhnzHCiWNgYG/hMBI2PABlJ5vsIhwLV4eXkztOZVgn2AulF+ote06\ntYu3trxldojCCWR6ZDczICGEAQkhXCytIHvXcbTWKNVAA66nD4x7AQ6ugA/Gw77F0HWicwIWzeaB\nHg/wQI8H7PZ/nPsxf177Zw6dO0SHNh1MiEw4izzhu6nMlAiKzpWy4+g5xwvFDgSfIBmB28rILJzu\nQxK+mxreOQKAxY1pwPWwQqcbbHPsVFUZFJlwtpjAGBKDEiXhuwGp0nFT4YHe9IoJYsay/WTlFvHn\nW3vStX2bhgsmjYHtn8Nbg6HLRMh81vhgheEyYjKYlTuLJ5c+aXcsyDuIp/s/jaeHpwmRieYkT/hu\n7LEbkhmUGErusfPM2XDEsUIp46HnXbYFUla+CeWXjQ1SOMWEjhNIDEpk16ldtbZtJ7YxO282qwpX\nmR2iaAbSLVNw3/tr2X/iIkueGN5wA+4Ve7Lgo9vgnn9B0ihjAxSmKa0sJe3TNCZ2nMh/Dfovs8MR\n9ZBumcJhmV3acejkJfafuOh4ofhh4Okna+C2ct4e3gyKGkROfg6u9HAoro8kfEFmiq0BNzu3EQ24\nnj6QOBx2fyfr37Zy6THpFF4sZO+ZvWaHIppIGm0F0cG+dG4XyNs5+1m29wTPT+xGQph/wwWTRttW\nxnp/PHSZAIN/ZXywwunSotMA+PPaP5MYlGh3fGLHifQM7+nssMR1kIQvAPjViI58tPowK/ae4LP1\nR3h6bErDhbpOsiX847mw5M+25RGlJ0er086/HaM7jGZ90Xr2ndlX69i5snMcPneYGaNnmBSdaAxp\ntBW13DVjFWculbPgsXTHC+V+BbN/DFO+hoQ044ITLucv6//Ch7kfsvyu5fh7OvBXoTCENNqK6zIy\npR27jp2n4EwjulsmDgeLp0yh7IbSY9KpqKpg9dHVZociHNDURczvUErtUEpVKaX6XXXsGaXUXqVU\nnlJqTNPCFM4y4koDbmNG4HoHQvxQmXLBDaVGpBLoGcjSfFk8vSVoah3+duBW4O2aO5VSXYG7gG5A\neyBLKZWsta5s4vWEwTqG+9Mh1I/pS/ax7sAppo1LoX2wb8MFk8fCgmnwyY8g5UbofY/xwQrTeVo8\nGRI9hMVHFvPGxjfsjrfza8fkzpMdH98hDNWkhK+1zgXq+mFOAj7VWpcCB5RSe4EBgAzXc3FKKX6e\nlsisVYf4autROoYH8OgNSQ0X7HqzbcqFI6vhyBrodbdt8RTR6k1InMCSI0t4f3vthe6rqKJKV9G3\nXV86te1kUnSiJqN66UQDNSv18qv3iRbgx4M68ONBHbj57yvIzjvuWMJvEwUPZMHWOfDFA3B0I8Q0\n2IYkWoGM2AzW/9i+s8Wxi8cY9a9R5BTkSMJ3EQ0+gimlspRS2+vYJl2rWB376uwOpJSaqpRar5Ra\nX1xc7GjcwglGpkSw5cgZis+XOl6o00hQFqnPF0T6R5ISksLSI1K/7yoaTPha6xu01t3r2OZeo1g+\nEFvj+xjgaD3vP0Nr3U9r3S88PLxx0QtDXWnAXZLXiAZcvxCIGSBTLgjANmhrS/EWzpaeNTsUgXFV\nOvOAj5VSr2JrtE0CZPXkFqZb+za0a+PN2zn72VZwlkdHJhEa4N1wweQxsOh5mPsQJI+DLjcZH6xw\nSRmxGbyz7R1e2/AaCUEJtY4pFGPix9DOv51J0bmfJiV8pdQtwJtAOPCNUmqz1nqM1nqHUuozYCdQ\nATwkPXRaHqUUU4bEM3PlQf656hDRwb48mNGx4YLdb4VNH8L2L23LIabcCNJLwy11D+1OXGAcn+/5\nvM7j+87u4/khzzs5KvclI22FQ8a+nkOQryezHxzseKGN/4R5j8AvVkBkd+OCEy6toqqC0kr7dqDf\nr/g9m49vJuuOLCxKenQ1hYy0Fc0qMyWC9YdOc/ZyueOFkkbbvsoIXLdmtVjx9/S32zJiMii+XEzu\nqVyzQ3QbkvCFQ0Z2iaCySrNsTyN6UgVGQlQv2xTKQlxlWPQwFErW0nUimS1TOCQ1ti1t/Tx5e+l+\ndh87z8/TEwn0cWBmzKQxsOwVWPisrTE3oRGTsolWLdQ3lB5hPVhwYAGRfpF2x2MCY+gf2d+EyFov\nSfjCIR4WxeT+scxceZBtBWcJDfBmypD4hgt2vw3WvQur37J11Xxkg+GxipZjdPxoXln/Cs+tfM7u\nmFVZWXLnEoK8g0yIrHWSRlvRaCNeWUJciB8zfzbA8UJrZsD8J+GRjRDqQE8f4Ra01hRdKrJbPnH3\n6d08nP0wL6e/zNiEsSZF13JIo60wzIjOEazaf5JLZRWOF0q+0oAr9fniB0opIv0jiQqIqrUNix5G\nW++2MgtnM5OELxptZJcIyiqqWLn3pOOF2sZDWGcZgSsc4mHxYGj0UJYXLKeySobwNBepwxeN1j8+\nhABvK+8s28/hU5f40cA4fDw9Gi6YPBpWT4dlf4GOmdC+t/HBihYrIyaDr/d/zce7PqZ9QPtax6zK\nysCogfhYfUyKrmWShC8azctq4aaeUXy67ghrDpwiwNvK5P6xDRfsdqst4S/6f7ZRuL9cbnywosUa\nEj0EHw8fXlr3Up3Hf9371/y858+dHFXLJo224rporSmtqGLEK0voFRPM9J/0daxgZTmsfMOW9H+z\nE4Jk1mxRv6KLRZwuPW23/7kVz+Fp8eSjGz8yISrXI422wlBKKXw8PRiREsGyPcWUVjhYz+rhCZ1v\ntL2WBlzRgHb+7UgJSbHbRsSNYNuJbZy83Ih2JCEJXzRNZucILpZVsu6A/VNYvcI7Q3CczJkvrltG\nTAYazfICqRZsDKnDF00ypFMoXlYL7684wPHzJUzs1R6rRwPPEUrZRuBu+hDW/wPi0yFMVkQSjusS\n0oVw33C+PfCtXYMuQExADFEBUSZE5tok4Ysm8fOyMjIlgvnbj7Fo13GsHhYm9rL/D2in282w7h34\n+je2BVMe+N74YEWroZRiROwIPtv9GSuPrrQ7HuYbRtbtWXhYHOg95kYk4Ysme/Pu3py4UMb4N5aR\nnVvkWMKPHwZPHYAVr8OKN+DiCfAPMz5Y0Wr8tt9vGZsw1m6U7qbjm/jb5r+x7cQ2UiNSTYrONUkd\nvmgyq4eFyCAfhieHs2R3MZVVDvb88guBrjcDGvZmGRqjaH38PP3oH9mfAVEDam13pdyFh/KQWTjr\n0KSEr5S6Qym1QylVpZTqV2P/KKXUBqXUtuqvmU0PVbi6zC4RnLlUzqbDjWjAjUoF/wgZgSuaTZB3\nEKkRqSwrWGZ2KC6nqVU624Fbgbev2n8CmKC1PqqU6g4sBKTDdSuXlhSOh0Uxc9UhzpWUMzw5Aoul\ngaUNLRbbCNydX8GOL231+dI3XzRRekw6r214jXXH1tnNtulp8SS+TTzKDZfdbJaBV0qpJcATWmu7\nUVPKdldPAO211vbrnNUgA69avnv/sZac3bZFUt65tx+jujqwQPXuhfDxZNvrxOFw71zD4hPuYf+Z\n/UyaO6ne461tFk5HB145o9H2NmBTQ8letA7Tf9yHw6cuccdbq1iUW+RYwk8aDb/eZGu83TQLSs6B\nTxvjgxWtVmJwIu+Peb/OUbp/WvMnsg5ntaqE76gGE75SKguwX44GntVaX/NRTCnVDXgRGH2Nc6YC\nUwHi4uIaCke4OD8vKymRbUhLDiN713G01g3/6awUhCRCz8mw4X3Yl23rtilEE/SLrPuBd3nBcr47\n+B3lVeV4WhxYta0VabDRVmt9g9a6ex1bQ8k+BvgSuFdrve8a7z9Da91Pa90vPDy88f8C4ZIyU9px\n/HwpO46ec7xQzADwCZYpF4Sh0qPTuVB+gc3HN5sditMZ0i1TKRUMfAM8o7VeYcQ1hGsb3jkcpeDj\ntYfZePi0XV/pOnlYodMNtjr9/Uvg0inD4xTuZ1D7QVgtVhYeXMixi8fstrLKMrNDNEyTGm2VUrcA\nbwLhwBlgs9Z6jFLqv4BngD01Th+ttT5+rfeTRtvW5fa3VrL+kK0O9aMHBjK0kwMDq7Z/Af+6z/a6\n841w98cGRijc1YPfP1jnCF2A1PBUZo2f5eSImsbRRluZHlkY5tTFMvKOnee+D9ZyV/84/jixW8OF\nqqrg6CZY8xbkfg1PHwBPX+ODFW4l/3w+awrX2O1fV7SOb/Z/Q9btWbTzd6DDgYtwpV46wk2F+Hsx\nuGMoQzraGnD/MKFrww24FgvE9IWSu2HbHDiw7If1cIVoJjGBMcQExtjt7xXei2/2f8OygmXcnny7\nCZEZS6ZWEIbLTIng8KlL7Cu+6Hih+GHg6S8jcIVTdQzuSHv/9q12WgZJ+MJwI1IiAMjeVeR4Ias3\ndBxha8A9uhnKLhkUnRA/UEqRFpPG6sLVnC09y+WKy7W28spys0NsEqnSEYaLDvYlJTKQ7F3HmZre\n0fGCncfBrq9hRgb0uANue9e4IIWolhGTwey82Qz7dJjdMauyMmv8LLqHdTchsqaThC+cYkRKBO/k\n7Ofs5XKCfB0c7NLzLgiMhLXv2p70K8ttSyQKYaAh7Yfw+0G/50L5hVr7tdb8bfPf+O7gd5LwhbiW\nkSkRvLVkH8v2FHNTTwfmy4cf+uWXl8Du+XB4FSSkGxuocHseFg8md55c57HVhatZmr+Ux/s97uSo\nmofU4Qun6B3XlmA/T7J3XXMoRt0Sh4OHl6yBK0yXHpPO/rP7OXL+iNmhXBdJ+MIpPCzKtkBKXjFH\nTl2ivLLK8cLeAbZeO7sXwOlDUFlhXKBCXENGTAZAi+3FIwOvhNN8teUoj3yyCYCx3SKZ/pO+jhde\nMwPmP2l73fenMOGvzR+gEA6Y8OUEDp8/jFXZ14jf3Olmfj/4906PSUbaCpdTUVnF/O3HmLv5KMv2\nFLP5udH4ejm4yHR5CeR+BRtnwvGd8MQekAWqhQnWHVvH8oLldvs3Fm0k73Qey+9ajpeHl1NjkpG2\nwuVYPSxM6NWetn5eZOUWsXLfCUZ2cXD4uqcP9LzDNpXy5/dDwUaI7W9swELUoX9kf/pH2n/2cvJz\neGjRQ6w/tp4h0UNMiKxhUocvnG5AQgj+Xh4sup4G3I6ZoCywRxpwhWvpH9kfbw9vluYvNTuUeknC\nF07nZbUwLCmM7Nzjjk2bXJNfCMQOhLwFUHIWXKhKUrg3X6svA6MGkpOf0/jPtZNIlY4wxciUdizc\nUURu4Xm6tm/kcobJYyDrj/BCHAz5NYz+b0NiFKKx0qPTycnPYcRnI+wmCgzwDOC9Me8R4RdhUnSS\n8IVJhqfYVjdbnHe88Qm/3/3gFQBbPoHtn8Oo/2er2xfCZOMTx3Pw3EFKKktq7S+vLGfuvrlkHcri\nR11+ZFJ0kvCFSSICfegRHcSi3CIeGtGpcYV92sCAn4PVB+Y9DEXbIbKHMYEK0QiBXoE8PeDpOo9t\nLt5MTkGOqQm/SXX4Sqk7lFI7lFJVSim7LkFKqTil1AWl1BNNuY5onTJTIth05AynLl7nknJJ1fPk\nywhc0QKkRaexrnAdl8rNm/m1qY2224FbgfqGnb0GzG/iNUQrlZkSgdawdPd19NYBCGwHUamy6Llo\nETJiMyirKqtzpS1naZaBV0qpJcATWuv1NfbdDAwFLgIXtNavNPQ+MvDKvVRVaQb8aREnLpRiUTDv\n4WF0jw5q3Jss/jMsfcHWVbPPvTICV7is8spy0man4e3hTVvvtnbHh0UP44n+11cZYurAK6WUP/A0\nMAqQ6hxRJ4tF8dqdvVh74BQA4YHejX+TAT+3Jfuqcmjfu5kjFKL5eHp48tt+v2XV0VV1HndG750G\nn/CVUllAZB2HntVaz60+Zwk1nvCVUq8Aa7XWnyml/sg1nvCVUlOBqQBxcXF9Dx06dJ3/FCGEcE/N\n9oSvtb7hOq4/ELhdKfUSEAxUKaVKtNZ/q+P9ZwAzwFalcx3XEkII4QBDqnS01mlXXtd4wrdL9kII\nIZynqd0yb1FK5QODgW+UUtI/TgghXFSTnvC11l8CXzZwzh+bcg0hhBDNQyZPE0IINyEJXwgh3IQk\nfCGEcBOS8IUQwk241Jq2SqlioCkjr8KAE80UTnORmBzjijGBa8YlMTnOFeMyIqYOWuvwhk5yqYTf\nVEqp9Y6MNnMmickxrhgTuGZcEpPjXDEuM2OSKh0hhHATkvCFEMJNtLaEP8PsAOogMTnGFWMC14xL\nYnKcK8ZlWkytqg5fCCFE/VrbE74QQoh6tIqEr5Qaq5TKU0rtVUpNMymGWKXUYqVUbvU6v49W7/+j\nUqpAKbW5ehtvQmwHlVLbqq9/Zc2CEKXU90qpPdVf7ZfgMS6ezjXux2al1Dml1GPOvldKqX8opY4r\npbbX2FfnfVE2b1R/xrYqpfo4Oa6XlVK7qq/9pVIquHp/vFLqco17Nt2JMdX781JKPVN9r/KUUmOc\nGNPsGvEcVEptrt7vrPtUXx4w/XMFgNa6RW+AB7APSAS8gC1AVxPiiAL6VL8OBHYDXYE/Ylscxsx7\ndBAIu2rfS8C06tfTgBdN/PkdAzo4+14B6UAfYHtD9wUYj219ZgUMAtY4Oa7RgLX69Ys14oqveZ6T\nY6rz51X9ud8CeAMJ1f8/PZwR01XH/wI85+T7VF8eMP1zpbVuFU/4A4C9Wuv9Wusy4FNgkrOD0FoX\naq03Vr8+D+QC0c6OoxEmATOrX88EbjYpjpHAPq2105c601rnAKeu2l3ffZkE/FPbrAaClVJRzopL\na/2d1rqi+tvVQIwR125MTNcwCfhUa12qtT4A7MX2/9RpMSmlFDAZ+KS5r9tATPXlAdM/V9A6qnSi\ngSM1vs/H5ESrlIoHegNXlqd/uPrPtX84s+qkBg18p5TaoGxLSgK001oXgu1DChi/oGbd7qL2f0qz\n71V998WVPmc/w/ZUeEWCUmqTUmqpUiqtvkIGqevn5Qr3Kg0o0lrvqbHPqffpqjzgEp+r1pDwVR37\nTOt6pJQKAD4HHtNanwPeAjoCqUAhtj8znW2o1roPMA54SCmVbkIMdpRSXsBEYE71Lle4V/Vxic+Z\nUupZoAL4qHpXIRCnte4NPA58rJRq46Rw6vt5ucK9upvaDxJOvU915IF6T61jn2H3qjUk/Hwgtsb3\nMcBRMwJRSnli+yF/pLX+AkBrXaS1rtRaVwHvYMCftg3RWh+t/noc24I1A4CiK386Vn897uy4sP0C\n2qi1LqqOz/R7Rf33xfTPmVJqCnATcI+urgCurjY5Wf16A7b68mRnxHONn5ep90opZQVuBWbXiNVp\n96muPICLfK5aQ8JfByQppRKqnxjvAuY5O4jqOsP3gFyt9as19tesj7sF2H51WYPj8ldKBV55ja3x\nbzu2ezSl+rQpwFxnxlWt1lOY2feqWn33ZR5wb3WvikHA2St/ojuDUmos8DQwUWt9qcb+cKWUR/Xr\nRCAJ2O+kmOr7ec0D7lJKeSulEqpjWuuMmKrdAOzSWudf2eGs+1RfHsBVPldGt1o7Y8PW0r0b22/t\nZ02KYRi2P8W2Apurt/HALGBb9f55QJST40rE1mNiC7Djyv0BQoFFwJ7qryFOjssPOAkE1djn1HuF\n7ZdNIVCO7Unr/vruC7Y/vf9e/RnbBvRzclx7sdX1XvlsTa8+97bqn+sWYCMwwYkx1fvzAp6tvld5\nwDhnxVS9/wPgF1ed66z7VF8eMP1zpbWWkbZCCOEuWkOVjhBCCAdIwhdCCDchCV8IIdyEJHwhhHAT\nkvCFEMJNSMIXQgg3IQlfCCHchCR8IYRwE/8fWODaiFQzKp8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fd8bf373518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for n in [2, 3, 10]:\n",
    "    print('n = %d' % n)\n",
    "    gaps = []\n",
    "    P, q, A, b = generate_problem(n)\n",
    "    P *= 2\n",
    "    x = np.array([0.0] * n)\n",
    "    t = 1.5\n",
    "    mu = 1.5\n",
    "    eps = 0.000001\n",
    "    while True:\n",
    "        x = barrier_step(x, t, eps)\n",
    "        if n / t < eps:\n",
    "            break\n",
    "        t *= mu\n",
    "    print(x, x.dot(P).dot(x)/2 + q.dot(x))\n",
    "    plt.plot(np.log(gaps), label=str(n))\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Задача 3.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим ограничение вида\n",
    "$$a^Tx \\le b, \\quad \\forall [a, b]\\in U.$$\n",
    "Здесь\n",
    "$$U = \\{[a, b]\\ \\mid\\ [a, b] = [a_0, b_0] + \\sum_{i=1}^m [a_i, b_i]\\zeta_i,\\ \\zeta\\in Z\\},$$\n",
    "где\n",
    "$$Z = \\{\\zeta\\in\\mathbb R^m\\ \\mid\\ 0\\le p\\le\\zeta_1\\le\\cdots\\le\\zeta_m;\\ \\lVert\\zeta\\rVert_2\\le\\Omega\\}.$$\n",
    "Используя метод, показанный на лекции, постройте эквивалентное конечное множество выпуклых ограничений."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Заметим, что $Z$ записывается в виде\n",
    "$$Z = \\{\\zeta\\in\\mathbb R^m\\ \\mid\\ G\\zeta + g\\in\\mathbb R^m_+,\\ Q\\zeta + q\\in\\mathbb L^{m+1}\\},$$\n",
    "где\n",
    "$$G_{ij} = \\begin{cases}\n",
    "1,&\\mbox{если $i=j$}\\\\\n",
    "-1,&\\mbox{если $i = j + 1$}\\\\\n",
    "0,&\\mbox{иначе}\n",
    "\\end{cases}\\in\\mathbb R^{m\\times m}, \\quad\n",
    "g = (-p, 0, \\ldots, 0) \\in\\mathbb R^m\n",
    "$$\n",
    "$$Q_{ij} = \\begin{cases}\n",
    "1,&\\mbox{если $i=j$}\\\\\n",
    "0,&\\mbox{иначе}\n",
    "\\end{cases}\\in\\mathbb R^{(m + 1)\\times m}, \\quad\n",
    "q = (0, \\ldots, 0, \\Omega) \\in\\mathbb R^{m + 1}\n",
    "$$\n",
    "$$\\mathbb L^{m+1} = \\left\\{(z, t)\\in\\mathbb R^{m+1}\\ \\mid\\ t \\ge \\sqrt{z_1^2 + \\cdots + z_m^2}\\right\\}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Учитывая, что двойственный конус к декартовому произведению конусов есть декартово произведение двойственных конусов и используя теорему из лекции об общем виде конечных множеств выпуклых ограничений, мы получаем, что $Z$ состоит из векторов $\\zeta\\in\\mathbb R^n$, для которых существуют вектора $y\\in\\mathbb R^m, z\\in\\mathbb R^m, t\\in\\mathbb R$ (это части вектора $y\\in\\mathbb R^{2m+1}$ из теоремы), что\n",
    "$$g^Ty + q^T(z, t) + a_0^Tx - b_0 \\le 0$$\n",
    "$$(G^Ty)_i + (Q^T(z, t))_i + a_i^Tx - b_i = 0, \\quad i=1\\ldots m$$\n",
    "$$y\\in (\\mathbb R_+^m)_*$$\n",
    "$$z\\in (\\mathbb L^m)_*$$\n",
    "где $(z, t) = (z_1, \\ldots, z_m, t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Учитывая, что оба рассматриваемых конуса самодвойственны, и подставляя формулы для матриц и векторов, мы получаем"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$-py_1 + \\Omega t + a_0^Tx - b_0 \\le 0$$\n",
    "$$y_i - y_{i+1} + z_i + a_i^Tx - b_i = 0,\\quad i=1\\ldots m-1$$\n",
    "$$y_m + z_m + a_m^Tx - b_m = 0$$\n",
    "$$y_i \\ge 0,\\quad i=1\\ldots m$$\n",
    "$$t\\ge\\sqrt{z_1^2 + \\cdots + z_m^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Учитывая, что $\\Omega > 0$, можно убрать последнее неравенство, а в первом заменить $t$ на $\\sqrt{z_1^2 + \\cdots + z_m^2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Таким образом, мы получили конечное число выпуклых неравенств, что и требовалось."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}