2017-11-13

KaggleのTitanicをやってみたよ

問題概要

Kaggleの入門用として位置づけられているコンテスト？で、タイタニック号に乗っていた客のデータ（名前、年齢、客室の等級etc）が与えられるのでその生死を予測してくださいねというもの。実際の正解まで書いてある訓練用のデータと、正解は伏せられてるテスト用のデータが与えられるので、テスト用のデータに予測値を書き込んで提出してね、その正解率でスコアが出るよ、という感じです。

やったこと

1. ランダムフォレストの実装

これまでこの問題をやろうと思ってやめるを数か月繰り返してきたので、ランダムフォレストなる手法が有効であることはぼんやり把握していた。でググっていると診断人さんのブログにたどり着いた。

アルゴリズムの懇切丁寧な説明とC++の実装例が載っており、とりあえず自分でも実装してみて雰囲気を掴もうという気持ちになった。同じ言語で書いてもコピペになってしまいそうだったので競プロでよく使っているD言語で書いてみたが、まあ立ち位置的には似てるので結局割とコピペ感が強くなってしまった。

2. 欠損値の補完

ランダムフォレストができたので早速データをぶちこもうとすると、データ欠損の壁が立ちふさがる。この問題では一部項目でたまにデータがない部分があり、これをどうにかしないとちゃんとアルゴリズムを回すことができない。こういう欠損値への対処に関しては一般的にこれをやっておけば全部OKというようなやり方はないようで、都度アドホックなやり方でいい感じにデータを補完していく必要があるらしい（要出典）。この問題では特に年齢情報の欠損が多かったが、これはなんか名前にくっついてる敬称である程度推測可能らしい、ということを以下の記事を見て知った。

今回はとにかくさっき実装したランダムフォレストがちゃんと回るかを確かめたかったので、上の記事をそのまま参考に敬称内の平均年齢で欠損を埋めることにする。あとこまごました欠損があるが雰囲気で補った。

3. 提出

これでデータの体裁が整ったので、ランダムフォレストに投入する。データ数の少なさゆえに一瞬で終わった。早速提出してみたところ、スコア0.79904で約9000チーム中1500位くらいになった。

f:id:fluffyowl:20171113195521p:plain

自分で生み出したものが何もないのでどういう感情を持てばいいのか微妙だが、とりあえずランダムフォレストがちゃんと動いてるっぽいと言えそうな結果になったのは嬉しかった。

感想

競プロ志向が強いせいか、ランダムフォレストを実装しているときが一番楽しかった。まあ正直kaggleにおいてその辺の作業はかなりどうでもいい部分ではあると思うので、今後参加するときはもっと順位に効くところをちゃんとやっていければな～という感じです。とはいえ何かほんとにたくさんの情報がKernelに書き込まれていて、一体上位陣の優位性はどこから生まれているんだ…という疑問はある。ていうかまずは統計をちゃんと勉強しようね…

コード

以下に置きました。

2017-11-09

yukicoder No.511 落ちゲー～手作業のぬくもり～

yukicoder 並列二分探索問題記録

https://yukicoder.me/problems/no/511

問題概要

無限行W列のグリッドで、2人ゲームが行われる。このゲームでは2人交互に長方形のブロックを上から落とし、ブロックの各列は一番下か他のブロックにつくまで落下を続ける。ブロックを落とすことである列の高さがHを超えたとき、そのブロックを落としたプレイヤーに1列につき1ポイント入る。すべての列がHを越えたとき多くのポイントを獲得している方が勝者になる。あるゲームでどこにブロックが置かれたかの記録が順番にN個与えられるので、ゲームの勝敗を求めよ。

W, N <= 10⁵

H <= 10¹⁵

解法

Binary Indexed Tree で各列の高さを管理しつつ、並列二分探索をする。

まず列の高さをシミュレートするために区間add, 一点取得が高速でできるデータ構造が欲しい。これは区間add, 区間和取得が可能な普通のBITで実現できる。具体的にはimos法のイメージで、[l, r]に一律nを足すときはtable[l]にnを足し、table[r+1]に-nを足す。点xの現在値が知りたいときには[0, x]の区間和を取る。こうすると区間add, 一点取得ともにO(logW)でできるようになる。

次にある1列のみにおいてどちらが得点を得たかを判定してみる。これは前述のBITがあればブロックのクエリを二分探索してO(NlogNlogW)で求めることができる。

この二分探索を全部の列で行うことを考えると、普通にやったらO(NWlogNlogW)になって全く間に合わない。しかし各二分探索で探索する対象は同じであることを考えると、これを全部並列に行うことで計算量を落とすことができる（並列二分探索）。具体的には以下のような手順で行う。

各列について「今、二分探索のhigh, low, middleはどこか」を覚えておく。
BITで普通にブロックの積みをひとつずつシミュレートしていく
途中で「今ここがmiddleである」という列があれば、普通の二分探索のようにhigh, lowの再計算を行う。
ブロックを最後まで積み終わったら、また最初に戻って積むことを繰り返す。
全部の二分探索が収束したら終わり

個別の二分探索はlogN回で収束するため、ブロックの積みは高々logN回までしか繰り返されない。各繰り返し内でBITによるブロックのシミュレートはO(NlogW)かかるため、このパートはO(NlogNlogW)になる。また各二分探索の保持を平衡二分木とかでやると一回の処理がO(logW)ででる。そしてW個の二分探索それぞれがlogN回行われることになるので二分探索を捌くのはO(WlogNlogW)になる。2つのパートを合わせてO((N+W)logNlogW)となる。（たぶん）

感想

AGC002Dの解説で並列二分探索を知り、よくわからずにいろいろググってたらこの問題が並列二分探索の例題として挙げられていたのを見つけたので解いてみました。並列二分探索は頭良すぎて理解した瞬間感動しました。実装面の細かい話として最初各にぶたんを溜めておくところで連想配列内の動的配列に入れてたんですが、偏った入力で同じところへのpushがたくさん発生すると動的配列の伸長が起こりまくってクソ遅くなってしまったところ、平衡二分木に切り替えたらまともな速度になりました。そして何気にyukicoderの問題記録を書くのこれが初めての気がします。yukicoderも結構解いてるんですがブログに問題記録を書く習慣をつける前に適正レベルの問題を解きつくしてしまった感じなのであまり書く機会がありませんでした。今後もちょうどいいやつがあれば書いていきたいですね。この記事だけ異常に感想が長いな

コード (D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop;


void main() {
    auto s = readln.split;
    auto N = s[0].to!int;
    auto W = s[1].to!int;
    auto H = s[2].to!long;

    auto A = new int[](N);
    auto B = new long[](N);
    auto X = new int[](N);

    foreach (i; 0..N) {
        s = readln.split;
        A[i] = s[0].to!int;
        B[i] = s[1].to!long;
        X[i] = s[2].to!int - 1;
    }


    auto ans = new int[](W);
    auto rbt = new RedBlackTree!(Tuple!(int, int, int, int), "a[3] < b[3]", true)();

    foreach (i; 0..W)
        rbt.insert(tuple(i, -1, N-1, (N-2)/2));


    for (int settled = 0; settled < W; ) {
        BIT height = new BIT(W+10);

        foreach (i; 0..N) {
            height.add(X[i], X[i]+A[i], B[i]);
            auto nodes = rbt.equalRange(tuple(0, 0, 0, i)).array;

            foreach (node; nodes) {
                auto col = node[0];
                auto lo = node[1];
                auto hi = node[2];
                rbt.removeKey(node);

                if (hi - lo <= 1) {
                    settled += 1;
                    ans[col] = hi % 2 ? -1 : 1;
                } else {
                    if (height.sum(col) >= H) hi = (hi + lo) / 2;
                    else lo = (hi + lo) / 2;
                    rbt.insert(tuple(col, lo, hi, (hi+lo)/2));
                }
            }
        }
    }

    if (ans.sum > 0) writeln("A");
    else if (ans.sum == 0) writeln("DRAW");
    else writeln("B");
}


class BIT{
    int n;
    long[] table;

    this(int n){
        this.n = n;
        table = new long[](n+1);
    }

    void add(int i, long x){
        i++;
        while (i <= n){
            table[i] += x;
            i += i & -i;
        }
    }

    // [l, r)
    void add(int l, int r, long x) {
        add(l, x);
        add(r, -x);
    }

    // [0, i]
    long sum(int i){
        i++;
        long ret = 0;
        while (i > 0){
            ret += table[i];
            i -= i & -i;
        }
        return ret;
    }

    // [l, r)
    long sum(int l, int r){
        if (l >= r) return 0;
        return sum(r-1) - sum(l-1);
    }
}

2017-11-08

AtCoder Grand Contest 002 D - Stamp Rally

AtCoder AGC AGC #002 問題記録永続UnionFind

http://agc002.contest.atcoder.jp/tasks/agc002_d

問題概要

N頂点M辺の無向グラフがあり、頂点には1～N, 辺には1～Mの番号がついている。またQ個のクエリがあり、ひとつのクエリでは2つの頂点番号x, yと整数zが与えられる。各クエリごとに、x, yの2つを出発点として合計でz個の頂点を訪れるとき、通る辺の最大の番号の最小値を求めよ。

N, M, Q <= 10⁵

解法

まずクエリが1回だけの場合を考える。最初グラフのすべての辺がなくなったものとして、ここに番号の小さい順から辺を追加していくことを考えると、この操作の過程でx, yから訪れることができる頂点数が初めてzを超えるタイミングが答えになることがわかる。具体的にはUnionFindを用いて、(1)もしxとyが同じ連結成分に入っていればその連結成分のサイズ (2)違う連結成分であればそれぞれのサイズの合計を求めれば「x, yから訪れることができる頂点数」が割り出せるので、UnionFindで辺をつなぐたびこの値を計算すればよい。

この方法はクエリが複数回のときにも通用しそうに見える。つまりUnionFindで辺をつないでいく操作は各クエリの内容に関係なく共通なので、前計算しておけばよさそうな気がする。ただしここで問題があり、普通のUnionFindでは辺をつなぐ操作は非可逆的に行われるため「i番目の辺をつないだときの状態」を取り出すことができない。もちろん愚直には辺を追加した時点でコピーを作って取っておけば後から見ることができるが、これは計算量がかかりすぎる（時間も空間も）。

これを解決するのが （部分）永続UnionFind である。このデータ構造では空間計算量(N+Q)で以下の3つの操作をいずれもO(logN)で行うことができる。

unite(u, v): 頂点u, vをつなぐ
find(u, t): t回目のuniteを行った時点での頂点uの親
size(u, t): t回目のuniteを行った時点での頂点uが属する連結成分のサイズ

（仕組みとか実装に関してはhttps://camypaper.bitbucket.io/2016/12/18/adc2016/ がとてもわかりやすかったのでこちらを見てください）

これで各辺を追加した段階での特定頂点が属する連結成分のサイズをO(logN)で計算できるようになった。あとは各クエリごとに2分探索でzを超える瞬間の辺番号を求めればよい。合わせてO(Q(logN)²)となる。

感想

初めて永続的データ構造の類を使ったけど思ったよりゴテゴテしてなくてよかった（？）想定解はこれじゃなくて並列二分探索なるものらしいので要理解

コード (D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop;
 
const int INF = 1 << 29;
 
void main() {
    auto s = readln.split.map!(to!int);
    auto N = s[0];
    auto M = s[1];
    auto uf = new PersistentUnionFind(N);
    foreach (i; 0..M) {
        s = readln.split.map!(to!int);
        uf.unite(s[0]-1, s[1]-1);
    }
 
    auto Q = readln.chomp.to!int;
    while (Q--) {
        s = readln.split.map!(to!int);
        auto x = s[0]-1;
        auto y = s[1]-1;
        auto z = s[2];
        int hi = M;
        int lo = 0;
        while (hi - lo > 1) {
            int size;
            int mid = (hi + lo) / 2;
            if (uf.find(x, mid) == uf.find(y, mid))
                size = uf.size(x, mid);
            else
                size = uf.size(x, mid) + uf.size(y, mid);
            if (size >= z)
                hi = mid;
            else
                lo = mid;
        }
        hi.writeln;
    }
}
 
class PersistentUnionFind {
    int[][] rank;
    int[][] time;
    int[][] parent;
    int n;
    int global_time;
 
    this(int n) {
        this.n = n;
        rank = new int[][](n);
        time = new int[][](n);
        parent = new int[][](n);
        foreach (i; 0..n) {
            rank[i] ~= 1;
            time[i] ~= 0;
            parent[i] ~= i;
        }
        global_time = 0;
    }
 
    void unite(int u, int v) {
        global_time += 1;
        u = find(u, global_time);
        v = find(v, global_time);
        if (u == v) return;
 
        if (rank[u] < rank[v]) swap(u, v);
 
        int r = rank[u].back + rank[v].back;
 
        rank[u] ~= r;
        time[u] ~= global_time;
        parent[u] ~= u;
 
        rank[v] ~= r;
        time[v] ~= global_time;
        parent[v] ~= u;
    }
 
    int find(int u, int t) {
        if (parent[u].back == u) return u;
        if (time[u].back > t) return u;
        return find(parent[u].back, t);
    }
 
    int size(int u, int t) {
        int v = find(u, t);
        int i = time[v].assumeSorted.lowerBound(t+1).length.to!int;
        return rank[v][i-1];
    }
}

2017-11-08

Typical DP Contest M - 家

AtCoder TDPC 問題記録

http://tdpc.contest.atcoder.jp/tasks/tdpc_house

問題概要

H階建ての家があり、各階には部屋が1番～R番までのR個存在する。部屋間を以下のルールに従って移動するとき、H階の部屋1から1階の部屋1に移動する経路は何通りあるか。ルール1. h階の部屋rから(h-1)階の部屋rへ移動できる。ルール2. 同じ階の部屋間は、与えられる行列Gの対応要素が1であれば相互に移動可能。ルール3. 同じ部屋を2度通ることはできない。

H <= 10⁹

R <= 16

解法

横方向の移動の仕方をbitDPで計算した後、縦方向の移動の仕方を行列累乗で求める。

「その階で最初に部屋xにいて最後に部屋yにいたときの移動経路の総数」をすべてのx, y(1 <= x, y <= R)の組合せで求めると、これをR行R列の行列の形で表すことができる。すべての階で部屋の構造は同じなので、あとはH階分の行列累乗を行えばH階→1遷への遷移を求めることができる。

問題は「その階で最初に部屋xにいて最後に部屋yにいたときの移動経路の総数」をどう求めるかということであるが、これはbitDPで行うことができる。具体的には

dp(i, j, mask): 「その階で最初にいた部屋がiで、部屋集合maskをこれまでに通過し、最後にjにいる」ときの移動経路の総数

となる。遷移の際に経由する部屋を総当たりする必要があるのでこのdpの計算はO(R³×2^R)になる。行列累乗パートも考慮に入れると計算量はO(R³×2^R + R³×logH)となる。R=16だと若干怪しくなってくる計算量だが実はTDPCでこの問題だけTLが8secもあるので余裕を持って間に合う。

感想

最初2^R×R³の桁を一個見間違えて無駄に迷走してしまった……

コード (D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop;
 
immutable long MOD = 10^^9+7;
 
void main() {
    auto s = readln.split.map!(to!int);
    auto H = s[0];
    auto R = s[1];
    auto G = R.iota.map!(_ => readln.split.map!(to!int).array).array;
 
    
    auto dp = new long[][][](R, R, 1 << R);
    foreach (i; 0..R) dp[i][i][1<<i] = 1;
    auto masks = (1<<R).iota.array.sort!((a, b) => a.popcnt < b.popcnt).array;
 
    foreach (i; 0..R) // 出発
        foreach (mask; masks)
            foreach (j; 0..R) // 経由
                foreach (k; 0..R) { // 到着
                    if (!(mask & (1 << i))) continue;
                    if (!(mask & (1 << j))) continue;
                    if ((mask & (1 << k))) continue;
                    if (!G[j][k]) continue;
                    (dp[i][k][mask|(1<<k)] += dp[i][j][mask]) %= MOD;
                }
 
    
    auto mat = new long[][](R, R);
    foreach (i; 0..R)
        foreach (j; 0..R)
            foreach (mask; 0..1<<R)
                (mat[i][j] += dp[i][j][mask]) %= MOD;
 
    matpow(R, H, mat)[0][0].writeln;
}
 
 
long[][] matmul(int N, long[][] m1, long[][] m2) {
    auto ret = new long[][](N, N);
    foreach (i; 0..N) {
        foreach (j; 0..N) {
            ret[i][j] = 0;
            foreach (k; 0..N) {
                (ret[i][j] += m1[i][k] * m2[k][j] % MOD) %= MOD;
            }
        }
    }
    return ret;
}
 
long[][] matpow(int N, long K, long[][] mat) {
    auto ret = new long[][](N, N);
    foreach (i; 0..N)
        foreach (j; 0..N)
            ret[i][j] = i == j ? 1 : 0;
    
    while (K > 0) {
        if (K % 2 == 1)
            ret = matmul(N, ret, mat);
        mat = matmul(N, mat, mat);
        K /= 2;
    }
 
    return ret;
}

2017-11-07

AtCoder Regular Contest 084 D - Small Multiple

AtCoder ARC ARC #084 問題記録

http://arc084.contest.atcoder.jp/tasks/arc084_b

問題概要

正整数Kのある正の倍数を10進法で表したときの各桁の和の最小値を求めよ。

2 <= K <= 10⁵

解法

以下のようなグラフを考え、最短経路問題を解く。

頂点は0～(K-1)のK個
頂点x → 頂点(x+1)%K にコスト1の辺が出る
頂点x → 頂点(x×10)%K にコスト0の辺が出る
コスト1を持って頂点1からスタートし、頂点0をゴールとする

K <= 10⁵ よりダイクストラで解いてもO(KlogK)で間に合う。

この最短経路問題で解いているのは「mod K が x であるすべての正整数の集合における各桁の和の最小値」である。つまり0～(K-1)の各頂点が「mod K で x である正整数の集合」を表し、頂点1から頂点xまでの最短距離が「各桁の和の最小値」を表す。イメージとしては、「1」という数字からスタートしてその数に＋1するか×10するかの遷移を繰り返している、という感じになる。頂点0が「modKが0であるすべての正整数の集合」、すなわちKの正の倍数ということになるので、頂点0までのコストが解となる。

感想

mod Kで0の集合 = Kの倍数の集合というの自明に見えるけど、なんか頭の中で結びついてなかったので認識しておきたいですね

コード (D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop;

const int INF = 1 << 29;

void main() {
    auto N = readln.chomp.to!int;

    auto dist = new int[](N);
    dist.fill(INF);

    auto pq = new BinaryHeap!(Array!(Tuple!(int, int)), "a[1] > b[1]");
    pq.insert(tuple(1, 1));

    while (!pq.empty) {
        auto n = pq.front[0];
        auto c = pq.front[1];
        pq.removeFront;
        if (dist[n] <= c) continue;
        dist[n] = c;
        if (dist[(n+1)%N] > c+1) pq.insert(tuple((n+1)%N, c+1));
        if (dist[(n*10)%N] > c) pq.insert(tuple((n*10)%N, c));
    }

    dist[0].writeln;
}

2017-11-04

DDCC2017本戦に参加しました

参加記

11/3に行われた「DISCO presents ディスカバリーチャンネルコードコンテスト2017」(DDCC2017)の本戦に参加しました。新卒枠100人＋それ以外100人という非常に大規模なオンサイトコンテストで、枠の大きさに助けられて本戦に行くことができました。社会人にとってここまで本戦チャンスが大きいオンサイトコンテストは他にないので本当にありがたい限りです。ありがとうございます。ありがとうございました。

コンテスト

A, Bで2完の75位という何ともいえない結果でした。始まってまず思ったのが、Aからむずない？ということで、最終的に角だけ見ればええやんというところにたどり着くまで結構時間がかかってしまい焦りました。Bはサンプル見てパッと出てくるGCD->LCMを書いて通したけど何でこれでいいのかはいまだによくわかってないです。C問題はあまりにもわからないのでサイクルの列挙方法についてググったり、窓の外のイトーヨーカドーを眺めたりしていました。電源が少ないということでバッテリーが不安だったんですが、1時間半その状態だったのでメチャクチャ余裕でした。

イベント

コンテスト終了後にもイベントがあり、Ponanzaの山本一成さん、プロ棋士の西尾明さん、司会としてAtCoderの高橋直大さん、の対談企画がありました。将棋やらないので普段聞く機会のないような新鮮な話を多く聞けて楽しかったですが、やはり個人的に一番興味深かったのは既にレーティング1000くらいの実力はあるという競プロAIの話でした。あと社内見学ということで、自分は初参加なのでAコースに参加して機械とかいろいろ見せてもらいました。結構な時間が福利厚生施設の紹介に当てられており温かい気持ちになりました。

交流

コンテストが終わったあと、kosakkunさんに声をかけてもらって色々喋りました（話してる途中で判明したが奇跡的に順位が一つ違いだった）。他にもiwashi31さん、kuusoさんなどツイッターでフォローしている方々やkosakkunさんのお友達とも少しお話ができて望外の楽しい時間を過ごすことができました。普段周りに競プロどころかプログラミングやってる人すらマジで全然存在しないので、声に出して競技プログラミングのことを喋れるのこんなに楽しいんかいという感じでした。皆さんありがとうございました！

2017-10-27

Topcoder SRM 714 Div.1 Easy - ParenthesisRemoval

問題記録 Topcoder SRM SRM #714

問題概要

開き括弧と閉じ括弧のみからなる文字列Sが与えられる。初めSの括弧は正しく対応している。Sの一番前にある開き括弧と任意の閉じ括弧のペアを括弧の対応が不正にならないよう取り除く、という操作をSが空になるまで繰り返すとき、括弧の選び方の順番が何通りあるか答えよ。すべての括弧は一番最初のインデックスによって区別されるものとする。

|S| <= 2500

解法

閉じ括弧の側から考える。Sを前から見ていって最初の閉じ括弧に当たったとき、その閉じ括弧と一緒に消せる可能性があるのはそれより前にある開き括弧だけである。よって以下の手順で計算を行える。Sを前から見ていき、その文字が開き括弧ならカウンターを1増やす、閉じ括弧ならばカウンターを1減らす。閉じ括弧に当たった時点でのカウンターの数を掛け合わせたものが答え。O(|S|)

感想

視点を変えると一瞬で答えが出て面白かった

コード (C++)

#include <bits/stdc++.h>
using namespace std;
#define REP(i,n) for (int i=0;i<(n);i++)
#define REP2(i,m,n) for (int i=m;i<(n);i++)
typedef long long ll;

const ll MOD = 1000000007;

class ParenthesisRemoval {
public:
    int countWays(string s) {
        int N = s.length();
        ll ans = 1;
        ll open = 0;
        
        REP(i, N) {
            if (s[i] == '(') {
                open += 1;
            } else {
                ans = ans * open % MOD;
                open -= 1;
            }
        }
        
        return (int)ans;
    }
};