#' 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.