#### Incorporating additional information into the formation of clusters.

[Covariate assisted spectral clustering (CASC)](http://arxiv.org/abs/1411.2158) is a way of including additional node level data to bias the clusters towards something that we suspect matters.  In this case, we want clusters that are "continuous in space".  

Below, I have not included the literal version of CASC.  Instead, the approach is as follows.  

1) Construct two graphs.   The first graph is the same as above.  The second graph is the 10 nearest neighbor graph for zip code locations.  
2) Compute the regularized $L_{sym}$ matrix for both graphs. Add these matrices to form a new similarity matrix. 
3) Compute the leading eigenvectors of this matrix and row normalize.  Run k-means.


```{r include = FALSE}

# this loads the data, the libraries, and a function to draw a convex hull.

library(Matrix)
library(igraph)
library(data.table)
library(rARPACK)
library(zipcode)
library(maps)
library(grDevices)

set.seed(1)
rm(list = ls())

data(zipcode); zipcode = as.data.table(zipcode); setkey(zipcode, zip)

# This is copy pasted from https://chitchatr.wordpress.com/2011/12/30/convex-hull-around-scatter-plot-in-r/
### Plotting function to plot convex hulls
### Filename: Plot_ConvexHull.R
### Notes:
############################################################################

# INPUTS:
# xcoords: x-coordinates of point data
# ycoords: y-coordinates of point data
# lcolor: line color

# OUTPUTS:
# convex hull around data points in a particular color (specified by lcolor)

# FUNCTION:
Plot_ConvexHull<-function(xcoord, ycoord, lcolor, thick=1){
  hpts <- chull(x = xcoord, y = ycoord)
  hpts <- c(hpts, hpts[1])
  lines(xcoord[hpts], ycoord[hpts], col = lcolor, lwd = thick)
  }  
# END OF FUNCTION
```


```{r}
load(url("http://pages.stat.wisc.edu/~karlrohe/netsci/data/zipA.RData"))
library(FNN)
dat = as.matrix(zipcode[rownames(A),c("longitude", "latitude"), with = F])
good = complete.cases(dat)
Knng = get.knn(dat[good,])
dict = rownames(A)[good]
cbind(dict[1], dict[Knng$nn.index[1,]])
knnEL = matrix("", nrow = prod(dim(Knng$nn.index)), ncol = 2)
for(i in 1:nrow(Knng$nn.index))  knnEL[1:10 + 10*(i-1),] = cbind(dict[i], dict[Knng$nn.index[i,]])
gknn=graph.edgelist(knnEL,directed = F) # this creates a graph.
Aknn = get.adjacency(gknn)
cross = match(rownames(A)[good], rownames(Aknn))
mean(rownames(Aknn)[cross] == rownames(A)[good])
Aknn = Aknn[cross,cross]

# this loads a function called "sumGraphs" that takes to adjacency matrices, converts them to symmetric Laplacians, and adds these matrices together.  
source(url("http://pages.stat.wisc.edu/~karlrohe/netsci/code/globalTransitivityFunctions.R"))
CASC = function(A, B, k){
  X = eigs(sumGraphs(A,B),k)$vec
  X = t(apply(X, 1, function(x) return(x/sqrt(sum(x^2)))))
  return(kmeans(X, k, nstart = 20))
  }

csc = CASC(A,Aknn, 10)
sort(table(csc$clust))
loc = zipcode[rownames(A), c("longitude", "latitude"), with  =F]
loc = as.matrix(loc)
# pdf(file = "newRegions2.pdf", height = 13, width = 20)
map('state', add = F) 
# draw a polygon around the convex hull of the first 10 clusters.
good = complete.cases(loc)
comploc = loc[good,]
points(rnorm(comploc[,1],comploc[,1],.1),rnorm(comploc[,2],comploc[,2],.1), col=as.factor(csc$clust[good]), pch  =19)

for(i in 1:10){
  i = i+1
  inclust = which(csc$cluster[good] ==i)
  Plot_ConvexHull(xcoord = comploc[inclust,1], ycoord = comploc[inclust,2],lcolor = "black", thick = 2)
  }
# dev.off()

```

```{r}
k = 250
csc = CASC(A,Aknn, k)
as.vector(sort(table(csc$clust)))
plot(as.vector(sort(table(csc$clust))))
loc = zipcode[rownames(A), c("longitude", "latitude"), with  =F]
loc = as.matrix(loc)
# pdf(file = "newRegions2.pdf", height = 13, width = 20)
plot(rnorm(loc[,1],loc[,1],.1),rnorm(loc[,2],loc[,2],.1), col=as.factor(csc$clust),
     xlim = c(-125, -65), ylim = c(24, 50), pch  =19)
map('state', add = T) 
# draw a polygon around the convex hull of the first 10 clusters.
good = complete.cases(loc)
comploc = loc[good,]
for(i in 1:k){
  i = i+1
  inclust = which(csc$cluster[good] ==i)
  Plot_ConvexHull(xcoord = comploc[inclust,1], ycoord = comploc[inclust,2],lcolor = "black", thick = 2)
  }
# dev.off()

map('state', add = F) 
# draw a polygon around the convex hull of the clusters
good = complete.cases(loc)
comploc = loc[good,]
for(i in 1:k){
  inclust = which(csc$cluster[good] ==i)
  Plot_ConvexHull(xcoord = comploc[inclust,1], ycoord = comploc[inclust,2],lcolor = "black", thick = 2)
  }

# focus on the North East
# pdf(file = "northEast.pdf", height = 6, width= 6)
NorthEast = c("Connecticut", "Maine", "Massachusetts", "New Hampshire", "Rhode Island", "Vermont", "New Jersey", "New York", "Pennsylvania", "Mary", "Del")
map('county', NorthEast, add = F) 
good = complete.cases(loc)
comploc = loc[good,]
points(rnorm(comploc[,1],comploc[,1],.03),rnorm(comploc[,2],comploc[,2],.03), col=as.factor(csc$clust[good]), pch = ".")
for(i in 1:k){
  i = i+1
  inclust = which(csc$cluster[good] ==i)
  Plot_ConvexHull(xcoord = comploc[inclust,1], ycoord = comploc[inclust,2],lcolor = "black", thick = 2)
  }
# dev.off()
```


What is going on with the crazy large clusters?  Let's find them.

```{r}
# given a set of 2d points, compute the sum of squares from the center
sd2 = function(x){
  x = x[complete.cases(x),]
  return( sqrt(sum(apply(scale(x, scale = F), 2, var))))
  }
# given a set of 2d points, compute the maximal distance from the center
max2 = function(x){
  x = x[complete.cases(x),]
  return( max(apply(scale(x, scale = F), 1, function(z) return(sqrt(sum(z^2))))))
  }
sizeSD = rep(NA, k)
sizeMAX = sizeSD
num = sizeSD
for(i in 1:k){
  inclust = which(csc$cluster[good] ==i)
  sizeSD[i] = sd2(comploc[inclust,])
  sizeMAX[i] = max2(comploc[inclust,])
  num[i]=length(inclust)
  }
plot(log(sizeSD,10), log(sizeMAX,10))
plot(num, log(sizeMAX))
```

Define the big clusters to be the ones with log(sizeSD,10) > .5. 

```{r}
bigClust = which(log(sizeMAX,10)>.5)

# draw a polygon around the big clusters
map('world', add = F) 
for(i in bigClust){
  inclust = which(csc$cluster[good] ==i)
  Plot_ConvexHull(xcoord = comploc[inclust,1], ycoord = comploc[inclust,2],lcolor = "red", thick = 2)
  }

map('state', add = F) 
for(i in bigClust){
  inclust = which(csc$cluster[good] ==i)
  Plot_ConvexHull(xcoord = comploc[inclust,1], ycoord = comploc[inclust,2],lcolor = "black", thick = 2)
  }
```


The really large cluster on the map is also is the cluster with the largest number of nodes.
```{r}
map('state', add = F) 
i = which.max(num)
inclust = which(csc$cluster[good] ==i)
Plot_ConvexHull(xcoord = comploc[inclust,1], ycoord = comploc[inclust,2],lcolor = "black", thick = 2)
points(rnorm(comploc[inclust,1],comploc[inclust,1],.03),rnorm(comploc[inclust,2],comploc[inclust,2],.03))

```
This cluster contains 165 zip codes; roughly 1% of the zip codes.  Note: not all zip codes that fall within the plotted region belong to the cluster!  

Perhaps we simply need to increase the number of clusters! 
```{r}
k = 350
csc = CASC(A,Aknn, k)
loc = zipcode[rownames(A), c("longitude", "latitude"), with  =F]
loc = as.matrix(loc)
# pdf(file = "newRegions2.pdf", height = 13, width = 20)
plot(rnorm(loc[,1],loc[,1],.1),rnorm(loc[,2],loc[,2],.1), col=as.factor(csc$clust),
     xlim = c(-125, -65), ylim = c(24, 50), pch  =19)
map('state', add = T) 
# draw a polygon around the convex hull of the first 10 clusters.
good = complete.cases(loc)
comploc = loc[good,]
for(i in 1:k){
  i = i+1
  inclust = which(csc$cluster[good] ==i)
  Plot_ConvexHull(xcoord = comploc[inclust,1], ycoord = comploc[inclust,2],lcolor = "black", thick = 2)
  }
```
Still big clusters!  One possibility is to "tune" the CASC weighting of the knn graph and the referral network.  The next lecture tries a vastly different algorithm. 

