LocalImprove算法

这里简单的实现了一下soda14-Flow-Based Algorithms for Local Graph Clustering上说到的LocalImprove算法。

距离的论文链接在这里：http://arxiv.org/pdf/1307.2855.pdf

算法的伪代码如下：

#include <bits/stdc++.h>
#define maxn 100000
#define maxm 3000000
#define eps 1e-6
#define INF 1<<30
using namespace std;

struct Edge{
int u;          //u表示边的起点
int v;          //v表示边的终点
int next;       //next表示下一条边编号
double w;       //w表示边的容量
}edge[maxm];        //网络流建模时用到的边

int e;                  //边的计数器
int S;                  //超级源点
int T;                  //超级汇点
int first[maxn];        //图邻接表头结点
int d[maxn];            //层次图距离标记
int work[maxn];         //dinic优化
int q[maxn];            //bfs队列
int deg[maxn];          //每个点的度数
set<int> A;             //A集合中的点
set<int> Bs;            //Bs集合中的点
vector<int> G[maxn];    //原始图上的边

int n;                  //图上的点数
int m;                  //图上的边数
int k;                  //A集合中点的数量
int vol_a;              //A集合的容量vol(A)
double epsilon_sigma;   //必要的参数ε(σ)即ε,文中定义ε=1/(3*(1/σ-1))
double sigma;           //必要的参数σ,可推导出σ=3ε/(3ε+1)

/*图的初始化*/
void init(){
e = 0;
memset(first,-1,sizeof(first));
}

/*加边*/
void add_edge(int a,int b,double c){
//printf("add_edge:from %d to %d,val = %.4f\n",a,b,c);
edge[e].u = a;
edge[e].v = b;
edge[e].next = first[a];
edge[e].w = c;
first[a] = e++;

edge[e].u = b;
edge[e].v = a;
edge[e].next = first[b];
edge[e].w = 0;
first[b] = e++;
}

/*bfs构造层次图*/
int bfs(){
int rear = 0;
memset(d,-1,sizeof(d));
d[S] = 0;q[rear++] = S;
for(int i = 0;i < rear;i++){
for(int j = first[q[i]];j != -1;j = edge[j].next){
int to = edge[j].v;
double val = edge[j].w;
if(abs(val) > eps && d[to] == -1){
d[to] = d[q[i]]+1;
q[rear++] = to;
if(to == T) return 1;
}
}
}
return 0;
}

/*dfs计算阻塞流*/
double dfs(int cur,double a){
if(cur == T)    return a;
for(int &i = work[cur];i != -1;i = edge[i].next){
int to = edge[i].v;
double val = edge[i].w;
if(abs(val) > eps && d[to] == d[cur]+1){
if(double t = dfs(to,min(a,val))){
edge[i].w -= t;edge[i^1].w += t;
return t;
}
}
}
return 0;
}

bool f1[maxn],f2[maxn];

void dfs1(int u){
f1[u] = 1;
for(int i = first[u];i != -1;i = edge[i].next){
int to = edge[i].v;
double val = edge[i].w;
if(!f1[to] && val > eps)    dfs1(to);
}
}

void dfs2(int u){
f2[u] = 1;
for(int i = first[u];i != -1;i = edge[i].next){
int to = edge[i].v;
double val = edge[i^1].w;
if(!f2[to] && val > eps)    dfs2(to);
}
}

//min_cut输出最小割，并给出当前划分的导率conductance
double min_cut(){
/*for(int i = 0;i < e;i += 2){
printf("edge:from %d to %d,val = %.4f\n",edge[i].u,edge[i].v,edge[i].w);
}*/
memset(f1,0,sizeof(f1));
memset(f2,0,sizeof(f2));
dfs1(S);
dfs2(T);

int cnt = 0;    //表示跨越集合的边数

//printf("min_cut edge:\n");
for(int i = 0;i < e;i += 2){
if(f1[edge[i].u] && f2[edge[i].v] && edge[i].w < eps){
//printf("from %d to %d\n",edge[i].u,edge[i].v);
if(edge[i].u == S || edge[i].u == T)    continue;
if(edge[i].v == S || edge[i].v == T)    continue;
cnt++;
}
}
printf("\n");

int vol1 = 0;
//set<int> S_new;//新的划分中S中的点
for(int i = S+1;i < T;i++){
if(f1[i]){
vol1 += deg[i];
//printf("flag1:%d\n",i);
}
}
/*set<int>::iterator iter;
for(iter = S_new.begin();iter != S_new.end();iter++){
int u = *iter;
for(int i = 0;i < G[u].size();i++){
int v = G[u][i];
if(S_new.find(v) != S_new.end())    continue;
cnt++;
}
}*/
printf("cnt = %d,vol1 = %d,total-vol1 = %d\n",cnt,vol1,2*m-vol1);
if(cnt == 0){
return 0;
}
double conductance = cnt*1.0/min(vol1,2*m-vol1);
printf("conductance = %.5f\n",conductance);
return conductance;
}

double LocalFlow(double alpha){

init();
Bs.clear();
double ans = 0;
printf("alpha=%.5f,vol_a=%d,epsilon_sigma=%.5f,sigma=%.5f\n",alpha,vol_a,epsilon_sigma,sigma);
int max_number_of_iterations = (int)(5.0/alpha*log(3*vol_a/sigma)/log(2));
printf("max_number_of_iterations=%d\n",max_number_of_iterations);
int number_of_iterations = 0;

set<int>::iterator iter;

//构造local graph G''=G'<Bs>
//edges(s,v) for all v ∈ A
for(iter = A.begin();iter != A.end();iter++){
int v = *iter;
add_edge(S,v,deg[v]);
}

set<int> old_to_T_set;
set<int> old_A_union_Bs;

while(++number_of_iterations < max_number_of_iterations){
set<int> A_union_Bs;
set_union(A.begin(),A.end(),Bs.begin(),Bs.end(),inserter(A_union_Bs, A_union_Bs.begin()));

set<int> to_T_set;
set<int> N_A_union_Bs;

for(iter = A_union_Bs.begin();iter != A_union_Bs.end();iter++){
int u = *iter;
for(int i = 0;i < G[u].size();i++){
int v = G[u][i];
if(A_union_Bs.find(v) != A_union_Bs.end())  continue;
N_A_union_Bs.insert(v);
}
}
set_union(N_A_union_Bs.begin(),N_A_union_Bs.end(),Bs.begin(),Bs.end(),inserter(to_T_set, to_T_set.begin()));
//edges(v,t) for all v ∈ Bs union N(A union Bs)
for(iter = to_T_set.begin();iter != to_T_set.end();iter++){
int v = *iter;
if(old_to_T_set.find(v) != old_to_T_set.end())  continue;
add_edge(v,T,epsilon_sigma*deg[v]);
}
//edges(u,v) for all (v ∈ A union Bs) and (v ∈ V)
for(iter = A_union_Bs.begin();iter != A_union_Bs.end();iter++){
int u = *iter;
if(old_A_union_Bs.find(u) != old_A_union_Bs.end())  continue;
for(int i = 0;i < G[u].size();i++){
int v = G[u][i];
if(old_A_union_Bs.find(v) != old_A_union_Bs.end())  continue;
add_edge(u,v,1/alpha);
add_edge(v,u,1/alpha);
}
}

/*
因为每次增广只需要在局部的图上进行，
因此我们维护两个集合old_to_T_set和old_A_union_Bs，
这样每次寻找阻塞流之前通过对比，
我们就可以把Bs集合更新后需要加入局部图的新边在图上建边即可。
*/
old_to_T_set = to_T_set;
old_A_union_Bs = A_union_Bs;

//bfs确定层次图
if(!bfs())  break;
//寻找阻塞流
memcpy(work,first,sizeof(first));
while(true){
double t = dfs(S,INF);
if(abs(t) < eps)    break;
ans += t;
}
//寻找A union Bs中指向超级汇点t中已经饱和的边，生成新的Bs集合
for(iter = N_A_union_Bs.begin();iter != N_A_union_Bs.end();iter++){
int u = *iter;
//printf("u = %d\n",u);
for(int i = first[u];i != -1;i = edge[i].next){
int to = edge[i].v;
double val = edge[i].w;
if(to != T)  continue;
//printf("u = %d,val = %.3f\n",u,val);
if(val < eps){
Bs.insert(u);
break;
}
}
}
//puts("----------------------");
}
/*printf("number_of_iterations=%d\n",number_of_iterations);
for(int i = 0;i < e;i += 2){
printf("edge:from %d to %d,val = %.4f\n",edge[i].u,edge[i].v,edge[i].w);
}*/
printf("LocalFlow(%.5f)=%.5f\n",alpha,ans);
return ans;
}

int main()
{
freopen("wing_nodal.txt","r",stdin);
//freopen("output.txt","w",stdout);

scanf("%d%d%d%lf",&n,&m,&k,&epsilon_sigma);
sigma = 3*epsilon_sigma/(3*epsilon_sigma+1);
S = 0,T = n+1;
vol_a = 0;

for(int i = 0;i < m;i++){
int a,b;
scanf("%d%d",&a,&b);
G[a].push_back(b);
G[b].push_back(a);
deg[a]++;deg[b]++;
}
for(int i = 0;i < k;i++){
int a;
scanf("%d",&a);
A.insert(a);
vol_a += deg[a];
}
set<int>::iterator iter;
int cut_a = 0;//统计E(A,V-A),即跨越A与V-A的边的数量
for(iter = A.begin();iter != A.end();iter++){
int u = *iter;
for(int i = 0;i < G[u].size();i++){
int v = G[u][i];
if(A.find(v) != A.end())    continue;
cut_a++;
}
}
double conductance = cut_a*1.0/min(vol_a,2*m-vol_a);
printf("init:cut_a = %d,vol_a = %d,vol_a/vol_v=%.5f\n",cut_a,vol_a,vol_a*1.0/(2*m));
printf("before algorithm,conductance = %.5f\n",conductance);

clock_t time_begin,time_end;
time_begin = clock();
double alpha_max = 1.0,alpha_min = 0.0;
while(alpha_max - alpha_min > 1e-5){
double alpha_mid = (alpha_max+alpha_min)*0.5;
double max_flow = LocalFlow(alpha_mid);
conductance = min_cut();
if(abs(conductance) < eps){
alpha_min = alpha_mid;
}else{
alpha_max = alpha_mid;
}
}
puts("------------------result-----------------");
printf("best alpha=%.5f\n",alpha_max);
LocalFlow(alpha_max);
min_cut();
time_end = clock();
double delta = (double)(time_end - time_begin) / CLOCKS_PER_SEC;
printf("\n,delta_time = %.5f s\n",delta);
return 0;
}

View Code