R代码：随机删除点线并计算网络指标

由于总有朋友来问如何做此类模拟，现奉上本人写的代码，如有疑问，欢迎探讨。不过由于本人事务繁忙，不一定保证能立即回答。 另外，此代码只是本人写论文的时候用的一部分内容，虽然各位可以随意使用，但不保证一定正确。请各位朋友自己检查。 **如果网络很大，这段代码的效率是不高的，算个几天也有可能，请各位注意啦。 ***这个代码写于多年以前，随着R软件包的改进，有些计算步骤只需要用一个函数就能够解决了，所以请各位自己查找R的教程。 ########################################################## #对输入的网络进行随机顶点删除或随机连线删除模拟； #删除后计算组元数量、主组元规模和割点数，以便进行曲线拟合； #删除的条件是： #~~~~~~~~~~~按照指定比例删除点或线，并执行1000次 #~~~~~~~~~~~对每次结果计算: #~~~~~~~~~~~~~~~~~~~~~~~~~average shortest path length,transitivity, #~~~~~~~~~~~~~~~~~~~~~~~~~~degree distribution,hierarchical clustering #~~~~~~~~~~~~~~~~~~~~~~~~~~assortativeness #***注意：此处建立了全局效率的函数global.efficiency(G,N); ########################################################### ################################################################### #以下建立函数，计算全局效率global efficiency #全局效率=1/N(N-1)乘以网络中任意两点的最短路径的倒数之和。 #注意，一定要调用library(igraph) global.efficiency<-function(G){ N<-vcount(G); geo.recip<-1/shortest.paths(G); #计算网络中两两顶点之间的最短距离，并且取倒数，形成一个矩阵。 geo.recipNA<-ifelse(geo.recip==Inf,NA,geo.recip); #由于矩阵对角线上为1/0，显示为Inf，将其指定为缺失值NA #当然原来距离为0的顶点对在矩阵中也显示为NA，不影响计算结果 globalEff<-sum(geo.recipNA,na.rm=TRUE)/(N*(N-1)); #计算非缺失值的所有数据之和，并除以顶点数。 return(globalEff) #返回全局效率。 } ##################################################################### #==============设置初始条件================================ library(igraph); library(intergraph); workPath<-" ";#工作路径，自己填啦。 originGraph<-" ";#原始Pajek格式的.net网络文件的名字，注意带.net后缀啦。 delPercent<-0.2; #要删除的百分比，删点和删线分别都用相同的百分比 forMax.V<-1000;#删点模拟的次数 forMax.E<-1000;#删线模拟的次数，通常与删点模拟的次数相同 #========================================================== #================导入网络，边权值忽略====================== originIgraph<-read.graph(paste(workPath,originGraph,sep=""),"pajek"); origin.unweighted<-remove.edge.attribute(originIgraph,"weight"); #========================================================== #=====================删点模拟============================= number.vertices<-vcount(origin.unweighted);#顶点总数 number.del.V<-round(number.vertices*delPercent);#将要删除的顶点数 cat("Total vertices:",number.vertices,".",spe="","\n", number.del.V,"vertices will be attacked.","\n"); #--------------准备数组------------------------------------- #compV.number<-numeric(); #删点后的组元数量 #compV.maxSize<-numeric(); #删点后的最大组元规模 #articuV.number<-numeric();#删点后割点数； a.V<-numeric();#删点后协调相关系数 transi.V<-numeric();#删点后的网络整体云集系数 avePathL.V<-numeric();#删点后平均最短路径 #gEf.V<-numeric();#删点后全局效率 #triad.V<-list(); #删点后三方谱 #centralization.degree.V<-list(); #centralization.betweenness.V<-list(); #centralization.closeness.V<-list(); #ern.compV.number<-numeric(); #ern.compV.maxSize<-numeric();#最大组元的规模 #ern.articuV.number<-numeric();#割点数 ern.a.V<-numeric();#协调相关系数 ern.transi.V<-numeric();#整体云集系数 ern.avePathL.V<-numeric(); #ern.gEf.V<-numeric(); #ern.triad.V<-list(); #ern.centralization.degree.V<-list(); #.centralization.betweenness.V<-list(); #ern.centralization.closeness.V<-list(); #----------------------------------------------------------- #-------------------实施删点计算-------------------------------- for(i in 1: forMax.V){ startTime<-Sys.time(); #~~~~~~~~~~~~~~~~~~~~~~~制备删点网络~~~~~~~~~~~~~~~~~~~~~~~~~~~ set.seed(i); verDelNets<-delete.vertices( origin.unweighted, igraph.sample(1,number.vertices,number.del.V)); #对网络进行删点,igraph.sample()是igraph包内的随机数发生器 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #~~~~~~~~~~~~~~~~~~~~计算实际网络删点的参数~~~~~~~~~~~~~~~~~~~~ comp<-clusters(verDelNets,mode="strong"); #计算删点后的网络的组元 #组元数量原为comp$no，因不考虑仅有一点的组元，故用以下代码 #如希望最小组元为3顶点，则comp$csize>2 #compV.number[i]<-sum(ifelse(comp$csize>1,1,0)); #compV.maxSize[i]<-max(comp$csize);#最大组元的规模 #if(compV.maxSize[i]=="-Inf") compV.maxSize[i]<-0; #articuV.number[i]<-length(articulation.points(verDelNets));#割点数 a.V[i]<-assortativity.degree(verDelNets,directed=FALSE);#协调相关系数 transi.V[i]<-transitivity(verDelNets,type="global");#整体云集系数 avePathL.V[i]<-average.path.length( verDelNets,unconnected=TRUE);#unconnected=TRUE只考虑现有路径 #gEf.V[i]<-global.efficiency(verDelNets); #triad.V[[i]]<-triad.census(verDelNets); #centralization.degree.V[[i]]<-centralization.degree(verDelNets,normalized = TRUE); #centralization.betweenness.V[[i]]<-centralization.betweenness(verDelNets,normalized = TRUE); #centralization.closeness.V[[i]]<-centralization.closeness(verDelNets,normalized = TRUE); #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #~~~~~~~~~~~~~~~~~~~~与删点网络点线相同的ER随机网络~~~~~~~~~~~~~~ ern.V<-erdos.renyi.game(vcount(verDelNets),ecount(verDelNets),type="gnm"); #ern.comp<-clusters(ern.V,mode="strong"); #计算删点后的网络的组元 #组元数量原为ern.comp$no，因不考虑仅有一点的组元，故用以下代码 #如希望最小组元为3顶点，则ern.comp$csize>2 #ern.compV.number[i]<-sum(ifelse(ern.comp$csize>1,1,0)); #ern.compV.maxSize[i]<-max(ern.comp$csize);#最大组元的规模 #if(ern.compV.maxSize[i]=="-Inf") ern.compV.maxSize[i]<-0; #ern.articuV.number[i]<-length(articulation.points(ern.V));#割点数 ern.a.V[i]<-assortativity.degree(ern.V,directed=FALSE);#协调相关系数 ern.transi.V[i]<-transitivity(ern.V,type="global");#整体云集系数 ern.avePathL.V[i]<-average.path.length( ern.V,unconnected=TRUE);#unconnected=TRUE只考虑现有路径 #ern.gEf.V[i]<-global.efficiency(ern.V); #ern.triad.V[[i]]<-triad.census(ern.V); #ern.centralization.degree.V[[i]]<-centralization.degree(ern.V,normalized = TRUE); #ern.centralization.betweenness.V[[i]]<-centralization.betweenness(ern.V,normalized = TRUE); #ern.centralization.closeness.V[[i]]<-centralization.closeness(ern.V,normalized = TRUE); #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ endTime<-Sys.time(); cat("The",i,"of",forMax.V,"total iteration need",endTime-startTime,"seconds","\n"); } #============================================================ #==================删线模拟========================================== number.edges<-ecount(origin.unweighted);#连线总数 number.del.E<-round(number.edges*delPercent);#将要删除的连线数 cat("Total edges:",number.edges,spe="",".\n", number.del.E,"edges will be attacked.","\n"); #--------------准备数组---------------------------------------- #compE.number<-numeric(); #删线后的组元数量 #compE.maxSize<-numeric(); #删线后的最大组元规模 #articuE.number<-numeric();#删线后割点数； a.E<-numeric();#删线后协调相关系数 transi.E<-numeric();#删线后的网络云集系数 avePathL.E<-numeric();#删线后平均最短路径 #gEf.E<-numeric(); #删线后全局效率 #triad.E<-list(); #centralization.degree.E<-list(); #centralization.betweenness.E<-list(); #centralization.closeness.E<-list(); #ern.compE.number<-numeric(); #ern.compE.maxSize<-numeric();#最大组元的规模 #ern.articuE.number<-numeric();#割点数 ern.a.E<-numeric();#协调相关系数 ern.transi.E<-numeric();#整体云集系数 ern.avePathL.E<-numeric();#unconnected=TRUE只考虑现有路径 #ern.gEf.E<-numeric(); #ern.triad.E<-list(); #ern.centralization.degree.E<-list(); #ern.centralization.betweenness.E<-list(); #ern.centralization.closeness.E<-list(); #--------------------------------------------------------------- #-------------------实施删线计算-------------------------------- for(j in 1: forMax.E){ startTime<-Sys.time(); #~~~~~~~~~~~~~~~~~~~~~制备删线网络~~~~~~~~~~~~~~~~~~~~~~~~~~~ set.seed(j); edgDelNets<-delete.edges(origin.unweighted, igraph.sample(1,number.edges,number.del.E)); #对网络进行删线,igraph.sample()是igraph包内的随机数发生器 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #~~~~~~~~~~~~~~~~~~~~计算实际网络删线的参数~~~~~~~~~~~~~~~~~~ #comp<-clusters(edgDelNets,mode="strong"); #计算删线后的网络的组元 #组元数量原为comp$no，因不考虑仅有一点的组元，故用以下代码 #如希望最小组元为3顶点，则comp$csize>2 #compE.number[j]<-sum(ifelse(comp$csize>1,1,0)); #compE.maxSize[j]<-max(comp$csize);#最大组元的规模 #articuE.number[j]<-length(articulation.points(edgDelNets));#割点数 a.E[j]<-assortativity.degree(edgDelNets,directed=FALSE); transi.E[j]<-transitivity(edgDelNets,type="global"); avePathL.E[j]<-average.path.length( edgDelNets,unconnected=TRUE);#unconnected=TRUE只考虑现有路径 #gEf.E[j]<-global.efficiency(edgDelNets); #triad.E[[j]]<-triad.census(edgDelNets); #centralization.degree.E[[j]]<-centralization.degree(edgDelNets,normalized = TRUE); #centralization.betweenness.E[[j]]<-centralization.betweenness(edgDelNets,normalized = TRUE); #centralization.closeness.E[[j]]<-centralization.closeness(edgDelNets,normalized = TRUE); #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #~~~~~~~~~~~~~~~~~与删线网络点线相同的ER随机网络~~~~~~~~~~~~~~ ern.E<-erdos.renyi.game(vcount(edgDelNets),ecount(edgDelNets),type="gnm"); #ern.comp<-clusters(ern.E,mode="strong"); #计算删点后的网络的组元 #组元数量原为ern.comp$no，因不考虑仅有一点的组元，故用以下代码 #如希望最小组元为3顶点，则ern.comp$csize>2 #ern.compE.number[j]<-sum(ifelse(ern.comp$csize>1,1,0)); #ern.compE.maxSize[j]<-max(ern.comp$csize);#最大组元的规模 # if(ern.compE.maxSize[j]=="-Inf") ern.compE.maxSize[j]<-0; #ern.articuE.number[j]<-length(articulation.points(ern.E));#割点数 ern.a.E[j]<-assortativity.degree(ern.E,directed=FALSE);#协调相关系数 ern.transi.E[j]<-transitivity(ern.E,type="global");#整体云集系数 ern.avePathL.E[j]<-average.path.length( ern.E,unconnected=TRUE);#unconnected=TRUE只考虑现有路径 #ern.gEf.E[j]<-global.efficiency(ern.E); #ern.triad.E[[j]]<-triad.census(ern.E); #ern.centralization.degree.E[[j]]<-centralization.degree(ern.E,normalized = TRUE); #ern.centralization.betweenness.E[[j]]<-centralization.betweenness(ern.E,normalized = TRUE); #ern.centralization.closeness.E[[j]]<-centralization.closeness(ern.E,normalized = TRUE); #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ endTime<-Sys.time(); cat("The",j,"of",forMax.E,"total iteration need",endTime-startTime,"seconds","\n"); } #============================================================