R Under development (unstable) (2024-02-11 r85891 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ####' Tests the causal inference algorithms for interventional data:
> ####' GIES, GES, DP
> ####'
> ####' @author Alain Hauser
> ####' $Id: test_gies.R 454 2017-07-11 15:03:05Z mkalisch $
> 
> cat("Testing the causal inference algorithms for interventional data:\n")
Testing the causal inference algorithms for interventional data:
> 
> library(pcalg)
> 
> source(system.file(package="Matrix", "test-tools-1.R", mustWork=TRUE))
Loading required package: tools
> ##--> showProc.time(), assertError(), relErrV(), ...
> 
> load("test_bicscore.rda") # in directory tests/ i.e., typically *not* installed
> # str(gauss.data)
> p <- ncol(gauss.data)
> 
> (doExtras <- pcalg:::doExtras())
[1] FALSE
> DBG <- if(doExtras) TRUE else FALSE # no debugging by default
> ## Tolerance for numerical comparison
> tol <- sqrt(.Machine$double.eps) # = default for all.equal()
> 
> ## Define all test settings
> settings <- expand.grid(
+     fcn = c("gies", "gds"),
+     cpp = c(FALSE, TRUE),
+     format = c("scatter", "raw"),
+     stringsAsFactors = FALSE)
> nreps <- 5
> 
> for (m in seq_along(settings)) {
+   cat(sprintf("Algorithm: %s, C++: %s, storage format: %s\n", 
+           settings$fcn[m], settings$cpp[m], settings$format[m]))
+   
+   for (i in 1:nreps) {
+     perm <- 1:nrow(gauss.data)
+     
+     ## Randomly permute data
+     if (i > 1) {
+       set.seed(i)
+       perm <- sample(perm)
+     }
+     
+     score <- new("GaussL0penIntScore", 
+                  targets = gauss.targets, 
+                  target.index = gauss.target.index[perm], 
+                  data = gauss.data[perm, ],
+                  format = settings$format[m],
+                  use.cpp = settings$cpp[m])
+     f <- get(settings$fcn[m])
+     est.graph <- f(score, verbose = DBG)
+     for (j in 1:p) {
+       if (!isTRUE(all.equal(est.graph$essgraph$.in.edges[[j]],
+                             gauss.parents[[j]], tolerance = tol)))
+         stop("Parents are not estimated correctly.")
+     }
+       showProc.time()
+   }
+   cat("[Ok]\n")
+ }
Algorithm: gies, C++: FALSE, storage format: scatter
Time (user system elapsed): 0.33 0 0.33 
Time (user system elapsed): 0.06 0 0.06 
Time (user system elapsed): 0.04 0.02 0.07 
Time (user system elapsed): 0.08 0 0.07 
Time (user system elapsed): 0.08 0 0.08 
[Ok]
Algorithm: gds, C++: FALSE, storage format: scatter
Time (user system elapsed): 0.06 0.01 0.08 
Time (user system elapsed): 0.07 0 0.06 
Time (user system elapsed): 0.07 0 0.08 
Time (user system elapsed): 0.08 0 0.08 
Time (user system elapsed): 0.1 0 0.09 
[Ok]
Algorithm: gies, C++: TRUE, storage format: scatter
Time (user system elapsed): 0 0 0 
Time (user system elapsed): 0.01 0 0.02 
Time (user system elapsed): 0 0.02 0.01 
Time (user system elapsed): 0.02 0 0.02 
Time (user system elapsed): 0.01 0 0.02 
[Ok]
> 
> ## Test compatibility with deprecated calling conventions
> cat("Compatibility with deprecated calling conventions... ")
Compatibility with deprecated calling conventions... > score <- new("GaussL0penIntScore", 
+     targets = gauss.targets, 
+     target.index = gauss.target.index, 
+     data = gauss.data)
> 
> warningIssued <- FALSE
> tryCatch(est.graph <- gies(p, gauss.targets, score),
+     warning = function(w) warningIssued <<- TRUE)
> if (!warningIssued) {
+   stop("No warning issued for old calling conventions.")
+ } else {
+   for (j in 1:p) {
+     if (!isTRUE(all.equal(est.graph$essgraph$.in.edges[[j]],
+             gauss.parents[[j]], tolerance = tol)))
+       stop("Parents are not estimated correctly.")
+   }
+ }
> warningIssued <- FALSE
> tryCatch(est.graph <- gies(p = p, targets = gauss.targets, score = score),
+     warning = function(w) warningIssued <<- TRUE)
> if (!warningIssued) {
+   stop("No warning issued for old calling conventions.")
+ }
> cat("[OK]\n")
[OK]
> 
> 
> ## Test stepwise execution of GIES
> cat(if(doExtras)"\n\n", "GIES stepwise", if(doExtras)":\n" else ": ... ",
+     if(doExtras) paste0(paste(rep("=", 14), collapse=""), "\n"),
+     sep = "")
GIES stepwise: ... > for (cpp in c(FALSE, TRUE)) {
+   ## Randomly permute data
+   for (i in 1:nreps) {
+     perm <- 1:nrow(gauss.data)
+     if (i > 1) {
+       set.seed(i)
+       perm <- sample(perm)
+     }
+     score <- new("GaussL0penIntScore",
+         targets = gauss.targets,
+         target.index = gauss.target.index[perm],
+         data = gauss.data[perm, ],
+         use.cpp = cpp)
+ 
+     ## Stepwise execution
+     essgraph <- new("EssGraph", nodes = as.character(1:p), score = score)
+     cont <- TRUE
+     while(cont) {
+       cont <- FALSE
+       while(essgraph$greedy.step("forward")) cont <- TRUE
+       while(essgraph$greedy.step("backward")) cont <- TRUE
+       while(essgraph$greedy.step("backward")) cont <- TRUE
+     }
+     for (i in 1:p) {
+       if(doExtras) cat("  use.cpp = ", cpp,"; i = ", i, "\n", sep="")
+       if (!isTRUE(all.equal(est.graph$essgraph$.in.edges[[i]],
+               gauss.parents[[i]], tolerance = tol)))
+         stop("Parents are not estimated correctly.")
+     }
+     showProc.time()
+   }
+ }
Time (user system elapsed): 0.13 0 0.12 
Time (user system elapsed): 0.08 0 0.08 
Time (user system elapsed): 0.07 0 0.08 
Time (user system elapsed): 0.08 0 0.08 
Time (user system elapsed): 0.1 0 0.09 
Time (user system elapsed): 0 0 0 
Time (user system elapsed): 0.01 0 0.01 
Time (user system elapsed): 0.02 0 0.02 
Time (user system elapsed): 0.01 0 0.02 
Time (user system elapsed): 0.02 0 0.01 
> 
> cat(if(doExtras) "\n", "Done.\n")
 Done.
> 
> ## add test to confirm bug-fix in 2.5-0
> new("GaussParDAG", "a") ## produced a warning prior to 2.5-0
Reference class object of class "GaussParDAG"
Field ".nodes":
[1] "a"
Field ".in.edges":
[[1]]
integer(0)

Field ".params":
[[1]]
[1] 0 0

> 
> proc.time()
   user  system elapsed 
   2.95    0.20    3.15