R/embed.2007CYT.R

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

#' Directed Graph Embedding by Chen, Yang, and Tang (2007)
#' 
#' @param graph an \code{igraph} object or \eqn{(N\times N)} adjacency matrix.
#' @param ndim an embedding dimension for a given graph.
#' @param alpha perturbation factor \eqn{\in (0,1)}. Default is 0.01
#' 
#' @return a named list containing
#' \describe{
#' \item{embed}{an \eqn{(N\times ndim)} matrix of embedded coordinates.}
#' }
#' 
#' @examples 
#' ## create a simple directed ring graph
#' library(igraph)
#' mygraph = graph.ring(10, directed=TRUE)
#' 
#' ## embed in R^2
#' embed2  = embed.2007CYT(mygraph)$embed
#' 
#' ## visualize the results
#' par(mfrow=c(1,2), pty="s")
#' plot(mygraph, main="directed ring graph")
#' plot(embed2[,1], embed2[,2], pch=19, main="embedded in R^2")
#' 
#' @references 
#' \insertRef{chen_directed_2007}{PNAS}
#' 
#' @author Kisung You
#' @export
embed.2007CYT <- function(graph, ndim=2, alpha=0.01){
  ##############################################################
  # Check Input & Transform into 'igraph' object if matrix
  if (!check_network(graph)){
    stop("* embed.2007CYT : an input graph should be either an 'igraph' object or adjacency matrix.")
  }
  if (is.matrix(graph)){
    if (isSymmetric(graph)){
      graph = igraph::graph_from_adjacency_matrix(graph, mode="undirected")  
    } else {
      graph = igraph::graph_from_adjacency_matrix(graph, mode="directed")  
    }
  }
  W = as.matrix(igraph::as_adjacency_matrix(graph))
  if ((ndim < 1)||(ndim >= nrow(W))){
    stop("* embed.2007CYT : target dimension 'ndim' should be in [1,#{nodes}).")
  }
  ndim = as.integer(ndim)
  if ((alpha <= 0)||(alpha >= 1)){
    stop("* embed.2007CYT : perturbation factor 'alpha' should be a small number in (0,1).")
  }
  nnode = nrow(W) # number of nodes; size of W
  
  ##############################################################
  # Preliminary Computation
  dout = rowSums(W)       # out-degree
  mu   = rep(0,nnode)     # 1 if i-th row of W is all zeros
  ee   = rep(1,nnode)
  for (i in 1:nnode){
    if (all(as.vector(W[i,])==0)){
      mu[i] = 1
    } 
  }
  
  ##############################################################
  # Main Computation
  #   1. teleport random walk
  doutinv = 1/dout
  doutinv[is.infinite(doutinv)] = 0
  P = alpha*(diag(doutinv)%*%W + (1/nnode)*base::outer(mu,ee)) + (1-alpha)*outer(ee,ee)/nnode
  
  #   2. stationary distributions
  vecpi = as.vector(base::Re(base::eigen(t(P))$vectors[,1]))
  vecpi = vecpi/sum(vecpi)
  
  #   3. graph Laplacian
  Phi = diag(vecpi)
  L   = Phi - (Phi%*%P + t(P)%*%Phi)/2
  
  #   4. generalized eigenvalue problem
  geigsol  = aux_geigen(L, Phi, decreasing = FALSE)
  solution = geigsol$vectors[,2:(2+ndim-1)]
  
  ##############################################################
  # Report
  return(list(embed=solution))
}

## I Love Ring Graph
## graph = graph.ring(10, directed=TRUE)