poj 2728 Desert King

点击打开链接poj 2728

思路：最优比例生成树

1问题描述：
  已知一个完全图，每条边有两个参数（dis和c），求一棵生成树，使(∑xi×ci)/(∑xi×disi)最小，其中xi当第i条边包含在生成树中时为1，否则为0。

2处理方法：
   1迭代法
    假设rate为当前比率，以ci-rate×disi作为各边的权重，使用Prim算法构造最小生成树，再对该最小生成树求(∑xi×ci)/(∑xi×disi)更新rate，可证明rate可收敛且收敛值即为所求。
   2二分法
    在一个精度范围内（以1e-6为例），二分查找[0,maxRate]之间的值rate，使(z=∑xi×ci-rate×∑xi×disi)==0，可证明该值即为所求。其中xi为以ci-rate×disi作为各边权重时，使用Prim算法所构造的最小生成树。在二分搜索过程中若z>0，则将区间上移(low=mid)，否则将区间下移(high=mid)。


3解题思路及相关证明：
  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值即为所求。

   2下面证明其存在性和唯一性：
      存在性：对于题目所求的最小比率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。
   
    3求解：有两种方法对该问题进行求解：迭代法和二分法。
         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算法，其复杂度只与节点数有关。

代码：

/*二分法*/
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
#define MAXN 1010
#define INF 0xFFFFFFF

int n;
struct Point{
   int x;
   int y;
   int h;
}p[MAXN];
double G[MAXN][MAXN];
double H[MAXN][MAXN];
double lowcost[MAXN];
int vis[MAXN];

/*初始化*/
void init(){
    for(int i = 1 ; i <= n ; i++){
      for(int j = i+1 ; j <= n ; j++){
         double tmp_x = (p[i].x-p[j].x)*(p[i].x-p[j].x);
         double tmp_y = (p[i].y-p[j].y)*(p[i].y-p[j].y);
         G[i][j] = G[j][i] = sqrt(tmp_x + tmp_y);
         H[i][j] = H[j][i] = abs(p[i].h-p[j].h);
      }
    }
}

/*prime算法*/
double prime(double tmp){
     int pos;
     double ans = 0;
     memset(vis , 0 , sizeof(vis));
     vis[1] = 1;
     for(int i = 1 ; i <= n ; i++)
         lowcost[i] = H[1][i]-tmp*G[1][i];
     for(int i = 1 ; i <= n ; i++){
         pos = -1;
         for(int j = 1 ; j <= n ; j++){
            if(!vis[j] && (pos == -1 || lowcost[j] < lowcost[pos]))
              pos = j;
         }
         if(pos == -1)
             break;
         vis[pos] = 1;
         ans += lowcost[pos];
         for(int j = 1 ; j <= n ; j++){
            /*TLE的写法
            if(!vis[j] && lowcost[j] > H[j][pos]-tmp*G[j][pos])
               lowcost[j] = H[j][pos] - tmp*G[j][pos];
            */            
            if(!vis[j] && lowcost[j] > H[pos][j]-tmp*G[pos][j])
               lowcost[j] = H[pos][j] - tmp*G[pos][j];
         }
     }
     return ans;
}

int main(){
    double left , right , mid , tmp , pre;
    while(scanf("%d" , &n) && n){
      for(int i = 1 ; i <= n ; i++)
          scanf("%d%d%d" , &p[i].x , &p[i].y , &p[i].h);
      init();
      left = 0;
      right = 100;/*二分上界*/
      while(left < right){
          mid = (left + right)/2;
          tmp = prime(mid);
          if(fabs(tmp-0.0) <= 0.0001)/*注意这里的精度问题*/
              break;
          else if(tmp > 0.0001)
              left = mid;
          else
              right = mid;
      }
      printf("%0.3lf\n" , mid);
    }
    return 0;
}

/*迭代法*/
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
#define MAXN 1010
#define INF 0xFFFFFFF

int n;
struct Point{
   int x;
   int y;
   int h;
}p[MAXN];
double G[MAXN][MAXN];
double H[MAXN][MAXN];
double lowcost[MAXN];
int pre[MAXN];/*记录前驱点的编号*/
int vis[MAXN];

void init(){
    for(int i = 1 ; i <= n ; i++){
      for(int j = i+1 ; j <= n ; j++){
         double tmp_x = (p[i].x-p[j].x)*(p[i].x-p[j].x);
         double tmp_y = (p[i].y-p[j].y)*(p[i].y-p[j].y);
         G[i][j] = G[j][i] = sqrt(tmp_x + tmp_y);
         H[i][j] = H[j][i] = abs(p[i].h-p[j].h);
      }
    }
}

double prime(double tmp){
     int pos;
     double ansHeight , ansLength;
     ansHeight = ansLength = 0;
     memset(vis , 0 , sizeof(vis));
     vis[1] = 1;
     for(int i = 1 ; i <= n ; i++){
         lowcost[i] = H[1][i]-tmp*G[1][i];
         pre[i] = 1;
     }
     for(int i = 1 ; i <= n ; i++){
         pos = -1;
         for(int j = 1 ; j <= n ; j++){
            if(!vis[j] && (pos == -1 || lowcost[j] < lowcost[pos]))
              pos = j;
         }
         if(pos == -1)
             break;
         vis[pos] = 1;
         ansHeight += H[pre[pos]][pos];
         ansLength += G[pre[pos]][pos];
         for(int j = 1 ; j <= n ; j++){       
            if(!vis[j] && lowcost[j] > H[pos][j]-tmp*G[pos][j]){
               lowcost[j] = H[pos][j] - tmp*G[pos][j];
               pre[j] = pos;
            }
         }
     }
     return ansHeight/ansLength;/*返回比例*/
}

int main(){
    double a , b;
    while(scanf("%d" , &n) && n){
      for(int i = 1 ; i <= n ; i++)
          scanf("%d%d%d" , &p[i].x , &p[i].y , &p[i].h);
      init();
      a = 0;/*初始化为0*/
      while(1){
          b = prime(a);/*求出b*/
          if(fabs(b-a) <= 1e-4)/*满足条件*/
              break;
          a = b;/*重新赋值*/
      }
      printf("%0.3lf\n" , b);/*输出最后的答案*/
    }
    return 0;
}