matplotlib / Drkcore

Machine Learning: An Algorithmic Perspective 12章

Thu, 18 Mar 2010 19:38:47 +0919

遺伝アルゴリズム

Machine Learning: An Algorithmic Perspective (Chapman & Hall/Crc Machine Learning & Patrtern Recognition)
Stephen Marsland
Chapman & Hall / ￥ 6,529 ()
通常2～3週間以内に発送

Four Peaks Problemってのがあるらしい。目的関数が

最初のビットからの連続した1の長さか最後のビットから連続した０の長さの大きいほう
但し、両端の1,0の連続した数が10以上の場合は100ポイント獲得する。

目的関数はこんな感じ。

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

from itertools import takewhile

def o(bits):
   return len(list(takewhile(lambda x: x == 1,bits)))

def z(bits):
       return len(list(takewhile(lambda x: x == 0,reversed(bits))))

def f(bits):
   reward = 100 if o(bits) > 10 and z(bits) > 10 else 0
   return max(o(bits),z(bits)) + reward


if __name__ == "__main__":
    bits1 = [1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    bits2 = [1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0]

    for bits in [bits1,bits2]:
        print "score: %d %s" % (f(bits),bits)

#score: 20  [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
#score: 114 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

要するに素直に探索していくとローカルミニマムに落ちるようになっていて、ピークの数が4つあるのでFour Peaks Problem

連続した1,0の長さでcontour plotを描いた。

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

def sim_score():
    for x in range(100):
        for y in range(100):
            if x + y > 100:
                yield 0
            else:
                score = four_peak(x,y)
                yield score

def four_peak(x,y):
    reward = 100 if (x > 10 and y > 10) else 0
    score = max(x,y) + reward
    return score

if __name__ == "__main__":
    from pylab import *

    delta = 1
    x = arange(0, 100, delta)
    y = arange(0, 100, delta)
    X, Y = meshgrid(x, y)
    Z = array([z for z in sim_score()])
    Z.shape = 100,100

    im = imshow(Z,origin='lower' ,alpha=.9)
    colorbar(im)

    cset = contour(X,Y,Z)   
    clabel(cset,inline=1,fmt='%1.1f',fontsize=10)

    hot()
    savefig('4peaks.png')

LeadOptimizationをシミュレーションしたい

Wed, 10 Mar 2010 21:56:46 +0919

ふと、つぶやいた。

例えば

プロジェクトで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")

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

matplotlibのバックエンドを切り替える

Tue, 12 Jan 2010 11:19:02 +0919

cronとか別にディスプレイ上に表示しなくていい場合は適当に引数を与えてやることでmatplotlibのバックエンドを切り替えることができる。

some_matplotlib.py -dAgg

別にこんなことしなくてもuseしても同じことが出来るということを最近知った。

import matplotlib
matplotlib.use('Agg')
from pylab import *

ただし、pylab呼び出しよりも先にやっておかないと

RuntimeError:  matplotlib.use() must be called *before* pylab

と怒られる。

SQLObjectとmatplotlibで血圧管理

Tue, 12 Jan 2010 11:16:45 +0919

前週末から風邪がひどくて、会社を休んで寝ている。今日は多少ボーっとするけど昨日よりはずっと楽なので蒲団の中でノートパソコンと戯れたり。で、この前の血圧管理のグラフ化をやってみた。はてなのグラフを二つ重ねるというのは僕も考えたんだけど、CUIから

mybp -b 142,88

とか

mybp -b 142,88 -d 2007-07-31

とタイプするだけのほうがなんとなく長続きしそうだし（今までの経験から）。とここまでエントリを書いたところで、はてなグラフでそんなのやってるの思い出した。

むむーですよ。ま、pythonで書きたかった+綺麗なグラフが好きだった、、、ということで、書いたのはココ。

O/RマッパはSQLObjectで。クラスを定義したら、

class BloodPressure(SQLObject):
    sbp = IntCol()
    dbp = IntCol()
    date = DateCol(default=datetime.date.today(), unique = True)

createTableメソッド呼ぶ

BloodPressure.createTable()

これで、実際にテーブルつくってくれるのでSQL文考えなくてよいし、テーブル作り直したいときにもdb消してもう一度createTableメソッド呼べばいいので心理的にも楽。

matplotlibは昨日書いたのを。データが少ないと目盛りがおかしいけど、1週間も続ければきちんとなるでしょう。

それにしても僕の血圧やたら高いなぁ。困った困った。

あとoptparseで固定引数の処理の仕方がよくわからなかったのだけど、パスワード生成スクリプトをきっかけにoptparseを使ってみましたを眺めたら解決した。

(options, args) = parser.parse_args()

ここのargsに固定引数が入ってんのね。

matplotlibを使う

Tue, 12 Jan 2010 11:16:30 +0919

pythonのグラフ描画ライブラリであるmatplotlibで日付ごとにプロットするサンプル。

from pylab import *
from matplotlib.dates import MONDAY, WeekdayLocator
import datetime

start_date = datetime.date(2007, 5, 1)
end_date   = datetime.date.today()
delta = datetime.timedelta(days=1)

mondays   = WeekdayLocator(MONDAY)

dates = drange(start_date,end_date,delta)

s = rand(len(dates)) + 85
t = rand(len(dates)) + 80

fig = figure()
ax = fig.add_subplot(111)
ax.plot_date(dates, s, 'r--')
ax.plot_date(dates, t, 'bo-')

ax.xaxis.set_major_locator(mondays)
ax.autoscale_view()
ax.grid(True)
fig.autofmt_xdate()

title('My Blood Pressure')
xlabel('Date')
ylabel('Blood Pressure')

fig.savefig('datetest')

matplotlibのimportエラー

Tue, 12 Jan 2010 11:15:50 +0919

matplotlibを0.91.1にあげたら

NameError: name 'gtk' is not defined

とかでて動かなくなった。仕方ないのでsvnのrev4730を入れたら動いたけど、それでも

  libpng: 1.2.22
 Tkinter: Tkinter: 50704, Tk: 8.4, Tcl: 8.4
wxPython: 2.8.4.0
         * WxAgg extension not required for wxPython >= 2.8
    Gtk+: no
         * pygtk present but import failed
      Qt: Qt: 3.3.8, PyQt: 3.17.3
     Qt4: no
   Cairo: 1.4.0

となっている。pygtkは入ってんのになー

pythonのdatetimeモジュール

Tue, 12 Jan 2010 11:04:14 +0919

6.10 datetime -- 基本的な日付型および時間型みながら。

>>> import datetime
>>> datetime.date.today()
datetime.date(2007, 7, 28)
>>> d = datetime.date.today()
>>> d.month
7
>>> d.day
28
>>> d.year
2007

で、matplotlibにはグラフを描きやすいようにmatplotlib.datesというモジュールがあって何日おきとか何年おきとかそんな感じの目盛りが簡単にうてる。

macbookにscipy,ipython,matplotlibなどをインストール

Tue, 12 Jan 2010 10:51:04 +0919

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のインストールが途中でこけるのだけど、再度同じコマンドを打てばいいらしい。

matplotlibでレーダーチャート

Tue, 12 Jan 2010 10:50:43 +0919

元ネタは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:多変量がある決まったレンジ内に収まっているかどうかを視覚的に捉えやすい表現手段を探す。

matplotlibで格子のパターンをつくる

Tue, 12 Jan 2010 10:37:24 +0919

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

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