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
