12012010 Python matplotlib

Beginning Python Visualization

matplotlibがメインなのかな？目次だけだとよくわからん。

ProductName Beginning Python Visualization: Crafting Visual Transformation Scripts (Books for Professionals by Professionals)
Shai Vaingast
Apress / ￥ 4,123 ()
在庫あり。

ビジュアライジング・データ読んどけば良い気もするが。

ProductName ビジュアライジング・データ ―Processingによる情報視覚化手法
Ben Fry
オライリージャパン / ￥ 3,780 ()
在庫あり。

12012010 Python llvm

macbookにllvm-pyを入れてみた

おもむろにllvmをインストール。

sudo port install llvm

でok。llvm-pyはサイトからダウンロードしてきていれた。

で、サンプルをやってみる。

#!/usr/bin/env python

from llvm import *
from llvm.core import *
from llvm.ee import *

my_module = Module.new('my_module')
ty_int = Type.int()   # by default 32 bits
ty_func = Type.function(ty_int, [ty_int, ty_int])
f_sum = my_module.add_function(ty_func, "sum")
f_sum.args[0].name = "a"
f_sum.args[1].name = "b"
bb = f_sum.append_basic_block("entry")
builder = Builder.new(bb)
tmp = builder.add(f_sum.args[0], f_sum.args[1], "tmp")
builder.ret(tmp)

mp = ModuleProvider.new(my_module)

ee = ExecutionEngine.new(mp)

arg1 = GenericValue.int(ty_int, 100)
arg2 = GenericValue.int(ty_int, 42)

print my_module
retval = ee.run_function(f_sum, [arg1, arg2])

print "returned", retval.as_int()

実行結果

; ModuleID = 'my_module'

define i32 @sum(i32 %a, i32 %b) {
entry:
    %tmp = add i32 %a, %b       ; <i32> [#uses=1]
    ret i32 %tmp
}

returned 142

12012010 chemoinformatics perl bioinformatics Python

オライリーからはやっぱケモインフォマティクスの本も出して欲しいと思う

バイオインフォだけだと片手落ちの気がする。表紙はケクレ構造にちなんで亀でお願いしたい。

ProductName バイオインフォマティクスのためのPerl入門
水島 洋
オライリー・ジャパン / ￥ 5,040 ()


ProductName 実践 バイオインフォマティクス -ゲノム研究のためのコンピュータスキル-
Cynthia Gibas,Per Jambeck
オライリー・ジャパン / ￥ 4,410 ()


今のCPUとかmemoryの状況でLLで量子化学とか、かなりアリだと思うんだけどなぁ。pyquanteとかやばいでしょう。

12012010 Python network

禁断探索法をpythonで

ヒューリスティクスは割と好きだ。スパッととけないものをなんとか現実と折り合いを付けながらよりよい解を見つけるみたいな、なんとも煮え切らない美しさのなさがたまらん。インフォマティクスも似たような香りはする。

特に、「シミュレーションできない学問は学問として未成熟である」という言葉に従うのであれば、創薬研究は（学際領域）と言われている割には未成熟な学問の組み合わせのために発見的な手法の割合が増えすぎるし、精度の高い予測法の誕生もまだまだ先であろう。というわけで、発見的な探索手法は当分有用だし、学際領域故に融合部分での応用が期待されるため、アルゴリズムとしてきちんと押さえておくと色々役に立つ。

メタヒューリスティクスの数理の4章はpythonのコードが載っているので、勉強になる。

ProductName メタヒューリスティクスの数理
久保 幹雄,J. P. ペドロソ
共立出版 / ￥ 3,675 ()
在庫あり。

  • 第4章 応用
  • 4.1 グラフ分割問題
  • 4.2 最大安定集合問題
  • 4.3 グラフ彩色問題
  • 4.4 巡回セールスマン問題
  • 4.5 2次割当問題
  • 4.6 多制約ナップサック問題
  • 4.7 数分割問題

4.1のグラフ分割問題をnetworkxで可視化するように書き換えた。

path

import math,random
from sys import maxint as Infinity

def rnd_graph_fast(nnodes, prob):
    nodes = range(nnodes)
    adj = [[] for i in nodes]
    i=1;j=1
    logp = math.log(1.0 - prob)
    while i<nnodes:
        logr = math.log(1.0 - random.random())
        j += 1+int(logr/logp)
        while j>=i and i<nnodes:
            j -= i
            i += 1
        if i<nnodes:
            adj[i].append(j)
            adj[j].append(i)
    return nodes, adj

def construct(nodes):
    sol = [0 for i in nodes]
    for i in range(len(nodes)/2):
        sol[i] = 1
    return sol

def evaluate(nodes, adj, sol):
    cost = 0
    s = [0 for i in nodes]
    d = [0 for i in nodes]
    for i in nodes:
        for j in adj[i]:
            if sol[i] == sol[j]:
                s[i] += 1
            else:
                d[i] += 1
    for i in nodes:
        cost += d[i]
    return cost/2, s, d\

def find_move(part, nodes, adj, sol, s, d, tabu, tabulen, iteration):
    mindelta = Infinity
    istar  = None
    for i in nodes:
        if sol[i] != part:
            if tabu[i] <= iteration:
                delta = s[i] - d[i]
                if delta < mindelta:
                    mindelta = delta 
                    istar = i
    if istar != None:
        return istar, mindelta
    tabu = [0 for i in nodes]
    return find_move(part, nodes , adj, sol, s, d, tabu, tabulen, iteration) 

def move(part, nodes, adj, sol, s, d, tabu, tabulen, iteration):
    i, delta = find_move(part, nodes, adj, sol, s, d, tabu, tabulen, iteration)
    sol[i] = part
    tabu[i] = iteration + tabulen
    s[i], d[i] = d[i], s[i]
    for j in adj[i]:
        if sol[j] != part:
            s[j] -= 1
            d[j] += 1
        else:
            s[j] += 1
            d[j] -= 1
    return delta

def tabu_search(nodes, adj, sol, max_iter, tabulen):
    cost, s, d = evaluate(nodes, adj, sol)
    tabu = [0 for i in nodes]
    bestcost = Infinity
    for it in range(max_iter):
        cost += move(1, nodes, adj, sol, s, d, tabu, tabulen, it)
        cost += move(0, nodes, adj, sol, s, d, tabu, tabulen, it)
        if cost < bestcost:
            bestcost = cost
            bestsol = list(sol)
    return bestsol, bestcost


if __name__ == "__main__":
    import networkx as nx
    import matplotlib.pyplot as plt
    num_nodes = 50
    nodes, adj = rnd_graph_fast(num_nodes,0.3)
    G=nx.Graph()
    for i in range(num_nodes):
        G.add_node(i)
    for i in range(num_nodes-1):
        for j in adj[i]:
            if i < j:
                G.add_edge(i,j)
    pos=nx.spring_layout(G)

    max_iter = 1000
    tabulen = 10
    sol = construct(nodes)
    sol, cost = tabu_search(nodes, adj, sol, max_iter, tabulen)

    rnodelist = []
    bnodelist = []
    for i in range(len(sol)):
        if sol[i] == 1:
            rnodelist.append(i)
        else:
            bnodelist.append(i)
    nx.draw(G,pos,nodelist=rnodelist,node_color='r')
    nx.draw(G,pos,nodelist=bnodelist,node_color='b')
    plt.savefig("path.png")

これをちょっと変えてSimulated Annealingにできる。

12012010 Python

最大クリーク問題を禁断探索法で

エッジをネガティブな力とした時に、エッジを構成しないノードの最大集合を求める問題を最大安定集合問題といい、これは補グラフの最大クリーク問題と同値である。

ProductName メタヒューリスティクスの数理
久保 幹雄,J. P. ペドロソ
共立出版 / ￥ 3,675 ()
在庫あり。

networkx用のコード

if __name__ == "__main__":
    import networkx as nx
    import matplotlib.pyplot as plt

    num_nodes = 100
    nodes,edges = rnd_graph(num_nodes,0.5)
    cedges = complement(nodes,edges)
    adj = adjacent(nodes, cedges)
    max_iterations = 1000
    tabulength = len(nodes)/10

    G=nx.Graph()
    for i in range(num_nodes):
        G.add_node(i)
    for i in range(num_nodes-1):
        for j in adj[i]:
            if i < j:
                G.add_edge(i,j)

    pos=nx.spring_layout(G)
    sol = construct(nodes,adj)
    xcard, xinfeas, xb = evaluate(nodes, adj, sol)
    sol, card = tabu_search(nodes, adj, sol, max_iterations, tabulength)

    rnodelist = []
    bnodelist = []
    for i in range(num_nodes):
        if i in sol:
            rnodelist.append(i)
        else:
            bnodelist.append(i)

    redges = []
    for i in range(len(rnodelist)-1):
        for j in range(i+1,len(rnodelist)):
                redges.append((rnodelist[i],rnodelist[j]))


    nx.draw(G,pos,nodelist=rnodelist,edgelist=redges,node_color='r')
    nx.draw(G,pos,nodelist=bnodelist,edgelist=[],node_color='b')
    plt.savefig("max_clique.png")

max_clique

例えばHTSヒットの化合物群をMCSベースでネットワークをつくり、最大安定集合問題として解けば、部分構造ベースで非類似度の高い化合物セットを抽出できる（はず）。まぁこういうのはむしろパテントを対象とするべきなのかも。

問題はエッジをどう定義するかなんだけどなぁ(MCSで完璧とも思えないし)。いい方法ないもんかな。

12012010 Python macbook Clojure

pythonでthreadingとmultiprocessing

Python vs Clojure – Evolving « Best in Classをみてて、「ほう」となったので、手元のmacbookで。

t_seq.py 

def count(n):
    while n > 0: n -= 1

count(100000000)
count(100000000)

t_thread.py 

from threading import Thread

def count(n):
    while n > 0: n -= 1

t1 = Thread(target=count,args=(100000000,))
t1.start()
t2 = Thread(target=count,args=(100000000,))
t2.start()
t1.join()
t2.join()

t_mprocessing.py 

from  multiprocessing import Process

def count(n):
    while n > 0:
        n -= 1

t1 = Process(target=count,args=(100000000,))
t1.start()
t2 = Process(target=count,args=(100000000,))
t2.start()
t1.join()
t2.join()

結果

$ time python t_seq.py;time python t_thread.py;time python t_mprocessing.py

real    0m17.032s
user    0m16.886s
sys 0m0.067s

real    0m37.139s
user    0m29.756s
sys 0m23.039s

real    0m8.777s
user    0m16.973s
sys 0m0.083s

Clojureのほう

user=> (defn countdown [n] (when (pos? n) (recur (dec n))))

#'user/countdown
user=> user=> (time (doall (pmap countdown (repeat 2 100000000))))
"Elapsed time: 8164.085 msecs"
(nil nil)

11012010 Python

Essential AQLAlchemy

読みやすかった。これ読めばoverviewが大体つかめる。というか、engine,session,metadataの関係がつかめると理解が進む。

ProductName Essential Sqlalchemy
Rick Copeland
Oreilly & Associates Inc / ￥ 3,356 ()
通常1～3週間以内に発送

  • Sessionクラスはsessionmakerで作る。
  • metadata,engineのbound,unbound
  • generative query interface
  • 1:N,N:Mrelationsの指定の仕方
  • sessionの使い方
  • Elixir,SQLSoup

11012010 life

TopHatenar

調べ物してたら、tophatenarの自分のランキングにぶつかった。

blogopolisの区画も割り当てられてた。

blogopolis

11012010 Haskell

「プログラミングHaskell」の読書会

今週末に「プログラミングHaskell」の読書会があります。

  • 日時： 2010年1月16日（土） 13:00～17:30
  • 場所： 静岡市産学交流センター 演習室 1 (http://www.hanjyou.jp/)
  • 地図： http://www.hanjyou.jp/map.html

ProductName プログラミングHaskell
Graham Hutton
オーム社 / ￥ 2,940 ()
在庫あり。

11012010 life macbook

LinkStation 1.0TB LS-CH1.0TL

macbookが二台になったので、Time Machineの環境を揃えることにした。

ProductName BUFFALO ネットワーク対応HDD LinkStation 1.0TB LS-CH1.0TL

バッファロー / ￥ 22,050 (2009-05-02)
在庫あり。

  • timemachineの設定をonにする
  • AFPのチェックボックスをonにする
  • ユーザーを作成する
  • macアドレスはethernetのものを
    • このmacについて -> 詳しい情報