Graph laplacian and reursive newman modularity

Schrick-Noah_CG-Analysis.R

@@ -11,8 +11,24 @@ setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
 source("./CG_Files/manual_import.R")
 
 car <- import_networks(1)
+car.netname <- "Vehicle Maintenance"
 hipaa <- import_networks(2)
+hipaa.netname <- "HIPAA Compliance"
 pci <- import_networks(3)
+pci.netname <- "PCI Compliance"
+
+# Get basic network attributes
+car.adj <- get.adjacency(car)
+car.deg <- rowSums(as.matrix(car.adj)) # degree
+car.n <- length(V(car))
+
+hipaa.adj <- get.adjacency(hipaa)
+hipaa.deg <- rowSums(as.matrix(hipaa.adj)) # degree
+hipaa.n <- length(V(hipaa))
+
+pci.adj <- get.adjacency(pci)
+pci.deg <- rowSums(as.matrix(pci.adj)) # degree
+pci.n <- length(V(pci))
 
 ################################ Centralities ################################
 source("centralities.R")
@@ -29,8 +45,73 @@ pci.pr <- page.rank(pci)
 
 ### K-path
 car.kpe <- geokpath(car, V(car), "out")
+hipaa.kpe <- geokpath(hipaa, V(hipaa), "out")
+pci.kpe <- geokpath(pci, V(pci), "out")
 
 
+############## Clustering
+source("self_newman_mod.R")
+
+### Laplacian
+car.Lap <- diag(car.deg) - car.adj # L = D-A
+hipaa.Lap <- diag(hipaa.deg) - hipaa.adj
+pci.Lap <- diag(pci.deg) - pci.adj
+
+# get eigvals and vecs
+car.eigs <- eigen(car.Lap)$vectors[,car.n-1]
+car.eig_val <- eigen(car.Lap)$values[car.n-1]
+names(car.eigs) <- names(V(car))
+car.l_clusters <- ifelse(car.eigs>0,1,-1)
+
+hipaa.eigs <- eigen(hipaa.Lap)$vectors[,hipaa.n-1]
+hipaa.eig_val <- eigen(hipaa.Lap)$values[hipaa.n-1]
+names(hipaa.eigs) <- names(V(hipaa))
+hipaa.l_clusters <- ifelse(hipaa.eigs>0,1,-1)
+
+pci.eigs <- eigen(pci.Lap)$vectors[,pci.n-1]
+pci.eig_val <- eigen(pci.Lap)$values[pci.n-1]
+names(pci.eigs) <- names(V(pci))
+pci.l_clusters <- ifelse(pci.eigs>0,1,-1)
+
+### Recursive Newmann
+car.modularity <- fastgreedy.community(car,merges=TRUE, modularity=TRUE, membership=TRUE)
+car.membership.ids <- unique(car.modularity$membership)
+cat(paste('Number of detected communities in the car network =',length(car.membership.ids)))
+cat("community sizes: ")
+sapply(membership.ids,function(x) {sum(x==car.modularity$membership)})
+cat("modularity: ")
+max(car.modularity$modularity)
+V(car)$color=car.modularity$membership
+plot(car,vertex.size=10,
+     vertex.label=NA,vertex.color=V(car)$color,
+     main=paste(car.netname, " Recursive Newman Modularity"))
+
+
+hipaa.modularity <- fastgreedy.community(hipaa,merges=TRUE, modularity=TRUE, membership=TRUE)
+hipaa.membership.ids <- unique(hipaa.modularity$membership)
+cat(paste('Number of detected communities in the HIPAA network =',length(hipaa.membership.ids)))
+cat("community sizes: ")
+sapply(membership.ids,function(x) {sum(x==hipaa.modularity$membership)})
+cat("modularity: ")
+max(hipaa.modularity$modularity)
+V(hipaa)$color=hipaa.modularity$membership
+plot(hipaa,vertex.size=10,
+     vertex.label=NA,vertex.color=V(hipaa)$color,
+     main=paste(hipaa.netname, " Recursive Newman Modularity"))
+
+
+pci.modularity <- fastgreedy.community(pci,merges=TRUE, modularity=TRUE, membership=TRUE)
+pci.membership.ids <- unique(pci.modularity$membership)
+cat(paste('Number of detected communities in the PCI network =',length(pci.membership.ids)))
+cat("community sizes: ")
+sapply(membership.ids,function(x) {sum(x==pci.modularity$membership)})
+cat("modularity: ")
+max(pci.modularity$modularity)
+V(pci)$color=pci.modularity$membership
+plot(pci,vertex.size=10,
+     vertex.label=NA,vertex.color=V(pci)$color,
+     main=paste(pci.netname, " Recursive Newman Modularity"))
+
 ############# Other- Tmp work
 min_cut(car,"0", "2490")
 min_cut(hipaa,"0","2320")

self_estrada.R

@@ -0,0 +1,81 @@
+estrada.index <- function(A, beta=NULL){
+    g <- A
+    if (class(A) == 'igraph'){
+        # Error checking. Turn into adj matrix if needed.
+        A <- get.adjacency(A)
+    }
+    
+    
+    if (is.null(beta)){
+        beta <- 1.0
+    }
+    
+    lam.dom <- eigen(A)$values[1] #dom eigenvec
+
+    A.eigs <- eigen(A) 
+    V <- A.eigs$vectors # where columns are the v_i terms
+    lams <- A.eigs$values
+    n <- length(lams)
+    
+    # Create subfunction to compute centrality for one node, then use sapply
+    # for all nodes
+    subg.node.i <- function(i){sum(V[i,]^2*exp(beta*lams))} 
+    subg.all <- sapply(1:n, subg.node.i)
+    EE <- sum(subg.all)
+    
+    return(EE)
+}
+
+microstate.prob <- function(A, beta=NULL){
+    EE <- estrada.index(A, beta)
+    g <- A
+    if (class(A) == 'igraph'){
+        # Error checking. Turn into adj matrix if needed.
+        A <- get.adjacency(A)
+    }
+    
+    if (is.null(beta)){
+        beta <- 1.0
+    }
+    
+    A.eigs <- eigen(A) 
+    lams <- A.eigs$values
+
+    probs <- (exp(beta*lams))/EE
+    
+    # Experiment with lambda being negative
+    #probs <- (exp(beta*-lams))/EE
+    
+    
+    # Add names to output
+    names(probs) <- V(g)$name
+    return(probs)
+}
+
+entropy <- function(A, beta=NULL, kb=NULL){
+    microstate_probs <- microstate.prob(A, beta)
+    EE <- estrada.index(A, beta)
+    g <- A
+    
+    if (class(A) == 'igraph'){
+        # Error checking. Turn into adj matrix if needed.
+        A <- get.adjacency(A)
+    }
+    
+    if (is.null(beta)){
+        beta <- 1.0
+    }
+    
+    if (is.null(kb)){
+        kb <- 1.0
+    }
+    
+    lam.dom <- eigen(A)$values[1] #dom eigenvec
+    A.eigs <- eigen(A) 
+    V <- A.eigs$vectors # where columns are the v_i terms
+    lams <- A.eigs$values
+    
+    S <- -kb*beta*sum(lams*microstate_probs)+kb*log(EE)*sum(microstate_probs)
+    return(S)
+}
+    

self_newman_mod.R

@@ -0,0 +1,43 @@
+newman_mod <- function(g, weights=NULL){
+  A <- get.adjacency(g) # adj
+  m <- ecount(g)
+  n <- vcount(g)
+  
+  if (is.null(weights)){
+    weights <- rep(1,n)
+  }
+  
+  # Obtain the modularity matrix
+  B.node.i <- function(i){degree(g)[i]*degree(g)}
+  B.node.all <- sapply(1:n, B.node.i)
+  B <- A - (B.node.all/(2*m))
+  # NOTE: This is identical to: modularity_matrix(g) ! Can verify with:
+  # modularity_matrix(g) == B
+  B.eigs <- eigen(B)
+  max.lam <- B.eigs$values[1]
+  s <- ifelse(B.eigs$vectors[,1]>0,1,-1)
+  
+  weights <- B.eigs$vectors[n]/B.eigs$vectors[,1]
+  # Plotting
+  #V(g)$color <- ifelse(B[1,]>0,"green","yellow")
+  V(g)$color <- ifelse(B.eigs$vectors[,1]>0,"green","yellow")
+  V(g)$size <- 10
+  plot(g, main=paste(g.netname, " Newman Modularity"))
+  clust1 = list()
+  clust2 = list()
+  clusters = list()
+  
+  # Make list of clusters
+  for(i in 1:n){
+      ifelse(V(g)[i]$color=="green",
+              clust1 <- append(clust1, V(g)[i]$name),
+              clust2 <- append(clust2, V(g)[i]$name))}
+
+  clusters <- list(clust1, clust2)
+  
+  Q.node.i <- function(i){sum(
+    (((B.eigs$vectors[i])*weights[i]*s)^2)*B.eigs[i]$values)}
+  Q <- (1/(4*m))*sapply(1:n, Q.node.i)
+ 
+  return(list(Q=Q,max.lam=max.lam,weights=weights,clusters=clusters))
+}