【usaco-Earthquake, 2001 Open】 0-1分数规划 & 最优比率生成树

题意：给定n个点m条边，一开始这些边全都是断的，要修一些边使得n个点全部联通。修完一共可以得到F元，修一条边有成本di和时间ti，要使得 得到的钱数 / 总时间 这个比值最大。



参考资料：

红线内的内容转载自http://www.cnblogs.com/scau20110726/archive/2012/10/19/2730896.html

------------------------------------------------------------------------------------------------------------------

****************************************************

解法之一 0-1分数规划

设x[i]等于1或0, 表示边e[i]是否属于生成树.

则我们所求的比率 r = ∑(benifit[i] * x[i]) / ∑(cost[i] * x[i]), 0≤i<m .

为了使 r 最大, 设计一个子问题---> 让 z = ∑(benifit[i] * x[i]) - l * ∑(cost[i] * x[i]) = ∑(d[i] * x[i]) 最大 (d[i] = benifit[i] - l * cost[i]) , 并记为z(l). 我们可以兴高采烈地把z(l)看做以d为边权的最大生成树的总权值.

然后明确两个性质:

　1. z单调递减

　　证明: 因为cost为正数, 所以z随l的减小而增大.

　2. z( max(r) ) = 0

　　证明: 若z( max(r) ) < 0, ∑(benifit[i] * x[i]) - max(r) * ∑(cost[i] * x[i]) < 0, 可化为 max(r) < max(r). 矛盾;

　　 若z( max(r) ) >= 0, 根据性质1, 当z = 0 时r最大.

到了这个地步, 七窍全已打通, 喜欢二分的上二分, 喜欢Dinkelbach的就Dinkelbach.

****************************************************

已知一个完全图，每条边有两个参数（dis和c），求一棵生成树，使(∑xi×ci)/(∑xi×disi)最小，其中xi当第i条边包含在生成树中时为1，否则为0。
迭代法：
假设rate为当前比率，以ci-rate×disi作为各边的权重，使用Prim算法构造最小生成树，再对该最小生成树求(∑xi×ci)/(∑xi×disi)更新rate，可证明rate可收敛且收敛值即为所求。
二分法：
在一个精度范围内（以1e-6为例），二分查找[0,maxRate]之间的值rate，使(z=∑xi×ci-rate×∑xi×disi)==0，可证明该值即为所求。其中xi为以ci-rate×disi作为各边权重时，使用Prim算法所构造的最小生成树。在二分搜索过程中若z>0，则将区间上移(low=mid+1)，否则将区间下移(high=mid-1)。
使用迭代法，以0作为初始迭代比率：188MS
使用二分法，固定查找范围为[0,31]，精度为1e-6：1422MS(不知二分法大家都有什么优化，分享一下吧^_^)

下面介绍一下该题目的解题思路及相关证明：

1、问题转化：

给定一个rate，z(rate)=∑xi×ci-rate*∑xi×disi，xi为一棵生成树使(∑xi×ci-rate*∑xi×disi)的值最小(下面会介绍求此生成树的方

法)，则rate=(∑xi×ci-z(rate))/( ∑xi×disi)，令rateNex=(∑xi×ci)/( ∑xi×disi)。

若z(rate)>0，则肯定不存在一棵生成树使rate=(∑yi×ci)/( ∑yi×disi)，即rate值无效，且有rateNex >rate；

若z(rate)<0，则我们可得到一个有效比率rateNex，且rateNex<rate；

若z(rate)==0，则rateNex==rate，且rate有效。

因此，如果有且仅有一个rate使z(rate)==0，又由题目所求的最小比率的存在性（有限节点的完全图其生成树只有有限个）可知，该rate值即为所求。

下面证明其存在性和唯一性：

存在性：对于题目所求的最小比率rateMin，由上面的分析必有z(rateMin)=0，又由该最小比率的存在性可知，至少存在一个有效的比率rate使z(rate)=0；

唯一性：只要证明z(rate)函数的单调性即可：

对于两个比率rate1<rate2，

z(rate1)-z(rate2)= (∑xi×ci-rate1*∑xi×disi)-(∑yi×ci-rate2*∑yi×disi)

<=(∑yi×ci-rate1*∑yi×disi)-(∑yi×ci-rate2*∑yi×disi)

=(rate2-rate1)* ∑yi×disi>0

因此，z(rate)为单调递减函数。从而也就证明了满足z(rate)==0的比率rate的唯一性。

由上面的分析，我们就将问题转化为求解比率rate使z(rate)==0。

2、求解：有两种方法对该问题进行求解：迭代法和二分法。

A、迭代法：由上面的分析可知：

当 z(rate)>0时，rate值无效，而rateNex有效且z(rateNex)<=0；

当z(rate)<0时，rateNex<rate，又z(rate)为单调递减函数，故0>=z(rateNex)>z(rate)。

故迭代过程向z(rate)==0收敛，又由有限节点的完全图其生成树只有有限个可知迭代次数必为有限值。

B、二分法：在定出一个搜索范围和精度之后(以[0,MAXRATE]为例)，我们就可以使用二分法进行搜索：对于当前的搜索区间[low,high]，mid=(low+high)/2，根据z(mid)的值及z(rate)的增减性，对搜索区间进行更新：
若z(mid)>0，上调区间，使low=mid+1；

若z(mid)<0，下调区间，使high=mid-1；

从而在经过有限次搜索之后便能找到所求比率。

3、求解生成树xi使(∑xi×ci-rate*∑xi×disi)的值最小的方法：

∑xi×ci-rate*∑xi×disi=∑xi(ci+rate×disi)，从而问题转化为以ci+rate×disi为边的权重，求解最小生成树，对于完全图，可使用

prim算法，其复杂度只与节点数有关。

******************************************************************

简单来说(迭代的不是二分)

1.先设比率r=0，对于每条边，我们计算一个新的权值，权值为w[i][j]=c[i][j]-r*d[i][j] (其中c[i][j]为第i个点和第j个点的垂直高度差，d[i][j]为水平距离,w[i][j]为计算出来的权值)

2.以这个权值去构建最小生成树，用prim算法（时间复杂度只与顶点数有关）去构建。最后统计这个MST的垂直高度差的和sumc，水平距离的和sumd，算出新的比率为R=sumc/sumd;

3.判断新的比率和R和旧的比率r是否相等（精度范围内，这里设为0.00001）,如果相等那么R就是答案，否则就r=R（迭代），然后再次去做（1），依次循环直到找到答案

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------

1 #include<cstdio>
2 #include<cstdlib>
3 #include<cstring>
4 #include<iostream>
5 #include<queue>
6 #include<algorithm>
7 using namespace std;
8
9 const int N=2*450,M=2*11000;
10 const double INF=(double)1e9;
11 int n,m,len,first
,fa
;
12 bool vis
;
13 double F;
14 struct node{
15     int x,y,next;
16     double w,d,t;
17     bool in;
18 }a[M];
19
20 struct cmpp
21 {
22     bool operator () (int &x,int &y)
23     {
24         return a[x].w<a[y].w;
25     }
26 };
27
28 priority_queue<int,vector<int>,cmpp> q;
29
30 void ins(int x,int y,double d,double t)
31 {
32     a[++len].x=x;a[len].y=y;a[len].d=d;a[len].t=t;
33     a[len].next=first[x];first[x]=len;
34 }
35
36 int findfa(int x)
37 {
38     if(fa[x]==x) return fa[x];
39     return findfa(fa[x]);
40 }
41
42 double myabs(double x){return x>0 ? x:-x;}
43
44 bool cmp(node x,node y){return x.w>y.w;}
45
46 bool check(double r)
47 {
48     int cnt=0;
49     double sum=F;
50     while(!q.empty()) q.pop();
51     memset(vis,0,sizeof(vis));
52     for(int i=1;i<=len;i++)
53     {
54         a[i].w=-a[i].d-r*a[i].t;
55         a[i].in=0;
56     }
57     vis[1]=1;
58     for(int i=first[1];i;i=a[i].next) a[i].in=1,q.push(i);
59
60     while(cnt!=n-1)
61     {
62         int ind;
63         while(1)
64         {
65             ind=q.top();q.pop();
66             int x=a[ind].x,y=a[ind].y;
67             if(!vis[x])
68             {
69                 for(int i=first[x];i;i=a[i].next) if(!a[i].in) a[i].in=1,q.push(i);
70                 break;
71             }
72             if(!vis[y])
73             {
74                 for(int i=first[y];i;i=a[i].next) if(!a[i].in) a[i].in=1,q.push(i);
75                 break;
76             }
77         }
78         cnt++;
79         vis[a[ind].x]=vis[a[ind].y]=1;
80         sum+=a[ind].w;
81     }
82     return (sum>=0);
83 }
84
85 int main()
86 {
87     freopen("a.in","r",stdin);
88     // freopen("quake.in","r",stdin);
89     // freopen("quake.out","w",stdout);
90
91     scanf("%d%d%lf",&n,&m,&F);
92     len=0;
93     memset(first,0,sizeof(first));
94     for(int i=1;i<=m;i++)
95     {
96         int x,y;
97         double d,t;
98         scanf("%d%d%lf%lf",&x,&y,&d,&t);
99         ins(x,y,d,t);
100         ins(y,x,d,t);
101     }
102
103     double l=0,r=INF,mid;
104     while(l<r)
105     {
106         mid=(l+r)/2;
107         if(check(mid)) l=mid;
108         else r=mid;
109         if(myabs(l-r)<=0.000001) break;
110     }
111     printf("%.4lf\n",r);
112     return 0;
113 }