2月後半

雑記

サークル

淡々とRWHやPRMLの読書会に参加している．PRMLの方は最後の章が担当になったので，そのうち資料を作らねばならないが，年度内には終わらなそうなのでまだまだまったり．

ちょっと変わったイベントとして皆でサバゲに行った．当日朝までルールを何も知らないド素人だったが，スタッフの方々が丁寧に解説してくれたおかげで楽しく参加出来た．また友人に誘われたら行ってみたい気もする．

毎年恒例の春合宿が近いので講座資料を考え中．機械学習の話でもするかなぁ，という気分．出来るだけ難しい説明を省いて，1回生などに楽しんでもらえる内容にできればと考え中．

自炊

電子書籍リーダーを買った．KindleとSony Readerで悩んだが結局Reader(PRS-650)にした．今の所自炊中心なのでkindleだとSDカード等が刺せずちょっとツラい．適当な場所に置いたファイルを見たくなったら3Gで通信して取ってくれば良い気もしたが，いつまで3Gが無料か分からないのでReaderにしてみた．

ScanSnapで読み取ったラノベはそのままで読めるが，ChainLPで変換した方が見易い．余白をカットすれば大体元の書籍と同じぐらいのサイズで文字が表示される．E Inkは読み易くて素晴らしいが，冷静に考えると一日中ディスプレイを眺めている日もあるので，液晶でも多分困らなかった気もする;-p

ラノベ以外では途中まで読みかけていたIntroduction to Information Retrievalを入れた．そのままhtmlをwget -rしてebook-convertでepubにするものの，ファイルが大きすぎるのか開くと途中でReaderが死ぬ．じゃぁ試しにLRFにでもしてみるか，と思ったらebook-convertが死ぬ．epubをsectionごとに適当に分割してくれるものも見つからなかったので，仕方無く自分でセクション毎にebook-convertした(ついでに邪魔なcopywriteやnavigationも消したり)．この手の投げ遣りなコード書くときはついbashとかsedとか使って書いてるけど何だか時代錯誤っぽい気もする．

大学

何も無い．他の学科の友人は研究室見学をしているようだが，うちの情報学科は4月頭に行うらしい．あいかわらずの適当さ．狙っている研究室は今のところ知能情報のY研かな，という気分．

その他

AmazonのAPIで返ってきたXMLをPythonのElementTreeのfindで検索したら上手く行かずちょっとだけハマった．名前空間の指定は find("//{namespace}hogefuga")とやれば良いらしい．

言語処理のための機械学習入門を読んだ．全体を通して例題が多く，理解の助けになって良いと思う．具体的な例があるとやっぱり楽．章末の問題にも解答が付いており独習しやすい．ただCRF以外はある程度既知だったので，いや何というか読む順番を間違えた気がするというか，PRMLを読み始めた後輩(1回生)に今度勧めてみるかなと思った．

自炊した本のpdfファイルからバーコード画像を読み取りISBNを取得してAmazonのAPIから得たタイトルのリンクを張りたい

Python sh

ScanSnapと裁断器で延々と刺身にタンポポを載せるような仕事をしていたのですが，タイトルを入れるのが面倒なので自動で頑張ることにしました．ScanSnapを使ったものの以下の作業はLinux上で行っています．

前提

自炊した本にISBNのバーコード画像が入っていること．裏表紙を捨てずにスキャンしていれば多分OK．ファイルはpdfを仮定していますが，そうで無い場合(zipに固めてみました等)は適当に頑張って下さい．

方針

自炊した本のファイル名は".pdf"に変換．日本語のファイル名は別途リンクを貼ります．直接ファイル名を変更しないのは，

  • あとでファイル名の形式を変えたくなるかもしれない．その際もう一度ISBNの番号を取得するのはバカっぽい．
  • utf8だけで無く，sjisのタイトルのリンクも作っておけばWinSCPでWindowsのノートPC等に転送する場合も楽……な気がしたので

準備

ZBar bar code reader
zbarimgで画像からバーコードのISBNを取得します．squeezeやUbuntuならaptitudeで入るようですが*1，我が家のDebianは未だLennyなのでソースから入れました．configureは --without-qt --without-gtk で良いかと．ちなみにconfigureでMagicWandがうんたらと怒られたのでMagicWandもインストール．
pdfimages, pdfinfo
Debianならxpdf-utilsをパッケージ管理システムで入れれば良いでしょう．
python
2.5以上3未満．

注意

ISBN13の先頭が978であると勝手に仮定しています．

ファイル名をISBNに変換

pdfの最初と最後の3ページのどれかにバーコード画像が入っていると仮定しています．pdfが置かれたディレクトリ内で

% ./isbnfilename.sh

でOK．ソースはgistに．bashに慣れちゃってshが怪しいゆとりですよ，ええ……．

Amazon APIでISBNからタイトル等を取得

AccessKeyIDとSecretAccessKeyとやらが必要らしいので，Amazon Web Serviesのページから登録．登録後は アカウント->セキュリティ証明書 から "アクセスキー ID" と "シークレットアクセスキー" とやらを取得しておきます．

ISBNからタイトルを得るコードはgistに．pdfがあるディレクトリに置いて

% ./isbn2title.py <AccessKeyID> <SecretAccessKey> <ISBN>

とすればタイトルが取得できます．以下は例．

% ./isbn2title.py ひみつ ひみつ 9784047270305 
[エンターブレイン][庵田 定夏][白身魚]ココロコネクト ミチランダム (ファミ通文庫).pdf

"[出版社][著者][イラスト]タイトル.pdf"な形式です．何というラノベ専用形式．気に入らなければ適当に変更して下さい．

日本語タイトルのリンクを作成

% for file in *.pdf; do ln $file "$(./isbn2title.py ひみつ ひみつ $(echo $file | sed -e 's/\.pdf//'))"; done

感想とか

毎回コマンドを叩くのが面倒．次は指定ディレクトリにpdfを放り込んだら勝手に変換するようにすると幸せなのかな．

*1:zbar-tools

「応用のための 確率論入門」を読んだ

トピックモデル楽しそうだなーと思って資料を漁っていたらディリクレ過程の説明で"測度"とか出てきて脳が拒否反応したでござる．「情報の計算機科学コースじゃ真面目な確率なんていらないだろwww」と思って測度は避けていたが，せっかく春休みで時間もあるので勉強しようと思って本を読んでみた．

集合や測度の話は1から丁寧に解説してくれるので，前提知識は殆どいらない．大学1回生とかでも余裕で読める．もう少し早く読んでおけば良かったなと思う．

2章の測度の説明は面積を定義する話から始まり，完全加法性を満たしながら面積を持つ部分集合を広げようという流れで話が展開されていく．具体的な例と丁寧な説明のおかげですんなりと理解できた．

分かりやすくて素晴らしい本だったが，証明はいくらか省かれており，測度やルベーグ積分については必要最低限だけを紹介しているように感じた．数学科の人は怒り出しそうだが，私みたいにてっとり早く測度論的確率論を知りたい人には丁度良いと思う．

測度論的確率論はこの本でしか学んでいないので，他の書籍と比較してどうなのかは分からないが，数学科出身のサークルの大先輩が情報ならこれで十分と仰っていたので，気にしないことにする．


しかし試験前日だとこのぐらいの本なら1日で読み終わるのに，春休みだと2，3日かかってしまうのは何なんだ．あれか，追い込まれないとダメなのか．

3回後期成績発表

大学

1科目(情報システム)未採点で+23単位．108->131．戦闘力*1(暫定)は 7780 ．

優

計算機科学実験及演習4
作業．出席が一番面倒．
工業数学A1
5問中レポートから1題，講義内容から1題だった．留数定理まで理解して，レポートをやっていれば取れそうな内容．
アルゴリズム論
今年はレポートの比重が大きかった．試験は問題が非可解であることを証明したりする感じ．
マルチメディア
大問3つのうち1つが「人間の眼の働きについて述べよ」「カメラと比較しながら説明せよ」という予想外の問題だった*2．
ヒューマンインターフェース
レポート(2回)が面倒．試験は毎年似た傾向なので簡単．
生命情報学
出席とレポートのみ．楽勝の部類．
パターン認識と機械学習
ありがとうPRML！

良

ソフトウェア工学
必修なので落とすと面倒．
数理統計
一夜漬け
集積システム入門
試験は持ち込みアリだが簡単だった(毎回課題をちゃんと出していれば)．

可

計算と論理
CAL100問で可．途中のHPropを飛ばしてさっさとLambdaTerm(or Arith)に行くと楽．

不受験

応用代数学
試験期間に過去のレポートの後出し提出を許可していた上，試験の採点も緩かった模様．受ければ良かったな．

*1:B群+専門の単位数に 可:60 良:70 優:80 をかけたもの

*2:そこは出ないと思ってたぞ……．適当に書けば点をもらえるんだろうけど．

scipy/numpyのufunc備忘録

Python

scipy/numpyの行列演算は詳しいことを知らなくても何となく使えて便利なのですが，慣れていないとついついループや内包表記を使ってしまいます．pythonのループで書くと遅くなるのは分かっているのに……．最近「ufuncのbroadcastingが分かっていれば大分マシになるのでは無いか」と考えるようになったので，備忘録も兼ねて適当なことを書いてみます．


ちなみに私の数値計算の知識はサッパリなので注意．またnumpy/scipyに関して何も知らない場合は他の記事を当たった方が良いと思います．
計算の速度やメモリのことは良く分からないので，難しいことを考えるのを止めて，目標を ループを使わずにとにかく1行で処理を書く にします．おいおい．

broadcasting

universal function(以下ufunc)はarray*1に作用する関数のことで，これの"array broadcasting"のおかげで行列とスカラーの足し算等が簡単に出来たりします． 詳しくは公式の解説を読め！で終わる話なのですが；

broadcastingのルールは以下の4つ．

  1. 各入力のarrayに対してndimが最大の入力より小さければ，1がshapeの先頭に追加される．
  2. 出力の各次元の大きさは，入力の対応する次元の最大値となる．
  3. 入力の各次元の大きさは，対応する出力の次元の大きさと等しいか1で無ければならない．
  4. 入力のある次元の大きさが1の時，その要素がその次元の計算で繰り返し利用される．

何のこっちゃさっぱり分かりません．ということで簡単な例で試してみます．
ちなみにndimは次元数，shapeは各次元の大きさのタプルです．array([ [1,2,3], [4,5,6] ]) の ndimは2，shapeは(2, 3) といった感じ．

例1

行列とベクトルの足し算を考えます．

import scipy as sp
a = sp.array([[1,2,3],[4,5,6],[7,8,9]])
b = sp.array([10,20,30])
print(a+b)

[[1, 2, 3],
 [4, 5, 6],  + [10, 20, 30]
 [7, 8, 9]]

a,bのndimが一致していないので，ルール1よりbのshapeの先頭に1が追加される(b.shape -> (1,3))

[[1, 2, 3],                  
 [4, 5, 6],  + [[10, 20, 30]]
 [7, 8, 9]]

a,bのshapeはそれぞれ (3,3), (1,3) なのでルール2より出力は(3,3)となる

[[1, 2, 3],                     [[?, ?, ?],
 [4, 5, 6],  + [[10, 20, 30]] =  [?, ?, ?],
 [7, 8, 9]]                      [?, ?, ?]]

ルール3は満たしているのでエラーは起きない．ルール4よりbの[10,20,30]が繰り返し使われる

[[1, 2, 3],    [[10, 20, 30],   [[?, ?, ?],
 [4, 5, 6],  +  [10, 20, 30], =  [?, ?, ?],
 [7, 8, 9]]     [10, 20, 30]]    [?, ?, ?]]

これより結果は

[[11, 22, 33], 
 [14, 15, 16], 
 [17, 28, 39]]

となります．

例2

a = [a0, a1, ..., aM] と b = [b0, b1, ..., bN] から

[[a0*b0  a0*b1 ... a0*bN ],
 [a1*b0  a1*b1 ... a1*bN ],
 ...
 [aM*b0            aM*bN ]]

を求めたい．共分散行列とかそんな時に．ちなみにnumpy.outer使えば終わる話なのですが;*2
aとbのshapeはそれぞれ(M,)と(N,)なので結果の(M,N)を求めるには，aのshapeを(M,1)にすれば上手く行きそうです．このような操作はnumpy.newaxisを利用してやればOK*3．具体的には

In [49]: import scipy as sp
In [50]: a = sp.array([1,2,3,4,5])
In [51]: b = sp.array([2,2,2,3,3])
In [52]: a[:, sp.newaxis]
Out[52]: 
array([[1],
       [2],
       [3],
       [4],
       [5]])
In [53]: a[:, sp.newaxis]*b
Out[53]: 
array([[ 2,  2,  2,  3,  3],
       [ 4,  4,  4,  6,  6],
       [ 6,  6,  6,  9,  9],
       [ 8,  8,  8, 12, 12],
       [10, 10, 10, 15, 15]])

具体的な添字の代わりに":"を指定すると「その次元の要素を全て」と認識しているが，実際どうなのかはよく知りません……．
ちなみに":"の個数が分からない場合は(入力の次元数が変化する)，"..."とかを利用すると良い感じ．(例)

例3

EMアルゴリズムで混合ガウス分布のパラメータを求める．PRMLの9章に載ってるやつですね．
以前書いた記事の使い回しですが……．

記事書くの疲れてきたので数式は省略します．下巻p154とか見ればいいんじゃないでしょうか．各変数のshapeは

  • gamma: (N,K)
  • Nk: (K,)
  • mu: (K,M)
  • sigma: (K,M,M)
  • pi: (K,)
  • X: (N,M)

とします．Mステップの更新だと
Nkはgamma(N,K)から(K,)を足し算で生成すれば良さそうです．ということで

Nk = sp.sum(gamma, 0)

sumの第二引数に数字を指定すると，その次元で足し算します．

sigmaはgamma(N,K)とX(N,M)とmu(K,M)とNk(K,)から(K,M,M)を作れば良さそうです．
まずは(N,K,M,M)を作り，そのあとsumでΣ相当のことをすると考えて

xmu = X[:, sp.newaxis, :] - mu
sigma = sp.sum(gamma[..., sp.newaxis, sp.newaxis] * xmu[:, :, sp.newaxis, :] * xmu[:, :, :, sp.newaxis], 0) / Nk[..., sp.newaxis, sp.newaxis]

PEP8(笑) *4

例4

以前しましまさんにコメントで尋ねられた例．

ところで numpy/scipy とかpythonの数値計算系はお詳しいですか？
私は numpy 修行中の身なのですが，演算の broadcasting とかがよく分からず，for文に手をだしてしまいます．

for x in xrange(m.nx):
for y in xrange(m.ny):
q[:, x, y] = m.pz * m.px[x, :] * m.py[y, :]

とかを1行で書こうとして挫折しました．scipy のサイトや，Python Scripting for Computational Science 本以外に良い資料をご存じでしたら教えていただけませんか？

http://d.hatena.ne.jp/se-kichi/20100519/1274275285#c1274496705

nxやnyやpxが何かよく分からないので「m.nxやm.px等はscipy/numpyとは関係無い」として考えます．
コメントのコードを見ると

  • qの次元数は3．
  • m.px と m.py の次元数は2．*5

と考えられる．m.nx, m.ny(それぞれX, Yと置く)でループをしているので，この値が添字で使われている次元の大きさに等しいと仮定すれば

  • q.shape -> (K, X, Y)
  • m.px.shape -> (X, L)
  • m.py.shape -> (Y, M)

となる．

  • q[:,x,y].shape -> (K,)
  • m.px[x, :].shape -> (L,)
  • m.py[y, :].shape -> (M,)

であるから K=L=M と分かる．m.pzは()*6と(K,)の可能性があるが，前者ならばループの外側で乗算を行うだろうから，後者だと考えて

  • q.shape -> (K, X, Y)
  • m.px.shape -> (X, K)
  • m.py.shape -> (Y, K)
  • m.pz.shape -> (K,)

となる．結局 m.pz, m.px, m.py から (K, X, Y) を乗算で作るので

  • m.px.shape -> (X, K) -> (K, X) -> (K, X, 1)
  • m.py.shape -> (Y, K) -> (K, Y) -> (K, 1, Y)
  • m.pz.shape -> (K,) -> (K, 1, 1)

と変換して計算すれば良さそう(というかこれ以外の選択肢無い……)．コードにすると

q = m.pz[:, np.newaxis, np.newaxis]*m.px.T[:, :, np.newaxix]*m.py.T[:, np.newaxix, :]

実際に対話環境で試してみると

In [84]: import scipy as sp
In [85]: K = 10
In [86]: X = 13
In [87]: Y = 19
In [88]: px = sp.random.randint(0, 100, (X, K))
In [89]: py = sp.random.randint(0, 100, (Y, K))
In [90]: pz = sp.random.randint(0, 100, (K,))
In [91]: q1 = pz[:, sp.newaxis, sp.newaxis]*px.T[:, :, sp.newaxis]*py.T[:, sp.newaxis, :]
In [97]: q2 = sp.zeros((K, X, Y))
In [98]: for x in xrange(X):
             for y in xrange(Y):
                 q2[:, x, y] = pz * px[x, :] * py[y, :]
In [105]: sp.all(q2 == q1)
Out[105]: True

*1:正確にはndarray

*2:ちなみにufuncにはouterメソッドがあるので，足し算で似たようなことがやりたければ numpy.add.outer(a, b)

*3:ちなみに newaxis は None

*4:Limit all lines to a maximum of 79 characters.

*5:m.px[x, :]から分かるのは「m.pxの次元数 >= 2」だが，q[:, x, y]の次元数は1なのでm.pxの次元数は2．

*6:スカラー

2月前半

雑記

はてダの「記事を書く」をクリックしたら "2010年を振り返って" "2011/01/01 10:02 にバックアップされた編集内容を挿入しました" と出てきて死にたくなった．しかもほぼ白紙じゃないですかー！

サークル

何故かフリークライミングに行った．日頃引きこもっているせいか，体を動かすとサッパリする．次にサークル関連で運動するのは3月のサバゲ*1．

3月には恒例の春合宿があるので，講座で何を喋るか考え中．クラスタリングとか文書分類とか何か喋ろうかなーと思いつつ予定は未定．

大学

試験は大体取れた(ハズ)．他学部の方々は既に研究室配属が決まっているようだが，京大情報学科は4月に研究室見学+配属決定．しかしどこにしようか……．*2

*1:一体何のサークルなんでしょうね

*2:2つぐらいまでは絞った気分

学園祭終わってました

雑記

シフト以外は大抵寝てた．KUSFAと漫トロと朽木の冊子は回収したけど，他は全く覗いてない．酷い．


サークルの方ではC#とBox2D.XNAでパズルのようなアクションのようなゲームを作ったりした．多分オリジナル．あとはwiiリモコンでガンシューティングっぽいのとかも作ったり．


締切前は「そもそもこれはゲームとして面白いのか」「簡単すぎる．難易度もうちょっと上げるか……と思ったけど客は大丈夫なのか」「死ぬ．講義なんて無かったレポートは捨てよう今日の講義は止めよう」と苦しんでいるのに，終わってみると楽しい思い出に感じる不思議．来年もやれる範囲で頑張りたい．研究室配属とか卒論とかで大変そうだが．


来年と言えば今回はスタートが異様に遅かったので，次はもっとちゃんと準備したいところ．*1

*1:3週間前にのろのろ書き始めたのは酷い．しかもゲームの内容決まってなかったし……