2015-01-05

2015年の目標

2014年は、特にTopCoderのMarathon MatchでTシャツに手が届いたのでめちゃくちゃ嬉しかったですが、目標が具体的じゃなかったので、あまり振り返ることが出来なかったような気がします。そこで、今年はTL上の方々を参考にして得点形式にすることにしました。

必要に迫られているので、ややC++の学習を重くしています。正直必要になるかはわからない感じですが、
ソースを読まないで勉強できるテクニック等は今のうちにひと通りやってみたいなぁという感じです。

できるだけモチベが保たれるように、できるだけ連続的に得点が動くようにはしてみました。

2013-12-31

2013年を振り返って

2013年は比較的色んなことに手を出せた、いい年だったと思います。特に2013年のはじめに目標を立ててなかったので、今年手をだせたことを振り返ってみます。

マラソンマッチ

普段コンテストには出てなかったものの、twitterのTL見てて面白そうだなぁと思ったのをきっかけに、マラソンに出ることにしました。
https://twitter.com/xyz600600/status/332127635557720064:twitter:detail:left
実際やってみると、Viewerとにらめっこしてあれこれ工夫を考えるのが凄く楽しくて、どっぷりはまってしましました。今年はTestMMも含めて４回出場して、1766までレートが上がりました。もう少しレートは上げられる感じがありますが、1900とかいくと現状どうにもならなそうな壁があるのは間違いないっぽいので、2014年はとにかくコンテストの数をこなして経験を積めたらいいなぁと思います。あと、TCOのTシャツがどうしても欲しいので、2014年は狙いに行きます!

レンダリング（パストレーシング）

始めたきっかけは、夏休みの終わりあたりに、同じ学校の友人が実装したいと言い出したことだったと思います。資料は、@h013さんが書いた資料を参考にしました（tweet内のURL）。正直、これがなかったらレンダリングに興味持ってた可能性は低かったと思います。
https://twitter.com/xyz600600/status/376557974489030656:twitter:detail:left
2013年はあれこれ論文を読みながら、BVHの実装、基礎的なパストレの実装など基本的な実装をしてました。2014年は双方向パストレを筆頭にもう少しリッチなアルゴリズムに手を出したいなぁと思います。あと、レンダリングをきっかけに準モンテカルロ法の理論に興味を持ったので、時間に余裕があれば手を出したいです。

プログラミング言語

毎年１つのプログラミング言語に手を出してみようと、修士入ってから考えていました。2013年はScalaのつもりだったんですが、見事に使う機会がなくて結局定着しませんでした… 少し前にD言語を初めてみたので、2014年はD言語を頑張ってみようと思います。あと、マラソンの上位陣がみんなC++なのを受けて、C++を少しずつ使ってみるつもりです。C++11が思ったより使いやすそうだったので、Javaから完全に離れる日が来るのかもしれません。

ジョギング

今年は何故かマラソン大会(5km〜20kmまで）に割と数多く出場しました。もう部活動みたいにガチでやることもないと思いますが、これからも定期的に大会に出て楽しく走りたいです。出場した大会は4つだったと思います。

どうでもいいですが、趣味のマラソンが２種類あるせいで割と適当に話してた時に誤解させた覚えが何度かあったので、2014年は頑張ってどちらか一意に特定できる話し方にします…w

2012-12-02

SVMを実装してみた

python

授業でSVMについて習ったけど、実際に実装したことなかったからやってみた。簡単って言われてるけど、制約付き2次計画問題の実装が結構大変だった（収束しないケースとかたくさんあったり、制約条件を遵守したり）

参考にしたのは、以下の本やらページやら
http://www.amazon.co.jp/%E3%82%B5%E3%83%9D%E3%83%BC%E3%83%88%E3%83%99%E3%82%AF%E3%82%BF%E3%83%BC%E3%83%9E%E3%82%B7%E3%83%B3%E5%85%A5%E9%96%80-%E3%83%8D%E3%83%AD-%E3%82%AF%E3%83%AA%E3%82%B9%E3%83%86%E3%82%A3%E3%82%A2%E3%83%8B%E3%83%BC%E3%83%8B/dp/4320121341/ref=sr_1_1?ie=UTF8&qid=1354439820&sr=8-1
http://home.hiroshima-u.ac.jp/tkurita/lecture/svm/index.html

制約条件の中でも厄介なのが、ラベルのベクトルyと重み係数ベクトルαの内積が0になる必要があるという制約で、普通にやると結構この制約を満たすだけでも一苦労…
そこで考えたのが、最初からyに直交する空間上で最適化するという案。具体的に基底を定めてやると、意外と簡単にαからβという別の重み係数の最適化問題に落とせたので、それで実装した。ちなみに最適化エンジンは単なる勾配法。
完成はしたけど学習率とかカーネル等を上手く調整しないと収束しない場合が結構ある。理由はよく分かってないけど、経験的にはサポートベクターの数が多いと収束しやすい傾向があると予想してる（性質としては良くないけどw）

コードはこんな感じになった。サポートベクター表示させてなかったり、サポートベクターじゃない奴も含めてクラス分類させてたりあるけど…
行列による一斉更新ってところは、βを逐次的に変えるより速いんだけど収束しなかったり収束が遅かったりするからコメントアウトしてある。

# -*- coding: utf-8 -*-

from pylab import *
from numpy.linalg import *

LABEL = 0
VECTOR = 1

# y = func(x)が境界線
def make_input_data(Nsample, func):
    data = []

    #正例のデータ作成
    while len(data) < Nsample / 2:
        xs = rand(2)
        if xs[1] >= func(xs[0]):
            data.append((1, xs))
            
    #負例のデータ作成
    while len(data) < Nsample:
        xs = rand(2)
        if xs[1] < func(xs[0]):
            data.append((-1, xs))

    return data

#係数alphaの最適化を行う
#要件: 最後のサンプルのラベルが-1であること
def optimize(samples, kernel):

    #bの次元数
    dim = len(samples) - 1

    #alphaの代わりにbを最適化する
    #最適化問題: min f(b) = b^TGb - c^Tb
    
    c = mat(zeros((dim, 1)))
    for i in xrange(dim):
        c[i] = 1 + samples[i][LABEL]

    #kernel(x_i, x_n)を求める
    gn = zeros(dim + 1)
    for i in xrange(dim + 1):
        gn[i] = kernel(samples[i][VECTOR], samples[dim][VECTOR])
    
    G = mat(zeros((dim, dim)))
    for y in xrange(dim):
        for x in xrange(y, dim):
            G[y, x] = G[x, y] = kernel(samples[x][VECTOR], samples[y][VECTOR]) \
                - gn[x] - gn[y] + gn[dim]

    #bの初期化
    b = mat(zeros((dim, 1)))
    for i in xrange(dim):
        b[i] = 0.5 * samples[i][LABEL]
    
    before_b = mat(ones((dim, 1)))
    dist_eps = 1.0e-4
    nu = 1.0e-4
    iteration = 0
    max_iteration = 1e5

    while norm(b - before_b) > dist_eps and iteration < max_iteration:

        if iteration % 1000 == 0:
            print norm(b - before_b), iteration

        before_b = b.copy()

        #bを行列による一斉更新
        #d = dot(G, b) - c
        #b -= nu * d

        #bを逐次更新
        for i in range(dim):
            d = dot(G[i,:], b) - c[i]
            b[i] -= nu * d
        
        #制約による修正 (bi * yi >= 0)
        for i in range(dim):
            if samples[i][LABEL] * b[i] < 0:
                b[i] = 0

        #制約による修正 (sum(b) >= 0)
        sum_b = sum(b)
        if sum_b < 0:
            Np = len(filter(lambda item: item[LABEL] == 1, samples))
            for i in xrange(dim):
                if samples[i][LABEL] == 1:
                    b[i] -= sum_b / Np
                    
        iteration += 1

    print iteration
        
    #biをalpha_iに直す
    alphas = zeros(dim + 1)
    for i in xrange(dim):
        alphas[i] = samples[i][LABEL] * b[i]
    alphas[dim] = sum(b)

    #alpha_iの最大値より十分小さい要素は0
    alpha_threashold = max(alphas) * 1.0e-5
    for i in xrange(len(alphas)):
        if alphas[i] < alpha_threashold:
            alphas[i] = 0.0
            
    return alphas

#最適化したalphaを用いてサンプルの最適化を行う
def make_classifier(kernel, alphas, samples):

    def weighted_sum_kernel(x):
        ans = 0.0
        for i in range(len(samples)):
            ans += kernel(samples[i][VECTOR], x) * alphas[i] * samples[i][LABEL]
        return ans
        
    psamples = filter(lambda sample: sample[LABEL] == 1, samples)
    min_psection = min(map(lambda sample: weighted_sum_kernel(sample[VECTOR]), psamples))
    nsamples = filter(lambda sample: sample[LABEL] == -1, samples)
    max_nsection = max(map(lambda sample: weighted_sum_kernel(sample[VECTOR]), nsamples))
    section = - (max_nsection + min_psection) / 2.0
    
    def classify(x):
        return weighted_sum_kernel(x) + section

    return classify    

#カーネルの一例
def gauss_kernel(sigma):
    def kernel(x1, x2):
        return exp(-(norm(x1 - x2)) / ((2 * sigma) ** 2))
    
    return kernel

def inner_product():
    def kernel(x1, x2):
        return dot(x1.T, x2)
    return kernel

#境界面の関数の一例
def target(x):
    return sin(4 * x) / 2.0 + 0.25

def target2(x):
    return 1 - x

#実際の使用例
Nsample = 200

func = target

samples = make_input_data(Nsample, func)
kernel = gauss_kernel(0.5)
alphas = optimize(samples, kernel)

classify = make_classifier(kernel, alphas, samples)

for i in range(len(samples)):
    print samples[i][LABEL], classify(samples[i][VECTOR]), alphas[i]

class1 = [item[VECTOR] for item in samples if item[LABEL] == 1]
class2 = [item[VECTOR] for item in samples if item[LABEL] == -1]

xs1 = [v[0] for v in class1]
ys1 = [v[1] for v in class1]
xs2 = [v[0] for v in class2]
ys2 = [v[1] for v in class2]

plot(xs1, ys1, 'ro')
plot(xs2, ys2, 'bo')

xs = arange(0, 1, 0.01)

size = len(xs)

mesh_z = zeros((size, size))
xy = zeros(2)

for y in range(size):
    for x in range(size):
        xy[0], xy[1] = xs[x], xs[y]
        mesh_z[y, x] = classify(xy)

mesh_x, mesh_y = meshgrid(xs, xs)
contourf(mesh_x, mesh_y, mesh_z, 0)
show()

これが結果の一例で、上手くいけばこうなった。

SVMはやっぱすごいけど、実際使ったり実験するならやっぱパッケージ使うよなーと実感した

2012-11-11

Mona Lisa TSPでモナリザの顔を表示させる -その2-

前回：
http://d.hatena.ne.jp/xyz600/20121028/1351428339

前回swapする方法だと上手く行かなかったからKD-Treeを作って近傍に線を引いていこうと思ってたんだけど、どうやらinverseの法を高速化して工夫すると十分良いらしい(研究室の先輩談)。

inverseさせる方法というのは、単にランダムに2点を選択して2点間の経路を反転させるというもので、２点をswapさせる方法より辺の変化数が少ない(inverse = 2, swap = 4)分優秀とのこと。
頂点の数をnとするとナイーブに計算すると経路反転にO(n)かかるから、データ構造を作ることになるけど、経路反転を高速に出来るデータ構造は知らない...
しばらく探してたら、こんなのが見つかった。

http://www.slideshare.net/iwiwi/2-12188757/1

Treapって聞いたことなかったけど、スライドがわかりやすいから実装はそんなにムズくなかった。多分反転の操作にO(log n)かかるように見えるんだけど、どうなんだろう。
コードはこんな感じになった。ずいぶん投げやりだけど。

import random, math

class Point:
    def __init__(self, y, x):
        self.x = x
        self.y = y
        #ランダムに割り振られる優先度
        self.priority = random.random()

    def toString(self):
        return "(%d, %d)" % (self.y, self.x)
        
class Node:
    def __init__(self, lst):
        self.reverse = False
        self.size = len(lst)
        if len(lst) == 0:
            pass
        else:
            self.size = len(lst)
            max_value = max(lst, key = (lambda x: x.priority))
            index = lst.index(max_value)
            self.location = lst[index]
            self.left = Node(lst[:index])
            self.right = Node(lst[index+1:])

    def get(self, index):
        self.push()
        if index < self.left.size:
            return self.left.get(index)
        elif self.left.size < index:
            return self.right.get(index - self.left.size - 1)
        else:
            return self.location
            
    def is_null(self):
        return self.size == 0
            
    def is_leaf(self):
        return self.right.is_null() and self.left.is_null()
            
    def update(self):
        if not self.is_null():
            self.size = self.left.size + self.right.size + 1
        return self

    def push(self):
        if self.reverse and not self.is_null():
            self.left, self.right = self.right, self.left
            self.reverse = False
            self.right.reverse ^= True
            self.left.reverse ^= True

    def toString(self):
        self.push()
        if self.is_null():
            return "(NULL)"
        else:
            return "[%s{%s}%s]" % (self.left.toString(), self.location.toString(), self.right.toString())

    def toList(self):
        self.push()
        if self.is_null():
            return []
        elif self.is_leaf():
            return [self.location]
        else:
            return self.left.update().toList() + [self.location] + self.right.update().toList()
            
def merge(node_l, node_r):
    node_l.push()
    node_r.push()
    if node_l.is_null():
        return node_r
    elif node_r.is_null():
        return node_l
    elif node_l.location.priority > node_r.location.priority:
        node_l.right = merge(node_l.right, node_r).update()
        return node_l.update()
    else:
        node_r.left = merge(node_l, node_r.left).update()
        return node_r.update()
    
# [0, n) -> [0, k) + [k, n)
def split(node, k):
    node.push()
    if node.is_null():
        return (node, node)
    if k == node.size:
        return (node, Node([]))
    elif k == 0:
        return (Node([]), node)
    elif node.left.size < k:
        (node_rl, node_rr) = split(node.right, k - node.left.size - 1)
        node.right = node_rl.update()
        return (node.update(), node_rr.update())
    else:
        (node_ll, node_lr) = split(node.left, k)
        node.left = node_lr.update()
        return (node_ll.update(), node.update())

def reverse(node, frm, to):
    node.push()
    if node.is_null():
        return node
    else:
        left, middle_tmp = split(node, frm)
        middle, right = split(middle_tmp, to - frm + 1)
        middle.reverse ^= True
        return merge(left, merge(middle, right))

updateで木のサイズを再計算していて(ただし再帰的に計算するわけではなく、左右の部分木のサイズが正しく計算されているという仮定が必要)、pushでinverseの操作を実際に行っている。
updateとpushのタイミングがよく分からなかったんだけど、とりあえずオブジェクトに破壊的な操作をする時は終了したタイミングでupdateをかけて、オブジェクトのデータにアクセスする時にpushしてるんだけど、もう少し上手くできないものか。
これをpypyで3時間程度動かすと、1700万〜1800万程度の経路長が出てくる。これで実際の既知最良解の3、４倍程度か。前回より10倍近く短くなったからまぁいいや。
今回は割とモナリザっぽくなった！

次はどうしようかな。最近傍ってどの程度精度出るのか興味あるし、やっぱりKD-treeかなぁ。工夫すればinverseでももっといけるっぽいけど、悩み中。

2012-11-06

matplotlibで２軸グラフと標準偏差付きのグラフを書く

python

研究で使用したので、自分用のメモで。ここらへんを参考にしてみた。

http://chick.g.hatena.ne.jp/allegro/20091010/p5
http://matplotlib.org/examples/api/fahrenheit_celcius_scales.html
http://matplotlib.org/examples/pylab_examples/errorbar_demo.html
http://www.comp.tmu.ac.jp/shintani/japanese/openUniversity/node31.html
http://w.livedoor.jp/met-python/d/matplotlib#content_5_14

個人的に軸やラベルのフォントサイズがどうにかしたかったので、試行錯誤してたらこんなんでうまく行った。

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt

parameter_list = [(dim, kids_rate) for dim in [2, 3, 5, 10, 20, 50] for kids_rate in [6, 10, 20, 50, 100]]

for dimension, kids_rate in parameter_list:

    #prepare for data
    data_set = [map(float, line.strip().split(",")) for line in open("../result_d%d_kr%d.txt" % (dimension, kids_rate)).readlines()]
    Crs = [item[0] for item in data_set]
    successRates = [item[1] for item in data_set]
    means = [item[2] for item in data_set]
    variances = [item[3] for item in data_set]

    #prepare for plot
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax2 = ax1.twinx()

    fs_num = 16
    ax1.errorbar(Crs, means, fmt='-r', yerr=[variances, variances], lw = 4)
    plt.setp(ax1.get_xticklabels(), fontsize=fs_num, visible=True)
    plt.setp(ax1.get_yticklabels(), fontsize=fs_num, visible=True)
    ax2.plot(Crs, successRates, '-b', lw = 4)
    plt.setp(ax2.get_yticklabels(), fontsize=fs_num, visible=True)
    
    fs_label = 20
    ax1.set_xlabel('Cr', fontsize = fs_label)
    ax1.set_ylabel('# of generation', fontsize = fs_label)
    ax1.set_xlim(0.0, 0.9)
    ax2.set_ylabel('success rate', fontsize = fs_label)
    ax2.set_ylim(0.8, 1.0)
    ax1.grid(True)
    plt.savefig('graph_d%d_kr%d.eps' % (dimension, kids_rate))
    plt.clf()

実際のグラフはこんな感じ。

重要なのは、ax1.set_xlabelとかplt.setpあたりで、ここでフォントサイズの調整をしてる。
かなり無理やりな感じもあるが、上手くいったからしばらくこれでやってみる。

2012-10-28

Mona Lisa TSPでモナリザの顔を表示させる -その１-

python TSP

少し前、研究室の先輩からこんな問題を教えてもらった。

http://www.tsp.gatech.edu/data/ml/monalisa.html

巡回セールスマン問題(TSP)の最適解を解くと、モナリザの顔が浮かび上がってくるらしい。段々モナリザの顔が…みたいな事をやりたくなったので、最適解は無理でもモナリザの顔が見えるくらいは頑張ってみることに。
とりあえず、簡単な最適化アルゴリズムで様子をみることにした。とりあえず簡単そうなアルゴリズムとして、二つ選んだ。最初に訪問順序をランダムに初期化して、経路順序{x_1,x_2,...,x_n}（nは頂点数）が得られたとすると、

1. ランダムに選択した2頂点x_i、x_jの訪問順序を入れ替え、回路長が短くなったら入れ替えを保持。これを連続max_fail回失敗するまで続ける。
2. ランダムに選択した2頂点x_iからx_jまでの経路をひっくり返す。回路長が短くなったら入れ替えを保持。これを連続max_fail回失敗するまで続ける。

実行速度が心配だったけどpythonで書いた。コードはこんな感じ

import random, math

def solve(max_fail, max_iteration):
    data = read_from("mona-lisa100K.tsp")
    min_answer = data
    for i in xrange(max_iteration):
        local_answer = opt_swap(list(data), max_fail)
        if dist_all(min_answer) > dist_all(local_answer):
            min_answer = local_answer
    print dist_all(min_answer)
    write_to("answer.txt", min_answer)

def choose_index_pair(size):
    index1 = random.randint(0, size - 1)
    index2 = random.randint(0, size - 1)
    while index1 == index2:
        index1 = random.randint(0, size - 1)
        index2 = random.randint(0, size - 1)
    return (index1, index2)

def dist_all(points):
    return sum(dist(points[i], points[i - 1]) for i in xrange(len(points) - 1, -1, -1))

def dist(p1, p2):
    return math.sqrt((p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2)

def opt_inverse(data, max_fail):
    random.shuffle(data)
    iteration = 0
    continuous_fail = 0
    while continuous_fail < max_fail:
        iteration += 1
        index1, index2 = choose_index_pair(len(data))
        if can_improve_inv(data, index1, index2):
            max_i, min_i = max(index1, index2), min(index1, index2)
            data[min_i:max_i] = data[max_i - 1:min_i - 1:-1] if min_i > 0 else data[max_i - 1::-1]
            continuous_fail = 0
        else:
            continuous_fail += 1
        if iteration % 50000 == 0:
            print iteration
    return data

def can_improve_inv(data, index1, index2):
    l = len(data)
    nIndex1 = index1 + 1 if index1 < l - 1 else 0
    nIndex2 = index2 + 1 if index2 < l - 1 else 0
    before_dist = dist(data[index1 - 1], data[index1]) + dist(data[index2], data[nIndex2])
    after_dist  = dist(data[index1 - 1], data[index2]) + dist(data[index1], data[nIndex2])
    return before_dist > after_dist

def opt_swap(data, max_fail):
    random.shuffle(data)
    iteration = 0
    continuous_fail = 0
    while continuous_fail < max_fail:
        iteration += 1
        index1, index2 = choose_index_pair(len(data))
        if can_improve_opt_swap(data, index1, index2):
            data[index1], data[index2] = data[index2], data[index1]
            continuous_fail = 0
        else:
            continuous_fail += 1
        if iteration % 5000000 == 0:
            print iteration
    return data

def can_improve_opt_swap(data, index1, index2):
    l = len(data)
    nIndex1 = index1 + 1 if index1 < l - 1 else 0
    nIndex2 = index2 + 1 if index2 < l - 1 else 0
    before_dist = dist(data[index1 - 1], data[index1]) + dist(data[index1], data[nIndex1]) + \
        dist(data[index2 - 1], data[index2]) + dist(data[index2], data[nIndex2])
    after_dist = dist(data[index1 - 1], data[index2]) + dist(data[index2], data[nIndex1]) + \
        dist(data[index2 - 1], data[index1]) + dist(data[index1], data[nIndex2])
    return before_dist > after_dist

def write_to(filename, data):
    fout = open(filename, 'w')
    fout.writelines([" ".join(map(str, item)) + "\n" for item in data])
    fout.close()    

def read_from(filename):
    lst = []
    for line in open(filename).readlines():
        coordinate = line.split()
        lst.append((int(coordinate[1]), int(coordinate[2])))
    return lst

if __name__ == "__main__":
    max_fail = 20000
    max_iteration = 1
    solve(max_fail, max_iteration)

標準の機能のみだったから、素のpythonじゃなくてずっと高速なpypyで実行。どうも数分程度で同一の実行時間ならswapの方が優秀らしい。とりあえず、最適化後の経路長はおおよそ1.1〜1.2億程度になる。最短経路はどんなもんかなーと見てみると何と570万...

嫌な予感しかしないが一応プロットしてみると、

何がなんだかwちなみに点だけを打った本物はこんな感じ

さすがに乱択のみでは厳しかったか。次は貪欲法っぽく探索してない点を結んで行くとか
を実装する予定。素で計算すると計算量がちょっとヤバそうなので、KD-Treeあたりを実装する予定。

2012-04-21

xyz600の日記