第十五篇 NumPy⾼级应⽤

在这篇中,将会深⼊NumPy库的数组计算。这会包括ndarray更内部的细节,和更⾼级的数组操作和算法。

一、ndarray对象的内部机理
NumPy的ndarray提供了⼀种将同质数据块(可以是连续或跨越)解释为多维数组对象的⽅式。正如你之前所看到的那样,数据类型(dtype)决定了数据的解释⽅式,⽐如浮点数、整数、布尔值等。

ndarray如此强⼤的部分原因是所有数组对象都是数据块的⼀个跨度视图(strided view)。你可能想知道数组视图arr[::2,::-1]不复制任何数据的原因是什么。简单地说,ndarray不只是⼀块内存和⼀个dtype,它还有跨度信息,这使得数组能以各种步幅(step size)在内存中移动。更准确地讲,ndarray内部由以下内容组成:
             ⼀个指向数据(内存或内存映射⽂件中的⼀块数据)的指针。
             数据类型或dtype,描述在数组中的固定⼤⼩值的格⼦。
             ⼀个表示数组形状(shape)的元组。
             ⼀个跨度元组(stride),其中的整数指的是为了前进到当前维度下⼀个元素需要“跨过”的字节数。

图15-1简单地说明了ndarray的内部结构。

图15-1 简单地说明了ndarray的内部结构

例如,⼀个10×5的数组,其形状为(10,5):
np.ones((10, 5)).shape              # 输出为数据的形状: (10, 5)

⼀个典型的(C顺序,稍后详细讲解)3×4×5的float64(8个字节)数组,其跨度为(160,40,8) —— 知道跨度是⾮常有⽤的,通常,跨度在⼀个轴上越⼤,沿这个轴进⾏计算的开销就越⼤:
np.ones((3, 4, 5), dtype=np.float64).strides    # 输出:(160, 40, 8)

虽然NumPy⽤户很少会对数组的跨度信息感兴趣,但它们却是构建⾮复制式数组视图的重要因素。跨度甚⾄可以是负数,这样会使数组在内存中后向移动,⽐如在切⽚obj[::-1]或obj[:,::-1]中就是这样的。

1、NumPy数据类型体系
你可能偶尔需要检查数组中所包含的是否是整数、浮点数、字符串或Python对象。因为浮点数的种类很多(从float16到float128),判断dtype是否属于某个⼤类的⼯作⾮常繁琐。所幸,dtype都有⼀个超类(⽐如np.integer和np.floating),它们可以跟np.issubdtype函数结合使⽤
ints = np.ones(10, dtype=np.uint16)
floats = np.ones(10, dtype=np.float32)
np.issubdtype(ints.dtype, np.integer)           # 输出:True
np.issubdtype(floats.dtype, np.floating)       # 输出:True

调⽤dtype的mro⽅法即可查看其所有的⽗类:
np.float64.mro()        # 输出如下:
[numpy.float64,
  numpy.floating,
  numpy.inexact,
  numpy.number,
  numpy.generic,
  float,
  object]
然后得到:
np.issubdtype(ints.dtype, np.number)            # 输出:True

⼤部分NumPy⽤户完全不需要了解这些知识,但是这些知识偶尔还是能派上⽤场。图15-2说明了dtype体系以及
⽗⼦类关系。

图15-2 NumPy的dtype体系
                                       图15-2  NumPy的dtype体系

二、⾼级数组操作
除花式索引、切⽚、布尔条件取⼦集等操作之外,数组的操作⽅式还有很多。虽然pandas中的⾼级函数可以处理数据分析⼯作中的许多重型任务,但有时你还是需要编写⼀些在现有库中找不到的数据算法。

1、数组重塑
多数情况下,你可以⽆需复制任何数据,就将数组从⼀个形状转换为另⼀个形状。只需向数组的实例⽅法reshape传⼊⼀个表示新形状的元组即可实现该⽬的。例如,假设有⼀个⼀维数组,希望将其重新排列为⼀个矩阵(结果⻅图15-3):
arr = np.arange(12)
arr.reshape((4, 3))     # 输出如下:默认order参数为 'C'
array([[ 0,  1,  2],
           [ 3,  4,  5],
           [ 6,  7,  8],
           [ 9, 10, 11]])
          
arr1 = np.arange(12)
arr1.reshape((4,3), order='F')      # 输出如下:修改order参数为 'F'
array([[ 0,  4,  8],
           [ 1,  5,  9],
           [ 2,  6, 10],
           [ 3,  7, 11]])

图15-3 按C顺序(按⾏)和按Fortran顺序(按列)进⾏重塑
                        图15-3   按C顺序(按⾏)和按Fortran顺序(按列)进⾏重塑

多维数组也能被重塑
arr.reshape((4, 3)).reshape((3, 4))   # 输出如下:
array([[ 0,  1,  2,  3],
           [ 4,  5,  6,  7],
           [ 8,  9, 10, 11]])

形状参数的其中⼀维可以是 -1,它表示该维度的⼤⼩由数据本身推断⽽来
arr = np.arange(15)
arr.reshape((5, -1))    # 列参数是-1,由数据本身推断而来,输出如下:
array([[ 0,  1,  2],
           [ 3,  4,  5],
           [ 6,  7,  8],
           [ 9, 10, 11],
           [12, 13, 14]])

与reshape将⼀维数组转换为多维数组的运算过程相反的运算通常称为扁平化(flattening)或散开(raveling)
arr.ravel()   # ravel() 方法将arr阵列扁平化,输出如下:
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])

如果结果中的值与原始数组相同,ravel不会产⽣源数据的副本。flatten⽅法的⾏为类似于ravel,只不过它总是返回数据的副本:
arr1 = np.arange(15).reshape((5, 3))       # arr1是5行3列的阵列,输出省略
arr1.flatten()          # flatten()方法不会修改原始的arr1阵列,输出如下:
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])

数组可以被重塑或散开为别的顺序。下⼀⼩节中专⻔讲解这个问题。

2、C和Fortran顺序
NumPy允许你更为灵活地控制数据在内存中的布局。默认情况下,NumPy数组是按⾏优先顺序创建的。在空间⽅⾯,这就意味着,对于⼀个⼆维数组,每⾏中的数据项是被存放在相邻内存位置上的。另⼀种顺序是列优先顺序,它意味着每列中的数据项是被存放在相邻内存位置上的

由于⼀些历史原因,⾏和列优先顺序⼜分别称为C和Fortran顺序。在FORTRAN 77中,矩阵全都是列优先的。

像reshape和reval这样的函数,都可以接受⼀个表示数组数据存放顺序的order参数。⼀般可以是'C'或'F'(还有'A'和'K'等不常⽤的选项,具体请参考NumPy的⽂档)。前面的图15-3对此进⾏了说明:
arr = np.arange(12).reshape((3, 4))
arr         # 输出如下:
array([[ 0,  1,  2,  3],
           [ 4,  5,  6,  7],
           [ 8,  9, 10, 11]])
arr.ravel()     # 默认行优先,输出如下:
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
arr.ravel('F')      # 传递'F'参数,以列优先,输出如下:
array([ 0,  4,  8,  1,  5,  9,  2,  6, 10,  3,  7, 11])

⼆维或更⾼维数组的重塑过程⽐较令⼈费解(⻅图15-3)。C和Fortran顺序的关键区别就是维度的⾏进顺序
             C/⾏优先顺序:先经过更⾼的维度(例如,轴1会先于轴0被处理)。
             Fortran/列优先顺序:后经过更⾼的维度(例如,轴0会先于轴1被处理)。

3、数组的合并和拆分
numpy.concatenate可以按指定轴将⼀个由数组组成的序列(如元组、列表等)连接到⼀起:
arr1 = np.array([[1, 2, 3], [4, 5, 6]])
arr2 = np.array([[7, 8, 9], [10, 11, 12]])
np.concatenate([arr1, arr2], axis=0)            # axis=0按行连接。输出如下:
array([[ 1,  2,  3],
           [ 4,  5,  6],
           [ 7,  8,  9],
           [10, 11, 12]])
np.concatenate([arr1, arr2], axis=1)            # axis=1按列连接。输出如下:
array([[ 1,  2,  3,  7,  8,  9],
           [ 4,  5,  6, 10, 11, 12]])

对于常⻅的连接操作,NumPy提供了⼀些⽐较⽅便的⽅法(如vstack和hstack)。因此,上⾯的运算还可以表达为:
np.vstack((arr1, arr2))     # 注意参数是元组,按行连接,输出如下:
array([[ 1,  2,  3],
           [ 4,  5,  6],
           [ 7,  8,  9],
           [10, 11, 12]])
np.hstack((arr1, arr2))    # 注意参数是元组,按列连接,输出如下:
array([[ 1,  2,  3,  7,  8,  9],
           [ 4,  5,  6, 10, 11, 12]])

与此相反,split⽤于将⼀个数组沿指定轴拆分为多个数组
arr = np.random.randn(5, 2)
arr         # 输出如下:
array([[ 1.4841,  1.0304],
           [ 2.6368, -0.8897],
           [ 0.5567,  1.1119],
           [-0.1207,  0.147 ],
           [ 2.2737,  0.3816]])
first, second, third = np.split(arr, [1, 3])    # 参数[1, 3]表示在行索引为1和3处拆分,索引从0开始
first            # 输出:array([[1.4841, 1.0304]])
second     # 输出如下:
array([[ 2.6368, -0.8897],
           [ 0.5567,  1.1119]])
third       # 输出如下:
array([[-0.1207,  0.147 ],
           [ 2.2737,  0.3816]])
传⼊到np.split的值[1,3]指示在哪个索引处分割数组。

表15-1中列出了所有关于数组连接和拆分的函数,其中有些是专⻔为了⽅便常⻅的连接运算⽽提供的。
函数                                  说明
concatenate                      最一般化的连接,沿一条轴连接一组数组
vstack, row_stack              以面向行的方式对数组进行堆叠(沿轴0)
hstack                               以面向列的方式对数组进行堆叠(沿轴1)
column_stack                    类似于hstack,但是会先将一维数组转换为二维列向量
dstack                               以面向“深度”的方式对数组进行堆叠(沿轴2)
split                                   沿指定轴在指定的位置拆分数组
hsplit/vsplit/dsplit            split的便捷化函数,分别沿轴0、轴1、轴2进行拆分

4、堆叠辅助类:r_和c_
NumPy命名空间中有两个特殊的对象 r_和c_,它们可以使数组的堆叠操作更为简洁:
arr = np.arange(6)
arr1 = arr.reshape((3, 2))
arr2 = np.random.randn(3, 2)
np.r_[arr1, arr2]       # 注意参数形式,按行连接,输出如下:
array([[ 0.    ,      1.    ],
           [ 2.    ,      3.    ],
           [ 4.    ,      5.    ],
           [-0.0714,  0.8247],
           [-2.0555, -0.8911],
           [-0.1912, -1.5544]])

np.c_[np.r_[arr1, arr2], arr]       # 注意参数形式,先按行连接,后按列连接,输出如下:
array([[ 0.    ,      1.    ,       0.    ],
           [ 2.    ,      3.    ,       1.    ],
           [ 4.    ,      5.    ,       2.    ],
           [-0.0714,  0.8247,  3.    ],
           [-2.0555, -0.8911,  4.    ],
           [-0.1912, -1.5544,  5.    ]])

它还可以将切⽚转换成数组
np.c_[1:6, -10:-5]      # 输出如下:
array([[  1, -10],
           [  2,  -9],
           [  3,  -8],
           [  4,  -7],
           [  5,  -6]])
r_和c_的具体功能请参考其⽂档。

5、元素的重复操作:tile和repeat
对数组进⾏重复以产⽣更⼤数组的⼯具主要是repeat和tile这两个函数。repeat会将数组中的各个元素重复⼀定次数,从⽽产⽣⼀个更⼤的数组
arr = np.arange(3)
arr         # 输出:array([0, 1, 2])
arr.repeat(3)           # 每个元素重复3次,输出如下:
array([0, 0, 0, 1, 1, 1, 2, 2, 2])
注意:跟其他流⾏的数组编程语⾔(如MATLAB)不同,NumPy中很少需要对数组进⾏重复(replicate)。这主要是因为⼴播能更好地满⾜该需求。

默认情况下,如果传⼊的是⼀个整数,则各元素就都会重复那么多次。如果传⼊的是⼀组整数,则各元素就可以重复不同的次数:
arr.repeat([2, 3, 4])   # 参数是列表,重复的次数会不一样
array([0, 0, 1, 1, 1, 2, 2, 2, 2])

对于多维数组,还可以让它们的元素沿指定轴重复:
arr = np.random.randn(2, 2)
arr         # 输出如下:
array([[ 0.6665, -1.0611],
           [ 2.1659, -0.5874]])
arr.repeat(2, axis=0)   # 重复行,输出如下:
array([[ 0.6665, -1.0611],
           [ 0.6665, -1.0611],
           [ 2.1659, -0.5874],
           [ 2.1659, -0.5874]])

注意,如果没有设置轴向,则数组会被扁平化,这可能不会是你想要的结果。同样,在对多维进⾏重复时,也可以传⼊⼀组整数,这样就会使各切⽚重复不同的次数:
arr.repeat([2, 3], axis=0)          # 输出如下:第一行重复2次,第二行重复3次
array([[ 0.6665, -1.0611],
           [ 0.6665, -1.0611],
           [ 2.1659, -0.5874],
           [ 2.1659, -0.5874],
           [ 2.1659, -0.5874]])
arr.repeat([2, 3], axis=1)          # 输出如下:
array([[ 0.6665,  0.6665, -1.0611, -1.0611, -1.0611],
           [ 2.1659,  2.1659, -0.5874, -0.5874, -0.5874]])

tile的功能是沿指定轴向堆叠数组的副本。你可以形象地将其想象成“铺瓷砖”:
arr         # 输出如下:
array([[ 0.6665, -1.0611],
           [ 2.1659, -0.5874]])
np.tile(arr, 2)         # 输出如下:
array([[ 0.6665, -1.0611,  0.6665, -1.0611],
           [ 2.1659, -0.5874,  2.1659, -0.5874]])

第⼆个参数是瓷砖的数量。对于标量,瓷砖是⽔平铺设的,⽽不是垂直铺设。它可以是⼀个表示“铺设”布局的元组:
np.tile(arr, (2, 1))    # 行重复2次,列重复1次,输出如下:
array([[ 0.6665, -1.0611],
           [ 2.1659, -0.5874],
           [ 0.6665, -1.0611],
           [ 2.1659, -0.5874]])
np.tile(arr, (3, 2))    # 行重复3次,列重复2次,输出如下:
array([[ 0.6665, -1.0611,  0.6665, -1.0611],
           [ 2.1659, -0.5874,  2.1659, -0.5874],
           [ 0.6665, -1.0611,  0.6665, -1.0611],
           [ 2.1659, -0.5874,  2.1659, -0.5874],
           [ 0.6665, -1.0611,  0.6665, -1.0611],
           [ 2.1659, -0.5874,  2.1659, -0.5874]])

6、花式索引的等价函数:take和put
在第四篇中说过,获取和设置数组⼦集的⼀个办法是通过整数数组使⽤花式索引:
arr = np.arange(10) * 100           # 每个元素乘以100
inds = [7, 1, 2, 6]
arr[inds]   # 根据inds列表的索引选取数组的元素,输出如下:
array([700, 100, 200, 600])

ndarray还有其它⽅法⽤于获取单个轴向上的选区:
arr.take(inds)          # 使用take方法根据inds列表所指的索引选取元素,输出如下:
array([700, 100, 200, 600])
arr.put(inds, 42)       # 使用put方法根据inds列表所指的索引元素替换为42
arr         # arr数组对应inds元素索引处的值已被替换,输出如下:
array([  0,  42,  42, 300, 400, 500,  42,  42, 800, 900])
arr.put(inds, [40, 41, 42, 43])     # 使用列表一一对应替换
arr         # 输出如下:
array([  0,  41,  42, 300, 400, 500,  43,  40, 800, 900])

要在其它轴上使⽤take,只需传⼊axis关键字即可:
inds = [2, 0, 2, 1]
arr = np.random.randn(2, 4)
arr         # 输出如下:
array([[-0.1307,  0.6228, -0.9549,  0.9602],
           [ 0.1764, -1.1677, -1.1503,  0.5795]])
arr.take(inds, axis=1)  # 根据inds的索引按列选取,输出如下:
array([[-0.9549, -0.1307, -0.9549,  0.6228],
           [-1.1503,  0.1764, -1.1503, -1.1677]])

put不接受axis参数,它只会在数组的扁平化版本(⼀维,C顺序)上进⾏索引。因此,在需要⽤其他轴向的索引设置元素时,最好还是使⽤花式索引。

三、⼴播
广播(broadcasting)指的是不同形状的数组之间的算术运算的执⾏⽅式。它是⼀种⾮常强⼤的功能,但也容易令⼈误解,即使是经验丰富的⽼⼿也是如此。将标量值跟数组合并时就会发⽣最简单的⼴播:
arr = np.arange(5)
arr         # 输出:array([0, 1, 2, 3, 4])
arr * 4     # 标量乘法广播到arr数组的每个元素,输出如下:
array([ 0,  4,  8, 12, 16])
这⾥我们说:在这个乘法运算中,标量值4被⼴播到了其他所有的元素上。

看⼀个例⼦,我们可以通过减去列平均值的⽅式对数组的每⼀列进⾏距平化处理。这个问题解决起来⾮常简单:
arr = np.random.randn(4, 3)
arr.mean(0)       # 参数0指定按列求平均值,输出如下:
array([0.3326, 0.2594, 0.1337])
demeaned = arr - arr.mean(0)        # 减去平均值
demeaned    # 输出如下:
array([[-0.0583,  0.0324, -1.2212],
           [-1.4145, -0.7383,  0.4903],
           [ 0.8041, -0.1439,  0.9819],
           [ 0.6687,  0.8498, -0.251 ]])
demeaned.mean(0)        # 输出:array([0., 0., 0.])

图15-4形象地展示了该过程。⽤⼴播的⽅式对⾏进⾏距平化处理会稍微麻烦⼀些。幸运的是,只要遵循⼀定的规则,低维度的值是可以被⼴播到数组的任意维度的(⽐如对⼆维数组各列减去⾏平均值)。
图15-4  ⼀维数组在轴0上的⼴播
                                      图15-4   ⼀维数组在轴0上的⼴播

于是就得到了广播的原则:
          如果两个数组的后缘维度(trailing dimension ,即从末尾开始算起的维度)的轴长相符或其中一方的长度为1,则认为它们是广播兼容的。广播会在缺失和(或)长度为1的维度上进行。

作为⼀名经验丰富的NumPy⽼⼿,经常还是得停下来画张图并想想⼴播的原则。再来看⼀下最后那个例⼦,假设你希望对各⾏减去那个平均值。由于arr.mean(0)的⻓度为3,所以它可以在0轴向上进⾏⼴播:因为arr的后缘维度是3,所以它们是兼容的。根据该原则,要在1轴向上做减法(即各⾏减去⾏平均值),较⼩的那个数组的形状必须是(4,1):
arr         # 有arr阵列如下:
array([[ 0.2551, -1.207 , -0.5938],
           [ 0.5544, -0.6003, -0.4302],
           [-2.1395, -0.3787, -1.5833],
           [ 0.7381,  0.5709, -0.6565]])
row_means = arr.mean(1)
row_means.shape         # 其形状是一维数组,输出是:(4,)
row_means.reshape((4, 1))           # 重设形状
array([[-0.5152],
           [-0.1587],
           [-1.3672],
           [ 0.2175]])
demeaned = arr - row_means.reshape((4, 1))
demeaned.mean(1)        # 极限是0,输出如下:
array([-3.7007e-17,  0.0000e+00, -7.4015e-17, -7.4015e-17])
图15-5说明了该运算的过程。

图15-5  ⼆维数组在轴1上的⼴播
                                 图15-5   ⼆维数组在轴1上的⼴播

图15-6展示了另外⼀种情况,这次是在⼀个三维数组上沿0轴向加上⼀个⼆维数组。

图15-6  三维数组在轴0上的⼴播
                                      图15-6   三维数组在轴0上的⼴播

1、沿其它轴向⼴播
⾼维度数组的⼴播似乎更难以理解,⽽实际上它也是遵循⼴播原则的。如果不然,你就会得到下⾯这样⼀个错误:
arr - arr.mean(1)       # 输出如下错误提示:
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-118-8b8ada26fac0> in <module>()
----> 1 arr - arr.mean(1)
ValueError: operands could not be broadcast together with shapes (4,3) (4,)

⼈们经常需要通过算术运算过程将较低维度的数组在除0轴以外的其他轴向上⼴播。根据⼴播的原则,较⼩数组的“⼴播维”必须为1。在上⾯那个⾏距平化的例⼦中,这就意味着要将⾏平均值的形状变成(4,1)⽽不是(4,):
arr - arr.mean(1).reshape((4, 1))   # 输出如下:
array([[ 0.7703, -0.6917, -0.0786],
           [ 0.7131, -0.4416, -0.2715],
           [-0.7723,  0.9885, -0.2161],
           [ 0.5206,  0.3534, -0.874 ]])

对于三维的情况,在三维中的任何⼀维上⼴播其实也就是将数据重塑为兼容的形状⽽已。图15-7说明了要在三维数组各维度上⼴播的形状需求。
图15-7   能在该三维数组上⼴播的⼆维数组的形状
                                           图15-7    能在该三维数组上⼴播的⼆维数组的形状

于是就有了⼀个⾮常普遍的问题(尤其是在通⽤算法中),即专⻔为了⼴播⽽添加⼀个⻓度为1的新轴。虽然
reshape是⼀个办法,但插⼊轴需要构造⼀个表示新形状的元组。这是⼀个很郁闷的过程。因此,NumPy数组提
供了⼀种通过索引机制插⼊轴的特殊语法。下⾯这段代码通过特殊的np.newaxis属性以及“全”切⽚来插⼊新轴
arr = np.zeros((4, 4))              # 4行4列阵列
arr_3d = arr[:, np.newaxis, :]      # 4个1行4列阵列
arr_3d.shape            # 输出:(4, 1, 4)
arr_1d = np.random.normal(size=3)   # 1行3列
arr_1d[:, np.newaxis]   # 输出如下:
array([[ 1.0429],
           [-0.5402],
           [-0.5034]])
arr_1d[np.newaxis, :]   # 输出:array([[ 1.0429, -0.5402, -0.5034]])

因此,如果我们有⼀个三维数组,并希望对轴2进⾏距平化,那么只需要编写下⾯这样的代码就可以了:
arr = np.random.randn(3, 4, 5)
depth_means = arr.mean(2)
depth_means             # 输出如下:
array([[-0.1045,  0.4648,  0.2302,  0.1905],
           [-0.63  ,  0.0221,  0.5197,  0.2534],
           [ 0.9099, -0.3619,  0.231 ,  0.2467]])
depth_means.shape       # 输出:(3, 4)
demeaned = arr - depth_means[:, :, np.newaxis]
demeaned.mean(2)        # 输出如下:计算有误差,极限是0
array([[-4.4409e-17,  4.4409e-17,  0.0000e+00,  0.0000e+00],
           [-8.8818e-17,  0.0000e+00, -6.6613e-17,  4.4409e-17],
           [-3.3307e-17,  0.0000e+00,  0.0000e+00,  2.2204e-17]])

你可能会想,在对指定轴进⾏距平化时,有没有⼀种既通⽤⼜不牺牲性能的⽅法呢?实际上是有的,但需要⼀些索引⽅⾯的技巧
def demean_axis(arr, axis=0):
       means = arr.mean(axis)
      
       # This generalizes thing like [:, :, np.newaxis] to N dimensions
       indexer = [slice(None)] * arr.ndim
       indexer[axis] = np.newaxis
       return arr - means[indexer]
      
slice(None)           # 该方法产生的结果不可迭代,输出如下:
slice(None, None, None)
在上面代码中:假设arr.ndim是3,即3维数组,由于参数axis默认是0,则means[indexer]表达式就是
means[np.newaxis, :, :],若axis=1,则是means[:, np.newaxis, :]。熟悉该函数的写法。

2、通过⼴播设置数组的值
算术运算所遵循的⼴播原则同样也适⽤于通过索引机制设置数组值的操作。对于最简单的情况,我们可以这样做:
arr = np.zeros((4, 3))
arr[:] = 10             # 广播到整个数据,不能使用arr = 10
arr         # 输出如下:
array([[10., 10., 10.],
           [10., 10., 10.],
           [10., 10., 10.],
           [10., 10., 10.]])

但是,假设我们想要⽤⼀个⼀维数组来设置⽬标数组的各列,只要保证形状兼容就可以了:
col = np.array([1.28, -0.42, 0.44, 1.6])
arr[:] = col[:, np.newaxis]         # col[:, np.newaxis]方法将col转换为4行1列的数组
arr         # 输出如下:
array([[ 1.28,  1.28,  1.28],
           [-0.42, -0.42, -0.42],
           [ 0.44,  0.44,  0.44],
           [ 1.6 ,  1.6 ,  1.6 ]])
arr[:2] = [[0], [1]]    # 对切片重复赋值,注意值是列表的列表
arr          # 输出如下:
array([[0.  , 0.  , 0.  ],
           [1.  , 1.  , 1.  ],
           [0.44, 0.44, 0.44],
           [1.6 , 1.6 , 1.6 ]])

四、ufunc⾼级应⽤
虽然许多NumPy⽤户只会⽤到通⽤函数所提供的快速的元素级运算,但通⽤函数实际上还有⼀些⾼级⽤法能使我们丢开循环⽽编写出更为简洁的代码。

1、ufunc实例⽅法
NumPy的各个⼆元ufunc都有⼀些⽤于执⾏特定⽮量化运算的特殊⽅法。表15-2汇总了这些⽅法,下⾯将通过⼏个具体的例⼦对它们进⾏说明。

reduce接受⼀个数组参数,并通过⼀系列的⼆元运算对其值进⾏聚合(可指明轴向)。例如,我们可以⽤np.add.reduce对数组中各个元素进⾏求和:
arr = np.arange(10)
np.add.reduce(arr)      # 输出:45
arr.sum()               # 输出:45

起始值取决于ufunc(对于add的情况,就是0)。如果设置了轴号,约简运算就会沿该轴向执⾏。这就使你能⽤⼀种⽐较简洁的⽅式得到某些问题的答案。在下⾯这个例⼦中,我们⽤np.logical_and检查数组各⾏中的值是否是有序的:
np.random.seed(12346)   # for reproducibility
arr = np.random.randn(5, 5)
arr[::2].sort(1)        #  排序,隔2行取一行,要修改原数组
arr[:, :-1] < arr[:, 1:]            # 输出如下:注意是5行4列
array([[ True,  True,  True,  True],
           [False,  True, False, False],
           [ True,  True,  True,  True],
           [ True, False,  True,  True],
           [ True,  True,  True,  True]])
下面一行代码的输出表示:从上面的输出可知,行全为True,下面的输出才为True,否则为False
np.logical_and.reduce(arr[:, :-1] < arr[:, 1:], axis=1)     # 输出如下:
array([ True, False,  True, False,  True])
注意,logical_and.reduce跟all⽅法是等价的。

accumulate跟reduce的关系就像cumsum跟sum的关系那样。它产⽣⼀个跟原数组⼤⼩相同的中间“累计”值数组
np.add.accumulate(arr, axis=1)      # 与arr.cumsum(1)的结果是一样的,输出如下:
array([[ 0,  1,  3,  6, 10],
           [ 5, 11, 18, 26, 35],
           [10, 21, 33, 46, 60]], dtype=int32)

outer⽤于计算两个数组的叉积
arr = np.arange(3).repeat([1, 2, 2])
arr         # 输出:array([0, 1, 1, 2, 2])
np.multiply.outer(arr, np.arange(5))            # 输出如下:
array([[0, 0, 0, 0, 0],
           [0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4],
           [0, 2, 4, 6, 8],
           [0, 2, 4, 6, 8]])

outer输出结果的维度是两个输⼊数据的维度之和:
x, y = np.random.randn(3, 4), np.random.randn(5)
result = np.subtract.outer(x, y)
result.shape            # 输出:(3, 4, 5)

最后⼀个⽅法reduceat⽤于计算“局部约简”,其实就是⼀个对数据各切⽚进⾏聚合的groupby运算。它接受⼀组⽤于指示如何对值进⾏拆分和聚合的“⾯元边界”:
arr = np.arange(10)
np.add.reduceat(arr, [0, 5, 8])     # 输出:array([10, 18, 17], dtype=int32)
最终结果是在arr[0:5]、arr[5:8]以及arr[8:]上执⾏的约简。

跟其他⽅法⼀样,这⾥也可以传⼊⼀个axis参数
arr = np.multiply.outer(np.arange(4), np.arange(5))
arr          # 输出如下:
array([[ 0,  0,  0,  0,  0],
           [ 0,  1,  2,  3,  4],
           [ 0,  2,  4,  6,  8],
           [ 0,  3,  6,  9, 12]])

np.add.reduceat(arr, [0, 2, 4], axis=1)          # 输出如下:
array([[ 0,  0,  0],
           [ 1,  5,  4],
           [ 2, 10,  8],
           [ 3, 15, 12]], dtype=int32)

表15-2总结了部分的ufunc⽅法。
表15-2  部分ufunc方法
方法                                 说明
reduce(x)                         通过连续执行原始运算的方式对值进行聚合
accumulate(x)                  聚合值,保留所有局部聚合结果
reduceat(x, bins)            “局部”约简(也就是groupby)。约简数据的各个切片以产生聚合型数组
outer(x, y)                         对x和y中的每对元素应用原始运算。结果数组的形状为x.shape + y.shape

2、编写新的ufunc
有多种⽅法可以让你编写⾃⼰的NumPy ufuncs。最常⻅的是使⽤NumPy C API。在这里,只说纯粹的Python ufunc。

numpy.frompyfunc接受⼀个Python函数以及两个分别表示输⼊输出参数数量的参数。例如,下⾯是⼀个能够实现元素级加法的简单函数:
def add_elements(x, y):
       return x + y
add_them = np.frompyfunc(add_elements, 2, 1)    # 2表示需要输入2个参数,1表示只有一个输出
add_them(np.arange(8), np.arange(8))            # 输出如下:
array([0, 2, 4, 6, 8, 10, 12, 14], dtype=object)

frompyfunc创建的函数总是返回Python对象数组,这⼀点很不⽅便。幸运的是,还有另⼀个办法,即
numpy.vectorize。虽然没有frompyfunc那么强⼤,但可以让你指定输出类型:
add_them = np.vectorize(add_elements, otypes=[np.float64])
add_them(np.arange(8), np.arange(8))            # 输出如下:
array([ 0.,  2.,  4.,  6.,  8., 10., 12., 14.])

虽然这两个函数提供了⼀种创建ufunc型函数的⼿段,但它们⾮常慢,因为它们在计算每个元素时都要执⾏⼀次Python函数调⽤,这就会⽐NumPy⾃带的基于C的ufunc慢很多:
arr = np.random.randn(10000)
%timeit add_them(arr, arr)
35.2 ms ± 433 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit np.add(arr, arr)
4.76 µs ± 243 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

本篇的后⾯,会介绍使⽤Numba(http://numba.pydata.org/),创建快速Python ufuncs。

五、结构化和记录式数组
到⽬前为⽌所讨论的ndarray都是⼀种同质数据容器,也就是说,在它所表示的内存块中,各元素占⽤的字节数相同(具体根据dtype⽽定)。从表⾯上看,它似乎不能⽤于表示异质或表格型的数据。结构化数组是⼀种特殊的ndarray,其中的各个元素可以被看做C语⾔中的结构体(struct,这就是“结构化”的由来)或SQL表中带有多个命名字段的⾏
dtype = [('x', np.float64), ('y', np.int32)]    # 元组列表定义结构化dtype
sarr = np.array([(1.5, 6), (np.pi, -2)], dtype=dtype)
sarr        # 输出如下
array([(1.5   ,  6), (3.1416, -2)], dtype=[('x', '<f8'), ('y', '<i4')])

定义结构化dtype(请参考NumPy的在线⽂档)的⽅式有很多。最典型的办法是元组列表,各元组的格式为(field_name,field_data_type)。这样,数组的元素就成了元组式的对象,该对象中各个元素可以像字典那样进⾏访问
sarr[0]     # 输出:(1.5, 6)
sarr[0]['y']            # 从输出可知sarr是按列方式与dtype对齐的
输出结果是:6

字段名保存在dtype.names属性中。在访问结构化数组的某个字段时,返回的是该数据的视图,所以不会发⽣数据复制:
sarr.dtype       # 输出:dtype([('x', '<f8'), ('y', '<i4')])
sarr.dtype.name         # 输出:'void96'
sarr['x']              # 输出:array([1.5   , 3.1416])

1、嵌套dtype和多维字段
在定义结构化dtype时,你可以再设置⼀个形状(可以是⼀个整数,也可以是⼀个元组)
dtype = [('x', np.int64, 3), ('y', np.int32)]
arr = np.zeros(4, dtype=dtype)
arr         # 输出如下:理解输出的含义
array([([0, 0, 0], 0), ([0, 0, 0], 0), ([0, 0, 0], 0), ([0, 0, 0], 0)],
          dtype=[('x', '<i8', (3,)), ('y', '<i4')])

在这种情况下,各个记录的x字段所表示的是⼀个⻓度为3的数组:
arr[0]['x']             # 输出:array([0, 0, 0], dtype=int64)

这样,访问arr['x']即可得到⼀个⼆维数组,⽽不是前⾯那个例⼦中的⼀维数组:
arr['x']       # 输出如下:
array([[0, 0, 0],
           [0, 0, 0],
           [0, 0, 0],
           [0, 0, 0]], dtype=int64)

这就使你能⽤单个数组的内存块存放复杂的嵌套结构。你还可以嵌套dtype,作出更复杂的结构。下⾯是⼀个简单的例⼦:
dtype = [('x', [('a', 'f8'), ('b', 'f4')]), ('y', np.int32)]
data = np.array([((1, 2), 5), ((3, 4), 6)], dtype=dtype)
data['x']      # 输出如下:
array([(1., 2.), (3., 4.)], dtype=[('a', '<f8'), ('b', '<f4')])
data['y']       # 输出:array([5, 6])
data['x']['a']          # 输出:array([1., 3.])

pandas的DataFrame并不直接⽀持该功能,但它的分层索引机制跟这个差不多。

2、为什么要⽤结构化数组
跟pandas的DataFrame相⽐,NumPy的结构化数组是⼀种相对较低级的⼯具。它可以将单个内存块解释为带有任意复杂嵌套列的表格型结构。由于数组中的每个元素在内存中都被表示为固定的字节数,所以结构化数组能够提供⾮常快速⾼效的磁盘数据读写(包括内存映像)、⽹络传输等功能

结构化数组的另⼀个常⻅⽤法是,将数据⽂件写成定⻓记录字节流,这是C和C++代码中常⻅的数据序列化⼿段(业界许多历史系统中都能找得到)。只要知道⽂件的格式(记录的⼤⼩、元素的顺序、字节数以及数据类型等),就可以⽤np.fromfile将数据读⼊内存。这种⽤法有点超范围,知道这点就可以了。

六、更多有关排序的话题
跟Python内置的列表⼀样,ndarray的sort实例⽅法也是就地排序。也就是说,数组内容的重新排列是不会产⽣新数组的:
arr = np.random.randn(6)
arr.sort()    # 默认升序排序,从小到大
arr         # 输出如下:
array([-1.1431, -0.3266,  0.0883,  0.1259,  0.2245,  1.6527])

在对数组进⾏就地排序时要注意⼀点,如果⽬标数组只是⼀个视图,则原始数组将会被修改
arr = np.random.randn(3,5)
arr         # 输出如下:
array([[ 0.1167, -1.8784, -2.1457,  0.9693, -0.2629],
           [-2.259 , -0.2965,  0.445 ,  0.0813, -1.8176],
           [-0.066 , -1.6586,  0.056 , -0.3691, -0.8012]])
arr[:, 0].sort()        # 对数组的第一列值就地(in-place)排序
arr         # 输出如下:
array([[-2.259 , -1.8784, -2.1457,  0.9693, -0.2629],
           [-0.066 , -0.2965,  0.445 ,  0.0813, -1.8176],
           [ 0.1167, -1.6586,  0.056 , -0.3691, -0.8012]])

相反,numpy.sort会为原数组创建⼀个已排序副本。另外,它所接受的参数(⽐如kind)跟ndarray.sort⼀样:
arr = np.random.randn(5)
arr         # 输出如下:
array([ 0.505 , -0.9533,  0.2046,  1.3725, -0.477 ])
np.sort(arr)        # 输出是arr数组已排序的一个副本,输出如下:
array([-0.9533, -0.477 ,  0.2046,  0.505 ,  1.3725])
arr         # 输出如下:原数组未变
array([ 0.505 , -0.9533,  0.2046,  1.3725, -0.477 ])

这两个排序⽅法都可以接受⼀个axis参数,以便沿指定轴向对各块数据进⾏单独排序:
arr = np.random.randn(3, 5)
arr         # 输出如下:
array([[-0.3618, -1.382 ,  0.768 , -0.2806, -0.8925],
           [-0.9492, -1.2928, -0.4233,  0.6935,  1.1833],
           [ 2.1808,  0.4717, -0.2035,  0.3922,  0.528 ]])
arr.sort(axis=1)        # axis=1在列方向广播
arr         # 输出如下:
array([[-1.382 , -0.8925, -0.3618, -0.2806,  0.768 ],
           [-1.2928, -0.9492, -0.4233,  0.6935,  1.1833],
           [-0.2035,  0.3922,  0.4717,  0.528 ,  2.1808]])

要注意,这两个排序⽅法都不可以被设置为降序。其实⽆所谓,因为数组切⽚会产⽣视图(也就是说,不会产⽣副本,也不需要任何其他的计算⼯作)。Python⼀个有关列表的⼩技巧:values[::-1]可以返回⼀个反序的列表。对ndarray也是如此:
arr[:, ::-1]            # 输出如下:在列方向按从大到小排序
array([[ 0.768 , -0.2806, -0.3618, -0.8925, -1.382 ],
           [ 1.1833,  0.6935, -0.4233, -0.9492, -1.2928],
           [ 2.1808,  0.528 ,  0.4717,  0.3922, -0.2035]])

1、间接排序:argsort和lexsort
在数据分析⼯作中,常常需要根据⼀个或多个键对数据集进⾏排序。例如,⼀个有关学⽣信息的数据表可能需要以姓和名进⾏排序(先姓后名)。这就是间接排序的⼀个例⼦,在pandas的章节,有不少⾼级例⼦。给定⼀个或多个键,你就可以得到⼀个由整数组成的索引数组(我亲切地称之为索引器),其中的索引值说明了数据在新顺序下的位置。argsort和numpy.lexsort就是实现该功能的两个主要⽅法。下⾯是⼀个简单的例⼦:
values = np.array([5, 0, 1, 3, 2])
indexer = values.argsort()          # 生成一个values已排序的索引器
indexer     # 注意是索引器,输出如下:
array([1, 2, 4, 3, 0], dtype=int64)
values[indexer]         # 通过已排序的索引器获取数组
array([0, 1, 2, 3, 5])

⼀个更复杂的例⼦,下⾯这段代码根据数组的第⼀⾏对其进⾏排序:
arr = np.random.randn(3, 5)
arr[0] = values
arr         # 输出如下:
array([[ 5.        ,  0.        ,  1.        ,  3.        ,  2.         ],
           [ 0.2341,  0.5446, -1.3734, -0.1121,  0.7825],
           [-0.9514, -2.1626,  0.0207,  0.634 , -1.3081]])
arr[:, arr[0].argsort()]            # 输出如下:
array([[ 0.        ,  1.        ,  2.        ,  3.         ,  5.        ],
           [ 0.5446, -1.3734,  0.7825, -0.1121,  0.2341],
           [-2.1626,  0.0207, -1.3081,  0.634 , -0.9514]])

lexsort跟argsort差不多,只不过它可以⼀次性对多个键数组执⾏间接排序(字典序)。假设我们想对⼀些以姓和名标识的数据进⾏排序:
first_name = np.array(['Bob', 'Jane', 'Steve', 'Bill', 'Barbara'])
last_name = np.array(['Jones', 'Arnold', 'Arnold', 'Jones', 'Walters'])
sorter = np.lexsort((first_name,last_name))
sorter      # 输出:array([1, 2, 3, 0, 4], dtype=int64)
zip(last_name[sorter], first_name[sorter])      # 输出是一个可迭代对象,输出如下:
<zip at 0x1ff2b8dba48>
for i,j in zip(last_name[sorter], first_name[sorter]):
       print((i,j))
上面的for循环输出如下:
('Arnold', 'Jane')
('Arnold', 'Steve')
('Jones', 'Bill')
('Jones', 'Bob')
('Walters', 'Barbara')

刚开始使⽤lexsort的时候可能会⽐较容易头晕,这是因为键的应⽤顺序是从最后⼀个传⼊的算起的。不难看出,last_name是先于first_name被应⽤的。

注意:Series和DataFrame的sort_index以及Series的order⽅法就是通过这些函数的变体(它们还必须考虑缺失值)实现的。

2、其他排序算法
稳定的(stable)排序算法会保持等价元素的相对位置。对于相对位置具有实际意义的那些间接排序⽽⾔,这⼀点⾮常重要:
values = np.array(['2:first', '2:second', '1:first', '1:second', '1:third'])
key = np.array([2, 2, 1, 1, 1])
indexer = key.argsort(kind='mergesort')
indexer     # 这里[2, 2]对应排序是[0, 1],输出如下:
array([2, 3, 4, 0, 1], dtype=int64)
values.take(indexer)    # 输出如下:注意理解这种关系
array(['1:first', '1:second', '1:third', '2:first', '2:second'],
          dtype='<U8')

mergesort(合并排序)是唯⼀的稳定排序,它保证有O(n log n)的性能(空间复杂度),但是其平均性能⽐默认的quicksort(快速排序)要差。表15-3列出了可⽤的排序算法及其相关的性能指标。⼤部分⽤户完全不需要知道这些东⻄,但了解⼀下总是好的。
表15-3  数组排序算法
kind                    速度       稳定性    工作空间        最坏的情况
'quicksort'           1           No          0                   o(n^2)
'mergesort'         2           Yes         n/2                o(n log n)
'heapsort'           3            No         0                   o(n log n)

3、部分排序数组
排序的⽬的之⼀可能是确定数组中最⼤或最⼩的元素。NumPy有两个优化⽅法,numpy.partition和np.argpartition,可以在第k个最⼩元素划分的数组:
np.random.seed(12345)
arr = np.random.randn(20)
arr         # 输出如下:
array([-0.2047,  0.4789, -0.5194, -0.5557,  1.9658,  1.3934,  0.0929,
            0.2817,  0.769 ,  1.2464,  1.0072, -1.2962,  0.275 ,  0.2289,
            1.3529,  0.8864, -2.0016, -0.3718,  1.669 , -0.4386])
  np.partition(arr, 3)   # 输出如下:头三个元素是最⼩的三个
array([-2.0016, -1.2962, -0.5557, -0.5194, -0.3718, -0.4386, -0.2047,
             0.2817,  0.769 ,  0.4789,  1.0072,  0.0929,  0.275 ,  0.2289,
            1.3529,  0.8864,  1.3934,  1.9658,  1.669 ,  1.2464])

当你调⽤partition(arr, 3),结果中的头三个元素是最⼩的三个,没有特定的顺序。numpy.argpartition与numpy.argsort相似,会返回索引,重排数据为等价的顺序:
indices = np.argpartition(arr, 3)
indices     # 输出如下:
array([16, 11,  3,  2, 17, 19,  0,  7,  8,  1, 10,  6, 12, 13, 14, 15,  5,
             4, 18,  9], dtype=int64)
arr.take(indices)       # 输出如下:
array([-2.0016, -1.2962, -0.5557, -0.5194, -0.3718, -0.4386, -0.2047,
             0.2817,  0.769 ,  0.4789,  1.0072,  0.0929,  0.275 ,  0.2289,
             1.3529,  0.8864,  1.3934,  1.9658,  1.669 ,  1.2464])

4、numpy.searchsorted:在有序数组中查找元素
searchsorted是⼀个在有序数组上执⾏⼆分查找的数组⽅法,只要将值插⼊到它返回的那个位置就能维持数组的有序性
arr = np.array([0, 1, 7, 12, 15])
arr.searchsorted(9)     # 表示9在数组arr中索引位置,输出:3

你可以传⼊⼀组值就能得到⼀组索引:
arr.searchsorted([0, 8, 11, 16])    # 输出如下:
array([0, 3, 3, 5], dtype=int64)

从上⾯的结果中可以看出,对于元素0,searchsorted会返回0。这是因为其默认⾏为是返回相等值组的左侧索引
arr = np.array([0, 0, 0, 1, 1, 1, 1])
arr.searchsorted([0, 1])             # 输出如下:
array([0, 3], dtype=int64)
arr.searchsorted([0, 1], side='right')           # 输出如下:
array([3, 7], dtype=int64)

再来看searchsorted的另⼀个⽤法,假设我们有⼀个数据数组(其中的值在0到10000之间),还有⼀个表示“⾯元边界”的数组,我们希望⽤它将数据数组拆分开:
data = np.floor(np.random.uniform(0, 10000, size=50))
bins = np.array([0, 100, 1000, 5000, 10000])
data        # 输出如下:
array([9940., 6768., 7908., 1709.,  268., 8003., 9037.,  246., 4917.,
           5262., 5963.,  519., 8950., 7282., 8183., 5002., 8101.,  959.,
           2189., 2587., 4681., 4593., 7095., 1780., 5314., 1677., 7688.,
           9281., 6094., 1501., 4896., 3773., 8486., 9110., 3838., 3154.,
           5683., 1878., 1258., 6875., 7996., 5735., 9732., 6340., 8884.,
           4954., 3516., 7142., 5039., 2256.])

然后,为了得到各数据点所属区间的编号(其中1表示⾯元[0,100)),我们可以直接使⽤searchsorted:
labels = bins.searchsorted(data)
labels      # 输出如下:
array([4, 4, 4, 3, 2, 4, 4, 2, 3, 4, 4, 2, 4, 4, 4, 4, 4, 2, 3, 3, 3, 3,
           4, 3, 4, 3, 4, 4, 4, 3, 3, 3, 4, 4, 3, 3, 4, 3, 3, 4, 4, 4, 4, 4,
           4, 3, 3, 4, 4, 3], dtype=int64)

通过pandas的groupby使⽤该结果即可⾮常轻松地对原数据集进⾏拆分:
pd.Series(data).groupby(labels).mean()          # 输出如下:
2     498.000000
3    3064.277778
4    7389.035714
dtype: float64

七、⽤Numba编写快速NumPy函数
Numba是⼀个开源项⽬,它可以利⽤CPUs、GPUs或其它硬件为类似NumPy的数据创建快速函数。它使⽤了LLVM项⽬(http://llvm.org/),将Python代码转换为机器代码。

为了介绍Numba,来考虑⼀个纯粹的Python函数,它使⽤for循环计算表达式(x - y).mean():
def mean_distance(x, y):
       nx = len(x)
       result = 0.0
       count = 0
       for i in range(nx):
           result += x[i] - y[i]
           count += 1
       return result / count
这个函数很慢:
x = np.random.randn(10000000)
y = np.random.randn(10000000)
%timeit mean_distance(x, y)         # 输出如下:
4.72 s ± 128 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit (x - y).mean()              # 输出如下:
70.7 ms ± 762 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

NumPy的版本要⽐它快过100倍。我们可以转换这个函数为编译的Numba函数,使⽤numba.jit函数
import numba as nb
numba_mean_distance = nb.jit(mean_distance)
也可以写成装饰器:
@nb.jit
def mean_distance(x, y):
       nx = len(x)
       result = 0.0
       count = 0
       for i in range(nx):
           result += x[i] - y[i]
           count += 1
       return result / count

它要⽐⽮量化的NumPy快:
%timeit numba_mean_distance(x, y)   # 输出如下:
12.6 ms ± 102 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Numba不能编译Python代码,但它⽀持纯Python写的⼀个部分,可以编写数值算法。

Numba是⼀个深厚的库,⽀持多种硬件、编译模式和⽤户插件。它还可以编译NumPy Python API的⼀部分,⽽不⽤for循环。Numba也可以识别可以编译为机器编码的结构体,但是若调⽤CPython API,它就不知道如何编译。Numba的jit函数有⼀个选项,nopython=True,它限制了可以被转换为Python代码的代码,这些代码可以编译为LLVM,但没有任何Python C API调⽤。jit(nopython=True)有⼀个简短的别名numba.njit

前⾯的例⼦,我们还可以这样写:
from numba import float64, njit
@njit(float64(float64[:], float64[:]))
  def mean_distance(x, y):
       return (x - y).mean()

建议学习Numba的线上⽂档(http://numba.pydata.org/)。下⼀节介绍⼀个创建⾃定义Numpy ufunc对象的例⼦。

1、⽤Numba创建⾃定义numpy.ufunc对象
numba.vectorize创建了⼀个编译的NumPy ufunc,它与内置的ufunc很像。考虑⼀个numpy.add的Python例⼦:
from numba import vectorize
@vectorize
def nb_add(x, y):
       return x + y
现在有:
x = np.arange(10)
nb_add(x, x)            # 输出如下:
array([ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18], dtype=int64)
nb_add.accumulate(x, axis=0)        # 输出如下:x数组的累加
array([ 0.,  1.,  3.,  6., 10., 15., 21., 28., 36., 45.])

八、⾼级数组输⼊输出
在NumPy基础中讲过,np.save和np.load可⽤于读写磁盘上以⼆进制格式存储的数组。其实还有⼀些⼯具可⽤于更为复杂的场景。尤其是内存映像(memory map),它使你能处理在内存中放不下的数据集。

1、内存映像⽂件
内存映像⽂件是⼀种将磁盘上的⾮常⼤的⼆进制数据⽂件当做内存中的数组进⾏处理的⽅式。NumPy实现了⼀个类似于ndarray的memmap对象,它允许将⼤⽂件分成⼩段进⾏读写,⽽不是⼀次性将整个数组读⼊内存。另外,memmap也拥有跟普通数组⼀样的⽅法,因此,基本上只要是能⽤于ndarray的算法就也能⽤于memmap。

要创建⼀个内存映像,可以使⽤函数np.memmap并传⼊⼀个⽂件路径、数据类型、形状以及⽂件模式
mmap = np.memmap('datasets/mymap', dtype='float64', mode='w+',
                                       shape=(100, 100))
mmap        # 输出如下:
memmap([[0., 0., 0., ..., 0., 0., 0.],
                   [0., 0., 0., ..., 0., 0., 0.],
                   [0., 0., 0., ..., 0., 0., 0.],
                         ...,
                   [0., 0., 0., ..., 0., 0., 0.],
                   [0., 0., 0., ..., 0., 0., 0.],
                   [0., 0., 0., ..., 0., 0., 0.]])

对memmap切⽚将会返回磁盘上的数据的视图:
section = mmap[:5]

如果将数据赋值给这些视图:数据会先被缓存在内存中(就像是Python的⽂件对象),调⽤flush即可将其写⼊磁盘
section[:] = np.random.randn(5, 100)
mmap.flush()            # 写入磁盘
mmap        # 输出如下:
memmap([[ 1.3714,  0.9313,  0.6057, ..., -0.0931,  0.3793, -0.0507],
                   [-0.9518,  0.3284, -1.9309, ..., -0.6463,  1.2956,  0.6098],
                   [ 0.4409, -1.678  ,  1.8496, ..., -0.6181, -0.0625,  1.0209],
                                                 ...,
                   [ 0.    ,  0.    ,  0.    , ...,  0.    ,  0.    ,  0.    ],
                   [ 0.    ,  0.    ,  0.    , ...,  0.    ,  0.    ,  0.    ],
                   [ 0.    ,  0.    ,  0.    , ...,  0.    ,  0.    ,  0.    ]])
del mmap    # 删除内存中mmap对象

只要某个内存映像超出了作⽤域,它就会被垃圾回收器回收,之前对其所做的任何修改都会被写⼊磁盘。当打开⼀
个已经存在的内存映像时,仍然需要指明数据类型和形状,因为磁盘上的那个⽂件只是⼀块⼆进制数据⽽已,没有
任何元数据:
mmap = np.memmap('datasets/mymap', dtype='float64', shape=(100, 100))
mmap        # 输出如下:
memmap([[ 1.3714,  0.9313,  0.6057, ..., -0.0931,  0.3793, -0.0507],
                   [-0.9518,  0.3284, -1.9309, ..., -0.6463,  1.2956,  0.6098],
                   [ 0.4409, -1.678  ,  1.8496, ..., -0.6181, -0.0625,  1.0209],
                                                 ...,
                   [ 0.    ,  0.    ,  0.    , ...,  0.    ,  0.    ,  0.    ],
                   [ 0.    ,  0.    ,  0.    , ...,  0.    ,  0.    ,  0.    ],
                   [ 0.    ,  0.    ,  0.    , ...,  0.    ,  0.    ,  0.    ]])

内存映像可以使⽤前⾯介绍的结构化或嵌套dtype。

2、HDF5及其他数组存储⽅式
PyTables和h5py这两个Python项⽬可以将NumPy的数组数据存储为⾼效且可压缩的HDF5格式(HDF意思是“层次化数据格式”)。你可以安全地将好⼏百GB甚⾄TB的数据存储为HDF5格式。要学习Python使⽤HDF5,请参考pandas线上⽂档。

九、性能建议
使⽤NumPy的代码的性能⼀般都很不错,因为数组运算⼀般都⽐纯Python循环快得多。下⾯⼤致列出了⼀些需要
注意的事项:
             将Python循环和条件逻辑转换为数组运算和布尔数组运算。
             尽量使⽤⼴播。
             避免复制数据,尽量使⽤数组视图(即切⽚)。
             利⽤ufunc及其各种⽅法。

如果单⽤NumPy⽆论如何都达不到所需的性能指标,可以考虑⼀下⽤C、Fortran或Cython(等下会稍微介绍⼀下)来编写代码。在⼯作中经常会⽤到Cython(http://cython.org),因为它不⽤花费太多精⼒就能得到C语⾔那样的性能。

1、连续内存的重要性
这个话题有点超范围,但还是要提⼀下,因为在某些应⽤场景中,数组的内存布局可以对计算速度造成极⼤的影响。这是因为性能差别在⼀定程度上跟CPU的⾼速缓存(cache)体系有关。运算过程中访问连续内存块(例如,对以C顺序存储的数组的⾏求和)⼀般是最快的,因为内存⼦系统会将适当的内存块缓存到超⾼速的L1或L2CPU Cache中。此外,NumPy的C语⾔基础代码(某些)对连续存储的情况进⾏了优化处理,这样就能避免⼀些跨越式的内存访问。

⼀个数组的内存布局是连续的,就是说元素是以它们在数组中出现的顺序(即Fortran型(列优先)或C型(⾏优先))存储在内存中的。默认情况下,NumPy数组是以C型连续的⽅式创建的。列优先的数组(⽐如C型连续数组的转置)也被称为Fortran型连续。通过ndarray的flags属性即可查看这些信息:
arr_c = np.ones((100, 100), order='c')
arr_f = np.ones((100, 100), order='F')
arr_c.flags             # 输出如下:
   C_CONTIGUOUS : True
   F_CONTIGUOUS : False
   OWNDATA : True
   WRITEABLE : True
   ALIGNED : True
   WRITEBACKIFCOPY : False
   UPDATEIFCOPY : False

arr_f.flags             # 输出如下:
   C_CONTIGUOUS : False
   F_CONTIGUOUS : True
   OWNDATA : True
   WRITEABLE : True
   ALIGNED : True
   WRITEBACKIFCOPY : False
   UPDATEIFCOPY : False

arr_f.flags.f_contiguous            # 输出:True

在这个例⼦中,对两个数组的⾏进⾏求和计算,理论上说,arr_c会⽐arr_f快,因为arr_c的⾏在内存中是连续的。我们可以在IPython中⽤%timeit来确认⼀下:
%timeit arr_c.sum(1)    # 运算速度够慢的,输出如下:
8.93 µs ± 43.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit arr_f.sum(1)    # 运算速度够慢的,输出如下:
10.3 µs ± 130 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

如果想从NumPy中提升性能,这⾥就应该是下⼿的地⽅。如果数组的内存顺序不符合要求使⽤copy并传⼊'C'或'F'即可解决该问题:
arr_f.copy('C').flags      # 输出如下:
   C_CONTIGUOUS : True
   F_CONTIGUOUS : False
   OWNDATA : True
   WRITEABLE : True
   ALIGNED : True
   WRITEBACKIFCOPY : False
   UPDATEIFCOPY : False

注意,在构造数组的视图时,其结果不⼀定是连续的
arr_c[:50].flags.contiguous         # 输出:True

arr_c[:, :50].flags     # 输出如下:
   C_CONTIGUOUS : False
   F_CONTIGUOUS : False
   OWNDATA : False
   WRITEABLE : True
   ALIGNED : True
   WRITEBACKIFCOPY : False
   UPDATEIFCOPY : False

原文地址:https://www.cnblogs.com/Micro0623/p/10266598.html