1 # encoding: utf-8
  2 import numpy as np
  3 import matplotlib.pyplot as plt
  4 import cPickle
  5 import gzip
  6
  7 class SVC(object):
  8     def __init__(self, c=1.0, delta=0.001):  # 初始化
  9         self.N = 0
 10         self.delta = delta
 11         self.X = None
 12         self.y = None
 13         self.w = None
 14         self.wn = 0
 15         self.K = np.zeros((self.N, self.N))
 16         self.a = np.zeros((self.N, 1))
 17         self.b = 0
 18         self.C = c
 19         self.stop=1
 20         self.k=0
 21         self.cls=0
 22         self.train_result=[]
 23
 24     def kernel_function(self,x1, x2):  # 核函数
 25         return np.dot(x1, x2)
 26
 27     def kernel_matrix(self, x):  # 核矩阵
 28         for i in range(0, len(x)):
 29             for j in range(i, len(x)):
 30                 self.K[j][i] = self.K[i][j] = self.kernel_function(self.X[i], self.X[j])
 31
 32     def get_w(self):  # 计算更新w
 33         ay = self.a * self.y
 34         w = np.zeros((1, self.wn))
 35         for i in range(0, self.N):
 36             w += self.X[i] * ay[i]
 37         return w
 38
 39     def get_b(self, a1, a2, a1_old, a2_old):  # 计算更新B
 40         y1 = self.y[a1]
 41         y2 = self.y[a2]
 42         a1_new = self.a[a1]
 43         a2_new = self.a[a2]
 44         b1_new = -self.E[a1] - y1 * self.K[a1][a1] * (a1_new - a1_old) - y2 * self.K[a2][a1] * (
 45             a2_new - a2_old) + self.b
 46         b2_new = -self.E[a2] - y1 * self.K[a1][a2] * (a1_new - a1_old) - y2 * self.K[a2][a2] * (
 47             a2_new - a2_old) + self.b
 48         if (0 < a1_new) and (a1_new < self.C) and (0 < a2_new) and (a2_new < self.C):
 49             return b1_new[0]
 50         else:
 51             return (b1_new[0] + b2_new[0]) / 2.0
 52
 53     def gx(self, x):  # 判别函数g(x)
 54         return np.dot(self.w, x) + self.b
 55
 56     def satisfy_kkt(self, a):  # 判断样本点是否满足kkt条件
 57         index = a[1]
 58         if a[0] == 0 and self.y[index] * self.gx(self.X[index]) > 1:
 59             return 1
 60         elif a[0] < self.C and self.y[index] * self.gx(self.X[index]) == 1:
 61             return 1
 62         elif a[0] == self.C and self.y[index] * self.gx(self.X[index]) < 1:
 63             return 1
 64         return 0
 65
 66     def clip_func(self, a_new, a1_old, a2_old, y1, y2):  # 拉格朗日乘子的裁剪函数
 67         if (y1 == y2):
 68             L = max(0, a1_old + a2_old - self.C)
 69             H = min(self.C, a1_old + a2_old)
 70         else:
 71             L = max(0, a2_old - a1_old)
 72             H = min(self.C, self.C + a2_old - a1_old)
 73         if a_new < L:
 74             a_new = L
 75         if a_new > H:
 76             a_new = H
 77         return a_new
 78
 79     def update_a(self, a1, a2):  # 更新a1,a2
 80         partial_a2 = self.K[a1][a1] + self.K[a2][a2] - 2 * self.K[a1][a2]
 81         if partial_a2 <= 1e-9:
 82             print "error:", partial_a2
 83         a2_new_unc = self.a[a2] + (self.y[a2] * ((self.E[a1] - self.E[a2]) / partial_a2))
 84         a2_new = self.clip_func(a2_new_unc, self.a[a1], self.a[a2], self.y[a1], self.y[a2])
 85         a1_new = self.a[a1] + self.y[a1] * self.y[a2] * (self.a[a2] - a2_new)
 86         if abs(a1_new - self.a[a1]) < self.delta:
 87             return 0
 88         self.a[a1] = a1_new
 89         self.a[a2] = a2_new
 90         self.is_update = 1
 91         return 1
 92
 93     def update(self, first_a):  # 更新拉格朗日乘子
 94         for second_a in range(0, self.N):
 95             if second_a == first_a:
 96                 continue
 97             a1_old = self.a[first_a]
 98             a2_old = self.a[second_a]
 99             if self.update_a(first_a, second_a) == 0:
100                 return
101             self.b= self.get_b(first_a, second_a, a1_old, a2_old)
102             self.w = self.get_w()
103             self.E = [self.gx(self.X[i]) - self.y[i] for i in range(0, self.N)]
104             self.stop=0
105
106     def train(self, x, y, max_iternum=100):  # SMO算法
107         x_len = len(x)
108         self.X = x
109         self.N = x_len
110         self.wn = len(x[0])
111         self.y = np.array(y).reshape((self.N, 1))
112         self.K = np.zeros((self.N, self.N))
113         self.kernel_matrix(self.X)
114         self.b = 0
115         self.a = np.zeros((self.N, 1))
116         self.w = self.get_w()
117         self.E = [self.gx(self.X[i]) - self.y[i] for i in range(0, self.N)]
118         self.is_update = 0
119         for i in range(0, max_iternum):
120             self.stop=1
121             data_on_bound = [[x,y] for x,y in zip(self.a, range(0, len(self.a))) if x > 0 and x< self.C]
122             if len(data_on_bound) == 0:
123                 data_on_bound = [[x,y] for x,y in zip(self.a, range(0, len(self.a)))]
124             for data in data_on_bound:
125                 if self.satisfy_kkt(data) != 1:
126                     self.update(data[1])
127             if self.is_update == 0:
128                 for data in [[x,y] for x,y in zip(self.a, range(0, len(self.a)))]:
129                     if self.satisfy_kkt(data) != 1:
130                         self.update(data[1])
131             if self.stop:
132                 break
133         return self.w, self.b
134
135     def fit(self,x, y):  # 训练模型, 一对一法k(k-1)/2个SVM进行多类分类
136         self.cls, y = np.unique(y, return_inverse=True)
137         self.k=len(self.cls)
138         for i in range(self.k):
139             for j in range(i):
140                 a,b=self.sub_data(x,y,i,j)
141                 self.train_result.append([i,j,self.train(a,b)])
142
143     def predict(self,x_new):  # 预测
144          p=np.zeros(self.k)
145          for i,j,w in self.train_result:
146              self.w=w[0]
147              self.b=w[1]
148              if self.classfy(x_new)==1:
149                  p[j]+=1
150              else:
151                  p[i]+=1
152          return self.cls[np.argmax(p)]
153
154     def sub_data(self,x,y,i,j):  # 数据分类
155         subx=[]
156         suby=[]
157         for a,b in zip(x,y):
158             if b==i:
159                  subx.append(a)
160                  suby.append(-1)
161             elif b==j:
162                  subx.append(a)
163                  suby.append(1)
164         return subx,suby
165
166     def classfy(self,x_new):  # 预测
167         y_new=self.gx(x_new)
168         cl = int(np.sign(y_new))
169         if cl == 0:
170             cl = 1
171         return cl
172
173
174 def load_data():
175     f = gzip.open('../data/mnist.pkl.gz', 'rb')
176     training_data, validation_data, test_data = cPickle.load(f)
177     f.close()
178     return (training_data, validation_data, test_data)
179
180 if __name__ == "__main__":
181     svc = SVC()
182     np.random.seed(0)
183     l=1000
184     training_data, validation_data, test_data = load_data()
185     svc.fit(training_data[0][:l],training_data[1][:l])
186     predictions = [svc.predict(a) for a in test_data[0][:l]]
187     num_correct = sum(int(a == y) for a, y in zip(predictions, test_data[1][:l]))
188     print "%s of %s values correct." % (num_correct, len(test_data[1][:l]))  #72/100  #808/1000  #8194/10000(较慢）