from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import random
from random import randint



##TASK 2.1
def OneStep(x,y,beta,b,f_values,i,j,C):
	delta=y[i]*((f_values[j]-y[j])-(f_values[i]-y[i]))
	s=y[i]*y[j]
	chi=np.dot(x[i],x[i])+np.dot(x[j],x[j])-2*np.dot(x[i],x[j])
	gamma=s*beta[i]+beta[j]

	if s==1:
		L=max(0,gamma-C)
		H=min(gamma,C)
	else:
		L=max(0,-gamma)
		H=min(C-gamma,C)

	if chi>0:
		betaI=min(max(beta[i]+delta/chi,L),H)
	elif delta>0:
		betaI=L
	else:
		betaI=H

	betaJ=gamma-s*betaI

	k=0
	for point in x:
		f_values[k]=f_values[k]+(betaI-beta[i])*y[i]*np.dot(point,x[i])+(betaJ-beta[j])*y[j]*np.dot(point,x[j])
		k=k+1

	bNew=b-0.5*(f_values[i]-y[i]+f_values[j]-y[j])
	f_values=f_values+np.ones(len(f_values))*(bNew-b)
	return betaI,betaJ,bNew,f_values


##TASK 2.2

coor11=np.random.exponential(4,20)
coor12=np.random.exponential(4,20)
x1=np.array([coor11,coor12]).T
coor21=np.random.exponential(0.5,20)
coor22=np.random.exponential(0.5,20)
x2=np.array([coor21,coor22]).T
x=np.concatenate((x1,x2))
y=np.concatenate((np.ones(20,int),-1*np.ones(20,int)))


#TASK 2.5

def KKT(i,C,beta,y,f_values):
	return C*max(0,1-y[i]*f_values[i])-beta[i]*(1-y[i]*f_values[i])

def SMO2(x,y,n,C,iterations):
	beta=np.zeros(n,int)
	b=0
	f_values=np.zeros(n,int)

	list_0bC=np.array([])
	i=0
	while i<n and iterations>0:
		if KKT(i,C,beta,y,f_values)>0:
			if len(list_0bC)!=0:
				w=np.where(list_0bC==i)[0]
				if len(w)!=0:
					j=random.choice(np.delete(list_0bC,w[0]))
				else:
					j=random.choice(list_0bC)

			else:
				j=random.choice(np.delete(range(n),i))

			beta[i],beta[j],b,f_values=OneStep(x,y,beta,b,f_values,i,j,C)

			if C>beta[i]>0:
				if not np.any(list_0bC[:]==i):
					list_0bC=np.concatenate((list_0bC,[i]))

			if C>beta[j]>0:
				if not np.any(list_0bC[:]==j):
					list_0bC=np.concatenate((list_0bC,[j]))

			i=0


		else:
			i+=1

		iterations-=1



	prov=np.array([])
	for k in xrange(n):
		if beta[k]>0:
			prov=np.concatenate((prov,[f_values[k]-y[k]]))

	if len(prov)!=0:
		b=b-np.median(prov)

	return beta,b

def LLS(x,y):
	X_matrix=np.vstack([np.ones(len(x)),x.T]).T
	return np.dot(  np.linalg.inv(np.dot(X_matrix.T,X_matrix))   ,  np.dot(X_matrix.T,y.T))


def RunSMO2(x,y,n,C,iterations):
	
	
	beta,b=SMO2(x,y,n,C,iterations)
	indSuppVec=np.array([0])
	indMargVec=np.array([0])

	for k in xrange(n):
		if beta[k]>0:
			indSuppVec=np.concatenate((indSuppVec,[k]))
			if beta[k]<C:
				indMargVec=np.concatenate((indMargVec,[k]))

	
	if len(indMargVec)!=1:
		plt.scatter(np.take(x.T[0],indMargVec[1:]),np.take(x.T[1],indMargVec[1:]),color='grey', marker='s',s=50)


	plt.plot(np.take(x.T[0],xrange(0,int(n/2))), np.take(x.T[1],xrange(0,int(n/2))),'go', np.take(x.T[0],xrange(int(n/2),n)), np.take(x.T[1],xrange(int(n/2),n)),'mo')

	if len(indSuppVec)!=1:
		plt.scatter(np.take(x.T[0],indSuppVec[1:]),np.take(x.T[1],indSuppVec[1:]), color='k',marker='x',s=100)


	alpha1=0
	alpha2=0
	for k in xrange(n):
		alpha1=alpha1+beta[k]*y[k]*x[k][0]
		alpha2=alpha2+beta[k]*y[k]*x[k][1]


		
	if alpha2!=0:
		plt.plot(x, (-alpha1*x-b)/alpha2, 'r')
		
	elif alpha1!=0:
		plt.axvline(x=-b/alpha1, color='r')
	

	delta0, delta1, delta2= LLS(x,y)

	
	plt.plot(x, (-delta0-delta1*x)/delta2, 'b')

	plt.show()

	return b,alpha1,alpha2


RunSMO2(x,y,40,0.01,10000)


coor11=np.random.exponential(4,1000)
coor12=np.random.exponential(4,1000)
x1=np.array([coor11,coor12]).T
coor21=np.random.exponential(0.5,1000)
coor22=np.random.exponential(0.5,1000)
x2=np.array([coor21,coor22]).T
xTest=np.concatenate((x1,x2))
yTest=np.concatenate((np.ones(1000,int),-1*np.ones(1000,int)))

def Accuracy(xTest,alpha0, alpha1,alpha2):
	Acc1, Acc2=0,0
	for i in xrange(0,1000):
		if (np.dot(xTest[i],[alpha1,alpha2])+alpha0)>0:
			Acc1+=1

	for i in xrange(1000,2000):
		if (np.dot(xTest[i],[alpha1,alpha2])+alpha0)<0:
			Acc2+=1

	return (Acc1+Acc2)/2000

print "Accuracy LLS:"
a0,a1,a2=LLS(xTest,yTest)
print Accuracy(xTest,a0,a1,a2)
print "Accuracy SMO 1:"
a0,a1,a2=RunSMO2(x,y,40,1,10000)
print Accuracy(xTest,a0,a1,a2)
print "Accuracy SMO 100:"
a0,a1,a2=RunSMO2(x,y,40,100,10000)
print Accuracy(xTest,a0,a1,a2)