This are a demo and instructions for using the graph module, Laplacians for Julia.
pwd()
push!(LOAD_PATH,"src")
using Laplacians
While you do not stricly need the following, it will allow you to plot without using Laplacians routines.
using PyPlot
VERSION
gr = grid2(4)
gr = grid2(4,3)
(x,y) = grid2coords(4,3)
p = plotGraph(gr,x,y;dots=false)
gr = randRegular(15,3)
spy(gr)
gr = grownGraph(200,2)
spy(gr)
gr = full(hyperCube(3))
a = completeBinaryTree(15)
spy(a)
a0 = completeBinaryTree(3)
a1 = completeBinaryTree(5)
a = productGraph(a0,a1)
[eig(full(a))[1]';
sort(kron(eig(full(a0))[1],ones(5)) + kron(ones(3),eig(full(a1))[1]))']
spectralDrawing(a)
The following code computes a 30-by-30 grid graph, samples edges with probability 1/2, and then computes the components.
gr = grid2(30);
grs = subsampleEdges(gr,.5);
co = components(grs)
comps = vecToComps(co)
(x,y) = grid2coords(30,30)
for i = 1:length(comps)
ind = comps[i]
plotGraph(grs[ind,ind],x[ind],y[ind],rand(3);dots=false,setaxis=false)
end
pm = plot(x,y,marker="o",linestyle="none",color="black",markersize=2)
We could read in an edge list, either in ij or ijv form as shown here. The following is for an ij, and will be used to test the speed of our connected components code. I generated this graph in matlab, and then saved it using
[asi,asj] = find(triu(as)); dlmwrite('testGraph3.txt',[asi,asj],'precision',9);Note that, without the precision, we won't get all the digits in the vertex ids.
The following is the steps that appear in the function that I have since named readIJ
edlist = readdlm("testGraph3.txt",',')
Let's try to make a sparse matrix with these as entries. Somehow, we will have to first convert to Integer types.
n = maximum(edlist)
m = size(edlist)
edlist = convert(Array{Int64,2}, edlist)
a = sparse(edlist[:,1],edlist[:,2],ones(m[1]),n,n)
a = a + a'
a = readIJ("testGraph3.txt")
@time co = components(a);
That was 0.23 seconds. In comparison, calling the Lapsolver routine from matlab for this took longer:
>> gr = a2g(a); >> gu = GraphUtils >> tic; co = gu.getComponents(gr); toc Elapsed time is 0.286240 seconds.And, matlab_bgl is even slower
>> addpath ~/progs/matlab_bgl/ >> tic; [comp,sz] = components(a); toc Elapsed time is 1.142663 seconds.
Since you probably don't have that graph, let me repeat this experiment with percolation on the grid, again in Matlab and in Julia.
>> a = grid2(1000); >> [ai,aj,av] = find(tril(a)); >> keep = (rand(size(av))>.5); >> asamp = sparse(ai(keep), aj(keep), av(keep), 1000000, 1000000); >> asamp = asamp + asamp'; >> gr = a2g(asamp); >> gu = GraphUtils; >> tic; co = gu.getComponents(gr); toc Elapsed time is 0.142004 seconds. >> addpath ~/progs/matlab_bgl/ >> tic; [comp,sz] = components(asamp); toc Elapsed time is 0.238816 seconds.
a = grid2(1000);
asamp = subsampleEdges(a,.5);
@time co = components(asamp);
In case that data disappears somehow, that was:
If we want to use less memory, we can convert the indices in the graph to Int32 from Int64. I would have expected this to speed up the code, but it does not seem to do so.
ashort = shortIntGraph(asamp)
@time co = components(ashort);
I will now compare the running time of two versions of shortest paths code that I wrote in Julia against matlab_bgl and my java code. The "slow" version uses Julia's own priority queues. The faster version uses my own priority queue called nodeHeap, and is an exact transcription of my java code. It is about 10x faster! It could be sped up some more, but this proves the point.
a = readIJ("testGraph3.txt");
@time dists, pArray = shortestPathsSlow(a,1)
@time dists2, pArray2 = shortestPaths(a,1)
Here are the results in Matlab.
>> dl = dlmread('testGraph3.txt'); a = sparse(dl(:,1),dl(:,2),1); n = max(max(dl)) a(n,n) = 0; a = a + a'; >> [ai,aj,av] = find(tril(a)); pg = graphs.PreGraph(ai,aj,av) >> tic; pa = pg.shortestPathTree(0); toc Elapsed time is 1.428569 seconds. tic; [d pa] = shortest_paths(a,1); toc Elapsed time is 1.795666 seconds.
This is a very clear win for Julia!
Let's fight against one other random graph distribution, just to check.
a = randRegular(1000000,3); [ai,aj,av] = find(a); ar = sparse(ai,aj,rand(size(av))); ar2 = sparse(ai,aj,1./rand(size(av))); tic; [d pa] = shortest_paths(ar,1); toc Elapsed time is 3.066359 seconds. tic; [d pa] = shortest_paths(ar2,1); toc Elapsed time is 3.124175 seconds.pg = graphs.PreGraph(ai,aj,av); tic; pa = pg.shortestPathTree(0); toc Elapsed time is 3.151916 seconds. [ai,aj,av] = find(tril(ar2)); pg = graphs.PreGraph(ai,aj,av); tic; pa = pg.shortestPathTree(0); toc Elapsed time is 2.916280 seconds. </pre> So, the time to beat is 3 seconds. Julia does it in just a little over 2, as we will now see.
a = randRegular(1000000,3);
(ai,aj,av) = findnz(triu(a));
ar = sparse(ai,aj,rand(length(av)),1000000,1000000);
ar = ar + ar';
ar2 = sparse(ai,aj,1./rand(length(av)),1000000,1000000);
ar2 = ar2 + ar2';
@time dists2, pArray2 = shortestPaths(ar,1)
@time dists2, pArray2 = shortestPaths(ar2,1);
We will also compare the running time of mst code. In this case, I will do it with a 2d grid with randomly chosen edge weights. We will see that the Julia code is over twice as fast as the matlab_bgl and java code that we have.
a = grid2(1000)
n = size(a)[1]
(ai,aj,av) = findnz(triu(a))
ar = sparse(ai,aj,rand(size(av)),n,n)
ar = ar + ar'
@time tree = kruskal(ar);
>> a = grid2(1000); >> n = length(a); >> [ai,aj,av] = find(triu(a)); >> ar = sparse(ai,aj,rand(length(av),1),n,n); ar = ar + ar'; >> tic; t = kruskal_mst(ar); toc Elapsed time is 6.501457 seconds. >> gr = a2g(ar); >> kt = KruskalTree; >> tic; tr = kt.getTree(gr); toc Elapsed time is 6.141985 seconds.
By default, this computes the minimum spanning tree. To get the max, do this.
@show sum(triu(tree).nzval)
maxTree = kruskal(ar,kind="max")
@show sum(triu(maxTree).nzval)
la = lap(grid2(3))
E = eigs(la, nev = 3, which=:SR)
V = E[2]
You would think that you should use
E = eigs(la, nev = 3, which=:SM)But, that gives horrible results.
plotGraph(la,V[:,2],V[:,3])
a = hyperCube(3)
la = lap(a)
E = eigs(la, nev = 3, which=:SR)
V = E[2]
plotGraph(a,V[:,2],V[:,3])
a = grid2(5)
typeof(a)
fieldnames(SparseMatrixCSC)