From d7e9098582597d27e11e76c001ced062f3483246 Mon Sep 17 00:00:00 2001
From: noah <noah@theschricks.com>
Date: Sat, 30 Apr 2022 15:29:07 -0500
Subject: [PATCH] Graph laplacian

---
 self_estrada.R    | 81 +++++++++++++++++++++++++++++++++++++++++++++++
 self_newman_mod.R | 43 +++++++++++++++++++++++++
 2 files changed, 124 insertions(+)
 create mode 100644 self_estrada.R
 create mode 100644 self_newman_mod.R

diff --git a/self_estrada.R b/self_estrada.R
new file mode 100644
index 0000000..b495b08
--- /dev/null
+++ b/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)
+}
+    
diff --git a/self_newman_mod.R b/self_newman_mod.R
new file mode 100644
index 0000000..26a0005
--- /dev/null
+++ b/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))
+}