0/1分数规划？我不会啊！

时间  2019-11-17
标签 1分 规划 不会 繁體版
原文   原文链接

最近入门了 01 分数规划，因而写篇博客记念html

(若是cnblogs上的格式您实在受不了，请转至个人洛谷blog，而且那里的讲解会稍微详细一些，另外，更新内容也在个人洛谷博客中)ios

分数规划是一项不经常使用到的（但还蛮实用的）算法，通常来说就是求一个最优比率。git

分数规划的简单介绍

分数规划顾名思义就是求一个分数表达式的最大（小）值。算法

好比说有 n 个物品，每一个物品有两个权值 a 和 b ，而后让你选出任意件数（但可能会有限制）的物品，使得两个权值和间的比值最大，即求 \(\dfrac{\sum_{i=1}^{k} a[i]}{\sum_{j=1}^{k} b[j]}\) （在这里 \(1-k\) 为挑出的 k 件物品）的最大值，而后对选择物品方面给出必定的限制条件，那么一道 0/1 分数规划的题目就完成了网络

分数规划的求解方法

分数规划有两种求解方法，一种是 二分法，而另外一种则是 Dinkelbach 算法，通常来说咱们选用第一种方法进行分数规划求解函数

1.二分法

问题同上，求 \(\dfrac{\sum_{i=1}^{k} a[i]}{\sum_{j=1}^{k} b[j]}\) 的最大值，而后加上一个限制：\(k>=W\)post

咱们令 \(sum=\sum_{i=1}^{k} a[i] ,\ tot=\sum_{j=1}^{k} b[j]\) ，\(k>=W\)优化

而后假设问题中的最优解为 \(ans\) ，那么必然有：
\[\dfrac{sum}{tot}<=ans\]ui

移项得：
\[sum<=ans\times tot\]spa

继续移就获得：

\[sum-ans\times tot<=0\]

这样转化有什么用呢？那咱们尝试将 \(sum\) 和 \(tot\) 带回去，就能够获得这么一个式子：

\[\sum_{i=1}^{k} (a[i]-b[i]\times ans) <=0\]

这个式子不难理解，就是把总体的贡献转化为了单件物品的贡献。

那么咱们只须要二分这个 ans， 计算出每件物品的 \(a-b\times ans\)，而后排个序，贪心取前 \(W\) 个加起来，看看最后的值是否 \(<=0\) ，而后就能够根据结果移动左右边界了。

2.Dinkelbach 算法

这个算法其实就是以二分函数图像的单调性为原理，利用了check中计算得出的结果，经过改变二分枚举的中间点在图像上的位置来优化二分。

而后先给一张图：

在这里咱们看到每次枚举的中间点都在一条直线上（但单一的直线并非二分的函数图像），那么对于枚举出的中间点其实能够不着急扔，先计算出当前点所在的直线，而后找到这条直线在 x 轴上的截距，用这个截距做为下一次二分的中间点，这样可能会更快逼近正确答案

尽管理解起来不难，但仍是配套图吧。下面是算法实现的图组

这里显然红点就是答案

这里咱们第一次能够直接二分中点，获得一个初始点

而后咱们过这个获得的点作一条直线，与函数图像相切，

接着咱们将该直线与 \(x\) 轴交点的 \(x\) 坐标做为新点的 \(x\) 坐标（好像离答案更远了啊）

而后新点重复上述过程，作直线

这时咱们发现已经获得答案了

然鹅这个算法的效果却不见的有多么优秀，由于二分的函数图像有多是坑坑洼洼的，因此有时这个算法跑数据的时间反而比二分大，可是若是图像是较为光滑的弧线，或许\(Dinkelbach\) 算法就能充分展示它的优点了

至于其余的地方（包括证实），和二分其实差很少，就不啰嗦了

分数规划的有机结合

分数规划通常来说不会单独成题，通常来说有如下几种形式：

0.不与任何算法结合，即分数规划裸题

1.与\(0/1\)背包结合，即最优比率背包

2.与生成树结合，即最优比率生成树

3.与负环断定结合，即最优比率环

4.与网络流结合，即最大密度子图

5.与费用流结合，即最优比率流（这个是我乱叫的）

6.与其余的各类带选择的算法乱套，即最优比率啥啥的...

对于上面的后三个结合图论算法，其实相似是把供选择的物品变为了图中的边，而后限制条件求解最优比率。

0.分数规划裸题

这个在介绍二分算法的时候讲过了

题目

Dropping tests

代码

//by Judge
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define eps 1e-4
using namespace std;
const int M=1005;
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
char buf[1<<21],*p1=buf,*p2=buf;
inline int read(){ int x=0,f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} int n,k; double a[M],b[M],c[M],l,r=1e6,mid,sum;
inline bool check(){
    for(int i=1;i<=n;++i)
        c[i]=a[i]-b[i]*mid;
    sort(c+1,c+1+n),sum=0;
    for(int i=k;i<=n;++i)
        sum+=c[i];
    return sum>=0;
}
int main(){
    while(1){
        n=read(),k=read();
        if(!n&&!k) return 0;
        l=0,r=1e6,++k;
        for(int i=1;i<=n;++i)
            a[i]=read();
        for(int i=1;i<=n;++i)
            b[i]=read();
        while(r-l>eps){
            mid=(l+r)/2;
            check()?l=mid:r=mid;
        }
        printf("%.0lf\n",100.0*l);
    }
    return 0;
}

1.最优比率背包

题目大意和裸题的差很少，也是求 n 个物品中，\(\dfrac{\sum_{i=1}^{k} a[i]}{\sum_{j=1}^{k} b[j]}\) 的最大值，可是条件有变化，此次咱们要对权值 b 进行限制，首先咱们仍是令 \(sum=\sum_{i=1}^{k} a[i] ,\ tot=\sum_{j=1}^{k} b[j]\) ，那么限制条件就是 \(tot>=W\)（\(W\)已给出）

那么这个时候咱们设完 \(ans\) 以后获得的式子一样是：
\[\sum_{i=1}^{k} (a[i]-b[i]\times ans) <=0\]

因而 \(ans\) 照常二分，对于条件成立的判断方式变一变就行了。

其实你应该也猜到了判断方式，原先咱们不是贪心取前面的 k 个数么，那此次咱们把贪心改为背包就行了啊

那样的话，对于每一个物品来讲，体积就是 b ，而价值就是上面的 \(a[i]-b[i]\times ans\)，而后咱们跑 0/1 背包，判断满背包价值的正负性。

题目

Talent Show

代码

//by Judge
#include<cstdio>
#include<cstring>
#include<iostream>
#define ll long long
using namespace std;
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
inline void cmax(ll& a,ll b){if(a<b)a=b;}
inline int Min(int a,int b){return a<b?a:b;}
char buf[1<<21],*p1=buf,*p2=buf;
inline int read(){ int x=0,f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} int n,l,r=1e6,mid,V,t[305]; ll w[305],f[10005];
inline bool check(int tim){
    memset(f,-0x3f,sizeof(f)),f[0]=0;
    ll tmp=f[V];
    for(int i=1,j,k;i<=n;++i)
        for(j=V;~j;--j) if(f[j]!=tmp){
            cmax(f[Min(j+w[i],V)],f[j]+t[i]-w[i]*tim);
        }
    return f[V]>=0;
}
inline int calc(){
    while(l<=r){
        mid=l+r>>1;
        check(mid)?l=mid+1:r=mid-1;
    }
    return r;
}
int main(){
    n=read(),V=read();
    for(int i=1;i<=n;++i){
        w[i]=read(),
        t[i]=read()*1000;
    }
    return !printf("%d\n",calc());
}

2.最优比率生成树

对于构造最优比率生成树的分数规划，其实差很少相似是将 \(n\) 物品转换为 \(m\) 条边，而后求解，只不过限制改了，选出的 \(n-1\) 条边要构成一棵树

具体点说，题目通常会给你一张图（有时候是彻底图），而后给你 n 个点 m 条边，求出原图包含 n 个点的一棵树，使得所选的边的两个权值和比值最大

那么求法呢，其实也很相似：

咱们考虑构造最小生成树的时候，边是从小到大加的

那么对于每一条边，咱们让他的边权为 \(a[i]-b[i]\times ans\)

而后和克鲁斯卡尔同样排个序一条一条尝试加就行了（然鹅有时候你得尝试普里姆算法）

若是加成功了就累计入 \(sum\)

最后根据 \(sum\) 的正负性来判断 \(check\) 的返回值

其实和前面两个没高明到哪里去，基本想法同样的

题目

Desert King

代码

（慢着，kruskal会炸？仍是我打开方式不对？只好开 \(Prim\) ）

//by Judge
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#define eps 1e-5
using namespace std;
const int M=1005;
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
char buf[1<<21],*p1=buf,*p2=buf;
inline bool cmin(double& a,double b){return a>b?a=b,1:0;}
inline int read(){ int x=0,f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} int n,pat,x[M],y[M],z[M],vis[M];
double f[M][M],g[M][M],minv[M];
double l,r,mid,sum,tmp;
inline bool check(){
    memset(vis,0,sizeof(vis));
    sum=0,vis[1]=1;
    for(int i=1;i<=n;++i)
        minv[i]=g[1][i]-f[1][i]*mid;
    for(int i=2,j,k;i<=n;++i){
        for(tmp=1e16,k=-1,j=2;j<=n;++j)
            if(!vis[j]&&cmin(tmp,minv[j])) k=j;
        if(k==-1) break;
        for(vis[k]=1,sum+=tmp,j=2;j<=n;++j) if(!vis[j])
            cmin(minv[j],g[k][j]-f[k][j]*mid);
    }
    return sum>=0;
}
inline double dis(int i,int j){
    static double dx,dy;
    dx=x[i]-x[j],dy=y[i]-y[j];
    return sqrt(dx*dx+dy*dy);
}
int main(){
    while(1){
        n=read();
        if(!n) return 0;
        for(int i=1;i<=n;++i){
            x[i]=read(),
            y[i]=read(),
            z[i]=read();
        }
        for(int i=1;i<=n;++i)
            for(int j=i+1;j<=n;++j)
                f[j][i]=f[i][j]=dis(i,j),
                g[j][i]=g[i][j]=abs(z[i]-z[j]);
        l=0,r=100;
        while(r-l>eps){
            mid=(l+r)/2,
            (check()?l:r)=mid;
        }
        printf("%.3lf\n",l);
    }
    return 0;
}

3.最优比率环

简而言之，最优比率环就是给你一个图，图上每条边有两个权值（固然，第二个权值可能恒为 1 ），而后让你在图中找到一个环，令环上边的个两个权值和的比值最大（或是最小）

而后就是要用上面生成树相似的思路，获得每条边边权为： \(a[i]-b[i]\times ans\) ，而后在图上找负环（有时是正环，可是正负环能够人工转化的嘛），找到就 \(check\) 成功

题目

[HNOI2009]最小圈

代码

//by Judge
#include<cstdio>
#include<cstring>
#include<iostream>
#define eps 1e-10
using namespace std;
const int M=1e5+7;
inline bool cmin(double& a,double b){return a>b?a=b,1:0;}
inline int read(){ int x=0,f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} int n,m,tim,pat,head[M],vis[M]; double l,r,mid,d[M];
struct Edge{ int to,next; double val; Edge(){}
    Edge(int a,double b,int c){ to=a,val=b,next=c; }
}e[M<<1];
inline void add(int u,int v,double z){
    e[++pat]=Edge(v,z,head[u]),head[u]=pat;
}
bool dfs(int u){ vis[u]=1;
    for(int i=head[u];i;i=e[i].next)
        if(cmin(d[e[i].to],d[u]+e[i].val-mid))
            if(vis[e[i].to]||dfs(e[i].to))
                return vis[u]=0,1;
    return vis[u]=0;
}
inline bool check(){
    memset(d,0,sizeof(d));
    for(int i=1;i<=n;++i)
        if(dfs(i)) return 1;
    return 0;
}
int main(){ n=read(),m=read(); double z;
    for(int i=1,x,y;i<=m;++i){
        x=read(),y=read();
        scanf("%lf",&z);
        add(x,y,z);
        if(z<0) l+=z; else r+=z;
    }
    while(r-l>eps){
        mid=(l+r)/2,
        (check()?r:l)=mid;
    }
    return !printf("%.8lf\n",l);
}

4.最大密度子图

题目通常就是给出 \(n\) 个点而后让你导出一个子图，让这个子图中边数与点数的比值最大

不过有时候题目中的边（甚至是点）可能会带上权值，那么具体问题具体分析就行了

题目

Hard Life

（话说这道题翻译照搬的\(bzoj\)，\(UVA\) 是要输出方案的啊...）

代码

（好像就这个代码最长了，题目讲得简单然鹅码量巨大）

//by Judge
#include<cmath>
#include<queue>
#include<cstdio>
#include<cstring>
#include<iostream>
#define db double
using namespace std;
const db eps=1e-8;const int N=1005;
inline int read(){ int x=0,f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} int n,m,S,T,pat,ans,h[N],q[N],head[N],cur[N];
int u[N],v[N],du[N],vis[N]; db U,l,r,mid;
struct Edge { int to,nxt; db flow; } e[N<<5];
void insert(int u,int v,db w) {
    e[++pat]=(Edge){v,head[u],w},head[u]=pat;
    e[++pat]=(Edge){u,head[v],0},head[v]=pat;
}
#define v e[i].to
inline bool bfs() { int hd=0,tl=1,u;
    memset(h,0,sizeof(h)),h[q[1]=S]=1;
    while(hd<tl) { u=q[++hd];
        for(int i=head[u]; i; i=e[i].nxt)
            if(e[i].flow>eps&&!h[v])
                h[v]=h[u]+1,q[++tl]=v;
    } return h[T];
}
db dfs(int x,db F) {
    if(x==T) return F; db w,used=0;
    for(int &i=cur[x]; i; i=e[i].nxt)
        if(h[x]+1==h[v]) {
            w=dfs(v,min(e[i].flow,F-used));
            e[i].flow-=w,e[i^1].flow+=w;
            used+=w; if(used==F) return F;
        }
    if(!used) h[x]=0;
    return used;
}
db dinic() { db res=0;
    for(;bfs();res+=dfs(S,1e16))
        memcpy(cur,head,sizeof(cur));
    return res;
}
void DFS(int x) { ++ans,vis[x]=1;
    for(int i=head[x]; i; i=e[i].nxt)
        if(e[i].flow>eps&&!vis[v]) DFS(v);
}
#undef v
inline bool check() {
    pat=1,S=0,T=n+1;
    memset(head,0,sizeof(head));
    for(int i=1; i<=n; ++i){
        insert(S,i,U),
        insert(i,T,U+mid+mid-du[i]);
    }
    for(int i=1; i<=m; ++i){
        insert(u[i],v[i],1),
        insert(v[i],u[i],1);
    }
    return U*n-dinic()>eps;
}
int main() {
    while(scanf("%d%d",&n,&m)!=EOF) {
        for(int i=0; i<=n; ++i) du[i]=0;
        memset(vis,0,sizeof(vis)),ans=-1,U=m;
        if(!m) {puts("1\n1\n\n");continue;}
        for(int i=1; i<=m; ++i)
            u[i]=read(),v[i]=read(),
            ++du[u[i]],++du[v[i]];
        for(l=0,r=m;r-l>1.0/n/n;){
            mid=(l+r)/2,
            (check()?l:r)=mid;
        }
        mid=l,check(),DFS(S);
        printf("%d\n",ans);
        for(int i=1; i<=n; ++i)
            if(vis[i]) printf("%d\n",i);
    } return 0;
}

5.最优比率流

仍是那句老话，和以前同样的，就是根据二分出来的答案赋边权，而后通常是每次建一边图跑费用流根据最小费用的正负性判断 \(ans\) 合法性（用 \(zkw\) 跑费用流有时候会挂，好比这个例题，固然多是我天生大常数那就 \(GG\) ）

题目

[SDOI2017]新生舞会

（这里是一道二分图最大费用最大流，可是代码实现中能够边权取反而后跑最小费用）

代码

//by Judge
#include<queue>
#include<cstdio>
#include<cstring>
#include<iostream>
#define eps 1e-8
#define ll long long
using namespace std;
const double inf=1e18+7;
const int M=205;
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
inline int Min(int a,int b){return a<b?a:b;}
inline bool cmax(double& a,double b){return a<b?a=b,1:0;}
char buf[1<<21],*p1=buf,*p2=buf;
inline int read(){ int x=0,f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} int n,S,T,pat=1,tot,head[M],cur[M]; ll a[M][M],b[M][M];
int inq[M],vis[M]; double l,r=1e4,mid,ans,dis[M]; queue<int> q;
struct Edge{ int to,flow,next; double val;
    Edge(int a,double b,int c,int d){
        to=a,val=b,flow=c,next=d;
    } Edge(){}
}e[M*M*M];
inline void add(int u,int v,double c,int w){
    e[++pat]=Edge(v,c,w,head[u]),head[u]=pat;
    e[++pat]=Edge(u,-c,0,head[v]),head[v]=pat;
}
#define v e[i].to
inline bool spfa(){
    for(int i=S;i<=T;++i)
        dis[i]=-inf,vis[i]=0;
    dis[S]=0,q.push(S);
    while(!q.empty()){
        int u=q.front();
        q.pop(),vis[u]=0;
        for(int i=head[u];i;i=e[i].next)
            if(e[i].flow&&cmax(dis[v],dis[u]+e[i].val))
                if(!vis[v]) vis[v]=1,q.push(v);
    } return dis[T]!=-inf;
}
int dfs(int u,int flow){
    if(u==T) return ans+=dis[T]*flow,flow;
    int res=0; vis[u]=1;
    for(int &i=cur[u],k;i;i=e[i].next)
        if(dis[v]==dis[u]+e[i].val)
            if(!vis[v]&&e[i].flow){
                k=dfs(v,Min(e[i].flow,flow-res));
                if(k){
                    res+=k;
                    e[i].flow-=k;
                    e[i^1].flow+=k;
                    if(res==flow)
                        break;
                }
            }
    return res;
}
inline bool check(){
    pat=1,ans=0,
    memset(head,0,sizeof(head));
    for(int i=1;i<=n;++i){
        add(S,i,0,1),
        add(n+i,T,0,1);
    }
    for(int i=1;i<=n;++i) 
        for(int j=1;j<=n;++j)
            add(i,n+j,a[i][j]-b[i][j]*mid,1);
    for(;spfa();dfs(S,1e9))
        memcpy(cur,head,sizeof(cur));
    return ans<=0;
}
int main(){
    n=read(),T=n+1<<1;
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
            a[i][j]=read();
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
            b[i][j]=read();
    for(tot=pat;r-l>eps;){
        mid=(l+r)/2,
        (check()?r:l)=mid;
    }
    return !printf("%.6lf\n",l);
}

6.各类其余结合

洛谷 blog 已添加

分数规划的一点感想

讲道理，其实分数规划这个东西哪儿都能套，然鹅如今关于分数规划的题目却并非不少，至少没有到泛滥的地步，不过相信在不（知道多）久的未来，分数规划会成为人尽皆知的一种算法

\[tHe\ EnD\]

哦，对了，参考文献？

01分数规划问题相关算法与题目讲解

poj 3155 (最大密度子图)

另外的另外，有时间的话会把本身出的分数规划的题目放到6这块内容上来

相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程