第十五篇 NumPy⾼級應⽤

在這篇中,將會深⼊NumPy庫的數組計算。這會包括ndarray更內部的細節,和更⾼級的數組操做和算法。python

1、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]中就是這樣的。
dom

一、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)       # 輸出:Trueide

調⽤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體系以及
⽗⼦類關係。oop

圖15-2 NumPy的dtype體系
                                       圖15-2  NumPy的dtype體系

2、⾼級數組操做
除花式索引、切⽚、布爾條件取⼦集等操做以外,數組的操做⽅式還有不少。雖然pandas中的⾼級函數能夠處理數據分析⼯做中的許多重型任務,但有時你仍是須要編寫⼀些在現有庫中找不到的數據算法。

一、數組重塑
多數狀況下,你能夠⽆需複製任何數據,就將數組從⼀個形狀轉換爲另⼀個形狀。只需向數組的實例⽅法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])

數組能夠被重塑或散開爲別的順序。下⼀⼩節中專⻔講解這個問題。

二、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被處理)。

三、數組的合併和拆分
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、軸一、軸2進行拆分

四、堆疊輔助類: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_的具體功能請參考其⽂檔。

五、元素的重複操做: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]])

六、花式索引的等價函數: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順序)上進⾏索引。所以,在須要⽤其餘軸向的索引設置元素時,最好仍是使⽤花式索引。

3、⼴播
廣播(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上的⼴播

一、沿其它軸向⼴播
⾼維度數組的⼴播彷佛更難以理解,⽽實際上它也是遵循⼴播原則的。若是否則,你就會獲得下⾯這樣⼀個錯誤:
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, :]。熟悉該函數的寫法。

二、經過⼴播設置數組的值
算術運算所遵循的⼴播原則一樣也適⽤於經過索引機制設置數組值的操做。對於最簡單的狀況,咱們能夠這樣作:
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 ]])

4、ufunc⾼級應⽤
雖然許多NumPy⽤戶只會⽤到通⽤函數所提供的快速的元素級運算,但通⽤函數實際上還有⼀些⾼級⽤法能使咱們丟開循環⽽編寫出更爲簡潔的代碼。

一、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

二、編寫新的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。

5、結構化和記錄式數組
到⽬前爲⽌所討論的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])

一、嵌套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並不直接⽀持該功能,但它的分層索引機制跟這個差很少。

二、爲何要⽤結構化數組
跟pandas的DataFrame相⽐,NumPy的結構化數組是⼀種相對較低級的⼯具。它能夠將單個內存塊解釋爲帶有任意複雜嵌套列的表格型結構。因爲數組中的每一個元素在內存中都被表示爲固定的字節數,因此結構化數組可以提供⾮常快速⾼效的磁盤數據讀寫(包括內存映像)、⽹絡傳輸等功能。

結構化數組的另⼀個常⻅⽤法是,將數據⽂件寫成定⻓記錄字節流,這是C和C++代碼中常⻅的數據序列化⼿段(業界許多歷史系統中都能找獲得)。只要知道⽂件的格式(記錄的⼤⼩、元素的順序、字節數以及數據類型等),就能夠⽤np.fromfile將數據讀⼊內存。這種⽤法有點超範圍,知道這點就能夠了。

6、更多有關排序的話題
跟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]])

一、間接排序: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⽅法就是經過這些函數的變體(它們還必須考慮缺失值)實現的。

二、其餘排序算法
穩定的(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)

三、部分排序數組
排序的⽬的之⼀多是肯定數組中最⼤或最⼩的元素。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])

四、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

7、⽤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對象的例⼦。

一、⽤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.])

8、⾼級數組輸⼊輸出
在NumPy基礎中講過,np.save和np.load可⽤於讀寫磁盤上以⼆進制格式存儲的數組。其實還有⼀些⼯具可⽤於更爲複雜的場景。尤爲是內存映像(memory map),它使你能處理在內存中放不下的數據集。

一、內存映像⽂件
內存映像⽂件是⼀種將磁盤上的⾮常⼤的⼆進制數據⽂件當作內存中的數組進⾏處理的⽅式。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。

二、HDF5及其餘數組存儲⽅式
PyTables和h5py這兩個Python項⽬能夠將NumPy的數組數據存儲爲⾼效且可壓縮的HDF5格式(HDF意思是「層次化數據格式」)。你能夠安全地將好⼏百GB甚⾄TB的數據存儲爲HDF5格式。要學習Python使⽤HDF5,請參考pandas線上⽂檔。

9、性能建議
使⽤NumPy的代碼的性能⼀般都很不錯,由於數組運算⼀般都⽐純Python循環快得多。下⾯⼤致列出了⼀些須要
注意的事項:
             將Python循環和條件邏輯轉換爲數組運算和布爾數組運算。
             儘可能使⽤⼴播。
             避免複製數據,儘可能使⽤數組視圖(即切⽚)。
             利⽤ufunc及其各類⽅法。

若是單⽤NumPy⽆論如何都達不到所需的性能指標,能夠考慮⼀下⽤C、Fortran或Cython(等下會稍微介紹⼀下)來編寫代碼。在⼯做中常常會⽤到Cython(http://cython.org),由於它不⽤花費太多精⼒就能獲得C語⾔那樣的性能。

一、連續內存的重要性
這個話題有點超範圍,但仍是要提⼀下,由於在某些應⽤場景中,數組的內存佈局能夠對計算速度形成極⼤的影響。這是由於性能差異在⼀定程度上跟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

相關文章
相關標籤/搜索