# Spectral methods and DBSCAN clustering

This workshop explores spectral clustering on cyber data using DBSCAN. Because of the peculiarities of this data as visualised in [Lecture 3.1](https://dsbristol.github.io/dst/assets/slides/03.1-LatentStructuresPCA.pdf) using Spectral decomposition and clustered in [Lecture 3.2](https://dsbristol.github.io/dst/assets/slides/03.2-Clustering.pdf), it generates strange structure that can only be handled with carefully chosen methods. 

Whilst not all data are as extreme as this example, the methods discussed here are robust and apply in many other contexts.

Remember to look at the code examples from the lectures for analogous data for other methods.

## 0. Requirements

```{r}
if(!require("purrr")) install.packages("purrr") # allows the nice %>% notation; part of the "tidyverse"
if(!require("dbscan")) install.packages("dbscan") # density based clustering
if(!require("readr")) install.packages("readr") # for read_csv which is improved over read.csv
library("purrr")
library("readr")
```

## 1. Data

We will continue with the KDD11 dataset from the first workshop. Recall that we downloaded it with:

```{r}
if(!dir.exists(file.path("..","data"))) dir.create(file.path("..","data"))
if(!file.exists(file.path("..","data","kddcup.data_10_percent.zip"))) download.file("http://kdd.org/cupfiles/KDDCupData/1999/kddcup.data_10_percent.zip", destfile=file.path("..","data","kddcup.data_10_percent.zip"))
if(!file.exists(file.path("..","data","kddcup.names"))) download.file("http://kdd.org/cupfiles/KDDCupData/1999/kddcup.names",destfile=file.path("..","data","kddcup.names"))
```

We'll read it in as before.

```{r}
kddata<-as.data.frame(read_csv(file.path("..","data","kddcup.data_10_percent.zip"),col_names=FALSE)) ## Ignore the warnings - there is a bug with the header
kddnames=read.table(file.path("..","data","kddcup.names"),sep=":",skip=1,as.is=T)
colnames(kddata)=c(kddnames[,1],"normal") # here we fix the bug with the header
goodcat=names(which(table(kddata[,"normal"])>1))
kddata=kddata[kddata[,"normal"]%in%goodcat,]
```

### 1.1 Data transformations

Now we make the data suitable for spectral embedding. There are several motivations here:

a) You can take an SVD of any matrix. However, any reasonable interpretation of SVD assumes **additive noise**. When the data are inherently multiplicative, we need to take a logarithm for this assumption hold. Here the data are heavy tailed, an indication of some sort of multiplicative noise.
b) Much of these data are categorical. Whilst we could use categorical data for SVD, it would need to be one-hot encoded and scaled appropriately. You should do that for assignments but its overkill here.
c) We did EDA on the data and found several columns to contain hard-to-interpret structure. These do too, but we'll be able to do the work of exploring them.

We will consider a few transformations here. You should rerun the analysis with different choices to get an understanding of how they effect the inference.

```{r}
testdata=kddata[,c("src_bytes","dst_bytes",
                   "count","srv_count",
                     "dst_host_count","dst_host_srv_count")]
#transformation=function(x)x # Identity
#transformation=function(x)(x-mean(x))/sd(x) # standardize and center
#transformation=function(x)x/sd(x) # standardize only
transformation=function(x){xx=log(1+x);(xx)/sd(xx)} # log transform and then standardize
# transformation=function(x)log(1+x) # log transform
testdatatrans=apply(testdata,2,transformation)
```

### 1.2 Spectral Embedding

Make a suitable "spectral embedding" (SVD or PCA) from the transformed test data, and plot the Eigenvalues to choose the number of PCs.

```{r}
testdatatrans.svd <- svd(testdatatrans)
```

Plotting the eigenvalues (a "scree plot"):
```{r}
plot(testdatatrans.svd$d,xlab="Eigenvalue index",ylab="Eigenvalue",log="")
```

#### See Q1 in Block 3 Portfolio

about exploring the different transformations that are possible and how they impact inference.

With this transformation choice it looks like we want the first 3 PCs, which we will record in a variable:

```{r}
npcs=3
```

Extensions: Try rerunning the following analysis using more or fewer PCs. How does the choice of representation effect our analysis - what if we scale differently? What if we normalise features?

#### See Q2 in Block 3 Portfolio

about the correct number of PCs to retain.

## 2.1 Visualising Eigenvectors

We make a plot of the data, annotated by the "normal" category that it has been assigned. 
Because there are very many examples of some classes of data, we retain only up to 5000 of each class.

First we will obtain class labels the indices of different classes of data, and order them by the number of each class:
```{r}
normaltab=table(kddata[,"normal"])
normallist=lapply(names(normaltab),function(x)which(kddata[,"normal"]==x))
normalsamples=lapply(normallist,function(x){
  if(length(x)<5000) return(x)
  return(sample(x,size = 5000))
})
names(normalsamples)=names(normaltab)
normalsizes=sapply(normalsamples,length)
normalsamples=normalsamples[order(normalsizes,decreasing=TRUE)]
normalsizes=normalsizes[order(normalsizes,decreasing=TRUE)]
normalsizes
```

From this we can see that we have downsampled `neptune`, `normal`, and `smurf` classes. We have very few cases of several attacks.

Plot the PCA:
```{r}
i=2;j=3
plot(testdatatrans.svd$u[,i],
     testdatatrans.svd$u[,j],type="n",xlab=paste0("PC",i),ylab=paste0("PC",j),
     col="#33333311",pch=19,cex=0.3) ## Plot where the points are omitted
legend("topleft",legend=names(normalsamples),text.col=1:length(normalsamples),col=1:length(normalsamples),pch=c(19,1,2)[ceiling((1:30)/10)],cex=0.5)
for(index in 1:length(normalsamples)){ ## Question: why would we plot the points in decreasing size order?
  points(testdatatrans.svd$u[normalsamples[[index]],c(i,j)],
         col=index+1,pch=c(19,1,2)[ceiling((1:30)/10)],cex=0.3)
}
```

This looks reassuring that there is structure to find for classification. The "line" structures are potentially interesting, and contrast to the "blob" structures in the larger clusters.


#### See Q3 in Block 3 Portfolio

about the interpretation of these results.

## 3.1 Clustering with DBscan

Here we perform parameter tuning for clustering using "dbscan". 

First, we use the function "kNNdist" which obtains the k nearest neighbours. We will explore this to test what sort of thresholds for the parameter "eps" might be useful.

```{r}
library("dbscan")
test=kNNdist(testdatatrans.svd$u[,1:npcs], k = 5,all=TRUE)
testmin=apply(test,1,min)
range(testmin)
```

Now we'll plot all of the minimum values. This gives the distance to the closest other point, for every point.

```{r}
plot(sort(testmin[testmin>1e-8]),log="y",type="l")
threshholds= c(0.01,0.001,0.0001,0.00001,0.000001)
abline(h=c(0.01,0.001,0.0001,0.00001,0.000001))
abline(h=0.0001,col="red") # we chose this number.
## abline(h=0.001) # would give bigger clusters
## abline(h=0.00001) # would give smaller clusters
```

From this analysis, we can see that nearly all points are within 1e-04 of another point.

Extension: Experiment with other thresholds and see how the clustering behaves.

#### See Q4 in Block 3 Portfolio

about running DBSCAN.

#### 3.1.1 Downsampling

Unfortunately, dbscan runs out of memory when applied to the full dataset. Here we have 

```{r}
dim(testdatatrans.svd$u)[1]
``` 

so a small amount of downsampling should solve the problem. Make a list of indices that can be used to downsample the main dataset, keeping all examples of small classes and up to 60000 of the larger classes. You should get 188752 entries. Then use this clustering to run dbscan.
```{r}
set.seed(1)
normalsamples2=lapply(normallist,function(x){
  if(length(x)<60000) return(x)
  return(sample(x,size = 60000))
})
myindices=unlist(normalsamples2) %>% sort
length(myindices)
```

### 3.1.2 DBScan clustering

Now do the clustering:
```{r}
dbscanres=dbscan(testdatatrans.svd$u[myindices,1:3],0.0001)
```

We can examine the results in a simple way using "table", to ask: How many classes are recovered? How do they associate with the labels?

```{r}
table(dbscanres$cluster)
```

Now we will display this with a heatmap:
```{r}
##pdf("test_kddspectralembedding_dbscan.pdf",height=10,width=10)
    plot(testdatatrans.svd$u[myindices,1],
         testdatatrans.svd$u[myindices,2],
         xlab="PC1",ylab="PC2",main=paste("K=41 + noise"),
         col=c("#66666666",rainbow(41))[dbscanres$cluster+1],pch=19,cex=0.5)
    plot(testdatatrans.svd$u[myindices,1],
         testdatatrans.svd$u[myindices,3],xlab="PC1",ylab="PC3",main=paste("K=41 + noise"),
         col=c("#66666666",rainbow(41))[dbscanres$cluster+1],pch=19,cex=0.5)
    plot(testdatatrans.svd$u[myindices,2],
         testdatatrans.svd$u[myindices,3],xlab="PC2",ylab="PC3",,main=paste("K=41 + noise"),
         col=c("#66666666",rainbow(41))[dbscanres$cluster+1],pch=19,cex=0.5)
##dev.off()
```

Examining these plots, it is clear that there are some important distinctions between the inferred clusters. We can make a table of the labels and clustering:

```{r}
mytable=data.frame(cluster=dbscanres$cluster,label=kddata[myindices,"normal"])
mytab=table(mytable)
myfreq=mytab/rowSums(mytab)
heatmap(t(myfreq),scale="none",main="Frequency of label given cluster",cexRow=0.7)
```


```{r}
myfreq2=t(t(mytab)/colSums(mytab))
heatmap(t(myfreq2),scale="none",main="Frequency of cluster given label",cexRow = 0.7)
```

### 3.1.3 Interpretation as a classifier

How well does the clustering describes each label?

```{r}
bestmatch=cbind(bestfreq=apply(myfreq,1,max),count=apply(mytab,1,sum))
print(bestmatch)
accuracy=sum(bestmatch[,1]*bestmatch[,2])/sum(bestmatch[,2])
print(paste("Accuracy in training data:",accuracy))
```

Overall:

* The mean accuracy is low, but the accuracy for every cluster except 0 and 3 is high.
  * Cluster 3 is very large.
  * There are very many clusters that are close to perfectly associated with a label. 
  * Even the "unclustered" data (0) is only very small. 
* Therefore, with any other assignment, we can be more confident.
  * "Cluster" is a good feature.
* If we were given new data, we might find which cluster it most resembled but this might be concerning as the data might not look like any of the currently observed data. 
* Also, a clustering offers little in the way of uncertainty assessment, and probably overfits to the small strands of data that we have available.
* for example, all those tiny clusters are unlikely to be "real".

#### See Q5 in Block 3 Portfolio

about interpreting the results DBSCAN.


Extensions: 

* Comment on what we can and cannot do with this representation of the data.
* Would we expect clustering to match labels? How would we do better?

## Mastery:

This worksheet applied dbscan and investigated its' clustering. What about the other clustering methods discussed? Which is best? Which even run? What is the right metric of performance? Can we use labels to choose the correct density threshold?

How would you go about making a real classifier from a clustering model, that can be deployed on a training dataset?