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

ふと、つぶやいた

例えば

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

project

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

Machine Learning: An Algorithmic Perspective 9章

Unsupervised Learning。なかなか面白かった。

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

K-meansの話から入って、K-Meansは一層のニューラルネットで表現することが出来ることを示していく。ニューラルネットで表現されたK-Meansはオンライン更新できるので、入力を一度に読み込まなくて良い(はず)。

そのあとSOMだったけど、SOMは知ってるので流した。

Machine Learning: An Algorithmic Perspective 8章

EM Algorithmとkd-Tree

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

内容はPRMLのほうが詳しい。Machine Leaningのほうはコードを読んで実装を理解するって感じだな。

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

kd-Treeは使ったことなかったけど、近傍探索はよく使うので、覚えておいて自由に使えるようにしようかな。

Machine Learning: An Algorithmic Perspective 7章

BoostingとBaggingの章。

AdaBoostのPython実装は参考になった。

でも、AdaBoostはRにあるのは知っているので、そっちを使ってみる。

library(ada)
data(iris)
iris[iris$Species!="setosa",]->iris
n<-dim(iris)[1]
trind<-sample(1:n,floor(.6*n),FALSE)
teind<-setdiff(1:n,trind)
iris[,5]<- as.factor((levels(iris[,5])[2:3])[as.numeric(iris[,5])-1])
gdis<-ada(Species~.,data=iris[trind,],iter=20,nu=1,type="discrete")
gdis=addtest(gdis,iris[teind,-5],iris[teind,5])
plot(gdis,TRUE,TRUE)

adaboost

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

Machine Learning: An Algorithmic Perspective 6章

決定木。ID3アルゴリズムの説明と実装をPythonで。

C4.5とかCARTの説明はあまりしていない、触りだけ。

CARTは統計的にテキスト解析(14) 〜テキストの分類分析1〜を読むのが良い。

C4.5とかはwekaの作者の本がよいらしい。

autoinsertでPythonスクリプトの定型文を自動的に埋め込む

普段Emacsを使うのだけど、autoinsertでスクリプト書くときの最初の二行は自動的に入るようにした。

;;auto insert
(require 'autoinsert)
(add-hook 'find-file-hooks 'auto-insert)
(setq auto-insert-directory "~/.emacs.d/insert/")
(setq auto-insert-query nil) ;自動的に挿入

(setq auto-insert-alist
      (append '(("\\.py$" . "pyheader.txt")
           ) auto-insert-alist))

pyheader.txt

1
2
3
4
#!/usr/bin/env python
# -*- encoding:utf-8 -*-

# kzfm <kerolinq@gmail.com>

これでOK

tweepyでOAuthを試してみた

PythonでOAuthを試したかったのでtweepyを使ってみた。

easy_installじゃなくてgitから

git clone git://github.com/joshthecoder/tweepy.git

まずはBasic認証

>>> import tweepy
>>> auth = tweepy.BasicAuthHandler("user", "password")
>>> api = tweepy.API(auth)
>>> api.update_status('hello from tweepy(Basic)')
<tweepy.models.Status object at 0x13c8db0>

http://twitter.com/kzfm/status/9330123865

続いてOAuth。consumer_token,consumer_secretはあらかじめtwitterから取得しておく。

>>> import tweepy
>>> consumer_token = "xxxxxxxxxxxxxx"
>>> consumer_secret = "xxxxxxxxxxxxxxxx"
>>> auth = tweepy.OAuthHandler(consumer_token, consumer_secret)
>>> redirect_url = auth.get_authorization_url()
>>> redirect_url

redirect_urlのURLにブラウザでアクセスして許可すると7桁の数字が表示されるのでそれを入力

>>> verifier = '1234567'
>>> auth.get_access_token(verifier)
<tweepy.oauth.OAuthToken object at 0x1393510>
>>> key = auth.access_token.key
>>> secret = auth.access_token.secret
>>> auth = tweepy.OAuthHandler(consumer_token, consumer_secret)
>>> auth.set_access_token(key, secret)
>>> oauthapi = tweepy.API(auth)
>>> oauthapi.update_status('hello from tweepy(OAuth)')
<tweepy.models.Status object at 0x1397b50>

http://twitter.com/kzfm/status/9330367673

wsgioauthもいいかもとか思い始めたでござる。

ProductName エキスパートPythonプログラミング
Tarek Ziade
アスキー・メディアワークス / 3780円 ( 2010-05-28 )


Scalaのzipが

なんか慣れない。

scala> Array(1,2,3) zip Array(4,5,6)
res1: Array[(Int, Int)] = Array((1,4), (2,5), (3,6))

Haskellだとこんな感じでしょ。

Prelude> [1,2,3] `zip` [4,5,6]
[(1,4),(2,5),(3,6)]

Pythonだと

[1,2,3].zip([4,5,6])

とはできない。

>>> zip([1,2,3],[4,5,6])
[(1, 4), (2, 5), (3, 6)]

PythonでK-means

Haskellな気分で。最後while入れたけど。

from itertools import groupby
from math import sqrt

zipWith = lambda f, xs, ys : [f(x, y) for x,y in zip(xs, ys)]
snd  = lambda x: x[1]
fst  = lambda x: x[0]
euclid  = lambda x, y : (x-y)**2

def distance(a,b): 
    return sqrt(sum(zipWith(euclid,a,b)))

def centroid(clusters): 
    sums = reduce(lambda x,y: zipWith(lambda a,b:a+b,x,y),clusters)
    return map(lambda x: x / float(len(clusters)), sums)

def closest(pts, pt):
    closest_ct = pts[0]
    for ct in pts[1:]:
        if distance(pt,closest_ct) > distance(pt,ct):
            closest_ct = ct
    return closest_ct

def recluster_(centroids,points):
    reclustered = [(closest(centroids,a), a) for a in points]
    reclustered.sort()
    return [map(snd,list(g)) for k, g in groupby(reclustered, fst)]

def recluster(clusters):
    centroids = map(centroid, clusters)
    concated_clusters = reduce(lambda a,b: a+b, clusters)
    return recluster_(centroids,concated_clusters)

def part(l,points):
    size = len(l)/points
    return [l[i:i+size] for i in range(0,len(l),size)]

def kmeans(k,points):
    cluster = part(k,points)
    newcluster = recluster(cluster)
    while(cluster !=  newcluster):
        cluster = newcluster
        newcluster = recluster(cluster)
    return newcluster

if __name__ == "__main__":
    pts = [[1,2,4],[1,3,3],[4,3,0],[2,5,1],[7,3,8],[0,0,0],[4,3,2],[6,1,8]]
    print kmeans(pts,3)

SQLAlchemyのDeclarative APIの使い方

SQLAlchemy 0.5から使えるようになったDeclarative APIはマッパーとテーブルとクラスを一度に定義するものらしい。

なのであとはengineをバインドすればよろしくやってくれる。

上の例だと

import datetime
from sqlalchemy import schema, types, orm

metadata = schema.MetaData()

def now():
    return datetime.datetime.now()

from sqlalchemy.ext.declarative import declarative_base

# Assign the same metadata object we created earlier.
Base = declarative_base(metadata=metadata)

# We still need the pagetag table because we don't want to explicitly define a
# Pagetag class but still
# need to specify the table in the relation between pages and tags.
pagetag_table = schema.Table('pagetag', metadata,
    schema.Column('id', types.Integer,
        schema.Sequence('pagetag_seq_id', optional=True), primary_key=True),
    schema.Column('pageid', types.Integer, schema.ForeignKey('page.id')),
    schema.Column('tagid', types.Integer, schema.ForeignKey('tag.id')),
)

class Page(Base):
    __tablename__ = 'page'

    id = schema.Column(types.Integer,
        schema.Sequence('page_seq_id', optional=True), primary_key=True)
    content = schema.Column(types.Text(), nullable=False)
    posted = schema.Column(types.DateTime(), default=now)
    title = schema.Column(types.Unicode(255), default=u'Untitled Page')
    heading = schema.Column(types.Unicode(255))
    comments = orm.relation("Comment", backref="page")
    tags = orm.relation("Tag", secondary=pagetag_table)

class Comment(Base):
    __tablename__ = 'comment'

    id = schema.Column(types.Integer,
        schema.Sequence('comment_seq_id', optional=True), primary_key=True)
    pageid = schema.Column(types.Integer,
        schema.ForeignKey('page.id'), nullable=False)
    content = schema.Column(types.Text(), default=u'')
    name = schema.Column(types.Unicode(255))
    email = schema.Column(types.Unicode(255), nullable=False)
    created = schema.Column(types.TIMESTAMP(), default=now())

class Tag(Base):
    __tablename__ = 'tag'

    id = schema.Column(types.Integer,
       schema.Sequence('tag_seq_id', optional=True), primary_key=True)
    name = schema.Column(types.Unicode(20), nullable=False, unique=True)

page_table = Page.__table__

これをmodel.pyとかしといてpython対話環境で

>>> execfile('model.py')
>>> from sqlalchemy.engine import create_engine
>>> engine = create_engine('sqlite:///test.db')
>>> metadata.bind = engine
>>> metadata.create_all()

とやってengineをバインドすればsqliteのファイルが作られる。metadataのところはテーブルからもアクセスできて

>>> page_table.metadata.bind = engine
>>> page_table.metadata.create_all()

でも同じことができる。

sqlite> .schema
CREATE TABLE comment (
    id INTEGER NOT NULL, 
    pageid INTEGER NOT NULL, 
    content TEXT, 
    name VARCHAR(255), 
    email VARCHAR(255) NOT NULL, 
    created TIMESTAMP, 
    PRIMARY KEY (id), 
     FOREIGN KEY(pageid) REFERENCES page (id)
);
CREATE TABLE page (
    id INTEGER NOT NULL, 
    content TEXT NOT NULL, 
    posted TIMESTAMP, 
    title VARCHAR(255), 
    heading VARCHAR(255), 
    PRIMARY KEY (id)
);
CREATE TABLE pagetag (
    id INTEGER NOT NULL, 
    pageid INTEGER, 
    tagid INTEGER, 
    PRIMARY KEY (id), 
     FOREIGN KEY(pageid) REFERENCES page (id), 
     FOREIGN KEY(tagid) REFERENCES tag (id)
);
CREATE TABLE tag (
    id INTEGER NOT NULL, 
    name VARCHAR(20) NOT NULL, 
    PRIMARY KEY (id), 
     UNIQUE (name)
);