diff --git a/.gitignore b/.gitignore deleted file mode 100644 index d0ab673..0000000 --- a/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -*~ -*.asv - diff --git a/Contents.m b/Contents.m deleted file mode 100644 index 0307020..0000000 --- a/Contents.m +++ /dev/null @@ -1,57 +0,0 @@ -%========================================= -% Graph Algorithms in Matlab Code (gaimc) -% Written by David Gleich -% Version 1.1 -% 2008-2009 -%========================================= -% -% Search algorithms -% dfs - depth first search -% bfs - breadth first search -% -% Shortest path algorithms -% dijkstra - Dijkstra's shortest path algorithm -% floydwarshall - Floyd-Warshall all 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 -% Update demo with floydwarshall -% Update airports with floydwarshall - - -% Future todos -% Implement geometric random graph generation -% Implement erdos reyni graph generatin -% Implement spectral graph partitioning -% Implement a max-flow solver -% Implement weighted core nums -% More testing diff --git a/bfs.m b/bfs.m deleted file mode 100644 index 93c7d59..0000000 --- a/bfs.m +++ /dev/null @@ -1,49 +0,0 @@ -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 diff --git a/bipartite_matching.m b/bipartite_matching.m deleted file mode 100644 index c35e5d7..0000000 --- a/bipartite_matching.m +++ /dev/null @@ -1,240 +0,0 @@ -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 - - - diff --git a/build/buildgaimc.m b/build/buildgaimc.m deleted file mode 100644 index e2819cb..0000000 --- a/build/buildgaimc.m +++ /dev/null @@ -1,87 +0,0 @@ -%% build the gaimc distribution -function buildgaimc(varargin) -if ~isempty(strmatch('all',varargin)) - buildlist = {'test','pages','zip'}; -else - buildlist = varargin; -end -tobuild = @(type,buidlist) ~isempty(strmatch(type,buildlist)); -gaimcversion = '1.0'; -mypath = fileparts(mfilename('fullpath')); -oldpath = pwd; -matlablpath=path; -try - if tobuild('test',buildlist) - fprintf('Running tests...\n\n'); - cd(mypath) % move to where we think we are - cd ../ % in main gaimc directory - addpath(pwd); % add it to the path - cd test % run tests - try - test_main - catch ME - fprintf('*************\n') - fprintf('Tests failed!\n'); - fprintf('*************\n') - fprintf('\n'); - fprintf('Build halted.\n'); - fprintf('\n'); - fprintf('Test error:\n'); - rethrow(ME); - end - end - - if tobuild('pages',buildlist) - fprintf('Building demos...\n\n'); - cd(mypath) % move to where we think we are - cd ../ % in main gaimc directory - addpath(pwd); % add it to the path - cd demo - try - publish('demo'); - publish('airports'); - publish('performance_comparison_simple'); - if tobuild('pages_fullperf') - publish('performance_comparison'); - end - catch ME - fprintf('******************\n') - fprintf('Publishing failed!\n'); - fprintf('******************\n') - fprintf('\n'); - fprintf('Build halted.\n'); - fprintf('\n'); - fprintf('Test error:\n'); - rethrow(ME); - end - end - - if tobuild('zip',buildlist) - cd(mypath) % move to where we think we are - try - cd ../../ - [status,result] = system(... - sprintf('zip -r gaimc/build/gaimc-%s.zip gaimc/*.m gaimc/demo gaimc/test gaimc/graphs',... - gaimcversion)); - result - catch ME - fprintf('******************\n') - fprintf('Zipping failed!\n'); - fprintf('******************\n') - fprintf('\n'); - fprintf('Build halted.\n'); - fprintf('\n'); - fprintf('Error:\n'); - rethrow(ME); - end - end - -catch ME - cd(oldpath); % move to where we think we are - path(matlabpath); - rethrow(ME) -end - -cd(oldpath); -path(matlabpath); - diff --git a/build/matlab_central_notes.txt b/build/matlab_central_notes.txt deleted file mode 100644 index bf4a4cb..0000000 --- a/build/matlab_central_notes.txt +++ /dev/null @@ -1,35 +0,0 @@ -Picture: minnesota graph (demo_08.png) - -Title : gaimc : Graph Algorithms In Matlab Code -Summary (100 chars) : -Efficient pure-Matlab implementations of graph algorithms to complement MatlabBGL's mex functions. -Description (5000 chars) : - -While MatlabBGL uses the Boost Graph Library for efficient graph routines, -gaimc implements everything in pure Matlab code. While the routines are -slower, they aren't as slow as I initially thought. Since people often -have problems getting MatlabBGL to compile on new versions of Matlab -or on new architectures, this library is then a complement to MatlabBGL. - -See the published M-files for a few examples of the capabilities. - -Functions - depth first search (dfs) - breadth first search (bfs) - connected components (scomponents) - maximum weight bipartite matching (bipartite_matching) - Dijkstra's shortest paths (dijkstra) - Prim's minimum spanning tree (mst_prim) - clustering coefficients (clustercoeffs) - directed clustering coefficients (dirclustercoeffs) - core numbers (corenums) - -The project is maintained at github : http://github.com/dgleich/gaimc/tree/master - -Tags: - -graph, network, dijkstra, core numbers, prim, dfs, bfs, connected components, -strongly connected components, clustering coefficients, directed clustering coefficients, -directed network, directed graph, depth first search, breadth first search, -shortest paths, minimum spanning tree, mst - diff --git a/clustercoeffs.m b/clustercoeffs.m deleted file mode 100644 index 19fc771..0000000 --- a/clustercoeffs.m +++ /dev/null @@ -1,74 +0,0 @@ -function cc=clustercoeffs(A,weighted,normalized) -% CLUSTERCOEFFS Compute undirected clustering coefficients for a graph -% -% ccfs=clustercoeffs(A) compute normalized, weighted clustering -% coefficients from a graph represented by a symmetric adjacency matrix A. -% -% ccfs=clustering(A,weighted,normalized) takes optional parameters to -% control normalization and weighted computation. If normalized=0 or -% false, then the computation is not normalized by d*(d-1) in the -% unweighted case. If weighted=0 or false, then the weights of the matrix -% A are ignored. Either parameter will assume it's default value if you -% specify an empty matrix. -% -% See also DIRCLUSTERCOEFFS -% -% Example: -% load_gaimc_graph('clique-10'); -% cc = clustercoeffs(A) % they are all equal! as we expect in a clique - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2009-05-15: First history comment, originally written in 2008 - -if ~exist('normalized','var') || isempty(normalized), normalized=true; end -if ~exist('weighted','var') || isempty(weighted), weighted=true; end -donorm=1; usew=1; -if ~normalized, donorm=0; end -if ~weighted, usew=0; end - -if isstruct(A) - rp=A.rp; ci=A.ci; %ofs=A.offset; - if usew, ai=A.ai; end -else - if ~isequal(A,A'), error('gaimc:clustercoeffs',... - 'only undirected (symmetric) inputs allowed: see dirclustercoeffs'); - end - if usew, [rp ci ai]=sparse_to_csr(A); - else [rp ci]=sparse_to_csr(A); - end - if any(ai)<0, error('gaimc:clustercoeffs',... - ['only positive edge weights allowed\n' ... - 'try clustercoeffs(A,0) for an unweighted comptuation']); - end -end -n=length(rp)-1; - -cc=zeros(n,1); ind=false(n,1); cache=zeros(n,1); -ew=1; ew2=1; -for v=1:n - for rpi=rp(v):rp(v+1)-1 - w=ci(rpi); if usew, ew=ai(rpi); end - if v~=w, ind(w)=1; cache(w)=ew^(1/3); end - end - curcc=0; d=rp(v+1)-rp(v); - % run two steps of bfs to try and find triangles. - for rpi=rp(v):rp(v+1)-1 - w=ci(rpi); if v==w, d=d-1; continue; end % discount self-loop - for rpi2=rp(w):rp(w+1)-1 - x=ci(rpi2); if x==w, continue; end - if ind(x) - if usew, ew=ai(rpi); ew2=ai(rpi2); end - curcc=curcc+ew^(1/3)*ew2^(1/3)*cache(x); - end - end - end - if donorm&&d>1, cc(v)=curcc/(d*(d-1)); - elseif d>1, cc(v)=curcc; - end - for rpi=rp(v):rp(v+1)-1, w=ci(rpi); ind(w)=0; end % reset indicator -end - - diff --git a/convert_sparse.m b/convert_sparse.m deleted file mode 100644 index f3fdcfb..0000000 --- a/convert_sparse.m +++ /dev/null @@ -1,22 +0,0 @@ -function As = convert_sparse(A) -% CONVERT_SPARSE Convert a sparse matrix to the native gaimc representation -% -% As = convert_sparse(A) returns a struct with the three arrays defining -% the compressed sparse row structure of A. -% -% Example: -% load('graphs/all_shortest_paths_example') -% As = convert_sparse(A) -% -% See also SPARSE_TO_CSR SPARSE - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2009-04-29: Initial coding - -[rp ci ai] = sparse_to_csr(A); -As.rp = rp; -As.ci = ci; -As.ai = ai; diff --git a/corenums.m b/corenums.m deleted file mode 100644 index df9ecba..0000000 --- a/corenums.m +++ /dev/null @@ -1,72 +0,0 @@ -function [d rt]=corenums(A) -% CORENUMS Compute the core number for each vertex in the graph. -% -% [cn rt]=corenums(A) returns the core numbers for each vertex of the graph -% A along with the removal order of the vertex. The core number is the -% largest integer c such that vertex v exists in a graph where all -% vertices have degree >= c. The vector rt returns the removal time -% for each vertex. That is, vertex vi was removed at step rt[vi]. -% -% This method works on directed graphs but gives the in-degree core number. -% To get the out-degree core numbers, call corenums(A'). -% -% The linear algorithm comes from: -% Vladimir Batagelj and Matjaz Zaversnik, "An O(m) Algorithm for Cores -% Decomposition of Networks." Sept. 1 2002. -% -% Example: -% load_gaimc_graph('cores_example'); % the graph A has three components -% corenums(A) -% - -% See also WCORENUMS - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2008-04-21: Initial Coding - -if isstruct(A), rp=A.rp; ci=A.ci; %ofs=A.offset; -else [rp ci]=sparse_to_csr(A); -end -n=length(rp)-1; -nz=length(ci); - -% the algorithm removes vertices by computing bucket sort on the degrees of -% all vertices and removing the smallest. - -% compute in-degrees and maximum indegree -d=zeros(n,1); maxd=0; rt=zeros(n,1); -for k=1:nz, newd=d(ci(k))+1; d(ci(k))=newd; if newd>maxd; maxd=newd; end, end - -% compute the bucket sort -dp=zeros(maxd+2,1); vs=zeros(n,1); vi=zeros(n,1); % degree position, vertices -for i=1:n, dp(d(i)+2)=dp(d(i)+2)+1; end % plus 2 because degrees start at 0 -dp=cumsum(dp); dp=dp+1; -for i=1:n, vs(dp(d(i)+1))=i; vi(i)=dp(d(i)+1); dp(d(i)+1)=dp(d(i)+1)+1; end -for i=maxd:-1:1, dp(i+1)=dp(i); end - -% start the algorithm -t=1; -for i=1:n - v = vs(i); dv = d(v); rt(v)=t; t=t+1; - for rpi=rp(v):rp(v+1)-1 - w=ci(rpi); dw=d(w); - if dw<=dv, % we already removed w - else % need to remove edge (v,w), which decreases d(w) - % swap w with the vertex at the head of its degree - pw=vi(w); % get the position of w - px=dp(dw+1); % get the pos of the vertex at the head of dw list - x=vs(px); - % swap w, x - vs(pw)=x; vs(px)=w; vi(w)=px; vi(x)=pw; - % decrement the degree of w and increment the start of dw - d(w)=dw-1; - dp(dw+1)=px+1; - end - end -end - - - diff --git a/cosineknn.m b/cosineknn.m deleted file mode 100644 index b0f7ed9..0000000 --- a/cosineknn.m +++ /dev/null @@ -1,94 +0,0 @@ -function S = cosineknn(A,K) -% COSINEKNN Compute k-nearest neighbors with a cosine distance metric -% -% S = cosineknn(A,k) computes a cosine similarity metric between the -% vertices of A or the upper half of a bipartite graph A -% - -% David F. Gleich -% Copyright, Stanford University, 2009 - -% History -% 2009-08-12: Initial version - -% TODO support struct input -% TODO support other similarity functions -% TODO allow option for diagonal similarity too -% TODO allow option for triplet output - - -if isstruct(A), - error('cosineknn:notSupported','struct input is not supported yet'); -else - [rp ci ai]=sparse_to_csr(A); check=1; - [rpt cit ait]=sparse_to_csr(A'); - [n,m] = size(A); -end - -% compute a row norm (which means we accumlate the column indices of At) -rn = accumarray(cit,ait.^2); -rn = sqrt(rn); -rn(rn>eps(1)) = 1./rn(rn>eps(1)); % do division once - -% the output is a n-by-n matrix, allocate storage for it -si = zeros(n*K,1); -sj = zeros(n*K,1); -sv = zeros(n*K,1); -nused = 0; - -dwork = zeros(n,1); -iwork = zeros(n,1); -iused = zeros(n,1); -for i=1:n - % for each row of A, compute an inner product against all rows of A. - % to do this efficiently, we know that the inner-product will - % only be non-zero if two rows of A overlap in a column. Thus, in - % row i of A, we look at all non-zero columns, and then look at all - % rows touched by those columns in the A' matrix. We compute - % the similarity metric, then sort and only store the top k. - - curused = 0; % track how many non-zeros we compute during this operation - - for ri=rp(i):rp(i+1)-1 - % find all columns j in row i - j = ci(ri); - aij = ai(ri)*rn(i); - for rit=rpt(j):rpt(j+1)-1 - % find all rows k in column j - k = cit(rit); - if k==i, continue; end % skip over the standard entries - akj = ait(rit)*rn(k); - if iwork(k)>0 - % we already have a non-zero between row i and row k - dwork(iwork(k)) = dwork(iwork(k)) + aij*akj; - else - % we need to add a non-zero betwen row i and row k - curused = curused + 1; - iwork(k) = curused; - dwork(curused) = aij*akj; - iused(curused) = k; - end - end - end - % don't sort if fewer than K elements - if curused < K, - sperm = 1:curused; - else - [sdwork sperm] = sort(dwork(1:curused),1,'descend'); - end - for k=1:min(K,length(sperm)) - nused = nused + 1; - si(nused) = i; - sj(nused) = iused(sperm(k)); - sv(nused) = dwork(sperm(k)); - end - % reset the iwork array, no need to reset dwork, as it's just - % overwritten - for k=1:curused - iwork(iused(k)) = 0; - end -end -si = si(1:nused); -sj = sj(1:nused); -sv = sv(1:nused); -S = sparse(si,sj,sv,n,n); \ No newline at end of file diff --git a/csr_to_sparse.m b/csr_to_sparse.m deleted file mode 100644 index 3c3e508..0000000 --- a/csr_to_sparse.m +++ /dev/null @@ -1,54 +0,0 @@ -function [nzi,nzj,nzv] = csr_to_sparse(rp,ci,ai,ncols) -% CSR_TO_SPARSE Convert from compressed row arrays to a sparse matrix -% -% A = csr_to_sparse(rp,ci,ai) returns the sparse matrix represented by the -% compressed sparse row representation rp, ci, and ai. The number of -% columns of the output sparse matrix is max(max(ci),nrows). See the call -% below. -% -% A = csr_to_sparse(rp,ci,ai,ncol) While we can infer the number of rows -% in the matrix from this expression, you may want a -% different number of -% -% [nzi,nzj,nzv] = csr_to_sparse(...) returns the arrays that feed the -% sparse call in matlab. You can use this to avoid the sparse call and -% customize the behavior. -% -% This command "inverts" the behavior of sparse_to_csr. -% Repeated entries in the matrix are summed, just like sparse does. -% -% See also SPARSE SPARSE_TO_CSR -% -% Example: -% A=sparse(6,6); A(1,1)=5; A(1,5)=2; A(2,3)=-1; A(4,1)=1; A(5,6)=1; -% [rp ci ai]=sparse_to_csr(A); -% A2 = csr_to_sparse(rp,ci,ai) - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2009-05-01: Initial version -% 2009-05-16: Documentation and example - -nrows = length(rp)-1; -nzi = zeros(length(ci),1); -for i=1:nrows - for j=rp(i):rp(i+1)-1 - nzi(j) = i; - end -end - -if nargout<2, - if nargin>3, - nzi = sparse(nzi,ci,ai,nrows,ncols); - else - % we make the matrix square unless there are more columns - ncols = max(max(ci),nrows); - if isempty(ncols), ncols=0; end - nzi = sparse(nzi,ci,ai,nrows,ncols); - end -else - nzj = ci; - nzv = ai; -end diff --git a/demo/airports.m b/demo/airports.m deleted file mode 100644 index af1970a..0000000 --- a/demo/airports.m +++ /dev/null @@ -1,119 +0,0 @@ -%% The US airport network -% THe North American airport network is an interesting graph to examine. -% The source for this data was a file on Brendan Frey's affinity -% propagation website. A(i,j) is the negative travel time between two -% airports. Although the data didn't include the airport locations, I used -% the Yahoo! Geocoding API to generate a latitude and longitude for each -% airport. - -%% The data -load_gaimc_graph('airports'); -%% -% Plot a histogram for all route time estimates -[si, ti, rt] = find(A); -hist(-rt,100); % times are stored as negative values - -%% -% Find the lengthiest route -[val,ind] = max(-rt) -{labels{si(ind)} labels{ti(ind)}} - -%% -% Some of the routes include stop overs, so it's probable that is what -% we find in this case. - -%% Graph analysis: connected? -% One of the first questions about any graph should be if it's connected or -% not. - -max(scomponents(A)) - -%% -% There is only one connected component, so the graph is connected. - -%% Distance instead of time -% Let's see how the edges correlate distance with estimated travel time. -[ai aj te] = find(A); -de = distance(xy(ai,:), xy(aj,:)); -plot(de,-te,'.'); -xlabel('distance (arclength)'); ylabel('time (?)'); - -%% -% Wow! It's all over the place, but there is a lower bound. Some of -% these routes can include stop - -%% Minimum spanning tree -% This section repeats and extends some analysis in the overall gaimc demo. -% First, let's recompute the minimum spanning tree based on travel time. -load_gaimc_graph('airports') -A = -A; % we store the negative travel time -A = max(A,A'); % travel time isn't symmetric -T = mst_prim(A); -clf; -gplot(T,xy); - - -% These next lines plot a map of the US with states colored. -ax = worldmap('USA'); -load coast -geoshow(ax, lat, long,... - 'DisplayType', 'polygon', 'FaceColor', [.45 .60 .30]) -states = shaperead('usastatelo', 'UseGeoCoords', true); -faceColors = makesymbolspec('Polygon',... - {'INDEX', [1 numel(states)], 'FaceColor', polcmap(numel(states))}); - geoshow(ax, states, 'DisplayType', 'polygon', 'SymbolSpec', faceColors) -set(gcf,'Position', [ 52 234 929 702]); -%% -% That's the US, now we need to plot our data on top of it. -[X,Y] = gplot(T,xy); % get the information to reproduce a gplot -plotm(Y,X,'k.-','LineWidth',1.5); % plot the lines on the map - -%% -% Let's just look at the continential US too. -clf; ax = usamap('conus'); -states = shaperead('usastatelo', 'UseGeoCoords', true, 'Selector',... - {@(name) ~any(strcmp(name,{'Alaska','Hawaii'})), 'Name'}); -faceColors = makesymbolspec('Polygon',... - {'INDEX', [1 numel(states)], 'FaceColor', polcmap(numel(states))}); -geoshow(ax, states, 'DisplayType', 'polygon', 'SymbolSpec', faceColors) -framem off; gridm off; mlabel off; plabel off -set(gcf,'Position', [ 52 234 929 702]); -plotm(Y,X,'k.-','LineWidth',1.5); % plot the lines on the map -%% -% One interesting aspect of this map is that major airline hubs (Chicago, -% New York, etc. are not well represented. One possible explaination is -% that they have larger delays than other regional airports. - -%% Honolulu to St. Johns? -% Before, we saw that the lengthiest route was between St. John's and -% Honolulu. But, does the network have a better path to follow between -% these cities? Let's check using Dijkstra's shortest path algorithm! - -% Reload the network to restore it. -load_gaimc_graph('airports') -A = -A; - -% find the longest route again -[si, ti, rt] = find(A); -[val,ind] = max(rt); % we've already negated above, so no need to redo it -start = si(ind); -dest = ti(ind); -[d pred] = dijkstra(A,start); % compute the distance to everywhere from St. Johns -d(dest) - -%% -% That value is considerably shorter than the direct time. How do we find -% this awesome route? -path =[]; u = dest; while (u ~= start) path=[u path]; u=pred(u); end -fprintf('%s',labels{start}); -for i=path; fprintf(' --> %s', labels{i}); end, fprintf('\n'); - -%% All pairs shortest paths -% At this point, the right thing to do would be to recompute each edge in -% the network using an all-pairs shortest path algorithm, and then look at -% how distance correlates with time in that network. However, I haven't -% had time to implement that algorithm yet. - -%% Conclusion -% Hopefully, you will agree that network algorithms are a powerful way to -% look at the relationships between airports! diff --git a/demo/demo.m b/demo/demo.m deleted file mode 100644 index 7778ae7..0000000 --- a/demo/demo.m +++ /dev/null @@ -1,310 +0,0 @@ -%% Demo of gaimc - 'Graph Algorithms In Matlab Code' -% Matlab includes great algorithms to work with sparse matrices but does -% provide a reasonable set of algorithms to work with sparse matrices as -% graph data structures. My other project -- MatlabBGL -- provides a -% high-performance solution to this problem by directly interfacing the -% Matlab sparse matrix data structure with the Boost Graph Library. -% That library, however, suffers from enormous complication because -% it must be compiled for each platform. The Boost Graph Library -% heavily uses advanced C++ features that impair easy portability -% between platforms. In contrast, the gaimc library is implemented in -% pure Matlab code, making it completely portable. -% -% The cost of the portability for this library is a 2-4x slowdown in -% the runtime of the algorithms as well as significantly fewer -% algorithms to choose from. - -%% Sparse matrices as graphs -% To store the connectivity structure of the graph, gaimc uses the -% adjacency matrix of a graph. -% -% A graph is represented by a set of vertices and a set of -% edges between the vertices. Often, we write $G = (V,E)$ to denote the -% graph, the set of vertices, and the set of edges, respectively. In -% gaimc, like in my other package MatlabBGL, we represent graphs with their -% adjacency matrices. This representation is handy in Matlab because -% Matlab is rather efficient at working with the large sparse matrices that -% typically arise as adjacency matrices. -% -% To convert from $G=(V,E)$ to an adjacency matrix, we identify each vertex -% with a row of the matrix via a bijective map. The adjacency matrix is -% then a $|V| \times |V|$ matrix called A. The entry A(i,j) = 1 for -% any edge between in $E$ and 0 otherwise. Let's look at an example. - -load_gaimc_graph('bfs_example'); -graph_draw(A,xy,'labels',labels) -full(A) -labels' - -%% -% This output means that vertex 'r' is row 1, vertex 's' is row 2 and -% because A(1,2) = 1, then there is an edge between them, just like in the -% picture. -% -% One funny property is that A(2,1) = 1 too! So we actually have to store -% each edge twice in the adjacency matrix. This might seem wasteful, but -% its hard to avoid as I've learned while working on graph algorithms. So -% don't worry about it! It also makes the generalization to directed -% graphs (below) easy. -% -% For more information about the adjacency matrix representation of a -% graph, see a standard book on graph algorithms. -% - -%% Weighted and directed graphs -% Our previous case handled the situation for undirected graphs only. To -% encode weighted and directed graphs, we use weighted and non-symmetric -% adjacency matrices. -% -% For a weighted matrix, A(i,j) = distance between i and j for most of the -% algorithms in gaimc. But A(i,j) = 0 means there is no edge, and so -% sometimes things can get a little tricky to get what you want. -% -% For a directed graph, just set A(i,j) ~= A(j,i). The adjacency matrix -% won't be symmetric, but that's what you want! -% -% To understand more, explore the examples or read up on adjacency matrices -% in graph theory books. - -%% Loading helper -% To make loading our sample graphs easy, gaimc defines it's own function -% to load graphs. - -load_gaimc_graph('dfs_example'); % loads one of our example graphs -whos - -%% -% This helps make our examples work regardless of where the current -% directory lies. - -%% Search algorithms -% The two standard graph search algorithms are depth first search and -% breadth first search. This library implements both. - -%% -% Load the example matrix from the Boost Graph Library -load_gaimc_graph('dfs_example'); -figure(1); graph_draw(A,xy,'labels',labels); - -%% -% Run a depth first search. The output records the distance to the other -% vertices, except where the vertices are not reachable starting from the -% first node A. -d=dfs(A,1) - -%% -% From this example, we see that vertices a-f are reachable from vertex a, -% but that verice g-i are not reachable. Given the of the edges, this -% makes sense. - -%% -% Let's look at breadth first search too, using a different example. -load_gaimc_graph('bfs_example'); -figure(1); clf; graph_draw(A,xy,'labels',labels); - -%% -% The breadth first search algorithm records the distance from the starting -% vertex to each vertex it visits in breadth first order. This means it -% visits all the vertices in order of their distance from the starting -% vertex. The d output records the distance, and the dt output records the -% step of the algorithm when the breadth first search saw the node. -[d dt] = bfs(A,2); -% draw the graph where the label is the "discovery time" of the vertex. -figure(1); clf; graph_draw(A,xy,'labels',num2str(dt)); - -%% -% Notice how the algorithm visits all vertices one edge away from the start -% vertex (0) before visiting those two edges away. - - -%% Shortest paths -% In the previous two examples, the distance between vertices was -% equivalent to the number of edges. Some graphs, however, have specific -% weights, such as the graph of flights between airports. We can use this -% information to build information about the _shortest path_ between two -% nodes in a network. - -% Find the minimum travel time between Los Angeles (LAX) and -% Rochester Minnesota (RST). -load_gaimc_graph('airports') -A = -A; % fix funny encoding of airport data -lax=247; rst=355; - -[d pred] = dijkstra(A,lax); % find all the shorest paths from Los Angeles. - -fprintf('Minimum time: %g\n',d(rst)); -% Print the path -fprintf('Path:\n'); -path =[]; u = rst; while (u ~= lax) path=[u path]; u=pred(u); end -fprintf('%s',labels{lax}); -for i=path; fprintf(' --> %s', labels{i}); end, fprintf('\n'); - -%% Minimum spanning trees -% A minimum spanning tree is a set of edges from a graph that ... -% -% This demo requires the mapping toolbox for maximum effect, but we'll do -% okay without it. - -%% -% Our data comes from a graph Brendan Frey prepared for his affinity -% propagation clustering tool. For 456 cities in the US, we have the mean -% travel time between airports in those cities, along with their latitude -% and longitude. -load_gaimc_graph('airports') - -%% -% For some reason, the data is stored with the negative travel time between -% cities. (I believe this is so that closer cities have larger edges -% between them.) But for a minimum spanning tree, we want the actual -% travel time between cities. -A = -A; - -%% -% Now, we just call MST and look at the result. -% T = mst_prim(A); -% This command means we can't run the demo, so it's commented out. - -%% -% Oops, travel time isn't symmetric! Let's just pick the longest possible -% time. -A = max(A,A'); -T = mst_prim(A); -sum(sum(T))/2 % total travel time in tree - -%% -% Well, the total weight isn't that helpful, let's _look_ at the data -% instead. -clf; -gplot(T,xy); - -%% -% Hey! That looks like the US! You can see regional airports and get some -% sense of the overall connectivity. - -%% Connected components -% The connected components of a network determine which parts of the -% network are reachable from other parts. One of your first questions -% about any network should generally be: is it connected? -% -% There are two types of connected components: components and strongly -% connected components. gaimc only implements an algorithm for the latter -% case, but that's okay! It turns out it computes exactly the right thing -% for connected components as well. The difference only occurs when the -% graph is undirected vs. directed. - -load_gaimc_graph('dfs_example') -graph_draw(A,xy) - -%% -% This picture shows there are 3 strongly connected components and 2 -% connected components - -% get the number of strongly connected components -max(scomponents(A)) - -%% - -% get the number of connected components -max(scomponents(A|A')) % we make the graph symmetric first by "or"ing each entry - -%% -% Let's look at the vertices in the strongly connected components -cc = scomponents(A) -%% -% The output tells us that vertices 1,2,3,5,6 are in one strong component, -% vertex 4 is it's own strong component, and vertices 7,8,9 are in another -% one. Remember that a strong component is all the vertices mutually -% reachable from a given vertex. If you start at vertex 4, you can't get -% anywhere else! That's why it is in a different component than vertices -% 1,2,3,5,6. - -%% -% We also have a largest_component function that makes it easy to just get -% the largest connected component. -clf; -[Acc,f] = largest_component(A); -graph_draw(Acc,xy(f,:)) -%% -% The filter variable f, tells us which vertices in the original graph made -% it into the largest strong component. We can just apply that filter to -% the coordinates xy and reuse them for drawing the graph! - -%% Statistics -% Graph statistics are just measures that indicate a property of the graph -% at every vertex, or at every edges. Arguably the simplest graph -% statistic would be the average vertex degree. Because such statistics -% are easy to compute with the adjaceny matrix in Matlab, they do not have -% special functions in gaimc. - -%% -% Load a road network to use for statistical computations -load_gaimc_graph('minnesota'); -gplot(A,xy); - -%% -% Average vertex degree -d = sum(A,2); -mean(d) - -%% -% So the average number of roads at any intersection is 2.5. My guess is -% that many roads have artificial intersections in the graph structure that -% do not correspond to real intersections. Try validating that hypothesis -% using the library! - -%% -% Average clustering coefficients -ccfs = clustercoeffs(A); -mean(ccfs) -% The average clustering coefficient is a measure of the edge density -% throughout the graph. A small value indicates that the network has few -% edges and they are well distributed throughout the graph. - -%% -% Average core numbers -cn = corenums(A); -mean(cn) - - - -%% Efficient repetition -% Every time a gaimc function runs, it converts the adjacency matrix into a -% set of compressed sparse row arrays. These arrays yield efficient access -% to the edges of the graph starting at a particular vertex. For many -% function calls on the same graph, this conversion process slows the -% algorithms. Hence, gaimc also accepts pre-converted input, in which case -% it skips it's conversion. -% -% Let's demonstrate how this works by calling Dijkstra's algorithm to -% compute the shortest paths between all vertices in the graph. The -% Floyd Warshall algorithm computes these same quantities more efficiently, -% but that would just be one more algorithm to implement and maintain. - -%% -% Load and convert the graph. -load_gaimc_graph('all_shortest_paths_example'); -A = spfun(@(x) x-min(min(A))+1,A); % remove the negative edges -As = convert_sparse(A); -%% -% Now, we'll run Dijkstra's algorithm for every vertex and save the result -% On my 2GHz laptop, this takes 0.000485 seconds. -n = size(A,1); -D = zeros(n,n); -tic -for i=1:n - D(i,:) = dijkstra(As,i); -end -toc -%% -% Let's try it without the conversion to see if we can notice the -% difference in speed. -% On my 2GHz laptop, this takes 0.001392 seconds. -D2 = zeros(n,n); -tic -for i=1:n - D2(i,:) = dijkstra(A,i); -end -toc -%% -% And just to check, let's make sure the output is the same. -isequal(D,D2) diff --git a/demo/performance_comparison.m b/demo/performance_comparison.m deleted file mode 100644 index f823fd5..0000000 --- a/demo/performance_comparison.m +++ /dev/null @@ -1,261 +0,0 @@ -%% Compare performance of gaimc to matlab_bgl -% While the |gaimc| library implements its graph routines in Matlab -% "m"-code, the |matlab_bgl| library uses graph algorithms from the Boost -% graph library in C++ through a mex interface. Folklore has it that -% Matlab code with for loops like those required in the |gaimc| library is -% considerably slower. This example examines this lore and shows that -% |gaimc| is typically within a factor of 2-4 of the mex code. - -%% -% *This demo is unlikely to work on your own computer.* -% *They depend on having the MatlabBGL routines in one spot.* - -%% Setup the environment -% We need MatlabBGL on the path -graphdir = '../graphs/'; -matlabbgldir = '~/dev/matlab-bgl/4.0'; -try - addpath(matlabbgldir); % change this to your matlab_bgl path - ci=components(sparse(ones(5))); -catch - error('gaimc:performance_comparison','Matlab BGL is not working, halting...'); -end - -%% -% Check to make sure we are in the correct directory -cwd = pwd; dirtail = ['gaimc' filesep 'demo']; -if strcmp(cwd(end-length(dirtail)+1:end),dirtail) == 0 - error('%s should be executed from %s\n',mfilename,dirtail); -end - -%% -% initalize the results structure -results=[]; -mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; - -%% Depth first search -% To compare these functions, we have to make a copy and then delete it -copyfile(['..' filesep 'dfs.m'],'dfstest.m'); -graphs = {'all_shortest_paths_example', 'clr-24-1', 'cs-stanford', ... - 'minnesota','tapir'}; -nrep=30; ntests=100; -fprintf(repmat(' ',1,76)); -for rep=1:nrep - % Matlab needs 1 iteration to compile the function - if nrep==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end - for gi=1:length(graphs) - load([graphdir graphs{gi} '.mat']); n=size(A,1); - At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - for ti=1:ntests - fprintf([repmat('\b',1,76) '%20s rep=%3i graph=%-30s trial=%4i'], ... - 'dfs',rep,graphs{gi},ti); - v=ceil(n*rand(1)); - tic; d1=dfs(A,v); mex_std=mex_std+toc; - tic; d2=dfs(At,v,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; d3=dfstest(A,v); mat_std=mat_std+toc; - tic; d4=dfstest(As,v); mat_fast=mat_fast+toc; - if any(d1 ~= d2) || any(d2 ~= d3) || any(d3 ~= d4) - error('gaimc:dfs','incorrect results from dijkstra'); - end - end - end -end -fprintf('\n'); -delete('dfstest.m'); -results(end+1).name='dfs'; -results(end).mex_fast = mex_fast; -results(end).mat_fast = mat_fast; -results(end).mex_std = mex_std; -results(end).mat_std = mat_std; - -%% Connected components -% To evaluate the performance of the connected components algorithm, we use -% sets of random graphs. -nrep=30; -szs=[1 10 100 5000 10000 50000]; -comp_results=[mex_fast mat_fast mex_std mat_std]; -for szi=1:length(szs) - fprintf('\n%20s size=%-5i ', 'scomponents', szs(szi)); - % Matlab needs 1 iteration to compile the function - if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end - for rep=1:nrep - fprintf('\b\b\b\b'); fprintf(' %3i', rep); - A=sprand(szs(szi),szs(szi),25/szs(szi)); - At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - tic; cc1=components(A); mex_std=mex_std+toc; - tic; cc2=components(At,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; cc3=scomponents(A); mat_std=mat_std+toc; - tic; cc4=scomponents(As); mat_fast=mat_fast+toc; - cs1=accumarray(cc1,1,[max(cc1) 1]); - cs2=accumarray(cc2,1,[max(cc2) 1]); - cs3=accumarray(cc3,1,[max(cc3) 1]); - cs4=accumarray(cc4,1,[max(cc4) 1]); - if any(cs1 ~= cs2) || any(cs2 ~= cs3) || any(cs2 ~= cs4) - error('gaimc:scomponents','incorrect results from scomponents'); - end - end - comp_results(end+1,:) = [mex_fast mat_fast mex_std mat_std]; -end -comp_results=diff(comp_results); -results(end+1).name='scomponents'; -results(end).mex_fast = mex_fast; -results(end).mat_fast = mat_fast; -results(end).mex_std = mex_std; -results(end).mat_std = mat_std; - -%% -% Plot the data for connected components -% This plot isn't too meaningful, so I've omitted it. - -% plot(szs, comp_results,'.-'); -% legend('mex fast','mat fast','mex std','mat std','Location','Northwest'); -% ylabel('time'); xlabel('graph size'); - -%% Dijkstra's algorithm -% To evaluate the performance of Dijkstra's algorithm, we pick -graphs = {'clr-25-2', 'clr-24-1', 'cs-stanford', ... - 'minnesota', 'tapir'}; -nrep=30; ntests=100; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; -for rep=1:nrep - for gi=1:length(graphs) - load([graphdir graphs{gi} '.mat']); n=size(A,1); - At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - for ti=1:ntests - fprintf([repmat('\b',1,66) '%20s rep=%3i graph=%-20s trial=%4i'], ... - 'dijkstra',rep,graphs{gi},ti); - v=ceil(n*rand(1)); - tic; d1=dijkstra_sp(A,v); mex_std=mex_std+toc; - tic; d2=dijkstra_sp(At,v,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; d3=dijkstra(A,v); mat_std=mat_std+toc; - tic; d4=dijkstra(As,v); mat_fast=mat_fast+toc; - if any(d1 ~= d2) || any(d2 ~= d3) || any(d3 ~= d4) - error('gaimc:dijkstra','incorrect results from dijkstra'); - end - end - end -end -fprintf('\n'); -results(end+1).name='dijkstra'; -results(end).mex_fast = mex_fast; -results(end).mat_fast = mat_fast; -results(end).mex_std = mex_std; -results(end).mat_std = mat_std; - -%% Directed Clustering coefficients -nrep=30; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; -comp_results=[mex_fast mat_fast mex_std mat_std]; -szs=[1 10 100 5000 10000 50000]; -for szi=1:length(szs) - fprintf('\n%20s size=%-5i ', 'dirclustercoeffs', szs(szi)); - % Matlab needs 1 iteration to compile the function - if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end - for rep=1:nrep - fprintf('\b\b\b\b'); fprintf(' %3i', rep); - A=sprand(szs(szi),szs(szi),25/szs(szi)); - At=A'; - [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - [cp ri ati]=sparse_to_csr(At); As.cp=cp; As.ri=ri; As.ati=ati; - tic; cc1=clustering_coefficients(A); mex_std=mex_std+toc; - tic; cc2=clustering_coefficients(At,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; cc3=dirclustercoeffs(A); mat_std=mat_std+toc; - tic; cc4=dirclustercoeffs(As); mat_fast=mat_fast+toc; - end - comp_results(end+1,:) = [mex_fast mat_fast mex_std mat_std]; -end -fprintf('\n'); -comp_results=diff(comp_results); -results(end+1).name='dirclustercoeffs'; -results(end).mex_fast = mex_fast; -results(end).mat_fast = mat_fast; -results(end).mex_std = mex_std; -results(end).mat_std = mat_std; - -%% Minimum spanning tree -nrep=30; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; -comp_results=[]; -szs=[10 100 5000 10000]; -for szi=1:length(szs) - fprintf('\n%20s size=%-5i ', 'mst_prim', szs(szi)); - % Matlab needs 1 iteration to compile the function - if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end - for rep=1:nrep - fprintf('\b\b\b\b'); fprintf(' %3i', rep); - A=abs(sprandsym(szs(szi),25/szs(szi))); - At=A'; - [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - tic; T1=prim_mst(A); mex_std=mex_std+toc; - tic; [t1i t1j t1v]=prim_mst(At,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; T2=mst_prim(A,0); mat_std=mat_std+toc; - tic; [t2i t2j t2v]=mst_prim(As,0); mat_fast=mat_fast+toc; - T1f=sparse(t1i,t1j,t1v,size(A,1),size(A,2)); - T2f=sparse(t2i,t2j,t2v,size(A,1),size(A,2)); - if ~isequal(T1,T2) || ~isequal(T1f+T1f',T2f+T2f') || ... - ~isequal(T1,T2f+T2f') - keyboard - warning('gaimc:mst_prim',... - 'incorrect results from mst_prim (%i,%i)', szi, rep); - end - end - comp_results(end+1,:) = [mex_fast mat_fast mex_std mat_std]; -end -fprintf('\n'); -comp_results=diff(comp_results); -results(end+1).name='mst_prim'; -results(end).mex_fast = mex_fast; -results(end).mat_fast = mat_fast; -results(end).mex_std = mex_std; -results(end).mat_std = mat_std; - -%% Undirected Clustering coefficients -nrep=30; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; -undircc_results=[mex_fast mat_fast mex_std mat_std]; -szs=[1 10 100 5000 10000 50000]; -for szi=1:length(szs) - fprintf('\n%20s size=%-5i ', 'clustercoeffs', szs(szi)); - % Matlab needs 1 iteration to compile the function - if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end - for rep=1:nrep - fprintf('\b\b\b\b'); fprintf(' %3i', rep); - A=abs(sprandsym(szs(szi),25/szs(szi))); - At=A'; - [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - tic; cc1=clustering_coefficients(A,struct('undirected',1)); mex_std=mex_std+toc; - tic; cc2=clustering_coefficients(At,struct('undirected',1,'istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; cc3=clustercoeffs(A); mat_std=mat_std+toc; - tic; cc4=clustercoeffs(As); mat_fast=mat_fast+toc; - end - undircc_results(end+1,:) = [mex_fast mat_fast mex_std mat_std]; -end -%% -fprintf('\n'); -undircc_results=diff(undircc_results); -results(end+1).name='clustercoeffs'; -results(end).mex_fast = mex_fast; -results(end).mat_fast = mat_fast; -results(end).mex_std = mex_std; -results(end).mat_std = mat_std; - -%% Summarize the results -% We are going to summarize the results in a bar plot based on the -% algorithm. Each algorithm is a single bar -graphs = {'all_shortest_paths_example', 'clr-24-1', 'cs-stanford', ... - 'minnesota','tapir'};, where the performance of the -% mex code is 1. -nresults=length(results); -Ystd = zeros(nresults,1); -Yfast = zeros(nresults,1); -for i=1:nresults - Ystd(i)=results(i).mat_std/results(i).mex_std; - Yfast(i)=results(i).mat_fast/results(i).mex_fast; -end -bar(1:nresults,[Ystd Yfast]); set(gca,'XTickLabel',{results.name}); -legend('Standard','Fast','Location','Northwest'); - - - diff --git a/demo/performance_comparison_simple.m b/demo/performance_comparison_simple.m deleted file mode 100644 index 54c8905..0000000 --- a/demo/performance_comparison_simple.m +++ /dev/null @@ -1,131 +0,0 @@ -%% Compare performance of gaimc to matlab_bgl -% While the |gaimc| library implements its graph routines in Matlab -% "m"-code, the |matlab_bgl| library uses graph algorithms from the Boost -% graph library in C++ through a mex interface. Folklore has it that -% Matlab code with for loops like those required in the |gaimc| library is -% considerably slower. This example examines this lore and shows that -% |gaimc| is typically within a factor of 2-4 of the mex code for one -% function. The full performance_comparison suite evaluates the rest, but -% it takes a while to run, so I only run it when I'm interested in a -% complete picture and can afford to run it overnight. - -%% -% *This demo is unlikely to work on your own computer.* -% *They depend on having the MatlabBGL routines in one spot.* - -%% Setup the environment -% We need MatlabBGL on the path -graphdir = '../graphs/'; -matlabbgldir = '~/dev/matlab-bgl/4.0'; -try - addpath(matlabbgldir); % change this to your matlab_bgl path - ci=components(sparse(ones(5))); -catch - error('gaimc:performance_comparison','Matlab BGL is not working, halting...'); -end - -%% -% Check to make sure we are in the correct directory -cwd = pwd; dirtail = ['gaimc' filesep 'demo']; -if strcmp(cwd(end-length(dirtail)+1:end),dirtail) == 0 - error('%s should be executed from %s\n',mfilename,dirtail); -end - -%% -% initalize the results structure -results=[]; -mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; - -%% Connected components -% To evaluate the performance of the connected components algorithm, we use -% sets of random graphs. -nrep=30; -szs=[1 10 100 5000 10000 50000]; -comp_results=[mex_fast mat_fast mex_std mat_std]; -for szi=1:length(szs) - fprintf('\n%20s size=%-5i ', 'scomponents', szs(szi)); - % Matlab needs 1 iteration to compile the function - if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end - for rep=1:nrep - fprintf('\b\b\b\b'); fprintf(' %3i', rep); - A=sprand(szs(szi),szs(szi),25/szs(szi)); - At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - tic; cc1=components(A); mex_std=mex_std+toc; - tic; cc2=components(At,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; cc3=scomponents(A); mat_std=mat_std+toc; - tic; cc4=scomponents(As); mat_fast=mat_fast+toc; - cs1=accumarray(cc1,1,[max(cc1) 1]); - cs2=accumarray(cc2,1,[max(cc2) 1]); - cs3=accumarray(cc3,1,[max(cc3) 1]); - cs4=accumarray(cc4,1,[max(cc4) 1]); - if any(cs1 ~= cs2) || any(cs2 ~= cs3) || any(cs2 ~= cs4) - error('gaimc:scomponents','incorrect results from scomponents'); - end - end - comp_results(end+1,:) = [mex_fast mat_fast mex_std mat_std]; -end -comp_results=diff(comp_results); -results(end+1).name='scomponents'; -results(end).mex_fast = mex_fast; -results(end).mat_fast = mat_fast; -results(end).mex_std = mex_std; -results(end).mat_std = mat_std; - -%% Dijkstra's algorithm -% To evaluate the performance of Dijkstra's algorithm, we pick -graphs = {'clr-25-2', 'clr-24-1', 'cs-stanford', ... - 'minnesota', 'tapir'}; -nrep=30; ntests=100; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; -for rep=1:nrep - for gi=1:length(graphs) - load([graphdir graphs{gi} '.mat']); n=size(A,1); - At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - for ti=1:ntests - fprintf([repmat('\b',1,66) '%20s rep=%3i graph=%-20s trial=%4i'], ... - 'dijkstra',rep,graphs{gi},ti); - v=ceil(n*rand(1)); - tic; d1=dijkstra_sp(A,v); mex_std=mex_std+toc; - tic; d2=dijkstra_sp(At,v,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; d3=dijkstra(A,v); mat_std=mat_std+toc; - tic; d4=dijkstra(As,v); mat_fast=mat_fast+toc; - if any(d1 ~= d2) || any(d2 ~= d3) || any(d3 ~= d4) - error('gaimc:dijkstra','incorrect results from dijkstra'); - end - end - end -end -fprintf('\n'); -results(end+1).name='dijkstra'; -results(end).mex_fast = mex_fast; -results(end).mat_fast = mat_fast; -results(end).mex_std = mex_std; -results(end).mat_std = mat_std; - - -%% Summarize the results -% We are going to summarize the results in a bar plot based on the -% algorithm. Each algorithm is a single bar where the performance of the -% mex code is 1. -nresults=length(results); -Ystd = zeros(nresults,1); -Yfast = zeros(nresults,1); -for i=1:nresults - Ystd(i)=results(i).mat_std/results(i).mex_std; - Yfast(i)=results(i).mat_fast/results(i).mex_fast; -end -bar(1:nresults,[Ystd Yfast]); set(gca,'XTickLabel',{results.name}); -legend('Standard','Fast','Location','Northwest'); - -%% -% From this, we see that the connected component codes are about half the -% speed of the Matlab BGL functions and Dijkstra's is about 1/4th the -% speed. This seems unideal, but the code is much more portable and -% flexible. -% -% The difference between the Standard and Fast results is that the fast -% results eliminate all data translation for both gaimc and MatlabBGL and -% are evaluating the actual algorithm implementation and not any of the -% data translation components. - diff --git a/dfs.m b/dfs.m deleted file mode 100644 index 53a8367..0000000 --- a/dfs.m +++ /dev/null @@ -1,94 +0,0 @@ -function [d dt ft pred] = dfs(A,u,full,target) -% DFS Compute depth first search distances, times, and tree for a graph -% -% [d dt ft pred] = dfs(A,u) returns the distance (d), the discover (dt) and -% finish time (ft) for each vertex in the graph in a depth first search -% starting from vertex u. -% d = dt(i) = ft(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. -% -% [...] = dfs(A,u,1) continues the dfs for all components of the graph, not -% just the connected component with u -% [...] = dfs(A,u,[],v) stops the dfs when it hits the vertex v -% -% Note 1: When target is specified, the finish time array only records the -% finish time for all vertices that actually finished. The array will then -% be undefined on a significant portion of the graph, but that doesn't -% indicate the vertices are unreachable; they just haven't been reached -% yet. -% -% Example: -% load_gaimc_graph('dfs_example.mat') % use the dfs example from Boost -% d = dfs(A,1) -% -% See also BFS - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2008-04-10: Initial coding - -if ~exist('full','var') || isempty(full), full=0; end -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); ft=-1*ones(n,1); pred=zeros(1,n); -rs=zeros(2*n,1); rss=0; % recursion stack holds two nums (v,ri) - -% start dfs at u -t=0; targethit=0; -for i=1:n - if i==1, v=u; - else v=mod(u+i-1,n)+1; if d(v)>0, continue; end, end - d(v)=0; dt(v)=t; t=t+1; ri=rp(v); - rss=rss+1; rs(2*rss-1)=v; rs(2*rss)=ri; % add v to the stack - while rss>0 - v=rs(2*rss-1); ri=rs(2*rss); rss=rss-1; % pop v from the stack - if v==target || targethit, - ri=rp(v+1); targethit=1; % end the algorithm if v is the target - end - while ri %s', labels{i}); end, fprintf('\n'); - -% David F. Gleich -% Copyright, Stanford University, 2008-200909 - -% History -% 2008-04-09: Initial coding -% 2009-05-15: Documentation - -if isstruct(A), - rp=A.rp; ci=A.ci; ai=A.ai; - check=0; -else - [rp ci ai]=sparse_to_csr(A); check=1; -end -if check && any(ai)<0, error('gaimc:dijkstra', ... - 'dijkstra''s algorithm cannot handle negative edge weights.'); end - -n=length(rp)-1; -d=Inf*ones(n,1); T=zeros(n,1); L=zeros(n,1); -pred=zeros(1,length(rp)-1); - -n=1; T(n)=u; L(u)=n; % oops, n is now the size of the heap - -% enter the main dijkstra loop -d(u) = 0; -while n>0 - v=T(1); ntop=T(n); T(1)=ntop; L(ntop)=1; n=n-1; % pop the head off the heap - k=1; kt=ntop; % move element T(1) down the heap - while 1, - i=2*k; - if i>n, break; end % end of heap - if i==n, it=T(i); % only one child, so skip - else % pick the smallest child - lc=T(i); rc=T(i+1); it=lc; - if d(rc)d(v)+ew - d(w)=d(v)+ew; pred(w)=v; - % check if w is in the heap - k=L(w); onlyup=0; - if k==0 - % element not in heap, only move the element up the heap - n=n+1; T(n)=w; L(w)=n; k=n; kt=w; onlyup=1; - else kt=T(k); - end - % update the heap, move the element down in the heap - while 1 && ~onlyup, - i=2*k; - if i>n, break; end % end of heap - if i==n, it=T(i); % only one child, so skip - else % pick the smallest child - lc=T(i); rc=T(i+1); it=lc; - if d(rc)1, % j==1 => element at top of heap - j2=floor(j/2); tj2=T(j2); % parent element - if d(tj2)1, cccyc=zeros(n,1); end -if nargout>2, ccmid=zeros(n,1); end -if nargout>3, ccin=zeros(n,1); end -if nargout>4, ccout=zeros(n,1); end -if nargout>5, nf=zeros(n,1); end -% precompute degrees -for v=1:n, - for rpi=rp(v):rp(v+1)-1, - w=ci(rpi); - if v==w, continue; else degs(w)=degs(w)+1; degs(v)=degs(v)+1; end - end -end -ew=1; ew2=1; -for v=1:n - % setup counts for the different cycle types - bilatedges=0; curcccyc=0; curccmid=0; curccin=0; curccout=0; - % 1. - % find triangles with out links as last step, so precompute the inlinks - % back to node v - for cpi=cp(v):cp(v+1)-1 - w=ri(cpi); if usew, ew=ati(cpi); end - if v~=w, ind(w)=1; cache(w)=ew^(1/3); end - end - % count cycles (cycles are out->out->out) - for rpi=rp(v):rp(v+1)-1 - w=ci(rpi); if v==w, continue; end % discount self-loop - if usew, ew=ai(rpi)^(1/3); end - for rpi2=rp(w):rp(w+1)-1 - x=ci(rpi2); if x==w, continue; end - if x==v, bilatedges=bilatedges+1; continue; end - if ind(x) - if usew, ew2=ai(rpi2); end - curcccyc=curcccyc+ew*ew2^(1/3)*cache(x); - end - end - end - % count middle-man circuits (out->in->out) - for rpi=rp(v):rp(v+1)-1 - w=ci(rpi); if v==w, continue; end % discount self-loop - if usew, ew=ai(rpi)^(1/3); end - for cpi=cp(w):cp(w+1)-1 - x=ri(cpi); if x==w, continue; end - if ind(x) - if usew, ew2=ati(cpi); end - curccmid=curccmid+ew*ew2^(1/3)*cache(x); - end - end - end - % count in-link circuits (in->out->out) - for cpi=cp(v):cp(v+1)-1 - w=ri(cpi); if v==w, continue; end % discount self-loop - if usew, ew=ati(cpi)^(1/3); end - for rpi2=rp(w):rp(w+1)-1 - x=ci(rpi2); if x==w, continue; end - if ind(x) - if usew, ew2=ai(rpi2); end - curccin=curccin+ew*ew2^(1/3)*cache(x); - end - end - end - % reset and reinit the cache for outlinks - for cpi=cp(v):cp(v+1)-1, w=ri(cpi); ind(w)=0; end % reset indicator - for rpi=rp(v):rp(v+1)-1 - w=ci(rpi); if usew, ew=ai(rpi); end - if v~=w, ind(w)=1; cache(w)=ew^(1/3); end - end - % count out-link circuits (out->out->in) - for rpi=rp(v):rp(v+1)-1 - w=ci(rpi); if v==w, continue; end % discount self-loop - if usew, ew=ai(rpi)^(1/3); end - for rpi2=rp(w):rp(w+1)-1 - x=ci(rpi2); if x==w, continue; end - if ind(x) - if usew, ew2=ai(rpi2); end - curccout=curccout+ew*ew2^(1/3)*cache(x); - end - end - end - for rpi=rp(v):rp(v+1)-1, w=ci(rpi); ind(w)=0; end % reset indicator - % store the values - nf=degs(v)*(degs(v)-1) - 2*bilatedges; - curcc=curcccyc+curccmid+curccin+curccout; - if nf>0 && donorm, curcc=curcc/nf; end - cc(v)=curcc; - if nargout>1, cccyc(v)=curcccyc; end - if nargout>2, ccmid(v)=curccmid; end - if nargout>3, ccin(v)=curccin; end - if nargout>4, ccout(v)=curccout; end - if nargout>5, nf(v)=nf; end -end diff --git a/floydwarshall.m b/floydwarshall.m deleted file mode 100644 index 376d2e2..0000000 --- a/floydwarshall.m +++ /dev/null @@ -1,78 +0,0 @@ -function [D,P] = floydwarshall(A) -% FLOYDWARSHALL Compute all shortest paths using the Floyd-Warshall algorithm -% -% D=floydwarshall(A) returns the shortest distance matrix between all pairs -% of nodes in the graph A. If A has a negative weight cycle, then this -% algorithm will throw an error. -% -% [D,P]=floydwarshall(A) also returns the matrix of predecessors. -% -% See also DIJKSTRA -% -% Example: -% load_gaimc_graph('all_shortest_paths_example') -% D = floydwarshall(A) - -% David F. Gleich -% Copyright, Stanford University, 2009 - -% History -% 2008-05-23: Initial coding - -if isstruct(A), - rp=A.rp; ci=A.ci; ai=A.ai; - [ri ci ai]=csr_to_sparse(rp,ci,ai); - n = length(rp)-1; -else - [ri ci ai]=find(A); - n = size(A,1); -end - -nz = length(ai); -computeP = nargout>1; -D = Inf*ones(n,n); - -if computeP - P = zeros(n,n); - % initialize the distance and predecessor matrix - for ei=1:nz - i=ri(ei); - j=ci(ei); - v=ai(ei); - if v -% 'shapes' - 1 if node is a box, 0 if oval -% -% h = graph_draw(A,xy,...) returns a handle for each object. h(i,1) is -% the text handle for vertex i, and h(i,2) is the circle handle for -% vertex i. -% -% Originally written by Erik A. Johnson, Ali Taylan Cemgil, and Leon Peskin -% Modified by David F. Gleich for gaimc package. -% -% See also GPLOT -% -% Example: -% load_gaimc_graph('dfs_example'); -% graph_draw(A,xy); - - -% 2009-02-26 interface modified by David Gleich -% to remove automatic layout -% 2009-05-15: Added example - -% 24 Feb 2004 cleaned up, optimized and corrected by Leon Peshkin pesha @ ai.mit.edu -% Apr-2000 draw_graph Ali Taylan Cemgil -% 1995-1997 arrow Erik A. Johnson - -linestyle = '-'; % -- -. -linewidth = .5; % 2 -linecolor = 'Black'; % Red -fontsize = 8; -N = size(adj,1); -color = ones(N, 3); % colors of elipses around text -labels = cellstr(int2str((1:N)')); % labels = cellstr(char(zeros(N,1)+double('+'))); -node_t = zeros(N,1); % -for i = 1:2:length(varargin) % get optional args - switch varargin{i} - case 'linestyle', linestyle = varargin{i+1}; - case 'linewidth', linewidth = varargin{i+1}; - case 'linecolor', linecolor = varargin{i+1}; - case 'labels', labels = varargin{i+1}; - case 'fontsize', fontsize = varargin{i+1}; - case 'shapes', node_t = varargin{i+1}; node_t = node_t(:); - end -end - -x = xy(:,1); -x = x - min(x); -y = xy(:,2); -y = y - min(y); -% scale the graph so it's between 0 and 1 -xrange = max(x); -yrange = max(y); -scalefactor = max(xrange,yrange); -x = x/scalefactor; -y = y/scalefactor; - -lp_ndx = find(diag(adj)); % recover from self-loops = diagonal ones -color(lp_ndx,:) = repmat([.8 .8 .8],length(lp_ndx),1); % makes self-looped nodes blue -adj = adj - diag(diag(adj)); % clean up the diagonal - -axis([-0.1 1.1 -0.1 1.1]); -axis off; -set(gcf,'Color',[1 1 1]); -set(gca,'XTick',[], 'YTick',[], 'box','on'); % axis('square'); %colormap(flipud(gray)); - -idx1 = find(node_t == 0); wd1 = []; % Draw nodes -if ~isempty(idx1), - [h1 wd1] = textoval(x(idx1), y(idx1), labels(idx1), fontsize, color); -end; - -idx2 = find(node_t ~= 0); wd2 = []; -if ~isempty(idx2), - [h2 wd2] = textbox(x(idx2), y(idx2), labels(idx2), color); -end; - -wd = zeros(size(wd1,1) + size(wd2,1),2); -if ~isempty(idx1), wd(idx1, :) = wd1; end; -if ~isempty(idx2), wd(idx2, :) = wd2; end; - -for node = 1:N % Draw edges - edges = find(adj(node,:) == 1); - for node2 = edges - sign = 1; - if ((x(node2) - x(node)) == 0) - if (y(node) > y(node2)), alpha = -pi/2; else alpha = pi/2; end; - else - alpha = atan((y(node2)-y(node))/(x(node2)-x(node))); - if (x(node2) <= x(node)), sign = -1; end; - end; - dy1 = sign.*wd(node,2).*sin(alpha); dx1 = sign.*wd(node,1).*cos(alpha); - dy2 = sign.*wd(node2,2).*sin(alpha); dx2 = sign.*wd(node2,1).*cos(alpha); - if (adj(node2,node) == 0) % if directed edge - my_arrow([x(node)+dx1 y(node)+dy1], [x(node2)-dx2 y(node2)-dy2]); - else - line([x(node)+dx1 x(node2)-dx2], [y(node)+dy1 y(node2)-dy2], ... - 'Color', linecolor, 'LineStyle', linestyle, 'LineWidth', linewidth); - adj(node2,node) = -1; % Prevent drawing lines twice - end; - end; -end; - -if nargout > 2 - h = zeros(length(wd),2); - if ~isempty(idx1), h(idx1,:) = h1; end; - if ~isempty(idx2), h(idx2,:) = h2; end; -end; - -function [t, wd] = textoval(x, y, str, fontsize, c) -% [t, wd] = textoval(x, y, str, fontsize) Draws an oval around text objects -% INPUT: x, y - Coordinates -% str - Strings -% c - colors -% OUTPUT: t - Object Handles -% width - x and y width of ovals -if ~isa(str,'cell'), str = cellstr(str); end; -N = length(str); -wd = zeros(N,2); -temp = zeros(N,2); -for i = 1:N, - tx = text(x(i),y(i),str{i},'HorizontalAlignment','center','VerticalAlign','middle','FontSize', fontsize); - sz = get(tx, 'Extent'); - wy = sz(4); - wx = max(2/3*sz(3), wy); - wx = 0.9 * wx; % might want to play with this .9 and .5 coefficients - wy = 0.5 * wy; - ptc = ellipse(x(i), y(i), wx, wy, c(i,:)); - set(ptc, 'FaceColor', c(i,:)); % 'w' - wd(i,:) = [wx wy]; - delete(tx); - tx = text(x(i),y(i),str{i},'HorizontalAlignment','center','VerticalAlign','middle', 'FontSize', fontsize); - temp(i,:) = [tx ptc]; -end; -t = temp; - -function [p] = ellipse(x, y, rx, ry, c) -% [p] = ellipse(x, y, rx, ry) Draws Ellipse shaped patch objects -% INPUT: x,y - N x 1 vectors of x and y coordinates -% Rx, Ry - Radii -% C - colors -% OUTPUT: p - Handles of Ellipse shaped path objects - - if length(rx)== 1, rx = ones(size(x)).*rx; end; - if length(ry)== 1, ry = ones(size(x)).*ry; end; -N = length(x); -p = zeros(size(x)); -t = 0:pi/30:2*pi; -for i = 1:N - px = rx(i) * cos(t) + x(i); py = ry(i) * sin(t) + y(i); - p(i) = patch(px, py, c(i,:)); -end; - -function [h, wd] = textbox(x,y,str,c) -% [h, wd] = textbox(x,y,str) draws a box around the text -% INPUT: x, y - Coordinates -% str - Strings -% OUTPUT: h - Object Handles -% wd - x and y Width of boxes - -if ~isa(str,'cell'), str=cellstr(str); end -N = length(str); -wd = zeros(N,2); -h = zeros(N,2); -for i = 1:N, - tx = text(x(i),y(i),str{i},'HorizontalAlignment','center','VerticalAlign','middle'); - sz = get(tx, 'Extent'); - wy = 2/3 * sz(4); wyB = y(i) - wy; wyT = y(i) + wy; - wx = max(2/3 * sz(3), wy); wxL = x(i) - wx; wxR = x(i) + wx; - ptc = patch([wxL wxR wxR wxL], [wyT wyT wyB wyB], c(i,:)); - set(ptc, 'FaceColor', c(i,:)); % 'w' - wd(i,:) = [wx wy]; - delete(tx); - tx = text(x(i),y(i),str{i},'HorizontalAlignment','center','VerticalAlign','middle'); - h(i,:) = [tx ptc]; -end; - -function [h,yy,zz] = my_arrow(varargin) -% [h,yy,zz] = my_arrow(varargin) Draw a line with an arrowhead. - -% A lot of the original code is removed and most of the remaining can probably go too -% since it comes from a general use function only being called inone context. - Leon Peshkin -% Copyright 1997, Erik A. Johnson , 8/14/97 - -ax = []; % set values to empty matrices -deflen = 12; % 16 -defbaseangle = 45; % 90 -deftipangle = 16; -defwid = 0; defpage = 0; defends = 1; -ArrowTag = 'Arrow'; % The 'Tag' we'll put on our arrows -start = varargin{1}; % fill empty arguments -stop = varargin{2}; -crossdir = [NaN NaN NaN]; -len = NaN; baseangle = NaN; tipangle = NaN; wid = NaN; -page = 0; ends = NaN; -start = [start NaN]; stop = [stop NaN]; -o = 1; % expand single-column arguments -ax = gca; -% set up the UserData data (here so not corrupted by log10's and such) -ud = [start stop len baseangle tipangle wid page crossdir ends]; -% Get axes limits, range, min; correct for aspect ratio and log scale -axm = zeros(3,1); axr = axm; axrev = axm; ap = zeros(2,1); -xyzlog = axm; limmin = ap; limrange = ap; oldaxlims = zeros(1,7); -oneax = 1; % all(ax==ax(1)); LPM -if (oneax), - T = zeros(4,4); invT = zeros(4,4); -else - T = zeros(16,1); invT = zeros(16,1); -end -axnotdone = 1; % logical(ones(size(ax))); LPM -while (any(axnotdone)) - ii = 1; % LPM min(find(axnotdone)); - curax = ax(ii); - curpage = page(ii); - % get axes limits and aspect ratio - axl = [get(curax,'XLim'); get(curax,'YLim'); get(curax,'ZLim')]; - oldaxlims(find(oldaxlims(:,1)==0, 1),:) = [curax reshape(axl',1,6)]; - % get axes size in pixels (points) - u = get(curax,'Units'); - axposoldunits = get(curax,'Position'); - really_curpage = curpage & strcmp(u,'normalized'); - if (really_curpage) - curfig = get(curax,'Parent'); pu = get(curfig,'PaperUnits'); - set(curfig,'PaperUnits','points'); pp = get(curfig,'PaperPosition'); - set(curfig,'PaperUnits',pu); set(curax,'Units','pixels'); - curapscreen = get(curax,'Position'); set(curax,'Units','normalized'); - curap = pp.*get(curax,'Position'); - else - set(curax,'Units','pixels'); - curapscreen = get(curax,'Position'); - curap = curapscreen; - end - set(curax,'Units',u); set(curax,'Position',axposoldunits); - % handle non-stretched axes position - str_stretch = {'DataAspectRatioMode'; 'PlotBoxAspectRatioMode' ; 'CameraViewAngleMode' }; - str_camera = {'CameraPositionMode' ; 'CameraTargetMode' ; ... - 'CameraViewAngleMode' ; 'CameraUpVectorMode'}; - notstretched = strcmp(get(curax,str_stretch),'manual'); - manualcamera = strcmp(get(curax,str_camera),'manual'); - if ~arrow_WarpToFill(notstretched,manualcamera,curax) - % find the true pixel size of the actual axes - texttmp = text(axl(1,[1 2 2 1 1 2 2 1]), ... - axl(2,[1 1 2 2 1 1 2 2]), axl(3,[1 1 1 1 2 2 2 2]),''); - set(texttmp,'Units','points'); - textpos = get(texttmp,'Position'); - delete(texttmp); - textpos = cat(1,textpos{:}); - textpos = max(textpos(:,1:2)) - min(textpos(:,1:2)); - % adjust the axes position - if (really_curpage) % adjust to printed size - textpos = textpos * min(curap(3:4)./textpos); - curap = [curap(1:2)+(curap(3:4)-textpos)/2 textpos]; - else % adjust for pixel roundoff - textpos = textpos * min(curapscreen(3:4)./textpos); - curap = [curap(1:2)+(curap(3:4)-textpos)/2 textpos]; - end - end - % adjust limits for log scale on axes - curxyzlog = [strcmp(get(curax,'XScale'),'log'); ... - strcmp(get(curax,'YScale'),'log'); strcmp(get(curax,'ZScale'),'log')]; - if (any(curxyzlog)) - ii = find([curxyzlog;curxyzlog]); - if (any(axl(ii)<=0)) - error([upper(mfilename) ' does not support non-positive limits on log-scaled axes.']); - else - axl(ii) = log10(axl(ii)); - end - end - % correct for 'reverse' direction on axes; - curreverse = [strcmp(get(curax,'XDir'),'reverse'); ... - strcmp(get(curax,'YDir'),'reverse'); strcmp(get(curax,'ZDir'),'reverse')]; - ii = find(curreverse); - if ~isempty(ii) - axl(ii,[1 2])=-axl(ii,[2 1]); - end - % compute the range of 2-D values - curT = get(curax,'Xform'); - lim = curT*[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1;0 0 0 0 1 1 1 1;1 1 1 1 1 1 1 1]; - lim = lim(1:2,:)./([1;1]*lim(4,:)); - curlimmin = min(lim,[],2); - curlimrange = max(lim,[],2) - curlimmin; - curinvT = inv(curT); - if ~oneax - curT = curT.'; curinvT = curinvT.'; curT = curT(:); curinvT = curinvT(:); - end - % check which arrows to which cur corresponds - ii = find((ax==curax)&(page==curpage)); - oo = ones(1,length(ii)); axr(:,ii) = diff(axl,1,2) * oo; - axm(:,ii) = axl(:,1) * oo; axrev(:,ii) = curreverse * oo; - ap(:,ii) = curap(3:4)' * oo; xyzlog(:,ii) = curxyzlog * oo; - limmin(:,ii) = curlimmin * oo; limrange(:,ii) = curlimrange * oo; - if (oneax), - T = curT; invT = curinvT; - else - T(:,ii) = curT * oo; invT(:,ii) = curinvT * oo; - end; - axnotdone(ii) = zeros(1,length(ii)); -end; -oldaxlims(oldaxlims(:,1)==0,:) = []; - -% correct for log scales -curxyzlog = xyzlog.'; ii = find(curxyzlog(:)); -if ~isempty(ii) - start(ii) = real(log10(start(ii))); stop(ii) = real(log10(stop(ii))); - if (all(imag(crossdir)==0)) % pulled (ii) subscript on crossdir, 12/5/96 eaj - crossdir(ii) = real(log10(crossdir(ii))); - end -end -ii = find(axrev.'); % correct for reverse directions -if ~isempty(ii) - start(ii) = -start(ii); stop(ii) = -stop(ii); crossdir(ii) = -crossdir(ii); -end -start = start.'; stop = stop.'; % transpose start/stop values -% take care of defaults, page was done above -ii = find(isnan(start(:))); if ~isempty(ii), start(ii) = axm(ii)+axr(ii)/2; end; -ii = find(isnan(stop(:))); if ~isempty(ii), stop(ii) = axm(ii)+axr(ii)/2; end; -ii = find(isnan(crossdir(:))); if ~isempty(ii), crossdir(ii) = zeros(length(ii),1); end; -ii = find(isnan(len)); if ~isempty(ii), len(ii) = ones(length(ii),1)*deflen; end; -baseangle(ii) = ones(length(ii),1)*defbaseangle; tipangle(ii) = ones(length(ii),1)*deftipangle; -wid(ii) = ones(length(ii),1) * defwid; ends(ii) = ones(length(ii),1) * defends; -% transpose rest of values -len = len.'; baseangle = baseangle.'; tipangle = tipangle.'; wid = wid.'; -page = page.'; crossdir = crossdir.'; ends = ends.'; ax = ax.'; - -% for all points with start==stop, start=stop-(verysmallvalue)*(up-direction); -ii = find(all(start==stop)); -if ~isempty(ii) - % find an arrowdir vertical on screen and perpendicular to viewer - % transform to 2-D - tmp1 = [(stop(:,ii)-axm(:,ii))./axr(:,ii);ones(1,length(ii))]; - if (oneax), twoD=T*tmp1; - else tmp1=[tmp1;tmp1;tmp1;tmp1]; tmp1=T(:,ii).*tmp1; - tmp2=zeros(4,4*length(ii)); tmp2(:)=tmp1(:); - twoD=zeros(4,length(ii)); twoD(:)=sum(tmp2)'; - end - twoD=twoD./(ones(4,1)*twoD(4,:)); - % move the start point down just slightly - tmp1 = twoD + [0;-1/1000;0;0]*(limrange(2,ii)./ap(2,ii)); - % transform back to 3-D - if (oneax), threeD=invT*tmp1; - else tmp1=[tmp1;tmp1;tmp1;tmp1]; tmp1=invT(:,ii).*tmp1; - tmp2=zeros(4,4*length(ii)); tmp2(:)=tmp1(:); - threeD=zeros(4,length(ii)); threeD(:)=sum(tmp2)'; - end - start(:,ii) = (threeD(1:3,:)./(ones(3,1)*threeD(4,:))).*axr(:,ii)+axm(:,ii); -end; -% compute along-arrow points -% transform Start points - tmp1 = [(start-axm)./axr; 1]; - if (oneax), X0=T*tmp1; - else tmp1 = [tmp1;tmp1;tmp1;tmp1]; tmp1=T.*tmp1; - tmp2 = zeros(4,4); tmp2(:)=tmp1(:); - X0=zeros(4,1); X0(:)=sum(tmp2)'; - end - X0=X0./(ones(4,1)*X0(4,:)); -% transform Stop points - tmp1=[(stop-axm)./axr; 1]; - if (oneax), Xf=T*tmp1; - else tmp1=[tmp1;tmp1;tmp1;tmp1]; tmp1=T.*tmp1; - tmp2=zeros(4,4); tmp2(:)=tmp1(:); - Xf=zeros(4,1); Xf(:)=sum(tmp2)'; - end - Xf=Xf./(ones(4,1)*Xf(4,:)); -% compute pixel distance between points - D = sqrt(sum(((Xf(1:2,:)-X0(1:2,:)).*(ap./limrange)).^2)); -% compute and modify along-arrow distances - len1 = len; - len2 = len - (len.*tan(tipangle/180*pi)-wid/2).*tan((90-baseangle)/180*pi); - slen0 = 0; slen1 = len1 .* ((ends==2)|(ends==3)); - slen2 = len2 .* ((ends==2)|(ends==3)); - len0 = 0; len1 = len1 .* ((ends==1)|(ends==3)); - len2 = len2 .* ((ends==1)|(ends==3)); - ii = find((ends==1)&(D -% check if we are in "WarpToFill" mode. - out = strcmp(get(curax,'WarpToFill'),'on'); - % 'WarpToFill' is undocumented, so may need to replace this by - % out = ~( any(notstretched) & any(manualcamera) ); diff --git a/graphs/airports.mat b/graphs/airports.mat deleted file mode 100644 index 6134dd0..0000000 Binary files a/graphs/airports.mat and /dev/null differ diff --git a/graphs/all_shortest_paths_example.mat b/graphs/all_shortest_paths_example.mat deleted file mode 100644 index edc728c..0000000 Binary files a/graphs/all_shortest_paths_example.mat and /dev/null differ diff --git a/graphs/bfs_example.mat b/graphs/bfs_example.mat deleted file mode 100644 index 7ee9b0b..0000000 Binary files a/graphs/bfs_example.mat and /dev/null differ diff --git a/graphs/bgl_cities.mat b/graphs/bgl_cities.mat deleted file mode 100644 index 819df81..0000000 Binary files a/graphs/bgl_cities.mat and /dev/null differ diff --git a/graphs/celegans.mat b/graphs/celegans.mat deleted file mode 100644 index 79d1b29..0000000 Binary files a/graphs/celegans.mat and /dev/null differ diff --git a/graphs/clique-10.mat b/graphs/clique-10.mat deleted file mode 100644 index 1b878b7..0000000 Binary files a/graphs/clique-10.mat and /dev/null differ diff --git a/graphs/clr-24-1.mat b/graphs/clr-24-1.mat deleted file mode 100644 index d01040f..0000000 Binary files a/graphs/clr-24-1.mat and /dev/null differ diff --git a/graphs/clr-25-2.mat b/graphs/clr-25-2.mat deleted file mode 100644 index a843a0c..0000000 Binary files a/graphs/clr-25-2.mat and /dev/null differ diff --git a/graphs/clr-26-1.mat b/graphs/clr-26-1.mat deleted file mode 100644 index edc728c..0000000 Binary files a/graphs/clr-26-1.mat and /dev/null differ diff --git a/graphs/clr-27-1.mat b/graphs/clr-27-1.mat deleted file mode 100644 index da42e0a..0000000 Binary files a/graphs/clr-27-1.mat and /dev/null differ diff --git a/graphs/cores_example.mat b/graphs/cores_example.mat deleted file mode 100644 index 4c0b964..0000000 Binary files a/graphs/cores_example.mat and /dev/null differ diff --git a/graphs/cs-stanford.mat b/graphs/cs-stanford.mat deleted file mode 100644 index 5a8955c..0000000 Binary files a/graphs/cs-stanford.mat and /dev/null differ diff --git a/graphs/dfs_example.mat b/graphs/dfs_example.mat deleted file mode 100644 index 2698814..0000000 Binary files a/graphs/dfs_example.mat and /dev/null differ diff --git a/graphs/dominator_tree_example.mat b/graphs/dominator_tree_example.mat deleted file mode 100644 index d8965ad..0000000 Binary files a/graphs/dominator_tree_example.mat and /dev/null differ diff --git a/graphs/kt-3-2.mat b/graphs/kt-3-2.mat deleted file mode 100644 index 7e77903..0000000 Binary files a/graphs/kt-3-2.mat and /dev/null differ diff --git a/graphs/kt-3-7.mat b/graphs/kt-3-7.mat deleted file mode 100644 index acb6c39..0000000 Binary files a/graphs/kt-3-7.mat and /dev/null differ diff --git a/graphs/kt-6-23.mat b/graphs/kt-6-23.mat deleted file mode 100644 index 04e752a..0000000 Binary files a/graphs/kt-6-23.mat and /dev/null differ diff --git a/graphs/kt-7-2.mat b/graphs/kt-7-2.mat deleted file mode 100644 index 0a498b5..0000000 Binary files a/graphs/kt-7-2.mat and /dev/null differ diff --git a/graphs/line-7.mat b/graphs/line-7.mat deleted file mode 100644 index b48a610..0000000 Binary files a/graphs/line-7.mat and /dev/null differ diff --git a/graphs/matching_example.mat b/graphs/matching_example.mat deleted file mode 100644 index df73c5c..0000000 Binary files a/graphs/matching_example.mat and /dev/null differ diff --git a/graphs/max_flow_example.mat b/graphs/max_flow_example.mat deleted file mode 100644 index b11b2e0..0000000 Binary files a/graphs/max_flow_example.mat and /dev/null differ diff --git a/graphs/minnesota.mat b/graphs/minnesota.mat deleted file mode 100644 index 708c960..0000000 Binary files a/graphs/minnesota.mat and /dev/null differ diff --git a/graphs/padgett-florentine.mat b/graphs/padgett-florentine.mat deleted file mode 100644 index c502886..0000000 Binary files a/graphs/padgett-florentine.mat and /dev/null differ diff --git a/graphs/tapir.mat b/graphs/tapir.mat deleted file mode 100644 index 6d20d72..0000000 Binary files a/graphs/tapir.mat and /dev/null differ diff --git a/graphs/tarjan-biconn.mat b/graphs/tarjan-biconn.mat deleted file mode 100644 index 2249c31..0000000 Binary files a/graphs/tarjan-biconn.mat and /dev/null differ diff --git a/index.md b/index.md new file mode 100644 index 0000000..6319fe5 --- /dev/null +++ b/index.md @@ -0,0 +1,11 @@ +--- +title: gaimc +--- + +[david]: http://www.stanford.edu/~dgleich + +gaimc +===== + +Written by [David F. Gleich][david] +Mathworks File Exchange Page *coming soon* diff --git a/largest_component.m b/largest_component.m deleted file mode 100644 index 2e70d02..0000000 --- a/largest_component.m +++ /dev/null @@ -1,46 +0,0 @@ -function [Acc,p] = largest_component(A,sym) -% LARGEST_COMPONENT Return the largest connected component of A -% -% Acc = largest_component(A) returns the largest connected component -% of the graph A. If A is directed, this returns the largest -% strongly connected component. -% -% Acc = largest_component(A,1) returns the largest connected piece of -% a directed graph where connectivity is undirected. Algorithmically, -% this takes A, drops the directions, then components the largest component -% and returns just this piece of the original _directed_ network. So the -% output Acc is directed in this case. -% -% [Acc,p] = largest_component(A,...) also returns a logical vector -% indicating which vertices in A were chosen. -% -% See also SCOMPONENTS -% -% Example: -% load_gaimc_graph('dfs_example') -% [Acc p] = largest_component(A); % compute the largest component -% xy2 = xy(p,:); labels2 = labels(p); % get component metadata -% % draw original graph -% subplot(1,2,1); graph_draw(A,xy,'labels',labels); title('Original'); -% % draw component -% subplot(1,2,2); graph_draw(Acc,xy2,'labels',labels2); title('Component'); - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2009-04-29: Initial coding - -if ~exist('sym','var') || isempty(sym), sym=0; end - -if sym - As = A|A'; - [ci sizes] = scomponents(As); -else - [ci sizes] = scomponents(A); -end -[csize cind] = max(sizes); -p = ci==cind; -Acc = A(p,p); - - diff --git a/license.txt b/license.txt deleted file mode 100644 index 2842d19..0000000 --- a/license.txt +++ /dev/null @@ -1,23 +0,0 @@ -Copyright (c) 2006-2015, David F. Gleich -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: - -* Redistributions of source code must retain the above copyright notice, this - list of conditions and the following disclaimer. - -* Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/load_gaimc_graph.m b/load_gaimc_graph.m deleted file mode 100644 index a19f27e..0000000 --- a/load_gaimc_graph.m +++ /dev/null @@ -1,35 +0,0 @@ -function varargout=load_gaimc_graph(graphname) -% LOAD_GAIMC_GRAPH Loads a graph from the gaimc library -% -% load_gaimc_graph is a helper function to load a graph provided with the -% library regardless of the current working directory. -% -% If it's called without any output arguments, it functions just like a -% load command executed on the .mat file with the graph. If it's called -% with an output arguemnt, it functions just like a load command with -% output arguments. It's somewhat complicated to explain because this is -% just a convinence function to make the examples work for any path, and -% not just from the gaimc root directory. -% -% Example: -% % equivalent to load('graphs/airports.mat') run from the gaimc directory -% load_gaimc_graph('airports') -% % equivalent to P=load('graphs/kt-7-2.mat') run from the gaimc directory -% P=load_gaimc_graph('kt-7-2.mat') -% % so you don't have to put the path in for examples! - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2009-04-27: Initial coding - - -path=fileparts(mfilename('fullpath')); -if nargout==0 - evalin('caller',['load(''' fullfile(path,'graphs',graphname) ''');']); -else - P = load(fullfile(path,'graphs',graphname)); - varargout{1} = P; -end - diff --git a/mst_prim.m b/mst_prim.m deleted file mode 100644 index 1532882..0000000 --- a/mst_prim.m +++ /dev/null @@ -1,142 +0,0 @@ -function varargout=mst_prim(A,full,u) -% MST_PRIM Compute a minimum spanning tree with Prim's algorithm -% -% T = mst_prim(A) computes a minimum spanning tree T using Prim's algorithm -% for the spanning tree of a graph with non-negative edge weights. -% -% T = mst_prim(A,0) produces an MST for just the component at A containing -% vertex 1. T = mst_prim(A,0,u) produces the MST for the component -% containing vertex u. -% -% [ti tj tv] = mst_prim(...) returns the edges from the matrix and does not -% convert to a sparse matrix structure. This saves a bit of work and is -% required when there are 0 edge weights. -% -% Example: -% load_gaimc_graph('airports'); % A(i,j) = negative travel time -% A = -A; % convert to travel time. -% A = max(A,A'); % make the travel times symmetric -% T = mst_prim(A); -% gplot(T,xy); % look at the minimum travel time tree in the US - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History: -% 2009-05-02: Added example -% 2009-11-25: Fixed bug with target - -% TODO: Add example - -if ~exist('full','var') || isempty(full), full=0; end -if ~exist('u','var') || isempty(u), u=1; end - -if isstruct(A), - rp=A.rp; ci=A.ci; ai=A.ai; - check=0; -else - [rp ci ai]=sparse_to_csr(A); - check=1; -end -if check && any(ai)<0, error('gaimc:prim', ... - 'prim''s algorithm cannot handle negative edge weights.'); -end -if check && ~isequal(A,A'), error('gaimc:prim', ... - 'prim''s algorithm requires an undirected graph.'); -end -nverts=length(rp)-1; -d=Inf*ones(nverts,1); T=zeros(nverts,1); L=zeros(nverts,1); -pred=zeros(1,length(rp)-1); - -% enter the main dijkstra loop -for iter=1:nverts - if iter==1, root=u; - else - root=mod(u+iter-1,nverts)+1; - if L(v)>0, continue; end - end - n=1; T(n)=root; L(root)=n; % oops, n is now the size of the heap - d(root) = 0; - while n>0 - v=T(1); L(v)=-1; ntop=T(n); T(1)=ntop; n=n-1; - if n>0, L(ntop)=1; end % pop the head off the heap - k=1; kt=ntop; % move element T(1) down the heap - while 1, - i=2*k; - if i>n, break; end % end of heap - if i==n, it=T(i); % only one child, so skip - else % pick the smallest child - lc=T(i); rc=T(i+1); it=lc; - if d(rc)ew - d(w)=ew; pred(w)=v; - % check if w is in the heap - k=L(w); onlyup=0; - if k==0 - % element not in heap, only move the element up the heap - n=n+1; T(n)=w; L(w)=n; k=n; kt=w; onlyup=1; - else kt=T(k); - end - % update the heap, move the element down in the heap - while 1 && ~onlyup, - i=2*k; - if i>n, break; end % end of heap - if i==n, it=T(i); % only one child, so skip - else % pick the smallest child - lc=T(i); rc=T(i+1); it=lc; - if d(rc)1, % j==1 => element at top of heap - j2=floor(j/2); tj2=T(j2); % parent element - if d(tj2)0, nmstedges=nmstedges+1; end -end -ti = zeros(nmstedges,1); tj=ti; tv = zeros(nmstedges,1); -k=1; -for i=1:nverts - if pred(i)>0, - j = pred(i); - ti(k)=i; tj(k)=j; - for rpi=rp(i):rp(i+1)-1 - if ci(rpi)==j, tv(k)=ai(rpi); break; end - end - k=k+1; - end -end -if nargout==1, - T = sparse(ti,tj,tv,nverts,nverts); - T = T + T'; - varargout{1} = T; -else - varargout = {ti, tj, tv}; -end - diff --git a/scomponents.m b/scomponents.m deleted file mode 100644 index 65a7cea..0000000 --- a/scomponents.m +++ /dev/null @@ -1,73 +0,0 @@ -function [sci sizes] = scomponents(A) -% SCOMPONENTS Compute the strongly connected components of a graph -% -% ci=scomponents(A) returns an index for the component number of every -% vertex in the graph A. The total number of components is max(ci). -% If the input is undirected, then this algorithm outputs just the -% connected components. Otherwise, it output the strongly connected -% components. -% -% The implementation is from Tarjan's 1972 paper: Depth-first search and -% linear graph algorithms. In SIAM's Journal of Computing, 1972, 1, -% pp.146-160. -% -% See also DMPERM -% -% Example: -% load_gaimc_graph('cores_example'); % the graph A has three components -% ci = scomponents(A) -% ncomp = max(ci) % should be 3 -% R = sparse(1:size(A,1),ci,1,size(A,1),ncomp); % create a restriction matrix -% CG = R'*A*R; % create the graph with each component -% % collapsed into a single node. - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2008-04-11: Initial coding - -if isstruct(A), rp=A.rp; ci=A.ci; %ofs=A.offset; -else [rp ci]=sparse_to_csr(A); -end - -n=length(rp)-1; sci=zeros(n,1); cn=1; -root=zeros(n,1); dt=zeros(n,1); t=0; -cs=zeros(n,1); css=0; % component stack -rs=zeros(2*n,1); rss=0; % recursion stack holds two nums (v,ri) -% start dfs at 1 -for sv=1:n - v=sv; if root(v)>0, continue; end - rss=rss+1; rs(2*rss-1)=v; rs(2*rss)=rp(v); % add v to the stack - root(v)=v; sci(v)=-1; dt(v)=t; t=t+1; - css=css+1; cs(css)=v; % add w to component stack - while rss>0 - v=rs(2*rss-1); ri=rs(2*rss); rss=rss-1; % pop v from the stack - while ridt(root(w)), root(v)=root(w); end - end - end - if root(v)==v - while css>0 - w=cs(css); css=css-1; sci(w)=cn; - if w==v, break; end - end - cn=cn+1; - end - end -end - -if nargout>1 - sizes=accumarray(sci,1,[max(sci) 1]); -end diff --git a/sparse_to_csr.m b/sparse_to_csr.m deleted file mode 100644 index d36abf0..0000000 --- a/sparse_to_csr.m +++ /dev/null @@ -1,101 +0,0 @@ -function [rp ci ai ncol]=sparse_to_csr(A,varargin) -% SPARSE_TO_CSR Convert a sparse matrix into compressed row storage arrays -% -% [rp ci ai] = sparse_to_csr(A) returns the row pointer (rp), column index -% (ci) and value index (ai) arrays of a compressed sparse representation of -% the matrix A. -% -% [rp ci ai] = sparse_to_csr(i,j,v,n) returns a csr representation of the -% index sets i,j,v with n rows. -% -% Example: -% A=sparse(6,6); A(1,1)=5; A(1,5)=2; A(2,3)=-1; A(4,1)=1; A(5,6)=1; -% [rp ci ai]=sparse_to_csr(A) -% -% See also CSR_TO_SPARSE, SPARSE - -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% History -% 2008-04-07: Initial version -% 2008-04-24: Added triple array input -% 2009-05-01: Added ncol output -% 2009-05-15: Fixed triplet input - -error(nargchk(1, 5, nargin, 'struct')) -retc = nargout>1; reta = nargout>2; - -if nargin>1 - if nargin>4, ncol = varargin{4}; end - nzi = A; nzj = varargin{1}; - if reta && length(varargin) > 2, nzv = varargin{2}; end - if nargin<4, n=max(nzi); else n=varargin{3}; end - nz = length(A); - if length(nzi) ~= length(nzj), error('gaimc:invalidInput',... - 'length of nzi (%i) not equal to length of nzj (%i)', nz, ... - length(nzj)); - end - if reta && length(varargin) < 3, error('gaimc:invalidInput',... - 'no value array passed for triplet input, see usage'); - end - if ~isscalar(n), error('gaimc:invalidInput',... - ['the 4th input to sparse_to_csr with triple input was not ' ... - 'a scalar']); - end - if nargin < 5, ncol = max(nzj); - elseif ~isscalar(ncol), error('gaimc:invalidInput',... - ['the 5th input to sparse_to_csr with triple input was not ' ... - 'a scalar']); - end -else - n = size(A,1); nz = nnz(A); ncol = size(A,2); - retc = nargout>1; reta = nargout>2; - if reta, [nzi nzj nzv] = find(A); - else [nzi nzj] = find(A); - end -end -if retc, ci = zeros(nz,1); end -if reta, ai = zeros(nz,1); end -rp = zeros(n+1,1); -for i=1:nz - rp(nzi(i)+1)=rp(nzi(i)+1)+1; -end -rp=cumsum(rp); -if ~retc && ~reta, rp=rp+1; return; end -for i=1:nz - if reta, ai(rp(nzi(i))+1)=nzv(i); end - ci(rp(nzi(i))+1)=nzj(i); - rp(nzi(i))=rp(nzi(i))+1; -end -for i=n:-1:1 - rp(i+1)=rp(i); -end -rp(1)=0; -rp=rp+1; - -% Copyright (c) 2006-2015, David F. Gleich -% All rights reserved. -% -% Redistribution and use in source and binary forms, with or without -% modification, are permitted provided that the following conditions are met: -% -% * Redistributions of source code must retain the above copyright notice, this -% list of conditions and the following disclaimer. -% -% * Redistributions in binary form must reproduce the above copyright notice, -% this list of conditions and the following disclaimer in the documentation -% and/or other materials provided with the distribution. -% -% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE -% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - diff --git a/test/dijkstra_perf.m b/test/dijkstra_perf.m deleted file mode 100644 index ea32211..0000000 --- a/test/dijkstra_perf.m +++ /dev/null @@ -1,33 +0,0 @@ -%% Dijkstra's algorithm -% To evaluate the performance of Dijkstra's algorithm, we pick -graphdir = '../graphs/'; -graphs = {'cs-stanford', 'tapir'}; -profile off; -if exist('prof','var') && prof, profile on; end -nrep=15; ntests=10; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; -for rep=1:nrep - for gi=1:length(graphs) - load([graphdir graphs{gi} '.mat']); n=size(A,1); - At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - for ti=1:ntests - fprintf([repmat('\b',1,66) '%20s rep=%3i graph=%-20s trial=%4i'], ... - 'dijkstra',rep,graphs{gi},ti); - v=ceil(n*rand(1)); - tic; d2=dijkstra_sp(At,v,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; d4=dijkstra2(As,v); mat_fast=mat_fast+toc; - if any(d2 ~= d4) - error('gaimc:dijkstra','incorrect results from dijkstra'); - end - end - end -end -fprintf('\n'); -fprintf('mex time: %f\n', mex_fast); -fprintf('mat time: %f\n', mat_fast); -fprintf('\n'); -if exist('prof','var') && prof, profile report; end -profile off; - - - diff --git a/test/prim_mst_perf.m b/test/prim_mst_perf.m deleted file mode 100644 index 7fd882e..0000000 --- a/test/prim_mst_perf.m +++ /dev/null @@ -1,38 +0,0 @@ -%% Help debug prim_mst performance -% -profile off; -if exist('prof','var') && prof, profile on; end -nrep=15; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; -comp_results=[]; -szs=[1 10000]; -for szi=1:length(szs) - fprintf('\n%20s size=%-5i ', 'mst_prim', szs(szi)); - % Matlab needs 1 iteration to compile the function - if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end - for rep=1:nrep - fprintf('\b\b\b\b'); fprintf(' %3i', rep); - A=abs(sprandsym(szs(szi),25/szs(szi))); - At=A'; - [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai; - tic; [t1i t1j t1v]=prim_mst(At,struct('istrans',1,'nocheck',1)); - mex_fast=mex_fast+toc; - tic; [t2i t2j t2v]=mst_prim2(As,0); mat_fast=mat_fast+toc; - T1f=sparse(t1i,t1j,t1v,size(A,1),size(A,2)); - T2f=sparse(t2i,t2j,t2v,size(A,1),size(A,2)); - if ~isequal(T1f+T1f',T2f+T2f') - warning('gaimc:mst_prim',... - 'incorrect results from mst_prim (%i,%i)', szi, rep); - fprintf('sum diff: %g\n', full(sum(sum(T1f))-sum(sum(T2f)))); - fprintf('cmp diff: %10g %10g\n', ... - max(scomponents(T1f+T1f')), max(scomponents(T2f+T2f'))); - keyboard - end - end - comp_results(end+1,:) = [mex_fast mat_fast]; -end -fprintf('\n'); -fprintf('mex time: %f\n', mex_fast); -fprintf('mat time: %f\n', mat_fast); -fprintf('\n'); -if exist('prof','var') && prof, profile report; end -profile off; \ No newline at end of file diff --git a/test/test_bfs.m b/test/test_bfs.m deleted file mode 100644 index c08e76d..0000000 --- a/test/test_bfs.m +++ /dev/null @@ -1,2 +0,0 @@ -function test_bfs -end diff --git a/test/test_bipartite_matching.m b/test/test_bipartite_matching.m deleted file mode 100644 index e029f56..0000000 --- a/test/test_bipartite_matching.m +++ /dev/null @@ -1,45 +0,0 @@ -function test_bipartite_matching -%% Info -% David F. Gleich -% Copyright, Stanford University, 2008-2009 - -% 2009-05-15: Initial version - -%% empty input -val = bipartite_matching([]); - -%% exact answer -A = [5 1 1 1 1; - 1 5 1 1 1; - 1 1 5 1 1; - 1 1 1 5 1; - 1 1 1 1 5; - 1 1 1 1 1]; -[val m1 m2] = bipartite_matching(A); -if val ~= 25 || any(m1 ~= m2) - error('gaimc:test_bipartite_matching','failed on known answer'); -end - -%% test triple input/output -A = ones(6,5); A = A-spdiags(diag(A),0,size(A,1),size(A,2)); -[ai aj av] = find(A); -ai = [(1:5)';ai]; -aj = [(1:5)';aj]; -av = [5*ones(5,1);av]; -[val m1 m2 mi] = bipartite_matching(av,ai,aj); -mitrue = zeros(length(av),1); mitrue(1:5)=1; -if val ~= 25 || any(m1 ~= m2) || sum(mi) ~= 5 || any(mi~=mitrue) - error('gaimc:test_bipartite_matching','failed on known answer'); -end - -%% 100 random trials against dmperm -for t=1:100 - A = sprand(500,400,5/500); % 5 nz per row - A = spones(A); - p = dmperm(A); - val = bipartite_matching(A); - if sum(p>0)~=val - error('gaimc:test_bipartite_matching','failed unweighted dmperm test'); - end -end - diff --git a/test/test_corenums.m b/test/test_corenums.m deleted file mode 100644 index 5d01412..0000000 --- a/test/test_corenums.m +++ /dev/null @@ -1,3 +0,0 @@ -function test_corenums - -end \ No newline at end of file diff --git a/test/test_csr_to_sparse.m b/test/test_csr_to_sparse.m deleted file mode 100644 index f0c1152..0000000 --- a/test/test_csr_to_sparse.m +++ /dev/null @@ -1,33 +0,0 @@ -function test_csr_to_sparse -%% empty arguments -A = csr_to_sparse(1,[],[]); -if ~isequal(A,sparse([])) - error('gaimc:csr_to_sparse','failed to convert empty matrix'); -end - -A = csr_to_sparse(1,[],[],50); -if size(A,2)~=50 - error('gaimc:csr_to_sparse','incorrect empty size output'); -end - -%% exact test -% clique to sparse -rp = 1:5:26; -ci = reshape(repmat(1:5,5,1)',25,1); -ai = ones(25,1); -A=csr_to_sparse(rp,ci,ai); -if ~isequal(A,sparse(full(ones(5)))) - error('gaimc:csr_to_sparse','failed to convert clique'); -end - -%% 100 random trials of round trips between sparse_to_csr and csr_to_sparse -for t=1:100 - A = sprand(100,80,0.01); - [rp ci ai]=sparse_to_csr(A); - A2 = csr_to_sparse(rp,ci,ai,80); - if ~isequal(A,A2) - error('gaimc:csr_to_sparse','random sparse test failed'); - end -end - - diff --git a/test/test_dfs.m b/test/test_dfs.m deleted file mode 100644 index 2b9eaea..0000000 --- a/test/test_dfs.m +++ /dev/null @@ -1,15 +0,0 @@ -function test_dfs - -%% Line graph test -for n=10:10:100 - A=spdiags(ones(n,2),[-1 1],n,n); - [d dt ft pred]=dfs(A,1); - if any(d~=(0:n-1)') || ... - any(dt~=(0:n-1)') || ... - any(ft~=(2*n-1:-1:n)') || ... - any(pred~=0:n-1) - error('gaimc:dfs','line graph test error'); - end -end - -%% \ No newline at end of file diff --git a/test/test_examples.m b/test/test_examples.m deleted file mode 100644 index 5e351b3..0000000 --- a/test/test_examples.m +++ /dev/null @@ -1,88 +0,0 @@ -%% bfs -load_gaimc_graph('bfs_example.mat') % use the dfs example from Boost -d = bfs(A,1) - -%% bipartite_matching -A = rand(10,8); % bipartite matching between random data -[val mi mj] = bipartite_matching(A); -val - -%% clustercoeffs -load_gaimc_graph('clique-10'); -cc = clustercoeffs(A) % they are all equal! as we expect in a clique - - -%% dfs -load_gaimc_graph('dfs_example.mat') % use the dfs example from Boost -d = dfs(A,1) - -%% dijkstra -% Find the minimum travel time between Los Angeles (LAX) and -% Rochester Minnesota (RST). -load_gaimc_graph('airports') -A = -A; % fix funny encoding of airport data -lax=247; rst=355; -[d pred] = dijkstra(A,lax); -fprintf('Minimum time: %g\n',d(rst)); -% Print the path -fprintf('Path:\n'); -path =[]; u = rst; while (u ~= lax) path=[u path]; u=pred(u); end -fprintf('%s',labels{lax}); -for i=path; fprintf(' --> %s', labels{i}); end, fprintf('\n'); - -%% dirclustercoeffs -load_gaimc_graph('celegans'); % load the C elegans nervous system network -cc=dirclustercoeffs(A); -[maxval maxind]=max(cc) -labels(maxind) % most clustered vertex in the nervous system - -%% graph_draw -load_gaimc_graph('dfs_example'); -graph_draw(A,xy); - - -%% mst_prim -load_gaimc_graph('airports'); % A(i,j) = negative travel time -A = -A; % convert to travel time. -A = max(A,A'); % make the travel times symmetric -T = mst_prim(A); -gplot(T,xy); % look at the minimum travel time tree in the US - -%% scomponents -% scomponents -load_gaimc_graph('cores_example'); % the graph A has three components -ci = scomponents(A) -ncomp = max(ci) % should be 3 -R = sparse(1:size(A,1),ci,1,size(A,1),ncomp); % create a restriction matrix -CG = R'*A*R; % create the graph with each component - % collapsed into a single node. - -%% load_gaimc_graph -% equivalent to load('graphs/airports.mat') run from the gaimc directory -load_gaimc_graph('airports') -% equivalent to P=load('graphs/kt-7-2.mat') run from the gaimc directory -P=load_gaimc_graph('kt-7-2.mat') -% so you don't have to put the path in for examples! - - -%% largest_component -load_gaimc_graph('dfs_example') -[Acc p] = largest_component(A); % compute the largest component -xy2 = xy(p,:); labels2 = labels(p); % get component metadata -% draw original graph -subplot(1,2,1); graph_draw(A,xy,'labels',labels); title('Original'); -% draw component -subplot(1,2,2); graph_draw(Acc,xy2,'labels',labels2); title('Component'); - -%% corenums -load_gaimc_graph('cores_example'); % the graph A has three components -corenums(A) - -%% sparse_to_csr -A=sparse(6,6); A(1,1)=5; A(1,5)=2; A(2,3)=-1; A(4,1)=1; A(5,6)=1; -[rp ci ai]=sparse_to_csr(A); - -%% csr_to_sparse -A=sparse(6,6); A(1,1)=5; A(1,5)=2; A(2,3)=-1; A(4,1)=1; A(5,6)=1; -[rp ci ai]=sparse_to_csr(A); -A2 = csr_to_sparse(rp,ci,ai); diff --git a/test/test_floydwarshall.m b/test/test_floydwarshall.m deleted file mode 100644 index 877d9e0..0000000 --- a/test/test_floydwarshall.m +++ /dev/null @@ -1,30 +0,0 @@ -function test_floydwarshall - -%% Compare against dijkstra -load_gaimc_graph('all_shortest_paths_example'); -A = spfun(@(x) x-min(min(A))+1,A); % remove the negative edges -n = size(A,1); -% Now, we'll run Dijkstra's algorithm for every vertex and save the result - -D2 = zeros(n,n); -for i=1:n - D2(i,:) = dijkstra(A,i); -end -D = floydwarshall(A); -if any(D-D2) - error('test:floydwarshall','floyd warshall returned incorrect distances'); -end - -%% Compare predecessors -load_gaimc_graph('all_shortest_paths_example'); -A = spfun(@(x) x-min(min(A))+1,A); % remove the negative edges -n = size(A,1); -P2 = zeros(n,n); -for i=1:n - [d,p]=dijkstra(A,i); - P2(i,:) = p; -end -[D,P] = floydwarshall(A); -if any(P-P2) - error('test:floydwarshall','floyd warshall returned incorrect predecessors'); -end diff --git a/test/test_largest_component.m b/test/test_largest_component.m deleted file mode 100644 index f55f0bc..0000000 --- a/test/test_largest_component.m +++ /dev/null @@ -1,33 +0,0 @@ -function test_largest_component - -fname = 'largest_component'; % function name -ntest = 0; % test number - -ntest = ntest+1; tid = sprintf('%32s test %3i : ', fname, ntest); -load_gaimc_graph('dfs_example'); -Acc = largest_component(A); -if size(Acc,1) ~= 5 - error('gaimc:test','%s failed on dfs_example', tid); -else - fprintf([tid 'passed\n']); -end - -ntest = ntest+1; tid = sprintf('%32s test %3i : ', fname, ntest); -load_gaimc_graph('dfs_example'); -Acc = largest_component(A,1); -if size(Acc,1) ~= 6 - error('gaimc:test','%s failed on sym=1 dfs_example', tid); -else - fprintf([tid 'passed\n']); -end - -ntest = ntest+1; tid = sprintf('%32s test %3i : ', fname, ntest); -load_gaimc_graph('cores_example'); % the graph A is symmetric -Acc1 = largest_component(A); -Acc2 = largest_component(A,1); -if ~isequal(Acc1,Acc2) - error('gaimc:test','%s failed on cores_example', tid); -else - fprintf([tid 'passed\n']); -end - \ No newline at end of file diff --git a/test/test_load_gaimc_graph.m b/test/test_load_gaimc_graph.m deleted file mode 100644 index b30d7e0..0000000 --- a/test/test_load_gaimc_graph.m +++ /dev/null @@ -1,31 +0,0 @@ -function load_test_gaimc_graph - -%% -P1=load('../graphs/kt-7-2.mat'); -P2=load_gaimc_graph('kt-7-2.mat'); -if ~isequal(P1,P2) - error('gaimc_graph failed on kt-7-2.mat'); -else - fprintf('gaimc_graph passed load test on kt-7-2.mat\n'); -end - -%% -P1=load('../graphs/kt-7-2'); -P2=load_gaimc_graph('kt-7-2'); -if ~isequal(P1,P2) - error('gaimc_graph failed on kt-7-2'); -else - fprintf('gaimc_graph passed load test on kt-7-2\n'); -end - -%% -load('../graphs/clr-24-1'); -P1 = struct('A',A,'labels',labels,'xy',xy); -load_gaimc_graph('clr-24-1'); -P2 = struct('A',A,'labels',labels,'xy',xy); - -if ~isequal(P1,P2) - error('gaimc_graph failed on clr-24-1'); -else - fprintf('gaimc_graph passed load test on clr-24-1\n'); -end \ No newline at end of file diff --git a/test/test_main.m b/test/test_main.m deleted file mode 100644 index dd12019..0000000 --- a/test/test_main.m +++ /dev/null @@ -1,15 +0,0 @@ -function test_main -% TODO Check the directory -test_examples -test_sparse_to_csr -test_csr_to_sparse % should always come after sparse_to_csr -test_bipartite_matching -test_dfs -test_load_gaimc_graph -test_largest_component -test_examples -% test_dijkstra % required for floydwarshall -test_floydwarshall -test_mst_prim - -end diff --git a/test/test_mst_prim.m b/test/test_mst_prim.m deleted file mode 100644 index 1c38ded..0000000 --- a/test/test_mst_prim.m +++ /dev/null @@ -1,21 +0,0 @@ -function test_mst_prim -msgid = 'gaimc:mst_prim'; -load_gaimc_graph('clr-24-1'); -A(2,3) = 9; A(3,2) = 9; -T = mst_prim(A); % root the tree at vertex -Ttrueijv = [ - 2 8 1 4 6 9 3 5 4 3 7 6 8 1 7 3 - 1 1 2 3 3 3 4 4 5 6 6 7 7 8 8 9 - 4 8 4 7 4 2 7 9 9 4 2 2 1 8 1 2 ]; -Ttrue = sparse(Ttrueijv(1,:), Ttrueijv(2,:), Ttrueijv(3,:), 9,9); -if nnz(T - Ttrue) ~= 0 - error(msgid, 'mst_prim failed'); -end - -load_gaimc_graph('clr-24-1'); -T1 = mst_prim(A); % root the tree at vertex -T2 = mst_prim(A,[],5); -T1T2diff = sparse(9,9); T1T2diff(2,3) = 8; T1T2diff(1,8) = -8; -if nnz(triu(T1-T2)-T1T2diff) ~= 0 - error(msgid, 'mst_prim failed rooted test'); -end diff --git a/test/test_sparse_to_csr.m b/test/test_sparse_to_csr.m deleted file mode 100644 index ff4f65c..0000000 --- a/test/test_sparse_to_csr.m +++ /dev/null @@ -1,24 +0,0 @@ -function test_sparse_to_csr -%% Previous failure -[ai,aj,av]=find(ones(5)); -sparse_to_csr(ai,aj,av); - -%% 100 random trials -for t=1:100 - A = sprand(100,80,0.01); - [rp ci ai]=sparse_to_csr(A); - i=zeros(length(ai),1); j=i; a=i; - n = length(rp)-1; nz=0; - for cr=1:n - for ri=rp(cr):rp(cr+1)-1 - nz=nz+1; i(nz)=cr; j(nz)=ci(ri); a(nz)=ai(ri); - end - end - A2 = sparse(i,j,a,n,80); - if ~isequal(A,A2) - error('gaimc:sparse_to_csr','random sparse test failed'); - end -end - -%% empty arguments -[rp ci ai] = sparse_to_csr([]);