hkrelax.py

"""
hkrelax.py

A demonstration of a relaxation method for computing a heat-kernel based
community that implements the algorith from "Heat-kernel based community
detection" by Kloster & Gleich.


Written by Kyle Kloster and David F. Gleich
"""

import collections
import sys
import math
 
# setup the graph
G = {
  1:set([ 2, 3, 5, 6,]),
  2:set([ 1, 4,]),
  3:set([ 1, 6, 7,]),
  4:set([ 2, 5, 7, 8,]),
  5:set([ 1, 4, 6, 8, 9, 10,]),
  6:set([ 1, 3, 5, 7,]),
  7:set([ 3, 4, 6, 9,]),
  8:set([ 4, 5, 9,]),
  9:set([ 5, 7, 8, 20,]),
  10:set([ 5, 11, 12, 14, 15,]),
  11:set([ 10, 12, 13, 14,]),
  12:set([ 10, 11, 13, 14, 15,]),
  13:set([ 11, 12, 15,]),
  14:set([ 10, 11, 12, 25,]),
  15:set([ 10, 12, 13,]),
  16:set([ 17, 19, 20, 21, 22,]),
  17:set([ 16, 18, 19, 20,]),
  18:set([ 17, 20, 21, 22,]),
  19:set([ 16, 17,]),
  20:set([ 9, 16, 17, 18,]),
  21:set([ 16, 18,]),
  22:set([ 16, 18, 23,]),
  23:set([ 22, 24, 25, 26, 27,]),
  24:set([ 23, 25, 26, 27,]),
  25:set([ 14, 23, 24, 26, 27,]),
  26:set([ 23, 24, 25,]),
  27:set([ 23, 24, 25,]),
}
Gvol = 102

def compute_psis(N,t):
    psis = {}
    psis[N] = 1.
    for i in xrange(N-1,0,-1):
        psis[i] = psis[i+1]*t/(float(i+1.))+1.
    return psis
    

## Setup parameters that can be computed automatically
N = 11 # see paper for how to set this automatically
t = 5.
eps = 0.01
seed = [1]
psis = compute_psis(N,t)
    
## Estimate hkpr vector 
# G is graph as dictionary-of-sets, 
# t, tol, N, psis are precomputed
x = {} # Store x, r as dictionaries
r = {} # initialize residual
Q = collections.deque() # initialize queue
for s in seed: 
  r[(s,0)] = 1./len(seed)
  Q.append((s,0))
while len(Q) > 0:
  (v,j) = Q.popleft() # v has r[(v,j)] ...
  rvj = r[(v,j)]
  # perform the hk-relax step
  if v not in x: x[v] = 0.
  x[v] += rvj 
  r[(v,j)] = 0.
  mass = (t*rvj/(float(j)+1.))/len(G[v])
  for u in G[v]:   # for neighbors of v
    next = (u,j+1) # in the next block
    if j+1 == N: 
      x[u] += rvj/len(G[v])
      continue
    if next not in r: r[next] = 0.
    thresh = math.exp(t)*eps*len(G[u])
    thresh = thresh/(N*psis[j+1])
    if r[next] < thresh and \
       r[next] + mass >= thresh:
       Q.append(next) # add u to queue
    r[next] = r[next] + mass

# Find cluster, first normalize by degree
for v in x: x[v] = x[v]/len(G[v])  
  
for v in xrange(1,len(G)+1):
    if v in x:
        print "hk[%2i] = %.16lf"%(v, x[v])
    else:
        print "hk[%2i] = -0."%(v)
print 
    
## Step 2 do a sweep cut based on this vector 
  
# now sort x's keys by value, decreasing
sv = sorted(x.iteritems(), key=lambda x: x[1], reverse=True)
S = set()
volS = 0.
cutS = 0.
bestcond = 1.
bestset = sv[0]
for p in sv:
  s = p[0] # get the vertex
  volS += len(G[s]) # add degree to volume
  for v in G[s]:
    if v in S:
      cutS -= 1
    else:
      cutS += 1
  print "v: %4i  cut: %4f  vol: %4f"%(s, cutS,volS)
  S.add(s)
  if cutS/min(volS,Gvol-volS) < bestcond:
    bestcond = cutS/min(volS,Gvol-volS)
    bestset = set(S) # make a copy
print "Best set conductance: %f"%(bestcond)
print "  set = ", str(bestset)