Ramp loss SVM.py

class SVM(object):
    def __init__(self,n_in,c0=1,c1=1,loss = RAMP_LOSS,gamma = 0.5):
        self.dim = n_in
        self.w = np.zeros(n_in)
        self.b = np.zeros(1)
        self.c0 = c0
        self.c1 = c1
        self.gamma = gamma
        self.x = None
        self.y = None

    def train(self,x,y):
        #setup
        self.x = x
        self.y = y
        x_2 = np.zeros(len(y))
        for i in range(self.dim):
            x_2 += x[:,i]**2
        x_2 = np.tile(x_2,(len(y),1))
        yy = np.tile(y,(len(y),1))
        yy = yy * yy.T
        K = np.exp(-self.gamma*(x_2 + x_2.T - 2*x.dot(x.T))) * yy
        '''
        conpute a^t by solving the following convex problem
        max(a - 1/2*a^T*K*a)
        subject to y*a = 0, -b^(t-1)<=a<=C-b^(t-1)
        '''
        self.w = np.ones(len(x))
        s = -1
        c = np.array([self.c0 if y_i == 1 else self.c1 for y_i in y])
        b_old = np.zeros(len(y))
        b_new = np.zeros(len(y))

        
        def f(x_l):
            return (y*np.asarray(alpha)).dot(np.exp((-self.gamma*(self.x-x_l)**2).sum(1)))+self.b
        

        while True:
            #set variables and constraint -b^(t-1)<=a<=C-b^(t-1)
            m = gurobi.Model("qp")
            a = {}
            for i in range(len(y)):
                a[i] = m.addVar(lb=-b_old[i],ub=c[i]-b_old[i])
                m.update()
                
            #set objective
            Ka = {}
            aKa = gurobi.QuadExpr()
            for i in range(len(y)):
                Ka[i] = gurobi.LinExpr()
                for j in range(len(y)):
                    Ka[i].add(a[i],K[i,j])
            
            for i in range(len(y)):
                aKa.add(a[i] * Ka[i])
                
            objective = gurobi.quicksum(a) - 1.0/2*aKa
            m.setObjective(objective,gurobi.GRB.MAXIMIZE)
            m.update()
            
            #add constraints
            ya = gurobi.LinExpr()
            for i in range(len(y)):
                ya += a[i]*y[i]

            m.addConstr(ya == 0,"c0")
            m.update()

            #m.params.logToConsole = 0
            m.optimize()
            alpha = {}
            i = 0
            for v in m.getVars():
                alpha[i] = v.x
                i += 1
            alpha = alpha.values()

            #compute b
            new_b = 0.
            count = 0
            for i in range(len(alpha)):
                if (0 < alpha[i]) & (alpha[i] < c[i]):
                    new_b +=  y[i]-(f(x[i])-self.b)
                    count += 1
            self.b = new_b/count
            print self.b
                
            #compute beta
            b_new = np.array([c[i] if y[i]*f(x[i]) < s else 0 for i in range(len(y))])
            if (b_old == b_new).all():
                break
            b_old = b_new
        self.w = y*np.array(alpha)#np.array(-np.asarray([1 if y[i]*f(x[i]) <= 0 else 0 for i in range(len(y))])*c*y)
        #print self.w,self.b

    def predict(self,x):
        return self.w.dot(np.exp((-self.gamma*(self.x-x)**2).sum(1)))+self.b