最大密度子图

最大密度子图

一、了解概念

首先要了解什么是最大密度子图，顾名思义，就是所有子图中密度最大的那一个。那么什么是一个子图的密度呢，这里规定子图的密度就是子图中边数和点数的比值。

二、怎么做

大体上是一个二分的01分数规划，但二分的判断函数要通过最小割模型来求，这里的求法是问题的关键，具体有如下两种求法。

1、方法一

我们选一个子图时，规定说选取原图的某条边时，那么这条边的两个顶点也一定要选，而选一个单独的点是可以选的。也就是说不能只选某条边而不选该边两个顶点，这不是一个子图的定义。在这个限制条件下，我们要尽可能的多选边，少选点，那么我们可以将边看作点，将点也看做点，从边变成的点向其两个端点做一条有向边，我们选这条边就必须要选这两个端点，这不就是闭合子图嘛，那么求最大密度子图的问题就变成了在新图中求最大权闭合子图的问题了（不懂的可以了解一下）
建图时我们将边变成的点的点权设为1，点变成的点的点权设为-g，这样最大权闭合子图求出来的答案就是我们要求的$Ne-g\times Nv$的最大值了。

2、方法二

但方法一需要建的点和边太多了，不方便，所以我们可以利用已有的性质来找到一个更好的方法，即方法二。
我们已知图中点数是n，边数是m，我们从源点s向每一个点建一条容量为m的边，原图内部的边容量都建成1，再从每个点向汇点t建一条容量为2g-d_v+m的边。我们求出该流网络中的最小割，然后我们要求的$Ne-g\times Nv$的最大值就是$\frac{n\times m-最小割}{2}$。

证明如下，可选择性的看，看不懂的话知道做法后可跳过。。。

三、证明

首先我们先写出密度的计算式：$\frac{Ne}{Nv}$，其中Ne是边的个数，Nv是点的个数。我们要求这个式子的最大值，见到分式和最大化，我们就想到了01分数规划。我们假设最大密度为g，那么就有$\frac{Ne}{Nv}\le g$，通过化简得到$Ne-g\times Nv\le 0$，根据01分数规划，我们不断二分g，求$Ne-g\times Nv$的最大值，使其最大值等于0，则这个g就是我们要求的答案了。

那么问题的关键就到了如何 计算$Ne-g\times Nv$的最大值？我们通过变形，原式的最大值就等于$g\times Nv-Ne$的最小值。所以我们要求的就是一个子图的点数和边数，对于确定的一个子图，其点数很好算，那么子图中的边数怎么算呢？我们可以将每一条边分成两半来看，一半给左端点，一半给右端点，那么子图中的边数就等于各点的度数之和$÷2$ 减去 该点集和点集外的边数$÷2$ ，这里的点集内和点集外的边数就等于这里的割集。
我们用公式来表示就是$g\times Nv-Ne=g\times Nv-(\frac{{\sum }_{v\in v_1}{d}_v}{2}-\frac{C[v_1,v_2]}{2})$。这里的$v_1$是选取子图的点集、$v_2$是$v_1$的补集、$d_v$是v点的度数。然后我们通过化简就可以得到$\frac{1}{2}({\sum }_{v\in v_1}(2g-d_v)+C[v_1,v_2])$。我们就从每个点到T建一条容量为$2g-d_v$的边来限制点权，为了使其是个正数，我们给他再加个大数m（边数[这里随便加一个大数最终是正数就行，我们选择用m]），保证是个正数，通过证明加上m后不会影响最小割的正确性。

然后我们我们要求按照方法二的建图，其最小割计算结果和我们要求的式子$Ne-g\times Nv$的最大值之间的关系。我们将原图分成两个集合S集合和T集合，所有的割集是S和T之间的边，S集合包含源点s和子图$v_1$，T集合包含汇点t和$v_1$的补集$v_2$。那么所有割集的情况如下图，共四种：

因为没有s到t的边，所以第一种情况不存在。那么最小割的容量$C[S,T]={\sum }_{v\in v_2}m+{\sum }_{v\in v_2}(2g-d_v)+{\sum }_{u\in v_1}{\sum }_{v\in v_2}{C}_{(u,v)}$，化简后就是$C[S,T]=m\times n+2g\times Nv-2Ne$，那么$Ne-g\times Nv=\frac{n\times m-C[S,T]}{2}$，其他都是定值，因此最小割就是原式的最大值。这样我们也自动就确定了子图的点集就是S集合减去源点s的集合。

经典例题

题目链接

传送门：POJ3155：Hard Life

思路

这就是一个裸的求最大密度子图的题，直接套方法二。这里要求输出子图的点集$v_1$，最后找到答案后，从残留网络中的源点s出发，沿着流量为正的边走，能走到的点就是S集合中的点，即$v_1$，输出即可。

代码

#include<iostream>
#include<cstring>
#include<queue>
using namespace std;
const int N = 110,M = (1000+2*N)*2,INF = 0x3f3f3f3f;
const double eps = 1e-8;

int n,m,s,t,dis[N],cur[N];	//dinic
int e[M],ne[M],h[N],idx;	//链式前向星存图
double w[M];			//容量存成double型
int res,dg[N];			//dg存点的度数
bool vis[N];
pair<int,int> es[M];

void add(int a,int b,double c)	//建图模板
{
	e[idx] = b;
	w[idx] = c;
	ne[idx] = h[a];
	h[a] = idx++;
	e[idx] = a;
	w[idx] = c;
	ne[idx] = h[b];
	h[b] = idx++;
}

void build(double g)	//建图函数，每次二分到一个答案时要重新建图
{
	memset(h,-1,sizeof h);
	idx = 0;
	for(int i = 0;i < m;i++)
		add(es[i].first,es[i].second,1);
	for(int i = 1;i <= n;i++)
	{
		add(s,i,m);
		add(i,t,m+g*2-dg[i]);
	}
}

bool bfs()		//dinic模板
{
	queue<int> q;
	memset(dis,-1,sizeof dis);
	q.push(s);
	dis[s] = 0;
	cur[s] = h[s];
	while(q.size())
	{
		int u = q.front();
		q.pop();
		for(int i = h[u];~i;i = ne[i])
		{
			int v = e[i];
			if(dis[v] == -1 && w[i])
			{
				dis[v] = dis[u] + 1;
				cur[v] = h[v];
				if(v == t)
					return 1;
				q.push(v);
			}
		}
	}
	return 0;
}

double dfs(int u,double limit)
{
	if(u == t)
		return limit;
	double flow = 0;
	for(int i = cur[u];~i;i = ne[i])
	{
		cur[u] = i;
		int v = e[i];
		if(dis[v] == dis[u]+1 && w[i])
		{
			double minf = dfs(v,min(w[i],limit-flow));
			w[i] -= minf;
			w[i^1] += minf;
			flow += minf;
			if(flow == limit)
				return flow;
		}
	}
	return flow;
}

double dinic()
{
	double ans = 0;
	while(bfs())
		ans += dfs(s,INF);
	return ans;
}

bool judge(double mid)		//二分的判断函数
{
	build(mid);
	double res = m*n-dinic();
	return res > 0;
}

void find(int u)		//一个dfs找方案
{
	vis[u] = 1;
	if(u != s)
		res++;
	for(int i = h[u];~i;i = ne[i])
	{
		int v = e[i];
		if(!vis[v] && w[i])
			find(v);
	}
}

int main()
{
	cin >> n >> m;
	s = 0,t = n+1;		//手动设置s和t
	for(int i = 0;i < m;i++)
	{
		int a,b;
		cin >> a >> b;
		dg[a]++,dg[b]++;	//a和b的度数+1
		es[i] = {a,b};		//暂时先将边存下
	}

	double l = 0,r = m;		//定义二分的左右端点
	while(r-l > eps)		//二分答案g
	{
		double mid = (r+l)/2;
		if(judge(mid))
			l = mid;
		else
			r = mid;
	}
	build(l);	//重新建一个图
	dinic();	//跑一边dinic
	find(s);	//从s出发找能到达的点
	if(res)
	{
		cout << res << endl;
		for(int i = 1;i <= n;i++)
			if(vis[i])
				cout << i << endl;
	}
	else
		cout << 1 << endl << 1 << endl;
 
	return 0;
}