Showing posts with label python. Show all posts
Showing posts with label python. Show all posts

2011-05-01

pypyで何も考えずにpythonを高速化する

Python公式のCPythonよりも高速という噂のpypyPython2.7.1相当の機能を実装したというので、その実行速度を測定しました。



■インストール(OSX)

pypyのインストールは非常に簡単で、本家からbzipをダウンロードして解凍したらできあがりです。


$ tar xjvf pypy-1.5-osx64.tar.bz2
$ cd pypy-c-jit-43780-b590cf6de419-osx64
$ ls
LICENSE README bin include lib-python lib_pypy site-packages

binの中にpypyというコマンドがありますので、それを実行すると起動します。


$ bin/pypy

Python 2.7.1 (b590cf6de419, Apr 30 2011, 03:30:00)
[PyPy 1.5.0-alpha0 with GCC 4.0.1] on darwin
Type "help", "copyright", "credits" or "license" for more information.
And now for something completely different: ``how to construct the blackhole
interpreter: we reuse the tracing one, add lots of ifs and pray''
>>>>

>が3つではなく4つ出てくるところが、公式CPythonよりも強そうですよね!



■実測

実際にPython標準のパフォーマンス測定ツールtimeitで、実行速度を比較してみました。

比較のため、OSX上のCPython2.6.1, 2.7,1とpypyを使いました。



まずは単純に配列を作る処理。


------- CPython 2.6.1
>>> import timeit
>>> t = timeit.Timer(stmt="[x*x for x in xrange(1000)]")
>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100000)/100000)
184.11 usec/pass

------- CPython 2.7.1
>>> import timeit
>>> t = timeit.Timer(stmt="[x*x for x in xrange(1000)]")
>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100000)/100000)
153.30 usec/pass

------- pypy 1.5
>>>> import timeit
>>>> t = timeit.Timer(stmt="[x*x for x in xrange(1000)]")
>>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100000)/100000)
120.95 usec/pass


CPython2.7は2.6よりもだいぶ速いですが、それをさらにpypyが上回る速度です。




次に、集合(set)の中の検索


------- CPython 2.6.1
>>> t = timeit.Timer(setup="a=set([x*x for x in xrange(100)])", stmt="[(x in a) for x in xrange(100)]")
>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100000)/100000)
23.41 usec/pass

------- CPython 2.7.1
>>> t = timeit.Timer(setup="a=set([x*x for x in xrange(100)])", stmt="[(x in a) for x in xrange(100)]")
>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100000)/100000)
19.59 usec/pass

------- pypy1.5
>>>> t = timeit.Timer(setup="a=set([x*x for x in xrange(100)])", stmt="[(x in a) for x in xrange(100)]")
>>>> print "%.2f usec/pass" % (1000000 * t.timeit(number=100000)/100000)
11.97 usec/pass

これも結果は同様に圧倒的にpypyが高速。(2.7の2倍に迫る速度!)



こんなに簡単に高速化できるなら、どこかで試しに使ってみたくなりますよね。

2011-02-08

Ajaxを補足するためのPythonによる手軽なマルチスレッド処理

Webアプリの開発で、重い処理はAjaxで別途処理して取得するという手段が取れない場合があります。

「画面の一部で矢印がぐるぐる回って待たせるのではなく、本当に一瞬で全部表示したい」

などの場合です。

そんなときにお手軽に使えるテクニックとして、サーバー処理のマルチスレッドによる並列化があります。
Pythonのマルチスレッドでは、高度なCPUの並列処理はできませんが(グローバルインタプリタロックと言います)、少なくともI/O待ちの間別の処理を流す程度の事は可能です。[参考]
その特徴を生かして、Webアプリのレスポンスを速めることができます。


これを使う場面ですが、たとえば、以下のようなものがあります。



この画面はうちの会社で作っている特許分析WEBアプリのBiz Cruncherというもので、特許の詳細な内容を表示しているところです。


特許の閲覧というのは企業知財部では古来から「手めくり」と呼ばれ、いかに速く大量の特許の内容を把握するかに心血を注ぎます。
そのため、図などは一瞬で見れるべきだし、電子的な画面であろうと、一瞬で次々表示されなければストレスを感じるものです。


そういう場合は付属情報はAjax化で別リクエストにするにしても、必須の要素だけは泣く泣く1リクエストに押し込んで、極限まで速度を向上させる必要があります。

たとえば、上の特許詳細画面では時間のかかる処理として、テキストマイニング、画像の生成、詳細テキスト(長文)の取得、引用情報・経過情報の取得 etc.をそれこそ一瞬で行う必要がありました。

改良前はこの表示に2.8秒かかっていたのですが、それでも不満ということで、冒頭のマルチスレッド化を行って0.7秒にまで縮めることができました。

具体的には以下の関数 call_async を使います。

import threading
def call_async(target, args=None, kwargs=None):
u"""
関数非同期呼び出しツール。
関数を別スレッドで呼び出して、その後結果を取得できます。

#使用例
r = call_async(f,args=(5,))
#ここで別の処理を片付けます。
r["thread"].join() #スレッドが終わるのを待ちます。
if "result" in r:
print "RESULT", r["result"]
else:
print "ERROR", str(r["exception"])
"""
res = {}
def f(*inargs, **inkwargs):
try:
inargs[-1]["result"] = target(*(inargs[:-1]), **inkwargs)
except Exception, e:
inargs[-1]["exception"] = e
res["thread"] = threading.Thread(target=f, args=(args+(res,)), kwargs=kwargs)
res["thread"].start()
return res



この関数を使うと、関数呼び出しを別スレッドにしながら、その返り値を取得することもできます。
使用方法は上記関数のコメント内にありますが、以下のように、先に遅い処理を実行開始しておき、平行してその他の処理を済ませたあとで、ゆっくりと結果を取得します。

#使用例 slow1~3という遅い処理の関数がある場合

#まず冒頭で遅い処理を別スレッドで実行開始しておきます。
r1 = call_async(slow1,args=(5,)) #遅い処理1を開始!
r2 = call_async(slow2,args=("ABC","DEF")) #遅い処理2を開始!
r3 = call_async(slow3,args=([1,2,3,4])) #遅い処理3を開始!

#ここでその他の処理を片付けます。

r1["thread"].join() #スレッドが終わるのを待ちます。
r2["thread"].join() #スレッドが終わるのを待ちます。
r3["thread"].join() #スレッドが終わるのを待ちます。

if "result" in r1: #slow1の結果を使う処理
print "RESULT 1", r1["result"]
else:
print "ERROR 1", str(r1["exception"])

if "result" in r2: #slow2の結果を使う処理
print "RESULT 2", r2["result"]
else:
print "ERROR 2", str(r2["exception"])

if "result" in r3: #slow3の結果を使う処理
print "RESULT 3", r3["result"]
else:
print "ERROR 3", str(r3["exception"])


こうすることで、I/O待ちの発生する処理をいくつも並列で実行でき、実行時間の効果的な短縮を実現できます。

こうした手段としてPython標準のmultiprocessingが検討にあがりますが、あれはより高度な分散処理用ですので、関数の結果を取得しようとすると、とたんに大仰なものになります。

また、もちろんですが、この秒数の短縮で効果的なのはJavascriptの最適化もあるわけです。
それをした後に、さらにどうしてもサーバー側も高速化しないといけない場合は、このお手軽なマルチスレッド化も考慮するとよいと思います。

2010-12-04

PythonのUnicode->UTF-8変換はSJISやEUCへの変換よりもちょっと速い

想像していた通りですが、UnicodeをUTF-8に変換するのは、SJISやEUCに変換するよりもちょっと速いことがわかった。




>>> import timeit
>>> def unicode_henkan(c):
... t = timeit.Timer("a.encode('%s')" % c, "a=u'\u3042'*10000")
... print t.timeit(number=1000)
...
>>> unicode_henkan("utf-8")
0.075012922287
>>> unicode_henkan("sjis")
0.19597196579
>>> unicode_henkan("eucjp")
0.146716833115
>>> unicode_henkan("cp932")
0.21554517746
>>>


ちなみに、timeitは与えられた計測対象のプログラムを内部でpythonのプログラムとして解釈する際に、文字コードのことをあまり考えてくれないようなので、

t = timeit.Timer("a.encode('%s')" % c, "a=u'あ'*10000")

のように書くとUnicodeErrorになります。

2010-09-14

関数呼び出しをタプル形式にして遅延させる関数

関数呼び出しをタプルの情報にして保存するしくみを作った。
これを使えばWEBアプリでリクエストをまたいで関数呼び出しを遅延させることができるようになります。
関数呼び出しをタプル形式にして、memcacheなどに保存しておき、次に来たときに本当に呼び出したり。

# coding: utf-8
import inspect
import imp

def freeze_func(func, *args, **kwargs):
u"""
関数呼び出しを保存可能なタプル形式に変換する。
"""
m = inspect.getmodule(func)
return m.__name__, func.__name__, args, kwargs

def __importer(m, path):
a = imp.find_module(m[0], path)
b = imp.load_module(m[0], *a)
if len(m) > 1:
return __importer(m[1:], b.__path__)
else:
return b

def _importer(modname):
modname = modname.split(".")
return __importer(modname, None)



def unfreeze_func(info):
u"""
タプル形式で保存された関数呼び出しを関数に復元して実行する。
"""
b = _importer(info[0])
f = getattr(b, info[1])
args = info[2]
kwargs = info[3]
return f(*args, **kwargs)

if __name__ == "__main__":
import re
info = freeze_func(re.sub, r"[A-Za-z]", "*", "5-3-1 Asakusabashi")
print info # -> ('re', 'sub', ('[A-Za-z]', '*', '5-3-1 Asakusabashi'), {})
print unfreeze_func(info) # -> "5-3-1 ************"

2008-05-23

はてブのホットエントリをなんとなーく分類するPythonスクリプトサンプル

昨日作った、Ward法のPythonバインディングを使って、はてなブックマークのホットエントリをなんとなーく分類するPythonスクリプトのサンプルを書いてみた。

サンプルプログラムは例によって、CodeReposに置きました。

試しに今の時点(5/23 深夜2:06)のホットエントリを分類した結果は以下のような感じです。
なんとなーく、分類されているような雰囲気。
もっと各エントリの情報を取得すれば精度は上がるんでしょうが、今の所はホットエントリのRSS一枚から取得できる情報(今の所タイトル中の名詞と、タグの単語)だけを使って分類しています。

----- 以下、その結果 -----

Group [ 50 ] --------------------
「トレードオフの概念は日本に無いのか」 三菱東京UFJ銀のシステム一本化報道に思う:NBonline(日経ビジネス オンライン)
「トレードオフの概念は日本に無いのか」 三菱東京UFJ銀のシステム一本化報道に思う:NBonline(日経ビジネス オンライン)
なぜか変換できない vs. なぜか変換できる:最近の「MS-IME」は目に余る――よろしい、ならば「ATOK」だ (1/4) - ITmedia +D PC USER
今日から始めるデジカメ撮影術:第97回 一眼レフとボケの関係 (1/3) - ITmedia D LifeStyle
NHK「クローズアップ現代」児ポ法単純所持処罰に絞って特集 - フィギュア萌え族(仮)犯行説問題ブログ版・サブカル叩き報道を追う
<ネットカフェ難民>「しんどい」と訴える妊婦まで 100人の実態調査(毎日新聞) - Yahoo!ニュース

Group [ other ] --------------------
Steve Jobs氏のようにプレゼンテーションをする方法 - builder by ZDNet Japan
調査結果に基づいて合成「人々が最も不愉快になる曲/最も聴きたい曲」:試聴可能 | WIRED VISION
グーグル、検索アルゴリズムを少しずつ明らかに:マーケティング - CNET Japan

Group [ 59 ] --------------------
愛し合うってすばらしい、セックスすると得られる10の効能 - GIGAZINE
はてなクラブ - はてな
十字軍はバカに勝てるか - モジモジ君の日記。みたいな。
「宅配寿司 銀のさら」の自虐CM | DIGITAL DJ
変な人があまりにも酷すぎる件について
高学歴のお嬢様をお持ちの方、教えて下さい - 発狂小町
目を開いているのに閉じていると感じる錯覚作成法 | 医学都市伝説

Group [ 60 ] --------------------
ケーキを売ればいいのに - 福耳コラム
落合博満ガンダム宣言
電撃速報!! Amazonを倉庫代わりにしていた転売厨終了のお知らせ
『らき☆すた』を見た非オタの友人が衝撃の一言 - GilCrowsのペネトレイト・トーク

Group [ 64 ] --------------------
世界のPhotoshopチュートリアル トップ50 - ITクオリティ
光源、光線、光跡を表現するPhotoshopのブラシ -2XGU.NET | コリス
歌舞伎町で遇ったギークなタクシー運転手の編み出した、群衆の英知を活かした馬券購入法 - 雑種路線でいこう

Group [ 68 ] --------------------
人間関係が苦手で社会に適合できない人向けのDVDソフト「ミテルだけ」が登場 - GIGAZINE
YouTubeから削除された動画の墓場「YouTomb」が登場、日本のアニメの削除件数は圧倒的 - GIGAZINE

Group [ 72 ] --------------------
AIR-users.jp - 日本の AIR ユーザのためのハブサイト
js-users.jp - 日本の JavaScript ユーザのためのハブサイト
perl-users.jp - 日本の Perl ユーザのためのハブサイト

Group [ 75 ] --------------------
痛いニュース(ノ∀`):「美少女ゲーム(エロゲ)・アニメは、人間性を破壊し少女殺害事件につながる。販売規制を」…民主党女性議員が請願
「美少女ゲーム・アニメをする人は心を破壊され、人間性を失っているので規制すべき」と主張するトンデモ請願が参議院に - GIGAZINE
痛いニュース(ノ∀`):「ジョジョの奇妙な冒険」、中東で「日本が侮辱」「テレビ局爆破しろ」と批判殺到→集英社、原作の第3部など出荷停止

Group [ 76 ] --------------------
人工知能の叛乱
Gmailの検索オプションまとめ | IDEA*IDEA
マックユーザがインストールしているアプリケーション - soundscape out
[okyuu.com] ソーシャルITメディア
透過をきれいに使ったウェブデザインいろいろ - DesignWalker
特集:Firefox 3とFirebugで始めるJavaScript開発|gihyo.jp … 技術評論社
PHP コード最適化 Best Practices 63+ - カタコト日記

Group [ 77 ] --------------------
「ジョジョの奇妙な冒険」が中東で非難 - 社会ニュース : nikkansports.com
「コーラン読み殺害指示」日本アニメ、中東で非難 - MSN産経ニュース
asahi.com:トヨタ、「カイゼン」に残業代 業務と認定、来月から - ビジネス

Group [ 78 ] --------------------
ヒトラー名言集:アルファルファモザイク
傷つきやすい人はどうすれば強くなる?:アルファルファモザイク
ウマーな牛スジカレーの作り方:アルファルファモザイク
ねとらぼ:米WIREDに「やる夫」の特大AA掲載 ひろゆき氏インタビューも - ITmedia News
1日500アクセス以上のブログは全ブログの●%:Garbagenews.com
@nifty:デイリーポータルZ:50年前の新聞広告

Group [ 79 ] --------------------
MOONGIFT: � ブラウザベースのSubversion管理ツール「USVN」:オープンソースを毎日紹介
CSSで背景画像をリサイズするテクニック:phpspot開発日誌
ナムコの成果主義、有能なゲーム開発者が報われない開発現場

2008-05-22

Ward法のPythonバインディングを作りました。

先日から作っていたクラスタリング(Ward法)のPython用バインディングを作りました。Pythonから一応つないで使えるようになりました。

ソースコードの取得とインストール。


$ svn export http://svn.coderepos.org/share/lang/cplusplus/misc/clustercpp
$ cd clustercpp/python/wardcluster
$ python setup.py build_ext --swig-cpp
$ sudo python setup.py install

サンプル実行。
引き続きディレクトリclustercpp/python/wardcluster内で以下を実行してください。

$ python sample.py
1 + 0 => 4 ( 0.244948974278 )
3 + 2 => 5 ( 0.412310562562 )
5 + 4 => 6 ( 4.80645253132 )
$

ちなみに、今上で動いたsample.pyは以下のようなスクリプトで、4本のベクトルをクラスタリングしています。
from wardcluster import Matrix, DoubleVec, Ward

m = Matrix()
m.append(DoubleVec([1.0, 2.0, 3.0]))
m.append(DoubleVec([1.1, 2.2, 2.9]))
m.append(DoubleVec([0.1, 1.2, 5.0]))
m.append(DoubleVec([0.3, 1.0, 5.3]))

w = Ward()

history = w.cluster(m, 4) #この4ってのは、ベクトルの本数

for h in history:
print h.index1, "+" , h.index2, "=>", h.newindex , "(", h.distance ,")"


メモ--------
* setuptoolsを使ったオシャレなモジュールにするのは断念。(swigの結果ファイルとの調整が難しかった)
* swigのtypemap?とか、そういうのはむずかしくて分からないので、めんどくさいインターフェース。
* あくまで実験的なものですよ。
* 相変わらずswigを理解できていない。
* 今回ベクトルの集合を表すために内部でstd::vector<std::vector<double> >を使うようにしたので、パフォーマンスが落ちています。(将来もとのstd::vector<double>*を使って渡す方式に戻したいものです。。。)

2008-05-16

Python -> C++への移植でどれぐらい速くなったか

先日、Pythonで実装したWard法をC++で再実装し、比較してみた。

まずは、前回のPython版の処理時間を再掲。

件数   処理時間(s)    そのうち初期行列作成時間
--------------------------------------------
196 0.677 0.334
392 2.158 1.356
588 4.636 3.073
1176 18.071 12.088
2352 79.09 49.507


これをC++に移植すると、以下の処理時間に短縮される。

件数 処理時間(s) そのうち初期行列作成時間
--------------------------------------------
196 0.216 0.04
392 0.842 0.17
588 1.627 0.39
1176 7.086 1.57
2352 33.432 6.39

今回は、まずは行列作成部分に、できるだけC++で可能なチューニングを施してある。
樹形図作成部分は、あまり最適化していない。
標準のstd::mapとかstd::vectorを使用しており、特別な数値計算用ライブラリなどは一切使っていない。
また、コンパイル時の様々なコンパイルオプションでの最適化もまだ行っていない。

C++(C)のいい所は、メモリの使い方などをダイレクトに想像しながらチューニングが行える所だろう。
今回も、3重ループの効率化を行うだけで、ここまで処理時間が短くなった。

しかし、作っていて気づいたのは、最適化を意識しないでC++を書くと、Pythonとたいして変わらない処理時間だったことだ。
その証拠に、上の表での処理時間の短縮分のほとんどは、今回意識して最適化を図った行列作成部分のものである。

次は、実際の樹形図作成部分がネックになっているので、この部分をより速くしていきたいと思う。

今回のソースはcodereposに置いておきました。

2008-05-15

Ward法さらにキャッシュを効かせて高速に

さらに昨日の続き。

行列の各行の最小距離成分をキャッシュするようにしてみた。
この改造によって、昨日O(N**2)でネックとなっていた、距離表からの最小値の探索がほぼO(N)になったことになる。
さらにパフォーマンスは向上し、以下のようになった。

件数   処理時間(s)    そのうち初期行列作成時間
--------------------------------------------
196 0.677 0.334
392 2.158 1.356
588 4.636 3.073
1176 18.071 12.088
2352 79.09 49.507

昨日、おとといに比べると、だいぶ速くなった。例えば昨日は588件の処理に19秒かかっていたのが、今日は4.6秒!
しかし、計算は速くなったものの、今度は実は最初の行列作成に大半の時間を使われていることがわかる。

ここが速くなれば、あるいは、行列の作成時間を気にしなくていい状況とかであれば、十分使い物になるパフォーマンスになってきた。
C/C++にはとうていかなわないでしょうけど、Pythonでも工夫次第でいろいろがんばれるものだ。

以下、今日の改造版のソース。
#!/usr/local/python25/bin/python
# _*_ coding: utf-8 _*_
# Copyright (C) 2007 Ayukawa Hiroshi
import sys
import time
from operator import itemgetter
from collections import defaultdict
from itertools import repeat, groupby
from math import sqrt

def distmatrix_by_dict(vectors):
"""
上三角部分のみの距離表を作ります(ベクトルがdict型のsparseベクトルの場合)
返り値 (距離, ベクトル番号1, ベクトル番号2) の列
"""
mat = defaultdict(lambda: defaultdict(float))
for i in range(len(vectors)):
dist = []
for j in range(i):
d = 0.0
for k in set(vectors[i].keys()) & set(vectors[j].keys()): #ユークリッド距離
d += (vectors[i].get(k, 0.0) - vectors[j].get(k, 0.0))**2
dist.append((j, sqrt(d)))
mat[i] = dafaultdict(float, dist)
return mat

def distmatrix_by_list(vectors):
"""
上三角部分のみの距離表を作ります(ベクトルが通常のリストの場合)
返り値 (距離, ベクトル番号1, ベクトル番号2) の列
"""
mat = defaultdict(lambda: defaultdict(float))
for i in range(1,len(vectors)):
v = vectors[i]
dist = []
for j in range(i):
dist.append((j, sqrt(sum([(x[0]-x[1])**2 for x in zip(v, vectors[j])])))) #ユークリッド距離
mat[i] = defaultdict(float, dist)
return mat

def _initcache(mat):
"""
距離表の各行ごとの最小値をキャッシュする配列を初期化します。
"""
cache = []
for i in sorted(mat.keys()):
cache.append((i, min(((v, i, j) for j,v in mat[i].iteritems()))))
return dict(cache)

class _Rap:
"""
パフォーマンス測定用クラス
"""
def __init__(self):
self.t = 0
def b(self):
self.start = time.time()
def e(self):
self.t += time.time() - self.start

def cluster(vectors, distfunc=distmatrix_by_list):
"""
クラスタリング(ウォード法)
返り値 クラスタ結合履歴配列
"""
rap0 = _Rap()
rap00 = _Rap()

#距離表の用意
rap0.b()
mat = distfunc(vectors)
rap0.e()

#距離表の最小値キャッシュを初期化する
rap00.b()
cache = _initcache(mat)
rap00.e()

size = len(vectors) #全体のサイズ
sizes = list(repeat(1, size)) # 各クラスタのサイズ

#ノードを結合した履歴を保存する配列(この列を最終的に返します)
hist = [dict(
joined=None, #joined: この回で結合したインデックス
dist=0.0, #dist: この回で結合したインデックス間の距離
newid=None, #newid: この回で新設されたインデックス
)]
rap1 = _Rap()
rap2 = _Rap()
nextnodenum = size #次に使うクラスタ番号
kieta = set() #消えた行番号の保存用

#size-1回繰り返して結合していきます。
while len(hist) < size:
#最小距離のものを選ぶ
rap1.b()
dist, q, p = min(cache.values()) # p < q
rap1.e()

#消えるもの集合を更新
kieta.add(p)
kieta.add(q)

#高速化のための事前計算
pairsize = sizes[p]+sizes[q]
r = mat[q][p]
obj = mat[nextnodenum]
sp = sizes[p]
sq = sizes[q]

#距離の計算しなおし
rap2.b()
#キャッシュから、今回消えるものをあらかじめ排除
if p in cache: cache.pop(p)
if q in cache: cache.pop(q)
tmpdist = [] #新ノードからの最小距離を求めるための一時記憶場所
for i in xrange(nextnodenum):
if i in kieta: continue #すでに消えているなら何もしない

#距離の計算 Ward法
if i < p:
obj[i] = ((sp+sizes[i])*mat[p][i]+(sq+sizes[i])*mat[q][i]-sizes[i]*r)/float(pairsize+sizes[i])
elif p < i and i < q:
obj[i] = ((sp+sizes[i])*mat[i][p]+(sq+sizes[i])*mat[q][i]-sizes[i]*r)/float(pairsize+sizes[i])
mat[i].pop(p) #今後使わなくなる成分の消去
else:
obj[i] = ((sp+sizes[i])*mat[i][p]+(sq+sizes[i])*mat[i][q]-sizes[i]*r)/float(pairsize+sizes[i])
mat[i].pop(p) #今後使わなくなる成分の消去
mat[i].pop(q) #今後使わなくなる成分の消去

#キャッシュ更新
if i in cache:
if obj[i] < cache[i][0]: #新ノードとi番目の距離がより近い場合はi行目のキャッシュを置き換えます。
cache[i] = (obj[i], nextnodenum, i)
elif mat[i] and (p in cache[i][1:] or q in cache[i][1:]):
cache[i] = min(((v, i, j) for j,v in mat[i].iteritems()))
if not mat[i]: #i行が空になってたら、iのキャッシュは不要。matからもけずる
cache.pop(i)
mat.pop(i)
elif not mat[i]: mat.pop(i)
tmpdist.append((obj[i], i))
if p in mat: mat.pop(p) #今後使わなくなる成分の消去
if q in mat: mat.pop(q) #今後使わなくなる成分の消去

#新ノードからの最小距離をキャッシュに記録
if tmpdist:
mindist = min(tmpdist)
cache[nextnodenum] = (mindist[0], nextnodenum, mindist[1])

sizes.append(pairsize) #新クラスタのサイズデータを格納
rap2.e()
#履歴情報に今回の結合結果を格納
hist.append(dict(
joined=(p,q),
dist=dist,
newid=nextnodenum,
))
nextnodenum += 1 #次のクラスタID

print "rap0", rap0.t
print "rap00", rap00.t
print "rap1", rap1.t
print "rap2", rap2.t

return hist

def main(nottest=False):
"""
引数 nottest パフォーマンステスト用フラグ Falseの場合はテストのために、結果出力しません。
"""
ids = [] # データラベル格納用
vectors = [] # ベクトル格納用
for line in sys.stdin: #標準入力からベクトルを読む
data = line.strip().split("\t")
ids.append(data[0])
vectors.append(map(float, data[1:]))

#クラスタリング実行
hist = cluster(vectors)

if nottest: #出力処理
for c, h in enumerate(hist):
j = h["joined"]
if not j: continue
i = "(%s+%s)" % (ids[j[0]], ids[j[1]])
ids.append(i)
print c, ":", h["dist"], ":", i

if __name__ == "__main__":
import optparse

parser = optparse.OptionParser(u"""
Ward法によるクラスタリング。
""")
parser.add_option("-t", "--test", dest="test", help=u"パフォーマンステスト実行", default=False, action="store_true")
(options, args) = parser.parse_args()

if options.test:
#パフォーマンステスト
import hotshot, hotshot.stats
prof = hotshot.Profile("cluster.prof")
prof.runcall(main)
prof.close()
stats = hotshot.stats.load("cluster.prof")
stats.strip_dirs()
stats.sort_stats('time', 'calls')
stats.print_stats(20)
else:
#通常実行
main(True)

2008-05-14

Ward法再実装で少し速くなった

PythonでWard法によるクラスタリング
Ward法のパフォーマンスと分散処理のための考察

このシリーズの続き。
プログラムを書き直し、真っ当な行列的データを用意したところ、ずいぶん速くなり、588件でも19秒程度となった。(1100件でも2分程度なので、ずいぶんと自分が想定している実用化に近づいたと思う。)

しかし、真っ当な行列データにすると、今度は最小距離の探索に9割の時間を取られるようになったため、ここを改造すべきと強く思う。

ここまでくれば、常に行列の成分のソート情報をなんらか記憶しておけばいいように思うのだが、そのソート情報の維持にはまたコストがかかるというジレンマ。

今は可変長の行列として、pythonのdefaultdict(lambda: defaultdict(float))を使っている。これで定義した簡易行列クラスは非常に便利なのでそんなにシビアでない計算の目的にはよく使っている。可変長&メモリ節約というのがとくにすばらしい。

これがC++であったらどうするべきであろうか。
こうした細かい計算アルゴリズムの王道としては最初から使用するメモリをまとめて最小限確保しておき、その上でメモリを上書きしながら計算して行くべきか。
それとも今回のPythonでの実装にならってSTLを使って map<size_t, map<size_t, double>*>* などとするべきか。ああ、こういう場合shared_ptrにするべきなんだっけ、、、。ここはいつも覚えてないので参考書を読み直さなくては。
はたまたboostになんかいいのないかと見ていたら、いつの間にかboostのコンテナ系が増えている!調査しなければ。
うーん、前途多難。

以下、今日の改造版。

#!/usr/local/python25/bin/python
# _*_ coding: utf-8 _*_
# Copyright (C) 2007 Ayukawa Hiroshi
import sys
import time
from operator import itemgetter
from collections import defaultdict
from itertools import repeat, groupby
from math import sqrt

def distmatrix_by_dict(vectors):
"""
上三角部分のみの距離表を作ります(ベクトルがdict型のsparseベクトルの場合)
返り値 (距離, ベクトル番号1, ベクトル番号2) の列
"""
for i in xrange(len(vectors)):
for j in xrange(i):
dist = 0.0
for k in set(vectors[i].keys()) & set(vectors[j].keys()): #ユークリッド距離
dist += (vectors[i].get(k, 0.0) - vectors[j].get(k, 0.0))**2
yield (sqrt(dist), i, j)

def distmatrix_by_list(vectors):
"""
上三角部分のみの距離表を作ります(ベクトルが通常のリストの場合)
返り値 (距離, ベクトル番号1, ベクトル番号2) の列
"""
for i in xrange(len(vectors)):
for j in xrange(i):
dist = sum([(x[0]-x[1])**2 for x in zip(vectors[i], vectors[j])]) #ユークリッド距離
yield (sqrt(dist), i, j)

def minmatrix(mat, size):
d = min((min(((v, j) for j,v in mat[i].iteritems())), i) for i in mat.iterkeys() if mat[i])
return (d[0][0], d[1], d[0][1])

def cluster(vectors, distfunc=distmatrix_by_list):
"""
クラスタリング(ウォード法)
返り値 クラスタ結合履歴配列
"""
#距離表の用意
mat = defaultdict(lambda : defaultdict(float))
for d in distfunc(vectors):
mat[d[1]][d[2]] = d[0]

size = len(vectors) #全体のサイズ
sizes = list(repeat(1, size)) # 各クラスタのサイズ

#ノードを結合した履歴を保存する配列(この列を最終的に返します)
hist = [dict(
joined=None, #joined: この回で結合したインデックス
dist=0.0, #dist: この回で結合したインデックス間の距離
newid=None, #newid: この回で新設されたインデックス
)]

class Rap:
def __init__(self):
self.t = 0
def b(self):
self.start = time.time()
def e(self):
self.t += time.time() - self.start
rap1 = Rap()
rap2 = Rap()
nextnodenum = size #次に使うクラスタ番号
#size-1回繰り返して結合していきます。
kieta = set()
while len(hist) < size:
#最小距離のものを選ぶ
rap1.b()
dist, q, p = minmatrix(mat, nextnodenum) # p < q
rap1.e()

#消えるもの集合を更新
kieta.add(p)
kieta.add(q)

#高速化のための事前計算
pairsize = sizes[p]+sizes[q]
r = mat[q][p]
obj = mat[nextnodenum]
sp = sizes[p]
sq = sizes[q]

#距離の計算しなおし
rap2.b()
for i in xrange(nextnodenum):
if i in kieta: continue #すでに消えているなら何もしない

#距離の計算 Ward法
if i < p:
obj[i] = ((sp+sizes[i])*mat[p][i]+(sq+sizes[i])*mat[q][i]-sizes[i]*r)/float(pairsize+sizes[i])
elif p < i and i < q:
obj[i] = ((sp+sizes[i])*mat[i][p]+(sq+sizes[i])*mat[q][i]-sizes[i]*r)/float(pairsize+sizes[i])
mat[i].pop(p) #今後使わなくなる成分の消去
else:
obj[i] = ((sp+sizes[i])*mat[i][p]+(sq+sizes[i])*mat[i][q]-sizes[i]*r)/float(pairsize+sizes[i])
mat[i].pop(p) #今後使わなくなる成分の消去
mat[i].pop(q) #今後使わなくなる成分の消去
if p in mat: mat.pop(p) #今後使わなくなる成分の消去
if q in mat: mat.pop(q) #今後使わなくなる成分の消去
sizes.append(pairsize) #新クラスタのサイズデータを格納
rap2.e()
#履歴情報に今回の結合結果を格納
hist.append(dict(
joined=(p,q),
dist=dist,
newid=nextnodenum,
))
nextnodenum += 1 #次のクラスタID
print "rap1", rap1.t
print "rap2", rap2.t

return hist

def main(nottest=False):
"""
引数 nottest パフォーマンステスト用フラグ Falseの場合はテストのために、結果出力しません。
"""
ids = [] # データラベル格納用
vectors = [] # ベクトル格納用
for line in sys.stdin: #標準入力からベクトルを読む
data = line.strip().split("\t")
ids.append(data[0])
vectors.append(map(float, data[1:]))

#クラスタリング実行
hist = cluster(vectors)

if nottest: #出力処理
for c, h in enumerate(hist):
j = h["joined"]
if not j: continue
i = "(%s+%s)" % (ids[j[0]], ids[j[1]])
ids.append(i)
print c, ":", h["dist"], ":", i

if __name__ == "__main__":
import optparse

parser = optparse.OptionParser(u"""
Ward法によるクラスタリング。
""")
parser.add_option("-t", "--test", dest="test", help=u"パフォーマンステスト実行", default=False, action="store_true")
(options, args) = parser.parse_args()

if options.test:
#パフォーマンステスト
import hotshot, hotshot.stats
prof = hotshot.Profile("cluster.prof")
prof.runcall(main)
prof.close()
stats = hotshot.stats.load("cluster.prof")
stats.strip_dirs()
stats.sort_stats('time', 'calls')
stats.print_stats(20)
else:
#通常実行
main(True)

2008-05-13

2ちゃんねる Pythonのお勉強 Part26よりParallel Python関係の議論抜粋

Pythonのお勉強 Part26
からの抜粋。なんともタイムリーなことに、Parallel Pythonが話題になっている。
この議論で288が指摘している事は、やはり的を得ているのではないだろうか。
Python->C/C++という移植で10倍速くなるなんてことは珍しくない。
(例えば、C++からMecabを利用した時の速度は驚異的ですらあるが、Pythonから利用した場合はちょっとのんびりとしたものになってしまう。)
numpyがうまくハマっているPythonスクリプトとかなら別ですけど。

で、以下その議論抜粋です。。。。
---------------------

169 :デフォルトの名無しさん:2008/05/06(火) 21:31:30
threadって計算速度の向上に効果ある?
a = b1 + b2の計算を(a,b1,b2はarray)
b1とb2をthreadで計算して,
それぞれ計算終了後に足すというプログラムを作ったんだけど,
普通にb1とb2をメインスレッドで順番に求めて足した方が早かった...

ちなみにb1とb2の計算はダミーのforループ10000回です.

244 :デフォルトの名無しさん:2008/05/10(土) 04:27:34
>>169
Parallel Python というパッケージがあるのを知ったので参考までに紹介しとく。
http://www.parallelpython.com/
スレッドじゃなく複数プロセスで並列実行する仕組みらしい。
Python のみで実装されていて非常にシンプルなパッケージ構成になっている。
マルチコア環境の場合はスレッドモジュール風に使える。
クラスタ環境では各ノードで計算サーバ(ppserver.py)をしておく仕組みになっている。
応用プログラムはどちらの環境でも同じにできるっぽい。

手元の共有メモリマシンで付属サンプル(sum_primes.py)を実行したら下のようになった。
素数の和を計算するプログラムなのでプロセス間通信はほとんどないけど、
並列化のオーバヘッドはあるわけで、なかなか良好な結果だと思う。

並列度,実行時間(秒),速度向上率
1,  17.54,  1.000
2,  8.781,  1.998
4,  4.424,  3.965
8,  2.261,  7.759

245 :169:2008/05/10(土) 13:53:14
>>244
おおお,サンクス!
早速いじってみまふ。

280 :244:2008/05/10(土) 23:51:28
引続き Parallel Python をいじってみたんだけども、このシステムはマスタワーカモデルを前提にしているみたいだ。
つまり、マスタが仕事の集合を持っていて、有限個のワーカに1つずつ仕事を割り当てる。仕事を終えたワーカは
マスタから新しい仕事をもらって処理する。これを仕事の集合が空になるまで続ける。

ワーカ間で通信をする機能は提供されていないし、特定のワーカ(たとえば特定のリモートホスト上のppserver.py)に
特定の仕事を割り当てることもできないっぽい。Parallel Python が有効かどうかは実現したい並列アルゴリズムが
マスタワーカモデルに適合するかどうかに依る希ガス。

283 :デフォルトの名無しさん:2008/05/11(日) 00:40:31
>>280
具体的には、どういう用途に使えそうなんですか?

285 :デフォルトの名無しさん:2008/05/11(日) 01:34:00
>>283
最も適しているのは並列処理の分野で embarrassingly parallel と呼ばれるカテゴリに属する種々の計算、
すなわち入力データを複数の独立した部分に分割できて、各部分について他から独立して計算が行なえる
タイプの用途に適している。
http://en.wikipedia.org/wiki/Embarrassingly_parallel
Wikipediaに例があがっている。フラクタル計算とか3DCGのレンダリング(例:レイトレーシング)とか力任せの
暗号解読とか色々ある。

並列処理って概して複雑で長時間計算し続ける必要のある応用が多い。そういう応用には C とか Fortran を
使うことがほとんどなんだけど、きちんと動くようになるまでの開発時間が長くてデバッグがたいへんだったりする。
そういう並列プログラムの試作(プロトタイピング)に Python が使えたらいいなーと個人的には思っている。

288 :デフォルトの名無しさん:2008/05/11(日) 09:43:58
>>285
PythonでCPU8個使って並列計算するよりCで書いたプログラムの方が速そうだったり。。。

289 :デフォルトの名無しさん:2008/05/11(日) 09:57:09
>>288
何言ってんだ?

290 :デフォルトの名無しさん:2008/05/11(日) 10:28:52
>>289
あ?ヤンのかコラ


295 :デフォルトの名無しさん:2008/05/11(日) 14:39:25
>>285
自分も並列化するときのプロトタイプとしてPythonの並列環境を使ってみたいと
考えているんだけれど、今のところプロトタイプとしての感触はどう?

データ分割の楽な問題ってCとかFORTRANでもOpenMPを注意深く使うだけでも早くなる事が
多いし、MPIでも慣れれば苦労せずに書けるけれど、そうでない、スレーブ間での通信が必要
とかうまく分割できないとかそういう問題はMPIでやろうとするとデバッグが(;゚д゚)
なことになりがちで、そういう用途にPythonのプロトタイプが役に立つなら素晴らしいと思う
実はC並列版とそれほど速度とメモリ使用量が変わらないとかなら最高
C並列版より速ければ神

296 :デフォルトの名無しさん:2008/05/11(日) 16:53:12
Parallel Python だけど pp.Server.__scheduler() を適当に書き直せばワーカと仕事の対応(割り当て)を
ユーザ側で決められそうな希ガス。問題はワーカ間のプロセス間通信。これが一番面倒なところなんで
自前で実装となると Parallel Python を使ううまみがほとんどない・・・。

やっぱ MPI の Python バインディングあたりが一番現実的な解かなー。でもなんか気が重いんだなー。
もっと Python らしく lightweight なソリューションがないかなー。

>>288
Python の並列プログラムより C の逐次プログラムの方が速そうってことだよね。
そういうことは十分(多分頻繁に)あると思われ。

>>295
残念ながらまだ並列プログラムのプロトタイプ用途には使えてなくて感触を得るところまでいってない。
理由は単純で、Python で手軽に並列プログラミングを実現できる道具を見つけられていないから。
>>209に書いたようにいろいろ試してるんだけどなかなか・・・。

ただ、個人的には Python を並列化のプロトタイピングに使うのは大いに有望だと思っている。
プロトタイピングの場合、欲しいのは並列度やデータ量を上げたときの実行時間の変化であって、
実行速度が多少遅くて実験に時間がかかるとしても知りたいことは分かるはずだから。
プロセス間通信にソケットを使うとすると、データ量が大きくなれば Python でも C 等と同じぐらいの
速度が出る(ハズ)。演算量が多い部分に numpy 等の C/Fortran で書かれた数値カーネルを使うことに
すれば、プロトタイピング用途には十分な速度とメモリ使用量が得られるのではと思う。

Ward法のパフォーマンスと分散処理のための考察

昨日実装したWard法の処理だが、200件程度なら2秒で終わるが、詳しく調べると、400件で10秒、700件で60秒と、O(n**2)どころか、「もしかしてO(2^n)?!」ぐらいのできであった(T_T)。。。

その99%は距離計算直し対象の発見と選り分けの部分であり、本来のWard法の距離計算の部分は微々たるもの。

今後はこの選り分け対象の部分でちゃんとインデックスの効いた選り分けができるように変更すべきか。

やはり、まずは王道に従って、きちんとした行列データとして距離を持つべきだろう。

作り直しだ。

2008-05-12

PythonでWard法によるクラスタリング

Pythonで、Ward法によるクラスタリング(デンドログラム作成まで)を実装してみた。
参考にしたのは、岡山大学が公開しているこの資料
(岡山大学の先生ありがとう!)

今回こんな車輪の再発明をしたのは、距離表と呼ばれる行列データを行列のイメージで保存せず、

[(ベクトル1とベクトル2の距離, ベクトル番号1, ベクトル番号2), ...]

という形式のリストにして、距離でソートしておき、そのソート結果を維持したままこのリストを分割、更新していくという方針で作ったらどうかと思ったからです。

196個の19次元の非スパースなデータのデンドログラム作成に2.5秒だからなかなかのパフォーマンスではないだろうか。
CもC++もNumpy系の数値計算ライブラリも使わないでこれだから、いいですねえ。

使い方
読み込むデータファイルをタブ区切りで用意する。
1カラム目はデータの識別用ID文字列、2カラム目以降にベクトルの数値データを記入。
あとはそのデータを標準入力から流し込むだけ。

以下のsample.txtを用意する。
A 1.0 1.1 0.2
B 1.1 0.5 0.3
C 3.1 0.5 0.2
D 4.0 1.1 1.2

以下を実行
ayu@~/work/cluster% ./cluster.py < sample.txt 
1 : (B+A)
2 : (D+C)
3 : ((D+C)+(B+A))

最初にBとAを結合し、次にDとCを結合し、最後に(D+C)と(B+A)を結合していることが分かる。

ソースcluster.py
#!/usr/local/Python25/bin/python
# _*_ coding: utf-8 _*_
# Copyright (C) 2007 Ayukawa Hiroshi
import sys
from collections import defaultdict
from itertools import repeat, groupby
from math import sqrt

def distmatrix_by_dict(vectors):
"""
上三角部分のみの距離表を作ります(ベクトルがdict型のsparseベクトルの場合)
返り値 (距離, ベクトル番号1, ベクトル番号2) の列
"""
for i in xrange(len(vectors)):
for j in xrange(i):
dist = 0.0
for k in set(vectors[i].keys()) & set(vectors[j].keys()): #ユークリッド距離
dist += (vectors[i].get(k, 0.0) - vectors[j].get(k, 0.0))**2
yield (sqrt(dist), i, j)

def distmatrix_by_list(vectors):
"""
上三角部分のみの距離表を作ります(ベクトルが通常のリストの場合)
返り値 (距離, ベクトル番号1, ベクトル番号2) の列
"""
for i in xrange(len(vectors)):
for j in xrange(i):
dist = sum([(x[0]-x[1])**2 for x in zip(vectors[i], vectors[j])]) #ユークリッド距離
yield (sqrt(dist), i, j)

def cluster(vectors, distfunc=distmatrix_by_list):
"""
クラスタリング(ウォード法)
返り値 クラスタ結合履歴配列
"""
#距離表の用意
data = list(distfunc(vectors))
data.sort()

size = len(vectors) #全体のサイズ
sizes = list(repeat(1, size)) # 各クラスタのサイズ

#ノードを結合した履歴を保存する配列(この列を最終的に返します)
hist = [dict(
data=data, #data: まだ結合されていない成分一覧(常に距離でソート済みの状態で格納)
joined=None, #joined: この回で結合したインデックス
dist=0.0 #dist: この回で結合したインデックス間の距離
)]

nextnodenum = size #次に使うクラスタ番号

#size-1回繰り返して結合していきます。
while len(hist) < size:
#最小距離のものを選ぶ
firstpair = min([h["data"][0] for h in hist if h["data"]])
#最小距離のもののインデックスペア
pairindex = set(firstpair[1:])

movedata = [] #新クラスタ形成のために計算し直し対象とするもの格納用
for h in hist: #履歴をさかのぼって、計算しなおしとそうでないものに分ける。
newdata = [] #計算しなおさない居残り組格納用
for d in h["data"]:
if d == firstpair: #今回結合対象のものは除外
continue
intersect = set(d[1:]) & pairindex
if intersect: #今回結合対象とインデックスがかぶるなら、計算しなおし対象へ
movedata.append(d)
else: #そうでないなら居残り組
newdata.append(d)
h["data"] = newdata #居残り組に更新

newdata = [] #計算しなおしデータ格納用

dpq, p, q = firstpair #改めて今回の結合ペアをp,qなどと名付ける
#計算し直しデータを点ごとにグループ化。この点ごとに新距離を計算しなおします。
movedata = groupby(sorted((((set(x[1:])-pairindex).pop(), x) for x in movedata)), key=lambda x: x[0])
for r, d in movedata:
newdist = 0.0 #新距離計算用
for w in d: #ward法計算式にしたがって計算
w = w[1]
newdist += (sizes[w[1]]+sizes[w[2]])*w[0]
newdist = (newdist-sizes[r] * dpq) / float(sizes[r]+sizes[p]+sizes[q])
newdata.append((newdist, nextnodenum, r)) #新距離を格納
newdata.sort() #距離データは常にソートして整理整頓
sizes.append(sizes[p]+sizes[q]) #新クラスタのサイズデータを格納
nextnodenum += 1 #次のクラスタID
#履歴情報に今回の結合結果を格納
hist.append(dict(
data=newdata,
joined=(p,q),
dist=firstpair[0]
))
return hist

def main(nottest=False):
"""
引数 nottest パフォーマンステスト用フラグ Falseの場合はテストのために、結果出力しません。
"""
ids = [] # データラベル格納用
vectors = [] # ベクトル格納用
for line in sys.stdin: #標準入力からベクトルを読む
data = line.strip().split("\t")
ids.append(data[0])
vectors.append(map(float, data[1:]))

#クラスタリング実行
hist = cluster(vectors)

if nottest: #出力処理
for c, h in enumerate(hist):
j = h["joined"]
if not j: continue
i = "(%s+%s)" % (ids[j[0]], ids[j[1]])
ids.append(i)
print c, ":", i


if __name__ == "__main__":
import optparse

parser = optparse.OptionParser(u"""
Ward法によるクラスタリング。
""")
parser.add_option("-t", "--test", dest="test", help=u"パフォーマンステスト実行", default=False, action="store_true")
(options, args) = parser.parse_args()

if options.test:
#パフォーマンステスト
import hotshot, hotshot.stats
prof = hotshot.Profile("cluster.prof")
prof.runcall(main)
prof.close()
stats = hotshot.stats.load("cluster.prof")
stats.strip_dirs()
stats.sort_stats('time', 'calls')
stats.print_stats(20)
else:
#通常実行
main(True)

2008-05-11

Parallel Pythonで分散処理

Erlangでなくて、Pythonで分散処理を書く意義はなんだろうかと考えた。
Erlangはたしかに分散処理が得意なんだろうけど、いろいろ調べた感じでは、複雑な数値計算などの分散処理には向いていないというウワサだ。良く知られているTwitterや通信の例のようなシンプルな処理を膨大な量さばくにはいいようだけど。

計算を分散で行う場合、本来はGoogleが採用しているように、C++をベースにすべきだろう。
(ただし、GoogleではSawzallという独自言語で記述し、それをC++に変換して実行するそうな。)

そうなると、「なぜPythonで分散処理?」というのが重要になる訳だが、おそらく以下のようなことだろうか。

  1. 既存の豊富なモジュール(しかも多くはCで書かれている)を使える。
  2. C、C++で書かれたルーチンをswig等でPythonに連結し、それを分散させれば、実質C、C++で実行しているようなものだろう。
  3. 簡単で読みやすい。(重要!)

そこで、まさに自分が今仕事で考えているものの要件でぴったりの問題があった。
問題
文の列(数十文程度)の各文どうしの総当たりのdiff(通常の行単位ではなく、
単語単位でのdiff)を取り、diffで異なった割合を数値にした総当たり表を
リアルタイムで返す。

都合のいいことに、Pythonには標準でdiffを取るためのライブラリdifflibがついているので、こむずかしいdiff計算はそいつに丸投げできる。
これが、文の数が10、20程度なら分散させずに普通に解いても速度上問題はない。
しかし、これが40、50になってくると、ちょっと「リアルタイム」とは言えなくなってくる。
(文の数をnとすると、対称であることと、自分自身とのdiffは計算しないとしても、n*(n-1)/2回diffを取る必要がある。)

ということで、これを分散処理させてみよう。

準備
Parallel Pythonを使う。インストールは通常のpython setup.py install 一発で完了。
実験に使ったソースdiffdist.pyは本記事末尾に添付した。
以下のように実行できる。
分散処理で働いてもらいたいマシン群で、Parallel Pythonのサーバーを起動しておく。
ppserver.pyコマンドはParallel Pythonをインストールした時に、pythonコマンドと同じディレクトリ(Linux, OSX)又は、c:¥Python25¥Scritpsフォルダ(Windows)に格納されているものです。
$ /PATH/TO/PYTHONBIN/ppserver.py

あとは、diffdist.pyのppservers変数(113行あたり)にサーバー群のアドレスを記入しておいて、
$ ./diffdist.py

で実行。
(※Parallel Python本家ドキュメントに従ってauto discovery機能を使えば、IPアドレスの設定記入を省く事もできる。)
今回は、140文(合計サイズ32kb程度)の英文の総当たり表を計算した。(今回このデータは添付しない。)
英文データファイルの各行どうしで、単語ごとのdiffを取るというサンプルになっている。

結果
まず普通に分散させないコードだとどれぐらいになるかというと、
ayu@~/work/difftest% time ./diffdist.py 
Starting pp with 2 workers
./diffdist.py 26.64s user 0.36s system 97% cpu 27.705 total

ていう感じ。余裕で20秒以上かかります。

これを、自分の手元マシン(MacBook白2GHz Intel Core2 Duo & メモリ2G)でppserverを起動して分散処理させると、
ayu@~/work/difftest% time ./diffdist.py 
Starting pp with 2 workers
./diffdist.py 14.22s user 0.62s system 90% cpu 16.412 total

ここに、さらに別マシン(DELL SC1420、Intel Xeon 2.8GHz & メモリ1G)を追加すると、、、
ayu@~/work/difftest% time ./diffdist.py 
Starting pp with 2 workers
./diffdist.py 9.62s user 0.59s system 89% cpu 11.394 total

順調に速くなっている。。
おまけだが、我が家の物置で眠っていたSharpのメビウス君(AMD Duron 897MHz & メモリ256M)を追加すると、、、
ayu@~/work/difftest% time ./diffdist.py 
Starting pp with 2 workers
./diffdist.py 8.03s user 0.58s system 85% cpu 10.100 total

ほんのちょっとだけど速くなるね!

この調子で、どんどんマシンを追加していけば、この処理はすごく速くなっていくだろう。
明日会社で本格的に試してみよう。

ソース diffdist.py
#!/usr/local/Python25/bin/python
# _*_ coding: utf-8 _*_
# Copyright (C) 2007 Ayukawa Hiroshi
import sys
import difflib
from itertools import islice
from collections import defaultdict

import pp

def getdiffrate(src1, src2):
"""
diffを取って、異なった割合(diff割合)を返す。
"""
SPC = " "
same = 0
diff = 0
for x in difflib.ndiff(src1, src2):
if x[0] == SPC:
same += 2 #一致部分カウントはsrc1とsrc2の分で+2
else:
diff += 1
return 100 * same / float(same+diff)

def diffrow(src):
"""
一つの文と、対比する複数の文のdiff割合をバッチ的に計算する。

引数 src
[(文番号, 原文, [(対比する文1の番号,対比する文1), (対比する文2の番号,対比する文), ..]), ...]

返り値
[(文番号, 対比する文の番号, 原文と対比文のdiff割合), ...]
"""
ans = []
for i1, s1, osrc in src:
for i2, s2 in osrc:
rate = getdiffrate(s1, s2)
ans.append((i1, i2, rate))
return ans

def d_diffmatrix(job_server, srcs, unit=140):
"""
文一覧を与えて、各文同士の間のdiff割合行列を返す。(分散処理バージョン)

引数
job_server 分散処理サーバー
srcs 文一覧(文字列の配列)
unit バッチ処理する単位(文の数)

返り値
行列データ(dict型データで、res[1][2]のようにdiff割合を格納しています。)
ただし、対角線成分上には、同じ文同士のdiff割合を格納すべきですが、計算を省略するため、
0.0を格納しています。
"""
#各タスク実行に必要なデータを作成する。
task = [] #タスクに与えるデータ格納用
subtask = [] #サブタスクとりわけ用
sublen = 0 #サブタスクの分量計算用
#全タスクを、文の数に従って、サブタスクに分割します。
#サブタスクごとに分散処理させます。
for i1, s1 in enumerate(srcs):
osrc = list(islice(enumerate(srcs), i1))
subtask.append((i1, s1, osrc))
sublen += len(osrc)
if sublen > unit: #文の数がunitを超えたら、一回分としてtaskに格納。
task.append(subtask)
subtask = []
sublen = 0
if subtask:
task.append(subtask)

#分散処理します。
jobs = [] #実行ジョブ格納用
for t in task:
jobs.append(job_server.submit(diffrow,(t, ), (getdiffrate,), ("difflib",)))

#結果を取得して格納
diffrates = defaultdict(lambda: defaultdict(float))
for job in jobs:
for i1, i2, rate in job():
#行列状に格納
diffrates[i1][i2] = rate
diffrates[i2][i1] = rate
return diffrates

def diffmatrix(srcs):
"""
文一覧を与えて、各文同士の間のdiff割合行列を返す。(比較のための通常処理バージョン)

引数
srcs 文一覧(文字列の配列)

返り値
行列データ(dict型データで、res[1][2]のようにdiff割合を格納しています。)
ただし、対角線成分上には、同じ文同士のdiff割合を格納すべきですが、計算を省略するため、
0.0を格納しています。
"""

#全タスクのジェネレーター
def gen():
for i1, s1 in enumerate(srcs):
osrc = list(islice(enumerate(srcs), i1))
yield (i1, s1, osrc)
#処理しながら、結果を行列に格納
diffrates = defaultdict(lambda: defaultdict(float))
for i1, i2, rate in diffrow(gen()):
diffrates[i1][i2] = rate
diffrates[i2][i1] = rate
return diffrates

if __name__ == "__main__":
ppservers = ("192.168.1.21","192.168.1.22","192.168.1.23")

if len(sys.argv) > 1:
ncpus = int(sys.argv[1])
job_server = pp.Server(ncpus, ppservers=ppservers)
else:
job_server = pp.Server(ppservers=ppservers)

print "Starting pp with", job_server.get_ncpus(), "workers"

src = file("samplepatent.txt").readlines()
res = d_diffmatrix(job_server, [x.split(" ") for x in src])

#比較用
#res = diffmatrix([x.split(" ") for x in src])

2008-05-08

pythonで英語の同義語などを取得する

nltkにはwordnetが組み込まれているので、nltkひとつをインストールするだけで、同義語、上位語、下位語などを自由に取得できる。

>>> from nltk import wordnet
>>> dog = wordnet.N["dog"]
>>> dog.synsets()
[{noun: dog, domestic_dog, Canis_familiaris}, {noun: frump, dog},
{noun: dog}, {noun: cad, bounder, blackguard, dog, hound, heel},
{noun: frank, frankfurter, hotdog, hot_dog, dog, wiener,
wienerwurst, weenie}, {noun: pawl, detent, click, dog},
{noun: andiron, firedog, dog, dog-iron}]


英語の検索エンジンの検索クエリ入力画面を作ったときなんかに、このデータでクエリーを膨らませて検索することとかも考えてみると使い手がありそうですね。

より詳しい使い方はnltkのデモでざっと眺める事ができます。

nltkのインストール方法は、自分の過去記事にあります。

2008-05-04

日本語テキストのトピック分割

先日からcodereposに置いている自動要約モジュールに、日本語テキストのトピック分割のソースをコミットしました。
(->そのソース)

このスクリプトでは、与えられた日本語のテキストを、トピックごとに分割する機能を提供しています。

基本的には論文"Advances in domain independent linear text segmentation"を参考にしています。
この論文では
1. 文ごとにTFベクトルを計算し、
2. そのベクトル間でcos正規化された内積を計算して文間の類似度を算出、
3. 近接する文同士の類似度の変化具合を見て、トピックの変わり目を決定。
という方式をとっています。
ですが、今回の実装では、上記を日本語にも適用するためにさらに以下の改良を加えています。

動詞も使うことにした。


日本語のテキスト、特にブログの記事などは、名詞だけでなく、「節約する」とか「怒る」とかの動詞も考慮にいれたほうがよかろうと思い、動詞も使っています。ただし、動詞は活用しますので、mecabでの処理結果の品詞情報の行から動詞の原型を取得して使うことにしました。
(->対応するソースの部分

表記揺れ対策


原論文では英語のporter stemmingをしていましたが、日本語では独自の表記揺れ対策をしなければいけません。
表記揺れ対策には以下の方針をとりました。

1. カタカナ語は地道に辞書を作る。
カタカナ語の表記揺れ辞書データは本職で半自動生成したことがありますので、近いうちに作ろうと思います。大量の生テキストデータが必要ですので、どこかからテキストを取得せねばなりません。wikipediaデータでもいいかもしれませんが、できればあんなきれいな文章ではないほうがいいので考えます。
(->対応するソースHyokiyureData.py)

2. 漢字の単語はunigramも生成する
名詞&未知語の場合のみ、漢字の単語の表記揺れは、その語のunigramもTFカウントすることで対応します。
例えば、通常「人」と「人々」は異なる語と認識されてしまいますが、「人々」の方を文字単位で分割し、「人」と「々」それぞれについてもカウントすることにより、「人」と「人々」が互いに多少は類似していると判定できるようになります。
(->対応するソース)

語の重みのバリエーション


語の重みを以下のように変えて計算しています。
* 名詞、未知語: 1語 = 重み1
* 名詞、未知語unigram: 1語(1文字) = 重み0.3
* 動詞: 1語 = 重み0.5
これは、一文でTFベクトルを作っても、たいていはTF=1の語ばかりで、あまりその文の特徴を表しているとも思えませんでしたので、unigramや動詞も考慮に入れた上で、それらの重みの間に差を付けることで、文の間により特色の違いをだせないかという試みです。検証してませんが、まあ、見たかんじ結果がよくなったように思うので、ひとまずこれで行きます。
(->対応するソース)

TFベクトルではなく、TFIDF的ベクトルを使用


ブログ記事のような中程度の長さの文を扱うことを想定していますので、より文の間の特徴を鋭く認識したほうがいいかと思いました。
そこで、文全体でまんべんなく頻出している単語の効果を間引くために語が現れた文の数の逆数をかけることで、ベクトルの重みを調整しています。
(->対応するソース


まだまだ先は長いですが、こつこつ作って行きます。

2008-05-01

Jython本を注文した

Jython本を注文しました。

まだAmazonのレビューはない見たいだけど、ブログでは話題になっているようだ。

随所で最高の評価を受けているようなので、とても楽しみです!

今度「Python禁止!JAVAだけで書け!」みたいなプロジェクトに無理矢理入れられそうなので、ちょっとJythonの準備しとくかな。。。

2008-04-30

Python用の自動要約モジュールをcodereposにimportしたよ

けっこう前からのんびり作っていた自動要約モジュールをCodereposにimportしました。

http://coderepos.org/share/browser/lang/python/yoyaku

この自動要約モジュールはまだあまり精度はよくないです。
が、どのぐらい良くないかを知るには、以下のサイトでいちおうサービスとして設置していますので、お試しいただけます。

自動要約サービス Civory

今の所は単純なMMRアルゴリズムで文を抜粋するだけのものですので、まあ、あまり、、、。

今取り組んでいるのは、英文の自動要約の際に、文をまずトピックごとに分割して、その上で要約をかけようとしています。
トピックに分けるアルゴリズムは以下の論文を参考にしています。
Advances in domain independent linear text segmentation
いちおうこれの実装はひとまず終わっていて、それに対応するソースは以下になります。
TopicSegmentationEN.py
今後これにバグがないか注意しながらいろんな例で実験を重ねていく予定です。

NLTKでstemming

nltkにはステミングの有名なアルゴリズムであるporterアルゴリズムが用意されている。

>>> import nltk
>>> nltk.PorterStemmer().stem("application")
'applic'
>>>

2008-04-20

NLTKで英文の文末判定

英文の文末を判定する簡易なルールベースのアルゴリズム。

Separating sentences

1年前にこの手のアルゴリズムを実装しようとしたが、この問題は非常にやっかいです。
たとえばこんな例:

CELLULAR COMMUNICATIONS INC. sold 1,550,000 common shares at $21.75 each
yesterday, according to lead underwriter L.F. Rothschild & Co.

"INC."の直後や"$21.75"、"L.F."などのピリオドを文末と認識しては大間違いになるのです。

この問題を解決するのに自分が1年前に着目していた論文は以下のもの。

Unsupervised Multilingual Sentence Boundary Detection

この論文では、特に言語を英語だけに限定しない方法を提案しています。
大規模な生のテキストデータから得られる統計情報のみで、文末判定を行えます。
難点はルールベースなどとは違って、事前の綿密な統計の作成、統計処理後の各種特別処理の実装が面倒くさいこと。
昨年途中まで実装したが、Python用の自然言語処理ライブラリNLTKでどうも実装予定との情報を見つけて半端でやめていました。

それで久しぶりに調べてみたら、このアルゴリズムがすでにNLTKで実装されて公開されているではないですか!

早速インストールして使ってみました。

インストール
% wget http://prdownloads.sourceforge.net/nltk/nltk-0.9.2.tar.gz
% tar zxvf nltk-0.9.2.tar.gz
% cd nltk-0.9.2
% sudo python setup.py install
% cd ..
% sudo mkdir /usr/share/nltk
% cd /usr/share/nltk
% sudo wget http://prdownloads.sourceforge.net/nltk/nltk-data-0.9.2.zip
% sudo unzip nltk-data-0.9.2.zip
% sudo chmod -R g+r data
% export NLTK_DATA=/usr/share/nltk/data
% python
>>> import nltk
>>> nltk.corpus.brown.words()
['The', 'Fulton', 'County', 'Grand', 'Jury', 'said', ...]
>>>

ここまでできればインストール完了。

生テキストデータを用意
例えば、Google News(英語版)などからリンクをたどって、ひたすらニュース本文をファイルnews.txtにコピペする。
実験でやるんであれば1000行程度で十分でした。

実験!!
まず生テキストを食わせて学習
>>> from nltk.tokenize import PunktSentenceTokenizer
>>> p = PunktSentenceTokenizer()
>>> fp = file("news.txt")
>>> p.train(fp.read())


次に実際に判定してみる。
判定に使うテキストは次のような少々意地の悪い例
The Finland-based company expects a weaker dollar and slower economic growth in the U.S. and parts of Europe to dampen the overall handset market this year. About half of Nokia's (NOK) sales are in dollars or currencies tied to it; a weaker dollar makes imports more expensive.

"What spooked us was its outlook for the industry in general," said Rick Franklin, equities analyst at Edward Jones.

Nokia reiterated projections that the industry shipments of handsets will grow 10% this year over last. In the first quarter, though, global shipments rose 17%, suggesting a slowdown in the remainder of the year.

For the quarter that ended March 31, Nokia earned $1.9 billion (1.2 euros), up 25% from the same quarter last year but short of an expected $2.3 billion. Overall sales rose 28% to $20.1 billion (12.6 billion euros), roughly in line with views.


>>> a=p.tokenize("""The Finland-based company expects
a weaker dollar and slower economic growth in the U.S. and parts of Europe to
dampen the overall handset market this year. About half of Nokia's (NOK) sales
are in dollars or currencies tied to it; a weaker dollar makes imports more
expensive.

"What spooked us was its outlook for the industry in general," said Rick Franklin,
equities analyst at Edward Jones.

Nokia reiterated projections that the industry shipments of handsets will grow 10%
this year over last. In the first quarter, though, global shipments rose 17%,
suggesting a slowdown in the remainder of the year.

For the quarter that ended March 31, Nokia earned $1.9 billion (1.2 euros), up 25%
from the same quarter last year but short of an expected $2.3 billion. Overall
sales rose 28% to $20.1 billion (12.6 billion euros), roughly in line with views.""")
>>> for x in a:
... print x
... print "-"*20
The Finland-based company expects a weaker dollar and slower economic growth
in the U.S. and parts of Europe to dampen the overall handset market this year.
--------------------
About half of Nokia's (NOK) sales are in dollars or currencies tied to it;
a weaker dollar makes imports more expensive.
--------------------
"What spooked us was its outlook for the industry in general," said Rick Franklin,
equities analyst at Edward Jones.
--------------------
Nokia reiterated projections that the industry shipments of handsets will
grow 10% this year over last.
--------------------
In the first quarter, though, global shipments rose 17%, suggesting a slowdown
in the remainder of the year.
--------------------
For the quarter that ended March 31, Nokia earned $1.9 billion (1.2 euros), up
25% from the same quarter last year but short of an expected $2.3 billion.
--------------------
Overall sales rose 28% to $20.1 billion (12.6 billion euros), roughly in line
with views.
--------------------
>>>

きちんと、ピリオドで文をわけつつも、"U.S."や、"1.2 euros"などで区切るのは避けていることが分かります。

精度を上げるには、もっともっと大量の生テキストを食わせる必要があります。

2008-04-10

Google App Engineが待ち遠しすぎる。。。

TekTekGoogle App Engine対応は完了したのだが、まだ登録の順番が回ってこない!!

なので、ちょっとメモなど

* Google App Engineではdjangoのauth、sessionは使えないようだ。sessionが使えないのは、djangoのsessionがdjangoのmodelに依存しているからだそうな。この制限は厳しいが誰かが簡単な解決を探してくれるんじゃないかな。
* TekTekでは、sessionをやめて、できる限りCookieに保存するようにした。