#' Compute Effective Resistance
#'
#' @param graph an \code{igraph} object or \eqn{(N\times N)} adjacency matrix.
#' @param pad.inf a logical; \code{TRUE} to take extremely large values as \code{Inf}, \code{FALSE} otherwise.
#'
#' @return an \eqn{(N\times N)} matrix containing effective resistance values on edge set.
#'
#' @examples
#' ## create a graph with disconnected and singleton components, adapted from
#' # https://stackoverflow.com/questions/29730624/how-to-split-an-igraph-into-connected-subgraphs
#' library(igraph)
#' g <- simplify(
#' graph.compose(graph.ring(10),
#' graph.star(5, mode = "undirected"))
#' ) + edge("7","8")
#'
#' # compute effective resistance
#' ermat = effective(g)
#'
#' # visualize
#' par(mfrow=c(1,2),pty="s")
#' plot(g, main="the graph 'g'")
#' image(ermat[,10:1],col=gray(0:32/32),
#' main="white for Inf values",xlab="",ylab="")
#'
#'
#' @references
#' \insertRef{young_new_2016}{PNAS}
#'
#' \insertRef{young_new_2016-1}{PNAS}
#'
#' @author Kisung You
#' @export
effective <- function(graph, pad.inf=FALSE){
##############################################################
# Check Input & Transform into 'igraph' object if matrix
input = graph
if (!check_network(input)){
stop("* effective : an input graph should be either an 'igraph' object or adjacency matrix.")
}
if (is.matrix(input)){
if (isSymmetric(input)){
input = igraph::graph_from_adjacency_matrix(input, mode="undirected")
} else {
input = igraph::graph_from_adjacency_matrix(input, mode="directed")
}
}
dflag = (igraph::is_directed(input)) # directed ?
cflag = (igraph::is_connected(input)) # connected ?
adjmat = (as.matrix(igraph::as_adjacency_matrix(input))) # adjacency matrix
NN = nrow(adjmat)
##############################################################
# compute and return : let's just ignore and go with directed case
gg = effective_cnn(adjmat)
if (pad.inf){ # pad inf values
gg[(gg>1e+10)] = Inf
} else {
if (!cflag){
warning("* effective : input graph is not connected. Larger values mean those edge sets are between non-connected components.")
}
}
return(gg)
# ## old codes
# cnnd = igraph::components(input)
# if (cnnd$no < 1){ # no connected component
# stop("* effective : there is no connected component in the given network.")
# } else if (cnnd$no < 2){ # case 1 : there is single one
# return(effective_cnn(adjmat))
# } else { # case 2 : multiple connected components
# # prepare basic computation
# output = array(Inf, c(NN,NN))
# # membership
# ulabel = unique(cnnd$membership)
# if (length(ulabel)!=(cnnd$no)){
# stop("* something wrong here..")
# }
# for (k in 1:cnnd$no){
# tgtid = which(cnnd$membership==ulabel[k])
# if (length(tgtid) > 1){ # there might exist singletons.. ignore !
# output[tgtid,tgtid] = effective_cnn(adjmat[tgtid,tgtid])
# }
# }
# # return the output
# return(output)
}
# auxiliary ---------------------------------------------------------------
#' @keywords internal
#' @noRd
effective_cnn <- function(A){ # adjacency matrix
# parameters and pre-compute L
if (inherits(A, "dgCMatrix")){
A = as.matrix(A)
}
N = nrow(A)
diag(A) = rep(0, as.integer(N))
L = base::diag(base::rowSums(A)) - A # graph laplacian
# 1. compute Q
tmp = diag(N)
tmp[,1] = rep(1,N)/sqrt(N)
qrq = base::qr.Q(base::qr(tmp))
Q = t(qrq[,2:N])
# 2. compute \bar{L}
Lbar = Q%*%L%*%t(Q)
# 3. solve Lyapunov equation
Sigma = maotai::lyapunov(Lbar, diag(rep(1,N-1)))
# 4. compute X
X = 2*(t(Q)%*%Sigma%*%Q)
# 5. compute
out.cpp = cpp_effective_arith(X)
if (!isSymmetric(out.cpp)){
out.cpp = (out.cpp + t(out.cpp))/2
}
return(out.cpp)
}
#' @keywords internal
#' @noRd
effectivesym <- function(A){
n = nrow(A)
L = diag(rowSums(A))-A
Linv = aux_pinv(L)
output = array(0,c(n,n))
for (i in 1:(n-1)){
for (j in (i+1):n){
est = rep(0,n)
est[i] = 1
est[j] = -1
output[i,j] <- output[j,i] <-sum(as.vector(Linv%*%est)*est)
}
}
return(output)
}
# # check some conditions for Q matrix
# N = ncol(Q)
# Q%*%rep(1,N) # Q1n = 0
# Q%*%t(Q) # QQt = I
# t(Q)%*%Q # QtQ = Projection
# ## compare sym and naive implementation
# # 1. build a fully connected graph
# karate <- igraph::make_graph("Zachary")
# A <- as.matrix(igraph::as_adjacency_matrix(karate))
# e1 = effective(A)
# e2 = effectivesym(A)
#
# # 2. visualize :: result -> same !
# graphics.off()
# par(mfrow=c(1,2), pty="s")
# image(e1, main="asymmetric")
# image(e2, main="original")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.