どんな本かな？と思い、目次を眺めてみたら、なかなか面白そうだった。

Python Algorithms: Mastering Basic Algorithms in the Python Language
Magnus Lie Hetland
Apress / ￥ 4,294 ()
近日発売　予約可

1-3章読んだ。

どちらかというと、リファレンスっぽい。

あと1章は正誤表がないと動かすことができないので、これを見ながら。

Kayはこれからというか、ソースコードをちょっと眺めはじめたくらい。

Google App Engineプログラミング入門
中居 良介,岡野 真也,船井 裕,松尾 貴史
アスキー・メディアワークス / ￥ 1,890 ()
在庫あり。

結局、自分でなんかつくってみないと覚えないわなという、当たり前の結論に達したので、なんか作ってみようと思ったのであった。

flask-gae-skeletonとかも試さないといけない。

最近寒いですね。昨晩は御殿場のうえのほうで雪が降って、ノーマルタイヤなのに道凍ってたらどうしようかなぁと。ちょっと悩んだけど、結局高速道路で帰った。((((;゜Д゜)))ガクガクブルブル

さて、これはチャーチ数で2^2=4、2^4=16を計算し、「w」を16個出力します。(grass.elのぺえじから)

((;゜Д゜))((;゜Д゜))(((;゜Д゜)))(((;゜Д゜)))((;゜Д゜))(((;゜Д゜)))(((;゜Д゜)))(((;゜Д゜)))((;゜Д゜))(;゜Д゜)((;゜Д゜))(((;゜Д゜)))(((;゜Д゜)))((;゜Д゜))((;゜Д゜))(((;゜Д゜)))((;゜Д゜))((;゜Д゜))((;゜Д゜))(((;゜Д゜)))((;゜Д゜))((;゜Д゜))((;゜Д゜))((;゜Д゜))((;゜Д゜))(((;゜Д゜)))((;゜Д゜))((;゜Д゜))((;゜Д゜))((;゜Д゜))((;゜Д゜))((;゜Д゜))((;゜Д゜))((;゜Д゜))

$ ./gkbr.py church.gkbr 
wwwwwwwwwwwwwwww

最近計算論を読んでいて、型なしラムダ計算の実装としてgrassを追っかけていたら思いついたので書いてみた。

Python実装のデバッグ出力には大変世話になった。先人のコードが読めるというのは素晴らしい。

あと数式でラムダ計算を眺めると外から見た感じになって理解を深めるのに良いですな。

計算論 計算可能性とラムダ計算 (コンピュータサイエンス大学講座)
高橋 正子
近代科学社 / ￥ 3,570 ()
在庫あり。

でもこの本は難しい。2章までは読んだけど、3章はわからん。

#!/usr/bin/env python
# -*- encoding:utf-8 -*-
"""
 gkbr.py - Gkbr interpreter

 2010-11-19
   - Modified by kzfm

---- Original Grass --
 grass.py - Grass interpreter
 http://www.blue.sky.or.jp/grass/

 Copyright (C) 2008 NISHIO Hirokazu. All rights reserved.

 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions
 are met:
 1. Redistributions of source code must retain the above copyright
    notice, this list of conditions and the following disclaimer.
 2. Redistributions in binary form must reproduce the above copyright
    notice, this list of conditions and the following disclaimer in
    the documentation and/or other materials provided with the
    distribution.

 THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS''
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS
 BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
 BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN
 IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 History:
 2008-06-06
   - Tlanslated to python by NISHIO
 2007-10-02
   - Follow the latest changes of the definition of Grass.
   - by UENO Katsuhiro
 2007-09-20
   - First version by UENO Katsuhiro.

"""
from copy import deepcopy
import sys
import re
import codecs

RELEASE = 0
ONLY_RETURN = 40
DEBUG = 50
DEBUG2 = 60
LOGLEVEL = RELEASE
NUMERICAL_OUTPUT = False
def log(level, *msg):
    if level <= LOGLEVEL:
        print "\t".join(map(str, msg))

def Struct(*keys):
    class _Struct(object):
        def __init__(self, *values):
            self.__dict__.update(zip(keys, values))
    return _Struct

Machine = Struct("code", "env", "dump")


class Value(object):
    def __repr__(self):
        return self.__class__.__name__

class Insn(object):
    pass

class App(Insn):
    def __init__(self, m, n):
        self.m = m
        self.n = n

    def eval(self, m):
        f, v = m.env[-self.m], m.env[-self.n]
        f.app(m, v)

    def __repr__(self):
        return "App(%(m)d, %(n)d)" % self.__dict__


class Abs(Insn):
    def __init__(self, body):
        self.body = body

    def eval(self, m):
        m.env.append(Fn(self.body, deepcopy(m.env)))

    def __repr__(self):
        return "Abs(%s)" % self.body

class Fn(Value):
    count = 0
    name = ""
    def __init__(self, code, env):
        self.code, self.env = code, env
        Fn.count += 1
        self.count = Fn.count

    def app(self, m, arg):
        m.dump.append((m.code, m.env))
        m.code, m.env = deepcopy(self.code), deepcopy(self.env)
        m.env.append(arg)

    def __repr__(self):
        if self.name:
            return self.name
        return "Fn%d" % self.count

ChurchTrue  = Fn([Abs([App(3,2)])], [Fn([],[])])
ChurchTrue.name = "CTrue"
ChurchFalse = Fn([Abs([])], [])
ChurchFalse.name = "CFalse"

class CharFn(Value):
    def __init__(self, char_code):
        self.char_code = char_code

    def app(self, m, arg):
        if self.char_code == arg.char_code:
            ret = ChurchTrue
        else:
            ret = ChurchFalse
        m.env.append(ret)

    def __repr__(self):
        return "Char(%s)" % self.char_code

class Succ(Value):
    def app(self, m, arg):
        m.env.append(CharFn((arg.char_code + 1) & 255))

class Out(Value):
    def app(self, m, arg):
        if NUMERICAL_OUTPUT:
            sys.stdout.write("%d(%c)" % (arg.char_code, arg.char_code))
        else:
            sys.stdout.write(chr(arg.char_code))
        m.env.append(arg)


class In(Value):
    def app(self, m, arg):
        ch = sys.stdin.read(1)
        if ch == "":
            ret = arg
        else:
            ret = CharFn(ord(ch))
        m.env.append(ret)

def eval(m):
    while True:
        if not m.code:
            if not m.dump: break
            ret = m.env[-1]
            m.code, m.env = m.dump.pop()
            m.env.append(ret)
            log(ONLY_RETURN, m.env)
        else:
            insn = m.code.pop(0)
            insn.eval(m)

    return m.env[0]

InitialEnv = [In(), CharFn(ord("w")), Succ(), Out()]
InitialDump = [[[], []], [[App(1, 1)], []]]

def start(code):
    return eval(Machine(code, deepcopy(InitialEnv), deepcopy(InitialDump)))

def parse(src):
    code = []
    src = src.replace(u'(((;゜Д゜)))','W').replace(u'((;゜Д゜))','w').replace(u'(;゜Д゜)','v')
    src = re.subn("[^w]*", "", src, 1)[0]
    src = re.sub("[^wWv]", "", src)
    for s in re.split("v+", src):
        if not s: continue
        a = re.findall(r"w+|W+", s)
        a = map(len, a)
        arity = 0
        if s[0] in "w":
            arity = a.pop(0)
        if len(a) % 2 != 0: raise RuntimeError("parse error at app")
        body = []

        for i in range(0, len(a) - 1, 2):
            body.append(App(a[i], a[i+1]))

        for i in range(arity):
            body = [Abs(body)]

        code += body

    return code

def run(src):
    start(parse(src))

def run_stdin():
    sys.stdin  = codecs.getreader('utf_8')(sys.stdin)
    src = sys.stdin.read()
    run(src) 

def run_from_file():
    with codecs.open(sys.argv[1], 'r','utf_8') as f:
        run(f.read())

if __name__ == "__main__":
    if len(sys.argv) > 1:
        run_from_file()
    else:
        run_stdin()

こう書いてもいいのか。

def Struct(*keys):
    class _Struct(object):
        def __init__(self, *values):
            self.__dict__.update(zip(keys, values))
    return _Struct

なるほど。

上巻読み終えたけど、この本はグラフ理論のさわりをおさえるのにちょうどいい難易度だ（難しすぎず、かといって平易すぎもしない）

離散数学への招待〈上〉
J. マトウシェク,J. ネシェトリル
シュプリンガー・フェアラーク東京 / ￥ 2,835 ()
在庫あり。

• 第2章 組合せ的数え上げ
  • 2.1 関数と部分集合
  • 2.2 置換と階乗
  • 2.3 二項係数
  • 2.4 評価 ―― 入門編
  • 2.5 評価 ―― 階乗関数
  • 2.6 評価 ―― 二項係数
  • 2.7 包除原理
  • 2.8 クローク係嬢の問題
• 第3章 グラフ理論入門
  • 3.1 グラフの概念 ―― 同型
  • 3.2 部分グラフ，連結成分，隣接行列
  • 3.3 次数列
  • 3.4 オイラー・グラフ
  • 3.5 オイラー回路を求めるアルゴリズム
  • 3.6 オイラー有向グラフ
  • 3.7 2-連結性
• 第4章 木
  • 4.1 木の定義と特徴づけ
  • 4.2 木の同型
  • 4.3 グラフの全域木
  • 4.4 最小全域木問題
  • 4.5 ヤルニークとボルーフカのアルゴリズム
• 第5章 グラフを平面に描く
  • 5.1 平面や曲面の上の描画
  • 5.2 平面的グラフの中の閉路
  • 5.3 オイラーの公式
  • 5.4 地図の色分け ―― 四色定理

1.4

X Yを任意の集合とする。直感的には、関数fとはXの各元にYの元を一つだけ割り当てる何かである

という記述とその次のページからの関数合成とp.36の関数合成の集合のイメージでHaskellの関数合成のイメージがかなり強く沸いた。

3.3

次数列定理って面白いなと調べたら、Havel-Hakimiアルゴリズムというらしい。perlで書いてあったので、pythonでかいてみた。

def hakimi(deg):
    deg.sort()

    n = len(deg)
    m = deg[-1]

    if m == 0:           return True
    if m < 0:            return False
    if m > 0 and n == 1: return False

    new_deg = []
    for i,v in enumerate(deg[:-1]):
        if i+1 < n-v:
            new_deg.append(v)
        else:
            new_deg.append(v-1)
    return hakimi(new_deg)


if __name__ == "__main__":
    print [1,1,1,2,2,3,4,5,5]
    print hakimi([1,1,1,2,2,3,4,5,5]) # True
    print [2,2,2]
    print hakimi([2,2,2]) # True
    print [3,3,3]
    print hakimi([3,3,3]) # False

4

木の章もおもしろかった(符号化アルゴリズムとか)。

クラスカルのアルゴリズムやヤルニークとボルーフカのアルゴリズム（プリムのアルゴリズム）は読みなおした。そっちに興味がわくようだったら、最短経路の本がより深いのでよいかも。

5

4色定理はいわんとしていることは分かるが証明の過程をきちんと理解しているとは言いがたい。確か図書館に四色定理の本が置いてあったので今度借りてみるかな

四色問題
ロビン・ウィルソン
新潮社 / ￥ 2,310 ()
在庫あり。

これは評判通りの良書。

Scheme手習い
Daniel P. Friedman,Matthias Felleisen
オーム社 / ￥ 2,940 ()
在庫あり。

最初のほうはさらさら読めて、僕は7章からが面白かった。

7章では関数とは何か？ということが集合から説明されていく。有限関数の表現の仕方とか一対一対応とか。

8章はlamdaから始まって継続の話でこれは面白いがさらっと流れる。継続とは何かということは良く考えないと出てこないかも

9章はコラッツの問題,アッカーマン関数が出てきて最後にY-combinator

pythonで書いてみた( fact(5) )。

(lambda g: (lambda f: g(lambda arg: f(f)(arg))) 
(lambda f: g(lambda arg: f(f)(arg))))
(lambda f: lambda n: (1 if n<2 else n*f(n-1)))(5)

10章はSchemeを実装する話で、Javaなんかで実装したことがあるが、こっちのほうが環境の動きというものがよくわかった。

8,9,10は色々分かっている人向けかな。でもちょっとは知っていればかなり理解が進むし、深まり度も高まってオススメかも。

Izu.R #0お疲れ様でした。いつもは22時前には寝るのだけど、2時半くらいまで起きてコード書いてた。たのしす。

今回が初めてで、Rらしく発表+LT的なほうがいいのか、ハッカソン的なもののほうがいいのか、わからなかったうえに初めての宿という不安はあったが、とりあえず集まってみて様子見しようという感じでしたが、蓋を開けてみれば、楽しい会になってくれて良かった。

次回もやる方向で考えているので、興味があれば。

今回のキーワード

Galaxy

pythonで書かれたワークフロー管理ツール。@bonohu さんに初期の導入のあたりを教えてもらいながら。みんな、ほとんどこれいじってたような。

• Galaxy ツールを作る

ちなみにmacportsのpythonだと依存関係トラブリまくってうまく動かないので、潔く/usr/binのを使えという結論に落ち着いた。

これはもう少しいじる。次回にはなんか成果物をもっていきたいかな。

別トラックで流れていたセッションの情報交換

これはよい。あとで内容をきちんと把握する。あの、ユーザー会は一定の確率で面白い内容出てくるからよいです。面白さを面白いと感じるためには、面白さに対して感度をあげる努力を怠るとダメですな。

ジャパネットムスカ

きっかけは覚えてないが、MADの方向に行ってきた。

ついでに、一人ふぁぼ界に引きずり込まれてた。

Python

• Dive into Python

みんなのPythonは良書ですよね。

みんなのPython 改訂版
柴田 淳
ソフトバンククリエイティブ / ￥ 2,940 ()
在庫あり。

そういえば、前日の二次会でもPython利用者率を上げていこうという結論で落ち着いた(Pythonistaの中で)

書いたもの

openbabelで荒く構造立ち上げて、その後量子化学計算で精密化して、Mulliken Chargeとか楽に求めたいこと多いので、サクサク出来るようにラッパー欲しかったので書いてた。 ファイルの読み込みとかよくわからんバグが出ててはまったが、なんとか動くとこまでいった。

あとはもう少しちゃんと書いてクックブックにでもあげとく。

import openbabel as ob
from tempfile import mkstemp, mkdtemp
from os import removedirs, unlink, system, environ
import re
import os.path
import string
from random import choice

def randstr(n):
    "ランダムなファイル名を生成するため"
    return u''.join(choice('abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ') for i in xrange(n))

class Gamess(object):

    def __init__(self):
        self.tempdir = mkdtemp()
        self.gamess  = "/Users/kzfm/gamess/rungms"
        self.jobname = ''
        self.cwd    = os.getcwd()

    def calc(self, mol):
        self.jobname = randstr(6)
        obc = ob.OBConversion()
        obc.SetInFormat("gamout")

        gamin = self.write_file(mol)
        gamout = self.tempdir + "/" + self.jobname + ".out"
        os.chdir(self.tempdir)
        os.system("%s %s> %s  2> /dev/null" % (self.gamess, self.jobname, gamout))        

        # エラーが出たのでstringを渡したらなおった
        # TypeError: in method 'OBConversion_ReadFile', argument 3 of type 'std::string'
        new_mol = ob.OBMol()
        s = open(gamout).read()
        next = obc.ReadString(new_mol, s)

        os.chdir(self.cwd)
        unlink(gamin)
        unlink(gamout)
        return new_mol

    def header(self):
        h = """ $BASIS GBASIS=STO NGAUSS=3 $END
 $CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE $END
 $STATPT OPTTOL=0.0001 NSTEP=20 $END"""

        return h

    def write_file(self,mol):
        obc = ob.OBConversion()
        obc.SetOutFormat("gamin")

        gamess_input_file = self.tempdir + "/" + self.jobname +".inp"
        gamess_input_str = obc.WriteString(mol)

        h = self.header()
        ng = gamess_input_str.replace(" $CONTRL COORD=CART UNITS=ANGS $END",h)

        with open(gamess_input_file, "w") as f:
            f.write(ng)
        return gamess_input_file

    def __del__(self):
        #print "remove " + self.tempdir
        removedirs(self.tempdir)

if __name__ == '__main__':

    g = Gamess()

    obc = ob.OBConversion()
    obc.SetInFormat("mol")

    mol = ob.OBMol()
    next = obc.ReadFile(mol,"test.mol")

    newmol = g.calc(mol)
    print [obatom.GetPartialCharge() for obatom in ob.OBMolAtomIter(newmol)]

GAMESSに関しては動く環境をもっているのが前提なので、そこをやりたい場合はこの本が良いです。

すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)

講談社 / ￥ 3,360 ()
在庫あり。

今日の朝、時計をもっていないことに気づいて。

import sys,os
import Growl
import threading

def timenotify(notify_string):
    image = Growl.Image.imageFromPath("/Users/kzfm/.tw_growl/kzfm.png")
    g = Growl.GrowlNotifier(applicationName='Timer', notifications=['Watch'])
    g.register()
    g.notify(
        noteType='Watch', 
        title="GrowlTimer", 
        description= notify_string,
        icon=image,
        sticky=False)

if __name__ == "__main__":
    t1 = threading.Timer(600, timenotify, args=["10分経過しました"])
    t2 = threading.Timer(1200, timenotify, args=["20分経過しました"])
    t3 = threading.Timer(1500, timenotify, args=["25分経過しました。そろそろまとめ"])
    t4 = threading.Timer(1800, timenotify, args=["30分経過しました。発表終了"])
    t1.start()
    t2.start()
    t3.start()
    t4.start()

やっつけスクリプト

SAR News No19に寄稿しました。この号はSVM,RFとかの統計手法を使ったSARから可視化とかグラフっぽい処理とか、ToolKitを使ったプログラミングとか、日本ではあまり見かけない内容なので面白いのではないでしょうか？

余談ですがハッカソンっぽい場もあるので、ちと書いてみるか（または飲んでみるか）という方がいれば連絡を下さい(まだ空きあります)。僕はOpenBabelを使ってGAMESSの量子化学計算で構造最適化をするラッパーを書く予定にしています。それかHaskellのchemoinformaticsライブラリのソース読んでます。

さて、今回の原稿を書くのにSphinxを使った。これ最高ですな、EPubも出力できるし、超便利。ただ、今回は原稿の提出がwordだったのでwordの行間調整したりとかよくわからない部分で難儀したけど(word嫌い)。でも履歴管理の仕組み（変更モードって言うの？）はメールでやり取りするときに便利だったけど。

あとバージョン管理はMercurialを使った。Gitに比べてwebアクセスが簡単に設定できたのでチョイス、あと文書管理だとそんなに複雑な操作しないし。

入門Mercurial Linux/Windows対応
藤原 克則
秀和システム / ￥ 2,310 ()
在庫あり。

以前、おもむろに思いついたLeadOptimizationのシミュレーターで、networkxはgml出力できてCytoscapeでimportすればいいじゃんということに気づいてヤル気が出た。

ついでにコードもちょっと直して、活性はpIC50とかそういうレンジにしてみて、μMオーダーから上げていくようにした。

エッジは世代を表していて、黄色から緑にむけてすすんでいく。ノードの色は活性を表していて黒から赤になるにしたがって活性が向上する。最初の方は活性が振れるけど段々向上していって安定するようにしてみた。

リードホッピングとは選択されなかったより良い未来を選択しなおすことであるとするならばそのような希望はLO初期のあたりに存在すると思うのだ。

以下、書いたコード

from random import random
import networkx as nx

class Branch(object):
    """LeadOptimization flow
    potency: pIC50 or such
    weight : range for activity
    """

    count = -1

    @classmethod
    def inc_count(cls):
        cls.count = cls.count + 1
        return cls.count

    @classmethod
    def get_count(cls): return cls.count

    def __init__(self,potency,weight):
        self.id = Branch.inc_count()
        self.potency = potency
        self.weight  = weight
        self.activity = self.potency + self.weight * random()

    def make_child(self,num_childs,potency,weight):
        return [Branch(self.potency + self.weight * (random()-0.5)*2,self.weight * 0.5) for i in range(num_childs)]

if __name__ == "__main__":
    max_comps        = 500 # total compounds
    initial_branches = 1   # number of lead compounds
    lead_potency     = 5   # lead activity
    generation       = 0

    G=nx.Graph()

    heads = [Branch(lead_potency,2) for i in range(initial_branches)]
    map(lambda b: G.add_node(b.id, activity=b.activity),heads)

    while Branch.get_count() < max_comps:
        new_heads = []
        generation += 1
        for branch in heads[:int(2+5*random())]:
            for new_comp in branch.make_child(int(40*random()),branch.potency,branch.weight):
                G.add_node(new_comp.id, activity=new_comp.activity)
                G.add_edge(branch.id, new_comp.id, genneration=generation)
                if new_comp.activity > 7:
                    new_heads.append(new_comp)
        heads = new_heads
        heads.sort(key=lambda x:x.activity)
        heads.reverse()

    nx.write_gml(G,"test.gml")