POJ 2728 Desert King（最优比率生成树）

Desert King
Time Limit: 3000MS
Memory Limit: 65536K
Total Submissions: 26364
Accepted: 7307
Description
David the Great has just become the king of a desert country. To win the respect of his people, he decided to build channels all over his country to bring water to every village. Villages which are connected to his capital village will be watered. As the dominate ruler and the symbol of wisdom in the country, he needs to build the channels in a most elegant way.

After days of study, he finally figured his plan out. He wanted the average cost of each mile of the channels to be minimized. In other words, the ratio of the overall cost of the channels to the total length must be minimized. He just needs to build the necessary channels to bring water to all the villages, which means there will be only one way to connect each village to the capital.

His engineers surveyed the country and recorded the position and altitude of each village. All the channels must go straight between two villages and be built horizontally. Since every two villages are at different altitudes, they concluded that each channel between two villages needed a vertical water lifter, which can lift water up or let water flow down. The length of the channel is the horizontal distance between the two villages. The cost of the channel is the height of the lifter. You should notice that each village is at a different altitude, and different channels can’t share a lifter. Channels can intersect safely and no three villages are on the same line.

As King David’s prime scientist and programmer, you are asked to find out the best solution to build the channels.
Input
There are several test cases. Each test case starts with a line containing a number N (2 <= N <= 1000), which is the number of villages. Each of the following N lines contains three integers, x, y and z (0 <= x, y < 10000, 0 <= z < 10000000). (x, y) is the position of the village and z is the altitude. The first village is the capital. A test case with N = 0 ends the input, and should not be processed.
Output
For each test case, output one line containing a decimal number, which is the minimum ratio of overall cost of the channels to the total length. This number should be rounded three digits after the decimal point.
Sample Input
4
0 0 0
0 1 1
1 1 2
1 0 3
0
Sample Output
1.000
Source
Beijing 2005

题目大意：

　　给你一张无向图，告诉你任意两点的距离和连接的花费，求一棵生成树使得花费／生成树上的边距离之和最小。

解题思路：

　　一道裸的最优比例生成树，也是01分数规划的一个经典应用。
　　
　　首先给出分数规划的一般形式：$\lambda =\frac{a(x)}{b(x)} (x\in S)\\s.t.\forall x\in S, b(x)>0$
其中，解向量x在解空间S内，a(x)与b(x)都是连续的实值函数。我们要求的就是$\lambda$的最小或最大值。

　　01分数规划：给定两个数组，a[i]表示选取i的收益，b[i]表示选取i的代价。如果选取i，定义x[i]=1否则x[i]=0。每一个物品只有选或者不选两种方案，求一个选择方案使得$R=\frac{\sum a[i] \times x[i]}{\sum b[i] \times x[i]}$取得最值，即所有选择物品的总收益/总代价的值最大或是最小。
　　
　　可以发现本题的最优比例生成树很明显是一个01分数规划问题。
　　
　　对于01分数规划的求解。我们通常将目标式变形，定义函数$f(R)=\sum (a[i]-R \times b[i]) \times x[i]$。由于a[i]与b[i]是已知的，那么$f(R)$只与R和x数组有关。可以证明$f(R)$是满足单调性的。当$f(R)$大于0的时候说明我们可以通过当前x数组的取值得到一个更大的R，当$f(x)$小于0的时候说明可以得到一个更小的R。所以我们就可以使用二分R来逐步逼近答案，对于当前R如果可以找到一个x的取值使得$f(R)$大于0就移动下界，否则移动上界。
　　不过我们可以发现，很多题目中在我们计算$f(R)$的过程中是可以得到一个x数组的取值的，如果$f(R)$大于0，那么我们利用组x求得的R一定比原先的R更大，$f(R)$小于0类似。比起漫无目的的二分我们可以直接把答案更新到这个新的R。这就是Dinkelbach算法的基本思想。他并不会去二分答案，而是先随便给定一个答案，然后根据更优的解不断移动答案，逼近最优解。由于他对每次判定使用的更加充分，所以它比二分会快上很多。但是有的时候判断一个解和计算新的解的复杂度是不一样的，所以在某些情况下我们还是要使用二分求解。

　　最后回到这道题的具体实现。我们可以得到$f(R)=\sum (dz[i]-dx[i] \times R)\times x[i]$其中$dz$为边的高度差，$dx$为边的水平距离，x的选择必须是一棵生成树。那么对于当前的L我们可以把$dz[i]-dx[i]\times R$作为一条边的权值，这样的到的生成树的$f(R)$当前R下所有x取值中最小的，当然由这组解得到的新的R也已定是可以有当前R得到的最小的。于是我们就可以用二分或迭代得到答案。

AC代码（二分）

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <vector>
#include <queue>
#include <stack>
#include <deque>
#include <string>
#include <map>
#include <set>
#include <list>
using namespace std;
#define INF 0x3f3f3f3f
#define LL long long
#define fi first
#define se second
#define mem(a,b) memset((a),(b),sizeof(a))

const double eps=1e-5;
const int MAXV=1000+3;
double y[MAXV], x[MAXV], z[MAXV];
double G[MAXV][MAXV];//顶点之间的距离
int V;//顶点数
double lowcost[MAXV];//到当前生成树的最小花费
int nearvex[MAXV];//当前生成树中最近的点

inline double get_dist(int i, int j)
{
    return sqrt((y[i]-y[j])*(y[i]-y[j])+(x[i]-x[j])*(x[i]-x[j]));
}

double prim(int s, double l)//起点，当前答案
{
    double sum=0;
    for(int i=0;i<V;++i)
    {
        nearvex[i]=s;
        lowcost[i]=fabs(z[s]-z[i])-G[s][i]*l;
    }
    nearvex[s]=-1;
    while(true)
    {
        double the_min=INF;
        int v=-1;
        for(int i=0;i<V;++i)
            if(nearvex[i]!=-1 && lowcost[i]<the_min)
            {
                v=i;
                the_min=lowcost[i];
            }
        if(v==-1)//没有可以加入生成树的点
            break;
        nearvex[v]=-1;
        sum+=lowcost[v];
        for(int i=0;i<V;++i)
        {
            double tmp=fabs(z[v]-z[i])-G[v][i]*l;
            if(nearvex[i]!=-1 && tmp<lowcost[i])
            {
                lowcost[i]=tmp;
                nearvex[i]=v;
            }
        }
    }
    return sum;
}

int main()
{
    while(~scanf("%d", &V) && V)
    {
        for(int i=0;i<V;++i)
            scanf("%lf%lf%lf", &y[i], &x[i], &z[i]);
        for(int i=0;i<V;++i)
            for(int j=0;j<V;++j)
                G[i][j]=get_dist(i, j);
        double l=0.0, r=1000000.0;
        while(r-l>eps)//二分查找最优比率
        {
            double mid=(l+r)/2;
            if(prim(0, mid)>=0)
                l=mid;
            else r=mid;
        }
        printf("%.3f\n", r);
    }

    return 0;
}

AC代码（Dinkelbach算法）

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <vector>
#include <queue>
#include <stack>
#include <deque>
#include <string>
#include <map>
#include <set>
#include <list>
using namespace std;
#define INF 0x3f3f3f3f
#define LL long long
#define fi first
#define se second
#define mem(a,b) memset((a),(b),sizeof(a))

const double eps=1e-5;
const int MAXV=1000+3;
double y[MAXV], x[MAXV], z[MAXV];
double G[MAXV][MAXV];//顶点之间的距离
int V;//顶点数
double lowcost[MAXV];//到当前生成树的最小花费
int nearvex[MAXV];//当前生成树中最近的点

inline double get_dist(int i, int j)
{
    return sqrt((y[i]-y[j])*(y[i]-y[j])+(x[i]-x[j])*(x[i]-x[j]));
}

double prim(int s, double l)//起点，当前答案
{
    double cost=0, len=0;//生成树的花费，长度
    for(int i=0;i<V;++i)
    {
        nearvex[i]=s;
        lowcost[i]=fabs(z[s]-z[i])-G[s][i]*l;
    }
    nearvex[s]=-1;
    while(true)
    {
        double the_min=INF;
        int v=-1;
        for(int i=0;i<V;++i)
            if(nearvex[i]!=-1 && lowcost[i]<the_min)
            {
                v=i;
                the_min=lowcost[i];
            }
        if(v==-1)//没有可以加入生成树的点
            break;
            cost+=fabs(z[nearvex[v]]-z[v]);
            len+=G[nearvex[v]][v];
        nearvex[v]=-1;
        for(int i=0;i<V;++i)
        {
            double tmp=fabs(z[v]-z[i])-G[v][i]*l;
            if(nearvex[i]!=-1 && tmp<lowcost[i])
            {
                lowcost[i]=tmp;
                nearvex[i]=v;
            }
        }
    }
    return cost/len;//得到当前最优比率
}

int main()
{
    while(~scanf("%d", &V) && V)
    {
        for(int i=0;i<V;++i)
            scanf("%lf%lf%lf", &y[i], &x[i], &z[i]);
        for(int i=0;i<V;++i)
            for(int j=0;j<V;++j)
                G[i][j]=get_dist(i, j);
        double last=0, now;
        while(true)//迭代
        {
            now=prim(0, last);
            if(fabs(last-now)<eps)
                break;
            last=now;
        }
        printf("%.3f\n", now);
    }

    return 0;
}