{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# COURS OPTIMISATION CONTINUE - Antonin Chambolle - 2020: Examples of stochastic algorithms : coordinate descent, gradient descent.\n",
    "\n",
    "## Classification of handwritten digits using a \"Support Vector Machine\" and stochastic optimization**\n",
    "\n",
    "Digits from the MNIST Database, [LeCun et al., 1998a: Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. \"Gradient-based learning applied to document recognition.\" Proceedings of the IEEE, 86(11):2278-2324, November 1998.]\n",
    "\n",
    "We will show in this notebook how to solve an easy convex problems either by stochastic gradient descent or (proximal) stochastic block coordinate descent. We will address a very simplify task which is distinguishing a subset of one or several digits from another one. \n",
    "\n",
    "### The problem:\n",
    "\n",
    "In practice, the learning set is made of a set of $N$ vectors $(X_i)$ (here, $28\\times 28=784$ small images representing handwritten digits) and of labels $y_i\\in \\{-1,1\\}$. (Initially, the labels are in $\\{0,\\dots,9\\}$ and indicate the value of the handwritten digits. We define two lists \"mylabelplus\" and \"mylabelminus\" of digits, the first list will be labelled as $+1$ and the second as $-1$.)\n",
    "\n",
    "One wants to learn a hyperplane best separating $\\{X_i : y_i=+1\\}$ and $\\{X_i : y_i=-1\\}$ in $\\mathbb{R}^{784}$. That is, we want to find $(w,b)\\in\\mathbb{R}^{784}\\times \\mathbb{R}$ such that $w\\cdot X_i>b$ for $y_i=1$\n",
    "and $<b$ for $y_i=-1$. Assuming this is possible, to maximize the margins around this hyperplane, one solves the convex problem\n",
    "$$\n",
    "\\min_{w,b} \\sum_{i=1}^N h(y_i(w\\cdot X_i-b))\n",
    "$$\n",
    "where $h(t) = (1-t)^+$ (\"Hinge loss\"). This function vanishes only if $w\\cdot X_i>b+1$ each time $y_i=+1$\n",
    "and $w\\cdot X_i<b-1$ each time $y_i=-1$. In case one cannot find a strip (here of width $2/|w|$) separating\n",
    "the positive and negative vectors, then it will try to find the hyperplane which separates, in some sense,\n",
    "most vectors.\n",
    "\n",
    "To simplify (losing a bit of precision...) here we will assume $b=0$ (in high dimension, it is not as\n",
    "restrictive as it would be in low dimension), in addition, to make sure the band in between the two\n",
    "sets is not too narrow, we will penalize also the size of $w$. We end up with the simpler problem:\n",
    "$$\n",
    "\\min_{w} \\frac{1}{2}|w|^2+\\lambda\\sum_{i=1}^N h(y_i(w\\cdot X_i)). \\quad (P)\n",
    "$$\n",
    "\n",
    "It is easy to find a dual of this problem, writing $h(y_i(w\\cdot X_i-b))= \\sup_{z_i} z_i y_i (w\\cdot X_i) - h^*(z_i)$. \n",
    "We get the problem\n",
    "$$\n",
    "\\sup_{(z_i)_{i=1}^N} \\min_w \\frac{1}{2}|w|^2+\\lambda w\\cdot \\sum_{i=1}^N (z_i y_i X_i) -\n",
    "\\lambda \\sum_{i=1}^N h^*(z_i)\n",
    "= \n",
    "\\sup_{(z_i)_{i=1}^N}  -\\frac{1}{2}\\left|\\lambda \\sum_{i=1}^N z_i y_i X_i\\right|^2 - \\lambda \\sum_{i=1}^N h^*(z_i).\n",
    "$$\n",
    "We find as dual problem \n",
    "$$\n",
    "- \\lambda \\min_{(z_i)_{i=1}^N} \\frac{\\lambda}{2}\\left| \\sum_{i=1}^N z_i y_i X_i\\right|^2 + \\sum_{i=1}^N h^*(z_i)\n",
    "\\quad (D)\n",
    "$$\n",
    "where at optimum, one recovers $w$ by letting $w  = -\\lambda\\sum_{i=1}^N z_i y_i X_i$. In practice,\n",
    "only a fraction of the vectors appear in this sum, called the \"support vectors\" (this means that\n",
    "most of the $(X_i)$ will live in the area where $h$ vanishes, and $z_i=0$).\n",
    "\n",
    "Once $w$ is computed, one can evaluate on the test set how many digits are correctly labeled, by comparing the sign of $y\\cdot X$ with the corresponding $y$, for $(X,y)$ from the test set.\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "In practice, $h(t) = (1-t)^+$ so that $h^*(s) = s$ if $s\\in [-1,0]$ and $+\\infty$ else. We will also used\n",
    "a regularized version, as follows: for $\\varepsilon>0$, we let\n",
    "$h_{\\varepsilon}^*(s) = s + \\varepsilon s^2/2$ for $s\\in [-1,0]$ and $+\\infty$ else. Show in this case that\n",
    "$h_\\varepsilon'$ is $1/\\varepsilon$ Lipschitz and that:\n",
    "$$\n",
    "h_\\varepsilon(t) = \\begin{cases} 1-t-\\frac{\\varepsilon}{2} & \\text{ if } t\\le 1-\\varepsilon \\\\\n",
    "\\frac{(t-1)^2}{2\\varepsilon} & \\text{ if } 1-\\varepsilon \\le t \\le 1 \\\\\n",
    "0 & \\text{ if } t\\ge 1\n",
    "\\end{cases}\\quad ,\\quad\n",
    "h'_\\varepsilon(t) = \\begin{cases} -1 & \\text{ if } t\\le 1-\\varepsilon \\\\\n",
    "\\frac{t-1}{\\varepsilon} & \\text{ if } 1-\\varepsilon \\le t \\le 1 \\\\\n",
    "0 & \\text{ if } t\\ge 1\n",
    "\\end{cases}\n",
    "$$\n",
    "We will normalize the images so that $|X_i|=1$. In the dual problem we find that the derivative of the quadratic term of energy with respect to $z_i$ is $1$-Lipschitz. \n",
    "In the primal problem, if we replace $h$ with $h_\\varepsilon$, we find that the Lipschitz constant of\n",
    "the gradient of the objective (which is now $C^1$) with respect to $w$ is at most $1+\\lambda/\\varepsilon$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import gzip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we download the MNIST database. We used the .csv files in https://www.python-course.eu/neural_network_mnist.php\n",
    "and then installed these in the directory data/mnist/: please adapt the code by setting the variable data_path below accordingly. The code in the next cell is copy-pased from the above link."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load MNIST data \n",
    "\n",
    "# training data\n",
    "mnist_path = \"data/mnist/\"\n",
    "train_images_path = mnist_path+\"train-images-idx3-ubyte.gz\"\n",
    "with gzip.open(train_images_path, 'rb') as train_imgpath:\n",
    "    train_imgs = np.frombuffer(train_imgpath.read(), dtype=np.uint8,\n",
    "                                 offset=16).reshape(-1, 784)/255.0\n",
    "    \n",
    "train_labels_path = mnist_path+\"train-labels-idx1-ubyte.gz\"\n",
    "with gzip.open(train_labels_path, 'rb') as train_lbpath:\n",
    "    train_labels = np.frombuffer(train_lbpath.read(), dtype=np.uint8,\n",
    "                                 offset=8)\n",
    "    \n",
    "# test data\n",
    "test_images_path = mnist_path+\"t10k-images-idx3-ubyte.gz\"\n",
    "with gzip.open(test_images_path, 'rb') as test_imgpath:\n",
    "    test_imgs = np.frombuffer(test_imgpath.read(), dtype=np.uint8,\n",
    "                                offset=16).reshape(-1, 784)/255.0\n",
    "    \n",
    "test_labels_path = mnist_path+\"t10k-labels-idx1-ubyte.gz\"\n",
    "with gzip.open(test_labels_path, 'rb') as test_lbpath:\n",
    "    test_labels = np.frombuffer(test_lbpath.read(), dtype=np.uint8,\n",
    "                                offset=8)\n",
    "# remove spurious dimensions\n",
    "train_labels=np.squeeze(train_labels)\n",
    "test_labels=np.squeeze(test_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reduce the size of train set if needed\n",
    "# this might be necessary for the random coordinate ascent which starts \n",
    "# by computing a N*N matrix of scalar products. We convert before this the\n",
    "# values to float32 (4 bytes), yet N=60000 values (the size of the training set)\n",
    "# require 14.4 Gb\n",
    "# depending on the number of figures selected below and your memory you should\n",
    "# adjust the value of r (the coordinate descent on the primal does not need this)\n",
    "# (the speed of your processor can also be a bottleneck, computing the matrix is best\n",
    "# with many cores)\n",
    "\n",
    "r = 60/100 # 30% \n",
    "\n",
    "indices = (np.random.random(train_labels.shape)<r)\n",
    "train_labels = train_labels[indices]\n",
    "train_imgs = train_imgs[indices,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we define the digits we want to label as +1 and the digit labelled as -1 (first arrays).\n",
    "Then we extract from the database the $X_i$ (array mytri), the $y_i$ (array mytrls). We do the same for\n",
    "the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mylabelsplus = np.array([9,1,6],dtype=float)\n",
    "#mylabelsminus = np.array([0,2,3,4,5,6,7,8,9],dtype=float)\n",
    "mylabelsminus = np.array([3,7,4,8],dtype=float)\n",
    "\n",
    "# extract the labels of interest\n",
    "indicesp = np.zeros(train_labels.shape,dtype=bool)\n",
    "indicesm = np.zeros(train_labels.shape,dtype=bool)\n",
    "for l in mylabelsplus:\n",
    "    indicesp = np.logical_or(indicesp,train_labels==l)\n",
    "for l in mylabelsminus:\n",
    "    indicesm = np.logical_or(indicesm,train_labels==l)\n",
    "indices = np.logical_or(indicesm,indicesp)\n",
    "#indices = np.logical_or(train_labels==mylabels[0],train_labels==mylabels[1])\n",
    "\n",
    "mytrl = train_labels[indices]\n",
    "mytri = train_imgs[indices,:]\n",
    "mytrls = np.empty(mytrl.shape, dtype=np.int8)\n",
    "mytrls[indicesp[indices]] = 1\n",
    "mytrls[indicesm[indices]] = -1\n",
    "\n",
    "# same for test\n",
    "indicesp = np.zeros(test_labels.shape,dtype=bool)\n",
    "indicesm = np.zeros(test_labels.shape,dtype=bool)\n",
    "for l in mylabelsplus:\n",
    "    indicesp = np.logical_or(indicesp,test_labels==l)\n",
    "for l in mylabelsminus:\n",
    "    indicesm = np.logical_or(indicesm,test_labels==l)\n",
    "indices = np.logical_or(indicesm,indicesp)\n",
    "#indices = np.logical_or(test_labels==mylabels[0],test_labels==mylabels[1])\n",
    "\n",
    "mytel = test_labels[indices]\n",
    "mytei = test_imgs[indices,:]\n",
    "mytels = np.empty(mytel.shape, dtype=np.int8)\n",
    "mytels[indicesp[indices]] = 1\n",
    "mytels[indicesm[indices]] = -1\n",
    "\n",
    "# we normalize the l2 norm of the X's to 1. This is interesting as it makes the quadratic\n",
    "# term in the dual problem have a lambda-Lipschitz derivative wr z_i\n",
    "mytri = np.divide(mytri,np.reshape(np.linalg.norm(mytri,axis=1),(-1,1)))\n",
    "mytei = np.divide(mytei,np.reshape(np.linalg.norm(mytei,axis=1),(-1,1)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 3 1 4 3 8 6 9 9 1]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAMe0lEQVR4nO3dYYhc9bnH8d8v3kZhUySa1YY0dNuqeKXYNAxB8FKUcqtGJFawNC9KFGkqKLRYQbEvEvBNuF5bCl4KiYampbUUWzWC3BsJBS1CySq5Gg1qKmmTGLITgtQg0uo+fbHHso07ZyZzzpkzyfP9wDIz5zlz/g+H/e2ZmXNm/44IATj7LWq7AQCjQdiBJAg7kARhB5Ig7EAS/zbKwZYtWxZTU1OjHBJI5eDBgzp+/LgXqlUKu+3rJf1E0jmSHo2ILWXrT01NaXp6usqQAEp0Op2etaFfxts+R9L/SLpB0hWS1tu+YtjtAWhWlffsayQdiIi3I+Jvkn4taV09bQGoW5Wwr5B0aN7jw8Wyf2F7o+1p29PdbrfCcACqqBL2hT4E+MS1txGxNSI6EdGZnJysMByAKqqE/bCklfMef1bSO9XaAdCUKmHfI+lS25+3vVjStyTtrKctAHUb+tRbRHxo+25J/6e5U2/bI+K12joDUKtK59kj4llJz9bUC4AGcbkskARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSVSaxRVo0rZt20rrd955Z2l9dna2Z+2NN94ofe5ll11WWj8TVQq77YOS3pP0kaQPI6JTR1MA6lfHkf3aiDhew3YANIj37EASVcMeknbZfsn2xoVWsL3R9rTt6W63W3E4AMOqGvarI2K1pBsk3WX7q6euEBFbI6ITEZ3JycmKwwEYVqWwR8Q7xe2MpCclramjKQD1Gzrstidsf/rj+5K+LmlfXY0BqFeVT+MvlvSk7Y+386uI+N9aukIKu3fvLq3fc889pfVFi4Z/YVr83qYydNgj4m1JX66xFwAN4tQbkARhB5Ig7EAShB1IgrADSfAVV7TmzTffLK1/8MEHI+okB47sQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AE59nRqNdff71nbfPmzZW2vXr16tL6rl27etYmJiYqjX0m4sgOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0lwnh2VHDhwoLS+du3anrUTJ05UGnvLli2l9fPPP7/S9s82HNmBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnOs6OSRx99tLR+6NChobd9yy23lNavvfbaobedUd8ju+3ttmds75u37ALbz9l+q7hd2mybAKoa5GX8zyRdf8qy+yXtjohLJe0uHgMYY33DHhHPSzr1usZ1knYU93dIurnmvgDUbNgP6C6OiKOSVNxe1GtF2xttT9ue7na7Qw4HoKrGP42PiK0R0YmIzuTkZNPDAehh2LAfs71ckorbmfpaAtCEYcO+U9KG4v4GSU/X0w6ApvQ9z277cUnXSFpm+7CkTZK2SPqN7Tsk/UXSrU02ifa8//77pfWHHnqotL5oUe/jyYUXXlj63AcffLC0jtPTN+wRsb5H6Ws19wKgQVwuCyRB2IEkCDuQBGEHkiDsQBJ8xTW5d999t7S+bt26xsbuN2Xz5Zdf3tjYGXFkB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkOM+e3AsvvFBaf/HFFytt/9Zbe3/7+bbbbqu0bZwejuxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kATn2c9ye/bsKa1v2LChtN7PTTfdVFrftm1bz9p5551XaWycHo7sQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AE59nPAmX/+/2qq65qdOxLLrmktD4xMdHo+Bhc3yO77e22Z2zvm7dss+0jtvcWP2ubbRNAVYO8jP+ZpOsXWP7jiFhV/Dxbb1sA6tY37BHxvKQTI+gFQIOqfEB3t+1Xipf5S3utZHuj7Wnb091ut8JwAKoYNuw/lfRFSaskHZX0cK8VI2JrRHQiojM5OTnkcACqGirsEXEsIj6KiFlJ2yStqbctAHUbKuy2l897+A1J+3qtC2A89D3PbvtxSddIWmb7sKRNkq6xvUpSSDoo6bsN9og+Hn6457soLVrU7HVT9913X6PbR336hj0i1i+w+LEGegHQIC6XBZIg7EAShB1IgrADSRB2IAm+4noGOHLkSGn9iSeeaGzs22+/vbTOVZFnDo7sQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AE59nPAJ1Op7R+/Pjxobd93XXXldYfeeSRobeN8cKRHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeS4Dz7GWBmZqa0XuXfRff7V9CLFy8eetsYLxzZgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJzrOPgXvvvbe0Pjs729jYV155ZWPbxnjpe2S3vdL2723vt/2a7e8Vyy+w/Zztt4rbpc23C2BYg7yM/1DSDyLi3yVdJeku21dIul/S7oi4VNLu4jGAMdU37BFxNCJeLu6/J2m/pBWS1knaUay2Q9LNTTUJoLrT+oDO9pSkr0j6o6SLI+KoNPcHQdJFPZ6z0fa07elut1utWwBDGzjstpdI+q2k70fEXwd9XkRsjYhORHSYBBBoz0Bht/0pzQX9lxHxu2LxMdvLi/pySeVfzQLQqr6n3mxb0mOS9kfEj+aVdkraIGlLcft0Ix2eBapOudzvK6znnntuz9qmTZtKnzsxMVFax9ljkPPsV0v6tqRXbe8tlj2guZD/xvYdkv4i6dZmWgRQh75hj4g/SHKP8tfqbQdAU7hcFkiCsANJEHYgCcIOJEHYgST4iusInDx5srTe7zx8P1NTUz1r/f5VNPLgyA4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJ8H32EVixYkVp/cYbbyytP/PMM3W2g6Q4sgNJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoPMz75S0s8lfUbSrKStEfET25slfUdSt1j1gYh4tqlGz2RLliwprT/11FMj6gSZDXJRzYeSfhARL9v+tKSXbD9X1H4cEf/dXHsA6jLI/OxHJR0t7r9ne7+k8kvCAIyd03rPbntK0lck/bFYdLftV2xvt720x3M22p62Pd3tdhdaBcAIDBx220sk/VbS9yPir5J+KumLklZp7sj/8ELPi4itEdGJiM7k5GQNLQMYxkBht/0pzQX9lxHxO0mKiGMR8VFEzEraJmlNc20CqKpv2G1b0mOS9kfEj+YtXz5vtW9I2ld/ewDqMsin8VdL+rakV23vLZY9IGm97VWSQtJBSd9tpEMAtRjk0/g/SPICJc6pA2cQrqADkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4k4YgY3WB2V9Kf5y1aJun4yBo4PePa27j2JdHbsOrs7XMRseD/fxtp2D8xuD0dEZ3WGigxrr2Na18SvQ1rVL3xMh5IgrADSbQd9q0tj19mXHsb174kehvWSHpr9T07gNFp+8gOYEQIO5BEK2G3fb3tN2wfsH1/Gz30Yvug7Vdt77U93XIv223P2N43b9kFtp+z/VZxu+Acey31ttn2kWLf7bW9tqXeVtr+ve39tl+z/b1ieav7rqSvkey3kb9nt32OpDcl/aekw5L2SFofEa+PtJEebB+U1ImI1i/AsP1VSScl/TwivlQs+y9JJyJiS/GHcmlE3DcmvW2WdLLtabyL2YqWz59mXNLNkm5Ti/uupK9vagT7rY0j+xpJByLi7Yj4m6RfS1rXQh9jLyKel3TilMXrJO0o7u/Q3C/LyPXobSxExNGIeLm4/56kj6cZb3XflfQ1Em2EfYWkQ/MeH9Z4zfceknbZfsn2xrabWcDFEXFUmvvlkXRRy/2cqu803qN0yjTjY7Pvhpn+vKo2wr7QVFLjdP7v6ohYLekGSXcVL1cxmIGm8R6VBaYZHwvDTn9eVRthPyxp5bzHn5X0Tgt9LCgi3iluZyQ9qfGbivrYxzPoFrczLffzT+M0jfdC04xrDPZdm9OftxH2PZIutf1524slfUvSzhb6+ATbE8UHJ7I9IenrGr+pqHdK2lDc3yDp6RZ7+RfjMo13r2nG1fK+a33684gY+Y+ktZr7RP5Pkn7YRg89+vqCpP8vfl5ruzdJj2vuZd3fNfeK6A5JF0raLemt4vaCMertF5JelfSK5oK1vKXe/kNzbw1fkbS3+Fnb9r4r6Wsk+43LZYEkuIIOSIKwA0kQdiAJwg4kQdiBJAg7kARhB5L4Byx6uPlgmCU1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAN9ElEQVR4nO3df6hcdXrH8c/HGBXNBmNztTf+ylYEjYVml0EK/iBladBgjBEs+oeoUaOg4EKUBvuHEUWkdA0KZfFuFZO6dVV2gwZCu0EWo/+YjCE1SUMbldTNemPuJcgqKjbJ0z/ucbnGO2euM2fmTPK8X3CZmfPMOefJ5H7umZnvmfk6IgTgxHdS3Q0A6A/CDiRB2IEkCDuQBGEHkji5nzubO3duzJ8/v5+7BFLZt2+fxsfHPVWtq7DbvkbS05JmSPqXiHiy7P7z589Xs9nsZpcASjQajZa1jp/G254h6Z8lXStpgaRbbC/odHsAequb1+yXS3o/Ij6MiK8l/UrSsmraAlC1bsJ+rqTfT7q9v1j2LbZX2m7abo6NjXWxOwDd6CbsU70J8J1zbyNiJCIaEdEYGhrqYncAutFN2PdLOn/S7fMkfdxdOwB6pZuwb5N0se0f2j5F0s2SXq+mLQBV63joLSIO275f0n9oYujt+YjYXVlnACrV1Th7RGyStKmiXgD0EKfLAkkQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kERXs7hiMIyPj7esHT58uHTdrVu3ltaXLVtWWj/ppME9Xtxxxx0ta88++2zpujNmzKi6ndp1FXbb+yR9JumIpMMR0aiiKQDVq+LI/jcR0frQAmAgDO5zMACV6jbsIem3tt+1vXKqO9heabtpuzk2Ntbl7gB0qtuwXxERP5Z0raT7bF997B0iYiQiGhHRGBoa6nJ3ADrVVdgj4uPi8qCkDZIur6IpANXrOOy2z7D9g2+uS1osaVdVjQGoVjfvxp8jaYPtb7bzbxHx75V0lcyBAwdK6+vXry+tj4yMtKwdPXq0dN2PPvqotN5uHL34/x9IL7zwQsvanDlzStd9/PHHS+unnnpqJy3VquOwR8SHkv6qwl4A9BBDb0AShB1IgrADSRB2IAnCDiTBR1wHwOrVq0vrL774Yp86yWPt2rWl9Xvvvbe0ftFFF1XZTl9wZAeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJBhnHwBLly4trXczzj5v3rzS+oMPPlhab/cR2W6+Svqtt94qrW/YsKHjbeO7OLIDSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKMsw+A5cuXl9YPHTrU8bbbjYPPmjWr421365577imtX3rppaX1dl+DXWbFihWl9QsvvLDjbQ8qjuxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kATj7AOg3Vj47Nmz+9RJf23fvr20Pj4+3rN9X3DBBaX1k08+8aLR9shu+3nbB23vmrTsLNubbe8tLssnuwZQu+k8jX9B0jXHLFst6Y2IuFjSG8VtAAOsbdgjYoukY8/XXCZpXXF9naQbKu4LQMU6fYPunIgYlaTi8uxWd7S90nbTdnNsbKzD3QHoVs/fjY+IkYhoRERjaGio17sD0EKnYf/E9rAkFZcHq2sJQC90GvbXJd1WXL9N0mvVtAOgV9oOJtp+SdIiSXNt75f0iKQnJb1i+05JH0m6qZdN4vj19ttvt6w9/fTTpet+8cUXVbfzJw899FDPtj2o2oY9Im5pUfpJxb0A6CFOlwWSIOxAEoQdSIKwA0kQdiCJE+9zfKjUli1bSuurVq0qre/evbtl7euvv+6op+m66qqrWta6mWr6eJXvXwwkRdiBJAg7kARhB5Ig7EAShB1IgrADSTDOPgA+/fTT0vorr7xSWt+0aVOV7XzLxo0bS+u2e7bvM888s7S+fv360vqVV17ZsjZz5syOejqecWQHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQYZ++D0dHR0vqiRYtK6x988EGF3Rw/li5dWlpfsmRJnzo5MXBkB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkGGcfABHRVb2Xjh49Wlrv5fevt/u8+gMPPFBaX7hwYZXtHPfa/k/Zft72Qdu7Ji1bY/sPtncUP5zdAAy46fxZfkHSNVMsXxsRC4uf3n1VCoBKtA17RGyRdKgPvQDooW5ecN1v+73iaf6cVneyvdJ203ZzbGysi90B6EanYf+5pIskLZQ0Kulnre4YESMR0YiIxtDQUIe7A9CtjsIeEZ9ExJGIOCrpF5Iur7YtAFXrKOy2hyfdXC5pV6v7AhgMbcfZbb8kaZGkubb3S3pE0iLbCyWFpH2S7ulhj8e94eHh0vq2bdtK66+++mppffHixS1rp5xySum6vfbcc8+1rD3yyCN97ARtwx4Rt0yxuPX/IICBxOmyQBKEHUiCsANJEHYgCcIOJOF+fnyy0WhEs9ns2/5Qv6+++qplbdasWV1tu93vUsaPuDYaDTWbzSnn0ebIDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJ8FXS6Knt27fX3QIKHNmBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnG2afpyJEjLWs7d+4sXfeyyy4rrc+cObOjngbB5s2bS+s33XRTnzpBOxzZgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJxtkLe/fuLa2vWbOmZe3ll18uXffQoUOl9TrH2b/88svS+tatW0vrN998c2n9888//949feP0008vrZ922mkdbzujtkd22+fb/p3tPbZ3236gWH6W7c229xaXc3rfLoBOTedp/GFJqyLiUkl/Lek+2wskrZb0RkRcLOmN4jaAAdU27BExGhHbi+ufSdoj6VxJyyStK+62TtINvWoSQPe+1xt0tudL+pGkdySdExGj0sQfBElnt1hnpe2m7ebY2Fh33QLo2LTDbnuWpF9L+mlE/HG660XESEQ0IqIxNDTUSY8AKjCtsNueqYmg/zIiflMs/sT2cFEflnSwNy0CqELboTfblvScpD0R8dSk0uuSbpP0ZHH5Wk867JPbb7+9tP7OO+90vO21a9eW1mfPnt3xtru1cePG0vqbb75ZWp/49ejMjTfeWFpftWpVaf2SSy7peN8ZTWec/QpJt0raaXtHsexhTYT8Fdt3SvpIEh9cBgZY27BHxNuSWv35/km17QDoFU6XBZIg7EAShB1IgrADSRB2IAk+4toHjz32WN0t9My8efNK67feemvL2qOPPlq67skn8+tZJY7sQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEA5mFdl8H/cwzz7SsPfXUUy1rdVuwYEFpvd1n6RcvXlxav/vuu0vrw8PDpXX0D0d2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCcfbCeeedV1p/4oknWtauvvrq0nXvuuuu0vr4+HhpfcWKFaX166+/vmVt0aJFpevOmjWrtI4TB0d2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUhiOvOzny9pvaQ/l3RU0khEPG17jaS7JY0Vd304Ijb1qtG6lX2H+XXXXVe67oEDB6puB/jepnNSzWFJqyJiu+0fSHrX9uaitjYi/ql37QGoynTmZx+VNFpc/8z2Hknn9roxANX6Xq/Zbc+X9CNJ7xSL7rf9nu3nbc9psc5K203bzbGxsanuAqAPph1227Mk/VrSTyPij5J+LukiSQs1ceT/2VTrRcRIRDQiojE0NFRBywA6Ma2w256piaD/MiJ+I0kR8UlEHImIo5J+Ieny3rUJoFttw27bkp6TtCcinpq0fPLXhi6XtKv69gBUZTrvxl8h6VZJO23vKJY9LOkW2wslhaR9ku7pSYcAKjGdd+PfluQpSifsmDpwIuIMOiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKOiP7tzB6T9L+TFs2VVD5fcX0GtbdB7Uuit05V2duFETHl97/1Nezf2bndjIhGbQ2UGNTeBrUvid461a/eeBoPJEHYgSTqDvtIzfsvM6i9DWpfEr11qi+91fqaHUD/1H1kB9AnhB1Iopaw277G9n/bft/26jp6aMX2Pts7be+w3ay5l+dtH7S9a9Kys2xvtr23uJxyjr2aeltj+w/FY7fD9pKaejvf9u9s77G92/YDxfJaH7uSvvryuPX9NbvtGZL+R9LfStovaZukWyLiv/raSAu290lqRETtJ2DYvlrS55LWR8RfFsv+UdKhiHiy+EM5JyL+fkB6WyPp87qn8S5mKxqePM24pBsk3a4aH7uSvv5OfXjc6jiyXy7p/Yj4MCK+lvQrSctq6GPgRcQWSYeOWbxM0rri+jpN/LL0XYveBkJEjEbE9uL6Z5K+mWa81seupK++qCPs50r6/aTb+zVY872HpN/aftf2yrqbmcI5ETEqTfzySDq75n6O1XYa7346ZprxgXnsOpn+vFt1hH2qqaQGafzvioj4saRrJd1XPF3F9ExrGu9+mWKa8YHQ6fTn3aoj7PslnT/p9nmSPq6hjylFxMfF5UFJGzR4U1F/8s0MusXlwZr7+ZNBmsZ7qmnGNQCPXZ3Tn9cR9m2SLrb9Q9unSLpZ0us19PEdts8o3jiR7TMkLdbgTUX9uqTbiuu3SXqtxl6+ZVCm8W41zbhqfuxqn/48Ivr+I2mJJt6R/0DSP9TRQ4u+/kLSfxY/u+vuTdJLmnha93+aeEZ0p6Q/k/SGpL3F5VkD1Nu/Stop6T1NBGu4pt6u1MRLw/ck7Sh+ltT92JX01ZfHjdNlgSQ4gw5IgrADSRB2IAnCDiRB2IEkCDuQBGEHkvh/S1oWzvuBodgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAL7ElEQVR4nO3dX4hc9RnG8edx2xAwRWIz2mBCkxTBSqFRhliwiKVUohBiBUtzUVKQphcKBryoWFDxSkptqCDCWoOb0ipCK+YitJUgSG+qY9hq7NL6L62rIZmgRnNjNL692JN2jTtnJnPOzBl9vx8YZub8zux5GPbZM3POzP4cEQLw+XdO0wEAjAdlB5Kg7EASlB1IgrIDSXxhnBtbtWpVrFu3bpybBFI5dOiQjh075qXGKpXd9mZJv5Y0Jek3EXFv2frr1q1Tp9OpskkAJdrtds+xoV/G256S9ICkayVdKmmb7UuH/XkARqvKe/ZNkl6JiNci4qSkxyRtrScWgLpVKftFkt5YdH++WPYJtnfY7tjudLvdCpsDUEWVsi91EOBTn72NiOmIaEdEu9VqVdgcgCqqlH1e0tpF99dIeqtaHACjUqXsz0m62PZ628sk/VDS3npiAajb0KfeIuIj27dI+rMWTr3tjoiXaksGoFaVzrNHxD5J+2rKAmCE+LgskARlB5Kg7EASlB1IgrIDSVB2IAnKDiRB2YEkKDuQBGUHkqDsQBKUHUiCsgNJUHYgCcoOJEHZgSQoO5AEZQeSoOxAEpQdSIKyA0lQdiAJyg4kQdmBJCg7kARlB5Kg7EASlB1IgrIDSVSaxRWT79SpU6Xjr776aun4zp07S8f37WMS38+KSmW3fUjS+5JOSfooItp1hAJQvzr27N+JiGM1/BwAI8R7diCJqmUPSX+x/bztHUutYHuH7Y7tTrfbrbg5AMOqWvYrI+JySddKutn2VWeuEBHTEdGOiHar1aq4OQDDqlT2iHiruD4q6QlJm+oIBaB+Q5fd9rm2v3T6tqRrJB2sKxiAelU5Gn+hpCdsn/45v4+IP9WSCrX54IMPSscvueSS0vE1a9aUjp84caJ0fMWKFaXjGJ+hyx4Rr0n6Zo1ZAIwQp96AJCg7kARlB5Kg7EASlB1Igq+4otT8/Hzp+PHjx0vHOfU2OdizA0lQdiAJyg4kQdmBJCg7kARlB5Kg7EASnGdHqYhoOgJqwp4dSIKyA0lQdiAJyg4kQdmBJCg7kARlB5LgPDtKFf8qvKd+/6oak4M9O5AEZQeSoOxAEpQdSIKyA0lQdiAJyg4kwXl2VDI7O1s6vmHDhjElQT999+y2d9s+avvgomXn237K9svF9crRxgRQ1SAv4x+RtPmMZbdL2h8RF0vaX9wHMMH6lj0inpH09hmLt0qaKW7PSLq+5lwAajbsAboLI+KwJBXXF/Ra0fYO2x3bnW63O+TmAFQ18qPxETEdEe2IaLdarVFvDkAPw5b9iO3VklRcH60vEoBRGLbseyVtL25vl/RkPXEAjErf8+y2H5V0taRVtucl3SXpXkmP275J0n8k3TjKkBjeOeeU/z1fubL8rOk777xTOj43N3fWmdCMvmWPiG09hr5bcxYAI8THZYEkKDuQBGUHkqDsQBKUHUiCr7h+zi1fvrx0fMuWLaXje/bsqTMOGsSeHUiCsgNJUHYgCcoOJEHZgSQoO5AEZQeSoOxAEpQdSIKyA0lQdiAJyg4kQdmBJCg7kARlB5Lg++yo5NixY01HwIDYswNJUHYgCcoOJEHZgSQoO5AEZQeSoOxAEpxnRyUzMzOl47t27RpTEvTTd89ue7fto7YPLlp2t+03bc8Wl+tGGxNAVYO8jH9E0uYllu+KiI3FZV+9sQDUrW/ZI+IZSW+PIQuAEapygO4W2y8UL/NX9lrJ9g7bHdudbrdbYXMAqhi27A9K+pqkjZIOS7qv14oRMR0R7Yhot1qtITcHoKqhyh4RRyLiVER8LOkhSZvqjQWgbkOV3fbqRXe/L+lgr3UBTIa+59ltPyrpakmrbM9LukvS1bY3SgpJhyT9dIQZMUKbNy91ouX/mJ/986Nv2SNi2xKLHx5BFgAjxMdlgSQoO5AEZQeSoOxAEpQdSIKvuCa3fv36So8/efJk6fjx48d7jp133nmVto2zw54dSIKyA0lQdiAJyg4kQdmBJCg7kARlB5LgPHtyU1NTlR4fEaXjH374YaWfj/qwZweSoOxAEpQdSIKyA0lQdiAJyg4kQdmBJDjPnly73S4d37hxY+n47Oxs6fj999/fc+yee+4pfSzqxZ4dSIKyA0lQdiAJyg4kQdmBJCg7kARlB5LgPDtK3XDDDaXjr7/+eun4nXfeWWccVNB3z257re2nbc/Zfsn2rcXy820/Zfvl4nrl6OMCGNYgL+M/knRbRHxd0rck3Wz7Ukm3S9ofERdL2l/cBzCh+pY9Ig5HxIHi9vuS5iRdJGmrpJlitRlJ148qJIDqzuoAne11ki6T9DdJF0bEYWnhD4KkC3o8Zoftju1Ot9utlhbA0AYuu+0Vkv4gaWdEvDfo4yJiOiLaEdFutVrDZARQg4HKbvuLWij67yLij8XiI7ZXF+OrJR0dTUQAdRjkaLwlPSxpLiJ+tWhor6Ttxe3tkp6sPx4mne3Sy9TUVM8LxmuQ8+xXSvqRpBdtn/7y8h2S7pX0uO2bJP1H0o2jiQigDn3LHhF/leQew9+tNw6AUeHjskASlB1IgrIDSVB2IAnKDiTBV1xRybvvvls6/uyzz/Ycu+KKK+qOgxLs2YEkKDuQBGUHkqDsQBKUHUiCsgNJUHYgCc6zo9T09HTp+PLly0vHN2zYUGccVMCeHUiCsgNJUHYgCcoOJEHZgSQoO5AEZQeS4Dw7Sm3ZsqV0/MCBA6Xjy5YtqzMOKmDPDiRB2YEkKDuQBGUHkqDsQBKUHUiCsgNJ9D3PbnutpD2SviLpY0nTEfFr23dL+omkbrHqHRGxb1RB0YwHHnig6QioySAfqvlI0m0RccD2lyQ9b/upYmxXRPxydPEA1GWQ+dkPSzpc3H7f9pyki0YdDEC9zuo9u+11ki6T9Ldi0S22X7C92/bKHo/ZYbtju9PtdpdaBcAYDFx22ysk/UHSzoh4T9KDkr4maaMW9vz3LfW4iJiOiHZEtFutVg2RAQxjoLLb/qIWiv67iPijJEXEkYg4FREfS3pI0qbRxQRQVd+y27akhyXNRcSvFi1fvWi170s6WH88AHUZ5Gj8lZJ+JOlF27PFsjskbbO9UVJIOiTppyNJCKAWgxyN/6skLzHEOXXgM4RP0AFJUHYgCcoOJEHZgSQoO5AEZQeSoOxAEpQdSIKyA0lQdiAJyg4kQdmBJCg7kARlB5JwRIxvY3ZX0r8XLVol6djYApydSc02qbkksg2rzmxfjYgl///bWMv+qY3bnYhoNxagxKRmm9RcEtmGNa5svIwHkqDsQBJNl3264e2XmdRsk5pLItuwxpKt0ffsAMan6T07gDGh7EASjZTd9mbb/7T9iu3bm8jQi+1Dtl+0PWu703CW3baP2j64aNn5tp+y/XJxveQcew1lu9v2m8VzN2v7uoayrbX9tO052y/ZvrVY3uhzV5JrLM/b2N+z256S9C9J35M0L+k5Sdsi4h9jDdKD7UOS2hHR+AcwbF8l6YSkPRHxjWLZLyS9HRH3Fn8oV0bEzyYk292STjQ9jXcxW9HqxdOMS7pe0o/V4HNXkusHGsPz1sSefZOkVyLitYg4KekxSVsbyDHxIuIZSW+fsXirpJni9owWflnGrke2iRARhyPiQHH7fUmnpxlv9LkryTUWTZT9IklvLLo/r8ma7z0k/cX287Z3NB1mCRdGxGFp4ZdH0gUN5zlT32m8x+mMacYn5rkbZvrzqpoo+1JTSU3S+b8rI+JySddKurl4uYrBDDSN97gsMc34RBh2+vOqmij7vKS1i+6vkfRWAzmWFBFvFddHJT2hyZuK+sjpGXSL66MN5/mfSZrGe6lpxjUBz12T0583UfbnJF1se73tZZJ+KGlvAzk+xfa5xYET2T5X0jWavKmo90raXtzeLunJBrN8wqRM491rmnE1/Nw1Pv15RIz9Iuk6LRyRf1XSz5vI0CPXBkl/Ly4vNZ1N0qNaeFn3oRZeEd0k6cuS9kt6ubg+f4Ky/VbSi5Je0EKxVjeU7dtaeGv4gqTZ4nJd089dSa6xPG98XBZIgk/QAUlQdiAJyg4kQdmBJCg7kARlB5Kg7EAS/wXieZ4s0PqxJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAANwklEQVR4nO3df6xU9ZnH8c+jC4laiLBc9Wp1Lwua9LrJQnNDGjSNmyr+DvSPq+WPhk1I8A9MAPtHtU2sRo1kkZJVNzV0JbAbtDa2RIzErSJRUNJ4MVRwye51DS3gDXeIUSQRWeXZP+5xc4F7vjPMOTNnLs/7lUxm5jxz5jwZ+Nwzc75n5mvuLgDnvvOqbgBAexB2IAjCDgRB2IEgCDsQxF+1c2PTpk3znp6edm4SCGX//v06cuSIjVUrFHYzu0XSP0s6X9K/uvvK1ON7eno0MDBQZJMAEvr6+nJrTb+NN7PzJf2LpFsl9UpaaGa9zT4fgNYq8pl9jqQP3f0jdz8h6TeS5pfTFoCyFQn7FZIOjLp/MFt2CjNbYmYDZjZQq9UKbA5AEUXCPtZBgDPOvXX3te7e5+59XV1dBTYHoIgiYT8o6cpR978t6eNi7QBolSJhf1fS1WY23cwmSvqRpM3ltAWgbE0Pvbn7V2Z2r6T/0MjQ2zp3/6C0zgCUqtA4u7tvkbSlpF4AtBCnywJBEHYgCMIOBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQRB2IAjCDgRB2IEgCDsQRFunbAbGi/7+/mTd/YzJj07x4osvltlOKdizA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQjLMjpMceeyxZf+WVV5L1FStWlNlOWxQKu5ntl/S5pK8lfeXufWU0BaB8ZezZ/8Hdj5TwPABaiM/sQBBFw+6S/mBmu8xsyVgPMLMlZjZgZgO1Wq3g5gA0q2jYr3P370q6VdJSM/v+6Q9w97Xu3ufufV1dXQU3B6BZhcLu7h9n18OSNkmaU0ZTAMrXdNjN7CIzm/TNbUnzJO0tqzEA5SpyNP5SSZvM7Jvnec7dXy2lK6AEq1evzq3VG2efOHFisn777bc31VOVmg67u38k6e9L7AVACzH0BgRB2IEgCDsQBGEHgiDsQBB8xRXnrB07duTWTpw4kVz3zjvvTNbnzp3bVE9VYs8OBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0Ewzn6OGxwcTNYffPDBZH3dunXJ+gUXXHDWPZVl+/btyfo777yTW+vt7U2uu2bNmqZ66mTs2YEgCDsQBGEHgiDsQBCEHQiCsANBEHYgCMbZz3H9/f3J+p49e5L1Rx55JFmfOXPmWfdUlvvuuy9ZHx4ezq29/PLLyXUvv/zypnrqZOzZgSAIOxAEYQeCIOxAEIQdCIKwA0EQdiAIxtnPcZMnT07Wsym3c9X7ffVWOnToULJe77v6552Xvy/78ssvm+ppPKu7ZzezdWY2bGZ7Ry2bamavmdlgdj2ltW0CKKqRt/HrJd1y2rL7JW1196slbc3uA+hgdcPu7m9J+uS0xfMlbchub5C0oOS+AJSs2QN0l7r7kCRl15fkPdDMlpjZgJkN1Gq1JjcHoKiWH41397Xu3ufufV1dXa3eHIAczYb9sJl1S1J2nf/1IgAdodmwb5a0KLu9SNJL5bQDoFXqjrOb2fOSbpA0zcwOSvqFpJWSfmtmiyX9RVL6S9Noqaeeeiq3tnPnzuS6s2fPTtZ7enqaaakh9cbwH3/88WT92LFjyfrNN9+cWxuP86sXVTfs7r4wp/SDknsB0EKcLgsEQdiBIAg7EARhB4Ig7EAQfMV1HDh69GiyvnLlytzahAkTkutu3LgxWb/wwguT9SIefvjhZP2ZZ55J1q+66qpkfcuWLWfd07mMPTsQBGEHgiDsQBCEHQiCsANBEHYgCMIOBME4ewcYGhpK1m+88cZk/fDhw7m1emPZ11xzTbJeVGoc/4knnij03E8++WSh9aNhzw4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQTDOXoKTJ08m69u2bUvW582bV+j5U1MTv/nmm8l1L7vssmR90aJFyfrx48eT9fXr1+fW3D257ooVK5L1O+64I1nHqdizA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQjLOXYPv27cl6aupgSTKzZD01ji5J1157bW7tjTfeSK5br/7CCy8k64ODg8n6gQMHcmv1xvhXrVqVrOPs1N2zm9k6Mxs2s72jlj1kZofMbHd2ua21bQIoqpG38esl3TLG8jXuPiu7MPUG0OHqht3d35L0SRt6AdBCRQ7Q3Wtm72dv86fkPcjMlpjZgJkN1Gq1ApsDUESzYf+VpBmSZkkakrQ674Huvtbd+9y9r6urq8nNASiqqbC7+2F3/9rdT0r6taQ55bYFoGxNhd3Mukfd/aGkvXmPBdAZ6o6zm9nzkm6QNM3MDkr6haQbzGyWJJe0X9I9LeyxI7z99tu5tXq/6z5x4sRkferUqcn666+/nqxPmjQpt7Z8+fLkups2bUrW643D1/tOeuocgtTv3UvS9OnTk/Vdu3Yl6/Ve12jqht3dF46x+NkW9AKghThdFgiCsANBEHYgCMIOBEHYgSD4imuD1qxZk1ubOXNmct16UwvfdNNNTfXUiKeffjpZ/+KLL5L1V199tcx2TlFv2G7BggXJOkNrZ4c9OxAEYQeCIOxAEIQdCIKwA0EQdiAIwg4EwTh7g+6+++7cWr2fip48eXLZ7TTs6NGjyfrOnTsLPX+9n9GeMWNG08998cUXN70uzsSeHQiCsANBEHYgCMIOBEHYgSAIOxAEYQeCYJy9Qf39/VW3kOv48eO5tY0bNybX/fTTT5P13t7eZH3u3LnJOjoHe3YgCMIOBEHYgSAIOxAEYQeCIOxAEIQdCIJx9nPAc889l1t79NFHk+t2d3cn66mpqjG+1N2zm9mVZrbNzPaZ2QdmtixbPtXMXjOzwex6SuvbBdCsRt7GfyXpJ+7+HUnfk7TUzHol3S9pq7tfLWlrdh9Ah6obdncfcvf3stufS9on6QpJ8yVtyB62QVJ6rh4AlTqrA3Rm1iNptqQ/SrrU3YekkT8Iki7JWWeJmQ2Y2UCtVivWLYCmNRx2M/uWpN9JWu7u6V8xHMXd17p7n7v3dXV1NdMjgBI0FHYzm6CRoG90999niw+bWXdW75Y03JoWAZSh7tCbmZmkZyXtc/dfjiptlrRI0srs+qWWdAh99tlnyfqqVatyayP/fPkeeOCBZL3Kn8FGuRoZZ79O0o8l7TGz3dmyn2kk5L81s8WS/iKpc7/wDaB+2N19h6S83cMPym0HQKtwuiwQBGEHgiDsQBCEHQiCsANB8BXXceD6669P1gcHB3Nry5YtS667dOnSpnrC+MOeHQiCsANBEHYgCMIOBEHYgSAIOxAEYQeCYJx9HFi+fHmyfs899+TW7rrrrrLbwTjFnh0IgrADQRB2IAjCDgRB2IEgCDsQBGEHgmCcfRxYvHhxoTogsWcHwiDsQBCEHQiCsANBEHYgCMIOBEHYgSDqht3MrjSzbWa2z8w+MLNl2fKHzOyQme3OLre1vl0AzWrkpJqvJP3E3d8zs0mSdpnZa1ltjbs/0br2AJSlkfnZhyQNZbc/N7N9kq5odWMAynVWn9nNrEfSbEl/zBbda2bvm9k6M5uSs84SMxsws4FarVaoWQDNazjsZvYtSb+TtNzdj0r6laQZkmZpZM+/eqz13H2tu/e5e19XV1cJLQNoRkNhN7MJGgn6Rnf/vSS5+2F3/9rdT0r6taQ5rWsTQFGNHI03Sc9K2ufuvxy1vHvUw34oaW/57QEoSyNH46+T9GNJe8xsd7bsZ5IWmtksSS5pv6T83zMGULlGjsbvkGRjlLaU3w6AVuEMOiAIwg4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQRB2IAjCDgRB2IEgCDsQBGEHgiDsQBDm7u3bmFlN0p9HLZom6UjbGjg7ndpbp/Yl0Vuzyuztb9x9zN9/a2vYz9i42YC791XWQEKn9tapfUn01qx29cbbeCAIwg4EUXXY11a8/ZRO7a1T+5LorVlt6a3Sz+wA2qfqPTuANiHsQBCVhN3MbjGz/zKzD83s/ip6yGNm+81sTzYN9UDFvawzs2Ez2ztq2VQze83MBrPrMefYq6i3jpjGOzHNeKWvXdXTn7f9M7uZnS/pvyXdJOmgpHclLXT3/2xrIznMbL+kPnev/AQMM/u+pGOS/s3d/y5b9k+SPnH3ldkfyinu/tMO6e0hSceqnsY7m62oe/Q045IWSPpHVfjaJfq6S2143arYs8+R9KG7f+TuJyT9RtL8CvroeO7+lqRPTls8X9KG7PYGjfxnabuc3jqCuw+5+3vZ7c8lfTPNeKWvXaKvtqgi7FdIOjDq/kF11nzvLukPZrbLzJZU3cwYLnX3IWnkP4+kSyru53R1p/Fup9OmGe+Y166Z6c+LqiLsY00l1Unjf9e5+3cl3SppafZ2FY1paBrvdhljmvGO0Oz050VVEfaDkq4cdf/bkj6uoI8xufvH2fWwpE3qvKmoD38zg252PVxxP/+vk6bxHmuacXXAa1fl9OdVhP1dSVeb2XQzmyjpR5I2V9DHGczsouzAiczsIknz1HlTUW+WtCi7vUjSSxX2copOmcY7b5pxVfzaVT79ubu3/SLpNo0ckf8fST+vooecvv5W0p+yywdV9ybpeY28rftfjbwjWizpryVtlTSYXU/toN7+XdIeSe9rJFjdFfV2vUY+Gr4vaXd2ua3q1y7RV1teN06XBYLgDDogCMIOBEHYgSAIOxAEYQeCIOxAEIQdCOL/ADloHJ7szM11AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAOhElEQVR4nO3db4xUZZbH8d9RMRqY8GdplNDt9uwEExGVgQquuhjNxImiBglhGV+MbGLsSYRkxswLDPsCX0hi1mUIMYaEER1GRwzCEEg07ighEY1BS8IqbrvqEmQYCTSS4L8EFM6+6Mukwaqnmrq36pac7yfpVPU9des5qfSvb1U9t+oxdxeA898FZTcAoD0IOxAEYQeCIOxAEIQdCOKidg42fvx47+3tbeeQQCj79u3TkSNHrFYtV9jN7HZJqyRdKOkpd38sdfve3l5Vq9U8QwJIqFQqdWtNP403swslPSnpDklTJN1rZlOavT8ArZXnNftMSZ+4+153PyHpBUlzimkLQNHyhH2SpL8O+f1Atu0MZtZnZlUzqw4MDOQYDkAeecJe602A75176+5r3L3i7pWurq4cwwHII0/YD0jqGfJ7t6TP8rUDoFXyhP0dSZPN7MdmdrGkX0jaWkxbAIrW9NSbu39nZosl/ZcGp96edvcPCusMQKFyzbO7+8uSXi6oFwAtxOmyQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQRB2IAjCDgRB2IEgCDsQBGEHgiDsQBCEHQiCsANBEHYgCMIOBJFrFddOcuLEiWR9+/btyfqll16arL/55pt1a8eOHUvu+8QTTyTrc+fOTda7u7uT9VaaNGlSsj5nzpxk/YorriiyHeSQK+xmtk/Sl5JOSvrO3StFNAWgeEUc2W919yMF3A+AFuI1OxBE3rC7pL+Y2btm1lfrBmbWZ2ZVM6sODAzkHA5As/KG/SZ3ny7pDkmLzOzms2/g7mvcveLula6urpzDAWhWrrC7+2fZ5WFJmyXNLKIpAMVrOuxmNtLMfnT6uqSfS9pTVGMAipXn3fjLJG02s9P387y7v1JIV01YtWpVsr5kyZI2dXLunn/++bJbaNpDDz2UrFcq9WdjH3jggeS+8+bNS9bHjBmTrONMTYfd3fdKuq7AXgC0EFNvQBCEHQiCsANBEHYgCMIOBHHefMT1mWeeKW3sCRMmJOuzZs1qUyffd9VVVyXr/f39yfrhw4eT9R07diTrb7/9dlM1SZoxY0ayPm3atGQdZ+LIDgRB2IEgCDsQBGEHgiDsQBCEHQiCsANBnDfz7G+88Uayvn///mQ9z1ceX3zxxcn6qFGjmr7vsh0/fjxZv/rqq5P1vXv3Nj32iy++mKwzz35uOLIDQRB2IAjCDgRB2IEgCDsQBGEHgiDsQBDnzTz7uHHjctVR286dO5P1PPPol1xySbLe11dzRTE0iSM7EARhB4Ig7EAQhB0IgrADQRB2IAjCDgRx3syzo7aTJ08m68uWLUvWV65cWWQ7Z/joo4+S9e7u7paNHVHDI7uZPW1mh81sz5Bt48zsVTP7OLsc29o2AeQ1nKfxf5B0+1nbHpa0zd0nS9qW/Q6ggzUMu7u/LunoWZvnSFqXXV8n6Z6C+wJQsGbfoLvM3Q9KUnZZd7EzM+szs6qZVQcGBpocDkBeLX833t3XuHvF3StdXV2tHg5AHc2G/ZCZTZSk7DK91CeA0jUb9q2SFmbXF0raUkw7AFql4Ty7ma2XdIuk8WZ2QNIySY9J2mBm90vaL2l+K5tE2ocffli3tnbt2uS+K1asyDX2iBEjkvVNmzbVrV1++eW5xsa5aRh2d7+3TulnBfcCoIU4XRYIgrADQRB2IAjCDgRB2IEg+IjrD8Cnn36arF9zzTV1a40+4prXBRekjxepj6maWdHtIIEjOxAEYQeCIOxAEIQdCIKwA0EQdiAIwg4EwTz7D8ALL7yQrLd6Lj3l+PHjyfr06dPr1m699dbkvgsWLEjW77777mR94sSJyXo0HNmBIAg7EARhB4Ig7EAQhB0IgrADQRB2IAjm2X8A5s9Pf1P3nj176tZee+215L6HDh1qqqcibN++PVf9wQcfTNaXL19et7Z48eLkviNHjkzWf4g4sgNBEHYgCMIOBEHYgSAIOxAEYQeCIOxAEObubRusUql4tVpt23iQjh07lqx/8cUXyfrRo0eT9fXr1yfrjz/+eN1aO//2zjZ37txkfePGjcl6p37nfaVSUbVardlcwyO7mT1tZofNbM+QbY+Y2d/MbHf2M7vIhgEUbzhP4/8g6fYa21e6+7Ts5+Vi2wJQtIZhd/fXJaWfywHoeHneoFtsZu9lT/PH1ruRmfWZWdXMqgMDAzmGA5BHs2FfLeknkqZJOihpRb0buvsad6+4e6Wrq6vJ4QDk1VTY3f2Qu59091OSfi9pZrFtAShaU2E3s6Hf0TtXUv3PWALoCA0/z25m6yXdImm8mR2QtEzSLWY2TZJL2ifpVy3sETmMHj06V72npydZv+6665L12bPrz8o++uijyX0bfRY/j82bNyfrzz77bLJ+3333FdlOWzQMu7vfW2Pz2hb0AqCFOF0WCIKwA0EQdiAIwg4EQdiBIPgqabTUzTffXLf2yiuvJPedN29esr5ly5amehqO/v7+lt13WTiyA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQzLOjNBdckD7WXH/99cl6K+fZp06d2rL7LgtHdiAIwg4EQdiBIAg7EARhB4Ig7EAQhB0Ignn2Nvjqq6+S9eeeey5Zv/baa5P1G2+88Zx76gSnTp1K1nft2tWysS+6KP2nP3Pm+bfuCUd2IAjCDgRB2IEgCDsQBGEHgiDsQBCEHQiCefYCNJpHv+2225L1nTt3JuvffPPNOffUKb7++uu6tdWrVyf33bhxY9Ht/N2MGTOS9cmTJ7ds7LI0PLKbWY+ZbTezfjP7wMx+nW0fZ2avmtnH2eXY1rcLoFnDeRr/naTfuvtVkv5Z0iIzmyLpYUnb3H2ypG3Z7wA6VMOwu/tBd9+VXf9SUr+kSZLmSFqX3WydpHta1SSA/M7pDToz65X0U0k7JV3m7gelwX8IkibU2afPzKpmVh0YGMjXLYCmDTvsZjZK0iZJv3H3L4a7n7uvcfeKu1e6urqa6RFAAYYVdjMbocGg/8nd/5xtPmRmE7P6REmHW9MigCI0nHozM5O0VlK/u/9uSGmrpIWSHssuW/e9vh1uyZIlyXqjqbVGPv/882R9woSar6AkSSNGjMg19rfffpusP/XUU8n60qVL69aOHTvWVE+nuXuyPnr06Lq1devW1a2dr4Yzz36TpF9Ket/Mdmfblmow5BvM7H5J+yXNb02LAIrQMOzu/oYkq1P+WbHtAGgVTpcFgiDsQBCEHQiCsANBEHYgCD7iWoC77rorWW/0Uc5Genp6kvVZs2bVreU9a7HRKc47duzIdf95pObRJemtt96qW7vyyiuLbqfjcWQHgiDsQBCEHQiCsANBEHYgCMIOBEHYgSCYZy/ADTfckKwvWrQoWX/yySdzjV/mXHcejZZNXr58ebI+f376U9W9vb3n2tJ5jSM7EARhB4Ig7EAQhB0IgrADQRB2IAjCDgTBPHsBxowZk6yvXLkyWV+wYEGy/tJLLyXrU6dOrVvbsGFDct9GpkyZkmv/O++8s26t0Tx4d3d3rrFxJo7sQBCEHQiCsANBEHYgCMIOBEHYgSAIOxCENVrj2sx6JP1R0uWSTkla4+6rzOwRSQ9IOv3F4kvd/eXUfVUqFa9Wq7mbBlBbpVJRtVqtuerycE6q+U7Sb919l5n9SNK7ZvZqVlvp7v9ZVKMAWmc467MflHQwu/6lmfVLmtTqxgAU65xes5tZr6SfStqZbVpsZu+Z2dNmNrbOPn1mVjWzaqOlhAC0zrDDbmajJG2S9Bt3/0LSakk/kTRNg0f+FbX2c/c17l5x90redccANG9YYTezERoM+p/c/c+S5O6H3P2ku5+S9HtJM1vXJoC8GobdzEzSWkn97v67IdsnDrnZXEl7im8PQFGG8278TZJ+Kel9M9udbVsq6V4zmybJJe2T9KuWdAigEMN5N/4NSbXm7ZJz6gA6C2fQAUEQdiAIwg4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQRB2IAjCDgRB2IEgCDsQBGEHgmj4VdKFDmY2IOnTIZvGSzrStgbOTaf21ql9SfTWrCJ7+0d3r/n9b20N+/cGN6u6e6W0BhI6tbdO7Uuit2a1qzeexgNBEHYgiLLDvqbk8VM6tbdO7Uuit2a1pbdSX7MDaJ+yj+wA2oSwA0GUEnYzu93M/tfMPjGzh8vooR4z22dm75vZbjMrdX3pbA29w2a2Z8i2cWb2qpl9nF3WXGOvpN4eMbO/ZY/dbjObXVJvPWa23cz6zewDM/t1tr3Uxy7RV1set7a/ZjezCyV9JOk2SQckvSPpXnf/n7Y2UoeZ7ZNUcffST8Aws5slfSXpj+4+Ndv2H5KOuvtj2T/Kse6+pEN6e0TSV2Uv452tVjRx6DLjku6R9G8q8bFL9PWvasPjVsaRfaakT9x9r7ufkPSCpDkl9NHx3P11SUfP2jxH0rrs+joN/rG0XZ3eOoK7H3T3Xdn1LyWdXma81Mcu0VdblBH2SZL+OuT3A+qs9d5d0l/M7F0z6yu7mRouc/eD0uAfj6QJJfdztobLeLfTWcuMd8xj18zy53mVEfZaS0l10vzfTe4+XdIdkhZlT1cxPMNaxrtdaiwz3hGaXf48rzLCfkBSz5DfuyV9VkIfNbn7Z9nlYUmb1XlLUR86vYJudnm45H7+rpOW8a61zLg64LErc/nzMsL+jqTJZvZjM7tY0i8kbS2hj+8xs5HZGycys5GSfq7OW4p6q6SF2fWFkraU2MsZOmUZ73rLjKvkx6705c/dve0/kmZr8B35/5P072X0UKevf5L039nPB2X3Jmm9Bp/WfavBZ0T3S/oHSdskfZxdjuug3p6V9L6k9zQYrIkl9fYvGnxp+J6k3dnP7LIfu0RfbXncOF0WCIIz6IAgCDsQBGEHgiDsQBCEHQiCsANBEHYgiP8HOXJIQ1I98wYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAN9ElEQVR4nO3dX4hcZZrH8d9PnSEaR6KbNgaVzazmYmU1Kh0RsowukiF600pUDDFmQTaiCTgYUHEheuFFEJ1BYRXjv8nIrOPIKOZCZyfoiIzIaCdEEze4GukdM2k7LREmQdA1efaiK0NP7HqrrXPqT/J8P1BU1XnqrfNQ9K9PVb1V9ToiBODYd1yvGwDQHYQdSIKwA0kQdiAJwg4kcUI3dzZ79uyYN29eN3cJpDIyMqLPP//cU9Uqhd32EkkPSzpe0pMRsb50+3nz5ml4eLjKLgEUDA4ONq21/TTe9vGS/kPSlZLOk7TM9nnt3h+Azqrymv0SSR9HxCcR8bWkX0kaqqctAHWrEvYzJX066fruxra/YXuV7WHbw+Pj4xV2B6CKKmGf6k2Ab332NiI2RMRgRAwODAxU2B2AKqqEfbeksyddP0vSnmrtAOiUKmF/V9J82z+0/X1JN0jaVE9bAOrW9tRbRHxje42k/9LE1NvTEfFBbZ0BqFWlefaIeEXSKzX1AqCD+LgskARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSVRaxRXopa+//rpYX7JkSdParl27imPfe++9Yn3WrFnFej+qFHbbI5L2Szoo6ZuIGKyjKQD1q+PI/i8R8XkN9wOgg3jNDiRRNewh6Xe2t9heNdUNbK+yPWx7eHx8vOLuALSratgXRcTFkq6UtNr2j468QURsiIjBiBgcGBiouDsA7aoU9ojY0zjfK+klSZfU0RSA+rUddtszbf/g8GVJP5a0o67GANSryrvxcyS9ZPvw/fxnRPy2lq5w1Ni/f3+lesnMmTOL9S1bthTrb7zxRtPaggULimNPPPHEYv1o1HbYI+ITSeVHDEDfYOoNSIKwA0kQdiAJwg4kQdiBJPiK6zFgdHS0ae2RRx4pjh0ZGam071bTX62+Slry0EMPFevbt28v1iOiaW3+/PnFsYcOHSrWj0Yc2YEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCebZjwFvvfVW09oDDzzQ0X3PmDGjWL/99tub1l588cXi2LVr17bV02GNr19PafXq1cWxx+JXXDmyA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EASzLMfBR599NFi/c4772z7vu+4445ifc6cOcX6bbfdVqyfdNJJTWut5tEXLlxYrI+NjRXrZ5xxRtPaokWLimOPRRzZgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJ5tmPAgcOHCjWv/zyy6a1c889tzj23nvvLdZPPvnkYr2Vffv2Na3df//9xbGfffZZsd5qSefHHnusae2EE/L96bc8stt+2vZe2zsmbTvN9mbbHzXOT+1smwCqms7T+J9LWnLEtrslvRYR8yW91rgOoI+1DHtEvCnpyOdiQ5I2Ni5vlHR1zX0BqFm7b9DNiYhRSWqcn97shrZX2R62PTw+Pt7m7gBU1fF34yNiQ0QMRsTgwMBAp3cHoIl2wz5me64kNc731tcSgE5oN+ybJK1sXF4p6eV62gHQKS0nG20/J+lySbNt75Z0r6T1kn5t+2ZJf5J0XSebzO76668v1l944YWmta1btxbHrlu3rlhfv359sf7VV18V66Xvyz/77LPFsa1e9j388MPF+tDQULGeTcuwR8SyJqUrau4FQAfxcVkgCcIOJEHYgSQIO5AEYQeSyPc9v6PQWWedVaxfcUXziZFWU2+tlk1etqzZZMyE5cuXF+u7du0q1kta/YT20qVL277vjDiyA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EASzLMfBVr97PGsWbPavu9PP/20WL/00kuL9Ygo1m03rbVaanrx4sXFOr4bjuxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kATz7MeAVssy99KNN97YtLZ27dri2FNOOaXudlLjyA4kQdiBJAg7kARhB5Ig7EAShB1IgrADSTDPfhQ4dOhQsb558+amtVbfN69qxYoVxfrGjRs7un9MX8sju+2nbe+1vWPStvts/9n2tsbpqs62CaCq6TyN/7mkJVNs/1lEXNg4vVJvWwDq1jLsEfGmpH1d6AVAB1V5g26N7fcbT/NPbXYj26tsD9seHh8fr7A7AFW0G/bHJJ0j6UJJo5IeanbDiNgQEYMRMTgwMNDm7gBU1VbYI2IsIg5GxCFJT0i6pN62ANStrbDbnjvp6jWSdjS7LYD+0HKe3fZzki6XNNv2bkn3Srrc9oWSQtKIpFs62GN6t956a7H+5JNPNq2Vfre9Dp2+f9SnZdgjYtkUm5/qQC8AOoiPywJJEHYgCcIOJEHYgSQIO5AEX3Htgv379xfrzz//fLH+xBNPFOul6a/LLrusOHbhwoXF+oMPPlis79mzp1hH/+DIDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJMM/eBVu2bCnWb7ml2jeES/Pwy5cvL459++23i/VW8+wLFiwo1tE/OLIDSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBLMs9fgww8/LNaXLl1a6f5bzdOff/75TWsHDhwojl29enVbPR12zjnnVBqP7uHIDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJMM9eg1dffbVY/+KLL4r1a665pli/6KKLivWDBw82rb3++uvFsfv27SvWI6JYnzt3brGO/tHyyG77bNu/t73T9ge2b29sP832ZtsfNc5P7Xy7ANo1nafx30haGxH/KOlSSattnyfpbkmvRcR8Sa81rgPoUy3DHhGjEbG1cXm/pJ2SzpQ0JGlj42YbJV3dqSYBVPed3qCzPU/SRZL+KGlORIxKE/8QJJ3eZMwq28O2h8fHx6t1C6Bt0w677ZMl/UbSTyLiL9MdFxEbImIwIgYHBgba6RFADaYVdtvf00TQfxkRLzY2j9me26jPlbS3My0CqEPLqTdPrAf8lKSdEfHTSaVNklZKWt84f7kjHR4Fjjuu/D+ztKTydOqlqTVJeuedd5rWrrvuuuLY2bNnF+t33XVXsT40NFSso39MZ559kaQVkrbb3tbYdo8mQv5r2zdL+pOk8l8VgJ5qGfaI+IOkZoeeK+ptB0Cn8HFZIAnCDiRB2IEkCDuQBGEHkuArrjUYGxurNP7006f8pPFfXXvttcX6pk2b2t53q6/nXnzxxW3fN/oLR3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIJ59hpccMEFlcY//vjjxXqrn3Mu/QLQunXrimNLyz3j2MKRHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSYJ69Bq1+O/2ZZ54p1tesWVOsL168uFgv/Tb8DTfcUByLPDiyA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAS01mf/WxJv5B0hqRDkjZExMO275P0b5LGGze9JyJe6VSj/WzGjBnF+k033VSpDtRhOh+q+UbS2ojYavsHkrbY3tyo/SwiHuxcewDqMp312UcljTYu77e9U9KZnW4MQL2+02t22/MkXSTpj41Na2y/b/tp26c2GbPK9rDt4fHx8aluAqALph122ydL+o2kn0TEXyQ9JukcSRdq4sj/0FTjImJDRAxGxGDpt9IAdNa0wm77e5oI+i8j4kVJioixiDgYEYckPSHpks61CaCqlmG3bUlPSdoZET+dtH3upJtdI2lH/e0BqMt03o1fJGmFpO22tzW23SNpme0LJYWkEUm3dKRDALWYzrvxf5DkKUop59SBoxWfoAOSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiThiOjezuxxSf87adNsSZ93rYHvpl9769e+JHprV529/X1ETPn7b10N+7d2bg9HxGDPGijo1976tS+J3trVrd54Gg8kQdiBJHod9g093n9Jv/bWr31J9NaurvTW09fsALqn10d2AF1C2IEkehJ220tsf2j7Y9t396KHZmyP2N5ue5vt4R738rTtvbZ3TNp2mu3Ntj9qnE+5xl6PervP9p8bj90221f1qLezbf/e9k7bH9i+vbG9p49doa+uPG5df81u+3hJ/yNpsaTdkt6VtCwi/rurjTRhe0TSYET0/AMYtn8k6YCkX0TEPzW2PSBpX0Ssb/yjPDUi7uqT3u6TdKDXy3g3ViuaO3mZcUlXS/pX9fCxK/R1vbrwuPXiyH6JpI8j4pOI+FrSryQN9aCPvhcRb0rad8TmIUkbG5c3auKPpeua9NYXImI0IrY2Lu+XdHiZ8Z4+doW+uqIXYT9T0qeTru9Wf633HpJ+Z3uL7VW9bmYKcyJiVJr445F0eo/7OVLLZby76YhlxvvmsWtn+fOqehH2qZaS6qf5v0URcbGkKyWtbjxdxfRMaxnvbplimfG+0O7y51X1Iuy7JZ096fpZkvb0oI8pRcSexvleSS+p/5aiHju8gm7jfG+P+/mrflrGe6plxtUHj10vlz/vRdjflTTf9g9tf1/SDZI29aCPb7E9s/HGiWzPlPRj9d9S1JskrWxcXinp5R728jf6ZRnvZsuMq8ePXc+XP4+Irp8kXaWJd+R3Sfr3XvTQpK9/kPRe4/RBr3uT9Jwmntb9nyaeEd0s6e8kvSbpo8b5aX3U27OStkt6XxPBmtuj3v5ZEy8N35e0rXG6qtePXaGvrjxufFwWSIJP0AFJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEv8PAtwpxrzodd0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAANR0lEQVR4nO3df6jVdZ7H8ddrzTFSM11vdWsuOTv0x8bSql0kMoaWYe0Hkdc/ZtM/BhcKJylQGHDDIKO/JHaUgmXIKRt3MYeBmcogckIEGYihW7ilK2UrrqP542qCTRCz13nvH/frcrV7vvd6vt/zw/t+PuBwzvm+z/d83xzvy+8538/3nI8jQgAmv7/qdAMA2oOwA0kQdiAJwg4kQdiBJK5p58bmzp0b8+bNa+cmgVSOHDmiM2fOeKxapbDbfkDSi5KmSHolIjaWPX7evHkaHBysskkAJfr7+xvWmn4bb3uKpH+T9KCkOyStsH1Hs88HoLWqfGZfJOnziDgcEX+W9CtJS+tpC0DdqoT9Vkl/HHX/WLHsErZX2R60PTg0NFRhcwCqqBL2sQ4CfOvc24jYEhH9EdHf09NTYXMAqqgS9mOS+kbd/66kL6q1A6BVqoT9A0m32/6e7e9IWi5pZz1tAahb00NvETFs+ylJuzQy9LY1Ig7U1hmAWlUaZ4+IdyS9U1MvAFqI02WBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiirT8ljXzOnTvXsPbMM8+UrvvKK6+U1o8fP15a55eRLsWeHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSYJwdlRw+fLi0vnDhwoa13t7e0nU3bNhQWp85c2ZpHZdizw4kQdiBJAg7kARhB5Ig7EAShB1IgrADSTDOjlL79u0rrd9zzz2l9bLvrK9bt6503alTp5bWcWUqhd32EUlfSbogaTgi+utoCkD96tiz/0NEnKnheQC0EJ/ZgSSqhj0k/c72h7ZXjfUA26tsD9oeHBoaqrg5AM2qGvbFEbFQ0oOSnrT9g8sfEBFbIqI/Ivr5AUCgcyqFPSK+KK5PS3pD0qI6mgJQv6bDbnu67ZkXb0taIml/XY0BqFeVo/E3SXrD9sXneT0i3q2lK7TN2bNnS+v33ntvaX1gYKC0vn79+oa14m8HbdJ02CPisKS/r7EXAC3E0BuQBGEHkiDsQBKEHUiCsANJ8BXXSW54eLi0/vjjj5fWx/sK62uvvVZaZ3ite7BnB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkGGef5DZt2lRa37VrV2n95MmTpfVp06ZdcU/oDPbsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AE4+yTQNl31l966aXSdZctW1Zav/7665vqCd2HPTuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJME4+ySwdevWhrXz58+Xrrtx48a620GXGnfPbnur7dO2949aNsf2e7YPFdezW9smgKom8jb+l5IeuGzZ05J2R8TtknYX9wF0sXHDHhF7JX152eKlkrYVt7dJGqi5LwA1a/YA3U0RcUKSiusbGz3Q9irbg7YHh4aGmtwcgKpafjQ+IrZERH9E9Pf09LR6cwAaaDbsp2z3SlJxfbq+lgC0QrNh3ylpZXF7paS36mkHQKuMO85ue4ek+yTNtX1M0gZJGyX92vZjko5K+lErm0S5t99+u2Ht4YcfLl23r6+v7nbQpcYNe0SsaFD6Yc29AGghTpcFkiDsQBKEHUiCsANJEHYgCb7iehU4dOhQaf3dd99tWDt+/Hjd7VyRTz/9tGFt1qxZpevefPPNdbeTGnt2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCcfarwLZt20rrd911V8PanDlzKm179+7dpfXly5eX1s+ePduwdu2115au+/rrr5fWBwb46cMrwZ4dSIKwA0kQdiAJwg4kQdiBJAg7kARhB5JgnP0q8MILL5TW9+7d27B2zTXl/8QXLlwora9evbq0XvYz1pK0YMGChrXBwcHSdZcsWVJaf//990vrd955Z2k9G/bsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AE4+xd4OTJk6X14eHh0vqUKVOa3vbRo0dL648++mhp/e67725624sXLy6tr1u3rrT+7LPPltbffPPNK+5pMht3z257q+3TtvePWvac7eO29xWXh1rbJoCqJvI2/peSHhhj+eaImF9c3qm3LQB1GzfsEbFX0pdt6AVAC1U5QPeU7Y+Lt/mzGz3I9irbg7YHh4aGKmwOQBXNhv3nkr4vab6kE5J+1uiBEbElIvojor+np6fJzQGoqqmwR8SpiLgQEX+R9AtJi+ptC0Ddmgq77d5Rd5dJ2t/osQC6w7jj7LZ3SLpP0lzbxyRtkHSf7fmSQtIRST9pYY+T3rlz5yqtf9tttzW97i233FJaX79+fdPPXdUTTzxRWh+vd1xq3LBHxIoxFr/agl4AtBCnywJJEHYgCcIOJEHYgSQIO5AEX3GdBGbNmtX0utOmTauxk3rdcMMNnW5hUmHPDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJMM7eBSKiUn2yOnDgQGm9yvkFGbFnB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkGGfvArYr1a9WFy5cKK1v3ry5tL5mzZo625n02LMDSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKMs3eBvr6+SvVdu3Y1rD3yyCNN9VSXsrH08aaDPnjwYGn95ZdfbqqnrMbds9vus73H9kHbB2yvKZbPsf2e7UPF9ezWtwugWRN5Gz8s6acR8beS7pb0pO07JD0taXdE3C5pd3EfQJcaN+wRcSIiPipufyXpoKRbJS2VtK142DZJA61qEkB1V3SAzvY8SQsk/UHSTRFxQhr5D0HSjQ3WWWV70Pbg0NBQtW4BNG3CYbc9Q9JvJK2NiPMTXS8itkREf0T09/T0NNMjgBpMKOy2p2ok6Nsj4rfF4lO2e4t6r6TTrWkRQB3GHXrzyPcrX5V0MCI2jSrtlLRS0sbi+q2WdJjAjBkzSusvvvhiaX3FihUNazt27Chd9/777y+tj/fR67PPPiutr127tmFt+vTppevu2bOntH7dddeV1nGpiYyzL5b0Y0mf2N5XLFuvkZD/2vZjko5K+lFrWgRQh3HDHhG/l9To1xN+WG87AFqF02WBJAg7kARhB5Ig7EAShB1Igq+4XgUGBsq/drB9+/aGtbIxeEn65ptvmurpovGmTX7++ecb1lavXl267pQpU5rqCWNjzw4kQdiBJAg7kARhB5Ig7EAShB1IgrADSTDOPgmUjcN//fXXbewE3Yw9O5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiQxbtht99neY/ug7QO21xTLn7N93Pa+4vJQ69sF0KyJ/HjFsKSfRsRHtmdK+tD2e0Vtc0T8a+vaA1CXiczPfkLSieL2V7YPSrq11Y0BqNcVfWa3PU/SAkl/KBY9Zftj21ttz26wzirbg7YHh4aGKjULoHkTDrvtGZJ+I2ltRJyX9HNJ35c0XyN7/p+NtV5EbImI/ojo7+npqaFlAM2YUNhtT9VI0LdHxG8lKSJORcSFiPiLpF9IWtS6NgFUNZGj8Zb0qqSDEbFp1PLeUQ9bJml//e0BqMtEjsYvlvRjSZ/Y3lcsWy9phe35kkLSEUk/aUmHAGoxkaPxv5fkMUrv1N8OgFbhDDogCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EASjoj2bcwekvQ/oxbNlXSmbQ1cmW7trVv7kuitWXX2dltEjPn7b20N+7c2bg9GRH/HGijRrb11a18SvTWrXb3xNh5IgrADSXQ67Fs6vP0y3dpbt/Yl0Vuz2tJbRz+zA2ifTu/ZAbQJYQeS6EjYbT9g+1Pbn9t+uhM9NGL7iO1PimmoBzvcy1bbp23vH7Vsju33bB8qrsecY69DvXXFNN4l04x39LXr9PTnbf/MbnuKpM8k/aOkY5I+kLQiIv6rrY00YPuIpP6I6PgJGLZ/IOlPkv49Iv6uWPaCpC8jYmPxH+XsiPiXLuntOUl/6vQ03sVsRb2jpxmXNCDpn9XB166kr39SG163TuzZF0n6PCIOR8SfJf1K0tIO9NH1ImKvpC8vW7xU0rbi9jaN/LG0XYPeukJEnIiIj4rbX0m6OM14R1+7kr7aohNhv1XSH0fdP6bumu89JP3O9oe2V3W6mTHcFBEnpJE/Hkk3drify407jXc7XTbNeNe8ds1Mf15VJ8I+1lRS3TT+tzgiFkp6UNKTxdtVTMyEpvFulzGmGe8KzU5/XlUnwn5MUt+o+9+V9EUH+hhTRHxRXJ+W9Ia6byrqUxdn0C2uT3e4n//XTdN4jzXNuLrgtevk9OedCPsHkm63/T3b35G0XNLODvTxLbanFwdOZHu6pCXqvqmod0paWdxeKemtDvZyiW6ZxrvRNOPq8GvX8enPI6LtF0kPaeSI/H9LeqYTPTTo628k/WdxOdDp3iTt0Mjbuv/VyDuixyT9taTdkg4V13O6qLf/kPSJpI81EqzeDvV2r0Y+Gn4saV9xeajTr11JX2153ThdFkiCM+iAJAg7kARhB5Ig7EAShB1IgrADSRB2IIn/A4zA71w4qxB9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAANPUlEQVR4nO3db6xU9Z3H8c9HaGMUEmG5AloCbOODmk2WNhPEsKlumm3EJ1ijGzAhrCF7ScSEJjxYdR8U4xOzsa0YNzWXlRTWrk1jq/KA7JYgiSEm1cGwipJdES+USi5D/FMbNRX47oN72FzxzpnLzJk/3O/7lUxm5nzPmfPN4X44c+c35/4cEQIw/V3R7wYA9AZhB5Ig7EAShB1IgrADSczs5c7mzZsXS5Ys6eUugVRGR0d15swZT1brKOy2b5O0TdIMSf8WEY+Wrb9kyRLV6/VOdgmgRK1Wa1pr+2287RmS/lXSKkk3Slpr+8Z2Xw9Ad3XyO/tySUcj4lhE/FnSLyWtrqYtAFXrJOzXS/r9hOcni2VfYnvYdt12vdFodLA7AJ3oJOyTfQjwle/eRsRIRNQiojY0NNTB7gB0opOwn5S0aMLzb0h6v7N2AHRLJ2F/TdINtpfa/rqkNZJ2V9MWgKq1PfQWEWdt3y/pvzQ+9LYjIt6qrDMAleponD0i9kjaU1EvALqIr8sCSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkuhoymbbo5I+kXRO0tmIqFXRFIDqdRT2wt9GxJkKXgdAF/E2Hkii07CHpN/aPmh7eLIVbA/brtuuNxqNDncHoF2dhn1lRHxH0ipJm2x/9+IVImIkImoRURsaGupwdwDa1VHYI+L94v60pOclLa+iKQDVazvstq+2PfvCY0nfl3S4qsYAVKuTT+PnS3re9oXX+Y+I+M9KugJQubbDHhHHJP11hb0A6CKG3oAkCDuQBGEHkiDsQBKEHUiiigthcBn77LPPSut79+4trW/evLm0fuLEiUvu6YKRkZHS+oYNG9p+7Yw4swNJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoyzT3PHjh0rrd93332l9Vbj7MUlzm3Xy2zcuLG0fuTIkdL6Y4891va+pyPO7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBOPs08Cnn37atLZ27drSbQ8ePFhaX7BgQWn93nvvLa3ffffdTWvbt28v3fapp54qre/bt6+0fu7cuaa1GTNmlG47HXFmB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkGGefBu66666mtVbj6GvWrCmtP/PMM231NBUPP/xwaf25554rrb/77rul9bGxsaa16667rnTb6ajlmd32DtunbR+esGyu7b223ynu53S3TQCdmsrb+J9Luu2iZQ9I2hcRN0jaVzwHMMBahj0iXpb0wUWLV0vaWTzeKemOivsCULF2P6CbHxGnJKm4v7bZiraHbddt1xuNRpu7A9Cprn8aHxEjEVGLiNrQ0FC3dwegiXbDPmZ7oSQV96erawlAN7Qb9t2S1heP10t6sZp2AHRLy3F2289KulXSPNsnJf1I0qOSfmV7g6QTkppftIyumzVrVtvbthpnH2Rz5pSP+M6ePbtHnVweWoY9Ipr99YPvVdwLgC7i67JAEoQdSIKwA0kQdiAJwg4kwSWu00BEtFWTpLlz55bWz549W1o/c+ZMaX3Xrl1Na6+++mrptosXLy6tv/DCC6V1ht6+jDM7kARhB5Ig7EAShB1IgrADSRB2IAnCDiTBOPs08MorrzSt2S7dduvWraX1VuP0+/fvL62XOXDgQGl9xYoVbb82voozO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kwTj7NLBw4cKmtY8//rh025deeqm03mqcvdU4ftmfuV66dGnptqgWZ3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIJx9mmgXq83rb333nul246NjZXWV65c2VZPF2zYsKFpbf78+R29Ni5NyzO77R22T9s+PGHZVtt/sH2ouN3e3TYBdGoqb+N/Lum2SZb/NCKWFbc91bYFoGotwx4RL0v6oAe9AOiiTj6gu9/2G8Xb/DnNVrI9bLtuu95oNDrYHYBOtBv2n0n6pqRlkk5J+nGzFSNiJCJqEVEbGhpqc3cAOtVW2CNiLCLORcR5SdslLa+2LQBVayvstideU/kDSYebrQtgMLQcZ7f9rKRbJc2zfVLSjyTdanuZpJA0KmljF3tEB1pdM37llVd2df8PPvhgV18fU9cy7BGxdpLFT3ehFwBdxNdlgSQIO5AEYQeSIOxAEoQdSIJLXJM7evRoaf38+fOl9Suu4HxxueBfCkiCsANJEHYgCcIOJEHYgSQIO5AEYQeSYJw9uVaXuLYaR7/zzjtL6zNn8iM2KDizA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EASDIJOc6dPny6tP/nkk6X1BQsWlNa3bNlSWr/mmmtK6+gdzuxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kATj7NPA559/3rR28803l257/Pjx0vquXbtK6ytWrCitY3C0PLPbXmR7v+0jtt+yvblYPtf2XtvvFPdzut8ugHZN5W38WUlbIuJbklZI2mT7RkkPSNoXETdI2lc8BzCgWoY9Ik5FxOvF408kHZF0vaTVknYWq+2UdEe3mgTQuUv6gM72EknflvQ7SfMj4pQ0/h+CpGubbDNsu2673mg0OusWQNumHHbbsyT9WtIPI+KPU90uIkYiohYRtaGhoXZ6BFCBKYXd9tc0HvRfRMRvisVjthcW9YWSyi+vAtBXLYfebFvS05KORMRPJpR2S1ov6dHi/sWudIiWHnnkkaa1VkNrmzZtKq3fc889bfWEwTOVcfaVktZJetP2oWLZQxoP+a9sb5B0QtLd3WkRQBVahj0iDkhyk/L3qm0HQLfwdVkgCcIOJEHYgSQIO5AEYQeS4BLXy8Dbb79dWt+2bVvT2lVXXVW67bp169rqCZcfzuxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kATj7APgo48+Kq3fcsstbb/23r17S+u1Wq3t18blhTM7kARhB5Ig7EAShB1IgrADSRB2IAnCDiTBOHsPfPHFF6X1xx9/vLT+4YcfltaHh4eb1m666abSbZEHZ3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSGIq87MvkrRL0gJJ5yWNRMQ221sl/aOkRrHqQxGxp1uNXs727Ck/LGXzq0vSqlWrSutPPPHEJfeEfKbypZqzkrZExOu2Z0s6aPvCX0T4aUQ81r32AFRlKvOzn5J0qnj8ie0jkq7vdmMAqnVJv7PbXiLp25J+Vyy63/YbtnfYntNkm2Hbddv1RqMx2SoAemDKYbc9S9KvJf0wIv4o6WeSvilpmcbP/D+ebLuIGImIWkTUhoaGKmgZQDumFHbbX9N40H8REb+RpIgYi4hzEXFe0nZJy7vXJoBOtQy7bUt6WtKRiPjJhOULJ6z2A0mHq28PQFWm8mn8SknrJL1p+1Cx7CFJa20vkxSSRiVt7EqHl4Hjx4+X1jdv3lxabzV01mpa5ZkzuVIZrU3l0/gDkjxJiTF14DLCN+iAJAg7kARhB5Ig7EAShB1IgrADSTBAW4HFixeX1kdHR3vTCFCCMzuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJOGI6N3O7IakiRd/z5N0pmcNXJpB7W1Q+5LorV1V9rY4Iib9+289DftXdm7XI6LWtwZKDGpvg9qXRG/t6lVvvI0HkiDsQBL9DvtIn/dfZlB7G9S+JHprV0966+vv7AB6p99ndgA9QtiBJPoSdtu32f4f20dtP9CPHpqxPWr7TduHbNf73MsO26dtH56wbK7tvbbfKe4nnWOvT71ttf2H4tgdsn17n3pbZHu/7SO237K9uVje12NX0ldPjlvPf2e3PUPS/0r6O0knJb0maW1EvN3TRpqwPSqpFhF9/wKG7e9K+pOkXRHxV8Wyf5H0QUQ8WvxHOSci/mlAetsq6U/9nsa7mK1o4cRpxiXdIekf1MdjV9LX36sHx60fZ/blko5GxLGI+LOkX0pa3Yc+Bl5EvCzpg4sWr5a0s3i8U+M/LD3XpLeBEBGnIuL14vEnki5MM97XY1fSV0/0I+zXS/r9hOcnNVjzvYek39o+aHu4381MYn5EnJLGf3gkXdvnfi7WchrvXrpomvGBOXbtTH/eqX6EfbKppAZp/G9lRHxH0ipJm4q3q5iaKU3j3SuTTDM+ENqd/rxT/Qj7SUmLJjz/hqT3+9DHpCLi/eL+tKTnNXhTUY9dmEG3uD/d537+3yBN4z3ZNOMagGPXz+nP+xH21yTdYHup7a9LWiNpdx/6+ArbVxcfnMj21ZK+r8Gbinq3pPXF4/WSXuxjL18yKNN4N5tmXH0+dn2f/jwien6TdLvGP5F/V9I/96OHJn39paT/Lm5v9bs3Sc9q/G3dFxp/R7RB0l9I2ifpneJ+7gD19u+S3pT0hsaDtbBPvf2Nxn81fEPSoeJ2e7+PXUlfPTlufF0WSIJv0AFJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEv8HhVr4ARAFvh8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAANcElEQVR4nO3df4hd9ZnH8c/H2BhJFWMzZgeNTjeIrC6YlkHWuBSlWn9CDGJIICUVcaIk0EL/WHUlHUQkyKYlyCqma2i6dFMCqRgkug0hEKpQHCWbxI2r2ZBNE5PJBMUaELtJnv1jTpZJnHvu5N5zf0ye9wuGe+95zpnzcJjPnDP3e+58HRECcOG7qNMNAGgPwg4kQdiBJAg7kARhB5K4uJ07mzlzZvT19bVzl0AqBw4c0PHjxz1eramw275H0hpJUyT9S0SsKlu/r69PQ0NDzewSQIn+/v6atYYv421PkfTPku6VdKOkxbZvbPT7AWitZv5mv0XSvojYHxF/kfRbSfOraQtA1ZoJ+9WS/jTm9aFi2VlsD9gesj00MjLSxO4ANKOZsI/3JsDX7r2NiLUR0R8R/T09PU3sDkAzmgn7IUmzx7y+RtInzbUDoFWaCfu7kq63/W3bUyUtkrS5mrYAVK3hobeIOGl7haR/1+jQ27qI+KCyzgBUqqlx9ojYImlLRb0AaCFulwWSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4k0dSUzbYPSPpC0ilJJyOiv4qmAFSvqbAX7oiI4xV8HwAtxGU8kESzYQ9Jv7f9nu2B8VawPWB7yPbQyMhIk7sD0Khmw35bRHxX0r2Sltv+3rkrRMTaiOiPiP6enp4mdwegUU2FPSI+KR6PSXpN0i1VNAWgeg2H3fZ025edeS7pB5L2VNUYgGo18278LEmv2T7zff4tIt6qpCtMGp9//nlp/cUXX6xZe/bZZ0u3PXnyZGl9+fLlDe87o4bDHhH7Jd1cYS8AWoihNyAJwg4kQdiBJAg7kARhB5Ko4oMwuIBt3bq1tL5w4cLS+ty5c2vWtmzZUrrtvn37Suv1ht5WrlxZs5bxbk7O7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBOPsF7hTp06V1tevX19aX7FiRWn9+eefL60//vjjNWvTpk0r3faGG24ordcbZ6/3/bPhzA4kQdiBJAg7kARhB5Ig7EAShB1IgrADSTDOfoF78803S+uPPfZYaX3jxo2l9Yceeui8e5qobdu2ldZ7e3tL65dddlmV7Ux6nNmBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnG2S8AX375Zc1avXH0ZcuWldYXLFjQUE8TUW+65xdeeKFl+86o7pnd9jrbx2zvGbPsSttbbX9cPM5obZsAmjWRy/hfSbrnnGVPStoWEddL2la8BtDF6oY9InZI+vScxfMlnfl/RuslPVhxXwAq1ugbdLMi4ogkFY9X1VrR9oDtIdtDIyMjDe4OQLNa/m58RKyNiP6I6M84mR7QLRoN+7DtXkkqHo9V1xKAVmg07JslLS2eL5X0ejXtAGiVuuPstjdIul3STNuHJP1M0ipJG20/KumgpIdb2WR2p0+fLq3PmzevZu2aa64p3Xb16tWl9Ysuat1fekuWLCmtf/jhh6X1VatWVdnOBa9u2CNicY3S9yvuBUALcbsskARhB5Ig7EAShB1IgrADSfAR10ngnXfeKa3v2rWrZm3v3r2l21566aUN9TRRO3bsqFl76623Sre9/PLLS+v1pmzG2TizA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EASjLNPAi+//HJp/eabb65ZmzNnTtXtnOXEiROl9UWLFtWs1fvo7uDgYGl9+vTppXWcjTM7kARhB5Ig7EAShB1IgrADSRB2IAnCDiTBOPsksGHDhtL6K6+8UrM2ZcqUpvb91VdfldYXLlxYWh8eHq5Ze+aZZ0q3feKJJ0rrOD+c2YEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcbZu0C9/+1ezwMPPNDwtrt37y6tL1iwoLS+f//+0vpNN91Us/bUU0+VbnvJJZeU1nF+6p7Zba+zfcz2njHLBm0ftr2z+LqvtW0CaNZELuN/JemecZb/IiLmFl9bqm0LQNXqhj0idkj6tA29AGihZt6gW2F7V3GZP6PWSrYHbA/ZHhoZGWlidwCa0WjYX5Y0R9JcSUckra61YkSsjYj+iOjv6elpcHcAmtVQ2CNiOCJORcRpSb+UdEu1bQGoWkNht9075uUCSXtqrQugO9QdZ7e9QdLtkmbaPiTpZ5Jutz1XUkg6IGlZC3u84F1xxRWl9XpzqN9///01a5999lnptkePHi2tT5s2rbQeEaX1lStXNvy9Ua26YY+IxeMsfrUFvQBoIW6XBZIg7EAShB1IgrADSRB2IAk+4toFent7S+tvvPFGaX3NmjU1a/PmzSvd9pFHHimtP/zww6X1O++8s7Re7yOyaB/O7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBOPsk8Add9zRcL3eR1Cfe+650vrhw4dL69u3by+tX3wxP2LdgjM7kARhB5Ig7EAShB1IgrADSRB2IAnCDiTBIOgF7qOPPiqtDw4OltbLPisvSdddd935toQO4cwOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kwzn6Bu+uuu0rr1157bWl9YGCgynbQQXXP7LZn295ue6/tD2z/uFh+pe2ttj8uHme0vl0AjZrIZfxJST+NiL+R9HeSltu+UdKTkrZFxPWSthWvAXSpumGPiCMR8X7x/AtJeyVdLWm+pPXFauslPdiqJgE077zeoLPdJ+k7kv4oaVZEHJFGfyFIuqrGNgO2h2wPjYyMNNctgIZNOOy2vylpk6SfRMSfJ7pdRKyNiP6I6O/p6WmkRwAVmFDYbX9Do0H/TUT8rlg8bLu3qPdKOtaaFgFUoe7Qm21LelXS3oj4+ZjSZklLJa0qHl9vSYeo6+DBgzVrR48eLd1206ZNpfWpU6c21BO6z0TG2W+T9ENJu23vLJY9rdGQb7T9qKSDkson8gbQUXXDHhF/kOQa5e9X2w6AVuF2WSAJwg4kQdiBJAg7kARhB5LgI66TwMmTJ0vrS5YsqVnr6+sr3fbuu+9upCVMQpzZgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJxtkngS1btpTW33777Zq1Q4cOlW7L59Xz4MwOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kwzj4JDA4OltZvvfXWmrVZs2ZV3A0mK87sQBKEHUiCsANJEHYgCcIOJEHYgSQIO5DEROZnny3p15L+StJpSWsjYo3tQUmPSRopVn06Iso/eI2GDA8Pl9ZfeumlmrWLLuL3OUZN5Kaak5J+GhHv275M0nu2txa1X0TEP7WuPQBVmcj87EckHSmef2F7r6SrW90YgGqd1zWe7T5J35H0x2LRCtu7bK+zPaPGNgO2h2wPjYyMjLcKgDaYcNhtf1PSJkk/iYg/S3pZ0hxJczV65l893nYRsTYi+iOiv6enp4KWATRiQmG3/Q2NBv03EfE7SYqI4Yg4FRGnJf1S0i2taxNAs+qG3bYlvSppb0T8fMzy3jGrLZC0p/r2AFRlIu/G3ybph5J2295ZLHta0mLbcyWFpAOSlrWkQ+jw4cOdbgEXgIm8G/8HSR6nxJg6MIlwxwWQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJR0T7dmaPSPqfMYtmSjretgbOT7f21q19SfTWqCp7uy4ixv3/b20N+9d2bg9FRH/HGijRrb11a18SvTWqXb1xGQ8kQdiBJDod9rUd3n+Zbu2tW/uS6K1Rbemto3+zA2ifTp/ZAbQJYQeS6EjYbd9j+79s77P9ZCd6qMX2Adu7be+0PdThXtbZPmZ7z5hlV9reavvj4nHcOfY61Nug7cPFsdtp+74O9Tbb9nbbe21/YPvHxfKOHruSvtpy3Nr+N7vtKZI+knSXpEOS3pW0OCL+s62N1GD7gKT+iOj4DRi2vyfphKRfR8TfFstekPRpRKwqflHOiIh/6JLeBiWd6PQ03sVsRb1jpxmX9KCkH6mDx66kr4Vqw3HrxJn9Fkn7ImJ/RPxF0m8lze9AH10vInZI+vScxfMlrS+er9foD0vb1eitK0TEkYh4v3j+haQz04x39NiV9NUWnQj71ZL+NOb1IXXXfO8h6fe237M90OlmxjErIo5Ioz88kq7qcD/nqjuNdzudM8141xy7RqY/b1Ynwj7eVFLdNP53W0R8V9K9kpYXl6uYmAlN490u40wz3hUanf68WZ0I+yFJs8e8vkbSJx3oY1wR8UnxeEzSa+q+qaiHz8ygWzwe63A//6+bpvEeb5pxdcGx6+T0550I+7uSrrf9bdtTJS2StLkDfXyN7enFGyeyPV3SD9R9U1FvlrS0eL5U0usd7OUs3TKNd61pxtXhY9fx6c8jou1fku7T6Dvy/y3pHzvRQ42+/lrSfxRfH3S6N0kbNHpZ978avSJ6VNK3JG2T9HHxeGUX9favknZL2qXRYPV2qLe/1+ifhrsk7Sy+7uv0sSvpqy3HjdtlgSS4gw5IgrADSRB2IAnCDiRB2IEkCDuQBGEHkvg/jbsAM0CS124AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAANIUlEQVR4nO3dYaxU9ZnH8d9voRgjVREugpQstfpizZqleINr3BCNbgP6AnjRTXmBbNDQF5K0SY1L6gtMjNGYbZuNWWvoQspuqk1ja8REuyhpJH1h45WwgJJd9IYtIIF7QYNAIgjPvrjHzRXvnBlmzswZeL6f5GZmzjNnzpMTfpyZ8z8zf0eEAFz+/qLuBgD0BmEHkiDsQBKEHUiCsANJTO7lxmbMmBHz5s3r5SaBVPbv36/R0VFPVOso7LYXS/oXSZMk/VtEPF32/Hnz5mloaKiTTQIoMTg42LDW9tt425Mk/aukJZJukbTC9i3tvh6A7urkM/tCSR9ExHBEnJH0a0lLq2kLQNU6CfscSQfGPT5YLPsS22tsD9keGhkZ6WBzADrRSdgnOgnwlWtvI2JDRAxGxODAwEAHmwPQiU7CflDS3HGPvyHpo87aAdAtnYT9HUk32/6m7SmSvidpSzVtAaha20NvEfG57bWS/lNjQ2+bIuK9yjoDUKmOxtkj4jVJr1XUC4Au4nJZIAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeS6OlPSePys3bt2tL61q1b26pJY79GjOpwZAeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJBhnR0c+/PDD0vq+ffsa1u6///7SdXft2lVanzRpUmkdX8aRHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSYJwdpU6cOFFa37ZtW9uv/f7775fWz58/X1pnnP3idBR22/slfSrpnKTPI2KwiqYAVK+KI/vdETFawesA6CI+swNJdBr2kLTV9ru210z0BNtrbA/ZHhoZGelwcwDa1WnY74yIBZKWSHrY9qILnxARGyJiMCIGBwYGOtwcgHZ1FPaI+Ki4PSrpZUkLq2gKQPXaDrvtq2x//Yv7kr4jaU9VjQGoVidn46+X9LLtL17nhYj4fSVdoW+cOnWqtH727Nm2X3v16tWl9cmTuQykSm3vzYgYlvQ3FfYCoIsYegOSIOxAEoQdSIKwA0kQdiAJxjaSa/Y10qeeeqpr237ooYdK68WwLirCkR1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkmCcPbknnniitP7ss8/2qBN0G0d2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCcfbL3Ouvv15af+aZZ3rUCerGkR1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkmCc/TLw5ptvNqwtW7asdN0zZ86U1hctWlRa3759e2kd/aPpkd32JttHbe8Zt+w622/Y3lfcTutumwA61crb+F9KWnzBsnWStkXEzZK2FY8B9LGmYY+I7ZKOX7B4qaTNxf3NksrfKwKoXbsn6K6PiMOSVNzObPRE22tsD9keGhkZaXNzADrV9bPxEbEhIgYjYnBgYKDbmwPQQLthP2J7tiQVt0erawlAN7Qb9i2SVhX3V0l6pZp2AHRL03F22y9KukvSDNsHJa2X9LSk39h+UNKfJX23m032g7Lx6AMHDpSuOzw8XFp/6623SuvPP/98af3YsWOl9TIvvfRSaf3uu+8urU+fPr3tbaO3moY9IlY0KN1TcS8AuojLZYEkCDuQBGEHkiDsQBKEHUiCr7i26JNPPmlYW758eem6u3fv7mjb06aVf6lw3brG30N65JFHOnrtEydOlNZx6eDIDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJMM7eopkzG/7ylnbs2FG67ujoaEfbvvLKK0vr11xzTUevX5fTp0/X3UIqHNmBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnG2SsweXL5bpw1a1aPOqneFVdcUVqfM2dOaf3QoUMNay+88ELpuvfcww8YV4kjO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kwTg7SjX7Lv2tt95aWi8bZ1+yZElbPaE9TY/stjfZPmp7z7hlj9s+ZHtn8Xdfd9sE0KlW3sb/UtLiCZb/LCLmF3+vVdsWgKo1DXtEbJd0vAe9AOiiTk7QrbW9q3ib33DCMNtrbA/ZHhoZGelgcwA60W7Yfy7pW5LmSzos6SeNnhgRGyJiMCIGBwYG2twcgE61FfaIOBIR5yLivKRfSFpYbVsAqtZW2G3PHvdwuaQ9jZ4LoD80HWe3/aKkuyTNsH1Q0npJd9meLykk7Zf0/S72iMtUs+/Co1pNwx4RKyZYvLELvQDoIi6XBZIg7EAShB1IgrADSRB2IAm+4oqust2wNn369B52Ao7sQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AE4+yXgGPHjpXWT5482fZrT5kypbQ+derU0vr69etL68PDww1rx4+X/7Th6dOnS+tnz54trZf1tnLlytJ1b7vtttL6pYgjO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kwTh7Bc6dO1dabzZOvnFj+Y/1Pvnkk6X1U6dOldbLNBtnv/rqq0vro6OjbW/7jjvuKK3PmjWrtP7ZZ5+V1j/++OOGtRtuuKF0XcbZAVyyCDuQBGEHkiDsQBKEHUiCsANJEHYgCcbZW1Q2lv3oo4+Wrvvcc89V3U7L5s6dW1ov+113qflY+IIFCy66p37wwAMP1N1CzzU9stuea/sPtvfafs/2D4rl19l+w/a+4nZa99sF0K5W3sZ/LulHEfFXkv5W0sO2b5G0TtK2iLhZ0rbiMYA+1TTsEXE4InYU9z+VtFfSHElLJW0unrZZ0rJuNQmgcxd1gs72PEnflvQnSddHxGFp7D8ESTMbrLPG9pDtoZGRkc66BdC2lsNue6qk30r6YUScaHW9iNgQEYMRMTgwMNBOjwAq0FLYbX9NY0H/VUT8rlh8xPbsoj5b0tHutAigCk2H3jw2NrNR0t6I+Om40hZJqyQ9Xdy+0pUO+8Srr77asNbtobVmP3v82GOPNazddNNNpetOmjSprZ5w6WllnP1OSSsl7ba9s1j2Y42F/De2H5T0Z0nf7U6LAKrQNOwR8UdJja68uKfadgB0C5fLAkkQdiAJwg4kQdiBJAg7kARfcW3R4sWLG9Zuv/320nVvvPHG0vrq1atL6/fee29pHWgFR3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIJx9hZde+21DWtvv/12DzsB2sORHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5JoGnbbc23/wfZe2+/Z/kGx/HHbh2zvLP7u6367ANrVyo9XfC7pRxGxw/bXJb1r+42i9rOI+OfutQegKq3Mz35Y0uHi/qe290qa0+3GAFTroj6z254n6duS/lQsWmt7l+1Ntqc1WGeN7SHbQyMjIx01C6B9LYfd9lRJv5X0w4g4Iennkr4lab7Gjvw/mWi9iNgQEYMRMTgwMFBBywDa0VLYbX9NY0H/VUT8TpIi4khEnIuI85J+IWlh99oE0KlWzsZb0kZJeyPip+OWzx73tOWS9lTfHoCqtHI2/k5JKyXttr2zWPZjSStsz5cUkvZL+n5XOgRQiVbOxv9RkicovVZ9OwC6hSvogCQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSTgiercxe0TS/45bNEPSaM8auDj92lu/9iXRW7uq7O0vI2LC33/radi/snF7KCIGa2ugRL/21q99SfTWrl71xtt4IAnCDiRRd9g31Lz9Mv3aW7/2JdFbu3rSW62f2QH0Tt1HdgA9QtiBJGoJu+3Ftv/b9ge219XRQyO299veXUxDPVRzL5tsH7W9Z9yy62y/YXtfcTvhHHs19dYX03iXTDNe676re/rznn9mtz1J0v9I+ntJByW9I2lFRLzf00YasL1f0mBE1H4Bhu1Fkk5K+veI+Oti2TOSjkfE08V/lNMi4p/6pLfHJZ2sexrvYrai2eOnGZe0TNI/qsZ9V9LXP6gH+62OI/tCSR9ExHBEnJH0a0lLa+ij70XEdknHL1i8VNLm4v5mjf1j6bkGvfWFiDgcETuK+59K+mKa8Vr3XUlfPVFH2OdIOjDu8UH113zvIWmr7Xdtr6m7mQlcHxGHpbF/PJJm1tzPhZpO491LF0wz3jf7rp3pzztVR9gnmkqqn8b/7oyIBZKWSHq4eLuK1rQ0jXevTDDNeF9od/rzTtUR9oOS5o57/A1JH9XQx4Qi4qPi9qikl9V/U1Ef+WIG3eL2aM39/L9+msZ7omnG1Qf7rs7pz+sI+zuSbrb9TdtTJH1P0pYa+vgK21cVJ05k+ypJ31H/TUW9RdKq4v4qSa/U2MuX9Ms03o2mGVfN+6726c8joud/ku7T2Bn5DyU9VkcPDfq6UdJ/FX/v1d2bpBc19rburMbeET0oabqkbZL2FbfX9VFv/yFpt6RdGgvW7Jp6+zuNfTTcJWln8Xdf3fuupK+e7DculwWS4Ao6IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUji/wAneu2dEjQM1AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# show some images\n",
    "print(mytrl[0:10])\n",
    "for i in range(10):\n",
    "    img = mytri[i].reshape((28,28))\n",
    "    plt.imshow(img, cmap=\"Greys\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider the primal problem, with regularized $h_\\varepsilon$, which we write:\n",
    "$$\n",
    "\\min_{w}\\sum_{i=1}^N \\frac{1}{2N}|w|^2+ \\lambda h_\\varepsilon(y_i(w\\cdot X_i)). \\quad (P)\n",
    "$$\n",
    "and letting $f_i(w):= |w|^2/(2N) +  \\lambda h_\\varepsilon(y_i(w\\cdot X_i))$ we will implement a randomized \n",
    "gradient descent, updating at each iteration $w^{k+1} = w^k -\\tau \\nabla f_i(w^k)$ for a randomly\n",
    "chosen $i$, and $\\tau \\to 0$ (or $\\tau \\sim 1/\\sqrt{\\mbox{number of iterations}}$).\n",
    "\n",
    "One has\n",
    "$$\n",
    "\\nabla f_i (w) = \\frac{w}{N} +  \\lambda h'_\\varepsilon(y_i(w\\cdot X_i)) y_i X_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solving a SVM by randomized gradient descent\n",
      "to classify digits  [9. 1. 6.]  and  [3. 7. 4. 8.]  from the MNIST database\n",
      "done\n"
     ]
    }
   ],
   "source": [
    "# randomized coordinate descent on the primal\n",
    "\n",
    "print('solving a SVM by randomized gradient descent')\n",
    "print('to classify digits ', mylabelsplus ,' and ',mylabelsminus ,' from the MNIST database')\n",
    "\n",
    "N= mytrl.shape[0] # number of elements\n",
    "epoch = 20 # number of \"epochs\" (= iterations / N)\n",
    "\n",
    "w = np.zeros(mytri.shape[1]) # initial w\n",
    "lbda = 1\n",
    "\n",
    "random.seed()\n",
    "\n",
    "# \"Huber\" type regularization for (1-t)^+ (epsilon 1 works)\n",
    "epsilon = .5 \n",
    "# step for gradient descent\n",
    "le = epsilon/(epsilon+lbda) * 10 / np.sqrt(N*epoch)\n",
    "\n",
    "def dh(x):\n",
    "    global epsilon\n",
    "    if x < 1-epsilon:\n",
    "        return -1\n",
    "    elif x < 1:\n",
    "        return (x-1)/epsilon\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "for it in range(epoch*N):\n",
    "    i = random.randint(0,N-1)\n",
    "    gra = w/N + mytrls[i]*mytri[i,:]*dh(mytrls[i]*np.dot(mytri[i,:],w))\n",
    "    w = w - le * gra\n",
    "    \n",
    "print('done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " train success =  [90.2] %\n",
      " fail =  [9.8] %\n",
      " test success =  [90.18] %\n",
      " fail =  [9.82] %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAUIUlEQVR4nO3dbYxc5XUH8P9/XvbFs/au1/YaYwwESquivjjVyo1EVVGljQhfIB9ShQ8plVBNpSAlalQV0UrhI6oaonyoojoFxUlTokgJAqmoiYWionwoYkEumDoUg3gxNrs2tvHu7O7s7Mzphx3axew9Z5m3O8vz/0nWrufMnXvmzpy5u3vu8zw0M4jIJ18h7wREpD9U7CKJULGLJELFLpIIFbtIIkr93FmxUrHyxGQ/dymSlPqlC2hUq9wo1lGxk7wNwLcBFAH8s5k95N2/PDGJ6+79q052KSKON//p4cxY2z/GkywC+EcAnwdwM4C7SN7c7uOJSG918jv7IQCnzOx1M1sB8CMAd3QnLRHptk6KfT+At9f9/3Trtg8heZjkDMmZRrXawe5EpBOdFPtGfwT4yLW3ZnbEzKbNbLpYqXSwOxHpRCfFfhrAgXX/vwbAmc7SEZFe6aTYnwNwE8lPkRwC8CUAT3YnLRHptrZbb2a2SvI+AD/DWuvtUTN7uWuZbSG2YVfzY2DvRh5GuXW8646ffO/08LBuSR312c3sKQBPdSkXEekhXS4rkggVu0giVOwiiVCxiyRCxS6SCBW7SCL6Op59K/PayVbwG7pWDB47ipeCx3c+sqNtO/24Z93vs7PhxFaDHn2Qeid9dDaDuJM3gI6vL8jjGgCd2UUSoWIXSYSKXSQRKnaRRKjYRRKhYhdJRDKtt05HYnrttU5ba83hoLUWtd6KTjzYlsN+j6lY8ntU5aFVNz5czo4z6D8t1YbceG2x7MZtJfvAs+af57jiv2GKNTcM+ocFaGY/fq/acjqziyRCxS6SCBW7SCJU7CKJULGLJELFLpIIFbtIIpLps4e9y2g4pdMXRdBT9YagAkAh2L5ZiuaDzg41tgVjORf9iwCawXFZ9nr8AJbK2fFCpe5uWwgee9euBTfecF6zhcVhd9t61e/x84Lf4y8Gr6mGuIpIz6jYRRKhYhdJhIpdJBEqdpFEqNhFEqFiF0nEJ6fPHvbJg3jQFy2uZMcKfrs4HM9eCKZjDsfiO3E2/J2PnvcPzNjbS2689NpZN96Yncvedt9V7rbzv3+tGz/3O9vceG2PN4+1uymGLvrnwfJC+9c+AEDTqTzrUVV29LAk3wAwD6ABYNXMpruRlIh0Xzc+Q/7IzM534XFEpIf0O7tIIjotdgPwc5LPkzy80R1IHiY5Q3KmUa12uDsRaVenP8bfYmZnSE4BOEbyV2b2zPo7mNkRAEcAYGT/gRwu/xcRoMMzu5mdaX2dA/A4gEPdSEpEuq/tYidZIbn9g+8BfA7AiW4lJiLd1cmP8XsBPE7yg8f5VzP7965klcX5JSBaYjfqhZf9odEYvpTdjx6q+r3q0oKfXGnJjzdG/F55dV/22OrqPv/z/PxBvyE8d2jUjVvhBjc+cu7XnI3dTdEY8ePDl/x4eSH7uEXzvkfXHxTrfnx5wj/uy5PZxz163u324dsudjN7HcDvtru9iPSXWm8iiVCxiyRCxS6SCBW7SCJU7CKJ+MQMcY2GsJYW/fjIBf8Bts1m9+6GZ/2+HZf8Pk9jvOLHh/3W29Lu7M/shRv9sbtje/3c/+ymZ934ntK8G/+PS7+RGTs+u9/ddv7VCTdeXPbPVd6w5E6H9kaMfv9sdVv2a9r0Z7GOOpaZdGYXSYSKXSQRKnaRRKjYRRKhYhdJhIpdJBEqdpFEJNNnL6743cnSsh8vLjv96qa/88ZOv49eH2+/Jwv4Sz5ve8t/iReG/OmYf3DKn4+k0fDPF4uz2c+9fNl/XiMX/eG3Y6eDYajOa145vexuW37rnBtv7tzh73vSXxKajezcjNHc4e3RmV0kESp2kUSo2EUSoWIXSYSKXSQRKnaRRKjYRRKxtfrsPZxKmub32ZennL7pHr+n2hjy+6bNsh+vV6Lts2MFZ0w3AJTmnI0BLF4cd+PD5/3zxe53s49rdG0E6L8mw5f9B2Aze/tCPXjDlIPSKPqvSWMkeM1K3jrb/q7dAe1OTGd2kUSo2EUSoWIXSYSKXSQRKnaRRKjYRRKhYhdJxGD12YMJsd2+bLCt+UOnsTzh36E+lt38XNrjP3Z9u98PHn4v+MwN+q7eNQbRUtSVM/6DL+31t2+M+gd+aCH7uXvj8AFgeWcwVn6P/5qZs/nqqD/HAK7341bwj1ttR9CH9y/N8HkP7cTCMzvJR0nOkTyx7rZJksdIvtr6uvNjJSsifbeZH+O/B+C2K267H8DTZnYTgKdb/xeRARYWu5k9A+DCFTffAeBo6/ujAO7scl4i0mXt/oFur5mdBYDW16msO5I8THKG5EyjWm1zdyLSqZ7/Nd7MjpjZtJlNFyvBH0VEpGfaLfZZkvsAoPV1rnspiUgvtFvsTwK4u/X93QCe6E46ItIrYZ+d5GMAbgWwm+RpAN8A8BCAH5O8B8BbAL7YlWw66bMHH1v1bX7fczX4DWPx2uym8Ohuf/F31vwx4zX6TddCzc9914nsA7f9bX9t+NqEn9v8Df6BtQl/ooDyf2Zvvzjlv/1WR/3nXd3vv2Hqe7JfMy76PfrSor/vKB7NrxBd99ELYbGb2V0Zoc92ORcR6SFdLiuSCBW7SCJU7CKJULGLJELFLpKIwRri2sFKtc2gldH0VybG8pQ/DHX7VfOZMbNgOONy0OYJWmulZT/u7Z6rfntqfr//FvjLP/6ZG/+3M7/txhen9mXGVrb7z2vhWj/3kRsvu/HRQvZrOn951N22Phy0BZ0hzwBQmvfPo8UVZ/toiu026cwukggVu0giVOwiiVCxiyRCxS6SCBW7SCJU7CKJGKw+e8AdFhgMj20O+XdojgdrOjuWFv0hqgz67NF0zI2K33idr3uP7/eT3//Mshv/68nX3Pi/nDrkxusHsvvJ5exLFwAAts/PbaXmv30bZ7Ofu9OCBwA0h6M3VHDtQ1BZ5rzd2MH1Jh6d2UUSoWIXSYSKXSQRKnaRRKjYRRKhYhdJhIpdJBGD1WcP+oveErxebDMK7/tTKi9eGM+McTWYVrgc9PgngrWLo8d3LkBY9dvs+PVrZt34383549Uvv7vdjW93VvxqRu++8/71CztO+Mel6MyiXdsZLFU9FbxmQe5Rq9ydFl19dhHphIpdJBEqdpFEqNhFEqFiF0mEil0kESp2kUQMVp+9A9ESuWz4zcuRuWAe8KXsWDRnPQrB0sNl/wFs1H9y3nNfDZaqPvXuHjf+yiv73fjQBT/3ocvZ/erq1X5uO4M++sRr/nLUi1NDmTHv9QSAUjVYCyCYgyCPJZkj4Zmd5KMk50ieWHfbgyTfIXm89e/23qYpIp3azI/x3wNw2wa3f8vMDrb+PdXdtESk28JiN7NnAFzoQy4i0kOd/IHuPpIvtn7M35l1J5KHSc6QnGlUnQulRaSn2i327wC4EcBBAGcBfDPrjmZ2xMymzWy6WKm0uTsR6VRbxW5ms2bWMLMmgO8C8KcYFZHctVXsJNevw/sFACey7isigyHss5N8DMCtAHaTPA3gGwBuJXkQa7O1vwHg3h7m+P+5OGOAGQwJL6748UIQ9/bdGPO3bWa3e1t38Hu65Yo/p/3i1dkv48h5//Pc3h1x4yVvHXEA284G/Wjn4Ycvupti5KI/uXt9zH/71sazc6sH66s3/aH04fwJFs3NUMzu0zN4P7QrLHYzu2uDmx/pQS4i0kO6XFYkESp2kUSo2EUSoWIXSYSKXSQRW2uIqzOqsBAMcS1V/SGJhWDF5vqO7HZII5iuuREs/8tgdeBGI2ifbc/uOy6X/LGWI+/68aK/anLMeW7FFf+J18b95x0NkV0dcx7fOpwqOljyORpy7bXmejSTtM7sIqlQsYskQsUukggVu0giVOwiiVCxiyRCxS6SiK3VZ/cEfc+oj15eCvquzorOxaFguKS/GjSawVTR5aL/5Ia8Pjv8sZqNIf/zvlkKetnBNQZe0zia3nu14r8m9algXHIt+7kVVvznzbqfWzTENeqz50FndpFEqNhFEqFiF0mEil0kESp2kUSo2EUSoWIXScSW6rO7U0lH44uDMePR2OoRZ9rjor9yMBrDwQjlIDyxY9GNXz12OTN2adxvhJ8emXDjjVqw9nAwZ3LhUvZbrFiLpqEOXrSIcyqL+uDRdRnR+ym67iPcvgd0ZhdJhIpdJBEqdpFEqNhFEqFiF0mEil0kESp2kURsqT67O9d20LeMxpRHhuazG7OloF+8ssM/zCuTfi/7vQv+mtAF58lXa/560ZPjVTdeq/u5L9f8A+tdvtAMxrMXFoIevzNeHQAKTrxUDfYdjWcPl2T243k02sMzO8kDJH9B8iTJl0l+tXX7JMljJF9tfd3Z+3RFpF2b+TF+FcDXzew3AXwGwFdI3gzgfgBPm9lNAJ5u/V9EBlRY7GZ21sxeaH0/D+AkgP0A7gBwtHW3owDu7FWSItK5j/UHOpLXA/g0gGcB7DWzs8DaBwKAqYxtDpOcITnTqPq/H4pI72y62EmOAfgJgK+ZWfbIiyuY2REzmzaz6WKl0k6OItIFmyp2kmWsFfoPzeynrZtnSe5rxfcBmOtNiiLSDWHrjSQBPALgpJk9vC70JIC7ATzU+vpETzJcz/loavgdprBXsjIWLIvsbB5Nt1xe8Nss478KhnoO+8NUF0az4ys7/H1Xx/yxmDbqx0sX/bdQydl91H0afi84LsE02J10t6Ih081o5G/Ueot6dz2wmT77LQC+DOAlksdbtz2AtSL/Mcl7ALwF4Iu9SVFEuiEsdjP7JbKnV/hsd9MRkV7R5bIiiVCxiyRCxS6SCBW7SCJU7CKJ2FpDXJ2PpmgIa7TEbjTdc9kZblmu+g3doSA+es6ft3j4jH/BIuezp5pu7h53t20O+Q3h5Sm/x1+o+7kvTmW/MKVa1Aj347Ud/mtWH8uON4PrMsI++paqnDU6s4skQsUukggVu0giVOwiiVCxiyRCxS6SCBW7SCK2VrfQaavGU/f64UbZ7+k2h7yebdDvDaYlbpb8iwQKK/5U0kOLy5kxnp71H7u24sYrb/p9dla2+Y+/uisz1hj2zzXLwRTbK0GfvTHixIbdTWFF//3AYBrsQaQzu0giVOwiiVCxiyRCxS6SCBW7SCJU7CKJULGLJGJr9dk7EE3TzVX/DsWl7Fih5j92MRi3XZvwP3MvXxeMKV89kL1vv40ezmnPVX/71W3BNQbbs2PLu/x9N7ZnL5MNAGAQd3rlrPk9/MJS8H7IvrRhbfvg/ZTDtPE6s4ukQsUukggVu0giVOwiiVCxiyRCxS6SCBW7SCI2sz77AQDfB3AVgCaAI2b2bZIPAvgLAOdad33AzJ7qVaK9Fq3l7c1L7/WSAaA2GfSigzXUG1PBmPPx7IsADuy86G57ueYM+gZQXfHH2tfrfnz3WDUzNjHsXLwAYCWYvP38YsWNv7/grFtvwYD2pWiCBF8effTIZi6qWQXwdTN7geR2AM+TPNaKfcvM/qF36YlIt2xmffazAM62vp8neRLA/l4nJiLd9bF+Zyd5PYBPA3i2ddN9JF8k+SjJnRnbHCY5Q3KmUc3+kU5EemvTxU5yDMBPAHzNzC4D+A6AGwEcxNqZ/5sbbWdmR8xs2symixX/dywR6Z1NFTvJMtYK/Ydm9lMAMLNZM2uYWRPAdwEc6l2aItKpsNhJEsAjAE6a2cPrbt+37m5fAHCi++mJSLds5q/xtwD4MoCXSB5v3fYAgLtIHsTaurpvALi3JxluVtDqiJZsjhYP9o5U2LaLpiUOhpGW3/HXF159PbuN9EolWLI5mEI7Oh0Uav6Bf3s4uy/5tv/QcfsqOK5wpnsuB0NYC8HQ4Ci56D2Rh838Nf6X2LiUtmxPXSRFuoJOJBEqdpFEqNhFEqFiF0mEil0kESp2kUQkM5V0nn3PQrC8rzU73UF2qLQYXYAQxKPTQXRcvf1H20bXRkR9eO9F73AM6iD20SM6s4skQsUukggVu0giVOwiiVCxiyRCxS6SCBW7SCJo1r+GIclzAN5cd9NuAOf7lsDHM6i5DWpegHJrVzdzu87M9mwU6Guxf2Tn5IyZTeeWgGNQcxvUvADl1q5+5aYf40USoWIXSUTexX4k5/17BjW3Qc0LUG7t6ktuuf7OLiL9k/eZXUT6RMUukohcip3kbSRfIXmK5P155JCF5BskXyJ5nORMzrk8SnKO5Il1t02SPEby1dbXDdfYyym3B0m+0zp2x0nenlNuB0j+guRJki+T/Grr9lyPnZNXX45b339nJ1kE8D8A/gTAaQDPAbjLzP67r4lkIPkGgGkzy/0CDJJ/CGABwPfN7Ldat/09gAtm9lDrg3Knmf3NgOT2IICFvJfxbq1WtG/9MuMA7gTw58jx2Dl5/Sn6cNzyOLMfAnDKzF43sxUAPwJwRw55DDwzewbAhStuvgPA0db3R7H2Zum7jNwGgpmdNbMXWt/PA/hgmfFcj52TV1/kUez78eGVf05jsNZ7NwA/J/k8ycN5J7OBvWZ2Flh78wCYyjmfK4XLePfTFcuMD8yxa2f5807lUewbTf41SP2/W8zs9wB8HsBXWj+uyuZsahnvftlgmfGB0O7y553Ko9hPAziw7v/XADiTQx4bMrMzra9zAB7H4C1FPfvBCrqtr3M55/N/BmkZ742WGccAHLs8lz/Po9ifA3ATyU+RHALwJQBP5pDHR5CstP5wApIVAJ/D4C1F/SSAu1vf3w3giRxz+ZBBWcY7a5lx5Hzscl/+3Mz6/g/A7Vj7i/xrAP42jxwy8roBwH+1/r2cd24AHsPaj3V1rP1EdA+AXQCeBvBq6+vkAOX2AwAvAXgRa4W1L6fc/gBrvxq+COB469/teR87J6++HDddLiuSCF1BJ5IIFbtIIlTsIolQsYskQsUukggVu0giVOwiifhfmTAZXnRxdDUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now we test the results\n",
    "succeed = (np.sum(mytrls * (mytri @ w) >= 0)) # if y_i and w.X_i have same sign then the classification is correct\n",
    "total = (mytrls.shape)\n",
    "\n",
    "print(' train success = ',np.round(100*succeed/total,2),'%\\n fail = ',np.round(100*(total-succeed)/total,2),'%')\n",
    "\n",
    "succeed = (np.sum(mytels * (mytei @ w) >= 0))\n",
    "total = (mytels.shape)\n",
    "\n",
    "print(' test success = ',np.round(100*succeed/total,2),'%\\n fail = ',np.round(100*(total-succeed)/total,2),'%')\n",
    "\n",
    "img = w.reshape((28,28))\n",
    "plt.imshow(img)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we try a different algorithm, that is, a stochastic coordinate ascent on the dual problem. We will even try a proximal coordinate ascent, which is shown to converge in [Peter Richtárik and Martin Takáč. Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function. Math. Program., 144(1-2, Ser. A):1–38, 2014.]\n",
    "\n",
    "The problem is:\n",
    "$$\n",
    "- \\lambda \\min_{(z_i)_{i=1}^N} \\frac{\\lambda}{2}\\left| \\sum_{i=1}^N z_i y_i X_i\\right|^2 + \\sum_{i=1}^N h^*(z_i)\n",
    "\\quad (D)\n",
    "$$\n",
    "Now, the gradient wr $z_i$ of the quadratic term is given by\n",
    "$$\n",
    "g_i(z):=\\lambda y_i\\sum_{j=1}^N z_j y_j (X_i\\cdot X_j)\n",
    "$$\n",
    "so that, in order to solve the problem, one should compute all the scalar products $X_i\\cdot X_j$ (on the\n",
    "other hand observe that it is the only information we need about the $X_i$'s to solve this problem, which \n",
    "will allow to replace this scalar product by arbitrary kernels $k(X_i,X_j)$, in more efficient \"kernel-SVM\"\n",
    "methods).\n",
    "\n",
    "The Lipschitz constant of this expression with respect to $z_i$ is $\\lambda y_i^2 |X_i|^2=\\lambda$. One\n",
    "should therefore take a step $\\tau=1/\\lambda$.\n",
    "Then, the optimization consists in picking randomly $i$ and updating $z_i^{k+1}=\n",
    "\\text{prox}_{\\tau h^*} (z^k_i- \\tau g_i(z^k))$. Here for $\\varepsilon \\ge 0$,\n",
    "$$\n",
    "\\text{prox}_{\\tau h_\\varepsilon^*} (t) = \\max\\left\\{\\min\\left\\{ \\frac{t-\\tau}{1+\\epsilon\\tau}, 0\\right\\}, -1\\right\\}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computing 648160681 scalar products\n",
      "done\n",
      "solving a SVM by randomized coordinate ascent on the dual\n",
      "to classify digits  [3. 7. 4. 8.]  and  [9. 1. 6.]  from the MNIST database\n",
      "done\n"
     ]
    }
   ],
   "source": [
    "N = mytrl.shape[0]\n",
    "epoch = 30\n",
    "\n",
    "z = np.zeros(N)\n",
    "\n",
    "print('computing',N*N,'scalar products')\n",
    "# reduce the size of array (float64 is too large for 60000 images, even 32)\n",
    "mytri = np.float32(mytri)\n",
    "scapro = mytri @ np.transpose(mytri) ## this is a bottleneck for too large arrays\n",
    "print('done')\n",
    "\n",
    "import random\n",
    "random.seed()\n",
    "\n",
    "# step for gradient descent\n",
    "lbda = 1\n",
    "tau = 1/lbda\n",
    "epsilon = 0\n",
    "\n",
    "# randomized coordinate descent on the dual\n",
    "\n",
    "print('solving a SVM by randomized coordinate ascent on the dual')\n",
    "print('to classify digits ',mylabelsminus,' and ',mylabelsplus,' from the MNIST database')\n",
    "\n",
    "for it in range(epoch*N):\n",
    "    i = random.randint(0,N-1)\n",
    "    gra = lbda*mytrls[i]*np.sum(z*scapro[i,:]*mytrls)\n",
    "    z[i] = np.minimum(np.maximum(z[i]-tau*gra - tau,-1),0)/(1+tau*epsilon) # the prox\n",
    "\n",
    "# compute finally w\n",
    "\n",
    "w = -np.sum(np.reshape(z*mytrls,(-1,1))*mytri,axis=0)\n",
    "\n",
    "print('done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " train success =  [91.3] %\n",
      " fail =  [8.7] %\n",
      " test success =  [90.81] %\n",
      " fail =  [9.19] %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAVmUlEQVR4nO3dW2yk5XkH8P8zJ3ttr9frPWF2l/OSsmmbJbXogaQijRoRIhVykSqrqqIS6eYiSImSiyJ6ES5R1RDlooq0ARpSpUSREgQXtA3dIlFuEOa8sMCyy4Y923uy1zP2eA5PLzwUB/z+XzMHzyzv/ydZtued7/ue+Wae+cZ+3oO5O0Tkky/T7QBEZHUo2UUSoWQXSYSSXSQRSnaRRORW82DZgUHPj4yu5iFFklK5cA61UtGWa2sp2c3sVgA/ApAF8KC738/unx8ZxVXf+G4rhxQR4siDDwTbmv4Yb2ZZAP8C4MsAdgLYbWY7m92fiHRWK3+z3wTgHXc/7O4LAH4B4Pb2hCUi7dZKsm8FcHTJ78cat/0OM9tjZhNmNlErFVs4nIi0opVkX+6fAB/pe+vue9193N3HswODLRxORFrRSrIfA7B9ye/bAJxoLRwR6ZRWkv15ADvM7GozKwD4OoAn2hOWiLRb06U3d6+a2d0A/guLpbeH3f31tkXWYzxL2iJvmZlqZN+R7esFPjLRaqStumzJ9YP2Oj92R/HQoufFs7HzEj4AO2efVC3V2d39SQBPtikWEekgdZcVSYSSXSQRSnaRRCjZRRKhZBdJhJJdJBGrOp69oyI124925P1427Oabq2fb1uJ1Mnra3jRNzPIC/WZDImtyt/PMzleaK9NF2h7dpZ0QADvY5Cd5yc92v8g8ur1vvB5scjrIVOO9E+I1Om72n8hQFd2kUQo2UUSoWQXSYSSXSQRSnaRRCjZRRLxiSm9eaz0FmmPlUpYmac6zOsw1h8prWX5wf0cL39VWWkvUlrrGyzT9lI/L631b+FTjRVy4drbNevP8mNX+eM+dXEtbb9wnsyMNJun2/p6ft4yc/y8WGRYc34mfJ3tVNlOV3aRRCjZRRKhZBdJhJJdJBFKdpFEKNlFEqFkF0nEJ6bO3ukhh2xIZGGK11w9z9vRamzksWeq/CmuHee17HyejwWdm+b16trWcB3+pSPbg20AUK/w85bt58Vsnws/9szwAt22f02Fts/N8HHNuSl+Xmi/j1ifkSbpyi6SCCW7SCKU7CKJULKLJELJLpIIJbtIIpTsIom4pOrsbGrhWB09NnVwjg/LRq4ULn4WLvKdz27jhdMML+mizkvhfMnmSP+DdYf5iRs+NEvbM0U+Hr524GCwLXclr7NP3bKNtpdH+YmZ+VS4Dm+neJ28nO+j7f3nIlN0R8azs+eltoZv22y/jJaS3cyOALgIoAag6u7jrexPRDqnHVf2L7j7mTbsR0Q6SH+ziySi1WR3AL8xsxfMbM9ydzCzPWY2YWYTtVLkD2MR6ZhWP8bf7O4nzGwzgKfM7E13f2bpHdx9L4C9ANB/+fbYimsi0iEtXdnd/UTj+ySAxwDc1I6gRKT9mk52Mxs0s7Xv/wzgSwD2tyswEWmvVj7GbwHwmJm9v59/d/f/bEtUAax2mZvj2+aK/C+IfKR9YCp88FqBv2eOHOKF9Pz+39L2+T+6hrbPXBEeOz3DN8Wpr/Bx3admyNzrALzAi8L9x/8s2FbeEOkEkInMeT/Jx7uveyP88s7O8+d78HSs4wZvn1vPY2N9J8qRAe31yFD5kKaT3d0PA/hMs9uLyOpS6U0kEUp2kUQo2UUSoWQXSYSSXSQRvTXEtYX+dRk+0jJamht+l98hf2r6Y0b0gcrYCG2fv/Fq2j47xmsts9vDpZrqWl7e8ln+EhjZzh/3Z7cco+2vbR0LtpXf2UC3zc3wElRs+G6VVAWzkddLdp6X1rILvD1X4u3VgfB1tp7lz8kCeTmxody6soskQskukgglu0gilOwiiVCyiyRCyS6SCCW7SCJ6q84ew2qIVV6k77/Ai7KZOT73b/nK0WCb1fmx6zn+npqd57GtO8T7AHgmXFCuH+fHnr6exz59bB1t/5+pIdpupfBQz/4zkWtNZJTpplf4czZ9TfjlbZGpnrPlyPDad6f4DiJKv7cl3Bi7BDfZH0VXdpFEKNlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXSURP1dljyyqz2Xtjyx7nZyN19M18SuQL14bn/h06xevk8yMtvqdGzkt5NDzuO3ZeCuf5mPG5DXzgd+40XzZ5LZkl2/mhseYsr3UXx/h0zdly+MTFllSuR6YHR5YfO4b1C4mdl6aP2ZndikivUbKLJELJLpIIJbtIIpTsIolQsoskQskukoieqrPHsLm+swu8GL2wjj/U0oZIzZbsf/oqvu3s1bwOP3g0UrON1NkrfEg5NfoGr2WXN/E56/PF2Nzu4eALRX7sC9dHlsI+yLdn9erKEI976g95/4FN2Ejb+04XaXtlbWt1+mZEr+xm9rCZTZrZ/iW3jZrZU2Z2sPF9fWfDFJFWreRj/E8B3Pqh2+4BsM/ddwDY1/hdRHpYNNnd/RkA5z508+0AHmn8/AiAO9ocl4i0WbP/oNvi7icBoPF9c+iOZrbHzCbMbKJW4n/HiEjndPy/8e6+193H3X08OzDY6cOJSECzyX7azMYAoPF9sn0hiUgnNJvsTwC4s/HznQAeb084ItIp0Tq7mT0K4BYAG83sGIDvA7gfwC/N7C4A7wH4WjuCia233XeBjE+ObFuO1DWLZI1zAKj1hY9d2cAHjecu8NM8P8oL6X2RMeebXwoPzq5n+baxdcRH3uCxz1zDY9/yXji2yRt5Db8eeXWe/TR/bL4j/D+i6jw/NmZiqcHr8Jkaby+Q17LHul2Qdta3IJrs7r470PTF2LYi0jvUXVYkEUp2kUQo2UUSoWQXSYSSXSQRl9QQVydvTWyaaSA+DLQyxHeQ2VYKH3uST0MdKynGSkx1XsVBaWO4FjNwhh98ahcvQf3N7n20/cGXbqbtF7f3Bdsqw7G5w3lzZRMvebJH5vXIzof4eZv9FH+9ZCPl1mp/+Pix1zLLA0ZXdpFEKNlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXSURP1dlj9eby+nBtcm4zr5vWI8MG11w7Tdtnzw0E27IVfuzKRr4+cO58ZAjsZl7zra4Jv2fXCpF67yCvdd+78S3a/q+n/oK2z5PlpPvO0k0xv4nHlh9c4Ds4HJ4ZydbxYnbfJH/BzG/lNX7PRmKfDZ+XavilBiBehw/RlV0kEUp2kUQo2UUSoWQXSYSSXSQRSnaRRCjZRRLRU3X22Pjlaj9pG+J1zcwC3/nC4WHa3keWJl4Y5YXPwml+mj3yLHieP7b+M6SWPc1jK35+jrbf+uZXaHt2PnJe14djX/su3RTVgchy0Pv5JAXr3w4/djYHAADMXcbPeXaGb2+R8fK18DD/aH+TZunKLpIIJbtIIpTsIolQsoskQskukgglu0gilOwiieitOnsMKV1m53hdM8OHH2PNFH/f6z8brtlOX8drrvlZfuwqn3Ye1aHma7bFrfxxjQyH58MHgKP/fSVtXxjjY+03PxeO/cxn+OMaPsRr3QNTvA/B+evDz4tFpqxnfRcAYH5DZP6EAj9AjfQZic4bH+mPEhK9spvZw2Y2aWb7l9x2n5kdN7OXG1+3NXd4EVktK/kY/1MAty5z+w/dfVfj68n2hiUi7RZNdnd/BsC5VYhFRDqolX/Q3W1mrzY+5q8P3cnM9pjZhJlN1ErFFg4nIq1oNtl/DOBaALsAnATwg9Ad3X2vu4+7+3h2IDwBoIh0VlPJ7u6n3b3m7nUAPwFwU3vDEpF2ayrZzWxsya9fBbA/dF8R6Q3ROruZPQrgFgAbzewYgO8DuMXMdgFwAEcAfLMt0cSW6yb1xwKf9h2ZKt+5Gy9ezm0Mvy8anxYeC2t5e3UgEluBF17LpObbf4a/n585w4PLr+GxrT3M+xjMbgu3DR/m+64MR+bjL/HHNnA6vH82nz0AzEXmrC9MR+rwG2kzrcPH5l5oVjTZ3X33Mjc/1IFYRKSD1F1WJBFKdpFEKNlFEqFkF0mEkl0kEb01xDVScciWw225uch0y+d5+Wp+hL/vVdaGg2PDFQEgx2drRp5MUw0ApSv4+NxajcQ2wx/X4Os8+FqBNkeHDjOxoZx95/lzGhtmOntF+ACZyGrPsTJwPc/bM5FyrEeWEO8EXdlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXSYSSXSQRvVVnjyFl1XjNlhc+c5HhkuduCBdWPRsZqhmZoGdhE5+OGcVIUTdHhnJeyQvKQwd4IX3ucn7e6gVeMGb1ZI9camLDSIs7+GOzYvjg1c28g0DmPD/ntcj03/VY/wQWeqTG32yNXld2kUQo2UUSoWQXSYSSXSQRSnaRRCjZRRKhZBdJxCVVZ2dTNsfGVdfzvGabm+O17rVHw++Lfef5e2bpssiSy0W+vVX49kM3nA+2zRb5ePXSLv6417zJC8rzW/j2hXPhonA+Ms6/tJXv20qRgjM7bWQOACDeb2NhA7+Dlfn+s/Os0wg/dmy56RBd2UUSoWQXSYSSXSQRSnaRRCjZRRKhZBdJhJJdJBGXVJ2djeO1Oi8+Vgd4TTbWzsbD52f5e2auzPc9uYk2o7YxMvY6E6751k/xOvvYzknafuK6yPWgyF9C5bFw7Atlvu9MpFYdHw8fvkM9zzeOjRnPlXhssXnlY7HzjZvbLHpIM9tuZk+b2QEze93Mvt24fdTMnjKzg43v65sLQURWw0reX6oAvufuNwD4EwDfMrOdAO4BsM/ddwDY1/hdRHpUNNnd/aS7v9j4+SKAAwC2ArgdwCONuz0C4I5OBSkirftYfzmY2VUAbgTwHIAt7n4SWHxDALA5sM0eM5sws4laqdhatCLStBUnu5kNAfgVgO+4+8xKt3P3ve4+7u7j2YHIzIsi0jErSnYzy2Mx0X/u7r9u3HzazMYa7WMA+L91RaSroqU3MzMADwE44O4PLGl6AsCdAO5vfH+8IxEuwcohlaHWhizGhg1e3NYXbOub5huXR3hsg8f5sWtTfF7i+cMbg235fh7b6Ve28IMP8BOXuxgpn1XCT1ps2eR8ibcj8pyy8ld1kD8nmcixY6+neiSznLTHynKxY4espM5+M4C/BfCamb3cuO1eLCb5L83sLgDvAfhacyGIyGqIJru7P4vwcPovtjccEekUdZcVSYSSXSQRSnaRRCjZRRKhZBdJxCU1xJWpFXjddGEt335hXWS6Z1Lq7pvm++4/xwujhVnePvT2BdpeXRee7rm8Idw/AOD1XgCYW8/Hetb4CFqUR8JthUg/zMoQbx99iy8nPXsZiz1SZ4+soh2LzSLbg7U3OYQ1Rld2kUQo2UUSoWQXSYSSXSQRSnaRRCjZRRKhZBdJxCVVZ2e17gwvucLmeXtpa5ODhBePTltzbHleAMUxvn09xyfuHXwvPPi6cJFPQ50/wwdu9w/zQnq2xAelz4+FC9KVwchS1ZF68/nr+cu3Rsazz2/mz7cP8RdU7iyfKzo3G5lfgdXZI0s2N0tXdpFEKNlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXScQlVWdn9cfYErnVNbx42T8ZaT8TLvp6hheEPcP3PR+e9h0AcOLzfPvhQ+HB+rFx1YWLkQHpEdPXDdD24cPhtnOf5vuurYvUuofKtL2vP9zHYH2e9z8oV/gLqliJ9K0oRVIrNt69A3RlF0mEkl0kEUp2kUQo2UUSoWQXSYSSXSQRSnaRRKxkffbtAH4G4DIsroi9191/ZGb3Afh7AFONu97r7k92KtDFYMJNdT69OSxSh4/N1b0wHD74mim+8cIw33e9L1Kn7+dF2dk/Dtebv3Dd23Tb4dwcbf+Pd3fS9p2bT/P2vzoVbPvToYN026MLG2j7sYVR2s72/+sz43Tbp9+6nrZni5HrZGx6hA6NWWdW0qmmCuB77v6ima0F8IKZPdVo+6G7/3PnwhORdlnJ+uwnAZxs/HzRzA4A2NrpwESkvT7W3+xmdhWAGwE817jpbjN71cweNrNl504ysz1mNmFmE7VSsaVgRaR5K052MxsC8CsA33H3GQA/BnAtgF1YvPL/YLnt3H2vu4+7+3h2YLANIYtIM1aU7GaWx2Ki/9zdfw0A7n7a3WvuXgfwEwA3dS5MEWlVNNnNzAA8BOCAuz+w5PaxJXf7KoD97Q9PRNrF3HnZx8w+B+B/AbyGDwoK9wLYjcWP8A7gCIBvNv6ZF9R/+Xa/6hvfbTHk5sSWJo4NU2VLE9fzsSGukWP38TpN7gKvKzpp9iyPLbvAa0Ajf3CGthef3UTba2vCx6+TqcEBoDLMS46Dl/H/AbGX9txJvuZyfjoyzXU1suQzH0HbMUcefADzJ44uG9xK/hv/LJavCna2pi4ibaUedCKJULKLJELJLpIIJbtIIpTsIolQsosk4tKaSroFFlvSOTLmMENXJm51vGJkyebI8FzWRyBb5vvOX+T7Lk7zOrpFhnLmiuFzkz3Ht+07y1+e9WPraHutP3xe+iJ18tgU3N0YotoqXdlFEqFkF0mEkl0kEUp2kUQo2UUSoWQXSYSSXSQR0fHsbT2Y2RSA3y65aSMAPmC6e3o1tl6NC1BszWpnbFe6+7KdI1Y12T9ycLMJd+cTeHdJr8bWq3EBiq1ZqxWbPsaLJELJLpKIbif73i4fn+nV2Ho1LkCxNWtVYuvq3+wisnq6fWUXkVWiZBdJRFeS3cxuNbO3zOwdM7unGzGEmNkRM3vNzF42s4kux/KwmU2a2f4lt42a2VNmdrDxfdk19roU231mdrxx7l42s9u6FNt2M3vazA6Y2etm9u3G7V09dySuVTlvq/43u5llAbwN4C8BHAPwPIDd7v7GqgYSYGZHAIy7e9c7YJjZnwOYBfAzd//9xm3/BOCcu9/feKNc7+7/0COx3QdgttvLeDdWKxpbusw4gDsA/B26eO5IXH+NVThv3biy3wTgHXc/7O4LAH4B4PYuxNHz3P0ZAB+ez+V2AI80fn4Eiy+WVReIrSe4+0l3f7Hx80UA7y8z3tVzR+JaFd1I9q0Aji75/Rh6a713B/AbM3vBzPZ0O5hlbHl/ma3G981djufDost4r6YPLTPeM+eumeXPW9WNZF9u9q5eqv/d7O6fBfBlAN9qfFyVlVnRMt6rZZllxntCs8uft6obyX4MwPYlv28DcKILcSzL3U80vk8CeAy9txT16fdX0G18n+xyPP+vl5bxXm6ZcfTAuevm8ufdSPbnAewws6vNrADg6wCe6EIcH2Fmg41/nMDMBgF8Cb23FPUTAO5s/HwngMe7GMvv6JVlvEPLjKPL567ry5+7+6p/AbgNi/+RPwTgH7sRQyCuawC80vh6vduxAXgUix/rKlj8RHQXgA0A9gE42Pg+2kOx/RsWl/Z+FYuJNdal2D6HxT8NXwXwcuPrtm6fOxLXqpw3dZcVSYR60IkkQskukgglu0gilOwiiVCyiyRCyS6SCCW7SCL+D38LiL89ZrQcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now we test the results\n",
    "succeed = (np.sum(mytrls * (mytri @ w) >= 0)) # if y_i and w.X_i have same sign then the classification is correct\n",
    "total = (mytrls.shape)\n",
    "\n",
    "print(' train success = ',np.round(100*succeed/total,2),'%\\n fail = ',np.round(100*(total-succeed)/total,2),'%')\n",
    "\n",
    "succeed = (np.sum(mytels * (mytei @ w) >= 0))\n",
    "total = (mytels.shape)\n",
    "\n",
    "print(' test success = ',np.round(100*succeed/total,2),'%\\n fail = ',np.round(100*(total-succeed)/total,2),'%')\n",
    "\n",
    "img = w.reshape((28,28))\n",
    "plt.imshow(img)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(17942, 784)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_imgs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}