#!/usr/bin/env python
# coding: utf-8

# ## Support Vector Machine

# In[790]:


import numpy as np
import math, random
import matplotlib.pyplot as plt
from scipy.optimize import minimize
global C, t, x, b, nonzero


# In[1000]:


# Generating the data sets
np.random.seed(100) # Will generate the same random data everytime you run the program

classA = np.concatenate(
    (np.random.randn(10 , 2) * 0.5 + [ 1.5 , 0.5] ,
     np.random.randn(10 , 2) * 0.5 + [ -1.5 , 0.5]))

classB = np.random.randn(20 , 2) * 0.5 + [0.0 , - 0.5]


inputs = np.concatenate((classA , classB))
targets = np.concatenate(
    (np.ones( classA.shape[0]) ,
     - np.ones( classB.shape[0])))

N = inputs.shape[0] #Number of rows (samples), .shape[checks number of rows at place 0]


permute=list(range (N))
random.shuffle(permute)
inputs = inputs[permute , : ] # List of 40 containing 2 elemens in each
targets = targets[ permute ]        # 1 or -1 Boolean value


# In[1001]:


def kernel(s,u):
    return np.dot(s,u) # Linear works

    #return np.power(np.dot(s, u) + 1, 5) # Polynomial
    #Radial Basis Function
    #e_d = - np.power(np.linalg.norm(s - u), 2)
    #return math.exp(e_d/(2*5**2))



# In[1002]:


def create_matrix(x,t):
    #x = np.array(N)
    #t = np.array(N)
    matrix = np.zeros((N,N))

    for n in range(N):
        for m in range(N):
            matrix[n][m] = targets[n] * targets[m] * kernel(inputs[n],inputs[m])

    return matrix

x=inputs
t=targets
matrix = create_matrix(x,t)


# In[ ]:





# In[ ]:





# In[ ]:





# In[1003]:


def objective(alpha):
    return 0.5*np.dot(alpha, np.dot(alpha, matrix)) - np.sum(alpha)


def zerofun(alpha):
    return np.dot(alpha, t)


def indicator(a, c):
    sum = 0
    for value in nonzero:
        sum += value[0] * value[2] * kernel([a,c],value[1])
    return sum - b


def main():
    # Initial guess
    alpha_zero = np.zeros(N)
    B =[(0, C) for b in range(N)] #Upper bound boundary
    XC = {'type':'eq', 'fun':zerofun}

    # Get solution
    sol = minimize( objective , alpha_zero , bounds = B, constraints = XC )
    return sol


b = 0
C = 1
sol = main()

alpha = sol['x']
not_zero = [(alpha[i], inputs[i], targets[i]) for i in range(N) if abs(alpha[i]) > 10e-5]


# In[1004]:


def b():
    sum = 0
    for j in nonzero:
        sum += j[0] * j[2] * kernel(j[1], not_zero[0][1])
    return sum - nonzero[0][2]

b = b()


# In[979]:


def plot():
    plt.plot([p [0] for p in classA ] ,
             [p [1] for p in classA ] ,
            'b.')
    plt.plot([p [0] for p in classB ] ,
             [p [1] for p in classB ] ,
            'r.' )

    plt.axis('equal') # Force same scale on both axes
    #plt.savefig('svmplot .pdf')

    # Decision boundaries
    xgrid=np.linspace( - 5, 5)
    ygrid=np.linspace ( - 4, 4)

    grid=np.array( [ [ indicator (h , g) for h in xgrid ] for g in ygrid ])

    plt.contour( xgrid , ygrid , grid , ( - 1.0, 0.0 , 1.0) , colors = ( 'r' , 'k' , 'b' ) , linewidths =(1 , 3 , 1))
    plt.show()



plot()


# In[809]:


plot()
# RBF


#

# In[821]:


plot()
# Linear Kernel


# In[826]:


plot()
# Polynomial p = 2


# In[832]:


# New generated data set - Polynomial
plot()


# ## Theory

# C \
# More slack allowed - lower c - Noisy data\
# Less slack allowed - higher c - Less noisy data
#

# In[932]:


plot()


# Linear - No solution obtained, high C so less slack allowed

# Non-linear kernels\
# - Sigma\
# Higher sigma, higher smoothness, less variance
#
# - p\
# Higher p, higher variance, less smoothness

# In[942]:


plot()


# Polynomial p = 2

# In[944]:


plot()


# Polynomial p = 5

# In[954]:


plot()


# sigma = 2

# In[961]:


plot()


# sigma = 5

# # Question

# Imagine that you are given data that is not easily separable. When
# should you opt for more slack rather than going for a more complex
# model (kernel) and vice versa?

# If the data are noisy (moving into each other) then it is maybe good to increase slack\
# but if it is more about how the spread of the data is then it is better to change the kernel
#

# In[1005]:


plot()


# In[991]:


plot()


# In[985]:


plot()


# In[ ]: