4bitにっき

ぽよ

SuperCon2016優勝


実際に書いたソースコードは最後に載せています。
感覚的で適当な表現が多いです。
誤りがあったら教えて下さい。

らてプロ@LatteMalta と 1年生のnarahara君(今年のJOI本選参加します!)と一緒にチームGhostDivとしてSuperCon2016に出た。
また久留米高専から、KO@sujinbemani_79c, つよっしー@TSUYOSSHI, ゲドリンク@gedorinku のチームKsanteam(久留米㌠)も出ていた。

呪術の甲斐あって夏季セミとSuperConの期間が丸かぶりしてくれたので(強い人などが夏季セミに吸われる)、SuperCon2016でついに優勝を成し遂げた。
ついでに3年前の先輩達を見習って、GPU縛りのレギュレーションを遵守した。(!)

0日目

阪大への移動日。
新幹線内でトランプに入れてもらえず悲しかった。
早速たこ焼きを食べたり、ジュンク堂梅田店で機械学習入門とかゲーム理論入門とか4冊買ってバッグをもらったりした。
ちなみにたこ焼きを自分で焼く店へ行ったけど、結構焼けないとうまくひっくり返せなくて、焼けるのを待っていたら凸凹したお好み焼きができかけた。
調節の仕方が分からないエアコンで激冷えのホテルの部屋のお風呂に浸かったりして、ブルジョワ生活を楽しむ。

1日目

問題が配られるまで魂が抜けていた。

問題概要

==============================================================
C(2~3)色に塗り分けられたN個(512~2048)の頂点集合が与えられるので、良さみの深いD次(4~16)正則グラフを作って、総最短経路長(全ての頂点のペアに対する最短経路長の総和)をなるべく最小化せよ。
ただし、両端の頂点の色がi,jである辺の重みはw_{ij}(1~9)である。
ちなみに正則グラフとは、全ての頂点の次数が等しい無向グラフのことらしい。
==============================================================

競技方法

==============================================================
プログラムを提出して運営側で実行するのではなく、事前に準備した上で2時間の競技時間中に自分たちで与えられた9つの問題への答えをそれぞれ出し、提出して順位を競う。
30秒に1回行える提出をすると、各問題の自分たちの順位は分かるが、他のチームについては何も分からないようになっていた。
最終的な全体の順位は、順位毎に決まった点数を総和して競う。
1位から順に 25-18-15-12-10-8-6-4-2-1 という配点で、1位の点が他よりかなり重いことが分かる。
==============================================================


僕はこの日は考察に徹していた。
このとき考えていた解法は、規則的に作ったグラフに対して、色の割り振りを最適化するというもの。(のちに却下する)
他の有力な解法との消去法で選んだ。
頂点が多い上にグラフの辺も多いので、どんな状態からでも割と良い解へ辿り着けるし、焼きなまし法では十分収束させられないと思ったので、最適化には山登り法を使うことにした。
他に考えた解法として以下の様なものがあった。
ちなみにいずれの解法でもグラフが非連結になってしまうことが気になるが、今回の問題ではD=4の時点で適当にグラフを作って非連結になることはまずあり得なさそうだった。

正則グラフを保ったまま辺の接続先をSwapする山登り法

評価でASPLを使う場合、辺が多すぎるので大してループできない今回の問題ではつらすぎでは?

頂点を1つずつ追加するビームサーチ

サイズNで条件を満たすグラフがサイズN+1の最適化に繋がりそうで良さそうに見えるが、ASPLが重いのでやはり計算が間に合わなさそう

頂点集合を等しい2つの集合に分けて分割統治

統治つらすぎやんけ


一方らて君がGPU並列化を学んでワーシャルフロイド法を書いたら300msくらいになったという報告を聞いて、予想よりひどいのでドン引きする。
GPUパワーをよく知らないので行列積でASPLを求める手法(配布された資料にワーシャルフロイド法などと共に書いてあった)はワーシャルフロイド法より早くできなさそうだし、これが最適だろうということになった。


しかし帰ってから、らてくんが「ダイクストラ法のソート部分を(辺の最大重み+1)個のqueueで代用するマイナーなテク」が使えてO(ND)と定数倍で1頂点からの全最短経路長が求まると言うので、希望が見えた。
(つまりO(N^2D)でASPLが求まる)
別にマイナーでもないらしい

2日目

らて君が上記のアルゴリズムを書くととGPU並列化なしでASPLが250msで求まった。なおGPUを使っても180msとかだった。
この分だとGPU使うよりもCPUの残る11スレッド(全部で12スレッドあるため)を使ったほうが早いのではないかと話し合った。

僕の方は、前述のアルゴリズムを使ったら全然良い結果が出なくて冷えていた。
ASPL以外の早くて有効な評価方法があれば強いのにな~と思って、木の直径を求めるO(N)アルゴリズムがグラフでも適用できないか、せめて近似できないかという話をした。
直径がO(N)で求まれば相当強いが、そんなうまくいくはずはなく、ちゃんと全点間の最短経路長を求めないとかすりもしないケースがあるようだった。

ただ、直径を求めた時に、N=2048,D=4のランダムなケースでも直径が11(確か)くらいしかなくて、つまりグラフがめちゃくちゃ密だと分かった。N=2048,D>=10ともなると直径は4くらいだった。
(人脈を7人たどると地球上のほとんどの人と知り合えるという話を聞いたことがあるが、僕の中での信ぴょう性が増した。)

この日は、こんなに頂点同士が近いなら頂点を何個か乱択して他の頂点への最短経路長を求めれば、ASPLを近似できるのではないかとか考えながらホテルに戻った。
そのような近似は、グラフの形が変わった時にちゃんと評価値が反応してくれると限らないのでやばいという結論に至った。

3日目

僕はまだ考察をしていた。
(心の声)
辺の接続先がSwapされてグラフが変形した時、他の頂点への最短経路長が一番変化しやすいのは変形箇所周辺の頂点なはずで。それに辺がめっちゃ多いので、いい感じに辺が助けあってくれて、頂点毎に最短経路長を改善しても最終的にはASPLがそこそこ下がるんじゃないか。
(心の声終わり)
ということで、辺の接続先をSwapした時に限り、その変形箇所に一番関係する4つの頂点を取ってきて、それらからの各最短経路長の2乗の総和を変更前と後で求めて差を取ることで、ASPLの差分を代用できるのではないかと考えた。
2乗したのは、他のより長めの最短経路長があると1乗程度よりもやばい気がしたから。

ともかく、これなら高速に評価が行えるので明るい未来を得た。
ただし、多分このような手法は辺の接続先がSwapされた時の相対評価としてしか使えないので、正則グラフを保ったまま辺の接続先をSwapする山登り法を行うことに決めた。

そして実装途中でホテルに帰る。初めて印刷をした。

この頃らて君は、並列化の高速化を目論んだり必要そうな関数をせっせと整えたりするのもとっくにやり終えていた。
そこでデバッグなどをかなり手伝ってもらっていた。というかほとんどやってもらった。

3日目の夜~4日目の朝

徹夜する。
帰ってすぐ実装途中の 山登り法+高速な評価 を完成させたところ、Graph Golf(今回の問題の背景となったコンテストらしい。このコンテストのサイトだけは閲覧を許されていた。今回の問題との相違点は、辺の重みが全て1であることくらい。)の一部のケースでほぼ1位、そうでなくとも大体上位みたいな点数が手元のPCで出せたりしたので調子に乗りまくった。
ちょっとした工夫程度のはずの、局所解を回避するために書いた「たまにグラフの辺を山登りとか考えずにSwapする」という部分がかなり効いた。
感覚的な話だけど、このアルゴリズムでの登る山が細かいギザギザだらけという感じだったんだと思う。
この時点でGPUを使う方針は完全に消えた。
並列化が下手なせいかうまみがあまり得られないし、それよりもCPUの眠っているスレッドを使って9つのプログラムを同時実行したほうがおいしかったから。

4日目 午前中

この日の13時から競技が行われたので、午前中しか純粋な作業時間がなかった。
ホテルから持ち込んだソースコードを必死に写して、本番の競技への準備などをした。
具体的には、プログラムを終了しても再起動時に前回の状態から最適化を始められる機能や、局所解回避が発生する確率をプログラムの実行中に変更できる機能を実装した。

4日目 競技中

ここでうっかりと言うにはあまりにひどいミスが発覚した。
深夜テンションでチューニングをしている時に、Graph Golfと同じケースしか試さず、つまり辺の全ての重みが1のケースだったので、本番ではチューニングした山登り法がうまくいかなかった。
ここで慌ててその場しのぎでパラメータを調整してほんの少しは改善する。
よく見ると良い解に辿り着いているのにASPLが突然落ちる時があるみたいなので、ネックにならない程度にたまにASPLを求め、極端に悪くなっていたら元に戻すという改善をしたSolver_.cなるファイルを生成して、確か終了40,50分前くらいからそれを動かした。
らて君がASPLを求める関数を整備していなかったら確実に死亡していた。
周りのチームがガンガン迫ってきていたのでものすごいプレッシャーを感じながら震える手で書いていたのを一番覚えている。


時間がないのでかなり雑な改善だった割に急に解がよくなり始めるのだが、書いた山登り法に、序盤の解の改善がゆっくりすぎるという弱点があったため、相変わらず他のチームに攻め立てられて怖かった。
今考えると競技中にせめて初期解の工夫とかを加えればよかったのだけど...。
この時はパラメータの調整で序盤の収束の遅さを本当に雀の涙程度改善した。

あと、6問目の順位があまりにひどかったのでプログラムの定数を調整して6Solver_.cを作った。

確か終了30分前くらいに、7問目で突然1位を剥奪され、左後ろのほうで仲の良い明石高専のチームが喜んでいるのが聞こえたので、慌てて7Solver_.cなるものを生成した。
7問目のケースが、辺の重みに1と9が入り乱れる超苦手ケースだったのだが、他のチームもさすがにこの極端なケースには苦戦しているみたいなので、このケースで点数の重い1位を取り返そうと考えていた。
ちなみに7Solver_.cは、ASPLが下がることをほとんど許さないように定数を調整したプログラムだった。
7Solver_.cのテストをしている間既存のプログラムのパラメータ調整でその場しのぎをして、明石と1位を取ったり取り返されたりしていた。

最終的に7問目は1位が取り返せそうなスコアが出たので、1位を保持しているケースのスコアを最後に念のため更新しつつ、時間ギリギリで7問目の答えを提出した。

5日目 結果発表

1問目から順に1,4,5,2,1,5,1,1,1という順位を取って、優勝できた。

感想

作業時間が3日間半、しかも1日に10時間しかスパコンが使えない、しかもC言語のみというコンテストなので、色々試しているうちに時間がほとんどなくなってしまい、無骨なソースコードしか出来上がらなかったのは後悔している。
しかし、きちんとインターフェイスを決めていたし、YellowCoderらて君の能力が異常に高かったので、僕はらて君の作った関数を呼ぶだけで細かい実装は考えずに最適化アルゴリズムの部分に集中することができて、仕事の分担がうまく出来たと思う。
デバッグもらて君にかなりやってもらったので、本当にありがとう。
ひどいミスなどもしてしまったが、結果オーライということで...。
あと懇談会で寿司がなかったのは罪深い。

ソースコード

僕の書いたSolver_.cと、らて君の書いた関数群を記事の末尾に貼っておきます。

最終日の写経でコメントを書き写すまでの気力がなかったため、コメントが消えてしまっています。
でも最終的に使ったソースコードだけなのででそんなに長くなく、読みやすいと思います。
また、Global.cという擬似乱数の関数や問題の入力などを担当する部分もありましたが、持ち帰り忘れました。
RandomManager.cという乱択を補助する関数群も(最終的には全く呼び出していないが)残っていたので一応貼っておきます。

GraphEvaluator.c(らてくん)

#ifndef _GLOBAL
#define _GLOBAL
#include "./Global.c"
#endif

#ifndef _GRAPHGENERATOR
#define _GRAPHGENERATOR
#include "./GraphGenerator.c"
#endif

#define INF (1001001001)

void GetDistance(int s,int N,int graph[MAX_N][MAX_D],int color[MAX_N],int dist[MAX_N]){
	for(int i=0;i<N;i++)dist[i]=INF;
	dist[s]=0;

	int queue[10][MAX_N],queue_size[10]={};
	queue_size[0]=1;
	queue[0][0]=s;
	
	int r=0;
	for(int i=0;i<=r;i++){
		int m=i%10;
		for(int j=0;j<queue_size[m];j++){
			int v=queue[m][j];
			if(dist[v]<i)continue;
			for(int k=0;k<D;k++){
				int u=graph[v][k];
				int d=dist[v]+W[color[v]][color[u]];
				if(dist[u]>d){
					dist[u]=d;
					queue[d%10][queue_size[d%10]++]=u;
					if(r<d)r=d;
				}
			}
		}
		queue_size[m]=0;
	}
}

void GetAllDistance(int N,int graph[MAX_N][MAX_D], int color[MAX_N], int dist[MAX_N][MAX_N]){

	for(int i=0;i<N;i++)GetDistance(i,N,graph,color,dist[i]);
}

int  ASPL(int N,int graph[MAX_N][MAX_D],int color[MAX_N]){
	int aspl=0;
	int dist[MAX_N];
	for(int i=0;i<N;i++){
		GetDistance(i,N,graph,color,dist);
		for(int j=0;j<N;j++)aspl+=dist[j];
	}
	return aspl/2;
}

//Union Find Tree
int uf_par[MAX_N];
int uf_find(int x){
	return x==uf_par[x]?x:(uf_par[x]=uf_find(uf_par[x]));
}
void uf_unite(int x,int y){
	uf_par[uf_find(y)]=uf_find(x);
}

int IsConnected(int N,int graph[MAX_N][MAX_D]){
	for(int i=0;i<N;i++)uf_par[i]=i;
	for(int i=0;i<N;i++){
		for(int j=0;j<D;j++){
			uf_unite(i,graph[i][j]);
		}
	}

	for(int i=0;i<N;i++)if(uf_find(i)!=uf_find(0))return 0;
	return 1;
}

int IsValid(int N,int graph[MAX_N][MAX_D],int color[MAX_N]){
	for(int i=0;i<N;i++){
		for(int j=0;j<D;j++){
			if(graph[i][j]<0||graph[i][j]>=N){
				puts("out of [0,N)");
				//abort();	
				return 0;
			
			}
			for(int k=0;k<D;k++){
				if(j!=k&&graph[i][j]==graph[i][k]){
					puts("not distinct");
					//abort();
					return 0;
				}
			}
		}
	}

	if(!IsConnected(N,graph)){
		puts("disconnected");
		//abort();	
		return 0;
	}

	int cnt[MAX_C]={0};
	for(int i=0;i<N;i++)cnt[color[i]]++;
	for(int i=0;i<C;i++){
		if(cnt[i]!=M[i]){
			puts("color error");
			//abort();
			return 0;
		}
	}

	return 1;//ok
}

Solver_.c(ぼく)

#define ANSWER_FILE "answer.txt"

#ifndef _GLOBAL
#define _GLOBAL
#include "./Global.c"
#endif

#ifndef _GRAPHEVALUATOR
#define _GRAPHEVALUATOR
#include "./GraphEvaluator.c"
#endif

#ifndef _RANDOMMANAGER
#define _RANDOMMANAGER
#include "./RandomManager.c"
#endif

int answerASPL = INT_MAX;
double destroyedProb = 0.188;

void SaveConfig()
{
	FILE *fp = fopen("config.txt", "w");
	if (fp == NULL)
	{
		puts("Config open failed.");
		return;
	}
	fprintf(fp, "%.16f\n", destroyedProb);
	fclose(fp);
}

void LoadConfig()
{
	FILE *fp = fopen("config.txt", "r");
	if (fp == NULL)
	{
		puts("Config is not found. file is made.");
		SaveConfig();
		return;
	}
	double temp = destroyedProb;
	if (fscanf(fp, "%lf", &destroyedProb) == EOF || destroyedProb < 0 || 1 < destroyedProb)
	{
		puts("Config is invalid.");
		destroyedProb = temp;
		fclose(fp);
		return;
	}
	fclose(fp);
	double abs = temp - destroyedProb;
	if (abs < 0) abs *= -1;
	if (abs > 1e-8) puts("Config is updated");
}

void Save(int graph[MAX_N][MAX_D])
{
	if (IsValid(N, graph, Color) == 0)
	{
		puts("Answer is invalid.");
		return;
	}
	int ret = ASPL(N, graph, Color);
	if (answerASPL > ret)
	{
		answerASPL = ret;
		printf("Answer is updated!!(ASPL = %d)\n", ret);
	}
	else return;
	char str[20] = ANSWER_FILE;
	
	FILE *fp = fopen(str, "w");
	if (fp == NULL)
	{
		puts("File open failed.");
		return;
	}
	int cnt = 0;
	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < D; j++)
		{
			if (i >= graph[i][j]) continue;
			fprintf(fp, "%d %d\n", i, graph[i][j]);
			cnt++;
		}
	}
	fclose(fp);
	if (cnt != N*D/2)
	{
		puts("few edges.");
		return;
	}
}

void Load(int graph[MAX_N][MAX_D])
{
	char str[20] = ANSWER_FILE;
	int size[MAX_N] = {};
	FILE *fp = fopen(str, "r");
	if (fp == NULL)
	{
		puts("File is not found. Graph is generated.");
		generate_random_graph(N, D, graph);
		return;
	}
	for (int i = 0; i < N*D/2; i++)
	{
		int a, b;
		if (fscanf(fp, "%d%d", &a, &b) == EOF)
		{
			puts("Answer is invalid. Graph is generated.");
			generate_random_graph(N, D, graph);
			fclose(fp);
			return;
		}
		graph[a][size[a]++] = b;
		graph[b][size[b]++] = a;
	}
	fclose(fp);
	puts("Answer was loaded.");
}

int TrySwapEdge(int graph[MAX_N][MAX_D], int u, int a, int v, int b)
{
	if (graph[u][a] == v || graph[v][b] == u || graph[u][a] == graph[v][b])
		return 0;
	for (int i = 0; i < D; i++)
	{
		if (graph[u][i] == graph[v][b]) return 0;
	}
	for (int i = 0; i < D; i++)
	{
		if (graph[v][i] == graph[u][a]) return 0;
	}
	int c, d;
	for (int i = 0; i < D; i++)
	{
		if (graph[graph[u][a]][i] == u) { c = i; break; }
	}
	for (int i = 0; i < D; i++)
	{
		if (graph[graph[v][b]][i] == v) { d = i; break; }
	}
	Swap(graph[graph[u][a]] + c, graph[graph[v][b]] + d);
	Swap(graph[u] + a, graph[v] + b);
	return 1;
}

long long EvaluateVertex(int n, int graph[MAX_N][MAX_D], int s)
{
	int dist[MAX_N];
	GetDistance(s, N, graph, Color, dist);
	long long res = 0;
	for (int i = 0; i < N; i++)
	{
		if (dist[i] > N) return LLONG_MAX/2;
		res += dist[i]*dist[i];
	}
	return res;
}

int Search(int graph[MAX_N][MAX_D])
{
	int loop = 10000;
	long long allCount = 0;
	int graphCopy[MAX_N][MAX_D] = {}, asplCopy = INT_MAX/2;
	Save(graph);
	while (1)
	{
		LoadConfig();
		for (int step = 0; step < 10; step++)
		{
			int ret = ASPL(N, graph, Color);
			if (ret - asplCopy > 200)
			{
				for (int i = 0; i < N; i++) for (int j = 0; j < D; j++)
					graph[i][j] = graphCopy[i][j];
				ret = asplCopy;
			}
			else
			{
				for (int i = 0; i < N; i++) for (int j = 0; j < D; j++)
					graphCopy[i][j] = graph[i][j];
				asplCopy = ret;
			}
			printf("score = %d, prob = %f, allCount = %lld\n",
				ret, destroyedProb, allCount);
			for (int cnt = 0; cnt < loop; cnt++, allCount++)
			{
				long long score = 0;
				int r = Rand(0, N*D), s = Rand(0, N*D);
				int u = r/D, a = r%D, v = s/D, b = r%D;
				score +=
					EvaluateVertex(N, graph, u) + EvaluateVertex(N, graph, v) +
					EvaluateVertex(N, graph, graph[u][a]) + EvaluateVertex(N, graph, graph[v][b]);
				
				int isAble = TrySwapEdge(graph, u, a, v, b);
				
				score -=
					EvaluateVertex(N, graph, u) + EvaluateVertex(N, graph, v) +
					EvaluateVertex(N, graph, graph[u][a]) + EvaluateVertex(N, graph, graph[v][b]);
				
				if (score <= 0)
				{
					if (isAble) TrySwapEdge(graph, u, a, v, b);
					else
					{
						if (DoubleRand() > destroyedProb) TrySwapEdge(graph, u, a, v, b);
					}
				}
			}
		}
		Save(graph);
	}
	return ASPL(N, graph, Color);
}

int main()
{
	Input();
	puts("FinishedInput");
	puts("l a t t e");
	int graph[MAX_N][MAX_D];
	Load(graph);
	printf("%d\n", Search(graph));
	return 0;
}

RandomManager.c

追記: この乱択管理には確率的に偏りがあります。特にShuffleArray()については同じ計算量で偏りのない実装があります。(とある偽シャッフルアルゴリズムとその分布 - 簡潔なQ)

//重複しない乱択を行うための配列を管理する関数群
//配列の中身はこれらの関数を使ってランダムに保たれるべき
//特に、要素を使用した後は必ずShuffleElementすること
#ifndef _GLOBAL
#define _GLOBAL
#include "Global.c"
#endif

//配列の要素をランダムに入れ替える
void ShuffleElement(int n, int a[], int pos)
{
	int r = Rand(0, n);
	if (pos == r) return;
	Swap(a + pos, a + r);
}

//配列全体をシャッフル
void ShuffleArray(int n, int a[])
{
	for (int i = 0; i < n; i++)
	{
		ShuffleElement(n, a, i);
	}
}

//配列に新たな要素を追加する
void AddToArray(int oldN, int a[], int value)
{
	a[oldN] = value;
	ShuffleElement(oldN + 1, a, oldN);
}