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

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

申明：文章来自百度用户carrot_hy

复杂网络的代码总共是三个m文件，复制如下：

第一个文件，CCM_ClusteringCoef.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: 1.0, Copywrite (c) 2007

error(nargchk(1,2,nargin,'struct'));

if(nargin < 2), Types = 'all'; end

N = length(gMatrix);

gMatrix(1:(N+1):end) = 0;%Clear self-edges

Cp_Nodal = zeros(N,1); %Preallocate

switch(upper(Types))

case 'BINARY'%Binary network

gMatrix = double(gMatrix > 0);%Ensure binary network

for i = 1:N

neighbor = (gMatrix(i,:) > 0);

Num = sum(neighbor);%number of neighbor nodes

temp = gMatrix(neighbor, neighbor);

if(Num > 1), Cp_Nodal(i) = sum(temp(:))/Num/(Num-1); end

end

case 'WEIGHTED'% Weighted network -- arithmetic mean

for i = 1:N

neighbor = (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, neighbor) > 0);

Cp_Nodal(i) = sum(n_weight(:))/(2*Si*(Num-1));

end

end

%case 'WEIGHTED'% Weighted network -- geometric mean

% A = (gMatrix~= 0);

% G3 = diag((gMatrix.^(1/3) )^3);)

% A(A == 0) = inf; %close-triplet no exist,let CpNode=0 (A=inf)

% CpNode = G3./(A.*(A-1));

case 'DIRECTED', % Directed network

for i = 1:N

inset = (gMatrix(:,i) > 0); %in-nodes set

outset = (gMatrix(i,:) > 0)'; %out-nodes set

if(any(inset & outset))

allset = and(inset, outset);

% Ensure aji*aik > 0,j belongs to inset,and k belongs to outset total = sum(inset)*sum(outset) - sum(allset);

tri = sum(sum(gMatrix(inset, outset)));

Cp_Nodal(i) = tri./total;

end

end

%case 'DIRECTED', % Directed network -- clarity format (from Mika Rubinov, UNSW) % G = gMatrix + gMatrix'; %symmetrized

% D = sum(G,2); %total degree

% g3 = diag(G^3)/2; %number of triplet

% D(g3 == 0) = inf; %3-cycles no exist,let Cp=0

% c3 = D.*(D-1) - 2*diag(gMatrix^2); %number of all possible 3-cycles

% Cp_Nodal = g3./c3;

%Note: Directed & weighted network (from Mika Rubinov)

case 'ALL',%All type

A = (gMatrix~= 0); %adjacency matrix

G = gMatrix.^(1/3) + (gMatrix.').^(1/3);

D = sum(A + A.',2); %total degree

g3 = diag(G^3)/2; %number of triplet

D(g3 == 0) = inf; %3-cycles no exist,let Cp=0

c3 = D.*(D-1) - 2*diag(A^2);

Cp_Nodal = g3./c3;

otherwise,%Eorr Msg

error('Type only four: "Binary","Weighted","Directed",and "All"');

end

Cp_Global =

sum(Cp_Nodal)/N;

%%

第二个文件：CCM_AvgShortestPath.m

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 connect matrix

% s source nodes, default is 1:N

% t target nodes, default is 1:N

% Usage:

% [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: 1.0, 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(:); end

if(nargin < 3 | isempty(t)), t = (1:N)';

else t = t(:); end

% Calculate the shortest-path from s to all node

D = dijk(gMatrix,s);%D(isinf(D)) = 0;

D = D(:,t); %To target nodes

D_Nodal = (sum(D,2)./sum(D>0,2));

% D_Nodal(isnan(D_Nodal)) = [];

D_Global = mean(D_Nodal);

第三个文件： dijk.m

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, Network Flows,

% Prentice-Hall, 1993, p. 109.)

% Copyright (c) 1998-2000 by Michael G. Kay

% Matlog Version 1.3 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(:); end

if ~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 ~= cA

error('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 A

D = 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 == 1

nLab = j - 1;

elseif isAcyclic == 2

nLab = n - j;

else

nLab = 0;

UnLab = 1:n;

isUnLab = logical(ones(n,1));

end

while nLab < n & ~all(isLab)

if isAcyclic

Dj = 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; end

P(i,jA(Dk < Di(jA))) = j;

Di(jA) = min(Di(jA),Dk);

if isAcyclic == 1 % Increment node index for upper triangular A

j = j + 1;

elseif isAcyclic == 2 % Decrement node index for lower triangular A j = j - 1;

end

%disp( num2str( nLab ));

end

D(i,:) = Di(t)';

end