2022-04-30 15:29:07 -05:00

82 lines
1.8 KiB
R

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)
}