tests/test_randDAG.R

In pcalg: Methods for Graphical Models and Causal Inference

library(pcalg)
suppressWarnings(RNGversion("3.5.0"))
## setwd("/sfs/u/kalischm/research/packages/unifDAGs/")
## source("aux_general.R")
## source("randDAG.R")

if (FALSE) { ## Tests currently don't work
  ### Check all methods: ----------------------------------------------
  
  ## MM hack: extract them from the randDAG() function definition
  body. <- body(randDAG)
  is.switch <- function(P) !is.symbol(P) && identical(as.symbol("switch"), P[[1]])
  switchCall <- body.[vapply(body., is.switch, NA)][[1]]
  stopifnot(identical(as.symbol("switch"), switchCall[[1]]))
  (rDAGmeths <- names(switchCall)[-c(1:2, length(switchCall))])
  rDAGall <- function(n, d, ...)
    sapply(rDAGmeths, function(meth) randDAG(n,d, method=meth, ...),
           simplify=FALSE)
  set.seed(37)
  rD.10.4 <- withCallingHandlers(
    rDAGall(10, 4) # "low-level warning" -- get it here and test it below
    , warning = function(w) {
      rDAG.warn <<- conditionMessage(w); invokeRestart("muffleWarning") })
  ## with a low-level warning:
  ## IGNORE_RDIFF_BEGIN
  rDAG.warn
  ## IGNORE_RDIFF_END
  stopifnot(grepl("graph_molloy_.*Cannot shuffle graph", rDAG.warn))
  ##
  rD.10.4 # looks ok
  ## Show, but ignore the package startup messages:
  ## IGNORE_RDIFF_BEGIN
  stopifnot( require("graph") )
  ## IGNORE_RDIFF_END
  
  stopifnot(vapply(rD.10.4, isDirected, NA),
            vapply(rD.10.4, inherits, NA, what="graph"))
  ## nice plot of all 8 :
  op <- par(mfrow=c(4,2))
  invisible(lapply(names(rD.10.4), function(nm) plot(rD.10.4[[nm]], main=nm)))
  par(op)
  
  
  stopifnot(identical(
    lapply(rD.10.4, leaves, "out"),
    list(er = "3", regular = c("1", "5", "6"), watts = c("3", "4", "6"),
         bipartite = c("1", "2", "5"), barabasi = c("4", "8"),
         geometric = c("4", "7"), power = c("4", "5", "9"),
         interEr = c("3", "7"))
  ))
  
  stopifnot(identical(
    lapply(rD.10.4, leaves, "in"),
    list(er = c("1", "4", "7"), regular = c("3", "7", "10"),
         watts = c("1", "8"), bipartite = c("4", "6"),
         barabasi = c("6", "7"), geometric = c("5", "10"),
         power = c("2", "7"), interEr = c("8", "10"))
  ))
  
  set.seed(47)
  rD.12.2 <- rDAGall(12, 2)
  stopifnot(vapply(rD.12.2, isDirected, NA),
            vapply(rD.12.2, numNodes, 1) == 12,
            identical(vapply(rD.12.2, numEdges, 1),
                      setNames(c(9, 12, 12, 11, 11, 11, 13, 8), rDAGmeths))
  )
  
  ## Use the output here
  require(Matrix)
  lapply(rD.10.4, function(g) as(as(g, "Matrix"),"nMatrix"))
  lapply(rD.12.2, function(g) as(as(g, "Matrix"),"nMatrix"))
}
##---------------------------------------------------------------------------

## Minimal checks on graphs generated via igraph
stopifnot( require("graph") )
set.seed(37)
dagList <- vector(mode = "list", length = 7)
dagList[[1]] <- randDAG(10, 4, "regular")
dagList[[2]] <- randDAG(10, 4, "watts")
dagList[[3]] <- randDAG(10, 4, "er")
dagList[[4]] <- randDAG(10, 4, "bipartite")
dagList[[5]] <- randDAG(10, 4, "barabasi")
dagList[[6]] <- randDAG(10, 4, "geometric")
dagList[[7]] <- randDAG(10, 4, "interEr", par2 = 0.5)

## number of nodes
stopifnot(all.equal(
  sapply(dagList, numNodes),
  rep(10,7)
))

## number of edges
stopifnot(all.equal(
  sapply(dagList, numEdges),
  c(20,20,16,15,0,10,19)
))

## check weights
set.seed(123)
n <- 100
g <- randDAG(n=n,d=3, wFUN=list(runif,min=0,max=1))
m <- wgtMatrix(g)
stopifnot(sum(m != 0) == 137)
v <- as.numeric(m)
v <- v[v!=0]
## dput(as.vector(summary(v, digits=7)))
stopifnot(all.equal(as.vector(summary(v, digits=7)),
                    c(0.008103577, 0.2589966, 0.5287397,
                      0.5232445, 0.8159941, 0.9915566)))
ct <- cut(x=v, breaks=seq(0,1,by=0.1))
stopifnot(all.equal(chisq.test(as.numeric(table(ct)), p = rep(0.1,10))$p.value,
                    0.3101796548))

## check generation of negative weights (fixed Bug)
set.seed(123) ; tmp1 <- randDAG(3,2, wFUN = list(runif, min =  2, max =  2))
set.seed(123) ; tmp2 <- randDAG(3,2, wFUN = list(runif, min = -2, max = -2))
stopifnot(unlist(tmp1@edgeData@data) ==  2,
	  unlist(tmp2@edgeData@data) == -2 )