解题报告 之 POJ2728 Desert King

解题报告 之 POJ2728 Desert King 

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

40 0 00 1 11 1 21 0 30

Sample Output

1.000

  题目大意是有对于三维空间中的一些点，距离len是x,y的欧几里德距离，代价cost是z之差的绝对值，问将所有点连起来的代价和/长度和的最小值。

  根据学长的提示，是一道最优比率生成树的题，学会这道题之后收获很多……特别是坚定了我将图论推给队长的决心……
于是就学习到了0-1分数规划。

 有带权图G, 对于图中每条边e[i], 都有benifit[i](收入)和cost[i](花费), 我们要求的是一棵生成树T, 它使得 ∑(benifit[i]) / ∑(cost[i]),i∈T最小.设x[i]等于1或0, 表示边e[i]是否属于生成树.所以所求比率r=∑(benifit[i]*x[i]) / ∑(cost[i]*x[i])。设答案为 r=∑(benifit[i]*x[i]) / ∑(cost[i]*x[i])。则对于每一颗生成树，有∑(benifit[i]*x[i])-r*∑(cost[i]*x[i])=0。所以对于所有生成树，有∑(benifit[i]*x[i])-rmin*∑(cost[i]*x[i])>=0。所以为了使得总权值能够==0，求最小生成树。之后二分r：
 

                                            当 r>rmin时，
                                            某序列x1,x2……xn使得  ∑(benifit[i]*x[i])-r* ∑(cost[i]*x[i])<0
                                            当 r<rmin时，
                                            某序列x1,x2……xn使得  ∑(benifit[i]*x[i])-r* ∑(cost[i]*x[i])>0
                                            当 r=rmin时，
                                            某序列x1,x2……xn使得  ∑(benifit[i]*x[i])-r* ∑(cost[i]*x[i])=0  
                                            即为所求。 
 

  所以我们可以令dis[i]= cost[i] - l*len[i] 作为图中边的权值，而自然z=∑(cost[i]*x[i]) - l* ∑(len[i]*x[i])，在生成最小生成树。于是我们可以用Prime算法求出这个奇葩的东西……最后用二分法rmin即可。有另一种Dinkelbach仍然不明觉厉，于是我就记住了一个模板呵呵呵。。。。。。Orz

#include<cstdio>

#include<cmath>

#include<cstring>

#define MAXN 1010

#define INF 10000000

#define bisection   //用于切换用二分法和Dinkelbach迭代法

using namespace std;

 

double dis[MAXN];   // dis[i]=cost[i] - l*len[i]  作为边的权值

int vis[MAXN], pre[MAXN]; 

int n;

 

struct NODE

{

      double x, y, z;

};  

 

NODE nodes[MAXN];

 

double distance(NODE a, NODE b)

{

      return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));    //两结点距离公式

}

 

double prime(double l)                      //变形Prime算法构建最小生成树，返回z(l)的值 

{

      memset(vis, 0, sizeof(vis));

      for (int i = 1; i <= n; i++)

      {

            dis[i] = fabs(nodes[1].z - nodes[i].z) - distance(nodes[1], nodes[i])*l ;

            pre[i] = 1;

      }

      double sum = 0.0;

      int pos = 1;

      vis[1] = 1;

      dis[1] = 0;

      double cost=0, len = 0;

      for (int i = 1; i < n; i++)

      {

            double MIN = INF;

            for (int j = 1; j <= n; j++)

                  if (vis[j] == 0 && MIN > dis[j])

                  {

                      MIN = dis[j];

                      pos = j;

                  }

            sum += MIN;

            vis[pos] = 1;

            cost += fabs(nodes[pre[pos]].z - nodes[pos].z);      //用于Dinkelbach算法
                                                                                        //由于题目特殊性，此处cost代表公式中的benifit

            len += distance(nodes[pre[pos]], nodes[pos]);          //len代表公式中的cost

            for (int j = 1; j <= n; j++)

            {

                  double val = fabs(nodes[pos].z - nodes[j].z) - distance(nodes[pos], nodes[j])*l;

                  if (vis[j] == 0 && dis[j] > val)

                  {

                        dis[j] = val;

                        pre[j] = pos;

                  }

            }

      }

#ifdef bisection

      return sum;

#else

      return cost / len;

#endif                              //用于切换二分法和Dinkelbach法

}

 

int main()

{

      while (scanf("%d", &n) == 1 && n)

      {

            for (int i = 1; i <= n; i++)

                  scanf("%lf%lf%lf", &nodes[i].x, &nodes[i].y, &nodes[i].z);

 

            double a = 0, b = 0;

#ifdef bisection                                             //用二分法求r的值

            double head = 0, tail = 10000.0;

            while (tail - head > 1e-5)                    // 注意精度

            {

                  double mid = (tail+head) / 2.0;

                  a = prime(mid);

                  if (a >= 0)

                        head=mid;

                  else

                        tail = mid;

 

            }

            printf("%.3lf\n", tail);

 

#else

            while (1)

            {

                  b = prime(a);

                  if (fabs(b - a) < 1e-4)   break;               //Dinkelbach迭代法……原理不明，求大神指点Orz

                  a = b;

            }

            printf("%.3lf\n", b);

#endif

      }

      return 0;

}

其实是一篇以前的解题报告。。发现自己看不懂。。囧