USACO Pollutant Control 解题报告

这道题是mincut问题。maxflow-mincut原理是说，一个图的maxflow（同一时刻从源到sink的流的capacity总和）和mincut（使源和目的分开需要切断的边的capacity之和的最小值）是一样大的。

在跑一遍maxflow,得到maxflow值后，如何得到mincut有两种方法。最快的方法是在residual graph上面从源跑一次DFS，看源能到那些节点。然后，对所有的边，如果它的一个节点可以从源到，另一个节点到不了，则说明这条边在最小割里面。方法二是依次隔断每条边，再求一次最大流。如果和原图最大流之差为这条边的capacity,则这条边在最大流里面。显然，方法二需要跑很多次最大流，没有方法一效率高。

这道题的test cases很tricky,每个测试点都需要改一遍程序。

首先，两个点之间可能会有多条边。这里的处理方法之一是新构建一个点，好处是这样每条边都可以单独处理，和普通的图一模一样了，缺点自然是节点可能会多很多，过不了有些测试点。处理方法之二是把两点之间所有边的capacity加到一起。这也是最有效的方法。

其次，maxflow可能会超过int范围。用过求capacity公约数的方法，但最简单的还是用unsigned long long或者double。这里用的是double。

再次，题目对输出的最小割有限定。首先得是边数最小的最小割，其次，得是index最小的。要求边数最小，可以把capacity * (M + 1) + 1。乘以(M + 1)是为了保证结果仍然是最小割，而不仅仅是因为边数少。+1是为了得到边数最小的最小割。详细解释见http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=maxFlow2。为了达到index最小，舍弃了求最小割的高效方法，方法一，而是采用了方法二。即按照index从小到大的顺序，一次判断该边是不是在最小割里面（通过再跑一次maxflow的方法）。需要注意的是，如果确定在最小割里面，可以直接在原图中把这条边去掉，再接着判断剩下的边。

/*
 ID: thestor1
 LANG: C++
 TASK: milk6
 */
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <climits>
#include <cassert>
#include <string>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <stack>
#include <algorithm>

using namespace std;

double find_path(vector<vector<double> > &capacity, const vector<vector<int> > &adjs, const int sink)
{
	int M = capacity.size();
	std::vector<bool> visited(M, false);
	// from[x] is the previous vertex visited in the shortest path from the source to x
	std::vector<int> from(M, -1);
	queue<int> Q;
	Q.push(0);
	visited[0] = true;
	int where = -1, next = -1, prev = -1;
	bool found = false;
	while (!Q.empty())
	{
		where = Q.front();
		Q.pop();
		for (int i = 0; i < adjs[where].size(); ++i)
		{
			next = adjs[where][i];
			if (!visited[next] && capacity[where][next] > 0)
			{
				Q.push(next);
				visited[next] = true;
				from[next] = where;
				if (next == sink)
				{
					found = true;
					break;
				}
			}
		}
		if (found)
		{
			break;
		}
	}
	
	// we compute the path capacity
	where = sink;
	double path_cap = -1;
	while (from[where] > -1)
	{
		prev = from[where];
		if (path_cap > -1)
		{
			path_cap = min(path_cap, capacity[prev][where]);
		}
		else
		{
			path_cap = capacity[prev][where];
		}
		where = prev;
	}
	
	// we update the residual network; if no path is found the while loop will not be entered
	where = sink;
	while (from[where] > -1)
	{
		prev = from[where];
		// if (capacity[prev][where] != INT_MAX)
		// {
        capacity[prev][where] -= path_cap;
		// }
		// if (capacity[where][prev] < INT_MAX - path_cap)
		// {
        capacity[where][prev] += path_cap;
		// }
		// else
		// {
        // assert(false);
        capacity[where][prev] = INT_MAX;
		// }
		where = prev;
	}
	
	// if no path is found, path_cap is infinity
	if (path_cap == -1)
	{
		return 0;
	}
	return path_cap;
}

void DFS(int u, vector<bool> &visited, const vector<vector<double> > &capacity, const vector<vector<int> > &adjs)
{
	visited[u] = true;
	for (int i = 0; i < adjs[u].size(); ++i)
	{
		int v = adjs[u][i];
		if (!visited[v] && capacity[u][v] > 0)
		{
			DFS(v, visited, capacity, adjs);
		}
	}
}

int gcd(int a, int b)
{
	if (b == 0)
	{
		return a;
	}
	else
	{
		return gcd(b, a % b);
	}
}

double maxflow(vector<vector<double> > capacity, const vector<vector<int> > &adjs, const int sink)
{
	double mf = 0;
	// vector<int> path_caps;
	while (true)
	{
		double path_cap = find_path(capacity, adjs, sink);
		if (path_cap == 0)
		{
			break;
		}
		else
		{
			// path_cap = (path_cap - 1) / (M + 1) * G;
			mf += path_cap;
			// path_caps.push_back(path_cap);
		}
	}
	return mf;
}

class Edge
{
public:
	int u, v, w;
    
	Edge(int u, int v, int w)
	{
		this->u = u;
		this->v = v;
		this->w = w;
	}
    
	Edge()
	{
		Edge(0, 0, 0);
	}
};

int main()
{
	ifstream fin("milk6.in");
	ofstream fout("milk6.out");
	
	int N, M;
	fin>>N>>M;
	vector<vector<double> > capacity(N, vector<double>(N, 0));
	vector<vector<int> > adjs(N, vector<int>());
	vector<Edge> edges(M);
    
	// int G = 1;
	for (int i = 0; i < M; ++i)
	{
		int u, v, w;
		fin>>u>>v>>w;
		// if (i == 0)
		// {
        // G = w;
		// }
		// else
		// {
        // G = gcd(w, G);
		// }
		
		if (capacity[u - 1][v - 1] == 0 && capacity[v - 1][u - 1] == 0)
		{
			adjs[u - 1].push_back(v - 1);
			adjs[v - 1].push_back(u - 1);
		}
		w = w * (M + 1) + 1;
		capacity[u - 1][v - 1] += w;
        
		edges[i] = Edge(u - 1, v - 1, w);
	}
    
	// for (int i = 0; i < N; ++i)
	// {
    // for (int j = 0; j < N; ++j)
    // {
    // if (capacity[i][j] > 0)
    // {
    // cout<<(INT_MAX - 1) / (M + 1) * G<<","<<capacity[i][j]<<endl;
    // assert(capacity[i][j] < (INT_MAX - 1) / (M + 1) * G);
    // capacity[i][j] = capacity[i][j] / G * (M + 1) + 1;
    // }
    // }
	// }
    
	double mf = maxflow(capacity, adjs, N - 1);
	// cout<<"[debug]mf:"<<mf<<endl;
    
	if (mf == 0)
	{
		fout<<"0 0"<<endl;
	}
	else
	{
		std::vector<int> mincut;
        double total = mf;
		for (int i = 0; i < M; ++i)
		{
			capacity[edges[i].u][edges[i].v] -= edges[i].w;
			double nmf = maxflow(capacity, adjs, N - 1);
			if (total - nmf == edges[i].w)
			{
				mincut.push_back(i + 1);
                total = nmf;
                if (total == 0)
                {
                    break;
                }
			}
			else
            {
                capacity[edges[i].u][edges[i].v] += edges[i].w;
            }
		}
		mf = (mf - mincut.size()) / (M + 1);
		fout<<(int)mf<<" "<<mincut.size()<<endl;
		for (int i = 0; i < mincut.size(); ++i)
		{
			fout<<mincut[i]<<endl;
		}
	}
    
	fin.close();
	fout.close();
	return 0;  
}