R Under development (unstable) (2024-09-06 r87103 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.

> library(pcalg)
> (doExtras <- pcalg:::doExtras())
[1] FALSE
> 
> showProc.time <- local({
+     pct <- proc.time()
+     function() { ## CPU elapsed __since last called__
+ 	ot <- pct ; pct <<- proc.time()
+ 	cat('Time elapsed: ', (pct - ot)[1:3],'\n')
+     }
+ })
> 
> 
> ##################################################
> ## Standard PC
> ##################################################
> ##library(RBGL)
> nreps <- 10
> 
> p <- 10
> set.seed(234)
> for (ii in 1:nreps) {
+   ## generate and draw random DAG :
+   myDAG <- randomDAG(p, prob = 0.3)
+   myCPDAG <- dag2cpdag(myDAG)
+   suffStat <- list(C = cov2cor(trueCov(myDAG)), n = 10^9)
+   res <- pc(suffStat, indepTest=gaussCItest, 0.99, p=p)
+   if( shd(res, myCPDAG) != 0)
+     stop("Test pc wrong: CPDAG ",ii," was not found correctly")
+ }
> 
> showProc.time()
Time elapsed:  0.52 0.05 0.56 
> 
> if (doExtras) {
+ ##################################################
+ ## Conservative PC
+ ##################################################
+ ##PC algorithm sample (compared with Tetrad)
+ ##__________________________________________________________________________________
+ ## Example 1
+ p <- 10
+ n <- 5000
+ set.seed(p*37+15673)
+ g <- randomDAG(p,2/(p-1))
+ 
+ ## generate n samples of DAG using standard normal error distribution
+ ##random data
+ 
+ (load(system.file("external", "test_conservative_pc_data1.rda", package = "pcalg")))
+ ## set.seed(67*37)
+ ## new.mat <- rmvnorm(n,mean=rep(0,p),sigma=trueCov(g))
+ ## save(new.mat, file = "/u/kalisch/research/packages/pcalg/inst/external/test_conservative_pc_data1.rda")
+ ## load(file = "/u/kalisch/research/packages/pcalg/inst/external/test_conservative_pc_data1.rda")
+ 
+ suffStat.data <- list(C=cor(new.mat1),n=n)
+ 
+ ##pcAlgo conservative sample
+ dag1 <- pc(suffStat.data, gaussCItest, alpha=0.005, p=p,
+ 	   u2pd="relaxed", conservative=TRUE)
+ ##adjacency matrix
+ dag1.amat <- as(dag1@graph,"matrix")
+ 
+ ##always save the transpose of the matrix to be loaded in Tetrad
+ ##write(t(new.mat1),file="test_conservative_pc_data1.txt",ncolumns=p)
+ 
+ ##check the output with Tetrad
+ 
+ amat.tetrad1 <- matrix(0,p,p)
+ amat.tetrad1[1,6] <- 1
+ amat.tetrad1[2,4] <- amat.tetrad1[2,6] <- 1
+ amat.tetrad1[3,4] <- 1
+ amat.tetrad1[4,6] <- amat.tetrad1[4,9] <- 1
+ 
+ correctEst1 <- all(dag1.amat == amat.tetrad1)
+ if (!correctEst1) stop("Test sample conservative PC wrong: example 1!")
+ showProc.time()
+ 
+ 
+ ## Example 2
+ p <- 12
+ n <- 5000
+ set.seed(p*67+15673)
+ g <- randomDAG(p,2/(p-1))
+ 
+ ## generate n samples of DAG using standard normal error distribution
+ ##random data
+ load(system.file("external", "test_conservative_pc_data2.rda", package = "pcalg"))
+ ## set.seed(67*67)
+ ## new.mat <- rmvnorm(n,mean=rep(0,p),sigma=trueCov(g))
+ ## save(new.mat, file = "/u/kalisch/research/packages/pcalg/inst/external/test_conservative_pc_data2.rda")
+ ## load(file = "/u/kalisch/research/packages/pcalg/inst/external/test_conservative_pc_data2.rda")
+ 
+ suffStat.data <- list(C=cor(new.mat2),n=n)
+ ## pcAlgo conservative sample
+ dag2 <- pc(suffStat.data, gaussCItest, alpha=0.005, p=p,
+            u2pd="relaxed", conservative=TRUE)
+ 
+ ##adjacency matrix
+ dag2.amat <- as(dag2@graph,"matrix")
+ 
+ ##always save the transpose of the matrix to be loaded in Tetrad
+ ##write(t(new.mat),file="test_conservative_pc_data2.txt",ncolumns=p)
+ 
+ ##check the output with Tetrad
+ 
+ amat.tetrad2 <- matrix(0,p,p)
+ amat.tetrad2[1,2] <- amat.tetrad2[1,7] <- 1
+ amat.tetrad2[2,1] <- amat.tetrad2[2,11] <- 1
+ amat.tetrad2[3,9] <- amat.tetrad2[3,11] <- amat.tetrad2[3,12] <- 1
+ amat.tetrad2[4,5] <- amat.tetrad2[4,11] <- amat.tetrad2[4,12] <- 1
+ amat.tetrad2[5,4] <- amat.tetrad2[5,8] <- amat.tetrad2[5,9] <- 1
+ amat.tetrad2[7,1] <- 1
+ amat.tetrad2[8,5] <- amat.tetrad2[8,12] <- 1
+ amat.tetrad2[10,11] <- 1
+ 
+ correctEst2 <- all(dag2.amat == amat.tetrad2)
+ if (!correctEst2) stop("Test sample conservative PC wrong: example 2!")
+ showProc.time()
+ 
+ 
+ ##PC algorithm population
+ ##_____________________________________________________________________________
+ ## Example 4
+ p <- 15
+ set.seed(15673)
+ g <- randomDAG(p, 2/(p-1))
+ 
+ ##population version
+ suffStat <- list(C=cov2cor(trueCov(g)),n=10^9)
+ dag4 <- pc(suffStat, gaussCItest, alpha=0.9999, p=p, u2pd="relaxed")
+ dag5 <- pc(suffStat, gaussCItest, alpha=0.9999, p=p, u2pd="relaxed",conservative=TRUE)
+ 
+ ##adjacency matrix
+ dag4.amat <- as(dag4@graph,"matrix")
+ dag5.amat <- as(dag5@graph,"matrix")
+ 
+ correctEst4 <- all(dag4.amat == dag5.amat)
+ if (!correctEst4) stop("Test population conservative PC wrong: example 4!")
+ showProc.time()
+ 
+ 
+ ## Example 5
+ p <- 25
+ set.seed(1589873)
+ g <- randomDAG(p,2/(p-1))
+ 
+ ##population version
+ suffStat <- list(C=cov2cor(trueCov(g)),n=10^9)
+ dag6 <- pc(suffStat, gaussCItest, alpha=0.9999, p=p, u2pd="relaxed")
+ dag7 <- pc(suffStat, gaussCItest, alpha=0.9999, p=p, u2pd="relaxed",conservative=TRUE)
+ 
+ ##adjacency matrix
+ dag6.amat <- as(dag6@graph,"matrix")
+ dag7.amat <- as(dag7@graph,"matrix")
+ 
+ correctEst5 <- all(dag6.amat == dag7.amat)
+ if (!correctEst5) stop("Test population conservative PC wrong: example 5!")
+ showProc.time()
+ 
+ ## Example 6
+ p <- 35
+ set.seed(78673)
+ g <- randomDAG(p,2/(p-1))
+ 
+ ##population version
+ suffStat <- list(C=cov2cor(trueCov(g)),n=10^9)
+ dag8 <- pc(suffStat, gaussCItest, alpha=0.9999, p=p, u2pd="relaxed")
+ dag9 <- pc(suffStat, gaussCItest, alpha=0.9999, p=p, u2pd="relaxed", conservative=TRUE)
+ 
+ ##adjacency matrix
+ dag8.amat <- as(dag8@graph,"matrix")
+ dag9.amat <- as(dag9@graph,"matrix")
+ 
+ correctEst6 <- all(dag8.amat == dag9.amat)
+ if (!correctEst6) stop("Test population conservative PC wrong: example 6!")
+ 
+ showProc.time()
+ }
> 
> ##################################################
> ## Test applyOrientationRules
> ##################################################
> amat <- matrix(c(0, 0, 1, 1,
+                  1, 0, 1, 1,
+                  1, 1, 0, 1,
+                  0, 0, 1, 0), 4, 4, byrow=TRUE)
> ### hier gibt es eigentlich nichts mehr zu orientieren; vor bugfix hat aber applyOrientationRules  Rule 3 angewendet
> amat2 <- t(pcalg:::applyOrientationRules(t(amat), verbose=TRUE))
> stopifnot( all(amat == amat2 ))
> 
> proc.time()
   user  system elapsed 
   1.29    0.17    1.46