MLMI_Exercise03_SVM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this exercise, the implementation is divided in four groups labeled as TM1, TM2, TM3 and TM4. Each team member will work in only one of the groups. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this first part, we will implement logistic regression using batch gradient descent. The next steps will guide you through the implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "First, we will load the numpy and matplotlib library. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we load the dataset that we will use for training and testing our algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xs = np.load('./data/ex03.reg.train.ft.npy') #training features\n",
    "ys = np.load('./data/ex03.reg.train.lb.npy') #training labels\n",
    "\n",
    "Xt = np.load('./data/ex03.reg.test.ft.npy') #testing features\n",
    "yt = np.load('./data/ex03.reg.test.lb.npy') #testing labels\n",
    "\n",
    "thi = np.load('./data/ex03.reg.w.npy') #initial parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The features Matrix (Xs and Xt) are $m\\times n$ numpy arrays where each row correspond to a feature vector and the labels (ys and yt) are $m\\times 1$ numpy arrays where each element can be either 0 or 1. Let's print the first 5 elements, to see how they look:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*Features Matrix\n",
      "Shape: (400, 2)\n",
      "First 5 elements:\n",
      "[[0.0257384  0.07356659]\n",
      " [0.05684302 0.0531371 ]\n",
      " [0.0428545  0.06111346]\n",
      " [0.10163394 0.09047114]\n",
      " [0.03650691 0.03184156]]\n",
      "*Labels\n",
      "Shape: (400, 1)\n",
      "First 5 elements:\n",
      "[[0.]\n",
      " [0.]\n",
      " [0.]\n",
      " [1.]\n",
      " [0.]]\n",
      "*Initial weights\n",
      "(3, 1)\n",
      "[[-4.24165502]\n",
      " [26.7368333 ]\n",
      " [29.77781363]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print('*Features Matrix')\n",
    "print('Shape: ' + str(Xs.shape))\n",
    "print('First 5 elements:\\n' + str(Xs[0:5]))\n",
    "print('*Labels')\n",
    "print('Shape: ' +  str(ys.shape))\n",
    "print('First 5 elements:\\n' + str(ys[0:5]))\n",
    "print('*Initial weights')\n",
    "print(thi.shape)\n",
    "print(thi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, $m$ is the number of samples and $n$ the number of features. We have $m=400$ for the training set and $m=200$ for the testing set. The number of features is $n=2$. We can obtain this values using the code below: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "400 2\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(3, 1)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = Xs.shape[0]\n",
    "n = Xs.shape[1]\n",
    "print(m,n)\n",
    "thi.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM1$ and $\\bf TM2$. Now let's begin with the implementation. As first point, we will implement a function for evaluating the sigmoid function given a set of parameters $g(\\theta^t {\\bf x}) = g(z) = 1 / (1 + exp(-z))$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters for the function will be:\n",
    "xs.- A numpy array of dimension $m \\times n$ with the features Matrix.\n",
    "th.- Is the vector of parameters or weights of the model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function should return a $m \\times 1$ numpy array $h$ where each element $h[i, 0] = g(\\theta^t {\\bf x^{(i)}})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-10-d8f3200e3ad9>, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-10-d8f3200e3ad9>\"\u001b[0;36m, line \u001b[0;32m2\u001b[0m\n\u001b[0;31m    sig + = 1.0/(1.0 + np.exp(-(np.dot(th.T,xs))))\u001b[0m\n\u001b[0m          ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "def eval(xs, th):\n",
    "    sig + = 1.0/(1.0 + np.exp(-(np.dot(th.T,xs))))\n",
    "    return sig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the function by running the next code: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 76, 100]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.0257384  0.07356659]\n",
      " [0.05684302 0.0531371 ]\n",
      " [0.0428545  0.06111346]]\n"
     ]
    }
   ],
   "source": [
    "Xaux = Xs[0:3]\n",
    "print(Xaux)\n",
    "gzaux = eval(xs=Xaux, th=thi)\n",
    "#print(gzaux)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should get an array similar to this: <br>\n",
    "[ [ 0.0257384 &nbsp;  0.07356659]<br>\n",
    " [ 0.05684302 &nbsp; 0.0531371 ]<br>\n",
    " [ 0.0428545 &nbsp;  0.06111346] ]<br>\n",
    " Maybe with some decimal variations due to possible differences in precision. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM1$ and $\\bf TM2$. Implement the loss function. Now, the parameters are: \n",
    "xs.- A numpy array of dimension $m \\times n$ with the features Matrix.\n",
    "ys.- A numpy array of dimensions $m \\times 1$ with the labels.\n",
    "th.- The vector of parameters o weights of the model. \n",
    "lamb.- The regularization parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function should return a single float. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss(xs, ys, th, lamb=0):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check we can run the next code: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "laux = loss(xs=Xs, ys=ys, th=thi)\n",
    "print(laux)\n",
    "laux = loss(xs=Xs, ys=ys, th=thi, lamb=1.0)\n",
    "print(laux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should get values around 102.3482 and 903.1364."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement the batch ($\\bf TM1$) and stochastic ($\\bf TM2$) gradient descent algorithm. The parameters for this function are: \n",
    "xs.- A numpy array of dimension $m \\times n$ with the features Matrix.\n",
    "ys.- A numpy array of dimensions $m \\times 1$ with the labels.\n",
    "th_init.- The initial vector of parameters o weights of the model. \n",
    "lt.- A float with the learning rate.\n",
    "T.- The number of running iterations.\n",
    "lamb.- The regularization parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It returns the optimized vector of parameters (th) and a python array with the loss for each iteration.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(xs, ys, th_init, lr, T, lamb=0):\n",
    "    losses = []\n",
    "    th = np.copy(th_init)\n",
    "    for i in range(T):\n",
    "        losses.append(loss(xs, ys, th, lamb))\n",
    "        #---Compute the gradient and update the weitghs:\n",
    "\n",
    "        #-----------------------------------------------\n",
    "    losses.append(loss(xs, ys, th, lamb))\n",
    "    return th, losses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, for testing, we set up the hyperparameters and load a $n \\times 1$ numpy array with the initial weights.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = 1e-3\n",
    "T = 100000\n",
    "lamb = 0 \n",
    "thi = np.load('./data/ex03.reg.w.npy') #initial parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run the algorithm and print the number of elements correctly classified. Also, we show the final weights and a graph of the loss. Running this part could take some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "op_th, ls = train(Xs, ys, thi, lr, T, lamb)\n",
    "gz = eval(Xt, op_th) #evaluating the logistic function\n",
    "y_predict = (gz >= 0.5).astype(np.int32)\n",
    "correct = np.sum(y_predict == yt)\n",
    "print(\"%d out of %d predictions correct\" % (correct, len(y_predict)))\n",
    "plt.plot(ls)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM1$ and $\\bf TM2$ Train again, but change the learning rate (lr) to 1e-5. Compare the graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = #TODO set the learning rate\n",
    "T = 100000\n",
    "lamb = 0 \n",
    "thi = np.load('./data/ex03.reg.w.npy') #initial parameters.\n",
    "op_th, ls = train(Xs, ys, thi, lr, T, lamb)\n",
    "gz = eval(Xt, op_th) #evaluating the logistic function\n",
    "y_predict = (gz >= 0.5).astype(np.int32)\n",
    "correct = np.sum(y_predict == yt)\n",
    "print(\"%d out of %d predictions correct\" % (correct, len(y_predict)))\n",
    "plt.plot(ls)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Support Vector Machine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset for this exercise is composed by 2000 images for training and 876 for testing. Each image is a $24\\times 24$ bounding box obtained from a blood sample and can contain either background or a blood parasite know as T. cruzi. The dataset was one of the result from a Mexican project for the diagnosis of Chagas disease (CONACYT/SALUD-2009-C01-113848, contact: Dr. Hugo Ruiz rpina@correo.uady.mx). The picture bellow show some negative (top) and positive (bottom) samples from the dataset. \n",
    "<img src=\"./data/ex03.dataset.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The features are $m\\times n$ numpy arrays and the labels are $m\\times 1$ numpy arrays with $m$ the number of samples and $n$ the number of features. We have $m=2000$ for the training set and $m=876$ for the testing set. The number of features is $n=12$. We save this values for the training set in two variables. As before, we load some libraries and the dataset for training and testing. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: you'll need the python module cvxopt. You can install via \"pip install cvxopt\". (For more detailed installations options see http://cvxopt.org/install/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxopt\n",
    "from cvxopt.solvers import qp\n",
    "from cvxopt.base import matrix\n",
    "import data.svm as svm\n",
    "import numpy as np\n",
    "\n",
    "Xs = np.load('./data/ex03.chagas.ft.train.npy') #training features\n",
    "ys = np.load('./data/ex03.chagas.label.1-1.train.npy') #training labels\n",
    "\n",
    "Xt = np.load('./data/ex03.chagas.ft.test.npy') #testing features\n",
    "yt = np.load('./data/ex03.chagas.label.1-1.test.npy') #testing labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we can see some of the samples and labels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*Features Matrix\n",
      "Shape: (2000, 12)\n",
      "First 5 elements:\n",
      "[[ 0.57777097  0.54774646  0.69579929  0.0460708   0.10126769  0.08421694\n",
      "   0.33381929  0.30002618  0.48413665  0.00212252  0.01025514  0.00709249]\n",
      " [ 0.52168437  0.52414216  0.66214597  0.02717556  0.02799127  0.02858336\n",
      "   0.27215458  0.274725    0.43843728  0.00073851  0.00078351  0.00081701]\n",
      " [ 0.53195125  0.45142293  0.62625272  0.0354268   0.11623582  0.08456892\n",
      "   0.28297214  0.20378266  0.39219247  0.00125506  0.01351077  0.0071519 ]\n",
      " [ 0.56395697  0.48425926  0.68728214  0.03398967  0.11100243  0.09859347\n",
      "   0.31804747  0.23450703  0.47235673  0.0011553   0.01232154  0.00972067]\n",
      " [ 0.52086057  0.48314951  0.62939815  0.0350372   0.09734438  0.07128154\n",
      "   0.27129573  0.23343345  0.39614203  0.00122761  0.00947593  0.00508106]]\n",
      "*Labels\n",
      "Shape: (2000, 1)\n",
      "First 5 elements:\n",
      "[[ 1.]\n",
      " [-1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]]\n"
     ]
    }
   ],
   "source": [
    "print('*Features Matrix')\n",
    "print('Shape: ' + str(Xs.shape))\n",
    "print('First 5 elements:\\n' + str(Xs[0:5]))\n",
    "print('*Labels')\n",
    "print('Shape: ' +  str(ys.shape))\n",
    "print('First 5 elements:\\n' + str(ys[0:5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the linear kernel as: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_kernel(x, z, kernel_params=None):\n",
    "    return np.dot(x, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where x and z are $n-$ dimensional numpy arrays (shape=[n,]) corresponding to particular features vectors. The function returns a single float. We can compute the kernel for the 0 and 1 features as: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.57777097  0.54774646  0.69579929  0.0460708   0.10126769  0.08421694\n",
      "   0.33381929  0.30002618  0.48413665  0.00212252  0.01025514  0.00709249]\n",
      " [ 0.52168437  0.52414216  0.66214597  0.02717556  0.02799127  0.02858336\n",
      "   0.27215458  0.274725    0.43843728  0.00073851  0.00078351  0.00081701]\n",
      " [ 0.53195125  0.45142293  0.62625272  0.0354268   0.11623582  0.08456892\n",
      "   0.28297214  0.20378266  0.39219247  0.00125506  0.01351077  0.0071519 ]]\n",
      "[[ 1.5734516   1.44127971  1.3565511 ]\n",
      " [ 1.44127971  1.32942692  1.24039027]\n",
      " [ 1.3565511   1.24039027  1.1765158 ]]\n"
     ]
    }
   ],
   "source": [
    "s = Xs[0:3,:]\n",
    "print(s)\n",
    "K_x = linear_kernel(s, s.T)\n",
    "print(K_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, following a similar strategy:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM3$ Compute the polynomial kernel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def polynomial_kernel(x, z, kernel_params={'a' : 1, 'c' : 1, 'd' : 3}):\n",
    "    a = kernel_params['a']\n",
    "    c = kernel_params['c']\n",
    "    d = kernel_params['d']\n",
    "    return ((np.dot(a,np.dot(x,z)) + c)**d)\n",
    "    #TODO Compute the polynomial kernel with the given parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the values for some pair of samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14.5496525513\n",
      "13.1720336949\n"
     ]
    }
   ],
   "source": [
    "K_x = polynomial_kernel(Xs[0], Xs[1])\n",
    "print(K_x)\n",
    "K_x = polynomial_kernel(Xs[3], Xs[20])\n",
    "print(K_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this, you should've obtain values around 14.5496 and 13.1720"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM4$ Compute the radial basis function (gaussian) kernel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from sklearn.svm import SVC\n",
    "def gaussian_kernel(x, z, kernel_params={'sigma' : 5.0}):\n",
    "    sigma = kernel_params['sigma']\n",
    "    return math.exp(-(np.linalg.norm(matrix(x)-matrix(z))**2)/(2*sigma**2))\n",
    "    #TODO Compute the gaussian kernel with the given parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we check the function: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9995937004876707\n",
      "0.9998377268103912\n"
     ]
    }
   ],
   "source": [
    "K_x = gaussian_kernel(Xs[0], Xs[1])\n",
    "print(K_x)\n",
    "K_x = gaussian_kernel(Xs[3], Xs[20])\n",
    "print(K_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we get values arround: 0.9996 and 0.9998"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this functions, the parameter \"kernel_params\" is a dictionary with the parameters for the kernel. \n",
    "In order to complete the svm code, we need to compute the kernel matrix. Also, once obtained the optimal Lagrange multipliers, we need to obtain the support vectors. (Note: if needed you can compute the norm of a vector with np.linalg.norm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM3$ and ${\\bf TM4}$ Given the features matrix (Xs) and a particular kernel, compute the kernel matrix (Gram matrix).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gram_matrix(X, kernel, kernel_params={'a' : 1, 'c' : 1, 'd' : 3, 'sigma' : 0.5}):\n",
    "    if kernel is gaussian_kernel:\n",
    "        return (kernel(X,X))\n",
    "    else:\n",
    "        return (kernel(X,X.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should return a $m\\times m$ numpy array. Note that the parameter \"kernel\" takes a function as input. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check, let's compute the kernel matrix for the first 3 features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.5734516   1.44127971  1.3565511 ]\n",
      " [ 1.44127971  1.32942692  1.24039027]\n",
      " [ 1.3565511   1.24039027  1.1765158 ]]\n"
     ]
    }
   ],
   "source": [
    "KM = gram_matrix(X=Xs[0:3], kernel=linear_kernel)\n",
    "print(KM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will  obtain an array similar to: <br>\n",
    "[ [ 1.5734516 &nbsp;  1.44127971 &nbsp; 1.3565511 ]<br>\n",
    " [ 1.44127971 &nbsp; 1.32942692 &nbsp; 1.24039027] <br>\n",
    " [ 1.3565511 &nbsp;  1.24039027 &nbsp; 1.1765158 ] ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM3$ and $\\bf TM4$ Given a $m$- dimensional numpy array with the lagrange multipliers (input variable $a$), the features matrix (input variable $X$), the labels array (input variable $y$), and the regularization parameter C (for softmargin svm), obtain the support vectors. (Hint: compare the multipliers with a small number instead of zero, e.g. 1e-5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def support_vectors(a,X,y,C):\n",
    "    sva, sv, svy = None, None, None\n",
    "    sv_indices = (a > 1e-5) & (a < C)\n",
    "    sva = a[sv_indices]\n",
    "    sv = X[sv_indices]\n",
    "    svy = y[sv_indices]\n",
    "    return sva,sv,svy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function should return 3 numpy arrays: a $k$- dimensional array with the lagrange multipliers for the support vectors, a $k \\times n$ numpy array with the support vectors, and a $k$ dimensional array with the labels for the support vectors. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try the function with the next code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  30.71    23.234  788.32 ]\n",
      "[[ 0.52168437  0.52414216  0.66214597  0.02717556  0.02799127  0.02858336\n",
      "   0.27215458  0.274725    0.43843728  0.00073851  0.00078351  0.00081701]\n",
      " [ 0.52086057  0.48314951  0.62939815  0.0350372   0.09734438  0.07128154\n",
      "   0.27129573  0.23343345  0.39614203  0.00122761  0.00947593  0.00508106]\n",
      " [ 0.54445806  0.55502451  0.68804466  0.03304641  0.09772327  0.07112868\n",
      "   0.29643458  0.30805221  0.47340546  0.00109207  0.00954984  0.00505929]]\n",
      "[[-1.]\n",
      " [ 1.]\n",
      " [ 1.]]\n"
     ]
    }
   ],
   "source": [
    "aaux = np.array([1.000e+03, 3.071e+01, 1.4039e+05, 1.00e+03, 23.234, 7.8832e+02, 1.000e+03, 1.232e+03])\n",
    "sa_aux, sv_aux, svy_aux = support_vectors(a=aaux, X=Xs[0:8], y=ys[0:8], C=900)\n",
    "# i,w = support_vectors(a=aaux, X=Xs[0:8], y=ys[0:8], C=900)\n",
    "# print(i)\n",
    "# b,w,sv,svy = support_vectors(a=aaux, X=Xs[0:8], y=ys[0:8], C=900)\n",
    "print(sa_aux)\n",
    "print(sv_aux)\n",
    "print(svy_aux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should have something similar to: <br>\n",
    "[  30.71 &nbsp;   23.234 &nbsp; 788.32 ]<br>\n",
    "[ [ 0.52168437 &nbsp; 0.52414216 &nbsp; 0.66214597 &nbsp; 0.02717556  &nbsp;0.02799127 &nbsp; 0.02858336<br>\n",
    "   0.27215458 &nbsp; 0.274725   &nbsp; 0.43843728 &nbsp; 0.00073851 &nbsp; 0.00078351 &nbsp; 0.00081701]<br>\n",
    " [ 0.52086057 &nbsp; 0.48314951 &nbsp; 0.62939815 &nbsp; 0.0350372 &nbsp;  0.09734438 &nbsp; 0.07128154<br>\n",
    "   0.27129573 &nbsp; 0.23343345 &nbsp; 0.39614203 &nbsp; 0.00122761 &nbsp; 0.00947593 &nbsp; 0.00508106]<br>\n",
    " [ 0.54445806 &nbsp; 0.55502451 &nbsp; 0.68804466 &nbsp; 0.03304641 &nbsp; 0.09772327 &nbsp; 0.07112868<br>\n",
    "   0.29643458 &nbsp; 0.30805221 &nbsp; 0.47340546 &nbsp; 0.00109207 &nbsp; 0.00954984 &nbsp; 0.00505929] ]<br>\n",
    "[ [-1.]<br>\n",
    " [ 1.]<br>\n",
    " [ 1.] ]<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's try the algorithm. Again, this part could take some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     pcost       dcost       gap    pres   dres\n",
      " 0: -1.9984e+03 -4.0528e+03  2e+03  2e-15  2e+00\n",
      " 1: -4.8986e+03 -4.9489e+03  5e+01  1e-13  1e+00\n",
      " 2: -4.8145e+05 -4.8150e+05  5e+01  9e-12  1e+00\n",
      " 3: -4.6359e+09 -4.6359e+09  5e+03  2e-08  1e+00\n",
      "Terminated (singular KKT matrix).\n",
      "2000 support vectors out of 2000 points\n",
      "667 out of 876 predictions correct\n"
     ]
    }
   ],
   "source": [
    "#Setting up the hyperparameters. \n",
    "C = None\n",
    "kernel_params = {\n",
    "    'a' : 1, \n",
    "    'c' : 1, \n",
    "    'd' : 3, \n",
    "    'sigma' : 0.5\n",
    "}\n",
    "#creating the model\n",
    "model = svm.SVM(kernel=gaussian_kernel, C=C)\n",
    "\n",
    "#finding the optimal margin\n",
    "#The function fit take as arguments:\n",
    "#the X feature matrix. \n",
    "#The y vector of labels. \n",
    "#The function to compute the kernel matrix\n",
    "#The function to compute the support vectors\n",
    "#A dictionary with the kernel parameters\n",
    "model.fit(X=Xs, y=ys, gram_matrix=gram_matrix, support_vectors=support_vectors, kernel_params=kernel_params)\n",
    "\n",
    "y_predict = model.predict(Xt)\n",
    "correct = np.sum(y_predict == yt[:,0])\n",
    "print(\"%d out of %d predictions correct\" % (correct, len(y_predict)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM3$ and $\\bf TM4$ Copy the last code in the next line and try with a value of C = 100 and sigma = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     pcost       dcost       gap    pres   dres\n",
      " 0: -1.0092e+05 -6.7488e+05  6e+05  1e-11  3e-12\n",
      " 1: -1.2990e+05 -2.0280e+05  7e+04  6e-12  2e-12\n",
      " 2: -1.9376e+05 -2.0028e+05  7e+03  4e-12  1e-12\n",
      " 3: -1.9435e+05 -1.9451e+05  2e+02  1e-13  2e-12\n",
      " 4: -1.9440e+05 -1.9440e+05  2e+00  5e-13  2e-11\n",
      " 5: -1.9440e+05 -1.9440e+05  2e-02  3e-11  7e-12\n",
      "Optimal solution found.\n",
      "2000 support vectors out of 2000 points\n",
      "667 out of 876 predictions correct\n"
     ]
    }
   ],
   "source": [
    "#Copy and change hyperparameters here\n",
    "C = 100\n",
    "kernel_params = {\n",
    "    'a' : 1, \n",
    "    'c' : 1, \n",
    "    'd' : 3, \n",
    "    'sigma' : 1.0\n",
    "}\n",
    "#creating the model\n",
    "model = svm.SVM(kernel=gaussian_kernel, C=C)\n",
    "\n",
    "#finding the optimal margin\n",
    "#The function fit take as arguments:\n",
    "#the X feature matrix. \n",
    "#The y vector of labels. \n",
    "#The function to compute the kernel matrix\n",
    "#The function to compute the support vectors\n",
    "#A dictionary with the kernel parameters\n",
    "model.fit(X=Xs, y=ys, gram_matrix=gram_matrix, support_vectors=support_vectors, kernel_params=kernel_params)\n",
    "\n",
    "y_predict = model.predict(Xt)\n",
    "correct = np.sum(y_predict == yt[:,0])\n",
    "print(\"%d out of %d predictions correct\" % (correct, len(y_predict)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\bf TM3$, $\\bf TM4$ Again copy and run the test with the polynomial kernel and C = 100. Keep the parameters a, c, and d to their default value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     pcost       dcost       gap    pres   dres\n",
      " 0:  2.0729e+05 -1.6044e+07  5e+07  1e+00  1e-11\n",
      " 1:  3.4324e+05 -5.5779e+06  8e+06  8e-02  1e-11\n",
      " 2:  1.8755e+05 -1.7660e+06  2e+06  2e-02  8e-12\n",
      " 3:  8.7444e+04 -8.0887e+05  1e+06  8e-03  7e-12\n",
      " 4:  3.4923e+04 -3.9424e+05  4e+05  3e-03  6e-12\n",
      " 5:  1.1596e+04 -2.0036e+05  2e+05  1e-03  6e-12\n",
      " 6:  1.7633e+03 -1.0651e+05  1e+05  5e-04  5e-12\n",
      " 7: -2.0712e+03 -7.4463e+04  7e+04  1e-04  6e-12\n",
      " 8: -2.5776e+03 -6.7242e+04  6e+04  9e-05  5e-12\n",
      " 9: -4.8111e+03 -4.5756e+04  4e+04  5e-05  5e-12\n",
      "10: -6.0293e+03 -3.2387e+04  3e+04  2e-05  5e-12\n",
      "11: -6.7371e+03 -2.8931e+04  2e+04  1e-05  6e-12\n",
      "12: -7.5910e+03 -2.4124e+04  2e+04  7e-06  6e-12\n",
      "13: -7.8689e+03 -2.2762e+04  1e+04  5e-06  6e-12\n",
      "14: -8.3269e+03 -1.9584e+04  1e+04  3e-06  7e-12\n",
      "15: -8.6954e+03 -1.8482e+04  1e+04  2e-06  7e-12\n",
      "16: -9.0228e+03 -1.7049e+04  8e+03  1e-06  7e-12\n",
      "17: -9.5818e+03 -1.5486e+04  6e+03  8e-07  7e-12\n",
      "18: -9.9215e+03 -1.4477e+04  5e+03  5e-07  7e-12\n",
      "19: -1.0033e+04 -1.3937e+04  4e+03  3e-07  7e-12\n",
      "20: -1.0287e+04 -1.3257e+04  3e+03  2e-07  8e-12\n",
      "21: -1.0375e+04 -1.2869e+04  2e+03  1e-07  7e-12\n",
      "22: -1.0629e+04 -1.2360e+04  2e+03  6e-08  7e-12\n",
      "23: -1.0766e+04 -1.2057e+04  1e+03  2e-08  7e-12\n",
      "24: -1.0942e+04 -1.1773e+04  8e+02  1e-08  7e-12\n",
      "25: -1.1034e+04 -1.1609e+04  6e+02  7e-09  7e-12\n",
      "26: -1.1126e+04 -1.1455e+04  3e+02  2e-09  8e-12\n",
      "27: -1.1199e+04 -1.1359e+04  2e+02  8e-10  8e-12\n",
      "28: -1.1235e+04 -1.1306e+04  7e+01  1e-12  1e-11\n",
      "29: -1.1261e+04 -1.1279e+04  2e+01  2e-12  9e-12\n",
      "30: -1.1268e+04 -1.1272e+04  3e+00  6e-14  9e-12\n",
      "31: -1.1270e+04 -1.1270e+04  4e-01  1e-12  1e-11\n",
      "32: -1.1270e+04 -1.1270e+04  2e-02  5e-13  1e-11\n",
      "33: -1.1270e+04 -1.1270e+04  2e-04  3e-13  1e-11\n",
      "Optimal solution found.\n",
      "134 support vectors out of 2000 points\n",
      "858 out of 876 predictions correct\n"
     ]
    }
   ],
   "source": [
    "#Copy and change hyperparameters here\n",
    "C = 100\n",
    "kernel_params = {\n",
    "    'a' : 1, \n",
    "    'c' : 1, \n",
    "    'd' : 3, \n",
    "    'sigma' : 0.5\n",
    "}\n",
    "#creating the model\n",
    "model = svm.SVM(kernel=polynomial_kernel, C=C)\n",
    "\n",
    "#finding the optimal margin\n",
    "#The function fit take as arguments:\n",
    "#the X feature matrix. \n",
    "#The y vector of labels. \n",
    "#The function to compute the kernel matrix\n",
    "#The function to compute the support vectors\n",
    "#A dictionary with the kernel parameters\n",
    "model.fit(X=Xs, y=ys, gram_matrix=gram_matrix, support_vectors=support_vectors, kernel_params=kernel_params)\n",
    "\n",
    "y_predict = model.predict(Xt)\n",
    "correct = np.sum(y_predict == yt[:,0])\n",
    "print(\"%d out of %d predictions correct\" % (correct, len(y_predict)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, you can test with different combinations of hyperparameters to improve the accuracy of the algorithm. The code presented here is a basic implementation of SVM. You can use a more robust library in your projects, e.g. LIBSVM (https://www.csie.ntu.edu.tw/~cjlin/libsvm/). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}