svm_IPmethods.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "CO Project.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/msdkhairi/svm/blob/master/svm_IPmethods.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qXtiqF4BJN15",
        "colab_type": "text"
      },
      "source": [
        "In this notebook **Interior Point Methods** are used to solve SVM. Two methods are implemented:\n",
        "\n",
        "*   Barrier method\n",
        "*   PDIP (Primal-Dual Interior Point method)\n",
        "\n",
        "For more information about these methods, you can refer to chapter 11 of Convex Optimization book by Stephen Boyd.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kfYaM5WPaybI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GuU1MzekFubS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# load data\n",
        "data = np.loadtxt('data.txt', delimiter=' ')\n",
        "X = data[:,:2]\n",
        "y = data[:,2]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gyWOyxMkpKYm",
        "colab_type": "text"
      },
      "source": [
        "# Barrier method\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AwtyEGubYg-B",
        "colab_type": "text"
      },
      "source": [
        "### RBF kernel"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hoLtf_oCmy7m",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# define functions to compute RBF kernel of data\n",
        "\n",
        "def rbf(x1, x2, gamma=500):\n",
        "  axis = None if x2.ndim==1 else 1\n",
        "  return np.exp(-gamma*np.linalg.norm(x1-x2, axis=axis)**2)\n",
        "\n",
        "\n",
        "def RBF_kernel(X, y):\n",
        "    \n",
        "  x_size = X.shape[0]\n",
        "  K = np.zeros((x_size, x_size))\n",
        "\n",
        "  for i in range(x_size):\n",
        "    for j in range(x_size):\n",
        "      K[i, j] = y[i] * y[j] * rbf(X[i], X[j])\n",
        "  return K\n",
        "\n",
        "def log(w):\n",
        "  \"\"\" return the log(w) if w > 0 \n",
        "      when w <= a very big negative number is returned\n",
        "  \"\"\"\n",
        "  w = (w > 0) * w + (w <= 0) * 1e-300\n",
        "  return np.log(w)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9HC31BkNqYQJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# compute RBF kernel\n",
        "K_hat= RBF_kernel(X, y)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kbICI8pdpHAO",
        "colab_type": "text"
      },
      "source": [
        "### Barrier Objective "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vsXTot2w_-WJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class Objective:\n",
        "    def __init__(self, w_0, t=10**4, K=K_hat, C=0.1):\n",
        "        self.w = w_0\n",
        "        self.t = t\n",
        "        self._K = K\n",
        "        self._C = C\n",
        "        \n",
        "        self._ones = np.ones_like(w_0)\n",
        "        self._update()\n",
        "    \n",
        "    def _update(self):\n",
        "        self.f = self._f(self.w)\n",
        "        self.grad = self.gradient(self.w)\n",
        "        self.hess = self.hessian(self.w)\n",
        "    \n",
        "    def _f(self, w):\n",
        "        f_hat = 0.5 * w.T @ self._K @ w - np.sum(w)\n",
        "        g = -np.sum(log(w)) - np.sum(log(self._C*self._ones - w))\n",
        "        return self.t*f_hat + g\n",
        "    \n",
        "    def gradient(self, w):\n",
        "        return self.t * self._K @ w - self.t * self._ones - (\n",
        "            np.divide(self._ones, w) + np.divide(self._ones, (w-self._C*self._ones)))\n",
        "    \n",
        "    def hessian(self, w):\n",
        "        return self.t * self._K + np.diag(np.divide(self._ones, w**2)) + (\n",
        "            np.diag(np.divide(self._ones, (w-self._C*self._ones)**2)))\n",
        "        \n",
        "    def setw(self, w):\n",
        "        self.w = w\n",
        "        self._update()\n",
        "        \n",
        "    def sett(self, t):\n",
        "        self.t = t\n",
        "        self._update()\n",
        "        \n",
        "    def __call__(self, w):\n",
        "        return self._f(w)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pHXYMgCypSsu",
        "colab_type": "text"
      },
      "source": [
        "### Newton's method"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rnXeEYJDCxZs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class Newton:\n",
        "    \n",
        "    def __init__(self, objective, A=y.reshape(-1, 1).T, eps=1e-16 ,alpha=0.01, beta=0.5):\n",
        "        \n",
        "        self.f = objective\n",
        "        self._A = A\n",
        "        self._epsilon = eps\n",
        "        \n",
        "        # params of backtracking line search\n",
        "        self._alpha = alpha\n",
        "        self._beta = beta\n",
        "        \n",
        "        self.dw = self._newton_step()\n",
        "        self.decrement = self._newton_decrement()\n",
        "    \n",
        "    def _newton_step(self):\n",
        "        # solving the KKT system to obtain newton step\n",
        "        \n",
        "        m = np.vstack((self.f.hess, self._A))\n",
        "        m = np.hstack((m, np.vstack((self._A.T, 0))))\n",
        "        \n",
        "        b = np.hstack((-self.f.grad, 0))\n",
        "        \n",
        "        dw_nu = np.linalg.inv(m) @ b\n",
        "        dw = dw_nu[:-1]\n",
        "        return dw\n",
        "    \n",
        "    def _newton_decrement(self):\n",
        "        # newton decrement ** 2\n",
        "        return self.dw.T @ self.f.hess @ self.dw\n",
        "    \n",
        "    def _line_search(self):\n",
        "        # back_tracking line search\n",
        "        t = 1\n",
        "        while self.f(self.f.w+t*self.dw) > self.f.f + (\n",
        "            self._alpha * t * self.f.grad.T @ self.dw):\n",
        "            t = self._beta * t\n",
        "        return t\n",
        "        \n",
        "    def run(self):\n",
        "        while self.decrement/2 > self._epsilon:\n",
        "            t = self._line_search()\n",
        "            self.f.setw(self.f.w + t*self.dw)            \n",
        "            self.dw = self._newton_step()\n",
        "            self.decrement = self._newton_decrement()\n",
        "        return self.f.w\n",
        "            \n",
        "    def setobj(self, obj):\n",
        "        self.objective = obj"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WFw8jyIhvXyt",
        "colab_type": "text"
      },
      "source": [
        "### Barrier method"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XsL1Y2jSqS7C",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class Barrier:\n",
        "    def __init__(self, objective, newton, t=10000, m=y.shape[0]*2, mu=4, eps=1e-8):\n",
        "        self.f = objective\n",
        "        self.newton = newton\n",
        "        self.t = t\n",
        "        self._m = m\n",
        "        self._mu = mu\n",
        "        self._eps = eps\n",
        "        \n",
        "        self.obj_values = []\n",
        "        \n",
        "    def update(self, w):\n",
        "        # update w, t, and f\n",
        "        self.f.setw(w)\n",
        "        self.f.sett(self.t)\n",
        "        self.newton.setobj(self.f)\n",
        "                \n",
        "    def report(self, plot=False):\n",
        "        print('Iter {}: f* = {}'.format(self.iteration+1, self.f.f))\n",
        "        \n",
        "    def run(self):\n",
        "        self.iteration = 0\n",
        "        while self._m/self.t > self._eps:\n",
        "            self.report()\n",
        "            self.obj_values.append(self.f.f)\n",
        "            \n",
        "            w_star = self.newton.run()\n",
        "            self.t = self._mu * self.t\n",
        "            self.update(w_star)\n",
        "            \n",
        "            self.iteration += 1\n",
        "            "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IyaxBs-aV6B_",
        "colab_type": "text"
      },
      "source": [
        "### SVM"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_eezeJkYV00C",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class SVMClassifier:\n",
        "    def __init__(self, w, K, X, y):\n",
        "        self.w = w\n",
        "        self.K = K\n",
        "        self.X = X\n",
        "        self.y = y\n",
        "        \n",
        "        self.b = self._beta_0()\n",
        "        \n",
        "    def _beta_0(self):\n",
        "        sup_vec = np.where(self.w > 1e-5)[0]\n",
        "        b = 0\n",
        "        for j in sup_vec:\n",
        "            b += self.y[j] - np.sum(self.w*self.y*self.K[:,j])\n",
        "        return b / sup_vec.size\n",
        "    \n",
        "    def predict(self, x):\n",
        "        size = x.shape[0]\n",
        "        y_pred = np.zeros(size)\n",
        "        if x.ndim==1:\n",
        "            return np.sign(np.sum(self.w*self.y*rbf(x, self.X)) + self.b)\n",
        "        else:\n",
        "            for i in range(size):\n",
        "                y_pred[i] = self.predict(x[i])\n",
        "            return y_pred\n",
        "        \n",
        "    def accuracy(self, x, y_true):\n",
        "        y_pred = self.predict(x)\n",
        "        return np.sum(y_pred==y_true) / y_true.size"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FdTHyOaZwh2h",
        "colab_type": "text"
      },
      "source": [
        "### RUN the optimiztion"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IRPM_9N3xX5M",
        "colab_type": "code",
        "outputId": "4a5ddf4c-a98c-4040-9382-8968636f2f43",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "source": [
        "# define initial point w_0\n",
        "w_0 = (y==-1) * np.count_nonzero(y==1) + (y==1) * np.count_nonzero(y==-1)\n",
        "w_0 = w_0 / 10000\n",
        "np.dot(w_0, y)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "-9.33975119465913e-15"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 210
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T7WF7ZjXThhw",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "objective = Objective(w_0, t=10**4)\n",
        "newton = Newton(objective, eps=1e-5)\n",
        "barr = Barrier(objective, newton)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oiLgahQ26e0q",
        "colab_type": "code",
        "outputId": "cc306816-d840-49c0-f7fb-efc067c4be73",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 254
        }
      },
      "source": [
        "barr.run()"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Iter 1: f* = -276649.1034325002\n",
            "Iter 2: f* = -362773805.3551815\n",
            "Iter 3: f* = -1452847330.3718767\n",
            "Iter 4: f* = -5813141430.438658\n",
            "Iter 5: f* = -23254317830.70578\n",
            "Iter 6: f* = -93019023431.77428\n",
            "Iter 7: f* = -372077845836.0482\n",
            "Iter 8: f* = -1488313135453.144\n",
            "Iter 9: f* = -5953254293921.527\n",
            "Iter 10: f* = -23813018927795.062\n",
            "Iter 11: f* = -95252077463289.2\n",
            "Iter 12: f* = -381008311605265.75\n",
            "Iter 13: f* = -1524033248173172.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kDbSO2ZH1ASq",
        "colab_type": "code",
        "outputId": "b636a0ad-ba35-4779-ccb6-f1eb914eb79d",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 281
        }
      },
      "source": [
        "plt.figure(), plt.title('f* in iterations')\n",
        "plt.plot(barr.obj_values);"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmUXOV95vHv06uWLiSEWtXaQGC6\n25YxW9oymDG2QSSAFwFjiJkskMSHOYmdOIkJg81kxuPM5OA48ZKD44wCGXDMiY0xMhjjYBDG2CdA\nLAEyErIWs2pvSYDUElKru3/zR90WJbk3dVXXreX5nNOn7731Vr2/a+R6uu59630VEZiZWe2pS7sA\nMzNLhwPAzKxGOQDMzGqUA8DMrEY5AMzMapQDwMysRjkArKxI6pT0jKS9kv7kGJ63RtL7xtlnj6RT\nxvPcYpD0Hknr0urfapf8PQArJ5JuA/ZExJ/lHVsAXBsRny1B/7cDmyLiv09gHwG0R8TGierDbCz8\nCcDKzUnAGgBJ50i6CWhI9s9P9suWpIa0azAbKweAlQ1JjwDvB26R1APsBlYDXwM+ClwCfGWY574o\naXGy/VlJd0n6enIpaY2krhH6DUmnSroO+C3ghuSy0PeSx+dI+o6kbkkv5F+aSvq6W9I3JO0BrpW0\nSNLjkl6TtFXSLZKakvaPJU9dlfTxm5LeJ2lT3mu+TdKjyfPXSPpw3mO3S/qqpO8n5/akpLckj0nS\nlyTtkLRH0rOSTjvW/w5WOxwAVjYi4gLgJ8AnIqIlItYD+dco+4/aH8mHgW8C04H7gFvG0P9S4E7g\nb5L+PySpDvgesAqYC1wI/Kmk38h76hLg7qSvO5M6/wyYCZybPOePkj7OT55zRtLHt/JrkNSY9PdD\nYBbwx8Cdkjrzmn0U+F/A8cBG4P8kx38dOB/oAKYBVwG7Rjtvq11lHwCS/jn5i2b1GNqeL+kpSX2S\nPnLUY/3JzcVnJN03cRVbsUg6Bzgd+ENyb+YPAZ8c49N/GhEPREQ/8C/AGeMs451Aa0R8LiJ6I+J5\n4J/IvQkPejwivhsRAxHxRkSsjIgnIqIvIl4E/i/w3jH2dw7QAtyc9PcIcD9wdV6bZRHxHxHRRy5w\nzkyOHwIywFvJ3d9bGxFbx3neVgMq4Xrl7eT+evv6GNq+DFwLXD/EY29ExJlDHLcyFRFPAE8kN4GJ\niB8DPx7j07flbe8HJklqSN40j8VJwBxJr+Udqyf3SWXQK/lPkNQBfBHoAqaQ+//ZyjH2Nwd4JSIG\n8o69RO7Tx6Cjz60FICIekXQL8FXgJEn3ANdHxJ4x9m01puw/AUTEY+SuBR8m6S2S/k3SSkk/kfTW\npO2LEfFzYGCo17LKlPx3/Wypujtq/xXghYiYnveTiYhLR3jO14BfkBvpcxzwGUBj7H8LMD+59DTo\nRGDzmIqP+PuI+DVgIblLQX8xxn6tBpV9AAxjKfDHyT/064F/GMNzJklaIekJSZdNbHlWwbYD+d8J\n+A9gr6T/JmmypHpJp0l65wivkQH2AD3JHyd/OEof+Z4k91f9DZIak+82fIjcJbARSXqnpHcl9xH2\nAQfwH0M2gooLAEktwLuBb0t6htz11dljeOpJEdEF/Bfgy4MjJ8yOchuwMBmB893kHsIHyV1nfwHY\nCdxK7ibrcK4n9+9sL7n7Bd866vHPAnckfVyV/0BE9JJ7w78k6esfgN+NiF+Mofbjkv5eJXfZaBfw\nhTE8z2pURXwRLLkGfH9EnCbpOGBdRAz7pp98mef+iLh7PI+bmdWCivsEkNzQekHSlXB47POIIzwk\nHS+pOdmeCZwHPDfhxZqZlbGy/wQg6V+B95EbU70d+J/AI+RutM0GGoFvRsTnkuuyy8iNjz4AbIuI\nt0t6N7lLRQPkQu/LEXFbqc/FzKyclH0AmJnZxKi4S0BmZlYcZf1FsJkzZ8aCBQvSLsPMrGKsXLly\nZ0S0jqVtWQfAggULWLFiRdplmJlVDEkvjbWtLwGZmdUoB4CZWY1yAJiZ1SgHgJlZjXIAmJnVqKIE\ngKSLJa2TtFHSjUM83izpW8njTw7O725mZukpOAAk1ZNbgOIScnOQXy1p4VHN/gB4NSJOBb4EfL7Q\nfs3MrDDF+B7AImBjslQekr5Jbo3U/MnWlpCbAhdya6feIkkxQfNQ/P3yDfT1exr0Y6KxrldSWuOp\nKv9UlLzCkceGaDfE+Q8eUl4Vbx7LbdclB+qkw/tS7vWUd3xwm8E2QF1d7rUH29cpt18nDh9rqq+j\nsb6OxnrR2FBHY10djQ2isb7uiMcaDu+L+joNeT5mRytGAMzlyCXxNgHvGq5NRPRJeh04gdx850eQ\ndB1wHcCJJ544roL+8ce/5I1D/eN6bi3ydFDVReJwQDTU54eFksA4cvv4qY3MbGmmtaWZmZnm3Ham\nmZktTcxsaWZSY33ap2QTpOy+CRwRS8mt+EVXV9e43pqe+9zFRa3JKkf+h8rBzRji8SOPDbb71ecy\nTLuI3GtEBAORe8GByL3CQCSPH7U/+JujjyWvk9vPvf7AQO6xvoHgUP8Ah/oG6O0f4FB/sn/Udm/f\nm/t9/QP0HtFugN6+oG/gze38563btpef7t3JngNDL5ecmdSQBEIuJPLDoTUvME5oaaK5wWFRSYoR\nAJuB+Xn78/jV9UsH22yS1EBuNaVdRejb7Aj5lz6GvgriSyPDOdjXz86eXnbuPcjOnoN0J7939vTS\nvfcg3T0HWbt1D49tOMjeYcLiuPywyPvdmc1wfkcrTQ0eeFhOihEAPwPaJZ1M7o3+o+SWw8t3H3AN\n8DjwEeCRibr+b2bj09xQz9zpk5k7ffKobQ8c6j8iHHb2HGRnEhK57V7WbNnDzr0H2XswFxbHT2nk\nQ2fM4Yqz53HGvGm+T1EGCg6A5Jr+J4AHgXrgnyNijaTPASsi4j5y66z+i6SNwG5yIWFmFWpSYz3z\njp/CvOOnjNr2jd5+nnh+F995ahPf+tkrfP3xlzhl5lSuOHsul501d0yvYROjrBeE6erqCs8GalY9\n9hw4xA+e3co9T23myRd2A/Cuk2dwxdlzueQdszluUmPKFVY+SSsjomtMbR0AZpaGV3bv595nNnPP\nU5t5fuc+mhvquGhhlivOnst72ltprPf9gvFwAJhZxYgIVm16nXue2sT3Vm3h1f2HmNnSxIfPmMsV\nZ8/l7XOO8/2CY+AAMLOK1Ns3wKPrdrDs6c0sX7uD3v4BOrItXH7WPC47aw6zp41+g7rWOQDMrOK9\nvv8Q9z+7hXue2szKl15Fgne/5QSuOGseF5/WxtTmsvsaU1lwAJhZVXlx5z6WPb2ZZU9v5uXd+5nc\nWM/Fp7Vx+VlzOe/UmdTX+RLRIAeAmVWliGDlS69yz9ObuX/VFvYc6CN7XDNLzpzLR35tHh3ZTNol\nps4BYGZV78Chfn70ix1856nNPLpuBwMR3P57izi/ozXt0lJ1LAHgcVZmVpEmNdZzyTtmc+s1XTz5\nmQs5pbWFG+7+Oa/vP5R2aRXDAWBmFe+Elma+dNWZ7Ow5yGe/tybtciqGA8DMqsI75k3jExecyrKn\nN/ODZ7emXU5FcACYWdX4+PtP5R1zp/GZZc/Svfdg2uWUPQeAmVWNxvo6vnjVGezr7efT9zxLOQ9y\nKQcOADOrKu3ZDDf8RicPr93O3Ss3pV1OWXMAmFnV+f3zTuZdJ8/gc997js2vvZF2OWXLAWBmVaeu\nTvztlWcwEMFffHsVAwO+FDQUB4CZVaX5M6bwlx9cyL//chd3PP5i2uWUJQeAmVWt33znfC546yxu\n/sEv+GV3T9rllJ2CAkDSDEkPSdqQ/D5+iDZnSnpc0hpJP5f0m4X0aWY2VpK4+Yp3MLmpnj+/axV9\n/QNpl1RWCv0EcCOwPCLageXJ/tH2A78bEW8HLga+LGl6gf2amY3JrOMm8b8vO41Vr7zGP/74l2mX\nU1YKDYAlwB3J9h3AZUc3iIj1EbEh2d4C7ABqe7YmMyupD54+hw+dMYcvP7yB1ZtfT7ucslFoAGQj\nYvA719uA7EiNJS0CmoBhY1jSdZJWSFrR3d1dYHlmZjl/teTtzJjaxKfuWsXBvv60yykLowaApIcl\nrR7iZ0l+u8h95W7YsVaSZgP/AvxeRAx7IS4ilkZEV0R0tbb6g4KZFcf0KU18/j+fzrrte/niQ+vT\nLqcsjLqmWkQsHu4xSdslzY6Irckb/I5h2h0HfB+4KSKeGHe1ZmYFeP9bZ3H1ovksfex5Lnpblq4F\nM9IuKVWFXgK6D7gm2b4GuPfoBpKagGXA1yPi7gL7MzMryE0fWMi84yfzqW+vYt/BvrTLSVWhAXAz\ncJGkDcDiZB9JXZJuTdpcBZwPXCvpmeTnzAL7NTMbl5bmBv72I2fw8u79/PUDa9MuJ1WjXgIaSUTs\nAi4c4vgK4GPJ9jeAbxTSj5lZMb3rlBP42H86mX/6yQv8+tvbeG+NLiPpbwKbWU361K930j6rhRvu\nXlWzy0g6AMysJk1qrOeLV53Jrp5e/ud9q9MuJxUOADOrWYPLSH73mS08UIPLSDoAzKymffz9p3L6\nvGnctOxZduw9kHY5JeUAMLOalr+M5GdqbBlJB4CZ1bxTZw0uI7mDb9fQMpIOADMzjlxGctOr+9Mu\npyQcAGZmvLmMZERwfY0sI+kAMDNLzJ8xhf/xoYU88fxubv/3F9MuZ8I5AMzM8lzVlVtG8vP/9gs2\n7qjuZSQdAGZmefKXkfzUt6t7GUkHgJnZUfKXkfzao9W7jKQDwMxsCIPLSH5lefUuI+kAMDMbxuAy\nkn9+1zMcOFR9y0g6AMzMhjF9ShOf/8jprN/ew5eqcBlJB4CZ2Qje3zmLqxedyNKfPM/PXtyddjlF\n5QAwMxvFTR94W24ZybuqaxnJggNA0gxJD0nakPw+foS2x0naJOmWQvs1MyuVluYG/u7KM3nl1f3c\n9tMX0i6naIrxCeBGYHlEtAPLk/3h/BXwWBH6NDMrqUUnz6BjVoafb3ot7VKKphgBsAS4I9m+A7hs\nqEaSfg3IAj8sQp9mZiXXnm1h3fa9aZdRNMUIgGxEDC6ls43cm/wRJNUBfwdcP9qLSbpO0gpJK7q7\nu4tQnplZcXRmM7yy+w3291bHfYAxBYCkhyWtHuJnSX67yK2kMNQUen8EPBARo060HRFLI6IrIrpa\nW1vHdBJmZqXQ0ZYBYMP26pgjqGEsjSJi8XCPSdouaXZEbJU0G9gxRLNzgfdI+iOgBWiS1BMRI90v\nMDMrK53ZXACs276XM+ZPT7mawo0pAEZxH3ANcHPy+96jG0TEbw1uS7oW6PKbv5lVmvkzptDcUMf6\nbdVxH6AY9wBuBi6StAFYnOwjqUvSrUV4fTOzslBfp6q6EVzwJ4CI2AVcOMTxFcDHhjh+O3B7of2a\nmaWhI5vh3zfuSruMovA3gc3MjkFnNsO2PQd4ff+htEspmAPAzOwYDI4EWr+j8i8DOQDMzI5Bx+BI\noCq4EewAMDM7BnOmTaKluYH1VXAj2AFgZnYMJNGRbXEAmJnVos62DOu27SU3+UHlcgCYmR2jjmyG\nV/cfYmdPb9qlFMQBYGZ2jAZvBFf6ZSAHgJnZMaqWkUAOADOzYzSzpYkZU5vYUOHfBXAAmJkdo8GR\nQP4EYGZWgzqzGdZv76nokUAOADOzcWjPZug52MeW1w+kXcq4OQDMzMahc3BOoAq+DOQAMDMbh45Z\nlT8U1AFgZjYO06Y00nbcpIpeHMYBYGY2Th1tmdr9BCBphqSHJG1Ifh8/TLsTJf1Q0lpJz0laUEi/\nZmbloGNWCxu299A/UJkjgQr9BHAjsDwi2oHlyf5Qvg58ISLeBiwCdhTYr5lZ6jraMhzsG+Dl3fvT\nLmVcCg2AJcAdyfYdwGVHN5C0EGiIiIcAIqInIirzfy0zszydFT4nUKEBkI2Ircn2NiA7RJsO4DVJ\n90h6WtIXJNUP94KSrpO0QtKK7u7uAsszM5s47dkWoHKHgjaM1kDSw0DbEA/dlL8TESFpqAthDcB7\ngLOAl4FvAdcCtw3VX0QsBZYCdHV1VeaFNTOrCVOaGjhxxpSKHQk0agBExOLhHpO0XdLsiNgqaTZD\nX9vfBDwTEc8nz/kucA7DBICZWSWp5NXBCr0EdB9wTbJ9DXDvEG1+BkyX1JrsXwA8V2C/ZmZloSOb\n4fnuffT2DaRdyjErNABuBi6StAFYnOwjqUvSrQAR0Q9cDyyX9Cwg4J8K7NfMrCx0tmXoGwhe3LUv\n7VKO2aiXgEYSEbuAC4c4vgL4WN7+Q8DphfRlZlaO8heHGdyuFP4msJlZAU5pnUp9nSryPoADwMys\nAM0N9Sw4YUpFLg7jADAzK1Bnhc4J5AAwMytQRzbDS7v3c+BQf9qlHBMHgJlZgTqzGSJg446etEs5\nJg4AM7MCteeNBKokDgAzswItOGEKTfV1FXcfwAFgZlaghvo63jKrpeLmBHIAmJkVQWe2peJmBXUA\nmJkVQUdbhi2vH2DvgUNplzJmDgAzsyLomDW4OEzljARyAJiZFUFnW+WtDuYAMDMrgrnTJzOlqb6i\nhoI6AMzMiqCuTrRnK2tKCAeAmVmRdGZbfA/AzKwWdWQz7Ow5yK6eg2mXMiYOADOzIhlcEKZSPgUU\nHACSZkh6SNKG5Pfxw7T7G0lrJK2V9PeSVGjfZmblpNJGAhXjE8CNwPKIaAeWJ/tHkPRu4Dxyy0Ke\nBrwTeG8R+jYzKxuzMs1Mm9xYMVNCFCMAlgB3JNt3AJcN0SaASUAT0Aw0AtuL0LeZWdmQRGc2w4Ya\nCoBsRGxNtrcB2aMbRMTjwI+ArcnPgxGxdqgXk3SdpBWSVnR3dxehPDOz0mnPtrBu214iIu1SRjWm\nAJD0sKTVQ/wsyW8XuTP+lbOWdCrwNmAeMBe4QNJ7huorIpZGRFdEdLW2th7zCZmZpamzLcOeA31s\n31P+I4EaxtIoIhYP95ik7ZJmR8RWSbOBHUM0uxx4IiJ6kuf8ADgX+Mk4ajYzK1uDI4HWbd9L27RJ\nKVczsmJcAroPuCbZvga4d4g2LwPvldQgqZHcDeAhLwGZmVWyw0NBK2BKiGIEwM3ARZI2AIuTfSR1\nSbo1aXM38EvgWWAVsCoivleEvs3MysqMqU20ZporYijomC4BjSQidgEXDnF8BfCxZLsf+K+F9mVm\nVgk6si0VEQD+JrCZWZF1ZDOs397DwEB5jwRyAJiZFVlnNsMbh/rZ9OobaZcyIgeAmVmRdbS9ORKo\nnDkAzMyKrH1WC1D+cwI5AMzMiiwzqZG50yc7AMzMalFHMiVEOXMAmJlNgI62DM937+NQ/0DapQzL\nAWBmNgE6sxl6+wd4ade+tEsZlgPAzGwCVMLqYA4AM7MJcOqsFupEWd8HcACYmU2ASY31nHTC1LIe\nCeQAMDObIB3ZlrL+MpgDwMxsgnRmM7y4cx8HDvWnXcqQHABmZhOkoy3DQMDz3eU5EsgBYGY2Qd4c\nCVSel4EcAGZmE2TBCVNprFfZ3gdwAJiZTZCmhjpOmdlStstDFhQAkq6UtEbSgKSuEdpdLGmdpI2S\nbiykTzOzStLRlqnaTwCrgSuAx4ZrIKke+CpwCbAQuFrSwgL7NTOrCB2zWtj06hvsO9iXdim/oqAA\niIi1EbFulGaLgI0R8XxE9ALfBJYU0q+ZWaUYXBxmw47ymxKiFPcA5gKv5O1vSo4NSdJ1klZIWtHd\n3T3hxZmZTaTOwZFAZXgfoGG0BpIeBtqGeOimiLi32AVFxFJgKUBXV1d5r6hsZjaK+TOmMKmxrizv\nA4waABGxuMA+NgPz8/bnJcfMzKpefZ1on5Upy+8ClOIS0M+AdkknS2oCPgrcV4J+zczKQnu2pfoC\nQNLlkjYB5wLfl/RgcnyOpAcAIqIP+ATwILAWuCsi1hRWtplZ5ejMZti+5yCv7e9Nu5QjjHoJaCQR\nsQxYNsTxLcClefsPAA8U0peZWaUaHAm0fnsPi06ekXI1b/I3gc3MJtjgSKByuxHsADAzm2Czp00i\n09xQdkNBHQBmZhNMUlneCHYAmJmVQGdbbihoRPl8vckBYGZWAh3ZDK/uP0R3z8G0SznMAWBmVgJv\nTglRPnMCOQDMzEpgcChoOY0EcgCYmZXAzJZmZkxtYoMDwMys9nRkW/wJwMysFnVmM6zfVj4jgRwA\nZmYl0tGWYV9vP5tfeyPtUgAHgJlZyRweCVQml4EcAGZmJdKefXNSuHLgADAzK5FpkxtpO25S2cwJ\n5AAwMyuhjrZM2YwEcgCYmZVQZ7aFDTt66B9IfySQA8DMrIQ6shl6+wZ4ade+tEspeEnIKyWtkTQg\nqWuYNvMl/UjSc0nbTxbSp5lZJesooxvBhX4CWA1cATw2Qps+4FMRsRA4B/i4pIUF9mtmVpHasy1A\neQwFLXRN4LWQW+xghDZbga3J9l5Ja4G5wHOF9G1mVommNDVw4owpZXEjuKT3ACQtAM4CnhyhzXWS\nVkha0d3dXarSzMxKpiOZEiJtowaApIclrR7iZ8mxdCSpBfgO8KcRsWe4dhGxNCK6IqKrtbX1WLow\nM6sInW0tvLBzH719A6nWMeoloIhYXGgnkhrJvfnfGRH3FPp6ZmaVrCOboW8geGHnPjqTdQLSMOGX\ngJS7QXAbsDYivjjR/ZmZlbvBkUBp3wcodBjo5ZI2AecC35f0YHJ8jqQHkmbnAb8DXCDpmeTn0oKq\nNjOrYKe0TqW+TqnfByh0FNAyYNkQx7cAlybbPwWGHyZkZlZjmhvqOXnm1Mr+BGBmZuPTkW1J/bsA\nDgAzsxR0ZDO8vHs/b/T2p1aDA8DMLAWd2QwRsHFHelNCOADMzFLQ0Zb+SCAHgJlZCk6aMYWmhrpU\n7wM4AMzMUtBQX8dbWltYl+JQUAeAmVlKOrMtbPAnADOz2tPRlmHL6wfYc+BQKv07AMzMUtKZTAmR\n1qcAB4CZWUoOzwm0LZ2hoA4AM7OUzJ0+mSlN9amNBHIAmJmlpK5OtGczDgAzs1rUmeKcQA4AM7MU\ndWQz7OzpZWfPwZL37QAwM0vR4IpgaXwKcACYmaVocCRQGovDOADMzFI0K9PMtMmNrE9hVtBCl4S8\nUtIaSQOSukZpWy/paUn3F9KnmVk1kURnNlORnwBWA1cAj42h7SeBtQX2Z2ZWdTraWli3fS8RUdJ+\nCwqAiFgbEetGaydpHvAB4NZC+jMzq0ad2Qx7D/Sxbc+BkvZbqnsAXwZuAAZGayjpOkkrJK3o7u6e\n+MrMzFLWfnhKiNJeBho1ACQ9LGn1ED9LxtKBpA8COyJi5VjaR8TSiOiKiK7W1taxPMXMrKJ1HJ4U\nrrQ3ghtGaxARiwvs4zzgw5IuBSYBx0n6RkT8doGva2ZWFWZMbaI101zy5SEn/BJQRHw6IuZFxALg\no8AjfvM3MztSZwpzAhU6DPRySZuAc4HvS3owOT5H0gPFKNDMrBZ0JAEwMFC6kUCjXgIaSUQsA5YN\ncXwLcOkQxx8FHi2kTzOzatSRbeHAoQFeeXU/J50wtSR9+pvAZmZloOPwnECluxHsADAzKwPts1qA\n0k4K5wAwMysDmUmNzJ0+uaTfBXAAmJmVic620o4EcgCYmZWJ9mwLv+zu4VD/qJMmFIUDwMysTHRm\nMxzqD17ata8k/TkAzMzKRMfhOYFKMxLIAWBmViZOndVCnSjZlBAOADOzMjGpsZ6TTphassVhHABm\nZmWkI9tSspFABU0FYWZmxfW+zlnMmNpERCBpQvtyAJiZlZGrF53I1YtOLElfvgRkZlajHABmZjXK\nAWBmVqMcAGZmNcoBYGZWowpdEvJKSWskDUjqGqHddEl3S/qFpLWSzi2kXzMzK1yhnwBWA1cAj43S\n7ivAv0XEW4EzgLUF9mtmZgUqdE3gtcCIX1aQNA04H7g2eU4v0FtIv2ZmVrhSfBHsZKAb+H+SzgBW\nAp+MiCHnO5V0HXBdstsjad04+50J7Bznc8tNtZxLtZwH+FzKUbWcBxR2LieNtaEiYuQG0sNA2xAP\n3RQR9yZtHgWuj4gVQzy/C3gCOC8inpT0FWBPRPzlWIscD0krImLY+xKVpFrOpVrOA3wu5ahazgNK\ndy6jfgKIiMUF9rEJ2BQRTyb7dwM3FviaZmZWoAkfBhoR24BXJHUmhy4Enpvofs3MbGSFDgO9XNIm\n4Fzg+5IeTI7PkfRAXtM/Bu6U9HPgTOCvC+l3jJaWoI9SqZZzqZbzAJ9LOaqW84ASncuo9wDMzKw6\n+ZvAZmY1ygFgZlajqi4AJF0saZ2kjZIqdrSRpPmSfiTpuWS6jU+mXVOhJNVLelrS/WnXUohqmdpE\n0p8l/7ZWS/pXSZPSrmmsJP2zpB2SVucdmyHpIUkbkt/Hp1njWA1zLl9I/n39XNIySdMnou+qCgBJ\n9cBXgUuAhcDVkhamW9W49QGfioiFwDnAxyv4XAZ9kuqYBqTipzaRNBf4E6ArIk4D6oGPplvVMbkd\nuPioYzcCyyOiHVhO5Qw3v51fPZeHgNMi4nRgPfDpiei4qgIAWARsjIjnkyknvgksSbmmcYmIrRHx\nVLK9l9ybzNx0qxo/SfOADwC3pl1LIfKmNrkNclObRMRr6VY1bg3AZEkNwBRgS8r1jFlEPAbsPurw\nEuCOZPsO4LKSFjVOQ51LRPwwIvqS3SeAeRPRd7UFwFzglbz9TVTwm+YgSQuAs4AnR25Z1r4M3AAM\npF1IgfKnNnla0q2SpqZd1LGKiM3A3wIvA1uB1yPih+lWVbBsRGxNtrcB2TSLKaLfB34wES9cbQFQ\ndSS1AN8B/jQi9qRdz3hI+iCwIyJWpl1LETQAZwNfi4izgH1UzqWGw5Lr40vIBdocYKqk3063quKJ\n3Pj2ih/jLukmcpeD75yI16+2ANgMzM/bn5ccq0iSGsm9+d8ZEfekXU8BzgM+LOlFcpflLpD0jXRL\nGrehpjY5O8V6xmsx8EJEdEfEIeAe4N0p11So7ZJmAyS/d6RcT0EkXQt8EPitmKAvbFVbAPwMaJd0\nsqQmcje17ku5pnFRbo7t24C1EfHFtOspRER8OiLmRcQCcv9NHomIivxrs4qmNnkZOEfSlOTf2oVU\n4M3so9wHXJNsXwPcm2ItBZF0MblLph+OiP0T1U9VBUBy0+QTwIPk/jHfFRFr0q1q3M4DfofcX8vP\nJD+Xpl2UAelMbVJUySeYu4HkzLHWAAAAaElEQVSngGfJvRdUzFQKkv4VeBzolLRJ0h8ANwMXSdpA\n7hPOzWnWOFbDnMstQAZ4KPn//j9OSN+eCsLMrDZV1ScAMzMbOweAmVmNcgCYmdUoB4CZWY1yAJiZ\n1SgHgJlZjXIAmJnVqP8Pcxf3agWen3AAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r_pCt_yOsmEz",
        "colab_type": "code",
        "outputId": "4cbf74f8-4240-4acc-d2b1-d232d771a2a5",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "source": [
        "barr_svm = SVMClassifier(barr.f.w, K_hat, X, y)\n",
        "print('Accuracy of SVM using Barrier method=', barr_svm.accuracy(X, y))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Accuracy of SVM using Barrier method= 0.5561993047508691\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3NMIfgxCFnud",
        "colab_type": "text"
      },
      "source": [
        "# PDIP (Primal - Dual Interior Point method)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TagCxWk8JyZq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class PDIP:\n",
        "  \n",
        "  def __init__(self, w_0, u_0, v_0, landa_0, \n",
        "               Y=y, K=K_hat, C=0.1, mu=1000, m=y.shape[0]*2, t_0=1, \n",
        "               alpha=0.01, beta=0.5, eps_feas=1e-6, eps=2e-6):\n",
        "    \n",
        "    self.w = w_0\n",
        "    self.u = u_0\n",
        "    self.v = v_0\n",
        "    self.landa = landa_0\n",
        "    self.z = np.hstack((self.w, self.u, self.v, self.landa))\n",
        "    self.dz = self.z.copy()\n",
        "    \n",
        "    \n",
        "    self.y = Y\n",
        "    self.K = K\n",
        "    self.C = C\n",
        "    self.mu = mu\n",
        "    self.m = m\n",
        "    self.t = t_0\n",
        "\n",
        "    self.alpha = alpha\n",
        "    self.beta = beta\n",
        "    self.eps = eps\n",
        "    self.eps_feas = eps_feas\n",
        "    \n",
        "    self.s = 0\n",
        "    \n",
        "    self.iteration = 0\n",
        "    self.all_obj = []\n",
        "    \n",
        "    self.update()\n",
        "  \n",
        "  def update(self):\n",
        "\n",
        "    self.eta = self._duality_gap()\n",
        "\n",
        "    self.r_dual, self.r_cent, self.r_prim = self.r(self.z)\n",
        "    \n",
        "    self.z = self.z + self.s*self.dz\n",
        "    \n",
        "    self.w, self.u, self.v, self.landa = self.unpack(self.z)\n",
        "    self.dw, self.du, self.dv, self.dlanda = self.unpack(self.dz)\n",
        "    \n",
        "    self.iteration += 1\n",
        "    \n",
        "    self.obj = self.objective()\n",
        "    self.all_obj.append(self.obj)\n",
        "    \n",
        "  def _duality_gap(self):\n",
        "    return np.dot(self.w, self.u) + np.dot((self.C - self.w), self.v)\n",
        "  \n",
        "  def unpack(self, z):\n",
        "    s = self.w.shape[0] # size\n",
        "    w = z[:s]\n",
        "    u = z[s:2*s]\n",
        "    v = z[2*s:3*s]\n",
        "    landa = z[-1]\n",
        "    return w, u, v, landa\n",
        "  \n",
        "  def r(self, z):\n",
        "    w, u, v, landa = self.unpack(z) \n",
        "    r_dual = self.K @ w -1 - u + v + self.y*landa\n",
        "    r_cent = np.diag(np.hstack((u, v))) @ np.hstack(\n",
        "        (w, (self.C - w))) - 1/self.t\n",
        "    r_prim = np.dot(self.y, w)\n",
        "    return r_dual, r_cent, r_prim\n",
        "  \n",
        "  def r2(self, z):\n",
        "    return(np.linalg.norm(np.hstack(self.r(z))))\n",
        "      \n",
        "  def pd_direction(self):\n",
        "    y = self.y.reshape(-1, 1)\n",
        "    n = self.w.shape[0]\n",
        "    I = np.eye(n)\n",
        "    Z = np.zeros((n, n))\n",
        "    z = np.zeros((n, 1))\n",
        "    \n",
        "    drz = np.block([[self.K, -I, I, y], \n",
        "                    [np.diag(self.u), np.diag(self.w), Z, z],\n",
        "                    [-np.diag(self.v), Z, np.diag(self.C - self.w), z],\n",
        "                    [y.T, z.T, z.T, 0]])\n",
        "    \n",
        "    rz = np.hstack((self.r_dual, self.r_cent, self.r_prim))\n",
        "    \n",
        "    \n",
        "    self.dz = -np.linalg.inv(drz) @ rz\n",
        "  \n",
        "  def _s_max(self):    \n",
        "    L = np.hstack((self.u, self.v))\n",
        "    dL = np.hstack((self.du, self.dv))\n",
        "    Ls = [-L[i] / dL[i] for i in range(self.m) if dL[i] < 0]\n",
        "    if not Ls:\n",
        "      return 1\n",
        "    return np.minimum(1, min(Ls))\n",
        "  \n",
        "  def line_search(self):\n",
        "    s_max = self._s_max()\n",
        "    self.s = 0.99 * s_max\n",
        "    \n",
        "    while (self.r2(self.z+self.dz) > \n",
        "        (1 - self.alpha*self.s) * self.r2(self.z)):\n",
        "      self.s = self.s * self.beta\n",
        "\n",
        "  \n",
        "  def objective(self):\n",
        "    return 0.5 * self.w.T @ self.K @ self.w - np.sum(self.w)\n",
        "    \n",
        "  def report(self):\n",
        "    print('Iter {}: f*={}'.format(self.iteration, self.obj))\n",
        "  \n",
        "  def run(self):\n",
        "    while True:\n",
        "      self.t = self.mu * self.m / self.eta\n",
        "      self.pd_direction()\n",
        "      self.line_search()\n",
        "      self.update()\n",
        "      self.report()\n",
        "    \n",
        "      if (np.sqrt(np.linalg.norm(self.r_prim)**2 + \n",
        "                  np.linalg.norm(self.r_dual)**2) <= self.eps_feas \n",
        "          and self.eta <= self.eps):\n",
        "        break\n",
        "        "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1TMRXRBK0Fnv",
        "colab_type": "text"
      },
      "source": [
        "### Run PDIP"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uZyVhE5-49iV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# PDIP initialization params\n",
        "w_0 = (y==-1) * np.count_nonzero(y==1) + (y==1) * np.count_nonzero(y==-1)\n",
        "w_0 = w_0 / 10000\n",
        "# w_0, u_0, v_0 are all strictly feasible\n",
        "u_0 = np.random.rand(w_0.shape[0])\n",
        "v_0 = np.random.rand(w_0.shape[0])\n",
        "landa_0 = 1"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pMjy6TGMBg6P",
        "colab_type": "code",
        "outputId": "e3b4c86d-4ebb-4518-d68d-66ff0c054c76",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        }
      },
      "source": [
        "pdip = PDIP(w_0, u_0, v_0, landa_0)\n",
        "pdip.run()"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Iter 1: f*=-37.755232350470536\n",
            "Iter 2: f*=-37.88678259812138\n",
            "Iter 3: f*=-37.94213445448706\n",
            "Iter 4: f*=-37.9978006201082\n",
            "Iter 5: f*=-38.062201236826695\n",
            "Iter 6: f*=-38.13625441300278\n",
            "Iter 7: f*=-38.21997794948979\n",
            "Iter 8: f*=-38.313220170589084\n",
            "Iter 9: f*=-38.415713859714046\n",
            "Iter 10: f*=-38.52705888695394\n",
            "Iter 11: f*=-38.64670519541937\n",
            "Iter 12: f*=-38.77394521413474\n",
            "Iter 13: f*=-38.907917408451326\n",
            "Iter 14: f*=-39.04762128908448\n",
            "Iter 15: f*=-39.19194308092098\n",
            "Iter 16: f*=-39.339689981850626\n",
            "Iter 17: f*=-39.48962980786405\n",
            "Iter 18: f*=-39.64053217296838\n",
            "Iter 19: f*=-39.79120738408049\n",
            "Iter 20: f*=-39.9405398747745\n",
            "Iter 21: f*=-40.08751394071307\n",
            "Iter 22: f*=-40.231230198571765\n",
            "Iter 23: f*=-40.37090999885564\n",
            "Iter 24: f*=-40.505871064751005\n",
            "Iter 25: f*=-40.635138775332436\n",
            "Iter 26: f*=-40.76011587100263\n",
            "Iter 27: f*=-40.87950157216804\n",
            "Iter 28: f*=-40.99153195587384\n",
            "Iter 29: f*=-41.0986636629861\n",
            "Iter 30: f*=-41.19929555091228\n",
            "Iter 31: f*=-41.29420300850996\n",
            "Iter 32: f*=-41.38466410973582\n",
            "Iter 33: f*=-41.46880359139838\n",
            "Iter 34: f*=-41.54738842589025\n",
            "Iter 35: f*=-41.62063573851381\n",
            "Iter 36: f*=-41.68876819257679\n",
            "Iter 37: f*=-41.75202427244648\n",
            "Iter 38: f*=-41.81065442106929\n",
            "Iter 39: f*=-41.86490966106038\n",
            "Iter 40: f*=-41.915035620440655\n",
            "Iter 41: f*=-41.96127138666931\n",
            "Iter 42: f*=-42.00384908677297\n",
            "Iter 43: f*=-42.042993931477426\n",
            "Iter 44: f*=-42.07892479741624\n",
            "Iter 45: f*=-42.111854584210874\n",
            "Iter 46: f*=-42.14198978472196\n",
            "Iter 47: f*=-42.169529412421525\n",
            "Iter 48: f*=-42.19466266149473\n",
            "Iter 49: f*=-42.2175896739513\n",
            "Iter 50: f*=-42.23845643167307\n",
            "Iter 51: f*=-42.25743536250456\n",
            "Iter 52: f*=-42.27439295096762\n",
            "Iter 53: f*=-42.29005265657916\n",
            "Iter 54: f*=-42.30427825355648\n",
            "Iter 55: f*=-42.3171703480118\n",
            "Iter 56: f*=-42.328841903387364\n",
            "Iter 57: f*=-42.33939961046147\n",
            "Iter 58: f*=-42.34894227374471\n",
            "Iter 59: f*=-42.35756101505492\n",
            "Iter 60: f*=-42.36533967310113\n",
            "Iter 61: f*=-42.372355230516476\n",
            "Iter 62: f*=-42.378678251924896\n",
            "Iter 63: f*=-42.38437332699289\n",
            "Iter 64: f*=-42.38949950903056\n",
            "Iter 65: f*=-42.39411073806656\n",
            "Iter 66: f*=-42.3982562388801\n",
            "Iter 67: f*=-42.40198088846409\n",
            "Iter 68: f*=-42.40532555266144\n",
            "Iter 69: f*=-42.40832739659252\n",
            "Iter 70: f*=-42.41102017593987\n",
            "Iter 71: f*=-42.41343451490509\n",
            "Iter 72: f*=-42.41559817247188\n",
            "Iter 73: f*=-42.41753629421184\n",
            "Iter 74: f*=-42.41927164504404\n",
            "Iter 75: f*=-42.420824819874795\n",
            "Iter 76: f*=-42.42221443240124\n",
            "Iter 77: f*=-42.42345728528397\n",
            "Iter 78: f*=-42.42456852599615\n",
            "Iter 79: f*=-42.42556179195804\n",
            "Iter 80: f*=-42.4264493469976\n",
            "Iter 81: f*=-42.42724220971582\n",
            "Iter 82: f*=-42.42795027347174\n",
            "Iter 83: f*=-42.42858241745786\n",
            "Iter 84: f*=-42.42914660851107\n",
            "Iter 85: f*=-42.429649993665734\n",
            "Iter 86: f*=-42.43009898381615\n",
            "Iter 87: f*=-42.430499329093976\n",
            "Iter 88: f*=-42.43085618665283\n",
            "Iter 89: f*=-42.43117418152578\n",
            "Iter 90: f*=-42.43145746114027\n",
            "Iter 91: f*=-42.43170974400442\n",
            "Iter 92: f*=-42.43193436307209\n",
            "Iter 93: f*=-42.432134304357504\n",
            "Iter 94: f*=-42.43231224144468\n",
            "Iter 95: f*=-42.432470566530384\n",
            "Iter 96: f*=-42.43261141852267\n",
            "Iter 97: f*=-42.43273670854907\n",
            "Iter 98: f*=-42.43284814310058\n",
            "Iter 99: f*=-42.43294724498887\n",
            "Iter 100: f*=-42.43303537229512\n",
            "Iter 101: f*=-42.43311373549812\n",
            "Iter 102: f*=-42.43318341295918\n",
            "Iter 103: f*=-42.433245364916644\n",
            "Iter 104: f*=-42.433300446116185\n",
            "Iter 105: f*=-42.43334941718457\n",
            "Iter 106: f*=-42.43339295484698\n",
            "Iter 107: f*=-42.43343166108341\n",
            "Iter 108: f*=-42.43346607131839\n",
            "Iter 109: f*=-42.43349666173573\n",
            "Iter 110: f*=-42.433523855804694\n",
            "Iter 111: f*=-42.43354803009825\n",
            "Iter 112: f*=-42.43356951947725\n",
            "Iter 113: f*=-42.43358862170711\n",
            "Iter 114: f*=-42.43360560156674\n",
            "Iter 115: f*=-42.43362069450431\n",
            "Iter 116: f*=-42.43363410988769\n",
            "Iter 117: f*=-42.43364603389419\n",
            "Iter 118: f*=-42.43365663207888\n",
            "Iter 119: f*=-42.43366605165693\n",
            "Iter 120: f*=-42.43367442353254\n",
            "Iter 121: f*=-42.43368186410249\n",
            "Iter 122: f*=-42.43368847685994\n",
            "Iter 123: f*=-42.43369435382077\n",
            "Iter 124: f*=-42.43369957679255\n",
            "Iter 125: f*=-42.43370421850359\n",
            "Iter 126: f*=-42.43370834360728\n",
            "Iter 127: f*=-42.43371200957606\n",
            "Iter 128: f*=-42.433715267496666\n",
            "Iter 129: f*=-42.433718162777836\n",
            "Iter 130: f*=-42.43372073577967\n",
            "Iter 131: f*=-42.43372302237396\n",
            "Iter 132: f*=-42.433725054442036\n",
            "Iter 133: f*=-42.433726860318174\n",
            "Iter 134: f*=-42.43372846518336\n",
            "Iter 135: f*=-42.433729891415965\n",
            "Iter 136: f*=-42.43373115890326\n",
            "Iter 137: f*=-42.43373228531843\n",
            "Iter 138: f*=-42.43373328636699\n",
            "Iter 139: f*=-42.4337341760056\n",
            "Iter 140: f*=-42.433734966637005\n",
            "Iter 141: f*=-42.43373566928268\n",
            "Iter 142: f*=-42.43373629373693\n",
            "Iter 143: f*=-42.4337368487034\n",
            "Iter 144: f*=-42.43373734191662\n",
            "Iter 145: f*=-42.43373778024981\n",
            "Iter 146: f*=-42.433738169810994\n",
            "Iter 147: f*=-42.43373851602817\n",
            "Iter 148: f*=-42.43373882372496\n",
            "Iter 149: f*=-42.43373909718805\n",
            "Iter 150: f*=-42.43373934022692\n",
            "Iter 151: f*=-42.433739556226996\n",
            "Iter 152: f*=-42.43373974819692\n",
            "Iter 153: f*=-42.433739918810446\n",
            "Iter 154: f*=-42.433740070443754\n",
            "Iter 155: f*=-42.43374020520859\n",
            "Iter 156: f*=-42.433740324981684\n",
            "Iter 157: f*=-42.43374043143098\n",
            "Iter 158: f*=-42.43374052603873\n",
            "Iter 159: f*=-42.433740610122335\n",
            "Iter 160: f*=-42.43374068485258\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kSxaGtqi_Jot",
        "colab_type": "code",
        "outputId": "1e6e3285-5c67-41d9-8e18-eaf857780017",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 269
        }
      },
      "source": [
        "plt.plot(pdip.all_obj);"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHONJREFUeJzt3Xl8FfW9//HXZ87JAmFfVLYIKq4g\noHG3VZS6VUV/2mvVSt1K9Wpvb+ttq+VRvfehXrX4uy5VW2nVa1u6qlWLSxWXWq0bUAQUUERBFktA\nDRCWkORz/zgTPOJZApNkksn7+XicB2e+czLzyZC8z+R7vvMdc3dERCT5grgLEBGRtqHAFxHpJBT4\nIiKdhAJfRKSTUOCLiHQSCnwRkU5CgS8i0kko8EVEOgkFvohIJ5GOu4Bs/fr186FDh8ZdhohIhzJz\n5szV7t6/2OvaVeAPHTqUGTNmxF2GiEiHYmZLmvM6demIiHQSCnwRkU4iUuCb2WQzW2Bmc8zsT2bW\nK2wvMbP7zWyumc03s6taplwREdlRUc/wnwZGuPv+wNtAU7B/BShz95HAgcA3zWxoxH2JiEgEkQLf\n3Z9y9/pw8RVgcNMqoMLM0kAXoA5YG2VfIiISTUv24V8IPBE+fwCoBVYCS4Gb3f2jFtyXiIhsp6LD\nMs1sOrBLjlWT3P2R8DWTgHpgarjuYKABGAj0Bv5mZtPdfXGO7U8EJgJUVlbuyPcgIiLNUDTw3X1c\nofVmdj5wMnCsf3q/xHOAJ919C7DKzF4CqoDPBb67TwGmAFRVVe3Q/RYXfriOaXNWcP7hQ+nbrWxH\nNiEiknhRR+mcAHwfONXdN2StWgocE76mAjgUWBBlX4UsWrWenzy7iDW1da21CxGRDi9qH/4dQHfg\naTObbWY/C9vvBLqZ2ZvA68B97j4n4r7ySgUGQH2DbsguIpJPpKkV3H2PPO3ryQzNbBNNgd/QqMAX\nEcknEVfappsC3xX4IiL5JCLwPz3Db4y5EhGR9itRga8+fBGR/BIV+OrSERHJLxGBn9aHtiIiRSUi\n8IOmLh0FvohIXokI/KYz/EYFvohIXokI/JTO8EVEikpU4KsPX0Qkv0QEvj60FREpLhGBnwoy34YC\nX0Qkv2QEvqkPX0SkmGQEfkqjdEREiklE4Kc1SkdEpKhEBH5gmjxNRKSYRAS+RumIiBSXiMBv6sNX\nl46ISH7JCHzTGb6ISDHJCHxNjywiUlSkwDeza81sTngD86fMbGDYbmZ2u5ktCtcf0DLl5ra1D183\nQBERySvqGf5kd9/f3UcD04Crw/YTgeHhYyLw04j7KUiTp4mIFBcp8N19bdZiBdCUuOOBX3rGK0Av\nMxsQZV+FmBmBQaO6dERE8kpH3YCZXQ9MAGqAsWHzIOCDrJctC9tW5vj6iWT+CqCysnKH60gHgc7w\nRUQKKHqGb2bTzWxejsd4AHef5O5DgKnA5dtbgLtPcfcqd6/q37//9n8HoSDQKB0RkUKKnuG7+7hm\nbmsq8DhwDbAcGJK1bnDY1mrSQUC9PrQVEckr6iid4VmL44EF4fNHgQnhaJ1DgRp3/1x3TktKBaY+\nfBGRAqL24d9oZnsBjcAS4JKw/XHgJGARsAG4IOJ+ikoFRr3m0hERyStS4Lv7GXnaHbgsyra3Vyow\n9eGLiBSQiCttIXPxlQJfRCS/xAR+YKZhmSIiBSQm8NMpneGLiBSSmMBXH76ISGHJCXxT4IuIFJKc\nwA/Uhy8iUkhiAj+dMhoV+CIieSUm8FMapSMiUlByAl8f2oqIFJSYwE8HgQJfRKSAxAS+pkcWESks\nMYGfuQGKJk8TEcknMYGfCgxNhy8ikl+yAl9n+CIieSUq8HXHKxGR/BIT+Gnd8UpEpKDEBH6gqRVE\nRApKTODrBigiIoUlJvB1pa2ISGGRAt/MrjWzOWY228yeMrOBYfu5YftcM/u7mY1qmXLz0/TIIiKF\nRT3Dn+zu+7v7aGAacHXY/h5wlLuPBK4FpkTcT1HplPrwRUQKSUf5Yndfm7VYAXjY/ves9leAwVH2\n0xypQNMji4gUEinwAczsemACUAOMzfGSi4AnCnz9RGAiQGVl5Q7XoemRRUQKK9qlY2bTzWxejsd4\nAHef5O5DgKnA5dt87Vgygf+DfNt39ynuXuXuVf3799/hbySl2TJFRAoqeobv7uOaua2pwOPANQBm\ntj/wC+BEd1+zwxU2UzqlD21FRAqJOkpneNbieGBB2F4JPASc5+5vR9lHc2lYpohIYVH78G80s72A\nRmAJcEnYfjXQF7jLzADq3b0q4r4KyvTha/I0EZF8oo7SOSNP+8XAxVG2vb1SgdHo4O6EbzIiIpIl\nMVfapoNMyKtbR0Qkt8QEfhAGvoZmiojklpjAbzrD1xTJIiK5JSbwUzrDFxEpKHGB36C7XomI5JSY\nwN/6oa26dEREckpM4KeCzLeiUToiIrklKPAz/6oPX0QktwQFfuZb0RTJIiK5JSbw0xqlIyJSUGIC\nP9h6pa3m0xERySUxgf/p1AoxFyIi0k4lJvA/vfBKiS8ikktyAt80eZqISCHJCfyUAl9EpJDEBL6m\nRxYRKSwxgd/UpaNhmSIiuSUn8JumR1bgi4jkFPUm5tea2Rwzm21mT5nZwG3WH2Rm9WZ2ZrQyi0un\ndIYvIlJI1DP8ye6+v7uPBqaRuXk5AGaWAm4Cnoq4j2YJNEpHRKSgSIHv7muzFiuA7LT9FvAgsCrK\nPporrdkyRUQKSkfdgJldD0wAaoCxYdsg4PRw+aCo+2gO3fFKRKSwomf4ZjbdzObleIwHcPdJ7j4E\nmApcHn7ZrcAP3L3oZa9mNtHMZpjZjOrq6h3+RlIalikiUlDRM3x3H9fMbU0FHgeuAaqA31mmX70f\ncJKZ1bv7wzm2PwWYAlBVVbXDaZ3SHa9ERAqK1KVjZsPd/Z1wcTywAMDdh2W95n+BabnCviWlNVum\niEhBUfvwbzSzvYBGYAlwSfSSdszWPnzdxFxEJKdIge/uZzTjNedH2Udzbb3wSl06IiI5JeZKW93x\nSkSksMQEfqBROiIiBSUm8DVbpohIYYkJfI3DFxEpLHGBrz58EZHcEhf4OsMXEcktMYGvydNERApL\nTOCHJ/jq0hERySMxgW9mpALTHa9ERPJITOBDph9fZ/giIrklK/DNNHmaiEgeiQr8dGA0KO9FRHJK\nVOCnUjrDFxHJJ1mBb+rDFxHJJ1mBH5jG4YuI5JGowE8r8EVE8kpU4AcKfBGRvBIV+GmNwxcRyStR\ngZ8KjAbd4lBEJKdIgW9m15rZHDObbWZPmdnArHVHh+1vmtlfo5daXCowGnQTcxGRnKKe4U929/3d\nfTQwDbgawMx6AXcBp7r7fsBXIu6nWVJBoC4dEZE8IgW+u6/NWqwAmtL2HOAhd18avm5VlP00Vzow\nGtWlIyKSUzrqBszsemACUAOMDZv3BErM7HmgO3Cbu/8y6r6KCfShrYhIXkXP8M1supnNy/EYD+Du\nk9x9CDAVuDz8sjRwIPBl4HjgR2a2Z57tTzSzGWY2o7q6OtI3kxmHr6kVRERyKXqG7+7jmrmtqcDj\nwDXAMmCNu9cCtWb2AjAKeDvH9qcAUwCqqqoinZ7rSlsRkfyijtIZnrU4HlgQPn8EONLM0mbWFTgE\nmB9lX82RmR5ZgS8ikkvUPvwbzWwvoBFYAlwC4O7zzexJYE647hfuPi/ivopKp4xN9Qp8EZFcIgW+\nu59RYN1kYHKU7W8v3eJQRCS/ZF1pq+mRRUTySlbg60NbEZG8EhX46ZQCX0Qkn0QFfqBROiIieSUq\n8DU9sohIfokK/FQQ6AxfRCSPhAU+CnwRkTwSFviaHllEJJ9EBb6mRxYRyS9RgZ8KjPoGzZYpIpJL\n4gJfffgiIrklKvDTuom5iEheiQr8QGf4IiJ5JSrwdeGViEh+iQr8VGC4oymSRURySFbgmwGoH19E\nJIdkBX4qDHyd4YuIfE6iAj8dKPBFRPJJVOAHYZdOzcYtMVciItL+RAp8M7vWzOaY2Wwze8rMBobt\nPc3sz2b2hpm9aWYXtEy5hR0yrC+l6YAL7nud6nWb22KXIiIdRtQz/Mnuvr+7jwamAVeH7ZcBb7n7\nKOBo4P+bWWnEfRU1cnBP7jv/IJZ+tIHxd7zI3X99lzXrFfwiIhAx8N19bdZiBdDUee5AdzMzoBvw\nEVAfZV/NdcQe/Zj6jUMY3LsrNzyxgMNueJZ/++0/mLnko7bYvYhIu2UecQijmV0PTABqgLHuXm1m\n3YFHgb2B7sBZ7v5Ynq+fCEwEqKysPHDJkiWR6sn29j/X8ZtXl/LgrGWs21TPQUN7861jhvOF4f2w\nsL9fRKSjM7OZ7l5V9HXFAt/MpgO75Fg1yd0fyXrdVUC5u19jZmcCRwDfBXYHngZGbfMXwedUVVX5\njBkzitW83TbU1fP71z/g5y8sZkXNJo7Yoy8/PGkf9hvYs8X3JSLS1los8Ldjh5XA4+4+wsweA250\n97+F654FrnT31wpto7UCv8nm+gamvrKUO55bxCcb6phw2FCuOG5PupeXtNo+RURaW3MDP+ooneFZ\ni+OBBeHzpcCx4Wt2BvYCFkfZV0soS6e48MhhPHfF0ZxzSCX3v/w+J93+N2Yu+Tju0kREWl3UUTo3\nmtk8M5sDHAd8O2y/FjjczOYCzwA/cPfVEffVYnp2LeG600byx28eRmMj/MvdL3Pb9Hd08xQRSbQW\n69JpCa3dpZPL2k1b+NHD83hk9gqqdu3NHeccwC49y9u0BhGRKNqkSycJepSXcNtXx3DrWaOZv3It\n4+98kbnLauIuS0SkxXX6wG9y2phBPHDp4aSDgH+5+2WenLcy7pJERFqUAj/LPgN68PBlR7D3gO5c\n8utZ3PX8ItpTl5eISBQK/G30717Gb79xKKeOGsiPn1zIdY/NV+iLSCKk4y6gPSovSXHrWaPpU1HK\nPS++x4a6Bq47bQSpQFfnikjHpcDPIwiMa07Zl4qyFHc+9y4b6+q5+SujSKf0R5GIdEwK/ALMjO8d\nvzddS9NM/stCNm1p5CfnjKFEoS8iHZCSqxkuG7sHV5+8L0+++SHf+f1s3VFLRDokneE304VHDmNL\nQyM3PLGA0nTAzWeOIlCfvoh0IAr87fDNo3Znc30j//P025SXpLj+tBGaZllEOgwF/nb61jF7sGlL\nA3c9/y7l6RQ/Onkfhb6IdAgK/O2U+SB3LzZuaeDel96jZ5cSvj1uePEvFBGJmQJ/B5gZP/ryvqzd\nWM8t09+mZ5c05x8xLO6yREQKUuDvoCAwbjpjJOs2beE///wWPbuWcPqYwXGXJSKSl4ZlRpBOBdx+\n9hgO370v//HHOUx/659xlyQikpcCP6LykhRTJlQxYmAP/vU3s3j53TVxlyQikpMCvwV0K0vzvxcc\nzK59uvKNX85gzrJP4i5JRORzFPgtpHdFKb+66BB6dS3h6/e+xqJV6+IuSUTkM1ok8M3sCjNzM+sX\nLpuZ3W5mi8xsjpkd0BL7ae926VnOry86hFQQcN49r7Hs4w1xlyQislXkwDezIWRuYL40q/lEYHj4\nmAj8NOp+Ooqh/Sr41UUHU7u5nvPueY3qdZvjLklEBGiZM/xbgO8D2TOKjQd+6RmvAL3MbEAL7KtD\n2GdAD+674CBW1mzk6/e+Rs3GLXGXJCISLfDNbDyw3N3f2GbVIOCDrOVlYVunceCufbj7vCreWbWO\ni+9/nY11DXGXJCKdXNHAN7PpZjYvx2M88EPg6igFmNlEM5thZjOqq6ujbKrdOWrP/tx61hhmLPmY\nS6fOpK6+Me6SRKQTKxr47j7O3Uds+wAWA8OAN8zsfWAwMMvMdgGWA0OyNjM4bMu1/SnuXuXuVf37\n94/6/bQ7X95/AP99+kieX1jNFX98Q3Ppi0hsdnhqBXefC+zUtByGfpW7rzazR4HLzex3wCFAjbuv\njFpsR3X2wZXUbNzCjU8soEd5mus0rbKIxKC15tJ5HDgJWARsAC5opf10GJcctTs1G7fw0+ffpVfX\nEr53/N5xlyQinUyLBb67D8167sBlLbXtpPj+8XvxyYYt3Pncu/QoL+GbR+0ed0ki0olotsw2ZGZc\nd9oI1m3awg1PLKDR4dKjFfoi0jYU+G0sFRi3nDWawIybnlzApi0N/Pu44erTF5FWp8CPQUkq4Jaz\nRlOWDrjtmXfYVN/AlSfsrdAXkValwI9JKjBuOmN/ykoC7v7rYjZsbuA/T92PVKDQF5HWocCPURAY\n144fQUVpmrtfWMzq9Zu55azRlJek4i5NRBJIgR8zM+Oqk/ahf/cyrntsPmvWv8bPJ1TRs2tJ3KWJ\nSMJoPvx24uIv7MZPzh7D7A8+4cyf/V1TK4tIi1PgtyOnjBrI/RcezIdrNzH+jpd0u0QRaVEK/Hbm\nsN378vBlR9Crawlfu+dV7nvpPTLXsYmIRKPAb4d279+Nhy87gmP23on/+vNbXPGHNzS9sohEpsBv\np7qXl3D31w7ku1/akz/NXs4pd7zImytq4i5LRDowBX47FgTGvx07nF9deAjrNm3htDtf4ucvLKZR\nUyyLyA5Q4HcARw7vx5Pf/iJj99qJ6x+fz4R7X+ODjzSKR0S2jwK/g+hdUcrd5x3IDf9vJLOWfsxx\nt7zAL/62WDdUEZFmU+B3IGbG2QdX8tR3vsihu/Xhusfmc/pdL6lvX0SaRYHfAQ3u3ZV7zz+I288e\nw4pPNnLKT17kygfnsGrtprhLE5F2TIHfQZkZp44ayPTvHsX5hw/jwVnLOPrm57nl6bep3Vwfd3ki\n0g5Ze7qop6qqymfMmBF3GR3S+6tr+fFfFvD43A/pU1HKRUcO47zDdqVHuebkEUk6M5vp7lVFX6fA\nT5ZZSz/m9mfe4fmF1XQvSzPh8F05//Bh9O9eFndpItJK2jTwzewK4Gagv7uvNrNzgR8ABqwDLnX3\nN4ptR4HfcuYtr+Gu5xfxxLwPSQfG8fvtwtcO3ZVDhvXRjVZEEqa5gR95emQzGwIcByzNan4POMrd\nPzazE4EpwCFR9yXNN2JQT+4690AWV6/n168s5YGZHzBtzkqG79SNMw8czCmjBjKwV5e4yxSRNhT5\nDN/MHgCuBR4Bqtx99TbrewPz3H1QsW3pDL/1bKxr4M9zVvCbV5cy+4NPADh4WB9OHTWQL+27Mzv3\nKI+5QhHZUW3SpWNm44Fj3P3bZvY+uQP/P4C93f3iYttT4LeNJWtqeXT2Ch6evZx3q2sBGDGoB8fu\nvTNj996JEQN7kE5pAJdIR9FigW9m04FdcqyaBPwQOM7da3IFvpmNBe4CjnT3nJO7m9lEYCJAZWXl\ngUuWLClWs7QQd+edVeuZPv+fPDt/FbOWfkyjQ7eyNAcN7c2hu/Wlamgf9h3Qgy6luu2iSHvV6mf4\nZjYSeAZomtRlMLACONjdPzSz/YE/ASe6+9vN2abO8OP1UW0dLy1azSuL1/Dy4jUsDs/+U4ExfKdu\njBzUk5GDezJyUE+G79ydbmW6Q6ZIe9DmwzKzz/DNrBJ4Fpjg7n9v7jYU+O3LqrWb+McHnzB3WQ1z\nl2ceH9XWbV2/c48yduvXjd13qmC3ft0Y1r+CQb26MLBXF70ZiLShNhulk8fVQF/grnAIYH1zipH2\nZace5Ry/3y4cv1+mR8/dWVGzibnLani3ej3vVq9ncXXm84C1mz57dW/38jQDe3ZhQK9ydulRTp+K\nUvpUlNK3Wym9u5bSt6KMPt1K6VtRSnmJuotE2kKLBb67D816fjFQ9ENa6VjMjEG9ujBom+Gc7s6a\n2jreX13L8k82srJmEys/2ciKmk2srNnImyvW8nFtHfV5ZvYsSRkVZWm6hY+K8NGtLEVFaZpu5Wkq\nStOUpgPK0kH4b+ozy03Py8J1qcA+fZh9djlsCwIjnaNNJKn0d7dEZmb061ZGv25l5Pszzt1Zu7Ge\njzbU8VHtZtasr+Oj2jrW1NaxfnM96zfVU7u5nvWb66mtq6dmQx3LP66ndnMDtWFbW80EnQ4ywR8Y\nGIZZ5grCpu/VAMI2s0/XN62z8AWftn92O9kXvpl9fn1ULXFhXYu87bXARtrL8YAWOiYFnHXQEC7+\nwm6tug8FvrQJM6Nn1xJ6di1hWL+KHdpGfUMjm+sbqav/9N+6hgY2bWmkrqGRzVv/baCh0Wlwz/yb\n/XCnsdGpz9PWGC7XNzrumTeqpo+5HDJtZLW5f669aZmm5W3WkfV6tr7eaYn3s5b4SK5l6oi+lRZ5\nf2+hk4SW+d8prF+31p/+RIEvHUY6FZBOBVRoWiCRHaKra0REOgkFvohIJ6HAFxHpJBT4IiKdhAJf\nRKSTUOCLiHQSCnwRkU5CgS8i0km0q5uYm1k1sKMT4vcDVhd9Vdtrr3VB+61NdW0f1bV9kljXru7e\nv9iL2lXgR2FmM9rjjJzttS5ov7Wpru2jurZPZ65LXToiIp2EAl9EpJNIUuBPibuAPNprXdB+a1Nd\n20d1bZ9OW1di+vBFRKSwJJ3hi4hIAYkIfDM7wcwWmtkiM7syxjqGmNlzZvaWmb1pZt8O2/uY2dNm\n9k74b++Y6kuZ2T/MbFq4PMzMXg2P2+/NrDSGmnqZ2QNmtsDM5pvZYe3heJnZd8L/w3lm9lszK4/j\neJnZvWa2yszmZbXlPD6WcXtY3xwzO6CN65oc/j/OMbM/mVmvrHVXhXUtNLPj27KurHVXmJmbWb9w\nOdbjFbZ/Kzxmb5rZj7PaW+d4Ze7o03EfQAp4F9gNKAXeAPaNqZYBwAHh8+7A28C+wI+BK8P2K4Gb\nYqrvu8BvgGnh8h+Ar4bPfwZcGkNN9wMXh89LgV5xHy9gEPAe0CXrOJ0fx/ECvggcAMzLast5fICT\ngCfI3I3vUODVNq7rOCAdPr8pq659w9/LMmBY+Puaaqu6wvYhwF/IXOfTr50cr7HAdKAsXN6ptY9X\nq/6wtsUDOAz4S9byVcBVcdcV1vII8CVgITAgbBsALIyhlsHAM8AxwLTwh3x11i/oZ45jG9XUMwxW\n26Y91uMVBv4HQB8yd4WbBhwf1/EChm4TFDmPD3A3cHau17VFXdusOx2YGj7/zO9kGLyHtWVdwAPA\nKOD9rMCP9XiROYEYl+N1rXa8ktCl0/TL2WRZ2BYrMxsKjAFeBXZ295Xhqg+BnWMo6Vbg+0BjuNwX\n+MTd68PlOI7bMKAauC/savqFmVUQ8/Fy9+XAzcBSYCVQA8wk/uPVJN/xaU+/CxeSOXuGmOsys/HA\ncnd/Y5tVcR+vPYEvhN2EfzWzg1q7riQEfrtjZt2AB4F/d/e12es885bdpkOjzOxkYJW7z2zL/TZD\nmsyfuT919zFALZkuiq1iOl69gfFk3pAGAhXACW1ZQ3PFcXyKMbNJQD0wtR3U0hX4IXB13LXkkCbz\nV+ShwPeAP5iZteYOkxD4y8n0zzUZHLbFwsxKyIT9VHd/KGz+p5kNCNcPAFa1cVlHAKea2fvA78h0\n69wG9DKzphvZx3HclgHL3P3VcPkBMm8AcR+vccB77l7t7luAh8gcw7iPV5N8xyf23wUzOx84GTg3\nfDOKu67dybxxvxH+/A8GZpnZLjHXBZmf/4c84zUyf333a826khD4rwPDwxEUpcBXgUfjKCR8d74H\nmO/u/5O16lHg6+Hzr5Pp228z7n6Vuw9296Fkjs+z7n4u8BxwZox1fQh8YGZ7hU3HAm8R8/Ei05Vz\nqJl1Df9Pm+qK9XhlyXd8HgUmhKNPDgVqsrp+Wp2ZnUCm2/BUd9+wTb1fNbMyMxsGDAdea4ua3H2u\nu+/k7kPDn/9lZAZWfEjMxwt4mMwHt5jZnmQGLaymNY9Xa31A0ZYPMp+2v03m0+xJMdZxJJk/r+cA\ns8PHSWT6y58B3iHzqXyfGGs8mk9H6ewW/iAtAv5IOFqgjesZDcwIj9nDQO/2cLyA/wIWAPOAX5EZ\nMdHmxwv4LZnPEbaQCauL8h0fMh/E3xn+HswFqtq4rkVk+p6bfvZ/lvX6SWFdC4ET27Kubda/z6cf\n2sZ9vEqBX4c/Y7OAY1r7eOlKWxGRTiIJXToiItIMCnwRkU5CgS8i0kko8EVEOgkFvohIJ6HAFxHp\nJBT4IiKdhAJfRKST+D9CfWj1T004tQAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_2Xc2gxMF5ye",
        "colab_type": "code",
        "outputId": "436d3a54-f214-4472-ec4a-57b678e698c1",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "source": [
        "pdip_svm = SVMClassifier(pdip.w, K_hat, X, y)\n",
        "print('Accuracy of SVM using PDIP=', pdip_svm.accuracy(X, y))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Accuracy of SVM using PDIP= 0.996523754345307\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vh91pLx-5HN6",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}