匈牙利算法的优化 - Hoblovski's Blog - 想拿Ag的蒟蒻.已经Ag滚出.
自动命名Python脚本

匈牙利算法的优化

Hoblovski posted @ 2014年7月12日 00:35 in Discussion with tags 匈牙利算法 网络流 二分图匹配 丧心病狂 , 1180 阅读
对于二分图匹配一类问题,我们有通用的网络流算法以及简单好写的匈牙利算法.
以下是匈牙利的核心代码(假定二分图|X|=|Y|)(直接txt中手打可能CE危险)
function augment(u:longint):boolean;
var t1:longint;
begin
v[u]:=true; t1:=g[u]; while t1<>0 do begin
	if (pair[mem[t1].n]=0)or((not v[pair[mem[t1].n]])and(augment(pair[mem[t1].n])) then begin
		pair[mem[t1].n]:=u; exit(true);
	end; t1:=mem[t1].next;
end; exit(false);
end;
function hungary:longint;
var i:longint;
begin
hungary:=0; fillchar(pair,sizeof(pair),0); for i:=1 to n do begin
        fillchar(v,sizeof(v),0); inc(hungary,byte(augment(i)));
end;  
end;

网络流算法代码略(参见文末),蒟蒻使用的是ISAP,静态内存邻接表(链式前向星).

 

不过令人惊异的是虽然匈牙利算法的常数很小,在小数据可以完胜网络流算法

但是在数据大了之后我们却要使用网络流算法(SCOI 游戏的匹配做法,如果要用匈牙利需要特殊优化).

我们估计,匈牙利算法在这种情况下的制限因素主要是
for i:=1 to n do begin
        fillchar(v,sizeof(v),0); inc(hungary,byte(augment(i)));
end;
中的fillchar,SCOI那道题也是需要根据题目特性改成 fillchar(v,i,0);
由于每次访问的点很少而又要全部fillchar一次显然很浪费,不妨把访问到的点插入队列之后手动改访问标志.
我试着改了改fillchar优化,但是
优化之后时间没有任何提高,反而有所下降,观察数据,
在40000个点160000条边随机图中平均匈牙利算法每次增广会访问1000个点左右
这时候的fillchar优化毫无作用.
 
蒟蒻不甘心,脑洞大开试着改了改匈牙利算法中语句的实行顺序,augment函数改为(同样txt手打小心CE)
 
function augment(u:longint):boolean;
var t1:longint;
begin
v[u]:=true; t1:=g[u]; while t1<>0 do begin
	if (pair[mem[t1].n]=0) then begin pair[mem[t1].n]:=u; exit(true); end; t1:=mem[t1].next;
end; t1:=g[u]; while t1<>0 do begin
	if (not v[pair[mem[t1].n]])and(augment(pair[mem[t1].n])) then begin
		pair[mem[t1].n]:=u; exit(true);
	end; t1:=mem[t1].next;
end;exit(false); 
匈牙利算法的速度提高了一倍!,每次增广只访问300个点左右了.
不过匈牙利在大数据上仍然跪给ISAP.
还能有更多优化么?随机打乱边?启发式安排边的顺序?
欢迎神牛讨(虐)论(场)

附 效率测试
生成数据无重边,全随机,匈牙利的数组大小统统是N
数据             ISAP   匈牙利(优化完全)
N=10000,M=10000  20ms   10ms
N=10000,M=100000 80ms   50ms
N=20000,M=60000  130ms  230ms
N=20000,M=100000 200ms  200ms
N=40000,M=80000  320ms  110ms
N=40000,M=160000 330ms  >1500ms
 
下附代码,Hungary自带增广访问结点平均计数(已测对速度影响很小)
time_stl是我自己写的一个库函数通过调用dos库中的gettime来达到时间计数
详情请自学Pascal gettime函数以及Pascal 库函数编写www
因为刚写完bzoj2150这道最小路径覆盖所以文件名什么的才会是bzoj2150
 
 
program bzoj2150;

uses time_stl;

type tnode=record
        n,c,next:longint;
     end;

const maxn=100017;
      maxm=1000017;
      maxint=longint($3f3f3f3f);

var g,d,dc,cur,pre,arc,dlt:array[0..maxn] of longint;
    mem:array[0..maxm] of tnode;
    n,m,memsize,st,trm:longint;

function min(i,j:longint):longint; inline; begin
if i<j then exit(i) else exit(j);  end;

procedure insnbs(u,v,ic:longint);
begin
inc(memsize); with mem[memsize] do begin
        n:=v; c:=ic; next:=g[u];
end; g[u]:=memsize;
inc(memsize); with mem[memsize] do begin
        n:=u; c:=0; next:=g[v];
end; g[v]:=memsize;
end;

function isap:longint;
var u,w,t1:longint;
begin
fillchar(d,sizeof(d),0); fillchar(dc,sizeof(dc),0);
dc[0]:=n; pre[st]:=st; dlt[st]:=maxint; cur:=g; u:=st; isap:=0;
while d[st]<n do begin
        while (cur[u]<>0)and((mem[cur[u]].c=0)or(d[u]<>d[mem[cur[u]].n]+1)) do cur[u]:=mem[cur[u]].next;
        if cur[u]=0 then begin
                dec(dc[d[u]]); if dc[d[u]]=0 then break; w:=n; t1:=g[u]; while t1<>0 do begin
                        if (mem[t1].c>0)and(d[mem[t1].n]<w) then w:=d[mem[t1].n]; t1:=mem[t1].next;
                end; d[u]:=w+1; inc(dc[d[u]]); cur[u]:=g[u]; u:=pre[u];
        end else begin
                pre[mem[cur[u]].n]:=u;
                arc[mem[cur[u]].n]:=cur[u];
                dlt[mem[cur[u]].n]:=min(dlt[u],mem[cur[u]].c);
                u:=mem[cur[u]].n; if u=trm then begin
                        w:=dlt[trm]; inc(isap,w); while u<>st do begin
                                dec(mem[arc[u]].c,w); inc(mem[arc[u]xor 1].c,w); u:=pre[u];
                        end;
                end;
        end
end;
end;

procedure init;
var i,j,k,u,v:longint;
    ch:char;
begin
fillchar(g,sizeof(g),0); memsize:=1; readln(n,m); st:=n<<1+1; trm:=st+1;
for i:=1 to m do begin readln(u,v);    insnbs(u,v+n,1);   end;
for i:=1 to n do begin insnbs(st,i,1); insnbs(i+n,trm,1); end;
n:=n<<1+2;
end;

begin
assign(input,'bzoj2150.in'); reset(input);
assign(output,''); rewrite(output);
init;
writeln(isap);
close(input); close(output);
end.
program bzoj2150;

uses time_stl;

type tnode=record
        n,next:longint;
     end;

const maxn=40017;
      maxm=1000017;

var g,pair,q:array[0..maxn] of longint;
    v:array[0..maxn] of boolean;
    mem:array[0..maxm] of tnode;
    n,m,memsize:longint;

procedure insnbs(u,v:longint);
begin
inc(memsize); with mem[memsize] do begin
        n:=v; next:=g[u];
end; g[u]:=memsize;
end;

function augment(u:longint):boolean;
var t1:longint;
begin
inc(q[0]); v[u]:=true; t1:=g[u]; while t1<>0 do begin
        if (pair[mem[t1].n]=0) then begin pair[mem[t1].n]:=u; exit(true); end; t1:=mem[t1].next;
end; t1:=g[u]; while t1<>0 do begin
        if (not v[pair[mem[t1].n]])and(augment(pair[mem[t1].n])) then begin pair[mem[t1].n]:=u; exit(true); end; t1:=mem[t1].next;
end; exit(false);
end;

function hungary:longint;
var i,j:longint;  t:extended;
begin
hungary:=0; fillchar(pair,sizeof(pair),0);  t:=0;
for i:=1 to n do begin
        fillchar(v,sizeof(v),0);   t:=t+q[0]; q[0]:=0;
        inc(hungary,byte(augment(i)));
end;  writeln((t/n):0:20);
end;

procedure init;
var i,j,u,v:longint;
begin
fillchar(g,sizeof(g),0); memsize:=0;
readln(n,m); for i:=1 to m do begin
        readln(u,v); insnbs(u,v);
end;
end;

begin
assign(input,'bzoj2150.in'); reset(input);
assign(output,''); rewrite(output);
init;
writeln(hungary);
close(input); close(output);
end.

 


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter