from cvxopt.solvers import qp
from cvxopt.base import matrix
import numpy as np, pylab, random, math


#kernel function
def linear_kernel(x, y):
    return np.dot(x,y) + 1

def radial_kernel(x, y):
    sigma=3
    return math.exp(-np.linalg.norm(((x[0] - y[0]), x[1] - y[1])) ** 2 / math.pow(2 * sigma, 2))

def polynomial_kernel(x, y):
    p = 5;
    return math.pow(np.dot(x, y) + 1, p)

def indicator(point, alpha, t, x, kernel):
    result = 0;
    for i in range(len(alpha)):
        result += alpha[i] * t[i] * kernel(x[i], point)
    return result

def createP(t, x, kernel):
    P = np.zeros( shape=(len(t), len(t)) )
    for i, ti in enumerate(t):
        for j, tj in enumerate(t):
            P[i][j] = ti * tj * kernel(x[i], x[j])

    return P

#create random dataset and store it in data
np.random.seed(100)
classA = [(random.normalvariate(-1.5, 1), random.normalvariate(0.5, 1), 1.0)
          for i in range(5)] + [(random.normalvariate(1.5, 1),
                                 random.normalvariate(0.5, 1), 1.0) for i in range(5)]
classB = [(random.normalvariate(0.0, 0.5),
           random.normalvariate(-1.5, 0.5), -1.0) for i in range(10)]
data = classA + classB
random.shuffle(data)

#plot dataset
pylab.hold(True)
pylab.plot([p[0] for p in classA], [p[1] for p in classA], 'bo')
pylab.plot([p[0] for p in classB], [p[1] for p in classB], 'ro')
pylab.show()

#calculate alpha for given dataset using given kernel
def find_alpha(data, kernel):
    tt = []
    Xx = []
    for d in data:
        tt.append(d[2])
        Xx.append((d[0], d[1]))

    N = len(data)

    q = np.full(N, -1)
    h = np.zeros(N)
    G = np.zeros((N, N))
    np.fill_diagonal(G, -1)
    P = createP(tt, Xx, kernel)
    q = q.astype(np.double)

    r = qp(matrix(P), matrix(q), matrix(G), matrix(h))
    tmp = list(r['x'])

    alpha = []
    t = []
    x = []
    for i in range(len(tmp)):
        x.append(Xx[i])
        t.append(tt[i])
        if tmp[i] < 0.00001:
            alpha.append(0)
        else:
            alpha.append(tmp[i])

    return (x, t, alpha)

(x, t, alpha) = find_alpha(data, linear_kernel)

#plot dataset with classifier
xrange = np.arange(-4, 4, 0.05)
yrange = np.arange(-4, 4, 0.05)
grid = matrix([[indicator([xd, yd], alpha, t, x, linear_kernel)
              for yd in yrange]
              for xd in xrange])
pylab.contour(xrange, yrange, grid,
              (-1.0, 0.0, 1.0),
              colors=('red', 'black', 'blue'),
              linewidths=(1, 3, 1))
pylab.plot([p[0] for p in classA], [p[1] for p in classA], 'bo')
pylab.plot([p[0] for p in classB], [p[1] for p in classB], 'ro')
pylab.show()

(x, t, alpha) = find_alpha(data, radial_kernel)

xrange = np.arange(-4, 4, 0.05)
yrange = np.arange(-4, 4, 0.05)
grid = matrix([[indicator([xd, yd], alpha, t, x, radial_kernel)
              for yd in yrange]
              for xd in xrange])
pylab.contour(xrange, yrange, grid,
              (-1.0, 0.0, 1.0),
              colors=('red', 'black', 'blue'),
              linewidths=(1, 3, 1))
pylab.plot([p[0] for p in classA], [p[1] for p in classA], 'bo')
pylab.plot([p[0] for p in classB], [p[1] for p in classB], 'ro')
pylab.show()

#calculate alpha with slack for given dataset using given kernel
def find_alpha_with_slack(data, kernel, C):
    tt = []
    Xx = []
    for d in data:
        tt.append(d[2])
        Xx.append((d[0], d[1]))

    N = len(data)

    q = np.full(N, -1)
    q = q.astype(np.double)

    h = np.concatenate((np.zeros(N), np.full(N, C)))
    G1 = np.zeros((N,N))
    G2 = np.zeros((N,N))
    np.fill_diagonal(G1, -1)
    np.fill_diagonal(G2, 1)
    G = np.concatenate((G1, G2))

    P = createP(tt, Xx, kernel)


    r = qp(matrix(P), matrix(q), matrix(G), matrix(h))
    tmp = list(r['x'])

    alpha = []
    t = []
    x = []
    for i in range(len(tmp)):
        x.append(Xx[i])
        t.append(tt[i])
        if tmp[i] < 0.00001 or tmp[i] > C:
            alpha.append(0)
        else:
            alpha.append(tmp[i])

    return (x, t, alpha)

(x, t, alpha) = find_alpha_with_slack(data, polynomial_kernel, 2)

xrange = np.arange(-4, 4, 0.05)
yrange = np.arange(-4, 4, 0.05)
grid = matrix([[indicator([xd, yd], alpha, t, x, polynomial_kernel)
              for yd in yrange]
              for xd in xrange])
pylab.contour(xrange, yrange, grid,
              (-1.0, 0.0, 1.0),
              colors=('red', 'black', 'blue'),
              linewidths=(1, 3, 1))
pylab.plot([p[0] for p in classA], [p[1] for p in classA], 'bo')
pylab.plot([p[0] for p in classB], [p[1] for p in classB], 'ro')
pylab.show()