NumPy-初学者指南中文第三版-四-

NumPy 初学者指南中文第三版（四）

原文：NumPy: Beginner’s Guide - Third Edition

协议：CC BY-NC-SA 4.0

十、当 NumPy 不够用时 - SciPy 及更多

SciPy 是建立在 NumPy 之上的世界著名的 Python 开源科学计算库。 它增加了一些功能，例如数值积分，优化，统计和特殊函数。

在本章中，我们将介绍以下主题：

• 文件 I/O
• 统计
• 信号处理
• 优化
• 插值
• 图像和音频处理

MATLAB 和 Octave

MATLAB 及其开源替代品 Octave 是流行的数学程序。 scipy.io包具有一些函数，可让您加载 MATLAB 或 Octave 矩阵，以及数字或 Python 程序中的字符串，反之亦然。 loadmat()函数加载.mat文件。 savemat()函数将名称和数组的字典保存到.mat文件中。

实战时间 – 保存并加载.mat文件

如果我们从 NumPy 数组开始并决定在 MATLAB 或 Octave 环境中使用所述数组，那么最简单的方法就是创建一个.mat文件。 然后，我们可以在 MATLAB 或 Octave 中加载文件。 让我们完成必要的步骤：

1. 创建一个 NumPy 数组，然后调用 savemat()函数来创建.mat文件。 该函数有两个参数：文件名和包含变量名和值的字典：

   a = np.arange(7)
   
   io.savemat("a.mat", {"array": a})
   

2. 在 MATLAB 或 Octave 环境中，加载.mat文件并检查存储的数组：

   octave-3.4.0:7> load a.mat
   octave-3.4.0:8> a
   
   octave-3.4.0:8> array
   array =
   
     0
     1
     2
     3
     4
     5
     6
   

刚刚发生了什么？

我们从 NumPy 代码创建了一个.mat文件，并将其加载到 Octave 中。 我们检查了创建的 NumPy 数组（请参见scipyio.py）：

import numpy as np
from scipy import io

a = np.arange(7)

io.savemat("a.mat", {"array": a})

小测验 - 加载.mat文件

Q1. 哪个函数加载.mat文件？

1. Loadmatlab
2. loadmat
3. loadoct
4. frommat

统计

SciPy 统计模块为 ，称为scipy.stats。 一类实现连续分布 ，一类实现离散分布。 同样，在此模块中，可以找到执行大量统计检验的函数。

实战时间 – 分析随机值

我们将生成模拟正态分布的随机值，并使用scipy.stats包中的统计函数分析生成的数据。

1. 使用scipy.stats包从正态分布生成随机值：

   generated = stats.norm.rvs(size=900)
   

2. 将生成的值拟合为正态分布。 这基本上给出了数据集的平均值和标准偏差：

   print("Mean", "Std", stats.norm.fit(generated))
   

   平均值和标准差如下所示：

   Mean Std (0.0071293257063200707, 0.95537708218972528)
   
   

3. 偏度告诉我们概率分布有多偏斜（不对称）。执行偏度检验。 该检验返回两个值。 第二个值是 p 值 -- 数据集的偏斜度不符合正态分布的概率。

   注意

   一般而言，p 值是结果与给定零假设所期望的结果不同的概率，在这种情况下，偏度与正态分布（由于对称而为 0）不同的概率。

   P 值的范围是0至1：

   print("Skewtest", "pvalue", stats.skewtest(generated))
   

   偏度检验的结果如下所示：

   Skewtest pvalue (-0.62120640688766893, 0.5344638245033837)
   
   

   因此，我们不处理正态分布的可能性为53% 。 观察如果我们生成更多的点会发生什么，这很有启发性，因为如果我们生成更多的点，我们应该具有更正态的分布。 对于 900,000 点，我们得到0.16的 p 值。 对于 20 个生成的值，p 值为0.50。

4. 峰度告诉我们概率分布的弯曲程度。 执行峰度检验。 此检验的设置与偏度检验类似，但当然适用于峰度：

   print("Kurtosistest", "pvalue", stats.kurtosistest(generated))
   

   峰度检验的结果显示如下：

   Kurtosistest pvalue (1.3065381019536981, 0.19136963054975586)
   
   

   900,000 个值的 p 值为0.028。 对于 20 个生成的值，p 值为0.88。

5. 正态检验告诉我们数据集符合正态分布的可能性。 执行正态性检验。 此检验还返回两个值，其中第二个是p值：

   print("Normaltest", "pvalue", stats.normaltest(generated))
   

   正态性检验的结果如下所示：

   Normaltest pvalue (2.09293921181506, 0.35117535059841687)
   
   

   900,000 个生成值的 p 值为0.035。 对于 20 个生成的值，p 值为0.79。

6. 我们可以使用 SciPy 轻松找到某个百分比的值：

   print("95 percentile", stats.scoreatpercentile(generated, 95))
   

   95th百分位的值显示如下：

   95 percentile 1.54048860252
   
   

7. 进行与上一步相反的操作，以找到 1 处的百分位数：

   print("Percentile at 1", stats.percentileofscore(generated, 1))
   

   1处的百分位数显示如下：

   Percentile at 1 85.5555555556
   
   

8. 使用matplotlib在直方图中绘制生成的值（有关matplotlib的更多信息，请参见前面的第 9 章， “Matplotlib 绘图”：

   plt.hist(generated)
   

   生成的随机值的直方图如下：

   

刚刚发生了什么？

我们从正态分布创建了一个数据集，并使用scipy.stats模块对其进行了分析（请参见statistics.py）：

from __future__ import print_function
from scipy import stats
import matplotlib.pyplot as plt

generated = stats.norm.rvs(size=900)
print("Mean", "Std", stats.norm.fit(generated))
print("Skewtest", "pvalue", stats.skewtest(generated))
print("Kurtosistest", "pvalue", stats.kurtosistest(generated))
print("Normaltest", "pvalue", stats.normaltest(generated))
print("95 percentile", stats.scoreatpercentile(generated, 95))
print("Percentile at 1", stats.percentileofscore(generated, 1))
plt.title('Histogram of 900 random normally distributed values')
plt.hist(generated)
plt.grid()
plt.show()

勇往直前 - 改善数据生成

从前面的“实战时间”部分中的直方图来看，在生成数据方面还有改进的余地。 尝试使用 NumPy 或scipy.stats.norm.rvs()函数的其他参数。

SciKits 样本比较

通常，我们有两个数据样本，可能来自不同的实验，它们之间存在某种关联。 存在可以比较样本的统计检验。 其中一些是在scipy.stats模块中实现的。

我喜欢的另一个统计检验是scikits.statsmodels.stattools的 Jarque-Bera 正态性检验。 SciKit 是小型实验 Python 软件工具箱。 它们不属于 SciPy。 还有 Pandas，这是scikits.statsmodels的分支。 可以在这个页面上找到 SciKit 的列表。 您可以使用安装工具通过以下工具安装statsmodels：

$ [sudo] easy_install statsmodels

实战时间 – 比较股票对数收益

我们将使用matplotlib下载两个追踪器的去年股票报价。 如先前的第 9 章，“matplotlib 绘图”，我们可以从 Yahoo Finance 检索报价。 我们将比较DIA和SPY的收盘价的对数回报（DIA 跟踪道琼斯指数； SPY 跟踪 S&P 500 指数）。 我们还将对返回值的差异执行 Jarque–Bera 检验。

1. 编写一个可以返回指定股票的收盘价的函数：

   def get_close(symbol):
      today = date.today()
      start = (today.year - 1, today.month, today.day)
   
      quotes = quotes_historical_yahoo(symbol, start, today)
      quotes = np.array(quotes)
   
      return quotes.T[4]
   

2. 计算 DIA 和 SPY 的对数返回。 通过采用收盘价的自然对数，然后采用连续值的差来计算对数收益：

   spy =  np.diff(np.log(get_close("SPY")))
   dia =  np.diff(np.log(get_close("DIA")))
   

3. 均值比较测试检查两个不同的样本是否可以具有相同的平均值。 返回两个值，第二个是从 0 到 1 的 p 值：

   print("Means comparison", stats.ttest_ind(spy, dia))
   

   均值比较检验的结果如下所示：

   Means comparison (-0.017995865641886155, 0.98564930169871368)
   
   

   因此，这两个样本有大约 98% 的机会具有相同的平均对数回报。 实际上，该文档的内容如下：

   注意

   如果我们观察到较大的 p 值（例如，大于 0.05 或 0.1），那么我们将无法拒绝具有相同平均分数的原假设。 如果 p 值小于阈值，例如 1% ，5% 或 10% ，则我们拒绝均值的零假设。

4. Kolmogorov–Smirnov 双样本检验告诉我们从同一分布中抽取两个样本的可能性：

   print("Kolmogorov smirnov test", stats.ks_2samp(spy, dia))
   

   再次返回两个值，其中第二个值为 p 值：

   Kolmogorov smirnov test (0.063492063492063516, 0.67615647616238039)
   
   

5. 对对数返回值的差异进行 Jarque–Bera 正态性检验：

   print("Jarque Bera test", jarque_bera(spy – dia)[1])
   

   Jarque-Bera 正态性检验的 p 值显示如下：

   Jarque Bera test 0.596125711042
   
   

6. 用matplotlib绘制对数收益的直方图及其差值：

   plt.hist(spy, histtype="step", lw=1, label="SPY")
   plt.hist(dia, histtype="step", lw=2, label="DIA")
   plt.hist(spy - dia, histtype="step", lw=3, label="Delta")
   plt.legend()
   plt.show()
   

   对数收益和差异的直方图如下所示：

   

刚刚发生了什么？

我们比较了 DIA 和 SPY 的对数回报样本。 另外，我们对对数返回值的差进行了 Jarque-Bera 检验（请参见pair.py）：

from __future__ import print_function
from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
from scipy import stats
from statsmodels.stats.stattools import jarque_bera
import matplotlib.pyplot as plt

def get_close(symbol):
   today = date.today()
   start = (today.year - 1, today.month, today.day)
   quotes = quotes_historical_yahoo(symbol, start, today)
   quotes = np.array(quotes)
   return quotes.T[4]

spy =  np.diff(np.log(get_close("SPY")))
dia =  np.diff(np.log(get_close("DIA")))

print("Means comparison", stats.ttest_ind(spy, dia))
print("Kolmogorov smirnov test", stats.ks_2samp(spy, dia))

print("Jarque Bera test", jarque_bera(spy - dia)[1])

plt.title('Log returns of SPY and DIA')
plt.hist(spy, histtype="step", lw=1, label="SPY")
plt.hist(dia, histtype="step", lw=2, label="DIA")
plt.hist(spy - dia, histtype="step", lw=3, label="Delta")
plt.xlabel('Log returns')
plt.ylabel('Counts')
plt.grid()
plt.legend(loc='best')
plt.show()

信号处理

scipy.signal模块包含过滤函数和 B 样条插值算法。

注意

样条插值使用称为样条的多项式进行插值 ）。 然后，插值尝试将样条线粘合在一起以拟合数据。 B 样条是样条的一种。

SciPy 信号定义为数字数组。 过滤器的一个示例是detrend()函数。 此函数接收信号并对其进行线性拟合。 然后从原始输入数据中减去该趋势。

实战时间 – 检测QQQ趋势

通常对数据样本的趋势比对其去趋势更感兴趣。 在下降趋势之后，我们仍然可以轻松地恢复趋势。 让我们以QQQ的年价格数据为例。

1. 编写获取QQQ收盘价和相应日期的代码：

   today = date.today()
   start = (today.year - 1, today.month, today.day)
   
   quotes = quotes_historical_yahoo("QQQ", start, today)
   quotes = np.array(quotes)
   
   dates = quotes.T[0]
   qqq = quotes.T[4]
   

2. 消除趋势：

   y = signal.detrend(qqq)
   

3. 为日期创建月和日定位器：

   alldays = DayLocator()
   months = MonthLocator()
   

4. 创建一个日期格式化器，该日期格式化器创建月份名称和年份的字符串：

   month_formatter = DateFormatter("%b %Y")
   

5. 创建图形和子图：

   fig = plt.figure()
   ax = fig.add_subplot(111)
   

6. 通过减去去趋势信号绘制数据和潜在趋势：

   plt.plot(dates, qqq, 'o', dates, qqq - y, '-')
   

7. 设置定位器和格式化器：

   ax.xaxis.set_minor_locator(alldays)
   ax.xaxis.set_major_locator(months)
   ax.xaxis.set_major_formatter(month_formatter)
   

8. 将 x 轴标签的格式设置为日期：

   fig.autofmt_xdate()
   plt.show()
   

   下图显示了带有趋势线的 QQQ 价格：

   

刚刚发生了什么？

我们用趋势线绘制了 QQQ 的收盘价（请参见trend.py）：

from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator

today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo("QQQ", start, today)
quotes = np.array(quotes)

dates = quotes.T[0]
qqq = quotes.T[4]

y = signal.detrend(qqq)

alldays = DayLocator()
months = MonthLocator()
month_formatter = DateFormatter("%b %Y")

fig = plt.figure()
ax = fig.add_subplot(111)

plt.title('QQQ close price with trend')
plt.ylabel('Close price')
plt.plot(dates, qqq, 'o', dates, qqq - y, '-')
ax.xaxis.set_minor_locator(alldays)
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(month_formatter)
fig.autofmt_xdate()
plt.grid()
plt.show()

傅立叶分析

现实世界中的信号通常具有周期性。 处理这些信号的常用工具是离散傅里叶变换。 离散傅立叶变换是从时域到频域的变换，即将周期信号线性分解为各种频率的正弦和余弦函数：

可以在scipy.fftpack模块中找到傅里叶变换的函数（NumPy 也有自己的傅里叶包numpy.fft）。 该包中包括快速傅立叶变换，微分和伪微分运算符，以及一些辅助函数。 MATLAB 用户将很高兴地知道scipy.fftpack模块中的许多函数与 MATLAB 的对应函数具有相同的名称，并且与 MATLAB 的等效函数具有相似的功能。

实战时间 – 过滤去趋势的信号

在前面的“实战时间”部分中，我们学习了如何使信号逆趋势。 该去趋势的信号可以具有循环分量。 让我们尝试将其可视化。 其中一些步骤是前面“实战时间”部分中的步骤的重复，例如下载数据和设置matplotlib对象。 这些步骤在此省略。

1. 应用傅立叶变换，得到频谱：

   amps = np.abs(fftpack.fftshift(fftpack.rfft(y)))
   

2. 滤除噪音。 假设，如果频率分量的幅度低于最强分量的10% ，则将其丢弃：

   amps[amps < 0.1 * amps.max()] = 0
   

3. 将过滤后的信号转换回原始域，并将其与去趋势的信号一起绘制：

   plt.plot(dates, y, 'o', label="detrended")
   plt.plot(dates, -fftpack.irfft(fftpack.ifftshift(amps)), label="filtered")
   

4. 将 x 轴标签格式化为日期，并添加具有超大尺寸的图例：

   fig.autofmt_xdate()
   plt.legend(prop={'size':'x-large'})
   

5. 添加第二个子图并在过滤后绘制频谱图：

   ax2 = fig.add_subplot(212)
   N = len(qqq)
   plt.plot(np.linspace(-N/2, N/2, N), amps, label="transformed")
   

6. 显示图例和图解：

   plt.legend(prop={'size':'x-large'})
   
   plt.show()
   

   下图是信号和频谱的图：

   

刚刚发生了什么？

我们对信号进行了去趋势处理，然后使用scipy.fftpack模块在其上应用了一个简单的过滤器（请参阅frequencies.py）：

from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from scipy import fftpack
from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator

today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo("QQQ", start, today)
quotes = np.array(quotes)

dates = quotes.T[0]
qqq = quotes.T[4]

y = signal.detrend(qqq)

alldays = DayLocator()
months = MonthLocator()
month_formatter = DateFormatter("%b %Y")

fig = plt.figure()
fig.subplots_adjust(hspace=.3)
ax = fig.add_subplot(211)

ax.xaxis.set_minor_locator(alldays)
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(month_formatter)

## make font size bigger
ax.tick_params(axis='both', which='major', labelsize='x-large')

amps = np.abs(fftpack.fftshift(fftpack.rfft(y)))
amps[amps < 0.1 * amps.max()] = 0

plt.title('Detrended and filtered signal')
plt.plot(dates, y, 'o', label="detrended")
plt.plot(dates, -fftpack.irfft(fftpack.ifftshift(amps)), label="filtered")
fig.autofmt_xdate()
plt.legend(prop={'size':'x-large'})
plt.grid()

ax2 = fig.add_subplot(212)
plt.title('Transformed signal')
ax2.tick_params(axis='both', which='major', labelsize='x-large')
N = len(qqq)
plt.plot(np.linspace(-N/2, N/2, N), amps, label="transformed")

plt.legend(prop={'size':'x-large'})
plt.grid()
plt.tight_layout()
plt.show()

数学优化

优化算法试图找到问题的最佳解决方案，例如，找到函数的最大值或最小值。 该函数可以是线性的或非线性的。 该解决方案也可能具有特殊的约束。 例如，可能不允许解决方案具有负值。 scipy.optimize模块提供了几种优化算法。 算法之一是最小二乘拟合函数leastsq()。 调用此函数时，我们提供了残差（错误项）函数。 此函数可将残差平方和最小化。 它对应于我们的解决方案数学模型。 还必须给算法一个起点。 这应该是一个最佳猜测-尽可能接近真实的解决方案。 否则，将在大约100 * (N+1)次迭代后停止执行，其中 N 是要优化的参数数量。

实战时间 – 正弦拟合

在前面的“实战时间”部分中，我们为脱趋势数据创建了一个简单的过滤器。 现在，让我们使用限制性更强的过滤器，该过滤器将只剩下主要频率分量。 我们将为其拟合正弦波模式并绘制结果。 该模型具有四个参数-幅度，频率，相位和垂直偏移。

1. 根据正弦波模型定义残差函数：

   def residuals(p, y, x):
      A,k,theta,b = p
      err = y-A * np.sin(2* np.pi* k * x + theta) + b
      return err
   

2. 将过滤后的信号转换回原始域：

   filtered = -fftpack.irfft(fftpack.ifftshift(amps))
   

3. 猜猜我们试图估计的从时域到频域的转换的参数值：

   N = len(qqq)
   f = np.linspace(-N/2, N/2, N)
   p0 = [filtered.max(), f[amps.argmax()]/(2*N), 0, 0]
   print("P0", p0)
   

   初始值如下所示：

   P0 [2.6679532410065212, 0.00099598469163686377, 0, 0]
   
   

4. 调用leastsq()函数：

   plsq = optimize.leastsq(residuals, p0, args=(filtered, dates))
   p = plsq[0]
   print("P", p)
   

   最终参数值如下：

   P [  2.67678014e+00   2.73033206e-03  -8.00007036e+03  -5.01260321e-03]
   
   

5. 用去趋势数据，过滤后的数据和过滤后的数据拟合完成第一个子图。 将日期格式用于水平轴并添加图例：

   plt.plot(dates, y, 'o', label="detrended")
   plt.plot(dates, filtered, label="filtered")
   plt.plot(dates, p[0] * np.sin(2 * np.pi * dates * p[1] + p[2]) + p[3], '^', label="fit")
   fig.autofmt_xdate()
   plt.legend(prop={'size':'x-large'})
   

6. 添加第二个子图，其中包含频谱主要成分的图例：

   ax2 = fig.add_subplot(212)
   plt.plot(f, amps, label="transformed")
   

   以下是结果图表：

   

刚刚发生了什么？

我们降低了 QQQ 一年价格数据的趋势。 然后对该信号进行过滤，直到仅剩下频谱的主要成分。 我们使用scipy.optimize模块（请参见optfit.py）将正弦拟合到过滤后的信号：

from __future__ import print_function
from matplotlib.finance import quotes_historical_yahoo
import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack
from scipy import signal
from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator
from scipy import optimize

start = (2010, 7, 25)
end = (2011, 7, 25)

quotes = quotes_historical_yahoo("QQQ", start, end)
quotes = np.array(quotes)

dates = quotes.T[0]
qqq = quotes.T[4]

y = signal.detrend(qqq)

alldays = DayLocator()
months = MonthLocator()
month_formatter = DateFormatter("%b %Y")

fig = plt.figure()
fig.subplots_adjust(hspace=.3)
ax = fig.add_subplot(211)

ax.xaxis.set_minor_locator(alldays)
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(month_formatter)
ax.tick_params(axis='both', which='major', labelsize='x-large')

amps = np.abs(fftpack.fftshift(fftpack.rfft(y)))
amps[amps < amps.max()] = 0

def residuals(p, y, x):
   A,k,theta,b = p
   err = y-A * np.sin(2* np.pi* k * x + theta) + b
   return err

filtered = -fftpack.irfft(fftpack.ifftshift(amps))
N = len(qqq)
f = np.linspace(-N/2, N/2, N)
p0 = [filtered.max(), f[amps.argmax()]/(2*N), 0, 0]
print("P0", p0)

plsq = optimize.leastsq(residuals, p0, args=(filtered, dates))
p = plsq[0]
print("P", p)
plt.title('Detrended and filtered signal')
plt.plot(dates, y, 'o', label="detrended")
plt.plot(dates, filtered, label="filtered")
plt.plot(dates, p[0] * np.sin(2 * np.pi * dates * p[1] + p[2]) + p[3], '^', label="fit")
fig.autofmt_xdate()
plt.legend(prop={'size':'x-large'})
plt.grid()

ax2 = fig.add_subplot(212)
plt.title('Tranformed signal')
ax2.tick_params(axis='both', which='major', labelsize='x-large')
plt.plot(f, amps, label="transformed")

plt.legend(prop={'size':'x-large'})
plt.grid()
plt.tight_layout()
plt.show()

数值积分

SciPy 具有数值积分包scipy.integrate，在 NumPy 中没有等效项。 quad()函数可以在两个点之间整合一个单变量函数。 这些点可以是无穷大。 该函数使用最简单的数值积分方法：梯形法则。

实战时间 – 计算高斯积分

高斯积分是error()函数相关（在数学上也称为erf），但没有限制。 计算结果为pi的平方根。

让我们用quad()函数计算积分（对于导入，请检查代码包中的文件）：

print("Gaussian integral", np.sqrt(np.pi),integrate.quad(lambda x: np.exp(-x**2), -np.inf, np.inf))

返回值是结果，其错误如下：

Gaussian integral 1.77245385091 (1.7724538509055159, 1.4202636780944923e-08)

刚刚发生了什么？

我们使用quad()函数计算了高斯积分。

勇往直前 – 更多实验

试用同一包中的其他集成函数。 只需替换一个函数调用即可。 我们应该得到相同的结果，因此您可能还需要阅读文档以了解更多信息。

插值

插值填充数据集中已知数据点之间的空白。 scipy.interpolate()函数根据实验数据对函数进行插值。 interp1d类可以创建线性或三次插值函数。 默认情况下，会创建线性插值函数，但是如果设置了kind参数，则会创建三次插值函数。 interp2d类的工作方式相同，但是是 2D 的。

实战时间 – 一维内插

我们将使用 sinc()函数创建数据点，并向其中添加一些随机噪声。 之后，我们将进行线性和三次插值并绘制结果。

1. 创建数据点并为其添加噪声：

   x = np.linspace(-18, 18, 36)
   noise = 0.1 * np.random.random(len(x))
   signal = np.sinc(x) + noise
   

2. 创建一个线性插值函数，并将其应用于具有五倍数据点的输入数组：

   interpreted = interpolate.interp1d(x, signal)
   x2 = np.linspace(-18, 18, 180)
   y = interpreted(x2)
   

3. 执行与上一步相同的操作，但使用三次插值：

   cubic = interpolate.interp1d(x, signal, kind="cubic")
   y2 = cubic(x2)
   

4. 用matplotlib绘制结果：

   plt.plot(x, signal, 'o', label="data")
   plt.plot(x2, y, '-', label="linear")
   plt.plot(x2, y2, '-', lw=2, label="cubic")
   plt.legend()
   plt.show()
   

   下图是数据，线性和三次插值的图形：

   

刚刚发生了什么？

我们通过sinc()函数创建了一个数据集，并添加了噪声。 然后，我们使用scipy.interpolate模块的interp1d类（请参见sincinterp.py）进行了线性和三次插值 ：

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt

x = np.linspace(-18, 18, 36)
noise = 0.1 * np.random.random(len(x))
signal = np.sinc(x) + noise

interpreted = interpolate.interp1d(x, signal)
x2 = np.linspace(-18, 18, 180)
y = interpreted(x2)

cubic = interpolate.interp1d(x, signal, kind="cubic")
y2 = cubic(x2)

plt.plot(x, signal, 'o', label="data")
plt.plot(x2, y, '-', label="linear")
plt.plot(x2, y2, '-', lw=2, label="cubic")

plt.title('Interpolated signal')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend(loc='best')
plt.show()

图像处理

使用 SciPy，我们可以使用scipy.ndimage包进行图像处理。 该模块包含各种图像过滤器和工具。

实战时间 - 操纵 Lena

scipy.misc模块是一个加载“Lena”图像的工具。 这是 Lena Soderberg 的图像，传统上用于图像处理示例。 我们将对该图像应用一些过滤器并旋转它。 执行以下步骤以执行 ：

1. 加载 Lena 图像并将其显示在带有灰度色图的子图中：

   image = misc.lena().astype(np.float32)
   plt.subplot(221)
   plt.title("Original Image")
   img = plt.imshow(image, cmap=plt.cm.gray)
   

   请注意，我们正在处理float32数组。

2. 中值过滤器扫描图像，并用相邻数据点的中值替换每个项目。 对图像应用中值过滤器，然后在第二个子图中显示它：

   plt.subplot(222)
   plt.title("Median Filter")
   filtered = ndimage.median_filter(image, size=(42,42))
   plt.imshow(filtered, cmap=plt.cm.gray)
   

3. 旋转图像并将其显示在第三个子图中：

   plt.subplot(223)
   plt.title("Rotated")
   rotated = ndimage.rotate(image, 90)
   plt.imshow(rotated, cmap=plt.cm.gray)
   

4. Prewitt 过滤器基于计算图像强度的梯度。 将 Prewitt 过滤器应用于图像，并在第四个子图中显示它：

   plt.subplot(224)
   plt.title("Prewitt Filter")
   filtered = ndimage.prewitt(image)
   plt.imshow(filtered, cmap=plt.cm.gray)
   plt.show()
   

   以下是生成的图像：

   

刚刚发生了什么？

我们使用 Lena 的图像：

from scipy import misc
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage

image = misc.lena().astype(np.float32)

plt.subplot(221)
plt.title("Original Image")
img = plt.imshow(image, cmap=plt.cm.gray)
plt.axis("off")

plt.subplot(222)
plt.title("Median Filter")
filtered = ndimage.median_filter(image, size=(42,42))
plt.imshow(filtered, cmap=plt.cm.gray)
plt.axis("off")

plt.subplot(223)
plt.title("Rotated")
rotated = ndimage.rotate(image, 90)
plt.imshow(rotated, cmap=plt.cm.gray)
plt.axis("off")

plt.subplot(224)
plt.title("Prewitt Filter")
filtered = ndimage.prewitt(image)
plt.imshow(filtered, cmap=plt.cm.gray)
plt.axis("off")
plt.show()

音频处理

既然我们已经完成了一些图像处理，那么您也可以使用 WAV 文件来完成令人兴奋的事情，您可能不会感到惊讶。 让我们下载一个 WAV 文件并重播几次。 我们将跳过下载部分的解释，该部分只是常规的 Python。

实战时间 – 重放音频片段

我们将下载 Austin Powers 的 WAV 文件，称为“Smashing baby”。 可以使用scipy.io.wavfile模块中的 read()函数将此文件转换为 NumPy 数组。 相同包中的 write()函数将在本节末尾用于创建新的 WAV 文件。 我们将进一步使用 tile()函数重播音频剪辑几次。

1. 使用read()函数读取文件：

   sample_rate, data = wavfile.read(WAV_FILE)
   

   这给了我们两项–采样率和音频数据。 对于本节，我们仅对音频数据感兴趣。

2. 应用tile()函数：

   repeated = np.tile(data, 4)
   

3. 使用write()函数编写一个新文件：

   wavfile.write("repeated_yababy.wav", sample_rate, repeated)
   

   下图显示了四次重复的原始音频数据和音频剪辑：

   

刚刚发生了什么？

我们读取一个音频剪辑，将其重复四次，然后使用新数组（请参见repeat_audio.py）创建一个新的 WAV 文件：

from __future__ import print_function
from scipy.io import wavfile
import matplotlib.pyplot as plt
import urllib.request
import numpy as np

response = urllib.request.urlopen('http://www.thesoundarchive.com/austinpowers/smashingbaby.wav')
print(response.info())
WAV_FILE = 'smashingbaby.wav'
filehandle = open(WAV_FILE, 'wb')
filehandle.write(response.read())
filehandle.close()
sample_rate, data = wavfile.read(WAV_FILE)
print("Data type", data.dtype, "Shape", data.shape)

plt.subplot(2, 1, 1)
plt.title("Original audio signal")
plt.plot(data)
plt.grid()

plt.subplot(2, 1, 2)

## Repeat the audio fragment
repeated = np.tile(data, 4)

## Plot the audio data
plt.title("Repeated 4 times")
plt.plot(repeated)
wavfile.write("repeated_yababy.wav",
    sample_rate, repeated)
plt.grid()
plt.tight_layout()
plt.show()

总结

在本章中，我们仅介绍了 SciPy 和 SciKits 可以实现的功能。 尽管如此，我们还是学到了一些有关文件 I/O，统计量，信号处理，优化，插值，音频和图像处理的知识。

在下一章中，我们将使用 Pygame（开源 Python 游戏库）创建一些简单而有趣的游戏。 在此过程中，我们将学习 NumPy 与 Pygame，Scikit 机器学习库，以及其他的集成。

十一、玩转 Pygame

本章适用于希望使用 NumPy 和 Pygame 快速轻松创建游戏的开发人员。 基本的游戏开发经验会有所帮助，但这不是必需的。

您将学到的东西如下：

• pygame 基础
• matplotlib 集成
• 表面像素数组
• 人工智能
• 动画
• OpenGL

Pygame

Pygame 是 Python 框架，最初由 Pete Shinners 编写， 顾名思义，可用于制作视频游戏。 自 2004 年以来，Pygame 是免费的开放源代码，并获得 GPL 许可，这意味着您基本上可以制作任何类型的游戏。 Pygame 构建在简单 DirectMedia 层（SDL）。 SDL 是一个 C 框架，可以访问各种操作系统（包括 Linux，MacOSX 和 Windows）上的图形，声音，键盘和其他输入设备。

实战时间 – 安装 Pygame

我们将在本节中安装 Pygame。 Pygame 应该与所有 Python 版本兼容。 在撰写时，Python3 存在一些不兼容问题，但很可能很快就会解决。

• 在 Debian 和 Ubuntu 上安装：Pygame 可以在 Debian 档案文件中找到。

• 在 Windows 上安装：从 Pygame 网站下载适用于您正在使用的版本的 Python 的二进制安装程序。

• 在 Mac 上安装 Pygame：适用于 MacOSX 10.3 及更高版本的二进制 Pygame 包可在这个页面中找到。

• 从源代码安装：Pygame 使用distutils系统进行编译和安装。 要开始使用默认选项安装 Pygame，只需运行以下命令：

  $ python setup.py
  
  

  如果您需要有关可用选项的更多信息，请键入以下内容：

  $ python setup.py help
  
  

• 要编译代码，您的操作系统需要一个编译器。 进行设置超出了本书的范围。 有关在 Windows 上编译 Pygame 的更多信息，可以在这个页面上找到。 有关在 MacOSX 上编译 Pygame 的更多信息，请参考这里。

HelloWorld

我们将创建一个简单的游戏，在本章中我们将进一步改进 。 与编程书籍中的传统方法一样，我们从Hello World!示例开始。

实战时间 – 创建一个简单的游戏

重要的是要注意所谓的主游戏循环，在该循环中所有动作都会发生，并使用Font模块渲染文本。 在此程序中，我们将操纵用于绘制的 Pygame Surface对象，并处理退出事件。

1. 首先，导入所需的 Pygame 模块。 如果正确安装了 Pygame，则不会出现任何错误，否则请返回安装“实战时间”：

   import pygame, sys
   from pygame.locals import *
   

2. 初始化 Pygame，按300像素创建400的显示，并将窗口标题设置为Hello world!：

   pygame.init()
   screen = pygame.display.set_mode((400, 300))
   
   pygame.display.set_caption('Hello World!')
   

3. 游戏通常会有一个游戏循环，该循环将一直运行直到发生退出事件为止。 在此示例中，仅在坐标(100, 100)上设置带有文本Hello world!的标签。 文字的字体大小为 19，颜色为红色：

   while True: 
      sysFont = pygame.font.SysFont("None", 19)
      rendered = sysFont.render('Hello World', 0, (255, 100, 100))
      screen.blit(rendered, (100, 100))
   
      for event in pygame.event.get():
         if event.type == QUIT:
            pygame.quit()
            sys.exit()
   
      pygame.display.update()
   

   我们得到以下屏幕截图作为最终结果：

   

   以下是 HelloWorld 的完整代码示例：

   import pygame, sys
   from pygame.locals import *
   
   pygame.init()
   screen = pygame.display.set_mode((400, 300))
   
   pygame.display.set_caption('Hello World!')
   
   while True: 
      sysFont = pygame.font.SysFont("None", 19)
      rendered = sysFont.render('Hello World', 0, (255, 100, 100))
      screen.blit(rendered, (100, 100))
   
      for event in pygame.event.get():
         if event.type == QUIT:
            pygame.quit()
            sys.exit()
   
      pygame.display.update()
   

刚刚发生了什么？

看起来似乎不多，但是我们在本节中学到了很多东西。 下表总结了通过审查的函数：

函数	描述
pygame.init()	此函数执行初始化，您必须在调用其他 Pygame 函数之前调用它。
pygame.display.set_mode((400, 300))	此函数创建一个要使用的所谓的Surface对象。 我们给这个函数一个表示表面尺寸的元组。
pygame.display.set_caption('Hello World!')	此函数将窗口标题设置为指定的字符串值。
pygame.font.SysFont("None", 19)	此函数根据逗号分隔的字体列表（在本例中为无）和整数字体大小参数创建系统字体。
sysFont.render('Hello World', 0, (255, 100, 100))	此函数在Surface上绘制文本。 最后一个参数是表示颜色的 RGB 值的元组。
screen.blit(rendered, (100, 100))	此函数使用Surface。
pygame.event.get()	此函数获取Event对象的列表。 事件表示系统中的特殊事件，例如用户退出游戏。
pygame.quit()	该函数清除由 Pygame 使用的资源。 退出游戏之前，请调用此函数。
pygame.display.update()	此函数刷新表面。

动画

大多数游戏，甚至是最静态的游戏，都有一定程度的动画效果。 从程序员的角度来看，动画就是 ，无非就是在不同的时间在不同的位置显示对象，从而模拟运动。

Pygame 提供了一个Clock对象，该对象管理每秒绘制多少帧。 这样可以确保动画与用户 CPU 的速度无关。

实战时间 – 使用 NumPy 和 Pygame 为对象设置动画

我们将加载图像，然后再次使用 NumPy 定义屏幕周围的顺时针路径。

1. 创建一个 Pygame 时钟，如下所示：

   clock = pygame.time.Clock()
   

2. 作为本书随附的源代码的一部分，应该有一张头像。 加载此图像并在屏幕上四处移动：

   img = pygame.image.load('head.jpg')
   

3. 定义一些数组来保存位置的坐标，我们希望在动画过程中将图像放置在这些位置。 由于我们将移动对象，因此路径有四个逻辑部分：right，down，left和up。 每个部分将具有40等距步长。 将0部分中的所有值初始化：

   steps = np.linspace(20, 360, 40).astype(int)
   right = np.zeros((2, len(steps)))
   down = np.zeros((2, len(steps)))
   left = np.zeros((2, len(steps)))
   up = np.zeros((2, len(steps)))
   

4. 设置图像位置的坐标很简单。 但是，需要注意一个棘手的问题-[::-1]表示法会导致数组元素的顺序颠倒：

   right[0] = steps
   right[1] = 20
   
   down[0] = 360
   down[1] = steps
   
   left[0] = steps[::-1]
   left[1] = 360
   
   up[0] = 20
   up[1] = steps[::-1]
   

5. 我们可以加入路径部分，但是在执行此操作之前，请使用T运算符转置数组，因为它们未正确对齐以进行连接：

   pos = np.concatenate((right.T, down.T, left.T, up.T))
   

6. 在主事件循环中，让时钟以每秒 30 帧的速度计时：

      clock.tick(30)
   

   摇头的屏幕截图如下：

   

   您应该能够观看此动画的电影。 它也是代码包（animation.mp4）的一部分。

   此示例的代码几乎使用了到目前为止我们学到的所有内容，但仍应足够简单以了解：

   import pygame, sys
   from pygame.locals import *
   import numpy as np
   
   pygame.init()
   clock = pygame.time.Clock()
   screen = pygame.display.set_mode((400, 400))
   
   pygame.display.set_caption('Animating Objects')
   img = pygame.image.load('head.jpg')
   
   steps = np.linspace(20, 360, 40).astype(int)
   right = np.zeros((2, len(steps)))
   down = np.zeros((2, len(steps)))
   left = np.zeros((2, len(steps)))
   up = np.zeros((2, len(steps)))
   
   right[0] = steps
   right[1] = 20
   
   down[0] = 360
   down[1] = steps
   
   left[0] = steps[::-1]
   left[1] = 360
   
   up[0] = 20
   up[1] = steps[::-1]
   
   pos = np.concatenate((right.T, down.T, left.T, up.T))
   i = 0
   
   while True: 
      # Erase screen
      screen.fill((255, 255, 255))
   
      if i >= len(pos):
         i = 0
   
      screen.blit(img, pos[i])
      i += 1
   
      for event in pygame.event.get():
         if event.type == QUIT:
            pygame.quit()
            sys.exit()
   
      pygame.display.update()
      clock.tick(30)
   

刚刚发生了什么？

在本节中，我们了解了一些有关动画的知识。 我们了解到的最重要的概念是时钟。 下表描述了我们使用的新函数：

函数	描述
pygame.time.Clock()	这将创建一个游戏时钟。
clock.tick(30)	此函数执行游戏时钟的刻度。 此处，30是每秒的帧数。

matplotlib

matplotlib是一个易于绘制的开源库，我们在第 9 章，“matplotlib 绘图”中了解到。 我们可以将 matplotlib 集成到 Pygame 游戏中并创建各种绘图。

实战时间 – 在 Pygame 中使用 matplotlib

在本秘籍中，我们采用上一节的位置坐标，并对其进行绘制。

1. 要将 matplotlib 与 Pygame 集成，我们需要使用非交互式后端； 否则，默认情况下，matplotlib 将为我们提供一个 GUI 窗口。 我们将导入主要的 matplotlib 模块并调用use()函数。 在导入主 matplotlib 模块之后以及在导入其他 matplotlib 模块之前，立即调用此函数：

   import matplotlib as mpl
   
   mpl.use("Agg")
   

2. 我们可以在 matplotlib 画布上绘制非交互式绘图。 创建此画布需要导入，创建图形和子图。 将数字指定为3乘3英寸大。 在此秘籍的末尾可以找到更多详细信息：

   import matplotlib.pyplot as plt
   import matplotlib.backends.backend_agg as agg
   
   fig = plt.figure(figsize=[3, 3])
   ax = fig.add_subplot(111)
   canvas = agg.FigureCanvasAgg(fig)
   

3. 在非交互模式下，绘制数据比在默认模式下复杂一些。 由于我们需要重复绘图，因此在函数中组织绘图代码是有意义的。 Pygame 最终在画布上绘制了绘图。 画布为我们的设置增加了一些复杂性。 在此示例的末尾，您可以找到有关这些函数的更多详细说明：

   def plot(data):
      ax.plot(data)
      canvas.draw()
      renderer = canvas.get_renderer()
   
      raw_data = renderer.tostring_rgb()
      size = canvas.get_width_height()
   
      return pygame.image.fromstring(raw_data, size, "RGB")
   

   下面的屏幕截图显示了正在运行的动画。 您还可以在代码包（matplotlib.mp4）和 YouTube 上查看截屏视频。

   

   更改后，我们将获得以下代码：

   import pygame, sys
   from pygame.locals import *
   import numpy as np
   import matplotlib as mpl
   
   mpl.use("Agg")
   
   import matplotlib.pyplot as plt
   import matplotlib.backends.backend_agg as agg
   
   fig = plt.figure(figsize=[3, 3])
   ax = fig.add_subplot(111)
   canvas = agg.FigureCanvasAgg(fig)
   
   def plot(data):
      ax.plot(data)
      canvas.draw()
      renderer = canvas.get_renderer()
   
      raw_data = renderer.tostring_rgb()
      size = canvas.get_width_height()
   
      return pygame.image.fromstring(raw_data, size, "RGB")
   
   pygame.init()
   clock = pygame.time.Clock()
   screen = pygame.display.set_mode((400, 400))
   
   pygame.display.set_caption('Animating Objects')
   img = pygame.image.load('head.jpg')
   
   steps = np.linspace(20, 360, 40).astype(int)
   right = np.zeros((2, len(steps)))
   down = np.zeros((2, len(steps)))
   left = np.zeros((2, len(steps)))
   up = np.zeros((2, len(steps)))
   
   right[0] = steps
   right[1] = 20
   
   down[0] = 360
   down[1] = steps
   
   left[0] = steps[::-1]
   left[1] = 360
   
   up[0] = 20
   up[1] = steps[::-1]
   
   pos = np.concatenate((right.T, down.T, left.T, up.T))
   i = 0
   history = np.array([])
   surf = plot(history)
   
   while True: 
      # Erase screen
      screen.fill((255, 255, 255))
   
      if i >= len(pos):
         i = 0
         surf = plot(history)
   
      screen.blit(img, pos[i])
      history = np.append(history, pos[i])
      screen.blit(surf, (100, 100))
   
      i += 1
   
      for event in pygame.event.get():
         if event.type == QUIT:
            pygame.quit()
            sys.exit()
   
      pygame.display.update()
      clock.tick(30)
   

刚刚发生了什么？

下表解释了绘图相关函数：

函数	描述
mpl.use("Agg")	此函数指定使用非交互式后端
plt.figure(figsize=[3, 3])	此函数创建一个3 x 3英寸的图形
agg.FigureCanvasAgg(fig)	此函数在非交互模式下创建画布
canvas.draw()	此函数在画布上绘制
canvas.get_renderer()	此函数为画布提供渲染器

表面像素

Pygame surfarray模块处理 Pygame Surface对象与 NumPy 数组之间的转换 。 您可能还记得，NumPy 可以快速有效地处理大型数组。

实战时间 – 用 NumPy 访问表面像素数据

在本节中，我们将平铺一个小图像以填充游戏屏幕。

1. array2d()函数将像素复制到二维数组中（对于三维数组也有类似的功能）。 将头像图像中的像素复制到数组中：

   pixels = pygame.surfarray.array2d(img)
   

2. 使用数组的shape属性从像素数组的形状创建游戏屏幕。 在两个方向上将屏幕放大七倍：

   X = pixels.shape[0] * 7
   Y = pixels.shape[1] * 7
   screen = pygame.display.set_mode((X, Y))
   

3. 使用 NumPy tile()函数可以轻松平铺图像。 数据需要转换为整数值，因为 Pygame 将颜色定义为整数：

   new_pixels = np.tile(pixels, (7, 7)).astype(int)
   

4. surfarray模块具有特殊函数blit_array()在屏幕上显示数组：

   pygame.surfarray.blit_array(screen, new_pixels)
   

   

   以下代码执行图像的平铺：

   import pygame, sys
   from pygame.locals import *
   import numpy as np
   
   pygame.init()
   img = pygame.image.load('head.jpg')
   pixels = pygame.surfarray.array2d(img)
   X = pixels.shape[0] * 7
   Y = pixels.shape[1] * 7
   screen = pygame.display.set_mode((X, Y))
   pygame.display.set_caption('Surfarray Demo')
   new_pixels = np.tile(pixels, (7, 7)).astype(int)
   
   while True: 
      screen.fill((255, 255, 255))
      pygame.surfarray.blit_array(screen, new_pixels)
   
      for event in pygame.event.get():
         if event.type == QUIT:
            pygame.quit()
            sys.exit()
   
      pygame.display.update()
   

刚刚发生了什么？

以下是我们使用的新函数和属性的简要说明：

函数	描述
pygame.surfarray.array2d(img)	此函数将像素数据复制到二维数组中
pygame.surfarray.blit_array(screen, new_pixels)	此函数在屏幕上显示数组值

人工智能

通常，我们需要模仿游戏中的智能行为。 scikit-learn项目旨在提供一种用于机器学习的 API，而我最喜欢的是其精美的文档。 我们可以使用操作系统的包管理器来安装scikit-learn，尽管此选项可能有效或无效，具体取决于您的操作系统，但这应该是最方便的方法。 Windows 用户只需从项目网站下载安装程序即可。 在 Debian 和 Ubuntu 上，该项目称为python-sklearn。 在 MacPorts 上，这些端口称为py26-scikits-learn和py27-scikits-learn。 我们也可以从源代码或使用easy_install安装。 PythonXY， Enthought 和 NetBSD。

我们可以通过在命令行中键入来安装scikit-learn：

$ [sudo] pip install -U scikit-learn

我们也可以键入以下内容而不是前一行：

$ [sudo] easy_install -U scikit-learn

由于权限的原因，这可能无法正常工作，因此您可能需要在命令前面放置sudo或以管理员身份登录。

实战时间 – 点的聚类

我们将生成一些随机点并将它们聚类，这意味着彼此靠近的点将放入同一簇中。 这只是scikit-learn可以应用的许多技术之一。聚类是一种机器学习算法，旨在基于相似度对项目进行分组。 接下来，我们将计算平方亲和矩阵。亲和度矩阵是包含亲和度值的矩阵：例如，点之间的距离。 最后，我们将这些点与[​​HTG2]中的AffinityPropagation类聚类。

1. 在400 x 400像素的正方形内生成 30 个随机点位置：

   positions = np.random.randint(0, 400, size=(30, 2))
   

2. 使用到原点的欧式距离作为亲和度度量来计算亲和度矩阵：

   positions_norms = np.sum(positions ** 2, axis=1)
   S = - positions_norms[:, np.newaxis] - positions_norms[np.newaxis, :] + 2 * np.dot(positions, positions.T)
   

3. 给AffinityPropagation类上一步的结果。 此类使用适当的群集编号标记点：

   aff_pro = sklearn.cluster.AffinityPropagation().fit(S)
   labels = aff_pro.labels_
   

4. 为每个群集绘制多边形。 涉及的函数需要点列表，颜色（将其绘制为红色）和表面：

   pygame.draw.polygon(screen, (255, 0, 0), polygon_points[i])
   

   结果是每个群集的一堆多边形，如下图所示：

   

   群集示例代码如下：

   import numpy as np
   import sklearn.cluster
   import pygame, sys
   from pygame.locals import *
   
   np.random.seed(42)
   positions = np.random.randint(0, 400, size=(30, 2))
   
   positions_norms = np.sum(positions ** 2, axis=1)
   S = - positions_norms[:, np.newaxis] - positions_norms[np.newaxis, :] + 2 * np.dot(positions, positions.T)
   
   aff_pro = sklearn.cluster.AffinityPropagation().fit(S)
   labels = aff_pro.labels_
   
   polygon_points = []
   
   for i in xrange(max(labels) + 1):
      polygon_points.append([])
   
   # Sorting points by cluster
   for label, position in zip(labels, positions):
      polygon_points[labels[i]].append(positions[i])
   
   pygame.init()
   screen = pygame.display.set_mode((400, 400))
   
   while True: 
      for point in polygon_points:
         pygame.draw.polygon(screen, (255, 0, 0), point)
   
      for event in pygame.event.get():
         if event.type == QUIT:
            pygame.quit()
            sys.exit()
   
      pygame.display.update()
   

刚刚发生了什么？

下表更详细地描述了人工智能示例中最重要的行 ：

函数	描述
sklearn.cluster.AffinityPropagation().fit(S)	此函数创建AffinityPropagation对象，并使用相似性矩阵执行拟合
pygame.draw.polygon(screen, (255, 0, 0), point)	给定表面，颜色（在这种情况下为红色）和点列表，此函数绘制多边形

OpenGL 和 Pygame

OpenGL 为二维和三维计算机图形指定了 API。 API 由函数和常量组成。 我们将专注于名为PyOpenGL的 Python 实现。 使用以下命令安装PyOpenGL：

$ [sudo] pip install PyOpenGL PyOpenGL_accelerate

您可能需要具有 root 访问权限才能执行此命令。 对应于easy_install的命令如下：

$ [sudo] easy_install PyOpenGL PyOpenGL_accelerate

实战时间 – 绘制 Sierpinski 地毯

为了演示的目的，我们将使用 OpenGL 绘制一个 Sierpinski 地毯，也称为 Sierpinski 三角形或 Sierpinski 筛子。 这是由数学家 Waclaw Sierpinski 创建的三角形形状的分形图案。 三角形是通过递归且原则上是无限的过程获得的。

1. 首先，首先初始化一些与 OpenGL 相关的原语。 这包括设置显示模式和背景颜色。 本节末尾提供逐行说明：

   def display_openGL(w, h):
     pygame.display.set_mode((w,h), pygame.OPENGL|pygame.DOUBLEBUF)
   
     glClearColor(0.0, 0.0, 0.0, 1.0)
     glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)
   
     gluOrtho2D(0, w, 0, h)
   

2. 该算法要求我们显示点，越多越好。 首先，我们将绘图颜色设置为红色。 其次，我们定义一个三角形的顶点（我称它们为点）。 然后，我们定义随机索引，该随机索引将用于选择三个三角形顶点之一。 我们在中间的某个地方随机选择一个点，实际上并不重要。 之后，在上一个点和随机选取的一个顶点之间的一半处绘制点。 最后，刷新结果：

       glColor3f(1.0, 0, 0)
       vertices = np.array([[0, 0], [DIM/2, DIM], [DIM, 0]])
       NPOINTS = 9000
       indices = np.random.random_integers(0, 2, NPOINTS)
       point = [175.0, 150.0]
   
       for i in xrange(NPOINTS):
          glBegin(GL_POINTS)
          point = (point + vertices[indices[i]])/2.0
          glVertex2fv(point)
          glEnd()
   
       glFlush()
   

   Sierpinski 三角形如下所示：

   

   带有所有导入的完整 Sierpinski 垫圈演示代码如下：

   import pygame
   from pygame.locals import *
   import numpy as np
   
   from OpenGL.GL import *
   from OpenGL.GLU import *
   
   def display_openGL(w, h):
     pygame.display.set_mode((w,h), pygame.OPENGL|pygame.DOUBLEBUF)
   
     glClearColor(0.0, 0.0, 0.0, 1.0)
     glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)
   
     gluOrtho2D(0, w, 0, h)
   
   def main():
       pygame.init()
       pygame.display.set_caption('OpenGL Demo')
       DIM = 400
       display_openGL(DIM, DIM)
       glColor3f(1.0, 0, 0)
       vertices = np.array([[0, 0], [DIM/2, DIM], [DIM, 0]])
       NPOINTS = 9000
       indices = np.random.random_integers(0, 2, NPOINTS)
       point = [175.0, 150.0]
   
       for i in xrange(NPOINTS):
          glBegin(GL_POINTS)
          point = (point + vertices[indices[i]])/2.0
          glVertex2fv(point)
          glEnd()
   
       glFlush()
       pygame.display.flip()
   
       while True:
           for event in pygame.event.get():
               if event.type == QUIT:
                   return
   
   if __name__ == '__main__':
     main()
   

刚刚发生了什么？

如所承诺的，以下是该示例最重要部分的逐行说明：

函数	描述
pygame.display.set_mode((w,h), pygame.OPENGL|pygame.DOUBLEBUF)	此函数将显示模式设置为所需的宽度，高度和 OpenGL 显示
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)	此函数使用遮罩清除缓冲区。 在这里，我们清除颜色缓冲区和深度缓冲区位
gluOrtho2D(0, w, 0, h)	此函数定义二维正交投影矩阵，其坐标为左，右，上和下剪切平面
glColor3f(1.0, 0, 0)	此函数使用 RGB 的三个浮点值（红色，绿色，蓝色）定义当前图形颜色。 在这种情况下，我们将以红色绘制
glBegin(GL_POINTS)	此函数定义了图元的顶点或图元的组。 这里的原语是点
glVertex2fv(point)	此函数在给定顶点的情况下渲染点
glEnd()	此函数关闭以glBegin()开头的一段代码
glFlush()	此函数强制执行 GL 命令

使用 Pygame 的模拟游戏

作为最后一个示例，我们将使用 Conway 的生命游戏模拟生命。 最初的生命游戏是基于一些基本规则。 我们从二维正方形网格上的随机配置开始。 网格中的每个单元可以是死的或活着的。 此状态取决于小区的邻居。 您可以在这个页面上详细了解规则。在每个时间步上，都会发生以下转换：

1. 少于两个活邻居的活细胞死亡。
2. 具有两个或三个活邻居的活细胞可以存活到下一代。
3. 具有三个以上活邻居的活细胞死亡。
4. 具有恰好三个活邻居的死细胞会成为活细胞。

卷积可用于求值游戏的基本规则。 卷积过程需要 SciPy 包。

实战时间 – 模拟生命

以下代码是生命游戏的实现 ，并进行了一些修改：

• 用鼠标单击一次会画一个十字，直到我们再次单击
• r键可将网格重置为随机状态
• b键根据鼠标位置创建块
• g键创建滑翔机

代码中最重要的数据结构是一个二维数组，其中包含游戏屏幕上像素的颜色值。 该数组用随机值初始化，然后针对游戏循环的每次迭代重新计算。 在下一部分中找到有关所涉及函数的更多信息。

1. 要求值规则，请使用卷积，如下所示：

   def get_pixar(arr, weights):
     states = ndimage.convolve(arr, weights, mode='wrap')
   
     bools = (states == 13) | (states == 12 ) | (states == 3)
   
     return bools.astype(int)
   

2. 使用我们在第 2 章，“从 NumPy 基础知识开始”中学习的基本索引技巧来画十字：

   def draw_cross(pixar):
      (posx, posy) = pygame.mouse.get_pos()
      pixar[posx, :] = 1
      pixar[:, posy] = 1
   

3. 用随机值初始化网格：

   def random_init(n):
      return np.random.random_integers(0, 1, (n, n))
   

   以下是完整的代码：

   from __future__ import print_function
   import os, pygame
   from pygame.locals import *
   import numpy as np
   from scipy import ndimage
   
   def get_pixar(arr, weights):
     states = ndimage.convolve(arr, weights, mode='wrap')
   
     bools = (states == 13) | (states == 12 ) | (states == 3)
   
     return bools.astype(int)
   
   def draw_cross(pixar):
      (posx, posy) = pygame.mouse.get_pos()
      pixar[posx, :] = 1
      pixar[:, posy] = 1
   
   def random_init(n):
      return np.random.random_integers(0, 1, (n, n))
   
   def draw_pattern(pixar, pattern):
        print(pattern)
   
        if pattern == 'glider':
         coords = [(0,1), (1,2), (2,0), (2,1), (2,2)]
        elif pattern == 'block':
         coords = [(3,3), (3,2), (2,3), (2,2)]
        elif pattern == 'exploder':
         coords = [(0,1), (1,2), (2,0), (2,1), (2,2), (3,3)]
        elif pattern == 'fpentomino':
         coords = [(2,3),(3,2),(4,2),(3,3),(3,4)]
   
        pos = pygame.mouse.get_pos()
   
        xs = np.arange(0, pos[0], 10)
        ys = np.arange(0, pos[1], 10)
   
        for x in xs:
           for y in ys:
              for i, j in coords:
                  pixar[x + i, y + j] = 1
   
   def main():
       pygame.init ()
       N = 400
       pygame.display.set_mode((N, N))
       pygame.display.set_caption("Life Demo")
   
       screen = pygame.display.get_surface()
   
       pixar = random_init(N)
       weights = np.array([[1,1,1], [1,10,1], [1,1,1]])
   
       cross_on = False
   
       while True:
          pixar = get_pixar(pixar, weights)
   
          if cross_on:
             draw_cross(pixar)
   
          pygame.surfarray.blit_array(screen, pixar * 255 ** 3)
          pygame.display.flip()
   
          for event in pygame.event.get():
            if event.type == QUIT:
                return
            if event.type == MOUSEBUTTONDOWN:
               cross_on = not cross_on
            if event.type == KEYDOWN:
               if event.key == ord('r'):
                  pixar = random_init(N)
                  print("Random init")
               if event.key == ord('g'):
                  draw_pattern(pixar, 'glider')
               if event.key == ord('b'):
                  draw_pattern(pixar, 'block')
               if event.key == ord('e'):
                  draw_pattern(pixar, 'exploder')
               if event.key == ord('f'):
                  draw_pattern(pixar, 'fpentomino')
   
   if __name__ == '__main__':
       main()
   

   您应该能够从代码包（life.mp4）或 YouTube 上观看截屏视频。 正在运行的游戏的屏幕截图如下：

   

刚刚发生了什么？

我们使用了一些 NumPy 和 SciPy 函数，这些函数需要说明：

函数	描述
ndimage.convolve(arr, weights, mode='wrap')	此函数在包装模式下使用权重将卷积运算应用于给定数组。 该模式与数组边界有关。
bools.astype(int)	此函数将布尔数组转换为整数。
np.arange(0, pos[0], 10)	此函数以 10 为步长创建一个从 0 到pos[0]的数组。因此，如果pos[0]等于 1000，我们将得到 0、10、20，… 990。

总结

您可能会发现本书中提到 Pygame 有点奇怪。 但是，阅读本章后，我希望您意识到 NumPy 和 Pygame 可以很好地结合在一起。 毕竟，游戏涉及大量计算，因此 NumPy 和 SciPy 是理想的选择，并且它们还需要scikit-learn中提供的人工智能功能。 无论如何，制作游戏都很有趣，我们希望最后一章相当于十道菜后的精美甜点或咖啡！ 如果您仍然渴望更多，请查看《NumPy Cookbook 第二版》， Ivan Idris，Packt Publishing，在本书的基础上以最小的重叠为基础。

附录 A：小测验答案

第 1 章，NumPy 快速入门

小测验 – arange()函数的功能

arange(5)做什么？	它创建一个 NumPy 数组，其值为从 0-4 创建的 NumPy 数组的值，0、1、2、3 和 4

第 2 章，从 NumPy 基本原理开始

小测验 – ndarray的形状

ndarray的形状如何存储？	它存储在一个元组中

第 3 章，熟悉常用函数

小测验 - 计算加权平均值

哪个函数返回数组的加权平均值？	average

第 4 章，为您带来便利的便利函数

小测验 - 计算协方差

哪个函数返回两个数组的协方差？	cov

第 5 章，使用矩阵和ufunc

小测验 – 使用字符串定义矩阵

mat和bmat函数接受的字符串中的行分隔符是什么？	分号

第 6 章，深入探索 NumPy 模块

小测验 - 创建矩阵

哪个函数可以创建矩阵？	mat

第 7 章，探索特殊例程

小测验 - 生成随机数

哪个 NumPy 模块处理随机数？	random

第 8 章，通过测试确保质量

小测验 - 指定小数精度

assert_almost_equal函数的哪个参数指定小数精度？	decimal

第 9 章，matplotlib 绘图

小测验 – plot()函数

plot函数有什么作用？	它既不执行 1、2 也不执行 3

第 10 章，当 NumPy 不够用时 – Scipy 和更多

小测验 - 加载.mat文件

哪个函数加载.mat文件？	loadmat

附录 B：其他在线资源

本附录包含指向相关网站的链接。

Python 教程

1. Think Python 中文第二版↗
2. 笨办法学 Python · 续 中文版
3. PythonSpot 中文系列教程
4. PythonBasics 中文系列教程
5. PythonGuru 中文系列教程
6. Python 分布式计算↗

数学教程

1. MIT 公开课课本/笔记
   1. MIT 18.06 线性代数笔记↗
   2. MIT 18.03 写给初学者的微积分

数据科学文档

1. Numpy 技术栈中文文档
   1. NumPy 1.11 中文文档⭐
   2. Pandas 0.19.2 中文文档⭐
   3. Matplotlib 2.0 中文文档⭐
   4. statsmodels 中文文档⭐
   5. seaborn 0.9 中文文档⭐
2. SimuPy 中文文档↗

数据科学教程

1. 斯坦福公开课课本/笔记
   1. 斯坦福 STATS60 课本：21 世纪的统计思维⭐
   2. 斯坦福博弈论中文笔记⭐
2. UCB 公开课课本/笔记
   1. UCB Data8 课本：计算与推断思维⭐
   2. UCB Prob140 课本：面向数据科学的概率论⭐
   3. UCB DS100 课本：数据科学的原理与技巧⭐
3. ApacheCN 数据科学译文集⭐
4. TutorialsPoint NumPy 教程
5. 复杂性思维 中文第二版
6. 利用 Python 进行数据分析 · 第 2 版
7. fast.ai 数值线性代数讲义 v2
8. Pandas Cookbook 带注释源码
9. 数据科学 IPython 笔记本
10. UCSD COGS108 数据科学实战中文笔记
11. USF MSDS501 计算数据科学中文讲义
12. 数据可视化的基础知识
13. Joyful Pandas↗

附录 C：NumPy 函数的参考

本附录包含有用的 NumPy 函数及其说明的列表。

• numpy.apply_along_axis（func1d, axis, arr, *args）：沿arr的一维切片应用函数func1d。

• numpy.arange([start,] stop[, step,], dtype=None)：创建一个 NumPy 数组，它在指定范围内均匀间隔。

• numpy.argsort(a, axis=-1, kind='quicksort', order=None)：返回对输入数组进行排序的索引。

• numpy.argmax(a, axis=None)：返回沿轴的最大值的索引。

• numpy.argmin(a, axis=None)：返回沿轴的最小值的索引。

• numpy.argwhere(a)：查找非零元素的索引。

• numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0)：从类似数组的序列（例如 Python 列表）创建 NumPy 数组。

• numpy.testing.assert_allclose((actual, desired, rtol=1e-07, atol=0, err_msg='', verbose=True)：如果两个对象在指定的精度下不相等，则引发错误。

• numpy.testing.assert_almost_equal()：如果两个数字在指定的精度下不相等，则引发异常。

• numpy.testing.assert_approx_equal()：如果两个数字在某个有效数字下不相等，则引发异常。

• numpy.testing.assert_array_almost_equal()：如果两个数组在指定的精度下不相等，则引发异常。

• numpy.testing.assert_array_almost_equal_nulp(x, y, nulp=1)：将数组与其最低精度单位（ULP）。

• numpy.testing.assert_array_equal()：如果两个数组不相等，则引发异常。

• numpy.testing.assert_array_less()：如果两个数组的形状不同，并且第一个数组的元素严格小于第二个数组的元素，则会引发异常。

• numpy.testing.assert_array_max_ulp(a, b, maxulp=1, dtype=None)：确定数组元素最多相差 ULP 的指定数量。

• numpy.testing.assert_equal()：测试两个 NumPy 数组是否相等。

• numpy.testing.assert_raises()：如果使用定义的参数调用的可调用对象未引发指定的异常，则失败。

• numpy.testing.assert_string_equal()：断言两个字符串相等。

• numpy.testing.assert_warns()：如果未引发指定的警告，则失败。

• numpy.bartlett(M)：返回带有M点的 Bartlett 窗口。 此窗口类似于三角形窗口。

• numpy.random.binomial(n, p, size=None)：从二项分布中抽取随机样本。

• numpy.bitwise_and(x1, x2[, out])：计算数组的按位AND。

• numpy.bitwise_xor(x1, x2[, out])：计算数组的按位XOR。

• numpy.blackman(M)：返回一个具有M点的布莱克曼窗口，该窗口接近最佳值，并且比凯撒窗口差。

• numpy.column_stack(tup)：堆叠以元组列形式提供的一维数组 。

• numpy.concatenate ((a1, a2, ...), axis=0)：将数组序列连接在一起。

• numpy.convolve(a, v, mode='full')：计算一维数组的线性卷积。

• numpy.dot(a, b, out=None)：计算两个数组的点积。

• numpy.diff(a, n=1, axis=-1)：计算给定轴的 N 阶差。

• numpy.dsplit(ary, indices_or_sections)：沿着第三轴将数组拆分为子数组。

• numpy.dstack(tup)：沿第三轴堆叠以元组形式给出的数组。

• numpy.eye(N, M=None, k=0, dtype=<type 'float'>)：返回单位矩阵。

• numpy.extract(condition, arr)：使用条件选择数组的元素。

• numpy.fft.fftshift(x, axes=None)：将信号的零频率分量移到频谱的中心。

• numpy.hamming(M)：返回带有M点的汉明窗口。

• numpy.hanning(M)：返回具有M点的汉宁窗口。

• numpy.hstack(tup)：水平堆叠以元组形式给出的数组。

• numpy.isreal(x)：返回一个布尔数组，其中True对应于输入数组的实数（而不是复数）元素。

• numpy.kaiser(M, beta)：对于给定的beta参数，返回带有M点的凯撒窗口。

• numpy.load(file, mmap_mode=None)：从.npy，.npz，或腌制中加载 NumPy 数组或腌制对象。 内存映射的数组存储在文件系统中，不必完全加载到内存中。 这对于大型数组尤其有用。

• numpy.loadtxt(fname, dtype=<type 'float'>, comments='#', delimiter=None, converters=None, skiprows=0, usecols=None, unpack=False, ndmin=0)：将文本文件中的数据加载到 NumPy 数组中。

• numpy.lexsort (keys, axis=-1)：使用多个键进行排序。

• numpy.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)：返回在间隔内均匀间隔的数字。

• numpy.max(a, axis=None, out=None, keepdims=False)：沿轴返回数组的最大值。

• numpy.mean(a, axis=None, dtype=None, out=None, keepdims=False)：沿给定轴计算算术平均值。

• numpy.median(a, axis=None, out=None, overwrite_input=False)：沿给定轴计算中位数 。

• numpy.meshgrid(*xi, **kwargs)：返回坐标向量的坐标矩阵。 例如：

  In: numpy.meshgrid([1, 2], [3, 4])
  Out:
  [array([[1, 2],
          [1, 2]]), array([[3, 3],
          [4, 4]])]
  

• numpy.min(a, axis=None, out=None, keepdims=False)：沿轴返回数组的最小值。

• numpy.msort(a)：返回沿第一轴排序的数组的副本。

• numpy.nanargmax(a, axis=None)：返回给定一个忽略 NaN 的轴的最大值的索引。

• numpy.nanargmin(a, axis=None)：返回给定的轴的最小值索引，忽略 NaN。

• numpy.nonzero(a)：返回非零数组元素的索引。

• numpy.ones(shape, dtype=None, order='C')：创建指定形状和数据类型的 NumPy 数组，包含 1s。

• numpy.piecewise(x, condlist, funclist, *args, **kw)：分段求值函数。

• numpy.polyder(p, m=1)：将多项式微分为给定阶数。

• numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)：执行最小二乘多项式拟合。

• numpy.polysub(a1, a2)：减去多项式。

• numpy.polyval(p, x)：以指定值求值多项式。

• numpy.prod(a, axis=None, dtype=None, out=None, keepdims=False)：返回指定轴上数组元素的乘积 。

• numpy.ravel(a, order='C')：展平数组，或在必要时返回副本。

• numpy.reshape(a, newshape, order='C')：更改 NumPy 数组的形状 。

• numpy.row_stack(tup)：逐行堆叠数组。

• numpy.save(file, arr)：以 NumPy .npy格式将 NumPy 数组保存到文件中。

• numpy.savetxt(fname, X, fmt='%.18e', delimiter=' ', newline='\n', header='', footer='', comments='# ')：将 NumPy 数组保存到文本文件。

• numpy.sinc(a)：计算sinc函数。

• numpy.sort_complex(a)：首先以实部，然后是虚部对数组元素进行排序。

• numpy.split(a, indices_or_sections, axis=0)：将数组拆分为子数组。

• numpy.std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False)：沿给定轴返回标准差。

• numpy.take(a, indices, axis=None, out=None, mode='raise')：使用指定的索引从数组中选择元素。

• numpy.vsplit(a, indices_or_sections)：将数组垂直拆分为子数组。

• numpy.vstack(tup)：垂直堆叠数组。

• numpy.where(condition, [x, y])：基于布尔条件从输入数组中选择数组元素。

• numpy.zeros(shape, dtype=float, order='C')：创建指定形状和数据类型的 NumPy 数组，其中包含零。