poj 2728 Desert King(最优比例生成树)

poj 2728 Desert King(最优比率生成树)

链接：

http://poj.org/problem?id=2728

题目：

Desert King

Time Limit: 3000MS		Memory Limit: 65536K
Total Submissions: 16802		Accepted: 4732

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

Sample Output

1.000

Source

Beijing 2005

题外话： 2005年北京赛区，楼教主（ACRush）第5分钟AC了一题之后，在第25分钟又AC了这题，结果误导了其他队伍也都也做这题，结果最终这题提交了200多次，但只有8个队伍AC了。膜拜Orz...

分析与总结：

通过这题学了最优比率生成树问题。

有二分和迭代两种方法，当年楼教主用的就是二分的，但是速度较慢，poj用1200ms，而迭代只用了250ms。

不过迭代的原理和证明还不太懂，而黑书干脆以“证明很复杂，这里略去”为由不讲了，囧rz....

参考资料：黑书 + 百度

若记代价/距离为r，则，其中Xi=0或1，表示取不取i这条边。

记，则Z(r)=0时，r为解（注意这个式子怎么变来的）。

那么我们将边权变为求生成树，如果Z(r)=0，则r就是我们的解。

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

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) /*由z(rate)的定义，xi为一棵生

成树使(∑xi×ci-rate*∑xi×disi)的值最小，故有z(rate1)= (∑xi×ci-rate1*∑xi×disi)<= (∑yi×ci-rate1*∑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.二分

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define Q(x) ((x)*(x))
using namespace std;

const int VN = 1005;
const int EN = VN*VN/2;
const double eps = 1e-4;

int n;
double cost[VN][VN],len[VN][VN],x[VN], y[VN],z[VN];
  
template<typename Type>
class Prim{
public:
    void init(int _n){
        n=_n; MSTsum = 0;
    }
    bool prim(double m){
        for(int i=1; i<n; ++i){ //初始化权值
            for(int j=i+1; j<=n; ++j){
                w[i][j]=w[j][i]=cost[i][j]-m*len[i][j];
            }
        }
        memset(vis, 0, sizeof(vis));
        vis[1] = true;
        MSTsum = 0;
        for(int i=1; i<=n; ++i){
            key[i] = w[1][i]; pre[i] = 1; 
        }
        for(int i=1; i<n; ++i){
            int u=-1;
            for(int j=1; j<=n; ++j)if(!vis[j]){
                if(u==-1||key[j]<key[u]) u=j;
            }
            MSTsum += w[pre[u]][u];
            vis[u] = true;
            for(int j=1; j<=n; ++j)if(!vis[j]){
                if(key[j]>w[u][j]){
                    key[j] = w[u][j]; pre[j] = u;
                }
            }
        }
        if(MSTsum>0.0) return true;
        return false;
    }
   
private:
    Type MSTsum;
    Type w[VN][VN], key[VN];
    int pre[VN], n;
    bool vis[VN];
};
Prim<double>G;
inline double getDist(double x1,double y1,double x2,double y2){
    return sqrt(Q(x1-x2)+Q(y1-y2));
}

int main(){
    while(~scanf("%d",&n)&&n){
        G.init(n);
        for(int i=1; i<=n; ++i)
            scanf("%lf%lf%lf",&x[i],&y[i],&z[i]);
        //初始化
        for(int i=1; i<n; ++i){
            cost[i][i]=len[i][i]=0;
            for(int j=i+1; j<=n; ++j){
                double d=getDist(x[i],y[i],x[j],y[j]);
                cost[i][j] = cost[j][i] = fabs(z[i]-z[j]);
                len[i][j] = len[j][i] = d;
            }
        }
        double left=0, right=100;
        while(right-left>=eps){
            double mid=(left+right)/2.0;
            if(G.prim(mid)){
                left=mid;
            }
            else{
                right=mid;
            }
        }
        printf("%.3f\n", left);
    }
    return 0;
}

2.迭代

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define Q(x) ((x)*(x))
using namespace std;

const int VN = 1005;
const int EN = VN*VN/2;
const double eps = 1e-4;

int n;
double cost[VN][VN],len[VN][VN],x[VN], y[VN],z[VN];
  
template<typename Type>
class Prim{
public:
    void init(int _n){
        n=_n; MSTsum = 0;
    }
    double prim(double m){
        for(int i=1; i<n; ++i){ //初始化权值
            for(int j=i+1; j<=n; ++j){
                w[i][j]=w[j][i]=cost[i][j]-m*len[i][j];
            }
        }
        memset(vis, 0, sizeof(vis));
        vis[1] = true;
        MSTsum = 0;
        double c=0, l=0;
        for(int i=1; i<=n; ++i){
            key[i] = w[1][i]; pre[i] = 1; 
        }
        for(int i=1; i<n; ++i){
            int u=-1;
            for(int j=1; j<=n; ++j)if(!vis[j]){
                if(u==-1||key[j]<key[u]) u=j;
            }
            MSTsum += w[pre[u]][u];
            vis[u] = true;

            c += cost[pre[u]][u];
            l += len[pre[u]][u];

            for(int j=1; j<=n; ++j)if(!vis[j]){
                if(key[j]>w[u][j]){
                    key[j] = w[u][j]; pre[j] = u;
                }
            }
        }
        return c/l;
    }
   
private:
    Type MSTsum;
    Type w[VN][VN], key[VN];
    int pre[VN], n;
    bool vis[VN];
};
Prim<double>G;
inline double getDist(double x1,double y1,double x2,double y2){
    return sqrt(Q(x1-x2)+Q(y1-y2));
}

int main(){
    while(~scanf("%d",&n)&&n){
        G.init(n);
        for(int i=1; i<=n; ++i)
            scanf("%lf%lf%lf",&x[i],&y[i],&z[i]);
        //初始化
        for(int i=1; i<n; ++i){
            cost[i][i]=len[i][i]=0;
            for(int j=i+1; j<=n; ++j){
                double d=getDist(x[i],y[i],x[j],y[j]);
                cost[i][j] = cost[j][i] = fabs(z[i]-z[j]);
                len[i][j] = len[j][i] = d;
            }
        }
        double rate = 0, r;
        do {
            r = rate;
            rate=G.prim(rate);
        } while (fabs(r - rate) > eps);
        printf("%.3f\n", r);
    }
    return 0;
}


——  生命的意义，在于赋予它意义。

          
     原创 http://blog.csdn.net/shuangde800 ， By   D_Double  (转载请标明)

poj 2728 Desert King(最优比例生成树)

相关推荐