code_for_samples.py

# -*- coding: utf-8 -*-
"""FourClass2.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1pB13hzKgS3zazyAnuDAU08NYutGvyMQS

#Questions

2.Z_out vs Z_0 - singular matrix - starting value of Z ## fixed

3. Newton2 - do we start with 5 iterations or start with the while loop

4. Do we set alpha set to 1 or do we use changed alpha in each iteration of while loop- just to make sure
"""

# Commented out IPython magic to ensure Python compatibility.
import sys
import numpy
import cvxpy as cp
from scipy.optimize import nnls
from pandas import *
import random
from scipy.optimize import minimize, rosen, rosen_der
from scipy.optimize import least_squares
import copy
import math
import matplotlib.pyplot as plt
from numpy.linalg import matrix_rank
from google.colab import files

!pip install ipython-autotime
# %load_ext autotime

def print_Xi_Si(Z,n,m):
  for i in range(n):
    print("Xi: ",i," : ",Z[m+i])
  for i in range(n):
    print("S: ",i," : ",Z[m+i+n])

def generate_X(X):
  new_x = []
  for i in range(len(X)):
    new_x.append(numpy.concatenate((X[i], [1])))
  return new_x

def generate_A_tilda(X,y):
  n= len(X)
  m = len(X[0])
  A_tilda = numpy.zeros(shape=(n,m))
  for i in range(n):
    for j in range(m):
      A_tilda[i][j]= y[i]*X[i][j]
  return A_tilda

def generate_A(X,y):
  n= len(X)
  m = len(X[0])
  A = numpy.zeros(shape=(n,m+2*n))
  for i in range(n):
    A[i][m+i] =1 #This is for \xi
    A[i][m+n+i] = -1 # this is for Si
    for j in range(m):
      A[i][j]= y[i]*X[i][j]
  return A

def generate_M(A):
  n = len(A)
  m = len(A[0])-2*n
  temp = numpy.zeros(shape=(m+2*n,m+2*n))
  for i in range(m-1):
    temp[i][i]=1
  gammas = numpy.zeros((m+2*n,n))
  mus = numpy.zeros((m+2*n,n))
  A_transpose = numpy.transpose(A)
  zero1 = numpy.zeros(shape = (n,2*n))
  zero2 = numpy.zeros(shape = (n,n))

  top = numpy.concatenate((temp,gammas,mus,A_transpose),axis=1) 
  bottom = numpy.concatenate((A,zero1,zero2),axis=1)
  result = numpy.concatenate((top,bottom),axis=0)
  return result

#The order of Z is W,b,X_i,S_i,Gamma_i,Mu_i, Lambda_i
def generate_F_Z(M,Z,n,m):
  Z_out = numpy.zeros(m+5*n)

  Xi=[]
  Gamma=[]
  S=[]
  Mu =[]
  for i in range(n):
    Xi.append(Z[i+m])
    S.append(Z[i+m+n])  
    Gamma.append(Z[i+m+2*n])
    Mu.append(Z[i+m+3*n])
  M_z = numpy.matmul(M,Z)
  for i in range(len(M_z)):
    Z_out[i] = M_z[i]
  for i in range(n):
    Z_out[i+m+3*n]= (Gamma[i]*Xi[i])
    Z_out[i+m+4*n] = (Mu[i]*S[i])
  return Z_out



def get_output_vector(n,m):
  output_vector = numpy.zeros(m+3*n)
  for i in range(n):
    output_vector[m+2*n+i]=1
  return output_vector
  
def generate_d(n,m,beta):
  d = numpy.zeros(m+5*n)
  for i in range(n):
    d[m+2*n+i]=1
    d[m+3*n+i] = beta
    d[m+4*n+i] = beta
  return d

def get_Z_0(M,output_vector):
  rows = len(M)
  cols = len(M[0])
  n = (int)((cols-rows)/2)
  m = cols-5*n
  def fn(Z):
    temp = numpy.matmul(M,Z)
    for i in range(len(temp)):
      temp[i] = temp[i]-output_vector[i] 
    return temp
  
  Z_0 = numpy.zeros(m+5*n)
  res = least_squares(fn, Z_0)
  Z_0 = res.x
  return Z_0

def get_null(M,index): # gives a null vector where index of Z_0 is 1
  rows = len(M)
  cols = len(M[0])
  n = (int)((cols-rows)/2)
  m = cols-5*n
  def fn(Z):
    for i in range(2*n):
      Z[m+i]=0
    Z[index]=1
    temp = numpy.matmul(M,Z)
    return temp
  
  Z_0 = numpy.zeros(m+5*n)
  res = least_squares(fn, Z_0)
  Z_0 = res.x
  for i in range(2*n):
      Z_0[m+i]=0
  Z_0[index]=1
  return Z_0

def convertToPositive(M,Z_0,n,m): 
  Z_0_changed = copy.deepcopy(Z_0)
  for i in range(m,m+2*n):
    if(Z_0[i]<0):
      k=math.ceil(abs(Z_0[i]))
      null_vector = list(map(lambda x: k*x ,get_null(M,i)))
      Z_0_changed = [Z_0_changed[j] + null_vector[j] for j in range(len(Z_0))]  
  return Z_0_changed

#The order of Z is W,b,X_i,S_i,Gamma_i,Mu_i, Lambda_i
def Jacobian( n,m,A_tilda,Z):
  Xi = []
  S=[]
  Gamma=[]
  Mu =[]
  Lambda = []
  for i in range(n):
    Xi.append(Z[i+m])
    S.append(Z[i+m+n])  
    Gamma.append(Z[i+m+2*n])
    Mu.append(Z[i+m+3*n])
    Lambda.append(Z[i+m+4*n])
  Onn = numpy.zeros((n,n))
  Omn = numpy.zeros((m,n))
  Onm = numpy.zeros((n,m))
  Inn = numpy.identity(n)
  Im0 = numpy.identity(m)
  Gamma_nn = numpy.identity(n)
  Xi_nn = numpy.identity(n)
  Mu_nn = numpy.identity(n)
  S_nn = numpy.identity(n)
  negative_Inn = numpy.identity(n)
  Im0[m-1][m-1]=0

  for i in range(n):
    for j in range(n):
      negative_Inn[i][i]=-1
      Gamma_nn[i][i] = Gamma[i]
      Xi_nn[i][i]=Xi[i]  
      Mu_nn[i][i] = Mu[i]
      S_nn[i][i]=S[i]
  negativeAT = numpy.transpose(A_tilda)
  for i in range(len(negativeAT)):
    for j in range(len(negativeAT[i])):
      negativeAT[i][j]= -negativeAT[i][j]
  rowBlocks = [[],[],[],[],[],[]]
  rowBlocks[0] = numpy.concatenate((Im0,Omn,Omn,Omn,Omn,negativeAT),axis=1)     # column wise addition
  rowBlocks[1] = numpy.concatenate((Onm,Onn,Onn,negative_Inn,Onn,Inn),axis=1)
  rowBlocks[2] = numpy.concatenate((Onm,Onn,Onn,Onn,negative_Inn,negative_Inn),axis=1)
  rowBlocks[3] = numpy.concatenate((A_tilda,Inn,negative_Inn,Onn,Onn,Onn),axis=1)
  rowBlocks[4] = numpy.concatenate((Onm,Gamma_nn,Onn,Xi_nn,Onn,Onn),axis=1)
  rowBlocks[5] = numpy.concatenate((Onm,Onn,Mu_nn,Onn,S_nn,Onn),axis=1)
  jacobian_result = numpy.concatenate((rowBlocks))
  b = numpy.zeros(n+m)
  for i in range(m,m+n):
    b[i]=-1
  #print (DataFrame(jacobian_result))
  return jacobian_result

def calculate_gap(Z,d,n,m,M):
  F_Z = generate_F_Z(M,Z,n,m)
  temp = numpy.subtract(F_Z,d)
  temp_sum =0
  for i in range(len(temp)):
    temp_sum+=temp[i]*temp[i]
  return temp_sum

def inverse(Jacobian):
  l = len(Jacobian)

  def fn(v):
    inv = v.reshape(l,l)
    I = numpy.identity(l)
    temp_result= numpy.subtract(numpy.matmul(Jacobian,inv),I)
    return temp_result.flatten()
  
  if(matrix_rank(Jacobian)==len(Jacobian)):
    return numpy.linalg.inv(Jacobian)
  v = numpy.zeros(l*l)
  res = least_squares(fn, v)
  inv = res.x.reshape(l,l)
  return inv

def check_xi_si_negative(Z,n,m):
  for i in range(m,m+2*n):
    if(Z[i]<0):
      return False
  return True

def newton(iterations,X,y,Z0,beta,C,printbool):
  Z = Z0
  A = generate_A(X,y)
  n = len(A)
  m = len(A[0])-2*n
  M = generate_M(A)
  d = generate_d(n,m,beta)
  A_tilda = generate_A_tilda(X,y)


  for i in range(iterations):
    jacobian = Jacobian(n,m,A_tilda,Z)
    jacobianInverse = inverse(jacobian)
    F_Z = generate_F_Z(M,Z,n,m)
    Z= numpy.add(Z,numpy.matmul(jacobianInverse,numpy.subtract(d,F_Z)))
    Z = convertToPositive(M,Z,n,m)
    if(printbool):
      print("newton's method iteration: ",i, " and d-F(Z_0) is:")
      print(numpy.subtract(d,F_Z))
  return Z

def newton2(iterations,X,y,Z0,beta,C,tolerance,printbool):
  Z = Z0
  A = generate_A(X,y)
  n = len(A)
  m = len(A[0])-2*n
  M = generate_M(A)
  d = generate_d(n,m,beta)
  A_tilda = generate_A_tilda(X,y)
  Z = convertToPositive(M,Z,n,m)

  i= 0
  gap = calculate_gap(Z,d,n,m,M)#changes
  while(gap>tolerance and i<20):
    jacobian = Jacobian(n,m,A_tilda,Z)
    jacobianInverse = inverse(jacobian)
    F_Z = generate_F_Z(M,Z,n,m)
    Z= numpy.add(Z,numpy.matmul(jacobianInverse,numpy.subtract(d,F_Z)))
    Z = convertToPositive(M,Z,n,m) 
    if(printbool):
      print("newton's method iteration: ",i, " and d-F(Z_0) is:", numpy.subtract(d,F_Z))
    gap = calculate_gap(Z,d,n,m,M) # changes
    i= i+1 

  
  alpha = 1

  d_beta_zero = generate_d(n,m,0) # make sure the initial gap is calculated in terms of d_beta_zero
  gap = calculate_gap(Z,d_beta_zero,n,m,M) # changes

  while(gap>tolerance and i<100):
    if(printbool):
      print("newton's method iteration: ",i, " and d-F(Z_0) is:", numpy.subtract(d,F_Z))
    beta = beta*0.9
    jacobian = Jacobian(n,m,A_tilda,Z)
    jacobianInverse = inverse(jacobian)
    F_Z = generate_F_Z(M,Z,n,m)
    #alpha =1 ????
    while(check_xi_si_negative(Z,n,m) == False):
      #do the multiplication with alpha
      right_hand = numpy.subtract(d_beta_zero,F_Z)*alpha
      Z= numpy.add(Z,numpy.matmul(jacobianInverse,right_hand)) #Z_new = Z_old + jacobianInverse(d-F(Z_old))
      alpha = alpha*0.9
    #Z = convertToPositive(M,Z,n,m) # you do not need to do this because you should always be in the feasability region
    gap = calculate_gap(Z,d_beta_zero,n,m,M)
    i=i+1
  return Z

def getSVM(X,y,printbool,iterations,beta,C,tolerance,method):  
  #The order of Z is W,b,X_i,S_i,Gamma_i,Mu_i, Lambda_i

  n = len(X)
  m = len(X[0])

  A_tilda= generate_A_tilda(X,y)
  A = generate_A(X,y)
  M = generate_M(A)
  op = get_output_vector(n,m)
  Z_0=get_Z_0(M,op)
  if(printbool):
    print(numpy.matmul(M,Z_0))


  Z_0=convertToPositive(M,Z_0,n,m) 

  if(printbool):
    print("Z_0 - changed")
    print(Z_0)
    print_Xi_Si(Z_0,n,m)
    print()
    print("M(Z_0_changed) is ")
    print(numpy.matmul(M,Z_0))


  #index of xi in Z_0 is m to m+n-1
  #index of Si in Z_0 is m+n to m+2n-1
  #find index of Z_0 in between (m, m+2n-1), where Z_0[index] are negative 


  if(method ==1):
    Z = newton(iterations,X,y,Z_0,beta,C,printbool)
  if(method == 2):
    Z = newton2(iterations,X,y,Z_0,beta,C,tolerance,printbool)


  if(printbool):
    
    print("*************************")
    print("Z is: ")
    print(Z)
  print_Xi_Si(Z,n,m)
  return Z

def plotSVM_Y(X,y,Z):

  x1 = [] 
  x2 = [] 
  color = []

  for i,point in enumerate(X):
    x1.append([point[0]])
    x2.append([point[1]])
    if(y[i]==1):
      color.append(['r'])
    else:
      color.append(['b'])

  for i,j,k in zip(x1,x2,color):
    plt.scatter(i,j,s=100,marker='x',color=k,linewidths=3)

  #w1x1+w2x2+b = 0
  #w1x1= -w2x2-b
  #x1 = (-w2x2-b)/w1
  def getX1Cord(w1,w2,b,x2Cord):
    x1Cord = []
    for i in range(len(x2Cord)):
      x1Cord.append((-w2*x2Cord[i]-b)/w1)
    return x1Cord



  w1 = Z[0]
  w2 = Z[1]
  b = Z[2]
  x2Cord= [0,200]
  x1Cord = getX1Cord(w1,w2,b,x2Cord)
  plt.plot(x1Cord, x2Cord)
  plt.show()

def plotSVM(X,y,Z): #scale X

  x1 = [] 
  x2 = [] 
  color = []

  for i,point in enumerate(X):
    x1.append([point[0]])
    x2.append([point[1]])
    if(y[i]==1):
      color.append(['r'])
    else:
      color.append(['b'])

  for i,j,k in zip(x1,x2,color):
    plt.scatter(i,j,s=100,marker='x',color=k,linewidths=3)

  #w1x1+w2x2+b = 0
  #w2x2= -w1x1-b
  #x2 = (-w1x1-b)/w2
  def getX2Cord(w1,w2,b,x1Cord):
    x2Cord = []
    for i in range(len(x1Cord)):
      x2Cord.append((-w1*x1Cord[i]-b)/w2)
    return x2Cord



  w1 = Z[0]
  w2 = Z[1]
  b = Z[2]
  x1Cord= [-2,0,2]
  x2Cord = getX2Cord(w1,w2,b,x1Cord)
  plt.plot(x1Cord, x2Cord)
  plt.show()

uploaded = files.upload()

file1 = open("fourclass.txt","r") 
lines = file1.readlines()
inputX = []
y = []
for i in range(7):
  classification = (int)(lines[i].split()[0])
  coordinate1 = (float)(lines[i].split()[1][2:])
  coordinate2 = (float)(lines[i].split()[2][2:])
  inputX.append([coordinate1,coordinate2])
  y.append(classification)

inputX = numpy.array(inputX)
y = numpy.array(y)
X = generate_X(inputX)
iterations = 100
beta = 20
C= 1
tolerance = 0.1
method =2
Z = getSVM(X,y,False,iterations,beta,C,tolerance, method)
plotSVM_Y(X,y,Z)