List::UtilとかList::MoreUtilsというものを教えてもらったので、少し読んでみた(というか自分用のメモ的な意味合いが強いかも)

あと、リスト操作の動作確認するために、対話型のperlshを使ってるが便利だ。

List::Util

リストの最大、最小の要素を返したりとか。reduceで同じことができるのでreduce覚えておけばよい感じ。日本語訳

perlmolのChemistry::InternalCoords ::Builderだとこんな感じで使ってる。

$ref = ${
    reduce { $a->[0] < $b->[0] ? $a : $b }
    map { [$atom->distance($_)] }
    grep { $_ ne $atom } @$atoms;
}[1];

grep,mapでリストを操作した後、reduceを使ってお望みの要素を取り出すと。
というわけで、これといった例が思い浮かばないが、シュワルツ変換っぽく使えばリストでなくて要素が返ってくると。

List::MoreUtils

anyとかallとかのリストに対する真偽判断用の関数はこっち。日本語訳

pairwise、each_array、zipはfor文を入れ子にしなくてよくなるので、コードの見通しがよくなる感じ。

あと、QSAR記述子選択用のコード書くのに便利そうなのが、minmax(min==maxの記述子はこれでさくっと落とせる)。さらに、uniq組み合わせて、さらに数値のバリエーションのない記述子も落としていけるし、any使えばnullデータの入ってる記述子も排除できるので、今度preparation用のコードをList::MoreUtils使って書いてみよう。

fp::functionals

リスト操作用だけでなく、関数操作用のモジュールもあるだろうと探したら、やっぱありました。

もうちょいhaskellまともに書けるようになったら使ってみたい（かな?）が、今のところはperlで関数型っぽくかけるのはリスト操作で十分かもという気がしてる。

WindowsXPでPoderosa3系使うとたまに落ちるのが困りもんだったが、4系は結構安定していていい感じ。

あとちょっとさくさく動くようになって快適感が増したような。

みんなのpythonを読んだら、もうちょっと実践的なコードを基に、他の言語と比較しながら説明してあるDive Into Pythonなど読むべし。pythonの初歩からhtml,xmlの処理、そしてSOAPによる通信といった割と基本的な（けど実用的な）章があって、終わりのほうは単体テスト、リファクタリングといったテスト技法や、関数型プログラミングに関しても触れてるので、全体を読んでpythonへの理解はかなり深まった(と思う)。わからないところも、みんなのpython片手にPython ライブラリリファレンスひけば、大概解決したし。

で、17章にフィボナッチ数を求める例があって、そこでジェネレータ使ってたのをみてジェネレータって便利だと理解した。

def fibonacci(max):
        a,b = 0,1
        while a < max:
                yield a
                a, b = b, a+b

for n in fibonacci(1000):
        print n

とfibonacci関数定義するだけでフィボナッチ数が。perlで書いたのと同じようなアルゴリズムなのにperlより簡単にかけているのはyield文のせい。

みんなのPython
柴田 淳
ソフトバンククリエイティブ / ?円 ( 2006-08-22 )

060914追記

perlでもクロージャ使えば同じようにかけるわ。ということに朝起きたときに気づいた

sub fibonacci {
    my $max = shift;
    my ($fib,$a,$b,) = (0,0,1);
    return sub {
        ($fib,$a, $b) = ($a,$b, $a+$b);
        if ($a < $max){
            return $fib;
        }
    }
}

my $f = fibonacci(100);

pythonにもクロージャってあるんだろうけど、ジェネレータとの使いわけってどうなるんだろうか?

ダブっていうか、ジャングルよりな気が。ラガな感じの、she's goneとかShout Ballout (Ray Keith Remix)なんかが気持ちヨサゲ。

Dub Dread 2
Ray Keith
Dread / 1594円 ( 2007-07-10 )

Blade Runnerの曲も古き良きっていうか、loving youなんかはアートコア目が強くて好き。なんといってもアートコアは浮遊感がね。あと1.5Kってのもお買い得感が高め。

抜天で送別会。2階が個室で畳の円卓になっているので、子連れもOK。しかも美味い。
子供OKの店ってあまり知らないので重宝する。

前菜4種盛りは、椎茸の含煮と砂肝が美味しかった。そしてふかひれスープも。

チキンと鮑の炒め物。鮑は柔らかく、メンバーの中でも好評だった。

なんかの揚げ物（忘れた）と海老のチリソース。

あと、炒め物がもう一品と、チャンポンと杏仁豆腐がデザートで。

紹興酒も一本空けて、満足。というわけで、また行きたい。

Statistics::RもRSPerlもイマイチ感が強かったので、rpyを使い始めているが、こんなに便利とは思わなかった。

というわけで、rpyで主成分分析を行い、できたモデルで新しいデータの主成分を求めてみるというサンプルを。 　今回はrpyのテストなのでデータはrnormで作った。普通バイオインフォとかケモインフォなパターンだと、遺伝子の発現レベルとか、化合物のディスクリプターを使うけど。

>>> from rpy import *
>>> r.library('stats')
>>> a = r.matrix(r.rnorm(15),5)
>>> b = r.matrix(r.rnorm(15),5)
>>> model = with_mode(NO_CONVERSION,r.prcomp)(a)
>>> model
<Robj object at 0xb7f9f1c0>
>>> r.predict(model,b)
array([[ 0.28817526, -2.20038126,  0.42937523],
      [ 0.76682094, -1.38050702,  0.24898564],
      [-0.28308137, -0.51609132,  0.12062971],
      [-0.68531169,  1.42599763,  0.03987597],
      [-1.35685327,  0.02394946, -1.03868517]])

と数行ほどで、予測したいデータセットの主成分がpythonのarray型で返ってくるので、このあとの処理が凄く楽チンになる。

rpyはデータのコンバージョンのやり方だけきちっと押さえておけばあとはRと一緒に扱えるのだが、コンバージョンのモードが色々あって、初めのほうは悩まされることが多かった。特にモデルをRで適用する場合、python形式にコンバートしてしまうとRで扱えなくなってしまうのでNO_CONVERSIONモードにしないといけない。

最近pythonに真面目に取り組んでいるので、Meadowにもpython-modeを入れた。

Emacs/Meadow2 - 2.00pre1

;;;Add python mode (setq auto-mode-alist (cons '("\.py$" . python-mode) auto-mode-alist)) (setq interpreter-mode-alist (cons '("python" . python-mode) interpreter-mode-alist)) (autoload 'python-mode "python-mode" "Python editing mode" t)

openbabelの分子重ね合わせツールのobfitをperlのモジュールにできんもんかな？とソース(C++)を眺めていたら、QTRFIT というアルゴリズム（実装はC++)を使っているようだ。

結局最大の重ねあわせを求めるのを固有値問題に落とし込んでるようだ。

QTRFIT

The eigenvector associated with the largest eigenvalue gives the desired quaternion. To generate a rotation matrix we consider the effects of the quaternion on the cartesian unit vectors

で、それをJacobi法で解いていると。感覚的にはそうなんだろなぁと思うが、ここら辺の理解度はちょっとあやしい。

Rをネットワークで使うには、CGIwithRで十分だろうと考えていたが、まぁいろいろあって（単にSOAPいじりたくなっただけともいう）RSOAPを入れてみることにした。

正直なところがパイソンいじったことがないので、ちょっと避けてたのもあるヨ。

用意するものは結構多い

• fpconst-0.7.2.tar.gz
• PyXML-0.8.4.tar.gz
• logging-0.4.9.5.tar.gz
• R-2.1.0.tar.gz
• session_1.0.1.tar.gz
• rpy-0.4.2.1.tar.gz
• SOAPpy-0.11.6.tar.gz
• RSOAP-1.1.3.tar.gz

まずRはrpyのために--enable-R-shlibオプションをつけてインストールしなければならない。ということで勢いtarボールから。

#./configure --enable-R-shlib #make #make install

いやー、Rのコンパイルは長い。RH4xあたりでのカーネルコンパイルくらい（古い？）こんなにかかるとは思わず、昨晩makeしたら、U隊長に大不評（僕のマシンは寝室においてあるので排気音とガリガリ音がうるさい）。で朝起きてからセッションいれとく。

#R CMD INSTALL session_1.0.1.tar.gz

tar xvfz して、パイソン関係のモジュールを順番に入れます。

#python setup.py build #su #python setup.py install

rpy だけちょっとメンドイ。まず、libR.soへのパスを通す。READMEだとbinに通してるけどlibでしょ。

#su #vi /etc/ld.so.conf ## /usr/local/lib/R/lib を記述 #ldconfig

ちなみに、setup.py も適当に書き換えないといけなかったりする場合もあるよ。-lR が見つからんぜよみたいなエラーの場合は、setup.py中の記述のbinをlibに変えるといい（まぁすぐわかる）。あとRのライブラリもtar xvfzする。僕の場合はR-2.1.0なので、

#tar xvfz R-2.1.0.tar.gz R-2.1.0/src/include #su #python setup.py install

最後にRSOAPのインストールは、普通

#python setup.py build #su #python setup.py install

ローカルでテストしたかったら、SimpleExample.pyというのを探してきて、ホスト名をlocalhostに書き換えておく。

#su #RSOAPManager

で、サーバーを立ち上げといて、別の端末から

#python SimpleExample.py

結果が戻ってくればOK

TCP/IPでやりたければ

#su #RSOAPManager　--tcp-ip

ここまではうまくいったヨ。今朝、作業したわりには記憶が曖昧なので順番が前後してるかもしれん。

みんなのpythonも読んだし、フレームワークでも試してみるかと、djangoとかTurbogearsに手をだしたら、火傷した。

• http://www-06.ibm.com/jp/developerworks/linux/060802/j_l-django.shtml
• http://www-06.ibm.com/jp/developerworks/linux/060825/j_l-turbogears.shtml

djangoは一通りサンプル動作をさせることができたけど、Turbogearsはコントローラのあたりからさっぱりわからん。
いきなりフレームワークは無謀すぎたかなということで、Biopythonでも入れて、もう少しpython慣れすることに。

easy_install numpy

wget http://www.egenix.com/files/python/egenix-mx-base-2.0.6.tar.gz
python setup.py install

svn co http://www.reportlab.co.uk/svn/public/reportlab/trunk
python setup.py install

wget http://biopython.org/DIST/biopython-1.42.tar.gz
python setup.py install

easy_installはperlでいうところのcpanコマンドみたいで楽チンです（ちゃんと動けば）。動かないのは地味にダウンロード、展開、setup.pyを行なった。

さて、インストールも無事に終了したところで、bioperlと比べてみた。比較したのはpdbファイルの処理。というのは、FMO用のインプット作るときに使っているPDBパーザーつまりBio::Structure::IOがやたらと遅く、不満だったからというありがちな理由。

コードはこんな感じで、1fatっていうpdbファイル読み込んでCalphaを出力してみた。

python

from Bio.PDB.PDBParser import PDBParser

parser=PDBParser(PERMISSIVE=1)
structure=parser.get_structure("1fat", "1fat.pdb")
for model in structure.get_list():
    for chain in model.get_list():
        for residue in chain.get_list():
            if residue.has_id("CA"):
                ca_atom=residue["CA"]
                print ca_atom.get_coord()

perl

use strict;
use Bio::Structure::IO;

my $pdb_file = "1fat.pdb";
my $structio =  Bio::Structure::IO->new(-file => "$pdb_file",
                            -format => 'PDB');

my $struc = $structio->next_structure;
for my $chain ($struc->get_chains) {
    for my $res ($struc->get_residues($chain)) {
        for my $atom ($struc->get_atoms($res)) {
            print join " ",$atom->xyz,"\n" if $atom->pdb_atomname eq " CA ";        
        }
    }
}

で、ベンチマーク

$ time python test.py
real    0m1.161s
user    0m0.996s
sys     0m0.116s

$ time perl pdbtest.pl
real    0m3.183s
user    0m3.056s
sys     0m0.044s

平均とってないけど、何回か実行してみてもpythonのほうがbioperlよりも3倍くらい速かった。bioperlのコードってなんだかなと思うことがままある。pdb_atomnameなんかも4文字で前後にスペース入ってるし。というわけで、SBDD周りのプログラミングはbiopythonのほうが扱いやすい感じがしてる。

みんなのPython
柴田 淳
ソフトバンククリエイティブ / ?円 ( 2006-08-22 )

余談だが、Turbogears普通に触れる程度までpython覚えてやる（意地でも）とか思ったので、Dive Into Pythonも読み中。