import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pylab as pl
import cvxopt
import cvxopt.solvers
from sklearn.decomposition import PCA
from sklearn import svm
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns

##DEFINE THE SVM AND USE cvxopt TO SOLEVE QUADRATIC PROBLEM
class SVM(object):

    #define the kernel
    def linear_kernel(x1, x2):
        return np.dot(x1, x2)
    
    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

        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),'d')
        b = cvxopt.matrix(0.0)

        tP1 = np.diag(np.ones(n_samples) * -1)
        tP2 = np.identity(n_samples)
        G = cvxopt.matrix(np.vstack((tP1, tP2)))
        tP1 = np.zeros(n_samples)
        tP2 = np.ones(n_samples) * self.C
        h = cvxopt.matrix(np.hstack((tP1, tP2)))

        # solve QP problem
        solution = cvxopt.solvers.qp(P, q, G, h, A, b)
        cvxopt.solvers.options['show_progress'] = False

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

        # Support vectors have nonzero lamnda
        sv = a > 1e-3
        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))
        print("The support vectors are:",self.sv,sep='\n')

        # 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
        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]

    def project(self, X):
        return np.dot(X, self.w) + self.b
        
    def predict(self, X):
        return np.sign(self.project(X))
    
    def support_vectors(self, X):
        return self.sv

##DEFINE THE TEST OF SOFT MARGIN
def test_soft(X_train, Y_train,X_test, Y_test):
        clf = SVM(C=1000)
        clf.fit(X_train, Y_train)

        y_predict = clf.predict(X_test)
        correct = np.sum(y_predict == Y_test)
        test = np.where(y_predict == Y_test,y_predict,0)
        accuracy = correct/(len(y_predict))
        #print("\n\nCompare differences between the test data and the solution through the SVM:",test,sep='\n' )
        print("%d out of %d predictions correct" % (correct, len(y_predict)))
        print("The accuracy of the SVM is:",accuracy)
        
##DATA IMPORT & CLEANING
BCW = pd.read_csv('breast-cancer-wisconsin.csv')
BCW.drop(['ID'], 1, inplace = True)
BCW_data = np.array(BCW.drop(['Class'],1))
BCW_result = np.array(BCW['Class'])

#For 3D
pca3 = PCA(n_components=3).fit(BCW_data)
BCW_3d = pca3.transform(BCW_data)
#separate into training and testing
X_train3 = BCW_3d[:120]
X_test3 = BCW_3d[120:]

#For 2D
pca2 = PCA(n_components=2).fit(BCW_data)
BCW_2d = pca2.transform(BCW_data)
#separate into training and testing
X_train2 = BCW_2d[:120]
X_test2 = BCW_2d[120:]

# separate into training and testing
Y_train = BCW_result[:120]
Y_test = BCW_result[120:]

print("The result of 2D:")
test_soft(X_train2, Y_train, X_test2, Y_test)
print("\nThe result of 3D:")
test_soft(X_train3, Y_train, X_test3, Y_test)

#hyperplane
clf = svm.SVC(kernel='linear', C=1000)
clf.fit(X_train2, Y_train)

# plt.scatter(X_train2[:, 0], X_train2[:, 1], c=Y_train, s=50, cmap=plt.cm.Paired)

for i in range(0, 120):
    if Y_train[i] == 1:
        d1=plt.scatter(X_train2[i,0], X_train2[i,1],  c='r', marker='o')
    else:
        d2=plt.scatter(X_train2[i,0], X_train2[i,1],  c='g', marker='o')
        
plt.legend([d1,d2],['benign','malignant'],loc='upper left')


# plot the decision function
ax = plt.gca()
xlim = ax.get_xlim()
ylim = ax.get_ylim()

# create grid to evaluate model
xx=np.linspace(xlim[0], xlim[1], 100)
yy=np.linspace(ylim[0], ylim[1], 100)
YY, XX = np.meshgrid(yy, xx)
xy = np.vstack([XX.ravel(), YY.ravel()]).T
Z = clf.decision_function(xy).reshape(XX.shape)

# plot decision boundary and margins
ax.contour(XX, YY, Z, colors='k', levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])

# plot support vectors
ax.scatter(clf.support_vectors_[:, 0], clf.support_vectors_[:, 1], s=60,
           linewidth=1, facecolors='none', edgecolors='k')

plt.show()

fig = plt.figure()
ax = fig.gca(projection='3d')

#Training
cls = svm.SVC(kernel='linear', C=1000)
cls.fit(X_train3, Y_train)

w = cls.coef_
b = cls.intercept_

#Plot Hyperplane
ax = plt.subplot(111, projection='3d')
x = np.arange(-10,10,0.01)
y = np.arange(-10,10,0.11)
x, y = np.meshgrid(x, y)
z = (w[0,0]*x + w[0,1]*y + b) / (-w[0,2])
surf = ax.plot_surface(x, y, z)

# Plot scatter
x_array = np.array(X_train3, dtype=float)
y_array = np.array(Y_train, dtype=int)
be = x_array[np.where(y_array==1)]
ma = x_array[np.where(y_array==-1)]
ax.scatter(be[:,0], be[:,1], be[:,2], c='r', label='benign')
ax.scatter(ma[:,0], ma[:,1], ma[:,2], c='g', label='malignant')

# Plot support vectors
X = np.array(X_train3,dtype=float)
for i in range(len(sv_idx)):
     ax.scatter(X[sv_idx[i],0], X[sv_idx[i],1], X[sv_idx[i],2],s=30,
                c='',marker='o', edgecolors='k')

# Set lables
ax.set_zlabel('Z')
ax.set_ylabel('Y')
ax.set_xlabel('X')
ax.set_zlim([-90, 90])
plt.legend(loc='upper left')

plt.show()

#Plot the heat map to see the correlation among attributes.
sns.set(style='ticks', color_codes=True)
plt.figure(figsize=(14, 12))
sns.heatmap(BCW.astype(float).corr(), 
            linewidths=0.1, 
            square=True, 
            linecolor='white', 
            annot=True)
plt.show()

#Bar plot 1 
fig = plt.figure()
ax = sns.countplot(x="Bland_Chromatin", 
                   hue='Class', 
                   #palette={0:'#EB434A', 1:'#61A98F'}, 
                   data=BCW)
ax.set(xlabel='Bland Chromatin', ylabel='No of cases')
fig.suptitle("Bland Chromatin w.r.t. Class", y=0.96);

#Bar plot 2 
fig = plt.figure()
ax = sns.countplot(x='Clump_Thickness', 
                   hue='Class', 
                   #palette={0:'#EB434A', 1:'#61A98F'}, 
                   data=BCW)
ax.set(xlabel='Thickness of Clump', ylabel='Total')
fig.suptitle("Clump_Thickness w.r.t. Class", y=0.96);