実験

並列の実験が終わった．LU分解やNQueenを逐次と並列で高速に解こうといった内容．並列の方はMPI等を使わずにpthreadで．(実験ページ)

LUの方はキャッシュミスをどれだけ減らすか，という感じだった．適当にブロック化してoprofileでプロファイルを取るとキャッシュミスが減ったよ！やったねたえちゃん！みたいな．
SIMDイントリンシック命令を加えるのは実験ページに載っていないのでズルい気もしたが，SSEについては言及されているので気にしないことにした．*1．まぁキャッシュミスを減らす方が性能に影響したんだけど．
ソースコードに溢れる_mm_add_pd等の文字列．可読性が落ちて素晴らしい！気持ちいいね！*2

課題はどちらもググれば色々事例が出てくるので，何だかなぁと言った感じ．特にNQueenの方はググって枝刈りの方法見つけちゃったがな．まぁぼちぼち面白かったけど．

ちなみに並列の実験は先週ぐらいに終了して，今はエージェントの実験に移っている．内容は「SVMでオークションの入札を予測して，最適な入札を考える」といったもの．

二次計画問題を解く部分はライブラリに丸投げ．実験のページではQuadProg++を紹介しているのだが，QuadProg++は半正定値行列の扱いが上手くない(?)ようで，ライブラリがサクサク死んでくれる．サークルの先輩からその話を聞いていたので，対角成分に誤差を加えるという投げ遣りな対処で乗り切った．酷い．実験wikiに解決方法が全く載っていないのが何とも．

最後にどうでもいい話だが，どちらの実験でもオプションの解析を手作業でやっている方をちらほら見かけた．悪いとは言わないが，getopt使った方が楽だと思う……．*3

*1:実際私もその記述からググって辿りついたし

*2:よくねえよ

*3:見えた範囲では突っ込みを入れたが;

11月前半

学園祭に向けてデスマーチしてる(現在進行形)．死ぬ．どうしてこうなった，誰か俺を殺せ．

サークル

NFで自作ゲームの展示とか販売やります．遊びに来てね! リンク

大学

応用代数学の課題を投げ捨てた．同じ講師だった工業数学A3は，最後の講義で「今までのレポートを今からでも受け取るよ」という発言をしていたので，それに賭ける方向で．次からのレポはちゃんと出す

実験

忙しかったこともあり，NQueenの課題は投げ遣りになった．最終レポートの提出がNFの後なので，そこで頑張るかも．

その他

俺妹7巻読んだ．黒猫ちゃんマジ天使

課題とか

初めてのTOEICを受けてきた．最後の方(15問ぐらい)が解ききれず，練習不足な感じ．

課題(括弧の中は締切)

• 応用代数学(11/8)
• 集積システム入門(11/1)
• 工業数学A1(11/4)
• 実験のレポート(11/11, 11/25)
• ヒューマンインターフェース(11/15)

サークル

• 11月祭(学園祭)に出すゲーム作り(内部締切は11/14ぐらい)
• PRML読書会のスライド作り(11/5)

大学の課題より下の2つがまずい．凄くまずい．
そして実験が現実逃避に．一人相撲感やばい

混合ガウスでiris datasetをクラスタリング

Rで学ぶクラスタ解析でやってることを何故かpythonでやってみたよという話．

コードはgistに置いてみた．gistyが何故か入れれなくてしょぼーん．混合ガウスはPRML9章から．
純度とかエントロピー
混合ガウス

irisのデータはここから取得．そのままだとちょっと微妙な気もしたので適当にcut

$ cut --delimiter=, -f1-4 iris.data | sed -e 's/,/ /g' > iris.input
$ cut --delimiter=, -f5 iris.data | cut -c6- > iris.ans

REPLで適当に試す

>>> import scipy as sp
>>> import mydefs
>>> from gmm import gmm
>>> x = sp.loadtxt("data/iris.input")
>>> ans = gmm(x, 3)
>>> cls = ans['classification']
>>> goldans = sp.loadtxt("data/iris.ans", dtype="|S5")
>>> mydefs.myeval(cls, goldans)
Entropy:  0.10167363274
Purity:  0.966666666667

初期値に応じて結果が結構変わる上，たまに潰れて発散する ^^;

3回後期履修科目

全部で18コマ(+α)になった．以下雑感．

数理統計

昨年履修ミスした科目．人多すぎ．教科書は指定されたので，人が減らないようなら出席はしない方向で．

応用代数学

群論とかやるらしい．初回講義の後半に3次方程式の解が云々いい始めたあたりで理解が怪しくなった．おい大丈夫か俺．

集積システム入門

N型とかP型とかの説明から入るあたり凄くゆとり講義でした．講師「情報学科の学生さんは電気回路とか嫌い」．すみません……

計算と論理

1限な上に，資料がpdfで配布されていてこれもまた眠気との闘(ry *1

ソフトウェア工学

初回講義は真面目に聞かずに別のことをしてしまった．というわけで感想は無し．つ，次からちゃんと聞きます．

生命情報学

初回の講義でプロジェクタが使えない→教卓にノートPCを開き，そこでスライドを写しながらプレゼン．立ち見が出るような講義だったので，「その発想は無かった」．生物とかさっぱりだが，出席とレポのみの科目なので頑張る…！

マルチメディア

講師は3人．コンピュータビジョン，自然言語処理，地図・文書画像処理等の色々な内容をやるらしい．あまり深入りせずに表面だけサラっと紹介していく感じなのだろうか．

パターン認識と機械学習

ここ半年ぐらい細々と勉強していた分野なので絶対に負けられない．教科書は「わかりやすいパターン認識」．せっかくなので読んでみる予定．

情報システム

情報検索とかそんな話．IIRは中途半端にしか読んでないので，ちゃんと読みたいなーと思ったりした(講義関係無い

ヒューマンインターフェース

講師の話が面白い．分野としては興味サッパリなのだが，「良いデザインを作れるようになるのではなく，デザイナーさんとコミュニケーションできるようになるのが目的」らしいので，まぁサボらんように頑張る．

アルゴリズム論

2回後期の言語オートマトンの続きの雰囲気．計算可能性とか計算量理論をやるらしい(KULASISからコピペ．正直よく分かってない．)

工業数学A1

昨年度敗北した．専門の再履は実は始めてのハズ．内容としては留数定理のあたりまで．工業数学A3で留数定理が必要になる→自習，というのを前期でやったので，今年は大丈夫だと信じたい．

計算機科学実験及演習(前半，6コマ)

前半は画像処理(QRコードの読み取り)，ロボットプログラミング，並列プログラミングの中から1つを選択．画像を希望したが抽選の結果並列に回された．
並列の内容は{逐次,並列}{LU分解,N-Queen}の高速化．逐次のLU分解は初期状態から3倍ぐらいの速度になったあたりで飽きた．今までキャッシュミス等を意識せずにコードを書いてきたんだなぁと身に染みて感じた．速度ランキングもあるので頑張って上位を狙いたいところ．

*1:せっかく大学生をやっているんだから，家で資料を読むのでは無く講義を聞きに行くのは大切なんだろうけど

zshのalias -sで解凍

zshのalias -sは拡張子に合わせてコマンドを実行できて便利なのだが*1

alias -s 'tar.gz'='tar xzvf'

としても期待通りには動いてくれない(何か間違っているのかな……)
仕方が無いので

function extract() {
    case $1 in
        *.tar.gz|*.tgz) tar xzvf $1 ;;
        *.tar.xz) tar Jxvf $1 ;;
        *.zip) unzip $1 ;;
        *.lzh) lha e $1 ;;
        *.tar.bz2|*.tbz) tar xjvf $1 ;;
        *.tar.Z) tar zxvf $1 ;;
        *.gz) gzip -dc $1 ;;
        *.bz2) bzip2 -dc $1 ;;
        *.Z) uncompress $1 ;;
        *.tar) tar xvf $1 ;;
        *.arj) unarj $1 ;;
    esac
}
alias -s {gz,tgz,zip,lzh,bz2,tbz,Z,tar,arj,xz}=extract

としましたとさ．

これで

$ ./hogefuga.tar.gz

等で解凍できるようになる．
ちなみに他の環境で解凍コマンドを忘れて困りそうだが，今迄もC-rの検索から解凍していたのであまり変わらない気もする．*2

追記(2010/10/21)

コメントでaunpackというのを教えてもらったので，aptitudeからatoolをインストールし，設定は

alias -s {gz,tgz,zip,lzh,bz2,tbz,Z,tar,arj,xz}=aunpack

に変更した．既にあるものを使うのが一番ですね！(おい

*1:参考: http://journal.mycom.co.jp/column/zsh/016/index.html

*2:*.zipと*.tar.gz以外の解凍はまともに覚えていない……

EMアルゴリズムで混合ガウス分布

numpyのufuncの説明を読んで，何かちょろっと書きたい気分になったので書いた．相変わらずの生産性0．

実行すると適当に乱数でデータを生成し，混合ガウス分布のパラメータを推定した後，負担率に応じて色分け．
ufuncのルールをある程度理解したので，3次以上のarrayの扱いがマシになり，ループが減ったものの，弊害として一行に170文字以上書くというPEP8(笑)なコードが出来あがった*1．だからどうしたんだという気もするが．

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

# PRML chapter 9
# Gaussian Mixture Model


import scipy as sp
from scipy.linalg import det, inv


def multivariate_normal_pdf(x, u, sigma):
    D = len(x)
    x, u = sp.asarray(x), sp.asarray(u)
    y = x-u
    return sp.exp(-(sp.dot(y, sp.dot(inv(sigma), y)))/2.0) / (((2*sp.pi)**(D/2.0)) * (det(sigma) ** 0.5))


def gmm(X, K, iter=1000, tol=1e-6):
    """
    Gaussian Mixture Model
    Arguments:
    - `X`: Input data (2D array, [[x11, x12, ..., x1D], ..., [xN1, ... xND]]).
    - `K`: Number of clusters.
    - `iter`: Number of iterations to run.
    - `tol`: Tolerance.
    """
    X = sp.asarray(X)
    N, D = X.shape
    pi = sp.ones(K) * 1.0/K
    mu = sp.rand(K, D)
    sigma = sp.array([sp.eye(D) for i in xrange(K)])
    L = sp.inf

    for i in xrange(iter):
        # E-step
        gamma = sp.apply_along_axis(lambda x: sp.fromiter((pi[k] * multivariate_normal_pdf(x, mu[k], sigma[k]) for k in xrange(K)), dtype=float), 1, X)
        gamma /= sp.sum(gamma, 1)[:, sp.newaxis]

        # M-step
        Nk = sp.sum(gamma, 0)
        mu = sp.sum(X*gamma.T[..., sp.newaxis], 1) / Nk[..., sp.newaxis]
        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]
        pi = Nk / N

        # Likelihood
        Lnew = sp.sum(sp.log2(sp.sum(sp.apply_along_axis(lambda x: sp.fromiter((pi[k] * multivariate_normal_pdf(x, mu[k], sigma[k]) for k in xrange(K)), dtype=float), 1, X), 1)))
        if abs(L-Lnew) < tol: break
        L = Lnew
        print "L=%s" % L

    return dict(pi=pi, mu=mu, sigma=sigma, gamma=gamma)


if __name__ == '__main__':
    data = sp.append(sp.random.multivariate_normal([-3.5, 5.0], sp.eye(2)*4, 50),
                     sp.random.multivariate_normal([-8.2, 10.0], sp.eye(2)*2, 70)).reshape(50+70, 2)
    K = 2
    d = gmm(data, K)
    print "π=%s\nμ=%s\nΣ=%s" % (d['pi'], d['mu'], d['sigma'])
    gamma = d['gamma']

    import matplotlib.pyplot as plt
    plt.scatter(data[:, 0][gamma[:, 0] >= 0.5], data[:, 1][gamma[:, 0] >= 0.5], color='r')
    plt.scatter(data[:, 0][gamma[:, 1] > 0.5 ], data[:, 1][gamma[:, 1] > 0.5 ], color='g')
    plt.show()

"Rで学ぶクラスタ解析"を今読んでおり，それに習ってirisとかのデータを試してエントロピーとか適当に指標とってみようかと思ったけど，やる気が無かったのでまた後日 ;-p

*1:まぁ適当に括弧の中で改行すればいいのだが