Schrick-Noah_CG-Analysis.R

@@ -11,24 +11,8 @@ 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")
@@ -45,73 +29,8 @@ 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

@@ -1,81 +0,0 @@
-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

@@ -1,43 +0,0 @@
-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))
-}