ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

保序回归问题

2021-06-16 23:34:56  阅读:233  来源: 互联网

标签:ll 更优 int 回归 mid 问题 theta 保序 ql


保序回归问题

基本形式

\[f(x)=\sum_{i=1}^nw_i|a_i-b_i|^k \]

有一些要求形如

\[b_x\le b_y \]

最小化\(f(x)\)

一般解法

我们可以整体二分,对于值域区间\([l,r]\),我们二分\(mid=\frac{l+r}{2}\),对于当前需要考虑的变量,判断它的取值在\([l,mid]\)更优还是\((mid,r]\)更优。

有结论:

取值在\([l,mid]\)更优还是\((mid,r]\)更优等价于取\(mid\)还是\(mid+\theta\)更优,其中\(\theta\)是根据题目要求的精度而确定的极小值

所以我们现在判断一个点取\(mid\)更优还是\(mid+\theta\)更优即可

考虑把这个问题转化为最小权闭合子图问题。我们先钦定所有点选择\(mid\),如果一个点变成了\(mid+\theta\),那么所有有要求\(b_x\le b_y\)的点的\(y\)也必须选择\(mid+\theta\)。也就是每个点的点权为\(f(mid+\theta)-f(mid)\),求最小权闭合子图。

根据我们网络流的结果判断每个点取值在\([l,mid]\)更优还是\((mid,r]\)更优,分治解决即可。

例题

问题的形式是

\[f(x)=\sum_{i=1}^nw_i(a_i-b_i)^2 \]

我们只需要找到偏序关系即可。

对\(A,B\)分别建线性基。如果我们将一个元素\(x\in A\)从线性基中删除,而另一个元素\(y\not\in A\)可以插入,那么用\(y\)替换\(x\)可以使得集合大小不变,所以\(b_x\le b_y\)。同理,如果将一个元素\(x\in B\)从线性基中删除,而另一个元素\(y\not\in B\)可以插入,要求\(b_x\ge b_y\)。

而这样就一定能覆盖所有的情况吗?答案是肯定的。因为\(A,B\)已经是满足条件的最大的集合了,所以必须一一替换。而因为线性基可以用唯一的方案表示一个数,所以如果同时删去线性基中一个集合\(S\),可以被插入的所有数一定需要有\(S\)里的元素来表示,而\(S\)里的所有元素能够替换的数都已经连好了边,所以可以插入的所有数一定是\(S\)里每个元素能够替换的并集。

因为取值必须为整数,\(mid\) 取$ \left\lfloor\frac{l+r}{2}\right\rfloor$, \(\theta=1\)。然后就分治即可

#include <bits/stdc++.h>
using namespace std;
const int N=1005,M=N*N*4,inf=1e9;
#define ull unsigned long long
int n,m,v[N],a[N],b[N],mx;
bool visa[N],visb[N];
ull buc[64],c[N];
inline int insert(ull x){
	for(int i=63;~i;--i){
		if(!((x>>i)&1))continue;
		if(buc[i])x^=buc[i];
		else{buc[i]=x;return i;}
	}
	return -1;
}
namespace Flow{
	struct Edge{int to,next,flow;}e[M<<1];
	int head[N],ecnt=1,cur[N];
	inline void adde(int u,int v,int flow){
		e[++ecnt]=(Edge){v,head[u],flow};head[u]=ecnt;
		e[++ecnt]=(Edge){u,head[v],0};head[v]=ecnt;
	}
	struct Queue{
		int q[N],head,tail;
		inline void init(){head=1;tail=0;}
		inline void push(int x){q[++tail]=x;}
		inline int front(){return q[head];}
		inline void pop(){++head;}
		inline bool empty(){return head>tail;}
	}q;
	int s,t,dep[N];
	bool bfs(){
		memset(dep+1,-1,t<<2);
		dep[s]=0;q.init();q.push(s);cur[s]=head[s];
		while(!q.empty()){
			int u=q.front();q.pop();
			for(int i=head[u],v;i;i=e[i].next){
				if(e[i].flow<=0)continue;v=e[i].to;
				if(dep[v]==-1){dep[v]=dep[u]+1;cur[v]=head[v];q.push(v);}
			}
		}
		return dep[t]>=0;
	}
	int dinic(int u,int flow){
		if(u==t)return flow;
		int ret=0,f;
		for(int i=cur[u],v;i;i=e[i].next){
			v=e[i].to;cur[u]=i;
			if(e[i].flow<=0||dep[v]!=dep[u]+1)continue;
			f=dinic(v,min(flow,e[i].flow));
			ret+=f;flow-=f;e[i].flow-=f;e[i^1].flow+=f;
			if(flow<=0)break;
		}
		return ret;
	}
	inline int Maxflow(int S,int T){
		s=S;t=T;
		int ans=0;
		while(bfs())ans+=dinic(s,inf);
		return ans;
	}
	inline void init(){
		memset(head,0,sizeof(head));
		ecnt=1;
	}
}
vector<int> to[N];
#define pb push_back
int f[N],q[N];
#define ll long long
inline ll calc(int x,int y){
	return 1ll*(y-v[x])*(y-v[x]);
}
int id[N],q1[N],q2[N],vis[N],inde;
void work(int l,int r,int ql,int qr){
	if(l==r){
		for(int i=ql;i<=qr;++i)f[q[i]]=l;
		return;
	}
	if(ql>qr)return;
	int mid=(l+r)>>1;++inde;
	for(int i=ql;i<=qr;++i)id[q[i]]=i-ql+1,vis[q[i]]=inde;
	int s=qr-ql+2,t=s+1;
	Flow::init();
	for(int i=ql,u;i<=qr;++i){
		u=q[i];
		for(int v:to[u])if(vis[v]==inde)
			Flow::adde(id[u],id[v],inf);
		ll x=calc(u,mid)-calc(u,mid+1);
		if(x>0)Flow::adde(s,id[u],x);
		else Flow::adde(id[u],t,-x);
	}
	Flow::Maxflow(s,t);
	int tot1=0,tot2=0;
	for(int i=ql;i<=qr;++i)
		if(Flow::dep[id[q[i]]]==-1)q1[++tot1]=q[i];
		else q2[++tot2]=q[i];
	memcpy(q+ql,q1+1,tot1<<2);
	memcpy(q+ql+tot1,q2+1,tot2<<2);
	work(l,mid,ql,ql+tot1-1);
	work(mid+1,r,ql+tot1,qr);
}
int main(){
	scanf("%d %d",&n,&m);
	for(int i=1;i<=n;++i)scanf("%llu",&c[i]);
	for(int i=1;i<=n;++i)scanf("%d",&v[i]);
	for(int i=1;i<=m;++i)scanf("%d",&a[i]),visa[a[i]]=1;
	for(int i=1;i<=m;++i)scanf("%d",&b[i]),visb[b[i]]=1;
	for(int i=1;i<=m;++i){
		memset(buc,0,sizeof(buc));
		for(int j=1;j<=m;++j)if(i!=j)insert(c[a[j]]);
		for(int j=1,p;j<=n;++j){
			if(visa[j])continue;
			p=insert(c[j]);
			if(p>=0)buc[p]=0,to[a[i]].pb(j);
		}
	}
	for(int i=1;i<=m;++i){
		memset(buc,0,sizeof(buc));
		for(int j=1;j<=m;++j)if(i!=j)insert(c[b[j]]);
		for(int j=1,p;j<=n;++j){
			if(visb[j])continue;
			p=insert(c[j]);
			if(p>=0)buc[p]=0,to[j].pb(b[i]);
		}
	}
	for(int i=1;i<=n;++i)q[i]=i,mx=max(mx,v[i]);
	work(0,mx,1,n);
	ll ans=0;
	for(int i=1;i<=n;++i)ans+=calc(i,f[i]);
	printf("%lld\n",ans);
	return 0;
}

如果不考虑图的特殊性(基环树)的话,直接用保序回归的一般解法即可。

但是图具有特殊性,我们考虑使用更高效的DP做法替换最小权闭合子图的网络流做法。

也就是说,我们需要在线性复杂度内,确定每个点选\(mid\)还是\(mid+1\)。

  • 如果图是一棵树

设\(f_{u,0/1}\)表示该点选择\(mid/mid+1\)的最小代价。

如果存在限制\(b_u\le b_v\),\(f_{u,0}\leftarrow f_{v,0}+f_{v,1},\ \ f_{u,1}\leftarrow f_{v,1}\)

否则\(b_u\ge b_v\),\(f_{u,1}\leftarrow f_{v,0}+f_{v,1},\ \ f_{u,0}\leftarrow f_{v,0}\)

随便选一个点为根DFS,确定了根的颜色之后染色。如果\(b_{fa}\le b_{u}\)且\(b_{fa}\)取了\(mid+1\),那么\(u\)必须取\(mid+1\),如果\(b_{fa}\ge b_{u}\)且\(b_{fa}\)取了\(mid\),那么\(u\)必须取\(mid\)。其它情况根据DP值取。

  • 图是一个基环树

设\(f_{u,0/1,0/1}\)表示根选择\(mid/mid+1\),该点选择\(mid/mid+1\)的最小代价。其中根为环上任意一点

转移与上面类似。\(f_{?,0,?}\)和\(f_{?,1,?}\)分别转移

但是因为基环树多一条边的限制,假设这一条环边是\((rt,u)\),那么需要根据\(b_{rt}\)和\(b_{u}\)的大小关系,将\(f_{u,0,1}\)或\(f_{u,1,0}\)设置为\(inf\),即该状态不合法。对于\(rt\),只有\(f_{rt,0,0}\)和\(f_{rt,1,1}\)合法。

染色和上面类似。

分治+树形DP即可

const int N=2e5+5;
#define ll long long
const ll inf=2e18;
inline void Min(ll &a,ll b){if(a>b)a=b;}
struct Edge{int to,next;}e[N<<1];
int head[N],ecnt=1;
inline void adde(int u,int v){e[++ecnt]=(Edge){v,head[u]};head[u]=ecnt;}
int inq[N],inde;
int n,m,k,a[N],root,dep[N];
void dfs1(int u,int pre){
	for(int i=head[u],v;i&&!root;i=e[i].next){
		v=e[i].to;if((i^1)==pre||inq[v]!=inde)continue;
		if(!dep[v])	dep[v]=dep[u]+1,dfs1(v,i);
		else if(dep[v]<dep[u])root=u;
	}
}
inline int find_rt(int x){
	if(m==n-1)return x;
	else{root=0;dep[x]=1;dfs1(x,0);return root?root:x;}
}
inline ll update(ll x){return min(x,inf);}
ll f[N][2][2],ans[N],w[N];
inline ll Pow(ll x){return k==1?x:x*x;}
inline ll calc(int x,int y){return w[x]*Pow(abs(a[x]-y));}
bool vis[N];
void dfs(int u,int rt,int pre,int mid){
	int flag=0;
	f[u][0][0]=f[u][1][0]=calc(u,mid);
	f[u][0][1]=f[u][1][1]=calc(u,mid+1);
	vis[u]=1;
	for(int i=head[u],v;i;i=e[i].next){
		v=e[i].to;if(inq[v]!=inde||pre==(i^1))continue;
		if(v==rt)	flag=i;
		else if(!vis[v]){
			dfs(v,rt,i,mid);
			bool cur=(i&1)^1;
			for(int j=0;j<2;++j){
				f[u][j][cur]=update(f[u][j][cur]+f[v][j][cur]);
				f[u][j][cur^1]=update(f[u][j][cur^1]+min(f[v][j][cur],f[v][j][cur^1]));
			}
		}
	}
	if(flag){
		if(flag&1)f[u][1][0]=inf;
		else f[u][0][1]=inf;
	}
	if(u==rt)	for(int i=0;i<2;++i)f[u][i][i^1]=inf;
}
bool col[N],coled[N];
void dfs2(int u,int pre,bool colu,bool c){
	col[u]=c;coled[u]=1;
	for(int i=head[u],v;i;i=e[i].next){
		v=e[i].to;if((i^1)==pre||coled[v])continue;
		if((i&1)!=c)dfs2(v,i,colu,c);
		else dfs2(v,i,colu,f[v][colu][0]>f[v][colu][1]);
	}
}
int q[N],q1[N],q2[N];
void work(int l,int r,int ql,int qr){
	if(ql>qr)return;
	if(l==r){for(int i=ql;i<=qr;++i)ans[q[i]]=l;return;}
	int mid=(l+r)>>1;++inde;
	for(int i=ql;i<=qr;++i)inq[q[i]]=inde,vis[q[i]]=0,coled[q[i]]=0,dep[q[i]]=0;
	for(int i=ql;i<=qr;++i)if(!vis[q[i]]){
		int rt=find_rt(q[i]);
		dfs(rt,rt,0,mid);
		if(f[rt][0][0]<f[rt][1][1])dfs2(rt,0,0,0);
		else dfs2(rt,0,1,1);
	}
	int tot1=0,tot2=0;
	for(int i=ql;i<=qr;++i)
		if(!col[q[i]])q1[++tot1]=q[i];
		else q2[++tot2]=q[i];
	memcpy(q+ql,q1+1,tot1<<2);	memcpy(q+ql+tot1,q2+1,tot2<<2);
	work(l,mid,ql,ql+tot1-1);	work(mid+1,r,ql+tot1,qr);
}
int main(){
	n=read();m=read();k=read();
	int mx=0;
	for(int i=1;i<=n;++i)a[i]=read(),mx=max(mx,a[i]);
	for(int i=1;i<=n;++i)w[i]=read();
	for(int i=1,u,v;i<=m;++i){
		u=read();v=read();bool val=readtype();
		if(val)adde(u,v),adde(v,u);
		else adde(v,u),adde(u,v);
	}
	for(int i=1;i<=n;++i)q[i]=i;
	work(0,mx,1,n);
	ll x=0;
	for(int i=1;i<=n;++i)x+=calc(i,ans[i]);
	printf("%lld\n",x);
	return 0;
}

标签:ll,更优,int,回归,mid,问题,theta,保序,ql
来源: https://www.cnblogs.com/harryzhr/p/14891588.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有