#' Directed Graph Embedding by Perrault-Joncas and Meila (2011)
#'
#' @param graph an \code{igraph} object or \eqn{(N\times N)} affinity matrix.
#' @param ndim an embedding dimension for a given graph.
#'
#' @return a named list containing
#' \describe{
#' \item{embed}{an \eqn{(N\times ndim)} matrix of embedded coordinates.}
#' \item{field}{an \eqn{(N\times ndim)} matrix of directional information on each one.}
#' }
#'
#' @examples
#' ## create a simple directed ring graph
#' library(igraph)
#' mygraph = graph.ring(10, directed=TRUE)
#'
#' ## embed in R^2
#' solution = embed.2011PJM(mygraph)
#' coords = solution$embed
#' fields = solution$field
#'
#' ## extract coordinates information
#' x0 = coords[,1]; y0 = coords[,2]
#' x1 = x0 + fields[,1]
#' y1 = y0 + fields[,2]
#'
#' ## visualize the results
#' par(mfrow=c(1,2), pty="s")
#' plot(mygraph, main="directed ring graph")
#' plot(x0,y0, pch=19, main="embedded in R^2",xlim=c(-0.5,0.5),ylim=c(-0.5,0.5))
#' text(x0+0.05,y0+0.05,labels=1:10)
#' arrows(x0,y0,x1,y1,length=0.1)
#'
#'
#' @export
embed.2011PJM <- function(graph, ndim=2){
##############################################################
# Check Input & Transform into 'igraph' object if matrix
if (!check_network(graph)){
stop("* embed.2011PJM : an input graph should be either an 'igraph' object or affinity 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 (isSymmetric(W)){
warning("* embed.2011PJM : directional vector field is 0 since graph is undirected.")
}
if ((ndim < 1)||(ndim >= nrow(W))){
stop("* embed.2011PJM : target dimension 'ndim' should be in [1,#{nodes}).")
}
ndim = as.integer(ndim)
nnode = nrow(W) # number of nodes; size of W
##############################################################
# COMPUTATION
# PART 1 : estimate embedding coordinates
S = (W+t(W))/2
Q = diag(rowSums(S)); Qinv = base::solve(Q)
V = Qinv%*%S%*%Qinv
Q1 = diag(rowSums(V))
Hss = base::solve(Q1,V)
eigHss = base::eigen(Hss)
if (nnode < (ndim+1)){
warning("* embed.2011PJM : we don't have enough data points for embedding you want.")
}
idstart = max(which.min((eigHss$values < 1)),2) # for spectral clustering, it's important to exclude non-null ones
if ((ndim+idstart-1)>nnode){
stop("* embed.2011PJM : invalid eigendecomposition. Try smaller 'ndim'.")
}
Phi = eigHss$vectors[,idstart:min((ndim+idstart-1),nnode)]
Lbd = eigHss$values[idstart:min((idstart+ndim+1),nnode)]
# PART 2 : Density
vecpi = as.vector(base::eigen(t(Hss))$vectors[,1])
vecpi = vecpi/sum(vecpi)
# PART 3 : Estimate Vector Field
PP = diag(rowSums(W)); PPinv = solve(PP)
TT = PPinv%*%W%*%PPinv
P1 = diag(rowSums(TT))
Haa = base::solve(P1,TT)
# PART 4 : compute vector field
# R = (Phi%*%diag(Lbd) - Haa%*%Phi)/2 # arXiv
R = (Haa-Hss)%*%Phi/2
##############################################################
# REPORT
output = list()
output$embed = Phi
output$field = R
return(output)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.