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