Python-数据分析秘籍-四-

布客飞龙
73 阅读46分钟

Python 数据分析秘籍（四）

原文：Mastering Python Data Analysis

协议：CC BY-NC-SA 4.0

十、评估分类器、回归器和聚类

在本章中，我们将介绍以下食谱:

  • 用混淆矩阵直接分类
  • 计算精度、召回率和 F1 分数
  • 检查接收器工作特性和曲线下的区域
  • 想象契合度
  • 计算均方误差和中值绝对误差
  • 用平均轮廓系数评价聚类
  • 将结果与虚拟分类器进行比较
  • 确定 MAPE 和多相萃取
  • 与虚拟回归器进行比较
  • 计算平均绝对误差和残差平方和
  • 检查分类的 kappa
  • 看看马修斯相关系数

简介

评估分类器、回归器和聚类是一个涉及多个方面的关键多维问题。纯粹从工程的角度来看，我们担心速度、内存和正确性。在某些情况下，速度就是一切。如果记忆不足，我们当然要优先考虑。这个世界是一个充满选择的巨大迷宫，你有时被迫选择一个模型而不是其他模型，而不是在一个集合中使用多个模型。当然，我们应该用适当的评估标准来指导我们的理性决策。

有这么多的评估标准，你需要多本书来描述它们。显然，许多指标非常相似。其中一些被接受和流行，在这些度量中，一些在 scikit-learn 中实现。

我们将评估第 9 章、集成学习和降维中的分类器和回归器。我们将这些估计应用于天气预报的样本问题。这不一定是人类擅长的问题。实现人的表现是一些问题的目标，例如人脸识别、字符识别、垃圾邮件分类和情感分析。作为击败的底线，我们经常选择某种形式的随机猜测。

用混淆矩阵直接分类

准确性是一个衡量模型在给定环境中表现如何的指标。准确性是 scikit-learn 分类器的默认评估指标。不幸的是，准确性是一维的，当类不平衡时，它没有帮助。我们在第 9 章、集成学习和降维中检查的降雨数据相当均衡。雨天的数量几乎等于不下雨的天数。在电子邮件垃圾邮件分类的情况下，至少对我来说，平衡转向了垃圾邮件。

A 混淆矩阵是一个表格，通常用来总结分类的结果。表的两个维度是预测类和目标类。在二元分类的背景下，我们谈论正类和负类。将一个类命名为负数是任意的——这并不一定意味着它在某种程度上是坏的。我们可以将任何多类问题简化为一类，而不是问题的其余部分；所以，当我们评估二元分类时，我们可以将框架扩展到多类分类。一个类可以被正确预测，也可以不被正确预测；我们相应地用“真”和“假”来标记这些实例。

我们有四种真、假、阳性和阴性组合，如下表所述:

|   |

预测类别

| | --- | --- | | 目标类 | 真阳性:下雨了，我们预测的没错。 | 误报:我们错误的预测了会下雨。 | | 假阴性:确实下雨了，但是我们预测不会下雨。 | 真否定:没下雨，我们预测的没错。 |

怎么做...

  1. 进口情况如下:

    import numpy as np
from sklearn.metrics import confusion_matrix
import seaborn as sns
import dautil as dl
from IPython.display import HTML
import ch10util

  2. 定义以下函数绘制混淆矩阵:

    def plot_cm(preds, y_test, title, cax):
    cm = confusion_matrix(preds.T, y_test)
    normalized_cm = cm/cm.sum().astype(float)
    sns.heatmap(normalized_cm, annot=True, fmt='.2f', vmin=0, vmax=1,
                xticklabels=['Rain', 'No Rain'],
                yticklabels=['Rain', 'No Rain'], ax=cax)
    cax.set_xlabel('Predicted class')
    cax.set_ylabel('Expected class')
    cax.set_title('Confusion Matrix for Rain Forecast | ' + title)

  3. 加载目标值并绘制随机森林分类器、打包分类器、投票和堆叠分类器的混淆矩阵:

    y_test = np.load('rain_y_test.npy')
sp = dl.plotting.Subplotter(2, 2, context)

plot_cm(y_test, np.load('rfc.npy'), 'Random Forest', sp.ax)

plot_cm(y_test, np.load('bagging.npy'), 'Bagging', sp.next_ax())

plot_cm(y_test, np.load('votes.npy'),'Votes', sp.next_ax())

plot_cm(y_test, np.load('stacking.npy'), 'Stacking', sp.next_ax())
sp.fig.text(0, 1, ch10util.classifiers())
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

源代码在本书代码包的 conf_matrix.ipynb文件中。

它是如何工作的

我们为四个分类器显示了四个混淆矩阵，每个矩阵的四个数字似乎在重复。当然，数字并不完全相等；然而，你必须考虑到一些随机变化。

另见

计算精度、召回率和 F1 分数

在用混淆矩阵配方直接分类中，您了解到我们可以将分类的样本标记为真阳性、假阳性、真阴性和假阴性。通过这些类别的计数，我们可以计算的许多评估指标，我们将在本食谱中涵盖四个，如下式所示:

Computing precision, recall, and F1-score

这些指标的范围从零到一，零是最差的理论分数，一是最好的。实际上，最差的分数是我们通过随机猜测得到的。实践中的最佳分数可能低于 1，因为在某些情况下，我们只能希望模仿人类的表现，对于正确的分类应该是什么可能存在歧义，例如，在情感分析的情况下(涵盖在 Python 数据分析一书中)。

  • 精确度(10.1)是正确预测的比率。
  • 精度 (10.2)将相关性度量为将阴性类别样本分类为阳性的可能性。选择哪个班是积极的多少有些武断，但我们可以说，未雨绸缪是积极的。高精度意味着我们将相对较少的非雨天(负)标记为雨天。对于搜索(网络、数据库或其他)，这意味着相关结果的数量相对较高。
  • 回忆 (10.3)是找到所有阳性样本的可能性。如果再来一次，雨天是我们的正课，雨天分类越正确，召回率越高。对于搜索，我们可以通过返回所有文档来获得完美的召回，因为这将自动返回所有相关文档。人脑有点像数据库，在这种情况下，回忆意味着记住某个 Python 函数如何工作的可能性。
  • F1 分数 (10.4)是精度和召回率的调和平均值(实际上，F1 分数有多个变化)。G 分数使用几何平均值；但是，据我所知，它不太受欢迎。F1 成绩背后的理念，相关的 F 成绩和 G 成绩，就是把精准性和召回率结合起来。这并不一定是最好的衡量标准。还有其他指标你可能更喜欢，比如马修斯相关系数(参考看一下马修斯相关系数配方)和科恩卡帕(参考检验卡帕的分类配方)。当我们开始选择这么多分类指标时，我们显然想要最好的指标。然而，你必须根据你的情况做出选择，因为没有一个标准适合所有人。

怎么做...

  1. 进口情况如下:

    import numpy as np
from sklearn import metrics
import ch10util
import dautil as dl
from IPython.display import HTML

  2. 加载目标值并计算指标:

    y_test = np.load('rain_y_test.npy')
accuracies = [metrics.accuracy_score(y_test, preds)
              for preds in ch10util.rain_preds()]
precisions = [metrics.precision_score(y_test, preds)
              for preds in ch10util.rain_preds()]
recalls = [metrics.recall_score(y_test, preds)
           for preds in ch10util.rain_preds()]
f1s = [metrics.f1_score(y_test, preds)
       for preds in ch10util.rain_preds()]

  3. 绘制降雨预报的指标:

    sp = dl.plotting.Subplotter(2, 2, context)
ch10util.plot_bars(sp.ax, accuracies)
sp.label()

ch10util.plot_bars(sp.next_ax(), precisions)
sp.label()

ch10util.plot_bars(sp.next_ax(), recalls)
sp.label()

ch10util.plot_bars(sp.next_ax(), f1s)
sp.label()
sp.fig.text(0, 1, ch10util.classifiers())
HTML(sp.exit())

最终结果见截图后的:

How to do it...

代码在本书代码包的precision_recall.ipynb文件中。

另见

检查接收器操作特性和曲线下的区域

接收器 操作特性 ( ROC )是二进制分类器的召回率(10.3)和假阳性率 ( FPR )的图。FPR 由下式给出:

Examining a receiver operating characteristic and the area under a curve

在本食谱中，我们将绘制我们在第 9 章、集成学习和降维中使用的各种分类器的 ROC。此外，我们将绘制与随机猜测相关的曲线和理想曲线。显然，我们想要超越基线，尽可能地接近理想曲线。

曲线下的 区域 ( AUC 、 ROC AUC 或 AUROC )是总结 ROC 的另一个评价指标。AUC 也可以用来比较模型，但它提供的信息比 ROC 少。

怎么做...

  1. 进口情况如下:

    from sklearn import metrics
import numpy as np
import ch10util
import dautil as dl
from IPython.display import HTML

  2. 加载数据并计算指标:

    y_test = np.load('rain_y_test.npy')
roc_aucs = [metrics.roc_auc_score(y_test, preds)
            for preds in ch10util.rain_preds()]

  3. 绘制雨预报器的 AUROC:

    sp = dl.plotting.Subplotter(2, 1, context)
ch10util.plot_bars(sp.ax, roc_aucs)
sp.label()

  4. 绘制降雨预报器的 ROC 曲线:

    cp = dl.plotting.CyclePlotter(sp.next_ax())

for preds, label in zip(ch10util.rain_preds(),
                        ch10util.rain_labels()):
    fpr, tpr, _ = metrics.roc_curve(y_test, preds,
                                    pos_label=True)
    cp.plot(fpr, tpr, label=label)

fpr, tpr, _ = metrics.roc_curve(y_test, y_test)
sp.ax.plot(fpr, tpr, 'k', lw=4, label='Ideal')
sp.ax.plot(np.linspace(0, 1), np.linspace(0, 1), 
           '--', label='Baseline')
sp.label()
sp.fig.text(0, 1, ch10util.classifiers())
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的roc_auc.ipynb文件中。

另见

可视化拟合优度

我们期望，或者至少希望，回归的残差只是随机噪声。如果不是这样，那么我们的回归者可能忽略了信息。我们期望残差是独立的并且是正态分布的。用直方图或者 QQ 图比较容易检查。一般来说，我们希望残差的均值尽可能接近零，并且希望残差的方差尽可能小。理想拟合会有零值残差。

怎么做...

  1. 进口情况如下:

    import numpy as np
import matplotlib.pyplot as plt
import dautil as dl
import seaborn as sns
from scipy.stats import probplot
from IPython.display import HTML

  2. 加载助推回归器的目标和预测:

    y_test = np.load('temp_y_test.npy')
preds = np.load('boosting.npy')

  3. 将实际值和预测值绘制如下:

    sp = dl.plotting.Subplotter(2, 2, context)
cp = dl.plotting.CyclePlotter(sp.ax)
cp.plot(y_test)
cp.plot(preds)
sp.ax.set_ylabel(dl.data.Weather.get_header('TEMP'))
sp.label()

  4. 将残差绘制如下:

    residuals = preds - y_test
sp.next_ax().plot(residuals)
sp.label()

  5. 绘制残差分布图:

    sns.distplot(residuals, ax=sp.next_ax())
sp.label()

  6. 绘制残差的 QQ 图:

    probplot(residuals, plot=sp.next_ax())
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的visualizing_goodness.ipynb文件中。

另见

计算均方误差和中值绝对误差

均方误差 ( 均方误差)和中值 绝对误差 ( 梅达)是流行的回归指标。它们由以下等式给出:

Computing MSE and median absolute error

均方误差(10.6)类似于总体方差。因此，均方误差的平方根( RMSE )类似于标准偏差。均方误差的单位与被分析的变量相同——在我们的例子中，是温度。理想拟合具有零值残差，因此其均方误差等于零。由于我们处理的是平方误差，均方误差的值大于或理想地等于零。

MedAE 类似于 MSE，但是我们从残差的绝对值开始，我们使用中值而不是平均值作为中心性的度量。MedAE 也类似于方差，理想情况下为零或非常小。取绝对值而不是平方有可能避免数值不稳定和速度问题，并且中值对于异常值比平均值更稳健。此外，取平方倾向于强调较大的误差。

在本食谱中，我们将为来自第 9 章、集成学习和降维的回归器绘制 MSE 和 MedAE 的自举种群。

怎么做...

  1. 进口情况如下:

    from sklearn import metrics
import ch10util
from IPython.display import HTML
import dautil as dl
from IPython.display import HTML

  2. 绘制温度预测指标的分布:

    sp = dl.plotting.Subplotter(3, 2, context)
ch10util.plot_bootstrap('boosting',
                        metrics.mean_squared_error, sp.ax)
sp.label()

ch10util.plot_bootstrap('boosting',
                        metrics.median_absolute_error, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('etr',
                        metrics.mean_squared_error, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('etr',
                        metrics.median_absolute_error, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('ransac',
                        metrics.mean_squared_error, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('ransac',
                        metrics.median_absolute_error, sp.next_ax())
sp.label()
sp.fig.text(0, 1, ch10util.regressors())
HTML(sp.exit())

最终结果参见以下截图:

How to do it...

代码在本书代码包的mse.ipynb文件中。

另见

用平均轮廓系数评价聚类

聚类是一种无监督的机器学习类型的分析。虽然我们一般不知道什么是最好的聚类，但我们仍然可以得到一个聚类结果有多好的概念。一种方法是计算轮廓系数，如下式所示:

Evaluating clusters with the mean silhouette coefficient

在上式中， a(i) 是样本 i 相对于同一聚类中其他样本的平均相异度。一个小的 a(i) 表示样品属于它的簇。 b(i) 是 i 与其他聚类的最低平均相异度。它为 i 指示下一个最佳集群。如果样本的轮廓系数 s(i) 接近 1，则表示样本分配正确。 s(i) 的值在-1 到 1 之间变化。所有样本的轮廓系数的平均值衡量聚类的质量。

我们可以使用平均轮廓系数来通知我们对 K-均值聚类算法的聚类数量的决定。在第 5 章、网络挖掘、数据库和大数据中的使用 Spark 配方对流式数据进行聚类中更详细地介绍了 K-means 聚类算法。

怎么做...

  1. 进口情况如下:

    import dautil as dl
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.metrics import silhouette_samples
from IPython.display import HTML

  2. 定义以下函数绘制轮廓样本:

    def plot_samples(ax, years, labels, i, avg):
    silhouette_values = silhouette_samples(X, labels)
    dl.plotting.plot_text(ax, years, silhouette_values,
                          labels, add_scatter=True)
    ax.set_title('KMeans k={0} Silhouette avg={1:.2f}'.format(i, avg))
    ax.set_xlabel('Year')
    ax.set_ylabel('Silhouette score')

  3. 加载数据并重新采样，如下所示:

    df = dl.data.Weather.load().resample('A').dropna()
years = [d.year for d in df.index]
X = df.values

  4. 绘制不同簇数的簇:

    sp = dl.plotting.Subplotter(2, 2, context)
avgs = []
rng = range(2, 9)

for i in rng:
    kmeans = KMeans(n_clusters=i, random_state=37)
    labels = kmeans.fit_predict(X)
    avg = silhouette_score(X, labels)
    avgs.append(avg)

    if i < 5:
        if i > 2:
            sp.next_ax()

        plot_samples(sp.ax, years, labels, i, avg)

sp.next_ax().plot(rng, avgs)
sp.label()
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的 evaluating_clusters.ipynb文件中。

另见

用虚拟分类器比较结果

scikit-learn DummyClassifier类实现了几种随机猜测的策略，可以作为分类器的基线。策略如下:

  • stratified:这使用了训练集类分布
  • most_frequent:这预示着最频繁的上课
  • prior:这在 scikit-learn 0.17 中提供，通过最大化课程优先级进行预测
  • uniform:这使用均匀分布对类进行随机采样
  • constant:这预测了用户指定的类

可以看到DummyClassifier类的一些策略总是预测同一个类。这可能会导致一些 scikit-learn 度量函数发出警告。我们将执行与在计算精度、召回率和 F1 评分配方中相同的分析，但添加了虚拟分类器。

怎么做...

  1. 进口情况如下:

    import numpy as np
from sklearn import metrics
import ch10util
from sklearn.dummy import DummyClassifier
from IPython.display import HTML
import dautil as dl

  2. 加载数据如下:

    y_test = np.load('rain_y_test.npy')
X_train = np.load('rain_X_train.npy')
X_test = np.load('rain_X_test.npy')
y_train = np.load('rain_y_train.npy')

  3. 创建虚拟分类器并使用它们进行预测:

    stratified = DummyClassifier(random_state=28)
frequent = DummyClassifier(strategy='most_frequent',
                           random_state=28)
prior = DummyClassifier(strategy='prior', random_state=29)
uniform = DummyClassifier(strategy='uniform',
                          random_state=29)
preds = ch10util.rain_preds()

for clf in [stratified, frequent, prior, uniform]:
    clf.fit(X_train, y_train)
    preds.append(clf.predict(X_test))

  4. 使用预测计算指标，如下所示:

    accuracies = [metrics.accuracy_score(y_test, p)
              for p in preds]
precisions = [metrics.precision_score(y_test, p)
              for p in preds]
recalls = [metrics.recall_score(y_test, p)
           for p in preds]
f1s = [metrics.f1_score(y_test, p)
       for p in preds]

  5. 绘制虚拟和常规分类器的度量:

    labels = ch10util.rain_labels()
labels.extend(['stratified', 'frequent',
               'prior', 'uniform'])

sp = dl.plotting.Subplotter(2, 2, context)
ch10util.plot_bars(sp.ax, accuracies, labels, rotate=True)
sp.label()

ch10util.plot_bars(sp.next_ax(), precisions, labels, rotate=True)
sp.label()

ch10util.plot_bars(sp.next_ax(), recalls, labels, rotate=True)
sp.label()

ch10util.plot_bars(sp.next_ax(), f1s, labels, rotate=True)
sp.label()
sp.fig.text(0, 1, ch10util.classifiers(), fontsize=10)
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的 dummy_clf.ipynb文件中。

另见

测定 MAPE 和多相萃取

平均百分比误差 (MPE )和平均绝对百分比误差 ( MAPE )将预测误差表示为比率，因此它们是无量纲的，易于解释。从下面的等式中可以看出，MPE 和 MAPE 的缺点是我们有被零除的风险:

Determining MAPE and MPE

目标变量等于零是完全有效的。对于温度，这恰好是冰点。冬天经常会结冰，所以我们要么忽略那些观测值，要么加入一个足够大的常数，避免被零值除。在下一节中，很明显，简单地忽略观察会导致奇怪的引导分布。

怎么做...

  1. 进口情况如下:

    import ch10util
import dautil as dl
from IPython.display import HTML

  2. 将自举指标绘制如下:

    sp = dl.plotting.Subplotter(3, 2, context)
ch10util.plot_bootstrap('boosting',
                        dl.stats.mape, sp.ax)
sp.label()

ch10util.plot_bootstrap('boosting',
                        dl.stats.mpe, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('etr',
                        dl.stats.mape, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('etr',
                        dl.stats.mpe, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('ransac',
                        dl.stats.mape, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('ransac',
                        dl.stats.mpe, sp.next_ax())
sp.label()
sp.fig.text(0, 1, ch10util.regressors())
HTML(sp.exit())

最终结果参见下面的截图:

How to do it...

代码在本书代码包的mape_mpe.ipynb文件中。

另见

与虚拟回归器进行比较

scikit-learn DummyRegressor类实现了随机猜测的几种策略，可以作为回归器的基线。这些战略如下:

  • mean:这预测了训练集的均值。
  • median:这预测了训练集的中位数。
  • quantile:当提供quantile参数时，这预测训练集的指定分位数。我们将通过指定第一个和第三个四分位数来应用这个策略。
  • constant:这预测了由用户提供的恒定值。

我们将使用 R 平方、均方误差、最大似然误差和最大似然误差将虚拟回归器与第 9 章、集成学习和降维中的回归器进行比较。

怎么做...

  1. 进口情况如下:

    import numpy as np
from sklearn.dummy import DummyRegressor
import ch10util
from sklearn import metrics
import dautil as dl
from IPython.display import HTML

  2. 加载温度数据如下:

    y_test = np.load('temp_y_test.npy')
X_train = np.load('temp_X_train.npy')
X_test = np.load('temp_X_test.npy')
y_train = np.load('temp_y_train.npy')

  3. 使用可用的策略创建虚拟回归，并用它们预测温度:

    mean = DummyRegressor()
median = DummyRegressor(strategy='median')
q1 = DummyRegressor(strategy='quantile', quantile=0.25)
q3 = DummyRegressor(strategy='quantile', quantile=0.75)

preds = ch10util.temp_preds()

for reg in [mean, median, q1, q3]:
    reg.fit(X_train, y_train)
    preds.append(reg.predict(X_test))

  4. 计算常规回归和虚拟回归的 R 平方、均方误差、中值绝对误差和平均百分比误差:

    r2s = [metrics.r2_score(p, y_test) for p in preds]
mses = [metrics.mean_squared_error(p, y_test)
        for p in preds]
maes = [metrics.median_absolute_error(p, y_test)
        for p in preds]
mpes = [dl.stats.mpe(y_test, p) for p in preds]

labels = ch10util.temp_labels()
labels.extend(['mean', 'median', 'q1', 'q3'])

  5. 将指标绘制如下:

    sp = dl.plotting.Subplotter(2, 2, context)
ch10util.plot_bars(sp.ax, r2s, labels)
sp.label()

ch10util.plot_bars(sp.next_ax(), mses, labels)
sp.label()

ch10util.plot_bars(sp.next_ax(), maes, labels)
sp.label()

ch10util.plot_bars(sp.next_ax(), mpes, labels)
sp.label()
sp.fig.text(0, 1, ch10util.regressors())
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的dummy_reg.ipynb文件中。

另见

计算平均绝对误差和残差平方和

平均绝对误差 ( 平均)和 残差平方和 ( RSS )是由以下公式给出的回归度量:

Calculating the mean absolute error and the residual sum of squares

平均绝对误差(10.11)类似于均方误差和最大误差，但在计算的一个步骤中有所不同。这些度量的共同特征是它们忽略了误差的符号，类似于方差。平均值大于或理想地等于零。

RSS (10.12)类似于 MSE，只是我们没有除以残差数。因此，使用 RSS 可以获得更大的值。然而，一个理想的适合给你一个零 RSS。

怎么做...

  1. 进口情况如下:

    import ch10util
import dautil as dl
from sklearn import metrics
from IPython.display import HTML

  2. 将自举指标绘制如下:

    sp = dl.plotting.Subplotter(3, 2, context)
ch10util.plot_bootstrap('boosting',
                        metrics.mean_absolute_error, sp.ax)
sp.label()

ch10util.plot_bootstrap('boosting',
                        dl.stats.rss, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('etr',
                        metrics.mean_absolute_error, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('etr',
                        dl.stats.rss, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('ransac',
                        metrics.mean_absolute_error, sp.next_ax())
sp.label()

ch10util.plot_bootstrap('ransac',
                        dl.stats.rss, sp.next_ax())
sp.label()
sp.fig.text(0, 1, ch10util.regressors())
HTML(sp.exit())

最终结果参见以下截图:

How to do it...

代码在本书代码包的mae_rss.ipynb文件中。

另见

检验分类的 kappa

Cohen 的 kappa 测量目标和预测类别之间的一致性，类似于精确度，但是它也考虑了获得预测的随机机会。科恩的 kappa 由下式给出:

Examining the kappa of classification

在这个方程中， p 0 是相对观测一致性， p e 是由数据推导出的随机一致性概率。Kappa 在负值和以下兰迪斯和科赫的粗略分类之间变化:

  • 不一致:kappa < 0
  • 略微一致:kappa = 0 至 0.2
  • 公平协议:kappa = 0.21 至 0.4
  • 中度一致:kappa = 0.41 至 0.6
  • 良好一致性:kappa = 0.61 至 0.8
  • 非常好的一致性:kappa = 0.81 到 1.0

我知道另外两个评定卡帕等级的方案，所以这些数字不是一成不变的。我想我们可以同意不接受低于 0.2 的 kappa。最合适的用例当然是对模型进行排名。科恩的 kappa 还有其他变体，但截至 2015 年 11 月，它们还没有在 scikit-learn 中实现。scikit-learn 0.17 通过cohen_kappa_score()功能增加了对 Cohen 的 kappa 的支持。

怎么做...

  1. 进口情况如下:

    import dautil as dl
from sklearn import metrics
import numpy as np
import ch10util
from IPython.display import HTML

  2. 计算雨预报器的准确度、精确度、召回率、F1 评分和 kappa:

    y_test = np.load('rain_y_test.npy')
accuracies = [metrics.accuracy_score(y_test, preds)
              for preds in ch10util.rain_preds()]
precisions = [metrics.precision_score(y_test, preds)
              for preds in ch10util.rain_preds()]
recalls = [metrics.recall_score(y_test, preds)
           for preds in ch10util.rain_preds()]
f1s = [metrics.f1_score(y_test, preds)
       for preds in ch10util.rain_preds()]
kappas = [metrics.cohen_kappa_score(y_test, preds)
          for preds in ch10util.rain_preds()]

  3. 散点图针对 kappa 的指标如下:

    sp = dl.plotting.Subplotter(2, 2, context)
dl.plotting.plot_text(sp.ax, accuracies, kappas,
                      ch10util.rain_labels(), add_scatter=True)
sp.label()

dl.plotting.plot_text(sp.next_ax(), precisions, kappas,
                      ch10util.rain_labels(), add_scatter=True)
sp.label()

dl.plotting.plot_text(sp.next_ax(), recalls, kappas,
                      ch10util.rain_labels(), add_scatter=True)
sp.label()

dl.plotting.plot_text(sp.next_ax(), f1s, kappas,
                      ch10util.rain_labels(), add_scatter=True)
sp.label()
sp.fig.text(0, 1, ch10util.classifiers())                     
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的kappa.ipynb文件中。

它是如何工作的

从前两个图中，我们可以得出结论，bagging 分类器具有最高的准确度、精确度和 kappa。所有分类器的 kappa 都在 0.2 以上，所以它们至少在一定程度上是可以接受的。

另见

看一下马修斯相关系数

马修斯相关系数 ( MCC )或φ系数是 Brian Matthews 在 1975 年发明的二进制分类的评价指标。监控中心是目标和预测的相关系数，在-1 和 1 之间变化(最佳一致性)。MCC 是总结混淆矩阵的一个非常好的方法(参考用混淆矩阵配方直接得到分类)，因为它使用了其中的所有四个数字。监控中心由以下等式给出:

Taking a look at the Matthews correlation coefficient

怎么做...

  1. 进口情况如下:

    import dautil as dl
from sklearn import metrics
import numpy as np
import ch10util
from IPython.display import HTML

  2. 计算降雨预报器的准确度、精确度、召回率、F1 分数和马修斯相关系数:

    y_test = np.load('rain_y_test.npy')
accuracies = [metrics.accuracy_score(y_test, preds)
              for preds in ch10util.rain_preds()]
precisions = [metrics.precision_score(y_test, preds)
              for preds in ch10util.rain_preds()]
recalls = [metrics.recall_score(y_test, preds)
           for preds in ch10util.rain_preds()]
f1s = [metrics.f1_score(y_test, preds)
       for preds in ch10util.rain_preds()]
mc = [metrics.matthews_corrcoef(y_test, preds)
      for preds in ch10util.rain_preds()]

  3. 将指标绘制如下:

    sp = dl.plotting.Subplotter(2, 2, context)
dl.plotting.plot_text(sp.ax, accuracies, mc,
                      ch10util.rain_labels(), add_scatter=True)
sp.label()

dl.plotting.plot_text(sp.next_ax(), precisions, mc,
                      ch10util.rain_labels(), add_scatter=True)
sp.label()

dl.plotting.plot_text(sp.next_ax(), recalls, mc,
                      ch10util.rain_labels(), add_scatter=True)
sp.label()

dl.plotting.plot_text(sp.next_ax(), f1s, mc,
                      ch10util.rain_labels(), add_scatter=True)
sp.label()
sp.fig.text(0, 1, ch10util.classifiers())
HTML(sp.exit())

最终结果参见下面的截图:

How to do it...

代码在本书的代码包matthews_correlation.ipynb文件中。

另见

十一、分析图像

在本章中，我们将介绍以下食谱:

  • 设置 OpenCV
  • 应用尺度不变特征变换
  • 使用 SURF 检测要素
  • 量化颜色
  • 图像去噪
  • 从图像中提取面片
  • 用哈尔级联检测人脸
  • 寻找明亮的星星
  • 从图像中提取元数据
  • 从图像中提取纹理特征
  • 在图像上应用分层聚类
  • 用谱聚类分割图像

简介

图像处理是一个非常大的研究领域。用于图像处理的技术通常也可以应用于视频分析。我们可以将图像处理视为一种特殊类型的信号处理。信号处理包含在第 6 章、信号处理和时间序列中。然而，图像带来了特殊的挑战，例如高维数(我们可以将每个图像像素定义为一个特征)和空间依赖性(像素位置很重要)。

与计算机相比，人类的视觉系统非常先进。我们能够识别物体、面部表情和物体运动。显然，这与食肉动物和它们食用人肉的倾向有关。在本章中，我们将集中精力寻找图像中的特征并对图像像素进行聚类(分割)，而不是试图理解人类视觉是如何工作的。

在这一章中，我们经常使用 OpenCV 库，由于它是一个相当大的库，我决定只为这一章创建一个特殊的 Docker 容器。你可能已经知道了，我制作了一个名为pydacbk的 Docker 图像。嗯，本章的 Docker 容器名为pydacbk11。

设置开放式课程

OpenCV ( 开源计算机视觉)是 2000 年创建的计算机视觉库，目前由 Itseez 维护。OpenCV 是用 C++编写的，但它也有针对 Python 和其他编程语言的绑定。支持多种操作系统和图形处理器。本章没有足够的篇幅来涵盖 OpenCV 的所有功能。即使是一本书可能也是不够的——对于皮德尼塔斯，我推荐约瑟夫·豪斯的《用 Python 打开计算机视觉》。

OpenCV 2.x.x 包中的一些第三方专利算法，如 SIFT 和 SURF(参考本章中的相关食谱)，已经被转移到一个特殊的 GitHub 存储库中。您仍然可以使用它们，但是您需要在安装过程中明确包含它们。

OpenCV 构建过程有许多选项。如果您不确定哪些选项最适合您，请阅读 OpenCV 文档或使用适合您的操作系统的软件包管理器。一般来说，你不应该使用太多的选项。尽管您可以灵活地关闭某些模块，但其他模块可能会依赖于它们，这可能会导致一连串的错误。

做好准备

如果你在 Windows 或者 Fedora 上，可以在http://docs . opencv . org/3 . 0 . 0/da/df6/tutorial _ py _ table _ contents _ setup . html阅读相应教程(2015 年 12 月检索)。对于基于 Ubuntu 的 Docker 容器，我需要用以下命令安装一些先决条件:

$ apt-get update
$ apt-get install -y cmake make git g++

怎么做...

以下说明作为示例，并对您的设置做了一些假设。例如，它假设您正在 Python 3 中使用 Anaconda。为了方便起见，我将基于 Ubuntu 的 Docker 容器的所有指令组织在一个 shell 脚本中；但是，如果您愿意，也可以在终端中分别键入每一行。

  1. 下载核心 OpenCV 项目的代码(如果没有 Git，也可以从 GitHub 网站下载代码):

    $ cd /opt
$ git clone https://github.com/Itseez/opencv.git
$ cd opencv
$ git checkout tags/3.0.0

  2. 下载 OpenCV(第三方)投稿的代码(如果没有 Git，也可以从 GitHub 网站下载代码):

    $ cd /opt
$ git clone https://github.com/Itseez/opencv_contrib
$ cd opencv_contrib
$ git checkout tags/3.0.0
$ cd /opt/opencv

  3. 创建一个构建目录并导航到它:

    $ mkdir build
$ cd build

  4. 此步骤显示了您可以使用的一些构建选项(您不必使用所有这些选项):

    $ ANACONDA=~/anaconda
$ cmake -D CMAKE_BUILD_TYPE=RELEASE \
 -D BUILD_PERF_TESTS=OFF \
 -D BUILD_opencv_core=ON \
 -D BUILD_opencv_python2=OFF \
 -D BUILD_opencv_python3=ON \
 -D BUILD_opencv_cuda=OFF \
 -D BUILD_opencv_java=OFF \
 -D BUILD_opencv_video=ON \
 -D BUILD_opencv_videoio=ON \
 -D BUILD_opencv_world=OFF \
 -D BUILD_opencv_viz=ON \
 -D WITH_CUBLAS=OFF \
 -D WITH_CUDA=OFF \
 -D WITH_CUFFT=OFF \
 -D WITH_FFMPEG=OFF \
 -D PYTHON3_EXECUTABLE=${ANACONDA}/bin/python3 \
 -D PYTHON3_LIBRARY=${ANACONDA}/lib/libpython3.4m.so \
 -D PYTHON3_INCLUDE_DIR=${ANACONDA}/include/python3.4m \
 -D PYTHON3_NUMPY_INCLUDE_DIRS=${ANACONDA}/lib/python3.4/site-packages/numpy/core/include \
 -D PYTHON3_PACKAGES_PATH=${ANACONDA}/lib/python3.4/site-packages \
 -D BUILD_opencv_latentsvm=OFF \
 -D BUILD_opencv_xphoto=OFF \
 -D BUILD_opencv_xfeatures2d=ON \
 -D OPENCV_EXTRA_MODULES_PATH=/opt/opencv_contrib/modules \
 /opt/opencv

  5. 运行make命令(8 芯)安装如下:

    $ make -j8
$ sudo make install

它是如何工作的

前面的说明假设您是第一次安装 OpenCV。如果您是从 OpenCV 2.x.x 升级，您将不得不采取额外的预防措施。另外，我假设您不想要某些选项，并且正在 Python 3 中使用 Anaconda。下表解释了我们使用的一些构建选项:

|

[计]选项

|

描述

| | --- | --- | | BUILD_opencv_python2 | 支持 Python 2 | | BUILD_opencv_python3 | 支持 Python 3 | | BUILD_opencv_java | 支持 OpenCV Java 绑定 | | PYTHON3_EXECUTABLE | Python 3 可执行文件的位置 | | PYTHON3_LIBRARY | Python 3 库的位置 | | PYTHON3_INCLUDE_DIR | Python 3 include目录的位置 | | PYTHON3_NUMPY_INCLUDE_DIRS | NumPy include目录的位置 | | PYTHON3_PACKAGES_PATH | Python 3 包的位置 | | BUILD_opencv_xfeatures2d | 支持某些第三方算法，如 SIFT 和 SURF | | OPENCV_EXTRA_MODULES_PATH | OpenCV 第三方贡献代码的位置 |

还有更多

如果您没有足够的空间，例如在 Docker 容器中，那么您可以使用以下命令进行清理:

$ rm -rf /opt/opencv
$ rm -rf /opt/opencv_contrib

应用尺度不变特征变换

SIFT 算法(1999)在图像或视频中发现特征，并获得了不列颠哥伦比亚大学的专利。通常，我们可以使用特征进行分类或聚类。SIFT 在平移、缩放和旋转方面是不变的。

算法步骤如下:

  1. 使用高斯模糊滤镜以不同的比例模糊图像。
  2. 一个八度 对应的是两倍的滤镜标准差。将模糊的图像按八度音分组，并加以区别。
  3. 为差分图像找到跨尺度的局部极值。
  4. 将与局部极值相关的每个像素与同一尺度和相邻尺度中的相邻像素进行比较。
  5. 从比较中选择最大值或最小值。
  6. 拒绝对比度低的点。
  7. 插值候选关键点(图像特征)以获得原始图像上的位置。

做好准备

按照中的说明设置 OpenCV 配方。

怎么做...

  1. 进口情况如下:

    import cv2
import matplotlib.pyplot as plt
import dautil as dl
from scipy.misc import face

  2. 将原始图像绘制如下:

    img = face()
plt.title('Original')
dl.plotting.img_show(plt.gca(), img)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

  3. 将灰度图像绘制如下:

    plt.figure()
plt.title('Gray')
dl.plotting.img_show(plt.gca(), gray, cmap=plt.cm.gray)

  4. 用关键点(蓝色)绘制图像，如下所示:

    sift = cv2.xfeatures2d.SIFT_create()
(kps, descs) = sift.detectAndCompute(gray, None)
img2 = cv2.drawKeypoints(gray, kps, None, (0, 0, 255))

plt.figure()
plt.title('With Keypoints')
dl.plotting.img_show(plt.gca(), img2)

有关最终结果，请参考以下屏幕截图:

How to do it...

程序在本书代码包的applying_sift.ipynb文件中。

另见

利用 SURF 检测特征

加速鲁棒特征 ( SURF )是一种类似于 SIFT 并受其启发的专利算法(参考应用比例不变特征变换配方)。SURF 于 2006 年推出，使用哈尔小波(参考应用离散小波变换配方)。SURF 最大的优点就是比 SIFT 快。

看看下面的等式:

Detecting features with SURF

算法步骤如下:

  1. 如有必要，变换图像以获得等效灰度。
  2. 计算不同尺度下的 积分图像，它是一个像素上面和左边的像素之和，如公式(11.1)所示。积分图像代替了 SIFT 中的高斯滤波器。
  3. 将包含灰度图像二阶导数的 黑森矩阵 (11.2)定义为像素位置 p 和比例 σ (11.3)的函数。
  4. 行列式是与方阵相关的值。黑森矩阵的行列式对应于一个点的局部变化。选择行列式最大的点。
  5. 音阶σ由 11.3 定义，就像 SIFT 一样，我们可以定义音阶八度。SURF 通过改变过滤器内核的大小来工作，而 SIFT 改变图像大小。在比例和图像空间中插入上一步的最大值。
  6. 将哈尔小波变换应用于关键点周围的圆。
  7. 使用滑动窗口对响应求和。
  8. 根据响应总和确定方向。

做好准备

按照中的说明设置 OpenCV 配方。

怎么做...

  1. 进口如下:

    import cv2
import matplotlib.pyplot as plt
import dautil as dl

  2. 将原始图像绘制如下:

    img = cv2.imread('covers.jpg')
plt.title('Original')
dl.plotting.img_show(plt.gca(), img)

  3. 将灰度图像绘制如下:

    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

plt.figure()
plt.title('Gray')
dl.plotting.img_show(plt.gca(), gray, cmap=plt.cm.gray)
surf = cv2.xfeatures2d.SURF_create()
(kps, descs) = surf.detectAndCompute(gray, None)
img2 = cv2.drawKeypoints(gray, kps, None, (0, 0, 255))

  4. 用关键点(蓝色)绘制图像，如下所示:

    plt.figure()
plt.title('With Keypoints')
dl.plotting.img_show(plt.gca(), img2)

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的 applying_surf.ipynb文件中。

另见

量化颜色

在古代，电脑游戏几乎是单色的。许多年后，互联网允许我们下载图像，但网络速度很慢，所以首选颜色很少的紧凑图像。我们可以得出结论，限制颜色的数量是传统的。颜色是图像的一个维度，所以如果我们从图像中去除颜色，我们就可以说是降维。实际过程叫做 颜色量化。

通常，我们在三维空间中为每个像素表示 RGB ( 红色、绿色、蓝色)值，然后对这些点进行聚类。对于每个集群，我们都有相应的平均颜色。在这个配方中，我们将使用 k-means 聚类(参考用 Spark 配方对流式数据进行聚类)，尽管这不一定是最好的算法。

做好准备

按照中的说明设置 OpenCV 配方。

怎么做...

代码在本书代码包的quantizing_colors.ipynb文件中:

  1. 进口情况如下:

    import numpy as np
import cv2
import matplotlib.pyplot as plt
import dautil as dl
from scipy.misc import face

  2. 将原始图像绘制如下:

    sp = dl.plotting.Subplotter(2, 2, context)
img = face()
dl.plotting.img_show(sp.ax, img)
sp.label()
Z = img.reshape((-1, 3))

Z = np.float32(Z)

  3. 应用 k-均值聚类并绘制结果:

    criteria = (cv2.TERM_CRITERIA_MAX_ITER, 7, 1.0)

for k in [2, 4, 8]:
    _, label, center = cv2.kmeans(Z, k, None, criteria, 7,
                                  cv2.KMEANS_RANDOM_CENTERS)

    center = np.uint8(center)
    res = center[label.flatten()]
    res2 = res.reshape((img.shape))

    dl.plotting.img_show(sp.next_ax(), res2)
    sp.label()

最终结果参见以下截图:

How to do it...

另见

图像去噪

噪声是数据和图像中的常见现象。当然，噪音是不可取的，因为它不会给我们的分析增加任何价值。我们通常假设噪声通常分布在零附近。我们认为像素值是真实值和噪声(如果有的话)的总和。我们还假设噪声值是独立的，即一个像素的噪声值独立于另一个像素。

一个简单的想法是在一个小窗口中平均像素，因为我们假设噪声的期望值为零。这是模糊背后的一般想法。我们可以将这个想法更进一步，在一个像素周围定义多个窗口，然后我们可以对相似的面片进行平均。

OpenCV 有几个去噪功能，通常我们需要指定过滤器的强度、搜索窗口的大小以及模板窗口的大小来进行相似性检查。您应该注意不要将滤镜强度设置得太高，因为这可能会使图像不仅更干净，而且有点模糊。

做好准备

按照中的说明设置 OpenCV 配方。

怎么做...

  1. 进口情况如下:

    import cv2
import matplotlib.pyplot as plt
from sklearn.datasets import load_sample_image
import numpy as np
import dautil as dl

  2. 将原始图像绘制如下:

    img = load_sample_image('china.jpg')
dl.plotting.img_show(plt.gca(), img) 
plt.title('Original')
Z = img.reshape((-1, 3))

  3. 向图像添加噪点并绘制噪点图像:

    np.random.seed(59)
noise = np.random.random(Z.shape) < 0.99

noisy = (Z * noise).reshape((img.shape))

plt.figure()
plt.title('Noisy')
dl.plotting.img_show(plt.gca(), noisy)

  4. 清洁图像并显示:

    cleaned = cv2.fastNlMeansDenoisingColored(noisy, None, 10, 10, 7, 21)
plt.figure()
plt.title('Cleaned')
dl.plotting.img_show(plt.gca(), cleaned)

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的denoising_images.ipynb文件中。

另见

从图像中提取斑块

图像分割是将图像分割成多个片段的过程。这些片段具有相似的颜色或强度。这些片段在医学、交通、天文学或其他方面通常也有的含义。

分割图像最简单的方法是使用阈值，这将产生两个片段(如果值等于阈值，我们将它们放在两个片段之一中)。大津阈值方法最小化两个片段的加权方差(参考以下等式):

Extracting patches from an image

如果我们分割图像，去除噪声或外来伪像是个好主意。借助膨胀(参见也参见部分)，我们可以找到图像中属于背景和前景的部分。然而，膨胀给我们留下了未识别的像素。

做好准备

按照中的说明设置 OpenCV 。

怎么做...

  1. 进口情况如下:

    import numpy as np
import cv2
from matplotlib import pyplot as plt
from sklearn.datasets import load_sample_image
import dautil as dl
from IPython.display import HTML

  2. 将原始图像绘制如下:

    sp = dl.plotting.Subplotter(2, 2, context)
img = load_sample_image('flower.jpg')
dl.plotting.img_show(sp.ax, img)
sp.label()

  3. 如下绘制大津阈值图像:

    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
_, thresh = cv2.threshold(gray, 0, 255,
                          cv2.THRESH_OTSU)

dl.plotting.img_show(sp.next_ax(), thresh)
sp.label()

  4. 绘制前景和背景分散的图像，如下所示:

    kernel = np.ones((3, 3), np.uint8)
opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN,
                           kernel, iterations=2)

bg = cv2.dilate(opening, kernel, iterations=3)

dist_transform = cv2.distanceTransform(opening, cv2.DIST_L2, 5)
_, fg = cv2.threshold(dist_transform, 0.7 * dist_transform.max(),
                      255, 0)

fg = np.uint8(fg)
rest = cv2.subtract(bg, fg)

dl.plotting.img_show(sp.next_ax(), rest)
sp.label()

  5. 用标记绘制图像，如下所示:

    _, markers = cv2.connectedComponents(fg)
markers += 1
markers[rest == 255] = 0

dl.plotting.img_show(sp.next_ax(), markers)
sp.label()

HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的 extracting_patches.ipynb文件中。

另见

用哈尔级联检测人脸

人脸是人体解剖学的一个识别特征。严格来说，许多动物也有脸，但这与大多数实际应用不太相关。人脸检测试图在图像中找到代表人脸的(矩形)区域。人脸检测是物体检测的一种，因为人脸是物体的一种。

大多数人脸检测算法都擅长检测干净的人脸，因为大多数训练图像都属于这一类。倾斜的面部、明亮的灯光或嘈杂的图像可能会导致面部检测出现问题。从一张脸推断年龄、性别或种族(例如内眦赘皮的存在)是可能的，这当然对营销有用。

一个可能的应用是分析社交媒体网站上的个人资料图片。OpenCV 使用基于哈尔特征的级联分类器系统来检测人脸。该系统还被命名为“紫百合-琼斯目标检测框架”T4，以其在 2001 年提出的“T5”清单命名。

该算法具有以下步骤:

  1. 哈尔特征选择:哈尔特征类似于哈尔小波(如第 6 章、信号处理和时间序列中应用离散小波变换配方所述)。
  2. 创建完整图像(参考使用 SURF 配方检测特征)。
  3. Adaboost 训练(参考第九章、集成学习和降维中 Boosting 更好的学习食谱)。
  4. 级联分类器。

当我们看人脸图像时，我们可以创建与亮度相关的试探法。

例如，鼻子区域比直接位于其左右两侧的区域更亮。因此，我们可以定义一个覆盖鼻子的白色矩形和覆盖邻近区域的黑色矩形。当然 Viola-Jones 系统并不知道鼻子的确切位置，但是通过定义不同大小的窗口并寻找相应的白色和黑色矩形，就有机会匹配到鼻子。实际的哈尔特征被定义为黑色矩形的亮度总和和相邻矩形的亮度总和。对于 24 x 24 的窗口，我们有超过 16 万个功能(大约是 24 的四次方)。

训练集由大量正面(有脸)图像和负面(没有脸)图像组成。只有大约 0.01%的窗口(大约 24×24 像素)实际包含人脸。级联分类器逐步过滤掉负像区域。在每个渐进阶段，分类器在更少的图像窗口上逐渐使用更多的特征。这个想法是把大部分时间花在包含人脸的图像块上。维奥拉和琼斯的原始论文有 38 个阶段，前五个阶段有 1、10、25、25 和 50 个特征。平均每个图像窗口评估 10 个特征。

在 OpenCV 中，可以自己训练一个级联分类器，如中所述。然而，OpenCV 有预先训练好的人脸、眼睛和其他特征的分类器。这些分类器的配置存储为 XML 文件，可以在安装 OpenCV 的文件夹中找到(在我的机器上，/usr/local/share/OpenCV/haarcascades /)。

做好准备

按照中的说明设置 OpenCV 。

怎么做...

  1. 进口情况如下:

    import cv2
from scipy.misc import lena
import matplotlib.pyplot as plt
import numpy as np
import dautil as dl
import os
from IPython.display import HTML

  2. 定义以下功能，绘制检测到人脸的图像(如果检测到):

    def plot_with_rect(ax, img):
    img2 = img.copy()

    for x, y, w, h in face_cascade.detectMultiScale(img2, 1.3, 5):
        cv2.rectangle(img2, (x, y), (x + w, y + h), (255, 0, 0), 2)

    dl.plotting.img_show(ax, img2, cmap=plt.cm.gray)

  3. 下载 XML 配置文件，创建分类器:

    # dir = '/usr/local/share/OpenCV/haarcascades/'
base = 'https://raw.githubusercontent.com/Itseez/opencv/master/data/'
url = base + 'haarcascades/haarcascade_frontalface_default.xml'
path = os.path.join(dl.data.get_data_dir(),
                    'haarcascade_frontalface_default.xml')

if not dl.conf.file_exists(path):
    dl.data.download(url, path)

face_cascade = cv2.CascadeClassifier(path)

  4. 用检测到的人脸绘制原始图像:

    sp = dl.plotting.Subplotter(2, 2, context)
img = lena().astype(np.uint8)
plot_with_rect(sp.ax, img)
sp.label()

  5. 绘制轻微旋转的图像(检测失败):

    rows, cols = img.shape
mat = cv2.getRotationMatrix2D((cols/2, rows/2), 21, 1)
rot = cv2.warpAffine(img, mat, (cols, rows))
plot_with_rect(sp.next_ax(), rot)
sp.label()

  6. 绘制添加了噪声的图像(检测失败):

    np.random.seed(36)
noisy = img * (np.random.random(img.shape) < 0.6)
plot_with_rect(sp.next_ax(), noisy)
sp.label()

  7. 用检测到的人脸绘制模糊图像:

    blur = cv2.blur(img, (9, 9))
plot_with_rect(sp.next_ax(), blur)
sp.label()

HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的detecting_faces.ipynb文件中。

另见

寻找明亮的星星

许多星星在晚上是可见的，甚至不用望远镜或任何其他光学设备。一般来说，恒星比行星地球大，但在它们进化的某些阶段，它们可以更小。由于距离很远，它们看起来像小点。通常，这些点由两颗(双星系统)或更多恒星组成。不是所有的恒星都发出可见光，也不是所有的星光都能到达我们这里。

我们有许多方法可以在星空图像中找到明亮的星星。在这个食谱中，我们将寻找亮度的局部最大值，它也高于一个阈值。为了确定亮度，我们将把图像转换到 HSV 颜色空间。在这个颜色空间中，三个维度是色调、饱和度和值(亮度)。OpenCV split()将颜色空间中的图像值转换为组成值，例如色调、饱和度和亮度。这是一个相对缓慢的操作。为了找到最大值，我们可以应用 SciPy argrelmax()函数。

做好准备

按照中的说明设置 OpenCV 配方。

怎么做...

  1. 进口情况如下:

    import dautil as dl
import os
import cv2
import matplotlib.pyplot as plt
from scipy.signal import argrelmax
import numpy as np
from IPython.display import HTML

  2. 定义以下功能来扫描水平或垂直轴上的局部亮度峰值:

    def scan_axis(v, axis):
    argmax = argrelmax(v, order=int(np.sqrt(v.shape[axis])),
                       axis=axis)

    return set([(i[0], i[1]) for i in np.column_stack(argmax)])

  3. 下载图片分析:

    dir = dl.data.get_data_dir()
path = os.path.join(dir, 'night-927168_640.jpg')
base = 'https://pixabay.com/static/uploads/photo/2015/09/06/10/19/'
url = base + 'night-927168_640.jpg'

if not dl.conf.file_exists(path):
    dl.data.download(url, path)

  4. 从图像中提取亮度值:

    img = cv2.imread(path)
hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)

h, s, v = cv2.split(hsv)

# Transform for normalization
v = v.astype(np.uint16) ** 2

  5. 绘制亮度值直方图:

    sp = dl.plotting.Subplotter(2, 2, context)
sp.ax.hist(v.ravel(), normed=True)
sp.label()

  6. 绘制轴 0 的亮度值直方图:

    dl.plotting.hist_norm_pdf(sp.next_ax(), v.mean(axis=0))
sp.label()

  7. 绘制轴 1 的亮度值直方图:

    dl.plotting.hist_norm_pdf(sp.next_ax(), v.mean(axis=1))
sp.label()

  8. 用我们认为包含亮星的点绘制图像:

    points = scan_axis(v, 0).intersection(scan_axis(v, 1))

limit = np.percentile(np.unique(v.ravel()), 95)

kp = [cv2.KeyPoint(p[1], p[0], 1) for p in points
      if v[p[0], p[1]] > limit]
with_kp = cv2.drawKeypoints(img, kp, None, (255, 0, 0))

dl.plotting.img_show(sp.next_ax(), with_kp)
sp.label()

HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的 searching_stars.ipynb文件中。

另见

从图像中提取元数据

数字照片经常包含额外的文本元数据，例如时间戳、曝光信息和地理位置。其中一些元数据可由相机所有者编辑。例如，在营销环境中，从社交媒体网站上的个人资料(或其他)图像中提取元数据可能会很有用。据称，告密者爱德华·斯诺登声称美国国家安全局正在从全球在线数据中收集 EXIF 元数据。

做好准备

在这个食谱中，我们将使用 ExifRead 提取 EXIF 元数据。

按照以下步骤安装 ExifRead:

$ pip install ExifRead

我用 ExifRead 2.1.2 测试了代码。

怎么做...

  1. 进口情况如下:

    import exifread
import pprint

  2. 打开图像如下:

    f = open('covers.jpg', 'rb')

  3. 打印标签和密钥如下:

    # Return Exif tags
tags = exifread.process_file(f)
print(tags.keys())
pprint.pprint(tags)
f.close()

参考以下最终结果:

dict_keys(['EXIF Flash', 'Image Make', 'EXIF Contrast', 'EXIF DateTimeOriginal', 'Image ResolutionUnit', 'EXIF ComponentsConfiguration', 'EXIF ISOSpeedRatings', 'Image ExifOffset', 'Image ImageDescription', 'EXIF MaxApertureValue', 'EXIF ExposureBiasValue', 'Image YResolution', 'Image Orientation', 'EXIF DateTimeDigitized', 'EXIF MeteringMode', 'EXIF Sharpness', 'EXIF WhiteBalance', 'EXIF ExposureTime', 'Image Model', 'EXIF SceneCaptureType', 'Image Software', 'EXIF SceneType', 'EXIF SubjectDistanceRange', 'EXIF LightSource', 'EXIF FocalLengthIn35mmFilm', 'Image XResolution', 'Image DateTime', 'EXIF FileSource', 'EXIF ExposureProgram', 'EXIF FocalLength', 'EXIF FNumber', 'EXIF Saturation', 'EXIF ExifImageWidth', 'EXIF ExposureMode', 'EXIF DigitalZoomRatio', 'EXIF FlashPixVersion', 'EXIF ExifVersion', 'EXIF ColorSpace', 'EXIF CustomRendered', 'EXIF GainControl', 'EXIF CompressedBitsPerPixel', 'EXIF ExifImageLength'])
{'EXIF ColorSpace': (0xA001) Short=sRGB @ 406,
 'EXIF ComponentsConfiguration': (0x9101) Undefined=YCbCr @ 298,
 'EXIF CompressedBitsPerPixel': (0x9102) Ratio=2 @ 650,
 'EXIF Contrast': (0xA408) Short=Normal @ 550,
 'EXIF CustomRendered': (0xA401) Short=Normal @ 466,
 'EXIF DateTimeDigitized': (0x9004) ASCII=0000:00:00 00:00:00 @ 630,
 'EXIF DateTimeOriginal': (0x9003) ASCII=0000:00:00 00:00:00 @ 610,
 'EXIF DigitalZoomRatio': (0xA404) Ratio=0 @ 682,
 'EXIF ExifImageLength': (0xA003) Long=240 @ 430,
 'EXIF ExifImageWidth': (0xA002) Long=940 @ 418,
 'EXIF ExifVersion': (0x9000) Undefined=0220 @ 262,
 'EXIF ExposureBiasValue': (0x9204) Signed Ratio=0 @ 658,
 'EXIF ExposureMode': (0xA402) Short=Auto Exposure @ 478,
 'EXIF ExposureProgram': (0x8822) Short=Program Normal @ 238,
 'EXIF ExposureTime': (0x829A) Ratio=10/601 @ 594,
 'EXIF FNumber': (0x829D) Ratio=14/5 @ 602,
 'EXIF FileSource': (0xA300) Undefined=Digital Camera @ 442,
 'EXIF Flash': (0x9209) Short=Flash fired, auto mode @ 370,
 'EXIF FlashPixVersion': (0xA000) Undefined=0100 @ 394,
 'EXIF FocalLength': (0x920A) Ratio=39/5 @ 674,
 'EXIF FocalLengthIn35mmFilm': (0xA405) Short=38 @ 514,
 'EXIF GainControl': (0xA407) Short=None @ 538,
 'EXIF ISOSpeedRatings': (0x8827) Short=50 @ 250,
 'EXIF LightSource': (0x9208) Short=Unknown @ 358,
 'EXIF MaxApertureValue': (0x9205) Ratio=3 @ 666,
 'EXIF MeteringMode': (0x9207) Short=Pattern @ 346,
 'EXIF Saturation': (0xA409) Short=Normal @ 562,
 'EXIF SceneCaptureType': (0xA406) Short=Standard @ 526,
 'EXIF SceneType': (0xA301) Undefined=Directly Photographed @ 454,
 'EXIF Sharpness': (0xA40A) Short=Normal @ 574,
 'EXIF SubjectDistanceRange': (0xA40C) Short=0 @ 586,
 'EXIF WhiteBalance': (0xA403) Short=Auto @ 490,
 'Image DateTime': (0x0132) ASCII=0000:00:00 00:00:00 @ 184,
 'Image ExifOffset': (0x8769) Long=204 @ 126,
 'Image ImageDescription': (0x010E) ASCII=           @ 134,
 'Image Make': (0x010F) ASCII=NIKON @ 146,
 'Image Model': (0x0110) ASCII=E7900 @ 152,
 'Image Orientation': (0x0112) Short=Horizontal (normal) @ 54,
 'Image ResolutionUnit': (0x0128) Short=Pixels/Inch @ 90,
 'Image Software': (0x0131) ASCII=E7900v1.1 @ 174,
 'Image XResolution': (0x011A) Ratio=300 @ 158,
 'Image YResolution': (0x011B) Ratio=300 @ 166}

代码在本书代码包的img_metadata.py文件中。

另见

从图像中提取纹理特征

纹理是图像的空间和视觉质量。在这个食谱中，我们将看看 哈拉利克纹理特征。这些特征基于如下定义的共现矩阵 (11.5) :

Extracting texture features from images

在等式 11.5 中， i 和 j 是强度，而 p 和 q 是位置。哈拉里克特征是从共现矩阵中导出的 13 个度量，其中一些在等式 11.6 中给出。更多完整列表，请参考http://murphylab . web . CMU . edu/publications/Boland/Boland _ node 26 . html(2015 年 12 月检索)。

我们将使用 mahotas API 计算 Haralick 特征，并将其应用于 scikit-learn 的手写数字数据集。

做好准备

按照以下步骤安装 maho tas:

$ pip install mahotas

我用 mahotas 1.4.0 测试了代码。

怎么做...

  1. 进口如下:

    import mahotas as mh
import numpy as np
from sklearn.datasets import load_digits
import matplotlib.pyplot as plt
from tpot import TPOT
from sklearn.cross_validation import train_test_split
import dautil as dl

  2. 加载 sci kit-学习数字数据如下:

    digits = load_digits()
X = digits.data.copy()

  3. 创建哈拉利克特征并添加它们:

    for i, img in enumerate(digits.images):
    np.append(X[i], mh.features.haralick(
        img.astype(np.uint8)).ravel())

  4. 用 TPOT(或我的叉子，如第 9 章、集成学习和降维所述)拟合模型并评分:

    X_train, X_test, y_train, y_test = train_test_split(
    X, digits.target, train_size=0.75)

tpot = TPOT(generations=6, population_size=101,
            random_state=46, verbosity=2)
tpot.fit(X_train, y_train)

print('Score {:.2f}'.format(tpot.score(X_train, y_train, X_test, y_test)))

  5. 将第一张原图绘制如下:

    dl.plotting.img_show(plt.gca(), digits.images[0])
plt.title('Original Image')

  6. 绘制该图像的核心特征:

    plt.figure()
dl.plotting.img_show(plt.gca(), digits.data[0].reshape((8, 8)))
plt.title('Core Features')

  7. 绘制该图像的哈拉利克特征:

    plt.figure()
dl.plotting.img_show(plt.gca(), mh.features.haralick(
    digits.images[0].astype(np.uint8)))
plt.title('Haralick Features')

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的extracting_texture.ipynb文件中。

另见

对图像应用层次聚类

我们在第 9 章、集成学习和降维中遇到了层次聚类的概念。在这个食谱中，我们将通过分层聚类来分割图像。我们将应用凝聚聚类 O(n 3 ，这是的一种层次聚类。

在聚集聚类中，每个项目在初始化时被分配其自己的聚类。随后，这些集群合并(聚集)并根据需要向上移动。显然，我们只合并在某种程度上相似的集群。

初始化之后，我们找到距离最近的一对，并将其合并。合并后的集群是由较低级别的集群组成的较高级别的集群。之后，我们再次找到最接近的一对并将其合并，以此类推。在此过程中，集群可以有任意数量的项目。当我们达到一定数量的集群后，或者当集群相距太远时，我们停止集群。

怎么做...

  1. 进口情况如下:

    import numpy as np
from scipy.misc import ascent
import matplotlib.pyplot as plt
from sklearn.feature_extraction.image import grid_to_graph
from sklearn.cluster import AgglomerativeClustering
import dautil as dl

  2. 加载图像并将其加载到数组中:

    img = ascent()
X = np.reshape(img, (-1, 1))

  3. 将图像聚类，聚类数设置为9(猜测):

    connectivity = grid_to_graph(*img.shape)
NCLUSTERS = 9
ac = AgglomerativeClustering(n_clusters=NCLUSTERS,
                             connectivity=connectivity)
ac.fit(X)
label = np.reshape(ac.labels_, img.shape)

  4. 绘制叠加了簇段的图像:

    for l in range(NCLUSTERS):
    plt.contour(label == l, contours=1,
                colors=[plt.cm.spectral(l/float(NCLUSTERS)), ])

dl.plotting.img_show(plt.gca(), img, cmap=plt.cm.gray)

最终结果参见下面的截图:

How to do it...

代码在本书代码包的 clustering_hierarchy.ipynb文件中。

另见

利用光谱聚类分割图像

光谱聚类是一种可以用来分割图像的聚类技术。scikit-learn spectral_clustering()函数实现了归一化图割谱聚类算法。该算法将图像表示为单位图。这里的“图”与第八章、文本挖掘和社交网络分析中的数学概念相同。该算法试图分割图像，同时最小化片段大小和沿切割的强度梯度比率。

怎么做...

  1. 进口情况如下:

    import numpy as np
import matplotlib.pyplot as plt
from sklearn.feature_extraction.image import img_to_graph
from sklearn.cluster import spectral_clustering
from sklearn.datasets import load_digits

  2. 加载数字数据集如下:

    digits = load_digits()
img = digits.images[0].astype(float)
mask = img.astype(bool)

  3. 根据图像创建图形:

    graph = img_to_graph(img, mask=mask)
graph.data = np.exp(-graph.data/graph.data.std())

  4. 应用谱聚类得到三个聚类:

    labels = spectral_clustering(graph, n_clusters=3)
label_im = -np.ones(mask.shape)
label_im[mask] = labels

  5. 将原始图像绘制如下:

    plt.matshow(img, False)
plt.gca().axis('off')
plt.title('Original')

  6. 将图像与三个簇绘制如下:

    plt.figure()
plt.matshow(label_im, False)
plt.gca().axis('off')
plt.title('Clustered')

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的clustering_spectral.ipynb文件中。

另见

十二、并行性和性能

在本章中，我们将介绍以下食谱:

  • 用 Numba 进行即时编译
  • 用 Numexpr 加速数值表达式
  • 使用threading模块运行多个线程
  • 使用concurrent.futures模块启动多个任务
  • 与asyncio模块异步访问资源
  • 带execnet的分布式处理
  • 分析内存使用情况
  • 即时计算平均值、方差、偏斜度和峰度
  • 使用最近最少使用的缓存进行缓存
  • 缓存 HTTP 请求
  • 使用最小计数草图进行流式计数
  • 利用 OpenCL 开发图形处理器的能力

简介

1943 年至 1946 年间建造的 ENIAC，用一万八千根管子填满了一个大房间，拥有 20 位内存。从那以后，我们走过了漫长的道路。摩尔定律也预测了这种指数增长。当然，我们面对的是一个自我实现的预言还是一个基本现象，这很难说。据称，增长开始减速。

鉴于我们目前对技术、热力学和量子力学的了解，我们可以为摩尔定律设定硬性限制。然而，我们的假设可能是错误的；例如，科学家和工程师可能会想出从根本上更好的技术来制造芯片。(其中一个发展就是量子计算，目前还远未普及。)最大的障碍是散热，散热通常以 kT 为单位测量，K 为玻尔兹曼常数(约 10-23 J/K)，而 T 以开尔文为单位(冰点为 273.15 K)。一个芯片的每比特散热量至少为kT(350K 时为 10-20 J)。90 年代的半导体至少消耗了十万 kT 。计算系统在运行过程中会经历能量水平的变化。能量上最小的容许差大概是 100 kT 。即使我们设法避免了这种限制，我们也很快会在接近原子水平的水平上运行，由于量子力学的原因，这是不切实际的(关于粒子的信息从根本上是有限的)，除非我们谈论的是量子计算机。目前的共识是，我们将在几十年内达到极限。另一个考虑是芯片的复杂布线。复杂的布线大大降低了芯片的预期寿命。

本章是关于软件性能的；然而，还有其他更重要的软件方面，例如可维护性、健壮性和可用性。押注摩尔定律是有风险的，也是不现实的，因为我们还有其他提高性能的可能性。第一种选择是使用多台机器、单台机器上的内核、GPU 或其他专用硬件(如)尽可能并行地完成工作。例如，我正在八核机器上测试代码。作为一名学生，我有幸参与了一个以创建网格为目标的项目。网格被认为是将大学计算机聚集到一个单一的计算环境中。在后期，也有计划连接其他计算机，有点像 SETI 项目。(众所周知，很多办公电脑在周末和晚上都是闲置的，为什么不让它们也工作呢？)

当然，目前有各种各样的商业云系统，比如亚马逊和谷歌提供的。我不会讨论这些，因为我觉得这些是更专业的主题，尽管我在 Python 数据分析中确实介绍了一些 Python 特有的云系统。

提高性能的第二种方法是应用缓存，从而避免不必要的函数调用。我在第 9 章、集成学习和降维中介绍了具有缓存功能的 joblib 库。Python 3 为我们带来了并行性和缓存的新特性。

第三种方法是靠近金属。众所周知，Python 是一种带有虚拟机和解释器的高级编程语言。Python 有一个额外的层，这是一种不同于 C 语言的语言。当我还是个学生的时候，我们被教导 C 是一种高级语言，以汇编和机器码为低级。据我所知，现在，几乎没有人用汇编语言编写代码。通过 Cython(包含在 Python 数据分析中)和类似的软件，我们可以编译我们的代码来获得与 C 和 C++相当的性能。编译是一件麻烦的事情，也是有问题的，因为它降低了平台依赖性带来的可移植性。一个常见的解决方案是使用 shell 脚本自动编译并生成文件。Numba 和其他类似的项目通过即时编译让生活变得更加容易，尽管有一些限制。

用 Numba 及时编译

Numba 软件使用特殊的函数装饰器执行即时编译。编译会自动生成本机机器代码。生成的代码可以在 CPU 和 GPU 上运行。Numba 的主要用例是使用 NumPy 数组的数学密集型代码。

我们可以用带有可选函数签名的@numba.jit装饰器(例如int32(int32))编译代码。这些类型对应于相似的 NumPy 类型。Numba 在nopython和object模式下运行。nopython模式速度更快，但限制更多。我们也可以通过nogil选项释放 全局解释器锁 ( GIL )。您可以通过使用cache参数请求文件缓存来缓存编译结果。

@vectorize装饰器将带有标量参数的函数转换成 NumPy ufuncs。向量化提供了额外的优势，例如自动广播，并且可以在单核、多核并行或 GPU 上使用。

做好准备

使用以下命令安装 Numba:

$ pip/conda install numba

我用 Numba 0.22.1 测试了代码。

怎么做...

  1. 进口情况如下:

    from numba import vectorize
from numba import jit
import numpy as np

  2. 定义以下功能来使用@vectorize装饰器:

    @vectorize
def vectorize_version(x, y, z):
    return x ** 2 + y ** 2 + z ** 2

  3. 定义以下功能来使用@jit装饰器:

    @jit(nopython=True)
def jit_version(x, y, z):
    return x ** 2 + y ** 2 + z ** 2

  4. 定义一些随机数组如下:

    np.random.seed(36)
x = np.random.random(1000)
y = np.random.random(1000)
z = np.random.random(1000)

  5. 测量数组平方和所需的时间:

    %timeit x ** 2 + y ** 2 + z ** 2
%timeit vectorize_version(x, y, z)
%timeit jit_version(x, y, z)
jit_version.inspect_types()

最终结果参见以下截图:

How to do it...

代码在本书代码包的compiling_numba.ipynb文件中。

它是如何工作的

在我的机器上测得的最佳时间是 1.82 微秒，这比正常 Python 代码测得的时间快得多。在截图的最后，我们看到了编译的结果，最后一部分被省略了，因为它太长了，很难阅读。我们会收到警告，这很可能是由 CPU 缓存引起的。我是故意留下的，但是你可以使用不适合缓存的更大的数组来删除它们。

另见

用 Numexpr 加速数值表达式

Numexpr 是一个数值数组表达式求值的软件包，在你安装 pandas 的时候也安装了，你可能在其他菜谱的水印里看到过公布(用 Numexpr 2.3.1 测试过)。Numexpr 试图通过避免创建临时变量来加快计算速度，因为读取变量可能是一个潜在的瓶颈。对于无法容纳在中央处理器高速缓存中的数组，预计加速比最大。

Numexpr 将大数组分割成块，这些块可以放入缓存中，并且它还在可能的情况下并行使用多个内核。它有一个evaluate()函数，该函数接受简单表达式并对其求值(有关支持特性的完整列表，请参考文档)。

怎么做...

  1. 进口情况如下:

    import numexpr as ne
import numpy as np

  2. 生成随机数组，数组应该太大，无法保存在缓存中:

    a = np.random.rand(1e6)
b = np.random.rand(1e6)

  3. 计算一个简单的算术表达式并测量执行时间:

    %timeit 2 * a ** 3 + 3 * b ** 9
%timeit ne.evaluate("2 * a ** 3 +3 * b ** 9 ")

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的speeding_numexpr.ipynb文件中。

它是如何工作的

我们生成了不应该放在缓存中的随机数据，以避免缓存效应，因为这是 Numexpr 的最佳用例。缓存的大小因机器而异，因此如有必要，请为数组使用更大或更小的大小。在示例中，我们放入了一个包含简单算术表达式的字符串，尽管我们可以使用稍微复杂一点的表达式。有关更多详细信息，请参考文档。我用一台有八个内核的机器测试了代码。加速比大于八倍，所以这显然是由于 Numexpr。

另见

用线程模块运行多个线程

计算机进程是正在运行的程序的一个实例。进程实际上是重量级的，所以我们可能更喜欢线程，线程更轻。事实上，线程通常只是进程的子单元。进程是相互分离的，而线程可以共享指令和数据。

操作系统通常为每个内核分配一个线程(如果有多个)，或者定期在线程之间切换；这叫 时间切片。线程作为进程可以有不同的优先级，操作系统有后台运行的守护线程，优先级非常低。

线程之间的切换比进程之间的切换更容易；但是，由于线程共享信息，因此使用起来更加危险。例如，如果多个线程能够同时增加一个计数器，这将使代码不确定，并且可能不正确。最小化风险的一种方法是确保一次只有一个线程可以访问共享变量或共享函数。这个策略在 Python 中实现为 GIL。

怎么做...

  1. 进口情况如下:

    import dautil as dl
import ch12util
from functools import partial
from queue import Queue
from threading import Thread
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import skew
from IPython.display import HTML

STATS = []

  2. 定义以下功能重新采样:

    def resample(arr):
    sample = ch12util.bootstrap(arr)
    STATS.append((sample.mean(), sample.std(), skew(sample)))

  3. 定义以下要引导的类:

    class Bootstrapper(Thread):
    def __init__(self, queue, data):
        Thread.__init__(self)
        self.queue = queue
        self.data = data
        self.log = dl.log_api.conf_logger(__name__)

    def run(self):
        while True:
            index = self.queue.get()

            if index % 10 == 0:
                self.log.debug('Bootstrap {}'.format(
                    index))

            resample(self.data)
            self.queue.task_done()

  4. 定义以下函数来执行串行重采样:

    def serial(arr, n):
    for i in range(n):
        resample(arr)

  5. 定义以下函数执行并行重采样:

    def threaded(arr, n):
    queue = Queue()

    for x in range(8):
        worker = Bootstrapper(queue, arr)
        worker.daemon = True
        worker.start()

    for i in range(n):
        queue.put(i)

    queue.join()

  6. 绘制时刻分布和执行时间:

    sp = dl.plotting.Subplotter(2, 2, context)
temp = dl.data.Weather.load()['TEMP'].dropna().values
np.random.seed(26)
threaded_times = ch12util.time_many(partial(threaded, temp))
serial_times = ch12util.time_many(partial(serial, temp))

ch12util.plot_times(sp.ax, serial_times, threaded_times)

stats_arr = np.array(STATS)
ch12util.plot_distro(sp.next_ax(), stats_arr.T[0], temp.mean())
sp.label()

ch12util.plot_distro(sp.next_ax(), stats_arr.T[1], temp.std())
sp.label()

ch12util.plot_distro(sp.next_ax(), stats_arr.T[2], skew(temp))
sp.label()

HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书的代码包中的 running_threads.ipynb文件中。

另见

用并发期货模块启动多个任务

concurrent.futures模块是一个 Python 模块，我们可以用它异步执行可调用的东西。如果您熟悉 Java 并浏览该模块，您会注意到与等效 Java API 的一些相似之处，例如类名和体系结构。根据 Python 文档，这不是巧合。

在这种情况下，任务是一个独立的工作单元。例如，打印文档可以被认为是一项任务，但通常我们会考虑更小的任务，例如添加两个数字。

怎么做...

  1. 进口情况如下:

    import dautil as dl
import ch12util
from functools import partial
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import skew
import concurrent.futures
from IPython.display import HTML

STATS = []

  2. 定义以下函数进行重采样:

    def resample(arr):
    sample = ch12util.bootstrap(arr)
    STATS.append((sample.mean(), sample.std(), skew(sample)))

  3. 定义以下要引导的类:

    class Bootstrapper():
    def __init__(self, data):
        self.data = data
        self.log = dl.log_api.conf_logger(__name__)

    def run(self, index):
        if index % 10 == 0:
            self.log.debug('Bootstrap {}'.format(
                index))

        resample(self.data)

  4. 定义以下函数来执行串行重采样:

    def serial(arr, n):
    for i in range(n):
        resample(arr)

  5. 定义以下函数执行并行重采样:

    def parallel(arr, n):
    executor = concurrent.futures.ThreadPoolExecutor(max_workers=8)
    bootstrapper = Bootstrapper(arr)

    for x in executor.map(bootstrapper.run, range(n)):
        pass

    executor.shutdown()

  6. 绘制时刻分布和执行时间:

    rain = dl.data.Weather.load()['RAIN'].dropna().values
np.random.seed(33)
parallel_times = ch12util.time_many(partial(parallel, rain))
serial_times = ch12util.time_many(partial(serial, rain))

sp = dl.plotting.Subplotter(2, 2, context)
ch12util.plot_times(sp.ax, serial_times, parallel_times)

STATS = np.array(STATS)
ch12util.plot_distro(sp.next_ax(), STATS.T[0], rain.mean())
sp.label()

ch12util.plot_distro(sp.next_ax(), STATS.T[1], rain.std())
sp.label()

ch12util.plot_distro(sp.next_ax(), STATS.T[2], skew(rain))
sp.label()
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书的代码包中的 launching_futures.ipynb文件中。

另见

与异步模块异步访问资源

输入/输出(例如，文件或数据库访问)缓慢是生活中的一个基本事实。I/O 不仅慢，而且不可预测。在一个常见的场景中，我们等待数据(来自 web 服务或传感器)并将数据写入文件系统或数据库。在这种情况下，我们会发现自己受到输入/输出的限制，花费更多的时间等待数据，而不是实际处理数据。我们可以定期轮询数据或触发事件(检查您的手表或设置警报)。图形用户界面通常有特殊的线程，在无限循环中等待用户输入。

异步输入输出的 Python asyncio模块使用了的概念，并通过相关的函数装饰器来协调。在第 5 章、网络挖掘、数据库和大数据的刮网配方中也给出了该模块的一个简单示例。子程序可以看作是协同程序的一个特例。一个子例程有一个开始和结束点，要么通过一个返回语句提前结束，要么到达子例程定义的末尾。相比之下，一个协同程序可以通过调用另一个协同程序，然后从该退出点恢复执行，来产生yield from语句。可以说，科罗廷正在让另一个科罗廷接管，并回到睡眠状态，直到它再次被激活。

子程序可以放在一个堆栈上。然而，协同程序需要多个堆栈，这使得理解代码和潜在的异常更加复杂。

怎么做...

代码在本书代码包的accessing_asyncio.ipynb文件中:

  1. 进口情况如下:

    import dautil as dl
import ch12util
from functools import partial
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import skew
import asyncio
import time
from IPython.display import HTML

STATS = []

  2. 定义以下函数进行重采样:

    def resample(arr):
    sample = ch12util.bootstrap(arr)
    STATS.append((sample.mean(), sample.std(), skew(sample)))

  3. 定义以下要引导的类:

    class Bootstrapper():
    def __init__(self, data, queue):
        self.data = data
        self.log = dl.log_api.conf_logger(__name__)
        self.queue = queue

    @asyncio.coroutine
    def run(self):
        while not self.queue.empty():
            index = yield from self.queue.get()

            if index % 10 == 0:
                self.log.debug('Bootstrap {}'.format(
                    index))

            resample(self.data)
            # simulates slow IO
            yield from asyncio.sleep(0.01)

  4. 定义以下功能执行串行重采样:

    def serial(arr, n):
    for i in range(n):
        resample(arr)
        # simulates slow IO
        time.sleep(0.01)

  5. 定义以下函数执行并行重采样:

    def parallel(arr, n):
    q = asyncio.Queue()

    for i in range(n):
        q.put_nowait(i)

    bootstrapper = Bootstrapper(arr, q)
    policy = asyncio.get_event_loop_policy()
    policy.set_event_loop(policy.new_event_loop())
    loop = asyncio.get_event_loop()

    tasks = [asyncio.async(bootstrapper.run())
             for i in range(n)]

    loop.run_until_complete(asyncio.wait(tasks))
    loop.close()

  6. 绘制力矩和执行时间的分布:

    pressure = dl.data.Weather.load()['PRESSURE'].dropna().values
np.random.seed(33)
parallel_times = ch12util.time_many(partial(parallel, pressure))
serial_times = ch12util.time_many(partial(serial, pressure))

dl.options.mimic_seaborn()
ch12util.plot_times(plt.gca(), serial_times, parallel_times)

sp = dl.plotting.Subplotter(2, 2, context)
ch12util.plot_times(sp.ax, serial_times, parallel_times)

STATS = np.array(STATS)
ch12util.plot_distro(sp.next_ax(), STATS.T[0], pressure.mean())
sp.label()

ch12util.plot_distro(sp.next_ax(), STATS.T[1], pressure.std())
sp.label()

ch12util.plot_distro(sp.next_ax(), STATS.T[2], skew(pressure))
sp.label()
HTML(sp.exit())

最终结果参见下面的截图:

How to do it...

另见

使用 execnet 进行分布式处理

execnet模块采用无共享模式，使用通道进行通信。这个上下文中的通道是用于在(分布式)计算机进程之间发送和接收消息的软件抽象。execnet最适用于将异构计算环境与不同的 Python 解释器和安装的软件相结合。这些环境可以有不同的操作系统和 Python 实现(CPython、Jython、pypypy 或其他)。

在 无共享架构中，计算节点不共享内存或文件。因此，该架构是完全分散的，具有完全独立的节点。明显的优势是我们不依赖任何一个节点。

做好准备

使用以下命令安装 execnet:

$ pip/conda install execnet

我用 execnet 1.3.0 测试了代码。

怎么做...

  1. 进口情况如下:

    import dautil as dl
import ch12util
from functools import partial
import matplotlib.pyplot as plt
import numpy as np
import execnet

STATS = []

  2. 定义以下助手函数:

    def run(channel, data=[]):
    while not channel.isclosed():
        index = channel.receive()

        if index % 10 == 0:
            print('Bootstrap {}'.format(
                index))

        total = 0

        for x in data:
            total += x

        channel.send((total - data[index])/(len(data) - 1))

  3. 定义以下函数来执行串行重采样:

    def serial(arr, n):
    for i in range(n):
        total = 0

        for x in arr:
            total += x

        STATS.append((total - arr[i])/(len(arr) - 1))

  4. 定义以下函数执行并行重采样:

    def parallel(arr, n):
    gw = execnet.makegateway()
    channel = gw.remote_exec(run, data=arr.tolist())

    for i in range(n):
        channel.send(i)
        STATS.append(channel.receive())

    gw.exit()

  5. 绘制平均值和执行时间的分布:

    ws = dl.data.Weather.load()['WIND_SPEED'].dropna().values
np.random.seed(33)
parallel_times = ch12util.time_many(partial(parallel, ws))
serial_times = ch12util.time_many(partial(serial, ws))

%matplotlib inline
dl.options.mimic_seaborn()
ch12util.plot_times(plt.gca(), serial_times, parallel_times)
plt.legend(loc='best')

plt.figure()
STATS = np.array(STATS)
ch12util.plot_distro(plt.gca(), STATS, ws.mean())
plt.title('Distribution of the Means')
plt.legend(loc='best')

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在书的代码包中的distributing_execnet.ipynb文件中。

另见

分析内存使用情况

在 Python 数据分析中，我们使用了各种分析工具。这些工具主要与测量执行时间有关。然而，记忆也很重要，尤其是当我们没有足够的记忆时。内存泄漏 是计算机程序的常见问题，我们可以通过执行内存分析来发现。当我们不释放不需要的内存时，就会发生泄漏。当我们使用需要比我们需要更多内存的数据类型时，也可能会出现问题，例如，整数数组就可以使用 NumPy float64数组。

Python memory_profiler模块可以逐行剖析代码的内存使用情况。安装后，您还可以通过各种神奇的命令在 IPython 笔记本中使用该模块。该模块通过与操作系统通信来工作。在 Windows 上，您将需要 Python psutil包进行通信。

做好准备

用以下命令安装memory_profiler:

$ pip install memory-profiler

我用 memory_profiler 0.39 测试了代码。

创建一个脚本进行分析(参见本书代码包中的mem_test.py文件):

import numpy as np

def test_me():
    a = np.random.random((999, 99))
    b = np.random.random((99, 99))
    a.ravel()
    b.tolist()

怎么做...

  1. 进口情况如下:

    import dautil as dl
from mem_test import test_me

  2. 按如下方式加载 IPython 扩展:

    %load_ext memory_profiler

  3. 使用以下命令逐行剖析测试脚本:

    %mprun -f test_me test_me()

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书的代码包中的 profiling_memory.ipynb文件中。

另见

即时计算平均值、方差、偏斜度和峰度

均值、方差、偏斜度和峰度是统计学中的重要量。一些计算涉及平方的和，对于大的值可能导致溢出。为了避免精度的损失，我们必须认识到方差在偏移某个常数时是不变的。

当我们在内存中有足够的空间时，我们可以直接计算力矩，必要时考虑数值问题。然而，我们可能不想将数据保存在内存中，因为它很多，或者因为它更便于计算飞行中的时刻。

一个在线和数值稳定的算法来计算方差已经提供了特里贝里(特里贝里，蒂莫西(2007)，计算高阶矩在线)。我们将把这个算法与LiveStats模块中的实现进行比较，尽管它不是最好的算法。如果你对改进的算法感兴趣，看看中列出的维基百科页面，也可以参见部分。

看看下面的等式:

Calculating the mean, variance, skewness, and kurtosis on the fly

偏斜度由 12.6 给出，峰度由 12.7 给出。

做好准备

使用以下命令安装 LiveStats:

$ pip install LiveStats

我用 LiveStats 1.0 测试了代码。

怎么做...

  1. 进口情况如下:

    from livestats import livestats
from math import sqrt
import dautil as dl
import numpy as np
from scipy.stats import skew
from scipy.stats import kurtosis
import matplotlib.pyplot as plt

  2. 定义以下函数来实现力矩计算的方程:

    # From https://en.wikipedia.org/wiki/
# Algorithms_for_calculating_variance
def online_kurtosis(data):
    n = 0
    mean = 0
    M2 = 0
    M3 = 0
    M4 = 0
    stats = []

    for x in data:
        n1 = n
        n = n + 1
        delta = x - mean
        delta_n = delta / n
        delta_n2 = delta_n ** 2
        term1 = delta * delta_n * n1
        mean = mean + delta_n
        M4 = M4 + term1 * delta_n2 * (n**2 - 3*n + 3) + \
            6 * delta_n2 * M2 - 4 * delta_n * M3
        M3 = M3 + term1 * delta_n * (n - 2) - 3 * delta_n * M2
        M2 = M2 + term1
        s = sqrt(n) * M3 / sqrt(M2 ** 3)
        k = (n*M4) / (M2**2) - 3
        stats.append((mean, sqrt(M2/(n - 1)), s, k))

    return np.array(stats)

  3. 初始化并加载数据如下:

    test = livestats.LiveStats([0.25, 0.5, 0.75])

data = dl.data.Weather.load()['TEMP'].\
    resample('M').dropna().values

  4. 用 LiveStats 计算各种统计，上节提到的算法，和对比一下我们一次对所有数据应用 NumPy 函数时的结果:

    ls = []
truth = []

test.add(data[0])

for i in range(1, len(data)):
    test.add(data[i])
    q1, q2, q3 = test.quantiles()

    ls.append((test.mean(), sqrt(test.variance()),
              test.skewness(), test.kurtosis(), q1[1], q2[1], q3[1]))
    slice = data[:i]
    truth.append((slice.mean(), slice.std(),
                  skew(slice), kurtosis(slice),
                  np.percentile(slice, 25), np.median(slice),
                  np.percentile(slice, 75)))

ls = np.array(ls)
truth = np.array(truth)
ok = online_kurtosis(data)

  5. 将结果绘制为如下:

    dl.options.mimic_seaborn()
cp = dl.plotting.CyclePlotter(plt.gca())
cp.plot(ls.T[0], label='LiveStats')
cp.plot(truth.T[0], label='Truth')
cp.plot(data)
plt.title('Live Stats Means')
plt.xlabel('# points')
plt.ylabel('Mean')
plt.legend(loc='best')

plt.figure()

mses = [dl.stats.mse(truth.T[i], ls.T[i])
        for i in range(7)]
mses.extend([dl.stats.mse(truth.T[i], ok[1:].T[i])
             for i in range(4)])
dl.plotting.bar(plt.gca(),
                ['mean', 'std', 'skew', 'kurt',
                 'q1', 'q2', 'q3',
                 'my_mean', 'my_std', 'my_skew', 'my_kurt'], mses)
plt.title('MSEs for Various Statistics')
plt.ylabel('MSE')

最终结果参见以下截图:

How to do it...

代码在本书的代码包中的 calculating_moments.ipynb文件中。

另见

使用最近最少使用的缓存进行缓存

缓存包括将结果(通常来自函数调用)存储在内存或磁盘上。如果操作正确，缓存有助于减少函数调用的次数。一般来说，出于空间原因，我们希望保持缓存较小。如果我们能够在缓存中找到项目，我们就谈论命中；否则，我们会有失误。显然，我们希望有尽可能多的命中和尽可能少的失误。这意味着我们希望最大化命中率。

缓存算法有很多，其中我们将介绍最近最少使用的**(LRU)算法。该算法跟踪缓存项何时被使用。如果缓存即将超过其最大指定大小，LRU 会删除最近最少使用的项目。理由是这些项目可能更古老，因此不再相关。LRU 有几种变体。其他算法则相反——删除最近的项目、最不常用的项目或随机项目。**

标准的 Python 库有一个 LRU 的实现，但也有一个专门的 Python 库，其中一些部分用 C 实现，因此它可能更快。我们将使用 NLTK lemmatize()方法比较两种实现(参考第 8 章、文本挖掘和社交网络分析中的词干、引理、过滤和 TF-IDF 评分方法)。

做好准备

按照以下步骤安装 fastcache:

$ pip/conda install fastcache

我用 fastcache 1.0.2 测试了代码。

怎么做...

  1. 进口情况如下:

    from fastcache import clru_cache
from functools import lru_cache
from nltk.corpus import brown
from nltk.stem import WordNetLemmatizer
import dautil as dl
import numpy as np
from IPython.display import HTML

  2. 定义以下函数进行缓存:

    def lemmatize(word, lemmatizer):
    return lemmatizer.lemmatize(word.lower())

  3. 定义以下函数来衡量缓存的效果:

    def measure(impl, words, lemmatizer):
    cache = dl.perf.LRUCache(impl, lemmatize)
    times = []
    hm = []

    for i in range(5, 12):
        cache.maxsize = 2 ** i
        cache.cache()
        with dl.perf.StopWatch() as sw:
            _ = [cache.cached(w, lemmatizer) for w in words]

        hm.append(cache.hits_miss())
        times.append(sw.elapsed)
        cache.clear()

    return (times, hm)

  4. 初始化一个单词列表和一个 NLTK WordNetLemmatizer对象:

    words = [w for w in brown.words()]
lemmatizer = WordNetLemmatizer()

  5. 测量执行时间如下:

    with dl.perf.StopWatch() as sw:
    _ = [lemmatizer.lemmatize(w.lower()) for w in words]

plain = sw.elapsed

times, hm = measure(clru_cache, words, lemmatizer)

  6. 绘制不同缓存大小的结果:

    sp = dl.plotting.Subplotter(2, 2, context)
sp.ax.plot(2 ** np.arange(5, 12), times)
sp.ax.axhline(plain, lw=2, label='Uncached')
sp.label()

sp.next_ax().plot(2 ** np.arange(5, 12), hm)
sp.label()

times, hm = measure(lru_cache, words, lemmatizer)
sp.next_ax().plot(2 ** np.arange(5, 12), times)
sp.ax.axhline(plain, lw=2, label='Uncached')
sp.label()

sp.next_ax().plot(2 ** np.arange(5, 12), hm)
sp.label()
HTML(sp.exit())

最终结果参见以下截图:

How to do it...

代码在本书代码包的caching_lru.ipynb文件中。

另见

缓存 HTTP 请求

有时，数据可以通过 HTTP 上的网络服务获得。这样做的好处是，我们不必太在意发送方使用的技术。例如，这与电子邮件的工作方式类似。然而，我们必须通过 HTTP GET(通常)或 HTTP POST(按照惯例大写)方法显式地请求信息。每当我们请求一个网页或下载一个文件时，我们通常会执行一个 GET 请求。另一端的 web 服务器必须处理该请求。如果有很多请求，我们可能会降低服务器的速度，因此组织通常会采取措施来防止这种情况。这可能意味着您的进一步请求将被阻止。

出于效率原因，避免多次发出相同的请求也是有利的。网络浏览器用缓存解决了这个问题，我们可以用requests-cache包做同样的事情。默认情况下，缓存存储在 SQLite 数据库中。

我们将不涉及的一个常见用例是用 HTTP 定期检索信息。显然，如果没有任何变化，我们不想检索内容。HTTP 协议提供了确定内容是否被修改的有效机制。但是，web 服务器不需要报告内容更改。

做好准备

使用以下命令安装请求缓存:

$ pip install --upgrade requests-cache

我用请求测试了代码——缓存 0.4.10。

怎么做...

  1. 进口情况如下:

    import requests
import requests_cache

  2. 安装缓存(这将默认创建一个 SQLite 数据库):

    requests_cache.install_cache()

  3. 请求建立缓存的网站:

    %time requests.get('http://google.com')

  4. 请求现在应该来自本地缓存的相同网站:

    %time requests.get('http://google.com')

  5. 清除缓存如下:

    requests_cache.clear()

  6. 再次请求网站(缓存现在应该是空的):

    %time requests.get('http://google.com')

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的 caching_requests.ipynb文件中。

另见

使用最小计数草图进行流式计数

流式或在线算法很有用，因为它们不需要像其他算法那样多的内存和处理能力。本章有一个在线计算统计矩的方法(参考即时计算平均值、方差、偏斜度和峰度)。

此外，在第 5 章、网络挖掘、数据库和大数据的用 Spark 配方聚类流数据中，我介绍了另一种流算法。

由于基本原因或舍入误差，流算法通常是近似的。因此，如果可能的话，你应该尝试使用其他算法。当然，在许多情况下，近似结果已经足够好了。例如，一个用户在社交媒体网站上有 500 或 501 个连接并不重要。如果你只是发送成千上万的邀请，你迟早会到达那里。

素描是你可能从绘画中知道的东西。在这种情况下，素描意味着勾勒出没有任何细节的物体的大致轮廓。类似的概念存在于流算法的世界中。

在这个食谱中，我涵盖了格雷厄姆·科尔莫和 s .穆图·穆图克里希南(S. Muthu Muthukrishnan)的计数-min 草图 (2003 年)，这在排名的上下文中很有用。例如，我们可能想知道新闻网站上浏览次数最多的文章、热门话题、点击次数最多的广告或联系最多的用户。天真的方法需要在字典或表格中记录每个条目的数量。字典使用散列函数来计算作为关键字的识别整数。出于理论上的原因，我们可能会有冲突——这意味着两个或多个项目有相同的密钥。Count-min 草图是一个二维表格数据结构，它很小，并且对每一行都使用散列函数。它容易发生碰撞，导致过度安装。

当事件发生时，例如有人观看广告，我们会执行以下操作:

  1. 对于草图中的每一行，我们应用相关的散列函数，例如，使用广告标识符来获得列索引。
  2. 递增与行和列对应的值。

每个事件都清晰地映射到草图中的每一行。当我们请求计数时，我们遵循相反的路径来获得多个计数。最低计数给出了该项目的估计计数。

这种设置背后的想法是，频繁项目可能会主导不太常见的项目。热门物品与冷门物品碰撞的概率大于热门物品之间的碰撞概率。

怎么做...

  1. 进口情况如下:

    from nltk.corpus import brown
from nltk.corpus import stopwords
import dautil as dl
from collections import defaultdict
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import HTML

  2. 将 NLTK Brown 和 stop 单词语料库的单词存储在列表中:

    words_dict = dl.collect.IdDict()
dd = defaultdict(int)
fid = brown.fileids(categories='news')[0]
words = brown.words(fid)
sw = set(stopwords.words('english'))

  3. 计算每个停止字的出现次数:

    for w in words:
    if w in sw:
        dd[w] += 1

  4. 绘制最小计数草图各种参数的计数误差分布:

    sp = dl.plotting.Subplotter(2, 2, context)
actual = np.array([dd[w] for w in sw])
errors = []

for i in range(1, 4):
    cm = dl.perf.CountMinSketch(depth=5 * 2 ** i,
                                width=20 * 2 ** i)

    for w in words:
        cm.add(words_dict.add_or_get(w.lower()))

    estimates = np.array([cm.estimate_count(words_dict.add_or_get(w))
                        for w in sw])
    diff = estimates - actual
    errors.append(diff)

    if i > 1:
        sp.next_ax()

    sp.ax.hist(diff, normed=True,
            bins=dl.stats.sqrt_bins(actual))
    sp.label()

sp.next_ax().boxplot(errors)
sp.label()
HTML(sp.exit())

有关最终结果，请参考以下屏幕截图:

How to do it...

代码在本书代码包的 stream_demo.py文件中。

另见

利用 OpenCL 利用图形处理器的能力

开放计算语言 ( OpenCL )最初由苹果公司(Apple Inc .)开发，是一个开放的程序技术标准，可以在各种设备上运行，包括在商品硬件上可用的 CPU 和 GPU，比如我用于这个配方的机器。自 2009 年以来，OpenCL 一直由 Khronos 计算工作组维护。许多硬件厂商，包括我比较喜欢的厂商，都有 OpenCL 的实现。

OpenCL 是一种类似 C 的语言(实际上有多种 C 语言方言或版本)，其功能称为 内核。内核可以在多个处理元件上并行运行。硬件供应商给出了处理元素的定义。OpenCL 程序是为了可移植性而在运行时编译的。

可移植性是 OpenCL 相对于类似技术(如作为 NVIDIA 产品的 CUDA)的最大优势。另一个优势是能够在中央处理器、图形处理器和其他设备之间共享工作。有人建议使用机器学习来优化分工。

皮奥尼塔斯可以用 PyOpenCL 包编写 OpenCL 程序。PyOpenCL 向 Python 异常添加了额外的功能，例如对象清理和错误转换。许多其他库使用并在某些方面增强了 pyoppencl(请参考 pyoppencl 文档)。

做好准备

使用以下命令安装pyopencl :

$ pip install pyopencl

我用 PyOpenCL 2015.2.3 测试了代码。更多信息请参考wiki.tiker.net/OpenCLHowTo。

怎么做...

代码在本书代码包的opencl_demo.ipynb文件中:

  1. 进口情况如下:

    import pyopencl as cl
from pyopencl import array
import numpy as np

  2. 定义以下函数来接受 NumPy 数组并执行简单的计算:

    def np_se(a, b):
    return (a - b) ** 2

  3. 定义以下函数，使用 OpenCL 进行与上一步相同的计算:

    def gpu_se(a, b, platform, device, context, program):

  4. 创建一个队列，启用概要分析(仅用于演示)和缓冲区，以便在

        queue = cl.CommandQueue(context,
                            properties=cl.command_queue_properties.
                            PROFILING_ENABLE)
    mem_flags = cl.mem_flags
    a_buf = cl.Buffer(context,
                      mem_flags.READ_ONLY | mem_flags.COPY_HOST_PTR,
                      hostbuf=a)
    b_buf = cl.Buffer(context,
                      mem_flags.READ_ONLY | mem_flags.COPY_HOST_PTR, hostbuf=b)
    error = np.empty_like(a)
    destination_buf = cl.Buffer(context,
                                mem_flags.WRITE_ONLY,
                                error.nbytes)

    周围洗牌

  5. 执行 OpenCL 程序并分析代码:

        exec_evt = program.mean_squared_error(queue, error.shape, None,
                                          a_buf, b_buf, destination_buf)
    exec_evt.wait()
    elapsed = 1e-9*(exec_evt.profile.end - exec_evt.profile.start)

    print("Execution time of OpenCL: %g s" % elapsed)

    cl.enqueue_copy(queue,
                    error, destination_buf)

    return error

  6. 生成随机数据如下:

    np.random.seed(51)
a = np.random.rand(4096).astype(np.float32)
b = np.random.rand(4096).astype(np.float32)

  7. 访问中央处理器和图形处理器。这部分依赖于硬件，所以您可能需要更改这些行:

    platform = cl.get_platforms()[0]
device = platform.get_devices()[2]
context = cl.Context([device])

  8. 用 OpenCL 语言定义一个内核:

    program = cl.Program(context, """
    __kernel void mean_squared_error(__global const float *a,
    __global const float *b, __global float *result)
    {
        int gid = get_global_id(0);
        float temp = a[gid] - b[gid];
        result[gid] =  temp * temp;
    }
        """).build()

  9. 用 NumPy 和 OpenCL (GPU)计算平方误差，并测量执行时间:

    gpu_error = gpu_se(a, b, platform, device, context, program)

np_error = np_se(a, b)
print('GPU error', np.mean(gpu_error))
print('NumPy error', np.mean(np_error))
%time np_se(a, b)

有关最终结果，请参考以下屏幕截图:

How to do it...

另见

avatar
文章
阅读
粉丝