其他基本用法
广播(Broadcasting)
广播规则解释NumPy如何处理不同维度数组之间的算数运算问题。在一定条件下,维度较小的数组广播为更大的数组,使他们具有相同的维度形状。广播提供了一种向量化数组操作方法,使得循环发生在C而不是Python。广播不会产生没必要的数据复制,通常会使算法运行更加有效。缺点是,广播会使内存的使用效率降低,并影响计算速度。
NumPy运算执行在数组的每个元素基础上。如下所示,两个数组必须有相同的形状。
>>> a = np.array([1.0, 2.0, 3.0])
>>> b = np.array([2.0, 2.0, 2.0])
>>> a * b
array([ 2., 4., 6.])
当数组满足一定的约束时,NumPy的广播规则会降低形状必须相同的约束。广播的一个简单例子就是一个数组与一个数值进行运算。
>>> a = np.array([1.0, 2.0, 3.0])
>>> b = 2.0
>>> a * b
array([ 2., 4., 6.])
结果与第一个例子b为数组时相同。我们可以认为数值b在与数组a进行算数运算时被拉伸为与数组a相同的形状。数组b中的元素是复制数值b所得。我们通过拉伸来做概念上的类比,NumPy更加有效地利用原始数值而不是复制,使广播尽可能保证内存和计算的效率。
通用广播规则
NumPy在数组运算时会自动对比两个数组的形状。首先对比数组最后一轴的大小,依次向前对比。以下两种情况是兼容可运算的:
- 形状相同
- 其中一个维度为1
如果不符合运算条件,则会抛出ValueError: frames are not aligned
异常,表明两个数组的形状不同不能运算。运算结果的数组形状每一维都是输入数组中的最大值。
数组不需要维度数相同。例如,你有一个256x256x3
的三维RGB(颜色)值,你想将图像中的每维值用不同的数值等比例扩展。你可以用一个包含3个值得一维数组乘以图像数组。根据广播规则,对准图像数组最后一轴的大小,显示两者之间是可运算的:
Image (3d array): 256 x 256 x 3
Scale (1d array): 3
Result (3d array): 256 x 256 x 3
在两个数组的两个维度比较时,任何一个为1,则另一个即可计算。可以认为,维度为1即可以被拉伸复制去对应另一个数组的这一维。
A (4d array): 8 x 1 x 6 x 1
B (3d array): 7 x 1 x 5
Result (4d array): 8 x 7 x 6 x 5
更多例子如下所示:
A (2d array): 5 x 4
B (1d array): 1
Result (2d array): 5 x 4
A (2d array): 5 x 4
B (1d array): 4
Result (2d array): 5 x 4
A (3d array): 15 x 3 x 5
B (3d array): 15 x 1 x 5
Result (3d array): 15 x 3 x 5
A (3d array): 15 x 3 x 5
B (2d array): 3 x 5
Result (3d array): 15 x 3 x 5
A (3d array): 15 x 3 x 5
B (2d array): 3 x 1
Result (3d array): 15 x 3 x 5
不能广播的例子:
A (1d array): 3
B (1d array): 4 # 最后一维的元素个数不匹配
A (2d array): 2 x 1
B (3d array): 8 x 4 x 3 # 倒数第二维不匹配
实际应用的例子:
>>> x = np.arange(4)
>>> xx = x.reshape(4,1)
>>> y = np.ones(5)
>>> z = np.ones((3,4))
>>> x.shape
(4,)
>>> y.shape
(5,)
>>> x + y
<type 'exceptions.ValueError'>: shape mismatch: objects cannot be broadcast to a single shape
>>> xx.shape
(4, 1)
>>> y.shape
(5,)
>>> (xx + y).shape
(4, 5)
>>> xx + y
array([[ 1., 1., 1., 1., 1.],
[ 2., 2., 2., 2., 2.],
[ 3., 3., 3., 3., 3.],
[ 4., 4., 4., 4., 4.]])
>>> x.shape
(4,)
>>> z.shape
(3, 4)
>>> (x + z).shape
(3, 4)
>>> x + z
array([[ 1., 2., 3., 4.],
[ 1., 2., 3., 4.],
[ 1., 2., 3., 4.]])
广播还提供了便利的计算两个数组的外积的方法。如下所示:
>>> a = np.array([0.0, 10.0, 20.0, 30.0])
>>> b = np.array([1.0, 2.0, 3.0])
>>> a[:, np.newaxis]
array[[ 0.]
[ 10.]
[ 20.]
[ 30.]]
>>> a[:, np.newaxis] + b # 将外积转化为求和操作
array([[ 1., 2., 3.],
[ 11., 12., 13.],
[ 21., 22., 23.],
[ 31., 32., 33.]])
其中newaxis
方法想数组a插入一个新的轴,使其成为一个两维的4×1数组。使其和shape为(3,)的数组b相加,得到一个4×3的数组。
索引使用技巧
NumPy提供了比Python更多的索引功能和技巧。除了使用整数和切片索引,数组可以由整数数组和布尔型数组构建索引。
整数数组构建索引
>>> a = np.arange(12)**2 # 数组a包含前12个数的平方数
>>> i = np.array( [ 1,1,3,8,5 ] ) # 数组的索引i
>>> a[i] # 数组a中i位置的元素,数组的索引i中的元素为索引值
array([ 1, 1, 9, 64, 25])
>>>
>>> j = np.array( [ [ 3, 4], [ 9, 7 ] ] ) # 二维数组的索引j
>>> a[j] # 返回结果与数组的索引j形状相同
array([[ 9, 16],
[81, 49]])
当索引数组a
是多维时,一个单一的数组的索引作用于a的第一维。如下所示,使用调色板将标签图像转换为彩色图像。
>>> palette = np.array( [ [0,0,0], # black
... [255,0,0], # red
... [0,255,0], # green
... [0,0,255], # blue
... [255,255,255] ] ) # white
>>> 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], # 数组a第一维的索引
... [1,2] ] )
>>> j = np.array( [ [2,1], # 数组a第二位的索引
... [3,3] ] )
>>>
>>> a[i,j] # 数组i、j的形状必须相同
array([[ 2, 5],
[ 7, 11]])
>>>
>>> a[i,2]
array([[ 2, 6],
[ 6, 10]])
>>>
>>> a[:,j]
array([[[ 2, 1],
[ 3, 3]],
[[ 6, 5],
[ 7, 7]],
[[10, 9],
[11, 11]]])
也可以把数组i、j放进一个列表中,然后用这个列表做索引。
>>> l = [i,j]
>>> a[l] # 等于a[i,j]
array([[ 2, 5],
[ 7, 11]])
但是,不能把数组i、j放进另一个数组,因为新生成的数组只会对索引数组的第一维进行构建索引。
>>> s = np.array( [i,j] )
>>> a[s] # 错误,不是我们要的
Traceback (most recent call last):
File "<stdin>", line 1, in ?
IndexError: index (3) out of range (0<=index<=2) in dimension 0
>>>
>>> a[tuple(s)] # 等于 a[i,j]
array([[ 2, 5],
[ 7, 11]])
利用数组构建索引的其他常见用法还有搜索时间依存序列的最大值。
>>>
>>> time = np.linspace(20, 145, 5) # 时标
>>> data = np.sin(np.arange(20)).reshape(5,4) # 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) # 每一系列最大值的索引,axis为0表示纵轴方向,为1表示横轴方向
>>> ind
array([2, 0, 3, 1])
>>>
>>> time_max = time[ ind] # 最大值对应的时间
>>>
>>> data_max = data[ind, xrange(data.shape[1])] # data.shape为(5,4),data.shape[1]为4,xrange(4)返回0-3的列表
>>>
>>> time_max
array([ 82.5 , 20. , 113.75, 51.25])
>>> 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在索引列表里出现两次,数组a中下标为0的元素只自加了一次。因为Python需要a+=1
等于a=a+1
。(不太懂)
布尔数组构建索引
当使用整数数组构建索引时,我们是提供了索引下标列表去获取数组中的数据。使用布尔数组构建索引则不同,目的是希望通过布尔数组索引明确指出我们需要数组中的哪些元素,不需要哪些。
最自然地理解布尔检索的方法是使用相同形状的布尔数组作为原始数组。
>>> a = np.arange(12).reshape(3,4)
>>> b = a > 4
>>> b # b是和数组a形状相同的布尔数组
array([[False, False, False, False],
[False, True, True, True],
[ True, True, True, True]], dtype=bool)
>>> a[b] # 生成只包含True值对应元素的一维数组
array([ 5, 6, 7, 8, 9, 10, 11])
这个属性对应赋值很有用:
>>> a[b] = 0 # 将所有大于4的值修改为0,说明a[b]是数组a的视图
>>> a
array([[0, 1, 2, 3],
[4, 0, 0, 0],
[0, 0, 0, 0]])
举个例子演示如何使用布尔检索生成一个Mandelbrot set(曼德勃罗集合)的图像。
import numpy as np
import matplotlib.pyplot as plt
def mandelbrot(h, w, maxit=20):
# 返回一个(h, w)大小的曼德尔布罗特图形
# ogrid方法用于构建一个多维的格点矩阵
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 # 分叉处
div_now = diverge & (divtime==maxit) # 当前分叉点
divtime[div_now] = i # 记录
z[diverge] = 2 # 避免分叉太频繁
return divtime
plt.imshow(mandelbrot(400, 400)) # imshow函数是有样板文件自动生成的绘图功能
plt.show()
布尔索引检索的第二种方法和整数索引相似。对于数组的每一维,使用一维的布尔数组获取我们想要的切片。
>>> a = np.arange(12).reshape(3,4)
>>> b1 = np.array([False,True,True]) # 第一维筛选
>>> b2 = np.array([True,False,True,False]) # 第二维筛选
>>>
>>> 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] # 奇怪的地方,不解为什么结果是这个,而不是[[4, 6],[8, 10]]
array([ 4, 10])
需要注意的是一维布尔数组的长度必须和你想获取切片的维度或轴的长度相同。在上个例子中,b1是长度为3的一维数组(和数组a的行数相同),b2是长度为4的一维数组(与数组a第二维的长度相等)。
ix_()函数
可以理解为对多个不同向量中的元素进行运算,通过ix函数修饰过的向量可以直接运算。例如,如果你想向量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[3,2,4]
17
>>> a[3]+b[2]*c[4]
17
也可以使用函数减少上述代码:
>>> 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的优势在于,使用广播避免产生一个输出大小乘以向量的大小的参数数组。
线性代数(Linear Algebra)
简单的数组操作
可以在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) # 2x2 单位矩阵; "eye" 代表 "I"
>>> u
array([[ 1., 0.],
[ 0., 1.]])
>>> j = np.array([[0.0, -1.0], [1.0, 0.0]])
>>> np.dot (j, j) # 矩阵积
array([[-1., 0.],
[ 0., -1.]])
>>> np.trace(u) # 迹,矩阵对角线之和
2.0
>>> y = np.array([[5.], [7.]])
>>> np.linalg.solve(a, y) # 求解线性方程的值,```1*x1+2*x2=5```,```3*x1+4*x2=7```,数组a为方程的参数
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]]))
技巧提示
自动变形(“Automatic” Reshaping)
想要改变数组的维度,可以忽略某个维度值,然后这个值会被自动推导出来。
>>> a = np.arange(30)
>>> a.shape = 2,-1,3 # -1 表示这个值需要被推导出来
>>> 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]]])
向量堆积(Vector Stacking)
如何从一个相同大小的行向量列表构建一个二维数组?在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])
更多关于MATLAB在NumPy中的实现可以参考Numpy for Matlab users
直方图(Histograms)
histogram 函数作用于一个数组,返回两个向量:数组的直方图和箱柜向量。注意:matplotlib
包中也有个构建直方图的函数hist
,与这里的histogram
不同。区别在于pylab.hist
函数自动绘制直方图,而numpy.histogram
仅生成数据。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> # 构建一个10000的向量,正态分布方差 0.5^2 ,平均值为 2
>>> mu, sigma = 2, 0.5
>>> v = np.random.normal(mu,sigma,10000)
>>> # 绘制一个有50个箱柜的正态分布图
>>> plt.hist(v, bins=50, normed=1)
>>> plt.show()
>>> # 计算数组的直方图,并绘制出来
>>> (n, bins) = np.histogram(v, bins=50, normed=True)
>>> plt.plot(.5*(bins[1:]+bins[:-1]), n)
>>> plt.show()
延伸阅读
Python3官方教程
NumPy参考目录
SciPy 教程
scipy讲义
A matlab, R, IDL, NumPy/SciPy dictionary