python2.5,R-2.7.1,rpy1.0.3という組み合わせで。rpy2-alphaはコンパイルが通らなかった。

Rのインストール

Rはソースからコンパイル。enable-R-shlibオプションをつけるのを忘れずに、というより配布されていたバイナリはこれつきでコンパイルされてないようなのでそのままではrpyが入らないはず。

R CMD config --ldflags

とか打つとOKかどうかわかる。　ダメげな場合は潔くソースから。

./configure --enable-R-shlib
make

ここでまた、make installするとshlibがなんかみつからなくなるという現象に悩まされた。そのため、ユーザー領域でコンパイルしたバイナリを使う事にした。/Users/kzfm/R-2.7.1/bin/Rが実体で/opt/local/bin/Rにシンボリックリンクを張って使っている。

Rpyのインストール

export R_HOME=/Users/kzfm/R-2.7.1
export DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH:$R_HOME/lib

これが重要。あとmacのgccでもコンパイルできなくてmacportのgcc使ってコンパイルする必要があるかもしれない。

あとは普通に

sudo python setup.py install

でいけるはず。

XYZでreadした分子(CC)もAddHydrogensできない

>>> from openbabel import *
>>> conv = OBConversion()
>>> conv.SetInFormat("xyz")
True
>>> mol = OBMol()
>>> conv.ReadFile(mol, "test.xyz")
True
>>> mol.NumAtoms()
2  
>>> mol.AddHydrogens()
True
>>> mol.NumAtoms()
2

でも分子中のそれぞれの原子を指定して

mol.AddHydrogens(atom)

としてやればきちんと水素が付加されることがわかったので、とりあえずスクリプト中でループをまわせばいいかな

pythonラッパーで。openbabelで気をつけないといけないところはAtomのインデックスは1から始まるのに対しBondのインデックスは0からはじまる。

from openbabel import *

mol = OBMol()

a = mol.NewAtom()
a.SetAtomicNum(6)
a.SetVector(0.0, 1.0, 2.0)

b = mol.NewAtom()
b.SetAtomicNum(6)
b.SetVector(1.0, 2.0, 2.0)

mol.AddBond(1, 2, 2) 
mol.SetTitle("Ethene")

obConversion = OBConversion()
obConversion.SetOutFormat("smi")

smiles_string = obConversion.WriteString(mol)

print smiles_string

あとAddHydrogensで水素を付与した場合の水素原子はSetVectorメソッドで座標が貼付けられないの。よくわからん。サンプルだときちんと水素がついてるので、なんかの属性をセットしとく必要があるのではなかろうか。

今日は休日だったので子供の面倒をみつつ、Project Eulerをはじめてみた。pythonで解いてる。今日一日で17問程解いたらjapan ranking 200位丁度になった。

ちょっと悩んだのはproblem 10で200万以下の素数の総和を求めるという問題。

最初は素数リストで割り算してみて割り切れるものがなかったら素数リストにappendという単純な方法でやってみたけど、10万くらいで遅くなってしまった。

というわけで篩で。200万までの奇数リストをつくって篩にかけていく。終了条件は最も大きい素数の二乗が200万を超えるかどうか。 で、できた素数のリストをreduceで畳んでprint

numlist = range(3,2000000,2)

def sieve(nums):
    prime_list = [2]
    while 1:
        prime = nums[0]
        prime_list.append(prime)
        if prime * prime > max(nums):
            prime_list = prime_list + nums[1:]
            break
        else:
            nums = [x for x in nums if x%prime != 0]
    return prime_list

print reduce(lambda x,y: x+y, sieve(numlist))

Project Eulerを解いていくには篩とか素因数分解とか組み合わせ発生用の関数を良く使う。

permutation用のジェネレータを見つけた。

def permutations(L):
    if len(L) == 1:
        yield [L[0]]
    elif len(L) >= 2:
        (a, b) = (L[0:1], L[1:])
        for p in permutations(b):
            for i in range(len(p)+1):
                yield p[:i] + a + p[i:]

再帰の中でyieldが使われている場合の動きのイメージがまだピンとこないなぁ。

それから、インデックスの範囲外は空リストが返ってくるからlispとかhaskellっぽく書けるのか。

>>> a = [1,2,3,4]
>>> a[1:]
[2, 3, 4]
>>> a[2:]
[3, 4]
>>> a[3:]
[4]
>>> a[4:]
[]
>>> a[5:]
[]
>>> a[:0]
[]

大きい数字の下位10桁を求める。桁溢れとかそんなんで数字があわなくなるので、欲しい桁以上の数字は捨てる。48(1^1+2^2...1000^1000)と97(大きい素数)はこれで解ける。

Problem 48

sum = 0
for i in range(1,1000):
    sum += reduce(lambda x,y: x*y%10000000000,[i] * i)

print sum % 10000000000

Problem 97

def fact2():
    num = 1
    for i in range(7830457):
        num = num * 2 % 10000000000
    return num

print (1 + 28433 * fact2() ) % 10000000000

25問解いてレベル1

多面体ゲット

pythonのyieldって呼ばれるとメインに戻ってくるようなイメージがあったので再帰の終了条件でyieldがあると気持ち悪かったのだけど、呼び出し側に返されるらしい。

というわけで、やっと理解したので、再帰でyieldをガンガンつかう。

あと、状態の凍結と継続の違いがちょっとあやふやかも。

def permute(item):
    if len(item) == 1: yield item
    for i in range(len(item)):
        for x in  permute(item[:i]+item[i+1:]):
            yield item[i] + x

c=0
for i in permute("0123456789"):
    c += 1
    if c == 1000000:
        print i
        break

mixiでpubchemからSMILESを抜き出すのは？みたいなエントリがたってたのだけど消えちゃったみたい。

なんか、ここは宿題まるなげちゃうわーみたいな厳しいコメントついてたからかな。まぁああいったコメント書くヒトはすっきりするのかもしれんけど、見るほうにとっては情報量ゼロのゴミなんだよなー。まだ、課題まるなげっていう情報のほうが情報量的に有意義。

ちょっと考えたのでここに書いておく。

XMLをBeautifulSoupで

import urllib2,sys
from BeautifulSoup import BeautifulSoup

cid = sys.argv[1]
url = 'http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid=%s&disopt=DisplayXML' % cid

opener = urllib2.build_opener()
html = opener.open(url).read()
soup = BeautifulSoup(html)

print soup.findAll('pc-infodata_value_sval')[-2].string 

xmlがなんかfindしにくいので配列の要素指定に-2とかやって良くない香り。この程度のだったらxmlじゃなくてSDFを正規表現でいじるな。

import urllib2,sys,re

cid = sys.argv[1]
url = 'http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid=%s&disopt=DisplaySDF' % cid

html = urllib2.build_opener().open(url).read()
p = re.compile('<PUBCHEM_OPENEYE_CAN_SMILES>\n(.+)\n')
m = p.search(html)
print  m.group(1)

biopythonとかでもいけそうな気がするし、urllib2+openbabelの組み合わせでも良いかもしれない。

1からnまでの連続した数をおのおの一つだけ桁にもつ素数のうち最大のものを探せ

987654321から数を減らしながら素直に見つけていくという安直なコードを書いてサブミット。

import math

def make_primes(n):
    primes = [2]
    nums = range(3,n+1,2)
    while 1:
        prime = nums[0]
        if prime ** 2 > n:
            primes += nums
            return primes
        else:
            primes.append(prime)
            nums = [x for x in nums if x%prime != 0]
    return list(set(primes))

def permute(item):
    if len(item) == 1: yield item
    for i in range(len(item)):
        for x in  permute(item[:i]+item[i+1:]):
            yield item[i] + x

def is_prime(n):
    for j in primes:
        if j*j > n: return True
        if n%j == 0: return False
    return True

if __name__ == "__main__":
    num = 7654321
    primes = make_primes(int(math.sqrt(num)))

    for i in permute(str(num)):
        x = int(i)
        if x%2 != 0:
            if is_prime(x):
                print x
                break

それよりもみずぴー日記のコメントが参考になった。

「各桁の数を足した和が3の倍数ならその数自身も3の倍数」って性質から、素数のPanditital数はあるとすれば1桁か4桁か7桁に絞られるんじゃないでしょうか。

えーなんでなんでー？？？？

３の倍数のANo3をみてなるほどなーと感心した。