svm-kernel.py

# %%
# Mathieu Blondel, September 2010
# License: BSD 3 clause
# http://www.mblondel.org/journal/2010/09/19/support-vector-machines-in-python/

# visualizing what translating to another dimension does
# and bringing back to 2D:
# https://www.youtube.com/watch?v=3liCbRZPrZA

# Docs: http://cvxopt.org/userguide/coneprog.html#quadratic-programming
# Docs qp example: http://cvxopt.org/examples/tutorial/qp.html

# Nice tutorial:
# https://courses.csail.mit.edu/6.867/wiki/images/a/a7/Qp-cvxopt.pdf


import timeit
import matplotlib.pyplot as plt
import random
import numpy as np
from numpy import linalg
import cvxopt
import cvxopt.solvers
from sklearn import datasets


def linear_kernel(x1, x2):
    return np.dot(x1, x2)


def polynomial_kernel(x, y, p=3):
    return (1 + np.dot(x, y)) ** p


def gaussian_kernel(x, y, sigma=5.0):
    return np.exp(-linalg.norm(x-y)**2 / (2 * (sigma ** 2)))


class SVM(object):

    def __init__(self, kernel=linear_kernel, C=None):
        self.kernel = kernel
        self.C = C
        if self.C is not None:
            self.C = float(self.C)

    def fit(self, X, y):
        n_samples, n_features = X.shape

        # Gram matrix
        K = np.zeros((n_samples, n_samples))
        for i in range(n_samples):
            for j in range(n_samples):
                K[i, j] = self.kernel(X[i], X[j])

        P = cvxopt.matrix(np.outer(y, y) * K)
        q = cvxopt.matrix(np.ones(n_samples) * -1)
        A = cvxopt.matrix(y, (1, n_samples))
        b = cvxopt.matrix(0.0)

        if self.C is None:
            G = cvxopt.matrix(np.diag(np.ones(n_samples) * -1))
            h = cvxopt.matrix(np.zeros(n_samples))
        else:
            tmp1 = np.diag(np.ones(n_samples) * -1)
            tmp2 = np.identity(n_samples)
            G = cvxopt.matrix(np.vstack((tmp1, tmp2)))
            tmp1 = np.zeros(n_samples)
            tmp2 = np.ones(n_samples) * self.C
            h = cvxopt.matrix(np.hstack((tmp1, tmp2)))

        # solve QP problem
        solution = cvxopt.solvers.qp(P, q, G, h, A, b)

        # Lagrange multipliers
        a = np.ravel(solution['x'])

        # Support vectors have non zero lagrange multipliers
        sv = a > 1e-5
        ind = np.arange(len(a))[sv]
        self.a = a[sv]
        self.sv = X[sv]
        self.sv_y = y[sv]
        print("%d support vectors out of %d points" % (len(self.a), n_samples))

        # Intercept
        self.b = 0
        for n in range(len(self.a)):
            self.b += self.sv_y[n]
            self.b -= np.sum(self.a * self.sv_y * K[ind[n], sv])
        self.b /= len(self.a)

        # Weight vector
        if self.kernel == linear_kernel:
            self.w = np.zeros(n_features)
            for n in range(len(self.a)):
                self.w += self.a[n] * self.sv_y[n] * self.sv[n]
        else:
            self.w = None

    def project(self, X):
        if self.w is not None:
            return np.dot(X, self.w) + self.b
        else:
            y_predict = np.zeros(len(X))
            for i in range(len(X)):
                s = 0
                for a, sv_y, sv in zip(self.a, self.sv_y, self.sv):
                    s += a * sv_y * self.kernel(X[i], sv)
                y_predict[i] = s
            return y_predict + self.b

    def predict(self, X):
        return (self.project(X))


if __name__ == "__main__":
    import pylab as pl

    def plot_margin(X1_train, X2_train, clf):
        def f(x, w, b, c=0):
            # given x, return y such that [x,y] in on the line
            # w.x + b = c
            return (-w[0] * x - b + c) / w[1]

        pl.plot(X1_train[:, 0], X1_train[:, 1], "ro")
        pl.plot(X2_train[:, 0], X2_train[:, 1], "bo")
        pl.scatter(clf.sv[:, 0], clf.sv[:, 1], s=100, c="g")

        # w.x + b = 0
        a0 = -4
        a1 = f(a0, clf.w, clf.b)
        b0 = 4
        b1 = f(b0, clf.w, clf.b)
        pl.plot([a0, b0], [a1, b1], "k")

        # w.x + b = 1
        a0 = -4
        a1 = f(a0, clf.w, clf.b, 1)
        b0 = 4
        b1 = f(b0, clf.w, clf.b, 1)
        pl.plot([a0, b0], [a1, b1], "k--")

        # w.x + b = -1
        a0 = -4
        a1 = f(a0, clf.w, clf.b, -1)
        b0 = 4
        b1 = f(b0, clf.w, clf.b, -1)
        pl.plot([a0, b0], [a1, b1], "k--")

        pl.axis("tight")
        pl.show()

    def plot_contour(X1_train, X2_train, clf):
        pl.plot(X1_train[:, 0], X1_train[:, 1], "ro")
        pl.plot(X2_train[:, 0], X2_train[:, 1], "bo")
        pl.scatter(clf.sv[:, 0], clf.sv[:, 1], s=100, c="g")

        X1, X2 = np.meshgrid(np.linspace(-6, 6, 50), np.linspace(-6, 6, 50))
        X = np.array([[x1, x2] for x1, x2 in zip(np.ravel(X1), np.ravel(X2))])
        Z = clf.project(X).reshape(X1.shape)
        pl.contour(X1, X2, Z, [0.0], colors='k', linewidths=1, origin='lower')
        pl.contour(X1, X2, Z + 1, [0.0], colors='grey',
                   linewidths=1, origin='lower')
        pl.contour(X1, X2, Z - 1, [0.0], colors='grey',
                   linewidths=1, origin='lower')

        pl.axis("tight")
        pl.show()


# %%
iris = datasets.load_iris()
X = iris.data[:, 1:3]
# print(X)
y = iris.target


def oneversusrest(X, testdata):

    setosa = X[:50]
    versicolor = X[50:100]
    virginica = X[100:]

    # Get 10 Random Samples From Each Class
    random_setosa = random.choices(setosa, k=10)
    random_versicolor = random.choices(versicolor, k=10)
    random_virginica = random.choices(virginica, k=10)

    # CLASSIFIER 1
    join1 = np.vstack((random_virginica, random_versicolor))
    X_train = np.vstack((random_setosa, join1))

    y1 = np.ones(len(random_setosa))
    y2 = np.ones(len(join1)) * -1
    y_train1 = np.hstack((y1, y2))
    clf1 = SVM(C=1000.1)
    clf1.fit(X_train, y_train1)

    print("prediction1 : ", clf1.predict(testdata))
    # plot_contour(X_train[y_train1 == 1], X_train[y_train1 == -1], clf1)
    # =============================================================

    # CLASSIFIER 2
    join2 = np.vstack((random_setosa, random_virginica))
    X_train2 = np.vstack((random_versicolor, join2))

    y1 = np.ones(len(random_versicolor))
    y2 = np.ones(len(join2)) * -1
    y_train2 = np.hstack((y1, y2))
    clf2 = SVM(gaussian_kernel)
    clf2.fit(X_train2, y_train2)
    print("prediction2 ,", np.average(clf2.predict(testdata)))
    # plot_contour(X_train2[y_train2 == 1], X_train2[y_train2 == -1], clf2)
    # =============================================================

    # CLASSIFIER 3
    join3 = np.vstack((random_setosa, random_versicolor))
    X_train3 = np.vstack((random_virginica, join3))

    y1 = np.ones(len(random_virginica))
    y2 = np.ones(len(join3)) * -1
    y_train3 = np.hstack((y1, y2))
    clf3 = SVM(C=1000.1)
    clf3.fit(X_train3, y_train3)
    print("Prediction3,", clf3.predict(testdata))
    # plot_contour(X_train3[y_train3 == 1], X_train3[y_train3 == -1], clf3)
    # =============================================================

    predictions = [clf1.predict(testdata), np.average(
        clf2.predict(testdata)), clf3.predict(testdata)]
    ovrpredict = predictions.index(max(predictions))
    print("END PREDICTION 1 VS Rest: Class =", ovrpredict+1)

    return ovrpredict
# =============================================================
# %%


def DAG(X, testData):
    # %%

    setosa = X[:50]
    versicolor = X[50:100]
    virginica = X[100:]

    # Get 10 Random Samples From Each Class
    random_setosa = random.choices(setosa, k=10)
    random_versicolor = random.choices(versicolor, k=10)
    random_virginica = random.choices(virginica, k=10)

    # CLASSIFIER 1
    X_train = np.vstack((random_setosa, random_virginica))

    y1 = np.ones(len(random_setosa))
    y2 = np.ones(len(random_virginica)) * -1
    y_train1 = np.hstack((y1, y2))
    clf1 = SVM(C=1000.1)
    clf1.fit(X_train, y_train1)

    # plot_contour(X_train[y_train1 == 1], X_train[y_train1 == -1], clf1)
    res1 = np.sign(clf1.predict(testData))

    def clf_2(data):
        X_train2 = np.vstack((random_setosa, random_versicolor))

        y1 = np.ones(len(random_setosa))
        y2 = np.ones(len(random_versicolor)) * -1
        y_train2 = np.hstack((y1, y2))
        clf2 = SVM(C=1000.1)
        clf2.fit(X_train2, y_train2)
        print("prediction2 ,", np.average(clf2.predict(data)))
        # plot_contour(X_train2[y_train2 == 1], X_train2[y_train2 == -1], clf2)
        predicted = clf2.predict(data)
        if np.sign(predicted) != 1:
            return 2
        else:
            return 1

    def clf_3(data):
        X_train3 = np.vstack((random_versicolor, random_virginica))

        y1 = np.ones(len(random_versicolor))
        y2 = np.ones(len(random_virginica)) * -1
        y_train3 = np.hstack((y1, y2))
        clf3 = SVM(C=1000.1)
        clf3.fit(X_train3, y_train3)
        predicted = clf3.predict(data)

        print("prediction2 ,", np.average(clf3.predict(data)))
        # plot_contour(X_train3[y_train3 == 1], X_train3[y_train3 == -1], clf3)
        if np.sign(predicted) != 1:
            return 3
        else:
            return 2
    if res1 != 1:
        print("HASIL AKHIR METODE DAG:", clf_3(testData))
    else:

        print("HASIL AKHIR METODE DAG :", clf_2(testData))
    # =============================================================

    # CLASSIFIER 2


# X = iris.data
# print(X[:, 0])
# f = plt.figure()
# s = plt.scatter(X[:, 0], X[:, 1], c='y')
# plt.show(s)
# Jalanin Satu Aja Kalo bingung
# oneversusrest(X, X[101])
# DAG(X, X[101])
def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped


short_list = range(10)
wrapped = wrapper(oneversusrest, X, X[101])
# print(timeit.timeit(wrapped, number=1))
wrapped = wrapper(oneversusrest, X, X[101])

print("Function took: ", timeit.timeit(wrapped, number=1))