diff --git a/MATLAB_Functions/HuskyTrack.mltbx b/MATLAB_Functions/HuskyTrack.mltbx
new file mode 100644
index 0000000..f03d778
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack.mltbx differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/Contents.m b/MATLAB_Functions/HuskyTrack/gaimc/Contents.m
new file mode 100644
index 0000000..6bb21b7
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/Contents.m
@@ -0,0 +1,52 @@
+%=========================================
+% Graph Algorithms in Matlab Code (gaimc)
+% Written by David Gleich
+% Version 1.0 (beta)
+% 2008-2009
+%=========================================
+%
+% Search algorithms
+% dfs - depth first search
+% bfs - breadth first search
+%
+% Shortest path algorithms
+% dijkstra - Dijkstra's shortest path algorithm
+%
+% Minimum spanning tree algorithms
+% mst_prim - Compute an MST using Prim's algorithm
+%
+% Matching
+% bipartite_matching - Compute a maximum weight bipartite matching
+%
+% Connected components
+% scomponents - Compute strongly connected components
+% largest_component - Selects only the largest component
+%
+% Statistics
+% clustercoeffs - Compute clustering coefficients
+% dirclustercoeffs - Compute directed clustering coefficients
+% corenums - Compute core numbers
+%
+% Drawing
+% graph_draw - Draw an adjacency matrix (from Leon Peshkin)
+%
+% Helper functions
+% sparse_to_csr - Compressed sparse row arrays from a matrix
+% csr_to_sparse - Convert back to Matlab sparse matrices
+% load_gaimc_graph - Loads a sample graph from the library
+
+% David F. Gleich
+% Copyright, Stanford University, 2008-2009
+
+% History
+% 2008-04-10: Initial version
+
+
+% TODO for release
+% Fix mlintrpt errors
+
+
+% Future todos
+% Implement weighted core nums
+% More testing
+% Implement all pairs shortest paths with Floyd Warshall
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/bfs.m b/MATLAB_Functions/HuskyTrack/gaimc/bfs.m
new file mode 100644
index 0000000..93c7d59
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/bfs.m
@@ -0,0 +1,49 @@
+function [d dt pred] = bfs(A,u,target)
+% BFS Compute breadth first search distances, times, and tree for a graph
+%
+% [d dt pred] = bfs(A,u) returns the distance (d) and the discover time
+% (dt) for each vertex in the graph in a breadth first search
+% starting from vertex u.
+% d = dt(i) = -1 if vertex i is not reachable from u
+% pred is the predecessor array. pred(i) = 0 if vertex (i)
+% is in a component not reachable from u and i != u.
+%
+% [...] = bfs(A,u,v) stops the bfs when it hits the vertex v
+%
+% Example:
+% load_gaimc_graph('bfs_example.mat') % use the dfs example from Boost
+% d = bfs(A,1)
+%
+% See also DFS
+
+% David F. Gleich
+% Copyright, Stanford University, 2008-20098
+
+% History
+% 2008-04-13: Initial coding
+
+if ~exist('target','var') || isempty(full), target=0; end
+
+if isstruct(A), rp=A.rp; ci=A.ci;
+else [rp ci]=sparse_to_csr(A);
+end
+
+n=length(rp)-1;
+d=-1*ones(n,1); dt=-1*ones(n,1); pred=zeros(1,n);
+sq=zeros(n,1); sqt=0; sqh=0; % search queue and search queue tail/head
+
+% start bfs at u
+sqt=sqt+1; sq(sqt)=u;
+t=0;
+d(u)=0; dt(u)=t; t=t+1; pred(u)=u;
+while sqt-sqh>0
+ sqh=sqh+1; v=sq(sqh); % pop v off the head of the queue
+ for ri=rp(v):rp(v+1)-1
+ w=ci(ri);
+ if d(w)<0
+ sqt=sqt+1; sq(sqt)=w;
+ d(w)=d(v)+1; dt(w)=t; t=t+1; pred(w)=v;
+ if w==target, return; end
+ end
+ end
+end
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/bipartite_matching.m b/MATLAB_Functions/HuskyTrack/gaimc/bipartite_matching.m
new file mode 100644
index 0000000..c35e5d7
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/bipartite_matching.m
@@ -0,0 +1,240 @@
+function [val m1 m2 mi]=bipartite_matching(varargin)
+% BIPARTITE_MATCHING Solve a maximum weight bipartite matching problem
+%
+% [val m1 m2]=bipartite_matching(A) for a rectangular matrix A
+% [val m1 m2 mi]=bipartite_matching(x,ei,ej,n,m) for a matrix stored
+% in triplet format. This call also returns a matching indicator mi so
+% that val = x'*mi.
+%
+% The maximum weight bipartite matching problem tries to pick out elements
+% from A such that each row and column get only a single non-zero but the
+% sum of all the chosen elements is as large as possible.
+%
+% This function is slightly atypical for a graph library, because it will
+% be primarily used on rectangular inputs. However, these rectangular
+% inputs model bipartite graphs and we take advantage of that stucture in
+% this code. The underlying graph adjency matrix is
+% G = spaugment(A,0);
+% where A is the rectangular input to the bipartite_matching function.
+%
+% Matlab already has the dmperm function that computes a maximum
+% cardinality matching between the rows and the columns. This function
+% gives us the maximum weight matching instead. For unweighted graphs, the
+% two functions are equivalent.
+%
+% Note: If ei and ej contain duplicate edges, the results of this function
+% are incorrect.
+%
+% See also DMPERM
+%
+% Example:
+% A = rand(10,8); % bipartite matching between random data
+% [val mi mj] = bipartite_matching(A);
+% val
+
+% David F. Gleich and Ying Wang
+% Copyright, Stanford University, 2008-2009
+% Computational Approaches to Digital Stewardship
+
+% 2008-04-24: Initial coding (copy from Ying Wang matching_sparse_mex.cpp)
+% 2008-11-15: Added triplet input/output
+% 2009-04-30: Modified for gaimc library
+% 2009-05-15: Fixed error with empty inputs and triple added example.
+
+[rp ci ai tripi n m] = bipartite_matching_setup(varargin{:});
+
+if isempty(tripi)
+ error(nargoutchk(0,3,nargout,'struct'));
+else
+ error(nargoutchk(0,4,nargout,'struct'));
+end
+
+
+if ~isempty(tripi) && nargout>3
+ [val m1 m2 mi] = bipartite_matching_primal_dual(rp, ci, ai, tripi, n, m);
+else
+ [val m1 m2] = bipartite_matching_primal_dual(rp, ci, ai, tripi, n, m);
+end
+
+function [rp ci ai tripi n m]= bipartite_matching_setup(A,ei,ej,n,m)
+% convert the input
+
+if nargin == 1
+ if isstruct(A)
+ [nzi nzj nzv]=csr_to_sparse(A.rp,A.ci,A.ai);
+ else
+ [nzi nzj nzv]=find(A);
+ end
+ [n m]=size(A);
+ triplet = 0;
+elseif nargin >= 3 && nargin <= 5
+ nzi = ei;
+ nzj = ej;
+ nzv = A;
+ if ~exist('n','var') || isempty(n), n = max(nzi); end
+ if ~exist('m','var') || isempty(m), m = max(nzj); end
+ triplet = 1;
+else
+ error(nargchk(3,5,nargin,'struct'));
+end
+nedges = length(nzi);
+
+rp = ones(n+1,1); % csr matrix with extra edges
+ci = zeros(nedges+n,1);
+ai = zeros(nedges+n,1);
+if triplet, tripi = zeros(nedges+n,1); % triplet index
+else tripi = [];
+end
+
+%
+% 1. build csr representation with a set of extra edges from vertex i to
+% vertex m+i
+%
+rp(1)=0;
+for i=1:nedges
+ rp(nzi(i)+1)=rp(nzi(i)+1)+1;
+end
+rp=cumsum(rp);
+for i=1:nedges
+ if triplet, tripi(rp(nzi(i))+1)=i; end % triplet index
+ ai(rp(nzi(i))+1)=nzv(i);
+ ci(rp(nzi(i))+1)=nzj(i);
+ rp(nzi(i))=rp(nzi(i))+1;
+end
+for i=1:n % add the extra edges
+ if triplet, tripi(rp(i)+1)=-1; end % triplet index
+ ai(rp(i)+1)=0;
+ ci(rp(i)+1)=m+i;
+ rp(i)=rp(i)+1;
+end
+% restore the row pointer array
+for i=n:-1:1
+ rp(i+1)=rp(i);
+end
+rp(1)=0;
+rp=rp+1;
+
+%
+% 1a. check for duplicates in the data
+%
+colind = false(m+n,1);
+for i=1:n
+ for rpi=rp(i):rp(i+1)-1
+ if colind(ci(rpi)), error('bipartite_matching:duplicateEdge',...
+ 'duplicate edge detected (%i,%i)',i,ci(rpi));
+ end
+ colind(ci(rpi))=1;
+ end
+ for rpi=rp(i):rp(i+1)-1, colind(ci(rpi))=0; end % reset indicator
+end
+
+
+function [val m1 m2 mi]=bipartite_matching_primal_dual(...
+ rp, ci, ai, tripi, n, m)
+% BIPARTITE_MATCHING_PRIMAL_DUAL
+
+alpha=zeros(n,1); % variables used for the primal-dual algorithm
+beta=zeros(n+m,1);
+queue=zeros(n,1);
+t=zeros(n+m,1);
+match1=zeros(n,1);
+match2=zeros(n+m,1);
+tmod = zeros(n+m,1);
+ntmod=0;
+
+
+%
+% initialize the primal and dual variables
+%
+for i=1:n
+ for rpi=rp(i):rp(i+1)-1
+ if ai(rpi) > alpha(i), alpha(i)=ai(rpi); end
+ end
+end
+% dual variables (beta) are initialized to 0 already
+% match1 and match2 are both 0, which indicates no matches
+i=1;
+while i<=n
+ % repeat the problem for n stages
+
+ % clear t(j)
+ for j=1:ntmod, t(tmod(j))=0; end
+ ntmod=0;
+
+
+ % add i to the stack
+ head=1; tail=1;
+ queue(head)=i; % add i to the head of the queue
+ while head <= tail && match1(i)==0
+ k=queue(head);
+ for rpi=rp(k):rp(k+1)-1
+ j = ci(rpi);
+ if ai(rpi) < alpha(k)+beta(j) - 1e-8, continue; end % skip if tight
+ if t(j)==0,
+ tail=tail+1; queue(tail)=match2(j);
+ t(j)=k;
+ ntmod=ntmod+1; tmod(ntmod)=j;
+ if match2(j)<1,
+ while j>0,
+ match2(j)=t(j);
+ k=t(j);
+ temp=match1(k);
+ match1(k)=j;
+ j=temp;
+ end
+ break; % we found an alternating path
+ end
+ end
+ end
+ head=head+1;
+ end
+
+ if match1(i) < 1, % still not matched, so update primal, dual and repeat
+ theta=inf;
+ for j=1:head-1
+ t1=queue(j);
+ for rpi=rp(t1):rp(t1+1)-1
+ t2=ci(rpi);
+ if t(t2) == 0 && alpha(t1) + beta(t2) - ai(rpi) < theta,
+ theta = alpha(t1) + beta(t2) - ai(rpi);
+ end
+ end
+ end
+
+ for j=1:head-1, alpha(queue(j)) = alpha(queue(j)) - theta; end
+
+ for j=1:ntmod, beta(tmod(j)) = beta(tmod(j)) + theta; end
+
+ continue;
+ end
+
+ i=i+1; % increment i
+end
+
+val=0;
+for i=1:n
+ for rpi=rp(i):rp(i+1)-1
+ if ci(rpi)==match1(i), val=val+ai(rpi); end
+ end
+end
+noute = 0; % count number of output edges
+for i=1:n
+ if match1(i)<=m, noute=noute+1; end
+end
+m1=zeros(noute,1); m2=m1; % copy over the 0 array
+noute=1;
+for i=1:n
+ if match1(i)<=m, m1(noute)=i; m2(noute)=match1(i);noute=noute+1; end
+end
+
+if nargout>3
+ mi= false(length(tripi)-n,1);
+ for i=1:n
+ for rpi=rp(i):rp(i+1)-1
+ if match1(i)<=m && ci(rpi)==match1(i), mi(tripi(rpi))=1; end
+ end
+ end
+end
+
+
+
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/clustercoeffs.m b/MATLAB_Functions/HuskyTrack/gaimc/clustercoeffs.m
new file mode 100644
index 0000000..19fc771
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/clustercoeffs.m
@@ -0,0 +1,74 @@
+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/MATLAB_Functions/HuskyTrack/gaimc/convert_sparse.m b/MATLAB_Functions/HuskyTrack/gaimc/convert_sparse.m
new file mode 100644
index 0000000..f3fdcfb
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/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 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/MATLAB_Functions/HuskyTrack/gaimc/corenums.m b/MATLAB_Functions/HuskyTrack/gaimc/corenums.m
new file mode 100644
index 0000000..df9ecba
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/corenums.m
@@ -0,0 +1,72 @@
+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/MATLAB_Functions/HuskyTrack/gaimc/csr_to_sparse.m b/MATLAB_Functions/HuskyTrack/gaimc/csr_to_sparse.m
new file mode 100644
index 0000000..3c3e508
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/csr_to_sparse.m
@@ -0,0 +1,54 @@
+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/MATLAB_Functions/HuskyTrack/gaimc/demo/airports.m b/MATLAB_Functions/HuskyTrack/gaimc/demo/airports.m
new file mode 100644
index 0000000..af1970a
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/demo/airports.m
@@ -0,0 +1,119 @@
+%% 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/MATLAB_Functions/HuskyTrack/gaimc/demo/demo.m b/MATLAB_Functions/HuskyTrack/gaimc/demo/demo.m
new file mode 100644
index 0000000..7778ae7
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/demo/demo.m
@@ -0,0 +1,310 @@
+%% 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/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports.html b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports.html
new file mode 100644
index 0000000..c8b66e2
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports.html
@@ -0,0 +1,321 @@
+
+
+
+
+
+
+
+ The US airport network
+
+
+
+
+
+
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.
+
+
+
Contents
+
+
The data
load_gaimc_graph('airports');
+
Plot a histogram for all route time estimates
[si, ti, rt] = find(A);
+hist(-rt,100);
+
Find the lengthiest route
[val,ind] = max(-rt)
+{labels{si(ind)} labels{ti(ind)}}
+
+val =
+
+ 2855
+
+
+ind =
+
+ 63478
+
+
+ans =
+
+ 'Honolulu, HI' 'St. Johns, NL'
+
+
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))
+
+ans =
+
+ 1
+
+
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;
+A = max(A,A');
+T = mst_prim(A);
+clf;
+gplot(T,xy);
+
+
+
+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);
+plotm(Y,X,'k.-','LineWidth',1.5);
+
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);
+
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!
+
+load_gaimc_graph('airports')
+A = -A;
+
+
+[si, ti, rt] = find(A);
+[val,ind] = max(rt);
+start = si(ind);
+dest = ti(ind);
+[d pred] = dijkstra(A,start);
+d(dest)
+
+ans =
+
+ 735
+
+
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');
+
Honolulu, HI --> Chicago, IL --> Montreal, QC --> St. Johns, NL
+
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!
+
+
+
+
+
\ No newline at end of file
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports.png
new file mode 100644
index 0000000..19847a6
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_01.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_01.png
new file mode 100644
index 0000000..b88c23c
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_01.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_02.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_02.png
new file mode 100644
index 0000000..6b650ab
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_02.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_03.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_03.png
new file mode 100644
index 0000000..6d38c92
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_03.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_04.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_04.png
new file mode 100644
index 0000000..3828376
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_04.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_05.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_05.png
new file mode 100644
index 0000000..6112d3f
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/airports_05.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo.html b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo.html
new file mode 100644
index 0000000..05f6e22
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo.html
@@ -0,0 +1,744 @@
+
+
+
+
+
+
+
+ Demo of gaimc - 'Graph Algorithms In Matlab Code'
+
+
+
+
+
+
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.
+
+
+
Contents
+
+
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 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 to an adjacency matrix, we identify each vertex with a row of the matrix via a bijective map. The adjacency matrix is then
+ a matrix called A. The entry A(i,j) = 1 for any edge between in and 0 otherwise. Let's look at an example.
+
load_gaimc_graph('bfs_example');
+graph_draw(A,xy,'labels',labels)
+full(A)
+labels'
+
+ans =
+
+ 0 1 0 0 1 0 0 0
+ 1 0 0 0 0 1 0 0
+ 0 0 0 1 0 1 1 0
+ 0 0 1 0 0 0 0 1
+ 1 0 0 0 0 0 0 0
+ 0 1 1 0 0 0 1 0
+ 0 0 1 0 0 1 0 1
+ 0 0 0 1 0 0 1 0
+
+
+ans =
+
+ 'r' 's' 't' 'u' 'v' 'w' 'x' 'y'
+
+
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');
+whos
+
Name Size Bytes Class Attributes
+
+ A 9x9 416 double sparse
+ Acc 5x5 176 double sparse
+ As 1x1 40792464 struct
+ At 50000x50000 20395704 double sparse
+ D 5x5 200 double
+ D2 5x5 200 double
+ T 456x456 18216 double sparse
+ X 2730x1 21840 double
+ Y 2730x1 21840 double
+ ai 1249731x1 9997848 double
+ aj 71959x1 575672 double
+ ans 1x8 912 cell
+ ati 1249731x1 9997848 double
+ ax 1x1 8 double
+ cc 9x1 72 double
+ cc1 50000x1 400000 double
+ cc2 50000x1 400000 double
+ cc3 50000x1 400000 double
+ cc4 50000x1 400000 double
+ ccfs 2642x1 21136 double
+ ci 1249731x1 9997848 double
+ cn 2642x1 21136 double
+ comp_results 6x4 192 double
+ cp 50001x1 400008 double
+ cs1 1x1 8 double
+ cs2 1x1 8 double
+ cs3 1x1 8 double
+ cs4 1x1 8 double
+ cwd 1x28 56 char
+ d 456x1 3648 double
+ d1 1024x1 8192 double
+ d2 1024x1 8192 double
+ d3 1024x1 8192 double
+ d4 1024x1 8192 double
+ de 71959x1 575672 double
+ dest 1x1 8 double
+ dirtail 1x10 20 char
+ dt 8x1 64 double
+ f 9x1 9 logical
+ faceColors 1x1 1904 struct
+ gi 1x1 8 double
+ graphdir 1x10 20 char
+ graphs 1x5 642 cell
+ i 1x1 8 double
+ ind 1x1 8 double
+ labels 9x1 1026 cell
+ lat 9865x1 78920 double
+ lax 1x1 8 double
+ long 9865x1 78920 double
+ mat_fast 1x1 8 double
+ mat_std 1x1 8 double
+ matlabbgldir 1x20 40 char
+ mex_fast 1x1 8 double
+ mex_std 1x1 8 double
+ n 1x1 8 double
+ nrep 1x1 8 double
+ ntests 1x1 8 double
+ path 1x3 24 double
+ pred 1x456 3648 double
+ rep 1x1 8 double
+ results 1x3 2140 struct
+ ri 1249731x1 9997848 double
+ rp 50001x1 400008 double
+ rst 1x1 8 double
+ rt 71959x1 575672 double
+ si 71959x1 575672 double
+ source 1x83 166 char
+ start 1x1 8 double
+ states 49x1 190442 struct
+ szi 1x1 8 double
+ szs 1x6 48 double
+ te 71959x1 575672 double
+ ti 1x1 8 double
+ u 1x1 8 double
+ v 1x1 8 double
+ val 1x1 8 double
+ xy 9x2 144 double
+
+
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)
+
+d =
+
+ 0
+ 1
+ 3
+ 3
+ 2
+ 4
+ -1
+ -1
+ -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);
+
+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.
+
+
+load_gaimc_graph('airports')
+A = -A;
+lax=247; rst=355;
+
+[d pred] = dijkstra(A,lax);
+
+fprintf('Minimum time: %g\n',d(rst));
+
+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 time: 244
+Path:
+Los Angeles, CA --> Minneapolis/St Paul, MN --> Rochester, MN
+
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
+
+ans =
+
+ (1,1) 24550
+
+
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
+max(scomponents(A))
+
+ans =
+
+ 3
+
+
+max(scomponents(A|A'))
+
+ans =
+
+ 2
+
+
Let's look at the vertices in the strongly connected components
cc = scomponents(A)
+
+cc =
+
+ 2
+ 2
+ 2
+ 1
+ 2
+ 2
+ 3
+ 3
+ 3
+
+
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)
+
+ans =
+
+ (1,1) 2.5034
+
+
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)
+
+
+
+
+ans =
+
+ 0.0160
+
+
Average core numbers
cn = corenums(A);
+mean(cn)
+
+ans =
+
+ 1.9463
+
+
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);
+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
+
Elapsed time is 0.000304 seconds.
+
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
+
Elapsed time is 0.000541 seconds.
+
And just to check, let's make sure the output is the same.
isequal(D,D2)
+
+ans =
+
+ 1
+
+
+
+
+
+
\ No newline at end of file
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo.png
new file mode 100644
index 0000000..6f0c93a
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_01.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_01.png
new file mode 100644
index 0000000..2dbcfe2
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_01.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_02.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_02.png
new file mode 100644
index 0000000..4fd7aee
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_02.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_03.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_03.png
new file mode 100644
index 0000000..9cad9bf
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_03.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_04.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_04.png
new file mode 100644
index 0000000..5e602d1
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_04.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_05.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_05.png
new file mode 100644
index 0000000..e0c2922
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_05.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_06.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_06.png
new file mode 100644
index 0000000..e08d0ee
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_06.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_07.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_07.png
new file mode 100644
index 0000000..d364e9f
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_07.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_08.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_08.png
new file mode 100644
index 0000000..1612ed5
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_08.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq02910.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq02910.png
new file mode 100644
index 0000000..d1fbfd3
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq02910.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq74756.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq74756.png
new file mode 100644
index 0000000..473e14e
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq74756.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq85525.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq85525.png
new file mode 100644
index 0000000..896d5f6
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq85525.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq91421.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq91421.png
new file mode 100644
index 0000000..d1fbfd3
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/demo_eq91421.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple.html b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple.html
new file mode 100644
index 0000000..fbc95aa
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple.html
@@ -0,0 +1,334 @@
+
+
+
+
+
+
+
+ Compare performance of gaimc to matlab_bgl
+
+
+
+
+
+
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.
+
+
+
Contents
+
+
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);
+ ci=components(sparse(ones(5)));
+catch
+ error('gaimc:performance_comparison','Matlab BGL is not working, halting...');
+end
+
Warning: Duplicate directory name: /home/dgleich/dev/matlab-bgl/4.0.
+
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));
+
+ 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;
+
+ scomponents size=1 30
+ scomponents size=10 30
+ scomponents size=100 30
+ scomponents size=5000 30
+ scomponents size=10000 30
+ scomponents size=50000 30
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;
+
dijkstra rep= 30 graph=tapir trial= 100
+
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.
+
+
+
+
+
+
\ No newline at end of file
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple.png
new file mode 100644
index 0000000..a2096d9
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple_01.png b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple_01.png
new file mode 100644
index 0000000..06386ec
Binary files /dev/null and b/MATLAB_Functions/HuskyTrack/gaimc/demo/html/performance_comparison_simple_01.png differ
diff --git a/MATLAB_Functions/HuskyTrack/gaimc/demo/performance_comparison.m b/MATLAB_Functions/HuskyTrack/gaimc/demo/performance_comparison.m
new file mode 100644
index 0000000..f823fd5
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/demo/performance_comparison.m
@@ -0,0 +1,261 @@
+%% 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/MATLAB_Functions/HuskyTrack/gaimc/demo/performance_comparison_simple.m b/MATLAB_Functions/HuskyTrack/gaimc/demo/performance_comparison_simple.m
new file mode 100644
index 0000000..54c8905
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/demo/performance_comparison_simple.m
@@ -0,0 +1,131 @@
+%% 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/MATLAB_Functions/HuskyTrack/gaimc/dfs.m b/MATLAB_Functions/HuskyTrack/gaimc/dfs.m
new file mode 100644
index 0000000..b8bae5c
--- /dev/null
+++ b/MATLAB_Functions/HuskyTrack/gaimc/dfs.m
@@ -0,0 +1,70 @@
+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