ふと、つぶやいた。

例えば

• プロジェクトで500検体合成
• ひとつの親あたり5-25化合物くらいの子を合成する
• プロジェクトの人員の制限により、一度に10ラインまでしか同時に走らない
  • 11番目以降はそれ以降のオプティマイゼーションは行われない

その時知りたいことは、

• {計画、合成、アッセイ、解析}というサイクルは何回まわるのか、
• 有望そうな化合物だけど、人的リソースの関係で選択されなかった歴史はどこにあるか？
• バックアップ候補の化合物はどこら辺で現れるか（初期？中期？）

あたり。

現実の系に近づけるために、親化合物にはポテンシャルがあって、近傍探索すると、ポテンシャルの近くで活性が変化するモデル。あと同じライン（ブランチ）を継続して合成すると、ポテンシャル曲線が底に近づくのと、合成もSARに関する知識がついてくるので、活性は上がりやすくなって、かつ変動幅も小さくなるようにweightはだんだん小さくなるようにしてみた。

#!/usr/bin/env python
# -*- encoding:utf-8 -*-

# kzfm <kerolinq@gmail.com>

from random import random
import networkx as nx
import matplotlib.pyplot as plt


class Branch(object):
    """LeadOptimization flow"""

    count = 0

    @classmethod
    def inc_count(cls):
        cls.count = cls.count + 1
        return cls.count

    @classmethod
    def get_count(cls): return cls.count

    def __cmp__(self, other):
        return cmp(self.act, other.act)

    def __init__(self,potency,weight):
        self.num = Branch.inc_count()
        self.potency = potency
        self.weight  = weight
        self.act = self.potency + self.weight * random()

    def make_child(self,num_childs,potency,weight):
        return [Branch(self.potency + self.weight * random(),self.weight * 0.7) for i in range(num_childs)]

if __name__ == "__main__":
    max_comps        = 500
    initial_branches = 3
    nodesize         = []
    nrange           = []
    erange           = []
    generation       = 1
    heads = [Branch(0.3,1) for i in range(initial_branches)]
    G=nx.Graph()

    for h in heads: 
        nodesize.append(h.act * 30)
        nrange.append(generation)

    while Branch.get_count() < max_comps:
        new_heads = []
        generation += 1
        for branch in heads[:10]:
            for new_comp in branch.make_child(int(5+20*random()),branch.potency,branch.weight):
                nodesize.append(new_comp.act * 30)
                nrange.append(generation)
                erange.append(generation)
                G.add_edge(branch.num,new_comp.num)
                if new_comp.act > 0.75:
                    new_heads.append(new_comp)
        heads = new_heads
        heads.sort()

    nx.draw_spring(G, node_size=nodesize, node_color=nrange, edge_color=erange,alpha=0.7, \
cmap=plt.cm.hot, edge_cmap=plt.cm.hot, with_labels=False)
    plt.savefig("proj.png")

んーなんかいまいち。一応活性の強弱にあわせてノードのサイズが変わるようにしてみたんだけど、大きさがわからん。

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

Matplotlib for Python Developers
Sandro Tosi
Packt Publishing / ￥ 4,216 ()
通常2～4週間以内に発送

欲しいのう

PRML11章。一様乱数から正規乱数を発生させる。

from math import sin,cos,pi,log,sqrt
from random import random

def box_muller():
    z1 = random()
    z2 = random()
    return sqrt(-2*log(z2)) * sin(2*pi*z1)

if __name__ == "__main__":
    from pylab import *
    x = [box_muller() for i in range(10000)]
    hist(x,20)
    savefig("box_muller.png")

パターン認識と機械学習 下 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ￥ 8,190 ()
在庫あり。

macbookでmatplotlibを使っていたら出たエラー。

LookupError: unknown encoding: X-MAC-JAPANESE

/opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/matplotlib/cbook.pyというファイルを書き換えた

def unicode_safe(s):
    if preferredencoding is None: return unicode(s)
    elif preferredencoding == 'X-MAC-JAPANESE': return unicode(s)
    else: return unicode(s, preferredencoding)

ローレンツのついでにgumowski

class Gumowski:
    "Gumowski-Mira Map"
    def __init__(self, x0, y0, a):
        self.x0 = x0; self.y0 = y0
        self.a  = a;

    def f(self,x):
        return self.a*x + 2*(1-self.a)*x**2/(1+x**2)

    def calc(self, n):
        dat = []
        x,y = self.x0, self.y0
        for i in range(n):
            nx = y  + self.f(x)
            ny = -x + self.f(nx)
            dat.append([nx,ny])
            x,y = nx,ny
        return dat

if __name__ == "__main__":
    from pylab import *
    from random import random
    n = 3
    for i in range(n*n):
        p = str(n) + str(n) + str(i+1)
        subplot(p)
        a = random() * 0.1 + 0.3
        gumow = Gumowski(20, 20, a)
        dat = gumow.calc(1000)
        x = [d[0] for d in dat]
        y = [d[1] for d in dat]
        plot(x,y,'bo', alpha=0.1,ms=5)

    savefig("gumow.png")

どう書くより。

python+matplotlibで。

class Lorenz:
    "Lorenz attractor"
    def __init__(self, x0, y0, z0, p, r, b, dt):
        self.x0 = x0; self.y0 = y0; self.z0 = z0
        self.p  = p;  self.r  = r;  self.b  = b
        self.dt = dt

    def calc(self, n):
        dat = []
        x,y,z = self.x0, self.y0, self.z0
        for i in range(n):
            dx = (-1 * self.p * x + self.p * y) * self.dt
            dy = (-x * z + self.r * x - y) * self.dt
            dz = (x * y - self.b * z) * self.dt
            x += dx; y += dy; z += dz
            dat.append([x,y,z])
        return dat

if __name__ == "__main__":
    from pylab import *
    lorenz = Lorenz(1.0, 1.0, 1.0, 10.0, 28.0, 8.0/3.0, 0.01)
    dat = lorenz.calc(5000)
    x = [d[0] for d in dat]
    y = [d[1] for d in dat]
    plot(x,y)
    savefig("lorenz.png")

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

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

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

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

ちょっと描きたいものがあったので調べたらpcolormeshを使えばいいだけだった。

from numpy import *
from pylab import *
a = rand(100,100)
pcolormesh(a)
colorbar()
savefig('colour.png')

元ネタはRadar / Spider Chars

五角形にしたかったのでrule of fiveにPolar Surface Areaを加えておいた。

#!/usr/bin/env python
from matplotlib.projections.polar import PolarAxes
from matplotlib.projections import register_projection
from pylab import *

class RadarAxes(PolarAxes):
    """Class for creating a radar chart (a.k.a. a spider or star chart)        
    http://en.wikipedia.org/wiki/Radar_chart
    """
    name = 'radar'
    # use 1 line segment to connect specified points
    RESOLUTION = 1

    def draw_frame(self, x0, y0, r):
        verts = [(r*cos(t) + x0, r*sin(t) + y0) for t in theta]
        return Polygon(verts, closed=True)

    def set_varlabels(self, labels):
        self.set_thetagrids(theta * 180/pi, labels)

    def get_axes_patch(self):
        x0, y0 = (0.5, 0.5)
        r = 0.5
        return self.draw_frame(x0, y0, r)

if __name__ == '__main__':

    register_projection(RadarAxes)
    N = 5

    theta = 2*pi * linspace(0, 1, N+1)[:-1]
    theta += pi/2
    labels = ['HBA', 'HBD', 'cLogP', 'MWT', 'PSA']
    rule_of_five = [10, 5, 5, 500, 140]
    desc = [12, 3, 3.6, 532, 160]
    desc_rate = [100*desc[i]/float(v) for (i,v) in enumerate(rule_of_five)]

    ax = subplot(111, projection='radar')

    ax.fill(theta, [100]*N)
    ax.fill(theta, desc_rate)

    for patch in ax.patches:
        patch.set_alpha(0.5)

    ax.set_varlabels(labels)
    rgrids((20, 40, 60, 80, 100))

    grid(True)
    show()

rule of fiveのようにある範囲内に収まっていること(超えるとリスク)というような指標を表すのにレーダーチャートは適しているんだろうか。つまり充足している事を示すような面積の表現はいいのかなぁ。あと、レンジが負になったりするのでそれもどうかと思う。

再考の余地はあるな。

TODO:多変量がある決まったレンジ内に収まっているかどうかを視覚的に捉えやすい表現手段を探す。

portで。gccのコンパイルに結構時間がかかった(3時間くらい?)。途中で放って寝たので知らんけど。

sudo port install py-f2py
sudo port install py25-numpy
sudo port install py25-scipy
sudo port install py25-matplotlib +tkinter
sudo port install py25-ipython

あと、scipy,matplotlibのインストールが途中でこけるのだけど、再度同じコマンドを打てばいいらしい。