SuperCon2010

先日、東京工業大学大阪大学の方でSuperConという大会がありました。
スーパーコンピュータを使わせていただくことの出来る太っ腹な大会です。

自分は、LovePlusというチーム名で大阪会場の方に参加させていただきました。

予選

2*n及び3*mの広場を1*2の大きさのタイルで敷き詰める通り数を求める問題。
普通に漸化式を導き行列累乗をバイナリ法を使ってlog(n)で求めているだけですが、興味のある方は下のスライドをどうぞ。

本選

n*m (n, m≦50)の広場を1*2の大きさのタイルで敷き詰める通り数を求める問題。障害物もあります。
これが100000問出題されるので、東工大のスーパーコンピュータ・TSUBAMEの9ノード*8コアを45分間使ってできるだけ多く解いたチームの勝ち。

  • 島ごとに分けて計算し後で乗算
  • ボトルネックになっている1lineを分割し再帰的に計算
  • ここに置かなければならないと確定する所は先に置く
  • 2個同時置き

以上を実装し、練習用の入力ファイル(10000問)では最終日では多分5000問弱位解けるようになっていました。
また、並列化に関しては

  • statistics.txt(問題の統計情報が書かれているファイル)を読み込み評価関数を通し簡単な順にソート
  • 100問ごとに同じファイルに分けて読み込む

の2点に重きを置きました。
これでまぁまぁ解けるのではないかなと思っていたのですが、本番で評価関数を通した際に何故か難しい問題が上位に混入してしまい
スーパーコンピュータさんがその問題を解こうと必死にうんうん唸り始めてしまったため残念なことになってしまいました。
健気で一途なのはいいのですが、タイムアウト機能を入れておけばよかったなと激しく思います。
あと最後の最後まで動的計画法を実装が面倒そうだからという思考実験だけで完結させてしまっていたことも心残りです。


そんな感じで余りいい結果は残せませんでしたが、スーパーコンピュータに触れるというだけでも非常に良い経験になったと思います。
SuperConは今年で最後でしたが、まだ高専プロコンやパソコン甲子園がありますし、来年以降はICPCなどでも活躍していけたらいいなと思っています。

最後にソースコードを貼っておきます。
コメントは基本的に嘘つきです。

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <math.h>
#include "sc10.h"

#define TEAMNAME "loveplus"

int RANK;		// 自分のRank
int PROCESS; 	// 全体のプロセス数
int TOTAL;		// 問題の数

void initialize(int* argc, char* argv[]);
void finalize();
void parent_arg(char* argv[]);
void parent_all(int argc, char* argv[]);
void child(char* argv[]);
int solve(int m, int n, int k, char table[MAXSIZE][MAXSIZE]);
int compute1(int m, int n, int k, char table[MAXSIZE][MAXSIZE]);
void unionfind(int m, int n, char table[MAXSIZE][MAXSIZE], int set[], int* setn);
void fill(int m, int n, char table[MAXSIZE][MAXSIZE], int y, int x, char c);
void animate(int m, int n, char table[MAXSIZE][MAXSIZE]);
void getseparate(int m, int n, char table[MAXSIZE][MAXSIZE], int* horizon, int* pos);
int separate(int m, int n, int k, char table[MAXSIZE][MAXSIZE], int horizon, int pos, int now);

typedef struct
{
	int s_filenum;
	int probnum;
	double score;
}sub_Level;

double serch_score2(int n, int m, int obs, double s1, double s2, double s3);
void fsort(char* statistics, int fnames[], int* qnum);
int compare2(const void* a, const void* b);


int main(int argc, char* argv[])
{
	int i;
	
	if(argc == 1){
		puts("Usage: loveplus (filename | -a dirname)");
		return 1;
	}
	
	initialize(&argc, argv);
	
	if(argc == 2)
		TOTAL = sc10_readfile(argv[1]);
	
	if(RANK == 0){
		if(argc == 2)
			parent_arg(argv);
		else
			parent_all(argc, argv);
	}else{
		child(argv);
	}
	
	finalize();
	return EXIT_SUCCESS;
}

/* ------------------------------------
	rank0 のプロセスの処理(引数1個, デバッグ用)
---------------------------------------*/
void parent_arg(char* argv[])
{
	int i, sent=0, solved=0;
	int problems[100000], pnum;
	
	for(i=1; i<PROCESS; ++i){
		MPI_Send(&sent, 1, MPI_INT, i, 99, MPI_COMM_WORLD);
		++sent;
	}
	
	while(solved<TOTAL){
		MPI_Status status;
		int result[2];
		
		MPI_Recv(result, 2, MPI_INT, MPI_ANY_SOURCE, 99, MPI_COMM_WORLD, &status);
		sc10_output(result[0], result[1]);
		++solved;
		
		if(sent < TOTAL){
			MPI_Send(&sent, 1, MPI_INT, status.MPI_SOURCE, 99, MPI_COMM_WORLD);
			++sent;
		}
	}
	
	MPI_Abort(MPI_COMM_WORLD, -1);
}

/* ------------------------------------
	rank0 のプロセスの処理(引数2個, 本番用)
---------------------------------------*/
void parent_all(int argc, char* argv[])
{
	int i, j, c=0;
	int problems[100000], qnum;
	
	char statistics[256];
	sprintf(statistics, "%s/statistics.txt", argv[2]);
	fsort(statistics, problems, &qnum);
	
	i=0;
	if(argc == 4)
		i = atoi(argv[3]);
	
	for(; i<qnum/100; i+=c){
		int sent = 0, solved = 0;
		char fname[16];
		
		for(c=1; (i+c)*100<qnum; ++c){
			if(problems[i*100]/1000 != problems[(i+c)*100]/1000)
				break;
		}
		
		for(j=1; j<PROCESS; ++j){
			int hoge = -(problems[i*100]/1000%100)-1;
			MPI_Send(&hoge, 1, MPI_INT, j, 99, MPI_COMM_WORLD);
		}
		
		sprintf(fname, "%s/prob%02d.in", argv[2], problems[i*100]/1000%100);
		fprintf(stderr, "%d: %s\n", i*100, fname);
		sc10_readfile(fname);
		
		for(j=1; j<PROCESS; ++j){
			int hoge = problems[i*100 + sent] % 1000;
			MPI_Send(&hoge, 1, MPI_INT, j, 99, MPI_COMM_WORLD);
			fprintf(stderr, "%d to %d\n", j, problems[i*100+sent]);
			++sent;
		}
		
		while(solved < 100*c){
			MPI_Status status;
			int result[2];
			
			MPI_Recv(result, 2, MPI_INT, MPI_ANY_SOURCE, 99, MPI_COMM_WORLD, &status);
			sc10_output(result[0], result[1]);
			++solved;
			
			if(sent < 100*c){
				int hoge = problems[i*100 + sent] % 1000;
				MPI_Send(&hoge, 1, MPI_INT, status.MPI_SOURCE, 99, MPI_COMM_WORLD);
				fprintf(stderr, "%d to %d\n", status.MPI_SOURCE, problems[i*100+sent]);
				++sent;
			}
		}
		
	}
}

/* ------------------------------------
	rank0 以外のプロセスの処理
---------------------------------------*/
void child(char* argv[])
{
	for(;;){
		MPI_Status status;
		int index, result[2];
		
		int probname, m, n, k;
		char table[MAXSIZE][MAXSIZE];
		
		MPI_Recv(&index, 1, MPI_INT, 0, 99, MPI_COMM_WORLD, &status);
		
		if(index < 0){
			char fname[16];
			sprintf(fname, "%s/prob%02d.in", argv[2], -index-1);
			
			sc10_readfile(fname);
			
		}else{
			sc10_getproblem(index, &probname, &m, &n, &k, table);
			result[0] = probname;
			result[1] = solve(m, n, k, table);
			
			MPI_Send(result, 2, MPI_INT, 0, 99, MPI_COMM_WORLD);
		}
	}
}

int solve(int m, int n, int k, char table[MAXSIZE][MAXSIZE])
{
	int i, j, o;
	int setn;
	int set[MAXSIZE * MAXSIZE];
	int result = 1;
	char buf[MAXSIZE][MAXSIZE], buf2[MAXSIZE][MAXSIZE];
	
//	fprintf(stderr, "hogehoge!\n");
//	animate(m, n, table);
//	sleep(1);
	
	unionfind(m, n, table, set, &setn);
	
	for(i=0; i<setn; ++i){
		int left, right, top, bottom;
		int pm, pn;
		int tmpresult;
		
		for(j=0; j<m; ++j){
			for(o=0; o<n; ++o){
				buf[j][o] = table[j][o] == '.' ? '#' : table[j][o];
			}
		}
		
		fill(m, n, buf, set[i]/n, set[i]%n, '.');
		
		
		left = top = MAXSIZE;
		right = bottom = 0;
		
		for(j=0; j<m; ++j){
			for(o=0; o<n; ++o){
				if(buf[j][o] == '.'){
					top = top < j ? top : j;
					bottom = j < bottom ? bottom : j;
					left = left < o ? left : o;
					right = o < right ? right : o;
				}
			}
		}
		
		pm = bottom - top + 1;
		pn = right - left + 1;
		
		for(j=0; j<MAXSIZE; ++j){
			for(o=0; o<MAXSIZE; ++o){
				buf2[j][o] = '*';
			}
		}
		
		for(j=top; j<=bottom; ++j){
			for(o=left; o<=right; ++o){
				buf2[j-top][o-left] = buf[j][o] == '#' ? '*' : buf[j][o];
			}
		}
		
		if(pm * pn <= 75){
			tmpresult = compute1(pm, pn, k, buf2);
			
		}else{
			int horizon, pos;
			getseparate(pm, pn, buf2, &horizon, &pos);
			//fprintf(stderr, "horizon=%d, pos=%d\n", horizon, pos);
			tmpresult = separate(pm, pn, k, buf2, horizon, pos, 0);
		}
		
		result = (result * tmpresult) % k;
	}
	
//	fprintf(stderr, "%d returned.\n", result);
//	sleep(1);
	
	return result;
}

/* ------------------------------------
	分割
---------------------------------------*/
void getseparate(int m, int n, char table[MAXSIZE][MAXSIZE], int* horizon, int* pos)
{
/**
	int hor[MAXSIZE], ver[MAXSIZE];
	int hsum[MAXSIZE] = {0}, vsum[MAXSIZE] = {0};
	int i, j;
	int best, bestdiff;
	
	for(i=0; i<m; ++i){
		for(j=0; j<n; ++j){
			hor[i] += table[i][j] == '.';
			ver[j] += table[i][j] == '.';
		}
	}
	
	for(i=1; i<=m; ++i){
		hsum[i] = hsum[i-1] + hor[i];
	}
	for(i=1; i<=n; ++i){
		vsum[i] = vsum[i-1] + ver[i];
	}
	
	best = MAXSIZE;
	bestdiff = MAXSIZE;
	
	for(i=0; i<m; ++i){
		if(hor[i] == 0)
			continue;
		
		if(best > hor[i]){
			*pos = i;
			best = hor[i];
			*horizon = 1;
			bestdiff = abs((hsum[i]) - (hsum[m]-hsum[i]));
		}else if(best == hor[i]){
			int diff = abs((hsum[i]) - (hsum[m]-hsum[i]));
			if(diff < bestdiff){
				*pos = i;
				best = hor[i];
				*horizon = 1;
				bestdiff = diff;
			}
		}
	}
	for(j=0; j<n; ++j){
		if(ver[j] == 0)
			continue;
		
		if(best > ver[j]){
			*pos = j;
			best = ver[j];
			*horizon = 0;
		}else if(best == ver[j]){
			int diff = abs((vsum[j]) - (vsum[n]-vsum[j]));
			if(diff < bestdiff){
				*pos = i;
				best = ver[j];
				*horizon = 0;
				bestdiff = diff;
			}
		}
	}
	
	
/*/
	*horizon = m > n;
	*pos = m > n ? m/2 : n/2;
	*pos-=2;
/**/
}

int separate(int m, int n, int k, char table[MAXSIZE][MAXSIZE], int horizon, int pos, int now)
{
	int result = 0;
	int i, j;
	
	int dy = horizon ? 0 : 1;
	int dx = horizon ? 1 : 0;
	
	int y = horizon ? pos : now;
	int x = horizon ? now : pos;
	
	int put = 0;
//	fprintf(stderr, "now = %d started.\n", now);
	
	for(; y<m && x<n; y+=dy, x+=dx){
//		fprintf(stderr, "(%d, %d)\n", x, y);
		
		if(table[y][x] != '.')
			continue;
		
		table[y][x] = '%';
		if(table[y+dy][x+dx] == '.'){
//			fprintf(stderr, "a in(now = %d)\n", now);
			table[y+dy][x+dx] = '*';
			result = (result + separate(m, n, k, table, horizon, pos, horizon ? x+2 : y+2)) % k;
			table[y+dy][x+dx] = '.';
//			fprintf(stderr, "a out(now = %d)\n", now);
			++put;
		}
		if(table[y+dx][x+dy] == '.'){
//			fprintf(stderr, "b in(now = %d)\n", now);
			table[y+dx][x+dy] = '%';
			result = (result + separate(m, n, k, table, horizon, pos, horizon ? x+1 : y+1)) % k;
			table[y+dx][x+dy] = '.';
//			fprintf(stderr, "b out(now = %d)\n", now);
			++put;
		}
		if(table[y-dx][x-dy] == '.'){
//			fprintf(stderr, "c in(now = %d)\n", now);
			table[y-dx][x-dy] = '%';
			result = (result + separate(m, n, k, table, horizon, pos, horizon ? x+1 : y+1)) % k;
			table[y-dx][x-dy] = '.';
//			fprintf(stderr, "c out(now = %d)\n", now);
			++put;
		}
		
		table[y][x] = '.';

		if(put != 0)
			return result;
	}
	result = (result + solve(m, n, k, table)) % k;
//	fprintf(stderr, "now = %d finished.\n", now);
	
	return result;
}

/* ------------------------------------
	unionfindかと思いきや
---------------------------------------*/
void unionfind(int m, int n, char table[MAXSIZE][MAXSIZE], int set[], int* setn)
{
	int i, j;
	char buf[MAXSIZE][MAXSIZE];
	
	*setn = 0;
	
	for(i=0; i<m; ++i){
		for(j=0; j<n; ++j){
			buf[i][j] = table[i][j] == '.' ? '#' : table[i][j];
		}
	}
	
	for(i=0; i<m; ++i){
		for(j=0; j<n; ++j){
			if(buf[i][j] == '#'){
				fill(m, n, buf, i, j, '@');
				set[*setn] = i*n + j;
				++*setn;
			}
		}
	}
}

void fill(int m, int n, char table[MAXSIZE][MAXSIZE], int y, int x, char c)
{
	int dy[] = {-1, 0, 0, 1};
	int dx[] = {0, 1, -1, 0};
	int i, j;
	
	table[y][x] = c;
	
	for(i=0; i<4; ++i){
		int py = dy[i] + y;
		int px = dx[i] + x;
		if(py<0 || m<=py || px<0 || n<=px)
			continue;
		
		if(table[py][px] != '#')
			continue;
		
		fill(m, n, table, py, px, c);
	}
}

/* ------------------------------------
	探索木による計算
---------------------------------------*/
void animate(int m, int n, char table[MAXSIZE][MAXSIZE])
{
	int i, j;
	for(i=0; i<m; ++i){
		for(j=0; j<n; ++j){
			putchar(table[i][j]);
		}
		puts("");
	}
	puts("");
}

int compute1(int m, int n, int k, char table[MAXSIZE][MAXSIZE])
{
	int i, j;
	int tmpresult = 0;
	int by=-1, bx=-1;
	
	for(i=0; i<m; i++){
		for(j=0; j<n; j++){
			int a, b, c, d;
			
			if(table[i][j] != '.')
				continue;
			
			a = table[i][j+1] == '.';
			b = table[i+1][j] == '.';
			c = j ? table[i][j-1] == '.' : 0;
			d = i ? table[i-1][j] == '.' : 0;
			
			if(a+b+c+d>=2 && by == -1){
				by = i;
				bx = j;
				
			}else if(a+b+c+d==1){
				table[i][j] = '*';
				if(a){
					table[i][j+1] = '*';
					tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
					table[i][j+1] = '.';
					
				}else if(b){
					table[i+1][j] = '*';
					tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
					table[i+1][j] = '.';
					
				}else if(c){
					table[i][j-1] = '*';
					tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
					table[i][j-1] = '.';
					
				}else if(d){
					table[i-1][j] = '*';
					tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
					table[i-1][j] = '.';
				}
				table[i][j] = '.';
				return tmpresult;
				
			}
			else if(a+b+c+d==0){
				return 0;
			}
		}
	}
	
	if(by == -1)
		return 1;
	
	/* It's mendokusai. */
	i = by;
	j = bx;
	
	table[i][j] = '*';
	
	/* horizonal */
	table[i][j+1] = '*';
	
	if(table[i+1][j] == '.'){
		table[i+1][j] = '*';
		
		if(table[i+1][j+1] == '.'){
			table[i+1][j+1] = '*';
			tmpresult = (tmpresult + compute1(m, n, k, table) * 2) % k;
			table[i+1][j+1] = '.';
		}
		if(table[i+2][j] == '.'){
			table[i+2][j] = '*';
			tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
			table[i+2][j] = '.';
		}
		if(j>0 && table[i+1][j-1] == '.'){
			table[i+1][j-1] = '*';
			tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
			table[i+1][j-1] = '.';
		}
		table[i+1][j] = '.';
		
	}else{
		tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
	}
	table[i][j+1] = '.';
	
	/* vertical */
	table[i+1][j] = '*';
	
	if(table[i+1][j+1] == '.'){
		table[i+1][j+1] = '*';
		
		if(table[i+2][j+1] == '.'){
			table[i+2][j+1] = '*';
			tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
			table[i+2][j+1] = '.';
		}
		if(table[i+1][j+2] == '.'){
			table[i+1][j+2] = '*';
			tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
			table[i+1][j+2] = '.';
		}
		table[i+1][j+1] = '.';
		
	}else{
		tmpresult = (tmpresult + compute1(m, n, k, table)) % k;
	}
	table[i+1][j] = '.';
	table[i][j] = '.';
	
	return tmpresult;
}

/* ------------------------------------
	初期化処理
---------------------------------------*/
void initialize(int* argc, char* argv[])
{
	MPI_Init(argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &RANK);
	MPI_Comm_size(MPI_COMM_WORLD, &PROCESS);
	sc10_init(TEAMNAME);
}

/* ------------------------------------
	終了時処理
---------------------------------------*/
void finalize()
{
	MPI_Finalize();
}

/* ------------------------------------
	サッ☆
---------------------------------------*/
void fsort(char* statistics, int fnames[], int* qnum)
{
	FILE *fp;
	int i, j,loop,l_count,cp_i[100]={0};
	int n, m, k;
	int obs;
	double s1, s2, s3;
	int pnum;
	char fname[16];
	int cp_filename[100][100];
	
	sub_Level* s_level = (sub_Level*)malloc(sizeof(sub_Level) * 100001);
	
	*qnum=0;
	
	fp = fopen(statistics, "r");
	
	while(fscanf(fp, "%d%s%d%d%d%d%lf%lf%lf", &pnum, fname, &n, &m, &k, &obs, &s1, &s2, &s3) != EOF){
		int f;
		/*問題番号*/
		s_level[*qnum].probnum = pnum;
		s_level[*qnum].s_filenum = pnum/1000%100;
		s_level[*qnum].score = serch_score2(n, m, obs,s1,s2,s3);
		++*qnum;
	}
	
	qsort(s_level, *qnum, sizeof(s_level[0]), compare2);
	
	loop=0;
	
	for(i=0;i<*qnum;i++)
	{
		cp_filename[ s_level[i].s_filenum ][ cp_i[s_level[i].s_filenum ]]=s_level[i].probnum;
		cp_i[ s_level[i].s_filenum ]++;
		if(cp_i[ s_level[i].s_filenum ]==100)
		{
			for(l_count=0;l_count<100;++loop,++l_count)
			{
				fnames[loop]=cp_filename[ s_level[i].s_filenum ][ l_count ];
			}
			cp_i[s_level[i].s_filenum]=0;
		}
	}
	
	free(s_level);
	fclose(fp);
}

double serch_score2(int n, int m, int obs, double s1, double s2, double s3)
{
	double temp = 1.0/(n*m-obs);//sqrt(s1*1.5+s2+s3*0.5)/10000.0;
	double _obs = ((double)obs/((double)n*(double)m));
//	double a = temp;
	int flg = 0;
	
	if(s1 < 60){
		temp *= 0.6;
	}
/*
	if(s1 > 40 && s3 > 200){
		temp *= 0.4;
		_obs = sqrt(_obs);
		flg = 1;
	}
*/
	else if(s1 < 10 && s3 < 70){
		temp *= 0.75;
		_obs = sqrt(_obs);
		flg = 1;
	}
	else if(s1+s2+s3 < 220){
		temp *=0.5;
		_obs = sqrt(_obs);
		flg = 1;
		if(s1+s2+s3 < 40){
			temp *= 0.75;
		}
	}
	if(s1 + s2 + s3 > 700){
		temp *= 3;
	}

	temp *= (_obs/5);
//	printf("%lf * %lf * srore >> %lf\n",a,flg == 0 ? 1.0 : 0.6,temp);
	return temp;
}

int compare2(const void* a, const void* b)
{
	double aa = ((const sub_Level*)a)->score;
	double bb = ((const sub_Level*)b)->score;
	
	return aa > bb ? -1 : (aa < bb ? 1 : 0);
}