import numpy as np
from sklearn import svm
from sklearn import datasets

def RBF(x,z,sigma=100):
    return np.exp(-np.dot(x-z,x-z)/2/sigma**2)
#    return (np.dot(x,z)+1)**3
#    return np.dot(x,z)

class SVM():
    def __init__(self,X_train,y_train,X_test,C,epsilon=10**-8,max_itera=100):
        self.X_train=X_train
        self.y_train=y_train
        self.X_test=X_test
        train_no,feature_no=np.shape(self.X_train)
        test_no,_=np.shape(self.X_test)
        self.train_no=train_no
        self.test_no=test_no
        self.C=C
        self.epsilon=epsilon
        self.max_itera=max_itera

    def g(self,x,b):
        gx=0
        for i in range(self.train_no):
            gx+=self.alpha[i]*self.y_train[i]*RBF(self.X_train[i],x)
        return gx+b

    def E(self,x,y,b):
        gx=0
        for i in range(self.train_no):
#            if self.alpha[i]>0 and self.alpha[i]<self.C:
#                gx+=self.alpha[i]*self.y_train[i]*RBF(self.X_train[i],x)
            gx+=self.alpha[i]*self.y_train[i]*RBF(self.X_train[i],x)
        result=gx+b-y
        return result
    
    def W(self,index_1,index_2):
        K11=RBF(X_train[index_1],X_train[index_1])
        print(X_train[index_1],X_train[index_2])
        K22=RBF(X_train[index_2],X_train[index_2])
        K12=RBF(X_train[index_1],X_train[index_2])
        K21=RBF(X_train[index_2],X_train[index_1])       
        v1=self.g(X_train[index_1],self.b)-self.alpha[index_1]*self.y_train[index_1]*K11-self.alpha[index_2]*self.y_train[index_2]*K12-self.b
        v2=self.g(X_train[index_2],self.b)-self.alpha[index_1]*self.y_train[index_1]*K21-self.alpha[index_2]*self.y_train[index_2]*K22-self.b
        w=0.5*K11*self.alpha[index_1]**2+0.5*K22*self.alpha[index_2]**2+self.y_train[index_1]*self.y_train[index_2]*K12*self.alpha[index_1]*self.alpha[index_2] \
            -(self.alpha[index_1]+self.alpha[index_2])+self.y_train[index_1]*v1*self.alpha[index_1]+self.y_train[index_2]*v2*self.alpha[index_2]
        print(K11,K12,K21,K22,v1,v2,w)                
        return w
    
    def update_index_2(self,index_1,w_old,e):
        if e[index_1]>0:
            index_2=np.argmin(e)                
        else:
            index_2=np.argmax(e)
        print('index_2: ', index_2)

        w=self.W(index_1,index_2)

        if (w_old-w)<10**-8:

            for i in range(self.train_no):
                if self.alpha[i]>0 and self.alpha[i]<self.C:
                    w_temp=self.W(index_1,i)
                    if w_temp<w:
                        w=w_temp
                        index_2=i
        return w, index_2

    def is_stop(self,itera):
         

        if np.abs((np.dot(self.alpha,self.y_train))) <= self.epsilon:
            print("stop1:",np.abs((np.dot(self.alpha,self.y_train))))
#            print(self.alpha)
            if (all(self.alpha>=0) and all(self.alpha<=self.C))==True:
                result=np.zeros((self.train_no))
                for i in range(self.train_no):                        
                    if self.alpha[i]==0:
                        if self.y_train[i]*self.g(self.X_train[i],self.b)>=1:
                            print('stop2==0:',i ,y_train[i], self.g(self.X_train[i],self.b))
                            result[i]=1
                    elif self.alpha[i]==self.C:
                        if self.y_train[i]*self.g(self.X_train[i],self.b)<=1:
                            print('stop2==C:',i ,y_train[i], self.g(self.X_train[i],self.b))
                            result[i]=1
                    else:
                        if np.abs(self.y_train[i]*self.g(self.X_train[i],self.b)-1) <= self.epsilon:
                            print('stop2:',i ,y_train[i], self.g(self.X_train[i],self.b))
#                            print(self.alpha[0]*self.y_train[0]*RBF(self.X_train[0],self.X_train[0]))
                            result[i]=1
                if all(result==1)==True:
                    return True
                    
        if itera>=self.max_itera:
            return True
        return False
        
    def train(self):
        
#        self.alpha=np.random.rand((self.train_no))
        self.alpha=np.zeros((self.train_no))

        self.b=0
        e=np.zeros((self.train_no))
        w=np.inf
        itera=0

        while True:
#            print(y_train[1],self.g(self.X_train[1],self.b))
            flag_1st_via=True
#            list1=[]
            for i in range(self.train_no):
                e[i]=self.E(self.X_train[i],self.y_train[i],self.b)
            print(e)
#            print(self.alpha)
                
            for i in range(self.train_no):
                if self.alpha[i]>0 and self.alpha[i]<self.C:
                    print('a',i)
                    temp=0
                    if self.y_train[i]*self.g(self.X_train[i],self.b)!= 1:
                        if abs(self.y_train[i]*self.g(self.X_train[i],self.b) - 1)>temp:
                            index_1=i
                            temp=abs(self.y_train[i]*self.g(self.X_train[i],self.b) - 1)
                        flag_1st_via=False
#                        break

            if flag_1st_via==True:
                temp=0                
                for i in range(self.train_no):
                    if self.alpha[i]==0:
#                        print('b', i)
                        if self.y_train[i]*self.g(self.X_train[i],self.b) < 1:
                            if abs(self.y_train[i]*self.g(self.X_train[i],self.b) - 1)>temp:
                                index_1=i
                                temp=abs(self.y_train[i]*self.g(self.X_train[i],self.b) - 1)
                                print('index_1:',index_1)
                                
                    if self.alpha[i]==self.C:
                        if self.y_train[i]*self.g(self.X_train[i],self.b) > 1:
                            if abs(self.y_train[i]*self.g(self.X_train[i],self.b) - 1)>temp:
                                index_1=i
                                temp=abs(self.y_train[i]*self.g(self.X_train[i],self.b) - 1)

            w_old=w
            w,index_2=self.update_index_2(index_1,w_old,e)
            print('w:',w)
                      
            i=0
            while (w_old-w)<10**-8:
                w, index_2=self.update_index_2(i,w_old,e)
                index_1=i
                i+=1
                if i>=self.train_no:
                    break
                
                
            if self.y_train[index_1]!=self.y_train[index_2]:
                L=max(0, self.alpha[index_2]-self.alpha[index_1])
                H=min(self.C, self.C+self.alpha[index_2]-self.alpha[index_1])
            else:
                L=max(0, self.alpha[index_2]+self.alpha[index_1]-self.C)
                H=min(self.C, self.alpha[index_2]+self.alpha[index_1])
            print('L,H:',L,H)
                
            K11=RBF(X_train[index_1],X_train[index_1])
            K22=RBF(X_train[index_2],X_train[index_2])
            K12=RBF(X_train[index_1],X_train[index_2])
            K21=RBF(X_train[index_2],X_train[index_1])
            eta=K11+K22-2*K12            
            alpha2_old=self.alpha[index_2]            
            self.alpha[index_2]=self.alpha[index_2]+self.y_train[index_2]*(e[index_1]-e[index_2])/eta
            print('alpha2',self.alpha[index_2])

            if self.alpha[index_2]>H:
                self.alpha[index_2]=H
            if self.alpha[index_2]<L:
                self.alpha[index_2]=L
            print('alpha2',self.alpha[index_2])
                
            alpha1_old=self.alpha[index_1]  
            self.alpha[index_1]=self.alpha[index_1]+self.y_train[index_1]*self.y_train[index_2]*(alpha2_old-self.alpha[index_2])
            print('alpha1',self.alpha[index_1])
            
            
            
            
            
            
            
            b1=-e[index_1]-self.y_train[index_1]*K11*(self.alpha[index_1]-alpha1_old)-self.y_train[index_2]*K21*(self.alpha[index_2]-alpha2_old)+self.b
            b2=-e[index_2]-self.y_train[index_1]*K12*(self.alpha[index_1]-alpha1_old)-self.y_train[index_2]*K22*(self.alpha[index_2]-alpha2_old)+self.b 
            if (0<self.alpha[index_1]) and (self.alpha[index_1]<self.C):
                self.b=b1
            elif (0<self.alpha[index_2]) and (self.alpha[index_2]<self.C):
                self.b=b2
            else:
                self.b=(b1+b2)/2
            self.b=(b1+b2)/2
            print(b1,b2)
            itera+=1            
            if self.is_stop(itera):
                break
#            print('b:', self.b)        
            
#            print(itera,':',self.alpha,self.b,w)
        for j in range(self.train_no):
            if self.alpha[j]>0 and self.alpha[j]<self.C:
                summ=0
                for i in range(self.train_no):
                    summ+= self.y_train[i]*self.alpha[i]*RBF(X_train[i],X_train[j])
                self.b=self.y_train[j]-summ
                break
        w_coef=np.sum(np.reshape(self.alpha,(self.train_no,1)).T*np.reshape(self.y_train,(self.train_no,1)).T*X_train.T,axis=1)
#        print(self.alpha)
#        print(w_coef)
        print('b:', self.b)                  
        return self.alpha, self.b
            
    def predict(self):
        summ=0
        y_test=np.ones((self.test_no))
        for j in range(self.test_no):        
            for i in range(self.train_no):
                summ+= self.y_train[i]*self.alpha[i]*RBF(X_train[i],X_test[j])
            y_test[j]=summ+self.b
        y_test=np.sign(y_test)

        return y_test

     

if __name__=='__main__':
#    X_train=np.array([[3,3],[4,3],[1,1]])
#    y_train=np.array([1,1,-1])
#    X_test=np.array([[2,1],[8,9]])
    index=np.random.permutation(150)
    iris=datasets.load_iris()
    X_train=iris['data'][:,(2,3)]
    y_train=(iris['target']==2).astype(np.float64)
    y_train[y_train==0]=-1
    y_train_1=y_train
    X_test=X_train
   
    X_train=X_train[index]
    y_train=y_train[index]


    model=SVM(X_train,y_train,X_test,10)
    model.train()
    y_test=model.predict()
#    print('SVM', y_test)
    print(sum(y_test==y_train_1))
    model =svm.SVC(C=1)
#    model =svm.LinearSVC(C=1)
    model.fit(X_train, y_train)
    y_test=model.predict(X_test)
#    print('sklearn', y_test)
    print(sum(y_test==y_train_1))
#    print(model.coef_)
    print(model.intercept_)
#    print(model.get_params)
#    print(model.dual_coef_)

