目录
1. 前言
1.1 关于作者
1.2 关于本书
1.3 许可
2. 简介
numpy 的一切都是围绕着向量化的，所以你要熟悉这些概念，「向量」（vectors）、「数组」（arrays）、「视图」（views）、「」（ufuncs）。
| 属性 | 描述 |
|---|---|
| ndim | 数组的维数（尺度） |
| shape | 数组的大小。这是一个整数元组，表示每个维度中数组的大小。对于具有 n 行和 m 列的矩阵，shape 将是(n, m) |
| size | 数组的元素总数。等于 shape 元素的乘积 |
| dtype | 数组中元素类型 |
| itemsize | 数组中每个元素的大小（以字节为单位） |
| data | 数组实际元素的缓冲区 |
| strides | 遍历数组时每个维度中需要步进的字节数，可提升效率 |
>>> import numpy as np
>>> Z = np.arange(9).reshape(3, 3).astype(np.float16)
>>> Z
array([[0., 1., 2.],
[3., 4., 5.],
[6., 7., 8.]], dtype=float16)
>>> Z.ndim
2
>>> Z.shape
(3, 3)
>>> Z.size
9
>>> Z.dtype
dtype('float16')
>>> Z.itemsize
2
>>> Z.data
<memory at 0x115385bb0>
>>> Z.strides
(6, 2)
3. 数组剖析
内容
3.1 简介
正如序言中所解释的，你应该对 numpy 有一个基本的经验来阅读这本书。如果不是这样的话，你最好在回来之前先从初学者教程开始。因此，我将在这里简要介绍 numpy 数组的基本结构，特别是关于内存布局、视图、拷贝和数据类型。如果你想让你的计算受益于 numpy 哲学，它们是需要理解的关键概念。
让我们考虑一个简单的例子，在这个例子中，我们希望从具有 dtype 为 np.float32 的数组中清除所有值。如何编写最快的操作？下面的语法非常简明（至少对于那些熟悉 numpy 的人来说是这样），但是上面的问题要求找到最快的操作。
>>> Z = np.ones(4*1000000, np.float32)
>>> Z[...] = 0
如果您更仔细地观察 dtype 和数组的大小，可以发现这个数组可以被转换（即查看）成许多其他“兼容”的数据类型。所谓兼容，我的意思是 Z.size*Z.itemsize 可以除以新的 dtype itemsize。
>>> timeit("Z.view(np.float16)[...] = 0", globals())
100 loops, best of 3: 2.72 msec per loop
>>> timeit("Z.view(np.int16)[...] = 0", globals())
100 loops, best of 3: 2.77 msec per loop
>>> timeit("Z.view(np.int32)[...] = 0", globals())
100 loops, best of 3: 1.29 msec per loop
>>> timeit("Z.view(np.float32)[...] = 0", globals())
100 loops, best of 3: 1.33 msec per loop
>>> timeit("Z.view(np.int64)[...] = 0", globals())
100 loops, best of 3: 874 usec per loop
>>> timeit("Z.view(np.float64)[...] = 0", globals())
100 loops, best of 3: 865 usec per loop
>>> timeit("Z.view(np.complex128)[...] = 0", globals())
100 loops, best of 3: 841 usec per loop
>>> timeit("Z.view(np.int8)[...] = 0", globals())
100 loops, best of 3: 630 usec per loop
有趣的是，清除所有值的简明方法并不是最快的。通过将数组转换为更大的数据类型，如 np.float64，我们得到了 25% 的速度系数。但是，通过将数组视为字节数组(np.int8)，我们得到了 50% 的速度因子。这种加速的原因可以在 numpy 内部机制和编译器优化中找到。这个简单的例子说明了 numpy 的哲学，我们将在下面的小节中看到。
3.2 内存布局
numpy 文档对类 ndarray 的定义非常明确：
类 ndarray 的实例由一个连续的一维计算机内存段（由数组或其他对象拥有）与索引方案（将N个整数映射到块中某个项的位置）组成。
换言之，数组主要是一个连续的内存块，其可以使用索引方案(indexing schema)访问。这样的索引方案又由 shape 和数据类型定义，这正是定义新数组时所需要的：
Z = np.arange(9).reshape(3, 3).astype(np.int16)
此外，由于 Z 不是视图，我们可以推导出数组的跨距（strides），它定义了遍历数组时每个维度的字节数。
>>> strides = Z.shape[1] * Z.itemsize, Z.itemsize
>>> strides
(6, 2)
>>> Z.strides
(6, 2)
有了这些信息，我们就知道如何访问特定的项（由索引元祖设计）更加准确的是，如何计算开始和结束的偏移量
offset_start = 0
for i in range(ndim):
offset_start += strides[i] * index[i]
offset_end = offset_start + Z.itemsize
可以用 tobytes 来验证一下上述偏移是否正确：
>>> Z = np.arange(9).reshape(3, 3).astype(np.int16)
>>> index = 1, 1
>>> Z[index].tobytes()
b'\x04\x00'
>>> offset = 0
>>> for i in range(Z.ndim):
... offset + = Z.strides[i]*index[i]
>>> print(Z.tobytes()[offset_start:offset_end]
b'\x04\x00'
这个数组实际上可以从不同的角度考虑（比如：布局）：
元素布局(Item layout)
shape[1]
(=3)
┌───────────┐
┌ ┌───┬───┬───┐ ┐
│ │ 0 │ 1 │ 2 │ │
│ ├───┼───┼───┤ │
shape[0] │ │ 3 │ 4 │ 5 │ │ len(Z)
(=3) │ ├───┼───┼───┤ │ (=3)
│ │ 6 │ 7 │ 8 │ │
└ └───┴───┴───┘ ┘
展平元素布局(Flattened item layout)
┌───┬───┬───┬───┬───┬───┬───┬───┬───┐
│ 0 │ 1 │ 2 │ 3 │ 4 │ 5 │ 6 │ 7 │ 8 │
└───┴───┴───┴───┴───┴───┴───┴───┴───┘
└───────────────────────────────────┘
Z.size
(=9)
内存布局(C order, big endian)
strides[1]
(=2)
┌─────────────────────┐
┌ ┌──────────┬──────────┐ ┐
│ p+00: │ 00000000 │ 00000000 │ │
│ ├──────────┼──────────┤ │
│ p+02: │ 00000000 │ 00000001 │ │ strides[0]
│ ├──────────┼──────────┤ │ (=2x3)
│ p+04 │ 00000000 │ 00000010 │ │
│ ├──────────┼──────────┤ ┘
│ p+06 │ 00000000 │ 00000011 │
│ ├──────────┼──────────┤
Z.nbytes │ p+08: │ 00000000 │ 00000100 │
(=3x3x2) │ ├──────────┼──────────┤
│ p+10: │ 00000000 │ 00000101 │
│ ├──────────┼──────────┤
│ p+12: │ 00000000 │ 00000110 │
│ ├──────────┼──────────┤
│ p+14: │ 00000000 │ 00000111 │
│ ├──────────┼──────────┤
│ p+16: │ 00000000 │ 00001000 │
└ └──────────┴──────────┘
└─────────────────────┘
Z.itemsize
Z.dtype.itemsize
(=2)
如果现在对 Z 取一个切片，那将得到一个基于数组 Z 的视图（view）：
V = Z[::2, ::2]
这样的视图是用 shape、dtype 和 strides 指定的，因为 strides （的存在）不能再仅从 dtype 和 shape 推断（）。
视图的元素布局
shape[1]
(=2)
┌───────────┐
┌ ┌───┬╌╌╌┬───┐ ┐
│ │ 0 │ │ 2 │ │ ┌───┬───┐
│ ├───┼╌╌╌┼───┤ │ │ 0 │ 2 │
shape[0] │ ╎ ╎ ╎ ╎ │ len(Z) → ├───┼───┤
(=2) │ ├───┼╌╌╌┼───┤ │ (=2) │ 6 │ 8 │
│ │ 6 │ │ 8 │ │ └───┴───┘
└ └───┴╌╌╌┴───┘ ┘
展平的元素布局
┌───┬╌╌╌┬───┬╌╌╌┬╌╌╌┬╌╌╌┬───┬╌╌╌┬───┐ ┌───┬───┬───┬───┐
│ 0 │ │ 2 │ ╎ ╎ │ 6 │ │ 8 │ → │ 0 │ 2 │ 6 │ 8 │
└───┴╌╌╌┴───┴╌╌╌┴╌╌╌┴╌╌╌┴───┴╌╌╌┴───┘ └───┴───┴───┴───┘
└─┬─┘ └─┬─┘ └─┬─┘ └─┬─┘
└───┬───┘ └───┬───┘
└───────────┬───────────┘
Z.size
(=4)
内存布局（C order, big endian）
┌ ┌──────────┬──────────┐ ┐ ┐
┌─┤ p+00: │ 00000000 │ 00000000 │ │ │
│ └ ├──────────┼──────────┤ │ strides[1] │
┌─┤ p+02: │ │ │ │ (=4) │
│ │ ┌ ├──────────┼──────────┤ ┘ │
│ └─┤ p+04 │ 00000000 │ 00000010 │ │
│ └ ├──────────┼──────────┤ │ strides[0]
│ p+06: │ │ │ │ (=12)
│ ├──────────┼──────────┤ │
Z.nbytes ─┤ p+08: │ │ │ │
(=8) │ ├──────────┼──────────┤ │
│ p+10: │ │ │ │
│ ┌ ├──────────┼──────────┤ ┘
│ ┌─┤ p+12: │ 00000000 │ 00000110 │
│ │ └ ├──────────┼──────────┤
└─┤ p+14: │ │ │
│ ┌ ├──────────┼──────────┤
└─┤ p+16: │ 00000000 │ 00001000 │
└ └──────────┴──────────┘
└─────────────────────┘
Z.itemsize
Z.dtype.itemsize
(=2)
3.3 视图和拷贝
视图和拷贝（Views and copies）对于优化数值计算是重要的概念。上一节中虽然进行了操作，但整个故事还是有点复杂。
3.3.1 直接和间接访问
首先，需要区分索引 indexing 和花式索引 fancy indexing。indexing 总是返回视图而 fancy indexing 返回拷贝。这种差异很重要，因为在第一种情况下，修改视图会修改源数组，而在第二种情况下则不会。
>>> Z = np.zeros(9)
>>> Z_view = Z[:3]
>>> Z_view[...] = 1
>>> Z
[ 1. 1. 1. 0. 0. 0. 0. 0. 0.]
>>> Z = np.zeros(9)
>>> Z_copy = Z[[0,1,2]]
>>> Z_copy[...] = 1
>>> Z
[ 0. 0. 0. 0. 0. 0. 0. 0. 0.]
因此，如果您需要花式索引，最好保留一份花式索引的副本（尤其是计算复杂的索引）并使用它：
>>> Z = np.zeros(9)
>>> index = [0,1,2]
>>> Z[index] = 1
>>> print(Z)
[ 1. 1. 1. 0. 0. 0. 0. 0. 0.]
如果你不确定索引的结果是视图还是拷贝，可以检查结果的 base 是什么，如果是 None 那么结果是拷贝：
>>> Z = np.random.uniform(0,1,(5,5))
>>> Z1 = Z[:3,:]
>>> Z2 = Z[[0,1,2], :]
>>> print(np.allclose(Z1,Z2))
True
>>> print(Z1.base is Z)
True
>>> print(Z2.base is Z)
False
>>> print(Z2.base is None)
True
请注意，有些 numpy 函数在可能的情况下会返回视图（例如：ravel），而有些总是返回拷贝（例如：flatten）
>>> Z = np.zeros((5, 5))
>>> Z.ravel().base is Z
True
>>> Z[::2. ::2].ravel().base is Z
False
>>> Z.flatten().base is Z
False
3.3.2 临时拷贝
可以像上一节那样显式地创建拷贝，但最常见的情况是隐式创建中间拷贝。在对数组执行某些运算时，会出现这种情况：
>>> X = np.ones(10, dtype=np.int)
>>> Y = np.ones(10, dtypr=np.int)
>>> A = 2*X + 2*Y
在上面的示例中，创建了三个中间数组。一个用于保存 2*X 的结果，一个用于保存 2*Y 的结果，最后一个用于保存 2*X + 2*Y 的结果。在这种特定情况下，数组足够小，这并没有真正的区别。但是，如果数组很大，那么您必须小心处理这些表达式，并想知道是否可以用不同的方法来处理。例如，如果只有最终结果很重要，而且之后不需要 X 或 Y，则另一种解决方案是：
>>> X = np.ones(10, dtype=np.int)
>>> Y = np.ones(10, dtype=np.int)
>>> np.multiply(X, 2, out=X)
>>> np.multiply(Y, 2, out=Y)
>>> np.add(X, Y, out=X)
使用此替代解决方案，并未创建临时数组。问题是，还有许多其他情况需要创建此类拷贝，这会影响下面示例中所示的性能：
>>> X = np.ones(1000000000, dtype=np.int)
>>> Y = np.ones(1000000000, dtype=np.int)
>>> timeit("X = X + 2.0*Y", globals())
100 loops, best of 3: 3.61 ms per loop
>>> timeit("X = X + 2*Y", globals())
100 loops, best of 3: 3.47 ms per loop
>>> timeit("X += 2*Y", globals())
100 loops, best of 3: 2.79 ms per loop
>>> timeit("np.add(X, Y, out=X); np.add(X, Y, out=X)", globals())
1000 loops, best of 3: 1.57 ms per loop
3.4 结论
最后，我们做一个练习。给定两个向量Z1和Z2。我们想知道Z2是否是Z1的视图，如果是，这个视图是什么？
>>> Z1 = np.arange(10)
>>> Z2 = Z1[1:-1:2]
╌╌╌┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬╌╌
Z1 │ 0 │ 1 │ 2 │ 3 │ 4 │ 5 │ 6 │ 7 │ 8 │ 9 │
╌╌╌┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴╌╌
╌╌╌╌╌╌╌┬───┬╌╌╌┬───┬╌╌╌┬───┬╌╌╌┬───┬╌╌╌╌╌╌╌╌╌╌
Z2 │ 1 │ │ 3 │ │ 5 │ │ 7 │
╌╌╌╌╌╌╌┴───┴╌╌╌┴───┴╌╌╌┴───┴╌╌╌┴───┴╌╌╌╌╌╌╌╌╌╌
首先，需要检查 Z2 的 base 是否是 Z1
>>> Z2.base is Z1
True
此时，我们知道 Z2 是 Z1 的视图，这意味着 Z2 可以表示为 Z1 [start:stop:step]。困难在于找到起点、终点和步长。对于步长，我们可以使用任何数组的 strides 属性，该属性给出每个维度中从一个元素到另一个元素所需的字节数。在我们的例子中，由于两个数组都是一维的，我们可以直接比较第一个 strides：
>>> step = Z2.strides[0] // Z1.strides[0]
>>> step
2
下一个困难是找到开始和结束位置。为此，我们可以利用 byte_bounds 方法返回指向数组端点的指针。
byte_bounds(Z1)[0] byte_bounds(Z1)[-1]
↓ ↓
╌╌╌┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬╌╌
Z1 │ 0 │ 1 │ 2 │ 3 │ 4 │ 5 │ 6 │ 7 │ 8 │ 9 │
╌╌╌┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴╌╌
byte_bounds(Z2)[0] byte_bounds(Z2)[-1]
↓ ↓
╌╌╌╌╌╌╌┬───┬╌╌╌┬───┬╌╌╌┬───┬╌╌╌┬───┬╌╌╌╌╌╌╌╌╌╌
Z2 │ 1 │ │ 3 │ │ 5 │ │ 7 │
╌╌╌╌╌╌╌┴───┴╌╌╌┴───┴╌╌╌┴───┴╌╌╌┴───┴╌╌╌╌╌╌╌╌╌╌
>>> offset_start = np.byte_bounds(Z2)[0] - np.byte_bounds(Z1)[0]
>>> offset_start
8
>>> offset_stop = np.byte_bounds(Z2)[-1] - np.byte_bounds(Z1)[-1]
offset_stop
-16
使用 itemsize 并考虑到 offset_stop 为负（Z2 的结束边界在逻辑上小于 Z1 数组的结束边界），将这些偏移转换为索引非常简单。因此，我们需要加上 Z1 的元素大小来获得右端索引。
>>> start = offset_start // Z1.itemsize
>>> stop = Z1.size + offset_stop // Z1.itemsize
>>> start, stop, step
(1, 8, 2)
最后我们测试一下结果：
>>> np.allclose(Z1[start:stop:step], Z2)
True
作为一个练习，您可以通过考虑以下因素来改进第一个简明的实现：
- 负步长
- 多维数组
这个练习的解决方案
4 代码向量化
4.1 简介
代码向量化意味着你要解决的问题本质上是可向量化的，只需要一些小技巧就能使它更快。当然，这并不意味着它是简单或直接的，但至少它不需要完全重新思考您的问题（正如问题向量化一章中的情况一样）。不过，它可能需要一些经验，看看代码是否可向量化。让我们通过一个简单的例子来说明这一点，在这个例子中，我们要对两个整数列表求和。使用纯 Python 的一种简单方法是：
def add_python(Z1, Z2):
return [z1 + z2 for z1, z2 in zip(Z1, Z2)]
这个简单的方案可以很容易的使用 numpy 进行向量化：
def add_numpy(Z1, Z2):
return np.add(Z1, Z2)
毫无疑问，对这两种方法进行基准测试表明，在同数量级下第二种方法是速度最快的一个。
>>> Z1 = random.sample(range(1000), 100)
>>> Z2 = random.sample(range(1000), 100)
>>> timeit("add_python(Z1, Z2)", globals())
1000 loops, best of 3: 68 usec per loop
>>> timeit("add_numpy(Z1, Z2)", globals())
10000 loops, best of 3: 1.14 usec per loop
(PS: 在我运行的时候结果刚好是相反的，黑人问号）
第二种方法不仅速度更快，而且很自然地适应了 Z1 和 Z2 的 shape。这就是我们没有写 Z1 + Z2 的原因，因为如果 Z1 和 Z2 都是列表，它将无法工作。在第一个 Python 方法中，根据两个对象的性质，对内置 + 的解释是不同的，因此，如果我们考虑两个嵌套列表，则会得到以下输出：
>>> Z1 = [[1, 2], [3, 4]]
>>> Z2 = [[5, 6], [7, 8]]
>>> Z1 + Z2
[[1, 2], [3, 4], [5, 6], [7, 8]]
>>> add_python(Z1, Z2)
[[1, 2, 5, 6], [3, 4, 7, 8]]
>>> add_numpy(Z1, Z2)
[[ 6 8]
[10 12]]
第一个方法将两个列表连接在一起，第二个方法将内部列表连接在一起，最后一个方法计算（数字）期望的值。作为练习，您可以重写Python版本，使其接受任何深度的嵌套列表。
4.2 Uniform 向量化
统一矢量化是矢量化的最简单形式，所有元素在每个时间步共享相同的计算，而不需要对任何元素进行特定处理。约翰最早发明的一个例子就是细胞自动机游戏。这些元胞自动机可以方便地看作是一个单元阵列，这些单元通过邻域的概念连接在一起，并且它们的矢量化很简单。首先让我来定义这个游戏，然后我们来看看如何将其矢量化。
图4.1: 织锦芋螺在它的壳上展示了一个细胞自动机模式。
4.2.1 生命游戏
生命游戏是由英国数学家约翰·霍顿·康威在1970年发明的一种细胞自动机。它是细胞自动机最著名的例子。“游戏”实际上是一个零玩家游戏，这意味着它的进化是由初始状态决定的，不需要人类玩家的输入。一个人通过创造一个初始形态并观察它如何演变来与生命游戏互动。
生命游戏的宇宙是一个无限的二维正交方格网格，每个方格都处于两种可能的状态之一，活的或死的。每一个细胞都与其八个相邻细胞相互作用，这八个邻居是水平、垂直或对角相邻的细胞。在每一步中，都会发生以下转变：
- 任何拥有少于两个活邻居的活细胞都会死亡，就好像由于人口不足导致的需要一样。
- 任何有超过三个活邻居的活细胞都会死亡，就像人满为患一样。
- 任何有两个或三个活邻居的活细胞都不会改变，直到下一代。
- 任何只有三个活邻居的死细胞都会成为活细胞。
初始模式构成了系统的“种子”。第一代是通过同时将上述规则应用于种子中的每个细胞而产生的——出生和死亡同时发生，发生这种情况的离散时刻有时称为滴答声。（换言之，每一代都是前一代的纯功能。）这些规则不断重复应用以创造更多的世代。
4.2.2 Python 实现
*注意：
我们本可以使用更高效的 python 数组接口，但使用熟悉的 list 对象更方便。*
在纯 Python 中，我们可以使用代表细胞应该进化的棋盘的列表的列表来编写生命游戏。这样的棋盘将配备 0 边界，通过避免在计算邻居数量时对边界进行特定测试，可以加快速度。
Z = [[0,0,0,0,0,0],
[0,0,0,1,0,0],
[0,1,0,1,0,0],
[0,0,1,1,0,0],
[0,0,0,0,0,0],
[0,0,0,0,0,0]]
考虑到边界，那么计算邻居就很简单了：
4.2.3 Numpy 的实现
4.2.4 练习
4.2.5 源代码
4.2.6 引用
4.3 临时向量化
4.3.1 Python 的实现
4.3.2 Numpy 的实现
4.3.3 更快的 Numpy 实现
4.3.4 练习
4.3.5 源代码
4.3.6 引用
4.4 空间向量化
4.4.1 Boids
4.4.2 Python 的实现
4.4.3 Numpy 的实现
4.4.4 练习
4.4.5 源代码
4.4.6 引用
4.5 结论
http://www.labri.fr/perso/nrougier/from-python-to-numpy
