最近一直沒時間寫，但 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 )

好，你看，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]]

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