1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768 |
- function [R,D] = reachdist(CIJ)
- %REACHDIST Reachability and distance matrices
- %
- % [R,D] = reachdist(CIJ);
- %
- % The binary reachability matrix describes reachability between all pairs
- % of nodes. An entry (u,v)=1 means that there exists a path from node u
- % to node v; alternatively (u,v)=0.
- %
- % The distance matrix contains lengths of shortest paths between all
- % pairs of nodes. An entry (u,v) represents the length of shortest path
- % from node u to node v. The average shortest path length is the
- % characteristic path length of the network.
- %
- % Input: CIJ, binary (directed/undirected) connection matrix
- %
- % Outputs: R, reachability matrix
- % D, distance matrix
- %
- % Note: faster but more memory intensive than "breadthdist.m".
- %
- % Algorithm: algebraic path count.
- %
- %
- % Olaf Sporns, Indiana University, 2002/2007/2008
- % initialize
- R = CIJ;
- D = CIJ;
- powr = 2;
- N = size(CIJ,1);
- CIJpwr = CIJ;
- % Check for vertices that have no incoming or outgoing connections.
- % These are "ignored" by 'reachdist'.
- id = sum(CIJ,1); % indegree = column sum of CIJ
- od = sum(CIJ,2)'; % outdegree = row sum of CIJ
- id_0 = find(id==0); % nothing goes in, so column(R) will be 0
- od_0 = find(od==0); % nothing comes out, so row(R) will be 0
- % Use these columns and rows to check for reachability:
- col = setxor(1:N,id_0);
- row = setxor(1:N,od_0);
- [R,D,powr] = reachdist2(CIJ,CIJpwr,R,D,N,powr,col,row);
- % "invert" CIJdist to get distances
- D = powr - D+1;
- % Put 'Inf' if no path found
- D(D==(N+2)) = Inf;
- D(:,id_0) = Inf;
- D(od_0,:) = Inf;
- %----------------------------------------------------------------------------
- function [R,D,powr] = reachdist2(CIJ,CIJpwr,R,D,N,powr,col,row)
- % Olaf Sporns, Indiana University, 2002/2008
- CIJpwr = CIJpwr*CIJ;
- R = double(R | ((CIJpwr)~=0));
- D = D+R;
- if ((powr<=N)&&(~isempty(nonzeros(R(row,col)==0))))
- powr = powr+1;
- [R,D,powr] = reachdist2(CIJ,CIJpwr,R,D,N,powr,col,row);
- end;
|