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

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

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

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

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

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

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にできる。

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

メタヒューリスティクスの数理
久保 幹雄,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")

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

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

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)

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

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

多数の化合物のファイルはSMILES形式でsqliteのデータベースにしておくと、あつかいやすいのだけど、たまにさくっとアクセスしてゴニョりたいときにsqliteでSQL文発行してとかめんどくさかったりする。

で、SQLAlchemy使うとpythonの対話環境で使えてそのままpybelとかに持っていけるので便利なんだけど、いちいちマッピングしなくても良いらしい。

create_engineでsqliteのファイルを指定して、reflectをTrueにすればよい。

from sqlalchemy import *

db = create_engine('sqlite:///test.db')
metadata = MetaData(bind=db, reflect=True)
table = metadata.tables['table_name']

このdrkcoreのファイルをつかって、タイトルを抜き出してみる。

>>> execfile("satest.py")
>>> stmt = entries.select()
>>> result = stmt.execute()
>>> for row in result:
...   print row.title
... 
blogを変えてみた
卓次郎商店でつけ麺
かど乃やで黒びしおラーメン

らくちん。

Essential Sqlalchemy
Rick Copeland
Oreilly & Associates Inc / 2719円 ( 2008-06 )

4章くらいまで読んだ。データの永続化用にsqlalchemy使えるようにしとくと仕事が非常に楽になりそうなので、ちゃんと読む。

HaskellのsplitAtはなんでタプルを返すんだろうか？

Prelude> :t splitAt
splitAt :: Int -> [a] -> ([a], [a])
Prelude> splitAt 3 [1,2,3,4,5]
([1,2,3],[4,5])

pythonで同じようなのを書いてみる。

>>> def splitAt (n,xs): return [xs[:n],xs[n:]]
... 
>>> splitAt(3,[1,2,3,4,5])
[[1, 2, 3], [4, 5]]
>>> def splitAt (n,xs): return (xs[:n],xs[n:])
... 
>>> splitAt(3,[1,2,3,4,5])
([1, 2, 3], [4, 5])

なんというか[[1, 2, 3], [4, 5]]のほうが見た目よろしい気がするのです。

追記09.12.25

• @ringtaro
• @mzp

あとで読む。

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

openbabelのringsearchがちょっと使いにくい気がする。部分構造をmolオブジェクトで返すメソッドがあっても良いんではないか？

from openbabel import *

obconversion = OBConversion()
obconversion.SetInFormat("smi")
obconversion.SetOutFormat("smi")
obmol = OBMol()

notatend = obconversion.ReadString(obmol,"c1ccccc1CC(=O)C2CCC2C(=O)O")

for r in obmol.GetSSSR():
    for a in  r._path:
        print a

これだと原子のIDしか返さない。

でも縮環チェックしたいんだったら、原子のリストがかぶるかどうかをリングの構成原子のIDをつかって調べればいいのでminimumの環からmaximumの環にするのは難しくはないか。

アトムのリストを渡したら部分構造検索を実行して、結果としてマッチした構造をmolとして返すメソッドがあればいいのか。

来年の冬に出るのか。

SQLAlchemy: Database Access Using Python (Developer's Library)
Mark Ramm,Michael Bayer
Addison-Wesley Professional / ￥ 4,874 ()
近日発売　予約可

Pythonの描画ライブラリであるmatplotlibの本が出てた。

Matplotlib for Python Developers: Build Remarkable Publication Quality Plots the Easy Way
Sandro Tosi
Packt Publishing / 4332円 ( 2009-10-30 )

欲しいのう