[心得] Scipy+Numpy,在Python裡面使用C語言

[painkiller (肚子餓~)]

最近開始用python進行一些數據分析...
想跟各位分享一下 使用Scipy + Numpy 將 C語言整合至python的心得:

python在某些情況下的速度實在是讓人無法忍受
(ex. 對大型List做兩三層迴圈以上的數學運算)

為了增強python在這方面的功能
有兩個好用的module: SciPy & Numpy

numpy 基本上是使用自訂的 "array"
相較於python list, numpy array裡的元素資料型態是一致的
所以在記憶體的使用上有較高效率
而且內建的 "universal function" (ufunc)功能
可以一次對整個numpy array做簡單的數學函數運算
可惜如果運算稍微複雜一點時,
numpy 提供的 frompython (將自訂函數轉成universal function)
相較於使用python本身去寫並沒有快多少
速度上還是不如C語言

這個時候Scipy就可以出場了...
Scipy本身的功能也相當強大,比如說scipy.matplotlib用來繪圖就蠻方便的
其中的scipy.weave更提供了直接在python 裡面使用C語言(and C++)的功能:

舉例來說, hello world:

from scipy import weave
code="""printf("Hello World! \\n");"""
weave.inline(code)

當然, 這看起來像脫褲子放屁 XD

不過如果有一個或好幾個很長的numpy array,
就可以用weave.blitz
結合weave.inline 直接傳到在python裡面寫的C code裡

舉例:
import numpy as np
from scipy import weave
a=np.arange(0,1000,dtype=float).reshape(500,2)
#產生((0,1),(2,3),(4,5)...(998,999))的500*2矩陣,並指定資料型態為float

code="""
        float k=0.0;
        for(int i=0;i<500;i++)
                {
                k+=a(i,0)+a(i,1); //直接在code裡面使用剛剛產生的矩陣
                }
        return Py_BuildValue("d", k); //傳回一個float給python
"""

print weave.inline(code,
                   ['a'],
                   type_converters=weave.converters.blitz,
                   compiler='gcc')
就會出現從0 加到 999的值囉


如果要傳回python list,
可以在C code裡宣告
ex.
PyObject *m = PyList_New(0);
對這個list 可以像在python 裡面一樣操作,
如: PyList_Append(m, PyFloat_FromDouble(k));
因為是在C語言裡面,所以資料型態要指定好

關於如何增加python的效率, 有興趣的人可以看看這邊
http://www.scipy.org/PerformancePython

個人是覺得簡單但是大量的運算
用weave.inline 搭配 numpy, weave.blitz是相當方便的組合
cython or SWIG速度比inline快上一點, 但是使用上稍微複雜一些(對我來說 XD)

如果是使用windows的話 python(x,y)裝了就可以直接用scipy.weave & numpy了
用linux的話就是手動裝scipy & numpy,
要使用weave的話得多裝python-dev
還有要更改這個
/usr/lib/python2.6/dist-packages/scipy/weave/blitz/blitz/blitz.h
在#include <stdio.h> 之前添加：
#include <cstdlib>

以上拋磚引玉
希望有更多高手可以分享使用python 做科學研究的心得 ~OTZ

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 129.170.26.186
※ 編輯: painkiller      來自: 129.170.26.186       (04/09 08:20)
※ 編輯: painkiller      來自: 129.170.26.186       (04/09 08:21)
→ juiz:Join us in http://conference.scipy.org/scipy2011/
推 pcedison:感謝分享
推 curist:推推
推 mjsg:感謝分享
推 cobrasgo:如果是我的話，重點應該會放在那兩三層的迴圈有沒有辦法
→ cobrasgo:改寫邏輯
→ painkiller:嗯...大概是我舉得例不好 光是取sin值100萬次的話
→ painkiller:python就比weave慢個10倍左右
推 guteres:請問有沒有辦法把code另外存成檔案再去call它?
→ painkiller:可能要用SWIG之類的套件喔... 那個我很不熟耶 :P
推 guteres:soga
推 cobrasgo:用c binding不就可以了？