实用算法实现-第 23 篇 最大流

23.1    Ford-Fullkerson方法

使用BFS来实现Ford-Fullkerson方法中的找增广路径的算法称为Edmods-Karp算法。Edmods-Karp算法是最短增广路算法，因为实用BFS找到的增广路径是所有可能的增广路径中最短的路径。它的复杂度是O(VE­­²­)，其中V是结点数，E是有向边数。

如果用使用DFS代替BFS，则Ford-Fullkerson方法退化成一般增广路算法。其复杂度是O(E| f* |)。其中f*是算法找出的最大流。

23.2    最大流的Edmods-Karp算法

23.2.1   实例

PKU JudgeOnline, 1273, Drainage Ditches.

23.2.2   问题描述

给出M个点然后给出N条边的最大流量，求最大流。

23.2.3   输入

54

12 40

14 20

24 20

23 30

3 4 10

23.2.4   输出

50

23.2.5   分析

基本的最大流算法Edmonds-Karp就可以完成。

PKU JudgeOnline, 3436, ACM Computer Factory也是基本的最大流问题，只是构图需要花一点时间。

23.2.6   程序

#include<stdio.h>
#include<string.h>
#include<iostream>
using namespace std;
#define MAX   0x7FFFFFFF
int route[201][201];
int c[201][201];
int f[201][201];
struct queue{
     int parent;
     int node;
};
queuequeue[2000];
int last;
int visited[201];
int M;
void EdmondsKarp()
{
     int i, j;
     int temp;
     int from,to;
     intpointer;
     int minC;
     int parent;
     while(1){
         memset(visited, 0, sizeof(visited));
         last = 0;
         queue[last].parent = 0;
         queue[last++].node = 1;
         pointer = -1;
         for(i =0; i < last; i++){
              //BFS，找增广路径
              temp = queue[i].node;
              for(j= 1; j <= M; j++){
                   if((visited[j]== 0)&&
                       (c[temp][j] != 0))
                   {
                       visited[j] = 1;
                       queue[last].node = j;
                       queue[last++].parent = i;
                       if(j == M){
                            //到达汇结点，找到一条增广路径
                            pointer = last - 1;
                            break;
                       }
                   }
              }
              if(pointer!= -1)
              {
                   break;
              }
         }
         if(pointer== -1)
         {
              //没有找到增广路径
              break;
         }
         minC = MAX;
         temp = pointer;
         while(temp!= 0){
              //计算增光路径的流
              parent = queue[temp].parent;
              from = queue[parent].node;
              to = queue[temp].node;
              if(minC> c[from][to]){
                   minC = c[from][to];
              }
              temp = parent;
         }
         temp = pointer;
         while(temp!= 0){
              //修改残留网络
              parent = queue[temp].parent;
              from = queue[parent].node;
              to = queue[temp].node;
              c[from][to] = c[from][to] - minC;
              c[to][from] = c[to][from] + minC;
              f[from][to] = f[from][to] + minC;
              f[to][from] = -f[from][to];
              temp = parent;
         }
     }
}
int main()
{
     int N;
     int i;
     int temp;
     int from,to;
     int flow;
     while(cin>> N >> M){
         memset(c, 0, sizeof(c));
         memset(route, 0, sizeof(route));
         memset(f, 0, sizeof(f));
         for(i =0; i < N; i++){
              cin >> from >> to;
              cin >> temp;
              route[from][to] += temp;
              c[from][to] = route[from][to];
         }
         EdmondsKarp();
         flow = 0;
         for(i =2; i <= M; i++)
         {
              flow += f[1][i];
         }
         cout << flow << endl;
     }
}

23.3    二分图的最小点权覆盖集算法

最小点权覆盖集（minimum weight vertexcovering set，MinWVCS）是在带点权无向图G中，点权之和最小的点覆盖集。

二分图的最小点权覆盖集问题可以转化成最小割问题。以下图为例：

图一

 

图二

    图一中的二分图的最小点权覆盖问题可以通过如上的方法转换为最小割问题。其中不和源或者汇连接的边的权值为∞。具体的证明参见《最小割模型在信息学竞赛中的应用》。

23.3.1   实例

PKU JudgeOnline, 3308, Paratroopers.

23.3.2   问题描述

有外星人入侵地球。已知外星人入侵的坐标，要将他们全部消灭。在横向（或者纵向）上可以布置武器，横向的武器能把整行的外星人消灭。横向（或者纵向）上每个坐标布置武器的花销已知。总代价是所有武器费用的乘积。求最小的总代价。

首先输入测试数据个数T。第二行输入行数、列数、外星人总数。接着一行输入每行代价。接着一行输入每列代价。接着几行输入外星人着陆的行号和列号。

输出总的代价，保留四位小数点。

23.3.3   输入

1

44 5

2.07.0 5.0 2.0

1.52.0 2.0 8.0

11

22

33

44

1 4

23.3.4   输出

16.0000

23.3.5   分析

这个问题和PKU JudgeOnline, 3041,Asteroids.比较相似。只不过3041是求最小点覆盖，这里每个点还带上了权。

这个问题很容易地归为二分图的最小点权覆盖问题。首先问题可以转化为对每个外星人可以从横向也可以从纵向上杀死，但是一定要杀死，这就是覆盖问题。然而需要注意的是，这里要求的是最小的武器费用的乘积，而最小点权覆盖问题求的是最小的权之和。将每个武器的费用取对数，问题就变成了要求的是最小的武器费用（取对数后）之和。

       构图完成之后直接采用最大流算法完成即可。

PKU JudgeOnline, 2125, Destroying TheGraph和这个题目很像。不过更困难的是2125还需要找出最小割的边。找最小割的边的方法是DFS遍历残余图。《最小割模型在信息学竞赛中的应用》中在介绍二分图的最小点权覆盖时进行了一些分析。由于最小割可能有多种情况，故此该问题需要Special Judge。

23.3.6   程序

#include<stdio.h>
#include<string.h>
#include<iostream>
#include<stdlib.h>
#include<math.h>
using namespace std;
#define MAX   0x7FFFFFFF
#define MIN   0.0000001
double route[201][201];
double c[201][201];
double f[201][201];
struct queue{
     int parent;
     int node;
};
queuequeue[2000];
int last;
int visited[201];
int M;
void EdmondsKarp()
{
     int i, j;
     int temp;
     int from,to;
     intpointer;
     doubleminC;
     int parent;
     while(1){
         memset(visited, 0, sizeof(visited));
         last = 0;
         queue[last].parent = 0;
         queue[last++].node = 0;
         pointer = -1;
         for(i =0; i < last; i++){
              temp = queue[i].node;
              for(j= 0; j <= M; j++){
                   if((visited[j]== 0)&&
                   (c[temp][j] > MIN))
                   //   (c[temp][j] != 0))
                   {
                       visited[j] = 1;
                       queue[last].node = j;
                       queue[last++].parent = i;
                       if(j== M){
                            pointer = last - 1;
                            break;
                       }
                   }
              }
              if(pointer!= -1)
              {
                   break;
              }
         }
         if(pointer== -1)
         {
              break;
         }
         minC = MAX;
         temp = pointer;
         while(temp!= 0){
              parent = queue[temp].parent;
              from = queue[parent].node;
              to = queue[temp].node;
              if(minC> c[from][to]){
                   minC = c[from][to];
              }
              temp = parent;
         }
         temp = pointer;
         while(temp!= 0){
              parent = queue[temp].parent;
              from = queue[parent].node;
              to = queue[temp].node;
              c[from][to] = c[from][to] - minC;
              c[to][from] = c[to][from] + minC;
              f[from][to] = f[from][to] + minC;
              f[to][from] = -f[from][to];
              temp = parent;
         }
     }
}
struct weight{
     doublecost;
     doublelogCost;
};
int main()
{
     int cases;
     int rowNum,columnNum;
     int i, j;
     int l;
     weight row[51];
     weight column[51];
     int from;
     int to; 
     doubleflow;
 
     cin >> cases;
     for(; cases> 0; cases--){
         memset(c, 0, sizeof(c));
         memset(route, 0, sizeof(route));
         memset(f, 0, sizeof(f));
         cin >> rowNum >> columnNum>> l;
         for(i =1; i <= rowNum; i++){
              cin >> row[i].cost;
              row[i].logCost =log10(row[i].cost);
 
              route[0][i] = row[i].logCost;
              c[0][i] = route[0][i];
         }
         for(i =1; i <= columnNum; i++){
              cin >> column[i].cost;
              column[i].logCost =log10(column[i].cost);
 
              from = i + rowNum;
              to = rowNum + columnNum + 1;
              route[from][to] =column[i].logCost;
              c[from][to] = route[from][to];
         }
         M = rowNum + columnNum + 1;
         for(i =0; i < l; i++){
              cin >> from >> to;
              //from =from;
              to = to + rowNum;
              route[from][to] = MAX;
              c[from][to] = route[from][to];
         }
         EdmondsKarp();
         flow = 0;
         for(i =1; i <= M; i++)
         {
              flow += f[0][i];
         }
 
         flow = pow(10.0, flow);
         printf("%.4f\n",flow);
     }
}
 

23.4    最大二分匹配

可以证明最大二分匹配问题可以转化为最大流问题。求最大二分匹配实际上就是求最大流。

二分图的最大匹配有两种求法，第一种是最大流（我在此假设读者已有网络流的知识）。第二种就是我现在要讲的匈牙利算法。这个算法说白了就是最大流的算法，但是它跟据二分图匹配这个问题的特点，把最大流算法做了简化，提高了效率。匈牙利算法其实很简单，但是网上搜不到什么说得清楚的文章。所以我决定要写一下。

最大流算法的核心问题就是找增广路径（augment path）。匈牙利算法也不例外，它的基本模式就是：

 

初始时最大匹配为空

while 找得到增广路径

do 把增广路径加入到最大匹配中去

 

可见和最大流算法是一样的。但是这里的增广路径就有它一定的特殊性，下面我来分析一下。

（注：匈牙利算法虽然根本上是最大流算法，但是它不需要建网络模型，所以图中不再需要源点和汇点，仅仅是一个二分图。每条边也不需要有方向。）

图1 图2

图1是我给出的二分图中的一个匹配：［1，5］和［2，6］。图2就是在这个匹配的基础上找到的一条增广路径：3->6->2->5->1->4。我们借由它来描述一下二分图中的增广路径的性质：

(1)有奇数条边。

(2)起点在二分图的左半边，终点在右半边。

(3)路径上的点一定是一个在左半边，一个在右半边，交替出现。（其实二分图的性质就决定了这一点，因为二分图同一边的点之间没有边相连，不要忘记哦。）

(4)整条路径上没有重复的点。

(5)起点和终点都是目前还没有配对的点，而其它所有点都是已经配好对的。（如图1、图2所示，［1，5］和［2，6］在图1中是两对已经配好对的点。而起点3和终点4目前还没有与其它点配对。）

(6)路径上的所有第奇数条边都不在原匹配中，所有第偶数条边都出现在原匹配中。（如图1、图2所示，原有的匹配是［1，5］和［2，6］，这两条配匹的边在图2给出的增广路径中分边是第2和第4条边。而增广路径的第1、3、5条边都没有出现在图1给出的匹配中。）

(7)最后，也是最重要的一条，把增广路径上的所有第奇数条边加入到原匹配中去，并把增广路径中的所有第偶数条边从原匹配中删除（这个操作称为增广路径的取反），则新的匹配数就比原匹配数增加了1个。（如图2所示，新的匹配就是所有蓝色的边，而所有红色的边则从原匹配中删除。则新的匹配数为3。）

不难想通，在最初始时，还没有任何匹配时，图1中的两条灰色的边本身也是增广路径。因此在这张二分图中寻找最大配匹的过程可能如下：

(1)找到增广路径1->5，把它取反，则匹配数增加到1。

(2)找到增广路径2->6，把它取反，则匹配数增加到2。

(3)找到增广路径3->6->2->5->1->4，把它取反，则匹配数增加到3。

(4)再也找不到增广路径，结束。

当然，这只是一种可能的流程。也可能有别的找增广路径的顺序，或者找到不同的增广路径，最终的匹配方案也可能不一样。但是最大匹配数一定都是相同的。

对于增广路径还可以用一个递归的方法来描述。这个描述不一定最准确，但是它揭示了寻找增广路径的一般方法：

“从点A出发的增广路径”一定首先连向一个在原匹配中没有与点A配对的点B。如果点B在原匹配中没有与任何点配对，则它就是这条增广路径的终点。反之，如果点B已与点C配对，那么这条增广路径就是从A到B，再从B到C，再加上“从点C出发的增广路径”。并且，这条从C出发的增广路径中不能与前半部分的增广路径有重复的点。

比如图2中，我们要寻找一条从3出发的增广路径，要做以下3步：

(1)首先从3出发，它能连到的点只有6，而6在图1中已经与2配对，所以目前的增广路径就是3->6->2再加上从2出发的增广路径。

(2)从2出发，它能连到的不与前半部分路径重复的点只有5，而且5确实在原匹配中没有与2配对。所以从2连到5。但5在图1中已经与1配对，所以目前的增广路径为3->6->2->5->1再加上从1出发的增广路径。

(3)从1出发，能连到的不与自已配对并且不与前半部分路径重复的点只有4。因为4在图1中没有与任何点配对，所以它就是终点。所以最终的增广路径是3->6->2->5->1->4。

但是严格地说，以上过程中从2出发的增广路径（2->5->1->4）和从1出发的增广路径（1->4）并不是真正的增广路径。因为它们不符合前面讲过的增广路径的第5条性质，它们的起点都是已经配过对的点。我们在这里称它们为“增广路径”只是为了方便说明整个搜寻的过程。而这两条路径本身只能算是两个不为外界所知的子过程的返回结果。

显然，从上面的例子可以看出，搜寻增广路径的方法就是DFS，可以写成一个递归函数。当然，用BFS也完全可以实现。

至此，理论基础部份讲完了。但是要完成匈牙利算法，还需要一个重要的结论：

定理：如果从一个点A出发，没有找到增广路径，那么无论再从别的点出发找到多少增广路径来改变现在的匹配，从A出发都永远找不到增广路径。

要用文字来证明这个定理很繁，话很难说，要么我还得多画一张图，我在此就省了。其实你自己画几个图，试图举两个反例，这个定理不难想通的。（给个提示。如果你试图举个反例来说明在找到了别的增广路径并改变了现有的匹配后，从A出发就能找到增广路径。那么，在这种情况下，肯定在找到别的增广路径之前，就能从A出发找到增广路径。这就与假设矛盾了。）

有了这个定理，匈牙利算法就成形了。如下：

 

初始时最大匹配为空

for 二分图左半边的每个点i

do 从点i出发寻找增广路径。如果找到，则把它取反（即增加了总了匹配数）。

 

如果二分图的左半边一共有n个点，那么最多找n条增广路径。如果图中共有m条边，那么每找一条增广路径（DFS或BFS）时最多把所有边遍历一遍，所花时间也就是m。所以总的时间大概就是O（n * m）。

可以转化为最大二分匹配的问题有以下几种：

1. 最小路径覆盖。有向图n个点，求最少的不相交（包括顶点相交）路径，覆盖所有n个结点，输出路径条数。

2. 最小点覆盖。最小覆盖要求用最少的点（X集合或Y集合的都行）让每条边都至少和其中一个点关联。

23.4.1   最小路径覆盖

在一个P*P的有向图中，路径覆盖就是在图中找一些路经，使之覆盖了图中的所有顶点，且任何一个顶点有且只有一条路径与之关联。（如果把这些路径中的每条路径从它的起始点走到它的终点，那么恰好可以经过图中的每个顶点一次且仅一次）。如果不考虑图中存在回路，那么每每条路径就是一个弱连通子集。

 

由上面可以得出：

1. 一个单独的顶点是一条路径。

2. 如果存在一路径p1,p2,......pk，其中p1 为起点，pk为终点，那么在覆盖图中，顶点p1,p2,......pk不再与其它的顶点之间存在有向边。

最小路径覆盖就是找出最小的路径条数，使之成为P的一个路径覆盖。

路径覆盖与二分图匹配的关系：

最小路径覆盖=|P|-最大匹配数。

其中最大匹配数的求法是把P中的每个顶点pi分成两个顶点pi'与pi''，如果在p中存在一条pi到pj的边，那么在二分图P＇中就有一条连接pi'与pj''的无向边。这里pi' 就是p中pi的出边，pj''就是p中pj 的一条入边。

对于公式：最小路径覆盖=|P|-最大匹配数。可以这么来理解：

如果匹配数为零，那么P中不存在有向边，于是显然有：

最小路径覆盖=|P|-最大匹配数=|P|-0=|P|。即P的最小路径覆盖数为|P|。

P＇中不在于匹配边时，路径覆盖数为|P|。

如果在P＇中增加一条匹配边pi'-->pj''，那么在图P的路径覆盖中就存在一条由pi连接pj的边，也就是说pi与pj 在一条路径上，于是路径覆盖数就可以减少一个。

如此继续增加匹配边，每增加一条，路径覆盖数就减少一条。直到匹配边不能继续增加时，路径覆盖数也不能再减少了，此时就有了前面的公式。但是这里只是说话了每条匹配边对应于路径覆盖中的一条路径上的一条连接两个点之间的有向边。下面来说明一个路径覆盖中的每条连接两个顶点之间的有向边对应于一条匹配边。

与前面类似，对于路径覆盖中的每条连接两个顶点之间的每条有向边pi--->pj，我们可以在匹配图中对应做一条连接pi'与pj''的边，显然这样做出来图的是一个匹配图（这一点用反证法很容易证明，如果得到的图不是一个匹配图，那么这个图中必定存在这样两条边  pi'---pj'' 及 pi' ----pk''，（j!= k），那么在路径覆盖图中就存在了两条边pi-->pj,pi--->pk ，那边从pi出发的路径就不止一条了，这与路径覆盖图是矛盾的。还有另外一种情况就是存在pi'---pj'',pk'---pj''，这种情况也类似可证）。

至此，就说明了匹配边与路径覆盖图中连接两顶点之间边的一一对应关系，那么也就说明了前面的公式成立。

1.4.2   最小点覆盖

König定理是一个二分图中很重要的定理。

König定理：一个二分图中的最大匹配数等于这个图中的最小点覆盖数。

可以证明：最少的点（即覆盖数）= 最大匹配数。

可以比较简单地想：如果有一条边没有被覆盖，也就是有一条边没有和一个覆盖点连接，则可以通过这条边流过流量，也就是这条边的两个端点可以匹配起来作为最大匹配的一对匹配。

 

23.4.3   实例

PKU JudgeOnline, 3041, Asteroids.

23.4.4   问题描述

在N*N的格子里有K个小行星。可以用激光炮把一行或者一列的小行星清除。问最少要发射多少次激光炮。

23.4.5   输入

34

11

13

22

32

23.4.6   输出

2

23.4.7   分析

本题就是要找最小点覆盖，根据König 定理就是求最大二分匹配了。PKU JudgeOnline, 2226, Muddy Fields是这个题目的威力加强版，构图的时候相对麻烦一点。

PKU JudgeOnline, 2771, Guardian of Decency 也是一个找最小点覆盖的问题，可以这样理解：将男女分成两队，能恋爱的连接起来。要求去掉最少的人，让所有的连接断掉。对比上面最小点覆盖的问题，可以发现这个问题就转换成了最小点覆盖问题。

PKU JudgeOnline, 1325,Machine Schedule也是一个找最小点覆盖的问题。

PKU JudgeOnline, 1422, Air Raid是一个最小路径覆盖问题。

23.4.8   程序

#include<string.h>
#include<stdio.h>
#include<iostream.h>
#define MAX_NUM 502
int used[MAX_NUM];   //记录R中结点是否使用
int link[MAX_NUM];   //记录当前与R结点相连的L的结点
int connected[MAX_NUM][MAX_NUM]; //记录连接L和R的边，如果i和j之间有边则为，否则为
//int gn,gm;    //二分图中L和R中点的数目
int numL, numR;
int augmentPath(int t)
{
     int i;
     for(i = 1;i <= numR;i++)
     {
         if(used[i]==0&& connected[t][i])
         {
              used[i]=1;
              if(link[i]== -1 || augmentPath(link[i]))
              {
                   link[i]=t;
                   return1;
              }
         }
     }
     return 0;
}
int MaxMatch()
{
     int i, num;
     num = 0;
     memset(link, 0xff, sizeof(link));
     for(i = 1;i <= numL; i++){
         memset(used, 0, sizeof(used));
         if(augmentPath(i))
         {
              num++;
         }
     }
     return num;
}
int main()
{
     int K, N;
     int i;
     int row,column;
     cin >> N >> K;
     numL = K;
     numR = K;
     memset(connected, 0, sizeof(connected));
     for(i = 1;i < K; i++){
         scanf("%d",&row, &column);
         connected[row][column] = 1;
     }
     cout << MaxMatch();
     return 1;
}

23.5    实例

23.5.1   最大流实例

PKU JudgeOnline, 1273, Drainage Ditches.

PKU JudgeOnline, 3436, ACM Computer Factory.

23.5.2   二分图的最小点权覆盖集算法实例

PKU JudgeOnline, 2125, Destroying TheGraph.

PKU JudgeOnline, 3308, Paratroopers.

 23.5.3   最大二分匹配实例

PKU JudgeOnline, 3041, Asteroids.

PKU JudgeOnline, 2771, Guardian of Decency.

PKU JudgeOnline, 2536, Gopher II.

PKU JudgeOnline, 1274, The Perfect Stall.

PKU JudgeOnline, 2239, Selecting Courses.

PKU JudgeOnline, 1469, COURSES.

PKU JudgeOnline, 1422, Air Raid.

PKU JudgeOnline, 1325, Machine Schedule.

PKU JudgeOnline, 2226, Muddy Fields.

本文章欢迎转载，请保留原始博客链接http://blog.youkuaiyun.com/fsdev/article