启发式规则找alphas SVM（二）

# coding=utf-8
from numpy import *
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
#加载数据集
def loadDataSet(fileName):
    dataMat=[]
    labelMat=[]
    fr=open(fileName)
    for line in fr.readlines():
        lineArr=line.strip().split('\t')
        dataMat.append([float(lineArr[0]),float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat
#i为alpha_1在样本集合中的下标 ，需要在找一个不同的alpha_2 其下表不能为i
def selectJrand(i,m):
    j=i
    while(j==i):
        j=int(random.uniform(0,m))
    return j
#alpha需要满足一定的约束条件   0<=alpha<=C sum(alpha*y)=0
def clipAlpha(aj,H,L):
    if aj>H:
        aj=H
    if L>aj:
        aj=L
    return aj
def calcWs(alphas,dataArr,classLabels):
    X=mat(dataArr);labelMat=mat(classLabels).transpose()
    m,n=shape(X)
    w=zeros((n,1))
    for i in range(m):
        w+=multiply(alphas[i]*labelMat[i],X[i,:].T)#对应的元素相乘
    return w
class optStruct:
    def __init__(self,dataMatIn,classLabels,C,toler):
        self.X=dataMatIn#输入矩阵 样本矩阵
        self.labelMat=classLabels#样本标签
        self.C=C
        self.tol=toler
        self.m=shape(dataMatIn)[0]#样本容量
        self.alphas=mat(zeros((self.m,1)))#生成一个全是0的alphas向量
        self.b=0
        self.eCache=mat(zeros((self.m,2)))#误差缓存 eCache第一列给出eCache是否有效的标志位，第二列给出的是实际值
def calcEk(oS,k):#计算给定alpha的时候的误差 k用来指定计算哪个样本的误差
    fXk=float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T))+oS.b
    Ek=fXk-float(oS.labelMat[k])
    return Ek
def selectJ(i,oS,Ei):
    #选择第二个/内循环的alpha值  选择合适的第二个alphas值 是的每次有优化的时候的步长都最大
    maxK=-1
    maxDeltaE=0
    Ej=0
    oS.eCache[i]=[1,Ei]#先设置为有效的
    validEcacheList=nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList))>1:
        for k in validEcacheList:
            if k==i:
                continue
            Ek=calcEk(oS,k)
            deltaE=abs(Ei-Ek)
            if (deltaE>maxDeltaE):
                maxK=k
                maxDeltaE=deltaE
                Ej=Ek
        return maxK,Ej
    else:#如果是第一次 就随机选择
        j=selectJrand(i,oS.m)
        Ej=calcEk(oS,j)
        return j,Ej
#更新缓存
def updateEk(oS,k):
    Ek=calcEk(oS,k)
    oS.eCache[k]=[1,Ek]

#SMO内部函数
def innerL(i,oS):
    Ei=calcEk(oS,i)#计算第一个alpha的误差
   # if ((oS.labelMat[i].Ei<-oS.tol) and (oS.alphas[i]<oS.C)) or \
    #   ((oS.labelMat[i].Ei>oS.tol) and( oS.alphas[i]>0)):
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        #SMO必须要满足KKT条件
        j,Ej=selectJ(i,oS,Ei)#选择具有最大步长的第二个alpha
        alphaIold=oS.alphas[i].copy()
        alphaJold=oS.alphas[j].copy()
        if (oS.labelMat[i]!=oS.labelMat[j]):
            L=max(0,oS.alphas[j]-oS.alphas[i])
            H=min(oS.C,oS.C+oS.alphas[j]-oS.alphas[i])
        else:
            L=max(0,oS.alphas[j]+oS.alphas[i]-oS.C)
            H=min(oS.C,oS.alphas[j]+oS.alphas[i])
        eta=2.0*oS.X[i,:]*oS.X[j,:].T-oS.X[i,:]*oS.X[i,:].T-oS.X[j,:]*oS.X[j,:].T
        oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta
        oS.alphas[j]=clipAlpha(oS.alphas[j],H,L)
        updateEk(oS,j)
        if (abs(oS.alphas[j]-alphaJold)<0.00001):
            return 0
        oS.alphas[i]+=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j])
        updateEk(oS,i)
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*(oS.X[i,:]*oS.X[i,:].T) - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*(oS.X[i,:]*oS.X[j,:].T)#oS.K[i,j]
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*(oS.X[i,:]*oS.X[j,:].T) - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*(oS.X[j,:]*oS.X[j,:].T)#oS.K[j,j]     
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0

#选择第一个alpha的策略，在所有的数据集上进行单边扫描，在非边界上进行扫描
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)#, kTup)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):# entireSet控制在边界上和整个数据集上进行选择alpha
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):        
                alphaPairsChanged += innerL(i,oS)
               # print "fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
            iter += 1
        else:#go over non-bound (railed) alphas
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                #print "non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
            iter += 1
        if entireSet: entireSet = False #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True  
        print "iteration number: %d" % iter
    return oS.b,oS.alphas,oS.C


dataArr,labelArr=loadDataSet("testSet.txt")
b,alphas,C=smoP(dataArr,labelArr,0.6,0.01,40)
ws=calcWs(alphas,dataArr,labelArr)
print ws

from matplotlib.patches import Circle
xcord0=[]
ycord0=[]
xcord1=[]
ycord1=[]
markers=[]
colors=[]
fr=open("testSet.txt")
for line in fr.readlines():
    lineSplit = line.strip().split('\t')
    xPt = float(lineSplit[0])
    yPt = float(lineSplit[1])
    label = int(lineSplit[2])
    if (label == -1):
        xcord0.append(xPt)
        ycord0.append(yPt)
    else:
        xcord1.append(xPt)
        ycord1.append(yPt)
fr.close()
xd=[]
yd=[]
m=shape(dataArr)[0]
for i in range(100):
    if alphas[i]>0.0 and alphas[i]!=C:
         xd.append(dataArr[i])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord0,ycord0, marker='s', s=50)
ax.scatter(xcord1,ycord1, marker='o', s=50, c='red')
plt.title('SVM')
totalm=shape(xd)[0];
for i in range(totalm):
    circle = Circle((xd[i][0], xd[i][1]), 0.3, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=2, alpha=0.5)
    ax.add_patch(circle)
w0=float(ws[0])
w1=float(ws[1])
x=arange(-2.0,12.0,0.1)
x=mat(x)
z=(-w0*x-b)/w1
z=mat(z)
ax.plot(x.tolist()[0],z.tolist()[0])
ax.axis([-2,12,-8,6])
plt.show()

速度特别快