2016-03-30

k平均法（k-meansクラスタリング）実装

C# 機械学習

Courseraの機械学習コースでk平均法を学んだので、理解を深めるために実装してみた。
わざわざ日記にするほど難しくはないが、メモということで。ソースコードはここ。

メインは次のGetKMeans_withCost()で、tryNumの数だけk平均法を試して、コスト(最寄りのセントロイドまでの距離の合計)が最小になった結果を返している。返り値は各セントロイドの位置と最終的なコストのタプル。多次元対応。

入力はクラスタ数kとMatrix型の変数x。xは行が１サンプル、列が特徴を表す。

アルゴリズムは単純で、　　

xからk個ランダムに選び出し、この位置をセントロイドの初期値にする
それぞれのxに、最も近いセントロイドを割り当てる
セントロイドの位置を、割り当てられているxの平均値にアップデート

を繰り返しているだけ。セントロイドの位置が変わらなくなるまでループする。

static Tuple<int[], double> GetKMeans_withCost(int tryNum, int K, Matrix x)
{
    var r = new Random();

    int[] ret = null;
    var centroids = new Matrix[K];

    var minCost = double.MaxValue;
    while (tryNum > 0)
    {
        Initialize(r, x, centroids);

        var J = 0.0;
        int[] c = null;
        while (true)
        {
            var cTuple = AssignC(centroids, x);
            c = cTuple.Item1;

            var prevJ = J;
            J = cTuple.Item2;
            if (DoubleUtil.Equals(prevJ, J)) break;

            UpdateCentroids(centroids, x, c);
        }

        if (J < minCost || ret == null)
        {
            ret = c;
            minCost = J;
        }

        tryNum--;
    }

    return new Tuple<int[], double>(ret, minCost);
}

こちらはkの値をいろいろ試してみるコード。kの候補値(Ks[])ごとの結果を返す。
うまい具合にエルボー（これ以上増やしてもあまりコストが下がらないようなkの値）があれば、それを使うのが良いらしい。

public static double[] FindK(int tryNum, int[] Ks, Matrix x)
{
    var ret = new double[Ks.Length];

    for (int i = 0; i < Ks.Length; i++)
    {
        var result = GetKMeans_withCost(tryNum, Ks[i], x);
        ret[i] = result.Item2;
    }

    return ret;
}

2016-03-25

ニューラルネットワーク実装

C# 機械学習

Courseraの機械学習コースでニューラルネットワークを学んだので、習作としてC#で実装してみた。

多層パーセプトロン対応。Classification専用（シグモイド関数をベタ書きしてるので）。
ソースコードはここにおいてある。

実装の概略
ネットワーク各層の関数を、行列型の配列Thetas[]で保持している。

public class NeuralNetwork
{
    public Matrix[] Thetas { get; set; }
    public int NumLayers { get; set; }
    public int[] NumNeurons { get; set; }

これに、（RandomizeThetas()で初期化したあと）関数Learn()に学習用データとパラメタを渡せば最急降下法で学習を開始する。

public void Learn(Matrix[] input, Matrix[] output, double alpha, double lambda, int maxItr)
{
    while (maxItr > 0)
    {
        var grad = GetGrad(Thetas, input, output, input.Length, lambda);
        UpdateThetas(input.GetLength(0), grad, alpha, lambda, input, output);
        maxItr--;
    }
}

上記のGetGrad()が、各ノードの偏微分を返す関数で、それを使ってUpdateThetas()内でThetas[]を更新する。

GetGrad()の実装も教科書通り。

public Matrix[] GetGrad(Matrix[] thetas, Matrix[] input, Matrix[] output, int m, double lambda)
{
    var L_deltas = InitializeLDeltas();

    for (int idx = 0; idx < input.GetLength(0); idx++)
    {
        Matrix[] z = null;
        var values = ForwardProp(input[idx], thetas, ref z);
        BackProp(L_deltas, output[idx], values, z);
    }

    var grad = GetGrad(thetas, L_deltas, input.Length, lambda);
    return grad;
}

ForwardProp()は名前の通りフォワードプロパゲーションを、BackProp()はバックプロパゲーションをする関数。
ForwardProp()は各層の中間値を返すので、それを使ってBackProp()で後ろの層から変化量もとめ、L_deltasを更新している。
全入力のループが終わったら、GetGrad()がL_deltasを使って偏微分を求める。

偏微分を求める実装がややこしいのだが、コースで推奨されているとおりに、近似値を求めて比較する方法でデバッグした。
${\theta}$ -eと ${\theta}$ +e（eは小さい数）のそれぞれでコストを求め、その傾きを使うというもの。

Matrix GetNumericGradient(NeuralNetwork nn, Matrix[] thetas, int idx, Matrix[] input, Matrix[] output)
{
    var t = thetas.ToArray();

    const double e = 1e-4;
    var ret = new Matrix(t[idx].RowNum, t[idx].ColNum);
    var diffMat = new Matrix(t[idx].RowNum, t[idx].ColNum);

    for (int row = 0; row < diffMat.RowNum; row++)
    {
        for (int col = 0; col < diffMat.ColNum; col++)
        {
            var orgValue = t[idx][row, col];

            t[idx][row, col] = orgValue - e;
            var loss1 = nn.J(thetas, input, output);

            thetas[idx][row, col] = orgValue + e;
            var loss2 = nn.J(thetas, input, output);

            ret[row, col] = (loss2 - loss1) / (2.0 * e);

            thetas[idx][row, col] = orgValue;
        }
    }

    return ret;
}

var numgrads = GetNumericGradients(nn, nn.Thetas, x, y);
var grads = nn.GetGrad(new Matrix[] { theta1, theta2 }, x, y, x.Length, lambda);

2通りの結果がほぼ同じなので、実装が合っていることを確認できる。

使ってみる

こちらのサイトにある手書き数字の学習と認識に挑戦してみる。

static void NeuralNetworkTest_Example()
{
    Console.WriteLine("**Neural Network test (test data from: http://archive.ics.uci.edu/ml/machine-learning-databases/pendigits/)**");

    var nn = new NeuralNetwork(3, new int[] { 16, 16, 10 });
    nn.RandomizeThetas();

    //train
    var trainingDataTuple = GetNNData(TraininData.Data);
    var inputs = trainingDataTuple.Item1;
    var outputs = trainingDataTuple.Item2;

    //test
    var testDataTuple = GetNNData(TestData.Data);
    var test_inputs = testDataTuple.Item1;
    var test_outputs = testDataTuple.Item2;

    var alphas = new double[] { 5, 7 };
    var lambdas = new double[] { 0.03, 0.1, 0.3, 1.0 };

    var parameters = nn.FindParameters(alphas, lambdas, inputs, outputs, test_inputs, test_outputs, 300);
            
    Console.WriteLine("Chosen params - alpha:{0}, lambda:{1}", parameters.Item1, parameters.Item2);

    nn.Learn(inputs, outputs, parameters.Item1, parameters.Item2, 300);

    var result = nn.GetResult(test_inputs, test_outputs);
    var err = result.Item1;
    var total = result.Item2;

    Console.WriteLine(string.Format("Total:{0}, Error:{1}, Success Rate: {2}", total, err, (double)(total - err) / (double)total));            
}

トレーニング用とテスト用のデータがすでに用意されていたので、そのまま利用した。
好ましいパラメタを選ぶために、本来ならクロスバリデーション用として別のデータを用意すべきだが、今回はテスト用を流用した。

3層ネットワークで、イテレーションは適当に300回。
合計5回の学習で1時間以上かかってしまった。

f:id:yambe2002:20160324232921p:plain

認識率は約90%だった。何も工夫しないとこんなものだろうか。または、イテレーション回数が少なくて収束しきっていないのかもしれない。
このあたりは、別にちゃんとしたライブラリを使っていろいろ試してみようと思う。

一から自分で組んでみてかなり勉強になったが、本格的に実装するならライブラリを使わないと難しいと改めて実感。（数学的に最適化されてないので実行がかなり遅いし、デバッグも難しい）

2016-03-03

HackerRank Week of Code - 19 参加日記

プログラミングコンテスト C#

HackerRankのWeek of Code - 19に参加したのでそのメモ。
結果は346位/3176人だった。

Fix the Cycles

4つのノードA,B,C,Dをもつ有向グラフを考える。D→Aのエッジをa、A→Bのエッジをb、B→Cのエッジをc、C→Dのエッジをd、A→Cのエッジをe、B→Dのエッジをfとする。
a~fの重みが与えられたとき、負の閉路をなくすためには、最小でいくら重みを足せばよいか。ただし、1つのエッジにのみ重みを足すことができる。
不可能な場合は-1を返こと。
制約：a~fの重みの値は-20から+20

全エッジに対し、足す重み全て(0~80)をチェックするだけでよい。
ちなみに、すべてのループがエッジaを通るため、答えが-1になることはない。

static bool Success(int a, int b, int c, int d, int e, int f, int val, int idx)
{
    switch (idx)
    {
        case 0: a += val; break;
        case 1: b += val; break;
        case 2: c += val; break;
        case 3: d += val; break;
        case 4: e += val; break;
        case 5: f += val; break;
    }

    return (b + f + a >= 0) && (b + c + d + a >= 0) && (e + d + a >= 0);
}

public static int Solve(int a, int b, int c, int d, int e, int f)
{
    for (int i = 0; i <= 80; i++)
    {
        for (int j = 0; j < 6; j++)
        {
            if (Success(a, b, c, d, e, f, i, j)) return i;
        }
    }

    return -1;  //should not be here since all loops have edge 'a'
}

Two Robots

無限個のキャンディーが入ったM個のコンテナが、１列に配置されている。そして2つのロボットがある。ロボットはコンテナ間を移動して、キャンディーを運ぶことができる。
キャンディー運ぶ順番（あるコンテナから別のコンテナへ）が指定されたとき、トータルで最小のロボット移動距離を求めよ（隣り合うコンテナ間の距離を距離1とする）。

キャンディーを運ぶ順番は次のように指定される。
N　#運ぶ回数
Ma1 Mb1　#コンテナMa1からMb1に運ぶ。必ず違う値
・・・
Man　Mbn

制約： ${ 1 \leq M \leq 1000, \, 1 \leq N \leq 1000 }$

次のDPを設定し、手順を後ろからループしていく。
DP[現在のロボットAの位置][現在のロボットBの位置] ＝以降の移動距離の最小（Aの位置＜Bの位置）

単純に実装するとメモリと実行時間が足りなくなるので、(1)DP配列の再利用(2)あり得ない条件は無視して枝刈り、を行った。

public class Solver
{
    public int n = 0;
    public int m = 0;
    public Tuple<int, int>[] queries;

    public int Solve()
    {
        int[][,] dp_all = new int[][,] { new int[m + 2, m + 2], new int[m + 2, m + 2] };

        int[] minIdx = new int[m + 1];
        for (int idx = 0; idx < m + 1; idx++)
        {
            minIdx[idx] = Int32.MaxValue;
        }

        for (int idx = 0; idx < queries.Length; idx++)
        {
            var val = queries[idx].Item2;
            if (minIdx[val] == Int32.MaxValue || idx < minIdx[val])
            {
                minIdx[val] = idx;
            }
        }

        int dpIdx = 0;
        for (int idx = queries.Length - 1; idx >= 0; idx--)
        {
            var q = queries[idx];

            var pos1 = idx == 0 ? 0 : queries[idx-1].Item2;

            for (int pos2 = 0; pos2 < m + 1; pos2++)
            {
                if (pos2 != 0 && idx < minIdx[pos2]) continue;  //cannot be

                var dp = dp_all[dpIdx];
                var dp_next = dp_all[dpIdx == 0 ? 1 : 0];

                dp[Math.Min(pos1, pos2), Math.Max(pos1, pos2)] = Math.Abs(q.Item1 - q.Item2);

                var movePos1 = (pos1 == 0 ? 0 : Math.Abs(pos1 - q.Item1))
                    + ((idx == queries.Length - 1) ? 0 : dp_next[Math.Min(q.Item2, pos2), Math.Max(q.Item2, pos2)]);

                var movePos2 = (pos2 == 0 ? 0 : Math.Abs(pos2 - q.Item1))
                    + ((idx == queries.Length - 1) ? 0 : dp_next[Math.Min(q.Item2, pos1), Math.Max(q.Item2, pos1)]);

                dp[Math.Min(pos1, pos2), Math.Max(pos1, pos2)] += Math.Min(movePos1, movePos2);
            }

            dpIdx = dpIdx == 0 ? 1 : 0;
        }

        return dp_all[dpIdx == 0 ? 1 : 0][0, 0];
    }
}

Scalar Products

次の数列を考える。
${ a_0 = 0, a_1 = C, a_{i+2} = (a_{i+1} + a_i) \% M,\,where\,0 \leq i \leq 2n}$
このとき、ベクトル ${v_1, v_2, ..., v_n}$ を次のように規定する。
${v_1 = (a_2, a_3), v_2 = (a_4, a_5), ..., v_n = (a_{2n}, a_{2n+1})}$
Sを ${v_i, v_j \, (v\leq i, j \leq n)}$ のスカラ積の集合としたとき、Sの剰余の個数をMでModしたものを求めよ。
制約： ${1 \leq C \leq 10^9, 1 \leq M \leq 10^9, 1 \leq n \leq 3 \times 10^5}$

${a_0 = 0,\, a_1 = C }$ で開始するフィボナッチ数列の変種と考えればよい。

フィボナッチ数列の特徴
${ a_x \times a_y + a_{x+1} \times a_{y+1} = a_{x+y+1} }$
を利用すると、この数列に対しては
${ a_x \times a_y + a_{x+1} \times a_{y+1} = C \times a_{x+y+1} }$
が導き出せる。

よって、 ${ C \times a_{x+y+1} }$ の部分を見ていけばよい。

public class Solver
{
    public long Solve(long c, long m, long n)
    {
        c %= m;
        var resultSet = new HashSet<long>();

        var prev = c;
        long prev2 = 0;
        for (int i = 2; i <= 4 * n - 1; i++)
        {
            var val = (prev2 + prev) % m;

            if (i % 2 == 1 && i >= 7) resultSet.Add((val * c) % m);

            prev2 = prev;
            prev = val;
        }

        return resultSet.Count() % m;
    }
}

これは解けなかったのが悔しかった。
「フィボナッチ数列」で検索してみるべきだった。

Coolguy and Two Subsequences

一次元配列A[]（長さN、Indexは1から始まる）が与えられたとき、次の疑似コードで得られるansを答えよ。
//f(a, b)はA[a, b]の最小値を返す
ans = 0
for a -> [1, N]
    for b -> [a, N]
        for c -> [b + 1, N]
            for d -> [c, N]
                ans = ans + min(f(a, b), f(c, d))

ans = ans % (10^9 + 1)
制約： ${ 1 \leq N \leq 2 \times 10^5, \, 1 \leq A_i \leq 10^9 }$

まった手が出なかった。
Editorialは理解できるのだが、最後が「これをSegmentTreeを使って解けばよい」で終わっていて困る・・・。このレベルなら自明であるべきなのだろうか。

仕方ないので、SegmentTreeの集中特訓をしてから、後でもう一度解いてみることにする。

Editorialの内容

配列Aにあるすべての数字xについて、 ${ min( f(a, \, b), \, f(c,\,d) ) = x }$ になる回数を求めればよい。
これを ${ cnt(x) }$ とすると、最終的な解答は ${ sum( x \times cnt(x) ) }$ になる。

ここで、
　 ${g(x) = min( f(a, \, b), \, f(c,\,d) ) \geq x になる回数}$
と定義すると、
${cnt(x) = cnt(x) - cnt(x + 1)}$ となる。

${g(x)}$ を求めるためには、別の配列Bを用意するとよい。
BはAと同じ長さで、 ${a_i \geq x}$ なら ${b_i = 1}$ 、それ以外なら ${b_i = 0}$ のもの。

すると、 ${g(x)}$ は
　 ${x1 \leq y1}$ かつ
　 ${y1 + 1 \leq x2}$ かつ ${ x2 \leq y2 }$
になる ${(x1,\,y1,\,x2,\,y2)}$ の組み合わせのうち、
${ B[x1:y1]}$ と ${ B[x2:y2] }$ の範囲内すべてが1になるものの数と等しくなる。

${g(x)}$ はSegmentTreeを使えば解くことができる。

2016-02-13

線形回帰とロジスティック回帰の実装

C# 機械学習

Courseraの機械学習コースで多変量線形回帰（最小二乗法）とロジスティック回帰を学んだので、理解を深めるためにC#で実装してみた。
特にロジスティック回帰は、コースで最急降下法の実装をやらなかった（Octave組み込みのfminuc()を使う）ため、練習として組んでみたかった。

行列クラス
まずは行列クラスを作成。
topcoder_template/MyLibrary_Matrix.cs at master · yambe2002/topcoder_template · GitHub

演算子のオーバーロードと、平均・標準偏差や合計を求める関数などもあり、ずいぶん大きくなってしまった。

逆行列を求める関数Inverse()は掃き出し法を使っているので、一般化逆行列を返すわけではない。
NormalizeFeature()は特徴量を標準化させるもので、行列クラスに入れるのは微妙だったかもしれない。

線形回帰

最急降下法（Gradient Descent）

各特徴量ごと、コストをその特徴量(j)にかかる係数 ${\theta_{j}}$ で偏微分し、コストが下がる方向にその係数 ${\theta_{j}}$ を変化させていく。返り値は ${\theta}$ のベクトルとコスト履歴のタプル。

public static Tuple<Matrix, double[]> GetResult_GradientDescent(Matrix x, Matrix y, double alpha, double lambda, int maxItr, bool logCosts=true)
{
    var theta = new Matrix(x.ColNum, 1);
    var costs = new List<double>();
    while (maxItr > 0)
    {
        var h = x * theta;

        if (logCosts)
        {
            var j = ((1.0 / (2.0 * (double)x.RowNum)) * ((h - y) ^ 2)).Sum();
            costs.Add(j);
        }
        var diff = (alpha / (double)x.RowNum) * (x.Transpose() * (h - y)) - ((alpha * lambda) / (double)x.RowNum) * theta;
        diff[0, 0] += ((alpha * lambda) / (double)x.RowNum) * theta[0, 0];
        theta = theta - diff;

        maxItr--;
    }

    return new Tuple<Matrix,double[]>(theta, costs.ToArray());
}

行列xがトレーニングデータ（データ数 ${ \times }$ 特徴数）、行列yが解のベクトル、alphaが降下パラメータ、lambdaが標準化係数。

正規方程式（Normal Equation）

一般解として計算する方法。 ${\theta}$ を求める係数のベクトルとすると、解は
${ \theta = (X ^{T} X) ^{-1} X^{T} Xy }$ 　になる。
なぜこれで求まるかは、線形回帰の Normal Equation（正規方程式）について - Qiitaが参考になった。
とても簡単で結果も厳密だが、データサイズが大きくなると計算に時間がかかる。

public static Tuple<Matrix, double[]> GetResult_NormalEquation(Matrix x, Matrix y, double lambda, bool logCosts = true)
{
    var xT = x.Transpose();
    var l = Matrix.Eye(x.ColNum);
    l[0, 0] = 0;

    var ret = (xT * x + lambda * l).Inverse() * xT * y;

    var costs = new List<double>();
    if (logCosts)
    {
        var h = x * ret;
        var j = ((1.0 / (2.0 * (double)x.RowNum)) * ((h - y) ^ 2)).Sum();
        costs.Add(j);
    }

    return new Tuple<Matrix,double[]>(ret, costs.ToArray());
}

ロジスティック回帰

線形回帰の最急降下法と似ているが、 ${h = g (\theta ^T x)}$ になる。
${g(z)}$ はシグモイド関数
${ g(z) = \frac{1}{ 1 + e^{-z} } }$
で、解が必ず0~1になる。
コストは ${ y = 0 }$ の場合と ${ y = 1 }$ の場合で場合わけして求める。

public static Tuple<Matrix, double[]> GetResult(Matrix x, Matrix y, double alpha, double lambda, int maxItr, bool logCosts=true)
{
    var theta = new Matrix(x.ColNum, 1);
    var costs = new List<double>();
    while (maxItr > 0)
    {
        var h = GetSigmoid(x * theta);
        var diff = h - y;

        if (logCosts)
        {
            var wktheta = theta.Clone();
            wktheta[0, 0] = 0;
            var j = ( 1.0 / (double)x.RowNum ) *
                    ( ( -1.0 * (y.Transpose() * h.Log())) - (1.0 - y).Transpose() * (1.0 - h).Log())
                    + ( lambda / ( 2.0 * (double)x.RowNum) ) * (wktheta^2).Sum();
            costs.Add(j[0, 0]);
        }

        var grad = ((alpha / (double)x.RowNum) * (diff.Transpose() * x).Transpose()) + ((alpha * lambda) / (double)x.RowNum) * theta;
        grad[0, 0] -= ((alpha * lambda) / (double)x.RowNum) * theta[0, 0];

        theta -= grad;

        maxItr--;
    }

    return new Tuple<Matrix, double[]>(theta, costs.ToArray());
}

static Matrix GetSigmoid(Matrix z)
{
    var a = Math.E ^ z;
    var b = 1.0 / a;
    var c = 1.0 + b;
    var d = 1.0 / c;

    var ret = 1.0 / ( 1.0 + (1.0 / (Math.E ^ ( z ))) );
    return ret;
}

検証
f:id:yambe2002:20160303041844p:plain
コースと同じデータを使って検証した。
Octaveの結果とほぼ一致。いずれも、それらしく動いているようだ。

メモ：当初、あるテストケースで最急降下法がうまく収束しなくて悩んでいたが、標準化を忘れていたのが原因だった。

2016-01-28

二分探索の練習

C# プログラミングコンテスト

二分探索の練習。
問題は「プログラミングコンテスト」「二分探索」あたりで適当にググったもの。

- UpperBoundとLowerBound
CのSTD関数をC#で実装。

public static int LowerBound(int[] ar, int val)
{
    var lb = -1;
    var ub = ar.Length;

    while (ub - lb > 1)
    {
        var mid = (lb + ub) / 2;
        if (ar[mid] >= val)
        {
            ub = mid;
        }
        else
        {
            lb = mid;
        }
    }

    return ub;
}

public static int UpperBound(int[] ar, int val)
{
    var lb = -1;
    var ub = ar.Length;

    while (ub - lb > 1)
    {
        var mid = (lb + ub) / 2;
        if (ar[mid] <= val)
        {
            lb = mid;
        }
        else
        {
            ub = mid;
        }
    }

    return lb + 1;
}

- Topcoder SRM676 Div1Easy(WaterTank)
https://community.topcoder.com/stat?c=problem_statement&pm=14019

Cリットル入る水がめと、配列tとxが与えられる。水がめには、t[i]の時間、x[i]リットルの水が注ぎこまれる。水がめが溢れないようになるべく小さい穴をあけたい。
この穴の単位時間あたりの水の排出量を求める。
配列の大きさは1～50、tとxは1～1,000,000、Cは1～10^9

LowerBoundと同じように二分探索する。
最初、1sごとにループするコードを書いてタイムアウトになってしまった。水の増え方が線形なので、t[i]単位でまとめて計算して問題ない。

public class WaterTank {
    double EPS = 1e-8;

    int[] t;
    int[] x;
    double C;

    bool NoOverflow(double rate)
    {
        var tank = 0.0;
        for (int n = 0; n < t.Length; n++)
        {
            //これだと間に合わない！
            //線形に増えるので1sごとにループする必要はない
            //var water = (double)x[n];
            //for (int sec = 0; sec < t[n]; sec++)
            //{
            //    tank += water;
            //    tank -= rate;
            //    if (tank > C) return false;
            //    if (tank < 0) tank = 0.0;
            //}

            tank += (double)x[n] * t[n];    //キャストを忘れない！（忘れてもサンプルテストは通る）
            tank -= rate * t[n];
            if (tank > C) return false;
            if (tank < 0) tank = 0.0;
        }

        return true;
    }

    public double minOutputRate(int[] pt, int[] px, int pC)
    {
        t = pt;
        x = px;
        C = (double)pC;

        double lb = 0.0, ub = 1000001.0;

        while (ub - lb > EPS)
        {
            var mid = (lb + ub) / 2;
            if (NoOverflow(mid))
                ub = mid;
            else
                lb = mid;
        }

        return ub;
    }
}

- Topcoder SRM456 Div2Hard(CutSticks)
https://community.topcoder.com/stat?c=problem_statement&pm=10563&rd=13909

いろいろな長さの棒sticks[]がある。そして最大でC回、棒を取り出して2つに割ることができる（割ってできた棒をさらに割ることもできる）。
最終的にできた棒を、増加しない順にならべて、K番目の棒を取り出す。このとき、その取り出した棒の最大長を求める。
なお、最低でもK本の棒ができるように、割っていく必要がある。
sticks[]の長さ：1～50
sticks[]の各要素：1～1000000000
C：1～1000000000
K:1～sticks[]の長さ+C

「ある長さにできるかどうか」を判定しながら二分探索する。
終了条件を上限と下限の差にすると無限ループになってしまうので、十分大きな固定回数のループにする。

public class CutSticks {
    double EPS = 1e-10;

    bool CanMake(double len, double[] sticks, int C, int K)
    {
        var cnt = 0;
        foreach (var stick in sticks)
        {
            if (stick - len > EPS)
            {
                var maxDevideNum = (int)(stick / len) - 1;
                var devideNum = Math.Min(C, maxDevideNum);

                cnt += (devideNum + 1);
                C -= devideNum;

                if (cnt >= K) return true;
            }
        }

        return false;
    }

    public double maxKth(int[] sticks, int C, int K) {
        double lb = 0;
        double ub = 1000000001;

        var dsticks = sticks.Select(a => (double)a).OrderByDescending(b => b).ToArray();

        //while(ub - lb > EPS)  //これだと (new int[] { 1000000000 }, 1000000000, 1) のケースで終了しない
        for (int i = 0; i < 500; i++)
        {
            var mid = (lb + ub) / 2;

            if (CanMake(mid, dsticks, C, K))
                lb = mid;
            else
                ub = mid;
        }

        return lb;
    }
}

- Topcoder SRM636 Div2Hard(ChocolateDividingHard)
https://community.topcoder.com/stat?c=problem_statement&pm=13388&rd=16079

チョコレートの板String[] chocolateがある。
chocolateの各char要素は、そのセルのクオリティを表している('0'～'9')。
縦に４つ、横に４つに割ったとき、１つの塊のクオリティの合計の最小値を最大にしたい。この値を答えよ。
chocolate[]の縦横の要素数は4～75

最大のクオリティを求める二分探索を行う。
その中で縦（または横）全パターンのループを回し、パターンごとにクオリティを満足する最小の横Indexを二分探索で求めていく。
横に4分割できればOK。

二分探索の中で二分探索をしていて、はまるとデバッグが大変になる問題だと思う。

public class ChocolateDividingHard {
    int[,] sums;
    int rowNum;
    int colNum;
    int sum;

    //row1 < row2, col1 < col2
    //(row1 + 1, col1 + 1) to (row2, col2)
    int GetArea(int row1, int row2, int col1, int col2)
    {
        var ret = sums[row2, col2] -
            (col1 >= 0 ? sums[row2, col1] : 0) -
            (row1 >= 0 ? sums[row1, col2] : 0) +
            (col1 >= 0 && row1 >= 0 ? sums[row1, col1] : 0);

        return ret;
    }

    bool CanDevide(int col1, int col2, int size, int[] rows)
    {
        for (int rowIdx = 0; rowIdx < rows.Length; rowIdx++)
        {
            var row1 = rowIdx > 0 ? rows[rowIdx - 1] : -1;
            var row2 = rows[rowIdx];
            if (GetArea(row1, row2, col1, col2) < size) return false;
        }
        return true;
    }

    int GetMinColIdx(int startIdx, int size, int[] rows)
    {
        //binary search (l, u]
        int l = startIdx - 1;   //ここを l = startIdx にすると、startIdxが答えの時に間違った解を返す！
        int u = colNum;

        while (l < u - 1)
        {
            var mid = (l + u) / 2;
            if (CanDevide(startIdx - 1, mid, size, rows))
                u = mid;
            else
                l = mid;
        }

        return (u == colNum) ? -1 : u;
    }

    bool CanDevide(int size, int[] rows)
    {
        var col = -1;
        for (int i = 0; i < 4; i++)
        {
            col = GetMinColIdx(col + 1, size, rows);
            if (col == -1) return false;
        }
        return true;
    }

    bool CanDevide(int size)
    {
        for (int row1 = 0; row1 < rowNum; row1++)
            for (int row2 = row1 + 1; row2 < rowNum; row2++)
                for (int row3 = row2 + 1; row3 < rowNum; row3++)
                    if (CanDevide(size, new int[] { row1, row2, row3, rowNum - 1 }))
                        return true;
        return false;
    }

    int findBest()
    {
        //binary search [l, u)
        int l = 0;
        int u = sum + 1;

        while (u > l + 1)
        {
            var mid = (l + u) / 2;

            if (CanDevide(mid))
                l = mid;
            else
                u = mid;
        }

        return l;
    }

    public int findBest(string[] chocolate) {
        rowNum = chocolate.Length;
        colNum = chocolate[0].Length;

        sums = new int[rowNum, colNum];
        for (int row = 0; row < rowNum; row++)
        {
            for (int col = 0; col < colNum; col++)
            {
                sums[row, col] = 
                    (row > 0 ? sums[row - 1, col] : 0) +
                    (col > 0 ? sums[row, col - 1] : 0) +
                    ((int)(chocolate[row][col] - '0')) -
                    (row > 0 && col > 0 ? sums[row - 1, col - 1] : 0);
                sum += (int)(chocolate[row][col] - '0');
            }
        }

        return findBest();
    }
}

2016-01-10

初等整数論の基本の復習

数学 C#

2016/8/20誤記修正

基本的な数学がいまいち使いこなせてないので復習しておく。

参考にしたのは蟻本と以下のサイト。
http://www2.cc.niigata-u.ac.jp/~takeuchi/tbasic/BackGround/

ユークリッドの互除法

整数aとbについて、a > b > 0 かつ a % b を r としたとき、 ${ GCD(a,\,b) = GCD(b,\,r) }$ が成り立つ。

    static int Gcd(int a, int b)
    {
        if (b == 0) return a;
        var p = a > b ? a : b;
        return Gcd(b, p % b);
    }

毎回、半分以下になっていくので計算量は ${ O\left ( log \, max (a, b) \right ) }$ 程度になる。

拡張ユークリッドの互除法

${ ax + by = GCD(a,b) }$ の整数解 ${(x, y)}$ は常に存在し、計算することができる。
この方程式の解を求める関数を int ExtGcd(int a, int b, ref int x, ref int y)とすると、次のように求められる。返り値はGCD(a, b)。

    static int ExtGcd(int a, int b, ref int x, ref int y)
    {
        int d = a;
        if (b != 0)
        {
            d = ExtGcd(b, a % b, ref y, ref x);
            y -= (a / b) * x;
        }
        else
        {
            x = 1;
            y = 0;
        }
        return d;
    }

手でやってみると、なぜこうなるかよく分かる。計算量はGcd()と同じくらい。

エラトステネスの篩

整数iが素数かどうか判定する。また、j番目の素数を求める。これは簡単。

    static void Sieve(int[] prime, bool[] isPrime)
    {
        for (int i = 0; i < prime.Length; i++) prime[i] = -1;
        for (int i = 0; i < isPrime.Length; i++) isPrime[i] = true;
        isPrime[0] = isPrime[1] = false;

        var idx = 0;
        for (int i = 2; i < isPrime.Length; i++)
        {
            if (isPrime[i])
            {
                prime[++idx] = i;
                for (int j = 2 * i; j < isPrime.Length; j += i) isPrime[j] = false;
            }
        }
    }

蟻本によると、計算量は ${ O(n \, log \, log \, n) }$ 程度で、通常の範囲ならほぼ線形と考えてよいらしい。

整数の合同

整数aとbのmで割った余りが等しい（a-bがmで割り切れる）とき、
${ a \equiv b \, (mod\: m) }$
と表し、aとbはmを法にして合同、という。

線形の合同には次のような性質がある。
${ a \equiv a \, (mod\: m) }$
${ a \equiv b \, (mod\: m) }$ 　ならば　 ${ b \equiv a \, (mod\: m) }$
${ a \equiv b \, (mod\: m) }$ 　かつ　 ${ b \equiv c \, (mod\: m) }$ 　ならば　 ${ a \equiv c \, (mod\: m) }$

${ a \equiv c \, (mod\: m) }$ 、 ${ b \equiv d \, (mod\: m) }$ 　とすると、
${ a + b \equiv c + d (mod\: m) }$
${ a - b \equiv c - d (mod\: m) }$
${ a \times b \equiv c \times d (mod\: m) }$
が成り立つ（除算は成り立たない）。

${ a \times c \equiv b \times c(mod\: m) }$ の場合、
${ a \equiv b(mod\:\, m\, /\, gcd(m,\, c)) }$
が成り立つ。

べき乗法

以下の単純なやり方で普通は十分。

    static ulong ModPow(ulong x, ulong n, ulong mod)
    {
        ulong ret = 1;
        while (n > 0)
        {
            if ((n & 1) == 1) ret = ret * x % mod;
            x = x * x % mod;
            n >>= 1;
        }
        return ret;
    }

計算量は ${ O(log\, n) }$ 。
modが大きいとオーバーフローする（途中でmod^2程度の数が表れる）ので、改善したアルゴリズムがここに紹介されていた。
http://www2.cc.niigata-u.ac.jp/~takeuchi/tbasic/BackGround/power.html

中国の剰余定理

${ m_{1}, m_{2} }$ が互いに素な整数なら、任意に与えられる整数 ${a_1, a_2}$ に対し
${ x \equiv a_{1} (mod\: m_{1}) }$
${ x \equiv a_{2} (mod\: m_{2}) }$
を満たす整数xが、 ${ m_{1}m_{2} }$ を法にしてただ一つ存在する。

これを拡張して、整数 ${ m_{1}, m_{2}, ...m_k }$ がどの二つも互いに素なら、任意に与えられる整数 ${a_1, a_2, ... a_k}$ に対し
${ x \equiv a_{1} (mod\: m_{1}) }$
${ x \equiv a_{2} (mod\: m_{2}) }$
...
${ x \equiv a_{k} (mod\: m_{k}) }$
を満たす整数xが、 ${ m_{1}m_{2}...m_k }$ を法にしてただ一つ存在する。

これを利用して、TopcoderのSRM657(Div2Hard)を解いてみる。

整数a、b、cが与えられたとき、 ${ ax^{2} + bx + c \equiv 0\: mod\: 10^{9} }$ を満たすxを求める。
（a、b、cは0～999,999,999）

xの範囲が0,...999999999と大きいので、単純にループさせては時間がかかりすぎる。

${10^{9}}$ は ${ 2^{9} }\times{ 5^{9} }$ かつ ${ 2^{9} }$ と ${ 5^{9} }$ は互いに素なので、中国の剰余定理を利用して

${ ax^{2} + bx + c \equiv 0\: mod\: 2^{9} }$
${ ax^{2} + bx + c \equiv 0\: mod\: 5^{9} }$

を満たすxを求めればよい。

それぞれの式から、

${ x \equiv x_{2}\: mod \: 2^{9} }$
${ x \equiv x_{5}\: mod \: 5^{9} }$

が導き出されるので、片方を満たす最小非負代表を求めてから、もう片方の式を満足できるまで法の数だけ足していけばよい。
下の例では、 ${ x \equiv x_{5}\: mod \: 5^{9} }$ を満たす最小値を算出し、 ${ x \equiv x_{2}\: mod \: 2^{9} }$ が満たされるまで ${ 5^{9} }$ を足していっている。

public class PolynomialRemainder {
    int POW2 = 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2;
    int POW5 = 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5;

    ulong GetAx2BxC(ulong mod, ulong x, ulong a, ulong b, ulong c)
    {
        return ((a * (x * x % mod) % mod) + (b * x) % mod + c) % mod;
    }

    int GetRoot(int mod, int a, int b, int c)
    {
        for (int x = 0; x < mod; x++)
        {
            var ret = GetAx2BxC((ulong)mod, (ulong)x, (ulong)a, (ulong)b, (ulong)c);
            if (ret % (ulong)mod == 0)
            {
                return x;
            }
        }

        return -1;
    }

    public int findRoot(int a, int b, int c)
    {
        var x2 = GetRoot(POW2, a, b, c);
        var x5 = GetRoot(POW5, a, b, c);

        if (x2 == -1 || x5 == -1) return -1;

        int x = x5;
        while (x % POW2 != x2)
            x += POW5;

        return x;
    }
}