Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Browse files
Browse the repository at this point in the history
Add files via upload
Commit of All MATLAB Files and Toolbox Part1
- Loading branch information
Showing
100 changed files
with
5,311 additions
and
0 deletions.
There are no files selected for viewing
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,52 @@ | ||
%========================================= | ||
% Graph Algorithms in Matlab Code (gaimc) | ||
% Written by David Gleich | ||
% Version 1.0 (beta) | ||
% 2008-2009 | ||
%========================================= | ||
% | ||
% Search algorithms | ||
% dfs - depth first search | ||
% bfs - breadth first search | ||
% | ||
% Shortest path algorithms | ||
% dijkstra - Dijkstra's shortest path algorithm | ||
% | ||
% Minimum spanning tree algorithms | ||
% mst_prim - Compute an MST using Prim's algorithm | ||
% | ||
% Matching | ||
% bipartite_matching - Compute a maximum weight bipartite matching | ||
% | ||
% Connected components | ||
% scomponents - Compute strongly connected components | ||
% largest_component - Selects only the largest component | ||
% | ||
% Statistics | ||
% clustercoeffs - Compute clustering coefficients | ||
% dirclustercoeffs - Compute directed clustering coefficients | ||
% corenums - Compute core numbers | ||
% | ||
% Drawing | ||
% graph_draw - Draw an adjacency matrix (from Leon Peshkin) | ||
% | ||
% Helper functions | ||
% sparse_to_csr - Compressed sparse row arrays from a matrix | ||
% csr_to_sparse - Convert back to Matlab sparse matrices | ||
% load_gaimc_graph - Loads a sample graph from the library | ||
|
||
% David F. Gleich | ||
% Copyright, Stanford University, 2008-2009 | ||
|
||
% History | ||
% 2008-04-10: Initial version | ||
|
||
|
||
% TODO for release | ||
% Fix mlintrpt errors | ||
|
||
|
||
% Future todos | ||
% Implement weighted core nums | ||
% More testing | ||
% Implement all pairs shortest paths with Floyd Warshall |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,49 @@ | ||
function [d dt pred] = bfs(A,u,target) | ||
% BFS Compute breadth first search distances, times, and tree for a graph | ||
% | ||
% [d dt pred] = bfs(A,u) returns the distance (d) and the discover time | ||
% (dt) for each vertex in the graph in a breadth first search | ||
% starting from vertex u. | ||
% d = dt(i) = -1 if vertex i is not reachable from u | ||
% pred is the predecessor array. pred(i) = 0 if vertex (i) | ||
% is in a component not reachable from u and i != u. | ||
% | ||
% [...] = bfs(A,u,v) stops the bfs when it hits the vertex v | ||
% | ||
% Example: | ||
% load_gaimc_graph('bfs_example.mat') % use the dfs example from Boost | ||
% d = bfs(A,1) | ||
% | ||
% See also DFS | ||
|
||
% David F. Gleich | ||
% Copyright, Stanford University, 2008-20098 | ||
|
||
% History | ||
% 2008-04-13: Initial coding | ||
|
||
if ~exist('target','var') || isempty(full), target=0; end | ||
|
||
if isstruct(A), rp=A.rp; ci=A.ci; | ||
else [rp ci]=sparse_to_csr(A); | ||
end | ||
|
||
n=length(rp)-1; | ||
d=-1*ones(n,1); dt=-1*ones(n,1); pred=zeros(1,n); | ||
sq=zeros(n,1); sqt=0; sqh=0; % search queue and search queue tail/head | ||
|
||
% start bfs at u | ||
sqt=sqt+1; sq(sqt)=u; | ||
t=0; | ||
d(u)=0; dt(u)=t; t=t+1; pred(u)=u; | ||
while sqt-sqh>0 | ||
sqh=sqh+1; v=sq(sqh); % pop v off the head of the queue | ||
for ri=rp(v):rp(v+1)-1 | ||
w=ci(ri); | ||
if d(w)<0 | ||
sqt=sqt+1; sq(sqt)=w; | ||
d(w)=d(v)+1; dt(w)=t; t=t+1; pred(w)=v; | ||
if w==target, return; end | ||
end | ||
end | ||
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,240 @@ | ||
function [val m1 m2 mi]=bipartite_matching(varargin) | ||
% BIPARTITE_MATCHING Solve a maximum weight bipartite matching problem | ||
% | ||
% [val m1 m2]=bipartite_matching(A) for a rectangular matrix A | ||
% [val m1 m2 mi]=bipartite_matching(x,ei,ej,n,m) for a matrix stored | ||
% in triplet format. This call also returns a matching indicator mi so | ||
% that val = x'*mi. | ||
% | ||
% The maximum weight bipartite matching problem tries to pick out elements | ||
% from A such that each row and column get only a single non-zero but the | ||
% sum of all the chosen elements is as large as possible. | ||
% | ||
% This function is slightly atypical for a graph library, because it will | ||
% be primarily used on rectangular inputs. However, these rectangular | ||
% inputs model bipartite graphs and we take advantage of that stucture in | ||
% this code. The underlying graph adjency matrix is | ||
% G = spaugment(A,0); | ||
% where A is the rectangular input to the bipartite_matching function. | ||
% | ||
% Matlab already has the dmperm function that computes a maximum | ||
% cardinality matching between the rows and the columns. This function | ||
% gives us the maximum weight matching instead. For unweighted graphs, the | ||
% two functions are equivalent. | ||
% | ||
% Note: If ei and ej contain duplicate edges, the results of this function | ||
% are incorrect. | ||
% | ||
% See also DMPERM | ||
% | ||
% Example: | ||
% A = rand(10,8); % bipartite matching between random data | ||
% [val mi mj] = bipartite_matching(A); | ||
% val | ||
|
||
% David F. Gleich and Ying Wang | ||
% Copyright, Stanford University, 2008-2009 | ||
% Computational Approaches to Digital Stewardship | ||
|
||
% 2008-04-24: Initial coding (copy from Ying Wang matching_sparse_mex.cpp) | ||
% 2008-11-15: Added triplet input/output | ||
% 2009-04-30: Modified for gaimc library | ||
% 2009-05-15: Fixed error with empty inputs and triple added example. | ||
|
||
[rp ci ai tripi n m] = bipartite_matching_setup(varargin{:}); | ||
|
||
if isempty(tripi) | ||
error(nargoutchk(0,3,nargout,'struct')); | ||
else | ||
error(nargoutchk(0,4,nargout,'struct')); | ||
end | ||
|
||
|
||
if ~isempty(tripi) && nargout>3 | ||
[val m1 m2 mi] = bipartite_matching_primal_dual(rp, ci, ai, tripi, n, m); | ||
else | ||
[val m1 m2] = bipartite_matching_primal_dual(rp, ci, ai, tripi, n, m); | ||
end | ||
|
||
function [rp ci ai tripi n m]= bipartite_matching_setup(A,ei,ej,n,m) | ||
% convert the input | ||
|
||
if nargin == 1 | ||
if isstruct(A) | ||
[nzi nzj nzv]=csr_to_sparse(A.rp,A.ci,A.ai); | ||
else | ||
[nzi nzj nzv]=find(A); | ||
end | ||
[n m]=size(A); | ||
triplet = 0; | ||
elseif nargin >= 3 && nargin <= 5 | ||
nzi = ei; | ||
nzj = ej; | ||
nzv = A; | ||
if ~exist('n','var') || isempty(n), n = max(nzi); end | ||
if ~exist('m','var') || isempty(m), m = max(nzj); end | ||
triplet = 1; | ||
else | ||
error(nargchk(3,5,nargin,'struct')); | ||
end | ||
nedges = length(nzi); | ||
|
||
rp = ones(n+1,1); % csr matrix with extra edges | ||
ci = zeros(nedges+n,1); | ||
ai = zeros(nedges+n,1); | ||
if triplet, tripi = zeros(nedges+n,1); % triplet index | ||
else tripi = []; | ||
end | ||
|
||
% | ||
% 1. build csr representation with a set of extra edges from vertex i to | ||
% vertex m+i | ||
% | ||
rp(1)=0; | ||
for i=1:nedges | ||
rp(nzi(i)+1)=rp(nzi(i)+1)+1; | ||
end | ||
rp=cumsum(rp); | ||
for i=1:nedges | ||
if triplet, tripi(rp(nzi(i))+1)=i; end % triplet index | ||
ai(rp(nzi(i))+1)=nzv(i); | ||
ci(rp(nzi(i))+1)=nzj(i); | ||
rp(nzi(i))=rp(nzi(i))+1; | ||
end | ||
for i=1:n % add the extra edges | ||
if triplet, tripi(rp(i)+1)=-1; end % triplet index | ||
ai(rp(i)+1)=0; | ||
ci(rp(i)+1)=m+i; | ||
rp(i)=rp(i)+1; | ||
end | ||
% restore the row pointer array | ||
for i=n:-1:1 | ||
rp(i+1)=rp(i); | ||
end | ||
rp(1)=0; | ||
rp=rp+1; | ||
|
||
% | ||
% 1a. check for duplicates in the data | ||
% | ||
colind = false(m+n,1); | ||
for i=1:n | ||
for rpi=rp(i):rp(i+1)-1 | ||
if colind(ci(rpi)), error('bipartite_matching:duplicateEdge',... | ||
'duplicate edge detected (%i,%i)',i,ci(rpi)); | ||
end | ||
colind(ci(rpi))=1; | ||
end | ||
for rpi=rp(i):rp(i+1)-1, colind(ci(rpi))=0; end % reset indicator | ||
end | ||
|
||
|
||
function [val m1 m2 mi]=bipartite_matching_primal_dual(... | ||
rp, ci, ai, tripi, n, m) | ||
% BIPARTITE_MATCHING_PRIMAL_DUAL | ||
|
||
alpha=zeros(n,1); % variables used for the primal-dual algorithm | ||
beta=zeros(n+m,1); | ||
queue=zeros(n,1); | ||
t=zeros(n+m,1); | ||
match1=zeros(n,1); | ||
match2=zeros(n+m,1); | ||
tmod = zeros(n+m,1); | ||
ntmod=0; | ||
|
||
|
||
% | ||
% initialize the primal and dual variables | ||
% | ||
for i=1:n | ||
for rpi=rp(i):rp(i+1)-1 | ||
if ai(rpi) > alpha(i), alpha(i)=ai(rpi); end | ||
end | ||
end | ||
% dual variables (beta) are initialized to 0 already | ||
% match1 and match2 are both 0, which indicates no matches | ||
i=1; | ||
while i<=n | ||
% repeat the problem for n stages | ||
|
||
% clear t(j) | ||
for j=1:ntmod, t(tmod(j))=0; end | ||
ntmod=0; | ||
|
||
|
||
% add i to the stack | ||
head=1; tail=1; | ||
queue(head)=i; % add i to the head of the queue | ||
while head <= tail && match1(i)==0 | ||
k=queue(head); | ||
for rpi=rp(k):rp(k+1)-1 | ||
j = ci(rpi); | ||
if ai(rpi) < alpha(k)+beta(j) - 1e-8, continue; end % skip if tight | ||
if t(j)==0, | ||
tail=tail+1; queue(tail)=match2(j); | ||
t(j)=k; | ||
ntmod=ntmod+1; tmod(ntmod)=j; | ||
if match2(j)<1, | ||
while j>0, | ||
match2(j)=t(j); | ||
k=t(j); | ||
temp=match1(k); | ||
match1(k)=j; | ||
j=temp; | ||
end | ||
break; % we found an alternating path | ||
end | ||
end | ||
end | ||
head=head+1; | ||
end | ||
|
||
if match1(i) < 1, % still not matched, so update primal, dual and repeat | ||
theta=inf; | ||
for j=1:head-1 | ||
t1=queue(j); | ||
for rpi=rp(t1):rp(t1+1)-1 | ||
t2=ci(rpi); | ||
if t(t2) == 0 && alpha(t1) + beta(t2) - ai(rpi) < theta, | ||
theta = alpha(t1) + beta(t2) - ai(rpi); | ||
end | ||
end | ||
end | ||
|
||
for j=1:head-1, alpha(queue(j)) = alpha(queue(j)) - theta; end | ||
|
||
for j=1:ntmod, beta(tmod(j)) = beta(tmod(j)) + theta; end | ||
|
||
continue; | ||
end | ||
|
||
i=i+1; % increment i | ||
end | ||
|
||
val=0; | ||
for i=1:n | ||
for rpi=rp(i):rp(i+1)-1 | ||
if ci(rpi)==match1(i), val=val+ai(rpi); end | ||
end | ||
end | ||
noute = 0; % count number of output edges | ||
for i=1:n | ||
if match1(i)<=m, noute=noute+1; end | ||
end | ||
m1=zeros(noute,1); m2=m1; % copy over the 0 array | ||
noute=1; | ||
for i=1:n | ||
if match1(i)<=m, m1(noute)=i; m2(noute)=match1(i);noute=noute+1; end | ||
end | ||
|
||
if nargout>3 | ||
mi= false(length(tripi)-n,1); | ||
for i=1:n | ||
for rpi=rp(i):rp(i+1)-1 | ||
if match1(i)<=m && ci(rpi)==match1(i), mi(tripi(rpi))=1; end | ||
end | ||
end | ||
end | ||
|
||
|
||
|
Oops, something went wrong.