SWERC 2011 Problem B Coin Collectingとマトロイド交差まとめメモ

問題

n組の封筒のペア(p[i], q[i])がある。
p[i], q[i]の封筒にはそれぞれ、コインが2枚ずつ入っている。
それぞれのコインは整数で表される種類がある。

n組のペアについて、高々片方の封筒を取っていく。（どちらも取らないことができる）
こうして取り終えた後、取った集合をSとする。
Sの任意の部分集合Tに含まれるコインについて、全てのコインが偶数枚あってはならない。
（すべてのTに対して、奇数枚のコインが一種類以上存在すればよい）

このようにして取れるコインの枚数の最大値を求めよ。

制約条件

n≦300
コインの種類の整数≦10000

方針

マトロイド交差。

このまずはこの問題をグラフに落としこむ。
頂点をコインの種類として、
封筒に2枚のコインa, bが入っているとするとき、a, bに辺を張ったグラフを考える。

問題の条件は二つあって、

任意の部分集合に対して全てのコインが偶数枚だとだめ
ペアになっている封筒は片方しか取れない

前者はグラフで考えると、グラフに閉路がないことが必要十分条件。
後者は、ペアの辺のうち片方しか取れない。

で、前者はグラフが無向森であることに相当するから、解の集合はマトロイド。
後者も、封筒の選び方の集合S, Tかつ|S| > |T|があったとして、
するとSには、Tで一個も使っていないペアがあるので、それを加えて有効な解にすることができるので、マトロイド。

この二つのマトロイド交差の最大独立集合が、求める答えである。

マトロイドの定義

集合EおよびF⊆2^Eがあったとき、次の３つが成り立つとき(E, F)はマトロイド。

{}∈F
X⊆Yについて、Y∈FならばX∈F
X, Y∈Fについて|Y| > |X|ならばx∈|Y| - |X|が存在しX∪{x} ∈ F

Eはグラフの辺で、Fは森の全体とかを考えるとマトロイドになっている。

マトロイド交差の最大独立集合問題

マトロイドにおける独立というのは、単にマトロイドの要素であること。
辞書的に言えば、マトロイド(E, F)においてX∈Fであるとき、集合Xは独立であるという。
最大独立問題は、要素が一番多いXを求めることなので、
たとえば森のマトロイドだったら、最大全域森を求めることになる。

二つのマトロイド"交差"の最大独立集合は、二つのマトロイド(E, F1), (E, F2)に対して、X∈F1∩F2であるようなもののうち、要素数が最大のもの。

マトロイド交差の最大独立集合は、２つのマトロイドまでなら多項式時間で求められる。
Edmonds' matroid intersection algorithmというアルゴリズムがある。
http://www.mathematik.hu-berlin.de/~wilhelm/greedy-ausarbeitung.pdf
３つ以上のマトロイド交差はNP-hardらしい。

二部グラフの最大マッチングみたいな感じなのを繰り返すアルゴリズム。
ただしここでの二部グラフの頂点は、E.
つまり森のマトロイドでいうと、辺が頂点になったグラフ。
片側が解の基に入れる要素で、もう片側が解に入れてない要素。

Xを空から始めて、
片方のマトロイドの基を一個ずつ追加していって、追加できるなら追加する。
もしXと追加する基がバッティングしたら、二部グラフ最大マッチングのDFSみたいのをやって、
E-Xの頂点とXの頂点をいったりきたりして増加パスを探して、Xを修正するってのを繰り返す。

交差させるマトロイドをA, Bとしたら、
E-X→Xへの辺は、それを加えるとマトロイドAで破綻する（サーキットができる）部分の辺。
X→E-Xへの辺は、それを加えるとマトロイドBで破綻する（サーキットができる）部分の辺。

サーキットは、独立（マトロイドの要素）でない集合のうち極小なもので、
森のマトロイドだとまさに閉路。自分を入れると閉路ができるような閉路の辺全てに対して辺を張る。

方針続き

Edmondsのマトロイド交差アルゴリズムなんだけど、
使用済み頂点から未使用頂点に行くときは、後者（ペアの片方しか使えない）の条件でサーキットになる頂点へ、
未使用頂点から使用済み頂点にいくときは、前者（閉路ができない）の条件でサーキットになる頂点へ辺を張るようにする。

証明とか

誰か僕にもわかるように教えて……

ソースコード

#define F first
#define S second

int n, N;
int a[1212], b[1212];

vector<vector<pi> > e;

int v[1212];
int st[1212], sz;

inline int rec(int c, int p, int depth = 1){
	v[c] = depth;
	rep(i, e[c].size()) if(e[c][i].S != p){
		
		st[sz++] = e[c][i].S;
		
		if(v[e[c][i].F]) return depth + 1 - v[e[c][i].F];
		
		int res = rec(e[c][i].F, e[c][i].S, depth + 1);
		if(res) return res;
	
		sz--;
	}
	return 0;
}

inline vi getcycle(int id){
	
	vi res;
	
	e[a[id]].pb(mp(b[id], id));
	e[b[id]].pb(mp(a[id], id));
	
	rep(i, N) v[i] = 0;
	rep(i, N) if(!v[i]){
		sz = 0;
		int len = rec(i, -1);
		if(len == 0) continue;
		
		rep(j, len) res.pb(st[sz - 1 - j]);
		break;
	}
	
	e[a[id]].pop_back();
	e[b[id]].pop_back();
	
	return res;
}

bool can(int id){
	
	vi cy = getcycle(id);
	
	if(cy.empty()){
		e[a[id]].pb(mp(b[id], id));
		e[b[id]].pb(mp(a[id], id));
		return 1;
	}
	
	
	int prev[1212];
	queue<pi> q;
	
	memset(prev, -1, sizeof(prev));
	prev[id] = id;
	q.push(mp(id, 0));
	
	while(!q.empty()){
		
		int cur = q.front().F, isX = q.front().S;
		q.pop();
		
		if(isX){
			if(prev[cur ^ 1] >= 0) continue;
			prev[cur ^ 1] = cur;
			q.push(mp(cur ^ 1, 0));
		}
		else{
			vi cy = getcycle(cur);
			
			if(cy.empty() || cy.size() == 1 && cy[0] == cur){
				
				while(1){
					bool found = 0;
					
					rep(i, N) rep(j, e[i].size()) if(e[i][j].S == cur){
						found = 1;
						break;
					}
					if(!found){
						e[a[cur]].pb(mp(b[cur], cur));
						e[b[cur]].pb(mp(a[cur], cur));
					}
					else{
						rep(it, N){
							int i = ii ? a[cur] : b[cur];
							
							for(int j = (int)e[i].size() - 1; j >= 0; j--){
								
								if(e[i][j].S != cur) continue;
								swap(e[i][j], e[i].back());
								e[i].resize(e[i].size() - 1);
							}
						}
					}
					if(cur == id) break;
					cur = prev[cur];
				}
				return 1;
			}
			rep(i, cy.size()) if(cy[i] != cur){
				
				if(prev[cy[i]] >= 0) continue;
				prev[cy[i]] = cur;
				q.push(mp(cy[i], 1));
			}
		}
	}
	return 0;
}

bool run(){
	
	scanf("%d", &n);
	if(n == 0) return 0;
	
	e.clear();
	e.resize(4*n);
	
	vi vs;
	
	rep(i, n){
		scanf("%d%d%d%d", a + 2*i, b + 2*i, a + 2*i+1, b + 2*i+1);
		vs.pb(a[i * 2]); vs.pb(a[i * 2 + 1]);
		vs.pb(b[i * 2]); vs.pb(b[i * 2 + 1]);
	}
	sort(all(vs));
	vs.erase(unique(all(vs)), vs.end());
	N = vs.size();
	
	rep(i, 2 * n){
		a[i] = lower_bound(all(vs), a[i]) - vs.begin();
		b[i] = lower_bound(all(vs), b[i]) - vs.begin();
	}
	
	int ans = 0;
	
	rep(i, n){
		if(can(i * 2) || can(i * 2 + 1)) ans += 2;
	}
	cout << ans << endl;
	
	return 1;
}
int main(){
	while(run());
	return 0;
}