最近在刷巨巨们放出来的专题,然后没做几题就卡住了,果然还是太弱了T U T...
这次做到了一题01分数规划求解的生成树问题。
题目大意是这样的:给你一个无向完全图,每条边i都有两个权值,长度a[ i ],花费b[ i ],需要选出其中的一些边构造一颗生成树,生成树需要满足条件:∑ b [ i ] / ∑ a [ i ]最小。
这样我还是先来介绍一下01分数规划吧~
给定一个上述的问题,我们可以设R = ∑x[ i ] * b[ i ] / ∑x[ i ] * a[ i ],其中x[ i ]为0或1表示边i取或者不取。
这里我们先假设对于所有的可行解x[ i ],a[ i ] * x[ i ]都是正数。
我们需要求R的最小值,那么不妨先设一个子问题Q(L) = ∑x[ i ] * b[ i ] - L * ∑x[ i ] * a[ i ] 的最小值,然后我们需要求满足条件(即它是一棵生成树)的L的最小值。
易知L就是R,所以问题等价于求R的最小值
首先我们先对Q(L)移项得:Q(L) = ∑x[ i ] * ( b[ i ] - R * a[ i ] )。
如果我们已知L,那么b[ i ] - L * a[ i ]也是已知的,所以子问题只和两个参数L,x[ i ]有关。并且对于一个L,x[ i ]的选取不同,将可能会导致Q(L)的取值不同。
现在我们假设Q(L) < 0 ,那么我们能得到什么信息?
根据Q(L) = ∑x[ i ] * ( b[ i ] - L * a[ i ] ) < 0 可以推出∑x[ i ] * b[ i ] / ∑x[ i ] * a[ i ] < L ,那么就说明我们找到了一个比R更优的解!为什么不选择更优一点的解呢~
如果Q(L) < 0,那么这时候比L小的解都是不是最优的,自然不需要了。(其实应该理解为它们都是意义的)
当然如果Q(L) = 0,那么我们找到了一个可行解R=L,但是不是最优解,暂时我们还不知道。
如果函数Q(L)满足单调性,则说明如果存在Q(L) = 0,则有且只有这一个解,并且这个解就是我们需要的答案。
证明:
假设L2 > L1,令X[ i ]是L2的一组解向量,Y[ i ]是L1的一组解向量。
那么Q(L2) = ∑X[ i ] * b[ i ] - L2 * ∑X[ i ] * a[ i ] , Q(L1) =∑Y[ i ] * b[ i ] - L1 * ∑Y[ i ] * a[ i ]。
易知∑Y[ i ] * b[ i ] - L2 * ∑Y[ i ] * a[ i ] >= ∑X[ i ] * b[ i ] - L2 * ∑X[ i ] * a[ i ](因为X[ i ],Y[ i ]分别是Q(L2),Q(L1)最小值对应的解向量,所以Y[ i ]放到L2中去一定不会比最小值还小。
所以有:Q(L2) - Q(L1) = ( ∑X[ i ] * b[ i ] - L2 * ∑X[ i ] * a[ i ] ) -( ∑Y[ i ] * b[ i ] - L1 * ∑Y[ i ] * a[ i ] ) <=
( ∑Y[ i ] * b[ i ] - L2 * ∑Y[ i ] * a[ i ] ) - ( ∑Y[ i ] * b[ i ] - L1 * ∑Y[ i ] * a[ i ] ) =
( L1 - L2 ) * ∑Y[ i ] * a[ i ] < 0 ( Y[ i ] > 0,a[ i ] > 0 )
因此我们得到了对于L2 > L1,Q(L2) < Q(L1),所以函数Q(L)单调递减。
所以根据Q(L)的单调性,我们很容易知道当Q(L) = 0的时候,R = L 取得最小值。
另外,假设L*是Q(L) = 0时的最优解:
Q(L) > 0
当且仅当
L < L*
Q(L) = 0
当且仅当
L = L*
Q(L) < 0
当且仅当
L > L*
现在我们可以通过二分法来寻找R值或者通过Dinkelbach算法迭代求解。
不同的情况下,两种算法会发挥不同的效果,视情况而定。
题目大意:给你一个无向完全图,每条边i都有两个权值,长度a[ i ],花费b[ i ],需要选出其中的一些边构造一颗生成树,生成树需要满足条件:∑ b [ i ] / ∑ a [ i ]最小。
题目分析:01分数规划——最优比率生成树。经典题!!!
本题二分明显没有迭代的效率高。
二分解法:
求出Q(L) = height[ i ] - R * length[ i ]
Q(L) <= 0 则移动上界,否则移动下界
二分代码:
//1454ms
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define REPF( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define REPV( i , a , b ) for ( int i = a ; i >= b ; -- i )
#define clear( a , x ) memset ( a , x , sizeof a )
const int MAXN = 1005 ;
const double INF = 1e18 ;
const double eps = 1e-5 ;
struct Point {
double x , y , h ;
void input () {
scanf ( "%lf%lf%lf" , &x , &y , &h ) ;
}
} ;
struct MST {
Point point[MAXN] ;
double G[MAXN][MAXN] ;
double length[MAXN][MAXN] ;
double height[MAXN][MAXN] ;
double d[MAXN] ;
bool vis[MAXN] ;
int fa[MAXN] ;
int n ;
double dist ( double x , double y ) {
return sqrt ( x * x + y * y ) ;
}
void input () {
REP ( i , n )
point[i].input () ;
REPF ( i , 0 , n - 1 )
REPF ( j , i + 1 , n - 1 ) {
length[i][j] = length[j][i] = dist ( point[i].x - point[j].x , point[i].y - point[j].y ) ;
height[i][j] = height[j][i] = fabs ( point[i].h - point[j].h ) ;
}
REP ( i , n )
height[i][i] = length[i][i] = G[i][i] = 0 ;
}
void cal ( double r ) {
REPF ( i , 0 , n - 1 )
REPF ( j , i + 1 , n - 1 )
G[i][j] = G[j][i] = height[i][j] - r * length[i][j] ;
}
double prim ( double r ) {
cal ( r ) ;
REP ( i , n )
d[i] = INF ;
clear ( vis , 0 ) ;
d[0] = 0 ;
fa[0] = 0 ;
double ans = 0 , m ;
int x ;
REP ( _ , n ) {
m = INF ;
REP ( i , n )
if ( !vis[i] && m > d[i] )
m = d[x = i] ;
vis[x] = 1 ;
ans += m ;
REP ( i , n )
if ( !vis[i] && d[i] > G[x][i] )
d[i] = G[x][i] , fa[i] = x ;
}
return ans ;
}
void solve () {
input () ;
double mid , l = 0 , r = 100 ;
while ( fabs ( r - l ) > eps ) {
mid = ( l + r ) / 2 ;
if ( prim ( mid ) <= eps )
r = mid ;
else
l = mid ;
}
printf ( "%.3f\n" , l ) ;
}
} ;
MST z ;
int main () {
while ( ~scanf ( "%d" , &z.n ) && z.n )
z.solve () ;
return 0 ;
}
Dinkelbach算法:
设初值0,第一次带入求出的Q(L)一定大于0,但是得到的L带入第二次得到的一定小于0,然后就不断的迭代直到Q(L) = 0。
迭代代码:
//282ms
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define REPF( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define REPV( i , a , b ) for ( int i = a ; i >= b ; -- i )
#define clear( a , x ) memset ( a , x , sizeof a )
const int MAXN = 1005 ;
const double INF = 1e18 ;
const double eps = 1e-5 ;
struct Point {
double x , y , h ;
void input () {
scanf ( "%lf%lf%lf" , &x , &y , &h ) ;
}
} ;
struct MST {
Point point[MAXN] ;
double G[MAXN][MAXN] ;
double length[MAXN][MAXN] ;
double height[MAXN][MAXN] ;
double d[MAXN] ;
bool vis[MAXN] ;
int fa[MAXN] ;
int n ;
double dist ( double x , double y ) {
return sqrt ( x * x + y * y ) ;
}
void input () {
REP ( i , n )
point[i].input () ;
REPF ( i , 0 , n - 1 )
REPF ( j , i + 1 , n - 1 ) {
length[i][j] = length[j][i] = dist ( point[i].x - point[j].x , point[i].y - point[j].y ) ;
height[i][j] = height[j][i] = fabs ( point[i].h - point[j].h ) ;
}
REP ( i , n )
height[i][i] = length[i][i] = G[i][i] = 0 ;
}
void cal ( double r ) {
REPF ( i , 0 , n - 1 )
REPF ( j , i + 1 , n - 1 )
G[i][j] = G[j][i] = height[i][j] - r * length[i][j] ;
}
double prim ( double r ) {
cal ( r ) ;
REP ( i , n )
d[i] = INF ;
clear ( vis , 0 ) ;
d[0] = 0 ;
fa[0] = 0 ;
double a = 0 , b = 0 , m ;
int x ;
REP ( _ , n ) {
m = INF ;
REP ( i , n )
if ( !vis[i] && m > d[i] )
m = d[x = i] ;
a += length[fa[x]][x] ;
b += height[fa[x]][x] ;
vis[x] = 1 ;
REP ( i , n )
if ( !vis[i] && d[i] > G[x][i] )
d[i] = G[x][i] , fa[i] = x ;
}
return b / a ;
}
void solve () {
input () ;
double res = 0 , tmp ;
while ( 1 ) {
tmp = prim ( res ) ;
if ( fabs ( tmp - res ) <= eps )
break ;
res = tmp ;
}
printf ( "%.3f\n" , res ) ;
}
} ;
MST z ;
int main () {
while ( ~scanf ( "%d" , &z.n ) && z.n )
z.solve () ;
return 0 ;
}