1# IGraph R package 2# Copyright (C) 2011-2012 Gabor Csardi <csardi.gabor@gmail.com> 3# 334 Harvard street, Cambridge, MA 02139 USA 4# 5# This program is free software; you can redistribute it and/or modify 6# it under the terms of the GNU General Public License as published by 7# the Free Software Foundation; either version 2 of the License, or 8# (at your option) any later version. 9# 10# This program is distributed in the hope that it will be useful, 11# but WITHOUT ANY WARRANTY; without even the implied warranty of 12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13# GNU General Public License for more details. 14# 15# You should have received a copy of the GNU General Public License 16# along with this program; if not, write to the Free Software 17# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 18# 02110-1301 USA 19# 20################################################################### 21 22#' Hierarchical random graphs 23#' 24#' Fitting and sampling hierarchical random graph models. 25#' 26#' A hierarchical random graph is an ensemble of undirected graphs with \eqn{n} 27#' vertices. It is defined via a binary tree with \eqn{n} leaf and \eqn{n-1} 28#' internal vertices, where the internal vertices are labeled with 29#' probabilities. The probability that two vertices are connected in the 30#' random graph is given by the probability label at their closest common 31#' ancestor. 32#' 33#' Please see references below for more about hierarchical random graphs. 34#' 35#' igraph contains functions for fitting HRG models to a given network 36#' (\code{fit_hrg}, for generating networks from a given HRG ensemble 37#' (\code{sample_hrg}), converting an igraph graph to a HRG and back 38#' (\code{hrg}, \code{hrg_tree}), for calculating a consensus tree from a set 39#' of sampled HRGs (\code{consensus_tree}) and for predicting missing edges in 40#' a network based on its HRG models (\code{predict_edges}). 41#' 42#' The igraph HRG implementation is heavily based on the code published by 43#' Aaron Clauset, at his website (not functional any more). 44#' 45#' @name hrg-methods 46#' @family hierarchical random graph functions 47NULL 48 49#' Fit a hierarchical random graph model 50#' 51#' \code{fit_hrg} fits a HRG to a given graph. It takes the specified 52#' \code{steps} number of MCMC steps to perform the fitting, or a convergence 53#' criteria if the specified number of steps is zero. \code{fit_hrg} can start 54#' from a given HRG, if this is given in the \code{hrg} argument and the 55#' \code{start} argument is \code{TRUE}. 56#' 57#' @aliases hrg.fit 58#' @param graph The graph to fit the model to. Edge directions are ignored in 59#' directed graphs. 60#' @param hrg A hierarchical random graph model, in the form of an 61#' \code{igraphHRG} object. \code{fit_hrg} allows this to be \code{NULL}, in 62#' which case a random starting point is used for the fitting. 63#' @param start Logical, whether to start the fitting/sampling from the 64#' supplied \code{igraphHRG} object, or from a random starting point. 65#' @param steps The number of MCMC steps to make. If this is zero, then the 66#' MCMC procedure is performed until convergence. 67#' @return \code{fit_hrg} returns an \code{igraphHRG} object. This is a list 68#' with the following members: 69#' \item{left}{Vector that contains the left children of the internal 70#' tree vertices. The first vertex is always the root vertex, so the 71#' first element of the vector is the left child of the root 72#' vertex. Internal vertices are denoted with negative numbers, starting 73#' from -1 and going down, i.e. the root vertex is -1. Leaf vertices 74#' are denoted by non-negative number, starting from zero and up.} 75#' \item{right}{Vector that contains the right children of the vertices, 76#' with the same encoding as the \code{left} vector.} 77#' \item{prob}{The connection probabilities attached to the internal 78#' vertices, the first number belongs to the root vertex (i.e. internal 79#' vertex -1), the second to internal vertex -2, etc.} 80#' \item{edges}{The number of edges in the subtree below the given 81#' internal vertex.} 82#' \item{vertices}{The number of vertices in the subtree below the 83#' given internal vertex, including itself.} 84#' @references A. Clauset, C. Moore, and M.E.J. Newman. Hierarchical structure 85#' and the prediction of missing links in networks. \emph{Nature} 453, 98--101 86#' (2008); 87#' 88#' A. Clauset, C. Moore, and M.E.J. Newman. Structural Inference of Hierarchies 89#' in Networks. In E. M. Airoldi et al. (Eds.): ICML 2006 Ws, \emph{Lecture 90#' Notes in Computer Science} 4503, 1--13. Springer-Verlag, Berlin Heidelberg 91#' (2007). 92#' @examples 93#' ## We are not running these examples any more, because they 94#' ## take a long time (~15 seconds) to run and this is against the CRAN 95#' ## repository policy. Copy and paste them by hand to your R prompt if 96#' ## you want to run them. 97#' 98#' \dontrun{ 99#' ## A graph with two dense groups 100#' g <- sample_gnp(10, p=1/2) + sample_gnp(10, p=1/2) 101#' hrg <- fit_hrg(g) 102#' hrg 103#' 104#' ## The consensus tree for it 105#' consensus_tree(g, hrg=hrg, start=TRUE) 106#' 107#' ## Prediction of missing edges 108#' g2 <- make_full_graph(4) + (make_full_graph(4) - path(1,2)) 109#' predict_edges(g2) 110#' } 111#' @export 112#' @family hierarchical random graph functions 113 114fit_hrg <- function(graph, hrg=NULL, start=FALSE, steps=0) { 115 # Argument checks 116 if (!is_igraph(graph)) { stop("Not a graph object") } 117 if (is.null(hrg)) { 118 hrg <- list(left=c(), right=c(), prob=c(), edges=c(), 119 vertices=c()) 120 } 121 hrg <- lapply(hrg[c("left","right","prob","edges","vertices")], 122 as.numeric) 123 start <- as.logical(start) 124 steps <- as.integer(steps) 125 126 on.exit( .Call(C_R_igraph_finalizer) ) 127 # Function call 128 res <- .Call(C_R_igraph_hrg_fit, graph, hrg, start, steps) 129 130 if (igraph_opt("add.vertex.names") && is_named(graph)) { 131 res$names <- V(graph)$name 132 } 133 134 class(res) <- "igraphHRG" 135 res 136} 137 138 139#' Create a consensus tree from several hierarchical random graph models 140#' 141#' \code{consensus_tree} creates a consensus tree from several fitted 142#' hierarchical random graph models, using phylogeny methods. If the \code{hrg} 143#' argument is given and \code{start} is set to \code{TRUE}, then it starts 144#' sampling from the given HRG. Otherwise it optimizes the HRG log-likelihood 145#' first, and then samples starting from the optimum. 146#' 147#' @aliases hrg.consensus 148#' @param graph The graph the models were fitted to. 149#' @param hrg A hierarchical random graph model, in the form of an 150#' \code{igraphHRG} object. \code{consensus_tree} allows this to be 151#' \code{NULL} as well, then a HRG is fitted to the graph first, from a 152#' random starting point. 153#' @param start Logical, whether to start the fitting/sampling from the 154#' supplied \code{igraphHRG} object, or from a random starting point. 155#' @param num.samples Number of samples to use for consensus generation or 156#' missing edge prediction. 157#' @return \code{consensus_tree} returns a list of two objects. The first 158#' is an \code{igraphHRGConsensus} object, the second is an 159#' \code{igraphHRG} object. The \code{igraphHRGConsensus} object has the 160#' following members: 161#' \item{parents}{For each vertex, the id of its parent vertex is stored, 162#' or zero, if the vertex is the root vertex in the tree. The first n 163#' vertex ids (from 0) refer to the original vertices of the graph, the 164#' other ids refer to vertex groups.} 165#' \item{weights}{Numeric vector, counts the number of times a given tree 166#' split occurred in the generated network samples, for each internal 167#' vertices. The order is the same as in the \code{parents} vector.} 168 169#' @include auto.R 170#' @family hierarchical random graph functions 171#' @export 172 173consensus_tree <- consensus_tree 174 175 176#' Create a hierarchical random graph from an igraph graph 177#' 178#' \code{hrg} creates a HRG from an igraph graph. The igraph graph must be 179#' a directed binary tree, with \eqn{n-1} internal and \eqn{n} leaf 180#' vertices. The \code{prob} argument contains the HRG probability labels 181#' for each vertex; these are ignored for leaf vertices. 182#' 183#' @aliases hrg.create 184#' @param graph The igraph graph to create the HRG from. 185#' @param prob A vector of probabilities, one for each vertex, in the order of 186#' vertex ids. 187#' @return \code{hrg} returns an \code{igraphHRG} object. 188#' 189#' @family hierarchical random graph functions 190#' @export 191 192hrg <- hrg 193 194 195#' Create an igraph graph from a hierarchical random graph model 196#' 197#' \code{hrg_tree} creates the corresponsing igraph tree of a hierarchical 198#' random graph model. 199#' 200#' @param hrg A hierarchical random graph model. 201#' @return An igraph graph. 202#' 203#' @family hierarchical random graph functions 204#' @export 205 206hrg_tree <- hrg_tree 207 208 209#' Sample from a hierarchical random graph model 210#' 211#' \code{sample_hrg} samples a graph from a given hierarchical random graph 212#' model. 213#' 214#' @aliases hrg.game 215#' @param hrg A hierarchical random graph model. 216#' @return An igraph graph. 217#' 218#' @family hierarchical random graph functions 219#' @export 220 221sample_hrg <- sample_hrg 222 223#' Predict edges based on a hierarchical random graph model 224#' 225#' \code{predict_edges} uses a hierarchical random graph model to predict 226#' missing edges from a network. This is done by sampling hierarchical models 227#' around the optimum model, proportionally to their likelihood. The MCMC 228#' sampling is stated from \code{hrg}, if it is given and the \code{start} 229#' argument is set to \code{TRUE}. Otherwise a HRG is fitted to the graph 230#' first. 231#' 232#' @aliases hrg.predict 233#' @param graph The graph to fit the model to. Edge directions are ignored in 234#' directed graphs. 235#' @param hrg A hierarchical random graph model, in the form of an 236#' \code{igraphHRG} object. \code{predict_edges}s allow this to be 237#' \code{NULL} as well, then a HRG is fitted to the graph first, from a 238#' random starting point. 239#' @param start Logical, whether to start the fitting/sampling from the 240#' supplied \code{igraphHRG} object, or from a random starting point. 241#' @param num.samples Number of samples to use for consensus generation or 242#' missing edge prediction. 243#' @param num.bins Number of bins for the edge probabilities. Give a higher 244#' number for a more accurate prediction. 245#' @return A list with entries: 246#' \item{edges}{The predicted edges, in a two-column matrix of vertex 247#' ids.} 248#' \item{prob}{Probabilities of these edges, according to the fitted 249#' model.} 250#' \item{hrg}{The (supplied or fitted) hierarchical random graph model.} 251#' 252#' @references A. Clauset, C. Moore, and M.E.J. Newman. Hierarchical structure 253#' and the prediction of missing links in networks. \emph{Nature} 453, 98--101 254#' (2008); 255#' 256#' A. Clauset, C. Moore, and M.E.J. Newman. Structural Inference of Hierarchies 257#' in Networks. In E. M. Airoldi et al. (Eds.): ICML 2006 Ws, \emph{Lecture 258#' Notes in Computer Science} 4503, 1--13. Springer-Verlag, Berlin Heidelberg 259#' (2007). 260#' @examples 261#' ## We are not running these examples any more, because they 262#' ## take a long time (~15 seconds) to run and this is against the CRAN 263#' ## repository policy. Copy and paste them by hand to your R prompt if 264#' ## you want to run them. 265#' 266#' \dontrun{ 267#' ## A graph with two dense groups 268#' g <- sample_gnp(10, p=1/2) + sample_gnp(10, p=1/2) 269#' hrg <- fit_hrg(g) 270#' hrg 271#' 272#' ## The consensus tree for it 273#' consensus_tree(g, hrg=hrg, start=TRUE) 274#' 275#' ## Prediction of missing edges 276#' g2 <- make_full_graph(4) + (make_full_graph(4) - path(1,2)) 277#' predict_edges(g2) 278#' } 279#' @export 280#' @family hierarchical random graph functions 281 282predict_edges <- function(graph, hrg=NULL, start=FALSE, num.samples=10000, 283 num.bins=25) { 284 285 # Argument checks 286 if (!is_igraph(graph)) { stop("Not a graph object") } 287 if (is.null(hrg)) { 288 hrg <- list(left=c(), right=c(), prob=c(), edges=c(), 289 vertices=c()) 290 } 291 hrg <- lapply(hrg[c("left","right","prob","edges","vertices")], 292 as.numeric) 293 start <- as.logical(start) 294 num.samples <- as.integer(num.samples) 295 num.bins <- as.integer(num.bins) 296 297 on.exit( .Call(C_R_igraph_finalizer) ) 298 # Function call 299 res <- .Call(C_R_igraph_hrg_predict, graph, hrg, start, num.samples, 300 num.bins) 301 res$edges <- matrix(res$edges, ncol=2, byrow=TRUE) 302 class(res$hrg) <- "igraphHRG" 303 res 304} 305 306 307 308#' Conversion to igraph 309#' 310#' These functions convert various objects to igraph graphs. 311#' 312#' You can use \code{as.igraph} to convert various objects to igraph graphs. 313#' Right now the following objects are supported: \itemize{ \item codeigraphHRG 314#' These objects are created by the \code{\link{fit_hrg}} and 315#' \code{\link{consensus_tree}} functions. } 316#' 317#' @aliases as.igraph as.igraph.igraphHRG 318#' @param x The object to convert. 319#' @param \dots Additional arguments. None currently. 320#' @return All these functions return an igraph graph. 321#' @export 322#' @author Gabor Csardi \email{csardi.gabor@@gmail.com}. 323#' @keywords graphs 324#' @examples 325#' 326#' g <- make_full_graph(5) + make_full_graph(5) 327#' hrg <- fit_hrg(g) 328#' as.igraph(hrg) 329#' 330as.igraph <- function(x, ...) 331 UseMethod("as.igraph") 332 333#' @method as.igraph igraphHRG 334#' @export 335 336as.igraph.igraphHRG <- function(x, ...) { 337 ovc <- length(x$left)+1L 338 ivc <- ovc-1L 339 ll <- ifelse(x$left < 0, -x$left + ovc, x$left + 1) 340 rr <- ifelse(x$right < 0, -x$right + ovc, x$right + 1) 341 edges <- c(rbind(seq_len(ivc)+ovc, ll), rbind(seq_len(ivc)+ovc, rr)) 342 res <- graph(edges) 343 344 V(res)$name <- c(if (!is.null(x$names)) x$names else as.character(1:ovc), 345 paste0("g", 1:ivc)) 346 V(res)$prob <- c(rep(NA, ovc), x$prob) 347 res$name <- "Fitted HRG" 348 res 349} 350 351buildMerges <- function(object) { 352 353 ## Build a merge matrix. This is done by a post-order 354 ## traversal of the tree. 355 356 S <- numeric() 357 vcount <- length(object$left)+1 358 nMerge <- vcount-1 359 merges <- matrix(0, nrow=vcount-1, ncol=3) 360 mptr <- 1 361 S[length(S)+1] <- -1 362 prev <- NULL 363 while (length(S) != 0) { 364 curr <- S[length(S)] 365 ## coming from parent? going left if possible. 366 if (is.null(prev) || 367 (prev < 0 && object$left[-prev] == curr) || 368 (prev < 0 && object$right[-prev] == curr)) { 369 if (curr < 0) { S <- c(S, object$left[-curr]) } 370 ## coming from left child? going right 371 } else if (curr < 0 && object$left[-curr] == prev) { 372 S <- c(S, object$right[-curr]) 373 ## coming from right child? going up 374 } else { 375 if (curr < 0) { 376 merges[mptr,] <- c(object$left[-curr], object$right[-curr], curr) 377 mptr <- mptr + 1 378 } 379 S <- S[-length(S)] 380 } 381 prev <- curr 382 } 383 merges 384} 385 386#' @method as.dendrogram igraphHRG 387 388as.dendrogram.igraphHRG <- function(object, hang=0.01, ...) { 389 390 nMerge <- length(object$left) 391 merges <- buildMerges(object) 392 393 .memberDend <- function(x) { 394 r <- attr(x,"x.member") 395 if(is.null(r)) { 396 r <- attr(x,"members") 397 if(is.null(r)) r <- 1:1 398 } 399 r 400 } 401 402 oHgt <- 1:nrow(merges) 403 hMax <- oHgt[length(oHgt)] 404 mynames <- if (is.null(object$names)) 1:(nMerge+1) else object$names 405 z <- list() 406 407 for (k in 1:nMerge) { 408 x <- merges[k,1:2] 409 if (any(neg <- x >= 0)) { 410 h0 <- if (hang < 0) 0 else max(0, oHgt[k] - hang * hMax) 411 } 412 if (all(neg)) { # two leaves 413 zk <- as.list(x+1) 414 attr(zk, "members") <- 2L 415 attr(zk, "midpoint") <- 1/2 # mean( c(0,1) ) 416 objlabels <- mynames[x+1] 417 attr(zk[[1]], "label") <- objlabels[1] 418 attr(zk[[2]], "label") <- objlabels[2] 419 attr(zk[[1]], "members") <- attr(zk[[2]], "members") <- 1L 420 attr(zk[[1]], "height") <- attr(zk[[2]], "height") <- h0 421 attr(zk[[1]], "leaf") <- attr(zk[[2]], "leaf") <- TRUE 422 } else if (any(neg)) { # one leaf, one node 423 X <- paste0("g", -x) 424 isL <- x[1] >= 0 425 zk <- if (isL) list(x[1]+1, z[[X[2]]]) else list(z[[X[1]]], x[2]+1) 426 attr(zk, "members") <- attr(z[[X[1+isL]]], "members") + 1L 427 attr(zk, "midpoint") <- 428 (.memberDend(zk[[1]]) + attr(z[[X[1+isL]]], "midpoint"))/2 429 attr(zk[[2 - isL]], "members") <- 1L 430 attr(zk[[2 - isL]], "height") <- h0 431 attr(zk[[2 - isL]], "label") <- mynames[x[2 - isL]+1] 432 attr(zk[[2 - isL]], "leaf") <- TRUE 433 } else { #two nodes 434 X <- paste0("g", -x) 435 zk <- list(z[[X[1]]], z[[X[2]]]) 436 attr(zk, "members") <- attr(z[[X[1]]], "members") + 437 attr(z[[X[2]]], "members") 438 attr(zk, "midpoint") <- (attr(z[[X[1]]], "members") + 439 attr(z[[X[1]]], "midpoint") + 440 attr(z[[X[2]]], "midpoint"))/2 441 } 442 attr(zk, "height") <- oHgt[k] 443 z[[k <- paste0("g", -merges[k,3])]] <- zk 444 } 445 z <- z[[k]] 446 class(z) <- "dendrogram" 447 z 448} 449 450#' @importFrom stats as.hclust 451#' @method as.hclust igraphHRG 452 453as.hclust.igraphHRG <- function(x, ...) { 454 merge3 <- buildMerges(x) 455 456 ## We need to rewrite the merge matrix, because hclust assumes 457 ## that group ids are assigned in the order of the merges 458 map <- order(-merge3[,3]) 459 460 merge <- merge3[,1:2] 461 gs <- which(merge < 0) 462 merge[ gs] <- map[ -merge[gs] ] 463 merge[-gs] <- -merge[-gs]-1 464 465 ## To get the ordering, we need to recode the merge matrix again, 466 ## without using group ids. Here the right node is merged _into_ 467 ## the left node. 468 map2 <- numeric(nrow(merge)) 469 mergeInto <- merge 470 for (i in 1:nrow(merge)) { 471 mr <- mergeInto[i,] 472 mr[mr > 0] <- -map2[mr[mr>0]] 473 mergeInto[i,] <- -mr 474 map2[i] <- -mr[1] 475 } 476 n <- nrow(merge)+1 477 hcass <- .C("igraphhcass2", n=as.integer(n), 478 ia=as.integer(mergeInto[,1]), 479 ib=as.integer(mergeInto[,2]), 480 order=integer(n), iia=integer(n), iib=integer(n)) 481 482 mynames <- if (is.null(x$names)) 1:n else x$names 483 res <- list(merge=merge, height=1:nrow(merge), order=hcass$order, 484 labels=mynames, method=NA_character_, 485 dist.method=NA_character_) 486 class(res) <- "hclust" 487 res 488} 489 490#' @method as_phylo igraphHRG 491#' @importFrom stats reorder 492 493as_phylo.igraphHRG <- function(x, ...) { 494 495 ovc <- length(x$left)+1L 496 ivc <- ovc-1L 497 ll <- ifelse(x$left < 0, -x$left + ovc, x$left + 1) 498 rr <- ifelse(x$right < 0, -x$right + ovc, x$right + 1) 499 edge <- matrix(rbind(seq_len(ivc)+ovc, ll, seq_len(ivc)+ovc, rr), 500 ncol=2, byrow=TRUE) 501 502 edge.length <- rep(0.5, nrow(edge)) 503 labels <- if (is.null(x$names)) 1:ovc else x$names 504 obj <- list(edge=edge, edge.length=edge.length/2, tip.label=labels, 505 Nnode=ivc) 506 class(obj) <- "phylo" 507 reorder(obj) 508} 509 510 511 512#' HRG dendrogram plot 513#' 514#' Plot a hierarchical random graph as a dendrogram. 515#' 516#' \code{plot_dendrogram} supports three different plotting functions, selected via 517#' the \code{mode} argument. By default the plotting function is taken from the 518#' \code{dend.plot.type} igraph option, and it has for possible values: 519#' \itemize{ \item \code{auto} Choose automatically between the plotting 520#' functions. As \code{plot.phylo} is the most sophisticated, that is choosen, 521#' whenever the \code{ape} package is available. Otherwise \code{plot.hclust} 522#' is used. \item \code{phylo} Use \code{plot.phylo} from the \code{ape} 523#' package. \item \code{hclust} Use \code{plot.hclust} from the \code{stats} 524#' package. \item \code{dendrogram} Use \code{plot.dendrogram} from the 525#' \code{stats} package. } 526#' 527#' The different plotting functions take different sets of arguments. When 528#' using \code{plot.phylo} (\code{mode="phylo"}), we have the following syntax: 529#' \preformatted{ 530#' plot_dendrogram(x, mode="phylo", colbar = rainbow(11, start=0.7, 531#' end=0.1), edge.color = NULL, use.edge.length = FALSE, \dots) 532#' } The extra arguments not documented above: \itemize{ 533#' \item \code{colbar} Color bar for the edges. 534#' \item \code{edge.color} Edge colors. If \code{NULL}, then the 535#' \code{colbar} argument is used. 536#' \item \code{use.edge.length} Passed to \code{plot.phylo}. 537#' \item \code{dots} Attitional arguments to pass to \code{plot.phylo}. 538#' } 539#' 540#' The syntax for \code{plot.hclust} (\code{mode="hclust"}): \preformatted{ 541#' plot_dendrogram(x, mode="hclust", rect = 0, colbar = rainbow(rect), 542#' hang = 0.01, ann = FALSE, main = "", sub = "", xlab = "", 543#' ylab = "", \dots) 544#' } The extra arguments not documented above: \itemize{ 545#' \item \code{rect} A numeric scalar, the number of groups to mark on 546#' the dendrogram. The dendrogram is cut into exactly \code{rect} 547#' groups and they are marked via the \code{rect.hclust} command. Set 548#' this to zero if you don't want to mark any groups. 549#' \item \code{colbar} The colors of the rectangles that mark the 550#' vertex groups via the \code{rect} argument. 551#' \item \code{hang} Where to put the leaf nodes, this corresponds to the 552#' \code{hang} argument of \code{plot.hclust}. 553#' \item \code{ann} Whether to annotate the plot, the \code{ann} argument 554#' of \code{plot.hclust}. 555#' \item \code{main} The main title of the plot, the \code{main} argument 556#' of \code{plot.hclust}. 557#' \item \code{sub} The sub-title of the plot, the \code{sub} argument of 558#' \code{plot.hclust}. 559#' \item \code{xlab} The label on the horizontal axis, passed to 560#' \code{plot.hclust}. 561#' \item \code{ylab} The label on the vertical axis, passed to 562#' \code{plot.hclust}. 563#' \item \code{dots} Attitional arguments to pass to \code{plot.hclust}. 564#' } 565#' 566#' The syntax for \code{plot.dendrogram} (\code{mode="dendrogram"}): 567#' \preformatted{ 568#' plot_dendrogram(x, \dots) 569#' } The extra arguments are simply passed to \code{as.dendrogram}. 570#' 571#' @aliases hrg.dendrogram 572#' @param x An \code{igraphHRG}, a hierarchical random graph, as returned by 573#' the \code{\link{fit_hrg}} function. 574#' @param mode Which dendrogram plotting function to use. See details below. 575#' @param \dots Additional arguments to supply to the dendrogram plotting 576#' function. 577#' @return Returns whatever the return value was from the plotting function, 578#' \code{plot.phylo}, \code{plot.dendrogram} or \code{plot.hclust}. 579#' @method plot_dendrogram igraphHRG 580#' @export 581#' @author Gabor Csardi \email{csardi.gabor@@gmail.com} 582#' @keywords graphs 583#' @examples 584#' 585#' g <- make_full_graph(5) + make_full_graph(5) 586#' hrg <- fit_hrg(g) 587#' plot_dendrogram(hrg) 588#' 589plot_dendrogram.igraphHRG <- function(x, mode=igraph_opt("dend.plot.type"), ...) { 590 591 if (mode=="auto") { 592 have_ape <- requireNamespace("ape", quietly = TRUE) 593 mode <- if (have_ape) "phylo" else "hclust" 594 } 595 596 if (mode=="hclust") { 597 hrgPlotHclust(x, ...) 598 } else if (mode=="dendrogram") { 599 hrgPlotDendrogram(x, ...) 600 } else if (mode=="phylo") { 601 hrgPlotPhylo(x, ...) 602 } 603} 604 605#' @importFrom graphics plot 606#' @importFrom grDevices rainbow 607#' @importFrom stats rect.hclust 608 609hrgPlotHclust <- function(x, rect=0, colbar=rainbow(rect), hang=.01, 610 ann=FALSE, main="", sub="", xlab="", ylab="", 611 ...) { 612 hc <- as.hclust(x) 613 ret <- plot(hc, hang=hang, ann=ann, main=main, sub=sub, xlab=xlab, 614 ylab=ylab, ...) 615 if (rect > 0) { 616 rect.hclust(hc, k=rect, border=colbar) 617 } 618 invisible(ret) 619} 620 621#' @importFrom graphics plot 622 623hrgPlotDendrogram <- function(x, ...) { 624 plot(as.dendrogram(x), ...) 625} 626 627#' @importFrom graphics plot 628#' @importFrom grDevices rainbow 629 630hrgPlotPhylo <- function(x, colbar=rainbow(11, start=.7, end=.1), 631 edge.color=NULL, use.edge.length=FALSE, ...) { 632 vc <- length(x$left)+1 633 phy <- as_phylo(x) 634 br <- seq(0,1,length=length(colbar)) ; br[1] <- -1 635 cc <- as.integer(cut(x$prob[phy$edge[,1] - vc], breaks=br)) 636 if (is.null(edge.color)) { edge.color <- colbar[cc] } 637 plot(phy, edge.color=edge.color, use.edge.length=use.edge.length, ...) 638} 639 640#' Print a hierarchical random graph model to the screen 641#' 642#' \code{igraphHRG} objects can be printed to the screen in two forms: as 643#' a tree or as a list, depending on the \code{type} argument of the 644#' print function. By default the \code{auto} type is used, which selects 645#' \code{tree} for small graphs and \code{simple} (=list) for bigger 646#' ones. The \code{tree} format looks like 647#' this: \preformatted{Hierarchical random graph, at level 3: 648#' g1 p= 0 649#' '- g15 p=0.33 1 650#' '- g13 p=0.88 6 3 9 4 2 10 7 5 8 651#' '- g8 p= 0.5 652#' '- g16 p= 0.2 20 14 17 19 11 15 16 13 653#' '- g5 p= 0 12 18 } 654#' This is a graph with 20 vertices, and the 655#' top three levels of the fitted hierarchical random graph are 656#' printed. The root node of the HRG is always vertex group #1 657#' (\sQuote{\code{g1}} in the the printout). Vertex pairs in the left 658#' subtree of \code{g1} connect to vertices in the right subtree with 659#' probability zero, according to the fitted model. \code{g1} has two 660#' subgroups, \code{g15} and \code{g8}. \code{g15} has a subgroup of a 661#' single vertex (vertex 1), and another larger subgroup that contains 662#' vertices 6, 3, etc. on lower levels, etc. 663#' The \code{plain} printing is simpler and faster to produce, but less 664#' visual: \preformatted{Hierarchical random graph: 665#' g1 p=0.0 -> g12 g10 g2 p=1.0 -> 7 10 g3 p=1.0 -> g18 14 666#' g4 p=1.0 -> g17 15 g5 p=0.4 -> g15 17 g6 p=0.0 -> 1 4 667#' g7 p=1.0 -> 11 16 g8 p=0.1 -> g9 3 g9 p=0.3 -> g11 g16 668#' g10 p=0.2 -> g4 g5 g11 p=1.0 -> g6 5 g12 p=0.8 -> g8 8 669#' g13 p=0.0 -> g14 9 g14 p=1.0 -> 2 6 g15 p=0.2 -> g19 18 670#' g16 p=1.0 -> g13 g2 g17 p=0.5 -> g7 13 g18 p=1.0 -> 12 19 671#' g19 p=0.7 -> g3 20} 672#' It lists the two subgroups of each internal node, in 673#' as many columns as the screen width allows. 674#' 675#' @param x \code{igraphHRG} object to print. 676#' @param type How to print the dendrogram, see details below. 677#' @param level The number of top levels to print from the dendrogram. 678#' @param ... Additional arguments, not used currently. 679#' @return The hierarchical random graph model itself, invisibly. 680#' 681#' @method print igraphHRG 682#' @export 683#' @family hierarchical random graph functions 684 685print.igraphHRG <- function(x, type=c("auto", "tree", "plain"), 686 level=3, ...) { 687 688 type <- igraph.match.arg(type) 689 if (type=="auto") { 690 type <- if (length(x$left <= 100)) "tree" else "plain" 691 } 692 if (type=="tree") { 693 return(print1.igraphHRG(x, level=level, ...)) 694 } else { 695 return(print2.igraphHRG(x, ...)) 696 } 697} 698 699print1.igraphHRG <- function(x, level=3, ...) { 700 cat(sep="", "Hierarchical random graph, at level ", level, ":\n") 701 702 ## Depth of printed top of the dendrogram 703 .depth <- function(b, l) { 704 l[2] <- max(l[2], nchar(format(x$prob[b], digits=2))) 705 if (l[1]==level) { return(l) } 706 if (x$left[b] < 0 && x$right[b] < 0) { 707 l1 <- .depth(-x$left[b], c(l[1]+1, l[2])) 708 l2 <- .depth(-x$right[b], c(l[1]+1, l[2])) 709 return(pmax(l1,l2)) 710 } 711 if (x$left[b] < 0) { return(.depth(-x$left[b], c(l[1]+1, l[2]))) } 712 if (x$right[b] < 0) { return(.depth(-x$right[b], c(l[1]+1, l[2]))) } 713 return(l) 714 } 715 cs <- .depth(1, c(1, 0)) 716 pw <- cs[2] 717 cs <- cs[1] * 3 718 vw <- nchar(as.character(length(x$left)+1)) 719 sp <- paste(collapse="", rep(" ", cs+pw+2+2)) 720 nn <- if (is.null(x$names)) seq_len(length(x$left)+1) else x$names 721 722 ## Function to collect all individual vertex children 723 .children <- function(b) { 724 res <- c() 725 if (x$left[b] < 0) { 726 res <- c(res, .children(-x$left[b])) 727 } else { 728 res <- c(x$left[b]+1, res) 729 } 730 if (x$right[b] < 0) { 731 res <- c(res, .children(-x$right[b])) 732 } else { 733 res <- c(x$right[b]+1, res) 734 } 735 return(res) 736 } 737 738 ## Recursive printing 739 .plot <- function(b, l, ind = "") { 740 if (b != 1) { 741 he <- format(paste(sep="", ind, "'- g", b), width=cs) 742 ind <- paste(" ", ind) 743 } else { 744 he <- format(paste(sep="", ind, "g", b), width=cs) 745 } 746 ## whether to go left and/or right 747 gol <- x$left[b] < 0 && l < level 748 gor <- x$right[b] < 0 && l < level 749 750 ## the children to print 751 ch1 <- character() 752 if (!gol && x$left[b] < 0) { 753 ch1 <- c(ch1, paste(sep="", "g", -x$left[b])) 754 } 755 if (!gor && x$right[b] < 0) { 756 ch1 <- c(ch1, paste(sep="", "g", -x$right[b])) 757 } 758 ch2 <- numeric() 759 if (!gol) { 760 if (x$left[b] < 0) { ch2 <- c(ch2, .children(-x$left[b])) } 761 if (x$left[b] >= 0) { ch2 <- c(ch2, x$left[b] + 1) } 762 } 763 if (!gor) { 764 if (x$right[b] < 0) { ch2 <- c(ch2, .children(-x$right[b])) } 765 if (x$right[b] >= 0) { ch2 <- c(ch2, x$right[b] + 1) } 766 } 767 768 ## print this line 769 ch2 <- as.character(nn[ch2]) 770 lf <- gsub(" ", "x", format(ch2, width=vw), fixed=TRUE) 771 lf <- paste(collapse=" ", lf) 772 lf <- strwrap(lf, width=getOption("width") - cs - pw - 3 - 2) 773 lf <- gsub("x", " ", lf, fixed=TRUE) 774 if (length(lf) > 1) { 775 lf <- c(lf[1], paste(sp, lf[-1])) 776 lf <- paste(collapse="\n", lf) 777 } 778 op <- paste(sep="", format(he, width=cs), 779 " p=", format(x$prob[b], digits=2, width=pw, justify="left"), 780 " ", paste(collapse=" ", lf)) 781 cat(op, fill=TRUE) 782 783 ## recursive call 784 if (x$left[b] < 0 && l < level) .plot(-x$left[b], l+1, ind) 785 if (x$right[b] < 0 && l < level) .plot(-x$right[b], l+1, ind) 786 } 787 788 ## Do it 789 if (length(x$left) > 0) .plot(b=1, l=1) 790 791 invisible(x) 792} 793 794print2.igraphHRG <- function(x, ...) { 795 cat("Hierarchical random graph:\n") 796 bw <- ceiling(log10(length(x$left)+1))+1 797 p <- format(x$prob, digits=1) 798 pw <- 4 + max(nchar(p)) 799 nn <- if (is.null(x$names)) seq_len(length(x$left)+1) else x$names 800 op <- sapply(seq_along(x$left), function(i) { 801 lc <- if (x$left[i] < 0) { 802 paste(sep="", "g", -x$left[i]) 803 } else { 804 nn[x$left[i]+1] 805 } 806 rc <- if (x$right[i] < 0) { 807 paste(sep="", "g", -x$right[i]) 808 } else { 809 nn[x$right[i]+1] 810 } 811 paste(sep="", format(paste(sep="", "g", i), width=bw), 812 format(paste(sep="", " p=", p[i]), width=pw), 813 "-> ", lc, " ", rc) 814 }) 815 op <- format(op, justify="left") 816 cat(op, sep=" ", fill=TRUE) 817 invisible(x) 818} 819 820## TODO: print as a tree 821 822#' Print a hierarchical random graph consensus tree to the screen 823#' 824#' Consensus dendrograms (\code{igraphHRGConsensus} objects) are printed 825#' simply by listing the children of each internal node of the 826#' dendrogram: \preformatted{HRG consensus tree: 827#' g1 -> 11 12 13 14 15 16 17 18 19 20 828#' g2 -> 1 2 3 4 5 6 7 8 9 10 829#' g3 -> g1 g2} 830#' The root of the dendrogram is \code{g3} (because it has no incoming 831#' edges), and it has two subgroups, \code{g1} and \code{g2}. 832#' 833#' @param x \code{igraphHRGConsensus} object to print. 834#' @param ... Ignored. 835#' @return The input object, invisibly, to allow method chaining. 836#' 837#' @method print igraphHRGConsensus 838#' @export 839#' @family hierarchical random graph functions 840 841print.igraphHRGConsensus <- function(x, ...) { 842 cat("HRG consensus tree:\n") 843 n <- length(x$parents) - length(x$weights) 844 mn <- if (is.null(x$names)) seq_len(n) else x$names 845 id <- c(mn, paste(sep="", "g", seq_along(x$weights))) 846 ch <- tapply(id, x$parents, c)[-1] # first is zero 847 bw <- nchar(as.character(length(x$weights))) 848 vw <- max(nchar(id)) 849 op <- sapply(seq_along(x$weights), function(i) { 850 mych <- format(ch[[i]], width=vw) 851 if (length(ch[[i]])*(vw+1) + bw + 4 > getOption("width")) { 852 mych <- gsub(" ", "x", mych, fixed=TRUE) 853 mych <- paste(collapse=" ", mych) 854 pref <- paste(collapse="", rep(" ", bw+5)) 855 mych <- strwrap(mych, width=getOption("width") - bw - 4, 856 initial="", prefix=pref) 857 mych <- gsub("x", " ", mych, fixed=TRUE) 858 mych <- paste(collapse="\n", mych) 859 } else { 860 mych <- paste(collapse=" ", mych) 861 } 862 paste(sep="", "g", format(i, width=bw), " -> ", mych) 863 }) 864 if (max(nchar(op)) < (getOption("width")-4)/2) { 865 op <- format(op, justify="left") 866 cat(op, sep=" ", fill=TRUE) 867 } else { 868 cat(op, sep="\n") 869 } 870 871 invisible(x) 872} 873 874" 875## How to print HRGs? 876 877B-1 p=0 878'- B-3 p=1 6 879 '- B-7 p=1 2 880 '- B-5 p=1 1 5 881'- B-6 p=1 7 882 '- B-2 p=1 4 883 '- B-4 p=1 3 8 884 885## The same at levels 1, 2 and 3: 886 887B-1 p=0 B-3 B-6 6 2 1 5 7 4 3 8 888 889B-1 p=0 890'+ B-3 p=1 B-7 6 2 1 5 891'+ B-6 p=1 B-2 7 4 3 8 892 893B-1 p=0 894'- B-3 p=1 6 895 '+ B-7 p=1 B-5 2 1 5 896'- B-6 p=1 7 897 '+ B-2 p=1 B-4 4 3 8 898 899## This can be tedious if the graph is big, as we always have n-1 900## internal nodes, we can restrict ourselves to (say) level 3 by default. 901 902## Another possibility is to order the lines according to the group ids. 903 904B-1 p=0 B-3 B-6 905B-2 p=1 B-4 4 906B-3 p=1 B-7 6 907B-4 p=1 3 8 908B-5 p=1 1 5 909B-6 p=1 B-2 7 910B-7 p=1 B-5 2 911 912" 913