机器学习算法实战项目—支持向量机（2）—完整版的SMO算法

2.完整版的SMO算法

在几百个点组成的小规模数据集上，简化版SMO算法的运行是没有什么问题的，但是在更大的数据集上的运行速度就会变慢。

刚才已经讨论了简化版SMO算法，下面我们就讨论完整版的Platt SMO算法。

在这两个版本中，实现alpha的更改和代数运算的优化环节一模一样。在优化过程中，唯一的不同就是选择alpha的方式。

完整版的Platt SMO算法应用了一些能够提速的启发方法。上一节的例子在执行时存在一定的时间提升空间。

Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间进行交替：一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单遍扫描。

而所谓非边界alpha指的就是那些不等于边界0或C的alpha值。

对整个数据集的扫描相当容易，而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的alpha值。

在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha值。在优化过程中，会通过最大化步长的方式来获得第二个alpha值。

在简化版SMO算法中，我们会在选择j之后计算错误率Ej。但在这里，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的alpha值。

辅助函数和上面的步骤一样，即下面的三个函数:

"""
函数说明：读取数据
Parameters:
    file_name - 文件名
Returns:
    data_mat - 数据矩阵
    label_mat - 数据标签
"""
def load_dataset(file_name):
    # 数据矩阵
    data_mat = []
    # 标签向量
    label_mat = []
    # 打开文件
    fr = open(file_name)
    # 逐行读取
    for line in fr.readlines():
        # 去掉每一行首尾的空白符，例如'\n','\r','\t',' '
        # 将每一行内容根据'\t'符进行切片
        line_array = line.strip().split('\t')
        # 添加数据(100个元素排成一行)
        data_mat.append([float(line_array[0]), float(line_array[1])])
        # 添加标签(100个元素排成一行)
        label_mat.append(float(line_array[2]))
    return data_mat, label_mat

"""
函数说明：随机选择alpha_j
Parameters:
    i - alpha_i的索引值
    m - alpha参数个数
Returns:
    j - alpha_j的索引值
"""
def select_j_random(i, m):
    j = i
    while(j == i):
        # uniform()方法将随机生成一个实数，它在[x, y)范围内
        j = int(random.uniform(0, m))
    return j

"""
函数说明：修剪alpha_j
Parameters:
    aj - alpha_j值
    H - alpha上限
    L - alpha下限 
Returns:
    aj - alpha_j值
"""
def clip_alpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

我们首先定义一个数据结构:

"""
类说明：维护所有需要操作的值
Parameters:
    dataMatIn - 数据矩阵
    classLabels - 数据标签
    C - 松弛变量
    toler - 容错率   
Returns:
    None
"""
"""
定义一个新的数据结构
"""
class opt_struct:
    def __init__(self, data_mat_in, class_labels, C, toler):
        # 数据矩阵
        self.X = data_mat_in #传进来的数据
        # 数据标签
        self.label_mat = class_labels
        # 松弛变量
        self.C = C
        # 容错率
        self.tol = toler
        # 矩阵的行数
        self.m = np.shape(data_mat_in)[0]
        # 根据矩阵行数初始化alphas矩阵，一个m行1列的全零列向量
        self.alphas = np.mat(np.zeros((self.m, 1)))
        # 初始化b参数为0
        self.b = 0
        # 根据矩阵行数初始化误差缓存矩阵，第一列为是否有效标志位，其中0无效，1有效;第二列为实际的误差Ei的值
        """
        我们之前的定义为:Ei=gxi-yi
        yi是标签的实际值。
        gx=alpha_i*y_i*x_i.x,就相当于是w.x+b
        因为误差值经常用到，所以希望每次计算后放到一个缓存当中，将ecache一分为二，第一列是标志位，取值为0或者1，为1时表示已经算出来
        """
        self.ecache = np.mat(np.zeros((self.m, 2)))

首要的事情就是建立一个数据结构来保存所有的重要值，而这个过程可以通过一个对象来完成。只是作为一个数据结构来使用对象。

在将值传给函数时，我们可以通过将所有数据移到一个结构中来实现，这样就可以省掉手工输入的麻烦了。而此时，数据就可以通过一个对象来进行传递。实际上，当完成其实现时，可以很容易通过Python的字典来完成。

该方法可以实现其成员变量的填充。除了增加了一个m×2的矩阵成
员变量ecache之外，这些做法和简化版SMO一模一样。ecache的第一列给出的是ecache是否有效的标志位，而第二列给出的是实际的E值。

对于给定的alpha值，第一个辅助函数cal-Ek()能够计算E值并返回。以前，该过程是采用内嵌的方式来完成的，但是由于该过程在这个版本的SMO算法中出现频繁，这里必须要将其单独拎出来。

"""
函数说明：计算误差
Parameters:
    os - 数据结构
    k - 标号为k的数据
Returns:
    Ek - 标号为k的数据误差
"""
def cal_Ek(os, k):
    # multiply(a,b)就是个乘法，如果a,b是两个数组，那么对应元素相乘
    # .T为转置
    fXk = float(np.multiply(os.alphas, os.label_mat).T * (os.X * os.X[k, :].T) + os.b)
    # 计算误差项
    Ek = fXk - float(os.label_mat[k])
    # 返回误差项
    return Ek

函数select-j()用于选择第二个alpha或者说内循环的alpha值。回想一下，这里的目标是选择合适的第二个alpha值以保证在每次优化中采用最大步长。该函数的误差值与第一个alpha值Ei和下标i有关。

首先将输入值Ei在缓存中设置成为有效的。这里的有效意味着它已经计算好了。在ecache中，代码nonzero(os.ecache[:,0].A)[0]构建出了一个非零表。

NumPy函数nonzero()返回了一个列表，而这个列表中包含以输入列表为目录的列表值，当然这里的值并非零。nonzero()语句返回的是非零E值所对应的alpha值，而不是E值本身。

程序会在所有的值上进行循环并选择其中使得改变最大的那个值。如果这是第一次循环的话，那么就随机选择一个alpha值。当然也存在有许多更复杂的方式来处理第一次循环的情况，而上述做法就能够满足我们的目的。

"""
函数说明：内循环启发方式2
选择第二个待优化的alpha_j，选择一个误差最大的alpha_j
即，我们在选择alpha_2的时候做了改进，选择误差最大的
Parameters:
    i - 标号为i的数据的索引值
    oS - 数据结构
    Ei - 标号为i的数据误差 
Returns:
    j - 标号为j的数据的索引值
    maxK - 标号为maxK的数据的索引值
    Ej - 标号为j的数据误差
"""
def select_j(i, os, Ei):
    # 初始化
    max_K = -1 #下标的索引值
    max_delta_E = 0
    Ej = 0
    # 根据Ei更新误差缓存，即先计算alpha_1以及E1值
    os.ecache[i] = [1, Ei] #放入缓存当中，设为有效
    # 对一个矩阵.A转换为Array类型
    # 返回误差不为0的数据的索引值
    valid_ecache_list = np.nonzero(os.ecache[:, 0].A)[0] #找出缓存中不为0
    # 有不为0的误差
    if(len(valid_ecache_list) > 1):
        # 遍历，找到最大的Ek
        for k in valid_ecache_list: #迭代所有有效的缓存，找到误差最大的E
            # 不计算k==i节省时间
            if k == i: #不选择和i相等的值
                continue
            # 计算Ek
            Ek = cal_Ek(os, k)
            # 计算|Ei - Ek|
            delta_E = abs(Ei - Ek)
            # 找到maxDeltaE
            if(delta_E > max_delta_E):
                max_K = k
                max_delta_E = delta_E
                Ej = Ek
        # 返回max_K，Ej
        return max_K, Ej #这样我们就得到误差最大的索引值和误差最大的值
    # 没有不为0的误差
    else: #第一次循环时是没有有效的缓存值的，所以随机选一个（仅会执行一次）
        # 随机选择alpha_j的索引值
        j = select_j_random(i, os.m)
        # 计算Ej
        Ej = cal_Ek(os, j)
    # 返回j，Ej
    return j, Ej

最后一个辅助函数是update-Ek()，它会计算误差值并存入缓存当中。在对alpha值进行优化之后会用到这个值。

"""
函数说明：计算Ek,并更新误差缓存
Parameters:
    os - 数据结构
    k - 标号为k的数据的索引值
Returns:
    None
"""
def update_Ek(os, k):
    # 计算Ek
    Ek = cal_Ek(os, k)
    # 更新误差缓存
    os.ecache[k] = [1, Ek]

接下来将简单介绍一下用于寻找决策边界的优化例程:

"""
函数说明：优化的SMO算法
Parameters:
    i - 标号为i的数据的索引值
    os - 数据结构
Returns:
    1 - 有任意一对alpha值发生变化
    0 - 没有任意一对alpha值发生变化或变化太小
"""
def innerL(i, os):
    # 步骤1：计算误差Ei
    Ei = cal_Ek(os, i)
    # 优化alpha,设定一定的容错率
    if((os.label_mat[i] * Ei < -os.tol) and (os.alphas[i] < os.C)) or ((os.label_mat[i] * Ei > os.tol) and (os.alphas[i] > 0)):
        # 使用内循环启发方式2选择alpha_j,并计算Ej
        j, Ej = select_j(i, os, Ei) #这里不再是随机选取了
        # 保存更新前的alpha值，使用深层拷贝
        alpha_i_old = os.alphas[i].copy()
        alpha_j_old = os.alphas[j].copy()
        # 步骤2：计算上界H和下界L
        if(os.label_mat[i] != os.label_mat[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])
        if L == H:
            print("L == H")
            return 0
        # 步骤3：计算eta
        eta = 2.0 * os.X[i, :] * os.X[j, :].T - os.X[i, :] * os.X[i, :].T - os.X[j, :] * os.X[j, :].T
        if eta >= 0:
            print("eta >= 0")
            return 0
        # 步骤4：更新alpha_j
        os.alphas[j] -= os.label_mat[j] * (Ei - Ej) / eta
        # 步骤5：修剪alpha_j
        os.alphas[j] = clip_alpha(os.alphas[j], H, L)
        # 更新Ej至误差缓存
        update_Ek(os, j)
        if(abs(os.alphas[j] - alpha_j_old) < 0.00001):
            print("alpha_j变化太小")
            return 0
        # 步骤6：更新alpha_i
        os.alphas[i] += os.label_mat[i] * os.label_mat[j] * (alpha_j_old - os.alphas[j])
        # 更新Ei至误差缓存
        update_Ek(os, i)
        # 步骤7：更新b_1和b_2:
        b1 = os.b - Ei - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[i, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[i, :].T
        b2 = os.b - Ej - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[j, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[j, :].T
        # 步骤8：根据b_1和b_2更新b
        if(0 < os.alphas[i] < os.C):
            os.b = b1
        elif(0 < os.alphas[j] < os.C):
            os.b = b2
        else:
            os.b = (b1 + b2) / 2.0
        return 1 #表示有更新
    else:
        return 0 #表示没有更新

代码几乎与上一节中给出的smo-simple函数一摸一样，但是这里的代码已经使用了自己的数据结构。该结构在参数os中传递。

完整版SMO算法的外循环代码:

"""
函数说明：完整的线性SMO算法
Parameters:
    dataMatIn - 数据矩阵
    classLabels - 数据标签
    C - 松弛变量
    toler - 容错率
    maxIter - 最大迭代次数
Returns:
    oS.b - SMO算法计算的b
    oS.alphas - SMO算法计算的alphas
"""
def smo_p(data_mat_in, class_labels, C, toler, max_iter):
    # 初始化数据结构
    os = opt_struct(np.mat(data_mat_in), np.mat(class_labels).transpose(), C, toler)
    # 初始化当前迭代次数
    iter = 0
    entrie_set = True #是否在全部数据集上迭代
    alpha_pairs_changed = 0
    # 遍历整个数据集alpha都没有更新或者超过最大迭代次数，则退出循环
    while(iter < max_iter) and ((alpha_pairs_changed > 0) or (entrie_set)):
        alpha_pairs_changed = 0
        if entrie_set: # 遍历整个数据集
            for i in range(os.m):
                # 使用优化的SMO算法
                alpha_pairs_changed += innerL(i, os) #innerL返回的值是0或者1
                print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
            iter += 1
        # 遍历非边界值
        else:
            # 遍历不在边界0和C的alpha
            non_bound_i_s = np.nonzero((os.alphas.A > 0) * (os.alphas.A < C))[0]
            for i in non_bound_i_s:
                alpha_pairs_changed += innerL(i, os)
                print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
            iter += 1
        # 遍历一次后改为非边界遍历
        if entrie_set:
            entrie_set = False #进行切换，遍历非边界数据集
        # 如果alpha没有更新，计算全样本遍历
        elif(alpha_pairs_changed == 0):
            entrie_set = True
        print("迭代次数:%d" % iter)
    # 返回SMO算法计算的b和alphas
    return os.b, os.alphas

其输入和函数smo-simple()完全一样。函数一开始构建一个数据结构来容纳所有的数据，然后需要对控制函数退出的一些变量进行初始化。整个代码的主体是while循环，这与smo-simple()有些类似，但是这里的循环退出条件更多一些。

当迭代次数超过指定的最大值，或者遍历整个集合都未对任意alpha对进行修改时，就退出循环。

这里的max-iter变量和函数smo-simple()中的作用有一点不同，后者当没有任何alpha发生改变时会将整个集合的一次遍历过程计成一次迭代，而这里的一次迭代定义为一次循环过程，而不管该循环具体做了什么事。此时，如果在优化过程中存在波动就会停止，因此这里的做法优于smo-simple()函数中的计数方法。

while循环的内部与smo-simple()中有所不同，一开始的for循环在数据集上遍历任意可能的alpha。

我们通过调用innerL()来选择第二个alpha，并在可能时对其进行优化处理。如果有任意一对alpha值发生改变，那么会返回
1。第二个for循环遍历所有的非边界alpha值，也就是不在边界0或C上的值。

接下来，我们对for循环在非边界循环和完整遍历之间进行切换，并打印出迭代次数。最后程序将会返回常数b和alpha值。

我们的部分结果如下:

全样本遍历:第0次迭代 样本:0, alpha优化次数:1
全样本遍历:第0次迭代 样本:1, alpha优化次数:1
全样本遍历:第0次迭代 样本:2, alpha优化次数:1
全样本遍历:第0次迭代 样本:3, alpha优化次数:1
全样本遍历:第0次迭代 样本:4, alpha优化次数:1
全样本遍历:第0次迭代 样本:5, alpha优化次数:1
全样本遍历:第0次迭代 样本:6, alpha优化次数:1
全样本遍历:第0次迭代 样本:7, alpha优化次数:1
全样本遍历:第0次迭代 样本:8, alpha优化次数:2
全样本遍历:第0次迭代 样本:9, alpha优化次数:2
......
非边界遍历:第1次迭代 样本:26, alpha优化次数:3
alpha_j变化太小
非边界遍历:第1次迭代 样本:29, alpha优化次数:3
alpha_j变化太小
非边界遍历:第1次迭代 样本:46, alpha优化次数:3
alpha_j变化太小
非边界遍历:第1次迭代 样本:54, alpha优化次数:3
非边界遍历:第1次迭代 样本:55, alpha优化次数:3
迭代次数:2
alpha_j变化太小
非边界遍历:第2次迭代 样本:8, alpha优化次数:0
alpha_j变化太小
非边界遍历:第2次迭代 样本:10, alpha优化次数:0
非边界遍历:第2次迭代 样本:23, alpha优化次数:0
alpha_j变化太小
非边界遍历:第2次迭代 样本:29, alpha优化次数:0
alpha_j变化太小
非边界遍历:第2次迭代 样本:46, alpha优化次数:0
alpha_j变化太小
非边界遍历:第2次迭代 样本:54, alpha优化次数:0
非边界遍历:第2次迭代 样本:55, alpha优化次数:0
......
全样本遍历:第3次迭代 样本:97, alpha优化次数:0
全样本遍历:第3次迭代 样本:98, alpha优化次数:0
全样本遍历:第3次迭代 样本:99, alpha优化次数:0
迭代次数:4

我们最终的决策边界和支持向量可视化为:

b= [[-3.06127512]]
w= [[ 0.69141217]
 [-0.3411598 ]]
alphas= [[0.06158536 0.08744873 0.09509032 0.07285032 0.01603246 0.15524162]]
[3.542485, 1.977398] -1.0
[3.125951, 0.293251] -1.0
[4.658191, 3.507396] -1.0
[8.197181, 1.545132] 1.0
[5.286862, -2.358286] 1.0
[6.080573, 0.418886] 1.0

完整版代码如下:

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import random
import matplotlib as mpl
"""
改进SMO算法以加快我们的SVM运行速度！
"""

"""
函数说明：读取数据
Parameters:
    file_name - 文件名
Returns:
    data_mat - 数据矩阵
    label_mat - 数据标签
"""
def load_dataset(file_name):
    # 数据矩阵
    data_mat = []
    # 标签向量
    label_mat = []
    # 打开文件
    fr = open(file_name)
    # 逐行读取
    for line in fr.readlines():
        # 去掉每一行首尾的空白符，例如'\n','\r','\t',' '
        # 将每一行内容根据'\t'符进行切片
        line_array = line.strip().split('\t')
        # 添加数据(100个元素排成一行)
        data_mat.append([float(line_array[0]), float(line_array[1])])
        # 添加标签(100个元素排成一行)
        label_mat.append(float(line_array[2]))
    return data_mat, label_mat

"""
函数说明：随机选择alpha_j
Parameters:
    i - alpha_i的索引值
    m - alpha参数个数
Returns:
    j - alpha_j的索引值
"""
def select_j_random(i, m):
    j = i
    while(j == i):
        # uniform()方法将随机生成一个实数，它在[x, y)范围内
        j = int(random.uniform(0, m))
    return j

"""
函数说明：修剪alpha_j
Parameters:
    aj - alpha_j值
    H - alpha上限
    L - alpha下限 
Returns:
    aj - alpha_j值
"""
def clip_alpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

"""
类说明：维护所有需要操作的值
Parameters:
    dataMatIn - 数据矩阵
    classLabels - 数据标签
    C - 松弛变量
    toler - 容错率   
Returns:
    None
"""
"""
定义一个新的数据结构
"""
class opt_struct:
    def __init__(self, data_mat_in, class_labels, C, toler):
        # 数据矩阵
        self.X = data_mat_in #传进来的数据
        # 数据标签
        self.label_mat = class_labels
        # 松弛变量
        self.C = C
        # 容错率
        self.tol = toler
        # 矩阵的行数
        self.m = np.shape(data_mat_in)[0]
        # 根据矩阵行数初始化alphas矩阵，一个m行1列的全零列向量
        self.alphas = np.mat(np.zeros((self.m, 1)))
        # 初始化b参数为0
        self.b = 0
        # 根据矩阵行数初始化误差缓存矩阵，第一列为是否有效标志位，其中0无效，1有效;第二列为实际的误差Ei的值
        """
        我们之前的定义为:Ei=gxi-yi
        yi是标签的实际值。
        gx=alpha_i*y_i*x_i.x,就相当于是w.x+b
        因为误差值经常用到，所以希望每次计算后放到一个缓存当中，将ecache一分为二，第一列是标志位，取值为0或者1，为1时表示已经算出来
        """
        self.ecache = np.mat(np.zeros((self.m, 2)))

"""
函数说明：计算误差
Parameters:
    os - 数据结构
    k - 标号为k的数据
Returns:
    Ek - 标号为k的数据误差
"""
def cal_Ek(os, k):
    # multiply(a,b)就是个乘法，如果a,b是两个数组，那么对应元素相乘
    # .T为转置
    fXk = float(np.multiply(os.alphas, os.label_mat).T * (os.X * os.X[k, :].T) + os.b)
    # 计算误差项
    Ek = fXk - float(os.label_mat[k])
    # 返回误差项
    return Ek


"""
函数说明：内循环启发方式2
选择第二个待优化的alpha_j，选择一个误差最大的alpha_j
即，我们在选择alpha_2的时候做了改进，选择误差最大的
Parameters:
    i - 标号为i的数据的索引值
    oS - 数据结构
    Ei - 标号为i的数据误差 
Returns:
    j - 标号为j的数据的索引值
    maxK - 标号为maxK的数据的索引值
    Ej - 标号为j的数据误差
"""
def select_j(i, os, Ei):
    # 初始化
    max_K = -1 #下标的索引值
    max_delta_E = 0
    Ej = 0
    # 根据Ei更新误差缓存，即先计算alpha_1以及E1值
    os.ecache[i] = [1, Ei] #放入缓存当中，设为有效
    # 对一个矩阵.A转换为Array类型
    # 返回误差不为0的数据的索引值
    valid_ecache_list = np.nonzero(os.ecache[:, 0].A)[0] #找出缓存中不为0
    # 有不为0的误差
    if(len(valid_ecache_list) > 1):
        # 遍历，找到最大的Ek
        for k in valid_ecache_list: #迭代所有有效的缓存，找到误差最大的E
            # 不计算k==i节省时间
            if k == i: #不选择和i相等的值
                continue
            # 计算Ek
            Ek = cal_Ek(os, k)
            # 计算|Ei - Ek|
            delta_E = abs(Ei - Ek)
            # 找到maxDeltaE
            if(delta_E > max_delta_E):
                max_K = k
                max_delta_E = delta_E
                Ej = Ek
        # 返回max_K，Ej
        return max_K, Ej #这样我们就得到误差最大的索引值和误差最大的值
    # 没有不为0的误差
    else: #第一次循环时是没有有效的缓存值的，所以随机选一个（仅会执行一次）
        # 随机选择alpha_j的索引值
        j = select_j_random(i, os.m)
        # 计算Ej
        Ej = cal_Ek(os, j)
    # 返回j，Ej
    return j, Ej

"""
函数说明：计算Ek,并更新误差缓存
Parameters:
    os - 数据结构
    k - 标号为k的数据的索引值
Returns:
    None
"""
def update_Ek(os, k):
    # 计算Ek
    Ek = cal_Ek(os, k)
    # 更新误差缓存
    os.ecache[k] = [1, Ek]

"""
函数说明：优化的SMO算法
Parameters:
    i - 标号为i的数据的索引值
    os - 数据结构
Returns:
    1 - 有任意一对alpha值发生变化
    0 - 没有任意一对alpha值发生变化或变化太小
"""
def innerL(i, os):
    # 步骤1：计算误差Ei
    Ei = cal_Ek(os, i)
    # 优化alpha,设定一定的容错率
    if((os.label_mat[i] * Ei < -os.tol) and (os.alphas[i] < os.C)) or ((os.label_mat[i] * Ei > os.tol) and (os.alphas[i] > 0)):
        # 使用内循环启发方式2选择alpha_j,并计算Ej
        j, Ej = select_j(i, os, Ei) #这里不再是随机选取了
        # 保存更新前的alpha值，使用深层拷贝
        alpha_i_old = os.alphas[i].copy()
        alpha_j_old = os.alphas[j].copy()
        # 步骤2：计算上界H和下界L
        if(os.label_mat[i] != os.label_mat[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])
        if L == H:
            print("L == H")
            return 0
        # 步骤3：计算eta
        eta = 2.0 * os.X[i, :] * os.X[j, :].T - os.X[i, :] * os.X[i, :].T - os.X[j, :] * os.X[j, :].T
        if eta >= 0:
            print("eta >= 0")
            return 0
        # 步骤4：更新alpha_j
        os.alphas[j] -= os.label_mat[j] * (Ei - Ej) / eta
        # 步骤5：修剪alpha_j
        os.alphas[j] = clip_alpha(os.alphas[j], H, L)
        # 更新Ej至误差缓存
        update_Ek(os, j)
        if(abs(os.alphas[j] - alpha_j_old) < 0.00001):
            print("alpha_j变化太小")
            return 0
        # 步骤6：更新alpha_i
        os.alphas[i] += os.label_mat[i] * os.label_mat[j] * (alpha_j_old - os.alphas[j])
        # 更新Ei至误差缓存
        update_Ek(os, i)
        # 步骤7：更新b_1和b_2:
        b1 = os.b - Ei - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[i, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[i, :].T
        b2 = os.b - Ej - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[j, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[j, :].T
        # 步骤8：根据b_1和b_2更新b
        if(0 < os.alphas[i] < os.C):
            os.b = b1
        elif(0 < os.alphas[j] < os.C):
            os.b = b2
        else:
            os.b = (b1 + b2) / 2.0
        return 1 #表示有更新
    else:
        return 0 #表示没有更新

"""
函数说明：完整的线性SMO算法
Parameters:
    dataMatIn - 数据矩阵
    classLabels - 数据标签
    C - 松弛变量
    toler - 容错率
    maxIter - 最大迭代次数
Returns:
    oS.b - SMO算法计算的b
    oS.alphas - SMO算法计算的alphas
"""
def smo_p(data_mat_in, class_labels, C, toler, max_iter):
    # 初始化数据结构
    os = opt_struct(np.mat(data_mat_in), np.mat(class_labels).transpose(), C, toler)
    # 初始化当前迭代次数
    iter = 0
    entrie_set = True #是否在全部数据集上迭代
    alpha_pairs_changed = 0
    # 遍历整个数据集alpha都没有更新或者超过最大迭代次数，则退出循环
    while(iter < max_iter) and ((alpha_pairs_changed > 0) or (entrie_set)):
        alpha_pairs_changed = 0
        if entrie_set: # 遍历整个数据集
            for i in range(os.m):
                # 使用优化的SMO算法
                alpha_pairs_changed += innerL(i, os) #innerL返回的值是0或者1
                print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
            iter += 1
        # 遍历非边界值
        else:
            # 遍历不在边界0和C的alpha
            non_bound_i_s = np.nonzero((os.alphas.A > 0) * (os.alphas.A < C))[0]
            for i in non_bound_i_s:
                alpha_pairs_changed += innerL(i, os)
                print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
            iter += 1
        # 遍历一次后改为非边界遍历
        if entrie_set:
            entrie_set = False #进行切换，遍历非边界数据集
        # 如果alpha没有更新，计算全样本遍历
        elif(alpha_pairs_changed == 0):
            entrie_set = True
        print("迭代次数:%d" % iter)
    # 返回SMO算法计算的b和alphas
    return os.b, os.alphas

"""
函数说明：分类结果可视化
Returns:
    dataMat - 数据矩阵
    classLabels - 数据标签
    w - 直线法向量
    b - 直线截距 
Returns:
    None
"""
def show_classifer(data_mat, class_labels, w, b):
    data_mat, label_mat = load_dataset('testSet.txt')
    fig = plt.figure()
    ax = fig.add_subplot(111)
    cm = mpl.colors.ListedColormap(['g', 'b'])
    ax.scatter(np.array(data_mat)[:, 0], np.array(data_mat)[:, 1], c=np.array(label_mat), cmap=cm, s=20)
    # 绘制直线
    x1 = max(data_mat)[0]
    x2 = min(data_mat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b - a1 * x1) / a2, (-b - a1 * x2) / a2
    plt.plot([x1, x2], [y1, y2])
    # 找出支持向量点
    # enumerate在字典上是枚举、列举的意思
    for i, alpha in enumerate(alphas):
        # 支持向量机的点
        if (abs(alpha) > 0):
            x, y = data_mat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolors='red')
    plt.show()
"""
函数说明：计算w
Returns:
    dataArr - 数据矩阵
    classLabels - 数据标签
    alphas - alphas值
Returns:
    w - 直线法向量
"""
def cal_w_s(alphas, data_array, class_labels):
    X = np.mat(data_array)
    label_mat = np.mat(class_labels).transpose()
    m, n = np.shape(X)
    w = np.zeros((n, 1))
    for i in range(m):
        w += np.multiply(alphas[i] * label_mat[i], X[i, :].T)
    return w


if __name__ == '__main__':
    data_array, class_labels = load_dataset('testSet.txt')
    b, alphas = smo_p(data_array, class_labels, 0.6, 0.001, 40)
    w = cal_w_s(alphas, data_array, class_labels)
    print('b=', b)
    print('w=', w)
    print('alphas=', alphas[alphas > 0])
    for i in range(100):
        if alphas[i] > 0.0:
            print(data_array[i], class_labels[i])
    show_classifer(data_array, class_labels, w, b)