デコンボリューション法が紹介されていたのできちんと理解しておきたいところ。

ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ￥ 6,300 ()
在庫あり。

で、ちょっと調べたら古い論文がヒットした。

An algorithm and computer program for deconvolution in linear pharmacokinetics

fortranのコードが載っているのでpythonに移植すればいいやと軽く考えてたんだけど、go toで混乱している最中。if else 的なgo toは読みやすいんだけど、あっちいったりこっちいったりするのはちょっと辛い。

あとは式の内容ちゃんと理解してないのでコード追えないってのもあるが。 移植しながら論文の式の内容を理解していくつもりだったのだけどなかなか疲れる。

fortranの組み込み関数にsignという符号付け替え関数が用意されている。

def sign(a1,a2):
    if a2 >= 0: return abs(a1)
    else:       return -abs(a1)

古いfortranのプログラムを読んでいたら

   EPS = 1.0
10 EPS = EPS/2.0
   S = 1.0 + EPS
   IF (S .GT. 1.0) GO TO 10

というコードをみつけて？となったがすぐに限りなく小さいの求め方か！と気がついた。

pythonだとこんな感じ

eps = 1.0
while eps+1.0 > 1.0: eps = eps/2.0 
# eps=1.1102230246251565e-16

レイアウトをもう少しちゃんとしたい。

import pygraphviz as pgv
G=pgv.AGraph(strict=False,directed=True)
G.add_node('V1,C1')
G.add_node('V2,C2')
G.add_node('X1')
G.add_node('X2')
G.get_node('V2,C2').attr['shape'] = 'box'
G.get_node('V1,C1').attr['shape'] = 'box'
G.add_edge('V1,C1','V2,C2')
G.add_edge('V2,C2','V1,C1')
G.add_edge('V1,C1','X1')
G.add_edge('V1,C1','X2')
G.get_edge('V1,C1','V2,C2').attr['label'] = 'k12'
G.get_edge('V2,C2','V1,C1').attr['label'] = 'k21'
G.get_edge('V1,C1','X1').attr['label'] = 'Km,Vmax'
G.get_edge('V1,C1','X2').attr['label'] = 'K'
G.layout(prog='dot')
G.draw('pk1-3.png')

コンパートメントモデルを書くのにいちいちInkspaceを立ち上げるのもめんどうなのでGraphvizにでも頼るかと思ったらmacbookには入ってなかった。networkxもついでに入れといた。

sudo port install graphviz
sudo easy_install-2.6 pygraphviz
sudo easy_install-2.6 networkx

使う

>>> import pygraphviz as pgv
>>> d={'1': {'2': None}, '2': {'1': None, '3': None}, '3': {'2': None}}
>>> A=pgv.AGraph(d)
>>> A.layout(prog='dot')
>>> A.draw('pyg.png')

pybelでのRECAP実装が投稿されたもようだが、実際はopenbabelバインディングでの実装なのかな。

あとで試す。

オープンソースで始めるゲノム・プロテオーム・メタボローム解析
樋口 千洋
オーム社 / ￥ 6,090 ()
在庫あり。

soaplibにはCherryPyが組み込まれているので簡単にウェブサービスがつくれる。さらにrpy2を使えばRのウェブサービスがお手軽に構築できて便利。

インストールの手順はwikiを参照

以下のコードをsoapr.pyというファイル名で保存する。これは与えた数字の数だけRがrnormして返すという簡単なサービス。

from soaplib.wsgi_soap import SimpleWSGISoapApp
from soaplib.service import soapmethod
from soaplib.serializers.primitive import String, Integer, Array, Float
import rpy2.robjects as robjects

class RService(SimpleWSGISoapApp):

    @soapmethod(Integer,_returns=Array(Float))
    def rnorm(self,num):
        r = robjects.r        
        y = r.rnorm(num)
        return [i for i in y]

if __name__=='__main__':
    from cherrypy.wsgiserver import CherryPyWSGIServer
    server = CherryPyWSGIServer(('localhost',7789),RService())
    server.start()

実行するとサーバが立ち上がりhttp://localhost:7789/soapr.wsdlにアクセスするとwsdlを出力する。

テスト用のクライアント

>>> from soaplib.client import make_service_client
>>> from soapr import RService
>>> client = make_service_client('http://localhost:7789/',RService())
>>> print client.rnorm(5)
[-0.40538156526399999, 0.45889158855099998, -0.67243766066699995, 2.6881017757299999, 0.113257266309]
>>> print client.rnorm(3)
[-0.16610789340500001, 0.66591497139099998, -1.74812066455]

参考

• Documentation for rpy2
• soaplib

CPANのアップデートは追いかけているのに、PyPIのアップデートを追いかけないのは如何なものかと思い、PyPI recent updatesをLDRで購読した。

が、他に誰も購読してないが、なんか購読するサイト間違ってる？

みんなのPython 改訂版
柴田 淳
ソフトバンククリエイティブ / ?円 ( 2009-04-11 )

09.04.11 追記

間違ってた。RSS autodiscoveryだとhttp://www.python.org/をさすようになってた。

気になる

Natural Language Processing with Python
Steven Bird
Oreilly & Associates Inc / 4502円 ( 2009-06-30 )

Common Lispのdocstringについて調べてたのだけど、pythonって関数にしかdocstringを入れられないんだっけ？と思ったのでやってみた。

まず関数。__doc__でdocstring。

>>> def sumab (a,b):
...   "test sum"
...   return a + b
... 
>>> print sumab.__doc__
test sum

まぁ、ここまではok

続いて変数

>>> a = 10
>>> dir(a)
['__abs__', '__add__', '__and__', '__class__', '__cmp__', '__coerce__', 
'__delattr__', '__div__', '__divmod__', '__doc__', '__float__', '__floordiv__', 
...
'conjugate', 'denominator', 'imag', 'numerator', 'real']

おー__doc__あるじゃん。

>>> print (a.__doc__)
int(x[, base]) -> integer

Convert a string or number to an integer, if possible.  A floating point
argument will be truncated towards zero (this does not include a string
representation of a floating point number!)  When converting a string, use
the optional base.  It is an error to supply a base when converting a
non-string.  If base is zero, the proper base is guessed based on the
string content.  If the argument is outside the integer range a
long object will be returned instead.

ふむ、これは任意のstringに置き換えられるのかな？

>>> a.__doc__ = "test int"
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'int' object attribute '__doc__' is read-only

リードオンリーだった。str,float,listも一緒だった。

プリミティブな型にもdocstringはあるけど、それは書き換えられないってことかな。