R/effective.R

In kyoustat/PNAS: 'PNAS' is for Network Analysis and Statistical methods

#' 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")