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 npx = [i * 0.001 for i in range(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 range(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'>

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函數對數組進行計算了。讓我們來看一個例子。

三角波的樣子如下圖所示：

很容易根據上圖所示寫出如下的計算三角波某點y坐標的函數:

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

顯然 triangle_wave 函數只能計算單個數值，不能對數組直接進行處理。可以用下面的方法先使用列表包容 (List comprehension) ，計算出一個 list ，然后用 array 函數將列表轉換為數組：

x = np.linspace(0, 2, 1000) y = np.array([triangle_wave(t, 0.6, 0.4, 1.0) for t in x])

這種做法每次都需要使用列表包容語法調用函數，對于多維數組是很麻煩的。frompyfunc 栗子：

triangle_ufunc = np.frompyfunc(lambda x: triangle_wave(x, 0.6, 0.4, 1.0), 1, 1) y2 = triangle_ufunc(x)

frompyfunc的調用格式為frompyfunc(func, nin, nout)，其中func是計算單個元素的函數，nin是此函數的輸入參數的個數，nout是此函數的返回值的個數。雖然triangle_wave函數有4個參數，但是由于后三個c, c0, hc在整個計算中值都是固定的，因此所產生的ufunc函數其實只有一個參數。為了滿足這個條件，我們用一個lambda函數對triangle_wave的參數進行一次包裝。這樣傳入frompyfunc的函數就只有一個參數了。這樣子做，效率并不是太高，另外還有一種方法：

def triangle_func(c, c0, hc):def trifunc(x):x = x - int(x) # 三角波的周期為1，因此只取x坐標的小數部分進行計算if x >= c: r = 0.0elif x < c0: r = x / c0 * hcelse: r = (c-x) / (c-c0) * hcreturn 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)將其轉換為雙精度浮點數組。

廣播

當使用 ufunc 函數對兩個數組進行計算時，ufunc 函數會對這兩個數組的對應元素進行計算，因此它要求這兩個數組有相同的大小(shape相同)。如果兩個數組的shape不同的話，會進行如下的廣播(broadcasting)處理：

讓所有輸入數組都向其中shape最長的數組看齊，shape 中不足的部分都通過在前面加1補齊
輸出數組的shape是輸入數組shape的各個軸上的最大值
如果輸入數組的某個軸和輸出數組的對應軸的長度相同或者其長度為1時，這個數組能夠用來計算，否則出錯
當輸入數組的某個軸的長度為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 是一個很有趣的對象，它像一個多維數組一樣，用切片組元作為下標進行存取，返回的是一組可以用來廣播計算的數組。其切片下標有兩種形式：

開始值:結束值:步長，和 np.arange (開始值, 結束值, 步長) 類似

開始值:結束值:長度 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 為什么不是函數

根據 Python 的語法，只有在中括號中才能使用用冒號隔開的切片語法，如果 ogrid 是函數的話，那么這些切片必須使用 slice 函數創建，這顯然會增加代碼的長度。

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

import numpy as np from enthought.mayavi import mlabx, 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曲面。

ufunc的方法

ufunc 函數本身還有些方法，這些方法只對兩個輸入一個輸出的 ufunc 函數有效，其它的ufunc對象調用這些方法時會拋出ValueError異常。

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

<op>.reduce (array=, axis=0, dtype=None)

例如：

>>> 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]])

reduceat 方法計算多組reduce的結果，通過indices參數指定一系列reduce的起始和終了位置。reduceat的計算有些特別，讓我們通過一個例子來解釋一下：

>>> a = np.array([1,2,3,4]) >>> result = np.add.reduceat(a,indices=[0,1,0,2,0,3,0]) >>> result array([ 1, 2, 3, 3, 6, 4, 10])

對于indices中的每個元素都會調用reduce函數計算出一個值來，因此最終計算結果的長度和indices的長度相同。結果result數組中除最后一個元素之外，都按照如下計算得出：

if indices[i] < indices[i+1]:result[i] = np.reduce(a[indices[i]:indices[i+1]]) else:result[i] = a[indices[i]]

而最后一個元素如下計算:

np.reduce(a[indices[-1]:])

因此上面例子中，結果的每個元素如下計算而得：

1 : a[0] = 1 2 : a[1] = 2 3 : a[0] + a[1] = 1 + 2 3 : a[2] = 3 6 : a[0] + a[1] + a[2] = 1 + 2 + 3 = 6 4 : a[3] = 4 10: a[0] + a[1] + a[2] + a[4] = 1+2+3+4 = 10

可以看出 result[::2] 和a相等，而 result[1::2] 和 np.add.accumulate(a) 相等。

outer 方法，<op>.outer(a,b) 方法的計算等同于如下程序：

>>> a.shape += (1,)*b.ndim >>> <op>(a,b) >>> a = a.squeeze()

其中squeeze的功能是剔除數組a中長度為1的軸。如果你看不太明白這個等同程序的話，栗子：

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

轉載于:https://www.cnblogs.com/oneTOinf/p/10502718.html

與50位技術專家面對面20年技術見證，附贈技術全景圖

總結

以上是生活随笔為你收集整理的07.Numpy广播和ufunc的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。