• 0/1分数规划?我不会啊!


    最近入门了 01 分数规划,于是写篇博客纪念

    (如果cnblogs上的格式您实在受不了,请转至我的洛谷blog,并且那里的讲解会稍微详细一些,另外,更新内容也在我的洛谷博客中)

    分数规划是一项不常用到的(但还蛮实用的)算法,一般来讲就是求一个最优比率。

    分数规划的简单介绍

    分数规划顾名思义就是求一个分数表达式的最大(小)值。

    比如说有 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)

    我们令 (sum=sum_{i=1}^{k} a[i] , tot=sum_{j=1}^{k} b[j])(k>=W)

    然后假设问题中的最优解为 (ans) ,那么必然有:

    [dfrac{sum}{tot}<=ans ]

    移项得:

    [sum<=ans imes tot ]

    继续移就得到:

    [sum-ans imes tot<=0 ]

    这样转化有什么用呢?那我们尝试将 (sum)(tot) 带回去,就可以得到这么一个式子:

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

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

    那么我们只需要二分这个 ans, 计算出每件物品的 (a-b imes 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
    ",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] imes ans) <=0 ]

    于是 (ans) 照常二分,对于条件成立的判断方式变一变就好了。

    其实你应该也猜到了判断方式,原先我们不是贪心取前面的 k 个数么,那这次我们把贪心改成背包就好了啊

    那样的话,对于每个物品来说,体积就是 b ,而价值就是上面的 (a[i]-b[i] imes 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
    ",calc());
    }
    

    2.最优比率生成树

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

    具体点说,题目一般会给你一张图(有时候是完全图),然后给你 n 个点 m 条边,求出原图包含 n 个点的一棵树,使得所选的边的两个权值和比值最大

    那么求法呢,其实也很类似:

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

    那么对于每一条边,我们让他的边权为 (a[i]-b[i] imes 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
    ",l);
    	}
    	return 0;
    }
    

    3.最优比率环

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

    然后就是要用上面生成树类似的思路,得到每条边边权为: (a[i]-b[i] imes 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
    ",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
    1
    
    ");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
    ",ans);
            for(int i=1; i<=n; ++i)
                if(vis[i]) printf("%d
    ",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
    ",l);
    }
    

    6.各种其他结合

    洛谷 blog 已添加

    分数规划的一点感想

    讲道理,其实分数规划这个东西哪儿都能套,然鹅现在关于分数规划的题目却并不是很多,至少没有到泛滥的地步,不过相信在不(知道多)久的将来,分数规划会成为人尽皆知的一种算法

    [tHe EnD ]

    哦,对了,参考文献?

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

    poj 3155 (最大密度子图)

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

  • 相关阅读:
    Verilog手绘FVH信号
    Verilog编码规范与时序收敛
    关于DDS的基础知识
    阅读ug949-vivado-design-methodology笔记
    在windows系统上使用pip命令安装python的第三方库
    pandas第一课
    视频外同步信号研究---fvh
    FPGA调试技巧
    关于FIFO异步复位的问题
    搭建一个microblaze的最小系统
  • 原文地址:https://www.cnblogs.com/Judge/p/10173795.html
Copyright © 2020-2023  润新知