hans

hans

【Python】二、Numpy——《用Python做科學計算》


標準安裝的 Python 中用列表 (list) 保存一組值,可以用來當作數組使用,不過由於列表的元素可以是任何對象,因此列表中所保存的是對象的指針。這樣為了保存一個簡單的 [1,2,3],需要有 3 個指針和三個整數對象。對於數值運算來說這種結構顯然比較浪費內存和 CPU 計算時間。
此外 Python 還提供了一個 array 模塊,array 對象和列表不同,它直接保存數值,和 C 語言的一維數組比較類似。但是由於它不支持多維,也沒有各種運算函數,因此也不適合做數值運算。
NumPy 的誕生彌補了這些不足,NumPy 提供了兩種基本的對象:ndarray(N-dimensional array object)和
ufunc(universal function
object)。ndarray (下文統一稱之為數組) 是存儲單一數據類型的多維數組,而 ufunc 則是能夠對數組進行處理的函數。#

這裡我想講一下列表(list),元組 or 組元(tuple)和數組(array)的區別:
列表: 元素可以是任何類型的對象,必須初始化,並且初始化後內容可以修改。當所有元素對象一樣的時候可以當數組使用。

>>> a = ['a', "hello", 1, 2.2]
>>> a
['a', 'hello', 1, 2.2]
>>> a[1] = "HELLO"
>>> a
['a', 'HELLO', 1, 2.2]

元組: 元素可以是任何類型的對象,必須初始化,並且初始化後內容不可以修改。

>>> b = ('a', "hello", 1, 2.2)
>>> b
('a', 'hello', 1, 2.2)
>>> b[2] = "HELLO"
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'tuple' object does not support item assignment

數組: 元素只可以是一種類型的對象,非必須初始化,並且初始化後內容可以修改。

>>> c = np.array([1,2,3])
>>> d = np.array(['a','b','c'])
>>> e = np.array([1,2,'c','d'])
>>> e
array(['1', '2', 'c', 'd'],
      dtype='|S21')

看到沒,數組 e 中將 1,2, 默認為了 string 類型與後面的‘c',’d' 保持一致。
善於觀察的話可以看到元組是用的小括號(),列表和數組用的中括號 []。#

1. 函數庫導入#

import numpy as np

導入 numpy 函數庫,並在後面用縮略 np 表示。

2. 創建#

a,b 是一維數組,c 是多維數組。

>>> a = np.array([1, 2, 3, 4])
>>> b = np.array((5, 6, 7, 8))
>>> c = np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]])
>>> b
array([5, 6, 7, 8])
>>> c
array([[1, 2, 3, 4],
       [4, 5, 6, 7],
       [7, 8, 9, 10]])
>>> c.dtype
dtype('int32')

數組的大小可以通過其 shape 屬性獲得:

>>> a.shape
(4,)
>>> c.shape
(3, 4)

數組 a 的 shape 只有一個元素,因此它是一維數組。而數組 c 的 shape 有兩個元素,因此它是二維數組,其中第 0 軸(3 個橫軸)的長度為 3,第 1 軸(4 個縱軸)的長度為 4。還可以通過修改數組的 shape 屬性,在保持數組元素個數不變的情況下,改變數組每個軸的長度。下面的例子將數組 c 的 shape 改為 (4,3),注意從 (3,4) 改為 (4,3) 並不是對數組進行轉置,而只是改變每個軸的大小,數組元素在內存中的位置並沒有改變:

>>> c.shape = 4,3
>>> c
array([[ 1,  2,  3],
       [ 4,  4,  5],
       [ 6,  7,  7],
       [ 8,  9, 10]])

看上面的數組 a 是一維數組,a.shape 只有一個數,這個數字表示數組裡面元素的個數。我們還可以用多維數組的表達方法表示 a。注意觀察區別:

>>> a = np.array([1, 2, 3, 4])
>>> a
array([1, 2, 3, 4])
>>> a.shape
(4,)
>>> a.shape = 1,4
>>> a
array([[1, 2, 3, 4]])
>>> a.shape
(1, 4)

當某個軸的元素為 - 1 時,將根據數組元素的個數自動計算此軸的長度,因此下面的程序將數組 c 的 shape 改為 (2,6):

>>> c.shape = 2,-1
>>> c
array([[ 1,  2,  3,  4,  4,  5],
       [ 6,  7,  7,  8,  9, 10]])

使用數組的 reshape 方法,可以創建一個改變了尺寸的新數組,原數組的 shape 保持不變:

>>> d = a.reshape((2,2))
>>> d
array([[1, 2],
       [3, 4]])
>>> a
array([1, 2, 3, 4])

數組 a 和 d 其實共享數據存儲內存區域,因此修改其中任意一個數組的元素都會同時修改另外一個數組的內容:

>>> a[1] = 100 # 將數組a的第一個元素改為100
>>> d # 注意數組d中的2也被改變了
array([[  1, 100],
       [  3,   4]])

數組的元素類型可以通過 dtype 屬性獲得。上面例子中的參數序列的元素都是整數,因此所創建的數組的元素類型也是整數,並且是 32bit 的長整型。可以通過 dtype 參數在創建時指定元素類型:

>>> np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]], dtype=np.float)
array([[  1.,   2.,   3.,   4.],
       [  4.,   5.,   6.,   7.],
       [  7.,   8.,   9.,  10.]])
>>> np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]], dtype=np.complex)
array([[  1.+0.j,   2.+0.j,   3.+0.j,   4.+0.j],
       [  4.+0.j,   5.+0.j,   6.+0.j,   7.+0.j],
       [  7.+0.j,   8.+0.j,   9.+0.j,  10.+0.j]])

上面都是通過 array 函數創建數組。下面是三種專門用來創建數組的函數。
arange 函數類似於 python 的 range 函數,通過指定開始值、終值和步長來創建一維數組, 注意數組不包括終值 :

>>> np.arange(0,1,0.1)
array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9])

linspace 函數通過指定開始值、終值和元素個數來創建一維數組,可以通過 endpoint 關鍵字指定是否包括終值,缺省(默認)設置是 包括終值 :

>>> np.linspace(0, 1, 12)
array([ 0.        ,  0.09090909,  0.18181818,  0.27272727,  0.36363636,
        0.45454545,  0.54545455,  0.63636364,  0.72727273,  0.81818182,
        0.90909091,  1.        ])

logspace 函數和 linspace 類似,不過它創建等比數列,下面的例子產生 1 (10^0) 到 100 (10^2)、有 20 個元素的等比數列:

>>> np.logspace(0, 2, 20)
array([   1.        ,    1.27427499,    1.62377674,    2.06913808,
          2.6366509 ,    3.35981829,    4.2813324 ,    5.45559478,
          6.95192796,    8.8586679 ,   11.28837892,   14.38449888,
         18.32980711,   23.35721469,   29.76351442,   37.92690191,
         48.32930239,   61.58482111,   78.47599704,  100.        ])

此外,使用 frombuffer, fromstring, fromfile 等函數可以從字節序列創建數組,下面以 fromstring 為例:

>>> s = "abcdefgh"

Python 的字符串實際上是字節序列,每個字符占一個字節,因此如果從字符串 s 創建一個 8bit 的整數數組的話,所得到的數組正好就是字符串中每個字符的 ASCII 編碼:

>>> np.fromstring(s, dtype=np.int8)
array([ 97,  98,  99, 100, 101, 102, 103, 104], dtype=int8)

如果從字符串 s 創建 16bit 的整數數組,那麼兩個相鄰的字節就表示一個整數,把字節 98 和字節 97 當作一個 16 位的整數,它的值就是 98*256+97 =
25185。可以看出內存中是以 little endian (低位字節在前) 方式保存數據的。

>>> np.fromstring(s, dtype=np.int16)
array([25185, 25699, 26213, 26727], dtype=int16)
>>> 98*256+97
25185

我們可以寫一個 Python 的函數,它將數組下標轉換為數組中對應的值,然後使用此函數創建數組:

>>> def func(i): #定義函數,函數名func,參數i
...   return i%4+1 #返回i除以4的余數加一
...
>>> np.fromfunction(func, (10,))
array([ 1.,  2.,  3.,  4.,  1.,  2.,  3.,  4.,  1.,  2.])

fromfunction 第二個參數 n 可以理解為從 0 開始一共有 n 個取值,賦值給前面的 function。
這個例子是 i 先取 0, 取余加一後返回數值 1
然後 i 取 1, 取余加一後返回數值 2
依次類推
最後 i 取 9,取余加一後返回數值 2#

fromfunction 函數的第一個參數為計算每個數組元素的函數,第二個參數為數組的大小 (shape),因為它支持多維數組,所以第二個參數必須是一個序列,本例中用 (10,) 創建一個 10 元素的一維數組。

下面的例子創建一個二維數組表示九九乘法表,輸出的數組 a 中的每個元素 a [i, j] 都等於 func2 (i, j):

>>> def func2(i, j): #定義函數,函數名func,參數i,j
...     return (i+1) * (j+1) #函數返回值
...
>>> a = np.fromfunction(func2, (9,9))
>>> a
array([[  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.],
       [  2.,   4.,   6.,   8.,  10.,  12.,  14.,  16.,  18.],
       [  3.,   6.,   9.,  12.,  15.,  18.,  21.,  24.,  27.],
       [  4.,   8.,  12.,  16.,  20.,  24.,  28.,  32.,  36.],
       [  5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.],
       [  6.,  12.,  18.,  24.,  30.,  36.,  42.,  48.,  54.],
       [  7.,  14.,  21.,  28.,  35.,  42.,  49.,  56.,  63.],
       [  8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,  72.],
       [  9.,  18.,  27.,  36.,  45.,  54.,  63.,  72.,  81.]]) #可以理解為每一點的下標相乘

3. 存取元素#

數組元素的存取方法和 Python 的標準方法相同:

>>> a = np.arange(10)
>>> a[5]    # 用整數作為下標可以獲取數組中的某個元素
5
>>> a[3:5]  # 用範圍作為下標獲取數組的一個切片,包括a[3]不包括a[5]
array([3, 4])
>>> a[:5]   # 省略開始下標,表示從a[0]開始
array([0, 1, 2, 3, 4])
>>> a[:-1]  # 下標可以使用負數,表示從數組最後一個元素往前數
array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> a[2:4] = 100,101    # 下標還可以用來修改元素的值
>>> a
array([  0,   1, 100, 101,   4,   5,   6,   7,   8,   9])
>>> a[1:-1:2]   # 範圍中的第三個參數表示步長,2表示隔一個元素取一個元素
array([  1, 101,   5,   7])
>>> a[::-1] # 省略範圍的開始下標和結束下標,步長為-1,整個數組頭尾顛倒
array([  9,   8,   7,   6,   5,   4, 101, 100,   1,   0])
>>> a[5:1:-2] # 步長為負數時,開始下標必須大於結束下標
array([  5, 101])

和 Python 的列表序列不同,通過下標範圍獲取的新數組是原始數組的一個視圖。它與原始數組共享同一塊數據空間:

>>> b = a[3:7] # 通過下標範圍產生一個新的數組b,b和a共享同一塊數據空間
>>> b
array([101,   4,   5,   6])
>>> b[2] = -10 # 將b的第2個元素修改為-10
>>> b
array([101,   4, -10,   6])
>>> a # a的第5個元素也被修改為10
array([  0,   1, 100, 101,   4, -10,   6,   7,   8,   9])

除了使用下標範圍存取元素之外,NumPy 還提供了兩種存取元素的高級方法。
使用整數序列
當使用整數序列對數組元素進行存取時,將使用整數序列中的每個元素作為下標,整數序列可以是列表或者數組。使用整數序列作為下標獲得的數組不和原始數組共享數據空間。

>>> x = np.arange(10,1,-1)
>>> x
array([10,  9,  8,  7,  6,  5,  4,  3,  2])
>>> x[[3, 3, 1, 8]] # 獲取x中的下標為3, 3, 1, 8的4個元素,組成一個新的數組
array([7, 7, 9, 2])
>>> b = x[np.array([3,3,-3,8])]  #下標可以是負數
>>> b[2] = 100
>>> b
array([7, 7, 100, 2])
>>> x   # 由於b和x不共享數據空間,因此x中的值並沒有改變
array([10,  9,  8,  7,  6,  5,  4,  3,  2])
>>> x[[3,5,1]] = -1, -2, -3 # 整數序列下標也可以用來修改元素的值
>>> x
array([10, -3,  8, -1,  6, -2,  4,  3,  2])

使用布爾數組
當使用布爾數組 b 作為下標存取數組 x 中的元素時,將收集數組 x 中所有在數組 b 中對應下標為 True 的元素。使用布爾數組作為下標獲得的數組不和原始數組共享數據空間,注意這種方式只對應於布爾數組,不能使用布爾列表。(列表,元組,數組區別還記得?)

>>> x = np.arange(5,0,-1)
>>> x
array([5, 4, 3, 2, 1])
>>> x[np.array([True, False, True, False, False])]
>>> # 布爾數組中下標為0,2的元素為True,因此獲取x中下標為0,2的元素
array([5, 3])
>>> x[[True, False, True, False, False]]
>>> # 如果是布爾列表,則把True當作1, False當作0,按照整數序列方式獲取x中的元素
array([4, 5, 4, 5, 5])
>>> x[np.array([True, False, True, True])]
>>> # 布爾數組的長度不夠時,不夠的部分都當作False
array([5, 3, 2])
>>> x[np.array([True, False, True, True])] = -1, -2, -3
>>> # 布爾數組下標也可以用來修改元素
>>> x
array([-1,  4, -2, -3,  1])

布爾數組一般不是手工產生,而是使用布爾運算的 ufunc 函數產生,關於 ufunc 函數請參照 ufunc 運算 一節。

>>> x = np.random.rand(10) # 產生一個長度為10,元素值為0-1的隨機數的數組
>>> x
array([ 0.72223939,  0.921226  ,  0.7770805 ,  0.2055047 ,  0.17567449,
        0.95799412,  0.12015178,  0.7627083 ,  0.43260184,  0.91379859])
>>> x>0.5
>>> # 數組x中的每個元素和0.5進行大小比較,得到一個布爾數組,True表示x中對應的值大於0.5
array([ True,  True,  True, False, False,  True, False,  True, False,  True], dtype=bool)
>>> x[x>0.5]
>>> # 使用x>0.5返回的布爾數組收集x中的元素,因此得到的結果是x中所有大於0.5的元素的數組
array([ 0.72223939,  0.921226  ,  0.7770805 ,  0.95799412,  0.7627083 ,
        0.91379859])

4. 多維數組#

多維數組的存取和一維數組類似,因為多維數組有多個軸,因此它的下標需要用多個值來表示,NumPy 採用組元 (tuple) 作為數組的下標。如圖所示,a 為一個 6x6 的數組,圖中用顏色區分了各個下標以及其對應的選擇區域。
1668631158233.jpg


a [2::2,::2] 這一句分解開這麼看,首先 “2:” 表示從索引為 2 的橫軸(0 軸)開始到末尾,然後 “:2” 表示隔兩行取一次。其實就是取 2,4 行。後面第一個 “:” 表示所有列(1 軸),“:2” 表示隔兩列取一次。就是取 0,2,4 三列。

比如 24 的下標是 (2,4),2 就是 0 軸上的索引,4 就是 1 軸上的索引。


數組 a 是一個加法表,縱軸的值為 0, 10, 20, 30, 40, 50;橫軸的值為 0, 1, 2, 3, 4, 5。縱軸的每個元素都和橫軸的每個元素求和。

>>> np.arange(0, 60, 10).reshape(-1, 1) + np.arange(0, 6)
array([[ 0,  1,  2,  3,  4,  5],
       [10, 11, 12, 13, 14, 15],
       [20, 21, 22, 23, 24, 25],
       [30, 31, 32, 33, 34, 35],
       [40, 41, 42, 43, 44, 45],
       [50, 51, 52, 53, 54, 55]])

多維數組同樣也可以使用整數序列和布爾數組進行存取。
1668631278744.jpg
a[(0,1,2,3,4),(1,2,3,4,5)] :
用於存取數組的下標和仍然是一個有兩個元素的組元,組元中的每個元素都是整數序列,分別對應數組的第 0 軸和第 1 軸。從兩個序列的對應位置取出兩個整數組成下標:
a[0,1], a[1,2], ..., a[4,5]。
a [3:, [0, 2, 5]] : 下標中的第 0 軸是一個範圍,它選取第 3 行之後的所有行;第 1 軸是整數序列,它選取第 0, 2, 5 三列。
a [mask, 2] : 下標的第 0 軸是一個布爾數組,它選取第 0,2,5 行;第 1 軸是一個整數,選取第 2 列。

5. 數組結構#

假設我們需要定義一個結構數組,它的每個元素都有 name, age 和 weight 字段。在 NumPy 中可以如下定義:

import numpy as np
persontype = np.dtype({
    'names':['name', 'age', 'weight'],
    'formats':['S32','i', 'f']})
a = np.array([("Zhang",32,75.5),("Wang",24,65.2)],
    dtype=persontype)

我們先創建一個 dtype 對象 persontype,通過其字典參數描述結構類型的各個字段。字典有兩個關鍵字:names,formats。每個關鍵字對應的值都是一個列表。names 定義結構中的每個字段名,而 formats 則定義每個字段的類型:
S32 : 32 個字節的字符串類型,由於結構中的每個元素的大小必須固定,因此需要指定字符串的長度
i : 32bit 的整數類型,相當於 np.int32
f : 32bit 的單精度浮點數類型,相當於 np.float32
然後我們調用 array 函數創建數組,通過關鍵字參數 dtype=persontype,
指定所創建的數組的元素類型為結構 persontype。運行上面程序之後,我們可以在 IPython 中執行如下的語句查看數組 a 的元素類型

>>> a.dtype
dtype([('name', '|S32'), ('age', '<i4'), ('weight', '<f4')])

結構數組的存取方式和一般數組相同,通過下標能夠取得其中的元素,注意元素的值看上去像是組元,實際上它是一個結構:

>>> a[0]
('Zhang', 32, 75.5)
>>> a[0].dtype
dtype([('name', '|S32'), ('age', '<i4'), ('weight', '<f4')])

a [0] 是一個結構元素,它和數組 a 共享內存數據,因此可以通過修改它的字段,改變原始數組中的對應字段:

>>> c = a[1]
>>> c["name"] = "Li"
>>> a[1]["name"]
"Li"

結構像字典一樣可以通過字符串下標獲取其對應的字段值:

>>> a[0]["name"]
'Zhang'

我們不但可以獲得結構元素的某個字段,還可以直接獲得結構數組的字段,它返回的是原始數組的視圖,因此可以通過修改 b [0] 改變 a [0]["age"]:

>>> b=a[:]["age"] # 或者a["age"]
>>> b
array([32, 24])
>>> b[0] = 40
>>> a[0]["age"]
40

通過調用 a.tostring 或者 a.tofile 方法,可以直接輸出數組 a 的二進制形式:

>>> a.tofile("test.bin")

6.ufunc 運算#

ufunc 是 universal
function 的縮寫,它是一種能對數組的每個元素進行操作的函數。NumPy 內置的許多 ufunc 函數都是在 C 語言級別實現的,因此它們的計算速度非常快。讓我們來看一個例子:

>>> x = np.linspace(0, 2*np.pi, 10)
# 對數組x中的每個元素進行正弦計算,返回一個同樣大小的新數組
>>> y = np.sin(x)
>>> y
array([  0.00000000e+00,   6.42787610e-01,   9.84807753e-01,
         8.66025404e-01,   3.42020143e-01,  -3.42020143e-01,
        -8.66025404e-01,  -9.84807753e-01,  -6.42787610e-01,
        -2.44921271e-16])

先用 linspace 產生一個從 0 到 2*PI 的等距離的 10 個數,然後將其傳遞給 sin 函數,由於 np.sin 是一個 ufunc 函數,因此它對 x 中的每個元素求正弦值,然後將結果返回,並且賦值給 y。計算之後 x 中的值並沒有改變,而是新創建了一個數組保存結果。如果我們希望將 sin 函數所計算的結果直接覆蓋到數組 x 上去的話,可以將要被覆蓋的數組作為第二個參數傳遞給 ufunc 函數。例如:

>>> t = np.sin(x,x)
>>> x
array([  0.00000000e+00,   6.42787610e-01,   9.84807753e-01,
         8.66025404e-01,   3.42020143e-01,  -3.42020143e-01,
        -8.66025404e-01,  -9.84807753e-01,  -6.42787610e-01,
        -2.44921271e-16])
>>> id(t) == id(x)
True

sin 函數的第二個參數也是 x,那麼它所做的事情就是對 x 中的每給值求正弦值,並且把結果保存到 x 中的對應的位置中。此時函數的返回值仍然是整個計算的結果,只不過它就是 x,因此兩個變量的 id 是相同的 (變量 t 和變量 x 指向同一塊內存區域)。
下面這個小程序,比較了一下 numpy.math 和 Python 標準庫的 math.sin 的計算速度:

import time
import math
import numpy as np

x = [i * 0.001 for i in xrange(1000000)]
start = time.clock()
for i, t in enumerate(x):
    x[i] = math.sin(t)
print "math.sin:", time.clock() - start

x = [i * 0.001 for i in xrange(1000000)]
x = np.array(x)
start = time.clock()
np.sin(x,x)
print "numpy.sin:", time.clock() - start

# 輸出
# math.sin: 1.15426932753
# numpy.sin: 0.0882399858083

計算 100 萬次正弦值,numpy.sin 比 math.sin 快 10 倍多。這得利於 numpy.sin 在 C 語言級別的循環計算。numpy.sin 同樣也支持對單個數值求正弦,例如:numpy.sin (0.5)。不過值得注意的是,對單個數的計算 math.sin 則比 numpy.sin 快得多了,讓我們看下面這個測試程序:

x = [i * 0.001 for i in xrange(1000000)]
start = time.clock()
for i, t in enumerate(x):
    x[i] = np.sin(t)
print "numpy.sin loop:", time.clock() - start

# 輸出
# numpy.sin loop: 5.72166965355

請注意 numpy.sin 的計算速度只有 math.sin 的 1/5。這是因為 numpy.sin 為了同時支持數組和單個值的計算,其 C 語言的內部實現要比 math.sin 複雜很多,如果我們同樣在 Python 級別進行循環的話,就會看出其中的差別了。此外,numpy.sin 返回的數的類型和 math.sin 返回的類型有所不同,math.sin 返回的是 Python 的標準 float 類型,而 numpy.sin 則返回一個 numpy.float64 類型:

>>> type(math.sin(0.5))
<type 'float'>
>>> type(np.sin(0.5))
<type 'numpy.float64'>

通過上面的例子我們了解了如何最有效率地使用 math 庫和 numpy 庫中的數學函數。
NumPy 中有眾多的 ufunc 函數為我們提供各式各樣的計算。除了 sin 這種單輸入函數之外,還有許多多個輸入的函數,add 函數就是一個最常用的例子。先來看一個例子:

>>> a = np.arange(0,4)
>>> a
array([0, 1, 2, 3])
>>> b = np.arange(1,5)
>>> b
array([1, 2, 3, 4])
>>> np.add(a,b)
array([1, 3, 5, 7])
>>> np.add(a,b,a)
array([1, 3, 5, 7])
>>> a
array([1, 3, 5, 7])

add 函數返回一個新的數組,此數組的每個元素都為兩個參數數組的對應元素之和。它接受第 3 個參數指定計算結果所要寫入的數組,如果指定的話,add 函數就不再產生新的數組。
由於 Python 的操作符重載功能,計算兩個數組相加可以簡單地寫為 a+b,而 np.add (a,b,a) 則可以用 a+=b 來表示。下面是數組的運算符和其對應的 ufunc 函數的一個列表,注意除號 "/" 的意義根據是否激活__future__.division 有所不同。
y = x1 + x2: add(x1, x2 [, y])
y = x1 - x2: subtract(x1, x2 [, y])
y = x1 * x2: multiply (x1, x2 [, y])
y = x1 /x2: divide (x1, x2 [, y]), 如果兩個數組的元素為整數,那麼用整數除法
y = x1 /x2: true divide (x1, x2 [, y]), 總是返回精確的商
y = x1 //x2: floor divide (x1, x2 [, y]), 總是對返回值取整
y = -x: negative(x [,y])
y = x1**x2: power(x1, x2 [, y])
y = x1 % x2: remainder(x1, x2 [, y]), mod(x1, x2, [, y])
數組對象支持這些操作符,極大地簡化了算式的編寫,不過要注意如果你的算式很複雜,並且要運算的數組很大的話,會因為產生大量的中間結果而降低程序的運算效率。例如:假設 a
b c 三個數組採用算式 x=a*b+c 計算,那麼它相當於:

t = a * b
x = t + c
del t

也就是說需要產生一個數組 t 保存乘法的計算結果,然後再產生最後的結果數組 x。我們可以通過手工將一個算式分解為 x = a*b; x +=
c,以減少一次內存分配。

通過組合標準的 ufunc 函數的調用,可以實現各種算式的數組計算。不過有些時候這種算式不易編寫,而針對每個元素的計算函數卻很容易用 Python 實現,這時可以用 frompyfunc 函數將一個計算單個元素的函數轉換成 ufunc 函數。這樣就可以方便地用所產生的 ufunc 函數對數組進行計算了。讓我們來看一個例子。

我們想用一個分段函數描述三角波,三角波的樣子如圖所示:
1668631406858.jpg
我們很容易根據上圖所示寫出如下的計算三角波某點 y 坐標的函數:

def triangle_wave(x, c, c0, hc):
    x = x - int(x) # 三角波的周期為1,因此只取x坐標的小數部分進行計算
    if x >= c: r = 0.0
    elif x < c0: r = x / c0 * hc
    else: r = (c-x) / (c-c0) * hc
    return r

顯然 triangle_wave 函數只能計算單個數值,不能對數組直接進行處理。
另外還有一種方法:

def triangle_func(c, c0, hc):
    def trifunc(x):
        x = x - int(x) # 三角波的周期為1,因此只取x坐標的小數部分進行計算
        if x >= c: r = 0.0
        elif x < c0: r = x / c0 * hc
        else: r = (c-x) / (c-c0) * hc
        return r

    # 用trifunc函數創建一個ufunc函數,可以直接對數組進行計算, 不過通過此函數
    # 計算得到的是一個Object數組,需要進行類型轉換
    return np.frompyfunc(trifunc, 1, 1)

y2 = triangle_func(0.6, 0.4, 1.0)(x)

我們通過函數 triangle_func 包裝三角波的三個參數,在其內部定義一個計算三角波的函數 trifunc,trifunc 函數在調用時會採用 triangle_func 的參數進行計算。最後 triangle_func 返回用 frompyfunc 轉換結果。

值得注意的是用 frompyfunc 得到的函數計算出的數組元素的類型為 object,因為 frompyfunc 函數無法保證 Python 函數返回的數據類型都完全一致。因此還需要再次
y2.astype (np.float64) 將其轉換為雙精度浮點數組。

7. 廣播#

當我們使用 ufunc 函數對兩個數組進行計算時,ufunc 函數會對這兩個數組的對應元素進行計算,因此它要求這兩個數組有相同的大小 (shape 相同)。如果兩個數組的 shape 不同的話,會進行如下的廣播 (broadcasting) 處理:
1. 讓所有輸入數組都向其中 shape 最長的數組看齊,shape 中不足的部分都通過在前面加 1 補齊
2. 輸出數組的 shape 是輸入數組 shape 的各個軸上的最大值
3. 如果輸入數組的某個軸和輸出數組的對應軸的長度相同或者其長度為 1 時,這個數組能夠用來計算,否則出錯
4. 當輸入數組的某個軸的長度為 1 時,沿著此軸運算時都用此軸上的第一組值
上述 4 條規則理解起來可能比較費勁,讓我們來看一個實際的例子。

先創建一個二維數組 a,其 shape 為 (6,1):

>>> a = np.arange(0, 60, 10).reshape(-1, 1)
>>> a
array([[ 0], [10], [20], [30], [40], [50]])
>>> a.shape
(6, 1)
再創建一維數組b,其shape為(5,):
>>> b = np.arange(0, 5)
>>> b
array([0, 1, 2, 3, 4])
>>> b.shape
(5,)

計算 a 和 b 的和,得到一個加法表,它相當於計算 a,b 中所有元素組的和,得到一個 shape 為 (6,5) 的數組:

>>> c = a + b
>>> c
array([[ 0,  1,  2,  3,  4],
       [10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24],
       [30, 31, 32, 33, 34],
       [40, 41, 42, 43, 44],
       [50, 51, 52, 53, 54]])
>>> c.shape
(6, 5)

由於 a 和 b 的 shape 長度 (也就是 ndim 屬性) 不同,根據規則 1,需要讓 b 的 shape 向 a 對齊,於是將 b 的 shape 前面加 1,補齊為 (1,5)。相當於做了如下計算:

>>> b.shape=1,5
>>> b
array([[0, 1, 2, 3, 4]])

這樣加法運算的兩個輸入數組的 shape 分別為 (6,1) 和 (1,5),根據規則 2,輸出數組的各個軸的長度為輸入數組各個軸上的長度的最大值,可知輸出數組的 shape 為 (6,5)。

由於 b 的第 0 軸上的長度為 1,而 a 的第 0 軸上的長度為 6,因此為了讓它們在第 0 軸上能夠相加,需要將 b 在第 0 軸上的長度擴展為 6,這相當於:

>>> b = b.repeat(6,axis=0)
>>> b
array([[0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4]])

由於 a 的第 1 軸的長度為 1,而 b 的第一軸長度為 5,因此為了讓它們在第 1 軸上能夠相加,需要將 a 在第 1 軸上的長度擴展為 5,這相當於:

>>> a = a.repeat(5, axis=1)
>>> a
array([[ 0,  0,  0,  0,  0],
       [10, 10, 10, 10, 10],
       [20, 20, 20, 20, 20],
       [30, 30, 30, 30, 30],
       [40, 40, 40, 40, 40],
       [50, 50, 50, 50, 50]])

經過上述處理之後,a 和 b 就可以按對應元素進行相加運算了。

當然,numpy 在執行 a+b 運算時,其內部並不會真正將長度為 1 的軸用 repeat 函數進行擴展,如果這樣做的話就太浪費空間了。

由於這種廣播計算很常用,因此 numpy 提供了一個快速產生如上面 a,b 數組的方法: ogrid 對象:

>>> x,y = np.ogrid[0:5,0:5]
>>> x
array([[0],
       [1],
       [2],
       [3],
       [4]])
>>> y
array([[0, 1, 2, 3, 4]])

ogrid 是一個很有趣的對象,它像一個多維數組一樣,用切片組元作為下標進行存取,返回的是一組可以用來廣播計算的數組。其切片下標有兩種形式:
1. 開始值:結束值:步長,和 np.arange (開始值,結束值,步長) 類似
2. 開始值:結束值:長度 j,當第三個參數為虛數時,它表示返回的數組的長度,和 np.linspace (開始值,結束值,長度) 類似:

>>> x, y = np.ogrid[0:1:4j, 0:1:3j]
>>> x
array([[ 0.        ],
       [ 0.33333333],
       [ 0.66666667],
       [ 1.        ]])
>>> y
array([[ 0. ,  0.5,  1. ]])

利用 ogrid 的返回值,我能很容易計算 x, y 網格面上各點的值,或者 x, y, z 網格體上各點的值。下面是繪製三維曲面 x * exp (x2 -
y
2) 的程序:

import numpy as np
from enthought.mayavi import mlab

x, y = np.ogrid[-2:2:20j, -2:2:20j]
z = x * np.exp( - x**2 - y**2)

pl = mlab.surf(x, y, z, warp_scale="auto")
mlab.axes(xlabel='x', ylabel='y', zlabel='z')
mlab.outline(pl)

此程序使用 mayavi 的 mlab 庫快速繪製如圖所示的 3D 曲面
1668631471880.jpg

8.ufunc 的方法#

reduce
方法和 Python 的 reduce 函數類似,它沿著 axis 軸對 array 進行操作,相當於將運算符插入到沿 axis 軸的所有子數組或者元素當中。
例如:

>>> np.add.reduce([1,2,3]) # 1 + 2 + 3
6
>>> np.add.reduce([[1,2,3],[4,5,6]], axis=1) # 1,4 + 2,5 + 3,6
array([ 6, 15])

accumulate 方法和 reduce 方法類似,只是它返回的數組和輸入的數組的 shape 相同,保存所有的中間計算結果:

>>> np.add.accumulate([1,2,3])
array([1, 3, 6])
>>> np.add.accumulate([[1,2,3],[4,5,6]], axis=1)
array([[ 1,  3,  6],
       [ 4,  9, 15]])

outer 方法,.outer (a,b) 方法:

>>> np.multiply.outer([1,2,3,4,5],[2,3,4])
array([[ 2,  3,  4],
       [ 4,  6,  8],
       [ 6,  9, 12],
       [ 8, 12, 16],
       [10, 15, 20]])

可以看出通過 outer 方法計算的結果是如下的乘法表:

#    2, 3, 4
# 1
# 2
# 3
# 4
# 5

如果將這兩個數組按照等同程序一步一步的計算的話,就會發現乘法表最終是通過廣播的方式計算出來的。

9. 矩陣運算#

矩陣的乘積可以使用 dot 函數進行計算。對於二維數組,它計算的是矩陣乘積,對於一維數組,它計算的是其點積。當需要將一維數組當作列矢量或者行矢量進行矩陣運算時,推薦先使用 reshape 函數將一維數組轉換為二維數組:

>>> a = array([1, 2, 3])
>>> a.reshape((-1,1))
array([[1],
       [2],
       [3]])
>>> a.reshape((1,-1))
array([[1, 2, 3]])

dot : 對於兩個一維的數組,計算的是這兩個數組對應下標元素的乘積和 (數學上稱之為內積);其實就是點乘!

10. 文件存儲#

NumPy 提供了多種文件操作函數方便我們存取數組內容。文件存取的格式分為兩類:二進制和文本。而二進制格式的文件又分為 NumPy 專用的格式化二進制類型和無格式類型。

使用數組的方法函數 tofile 可以方便地將數組中數據以二進制的格式寫進文件。tofile 輸出的數據沒有格式,因此用 numpy.fromfile 讀回來的時候需要自己格式化數據:

>>> a = np.arange(0,12)
>>> a.shape = 3,4
>>> a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> a.tofile("a.bin")
>>> b = np.fromfile("a.bin", dtype=np.float) # 按照float類型讀入數據
>>> b # 讀入的數據是錯誤的
array([  2.12199579e-314,   6.36598737e-314,   1.06099790e-313,
         1.48539705e-313,   1.90979621e-313,   2.33419537e-313])
>>> a.dtype # 查看a的dtype
dtype('int32')
>>> b = np.fromfile("a.bin", dtype=np.int32) # 按照int32類型讀入數據
>>> b # 數據是一維的
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
>>> b.shape = 3, 4 # 按照a的shape修改b的shape
>>> b # 這次終於正確了
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

從上面的例子可以看出,需要在讀入的時候設置正確的 dtype 和 shape 才能保證數據一致。並且 tofile 函數不管數組的排列順序是 C 語言格式的還是 Fortran 語言格式的,統一使用 C 語言格式輸出。

此外如果 fromfile 和 tofile 函數調用時指定了 sep 關鍵字參數的話,數組將以文本格式輸入輸出。

numpy.load 和 numpy.save 函數以 NumPy 專用的二進制類型保存數據,這兩個函數會自動處理元素類型和 shape 等信息,使用它們讀寫數組就方便多了,但是 numpy.save 輸出的文件很難和其它語言編寫的程序讀入:

>>> np.save("a.npy", a)
>>> c = np.load( "a.npy" )
>>> c
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

如果你想將多個數組保存到一個文件中的話,可以使用 numpy.savez 函數。savez 函數的第一個參數是文件名,其後的參數都是需要保存的數組,也可以使用關鍵字參數為數組起一個名字,非關鍵字參數傳遞的數組會自動起名為 arr_0,
arr_1,
...。savez 函數輸出的是一個壓縮文件 (擴展名為 npz),其中每個文件都是一個 save 函數保存的 npy 文件,文件名對應於數組名。load 函數自動識別 npz 文件,並且返回一個類似於字典的對象,可以通過數組名作為關鍵字獲取數組的內容:

>>> a = np.array([[1,2,3],[4,5,6]])
>>> b = np.arange(0, 1.0, 0.1)
>>> c = np.sin(b)
>>> np.savez("result.npz", a, b, sin_array = c)
>>> r = np.load("result.npz")
>>> r["arr_0"] # 數組a
array([[1, 2, 3],
       [4, 5, 6]])
>>> r["arr_1"] # 數組b
array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9])
>>> r["sin_array"] # 數組c
array([ 0.        ,  0.09983342,  0.19866933,  0.29552021,  0.38941834,
        0.47942554,  0.56464247,  0.64421769,  0.71735609,  0.78332691])

如果你用解壓軟件打開 result.npz 文件的話,會發現其中有三個文件:arr_0.npy, arr_1.npy,
sin_array.npy,其中分別保存著數組 a, b, c 的內容。

使用 numpy.savetxt 和 numpy.loadtxt 可以讀寫 1 維和 2 維的數組:

>>> a = np.arange(0,12,0.5).reshape(4,-1)
>>> np.savetxt("a.txt", a) # 缺省按照'%.18e'格式保存數據,以空格分隔
>>> np.loadtxt("a.txt")
array([[  0. ,   0.5,   1. ,   1.5,   2. ,   2.5],
       [  3. ,   3.5,   4. ,   4.5,   5. ,   5.5],
       [  6. ,   6.5,   7. ,   7.5,   8. ,   8.5],
       [  9. ,   9.5,  10. ,  10.5,  11. ,  11.5]])
>>> np.savetxt("a.txt", a, fmt="%d", delimiter=",") #改為保存為整數,以逗號分隔
>>> np.loadtxt("a.txt",delimiter=",") # 讀入的時候也需要指定逗號分隔
array([[  0.,   0.,   1.,   1.,   2.,   2.],
       [  3.,   3.,   4.,   4.,   5.,   5.],
       [  6.,   6.,   7.,   7.,   8.,   8.],
       [  9.,   9.,  10.,  10.,  11.,  11.]])

本文參考自: http://old.sebug.net/paper/books/scipydoc/numpy_intro.html

載入中......
此文章數據所有權由區塊鏈加密技術和智能合約保障僅歸創作者所有。