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;
}