复杂网络聚类系数和平均路径长度计算的MATLAB源代码

申明：文章来自百度用户carrot_hy复杂网络的代码总共是三个m文件，复制如下:第一个文件，function [Cp_Global, Cp_Nodal] = CCM_ClusteringCoef(gMatrix, Types)% CCM_ClusteringCoef calculates clustering coefficients.% Input:% gMatrix adjacency matrix% Types type of graph:'binary','weighted','directed','all'(default).% Usage:% [Cp_Global, Cp_Nodal] = CCM_ClusteringCoef(gMatrix, Types) returns% clustering coefficients for all nodes "Cp_Nodal" and average clustering % coefficient of network "Cp_Global".% Example:% G = CCM_TestGraph1('nograph');% [Cp_Global, Cp_Nodal] = CCM_ClusteringCoef(G);% Note:% 1) one node have vaule 0, while which only has a neighbour or none.% 2) The dircted network termed triplets that fulfill the follow condition % as non-vacuous: j->i->k and k->i-j,if don't satisfy with that as% vacuous, just like: j->i,k->i and i->j,i->k. and the closed triplets% only j->i->k == j->k and k->i->j == k->j.% 3) 'ALL' type network code from Mika Rubinov's BCT toolkit.% Refer:% [1] Barrat et al. (2004) The architecture of the complex weighted networks. % [2] Wasserman,S.,Faust,K.(1994) Social Network Analysis: Methods and %Applications.% [3] Tore Opsahl and Pietro Panzarasa (2009). "Clustering in Weighted% Networks". Social Networks31(2).% See also CCM_Transitivity% Written by Yong Liu, Oct,2007% Center for Computational Medicine (CCM),% National Laboratory of Pattern Recognition (NLPR),% Institute of Automation,Chinese Academy of Sciences (IACAS), China.% Revise by Hu Yong, Nov, 2010% E-mail:% based on Matlab 2006a% $Revision: , Copywrite (c) 2007error(nargchk(1,2,nargin,'struct'));if(nargin < 2), Types = 'all'; endN = length(gMatrix);gMatrix(1:(N+1):end) = 0;%Clear self-edgesCp_Nodal = zeros(N,1); %Preallocateswitch(upper(Types))case 'BINARY'%Binary networkgMatrix = double(gMatrix > 0);%Ensure binary networkfor i = 1:Nneighbor = (gMatrix(i,:) > 0);Num = sum(neighbor);%number of neighbor nodestemp = gMatrix(neighbor, neighbor);if(Num > 1), Cp_Nodal(i) = sum(temp(:))/Num/(Num-1);endcase 'WEIGHTED'% Weighted network -- arithmetic mean for i = 1:Nneighbor = (gMatrix(i,:) > 0); n_weight =gMatrix(i,neighbor); Si = sum(n_weight);Num = sum(neighbor); if(Num > 1),n_weight = ones(Num,1)*n_weight;n_weight = n_weight + n_weight';n_weight = n_weight.*(gMatrix(neighbor,Cp_Nodal(i) = sum(n_weight(:))/(2*Si*(Num-1)); end end %case 'WEIGHTED'% Weighted network -- geometric mean% A = (gMatrix~= 0);% G3 = diag((gMatrix.A(1 ⑶)人3);)% A(A == 0) = inf; %close-triplet no exist,let CpNode=0 (A=inf)% CpNode = G3./(A.*(A-1));case 'DIRECTED', % Directed networkfor i = 1:Ninset = (gMatrix(:,i) > 0);outset = (gMatrix(i,:) > 0)'; %out-nodes setif(any(inset & outset)) allset = and(inset, outset);% Ensure aji*aik > 0,j belongs to inset,and kbelongs tooutsetend neighbor) > 0);%in-nodes settotal = sum(inset)*sum(outset) - sum(allset);tri = sum(sum(gMatrix(inset, outset))); Cp_Nodal(i)= tri./total;endend%Note: Directed & weighted network (from Mika Rubinov) case 'ALL',%All typeA = (gMatrix~= 0);G = gMatrix.A(1/3) + (gMatrix.')4(1/3); D = sum(A + A.',2);g3 = diag(G A 3)/2; D(g3 == 0) = inf; Cp=0 c3 = D.*(D-1) - 2*diag(AA2);Cp_Nodal = g3./c3;otherwise,%Eorr Msgerror('Type only four: "Binary","Weighted","Directed",and "All"');endCp_Global = sum(Cp_Nodal)/N;%%第二个文件：function [D_Global, D_Nodal] = CCM_AvgShortestPath(gMatrix, s, t)% CCM_AvgShortestPath generates the shortest distance matrix of source nodes% indice s to the target nodes indice t.% Input:% gMatrix symmetry binary connect matrix or weighted connectmatrix source nodes, default is 1:N% t target nodes, default is 1:N% Usage: % G = gMatrix + gMatrix'; %symmetrized% D = sum(G,2); %total degree % g3 = diag(GA3)/2; %number of triplet % D(g3 == 0) = inf; %3-cycles no exist,let % c3 = D.*(D-1) - 2*diag(gMatrixA2); %number of all possible 3-cycles% Cp_Nodal = g3./c3;%case 'DIRECTED', % Directed network -- clarity format (from Mika Rubinov, UNSW) Cp=0%adjacency matrix%total degree %number of triplet%3-cycles no% [D_Global, D_Nodal] = CCM_AvgShortestPath(gMatrix) returns the mean% shortest-path length of whole network D_Global,and the mean shortest-path % length of each node in the network% Example:% G = CCM_TestGraph1('nograph');% [D_Global, D_Nodal] = CCM_AvgShortestPath(G);% See also dijk, MEAN, SUM% Written by Yong Liu, Oct,2007% Modified by Hu Yong, Nov 2010% Center for Computational Medicine (CCM),% Based on Matlab 2008a% $Revision: , Copywrite (c) 2007% ###### Input check ######### error(nargchk(1,3,nargin,'struct'));N = length(gMatrix); if(nargin < 2 | isempty(s)), s = (1:N)';else s = s(:); endif(nargin < 3 | isempty(t)), t = (1:N)';else t = t(:); end % Calculate the shortest-path from s to all nodeD = dijk(gMatrix,s);%D(isinf(D)) = 0;D = D(:,t); %To target nodesD_Nodal = (sum(D,2)./sum(D>0,2)); % D_Nodal(isnan(D_Nodal)) = []; D_Global = mean(D_Nodal);第三个文件：function D = dijk(A,s,t)%DIJK Shortest paths from nodes 's' to nodes 't' using Dijkstra algorithm.% D = dijk(A,s,t)% A = n x n node-node weighted adjacency matrix of arc lengths% (Note: A(i,j) = 0 => Arc (i,j) does not exist;A(i,j) = NaN => Arc (i,j)exists with 0 weight)% s = FROM node indices% = [] (default), paths from all nodes% t = TO node indices% = [] (default), paths to all nodes% D = |s| x |t| matrix of shortest path distances from 's' to 't'% = [D(i,j)], where D(i,j) = distance from node 'i' to node 'j'% % (If A is a triangular matrix, then computationally intensive node % selection step not needed since graph is acyclic (triangularity is a% sufficient, but not a necessary, condition for a graph to be acyclic)% and A can have non-negative elements)% % (If |s| >> |t|, then DIJK is faster if DIJK(A',t,s) used, where D is now% transposed and P now represents successor indices)% % (Based on Fig. 4.6 in Ahuja, Magnanti, and Orlin, NetworkFlows, % Prentice-Hall, 1993, p. 109.) % Copyright (c) 1998-2000 by Michael G. Kay% Matlog Version 29-Aug-2000%% Modified by JBT, Dec 2000, to delete paths% Input Error Checking ******************************************************error(nargchk(1,3,nargin,'struct'));[n,cA] = size(A);if nargin < 2 | isempty(s), s = (1:n)'; else s = s(:); end if nargin < 3 | isempty(t), t = (1:n)'; else t = t(:); endif ~any(any(tril(A) ~= 0)) % A is upper triangular isAcyclic = 1;elseif ~any(any(triu(A) ~= 0)) % A is lower triangular isAcyclic = 2;else % Graph may not be acyclic isAcyclic = 0;end if n ~= cAerror('A must be a square matrix'); elseif ~isAcyclic &any(any(A < 0))error('A must be non-negative'); elseif any(s < 1 | s > n)error(['''s'' must be an integer between 1 and ',num2str(n)]); elseif any(t < 1 | t > n)error(['''t'' must be an integer between 1 and ',num2str(n)]);end % End (Input Error Checking)************************************************A = A'; % Use transpose to speed-up FIND for sparse AD = zeros(length(s),length(t));P = zeros(length(s),n);for i = 1:length(s) j = s(i);Di = Inf*ones(n,1); Di(j) = 0;isLab = logical(zeros(length(t),1)); if isAcyclic == 1nLab = j - 1;elseif isAcyclic == 2nLab = n - j;elsenLab = 0; UnLab = 1:n;isUnLab = logical(ones(n,1)); endwhile nLab < n & ~all(isLab)if isAcyclicDj = Di(j);else % Node selection[Dj,jj] = min(Di(isUnLab)); j = UnLab(jj);UnLab(jj) = []; isUnLab(j) = 0; end nLab = nLab + 1;if length(t) < n, isLab = isLab | (j == t); end[jA,kA,Aj] = find(A(:,j));Aj(isnan(Aj)) = 0;if isempty(Aj), Dk = Inf; else Dk = Dj + Aj; endP(i,jA(Dk < Di(jA))) = j; Di(jA) = min(Di(jA),Dk); if isAcyclic == 1 j = j +1;elseif isAcyclic == 2triangular A j = j - 1;end %disp( num2str( nLab )); endD(i,:) = Di(t)'; end % Increment node index for upper triangular A% Decrement node index for lower。