Merge pull request #2 from mtb14006/mtb14006-patch-1

Add files via upload

MATLAB_Functions/HuskyTrack/gaimc/Contents.m

@@ -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

MATLAB_Functions/HuskyTrack/gaimc/bfs.m

@@ -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

MATLAB_Functions/HuskyTrack/gaimc/bipartite_matching.m

@@ -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
+
+
+