MVA-2021 / convex_optim / exam_svm.ipynb
exam_svm.ipynb
Raw
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 6 - Pierre Fernandez"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given $m$ data points $x_i \\in \\mathbb{R}^n$ with labels $y_i \\in \\{-1,1\\}$. The goal of\n",
    "this exercise is to write a function to solve the classification problem:\n",
    "\n",
    "$$ \n",
    "\\begin{aligned} \n",
    "&\\textrm{min}_{w\\in \\mathbb{R}^n, z\\in \\mathbb{R}^m} \\quad && \\frac{1}{2} \\| w\\|^2_2 + C\\mathbf{1}^Tz \\\\ \n",
    "&\\textrm{subject to} && y_i(w^T x_i) \\geq 1 - z_i, \\quad i=1,...,m \\\\\n",
    "&                    && z \\geq 0\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1\n",
    "\n",
    "#### Derive the dual problem.\n",
    "\n",
    "The associated Lagragian is ($\\alpha, \\pi \\in \\mathbb{R}^m$): \n",
    "$$ L(w,z,\\alpha, \\pi) = \\frac{1}{2} \\| w \\|^2_2 + C\\mathbb{1}^Tz + \\sum_{i=1}^{m} \\alpha_i (1 - y_i(w^Tx_i)-z_i) -\\pi^Tz $$\n",
    "\n",
    "Therefore $$ g(\\alpha ,\\pi) = \\inf_{w,z} L(w,z, \\alpha, \\pi) = \\sum_{i=1}^{m} \\alpha_i + \\inf_{w} \\left( \\frac{1}{2} \\| w \\|^2_2 - \\left(\\sum_{i=1}^{m} \\alpha_i y_i x_i\\right)^T w  \\right) + \\inf_{z} \\left( C\\mathbb{1}^Tz - \\alpha^Tz - \\pi^Tz \\right) $$\n",
    "\n",
    "Now, recall (from example 3.27 from the book and the fact that $L_2$-norm is self-dual) that the conjugate of $f=\\frac{1}{2}\\|.\\|_2^2$ is $f^\\star=\\frac{1}{2}\\|.\\|_2^2$.\n",
    "\n",
    "Hence:\n",
    "- if $C\\mathbb{1} -\\alpha -\\pi \\neq 0$, $g(\\alpha ,\\pi) = -\\infty$\n",
    "- otherwise, $g(\\alpha ,\\pi) = \\mathbb{1}^T \\alpha - \\frac{1}{2} \\| \\sum_{i=1}^{m} \\alpha_i y_i x_i \\|_2^2$\n",
    "\n",
    "Therefore the dual problem is:\n",
    "$$ \n",
    "\\begin{aligned} \n",
    "&\\textrm{max}_{\\alpha ,\\pi} \\quad &&\\mathbb{1}^T \\alpha - \\frac{1}{2} \\left\\| \\sum_{i=1}^{m} \\alpha_i y_i x_i \\right\\|_2^2 \\\\\n",
    "&\\textrm{subject to} &&C\\mathbb{1} -\\alpha -\\pi = 0 \\\\\n",
    "& && \\alpha \\geq 0 \\\\\n",
    "& && \\pi \\geq 0 \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Which is equivalent to:\n",
    "\n",
    "$$ \n",
    "\\begin{aligned} \n",
    "&\\textrm{max}_{\\alpha} \\quad &&\\mathbb{1}^T \\alpha - \\frac{1}{2} \\left\\| \\sum_{i=1}^{m} \\alpha_i y_i x_i \\right\\|_2^2 \\\\\n",
    "&\\textrm{subject to} && 0 \\leq \\alpha \\leq C \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "\n",
    "Now let us rewrite this as a general Quadratic Problem.\n",
    "\n",
    "- Let $ A = \\begin{pmatrix}I_m \\\\ -I_m \\end{pmatrix} $ and $ b = \\begin{pmatrix}C.\\mathbf{1}_m \\\\ 0.\\mathbf{1}_m \\end{pmatrix}$. \n",
    "We have that $0 \\leq \\alpha \\leq C \\Leftrightarrow A\\alpha \\leq b$.\n",
    "\n",
    "- Let $X = (x^T_1,..., x^T_m) \\in \\mathbb{R}^{nxm}$ and $\\circ y = diag(y) \\in \\mathbb{R}^{mxm}$. Let $Q = \\frac{1}{2} (X\\circ y)^T (X\\circ y)$ and $p=-\\mathbb{1}_m$. We have that \n",
    "$\\frac{1}{2} \\left\\| \\sum_{i=1}^{m} \\alpha_i y_i x_i \\right\\|_2^2 = \\alpha^T Q \\alpha$, \n",
    "therefore $\\frac{1}{2} \\| \\nu \\|_2^2 + \\nu^Ty = \\nu^TQ\\nu + p^T\\nu$.\n",
    "\n",
    "With these notations, the dual problem becomes a QP: \n",
    "\n",
    "$$ \n",
    "\\begin{aligned} \n",
    "&\\textrm{min}_{\\alpha} \\quad &&\\alpha^TQ\\alpha + p^T\\alpha \\\\ &\\textrm{subject to} && A\\alpha \\leq b\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2 & 3\n",
    "\n",
    "#### Use the barrier method to solve both primal and dual problems up to arbitrary precision $\\epsilon$. Test your code on random clouds of points (e.g. generate two classes of data points by picking two bivariate Gaussian samples with different moments). Try various values of $C > 0$, e.g. $0.1, 1, 10, ...$\n",
    "\n",
    "#### Plot convergence of the objective function in Newton iterations for the barrier problem (i.e. inner iterations), with respect to the best value found, in semilog scale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\Phi_t(v) = t(v^TQv+p^Tv) + \\phi(v)$, where $\\phi(v) = -\\sum_{i=1}^{2m} \\log(b_i-A_i v)$ is the logarithmic barrier function (and $A_1, ..., A_{2m}$ are the rows of $A$)\n",
    "\n",
    "In order to perform the backtracking line search we first have to derive the gradient and Hessian of $\\Phi_t$.\n",
    "\n",
    "We have, with the notation $d(v)_i = \\frac{1}{b_i-A_iv}$, that $\\nabla\\phi(v) =  A^Td(v)$ and $\\nabla^2\\phi(v) = A^Tdiag(d(v))^2A$. Therefore $$ \\nabla\\Phi_t(v) = t(2Qv+p)+A^Td(v) $$ $$ \\nabla^2\\Phi_t(v) = 2tQ + A^Tdiag(d(v))^2A $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also want to retrieve the primal solution from the dual solution.\n",
    "\n",
    "First, note that the primal problem is convex since all the functions are either a sum of a square norm that is convex with regards to $w$, and an affine function with regards to $w$ or $z$. Then the point $(z=2, w=0)$ is strictly feasible so Slater's condition holds. \n",
    "- The KKT stationarity condition implies that $\\nabla_w L = 0 \\Rightarrow w^\\star = \\sum_{i=1}^m\\alpha^\\star_iy_ix_i$\n",
    "- From the complementary slackness, one has for $i=1,...,m$ that:  \n",
    "$\\alpha^\\star_i f_i(w^\\star, z^\\star)= 0 \\quad \\Rightarrow \\quad \\alpha^\\star_i (1 - y_i( (w^\\star)^Tx_i)-z^\\star_i ) = 0 \\quad$ and $\\quad (C-\\alpha^\\star_i)z^\\star_i= 0$  \n",
    "The last condition gives that $\\alpha_i^\\star = C$ if $z_i^\\star > 0$, and the first one gives $z^\\star_i = 1 - y_i( (w^\\star)^Tx_i)$ for $\\alpha^\\star_i \\neq 0$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [],
   "source": [
    "# backtracking line search parameters, chosen as they are in Boyd's book\n",
    "alpha = 0.01\n",
    "beta = 0.5\n",
    "assert 0 < alpha and alpha < 0.5 and 0 < beta and beta < 1.0\n",
    "\n",
    "def f(Q,p,v):\n",
    "    return v.T@Q@v + p.T@v\n",
    "\n",
    "def phi(A,b,v):\n",
    "    assert (b-A.dot(v)>0).all()\n",
    "    return -np.sum(np.log(b-A.dot(v)))\n",
    "\n",
    "# Phi_t at v\n",
    "def phi_t(Q, p, A, b, t, v):\n",
    "    return t * f(Q,p,v) + phi(A,b,v)\n",
    "\n",
    "def centering_step(Q,p,A,b,t,v0,eps):    \n",
    "    v_seq = [v0]\n",
    "    while True:\n",
    "        v = v_seq[-1]\n",
    "        # compute Newton step and decrement\n",
    "        d = 1/(b-A.dot(v))\n",
    "        D2 = np.diag(d.flatten()**2)\n",
    "        nabla = t*(2*Q@v+p) + A.T.dot(d) # gradient of Phi_t at v\n",
    "        nabla2_inv = np.linalg.inv(2*t*Q + A.T @ D2 @ A) # hessian of Phi_t at v\n",
    "        delta_v_nt = - nabla2_inv @ nabla # newton step\n",
    "        lambda2 = nabla.T @ nabla2_inv @ nabla # newton decrement\n",
    "        if lambda2/2 < eps: # stopping criterion\n",
    "            break \n",
    "        # backtracking line search to find t\n",
    "        gamma = 1\n",
    "        while not((b-A.dot(v+gamma*delta_v_nt)>0).all()) or phi_t(Q,p,A,b,t, v + gamma*delta_v_nt) > phi_t(Q,p,A,b,t,v) + alpha*gamma*nabla.T@delta_v_nt:\n",
    "            gamma *= beta\n",
    "        # update v\n",
    "        v_seq.append(v + gamma*delta_v_nt) \n",
    "    return v_seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {},
   "outputs": [],
   "source": [
    "def barr_method(Q,p,A,b,v0,eps):\n",
    "    # barrier method parameter\n",
    "    assert mu>1\n",
    "    v_seq = [v0]\n",
    "    m_A = A.shape[0]\n",
    "    eps_cs = 1e-6\n",
    "    t = 1\n",
    "    while True:\n",
    "        # centering step\n",
    "        v_seq_cs = centering_step(Q,p,A,b,t,v_seq[-1],eps_cs)\n",
    "        # update\n",
    "        v_seq += v_seq_cs[1:]\n",
    "        if m_A/t < eps: # stopping criterion\n",
    "            break\n",
    "        # increase\n",
    "        t *= mu\n",
    "    return v_seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of the X :  (50, 2)\n",
      "Shape of the y :  (50, 1)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x26162e3eb00>"
      ]
     },
     "execution_count": 261,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd5iU1fXA8e+Z2SlbWJZlURCkqAQExAJYMSJIsANGbNGgJkGNGkvUnyUmqNHYYxJb7BgLYCcqKGKwF7ADFlAEUXrbOv38/pgBlt3ZpezMvDOz5/M8++zOvTNzz77i2Xfue99zRVUxxhiTn1xOB2CMMSZ9LMkbY0wesyRvjDF5zJK8McbkMUvyxhiTxwqcDqC+iooK7d69u9NhGGNMTvnoo49WqWqHZH1ZleS7d+/O7NmznQ7DGGNyiogsaqrPpmuMMSaPWZI3xpg8ZkneGGPymCV5Y4zJY5bkjdkKS79bzqIvlxCLxZwOxZhtklWra4zJNku++Ynxx93CsoUrEJdQ3LaIK5+4kP4/7+N0aMZsFTuTN6YJ4VCYiw/5M4u//JFgXYhATZDVP63lqqNuYPXStU6HZ8xWsSRvTBNmTf2UYG2IhuW4Y9EYr074n0NRGbNtLMkb04Q1y9YRjUQbtYcCYVYsXu1ARMZsO0vyxjSh74E/S9ruL/Gz15C+GY7GmO1jSd6YJvTYoxv7HT0AX5FvY5vX76Xzbh05aPS+DkZmzNaz1TXGNOPKJy5g2oOv89J90wkFIww7ZTCj/nAkBR77X8fkBsmmPV4HDhyo+VCgLBQM8+nrc4iEIuw1tB9FbQqdDskYk8dE5CNVHZisz05HUuyzmXP5y+ibN67IiIajXHT/2Qw75WCHIzPGtEY2J59CtVV1XH3sjdSsr6W2so7ayjqCdSFu/929/PTtMqfDM8a0QpbkU+i9KcmnmmKRKK899maGozHGGEvyKVVbVUc02ri2SSQcpWZdjQMRGWNaO0vyKTRgeH9IciHbX+xj/2OSXhMxxpi0siSfQjvt2pHjLjgKf7EPkXibv9jHwBF7sdeh/ZwNzhjTKtnqmhT7zd9+xcARezHt4dcJB8IMPeVg9j9mALIh6xtjTAZZkk+DPYf0ZU+77d0YkwVaPF0jIjuLyP9E5EsRmSsiFyTay0VkuojMT3xv1/JwjTHGbItUzMlHgD+q6u7A/sC5ItIHuByYoao9gRmJx8YYYzKoxUleVZeq6seJn6uAL4HOwEhgQuJpE4BRLR3LGGPMtknp6hoR6Q7sDXwA7KiqSyH+hwDYoYnXjBOR2SIye+XKlakMxxhjWr2UJXkRKQGeAS5U1cqtfZ2q3qeqA1V1YIcOHVIVjjE5IxQM89H0z5j96meEAiGnwzF5JiWra0TEQzzBP66qzyaal4tIJ1VdKiKdgBWpGMvkt9qqOh69ZjKvP/4WqjDkxAM5/doTKW5b7HRoafHxa59z7ZjbNha0U1WuevIi9jtyH4cjM/mixaWGJb4AfAKwRlUvrNd+C7BaVW8UkcuBclW9rLn3ypdSw2b7xGIxzh10OYvm/UA4GAGgwFtA556d+Pcnt+AucDscYWpVrq7ilG7nEKwNbtbuK/Ty6Ld3Ut7RFqSZrdNcqeFUTNccBJwGDBWRTxNfRwI3AsNFZD4wPPHYmCZ9NP1zfpy/dGOCB4iEIqxYtJIPXvrYwcjS461n3k/arqrMnPRuhqMx+arF0zWq+jbQ1O2cw1r6/qb1WPDxQoJ1jeek66oDLPhkIQeOHORAVOlTs76WaDjSqD0cjFCzvtaBiEw+sto1Jmt07LEDviJvo3Z/iZ+OPZIuzspp+wzvj9vTeArKV+RlwC/2dCAik48syZuscdCoQRQW+3G5Nn0wFJfgK/Ty8zEHOBhZeuy2Vw8OPWkw/uJNG4X7i30ccOxAdt+vp4ORmXxie7yarLJ04XJuPv1Ovnx/PgC9Bu3KZY+cR+fdOjkc2darWV+Dy+2isGTLe/uqKu+/+BGvPPI/NKb8YuwQDjh2IC6XnX+ZrdfchVdL8iYr1VbVoaoUlxY5HcpWW/TlEm4eeyffffY9CuwxuDeXPnIeO+xc4XRoJs+le3WNMSlX1KYwpxJ8zfoaLhz8J+Z/9C2RcJRoOMrnb37JRQdfTTQSdTo804pZkjcmBV5/4m0iwchmG4PFojGq1tbwwcv5t/zT5A5L8sakwI8LlhJocFMTQCQUZtlCu9nbOMeSvDEp0GtQTwpL/I3a3Z4Cdtu7hwMRGRNnSd6YFDho9L6Ud2pHQb117x6/h136d2OPg3d3MDLT2lmSNyYFvD4P/3zveo783WGUdSilvFMZx11wJDe9erXt72scZUsojTEmx9kSSmOMaaUsyRtjTB6zJG+MMXksJTtDmfQJBcO8MfldPn9jHh17dGDEGUOp2Knc6bCMMTnCknwWq1lfw/kHXMXKH1YRqAni8XmYeOPz/G3an+h3UG+nwzPG5ACbrsliE296nmULVxCoid9JGQ6GCdQEufHUf5JNq6KMMdkrL5L8mmVrmfPOV6xdsd7pUFLqjcnvEg6GG7WvW7GeZd/brfLGmC3L6emacCjMrWfew9vPvo/H5yEUCDPsVwdz4b3j8mLTZ6+/8S5JALGYNtlnjDH15fSZ/CN/msg7z31AKBCmZn0t4WCY/018m8evf8bp0FLiqLOGN9oOz+USdunflfad2jkUlTEml+RskldV/nvvq402fg7WhnjhzqkORZVax54zgn2P2AdfoRdfkY/CNn7ady7nT5Mudjo0Y0yOyNnpGlXdeEGyoZr1dRmOJj3cBW7+/NQfWfjFIr76cAEVncvjmz+7c38qyhiTGTmb5F0uF7vu1Z0Fnyxs1Nc7zzZB7rFHN3rs0c3pMIwxOShnp2sAzr/zN/iLfLjc8V/D5XbhL/Fz7j/OcDgyY4zJDjl7Jg/Q54Be3DX7Jibf/DzffraIngN24YRLR9KlZyenQzMN/PD1j7xw1zSWfbeCvQ/bg8PPHJpTe7hmg5r1Ncx+5TNUlYEj9qKkrNjpkEwOsFLDJu1mTfuEa46/lUgoQjQSw1fkpW1FKfd8dDOl7ds0ev63n33Pk397lu8+X8xue3Xn5CuPo0e/rg5Enj3eeuZ9bvr1vzYuDY5Eolzy0O859MSDHI7MZIPmSg1bkjdpFYvFOKnzONYu3/xGtQJvASPPO5yzbx27WfsXb33JFUf8lVAgjMYUl0vw+L3cMuMv7N7MtZa5737NzEnv4C5wM/SUwfxswK5p+X2csGbZWk7b9TxCDVaSeQu9TPjmn1R0bu9QZCZbWD1545ifvl1OXXWgUXskFOHd52c1ar/zDw8SrA2hsfjJRyymBGuD3HXBQ02OcfeFD3P5iOt44c5pPPuPl7j4kD/z6DWTU/dLOOzNp99P2q4x5Y3J72U4GpNrLMmbtCpq4ycaiSXvKy3c7LGq8t1ni5I+d8HH3yVtn//xd7z8wGsEaoKoKhpTgrUhJt30PD8uWNqy4LNEsDZELBJt1B6NRJtcRmzMBpbkTVqVd2xHr0G74i7Y/J+av9jH6D8cuVmbiDRK/Bs0dZHxvSmzCAca1/dRhQ9e/Hg7o84u+x21T9IyHR5fAfsdvY8DEZlcYknepN2fJl1Ml1474S/xU1RaiMfnYfivD+EXY4c0eu6o849sVMrBV+Rj9AVHJX1vj9+7cQltfS634PF7UhK/07r33Zmjz/kFviIfIiAS/yN5+JlD2W2vHk6HZ7KcXXg1223dyvU8fdt/+eCljynboZTjLz6G/Y4akPS5qsrXsxaw6sc19Bq0Gx26JL9YGI1E+ee5D/Daf96gwFtAJBTh8N8M5fd3nJH0Tt+l3y3nt3tcnPSi5GML76bdDm1b/otmiTlvf8mMx99CVTn05MH0/3kfRMTpsEwWsNU1JuUqV1cxbs8/UrmqinAoAsTPLk+9+nhOvGxUy99/TRXLv19Jxx470KZdSbPPnfrgDO48/0FcBW4EiEZjXPbIeRwy5oBmX7dk/lIeHT+JL976ioqd2nHyFcdx4MhBLY7dmEyzJG9SbsL4SUy++QVCDebDvYVeJi+9P+M3Oq1buZ4PXvoYl9vF/kcP2OIfhp++XcY5Ay4jUB0glljJ4yvy8dsbT2HUeUc2+1pjso0toTQpN2vqp40SPIDHW8C3n36f8XjKOrRlxOmHMvy0Q7aY4AH+c+1TmyV4gGBtkIeuepJQko1ajMlVluTNdumwc3uSTQdHwhHKO5ZlPqBtNOftrzZL8PUtW2i7bpn8kZIkLyIPicgKEZlTr61cRKaLyPzEd9vlIo/88qKj8Rb6Nmtze9z06NeVLj/byaGott4OXSuStkdCUcp2KM1wNMakT6rO5B8BDm/QdjkwQ1V7AjMSj02e6HdQb86/8zcUtSmkqLQQr99D731349opufGf+eQrjmu0VNPr93DgqEGUljeup2NMrkrZhVcR6Q68qKr9Eo+/Boao6lIR6QTMVNVezb2HXXjNPaFgmMXzltCmvIQdu3VwOpxtMvWhGdx3yX+IhOOF0w7+5X5cdN/Z+It8W36xMVkkI6trkiT5dapaVq9/rao2mrIRkXHAOICuXbsOWLQo+W3txqRDJBxh5Q+rKW1fQnHbxnfVLv7qR2ZP+5TCNn4GH7ffVl3UNSbTsjrJ12dn8q2PqvLVhwv49PU5lLYv4edjDsiKRKqq3HPRI7x0/2toLLaxrMD4Zy9lwPA9Uz6W3dRkWqK5JJ/OTUOWi0inetM1rXLJQjQaZeoDr/Py/dOJhKMM+9XBjDr/CHyFNiUQjUa5/qQ7mDXtE0KBMF6fh3sveZQbXrqSPQ7e3dHYPn7tc6Y+OGPjnbThYPyGr2uOv42nlz+A1+9t7uVbFIvFmHTzCzx923+pWlNNtz5dOOeOM9hn2B4tjt2Y+tK5hHIKsKFY+FjghTSOlbX+euLfufePE5j/8UIWfrGYR695iksOHU802riqYGszc+K7zJr2CYGaILFojEBtkEB1gGt+eYvjx+fVCW8krfAoAp+8PifJK7bNg1c+weN/fYbK1VWoKt/P/YE/j7yRee9/0+L3Nqa+VC2hfBJ4D+glIktE5DfAjcBwEZkPDE88blUWfLKQWdM+JVi7KVmE6kIsmrck5ysk1lXX8dazHzBz0jtUra3erveY9vDrSRNpKBjm61nftjTEFmnuj0wsmrx08tYK1AZ54V9TN/t3AfGSwo+Oz586+CY7pGS6RlVPbqJrWCreP1fNefsrNNY4IdRVB/hs5pycrZMya9onXDvmNlxuF6pKNBzjD3f/lhGnH7ptb9Tc5SCHy20MPXkwH7z4UaM/QtFIjL0O7dui91790xokSeVMgEXzlrTovY1pyO54TaN2HctwexpXTvT6PVQ0UYUx21Wvq+Ga428jUBOktrKOuqoAoUCIf577wDZv0vGL04fgL258bcLj9dBr0G6pCnm7HHDMQA4cOSgen8TLNXgLvVz2yLkUliSveb+12u9UvnHnq4a69925Re9tTEOW5NPogGMG4PF5Gt3+73K7Oey0Q5wJqoXeef7DpOUMYpEorz/59ja919BTBrPPYf3xF/sQl+Ar8uIv9vHnp/6YdJOMTBIRLv/PH7jxlas56f9GM/baE3n4q3/w8+Obr2y5NfxFPkb/4Qh8Ddbj+wq9/Hr8CS1+f2PqS+fqmlbP6/dy+8xrGH/craxcsgoRoaSsmKuevDBn6pzX1QSY9uAM3n1hNu06ltG+U7ukNV8ikSh1VY33cm2O2+1m/LOXMvfdrxNLKNsw5MQDKW2fHXecigh9D+xF3wObvYdvu5zx15Mpbd+Gybe8QOXqarr33Zlz/n56s5uVG7M9rNRwBqgqPy5YRjQcoevuXXJmTXRddR3n7nsFKxavJFgbQiS+21I0cYdoff5iH3+behX9Bju79NGY1sipdfImQUTo0rOT02Fssxf/PZ0Vi1YSTKwVV1VCdSHcHjfeQi/hQBhVxV/sY/Bx+9H3oN4OR2yMaciSvGnSO89/uDHB1+cr9HLm9Sez8IsfiIQjDDnxIAYM758zn1CMaU0syZsmta1IXnI3Fo3Rb/DujDz3iAxHZIzZVra6xjRp1PmNV4C4XEKHLu3ZpX83h6IyxmwLS/KmSXsP3YNfjz8Br99DUWkR/hI/nXbdketfvtKmZozJEba6xmxR9boavvxgPqXt2/CzAbtYgjcmy9jqGtMiJWXFDBqxl9NhZMT6VZW89cwHBGoCDDp8L7r1sTtQTW6zJG9Mwgcvf8x1J9wGCLFIlIevnsjR4w7j7NtPt08vJmfZnLwxxO/s/euJtxOsDRGsDRIORQjVhXj5gRl8NnOu0+EZs90syZucEgqGefPp93juny/z9ezUlSP+5LUvcCWpDBmsDTL90TdSNo4xmWbTNSZnLP7qR/54yJ8JBsJEQhHcBS72GtqP8c9c2uKCZrEkJaEhXvG4qT5jcoGdyZucce3xt7J+VSV1VXWEg2ECNUE+mfEFU+55pcXvvc9h/YlGGm8U4i/2MfSUg1v8/sY4xZK8yQnLvl/B0oUrGu0lEqwNMfWBGS1+/6I2hVz68Hl4C714fAWICL4iH0NOOJCBv0jtxt2pFA6FWb10LZFwxOlQTJay6RqTEyLhaJMrXCLh1OwHe8iYA+hzwM+YOfEd6qrr2PfIfei9b3aW/lVVHrvuaSbfOoVYNEZBgZuTrxzNiZeNspVAZjOW5E1O6LxbR9pWtGHF4s234/P6PRx2auqmUzp0ac+YS45N2fuly9O3/5fJN79AILFPbAh4/LpnKCot4thzRjgbnMkqNl1jcoKIcNXEiygs8eMt9AJQWOKnW9+dOe7Cox2OLvMm3vT8xgS/QaA2yBM3POtQRCZb2Zm8yRl99v8Z//nuLl577E1W/rCKPQ7uw/5HD3B8q8BMi8ViVK6qStq3bvm6DEdjsp0leZNT2laU8stWeOZen8vlonPPTvw4v/HG6V137+JARCab2XSNMTno7NvG4ktMW23gK/Ry9m1jHYrIZCtL8sbkoP2PHsB1/72cPgf+jNL2beh38O7cMPUq9jmsv9OhmSxj0zXG5Ki9h+7B3kP3cDoMk+XsTN4YY/KYJXljjMljluSNMSaPWZLPE6FAiBWLVxIKhp0OxRiTRezCa46LxWJMGD+ZZ25/EQBxCSdceiyn/ul4q2Fitko0GmXxlz/iL/bRqceOTofTqqjWobXPQ+gNcO2EFJ+CFOyW0jEsyee4p26dwjO3v0iw3i3uk256gZKyYkaff6SDkZlc8OHUT7hp7L8IB8JEozF27rUT45+9lI7dd3A6tLynsWp09S8hugyoA9xo3dNo29twFQ5P2Tg2XZPjJt/ywmYJHuK7GU288XmHIjK5Ysn8pVw75lYqV1VRVx0gVBdi4eeLuGToeNsoJQO05mGI/kQ8wQNEgQBUXoFq6qZdLcnnMFWlcnV10r71K9dnOBqTa17896uNyjTHYkrl6iq+ePPLjMSgGkbrniO25nfE1l2EBj/IyLhZIfgqEEzSEYXI/JQNY9M1OUxE6NJrJ5Z8/VOjvm59d3YgIpNLVv6wmmgTtfjXLEt/oTPVCLpmLITnsuFsVoOvo8Vn4Sr5fdrHd5yUJG/XKEhxyoaxM/kc9/u/n964hkmR1TAxWzZg+J74i32N2iPhKLvvn4HNUgKvQmRTggdA66D6bjS6Mv3jO0yKTgMpbNDqgoLuSEG3lI1jST7HDTp8b26YehV7DulLux3bsvewftz06p/tdnezRcN+NZiKLu3x+j0b2/zFPkacfmhGLrxq8LV4Um9IPBBqBdM2/iOg8CTAFz+rlyJwd0ba3ZPSYUQbbpqZYiJyOPAPwA08oKo3NvXcgQMH6uzZs9MajzFmk5rKWp694yXeeOpditoUMvLcIxh6yuCMLL+Nrb8W6p4AGlzklWKk7A7Ed0jaY8gGGl0O4U/BVQGefbbr2IvIR6o6MGlfOpO8iLiBb4DhwBJgFnCyqs5L9nxL8sa0Hhr+Cl19AhDYvEPaITu8jYgn6etMY80l+XRP1+wLLFDV71Q1BEwERqZ5TGNMDhBPbyi9mk3TFcXgao+UP2wJPoXSvbqmM/BDvcdLgP3qP0FExgHjALp27ZrmcIwx2cRVNAb1HwHh2fE5ac8A4hMAJlXSfSafbHJps/khVb1PVQeq6sAOHTqkORxjTLYRVwniG4J497UEnwbpTvJLgPoLtrsAjRd1G2OMSYt0J/lZQE8R6SEiXuAkYEqaxzTGGJOQ1jl5VY2IyHnAK8SXUD6kqnPTOaYxxphN0l7WQFVfBl5O9zgmuRmPv8WE8ZNYtWQNXXt35nc3n8qA4Xs6HZYxJkPsjtc89tL90/n7Wf9m6bfLCQfDfPvZ9/xl1M18POMLp0MzxmSIJfk8FYvFeOiqJxuXIa4L8eDljzkUlTEm06wKZZ6qqw5Qs742ad/iJFUrjdlWGl0FwdfjD3xDEXeFswGZpOxMPk8VlvjxFzWuMAjQsZvdj2BaJlb7FLryULTyerTqenTlocRqn3E6LJOEJfk85XK5OPmK0Y0Sva/Iy+nXneRQVCYfaPRHqLyW+IYXdYlKkkGoHI9GlzkcnWnIknweO+HSkfz62hNoU16CyyVUdC7n4vvP5qBR+zodmsllgWk0uHG9Xt8rGQ3FbJnNyecxEWHMxcdy/EXHEA5F8HgLMlJC1uQ5jRDfj7ShGJC6vUlNatiZfCsgInh9HkvweUo1Qqz6fmIrhhJbvj+x9Vei0RXpG9A/DEhWJdIFvmHpG9dsF0vyxuQ4XXcJVP8LYktA10Dd8+jq0WisKi3jScFuUHwm4CeeQlzxn0vGIQU90jKm2X42XWNMDtPI9xCcQfwi6AYRiFWhtU8hJWemZVxXmwtR/y/QwEuAIP4jEU+ftIyloc/Qmn9DZCF490aKz0rpHqj5zpK8MbksPBekADTYoCMA4Y+A9CR5APH0SVti30AD/0PXXUD8j5hC3fdoYCq0fyr+icJskU3XGJNjVBUNTCO25gy05j7QZBc7PVCwS8ZjSyVVRSvHE98ecMNqnihoLVp5s3OB5Rg7kzcmx2jl1VD3X6Au0ZLkgroUIEUnZzKs1NP1EFuVrCPxKcVsDTuTNyaHaGQB1L3ApgQP8bNcF/Fq3h5w90DaPYy4d3IkxpSRIpJvLge42mU0lFxmSd6YXBL6kOSJLwaFY5Ad3sLV4RXEu0+mI0s5ES8UjgIalucohKLfORFSTrIkb0wukTJIug+qF9w7Iq7yjIeUTlJ6dWJdvhekBPBB0a+RohOcDi1n2Jy8MbnEPxQqkyV5F1J4XMbDSTcRH1J2BxpdDbGl4O6OuEo29mt4HoTei/8B8B+OuNo6GG12siRvTA4R8UP5BHTtWaA1xKduBCm7DXF3cjq8tBF3e3C33/g4vvLmSqh7CYiAeKDqBii7F/Ed4FygWciSvDE5Rjx9ocObEJkTXz7p6Y9IsjIDeSz4OgReJr68kkQ9HdB158EO78Xn8w1gSd6YnCTiAk9/p8NwjNY9myhx3KgHQrPBd2DGY8pWduHVGJODYs30NVEGuZWyJG+MyTlSOAooStKj4B2Y6XCymiV5Y0zu8Q0H3xCgkPjFZy/gR8r+jkjybS9bK5uTNybNNFaN1tydKEUgUHgcUnIWIoVOh5azRFxQ9ncIf4IG30FcbcB/FOK2/YsbsiRvTBqpRtA1J8fL5BKKN9Y8iIbehfJJrXYjF9UwWn0v1D0BsVrw7o+UXoEUdE/0RyDyJeCDgp5Jj5OIgHefvLi7N50syRuTTsGZEP2BjQk+3giRbyD0PmRgTbdGl6LVd0PoHXBVIMW/Q/zD0z5uszGtvwwCM9i4BDI0E139EVS8DJF56LpLgQgQA1cHaHevlRbeTjYnb0waafgL0NokHSEIf5H+8aPL0VUjoe4ZiC6B8Kfo+kuIVd+f9rGbjCmyBAKvsTHBx1tBA2j1XejaP8QrUGpNfJlk9Ad0zWlo0pLKZkssyRuTRuLuTPziYMMOH7g7p318rbkftJr4WfGGxjqo/hcaq0n7+ElF5kPSm5VCiTtYQw3a438ACL6dgeDyjyV5Y9LJf2QiodWfU3aB+MF/WPrHD73PZgl+AymA6LfpHz+Zgu5NbHQCUEnyNfAKsTXpiymPWZI3Jo3EVYK0fxIKdgc88S/PHkj5pMws9XM1Uc9Gw+DaIf3jJyEFPcA7gPiyx62kUVv/vp3swqsxaSYFuyEVz6OxNYAgGdzwQkrGoWs+YPP5by94ByHujhmLoyEpuwutui6xrHTDWX1Td6oWQuFI27x7O9mZvDEZIq7yjCZ4APEOgtJrQdokdlrygu9ApOwfGY2jUVyuIlxt/4bs+Bm0fx7wN/HMYqTtDUjpNZkML6/Ymbwxec5VNAotPAqii8HVLqs2FhFxQ0Fv1F2RWGpaXyG0uRwpPCotY2usCq35d/xir3ih8ESk+LS8q+hpSd6YVkDEAwW7Oh1GUiICZXeha04DIvHlpRSA7yCk6Pi0jKkaQlePiS8r3bCap/oONPwh0u7etIzpFEvyxuQ5Dc5EayZAbC34hiHFYxFXaWrHiFWh1XdB4CXADYW/RErGbfXFZfH0hh3eiq+fj60E70DEs0dKY9xMYBrElrH5cs0ABN9Fw/MQT5/0jZ1hluSNyWOx6nug5t5Ntdcj36KB56D9lM220WsJ1TC6+sT4dNDG0g33oaH3ofyxrS7dIOKHwqNTEtOWaGh28pvUAMKfQx4l+RZdeBWRMSIyV0RiIjKwQd8VIrJARL4WkREtC9MYs600th6q726wuUYQoqvQ2kmpGyg4A2I/0bh0w1wIf5y6cVLJ3QVI8ilD3ODgqqN0aOnqmjnAccCb9RtFpA9wEtAXOBy4WyTpFvPGmHQJz4nvfdpIIL59Xopo6NMmSjeEM1K6YXtI4XHxhL4ZV3wVknewIzGlS4uSvKp+qapfJ+kaCUxU1aCqLgQWAPu2ZCxjzDZylQPRJB0C7h1TNowU7Ezy0g1ecO+UsnFSSdwVSLsJ4O5K/Izem7hJ7QlE8msWO12/TWfg/XqPlyTaGhGRccA4gK5du6YpHGNaoYLe8WmJyHdsnux9SNHY1I3jPwaqbm9wL5MLpBh8h6ZunBQT755QMR1iSwFP3tQrpH4AAAr4SURBVNai3+KZvIi8JiJzknyNbO5lSdqS3s6mqvep6kBVHdihQ34eZGOcICJIuwehoBfgBymJJ97Sv8QTXKrGcZUi5Y8nxvESL92wJ1L+ZNavORcRxL1TShK8xmrRuufR6nvR4HuoZsdes1s8k1fV7amitATYud7jLsBP2/E+xpgWEHfHeEmFyEKIVYKnd1pq5oinN1LxXzS6GsSNuMpSPkY208gCdPXJQDheMVP88U9S5RMc344wXWUNpgAniYhPRHoAPYEP0zSWMWYLpKAH4t0z7QlH3O1bXYIH0LUXgFYmLkDH4t/D89Cah5wOrcVLKEeLyBLgAOAlEXkFQFXnApOBecA04FxVTXYFyBhjcppGlyXuEWg4PROAumedCGkzLbrwqqrPAc810Xc9cH1L3t8YY3Kb8/Py+bVWyBiz3TT6U/wmqfAcwAUFvZDCYxBPL6dDy2ri7oi6uyTZhMUHhaMdiak+S/LGGDQ0C137W9AgG3dmCr2B1k5AS87CVXKeM3FFl0PgFSAMvkORgl0ciWNLpOwOdM2vQCNAbbysc0FPpPg3TodmSd6Y1k5V0XWXNSh/sEEQqu9D/UcgGa5iGaudApVXEV+RHYWqO9Di3+Jqc0FG49ga4ukFHWZuKnzm2RO8ByHi/JYdluSNae1iyyC2qpknRCAwA0oyl+Q1tiaR4IObd9Q8iPqHIZ5+GYtla4mrBNJUGrklnP8zY4xxlvho/gKhK0mdlzQLzgSSjRlC617MbCw5zpK8Ma2cuMrj0wtNpgMBX4YLySrJ75tHyYYVK7nEkrwxrYBqHbGax4mtOYPYuj+ioU8265ey2xPFurz1Wl2AD9pcgRR0yWS44B8CSW+t8SH+IzMbS46zOXlj8pxqXXyru8gPQB0gaGA62ub/cBX/CgBx7wgVr0B4Nhr+CmJrEXf7+E5SDtRXF1c5WjoeKscTP3OPAh4oOjWldXdaA0vyxuQ5rX0KIouBwIaW+M9VN6GFIzfuECUi4B2EeAc5FepmXEW/RH37x1esaAh8Q23N/nawJG9Mvgu8yqYEX48UQPgz8B2U8ZC2lrg7wxbWmmt0NVpzDwT/B1KKFJ8B/mO2etvBfGdJ3ph85ypvoiMGKd7QO9M0th5dPTK+STnheNv6qyE8Dym93NngsoRdeDUmz0nxqTTeuUnAVQEF2bfefFto7RMQW8+GBB9XB7WPxcseG0vyxuQ78e4LbS4AfImNQ4rA3Rlp92DuT2kE36XRDVMQX/sfmZvxcLKRTdcY0wq4is9EC4+H8CcgZeDpn/sJHqCgC4RnsbHezgYaAVfq9rHNZXYmb0wrIa5SxHdIYvOQPEjwkNir1tugtQAKdrGVOAmW5I0xOUs8veM3ckl5fBoKL3gGIO0ecDq0rGHTNcaYnCb+w8B3KEQXgbRJyabc+cSSvDEm54m4IUtrzTvNpmuMMSaPWZI3xpg8ZkneGGPymCV5Y4zJY3bh1RizRaoRCL0N0eXg2RPx9HY6JLOVLMkbY5qlkR/QNb8Crdq4kYf6BiNl/0TEUki2s+kaY0yzdN2FEFsBWkO8ZHEAgm+jtY87HVpW0NCnxFafSGxZf2IrhxKrfQrV7Nmi0JK8MaZJGl0Bka9pVBuGANROdCKkrKLhOeiaX8drAhGA6BKo/Ctac7/ToW1kSd4Y0zQN0XSaCGUykqykVXfQuApmHdTcjWp2HB9L8saYprk7g7siSYcXbENtiHxFfDvFBhSILst0NElZkjfGNElEkLa3bSr+BYl69F2Q4rMcjS0ruLs10RGDLKmhY5fGjTHNEu/eUDEdrXsGoj8g3v3AfzgiDUv8tj5Scj66dhyb76Hrh6KTEGm4G5czLMkbY7ZI3B2QkrOdDiPriG9/tO1tUHU9xJaD+KHoNKTkAqdD28iSvDHGtICrcDjqP4z42bwPkeyaBbckb4zJGxpZgFbfA+F5UNATKTkb8fRJ+7jxnbayY3qmIUvyxpi8oOEv0DWngQaAGES/Q4Mzod39iG8/p8NzTHZ9rjDGmO2klTeA1rLpxi0FAmjlNQ5G5TxL8saY/BD+Inl79Nt4gbVWqkVJXkRuEZGvRORzEXlORMrq9V0hIgtE5GsRGdHyUI0xphmu0uTtUgi4MxpKNmnpmfx0oJ+q9ge+Aa4AEJE+wElAX+Bw4G4Rab1H2ZgsopFviVVeT2zteWjtZFQDW35RLig6HfA3aPRD4SmJC6OtU4uSvKq+qps+B70PdEn8PBKYqKpBVV0ILAD2bclYxpiW08Br6KrRUPsYBF9FK69HV41GY9VOh9ZiUvxbKBoD+EBK4t8Lj0LaXOR0aI5K5eqaM4FJiZ87E0/6GyxJtDUiIuOAcQBdu3ZNYTjGmPpUw+j6y9n87sw6iC5Ba/+DlJzjVGgpIeJCSq9GSy6A6GJwd0Zc7ZwOy3FbPJMXkddEZE6Sr5H1nnMVEAE2FJhO9tkoaYFlVb1PVQeq6sAOHbKj1oMxeSnyDRBN0hGEwNRMR5M24ipFPP0swSds8UxeVQ9rrl9ExgJHA8N0U6X8JcDO9Z7WBfhpe4M0xqSAFG3c2alxX0lmYzEZ09LVNYcD/wccq6q19bqmACeJiE9EegA9gQ9bMpYxpmWkoAe4d6bx//aFSNGpToRkMqClq2vuBNoA00XkUxG5F0BV5wKTgXnANOBc1aZOIYwxmSLt7gFXJ5Di+BdeKDoB/Ec4HZpJkxZdeFXV3Zrpux64viXvb4xJLSnoCh1mQHg2RFeBd2/E3cnpsEwaWe0aY1oZERd4bUVza2FlDYwxJo9ZkjfGmDxmSd4YY/KYJXljjMljluSNMSaPyaabVJ0nIiuBRU7HkQYVwCqng8hidnyaZ8eneXZ8oJuqJq0Lk1VJPl+JyGxVHeh0HNnKjk/z7Pg0z45P82y6xhhj8pgleWOMyWOW5DPjPqcDyHJ2fJpnx6d5dnyaYXPyxhiTx+xM3hhj8pgleWOMyWOW5NNERG4Rka9E5HMReU5Eyur1XSEiC0TkaxEZ4WScThGRMSIyV0RiIjKwQV+rPz4Q35QncQwWiMjlTseTDUTkIRFZISJz6rWVi8h0EZmf+G77/tVjST59pgP9VLU/8A1wBYCI9AFOAvoChwN3i4jbsSidMwc4DnizfqMdn7jE73wXcATQBzg5cWxau0eI/7uo73Jghqr2BGYkHpsES/Jpoqqvqmok8fB94vvcAowEJqpqUFUXAguAVlfcW1W/VNWvk3TZ8YnbF1igqt+pagiYSPzYtGqq+iawpkHzSGBC4ucJwKiMBpXlLMlnxpnA1MTPnYEf6vUtSbSZODs+cXYctt6OqroUIPF9B4fjySq2M1QLiMhrQMckXVep6guJ51wFRIDHN7wsyfPzch3r1hyfZC9L0paXx2cL7DiYlLAk3wKqelhz/SIyFjgaGKabbkhYAuxc72ldgJ/SE6GztnR8mtBqjs8W2HHYestFpJOqLhWRTsAKpwPKJjZdkyYicjjwf8Cxqlpbr2sKcJKI+ESkB9AT+NCJGLOUHZ+4WUBPEekhIl7iF6OnOBxTtpoCjE38PBZo6lNiq2Rn8ulzJ+ADposIwPuqeraqzhWRycA84tM456pq1ME4HSEio4F/AR2Al0TkU1UdYccnTlUjInIe8ArgBh5S1bkOh+U4EXkSGAJUiMgS4C/AjcBkEfkNsBgY41yE2cfKGhhjTB6z6RpjjMljluSNMSaPWZI3xpg8ZkneGGPymCV5Y4zJY5bkjTEmj1mSN8aYPPb/YJJPs1ik4B8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate gaussian data\n",
    "from sklearn.datasets import make_blobs\n",
    "m = 50\n",
    "n = 2\n",
    "centers = [[-10.0,10.0], [8.0,-10.0]]\n",
    "\n",
    "np.random.seed(0)\n",
    "X, y = make_blobs(n_samples=m, centers=centers, n_features=n, cluster_std = 6.0)\n",
    "y = y.reshape((-1,1))\n",
    "y = 2*(y-1/2)\n",
    "\n",
    "print(\"Shape of the X : \", X.shape)\n",
    "print(\"Shape of the y : \", y.shape)\n",
    "\n",
    "plt.scatter(X[:,0], X[:,1], c=y[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the dual QP parameters\n",
    "def initialize(c=1):\n",
    "    global C, A, b, Q_sqrt, Q, p, v0, eps\n",
    "    C = c\n",
    "    Im = np.eye(m)\n",
    "    A = np.vstack((Im, -Im))\n",
    "    b = C * np.vstack( (np.ones((m,1)), np.zeros((m,1))) )\n",
    "    Q_sqrt = X * y\n",
    "    Q = Q_sqrt @ Q_sqrt.T\n",
    "    p = - np.ones((m,1))\n",
    "    v0 = 0.5 * C * np.ones((m,1)) # stricly feasible\n",
    "    eps = 1e-6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_w(v):\n",
    "    return np.sum(Q_sqrt * v, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can plot the gap $f(v_t)-f^*$ for different values of $C$ and check the impact on $w$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Barrier method for C = 0.001\n",
      "Barrier method for C = 0.01\n",
      "Barrier method for C = 0.1\n",
      "Barrier method for C = 0.5\n",
      "Barrier method for C = 1\n",
      "Barrier method for C = 5\n",
      "Barrier method for C = 10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAE9CAYAAABdmIXpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de3zV5Znv/c8FQZEE5JgoLE4pcoZmZhLS2sJmurEwKdXap0qwxgEUnOfRjhtm2MN0+iqwZ6bFZzN72orblscT9Cmg6IDWQdS2OlJrIdhBOYsFLAGUkzaEhEPCvf9YKxgha607yS/r+H2/Xrxg/XKt+3dl4Stc3r/7vm5zziEiIiIi6aFDshMQEREREX8q3kRERETSiIo3ERERkTSi4k1EREQkjah4ExEREUkjKt5ERERE0khOshNIlN69e7tBgwYlOw0RERGRuN56660Tzrk+zX0ta4q3QYMGsXXr1mSnISIiIhKXmb0f7Wt6bCoiIiKSRlS8iYiIiKQRFW8iIiIiaSRr1ryJiIhI6124cIGqqirOnj2b7FQySufOnQmFQnTq1Mn7PSreREREJK6qqiq6du3KoEGDMLNkp5MRnHOcPHmSqqoqBg8e7P0+PTYVERGRuM6ePUuvXr1UuAXIzOjVq1eLZzNVvImIiIgXFW7Ba81nquJNREREJI2oeBMRERFJIyreAvTxz3/Pxz//fbLTEBERyWhr166ltLSUsWPHMmTIEBYvXtym8TZu3MiwYcMYMmQIS5YsaXFctOuzZs0iPz+f0aNHtym/y6l4C9D5I2c4f+RMstMQERHJWCtWrODBBx/k2Wef5Z133mHbtm106dKl1eM1NDRw33338eKLL7Jr1y5Wr17Nrl27vONivX/GjBls3Lix1blFo1YhAbtwtIZjP3knZkyXoj7klV6foIxEREQyQ3V1NfPmzaOyspJQKARAXl4e8+fPb/WYW7ZsYciQIRQWFgJQXl7Oc889x8iRI73iJk6cGPX9EyZM4ODBg63OLRoVbwHqUtSH2jgxF47WUAsq3kRERFpo3bp1lJaWXiqUYhk/fjynT5++4vrSpUuZNGnSpdeHDx+mf//+l16HQiE2b958xfuixfm+P0gq3gKUV3p93KIs3qyciIhIqlv8853sOlId6Jgj+3Zj4VdHxYzZuXMnRUVFXuNt2rTJK845d8W15tp3RIvzfX+QVLyJiIhIWsjNzaWurs4r1nfmLRQKcejQoUuvq6qq6Nu37xXvixbn+/4gqXhLAp91caC1cSIikprizZC1l7KyMsrLy5k7dy4FBQWcO3eOlStXMnv27CtifWfeSkpK2LdvHwcOHKBfv36sWbOGVatWeccNGzbM6/1B0m7TBOtS1IdO1+fFjbtwtIbabccTkJGIiEh6KCkpYdGiRUyePJkxY8ZQVFTEsWPH2jRmTk4Oy5YtY/LkyYwYMYLbb7+dUaM+KU7Lyso4cuRI1LhY758+fTqf//zn2bt3L6FQiMcee6xNuTay5p7VZqLi4mK3devWZKfhrXFmLv/esUnOREREBHbv3s2IESOSnUZGau6zNbO3nHPFzcVr5k1EREQkjah4C9CrTy7n1SeXJzsNERERyWDasBCgY+/vT3YKIiIikuE08yYiIiKSRtJy5s3MOgD/CHQDtjrnViQ5JREREZGESJnizcweB6YCx5xzo5tcnwL8EOgIPOqcWwLcAvQDTgFVSUg3quMHD/DU4gUxY0Z8YSJjJ01JUEYiIiKSSVKmeAOeBJYBKxsvmFlH4GHgJsJFWqWZPQ8MA950zv3EzJ4Bfpn4dK804gsT48YcP3gAwKt40yH3IiIicrmUKd6cc6+b2aDLLo8D3nPO7QcwszWEZ90OAecjMQ2JyjGesZOmxC3K4s3KNdIh9yIiItKcVN+w0I9wodaoKnLt34DJZvYQ8Hq0N5vZHDPbamZbjx9Pr9MK8kqvJ//esTF/+ZzUICIikmnWrl1LaWkpY8eOZciQISxevLhN423cuJFhw4YxZMgQlixZ0uK4WbNmkZ+fz+jRo6O+N0ipXrxZM9ecc67WOXe3c+5bzrmHo73ZObfcOVfsnCvu06dPO6YpIiIiibBixQoefPBBnn32Wd555x22bdtGly5dWj1eQ0MD9913Hy+++CK7du1i9erV7Nq1q0VxM2bMYOPGja3OoaVS5rFpFFVA/yavQ8CRJOWSknTIvYiIZIvq6mrmzZtHZWUloVAIgLy8PObPn9/qMbds2cKQIUMoLCwEoLy8nOeee46RI0d6x02YMIGDBw+2OoeWSvXirRK4wcwGA4eBcuCO5KaUOnzWxYHWxomISGZYt24dpaWllwqoWMaPH8/p06evuL506VImTZp06fXhw4fp3/+TeaJQKMTmzZuveJ9vXCKkTPFmZquBiUBvM6sCFjrnHjOz+4GXCLcKedw5tzOJaaaUvNLrvQoyn5k5ERERby8ugA+2BzvmdWPgL6KvNwPYuXMnRUVFXsNt2rTJK845d8U1sytXbfnGJULKFG/OuelRrm8ANiQ4HREREUkxubm51NXVecX6zryFQiEOHfpkb2RVVRV9+/a94n2+cYmQMsVbNvFp5AvBNvPV2jgREQlMnBmy9lJWVkZ5eTlz586loKCAc+fOsXLlSmbPnn1FrO/MW0lJCfv27ePAgQP069ePNWvWsGrVqlbHJUKq7zbNOCO+MJE+gwbHjTt+8AC733gtkHt2Kerj1VbkwtEaarelV0sVERHJHiUlJSxatIjJkyczZswYioqKOHbsWJvGzMnJYdmyZUyePJkRI0Zw++23M2rUqEtfLysr48iRIzHjpk+fzuc//3n27t1LKBTisccea1NO8Vhzz3AzUXFxsdu6dWuy0/DWODM3bWHi/u+mcWYu/96xCbuniIikh927dzNixIhkp5GRmvtszewt51xxc/GaeRMRERFJI1rzJp+i81RFRERSm4o3uUTnqYqIiKQ+FW9yiU/fOPWMExERSS6teRMRERFJI5p5C9Cmp98FYPztQ5OcSftSzzgREZHkUfEWoBOHagIdz6eZb5CNfH3oPFUREZHkUvGWokZ8YWLcmOMHDwAktHjTeaoiIiLJpeItRY2dNCVuUeZzxJaIiEimWbt2LUuXLqWuro7a2loqKipYuHBhq8fbuHEjDzzwAA0NDdxzzz0sWND8v6+zZs3ihRdeID8/nx07drT6fm2l4i1gJ6pqWPcvv4sZM3RcAaPG90tQRsmjnnEiIhK0FStW8NBDD7F+/XpCoRA1NTU88sgjrR6voaGB++67j1deeYVQKERJSQk333wzI0eOvCJ2xowZ3H///dx1111t+RbaTMVbgIaOK4gbc6IqvC4uqOItGYfc+1DPOBERCVp1dTXz5s2jsrKSUCgEQF5eHvPnz2/1mFu2bGHIkCEUFhYCUF5eznPPPdds8TZhwgQOHjzY6nsFRcVbgEaN7xe3KIs3K9cSPuviIHXXxmldnIiItMS6desoLS29VGjFMn78eE6fPn3F9aVLlzJp0qRLrw8fPkz//v0vvQ6FQmzevDmYhNuJirc05rMuDrQ2TkREgvXglgfZc2pPoGMO7zmcvxv3dzFjdu7cSVFRkdd4mzZt8opzzl1xzcy83pssKt4kqdQzTkREfOXm5lJXV+cV6zvzFgqFOHTo0KXXVVVV9O3bt+3JtiMVb1kiFdfGqWeciEh6ijdD1l7KysooLy9n7ty5FBQUcO7cOVauXMns2bOviPWdeSspKWHfvn0cOHCAfv36sWbNGlatWhV06oFS8RagD773PQCu+/a3k5zJp6Xq2jj1jBMRkZYoKSlh0aJFTJ48mYaGBurr67nzzjvbNGZOTg7Lli27NOasWbMYNWrUpa+XlZXx6KOP0rdvX6ZPn85rr73GiRMnCIVCLF68mLvvvrut31bLc074HTPYud3BPv8PSiasjVPbERERAaioqKCioiLQMcvKyigrK2v2axs2bLj059WrVwd639ZS8SYpT21HREREPpG2xZuZ5QKvAwudcy8kOx9pP2o7IiIi8omUKd7M7HFgKnDMOTe6yfUpwA+BjsCjzrklkS/9HfB0whMNgM8pDJA9JzGIiIiIv5Qp3oAngWXAysYLZtYReBi4CagCKs3seaAvsAvonPg028bnFAYI/iSGbKC2IyIikg1Spnhzzr1uZoMuuzwOeM85tx/AzNYAtwB5QC4wEqgzsw3OuYsJTLfVfE5hgGBPYsgGajsiIiLZImWKtyj6AYeavK4CSp1z9wOY2QzgRLTCzczmAHMABgwY0L6ZSlKp7YiIiGSLVC/emjuf4tI5Fs65J2O92Tm3HFgOUFxcfOX5F3IFn2a+iT7kXkRERD6R6sVbFdC/yesQcCRJuXg5u2cP71fcFTOm29Sp9Jh2e4Iy8ufTzDcZh9yLiIjIJ1K9eKsEbjCzwcBhoBy4I7kpRddt6tS4MWf3hBv5pmLx5tPMN5Ub+YqISHZYu3YtS5cupa6ujtraWioqKli4cGGrx9u4cSMPPPAADQ0N3HPPPSxY0Py/dYMGDaJr16507NiRnJwctm7d2up7tkXKFG9mthqYCPQ2syrC/dseM7P7gZcItwp53Dm3M4lpxtRj2u1xi7J4s3IiIiIS3YoVK3jooYdYv349oVCImpoaHnnkkVaP19DQwH333ccrr7xCKBSipKSEm2++mZEjRzYb/+qrr9K7d+9W3y8IKVO8OeemR7m+AdjQ3NdEREQke1RXVzNv3jwqKysJhUIA5OXlMX/+/FaPuWXLFoYMGUJhYSEA5eXlPPfcc1GLt1SQMsWbXMmnma8a+YqISLZYt24dpaWllwqtWMaPH8/p06evuL506VImTZp06fXhw4fp3/+T5fWhUIjNmzc3O6aZ8eUvfxkz495772XOnDmt+C7aTsVbivJp5qtGviIikgwffO97nNu9J9Axrx4xnOu+/e2YMTt37qSoqMhrvE2bNnnFOXdlMwqz5ppdwBtvvEHfvn05duwYN910E8OHD2fChAle9wmSircU5dPMV418W0cnMYiIpKfc3Fzq6uq8Yn1n3kKhEIcOfdJStqqqir59+zY7ZuP1/Px8br31VrZs2aLiTaS96SQGEZG2izdD1l7KysooLy9n7ty5FBQUcO7cOVauXMns2bOviPWdeSspKWHfvn0cOHCAfv36sWbNGlatWnVF3JkzZ7h48SJdu3blzJkzvPzyy3z3u99t8/fUGirepMV8GvlCajbz1UkMIiLpq6SkhEWLFjF58mQaGhqor6/nzjvvbNOYOTk5LFu27NKYs2bNYtSoUZe+XlZWxqOPPsrZs2e59dZbAaivr+eOO+5gypTk/Bun4i0JfBr5Qmo28/Vp5Atq5isiIu2joqKCioqKQMcsKyujrKys2a9t2PBJw4u333470Pu2loq3BPNp5Aup28zXp5EvqJmviIhIe1HxlmA+jXxBzXxFRESkeR2SnUAmeXDLgzy45cFkpyEiIiIZTDNvAdpzKtieNz58GvmCmvmKiIhkChVvacynkS+oma+IiEgmUfGWxnwa+YKa+baWTzNfNfIVEZFEU/Em0gyfZr5q5CsiIsmg4k3ajU8z31Rs5At+zXzVyFdERJJBxVsK82nmm4qNfMGvma8a+YqISGusXbuWpUuXUldXR21tLRUVFSxcuLDV482aNYsXXniB/Px8duzYEWCm7UPFW4ryaeabqo18wa+Zrxr5iohIS61YsYKHHnqI9evXEwqFqKmp4ZFHHmnTmDNmzOD+++/nrrvSo8eqircU5dPMV418RUQkm1RXVzNv3jwqKysJhUIA5OXlMX/+/DaNO2HCBA4ePBhAhomh4i1ge0/tZebGmTFjygrLuG3obQnKSEREJDOsW7eO0tJSCgsL48aOHz+e06dPX3F96dKlTJo0qT3SSxgVbwEqK2z+UNum9p7aC5Dw4k3NfEVEJCibnn6XE4dqAh2zd/88xt8+NGbMzp07KSoq8hpv06ZNQaSVklS8Bei2obfFLcrizcq1BzXzbT8+veBA/eBERIKQm5tLXV2dV6xm3iStqZlv+/DpBQfqBycimSfeDFl7KSsro7y8nLlz51JQUMC5c+dYuXIls2fPviJWM28pyMy+BnwFyAceds69nOSUpBV8esFBavaD8+kFB+oHJyISlJKSEhYtWsTkyZNpaGigvr6eO++8s83jTp8+nddee40TJ04QCoVYvHgxd999dwAZt4+UKt7M7HFgKnDMOTe6yfUpwA+BjsCjzrklzrn1wHoz6wEsBbKyePPpBQep2Q/OpxccZEY/OB21JSISjIqKCioqKgIdc/Xq1YGO195SqngDngSWASsbL5hZR+Bh4CagCqg0s+edc7siId+JfD3r+PSCg9TtB+fTCw7Svx+cjtoSEZEgpVTx5px73cwGXXZ5HPCec24/gJmtAW4xs93AEuBF51xWLtby6QUHmdEPTkdtiYiIhKVU8RZFP+BQk9dVQCnwLWAScK2ZDXHO/fjyN5rZHGAOwIABAxKQqh+fXnCQnH5wPi1FEt1OREdtiYiIfCIdijdr5ppzzv0I+FGsNzrnlgPLAYqLi1075NZiPr3gIDn94HxaiiSjnYiO2hIREflEOhRvVUD/Jq9DwJEk5dJmPr3gIDn94HxaiqidiIiISHKlQ/FWCdxgZoOBw0A5cEdyUxIJnhr+ioiIjw7JTqApM1sNvAkMM7MqM7vbOVcP3A+8BOwGnnbO7UxmniJB61LUh07X58WNu3C0htptxxOQkYiIpKqUmnlzzk2Pcn0DsCHB6YgkjBr+ioiIr5SaeRMRERGJZ+3atZSWljJ27FiGDBnC4sWL2zTeoEGDGDNmDEVFRRQXFweUZftJqZk3ERERkVhWrFjBQw89xPr16wmFQtTU1PDII4+0edxXX32V3r17B5Bh+1PxJtJK7/xiI7vfeC2h9xx97nMAvLp4VdSYgvr+9Gnwa+XSpVt3cnv0jB+nTRIikgKqq6uZN28elZWVhEIhAPLy8pg/f36SM0ssFW9Z4HjdcU7WnWRRAO1H/vTgRK6tdbxf8YOYce/nDOVwTmGb79cop1cvcvr0ifr1E4dquKbrVXHHCbLgqtq1A4DQyNFxIoOVe7HbpSKuOdde7AXAHzucjDnOhbNnqeXjuMWbju4SkVSxbt06SktLKSyM/+/L+PHjOX369BXXly5dyqRJkz51zcz48pe/jJlx7733MmfOnMBybg8q3tLY2nfXsmF//H0cX6k+6DXeZ3/zISPeOhEz5kTX0vAfOsUe63BOIdUdetLt4imve8dysbaWeohZvF041wCcjzvW7jde4/jBA/QZNLjNeYVGjk74kVw1m49Su+04XYg9td+lqA+h0vExYxobG0+792sx47RJQkQu9+qTyzn2/v5Ax8wfWMifz4hdNO3cuZOioiKv8TZt2uR97zfeeIO+ffty7NgxbrrpJoYPH86ECRO8359oKt7S2Ib9G9h7ai/Deg6LGde1U1d6XdOLJ6Y8ETPu/Z/dxdljJ+g8fHjUmF90zKHumr78bsi4mGPVVNWQH8rj1r9pe2HzfsVdcBYG/k309n4P39OR87UfxD1pobFwm7ZwSZvzSgbfXakiIpkoNzeXuro6r9iWzLz17dsXgPz8fG699Va2bNmi4k3az7Cew7yKMl+dhw9n4E9XRv36nh9/j16HugDXxRyndyjP67itoOT1GstH5xo4cagmZpx17MO1+Z9NUFap7/jBA3EL3tHnPkfXDj3jzsBpXZxI9og3Q9ZeysrKKC8vZ+7cuRQUFHDu3DlWrlzJ7Nmzr4j1nXk7c+YMFy9epGvXrpw5c4aXX36Z7373u0GnHigVbyks3gH2PrNujc7u2ROewYoTE2vWDeDEoH2cGLSPBVMSe8hFvPwH5wylS7+vxXy0CuGzWetq4zfDzQYjvjDRK27/x+9Q2H0snekeNUbr4kQkEUpKSli0aBGTJ0+moaGB+vp67rzzzjaN+eGHH3LrrbcCUF9fzx133MGUKYlbDtMaKt6C9GJkBuMv2v5IzucA+2E9h3nFdZs61euenYcP945NJJ+cCnY8z8D6dxm4JPqsIehs1qbGTpritV7vqcUL2MFvY66N07o4EUmUiooKKioqAhuvsLCQt99+O7DxEkHFW5A+2B7YUL4H2PvoMe12eky7PZCxksEn/3iziiIiIplCJyyIiIiIpBHNvAXtg+3wxFcSd78x34DitvdvE4kl3sYG300NoI0NIiJtpeItSGO+kdj7NT6mTXDxFm8jRaOywrLAHv1K8vhsbPDZ1ADa2CAiEgQVb0EqnpnYQiqRM3wRPhskIFzgASreMoDPxgafTQ2gjQ0iIkFQ8SYt4ruRwmdmTjKLesaJiCSGijfJGF697DpPIadXrwRllD3UM05EJHFUvElG8O1Pd7G2llMdenr1exs6roBR4/u1NbWsoJ5xIiKJo+JNMoJvL7v3Z/4Thy92AfJjxp2oCh+zpeJNRCT1rF27lqVLl1JXV0dtbS0VFRUsXLiw1ePNmjWLF154gfz8fHbs2PGpr23cuJEHHniAhoYG7rnnHhYsiL08JBFUvElWGVj/bvhXjEPuQScxiIikqhUrVvDQQw+xfv16QqEQNTU1PPLII20ac8aMGdx///3cddenl940NDRw33338corrxAKhSgpKeHmm29m5MiRbbpfW6l4S3e+feXUD05ERNJcdXU18+bNo7KyklAoBEBeXh7z589v07gTJkzg4MGDV1zfsmULQ4YMobCwEIDy8nKee+651C/ezCzXOXfGzPKcczWJSEo8+faVS+F+cOoFl33U8FdEWmvdunWUlpZeKqZiGT9+PKdPn77i+tKlS5k0aZLX/Q4fPkz//v0vvQ6FQmzevNk/4XbiM/PWw8xmAu8BG9s5Hy9mlgv8b+A88Jpz7mdJTik5fPvKpWg/uGT1gtOu1ORRw1+RzPDxz3/P+SNnAh3zqr65dP/qZ2LG7Ny5k6KiIq/xNm3a1OacnHNXXDOzNo/bVj7F238FZgCPm1m+c+5YeyRiZo8DU4FjzrnRTa5PAX4IdAQedc4tAb4OPOOc+7mZPQVkZ/GWwnz6wSWjF1zQu1K1I7Vl1PBXRNoiNzeXuro6r9ggZt5CoRCHDh269Lqqqoq+ffv6JduOfIq3LcAsoH97FW4RTwLLgJWNF8ysI/AwcBNQBVSa2fNACIg8C6ShHXOSDBPkrlTtSG0/avgrktrizZC1l7KyMsrLy5k7dy4FBQWcO3eOlStXMnv27Ctig5h5KykpYd++fRw4cIB+/fqxZs0aVq1a1eZx2ypu8eac2x35Y7v+b65z7nUzG3TZ5XHAe865/QBmtga4hXAhFwK2AR3aMy/JTj67UrUjtX2o4a+IRFNSUsKiRYuYPHkyDQ0N1NfXc+edd7Z53OnTp/Paa69x4sQJQqEQixcv5u677yYnJ4dly5Zdut+sWbMYNWpUAN9J23jtNjWzHzrnHjCza5xzfvOVwegHHGryugooBX4ELDOzrwA/j/ZmM5sDzAEYMGBAO6Ypkhg7Nx3m3S0fJjuNZgX1CLklDX9fPbiaPlcPjhoz+vznyD1wlp3fXt/mvFqiS7fu5PboGT9Os4IiLVZRUUFFRUWgY65evTrq18rKyigr8zvXO1F8W4X818jvvwb+rJ1yaU5zqwKdc+4MEHfBlHNuObAcoLi4+MpVhyJtUH/8OPUnT/J+xQ9ixr2fM5TDOfF3Rvk41fE6APreEHsxf6Il4xGyzwzd8Y6H2z+Ry1w4e5ZaPo5bvGlWUERay7d422hmbwLXmdks4G1gp3PubPulBoRn2vo3eR0CjrTzPUW81J88SbW7lt90jj1L1Fhw9Wz4oM337P7xuxR8uJUbzpyPGddt6lSvtX1BScYjZN8ZukRrXKunDRci0l68ijfn3N+aWSHwGjAYuBkYZWbngR3OuWntlF8lcIOZDQYOA+VA7Nb4IgGJ11Ik/0Q3CN1I5+HDY47Tl8ZHil9qc04fPfU01S/ELtxqKyuprayk+oUX4o6X6CJPRETazvuEBefcfjOb5Jx7t/GameUBo2O8zZuZrQYmAr3NrApY6Jx7zMzuB14i3CrkcefcziDuJ+lj7btr2bB/Q0Lv+dnBJxhxCji1J3pQB+gzJpdJf524/5/w2S0bLvDiF24q8tqPdstKpnLOpUSfs0zSXC+5eFp0PFbTwi3yugb4bYvv2vzY06Nc3wAk9l9uiW7rE7D9mWDGssjC+zhNhDfYh+ztCMP6jAnmvh7evrGAt28siBmz99RehvU8hl+3oMTxbYcSZJGnpsaf0G5ZyVSdO3fm5MmT9OrVSwVcQJxznDx5ks6dO7fofTrbVD7hU5i9/+vw7wO/GMgt93KemRZ79+Red5ZhDZ15YsoTgdwzKMloMhykIIu8i7W11AeVWJpryW7ZeM2ItS5OUkkoFKKqqorjx48nO5WM0rlz50vntPpqVfFmZl91zkVt0SFpavsz4XNQr4sxwzXwi4Edcl/27lrweBw67IPtlP3xo/jHfAWUl3yaT5H3m9nJb1opIu2rU6dODB4cvTWPJE5rZ97+mRj91SSNXTcGZv57Qm7lc4QW4Dcj+EHkwA0Vb0lzsbY27pmxoPVzIiJt1driTQ+7080H2+PPXMWbdUuW4pnxi7J431s72Xtqr9fj07LCMr9CNU3l9OoVfmwap3nQ2T3hDSAq3kREWq+1xZsa3qaTMd/wi7tujH+sUFbo13F776m9AJldvPXpQ06fPjGPEwO8ZuZERCQ2bVjIBj4zV5nAZ3YRAlsb5/vYd+bGmV4zdJk+OydXitdSxLedCKiliEg2UfEmmcF3xjAJa+N8ZuiyYXZOPs2npYhPOxFQSxGRbNPa4i01T8aW7OU7u5iEtXE+M3S+s3OgGbpM4dNSxKedCKiliEi2aVXx5py7KehERLKZ1s+JiIgvPTYN0OKf72TXkeq4cbcU9eOO0gEJyEjSRUvWz6WqE1U1cQ+oP9t5CvlVv4E4GxfUTuQTOmpLRC6n4i3Bdh0NF3cq3iSTDB0X+yixRqc6XsepgV/nWMMHUWMunq6m4Mcvc4POXNVRWyLSrNaesJALnHXONQScT1pb+NVRcWOm/eRNdh2tZl5DxfQAAB7vSURBVNpP3gzknprFk1Qwanw/Ro3vFzdu56bDvLvlQ4hRZBw/cIrjXbtxw9mNMcfKhp5xOmpLRJrjVbyZWQegHPgmUAKcA642s+OED41f7pzb125ZZpBbiuL/A+dLs3iSbnyKvPCj155ePePO7tkTt3dcps/OiUj28Z15exX4BfD3wA7n3EUAM+sJ/DmwxMzWOef+//ZJM3PcUTogsGIrqNm7rJPgfnDScj7r5+r7fZOCnEoG1r8bNSYbZudEJPv4Fm+TnHMXLr/onDsFPAs8a2adAs1MvPg+gtXj1YgU7gcnYb7r5z4+dw05o8uY8DffiRqjEx1EJBN5FW/NFW4AZjbPOfe/Ii8Lgb1BJSbx+T6C1ePVJlrSD85nhi4Js3OZflqD7/q5eDNzIiKZqrUbFroD/woMM7OzwDvA3YCmKRLI9xGsHq+2gs8MnU5rSAs+6+JAa+NEJH20tknvx8BMM5sMnADGAv8WZGIiSeUzQ5fCpzVki3hr43zWxYHWxolIemlR8WZmvwH+wTn3KoBz7qXIl94KOjERab1sOGrLZ22cz7o40No4EUkvLZ15mwMsNrPvAN9xzul5nGS3FNy56nvU1tYPt7L1w61s2L/Ba8xUK/L824740eNVEUkXLSrenHM7gP/LzP4U+B9mBuEiblt7JCfB0a7UdpCiO1d9j9pa++5ar8KtJUVeqhl2KnwM88yND8WMK7quLz2u+hKdOlwVM+7i6WpYd4oOv1gVWI7x9KvfH/+x74Vqcnr1CuR+7/xiI7vfeC1uXEF9f/o0BNe3sku37uT26Bk7Rsd7iQBgzrmWv8msGzCC8AaFe5xzCT9my8y+BnwFyAceds69HCu+uLjYbd26NSG5pZpVm//Ac9sOx43bdbSakdd346l7P5+ArLJI487V68YkO5MrecwI+hZ5qWjYppu45o89qbv2VMy4rievA6DvDdFPfgCoP36c+pMnA8svnuoOPel28RQ3xjlt4lenDlN9zVV0z7k6asyfXnc7Xa8u4EzH0zHHOnfmDABX5+bGjLv2YrhY/GOHtn8eF86epVPnzvQZWBg95mgNna7PI//esW2+n0g6MLO3nHPFzX2tpWvefgXcAJwFdkV+zWhFQo8DU4FjzrnRTa5PAX4IdAQedc4tiTaGc249sN7MegBLgZjFWzbTrtQk852hSzTPGUHfmbxUtDO38Tiu/Jhxe9nDyf4Hue+vvp2YxDyFH/vmxz1tYviSf+LdbbH/5/Toqbeh52fp0KVLzLirc3O9ZsEgPBMWKh0fNy6epxYvANDxXiKeWjpj9rfAbudcXRvv+ySwDFjZeMHMOgIPAzcBVUClmT1PuJD7/mXvn+WcOxb583ci7xNJTb695RLNt5cdpO1pE7494+I9Vk11n1/wHeLNl79fcRecPMrAH6yMEykiqc73bFNzYVFX/zbG+IznnHvdzAZddnkc8J5zbn9kvDXALc657xOepbvifsAS4MVYeYlIFCm6Zi9ZsmGHrohkBt+Zt1+Z2b8Bzznn/tB40cyuAr4I/CXh80+fbEMu/YBDTV5XAaUx4r8FTAKuNbMhzrkfXx5gZnMI75BlwAAtwg+C7/o5X9ogkUQtOW0iw/nu0FUDZBFJBb7F2z6gAVhnZtcDHwOdCT/SfBn41wB2nFoz16LO5DnnfgT8KNaAzrnlwHIIb1hoU3ZZIt6u1M0Hwgu/SwfHXw/jcy/QsV1pIUWPCguK77q+dG+ArHYoIpnBt3i70Tk3x8zuAQYAfYC6yEkLQakC+jd5HQKOBDi+xOFzVmrp4J6BzZZN+8mbamGSDnwer77/6/Cv7c/4jZemRR6k79my3aZesfqkWTptQiT1+RZvL5nZm0ABcBfwNrAz4FwqgRvMbDBwGCgHYm+xkkD57koNik+xCJqhSzqfx6tbn/Ar3NJ8/Vw6ny3bY9rtXgWZTpsQSX1exZtz7m/MrBB4DRgM3AyMMrPzwA7n3LSW3NTMVgMTgd5mVgUsdM49Zmb3Ay8Rfhz7uHMu6AJRUohamGSQLFk/l4yzZeOd39po6LgCr521qer4wQOXWoY0Z/S5z9G1Q0+vliFq5iuZzrtViHNuv5lNcs5davVtZnnA6BhvizbW9CjXNwDp2Q1U2pXP41U9Wk0TGb5+Lkg+57dCuMADAiveEr02bsQXJsaN2f/xOxR2H0tnYjdSvnC0hlpQ8SYZraXHY7172esa4LeBZiRyGZ/Hq3q0miZ81s+l+aNVCK7tiG+fupac4RpPMtbGjZ00hbGTpsSMeWrxAnbw25iNfEHNfCU7JPxYK5GW8nm8qkeracLn8WqaP1pN97YjWhsnkvpUvEnG0M7VDJLGJz9kS9sREUkeFW+SEbRzNYPo5AcRkZhUvElG0M7VDJIlO1czgc/GBjX8FQmeijcREWkxn40Navgr0j5UvEnW0do4yUQ+/eCC7AXns7Eh6E0N8XrBgX8/OPWCk3Sm4k2yitbGSSby6QcXdC+4RPPpBQd+/eDUC07SnYo3ySpaGyeZyKcfXJC94JLBpxcc+PWDUy84SXcdkp2AiIiIiPhT8SYiIiKSRvTYVETSl85JFZEspOJNJAqfXanakZpEWXJOqojI5VS8iTTDZ1eqdqQmme85qWl81JaISHNUvIk0w2dX6rSfvKmecakuxY/a2ntqb9wzTssKywI7vN6nFxwE2w8uGeL1g/PtBQfqByepyZxzyc4hIYqLi93WrVuTnYZkkFWb/8Bz2w7Hjdt84BQApYN7tndKLaaiMqJxZm7mv7d9rK1PwPZn4oatpYYNdiZmzF7OM4yreMLF7+MWz8pjQzlyYmTcuK6nPwPA6a6/b/M9+xy7SKcOOXQd/SdxY4MqGN9c8k+8uy32z/pueaMYmDuSTmYx43p07g/AR2cPxb2v5XTCruoUM6ZLt+7k9oj/c0AFowCY2VvOueLmvqaZN5FW8u0Z51vkJZoe+17G9/FqPO//Ovz7wC/GDLuNPG5zeTFjZtqH7OU8M+3DNqe19bo/wHW/oNhdHTOu94el9DpR1Ob7AVzEceFifdw4nwbCHz31NNUvvBB3rO6VlYwDupSURI3ZX7uP7TU7447VN28M1+WNiBvnGhoAYhZvF86epZaP4xZvaiAsPlS8ibQz3yIv0dSIuAnfx6s+Bn4xsPVzZe+uhf0bAkgKign2EayPl24aBcDkFbELJZ9HudUvvMDZPXvoPHx4zLguJSV0mzo15tFdA+PerWUajwEb+IOVUWMaH+PGah4MaiAsflS8iWQxrdmL8Nn8kAS3Db0tocVWu3AX485o1r/3derPdeT9ih9EjWks3Ab+NHqBJJItVLyJZCmd8yrtrmMnaLgQN6y+5gIX6+tj/ovUefhwuk2dGmByIulLxZtIltI5r9LuOl5F/gf1vP+rXjHDLtYbHa4yBq7QrJqIj7Qu3swsF3gdWOici7+SVURaRQ2LpTV2/1lv4ATx9ld2uMrI6dIxESmJZISkFG9m9jgwFTjmnBvd5PoU4IdAR+BR59ySOEP9HfB0uyUqImpYLK329o0FvH1jAU9MeSJm3O/mr0hQRiKZIVkzb08Cy4BLc+Rm1hF4GLgJqAIqzex5woXc9y97/yxgLLAL6JyAfEWylm/DYhERSYykFG/OudfNbNBll8cB7znn9gOY2RrgFufc9wnP0n2Kmf05kAuMBOrMbINz7mK7Ji4iItKMs3v2XGoZ0uzXL1ST0yv22j8RX6m05q0f0LSNdRVQGi3YOfcPAGY2AzjRXOFmZnOAOQADBuhxjohIovkcATaMG+lW2y9uv7dUPbbLZxfsxdpaTrqGmMd2gf/RXTqFIbulUvHW3Dklcc/ucs49GeNry4HlED4eq9WZiYhIi5UVlnnF7eq9lZEn4HquixrjcwpDsvSYdnvMpsAAg+64nT9cPB93rP0fv0Nh97F0pnvUGJ3CIKlUvFUB/Zu8DgFHkpSLiIi0kW+T4ZkfFLO34C0WzIh+JqnPKQyprLBjZwo7dmbgwtj78J5avIAd/DbmSQw6hUE6JDuBJiqBG8xssJldBZQDzyc5JxEREZGUkqxWIauBiUBvM6si3KftMTO7H3iJ8A7Tx51z8U8OFhGR9Hf+TOxjtD64HXL7JC4fkRSWrN2m06Nc3wAEcwqziIikh8aiLNbK5PNnEpKKSDpIpTVvIiKSjbpeF/4Vq5mvGvmKXKLiTUQC4XOEFugYLWlevJYiw7iRbmeuZ51HETf0s10YdWf8jRKp6vjBAzFbivi2EwG1FMlUKt5EpM18jtAC2HzgFJsPnOK5bYfjjqcCL3v4tBTZlb+Nkcfgeq6KGXeiphe8fZJRdwaVXXDiNfIFyG84S32PrjFjfNqJgFqKZDIVbyLSZj5HaAGs2vyHuIWbzknNPj4tRWYyk738hgVxzkldN39F/M0PAGO+AcWxmwcHyaeRL0DffQcpHD6cgQ9Hbyni004E1FIkk6l4E5GE0Tmp0u58dqS+/+vwr+3PxI8NqMjzaeQLxJ2ZEwEVbyIikkkaNz/M/MvoMVuf8CvcPtge/j2BM3QiPlS8iUjK0eYHaVfFM/0Ksie+Ei7g4j2ChYQ/hpXspuJNRFKK7+YHrY2TdjfmG35xmqGTBFPxJiIpxXfzg9bGSbtryQxdgOLtSj17oZqcXr0CvaekFxVvIpK2fB6v6tGqpBOfXakXa2upT0AukrpUvIlIWvJ5vBr0o1WfVifJkg1FarxGvgDDTt3ENX/syZJ/WBXIPXte04s+18TZwfrB7Qzt9AtGBTAD1wPo8aXYMR1+7vhjw3me+qvY9xt1zXRycwrY/tc/jT2gGViHuLldjXGNT9yQLuTf99W4cdJ6Kt5EJC35th3x3fzgY/OBUwCUDu4ZyHhByYb1fz6NfAFO9j9IUA8Ua+vroO5k3OLtxLl+wCRG8XRAd45twPla/nBVF4jTsPjox9spyHOYxYpy4TNlY8ZAvXNgcE285Dr04Nx7H8WLkjZS8SYiGct384Ov0sE9U3KGKxvW//k08g1a4yzfgil3xIxb9y+/A3rGbk8SoMJf3UXfPXvonDs8ZtzZPa/SefhRBv50ZfSgxtnCmf8ec6zG47qmLYzePBjg0NynYn5dgqHiTUQylu/mh0yg9X/Zw/e0hs7Dh3vHSnpR8SYikuaSsf5Pksf3tAbJXCreRETSnI4dE8ku8beNiIiIiEjKUPEmIiIikkZUvImIiIikERVvIiIiImlEGxZERESy2Qfb45/PehTIi3PSRMTF2tqYZ7NCuN2Jdsy2XtoWb2bWAfhHoBuw1Tm3IskpiYhIhvE9kqtbTX6kWW90Q8cVMGp8sI2j22zMN/zizp/h+NFPmvVGM6rjn9H12v506BK9GLxYW8vpX71Hj2ktSVSaSkrxZmaPA1OBY8650U2uTwF+CHQEHnXOxWrlfAvQDzgFVLVjuiIikoV8j+Ta1fO3jORzXE/0Y9NOVNUApF7xVjwz/CuOEYe+AkfiD3egdjeF3cfS5zMjosac+/0xuNClJVnKZZI18/YksAy4dGaHmXUEHgZuIlyMVZrZ84QLue9f9v5ZwDDgTefcT8zsGeCXCchbRESyhO+RXDOZyV5eiXmMVrxZuVQ3tj+M7bQdrosd99RR2HF+P9PufTJqjI7QarukFG/OudfNbNBll8cB7znn9gOY2RrgFufc9wnP0n2KmVUB5yMvG9ovWxERkSzXgser1LRvKpJaa976AYeavK4CSmPE/xvwkJmNB15vLsDM5gBzAAYM0JEwIiIireL5eJXNcTY+SCBSqXizZq65aMHOuVrg7lgDOueWA8sBiouLo44lIiIiki5Sqc9bFdC/yesQXssjRURERLJHKhVvlcANZjbYzK4CyoHnk5yTiIiISEpJVquQ1cBEoHdk48FC59xjZnY/8BLhHaaPO+d2JiM/ERGRlojXD863FxykaD+4Fjh+OnY/uJFuLN1yenvtOr16SBfy7/tqkOllhGTtNp0e5foGYEOC0xEREWk1n35wPr3gIIX7wXka0Ze4C54O1u1h0DXD6W4dYwd26MG59z4KLLdMkkobFkRERNKOTz84n15wkCH94PoDM6P32H9q8QJ28Q7TFsbqw69+cLGoeBMRyRK7jlYz7Sdvxo27pagfd5SqvZJIqlLxJiKSBW4p8nsMt+toNYCKN5EUpuJNRCQL3FE6wKsg85mZE5HkSqVWISIiIiISh2beRETkU7Q2TiS1qXgTEZFLtDZOJPWpeBMRkUu0Nk4k9WnNm4iIiEgaUfEmIiIikkZUvImIiIikEa15ExGRVvHZlaodqSLBU/EmIiIt5rMrVTtSRdqHijcREWkxn12p2pEq0j605k1EREQkjah4ExEREUkjKt5ERERE0oiKNxEREZE0ouJNREREJI2oeBMRERFJIyreRERERNJI2vZ5M7MBwDLgBPCuc25JklMSERGJau+pvczcODNmzLBTN9GtJp91//K7mHFDxxUwanz8Rsmp6vjBAzy1eEHMmJFuLN1yenNo7lMx464e0oX8+74aZHopLynFm5k9DkwFjjnnRje5PgX4IdAReDROQTYU+Hfn3E/MbGW7JiwiItIGZYVlXnG7ev6WkXyO6+kZNeZEVQ1A6hZvH2yHJ74S9csjDOjVJ+4wB+v2MOia4XS3jtGDOvTg3HsftSLJ9JasmbcnCc+aXSq6zKwj8DBwE1AFVJrZ84QLue9f9v5ZwH8C/2Bm04CfJiBnERGRVrlt6G3cNvS2uHEzmcleXmHBlDuixsSblUuqMd+IGzK203bGjhoDM5+MGffU4gXs4h2mLYw+jxNvVi5TJaV4c869bmaDLrs8DnjPObcfwMzWALc4575PeJbuU8zsb4GFkbGeAZ5o36xFREQkpuKZ4V+xxJiVEz+ptGGhH3CoyeuqyLVoNgJ/bWY/Bg42F2Bmc8xsq5ltPX78eGCJioiIiCRLKm1YsGauuWjBzrkdQMz5WefccmA5QHFxcdSxRERERNJFKs28VQH9m7wOAUeSlIuIiIhISkql4q0SuMHMBpvZVUA58HyScxIRERFJKUkp3sxsNfAmMMzMqszsbudcPXA/8BKwG3jaObczGfmJiIiIpKpk7TadHuX6BmBDgtMRERERSRup9NhUREREROJQ8SYiIiKSRlS8iYiIiKQRFW8iIiIiaUTFm4iIiEgaUfEmIiIikkZUvImIiIikkVQ621RERESywQfb4YmvxI45CuT1SUg66UbFm4iIiCTOmG/4xZ0/AzXtm0q6UvEmIiIiiVM8M/wrns1xZuaymNa8iYiIiKQRFW8iIiIiaUTFm4iIiEgaUfEmIiIikkZUvImIiIikERVvIiIiImlExZuIiIhIGlHxJiIiIpJGVLyJiIiIpBEVbyIiIiJpRMWbiIiISBox51yyc0gIMzsOvJ+AW/UGTiTgPnIlffbJpc8/ufT5J5c+/+TJ1M9+oHOuT3NfyJriLVHMbKtzrjjZeWQjffbJpc8/ufT5J5c+/+TJxs9ej01FRERE0oiKNxEREZE0ouIteMuTnUAW02efXPr8k0uff3Lp80+erPvsteZNREREJI1o5k1EREQkjah4C4iZTTGzvWb2npktSHY+mc7MHjezY2a2o8m1nmb2ipnti/zeI5k5ZjIz629mr5rZbjPbaWYPRK7r76CdmVlnM9tiZm9HPvvFkeuDzWxz5LN/ysyuSnaumczMOprZf5rZC5HX+vwTxMwOmtl2M9tmZlsj17LqZ4+KtwCYWUfgYeAvgJHAdDMbmdysMt6TwJTLri0AfumcuwH4ZeS1tI964G+ccyOAzwH3Rf6b199B+zsHfMk591mgCJhiZp8DHgT+NfLZfwTcncQcs8EDwO4mr/X5J9afO+eKmrQIyaqfPSregjEOeM85t985dx5YA9yS5JwymnPudeDUZZdvAVZE/rwC+FpCk8oizrmjzrnfRf58mvA/Yv3Q30G7c2E1kZedIr8c8CXgmch1ffbtyMxCwFeARyOvDX3+yZZVP3tUvAWjH3CoyeuqyDVJrALn3FEIFxdAfpLzyQpmNgj4E2Az+jtIiMgju23AMeAV4PfAx865+kiIfga1rx8A/x24GHndC33+ieSAl83sLTObE7mWVT97cpKdQIawZq5pG69kPDPLA54F/ptzrjo8ASHtzTnXABSZWXdgHTCiubDEZpUdzGwqcMw595aZTWy83EyoPv/28wXn3BEzywdeMbM9yU4o0TTzFowqoH+T1yHgSJJyyWYfmtn1AJHfjyU5n4xmZp0IF24/c879W+Sy/g4SyDn3MfAa4XWH3c2s8X/I9TOo/XwBuNnMDhJeIvMlwjNx+vwTxDl3JPL7McL/8zKOLPvZo+ItGJXADZHdRlcB5cDzSc4pGz0P/GXkz38JPJfEXDJaZI3PY8Bu59z/avIl/R20MzPrE5lxw8yuASYRXnP4KvCNSJg++3binPt751zIOTeI8M/6Xznnvok+/4Qws1wz69r4Z+DLwA6y7GePmvQGxMzKCP/fV0fgcefcPyc5pYxmZquBiUBv4ENgIbAeeBoYAPwBuM05d/mmBgmAmX0R2ARs55N1P98mvO5NfwftyMzGEl6Q3ZHw/4A/7Zz7H2ZWSHgmqCfwn8Cdzrlzycs080Uem/6tc26qPv/EiHzO6yIvc4BVzrl/NrNeZNHPHhVvIiIiImlEj01FRERE0oiKNxEREZE0ouJNREREJI2oeBMRERFJIyreRERERNKIijcRicvMnJn9S5PXf2tmiwIa+0kz+0b8yDbf5zYz221mr152fVDk+/tWk2vLzGxGG+830cxubMsYl433m8jvg8zsjqDGjYz57ebuJSKpScWbiPg4B3zdzHonO5GmzKxjC8LvBv4f59yfN/O1Y8ADkSbbQZkIBFa8OecaxxoEtKh48/icPlW8NbmXiKQgFW8i4qMeWA7MvfwLl8+cmVlN5PeJZvYfZva0mb1rZkvM7JtmtsXMtpvZZ5oMM8nMNkXipkbe39HM/qeZVZrZO2Z2b5NxXzWzVYSbBF+ez/TI+DvM7MHIte8CXwR+bGb/s5nv7zjwSz7p0N50vM+Y2cbIIdibzGx4JLf9FtbdzC6a2YRI/CYzGwL8FTDXzLaZ2XgzG2hmv4x8L780swFNPr8fmdlvImM2OwvZ+LkCS4DxkXHntuRzMrP1ke9jp0UO9DazJcA1kfF+dtnfoUXG3hH5TKc1Gfs1M3vGzPaY2c/MwgfbRv6ed0VyWdrc9yIibaOD6UXE18PAO2b2/7bgPZ8lfGj6KWA/8KhzbpyZPQB8C/hvkbhBwH8BPgO8Gil+7gL+6JwrMbOrgTfM7OVI/DhgtHPuQNObmVlf4EHgz4CPgJfN7GuREwi+RLgb/tYouS4BXjSzxy+7vhz4K+fcPjMrBf63c+5LZvYuMBIYDLxFuKDaDIScc++Z2Y+BGufc0khuPwdWOudWmNks4EfA1yL3uJ5wcTmc8DE/z8T4TBdEvo/GIndOCz6nWc65UxY+VqvSzJ51zi0ws/udc0XN3OvrQBHhv8fekfe8HvnanwCjCJ/h+QbwBTPbBdwKDHfOOYsc4yUiwdLMm4h4cc5VAyuBv27B2yqdc0cjxwT9HmgsKrYTLtgaPe2cu+ic20e4yBtO+MzCu8xsG+Fjt3oBN0Tit1xeuEWUAK8554475+qBnwETPL+/A8AWmjySNLM8wo8+10by+AnhQgvCx4NNiPz6PuHiq4TwWcfN+TywKvLnn0biG62PfP+7gAKffJtoyef012b2NvBboH+TuGi+CKx2zjU45z4E/oPw99g4dpVz7iKwjfDfZzVwFnjUzL4O1LbwexERDyreRKQlfkB47Vhuk2v1RH6WRB6dNV031vRsx4tNXl/k0zP/l5/T5wADvuWcK4r8Guycayz+zkTJz3y/kSi+B/wdn/xs7AB83CSHIufciMjXNgHjCc9ubQC6E17n9jp+mn7PTT+nln4PXp+Thc/hnAR83jn3WcLnb3b2GDuapjk3ADmRgnkc8CzhWcWNLfpORMSLijcR8RY56PlpwgVco4OEH1MC3AJ0asXQt5lZh8g6uEJgL/AS8H+bWScAMxtqZrmxBiE88/RfzKy3hRfpTyc8W+TFObcH2AVMjbyuBg6Y2W2RHMzMPtvkXjcCF51zZwnPPt1LuKgDOA10bTL8b4DyyJ+/CfzaN6/LXD6u7+d0LfCRc67WzIYDn2vytQuN77/M68C0yLq6PoRnGbdESywyU3mtc24D4UfizT2KFZE2UvEmIi31L4TXPzX6/wgXTFuAUqLPisWyl3CR9SLh9WVngUcJF1K/M7MdhB9Zxlyn65w7Cvw98CrwNvA759xzLczln4FQk9ffBO6OPG7cSbhAJfIo+BDhR5AQLtq68skmip8DtzZuWCD8uHmmmb0DVAAPtDCvRu8A9Wb2tpnNxf9z2gjkRO7/j03yhvC6vncaNyw0sS5yv7eBXwH/3Tn3QYzcugIvRO7xHzSzwUVE2s6cu/xphYiIiIikKs28iYiIiKQRFW8iIiIiaUTFm4iIiEgaUfEmIiIikkZUvImIiIikERVvIiIiImlExZuIiIhIGlHxJiIiIpJG/g8HHMNSXW7ZgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5)) \n",
    "plt.xlabel('Number of Newton iterations')\n",
    "plt.ylabel('$f(v_t)-f^*$')\n",
    "\n",
    "cs = [0.001, 0.01, 0.1, 0.5, 1, 5, 10]\n",
    "w_seq_C = {}\n",
    "\n",
    "for c in cs:\n",
    "    initialize(c)\n",
    "    print('Barrier method for C =',C)\n",
    "    v_seq = barr_method(Q,p,A,b,v0,eps)\n",
    "    w_seq_C[C] = [get_w(v) for v in v_seq]\n",
    "    y_axis = [(f(Q,p,v) - f(Q,p,v_seq[-1])).flatten() for v in v_seq]\n",
    "    x_axis = np.arange(len(y_axis))\n",
    "    plt.step(x_axis, y_axis, label='$C$ = '+str(C))\n",
    "\n",
    "plt.semilogy()\n",
    "plt.legend(loc = 'upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We chose $\\mu=50$ which generally is a good trade-off between the number of inner iterations and the number of outer iterations, and generally gives a quicker convergence. We can see that in each case, the algorithm converge at a $1e-7$ solution in about $50$ iterations.\n",
    "\n",
    "Now, we want to check the impact of $C$ on $w$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w obtained for C= 0.001  :  [ 0.03515183 -0.04399957]\n",
      "w obtained for C= 0.01  :  [ 0.05119681 -0.09756987]\n",
      "w obtained for C= 0.1  :  [ 0.07466423 -0.14263758]\n",
      "w obtained for C= 0.5  :  [ 0.10345755 -0.46723647]\n",
      "w obtained for C= 1  :  [ 0.12944951 -0.76025446]\n",
      "w obtained for C= 5  :  [ 0.12944951 -0.76025447]\n",
      "w obtained for C= 10  :  [ 0.12944952 -0.76025447]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb8AAAGpCAYAAAAdnDj0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3iV5fnA8e979sggg4Qs9p5BtgoIiiLuqggqarWASq0/rKO17lon1bqqolVBZJa6QQUEQVQgQAh7j4SRhOycnH2e3x/BACZAIOOcJPfnurgg73nHfQKcO8+8NaUUQgghRFOiC3YAQgghRH2T5CeEEKLJkeQnhBCiyZHkJ4QQosmR5CeEEKLJMQQ7gNoQGxurWrduHewwhBBChJC1a9ceVUo1r+q1RpH8WrduTVpaWrDDEEIIEUI0Tdt/qtek21MIIUSTI8lPCCFEkyPJTwghRJPTKMb8hBAi1Hm9XrKysnC5XMEOpdGxWCwkJydjNBqrfY0kPyGEqAdZWVmEh4fTunVrNE0LdjiNhlKKvLw8srKyaNOmTbWvk25PIYSoBy6Xi5iYGEl8tUzTNGJiYs66RS3JTwgh6okkvrpxLt9XSX5CCCGaHEl+QgghmhxJfkIIIZocSX5CCNHEzJs3jwEDBtCzZ0/at2/P008/XaP7ffPNN3Tq1In27dvzwgsvnPV5pzp+5513EhcXR/fu3WsUX1Uk+QkhRBMybdo0XnzxRebPn09GRgbp6enYbLZzvp/f72fSpEksXLiQLVu2MGvWLLZs2VLt8053/R133ME333xzzrGdjiQ/IUKQx+Vhxt/nMa7dJG5tcy8fPjYLp0MWR4uaKS4u5oEHHmDu3LkkJycDEBYWxkMPPXTO91y9ejXt27enbdu2mEwmxowZw+eff17t8053/ZAhQ4iOjj7n2E5HFrkLEWKUUjxy6d/ZkbYbj8sLwH9f+ZLVC9fz5urn0ev1QY5Q1NTTX25my6HiWr1n18QInryq22nP+fTTTxkwYABt27Y94/0GDx5MSUlJpeNTpkzhkksuqfj64MGDpKSkVHydnJzMqlWrKl13qvOqe31tk+QnRIjZsGwzu9L3VSQ+AI/Ly8Gdh1mzMJ2BV/YJYnSiIdu8eTOpqanVOnfFihXVOk8pVelYVevuTnVeda+vbZL8hAgx29fsxntC4vuVs9TFttU7Jfk1AmdqodUVu92O0+ms1rnVbfklJyeTmZlZ8XVWVhaJiYmVrjvVedW9vrZJ8hMixDRPicFkMeIs9Z903GI3E5cSG6SoRGMwatQoxowZw+TJk4mPj8ftdjN9+nTGjx9f6dzqtvz69evHzp072bt3L0lJScyePZuZM2dW+7xOnTpV6/raJhNehAgxF17XH7PVVKnrx2A0MPSm84MUlWgM+vXrx1NPPcVll11Gjx49SE1NJScnp0b3NBgMvPnmm1x22WV06dKF0aNH063b8ZbtqFGjOHTo0CnPO931Y8eOZdCgQWzfvp3k5GT+85//1CjWE2lV9bc2NH379lVpaWnBDkOIWpO14xD/uPlf7N+cCWgkdUjg0U/+RJserYIdmjhHW7dupUuXLsEOo9Gq6vuradpapVTfqs6Xbk8hQlByx0TeTnuJguxCAgFFTEJUsEMSolGR5CdECIuKbxbsEIRolGTMTwghRJMjyU8IIUSTI8lPCCFEkyPJTwghRJMjya8BUEqRf6QAR3FZsEMRQohGQWZ7hrj0pZuYcte/yT9SiAoo+l7ai4c+mkREdHiwQxNCiAZLWn4hLGvHIR676gWy9+XidXnxeXykfZvO30Y9F+zQhBANWKgVs63LorWnIskvhH36xkJ8npM3OPZ5/ezdlMmejP1BikoI0ZCFWjFbqNuitaciyS+EHdxxCL8vUOm43qAje39uECISQjRkoVjMFuq2aO2pyJhfCOs5tCsbV2w9qa4bgM/jo11q6+AEJYSouYV/gSMba/eeLXrA5afucoTQLGYbLJL8QtiVEy/l09cX4veVVLQAzTYzw8deIKVthBBnLRSL2QaLJL8QFhETzttrX2TaE3NYtWA9tggr1953OVffe1mwQxNC1MQZWmh1JRSL2QaNUqrB/+rTp48SQohQtmXLlmCHoFavXq3atm2rjhw5opRSyuVyqalTp9bonl6vV7Vp00bt2bNHud1u1bNnT7Vp06azPm/v3r2qW7du5xxHVd9fIE2dIm/IhBchhGgiQrGYLdRt0dpTkWK2QghRD6SYbd0622K20vITQgjR5EjyE0II0eRI8hNCCNHkSPITQgjR5AQt+WmalqJp2lJN07ZqmrZZ07T7jx2P1jRtkaZpO4/9HhWsGIUQQjROwWz5+YA/K6W6AAOBSZqmdQX+AixRSnUAlhz7WgghhKg1QUt+SqnDSql1x/5cAmwFkoBrgGnHTpsGXBucCIUQQjRWITHmp2laa6A3sAqIV0odhvIECcSd4poJmqalaZqWlpsrFQ6EEEJUX9CTn6ZpYcB84P+UUsXVvU4pNVUp1Vcp1bd58+Z1F6AQQjQywSpmG4yitacS1OSnaZqR8sT3iVLqf8cOZ2ualnDs9QSgZnvvCCGEqBCsYrYQnKK1pxLM2Z4a8B9gq1LqlRNe+gK4/difbwcqV0UUQghx1oJZzBaCU7T2VIJZ0ugCYBywUdO09GPHHgVeAOZqmnYXcAC4MUjxCSFEnXhx9Ytsy99Wq/fsHN2ZR/o/ctpzglnMNtQELfkppX4ETlXJ8OL6jEWIpszn9bFxxVZ8Hh89hnTFYjMHOyRRR4JZzDbUSDFb0WT5/X5+/iKNVV+vJTI2gpF3XUxyh4Rgh1WvNq3cxhPXvIjf5wcg4A/w0IeTGHLDoCBH1ridqYVWV4JZzDbUSPITTZLP6+OREX9nx7o9uEpd6I16PntjIQ99OImho88Pdnj1wlnq5NFR/8BZ4jrp+Eu3v0mHPm1JaBMfpMhEXRk1ahRjxoxh8uTJxMfH43a7mT59OuPHj690bnVbfv369WPnzp3s3buXpKQkZs+ezcyZM2s79FoX9KUOQgTD4hkr2J62G1dp+Qe/3+vH7fTwzz+8jdvpDnJ09ePnL9KginKefr+fxR8vr/+ARJ0LZjFbCE7R2lORlp9okr6fuQJ3WeUkp9NpbPl5B72H9whCVPXLUezE7w9UOu7z+CkpKA1CRKI+jBs3jnHjxtXqPUeNGsWoUaOqfG3BggUVf541a1atPrcmpOUnmiSLvepJHQGlMFtN9RxNcPS+uAdUMVnBYjczYNR5QYhIiPojyU80SVdMGFFlArSGWeg8oEMQIqp/yR0SuOruS0/6PljsZlKHdS9PjEI0YtLtKZqk/pf35qp7LuPzNxei0+vQ6XTojXr+8dWj6HRN52fCif+8nT6XpfLNf5bgcXu5+ObBDL5hYJP6HoimSZKfaJI0TWPCS+O4+t7LSF+6mYjoMPqOTMVkNgY7tHqlaRr9Lkul32XVW/slRGMhyU80aS1axzHy91UWDhFCNGLStyGEEKLJkeQnhBCiyZHkJ4QQosmR5CdEI+Dz+lg8YzlPXPsiL93xJpt/2h7skEQIC1Yx29atW1fsLNO3b98aPbOmZMKLEA2cz+vjoYufZtf6vbgcbjRNY/l/f+G2p0Yz+sGrgx2eCDHTpk3jjTfe4LPPPiM5OZnS0lLefvvtc77fr8VsFy1aRHJyMv369ePqq6+ma9euVZ6/dOlSYmNjz/l5tUVafkI0cMvn/VyR+KC8xIy7zM20J2ZTdLQ4yNGJUBLsYrahRFp+QjRwP366uiLxnchgMpDxwxYGXz8wCFGJ0zny3HO4t9ZuMVtzl860ePTR054T7GK2mqZx6aWXomkaEydOZMKECWeMo65I8hOigQuPDkPTaahA5X06bRHWIEQkQlWwi9muXLmSxMREcnJyGDFiBJ07d2bIkCHVek5tk+QnRAM3avwlLPmkcpUKk9lI6rDuQYpKnM6ZWmh1JdjFbH89HhcXx3XXXcfq1auDlvxkzE+IBq5T33ZMeOlWTBYTtggrtnArzeIief7bx9Ab9MEOT4SQUaNGMW/ePLKzswFwu9289957VZ67YsUK0tPTK/06MfHBycVsPR4Ps2fP5uqrK0+0cjgcFcnU4XDw3Xff0b178H44k5afEI3A1feOZPjNg8lYvgVbuJUeg7tI4hOVnFjM1u/34/P5uPXWW2t0zxOL2fr9fu68885KxWzff/99XC4X1113HQA+n4+bb76ZkSNH1ujZNaFV1V/b0PTt21elpaUFOwwhhDilrVu30qVLl2CH0WhV9f3VNG2tUqrKBYXS7SmEEKLJkeQnhBCiyZHkJ4QQosmR5CeEECJk+PwBXF5/nT9Hkp8QQoiQ4PT42ZVTyv68sioXz9cmWeoghBAi6ArLPGQVONHrNFpGWU+5S0xtkeQnhBAiaJRSHCl2kVvixmYy0CrGhlFf952SkvyEEEIEhc8fILPASYnLS7TdRGIzK7o6bvH9Ssb8hBCiiantYrZ33nkncXFxZ7VdmcvrZ3duKaUuH0nNrCRH2eot8YEkPyGEaFKmTZvGiy++yPz588nIyCA9PR2bzVaje95xxx1888031T6/yOllV04p/gC0bW4nJsxco+efC+n2FEKIJuLXYrZr1qyptWK2AEOGDGHfvn1nPE8pRU6Jm+xiF1aTnlbRdkyG4LTBJPk1UYf3ZvP9zB9xOdwMvLIPXQd1rPPZVUKIcivm7uBoZmmt3jM2JYzBozue9py6KGZbXf5AgMx8J8UuL1E2E0nNrOh0wfvMkeTXBC3+ZDmvTniXgM+Pz+fnszcWMPj6gTz04SRJgEI0YnVRzLY63F4/+/LK8PgCJDazEmM3Bf2zRpJfE1Na6ODVCe/icXoqjrkcblbM/4VhYy+k32XV+48hhDh3Z2qh1ZW6KGZ7JsVOL5kFZWhAm1gbYRZjta+tS5L8mph1izMwGPR4fnPc5XCzbPZKSX5CNGKjRo1izJgxTJ48mfj4eNxuN9OnT2f8+PGVzq1py08pRW6JmyPFLqxGPa1ibJhCqMakzPZsYnR6HVTR26BpoDeGzj9MIUTtO7GYbY8ePUhNTSUnJ6fG9x07diyDBg1i+/btJCcn895773Mgv4wjxS6aWU20ax4WUokPpOXX5PS5tBcBf6DScZPVzIhxQ4MQkRCiPo0bN45x48bV6j1nzZpV8We3z8/+vDKKnV4SIi3EhpmDPr5XlSbf8lNK8ctXa3nmxik8fcPL/PT5GgKBysmhsbDaLTw+98+YbWYsdjMmixGTxcR1911Oj8FSZVoIce5KXOXr97z+AK1j7TQPt4Rk4gNp+fHK+HdYNmclLocbgLRvN3D+tf35y/T7QvYvrab6X96bWZnvsPLT1bgcbvpdnkpS+4RghyWEaKCUUhwt9XCkyIn52PieOcS6OX+rSSe/Xel7WTr7R9xlJ898XPnparav2UXn/h2CGF3dCo8KY+Sdw4MdhhCigQsEFFmFTgrLPERajSRH2dAHcf1edTXpbs+132Xgq6JoosflIe3bDUGISAghGg6Pr3x/zsIyDy0iLLSMbhiJD5p4y88eacNgNOD/TQI0mgzYIqxBikoIIUJfqcvHgfzyorOtY+xEWENj/V51NemW35AbBlZ5XNM0ho25oJ6jEUKI0Fc+vudm71EHep1Gu7iwBpf4oIknv4iYcJ7630PYIqwVv6zhFh6f+wBR8c2CHZ4QQoSUQECRVeDkUKGTcIuB9nF2LA10fXCT7vYE6HtpL+Zl/4eMH7aglKLX0K6YLKZghyWEECHF6wuwP7+MMo+PuAgL8eGhuX6vupp88gMwmY30vbRXsMMQQoh6MW/ePKZMmYLT6aSsrIxx48bx5JNPnvJ8h9vH/rwyAkrRKsZO5G+6OVu3bk14eDh6vR6DwUBaWlpdv4Uak+QnhBBNyLRp03jjjTf47LPPSE5OprS0lLfffvuU5+eVujlU5MKk12gbE3bKbs6lS5cSGxtbV2HXOkl+QgjRRJxNMduAUhwqdJLv8BBuMZISZcWgbzzTRCT5iUYn/0gBzlIXCW3j0ekaz39W0Xgs/WgqOfv31Oo941q1ZdgdE057TnWL2Xr9Ac6/4EJKSkrQ6zSMJyS9qkoaaZrGpZdeiqZpTJw4kQkTTh9HKJDkJxqNvMMFPHvTK2xfsxu9Xoc1wspDH06SMk21JBAIsPLT1Xw/cwUGs5GRvx/GeZf0bNCTHpqa6hSzLfOUj+99OH8hyVFWmtnOPAFw5cqVJCYmkpOTw4gRI+jcuTNDhgyprbDrhCQ/0SgopXjk0mfI2n4Ivy+AF3CVuXn6+im8s+4lkjsmVnndod1HWPDeYrL359JnRC+Gjb0As9Vcv8E3AEopnrlhCmsXZVTsg/vLl2lcMXEEd0+5PcjRNTxnaqHVlTMVs813eDhY6MSo05h40xU4SksrnVNVyy8xsfz/V1xcHNdddx2rV68O+eQnfUKiUdi+ZhfZ+4/i951ckcPn8fL5W99Uec2ab9OZ0OtB5r/6Fcvm/MRb93/APX0ewVFcVh8hNyjrv9/E2kUbKxIflO+D++W/v+XgrsNBjEycjVGjRjFv3jyys7MBcLvdvPfeexXje1kFZdhNetrHhbHyxx9JT0+v9Ou3ic/hcFRUfHc4HHz33Xd079693t/b2ZKWn2gUjh7MR1fFnoJ+X4DDe7IrH/f7eem2N3CXnfxhnr0vh/mvfsVtT44+7fN2rN3Nqq/XYbGZGTp6EHEtm9f8TYSw1QvW4XK4Kh3XNI2132VIVZAG4sRitn6/H5/Px80338Leow4cbh+xYWYSIs+uDFF2djbXXXcdwLH73czIkSPr6i3UGkl+olHo2LcdPo+v0nGzzcR5l/SsdDxz2yFcTnel4x6Xl+Xzfj5l8lNK8fqk91k0/Qc8Lg96g56PnpzDn9+/m+FjB9f8jYSosGZ2DCY9Ps/J++Dq9DrskbYgRSXOxYnFbJ3HxvfKPH5SomxE2c9+g4+2bduyYUPDKwQg3Z6iUYhLieXS2y/CYjs+XmcwGYiMjeCy3w+rdL7ZZiLgV1XeyxJmOeVzNizbzOKPf8Bd5kYFFD6PD4/Twz//8A6lhY6av5EQdfGtg9Hrq17fNejqvvUcjagNhWUeduc6UEC75vZzSnwNmSQ/0Wj86d/jmfTGXbRLbU2LNnFc88eR/DvtRewRlVsmCW3iSemUWKmr1GI3c829p+6y+X7Wjyd1lf5Kb9CR9m16zd9EiEpoE89DH07CbDNX7INrj7Tx7Fd/xRYuFVAaEqUUh4ucHMgvw2osH9+zmZpeJ2DTe8ei0dI0jZG/H8bIKlp6VXly/oM8NPxpivNLQIHf52f4zRdyybhTz1IrXzeoAb9tNWpojXxN4dDR59P/ivPIWLYZg8lAz6FdMZoa3m7+TZnPH+BAfhmlbh8xx8b3dE10qYokP9FkJbSJZ/ruN9mwbDP5hwvpen5HEtrEn/aaS24dzOIZyyu1/gJ+P/1GNv71hFa7hQFX9Al2GOIcOL1+9uc58PoVyVFWou1Ne0mPJD/RpOl0OnoP71Ht87tf2IVrJl3GZ29+g/IH0Bl0KAV/nXG/dP+JkFVU5iGzwIlep9E21o7dLB/98h0Q4iyNf3Ecl/1+ePlSB7uZwdcPoFnzyGCHJUQlSimyi13klLixmQy0irGdtFVZUybJT4Q8v99PwZFC7M3sWO2nnolZn1p2TqJl56RghyHEKfkCAbLynRS7vETbTCRGWZvs+F5VJPmJkLbo4x9454FpuMrcoBTDbx7MfW/9AZNZJloIcSour5/9eWV4fAGSmlmJtptkD9bfkPavCFnrFmfw2j1TKc4rweP04HF5+X7Wj/xr4rvBDq1Gdq7bw5yXPmfBe4sb9dpAERzFTi+7c0rxBxRtm9uJCatccX3evHkMGDCAnj170r59e55++ukaPfPOO+8kLi6uym3NvvnmGzp16kT79u154YUXavSc2iTJT4SsT56dj7vMc9Ixj9PDsjk/NcikEQgEeOG2N5g85Ak+fHwWbz/wEWNb3s2GHzYHOzTRCPw6vrcvz4HJoKN9XFiVE1umTZvGiy++yPz588nIyCA9PR2brWa79Nxxxx18803lPXT9fj+TJk1i4cKFbNmyhVmzZrFly5YaPau2SPITIevIvpwqjxuMegqyC+s5mpr78X+rWPnpKtxlbvxePy6HG1epi6evn4LPW3lrNiGqyx9Q7M8rI7vYRZTNRLvmYZgMlT/efy1mO3fu3DMWsz0bQ4YMITo6utLx1atX0759e9q2bYvJZGLMmDF8/vnnNXpWbZExPxGyug7qyNGsPAKBytuQxbeOC0JENfPth0tPqorwK7/Xz5afd9BzSNcgRCWCofDL3XgO1U7vhVIKl9ePL9ZK4pVtiQk79fhedYvZAgwePLiiWsOJqippdCoHDx4kJSWl4uvk5GRWrVpVrWvrmiQ/EbJue2o0q75eh8vhRqnyBGixmbn9mZsa5ISXX99DJdppXhPiNPyBAO5jZbwiLAZiw0+/cL06xWx/tWLFihrHV9W/62pNvFEB0Oq2Y1KSnwhZKZ2SeOOX5/jw8dlsXrmdmMQobn70dwy5YVCwQzsnI267iI0rtlZq/WmaRtdBHYMUlQiGZle1q9H1SilyS93kFrmwGPW0jrFhMlS98fiJzlTM9kS10fJLTk4mMzOz4uusrKyKwrdV8XvKyC/JJADEx3Sq1jPOlSQ/EdJadU3hqfk1G48IFUNHD2L5f38m7dt0XGVuTBYTmgaPzXlA9sgU1eYPKLIKyihyemlmNZIUZUNfRS3LqowaNYoxY8YwefJk4uPjcbvdTJ8+nfHjx1c6tzZafv369WPnzp3s3buXpKQkZs+ezcyZMyudF/C6KCjJIjfgwo9GhM6EUqpOl2cENflpmvYBcCWQo5TqfuxYNDAHaA3sA0YrpQqCFaMQtUWn0/HEvD+z+aftrF20gYjocC4acwFRcbI7jKgej8/PvrwyXF4/LSItNK9iGcPpVFXM9tZbb61xXGPHjmXZsmUcPXqU5ORknn76ae666y4MBgNvvvlmxfPuvPNOunXrVnGd8nspLM4k11eGV9Owa0biwhOxmcJrHNOZaMEca9A0bQhQCkw/Ifm9BOQrpV7QNO0vQJRS6pHT3adv374qLS2t7gMWQohztHXrVrp06XLO15e6vOzPLwOgZbSNcEvD7S1QAR8lxQfJ8Zbg1jQsmo54ewJ2c+Q5t/aq+v5qmrZWKVVlwcmgtvyUUss1TWv9m8PXABcd+/M0YBlw2uQnhBCNlVKKo6UejhQ5MRv1tIq2YTaeeXwvJAUClJYeIsddiFPTMGk6UmzxhFuj630HmlAc84tXSh0GUEod1jStyjntmqZNACYAtGzZsh7DE0KI+hEIKA4WOiko8xBpNZJ8FuN7IUUpykqPkOPKw6FpGDUdidYYmtnigrbtWigmv2pRSk0FpkJ5t2eQwxFCiDM6m0kcHl+A/XkOnF4/8REW4sLPbnwvJCiFqyyXnLJcSjTQaxotzFFEhbVAV4tLGc5l+C4Uk1+2pmkJx1p9CUDV23wIIUQDYrFYyMvLIyYm5oxJrNTt40BeGUopWsfYibA2sPE9pfA4C8h1HKFQU+g0aG6KICYsEb2udrtslVLk5eVhsZxdxZdQTH5fALcDLxz7PTT2whGinvh9flZ9vY6NP26leUoMF988mMjYiGCHJWooOTmZrKwscnNzT3teqdtHUZkXg14j2m7iYLGOg/UUY23we52Uuosoo3zxvV1vJszcjDxdKXnsqJNnWiyWiu3aqivYsz1nUT65JRbIBp4EPgPmAi2BA8CNSqn8091HZnuKxsJV5uaBoU+Qtf0QzlIXJqsJvUHHS4ueoHP/DsEOT9Qht8/P459tYm5aFsM7x/GvMalENKAZncVZa/ho6SPM8GXj0TSujerO3cNepkVEypkvriOhPNtz7CleurheAxGNXiAQYOH7S/jsjYU4issYdHU/bn38hpBbYzf/1S/ZvzkTj8sLlFexAPjHmFeZvvuthjfmI6olu9jFxI/Xkp5ZyH3D2zP5ko7oGsjEFmfONmYu+TMfOPdSrNczMrwtky56idax576soz6EYrenOIMda3ez4L0llBY6GHz9QC68rj/6amxt1JS9dvdUvp/5Y3lRXGDB1EX89Nlq3t/0CvZIe5CjO27JjBUVie9EBTlFHNp9hKT2CbXyHJ/Xx7Sn5vLlv7+lrMRJp37t+eMbd9Gpb8223RJnb+3+Au6esRaH28fbt5zH5T1q5++4rnmLDvLpogd4p2gjuQY9F9oS+dPgZ+mSNDDYoVWLJL8G5vO3FvLeIzPwurwEAopVX6/jq3e+44VvH5MEeAo5B3JZNGM53hOSis/rpyS/lAXvL+HGP18dxOhOdqq/QxVQGIy199/1n3e9zYr5v+A+1rLctmonDw57krfXvUxyh4bx4dsYzF59gMc/30RCpJUZdw2gU4u639mkpgJl+Sxc/BBv5f5MpkFPb0s0Lw96gj5tLw12aGdF6vk1IMX5JUx96GPcZZ6KMj8uh4ttq3eyYv4vQY4udO1ctxejqXLicDs9bFgaWoVkR42/GLPNdNIxTYOEdi2Ib9W8Vp6Rd7iAH/77c0Xi+5XH5WXuyzK/rD54fAEe+2wjf/nfRga2jeGLP14Q8olPuR388O1kbpx5IX8pWI3VFM5b/R5j2s0rGlziA2n5NSgZP2zBYDJU6hZzOdwsn/czF910QZAiC22xyTEE/IFKx/UGPQnt4oMQ0alddc9lpH23gYwftuD3+TGYDJgsJp6Y9+dae8ahXUcwmY0ntYQBAv4Au9bvrbXniKrllri595O1rNlXwMShbXn4ss6hvXDd7yVtxXO8tmsu6UYdKUYLL/acyMjU8bW6Vq++SfJrQKxhFqhicq6madgjbfUfUB3LOZDL+u83YY+00f/y3pgspjNfVIWOfdqS0DaeA1sP4vf5K44bTAaumTSytsKtFQajgWe//Cvb1+xiy887iEmMZtDVfWu1fmFi+xZVjivq9Dra925Ta88RlW3ILGTix2spdHp4fWxvru516vI+QRcIsHXNW7y2cSorjRBnMPB4x7FcN+DPGHUNZxbqqUjya0B6Xd8ZPeYAACAASURBVNQNg9kAvymxZbKauHx89eprNRQf/G0m81/9Cp1eh06nQ6fX8cK3j9GpX/uzvpemabz43eM8d8trbF65DZ1OR3h0GA9+cC/JHUPvw0fTNDr371BnSxtiEqIYcsNAfvzfqpO6Pk0WE6MfDJ3xz8Zm/tos/vrpRpqHmZl/z/l0SwytmcYVlGJfxie8mfYK3xq8RBjggZZXMvbCJ7AYrcGOrtYEdZ1fbWlK6/y2p+3mryOfLW/BKPB6fPz+2TEhNWmjptYtzuDJ616qVPQ1snkEcw5OrdHEnqKjxThLXcS3at6klw34vD6mPTmXL98+Ybbn63ee0w8X4vS8/gDPLdjKhyv3MahtDG/dch7R9nPrxahrR7Yv4J2fnuYzzYEJjXEJg7lj6POEW0I0UZ/B6db5SfJrgHxeH+uXbKSs2EmvYd1o1rxh/sM8lb+P/ifL/1t5Ao8t3MozXzxCr6HdqrhKiNCT7/Aw6ZN1/Lwnj99f0JpHR3XBqA+9cbKCAyt5f9lfme3PR2kwOqY3fxj2MrFhLYIdWo2E7CJ3cW4MRgP9RvYOdhh1xlXmqfoF7fiibyFC3eZDRUyYvpbcUjdTbuzFDX3Obvut+uDI3sT0JQ8yzZ2JU9O4KrIj9w57mcSoxr/eU5KfCDnDx15Ixg+bK3V7+n0Bug8O7V0jhAD4YsMhHv7vBppZTcybOIheKc2CHdJJ3AX7mbN4Mu+XbKNAr+cSewp/HPI87Vo03h+qf0uSnwg5F910Pt9OW8rWX3biKnWhN+jRG/VMnjoRq/3sdm4Xoj75A4qXvt3Guz/soV/rKP59Sx+ah5uDHVYFnyOXLxY9yNt5azhi0DPQEsf9FzxN91ZDgx1avZPkJ4LKUVzGew/PYOmsH/H7Awy8sg/3vHoHzy/8G6u+XsfPX6YRER3OyDuHkdIpKdjhCnFKhWUe7pu1nhU7j3LrwJY8cWU3TIbQGN8LuEtY9P1fefPg9+wz6ulhiuTv/R9hYKdrgx1a0MiEFxE0Sinu7fcI+zZl4vP4ANAbdETFN+PD7a9jsYXOT8xCnM72IyVM+DiNQ4VOnrmmO2P7twx2SAAor5ufVvyd1/Z8ylajjnaYua/3JIb3uKNJzHaWCS8iJG1YtpmDOw5XJD4oH9crLXTww9yfuOyOYUGMTojq+WbTYR6YuwG72cDsCQPp0yo62CFBwM+GX/7Fa1s+Yo0REg1Gnu18G1f2u7/Wi8k2VJL8RNDs25R50o4rv3I53Oxct1eSnwhpgYDi1cU7eOP7XfRKaca7t/ahRWSQx6SVYkf6h7yx7nWWGfxEGzT+0uZabjz/b5gM0pNyIkl+ImiSOiagN+rhN1ttWexmWncLXgFMIc6k2OVl8ux0lmzL4cY+yfz92u5YjMFtUWVu/ZS3fnmOBZoTu17jvsTh3DrkH9jMYUGNK1RJ8hNB02dET2ISozm8Jxu/t7wFqNNpmG1mht98YZCjE6Jqu3JKmfBxGgfyynjmmm6MG9gqqONnufuW8e7yx5gfKESvwR3N+3PXsJeItMUGLaaGQJKfCBqdTse/Vvyd1+59j58+W4NSil4XdWPyuxOxhTeePQRF47F4Szb/Nycdk0HHjD8MYGDbmKDFUnQ4nQ+/f4hPPIfwaRq/i+zKxOFTiIsMjck2oU6SnwiqyNgInpj7ZwKBAEop9HoZjBehJxBQvLl0F68s2kG3xAim3taXpGbB+QGtLG83Mxc/wAeOXZTqNC4Pa8Oki16kZXPZ9u9sSPITIUGnC431UEL8Vqnbx4NzN/DN5iNcm5rI87/ridVU/z+keUuz+e+iybxbkE6eXs9QawvuG/IsnZIG1XssjYEkPyGEOIV9Rx1M+DiNXTmlPHZFF+66sE29j+/5nUUs+P4R3jqynIMGPX1MUbw68DF6t7+8XuNobCT5CSFEFX7Ykct9M9eh02lMv3MAF3ao3wkkyuti6bIneOPA1+wy6OhitPNYnwe4oOuYJrFAva5J8hNCiBMopXh3+R5e+mYbHePDmTquLy1jbPUXQMDP6p9e4rXtM8kwQGu9iZe7j+fS8+5Gp8nwQG2R5CeEEMeUeXw8/N8Mvso4zBU9Enj5xp7YTPX0MakUm9Pe5bWMt/nZECBer+Op9qO5ZuDDGPTG+omhCZHkJ4QQQGZ+GRM+Xsu2I8U8PLIT9wxtV2/di3s2zeHNNS+zSOemmQ4eTLmcMUOewWyQKiZ1RZKfEEHm8/pY+elqNq3cRnzr5owYN5TI2Ihgh9Wk/LTrKJNmrsMXUHxwRz+GdYqrl+ce3r2If//4JF+oYiwa3BN3Abdd9CJh1tCq/9cYSfITIojKSpzcf8HfyN6Xi7PUhclq4uOn5vHSkifp1LfxV9MONqUUH6zcx3MLttIm1s57t/WlTay9zp+bl7WG95c9zBxfLgC3RPXiD8OnEB2eUOfPFuUk+QkRRLNf+JSDO4/gdZfvb+pxevAAz9/yLz7c9rrM6qtDLq+fRz/dyP/WHWRE13heGd2LcEvdjq2V5G5j2pI/87FzHy5N45qwdtwz7GUSYjrW6XNFZZL8hAiipbNXViS+E+Vm5pGblUdciuzPWBcOFTq5e8ZaMrKK+L9LOvCn4R3Q6eruBw1X8UHmLHqA94o2UaTXMcKaxB+HvkDbhPPq7Jni9CT5CRFEhlNUAlDq1K+Jmlm9N597P1mL0+Nn6rg+XNqtRZ09y+ss4LPFD/JO7i/k6HWcb47hTxc8RbfWw+vsmaJ6JPkJcRplJU52rd9LZPMIWnVJrvX7X/6HS5j+5BzcTk/FMU2n0aZ7CtEtomr9eU2ZUooZqw7w9BebSYm2MWv8QDrEh9fJswKeMr5d+jfeyvqO/QYdPQ1hvND/Yfp1vr5OnifOniQ/IU7hv698yYePz8Zg1OP3+WnZOYlnv/prrSal6/50OenfbyRj+VZUIIDBaMASZuGxOQ/U2jMEuH1+nvx8M7PXZHJRp+a8NqY3kdbaH99TPi8/rnyO13fNY5tBo73BzOu97uWiXnfJ+G2I0ZRSwY6hxvr27avS0tKCHYZoRNK+28BTv3sZd5m74pjeoKN97za8ueqFWn/e9jW72LpqJ82TYxhwxXkYjPJzaW3JLnZx94y1rD9QyKRh7XhgRCf0tT2+pxTrV7/Bvza9zzqDIknpmNT5Vkb1fwC9Trqvg0XTtLVKqb5VvSb/w4SowvxXvzwp8QH4fQH2bcrk4K7DJLWv3Snpnfq1p1O/9rV6TwHrDhRw98drKXH5eOvm87iiZ+0vJdi+YQavr32F5XovsTr4W6uruf7CJzEaTLX+LFF7JPkJUYWC7KIqj+uNeorzSkmSPBXy5qw5wOOfbSY+0sy0O8+nS0LtbhxwYMcC3vz5GRbiIFxT3N9iGDdf9Dw2c1itPkfUDUl+QlRh4FV9ObA1C6/bd9LxgD9A255SKTuUeXwB/v7VFj7+ZT8Xto/ljbG9ibLXXiss+8BK3v3hUT7152EE/hDThzuGTyHS3rzWniHqniQ/Iarwu/tH8d1HSynKLcbj8qJpYLKaufuV2zFbzcEOL6Rk7TzM+iUbCY+yM/Cqvlhswfv+HC11c++Mdazel8/4wW14ZGRnDPraqYRQmL2RD75/iJmuLPwa3BDRmYnDpxDbrHWt3F8cF3A60VmtdfoMSX5CVCEiOpx306fw+ZsLWbVgPTGJUVz/f1fSY3CXYIcWMpRSvD35I76eughN09AZdGiaxvPfPEbXgfW/Y8nGrCImfJxGvsPDv25K5dreSbVy37LCA3y8eDIflWzDoWlcaU/hnoteJiWue63cv6lTHg+urVtxpqdTlp6OM30DpqQkWs34uE6fK7M9hRDnZNWCdTx70yu4HCdPDIqIDWfuoffQG+pvluP/1mXx1/9tJMZuYuptfemeFFnje3ocucxb/CBT89LI1+sYZozlvsHP0iHlglqIuOnyZufgTE+v+OXavBnlKV/nakhMwJaaiq3/AKLG3FTjZ8lsTyFErVv4nyWVEh+Az+1j80/b6Tmka53H4PMHeH7hNv7z414GtInmrVvOIzasZt2ufncpX37/F94+vJRDeh39jBG8PvAxenW4opaibjqU14tr2zac648nO++hQwBoRiOWbt2IuvlmrKmpWHunYoyPr7fYJPkJIc6Jx1V5T1IANKrcr7S25Ts8/HHmOn7anccd57fmb1d0wViD8T3l87Jk+VO8sfdz9hg0uuqtPHne/QzqfqssUK8mX27usa7L8u5L16ZNKHf5D0iGFi2wpqYSdds4bKmpmLt2RWcK3nIQSX5CiHNy8c2D2bh8S6XWnwooul/YuU6fveVQMRM+TiOn2M1LN/RkdN+Uc79ZIMAvv7zCa9ums0mvaKM38ErXO7mk732S9E5Deb24tu84qQvTm5VV/qLRiKVrF6LG3FTeqktNxZgQWuWaJPkJIc7JRTedz5JPlrPxx224Sl0YTHr0ej0PfTipTmfEfpVxiIfmZRBhNTBn4kB6tzzH7eaUYuP6D3gt/U1W6X200DSeaXMDV13wKAZ93ZY2aoh8+fnlSW79+vJuzE2bUC4XAIa4uPJW3bEuTEu3rujMoT0rWpKfEOKc6A16nv3qr6xdlMGqr9cSERPOiNuGktCmbsZt/AHFlO+28/ay3fRpFcXbt55HXLjlnO61a+unvLHqeb7XnERp8HDiCEZf9Cxmo62Wo26YlM+He8eOk7owvQcOlL9oMGDp0oVmN96INbUXttRUDImJDa6VLLM9hRAhr6jMy59mr+eHHbmM7d+Sp6/uhslw9uN7B/cu498/Ps5X/gKsCm6PG8htw1/Gbm3aFTR8BQUVSc65fn15q66sDAB981hsx7ouy1t13dBZzu2Hjury+/3o9TWfLSyzPYUQDdaO7BImTE/jYKGTf1zXnVsGtDrrexw9vI73lj7CXM9hdMC4yG7cdckrRIXXzlrAhkT5/bh37jzWhVnesvPs31/+ol6PpXNnml133fEZmElJddaq83q8HF67l6MbD+A+WIjB4cemrLjxcN4/b6iTZ/5Kkp8QImR9u/kID8xJx2oyMHP8QPq1jj6r64vzdvPRkgeY4diFR9O41t6Gu4e9TIvYup2QE0r8hYU4N2yo6MJ0bcgg8GurLjoaa+/eRN5wPbbUVCzdu9fJzip+v5/sjZnkpO/DlZmPVuzFFjARro/EoDMShxmIp5RiHDjwRNT9GlFJfkKIkBMIKP61ZCevL9lJr+RI3hnXh4TI6n8oO0uOMGvJA/ynYAPFOh0jLQlMGvI8rZP61WHUwaf8fty7dp80A9Ozd2/5i3o95k4dibz2muMzMFNSar1Vl7frMAdX76Zs31G0IjcWr5FwQwQmnZlY9EBznMpBKaUc0R9Fi7ES2bEFSQM7kBzXrFZjOR1JfkKIkFLi8jJ5zgYWb83m+vOS+cd13bEYq9cS8LqK+XTJQ7yT/SO5eh0XGqP40/lP0qXtiDqOOjj8RUU4MzKOLyLPyCBQWgqAvlmz8lbdtdeWJ7se3dHZam9CT8nhAjJ/2UHJrmxUnhOzV0+4LhyL3kY0EE0M7oCLEq2YHC0PmpkJaxNHUv92JNfRpKizIclPCBEy9uSWMn56Gvvyynjqqq7cfn7rarVMAj43C5c9xlsHFpKp1+itt/Nyvwfp03V0PURdP1QggGfPHsrWr6+YnOLZvbv8RZ0Oc8eORFx5BdbUVGypqRhbtaqVVp2zoJTMn3dSuOMQ/hwHJrdGmGbHbggnAoggCm/ATgnF5GmFqDAH1laxtOjTloROCbUycaUuSPITQoSE77dlc/+sdIwGHTPuGsCgdjFnvEb5/Sz/+UVe3zGLHXroqDPyVo+7Gdx7QoObev9b/pISnBsyjndhZmQQKC4GQB8ZiSW1F5FXXVk+A7N7D/Rh9ho9z13m4uCq3eRvycJ7pBhjmcKODbshHJumw0YEfmWjhGKKKCHfWoYlKYrmvVrRsldr9KaGlU4aVrRCiEZHKcW/l+1mynfb6ZoQwbvj+pAcdYbuOaVIW/sOr2e8y3q9nxRN48UONzNy4EPodKHZ0jgdFQjg2bevfJnBsWTn3rUblAJNw9yhAxEjR1aM1ZnaVK9FXBW/x8eh9XvI2XAAz6FC9KUB7MpMmCESk6anBTYCylI++UQro9DkxJQQSXSPZJL7tsdobRwV6iX5CSGCxuH28dB/N7Bg4xGu7pXIi9f3xGo6ffLaumkOr6W9zErNTXMNHk+5guuGPI3RENo7ipzIX1qKKyPj+CLyDRkEiooA0EVEYO3Vi/CRI7H17o2lZ0/0YWdfHd7v95O75SDZ6/ZQdiAf3W9mWMYfm2HpoIRSHBQbcjA0D6NZ5wSSB3WkZWTNWpKhTpKfECIo9uc5mDB9LTtzSvjbqC78YXCb07Zm9u3+jjdXPsW3qoQIpXggfjBjhr+I1RxRj1GfPaVUeasufUPF9mDunTuPt+ratyPi0hEntOraoOnObgF/3p4jHFq1G8f+XChwY/EZiNBFYNJbiEFPDM1x8usMy1y0GBsRHeJJHtCR5Bb1N8MylEjyE0LUuxU7c/njzPUAfPT7/gzp2PyU5x7JWsU7P/yVz7w5mBRMiEnljuGvEB4WV1/hnpWAw4Fz48bji8g3bMBfWAiALjwca8+ehI8YgbV3b6w9e6CPqH7yLjlSSNYvOyjelY3KK8Ps0ROmC8OqtxMFRBGDJ+CimGJydPnlMyxbNyehfzuS27aoo3fcMEnyE0LUG6UU763YwwsLt9EhLpypt/WhVUzV3WsFudt4//s/M9u5nwBwU1gHxl88hdiodvUb9GkopfAeOHC8Cvn6dNw7dkAgAICpXTvCLh5eMQPT1K5dtVp1ziIHWT/voHDbYXy5pZhcx2dYhgPhNMMXsFNMEflaESrMiTUlivjz2pLQNSlkZ1iGEkl+Qoh64fT4+cv/Mvg8/RCXd2/BlBt7YTdX/ghyFB9k+uLJTCvajFPTuNKSxL0XvURSi15BiPpkgbIynBs3HZ+BuWED/vx8AHR2O9ZePQm/e2J5F2bPnuibnb5L0ev0kLl6J/mbs/AeLsboVNiVFbshAqumw0o4fmWjlCKKKCXf4sScFEnzXq1JSW14MyxDiXznxGntXLeHjcu30iwugvOv7Y/F1nAmFYjQkVVQxsSP17LlcDEPXdaJey9qV2l8z+3MZ87iB3k/dxUFeh0Xm2K578JnaNdqaFBiVkrhzco6aQ9M1/bt4PcDYGrThrChQyvG6szt26GdosXl9/g4vGEfuRv24zpYgL7Ej01Zync+0QwVMywdqhgHTgpNLowtIojumkzSgHaYbXW7kXRTJFUdRJX8fj/PjX2NVQvWEfD7MRgN6A16Xl7yJO17twl2eKIB+Xl3HpNmrsPrC/Da2FSGdz55dw+fp4wvlj3K21mLOaLXGKDZuX/AX+nR6Zp6jTPgdOLatOnYDMzyySn+vDwAdDYblp49sab2Kk92vXphiKpcCcLv93N02yGOrN2N80ABWrEHq99EuD4Co+74EgGHr4RS5cBjURia24nslEjKwA5Yo85+Vmd9UUqB52eU60tAQ7NcDaYBIb2eskZVHTRN+yPwiVKqoNYjEyFr8cfLWb1wHe6y8irdXrcPgCeve4kZe/8d0v/gRWhQSjHtp338/euttI6xMfW2vrRrfvzDXfl9LPrxWd7YPZ99euimGRm2rQ9b55h5P2E1Nz2cyPnX1M1enEopvAcPnbQHpmvbNvCV/zs3tWpF2IUXYu19rFXXoUOlVl3B3mwOrt5N6d6c8hmWXgPh+gjMx2ZYQmzFHpbZuvI9LMPbx5M0oAPJiWe3QXcoUMVPgeszUE5AQzm/BtuNaBGPBTu0c1Kdbs8WwBpN09YBHwDfqsbQXBSnteD9xbgc7krHi/NK2LfpAG16nH1ZGdF0uLx+HvtsE/9dm8UlXeJ59aZehFvKq6OrQICf17zBa5s/YIs+QDtNxz8SRvPh2B0sLynB7yviyL4cnr/lNW5/5iZueOCqGscTcLtxbd5csYi8LD0df+5RADSrFWuPHsTceeexLsxeGKKPJ6eSI4Xs/WotxTuP7WHp0RGmhWE12GkGNCMGj3JTohWTq8uHSDP21s1J7NeW5PYJNY49FCjvZnB+Crh+PQI4oWwuyjoazdgxiNGdmzMmP6XUY5qmPQ5cCvweeFPTtLnAf5RSu+s6QBEcfm+gyuOaTsPvq/o1IQAOFzm5++O1bMgq4v6LO3D/xR3Q6cp7CjZkTOe1ta+xRuchUYNnW1/LlRc8wb//bxrOko0n/dtylbmZ9uQcrrz70rMaa1ZK4Tt8+PgMzPQNuLZuBa8XAGNKCvZBgypmYJo7dkQzGHAVO8n8ZTsFH/+CP+fXGZY27IaIk2ZYlvw6w9JehiUlmvjz2tCqW3LjnmHpXg54q3jBB+4foDEmPwCllNI07QhwBPABUcB/NU1bpJR6uC4DFMEx4rYh7Nt8AHeZ56TjZpuZtr2k1SeqlrYvn7tnrMPp8fHuuD5c1q18bdnOHV/x+s/PsgwH0Sj+kjCcGy96HpOpfJlD+tLN+H3+SvfT6XVkbT902nHmgMdzrFV3vAvTl5MDgGaxYO3enZg7bq+YmKJs4RxM282RTZl4Zm/GULYJu7ISdsIMy4CyUUIRxTjItzgxJTYjrldLknq3wWgy1sF3LsRpFsrTxW//jgyg1X79v/pQnTG/PwG3A0eB94GHlFJeTdN0wE5Akl8jNGr8JSyf9ws71u3BVerCZDGi0+t4bPZkdGe5+4RoGj5ZtZ+nvthMUjMrM8cPoGN8OJkHVvLv5Y/ytS8Pu4L7Yvtx68VTsNlO3rQ6rmUMB7ZmVbqnz+OjWXzkSce8R46cPANzyxbUr626pCRs/fph7d0bU4/uFLqs5G7KwpVVgP6bo9gWrDi2vZeBeKwoZaFUFeOgjCKTC2OLcKK6JJHUvx0ppi3g2wn6lmDqT/lHXhNluRxKXqniBQWWkfUeTm2oTssvFvidUmr/iQeVUgFN066sm7BEsBlNRl7+/knSvt3AhmWbiGrRjItvGUJUXOSZLxZNiscX4MkvNjNr9QGGdmzO62N64ynZyrOzH2a+Kws9cEdEF+66+BUiI1OqvMdND1/LxhVbT+ppMJqNnDesC/bcTPIXfVmxiNx35AgAmtmMpXt3Im65FU9KJwrLrBTlONCKPFjXGQnPOIpRZ6I5RiCOMlVCCQ5KDLno4+w065hI0sAOpESfPMNSBUpR+eNQjt2gAqDpQJ8E0Z+g6ZrmVmCavgUq8mUoehi0Y927KgCRU9D0scEN7hzJUgchxDnLKXZxzyfrWLu/gHsuaseE/hamff8An5TswKfB76wtmTh8CnHNu57xXt9+tJQZf55Kir6EFH0J3VpAc38heMoToiExgUCPAZRGt8ej7GhODYvPQLguHLP+eNeby++kJFCM2+hHiz42w3JgByKrOcMyUPQkOOcDJ3b5G8FyKbpmr57Nt6fRUYFS8Kws/8J0IZoutDe/rtFSByGEqEp6ZiETP06j2OnjtRtaknvwOa74fC2lmsblpuZMGvIcLVMGnfJ65fXi2ratovuybfp6Hm1+GAB3WDTu7kPIi2qFTheGRZkJ04VhN4Rhd5Zf71FuSlQxR3UFqAgHtlaxJPRrR3LHxJq9MdcXnJz4ALzg+halAk26+1PThYHlsmCHUSsk+QkhztrctEwe+3QTieEB7jzvc17N+Ik8vcZQQyT3DXqCTu0rjwP5cnOPl/BJ34Br0ya8fo2SFl3wJnRG33Uslp7hhBnCiDUc7173BbyU+Isp1IrJs5VhTYkmrncbWvVIqZsZlsp3ihcClE/xF42BJD8hRLV5/QH+8fVWpv+0i2tafsFWyyreLNboo7Pw6nkP0LvHzcCxVt32HRWzLx3pGeR5w3A3b48uMgFL9GWEXTKaSEMkUcdaUgHlp8RXTAllFJhdmBMiienRkqR+bet3hqX5InAvojzZ/Up3bDeTRrycoYmRMT8hRLXklbq5Z0YagfxpeOJ+Zq8RugT0/KnHeAa0Go0rI4PStevJ2ZxJocMIYS0w22Kwm5sRbiwvoArl6/Ac/mJKVRleq4axRRhRXZJJHtAec1jwp80r/xFU3vUQKAWcx6byW9Bi5qIZZJlPQ3K6Mb+QTX6apo0EXgP0wPtKqRdOda4kPyHq1qaDRbww6zkc4QvZYwow4KDGTWWDsZW0wevUYzBGYDc3I8LUDKPu+IL0Ml8ppcqB2xxAH2sjsmMCyQM7YI8N8QK0AQfK+RX4NoOhI5r1GjRdeLDDEmepwU140cr7Ft4CRgBZlG+v9oVSaktwIxOi6Zn5weuUrN/K9d4kYrRJhJuaEW5qhkVvAztgB7ffSbG/mGxdHlqUlfD2ceV7WCbHnPH+oUjT2dHsNwU7DFGHQjL5Af2BXUqpPQCaps0GrgEk+QlRRxz5pWT+uJXCjH0E8p2Y/UbCjZEMMfQGe28AvAEPxd5CjvqPosKM2Nu1IKF/exLaxTfu7b1EoxOqyS8JyDzh6yxgwIknaJo2AZgA0LJly/qLTIgGzl3qJGvVLgq2ZuE9VF5AtXwZQQRhmkYYLfAbfBSrAnJd2ZQEthPZMo6U4X1p2ae9JDnRKIRq8quqXs5Jg5NKqanAVCgf86uPoIRoSLweL4fW7CFvcybug4UYHH5sykq4IQKzpqcFYQSUjVKtkGJ3PjmOfZT6j5Bu28byltuJ10dxae9nuWXokGC/FSFqXagmvyzgxH2QkoFDQYpFiJDm9/vJ3phJzvq9ODPz0Zf4sAbMhOsjMOiMxGEG4ilVRZR4Cskr2k/AmU9YmJ+EHok0796KH3K+5iPHZlwa9CqJQCudxP233MR5LSsXbBWiMQjV5LcG6KBpWhvgIDAGuDm4IQkRfHm7DnNw9W7K9h1FK3Jj8RoJN0Rg0pmJxUD5HpallPiKyCo5jK80B0P+fqJtLqJ6dSIhNRVr6mCMLVvidhUx6DxogAAAIABJREFU5/uHeW/vmxTpNAb6w9iReS2l8YOZNqkPcRGWYL9dIepMSCY/pZTvWAX5bylf6vCBUmpzkMMSot4UHcrn4C87KdmVjcp3YvbqCdeFY9HbiAaiicEdcFKiijhSth9vSS66o/uwZ28hTO8hLjUV63mpWFMvw9K9B/qw43sw+rwu5i95iLcPfEOOXmOgZiPMdSOf7u3KmH4pPH1NN8wGGdcTjVtIJj8ApdQCYEGw4xCiLjkLSsn8eSeFOw7hz3FgcmuEa3ZshnAigAii8AbslKgijvpy8bpK0QozMe/4//buO0yq8uzj+PeZtjuzFZalLr0XYSkC9oIaYtdExW4kArFFo8b2GqPGFpSoqFGMiYoi9t4bYkOlLL2D9LLAsm36nPv9Y1akLEjZnTOzc3+uyyvLOTt7bk7Y+c156g/4yleTYQy5nTvH96o7fSje4hvxtGuHMbt2m1tWjI+/uZdHl7zMCofQ27i5rv1wRv/Qi1Vb/Pzj9J6cP6hNra9VqqFJ2vBTqiEJ+YOsnrKYsvlriKyvxO0XsvCR7crFZww+colZ8Q1Ut1JJKVtwUIV3yzIyZn6Nu7yMPMCRl4e3T2+8F5+Br7iYzN69cWZn7/HaIsLXUx9n7JynmO+I0UkMj3S+kHDOBVz36mwy3REmXDaYge33btcDpRoCDT+l6lAsHGXtjGVsnLmS8NqtuKosfJJBtiuvZoRlFpZ4qaKCKuNna0YAV56bbLef7A0L8S0owbl4MYiAMWR06oj3hCF4i/vi7Vscf6rbh82EZ8yZyEPTHmA6IVqJcE/rUxh6+N95bPIKHnp7Br2L8njigv60zLd/WTGlEknDT6n9EIvFKJ23hg3Tl+FfuQVHRQSf5anZJdxNs5oRltVUUIWfCtdGnE2zye/YhMbeEL7Fa+Jb+Xw9k9jWrQD4c3Lw9ulDzgknxJsx+/TGmbN/S2otXPIhj3x3J5OtSprEhFtbHMXvjr2foGRw+cSZfDJvA2f2a8U9ZxxEplv791T60fBT6ldsXraetd8vpfqnUtgawht1k+OMj7AswEkBhQSkmiqqWO/chCnwktelOa0GdaJpsIJASQn+GTMIfDGT0JOLKLXiuwV4OnYke8ix+Pr2xVtcjKdDh316qqvNytVTeHTyLXwY3ki2CH9u3J/zjnsAX1ZTlpVWMWL8tyzfVM3fTu7BHw6rvW9QqXSg4adUjcp1Zaz+fjEVSzYgm/1khJ1kO7LxOrNoBDSigJAVpFIq2Gg2Q34G2e2b0mpgR4raN8Py+wnMnhPfxueDD1h3bwmxsjIAHNnZeHv3JmfUKLx9i/H27o0zL2/PBe2DjaXzeeLz63gjsBK3wPCcLlwyZAx5jdoB8MWCjVw9cQYuh2H8pQM5tFOTOru2UqlIw0/tFxHh42cm8dLotygvreCgI7sz/J7zaN21ld2l/apAWRWrpiymfOFaoqXVeIKGbJNFliuHHCCHfCKWj0oq2GLKkWw/3jaNad6/Iy26tcTpdCIiRFavjgfd+M9YVjKD0MJFEIsB4Gnfnuyjj443XxYXk9GpI6YelgUrL1/F05/9hQkV84kBv89sxchjRtOkWW8g/v/Tv79cyuiPFtKteS7jLuxP68a+Oq9DqVSTtFsa7Qvd0ijxnrl9Iq89+C5BfwgAYwzenEyemDGaFu2b2VxdXMgfZM2PS9kydzWRdfE1LLPES5YrF0fNBqoxiVIZrcBvgsRynGS2akRhn7a06NMOp+eXz4ZWIEBwzpyanchnEigpIbZ5MwAOn4/MPr3xFhfjKy7G26cPzvz8ev27+f2bGf/59TxT+iPVBk5yFXD5kf+gdZsjfvmecJQbXpnFe7PXcUqflvzzd73xerR/T6WPlNvSSCU3f2WAV0a/QzgY3nZMRAj5Q7x47xv8ZdyohNYTC0dZW/ITpTN/IrSmHGdVFJ9kxlc+MS6a48OSTKqkgmoTYKsniLtFLo17FtHq4I609e24komIEFmzlqqaXcgDM2YQXLgQolEAPG3bkn344fHmy+JiMjp3rpenutqEw35emXQT49Z8zhaH4RhHNlcOvoUuXU/d4ftWbvYzYvxUFm2o5ObfdmPEkR20f0+p7Wj4qX22etFaXB4n4eCOx2NRi3nfLqy368ZiMTYtWMv6aUsJrCzDVITxxjzk1oywjK9h2ZRqKqmimkpXKa7CbPK7taDokC60ycuq9edawSDBuXPjQVdSgr+khFjpJgCMz4f3oIMoGD4cb3EfvMXFuBolfr3LWDTCu9/cxePL3mCtAw7GwyPF19Cnz0W7fO/Xizdx5YvTsSzhf38YyFFdChNer1LJTsNP7bPCogKi4egux42Blh2b18k1ypZvYM0PS6lavhHKQmRGXOQ4c8lwZlKAE2iy3QjLUkyBj9xOzSga1JmiFrsPJxEhum7dtpALlMwkOH8+RCIAuNu0IeuQQ7Y1YWZ06YJx2fdrIpbFZz8+xNh5z7LMYdFDHNze9Q8cMvDqXUaGighPf72ce96fT6em2Yy7cADtmtQe+EqlOw0/tc8aNctn0Mn9+f7daYSDkW3HPV4Pw24+Y59+VuX6raz+fhEVi3cdYZkP5FNA2ApSYSoodWyJj7BsV0iLgR0p6vDrQWuFQgTnztv2VBcoKSG6cSMAJjMTb69eFFxySbwJs08fXAXJs/P4lJnP8PCMscwxYdoJPNjuTI4//DaMc9df22Akxk2vzeLNkrUM7dmcB87uQ3aG/nortTv626H2y43PXskjV/yHL178BoDcxtlc9dgf6TG4S63fHyivZvX3iymbv5bYxqqaEZY+sly520ZYRq0sKilnC+VYWX58bRrTrF8HWvRotdcbqEbWr6/pp4sHXXDePOTnp7qiInwDB24bgZnZtQvG7a6T+1GXZi98m4e/v4fvpZrmlnBnq+M45Zh7cHlqH6W5ZmuAkeOnMndtBdcd34UrjumEw5Fe/XsiAqGPEf9LICHIPBXjOwNjPHaXppKUjvZUByQUCOGvCJDfNA9jDJFAmNVTl7Bl9mrC68pxBwSfeMneYYRljKpoOdUmRCzbSUarPAr7tKNl8Y4jLH+NhMME58+PTyCvGYEZXb8eAJORQWavXvhqBqV4+/TBVZjcfV9LV05m7Fe38Vl0C41iFpcVDuTsIQ+Q4dv90+iUZZu54oXphKMW/zqnmON6JMdI20Szyv8OwTdAAvEDxguuXpjGz2GMjnBNVzraU9W5WDjKulkrKJ25guDqsvgISys+wtJtXDTDi/w8whI/5Z4g7ua5NO5RRKtBHcnw7ftecZENG3dovgzOnYuE4yNO3S1b4uvXLx50fYvJ7NoV40mNT/1r1s/k8Ul/5d3gGrwiXJ7Xk4uGjCErr2i3rxERnvtuBXe9O482BT7GXTiATk33vMB1QyXR5RB4DQhtdzAA0bkQmgSZQ+wqTSUxDT+1R7FYjE0L17F+2jICKzZhKiLx5b1cebgdHgpxA02plu1HWGaR17UlrQd3pnWj/XtDlnCY4MKFBGbM2DY4Jbp2HQDG4yGzZ08anX9+TV9dMe5mTevwb50Ym8qW8dTn1/Fy5WIcwIXetgw/9gEaFXbf4+uCkRh/e2sOL09dzZBuTfnXsGJyM5Ov+TZhwj8AtTTzih8JTcZo+KlaaPipbcpWbGTN90uoWl4KZcH4CEtHDhlOLwUYoJCg+Kk0lWxwxNewzOnUjFaDOlPU8sC2w4mWltaMvqwZgTlnDhKKf5J3tWiBt7gPvosvjvfVde+eMk91tamo2sAzn13H81tKCBs43d2UUUfdS/OiQb/62vXlQUY9P42SVVu5+thOXHNcl7Tr39uFIx+ME3bpwXGDI3kGMKnkouGXhio3bmXNlMWUL96AtdlPRshBtsnC58quGWHZmLAVopIKSh1lkOfH17YJLQZ0oKhLywO+vkQiBBcu+qUJc8YMImvWAGDcbjJ79KDRsGHbJpG7m9fN9Am7BUIVvPjFjTy97isqHIahzlyuOPR22nX6zV69ftqKLYx6fjrVoShPXNCPob1a1HPFKSLjaGp/K3NivGcmuBiVKjT8GrBgRYBV3y9i6/y1RDdW4glAlski25VLNpBNHtGaDVTLTAVbsgJktm5Ms37taduzaK9HWP6a6ObN2wVdCYE5c5BgfIa8q1kzvMXFNLrgArzFfcjs2RNHCj/V1SYSDfHG5Nt5YsV7lDrgcDK5uv/1dO81bK9/xos/rORvb82hZb6X54cPomvz/dvqqCEyJgMaP4uUjQSp4ucmUJP3AMa1+35Tld40/BqASDjC2h+Xsmn2SsLrKnBVW2RJJtmuPLzGgZdsLPFSSTmVVFGWESCjVT5NDmpDq/7tcXvqrr9IolFCixb90oQ5o4TIqlXxk243md27k3/2WfE1MPv2xd2i4T69WFaMD74bzWOLXmSVw6KvOBndcyT9+4+KrwiwF8JRizvemcsL36/kyC6FjB3WlzxfGvfv7YZx94DCLyE6Jz7Vwd1HpzmoPdLwSyGxWIx1M1dQWvJTfIRlZRSflUGOMw+no7YRliHczXNo1L0VRQM70Sa77nfrjpaVbZtTFyipearz+wFwFRbGn+pqmjAze/TAkbnvozxTjYgwefo4Hpn9BItMlC4Cj3UYxhGH3rRPa4BurAxy+fPTmbqijFFHdeSG33TFmcT9exLbhFT/B8LfgKMZJuuPmIzBCbu+MQ5w907Y9VRq0/BLIaXz1sDLa7aNsPRLFVVUsdZVirMwi7wuzSk6pCutG9fPkHeJxQgtXrzDJPLwihXxky4Xmd26kX/mmfGlwfoW42rZMu0WU54672Ue+fEBZhCgtSXc3+a3DD3qLhzufQv9mau2MnL8NLYGwjxybl9O7XPgfa31SWKlyKZTQCqBCLAQCf+I5N6Cw3eO3eUptQsNvxRS2KMVU5vMJKdDU1oN7ERRUf2OZIuWlRGYOfOXEZizZmHVPNU5mzTBW9yH/LN+Hx+B2bMnDm/dP1mmivnLP+Xhb/7ON7FyCmMWtzU7jDOGjMadue8b1r46bTW3vDGbwuwMXvvTofRsWXeb3tYXqR4HUgFsv+ZrACrvQ7y60opKPhp+KcTpdDLo+lPq5WdLLEZoydIdJpGHly//+cJkdu1K3umn/zICs6go7Z7qavPT2qk8+uVNfBTeQG7M4tpGfTj3uAfx5ux7X2YkZnH3e/N55tufOKRDAY+d34/GWSkSGqGv2TH4thNdAu4eCS0nkcSqAOPBmIbfpN+QaPilqVh5OYFZs2qaL2cQmDkLq7oaAGejRnj79iXvjDPi2/j06oXDp7t/b2/95oU88fn1vFm9HI8Il2V15JIhY8gt6LRfP29zVYgrJ8zgu2WbufSw9txyYjdcTsevvzBZOAshtnTX4xIFx4HNAU1WEp6OlN8CsZWAQTKOxeTdjXHk2l2a2gsafmlALIvw0qU7TCIPL615o3I4yOjaldxTT4mPwCwuxt2mjT7V7UZZ5Rr+89l1TNw6Bws4x9OCy465nyYt+u33z5yzppyR46dRWhXiwbP68Lv+qTc832QNR8IzgcB2R93g6YdxNox5mtuT6Cqk7FIQ/y8HQ58jZX/EFLxsX2Fqr2n4NUCxykoCM2dtm0AemDULq7ISAGd+Pt7iYvJOOSW+DuZBvXBk6Z5vv6Y6UMZzn9/AsxunEDBwsiOfy4+4i1btjzmgn/tWyRpufG0WjXweXh11CL2L8uuo4sQyGUchOddC5b9qVluJgLsYk/+w3aXVC/GPj/8ddxCByEIksgDj7mZLXWrvafilOLEswsuX79BXF1qyFETiT3WdO5N74ok12/j0wdOunT7V7YNQxM/Lk/6Pp1Z/QpkDhhgvVw26mY7dD2zlkJgl/PPDBTw5eRkHt2vE4+f3pzAno46qTgyxKiH4PlhbwHMwxncxeM+G6GJwNsE4W9ldYv2JLiU+qnUnxgmx1aDhl/Q0/FJMrKqK4KxZvzRhzpyFVV4OgCMvD2+f3tvCLvOg3jiz9aluf0RjEd759l4eX/Iq6x3CIHHy54Ou5KC+w/d6gvrubPWHuerFGXy1eBMXDm7LbSf3wONKof49QMIlSNkfACs+qdxkgHsgptG/MZ4+dpdX/zz9axbUDu14XCLg6mpLSWrfaPilkMi6dSw5dkj8qc4YMjp1IveEE7Zt4+Np1w7jSK030WQjInwydSxj5/6Xn0yMXpbhrq4XMXjwdVAH93bB+gpGPDeNdeUB7jvzIIYNbFMHVSeWiIVsvRKkeruDgXgYBF4H39n2FZcgxncuUv1sTdOnVXM0EzKPw7ha21ma2ksafinE1bw5hddeS2bPHnh798aZo+s71hUR4bs5L/Dw9IeYR4iOMeGh9qdy7JF/x7jqZrrBB7PXcd0rM8nOcDFxxCH0b9uoTn5uwkUX1KyhubMAEngVkw7h52gETd5AKh+E0GQwPvCdj8m61O7S1F7S8EshxhiajLjM7jIanJlL3ufh7/7Bj1YlLaMW/2h5NCcfez/OjLpZKceyhDGfLOLRL5ZQ3DqfJy/sT7NcnROW6oyzJSb/QbvLUPtJw0+lrcWrv+WRybcyKbKJxjGLmwr6c9aQB/Bk193GuBXBCNdMLOHzBRs5e0ARd53eiwxX3eyWYRtXNzDZOw7zB8ALmb+zpSSl9pWGn0o7qzbO4fFJN/CefxVZIlyV3ZULjv8Xvvy2dXqdJRurGPHcVFZu8XPXaT25YHDbBjHS1hgH5I+Nz3PbNuAlE9wDMD4NP5UaNPxU2ijd+hNPfn49r1UswCnCJZmtGX7saPKaHVTn1/p03gaueamEDJeDF/44iEEdGtaO4sbTN76F0HZTHXAPaBDhrtKDhp9q8MqrS/nf59fzwuZpRIEznQWMPPJumrY9vM6vZVnCo18sYcwnizioVR5PXtiflvkNc8Fv48gF395vyKtUMtHwUw2WP1TFhEk389+1k6gywm/J4opDbqNN15Pr5XpVoSjXvVzCR3M3cEbfVtx75kFkulO8f0/tN7EqITITHLngOkifipOMhp9qcCLRMK9+fQfjlr/DJodwlLi5qt+1dO19wQFPUN+dnzZVc9lzU1m2qZrbTu7BpYfpSjrpzKp+BiofBOMGLHA0gUb/xbhSb15nQ6XhpxqMmBXj/R/G8NiC51ljLPpZhjHd/0jfQVfXW+gBTFq4katfnIHDYXju0oEc1qlJvV1LJT8J/xhf45RQfDAQQGw1UjYcmnysH4qShIafSnkiwhclTzN21r9ZQphuMeHxTr/n8MNuxbjc9XrdJ75cxj8/WkDXZjk8ddEAWjfWrZ/SnVSPZ8fdLQAssEohOg/cPe0oS+1Ew0+ltB8XvM5DP9zPLPHTNmoxuuh4TjjmHhye+g0hfzjKDa/O4r1Z6zipdwtG/743Po/+Oinio19r5QBra0JLUbunv60qJc1d8QWPfH0730bLaBqzuL1wEKcNeQC3r/43Tl21xc9lz01l4YZKbhzajVFHddCmLPWLzOMhMgsI7nhcouBOg0W/U4SGn0opy9ZN59Evb+KT0DryYxbX5/XknOPGkJmXmA1gv1myiSsmTMeyhP9dcjBHd6271WBUw2C8ZyH+l+JbGxEEDJAJOTdgHHWzZJ46cBp+KiWs27KYf39+A29VLSFThFHedlw85EGyCxOzb5qI8N9vfuKe9+fToUkW4y4aQPsmul2U2pVx+KDgVSTwGgQ/AWcBxncBxtPf7tLUdjT8VFLbUrmOpz6/npfKZgJwnquQy46+n8ZFAxNWQzAS45bXZ/P6jDWc0KMZY84pJjtDf3XU7hmHD5N1IWRdaHcpajf0N1glparAVp6ddCPPbfiGIHCayeFPh91Bi04nJLSOtVsDjBw/jdlryvnL8V248phOOBzav6dUqtPwU0klGAnw0ld/4z8rPmSrA463Mrhy4A106JX4ZbS+X7aZy1+YTihq8dRFAzi+R7OE16CUqh8afiopRK0ob353P/9e/DIbjcWhloOre42i54BR9TpBvTYiwvNTVnDHO/No09jHuIsG0KmpDlRQqiHR8FO2ssTi4+lP8Oicp1hBlN4x4b4u53HwoX8FZ+L/eYaiMf725lxemrqKY7s15aFhxeRm1t9E+USQyHyk+j8QXQ6evpisP2KcLewuSylbafgpW4gIX8+byNhpY5gvQTpFLB5pdxJHH3Unxm3PLucbKoKMen4aM1Zu5cpjOvGX47ukfP+ehL5Gyi4HwoAF0QVI4E0oeA3jamdzdclJxAJrE5js+MhN1SBp+KmEm7H0Qx767i6mxypoFY1xT7MjOHHI/Ti9+bbVNG1FGaOen0Z1KMrj5/fjxINS/8lIRJDy29hxsnUUpBqpfADT6FG7SktaVuAjqLwDrEpAkMwTMXl3YEzD3JYqnWn4qYRZuHoKj3x1C5PDpRTEYtySX8zvjx+DO6e5rXVN/GElt701hxZ5XsYPH0i35rm21lNnZCtYG2s5YUH4+4SXA/FAxloHuDDO5FogQMLTofwGdviwEPwAkQCm0VjE8iPV4yDwVvyc9zRM9kgNxhSl4afq3crSuTw66a98WL2CbEv4c1ZHzjtuDL6CTrbWFY5a3PXuPMZPWcERnZsw9ty+5Ps8ttZUp4yP+OoitXAkPuAlMhvZ+heIrQcEcXXC5D+McbVNeC21kaon2WVJMkIQmoQV3Qhb/wTRRfFjANVPI6GvoOAVjHEkuFp1oDT8VL3ZWL6SJz6/njfK5+ES4VJPC/5w7GjyWvS1uzRKK0Nc8cJ0fvhpCyOP7MANv+mKy9mw3sCMyUAyT4TgB2x7wwbAC75LE1qLWGXIlotAqn85GF2AbDkXCidhTBJ86IitrP24cUPoY4gtZcf7GIofC38DGUckokJVhzT8VJ0r92/m6S9uYELpD8SA35l8Rh51F4Xtj7G7NABmry5nxPiplPnDPDysmNOKW9ldUr0xuX9HrK0Q/i7+Ji4R8J2N8Z2X0Dok8BZIbKejFkgAQpMgM7GLF9TK0w8CPwE71SnR+HZEsvM2RcSPRWZr+KUgDT9VZ/yhKp6f/H/8b81nVCOcJF4uH3QLrXucYXdp27w+fTU3vz6bJtkZvDrqUHq1yrO7pHplHD5M43FIbC3E1oKrA8ZR/ztf7CK2ll2bFImHcWx9wsupjckahQQ/APEDVs1RL2RfBiYPyGDXv4MXdNpIStLwUwcsHA3xyrd3M27Zm2wxwtExJ1f1vYoufS9N+AT13YnGLO79YAFPf72cwR0a89h5/SjIzrC7rIQxzpbgbGnf9T39kMDLNcGyPWfSbPNjXK2h4DWk8l8Q/gGcBeAdBsEPIVICRHZ+BRgPZPzGjnLVAdLwU/stZsV4d+ojPD7vWdaaGAdH4eEeF1M8+FpwOO0ub5st1WGunDCdb5du5g+HteOWE7vjbmD9e0kvYwg4W8cn2hOuOZgJnv7g7m1nZTswrvaYRo8ANVNFNp9UU/POTbZucHXG5D+ocwFTlIaf2mciwuezn2VsyaMslRA9oha3dziDQ474P9smqO/OvLUVjBg/lY2VIUb/vjdnDWhtd0lpyRg3NJ6IVD8FwXcAF3jPwmRdlLwbAUdm1TTX7hx8TvCeiSPvLjuqUnVEw0/tkymL3uKR7+9jtlVFu0iMB1sdy/HH3IvJzLG7tF28M3MtN7w6k3yvh5dHHkJxa/sm0SswjixMzjWQc43dpewdaz1QWwtBLD4ARqU0DT+1V2avnMzDX/+N7yObaR6NcWfBAE457gFc2ck1URkgZgmjP1rIE18uZUDbRjx+QT+a5iTXE6lKAe5eIOFaTmSCZ3DCy1F1S8NP7dHS9TMY++XNfBZcQ6NYjL/mdOPs48aQ0aid3aXVqtwf4aqJM5i8qJTzB7Xh9lN64nFp/57ad8bZCvGeBoF3gZ+nObjBkY/x/t7O0lQd0PBTtVqzZQmPT/or71YswivC5RlFXDT0QbKa9bK7tN1atKGSy56bytqtAe454yDOG9TG7pJUijO5dyGu3hAYD1YVZB6Pyf4TxqFbXKU6DT+1g01V63nqixt4efMMHAIXOhsz/Kh7aNT2cLtL26MP56znupdL8GW4mDhiMP3b2jCXTTU4xjgwWedA1jl2l6LqmIafAqAyWM7/vryZ59d9RRjhdPEx6tC/0bzryXaXtkeWJTz02WIe+WwxfVrn8+QF/Wmep/17KvWJCOKfANVPgZSBuw8m5yaMu4fdpTUIGn5pLhDx8+LXd/L0ivepMMLQmIsr+v+Fdn0uSJoJ6rtTGYxw7UslfDp/I2f1L+Ku03uR6U6e+YVKHQip+hdUP8u2/sbwFGTLeVDwKsZl76LwDYGGX5qKWBHe+P5Bnlj4IqXG4rAo/LnXH+k+8CpwJP8AkaWlVYx4biorNvu587SeXDi4bfLOF1NqH4lVBdX/Y8eFtAEJIlX/xuQ/aEtdDYmGX5qxxOKDGU/x2OwnWUWE4qjFPzudxYDDbwFXEqysvxc+X7CBP79Ygtvl4Pk/DmJwhwK7S1J7SaLLatYY7YpxFtbvtSQEwY8htjo+bcFzWOpsPRRbDcYFslP4YcUn36sDpuGXJkSEyQte5pGpY1hk+ekSifFY699wxNF3YTJSY+SaiPDYF0t48JNF9GiRy7iLBtAqXzcSTQViVSJlo+I7IBg3SAjx/g6Te3u9BJJEVyBbhsV3XZAgmExwdoDGz6fGcmTO5vFFv3dh4n8PdcA0/NLA1OUf88i3dzEjupXWkRj3NR3Mb48bjcOXOk9M1aEo178ykw/mrOf04pbce2ZvvB7t30sVUn7LL4tDS83OCIE3EVcXTNb59XC968EqY9vuDOKH6GKk+t+YnOvq/Hp1zTjyEe/JEHifHXeSyMBkX25XWQ2Khl8DNn/t9zw8+Va+CW2gMBrjtrxenHH8GNx5RXaXtk9WbK5mxHPTWLyxkv87qTvDD2+v/XspRKxqCH3OrrsiBMD/LNRx+IlVDpF5/LIt0c9CEHgTUiD8AEy0LVGVAAAUAklEQVTunYjJBf9EIALOlpjcv2E8ybELRqqzJfyMMWcBfwe6AwNFZOp2524GhhNfTfZqEfnIjhpT2U+b5vPYpL/yYfVP5MZiXOttz7m/fQBv0+52l7bPJi8q5aoXZ2AMPHfpIA7v3MTuktS+Ej+wmw8rVkV9XHA/zyUXYzyY3FuQnL/G+/6MTz/01SG7nvzmAGcCT25/0BjTAxgG9ARaAp8aY7qI7LIFtKrF+opVPPHFDbxZNgePCJe5mnLJsfeR23qQ3aXtMxHhqa+Wcd8HC+jSLIdxFw6gTUEK9NWoXTmaxP+z1u58ol52QDeOfMTVDaJz2DHsPOA9tc6vV9+MccUHv6g6ZcsdFZH5QG2fYk4DJopICFhujFkCDAS+S2yFqaXMv5mnv7yJFzdMwUI4R3K47Ig7aNLpBLtL2y+BcIwbX5vF2zPXctJBLRh9Vm98Hv3lT1XGGMi7Gym7nPhefjHAE3+Sybm2fq6ZPxrZfG7NaEk/GB8422CytL9MxSXbO0orYMp2f15dc2wXxpgRwAiANm3Scw3H6nAVz311O8+u+pgAwskxD5cPvJFWvc5O+gnqu7O6zM+I56Yxf30Ffx3alT8d1VGbehoAk3EYNHkNqf5ffHNYz0CM7wKMs36asY2rAxR+Ed+F3VoDrl6QcSTG6CApFVdv4WeM+RRoXsupW0Xkrd29rJZjtTbSi8g4YBzAgAEDUqchvw6EYiFe/u5+nlryGmXGYkjUcFXvy+k4YGRKTFDfnW+XbuLKCTOIxCz+e8nBHNM1+bZLUvvPuDph8u5O3PUcPvCdmbDrqdRSb+EnIsftx8tWA9tvtV0E7NxRkLaiVpR3pj/O43P/x3qiDIpY/Lnr+Rx02A3gdNtd3n4TEZ759if+8d582jfJ4qmLBtC+SZbdZSmlGrBka/Z8G5hgjBlDfMBLZ+AHe0tKHhvW/MCdc8bRLWJxV7uTGXz0HeBO7UnewUiMW9+Yw2vTV3N8j2aMObsPOZmpG+RKqdRg11SHM4CxQCHwnjGmRER+IyJzjTEvA/OAKHCFjvT8RavWhzKx22V06X0RxtfI7nIO2LryAKPGT2Pm6nKuOa4zVx/bGYdD+/eUqm8SXQbRFeDugnHWOqyiwTMiqd9dNmDAAJk6deqvf6NKGj/+tIU/PT+NQDjGv84p5oSetXUPK6XqkljVyNYrIDy9Zpm5cHyD3rx/xqdUNDDGmGkiMqC2c6k7OkKlrOenrODccVPIyXTz5hWHafAplSBSeSeEpwJBkEogBMFPkeon7C4t4TT8VMKEojFufn02//fmHA7v3IQ3rziMzs1y7C5LqbQgEoXAe8TnWm4vCP4JdpRkq4b3nKuS0saKIH96YTrTVpRxxTEd+cvxXXFq/55SCRSt+a8WUp3QSpKBhp+qdzNWljHq+WlUBKI8dl4/Turdwu6SlEo7xmQirs4QXbjzGfAcaktNdtJmT1WvXp66inOenILH5eD1yw/V4FPKRib3H/Gl3vh5OpEHTA4m5yY7y7KFPvmpehGJWfzj3Xk8+90KDu/UhLHn9qVRVmrsFK9UQ2U8faDgHcT/LEQWgadvzTJzhXaXlnAafqrObaoKcfkL0/lh+RYuO6I9Nw7thsupjQxKJQPjao3J/T+7y7Cdhp+qU7NXlzNy/FQ2V4d56JxiTu+bnhNo1YGR8Eyk+t8QWQjGA85WGO+pkHkixmgLgjpwGn6qzrw5Yw03vjaLgiwPr/3pUHq1yrO7JJWCJDQJKbsaCP5yMLYcCf8I1c9BwYsYk2Fbfaph0PBTBywas7jvgwX85+vlDGrfmMfO70eTbH1zUvtORJDyO9gh+LYJQXQp4n8dk3VuokvbgYSnIlVPQGw1ePpjskZhXK1//YUqaWj4qQNSVh3mqhdn8PWSTVxyaDtuPak7bu3fU/tLqsDasIdvCEDwPbAx/KzA+1B+E9sCOrACCX4ABa9jXO1sq0vtGw0/td/mr6tgxPipbCgP8c/f9+bsAfrJVx0gk0n8bWk3k7EBHLmJqmYXIhZU3sWOT6YxED9S+S9Mo4ftKk3tI/2IrvbLe7PWcebj3xKOWrw0crAGn6oTxrjBeyawu2ZzL8Z3XiJL2pG1EazaVkOxIKK7r6USffJT+yRmCQ9+vJDHJy2lf9tG/PuCfjTNybS7LNWAmNxbEKmC4AfEnwCF+KRsA1nDMRmH21hcDmDVfs7RJKGlqAOj4af2Wnkgwp8nzmDSwlLOHdiGO07ticeljQeqbhnjweQ/gFi3INE1YG3EEAb3ANsnYxtHFpI5FIIfAaHtTngxWSNtq0vtOw0/tVcWb6hkxPhprC7zc/cZvTh/UFu7S1INnHE0xnga213GLkzeXYiEIDSpZk+8KGSNhMyT7C5N7QMNP/WrPpq7nr+8VILX42LCZYM5uF3yvSEplSjGeDGNxiKxzfGRqc52GIfP7rLUPtLwU7tlWcLDny3m4c8W07sojycv7E+LPK/dZakUJ1ZV/KlJwpBxOMbZ1O6S9otxFoCzwO4y1H7S8FO1qgxGuPalmXw6fwO/61fE3Wf0ItPttLssleIk9DWy9UrAABZUxJCca3FkDbe7NJVmNPzULpaVVjFi/DSWb6rm9lN6cMmh7TBGN55VB0as6njwiX/HE5UPI55DMO4e9hSm0pKGn9rBFws2cvXEGbgchvHDB3JoRx2+repI6EviT3w7CyOBNzT8VEJp+Ckgvqbi45OW8sDHC+nePJcnL+xP68baia/qUhhEajlugdS2lqdS9UfDT1EdinLDqzN5f/Z6Tu3Tkvt/1xuvR/v3VB3zHA7Edj1ufJjMoQkvJ1HEKoPQ14ATMo7EOLLtLkmh4Zf2Vm72M2L8VBZtqOSWE7tx2REdtH9P1QvjbILk3ACVDwARIAbGCxnHgudQu8urF5b/Fai4E3DFW3zFgvyHMJnH2F1a2tPwS2NfLS7lygkzAHjmDwM5sou9q2eohs+RdRHiGYQE3gQJYDJPAM8hDfIDl0RXQMVdxFeCCcVXaQNk65+h6VcYh+53aScNvzQkIvznq+Xc+8F8OjfNYdxF/WlbkGV3WSpNGHdXjPtGu8uodxJ4m9p3p3BA8FPw/S7RJantaPilmUA4xk2vz+KtkrUM7dmcB8/uQ1aG/jNQqs5JgFr7OLFqzik76bteGlld5mfk+GnMW1fB9Sd04YpjOjXI5ialkoHJPA7xvwDsHHQCGUfZUZLajoZfmvhu6WaumDCdSNTiPxcNYEj3ZnaXpFTD5u4L3hPjWzOJn/j2qR7Ivgzj0v0v7abh18CJCM99t4I7351H2wIfT100gI6FOtRaqfpmjIHceyDzVCT4PuDCeE/HePrYXZpCw69BC0Zi3PbmHF6ZtprjujdlzDnF5Ga67S5LqbRhjIGMQzAZh9hditqJhl8Dtb48yMjnpzFz1VauHtKZa4Z0xuHQ/j2llAINvwZp6k9bGPX8dALhKE9c0J+hvZrbXZJSSiUVDb8GZsL3K7n97Tm0yvcy4bJBdGmWY3dJSimVdDT8Gohw1OLv78xlwvcrOapLIY8M60ueT/v3lFKqNhp+DcDGyiCXPz+dqSvKGHVUR274TVec2r+nlFK7peGX4kpWbWXU+GlsDYQZe25fTunT0u6SlFIq6Wn4pbBXpq7i1jfn0DQng9f/dBg9WubaXZJSSqUEDb8UFIlZ3P3efJ759icO7VjAo+f1o3GWx+6ylFIqZWj4pZjNVSGumDCdKcu2MPzw9tz82264nA67y1JKqZSi4ZdCVm3xM2zcFEqrQow5uw9n9iuyuySllEpJGn4ppFluJv3aNuKyI9rTuyjf7nKUUiplafilEI/Lwdhz+9pdhlJKpTztLFJKKZV2NPyUUkqlHQ0/pZRSaUf7/JRSKUWsaqT6MQi8CQhknojJvgbj0EXc1d7T8FNKpQwRC9lyIUQXAeH4Qf9LSOg7aPI2xuhbmto72uyplEod4e8gtoxtwRc/CNY6CH1uV1UpQ6IrkOBHSGQuImJ3ObbSj0lKqdQRmQsS3vW4VCOReZjMExJfUwoQiSJbr4fQZ2DcQAycnaDxfzGOPLvLs4U++SmlUoerNZiMWk74MM7WCS8nVUj10zVPxiGQKpAARBcg5bfYXZptNPyUUqkjYwgYHzu+dZl4IGYOtauq5OefAAR3OhiB0CREAnZUZDsNP6VUyjDGg2n8Erj7Ee+1cYG7N6ZgIsaRZXd5yUv8uztRezNyGtA+P6VUSjGuIkzBBMSqAkSnOOyNjKMh+C4Q2/G4s632+SmlVCoxjmwNvr1kcq4Dkwf83F/qBuPD5N1jZ1m20ic/pZRq4IyzORR+iPhfgsh0cHXE+M7HOFvZXZptNPyUUioNGEc+Jnuk3WUkDW32VEoplXY0/JRSSqUdDT+llFJpR8NPKaVU2tHwU0oplXY0/JRSKoFEQkiarqqSTDT8lFIqASS6FGvzOciGYmRDMVbZKCS22e6y0paGn1JK1TOxKpDNwyBSQnyJsSiEJiNbzkPEsru8tKThp5RS9UwCb9YsIL39BrJRsDbGN+hVCafhp5RS9S26FKhl6yCJQWxFwstRGn5KKVXvjPsgwFfLCQOurgmvR2n4KaVU/fOeBI5cdlxO2RMPPnc/u6pKaxp+SilVz4zxYgpeg8yTwWSDyQff+ZhGz2CMsbu8tGTLrg7GmNHAKUAYWAr8QUS21py7GRhOfEjU1SLykR01KqVUXTLOQkz+P+0uQ9Ww68nvE6CXiPQGFgE3AxhjegDDgJ7AUOBxY4zTphqVUgkgEkAkYncZKs3YEn4i8rGIRGv+OAUoqvn6NGCiiIREZDmwBBhoR41Kqfol4ZlYm05BNvRDNvTF2vpXxKq2uyyVJpKhz+9S4IOar1sBq7Y7t7rm2C6MMSOMMVONMVNLS0vruUSlVF2S6Cqk7GKILiTewxGG4PtI2Si7S1Npot7CzxjzqTFmTi3/nbbd99wKRIEXfj5Uy4+SWo4hIuNEZICIDCgsLKz7v4BSqt6Ifzzs0tQZhshMJLrElppUeqm3AS8ictyezhtjLgZOBoaIyM8Btxpovd23FQFr66dCpZRtoouAWvr5jAuiK8HVKeElqfRiS7OnMWYocCNwqoj4tzv1NjDMGJNhjGkPdAZ+sKNGpVQ9chcDnl2PSwRcXRJejko/dvX5PQrkAJ8YY0qMMU8AiMhc4GVgHvAhcIWIxGyqUSlVT4zvAjCZ7PgWlAmZQzCuot29TKk6Y8s8PxHZbZuGiNwN3J3AcpRSCWacTaDgdaTyPgh/C8YHvvMwWTrgRSWGLeGnlFLG1QbT6HG7y1BpKhmmOiillFIJpeGnlFIq7Wj4KaWUSjsafkoppdKOhp9SSqm0o+GnlFIq7Wj4KaWUSjsafkoppdKOhp9SSqm0o+GnlFIq7Wj4KaWUSjsafkoppdKOhp9SSqm0Y37ZRD11GWNKgRV215FATYBNdheRgvS+7R+9b/tP793+qav71lZECms70SDCL90YY6aKyAC760g1et/2j963/af3bv8k4r5ps6dSSqm0o+GnlFIq7Wj4paZxdheQovS+7R+9b/tP793+qff7pn1+Siml0o4++SmllEo7Gn5KKaXSjoZfijDGjDbGLDDGzDLGvGGMyd/u3M3GmCXGmIXGmN/YWWcyMsacZYyZa4yxjDEDdjqn924PjDFDa+7NEmPMTXbXk8yMMf81xmw0xszZ7lhjY8wnxpjFNf/byM4ak40xprUx5gtjzPya39E/1xyv9/um4Zc6PgF6iUhvYBFwM4AxpgcwDOgJDAUeN8Y4basyOc0BzgQmb39Q792e1dyLx4DfAj2Ac2vumardM8T/HW3vJuAzEekMfFbzZ/WLKHCdiHQHBgNX1Pwbq/f7puGXIkTkYxGJ1vxxClBU8/VpwEQRCYnIcmAJMNCOGpOViMwXkYW1nNJ7t2cDgSUiskxEwsBE4vdM1UJEJgNbdjp8GvBszdfPAqcntKgkJyLrRGR6zdeVwHygFQm4bxp+qelS4IOar1sBq7Y7t7rmmPp1eu/2TO/PgWsmIusg/kYPNLW5nqRljGkH9AW+JwH3zVXXP1DtP2PMp0DzWk7dKiJv1XzPrcSbCl74+WW1fH/azV/Zm3tX28tqOZZ2924P9P6ohDDGZAOvAdeISIUxtf3Tq1safklERI7b03ljzMXAycAQ+WWC5mqg9XbfVgSsrZ8Kk9ev3bvd0Hu3Z3p/DtwGY0wLEVlnjGkBbLS7oGRjjHETD74XROT1msP1ft+02TNFGGOGAjcCp4qIf7tTbwPDjDEZxpj2QGfgBztqTEF67/bsR6CzMaa9McZDfHDQ2zbXlGreBi6u+fpiYHetEGnJxB/xngbmi8iY7U7V+33TFV5ShDFmCZABbK45NEVERtWcu5V4P2CUeLPBB7X/lPRkjDkDGAsUAluBEhH5Tc05vXd7YIw5EXgIcAL/FZG7bS4paRljXgSOJr4dzwbgduBN4GWgDbASOEtEdh4Uk7aMMYcDXwGzAavm8C3E+/3q9b5p+CmllEo72uyplFIq7Wj4KaWUSjsafkoppdKOhp9SSqm0o+GnlFIq7Wj4KaWUSjsafkoppdKOhp9SDYAx5uCavR4zjTFZNXuj9bK7LqWSlU5yV6qBMMb8A8gEvMBqEbnX5pKUSloafko1EDXrb/4IBIFDRSRmc0lKJS1t9lSq4WgMZAM5xJ8AlVK7oU9+SjUQxpi3ie+23h5oISJX2lySUklL9/NTqgEwxlwEREVkgjHGCXxrjDlWRD63uzalkpE++SmllEo72uenlFIq7Wj4KaWUSjsafkoppdKOhp9SSqm0o+GnlFIq7Wj4KaWUSjsafkoppdLO/wOECKDdM8zjJQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7,7)) \n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "\n",
    "plt.scatter(X[:,0], X[:,1], c=y[:,0])\n",
    "\n",
    "for c in cs:\n",
    "    w = w_seq_C[c][-1]\n",
    "    print(\"w obtained for C=\",c,\" : \", w)\n",
    "    x_line=[-20.0, 20.0]\n",
    "    y_line=[-w[0]/w[1]*xi for xi in x_line]\n",
    "    plt.plot(x_line,y_line, label='$C$ = '+str(c))\n",
    "\n",
    "plt.legend(loc = 'upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(in the above graph, some hyperplanes are on top of each other so we can not distinguish them, namely $C=0.01$, $C=0.1$ and $C=1$, $C=5$, $C=10$.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this simple case, we can analyse the influence of $C$. It tells the SVM optimization how much we want to avoid misclassifying each instance and represents a trade-off between:\n",
    "- a hyperplane with the largest minimum margin (low $C$)\n",
    "- a hyperplane that correctly separates as many instances as possible (large $C$)\n",
    "\n",
    "For large values of $C$, the optimization will choose a smaller-margin hyperplane if that hyperplane classifies better the points in the dataset. It is the case of the hyperplanes created for $C>0.5$ in the above figure, where the margin is close to 0 but the hyperplane accurately classifies all the points. \n",
    "\n",
    "Conversely, for small values of $C$ ($C<0.1$), we observe larger-margin separating hyperplane, even if that hyperplane misclassifies more points. Here it only misclassifies one single point, and the margin is bigger than the one obtained with large values of $C$."
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}