2008-05-24

Djangoでblogっぽいのを作ったので移動します

Djangoでblogっぽいのを作ったので移動します。RSSも吐いておりますので宜しくお願いします。
新TekTekBlog

移動の理由:
1. このBloggerでは一切のJavascriptなどが貼れないので不便。アクセス解析の類も一切無理。
2. データのエクスポート機能もないので、書き溜めても外部に効率的に持ち出せない。
3. 自分で作って自分で設置するほうが面白い。

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-20

C++のWard法->g++のオプションで驚異的な高速化

先日から作っているクラスタリングアルゴリズムの一種Ward法のプログラム今日もいじっていた。

こまごまとしたことをしたのだが、コンパイルオプションをいじると、恐ろしく速くなった。

まず、比較のために、前回の結果。

件数   処理時間(s)    そのうち初期行列作成時間
--------------------------------------------
2352 33.432 6.39

そこに、今回以下のコンパイルオプションをつけてみた。(ちなみにOSX10.4&Intelな白MacBookの環境です。)
g++ -O3 -ffast-math -funroll-loops -Wall -fno-common -march=nocona -o ward -I../include ward.cpp

その結果、以下のように高速化された。
件数   処理時間(s)    そのうち初期行列作成時間
--------------------------------------------
2352 12.628 1.06

シャアではないが、通常の3倍以上の速さである。
C++恐るべし。
分散処理とか考えるのはやはりC++で肝心の部分を書いて、その後なんだろうなと実感しました。

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法をここ数日試してみて、十分速くなったのだが、ここまでくるといろいろ思う事がある。

1. 最終的に今ネックになっているのは、一番最初の行列×行列の計算である。
2. その演算であれば、分散処理での高速化が十分研究されている。
3. 最初は樹形図をつくるところがネックかと思っていたが、そうでないのは驚きである。
4. dict使いまくりの富豪プログラミングはアルゴリズムの試作に便利。

ということで、例えば分散処理で行列×行列を一瞬で計算できるようになれば、直前のエントリのクラスタリングアルゴリズムで1000文書であろうとものの数秒で分類できるということになる。

このクラスタリングの計算方法をC++で書き直せば、もっともっと速くなるはずで、たとえばweb経由でアクセスがあった瞬間に分類アルゴリズムを走らせることもできるだろう。

樹形図をどこでちょんぎるかのさじ加減は別途研究の必要があるが。

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-10

Erlangでベクトルの集合をクラスタリング(k-means)

Erlang本を買ったので、k-means法によるクラスタリングをサンプル的に作ってみた。
Erlang簡潔でいいですね。

使い方
データファイルをタブ区切りで用意する。
1カラム目はデータの識別用ID文字列、2カラム目以降にベクトルの数値データを記入。改行コードは¥nで。

ayu@~/work% /usr/local/erlang_R12B_2/bin/erl
Erlang (BEAM) emulator version 5.6.2 [source] [smp:2] [async-threads:0] [kernel-poll:false]
Eshell V5.6.2 (abort with ^G)
1> c(kmeans).
{ok,kmeans}
2> Src = kmeans:read_line("sampledata2.txt").
[{"1",[1,9,2,4,4,163.5,54,3,1,4,4,2.5,1,1,1,4,3,3,3,2,5]},
{"2",[2,7,1,5,5,154.2,42,2,2,4,3,24.5,1,1,4,3,4,3,5,2,4]},
{"3",[2,8,2,5,3,153.8,44.8,3.5,1,2,5,24,1,1,5,3,2,5,4,3,2]},
.....
3> {L, R} = kmeans:splitvectorset(Src).
{[{"210",[2,4,2,3,3,156,42,3,2,3,4,2,2,2,3,2,3,5,4,1,3]},
{"205",[1,6,1,3,5,170,54,2.5,2,3,3,1.5,1,1,4,3,4,1,2,3,5]},
{"203",[1,2,1,2,3,163,53,2,2,4,4,24,2,1,3,2,2,1,3,1,3]},
....
4>

以上で、LとRに二分割された集合が格納される。

ソース
-module(kmeans).
-import(lists, [map/2, zip/2, sum/1]).
-import(math, [sqrt/1]).
-export([splitvectorset/1, read_line/1]).

%ベクトルのノルム
norm(X)->
sqrt(sum(map(fun(Y)->Y*Y end, X))).

%ベクトル集合の重心
median([],_)->
[];
median(X, L)->
First = lists:nth(1,X),
if
First == [] -> [];
true -> [sum(map(fun erlang:hd/1, X))/L | median(map(fun erlang:tl/1, X), L)]
end.
median(X)->
median(X, length(X)).

%ベクトルの差
sub([], _)->
[];
sub(X, Y) ->
[hd(X)-hd(Y)|sub(tl(X), tl(Y))].

%二つのベクトル間の距離
distance(X, Y) ->
norm(sub(X,Y)).

%与えられた二点X,YのどちらにベクトルZは近いか?
% Xに近ければ'Right', Yに近ければ'Left'を返す。
which(_,_,[])->
ok;
which(X, Y, Z) ->
DX = distance(X,Z),
DY = distance(Y,Z),
if
DX > DY -> 'Left';
true -> 'Right'
end.

%[{'Right/Left', ベクトル}, ...]というリストを'Right'/'Left'に従って二つの
%リストに分けて返す。
splitvector([], R, L)->
{R, L};
splitvector([{D,Z}|T],R,L) ->
case D of
'Right' -> splitvector(T, [Z|R], L);
'Left' -> splitvector(T, R, [Z|L])
end.

%集合ZSの2元からなる部分集合列を返す
pairs([])->
[];
pairs([Z|T]) ->
lists:append(map(fun(X)->{Z,X} end, T), pairs(T)).
%ベクトル集合の距離&ベクトルペアの列を作る
pairdistance([])->
[];
pairdistance(ZS) ->
map(fun({X,Y})->{distance(X,Y), X, Y} end, pairs(ZS)).

%ベクトル集合の最大距離を持つ2点を選ぶ
maxdistance(ZS)->
Pairs = pairdistance(ZS),
{_, MX, MY} = maxdistance(Pairs, hd(Pairs)),
{MX, MY}.
maxdistance([], M)-> M;
maxdistance([Z|T], M)->
{D1,_,_} = M,
{D2,_,_} = Z,
if
D2 > D1 -> maxdistance(T, Z);
true -> maxdistance(T, M)
end.


%二点X,Yのどっちに近いかで、ベクトル集合ZSを二分割する。
splitvectorset(X,Y,ZS, Thre)->
{R,L} = splitvector(zip(map(fun({_,W})->which(X,Y,W) end, ZS), ZS), [], []),
IdRmv = fun ({_, V})->V end,
CR = median(map(IdRmv, R)),
CL = median(map(IdRmv, L)),
MinDist = lists:min([distance(X,CR)+distance(Y, CL), distance(Y, CR)+distance(X, CL)]),
if
MinDist > Thre ->
splitvectorset(CR, CL, ZS, Thre);
true -> {R, L}
end.

%公開関数その1
%ベクトル集合ZSを二分割する。
%集合ZSは、もう一つの公開関数read_lineで読み込んだデータです。
splitvectorset(ZS)->
Threshold = 0.1,
case length(ZS) of
0 -> [[],[]];
1 -> [ZS|[[]]];
2 -> [[lists:nth(1, ZS)], [lists:nth(2, ZS)]];
_ ->
{X, Y} = maxdistance(map(fun ({_, X})->X end, ZS)),
splitvectorset(X, Y, ZS, Threshold)
end.

%分析するデータファイルを読んで{ID, ベクトル}のリストとして返します。
% データファイル仕様
% データファイルはタブ区切り、先頭カラムはレコードID文字列、次カラム以降が数値データ
% 改行コードは\n
read_line(File) ->
{ok, IoDevice} = file:open(File, read),
ANS = read_line(IoDevice, 1, []),
file:close(IoDevice),
ANS.

read_line(IoDevice, LineNumber, Buf) ->
case io:get_line(IoDevice, "") of
eof ->
Buf;
Line ->
Data = string:tokens(lists:delete(10, Line), "\t"), % 10 = "\n"
if
length(Data) == 0 -> ok;
true -> Vec = {hd(Data), map(fun(X)-> parse_num(X) end, tl(Data))},
read_line(IoDevice, LineNumber + 1, lists:append(Buf, [Vec]))
end
end.

%数値を表す文字列をinteger又はfloatに変換します。
parse_num(X)->
case string:str(X, ".") of
0 -> list_to_integer(X);
_ -> list_to_float(X)
end.

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的ベクトルを使用


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


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

coderepos

codereposに自分のページを作った。今後ここに開発中のライブラリのマニュアルなどを書いていこうと思います。

http://coderepos.org/share/wiki/Committers/ayu

2008-05-01

Jython本を注文した

Jython本を注文しました。

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

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

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