macbookにOpenCVを入れると楽しいことに気づいた

portでインストールすればよい

sudo port opencv

途中でこけたので、メッセージ通りに

Error: Some libs are missing from your X11 installation. Please run this command:
Error: sudo ln -s libXrandr.2.dylib /usr/X11/lib/libXrandr.2.0.0.dylib
Error: Target org.macports.fetch returned: missing /usr/X11/lib/libXrandr.2.0.0.dylib

あとサンプルがみつからなかったので結局ソースを落としてきてsample/pythonのスクリプトで遊んでみた。

motempl.pyってのが面白い。

motempl

pythonのサンプルが結構あるので、暇を見て読んでみる事にする。

ProductName Learning OpenCV
Gary Bradski
Oreilly & Associates Inc / ?円 ( 2008-10-03 )


これもやっぱ欲しくなってきた。

WassrTerm

職場のプロキシーがあれで、ブラウザからだとリロードしまくらないとちゃんと表示してくれないので、TwitTermベースでつくってみた。

urllib2でBasic認証のさせ方を覚えた。

#!/usr/bin/env python
# -*- encoding: utf-8 -*-

username = "user"
password = "password"
interval = 240

from datetime import datetime, timedelta
import os
import sys
import re
import time
import types
from threading import Thread
import urllib
import urllib2
import webbrowser
from xml.dom.minidom import parse
import thread
import warnings
warnings.filterwarnings("ignore")

try:
  import readline
except ImportError:
  sys.stderr.write('INFO: For better line editing capability, '
      'install "readline" module.\n')

passman = urllib2.HTTPPasswordMgrWithDefaultRealm()
passman.add_password(None, 'http://api.wassr.jp/', username, password)
authhandler = urllib2.HTTPBasicAuthHandler(passman)
opener = urllib2.build_opener(authhandler)
urllib2.install_opener(opener)

tagText = lambda node, tagName: \
  node.getElementsByTagName(tagName)[0].firstChild.nodeValue

class TwitterReader(Thread):
  target_url ='http://api.wassr.jp/statuses/friends_timeline.rss'
  def __init__(self):
    Thread.__init__(self)
    self.last_id = 0

  def run(self):
    while not TwitterWriter.already_exit:
      req = urllib2.Request(self.target_url)
      req.add_data("since_id=" + str(self.last_id))
      e = None
      try:
        e = parse(file=urllib2.urlopen(req))
      except Exception, ex: # HTTPError or ExpatError
        sys.stderr.write(str(ex) + " will retry after %s sec\n" % interval)

      if type(e) != types.NoneType:
        for status in reversed(e.getElementsByTagName("item")):
          screen_name = tagText(status, "author")
          text = tagText(status, "description")

          print_status(screen_name, text,
                       time_created_at(tagText(status, "dcterms:modified")))

      for i in range(interval):
        if TwitterWriter.already_exit:
          break
        time.sleep(1)


def print_status(screen_name, text, (rel_time, created_at)):
  print "[%s] %s (%s)" % (screen_name, text, rel_time)

plural = lambda n: n > 1 and "s" or ""

def time_created_at(s):
  """
  recieving text element of 'created_at' in the response of Twitter API,
  returns relative time string from now.
  """

  try:
    date = time.strptime(s, "%Y-%m-%dT%H:%M:%S+09:00")[:-2]
  except ValueError:
    return "", ""
  created_at = datetime(*date)
  d = datetime.now() - created_at

  if d.days:
    rel_time = "%s days ago" % d.days
  elif d.seconds > 3600:
    hours = d.seconds / 3600
    rel_time = "%s hour%s ago" % (hours, plural(hours))
  elif 60 <= d.seconds < 3600:
    minutes = d.seconds / 60
    rel_time = "%s minute%s ago" % (minutes, plural(minutes))
  elif 30 < d.seconds < 60:
    rel_time = "less than a minute ago"
  else:
    rel_time = "less than %s second%s ago" % (d.seconds, plural(d.seconds))

  return  rel_time, created_at.strftime("%H:%M:%S")

class TwitterWriter(Thread):
  target_url = "http://twitter.com/statuses/update.xml"
  already_exit = False
  def __init__(self):
    Thread.__init__(self)

  def run(self):
    while 1:
      try:
        update_text = raw_input().strip()
        if hasattr(sys, "winver"):
          update_text = update_text.decode("shift_jis", "replace").encode("utf-8")
      except EOFError: # EOFError
        TwitterWriter.already_exit = True
        break
      if len(update_text):
        params = {}
        params['status'] = update_text       
        up = urllib2.urlopen('http://api.wassr.jp/statuses/update.json', urllib.urlencode(params))
        if up.code == 200:
          print "(update) [%s] %s " % (username,update_text)

    TwitterWriter.already_exit = True

def main():
  if username == "" or password == "":
    print "Error: Please specify your twitter account info in", sys.argv[0]
    sys.exit(1)

  print "Note: CTRL-D (or CTRL-C) to exit"
  reader = TwitterReader()
  writer = TwitterWriter()

  reader.start()
  writer.start()

  reader.join()
  writer.join()

if __name__ == '__main__':
  main()

まだリプライとかできないけど快適だ。

macbookでRpyを利用できるようにしたのでメモ

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

でいけるはず。

openbabelのAddHydrogens

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)

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

openbabelで分子を構築するサンプル

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をはじめた

今日は休日だったので子供の面倒をみつつ、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))

pythonでpermutation

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]
[]

Project Euler Problem 48&97

大きい数字の下位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

Project Euler lv.1

25問解いてレベル1

euler lv.1

多面体ゲット

pythonのyield (project euler problem 24)

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