1 import tt 2 import numpy as np 3 import random 4 import time 5 st = time.time() 6 7 def wTx(alpha, trainx, trainy): return np.dot(trainx, np.dot(trainx.T, alpha * trainy)) 8 9 def EI(svmst, i): 10 wTtrainx = wTx(svmst.alpha, svmst.trainx, svmst.trainy) 11 return wTtrainx[i] - svmst.trainy[i] 12 13 def UI(svmst, i): 14 wTtrainx = wTx(svmst.alpha, svmst.trainx, svmst.trainy) 15 return float(svmst.trainy[i] * (wTtrainx[i]+ svmst.b) - 1) 16 17 def randomchoosej(i, m): 18 j = i 19 while(j==i): 20 j = random.randint(0, m-1) 21 return j 22 23 def deltacal(v1, v2): return 2.0*np.dot(v1, v2)-np.dot(v1, v1)-np.dot(v2,v2) 24 25 def LH(y1, y2, a1, a2, C): 26 if y1 == y2: L, H = max(0, a1+a2-C), min(a1+a2, C) 27 else: L, H = max(0, a2-a1), min(C+a2-a1, C) 28 return L,H 29 30 def newalpha(L,H,a): 31 if a>H: a = H 32 elif a<L: a = L 33 return a 34 35 def svm(trainx, trainy, C=0.6, toler=0.001, Maxiters=40): 36 st1 = time.time() 37 print "step 2: training..." # step 2: training... 38 iter = 0; p=0; m, n = trainx.shape; alpha = np.zeros((m, 1)); b = 0.0; wTtrainx = wTx(alpha, trainx, trainy) # m*1 wTxi 39 while (iter < Maxiters): 40 alphachange = 0 41 for i in xrange(m): 42 ui = float(trainy[i] * (wTtrainx[i]+ b) - 1) # yi(w.Txi+b) 43 if (alpha[i] < C and ui < -toler) or (alpha[i] > 0 and ui > toler): # choose i 44 j = randomchoosej(i, m) # randon choose j 45 alphajo = alpha[j].copy(); alphaio = alpha[i].copy() 46 L, H = LH(trainy[i], trainy[j], alphaio, alphajo, C)# update alpha or quit 47 delta = deltacal(trainx[i], trainx[j]) 48 alpha[j] += trainy[j]*(wTtrainx[j]-trainy[j]-wTtrainx[i]+trainy[i])/delta; alpha[j] = newalpha(L,H,alpha[j]) 49 if L == H: print "L=H,quit."; continue 50 if delta >= 0.0: print "delta>=0, quit."; continue 51 if (abs(alpha[j] - alphajo)[0] < 0.00001):print "j not moving enough, quit." ; continue 52 alpha[i] += trainy[i] * trainy[j]*(alphajo - alpha[j]) 53 wTtrainx = wTx(alpha, trainx, trainy) # update b 54 b1 = trainy[i] - wTtrainx[i]; b2 = trainy[j] - wTtrainx[j] 55 if alpha[i]>0 and alpha[i]<C: b = b1 56 elif alpha[j]>0 and alpha[j]<C: b = b2 57 else: b = (b1+b2)/2.0 58 alphachange += 1;p+=1 59 #print "iter: %d i: %d, j: %d pairs changed %d" % (iter, i, j, alphachange) 60 if alphachange == 0: iter += 1 61 #print "iteration number: ", iter 62 print "p",p 63 for i in alpha[alpha>0]: 64 print "%.8f" % i 65 print 'b', b 66 return alpha,b,(time.time()-st1) 67 68 def suportvector(alpha, trainx, trainy): 69 return trainx[np.nonzero(alpha>0)], trainy[np.nonzero(alpha>0)] 70 71 class svmsture(object): 72 def __init__(self, trainx, trainy, C, toler): 73 self.trainx = trainx 74 self.m = len(trainx) 75 self.trainy = trainy 76 self.C = C 77 self.toler = toler 78 self.alpha = np.zeros((self.m,1)) 79 self.b = 0.0 80 self.ek = np.zeros((self.m,2)) 81 82 def choosej(svmst, i, ei): 83 svmst.ek[i] = [1,ei] 84 nzeiid = np.nonzero(svmst.ek)[0] 85 if len(nzeiid) > 1: 86 maxid = 0; maxediff = 0.0 87 for id in nzeiid: 88 if id == i: continue 89 eiid = EI(svmst, id) 90 absdiff = abs(eiid - ei) 91 if absdiff > maxediff: 92 maxid = id; maxediff = absdiff 93 else: maxid = randomchoosej(i, svmst.m) 94 return maxid 95 96 def update(svmst, i): 97 ei = EI(svmst, i) 98 svmst.ek[i] = [1, ei] 99 100 def svminner(svmst, i): 101 ei = EI(svmst, i); ui = UI(svmst, i); wTtrainx = wTx(svmst.alpha, svmst.trainx, svmst.trainy) 102 if (svmst.alpha[i] < svmst.C and ui < -svmst.toler) or (svmst.alpha[i] > 0 and ui > svmst.toler): # choose i 103 j = choosej(svmst, i, ei) 104 alphajo = svmst.alpha[j].copy(); alphaio = svmst.alpha[i].copy() 105 L, H = LH(svmst.trainy[i], svmst.trainy[j], alphaio, alphajo, svmst.C) 106 delta = deltacal(svmst.trainx[i], svmst.trainx[j]) 107 if L == H: print "L=H,quit."; return 0 108 if delta >= 0.0: print "delta>=0, quit."; return 0 109 svmst.alpha[j] += svmst.trainy[j]*(wTtrainx[j]-svmst.trainy[j]-wTtrainx[i]+svmst.trainy[i])/delta; 110 svmst.alpha[j] = newalpha(L,H,svmst.alpha[j]); update(svmst, j) 111 if (abs(svmst.alpha[j] - alphajo)[0] < 0.00001):print "j not moving enough, quit." ; return 0 112 svmst.alpha[i] += svmst.trainy[i] * svmst.trainy[j]*(alphajo - svmst.alpha[j]) 113 update(svmst, i); wTtrainx = wTx(svmst.alpha, svmst.trainx, svmst.trainy) 114 b1 = svmst.trainy[i] - wTtrainx[i]; b2 = svmst.trainy[j] - wTtrainx[j] 115 if svmst.alpha[i]>0 and svmst.alpha[i]<svmst.C: svmst.b = b1 116 elif svmst.alpha[j]>0 and svmst.alpha[j]<svmst.C: svmst.b = b2 117 else: svmst.b = (b1+b2)/2.0 118 return 1 119 else: return 0 120 121 def svmouter(trainx, trainy, C = 0.6, toler = 0.001, Maxiters = 40): 122 svmst = svmsture(trainx, trainy, C, toler) 123 alphachange = 0; iter = 0 124 while (alphachange == 0 and iter < Maxiters): 125 for i in xrange(svmst.m): 126 alphachange += svminner(svmst, i) 127 print "FullSet, iter: %d, pairs changed %d" % (iter, alphachange) 128 iter += 1; alphachange = 1 129 while (alphachange > 0): 130 nonboundid = np.nonzero((svmst.alpha>0) * (svmst.alpha<svmst.C))[0]; alphachange = 0 ### 131 for i in nonboundid: 132 alphachange += svminner(svmst, i) 133 print "non-bound, iter: %d, pairs changed %d" % (iter, alphachange) 134 iter += 1 135 print "iteration number: ", iter 136 return svmst.b, svmst.alpha 137 138 def svmclassify(testx): 139 trainx, trainy = tt.loaddata('testSet') 140 b, alpha = svmouter(trainx, trainy) 141 wTtrainx = wTx(alpha, trainx, trainy) 142 testy = np.dot(testx, np.dot(trainx.T, alpha * trainy)) + b 143 print 'b',b 144 print 'alpha', alpha[alpha>0] 145 print 'w',np.dot(trainx.T, alpha * trainy) 146 print 'testy=wTtest + b',testy 147 148 trainx, trainy = tt.loaddata('testSet') 149 i=0 150 svmclassify(trainx[i]) 151 print 'trainy[%d] %d' % (i, trainy[i]) 152 """ 153 svmst = svmsture(trainx, trainy, 0.6, 0.001) 154 for i in xrange(svmst.m): 155 print i 156 print "svminner",svminner(svmst, i) 157 print "type(svm)",type(svminner(svmst, i)) 158 """ 159 160 161 162 print "cost time: ", time.time()-st