{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "name": "Machine learning classification as a constrained quadratic optimization problem.ipynb",
   "provenance": [],
   "authorship_tag": "ABX9TyM4mfQpET0ZUlPleqa8G3Zk"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "source": [],
    "metadata": {
     "collapsed": false
    }
   }
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Xh0psy4RWn_J"
   },
   "source": [
    "# Machine learning classification as a constrained quadratic optimization problem\n",
    "\n",
    "## Try me\n",
    " [![Open In Colab](../../_static/colabs_badge.png)](https://colab.research.google.com/github/ffraile/operations-research-notebooks/blob/main/docs/source/NLP/tutorials/Machine_learning_classification_as_a_constrained_quadratic_optimization_problem.ipynb)[![Binder](../../_static/binder_badge.png)](https://mybinder.org/v2/gh/ffraile/operations-research-notebooks/main?labpath=docs%2Fsource%2FNLP%2Ftutorials%2FMachine_learning_classification_as_a_constrained_quadratic_optimization_problem.ipynb)\n",
    "\n",
    "## Problem definition\n",
    "In this notebook we are to provide a closed form to reduce a machine learning linear classification problem to a constrained quadratic optimization problem.\n",
    "\n",
    "In this example, let us imagine that we want to predict the status of an engine by reading sensor data. We do not know what is the mathematical expression that models the status of the machine as a function of the sensor data reading, but we have collected data in different conditions, so we have a **dataset** with the data readings that have been generated by sensors measuring the engine in different states: normal operation, and failure state.\n",
    "\n",
    "In linear classification problems, we find the following definitions: \n",
    "\n",
    "- $x^{(t)}$: feature vector of length t of independent variables (e.g. in our example this is going to be the different sensor data readings like temperature, rotation speed, torque) x is a vector of length d i.e. $x^{(t)} \\in \\mathbb{R}^d$, and each element $x_1, x_2, ...$ is going to represent the readings of the different sensors.\n",
    "\n",
    "- $y^{(t)}$: Label t representation the variable that we want to explain from the feature vector (e.g. in our case, whether the engine is in a failure state or not)\n",
    "\n",
    "The objective is to learn the parameters of a **maximum margin linear separator** from a training dataset. A maximum margin linear separator is a [hyperplane](https://en.wikipedia.org/wiki/Hyperplane) in the vector space X. Think of a hyperplane as a plane that divides the vector space into two regions. For instance, a line, which is a 1 dimensional plane, divides a 2-dimensional space into two regions, the area above the line and the area below the line. Consider that a hyperplane in the vector space of size $d$ can be represented by a vector of size $d$. Let us note as $\\hat{\\theta} \\in \\mathbb{R}^{d}$ the vector that represents our linear separator.\n",
    "\n",
    "To clearly visualize the results, we are going to work with a linear separator of size $d=2$, so that we can visualize the results in a coordinate system. Therefore, we are only going to have data from two sensors, and our data points are noted as:\n",
    "\n",
    "$x = [x_1, x_2]$\n",
    "\n",
    "$x_1$: temperature sensor data readings \n",
    "\n",
    "$x_2$: Encoder measuring the engine speed in rpms \n",
    "\n",
    "## Dataset generation\n",
    "Now, since we do not really have a dataset, we are going to generate a synthetic one. The following script generates random points and represents them in the bi-dimensional containing all data points.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "kZqxNmsyc4Ul"
   },
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "## This parameter determines (half) the number of samples used\n",
    "no_samples = 30\n",
    "\n",
    "## This array represents our temperature readings in failure state\n",
    "temp_readings_failure = 7*np.random.random_sample((no_samples,1))+29\n",
    "\n",
    "## This array represents our temperature readings in normal state\n",
    "temp_readings_normal = 7*np.random.random_sample((no_samples,1))+22\n",
    "\n",
    "#This array represents our rpms in failure state\n",
    "rpms_failure = 400*np.random.random_sample((no_samples,1))+2900\n",
    "\n",
    "#This array represents our rpms in normal state\n",
    "rpms_normal = 400*np.random.random_sample((no_samples,1))+2500"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Problem representation\n",
    "Now, let us plot our dataset, together with a hyper-plane which is not necessarily a separator, to illustrate how a line (1-dimensional hyper-plane) divides the 2-dimensional space into two regions:"
   ],
   "metadata": {
    "id": "OSYKXDTADnoo"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "#Prepare the figure \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "#Plot the data in failure state in red\n",
    "ax.scatter(temp_readings_failure, rpms_failure, color='red')\n",
    "#Plot the data in normal operation in blue\n",
    "ax.scatter(temp_readings_normal, rpms_normal, color='blue')\n",
    "\n",
    "# Make a linear space to plot the temperature\n",
    "t = np.linspace(22,37)\n",
    "\n",
    "# Initial separator\n",
    "thetahat_1 =-57\n",
    "thetahat_0 = 4452\n",
    "\n",
    "# \n",
    "r = thetahat_1*t + thetahat_0\n",
    "ax.plot(t,r)\n",
    "plt.xlabel('$x_1$ (temperature in celsius degrees)')\n",
    "plt.ylabel('$x_2$ (engine speed in rpms)')\n",
    "\n"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 300
    },
    "id": "6bvtUmvKDnNs",
    "outputId": "723685b9-4986-4389-d845-09e3cd6521d1"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$x_2$ (engine speed in rpms)')"
      ]
     },
     "metadata": {},
     "execution_count": 27
    },
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUddbA8e9JSIAAUiMgJfQmFooIChJAERFFXV0Lort2Ze2sq2Iv6+4KunbF1w7YsYCIIBCKSpfepYogvYaW5Lx/3Jt1xEySm8zMnXI+z3Ofmblzy5kLmTP3V0VVMcYYY7xI8jsAY4wxsceShzHGGM8seRhjjPHMkocxxhjPLHkYY4zxrIzfAURCjRo1tEGDBn6HYYwxMWXOnDnbVDW9oPcSInk0aNCA2bNn+x2GMcbEFBFZF+w9K7YyxhjjmSUPY4wxnlnyMMYY45klD2OMMZ5Z8jDGGOOZJQ9jTGgMHw4NGkBSkvM4fLjfEZkwSoimusaYMBs+HG64AbKzndfr1jmvAfr18y8uEzZ252GMKb1Bg35LHPmys531Ji5Z8jDGlN769d7Wm5hnycMYU3r163tbb2KeJQ9jEkG4K7OffBLS0n6/Li3NWW/ikiUPY+JdfmX2unWg+ltldigTSL9+MHQoZGSAiPM4dKhVlscxSYQ5zNu3b682MKJJWA0aOAnjaBkZsHZtpKMxMURE5qhq+4LeszsPY+KdVWabMLDkYUy8s8psEwaWPIyJd9FYmW290WOeJQ9j4l20VWZHogLfhJ1VmBtjIssq8GOGVZgbY6KHVeDHBUsexpjIsgr8uGDJwxgTWdFYgW88s+RhjImsaKvANyVi83kYYyKvXz9LFjHO7jyMMcZ4ZsnDGGOMZxFLHiJSTkRmish8EVksIo+664eLyHIRWSQib4pIirteROR5EVklIgtEpG3Asa4WkZXucnWkPoMxxhhHJO88DgHdVfUk4GSgl4h0BIYDLYATgPLAde725wBN3eUG4BUAEakGPAycCnQAHhaRqhH8HMYYk/AiljzUsc99meIuqqpj3PcUmAnUdbfpC7zrvjUdqCIitYGzgfGqukNVdwLjgV7hinvVln0kQi98Y4zxIqJ1HiKSLCLzgC04CWBGwHspQH9grLuqDrAhYPef3XXB1h99rhtEZLaIzN66dWuJ4l2/PZvez03lqjdnsvLXvSU6hjHGxKOIJg9VzVXVk3HuLjqISOuAt18Gpqjq1BCda6iqtlfV9unp6SU6Ru0q5bivdwvmb9hFr+em8uioxew+cCQU4RljTEzzpbWVqu4CJuEWN4nIw0A6cFfAZhuBegGv67rrgq0PuZTkJP56ekMmDczk0lPq8fb3a+k2OIsRM9aTm2dFWcaYxBXJ1lbpIlLFfV4eOAtYJiLX4dRjXK6qeQG7fAlc5ba66gjsVtVNwDdATxGp6laU93TXhU31imX554UnMPrWzjRJr8j9ny3k/BenMWvtjnCe1hhjolYk7zxqA5NEZAEwC6fOYzTwKlAT+EFE5onIQ+72Y4DVwCrgdeAWAFXdATzuHmMW8Ji7LuyOP64yH97YkRcub8OO/Ye55NUfuO39H9m0+0AkTm+MMVHD5vMooezDObya9ROvTVlNkgi3ZDbm+jMaUS4lOaTnMcYYv9h8HmGQllqGu3o259u7upLZPJ0h41dw5jOTGbtokzXtNaVjU7SaGGDJo5TqVUvjlSvbMeK6U6mQWoabhs2l3//NYPlma9prSsCmaDUxwoqtQignN48RM9czZNwK9h3K4cpT63PnWc2okpYa9nObOGFTtJooUlixlSWPMNi5/zDPjF/B8BnrqFw+hbt6NueKDvVJTpKIxWBiVFKSc8dxNBHIy/vjemPCyOo8IqxqhVQev6A1X93Whea1KvHg54vo88I0pq/e7ndoJtrZFK2xI8Hrpix5hFHL2sfw/vUdeblfW/YcOMJlQ6czYMRcNu6ypr0JxcuXTGmmaE3wL7OIsropUNW4X9q1a6d+yz6Uo8+OX67NHxijzR8Yo8+OX67Zh3L8DsuE27Bhqmlpqs5XjLOkpTnrC9snI0NVxHksbNvSnMeUXEbG7691/pKR4XdkIQXM1iDfq1bnEWEbdx3gn2OW8tWCTdSpUp77e7ek9wm1ELH6kLgUqQpwq2iPrASpm7I6jyhSp0p5XrqiLR/c0JFjyqcwYMRcLn99Oks37fE7NBMO69d7Wx/t5zEOq5uy5OGXjo2qM/rWzjxxQWuWbd7Luc9P5YHPF7Jz/2G/QzOhFKkvGfsyi6zS1E3FCc/JQ0QqiIiNwRECyUnClR0zyBqYSf+OGbw/cwOZg7N494e15OTGz61vQovUl4x9mUVWv34wdKhTLCjiPA4d6qxPEEXWeYhIEnAZ0A84BWc62bLANuAr4DVVXRXmOEslmuo8CrN8814eHbWY73/aTvOalXj4vFac1qSG32GZ0ho+HAYNcoqQ6td3vtDD8SUTqfOYhFGqToIiMhn4FvgCWKTusOnuXOLdgCuAz1R1WEijDqFYSR7gtH77ZvFmnvhqKT/vPMA5rWtxf++W1KuWVvTOxhgTQqWtMD9TVR9X1QUaMN+GOnOIf6qqfwI+DFWwiU5E6NW6Nt/e1ZW7z2pG1vKtnPnMZJ4Zt5zswzl+h2dCKZb6ZcRSrCYiikweqnoEQEQuEZFK7vMHRWSkiLQN3MaETrmUZG7t0ZSJA7ty9vG1eH7iKnoMmcyX83+xUXvjQSx1MoulWE3EeKkwf1BV94pIZ6AH8AbwSnjCMvlqVy7P85e34eObOlGtQiq3vf8jf37tBxZt3O13aKY0Bg2C7Ozfr8vOdtZHm1iK1S8JeGdW7E6CIvKjqrYRkaeAhao6In9deEMsvViq8yhMbp7y0ewNPP3NcnZmH+ayU+ozsGczqlcs63doxqtY6mQWS7H6If/OLDDBpqXFReurUHUS3CgirwGXAmNEpKzH/U0pJScJl3eoz6S7M/nLaQ34aPYGug3O4s1pazhiTXtjSyz1y4ilWP2QoHdmXr78/wx8A5ytqruAasDfwxKVKVTltBQePu94xt7ehZPqVeGx0Uvo/dxUpq7c6ndopiAFFWnEUr+MWIrVDwnau7/YyUNVs1V1pKqudF9vUtVx4QvNFKVpzUq8e00HhvZvx6GcPPq/MZPr353N+u3ZRe9sIiNYZTPETicz6xBXuAS9M/NS59EeGARkAGUAAVRVTwxfeKERL3UehTmUk8sb09bw4sRV5OQq15/RkFsym1ChbBm/Q0tsNmBh/LM6jyINB94C/gScB/RxH00UKFsmmVsymzBpYCbnnliblyb9RPchWXz+40Zr2uunBC3SiCmlbSmVoHdmXu48pqlq5zDHExaJcOdxtDnrdvLoqMUs+Hk37TKq8sh5x3NC3cp+h5V47M4jusXxXUMohGQOcxHpAVwOTMAZ3woAVR0ZiiDDKRGTB0BenvLJnJ/5zzfL2L7/MH9uV4+/92pODWvaGzn25RTdLLkXKlTJYxjQAlgM5LcLVVW9JiRRhlGiJo98ew4e4YUJK3nru7WUT0nm9jObclWnBqSWsZbWEWEDFkYv68NSqFAlj+Wq2jykkUVIoiePfD9t3cdjo5YwecVWGqVX4KE+rchsfqzfYcUnSxixwe48ChWqCvPvRaRViGIyPmicXpG3/3oKb1zdnrw85S9vzeLat2exdtt+v0OLLzYWlP+KWwlufVhKzMudx1KgCbAap87DmurGsEM5ubz93Vqen7CSI7nKNZ0b8rfuTahoTXtLz37N+strPZPdJQYVqmKrjILWq2oBfyXRxZJHcFv2HOTfY5fz6dyfObZSWf7RqwUXtqlDUpL4HVrssnJ0f1nyDplQFVv9itPH41ngGeAid52JYcceU44hfz6Jz245jdpVynP3x/P506vfM3/DLr9Di10J2uM4akRT35o4Hm3XS/J4FzgeeAF4EWgFvBeOoEzktalflc9uPo3Bl5zEzzsP0Pel7/j7x/PZsveg36HFHitH91e0JO84r/vykjxaq+q1qjrJXa7HSSYmTiQlCRe3q8vEu7tyY9dGfD5vI90HT2bolJ84nGPFLcWWoD2OS6w4v869/IIvKHkD7NsX2S/ueB9tV1WLtQDDgI4Br08F3i3u/n4u7dq1U+Pd6q379Jq3ZmrGP0Zrt6cn6cSlv/odkok3w4appqWpOr/NnSUtzVnvZZuCjlu9+u/3Kc5+oSTyx/PnLyKqGRmRi6WEgNka5HvVa2ur5kB+wWF9YDmQQzFaXYlIOWAKUBZnYMVPVPVhEWkIfABUB+YA/VX1sDtfyLtAO2A7cKmqrnWPdR9wLZAL3Kaq3xR2bqswL51Jy7fw+KglrN62n27N03mwTysapVf0OywTD4pTuV3SCnC/K86DnT9QlI82ENbWVvm0iFZXIiJABVXdJyIpwDTgduAuYKSqfiAirwLzVfUVEbkFOFFVbxKRy4ALVfVSt6/J+0AH4DjgW6CZquYGO7clj9I7nJPHO987TXsP5uRyzelO095K5VL8Ds3EsuK0TCtp6zW/W70V1GS4IFHcCqzUra3cL/48VV0XbCnqGO5d0D73ZYq7KNAd+MRd/w5wgfu8r/sa9/0ebhx9gQ9U9ZCqrgFW4SQSE0apZZK4/oxGTByYyYVt6jB06mq6DZ7MR7M3kJdno/aaEipO5XZJK8D9rjg/uu4rmBgdYblYycMt+xpT2pOJSLKIzAO2AOOBn4BdqprjbvIzUMd9XgfY4J4/B9iNU7T1v/UF7BN4rhtEZLaIzN661WbYC5X0SmX5z8Un8cWA06lfrTz3fLKAC1/+jrnrd/odmolFxWmZVtLWa9HQ6q1fP+euIi/PSSIFidEm3F5aW80VkVNKczJVzVXVk4G6OHcLLUpzvCLONVRV26tq+/T09HCdJmGdWLcKn9x0Gs9eehKbdh/kope/564P5/HrHmvaazwoTsu0krZeK2i/q692WjsVp9VWqPtoREMyC6VgNelHL8AynMrxn4AFwEJgQXH3L+B4D+HMgb4NKOOu6wR84z7/BujkPi/jbifAfcB9Acf533bBFmttFV57Dx7Rf329VJveP0ZbPfi1vjRppR48kuN3WMb8npdWWyVp4VXcGDIyft/aqqB1UYIQtbYq1fAkIpIOHFHVXSJSHhgH/Bu4GvhUf6swX6CqL4vIAOAE/a3C/CJV/bOIHA+M4LcK8wlAU7UKc9+t3bafJ75ayrdLfyWjehoPntuKHi2PRQor7zUmUry0vopUS60on+8lJK2tQhDEiTgV4Mk4xWUfqepjItIIp6luNeBH4EpVPeQ27X0PaAPsAC5T1dXusQYB1+DcCd2hql8Xdm5LHpE1ZcVWHh21mJ+27ueMZuk81KclTY6t5HdYJtF5aX0VqZZafjcnLkJUJA8/WfKIvCO5ebz7wzr+++0KDhzO5erTGnD7mU05xpr2Gr9E452H382JixCqgRGNKbaU5CSu7dyQSQMzuaR9Xd78bg3dns7ig5nrybWmvcYPXiqsI1W57Xdz4lIodvIQkbIicoWI3C8iD+Uv4QzOxL4aFcvy1EUnMupvnWlYowL3jlxI35emMXvtDr9DM4nGS6utSI1PFsMtsLzceXyB00EvB9gfsBhTpNZ1KvPxTZ147rKT2bb3MBe/+gN3fPAjm3db015fxPFQ4YUK7Hexdu1vyaCg6xFs21DHE6uDaAZrhnX0Aiwq7rbRtlhT3eiy7+ARfXrsMm06aIy2fPBrfXHiSj1w2Jr2Rky4mqH6pbRNXYu6HlHclDbcKKSprpfkMRSn6azvycDrYskjOq3btl+vf2eWZvxjtHb590Qdu2iT5uXl+R1W/MvI+P0XZf6SkeF3ZN4F++K/+ebif+EXdj3iLdF6VFjy8NLPYwnOHOZrsDnMTQhNW7mNx0YvZsWv++jcpAYPn9eKpjWtaW/YRHkLH0+CtYoS+f1nLKzvRLDrUZgoaUobbjaHuSWPqJeTm8ew6et4ZvwK9h/OpX/HDO48sxmV06xpb8hFed8CT7x88Qf7fMUZOv1osZhoSyAkTXW1hKPpGlMcZZKT+MvpDcn6ezcuPaUe7/ywlm5Dshgxw5r2hlwMt/D5g2rVir9tsNFrg808WJgYaEobbkUmDxGZ5j7uFZE9ActeEdkT/hBNIqlWIZV/XngCo2/tTJP0itz/2ULOf3EaM9dY096QieUWPoGGD4c9Hr6Cgn3hF3fo9HyxmmhDzHqYm6ilqoxesImnxizll90HOe+k47jvnBYcV6W836GZaBCsuKlCBacoq6TjRRVWjJWR4SSOWEu0JWQ9zE1MEhHOO+k4JtydyW09mjJu8Wa6D8lyZjM8EnQczOiRqH0pIiVYMVR2dunurIIV6w0bFr7+HjHI7jxMzNiwI5unvl7KmIWbqVu1PIN6t6RX61rROWpvlI+WGhfCWfE/fLgz78f69U5xVwLdbQSyOw8TF+pVS+Plfu0Ycf2pVCxbhpuHz+WK12ewbHPh5d6+3AAMGvTHuauzs531JjTCWfEfid7lMa5MURuIyF2Fva+qz4QuHGOKtmZ6DZY915kdVdfzQ9cVnLN6Kv07ZXDXWc2okpb6u22PvgFYt855DWH+PghWpBKj81VHpfx/QLtD8EWRxVYi8rD7tDlwCvCl+/o8YKaqXhm+8ELDiq3ix9HJIKncYWp0W0HaieuonJbC3T2bc0WH+iQnOUVZvnVpiKe+FCZhhaqT4BTgXFXd676uBHylqmeELNIwseQRP4J+J5+4h44DFjN99Q5a1KrEI+cfT8dG1f3rTG11HiYOhKrOoyZwOOD1YXddwoq2xjTRFk9xeI05aGnQwmN4//qOvNyvLXsP5nDZ0OkMGD6X+i2yC9w+7H284qUvhTHBBBv06ugFGATMBx5xl3nA/cXd388lHAMjRtt4adEWT3GUJObijOl34HCO/nf8Cm3+wBhtfO8YrZG5XKVMTsxcF2OiBaEYGBFARNoCXdyXU1T1xxDmsbAJR7FVtBVpR1s8xVGSmL2UBm3cdYCnxixl9IJNsL88275tQY0DtXnySbEbAGOKISTFVuI0pm8FVFbV54DtItIhRDHGnOI0polkMVIsNu4pScxeSoPqVCnPi1e05cMbOtKySQo1+v7IqfdNp003G1XHmNLyUmH+CpAHdFfVliJSFRinqqeEM8BQ8OPOI9L1pYly51FSuXnK+zPXM2TccnYfOMLlHepzd8/mVKuQWvTOxiSoUFWYn6qqA4CDAKq6E0jYv7yi+idFuo9YLA6UGsmYk5OEKztmkDWwG1d1asAHszbQbXAWb3+3hpzc+B9a25hQ85I8johIMqAAIpKOcyeSkIoqPol0MVIsNu7xI+bKaSk8cv7xfH17F1rXOYZHRi3h3Oen8f2qbeE7qTFxyEuxVT/gUqAt8A5wMfCAqn4cvvBCw49+HrFYjJRoVJVvFv/Kk2OWsGHHAXodX4tB57akXjWPczsYE6dCNRnUcOAe4ClgE3BBLCQOv0SySKagivlY7PMRaSJCr9a1GH9nVwb2bMbkFVvp8cxkhoxbTvbhHL/DMyaq2ai6YRSJgTkLqphPTXV6NBw58ts669xctE27D/Cvr5fxxbxfqF25HPee04LzTzouOkftNSYCQjU8iQD9gEaq+piI1AdqqerM0IUaHvE8PImX6ZejucgsmkbAnr12B4+MWsyijXs4pUFVHj7veFrXqexPMMb4KFTJw5rqRqFgYzcVJOzjOZVQNA4DlZunfDx7A09/s5wd2Ye57JR6DOzZnOoVy/oTkDE+sKa6McZLfYWXMZrCPp5TCUXj1BfJScJlHeozcWAmfz2tIR/P/pnMwVm8MW0NR6xprzHWVDfa5P8KX7fOuaPIn38iWAIpqGI+NRVSUn6/Lpr7fERz7/jK5VN46LxWjL2jCyfXq8Ljo5dwznNTmbpyq9+hGeMrL8njeeAzoKaIPAlMA/4ZlqgSmNdf4QX1lXjzTXjrrdjp8xHsjiia7pSaHFuJd6/pwOtXtedIbh7935jJ9e/OZv32gkftNSbeeR0YsQXQw305UVWXhiWqEIulOg/f5p/wUTTWeRTmUE4ub0xbw4sTV5GTq1x/RkNuyWxChbJFTsxpTEwJ1cCI5YDewJlAd6CXu664+9cTkUkiskREFovI7e76k0VkuojME5HZ+YMtiuN5EVklIgvcEX3zj3W1iKx0l6uLG0MsiIVf4aEWa73jy5ZJ5pbMJkwamEmfE2vz0qSf6D4ki89/3EgiNH03BvA0n8dHwBtAN3d5HfjYw/61gbbu80rACpxRescB57jrewNZAc+/BgToCMxw11cDVruPVd3nVQs7dzjm8wiXaJ2XY9gwZ84MEeexuPGUdL9YMmfdDj3vhama8Y/RetHL3+mCDbv8DsmYkKCQ+Ty8JI8lxVnn4XhfAGcB3wCXuusuB0a4z18DLg/YfrmbgC4HXgtY/7vtClpiKXmoRt8XbkkTWrQmwnDIzc3TD2et13aPj9cG947Wez6er1v3HvQ7LGNKpbDk4aWfxzDgRVWd7r4+FRigqlcV6wC/P1YDYArQGqjjJhDBKUY7TVXXicho4F+qOs3dZwLwDyATKKeqT7jrHwQOqOrgo85xA3ADQP369dutK25POvMHJR2nKxHH99p78AgvTFzFm9PWUD4lmdvPbMpVnRqQWsZL2xRjokOo+nm0A74XkbUishb4AThFRBaKyAIPwVQEPgXuUNU9wM3AnapaD7gTp2is1FR1qKq2V9X26enpoThkwippU9poboIbLpXKpXB/75Z8c+cZtGtQlSe+Wkqv56aQtXyL36EZE1JekkcvoCHQ1V0auuv6AOcV5wAikoKTOIar6kh39dVA/vOPgfzZCTcC9QJ2r+uuC7behElJK/ETsfI/X+P0irz91w68+Zf2qMJf3prFde/MYu22/X6HZkxIeEkeHYAdqroO6A88C1RX1XXuukK5Y2O9ASxV1WcC3voFJxmB04prpfv8S+Aqt9VVR2C3qm7CKeLqKSJV3SFSerrrTJiUdITgWJygKtS6t6jJN3ecwX3ntOCHn7Zz1rOTeerrpew7ZKP2mhgXrDLk6AVY4D52BrKAc3FbQBVz/844vdMXAPPcpbe7fg4wH5gBtHO3F+Al4CdgIdA+4FjXAKvc5a9FnTvWKsyjkbW2Kr1fdx/Quz6cpxn/GK2nPDFeP5m9QXNz8/wOy5igCFGF+Y+q2kZEngIWquqI/HWlS1/hF0udBE38+3H9Th4ZtYT5G3Zxcr0qPHr+8ZxUr4rfYRnzB6GqMN8oIq/hzCY4RkTKetw/bthES6Y02tSvymc3n8bgS05i464D9H3pOwZ+PJ8tew/6HZoxxeblziMNp4J8oaquFJHawAmqOi6cAYZCKO88Ym0oDRPd9h48wouTnKa9Zcskc2v3Jvz19IbWtNdEhVBNQ5utqiNVdaX7elMsJI5Q82v4cLvb+aN4uCaVyqVw3zktGXdnVzo0rMZTXy/j7P9OYdIya9proptNQ+uRHwMX2t3OH8XrNZm0fAuPj1rC6m376dY8nQf7tKJRekW/wzIJKlR1HoaS910oza/kaJwsyW/xek26NT+WsXecwaDeLZm1didn/3cK/xyzlL0HjxS9szERZMnDo5L0XfA6wdPRErGndlHi+Zqklkni+jMaMWlgJhe2qcPrU1fTbfBkPpq9gby8+C8pMLHBy5DsIiJXishD7uv6+cOnJ5KSDB9e2l/JidxTO5hEuCbplcryn4tP4osBp1O/Wnnu+WQBF778HXPX7/Q7NGM83Xm8DHTCGdUWYC9OJ76E06+fM7BfXp7zWFQZe2l/JVtP7T9KpGtyYt0qfHLTaTx76Uls2n2Qi17+nrs+nMeve6xpr/GPl+RxqqoOAA4CqOpOIDUsUcWZ0v5KjrXJkiIh0a5JUpJwYZu6TBqYyS2ZjRm9YBPdBmfxctYqDuXk+h2eSUBe+nnMAE4DZqlqWxFJB8ZZD/OixWvLIOOfddv388RXSxm/5FcyqqfxwLmtOLPlsThDyBkTGqFqbfU88BlQU0SeBKYB/wxBfHEv0X4lm/DLqF6B169qz7vXdCAlOYnr353NVW/OZNWWvX6HZhKEp34eItIC6OG+nKiqS8MSVYj5fedhTDgdyc3j3R/W8d9vV3DgcC5XdWrA7Wc2pXL5FL9DMzEuJHce7lhWbYHKQHXgkvyWV8YY/6QkJ3Ft54ZMGpjJJe3r8tb3a+g+OIv3Z64n15r2mjDxUmz1BdAXyAH2ByzG+CIehicJpRoVy/LURScy6m+daVijAveNXEjfl6Yxe+0Ov0MzcchLhfkiVW0d5njCwoqt4o81QiicqvLl/F94aswyNu85SN+Tj+Pec1pQu3J5v0MzMaSwYisvyWMo8IKqLgxlcJFgySP+NGjg9NQ/WkaG0/fGOLIP5/BK1k+8NmU1ySIM6NaY67o0olxKst+hmRgQqtZWnYE5IrJcRBaIyEIRWRCaEGOXFZ38JpLXIp6HJwmltNQy3N2zORPu6krXZukMHreCs56dzNhFm0mEQVFN+HhJHucATXHmDD8P6OM+JqzSjlkVTyJ9LcI5PEl+EhSBMmWcx1j/YVCvWhqv9m/H8OtOpXxKMjcNm8OVb8xgxa/WtNeUjA3JXgpWdPKbSF+LcNV5FHTcUB4/GuTk5jF8xnqGjFvO/sO59O+YwZ1nNqNymjXtNb9XqjoPEZmmqp1FZC+ggAQ+quoxoQ441MKVPPyY2yNa+TXPyaBBTlFV/frOuFal/WIPlgTzxdMPgx37DzNk3HLen7meyuVTGHh2cy47pT7JSdZL3ThCUmEey+zOI/zi5VoES4L54vGHwZJf9vDIqMXMXLODVrWP4ZHzj6dDw2p+h2WiQKg6Cd5VwHKtiJwculBjSyKN7FqUeLkWRdWZxNOQ7/laHXcMH97QkRevaMOu7MP8+bUf+NuIufyy64DfoZko5qXCvD1wE1DHXW4EegGvi8g9YYgt6tmYVb+Jl2tRUBLMF4vJsLhEhD4nHseEuzO5vUdTxi/5le5Dsnju25UcPGKj9po/8tLPYwrQW1X3ua8rAl/hJJA5qtoqbFGWkvXzMF7k16WsWwfJyZCb6yTDUNSpxIqfd6jEfgIAABljSURBVGbz1JhlfLVwE3WqlGfQuS05p3UtG7U3wYSqk+Ay4ARVPeK+LgvMV9UWIvJjNA/NbsnDmJL54aftPDpqMcs276VTo+o8fH4rWtSK+jYyJkRC1UlwODBDRB4WkUeA74ERIlIBWFL6MI0x0aZT4+qMvrUzj1/QmqWb99D7uak89MUidmUf9js047NiJw9VfRy4AdgF7ABuVNXHVHW/qibIzXzisJ7zJl+Z5CT6d8wga2AmV3bMYNj0dWQOzuK9H9aSkxtnTc9MsXkdkr0ZUAGoAvS2IdlLLxq/pK3nvClIlbRUHuvbmjG3d6FlrWN48IvF9HlhGt//tM3v0IwPvNR5jAV2A3OA/zW/UNUh4QktdKK1ziNaR4aNlz4bJnxUlbGLNvPEV0vZuOsAvU+oxf29W1K3apCmaiYmharC3IZkD7Fo/ZK2nvPehKOne6w4eCSXoVNW83LWKlThxq6NublrY8qn2qi98SBUFebfi8gJIYrJEL0jw4Zz0MF4k+hFfOVSkrmtR1Mm3J3JWa1q8vyElfQYksXoBb/YqL1xzuuQ7HNtSPbQidYv6XjpLR4Jgwb9cRDF7GxnfSKpU6U8L17Rlg9v6EjltFT+NuJHLh06ncW/7PY7NBMmXodkb0IJh2QXkXoiMklElojIYhG5PeC9W0Vkmbv+PwHr7xORVW7COjtgfS933SoRudfDZ4gq0folHS+9xSMhWu8e/XJqI6dp75MXtmblr3s574Vp3P/ZQnbst6a98cZLnYcA/YBGqvqYiNQHaqnqzGLuXxuorapzRaQSTsX7BUBNYBBwrqoeEpFjVXWLiLQC3gc6AMcB3+K09gJYAZwF/AzMAi5X1aB9TaK1zgMSu7w8HkRrvVU02J19hGe/XcF709dRITWZO89qxpUdM0hJ9vKb1fgpVHUeLwOdgMvd13uBl4q7s6puUtW57vO9wFKcMbJuBv6lqofc97a4u/QFPlDVQ6q6BliFk0g6AKtUdbWqHgY+cLeNSf36OV8yeXnOoyWOyCtNc+lovXuMBpXTUnjk/OP5+vYunFi3Co+OWkLv56YybaU17Y0HXpLHqao6ADgIoKo7gdSSnFREGgBtgBk4dxNdRGSGiEwWkVPczeoAGwJ2+5nfBmUsaL0xnhVU4X3NNVCjRvGSiRXxFa1ZzUq8d20HXuvfjoM5uVz5xgxufG82G3YUMOOWiRllPGx7RESScSaCQkTSAc8NN90BFT8F7lDVPSJSBqgGdAROAT4SkUZej1vAeW7A6RFPfb9roE3UKqjC+/Bh2L7deZ7fegqCJ4R+/SxZFEVEOPv4WnRtls4b09bw4sRV9Fg+mRvPaMTNmY1JS/XyVWSigZc7j+eBz4BjReRJYBrwTy8nE5EUnMQxXFVHuqt/BkaqYyZOQqoBbATqBexe110XbP3vqOpQVW2vqu3T09O9hGkSSHEqthOx9VS4lEtJZkC3Jkwc2JVzWtfihYmr6DFkMl/M22hNe2OMl7GthgP3AE8Bm4ALVPXj4u7vVri/ASxV1WcC3voc6OZu0wynKGwb8CVwmYiUFZGGQFNgJk4FeVMRaSgiqcBl7rbGeFbcm9JEbT0VLrUrl+e5y9rwyU2dqF4xlds/mMefX/uBRRutaW+sKDJ5SMAA/qq6TFVfUtUXVXVpQdsU4nSgP9BdROa5S2/gTaCRiCzCqfy+2r0LWQx8hDNi71hggKrmqmoO8DfgG5xK94/cbY3xrLDJnwJZyWd4tG9QjS8GdOZfF53A6q37Oe/Fadw3cgHb9x3yOzRThCKb6opIFk5R0xequj5gfSpOx8GrgUmq+nb4wiydaG6qa/wX2Fy6WjXYsweOHPnt/WgYbywR7D5whOcnrOSd79dSPjWZO85sxlWdrGmvn0rbVLcXzkCI74vIL24nvzXASpxmu/+N5sRhEltxmuEGNpfetg3eestaT/mhcvkUHuzTirF3dKFN/ao8PnoJ5zw3lSkrtvodmilAsTsJwv8qvGsAB1R1V9iiCjG780hMkRq12Dp6hp6qMmHpFh7/agnrtmdzZsuaPNinJRnVK/gdWkIJyai6scySR2KKRO/vaB1WP14cysnlzWlreXHiSo7kKtd1aciAbk2oUNaa9kaCJQ9LHgkpEkPL2/AkkfHrnoP8e+wyRs7dSM1jynLvOS3oe1IdkpKK01bHlFSohicxJqZEYtRiGxgxMmoeU45n/nwyI285jVrHlOPOD+dz8avfs+DnmCk9jzuWPEzcisS4U9E6rH68alu/Kp/dcjpPX3wi63ccoO9L33HPJ/PZutea9kaa5+QhImeJyOsicrL7+obQh2VM6UVi3CkbGDHykpKES9rXY9LArlzfpRGf/biR7oOzeH3Kag7n2FSXkeK5zkNE3scZCfcBYAxwsareEobYQsbqPEw4WWsrf63euo/HRy9h0vKtNEqvwIN9WtGt+bF+hxUXQl3nsVdVd6nqQJyJoU4pagdjYl1h/UVsWH1/NUqvyFt/7cBbfzkFVfjrW7O45u1ZrNm23+/Q4lpJksdX+U9U9V7g3dCFY0z0SfR5ymNFtxbH8s0dZ3B/7xbMXLODns9O5qkxS9l78EjROxvPvMwk+BzOMOox17bXiq1MaVhz3NizZe9Bnh67nI/n/Ex6pbL8o1cLLmpjTXu9ClWx1V7gSxFJcw96toh8F4oAjYlm1hw39hxbqRxPX3ISnw84nTpVyjPw4/lc+Mr3zNtgTXtDxcuQ7A/gzCk+2U0adwH3hiuweFeaqU9NZFlz3Nh1cr0qjLz5NIZcchK/7DrABS99x90fzWfLnoN+hxbzip08RKQHcD2wH2d8q9tUdWq4AotnVoYeW6w5bmxLShL+1K4ukwZmclPXxnw5fyPdBmfx2uSfrGlvKXip85gIPKSq00TkBOA94C5VnRjOAEMh2uo8rAw99lhz3PixZtt+nhi9hAnLttCwRgUe7NOS7i1q+h1WVArL2FYiUhv4VFVPK01wkRBtySMSYy4ZYwqXtXwLj41ewuqt+8lsns6DfVrROL2i32FFlbCMbaWqm4AeJY4qgVkZujH+y2x+LGNvP4MHzm3JnLU7OfvZKTz51RL2WNPeYinV2FaqeiBUgSQSK0M3Jjqklkniui6NmDgwkz+1rcv/TVtD98FZfDRrA3l5MdcrIaJsYEQfRGLMJWNM8aVXKsu/Lz6RLwacTv1qadzz6QIuePk75qzb4XdoUcvm8zDGmACqyhfzfuGpr5fy655DXNimDvee04Kax5TzO7SIs/k8jDGmmESEC9rUYeLdmdyS2ZivFmyi2+AsXpq0ioNHcv0OL2pY8jDGmAJUKFuGe3q1YPxdZ3B6kxo8/c1yej47hXGLN5MIJTZFseRhjDGFyKhegdevas9713YgtUwSN7w3h6venMmqLXv9Ds1XljyMMaYYujRN5+vbu/BQn1bM27CLs/87lUdHLWb3gcRs2mvJwxhjiiklOYlrOjcka2Amf25fj7e/X0u3wVmMmLGe3ARr2mvJwxhjPKpesSxPXXQCo/7WmcbpFbj/s4Wc/+I0Zq1NnKa9ljyMMaaEWtepzEc3duL5y9uwY/9hLnn1B257/0c27Y7//tOWPIwxphREhPNPOo4Jd3fltu5NGLt4M90HT+aFCSvjummvJQ9jjAmBtNQy3NWzORPu6kpm83SGjF/Bmc9MZuyiTXHZtNeShzHGhFC9amm8cmU7Rlx3KhVSy3DTsLlc+cYMlm+Or6a9ljyMMSYMTmtSg69u68yj5x/Poo176P38VB75cjG7s+Ojaa8lD2OMCZMyyUlcfVoDJg3M5PIO9Xj3h7VkDp7EsOnrYr5pryUPY4wJs2oVUnnighMYfWsXmtWsxAOfL6LPC9OYsXq736GVWMSSh4jUE5FJIrJERBaLyO1HvX+3iKiI1HBfi4g8LyKrRGSBiLQN2PZqEVnpLldH6jMYY0xptDruGD64oSMvXtGG3dmHuXTodAaMmMvGXbHXtLdMBM+VA9ytqnNFpBIwR0TGq+oSEakH9ATWB2x/DtDUXU4FXgFOFZFqwMNAe0Dd43ypqjsj+FmMMaZERIQ+Jx5HjxY1eW3KT7yS9RMTlv7KTV0bc1PXxpRLSfY7xGKJ2J2Hqm5S1bnu873AUqCO+/azwD04ySBfX+BddUwHqrjzpp8NjFfVHW7CGA/0itTnMMaYUCifmswdZzZjwt1d6dGiJv/9diU9hkzmqwWx0bTXlzoPEWkAtAFmiEhfYKOqzj9qszrAhoDXP7vrgq0/+hw3iMhsEZm9devWEEZvjDGhU7dqGi/1a8v713ekUrkyDBgxl8tfn87STXv8Dq1QEU8eIlIR+BS4A6co637goVCfR1WHqmp7VW2fnp4e6sMbY0xIdWpcndG3dubxC1qzbPNezn1+Kg9+void+w/7HVqBIpo8RCQFJ3EMV9WRQGOgITBfRNYCdYG5IlIL2AjUC9i9rrsu2HpjjIlpZZKT6N8xg6yBmfTvmMGImevpNiSLd39YS05unt/h/U4kW1sJ8AawVFWfAVDVhap6rKo2UNUGOEVQbVV1M/AlcJXb6qojsFtVNwHfAD1FpKqIVMWpaP8mUp/DGGPCrUpaKo/2bc2Y27rQqvYxPPTFYvq8MI3vf9rmd2j/E8k7j9OB/kB3EZnnLr0L2X4MsBpYBbwO3AKgqjuAx4FZ7vKYu84YY+JK81qVGH7dqbx6ZVv2HcrhitdncMvwOfy8M9vv0JBYqNUvrfbt2+vs2bP9DsMYY0rs4JFchk5ZzctZq1CFG7s25uaujSmfGr6mvSIyR1XbF/Se9TA3xpgYUC4lmdt6NGXi3Zn0PL4Wz09YSY8hWYya/4svTXsteRhjTAw5rkp5Xri8DR/d2Ikqaanc+v6PXDp0Oot/2R3ROCx5GGNMDOrQsBqjbu3MPy88gZW/7uW8F6Zx/2cL2b7vUETOb8nDGGNiVHKScMWp9cka2I2rOjXgw1kb6DY4i7e+W8ORMDftteRhjDExrnJaCo+cfzxf396FE+tW4dFRS+j93FSmrQxf015LHsYYEyea1azEe9d24LX+7TiYk8uVb8xgwPC5YalQj+SousYYY8JMRDj7+Fp0bZbOG9PWcOBwLk4f7dCy5GGMMXGoXEoyA7o1CdvxrdjKGGOMZ5Y8jDHGeGbJwxhjjGeWPIwxxnhmycMYY4xnljyMMcZ4ZsnDGGOMZ5Y8jDHGeJYQk0GJyFZgXSkOUQOInvkf/yja44PojzHa4wOLMRSiPT6IrhgzVDW9oDcSInmUlojMDjabVjSI9vgg+mOM9vjAYgyFaI8PYiNGsGIrY4wxJWDJwxhjjGeWPIpnqN8BFCHa44PojzHa4wOLMRSiPT6IjRitzsMYY4x3dudhjDHGM0sexhhjPLPkEUBE6onIJBFZIiKLReR2d/3TIrJMRBaIyGciUiXaYgx4/24RURGpEW3xicit7nVcLCL/8SO+wmIUkZNFZLqIzBOR2SLSwaf4yonITBGZ78b3qLu+oYjMEJFVIvKhiKT6EV8RMQ4XkeUiskhE3hSRlGiLMeD950VkX7TFJ44nRWSFiCwVkdv8irFQqmqLuwC1gbbu80rACqAV0BMo467/N/DvaIvRfV0P+AanQ2SNaIoP6AZ8C5R13zs22q4hMA44x13fG8jyKT4BKrrPU4AZQEfgI+Ayd/2rwM0+XsNgMfZ23xPg/WiM0X3dHngP2Bdt8QF/Bd4Fktz3fPtbKWyxO48AqrpJVee6z/cCS4E6qjpOVXPczaYDdaMtRvftZ4F7AN9aQRQS383Av1T1kPveliiMUYFj3M0qA7/4FJ+qav4v4hR3UaA78Im7/h3gAh/CA4LHqKpj3PcUmIm/fysFxigiycDTOH8rvink3/lm4DFVzXO38+1vpTCWPIIQkQZAG5xfA4GuAb6OdDwFCYxRRPoCG1V1vq9BBTjqGjYDurjFLpNF5BQ/Y8t3VIx3AE+LyAZgMHCfj3Eli8g8YAswHvgJ2BXwI+ZnfvvR4IujY1TVGQHvpQD9gbF+xefGUVCMfwO+VNVNfsYGQeNrDFzqFp1+LSJN/Y2yYJY8CiAiFYFPgTtUdU/A+kFADjDcr9gCYvlfjDgx3Q885GtQAQq4hmWAaji35X8HPhIR8THEgmK8GbhTVesBdwJv+BWbquaq6sk4v9w7AC38iiWYo2MUkdYBb78MTFHVqf5E5yggxjOAS4AX/IwrX5BrWBY4qM4QJa8Db/oZYzCWPI7i/mL6FBiuqiMD1v8F6AP0c2/JfVNAjI2BhsB8EVmL8x9xrojUipL4wPmlPNK9VZ8J5OEMAOeLIDFeDeQ//xjnS9tXqroLmAR0AqqISBn3rbrARt8CCxAQYy8AEXkYSAfu8jOuQAExdgOaAKvcv5U0EVnlZ2zwh2v4M7/9P/wMONGvuApjySOA+0v4DWCpqj4TsL4XTvno+aqa7Vd8bix/iFFVF6rqsaraQFUb4Pzna6uqm6MhPtfnOH+4iEgzIBWfRg4tJMZfgK7u8+7AykjHBiAi6fkt+kSkPHAWTr3MJOBid7OrgS/8iM+Nq6AYl4nIdcDZwOX5ZfZRFuMcVa0V8LeSrapNoii+ZQT8reD8f1zhR3xFsR7mAUSkMzAVWIjzyxic4qDncW4lt7vrpqvqTZGPMHiMqjomYJu1QHtVjfiXcyHX8Fuc2++TgcPAQFWdGOn4iohxD/AcThHbQeAWVZ3jQ3wn4lSIJ+P8wPtIVR8TkUbABzjFfz8CV+Y3QIiiGHNwWvvtdTcdqaqPRVOMR22zT1UrRlN8bkIZDtQH9gE3RVNdZj5LHsYYYzyzYitjjDGeWfIwxhjjmSUPY4wxnlnyMMYY45klD2OMMZ5Z8jDGGOOZJQ9jjDGeWfIwYSMi5d1BEJNFpIqI3OJ3TEWJRJwi8n2Yj1/oHBWhPr+IPCIiA0N5zNISkVQRmRIwnIsJMUseJpyuwelhnAtUAaIiebiT7QT7v+85ziKO9weqepqX44ea3+cPxut1LIyqHgYmAJeG4njmjyx5mCKJM+veWe7zJ0SkuCOS9uO38Zf+BTQWZ5a+p0XkSnFmUZsnIq+5cywgIg3EmW3wbXcmteEicqaIfCciK0WkQ8A2w8WZae0TEUkLiPcPx3b3WS4i7wKLgHoi8rmIzBFnFrcbgsTZQEQWBRx7oPtLu6DjFfiZCrie+wI+61IRed2NYZw7xtHR218lziyW80XkvcI+Z8B7FUTkK3efRSJyacB7+4J9rqL2Ddh+kPvvMw1oftR7wf5tH3Sv2TQRed89Z7GvY5B/18Ji/Rzn/6AJh5LMIGVLYi3AGUAWzh/iV0ByMfZJBTYHvG4ALHKftwRGASnu65eBqwK2ywFOwPlxMwdnTCwB+uJ8ITTAmTTndHefN3HGygp6bHefPNyZ5Nz3qrmP5XG+uKoHxnl03O7rgcAjRx+vsM9UwLXZd9RnPdl9/RHOeFWB2x6PMzBejaNiLuwa7gP+BLwecJzKgecP9rnc50H3dV+3wxkXLA1n8qxVxbj+pwDzgHI4szeudM9ZrOtYyPrCPmcysNXvv594Xaw80BRJVaeIiOAMsZ2pqrniDNI3COeP9eICdqsB7ApyyB44X0CznMNSHmcynHxrVHUhgIgsBiaoqorIQpwvG4ANqvqd+3wYcBvOBE7Bjj0FWKeq0wPOc5uIXOg+rwc0BbyMRBx4vKI+UzBrVHWe+3xOwOfL1x34WN1BLlV1RzHPtxAYIiL/Bkart3k1itq3C/CZuiNMi8iXAe8Fi6sa8IWqHgQOisiogH2Kcx2PCbJ+RLBY3f+nh0WkkjozRpoQsuRhiiQiJ+DM+709/49QVVcD14rIJ0F2O4DzK7PAQwLvqGqwmfoCR4rNC3idx2//Z48e0TP/dYHHFmfGwP0BrzOBM4FOqpotIllB4s3h98W7gdvsD3he1GcKJvCz5uJ8KRZHoedT1RUi0hZnTvEnRGSC/n5E2aCfqxj7eo5LRO4oZJ8ir6OI3FrQeve9wmItizNCsgkxq/MwhRKR2jjDQ/cF9okzt0mRVHUnkCwi+V9Ke3GKK8CpyLxYRI51z1FNRDI8hlZfRDq5z68Apnk8dmVgp5s4WuDMcHh0nAC/AseKSHURKYszIVhBQvGZCjIRuEREqucftzjnE5HjcOaqGIYzX3fbo44b9HMVY98pwAXitKarBJwX8F6wuL4DzhORcuLM4Oj1Oha4vrBY3Wu2TVWPBDmXKQW78zBBiVMJPRK4W1WXisjjwL8p/rzU44DOwLequl2cSu9FOHPAPwCME6d1zRFgAM48EMW1HBggIm8CS4BXAFR1iYgUdOyji6PGAjeJyFL3WNPd/X8Xp6r+XUQeA2bizNy3rKBgCjmvl89U0HEXi8iTwGQRycWZx+MvxTjfCTjzsee579181HGPFPK5itp3roh8CMzHKTqaVdR1UNXpbvHWApzEtRDYXcDnLWz/gj5v5UJi7YZTR2fCwObzMCXi/qp7Emf2s/9T1acK2KYtzpzg/UN87gY45duti9jURBERqaiq+9wfJVOAG1R1bhjPNxK4V1Wjcia+WGd3HqZEVHU7UOhsiu4v1EkikqxOXw+T2IaKSCuc+pV3wpw4UoHPLXGEj915GGOM8cwqzI0xxnhmycMYY4xnljyMMcZ4ZsnDGGOMZ5Y8jDHGeGbJwxhjjGeWPIwxxnj2/wxSdxcKTl6kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lBIQo42pc4x6"
   },
   "source": [
    "### How linear separators work (Support Vector Machines)\n",
    "The line represents a linear separator with parameters $\\hat{\\theta} = [57, 1], \\hat{\\theta}_0 =4452$, that is, all points in the line satisfy the following equation: \n",
    "\n",
    "$\\hat{\\theta}·x + \\hat{\\theta}_0 = 0$\n",
    "\n",
    "conform a line (1 dimensional hyperplane) that *classifies* all points of the 2 dimensional space into two categories. The objective is to configure the parameters of this line so that all data points corresponding to the failure state fall into one category, and that all data points corresponding to the normal operation state fall into the other category. This way, we can use the expression $\\hat{\\theta}·x + \\hat{\\theta}_0$ to explain the status of the engine from a reading x:\n",
    "\n",
    "$\\hat{\\theta}·x + \\hat{\\theta}_0 \\geq 0 \\rightarrow$ the engine is in failure state\n",
    "\n",
    "$\\hat{\\theta}·x + \\hat{\\theta}_0 \\leq 0 \\rightarrow$ the engine is in normal state\n",
    "\n",
    "We can use the sign operator to provide a numeric value to the labels that represent both status: \n",
    "\n",
    "$h(x) = \\text {sign(}  \\hat{\\theta} \\cdot x +  \\hat{\\theta}_0 \\text {)}$\n",
    "\n",
    "If we read new sensor data and multiply it by the coefficients of our linear separator, we can predict the status of the machine as a failure if the sign is positive (the result is greater than zero) and we can predict the status of the machine as normal state if the result is negative. This type of classifiers are also known as **Support Vector Machines**.\n",
    "\n",
    "## Learning the parameters of an SVM\n",
    "### The training dataset\n",
    "Now, the objective is to learn the parameters of our separator from the dataset, a set of historic readings from which we know the status of the machine. Let a sample t of our data set be noted as:\n",
    "\n",
    "- $x^{(t)} = [x_1^{(t)}, x_2^{(t)}]$ temperature sensor and rpm readings at sample t\n",
    "- $y^{(t)}$ status at sample t {-1 if machine is in normal state, +1 if machine is in failure state}\n",
    "\n",
    "We are going to have many data readings, some of them in normal state, and some in failure state, and we are going to use these data to find the parameters of a separator that divides the space into two regions based on the information that we have in the database we have collected. \n",
    "\n",
    "Let us note our training data set as:\n",
    "\n",
    "$S_n = [x^{(t)}, y^{(t)}] \\forall t \\in [1, ..., n]$\n",
    "\n",
    "### Reducing the learning problem to an optimization problem\n",
    "Now, given that we are going to use the sign of the result to determine the status of the machine, every separator must ensure that: \n",
    "\n",
    "$\\hat{\\theta}·x^{(t)} + \\hat{\\theta}_0 \\geq 0 \\quad \\forall y^{(t)} = 1$\n",
    "\n",
    "$\\hat{\\theta}·x^{(t)} + \\hat{\\theta}_0 \\leq 0 \\quad \\forall y^{(t)} = -1$\n",
    "\n",
    "There might be many lines that satisfy these constraints, or there might be actually none, that is, the space might not be separable by a line. If the space can be separated by a line, the objective is to find a hard margin separator.\n",
    "Note that it is beneficial to maximise the **margin** of our separator, that is, the distance between the line and the closest points above and below it. Given that there is uncertainty on the information we have collected (for instance, we might not have gathered data in every possible situation), it is better to find a margin that maximises the distance between the distance between the data clusters in the two regions. In fact, this is what a human will solve the problem. \n",
    "This is known as the maximum margin SVM, i.e., to find the linear hyperplane with the maximum margin, such that all points are classified correctly. \n",
    "Furthermore, the parameters  $\\hat\\theta$  and  $\\hat\\theta_0$  of the maximum margin separator satisfy the following constraint: \n",
    "\n",
    "$y^{(t)}*(\\hat\\theta·x^{(t)}+\\hat\\theta_0) \\geq 1 \\quad \\forall t \\in [1,...,n]$ \n",
    "\n",
    "For convenience, we have set up the minimum distance to the separator to 1. This yields that the margin of our hyperplane is given by $\\frac{2}{\\left\\|{\\hat{\\theta}}\\right\\|}$. \n",
    "For instance, in our example, imagine that point $x_A$ is the closest point above the line ($y^A=1$) and point $x_B$ is the closest point below the line($y^B=1$). Let us now define a line that is parallel to the separator and passes through $x^A$:\n",
    "\n",
    "$\\hat\\theta·x^A+\\hat\\theta_0 = 1$\n",
    "\n",
    "$\\hat\\theta_0·x^A_0 + \\hat\\theta_1·x^A_1 +\\hat\\theta_0 - 1 = 0$\n",
    "\n",
    " and another one that is also parallel to the separator and passes through $x^B=-1$: \n",
    " \n",
    "$-\\hat\\theta·x^B-\\hat\\theta_0 = 1$\n",
    "\n",
    "$\\hat\\theta_1·x^B_1 + \\hat\\theta_1·x^B_2 +\\hat\\theta_0 + 1 = 0$\n",
    " \n",
    "The [distance between both parallel lines](https://en.wikipedia.org/wiki/Distance_between_two_parallel_lines) is noted as $\\omega$ and is given by:\n",
    "\n",
    "$\\omega = \\frac{|(\\hat\\theta_0+1)-(\\hat\\theta_0-1)|}{\\sqrt{\\hat\\theta_1^2 + \\hat\\theta_2^2}} = \\frac{2}{\\left\\|{\\hat{\\theta}}\\right\\|}$  \n",
    "\n",
    "Therefore, our objective function is to maximise this expression, or equivalently, to minimise $\\left\\|{\\hat{\\theta}}\\right\\|/2$, which is a function of the components of the vector $\\hat{\\theta}$ and $\\hat{\\theta}_0$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "klRkt4Vr72Wb"
   },
   "source": [
    "With this, our non-linear problem becomes:\n",
    "\n",
    "$\\min \\frac{\\sqrt{\\hat\\theta_1^2 + \\hat\\theta_2^2}}{2}$\n",
    "\n",
    "subject to:\n",
    "\n",
    "$y^{(t)}*(\\hat\\theta·x^{(t)}+\\hat\\theta_0) \\geq 1 \\quad \\forall t \\in [1,...,n]$ \n",
    "\n",
    "Or equivalently\n",
    "\n",
    "$\\min \\frac{\\sqrt{\\hat\\theta_1^2 + \\hat\\theta_2^2}}{2}$\n",
    "\n",
    "subject to:\n",
    "\n",
    "$\\hat\\theta·x^{(t)}+\\hat\\theta_0 \\geq 1 \\quad \\forall y^{(t)} =1$\n",
    "\n",
    "$\\hat\\theta·x^{(t)}+\\hat\\theta_0 \\leq -1 \\quad \\forall y^{(t)} =-1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sQhtcjJ-9tLb"
   },
   "source": [
    "## SVMs in Python\n",
    "Now, we are going to use two different approaches in Python to find the maximum margin separator for our dataset, using the package for data science [Scipy](https://www.scipy.org/).\n",
    "\n",
    "### Using non-linear programming\n",
    "In the first approach we are going to use the package **optimize** for optimization. More specifically, we are going to use the [LinearConstraint](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.LinearConstraint.html), the [Bounds](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.Bounds.html) modules, and the [minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) function. The script below imports the dependencies and packs the features into a matrix."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mtwRULBIWouV"
   },
   "source": [
    "from scipy.optimize import LinearConstraint, Bounds, minimize\n",
    "\n",
    "# Pack the 2 sample arrays into 1 array using np.append\n",
    "failure_features = np.append(temp_readings_failure, rpms_failure,1)\n",
    "normal_features =np.append(temp_readings_normal, rpms_normal, 1)\n"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "S1ZO_qImBFO_"
   },
   "source": [
    "#### Minimize\n",
    "The ```scipy.optimize.minimize``` function takes the following parameters:\n",
    "\n",
    "- ```fun```: A callable function that provides the objective function as a function of a one dimensional array, where each element of the array is a decision variable. \n",
    "\n",
    "- ```x0```:  An initial guess, or initial values for our decision variables\n",
    "\n",
    "- ```method```: (Optional) The algorithm used for optimization\n",
    "\n",
    "- ```constraints```: List or dictionary containing the constraints for the optimization problem. \n",
    "\n",
    "Given that our problem has constraints, the function will select the [L-BFGS-M](https://en.wikipedia.org/wiki/Limited-memory_BFGS) algorithm. \n",
    "\n",
    "Let us first define the objective function to minimize:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "qBW_3ixVBiJu"
   },
   "source": [
    "def objective_func(theta):\n",
    "  return 1/2*(theta[0]**2 + theta[1]**2)**(1/2)"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yWrJjarZByBj"
   },
   "source": [
    "#### Constraints\n",
    "Now, for the constraints, we can either introduce linear constraints or non-linear constraints. The linear constraints are instances of ```scipy.optimize.LinearConstraint```. This function has the following parameters:\n",
    "\n",
    "- ```A```: Matrix defining the coefficients in the left-hand-sides of the constraints\n",
    "- ```ub```: A vector defining the upper bound of the constraints or a scalar defining the upper bound for all the constraints\n",
    "- ```lb```: A vector defining the lower bound of the constraints or a scalar defining the lower bound for all the constraints\n",
    "\n",
    "Note that, playing with the values provided to ```ub``` and ```lb``` it is possible to defines constraints of type less or equal, greater or equal, or equal. \n",
    "With this function, the following code snippet defines the set of constraints corresponding to the failure state data set samples"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "Z6ZzOLaIBe-l"
   },
   "source": [
    "# The right hand sides are x_1^(t), x_2^(t), 1, so we append 1s to the failure features to obtain all right hand sides:\n",
    "failure_rhs_coefs = np.append(failure_features,np.array(np.ones((no_samples,1))),1)\n",
    "# For the label +1, we do not need to make any change to the RHS, the lower bound is 1, i.e. 1 <= x_1(i)*theta_1 + x_2(i)*theta_2 + theta_0 <= inf\n",
    "failure_linear_constraints = LinearConstraint(failure_rhs_coefs, lb=1, ub=np.inf)"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lF62FFkrJtD5"
   },
   "source": [
    "And the following code snippet to define the set of constraints for the normal state condition"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "KcnA2L43EnhT"
   },
   "source": [
    "normal_rhs_coefs = np.append(normal_features, np.array(np.ones((no_samples,1))),1)\n",
    "normal_rhs_coefs = normal_rhs_coefs\n",
    "# For the label +1, we do not need to make any change to the RHS, the lower bound is 1, i.e. -np.inf <= x_1(i)*theta_1 + x_2(i)*theta_2 + theta_0 <= -1\n",
    "normal_linear_constraints = LinearConstraint(normal_rhs_coefs, lb=-np.inf, ub=-1)\n",
    "\n"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PsVTL9CfLuru"
   },
   "source": [
    "#### Example\n",
    "Now, we have everything that we need to build our separator! let us use the initial separator at the beginning of the function as the initial guess:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 577
    },
    "id": "aVROIj-wLvFh",
    "outputId": "606ac655-4996-40c8-d663-ed126461d59d"
   },
   "source": [
    "theta_hat_0 = np.array([57, 1, -4452])\n",
    "linear_constraint = [failure_linear_constraints, normal_linear_constraints]\n",
    "print(len(linear_constraint))\n",
    "# x_0 = [ 1.13792022e+00,  5.42842896e-03, -4.83177976e+01]\n",
    "res = minimize(objective_func, theta_hat_0, constraints=linear_constraint, \n",
    "               options={\"maxiter\":1000, \"disp\": True})\n",
    "print(res)\n",
    "\n",
    "theta = res.x\n",
    "\n",
    "\n",
    "#Prepare the figure \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "#Plot the data in failure state in red\n",
    "ax.scatter(temp_readings_failure, rpms_failure, color='red')\n",
    "#Plot the data in normal operation in blue\n",
    "ax.scatter(temp_readings_normal, rpms_normal, color='blue')\n",
    "\n",
    "# Make a linear space to plot the temperature\n",
    "t = np.linspace(22,37)\n",
    "\n",
    "\n",
    "plt.xlabel('$x_1$ (temperature in celsius degrees)')\n",
    "plt.ylabel('$x_2$ (engine speed in rpms)')\n",
    "\n",
    "r_mms = (-theta[0] /theta[1])*t - (theta[2] / theta[1])\n",
    "ax.plot(t, r_mms)\n",
    "\n",
    "mod = (abs(theta[0]**2 + theta[1]**2))**-1/2\n",
    "print(mod)\n",
    "r_mms_max = (-theta[0] /theta[1])*t + (-theta[2] + 1)/ theta[1]\n",
    "\n",
    "r_mms_min = (-theta[0] /theta[1])*t + (-theta[2] - 1)/ theta[1]\n",
    "\n",
    "ax.plot(t, r_mms_min, color='red')\n",
    "ax.plot(t, r_mms_max, color='black')\n",
    "ax.set_xlim(21, 37)\n",
    "ax.set_ylim(2300, 3300)\n",
    "print(np.min(np.sum(failure_rhs_coefs*theta,axis=1)))"
   ],
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "2\n",
      "Optimization terminated successfully.    (Exit mode 0)\n",
      "            Current function value: 0.7663190555426902\n",
      "            Iterations: 10\n",
      "            Function evaluations: 51\n",
      "            Gradient evaluations: 10\n",
      "     fun: 0.7663190555426902\n",
      "     jac: array([0.00503157, 0.49997468, 0.        ])\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 51\n",
      "     nit: 10\n",
      "    njev: 10\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([ 1.54231359e-02,  1.53256051e+00, -4.45199915e+03])\n",
      "0.2128583851271692\n",
      "1.0\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de7RcdX338feHEC4RJAIBLSQ5qFiKQjFGwIq2gnKrPtiLVRsw1dYskVYRtI9KBUXTeml1SStSLKxySaUoQamggIhi9OGSIBAgIKkSAaOgEA2iXJLv88f+jZmczJwzvzl7Zvae+bzWmnVm9vU7+5wz39m/qyICMzOzHFsNOgAzM6sfJw8zM8vm5GFmZtmcPMzMLJuTh5mZZdt60AH0w6677hpjY2ODDsPMGlauhCee2HL5NtvAfvv1Px5racWKFT+LiFmt1o1E8hgbG2P58uWDDsPMGrZqU+jx5JPg/9XKkLSm3ToXW5lZ/82Zk7fcKsfJw8y2tGQJjI0VdwhjY8XrMi1eDDNmbL5sxoxiudWCk4eZbW7JEli0CNasgYji56JF5SaQBQvg7LNh7lyQip9nn10st1rQKAxPMn/+/HCdh1mHxsaKhDHe3Llw7739jsYGSNKKiJjfap3vPMxscz/6Ud5yG0lOHma2OVdmWwecPMxsc1WszO51Bb5lc/Iws81VrTK7HxX4ls0V5mZWba7AHxhXmJtZfbkCv5KcPMys2lyBX0lOHmZWbVWswDcnDzOruKpV4BswIqPqmlnNLVjgZFExvvMwM7NsTh5mZpbNycPMzLL1LXlI2k7SjZJulXSHpA+l5Usk3S3pdknnSpqelkvSGZJWS7pN0rymYy2UdE96LOzXezAzs0I/7zweBw6NiN8HDgCOlHQwsATYB9gP2B74m7T9UcDe6bEI+CyApJ2B04CDgAOB0yQ9o4/vw8xs5PUteUTh0fRyenpERFyR1gVwI7Bn2uYY4Py06npgpqRnAUcAV0fEwxHxCHA1cGS/3oeZmfW5zkPSNEm3AA9SJIAbmtZNB44DvpYW7QHc17T7/WlZu+Xjz7VI0nJJyx966KFy34iZ2Yjra/KIiA0RcQDF3cWBkl7QtPpM4LqI+HZJ5zo7IuZHxPxZs2aVcUgzM0sG0toqItYB15KKmySdBswCTmra7AFgdtPrPdOydsvNzKxP+tnaapakmen59sCrgLsk/Q1FPcYbI2Jj0y6XAW9Kra4OBn4REWuBK4HDJT0jVZQfnpaZmVmf9HN4kmcB50maRpG0Lo6Ir0h6ClgD/D9JAEsj4nTgCuBoYDXwGPBmgIh4WNKHgZvScU+PiIf7+D7MzEZe35JHRNwGvLDF8pYxpNZXJ7RZdy5wbqkBmplZx9zD3KyuPK+3DZBH1TWro8a83o89VrxuzOsNHn3W+sJ3HmZ1dMopmxJHw2OPFcvN+sDJw6yOPK+3DZiTh1kdeV7v+hqSuionD7OqyPlQ8bze9dSoq1qzBiI21VXVMIE4eZhVQe6HylTn9R6Sb7+1M0R1VSq6Uwy3+fPnx/Llywcdhll7Y2NFwhhv7ly4995yzzW+pRYUdy05yce6s9VWxZeD8STYuHHL5QMmaUVEzG+1znceZlXQzwrwIfr2WztDVFfl5GFWBf38UHFLrcEZoroqJw+zKujnh8oQffutnanWVVWIk4dZFfTzQ2WIvv3W0oIFRT3Wxo3FzxomDvDwJGbVsWBBfz5IGuc45ZSiqGrOnCJx1PRDzAbDycNsFPUrUdnQcrGVWdXVrU9G3eK1rvjOw6zK6jZ6bt3ita75zsOsyurWJ6Nu8Q5aje/SfOdhVmV165NRt3gHqeZ3adl3HpKeluYhN7Neq1ufjLrFO0g1v0ubNHlI2krSX0q6XNKDwF3AWkl3SvqEpOf2PkyzEVW3Phl1i3eQan6X1smdx7XAc4D3Ac+MiNkRsRtwCHA98DFJx/YwRrPR0Kr8u249kusW7yDV/C5t0lF1JU2PiCenus0geVRdqzyPdDt6avA7n9Kouo2kIOl1knZMzz8gaamkec3bmFmXal7+PXLKaCVV87u0jufzkHRbROwv6RDgI8AngFMj4qBeBlgG33lY5dVsnoeRVoM7hrKUNZ/HhvTzj4GzI+JyYJupBmdm1L78e6T4LhHISx4PSPp34PXAFZK2zdzfzNpxK6X6qHkrqbLkfPj/BXAlcERErAN2Bt7Tk6jMRsmSJZu+zU5LXahqVv49UnyXCGQkj4h4LCKWRsQ96fXaiLiqd6GZjYBG+Xlj/vINGzbdcThxVJPvEoGM5CFpvqRLJd0s6TZJKyXd1svgzIaey8+ro9MWVDVvJVWWnNZWd1MUU60Eftv8IyLW9Ca08ri1lVWWW1lVwwi1oMpRVmurhyLisoj4YUSsaTxKitFsNLn8vBp8B5gtJ3mcJuk/JL1R0p82Hj2LzGwUuPy8GqrWgqoGQ7XnJI83AwcARwKvSY9X9yIos5Hh8vPeyP3wrdIdYHMjiohNQ7VXLYFEREcP4O5Ot63a40UvelGY2RC58MKIuXMjpOLnhRduvm7GjIjio7d4zJix+Tatjjd+H6n4Of74vTZ37uZxNB5z5/YvhgRYHm0+V3PuPL4rad9uk5Sk7STdKOlWSXdI+lBavpekGyStlvTfkrZJy7dNr1en9WNNx3pfWn63pCO6jcnMamiyb+bd1F803wFCcRfYaMjQ72/+7YrK1qypVFFWTmurVcBzgR8AjwMCIiL273B/AU+LiEclTQeWAe8ETgKWRsRFks4Cbo2Iz0p6O7B/RLxN0huAP4mI16cE9nngQOB3gK8Dz4uIDS1PjFtbmQ2VsbFN/WKazZ0L99479RZskx2/19qdvzmhQV9ag5XV2upIiuRxOJvqO17T6c7pLujR9HJ6egRwKPDFtPw84LXp+THpNWn9YSkBHQNcFBGPR8QPgdUUicTMRsFkldtTrb8YdOV5q0YU4xMHDLw1WE7y+CnwZ8CngE8Cf5qWdUzSNEm3AA8CVwP/C6yLiKfSJvcDe6TnewD3AaT1vwB2aV7eYh8zG3aTJYeptmAbdOV5q0YU7UqIBjieVk7yOB94PvCvwL8B+wIX5JwsIjZExAHAnhR3C/vk7J9D0iJJyyUtf+ihh3p1GjPrt8mSw1RbsFWh+fSCBUUR2caNxc9GXcx4A+wPlJM8XhARfx0R16bHWymSSbYoBla8FngJMFPS1mnVnsAD6fkDwGyAtH4n4OfNy1vs03yOsyNifkTMnzVrVjdhmlkVdZIcxn/45tQLjD/+LrvA9tvDccd1VlHdiz4aVUho47VrhjX+AVwIHNz0+iDg/Iz9ZwEz0/PtgW9T1Jt8AXhDWn4W8Pb0/ATgrPT8DcDF6fnzgVuBbYG9KCrwp010bjfVNbOu5Db77aaZcE4szc2Tjz++fXPlkjBBU93c1la/CzQK2eYAdwNP0UGrK0n7U1SAT6O447k4Ik6X9GzgIooh3r8HHBsRj0vajqJY7IXAwynB/CAd6xTgLencJ0bEVyc6t1tbmVlXclte9aulVp/G4pqotVVO8mhT6FaICo9z5eRhZl3Jbfbbr4Eu+5SkJkoeW7da2OIAAjZGxH2TbmxmNizmzGn9IT1Ri6yc7bs16ObEdFhhnsq+ruhxLGZm1ZJbUd2viu1BNycmr7XVzZJe3LNIzMyqJrfZb78GuqxA66ucOo+7KHqYrwF+RebwJIPkOg8z64nG/PM/+lHxrb+f0wf34dxlDU9yBPAciuFEsocnMbMaq8H8ElnKeD+TDdDY62s2lb4sJeiowhyq3ZrKzHpofLPQxock1HPekbLez2Sj9w7TNWuh42KrOnOxldkUDHqU2bJN9H4WL+68KGiiZrntWl3V7JqV0s+jzpw8zKagX30X+qXd+4Gi0rnTjnftktBEanbNSqnzSJMz/aWk90s6tfEoL0wzq6QKNAstVbu4p03Lm0SqVYunbs9dQzkV5l+mmEvjKYrWVo2HmQ2zCjQLLdXRRxd3AM1mzIANbeaTa9fxbvzsg5Op8zVroeMKc2DPiDiyZ5GYWTU1imwG1SS1TEuWwHnnbV5sJcHChXDFFfm9wxcsKB4TFYU16kDqes3ayEke35W0X0Ss7Fk0ZlZNjQ/JumvVQiqiSByLF7cebLCTu4UhqSDPkVNsdQiwQtLdkm6TtFLSbb0KzGykDFs/iqqaaEyoqfQOH7aivQ7k3Hkc1bMozEbZsPWjqLLJBi7s9g5rmIr2OuSmul245LyL+M+z/r2rfTX5JuN2yN6jb9Sn2Nqfprzzr1sHP/0pPPkkTJ8Ou+8OM2d2d6zs63LXXcWJx5s+HfaZYKbmLq5/Vf+auvlb6uaTS+vWwX33bd5cdqutYPbstr/wvvydd3WO3sd1yTe/2v2Q7JKWRcQhktaz+e+rMbbV00uKszbuvOwrXH/9N3t+niqn9UHH1uvz//Rn3e1XelzL7i/tUIP+nbVTibgeWdtycSVia6EKcfnOowtP/ORBnvjfH+Tv2I9LXeXfZxexRRcX7Zpr4FOfgscf37Rs223hXe+Cww7bcvvjjg1++tMtl+++O5x/fvbp8x13HDz44JbLd9utbQD9+i0r83fWt8+Trv6W+iT3mvUojC1PlHumYOcjXuke5u5hPjpyR9MYeAfqPk0paparrFF1rUNuOFOu3OuZO8nawDtQ92sOCLMSOXmUbLJRmgcdW92SWjfXMzcZVKKV5YCH1zbLNWmxlaSTJlofEZ8sNaIe6GexVVUHIK1ryUg317Ob9zrIOX3MqmqqxVY7psd84Hhgj/R4GzCvrCCHRSdFJoO4A5hs6oGqyi2Cgu5KgfzF3yxPzjS01wF/HBHr0+sdgcsj4uU9jK8UVbrzGNQdwMArhbtU1Ts5s1FQVoX57sATTa+fSMusyWTl54O6Axh4pXCXKlEfYWZbyEke5wM3SvqgpA8CNwDn9SSqGpusyKSbYpgy1PVD2A2RzKopq5+HpHnAy9LL6yLiez2JqmRV6ucxyGIYVwqbWY6yZhIUsC+wU0R8Gvi5pANLinFk9PsOoLly/pRTivO4UtjMpiqn2OpM4CXAG9Pr9cBnSo9oyPWzGKZdH4m3v71+/T3MrFpyksdBEXEC8BuAiHgE2KYnUQ25fjULbVc5f9ZZ1ezEOJE6dnA0G2Y5yeNJSdNI43hJmgVUuJGntauEH1/NVfX+HlXutW82qnKSxxnApcBukhYDy4B/7ElUVoqcZri9bu01FXXt4Gg2zDpOHhGxBPh74J+AtcBrI+ILvQrMWsspvmlVOd9uzpkq9/cYVPNmM2svZxpaIuIu4K4exWKTyJ2ttNXMmEcfDeedt2UP9yr395hs5lAz67+sprqSjpV0ano9x011+6ub4pvxlfNnnlm/Tnd17eBoNsxyxrb6LEUF+aER8XuSngFcFREv7mWAZahSJ8GpqOv4VGVwB0ez/itrbKspNdWVNFvStZLulHSHpHem5QdIul7SLZKWN+5m0p3OGZJWS7ot9W5vHGuhpHvSY2HGe6i1uo5PVQaPemtWLf1sqvsUcHJE7AscDJwgaV/g48CHIuIA4NT0GuAoYO/0WAR8Np13Z+A04CDgQOC0dBc09KpcfNNtPwz33zCrp26a6u7eTVPdiFgbETen5+uBVRTzggTw9LTZTsCP0/NjgPOjcD0wU9KzgCOAqyPi4XT3czVwZMb7qK2qDhLYbT8M998wq6/cgRH3AQ5LL78REau6Oqk0BlwHvIAigVwJiCKZ/UFErJH0FeCjEbEs7XMN8H+BPwK2i4iPpOUfAH4dEf887hyLKO5YmDNnzovWtGquY6XodrBHz9VhVm1lDYy4HXA08ErgUODItCw3mB2AS4ATI+KXFLMTvisiZgPvAs7JPWYrEXF2RMyPiPmzZs0q45DWRrf9MNx/w6y+cufzeD5F8dW/UYywe0HOySRNp0gcSyJiaVq8EGg8/wJFPQbAA8Dspt33TMvaLbcB6bYif5QbAJjVXU7yeEFE/HVEXJseb6VIJh1JQ7qfA6yKiE82rfox8Ifp+aHAPen5ZcCbUqurg4FfRMRaiiKuwyU9I1WUH56W2YB0W5Ff5QYAZjaxnB7mN0s6OFVeI+kgIKfzxEuB44CVkm5Jy94PvBX4tKStKZoBpz7TXEFRTLYaeAx4M0BEPCzpw8BNabvTI+LhjDisZK16snfSD6Pb/cxs8HI6Ca4CfhdolEjPAe6maIIbEbF/TyIswbB0EjQz66eJKsxz7jxGojmsmZlNLqfO40Dg4YhYQ1H89Clgl4hYk5aNPHd4M7NRkZM8PhAR6yUdQtFc9xxSr28bbIc3J632fG3MeiMneWxIP/8YODsiLsfT0P7WoCYsci/t9nxtzHonp8L8KxT9KV4FzAN+DdwYEb/fu/DK0Y8K80GNeOte2u352phNTVmj6v4FRX+KIyJiHbAz8J4S4hsK3XZ4m2qxintpt+drY9Y7OdPQPhYRSyPinvR6bURc1bvQ6qWbDm9lFKu4l3Z7vjZmvZNz52ET6GbE2zLqSdxLuz1fG7PeyRpVt66q2kmwrHoSz7LXnq+NWfcmqvNw8hggV+iaWZWVNSS7JB0r6dT0ek5jyljrjotVzKyucuo8zgReArwxvV4PfKb0iEZIVWcGNDObTM7YVgdFxDxJ3wOIiEckuZPgFC1Y4GRhZvWTc+fxpKRpFHOOI2kW0MPub2ZmVlU5yeMM4FJgd0mLgWXAP/YkKrOSeGwrs97ouNgqIpZIWgEclha9NiJW9SYss6lrdMJs9KVpdMIEFxWaTVVOa6ttKca02gnYBXhdo+WVWRUNarBKs1GQU2z1ZeAYipkDf9X0MOtIv4uQPLaVWe/kJI89I+L1EfHxiPiXxqNnkdWcy9o3N4jh0Xs5tpV/vzbqcpLHdyXt17NIhojnkdjSIIqQetUJs9Xv97jjir46TiQ2KnLm87gTeC7wQ+BxQEBExP69C68c/R6exMOObGlQ8530Ymyrdr/fhhkz3NnThkMpY1tJmttqeR3mL+938hjUB2WVDVNCbff7bVbH92U2XiljW0XEmlaP8sIcHp5HYkvDNI5XJ79HV8rbsJs0eUhaln6ul/TL8T97H2L9DNMHZVmGaRyvVr/f8Ub5i4KNhkk7CUbEIennjr0PZzg0PhA9j8TmhmUcr+bf75o1RTJsLsYa9S8KNhpy6jxOarH4F8CKiLil1KhKVtX5PGw4eMIpG1YT1XnkjKo7Pz3+J71+NXAb8DZJX4iIj08tTLN6GpY7KrMcWZ0EgXkRcXJEnAy8CNgNeDnwVz2IzSrMneTMRltO8tiNon9Hw5PA7hHx63HLrQRV/nB2J0gzyym2WgLcIOnLFB0EXwP8l6SnAXf2IrhRVfXRYCfqLV6F+Mys9zquMAeQNB94KcWEUN+NiFrUQtetwrzqHercCdJsNJTSSTANyf484GnATOBoD8neG1UfDdadILtT5aJIs1wekr2Cqv7h7E6Q+VxPZMMmp85jz4g4smeR2G8tXrx5nQdU68PZnSDzuZ7Ihk3fhmSXNFvStZLulHSHpHc2rfs7SXel5R9vWv4+Sasl3S3piKblR6ZlqyW9t9uYqqoOQ3ksWFDUv2zcWPysUmxVVPWiSLNcOXcehwBvlvQDuhuS/Sng5Ii4WdKOwApJVwO7UxSH/X5EPC5pNwBJ+wJvAJ4P/A7wdUnPS8f6DPAq4H7gJkmXRcRQtfhyx7PhMmdO60YQVSmKNMuVc+dxFMV8HodTNNN9dfrZkYhYGxE3p+frgVXAHsDxwEcj4vG07sG0yzHARRHxeET8EFgNHJgeqyPiBxHxBHBR2tasp6ZS4e16Ihs2OcnjR8DLgIVpKPaguGvIJmkMeCFwA0ULrpdJukHStyS9OG22B3Bf0273p2Xtlpv1zFRnD6xDUaRZjpxiqzOBjcChwOnAeuAS4MUT7TSepB3SfidGxC8lbQ3sDBycjnWxpGfnHLPNeRYBiwDmuGzApqhVhXejr0unnThdFGnDJOfO46CIOAH4DUBEPAJsk3MySdMpEseSiFiaFt8PLI3CjRQJalfgAWB20+57pmXtlm8mIs6OiPkRMX/WrFk5YZptYbKK7V7Px25WNTnJ40lJ0yiKq5A0i+KDviOSBJwDrIqITzat+hLwirTN8ygS0s+Ay4A3SNpW0l7A3sCNwE3A3pL2krQNRaX6ZRnvwyybZw8021xO8jgDuBTYTdJiYBnwjxn7vxQ4DjhU0i3pcTRwLvBsSbdTVH4vTHchdwAXU4yb9TXghIjYEBFPAX8LXElR6X5x2tasZzx7oNnmcse22gc4jKKZ7jURsapXgZWpbmNbWTU1Jn1qN3ugK8Bt2ExpbKtU3ARARNwVEZ+JiH9rThzN25jVTadNcBsdIyPgggvccspGWyfFVtemHuCb3ZRL2kbSoZLOAxb2Jjyz3up2zKlueth7YEQbJpMWW0naDngLsADYC1gHbE+ReK4CzoyI7/U4zilxsZW106/h78fP0QIu6rLqm6jYKrfOYzpFM9pfR8S6kuLrOScPa6dfc5NUfY4Ws1ZKmc8DICKeTMOM1CZxmE2kX8Pfe2BEGzZZycNs2PRrzKmqz9FilsvJw0Zav8ac8sCINmxyxrYyG0r9GHPKE2jZsMm+85D0Kkmfk3RAer2o/LDMho8n0LJh0k2x1VuA9wDHSjoUOKDckMzqy305bFR0kzzWR8S6iHg3xcRQWUOymw2rbjscmtVRN8nj8saTiHgvcH554ZjVV6s5PzxUuw2rjpOHpE9LUkR8uXl5RPxr+WGZ1Y/7ctgoybnzWA9cJmkGgKQjJH2nN2GNNpeb15P7ctgo6Th5RMQ/AJ8HvpWSxknAe3sV2KhyuXl9uS+HjZKcYqvDgLcCv6IY3+odEfHtXgU2qlxuXl/96nBoVgUdD4wo6RvAqRGxTNJ+wAXASRHxjV4GWIY6DYzYr4H6zMwmU8rAiBFxaEQsS89XAkcBHyknRGtwubmZ1UHXY1tFxFqKKWmtRC43N7M6mNLAiBHx67ICsYLLzc2sDjwwYgX1Y6A+M7Op8JDsZmaWzcnDzMyyOXmYmVk2Jw8zM8vm5GFmZtmcPMzMLJuTh5mZZXPyMDOzbE4eZmaWzcnDzMyyOXmYmVk2Jw8zM8vm5GFmZtmcPMzMLFvfkoek2ZKulXSnpDskvXPc+pMlhaRd02tJOkPSakm3SZrXtO1CSfekx8J+vQczMyv0cz6Pp4CTI+JmSTsCKyRdHRF3SpoNHA78qGn7o4C90+Mg4LPAQZJ2Bk4D5gORjnNZRDzSx/diZjbS+nbnERFrI+Lm9Hw9sArYI63+FPD3FMmg4Rjg/ChcD8yU9CzgCODqiHg4JYyrgSP79T7MzGxAdR6SxoAXAjdIOgZ4ICJuHbfZHsB9Ta/vT8vaLR9/jkWSlkta/tBDD5UYvZmZ9T15SNoBuAQ4kaIo6/3AqWWfJyLOjoj5ETF/1qxZZR/ezGyk9TV5SJpOkTiWRMRS4DnAXsCtku4F9gRulvRM4AFgdtPue6Zl7ZabmVmf9LO1lYBzgFUR8UmAiFgZEbtFxFhEjFEUQc2LiJ8AlwFvSq2uDgZ+ERFrgSuBwyU9Q9IzKCrar+zX+zAzs/62tnopcBywUtItadn7I+KKNttfARwNrAYeA94MEBEPS/owcFPa7vSIeLh3YZuZ2Xh9Sx4RsQzQJNuMNT0P4IQ2250LnFtmfGZm1jn3MDczs2xOHmZmls3Jw8zMsjl5mJlZNicPMzPL5uRhZmbZnDzMzCybk4eZmWVz8jAzs2xOHmZmls3Jw8zMsjl5mJlZNicPMzPL5uRhZmbZnDzMzCybk4eZmWVTMefScJP0ELCm5MPuCvys5GP2guMsl+MsVx3irEOM0Js450bErFYrRiJ59IKk5RExf9BxTMZxlstxlqsOcdYhRuh/nC62MjOzbE4eZmaWzcmje2cPOoAOOc5yOc5y1SHOOsQIfY7TdR5mZpbNdx5mZpbNycPMzLI5eXRA0mxJ10q6U9Idkt6Zln9C0l2SbpN0qaSZVYyzaf3JkkLSrlWMUdLfpet5h6SPDyrGieKUdICk6yXdImm5pAMHHOd2km6UdGuK80Np+V6SbpC0WtJ/S9qmonEukXS3pNslnStpehXjbFp/hqRHBxVfUxztrqckLZb0fUmrJL2jZ0FEhB+TPIBnAfPS8x2B7wP7AocDW6flHwM+VsU40+vZwJUUnSV3rVqMwCuArwPbpnW7VfFaAlcBR6XlRwPfHHCcAnZIz6cDNwAHAxcDb0jLzwKOr2icR6d1Aj5f1TjT6/nABcCjg4xxkuv5ZuB8YKu0rmf/R77z6EBErI2Im9Pz9cAqYI+IuCoinkqbXQ/sOagYoX2cafWngL8HBtpCYoIYjwc+GhGPp3UPDi7KCeMM4Olps52AHw8mwkIUGt+Ep6dHAIcCX0zLzwNeO4DwfqtdnBFxRVoXwI0M/n+oZZySpgGfoPgfGrgJfu/HA6dHxMa0Xc/+j5w8MkkaA15IkembvQX4ar/jaac5TknHAA9ExK0DDWqccdfyecDLUlHLtyS9eJCxNRsX54nAJyTdB/wz8L7BRVaQNE3SLcCDwNXA/wLrmr7Y3M+mLxEDMz7OiLihad104Djga4OKrymWVnH+LXBZRKwdbHSbtInzOcDrU5HqVyXt3avzO3lkkLQDcAlwYkT8smn5KcBTwJJBxdasOU6KuN4PnDrQoMZpcS23BnamuPV+D3CxJA0wRKBlnMcD74qI2cC7gHMGGR9ARGyIiAMovrUfCOwz4JBaGh+npBc0rT4TuC4ivj2Y6DZpEefLgdcB/zrYyDbX5npuC/wmimFKPgec26vzO3l0KH0zugRYEhFLm5b/FfBqYEG69R6oFnE+B9gLuFXSvRR/aDdLemaFYoTi2/HSdDt+I7CRYqC3gWkT50Kg8fwLFB/WlRAR64BrgZcAMyVtnVbtCTwwsMDGaYrzSABJpwGzgJMGGdd4TXG+AngusDr9D82QtHqQsTUbdz3vZ9Pf56XA/r06r5NHB9I34HOAVRHxyablR1KUgf6fiHhsUPE1xbNFnBGxMiJ2i4ixiBij+OOaFxE/qUqMyZco/kmR9DxgGwY4kukEcf4Y+MP0/FDgnj/El0EAAAZTSURBVH7H1kzSrEYrP0nbA6+iqJ+5FvjztNlC4MuDibDQJs67JP0NcATwxkY5/SC1iXNFRDyz6X/osYh4bgXjvIum/yOKv9Pv9yyGCnxZrjxJhwDfBlZSfCOGoijoDIrbxJ+nZddHxNv6H2GhXZwRcUXTNvcC8yNiIB/ME1zLr1PcYh8APAG8OyK+MYgYYcI4fwl8mqKY7TfA2yNixUCCBCTtT1EhPo3iy+DFEXG6pGcDF1EUBX4POLbRGKFicT5F0QJwfdp0aUScPqAw28Y5bptHI2KHQcTXFEO76zmTovh8DvAo8LZe1XU6eZiZWTYXW5mZWTYnDzMzy+bkYWZm2Zw8zMwsm5OHmZllc/IwM7NsTh5mZpbNycN6RtL2aZDDaZJmSnr7oGOaTD/ilPTdHh9/wvkmyj6/pA9KeneZx5wqSdtIuq5piBYrmZOH9dJbKHoMbwBmApVIHmnCnHZ/+9lxTnK8LUTEH+Qcv2yDPn87uddxIhHxBHAN8PoyjmdbcvKwSamYUe9V6flHJHU6uugCNo2p9FHgOSpm4PuEpGNVzIR2i6R/T/MlIGlMxWyC/5lmQ1si6ZWSviPpHkkHNm2zRMVsaV+UNKMp3i2Onfa5W9L5wO3AbElfkrRCxUxsi9rEOSbp9qZjvzt90251vJbvqcX1fLTpva6S9LkUw1VpnKLx279JxWyVt0q6YKL32bTuaZIuT/vcLun1Tesebfe+Jtu3aftT0u9nGfC749a1+91+IF2zZZI+n87Z8XVs83udKNYvUfwNWi90M4OUH6P1AF4OfJPiH/FyYFoH+2wD/KTp9Rhwe3r+e8D/ANPT6zOBNzVt9xSwH8WXmxUUY14JOIbiA2GMYuKbl6Z9zqUYC6vtsdM+G0mzwqV1O6ef21N8cO3SHOf4uNPrdwMfHH+8id5Ti2vz6Lj3ekB6fTHFGFTN2z6fYnC7XcfFPNE1fBT4M+BzTcfZqfn87d5Xet523/T6RRRjfs2gmBhrdQfX/8XALcB2FDMz3pPO2dF1nGD5RO9zGvDQoP9/hvXh8kCbVERcJ0kUQ2b/UURsUDHw3ikU/6x/3mK3XYF1bQ55GMUH0E3FYdmeYkKbhh9GxEoASXcA10RESFpJ8WEDcF9EfCc9vxB4B8XkTO2OfR2wJiKubzrPOyT9SXo+G9gbyBltuPl4k72ndn4YEbek5yua3l/DocAXIg1kGREPd3i+lcC/SPoY8JXImydjsn1fBlwaaSRpSZc1rWsX187AlyPiN8BvJP1P0z6dXMent1n+X+1iTX+nT0jaMYrZIK1ETh42KUn7Uczp/fPGP2FE/AD4a0lfbLPbrym+ZbY8JHBeRLSbha959NeNTa83sulvdvyIno3XLY+tYjbAXzW9/iPglcBLIuIxSd9sE+9TbF6827zNr5qeT/ae2ml+rxsoPhQ7MeH5IuL7kuZRzBH+EUnXxOajw7Z9Xx3smx2XpBMn2GfS6yjp71otT+sminVbitGPrWSu87AJSXoWxRDPxwCPqpjDZFIR8QgwTVLjQ2k9RXEFFBWZfy5pt3SOnSXNzQxtjqSXpOd/CSzLPPZOwCMpcexDMYPh+DgBfgrsJmkXSdtSTPzVShnvqZVvAK+TtEvjuJ2cT9LvUMw7cSHF3Nvzxh237fvqYN/rgNeqaE23I/CapnXt4voO8BpJ26mYnTH3OrZcPlGs6Zr9LCKebHMumwLfeVhbKiqhlwInR8QqSR8GPkbn80xfBRwCfD0ifq6i0vt2irne/wG4SkXrmieBEyjmdejU3cAJks4F7gQ+CxARd0pqdezxxVFfA94maVU61vVp/83ijIj3SDoduJFiNr67WgUzwXlz3lOr494haTHwLUkbKObm+KsOzrcfxVzrG9O648cd98kJ3tdk+94s6b+BWymKjm6a7DpExPWpeOs2isS1EvhFi/c70f6t3u9OE8T6Coo6OusBz+dhXUnf6hZTzGD2HxHxTy22mUcx3/dxJZ97jKJ8+wWTbGoVImmHiHg0fSm5DlgUETf38HxLgfdGRM9m0xtlvvOwrkTEz4EJZ01M31CvlTQtir4eNtrOlrQvRf3KeT1OHNsAX3Li6B3feZiZWTZXmJuZWTYnDzMzy+bkYWZm2Zw8zMwsm5OHmZllc/IwM7NsTh5mZpbt/wOUidLpJeq7RAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Alternative solution: Adding a regulation parameter\n",
    "Note that the separator might not be the same as a human would have selected when solving this task. This is because the algorithm is quite sensitive to the initial guess passed as argument, and it can converge to a local minimum. To make the algorithm more resilient, we can modify slightly the objective function to factor in $\\theta_0$:"
   ],
   "metadata": {
    "id": "bzZTKqdRl7g0"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def objective_func(theta):\n",
    "  return 1/2*(theta[0]**2 + theta[1]**2 + theta[2]**2)**(1/2)"
   ],
   "metadata": {
    "id": "otiHb-J1maoj"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "With this regularization, we assume by looking at the data that the separator we are looking for crosses the y-axis at a high value, and therefore, we include it in the objective function.\n",
    "Now, if we try again with this objective function, we get a support vector machine which is closer to the solution that a human would provide:"
   ],
   "metadata": {
    "id": "hxCFOKkImeZH"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "theta_hat_0 = np.array([57, 1, -4452])\n",
    "linear_constraint = [failure_linear_constraints, normal_linear_constraints]\n",
    "print(len(linear_constraint))\n",
    "# x_0 = [ 1.13792022e+00,  5.42842896e-03, -4.83177976e+01]\n",
    "res = minimize(objective_func, theta_hat_0, constraints=linear_constraint, \n",
    "               options={\"maxiter\":1000, \"disp\": True})\n",
    "print(res)\n",
    "\n",
    "theta_prime = res.x\n",
    "\n",
    "\n",
    "#Prepare the figure \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "#Plot the data in failure state in red\n",
    "ax.scatter(temp_readings_failure, rpms_failure, color='red')\n",
    "#Plot the data in normal operation in blue\n",
    "ax.scatter(temp_readings_normal, rpms_normal, color='blue')\n",
    "\n",
    "# Make a linear space to plot the temperature\n",
    "t = np.linspace(22,37)\n",
    "\n",
    "\n",
    "plt.xlabel('$x_1$ (temperature in celsius degrees)')\n",
    "plt.ylabel('$x_2$ (engine speed in rpms)')\n",
    "\n",
    "r_mms = (-theta_prime[0] /theta_prime[1])*t - (theta_prime[2] / theta_prime[1])\n",
    "ax.plot(t, r_mms)\n",
    "\n",
    "mod = (abs(theta[0]**2 + theta[1]**2))**-1/2\n",
    "print(mod)\n",
    "r_mms_max = (-theta_prime[0] /theta_prime[1])*t + (-theta_prime[2] + 1)/ theta_prime[1]\n",
    "\n",
    "r_mms_min = (-theta_prime[0] /theta_prime[1])*t + (-theta_prime[2] - 1)/ theta_prime[1]\n",
    "\n",
    "ax.plot(t, r_mms_min, color='red')\n",
    "ax.plot(t, r_mms_max, color='black')\n",
    "ax.set_xlim(21, 37)\n",
    "ax.set_ylim(2300, 3300)\n",
    "print(np.min(np.sum(failure_rhs_coefs*theta,axis=1)))"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 577
    },
    "id": "fpQtk36Smuee",
    "outputId": "3029c66a-779f-49a5-ca7c-08e21586f634"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "2\n",
      "Optimization terminated successfully.    (Exit mode 0)\n",
      "            Current function value: 26.689281097128372\n",
      "            Iterations: 8\n",
      "            Function evaluations: 40\n",
      "            Gradient evaluations: 8\n",
      "     fun: 26.689281097128372\n",
      "     jac: array([ 1.22678280e-02,  5.22136688e-05, -4.99849558e-01])\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 40\n",
      "     nit: 8\n",
      "    njev: 8\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([ 1.30967158e+00,  5.57996292e-03, -5.33624927e+01])\n",
      "0.2128583851271692\n",
      "1.0\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3hV15W336WCQIBpFgYEEjbgjo1BILnFBWJc4t7AyGBsim4yk8lMkskk/hJnMuPMJPNN5pvMJBQDbshgTDDuxg0bN3rHGAO26aY3AUJtfX+co/haVjubW6X1Ps95dO8+53fOukfSXWfvtfdaoqoYhmEYRhBS4m2AYRiGkXyY8zAMwzACY87DMAzDCIw5D8MwDCMw5jwMwzCMwKTF24BYcPrpp2vPnj3jbYZRk4MH4fPPoXdvaNcusPzQ8XK2HTxOz06tadsy+J/y/v37+fLLLzn77LNp27ZtYL1xCqxZA2Vl325v0QL69o29PUatLFu2bJ+qZtW2r1k4j549e7J06dJ4m2HUpLwccnPh7LPhlVcCy09WVHLZv71D/9wOPDYyL7D+xIkTZGdnc/HFFzNr1qzAeuMUSKlj0KO8HOx/NWEQkS117bNhKyN+pKfDmDHw2mvwxReB5RlpqdwzsAdvr9/NzkMnAutbtWrF6NGjef7559m1a1dgvXEK5OQEazcSDnMeRnwZOxZEYPJkJ/l9g3JQYObirU76oqIiKioqmDp1qpO+yVJcDD17ej2Enj2995Hk0UchM/ObbZmZXruRFJjzMOJLjx5w880wdSqcPBlc3jGTq8/OYuaSbZRXVgXW9+nThyFDhjB58mQqKysD65skxcUwbhxs2QKq3s9x4yLrQEaM8B4YcnO9h4fcXO/9iBGRu4YRVcx5GPEnFIK9e2HOHCd5YUEue46e5M1PdjtePsS2bdt49dVXnfRNjocfhuPHv9l2/LjXHklGjIAvv4SqKu+nOY6kwpyHEX+++13o1QsmTHCSX31OZ7Lbt2L6wjpje/Vyyy230K1bNyY4Xr/JsbWOIcC62o1miTkPI/6kpMD48fD++7B2bWB5aopwX34OH23ez+a9JYH1aWlpjB07ltdff50vHAL3TQ4LZhuNwJyHkRiMHg0ZGTBxopP8nrwepKcKxQvdno7Hjh1LSkoKkyZNctI3KRIxmB3tAL4RGHMeRmJw+ulw993w1FNQErz3kNU2g6EXdGH2sm2cKAse+M7OzuaWW25h6tSpnHQI3DcpEi2YHYsAvhEYcx5G4hAKwdGj8MwzTvLCglyOlFbw0uqdjpcPsW/fPmbPnu2kb1IkUjA7VgF8IxDmPIzE4dJL4aKLvMC5Q5Gy/DM70qdzG4odA+eDBw+md+/eFjhPNCyAn5CY8zASBxGv97FyJSxa5CAXRuTnsGr7YdZsPxxYn5KSQlFRER9++CFr1qwJrDeihAXwExJzHkZiMWIEtGnjPG33jgHdaZWeSvEit97HAw88QEZGhvU+EolEDOAb5jyMBKNtW7j/fnj2Wdi/P7D8tJbp3NqvGy+s3MmR0vLA+k6dOnHPPffw9NNPc/To0cB6IwokWgDfAMx5GIlIKOSlKnniCSd5YUEuJ8ormbNsu+PlQ5SUlFBss3kSh0QK4BuAOQ8jEenbFy6/3FvzURU8X9WF2e24uEd7pi/aijoE3gsKCrj44ouZMGGCk94wmgPmPIzEJBSCTZvg7bed5IX5OWzaU8KiLw4E1ooIoVCI1atX8/HHHztd3zCaOuY8jMTkrru8hYOOgeubL+5Gu1bpzvmuRowYQdu2bS1wbhh1EDPnISItRWSxiKwSkXUi8s9+e7GIbBCRtSIyTUTS/XYRkT+KyCYRWS0i/cPONUpENvrbqFh9BiOGZGTAgw/Ciy/Cjh2B5S3TU7lrQHfmrfuKvUeDrxhv06YN999/P7NmzWLfvn2B9YbR1Illz+MkcK2qXgz0A64XkQKgGDgX6Au0Asb4x98A9PG3ccAEABHpCDwC5AODgEdEpEMMP4cRK8aP92Iejz3mJB+Rn0N5pTJr6TYnfSgUoqysjMcff9xJbxhNmZg5D/WoTlqU7m+qqq/6+xRYDHT3j7kVeMrftRBoLyJdgaHAm6p6QFUPAm8C18fqcxgx5KyzYOhQz3mUB592e1ZWGy7v3YlnFm2lsip44PvCCy/kiiuuYNKkSVQ5BO4NoykT05iHiKSKyEpgD54DWBS2Lx24H3jdb8oGwh8Zt/ttdbXXvNY4EVkqIkv37t0b2Q9ixI5QCHbuhJdecpIX5uey49AJ3t2wx/HyITZv3sybb77ppDeMpkpMnYeqVqpqP7zexSARuTBs95+BBar6foSuNVlV81Q1LysrKxKnNOLBTTd5pWodA9dDzj+Dzm0znAPnd955J1lZWRY4N4waxGW2laoeAubjDzeJyCNAFvAPYYftAHqEve/ut9XVbjRFUlO99NtvvQUbNwaWp6emMGxQDu9+tpdtB443LKhBRkYGDz74IC+99BLbtrnFTgyjKRLL2VZZItLef90K+C7wqYiMwYtjDFfV8IHlF4GR/qyrAuCwqu4C5gHXiUgHP1B+nd9mNFXGjIG0NOdCUcMG9kCAZxa7ZWEdP348qspjjoF7w2iKxLLn0RWYLyKrgSV4MY+XgYnAGcDHIrJSRH7lH/8q8DmwCXgM+D6Aqh4A/sU/xxLgN36b0VTp0gVuvx0efxxOnAgs79a+FYPPO4NZS7ZxsiJ4oagzzzyT66+/nilTplDuELg3jKZILGdbrVbVS1T1IlW9UFV/47enqWovVe3nb9Xtqqo/8Pf1VdWlYeeapqq9/c3mUTYHQiE4eBBmzXKSFxbksv9YGa+v/crx8iF27drFCy+84KQ3jKaGrTA3koOrr4Zzz3UOnF/Z+3RyO2U61zi/8cYbycnJSazAudX1NuKIOQ8jORCBoiKvSNSKFYHlKSnCfYNyWPzlATZ8FTzVempqKuPGjeOdd95hw4YNgfURx+p6G3HGnIeRPIwaBa1aOfc+7s7rQYu0FOdCUQ899BBpaWlMmjTJSR9RrK63EWfMeRjJQ/v2MHw4PPMMHA5eZrZj6xbc1Lcrc5bv4NjJisD6Ll26cMcdd/DEE09wwiFwH1GsrrcRZ8x5GMlFKATHjsH06U7ywoIcSk5W8MLKnY6XD3Hw4EGeffZZJ33EsLreyUsTiVWZ8zCSi7w8b5swwRvrD0j/nA6c26Ut0xducSr0dNVVV3HeeedFJ3Ae5EvF6nonJ00oVmXOw0g+QiFYtw4++CCwVEQoLMjlk11HWLHtkJO+qKiIxYsXs3z58sD6Ogn6pXKqdb2byNNv0tGEYlXSHMps5uXl6dKlSxs+0EgOjh+H7Gy44QYv/hGQkpMV5D/6FkMv7MIf7ukXWH/o0CGys7O57777IrfqvGdPz2HUJDfXq9kdSaodVfiXWGZmMOdjuJGSUnuPWcSp5HK0EZFlqppX2z7reRjJR2amN/Nq9mzYEzxbbpuMNG7vn83Lq3dx8FhZYH379u0ZPnw4zzzzDIcOBe+91EosA+BN6Ok36WhCsSpzHkZyUlTk1fiYNs1JXliQS1lFFbOXbXfSh0Ihjh8/ztNPP+2k/xax/FKxmVrxownFqsx5GMnJuefCNdfApElQGTxf1bldTiMvtwPFi7ZQ5VAoasCAAQwcOJAJEyY4Bd6/RSy/VJrQ02/ScaqxqgTCnIeRvBQVefGAeW5JlQsLcvly/3E+3OxWozwUCrF+/XoWLFjgpP8GsfxSaUJPv0nJiBHe321VlfczCR0HmPMwkpnbboMzznBecX5D3y50bN3CuVDUvffeS/v27SM3bTdWXypN6OnXiB/mPIzkpUULr9bHK6/UPlOpATLSUrk7rztvrd/DV4dLA+szMzN54IEHmDNnDrt37w6sjytN5OnXiB/mPIzkZtw47+l58mQn+YhBuVSpMsOxUFRRURHl5eVMnTrVSd8okm1NRrLZazhhzsNIbnJyvDrnU6ZAWfBptzmdMvlOnyxmLtlKeWXwefbnnHMO1157LZMmTaLSIXDfIMm2IjnZ7DWcMedhJD+hkLfe4/nnneQj8nPYfeQkb693G3oKhUJs3bqV119/3UlfL8m2JiPZ7I03SdxLM+dhJD9Dh8KZZzoHzq89tzNd27VkumOhqFtvvZWuXbtGJ99Vsq3JSDZ740mS99ICOw8RaS0iqdEwxjCcSEmB8ePhvfdg/frA8rTUFIYPyuGDTfv4Yt+xwPr09HTGjBnDq6++ypeRTiWSbGsyks3eeJLkvbQGnYeIpIjIfSLyiojsAT4FdonIJyLyHyLSO/pmGkYDPPigN/tq4kQn+bCBPUhLEZ5xLBQ1duxYRITJjoH7Okm2NRnJZm88SfJeWmN6HvOBXsDPgS6q2kNVOwNXAAuB34lIYRRtNIyGycqCu+6CJ5/06n0EpPNpLbnugjN4btl2SsuDB7579OjBzTffzNSpUylzCNwDtY9/J9uajGSzN54key9NVevdgPRIHBPPbcCAAWo0A95/XxVUp0xxkn+4ca/m/uxlnb10m5P+9ddfV0BnzJgRXDx9umpmpmd/9ZaZ6bUbTZMk+J0DS7WO79UGex6qWg4gIneLSFv/9S9FZI6I9A8/xjDiyuWXQ9++zoHzS3t14qys1kx3HLr67ne/y1lnneUWOE/y8e9mRyRmSSV5Ly1IwPyXqnpURK4ABgNTgShMLzEMR0S8fFfLlsGSJQ5yYUR+Liu2HmLdzuA10lNSUhg/fjwLFixg3bp1wcRJPv7drIjkLKkkXukfxHlUDwTfBExW1VeAFpE3yTBOgcJCaN3aufdxV//utExPcZ62O3r0aFq0aMHEoIH7ZB//bk5YLxEI5jx2iMgk4F7gVRHJCKg3jOhz2mmeA5k5Ew4eDCxvl5nOzRd144WVOzhaGnw0Nisri7vvvpunnnqKkpKSxgttllLyYL1EINiX/z3APGCoqh4COgI/jYpVhnEqhEJw4oQ388qBwoJcjpdV8vyKHY6XD3HkyBFmzJjROEFx8ddPs6n+EqokG/9uVlgvEQjgPFT1uKrOUdWN/vtdqvpG9EwzDEcuvhguvdRb8+FQqOniHu3pm92O6Qu3OBV6uuyyy+jbt2/jCkWFj5+DV9iqusdhjiMxsV4iEMB5iEieiDwvIstFZLWIrBGR1dE0zjCcCYVgwwaYP99JXliQw2e7S1jyZfChLxEhFAqxYsUKFi9eXP/BNn6eODR2BlWSz5KKFNLYJysR2YA3TLUG+Gv6UVV1m9cYQ/Ly8nTp0qXxNsOIJaWl0L27V6r2uecCy4+XVZD/27e55pzO/HH4JYH1R48epVu3btx555088cQTdR+YklJ770jEm4FjxIbqHmC4I8/MbJZOIRwRWaaqebXtCxLz2KuqL6rqF6q6pXqLkI2GEVlatoTRo2HuXNi1K7A8s0Uad/bvzmtrd7Gv5GRgfdu2bSksLOTZZ5/lwIEDdR9o4+eJgfUAAxPEeTwiIlNEZLiI3FG9Rc0ywzhVxo+Higqv1ocDI/JzKK9UZi3d5qQPhUKUlpbW3/Ow8fPEINFmUCVBqvYgzmM00A+4HrjZ374XDaMMIyL07g3XXecNPVRUBJb3OaMt+Wd25JlFW6msCh44v+iii7jsssuYOHFi3YFzGz+PDkG/fBOpB5gsqdrryltScwM2NPbYRNsst1Uz5vnnvZxBc+c6yV9cuUNzf/ayvrN+t5P+6aefVkDfeustJ71RB9Onq+bmqop4P8PzQbnkjKpNI+L9rHn+aJOb+007qrfc3NjZ4EM9ua2COI/HgfMbe3wt+pbAYmAVsA74Z7/9TGARsAl4Fmjht2f47zf5+3uGnevnfvsGvHUn5jyM2ikvV+3eXXXoUCf5yfJKHfAvb+qDjy920p84cUI7deqkd955p5PeqIWGnIPrl2+1Qwp3HPFIWFjz2jU/Q20OM0pEynmsB8r9L+zVeLOuVgfQC9DGf53uO4QCYBYwzG+fCIT8198HJvqvhwHP+q/P9x1Qhu94NgOp9V3bnEcz55//2ftT37TJSf7719drz396WbcdOOak/+lPf6qpqam6Y8cOJ71Rg4acQ11fviKROX+0qev6cXBo9TmPIDGP64HewHV8He+4ubFi35bqfA3p/qbAtcBsv/1J4Db/9a3+e/z9g0VE/PaZqnpSVb/A64EMCvA5jObGmDHeyu1Jk5zkwwd5494zFrsFT8ePH09lZSVTHAP3Rg0aCm6favwi3sHz2iZRiHx7SnecZ4MFcR67gTuB/wL+ANzhtzUaEUkVkZXAHuBNvF7DIVWtjmZuB7L919nANgB//2GgU3h7LRrD+DbdusFtt8G0ad76j4B075DJted05tkl2yirCL72olevXgwdOpTJkydT4RC4N2rQkHM41Rls8Q6e1zaJoqbjqCaO+bSCOI+ngAuA/wH+F2/46OkgF1PVSlXtB3TH6y2cG0QfBBEZJyJLRWTp3r17o3UZI1kIhWD/fpg9u+Fja6GwIJd9JWXMW/eV4+VD7Nixg5dfftlJb4TRkHM41RlsiTB9umaq9tzc2o+L53qgusazam7AJ41pC3C+X+GtWN8HpPltlwLz/NfzgEv912n+cYIXLP952Hn+elxdm8U8DK2qUj37bNXLLnOSV1RW6eX//rbeM/EjJ315ebl2795dr7vuOie9UYP6ZltF+vydOnlbY68VDdviVHWQCAXMpwMFYe/zgacC6LOA9v7rVsD7eHGT5/hmwPz7/usf8M2A+Sz/9QV8M2D+ORYwNxrDH/7g/cmvWuUk/9P8jZr7s5d14+4jTvrf/OY3CujGjRud9EYcCPqlHc0v+ZpOKRSK+uyrSDmP9Xg5rb70tyq/rVGzroCLgBV4M7XWAr/y28/Cm8K7yXckGX57S//9Jn//WWHnehgvXrIBuKGha5vzMFRVdf9+1ZYtVYuKnOR7j5Zq71+8oo+8sNZJv3PnTk1LS9Of/OQnTnojDgSdeRWrmVox6onU5zyCJEasY9DNQxM4z5UlRjT+ygMPwF/+Ajt3Qtu2geU/nLGC+Rv2sOgXg8lskRZYf/fddzN//ny2b99Oy5YtA+uNGBM0cWWsEl327Pl1Gv9wcnO9GEmEOOXEiP4U2SoNS4hYc4uYtYYRTUIhKCmB6dOd5IUFuRwtreClVTud9EVFRezfv5/nHDL9GnEg6MyrWM3Uivd0YhrpPPzuy6tRtsUwos+gQXDJJV6N80b2usMZ2LMDZ5/RxrnG+bXXXsvZZ5/NBMca60aMCTrzKlYzteI9nZhgU3WXi8jAqFliGLFAxOt9rFkDH33kIBcKC3JZs+Mwq7YdctIXFRXx8ccfs2rVqsB6I8YEnfYbq0SXiTCduK5gSM0N+BSowAtUB05PEs/NAubGNygpUT3tNNURI5zkR06U6Xm/fE1/Mmulk37//v3asmVLHT9+vJPeSCCiPWU4ztcmQulJhgK98NKJBE5PYhgJQ+vWMHKkV2HQYQFp25bp3NovmxdX7eTw8fLA+o4dOzJs2DCmT5/OkSNHAuvjQhLUlwhEJD5PQ6nTo33Pai4kjHUa/7q8SlParOdhfIu1a1VB9Xe/c5PvOKS5P3tZp7z/uZN+0aJFCuif/vQnJ31MidMCtagRqc9T37TcJnLPiMRU3WTGpuoatXLVVbB9O2zc6D0dBuT2P3/I4RPlvP0PV+FNSGw8qkpeXh5lZWWsXr06sD6mxGhaaMyo7/M8+qiXbHDrVi/4/OijdT/R1zctNyenSdyzSNUwN4ymRSgEn38Ob7zhJC/Mz+Xzvcf4ePP+wFoRIRQKsXbtWj788EOn68eMBJgWGlHqsrt62KmxFfzqmtlUrQ1y7SSk0c5DRDJE5D4R+YWI/Kp6i6ZxhhFVbr8dsrK8absO3HRRV9pnpjN9kdsyp+HDh9OuXbvEn7abANNCI0pddqememnOw6kv7XltM55cr52EBOl5vIBXS6MCOBa2GUZykpEBDz0EL7/s9ETYMj2Vuwd05411u9lzJHiq99atWzNy5Ehmz55NQmd+ToRpoZHkxhu9oaVwMjOhsrL24+v62wifltsYkvme1UIQ59FdVe9V1d+r6n9Wb1GzzDBiwfjx3jDDY485ye/Lz6WiSpm5ZFvDB9dCUVERZWVlTJs2zUkfE2K1diEWFBfDk09+M1YhAqNGuaU9r57xVF/MKtnvWR0EyW01GfgfVV0TXZMijwXMjXq56SZYvtx7wkxPDyy/f+oiNu0p4f1/vIa01OBhxKuvvpqtW7eyadMmUhwC90YAGgqWjxv3zaGrzMzGfek3tUkFPpEKmF8BLBORDSKyWkTWiMjqyJhoGHEkFIKvvoK5c53kI/Jz2XW4lHc+3eN4+RBffPEF87p2bTrrKBKV+oL/p9LDampDe40giPO4AeiDYw1zw0hYbrjB+6JwDFwPOa8zXU5ryfRFbjNpbi8t5Qxgwp49jZvlY7jTUPDfdeFdUxraaySNdh5q2XSNpkpqKisHjYP58zlPPg384J+WmsKwQT1Y8NletuwPPoekxSOP8BDwCvBX91PfLB/DnWj2EOK94jvGNOg8ROQD/+dRETkSth0VkSTJrWA0N4Jkhiguhttffogy0hnPRKcH/+GDckhNEZ5x6X1s3cpYQIHJNdqNCNMMewjRwlaYG02O6pRDjY17Vsc6ZzCMocwjmx2cIDNwrLPo6WUs+mI/H/98MC3TUxsv9A34HrAUr/fRApI+2GokP7bC3GhWPPxwsLVe1Q/4EwjRgUMMY+Y32htLYUEuB4+X89raXcGE/lBKCNgNzIUmH2w1kh9zHlGgqSUgjTdB72fQbBrVsdIFfId1nE+ICd9obyyX9erEmae3Dl4oyh9KuT4nh1xgQkaGDaUYCY85jwjTUJbmeNuWbE7N5X4GzabxdQxVmEgRA1nK5RlLAz/4p6QII/JzWLblIOt3BQwHjhhB6pYtjP/tb3n35EnW9+8fTG8YsaaudLvVG/AP9W0N6RNhi2VK9vqyNMeTZM0Q7XI/XT5rdV2ddhzSY5KpG6960Mneg8dOap+HX9WHn1/tpN+9e7emp6frD3/4Qye9YUQSTrEYVFt/ywNCQLa/FQH2eFSDxgyZxKMHEDQOkCi4JHR1mVBTPcvykLYj86H76L14Bhw8GNje9pkt+N5FXXl++Q5KTlYE1nfu3Jk777yTJ598kmPHLHWckcDU5VVqbsACoG3Y+7bAgsbq47klUs8jXj0AkdrtEonudU+VuPTkli3zLvL//p+bfMsBzf3Zy/r0x1866d977z0FdMqUKU56w4gURKgM7RlAWdj7Mr/NCKOhNUjx6gEka1btuGR96N8fBg2CiRNrL/bTAJf0aM/5XU9j+sIt1Q9agbjyyiu54IILEj9Vu9GsCeI8ngIWi8ivReTXwCLgyahYlcQ0NGQSr7o6yZp6J25rukIh+PRTePfdwFIRobAgl0+/OsryrcGHvkSEoqIili1bxpIlSwLrDSMWBFokKCL9gSv9twtUdUVUrIowibRIMJ7JN4uLG19hs9lz4gRkZ8OQITBrVmD5sZMV5P/2bb57/hn81739AusPHz5Mt27duPfeexM7XbvRpInIIkHxiiyfD7RT1f8G9ovIoAjZ2GyIdQ8gPDj/8MPedZpJ6p1To1UreOABeP55L+NuQFpnpHFH/2xeWb2LA8fKGhbUoF27dowYMYKZM2dy0CFwbxjRJsiw1Z+BS4Hh/vujwJ8iblETJ5bDMHWtkfj+95NvvUdcKCqCigqYOtVJXliQS1llFc8tdSsUFQqFOHHiBE8+aaPDRuIRpBjUclXtLyIrVPUSv22Vql4cVQsjQCINW8WSuobIRL4ZB25svZt4ErchtyFDYONG+Pxzr8Z1QO6Z+DG7j5Yy/8dXk5JST7W5OigoKODQoUOsX78eqa9anWFEgUjltioXkVS85J+ISBZQFQH7jChRVxC+5vNCoq/3iOuq/VDIu5GvvuokH1GQw5b9x3l/0z7Hy4fYsGED8+fPd9IbRrQI4jz+CDwPdBaRR4EPgN9GxSojIgSZhpvI2b/jusDxlluga1fnQlHXX9iFTq1bMH2hW+mbe+65hw4dOti0XSPhCFIMqhj4R+DfgF3Abar6XLQMM2onyOr02oLzdY18JPJ6j3hNbwa8muZjx8Lrr8MXXwSWZ6Slcs/AHry9fjc7D50IrG/VqhWjR49m7ty57NoVMFuvYUSRQIkRVfVTVf2Tqv6vqq6PllFG7QQdvqktOF9UlHzrPeK+wHHsWM9bT5rkJL9vUA4KzFzs5u2KioqoqKhgypQpTnrDiAp1LT2vuQECFAK/8t/nAIMaq4/nFsv0JNEkUqk6qpMAing/Ez05YkIkdbztNtXTT1ctLXWSPzBtkQ781ze1rKLSST9kyBDt3r27lpeXO+kNwwUilJ7EpurGmUgN3yRbqeWEqBwaCsG+ffCXvzjJCwty2XP0JG9+stvx8iG2b9/Oq46Be8OINEGcR76q/gAoBVDVg/jVMhuDiPQQkfki8omIrBORv/Pb+4nIQhFZKSJLqxceiscfRWSTiKz2V7dXn2uUiGz0t1EBPkNSE/fhmzgSd4c3ZAj06uUcOL/6nM5kt2/lHDi/5ZZb6NatmwXOjYQhllN1K4Afq+r5QAHwAxE5H/g98M+q2g/4lf8e4Aagj7+NA6+8m4h0BB4B8oFBwCMi0iGAHUlLIuenck0znzQFqlJSYPx4+OADWLMmsDw1RbgvP4ePNu9n896SwPq0tDTGjh3LvHnz+PzzzwPrDSPSuEzVPcNlqq6q7lLV5f7ro8B6vLogCpzmH9YO2Om/vhV4yh96Wwi0F5GuwFDgTVU94Pd+3gSuD/A5kpaEGL6pBdd1GIlcdbFWRo+GjAwv264D9+T1ID1VKA5aptZn7NixpKSkMMkxcG8YkSRoYsRzgcH+23fUccaViPTEqw9yIZ4DmYcXkE8BLlPVLSLyMvDvqvqBr3kb+BlwNdBSVf/Vb/8lcEJV/2+Na4zD67GQk5MzYEttS62NiOCa7DGeSSKduf9+eOEF2LkT2rQJLP+bZ5az4LO9LPrFEFq1CL5i/Y477uD9999n+/btZGRkBNYbRhAilRixJXAjMAS4FrjebwtqTIoefGUAACAASURBVBvgL8CPVPUIXnXCv1fVHsDfA26JhGqgqpNVNU9V87KysiJxSqMOXAP5cV2/4UooBEePOnePCgtyOVJawUurdzZ8cK2XD7Fv3z5mz57tpDeMSBG0nscFeMNX/4uXYffpIBcTkXQ8x1GsqnP85lFA9evn8OIYADuAHmHy7n5bXe1GnHAN5CflBIBLL4WLL/YC5w6FnvLP7Ejvzm0odgycDx48mN69e1vg3Ig7QZzHhar6kKrO97exeM6kUfgp3acC61X1D2G7dgJX+a+vBTb6r18ERvqzrgqAw6q6C2+I6zoR6eAHyq/z24w44RrIT+QJAHUi4q20XLUKFi1ykAsj8nNYtf0wa7YfDqxPSUmhqKiIDz/8kNWrVwfWG0akCOI8lvtf4gCISD4QJFXt5cD9wLX+tNyVInIjMBb4TxFZhReAH+cf/yrwObAJeAz4PoCqHgD+BVjib7/x24w44RrIT9QJAA0yYoQX73B8+r+jf3dapac6T9t94IEHyMjIYKJj4N4wIkGQlOzrgXOA6hHpHGAD3hRcVdWLomJhBGiuKdmNKPL978O0abBjB3TqFFj+s9mreWHVDhb9YgjtWqUH1o8aNYo5c+awc+dO2rZtG1hvGI0hUinZrwfOxBtiusp/fT3wPeDmUzXSMJKKUAhOnoQnnnCSFxbkUlpexZzl2x0vH6KkpITihJ3XbDR1gjiPQcABVd2CN/z0X0AnVd3itzV7kmbBm3Hq9O0Ll1/urfmoCl7Wpm/3dlzcvR3Fi7YSZLp8Nfn5+fTr148JEyY46Q3jVAniPH6pqkdF5Aq86bpT8Vd9G/Fd8GZOq26iem9CIdi0Cd5+20k+oiCXTXtKWPh58JCdiBAKhVi9ejUff/yx0/UN45SoK2NizQ1Y4f/8N+C+8LZE32KRVTdSGW+DkhAZZxOUqN+b0lIv0+7ttzvJj5+s0L6PvK7fL17mpD969Ki2bdtWCwsLnfSG0RBEKKvuDhGZBNwLvCoiGQSsB9KUideCt7hW2Utwon5vMjLgwQfhxRe9wHlAWrVI5a4BPZi39iv2HC0NrG/Tpg0jR45k1qxZ7NvnVubWMFwJ8uV/D956iqGqegjoCPw0KlYlIa4L3k51WCUpV2nHiJjcm/HjvZjHY485yUcU5FBRpcxass1JHwqFKCsr4/HHH3fSG4YrQcrQHlfVOaq60X+/S1XfiJ5pyYXLgrdIxEmScpV2jIjJvTnrLBg61HMe5eWB5b2y2nBZr07MWLyNyqrgge8LLriAK6+8kkmTJlHlELg3DFds2ClCuCx4i8SwSlKu0o4RMbs3oZCXKPGll5zkhQW57Dh0gvmf7nG8fIjNmzfz5ptvOukNw4m6giFNaUvUMrQitQfZRYKdJ9nKysaSmNybigrVHj1UhwxxkpdVVOrAf31TR01b5KQvLS3VrKwsvfXWW530hlEXRChgbkSYSA2rxL3KXgITk3uTmuqNN771Fmzc2PDxNUhPTWHYwB6899leth043rCgBhkZGTz00EO89NJLbNvmFjsxjKAESckuIlIoIr/y3+dUl4w13LAhpybEmDGQluZcKGrYoBwEKF7kFs0fP348qspjjoF7wwhKkJ7Hn4FLgeH++6PAnyJuUTMiaRMDGt+mSxe4/XZ4/HE4cSKwvFv7Vgw+7wxmLd3GyYrKwPqePXtyww03MGXKFModAveGEZQgziNfVX8AlAKoVwK2RVSsakbYkFMTIhSCgwdh1iwneWFBLgeOlfH62q+c9EVFRezatYsXXnjBSW8YQQjiPMpFJBWv5jgikgXY3EDDqObqq+Hcc51TtV/Z+3RyO2U6p2q/8cYbycnJsVTtRkwI4jz+CDwPnCEijwIf4NXfMIyEJaZ5v6oLRS1aBCtWBJanpAj3DcphyZcH+fSrI4H1qampjBs3jrfffpvPPvsssN4wghBkkWAx8I94DmMncJuqPhctwwzjVIlLsspRo6BVK+fex915PWiRlsIzjoHzhx56iLS0NOt9GFEnyGyrDKA/0A7oBNxdPfPKMBKRuOT9at8ehg/3PNTh4GVmO7ZuwU19uzJn+Q6OnawIrO/SpQt33HEHTzzxBCccAveG0ViCDFu9ANyKVznwWNhmGI0i1qnj45b3KxTyvNTTTzvJCwtyKDlZwQsrdzpePsTBgwd59tlnnfSG0RiCOI/uqnqvqv5eVf+zeouaZUmO1dj4JvEYQopmbqt6f795ed42YYL3YQPSP6cD53Zpy/SFW1AH/VVXXcV5553HBMehM8NoDEGcx0ci0jdqljQh4lkYKlGJxxBStBZh1vb7vf9+L17+V0cSCsEnn8D77wc+v4hQWJDLJ7uOsGLbISd9UVERixcvZvny5YH1htEYgjiPK4BlIrJBRFaLyBoRWR0tw5IZq7HxbeIxhBStRZi1/X6rOwjVDwozGebFPxyf/m+7JJvWLVKdp+2OHDmSzMxM630YUUMa2y0Wkdza2jUJ6pfn5eXp0qVLY3a9lJTaRytEnMpdNwl69vS+WGuSm+stjkwm6vr9hpObC1/e9iP4859h2zY444zA13n4+TU8t2w7i34+mA6tg6/HHTNmDDNmzGDnzp20a9cusN4wRGSZqubVti/IVN0ttW2RM7PpYDU2vk1TyuPVmN/j1q14az7Ky2HaNKfrFBbkUlZRxexl2530oVCI48eP89RTTznpDaM+GnQeIvKB//OoiByp+TP6JiYfTemLMlI0pTxetf1+a5KTg7fa/OqrYdIkqAyer+q8rqcxILcDxYu2UOVQKGrAgAEMHDiQCRMmOAXeDaM+GnQeqnqF/7Otqp5W82f0TUw+mtIXZSRpKnm8wn+/4P2Ow/nGg0Io5I3Xvf6607UKC3L4cv9xPtzsVqM8FAqxfv16FixY4KQ3jLoIEvP4h1qaDwPLVHVlRK2KMLGOeRjNi+JiL4i+davX43j00TDHWFbmNeblwcsvBz53aXkll/7b2ww6syOT7q916Llejh8/TnZ2NkOHDmXmzJmB9UbzJiIxDyAPKAKy/W08cD3wmIj84ylbaRhJSr09qhYtvFofr77qNDOgZXoq9+T14K31e/jqcGlgfWZmJqNGjWLOnDns3r07sN4w6iLQIkGgv6r+WFV/DAwAOgPfAR6Igm1GAmOLIAMwbpw3tjV5spP8vvwcKquUGYvd5jUXFRVRXl7O1KlTnfSGURtBnEdn4GTY+3LgDFU9UaPdiACJ/OVsiyADkpMDN90EU6d6w1gBye3Umu+cncXMJVsprww+1/vcc8/lmmuuYdKkSVQ6BO4NozaCOI9iYJGIPCIivwY+Ap4RkdbAJ9EwrrmS6F/OtgjSgVAI9uyBOXOc5IX5Oew+cpK317sNPYVCIbZu3cprr73mpDeMmjQ6YA4gInnA5XgFoT5S1aSIQidbwDzRF9TZIkgHqqqgd2/o0QPeey+wvKKyiit/P59eWW2YPiY/sL68vJycnBz69+/PK6+8ElhvNE8iEjD3U7KfDbQG2gM3Wkr26BC3bLCNxBZBOpCSwopB42HBAi6UdYGHItNSUxg+KIcPNu3ji33Bk1mnp6czZswYXnvtNb744ovAesOoiaVkT0AS/cvZFkEGp7gYbnvxQU7SgvFMdBqKHDawB2kpQrFjvqtx48YhIkx2DNwbxjdQ1UZtwNrGHpto24ABAzSZmD5dNTNT1Rsc8rbMTK89UZg+XTU3V1XE+5lItiUiubne73E69+khTtNMShS89iCEpi/Vi349T0+UVTjZccstt2hWVpaWlpY66Y3mBbBU6/hejVlKdhHpISLzReQTEVknIn8Xtu9vReRTv/33Ye0/F5FNfibfoWHt1/ttm0Tkn1xtSlSSYYV6U1ktHiuqhxwnEKIdRxjOjG+0N5bC/FwOnyjn5dW7nOwIhULs3buXOY6Be8OoJsgK80+APsDneFNzBVBVvaiR+q5AV1VdLiJtgWXAbcAZwMPATap6UkQ6q+oeETkfmAEMAroBb+HFXAA+A74LbAeWAMNVtc4ZX8kWMDeaHl9PglBWcxFltCCPpeTmSqBJEKrK4D+8x2kt05n7g8sD21FVVUWfPn3Izs62lCVGg0RqhfkNQG/gOuBm4Hv+z0ahqrtUdbn/+iiwHm+legj4d1U96e/b40tuBWaq6klV/QLYhOdIBgGbVPVzVS0DZvrHGkZUOZW1N1/HiYQJhBjAcq7MWBI4TiQijMjPZeW2Q6zdEbxGekpKCuPHj+f9999n7dq1gfWGUU0Q57EVuBIYpV4qdsXrNQRGRHoClwCL8HoTV4rIIhF5T0QG+odlA9vCZNv5OjVKbe2GETUaVT2wHsKHIosp5Ji05vFBE5yG++7q352W6SkUL3ILnI8ePZoWLVowadIkJ71hQDDn8WfgUmC4//4o8KegFxSRNsBfgB+p6hEgDegIFAA/BWaJ1MxTGhwRGSciS0Vk6d69e0/1dEYzpzHVAxvjQL78Eg7rabQeV0ivJTPhwIHAtrTLTOfmi7rxwsqdHCktD6zPysri7rvv5qmnnqKkpCSw3jAgmPPIV9UfAKUAqnoQCFTeTETS8RxHsapWR+y2A3P84P5ioAo4HdgB9AiTd/fb6mr/Bqo6WVXzVDUvKysriJmG8S0aCmwHXmEfCkFpKTz5pJM9hQW5HC+rZO6Kb/3pN/LyIY4cOcKMGTOc9IYRxHmUi0gq3nAVIpKF90XfKPzexFRgvar+IWzXXOAa/5iz8RzSPuBFYJiIZIjImXjB+sV4AfI+InKmiLQAhvnHGkbUaHT1wMZy8cVw6aUwcWLDNW1rk/doT9/sdkxfuMWp0NNll11G3759rVCU4UwQ5/FH4Hmgs4g8CnwA/DaA/nLgfuBaEVnpbzcC04CzRGQtXvB7lN8LWQfMwsub9TrwA1WtVNUK4G+AeXhB91n+sYYRNRpdPTAIoRB89hm8846TTYUFOXy2u4QlXx4MrBURQqEQK1asYPHixU7XN5o3QXNbnQsMxpum+7aqro+WYZHEpuoakaC66NOWLV6gPPxfJzPTYS1OaSl07+6Vqp09O7A9x8sqyP/t21xzTmf+OPySwPqjR4/SrVs37rzzTp544onAeqPpc0pTdcOD16r6qar+SVX/N9xxRCLAbRjxorFTcKsD3qrw9NMRWMTZsiWMHg1z58LOnYHtzmyRxp39u/Pa2l3sKwleFaFt27YUFhby7LPPcsAhcG80bxozbDXfXwH+jU65iLQQkWtF5ElgVHTMM4zo4pr+3mWFfa1Oavx4qKyEKVOc7C8syKG8Upm1dFvDB9dCKBSitLTUeh5GYBocthKRlsCDwAjgTOAQ0ArP8bwB/FlVV0TZzlPChq2MuohV+vtqJxU+3fevQ11PDYV167wLpqUFPvewyR+z/eAJ3vvpNaSmBB8EuPzyy9m7dy+ffvopKSlBwqBGU+eUhq1UtVRV/6yqlwO5eDGPS1Q1V1XHJrrjMIz6iFX6+3oLaBUVwY4d4Fhno7Agl+0HT7DgM7f1TKFQiI0bN/KOY+DeaJ4EesxQ1XI/zcihaBlkGLEkVunv63VSN98M2dkwYYLTua87vwunt8lgumOq9rvuuotOnToxwfH6RvPE+qhGsyZWtUnqdVJpaTB2LMybB5s3Bz53i7QUhg3swTsb9rD94PGGBTVo2bIlDz74IC+88AI7drgtOjSaH+Y8jGZNrNLfN+ikxoyB1FRwzDc1PD8HAWYsdhtvGz9+PJWVlUxxDNwbzQ9zHkazJxa1SRp0UtnZcOutMG2at/4jINntW3HtuZ15dsk2yiqCF5Lv1asXQ4cO5bHHHqOioiKw3mh+BHYeIvJdEXlMRPr578dF3izDaHo06KRCIdi/32nBIMCIglz2lZTxxidfOelDoRA7duzgpZdectIbzQuXnseDeNlvC0XkWqBfZE0yjOTlVGp+cO210KePc+D8qj5Z9OjYyjlwftNNN9G9e3cLnBuNwsV5HFXVQ6r6E7zCUAMbEhhGc8B1weFfSUnxpu1+9BGsXh34+ikpwn2Dcln4+QE27TkaWJ+Wlsa4ceN488032bhxY2C90bxwcR5/nYyuqv8EPBU5cwwjeal3LUdjeeABL22J49P/PXndSU8Vpi90C5yPGTOGtLQ0KxRlNEijnYeI/LeIiKq+EN6uqv8TebMMI/mIyILDjh3h3nth+nQ4Grz30KlNBjdc2JW/LN/O8bLgge+uXbty22238fjjj3PixInAeqP5EKTncRR4UUS8SswiQ0Xkw+iY1bw5pXFzI25EbMFhKAQlJZ4DcaCwIJejpRW8tCp4skXv8iEOHDjAc88956Q3mgeNdh6q+n+AGcB7vtP4B+CfomVYc+WUx82NuBGxBYeDBsEll3hDVw6Fmgb27MDZZ7RxHrq65pprOOeccyxwbtRLkGGrwcBY4Bhemdgfqur70TKsuRKRcXMjLkRswaGI1/tYswY+DN65FxEKC3JZs+Mwq7YFzyQkIhQVFbFw4UJWrlwZWG80D4IMWz0M/FJVrwbuAp71p+oaESRWifqM6BCxBYf33QenneYcOL/9kmwyW6Q6T9sdNWoUrVq1st6HUSdBhq2uVdUP/NdrgBuAf42WYc2VWCXqMxKc1q1h5EhvweDe4Nly27ZM59Z+2by0eieHj5cH1nfo0IFhw4ZRXFzMkSNHAuuNpo9zehJV3YWXnt2IILFK1GckAUVFUFYGjz/uJC8syKG0vIrZy7c76UOhEMeOHePpp5920htNm1PKbaWqNpcvwsQqUZ+RBFxwAXznO16yxKrg+aou6NaOS3LaU7xoCw0VfauNgQMHMmDAACZMmOCkN5o2lhgxAYlFoj4jSQiF4PPP4Y03nOSF+bl8vvcYH2/e73j5EOvWreODDz5w0htNF3MehpHI3HEHdO7sHDi/6aKutM9MZ/oit8D5sGHDaNeunQXOjW9hzsMwEpkWLeChh+Dll52m3LVMT+XuAd15Y91u9hwJnuq9devWjBw5ktmzZ7Nnz57AeqPpYs7DMBKdceO8xYKPPeYkvy8/l4oqZeaSbU76oqIiysvLedwxcG80Tcx5GEai07Mn3HgjTJkC5cGn3Z55emuu7HM6MxZvpaIyeOD9/PPP56qrrmLSpElUOQTujaaJOQ/DSAZCIfjqK5g710k+Ij+XXYdLeedTt6GnUCjEF198wbx585z0RtPDnIdhJAPXX+/N23YMXA85rzNdTmvJ9EVuqQpuv/12zjjjDAucG3/FnIdhJAOpqTB+PMyfD59+GlielprCsEE9WPDZXrbsPxZY36JFCx566CFeeeUVtlquHANzHoaRPDz0EKSnw8SJTvJhA3NITRGecex9jBs3DlVl8uTJTnqjaWHOwzCShc6d4c474cknv516uRF0adeS7553BrOWbqO0vDKwPjc3l5tuuokpU6ZQVlYWWG80Lcx5GEYyEQrBoUMwc6aTvLAgl4PHy3lt7S7Hy4fYvXs3cx0D90bTwZyHYSQTV14J55/vHDi/rFcnzjy9tXOhqKFDh9KzZ08LnBvmPAwjqRDxsu0uXeptAUlJEUbk57Bsy0HW7wqeaj01NZXx48fz7rvvsn79+sB6o+lgzsMwko2RI708/Y5P/3cN6E5GWopzoagHH3yQ9PR0JjoG7o2mgTkPw0g22rXzKg3OmAEHDwaWt89swfcu6sbcFTsoOVkRWN+5c2fuuusunnzySY4dCz7t12gaxMx5iEgPEZkvIp+IyDoR+bsa+38sIioip/vvRUT+KCKbRGS1iPQPO3aUiGz0t1Gx+gyGkTCEQnDiBDz1lJO8sCCHY2WVPL9ih+PlQxw+fJiZjoF7I/mJZc+jAvixqp4PFAA/EJHzwXMswHVAeBTvBqCPv40DJvjHdgQeAfKBQcAjItIhVh/CMBKC/v1h0CBvzYdDoaZ+PdpzQbfTKF7oVijqiiuu4IILLrDAeTMmZs5DVXep6nL/9VFgPZDt7/4v4B+B8L/iW4Gn1GMh0F5EugJDgTdV9YCqHgTeBK6P1ecwjIQhFPJWm7/7bmCpiFBYkMunXx1l2ZbgQ18iQlFREcuWLWPJkiWB9UbyE5eYh4j0BC4BFonIrcAOVV1V47BsIDyH9Ha/ra72mtcYJyJLRWTp3r17I2i9YSQI994LHTo4B85v7deNthlpzoHz+++/n8zMTOt9NFNi7jxEpA3wF+BHeENZvwB+FenrqOpkVc1T1bysrKxIn94w4k+rVvDAA/D887Ar+KK/zBZp3NE/m1fXfMX+kpOB9e3atWPEiBHMnDmTgw6BeyO5ianzEJF0PMdRrKpzgF7AmcAqEfkS6A4sF5EuwA6gR5i8u99WV7thND+KiqCiAqZOdZIXFuRSVlnFc8u2O+lDoRAnTpzgySefdNIbyUssZ1sJMBVYr6p/AFDVNaraWVV7qmpPvCGo/qr6FfAiMNKfdVUAHFbVXcA84DoR6eAHyq/z2wyj+XH22TB4MEyeDJXB81X1OaMtg87syDOLtlJVFTxwfskll5Cfn8/EiROdAu9G8hLLnsflwP3AtSKy0t9urOf4V4HPgU3AY8D3AVT1APAvwBJ/+43fZhjNk1AItm2DV15xkhcW5LL1wHEWbHSLDYZCITZs2MD8+fOd9EZyIs3haSEvL0+XOqRyMIykoLzcKxTVrx+8+mpgeVlFFZf9+9tcktOBx0bmBdafOHGC7OxsBg8ezHPPPRdYbyQuIrJMVWv9o7AV5oaR7KSnw9ix8Prr8MUXgeUt0lK4J68Hb6/fzc5DJwLrW7VqxejRo5k7dy67HAL3RnJizsMwmgJjx0JKCkya5CQfPigHBWYudsu2W1RUREVFBVOmTHHSG8mHOQ/DaAp07w433+zNujoZfNptj46ZXHNOZ2Yu2UZ5ZVVgfZ8+fRgyZAiTJ0+moiJ4viwj+TDnYRhNhVAI9u2Dv/zFST4iP4c9R0/y5ie7HS8fYvv27bziGLg3kgtzHobRVBgyBHr1cl5xfvU5nclu38p5xfktt9xCt27dbMV5M8Gch2E0FVJSvEWDH3wAa9YElqemCPfl5/DR5v1s3lsSWJ+WlsbYsWOZN28emzdvDqw3kgtzHobRlBg9GjIyvGy7DtyT14P0VKHYsUzt2LFjSU1NZZJj4N5IHsx5GEZTolMnuOceePppKAnee8hqm8HQC7owe9k2TpQFX7GenZ3NLbfcwrRp0ygtLQ2sN5IHcx6G0dQIheDoUSgudpIXFuRypLSCl1bvdLx8iP379zN79mwnvZEcmPMwjKZGQQFcfLEXOHfIIJF/Zkf6dG5DsWPgfPDgwfTu3dsC500ccx6G0dQQ8Xofq1bBwoUOcmFEfg6rth9mzfbDgfUpKSkUFRXx0UcfsXr16sB6Izkw52EYTZERI6BtW+dpu3cM6E6r9FTnabsPPPAAGRkZ1vtowpjzMIymSJs2cP/9MGsW7N8fWH5ay3Ru7deNF1bt4PCJ8sD6Tp06ce+99zJ9+nSOHj0aWG8kPuY8DKOpEgp5qUoef9xJXliQS2l5FXOWuxeKKikpYfr06U56I7Ex52EYTZULL4QrrvDWfFQFz1d1YXY7Lu7RnuJFW50KPeXn59OvXz8mTJhghaKaIOY8DKMpEwrB5s3w1ltO8sL8HDbtKWHh58HrrYkIoVCINWvW8NFHHzld30hcmkUxKBHZC7hF/urmdGBfhM8ZDczOyGJ2RpZksDMZbITo2Jmrqlm17WgWziMaiMjSuipsJRJmZ2QxOyNLMtiZDDZC7O20YSvDMAwjMOY8DMMwjMCY83BncrwNaCRmZ2QxOyNLMtiZDDZCjO20mIdhGIYRGOt5GIZhGIEx52EYhmEExpxHIxCRHiIyX0Q+EZF1IvJ3fvt/iMinIrJaRJ4XkfaJaGfY/h+LiIrI6Yloo4j8rX8/14nI7+NlY312ikg/EVkoIitFZKmIDIqznS1FZLGIrPLt/Ge//UwRWSQim0TkWRFpkaB2FovIBhFZKyLTRCQ9Ee0M2/9HEQleZSvC1HM/RUQeFZHPRGS9iPwwakaoqm0NbEBXoL//ui3wGXA+cB2Q5rf/DvhdItrpv+8BzMNbLHl6otkIXAO8BWT4+zon4r0E3gBu8NtvBN6Ns50CtPFfpwOLgAJgFjDMb58IhBLUzhv9fQLMSFQ7/fd5wNNASTxtbOB+jgaeAlL8fVH7P7KeRyNQ1V2qutx/fRRYD2Sr6huqWuEfthDoHi8boW47/d3/BfwjENcZEvXYGAL+XVVP+vv2xM/Keu1U4DT/sHaAW7m9CKEe1U/C6f6mwLVAdSm/J4Hb4mDeX6nLTlV91d+nwGLi/z9Uq50ikgr8B97/UNyp5/ceAn6jqlX+cVH7PzLnERAR6Qlcgufpw3kQeC3W9tRFuJ0iciuwQ1VXxdWoGtS4l2cDV/pDLe+JyMB42hZODTt/BPyHiGwD/i/w8/hZ5iEiqSKyEtgDvAlsBg6FPdhs5+uHiLhR005VXRS2Lx24H3g9XvaF2VKbnX8DvKiqu+Jr3dfUYWcv4F5/SPU1EekTreub8wiAiLQB/gL8SFWPhLU/DFQAbkWjI0y4nXh2/QL4VVyNqkEt9zIN6IjX9f4pMEtEJI4mArXaGQL+XlV7AH8PTI2nfQCqWqmq/fCe2gcB58bZpFqpaaeIXBi2+8/AAlV9Pz7WfU0tdn4HuBv4n/ha9k3quJ8ZQKl6aUoeA6ZF6/rmPBqJ/2T0F6BYVeeEtT8AfA8Y4Xe940otdvYCzgRWiciXeH9oy0WkSwLZCN7T8Ry/O74YqMJL9BY36rBzFFD9+jm8L+uEQFUPAfOBS4H2IpLm7+oO7IibYTUIs/N6ABF5BMgC/iGedtUkzM5rgN7AJv9/KFNENsXTtnBq3M/tfP33+TxwUbSua86jEfhPwFOB9ar6h7D26/HGQG9R1ePxsi/Mnm/ZqaprVLWzqvZU1Z54f1z9VfWrRLHRZy7ePykicjbQgjhmMq3Hzp3AkIg6hwAABnJJREFUVf7ra4GNsbYtHBHJqp7lJyKtgO/ixWfmA3f5h40CXoiPhR512PmpiIwBhgLDq8fp40kddi5T1S5h/0PHVbV3Atr5KWH/R3h/p59FzYYEeFhOeETkCuB9YA3eEzF4Q0F/xOsmVtf5XKiqRbG30KMuO1X11bBjvgTyVDUuX8z13Mu38LrY/YAy4Ceq+k48bIR67TwC/DfeMFsp8H1VXRYXIwERuQgvIJ6K9zA4S1V/IyJnATPxhgJXAIXVkxESzM4KvBmA1bVq56jqb+JkZp121jimRFXbxMO+MBvqup/t8YbPc4ASoChasU5zHoZhGEZgbNjKMAzDCIw5D8MwDCMw5jwMwzCMwJjzMAzDMAJjzsMwDMMIjDkPwzAMIzDmPAzDMIzAmPMwooaItPKTHKaKSHsR+X68bWqIWNgpIh9F+fz11puI9PVF5Nci8pNInvNUEZEWIrIgLEWLEWHMeRjR5EG8FcOVQHsgIZyHXzCnrr/9wHY2cL5voaqXBTl/pIn39esi6H2sD1UtA94G7o3E+YxvY87DaBDxKup913/9ryLS2OyiI/g6p9K/A73Eq8D3HyJSKF4ltJUiMsmvl4CI9BSvmuATfjW0YhEZIiIfishGERkUdkyxeNXSZotIZpi93zq3r9kgIk8Ba4EeIjJXRJaJV4ltXB129hSRtWHn/on/pF3b+Wr9TLXcz5Kwz7peRB7zbXjDz1NU8/iR4lWrXCUiT9f3OcP2tRaRV3zNWhG5N2xfSV2fqyFt2PEP+7+fD4Bzauyr63f7S/+efSAiM/xrNvo+1vF7rc/WuXh/g0Y0cKkgZVvz2oDvAO/i/SO+AqQ2QtMC+CrsfU9grf/6POAlIN1//2dgZNhxFUBfvIebZXg5rwS4Fe8LoSde4ZvLfc00vFxYdZ7b11ThV4Xz93X0f7bC++LqFG5nTbv99z8Bfl3zfPV9plruTUmNz9rPfz8LLwdV+LEX4CW3O72GzfXdwxLgTuCxsPO0C79+XZ/Lf12n1n8/AC/nVyZeYaxNjbj/A4GVQEu8yowb/Ws26j7W017f50wF9sb7/6epbjYeaDSIqi4QEcFLmX21qlaKl3jvYbx/1rtqkZ0OHKrjlIPxvoCWeKelFV5Bm2q+UNU1ACKyDnhbVVVE1uB92QBsU9UP/dfTgR/iFWeq69wLgC2qujDsOj8Ukdv91z2APkCQbMPh52voM9XFF6q60n+9LOzzVXMt8Jz6iSxV9UAjr7cG+E8R+R3wsgark9GQ9krgefUzSYvIi2H76rKrI/CCqpYCpSLyUpimMffxtDran6nLVv/vtExE2qpXDdKIIOY8jAYRkb54Nb33V/8TqurnwEMiMrsO2Qm8p8xaTwk8qap1VeELz/5aFfa+iq//Zmtm9Kx+X+u5xasGeCzs/dXAEOBSVT0uIu/WYW8F3xzeDT/mWNjrhj5TXYR/1kq8L8XGUO/1VPUzEemPVyP8X0Xkbf1mdtg6P1cjtIHtEpEf1aNp8D6KyN/W1u7vq8/WDLzsx0aEsZiHUS8i0hUvxfOtQIl4NUwaRFUPAqkiUv2ldBRvuAK8QOZdItLZv0ZHEckNaFqOiFzqv74P+CDgudsBB33HcS5eBcOadgLsBjqLSCcRycAr/FUbkfhMtfEOcLeIdKo+b2OuJyLd8OpOTMervd2/xnnr/FyN0C4AbhNvNl1b4OawfXXZ9SFws4i0FK86Y9D7WGt7fbb692yfqpbXcS3jFLCeh1En4gWh5wA/VtX1IvIvwO9ofJ3pN4ArgLdUdb94Qe+1eLXe/w/whniza8qBH+DVdWgsG4AfiMg04BNgAoCqfiIitZ275nDU60CRiKz3z7XQ13/DTlX9qYj8BliMV43v09qMqee6QT5TbeddJyKPAu+JSCVebY4HGnG9vni11qv8faEa5y2v53M1pF0uIs8Cq/CGjpY0dB9UdaE/vLUaz3GtAQ7X8nnr09f2edvVY+s1eDE6IwpYPQ/DCf+p7lG8CmZTVPXfajmmP1697/sjfO2eeOPbFzZwqJFAiEgbVS3xH0oWAONUdXkUrzcH+CdVjVo1veaM9TwMJ1R1P1Bv1UT/CXW+iKT+//bt2AhgEAaCoJRQjEt25g7oEAIq+IDIu03caH60zq8H//Z291NnX/kuh2NU1RSOe1weAMQM5gDExAOAmHgAEBMPAGLiAUBMPACIiQcAsQ0kTMPbOo+h7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Using Scipy\n",
    "#### SVM\n",
    "Luckily for us, we do not need to worry about in practice because Scipy packs an SVM model"
   ],
   "metadata": {
    "id": "5-PVaBijm26V"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "from matplotlib.axis import YTick\n",
    "from sklearn import svm\n",
    "\n",
    "clf_svm = svm.SVC(kernel=\"linear\")\n",
    "X = np.append(failure_features, normal_features, axis=0)\n",
    "#y is a vector of ones\n",
    "y = np.ones((2*no_samples,1))\n",
    "# make the normal features negative:\n",
    "y[no_samples:]=-1\n",
    "\n",
    "clf_svm.fit(X,y)\n",
    "\n",
    "xx = np.linspace(21, 37, no_samples)\n",
    "yy = np.linspace(2300, 3300, no_samples)\n",
    "YY, XX = np.meshgrid(yy,xx)\n",
    "xy = np.vstack([XX.ravel(), YY.ravel()]).T\n",
    "Z = clf_svm.decision_function(xy).reshape(XX.shape)\n",
    "\n",
    "#Prepare the figure \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "#Plot the data in failure state in red\n",
    "ax.scatter(temp_readings_failure, rpms_failure, color='red')\n",
    "#Plot the data in normal operation in blue\n",
    "ax.scatter(temp_readings_normal, rpms_normal, color='blue')\n",
    "\n",
    "ax.contour( XX, YY, Z, colors=\"k\", levels=[-1, 0, 1], alpha=0.5, linestyles=[\"--\", \"-\", \"--\"])\n",
    "\n",
    "plt.show()"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 300
    },
    "id": "iRX391FnnGTQ",
    "outputId": "7d97b7da-0bbe-47db-8d36-ffc9a8eb39d7"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "/usr/local/lib/python3.7/dist-packages/sklearn/utils/validation.py:993: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXRc5Znn8e9jWbIkL7ItC+NNls2S2AHjGLVxICE0Sdh7TE4PHWhDO8lh6EB6spDpTAdONx166DUnOck0JM0E5wTiDksgBNJuwBCzdbDB8oLxiowt71h4lS3LsqVn/ri3iiqVSmvt9fucU0el+96qeupW1XPf+77vfa+5OyIiUhyGZDsAERHJHCV9EZEioqQvIlJElPRFRIqIkr6ISBEZmu0AejJu3Divq6vLdhgiErFuHbS3Jy4vK4Pzz898PNKthoaGD9y9pruynE76dXV1rFy5MtthiEjEkCSNA6dOgX6rOcPMmpKVqXlHRPqutrZ/yyXnKOmLFJLFi6GuLqiR19UF/6fSffdBZWX8ssrKYLnkBSV9kUKxeDHcdhs0NYF78Pe221Kb+BcsgAcfhKlTwSz4++CDwXLJC5bL0zDU19e72vRF+qiuLkj0XU2dCtu3ZzoaySIza3D3+u7KVNMXKRQ7dvRvuRQlJX2RQqFOVukDJX2RQpGLnazp7liWflPSFykUudbJmomOZek3deSKSHqoYzlr1JErIpmnjuWcpKQvIumhjuWcpKQvIumRix3LoqQvImmSax3LAuT4LJsikucWLFCSzzGq6YuIFBElfRGRIqKkLyJSRHpN+mZWbmZvmtlaM1tvZt8Nly82s81m9o6ZLTKz0nC5mdmPzKzRzN42szkxz7XQzN4NbwvT97ZERKQ7fanpnwQud/cLgNnAVWY2D1gMfBQ4H6gAbg3Xvxo4J7zdBvwYwMzGAvcAFwFzgXvMbEzq3oqIiPSm16TvgWPhv6Xhzd19SVjmwJvA5HCd+cDDYdFyYLSZTQCuBJa6+0F3PwQsBa5K9RsSEZHk+tSmb2YlZrYG2E+QuFfElJUCtwDPhYsmATtjHr4rXJZsedfXus3MVprZyubm5v68FxER6UWfkr67d7j7bILa/FwzOy+m+AHgVXd/LRUBufuD7l7v7vU1NTWpeEoREQn1a/SOux8GlhE2y5jZPUANcGfMaruBKTH/Tw6XJVsuIiIZ0pfROzVmNjq8XwF8DthkZrcStNPf5O6dMQ95BvizcBTPPOCIu+8FngeuMLMxYQfuFeEyERHJkL5MwzAB+LmZlRDsJB5399+a2WmgCXjDzACecvd7gSXANUAj0Ap8CcDdD5rZ3wFvhc97r7sfTOm7ERGRHvWa9N39beDj3Szv9rHhaJ6vJilbBCzqZ4wiIpIiOiNXJNN03VjJIs2yKZJJkevGtrYG/0euGwuajVIyQjV9kUy6++4PE35Ea2uwXCQDlPRFMknXjZUsU9IXySRdNzZ/FUhfjJK+yGD1JxnourH5KdIX09QE7h/2xeRh4lfSFxmM/iaDwV43tkBqm3mngPpiLBhWn5vq6+t95cqV2Q5DJLm6uiDRdzV1KmzfntrX6jryB4KjBF1sPP2GDAl26l2ZQWdn4vIsM7MGd6/vrkw1fZHByGTHbAHVNvNOAfXFKOmLDEYmk4FG/mRPAfXFKOmLDEYmk0EB1TbzzmD7YnKIkr7IYGQyGRRQbTMvLVgQ9NN0dgZ/8zDhg6ZhEBm8BQsykwAir3H33UGTTm1tkPDzNPlIdijpi+STTO1gpGCpeUckXfJtTH2+xSsDopq+SDrk22ya+RavDJhq+iLpkG9j6vMt3mzL46Mi1fRF0iHfxtTnW7zZlOdHRarpi6RDvo2pz7d4synPj4qU9EXSId/G1OdbvNmU50dFSvoig9Vd+26+ncGZb/FmU54fFWmWTZHB0MyXxScPPnPNsimSLnnevlt0UjHqJs+PilTTFxmMPJtnvajlQQ09VVTTF0mXPG/fLSo6KgOU9EUGR6Ne8keej7pJFSV9kYFavPjD2mNJSbAsz9p3i4qOygAlfZGBib0gOkBHx4c1fCX83KSjMkBJX2Rg1D6cO/o6IifPR92kikbviAyERu3khiIakdMfGr0jkmpqH84NOuLqNyV9kYFQ+3BuyLUROXkw5bKSvshAqH04PfqbNHPpiCu2c9/9wymXcyzxK+mLDNSCBbB9e9CGv327En5f9JTUB5I0uzviMgsem+madp40NfWa9M2s3MzeNLO1ZrbezL4bLp9mZivMrNHMHjOzsnD5sPD/xrC8Lua5vhMu32xmV6brTYlIDuotqQ8kacYecUGQ8CMd7JmuaSdrUorsgHKkyafX0TtmZsBwdz9mZqXA68DXgTuBp9z9UTP7CbDW3X9sZncAs9z9K2Z2I/B5d/+Cmc0EfgnMBSYCLwLnuntHstfW6B2RAlJX9+F5DbGmTg2OlAY7Iqq350+3ZK8fuyOCjIwuGtToHQ8cC/8tDW8OXA78Klz+c+D68P788H/C8s+EO475wKPuftLdtwGNBDsAESkGvXW6DrZ9PtudusmamrruyLLc5NOnNn0zKzGzNcB+YCmwFTjs7qfDVXYBk8L7k4CdAGH5EaA6dnk3jxGRQtdbUh/siKhsd+p217mfrCUli/P99Cnpu3uHu88GJhPUzj+aroDM7DYzW2lmK5ubm9P1MiKSab0l9cGOiMqFYbRdO/cjfQ1dZfF8jn6N3nH3w8Ay4BPAaDMbGhZNBnaH93cDUwDC8irgQOzybh4T+xoPunu9u9fX1NT0JzwRyWV9SeqDGRHV9fmrq6GiAm65pW8dqOkYY58LO6Ku3L3HG1ADjA7vVwCvAdcBTwA3hst/AtwR3v8q8JPw/o3A4+H9jwFrgWHANOA9oKSn177wwgtdRKTffvEL98pK96CBJbhVVgbLU7F+f2OZOtXdLPh7++3x/6fiNboAVnqSvNqX0TuzCDpmSwiODB5393vNbDrwKDAWWA3c7O4nzawceAT4OHAw3DG8Fz7X3cCXgdPAN9z9P3t6bY3eEZEB6e9InkyN/MnQXEE9jd7RhGsiUnj6O/wzUxPoZWjnognXRKS49HckT6ZG/mR7WClK+iJSiPrbgZqpDtdsDytFSV9EClF/h39magK9HBjNozZ9ESk+kesb79gR1LIzeZnLDLy22vRFilUezO/eL6l4P71N/JbubZbl2VmH9r6KiOSlrsMDI8kN8nMa6FS9n95m8yykbdYNNe+IFKpszzqZaj29n/vu63uTSU/DM2trC2KbaZy+SDEqtIu3J3s/EHSG9vWEp2Q7j57k2TZTm75IMcqB4YEplSzukpL+XXyluxE0A33tPKSkL1KocmB4YEpdc01Q445VWQkdSa7DlOyEp65X2+pNPm+zbijpixSqQrp4++LF8POfxzfvmMHChQObvjgygqbrTiRWvm+zJNSmLyK5r7dO3IFOYlZond0htemLpEuhjYPPVT3NWTOYI5pCawLrA43TFxmoQhsHn8uSDaWMNOEsWDCwbR55TLbOzs2CnE76x44dY/PmzXzkIx9h165dPPnkk3HlZsZ1113H9OnT2bZtG7/97W8TnuPzn/88kydPZsuWLbzwwgvRx0XccMMNnHHGGaxfv55XXnkl7rkBbrrpJkaPHs3atWt54403Espvvvlmhg8fTkNDAw0NDQnlCxcupKysjBUrVrBu3bqE+G699VYAXn/9dTZv3hz32NLSUm655RYAXnnlFd5777248oqKCr7whS8A8NJLL7Fr1664566qquL664Pr1b/wwgu8//77ceXV1dVcc801ACxZsoSDBw/GbZvx48fz2c9+FoBnn32WlpaWuNefOHEin/70pwH4zW9+w4kTJ+LKa2tr+cQnPgHAU089xenTp6PPbWZMnz6dCy+8EIBf/epXCdvm3HPPZdasWZw6dYpnnnkm7rkBZsyYwapVM/jOd06wc+dz1NTAn/4pXHZZsM55553H2WefzbFjx3jppZfintvMmDVrFnV1dRw+fJjXXnst4fXnzJnDpEmTOHDgQPSzj339+r/6K8a3trIfiDZCtrZid94JY8cyd+5cqqur2bt3L2vXro17bYB58+ZRVVXFzp072bhxY0J8F198McOHD6epqYl33303Ib5PfepTDBs2jPfee49t27YlxHfppZcydOhQGhsb2blzZ9xzA1x22WUAbNmyhb1798aVl5SUcMkllwCwadMmul66tKysjIsuugiADRs2cPDgwbjnrqioYM6cOQC88847HD16NK58xIgRnH/++QCsW7eO48ePx5WPGjWKGTNmAPD222/T9uUvw9//PZw8iQGjgXPCGvnatWs5depUXHxjx45l+vTpAKxevZrOzs64bTNu3Dhqwx3G6pkz4de/jnv9mt27mTRpEh0dHbzzzjtxz21mjB8/nvHjx3Pq1Ck2bdpEVxMmTGDcuHGcPHmSxsbGhPKJEycyZswYTpw40e1nN3HiRKqqqmhtbWVHeJQTWz5p0iRGjBjBsWPH2LNnT0J8PcnppN/W1sa///u/M3r0aI4fPx5NerEOHz5MVVUVR48ejW68WC0tLYwcOZJDhw7RFFNTcHfMjOPHjzN8+HCam5vZsWNHwgZra2ujvLycffv2RX84seucPHmSsrIy9uzZw+7duyNXG4uuc/r0aUpKSmhqamLfvn3R1+76PFu3buX999+PLnN3hg4dSnt7O2bG5s2boz+8yDplZWW0trZiZqxfv56DBw9G3xcEP7wjR45gZqxZs4YjR47EvbeRI0dy4MABABoaGqJJPRLfmDFjojuKN954I5rUI6qrq6NfuFdffZWTJ08m7DSampowM1566aWEpD9p0iS2bt0KwPPPPw9AY6PT0ADHj0NV1UYWLtzCvHmnWLp0KV09+uhmfvObjbS3nwBeorkZ/vVfYfVqOOssaGxs5Oyzz6alpYVly5bFbRuAbdu2MW3aNA4fPszLL78cfe+RdXbu3MmUKVNobm7m9ddfjz4usn327NrFRGAP8PuYuHz/fnjsMfbv388ZZ5xBU1MTy5cvj/vsAA4cOMDYsWNpbGykoaGBrv1rR44cYdSoUWzcuJHVq1cnlLe2tlJZWcnatWt5++234747EHw3hw0bxltvvcWGDRsStl/E73//e7Zs2RK3rLS0lNOnT2NmvPzyywm/rcrKStrb2wF48cUX43YqEFQ4It+XJUuWRL/7EePGjYtWSJ5++mk++OCDuO0zYcIErr32WsyMxx9/PPjuzpsHK1dix49TW1XFlXfeCeecw+JFi2jtMmTznHPO4fLLLwfgZz/7WcJOYcaMGdEKy7/9278lbJMLLriAiy++mPb2dhYtWpRQXl9fz9y5czl+/DgPP/xwQvnFF1/M7NmzOXToEI8++mhC+WWXXcbHPvYx9u/fzxNPPJGQdz73uc/xkY98hN27d/PrLjskgGuvvZZp06axbds2/uM//iO6vLeEDzme9AFqamqYPHkyH3zwQbdJv6amhvHjx7Nv375uk351dTXV1dW4O9u76ZgZM2YMVVVVtLW1JfxozIyqqiqGDx8eTZ4RkXVGjhxJeXk5Bw4cSEjm7k5lZSVDhw6lrKws4bkBysvLcXdKSkoSXh+CH19k/e5ev7uYY//vCIezRS6VFrtOZ2cnJ0+eTCiPrNPZ2RmtgcWWRRJjR0cHR48e/fAybF3iO3XqFIcOHYp7TKy2tra4GmRjI/zXf0FHR7DekSOt/PjHezlwoINhwxLf+/PPt9Devgtoj5Z1dEBDQ5D0Dx06xPbt2xOOQCKPb25upqOjgxMnTnS7/fbs2UNra2vCNojuFKqqOHTkCC3Es+HDwYz33nuPffv2Jf3ubN68mfLycg4fPtxtRWDdunWUlZVFt2HX7bdq1SpKSkrijtBin2fFihUMGTKEDz74IOH13T26o2tubk7YobS3t7Ns2TIA9u/fH/0eRZ7nxIkT0aOn999/P+57Zma0tLREd9TNzc1xO3yAgwcPRnf0hw4dSijfv38/zz33HBDs/E6dOhW01YcjdXZXVLCksxOWLKGlpYXOzs6497B9+/bokX9ra2vC+9u6dWu0EhT5DcTatGlT9H2fPHky4fHr1q1j165dnD59Ovr9iv2ONzQ0sHXrVtrb2xN2SADLly9nw4YNtLW1JZSbGa+99hpr1qyJ+/7F+t3vfsfw4cOj5f0ZkJPTSX/48OFcf/310aQ/btw4ID4BzZs3jwkTJrB3717OPPPM6PLIOpdeeik1NTXs2LGDiRMnJrzGZz7zGcaMGcPWrVupra1N2HhXX301I0eOZNOmTdFD9Nh15s+fT0VFBevWrYseBsaW33DDDZSWlrJq1aqEQ3h35+abbwaCL0FjY2PcY4cOHcpNN90EBDXp7du3x5VXVFTwJ3/yJ0BQ24o070TWGTVqFH/8x38MxNe2IuXV1dXd1rYizjzzTK677joAHnvssYQjhdraWq666ioAHnnkkYQf19lnnx1tHvrpT38arW1F1pk5c2a0ieH+++/n6adjh1w7MJtTpz7Ja6+185d/+WDCZ/Ozn80FLgKOAYuijzt+HC64IGj+mDNnDgcPHoyrjUWe57Of/Sznn38+77//PotjOmAj5ddeey0f/ehH2blzJ4899hhdff688zjrr/+axtZWfh1GTGkp/NEfwYwZ3HjjjdTW1rJhwwaeffbZuOcG+OIXv8iZZ57J6tWrowkudp0///M/p7q6muXLlyc0TwF87WtfY+TIkbz66qu8+uqrCdvn29/+NsOGDeOFF15g+fLlCTvuv/3bv8XdefbZZ+OaJoO3Ucpdd90FwBNPPMH69evjykeMGMGdd94JwC9/+cuEI4Xq6mruuOMOAB5++OHoUXYkhgkTJnDrrbfi7vz0pz+NOxJwd6ZOnRpt2nzggQeizUcRZ511VrRp84c//CHHjh2LK585cybz588H4Hvf+170qCRi9uzZXH311bg7//AP/0BXc+fO5fLLL6etrY0f/OAHCeWXXHIJn/zkJ2lpaeH+++9PKL/sssuYO3cuBw4c4KGHHkr4bK688kpmzZrFvn37eOSRR+K2DcB1113HjBkz2LFjR9x3L7LO/PnzOeuss2hsbOTpp5+Oe+5eL4GrIZsfyuZsq4Wov9uzv7MG5MRouwL70nQ92oss65onIkenEBxNdpdHhg4N6pSnT5+ms8sHaGbRo9j29vaE1xgyZEj06LitrS1aFvlbUlLCsGHDgO5r8iUlJZSXlwNB32DX8tLS0mh55Gg1VllZGRUVFbh7tLITu055eTkVFRV0dnZy+PDhhG1TWVlJRUUFHR0d0SO1WMOHD6eiooJTp07FHelFjBo1ivLyctrb2zl06FBC+ejRoykvL6etrS3u+SPrTZo0KemQzW6vlp4rtwsvvNAz5Re/cK+sdA/STnCrrEzLheoHFNvUqe5mwd9ciKk3A9meU6fGrx+5TZ2autcQKQbASk+SV7Oe2Hu6ZTLp9zfhZEq+JraBbM+BvNd83CGKpFtPSV/NO6G+NC1k40g+J5owBmCgEzwWWGuJSFb0dEZuTnfkZlJv535k6zycnk5EzGW9bc9kBnqOjYj0jaZhCPV2NnZvF9tJl3ydHbcIz24XyQtK+qHepu/IVo07X5NnIU3wKFJI1KbfR9lsW1c7t4j0h2bZTIFM17hjJ2+8++7gdTo7gx2MEr6IDJSSfh9lsrki0mnc1BSMgIl0Gt9xh2bxFZHBUfNODkrWlGQWPwyyr9eJyCY1TYlknpp38kyyzuGu++dMjB4ajGRHLDpCEckeJf0c1J/hmLk8Xj9bw1xFJDkl/Qzpz1X1uus0TjZNdi6P18/XE8tECpmSfgb0t5mju07jr3wl/8br5+uJZSKFTEk/AwbSzLFgQTA8MzJM84EH8u9kp3w9sUykkGn0TgYMdPKxQqDROyKZN6jRO2Y2xcyWmdkGM1tvZl8Pl882s+VmtsbMVprZ3HC5mdmPzKzRzN42szkxz7XQzN4NbwtT9QZzXTE3c3Q9YlHCF8muvjTvnAa+5e4zgXnAV81sJvDPwHfdfTbwN+H/AFcD54S324AfA5jZWOAeguvbzQXuMbMxKXwvOSuXmzn608GciseJSHb1mvTdfa+7rwrvtwAbgUkElwQdFa5WBewJ788HHg7n8l8OjDazCcCVwFJ3P+juh4ClwFUpfTc5KlcnHxvoOHqNvxfJX/1q0zezOuBV4DyCxP88YAQ7j4vdvcnMfgv8o7u/Hj7mJeB/A5cB5e7+f8Llfw2ccPfvdXmN2wiOEKitrb2wqbtTUyUlBjqJXL5e2EWkWKTkjFwzGwE8CXzD3Y8CtwPfdPcpwDeBh1IRrLs/6O717l5fU1OTiqeUJAY6jl7j70XyV5+SvpmVEiT8xe7+VLh4IRC5/wRBOz3AbmBKzMMnh8uSLZcsGWgHczF3TIvku76M3jGCWvxGd/9+TNEe4NPh/cuBd8P7zwB/Fo7imQcccfe9BE1BV5jZmLAD94pwmWTJQDuYc7ljWkR61pdr5F4C3AKsM7M14bK7gP8B/NDMhgJthO3wwBLgGqARaAW+BODuB83s74C3wvXudfeDKXkXMiCRjuT+jqMf6ONEJPt0cpaISIHR1MoiIgIo6QM60UhEikfRJ/1snmiknU1y2jYi6VH0bfrZOtEosrOJnX0zHy5/mAnaNiKD01ObftEn/WzNgKmzWpPTthEZHHXk9mCgJxoNtvlBZ7Ump20jkj5Fn/QHcqJRKvoBdFZrcto2IulT9El/IDNgpuKC3zqrNTltG5H0Kfo2/YFIVT+AriqVnLaNyMCpIzfF1NEoIrlMHbkppuYHEclXSvoDkKtXwhIR6U1fZtmUbixYoCQvIvlHNX0RkSKipC85SXPviKSHmnck53Sdeydy8huoSU1ksFTTl5yTipPfRKR7SvrSq0w3tWjuHZH0UdLvQm3J8bJxvYF0zr2jz1eKnZJ+jGxeUCVXZaOpJV0nv3X3+d5yS3CuhXYAUiw0DUMMTa+QKFvXG0jH3DvJPt8IXahFCoXm3umjbCW4XFZIO8Jkn2+sfHxfIl1p7p0+0jzuiQppnqG+fI7qLJZCp6Qfo5ASXKoU0jxD3X2+XRXzDl6Kg07OihFJZJrHPV6hzDMU+/k2NQU7sdjmnmLfwUtxUJu+FC1dqEUKVU9t+qrpS9EqlCMYkf5Qm36R0clJIsVNST/Fcjmp6uQzEVHST6FcT6qayExElPRTKNeTqiYyExEl/RTK9aSqk88GJpeb7ET6S0k/hXI9qerks/7L9SY7kf5S0k+hXE+qhXR2babkepOdSH/1mvTNbIqZLTOzDWa23sy+HlP2P81sU7j8n2OWf8fMGs1ss5ldGbP8qnBZo5n9VerfTnblQ1JdsCCYUKyzM/ibS7HlolxvshPpr76cnHUa+Ja7rzKzkUCDmS0FxgPzgQvc/aSZnQFgZjOBG4GPAROBF83s3PC57gc+B+wC3jKzZ9x9Q2rfUnbphJ/CUlvb/SyjudJkJ9Jfvdb03X2vu68K77cAG4FJwO3AP7r7ybBsf/iQ+cCj7n7S3bcBjcDc8Nbo7u+5ezvwaLiuSFoNpiM215vsRPqrX236ZlYHfBxYAZwLfMrMVpjZK2b2B+Fqk4CdMQ/bFS5LtlwkbQZ7tax8aLIT6Y8+z71jZiOAJ4FvuPtRMxsKjAXmAX8APG5m0wcbkJndBtwGUKtjaBmk7jpiI3MMRkbiQM9JXE12Ukj6VNM3s1KChL/Y3Z8KF+8CnvLAm0AnMA7YDUyJefjkcFmy5XHc/UF3r3f3+pqamv6+H5E4vXW4aiSOFJu+jN4x4CFgo7t/P6boaeAPw3XOBcqAD4BngBvNbJiZTQPOAd4E3gLOMbNpZlZG0Nn7TCrfjEhXulqWSLy+1PQvAW4BLjezNeHtGmARMN3M3iHolF0Y1vrXA48DG4DngK+6e4e7nwb+AnieoDP48XBdkbTR1bJE4vXapu/urwOWpPjmJI+5D0gY3+DuS4Al/QlQZDB0tSyReDojV/JWX4diRk5Ic4dHHtFIHCluSvqSlwY6J85AzkjWhGtSSJT0JS9lak4cTbgmhUZJX/JSpubE0YRrUmiU9CUvZWoaa024JoVGSV/yUqbmxMn1aySI9JeSvuSlTM2JownXpND0ee4dkVyTiTlxYsf579gR1PDvu0/DPCV/KemL9EITrkkhUfOOCBqLL8VDNX0pepGx+JGhmX2dclkkH6mmL0VPY/GlmCjpS9HTWHwpJkr6KaR24fyksfhSTJT0U0RztOQvjcWXYqKknyJqF85fuvi5FBPz2CtK5Jj6+npfuXJltsPokyFD4i/OEWEWTOMrIpIpZtbg7vXdlammnyJqFxaRfKCknyJqFxaRfKCknyJqFxaRfKAzclNIc7SISK5TTV9EpIgo6YuIFBElfRGRIqKkLyJSRJT0RUSKiJK+iEgRUdIXESkiSvoiIkVESV9EpIgo6YuIFBElfRGRIqKkLyJSRJT0RUSKiJK+iEgR6TXpm9kUM1tmZhvMbL2Zfb1L+bfMzM1sXPi/mdmPzKzRzN42szkx6y40s3fD28LUvx0REelJX+bTPw18y91XmdlIoMHMlrr7BjObAlwB7IhZ/2rgnPB2EfBj4CIzGwvcA9QDHj7PM+5+KIXvR0REetBrTd/d97r7qvB+C7ARmBQW/wD4NkESj5gPPOyB5cBoM5sAXAksdfeDYaJfClyVurciIiK96VebvpnVAR8HVpjZfGC3u6/tstokYGfM/7vCZcmWd32N28xspZmtbG5u7k94IiLSiz4nfTMbATwJfIOgyecu4G9SHZC7P+ju9e5eX1NTk+qnFxEpan1K+mZWSpDwF7v7U8BZwDRgrZltByYDq8zsTGA3MCXm4ZPDZcmWi4hIhvRl9I4BDwEb3f37AO6+zt3PcPc6d68jaKqZ4+77gGeAPwtH8cwDjrj7XuB54AozG2NmYwg6gJ9Pz9sSEZHu9GX0ziXALcA6M1sTLrvL3ZckWX8JcA3QCLQCXwJw94Nm9nfAW+F697r7wQFHLiIi/dZr0nf31wHrZZ26mPsOfDXJeouARf0LUUREUkVn5IqIFBElfRGRIqKkLyJSROTAk4wAAAU0SURBVJT0RUSKiJK+iEgRUdIXESkiSvoiIkVESV9EpIgo6YuIFBElfRGRIqKkLyJSRJT0RUSKiJK+iEgRUdIXESkiSvoiIkVESV9EpIhYcM2T3GRmzUBTip92HPBBip8zHRRnainO1MqHOPMhRkhPnFPdvaa7gpxO+ulgZivdvT7bcfRGcaaW4kytfIgzH2KEzMep5h0RkSKipC8iUkSKMek/mO0A+khxppbiTK18iDMfYoQMx1l0bfoiIsWsGGv6IiJFS0lfRKSIFHTSN7MpZrbMzDaY2Xoz+3q4/F/MbJOZvW1mvzaz0bkYZ0z5t8zMzWxcLsZoZv8z3J7rzeyfsxVjT3Ga2WwzW25ma8xspZnNzXKc5Wb2ppmtDeP8brh8mpmtMLNGM3vMzMpyNM7FZrbZzN4xs0VmVpqLccaU/8jMjmUrvpg4km1PM7P7zGyLmW00s6+lLQh3L9gbMAGYE94fCWwBZgJXAEPD5f8E/FMuxhn+PwV4nuAktXG5FiPwh8CLwLCw7Ixc3JbAC8DV4fJrgJezHKcBI8L7pcAKYB7wOHBjuPwnwO05Guc1YZkBv8zVOMP/64FHgGPZjLGX7fkl4GFgSFiWtt9RQdf03X2vu68K77cAG4FJ7v6Cu58OV1sOTM5WjJA8zrD4B8C3gaz2uPcQ4+3AP7r7ybBsf/ai7DFOB0aFq1UBe7ITYcADkZpnaXhz4HLgV+HynwPXZyG8qGRxuvuSsMyBN8n+b6jbOM2sBPgXgt9Q1vXwud8O3OvuneF6afsdFXTSj2VmdcDHCfassb4M/Gem40kmNk4zmw/sdve1WQ2qiy7b8lzgU2GTxCtm9gfZjC1Wlzi/AfyLme0Evgd8J3uRBcysxMzWAPuBpcBW4HBMhWQXH+78s6ZrnO6+IqasFLgFeC5b8cXE0l2cfwE84+57sxvdh5LEeRbwhbDp8T/N7Jx0vX5RJH0zGwE8CXzD3Y/GLL8bOA0szlZssWLjJIjrLuBvshpUF91sy6HAWIJD1L8EHjczy2KIQLdx3g58092nAN8EHspmfADu3uHuswlqyXOBj2Y5pG51jdPMzospfgB41d1fy050H+omzkuBG4D/m93I4iXZnsOANg+mY/h/wKJ0vX7BJ/2wJvIksNjdn4pZ/kXgOmBBeIiaVd3EeRYwDVhrZtsJviCrzOzMHIoRgtroU+Fh65tAJ8EEUlmTJM6FQOT+EwRJNie4+2FgGfAJYLSZDQ2LJgO7sxZYFzFxXgVgZvcANcCd2Yyrq5g4/xA4G2gMf0OVZtaYzdhiddmeu/jw+/lrYFa6Xregk35Y43wI2Oju349ZfhVBG99/c/fWbMUXE09CnO6+zt3PcPc6d68j+FLMcfd9uRJj6GmCHxdmdi5QRhZnNuwhzj3Ap8P7lwPvZjq2WGZWExk1ZmYVwOcI+h+WAf89XG0h8JvsRBhIEucmM7sVuBK4KdIOnU1J4mxw9zNjfkOt7n52Dsa5iZjfEcH3dEvaYsiBSm7amNkngdeAdQQ1UAiaTH5EcDh1IFy23N2/kvkIA8nidPclMetsB+rdPSsJtYdt+SLBoehsoB34X+7+u2zECD3GeRT4IUFzVBtwh7s3ZCVIwMxmEXTUlhBUvh5393vNbDrwKEGT2Wrg5kgneY7FeZpgRFlLuOpT7n5vlsJMGmeXdY65+4hsxBcTQ7LtOZqgmbkWOAZ8JV19eQWd9EVEJF5BN++IiEg8JX0RkSKipC8iUkSU9EVEioiSvohIEVHSFxEpIkr6IiJF5P8DmfB6lFfS17oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Stochastic Gradient Descent\n",
    "Also, the library implements a Stochastic Gradient Descent model:"
   ],
   "metadata": {
    "id": "dGMfUsn9e6Zz"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "from sklearn.linear_model import SGDClassifier\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "\n",
    "X = np.append(failure_features, normal_features, axis=0)\n",
    "#y is a vector of ones\n",
    "y = np.ones((2*no_samples,1))\n",
    "# make the normal features negative:\n",
    "y[no_samples:]=-1\n",
    "\n",
    "clf_sgd = make_pipeline(StandardScaler(), SGDClassifier())\n",
    "\n",
    "clf_sgd.fit(X, y)\n",
    "\n",
    "xx = np.linspace(21, 37, no_samples)\n",
    "yy = np.linspace(2300, 3300, no_samples)\n",
    "YY, XX = np.meshgrid(yy,xx)\n",
    "xy = np.vstack([XX.ravel(), YY.ravel()]).T\n",
    "Z = clf_sgd.decision_function(xy).reshape(XX.shape)\n",
    "#Prepare the figure \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "#Plot the data in failure state in red\n",
    "ax.scatter(temp_readings_failure, rpms_failure, color='red')\n",
    "#Plot the data in normal operation in blue\n",
    "ax.scatter(temp_readings_normal, rpms_normal, color='blue')\n",
    "\n",
    "ax.contour( XX, YY, Z, colors=\"k\", levels=[-1, 0, 1], alpha=0.5, linestyles=[\"--\", \"-\", \"--\"])\n",
    "\n",
    "plt.show()"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 300
    },
    "id": "yanQIENde63K",
    "outputId": "c29abee4-9b58-408b-850f-0fc957b4ea3d"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "/usr/local/lib/python3.7/dist-packages/sklearn/utils/validation.py:993: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhcV53n//cpqbSv1mpL1mZLlmzZkuo6cdZusjnB6QdCQwaSdCb0giFAAwOZ7qb5zdDDTJ7pBciEngAJTXqAGEI20ul0EjsOBJI4seNb2hdLli3JshZL1r5YS9X5/VGli+RIsmRLqirV9/U8ely6Vao6p2R97qnvOfdepbVGCCFEcLD5ugFCCCHWjoS+EEIEEQl9IYQIIhL6QggRRCT0hRAiiIT6ugGLSU5O1jk5Ob5uhhBiRnU1TE5+cHtYGOzcufbtEfMyTbNXa50y331+Hfo5OTkcP37c180QQsywLVAcmJoC+Vv1G0qp1oXuk/KOEGLpsrKWt134HQl9IdaTAwcgJ8czIs/J8Xy/kh5+GKKi5m6LivJsFwFBQl+I9eLAAdi/H1pbQWvPv/v3r2zw33cfPPEEZGeDUp5/n3jCs10EBOXPp2HYvXu3lpq+EEuUk+MJ+otlZ0NLy1q3RviQUsrUWu+e7z4Z6QuxXrS1LW+7CEoS+kKsFzLJKpZAQl+I9cIfJ1lXe2JZLJuEvhDrhb9Nsq7FxLJYNpnIFUKsDplY9hmZyBVCrD2ZWPZLEvpCiNUhE8t+SUJfCLE6/HFiWUjoCyFWib9NLAvAz8+yKYQIcPfdJyHvZ2SkL4QQQURCXwghgohfh/7U1BRnz57Fn48lEEKIQHLJ0FdKRSiljimlKpVStUqp/+HdfkApdUIpVaOUelIpZfduV0qp7ymlTiqlqpRSjlnP9YBSqsn79cClXntkZIQf/ehH/Ou//qsEvxBCrICljPQngJu11iVAKXCHUuoa4ABQCOwEIoG/8D7+w0C+92s/8AMApdQG4JvAHuBq4JtKqcTFXjghIYE/+qM/oqCgAKUUWmt+/etfc+bMGdkJCCHEZbjk6h3tSdcR77d275fWWr8y8xil1DEg0/vtR4Gfen/uPaVUglJqI/Ah4HWtdZ/3Z14H7gB+sdBrT09Ps2nTJjZu3AjA4OAg7733Hr/73e9ITU3FMAx27dpFZGTk8nothBBBakk1faVUiFKqAjiHJ7iPzrrPDtwPvObdlAGcmfXj7d5tC22/+LX2K6WOK6WOd3R08MQTT/D444/z/vvvExERwUMPPcRHPvIR7HY7r776Kt/5zndok8O6hRBiSZa0Tl9r7QJKlVIJwK+UUsVa6xrv3d8Hfqe1fmslGqS1fgJ4AsAwDH3nnXdimib/8R//waFDhyguLsYwDP7iL/6C7u5uKioqrE8ClZWVjI2NUVJSQtTFRwIKIYRY3sFZWusBpdRv8JRlapRS3wRSgM/OethZYPOs7zO9287iKfHM3v7mYq83PT1NRkYGu3fvprOzE9M0qa6upry83Crv/OEf/iF2ux2A5uZmqqqqOHz4MEVFRRiGQU5ODkqp5XRTCCHWrUueWlkplQJMeQM/EjgE/AOQDvwZcIvWenzW4+8EvgjswzNp+z2t9dXeiVwTmFnN4wSMmRr/fLZs2aLvv/9+srOz+fSnP41SiomJCWpqajBNk46ODkJDQ9mxYweGYbB582Z6enowTZPKykouXLhAaWkpd91112W/QUIIEWgWO7XyUkb6G4GfKKVC8MwBPKO1flkpNQ20Au96R9IvaK2/BbyCJ/BPAmPAnwJorfuUUv8TeN/7vN9aLPDBs3pn3759TExMWKt3jhw5QkFBAZ/5zGfo6uqyRv+VlZWkpKRYo/9bb72Vuro64uPjARgaGuLgwYM4HA7y8vJk9C+ECEoBdRGVgYEBHnvsMaampkhLS7NW79hsNmpqanA6nbS3txMaGmqVd7Kzs1FK0dzczHPPPcf4+DiJiYk4HA7KysqIiYnxYQ+FEGLlLTbSD6jQB5iYmKC6uhrTNOns7MRut3P//feT5T1Hd3d3N6ZpUlVVxYULF0hKSsIwDEpKSggPD6e+vh6n08np06cJDQ3loYceIiIiwhfdE8HqwAH4xjc8FxPJyvKcalhOSiZW0LoK/dk6OjqorKzk1ltvxW63U1VVZa3eCQ0Npa6uDtM0aWtrIyQkhMLCQgzDIDc3l76+Ptra2igrKwPgxRdfJDExkbKyMuLi4taqiyLYzFw3dmzs99uiouSUw2JFrdvQv9jzzz9PdXU1oaGhbN++HcMwyMrKore315rcHR8fZ8OGDTgcDkpLS4mJicHlcvHzn/+c5uZmlFIUFBRgGAZbt27FZvPr0xOJQCPXjRVrIGhCH6Crqwun00llZSUTExOUlZXx0Y9+FPAsAa2vr8c0TVpaWrDZbNboPy8vj4GBAZxOJ+Xl5YyMjLBv3z6uvvrq1eiaCFY2G8z3N6cUuN1r3x6xLgVV6M+YmpqitraW+Ph4cnNzGRoa4tChQ9ba/fPnz+N0OqmoqGBsbIyEhARrcjcqKorGxkays7OJioqiqqqK6upqDMOgoKBARv/i8slIP3AF0FzMlS7ZDEh2u53S0lLr+3PnztHc3ExNTQ0bNmzAMAyuv/56br75ZhoaGnA6nfz617/mzTfftMo7MxO8LpeLrq4unn76aWJjYykrK6OsrIzExEXPFyeCxXLC4OGH56/py3Vj/dvFczGtrZ7vwW+DfyHrdqQ/n6mpKau809ra+oHVO319fVZ5Z3R0lPj4eCvgY2NjaWpqwjRNmpqaSE9P57Of9RyIrLWWdf/B6nImZq9kxBhAo811JcA+oQVleedSent7aWtrw+HwHCD80ksvkZiYSGlpKVFRUZw4cQLTNK3J3fz8fAzDID8/n+HhYUZGRsjIyGBiYoLHH3+c7du343A42LBhw6q0V/iptQwDWfnjOwE2FyOhfwkul4sDBw5w6tQpbDYb27Ztw+FwsGXLFgYHB+dM7sbGxlq1/4SEBAYGBnjttddobGzE7XaTl5eHw+GgsLCQ0NB1Wz0TM9YyDAJstLmuBNh7L6G/RLMnd0dHR+es3nG5XFZ55+TJkwBs2bLFmtwdHR2loqICp9PJwMAADz74IGlpabhcLkJCQtasD2KNrWUYBNhoc10JsE9ZEvrL5HK5OHHiBNnZ2URHR1NdXU1NTQ0Oh8Mq78yM/oeGhoiJiaGsrAyHw0F8fDzt7e3WEcIvvPACQ0NDGIZBUVGRjP7Xm7UMgwAbba47ATSfIqF/hWZW9oyMjBAXF2cFfGxsLCdPnsQ0TRobG9Fak5eXh2EYFBYWEhISwrvvvsuxY8fo7+8nMjKSkpISDMMgJSXF190SK2WtwiDARpvCdyT0V4DL5aKxsdGa3L149c7w8PCc8k50dDSlpaXW5O7p06cxTZOGhgauueYabrvtNtxuNy6Xy7oegBCXFECjTeE7EvorbGBggJGRETIzM5mYmOCJJ55gx44dlJWVER8fz6lTpzBNkxMnTuB2u8nNzcXhcFBUVMTExAQA0dHRnDx5kueee45du3ZhGAZpaWk+7pkQYj2Q0F9FAwMDvPLKKzQ1NQFY5Z1t27YxNjZmjf77+/uJioqipKQEh8NBSkoK3d3dvP3229TX1zM9PU1mZqZ1umiZ/F0HAm1UHmjtFQuS0F8DQ0NDlJeX43Q6GRwctFbvuN1ulFJzyjsul4vs7GwcDgfbt29namqKqqoqTNNkYmKCr3zlK9hsNkZHR4mOjvZ118TlCLT6e6C1VyxKQn8Nud1uzpw5Q3Z2NgC/+tWvrNU7hYWFTExMUFFRgWma9PX1ERERMWdyd3h4mLi4OFwuF4888ghxcXEYhkFxcTHh4eE+7p1YskBbaRNo7fU1P/9UJKHvQ0eOHOHYsWMMDAwQFRVlTe4mJSXR0tKCaZrU19fjcrnYvHkzhmGwY8cOlFI4nU5M06S7u5uwsDB27tzJtddeS3Jysq+7JS4l0NbUB1p7fSkAPhVJ6PuY1prm5macTicNDQ1ce+213HbbbWitcblcTE5OUllZiWma9Pb2Eh4ePmdyt729HafTSU1NDffccw95eXmMjo4SEhIiV/3yV4E2cg609vpSALxXEvp+ZGRkBICYmBhOnjzJ888/P2dyt62tDdM0qaurY3p6moyMDKu843a7CQ8PRynFa6+9hmmaFBcXYxgGGRkZctI3fxIAo8E5Aq29vhQAn4ok9P1UV1cXb731ljW5O1Pe2blz55zRf09Pj1XeMQyDTZs20dnZyfHjx6murmZycpK0tDT27NljnUBOrKGF6rt+Xvf9gEBrr6/ISH/1rPfQnzE6OmoF/NTUlLV6Z2xsjMjISNrb2zFNk9raWqampti4caO1cwCoqanBNE1SU1O56667AM/1gzdu3Cij/9UmI+TgEwC/cwn9AKG1ZmhoiPj4eNxuN4888gjx8fHW5K7b7baWdnZ3d2O329m5cycOh4OMjAxcLhehoaF0dnby+OOPk5KSYq37j4qK8nX31qcAGPWJWVbq04yffyqS0A9A09PTHD9+3CrvhIeHW6t3NmzYQEdHB6ZpUl1dzdTUFGlpaVbA22w2amtrMU2T9vZ2QkNDKSoqYu/evcTGxvq6a+tLANR3hVcAjNBXioR+ANNac+bMGWv1zr333kteXh5jY2PWUbvV1dWYpklnZyd2u50dO3ZgGAaZmZmcO3fOOiXEF77wBcLCwujs7CQuLk4O/FoJMtIPHEH0u5LQXycuXLhgrd45ePDgnNU7M5O7M6P/yclJq7xTUlJCeHg4NpsNrTU/+MEPOH/+PIWFhRiGQW5urtT+L1cQjR4DXhB9KpPQX4c6Ojp4//33qampYWpqivT0dPbs2UNZWRmTk5PW5O7Zs2cJDQ1l+/btGIZBVlYWPT09OJ1OKisrGR8fZ8OGDdx8880UFxf7uluBZaau29oKISHgcnlGjX5W3xVeMtL33CehH9guXLhglXfS0tL42Mc+BniWg6alpdHd3Y1pmlRVVTExMUFycjIOh4PS0lLCwsKoq6vD6XSye/duiouLGR0dpbOzky1btsjofzEywg88QfQ7k9APAlprpqensdvt1uqd1NTUOWftrKurwzRNzpw5Q0hICEVFRRiGQU5ODgBKKY4cOcKhQ4dISEiwrgUsk7/zCKJRo99bzkoaP191s1Ik9IPMTHnn+PHjdHR0EBoayo4dO7j11luJjY21Jnerqqqs8o5hGJSWlhIeHk5DQwNOp3POheLvvvtubDabr7vmP4KoPuzXgmj0vhwS+kGss7MTp9PJiRMn+OIXv0hYWBhdXV3ExcVht9upr6/HNE1aW1ux2WzW5G5eXh79/f04nU5GRkasg76qqqrIzs4mPj7exz3zMRnp+wf5PcxLQl/gdrs/sHpnZnI3Ozub8+fPY5omlZWVjI2NkZiYSFlZ2ZzyzujoKN/5znfQWpOfn49hGOTn5wfnJwAZYfoHf/vE5SflIwl9MUd3d7e1eufChQskJSVx8803s2PHDqanp2loaMA0TU6fPo3NZqOgoADDMNiyZQuDg4OUl5dTXl7O8PAwsbGxfOITn7CuHxBU/OQPfF1Z7nvqTyN9PxoISOiLeU1NTVmTu1dddRU7d+5kdHSU7u5ucnNz6evrw+l0UlFRwejoKPHx8dbkbnR0NE1NTTidTj7ykY8QExNDS0sL4+PjFBQUyOUexfwWC/XLCc35fkYpz+h/rZfP+tEO6IpCXykVAfwOCAdCgee01t9USuUCTwNJgAncr7WeVEqFAz8FDOA88EmtdYv3ub4O/DngAr6ktT642GtL6K8drfWc1TuJiYlWwEdGRnLixAlM06S5uRmlFAUFBTgcjjnlnWeffZba2lpiYmIoKyvD4XCQmJjo454Jv3GpUL/c0Jx9vMRM4M/3/KttoVITePqwhp8IrzT0FRCttR5RStmBt4EvA18FXtBaP62U+iFQqbX+gVLq88AurfXnlFKfAj6mtf6kUmo78AvgamATcBgo0Fq7FnptCf21Nz09bU3utrS0WJO7n/jEJ7DZbNbkbnl5OSMjI8TFxVkBHxsby8mTJzFNk8bGRrTWXHXVVdx5552+7pbwB5cK9Sutz/t6pL3Q6/tgR7Ri5R2lVBSe0H8Q+A8gXWs9rZS6Fvg7rfXtSqmD3tvvKqVCgS4gBfgbAK31//Y+l/W4hV5PQt+3ent7rdU7f/zHfwx4zvOTnZ1NdHQ0jY2N1ugfYOvWrdbk7ujoKOXl5SQkJFBSUsLk5CRvvfUWpaWlJCUl+bJbwlcuFepXGtq+ntRdrNR0sVXeES0W+qFLfIIQPCWcrcBjQDMwoLWe9j6kHcjw3s4AzgB4dwiDeEpAGcB7s5529s8IP5ScnMzevXut70dHR/nVr36F1tqa3L333nsZGhqyJneffvrpecs77e3tvPPOO7z11lvk5ubicDgoKioiNHRJ/wXFepCVNX+oZ2V5/n344fnLPw8/vDLPv9pmRu6z5yzmaw947veR5Y70E4BfAf8N+H9a663e7ZuBV7XWxUqpGuAOrXW7975mYA/wd8B7WuunvNt/7P2Z5y56jf3AfoCsrCyjdaE3TfjE7MndmfLOxz/+cbKzs3G73TQ1NWGaJk1NTWit2bJlC4ZhsG3bNsbGxqioqMDpdNLf309UVBSf//zniYmJ8XW3xFpYykTtlayI8qPVMxYflZwWG+mjtV7WF/Dfgf8K9AKh3m3XAge9tw8C13pvh3ofp4CvA1+f9TzW4xb6MgxDC/80PT2t6+rq9FNPPaWHhoa01lq3tLTo+vp67XK59ODgoH7zzTf1d7/7Xf3Nb35T/+M//qM+dOiQ7u3t1W63Wzc3N+tDhw5Zz/f222/ryspKPTk56asuibXw1FNaZ2drrZTn36eeWr3nT0ryfC31tVajbU89pXVUlNaeIo/nKypq5ft9EeC4XiBXlzKRmwJMaa0HlFKRwCHgH4AHgOf17ydyq7TW31dKfQHYqX8/kfvHWuv/pJTaAfyc30/kvgHka5nIXTdmVu/ExsZa5Z24uDiam5utyV23201ubi6GYVBYWEhoaChaax5//HG6urqIiIigpKQEwzBITU31dZdEoFruqH81PyVc/Oll3z545ZVVXc1zpat3dgE/AUIAG/CM1vpbSqk8PEs2NwDlwJ9orSe8Szx/BpQBfcCntNanvM/1DeDPgGngK1rrVxd7bQn9wOJ2u2lsbMTpdNLU1ATA7t27rdU7w8PDVFRUYJomAwMDREVFUVpaisPhICkpiZaWFkzTpL6+HpfLxe233861117ryy6JQLXcsspalWHWqAQlB2eJNTdz5G5CQgKlpaVMTk7y9ttvU1paSmJiIqdOncI0TRoaGnC73WRnZ2MYBtu3b2dycpLKykry8/NJTk6mtbWVmpoaDMMgPT3d110TgWC5K3nWauXPGu1cJPSFzzU3N/PUU0+htSYvL88q74yPj1NZWYlpmvT19REZGUlJSQkOh8Mq77z//vscPHiQ6elpMjIyMAyD4uJiwsLCfNwr4bf8daS/RjsXCX3hF4aGhqzVOwMDA0RHR/Pggw8SExOD1voD5Z3NmzdjGIZ1TqCqqipM0+TcuXMkJSXxxS9+US70IubnTzX92WSkvzgJ/fXJ7XZz6tQpTp06ZR0HcOTIEWJjYykqKmJiYoLKykqcTie9vb1ERESwa9cua3K3vb2d4eFhtm/fjtvt5umnn6agoICdO3cSHh7u494Jv7Hc5Z9rcQI9qekvTkI/OGit+eEPf0h3dzeRkZHW5G5ycjJtbW2YpkldXd285Z3BwUF+/vOf093djd1uty4Un5GRIZ8CxMJ8eYbUNXhtCX3h97TWnD592prcdblc3HHHHVxzzTUAc2r/PT09hIeHs3PnTmtyt6OjA9M0qa6uZmpqik9/+tPWZSCD2no7/fNK9OdSo+118J5J6IuAMjo6SkVFBdu2bbNG+7W1tdbk7pkzZzBNk9raWqanp9m4cSOGYbBz504A6urqKCkpwWaz8eabbzIwMIDD4WDz5s3BNfr3xyNUr8RK9WexuvpCp4IIsPdMQl8EtGPHjnHw4EFcLheZmZnW5K7L5aK6uhrTNOnu7iYsLMwq72zatAmlFG+88QZHjx5lcnKSlJQUDMOgpKSEyMhIX3dr9fn6rJMr7VJhvdTR+WIraBY6X06AvWcS+iLgjY2NWeWd3t7eOat3tNacPXsW0zSpqalhamqK9PR0a/Rvs9moqanBNE3Onj3Lrl27rLOGau91BNYlX591cqUtdr76qKilj84X2nksJsDeMwl9sW5orTlz5gzDw8Ps2LEDt9vNL3/5SwoKCiguLkZrTU1NDcePH6erqwu73c6OHTswDIPMzEy6u7sJCQkhJSWFc+fO8cwzz1ij/6ioKF93b2UFy0g/JARc85zNZaF+zlcmupQAe8+u+NTKQvgLpRRZs06VOzw8TH9/P//+7//OwYMHrcnd/fv309XVZU3uVlRUkJqaimEY7Nq1C4DJyUkiIyM5ePAghw8fpqioCMMwyMnJWR+j/ys9VbG/2bcPfvjDD16QZKHwXuj0xbNPgbyUEX8gv2fzkJG+CHhaa9rb263J3YtX70xMTFjlnY6ODkJDQ9mxYwcOh4OsrCx6enqs00W7XC4eeughwsPDcbvd1qUgA9Y6WIkCLHyBks99znPyssv9RLNYyWimxh+A75mUd0TQuHDhAnV1dZSWlmKz2fjtb3/LwMCAtXa/q6sLp9NJVVUVExMTJCcnW+Udu91OV1cXmzdvRmvNv/zLvxAfH49hGOTl5a2P0X+gWq0VN+utBOYloS+C1uzVO2lpaVZ5x2azUVtbi2matLe3ExISwvbt2zEMw7ogzOHDh6msrGRsbGzOheLnXPRlvYyk/d2lJqUv9/ew3pa1eknoi6B2cXnn4tU7586dw+l0UllZyYULF0hKSsLhcFBaWkp4eDgNDQ2Ypsnp06e56667KC0tZXp6GtsvfoHtc59bd4Hhl1ZzRL4Od9wS+mLdWu7fa2dnJ6GhodbqnWeffdYa/dvtdurq6jBNk7a2NkJCQigsLMQwDHJzc+nr6yMuLg673c67777Le3feiaO/nzIgbvaLBHhpwC+t0xH5apHQFwFjOSF+pTnQ3t7Oa6+9Rnt7O6GhoWzfvh2Hw0F2dja9vb2YpkllZSXj4+MkJiZiGAalpaXExMTQ3NzMka1bacZzLdB8YDdQAAG3pjtgrMMR+WqR0BcBYbkhvlKf+Lu7uzFNk6qqKlwuF1/72teIiIjA7Xbjdrupr6/HNE1aWlqw2Wxs27YNwzDYcsstDLS14cRz6bhU4D97GzBWV7f+1v2LgCGhLwLCckN8pQ84nZqaorOzk6ysLLTW/PjHPyYhIcFau3/+/HlraefY2BgJzc04nn2WsokJooAxIDYqiuFHH+X/dHZaF4vJz88nJCRk+Q0S4jItFvoBvgh5ZR044Akem83z74EDvm5RYFvu+7nQsTQLbZ91jNaStl+K3W63Dvxyu91kZmbS3NzMT37yE/75n/+ZhoYGrrvuOr761a9y9913s+G22/j1vn08Eh/PM0Dnpk24f/hDQu69lxtuuIHu7m6efvppHnnkEd544w1GR0cvr2FCrCAZ6Xv58zxRIJYyL+f9XO5Ify1+Z1NTU1Z5p7W1dc7qnZCQEPr7+3E6nZSXlzM6OkpcXJy1tDM2NpampiZM06S5uZkvf/nLxMXFMTQ0RHR0tIz+xaqR8s4S+OsxGv68M1rM5byfl9PXtdwh9vb2Eh8fb63eOXbsGGVlZZSVlREVFcWJEycwTZNTp04BkJ+fj8PhoKCggImJCevMnj/72c/o6uqyLhaTlJS0Og0WQUtCfwmWUh/2xYjbX3dGl3K59fZA+VRz8uRJ3nnnHU6fPo3NZqOgoMCq3w8MDFij/+HhYWJjYykrK8PhcJCQkMDJkyc5fvw4jY2NuN1ucnNzue6668jPz/d1t8Q6IaG/BJcKV1+NuAP17LiBurNartmTu6mpqTzwwAOA50pf4eHhNDY24nQ6aWpqAmDLli0YhkFBQQFjY2OUl5fjdDq56qqruP7665menmZgYIDk5GRfdksEOAn9JbhUqPsqxAI1PAO1LHW5XC6XVdMfHh7m0UcfnbN6Z3h42Ar4oaEhYmJirPJOYmIiLpeL0NBQqqqqeOGFF8jOzsYwDLZv305oqJwMVyyPhP4SLVZa8NWIO5DDM1BKNSttdHSU9957j/LyckZGRoiNjcXhcHD11VcTGRnJyZMnMU2TpqYm3G63tXMoLCxkfHyciooKnE4nfX19REZGUlJSwq233irhL5ZMQn8F+HLEHazhGehcLtcHVu/Ex8czPDxMVFSUdS1gp9PJwMAA0dHRlJSUYBgGGzZsoKWlBdM06evr4zOf+QxKKc6cOUN6ejp2u93X3RN+TEJ/Baz1iFuCfn0ZGxuzjtB96qmn6OrqsiZ34+PjOXXqFKZpcuLECdxuNzk5ORiGQVFRETabDZvNxuTkJN/+9rex2Wzs2rULwzBIS0vzcc+EP5LQXyFrFcQL7WAeeMBzvQjZEQS2xsZGTNOksbERrTVbtmzh2muvZevWrYyMjFi1//7+fiIjI63af3JyMq2trZimSX19PdPT02RmZnL77bezefNmX3dL+BEJ/QCzUClJqQ9eKc7fa/vyiWVhQ0NDVsBfffXV1uqdwcFBNmzYwOnTpzFNk4aGBlwuF1lZWdbk7tTUFFVVVZimycc//nHS09Pp7e1lamqKjRs3+rprwsck9APMYldwu5g/r+IJ5EnotTRzYrfZq3dyc3Otyd2JiQmr9n/+/HkiIiIoKSnB4XCQmppqXdHrxRdfpKKigk2bNuFwONi5cyfh4eE+7p3wBQn9ALPQSH8+/rxeP1CXm/rSxeWdqKgoSkpKuOWWWwgJCbHKO3V1dbhcLjIzMzEMgx07duByuaiursY0Tbq7uwkLC+Oqq67itttu83W3xBqT0PcDV3qe+ItLOzP8OUAD9cAyf6C1tiZ3+/v72b9/P0op2tvbSUtLY2pqisrKSkzTpLe3l/Dw8DmTu2fPnsXpdBIbG8tNN92E1pqKigqKirxn6VkAABypSURBVIqIiIjwdffEKpPQ97GVOKfMvn3wk58EVqlERvorw+12f2D1zszSzpSUFNra2nA6ndTW1jI9PU1GRgYOh4Pi4mKrvNPW1saTTz6J3W5nx44dGIZBZmamXOx9nZLQ97GVCr9AmxSVmv7K0lpba/fr6+txuVxs3ryZvXv3snnzZsbHx63J3XPnzhEWFsbOnTsxDIP09HS6urowTZPq6momJydJTU3l3nvvJSEhwdddEytMQt/HgrnMEWg7qkAxOjpqlXfuvvtu0tPTOX/+PFNTU6SlpdHe3o7T6aSmpsZa0TMzuauUora2loaGBj71qU9hs9mor68nKiqKrKwsGf2vA1cU+kqpzcBPgTRAA09orR9VSpUCPwQigGng81rrY8rzP+ZRYB+eiwl9Wmvt9D7XA8D/533q/6W1/slir71eQl/KHGK1aK0/sHonIyMDwzAoLi7G7XZbk7tdXV3Y7XaKi4sxDIOMjAyUUmit+f73v09PTw/JyckYhkFJSYlc7jGAXWnobwQ2aq2dSqlYwATuAv4P8IjW+lWl1D7gr7TWH/Le/ks8ob8HeFRrvUcptQE4juf60dr7PIbWun+h114voe/PZY7LHYnLCN7/jI+PW6P/np4ewsLC2LNnD7fccgtaazo6OnA6nVZ5Jy0tDcMw2LVrFzabjdraWkzTpL29nZCQEPbu3cuePXt83S1xGRYL/UuewUlr3Ql0em8PK6XqgQw8wR3nfVg80OG9/VHgp9qzN3lPKZXg3XF8CHhda93nbdTrwB3ALy63Y4FiJgz9LSQv3hm1tnq+h8Xbdrk/J1ZXZGQk11xzDXv27KG9vR3TNOeUanp6eti7dy979+6lpqYG0zR55ZVXeP3119m+fTuGYfDnf/7nnDt3DqfTaZ3i4fz58zQ0NFBaWkp0dLSvuidWyLJq+kqpHOB3QDGe4D8IKDzX2r1Oa92qlHoZ+Hut9dven3kD+Gs8oR+htf5f3u3/DRjXWn/7otfYD+wHyMrKMlqXumBdLNvllp2kXBV4ZlbvhIWFWeWdTZs2zZncnZiYICUlxRr9z5R3jh49yquvvkpISAiFhYU4HA7y8vKk9u/HVmQiVykVA/wWeFhr/YJS6nvAb7XWzyul/hOwX2t965WG/mzrpbzjry53gjmYJ6YD1Ux5ZybgZyZ877nnHhISEpicnKSmpgan00l7ezuhoaEUFRVhGAbZ2dn09vZaF4sZHx8nNTWVz33uc9hsNl93Tczjiso73iewA88DB7TWL3g3PwB82Xv7WeBfvLfPArPP/pTp3XYWT/DP3v7mUl5frI6srPlH7FlZq/NzwneUUmRkZJCRkcHtt99OdXU1DQ0NxMV5KrSnTp0iJSXFKu+YpklVVRXV1dUkJSVhGAY33HADt9xyC/X19QwODlqBf/jwYXJyctiyZYuM/gPAUiZyFfAToE9r/ZVZ2+uBB7XWbyqlbgH+UWttKKXuBL7I7ydyv6e1vto7kWsCDu9TOPFM5PYt9Noy0l9dlzvB7M8T02L5Zq/eSU1NxeFwUFJSQmhoqDW5e+bMGUJCQigqKsLhcJCbm4tSitHRUb7//e8zOjpKQkICDoeDsrIyYmNjfd2toHalq3duAN4CqoGZD+9/CwzhWZoZClzAs2TT9O4k/i+eSdox4E+11se9z/Vn3p8FT5noXxd7bQn91SerdwRglXdM0+Ts2bOEhoZy2223Wat3Zo/+x8fH2bBhAw6Hg9LSUiIjI2loaMA0TU6dOoXNZuPee+9l69atPu5V8JKDs4QQS9bV1YXT6aSoqIjc3Fz6+vo4ceIEJSUl2O126uvrMU2T1tZWbDYbhYWFGIZBXl4e/f39lJeXc8MNNxAeHk5VVRV9fX2UlZURHx/v664FDQl9IcRlm716Z2ZyNycnh/Pnz2OaJpWVlYyNjc1b3nn11Vc5duwYAPn5+TgcDgoKCmQCeJVJ6F+ClCqEWNxMeaeyspILFy6QlpbGZz/7WWw2G9PT01Z55/Tp09hsNgoKCjAMgy1btjA0NITT6aS8vJzh4WGKior45Cc/6esurWsS+ovw5aSk7GwWJu+Nf5qamrJW79x4440AvPHGG+Tk5JCXl0dfX5+1tHN0dJT4+HjrWsAxMTE0NjYSHh5Obm4uo6OjvPjii5SVlbFt2zZCQkJ83Lv1Q0J/Eb460EhWwCxM3pvAMTo6ymOPPcbY2BiJiYlWeScyMpITJ05gmibNzc0opcjPz8cwDPLz87HZbLS1tfHcc88xNDRETEyMdS3gDRs2+LpbAU9CfxG+OtBIjmpdmLw3gWV6epr6+nqcTqdV3pm9eqe/v98q74yMjBAbG2uN/uPi4jh58iSmadLU1ITWmq997WvExMT4uFeBTUJ/EZcbMFdafpCjWhcm703gOn/+POXl5dx4442Eh4dTXV1trd6Jjo6mqakJ0zQ5efIkAFu2bMEwDAoKChgdHaWlpYVdu3YB8G//9m9ERkZiGAZJSUm+7FbAueIjctezhx+ev5Tw8MML/8xKnHBMjmpdmLw3gSspKYlbb73V+v7MmTMcO3aMN99805rcveeeexgaGqK8vJzy8nJ++ctfEhMTY43+wXPA2OTkJJWVlRw5coScnBwMw6CoqIjQ0KCPrSsS9CN9WP6ofSXKD1K3Xpi8N+tLX1+fFfAjIyNzVu+43W6rvNPY2IjWmi1btuBwOCgsLGR8fJyKigrrWsG33norN9xww5zrCIgPkvLOClup8oOsUFmYvDfrj8vlslbv5OXlMTo6yksvvURZWRn5+fmMjo5SXl6O0+lkcHCQ6OjoOZO7p0+fJjU1lZiYGGprazl69CiGYbB9+3bsdruvu+dXJPRXmEw0CnHl2traePbZZxkeHrYmd2eO3G1ubrZG/263m9zcXAzDoLCwkNDQUOrq6jh8+DB9fX1ERERYF4pPTU31dbf8goT+CpPygxArw+12W5O7TU1NAHNW7wwPD1vlnYGBAaKioqyAT0pKorW1FdM0qaurIzExkS984QsopXC73UF91K+E/iqQ8oMQK2twcJDTp09TWloKwEsvvURUVBQOh4PExEROnTqFaZo0NDTgdrvJzs62JnenpqYYGBhg06ZNTE1N8dhjj1nHBaSnp/u4Z2tPQl8IEVC01jz77LNWwM8u71y4cMG6FnBfXx+RkZHs2rXLKu+MjIzw+uuvU1tby/T0NJs2bcIwDHbu3ElYWJivu7YmJPSFEAFpeHjYmtwdGBiwVu+AZ8fQ0tKCaZrU19fjcrnYvHkzhmGwY8cOpqenqaqqwjRNzp07x2c+8xkyMjKYnp5e98s+JfRFwJHymZhNa01zczPp6enExMRQV1c3Z/XOzJp+0zTp7e0lIiKCXbt24XA4SEtLo6uri/T0dJRSvPTSS3R2dlqj//DwcF93b8VJ6IuAIhPl4lJqa2s5fPgw/f39REZGWpO7ycnJtLW1WZO709PTZGRkYBgGxcXFhIWFUV5eztGjR+nq6sJut1NcXMxVV13Fpk2bfN2tFSOhLwKKLIkVS6G15vTp09bk7uzVO1rrObX/np4ewsPD2blzpzW5O3Oh+JqaGsrKyvjwhz+M1pqJiQkiIiJ83b0rIqEvrshal1rk3DtiuUZHRxkcHJyzemfmtA+pqam0t7djmia1tbVMTU2xceNGq7wDnpPGRUdH09LSwoEDB9i+fTuGYbB58+aAPPJXQn8ZpJY8ly9KLas50pff7/o3MjLCoUOHrPJOZmamNbnrdrutyd3u7m7CwsIoLi7GMAw2bdpEX18f7777LtXV1UxMTJCSkoJhGBiGEVBH/UroL5HUkj/IF6WW1fo9zPe8Snk+VWRnyw5gvRkbG7MCvqenZ87qnZCQEKu8U11dzdTUFGlpaRiGwa5du7DZbNTU1OB0Ounv7+erX/0qISEhDA0NERsb6/ejfwn9JZJa8gf5qtSyGiPyhX6/M4J9B79eaa3p6Ohg06ZN1uqdrq4ua3IXoLq6GtM06ezsxG63s2PHDgzDIDMzk7GxMaKjo3G73Xzve98jJCQEwzAoKSkhOjrax72bn4T+Ekkt+YPW045wod/vbIHYL7E8TqeTo0ePWuWdnTt3snv3bjZu3EhHRwdOp5OqqiomJydJTU3F4XBQUlJCWFgY1dXVOJ1O2traCAkJobCwkBtuuIGNGzf6ultzSOgv0XoKuJWynkpelxrpQ3Dv4IOJ1pqzZ8/OWb2zb98+6zz+SilqamowTZOzZ88SGhpqTe5mZWXR09OD0+mksrKSj33sY9ZFYLTWfnHVLwn9JVpPAbeS1svk53y/34sF8w4+WF24cAGXy0V0dDStra089dRTFBcX43A4yMzMpLu7G9M0qaqqYmJiguTkZKu8ExYWhs1mw2az8cYbb/DOO+9QWFiIYRjk5eX5rPYvob8M6yXgxPxmfr+trb+fxJ0hO3hx/vx5jhw5QnV1tVXeMQwDh8OB1pq6ujpM0+TMmTOEhIRQVFSEYRjk5ORw/vx5nE4nFRUVjI2NkZCQwNVXX81111235v2Q0BdiHrKDFwuZmJiwyjuDg4PW6p3h4WFiYmLo6enBNE0qKyu5cOECSUlJOBwOSktLCQ8Pp6GhAdM0iY6O5hOf+ATguX5AZmbmmpzyWUJfCCEu08jICDExMWitefTRR7Hb7dbSTrvdTn19PaZp0traak3uOhwO8vLycLvdhISE0NPTw2OPPUZ8fDwOh4OysjLi4uJWrc0S+sIio1shLo/L5bLW/be3txMaGkpRURHXX3896enp1uRuRUUF4+PjJCYmWgEfGRnJiRMnME2T5uZmlFLk5+fz4Q9/mMTExBVvq4T+GvLnUJWJaiFWRnd395zVO9u2bWNsbAytNeHh4dbov6WlBZvNxrZt26zJ3cHBQZxOJ7W1tezfv5+IiAjOnj1LdHQ0CQkJK9I+Cf014u+hKktShVhZU1NThISEWKt3jhw5Yq3eyc3Npa+vD6fTSXl5OWNjY3PKO7OP7P3Rj35ER0cHW7duxeFwUFBQQEhIyGW3S0J/jfh7qMrBZ0KsntmTuzPlnauuuorrrrsOl8tlTe6eOnUKpZR1QritW7cyNDREeXk55eXlDA0NERMTw0033YRhGJfVFgn9NeLvoervOyV/5c8lO+F/pqenrfJOdHQ0d999NwBnzpwhMzOT/v5+q/Y/MjJCXFwcZWVlOBwOYmNjaWpqwul0UlRURGlpKePj45w6dYrCwsIlj/4l9NeIv4eqv5ef/JG8Z+JKuFyuOat3EhISrPJOVFQUjY2N1uQuwNatWzEMg/z8fCvgjx8/zssvv0x0dDSlpaU4HA6SkpIWfV0J/TUSCAEho9bl8fcduQgM09PT1uqdU6dOYbPZKCgo4PbbbycxMZGBgQGr9j88PExsbCxlZWWUlZURHx9Pc3MzTqeTEydOWBeKv++++xa81u8Vhb5SajPwUyAN0MATWutHvff9JfAFwAX8h9b6r7zbvw78uXf7l7TWB73b7wAeBUKAf9Fa//1irx1ooQ8SquuNv5fsROCZmdytqanhs5/9LJGRkXR0dBAdHW2Vd0zTpKmpCYC8vDwMw7BWCFVUVHD+/HnuuusuACoqKsjMzCQ5Odl6jSsN/Y3ARq21UykVC5jAXXh2At8A7tRaTyilUrXW55RS24FfAFcDm4DDQIH36RqB24B24H3gHq113UKvHYihL9YXGemL1aK1/sDqnfz8fKu8Mzw8bE3uDg4OEh0dbdX+N2zYAHiOHP72t7/N1NQU2dnZGIZBUVERYWFhK1feUUr9G/B/gc/gGfUfvuj+r3s79L+93x8E/s57999prW+f73HzkdAXK+FKPn0FQslOBL7+/n4r4GfKOzfddBMOhwO3201zczOmadLY2GiVdwzDoLCwcM61gPv6+oiMjORv/uZvFgz9+QtCC1BK5QBlwFHgn4AblVIPAxeAh7TW7wMZwHuzfqzduw3gzEXb9yzn9YVYrotDu7UV7r8f/uRPlna1rJn7pGQnVlNiYiI333wzH/rQh6zJ3ZlPARMTE0xPT3P33XczNjZGeXk5TqeT5557jqioKGty97rrrqOlpQXTNBd9rSWHvlIqBnge+IrWekgpFQpsAK4BrgKeUUrlXW6nZ73OfmA/QFZW1pU+nQhy3/jGB0+lPPPhtrXVs0OASwe/hLxYCzabjcLCQgoLC61ttbW1vPzyy8TExFgBf+ONN1qTu++99x5Hjhyxyjsztf6FLCn0lVJ2PIF/QGv9gndzO/CC9tSHjiml3EAycBbYPOvHM73bWGS7RWv9BPAEeMo7S2mfEAtpa1v8/rExz45BQl34K4fDQVxcHKZp8s477/D222+Tl5fHvffey9atWxkZGaGiogLTNHnhhRd49dVXF32+S4a+8nzG+DFQr7X+7qy7XgRuAn6jlCoAwoBe4CXg50qp7+KZyM0HjgEKyFdK5eIJ+08B9y73DRBiObKyLn21rEvtGITwpZnlnQUFBQwNDVFRUUFvb6+1XLO5udk68dvp06dXpLxzPXA/UK2UqvBu+1vgSeBJpVQNMAk84B311yqlngHqgGngC1prF4BS6ovAQTxLNp/UWtcus/9CLMvDD1/6allSRRSBIi4ujj/4gz+wvp+YmODll19mamqKnJwcDMPgYx/72KLPIQdniXVPrpYl1rPh4WEqKipwOp309/cTFRXFX//1Xy+4emf1L+EixCo5cMCzjt5m8/x74MD8j7vvPs+aeq3hZz/zrNpRyvOvBL4IdLGxsdx444186Utf4v777yc3N3fRx8tIXwSktVw/L0dZi0Cz2BG5MtIXAWm+pZgzK3FW0szOpbXV80lhZpnnQp8qhPB3EvoiIC204malV+Ks1c5FiLUioS8C0kIrblZ6Jc5a7VyEWCsS+iIgPfywp4Y/W1SUZ/tKWqudixBrRUJfBKT77vNM2q72Spy12rkIsVaWdcI1IfzJWpwTR064JtYbCX0hLkFOuCbWEynvCMHSD/QSItDJSF8EvfnOub+UUy4LEYhkpC+CnqzFF8FEQl8EPVmLL4KJhP4KkrpwYJK1+CKYSOivEDlHS+CStfgimEjorxCpCweutTrQSwh/IKdWXiE229yLc8xQCtzutW+PECJ4yamV14DUhYUQgUBCf4VIXVgIEQgk9FeI1IWFEIFAjshdQXKOFiGEv5ORvhBCBBEJfSGECCIS+kIIEUQk9IUQIohI6AshRBCR0BdCiCAioS+EEEFEQl8IIYKIhL4QQgQRCX0hhAgiEvpCCBFEJPSFECKISOgLIUQQkdAXQoggcsnQV0ptVkr9RilVp5SqVUp9+aL7v6aU0kqpZO/3Sin1PaXUSaVUlVLKMeuxDyilmrxfD6x8d4QQQixmKefTnwa+prV2KqViAVMp9brWuk4ptRnYC7TNevyHgXzv1x7gB8AepdQG4JvAbkB7n+clrXX/CvZHCCHEIi450tdad2qtnd7bw0A9kOG9+xHgr/CE+IyPAj/VHu8BCUqpjcDtwOta6z5v0L8O3LFyXRFCCHEpy6rpK6VygDLgqFLqo8BZrXXlRQ/LAM7M+r7du22h7Re/xn6l1HGl1PGenp7lNE8IIcQlLDn0lVIxwPPAV/CUfP4W+O8r3SCt9RNa691a690pKSkr/fRCCBHUlhT6Sik7nsA/oLV+AdgC5AKVSqkWIBNwKqXSgbPA5lk/nundttB2IYQQa2Qpq3cU8GOgXmv9XQCtdbXWOlVrnaO1zsFTqnForbuAl4D/7F3Fcw0wqLXuBA4Ce5VSiUqpRDwTwAdXp1tCCCHms5TVO9cD9wPVSqkK77a/1Vq/ssDjXwH2ASeBMeBPAbTWfUqp/wm8733ct7TWfZfdciGEEMt2ydDXWr8NqEs8JmfWbQ18YYHHPQk8ubwmCiGEWClyRK4QQgQRCX0hhAgiEvpCCBFEJPSFECKISOgLIUQQkdAXQoggIqEvhBBBREJfCCGCiIS+EEIEEQl9IYQIIhL6QggRRCT0hRAiiEjoCyFEEJHQF0KIICKhL4QQQURCXwghgojyXPPEPymleoDWFX7aZKB3hZ9zNUg7V5a0c2UFQjsDoY2wOu3M1lqnzHeHX4f+alBKHdda7/Z1Oy5F2rmypJ0rKxDaGQhthLVvp5R3hBAiiEjoCyFEEAnG0H/C1w1YImnnypJ2rqxAaGcgtBHWuJ1BV9MXQohgFowjfSGECFoS+kIIEUTWdegrpTYrpX6jlKpTStUqpb7s3f5PSqkGpVSVUupXSqkEf2znrPu/ppTSSqlkf2yjUuovve9nrVLqH33VxsXaqZQqVUq9p5SqUEodV0pd7eN2RiiljimlKr3t/B/e7blKqaNKqZNKqV8qpcL8tJ0HlFInlFI1SqknlVJ2f2znrPu/p5Qa8VX7ZrVjofdTKaUeVko1KqXqlVJfWrVGaK3X7RewEXB4b8cCjcB2YC8Q6t3+D8A/+GM7vd9vBg7iOUgt2d/aCNwEHAbCvfel+uN7CRwCPuzdvg9408ftVECM97YdOApcAzwDfMq7/YfAg37azn3e+xTwC39tp/f73cDPgBFftvES7+efAj8FbN77Vu3vaF2P9LXWnVprp/f2MFAPZGitD2mtp70Pew/I9FUbYeF2eu9+BPgrwKcz7ou08UHg77XWE977zvmulYu2UwNx3ofFAx2+aaGH9pgZedq9Xxq4GXjOu/0nwF0+aJ5loXZqrV/x3qeBY/j+b2jediqlQoB/wvM35HOL/N4fBL6ltXZ7H7dqf0frOvRnU0rlAGV49qyz/Rnw6lq3ZyGz26mU+ihwVmtd6dNGXeSi97IAuNFbkvitUuoqX7Zttova+RXgn5RSZ4BvA1/3Xcs8lFIhSqkK4BzwOtAMDMwakLTz+52/z1zcTq310Vn32YH7gdd81b5ZbZmvnV8EXtJad/q2db+3QDu3AJ/0lh5fVUrlr9brB0XoK6VigOeBr2ith2Zt/wYwDRzwVdtmm91OPO36W+C/+7RRF5nnvQwFNuD5iPpfgWeUUsqHTQTmbeeDwH/RWm8G/gvwY1+2D0Br7dJal+IZJV8NFPq4SfO6uJ1KqeJZd38f+J3W+i3ftO735mnnHwB3A//s25bNtcD7GQ5c0J7TMfwIeHK1Xn/dh753JPI8cEBr/cKs7Z8G/gi4z/sR1afmaecWIBeoVEq14PkP4lRKpftRG8EzGn3B+7H1GODGcwIpn1mgnQ8AM7efxROyfkFrPQD8BrgWSFBKhXrvygTO+qxhF5nVzjsAlFLfBFKAr/qyXReb1c6bgK3ASe/fUJRS6qQv2zbbRe9nO7////krYNdqve66Dn3viPPHQL3W+ruztt+Bp8b3Ea31mK/aN6s9H2in1rpaa52qtc7RWufg+U/h0Fp3+UsbvV7E88eFUqoACMOHZzZcpJ0dwB96b98MNK1122ZTSqXMrBpTSkUCt+GZf/gN8Anvwx4A/s03LfRYoJ0NSqm/AG4H7pmpQ/vSAu00tdbps/6GxrTWW/2wnQ3M+jvC8/+0cdXa4AeD3FWjlLoBeAuoxjMCBU/J5Ht4Pk6d9257T2v9ubVvocdC7dRavzLrMS3Abq21TwJ1kffyMJ6PoqXAJPCQ1vrXvmgjLNrOIeBRPOWoC8DntdamTxoJKKV24ZmoDcEz+HpGa/0tpVQe8DSeklk58Cczk+R+1s5pPCvKhr0PfUFr/S0fNXPBdl70mBGtdYwv2jerDQu9nwl4ysxZwAjwudWay1vXoS+EEGKudV3eEUIIMZeEvhBCBBEJfSGECCIS+kIIEUQk9IUQIohI6AshRBCR0BdCiCDy/wMxjo/8BnmMTQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "1. Looking at the shapes of the different separators, which one do you think performed best given the training data? Why?"
   ],
   "metadata": {
    "id": "vA8kUFalWfmd"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "2. Note that the variable ```no_samples``` defined in the first cell defines the number of samples. Increase the value of samples to 3000 and run all the script again to see what happens. What are your main conclusions of this experiment?"
   ],
   "metadata": {
    "id": "DpVVL_BRU22g"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Set again the number of samples to 30 to answer the rest of the questions.\n",
    "3. Note that in the **Using Linear Programming Section**, the variable ```theta``` has the parameters of the first linear separator, and the variable ```theta_prime``` the parameter of the second linear separator. Recall that the margin is defined as $\\omega = \\frac{|(\\hat\\theta_0+1)-(\\hat\\theta_0-1)|}{\\sqrt{\\hat\\theta_1^2 + \\hat\\theta_2^2}} = \\frac{2}{\\left\\|{\\hat{\\theta}}\\right\\|}$º. Calculate the margin of the two linear separators. Which one has the largest margin? Why?"
   ],
   "metadata": {
    "id": "kqGzHFvyW6l8"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "print(\"Parameters of NLP estimator:\")\n",
    "print(theta)\n",
    "print(\"Parameters of NLP estimator with regularization:\")\n",
    "print(theta_prime)\n",
    "\n"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "gMybjiyIV3vd",
    "outputId": "8293952e-cf38-4df0-f2d6-f82c39c047ef",
    "pycharm": {
     "is_executing": true
    }
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "Lcxcvh30dliS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "4. Now, imagine we want to predict the status of the machine for the following points, given as pairs of temperature, rpms  (temperature, rpm): \n",
    "- (32, 2600)\n",
    "- (34, 2600)\n",
    "- (36, 2600)\n",
    "\n",
    "Note that to get the results with the two NLP linear separators, you need to calculate if the point is over or below the separators as, as: \n",
    "\n",
    "- ```theta[0]*temp + theta[1]*rpm + theta[2] >= 1``` → predict failure\n",
    "\n",
    "- ```theta[0]*temp + theta[1]*rpm + theta[2] <= -1``` → predict normal state$\n",
    "\n",
    "Test both linear separators in these provided points. Looking at the data, which one do you think provides more reliable results?"
   ],
   "metadata": {
    "id": "49toM4WrbvqC"
   }
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "QTH8bX6udmzC",
    "outputId": "3213c490-7fb8-4442-b595-88ac0f5ba960",
    "pycharm": {
     "is_executing": true
    }
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "5. Now, to predict values with the Scipy separators, you need to use the method ```predict()``` of the predictor. For instance, the following script uses the SVM predictor to predict the values at 24 degrees and 2200 rpms and 28 degrees and 3300 rpms:"
   ],
   "metadata": {
    "id": "zNioydGXeyey"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "res = clf_svm.predict([[24, 2200], [28, 3300]])\n",
    "print(res)"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "2gA_WRsSe_vP",
    "outputId": "cfb42b6d-b358-4554-c2a1-260464b37b0e"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "[-1.  1.]\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "Write down a script to predict the values at the same points as in exercise 4, but using the Scipy models you have trained. What are the main conclusions?"
   ],
   "metadata": {
    "id": "7IBih8-tf2TB"
   }
  }
 ]
}