ラベル python の投稿を表示しています。 すべての投稿を表示
ラベル python の投稿を表示しています。 すべての投稿を表示

2011/05/16

Google Art Projectを利用して高画質の絵画を手に入れる

"The Starry Night" by Vincent van Gogh
2011年2月2日、GoogleさんはGoogle Art Projectというサイトをリリースしました。どういうサイトかっていうとMoMAニューヨーク近代美術館Tate Britainなど、世界中の美術館にある有名な絵画を渡り歩くことができるという趣旨のサイトで、誰もが一度は見たことがある絵画がー家に居ながらー非常にリアリスティックな解像度で見ることができる。恐ろしい時代になったもんです。

ところがこのサイト、高解像度で見られるのは非常にいいことなのですが、肝心のダウンロードができない。せっかく高解像度なんですから保存できてもいいはずなのに、だけどできない。Google Mapなんかは常に最新の情報に更新しなければならないので、ダウンロードに制限をかけているのは合理性があるんですけど、絵画とか全く更新する必要がないでしょう。著作権も切れてるのに。よく分からないです。

ということで、なければ自分でなんとかしようという精神で、1時間くらいでちゃっちゃとPythonでスクリプト組んでダウンロードするようにしてみました。即興で組んだのでいろいろ粗い面はありますけど、そこら辺は勘弁してください。

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

import sys
import re
import urllib2
import random
from PIL import Image

def usage():
    print "%s [URL]" % sys.argv[0]
    return 1

def get_id(url):
    response = urllib2.urlopen(url)
    if response.code != 200:
        raise TypeError
    data = response.read()
    return re.search(
        r'<div id="microscope" data-microId="(.+)">', data).group(1)

def get_filename(x, y):
    return "%i_%i.jpg" % (x, y)

def download_image(target_id, x, y, z):
    urllist = [
        "http://lh3.ggpht.com/",
        "http://lh4.ggpht.com/",
        "http://lh5.ggpht.com/",
        "http://lh6.ggpht.com/"
    ]
    base_url = random.choice(urllist)
    url = "%s%s=x%i-y%i-z%i" % (base_url, target_id, x, y, z)

    try:
        response = urllib2.urlopen(url)
    except urllib2.HTTPError:
        return False
    with file(get_filename(x, y), "wb") as fp:
        fp.write(response.read())
    return True

def download_and_get_max_x(target_id, x, z):
    result = download_image(target_id, x, 0, z)
    if result:
        return download_and_get_max_x(target_id, x+1, z)
    else:
        return x

def download_and_get_max_y(target_id, y, z):
    result = download_image(target_id, 0, y, z)
    if result:
        return download_and_get_max_y(target_id, y+1, z)
    else:
        return y

def combine_images(name, max_x, max_y):
    x_size = 512
    y_size = 512

    result = Image.new("RGB", (max_x*x_size, max_y*y_size))
    for y in range(max_y):
        for x in range(max_x):
            image = Image.open(get_filename(x, y))
            result.paste(image, (x*x_size, y*y_size))
    result.save(name)

def main():
    if len(sys.argv) == 1:
        return usage()
    depth = 3
    target_addr = sys.argv[1]
    target_id = get_id(target_addr)
    max_x = download_and_get_max_x(target_id, 0, depth)
    max_y = download_and_get_max_y(target_id, 1, depth)
    for y in range(1, max_y):
        for x in range(1, max_x):
            download_image(target_id, x, y, depth)
    combine_images("result.png", max_x, max_y)

if __name__ == "__main__":
    sys.exit(main())

使い方は簡単で、
$ python artproject.py http://www.googleartproject.com/museums/tate/mariana-248
のようにURLを指定してあげると勝手に取得して結合してくれます。こんな感じで:
"Mariana" by John Everett Millais
あぁそれにしてもお美しい……巷ではオフィーリアちゃんが有名だけど、やっぱりマリアナさんのほうが綺麗だよ。絶対。

余談
  • もともとの動機としては高解像度な絵画の壁紙が欲しくて、いろいろと探し回ったわけですが、どうもGoogle Art Projectしか高解像度の絵画が載っていないので仕方なくといった次第です。
  • 利用はどうやらGoogle's Terms of Serviceに準じているようです。ざっと見した感じだとダウンロードを禁止する条項はないですし著作権も切れているので全く問題ないように思えるんですが、もし警告とかきたらチキンですから遠慮無く取り下げます。ご了承ください。
  • このネタを題材に時間があればHTMLのスクレイピング、パケット解析についての入門記事を書いてみたいです。
  • なんだかんだ言ってちゃんと実物見たほうが綺麗。ミレイの作品は以前東京に行ったときにBunkamuraで見ました。Bunkamuraはそういう催し物よくやっていて、そこらへんが東京羨ましかったり

2011/04/19

並列環境におけるreduceとscanアルゴリズム

"fireking" by cmyk1219
1. 概要
reducescanはGPGPUなどの並列環境下において重要な役割を果たす関数だけど、あまりこれらの関数、特にscanについて言及している記事はあまり見かけない。なので今回はreducescanについて調べた結果をまとめていきたいと思う。

2. reduce, scan関数の概要
2.1. reduce
scanについてはreduceが分かればより理解しやすくなるので、ひとまずはreduceから始めることにする。

reduceとは、端的に言うと、対象のリストを与えられた二項演算子で纏め上げる関数だ。とは言っても、いきなりそんなことを言われてもよく分からない。まずは例を見ていこう:
>>> import operator
>>> reduce(operator.add, [1, 2, 3, 4, 5])
15
>>> reduce(operator.sub, [1, 2, 3, 4, 5])
-13
>>> reduce(operator.mul, [1, 2, 3, 4, 5])
120
reduceにおいて重要なのは二項演算子であり、上記の場合、reduceは以下のような演算を行う。
(((1 + 2) + 3) + 4) + 5 = 15
(((1 - 2) - 3) - 4) - 5 = -13
(((1 * 2) * 3) * 4) * 5 = 120
これにより、最初の演算はsum関数、3番目の関数はprod関数と等価であることが分かる。
リストの型は何も数値に限られている訳ではない:
>>> reduce(operator.and_, [True, False, True])
False
>>> reduce(operator.or_, [True, False, True])
True
これより、上記2つの演算はall, any関数と機能的に等価であることが分かるのだけれど、all, any関数のほうが効率的なので、実際はちゃんとall, any関数を使うべきだ。

他にもreduceを用いて様々な関数、アルゴリズムを記述できるが、あまりにも多いので今回は後回しにして、scanを重点的に解説していく。

2.2. scan
scan関数はreduceと比べて文献が少ない。言及している文章もあまり見かけないが、並列処理においてscanは重要な関数だ。残念ながらpythonで実装されていないけれど、仮に実装されていたとしたら以下のような働きをする:
>>> scan(operator.add, [1, 2, 3, 4, 5])
[1, 3, 6, 10, 15]
>>> scan(operator.sub, [1, 2, 3, 4, 5])
[1, -1, -4, -8, -13]
>>> scan(operator.mul, [1, 2, 3, 4, 5])
[1, 2, 6, 24, 120]
きちんと書き下せばどのような働きをしているのか分かりやすい:
[1, 1+2, (1+2)+3, ((1+2)+3)+4, (((1+2)+3)+4)+5]
結局のところ、scanは先頭からreduceしていった結果をリストにまとめ、それを返す関数となる。scanは別名prefix sumとも呼ばれている。
ちなみに、haskellではscanl1(scanl)という名前で実装されている。
Prelude> :type scanl1
scanl1 :: (a -> a -> a) -> [a] -> [a]
Prelude> scanl1 (+) [1,2,3,4,5]
[1,3,6,10,15]
また、CUDAライブラリであるthrustにはinclusive_scan(exclusive_scan)という名前で実装されている。
#include <thrust/scan.h>
int data[6] = {1, 0, 2, 2, 1, 3};
thrust::inclusive_scan(data, data + 6, data); // in-place scan
// data is now {1, 1, 3, 5, 6, 9}

3. 並列による制限
何故並列下においてreduce, scanが重要になってくるのか、それは、両方ともメニーコアな環境において効率的に動作するアルゴリズムが複数提唱されているからだ。reducescanは並列と相性がいい。つまり、既存のアルゴリズムをどうにかしてreducescanに落とし込むことができれば、並列化の恩恵を受けれられると。

ただし、並列環境においては、演算にいくつかの制限が生じてくる。具体的には、集合Gと演算子"・"の組において、
  • 演算子"・"は交換則"a・b = b・a"を満たしていなければならない
  • 演算子"・"は結合則"(a・b)・c = a・(b・c)"を満たしていなければならない
後者に関してはあたり前の話で、これが成り立っていなければ並列に計算を行うことができない。前者については多少曖昧で、実際には満たさなくてもよいアルゴリズムも作ることはできるけど、コアレッシングなどの問題を考えると満たしていることが強く求められる。実際、thrustのreduceは交換則を満たしていることを前提としている。

reduceの場合はこの2つの制限で十分だけれど、scanの場合は上記2つの制限に加えてさらに2つの制限を満たしておくと都合がいい。具体的には、集合Gと演算子"・"の組において、
  • ただ一つの単位元"e"が存在する
  • 演算子"・"は任意の元"a"に対して"a・a-1 = a-1・a = e"を満たす、ただ一つの逆元"a-1"が存在する
これら4つの制限を満たしている集合Gと演算子"・"の組を『アーベル群』という。要するに、scan関数は、対象の型から成る集合と二項演算子の組がアーベル群を成していれば都合がいいということになる。

具体例を出そう。整数と和の組(Integer, (+))は明らかに上記4つの制限を満たす。よって、この組は並列にreduce, scanオペレーションを適用でき、かつ都合が良い。
1 + 2 == 2 + 1
(1 + 2) + 3 == 1 + (2 + 3)
1 + (-1) == (-1) + 1 == 0

真偽値と論理積の組(Bool, (&&))は上記2つの制限を満たすため並列にreducescanオペレーションを実行可能だが、後半2つの制限は満たしていないため都合は良くない。

3次元の回転行列と積の組(Matrix, (*))は交換則を満たしていないため―少なくともthrustでは―並列にreduce, scanオペレーションを実行することはできない。

3.1. 都合が良いとは?
とは言っても、一体アーベル群であればどこがどう都合がいいのか?それは、リスト中において局所的にreduceを行いたい場合に役に立つ。

リスト[a0, ... , an]にscanオペレーションを実行して

を手に入れたとする。この場合、i < jにおいてbj+bi-1を計算すると、

となる。つまり、i+1からjまでの要素についてのreduceを僅かO(1)で求めることができる。

これには大きな利点がある。例えば、大規模な時系列データを、その周りのデータから平均化(平滑化)したいとする。各要素ごとに計算を行うのは明らかに非効率であるので、まずscanオペレーションを(+)に対して適用し、その後で差分を計算し、正規化を行ってやればよい。結局、scan (+)のやっていることはf(x)の不定積分F(x)を求めることと本質的に等価となる。この考え方はCVでは、積分画像として様々な箇所で用いられる。

そしてさらに、この操作はアーベル群であれば何でも良い……とは言っても、アーベル群であることは結構厳しい制限となる。

例を出すと、リスト中の任意の範囲ーiからjの成分がすべてある条件式に適合しているかどうか確かめたいとする。
関数型の考え方だと、まずはmap関数を用いて(C++でならtransform_iteratorを用いて)リストの成分をBoolに変え、その後でscanを利用することを思いつくかもしれない。しかし、残念ながら(Bool, (&&))はアーベル群でないので、そのまま愚直に適用することができない。
Prelude> scanl1 (&&) $ map (<3) [1,2,3,2,1]
[True,True,False,False,False]
-- 最後の2要素については条件を満たしているのに、Falseで埋もれてしまって分からない!
この場合は、単純にアーベル群である(Integer, (+))scanに適用してやればよい。
Prelude> scanl1 (+) $ map (\x->if x<3 then 1 else 0) [1,2,3,2,1]
[1,2,2,3,4]
-- きちんと第3要素のみ条件から外れていることが分かる
Pythonっぽく書くとたぶんこんな感じ:
>>> import operator
>>> f = lambda x: 1 if x < 3 else 0
>>> scan(operator.add, [f(x) for x in [1,2,3,2,1]])

4. より詳しく知りたい方は
scanアルゴリズムのGPGPU実装に関する詳細については、現在のところ『Parallel Prefix Sum (Scan) with CUDA(Mark Harris, NVIDIA)』が一番分かりやすい。より詳しく知りたい方はそちらをどうぞ。

というより、基本NVIDIAの出しているテキストは分かりやすい。すごいことです。

2010/09/28

マルコフ連鎖モンテカルロ法を実装してみよう

"Chain Bridge, Budapest" by szeke

約2ヶ月ぶりということで、ご無沙汰しております。書きたいネタというのは結構あって書こう書こうとは思っているんですが、なにより書くとなるとあれもこれも伝えなきゃいけないみたいな思いになって、結局分量の多さから諦めてしまうというのが結構続いています。もう少し気を張らずに更新していきたいものです。

さて、最近の自分はマルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo, MCMC)についておりました。何にも知らない方ににとってはよく分からないかもしれませんが、この手法はマルコフ連鎖が持つ簡明さとモンテカルロ法が持つ実用性が合わさった凄まじい手法なんです。そしてなによりとてもエレガント。

とりあえず知らない人のためにてきとーに解説しますと、与えられた適当な関数から確率変数をサンプリングするための公式です。べつにこれじゃなくても確率変数のサンプリングにはいろんなアルゴリズムがあるわけですが、これを使うと関数の計算コストがやけに高い場合だとか、ピーク部分を集中的にサンプリングしてくれたりとかで結構便利で、金融や工学、物理の分野でも非常に役立っております。計算物理をやる人にとっては必須に近い手法といっても過言ではありません。

そんなMCMCですが、実際に提案されている手法には様々なものがあります。不変分布から遷移確率を逆算するため、自由度が増えるから遷移確率は一意じゃなくなってしまうんです。今回は最も有名なMetropolis-Hastingsアルゴリズム(*1)を用いて、任意の関数に従って乱数のサンプリングを行うコードを練習がてらpythonで書いてみました。

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

import sys
import math
import random

def main():
   f = lambda x: math.exp(-x*x / 10.0)
   for x, fx in metropolis(f):
       print x, fx

def metropolis(func, start=0.0, delta=2.0, burn_in=100, seed=1):
   metro_iter = _metropolis(func, start, delta, seed)
   for i, value in enumerate(metro_iter):
       if i > burn_in:
           yield value
           break
   for value in metro_iter:
       yield value

def _metropolis(func, start, delta, seed):
   random.seed(seed)
   initial_x = start
   initial_fx = func(initial_x)
   proposed_x, proposed_fx = 0.0, 0.0
   while True:
       proposed_x = initial_x + random.uniform(-delta, delta)
       proposed_fx = func(proposed_x)
       ratio = proposed_fx / initial_fx
       if ratio > 1.0 or ratio > random.random():
           initial_x, initial_fx = proposed_x, proposed_fx
           yield proposed_x, proposed_fx
       else:
           yield initial_x, initial_fx

if __name__ == '__main__':
   sys.exit(main())
うーん・・・きれいだ(*2)・・・

こいつのすごいところは、アルゴリズム自体は非常にシンプルで実装しやすいというところにあります。行数も上を見ると分かるのですがとっても少ないし、コア部分はたった25行くらいで構成されています。

なのに使っている学問は結構複雑で、定常分布から遷移行列を導出するために詳細釣り合いの条件からMetropolis的制約を用いて導出、サンプリングを行っております。このアルゴリズム自体は結構知ってる人が多いかもしれませんが、実際に理論方面で理解している人はどうなんでしょう。

んで、正規分布のガウス関数からサンプリングを行った結果がこれ。


赤線が元の確率分布。きちんとサンプリングが行われているようです。今回は正規分布で計算を行いましたが、非負の関数であればなんでも構いません。

あとはこれを研究テーマへと落とし込むことができれば、とりあえず一区切りつけるかなって感じです。

*1: 提案確率関数が対称であるという制約を持つのがMetropolis。それを非対称の場合にも拡張したのがHastings。個人的にはMetropolis-HastingよりMetropolisのほうが綺麗で好きです。というか提案確率密度関数が非対称じゃないといけないようなシチュエーションがよく分からない
*2: yieldとジェネレータって便利ですよね。計算速度も効率的だし、慣れてくるとすべて遅延評価で行わせたいくらい

2010/04/06

マルコフ連鎖を使って簡単な問題を解いてみよう

"Cute Gyoza / Meat Dumpling Key Chain" by ♥ Rainbowcatz ♥

はじめに

人工無能、一時期流行りましたねー。

難そうなわりに以外と簡単で、しかも結果が面白いということでRubyやなんかで実装している記事をよく見かけました。そういった中でマルコフ連鎖が注目されていた時期が一時期ありましたが、正直あれはマルコフ連鎖を利用しているというよりは、別のアルゴリズムにマルコフ連鎖という名前をつけたというか、あんまり関連性がないなーという印象でした。

マルコフ連鎖の本質はそういうことではなく、状態の遷移確率を行列によって求めたり、定常分布を固有値問題で計算するところにある、と自分は考えています。マルコフ連鎖についてのより詳しい説明は他のブログにお任せするとして(自分がやるとすぐにボロが出てしまう)、今回はマルコフ連鎖を使って簡単な問題を解き、実際にマルコフ連鎖がどんな感じなのか確かめてみることにします。

もんだい

二つの箱A,Bが用意されている。2つの箱の内部は完全に遮蔽されており、外部からは確認することができない。また、Aの箱には赤玉がN個、Bの箱には白玉がN個入っている。

ここで、両方の箱に腕を突っ込み、それぞれの箱から玉を一個づつ取り出し、交換し、反対側の箱のなかに入れることにしよう。今回の場合、Aの箱には赤玉がN-1、白玉が1個入ることになる。もちろん、Bの箱には赤玉が1個、白玉がN-1個入っている。

この試行をひたすら続けて、そう、例えば無限回行ってみよう。この場合、Aの箱に入っている赤玉の数はどれくらいであろうか?個数に対応した確率分布を求めよ。

回答

1回、2回ならまだしも、無限回なんて計算できるかクソが!そんなふうに考えていた時期が俺にもありました。

まず、Aの箱における赤玉の個数をmとおきます。すると、これは状態mの遷移問題と考えることができ、マルコフ連鎖が使えます。

また、1期間後にとりうることのできるmの遷移はm-1, m, m+1の3通りしか存在しません。以上を利用して、各mにおける状態の遷移確率を求めることができます。

まず、mがm+1に遷移する確率(m->m+1)を求めてみましょう。この場合、確率はAの箱の白玉が選ばれる確率とBの箱の赤玉が選ばれる確率の積に等しいので、

と表せます。

同様にして、m->m, m->m-1 についても表していきます。
  • mがmに遷移する確率は(i)Aの箱の赤玉が選ばれる確率とBの箱の赤玉が選ばれる確率の積と(ii)Aの箱の白玉が選ばれる確率とBの箱の白玉が選ばれる確率の積 の和((i)+(ii))に等しい。
  • mがm-1に遷移する確率は確率はAの箱の赤玉が選ばれる確率とBの箱の白玉が選ばれる確率に等しい。
以上の議論をまとめると、i->jへの遷移行列pijは以下のように表せます。

ここで、遷移行列はN×Nの次元を持っているものとします。

以上を利用して遷移行列を求め、pythonを用いて定常分布を算出しました。プログラムのソースコードを以下に示します。
#!/usr/bin/env python
#-*- coding:utf-8 -*-
"""
eig.py [N]

定常確率分布を記述したファイルを新しく書き込むプログラムです。
引数には箱に入っている玉の数を指定します。引数が指定されなかった場合には、N=10が用いられます。
"""

import sys

def main():
    number = int(sys.argv[1]) if len(sys.argv) == 2 else 10
    P = getStateMatrix(number)
    result = getFiniteStateSpace(P)
    fp = file("%i.dat" % number, "w")
    fp.write("\n".join([str(abs(num)) for num in result]))
    fp.close()

def getStateMatrix(n):
    """n個における遷移確率行列を返します"""
    from scipy import matrix
    P = [[getState(i, j, n) for j in xrange(n+1)] for i in xrange(n+1)]
    return matrix(P)

def getState(i, j, N):
    if   j == i-1: return float(i*i)/(N*N)
    elif j == i:   return float(2*i*(N-i))/(N*N)
    elif j == i+1: return float((N-i)**2)/(N*N)
    else: return 0.0

def getFiniteStateSpace(mat):
    """
    i->jに遷移する確率を行列(i,j)で表現した遷移確率行列から、定常分布を求めます。
    この関数は行列を転置する必要がない点に注意してください。
    
    mat : N×Nの遷移確率行列(numpyのmatrixでなければならない)
    戻り値 : N次元の定常分布ベクトル
    """
    from scipy import linalg
    la, v = linalg.eig(mat.T)
    result = zip(la, v.T)
    
    # 定常分布を求めるため、固有値が1に近い順で固有ベクトルをソートする
    result_sorted = sorted(result,
                           cmp=lambda x, y: cmp(abs(1.0-x), abs(1.0-y)),
                           key=lambda x: x[0])
    
    eig_value   = result_sorted[0][0]
    # 固有値が1近傍でなかったら例外を送出
    assert abs(eig_value - 1.0) < 1.0e-6, "the eigen value of FSS is apart from 1.0"
    
    fss = result_sorted[0][1]    
    return fss/sum(fss) # 正規化

if __name__ == "__main__":
 sys.exit(main())

実行は対象のソースコードを"eig.py"として保存してから、カレントディレクトリを適切なディレクトリに移動した後で、コマンドライン上から"python eig.py 10"のように指定してあげます。この場合、N=10(赤玉10個, 白玉10個)の際に遷移するであろう定常分布のデータが"10.dat"として書き出されます。

ここで、試しにN=6, 10, 50, 100における確率分布を載せておきます。
N = 6 のとき。離散的でよく分からないが、N=3付近の確率が一番高くなっている。
N = 10, 50 のとき。だんだん滑らかになっていく。
N = 100 のとき。N=50付近のピークが非常に鋭くなっている。

以上より、Nを増やせば増やすほど、N/2近傍により鋭いピークが表れるようになることが確認できます。

どーいうことが言いたいのか

結局のところ、お手玉の数を増やせば増やすほどmは中央付近で安定するわけです。「いや、最初に赤玉と白玉分けられているじゃん」それはそうなんですけど、今回の場合は無限回試行しているので、最初どのように分布してようが全く関係ないわけで。

最初は偏っていた分布がごちゃごちゃやっている内に中央付近に移動し、そこから多少のゆらぎはあるかもしれないが、結果的に中央へと帰ってしまうような、そういうイメージです。

まとめ

確率を無限回にすっ飛ばしたり、ある時点での確率分布を求めたり、モンテカルロ法にも使われているマルコフ連鎖はとっても賢くてすごいやつという記事なのでした

2010/02/08

pythonでunzip()関数を実現するためには

pythonには(当たり前だけど)リストをまとめるためのzip()関数が存在する。

>>> A = [1,2,3]
>>> B = [4,5,6]
>>> zip(A, B)
[(1, 4), (2, 5), (3, 6)]
これだけじゃちょっとありがたみが分かんないけど、配列っていうのは『オブジェクトjの状態』などを格納してある場合が結構あるので、jに関してループを回したいとかのような場合に便利だったりする。

for (state1, state2) in zip(A, B):
  # do something
で、これはいいんだけど、なぜかpythonにはzip関数でまとめられた配列を展開する"unzip()"関数が存在しない。なんでだろーとか思って調べたら、どーやらこれでできるらしい。

>>> C = zip(A, B)
>>> zip(*C)
[(1, 2, 3), (4, 5, 6)]
なるほどー…pythonらしい。

最初『…ですよねー!』ってなったけど、こーゆう解決案をぱっと思いつくのが賢い人間なんだろーなーとか思った。…うん、そりゃそうだ。

2009/10/07

(追記)郵便局で作る日本全国の光の地図

これは前回の記事
住所一覧からマックと吉野家の光の地図を作ってみた
の続きです。今回はデータとして、郵便局の住所一覧を用いました。
また、データは『全国郵便局名一覧』からダウンロードしました。ありがとうございます。


















当たり前ですけど、郵便局はやっぱり多いですね。解析にまる3,4日くらいかかったような気がします。
しかも日本全国に散らばっているあたり、郵政事業がいかに重要なのか思い知らされます。

2009/10/04

住所一覧からマックと吉野家の光の地図を作ってみた

おひさしぶりです。

お久しぶりです。どうやら巷ではアクセスアップのコツとして「内容よりも一日の更新頻度を上げたほうが人気が出る」らしいですが、ここはそんなこと全く気にせずにのんびりと更新していきます。(かと言って質が高いかと言えば疑問ですが…)

さて、更新してみようと思い立ったきっかけはこんな記事からでした。
「マックまでの近さ」が光で表された米国地図 | WIRED JAPAN
この地図を作成したのは、Steven Von Worley氏。ロサンゼルス盆地の、本当に何もないようなところで目にしたマクドナルドに刺激されてこの地図の作成を思い立ったという。
で、ブックマークのコメントを見たらこんなものが。
はてなブックマーク - 「マックまでの近さ」が光で表された米国地図 | WIRED VISION
Layzie food, ネタ 東と西でクッキリ分れてるねえ。誰か、日本で「吉野家」バージョンで作らないだろうか。
やってみようじゃないの。

というわけでちょっとした思いつきで吉野家地図の日本版を作ってみました。ちなみに本家ではきちんと距離を計算してそこからの最短距離で色を作っているらしいですが、そんなことしても面倒くさい割にそんな役に立たないので単純な加算処理にしてあります。

工程は後述するとして、画像は以下。

これが吉野家の光…

東京や大阪、京都の大都市は吉野家が密集していることが伺えます。また各種国道や高速自動車道をなぞるように点在しているようですね(*1)。しかし海岸沿いなどにはほとんどなく、東北地方、特に札幌を除く北海道は悲惨です。これが地方格差ってやつか。

ついでにマクドナルドの地図も作って見ました。結果の画像は以下。

やっぱり吉野家と比べると明るいですね。でも吉野家の光の密度が多くなったような感じで、全体的な傾向としてはあまり変わらないような気がします。
ちなみに2つとも画像映えを良くするために少しのグロー効果も入れてあります。

メイキング

作った流れとしては以下のような感じです。
  1. まずマックのすべての店舗の住所一覧を取得する。
  2. 住所から経度と緯度のペアに変換する。
  3. Proce55ingにデータを読み込ませ、パーティクルとして処理する。
店舗一覧はマックの場合は47都道府県すべてを検索し、電話番号などの不要なデータを消去し住所だけを抜き出すpythonスクリプトを組んで対応。吉野家は件数が少ないのでコピペと手動で対応しました。

住所から経度、緯度のペアに変換する手順では、今回はGoogleのウェブAPIサービスを利用しました。住所のデータを投げればXML形式で返ってくるので、ペアから経度と緯度部分だけを抜きだし、csv形式で保存するスクリプトを組みました(*2)。

最後にProce55ingを用いてデータを読み込ませ、経度、緯度を単なるXY座標のデータとして処理させました。ただし、座標軸の違いによりこのままでは日本が上下反対になってしまうので、単純に-1をかけて反転させてあります。
あとは見やすいように座標変換を行って、完成です。

ソースコード

今回のスクリプトをcodereposのリポジトリ上にアップしました。
svn checkout http://svn.coderepos.org/share/lang/java/misc/light
でチェックアウトしてください。

外部ライブラリに依存していないので単純にproce55ingを使えば動くと思います。ついでにpythonスクリプトも同様にアップしてみました(trans.py)。参考までに。
ちなみに使い方ですが、座標が書いてあるCSVファイルを指定した後はProce55ing側が勝手に読み取ってプロットしてくれます。

ドラッグで見たい場所を移動できます(Google Mapsと同じ)。
'f'キーでズームイン、'g'キーでズームアウトを行います。マウスポインタを中心としてズームを行います。
'j'キーでパーティクルサイズを大きく、'h'キーでパーティクルサイズを小さくします。
's'キーでフレームを保存します。

(ズームして分かる、東京都のマックの様子)

'b'キーで日本地図(Thanks http://www.freemap.jp/)をバックに表示します。

まぁ、おまけ機能です。なんで位置合わせしないかというと、そこまでするときちんと経度、緯度を座標変換しなきゃなんなくて面倒だからです。

終わりに

作業時間としては大体3日間くらいかかりました。そんなに難しくありませんが、いろんな知識が複合的に絡んでいたのでそれなりに面白かったです。
ちなみに、郵便局バージョンでも作れないかと今住所一覧を経度、緯度に変換しているんですがこれ3万件くらいあってすごい時間かかるんですよね…現在2日間くらい回してるんですがまだ1万件しか処理できていません。できたらちゃんと報告します。

*1 個人的に、国道にはよく吉野家が点在しているイメージがあります。
*2 もちろんこのままではサーバに過負荷(単なるDoS攻撃になってしまう)がかかってしまうので5秒のスリープを入れてあります。

追記
typo ×吉野屋 ○吉野家 

2009/07/15

Pythonistaのための2chライブラリ"twopy"


"Beautiful Kuwait 2" by creativesam

概要

さて、実は以前からちょっと2chの文章を利用して、いろいろ弄くることができないかなぁとは思っていたんですが、残念な事に2ch用のpythonライブラリなんてマニアックなものは作っている人は誰もいない。
それじゃあ自分でなんとかしようの精神で作ってみようと。正直使う人が居るのかどうかは甚だ疑問なライブラリですが、損するわけじゃないので一応公開してみます。

というわけで、Pythonistaのための2chライブラリ"twopy"をリリースしました。
現在のところ単純なスレッド一覧の取得、スレッド上のコメント取得、コメント書き込みや新規スレッド機能に対応しておりますが、まだ完全に動作確認しているわけではないので使用は自己責任で。ちなみにライセンスはMITライセンスなんで、比較的自由に扱えると思います。改良はこれから始まる期末テストに飽きたらやると思います。

インストール方法

インストールに必要な環境は以下の通りです。
Python 2.5 or higher
2.6.2での環境で動作を確認。外部ライブラリに依存していないので、もしかしたらそれより前のヴァージョンでも動作するかもしれません。

インストール方法はCodeReposに上げておきますので、

$ svn checkout http://svn.coderepos.org/share/lang/python/twopy/trunk

なりなんなりでチェックアウトしてください。

簡単な使い方

スレッド一覧の取得
import twopy
b = twopy.Board("http://takeshima.2ch.net/news4vip/")
b.retrieve()

スレッドごとのタイトル、レス数、速度(res/h)を取得することができます。
for t in b:
print u"%s (%i) %.4f" % (t.title, t.res, t.velocity)

配列のように呼び出すことも可能です。
print b[0].title
print b[0].since


スレッドの内容の取得
t = b[0]
t.retrieve()

で一番上にあるスレッドの内容を取得します。
t = twopy.Thread(b, "1247375844.dat")
t = twopy.Thread.initWithURL("http://takeshima.2ch.net/test/read.cgi/news4vip/1247375844/")

として、datファイルあるいはURLから直接初期化することも可能です。

コメントを出力するには以下のようにループさせます。また、特定のコメントのみを取り出すこともできます。
for c in t: print c.render()
print t[1].render()

ただし、その場合の引数は実際のレス番号を指定する必要があります(0から始まるわけではない)。

コメントの抽出
特定のコメントのURLやレスポンスを抜き出すことができます。
for c in t: print c.extractResponses(returnType="comment")
for c in t: print c.extractUrls()

追記: 10/07/22
returnType引数は廃止され、メソッドとして分割されることになりました。
for c in t:
    print c.extractResponses()
    print c.extractResponsesAsInteger()
    print c.extractResponsesAsComment()

コメントの書き込み
特定のスレッドに対して、コメントを書き込むことができます。

t.post(name=u"", mailaddr=u"", message=u"")

ただし、2ch側の仕様により必ずしも書き込まれるとは限らず、確認が求められる場合がありますので、実際には以下のようにします。
r = t.post(name=u"", mailaddr=u"", message=u"")
if r[0] == twopy.STATUS_COOKIE:
    t.post(name=u"", mailaddr=u"", message=u"", hidden=r[2])

すべての確認作業を自動的に行うメソッドも存在します。

t.autopost(name=u"", mailaddr=u"", message=u"")

ただし、上2つのコードは2ch側の確認をユーザに見せることなく書き込みを行っていますので、実際にはユーザーに文章を見せてから、同意を行なわせる必要があるでしょう。

他にもスレッド作成機能などのいろいろ細々したのもありますが、それはソースコード見ればだいたい理解できると思います。行数そんな多くないですし。
もし正常に動かないなどのバグがありましたらご連絡お願いします。暇があれば対処したいと思います。それでは。

2009/04/12

辞書と組み合わせと再帰関数


Photo by joka2000, Double-flowered Cherry Blossoms

辞書が与えられて、そのキーと値の群の中でタプルやリストを持っている値を、任意の値で組み合わせたリストを返すという処理はどう実装すべきか悩んでいました。
言葉にすると難しく聞こえますが、要は
{"a":(1,2), "b":4, "c":(3,4)}という辞書が与えられた場合、返すリストは
[{"a":1, "b":4, "c":3},
{"a":1, "b":4, "c":4},
{"a":2, "b":4, "c":3},
{"a":2, "b":4, "c":4}]
となるわけです。

通常の行列のような処理は単純にループを2回回せばいいわけですが、今回の問題の場合ループ回数は変動するわけで、単純にループを回すわけにはいかない。どうすればいいのかちょっと考えて見たら、再帰関数を使えば案外スマートに実装できそうな気がしたので、実際にやってみました。


from copy import copy

def makeKwargsList(kwargs, exceptions=[]):
"""
タプルやリストを含んだ辞書から、それらの値を含む任意の組み合わせを
すべてリストとして返します.
例:
kwargs = {"a":(1,2), "b":3, "c":(True,False), "d":(1,2)}
exceptions = "d"の場合、返される辞書のリストは、
[{"a":1, "b":3, "c":True, "d":(1,2)},
{"a":1, "b":3, "c":False, "d":(1,2)},
{"a":2, "b":3, "c":True, "d":(1,2)},
{"a":2, "b":3, "c":False, "d":(1,2)}] となります.

kwargs : タプルやリストを含んだ辞書
exceptions : 組み合わせとして出力してほしくないキーのリスト
"""
if type(kwargs) != dict: raise TypeError
# キーワードの抽出
kwdict = {}
base = {}
for (key, value) in kwargs.items():
if (type(value) == list or type(value) == tuple) and \
(key not in exceptions):
kwdict[key] = value
else: base[key] = value
# key, kw の分離
keys = kwdict.keys()
values = kwdict.values()
kwargs_list = []
__makeKw(values, [], kwargs_list)
# 結合
r = [dict(zip(base.keys() + keys, base.values() + ag)) for ag in kwargs_list]
return r

def __makeKw(kwlist=[], base=[], kwargs_list=[]):
argslist = kwlist[0]
if len(kwlist) > 1:
newKwlist = kwlist[1:]
for kw in argslist:
new_base = copy(base)
new_base.append(kw)
__makeKw(newKwlist, new_base, kwargs_list)
else:
for kw in argslist:
new_base = copy(base)
new_base.append(kw)
kwargs_list.append(new_base)

__makeKw()が再帰関数で、この場合[(1,2),(3,4)]のリストから、[[1,3], [1,4], [2,3], [2,4]]をkwargs_listに代入させていく関数となっています。

本来のプログラマがどうやってこのようなアルゴリズムを組んでいるのかは残念ながら勉強不足で分かりませんが、とりあえずきちんと動作しますし、そこまで速度が重要となるような処理ではないですので(せいぜい100〜10000程度)、べつにこれでいいかなと。

それにしても、1年前は思いつかずに挫折していたであろうアルゴリズムが実際に組めるようになってくると、ちょっと成長している感がありますね。

2009/03/20

MeCabを用いてスパムフィルタを作ってみよう


Photo by vsz, night glow

以前このブログでMeCabによる形態素解析を紹介しました。正直その後すっかり取り上げたことを忘れてのほほんと過ごしていたわけです(*1)が、ふとしたことでベイジアンフィルタに関するアルゴリズムの記事を見つけ、日本語でこのような記事があるなんて珍しいなということで、ちょっくら実装してみようと思い立ったわけです。

形態素解析部分はMeCabくんがやってくれるので、こっちがするのは名詞を抜き出してデータベース辞書を作り、Graham方式を用いて実装したくらいです。正直ただ単純に実装しただけなのでそこまで参考にならないと思いますが、一応coderepos上に公開してみます。

svn checkout http://svn.coderepos.org/share/lang/python/spam Somewhere

でチェックアウトしてください。

テストにつきましては正直自分のメールアドレスにスパムが届かなかったので(*2)、『Google ニュース』を使って『政治』分野をスパム、『スポーツ』分野を非スパム文章として40件サンプルを手作業で作り(*3)、2つのセンテンス

大阪府の橋下徹知事が進める、府庁を大阪市の第三セクター「大阪ワールドトレードセンタービルディング」(WTC)に移転する構想について、公明党府議団は19日、「原案には賛成できない」として府議会で審議中の移転関連予算案と条例案に反対する方針を決めた。 公明党は自民党とともに知事与党の立場だが「議会として十分な議論ができていない。あまりに拙速だ」としている。条例案の可決には出席議員の3分の2以上の賛成が必要で、原案通りの可決は極めて厳しい状況となった。橋下知事は同日、高齢の障害者らの医療費自己負担を引き上げる案などを、議会側の求めに応じて修正する意向を議会側に伝達。府庁移転問題の好転を図る狙いがあるとみられるが、先行きは不透明だ。橋下知事は、府議会総務常任委員会であらためて移転への意欲を示し、記者団に「(反対は)公明の政治的決断。自民、民主、諸派に納得していただくよう精いっぱい頑張る。可能性はゼロではない」と述べた。

フィギュアスケートの金妍児(キムヨナ)(韓国)が過去の国際大会の練習中、日本選手に進路を妨害されたなどとする韓国メディアの報道について、日本スケート連盟は19日、韓国スケート連盟に対し、報道内容についての調査を求める要望書を提出した。また、国際スケート連合(ISU)に対しても同日、報道への見解を問う文書を提出した。日本スケート連盟によると、韓国側からの直接の抗議や、ISUからの警告はこれまでないという。同連盟の常山正雄専務理事は「日本選手が意図的に妨害した事実はなく、今回の報道には大変困惑しているし、遺憾」としている。


がきちんと判定されているか確かめました。これらの記事を選んだ深い理由は得にないです。
結果は、

$ python test.py
センテンス1がスパムである確率 : 1.0
センテンス2がスパムである確率 : 6.07068879794e-100
スパム文章のサンプル数 : 20
非スパム文章のサンプル数 : 20
スパム辞書が学習した単語数 : 1040
非スパム辞書が学習した単語数 : 854

大体大丈夫そうです。よかったよかった。

学んだこと

  • ベイズ理論を用いたフィルタリングは結構面倒かと思ってましたが、意外と楽でした。まぁ面倒な部分はMeCabがやってくれたっていうのもあるんでしょうけど。
  • pythonのreduce()はとてもエレガントに扱えて便利。たった一行ですむっていうのは気分いいですね。
  • どうでもいいですけどフレッシュネスバーガーのスパムバーガーはおいしいのでお勧めです。

うーんおいしそう!明日食べに行ってきます。

*1 どれくらい忘れていたかというと、"MeCab python"で検索したら自分の記事がトップに出て、「そんなん書いたっけ?」と驚いたくらい。自分の欠点は熱しやすく冷めやすいところだと思う
*2 出会い系とかに登録していればたくさん来ていたかもしれない
*3 本来は4000件くらい作るらしいんですが、さすがに勘弁してください

追記

Robinson方式も扱えるようにレポジトリを更新しました。Robinson方式は未知の単語が出てきた際に、Graham方式よりも適切に判定が行えるアルゴリズムです。

$ python test.py
=====Graham方式による判定=====
センテンス1がスパムである確率 : 1.0
センテンス2がスパムである確率 : 6.07068879794e-100
=====Robinson方式による判定=====
センテンス1がスパムである確率 : 0.999957539928
センテンス2がスパムである確率 : 0.0234555720353
=====辞書に関する情報=====
スパム文章のサンプル数 : 20
非スパム文章のサンプル数 : 20
スパム辞書が学習した単語数 : 1040
非スパム辞書が学習した単語数 : 854

Robinson方式に対応させるために逆カイ二乗関数が必要になったので、今回はchi2pモジュール(http://garyrob.blogs.com/chi2p.py)を利用することにしました。chi2pモジュールはMIT Licenseの元で配布されています。

2009/03/05

PyBrain - a modular Machine Learning Library for Python


今、python界でPyBrainが熱い!…わけじゃないですけど、個人的にけっこう注目しているライブラリ。機械学習ライブラリにおける、期待の新人が出てきなような気持ちです。

0.PyBrainとは?

PyBrainっていうのはPythonによって動く、モジュール式の機械学習ライブラリです。python界ではいままでにもニューラルネットワークとかSVMなどを扱うライブラリが存在していましたが、PyBrainではそれらをより包括的に扱う、一種の環境としての機械学習ライブラリを目指しているようです。
PyBrainが優れているのはその思想もさることながら、扱っているアルゴリズムの多さにもあります。例えばFeaturesの欄を見てみると、
  • Backprop
  • Rprop
  • Policy Gradients
  • Support Vector Machines
  • Evolution Strategies
  • CMA-ES
  • Competitive Coevolution
  • Natural ES
  • Natural Actor-Critic
  • Reward-weighted regression
  • Evolino
  • Fitness Expectation Maximization
  • SPLA
  • PCA/pPCA
  • LSH for Hamming and Euclidean Spaces
と、比較的新しいアルゴリズムについても対応する予定のよう。(リファレンスをざっと見ると、まだ対応していないアルゴリズムはけっこう多いようです)とりあえずアルゴリズムについて確認したり、テストとして組んでみる分には申し分ない仕様です。

1.ダウンロード〜インストールまで

まず、PyBrainでは以下のツールが必要となるので、足りない場合は適宜apt-getやらeasy_installやらでインストールを行います。
  • g++
  • scipy
  • matplotlib
ちなみにホームページには書かれていませんが、自分はビルドした際にcblasが足りないと怒られてしまったので、ついでにそれもインストールします。
次に、PyBrainのサイトにアクセスして、Stable versionのソースコードを入手。圧縮ファイルを展開します。

~$ cd PyBrain-0.2
~/PyBrain-0.2$ ls
LICENSE arac docs examples pybrain setup.py

setuptoolsの慣例通りにビルド。

~/PyBrain-0.2$ sudo python setup.py build
arac/src/c/layers/common.c: In function ‘void make_layer(Layer*, int, int)’:
arac/src/c/layers/common.c:12: error: ‘malloc’ was not declared in this scope
arac/src/c/layers/common.c: In function ‘Layer* make_layer(int, int)’:
arac/src/c/layers/common.c:40: error: ‘malloc’ was not declared in this scope
Traceback (most recent call last):
File "setup.py", line 87, in <module>
compileArac()
File "setup.py", line 74, in compileArac
extra_postargs=['-O3', '-g0', '-DNDEBUG'])
File "/usr/lib/python2.6/distutils/ccompiler.py", line 697, in compile
self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
File "/usr/lib/python2.6/distutils/unixccompiler.py", line 180, in _compile
raise CompileError, msg
distutils.errors.CompileError: command 'g++' failed with exit status 1

怒られてしまいました。どうやらmallocが定義されていないみたいなので、暫定的な処置として"arac/src/c/common.h"に以下の文を付け加えます。

#include <stdlib.h>

再度ビルドを行うと、ビルドがうまくいったみたいです。

~/PyBrain-0.2$ sudo python setup.py install

でインストール。無事インストールが完了しました。
正常に動くかどうか確かめてみます。

~/PyBrain-0.2$ python
Python 2.6.1+ (r261:67515, Feb 24 2009, 20:00:00)
[GCC 4.3.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import pybrain
>>>

何事もなくインポートすることができました。

2.試しに使ってみる〜ニューラルネットワークの構築から学習

Pybrainのチュートリアル通りに、単純なFeed Forward型のニューラルネットワークを組んでみます。Pybrainではレイヤを個別に定義し、それらを繋げてネットワークを構成することもできるみたいですが、もっと簡単にネットワークを構成できるショートカットがあるみたいですので、それを利用します。

from pybrain.tools.shortcuts import buildNetwork
net = buildNetwork(2, 3, 1)


これで2つのインプット、3つの隠れユニット、1つのアウトプットから成るニューラルネットワークが出来上がりました。とても簡単ですね。
入力したインプットからアウトプットを出すには、activate()メソッドを用います。

net.activate([2,1])

buildNetworkの時点で各ニューロンはランダムな数値で初期化されていますので、既にactivate()ができる環境になっています。

ちなみに通常、各ニューロンに用いられる活性化関数はシグモイド関数が用いられますが、違う関数を用いたい場合は適宜指定することができるみたいです。まぁ好みの問題ですね。

from pybrain.structure import SoftmaxLayer
from pybrain.structure import TanhLayer
net = buildNetwork(2, 3, 1, hiddenclass=TanhLayer, outclass=SoftmaxLayer)


次に、教師あり学習を行うためにデータセットを用意しておきます。

from pybrain.datasets import SupervisedDataSet
ds = SupervisedDataSet(2, 1)
ds.addSample((0, 0), (0,))
ds.addSample((0, 1), (1,))
ds.addSample((1, 0), (1,))
ds.addSample((1, 1), (0,))


こんな感じでインプットとアウトプットのサンプルを増やしていきます。
そして、構成されたネットワークに誤差逆伝播法で学習を行わせるため、trainerを作成します。

from pybrain.supervised.trainers import BackpropTrainer
trainer = BackpropTrainer(net, ds)
trainer.train()


train()メソッドは返り値として誤り率を返しますが、もしある一定の割合に収束するまで学習を行わせたい場合、trainUntilConvergence()メソッドが有効です。そのまんまですね。

trainer.trainUntilConvergence(validationProportion=0.25)

だいたいニューラルネットワークの構築方法としてはそんな感じですが、他にもSVMを選択することもできるようです。詳しいことはPyBrainのドキュメントを参照してください。

3.使ってみて

とりあえずざっと使ってみた感想ですが、非常に組みやすく、かゆいところに手が届くような設計がなされているライブラリだと思います。まだバージョンは0.2ですし未完成な感じは否めないですけれど、将来が非常に楽しみなライブラリです。(しかも名前がカッコいい!!)

2009/02/12

Ambient Occlusion with Python


前々からやってみたかった、syoyoさんのAO benchでのコードをPython用に移植させてみました。上の画像はその出力結果です。ちなみに、出力を楽にするためにPython Imaging Library(PIL)を利用しています。

なんで前々からやってみたかったかというと、3Dでのレンダリングアルゴリズムの学習がしたかったから。ray tracing法やらAmbient Occlusionなどの手法は頭の中では分かってましたが、実際にどうやって動かしているのかは全然分かんなかったですから…今回の移植を通して、レンダラの実際の動きが少しだけですが掴めたような気がします。

それで気になる速度はというと、

  • サブサンプリング数:2
  • AOのサンプル数:8
  • 画像の長さと高さ:256px

という標準的な設定でレンダリングしたところ、
$ python sample.py
time:463.692063 sec
という結果が出ました(AMD Athlon 64 x2 3800+)。

まぁPythonですから大体予想してた通りの結果となりましたが、それにしてもものっそい遅いです。C言語での2.6secという結果を比べてみると、環境が違いますからそのまま評価することはできないんですがだいたい180倍くらい時間かかってます。もう少しなんとかできなかったかなー…

2008/12/08

Picasaのウェブアルバムに写真を自動的に載せるpythonスクリプトを組んでみた。

これはホント前々から作っておきたかったスクリプトの一つです。

Picasaのウェブアルバムは無料ながらにして容量が1GBととても大容量なので、クローズドに写真をアップロードする場所としてとても重宝してきました。だけどどうしても納得できないことが一つだけあって、それが写真のアップロード方法。


なんで一度に5枚しかアップロードできんのじゃい

普通アルバムに画像をアップロードしたいと思ったら2〜30枚くらいは軽くこえるだろうーどうして5枚だけ、しかもいちいち選択し直さなきゃならないのかとっても謎。

いや、「そんなことしないでPicasa使えば楽にアップロードできるのに」って思う人がいるかもしれません。そりゃそうなんですけど、個人的にあのアプリケーションは好きになれないのです。べつにキャッシュしなくてもいい画像ファイルもキャッシュしてしまうし。そんくらい自分でちゃんと整理します。

というわけで不自由なら自分でなんとかするの精神で、Picasaのウェブアルバムに指定したフォルダーの写真をすべて自動的にアップロードするpythonスクリプトを作りましたので公開します。
まず、使用するには以下の環境が必要です。

ソースみれば大体何やってるかわかると思いますが、画像が大量に入っているフォルダを用意しておいて、
$ python picasa.py new_album ./photos/
とか指定してやれば"new_album"っていうアルバムを新しく作って、./photos/ フォルダの画像すべてを自動的にアップロードしてくれます。すごく楽ちん。

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

import sys, os
import gdata.photos.service

USER_EMAIL = "ログインに用いるeメールアドレス"
USER_PASSWORD = "パスワード"

def main():
# ログイン開始
gd_client = gdata.photos.service.PhotosService()
gd_client.email = USER_EMAIL
gd_client.password = USER_PASSWORD
gd_client.source = "Picasa-AutoUploadApp"
print u"Picasaへのログインを開始します。"
gd_client.ProgrammaticLogin()

# アルバムを新しく追加
target_album = gd_client.InsertAlbum(title=sys.argv[1], summary='')
print u"新しいアルバム \"%s\" を追加しました。" % sys.argv[1]

# 指定されたディレクトリから写真の読み込み
filepath = []
for root, dirs, files in os.walk(sys.argv[2]):
for f in files:
if ( f.endswith(".JPG") or f.endswith(".jpg") ):
filepath.append( (os.path.join(root, f), f) )

# 写真の追加
print u"写真の読み込み完了。アップロードを開始します。"

album_url = "/data/feed/api/user/default/albumid/%s" % target_album.gphoto_id.text
for ( i , (photo_path, photo_name) ) in enumerate(filepath):
print u"(%i/%i) \"%s\"をアップロードしています..." % (i+1, len(filepath), photo_name),
gd_client.InsertPhotoSimple(album_url, photo_name, '', photo_path, content_type="image/jpeg")
print u"\t完了"

print u"すべての写真のアップロードが正常に行われました。"

if __name__ == '__main__':
sys.exit(main())

これでやっとひたすらボタンとにらめっこする重荷から開放される…

追記

これを書いている途中で気になったんですが、pythonってprint文で自動的に改行させないようにしたい場合、
print "hogehoge",
ってコンマを最後に付けるんですね。これなんとかなんないでしょうか?python 3.0ではさすがに修正されているとは思いますが、すんごくsyntax errorっぽい…

2008/09/20

エッシャーっぽい絵を生成する「エッシャーくん」を作ってみた。

エッシャーくん

エッシャーっていう画家は知っていますか?分かんない人のために説明しますと、こんな感じのふしぎーな絵を書いている人です。名前は知らなくても一度は見たことがあるのではないでしょうか。

それでなんですが、適当な画像からなんかエッシャーっぽいふしぎな画像を生成するフィルタ「エッシャーくん(仮称)」をPython Imaging Library(PIL)で作ってみました。これを使えばどんな画像もエッシャーっぽい世界にご招待です。ソースは近々公開します。

追記(09/24)

ソースコードをアップロードしました。subversionで管理されてますので、
svn checkout http://svn.coderepos.org/share/lang/python/escher Somewhere
でチェックアウトしてください。


たとえば、こんな感じのイラストにエッシャーくんを適用させてみると…


こんな感じの絵が出来上がりました。これはどこまでズームさせても終わりがない、無限に続いている絵になっています。


この、なんかごっついリングもエッシャーくんの手によれば…
(Photo by the justified sinner, "Roofing Nail Ring 1")


どこまでも続く絵に早変わり!なんかうにょんうにょんしています。


この画像はどーなるんだろ…


なんかコーヒーっぽい…

どーやって作っているのか

エッシャーくんは"A logarithmic image transformation"っていう論文を元に作られているので、気になる人はこれを見れば大体のことが分かると思いますが、一応大まかな解説をしてみると、複素平面をlog(z)で変形して、絵が循環するように座標をスケール、回転させてからexp(z)で繋げているようです…

あーやっぱ説明が下手くそなので、見てもらったほうが早いです。

余談

本当はprocessingで作る予定だったんですけど、あまりにもprocessingの数学ライブラリがアレだったんで途中からpythonで作り直しました。(対数関数すらないのにはビックリした。そりゃああんま使わないかもしれないけどさー)

それにしてもpythonの開発速度の早さは半端ない。まさか100行未満で出来るとは思わんかったです。

2008/08/13

どうやらDjangoの登録方法が変更になったみたい

どうやらDjangoはバージョン0.97から管理画面にModelを登録する方法が変更されているみたいで。今までは対象のModelクラスにAdminクラスメソッドを作る方法を採用していたんですが、Modelクラスに管理画面の設定が混在せざるを得ない仕様はやっぱり完璧主義者には馴染めないものだったんですね。

で、具体的にどうするかっていうと対象のアプリケーションディレクトリに新しくadmin.pyファイルを作成して、以下のように記述します。

#!/usr/bin/env python
#-*- coding:utf-8 -*-
from blogs.models import Entry,Blog,Comment
from django.contrib import admin

##
## The Configuration for Admin.
##
class EntryInline(admin.StackedInline):
model = Entry
extra = 1

class BlogAdmin(admin.ModelAdmin):
inlines = [EntryInline]
list_display = ('title','pub_date','user')
list_filter = ['pub_date']

# registration
admin.site.register(Blog, BlogAdmin)

admin.site.register()の部分で管理画面の登録を行ない、具体的な設定は専用のクラスを新しく作成することで解決します。
しかしこのままではadmin.pyが呼び出されないので、プロジェクトのurls.pyの先頭に以下の2行を加えます。

from django.contrib import admin
admin.autodiscover()

これでadmin.pyが呼び出され、無事管理画面の登録が行なわれるようになりました。

Djangoのオンラインドキュメントではmodels.pyに直接Adminの記述を行なっているんですが、これだとviews.pyとかがmodels.pyをインポートしたときに再びadmin.site.register()が実行されてしまいますので、二重登録のエラーが生じてしまいます。モデルと管理画面の設定を分離するという意味でも、こちらのほうがより便利じゃないかなーと思います。

2008/07/14

Black-Scholes Equation in Python 改良版

前回のコードをある程度改良しましたので、載せておきたいと思います。主な改良点は以下。

  • ヨーロピアン(コールorプット)オプションのデルタ、ガンマ、ベガ、シータ、ローを求める関数を追加しました。
  • ヨーロピアン(コールorプット)オプションのインプライド・ボラティリティをニュートン法より求める関数を追加しました。
  • その他若干の名称変更や、引数に関する説明を加えました。

数式は"Black-Scholes (1973) Option Pricing Formula"を参考にしました。

インプライド・ボラティリティに関しては適切な引数を加えないと、負の値や無限大に発散することがあるので注意が必要です。とはいえ実際の市場価格と理論価格がそこまで乖離することもないでしょうから、大抵正常な値を返すはずです。

商用、非商用に関わらず自由に使用、改変していいですし連絡する必要もないですが、無保証です。使う人はあんまいないと思いますが、一応建前として。

#!/usr/bin/env python
#-*- coding:utf-8 -*-
"""
Functions for calculating the Black-Scholes equation.
"""
import math

A1 = 0.319381530
A2 = -0.35653782
A3 = 1.781477937
A4 = -1.821255978
A5 = 1.330274429
RSQRT2PI = 0.3989422804
GAMMA = 0.2316419

def BS_Call(S,K,T,R,V):
"""Calculate the Black-Scholes formula for a call option.

Arguments:
S:the price of the underlying stock
K:the strike price
R:the continuously compounded risk free interest rate
T:the time in years until the expiration of the option
V:the volatility for the underlying stock
"""
d1,d2 = _getD1D2(S,K,T,R,V)
ret = S*_normalDist(d1)-K*math.exp(-R*T)*_normalDist(d2)
return ret

def BS_Put(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
ret = -S*_normalDist(-d1)+K*math.exp(-R*T)*_normalDist(-d2)
return ret

def BS_Call_IV(S,K,T,R,Pr,HV,maxNum=100,i=1):
"""Calculate the Implied Volatility using Newton's method.

Arguments:
Pr:the real price of call option
HV:the Historical Volatility
maxNum:the number of repeating calculation
i:index(Do not specify this parameter)
"""
if i > maxNum:
return HV
else:
newV = HV - (BS_Call(S,K,T,R,HV)-Pr)/BS_Vega(S,K,T,R,HV)
return BS_Call_IV(S,K,T,R,Pr,newV,maxNum,i+1)

def BS_Put_IV(S,K,T,R,Pr,HV,maxNum=10,i=1):
if i > maxNum:
return HV
else:
newV = HV - (BS_Put(S,K,T,R,HV)-Pr)/BS_Vega(S,K,T,R,HV)
return BS_Put_IV(S,K,T,R,Pr,newV,maxNum,i+1)

def BS_Call_Delta(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return _normalDist(d1)

def BS_Put_Delta(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return _normalDist(d1) - 1.0

def BS_Gamma(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return _normalDist(d1)/(S*V*math.sqrt(T))

def BS_Vega(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return S*_normalDist(d1)*math.sqrt(T)

def BS_Call_Theta(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return -S*_normalDist(d1)*V/(2*math.sqrt(T))-R*K*math.exp(-R*T)*_normalDist(d2)

def BS_Put_Theta(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return -S*_normalDist(d1)*V/(2*math.sqrt(T))+R*K*math.exp(-R*T)*_normalDist(-d2)

def BS_Call_Rho(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return K*T*math.exp(-R*T)*_normalDist(d2)

def BS_Put_Rho(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return -K*T*math.exp(-R*T)*_normalDist(-d2)

def _getD1D2(S,K,T,R,V):
vSqrtT = V*math.sqrt(T)
d1 = (math.log(float(S)/K)+(R+(0.5*V*V))*T)/(vSqrtT)
d2 = d1 - vSqrtT
return (d1,d2)

def _normalDist(x):
if x >= 0:
k = 1.0/(1.0 + GAMMA*x)
cnd = RSQRT2PI*math.exp(-0.5*x*x)*(k *(A1+ k*(A2+ k*(A3 + k*(A4 + k*A5)))))
ret = 1.0 - cnd
return ret
else:
return 1.0 - _normalDist(-x)

て、そんなことやってる場合じゃない。早くテスト対策行なわないと。

2008/06/16

Pythonのソートについて: sort()とsorted()


Photo by Fort Photo, "Fountains of Light"
自分のためにメモ。pythonでソートを行おうとする場合、リストに組み込まれている関数list.sort()と、単独で使える関数sorted()の2つがある。基本的に取る引数は同じなので機能自体はあまり変わりないのだが、以下の2点が異なる。
  • sort()関数は返り値を出さずに、リストの本体自体を変更してしまう。これはreverse()関数と同じ。
  • sort()関数はソートの対象がリストであるのに対し、sorted()関数は対象がイテレータとなる。
試しにpython上で確認してみる。
>>> x = [1,6,3,8,4]
>>> y = [1,6,3,8,4]
>>> y.sort()
>>> y
[1, 3, 4, 6, 8]
>>> sorted(x)
[1, 3, 4, 6, 8]
きちんと両方ともソートが行われているようだ。それでは次に引数cmp, key, reverseについて調べてみる。
"cmp"キーワードには、二つのオブジェクトxとyがあったときに、どちらが大きいのかを示す関数、あるいはλ式を指定する。普通はそのままで十分であるが、比較用のオブジェクトがちょっと特殊な場合に役立つ。
もしcmpキーワードを指定しなければならないような場合、通常ならばλ式の中にcmp(x,y)関数を導入する。cmp(x,y)関数はxとyを比較して、指定すべき適切な値を返してくれる関数。これを導入すれば、
あれ?cmpキーワードってx < yの場合にはx - yを指定すればいいんだっけ?それともy - xだっけ?』とかいうような余計な心配をしなくて済む。
"key"キーワードには、比較対象とするオブジェクトはどのオブジェクトであるのかを示す関数、あるいはλ式を指定する。何も指定しなかった場合、配列やタプルであれば一番目(インデックスが0)のオブジェクトが比較対象であると判別されるようだ。
"reverse"キーワードはソート結果が昇順であるか降順であるかを指定する。Trueを指定すると降順になる。sort()してからreverse()などのように2つメソッドをくっつける必要はない。
それでは今までの流れを踏まえて試しに簡単なプログラムを組んでみる。辞書には{"名前":"点数"}などのように名前と点数のペアが保存されている。これを点数が高い順に並び替えるというプログラムは以下のように記述できる。
>>> dic = {"Hanako":"50", "Taroh":"80", "Ichiroh":"30", "Yoshio":"60"}
>>> sorted(dic.items(), cmp=lambda x, y:cmp(int(x), int(y)), key=lambda x:x[1], reverse=True)
[('Taroh', '80'), ('Yoshio', '60'), ('Hanako', '50'), ('Ichiroh', '30')]
ちなみに上の式は以下のようにも記述できる(というかそっちのほうが普通)。
>>> sorted(dic.items(), key=lambda x:int(x[1]), reverse=True)
[('Taroh', '80'), ('Yoshio', '60'), ('Hanako', '50'), ('Ichiroh', '30')]
6/30 追記
最初のソート部分が間違っていたので修正しました。以前の書き方だとxとyは同一の配列オブジェクトを参照していることになりますね。すいません。

2008/06/01

Ubuntu 8.04 (Hardy Heron)にRpyをインストールする。

pythonの上からRを扱いたかったので、RPyをUbuntu 8.04(Hardy Heron)にインストールしようと試みたんですが、どうやらこのまま
$ sudo apt-get install python-rpy
とかでインストールすると、インポートの段階でRPyモジュールが足りないなどのエラーがでるようです。

えーそれじゃあソースからコンパイルするしかないのかなぁ面倒くさいなぁとか思っていろいろ調べてみたら、どうやら8.04だけではなく、Ubuntu全般で起こる問題なようでいろんな人が書き込みを行っていました→(rpy import error at Gusty)。自分の他に躓いた人もいるかもしれないのでメモしておきます。

いろんな方法がありますが、いまのとこ一番楽な解決方法は以下のとおり。

  1. "http://cran.r-project.org/bin/linux/ubuntu/"にアクセス。自分が使っているヴァージョンを選択(自分の場合は"hardy")
  2. "r-base-core_2.7.0-1hardy0_i386.deb"と"python-rpy_1.0.3-1hardy0_i386.deb"の両方をダウンロードしてインストールする(ただし自分の場合AMD64を使っていたので、amd64をインストールしました)。警告がでるかもしれないけど気にしない。
  3. 無事インストールが完了したら、
    $ python
    >>> import rpy
    できちんと動作するか確認する。

本来こんなどーでもいいことで時間が取られるのはあんま好きじゃないんですが、いやぁーとりあえず解決できてよかったよかった。

2008/05/27

Pythonでiniファイルとかの設定ファイルを読み込んでみる

Python上でiniファイルなどの設定ファイルを読み込むための方法。Pythonは"Battery Included"の発想があるため、特に何も用意しなくても設定ファイルを読み込むライブラリ"ConfigParser"が標準でついてくる。このConfigParserがとっても便利で感動したのでおもわず投稿。

ConfigParserを用いて設定ファイルを読み込むには、ただ以下の数行を付け加えるだけでよい。



import ConfigParser

config = ConfigParser.RawConfigParser()
config.read('settings.conf')
directory = config.items('Directory') #Directoryセクションの読み込み
print directory # [('item1','value1'),('item2','value2'),...]
print config.get('Directory','item1') # 'value1'

実質設定ファイルの読み込みが3行ですんでる!すごい!

これをC#で行うとなると結構大変で、StringBuilderなどを使って構文解析をする必要がある。
実際はもうiniファイルではなく、XMLなどを用いて設定ファイルを読み込むのが主流になっているはずだからフェアな比較とは言えないけど、ちょっとしたアプリケーションを作るときにXMLは少し気合が入りすぎというか、そこまでしなくても・・・という感じはある。


なのでこうやってたった3行、しかもリファレンスを見ただけで動作がきっちり理解できるpythonはやっぱすごいなと思ったのでした。

2008/05/08

python + BeautifulSoupでHTML解析を行ってみた。

諸事情により株価について整理されたデータを作ることになったんで、pythonとBeautifulSoupを使ってHTMLスクレイピングを行うことにしました。

最初は「えーpythonってインデントで識別するんでしょー慣れなさそー気持ちわるそー」って思ってたんですが、実際やってみると簡単で数時間後にはきちんと動作するスクリプトが組めたのでびっくり。すいません侮ってました。これは便利!

いやぁーそれにしてもすごいっすよBeautifulSoup。なにしろ名前が綺麗っすよね。混沌としたHTMLのSoupからBeautifulな部分を取り出すぜ!っていう思いが伝わってきて凄くいいです。そしてなによりHTML解析が便利すぎる。他の言語はどーやって行っているのかよく分かりませんが、トップレベルで使いやすいことは確かです。

これがHTMLスクレイピングの主要部分↓

soup = BeautifulSoup(rawdata)      
target = soup("tr",{'bgcolor':'#ffffff'},{'align':'right'})[0]('small')

こんな気持ち悪い文法で動くのはBeautifulSoupだけ!:)

こんな文初めて書きました。いやぁーすごいっす。