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)
> 
> ## Create a weighted DAG --- Example from ../man/jointIda.Rd  (keep in sync !)
> p <- 6
> V <- as.character(1:p)
> edL <- list(
+   "1" = list(edges=c(3,4), weights=c(1.1,0.3)),
+   "2" = list(edges=c(6),  weights=c(0.4)),
+   "3" = list(edges=c(2,4,6),weights=c(0.6,0.8,0.9)),
+   "4" = list(edges=c(2),weights=c(0.5)),
+   "5" = list(edges=c(1,4),weights=c(0.2,0.7)),
+   "6" = NULL)
> myDAG <- new("graphNEL", nodes=V, edgeL=edL, edgemode="directed") ## true DAG
> myCPDAG <- dag2cpdag(myDAG) ## true CPDAG
> covTrue <- trueCov(myDAG) ## true covariance matrix
> 
> n <- 1000
> ## simulate Gaussian data from the true DAG
> dat <- if (require("mvtnorm")) {
+   set.seed(123)
+   rmvnorm(n, mean=rep(0,p), sigma=covTrue)
+ } else readRDS(system.file(package="pcalg", "external", "N_6_1000.rds"))
Loading required package: mvtnorm
> 
> ## estimate CPDAG -- see  help(pc)
> suffStat <- list(C = cor(dat), n = n)
> pc.fit <- pc(suffStat, indepTest = gaussCItest, p = p, alpha = 0.01, u2pd="relaxed")
> 
> ## Suppose that we know the true CPDAG and covariance matrix
> m1 <- jointIda(x.pos=c(1,2), y.pos=6, covTrue, graphEst=myCPDAG, technique="RRC")
> m2 <- jointIda(x.pos=c(1,2), y.pos=6, covTrue, graphEst=myCPDAG, technique="MCD") ### MCD needed a bugfix
> 
> ##sorting for comparison
> order.m1 <- m1[,order(m1[1,])]
> order.m2 <- m2[,order(m2[1,])]
> tempres1 <- all( all.equal(order.m1, order.m2) )
> if (!tempres1) stop("There is a mismatch with jointIda RRC and MCD!")
> 
> ## Instead of knowing the true CPDAG, it is enough to know only
> ## all jointly valid parent sets of the intervention variables ## to use RRC or MCD
> ajv.pasets <- list(list(5, c(3,4)),
+                    list(integer(0), c(3,4)),
+                    list(3, c(3,4)))
> m3 <- jointIda(x.pos=c(1,2), y.pos=6, covTrue, all.pasets=ajv.pasets, technique="RRC")
> m4 <- jointIda(x.pos=c(1,2), y.pos=6, covTrue, all.pasets=ajv.pasets, technique="MCD") 
> 
> ##sorting for comparison
> order.m3 <- m3[,order(m3[1,])]
> order.m4 <- m4[,order(m4[1,])]
> tempres2 <- all( all.equal(order.m3, order.m4) )
> if (!tempres2) stop("There is a mismatch with jointIda RRC and MCD!")
> 
> ## From the true DAG, we can compute the true total joint effects
> ## using RRC or MCD
> m5 <- jointIda(x.pos=c(1,2),y.pos=6,covTrue,graphEst=myDAG,technique="RRC")
> m6 <- jointIda(x.pos=c(1,2),y.pos=6,covTrue,graphEst=myDAG,technique="MCD")
> 
> ##sorting for comparison
> order.m5 <- m5[,order(m5[1,])]
> order.m6 <- m6[,order(m6[1,])]
> tempres3 <- all( all.equal(order.m5, order.m6))
> if (!tempres3) stop("There is a mismatch with jointIda  RRC and MCD!")
> 
> ################################################################
> 
> mTrue <- rbind(c(0,0.99,0.99), c(0.4,0.4,0.4))
> Rnd <- 7
> res1 <- (all(round(order.m1,Rnd) == mTrue) &
+          all(round(order.m2,Rnd) == mTrue) &
+          all(round(order.m3,Rnd) == mTrue) &
+          all(round(order.m4,Rnd) == mTrue) )
> 
> if(!res1) stop("Test in jointIda: True causal effects were not recovered!")
> 
> 
> #############################################################################
> 
> ## jointIda also works when x.pos has length 1 and in the following example
> ## it gives the same result as ida() (see Note)
> ##
> ## When the CPDAG is known
> v1 <- sort(round(jointIda(x.pos=1,y.pos=6,covTrue,graphEst=myCPDAG,
+                      technique="RRC"), Rnd))
> v2 <- sort(round(jointIda(x.pos=1,y.pos=6,covTrue,graphEst=myCPDAG,
+                technique="MCD"), Rnd))    ###EMA MCD problem as above
> v3 <- sort(round(ida(x.pos=1,y.pos=6,covTrue,graphEst=myCPDAG,
+           method="global"), Rnd))   ###EMA there is an issue in v3 with lm.cov ??? i get the error:
> 
> res2 <- (all(v1==v2) & all(v2==v3) & all(v3==v1))
> 
> if(!res2) stop("Test in jointIda: ida() and jointIda() don't produce the same result!")
> 
> #########################################
> ## test extract parent set
> #########################################
> # library(digest)
> # 
> # 
> # ## compare outputs
> # jidaCompare <- function(a,b) {
> #   digest(a) ## has "integer (0)"
> #   digest(b) ## has "named integer(0)"
> #   for (i in 1:length(b)) {
> #     for (j in 1:length(b[[i]])) {
> #       if (length(b[[i]][[j]]) == 0) b[[i]][[j]] <- integer(0)
> #     }
> #   }
> #   
> #   for (i in 1:length(a)) {
> #     for (j in 1:length(a[[i]])) {
> #       if (length(a[[i]][[j]]) == 0) a[[i]][[j]] <- integer(0)
> #     }
> #   }
> #   
> #   all(unique(sort(sapply(a,digest))) == unique(sort(sapply(b,digest))))
> # }
> # 
> # 
> # ##testing
> # ##b[[5]]
> # ##a[[14]]
> # 
> # ##c <- a[[14]]
> # ##c2 <- b[[5]]
> # 
> # ###testing
> # ## possible_p <- c(seq(20,100,by=10))
> # possible_p <- c(seq(20,50,by=5))
> # 
> # ## possible_neighb_size <- c(3:15)
> # ## possible_neighb_size <- c(3:10)
> # possible_neighb_size <- c(3:8)
> # 
> # ## or choose prob and generate neighborhood size 
> # possible_prob <- seq(0.1,0.6,by=0.1)
> # 
> # ##for function msep
> # 
> # ##for joint IDA MCD:
> # ##if you want to run locally set nc to 1!!
> # nc <-  1
> # registerDoParallel(nc) ## Request nc processors in parallel
> # num_setings <-20
> # num_rep <- 5
> # 
> # ##amount of bg knowledge in percent
> # p_bg <- seq(0,1,by=0.1)
> # 
> # ##  start.time <- Sys.time()
> # 
> # seed1 <- 1 
> # resFin <- foreach(r=seed1: num_setings, .packages = 'pcalg') %dopar%{
> #   set.seed(r)
> #   ##Then we can choose one setting:
> #   
> #   p <-  sample(possible_p,1)
> #   neigh <- sample(possible_neighb_size,1)
> #   prob <- round(neigh/(p-1),3)
> #   res <- rep(FALSE, num_rep)
> #   
> #   cat("r=",r,", p=",p,", nb=",neigh,"\n")     
> #   
> #   ## then for every setting selected we can generate i.e. 20 graph
> #   for (d in 1: num_rep){
> #     ##i left the seed inside for now
> #     ## i am not sure which seed option makes it easier for us
> #     ## feel free to move this around
> #     #seed <- sample(possible_seed,1)
> #     #set.seed(seed)
> #     
> #     ##get graph
> #     g <- pcalg:::randomDAG(p, prob)  ## true DAG
> #     cpdag <- dag2cpdag(g) ## true CPDAG
> #     
> #     ## get adjacency matrix
> #     cpdag.amat <- t(as(cpdag,"matrix"))
> #     dag.amat <- t(as(g,"matrix"))
> #     dag.amat[dag.amat != 0] <- 1
> #     
> #     #############################
> #     ##                         ##    
> #     ## NEW CODE STARTS HERE!!! ##
> #     ##                         ##    
> #     #############################
> #     
> #     amat.cpdag <- cpdag.amat
> #     #plot(cpdag)
> #     x <- sample(p,3)
> #     
> #     a <- pcalg:::extract.parent.sets(x,t(amat.cpdag)) ## extract valid joint parent sets using pcalg
> #     ##      b <- extract.parent.sets(x,amat.cpdag)          ## and using my own function
> #     b <- b.extract.parent.sets(x,amat.cpdag)          ## and using my own function
> #     ## cat("l.a=",length(a), " - ","l.b=",length(b),"\n")
> #     ## cat("jidaCompare:",jidaCompare(a,b),"\n")
> #     res[d] <- jidaCompare(a,b)
> #     ##cat(res[d],"\n")
> #   }
> #   all(res)
> # }
> # 
> # 
> # stopifnot(all(unlist(resFin)))
> # 
> # print(resFin)
> 
> proc.time()
   user  system elapsed 
   1.31    0.15    1.46