Yukuan's Blog: scientific

Showing posts with label scientific. Show all posts

Sunday, September 24, 2006

Dynamic Plot on Python

之前介紹 Python 在複雜網路模擬實驗的應用時，曾提到用 matplotlib 來繪製該實驗的圖表。其模仿 Matlab 的繪圖功能，用起來方便，所繪製的圖也在水準之上。有圖有真相，這就秀秀用 matplotlib 為該實驗繪製的兩張圖表：

學長看了後提議：某些實驗邊跑模擬的同時，也把產生的數據動態繪出，效果會更好。

憑著殘留印象，我把相關用法的說明再次翻出來研讀，寫了下列的 Python code 來測試：

"""
Interactive random dot plotting.
 
This module demonstrates the usage of the interactive mode of matplotlib.
It is suggested to run this code under ipython. Details please see
http://matplotlib.sourceforge.net/interactive.html
 
 
Usage
=====
 
You can run this module in three ways:
 
    1. As a stand-alone application with GUI
    2. Run it on command line with a -d parameter to just generate the
       image with a specific format. (without GUI)
        * ref. http://matplotlib.sourceforge.net/backends.html
    3. Run it on ipython to call loop() many times
 
Examples
--------
You can simply run this module as a stand-alone application with GUI
 
    > irandot.py
 
Or just generate a image file with a specific format:
 
    > irandot.py -dAGG     # generate PNG file
    > irandot.py -dPS      # generate PS file
    > irandot.py -dSVG     # generate SVG file
    > irandot.py -dEMF     # generate EMF file
 
Running it with ipython is strongly suggested.
 
    In [1]: cd I:\trial\python  # suppose the module is put in I:\trial\python
    ...
    In [2]: run irandot.py
    ...
    In [3]: loop(50)
    ...
    In [4]: loop(20)
    ...
 
"""
__author__ = "Jiang Yu-Kuan, yukuan.jiang(at)gmail.com"
__date__ = "2006/09/24~2006/9/25"
__revision__ = "1.3"
 
from pylab import *
 
 
def dot_gen(dots=2):
    """A random dot generator.
 
    - dots: number of dots with plotting each loop
    """
    i = 0
    while 1:
        l, = plot(rand(dots), rand(dots), 'go')
        #setp(l, alpha=1,
        #    markerfacecolor='w', markeredgecolor='g', markersize=10)
        i += dots
        yield i
 
 
def loop(n=1):
    """Repeat dot plotting of n times."""
    g= dot_gen()  # gets a dot-plotting generator
    for i in xrange(n):
        print g.next(),
 
 
if __name__ == '__main__':
    if not isinteractive():
        ion()  # turns interactive on
 
    loop(50)
    savefig('irandot')
 
    try:
        __IPYTHON__
    except NameError:
        print "\nNOT on IPython"
        show()
    else:
        print "\nOn IPython"

行 59 是真正畫綠點的指令。
行 75 是關鍵所在，它切換 matplotlib 到互動模式。
行 80~86 是為了在 IPython 下不要執行 show()

這個程式取名為 irandot.py ，可以當作一般的執行檔，點兩下檔名執行。執行時我們可以看到隨機產生的座標，一點一點地在座標圖上冒出。以下是其 screenshot：

此外，強烈建議在 IPython shell 下執行 irandot.py：

開啟 IPython shell
在 IPython shell 敲進 "cd I:/trial/python" 來切換到存放 irandot.py 的目錄。
鍵入 "run irandot" 後就會繪出類似上面那幅 2*50=100 個綠點的隨機圖。
再鍵入 "loop(100)" ，程式會陸續再冒出 2*100 = 200 個綠點。
高興的話，可以一直下 loop 、改變參數，讓綠點繼續冒下去。

後記：

之前發現無論在 IDLE 或 PythonWin 下，圖表反覆繪製幾次，程式就會掛點的問題。碰巧藉這次測試，找到解法了。作法很簡單，只要改用 IPython shell 即可。
如果你也是裝 Python Enthought 來用，那 IPython 已經在你電腦了。
目前的經驗，開發過程使用 IDLE 比 PythonWin 順手，測試寫好的程式則用 IPython 較方便
如果嫌上述開發環境太陽春，我推薦 SPE 。
Python Enthought 附的 matplotlib 有些小 bug 。可以到這裡下載新版 matplotlib 安裝程式來更新。其底層還用到 NumPy ，可到這裡下載安裝。

Tuesday, August 29, 2006

Graph-based Modeling on Python

Agent-based modeling 的電腦實驗，最核心的架構不外乎一個大迴圈（super loop）包裹著一群規則。大迴圈每跑一輪，系統就更新一次狀態，就如同時鐘的滴答（tick）聲般。通常系統每次滴答都會收集一次統計資料。這類實驗，有許多現成的 famework 可用，如最經典的 Swarm 及其後進 Repast ，還有我模仿 Repast ，自己搞的一個，它們都提供了 start, pause, stop 等流程控制的介面。

模擬複雜網路，也可以套用 Agent-based modeling 架構。不過諸如網路的群聚度（clustering coefficient, C）及網路特徵的路經長度（characteristic path length, L）等統計數據計算需耗費的時間，隨著網路的規模成長很快，所以不適合運作太頻繁。但我們卻得靠這些統計數據來判斷網路是否收斂、試驗是否該終止了，想試驗各個參數排列組合時，更是雪上加霜。

透過 Intermediate File 是很直覺的解法。程式每次執行都餵入一個輸入檔來初始網路結構，並以參數決定 super loop 要跑幾次。統計數據只在程式要結束時計算，然後將網路狀態及統計數據吐給輸出檔。下次要執行，就以這個輸出檔當作新的輸入，然後一樣決定要跑幾圈，最後得到另一個輸出檔。採取接力的方式，不用每次都從頭跑，不會浪費先前程式執行的時間。

Multithread 是另一個可行作法。採用 Agent-based modeling framework 的方式，一個 thread 跑核心的規則。另一個 thread 則控制試驗的流程，讓 user 決定何時該暫停下來，計算統計數據，然後決定是否繼續執行。

Dynamic Programming Language 是這次的正解。利用動態程式語言，如 Python ，可以很方便於執行時期控制模擬試驗的程式流程。我就是採用 Python ，搭配 NetworkX 來處理 Graph ，並以 matplotlib 來繪製圖表。接下來就看看如何架設這個試驗平台：

下載 Python Enthought Edition 的 installer 。它是 Python 的加強版，除了 Python 基本工具外，還整合了許多有用的工具，例如 matplotlib 。
到這裡下載 NetworkX 的 installer 。
依序執行這兩個 installer 。

裝好了後，就以 Emergence of a small world from local interactions: modeling acquaintance networks 裡描述的實驗來操刀，試寫第一個 Python 程式：

"""
Simulation for the verifying of "Emergence of a Small World from Local
Interactions: Modeling Acquaintance Network", Davidsen J. et al. 2002.
 
Public Variables:
    - Iters: the interations of loop, used in loop()
 
Private Variables:
    - _G: graph of the acquaintance network
    - _N: the total of nodes
    - _p: the death-and-birth probability, used in _step2()
    - _run: a generator returned by loop()
 
Usage
=====
 
This module must run on a python shell with interative mode; and call
public functions e.g. runnext().
 
Example
-------
You can simply run this module with default settings.
 
    > python -i acq2002.py
    >>> runnext(10)         # loops 10 times of Iters time-steps.
    ...
 
Setting parameters is allowed.
 
    >>> start(1000, 0.04)   # starts the simulation with new N, and p
    >>> Iters=10000         # sets the iterations of each loop
    >>> runnext(5)          # loops 5 times of 10000 time-steps.
    ...
    >>> Iters=100000        # sets the iterations of each loop
    >>> runnext(2)          # loops another 2 times of 100000 time-steps
    ...
 
After above, you may want to save the result.
 
    >>> savepkfig()
    >>> savenet()
 
"""
__author__ = "Jiang Yu-Kuan, yukuan.jiang(at)gmail.com"
__date__ = "2006/08/22~2006/10/02"
__revision__ = "1.10"
 
import random as rnd
import networkx as nx # Must import this as a name to avoid namespace collision!
 
 
def _create():
    """Create an empty acquaintance network.
 
    Adds _N nodes to _G, and _G contains no any edges.
    """
    global _G
    _G = nx.Graph()
    _G.add_nodes_from(xrange(_N))  # adds n nodes into _G
 
 
def _step1():
    """Friend making of two persons.
 
    One randomly chosen person picks any two his acquaintances and
    introduces them to each another. If they have not met before, a new
    link between them is formed. In case the person chosen has less than
    two acquaintances, he introduces himself to one other random person.
    """
    u, v = rnd.sample(xrange(_N), 2)
    nb = _G.neighbors(u)
    if len(nb) > 1:
        u, v = rnd.sample(nb, 2)
    _G.add_edge(u, v)
 
 
def _step2(p):
    """Death and Birth of a chosen person.
 
    With probability p, one randomly chosen person is removed from the
    network, including all links connected to this node, and replaced by
    a new person with one randomly chosen acquaintance.
    """
    if rnd.random() < p:
        # sol. 1
        v = rnd.randrange(_N)
        _G.delete_node(v)
        _G.add_node(v)
 
        # sol. 2: slower than sol. 1.
        #_G.delete_edges_from(_G.edges(rnd.randrange(_N)))
 
 
def _avgK2():
    """Return the average square of degree of all nodes."""
    kss = [k*k for k in _G.degree(_G.nodes())]
    return sum(kss)/float(_N)
 
 
def _avgSPL():
    """Return the average shortest path length of Graph _G."""
    pathlengths=[]
    for v in _G.nodes():
        # sol. 1
        #spl = nx.single_source_shortest_path_length(_G,v)  # including v to v
        #for pl in spl.values():
        #    pathlengths.append(pl)
 
        # sol. 2
        pathlengths += nx.single_source_shortest_path_length(_G,v).values()
 
    return sum(pathlengths) / float(len(pathlengths)-_N)
 
 
def _stat():
    """Gather statistics for Graph _G."""
    print "N, p = %d, %f" % (_N, _p)
    print "Degree histogram:", nx.degree_histogram(_G)
    print "Avg. degree, <k>:", 2.*_G.number_of_edges()/_N
    print "Avg. square of degree, <k^2>:", _avgK2()
    print "Avg. clustering coefficient, C:", nx.average_clustering(_G)
    print "Avg. shortest path length, L:", _avgSPL()  # spends huge time
 
 
def savepkfig(fn='acq2002pk.png'):
    """Save the p(k) figure.
 
    - fn: file name
    """
    import pylab as mpl
    dh = nx.degree_histogram(_G)
    pk = mpl.array(dh)/float(_N)
 
    mpl.loglog(pk, 'r--')
    mpl.grid(True)
    mpl.gca().xaxis.grid(True, which='minor')
    mpl.xlabel('k')
    mpl.ylabel('p(k)')
    mpl.title('N=%d, p=%d' % (_N, _p) )
    mpl.savefig(fn)
 
 
def savenet(fn='acq2002', fm='GPickle'):
    """Save the acquaintance network.
 
    - fn: file name
    - fm: file format
    """
    save = {'edgelist':nx.write_edgelist,   # edge list
            'adjlist':nx.write_adjlist,     # node adjacency-list
            'yaml':nx.write_yaml,           # YAML text format
            'gpickle':nx.write_gpickle      # Python pickle format.
            }
    fm = fm.lower()
    if fm not in save.keys():
        fm = 'edgelist'
 
    import os
    fn, ext = os.path.splitext(fn)
    ext = ext[1:].lower()
    if ext in save.keys():
        fm = ext
    save[fm](_G, '.'.join([fn, fm]))
 
 
def _loop():
    """Loop _step1 and _step2 of Iters times and gather statistics.
 
    This function uses yield statement and returns the total iterations of loop
    """
    import time
    t_start = time.time()  # records the start time
    t_ttl = 0  # clears the total time spent
    i_cnt = i_ttl = 0  # clears the loop counter and total loops.
 
    _create()  # create an empty network
 
    while 1:
        i_cnt += 1
 
        _step1()
        _step2(_p)
 
        if i_cnt == Iters:
            i_ttl += i_cnt
            i_cnt = 0
 
            _stat()
 
            t_inc = time.time() - t_start
            t_ttl += t_inc
            print "Time spent (increment,total): (%f,%f)" % (t_inc, t_ttl)
            yield i_ttl
            t_start = time.time()
 
 
def runnext(n=1):
    """Call _run.next() of n times."""
    for i in xrange(n):
        print _run.next()
 
 
def start(N=100, p=0.04):
    """Start the module.
 
    - N: the total of nodes of the network
    - p: the death-and-birth probability
    """
    global _N, _p, _run
    _N, _p = N, p
    _run = _loop()  # gets a generator
 
 
Iters = 100000
 
if __name__ == "__main__":
    # Uses Psyco if available
    try:
        import psyco
        psyco.full()
        print "Psyco: OK"
    except ImportError:
        print "Psyco: FAIL"
 
    start(100, .04)
 
    """
    import profile, pstats
    profile.run('runnext()', 'acq2002.prof')
 
    s = pstats.Stats('acq2002.prof')
    s.sort_stats('time','name').print_stats()
    """

行 203~211 的 start() 除了測試這個模組外，也用作主程式。
行 211 讓我們取得 generator 給 _run 。
行 214 設定每次 _run.next() 要跑幾輪 super loop ，每 _run.next() 一次，才統計一次數據。
行 225 執行 start() ，並指定不同的 node 數 N 及死生發生的機率 p 。
這個模組執行後，可以在 Python shell 下 _run.next() ，每 _run.next() 一次，大迴圈就跑 Iters 輪。此外，還可以在 Python shell 中改變 Iters 的次數。
行 197~200 的 runnext() 讓 Python shell 下，調用多次 _run.next() 更方便。
行 217~223 使用 Psyco 來加速。

註：

NetworkX 的用法可以參考其 Tutorial 及 API Reference 。
Python 的 random 模組說明可以參考這裡。
要以 matplotlib 繪製圖表，可以參考其 Tutorial 及 Screenshots 。其刻意模仿 Matlab 的用法，用來倍感親切。
其他 Python 的說明文件可以參考 Python Tutorials 。

Saturday, January 07, 2006

Scientific Computation Tools

作科學研究常需要作些計算，可以幫我們作這些運算的軟體很多，其中最讓我印象深刻的是 Mathematica 及 Matlab 。它們用起來方便，功能更是不在話下，但就是要花錢買，而且還不便宜。所以我就找了一些功能雖沒那麼完整，但還算堪用的軟體來試試，重點是都是免費的。例如用於代數、符號運算的 Axiom, Maxima, Yacas 等。方便用於數值、矩陣運算的 FreeMat, Octave, RLaBPlus, SciLab 等。

Mathematica 是一套電腦代數系統（computer algebra system, CAS）。我們可以利用它來作因式分解，及微積分運算，值得注意的是，這些運算單元可以用代數符號表示就好，不用非得要是數值，因為它有精湛的符號運算能力。例如：我們可以要求它作 a^2 - b^2 的因式分解，而得到 (a-b) * (a+b) 這個答案。

Matlab 則對數值運算比較專門。其名字就是 Matrices Lab 的簡稱。可見它擅長在大量的矩陣運算場合發揮。我們把這類的軟體稱為數值運算環境（numerical computing environment, NCE）。

免費的 CAS ，這裡推薦 Maxima ，安裝之前，可以先到這裡試用一下，為了用起來更賞心悅目，記得順便把 wxMaxima 裝上。如果你想玩玩 Java 版的 CAS ，那就試試 Yacas ，介面雖陽春了些，但語法接近 C 語言，可能會讓你感到親切吧，可以偷喵一下 Yacas 的使用範例。此外 Axiom 也是不錯的選擇（它的 Windows 版本沒有繪圖模式，比較難看），你也可以先到 Axiom wiki's Sandbox 試用看看，再決定要不要裝來用。 Maxima 及 Yacas 的 license 是 GPL 的； Axiom 則是 Modified BSD License。

在逛 Axiom 網站時，還意外發現了 GNU TeXmacs 這個 wysiwyw 的科學文件編輯平台，可惜其 Windows 版的 porting 不大穩定，且無法正常處理中文。也許可以改而安裝其掛在 Cygwin 下的版本，有興趣的話可以試試。

要作數值及矩陣方面的科學計算，前述的 CAS 為主的軟體雖然也有支援部份，但總不是那麼方便。所以我另外找了提供專門數值運算環境（numerical computation enviroment）的免費軟體。 Octave 及 SciLab 是其中較有名： Octave 雖然是 GPL 授權的，但其 Windows 版要跟 Cygwin 綁在一起執行，不是那麼方便； SciLab 雖個人使用也是免費的，但它有專屬的授權方式，在商業用途有些限制，且說明文件的安排不是那麼好。我還另外試了一套叫 FreeMat 的，它是 MIT Type licence 的，說明文件淺顯易懂， Windows 的支援也不錯。另外還有 RLaBPlus 這套 RLaB 的後繼者，其 script language 的語言正交性較好，很有發展潛力，目前只有 Linux 版本。