1. 程式人生 > 程式設計 >python科學計算之numpy——ufunc函式用法

python科學計算之numpy——ufunc函式用法

寫在前面

ufunc是universal function的縮寫,意思是這些函式能夠作用於narray物件的每一個元素上,而不是針對narray物件操作,numpy提供了大量的ufunc的函式。這些函式在對narray進行運算的速度比使用迴圈或者列表推導式要快很多,但請注意,在對單個數值進行運算時,python提供的運算要比numpy效率高。

四則運算

numpy提供的四則ufunc有如下一些:

numpy提供的四則運算unfunc能夠大大的提高計算效率,但如果運算式複雜,且參與運算的narray過大,會產生大量的中間結果,從而降低計算效率。例如:計算x=a*b+c時,實際上會按照如下方式計算:

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

這會產生兩次記憶體分配,一次時x,一次時t,所以按照

x = a*b 
x = x+c

會節省一次記憶體的分配,從而提高效率。

比較運算 & bool運算

numpy同時提供了=、<、>等這些比較運算子,這些運算子的結果是bool值或者bool陣列。

np.array([1,2,3]) < np.array([0,3,4])
array([False,True,True],dtype=bool)

邏輯運算and、or、not、xor等由於python提供了,所以numpy的這些邏輯運算子是以logical_開頭的函式

a = np.arange(5)
b= np.arange(4,-1,-1)
## print a==b
[False False True False False]
## print a>b
[False False False True True]
## print np.logical_or(a == b,a > b)
[False False True True True]

對兩個bool陣列進行邏輯運算時,將發生ValueError異常,因為bool值也是True和False,numpy無法確定運算目的,可以使用numpy.any()和numpy.all()函式,他們的使用方法和python的any()、all()函式用法相同。以bitwise_開頭的函式時用於位運算,如(bitwise_and、bitwise_or)等,也可以使用&、|、~和^來進行運算。

除了numpy提供的內建ufunc函式,使用者也可以編寫自定義的ufunc函式,方式是:

1. 編寫對單個數值計算的目的函式;

2. 利用np.frompyfunc(func,nin,nout)將其轉換為ufunc函式,其中func是上面編寫的目的函式,nin是輸入的引數個數,nout是返回值的個數。

## 基本形式
u_func = np.frompyfunc(func,nout)
ret = u_func(narray_obj,param1,param2..)

這裡返回的ret是object型別,所以實際上需要用astype()轉換為目的型別。numpy.vectorize()也實現了和numpy.frompyfunc()一樣的功能,區別是前者可以t通過otypes指定返回值的型別,不用再用astype()進行轉換。

## 基本形式
u_func = np.frompyfunc(func,otypes=[dtype1,dtype2..]
ret = u_func(narray_object,param2..)

廣播

先看個例子:

a = np.arange(0,60,10).reshape(-1,1)
b = np.arange(0,5)
#a
array([[ 0],[10],[20],[30],[40],[50]])
#b
array([0,1,4])

ok,現在計算a+b,不過現在有一個問題,a和b的維度不一樣,那應該怎麼加?先看看結果吧

# a+b
array([[ 0,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]])

結果來看,是用a的每一行的元素(這裡a為列向量,每一行只有一個元素)與b的每一個元素相加,相當於:

a=array([[ 0,0],10,10],20,20],30,30],40,40],50,50]])

而b是一個行向量,現在我們將這一個行向量重複6次,和a的第0軸長度相同,構成一個二維陣列,相當於:

b=array([[0,[0,4]])

現在,再進行相加,自然就是對用元素相加了,也就是上面的結果,這就是numpy中的廣播,對進行運算的兩個narray物件shape不一樣時,進行的維度補齊。總的來說,numpy的廣播規則基於下面4個規則:

讓所有輸入數姐都向其中維數最多的陣列看齊,shape屬性中不足的部分都通過在前面加1補齊; 如上面的輸入中,a.shape=(6,1),b.shape=(,5),a的維數是2,b的維數是1,所以b向a看齊,並且用1補齊,那麼b.shape=(1,5)。

輸出陣列的shape屬性是輸入陣列的shape屬性的各個軸上的最大值 輸出是各軸上最大值,所以a+b的輸出的shape應該是(6,5);

如果輸入陣列的某個軸的長度 為 1或與輸出陣列的對應軸的長度相同,這個陣列能夠用來計算,否則出錯

當輸入陣列的某個軸的長度為1吋,沿著此軸運算時都用此軸上的第一組值

由於廣播在numpy計算中比較常見,所以numpy提供了ogrid和mgrid來建立廣播計算的陣列,前者返回的是兩個向量,後者返回的是進行廣播運算的陣列。

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

ogrid[]引數是有兩種形式(1):[start:end:step]即起點、終點和步長;(2):[start:end:len]即起點、終點和陣列長度,這裡的長度為了和步長區分,在傳入時時以虛數方式傳入的,例如[1,5,4j]將產生[1,4]這樣的陣列。

另外廣播支援特殊的下標None,意義時在對應位置上產生1的新軸,例如a[None,:]等效於a.reshape((1,-1)),a[:,None]等效於a.reshape((-1,1)),另外,np.ix_()能將一維陣列轉換為能進行廣播的二維陣列:

x=np.array([0,10])
y=np.array([2,8])
gy,gx=np.ix_(y,x)
#print gy
[[2]
[3]
[8]]
#print gx
[[ 0 1 4 10]]
# gx+gy
array([[ 2,6,12],[ 3,7,13],[ 8,9,18]])

np.ix_()支援多個引數,從而支援n維空間的廣播計算。

ufunc方法

ufunc函式物件本身還有一些方法函式,這些方法只對兩個輸入、一個輸出的ufunc 函式有效,其他的ufunc物件呼叫這些方法時會丟擲ValueError異常。

(1). reduce(),沿著指定軸對陣列進行操作,相當於將相應的操作放到該軸元素之間。

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

(2). accumulate()和reduce()類似,區別時是前者會保留中間結果:

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

(3). reduceat()方法計算多紐reduce()的結果,通 過 indices引數指定一系列的起始和終止位置。它的計算有些特別,,計算的方法如下:

if indices[i] < indices[i+1]:
 result[i] = <op>.reduce(a[indices[i]:indices[i+1]])
else:
 result[i] = a[indices[i]]
#result[-1]的計算如下:
<op>.reduce(a[indices[-1]:])

例:

a = np.array([1,4])
result = np.add.reduceat(a,indices=[0,0])
## result
array([1,10])

## 計算過程如下:
 : a[0] -> 1
 : a[1] -> 2
 : a[0] + a[1] -> 1 + 2
 : a[2] -> 3
 : a[0] + a[1] + a[2] -> 1 + 2 + 3 = 6
 : a[3] -> 4
 :a[0] + a[1] + a[2] + a[4] - > 1 + 2 + 3 + 4 = 1 0

再看多維陣列

在前一篇文章中我們提到過多維陣列,現在我們回頭再看看多維陣列的下標取資料。首先,多維陣列的下標應該是一個長度和陣列的維數相同的元組。如慄下標元組的長度比陣列的維數大,就會出錯;如果小,就 會 在 下 標 元 組 的 後 而 補 ,使得它的長度與陣列維數相同。如果下標物件不是元組,則 numpy會首先把它轉換為元組。這種轉換可能會和使用者所希望的不一致,因此為了避免出現問題,請 “顯式”地使用元組作為下標。

a = np.arange(3*4*5).reshape(3,5)
array([[[ 0,[ 5,8,9],[15,16,17,18,19]],[[20,[25,26,27,28,29],[35,36,37,38,39]],[[40,[45,46,47,48,49],54],[55,56,57,58,59]]])

lidx=[[0],[1]]
#a[lidx]

aidx = np.array(lidx)
#a[aidx]
array([[[[ 0,19]]],[[[20,39]]]])

可以看到,numpy把列表[[0],[1]]轉換成元組([0],[1]),而把陣列物件轉換成(aidx,:,:)。

整數陣列和切片做為下標

如果利用陣列作為下標,但陣列的維度都不一樣,那麼這個時候這些陣列將會應用廣播規則把這些陣列轉換成一樣,轉換的規則是上面說到的先用1補足然後取各個維度的最大值。

i0 = np.array([[1,0]])
i1 = np.array([[[0]],[[1]]])
i2 = np.array([[[2,2]]])
#print i0.shape
(2,3)
#print i1.shape
(2,1)
#print i2.shape
(1,3)
# i0補齊之後為 (1,2,3),然後取最大值為
(2,2,3)

這裡將所有的陣列進行廣播,把所有的陣列轉換成shape=(2,3),利用numpy.broadcast_arrays()可以檢視廣播後的陣列:

id0,id1,id2=np.broadcast_arrays(i0,i1,i2)
#id0
[[[1 2 1]
 [0 1 0]]

 [[1 2 1]
 [0 1 0]]
#id1
[[[0 0 0]
 [0 0 0]]

 [[1 1 1]
 [1 1 1]]]
#id2
[[[2 3 2]
 [2 3 2]]

 [[2 3 2]
 [2 3 2]]]

然後利用下標取元素時,例如i,j,k=1,1時,取得的元素應該是a[id0[i,k],id1[i,id2[i,k]]=a[1,3]=28。下標中含有切片時,首先來看第一種,就是整數陣列時連續的,也就是陣列之間沒有切片,例如a[1:3,i0,i1],這種情況下,首先會把整數陣列(i0,i1)廣播,然後將切片的長度放到相應的位置構成一個維度,i0,i1廣播後的shape=(2,3),切片的長度為2,所以最後的shape=(2,3)。最後的結果是a[1:3,id0[i,k]]。再者,如果切片是再整數陣列之間,那麼同樣會將陣列廣播,然後把切片位置替換為陣列的維度或者切片的長度新增到索引的最後面,如a[i0,i1],i0,i1廣播後的shape=(2,3),陣列a的第二個軸的長度為4,所以最後的shape=(2,4)。

bool陣列做為下標

用bool陣列作為下標時,會將bool陣列中為True的索引組成一個整數陣列,然後以此陣列作為下標。相當於用numpy.nonzero(bool_array)的結果,例如:

b2 = np.array([[True,False,[True,False]])
np.nonzero(b2)
##
(array([0,1]),array([0,0]))
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
## 相當於取[0,0]

除此之外,如果bool陣列中有切片,那麼相當於將bool值轉換成nonzero()後的整數陣列,切片位置不變。

以上這篇python科學計算之numpy——ufunc函式用法就是小編分享給大家的全部內容了,希望能給大家一個參考,也希望大家多多支援我們。