ufunc運(yùn)算

ufunc 是 universal function 的縮寫，它是一種能對(duì)數(shù)組的每個(gè)元素進(jìn)行操作的函數(shù)。numpy 內(nèi)置的許多 ufunc 函數(shù)都是在C語(yǔ)言級(jí)別實(shí)現(xiàn)的，因此它們的計(jì)算速度非常快。栗子：

>>> x = np.linspace(0, 2*np.pi, 10) # 對(duì)數(shù)組x中的每個(gè)元素進(jìn)行正弦計(jì)算，返回一個(gè)同樣大小的新數(shù)組 >>> 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 產(chǎn)生一個(gè)從 0到2*PI 的等距離的10個(gè)數(shù)，然后將其傳遞給 sin 函數(shù)，由于 np.sin 是一個(gè) ufunc 函數(shù)，因此它對(duì) x 中的每個(gè)元素求正弦值，然后將結(jié)果返回，并且賦值給 y 。計(jì)算之后 x 中的值并沒(méi)有改變，而是新創(chuàng)建了一個(gè)數(shù)組保存結(jié)果。如果希望將 sin 函數(shù)所計(jì)算的結(jié)果直接覆蓋到數(shù)組 x 上去的話，可以將要被覆蓋的數(shù)組作為第二個(gè)參數(shù)傳遞給 ufunc 函數(shù)。例如：:

>>> 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 函數(shù)的第二個(gè)參數(shù)也是 x ，那么它所做的事情就是對(duì) x 中的每給值求正弦值，并且把結(jié)果保存到x 中的對(duì)應(yīng)的位置中。此時(shí)函數(shù)的返回值仍然是整個(gè)計(jì)算的結(jié)果，只不過(guò)它就是 x ，因此兩個(gè)變量的 id 是相同的(變量 t 和變量 x 指向同一塊內(nèi)存區(qū)域)。

比較了一下 numpy.math 和 Python標(biāo)準(zhǔn)庫(kù)的 math.sin 的計(jì)算速度：:

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

計(jì)算100萬(wàn)次正弦值，numpy.sin 比 math.sin 快10倍多。這得利于 numpy.sin在 C 語(yǔ)言級(jí)別的循環(huán)計(jì)算。numpy.sin 同樣也支持對(duì)單個(gè)數(shù)值求正弦，例如：numpy.sin(0.5)。不過(guò)值得注意的是，對(duì)單個(gè)數(shù)的計(jì)算 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

請(qǐng)注意 numpy.sin 的計(jì)算速度只有 math.sin 的1/5。這是因?yàn)?numpy.sin 為了同時(shí)支持?jǐn)?shù)組和單個(gè)值的計(jì)算，其 C 語(yǔ)言的內(nèi)部實(shí)現(xiàn)要比math.sin復(fù)雜很多，同樣在 Python 級(jí)別進(jìn)行循環(huán)的話，就會(huì)看出其中的差別了。此外，numpy.sin返回的數(shù)的類型和math.sin返回的類型有所不同，math.sin返回的是Python的標(biāo)準(zhǔn)float類型，而numpy.sin則返回一個(gè)numpy.float64類型：

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

numpy 中有眾多的 ufunc 函數(shù)提供各式各樣的計(jì)算。除了sin這種單輸入函數(shù)之外，還有許多多個(gè)輸入的函數(shù)，add函數(shù)就是一個(gè)最常用的例子。栗子:

>>> 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 函數(shù)返回一個(gè)新的數(shù)組，此數(shù)組的每個(gè)元素都為兩個(gè)參數(shù)數(shù)組的對(duì)應(yīng)元素之和。它接受第3個(gè)參數(shù)指定計(jì)算結(jié)果所要寫入的數(shù)組，如果指定的話，add 函數(shù)就不再產(chǎn)生新的數(shù)組。

由于 Python 的操作符重載功能，計(jì)算兩個(gè)數(shù)組相加可以簡(jiǎn)單地寫為 a+b，而np.add(a,b,a) 則可以用 a+=b 來(lái)表示。下面是數(shù)組的運(yùn)算符和其對(duì)應(yīng)的 ufunc 函數(shù)的一個(gè)列表，注意除號(hào)"/"的意義根據(jù)是否激活__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]), 如果兩個(gè)數(shù)組的元素為整數(shù)，那么用整數(shù)除法 y = x1 / x2: true divide (x1, x2 [, y]), 總是返回精確的商 y = x1 // x2: floor divide (x1, x2 [, y]), 總是對(duì)返回值取整 y = -x: negative(x [,y]) y = x1**x2: power(x1, x2 [, y]) y = x1 % x2: remainder(x1, x2 [, y]), mod(x1, x2, [, y])

數(shù)組對(duì)象支持這些操作符，極大地簡(jiǎn)化了算式的編寫，不過(guò)要注意如果你的算式很復(fù)雜，并且要運(yùn)算的數(shù)組很大的話，會(huì)因?yàn)楫a(chǎn)生大量的中間結(jié)果而降低程序的運(yùn)算效率。例如：假設(shè)a b c三個(gè)數(shù)組采用算式x=a*b+c計(jì)算，那么它相當(dāng)于:

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

也就是說(shuō)需要產(chǎn)生一個(gè)數(shù)組t保存乘法的計(jì)算結(jié)果，然后再產(chǎn)生最后的結(jié)果數(shù)組x。可以通過(guò)手工將一個(gè)算式分解為 x = a*b; x += c，以減少一次內(nèi)存分配。

通過(guò)組合標(biāo)準(zhǔn)的 ufunc 函數(shù)的調(diào)用，可以實(shí)現(xiàn)各種算式的數(shù)組計(jì)算。不過(guò)有些時(shí)候這種算式不易編寫，而針對(duì)每個(gè)元素的計(jì)算函數(shù)卻很容易用Python實(shí)現(xiàn)，這時(shí)可以用frompyfunc函數(shù)將一個(gè)計(jì)算單個(gè)元素的函數(shù)轉(zhuǎn)換成ufunc函數(shù)。這樣就可以方便地用所產(chǎn)生的ufunc函數(shù)對(duì)數(shù)組進(jìn)行計(jì)算了。讓我們來(lái)看一個(gè)例子。

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

很容易根據(jù)上圖所示寫出如下的計(jì)算三角波某點(diǎn)y坐標(biāo)的函數(shù):

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

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

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

這種做法每次都需要使用列表包容語(yǔ)法調(diào)用函數(shù)，對(duì)于多維數(shù)組是很麻煩的。frompyfunc 栗子：

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

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

def triangle_func(c, c0, hc):def trifunc(x):x = x - int(x) # 三角波的周期為1，因此只取x坐標(biāo)的小數(shù)部分進(jìn)行計(jì)算if x >= c: r = 0.0elif x < c0: r = x / c0 * hcelse: r = (c-x) / (c-c0) * hcreturn r# 用trifunc函數(shù)創(chuàng)建一個(gè)ufunc函數(shù)，可以直接對(duì)數(shù)組進(jìn)行計(jì)算, 不過(guò)通過(guò)此函數(shù)# 計(jì)算得到的是一個(gè)Object數(shù)組，需要進(jìn)行類型轉(zhuǎn)換return np.frompyfunc(trifunc, 1, 1)y2 = triangle_func(0.6, 0.4, 1.0)(x)

通過(guò)函數(shù)triangle_func包裝三角波的三個(gè)參數(shù)，在其內(nèi)部定義一個(gè)計(jì)算三角波的函數(shù)trifunc，trifunc函數(shù)在調(diào)用時(shí)會(huì)采用triangle_func的參數(shù)進(jìn)行計(jì)算。最后triangle_func返回用frompyfunc轉(zhuǎn)換結(jié)果。

值得注意的是用 frompyfunc 得到的函數(shù)計(jì)算出的數(shù)組元素的類型為object，因?yàn)閒rompyfunc函數(shù)無(wú)法保證Python函數(shù)返回的數(shù)據(jù)類型都完全一致。因此還需要再次y2.astype(np.float64)將其轉(zhuǎn)換為雙精度浮點(diǎn)數(shù)組。

廣播

當(dāng)使用 ufunc 函數(shù)對(duì)兩個(gè)數(shù)組進(jìn)行計(jì)算時(shí)，ufunc 函數(shù)會(huì)對(duì)這兩個(gè)數(shù)組的對(duì)應(yīng)元素進(jìn)行計(jì)算，因此它要求這兩個(gè)數(shù)組有相同的大小(shape相同)。如果兩個(gè)數(shù)組的shape不同的話，會(huì)進(jìn)行如下的廣播(broadcasting)處理：

讓所有輸入數(shù)組都向其中shape最長(zhǎng)的數(shù)組看齊，shape 中不足的部分都通過(guò)在前面加1補(bǔ)齊
輸出數(shù)組的shape是輸入數(shù)組shape的各個(gè)軸上的最大值
如果輸入數(shù)組的某個(gè)軸和輸出數(shù)組的對(duì)應(yīng)軸的長(zhǎng)度相同或者其長(zhǎng)度為1時(shí)，這個(gè)數(shù)組能夠用來(lái)計(jì)算，否則出錯(cuò)
當(dāng)輸入數(shù)組的某個(gè)軸的長(zhǎng)度為1時(shí)，沿著此軸運(yùn)算時(shí)都用此軸上的第一組值
上述4條規(guī)則理解起來(lái)可能比較費(fèi)勁，讓我們來(lái)看一個(gè)實(shí)際的例子。

先創(chuàng)建一個(gè)二維數(shù)組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)

再創(chuàng)建一維數(shù)組b，其shape為(5,)：

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

計(jì)算 a 和 b 的和，得到一個(gè)加法表，它相當(dāng)于計(jì)算 a,b 中所有元素組的和，得到一個(gè)shape為 (6,5)的數(shù)組：

>>> 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長(zhǎng)度(也就是ndim屬性)不同，根據(jù)規(guī)則1，需要讓b的shape向 a 對(duì)齊，于是將b的shape前面加1，補(bǔ)齊為(1,5)。相當(dāng)于做了如下計(jì)算：

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

這樣加法運(yùn)算的兩個(gè)輸入數(shù)組的 shape 分別為 (6,1) 和 (1,5) ，根據(jù)規(guī)則2，輸出數(shù)組的各個(gè)軸的長(zhǎng)度為輸入數(shù)組各個(gè)軸上的長(zhǎng)度的最大值，可知輸出數(shù)組的shape為 (6,5)。

由于 b 的第0軸上的長(zhǎng)度為1，而 a 的第0軸上的長(zhǎng)度為6，因此為了讓它們?cè)诘?軸上能夠相加，需要將 b 在第0軸上的長(zhǎng)度擴(kuò)展為6，這相當(dāng)于：

>>> 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軸的長(zhǎng)度為1，而b的第一軸長(zhǎng)度為5，因此為了讓它們?cè)诘?軸上能夠相加，需要將a在第1軸上的長(zhǎng)度擴(kuò)展為5，這相當(dāng)于：

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

經(jīng)過(guò)上述處理之后，a 和 b 就可以按對(duì)應(yīng)元素進(jìn)行相加運(yùn)算了。

當(dāng)然，numpy 在執(zhí)行 a+b 運(yùn)算時(shí)，其內(nèi)部并不會(huì)真正將長(zhǎng)度為1的軸用repeat函數(shù)進(jìn)行擴(kuò)展，如果這樣做的話就太浪費(fèi)空間了。

由于這種廣播計(jì)算很常用，因此 numpy 提供了一個(gè)快速產(chǎn)生如上面 a,b 數(shù)組的方法： ogrid對(duì)象：

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

ogrid 是一個(gè)很有趣的對(duì)象，它像一個(gè)多維數(shù)組一樣，用切片組元作為下標(biāo)進(jìn)行存取，返回的是一組可以用來(lái)廣播計(jì)算的數(shù)組。其切片下標(biāo)有兩種形式：

開始值:結(jié)束值:步長(zhǎng)，和 np.arange (開始值, 結(jié)束值, 步長(zhǎng)) 類似

開始值:結(jié)束值:長(zhǎng)度 j，當(dāng)?shù)谌齻€(gè)參數(shù)為虛數(shù)時(shí)，它表示返回的數(shù)組的長(zhǎng)度，和 np.linspace(開始值, 結(jié)束值, 長(zhǎng)度) 類似：

>>> 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 為什么不是函數(shù)

根據(jù) Python 的語(yǔ)法，只有在中括號(hào)中才能使用用冒號(hào)隔開的切片語(yǔ)法，如果 ogrid 是函數(shù)的話，那么這些切片必須使用 slice 函數(shù)創(chuàng)建，這顯然會(huì)增加代碼的長(zhǎng)度。

利用 ogrid 的返回值，我能很容易計(jì)算 x, y 網(wǎng)格面上各點(diǎn)的值，或者 x, y, z 網(wǎng)格體上各點(diǎn)的值。下面是繪制三維曲面 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庫(kù)快速繪制如下圖所示的3D曲面。

ufunc的方法

ufunc 函數(shù)本身還有些方法，這些方法只對(duì)兩個(gè)輸入一個(gè)輸出的 ufunc 函數(shù)有效，其它的ufunc對(duì)象調(diào)用這些方法時(shí)會(huì)拋出ValueError異常。

reduce 方法和Python的reduce函數(shù)類似，它沿著axis軸對(duì)array進(jìn)行操作，相當(dāng)于將<op>運(yùn)算符插入到沿axis軸的所有子數(shù)組或者元素當(dāng)中。

<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方法類似，只是它返回的數(shù)組和輸入的數(shù)組的shape相同，保存所有的中間計(jì)算結(jié)果：>>> 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 方法計(jì)算多組reduce的結(jié)果，通過(guò)indices參數(shù)指定一系列reduce的起始和終了位置。reduceat的計(jì)算有些特別，讓我們通過(guò)一個(gè)例子來(lái)解釋一下：

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

對(duì)于indices中的每個(gè)元素都會(huì)調(diào)用reduce函數(shù)計(jì)算出一個(gè)值來(lái)，因此最終計(jì)算結(jié)果的長(zhǎng)度和indices的長(zhǎng)度相同。結(jié)果result數(shù)組中除最后一個(gè)元素之外，都按照如下計(jì)算得出：

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

而最后一個(gè)元素如下計(jì)算:

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

因此上面例子中，結(jié)果的每個(gè)元素如下計(jì)算而得：

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) 方法的計(jì)算等同于如下程序：

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

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

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

可以看出通過(guò)outer方法計(jì)算的結(jié)果是如下的乘法表：

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

轉(zhuǎn)載于:https://www.cnblogs.com/oneTOinf/p/10502718.html

與50位技術(shù)專家面對(duì)面20年技術(shù)見證，附贈(zèng)技術(shù)全景圖

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。