二分图匹配：匈牙利算法和最大流建模

话说一个半月以前纠结二分图匹配的匈牙利算法时想了好几天终于弄懂了，这种递归的算法设计得很巧妙，核心程序只有几行，但却有些难以理解。查了不少资料，画了不少图，终于模模糊糊懂一些了。

今天想重新热热手，于是复习一下，顺便用SAP写了一下。

一个半月前写的USACO：stall4，裸的二分图最大匹配

{
ID:lxyweb1
PROB:stall4
LANG:PASCAL
}

var
  a:array[1..1000,1..1000] of boolean;
  m,n,ans,i,j,k,x,y:longint;
  link:array[1..1000] of longint;
  flag:array[0..1000] of boolean;



function dfs(x:longint):boolean;
var
  i,j,k:longint;
begin
  for i:=1 to m do
    if (a[x,i]) and (not flag[i]) then
    begin
      flag[i]:=true;
      if (link[i]=0) or dfs(link[i]) then
      begin
        link[i]:=x;
        exit(true);
      end;
    end;
  exit(false);
end;

begin
  assign(input,'stall4.in');
  reset(input);
  assign(output,'stall4.out');
  rewrite(output);

  readln(n,m);
  for i:=1 to n do
  begin
    read(x);
    for j:=1 to x do
    begin
      read(y);
      a[i,y]:=true;
    end;
  end;
  for i:=1 to n do
  begin
    fillchar(flag,sizeof(flag),false);
    if dfs(i) then
      inc(ans);
  end;
  writeln(ans);

  close(input);
  close(output);

end.


然后想去纠结一下线性规划和网络流24题，不过好多题目都要输出方案数，是不是要Special Judge呢？然后又有些知难而退的想法。还是水平高了慢慢做吧。

01飞行员配对方案问题

var
  a:array[0..1000,0..1000] of boolean;
  i,j,k,ans,n,m,x,y:longint;
  flag:array[0..1000] of boolean;
  link:array[0..1000] of longint; 

function dfs(k:longint):boolean;
var
  i:longint;
begin
  for i:=1 to n+m do
    if a[k,i] and (not flag[i]) then
    begin
      flag[i]:=true;
      if (link[i]=0) or dfs(link[i]) then
      begin
        link[i]:=k;
        exit(true);
      end;
    end;
  exit(false);
end; 

begin
  assign(input,'air.in');
  reset(input);
  assign(output,'air.out');
  rewrite(output); 

  readln(n,m);
  while not eof do
  begin
    readln(x,y);
    if x=-1 then break;
    a[x,y]:=true;
  end;
  for i:=1 to n do
  begin
    fillchar(Flag,sizeof(flag),0);
    if dfs(i) then
      inc(ans);
  end;
  writeln(ans);
  for i:=1 to n do
    for j:=n+1 to n+m do
      if link[j]=i then
      writeln(i,' ',j); 

  close(input);
  close(output);
end.

我觉得link数组对应的就是相连的外国飞行员吧，没有测试，也不知道对不对。

另外，还用SAP写了一下stall4，发现还不是很熟悉，4天前刚学现在又要看着源代码敲了。不过基本思路是有的。很傻得敲错了一个变量，又对照着DEBUG好久。

关于建模：

在二分图的基础上增加源S和汇T。
1、S向X集合中每个顶点连一条容量为1的有向边。
2、Y集合中每个顶点向T连一条容量为1的有向边。
3、XY集合之间的边都设为从A集合中的点到B集合之中的点，容量为1的有向边。

求网络最大流，流量就是匹配数，所有满流边是一组可行解。

{
  ID:lxyweb1
  PROB:stall4
  LANG:PASCAL
}
const
  oo=19930508;
var
  size,i,j,k,m,n,flow,x,y:longint;
  e,next,opp,c:array[0..1000000] of longint;
  g,h,vh:Array[0..10000] of longint;
  s,t,tot:longint;

procedure insert(x,y,z:longint);
begin
  inc(size);  e[size]:=y; next[size]:=g[x];
  c[size]:=z; g[x]:=size; opp[size]:=size+1;
  inc(size);  e[size]:=x; next[size]:=g[y];
  c[size]:=0; g[y]:=size; opp[size]:=size-1;
end;

function sap(i,now:longint):longint;
var
  j,p,tmp,minh:longint;
begin
  minh:=tot-1;
  sap:=0;
  if i=t then
  begin
    inc(flow,now);
    exit(now);
  end;
  p:=g[i];
  while p>0 do
  begin
    j:=e[p];
    if c[p]>0 then
    begin
      if h[j]+1=h[i] then
      begin
        if c[p]<now then
          tmp:=sap(j,c[p])
        else tmp:=sap(j,now);
        dec(now,tmp);
        inc(sap,tmp);
        dec(c[p],tmp);
        inc(c[opp[p]],tmp);
        if h[s]>=tot then
          exit;
        if now=0 then break;
      end;
      if h[j]<minh then minh:=h[j];
    end;
    p:=next[p];
  end;
  if sap=0 then
  begin
    dec(vh[h[i]]);
    if vh[h[i]]=0 then h[s]:=tot;
    h[i]:=minh+1;
    inc(vh[h[i]]);
  end;
end;


procedure init;
begin
  readln(n,m);
  for i:=1 to n do
  begin
    read(x);
    for j:=1 to x do
    begin
      read(y);
      insert(i,y+n,1);
    end;
  end;
  s:=n+m+1;
  t:=s+1;
  tot:=t;
  for i:=1 to n do
    insert(s,i,1);
  for i:=n+1 to n+m do
    insert(i,t,1);
end;
procedure main;
begin
  vh[0]:=tot;
  while h[s]<tot do
    sap(s,oo);
end;
procedure print;
begin
  writeln(flow);
end;


begin
  assign(input,'stall4.in');
  reset(input);
  assign(output,'stall4.out');
  rewrite(output);

  init;
  main;
  print;

  close(input);
  close(output);
end.

总结一下，最大流的步骤：

1.输入，构图，一般用边集数组的话AddEdge（包括加反向边）有12个语句。

2.SAP。while h[s]<tot do sap(s,oo)

procedure sap(i,now)

minh=tot-1

sap=0

if i=t更新，返回

对每条边：

  如果是可行的，那么DFS，更新当前流量和调整边，如果h[s]>=tot则返回，当前流为0则跳出

  根据每条邻接的边调整minh

如果没有可行流则调整，4条语句