jythonがmacbookにインストールできなかったのでとりあえずlinuxで。

といってもgcjだとエラーを吐くので、sunのjavaをインストール。 ここみて設定。alternativeコマンドを使うとjavaの共存ができるのね。いままで、シンボリックリンクを上書きしてたのでめんどいなーと思ってたけど、これだとらくちん。

jythonをインストールしたら、あとはopsinのjarを落としてきてクラスパスに通す。

>>> import uk.ac.cam.ch.wwmm.opsin as opsin
>>> opsin.NameToStructure().parseToCML("4-iodobenzoic acid").toXML() 

二行でIUPAC名がCMLに。 すばらしい。

jrubyの例もある。

ちなみにjrubyはmacbookにさくっと入って、この例の通りにやれば動いた。

perlで同じ事をやる場合にはInline::Javaを使ってやればいいけど、Javaのライブラリを有効に利用するのはJavaで実装された言語処理系がやっぱ楽だ。

JRubyとかjythonとかはCGI書くのに苦労するよなとかつぶやいてみたけどjythonではSimpleHTTPSeverが使えることに気づいた。

% jython
Jython 2.2.1 on java1.6.0_06
Type "copyright", "credits" or "license" for more information.
>>> import SimpleHTTPServer
>>> SimpleHTTPServer.test()
Serving HTTP on 0:0:0:0:0:0:0:0 port 8000 ...
192.168.11.xx - - [26/Jun/2008 21:34:28] "GET / HTTP/1.1" 200 -
192.168.11.xx - - [26/Jun/2008 21:34:28] code 404, message File not found
192.168.11.xx - - [26/Jun/2008 21:34:28] "GET /favicon.ico HTTP/1.1" 404 -

あーでも、これだと、結局jythonスクリプトでCGI実行しないとあかんからうまくいかんわ。

というわけで、明日は会社にWebアプリ編を持ってってjythonでwebappサーバーつくるのにチャレンジ。

みんなのPython Webアプリ編 [みんなのシリーズ]
柴田 淳
ソフトバンククリエイティブ / ￥ 2,604 ()
在庫あり。

なんかやる気出てきた。

Programming Collective Intelligenceを読み始めた。

Programming Collective Intelligence: Building Smart Web 2.0 Applications
Toby Segaran
Oreilly & Associates Inc / 3446円 ( 2007-08 )

英語だとスペースで分割すれば単語に分けられるのだけど、日本語は品詞分解できないと類似度を計れないしクラスタリングもできないので、まずMeCabを使えるようにしておく。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys,re,feedparser
import MeCab

d = feedparser.parse('http://blog.kzfmix.com/rss/')
txt = ''
for entry in d.entries:
    txt += re.compile(r'<[^>]+>').sub('',entry.summary_detail.value)

try:
    t = MeCab.Tagger()
    m = t.parseToNode(txt.encode('utf-8'))
    while m:
        if m.stat < 2:
            if re.match('名詞',m.feature): print m.surface
        m = m.next
except RuntimeError, e:
    print "RuntimeError:", e;

1章はcollective intelligenceの紹介。

2章は類似性を評価するための距離の計りかた。ユークリッド距離とピアソン距離。あとエクササイズにタニモト距離。

Programming Collective Intelligence: Building Smart Web 2.0 Applications
Toby Segaran
Oreilly & Associates Inc / 3446円 ( 2007-08 )

無難に読み流す。

LDRのOPMLからRSSのURLを抜き出す

from xml.dom.minidom import parse, parseString
urls = []
dom = parse("export.xml")
for outline in dom.getElementsByTagName("outline"):
    url = outline.getAttribute("xmlUrl")
    if url: urls.append(url)

楽ちん

LDRで購読しているフィードから単語セットを抽出して遊びたい。 データセットは一度取っておけばいいので、永続化をしておく。入力はLDRから吐き出したOPMLファイル(export.xml)

集合知プログラミング
Toby Segaran
オライリージャパン / ￥ 3,570 ()
在庫あり。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys,re,feedparser,shelve,MeCab,sgmllib
from xml.dom.minidom import parse, parseString

opmlfile = "/Users/kzfm/export.xml"

def getwordcounts(i,url):
    print "#%d feedparser parse: %s" % (i,url)
    try:
        d = feedparser.parse(url)
    except UnicodeDecodeError,e:
        print "UnicodeDecodeError ", e
        return {}
    except sgmllib.SGMLParseError,e:
        print "sgmllib.SGMLParseError", e
        return {}

    txt = ''
    tage = re.compile(r'<[^>]+>');
    for entry in d.entries:
            if hasattr(entry,"summary_detail") : txt += tage.sub('',entry.summary_detail.value)

    try:
        t = MeCab.Tagger()
        m = t.parseToNode(txt.encode('utf-8'))
        wc = {}
        while m:
            if m.stat < 2:
                if re.match('名詞',m.feature): wc[m.surface] = wc.get(m.surface,0)+1
            m = m.next
        return  wc
    except RuntimeError, e:
        print "RuntimeError:", e;

def parse_opml(file):
    urls = []
    dom = parse(file)
    for outline in dom.getElementsByTagName("outline"):
        url = outline.getAttribute("xmlUrl")
        if url: urls.append(url)
    return urls;

if __name__ == "__main__":
    url_lists = parse_opml(opmlfile)
    tag_data = shelve.open("myldrwordsets")

    for i,url in enumerate(url_lists):
        tag_data[url] = getwordcounts(i,url)

    tag_data.close()

• 1000件ぐらい取ってくるのに30分ぐらいかかったが、ダウンしているサイトをずっと待っているのはよくない。feedparserのタイムアウトの設定ってどうやんのかな。
• 今回はshelveを使ってみたがpickleとの使い分けの基準がいまいちわからん

できたファイルのサイズは15Mくらいだった。これで類似度を計算したり、クラスタリングするためのデータは揃った。

追記 08.07.13

shelveだとなぜか辞書の復元がうまくいかなかったのでcPickleに変更した

if __name__ == "__main__":
    url_lists = parse_opml(opmlfile)
    tag_data = {}

    for i,url in enumerate(url_lists):
        tag_data[url] = getwordcounts(i,url)

    f=file("myldrwordsets","wb")
    pickle.dump(tag_data,f)
    f.close()

RSSから単語セットを抽出したら続いて、頻度のマトリックスを作成する。

集合知プログラミング
Toby Segaran
オライリージャパン / ￥ 3,570 ()
在庫あり。

流れとしては、全単語セットからある頻度以上の単語リストを作成（これが行のセット）し、それぞれのURL毎に何個含まれているかを調べてマトリックスにする。この時、あまりにもマイナーな単語とか、逆にあまりにもよく出る単語は類似性をはかるときに情報量ゼロにしかならないので除く。 今回1100件のフィードからなる（かなり多様性の高い）セットだったので、minの頻度はかなり低く設定した。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import cPickle,sys

def get_tagset(tag_data,min=0.1,max=0.7):
    all_tags = []
    tag_counts = {}
    size = float(len(tag_data))

    for tags in tag_data.values():
        for tag in tags:
            tag_counts[tag] = tag_counts.get(tag,0)+1

    for (tag,count) in tag_counts.items():
        if min < count/size < max: all_tags.append(tag)
    return all_tags

if __name__ == "__main__":
    out=file('blogdata.txt','w')
    f=file('myldrwordsets','rb')
    tag_data = cPickle.load(f)
    wordlist = get_tagset(tag_data,min=0.1,max=0.5)
    totalwords = len(wordlist)
    out.write('Blog')
    for word in wordlist: out.write('\t%s' % word)
    out.write('\n')

    for url,wordcount in tag_data.items():
        out.write(url)
        count = 0
        txt = ""
        for word in wordlist:
            c = wordcount.get(word,0)
            txt += '\t%s' % c
            if (c != 0): count += 1
        if float(count)/totalwords > 0.000000001:
            out.write(txt + '\n')
            print "%s: %2.8f" % (url,float(count)/totalwords)

実際に作成してみたデータを眺めた。というより、ケモインフォマティクス系はタニモト距離をよく使う関係上、密度(density)には非常に敏感。というわけで、url毎に非0の単語がどのくらいの割合含まれているかをチェック。やたらとビット密度の低いフィードは距離を算出しても仕方ないので除く事にした。

        if float(count)/totalwords > 0.000000001:
            out.write(txt + '\n')
            print "%s: %2.8f" % (url,float(count)/totalwords)

結局、915x300くらいのマトリックスが作成された。

類似度はかるのに尺度はなにを使おうか、、、とりあえず本にある通りにピアソン使おうかな。それとも今回の場合は単語の出現頻度よりも、ある文書に単語リスト中の文字が含まれるかどうかのほうが重要そうだからタニモトのほうがよさそうな気もする。

pythonで01の文字列を10進数に変換するには

int('101',2) # 5

とできるので、あとは、C=(A&B)としてC/(A+B-C)という式で1のビットを数えあげればタニモト係数は求まる。

ビットの数え上げはシフト演算子を使ってやればよさげ。

>>> def cb(x):
...   count = 0
...   while(x!=0):
...     count += x & 1
...     x >>= 1
...   return count

こんな感じでOKと思うんだけど大きい数字を与えた場合に11LとかLがついてくるのが気持ち悪い。

タニモト係数はさっきの関数を利用して

>>> a = int('10111',2)
>>> b = int('11010',2)
>>> float(cb(a&b))/(cb(a)+cb(b)-cb(a&b))
0.40000000000000002

さくっと求まりそうだ。

portで。gccのコンパイルに結構時間がかかった(3時間くらい?)。途中で放って寝たので知らんけど。

sudo port install py-f2py
sudo port install py25-numpy
sudo port install py25-scipy
sudo port install py25-matplotlib +tkinter
sudo port install py25-ipython

あと、scipy,matplotlibのインストールが途中でこけるのだけど、再度同じコマンドを打てばいいらしい。

元ネタはRadar / Spider Chars

五角形にしたかったのでrule of fiveにPolar Surface Areaを加えておいた。

#!/usr/bin/env python
from matplotlib.projections.polar import PolarAxes
from matplotlib.projections import register_projection
from pylab import *

class RadarAxes(PolarAxes):
    """Class for creating a radar chart (a.k.a. a spider or star chart)        
    http://en.wikipedia.org/wiki/Radar_chart
    """
    name = 'radar'
    # use 1 line segment to connect specified points
    RESOLUTION = 1

    def draw_frame(self, x0, y0, r):
        verts = [(r*cos(t) + x0, r*sin(t) + y0) for t in theta]
        return Polygon(verts, closed=True)

    def set_varlabels(self, labels):
        self.set_thetagrids(theta * 180/pi, labels)

    def get_axes_patch(self):
        x0, y0 = (0.5, 0.5)
        r = 0.5
        return self.draw_frame(x0, y0, r)

if __name__ == '__main__':

    register_projection(RadarAxes)
    N = 5

    theta = 2*pi * linspace(0, 1, N+1)[:-1]
    theta += pi/2
    labels = ['HBA', 'HBD', 'cLogP', 'MWT', 'PSA']
    rule_of_five = [10, 5, 5, 500, 140]
    desc = [12, 3, 3.6, 532, 160]
    desc_rate = [100*desc[i]/float(v) for (i,v) in enumerate(rule_of_five)]

    ax = subplot(111, projection='radar')

    ax.fill(theta, [100]*N)
    ax.fill(theta, desc_rate)

    for patch in ax.patches:
        patch.set_alpha(0.5)

    ax.set_varlabels(labels)
    rgrids((20, 40, 60, 80, 100))

    grid(True)
    show()

rule of fiveのようにある範囲内に収まっていること(超えるとリスク)というような指標を表すのにレーダーチャートは適しているんだろうか。つまり充足している事を示すような面積の表現はいいのかなぁ。あと、レンジが負になったりするのでそれもどうかと思う。

再考の余地はあるな。

TODO:多変量がある決まったレンジ内に収まっているかどうかを視覚的に捉えやすい表現手段を探す。