基于BGL的社团结构检测

一、功能描述

通过在控制台界面输入已有数据文件的名称，读取数据文件内容并建立图结构，进而利用fastQ算法对所得图进行社团检测，得到并输出最优划分方式。

二、设计

1. 数据类型的设计和说明

数据类型为一个结点无权，但边可赋权值并允许自回路的无向图：

① 边赋权：随着社团的合并，社团数量减少，社团间会出现大量平行边，在BGL里邻接矩阵的存储方式不支持平行边，所以改而采用边的权重增加的方式体现平行边的存在。

② 自回路：随着合并，社团内部的结点之间会有边相连，为了体现这种边的存在，采用了自回路等效表示的方法。同样，也是通过给自回路赋权来体现平行边的概念。

③ 无向：在该程序中，fastQ算法所讨论的图均无方向。如果数据文件出现重复数据，程序运行过程中会自动忽略。

④ 结点无权：因为是无向图并且已经用自回路赋权的方式表示社团内部结点之间的边，所以结点无需赋权。

2. 关键算法的设计

算法分为以下几步：

⑴ 准备工作：

①  浏览文件内容，得到结点数目n和边的数目m；

②  建立具有n个结点的图；

③  读取文件内容，并同时为图逐条添加边（包括边权值初始化）；

⑵ 图操作：

① 穷举所有社团两两合并的情况，计算使得ΔQ最大的合并方式；

② 记录最佳合并方式（在数组order内），对①所得的ΔQ最大值做累计；

③ 根据最佳合并方式更新图，对图的边进行修改；

④ 重复第①步，直到所有社团合并为一个社团；

⑶ 输出划分结果：

① 找到使得ΔQ的累加值最大的位置；

② 根据①中所得最大位置，从order中读取相应数量的数据，对社团编号进行归组；

 

3. 模块结构图及各模块的功能

关键在于图的合并过程：

for(jNghbrVrtx=adjacent_vertices(*j,ug);jNghbrVrtx.first!=jNghbrVrtx.second;++jNghbrVrtx.first){i2jNghbrEdge=edge(*i,*jNghbrVrtx.first,ug);double i2jNghbrEdgeWeight=get(edge_weight_t(),ug,i2jNghbrEdge.first);j2jNghbrEdge=edge(*j,*jNghbrVrtx.first,ug);double j2jNghbrEdgeWeight=get(edge_weight_t(),ug,j2jNghbrEdge.first);remove_edge(*i,*jNghbrVrtx.first,ug);if(*jNghbrVrtx.first!=*i){remove_edge(*i,*jNghbrVrtx.first,ug);add_edge(*i,*jNghbrVrtx.first,i2jNghbrEdgeWeight+j2jNghbrEdgeWeight,ug);}else{j2jEdge=edge(*j,*j,ug);double j2jEdgeWeight=get(edge_weight_t(),ug,j2jEdge.first);remove_edge(*i,*jNghbrVrtx.first,ug);add_edge(*i,*i,i2jNghbrEdgeWeight+j2jNghbrEdgeWeight+j2jEdgeWeight,ug);}}

下面就把代码贴出来，主要是一个main，三个头文件：

main.cpp

#include"getDeltaQ.h"#include"combineIJ.h"#include"tip.h"#include<vector>#include<string>#include<math.h>#include<fstream>#include<iostream>#include<boost\graph\graph_utility.hpp>//为了调用 print_vertices()和print_edges()#include<boost\graph\adjacency_matrix.hpp>using namespace std;using namespace boost;//边的权重属性类型---边有权，权重为double类型typedef property<edge_weight_t,double> edgeWeightProperty;//无向图中的结点类型:无向图、结点无权、边有权，图无权【是否允许平行边的概念只在邻接表里面有，邻接矩阵忽略此参数】typedef adjacency_matrix<undirectedS,no_property,edgeWeightProperty> UGraph;//边迭代器类型typedef graph_traits<UGraph>::edge_iterator edge_iter;//结点迭代器类型typedef graph_traits<UGraph>::vertex_iterator vertex_iter;//边的描述符类型typedef graph_traits<UGraph>::edge_descriptor edge_dscr;//入边迭代器类型typedef graph_traits<UGraph>::in_edge_iterator  in_edge_iter;int main(int,char*[]){/****************************************************************************///****************************************//*                                      *//*       构       造       图        *//*                                      *//****************************************/***************************************************************************/tip();int m=0,n=0;ifstream fBaseValue;string fName;do{cout<<"请输入文件名：";cin>>fName;fBaseValue.open(fName,ios::binary);} while (!fBaseValue);//得到m和n的值int node;do{fBaseValue>>node;m++;if(node>n){n=node;}fBaseValue.seekg(sizeof(char),ios::cur);} while (node!=0);m/=2;edgeWeightProperty eWP=1./m;UGraph ug(n);fBaseValue.seekg(0,ios::beg);//开始给图添加边int i,j;while(1){fBaseValue>>i;fBaseValue.seekg(sizeof(char),ios::cur);fBaseValue>>j;fBaseValue.seekg(sizeof(char),ios::cur);if(i==0)break;add_edge(i-1,j-1,eWP,ug);}fBaseValue.close();/****************************************************************************///****************************************//*                                      *//*     检 查 图 是 否 正 确 建 立        *//*                                      *//****************************************/***************************************************************************///检查边设置是否正确---正确//print_vertices(ug,"123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ");//cout<<endl;//print_edges(ug,"123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ");//cout<<endl<<endl;//print_graph(ug,"123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ");//检查边权重设置是否正确---正确//double sum=0;//property_map<UGraph,edge_weight_t>::type EdgeWeightMap=get(edge_weight_t(),ug);//std::pair<edge_iter,edge_iter> edgePair;//for(edgePair=edges(ug);edgePair.first!=edgePair.second;++edgePair.first)//{////std::cout<<EdgeWeightMap[*edgePair.first]<<" ";//这样写也是可以输出正确结果的，与下面的写法的区别：尚不明确//cout<<get(EdgeWeightMap,*edgePair.first);//sum+=get(EdgeWeightMap,*edgePair.first);//}//cout<<"sum="<<sum<<endl;//输出所有的结点的度数---正确//int sum=0;//std::pair<vertex_iter,vertex_iter> vertexPair;//for(vertexPair=vertices(ug);vertexPair.first!=vertexPair.second;++vertexPair.first)//{//cout<<in_degree(*vertexPair.first,ug)<<endl;//sum+=in_degree(*vertexPair.first,ug);//}//cout<<"sum="<<sum<<endl;//sum=156//说明：如果只是计算deltaQ的话，结点不需要具有权重属性!/****************************************************************************///****************************************//*                                      *//*    图 操 作--主 算 法 实 现 部 分      *//*                                      *//****************************************/***************************************************************************/vector<double> deltaQ(n,-10);//记录ΔQ的变化情况，为了防止迭代器越界，故意多申请一个位置，用到的是n-1个位置vector<double>::iterator pDeltaQ=deltaQ.begin();//deltaQ的迭代器vector<double>::iterator maxP_DeltaQ=deltaQ.begin();//指向deltaQ积分最大的位置vector<vector<int>> tree(n,vector<int>(n,0));vector<vector<int>::iterator> pTree(n);//迭代器数组for(int i=0;i<n;i++){tree[i][0]=i+1;pTree[i]=++tree[i].begin();}//检查树tree的设定是否正确---正确//cout<<endl<<endl;//for(int i=0;i<n;i++)//{//for(int j=0;j<n;j++)//{//cout<<tree[i][j]<<" ";//}//cout<<endl;//}//采用Cn2的方式进行合并，并且算各种情况下的ΔQ和Qstd::pair<vertex_iter,vertex_iter> pVertexI,pVertexJ;pVertexJ=vertices(ug);//如果不给它初始化，也不会报告，结果导致了后面循环里的内层循环出现死循环。吃一堑长一智啊！vertex_iter maxI,maxJ;int nNum=n;//n的替代品double sumDeltaQ=0;double tempSumDeltaQ=0;vector<vector<int>> order(n-1,vector<int>(2));int maxPosition=0;while(nNum>1){for(pVertexI=vertices(ug);pVertexI.second-pVertexI.first>=2;++pVertexI.first){for(pVertexJ.first=pVertexI.first+1;pVertexJ.first!=pVertexJ.second;++pVertexJ.first){if((in_degree(*pVertexI.first,ug)==0)||(in_degree(*pVertexJ.first,ug)==0))//这里必须与0做比较，返回值为0的时候并不等价于逻辑假{continue;}double tempDeltaQ=getDeltaQ(ug,pVertexI.first,pVertexJ.first);if(tempDeltaQ>*pDeltaQ)//打擂台,得到最优的i，j并且记录此时的ΔQ{*pDeltaQ=tempDeltaQ;maxI=pVertexI.first;maxJ=pVertexJ.first;}}}//为下一次循环做准备：修改ug，体现抱团关系，孤立被合并的社团。combineIJ(ug,maxI,maxJ);order[n-nNum][0]=*maxI;order[n-nNum][1]=*maxJ;//记录deltaQ积分最大的地方，等价于记录Q值最大的地方。sumDeltaQ+=*pDeltaQ;if(sumDeltaQ>tempSumDeltaQ){tempSumDeltaQ=sumDeltaQ;maxPosition++;}pDeltaQ++;nNum--;}for(int i=0;i<maxPosition;i++){vector<int>::iterator maxJtempPTree=tree[order[i][1]].begin();while (maxJtempPTree<pTree[order[i][1]])//指向的位置是一个待写的位置，当前值为0，不应复制到i位置上去{*pTree[order[i][0]]=*maxJtempPTree;*maxJtempPTree=0;maxJtempPTree++;pTree[order[i][0]]++;}}int groupNum=1;for(int i=0;i<n;i++){if(tree[i][0]==0)continue;cout<<"第"<<groupNum<<"组： ";for(int j=0;j<n;j++){if(tree[i][j]!=0)cout<<tree[i][j]<<" ";}groupNum++;cout<<endl;}return 0;}

tip.h

#include<iostream>#include<string>using std::cout;using std::cin;using std::endl;using std::string;void tip(){    cout<<endl<<"\t\t★---您好，欢迎进入fastQ算法计算及演示系统！---★"<<endl;cout<<endl<<endl<<endl;cout<<"\t\t╭----------------------------------------------╮"<<endl;cout<<"\t\t｜　　     注       意       事       项 　      ｜"<<endl;cout<<"\t\t｜-----------------------------------------------｜"<<endl;cout<<"\t\t｜　　　　1.数据文件是二进制文件                 ｜"<<endl;cout<<"\t\t｜　　　　　　　　　　　　　　　　　　　      　 ｜"<<endl;cout<<"\t\t｜　　　　2.数据文件的形式是序偶形式　　　 　    ｜"<<endl;cout<<"\t\t｜　　　　　　　　　　　　　　　　　　　      　 ｜"<<endl;cout<<"\t\t｜　　　　3.数据文件放在演示文件所在目录下　     ｜"<<endl;cout<<"\t\t｜　　　　　　　　　　　　　　　　　　　      　 ｜"<<endl;cout<<"\t\t｜　　　　4.序偶对之间以长为1位的任意字符间隔    ｜"<<endl;cout<<"\t\t｜　　　　　　　　　　　　　　　　　　　      　 ｜"<<endl;cout<<"\t\t｜　　　　5.在接下来要输入的文件名需包含扩展名   ｜"<<endl;cout<<"\t\t╰-----------------------------------------------╯"<<endl;}

getDeltaQ.h

#include<iostream>#include<boost\graph\adjacency_matrix.hpp>using namespace std;using namespace boost;//边的权重属性类型---边有权，权重为double类型typedef property<edge_weight_t, double> edgeWeightProperty;//无向图中的结点类型:无向图、结点无权、边有权，图无权【是否允许平行边的概念只在邻接表里面有，邻接矩阵忽略此参数】typedef adjacency_matrix<undirectedS,no_property,edgeWeightProperty> UGraph;//边迭代器类型typedef graph_traits<UGraph>::edge_iterator edge_iter;//结点迭代器类型typedef graph_traits<UGraph>::vertex_iterator vertex_iter;//边的描述符类型typedef graph_traits<UGraph>::edge_descriptor edge_dscr;double getDeltaQ(const UGraph &ug,vertex_iter i,vertex_iter j){double deltaQ=0;//入边迭代器类型typedef graph_traits<UGraph>::in_edge_iterator  in_edge_iter;//获取i，j的所有邻接边，判断ij之间是否有边。//deltaQ=2*(ij之间边的权重[有一条边，或者没有边]-（i的所有邻边权重之和[包括ij之间的边]*j的邻边权重之和[包括ij之间的边]）)//计算得到Ai,Aj,Eij的值double Ai=0,Aj=0,Eij=0;double temp;std::pair<edge_dscr,bool> i2jEdge;std::pair<in_edge_iter,in_edge_iter> iNghbrEdge,jNghbrEdge;//分别迭代i，j的邻接边for(iNghbrEdge=in_edges(*i,ug);iNghbrEdge.first!=iNghbrEdge.second;iNghbrEdge.first++){temp=get(edge_weight_t(),ug,*iNghbrEdge.first);if(target(*iNghbrEdge.first,ug)!=source(*iNghbrEdge.first,ug))//如果不是从自身到自身的边，要减半{temp/=2.0;}Ai+=temp;}for(jNghbrEdge=in_edges(*j,ug);jNghbrEdge.first!=jNghbrEdge.second;jNghbrEdge.first++){temp=get(edge_weight_t(),ug,*jNghbrEdge.first);if(target(*jNghbrEdge.first,ug)!=source(*jNghbrEdge.first,ug)){temp/=2.0;}Aj+=temp;}i2jEdge=edge(*i,*j,ug);Eij=get(edge_weight_t(),ug,i2jEdge.first)/2.;//计算deltaQ的值deltaQ=2*(Eij-Ai*Aj);return deltaQ;}

combineIJ.h

#include<iostream>#include<stdlib.h>#include<boost\graph\adjacency_matrix.hpp>#include<boost\graph\graph_traits.hpp>using namespace std;using namespace boost;void combineIJ(UGraph &ug,vertex_iter i,vertex_iter j){//边的权重属性类型---边有权，权重为double类型typedef property<edge_weight_t, double> edgeWeightProperty;//无向图中的结点类型:无向图、结点无权、边有权，图无权【是否允许平行边的概念只在邻接表里面有，邻接矩阵忽略此参数】typedef adjacency_matrix<undirectedS,no_property,edgeWeightProperty> UGraph;//边迭代器类型typedef graph_traits<UGraph>::edge_iterator edge_iter;//结点迭代器类型typedef graph_traits<UGraph>::vertex_iterator vertex_iter;//边的描述符类型typedef graph_traits<UGraph>::edge_descriptor edge_dscr;//入边迭代器类型typedef graph_traits<UGraph>::in_edge_iterator  in_edge_iter;std::pair<edge_dscr, bool> i2jNghbrEdge,j2jNghbrEdge;std::pair<edge_dscr, bool> j2jEdge;typedef graph_traits<UGraph>::adjacency_iterator adjacency_iter;std::pair<adjacency_iter, adjacency_iter> jNghbrVrtx;//迭代与第j个结点相邻的结点for(jNghbrVrtx=adjacent_vertices(*j,ug);jNghbrVrtx.first!=jNghbrVrtx.second;++jNghbrVrtx.first){i2jNghbrEdge=edge(*i,*jNghbrVrtx.first,ug);double i2jNghbrEdgeWeight=get(edge_weight_t(),ug,i2jNghbrEdge.first);j2jNghbrEdge=edge(*j,*jNghbrVrtx.first,ug);double j2jNghbrEdgeWeight=get(edge_weight_t(),ug,j2jNghbrEdge.first);remove_edge(*i,*jNghbrVrtx.first,ug);if(*jNghbrVrtx.first!=*i){remove_edge(*i,*jNghbrVrtx.first,ug);add_edge(*i,*jNghbrVrtx.first,i2jNghbrEdgeWeight+j2jNghbrEdgeWeight,ug);}else{j2jEdge=edge(*j,*j,ug);double j2jEdgeWeight=get(edge_weight_t(),ug,j2jEdge.first);remove_edge(*i,*jNghbrVrtx.first,ug);add_edge(*i,*i,i2jNghbrEdgeWeight+j2jNghbrEdgeWeight+j2jEdgeWeight,ug);}}//删除顶点j的所有出边和入边。该顶点则仍然保留在图的顶点集中。clear_vertex(*j,ug);}

有时间应该多看看侯捷的《STL源码剖析》这本书

另外顺便把刚开始的时候用MATLAB实现的代码也贴上来：

clear;clc%%%%%%%%%%%%n、m、baseValue的值都是一旦确定就不再改变的%%%%%%%%%%%%%%%%%%%%%%%%计算得到E和m值，另外给n赋初始值%%%%%%%%%%%%load baseValuen=size(baseValue,1);  %n是图中结点的数目，之后它的复制品nNum会在代码里逐一递减。m=sum(sum(baseValue))./2;   %m是初始时刻图中边的数目，经过计算后保证值不变。%记录每次输出的被合并的结点的编号，即maxU和maxV，每次输出的两个数在同一行。order=zeros(n-1,2);%记录使得△Q最大的一种合并方式是把编号为几的两个结点合并。maxI=0;maxJ=0;deltaQ=-100.*ones(1,n-1);%deltaQ用来记录每次合并的最大△Q值，从这里面找出最后一个正数△Q值就知道应该从哪儿切割树了。Qmax=zeros(1,n);%记录34个Q值，Q值最大的时候就最应该在哪儿切割树。E=baseValue./(2*m); %现在的E是n*n的矩阵clear baseValue i j%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%nNum=n;while nNum>=2 %结点数目，逐一递减        A=zeros(nNum,1);    for i=1:nNum        A(i,1)=sum(E(i,:));    end    t=1; %tempQ的下标值。    tempQ=zeros(1,(nNum.*(nNum-1))./2);        %Qmax的第一个值是有34个社团时候的模块度    for i=1:nNum        Qmax(1,n-nNum+1)=Qmax(1,n-nNum+1)+(E(i,i)-A(i,1).*A(i,1));    end        %将第i个社团和第j个社团合并并计算此时的Q值，必须满足：[i<j]才不会重复计算。    for i=1:(nNum-1)        for j=(i+1):nNum                        %%%%%%%%%%%%计算合并后的社团结构的△Q值%%%%%%%%%%%%            %计算Q值，所有的临时Q值都要存储于tempQ数组里。            tempQ(1,t)=2*(E(i,j)-A(i,1).*A(j,1));                        if(tempQ(1,t)>deltaQ(1,n-nNum+1))                deltaQ(1,n-nNum+1)=tempQ(1,t);                maxI=i;                maxJ=j;%记录下最佳的方式，合并当前编号体制下的第i和第j个结点能够使得ΔQ最大。            end                        t=t+1;                    end    end            %%%%%%%%%%%%记录有用的数据%%%%%%%%%%%%    order(n-nNum+1,1)=maxI;    order(n-nNum+1,2)=maxJ;            %%%%%%%%%%%%为下一次循环做预处理工作%%%%%%%%%%%%    temp=E(maxI,:)+E(maxJ,:);    temp(:,maxI)=[];    temp(:,maxJ-1)=[];        inE=E(maxI,maxJ);    inEii=E(maxI,maxI);    inEjj=E(maxJ,maxJ);        E(maxI,:)=[];    E(maxJ-1,:)=[];    E(:,maxI)=[];    E(:,maxJ-1)=[];        %temp是两个合并了的社团数据加和之后暂时形成的一行数据，用于加在testData上        %合并之后的社团数目，递减了1    nNum=nNum-1;        E=[E;temp];    temp(:,nNum)=2*inE+inEii+inEjj;%注意，这里应该是2倍的。    %%%%%%%%%%%社团内部的边所占比例    E=[E temp'];    end%%%%%%%%%%%%清除无用的变量，减少干扰%%%%%%%%%%%%clear i j k     m n      maxI maxJclear t E nNum temp Aclear inE inEii inEjj%在同一张图上做两条曲线，红色为Q的变化过程，蓝色为ΔQ的变化过程plot(Qmax,'red');hold onplot(Qmax,'*');plot(deltaQ,'blue');plot(deltaQ,'*');hold offgrid on