2021-03-25

平衡二分探索木を使った2-optの効率化について

こんにちは。平衡二分探索木で2-optを高速化する方法について研究したので、本記事ではそれについて書こうと思います。

2-optとは

2-optは巡回セールスマン問題(TSP)の近似解を求めるヒューリスティックアルゴリズムです。
下図のように、サイクルに対してランダムに二つのパスを選び、パスを組み替えた方がサイクルの長さが短くなる場合、パスを組み替えます。具体的には、A→B + C→D > A→C + B→D の時、辺ABと辺CDを辺ACと辺BDに組み替えます。
f:id:aspirator:20210325193100j:plain しかし、辺を張り直すだけではサイクルにならないので、B~CまたはD~Aの間の順序を逆にする必要があります。
2-opt法では順路を配列などで管理することでこれを実現しています。
f:id:aspirator:20210325193138j:plain そしてこれが`測定で用いた2-optの実装です。

#include<iostream>
#include<algorithm>
#include<vector>
#include<iterator>
#include<random>
#include<chrono>
using namespace std;
using time_point = chrono::system_clock::time_point;


time_point get_time() {
    return chrono::system_clock::now();
}
int time_diff(chrono::system_clock::time_point start) {
    return static_cast<int>(chrono::duration_cast<chrono::milliseconds>(chrono::system_clock::now() - start).count());
}

std::random_device seed;
std::mt19937 random_number(seed());

struct point {
    int x, y;
};
vector<point>points;
double dist(point a, point b) { //点a, b間の距離
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
int loop_count = 0;
void two_opt(vector<int>& route, int limit_ms) {
    time_point start_time = get_time();
    int a, b, size = route.size();
    while (time_diff(start_time) < limit_ms) {
        loop_count++;
        a = random_number() % size;
        b = random_number() % size;
        if (a == b)continue;
        if (b < a) swap(a, b);
        if (dist(points[route[a]], points[route[(a + 1) % size]]) + dist(points[route[b]], points[route[(b + 1) % size]]) >
            dist(points[route[a]], points[route[b]]) + dist(points[route[(a + 1) % size]], points[route[(b + 1) % size]])) {
            reverse(route.begin() + a + 1, route.begin() + b + 1);
        }
    }
}

int main() {
    const int N = 10000;
    points.resize(N);
    std::vector<int>route(N);
    for (int i = 0; i < N; i++) {
        points[i] = { random_number() % 10000,random_number() % 10000 };
        route[i] = i;
    }
    two_opt(route, 2000);
    cout << loop_count << endl;
}

しかし、順序を反転させる操作には $O(N)$ かかってしまいます。
この記事では平衡二分探索木を用いることで反転操作を $O(\log N)$ にします。
今回は平衡二分探索木としてtreapを使いました。
長いので別に載せておきます。
treapのライブラリ
そして測定で用いたコードがこちらです。

#include <iostream>
#include <vector>
#include <utility>
#include <cmath>
#include <random>
#include <chrono>
#include "reversibleArray.cpp"
using namespace std;
using time_point = chrono::system_clock::time_point;


time_point get_time() {
    return chrono::system_clock::now();
}
int time_diff(chrono::system_clock::time_point start) {
    return static_cast<int>(chrono::duration_cast<chrono::milliseconds>(chrono::system_clock::now() - start).count());
}

std::random_device seed;
std::mt19937 random_number(seed());

struct point {
    int x, y;
};
vector<point>points;
double dist(point a, point b) { //点a, b間の距離
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
int loop_count = 0;
void original_two_opt(vector<int>& route, int limit_ms) {
    time_point start_time = get_time();
    int l, r, size = route.size();
    reversible_array opt(route);
    while (time_diff(start_time) < limit_ms) { 
        loop_count++;
        l = random_number() % size;
        r = random_number() % size;
        if (l == r)continue;
        if (r < l) swap(l, r);
        int c = opt[l], d = opt[(l + 1) % size],
            e = opt[r], f = opt[(r + 1) % size];
        if (dist(points[c], points[d]) + dist(points[e], points[f]) >
            dist(points[c], points[e]) + dist(points[d], points[f])) {/*組み替えた方が短くなる場合反転*/
            opt.reverse(l, r);
        }
    }
    for (int i = 0; i < size; i++) { /*平衡二分探索木で出た結果をrouteに戻す*/
        route[i] = opt[i];
    }
}

int main() {
    const int N=10000;
    points.resize(N);
    std::vector<int>route(N);
    for (int i = 0; i < N; i++) {
        points[i] = { random_number() % 10000,random_number() % 10000 };
        route[i] = i;
    }
    original_two_opt(route, 2000);
    cout << loop_count << endl;
}

計測

計測方法

計測はwhileループを2000ms回してループした回数を数えます。
それを5回繰り返して平均値をとります。

計測結果

頂点数とループ回数

頂点数	$10^{2}$	$10^{2.5}$	$10^{3}$	$10^{3.5}$	$10^{4}$	$10^{4.5}$	$10^{5}$	$10^{5.5}$	$10^{6}$
2-opt	$1.07\cdot 10^{6}$	$1.20\cdot 10^{6}$	$1.17\cdot 10^{6}$	$1.13\cdot 10^{6}$	$9.56\cdot 10^{5}$	$4.67\cdot 10^{5}$	$9.68\cdot 10^{4}$	$3.04\cdot 10^{4}$	$9.84\cdot 10^{3}$
オリジナル	$1.69\cdot 10^{5}$	$1.37\cdot 10^{5}$	$1.12\cdot 10^{5}$	$9.29\cdot 10^{4}$	$7.61\cdot 10^{4}$	$6.41\cdot 10^{4}$	$5.80\cdot 10^{4}$	$5.14\cdot 10^{4}$	$4.48\cdot 10^{4}$

f:id:aspirator:20210830120559j:plain

考察

通常の2-optでは頂点数が $10^{2}$ ~ $10^{4}$ の範囲でループ回数が $1.1\cdot 10^{6}$ に固まっています。
これはある程度解が改善するとそれ以上改善されにくくなり、ボトルネックであるrevesrseパートが行われにくくなるからだと考えられます。
また、頂点数が $10^{5}$ 未満の時、通常の2-optにかなわなかったのも、同様の理由から取得 $O(\log N)$ の平衡二分探索木より取得 $O(1)$ のvectorの方が高速になったからだと考えられます。
結果の表やグラフから、頂点数が多い場合やタイムリミットが短い場合は平衡二分探索木を用いた2-optを使った方が効率が良く、頂点数が少ない場合やタイムリミットが長い場合は普通の2-optを使う方が効率が良いことが分かりました。

途中まで平衡二分探索木を用いて、ある程度解が良くなったところで通常の配列を用いた2-optを使うと、さらに効率が良くなるかもしれません。

感想

自分からこんな研究をするなんて自分でもびっくりしました。
マラソンコンテストへの採用は出来なさそうだったので、またいつか細かいところの定数倍高速化をしたり、続きの研究をしてリベンジしたいと思います。

夏休みの自由研究の課題としてこれを提出しました。

2021-02-21

自己参照構造体でセグ木を書く

これを作るに至った動機

平衡二分探索木を実装したい！
でもいきなり作るにはハードルが高い...
そうだ！自己参照構造体でセグ木を書いて練習しよう！
自己参照構造体で書くセグ木は、他にも発展性が高いらしいので一石二鳥！
（これが書けるようになると、必要なところだけ作るセグ木や、永続セグ木などをすっとかけるようになると赤コーダーの人がツイートしていた）
※始めに断っておくとこのセグ木は再帰で書いているため遅く、あまり実用的ではありません。

そもそも自己参照構造体とは

メンバに自分自身と同じ型の構造体へのポインタを持つ構造体のこと

struct hoge{
    hoge* next; //二分木だとポインタが二つになる
};

めっちゃ簡略化するとこんな感じ
std::listや大体の二分探索木がこれ

実装

C++での実装を乗せます。
折角なので軽く解説します。

constexpr long long LINF = 4611686015206162431;
int op(int x, int y) {
    return min(x, y);
}
int e() {
    return LINF;
}
namespace Library {
    template<class T>class slow_segment_tree {
        struct node {
            T sum = e();
            node* child[2];
        };
        node* root = NULL; //根ノード
        int depth = 0; //セグ木の深さ
        void make_node(node*& np, int d) {
            np = new node;
            if (d > 0) {
                make_node(np->child[0], d - 1);
                make_node(np->child[1], d - 1);
            }
            return;
        }
        void set_all(node*& np, int pos, int d, vector<int>& v) {
            np = new node;
            if (d == -1) {
                if (pos < v.size())np->sum = v[pos];
                return;
            }
            set_all(np->child[0], pos, d - 1, v);
            set_all(np->child[1], pos + (1 << d), d - 1, v);
            np->sum = op(np->child[0]->sum, np->child[1]->sum);
            return;
        }
        void set_sub(node* np, int pos, T n, int p) {//ノード 変更する位置 変更後の数値 深さ
            if (p == -1) {
                np->sum = n;
                return;
            }
            set_sub(np->child[(bool)(pos & (1 << p))], pos, n, p - 1);
            np->sum = op(np->child[0]->sum, np->child[1]->sum);
        }
        T prod_sub(node* np, int l, int r, int a, int b) {//今いるノード, ノードの範囲[l, r), 欲しい範囲[a, b)
            if (r <= a || b <= l) {
                return e();
            }else if (a <= l && r <= b) {
                return np->sum;
            }else {
                return op(prod_sub(np->child[0], l, (l + r) / 2, a, b), prod_sub(np->child[1], (l + r) / 2, r, a, b));
            }
        }
    public:
        slow_segment_tree(int size) {
            assert(size > 0);
            while ((1U << depth) < (unsigned int)(size)) depth++;
            make_node(root, depth);
        }
        slow_segment_tree(vector<int>& v) {
            int size = v.size();
            assert(size > 0);
            while ((1U << depth) < (unsigned int)(size)) depth++;
            set_all(root, 0, depth - 1, v);
        }
        void set(int pos, T n) {//pos番目をnにする
            assert(pos < (1 << depth));
            set_sub(root, pos, n, depth - 1);
        }
        T get(int pos) {//pos番目を返す
            assert(pos < (1 << depth));
            node* n = root;
            for (int d = depth - 1; d >= 0; d--) {
                n = n->child((bool)(pos & (1 << d)));
            }
            return n->sum;
        }
        T all() {
            return root->sum;
        }
        T prod(int l, int r) {// [l, r) の和
            assert(l < (1 << depth) && r <= (1 << depth));
            if (r <= l) {
                return e();
            }
            return prod_sub(root, 0, 1 << depth, l, r);
        }
    };
}

抽象化セグ木なので二項演算(op)と単位元(e)を変えることで、様々な機能を持たせることができます。
上ではrange max queryです。

コンストラクタ

コンストラクタは、深さを持たせてdfsしながらノードを作っていきます。
rootをNULLで初期化することで簡潔に実装出来ます。

値の取得

深さをfor文で回し、npを更新します。
左右どちらのノードに進むか判定するためにbit演算を使っています。
これは非再帰で書くことができました。

値の更新

深さを持たせ再帰します。
左右どちらのノードに進むか判定するためにbit演算を使っています。
葉まで行ったら値を書き換え、帰りがけに更新します。
非再帰で書くのは厳しそうでした。

最小値の取得

親から子へとたどっていき、区間に収まるものを見つけます。これは普通のトップダウン型セグ木とあまり変わりません。
またこれもボトムアップで非再帰にするのは厳しそうでした。

速さ

AOJのRMQの問題でACLのセグ木と比較しました。
ACLのセグ木⇒0.02s
自己参照構造体で作るセグ木⇒0.04s
このセグ木はどうしても再帰で書くことが多くなるため、速さを出すことができず、速いとは言えません。

2021-02-15

JOI2021 参加記

初めに

みんな書いてたので参加記を書こうと思いました。

こういうのを書くのは初めてなので色々おかしくても多めに見てください。

予選

A問題

A - 往復すごろく (Round Sugoroku)

setで殴ることができ、 $O(NlogN)$ 。

予選A問題にしては割と難しかったように思いました。

B問題

まだB問題なのにパッと方針が出ずかなり焦りました。

3進法にパンケーキの状態を3進法に直し、BFSで前計算して $3^{N}N+Q$ という満点解法が思いつき、小課題と実装量が変わらなそうだったのでいきなり満点解法の実装を始めました。

しかし、実行したところバグったためBはいったん飛ばしてCを解くことにしました。

Cを解いた後

C問題

見た瞬間にダイクストラ解を思いつき、そのまま実装しました。

多少バグらせて時間を掛けるも満点を得ることができました。

B問題へ戻る

一回D問題に目を通すも解法が思いつかなかったためB問題に戻りました。

競技時間終了までデバッグするも部分点すら取れず競技時間終了。

競技終了してから今更になって部分点を取っていなかったことを後悔しました。

終了後

kichiさんによるスプレッドシートを見ると、僕と同様に200点を取っている人は多く、本選進出はかなり厳しいという印象を受けました。

しかし、ふたを開けてみるとボーダーは186点だったようで、本選進出することができました！yay!

また、双子によれば、200点は高確率で通ると予想していたようで、経験者の勘の強さに驚きました。

本選

開会式

動画のクオリティが凄かったです。（小並感）

自己紹介動画

twitterでしか知らない人の顔を見ることができ、おもしろかったです。

またdefineさんに絶起する呪いを掛けられました。

本選競技

前日は11時半に就寝し、無事絶起することなく7時半に起床。

前日の地震の影響で協議開始時刻が10時に延期されたため、開始までライブラリの整備などをしていました。

A問題

家庭菜園が趣味のビ太郎は自宅の庭でビバ草という植物を育てている.

庭には N 株のビバ草が東西方向に一列に植えられており，西側から順に $1$ から $N$ までの番号が付いている.

現在，ビバ草 $i (1 ≦ i ≦ N)$ の背丈は Ai である.
育てているビバ草は特別な品種改良の結果，水を与えるたびに背丈が $1$ 伸びる.

ビ太郎は庭の見栄えを良
くするために水やりを複数回行い，以下の条件を満たすようにしたいと考えている．

すべての水やりを行った後のビバ草 $i (1 ≦ i ≦ N)$ の背丈を $Bi$ とする. このとき，「 $1 \leqq j \leqq k-1$ に対しB_j < B_j+1]」かつ「 $k \leqq N-1$ に対し B_j > B_j+1」を満たすような整数 $k(1 \leqq k \leqq N)$ が存在する.

ただし，ビ太郎は不器用なため，1 回の水やりにおいて，ある区間上のビバ草に一斉に水を与えることしかできない.

すなわち，水やりを行うたびにある整数 $L, R (1 ≦ L ≦ R ≦ N)$ を選び，ビバ草 $L, L + 1, . . . , R$ に水を与える.
ビ太郎は水やりの回数をできるだけ少なくしたい．
ビバ草の数と現在の背丈の情報が与えられたとき，条件を満たすのに必要な水やりの回数の最小値を求めるプログラムを作成せよ.

制約

$2 ≦ N ≦ 200 000$ .

$1 ≦ Ai ≦ 1 000 000 000 (1 ≦ i ≦ N)$ .

左右からシミュレーションしていくのが一番最初に思いつきました。

しかし、やや実装が面倒くさいと思ったため、「端からi番目まで山なりにするのにかかる水やりの回数」を左右から累積和のように取っていき、水やりの回数を二分探索しました。

しかし、日ごろ二分探索になれていなかったせいか二分探索をバグらせまくって1時間半かけても0点しか取れませんでした。

その後、A問題を飛ばし、Bを解きに行きます。

B問題

JOI 平原は東西方向に広がるとても大きな平原である．この平原は数直線と見なすことができ，各地点は
東向きを正とする座標であらわされる．JOI 平原は冬を迎え $N$ 個の雪玉が異なる座標に作られた．雪玉に
は西から順に $1$ から $N$ までの番号が付いている．最初，雪玉 $i (1 ≦ i ≦ N)$ の座標は整数 $X_i$ であった．
また JOI 平原は冬に強い風が吹くことで知られている．あなたは $Q$ 日間の風の観測データを入手した．

$j$ 日目 $(1 ≦ j ≦ Q)$ の風のデータは整数 $W_j$ であらわされる．

$W_j$ が負のときは西向きに， $W_j$ が負でないときは東向きに，強さ $| W_{j} |$ の風が吹いたことを意味する．
風が吹くと，雪玉は風と同じ向きに，風の強さと同じ距離だけ転がる．

すなわち $j$ 日目 $(1 ≦ j ≦ Q)$ の始
めに座標 $x$ に雪玉があったとき，その雪玉は座標 $x$ から座標 $x + W_j$ まで転がる．

$j$ 日目の終わりには，その雪玉の座標は $x + Wj$ になる．ただし，各日においてすべての雪玉が同時に，同じ速さで転がる．
最初，JOI 平原全体に雪が積もっていた．雪が積もっている範囲を雪玉が転がると，雪が付着し，雪玉の重さが増え，その範囲の雪はなくなる.

すなわち， $a$ を整数とし，座標 $a$ から座標 $a + 1$ までの範囲に雪が残っているとするこのとき，雪玉がこの範囲を転がると，その雪玉の重さが $1$ 増えて，座標 $a$ から座標 $a + 1$ までの範囲の雪がなくなる.

ただし，雪が残っていない範囲を雪玉が転がったとしても，雪玉の重さは変わらない．
最初，すべての雪玉の重さは $0$ であった．

また，観測した $Q$ 日間に新たに雪は降らなかった．
あなたは $Q$ 日目の終わりにおける雪玉の重さを知りたい．
雪玉の最初の座標， $Q$ 日間の風の観測データが与えられたとき， $Q$ 日目の終わりにおける，それぞれの雪
玉の重さを求めるプログラムを作成せよ．

制約

$1 ≦ N ≦ 200 000$ .

$1 ≦ Q ≦ 200 000$ .

$|X_i| ≦ 1 000 000 000 000 (= 10^{12}) (1 ≦ i ≦ N)$ .

X_{i} < X_{i+1} (1 ≦ i ≦ N − 1).

$|W_{j}| ≦ 1 000 000 000 000 (= 10^{12}) (1 ≦ j ≦ Q)$ .

小課題

1.(33点) $N\leqq 2000, Q\leqq 2000$ .