最近一直沒時間寫,但 numpy 實在值得大書特書。目前 debian unstable 裡已經有了 numpy 1.0rc1 (released at 2006-09-21) 的包裝,看完本文後強烈建議 apt-get 下來看一看。用 Windows 的人可以抓 win32 prebuild binary 或直接抓 enthon

再簡單介紹一下:numpy 是融合了 numeric python (python numerical extension) 和 numarray 能力,並在介面上大幅翻新的 python numerical extension。如果你用過 numeric 和 numarray,可以改換到 numpy 了。全新的 numpy 讓 coding 變成一件很快樂的事。

搞計算的人常需要處理多維陣列。Fortran 90 提供很方便的功能來進行行列操作:

real*8, dimension(:,;,:), allocatable :: ary
allocate( ary(12,11,10) )
ary(:,:,:) = 0.d0
ary(:,1,1) = 1.d0

在 Fortran 77 和 C 程式裡,最後兩行非得寫成迴圈不可。numpy 提供相同的方便性:

from numpy import empty
ary = empty( (10,11,12), dtype=float )
ary[:] = 0.0
ary[0,0,:] = 1.0

對了,後兩行還可以寫成:

ary.fill( 0.0 )
ary[0,0,:].fill( 1.d0 )

Note

很多時候呼叫 fill method 來把純量設裡陣列會比用 assign operator 來得快。

好,你看,numpy 很適合用來處理陣列資料 (笑)。

以上是很基本的功能。在處理多維陣列的時候我們還常會遇到另一個問題:依序處理陣列中不同維度 (軸) 裡的資料。在 Fortran 裡,好像非寫三次三重迴圈不可:

do i = 1, 10
  do j = 1, 11
    do k = 1, 11
      ary(k+1,j,i) = ary(k+1,j,i) + ary(k,j,i)
    end do
  end do
end do

do i = 1, 10
  do k = 1, 12
    do j = 1, 10
      ary(k,j+1,i) = ary(k,j+1,i) + ary(k,j,i)
    end do
  end do
end do

do k = 1, 12
  do j = 1, 11
    do i = 1, 9
      ary(k,j,i+1) = ary(k,j,i+1) + ary(k,j,i)
    end do
  end do
end do

啊,這裡的處理太簡單了,不必用三重迴圈,但同樣的碼還是得加三次:

do k = 1, 11
  ary(k+1,:,:) = ary(k+1,:,:) + ary(k,:,:)
end do
do j = 1, 10
  ary(:,j+1,:) = ary(;,j+1,:) + ary(:,j,:)
end do
do j = 1, 9
  ary(:,:,i+1) = ary(:,:,i+1) + ary(:,:,i)
end do

暫時別讓我對索引值最佳化,煩死人了,連檢查索引值上下限都覺得有點眼花。現在換 numpy:

c0sc = slice(0,None)
p1sc = slice(1,None)
m1sc = slice(0,-1)
for jcrd in range(3):
  sct = [c0sc,c0sc,c0sc]
  sct[jcrd] = p1sc
  p1sct = tuple(sct)
  sct[jcrd] = m1sc
  m1sct = tuple(sct)
  bry = ary[m1sct].copy()
  ary[m1sc] += bry

slice 是 python 內建的函式,它可以建出 python sequence 可接受的 slicing 參數;上面程式的 p1sc 等於是把 "1:" 這種寫法參數化,而 m1sc 則代表 ":-1",c1sc 代表 ":"。因此,善用 python 的 slice 函式,我們可以如上把維度的變化參數化,不必寫三個三重迴圈!

我知道,用上面的例子來看 slice 的用法是太抽象了點。那麼換個實際的陣列來看看:

>>> from numpy import arange
>>> a = arange(16,dtype=int).reshape(4,4)
>>> print a
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]
>>> s1 = slice(1,None); s2 = slice(0,-1)
>>> print a[(s1,s2)]
[[ 4  5  6]
 [ 8  9 10]
 [12 13 14]]
>>> print a[1:,:-1]
[[ 4  5  6]
 [ 8  9 10]
 [12 13 14]]

自由來去,縱橫左右。不錯吧!如果有人肯花時間弄懂我亂寫的程式的話 (笑)。總之謝謝收看。

Posted by yungyuc at 21:35, 0 comment, 0 trackback.
Navigate
Add a trackback
Add a comment

Your name. (required)

Your personal website. (optional)

Your email address. Will not show in page. (suggested, but optional)

Text format is "Plain Text".

Enter "Zthgb"
© hover year to navigate month: powered by django