first cut

Author: sbromberger

REQUIRE

@@ -1 +1,2 @@
 julia 0.6
+Clustering

src/CommunityDetection.jl

@@ -1,5 +1,103 @@
+__precompile__(true)
 module CommunityDetection
+using LightGraphs
+import Clustering: kmeans
 
-# package code goes here
+export community_detection_nback, community_detection_bethe
 
-end # module
+"""
+    community_detection_nback(g::AbstractGraph, k::Int)
+
+Return an array, indexed by vertex, containing commmunity assignments for
+graph `g` detecting `k` communities.
+Community detection is performed using the spectral properties of the 
+non-backtracking matrix of `g`.
+
+### References
+- [Krzakala et al.](http://www.pnas.org/content/110/52/20935.short)
+"""
+function community_detection_nback(g::AbstractGraph, k::Int)
+    #TODO insert check on connected_components
+    ϕ = real(nonbacktrack_embedding(g, k))
+    if k==2
+        c = community_detection_threshold(g, ϕ[1,:])
+    else
+        c = kmeans(ϕ, k).assignments
+    end
+    return c
+end
+
+function community_detection_threshold(g::AbstractGraph, coords::AbstractArray)
+    # TODO use a more intelligent method to set the threshold
+    # 0 based thresholds are highly sensitive to errors.
+    c = ones(Int, nv(g))
+    # idx = sortperm(λ, lt=(x,y)-> abs(x) > abs(y))[2:k] #the second eigenvector is the relevant one
+    for i=1:nv(g)
+        c[i] = coords[i] > 0 ?  1 : 2
+    end
+    return c
+end
+
+
+"""
+	nonbacktrack_embedding(g::AbstractGraph, k::Int)
+
+Perform spectral embedding of the non-backtracking matrix of `g`. Return
+a matrix ϕ where ϕ[:,i] are the coordinates for vertex i.
+
+### Implementation Notes
+Does not explicitly construct the `non_backtracking_matrix`.
+See `Nonbacktracking` for details.
+
+### References
+- [Krzakala et al.](http://www.pnas.org/content/110/52/20935.short).
+"""
+function nonbacktrack_embedding(g::AbstractGraph, k::Int)
+    B = Nonbacktracking(g)
+    λ, eigv, conv = eigs(B, nev=k+1, v0=ones(Float64, B.m))
+    ϕ = zeros(Complex64, nv(g), k-1)
+    # TODO decide what to do with the stationary distribution ϕ[:,1]
+    # this code just throws it away in favor of eigv[:,2:k+1].
+    # we might also use the degree distribution to scale these vectors as is
+    # common with the laplacian/adjacency methods.
+    for n=1:k-1
+        v= eigv[:,n+1]
+        ϕ[:,n] = contract(B, v)
+    end
+    return ϕ'
+end
+
+
+
+"""
+    community_detection_bethe(g::AbstractGraph, k=-1; kmax=15)
+
+Perform detection for `k` communities using the spectral properties of the 
+Bethe Hessian matrix associated to `g`.
+If `k` is omitted or less than `1`, the optimal number of communities
+will be automatically selected. In this case the maximum number of
+detectable communities is given by `kmax`.
+Return a vector containing the vertex assignments.
+
+### References
+- [Saade et al.](http://papers.nips.cc/paper/5520-spectral-clustering-of-graphs-with-the-bethe-hessian)
+"""
+function community_detection_bethe(g::AbstractGraph, k::Int=-1; kmax::Int=15)
+    A = adjacency_matrix(g)
+    D = diagm(degree(g))
+    r = (sum(degree(g)) / nv(g))^0.5
+
+    Hr = (r^2-1)*eye(nv(g))-r*A+D;
+    # Hmr = (r^2-1)*eye(nv(g))+r*A+D;
+    k >= 1 && (kmax = k)
+    λ, eigv = eigs(Hr, which=:SR, nev=min(kmax, nv(g)))
+    q = findlast(x -> x<0, λ)
+    k > q && warn("Using eigenvectors with positive eigenvalues,
+                    some communities could be meaningless. Try to reduce `k`.")
+    k < 1 && (k = q)
+    k < 1 && return fill(1, nv(g))
+    labels = kmeans(eigv[:,2:k]', k).assignments
+    return labels
+end
+
+end #module

test/runtests.jl

@@ -1,5 +1,93 @@
 using CommunityDetection
+using LightGraphs
 using Base.Test
 
-# write your own tests here
-@test 1 == 2
+""" Spectral embedding of the non-backtracking matrix of `g`
+(see [Krzakala et al.](http://www.pnas.org/content/110/52/20935.short)).
+
+`g`: input Graph
+`k`: number of dimensions in which to embed
+
+return : a matrix ϕ where ϕ[:,i] are the coordinates for vertex i.
+"""
+function nonbacktrack_embedding_dense(g::AbstractGraph, k::Int)
+    B, edgeid = non_backtracking_matrix(g)
+    λ,eigv,conv = eigs(B, nev=k+1, v0=ones(Float64, size(B,1)))
+    ϕ = zeros(Complex64, k-1, nv(g))
+    # TODO decide what to do with the stationary distribution ϕ[:,1]
+    # this code just throws it away in favor of eigv[:,2:k+1].
+    # we might also use the degree distribution to scale these vectors as is
+    # common with the laplacian/adjacency methods.
+    for n=1:k-1
+        v= eigv[:,n+1]
+        for i=1:nv(g)
+            for j in neighbors(g, i)
+                u = edgeid[Edge(j,i)]
+                ϕ[n,i] += v[u]
+            end
+        end
+    end
+    return ϕ
+end
+
+n = 10; k = 5
+pg = PathGraph(n)
+ϕ1 = CommunityDetection.nonbacktrack_embedding(pg, k)'
+
+nbt = Nonbacktracking(pg)
+B, emap = non_backtracking_matrix(pg)
+Bs = sparse(nbt)
+@test sparse(B) == Bs
+
+# check that matvec works
+x = ones(Float64, nbt.m)
+y = nbt * x
+z = B * x
+@test norm(y-z) < 1e-8
+
+#check that matmat works and full(nbt) == B
+@test norm(nbt*eye(nbt.m) - B) < 1e-8
+
+#check that we can use the implicit matvec in nonbacktrack_embedding
+@test size(y) == size(x)
+ϕ2 = nonbacktrack_embedding_dense(pg, k)'
+@test size(ϕ2) == size(ϕ1)
+
+#check that this recovers communities in the path of cliques
+n=10
+g10 = CompleteGraph(n)
+z = copy(g10)
+for k=2:5
+    z = blkdiag(z, g10)
+    add_edge!(z, (k-1)*n, k*n)
+
+    c = community_detection_nback(z, k)
+    @test sort(union(c)) == [1:k;]
+    a = collect(n:n:k*n)
+    @test length(c[a]) == length(unique(c[a]))
+    for i=1:k
+        for j=(i-1)*n+1:i*n
+            @test c[j] == c[i*n]
+        end
+    end
+
+    c = community_detection_bethe(z, k)
+    @test sort(union(c)) == [1:k;]
+    a = collect(n:n:k*n)
+    @test length(c[a]) == length(unique(c[a]))
+    for i=1:k
+        for j=(i-1)*n+1:i*n
+            @test c[j] == c[i*n]
+        end
+    end
+
+    c = community_detection_bethe(z)
+    @test sort(union(c)) == [1:k;]
+    a = collect(n:n:k*n)
+    @test length(c[a]) == length(unique(c[a]))
+    for i=1:k
+        for j=(i-1)*n+1:i*n
+            @test c[j] == c[i*n]
+        end
+    end
+end