保序迴歸問題
保序迴歸問題
基本形式
\[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\)，如果一個點變成了\(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\)
因為取值必須為整數，\(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;
}