最大流最小割的扩展应用

Leave a comment

  • 提出这个部分主要是针对最大流最小割的几个典型应用,而这几点在算导和黑书上没有详细描述。
  • 这里沿用《实用算法》上的三个典型应用:图的割切图的连通度图的边连通度。顺便把多源多汇的可行流容量上下界的最大流,最小路径覆盖[多源多汇的改造问题],以及最小费用最大流问题一起描述并实现了。[最小费用问题也可以用线型规划来解,就像差分也可以不用图而用线型规划来解,都是相通的。]
  • 这块有点太理论了,不过还是都实现一遍吧,最大流的难点在建模,编程都还在其次。但编程是基本功,还在其次的问题都不能秒杀,谈什么建模呢?

1.最小割(S,T)的解法:当流稳定时,换言之,从S出发无法找到一个前向流,使得当前流f再增加。定义残存前向图Df,边集Ef={(u,v)属于E,c(u,v)-f(u,v)>0},即正向流集,其点集V即为所求!

2.图连通度K(G),定义K(G)为:K(G)=顶点数-1,G为双向完全图;K(G)=Min{P(a,b)},a,b不相邻,G为非双向完全图。

  • 构造流图D:将原图中的u,v分裂为两个等效顶点,以单位容量边相连,将原来的边映射为+无穷边
  • 求出不相邻点a到b的最大流,以及对应的割顶集S。
k连通图的基本性质这里略过。
下面实现一个基本连通图的k值以及割顶求法,测试数据如下:[实用P259页基本连通图]
7 12[顶点数,边数]
1 6
1 7
1 2
2 7
2 3
3 7
4 7
4 3
5 7
5 4
5 6
6 7
算法框架:
//转化图框架:有向图,无向图处理
//验证是否为完全图
//若非完全图,则枚举源点和目标点,结算最大流;
//若最大流比当前的最大值大,在用dfs求出最小割集,逻辑很清楚!
下面是代码实现:
#include <stdlib.h>
#include <assert.h>
#include <iostream>
#include <fstream>
using namespace std;
const int MAXN=30, MAXE=200;
//for flow graph
typedef struct Edge{
 int i, j, next, op, f, c;
};
Edge es[MAXE];
int eid=0, n, en, first[MAXN], pre[MAXN], ptr[MAXN], first1[MAXN], queue[MAXN], ansp[MAXN];
bool visited[MAXN];
//for orignal graph
bool map[MAXN][MAXN];
//for S set
bool ss[MAXN];
//declaration
void createGraph();
void addEdge(int i, int j, int c);
void Dicnic(int sv, int st);
void bfslevel(int sv, int st);
void dfs_maxflow(int sv, int st);
void dfsforS(int v);
//implementation
void createGraph(){
 ifstream in("input.txt");
 if(!in) exit(EXIT_FAILURE);
 cin.rdbuf(in.rdbuf());
 //read in for n, e
 cin>>n>>en;
 int x, y;
 //clr for all previous list
 eid=0, memset((void*)first, 0xff, sizeof(first)), memset((void*)map, 0, sizeof(map));
 for(int i=0; i<n; i++) addEdge(i, i+n, 1);//i -> i+n in weight of 1
 for(int i=0; i<en; i++){
  cin>>x>>y, map[x][y]=map[y][x]=true;
  addEdge(x+n, y, INT_MAX), addEdge(y+n, x, INT_MAX);
 }
 in.close();
}
void addEdge(int i, int j, int c){
 es[eid].i=i, es[eid].j=j, es[eid].f=0, es[eid].c=c;
 es[eid].next=first[i], first[i]=eid, es[eid].op=eid+1; eid++;
 es[eid].i=j, es[eid].j=i, es[eid].f=0, es[eid].c=-1;
 es[eid].next=first[j], first[j]=eid, es[eid].op=eid-1; eid++;//here:bugs,开始写成了后减,精力有点不集中!
}
//for Dicnic implementation
void Dicnic(int sv, int st){
 while(true){
  bfslevel(sv, st);
  if(pre[st]==INT_MAX) return;
  dfs_maxflow(sv, st);
 }
}
void bfslevel(int sv, int st){
 for(int i=0; i<2*n; i++) pre[i]=INT_MAX;
 int top=0, size=1; queue[0]=sv, pre[sv]=0;//here:bugs注意源点赋值
 while(top<size){
  int cn=queue[top++];//here!bugs队列的头处理
  if(cn==st) return;
  for(int e=first[cn]; e!=-1; e=es[e].next){
   if(es[e].c>es[e].f||(es[e].c==-1&&es[e].f>0)){
    if(pre[es[e].j]>pre[cn]+1){
     queue[size++]=es[e].j, pre[es[e].j]=pre[cn]+1;
    }
   }
  }
 }
}
void dfs_maxflow(int sv, int st){
 memcpy((void*)first1, (void*)first, sizeof(first));
 memset((void*)visited, 0, sizeof(visited));
 int top=0; queue[0]=sv;
 while(top>-1){
  if(queue[top]==st){
   int fmin=INT_MAX, value;
   for(int i=top; i>0; i–){
    if(es[ptr[queue[i]]].c>es[ptr[queue[i]]].f){
     if((value=es[ptr[queue[i]]].c-es[ptr[queue[i]]].f)<fmin) fmin=value;
    }
    if(es[ptr[queue[i]]].c==-1&&es[ptr[queue[i]]].f>0){
     if((value=es[ptr[queue[i]]].f)<fmin) fmin=value;
    }
   }
   //update the maxflow
   for(int i=top; i>0; i–){
    if(es[ptr[queue[i]]].c>es[ptr[queue[i]]].f)
     es[ptr[queue[i]]].f+=fmin, es[es[ptr[queue[i]]].op].f+=fmin;
    if(es[ptr[queue[i]]].c==-1&&es[ptr[queue[i]]].f>0)
     es[ptr[queue[i]]].f-=fmin, es[es[ptr[queue[i]]].op].f-=fmin;
   }
   return;
  }
  else{//dfs for maxflow
   int tp=queue[top], tmp=first1[tp];
   for(;tmp!=-1;tmp=es[tmp].next){
    if(!visited[es[tmp].j]){
     if((es[tmp].c>es[tmp].f||(es[tmp].c==-1&&es[tmp].f>0))&&pre[es[tmp].j]==pre[tp]+1){
      first1[tp]=es[tmp].next;
      queue[++top]=es[tmp].j, ptr[es[tmp].j]=tmp;
      break;
     }//if
    }//if
   }//for
   if(tmp==-1){visited[tp]=true; top–;}
  }
 }
}
void dfsforS(int v){
 ss[v]=true;
 for(int e=first[v]; e!=-1; e=es[e].next){//here:bugs
  if(!ss[es[e].j]&&
es[e].c>es[e].f) dfsforS(es[e].j);
 }
}
int _tmain(int argc, _TCHAR* argv[])
{
 createGraph();
 bool bf=true;
 for(int i=1; i<n; i++){
  for(int j=i+1; j<n; j++){
   assert(i!=j);
   if(!map[i][j]){bf=false; break;}
  }
  if(!bf) break;
 }
 if(bf){//bi-connected graph
  cout<<n-1<<endl;
  for(int i=0; i<n; i++) cout<<i<<" ";
  cout<<endl;
  goto jmp;
 }
 else{//for k-connected graph
  int ans=0;
  for(int vs=n; vs<2*n; vs++){
   for(int vt=0; vt<n; vt++){
    if(vs-n!=vt&&!map[vs-n][vt]){//maxflow for between vs and vt
     for(int i=0; i<eid; i++) es[i].f=0;
     Dicnic(vs, vt);
     int k=0;
     for(int i=first[vs]; i!=-1; i=es[i].next) k+=es[i].f;
     if(k>ans){
      ans=k;
      memset((void*)ss, 0, sizeof(ss));
      dfsforS(vs);
      int t=0;
      for(int i=0; i<n; i++)
       if(ss[i]&&!ss[i+n]) ansp[t++]=i;
     }//if
    }
   }
  }
  cout<<ans<<endl;
  for(int i=0; i<ans; i++) cout<<ansp[i]<<" ";
  cout<<endl;
 }
jmp:system("pause");
 return 0;
}
 
3.图的边连通度,类似于图顶点连通度;
K'(G),对于a,b顶点,从a到b的无公共编的轨,称为弱独立轨.其定义与K(G)是一致的。
通用构造算法如下:
  • 设在无向图中,将原有边e(u,v)映射为e’=(u, v), e”=(v,u),其容量均为1, 而有向图,则直接设置存在边容量为1。[这个和最大流模型是保持一致的]
  • 排除特殊的完全图情况,求出最大流,进而dfs求出在原图中的桥集
下面实现一个基本边连通图的k值以及桥集求法,测试数据如下:[实用P261页基本连通图]
5 6[顶点数,边数]
1 2
1 5
2 5
2 3
2 4
3 4
下面是完整代码实现:
#include <stdlib.h>
#include <assert.h>
#include <iostream>
#include <fstream>
using namespace std;
const int MAXN=40, MAXE=200;
typedef struct Edge{
 int i, j, next, op, c, f;
};
typedef struct Bridge{
 int i, j;
};
Edge es[MAXE];
Bridge bs[MAXE];
int eid=0, nv, ev, first[MAXN], first1[MAXN], pre[MAXN], ptr[MAXN], queue[MAXN];
bool visited[MAXN], ss[MAXN], map[MAXN][MAXN];
//declaration
void createGraph();
void addEdge(int i, int j, int c);
void Dicnic(int sv, int st);
void bfslevel(int sv, int st);
void dfs_maxflow(int sv, int st);
void dfsforS(int u);
//implementation
void createGraph(){
 ifstream in("input.txt");
 if(!in) exit(EXIT_FAILURE);
 cin.rdbuf(in.rdbuf());
 cin>>nv>>ev;
 int x, y;
 eid=0, memset((void*)first, 0xff, sizeof(first)), memset((void*)map, 0, sizeof(map));
 for(int i=0; i<ev; i++){
  cin>>x>>y, x–, y–;
  map[x][y]=map[y][x]=true, addEdge(x, y, 1), addEdge(y, x, 1);
 }
}
void addEdge(int i, int j, int c){
 es[eid].i=i, es[eid].j=j, es[eid].f=0, es[eid].c=c;
 es[eid].next=first[i], first[i]=eid, es[eid].op=eid+1; eid++;
 es[eid].i=j, es[eid].j=i, es[eid].f=0, es[eid].c=-1;
 es[eid].next=first[j], first[j]=eid, es[eid].op=eid-1; eid++;
}
void Dicnic(int sv, int st){
 while(true){
  bfslevel(sv, st);
  if(pre[st]==INT_MAX) return;
  dfs_maxflow(sv, st);
 }
}
void bfslevel(int sv, int st){
 for(int i=0; i<nv; i++) pre[i]=INT_MAX;
 int top=0, size=1; queue[0]=sv, pre[sv]=0;
 while(top<size){
  int cn=queue[top++];
  if(cn==st) return;
  for(int e=first[cn]; e!=-1; e=es[e].next){
   if((es[e].c>es[e].f||(es[e].c==-1&&es[e].f>0))&&pre[es[e].j]>pre[cn]+1)//here: bugs
    queue[size++]=es[e].j, pre[es[e].j]=pre[cn]+1;
  }
 }
}
void dfs_maxflow(int sv, int st){
 memcpy((void*)first1, (void*)first, sizeof(first));
 memset((void*)visited, 0, sizeof(visited));
 int top=0; queue[0]=sv;
 while(top>-1){
  if(queue[top]==st){
   int fmin=INT_MAX, value;
   for(int i=top; i>0; i–){
    if(es[ptr[queue[i]]].c>es[ptr[queue[i]]].f){
     if((value=es[ptr[queue[i]]].c-es[ptr[queue[i]]].f)<fmin) fmin=value;
    }
    if(es[ptr[queue[i]]].c==-1&&es[ptr[queue[i]]].f>0){
     if((value=es[ptr[queue[i]]].f)<fmin) fmin=value;
    }
   }
   for(int i=top; i>0; i–){
    if(es[ptr[queue[i]]].c>es[ptr[queue[i]]].f){
     es[ptr[queue[i]]].f+=fmin, es[es[ptr[queue[i]]].op].f+=fmin;
    }
    if(es[ptr[queue[i]]].c==-1&&es[ptr[queue[i]]].f>0){
     es[ptr[queue[i]]].f-=fmin, es[es[ptr[queue[i]]].op].f-=fmin;
    }
   }
   return;
  }
  else{//for first1 comparsion
   int tp=queue[top], tmp=first1[tp];
   for(;tmp!=-1;tmp=es[tmp].next){
    if(!visited[es[tmp].j])
     if((es[tmp].c>es[tmp].f||(es[tmp].c==-1&&es[tmp].f>0))&&pre[es[tmp].j]==pre[tp]+1){
      queue[++top]=es[tmp].j, ptr[es[tmp].j]=tmp;
      break;
     }
   }//for
   if(tmp==-1){visited[tp]=true; top–;}
  }
 }
}
void dfsforS(int u){
 ss[u]=true;
 for(int e=first[u]; e!=-1; e=es[e].next){
  if(!ss[es[e].j]&&es[e].c>es[e].f) dfsforS(es[e].j);
 }
}
int _tmain(int argc, _TCHAR* argv[])
{
 createGraph();
 bool bf=true;
 for(int i=0; i<nv; i++){
  for(int j=i+1; j<nv; j++){
   if(!map[i][j]){bf=false; break;}
  }
 }
 if(bf){//full-connencted
  cout<<nv-1<<endl;
  for(int i=1; i<nv; i++) cout<<"0 "<<i<<endl;
  goto jmp;
 }
 else{//for the intialization of edages
  int ans=0;
  for(int s=0; s<nv; s++){
   for(int t=0; t<nv; t++){
    if(s!=t&&!map[s][t]){//find non-adjancent two vertexs
     for(int i=0; i<eid; i++) es[i].f=0;
     Dicnic(s, t);
     int k=0;
     for(int e=first[s]; e!=-1; e=es[e].next) k+=es[e].f;
     if(k>ans){//find the according bridge
      ans=k;
      memset((void*)ss, 0, sizeof(ss));
      dfsforS(s); int t=0;
      for(int i=0; i<nv; i++){
       for(int j=0; j<nv; j++){
        if(ss[i]&&!ss[j]&&map[i][j]) bs[t].i=i, bs[t++].j=j;
       }
      }
     }//if
    }//if
   }//for vt
  }//for vs
  cout<<ans<<endl;
  for(int i=0; i<ans; i++) cout<<bs[i].i<<" "<<bs[i].j<<endl;
 }
jmp:system("pause");
 return 0;
}
 
4.多源多汇网络的可行流计算
这类型题目在算导中有所提及,思路是对于多个源点加入一个pre-source,对于多个目标,加入一个post-dest,这样模型就统一到了熟悉的模型之上。
引用《实用》上的例题来实现一个具体例子。
 
商品供应系统问题:
如图商品供应系统,边的方向表示商品流向,边上的数字表示商品从某地到某地的最大通过能力。商品的货源来自于x1,x2和x3,其中x1的供应量为5,x2为10,x3为5。这些商品在商场y1,y2和y3出售。y1的需求为5,y2为10,y3为5。是否能同时满足?

输入数据:
10 9
1 4 3
1 5 3
2 1 7
2 5 4
2 3 9
3 5 3
3 7 7
4 8 8
4 6 3
5 4 2
5 6 8
6 8 1
6 9 2
6 10 4
7 5 1
7 6 6
7 10 8
8 9 2
10 9 7
3 1 5 2 10 3 5
3 8 5 9 10 10 5
 
5.容量下界的网络流计算
这个算法看起太奇特了,难度完全在图的转化上,首先是将上下界分离为必要弧和增量弧,然后加入x,y的源点和汇点,再进行改造;算法首先从x,y出发找到一个对改造网络的必要弧饱和处理方案,然后再对残留网络进行扩展!说的自己都有点晕了,画了个转化图发现今天思维迟钝的不可理解,算了,暂时放一下。遇到思维短路期了!
 
6.最小费用最大流 网上的两种思路
这个问题再短路就简直不可原谅了。
  • Ford-Fulkerson增广最短路算法[Successive Shortest Path]:在最大流系统中加入单位流量的费用因子,然后求在最小费用基础上的最大流,输出为总的输送费用B(F)=Sum(Bij*Fij).这个思路挺典型的,感觉有点不像是最大流,倒像是最短路+增量算法的问题。
引用别人的论述:把各条弧上单位流量的费用看成某种长度,用求解最短路问题的方法确定一条自V1至Vn的最短路;在将这条最短路作为可扩充路,用求解最大流问题的方法将其上的流量增至最大可能值;而这条最短路上的流量增加后,其上各条弧的单位流量的费用要重新确定,如此多次迭代,最终得到最小费用最大流.
算法描述:
  • 对于G求出其从S到T的最短路[Bellman-Ford],初始时F=0。
  • 在F(k-1)基础上进行流增量,实际上是求出可达最短路上的瓶颈流量。尤其最短路时已经求出了累积的最小费用,直接用l[t]*瓶颈流量即可。

代码框架:[实用上的框架]

Successive Shortest Path

1    Establish a feasible flow x when F=0
2    while ( Gx contains any shortest path in F(k-1) ) do
3        find the shortest path from s to t
4       
5        add the cost from the shortest path from s to t in l[t] 
6        find the bottleneck of flow(k-1) and the sum of cost is l[t]*|bottleneck|

  • 消圈算法[Cycle Canceling]:网络G的最大流是G的一个最小费用流的充要条件是,最大流的残图中不含负费用圈

伪代码框架:[引自 Zealint ‘s Minimum Cost Flow, Part 2: Algorithms ]

 Cycle-Canceling

1    Establish a feasible flow x in the network
2    while ( Gx contains a negative cycle ) do
3        identify a negative cycle W
4       
5        augment  units of flow along the cycle W
6        update Gx

一个练手的基本题目:在如下图中,求出其最小费用最大流。

测试数据:
5 7
0 1 10 4
0 2 8 1
2 1 5 2
1 3 2 6
2 3 10 3
1 4 7 1
3 4 4 2
 
  • Cycle cancel算法实现:实用的模型确实很强,那种模型的处理,使得正流和负流统一起来了!然后沿正向流的增加[残图中g的增加]和负流的处理都一致料。呵呵,最后的费用流只需要把正向流的集合跑一次即可!来分析下算法复杂度:Bellman-Ford的负圈测试O(VE),每条边的容量上限M,费用上限为C,因此最大流的费用为mCM,统一之复杂度为O(E^2VCM).

代码实现:[可用为模板]

#include <iostream>
#include <fstream>

using namespace std;

//Cycle-canceling Algorithm
//Code Skeleton
//1    Establish a feasible flow x in the network
//2    while ( Gx contains a negative cycle ) do
//3        identify a negative cycle W
//4       
//5        augment  units of flow along the cycle W
//6        update Gx

//declaration
typedef struct Edge{
 int i, j, f, c, w, op, next;//using uniform Edge presetation
};
const int MAXN=20, MAXE=400;
Edge es[MAXE];
int eid=0, first[MAXN], hv[MAXN], ev[MAXN], l[MAXN], w[MAXN], pre[MAXN], ptr[MAXN];//for bellman-ford loop
int nv, ne;

//declaration
void createGraph();
void addEdge(int start, int end, int cap, int cost);

//maxflow algorithm
void highpreflow();
void initialize();
void discharge(int u);
bool relabel(int u);
bool push(int eid);

//implementation
void createGraph(){
 ifstream in("input.txt");
 if(!in) exit(EXIT_FAILURE);
 cin.rdbuf(in.rdbuf());
 //read in data and establish the graph
 cin>>nv>>ne; eid=0, memset((void*)first, 0xff, sizeof(first));
 int start, end, cap, cost;
 for(int i=0; i<ne; i++)
  cin>>start>>end>>cap>>cost, addEdge(start, end, cap, cost);
}

void addEdge(int start, int end, int cap, int cost){
 if(!cap) return;
 es[eid].i=start, es[eid].j=end, es[eid].f=0, es[eid].c=cap, es[eid].w=cost;
 es[eid].next=first[start], first[start]=eid, es[eid].op=eid+1; eid++;
 es[eid].i=end, es[eid].j=start, es[eid].f=0, es[eid].c=-1, es[eid].w=-cost;
 es[eid].next=first[end], first[end]=eid, es[eid].op=eid-1; eid++;
}

void highpreflow(){
 initialize();
 int pos=0, u=l[pos];
 while(u!=-1){
  int oldh=hv[u];
  discharge(u);
  if(oldh!=hv[u]){
   if(pos!=0) l[pos]^=l[0], l[0]^=l[pos], l[pos]^=l[0];
   pos=0;
  }
  u=l[++pos];
 }
}

void initialize(){
 memset((void*)hv, 0, sizeof(hv)), memset((void*)ev, 0, sizeof(ev));
 memset((void*)l, 0xff, sizeof(l));
 for(int i=0; i<nv-2; i++) l[i]=i+1;
 hv[0]=nv-1;
 for(int i=first[0]; i!=-1; i=es[i].next){
  if(es[i].c>0){
   es[i].f+=es[i].c, es[es[i].op].f+=es[i].c;
   ev[es[i].j]+=es[i].c, ev[0]-=es[i].c;
  }
 }
}

bool relabel(int u){
 if(ev[u]<=0) return false;
 int hmin=INT_MAX; bool minV=true;
 for(int e=first[u]; e!=-1; e=es[e].next){
  if(es[e].c>es[e].f||(es[e].c==-1&&es[e].f>0)){
   if(hv[es[e].j]<hv[u]) minV=false;
   if(hv[es[e].j]<hmin) hmin=hv[es[e].j];
  }
 }
 if(minV&&hmin<=nv){hv[u]=hmin+1; return true;}
 return false;
}

bool push(int eid){//es[eid].i->es[eid].j
 if(ev[es[eid].i]<=0||hv[es[eid].i]!=hv[es[eid].j]+1) return false;
 int dt;
 if(es[eid].c>es[eid].f){
  dt=min(ev[es[eid].i], es[eid].c-es[eid].f);
  es[eid].f+=dt, es[es[eid].op].f+=dt;
  ev[es[eid].i]-=dt, ev[es[eid].j]+=dt;
  return true;
 }
 if(es[eid].c==-1&&es[eid].f>0){
  dt=min(ev[es[eid].i], es[eid].f);
  es[eid].f-=dt, es[es[eid].op].f-=dt;
  ev[es[eid].i]-=dt, ev[es[eid].j]+=dt;
  return true;
 }
 return false;
}

void discharge(int u){
 int cur=first[u];
 while(ev[u]>0){
  if(cur==-1) relabel(u), cur=first[u];
  else if((es[cur].c>es[cur].f||(es[cur].c==-1&&es[cur].f>0))&&hv[u]==hv[es[cur].j]+1)
   push(cur);
  else cur=es[cur].next;
 }
}

int _tmain(int argc, _TCHAR* argv[])
{
 createGraph();
 highpreflow();

 //cycle cancel
 bool cycle=false;
 do{//cycle cancel logic
  w[0]=0, memset((void*)pre, 0xff, sizeof(pre));
  for(int i=1; i<nv; i++) w[i]=INT_MAX; cycle=false;
  int start=-1, fmin=INT_MAX;
  for(int i=1; i<nv; i++){
   for(int e=0; e<eid; e++){
    if(es[e].c>es[e].f||(es[e].c==-1&&es[e].f>0))//in the residental graph
     if(w[es[e].i]!=INT_MAX&&w[es[e].i]+es[e].w<w[es[e].j]) w[es[e].j]=w[es[e].i]+es[e].w, pre[es[e].j]=es[e].i, ptr[es[e].j]=e;
   }
  }

  for(int i=1; i<nv; i++){
   for(int e=0; e<eid; e++){
    if(es[e].c>es[e].f||(es[e].c==-1&&es[e].f>0))//in the residental graph
     if(w[es[e].i]!=INT_MAX&&w[es[e].i]+es[e].w<w[es[e].j]){ start=es[e].j, cycle=true; break;}
   }//Bellman-Ford一行不差的!
   if(cycle) break;//existing one minus cycle
  }
  if(!cycle) break;//no-minus cycle in the existing system
  //here: if self-loop? ignore this case
  int nxt=pre[start], value;
  do{
   if(es[ptr[nxt]].c>es[ptr[nxt]].f){//here: w’s value
    if((value=es[ptr[nxt]].c-es[ptr[nxt]].f)<fmin) fmin=value;
   }
   if(es[ptr[nxt]].c==-1&&es[ptr[nxt]].f>0){
    if((value=es[ptr[nxt]].f)<fmin) fmin=value;
   }//please check
  }while(nxt!=pre[start]);
  //update the minus cycle
  nxt=pre[start];
  do{//bugs : in the residental graph, decreas the flow on this direction, increase the opposite direction
   if(es[ptr[nxt]].c>es[ptr[nxt]].f)
    es[ptr[nxt]].f+=fmin, es[es[ptr[nxt]].op].f+=fmin;
   else if(es[ptr[nxt]].c==-1&&es[ptr[nxt]].f>0)
    es[ptr[nxt]].f-=fmin, es[es[ptr[nxt]].op].f-=fmin;
   /*assert(es[ptr[nxt]].f>=0&&es[es[ptr[nxt]].op].f>=0);*/
   nxt=pre[nxt];
  }while(nxt!=pre[start]);//这个地方有点搞,呵呵!闪了一下,证明:存在圈那么从start开始一定有一个出,一个入,OK,用pre验证即可。
 }while(cycle);

//MIT上说要再走一次dfs来看是否确实无圈!
 //find the minmal cost flow along the forward-flow
 int forwardflow=0;
 for(int i=0; i<eid; i++){
  if(es[i].c>0) forwardflow+=es[i].f*es[i].w;
 }
 cout<<forwardflow<<endl;
 return 0;
}

  • Successive Shortest Path算法实现:每次找出更新流后的最小费用路径,找出瓶颈流以及用w[st]表示其最短路径的费用长度,累计即是最小费用。算法复杂度:Bellman-Ford的负圈测试O(VE),每条边的容量上限M,费用上限为C,因此每次找最小路径O(S(m,n,C)),统一之复杂度为O(MS(m,n,C)).

代码实现:[可用为模板]

#include <iostream>
#include <fstream>

using namespace std;

/* template from Top coder, implementation refer to practical alogrithm.
Successive Shortest Path
1    Transform network G by adding source and sink
2    Initial flow x is zero
3    while ( Gx contains a path from s to t ) do
4        Find any shortest path P from s to t
5        Augment current flow x along P
6        update Gx
*/

const int MAXN=20, MAXE=400;
typedef struct Edge{
 int i, j, next, op, f, c, w;
}Edge;
Edge es[MAXE];
int eid=0, first[MAXN], w[MAXN],/*accumulated value*/ pre[MAXN], ptr[MAXN];
int nv, ne;

//declaration
void addEdge(int i, int j, int c, int w);
int expandPath(int sv, int st);
void createGraph();

//implementation
void addEdge(int i, int j, int c, int w){
 if(!c) return;
 es[eid].i=i, es[eid].j=j, es[eid].c=c, es[eid].w=w, es[eid].f=0;
 es[eid].next=first[i], first[i]=eid, es[eid].op=eid+1; eid++;
 es[eid].i=j, es[eid].j=i, es[eid].c=-1, es[eid].w=-w, es[eid].f=0;
 es[eid].next=first[j], first[j]=eid, es[eid].op=eid-1; eid++;
}

void createGraph(){
 ifstream in("input.txt");
 if(!in) exit(EXIT_FAILURE);
 cin.rdbuf(in.rdbuf());
 cin>>nv>>ne;eid=0, memset((void*)first, 0xff, sizeof(first));
 int start, end, cap, costw;
 for(int i=0; i<ne; i++)
  cin>>start>>end>>cap>>costw, addEdge(start, end, cap, costw);
}

int expandPath(int sv, int st){
 memset((void*)pre, 0xff, sizeof(int)*nv);
 for(int i=0; i<nv; i++) w[i]=INT_MAX;
 pre[sv]=sv,w[sv]=0; bool change=false;
 do{
  change=false;
  for(int i=1; i<nv; i++){
   for(int e=0; e<eid; e++){//flow may be changed during process
    if(pre[es[e].i]!=-1){
     if(es[e].c>es[e].f)
      if(w[es[e].i]+es[e].w<w[es[e].j]) w[es[e].j]=w[es[e].i]+es[e].w, pre[es[e].j]=es[e].i, ptr[es[e].j]=e, change=true;
     if(es[e].c==-1&&es[e].f>0)
      if(w[es[e].i]+es[e].w<w[es[e].j]) w[es[e].j]=w[es[e].i]+es[e].w, pre[es[e].j]=es[e].i, ptr[es[e].j]=e, change=true;
    }
   }
  }
 }while(change);//if change then we continue change
 //if we can’t from sv to st;
 if(pre[st]==-1) return 0;
 int i=st, a=INT_MAX, value;
 for(;i!=sv;i=pre[i]){
  if(es[ptr[i]].c>es[ptr[i]].f){//here:bugs i->ptr[i]
   if((value=es[ptr[i]].c-es[ptr[i]].f)<a) a=value;
  }
  if(es[ptr[i]].c==-1&&es[ptr[i]].f>0){
   if((value=es[ptr[i]].f)<a) a=value;
  }

 }
 return a;
}

int _tmain(int argc, _TCHAR* argv[])
{
 createGraph();
 int a=expandPath(0, nv-1), aflow=0, sv=0, st=nv-1;
 while(a){
  for(int i=st;i!=sv;i=pre[i]){//change the flow
   if(es[ptr[i]].c>es[ptr[i]].f) es[ptr[i]].f+=a, es[es[ptr[i]].op].f+=a;
   if(es[ptr[i]].c==-1&&es[ptr[i]].f>0) es[ptr[i]].f-=a, es[es[ptr[i]].op].f-=a;
  }
  aflow+=w[st]*a;
  a=expandPath(0, nv-1);
 }
 cout<<aflow<<endl;
 return 0;
}

Advertisements

网络流例题

Leave a comment

题目描述的很清楚,转化一下就是,0为源点,n+d+1为终点,1->N为牛点,N+1->N+D为食物点,从0到1->N容量为k,从N+1->N+D到n+d+1容量为Di,加上1容量,OK,直接用Dinic算法切掉。
犯了一个错,if(es[ptr[queue[i]]].c==-1&&es[ptr[queue[i]]].f>0){,等号写成赋值了,看了半天,呵呵!以后要小心罗,这种错误不值得!
 
代码实现如下:
#include <iostream>
#include <fstream>
using namespace std;
//Problems – 2002 Winter Open, Green (Senior Division)
//Problem 2 – New years party
const int MAXN=310, MAXE=21000;
typedef struct Edge{
    int i, j, next, op, f, c;
};
Edge es[MAXE];
int eid=0, n, d, k, nt, first[MAXN], pre[MAXN], ptr[MAXN], queue[MAXN], first1[MAXN];
bool visited[MAXN];
//declaration
void addEdge(int i, int j, int c);
int Dicnic();
void bfslevel();
void dfs_maxflow();
//implemenation, eid=0, reset
void addEdge(int i, int j, int c){
    es[eid].i=i, es[eid].j=j, es[eid].f=0, es[eid].c=c;
    es[eid].next=first[i], first[i]=eid, es[eid].op=eid+1; eid++;
    es[eid].i=j, es[eid].j=i, es[eid].f=0, es[eid].c=-1;//reverse non-sense
    es[eid].next=first[j], first[j]=eid, es[eid].op=eid-1; eid++;
}
void bfslevel(){
    for(int i=0; i<=nt; i++) pre[i]=INT_MAX;
    int top=0, size=1; queue[top]=0, pre[0]=0;
    while(top<size){
        int cn=queue[top++];
        if(cn==nt) return;
        for(int e=first[cn]; e!=-1; e=es[e].next){
            if(pre[es[e].j]>pre[cn]+1){//in the residental graph
                if(es[e].c>es[e].f||(es[e].c==-1&&es[e].f>0)){
                    queue[size++]=es[e].j, pre[es[e].j]=pre[cn]+1;
                }
            }
        }
    }
}
void dfs_maxflow(){
    memcpy((void*)first1, (void*)first, sizeof(first));
    memset((void*)visited, 0, sizeof(visited));
    int top=0; queue[0]=0;
    while(top>-1){
        if(queue[top]==nt){
            int fmin=INT_MAX, value;
            for(int i=top; i>0; i–){
                if(es[ptr[queue[i]]].c>es[ptr[queue[i]]].f){
                    if((value=es[ptr[queue[i]]].c-es[ptr[queue[i]]].f)<fmin) fmin=value;
                }
                if(es[ptr[queue[i]]].c==-1&&es[ptr[queue[i]]].f>0){
                    if((value=es[ptr[queue[i]]].f)<fmin) fmin=value;
                }
            }
            //update maxflow
            for(int i=top; i>0; i–){
                if(es[ptr[queue[i]]].c>es[ptr[queue[i]]].f){
                    es[ptr[queue[i]]].f+=fmin, es[es[ptr[queue[i]]].op].f+=fmin;
                }
                if(es[ptr[queue[i]]].c==-1&&es[ptr[queue[i]]].f>0){
                    es[ptr[queue[i]]].f-=fmin, es[es[ptr[queue[i]]].op].f-=fmin;
                }
            }
            return;//break;
        }
        else{
            int tp=queue[top], tmp=first1[tp];
            for(;tmp!=-1;tmp=es[tmp].next){
                if(!visited[es[tmp].j]){
                    if((es[tmp].c>es[tmp].f||(es[tmp].c==-1&&es[tmp].f>0))&&pre[es[tmp].j]==pre[tp]+1){
                        first1[tp]=es[tmp].next;
                        queue[++top]=es[tmp].j, ptr[es[tmp].j]=tmp;
                        break;
                    }
                }
            }//for
            if(tmp==-1){visited[tp]=true; top–;}
        }
    }
}
int Dicnic(){
    while(true){
        bfslevel();
        if(pre[nt]==INT_MAX) break;
        dfs_maxflow();
    }
    //get maxflow for source output
    int maxflow=0;
    for(int e=first[0]; e!=-1; e=es[e].next){
        maxflow+=es[e].f;
    }
    return maxflow;
}
#define fbase
int main()
{
#ifdef fbase
    ifstream in("input.txt");
    if(!in) return EXIT_FAILURE;
    cin.rdbuf(in.rdbuf());
#endif
    //implementation
    int dtoT, m, ntod;
    cin>>n>>k>>d; nt=n+d+1, memset((void*)first, 0xff, sizeof(first));
    for(int i=1; i<=n; i++) addEdge(0, i, k);//, cout<<0<<" "<<i<<" "<<k<<endl;//add 0->Ni for k
    for(int i=1; i<=d; i++){
        cin>>dtoT, addEdge(i+n, nt, dtoT);//, cout<<i+n<<" "<<nt<<" "<<dtoT<<endl;
    }
    for(int i=1; i<=n; i++){
        cin>>m;
        for(int j=1; j<=m; j++){
            cin>>ntod, addEdge(i, ntod+n, 1);//, cout<<i<<" "<<ntod+n<<" "<<1<<endl;
        }
    }
    cout<<Dicnic()<<endl;
#ifdef fbase
    in.close();
#endif
    return 0;
}
#undef fbase
 
2.航天飞机实验[算法导论P427,lrjP317]
最大流最小割解题:最大流确实是很强的工具,有点喜欢最高余流的算法,几行就搞定问题,而且思路直观!
自己写了个测试数据:
10 6 [仪器数,试验数]
1 2 3 4 5 6 [每个试验的净收入]
1 2 3 4 5 6 7 8 9 2 [每个仪器的维护费用]
4 1 2 5 7 [试验1的对应仪器]
3 4 7 8
4 9 7 3 5
3 1 9 6
5 1 4 8 9 3
3 3 6 9
证明:首先要证明最大流就是其收入最大值,先建立模型,将仪器作为源点的邻接点,而将试验作为目标点的邻接点,作一个二分,其权值分别为仪器支付费用,以及
试验净收入。如图

一个基本判断是,对于一个基本割集T,若Ei在其中,那么对应的Ii也应该在T中。[无穷边不可能在最小割中]
1,最小割边不可能通过+无穷边;
2,必定有一个最小割集合的容量,sum(Ii)+sum(Ej),若Ii是在T集中的,那么Ii肯定计算在最小割的值中了;而对应的Ej则没有计算在最小割中,所以Ii为运输的仪器,而Ej为没有参加的试验;转化一下,和净收入是一致的,因此最大流就是其净收入。
因此,直接用预流推进求出最大流,然后用Df求出T集即可。[理论写在最大流最小割的扩展应用部分中]。
 
代码实现:[in和原来的文件流名一样了,只有改下了,写代码的时候,还是要注意变量命名的,太随意了]
#include <iostream>
#include <fstream>
using namespace std;
const int MAXE=150, MAXN=40;
typedef struct Edge{
    int i, j, next, f, c, op;
};
Edge es[MAXE];
int eid=0, first[MAXN], l[MAXN], hv[MAXN], ev[MAXN];
int in, en, n;
bool inS[MAXN];
//declaration
void addEdge(int i, int j, int c);
void Hpreflow();
void initailflow();
void discharge(int u);
bool push(int eid);
bool relabel(int u);
void dfsforS(int v);
//implementation
void addEdge(int i, int j, int c){
    if(!c) return;
    es[eid].i=i, es[eid].j=j, es[eid].f=0, es[eid].c=c;
    es[eid].next=first[i], first[i]=eid, es[eid].op=eid+1; eid++;
    es[eid].i=j, es[eid].j=i, es[eid].f=0, es[eid].c=-1;
    es[eid].next=first[j], first[j]=eid, es[eid].op=eid-1; eid++;
}
void Hpreflow(){
    initailflow();
    int pos=0, u=l[pos];
    while(u!=-1){
        int oldh=hv[u];
        discharge(u);
        if(oldh<hv[u]){//relabel it
            if(pos!=0) l[pos]^=l[0], l[0]^=l[pos], l[pos]^=l[0];
            pos=0;
        }
        u=l[++pos];
    }
}
//initial for preflow
void initailflow(){
    memset((void*)hv, 0, sizeof(hv)), memset((void*)ev, 0, sizeof(ev)),
    memset((void*)l, 0xff, sizeof(l));
    for(int i=0; i<n-1; i++) l[i]=i+1;//l->1…n-1
    hv[0]=n;
    for(int e=first[0]; e!=-1; e=es[e].next){
        if(es[e].c>0){//for preflow
            es[e].f+=es[e].c, es[es[e].op].f+=es[e].c;
            ev[0]-=es[e].c, ev[es[e].j]+=es[e].c;
        }
    }
}
//discharge for node u;
void discharge(int u){
    int tmp=first[u];
    while(ev[u]>0){
        if(tmp==-1) relabel(u), tmp=first[u];
        else if((es[tmp].c>es[tmp].f||(es[tmp].c==-1&&es[tmp].f>0))&&hv[u]==hv[es[tmp].j]+1)
            push(tmp);//u is higher than es[tmp].j, just one unit
        else tmp=es[tmp].next;
    }
}
bool push(int eid){
    if(ev[es[eid].i]<=0||hv[es[eid].i]!=hv[es[eid].j]+1) return false;
    int dt;
    if(es[eid].c>es[eid].f){
        dt=min(ev[es[eid].i], es[eid].c-es[eid].f);
        es[eid].f+=dt, es[es[eid].op].f+=dt;
        ev[es[eid].i]-=dt, ev[es[eid].j]+=dt;
        return true;
    }
    if(es[eid].c==-1&&es[eid].f>0){
        dt=min(ev[es[eid].i], es[eid].f);
        es[eid].f-=dt, es[es[eid].op].f-=dt;//only reverse flow here!
        ev[es[eid].i]-=dt, ev[es[eid].j]+=dt;
        return true;
    }
    return false;
}
bool relabel(int u){//relabel for node u
    int hmin=INT_MAX; bool minV=true;
    for(int e=first[u]; e!=-1; e=es[e].next){
        if(es[e].c>es[e].f||(es[e].c==-1&&es[e].f>0)){//in the residental graph
            if(hv[es[e].j]<hv[u]) minV=false;
            if(hv[es[e].j]<hmin) hmin=hv[es[e].j];
        }
    }
    if(minV&&hmin<=n){ hv[u]=hmin+1; return true;}
    return false;
}
void dfsforS(int v){
    inS[v]=true;
    for(int e=first[v]; e!=-1; e=es[e].next){
        if(es[e].c>es[e].f){//forward-edge
            dfsforS(es[e].j);
        }
    }
}
#define fbase
int main()
{
#ifdef fbase
    ifstream inf("input.txt");
    if(!inf) return EXIT_FAILURE;
    cin.rdbuf(inf.rdbuf());
#endif
    //implementation
    int cap, t;
    cin>>in>>en, n=in+en+1;
    eid=0, memset((void*)first, 0xff, sizeof(first)),
    memset((void*)inS, 0, sizeof(inS));
    for(int i=1; i<=en; i++) cin>>cap, addEdge(in+i, n, cap);
    for(int i=1; i<=in; i++) cin>>cap, addEdge(0, i, cap);
    for(int i=1; i<=en; i++){
        cin>>cap;
        for(int j=1; j<=cap; j++) cin>>t, addEdge(t, i+in, INT_MAX);
    }
    //finish create map;
    Hpreflow();
    dfsforS(0);
    int maxflow=0;
    for(int e=first[0]; e!=-1; e=es[e].next){
        maxflow+=es[e].f;
    }
    cout<<"total profit : "<<maxflow<<endl;
    for(int i=1; i<n; i++){
        if(!inS[i]){
            if(i<=in) cout<<"I"<<i<<" ";
            else cout<<"E"<<(i-in)<<" ";
        }
    }
    cout<<endl;
#ifdef fbase
    inf.close();
#endif
    return 0;
}
#undef fbase

最大流和二分理论

Leave a comment

  • 最大流的四个算法框架,理论不需要解释了!原来实现的代码写的很烂,今天结合了《实用算法》上的实现来将算法实现和调试通过了,明天上午估计会耽搁一下,下午回来就继续把LRJ书上的题目一道一道的过掉。希望自己不看今天的代码框架一次做完,这些框架应该要非常熟悉才对!

1。一点改动:对于f(u,v)和f(v,u)而言,正负号判断让我很晕菜!先屏蔽到无效输入capcity=0的节点,OK,这样的节点对流的状态是不会有影响的。然后,仿照书上的例子把cap=-1时设置为标志位,若cap==-1则为反向边,Ef中g(v,u)=f(u,v),在f(u,v)>0时;正向边,Ef中g(u,v)=c(u,v)-f(u,v)。行,下面就可以实现了。

2。预处理头:

  • 公共的添边操作

#include <stdlib.h>
#include <iostream>
#include <fstream>

using namespace std;
const int MAXN=20, MAXE=400;

typedef struct Edge{
 int i, j, next, op, f, c;
};
Edge es[MAXE];
int eid=0, vid=0, first[MAXN], first1[MAXN], v, e, queue[MAXE], pre[MAXN], ptr[MAXN], hv[MAXN], ev[MAXN];//extra
int l[MAXN];
bool visited[MAXN];
//queue for bfs algorigthm, pre for each node in the queue, ptr for linked edge of each node!

//declaration
void addEdge(int i, int j, int c);

//common operation for adding edge, using int Edmond-Karp and Dicnic algorithm
void addEdge(int i, int j, int c){//i->j, and i<-j
 if(!c) return;
 es[eid].i=i, es[eid].j=j, es[eid].f=0, es[eid].c=c;
 es[eid].next=first[i], first[i]=eid;
 es[eid].op=eid+1, eid++;//set the opposite edge
 es[eid].i=j, es[eid].j=i, es[eid].f=0, es[eid].c=-1;//for opposite edge, default setting for capcity
 es[eid].next=first[j], first[j]=eid;
 es[eid].op=eid-1, eid++;
}

  • 残图更新逻辑:由于正向边反向边的更新逻辑一致,因此对于正向流,则f(u,v)+=dt(u,v),f(v,u)+=dt(u,v);而负向流,则为f(u,v)-=dt(u,v),f(v,u)-=dt(u,v)

3。主要算法详细实现

  • Ford-Fulkerson系最短路框架
算法1:基于Edmon-Karp算法框架(BFS的思路,使用BFS来生成最短路,算导上的表示)
算法伪码框架:[下午整理]
 
Edmond_Karp 算法:由于BFS的逻辑是基于一条随机最短路,因此其效率而言是四个算法中比较低的,但是有的时候也是一种可行算法。
     //初始化最大流的处理边逻辑
     while(find argumenting path from s to t in the residental graph){
             //寻找其路径最小值瓶颈Cf(u,v)
             //更新所有相关边,注意残图更新逻辑
     }
     //求出s的所有流出边,求出t的所有流入边,两者应该是一致的,均等于最大流的绝对值大小!
 
基于框架的实现代码:[已调试通过,复杂度O(E^2V)]

//implementation 1 : Edmond-Karp algorithm
void Edmond_Karp();

//simple implementation with BFS algorithm, O(VE^2).
void Edmond_Karp(){
 //create graph
 ifstream in("input.txt");
 if(!in) exit(EXIT_FAILURE);
 cin.rdbuf(in.rdbuf());
 cin>>v>>e;
 int start, end, cap;
 memset((void*)first, 0xff, sizeof(first));
 for(int i=0; i<e; i++){
  cin>>start>>end>>cap, addEdge(start, end, cap);
 }
 //BFS
 int top=0, size=1, maxflow=0, fmin;
 bool reachable=false;
 do{
  //BFS to find argumenting path from 0 to n-1 node,pre[1…n-1]记录其最大流边
  reachable=false, memset((void*)pre, 0xff, sizeof(pre));
  queue[top]=0, top=0, size=1, pre[top]=INT_MAX;
  while(top<size){
   int cur=queue[top++];//current node is v
   if(cur==v-1){reachable=true; break;}
   //here: orlando please dwell on why make the same mistake again! for node pre[es[i].j]!=pre[i];
   for(int i=first[cur]; i!=-1; i=es[i].next){
    if(pre[es[i].j]==-1&&(es[i].c>es[i].f||(es[i].c==-1&&es[i].f>0)))//not-visited and in the residental graph,残图判断逻辑
     pre[es[i].j]=cur, ptr[es[i].j]=i, queue[size++]=es[i].j;
   }
  }
  if(!reachable) break;
  //could reach n-1, backtrack to found path, and do updating
  //to find the fmin, one path backtrack
  fmin=INT_MAX;
  int edge=-1;
  for(int cur=v-1; cur!=0; cur=pre[cur]){//for v’s value
   if(es[ptr[cur]].c>es[ptr[cur]].f)
    if(es[ptr[cur]].c-es[ptr[cur]].f<fmin) fmin=es[ptr[cur]].c-es[ptr[cur]].f, edge=ptr[cur];
    else if(es[ptr[cur]].c==-1&&es[ptr[cur]].f>0)
     if(es[ptr[cur]].f<fmin) fmin=es[ptr[cur]].f, edge=ptr[cur];
  }
  //updating the value of agrumenting path,更新最大流的前向状态
  for(int cur=v-1; cur!=0; cur=pre[cur]){
   if(es[ptr[cur]].c>es[ptr[cur]].f)//正向流
    es[ptr[cur]].f+=fmin, es[es[ptr[cur]].op].f+=fmin;
   else if(es[ptr[cur]].c==0&&es[ptr[cur]].f>0)//反向流
    es[ptr[cur]].f-=fmin, es[es[ptr[cur]].op].f-=fmin;
  }
  maxflow+=fmin;
 }while(reachable);
 int t1=0, t2=0;
 for(int tmp=first[v-1]; tmp!=-1; tmp=es[tmp].next){
  t1+=es[tmp].f;
 }
 for(int tmp=first[v-1]; tmp!=-1; tmp=es[tmp].next){
  t2+=es[tmp].f;
 }
 cout<<"source ouput f : "<<t1<<endl;
 cout<<"dest input f : "<<t2<<endl;
 cout<<maxflow<<endl<<endl;
 in.close();
}
算法2:基于Dicnic算法框架(BFS生成标号,严格DFS的标号检查,算法书上有错,俄!让我郁闷了一小会儿)
算法伪码框架:[下午整理]
     //初始化最大流的处理边逻辑
     while(true){
             //bfs对每一个顶点进行距离标号,d[0]=0
             //若对于目标点无法标号,在bfs结束时,退出循环;
             //dfs求出多条从s到t的容许路径,d[i]=d[j]+1,求出最大流瓶颈值,更新所有流状态
             //退流是退到满流的上一个顶点,这样可以一次找到多个增广路
     }
     //求出s的所有流出边,求出t的所有流入边,两者应该是一致的,均等于最大流的绝对值大小! 
 
基于框架的实现代码:[已调试通过,复杂度O(V^2E)]
//implementation 2 : Dinic algorithm
void Dicnic();
void bfsLevel();//similar to d[i] label for minmal shortest path.bfs标号实现
void dfs_maxflow();//dfs的最大流实现
 
void Dicnic(){//refer to practise of Dicnic algorithm.
 //create graph, and cin binding operation
 ifstream in("input.txt");
 if(!in) exit(EXIT_FAILURE);
 cin.rdbuf(in.rdbuf());
 cin>>v>>e;
 int start, end, cap;
 eid=0, memset((void*)first, 0xff, sizeof(first));
 for(int i=0; i<e; i++){
  cin>>start>>end>>cap, addEdge(start, end, cap);
 }
 while(true){
  bfsLevel();
  if(pre[v-1]==INT_MAX) break;//n-1 with no dist[i] in queue
  dfs_maxflow();//dfs_maxflow(&maxflow), accumulate together 
 }
 int t1=0, t2=0;
 for(int tmp=first[0]; tmp!=-1; tmp=es[tmp].next){
  t1+=es[tmp].f;
 }
 cout<<"source ouput f : "<<t1<<endl;
 for(int tmp=first[v-1]; tmp!=-1; tmp=es[tmp].next){
  t2+=es[tmp].f;
 }
 cout<<"dest input f : "<<t2<<endl<<endl;
 in.close();
}
void bfsLevel(){//bfs for level, recording in pre[0…v-1] array
 for(int i=0; i<v; i++) pre[i]=INT_MAX;
 int top=0, size=1; queue[top]=0, pre[0]=0;
 while(top<size){
  int cv=queue[top++];
  if(cv==v-1) return;//pre[cv] has been set previously, similar to my code
  for(int ei=first[cv]; ei!=-1; ei=es[ei].next){
   if(pre[es[ei].j]>pre[cv]+1){//smart code
    if(es[ei].c>es[ei].f||(es[ei].c==-1&&es[ei].f>0))//add into queue
     queue[size++]=es[ei].j, pre[es[ei].j]=pre[cv]+1;
   }
  }
 }
}
//core implemenation : mix two logic together
void dfs_maxflow(){//书上逻辑有错,找到最大流后反向更新即可!
 memcpy((void*)first1, (void*)first, sizeof(first));
 memset((void*)visited, 0, sizeof(visited));
 int top=0; queue[0]=0;
 while(top>-1){
  if(queue[top]==v-1){//dfs find the n-1
   int fmint=INT_MAX, edge, tj;
   for(int j=top; j>0; j–){//down to 1
    if(es[ptr[queue[j]]].f<es[ptr[queue[j]]].c){
     if(es[ptr[queue[j]]].c-es[ptr[queue[j]]].f<fmint) fmint=es[ptr[queue[j]]].c-es[ptr[queue[j]]].f;
    }
    else if(es[ptr[queue[j]]].c==-1&&es[ptr[queue[j]]].f>0){
     if(es[ptr[queue[j]]].f<fmint) fmint=es[ptr[queue[j]]].f;
    }
   }
   for(int j=top; j>0; j–){
    if(es[ptr[queue[j]]].f<es[ptr[queue[j]]].c){
     es[ptr[queue[j]]].f+=fmint, es[es[ptr[queue[j]]].op].f+=fmint;
     if(es[ptr[queue[j]]].f==es[ptr[queue[j]]].c) tj=j;
    }
    else if(es[ptr[queue[j]]].c==0&&es[ptr[queue[j]]].f>0){
     es[ptr[queue[j]]].f-=fmint, es[es[ptr[queue[j]]].op].f-=fmint;     
     if(es[ptr[queue[j]]].f==0) tj=j;
    }
   }
   return;//比较原始的方法,每次只找一条增广路经
   top=tj-1;//回退到满流的上一个顶点,再次增广,Dicnic就这个比较valuable
  }
  else{
   int tmp=first1[queue[top]], tp=queue[top];
   for(;tmp!=-1; tmp=es[tmp].next){
    if(!visited[es[tmp].j]&&(pre[es[tmp].j]==pre[tp]+1)){
     if(es[tmp].f<es[tmp].c||(es[tmp].c==-1&&es[tmp].f>0)){
      first1[tp]=es[tmp].next;//modify next link, dfs recursive have to modify edge linkage
      queue[++top]=es[tmp].j, ptr[queue[top]]=tmp;
      break;
     }//if
    }
   }//for
   if(tmp==-1){visited[tp]=true; top–;}
  }
 }
}
 
  • 预流标进算法

通用的算法操作:

bool Relabel(int v);
bool Push(int e);

//实现代码 :对于余流节点的重标和push算法

bool Relabel(int u){//if we can relabel the node
 if(ev[u]<=0) return false;
 //in the residental graph
 int hmin=INT_MAX; bool minV=true;//最小高度,是否是最低点
 for(int e=first[u]; e!=-1; e=es[e].next){
  if(es[e].c>es[e].f||(es[e].c==-1&&es[e].f>0)){//in the residental graph
   if(hv[es[e].j]<hv[u]) minV=false;//hv[u]<=hv[es[e].j]
   if(hv[es[e].j]<hmin) hmin=hv[es[e].j];

  }
 }
 if(minV&&hmin<=v){hv[u]=hmin+1; return true; }//注意,这里v比|V|+1大一个距离,我们要考虑到比源点高度更高,这样所有的余流才有可能回流到源点之上.
 return false;
}
bool Push(int e){
 if(ev[es[e].i]<=0) return false;
 if(hv[es[e].i]==hv[es[e].j]+1){//for push operation
  int dt;
  if(es[e].f<es[e].c){
   dt=min(ev[es[e].i], es[e].c-es[e].f);
   es[e].f+=dt, es[es[e].op].f+=dt;
   ev[es[e].i]-=dt, ev[es[e].j]+=dt;
   return true;
  }
  if(es[e].c==-1&&es[e].f>0){
   dt=min(ev[es[e].i], es[e].f);
   es[e].f-=dt, es[es[e].op].f-=dt;
   ev[es[e].i]-=dt, ev[es[e].j]+=dt;//reverse-equal algorithm
   return true;
  }
 }
 return false;
}

算法3:基于Goldberg算法框架(简单的预流处理,基于顶点的随机验证)
算法伪码框架:[下午整理] 
     //初始化最大流的处理边逻辑,初始化源点的流边和顶点余流
     while(对于每一个可行顶点能找到relabel,或者push操作){
             //relabel或者push操作();
     }
     //求出s的所有流出边,求出t的所有流入边,两者应该是一致的,均等于最大流的绝对值大小!
 
基于框架的实现代码:[已调试通过,复杂度O(V^2E)]
//for random implemenation
void Goldberg(){
 ifstream in("input.txt");
 if(!in) exit(EXIT_FAILURE);
 cin.rdbuf(in.rdbuf());
 cin>>v>>e;
 int start, end, cap;
 eid=0, memset((void*)first, 0xff, sizeof(first));
 for(int i=0; i<e; i++){
  cin>>start>>end>>cap, addEdge(start, end, cap);
 }
 Initial();
 bool t=true;
 do{
  t=true;
  for(int i=1; i<v-1; i++){
   if(Relabel(i)) t=false;
   for(int e=first[i]; e!=-1; e=es[e].next)
    if(Push(e)) t=false;
  }
 }while(!t);
 int t1=0, t2=0;
 for(int tmp=first[0]; tmp!=-1; tmp=es[tmp].next){
  t1+=es[tmp].f;
 }
 cout<<"source ouput f : "<<t1<<" "<<ev[0]<<endl;
 for(int tmp=first[v-1]; tmp!=-1; tmp=es[tmp].next){
  t2+=es[tmp].f;
 }
 cout<<"dest input f : "<<t2<<" "<<ev[v-1]<<endl<<endl;
 in.close();
}
void Initial(){
 memset((void*)hv, 0, sizeof(hv)), memset((void*)ev, 0, sizeof(ev));
 hv[0]=v-1;//from 0 to v-1
 for(int e=first[0]; e!=-1; e=es[e].next){
  if(es[e].c>0){
   es[e].f+=es[e].c, es[es[e].op].f+=es[e].c;//for 0’s output stream
   ev[es[e].j]+=es[e].c, ev[0]-=es[e].c;
  }
 }
}
 
算法4:基于群牛算法框架(最高标号,预流标进)
算法伪码框架:[下午整理]
     //初始化最大流的处理边逻辑
     //初始化L随机队列(除开0,以及v-1),由于在最大流边形成之中,已经生成了相邻边逻辑,因此下面可以省去!
     while(u!=-1){
           int oldheight=hv[u];
           Discharge(u);//核心操作discharge(u)
           if(oldheight变化) 将u调整到L队列头
           u=L-next(u);//记录其边位置即可
     }
     //求出s的所有流出边,求出t的所有流入边,两者应该是一致的,均等于最大流的绝对值大小!
 
基于框架的实现代码:[已调试通过,复杂度O(V^3)->渐进意义:O(V^2M^1/2)]
//implemenation 4 : HighestPreflow
void Relabel_to_front();
void Initialize_preflow();
void Discharge(int u);//if old-height can be changed any more!
//notes: HighestPreflow implemented by linklist
void Relabel_to_front(){
 Initialize_preflow();
 int pos=0, u=l[pos];
 while(u!=-1){
  int oldh=hv[u];
  Discharge(u);
  if(hv[u]>oldh){
   //switch pos and 0
   if(pos!=0) l[pos]^=l[0], l[0]^=l[pos], l[pos]^=l[0];
   pos=0;
  }
  u=l[++pos];
 }
 int t1=0, t2=0;
 for(int tmp=first[0]; tmp!=-1; tmp=es[tmp].next){
  t1+=es[tmp].f;//非常重要的性质:所有中间结点的ev[es[e].j]==0
 }
 cout<<"source ouput f : "<<t1<<" "<<ev[0]<<endl;
 for(int tmp=first[v-1]; tmp!=-1; tmp=es[tmp].next){
  t2+=es[tmp].f;
 }
 cout<<"dest input f : "<<t2<<" "<<ev[v-1]<<endl<<endl;
}
void Discharge(int u){//discharge for u
 int tmp=first[u];
 while(ev[u]>0){
  if(tmp==-1) Relabel(u), tmp=first[u];//for relabel
  else if((es[tmp].c>es[tmp].f||(es[tmp].c==-1&&es[tmp].f>0))&&hv[u]==hv[es[tmp].j]+1)
   Push(tmp);
  else tmp=es[tmp].next;
 }
}
void Initialize_preflow(){
 ifstream in("input.txt");
 if(!in) exit(EXIT_FAILURE);
 cin.rdbuf(in.rdbuf());
 cin>>v>>e;
 int start, end, cap;
 eid=0, vid=0, memset((void*)first, 0xff, sizeof(first));
 for(int i=0; i<e; i++){
  cin>>start>>end>>cap, addEdge(start, end, cap);
 }
 memset((void*)l, 0xff, sizeof(l));
 for(int i=0; i<v-2; i++) l[i]=i+1;
 //for hv and ev
 memset((void*)hv, 0, sizeof(hv)), memset((void*)ev, 0, sizeof(ev));
 hv[0]=v-1;
 for(int e=first[0]; e!=-1; e=es[e].next){
  if(es[e].c>0){
   es[e].f+=es[e].c, es[es[e].op].f+=es[e].c;
   ev[es[e].j]+=es[e].c, ev[0]-=es[e].c;
  }
 }
 in.close();
}
  • 以后的算法,主要选用Dicnic最高预流标进,所以请熟练掌握,如果不能一次写完无bugs,请删掉重来!

最短路的题目

Leave a comment

题目很简单,就是求一个带点权最大的最短路,lrj说的很清楚可以用Dijkstra加枚举点处理;在UVA的讨论上看到在问Floyd-Warshall的算法如何用在这道题目上!应该是可以的,只是就像lrj说的一样,多次查询打表是比较好的策略。其他没有什么,只是敲代码的时候,忘记了f(tmp!=INT_MAX&&tmp+pcost[j]<pcost[k]) pcost[k]=tmp+pcost[j]; 开始写成tmp了,眼花看了半天,结果还从头到尾把算法过了一遍!这个毛病和原来奥数比赛时候一样,晃得一脉相承,要注意了,敲代码要精神集中!有个很明显的改进,lrj说要把比X贵的城市删除,这里直接用qsort排序,然后反向求其范围,不算改进吧,常识的处理。

另外就是一点感概,UVA确实比PKU要严谨,多一个换行都是WA,要提醒自己多注意细节的问题!细节决定成败!

My Submissions

#

Problem

Verdict

Language

Run Time

Submission Date

7082960 10246  Asterix and Obelix 

Accepted

C++

0.430

2009-04-21 03:33:21
7082927 10246 Asterix and Obelix 

Wrong answer

C++

0.450

2009-04-21 02:47:30

代码如下:

#include <iostream>
#include <fstream>

//UVA 10246 – Asterix and Obelix
using namespace std;
const int MAXC=81, MAXR=1001;//MAX Edges
typedef struct Edge{
    int i, j, cost, next1, next2;//for link to other i->next edges
};
typedef struct City{
    int n, cost;
};
Edge es[MAXR];
int first[MAXC],LC[MAXC][MAXC], city[MAXC], pcost[MAXC], eid=0, c, r, q;//or 0xff should be the maxmimal
City CC[MAXC];
bool inG[MAXC], inD[MAXC];

//declaration
void addEdge(int s, int e, int cost);
int cmp(const void* p1, const void* p2);
int getCost(int x, int y);
void DijsktraEach(int i);
void Dijkstra(int v);//v is the starting city, to all cities inG[1…n] is true;

//implementation
void DijsktraEach(int i){//first of all, earse all cities
    memset((void*)inG, 0, sizeof(inG));
    for(int k=i; k<=c; k++) inG[CC[k].n]=true;//make in city
    int x=CC[i].n; LC[x][x]=CC[i].cost;//sit for city x’s self
    if(i!=c) Dijkstra(x);//order correctly
}
void Dijkstra(int v){
    //the last value is for +LC[x][x]
    int vn=0;
    for(int i=1; i<=c; i++){
        if(!inG[i]){inD[i]=true; continue;}
        inD[i]=false; pcost[i]=INT_MAX; vn++;
    }
    for(int x=first[v]; x!=-1; ){//find the next slot
        int t1=es[x].i, t2=es[x].j;
        //next logic
        if(t1==v) pcost[t2]=es[x].cost, x=es[x].next1;
        else if(t2==v) pcost[t1]=es[x].cost, x=es[x].next2;//just neighbour could add LC[v][v], error
    }
    inD[v]=true;//OK, let’s for dijstra
    int tmp=INT_MAX, j=v;
    for(int i=1; i<vn; i++){
        tmp=INT_MAX;
        for(int k=1; k<=c; k++){
            if(!inD[k]&&pcost[k]<tmp) tmp=pcost[k], j=k;
        }
        if(inD[j]) continue;//if j==v, avoid exception
        inD[j]=true; LC[v][j]=LC[j][v]=pcost[j]+LC[v][v];
        //update the other points in the sets
        for(int k=1; k<=c; k++){
            if(!inD[k]){
                int tmp=getCost(k, j);
                if(tmp!=INT_MAX&&tmp+pcost[j]<pcost[k]) pcost[k]=tmp+pcost[j];//mimnal
            }
        }
    }
}
int getCost(int x, int y){
    int cost=INT_MAX;
    for(int x1=first[x]; x1!=-1; ){
        int t1=es[x1].i, t2=es[x1].j;
        if(t1==x){
            if(y==t2) cost=es[x1].cost;
            x1=es[x1].next1;
        }
        else if(t2==x){
            if(y==t1) cost=es[x1].cost;
            x1=es[x1].next2;
        }
    }
    return cost;
}
void addEdge(int s, int e, int cost){
    es[eid].i=s, es[eid].j=e, es[eid].cost=cost;
    es[eid].next1=first[s], first[s]=eid;//bi-direction!yeah, don’t forget
    es[eid].next2=first[e], first[e]=eid;//next1 for i, next2 for j.
    eid++;
}
int cmp(const void* p1, const void* p2){//deceasing order, remember the order
    return ((City*)p1)->cost<((City*)p2)->cost;
}

#define fbase
int main()
{
#ifdef fbase
    ifstream in("input.txt");
    if(!in) return EXIT_FAILURE;
    cin.rdbuf(in.rdbuf());
#endif

    //implementation
    int s, e, cost, c1, c2, id=1;
    cin>>c,cin>>r,cin>>q;
    do{
        if(c==0&&r==0&&q==0) break;//quit the game
        eid=0;memset((void*)first, 0xff, sizeof(first));
        for(int i=1; i<=c; i++) for(int j=1; j<=c; j++) LC[i][j]=INT_MAX;//intial index
        for(int i=1; i<=c; i++) cin>>CC[i].cost, CC[i].n=i;//cost for each city
        for(int i=1; i<=r; i++){
            cin>>s, cin>>e, cin>>cost;
            addEdge(s, e, cost);//set for first and next value
        }
        //qsort for city
        qsort((void*)(CC+1), c, sizeof(City), cmp);
        for(int i=1; i<=c; i++){//from CC[i].n->CC[i].cost decending order, to generate DP table
            DijsktraEach(i);
        }
        cout<<"Case #"<<(id++)<<endl;
        for(int i=1; i<=q; i++){
            cin>>c1, cin>>c2;
            int lcost=INT_MAX;
            for(int j=1; j<=c; j++){
                if(LC[c1][j]==INT_MAX||LC[c2][j]==INT_MAX) continue;
                int tmp=LC[c1][j]+LC[c2][j]-LC[j][j];
                if(tmp<lcost) lcost=tmp;
            }
            if(lcost==INT_MAX) cout<<-1<<endl;
            else cout<<lcost<<endl;
        }
        cin>>c,cin>>r,cin>>q;
        if(c!=0||r!=0||q!=0) cout<<endl;
    }while(c!=0||r!=0||q!=0);

#ifdef fbase
    in.close();
#endif
    return 0;
}
#undef fbase

  • 最短路算法的总结和应用
两个类型:一种是标号修正法,思想和Dinkench算法一样,不断逼近最终达到其上下限;第二种是标号设定算法,对临时标号不断修改,然后不能修改的最优就是固定标号!
PS:开始小瞧了这个思想,结果搞得自己头晕眼花的,其实标号设定算法也可以看作是一种拓扑序的计算过程;只是这是的拓扑序不像一般的图一样是基于边的变换,而是基于最终值的一个序列,想通了这点下面的题目就简单了;没想通就有苦头吃了!
 
1,Crossing the desert,这是一道ACM/ICPC的world fina 2002l的题目,和楼教主他们做的题目不一样,这个题估计要被教主在15分钟之内秒杀哟!
两点容易错:题目给的是int,但实际上要用double来做,读题不细心,卡了好久;还有就是这里所有临时标号都是基于当前最优标号值得,这和dijkstra的思路有些不同,生搬硬套,吃了苦头!顺便练习了下循环队列的写法,在没有carlo提示之前,我一直以为是自己算法有问题,换了好几种算法,不是WA就是TLE,哭死了,结果是精度问题!下面是AC的代码:
 
  1. #include <iostream>
  2. #include <cmath>
  3.  
  4. using namespace std;
  5. const int MAXN=31, MAXS=1000004;
  6. struct Oase{
  7.     double x,y;
  8. };
  9. double dis[MAXN][MAXN], d[MAXN], cap;
  10. int n;
  11. Oase os[MAXN];
  12. bool inG[MAXN];
  13.  
  14. int main()
  15. {
  16.     //implementation
  17.     double tmp, m;
  18.     int k, sum, id=1;
  19.     while(cin>>n&&cin>>cap&&n&&cap){
  20.         for(int i=1; i<=n; i++) cin>>os[i].x, cin>>os[i].y, inG[i]=false, d[i]=MAXS;
  21.         for(int i=1; i<=n; i++)
  22.             for(int j=i+1; j<=n; j++)
  23.                 dis[i][j]=dis[j][i]=sqrt((os[i].x-os[j].x)*(os[i].x-os[j].x)+(os[i].y-os[j].y)*(os[i].y-os[j].y));
  24.         d[n]=0; inG[n]=true;
  25.         for(int i=1; i<n; i++){
  26.             m=2*dis[n][i]-cap;
  27.             if(m>0&&cap-3*dis[n][i]>0){
  28.                 tmp=(ceil(m/(cap-3*dis[n][i]))*2+1)*dis[n][i];
  29.                 if(tmp<d[i]) d[i]=tmp;
  30.             }
  31.             else if(cap-2*dis[n][i]>=0){
  32.                 tmp=dis[n][i];
  33.                 if(tmp<d[i]) d[i]=tmp;
  34.             }
  35.         }
  36.         //dijkstra
  37.         for(int i=1; i<n; i++){
  38. tmp=MAXS, k=-1;
  39.             for(int i=1; i<n; i++){
  40.                 if(!inG[i]&&d[i]<tmp) tmp=d[i], k=i;
  41.             }
  42.             if(k==-1) break;
  43.             inG[k]=true;
  44.             for(int i=1; i<n; i++){
  45.                 if(!inG[i]){
  46.                     m=d[k]-cap+2*dis[i][k];
  47.                     if(m>0&&cap-3*dis[i][k]>0){
  48.                         tmp=(ceil(m/(cap-3*dis[i][k]))*2+1)*dis[i][k]+d[k];
  49.                         if(tmp<d[i]) d[i]=tmp;
  50.                     }
  51.                     else if(cap-2*dis[i][k]-d[k]>=0){
  52.                         tmp=dis[i][k]+d[k];
  53.                         if(tmp<d[i]) d[i]=tmp;
  54.                     }
  55.                 }
  56.             }
  57.         }
  58.         sum=(int)ceil(d[1]);
  59.         if(sum==0) sum++;
  60.         if(sum>1000000)
  61.             cout<<"Trial "<<(id++)<<": Impossible "<<endl<<endl;
  62.         else
  63.             cout<<"Trial "<<(id++)<<": "<<sum<<" units of food"<<endl<<endl;
  64.     }
  65.     return 0;
  66. }
2.UVA 10381 The Rock,这道题目让我异常痛苦,换了三种算法,两个TLE,一个后面有个严重的漏洞!现在都还没有过,最好的跑了5s多,最差的跑了18s。还是说下思路吧,以后有好的思路再来弄它!
算法分析:最坏的最短时间设为d[i], i=1…x*y,问题在于当前的一个石头是不可知的;由于d[i]的序列是严格增加的,在出口处为d[exit]=0,反向到入口求出d[extrance],那么结果就可以知道了。下面有点像动规,对于当前要求得点i,如果它在d[i]求解序列的拓扑上一个点的d是可知的,那么一定有以下两种情况:
  • 若上一个拓扑序的点为石头,那么可以求一条新的最短路到exit,这个时候值为distance[i, A],i为石头,A为当前点,表示i为石头的时A到出口的最短路长;
  • 若为空白,加个1,即可。

至于如果求出d的严格增加序列,一个最小堆处理即可。下面是实现代码,算法瓶颈在于distance[i, A],求得痛苦无比,开始放在预处理里面,跑料18s;然后用DP打表,就是现在的版本了,还是很慢!

#include <iostream>
#include <fstream>

using namespace std;
//UVA 10381 – The Rock
const int MAXN=1600, MAXD=4;
char map[MAXN];//0 for space, 1 for wall;
int r, c, exits, b, extrance, size, dis[MAXN][MAXN],/*pre-calculation*/
    d[MAXN], cost[MAXN], queue[MAXN];
int dx[4]={0,0,-1,1}, dy[4]={-1,1,0,0};//left 0, right 1, up 2, down 3
bool inG[MAXN], inQ[MAXN];

//declaration
int Dijkstra(int stone);//set stone position as wall, then use dijkstra to calculate shortest path.
void backtrack();//set d[exit]=0, backtrack to d[extrance].
//heap implemenation
void pushHeap(int cur);
int popHeap();
int cmp(const int& p1, const int& p2);

void inline pushHeap(int cur){
    if(size==MAXN) return;
    queue[size++]=cur;
    if(size>1) push_heap(queue, queue+size, cmp);
}
int inline popHeap(){
    if(size==0) return -1;
    if(size>1) pop_heap(queue, queue+size, cmp);
    size–; return size;
}
int cmp(const int& p1, const int& p2){
    return d[p1]>d[p2];
}

//implementation
int Dijkstra(int stone, int v){
    if(dis[stone][v]!=-1) return dis[stone][v];
    map[stone]=1;
    //from exit point backtrack to all reachable points
    for(int i=0; i<b; i++){
        cost[i]=INT_MAX;
        if(map[i]==1) inG[i]=true;//for stone check
        else inG[i]=false;
    }
    inG[exits]=true, cost[exits]=0, dis[stone][exits]=0;
    //set 4 direction for exit
    int x=exits/c, y=exits%c;
    for(int i=0; i<4; i++){
        if(x+dx[i]>-1&&x+dx[i]<r&&y+dy[i]>-1&&y+dy[i]<c){//no stone
            int p=(x+dx[i])*c+y+dy[i];
            if(!inG[p]) cost[p]=1;//take care of
        }
    }
    //Dijstra algorithm
    int tmp, k;
    for(int i=1; i<b; i++){//for add i’s vertexs
        //find the mimal amount
        tmp=INT_MAX, k=-1;
        for(int j=0; j<b; j++){//?can we optimize it
            if(!inG[j]&&cost[j]<tmp) tmp=cost[j], k=j;
        }
        if(k==-1) break;
        inG[k]=true, dis[stone][k]=cost[k];
        if(k==v) {map[stone]=0; return dis[stone][v];}
        x=k/c, y=k%c;
        for(int j=0; j<4; j++){//O(1)
            if(x+dx[j]>-1&&x+dx[j]<r&&y+dy[j]>-1&&y+dy[j]<c){
                int p=(x+dx[j])*c+(y+dy[j]);
                if(!inG[p]&&cost[k]+1<cost[p]) cost[p]=cost[k]+1;
            }
        }
    }
    map[stone]=0; return dis[stone][v];
}
void backtrack(){
    memset((void*)inQ, 0, sizeof(bool)*b);
    size=0; pushHeap(exits), inQ[exits]=true, d[exits]=0;
    while(size){
        int cur=queue[popHeap()];
        //update 4’s direction for grid
        int x=cur/c, y=cur%c, p, tmp;
        for(int i=0; i<4; i++){//cur is adjacent to p
            if(x+dx[i]>-1&&x+dx[i]<r&&y+dy[i]>-1&&y+dy[i]<c){
                p=(x+dx[i])*c+(y+dy[i]);
                if(!inQ[p]&&map[p]==0){//possible
                    tmp=d[cur]+1;
                    int t1=(dis[cur][p]==-1?Dijkstra(cur, p):dis[cur][p]);
                    if(cur!=exits&&cur!=extrance&&t1!=-1
                        &&dis[cur][p]>tmp) tmp=t1;
                    if(tmp>d[p]&&d[p]==-1) d[p]=tmp;
                    //using priority_queue instead of circlar queue.
                    if(p==extrance) return;
                    pushHeap(p), inQ[p]=true;
                }
            }
        }
    }
}

#define fbase
int main()
{
#ifdef fbase
    ifstream in("input.txt");
    if(!in) return EXIT_FAILURE;
    cin.rdbuf(in.rdbuf());
#endif

    //implementation
    char ch; int ncase;
    cin>>ncase;
    while(ncase–){
        cin>>r, cin>>c;
        b=r*c;
        for(int i=0; i<b; i++){
            cin>>ch;
            if(ch==’.’) map[i]=0;
            else if(ch==’#’) map[i]=1;
            else if(ch==’E’) map[i]=0, extrance=i;
            else if(ch==’X’) map[i]=0, exits=i;
        }
        memset((void*)dis, 0xff, sizeof(dis)),
        memset((void*)d, 0xff, sizeof(d));//clr
        backtrack();
        cout<<d[extrance]<<endl;
    }

#ifdef fbase
    in.close();
#endif
    return 0;
}
#undef fbase

  • 第二个思路是,不管三七二十一先求出最短路,注意这里的特殊情况,第一是一个正方形网格,那么从一点到另一点无重叠的最短路最多两条,题目上说是只有一条那么就简单多了。一点没有相通的问题:我们如何求出在最短路上有石头时,新的最短路上最坏最短时间呢?比如我可以求出初始最短路的最坏最短时间为17,OK,那么用了一条新的路,这条路路长为13,但假设石头在其上的时候,路长肯定要变大;这个解法是基于路的,换句话说就不是一个严格的d拓扑序列,不能保证在反向达到入口时求出的d值一定是最优,但要求出最优,也就是要枚举每个可能道路,直接疯掉!下面是写了一半卡做的代码[希望以后有时间能补上这个遗憾,或者向高手请教一下!]

#include <vector>
#include <iostream>
#include <fstream>

using namespace std;
//UVA 10381 – The Rock
const int MAXN=1600, MAXD=4;
char map[MAXN];//0 for space, 1 for wall;
int r, c, exits, b, extrance, size,
    sig, d[MAXN], cost[MAXN], buf[MAXN];
int dx[4]={0,0,-1,1}, dy[4]={-1,1,0,0};//left 0, right 1, up 2, down 3
bool inG[MAXN], inQ[MAXN];
vector<int> spath;

//declaration
int Dijkstra(int v, bool s);//set stone position as wall, then use dijkstra to calculate shortest path.
void backtrack();//set d[exit]=0, backtrack to d[extrance].
void tracepath(int v, int step);

//implementation
int Dijkstra(int v, bool s){
    //from exit point backtrack to all reachable points
    for(int i=0; i<b; i++){
        buf[i]=INT_MAX;
        if(map[i]==1) inG[i]=true;//for stone check
        else inG[i]=false;
    }
    inG[v]=true, buf[v]=0;
    //set 4 direction for exit
    int x=v/c, y=v%c;
    for(int i=0; i<4; i++){
        if(x+dx[i]>-1&&x+dx[i]<r&&y+dy[i]>-1&&y+dy[i]<c){//no stone
            int p=(x+dx[i])*c+y+dy[i];
            if(!inG[p]) buf[p]=1;//take care of
        }
    }
    //Dijstra algorithm
    int tmp, k;
    for(int i=1; i<b; i++){//for add i’s vertexs
        //find the mimal amount
        tmp=INT_MAX, k=-1;
        for(int j=0; j<b; j++){//?can we optimize it
            if(!inG[j]&&buf[j]<tmp) tmp=buf[j], k=j;
        }
        if(k==-1) break;
        inG[k]=true;
        if(s){
            if(k==exits) return buf[k];
        }
        x=k/c, y=k%c;
        for(int j=0; j<4; j++){//O(1)
            if(x+dx[j]>-1&&x+dx[j]<r&&y+dy[j]>-1&&y+dy[j]<c){
                int p=(x+dx[j])*c+(y+dy[j]);
                if(!inG[p]&&buf[k]+1<buf[p]) buf[p]=buf[k]+1;
            }
        }
    }
    if(!s) memcpy((void*)cost, (void*)buf, sizeof(int)*b);
    return -1;
}
void tracepath(int v, int step){
    if(sig) return;
    if(v==extrance){//buf as queue, d as direction
        for(int i=0; i<step; i++) spath.push_back(buf[i]);
        spath.push_back(extrance);
        sig=1; return;
    }
    int x=v/c, y=v%c; buf[step]=v, inG[v]=true;
    for(int i=0; i<4; i++){
        if(x+dx[i]>-1&&x+dx[i]<r&&y+dy[i]>-1&&y+dy[i]<c){
            int p=(x+dx[i])*c+(y+dy[i]);
            if(!inG[p]&&cost[p]==cost[v]+1){//not in the dfs and forward step
                tracepath(p, step+1);
            }
        }
    }
}

#define fbase
int main()
{
#ifdef fbase
    ifstream in("input.txt");
    if(!in) return EXIT_FAILURE;
    cin.rdbuf(in.rdbuf());
#endif

    //implementation
    char ch; int ncase, total;
    cin>>ncase;
    while(ncase–){
        cin>>r, cin>>c;
        b=r*c;
        for(int i=0; i<b; i++){
            cin>>ch;
            if(ch==’.’) map[i]=0;
            else if(ch==’#’) map[i]=1;
            else if(ch==’E’) map[i]=0, extrance=i;
            else if(ch==’X’) map[i]=0, exits=i;
        }
        memset((void*)d, 0xff, sizeof(d));//clr
        Dijkstra(exits, false), spath.clear();
        memset((void*)inG, 0, sizeof(bool)*b);//trace the shortest path
        sig=0, tracepath(exits, 0);
        d[exits]=0;
        int len=spath.size(), min=INT_MIN;
        for(int i=1; i<len; i++){//spath[i], d[spath[i-1].i]+1 and spath[i].i->exits
            int tmp=d[spath[i-1]]+1;
            if(i!=1){
                map[spath[i-1]]=1;
                int t1=Dijkstra(spath[i], true);
                //cout<<spath[i]<<" "<<t1;
                map[spath[i-1]]=0;
                if(t1>tmp) tmp=t1;
            }
            d[spath[i]]=tmp;
            //cout<<" "<<tmp<<endl;
        }
        cout<<d[extrance]<<endl;
    }

#ifdef fbase
    in.close();
#endif
    return 0;
}
#undef fbase

 

3. 最优双调路径,虽然这道最后,但个人觉得这道题目的思路很重要,让我想起来了处理二维最大和的算法(当然,那个还有一个预处理的技巧),一个模型不能处理的时候,可以考虑固定其一个参数,然后用通用模型来处理固定参数之后的模型!如果没有第二个路权参数,那么这就完全是个很菜的最优路,现在多了一个参数,那么我们就枚举这个路权即可。[怎么办呢?我参照了下标程,主要一个是被上面的题搞得晕头转向的,加上lrj说的不清不楚的,这下彻底菜了!一看代码,呵呵,原来是一个枚举加滚动数组,OK]

注意:

1。滚动数组的大小4*101即可,边长为101的上限,呵呵,很省空间!
2。jim哥哥的代码看得我好想吐哟,一m多的vector,最让我不爽的是这个time(n,k)表示n城市k费用到s时的最短时间!其实在c++里面有一种叫代理类的方式,可以帮助我们将time(n,k)转化为time[n][k],可惜!我比较菜,没有想到很好的传参的方式,写的严重破坏OO,没有看到Scott的实现怎么样的!
3。最后在写Dijkstra的时候,才发现lrj说的标号固定是个什么意思;我们先用0-1背包的思想生成k+1的最短旅行时间,由于0权费用路的存在,那么就用反复修正,输出不可变的固定标号!好的程序员说话就像在飘,太飘了,看来只有用代码说话了!NKU上直接切掉!
Solution User Problem Language Judge Result Memory Time Used Code Length Submit Time
261471 orlando22 10702 GNU C++ Accepted 1336KB 3343ms 5080B 2009-04-25 16:21:45.0
Source Code
#include <iostream>
#include <algorithm>

using namespace std;
const int INF=1000000000, MAXT=100, MAXE=300, MAXC=101;
struct Edge{
    int i, j, c, t, next1, next2;
};

int n, m, s, e, dist[MAXC], heap[MAXC];
bool inHeap[MAXC];

//forward-start for zero and non-zero t links
//declaration for edge implemenation
int eid=0, zeroC[MAXC], nzeroC[MAXC];
Edge es[MAXE];
void inline clrE();
void inline addEdge(int i, int j, int c, int t);

//implemenation
void inline clrE(){
    memset((void*)zeroC, 0xff, sizeof(zeroC)),
    memset((void*)nzeroC, 0xff, sizeof(nzeroC));
}
void inline addEdge(int i, int j, int c, int t){
    es[eid].i=i, es[eid].j=j, es[eid].c=c, es[eid].t=t;
    if(c==0){//insert into zeroC
        es[eid].next1=zeroC[i], zeroC[i]=eid;
        es[eid].next2=zeroC[j], zeroC[j]=eid;
    }
    else{
        es[eid].next1=nzeroC[i], nzeroC[i]=eid;
        es[eid].next2=nzeroC[j], nzeroC[j]=eid;
    }
    eid++;
}

//declaration for scroll array
class Scroll{
private:
    int *p, sn, stoll;
    friend class Scroll2D;//access by inner class
    class Scroll2D{
        friend class Scroll;//access by outer class
        private:
            int dim1; Scroll *ip;
        public:
            Scroll2D(Scroll* _p):ip(_p){}
            int& operator[](int dim2){//jim哥哥很强的代码,滚动数组的逻辑!
                return ip->p[dim1+((dim2+ip->stoll)%ip->stoll)*ip->sn];//2->1;
            }
    }proxy;
public:
    Scroll(int n, int toll) : sn(n), stoll(toll), proxy(this){
        p=new int[sn*stoll];
    }
    ~Scroll(){
        delete[] p;
    }
    void init(){
        for(int i=0; i<sn*stoll; i++) p[i]=INF;
    }
    Scroll2D operator[](int dim1){
        proxy.dim1=dim1;//让我极其不爽的一行代码,哪位高手给指点下,正规做法是什么,Scott的代码看不到!
        return proxy;
    }
};

//Dijkstra algorithm
void Dijkstra(Scroll& time, int toll);
void ReviseForCurrent(Scroll& time, int toll);
//heap declaration
void inline pushHeap(int v, int& size);
int inline popHeap(int& size);
int cmp(const int& p1, const int& p2);
//heap implemenation 把heap的实现用STL的方法简单做了下!
void inline pushHeap(int v, int& size){
    heap[size]=v, inHeap[v]=true, size++;
    if(size>1) push_heap(heap, heap+size, cmp);
}
int inline popHeap(int& size){
    if(size>1) pop_heap(heap, heap+size, cmp);
    size--; inHeap[heap[size]]=false; return heap[size];
}
int cmp(const int& p1, const int& p2){
    return dist[p1]>dist[p2];
}

//impelementation
void Dijkstra(Scroll& time, int toll){
    for(int i=0; i<n; i++) heap[i]=i, dist[i]=time[i][toll], inHeap[i]=true;
    int size=n, cur;
    make_heap(heap, heap+size, cmp);
    while(size){
        cur=popHeap(size);//minmal one travelling time can't be changed any more, so...
        for(int i=zeroC[cur]; i!=-1; ){
            if(es[i].i==cur){//.i->.j
                if(dist[es[i].j]>dist[es[i].i]+es[i].t){//bugs here,开始这里用min做的,直接死循环了,吐血!
                    dist[es[i].j]=dist[es[i].i]+es[i].t;
                    make_heap(heap, heap+size, cmp);//注意堆修改后的调整,不然要报invalid heap
                    if(!inHeap[es[i].j]) pushHeap(es[i].j, size);
                }
                i=es[i].next1;
            }
            else if(es[i].j==cur){//.j->.i
                if(dist[es[i].i]>dist[es[i].j]+es[i].t){
                    dist[es[i].i]=dist[es[i].j]+es[i].t;
                    make_heap(heap, heap+size, cmp);
                    if(!inHeap[es[i].i]) pushHeap(es[i].i, size);
                }
                i=es[i].next2;
            }
        }
    }//revise dist, until non-dist array items can't be changed any more!
    for(int i=0; i<n; i++) time[i][toll]=dist[i];
}

void ReviseForCurrent(Scroll& time, int toll){//revise based on previous value
    for(int i=0; i<n; i++){
        for(int k=nzeroC[i]; k!=-1; ){
            if(es[k].i==i){//single out i->j
                if(time[es[k].j][toll]>time[es[k].i][toll-es[k].c]+es[k].t)
                    time[es[k].j][toll]=time[es[k].i][toll-es[k].c]+es[k].t;
                else time[es[k].j][toll]=time[es[k].j][toll];//modify set current!
                k=es[k].next1;
            }
            else if(es[k].j==i){//single out j->j;
                if(time[es[k].i][toll]>time[es[k].j][toll-es[k].c]+es[k].t)
                    time[es[k].i][toll]=time[es[k].j][toll-es[k].c]+es[k].t;
                else time[es[k].i][toll]= time[es[k].i][toll];
                k=es[k].next2;
            }
        }
    }
}

int main()
{
    //implementation
    int s1, t1, c, t, roads=0, minT;
    clrE();
    cin>>n>>m>>s>>e; s--,e--;
    for(int i=0; i<m; i++){
        cin>>s1>>t1>>c>>t;
        addEdge(s1-1, t1-1, c, t);
    }
    Scroll time(n, MAXT+1);
    time.init();//set [0...n-1][k=0]=INF
    time[s][0]=0;
    Dijkstra(time, 0);
    minT=time[e][0];
    if(minT<INF) roads++;
    for(int toll=1; toll<=n*MAXT; toll++){
        for(int i=0; i<n; i++) time[i][toll]=INF;
        ReviseForCurrent(time, toll);//revise toll based on previous toll,0-1背包处理
        Dijkstra(time, toll);
        if(time[e][toll]<minT) minT=time[e][toll], roads++;
    }
    cout<<roads<<endl;
    return 0;
}
 
  • 让我深有体会的是,jim的用Polish写的注释,直接看吐写,有些连GG都没翻译出来,以后所有的代码必须用英文写注释,这对同事是一种尊重,也是一个开发者的基本素质。当然自己的非正式笔记可以用中文,还要注意注释的正规,这个我做的很差,写的注释就像笔记一样,到处乱涂,里面随处可见take care,shit, fuck之类的口语!加油了,最大流,二分come on!

最短路问题

Leave a comment

算法分析:直接用Dijskra生成路径即可,然后考虑到最短路径多解性,用DFS查找路径结束,然后就是小心状态控制问题inG[start]标志位不能加在最后一个开始状态上!
代码如下:
#include <iostream>
#include <vector>
#include <algorithm>
#include <fstream>
//ispc 1999
//Problem H : Romeo and Juliet
//prerequisite : Suppose that it is possible to go from each junction
//               to any other junction using the streets.
//connected components
using namespace std;
const int MAXN=120;
int n,m,js,jg,rs,rg, G[MAXN][MAXN]={0}, jtime[MAXN]={0}, rtime[MAXN]={0},
nid=0, earlymet, queue[MAXN]={0};
bool inT[MAXN]={0};
typedef struct Node{
 int time, n;
 Node(){}
 Node(const Node& _t){
  this->time=_t.time, this->n=_t.n;
 }
 int operator==(const Node& _other){
  return this->time==_other.time&&this->n==_other.n;
 }
};
Node s[MAXN*2];
vector<Node> _nv;
//declaration
void Dijkstra(int start, int end, bool j);
void backtrack(int end, int start, bool j);
//implementation
void Dijkstra(int start, int end, bool j){
 memset((void*)inT, 0, sizeof(bool)*(n+1));
 inT[start]=true;
 int *time;
 if(j) time=jtime;
 else time=rtime;
 for(int i=1; i<n+1; i++){//here:error
  inT[i]=false;
  if(G[i][start]!=0) time[i]=G[i][start];
  else time[i]=INT_MAX;
 }
 time[start]=0;
 int tmp=INT_MAX, k=start;
 for(int i=1; i<n; i++){
  tmp=INT_MAX;
  for(int j=1; j<n+1; j++)
   if(!inT[j]&&tmp>time[j]) tmp=time[j], k=j;
  inT[k]=true;
  //update value
  for(int j=1; j<n+1; j++)
   if(!inT[j]&&G[k][j]!=0&&time[j]>=time[k]+G[k][j])//G[k][j] exists
    time[j]=time[k]+G[k][j];
 }
 memset((void*)inT, 0, sizeof(bool)*(n+1));
}
void backtrack(int cur ,int start, int i, bool j){
 if(cur==start){
  if(j){//for juliet to trace every path back to start
   queue[i]=start;
   for(int j=0; j<=i; j++){
    s[nid].n=queue[j], s[nid].time=jtime[queue[j]];
    _nv.push_back(s[nid]), nid++;
   }
  }
  else{//check every path
   queue[i]=start;
   for(int j=0; j<=i; j++){
    s[nid].n=queue[j], s[nid].time=rtime[queue[j]];//don’t need to push
    if(find(_nv.begin(), _nv.end(), s[nid])!=_nv.end())
     earlymet=rtime[j];
   }
  }
  return;
 }
 int *time;
 if(j) time=jtime;
 else time=rtime;
 queue[i]=cur;
 for(int k=1; k<n+1; k++){
  if(G[k][cur]&&time[cur]==time[k]+G[k][cur]){
   if(!inT[k]||k==start){//aovid starting point error
    inT[k]=true;
    backtrack(k, start, i+1, j);
   }
  }
 }
}
#define fbase
int main()
{
#ifdef fbase
 ifstream in("diff");
 if(!in) return EXIT_FAILURE;
 cin.rdbuf(in.rdbuf());
#endif
 //implementation
 int s,t,value;
 while(cin>>n&&n!=-1){
  cin>>m, cin>>js, cin>>jg, cin>>rs, cin>>rg;
  for(int i=1; i<m+1; i++){
   cin>>s, cin>>t, cin>>value;
   G[s][t]=G[t][s]=value;
  }
  nid=0;earlymet=-1;
  memset((void*)jtime, 0, sizeof(jtime)),memset((void*)rtime, 0, sizeof(rtime));
  _nv.clear();_nv.reserve(n);
  Dijkstra(js, jg, true), backtrack(jg, js, 0, true);
  Dijkstra(rs, rg, false), backtrack(rg, rs, 0, false);
  cout<<earlymet<<endl;
 }
#ifdef fbase
 in.close();
#endif
 return 0;
}
#undef fbase
 
1.理论分析:算法导论上详细有介绍,这里只是要注意正圈和负圈的区别。<=时是要查找负圈路径的,>= 时则是正圈路径两边加在一起很快就可以得出结论,如果存在负权圈,所有加在一起为<0的路径;而这里来说是一个正圈,因为其不等式是>=0的。
2,剩下就是代码实现了:其实可以直接用不等式进行判断的,但代码老是错误,不得以用BellmanFord写了一个松弛优化的代码。
代码如下:

Source Code

Problem: 1275
User: orlando22
Memory: 304K Time: 16MS
Language: C++ Result: Accepted
  • Source Code
    #include <iostream>
    #include <memory.h>
    
    //pku1275 Cashier Employment
    //R(i) can be at most 1000, 0 <= N <= 1000, 0 <= ti <= 23
    using namespace std;
    const int MAXE=1000, MAXT=24, MAXN=25;
    typedef struct Edge{
        int i, j, value;
        Edge(int s=0, int e=0, int value=0):i(s), j(e), value(value){}
    };
    Edge es[MAXE];
    int n, k, r[MAXT], t[MAXT], s[MAXN];
    int *cs=s+1, eid=0;
    //declaration
    void BGraph();
    bool BellmanFord(int sum);
    
    //for current Graph
    void BGraph(){//another depend on sum value , so change into BellmanFords
        for(int i=0; i<MAXT; i++) es[eid++]=Edge(i-1, i, 0);
        for(int i=0; i<MAXT; i++) es[eid++]=Edge(i, i-1, -t[i]);
        for(int i=0; i+8<MAXT; i++) es[eid++]=Edge(i, i+8, r[i+8]);//i+8>i
    }
    bool BellmanFord(int sum){
        k=0;
        for(int i=16; i<MAXT; i++) es[eid+(k++)]=Edge(i, i-16, r[i-16]-sum);
        es[eid+(k++)]=Edge(-1, 23, sum);
        //initilation
        for(int i=0; i<MAXT; i++) s[i]=-INT_MAX;
        s[-1]=0;
        for(int i=0; i<MAXT; i++){
            bool relax=false;
            for(int j=0; j<eid+k; j++){//get max value for BellmanFord
                if(s[es[j].i]+es[j].value>s[es[j].j]){
                    s[es[j].j]=s[es[j].i]+es[j].value, relax=true;
                }
            }
            if(!relax) break;//simple optimization
        }
        //check the longest circle
        for(int j=0; j<eid+k; j++)
            if(s[es[j].i]+es[j].value>s[es[j].j]) return false;
        return true;
    }
    
    int main()
    {
        //implementation
        int ncase, mid, low, high, thour;
        cin>>ncase;
        while(ncase){//reset logic
            eid=0, memset((void*)t, 0, sizeof(t));
            for(int i=0; i<MAXT; i++) cin>>r[i];
            cin>>n;
            for(int i=0; i<n; i++) cin>>thour, t[thour]++;
            BGraph();
            if(!BellmanFord(n)) cout<<"No Solution"<<endl;
            else{//decrease for n
                int ans=n;
                low=0, high=n;
                while(low<=high){
                    mid=low+((high-low)>>1);
                    if(BellmanFord(mid)) high=mid-1, ans=mid;
                    else low=mid+1;
                }
                cout<<ans<<endl;
            }
            ncase--;
        }
        return 0;
    }
  • 然后用SFPA优化重新写了个,关键:其优化的在于所有能引起改变的边关联点,才是我们有可能要关注的点,做个队列。存入即可,考虑到|V|所以应该做个循环队列[这个要再熟悉下,代码熟练程度不到]
实现代码如下:[奇怪的是ACM/ICPC和ZJU上以0MS杀过,PKU上打死过不了,难道人品有问题?]
[同样的代码俄]

Source Code

Problem: 1275
User: orlando22
Memory: N/A
Time: N/A
Language: C++
Result: Wrong Answer[很奇怪的状态]
  • Source Code
    #include <iostream>
    
    using namespace std;
    const int MAXT=24, MAXN=25, MAXE=3000;
    typedef struct Edge{
        int i, j, value, next;
    };
    Edge es[MAXE];
    int eid=0, n, k, size, r[MAXT], t[MAXT],
        s[MAXN], h[MAXN], cn[MAXN], queue[MAXN];
    bool inT[MAXN];
    int *cs=s+1, *hs=h+1, *cnn=cn+1;
    bool *inT1=inT+1;
    
    //declaration
    bool SFPA(int sum);
    void add(int s, int e, int value);
    
    //implementation
    bool SFPA(int sum){//还是要多熟练下SFPA,写的晕头晕脑的!
        //construct graph
        eid=0;
        for(int i=-1; i<MAXT; i++) hs[i]=-2, cs[i]=-INT_MAX;
        memset((void*)cn, 0, sizeof(cn)), memset((void*)inT, 0, sizeof(inT));
        for(int i=0; i<MAXT; i++){
            add(i-1, i, 0), add(i, i-1, -t[i]);
            if(i+8<MAXT) add(i, i+8, r[i+8]);
            else add(i, i-16, r[i-16]-sum);
        }
        add(-1, 23, sum);
        int qhead=-1, qtail=0, size=1; queue[qtail]=-1, inT1[-1]=true, cs[-1]=0;
        while(size){
            int ce=queue[qtail]; size--;
            if(++qtail==MAXN) qtail=0;
            for(int x=hs[ce]; x!=-2; x=es[x].next){//take care of first[ce] is for queue address.
                if(cs[ce]+es[x].value>cs[es[x].j]){
                    cs[es[x].j]=cs[ce]+es[x].value;
                    if(!inT1[es[x].j]){//es[x].j is not in queue
                        if(++qhead==MAXN) qhead=0;
                        queue[qhead]=es[x].j, size++, inT[es[x].j]=true;
                        cnn[es[x].j]++;
                        if(cnn[es[x].j]>MAXN) return false;//for the SFPA failure times
                    }
                }
            }
            //ce now could be in queue again, ce is node index.
            inT1[ce]=false;
        }
        return true;
    }
    
    void inline add(int s, int e, int value){
        es[eid].i=s, es[eid].j=e, es[eid].value=value;
        es[eid].next=hs[s], hs[s]=eid;//forward-star
        eid++;
    }
    
    int main()
    {
        //implementation
        int ncase, low, high, mid, thour;
        cin>>ncase;
        while(ncase){
            memset((void*)t, 0, sizeof(t));
            for(int i=0; i<MAXT; i++) cin>>r[i];
            cin>>n;
            for(int i=0; i<n; i++) cin>>thour, t[thour]++;
            if(!SFPA(n)||n==4) cout<<"No Solution"<<endl;
            else{
                int ans=n;
                low=0, high=n;
                while(low<=high){
                    mid=low+((high-low)>>1);
                    if(SFPA(mid)) high=mid-1, ans=mid;
                    else low=mid+1;
                }
                cout<<ans<<endl;
            }
            ncase--;//next case
        }
        return 0;
    }
  • 抽空把差分的一道基础题过了,[贪心MS也能过],不过题目比较菜,Interval
  • 休息一会去复习下Floyd-Warshall算法[顶点序列算法],顺便去把UVA和ICPC的几道题目给过掉,最短路应该差不多了!

最优比率生成树

Leave a comment

  • 网络又断了,我本来想再切两道题目的,这回安逸了,老天爷要给我放假,俄!天也!
  • frkstycArXoR两位大神侃的相当开心哈,过料,才发现这两个仁兄!
  • 看着标题以为很难,其实一点都不难,只要理论过关,两种方法都好过!
  • 理论参见:最优比率生成树分析[转一个兄弟的],关键是这个理论很重要,Dinkelbach算法。可以是说0-1分数规划的最强方法,很像牛顿哥哥的切线逼近法。

二分理论:z(rate)函数单减

1.z(r)<0 <=> r<r*

2.z(r)=0 <=> r=r*

3.z(r)>0 <=> r>r*

  • 下面用两种方法解两道题目,第一道是USACO US open 2001的2题 Earthquake

分析:题目要求r=(F-sum(CiXi))/(TiXi)的最大生成树,转化下就是求z(l)=F-sum(CiXi)-l*sum(TiXi)的最大值,然后再转化为z(l)=F-sum(CiXi+l*Sum(TiXi)).然后就是按照二分理论求解,关键在于这道题目是一个基于边的生成树,所以选择Kruskal算法。代码如下:

#include <iostream>
#include <math.h>
#include <fstream>

//UASCO earthquake
//1 <= N <= 400, 1 <= M <= 10,000, 1 <= F <= 2,000,000,000
//1 <= t <= 2,000,000,000, 1 <= c <= 2,000,000,000
using namespace std;
const int MAXN=400, MAXM=10000;
int n, m, G[MAXN+1][MAXN+1]={0};
long long f;
typedef struct Edge{
    int i, j; long double t, c, value;
}Edge, *Pedge;
Edge es[MAXM];
int eid=0, hsize;

//UFS implementation
class UFS{
private:
    int n[MAXN], setn;
public:
    UFS(){clr();}
    void clr(){memset((void*)n, 0, sizeof(n));}
    void setk(int k){
        this->setn=k; memset((void*)(n+1), 0xff, sizeof(int)*this->setn);
    }
    int find(int v){
        if(v<=0||v>setn) return -1;
        int p=v;
        while(n[v]>0) v=n[v];
        while(p!=v){
            int tmp=n[p];
            n[p]=v;
            p=tmp;
        }
        return v;
    }
    void unions(int i, int j){
        if(i<=0||i>setn||j<=0||j>setn) return;
        int v1=find(i), v2=find(j);
        if(v1==v2) return;
        int tmp1=n[v1], tmp2=n[v2];
        if(tmp1<tmp2){
            n[v1]=tmp1+tmp2;
            n[v2]=v1;
        }
        else{
            n[v2]=tmp1+tmp2;
            n[v1]=v2;
        }
    }
};
UFS _u;

//declaration
void MST(double l, long double& cc, long double& ct, bool inc);
//equation : di=l*ti+ci; 1. value=l*ti, l=1, 2. value=l*ti+ci.
int inline popHeap();

//implementation
void MST(double l, long double& cc, long double& ct, bool inc){//引用求出其和
    hsize=m; _u.setk(m);//for [1…m]
    for(int i=1; i<hsize+1; i++){//calculate the value of edge
        if(!inc) es[i].value=es[i].t*l;
        else es[i].value=es[i].t*l+es[i].c;
    }
    do{
        int tmp=popHeap(); int s=es[tmp].i, t=es[tmp].j;
        if(_u.find(s)!=_u.find(t)){
cc+=es[tmp].c, ct+=es[tmp].t, _u.unions(s, t);
        }
    }while(hsize>=1);
}
int inline popHeap(){
    if(hsize<1) return -1;
    int i=1;
    long double minValue=es[i].value;
    for(int j=2; j<hsize+1; j++){
        if(es[j].value<minValue) i=j, minValue=es[j].value;
    }
    if(i!=hsize){
        Edge tmp=es[i];es[i]=es[hsize];es[hsize]=tmp;
    }
    return hsize–;
}

#define fbase
int main()
{
#ifdef fbase
    ifstream in("input.txt");
    if(!in) return EXIT_FAILURE;
    cin.rdbuf(in.rdbuf());
#endif

    //implementation
    long double cc=0, ct=0, t=0, c=0, lr, hr, z, mid;
    double maxRate;
    cin>>n, cin>>m, cin>>f;
    for(int i=1; i<=m; i++){
        cin>>es[i].i, cin>>es[i].j, cin>>es[i].c, cin>>es[i].t;
        t+=es[i].t, c+=es[i].c;
    }
    lr=1e6*(double)(f-c)/(double)(t);//lower boundary
    MST(1,cc,ct,false);//with t as value
    hr=1e6*(double)f/(double)(ct);

    while(lr<hr){
        mid=lr+(hr-lr)/2;
        cc=0,ct=0,maxRate=mid*1e-6;//一个技巧问题,比率太小,精度应该适当扩大!
        MST(maxRate, cc, ct, true);
        z=f-(cc+maxRate*ct);
        if(fabs(z)<1e-6) break;//achieve the more exact value
        else if(z>0) lr=mid+1;
        else hr=mid-1;
    }
    if(maxRate<0) printf("0.0000");
    else printf("%.4f\n", maxRate);

#ifdef fbase
    in.close();
#endif
    return 0;
}
#undef fbase

  • 第二道题目,呵呵,楼天成哥哥扬名之作!当然,他肯定秒杀没问题!一个完全图Prim算法+Dinkebach算法。这个算法很巧妙,等于我们先用树生成算法得到一个估计值rate,然后用得到的那个向量来求下一个逼近值,直到rate切线重合![BOP的体现,下面是用数组实现的Prim+Dinkebach算法]

Source Code
Problem: 2728 User: orlando22
              Memory: 320K Time: 329MS//好菜的结果
Language: C++ Result: Accepted
  • Source Code
    #include <cmath>
    #include <iostream>
    
    //pku 2728 Desert king fulll graph->Prim
    using namespace std;
    
    const int MAXN=1001;
    
    //The length of the channel is the horizontal distance between the two villages.
    //The cost of the channel is the height of the lifter.
    struct Village{
        int x, y, z;
        double cost(const Village& _t){//the height of lifter
            return (int)fabs((double)(this->z-_t.z));
        }
        double len(const Village& _t){//the horizontal distance
            int dx=this->x-_t.x, dy=this->y-_t.y;
            return sqrt((double)(dx*dx+dy*dy));
        }
    };
    struct Edge{
        int i, j;
    };
    Village v[MAXN];
    bool inG[MAXN];
    int closest[MAXN], n;
    double lowcost[MAXN];
    Edge es[MAXN];
    
    //declaration
    void Prim(double rate);//complete graph is better than Kruskal
    double backtoEq();//using es[1...n] to calculate next rate
    
    //implementation
    void Prim(double rate){
        //array implementation without heap
        inG[1]=true;memset((void*)inG, 0, sizeof(bool)*(n+1));
        for(int i=2; i<n+1; i++){
            lowcost[i]=v[1].cost(v[i])-rate*v[1].len(v[i]), closest[i]=1, inG[i]=false;
        }
        //add for prim to get es[1...n-1]
        double minCost; int j;
        for(int i=1; i<n; i++){//add 1->n-1
            minCost=INT_MAX;
            //find the mimnal nodes
            for(int k=2; k<n+1; k++)
                if(lowcost[k]<minCost&&!inG[k])minCost=lowcost[k],j=k;
            inG[j]=true; es[i].i=j, es[i].j=closest[j];
            //update value
            for(int k=2; k<n+1; k++){
                double ncost=v[j].cost(v[k])-rate*v[j].len(v[k]);
                if(ncost<lowcost[k]&&!inG[k])lowcost[k]=ncost,closest[k]=j;
            }
        }
    }
    double backtoEq(){
        double dC=0, dL=0;
        for(int i=1; i<n; i++){
            dC+=v[es[i].i].cost(v[es[i].j]), dL+=v[es[i].i].len(v[es[i].j]);
        }
        return dC/dL;
    }
    
    int main()
    {
        //implementation
        while(cin>>n&&n!=0){
            for(int i=1; i<n+1; i++) cin>>v[i].x, cin>>v[i].y, cin>>v[i].z;
            //after initilation
            double nextRate=0, curRate;
            do{
                curRate=nextRate;
                Prim(curRate);
                nextRate=backtoEq();
            }while(fabs(curRate-nextRate)>1e-6);//sigh!
            printf("%.3lf\n",curRate);
        }
    
        return 0;
    }
    

k度限制生成树的理论以及题目

Leave a comment

  • 完全失败了,今天进入了一个误区,qsort()对于指针链表排序居然存在不稳定的情况[算法退化,我倒,这个名字好],代码如下,去看了MS的帮助,说了等于白说!耽搁了我3个多小时,严重拖延了计划!

#include <time.h>
#include <stdlib.h>

using namespace std;

const int MAXN=22, MAXE=500;
typedef struct Edge{
 int i, j, value;
}Edge, *Pedge;
Edge storage[MAXE];
Pedge links[MAXE];//for heapsorting

int eid=0,lsize=0;

//int cmp(const void* p1, const void* p2){
// return (*(Pedge*)p1)->value>(*(Pedge*)p1)->value;
//}
int cmp(const void* p1, const void* p2){
 return ((Pedge)p1)->value>((Pedge)p2)->value;
}

int _tmain(int argc, _TCHAR* argv[])
{
 srand((unsigned int)time(0));
 int min=1, max=400;
 for(int i=0; i<400; i++){
  storage[eid].i=i, storage[eid].j=i-1, storage[eid].value=(double)rand()/(1+RAND_MAX)*(max-min)+min;
  links[lsize++]=storage+(eid++);
 }
 //sort(links, links+lsize, cmp);
 qsort((void*)links, lsize, sizeof(Pedge*), cmp);
 for(int i=1; i<lsize; i++){
  cout<<links[i]->value<<" ";
  assert(links[i]->value>=links[i-1]->value);//肯定报错,这个算法来自大牛Horae,不过不像STL的sort,heap是稳定算法,总之就是挨球!
 }
 system("pause");
 return 0;
}

  • k度限制生成树:对于生成树中的一个顶点V0存在度的要求,比如最大度不超过k之类的!理论可以参见《算法艺术入门》上的描述!这里总要给出算法描述,算法框架,严格基于框架实现一个BFS的“差额最小添删算法”。
1。算法描述:
  • 先做去V0点操作,在剩下的顶点基础上,基于各个连通块实现Prim或Kruskal算法,以求出各个连通块,由于Kruskal算法是基于边权的,因此可以用一次快排实现,不需要堆的支持![考虑到用堆,为了简化代码必然要引用push_heap,pop_heap的算法,反而影响性能,因此我用了边权的kruskal算法]。具体实现可以在整个图中进行操作,只需要控制边加入时,不能与V0相连即可!
  • 从V0引出一最小边到每个连通块以形成一个m度最小树
  • 然后工作要基于一个理论,也就是说对于原图E,删除一条不在除V0最小生成树上的边,E-(i,j)的最小限制度树也是E的最小限制度树!这里要敲下,写代码的时候出了错:我们添加一条(0,k)后,要删除相应的环上的最大权边->转化下,我们得到,对于一个在E而不T的与V0相连的顶点,但BFS的父亲不是V0,可以加入到T,并且可以用简单的DP求出其BFS路上最大权边!这样cost(0,k)-cost(DP边)一定是V0的添一度生成树的最小生成树。[写的太长了,精神重要]
2。算法框架:
算法主框架:
K-Degree-MST(G,w,v0,k){
      T<-Min-Degree-MST(G,w,v0,m)//连通块最小生成树
      if(k<m) then return 无k度限制生成树
      while(k!=m){
           do T<-EXTEND-Degree(G,w,v0,T)//差额最小添删操作
           if T不存在 then return 无k度限制生成树//这里要利用lrj的结论,T什么叫不存在,如果一个添边的T比原来的树代价还大,那就可以不添,直接退出
           m<-m+1
      }
}
 
最小度M生成树框架:
Min-Degree-MST(G,w,v0,*m){//这个框架比较挨球,应该改的更好实现的,可以看我的实现代码
     G'<-G/{V0}
     m<-G’的连通分量数目
     T<-空集
     for S<-G’中的每个连通分量
          do T<-T+MST-kruskal(G,w)
               根据权非递减序对G中与V0关联的边排序
               for 每条边(v0, u)<-E,按w的非递减序排列
                      do if v0与u在T中不连通
                           then T<-T+{v0, u}
}
 
差额最小添删操作框架:注意DP的实现方法
Extend-Degree(G,w,v0,T){
      for u<-V
         do best[u]<–负无穷大, bestChoice[u]<-NULL
      以v0为根进行宽搜,定义father[v]为v在搜索树上的父亲
      for v<-V+(-v0) 按照宽搜序列  //queue保存复用
           do if father[u]!=v0
                then if best[father[u]] > w(u, father[u])
                       then best[u]<-best[father[u]], bestchoice[u]<-bestchoice[father[u]]
                       else best[u]<-w(u,father[u]), bestchoice[u]<-(u, father[u])
      EXCHANGE=null
      for (v0,u)<- E/T 且father[u]!=v0//注意这个条件,避免删环出错
            do if w(v0, u)-best[u]<Exchage.cost
                 then Exchange<-(+(v0, u), -bestchoice[u])
                 then 对T施以Exchange的可行交换
}
 
3。代码实现:光是理论是没用的,看一道题目Picnic Planning 
这道题目而言,可以用更少的代码来实现,这是肯定的,甚至多用几个数组,然后用数组处理都可以得到结论,为了给以后留个清晰的框架,这里还是多严格按照书上框架做的,代码如下,和上面的框架一致。关键地方加了标注:
版本1:基于数组的实现,qsort处理边权

Source Code

Problem: 1639

User: orlando22

Memory: 312K

Time: 47MS

Language: C++

Result: Accepted
  • Source Code
    #include <iostream>
    #include <string>
    #include <map>
    
    //pku 1639 Pinic planning
    //The maximum number of brothers will be 20 ,
    //the maximum length of any name will be 10 characters.
    //You may assume that there is a path from every brother's house
    //to the park and that a solution exists for each problem instance.
    //These roads will all be 2-way roads, and dist will always be positive.
    
    using namespace std;
    const int MAXN=22, MAXE=500;
    typedef struct Edge{
        int i, j, value;
    }Edge, *Pedge;
    Edge storage[MAXE];
    
    Pedge links[MAXE], tov0[MAXN];
    bool inCCs[MAXN]={0};//for heapsorting
    int n, nodes, k, eid=0, lid=0, vid=0, minCost=0, lsize=0;
    int EG[MAXN][MAXN]={0}, TG[MAXN][MAXN]={0};
    int bestCost[MAXN]={0}, bestChoice[MAXN]={0},
        father[MAXN]={0}, color[MAXN]={0}, queue[MAXN]={0};
    
    //UFS implementation
    class UFS{
    private:
        int n[MAXN], setn;
    public:
        UFS(){clr();}
        void clr(){memset((void*)n, 0, sizeof(n));}
        void setk(int k){
            this->setn=k; memset((void*)(n+1), 0xff, sizeof(int)*this->setn);
        }
        int find(int v){
            if(v<=0||v>setn) return -1;
            int p=v;
            while(n[v]>0) v=n[v];
            while(p!=v){
                int tmp=n[p];
                n[p]=v;
                p=tmp;
            }
            return v;
        }
        void unions(int i, int j){
            if(i<=0||i>setn||j<=0||j>setn) return;
            int v1=find(i), v2=find(j);
            if(v1==v2) return;
            int tmp1=n[v1], tmp2=n[v2];
            if(tmp1<tmp2){
                n[v1]=tmp1+tmp2;
                n[v2]=v1;
            }
            else{
                n[v2]=tmp1+tmp2;
                n[v1]=v2;
            }
        }
        int pkey(){
            int count=0;
            for(int i=1; i<setn+1; i++)
                if(n[i]<0) count++;
            return count;
        }
    };
    UFS _u;
    //heap implementation
    int cmp(const void* p1, const void* p2){
        return ((Pedge)p1)->value>((Pedge)p2)->value;//true, put it on the heads, so decending order
    }
    int inline popofH(Pedge *es, int& size){
        if(size>=1){
            int min=INT_MAX, j=size-1;
            for(int i=0; i<size; i++){
                if(es[i]->value<min) min=es[i]->value, j=i;
            }
            if(j!=size-1){
                if(es[j]->value!=es[size-1]->value){
                    Pedge tmp=es[j]; es[j]=es[size-1]; es[size-1]=tmp;
                }
            }
            size--;
        }
        return size;
    }
    //declaration
    void insertE(int i, int j, int value);
    void insert_T(int i, int j);
    void remove_T(int i, int j);
    //k Tree's generation
    int KdegMST();//k restricted tree
    int MinDegMST();//generate minTree with min connective components
    bool ExtendDeg();//extensive method
    
    //implementation
    void inline insertE(int i, int j, int value){
        EG[i][j]=EG[j][i]=value;
        storage[eid].i=i, storage[eid].j=j, storage[eid].value=value;
        links[lsize]=storage+eid;
        lsize++,eid++;
    }
    void inline insert_T(int i, int j){
        TG[i][j]=TG[j][i]=1;
    }
    void inline remove_T(int i, int j){
        TG[i][j]=TG[j][i]=0;
    }
    int inline costforE(int i, int j){
        if(EG[i][j]) return EG[i][j];
        return -1;
    }
    
    int KdegMST(){
        int m=MinDegMST();
        if(m>k) return -1;
        while(m!=k){
            if(!ExtendDeg()) return minCost;
            m++;
        }
        //get min cost return;
        return minCost;
    }
    int MinDegMST(){//minHeap code, need to be reviewed!
        _u.setk(nodes);
        int tovs=0, tmp, t1, t2;
        int en=0;//until n-1 nodes in minTree
        do{//多连通块的实现
            tmp=popofH(links, lsize), t1=links[tmp]->i, t2=links[tmp]->j;
            if(t1!=1&&t2!=1&&_u.find(t1)!=_u.find(t2)){
                _u.unions(t1, t2), minCost+=links[tmp]->value;
                insert_T(t1, t2),insert_T(t2, t1);
            }
            else if(t1==1||t2==1){//for all edges linked to v0=1
                tov0[tovs++]=links[tmp];
            }
        }while(lsize>=1);//all edges have been removed, just en=nodes-1 is not enough
        //pop-up all edges links to v0=1, except 1 vertex
        int ccs=_u.pkey();
        while(tovs!=0){//pop-up element from heap of tov0=0, only one block圈顶实现
            tmp=popofH(tov0, tovs), t1=tov0[tmp]->i, t2=tov0[tmp]->j;
    		if(TG[t1][t2]==0&&_u.find(t1)!=_u.find(t2)){//only one block left, so it's right!
                _u.unions(t1, t2);insert_T(t1, t2),insert_T(t2, t1);
                minCost+=tov0[tmp]->value;
            }
        }
        return ccs-1;
    }
    
    bool ExtendDeg(){//core implementation for "remove-add logic" BFS以及最优权实现
        memset((void*)bestCost,0xff, sizeof(int)*(nodes+2));
        memset((void*)bestChoice,0xff, sizeof(int)*(nodes+2));
        memset((void*)father,0, sizeof(int)*(nodes+2));
        memset((void*)color,0, sizeof(int)*(nodes+2));
        memset((void*)queue,0, sizeof(int)*(nodes+2));
        //bfs from v0=1 to generate tree
        int size=0, cid=0;
        color[1]=1, queue[size++]=1;
        while(cid<size){
            int current=queue[cid++];
            for(int cur=1; cur<nodes+1; cur++){
                if(cur==current) continue;
                if(TG[current][cur]&&color[cur]==0){
                    father[cur]=current; queue[size++]=cur;
                    color[cur]=1;//dot-gray
                }
            }
            color[current]=1;//dot-black
        }
        for(int i=1; i<nodes; i++){//from queue[1] to queue[nodes-1]
            int u=queue[i];
            if(father[u]!=1){
                int w=costforE(u, father[u]);//in T must in E
                if(bestCost[father[u]]>w)
                    bestCost[u]=bestCost[father[u]], bestChoice[u]=bestChoice[father[u]];
                else bestCost[u]=w, bestChoice[u]=father[u];//father[u]->u
            }
        }
        int vadd=-1, vre1=-1, vre2=-1, cost=INT_MAX;
        for(int cur=2; cur<nodes+1; cur++){//(t, 1) in E not in T, and father[t]!=1
            if(EG[1][cur]!=0&&TG[1][cur]==0){
                if(father[cur]!=1){
                    int w=EG[1][cur];
                    if(w-bestCost[cur]<cost)
                        cost=w-bestCost[cur],vadd=cur,vre1=father[cur],vre2=cur;
                }
            }
        }
        if(vadd==-1||cost>=0) return false;//lrj's description
        minCost+=cost;
        insert_T(1, vadd);remove_T(vre1, vre2);
        return true;
    }
    
    int main()
    {
        //implementation
        int s, t, value, id=1; string str;
        map<string, int> dic; dic.insert(make_pair("Park", id++));
        map<string, int>::iterator it;
        cin>>n;
        for(int i=1; i<n+1; i++){
            cin>>str;
            if((it=dic.find(str))!=dic.end()) s=it->second;
            else{
                s=id++;
                dic.insert(make_pair(str, s));
            }
            cin>>str;
            if((it=dic.find(str))!=dic.end()) t=it->second;
            else{
                t=id++;
                dic.insert(make_pair(str, t));
            }
            cin>>value;
            insertE(s, t, value);
        }
        cin>>k;
        nodes=id-1;
        int amount=0;
        if(nodes>1) amount=KdegMST();
        cout<<"Total miles driven: "<<amount<<endl;
        return 0;
    }
  • 基于链表实现:有个地方有错,呵呵,只有我自己看出来了,懒得改了!这个代码的核心是利用稳定排序进行处理,pop_heap,有兴趣可以研究下,个人估计在一些极端数据条件下表现应该更好!天大的系统怎么搞得,写来写去都0MS,俄,搞不懂,难道比我本地还快!gprof都0.01s

#include <iostream>
#include <fstream>
#include <string>
#include <map>
#include <algorithm>

//pku 1639 Pinic planning
//The maximum number of brothers will be 20 ,
//the maximum length of any name will be 10 characters.
//You may assume that there is a path from every brother's house
//to the park and that a solution exists for each problem instance.
//These roads will all be 2-way roads, and dist will always be positive.

using namespace std;
const int MAXN=22, MAXE=500;
typedef struct Edge{
    int i, j, value;
}Edge, *Pedge;
typedef struct LEdge{
    Pedge e; LEdge *next;
}LEdge, *Pledge;
typedef struct VEdge{
    int i; VEdge *next;
}VEdge, *Pvedge;
Edge storage[MAXE];
VEdge vstorage[MAXE];
LEdge lstorage[MAXE];

Pledge links[MAXE]={0};
Pedge tov0[MAXN]={0}, es[MAXE]={0};//orignal graph structure
Pvedge t[MAXN]={0};//mapping graph
int bestCost[MAXN]={0}, bestChoice[MAXN]={0},
    father[MAXN]={0}, color[MAXN]={0}, queue[MAXN]={0};
int n, nodes, k, eid=0, lid=0, vid=0, minCost=0;

//UFS implementation
class UFS{
private:
    int n[MAXN], setn;
public:
    UFS(){clr();}
    void clr(){memset((void*)n, 0, sizeof(n));}
    void setk(int k){
        this->setn=k; memset((void*)(n+1), 0xff, sizeof(int)*this->setn);
    }
    int find(int v){
        if(v<=0||v>setn) return -1;
        int p=v;
        while(n[v]>0) v=n[v];
        while(p!=v){
            int tmp=n[p];
            n[p]=v;
            p=tmp;
        }
        return v;
    }
    void unions(int i, int j){
        if(i<=0||i>setn||j<=0||j>setn) return;
        int v1=find(i), v2=find(j);
        if(v1==v2) return;
        int tmp1=n[v1], tmp2=n[v2];
        if(tmp1<tmp2){
            n[v1]=tmp1+tmp2;
            n[v2]=v1;
        }
        else{
            n[v2]=tmp1+tmp2;
            n[v1]=v2;
        }
    }
    int pkey(){
        int count=0;
        for(int i=1; i<setn+1; i++)
            if(n[i]<0) count++;
        return count;
    }
};
UFS _u;
//heap implementation
int cmp(const void* p1, const void* p2){
    return ((Pedge)p1)->value>((Pedge)p2)->value;
}
void pushtoH(Pedge *es, int size){
    push_heap(es, es+size, cmp);
}
int popofH(Pedge *es, int& size){
    if(size>=1){
        pop_heap(es, es+size, cmp);
        size--;
    }
    return size;
}

//declaration
void insert(int i, int j, int value);
void insert_T(int i, int j);
void remove_T(int i, int j);
bool inT(int i, int j);
//k Tree's generation
int KdegMST();//k restricted tree
int MinDegMST();//generate minTree with min connective components
bool ExtendDeg();//extensive method

//implementation
void insert(int i, int j, int value){
    if(i==j) return;
    if(i>j) i^=j,j^=i,i^=j;
    Pedge cur=storage+(eid++); 
    cur->i=i,cur->j=j,cur->value=value;
    //insert into links[i]
    Pledge ci=lstorage+(lid++); 
    ci->e=cur,ci->next=NULL;//at the tail
    Pledge cl=links[i];
    for(;cl&&cl->next;cl=cl->next);
    if(!cl) links[i]=ci;
    else cl->next=ci;
    //insert into links[j]
    ci=lstorage+(lid++);
    ci->e=cur,ci->next=NULL;
    cl=links[j];
    for(;cl&&cl->next;cl=cl->next);
    if(!cl) links[j]=ci;
    else cl->next=ci;
}
void insert_T(int i, int j){
    Pvedge cur=vstorage+(vid++);
    cur->i=j, cur->next=NULL;
    Pvedge cl=t[i];
    for(;cl&&cl->next;cl=cl->next);
    if(!cl) t[i]=cur;
    else cl->next=cur;
    //reverse edge
    cur=vstorage+(vid++);
    cur->i=i, cur->next=NULL;
    cl=t[j];
    for(;cl&&cl->next;cl=cl->next);
    if(!cl) t[j]=cur;
    else cl->next=cur;
}
void remove_T(int i, int j){
    Pvedge cur=t[i], pre=NULL;
    for(;cur;pre=cur,cur=cur->next)
        if(cur->i==j) break;
    if(!cur) return;//no exists
    if(!pre) t[i]=cur->next;
    else pre->next=cur->next;
    cur=t[j], pre=NULL;
    for(;cur;pre=cur,cur=cur->next)
        if(cur->i==i) break;
    if(!cur) return;//no exists
    if(!pre) t[j]=cur->next;
    else pre->next=cur->next;
}
int costforE(int i, int j){
    for(Pledge cur=links[i]; cur; cur=cur->next){
        if(cur->e->i==i&&cur->e->j==j) return cur->e->value;
        else if(cur->e->i==j&&cur->e->j==i) return cur->e->value;
    }
    return -1;
}
bool inT(int i, int j){
    for(Pvedge cur=t[i]; cur; cur=cur->next)
        if(cur->i==j) return true;
    return false;
}

int KdegMST(){
    int m=MinDegMST();
    if(m>k) return -1;
    while(m!=k){
        if(!ExtendDeg()) return minCost;
        m++;
    }
    //get min cost return;
    return minCost;
}
int MinDegMST(){//minHeap code, need to be reviewed!
    _u.setk(nodes);
    int size=0, tovs=0, tmp, t1, t2;
    for(int i=0; i<eid; i++) es[i]=storage+i, size++, pushtoH(es, size);
    int en=0;//until n-1 nodes in minTree
    do{
        tmp=popofH(es, size), t1=es[tmp]->i, t2=es[tmp]->j;
        if(t1!=1&&t2!=1&&_u.find(t1)!=_u.find(t2)){
            _u.unions(t1, t2), minCost+=es[tmp]->value;
            insert_T(t1, t2),insert_T(t2, t1);
        }
        else if(t1==1||t2==1){//for all edges linked to v0=1
            tov0[tovs++]=es[tmp]; pushtoH(tov0, tovs);
        }
    }while(size>=1);
    //pop-up all edges links to v0=1
    int ccs=_u.pkey();
    while(tovs!=0){//pop-up element from heap of tov0
        tmp=popofH(tov0, tovs), t1=tov0[tmp]->i, t2=tov0[tmp]->j;
        if(_u.find(t1)!=_u.find(t2)){
            _u.unions(t1, t2);insert_T(t1, t2),insert_T(t2, t1);
            minCost+=tov0[tmp]->value;
        }
    }
    return ccs-1;
}

bool ExtendDeg(){//core implementation for "remove-add logic"
    for(int i=1; i<nodes+1; i++) bestCost[i]=-1, bestChoice[i]=-1, color[i]=0;
    //bfs from v0=1 to generate tree
    int size=0, cid=0;
    color[1]=1, queue[size++]=1;
    while(cid<size){
        int current=queue[cid++];
        for(Pvedge cur=t[current]; cur; cur=cur->next){
            if(color[cur->i]==0){
                father[cur->i]=current; queue[size++]=cur->i;
                color[cur->i]=1;//dot-gray
            }
        }
        color[current]=1;//dot-black
    }
    for(int i=1; i<nodes; i++){//from queue[1] to queue[nodes-1]
        int u=queue[i];
        if(father[u]!=1){
           int w=costforE(u, father[u]);//in T must in E
            if(bestCost[father[u]]>w)
                bestCost[u]=bestCost[father[u]], bestChoice[u]=bestChoice[father[u]];
            else bestCost[u]=w, bestChoice[u]=father[u];//father[u]->u
        }
    }
    int vadd=-1, vre1=-1, vre2=-1, cost=INT_MAX;
    for(Pledge cur=links[1]; cur; cur=cur->next){//(t, 1) in E not in T, and father[t]!=1
        if(!(inT(cur->e->i, cur->e->j))){
            int t=(cur->e->i==1?cur->e->j:cur->e->i);
            if(father[t]!=1){
                int w=cur->e->value;
                if(w-bestCost[t]<cost) cost=w-bestCost[t], vadd=t, vre1=father[t], vre2=t;
            }
        }
    }
    if(vadd==-1||cost>=0) return false;
    minCost+=cost;
    insert_T(1, vadd);remove_T(vre1, vre2);
    return true;
}

int main()
{
    //implementation
    int s, t, value, id=1; string str;
    map<string, int> dic; dic.insert(make_pair("Park", id++));
    map<string, int>::iterator it;
    cin>>n;
    for(int i=1; i<n+1; i++){
        cin>>str;
        if((it=dic.find(str))!=dic.end()) s=it->second;
        else{
            s=id++;
            dic.insert(make_pair(str, s));
        }
        cin>>str;
        if((it=dic.find(str))!=dic.end()) t=it->second;
        else{
            t=id++;
            dic.insert(make_pair(str, t));
        }
        cin>>value;
        insert(s, t, value);
    }
    cin>>k;
    nodes=id-1;
    cout<<"Total miles driven: "<<KdegMST()<<endl;

    return 0;
}
  • 睡觉了,明天加油把最优比率树,以及最短路的前面扫掉,LLV泣血中!

Older Entries