Updated custom codes

custom/matching_dmperm.m

@@ -0,0 +1,45 @@
+function [p,q,r,s]=matching_dmperm(A)
+% MATCHING_DMPERM Implement an equivalent dmperm function with matching
+% 
+% The Matlab DMPERM function (as of R2007a) uses a recursive implementation
+% which can cause stack overflows for large matrices.  This function
+% implements the DMPERM function in terms of the non-recursive MATCHING
+% function in MatlabBGL.
+%
+% [p,q,r,s]=matching_dmperm(A) computes equivalent output to Matlab's
+% DMPERM function.
+%
+% See also DMPERM, SPRANK
+
+[nr nc]=size(A);
+GA=spaugment(A,0);
+rv=1:nr; cv=nr+1:nr+nc;
+m=matching(GA);
+cbar=find(m(cv)==0);
+rbar=find(m(rv)==0);
+d_cbar=setbfs(GA,cbar+nr);
+r1=find(mod(d_cbar(rv),2)==1 & d_cbar(rv)>0);
+c1=m(r1)-nr;
+d_rbar=setbfs(GA,rbar);
+c3=find(mod(d_rbar(cv),2)==1 & d_rbar(cv)>0);
+r3=m(c3+nr);
+rused=false(nr,1); rused(rbar)=1; rused(r1)=1; rused(r3)=1; 
+r2=find(~rused);
+cused=false(nc,1); cused(cbar)=1; cused(c1)=1; cused(c3)=1; 
+c2=find(~cused);
+B=A(r2,c2);
+ccB=components(B);
+[ignore ccp]=sort(ccB);
+p=[r1; r2(ccp); r3; rbar];
+q=[cbar; c1; c2(ccp); c3];
+end
+
+function d=setbfs(A,S)
+% SETBFS Compute BFS from a set of vertices
+n=size(A,1);
+A(n+1,n+1)=0; % extend A
+A(n+1,S)=1;
+d=bfs(A,n+1);
+d(end)=[];
+d(d>0)=d(d>0)-1;
+end

custom/path_histogram.m

@@ -1,4 +1,4 @@
-function [l c] = path_histogram(G,varargin)
+function [l c d deff] = path_histogram(G,varargin)
 % PATH_HISTOGRAM Compute a histogram of all shortest paths in graph G
 %
 % [l c] = path_histogram(G) computes all shortest paths in A one at a time
@@ -7,13 +7,18 @@
 % [l c] = path_histogram(G,struct('sample',N)) uses N random samples instead of an
 % exact solution.  (N=0 uses exact solution)
 %
+% [l,c,d,deff] = path_histogram(G,...) also returns the longest geodesics and
+%   longest effective geodesics ( for all nodes 
+%
 % Options:
 %   options.weighted: use Dijkstra algorithm for weighted paths [{0} | 1]
 %   options.intcount: use integer histogram bins [0 | {1}]
 %   options.nbins: number of bins when intcount=0 [{50} | positive integer]
 %   options.sample: use N iid samples instead of all sources
 %     [{0} | positive integer]
 %   options.verbose: output status for long runs [{0} | 1]
+%   options.effective_threshold: the threshold for effective geodesics 
+%     [{0.9} | 0 < float < 1]
 %
 % Example:
 %  load('graphs/cs-stanford.mat');
@@ -22,14 +27,16 @@
 
 %
 % David Gleich
-% Copyright, Stanford University, 2008
+% Copyright, Stanford University, 2008-2010
 %
 
 % History
 % 2008-03-14: Initial coding by David
+% :2010-05-28: Switched to including self in path count.
+%    Added option for geodesics and longest effective geodesics counts
 
 options=struct('weighted',0,'intcount',1,'nbins',50','sample',0,'istrans',0,...
-    'verbose',0);
+    'verbose',0,'effective_threshold',0.9);
 if ~isempty(varargin), options = merge_structs(varargin{1}, options); end
 
 if ~options.istrans, G = G'; end
@@ -42,6 +49,11 @@
 else vs=ceil(n*rand(options.sample,1)); 
 end
 
+if nargout>2
+  d=-1*ones(n,1);
+  deff=-1*ones(n,1);
+end
+
 % sample from all vertices
 H=sparse(n,1); t0=clock; p=1; Np=min(N,100);
 for vi=1:N
@@ -51,11 +63,21 @@
             100*(p-1)/Np, vi, N, N*dt/vi-dt); 
     end
     di=bfs(G,vi,bfsoptions);
-    di=di(di>0);             % only use reachable points
-    H=H+accumarray(di,1,[n 1],@sum,0,true);
-    bfsoptions.nocheck=1;    % don't repeat options check
+    di=di(di>=0);             % only use reachable points
+    H=H+accumarray(di+1,1,[n 1],@sum,0,true);
+    bfsoptions.nocheck=1;     % don't repeat options check
+    if nargout>2
+        d(vi) = max(di);
+        % compute the effective longest geodesic, which occurs
+        % when we can reach 90% of the vertices we can ever reach
+        reachable = length(di); 
+        reachable_at_d = cumsum(accumarray(di+1,1,[n,1],@sum,0,false));
+        deff(vi) = find(...
+            reachable_at_d>=options.effective_threshold*reachable,...
+            1,'first')-1;
+    end
 end
-l = find(H); c = nonzeros(H);
+l = find(H)-1; c = nonzeros(H);
 
 end