supercon参加記 その2

優勝しました!


最初は頭が全く回ってなかったせいもあり、サンプルに騙されてたけど、頭回り出してからは順調に行った。
まずはCPUで正しく動くフローを書こうと思い、ポテンシャルを使ったダイクストラ解法で解いた。
最初に出来たプログラムで既に0.04sくらいで、一番遅いダイクストラ部分が自明に一回減るのでその辺の高速化とかをしたら0.02sくらいになった。
3日目の午前中くらいにはもうデバッグも終わっていたので、GPUによる高速化を考えていた。
グラフに特徴があり、ベルマンフォードでも更新回数が少なくて済むことが分かったので、これを書いてGPU化することにした。
ベルマンフォードは前回更新があった頂点だけ調べるような枝狩りみたいなのもしてCPUだと0.15sくらいだったので、うまく並列化したらいい感じになりそうとか思ってたが、最後までバグが取れずよく分からない。
これと同時くらいに仲間がテストしてくれてて、サンプルと答えがずれるケースをいくつか検出してて焦った。
結局はテストケースジェネレータのミスだったので安心した。
以下コード。
こういうグラフの持ち方すると、逆辺がxor1するだけでわかるので素晴らしい。
詳しくはここにあるスライドにあります

/*******************************************************************
  SuperCon12: Template for input output procedures
  file: sc12_template2E.cu
  by S.Kishimoto and O.Watanabe Aug. 17, 2012
  ver.2 by T.Endo and A.Nukada Aug. 22, 2012

Remark: @@@ for test/tentative statements
*********************************************************************/

/********* Template Part (1) NO CHAGNGE FROM HERE **********/
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<sys/time.h>
// @@@ two include statements for CUDA
#include<cuda.h>
#include<cuda_runtime.h>

// @@@
#define MAXlid  5
#define MAXst   40
#define MAXntr  180
#define MAXtd   7
#define MAXtid  6000
// @@@
//#define MAXlid  3
//#define MAXst   20
//#define MAXntr  90
//#define MAXtd   4
//#define MAXtid  4000

#define MAXtm4test 1200
//#define MAXtm4test 120

#define MAXtm   1200
#define MAXtvtm 7
#define MAXstPline 20

//------------
int Nst;
int Ntid;
int maxDist;
short int trainTBnst[MAXtid];
short int trainTB[MAXtid][MAXstPline][3];
int outroute[MAXtm][3][2];

//------------
void input() {
        int i,j,k;
	scanf("%d", &Nst);

	int nlid;
	scanf("%d", &nlid);
	int yobi;
	for(i=0; i<nlid; i++) {
		int m;
		scanf("%d", &m);
		for(j=0; j<m; j++) scanf("%d", &yobi);
	}
	for(i=0; i<Nst; i++) for(j=0; j<Nst; j++) scanf("%d", &yobi);

	scanf("%d", &Ntid);
	int tid;
	for(tid=0; tid<Ntid; tid++) {
		int m;
		scanf("%d", &m);
		trainTBnst[tid] = m;
		int j;
		for(j=0; j<m; j++) {
			int st , tm, dist;
			scanf("%d %d %d", &st, &tm, &dist);
			trainTB[tid][j][0] = st;
			trainTB[tid][j][1] = tm;
			trainTB[tid][j][2] = dist;
		}
	}

	for(i=0; i<MAXtm; i++) for(j=0; j<3; j++) for(k=0; k<2; k++) outroute[i][j][k] = -1;

	// pay "CUDA tax" now
	{
		void *p;
		cudaError_t err;
		err = cudaMalloc(&p, 4096);
		if (err != cudaSuccess) {
			fprintf(stderr, "cudaMalloc in input() failed!\n");
			exit(1);
		}
		err = cudaFree(p);
		if (err != cudaSuccess) {
			fprintf(stderr, "cudaFree in input() failed!\n");
			exit(1);
		}
	}
}

void output() {
	int i,j;
	printf("##########\n");
	printf("maxDist: %d\n",maxDist);
	for(i=0; i<MAXtm; i++) {
		for(j=0; j<3; j++)  printf("%2d %2d  ", outroute[i][j][0], outroute[i][j][1]);
		printf("\n");
	}
}
/********* Template Part (1) NO CHANGE UP TO HERE **********/

/********* Your Program Part (1) FROM HERE **********/

#define MAX_V 100000
#define MAX_E 500000
#define MAX_T 150000
#define INF 10000000

struct edge{ int to, cost, next; short cap;};

int V;
int head[MAX_V], graph_i = 0;
struct edge graph[MAX_E];
int trainNum[MAX_T];

int h[MAX_V];
int dist[MAX_V];
int preve[MAX_V];

void add_edge(int from, int to, int cost){
    graph[graph_i] = (struct edge){to,-cost,head[from],1}; head[from] = graph_i++; 
    graph[graph_i] = (struct edge){from,cost,head[to],0}; head[to] = graph_i++; 
}

int calc_V(int i, int j){
    return (int)trainTB[i][j][0]+(Nst*(int)trainTB[i][j][1]<<1);
}



#define MAX_NODE 1000000

struct node{ int val,v;};

struct node heap[MAX_NODE];
int heap_size = 0;

void push(struct node n) {
    int i = heap_size++;
    while(i > 0) {
        int p = (i - 1) >> 1;
        if(heap[p].val <= n.val) break;
        heap[i] = heap[p];
        i = p;
    }
    heap[i] = n;
}

struct node pop() {
    struct node ret = heap[0];
    struct node n = heap[--heap_size];
    int i = 0;
    while((i<<1)+1 < heap_size) {
        int a = (i<<1)+1, b = a+1;
        if(b < heap_size && heap[b].val < heap[a].val) a = b;
        if(heap[a].val >= n.val) break;
        heap[i] = heap[a];
        i = a;
    }
    heap[i] = n;
    return ret;
}


/********* Your Program Part (1) UP TO HERE **********/

/********* Template Part (2) NO CHANGE FROM HERE **********/
int main() {
	struct timeval tstart, tlast;
	input();
	gettimeofday(&tstart, NULL);
/********* Template Part (2) NO CHANGE UP TO HERE **********/

/********* Your Program Part (2) FROM HERE **********/
    
    int i, j;
    int SV, TV;
    int nv, nto, ndist;
    
    V = (Nst*MAXtm<<1)+2;
    SV = V-2; TV = V-1;
    
    //init graph
    for(i = 0; i < V; i++) head[i] = -1;
    
    for(i = 0; i < Ntid; i++){
        for(j = 1; j < trainTBnst[i]; j++){
            trainNum[graph_i] = i;
            add_edge(calc_V(i,j-1)+Nst,calc_V(i,j),trainTB[i][j][2]);
        }
    }
    
    for(i = 1; i < (MAXtm<<1); i++) for(j = 0; j < Nst; j++) add_edge((i-1)*Nst+j,i*Nst+j,0);
    for(i = 0; i < Nst; i++){
        add_edge(SV,i,0);
        add_edge(Nst*((MAXtm<<1)-1)+i,TV,0);
    }   
    //
    
    //init potential
    for(i = 0; i < V; i++) h[i] = INF;
    h[SV] = 0;
    for(i = 0; i < V; i++){
        if(i){
            if(i == TV) nv = TV; else nv = i-1;
        } else nv = SV; // SV 0 1 2 3 ... V-3 TV
        for(j = head[nv]; j != -1; j = graph[j].next){
            if((j&1) == 0 && h[graph[j].to] > h[nv]+graph[j].cost){
                h[graph[j].to] = h[nv]+graph[j].cost;
                preve[graph[j].to] = j;
            }
        }
    }
    //
    
    //minimum cost flow
    int flow = 3, ans = 0;
    struct node xnode;

    while(flow){
        if(flow != 3){
            for(i = 0; i < V; i++) dist[i] = INF;
            dist[SV] = 0;
            heap_size = 0;
            push((struct node){0,SV});

            while(heap_size){
                xnode = pop();
                nv = xnode.v;
                if(dist[nv] < xnode.val) continue;

                for(i = head[nv]; i != -1; i = graph[i].next){
                    if(graph[i].cap == 0) continue;
                    nto = graph[i].to;
                    ndist = dist[nv] + graph[i].cost + h[nv]-h[nto];
                    if(dist[nto] > ndist){
                        dist[nto] = ndist;
                        preve[nto] = i;
                        push((struct node){ndist,nto});
                    }
                }
            }

            for(i = 0; i < V; i++) h[i] += dist[i];
        }
        
        ans += h[TV];
        for(nv = TV; nv != SV; nv = graph[i^1].to){
            i = preve[nv];
            graph[i].cap = 0;
            graph[i^1].cap = 1;
        }
        
        flow--;
    }
    //
    
    maxDist = -ans;
    
    //path
        int ntm;
        for(flow = 0; flow < 3; flow++){
            nv = SV; ntm = 0;
            while(nv != TV){
                for(i = head[nv]; i != -1; i = graph[i].next){
                    if((i&1) == 1 || graph[i].cap) continue;
                    graph[i].cap = 1;
                    nto = graph[i].to;
                    
                    if((nv/Nst&1) == 1){    
                        outroute[ntm][flow][0] = nv%Nst;
                        j = ntm;
                        ntm = nto/Nst>>1;
                        if(graph[i].cost){
                            for(; j < ntm; j++){
                                outroute[j][flow][1] = trainNum[i];
                            }
                        }
                    }
                    
                    nv = nto;
                    break;
                }
            }
        }
    //*/

/********* Your Program Part (2) UP TO HERE **********/

/********* Template Part (3) NO CHANGE FROM HERE TO THE END **********/
	gettimeofday(&tlast, NULL);
	printf("...end of execution... time %f\n", 
		   (double)(tlast.tv_sec - tstart.tv_sec)+
		   (double)(tlast.tv_usec - tstart.tv_usec)/ 1000000);
	output();
	return 0;
}

フィボナッチヒープをダイクストラと組み合わせると1桁ms位になることを後で聞いたけど、思いつかなかったなぁ・・・
準急がさらにヤバい解法をつぶやいてて怖かった。
こんな普通めな解法で優勝できるとは思ってなかった・・・
PANAIは全員が別々の方針を実装していたりすれば負けていたかも。


まとめ
高校最後の大会で優勝できて本当に嬉しい。
「優勝は、イマセン!」とかなると面白いなぁーとか思いながら付けたチーム名なので、当初の目的を達成できた。
一人では絶対に優勝できなかった。チームメイトにまじで感謝してる。