diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d0ab673 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*~ +*.asv + diff --git a/Contents.m b/Contents.m index 145ade2..faf8c88 100644 --- a/Contents.m +++ b/Contents.m @@ -2,7 +2,7 @@ % Graph Algorithms in Matlab Code (gaimc) % Written by David Gleich % Version 1.0 (beta) -% 2008-04-21 +% 2008-2009 %========================================= % % Search algorithms @@ -17,17 +17,34 @@ % % 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 +% load_gaimc_graph - Loads a sample graph from the library % David Gleich -% Copyright, Stanford University, 2008 +% Copyright, Stanford University, 2008-2009 % History % 2008-04-10: Initial version + + +% TODO for release +% Fix mlintrpt errors +% Update copyright info everywhere +% Implement bipartite matching + + +% Future todos +% Implement weighted core nums +% More testing +% Implement all pairs shortest paths with Floyd Warshall \ No newline at end of file diff --git a/bfs.m b/bfs.m index f069647..e6033d7 100644 --- a/bfs.m +++ b/bfs.m @@ -11,6 +11,8 @@ % [...] = 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 diff --git a/clustercoeffs.m b/clustercoeffs.m index 115868a..6b33c96 100644 --- a/clustercoeffs.m +++ b/clustercoeffs.m @@ -8,11 +8,14 @@ % 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. +% 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 if ~exist('normalized','var') || isempty(normalized), normalized=true; end if ~exist('weighted','var') || isempty(weighted), weighted=true; end @@ -37,7 +40,7 @@ end n=length(rp)-1; -cc=zeros(n,1); ind=false(n,1); cache=zeros(n,1); degs=zeros(n,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 diff --git a/convert_sparse.m b/convert_sparse.m new file mode 100644 index 0000000..4d18d11 --- /dev/null +++ b/convert_sparse.m @@ -0,0 +1,22 @@ +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 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; \ No newline at end of file diff --git a/corenums.m b/corenums.m index 5869154..83e0582 100644 --- a/corenums.m +++ b/corenums.m @@ -15,7 +15,7 @@ % Decomposition of Networks." Sept. 1 2002. % % Example: -% load('graphs/cores_example.mat') +% load_gaimc_graph('cores_example'); % the graph A has three components % corenums(A) % diff --git a/demo/airports.m b/demo/airports.m new file mode 100644 index 0000000..0f2aad4 --- /dev/null +++ b/demo/airports.m @@ -0,0 +1,124 @@ +%% 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 +%% +% We need to clear the axes after the mapping toolbox +clf; +%% +% Let's just look at the continential US too. +figure; 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. +% +% Clear the figure again before proceeding. +clf; + +%% 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 new file mode 100644 index 0000000..68fe5da --- /dev/null +++ b/demo/demo.m @@ -0,0 +1,308 @@ +%% 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('graphs/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('graphs/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); + +%% +% 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. +[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('graphs/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 ('graphs/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) \ No newline at end of file diff --git a/dfs.m b/dfs.m index dc5e258..5b61b78 100644 --- a/dfs.m +++ b/dfs.m @@ -19,6 +19,8 @@ % yet. % % Example: +% load_gaimc_graph('dfs_example.mat') % use the dfs example from Boost +% d = dfs(A,1) % % See also BFS diff --git a/dijkstra.m b/dijkstra.m index 52a050f..2a840b2 100644 --- a/dijkstra.m +++ b/dijkstra.m @@ -1,10 +1,42 @@ function [d pred]=dijkstra(A,u) +% DIJKSTRA Compute shortest paths using Dijkstra's algorithm +% +% d=dijkstra(A,u) computes the shortest path from vertex u to all nodes +% reachable from vertex u using Dijkstra's algorithm for the problem. +% The graph is given by the weighted sparse matrix A, where A(i,j) is +% the distance between vertex i and j. In the output vector d, +% the entry d(v) is the minimum distance between vertex u and vertex v. +% A vertex w unreachable from u has d(w)=Inf. +% +% [d pred]=dijkstra(A,u) also returns the predecessor tree to generate +% the actual shorest paths. In the predecessor tree pred(v) is the +% vertex preceeding v in the shortest path and pred(u)=0. Any +% unreachable vertex has pred(w)=0 as well. +% +% If your network is unweighted, then use bfs instead. +% +% See also BFS +% +% Example: +% % 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'); % David Gleich -% Copyright, Stanford University, 2008 +% Copyright, Stanford University, 2008-2009 % History % 2008-04-09: Initial coding +% 2009-05-15: Documentation if isstruct(A), rp=A.rp; ci=A.ci; ai=A.ai; diff --git a/dirclustercoeffs.m b/dirclustercoeffs.m index 32f85ee..911fea5 100644 --- a/dirclustercoeffs.m +++ b/dirclustercoeffs.m @@ -2,11 +2,30 @@ % DIRCLUSTERCOEFFS Compute clustering coefficients for a directed graph % % cc=dirclustercoeffs(A) returns the directed clustering coefficients -% (which are identical to the clustering coefficients on an undirected -% graph. +% (which generalize the clustering coefficients of an undirected graph, +% and so calling this function on an undirected graph will produce the same +% answer as clustercoeffs, but less efficiently.) +% +% This function implements the algorithm from Fagiolo, Phys Rev. E. 76 +% 026107 (doi:10:1103/PhysRevE.76.026107). +% +% [cc,cccyc,ccmid,ccin,ccout,nf]=dirclusteringcoeffs(A) returns different +% components of the clustering coefficients corresponding to cycles, +% middles, in triangles, and out triangles. See the manuscript for a +% description of the various types of triangles counted in the above +% metrics. +% +% See also CLUSTERCOEFFS +% +% Example: +% 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 % History % 2008-04-22: Initial coding +% 2009-05-15: Documentation and example if ~exist('normalized','var') || isempty(normalized), normalized=true; end if ~exist('weighted','var') || isempty(weighted), weighted=true; end @@ -37,9 +56,12 @@ 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); +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 + end +end ew=1; ew2=1; for v=1:n % setup counts for the different cycle types diff --git a/graph_draw.m b/graph_draw.m new file mode 100644 index 0000000..dc5b69e --- /dev/null +++ b/graph_draw.m @@ -0,0 +1,588 @@ +function h = graph_draw(adj, xy, varargin) +% GRAPH_DRAW Draw a picture of a graph when the coordinates are known +% +% graph_draw(A, xy) draws a picture of graph A where node i is placed +% at x = xy(i,1), y = xy(i,2). In the drawing, shaded nodes have +% self loops. +% +% Some of the parameters of the drawing are controlled by specifying +% optional parameters in the call graph_draw(A, xy, key, value). The keys +% and default values are +% 'linestyle' - default '-' +% 'linewidth' - default .5 +% 'linecolor' - default Black +% 'fontsize' - fontsize for labels, default 8 +% 'labels' - Cell array containing labels +% '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 new file mode 100644 index 0000000..6134dd0 Binary files /dev/null and b/graphs/airports.mat differ diff --git a/graphs/bfs_example.mat b/graphs/bfs_example.mat index c3ef97c..7ee9b0b 100644 Binary files a/graphs/bfs_example.mat and b/graphs/bfs_example.mat differ diff --git a/graphs/celegans.mat b/graphs/celegans.mat new file mode 100644 index 0000000..79d1b29 Binary files /dev/null and b/graphs/celegans.mat differ diff --git a/graphs/dfs_example.mat b/graphs/dfs_example.mat index 7d85743..2698814 100644 Binary files a/graphs/dfs_example.mat and b/graphs/dfs_example.mat differ diff --git a/largest_component.m b/largest_component.m new file mode 100644 index 0000000..21878dd --- /dev/null +++ b/largest_component.m @@ -0,0 +1,46 @@ +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 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/load_gaimc_graph.m b/load_gaimc_graph.m new file mode 100644 index 0000000..82430a0 --- /dev/null +++ b/load_gaimc_graph.m @@ -0,0 +1,35 @@ +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 Gleich +% Copyright, Stanford University, 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 + \ No newline at end of file diff --git a/mst_prim.m b/mst_prim.m index b24a940..96219a0 100644 --- a/mst_prim.m +++ b/mst_prim.m @@ -11,6 +11,21 @@ % [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 Gleich +% Copyright, Stanford University, 2008-2009 + +% History: +% 2009-05-02: Added example + +% TODO: Add example if ~exist('full','var') || isempty(full), full=0; end if ~exist('target','var') || isempty(full), u=1; end @@ -23,9 +38,11 @@ check=1; end if check && any(ai)<0, error('gaimc:prim', ... - 'prim''s algorithm cannot handle negative edge weights.'); end + '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 + '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); @@ -33,7 +50,10 @@ % 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 + 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 @@ -118,4 +138,4 @@ else varargout = {ti, tj, tv}; end - \ No newline at end of file + diff --git a/scomponents.m b/scomponents.m index 92be09b..3cf27ac 100644 --- a/scomponents.m +++ b/scomponents.m @@ -1,14 +1,25 @@ function [sci sizes] = scomponents(A) % SCOMPONENTS Compute the strongly connected components of a graph % -% ci=scomponents(A) returns the +% 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 implement is from Tarjan's 1972 paper: Depth-first search and linear -% graph algorithms. In SIAM's Journal of Computing, 1972, 1, pp.146-160. -% -% Example: +% 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('graphs/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 Gleich % Copyright, Stanford University, 2008 diff --git a/sparse_to_csr.m b/sparse_to_csr.m index 5a746ce..e617ece 100644 --- a/sparse_to_csr.m +++ b/sparse_to_csr.m @@ -12,7 +12,7 @@ % 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 +% See also CSR_TO_SPARSE, SPARSE % David Gleich % Copyright, Stanford University, 2008 @@ -33,9 +33,11 @@ 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 + length(nzj)); + end if reta && length(varargin) < 3, error('gaimc:invalidInput',... - 'no value array passed for triplet input, see usage'); end + '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 @@ -69,3 +71,4 @@ end rp(1)=0; rp=rp+1; + diff --git a/test/test_bfs.m b/test/test_bfs.m index e69de29..c08e76d 100644 --- a/test/test_bfs.m +++ b/test/test_bfs.m @@ -0,0 +1,2 @@ +function test_bfs +end diff --git a/test/test_examples.m b/test/test_examples.m new file mode 100644 index 0000000..1b268f8 --- /dev/null +++ b/test/test_examples.m @@ -0,0 +1,75 @@ +%% bfs +load_gaimc_graph('bfs_example.mat') % use the dfs example from Boost +d = bfs(A,1) + +%% 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('graphs/cores_example'); % the graph A has three components +ci = scomponents(A) +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) + diff --git a/test/test_largest_component.m b/test/test_largest_component.m new file mode 100644 index 0000000..f55f0bc --- /dev/null +++ b/test/test_largest_component.m @@ -0,0 +1,33 @@ +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 new file mode 100644 index 0000000..b30d7e0 --- /dev/null +++ b/test/test_load_gaimc_graph.m @@ -0,0 +1,31 @@ +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 index 729c2a2..f9c58b8 100644 --- a/test/test_main.m +++ b/test/test_main.m @@ -1,4 +1,8 @@ function test_main +% TODO Check the directory test_sparse_to_csr test_dfs +test_load_gaimc_graph +test_largest_component + end \ No newline at end of file