紧接着上一篇的翻译，开始最后一部分的内容。
高阶基础
广播（broadcasting）规则
Broadcasting 允许通用函数对不完全相同形状的输入以有意义的方式处理。
Broadcasting的第一个规则是，如果所有输入数组不具有相同数量的维度，则“1”将被重复地添加到较小数组的形状，直到所有数组具有相同数量的维度。Broadcasting的第二个规则确保沿着特定维度具有大小为1的数组表现得好像它们具有沿着该维度具有最大形状的数组的大小。假定数组元素的值沿“Broadcasting”数组的该维度相同。
在应用广播规则之后，所有阵列的大小必须匹配。更多细节可以在 Broadcasting 中找到。
花式索引和索引技巧
Numpy 提供了比常规 Python 序列更多的索引方法。除了整形和切片索引外，正如我们前面所看到的，数组可以被整形数组和布尔数组索引。
使用索引数组索引
>>> import numpy as np
>>> a=np.arange(12)**2
>>> i=np.array([1,1,3,8,5])
>>> a[i]
array([ 1, 1, 9, 64, 25], dtype=int32)
>>> j = np.array([[3,4],[9,7]])
>>> a[j]
array([[ 9, 16],
[81, 49]], dtype=int32)
>>> a
array([ 0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121],
dtype=int32)当被索引数组 a 是一个多维数组，单个索引数组指向 a 的第一个维度。以下示例通过使用调色板将标签图像转换为彩色图像来作为举例。
>>> palette=np.array([[0,0,0],[255,0,0],[0,255,0],[0,0,255],[255,255,255]])
>>> image=np.array([[0,1,2,0],[0,3,4,0]])
>>> palette[image] #数组为（2，4，3）为彩色图像
array([[[ 0, 0, 0],
[255, 0, 0],
[ 0, 255, 0],
[ 0, 0, 0]],
[[ 0, 0, 0],
[ 0, 0, 255],
[255, 255, 255],
[ 0, 0, 0]]])也可以给出索引数组是超过一维的。索引数组的每一维必须是相同的形状。
>>> a=np.arange(12).reshape(3,4)
>>> a
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> i=np.array([[0,1],[1,2]]) # a 的第一个维度索引
>>> j=np.array([[2,1],[3,3]]) # a 的第二个维度索引
>>> a[i,j] #i，j 的形状必须相同
array([[ 2, 5],
[ 7, 11]])
>>> a[i,2] #第二个维度只有 2，可以将其充满和第一个维度形状相同
array([[ 2, 6],
[ 6, 10]])
>>> a[:,j] #对于 a 的每一行做同样的操作
array([[[ 2, 1],
[ 3, 3]],
[[ 6, 5],
[ 7, 7]],
[[10, 9],
[11, 11]]])当然，我们可以把 i 和 j 放在一个序列中（列表），然后用列表进行索引。
>>> m=[i,j]
>>> a[m]
array([[ 2, 5],
[ 7, 11]])然而，我们不能将 i 和 j 放入一个数组中，因为这个数组将被解释为索引第一个维度。
>>> s=np.array([i,j])
>>> a[s]
Traceback (most recent call last):
File "<pyshell#35>", line 1, in <module>
a[s]
IndexError: index 3 is out of bounds for axis 0 with size 3
>>>
>>> a[tuple(s)]
array([[ 2, 5],
[ 7, 11]])索引数组的另一个常见用途是搜索时间相关序列的最大值：
>>> time = np.linspace(20,145,5)
>>> data=np.sin(np.arange(20).reshape(5,4)) #四个时间独立的序列
>>> time
array([ 20. , 51.25, 82.5 , 113.75, 145. ])
>>> data
array([[ 0. , 0.84147098, 0.90929743, 0.14112001],
[-0.7568025 , -0.95892427, -0.2794155 , 0.6569866 ],
[ 0.98935825, 0.41211849, -0.54402111, -0.99999021],
[-0.53657292, 0.42016704, 0.99060736, 0.65028784],
[-0.28790332, -0.96139749, -0.75098725, 0.14987721]])
>>> ind=data.argmax(axis=0) #对每一列找出最大值的索引
>>> ind
array([2, 0, 3, 1], dtype=int64)
>>> time_max=time[ind] #索引对应着最大值
>>> time_max
array([ 82.5 , 20. , 113.75, 51.25])
>>> data_max=data[ind,range(data.shape[1])] #得到每一列的最大值
>>> data_max
array([0.98935825, 0.84147098, 0.99060736, 0.6569866 ])
>>> np.all(data_max==data.max(axis=0))
True
还可以使用数组索引作为目标来赋值：
>>> a=np.arange(5)
>>> a
array([0, 1, 2, 3, 4])
>>> a[[1,3,4]]=0
>>> a
array([0, 0, 2, 0, 0])然而，当索引列表包含重复时，赋值完成多次，留下最后一个值：
>>> a=np.arange(5)
>>> a[[0,0,2]]=[1,2,3]
>>> a
array([2, 1, 3, 3, 4])这是合理的，但是观察以下如果想要使用Python +=结构，可能就不尽人意了：
>>> a=np.arange(5)
>>> a[[0,0,2]]+=1
>>> a
array([1, 1, 3, 3, 4])即使 0 在列表索引中出现两次，但是第 0 个元素仅仅增加一次。这是因为Python 要求 “a+=1” 等价于 “a = a+1”.
使用布尔数组索引
当我们用（整数）索引数组索引数组时，我们提供了要选择的索引列表。使用布尔值作为索引时，方法是不同的；我们明确地选择数组中的哪些元素我们想要的，哪些不是。
我们可以想到的布尔索引最自然的方式是使用与原始数组具有相同形状的布尔数组：
>>> a=np.arange(12).reshape(3,4)
>>> b=a>4
>>> b
array([[False, False, False, False],
[False, True, True, True],
[ True, True, True, True]])
>>> a[b]
array([ 5, 6, 7, 8, 9, 10, 11])这个属性在赋值时非常有用：
>>> a[b]=0
>>> a
array([[0, 1, 2, 3],
[4, 0, 0, 0],
[0, 0, 0, 0]])接下来的例子是通过布尔索引来产生图像：
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> def mandelbrot( h,w, maxit=20 ):
... """Returns an image of the Mandelbrot fractal of size (h,w)."""
... y,x = np.ogrid[ -1.4:1.4:h*1j, -2:0.8:w*1j ]
... c = x+y*1j
... z = c
... divtime = maxit + np.zeros(z.shape, dtype=int)
...
... for i in range(maxit):
... z = z**2 + c
... diverge = z*np.conj(z) > 2**2 # who is diverging
... div_now = diverge & (divtime==maxit) # who is diverging now
... divtime[div_now] = i # note when
... z[diverge] = 2 # avoid diverging too much
...
... return divtime
>>> plt.imshow(mandelbrot(400,400))
>>> plt.show()第二种用布尔来索引的方式类似于整形索引；对于数组中的每一个维度，给出 1D 布尔数组，选择想要的切片：
>>> a=np.arange(12).reshape(3,4)
>>> b1=np.array([False,True,True])
>>> b2=np.array([True,False,True,False])
>>> a
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> a[b1,:]
array([[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> a[b1]
array([[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> a[:,b2]
array([[ 0, 2],
[ 4, 6],
[ 8, 10]])
>>> a[b1,b2]
array([ 4, 10])请注意，1D布尔数组的长度必须与你要切片的维度（或轴）的长度一致。在前面的示例中， b1 是rank为1的数组，其长度为3（ a 中行的数量）， b2 （长度4）适合于索引 a 的第二个rank（列）。
ix_() 函数
可以使用 ix_ 函数来组合不同的向量以获得每个n-uplet的结果。例如，如果要计算从向量a、b和c中的取得的所有三元组的所有a + b * c：
>>> a=np.array([2,3,4,5])
>>> b=np.array([8,5,4])
>>> c=np.array([5,4,6,8,3])
>>> ax,bx,cx=np.ix_(a,b,c)
>>> ax
array([[[2]],
[[3]],
[[4]],
[[5]]])
>>> bx
array([[[8],
[5],
[4]]])
>>> cx
array([[[5, 4, 6, 8, 3]]])
>>> ax.shape,bx.shape,cx.shape
((4, 1, 1), (1, 3, 1), (1, 1, 5))
>>> result=ax+bx*cx
>>> result
array([[[42, 34, 50, 66, 26],
[27, 22, 32, 42, 17],
[22, 18, 26, 34, 14]],
[[43, 35, 51, 67, 27],
[28, 23, 33, 43, 18],
[23, 19, 27, 35, 15]],
[[44, 36, 52, 68, 28],
[29, 24, 34, 44, 19],
[24, 20, 28, 36, 16]],
[[45, 37, 53, 69, 29],
[30, 25, 35, 45, 20],
[25, 21, 29, 37, 17]]])
>>> result.shape
(4, 3, 5)
>>> result[3,2,4]
17
>>> a[3]+b[2]*c[4]
17还可以如下实现 reduce ：
>>> def ufunc_reduce(ufct, *vectors):
... vs = np.ix_(*vectors)
... r = ufct.identity
... for v in vs:
... r = ufct(r,v)
... return r然后将其用作：
>>> ufunc_reduce(np.add,a,b,c)
array([[[15, 14, 16, 18, 13],
[12, 11, 13, 15, 10],
[11, 10, 12, 14, 9]],
[[16, 15, 17, 19, 14],
[13, 12, 14, 16, 11],
[12, 11, 13, 15, 10]],
[[17, 16, 18, 20, 15],
[14, 13, 15, 17, 12],
[13, 12, 14, 16, 11]],
[[18, 17, 19, 21, 16],
[15, 14, 16, 18, 13],
[14, 13, 15, 17, 12]]])与正常的ufunc.reduce相比，这个版本的reduce的优点是它使用Broadcasting规则，以避免创建参数数组输出的大小乘以向量的数量。
线性代数
让我们继续。基本的线性代数就在这儿。
简单的数组操作
有关更多信息，请参阅numpy目录中的linalg.py
>>> import numpy as np
>>> a = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> print(a)
[[ 1. 2.]
[ 3. 4.]]
>>> a.transpose() # 矩阵转置
array([[ 1., 3.],
[ 2., 4.]])
>>> np.linalg.inv(a) #矩阵求逆
array([[-2. , 1. ],
[ 1.5, -0.5]])
>>> u = np.eye(2) # 单位矩阵，对角线为1，其余为0
>>> u
array([[ 1., 0.],
[ 0., 1.]])
>>> j = np.array([[0.0, -1.0], [1.0, 0.0]])
>>> j @ j # 矩阵乘积
array([[-1., 0.],
[ 0., -1.]])
>>> np.trace(u) # 求迹
2.0
>>> y = np.array([[5.], [7.]])
>>> np.linalg.solve(a, y) #求解方程组
array([[-3.],
[ 4.]])
>>> np.linalg.eig(j) #求特征值和特征向量，前一个是特征向量，后一个是特征值
(array([ 0.+1.j, 0.-1.j]), array([[ 0.70710678+0.j , 0.70710678-0.j ],
[ 0.00000000-0.70710678j, 0.00000000+0.70710678j]]))技巧和提示
这里给出一些简短而有用的技巧。
“自动”重定义数组形状
要更改数组的大小，你可以省略其中一个size，它将被自动推导出来：
>>> a = np.arange(30)
>>> a.shape = 2,-1,3 # -1 means "whatever is needed"
>>> a.shape
(2, 5, 3)
>>> a
array([[[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[ 9, 10, 11],
[12, 13, 14]],
[[15, 16, 17],
[18, 19, 20],
[21, 22, 23],
[24, 25, 26],
[27, 28, 29]]])向量堆叠
我们如何从一个相同大小的行向量列表构造一个二维数组？在MATLAB中，这很容易：如果x和y是两个长度相同的向量，那么只需要 m=[x;y] 。在NumPy中，这通过函数 column_stack ，dstack ，hstack 和 vstack 工作，具体取决于要做什么堆叠。例如：
x = np.arange(0,10,2) # x=([0,2,4,6,8])
y = np.arange(5) # y=([0,1,2,3,4])
m = np.vstack([x,y]) # m=([[0,2,4,6,8],
# [0,1,2,3,4]])
xy = np.hstack([x,y]) # xy =([0,2,4,6,8,0,1,2,3,4])这些函数在超过两维时背后的逻辑可能很奇怪。
另见：
NumPy for Matlab users直方图
NumPy的 histogram 函数应用于一个数组，并返回一对向量：数组的histogram和向量的bin。注意： matplotlib 也具有构建histograms的函数（在Matlab中称为 hist ），它与NumPy中的不同。主要区别是 pylab.hist 自动绘制histogram，而 numpy.histogram 仅生成数据。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> # 构建一个向量，长度为10000，分布是均值为2，方差为0.5**2的正太分布
>>> mu, sigma = 2, 0.5
>>> v = np.random.normal(mu,sigma,10000)
>>> # 用50个bin来绘画这个分布
>>> plt.hist(v, bins=50, normed=1) # matplotlib 版本 (plot)
>>> plt.show()>>> # Compute the histogram with numpy and then plot it
>>> (n, bins) = np.histogram(v, bins=50, normed=True) # NumPy 版本 (no plot)
>>> plt.plot(.5*(bins[1:]+bins[:-1]), n)
>>> plt.show()更进一步阅读
- The Python tutorial
- NumPy Reference
- SciPy Tutorial
- SciPy Lecture Notes
- A matlab, R, IDL, NumPy/SciPy dictionary