2010/03/18 19:38:47
遺伝アルゴリズム
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')
2010/03/10 21:56:46
ふと、つぶやいた。
例えば
- プロジェクトで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")
んーなんかいまいち。一応活性の強弱にあわせてノードのサイズが変わるようにしてみたんだけど、大きさがわからん。
2010/01/12 11:19:02
cronとか別にディスプレイ上に表示しなくていい場合は適当に引数を与えてやることでmatplotlibのバックエンドを切り替えることができる。
some_matplotlib.py -dAgg
別にこんなことしなくてもuseしても同じことが出来るということを最近知った。
import matplotlib
matplotlib.use('Agg')
from pylab import *
ただし、pylab呼び出しよりも先にやっておかないと
RuntimeError: matplotlib.use() must be called *before* pylab
と怒られる。
2010/01/12 11:16:45
前週末から風邪がひどくて、会社を休んで寝ている。今日は多少ボーっとするけど昨日よりはずっと楽なので蒲団の中でノートパソコンと戯れたり。で、この前の血圧管理のグラフ化をやってみた。はてなのグラフを二つ重ねるというのは僕も考えたんだけど、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に固定引数が入ってんのね。
2010/01/12 11:16:30
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')
2010/01/12 11:15:50
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は入ってんのになー
2010/01/12 11:04:14
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というモジュールがあって何日おきとか何年おきとかそんな感じの目盛りが簡単にうてる。
2010/01/12 10:51:04
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のインストールが途中でこけるのだけど、再度同じコマンドを打てばいいらしい。
2010/01/12 10:50:43
元ネタは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:多変量がある決まったレンジ内に収まっているかどうかを視覚的に捉えやすい表現手段を探す。
Machine Learning: An Algorithmic Perspective (Chapman & Hall/Crc Machine Learning & Patrtern Recognition)