12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364 |
- function [CIJtree,CIJclus] = backbone_wu(CIJ,avgdeg)
- %BACKBONE_WU Backbone
- %
- % [CIJtree,CIJclus] = backbone_wu(CIJ,avgdeg)
- %
- % The network backbone contains the dominant connections in the network
- % and may be used to aid network visualization. This function computes
- % the backbone of a given weighted and undirected connection matrix CIJ,
- % using a minimum-spanning-tree based algorithm.
- %
- % input: CIJ, connection/adjacency matrix (weighted, undirected)
- % avgdeg, desired average degree of backbone
- % output:
- % CIJtree, connection matrix of the minimum spanning tree of CIJ
- % CIJclus, connection matrix of the minimum spanning tree plus
- % strongest connections up to an average degree 'avgdeg'
- %
- % NOTE: nodes with zero strength are discarded.
- % NOTE: CIJclus will have a total average degree exactly equal to
- % (or very close to) 'avgdeg'.
- % NOTE: 'avgdeg' backfill is handled slightly differently than in Hagmann
- % et al 2008.
- %
- % Reference: Hidalgo et al. (2007) Science 317, 482.
- % Hagmann et al. (2008) PLoS Biol
- %
- % Olaf Sporns, Indiana University, 2007/2008/2010/2012
- N = size(CIJ,1);
- CIJtree = zeros(N);
- % find strongest edge (note if multiple edges are tied, only use first one)
- [i,j,s] = find(max(max(CIJ))==CIJ); %#ok<*ASGLU>
- im = [i(1) i(2)];
- jm = [j(1) j(2)];
- % copy into tree graph
- CIJtree(im,jm) = CIJ(im,jm);
- in = im;
- out = setdiff(1:N,in);
- % repeat N-2 times
- for n=1:N-2
-
- % find strongest link between 'in' and 'out',ignore tied ranks
- [i,j,s] = find(max(max(CIJ(in,out)))==CIJ(in,out));
- im = in(i(1));
- jm = out(j(1));
-
- % copy into tree graph
- CIJtree(im,jm) = CIJ(im,jm); CIJtree(jm,im) = CIJ(jm,im);
- in = [in jm]; %#ok<AGROW>
- out = setdiff(1:N,in);
- end;
- % now add connections back, with the total number of added connections
- % determined by the desired 'avgdeg'
- CIJnotintree = CIJ.*~CIJtree;
- [a,b] = sort(nonzeros(CIJnotintree),'descend');
- cutoff = avgdeg*N - 2*(N-1);
- thr = a(cutoff);
- CIJclus = CIJtree + CIJnotintree.*(CIJnotintree>=thr);
|