Running 'test_LINGAM.R' [172s]
  Running 'test_addBgKnowledge.R' [159s]
  Running 'test_adjustment.R' [3s]
  Running 'test_ages.R' [164s]
  Running 'test_amat2dag.R' [2s]
  Running 'test_arges.R' [2s]
  Running 'test_backdoor.R' [6s]
  Comparing 'test_backdoor.Rout' to 'test_backdoor.Rout.save' ... OK
  Running 'test_bicscore.R' [157s]
  Running 'test_causalEffect.R' [1s]
  Running 'test_coercion.R' [2s]
  Running 'test_compareGraphs.R' [1s]
  Running 'test_dag2cpdag.R' [2s]
  Running 'test_dag2essgraph.R' [2s]
  Running 'test_displayAmat.R' [151s]
  Running 'test_dsep.R' [150s]
  Running 'test_dsepAM.R' [1s]
  Running 'test_dsepAMTest.R' [2s]
  Running 'test_fci.R' [2s]
  Running 'test_fciPlus.R' [158s]
  Running 'test_gSquareBin.R' [159s]
  Running 'test_gSquareDis.R' [167s]
  Running 'test_gac.R' [3s]
  Running 'test_getNextSet.R' [2s]
  Running 'test_gies.R' [3s]
  Running 'test_ida.R' [165s]
  Running 'test_idaFast.R' [4s]
  Running 'test_isValidGraph.R' [156s]
  Running 'test_jointIda.R' [2s]
  Running 'test_mat2targets.R' [2s]
  Running 'test_optAdjSet.R' [162s]
  Running 'test_opttarget.R' [160s]
  Running 'test_pag2anc.R' [1s]
  Running 'test_pag2conf.R' [1s]
  Running 'test_pag2edge.R' [1s]
  Running 'test_pc.R' [2s]
  Running 'test_pcSelect.R' [3s]
  Running 'test_pcalg2dagitty.R' [2s]
  Running 'test_pcorOrder.R' [1s]
  Running 'test_pdag2allDags.R' [173s]
  Running 'test_pdag2dag.R' [2s]
  Running 'test_possDeAn.R' [171s]
  Running 'test_randDAG.R' [4s]
  Comparing 'test_randDAG.Rout' to 'test_randDAG.Rout.save' ...1a2,16
> R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
> Copyright (C) 2023 The R Foundation for Statistical Computing
> Platform: x86_64-pc-linux-gnu (64-bit)
> 
> R ist freie Software und kommt OHNE JEGLICHE GARANTIE.
> Sie sind eingeladen, es unter bestimmten Bedingungen weiter zu verbreiten.
> Tippen Sie 'license()' or 'licence()' für Details dazu.
> 
> R ist ein Gemeinschaftsprojekt mit vielen Beitragenden.
> Tippen Sie 'contributors()' für mehr Information und 'citation()',
> um zu erfahren, wie R oder R packages in Publikationen zitiert werden können.
> 
> Tippen Sie 'demo()' für einige Demos, 'help()' für on-line Hilfe, oder
> 'help.start()' für eine HTML Browserschnittstelle zur Hilfe.
> Tippen Sie 'q()', um R zu verlassen.
56c71
< Number of Edges = 21 
---
> Number of Edges = 0 
61c76
< Number of Edges = 13 
---
> Number of Edges = 12 
66c81
< Number of Edges = 16 
---
> Number of Edges = 13 
71c86
< Number of Edges = 19 
---
> Number of Edges = 16 
90,93c105,108
<  $ barabasi : chr "4"
<  $ geometric: chr [1:3] "3" "4" "7"
<  $ power    : chr [1:4] "2" "7" "8" "9"
<  $ interEr  : chr [1:2] "3" "7"
---
>  $ barabasi : chr [1:10] "1" "2" "3" "4" ...
>  $ geometric: chr [1:3] "2" "6" "10"
>  $ power    : chr [1:4] "2" "4" "5" "10"
>  $ interEr  : chr [1:2] "5" "6"
108,111c123,126
<  $ barabasi : chr [1:3] "1" "6" "7"
<  $ geometric: chr [1:3] "5" "6" "9"
<  $ power    : chr [1:3] "4" "5" "10"
<  $ interEr  : chr [1:2] "2" "4"
---
>  $ barabasi : chr [1:10] "1" "2" "3" "4" ...
>  $ geometric: chr [1:4] "1" "3" "5" "10"
>  $ power    : chr [1:4] "1" "6" "7" "8"
>  $ interEr  : chr [1:3] "1" "3" "8"
145c160
< Number of Edges = 13 
---
> Number of Edges = 0 
150c165
< Number of Edges = 6 
---
> Number of Edges = 12 
155c170
< Number of Edges = 11 
---
> Number of Edges = 13 
160c175
< Number of Edges = 10 
---
> Number of Edges = 11 
174a190
> Lade nötiges Paket: Matrix
240,242c256,258
< 1  . . . . | . . . | .
< 2  . . . | | . . | | .
< 3  . | . | . . . | | .
---
> 1  . . . . . . . . . .
> 2  . . . . . . . . . .
> 3  . . . . . . . . . .
244,249c260,265
< 5  . . . . . . . | . .
< 6  . | | . | . . . . |
< 7  . | . | . . . . . .
< 8  . . . | . . . . . .
< 9  . . . | . . . . . .
< 10 . . | . | . . . . .
---
> 5  . . . . . . . . . .
> 6  . . . . . . . . . .
> 7  . . . . . . . . . .
> 8  . . . . . . . . . .
> 9  . . . . . . . . . .
> 10 . . . . . . . . . .
255,264c271,280
< 1  . . | . . . . . . .
< 2  . . . . . . | . . .
< 3  . . . . . . . . . .
< 4  . . . . . . . . . .
< 5  . | . . . . | . . |
< 6  . . . . . . . . . |
< 7  . . . . . . . . . .
< 8  | . | | . . . . . .
< 9  | . | . . . . | . .
< 10 . . . . . . | . . .
---
> 1  . . . . . | . | . .
> 2  . . . . . . . . . .
> 3  . | . | . . . . . .
> 4  . | . . . . . . . .
> 5  . . . . . . . . | .
> 6  . . . . . . . . . .
> 7  . | . . . | . . . .
> 8  . . . . . | | . . .
> 9  . . . | . . | . . .
> 10 . . . . . . . . . .
270c286
< 1  . . . . . . | . . .
---
> 1  . . | . | . . . | .
272,279c288,295
< 3  | | . . . . | | | .
< 4  | . | . . | | . . .
< 5  . . | . . . . . . .
< 6  . . . . . . | . . .
< 7  . . . . . . . . . .
< 8  . . . . . . . . . .
< 9  . . . . . . . . . .
< 10 | . | . . . | . | .
---
> 3  . . . | . . . . . .
> 4  . . . . . . . . . .
> 5  . . . . . . . . . .
> 6  . | . . . . . . | .
> 7  . . . | . . . . | |
> 8  . . . | . . . . | .
> 9  . . | . . . . . . |
> 10 . . . . . . . . . .
285,294c301,310
< 1  . . . . . . . . . |
< 2  | . | . . . . . | |
< 3  . . . . . . . . . .
< 4  . . . . | | | | . .
< 5  . . . . . . | . . .
< 6  . . . . | . | | . .
< 7  . . . . . . . . . .
< 8  | . . . | . . . . .
< 9  | . | . . . . . . |
< 10 . . | . . . . . . .
---
> 1  . . . . | . . . . |
> 2  . . . . . | . . | .
> 3  . | . . . | | . | .
> 4  . . . . | . . . . .
> 5  . . . . . . . . . .
> 6  . . . . . . . . . .
> 7  . . . . . | . . . .
> 8  . . . | | | . . . |
> 9  . . . . . | . . . .
> 10 . . . . | . . . . .
369,372c385,388
< 1  . . . . . | . . . . . .
< 2  | . . . . . . . . . | .
< 3  | . . . . . . . . . . |
< 4  | . . . . . . . . . . .
---
> 1  . . . . . . . . . . . .
> 2  . . . . . . . . . . . .
> 3  . . . . . . . . . . . .
> 4  . . . . . . . . . . . .
375,379c391,395
< 7  | . . . . . . . . . . .
< 8  | . . . . . . . . . . .
< 9  . . . . | | . . . . . .
< 10 . . . | . . . . | . . .
< 11 . . | . . . . . . . . .
---
> 7  . . . . . . . . . . . .
> 8  . . . . . . . . . . . .
> 9  . . . . . . . . . . . .
> 10 . . . . . . . . . . . .
> 11 . . . . . . . . . . . .
386,387c402,403
< 1  . . . . . . . . | . . .
< 2  . . . | . . . . . | . .
---
> 1  . . . . . . . . . . . .
> 2  . . . . . . . . . . . .
389c405
< 4  . . . . . . . . . | . .
---
> 4  . . | . . . . . . | . .
391c407
< 6  . . . . . . . . . . . |
---
> 6  . | . . . . . | . | . .
393,394c409,410
< 8  . . . . . . . . . . . .
< 9  . . . . . . . . . . . .
---
> 8  . | . . | . . . . | . .
> 9  . . | | . . . . . . . .
396,397c412,413
< 11 . . | . . . . . . . . .
< 12 . . . . . . . . . . . .
---
> 11 . . . . . . . . . . . .
> 12 . . | . . . . . | . . .
403,407c419,423
< 1  . . . . . | | . . . . |
< 2  . . . . . . | . . . . .
< 3  | . . . . . . . . . . .
< 4  . . . . . . | . . . . .
< 5  | . . . . . . . . . . .
---
> 1  . . . . . . . . . . . .
> 2  . . . . . . . . . . . .
> 3  . . . . . . . . . . . .
> 4  . . . . . . . . . . . .
> 5  . . . . . . . . . . . .
410,414c426,430
< 8  . . . . . . | . | . . .
< 9  . . . . . . . . . . . .
< 10 . . . . . . | . . . . .
< 11 . . . . . . | . . . . .
< 12 . . . . . . . . . . . .
---
> 8  | | . | | . . . . . . |
> 9  . . | . . . . . . . . .
> 10 . | . . . . . | . . . .
> 11 . . . . . . . . | . . .
> 12 . | . . . | | . | . . .
420,422c436,438
< 1  . . . . . . . . . . . .
< 2  . . . . | . . . . . . .
< 3  . . . . . . | . . . | .
---
> 1  . . . . . . . . . . . |
> 2  | . . . . . . | . . . .
> 3  . . . . . . . . . . . |
425,431c441,447
< 6  . . . . . . . . . . . .
< 7  . . . . | . . . . . . .
< 8  . . . . . . . . . . . .
< 9  | . . . . . . . . . . .
< 10 . . . . . . | | . . . .
< 11 . . . . | . . . . . . .
< 12 | . . . . | . . . . . .
---
> 6  | . . . . . . . . . . |
> 7  . . . . . . . . . . . .
> 8  | . . . . . . . . . . |
> 9  . . . | | . . . . . . .
> 10 | . . . . . . . . . . .
> 11 . . . . . . . . . . . .
> 12 . . . . . . . . . . . .
523,525c539,541
< 1  . . . . | . . . | .
< 2  . . . | | . . | | .
< 3  . | . | . . . | | .
---
> 1  . . . . . . . . . .
> 2  . . . . . . . . . .
> 3  . . . . . . . . . .
527,532c543,548
< 5  . . . . . . . | . .
< 6  . | | . | . . . . |
< 7  . | . | . . . . . .
< 8  . . . | . . . . . .
< 9  . . . | . . . . . .
< 10 . . | . | . . . . .
---
> 5  . . . . . . . . . .
> 6  . . . . . . . . . .
> 7  . . . . . . . . . .
> 8  . . . . . . . . . .
> 9  . . . . . . . . . .
> 10 . . . . . . . . . .
538,547c554,563
< 1  . . | . . . . . . .
< 2  . . . . . . | . . .
< 3  . . . . . . . . . .
< 4  . . . . . . . . . .
< 5  . | . . . . | . . |
< 6  . . . . . . . . . |
< 7  . . . . . . . . . .
< 8  | . | | . . . . . .
< 9  | . | . . . . | . .
< 10 . . . . . . | . . .
---
> 1  . . . . . | . | . .
> 2  . . . . . . . . . .
> 3  . | . | . . . . . .
> 4  . | . . . . . . . .
> 5  . . . . . . . . | .
> 6  . . . . . . . . . .
> 7  . | . . . | . . . .
> 8  . . . . . | | . . .
> 9  . . . | . . | . . .
> 10 . . . . . . . . . .
553c569
< 1  . . . . . . | . . .
---
> 1  . . | . | . . . | .
555,562c571,578
< 3  | | . . . . | | | .
< 4  | . | . . | | . . .
< 5  . . | . . . . . . .
< 6  . . . . . . | . . .
< 7  . . . . . . . . . .
< 8  . . . . . . . . . .
< 9  . . . . . . . . . .
< 10 | . | . . . | . | .
---
> 3  . . . | . . . . . .
> 4  . . . . . . . . . .
> 5  . . . . . . . . . .
> 6  . | . . . . . . | .
> 7  . . . | . . . . | |
> 8  . . . | . . . . | .
> 9  . . | . . . . . . |
> 10 . . . . . . . . . .
568,577c584,593
< 1  . . . . . . . . . |
< 2  | . | . . . . . | |
< 3  . . . . . . . . . .
< 4  . . . . | | | | . .
< 5  . . . . . . | . . .
< 6  . . . . | . | | . .
< 7  . . . . . . . . . .
< 8  | . . . | . . . . .
< 9  | . | . . . . . . |
< 10 . . | . . . . . . .
---
> 1  . . . . | . . . . |
> 2  . . . . . | . . | .
> 3  . | . . . | | . | .
> 4  . . . . | . . . . .
> 5  . . . . . . . . . .
> 6  . . . . . . . . . .
> 7  . . . . . | . . . .
> 8  . . . | | | . . . |
> 9  . . . . . | . . . .
> 10 . . . . | . . . . .
652,655c668,671
< 1  . . . . . | . . . . . .
< 2  | . . . . . . . . . | .
< 3  | . . . . . . . . . . |
< 4  | . . . . . . . . . . .
---
> 1  . . . . . . . . . . . .
> 2  . . . . . . . . . . . .
> 3  . . . . . . . . . . . .
> 4  . . . . . . . . . . . .
658,662c674,678
< 7  | . . . . . . . . . . .
< 8  | . . . . . . . . . . .
< 9  . . . . | | . . . . . .
< 10 . . . | . . . . | . . .
< 11 . . | . . . . . . . . .
---
> 7  . . . . . . . . . . . .
> 8  . . . . . . . . . . . .
> 9  . . . . . . . . . . . .
> 10 . . . . . . . . . . . .
> 11 . . . . . . . . . . . .
669,670c685,686
< 1  . . . . . . . . | . . .
< 2  . . . | . . . . . | . .
---
> 1  . . . . . . . . . . . .
> 2  . . . . . . . . . . . .
672c688
< 4  . . . . . . . . . | . .
---
> 4  . . | . . . . . . | . .
674c690
< 6  . . . . . . . . . . . |
---
> 6  . | . . . . . | . | . .
676,677c692,693
< 8  . . . . . . . . . . . .
< 9  . . . . . . . . . . . .
---
> 8  . | . . | . . . . | . .
> 9  . . | | . . . . . . . .
679,680c695,696
< 11 . . | . . . . . . . . .
< 12 . . . . . . . . . . . .
---
> 11 . . . . . . . . . . . .
> 12 . . | . . . . . | . . .
686,690c702,706
< 1  . . . . . | | . . . . |
< 2  . . . . . . | . . . . .
< 3  | . . . . . . . . . . .
< 4  . . . . . . | . . . . .
< 5  | . . . . . . . . . . .
---
> 1  . . . . . . . . . . . .
> 2  . . . . . . . . . . . .
> 3  . . . . . . . . . . . .
> 4  . . . . . . . . . . . .
> 5  . . . . . . . . . . . .
693,697c709,713
< 8  . . . . . . | . | . . .
< 9  . . . . . . . . . . . .
< 10 . . . . . . | . . . . .
< 11 . . . . . . | . . . . .
< 12 . . . . . . . . . . . .
---
> 8  | | . | | . . . . . . |
> 9  . . | . . . . . . . . .
> 10 . | . . . . . | . . . .
> 11 . . . . . . . . | . . .
> 12 . | . . . | | . | . . .
703,705c719,721
< 1  . . . . . . . . . . . .
< 2  . . . . | . . . . . . .
< 3  . . . . . . | . . . | .
---
> 1  . . . . . . . . . . . |
> 2  | . . . . . . | . . . .
> 3  . . . . . . . . . . . |
708,714c724,730
< 6  . . . . . . . . . . . .
< 7  . . . . | . . . . . . .
< 8  . . . . . . . . . . . .
< 9  | . . . . . . . . . . .
< 10 . . . . . . | | . . . .
< 11 . . . . | . . . . . . .
< 12 | . . . . | . . . . . .
---
> 6  | . . . . . . . . . . |
> 7  . . . . . . . . . . . .
> 8  | . . . . . . . . . . |
> 9  . . . | | . . . . . . .
> 10 | . . . . . . . . . . .
> 11 . . . . . . . . . . . .
> 12 . . . . . . . . . . . .
  Running 'test_randomDAG.R' [169s]
  Running 'test_rfci.R' [165s]
  Running 'test_rmvDAG.R' [2s]
  Running 'test_searchAM.R' [177s]
  Running 'test_shd.R' [167s]
  Running 'test_skeleton.R' [3s]
  Running 'test_udag2pag.R' [1s]
  Running 'test_udag2pdag.R' [169s]
  Running 'test_wgtMatrix.R' [157s]
Running the tests in 'tests/test_LINGAM.R' failed.
Complete output:
  > library(pcalg)
  > ## source("/u/kalischm/research/packages/LINGAM/R/lingamFuns.R")
  > 
  > ##--> showProc.time(), assertError(), relErrV(), ...
  > R.home(); sessionInfo() # helping package maintainers to debug ...
  [1] "D:/RCompile/recent/R"
  R Under development (unstable) (2024-03-22 r86169 ucrt)
  Platform: x86_64-w64-mingw32/x64
  Running under: Windows Server 2022 x64 (build 20348)
  
  Matrix products: default
  
  
  locale:
  [1] LC_COLLATE=C                 LC_CTYPE=German_Germany.utf8
  [3] LC_MONETARY=C                LC_NUMERIC=C                
  [5] LC_TIME=C                   
  
  time zone: Europe/Berlin
  tzcode source: internal
  
  attached base packages:
  [1] stats     graphics  grDevices utils     datasets  methods   base     
  
  other attached packages:
  [1] pcalg_2.7-11
  
  loaded via a namespace (and not attached):
   [1] graph_1.80.0        igraph_2.0.3        DEoptimR_1.1-3     
   [4] clue_0.3-65         magrittr_2.0.3      fastICA_1.2-4      
   [7] abind_1.4-5         BiocGenerics_0.48.1 pkgconfig_2.0.3    
  [10] stats4_4.4.0        lifecycle_1.0.4     cli_3.6.2          
  [13] sfsmisc_1.1-17      ggm_2.5.1           corpcor_1.6.10     
  [16] robustbase_0.99-2   compiler_4.4.0      RBGL_1.78.0        
  [19] tools_4.4.0         cluster_2.1.6       bdsmatrix_1.3-7    
  [22] Rcpp_1.0.12         BiocManager_1.30.22 rlang_1.1.3        
  > .libPaths()
  [1] "D:/temp/RtmpmCFgj4/RLIBS_1177073ea38fa"
  [2] "D:/RCompile/recent/R/library"          
  > packageDescription("pcalg")
  Package: pcalg
  Version: 2.7-11
  Date: 2024-02-12
  Title: Methods for Graphical Models and Causal Inference
  Description: Functions for causal structure learning and causal
          inference using graphical models. The main algorithms for
          causal structure learning are PC (for observational data
          without hidden variables), FCI and RFCI (for observational data
          with hidden variables), and GIES (for a mix of data from
          observational studies (i.e. observational data) and data from
          experiments involving interventions (i.e. interventional data)
          without hidden variables). For causal inference the IDA
          algorithm, the Generalized Backdoor Criterion (GBC), the
          Generalized Adjustment Criterion (GAC) and some related
          functions are implemented. Functions for incorporating
          background knowledge are provided.
  Maintainer: Markus Kalisch <kalisch@stat.math.ethz.ch>
  Authors@R: c(person("Markus","Kalisch",
          email="kalisch@stat.math.ethz.ch", role=c("aut","cre")),
          person("Alain", "Hauser", role="aut"),
          person("Martin","Maechler", role="aut"), person("Diego",
          "Colombo", role="ctb"), person("Doris", "Entner", role="ctb"),
          person("Patrik", "Hoyer", role="ctb"), person("Antti",
          "Hyttinen", role="ctb"), person("Jonas", "Peters", role="ctb"),
          person("Nicoletta", "Andri", role="ctb"), person("Emilija",
          "Perkovic", role="ctb"), person("Preetam", "Nandy",
          role="ctb"), person("Philipp", "Ruetimann",
          role="ctb"),person("Daniel", "Stekhoven", role="ctb"),
          person("Manuel", "Schuerch", role="ctb"), person("Marco",
          "Eigenmann", role="ctb"), person("Leonard", "Henckel", role =
          "ctb"),person("Joris", "Mooij", role = "ctb"))
  Depends: R (>= 3.5.0)
  LinkingTo: Rcpp (>= 0.11.0), RcppArmadillo, BH
  Imports: stats, graphics, utils, methods, abind, graph, RBGL, igraph,
          ggm, corpcor, robustbase, vcd, Rcpp, bdsmatrix, sfsmisc (>=
          1.0-26), fastICA, clue
  Suggests: MASS, Matrix, Rgraphviz, mvtnorm, huge, ggplot2, dagitty
  NeedsCompilation: yes
  Encoding: UTF-8
  License: GPL (>= 2)
  URL: https://pcalg.r-forge.r-project.org/
  Packaged: 2024-02-12 10:30:09 UTC; kalischm
  Author: Markus Kalisch [aut, cre], Alain Hauser [aut], Martin Maechler
          [aut], Diego Colombo [ctb], Doris Entner [ctb], Patrik Hoyer
          [ctb], Antti Hyttinen [ctb], Jonas Peters [ctb], Nicoletta
          Andri [ctb], Emilija Perkovic [ctb], Preetam Nandy [ctb],
          Philipp Ruetimann [ctb], Daniel Stekhoven [ctb], Manuel
          Schuerch [ctb], Marco Eigenmann [ctb], Leonard Henckel [ctb],
          Joris Mooij [ctb]
  Repository: CRAN
  Date/Publication: 2024-02-12 12:00:05 UTC
  Built: R 4.4.0; x86_64-w64-mingw32; 2024-03-23 03:42:08 UTC; windows
  Archs: x64
  
  -- File: D:/RCompile/CRANpkg/lib/4.4/pcalg/Meta/package.rds 
  > packageDescription("Matrix")
  Package: Matrix
  Version: 1.7-0
  VersionNote: do also bump src/version.h, inst/include/Matrix/version.h
  Date: 2024-03-16
  Priority: recommended
  Title: Sparse and Dense Matrix Classes and Methods
  Description: A rich hierarchy of sparse and dense matrix classes,
          including general, symmetric, triangular, and diagonal matrices
          with numeric, logical, or pattern entries.  Efficient methods
          for operating on such matrices, often wrapping the 'BLAS',
          'LAPACK', and 'SuiteSparse' libraries.
  License: GPL (>= 2) | file LICENCE
  URL: https://Matrix.R-forge.R-project.org
  BugReports: https://R-forge.R-project.org/tracker/?atid=294&group_id=61
  Contact: Matrix-authors@R-project.org
  Authors@R: c(person("Douglas", "Bates", role = "aut", comment = c(ORCID
          = "0000-0001-8316-9503")), person("Martin", "Maechler", role =
          c("aut", "cre"), email = "mmaechler+Matrix@gmail.com", comment
          = c(ORCID = "0000-0002-8685-9910")), person("Mikael", "Jagan",
          role = "aut", comment = c(ORCID = "0000-0002-3542-2938")),
          person("Timothy A.", "Davis", role = "ctb", comment = c(ORCID =
          "0000-0001-7614-6899", "SuiteSparse libraries", "collaborators
          listed in dir(system.file(\"doc\", \"SuiteSparse\",
          package=\"Matrix\"), pattern=\"License\", full.names=TRUE,
          recursive=TRUE)")), person("George", "Karypis", role = "ctb",
          comment = c(ORCID = "0000-0003-2753-1437", "METIS library",
          "Copyright: Regents of the University of Minnesota")),
          person("Jason", "Riedy", role = "ctb", comment = c(ORCID =
          "0000-0002-4345-4200", "GNU Octave's condest() and
          onenormest()", "Copyright: Regents of the University of
          California")), person("Jens", "Oehlschlägel", role = "ctb",
          comment = "initial nearPD()"), person("R Core Team", role =
          "ctb", comment = "base R's matrix implementation"))
  Depends: R (>= 4.4.0), methods
  Imports: grDevices, graphics, grid, lattice, stats, utils
  Suggests: MASS, datasets, sfsmisc, tools
  Enhances: SparseM, graph
  LazyData: no
  LazyDataNote: not possible, since we use data/*.R and our S4 classes
  BuildResaveData: no
  Encoding: UTF-8
  NeedsCompilation: yes
  Packaged: 2024-03-19 17:15:14 UTC; maechler
  Author: Douglas Bates [aut] (<https://orcid.org/0000-0001-8316-9503>),
          Martin Maechler [aut, cre]
          (<https://orcid.org/0000-0002-8685-9910>), Mikael Jagan [aut]
          (<https://orcid.org/0000-0002-3542-2938>), Timothy A. Davis
          [ctb] (<https://orcid.org/0000-0001-7614-6899>, SuiteSparse
          libraries, collaborators listed in dir(system.file("doc",
          "SuiteSparse", package="Matrix"), pattern="License",
          full.names=TRUE, recursive=TRUE)), George Karypis [ctb]
          (<https://orcid.org/0000-0003-2753-1437>, METIS library,
          Copyright: Regents of the University of Minnesota), Jason Riedy
          [ctb] (<https://orcid.org/0000-0002-4345-4200>, GNU Octave's
          condest() and onenormest(), Copyright: Regents of the
          University of California), Jens Oehlschlägel [ctb] (initial
          nearPD()), R Core Team [ctb] (base R's matrix implementation)
  Maintainer: Martin Maechler <mmaechler+Matrix@gmail.com>
  Repository: CRAN
  Date/Publication: 2024-03-22 09:08:25 UTC
  Built: R 4.4.0; x86_64-w64-mingw32; 2024-03-23 02:52:59 UTC; windows
  Archs: x64
  
  -- File: D:/temp/RtmpmCFgj4/RLIBS_1177073ea38fa/Matrix/Meta/package.rds 
  > cat("doExtras:", (doExtras <- pcalg:::doExtras()), "\n")
  doExtras: FALSE 
  > 
  > 
  > ##################################################
  > ## Exp 1
  > ##################################################
  > set.seed(123)
  > n <- 500
  > eps1 <- sign(rnorm(n)) * sqrt(abs(rnorm(n)))
  > eps2 <- runif(n) - 0.5 #  ~ U[-1/2, 1/2]
  > 
  > X <- cbind(A = eps1 + 0.9*eps2,
  +            B =            eps2)
  > 
  > ## x1 <- x2
  > ## adjacency matrix:
  > ## 0 0
  > ## 1 0
  > (trueDAG <- rbind(c(0,0),
  +                   c(1,0)))
       [,1] [,2]
  [1,]    0    0
  [2,]    1    0
  > 
  > estDAG <- LINGAM(X, verbose = TRUE)
  Performing row permutation, nzdiag*() ...
    (Small dimensionality, using brute-force method.): Done!
  Performing permutation for causal order...
    (Small dimensionality, using brute-force method.): Done!
            [,1]        [,2]
  [1,] 0.0000000 0.001722255
  [2,] 0.8506492 0.000000000
  Causal B nicely triangular. No problems to report here.
  Pruning the network connections...
  Done!
  Warning message:
  In LINGAM(X, verbose = TRUE) : 'LINGAM' is deprecated.
  Use 'lingam' instead.
  See help("Deprecated")
  > 
  > stopifnot(identical(estDAG, LINGAM(X)),
  +           as.integer(estDAG$ Adj) == trueDAG,
  +           all.equal (estDAG$ B, cbind(0, c(0.878188262685122, 0))))
  Warning message:
  In LINGAM(X) : 'LINGAM' is deprecated.
  Use 'lingam' instead.
  See help("Deprecated")
  > 
  > if(doExtras) {
  +     ## using pcalg
  +     n <- nrow(X)
  +     V <- LETTERS[1:ncol(X)] # labels aka node names
  + 
  +     ## estimate CPDAG
  +     pc.fit <- pc(suffStat = list(C = cor(X), n = n),
  +                  indepTest = gaussCItest, ## indep.test: partial correlations
  +                  alpha = 0.01, labels = colnames(X))
  +     if (require(Rgraphviz)) {
  +         plot(pc.fit, main = "Estimated CPDAG")
  +     }
  + }
  > 
  > ##################################################
  > ## Exp 2
  > ##################################################
  > set.seed(123)
  > n <- 500
  > eps1 <- sign(rnorm(n)) * sqrt(abs(rnorm(n)))
  > eps2 <- runif(n) - 0.5
  > eps3 <- sign(rnorm(n)) * abs(rnorm(n))^(1/3)
  > eps4 <- rnorm(n)^2
  > 
  > x1 <- eps1 + 0.9*eps2
  > x2 <-            eps2
  > x3 <-        0.8*eps2 + eps3
  > x4 <- -0.9*x3 - x1 + eps4
  > 
  > X <- cbind(U = x1, V = x2, W = x3, Y = x4)
  > 
  > trueDAG <- cbind(c(0,1,0,0),c(0,0,0,0),c(0,1,0,0),c(1,0,1,0))
  > ## x4 <- x3 <- x2 -> x1 -> x4
  > ## adjacency matrix:
  > ## 0 0 0 1
  > ## 1 0 1 0
  > ## 0 0 0 1
  > ## 0 0 0 0
  > 
  > estDAG <- LINGAM(X, verbose = TRUE)
  Performing row permutation, nzdiag*() ...
    (Small dimensionality, using brute-force method.): Done!
  Performing permutation for causal order...
    (Small dimensionality, using brute-force method.): Done!
              [,1]        [,2]         [,3]        [,4]
  [1,]  0.00000000  0.01185382  0.009724318 0.010054538
  [2,]  0.89993102  0.00000000 -0.037587250 0.003634021
  [3,]  0.94576470 -0.04663875  0.000000000 0.018116176
  [4,] -0.05985199 -0.88907479 -1.057332765 0.000000000
  Causal B nicely triangular. No problems to report here.
  Pruning the network connections...
  Done!
  Warning message:
  In LINGAM(X, verbose = TRUE) : 'LINGAM' is deprecated.
  Use 'lingam' instead.
  See help("Deprecated")
  > 
  > B.est <- rbind(c(0, 0.986119553, 0, 0),
  +                c(0, 0, 0, 0),
  +                c(0, 0.89198226, 0, 0),
  +                c(-0.987301824, 0, -0.890961952, 0))
  > 
  > stopifnot(as.integer(estDAG$Adj) == trueDAG,
  +           all.equal(estDAG$B, B.est, tol=1e-9))
  > 
  > if(doExtras) {
  +     ## using pcalg
  +     n <- nrow(X)
  +     V <- colnames(X) # labels aka node names
  + 
  +     ## estimate CPDAG
  +     pc.fit <- pc(suffStat = list(C = cor(X), n = n),
  +                  indepTest = gaussCItest, ## indep.test: partial correlations
  +                  alpha=0.01, labels = V, verbose = FALSE)
  +     if (require(Rgraphviz)) {
  +         plot(pc.fit, main = "Estimated CPDAG")
  +     }
  + }
  > 
  > ## if(!doExtras && !interactive()) quit("no")
  > 
  > ### More tests for higher dimensions
  > 
  > ### p = 8 -------- Example 3 -----------
  > 
  > set.seed(127)
  > n <- 2000
  > x1 <- eps1 <- sign(rnorm(n)) * sqrt(abs(rnorm(n)))
  > x2 <- eps2 <- runif(n) - 0.5
  > x3 <- eps3 <- sign(rnorm(n)) * abs(rnorm(n))^(1/3)
  > x4 <- eps4 <- rnorm(n)^2
  > Z <- rnorm(n); eps5 <- sign(Z) * sqrt(abs(Z))
  > Z <- rnorm(n); eps6 <- sign(Z) * sqrt(abs(Z))
  > Z <- rnorm(n); eps7 <- sign(Z) * sqrt(abs(Z))
  > Z <- rnorm(n); eps8 <- sign(Z) * sqrt(abs(Z))
  > 
  > x5 <- 7/8*x1 - 3/4*x2 + 3/4*x3 + eps5
  > x6 <- 0.8*x4                   + eps6
  > x7 <- 3/4*x5 - 7/8*x6          + eps7
  > x8 <- .9*x7                    + eps8
  > 
  > X <- cbind(x1,x2,x3,x4,x5,x6,x7,x8, deparse.level = 2)
  > ## (x1, x2, x3) -> x5 -> x7 <- x6 <- x4;  x7 -> x8
  > ## adjacency matrix:
  > ##     1 2 3 4 5 6 7 8
  > ## x1  . . . . 1 . . .
  > ## x2  . . . . 1 . . .
  > ## x3  . . . . 1 . . .
  > ## x4  . . . . . 1 . .
  > ## x5  . . . . . . 1 .
  > ## x6  . . . . . . 1 .
  > ## x7  . . . . . . . 1
  > ## x8  . . . . . . . .
  > 
  > ## true DAG :
  > . <- 0
  > trDAG3 <- rbind(
  +     c(., ., ., ., 1, ., ., .),
  +     c(., ., ., ., 1, ., ., .),
  +     c(., ., ., ., 1, ., ., .),
  +     c(., ., ., ., ., 1, ., .),
  +     c(., ., ., ., ., ., 1, .),
  +     c(., ., ., ., ., ., 1, .),
  +     c(., ., ., ., ., ., ., 1),
  +     c(., ., ., ., ., ., ., .))
  > 
  > estB.3 <- rbind(
  +     c(., ., ., ., ., ., ., .),
  +     c(., ., ., ., ., ., ., .),
  +     c(., ., ., ., ., ., ., .),
  +     c(., ., ., ., ., ., ., .),
  +     c(.831899433, -.737954104, 0.725137273, ., ., ., ., .),
  +     c(., ., ., 0.788185348, ., ., ., .),
  +     c(., ., ., ., 0.774490692, -0.886143314, ., .),
  +     c(., ., ., ., ., ., 0.900617843, .))
  > 
  > eDAG3 <- LINGAM(X, verbose = 2)
  Centering
  Whitening
  Symmetric FastICA using logcosh approx. to neg-entropy function
  Iteration 1 tol = 0.7545483
  Iteration 2 tol = 0.1655535
  Iteration 3 tol = 0.04065082
  Iteration 4 tol = 0.0001496655
  Iteration 5 tol = 4.919493e-07
  Iteration 6 tol = 3.187801e-09
  Iteration 7 tol = 2.069445e-11
  Iteration 8 tol = 2.204903e-13
  Iteration 9 tol = 4.440892e-15
  Performing row permutation, nzdiag*() ...
    (Small dimensionality, using brute-force method.): Done!
  Performing permutation for causal order...
    (Small dimensionality, using brute-force method.): Done!
              [,1]         [,2]         [,3]         [,4]          [,5]
  [1,]  0.00000000  0.002919480 -0.006553738 -0.003062375  0.0008261561
  [2,]  0.02411382  0.000000000 -0.005439635  0.018401871 -0.0080918882
  [3,] -0.08248625  0.007102494  0.000000000 -0.032518681  0.0119294275
  [4,] -0.74202632  0.710888315  0.841560130  0.000000000  0.0104074191
  [5,] -0.13439282 -0.029284805 -0.044231816  0.039214627  0.0000000000
  [6,] -0.06915424  0.012373744 -0.013669545  0.002144062  0.7893267569
  [7,]  0.04342903 -0.043525297 -0.005390509  0.751614899  0.0050000676
  [8,]  0.01872415 -0.010379778 -0.003082474  0.024731302  0.0074836799
               [,6]          [,7]         [,8]
  [1,]  0.003732972  0.0099985514 -0.005287321
  [2,] -0.001156025 -0.0141942364  0.001135602
  [3,]  0.012845090  0.0307525826 -0.018495408
  [4,] -0.001767195 -0.0005324746  0.013116991
  [5,] -0.032577376 -0.0289342038  0.020539625
  [6,]  0.000000000 -0.0182766705  0.003333580
  [7,] -0.864727090  0.0000000000  0.011545886
  [8,] -0.018094416  0.8841932792  0.000000000
  Causal B nicely triangular. No problems to report here.
  Pruning the network connections...
  Done!
  Warning message:
  In LINGAM(X, verbose = 2) : 'LINGAM' is deprecated.
  Use 'lingam' instead.
  See help("Deprecated")
  > 
  > stopifnot(trDAG3 == eDAG3$Adj,
  +           with(eDAG3, all(t(B != 0) == Adj)),
  +           all.equal(eDAG3$B, estB.3, tol=1e-9))
  > 
  > 
  > ### p = 10 -------- Example 4 -----------
  > 
  > ### using same x1..,x4, and, eps5 .. eps8  as in Ex. 3
  > 
  > Z <- rnorm(n); eps9  <- sign(Z) * abs(Z)^(1/3)
  > Z <- rnorm(n); eps10 <- sign(Z) * abs(Z)^0.25
  > 
  > x5 <-  7/8*x1 - 3/4*x2          + eps5
  > x6 <-  0.8*x2 - 7/8*x3          + eps6
  > x7 <- -7/8*x4                   + eps7
  > x8 <-  0.9*x2 - 0.8*x5          + eps8
  > x9 <- -3/4*x6 + 7/8*x7          + eps9
  > x10 <- 3/4*x6 + 0.5*x8 +0.9*x9  + eps10
  > 
  > X <- cbind(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10, deparse.level = 2)
  > 
  > 
  > ## true DAG :
  > . <- 0
  > trDAG4 <- rbind(
  +     # 1  2  3  4  5  6  7  8  9  0
  +     c(., ., ., ., 1, ., ., ., ., .), #  1
  +     c(., ., ., ., 1, 1, ., 1, ., .), #  2
  +     c(., ., ., ., ., 1, ., ., ., .), #  3
  +     c(., ., ., ., ., ., 1, ., ., .), #  4
  +     c(., ., ., ., ., ., ., 1, ., .), #  5
  +     c(., ., ., ., ., ., ., ., 1, 1), #  6
  +     c(., ., ., ., ., ., ., ., 1, .), #  7
  +     c(., ., ., ., ., ., ., ., ., 1), #  8
  +     c(., ., ., ., ., ., ., ., ., 1), #  9
  +     c(., ., ., ., ., ., ., ., ., .)) # 10
  > 
  > estB.4 <- rbind(
  +     c(., ., ., ., .,  ., ., ., ., .),
  +     c(., ., ., ., .,  ., ., ., ., .),
  +     c(., ., ., ., .,  ., ., ., ., .),
  +     c(., ., ., ., .,  ., ., ., ., .),
  +     c(0.831899433, -0.737954104, ., ., .,  ., ., ., ., .),
  +     c(.,            0.687919607, -0.863361084, ., .,  ., ., ., ., .),
  +     c(.,             .,            .,        -0.878305407, .,  ., ., ., ., .),
  +     c(.,            0.864092647,   .,          .,        -0.77902333,  ., ., ., ., .),
  +     c(.,             .,            .,          .,          .,        -0.780116888, 0.929828083, ., ., .),
  +     c(.,             .,            .,          .,          .,         0.72436897,   .,     0.502210828, 0.913644804, .))
  > 
  > 
  > eDAG4 <- LINGAM(X, verbose = TRUE)
  Performing row permutation, nzdiag*() ...
    (Using the Hungarian algorithm.): Done!
  Performing permutation for causal order...
    (Using pruning algorithm.): Done!
               [,1]         [,2]          [,3]         [,4]         [,5]
   [1,]  0.00000000 -0.006326785  0.0029431058 -0.003816060  0.004737111
   [2,] -0.06430368  0.000000000 -0.0088146934 -0.035995357  0.025808763
   [3,]  0.01297078 -0.006196360  0.0000000000  0.009575279 -0.012113841
   [4,] -0.75694409  0.856720301 -0.0445599745  0.000000000  0.010837351
   [5,] -0.12624203 -0.044575737 -0.0254719709  0.047144333  0.000000000
   [6,]  0.72762061 -0.014204699 -0.8705396282 -0.006159921 -0.013464525
   [7,]  0.02750857 -0.005321567 -0.0327098113  0.019033795 -0.842958526
   [8,]  0.92330050 -0.003411800 -0.0048993505 -0.787883898 -0.010066586
   [9,] -0.03478589  0.004801767 -0.0113849437 -0.011812608  0.018412819
  [10,]  0.09828179 -0.004313719 -0.0005325424  0.016366714  0.020698254
                [,6]         [,7]         [,8]          [,9]         [,10]
   [1,] -0.006504593  0.005701682 -0.009052809 -0.0072221244  0.0076433826
   [2,] -0.009602760  0.031763460 -0.015247794 -0.0133862290 -0.0062761806
   [3,] -0.008791438 -0.002404077 -0.005191894 -0.0244044833  0.0126367391
   [4,] -0.026971906  0.019035635  0.007377126 -0.0190585825  0.0116252362
   [5,] -0.019736049  0.003906628  0.029344715  0.0032010324 -0.0192788022
   [6,]  0.000000000 -0.006334439  0.010147745  0.0009835252 -0.0126648632
   [7,]  0.018792555  0.000000000  0.005190981  0.0225601097  0.0134526632
   [8,] -0.003477516 -0.015750480  0.000000000  0.0015438015 -0.0024023179
   [9,] -0.766570434  0.896089525 -0.002912303  0.0000000000  0.0007192613
  [10,]  0.735553224 -0.001853620  0.501052022  0.9124372652  0.0000000000
  Causal B nicely triangular. No problems to report here.
  Pruning the network connections...
  Done!
  Warning message:
  In LINGAM(X, verbose = TRUE) : 'LINGAM' is deprecated.
  Use 'lingam' instead.
  See help("Deprecated")
  > 
  > stopifnot(trDAG4 == eDAG4$Adj,
  +           with(eDAG4, all(t(B != 0) == Adj)),
  +           all.equal(eDAG4$B, estB.4, tol=1e-9))
  > 
  > proc.time()
     user  system elapsed 
     2.21    0.31    2.53 
Running the tests in 'tests/test_addBgKnowledge.R' failed.
Complete output:
  > library(pcalg)
  > 
  > res <- rep(FALSE, 10)
  > set.seed(123)
  > g <- pcalg::randomDAG(n = 7, prob = 0.3)
  > plot(g)
  
  Attaching package: 'BiocGenerics'
  
  The following objects are masked from 'package:stats':
  
      IQR, mad, sd, var, xtabs
  
  The following objects are masked from 'package:base':
  
      Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
      as.data.frame, basename, cbind, colnames, dirname, do.call,
      duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
      lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
      pmin.int, rank, rbind, rownames, sapply, setdiff, table, tapply,
      union, unique, unsplit, which.max, which.min
  
  > cpdag <- dag2cpdag(g)
  > ## plot(cpdag)
  > cpdag.mat <- t(as(cpdag,"matrix")) ## has correct encoding
  > 
  > ## test 1: using graph, valid
  > g1 <- addBgKnowledge(gInput = cpdag, x = 3, y = 5)
  > m1 <- t(as(g1, "matrix"))
  > res[1] <- (m1[3,5]==0 & m1[5,3]==1) ## should be true
  > 
  > ## test 2: using matrix, valid
  > m2 <- addBgKnowledge(gInput = cpdag.mat, x = 3, y = 5, verbose = TRUE)
  Added orientation 3 -> 5 to the PDAG.
  > res[2] <- (m2[3,5]==0 & m2[5,3]==1) ## should be true
  > 
  > ## test 3: using matrix, invalid
  > m3 <- addBgKnowledge(gInput = cpdag.mat, x = 6, y = 3, verbose = TRUE)
  Invalid bg knowledge! Cannot add orientation 6 -> 3  because there is no edge between 6 and 3 in the PDAG. 
  > res[3] <- is.null(m3)
  > 
  > ## test 4: using graph, invalid
  > g4 <- addBgKnowledge(gInput = cpdag, x = 6, y = 3, verbose = TRUE)
  Invalid bg knowledge! Cannot add orientation 6 -> 3  because there is no edge between 6 and 3 in the PDAG. 
  > res[4] <- is.null(g4)
  > 
  > ## test 5: using matrix, invalid
  > m5 <- addBgKnowledge(gInput = cpdag.mat, x = 1, y = 3, verbose = TRUE)
  Invalid bg knowledge! Cannot add orientation 1 -> 3  because there is no edge between 1 and 3 in the PDAG. 
  > res[5] <- is.null(m5)
  > 
  > ## test 6: using graph, invalid
  > g6 <- addBgKnowledge(gInput = cpdag, x = 1, y = 3, verbose = TRUE)
  Invalid bg knowledge! Cannot add orientation 1 -> 3  because there is no edge between 1 and 3 in the PDAG. 
  > res[6] <- is.null(g6)
  > 
  > ## test 7: empty background knowledge: Meek rule 1
  > m7 <- matrix(0, 3,3)
  > colnames(m7) <- rownames(m7) <- as.character(1:ncol(m7))
  > r7T <- m7
  > m7[2,1] <- m7[3,2] <- m7[2,3] <- 1
  > r7 <- addBgKnowledge(gInput = m7, x = c(), y = c(), verbose = TRUE,
  +                      checkInput = FALSE)
  
  Rule 1: 1 -> 2  and  2 - 3  where  1  and  3  not connected:  2 -> 3 
  > r7T[2,1] <- r7T[3,2] <- 1
  > res[7] <- identical(r7,r7T)
  > 
  > ## test 8: empty background knowledge: Meek rule 2
  > m8 <- matrix(0, 3,3)
  > colnames(m8) <- rownames(m8) <- as.character(1:ncol(m8))
  > r8T <- m8
  > m8[1,2] <- m8[2,3] <- m8[3,1] <- m8[3,2] <- 1
  > r8 <- addBgKnowledge(gInput = m8, x = c(), y = c(), verbose = TRUE,
  +                      checkInput = FALSE)
  
  Rule 2: Kette  2 -> 1 -> 3 : 2 -> 3 
  > r8T[1,2] <- r8T[3,1] <- r8T[3,2] <- 1
  > res[8] <- identical(r8,r8T)
  > 
  > ## test 9: empty background knowledge: Meek rule 3
  > m9 <- matrix(0, 4,4)
  > colnames(m9) <- rownames(m9) <- as.character(1:ncol(m9))
  > r9T <- m9
  > m9[1,2:4] <- m9[2,1] <- m9[3,c(1,2,4)] <- m9[4,1] <- 1
  > r9 <- addBgKnowledge(gInput = m9, x = c(), y = c(), verbose = TRUE,
  +                      checkInput = FALSE)
  
  Rule 3: 1 -> 3 
  > r9T[1,c(2,4)] <- r9T[2,1] <- r9T[3,c(1,2,4)] <- r9T[4,1] <- 1
  > res[9] <- identical(r9,r9T)
  > 
  > ## test 10: empty background knowledge: Meek rule 4
  > m10 <- matrix(0, 4,4)
  > colnames(m10) <- rownames(m10) <- as.character(1:ncol(m10))
  > r10T <- m10
  > m10[1,2:4] <- m10[2,c(1,3)] <- m10[3,c(1,4)] <- m10[4,1] <- 1
  > r10 <- addBgKnowledge(gInput = m10, x = c(), y = c(), verbose = TRUE,
  +                       checkInput = FALSE)
  Rule 4 applied 
  > r10T[1,c(3,4)] <- r10T[2,c(1,3)] <- r10T[3,c(1,4)] <- r10T[4,1] <- 1
  > res[10] <- identical(r9,r9T)
  > 
  > ## final result
  > stopifnot(all(res))
  > 
  > proc.time()
     user  system elapsed 
     1.59    0.25    1.82 
Running the tests in 'tests/test_ages.R' failed.
Complete output:
  > ## Test file ages
  > library(pcalg)
  > suppressWarnings(RNGversion("3.5.0"))
  > (doExtras <- pcalg:::doExtras())
  [1] FALSE
  > 
  > ## Known example where ges and ages output a different result
  > bool3 <- TRUE
  > set.seed(77)
  > 
  > p <- 8
  > n <- 5000
  > ## true DAG:
  > vars <- c("Author", "Bar", "Ctrl", "Goal", paste0("V",5:8))
  > gGtrue <- randomDAG(p, prob = 0.3, V = vars)
  > data = rmvDAG(n, gGtrue)
  > 
  > ## Estimate the aggregated PDAG with ages
  > ages.fit <- ages(data = data)
  > 
  > 
  > ## Estimate the essential graph with ges
  > score <- new("GaussL0penObsScore", data)
  > ges.fit <- ges(score)
  > 
  > diff <- as(ges.fit$essgraph,"matrix") - as(ages.fit$essgraph,"matrix")
  > 
  > bool3 <- ( bool3 & (diff[6,2]==1) * (diff[8,2]==1) *(sum(abs( as(ges.fit$essgraph,"matrix") - as(ages.fit$essgraph,"matrix") ))==2) )
  > stopifnot(bool3)
  > 
  > if (doExtras) {
  + ## Test 1: Need to make sure that the skeleton of ges and ages are the same
  + 
  + bool1 <- TRUE
  + for(i in 1:10){
  +   
  +   p <- 20
  +   n <- 5000
  +   ## true DAG:
  +   gGtrue <- randomDAG(p, prob = 0.4)
  +   data <- rmvDAG(n, gGtrue)
  +   
  +   
  +   score <- new("GaussL0penObsScore", data=data, lambda = 200)
  +   ges.fit <- ges(score = score, phase = c("forward", "backward"), iterate = F)
  +   ages.fit <- ages(data = data, lambda_min = 200)
  +   
  +   
  +   ges.mat <- as(ges.fit$essgraph,"matrix")
  +   ages.mat <- as(ages.fit$essgraph,"matrix")
  +   # 
  +   # ages1mat2 <- ages1[[3]][[length(ages1[[3]])]]
  +   
  +   
  +   bool1 <- ( bool1 & all((((ges.mat + t(ges.mat))!=0)*1)==(((ages.mat + t(ages.mat))!=0)*1)) )
  + }
  + stopifnot(bool1)
  + 
  + ## Test 2: Number of CPDAGs used and number of penalty parameter has to be the same
  + bool2 <- TRUE
  + for(i in 1:10){
  +   
  +   p <- 20
  +   n <- 5000
  +   ## true DAG:
  +   gGtrue <- randomDAG(p, prob = 0.4)
  +   data = rmvDAG(n, gGtrue)
  +   
  +   agesfit <- ages(data = data, lambda_min = 100)
  +   
  +   bool2 <- ( bool2 & (length(agesfit$CPDAGsList) == length(agesfit$lambda)) )
  + }
  + stopifnot(bool2)
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.71    0.23    1.93 
Running the tests in 'tests/test_bicscore.R' failed.
Complete output:
  > #' Tests the calculation of BIC and MLE, as well as basic functions
  > #' the corresponding score class
  > #'
  > #' @author Alain Hauser
  > #' $Id: test_bicscore.R 409 2017-02-05 15:41:51Z alhauser $
  > 
  > cat("Testing the calculation of the BIC score...\n")
  Testing the calculation of the BIC score...
  > 
  > library(pcalg)
  > 
  > load("test_bicscore.rda") # in directory tests/ i.e., typically *not* installed
  > 
  > ## Tolerance for numerical comparison
  > tol <- sqrt(.Machine$double.eps)
  > 
  > ## Define all test settings
  > settings <- expand.grid(format = c("scatter", "raw"), cpp = c(FALSE, TRUE),
  +                         stringsAsFactors = FALSE)
  > nreps <- 5
  > 
  > ## Check data storage and calculation of scores
  > for (m in 1:nrow(settings)) {
  +   cat(sprintf("Setting: storage format = %s, C++ library = %s\n",
  +               settings$format[m], settings$cpp[m]))
  + 
  +   for (i in 1:nreps) {
  +     perm <- 1:nrow(gauss.data)
  + 
  +     ## Randomly permute data
  +     if (i > 1) {
  +       set.seed(i)
  +       perm <- sample(perm)
  +     }
  + 
  +     ## Create the score object with valid data
  +     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],
  +                  intercept = FALSE)
  + 
  +     if (any(score$pp.dat$data.count != 1000)) {
  +       stop("The number of non-interventions are not calculated correctly.")
  +     }
  + 
  +     stopifnot(is.list(score$pp.dat$non.int), length(score$pp.dat$non.int) == 5)
  +     stopifnot(all(sapply(1:5, function(j) {
  +       isTRUE(all.equal(score$pp.dat$non.int[[j]],
  +                        setdiff(1:1200, 200*j + 1:200)))
  +     })))
  + 
  +     if (settings$format[m] == "scatter") {
  +       if (!isTRUE(all.equal(score$pp.dat$scatter.index, 1:5))) {
  +         stop("The indices of the scatter matrices are not calculated correctly.")
  +       }
  + 
  +       for (j in 1:5) {
  +         if (!isTRUE(all.equal(score$pp.dat$scatter[[j]][1:5, 1:5],
  +                               gauss.scatter[[j]],
  +                               tolerance = tol))) {
  +           stop("The scatter matrices are not calculated correctly.")
  +         }
  +       }
  +     } # IF "scatter"
  + 
  +     for (j in 1:5) {
  +       if (!isTRUE(all.equal(gauss.loc.score[[j]],
  +                             score$local.score(j, gauss.parents[[j]]),
  +                             tolerance = tol))) {
  +         stop("The local score is not calculated correctly.")
  +       }
  +     }
  + 
  +     # print(lapply(1:5, function(i) score$local.fit(i, gauss.parents[[i]])))
  + 
  +     for (j in 1:5) {
  +       local.mle <- score$local.fit(j, gauss.parents[[j]])
  +       if (length(local.mle) != length(gauss.mle[[j]]) ||
  +           !isTRUE(all.equal(gauss.mle[[j]], local.mle, tolerance = tol))) {
  +         stop("The local MLE is not calculated correctly.")
  +       }
  +     }
  +   }
  + }
  Setting: storage format = scatter, C++ library = FALSE
  Setting: storage format = raw, C++ library = FALSE
  Setting: storage format = scatter, C++ library = TRUE
  Setting: storage format = raw, C++ library = TRUE
  > 
  > ## List targets in a non-unique way, check if representation is corrected for
  > ## internal storage
  > temp.targets <- gauss.targets
  > temp.targets[[2]] <- rep(temp.targets[[2]], 4)
  > score <- new("GaussL0penIntScore",
  +              targets = temp.targets,
  +              target.index = gauss.target.index,
  +              data = gauss.data,
  +              format = "scatter",
  +              use.cpp = FALSE,
  +              intercept = FALSE)
  > stopifnot(isTRUE(all.equal(score$pp.dat$targets, gauss.targets)))
  > stopifnot(isTRUE(all.equal(score$pp.dat$target.index, gauss.target.index)))
  > 
  > ## Try to create the score object with non-valid data,
  > ## check if error is thrown
  > stopifnot(isTRUE(
  +     tryCatch(
  +         score <- new("GaussL0penIntScore",
  +             targets = gauss.targets,
  +             target.index = gauss.target.index),
  +         error = function(e) {
  +           cat(paste("  Error caught:", e$message, "\n", sep = " "))
  +           TRUE
  +         }
  +     )))
    Error caught: Targets are out of range. 
  > 
  > set.seed(307)
  > temp.targets <- gauss.targets
  > temp.targets <- c(temp.targets, temp.targets[[6]])
  > temp.target.index <- gauss.target.index
  > temp.target.index[sample(which(gauss.target.index == 6), size = 20)] <- length(temp.targets)
  > stopifnot(isTRUE(
  +     tryCatch(
  +         score <- new("GaussL0penIntScore",
  +             targets = temp.targets,
  +             target.index = temp.target.index,
  +             data = gauss.data),
  +         error = function(e) {
  +           cat(paste("  Error caught:", e$message, "\n", sep = " "))
  +           TRUE
  +         }
  +     )))
    Error caught: List of targets must be unique. 
  > 
  > temp.targets <- gauss.targets
  > temp.targets[[2]] <- c(temp.targets[[2]], 9)
  > stopifnot(isTRUE(
  +     tryCatch(
  +         score <- new("GaussL0penIntScore",
  +             targets = temp.targets,
  +             target.index = gauss.target.index,
  +             data = gauss.data),
  +         error = function(e) {
  +           cat(paste("  Error caught:", e$message, "\n", sep = " "))
  +           TRUE
  +         }
  +     )))
    Error caught: Targets are out of range. 
  > 
  > temp.target.index <- gauss.target.index
  > temp.target.index[1] <- length(gauss.targets) + 1
  > stopifnot(isTRUE(
  +     tryCatch(score <- new("GaussL0penIntScore",
  +             targets = gauss.targets,
  +             target.index = temp.target.index,
  +             data = gauss.data),
  +         error = function(e) {
  +           cat(paste("  Error caught:", e$message, "\n", sep = " "))
  +           TRUE
  +         }
  +     )))
    Error caught: Target index is out of range. 
  > 
  > cat("Done.\n")
  Done.
  > 
  > proc.time()
     user  system elapsed 
     1.84    0.29    2.07 
Running the tests in 'tests/test_displayAmat.R' failed.
Complete output:
  > library(pcalg)
  > suppressWarnings(RNGversion("3.5.0"))
  > ##################################################
  > ## pcAlgo object
  > ##################################################
  > ## Load predefined data
  > data(gmG)
  > n <- nrow    (gmG8$x)
  > V <- colnames(gmG8$x)
  > 
  > ## define sufficient statistics
  > suffStat <- list(C = cor(gmG8$x), n = n)
  > ## estimate CPDAG
  > skel.fit <- skeleton(suffStat, indepTest = gaussCItest,
  +                      alpha = 0.01, labels = V)
  > (amSkel <- as(skel.fit, "amat"))
  Adjacency Matrix 'amat' (8 x 8) of type 'cpdag':
         Author Bar Ctrl Goal V5 V6 V7 V8
  Author      .   1    .    .  .  1  .  1
  Bar         1   .    1    .  1  .  .  .
  Ctrl        .   1    .    .  .  .  .  .
  Goal        .   .    .    .  .  .  .  .
  V5          .   1    .    .  .  1  .  1
  V6          1   .    .    .  1  .  1  .
  V7          .   .    .    .  .  1  .  .
  V8          1   .    .    .  1  .  .  .
  > str(amSkel)
   'amat' num [1:8, 1:8] 0 1 0 0 0 1 0 1 1 0 ...
   - attr(*, "dimnames")=List of 2
    ..$ : chr [1:8] "Author" "Bar" "Ctrl" "Goal" ...
    ..$ : chr [1:8] "Author" "Bar" "Ctrl" "Goal" ...
   - attr(*, "type")= chr "cpdag"
  > stopifnot(attr(amSkel, "type") == "cpdag",
  +           amSkel["Author", "Bar"] == 1,
  +           amSkel["Bar", "Author"] == 1,
  +           amSkel["Ctrl","Author"] == 0)
  > 
  > pc.fit <- pc(suffStat, indepTest = gaussCItest,
  +              alpha = 0.01, labels = V)
  > (amPC <- as(pc.fit, "amat"))
  Adjacency Matrix 'amat' (8 x 8) of type 'cpdag':
         Author Bar Ctrl Goal V5 V6 V7 V8
  Author      .   1    .    .  .  .  .  .
  Bar         1   .    1    .  .  .  .  .
  Ctrl        .   1    .    .  .  .  .  .
  Goal        .   .    .    .  .  .  .  .
  V5          .   1    .    .  .  .  .  .
  V6          1   .    .    .  1  .  .  .
  V7          .   .    .    .  .  1  .  .
  V8          1   .    .    .  1  .  .  .
  > stopifnot(attr(amPC, "type") == "cpdag",
  +           amPC["V5", "V8"] == 0,
  +           amPC["V8", "V5"] == 1,
  +           amPC["Goal","Author"] == 0)
  > 
  > ##################################################
  > ## fciAlgo object
  > ##################################################
  > set.seed(42)
  > p <- 7
  > ## generate and draw random DAG :
  > myDAG <- randomDAG(p, prob = 0.4)
  > 
  > ## find PAG using the FCI algorithm
  > myC <- cov2cor(trueCov(myDAG))
  > suffStat <- list(C = myC, n = 10^9)
  > V <- LETTERS[1:p] ## labels of nodes
  > 
  > fmFCI <- fci(suffStat, indepTest=gaussCItest, labels = V,
  +              alpha = 0.9999, doPdsep = FALSE)
  > (amFCI <- as(fmFCI, "amat"))
  Adjacency Matrix 'amat' (7 x 7) of type 'pag':
    A B C D E F G
  A . . 1 . 2 2 2
  B . . . . 2 2 2
  C 1 . . 1 2 2 2
  D . . 1 . . 2 .
  E 1 1 1 . . . 2
  F 1 1 1 1 . . 2
  G 3 3 3 . 1 3 .
  > stopifnot(attr(amFCI, "type") == "pag",
  +           amFCI["B","E"] == 2,
  +           amFCI["C","D"] == 1,
  +           amFCI["G","A"] == 3)
  > 
  > proc.time()
     user  system elapsed 
     1.67    0.17    1.84 
Running the tests in 'tests/test_dsep.R' failed.
Complete output:
  > #### Compare our dsep() with  dSep() from package "ggm" :
  > library(pcalg)
  > 
  > set.seed(22)
  > p <- 8
  > nreps <- 10
  > ok <- rep(FALSE,nreps)
  > for (i in 1:nreps) {
  +   myDAG <- randomDAG(p, prob = 0.3)
  +   amat <- as(myDAG,"matrix")
  +   amat[amat!=0] <- 1
  + 
  +   x <- sample(1:p,1)
  +   y <- sample(setdiff(1:p,x),1)
  +   S <- sample(setdiff(1:p,c(x,y)),sample(1:5,1))
  + 
  +   dsepOld <- ggm::dSep(amat,as.character(x),as.character(y),as.character(S))
  +   dsepRes <- dsep	   (as.character(x),as.character(y),as.character(S),
  + 			    myDAG)
  +   ok[i] <- (dsepRes == dsepOld)
  + }
  > 
  > if (!all(ok)) stop("Test dsep wrong: dsep oracle made a mistake!")
  > 
  > #### Test with a graph that is NOT top. sorted (need not be)
  > amat <- rbind(c(0,1,1,1,0,1,0),
  +               c(0,0,0,0,0,0,0),
  +               c(0,1,0,1,0,0,0),
  +               c(0,0,0,0,0,0,0),
  +               c(1,0,1,0,0,1,0),
  +               c(0,0,0,0,0,0,0),
  +               c(1,0,0,0,1,0,0))
  >                
  > colnames(amat) <- rownames(amat) <- as.character(1:7)
  > g <- as(amat,"graphNEL")
  > 
  > ok <- rep(FALSE, 2)
  > ok[1] <- ( dsep("7","2", c("1","5"),g) == TRUE ) ## should be TRUE
  > ok[2] <- ( dsep("7","2", c("3","5"),g) == FALSE ) ## should be FALSE
  > 
  > if (!all(ok)) stop("Test dsep wrong: dsep oracle made a mistake!")
  > 
  > proc.time()
     user  system elapsed 
     1.60    0.17    1.70 
Running the tests in 'tests/test_fciPlus.R' failed.
Complete output:
  > library(pcalg)
  > 
  > source(system.file(package="Matrix", "test-tools-1.R", mustWork=TRUE))
  Loading required package: tools
  > ##--> showProc.time(), assertError(), relErrV(), ...
  > cat("doExtras:", (doExtras <- pcalg:::doExtras()), "\n")
  doExtras: FALSE 
  > 
  > ##################################################
  > ## simulate graph; remove confounders
  > ##################################################
  > 
  > ##' (utility)
  > ##' for (i in 1:p),
  > ##'   vect[i] := new name of node i in the graph when series nodes are removed
  > conf2vec <- function(confs, p) {
  +   stopifnot(p >= 1) # confs may be *empty*
  +   vect <- sapply(1:p, function(r) r - sum(1 <= confs & confs <= r))
  +   vect[confs] <- NA_integer_
  +   vect
  + }
  > 
  > 
  > ##' random DAG  with confounders removed
  > rDAG_noC <- function(p, prob, max.no.conf = 3,
  +                      V = c(LETTERS,letters)[1:p], verbose=TRUE) {
  +   stopifnot(p == as.integer(p), length(p) >= 1, p >= 2, prob > 0,
  +             !anyDuplicated(V))
  + 
  +   g <- randomDAG(p, prob, V=V)  ## true DAG
  +   ## get the adjacency matrix
  +   A <- as(g, "matrix")
  +   A[A != 0] <- 1
  + 
  +   ## first three confounders should be removed
  +   ## confounders - vector of confounding nodes to be removed
  +   n <- ncol(A) ## number of nodes
  +   confounders <- numeric(max.no.conf)
  +   counter <- 1L
  +   for(k in 1:n) {
  +     if (sum(A[k,]) > 1) { ## k has more than one child
  +       confounders[counter] <- k ## mark node k; leave out later
  +       counter <- counter +1L
  +       if(counter > max.no.conf)
  +         break
  +     }
  +   }
  +   confounders <- confounders[is.c <- confounders > 0]
  +   if(verbose) {
  +     cat("confounders:\n"); print(confounders)
  +     vect <- conf2vec(confounders, p)
  +     cat("-> nodes after removing confounders:\n"); print(vect)
  +   }
  + 
  +   ## define sufficient statistics (d-separation oracle)
  +   trueC <- trueCov(g)
  +   ## delete rows and columns belonging to confounder variables in confounders
  +   C.1 <- if(any(is.c)) trueC[-confounders, -confounders] else trueC
  +   list(g = g, trueC = trueC, trueCov.1 = C.1, confounders = confounders, p. = p - sum(is.c))
  + } ## {rDAG_noC}
  > 
  > chk_rDAG <- function(rD) {
  +   stopifnot(is.list(rD), c("trueC", "trueCov.1", "p.","confounders") %in% names(rD))
  +   p <- nrow(rD$trueC)
  +   p. <- rD$p.
  +   conf <- rD$confounders
  +   no.c <- sum(conf > 0)
  +   stopifnot(is(rD$g, "graph"))
  +   V <- rownames(rD$trueC)
  +   stopifnot(p. == p - no.c, dim(rD$trueCov.1) == c(p., p.),
  +             isSymmetric(rD$trueC),
  +             isSymmetric(rD$trueCov.1),
  +             identical(V, colnames(as(rD$g, "matrix"))),
  +             identical(if(no.c) V[-conf] else V, rownames(rD$trueCov.1)))
  + }
  > 
  > seed <- 547 ## Seed 547: Dsep Link wird benoetigt
  > set.seed(seed)
  > 
  > ## will get slow if p > 14; these settings take already a few minutes
  > if (doExtras) {
  +   p <- 8
  + } else {
  +   p <- 14
  + }
  > r1 <- rDAG_noC(p, prob = 0.2)
  confounders:
  [1] 1 2 4
  -> nodes after removing confounders:
   [1] NA NA  1 NA  2  3  4  5  6  7  8  9 10 11
  > c(p = p, p. = r1$p.)
   p p. 
  14 11 
  > chk_rDAG(r1)
  > ###
  > 
  > indepTest1 <- gaussCItest
  > 
  > ## transform covariance matrix into a correlation matrix
  > p. <- r1$p.
  > suffStat1 <- list(C = cov2cor(r1$trueCov.1), n = 10^9)
  > 
  > ## FCI+: neue Variante
  > showSys.time(
  +   fciplus.fit <- fciPlus(suffStat1, indepTest = indepTest1, alpha = 0.99, p = p.)
  + )
  
   7 *-> 10 <-* 8 
  Sxz= and Szx=  
  
   8 *-> 10 <-* 7 
  Sxz=  and Szx= 
  
   4 *-> 11 <-* 5 
  Sxz= and Szx=  
  
   4 *-> 11 <-* 7 
  Sxz= and Szx=  
  
   5 *-> 11 <-* 4 
  Sxz=  and Szx= 
  
   5 *-> 11 <-* 7 
  Sxz= and Szx=  
  
   7 *-> 11 <-* 4 
  Sxz=  and Szx= 
  
   7 *-> 11 <-* 5 
  Sxz=  and Szx= 
  Time    user  system elapsed 
  Time    0.08    0.00    0.08 
  > 
  > ## Fit FCI
  > showSys.time(
  +   fci.fit <- fci(suffStat1, indepTest = indepTest1, alpha = 0.99, p = p.,
  +                  rules = rep(TRUE, 10), verbose = FALSE)
  + )
  Time    user  system elapsed 
  Time    0.08    0.00    0.08 
  > 
  > print.table(fciplus.fit@amat, zero.print = ".")
     1 2 3 4 5 6 7 8 9 10 11
  1  . . . . 1 . . . .  .  .
  2  . . 1 1 . 1 . . 1  .  .
  3  . 1 . 1 . 1 . . 1  .  .
  4  . 1 1 . . 1 . . 1  .  2
  5  1 . . . . . . 1 .  .  2
  6  . 1 1 1 . . . . 1  .  .
  7  . . . . . . . . .  2  2
  8  . . . . 1 . . . .  2  .
  9  . 1 1 1 . 1 . . .  .  .
  10 . . . . . . 1 1 .  .  .
  11 . . . 1 1 . 1 . .  .  .
  > ## MM{FIXME}: want to see node correct names !!
  > 
  > if(!all(fciplus.fit@amat == fci.fit@amat))
  +     stop("Test fciPlus() differing from fci(): PAG not found correctly?")
  > showProc.time()
  Time (user system elapsed): 0.45 0 0.45 
  > 
  > if(doExtras) { ## more tests
  +  for(p in 2:20) { ##     p = 10
  +    cat("\n\np = ", p, ":\n~~~~~~\n")
  +    set.seed(p)
  +    for(pr in c(0.75, 1.5) / p) {
  +      seed <- 10*p + ceiling(2*p*pr)
  +      cat("\nprob = ", pr, ", seed = ", seed,":\n")
  +      set.seed(seed)
  +      rr <- rDAG_noC(p, prob = pr) # quite sparse
  +      chk_rDAG(rr)
  + 
  +      ## transform covariance matrix into a correlation matrix
  +      p. <- rr$p.
  +      suffStat1 <- list(C = cov2cor(rr$trueCov.1), n = 10^9)
  + 
  +      ## FCI+: neue Variante
  +      showSys.time(
  +        fciplus.fit <- fciPlus(suffStat1, indepTest = indepTest1, alpha = 0.99, p = p.)
  +      )
  + 
  +      ## Fit FCI
  +      showSys.time(
  +        fci.fit <- fci(suffStat1, indepTest = indepTest1, alpha = 0.99, p = p.,
  +                       rules = rep(TRUE, 10), verbose = FALSE)
  +      )
  + 
  +      cat("Adjacency matrix:\n")
  +      print.table(fciplus.fit@amat, zero.print = ".")
  +      ## MM{FIXME}: want to see node correct names !!
  + 
  +      if(!all(fciplus.fit@amat == fci.fit@amat))
  +        stop("Test fciPlus() differing from fci(): PAG not found correctly?")
  +      else cat("=================================================\n")
  +    }# for different 'pr'
  + 
  +  }# for a set of 'p'
  + 
  + } ## if (doExtras)
  > 
  > proc.time()
     user  system elapsed 
     1.79    0.15    1.95 
Running the tests in 'tests/test_gSquareBin.R' failed.
Complete output:
  > library(pcalg)
  > 
  > cat("doExtras:", (doExtras <- pcalg:::doExtras()), "\n")
  doExtras: FALSE 
  > 
  > ##################################################
  > ## Simple 2 Variables  x1 -> x2
  > ##################################################
  > set.seed(123)
  > 
  > ## only 2 Variables
  > b0 <- 0
  > b1 <- 1
  > n <- 1000
  > x1 <- rbinom(n, 1, p = 1/2) ## == sample(c(0,1), n, replace=TRUE)
  > x2 <- rbinom(n, 1, p <- plogis(b0 + b1* x1))
  > 
  > res1 <- fisher.test(x1,x2)$p.value         # 6.14 e-21
  > res2 <- gSquareBin(1,2,NULL, cbind(x1,x2)) # 3.74 e-21
  > 
  > (eq1 <- all.equal(c(res1,     res2),
  +                   c(6.142e-21, 3.7401e-21), tol = 1e-5))
  [1] TRUE
  > ok1 <- isTRUE(eq1)
  > 
  > ##################################################
  > ## chain of 3 variables: x1 -> x2 -> x3
  > ##################################################
  > b0 <- 0
  > b1 <- 1
  > b2 <- 1
  > n <- 10000
  > 
  > set.seed(12)
  > x1 <- rbinom(n, 1, p=1/2)
  > x2 <- rbinom(n, 1, p2 <- plogis(b0 + b1* x1))
  > x3 <- rbinom(n, 1, p3 <- plogis(b0 + b2* x2))
  > ##
  > dat <- cbind(x1,x2,x3)
  > ##
  > system.time(
  + pv.2 <- c(gSquareBin(3,1, 2, dat),
  +           gSquareBin(1,3, 2, dat),
  +           gSquareBin(1,2, 3, dat),
  +           gSquareBin(3,2, 1, dat))
  + )
     user  system elapsed 
     0.01    0.00    0.02 
  > 
  > (eq2 <- c(all.equal(pv.2, c(rep(0.89393365, 2), 0, 0), tol=1e-8),
  +           all.equal(pv.2[1], pv.2[2], tol=1e-6)))
  [1] TRUE TRUE
  > 
  > ok2 <- is.logical(eq2) && all(eq2)
  > 
  > ok3 <- ok4 <- ok5 <- ok6 <- TRUE
  > if (doExtras) {
  +     suppressWarnings(RNGversion("3.5.0"))
  + ##################################################
  + ## collider: x1 -> x3 <- x2
  + ##################################################
  + b0 <- 0
  + b1 <- 1.3
  + b2 <- 1.7
  + n <- 10000
  + set.seed(13)
  + x1 <- rbinom(n, 1, p=1/2)
  + x2 <- rbinom(n, 1, p=1/2)
  + x3 <- rbinom(n, 1, p3 <- plogis(b0 + b1*x1 + b2*x2))
  + ##
  + dat <- cbind(x1,x2,x3)
  + ##
  + system.time(
  + pv3 <- c(gSquareBin(1,2, 3,  dat),
  +          gSquareBin(1,2,NULL,dat),
  +          gSquareBin(2,1, 3,  dat),
  +          gSquareBin(1,3, 2,  dat))
  + )
  + (eq3 <- all.equal(pv3, c(0, 0.9054922, 0,0), tol=1e-7))
  + ok3 <- isTRUE(eq3)
  + 
  + 
  + ##################################################
  + ## add noise variables to collider model
  + ##################################################
  + x4 <- sample(c(0,1),n, replace=TRUE)
  + x5 <- sample(c(0,1),n, replace=TRUE)
  + x6 <- sample(c(0,1),n, replace=TRUE)
  + 
  + dat <- cbind(x1,x2,x3,x4,x5,x6)
  + 
  + system.time(
  + pv4 <- c(gSquareBin(1,2, c(3,4,5,6), dat),
  +          gSquareBin(1,3, c(2,4,5,6), dat),
  +          gSquareBin(2,3, c(1,4,5,6), dat),
  +          gSquareBin(2,4, c(1,3,5,6), dat),
  +          gSquareBin(4,2, c(1,3,5,6), dat), # same
  +          gSquareBin(3,4, c(1,2,5,6), dat),
  +          ##
  +          gSquareBin(4,2, c(3,5,6),   dat))
  + )
  + 
  + (eq4 <- c(all.equal(pv4, c(0, 0, 0, 0.0269101, 0.0269101, 0.7260011, 0.1121826),
  +                    tol = 2e-7),
  +           all.equal(pv4[4], pv4[5], tol=1e-6)))
  + 
  + ok4 <- is.logical(eq4) && all(eq4)
  + 
  + 
  + ##################################################
  + ## rectangle model
  + ##################################################
  + b0 <- 0
  + b1 <- 1
  + b2 <- 1
  + b3 <- 0.5
  + b4 <- 0.5
  + n <- 10000
  + set.seed(31)
  + x1 <- rbinom(n, size=1, prob = 1/2)# = sample(c(0,1), n, replace=TRUE)
  + x2 <- rbinom(n, size=1, prob = (p2 <- plogis(b0 + b1* x1)))
  + x3 <- rbinom(n, size=1, prob = (p3 <- plogis(b0 + b2* x1)))
  + x4 <- rbinom(n, size=1, prob = (p4 <- plogis(b0 + b3* x2 + b4* x3)))
  + ##
  + dat <- cbind(x1,x2,x3,x4)
  + 
  + system.time(
  + pv5 <- c(gSquareBin(1,4,c(2,3),dat),
  +          gSquareBin(1,4,  2,   dat),
  +          gSquareBin(1,4,  3,   dat),
  +          gSquareBin(1,4,NULL,  dat),
  +          gSquareBin(2,3,  4,   dat),
  +          gSquareBin(2,3,  1,   dat),
  +          gSquareBin(2,3,c(1,4),dat))
  + )
  + 
  + (eq5 <- all.equal(pv5, c(0.38484, 0.13601, 0.10974, 0, 0.00023, 0.56828, 0.43436), 4e-4))
  + ok5 <- isTRUE(eq5)
  + 
  + ##################################################
  + ## rectangle model -- extended
  + ##################################################
  + x5 <- rbinom(n, size=1, prob = (p5 <- plogis(b0 + x2)))
  + x6 <- rbinom(n, size=1, prob = (p6 <- plogis(b0 + x3/2 - x4)))
  + x7 <- rbinom(n, size=1, prob = (p7 <- plogis(b0 + x5 - x6)))
  + x8 <- rbinom(n, size=1, prob = (p8 <- plogis(b0 - x3 + x6/2)))
  + x9 <- rbinom(n, size=1, prob = (p9 <- plogis(b0 + x2 - x7)))
  + ##
  + dat <- cbind(x1,x2,x3,x4,x5,x6,x7,x8,x9)
  + 
  + ## do many tests:
  + p <- ncol(dat) # = 9
  + ## all (x,y) pairs:
  + xy <- combn(p,2)
  + "%w/o%" <- function(x, y) x[!x %in% y] #--  x without y
  + ip <- 1:p
  + str(rest <- apply(xy, 2, function(e) ip %w/o% e)) # 7 x 36
  + set.seed(101)
  + in2 <- seq_len(ncol(xy))
  + args <- c(lapply(in2, function(i) list(x=xy[1,i], y=xy[2,i], S=rest[,i])),
  +           lapply(in2, function(i) list(x=xy[1,i], y=xy[2,i], S=sample(rest[,i], size= p-2-1))),
  +           lapply(in2, function(i) list(x=xy[1,i], y=xy[2,i], S=sample(rest[,i], size= sample(0:(p-2), 1)))),
  +           lapply(in2, function(i) list(x=xy[1,i], y=xy[2,i], S=sample(rest[,i], size= sample(2:(p-2), 1)))))
  + 
  + ia <- TRUE
  + 
  + system.time(
  + pv6 <- vapply(args, function(L) do.call(gSquareBin, c(L, list(dm=dat))), 1)
  + )# if(doExtras) about one minute
  + 
  + if(FALSE)
  +     dput(signif(zapsmall(pv6), 4))
  + pv6.ok <- c(0, 0, 0.6427, 0.4745, 0.9698, 0.4586, 0.5739, 0.6701, 0.1889,
  +             0, 0, 0.6476, 0.00601, 0.3586, 0, 0, 0.08916, 0, 0.7642, 0, 0.04021,
  +             0.4694, 0, 0.8799, 0.3966, 0.358, 0.01514, 0, 0.3615, 0.3244,
  +             0, 0, 0.6647, 0.6839, 0, 0.4136, 0, 0, 0.7023, 0.303, 0.9189,
  +             0.5501, 0.7896, 0.8239, 0.04095, 0, 0, 0.5531, 0.6374, 0.755,
  +             0, 0, 0.629, 0, 0.8377, 0, 0.5749, 0.6346, 0, 0.9513, 0.3138,
  +             0.6293, 0.000153, 0, 0.7855, 0.1503, 0, 0, 0.2033, 0.9071, 0,
  +             0.2587, 0, 0, 1.08e-05, 0.0002067, 0.03795, 0.5501, 3.5e-06,
  +             0.001827, 0.5683, 0, 0, 0.005468, 0.4085, 0.1476, 0, 0, 0.6822,
  +             0, 0.7642, 0, 0.7184, 0.4596, 0, 0.8799, 0.02051, 0.03413, 0.01689,
  +             0, 0.2786, 0.07084, 0, 0, 5.17e-05, 0.4915, 0, 0.5029, 0, 0,
  +             0.1891, 0.09992, 0.9698, 0.3937, 5.94e-05, 0.5473, 0.3894, 0,
  +             0, 0.2236, 2.3e-06, 0.784, 0, 0, 0.3504, 0, 0.7642, 0, 0.3572,
  +             0.4259, 0, 0.8139, 0.3326, 0.01518, 0.01514, 0, 0.3541, 0.021,
  +             0, 0, 0.6647, 0.6167, 0, 0.4136)
  + 
  + (eq6 <- all.equal(pv6, pv6.ok[ia], tol= 1e-4))
  + ok6 <- isTRUE(eq6)
  + 
  + }
  > if (!all(ok1,ok2,ok3,ok4,ok5,ok6))
  +     stop("Test gSquareBin wrong: Some dependence was not estimated correctly!")
  > 
  > proc.time()
     user  system elapsed 
     1.34    0.25    1.57 
Running the tests in 'tests/test_gSquareDis.R' failed.
Complete output:
  > library(pcalg)
  > 
  > cat("doExtras:", (doExtras <- pcalg:::doExtras()), "\n")
  doExtras: 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')
  +     }
  + })
  > options(warn = 1)# show warnings immediately
  > 
  > load("discr100k.rda")# in directory tests/ i.e., typically *not* installed
  > ## Note: --> ../man/gmD.Rd  is _same_ data, only just the first 10'000 obs.
  > str(dat)# 100'000 x 5 data frame of integers
  'data.frame':	100000 obs. of  5 variables:
   $ X1: int  2 2 1 1 1 2 2 0 2 0 ...
   $ X2: int  0 1 0 0 0 1 1 0 0 0 ...
   $ X3: int  0 1 1 2 0 1 0 1 2 0 ...
   $ X4: int  2 2 3 2 2 0 2 2 2 0 ...
   $ X5: int  1 1 0 1 1 1 1 1 1 0 ...
  > sapply(dat, table)
  $X1
  
      0     1     2 
  19771 30138 50091 
  
  $X2
  
      0     1 
  79875 20125 
  
  $X3
  
      0     1     2 
  27282 47715 25003 
  
  $X4
  
      0     1     2     3 
  24899 10715 47263 17123 
  
  $X5
  
      0     1 
  20120 79880 
  
  > nlev.dat <- sapply(dat, function(v) nlevels(factor(v)))
  > suffStat <- list(dm = dat, nlev = nlev.dat, adaptDF = TRUE)
  > dat. <- data.matrix(dat)
  > 
  > showProc.time()
  Time elapsed:  0.09 0.03 0.12 
  > 
  > ## Check for reasonable error message:
  > Sys.setlocale("LC_MESSAGES", "C")
  Warning in Sys.setlocale("LC_MESSAGES", "C") :
    LC_MESSAGES does not exist on Windows
  [1] ""
  > (er <- tryCatch(gSquareDis(1,3, S=20, dm=dat.),  error=conditionMessage))
  [1] "x, y, and S must all be in {1,..,p}, p=5"
  > ##  "x, y, and S must all be in {1,..,p}, p=5"
  > stopifnot(grepl("must all be in {1,..,p}", er, fixed=TRUE))
  > 
  > 
  > ## Tests for   |S| = 0  and  |S| = 1
  > 
  > pv.5 <- c(
  +     ## Collider
  +     ##        x y  S
  +     disCItest(1,2,NULL,suffStat),# 1 : Must be !=0
  +     disCItest(1,2, 3,  suffStat),# 2 : 0
  +     disCItest(5,2, 4,  suffStat),# 3 : 0
  +     disCItest(5,2,NULL,suffStat),# 4 : != 0
  +     ## Non-collider
  +     disCItest(4,3, 2,  suffStat),# 5 : !=0
  +     disCItest(4,3,NULL,suffStat),# 6 : 0
  +     ## Neighbours
  +     disCItest(4,2,NULL,suffStat),# 7 : 0
  +     disCItest(1,3,NULL,suffStat),# 8 : 0
  +     disCItest(1,5,NULL,suffStat))# 9 : !=0
  > 
  > stopifnot(all.equal(pv.5,
  +                     c(0.5150866,
  +                       0,
  +                       0,
  +                       0.289990516,
  +                       ## Non-collider
  +                       0.818905483,
  +                       0,
  +                       # Neighbours
  +                       0,
  +                       0,
  +                       0.774158263), tol = 4e-9))
  > 
  > showProc.time()
  Time elapsed:  0.39 0.08 0.47 
  > 
  > if(!doExtras) ## smaller data set
  +     dat. <- dat.[1:4096,]
  > 
  > ##  |S| = 2 :
  > pv2.5 <- c(
  +     gSquareDis(3,4, S = c(1,2), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,4, S = c(1,3), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,3, S = c(1,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,3, S = c(1,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,4, S = c(2,3), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,3, S = c(2,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,3, S = c(2,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,2, S = c(3,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,2, S = c(3,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,2, S = c(4,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  + 
  +     gSquareDis(3,5, S = c(1,2), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,5, S = c(1,3), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,5, S = c(1,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,4, S = c(1,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,5, S = c(2,3), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,5, S = c(2,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,4, S = c(2,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,5, S = c(3,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,4, S = c(3,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,3, S = c(4,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE))
  Edge  3 -- 4  with subset S = 1 2 
  Edge  2 -- 4  with subset S = 1 3 
  Edge  2 -- 3  with subset S = 1 4 
  Edge  2 -- 3  with subset S = 1 5 
  Edge  1 -- 4  with subset S = 2 3 
  Edge  1 -- 3  with subset S = 2 4 
  Edge  1 -- 3  with subset S = 2 5 
  Edge  1 -- 2  with subset S = 3 4 
  Edge  1 -- 2  with subset S = 3 5 
  Edge  1 -- 2  with subset S = 4 5 
  Edge  3 -- 5  with subset S = 1 2 
  Edge  2 -- 5  with subset S = 1 3 
  Edge  2 -- 5  with subset S = 1 4 
  Edge  2 -- 4  with subset S = 1 5 
  Edge  1 -- 5  with subset S = 2 3 
  Edge  1 -- 5  with subset S = 2 4 
  Edge  1 -- 4  with subset S = 2 5 
  Edge  1 -- 5  with subset S = 3 4 
  Edge  1 -- 4  with subset S = 3 5 
  Edge  1 -- 3  with subset S = 4 5 
  > 
  > showProc.time()
  Time elapsed:  0.06 0.07 0.13 
  > 
  > pv2.5.EXP <- {
  +     if(doExtras)
  +         c(0.98726, 0,       0, 0, 0.77097, 0,       0,       0, 0, 0.98974,
  +           0.95964, 0.99966, 0, 0, 0.99685, 0.99999, 0.42513, 1, 0.11179, 0)
  +     else
  +         c(0.73219, 0, 0.0003349, 0, 0.90758, 0, 0, 0.98106, 0.46011, 1,
  +           0.93245, 0.99998, 8e-07, 0, 0.99797, 0.99933, 0.86036, 0.99988, 0.84859, 0)
  + }
  > 
  > 
  > (eq52 <- all.equal(pv2.5, pv2.5.EXP, tol = 1e-5))
  [1] TRUE
  > if(!isTRUE(eq52)) stop("gSquareDis(|S| = 2) not ok:", eq52)
  > 
  > 
  > ##  |S| = 3 :
  > pv3.5 <- c(
  +     gSquareDis(1,2, S = c(3,4,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,3, S = c(2,4,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,4, S = c(2,3,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(1,5, S = c(2,3,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,3, S = c(1,4,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,4, S = c(1,3,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(2,5, S = c(1,3,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(3,4, S = c(1,2,5), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(3,5, S = c(1,2,4), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE),
  +     gSquareDis(4,5, S = c(1,2,3), dm = dat., nlev = nlev.dat, adaptDF = TRUE, verbose=TRUE))
  Edge  1 -- 2  with subset S = 3 4 5 
  Edge  1 -- 3  with subset S = 2 4 5 
  Edge  1 -- 4  with subset S = 2 3 5 
  Edge  1 -- 5  with subset S = 2 3 4 
  Edge  2 -- 3  with subset S = 1 4 5 
  Edge  2 -- 4  with subset S = 1 3 5 
  Edge  2 -- 5  with subset S = 1 3 4 
  Edge  3 -- 4  with subset S = 1 2 5 
  Edge  3 -- 5  with subset S = 1 2 4 
  Edge  4 -- 5  with subset S = 1 2 3 
  > showProc.time()
  Time elapsed:  0.08 0.09 0.17 
  > 
  > pv3.5.EXP <- {
  +     if(doExtras)
  +         c(0,       0, 0.95861, 1, 0,       0, 0,       0.99204, 1, 0)
  +     else
  +         c(0.99999, 0, 0.84316, 1, 0.15977, 0, 0.83432, 0.66329, 1, 0)
  + }
  > 
  > (eq53 <- all.equal(pv3.5, pv3.5.EXP, tol=6e-6))
  [1] TRUE
  > if(!isTRUE(eq53)) stop("gSquareDis(|S| = 3) not ok:", eq53)
  > 
  > ###  Extend data, and do (few) tests with  |S| in {4, 5, 6} !!
  > ###  ---------------------------------------------------------
  > 
  > ## Idee:  1) marginal-table fuer die no-incoming (aka 'root') notes
  > ##        2) bivariate table -->  X_j | X_i  for all edges  i -> j
  > ## ==> randomly generate generate the discrete data from these.
  > 
  > ###--> Look at  gR  CRAN task view -- package 'gRain'
  > ###--- can easily construct such objects
  > 
  > ##' Generate discrete (integer, in 0:k) random variable, given one parent (in a DAG)
  > ##'
  > ##' @title Generate Discrete Random Variable Given one parent
  > ##' @param pa.var
  > ##' @param Pr Markov / conditional probabability matrix Pr[i,j] = Prob{newX = i, pa.X = j}
  > ##' @return a discrete random variable \code{newX} of the same length as \code{pa.var}
  > ##' @author Martin Maechler
  > rDiscr1 <- function(pa.var, Pr, zero.based = TRUE) {
  +     stopifnot(pa.var %% 1 == 0, # <- fast test for integer valued
  +               is.matrix(Pr), is.numeric(Pr),
  +               (nlev <- ncol(Pr)) >= 1, ## #{levels} of the (discrete) parent variable
  +               nlev == nlevels(pa.f <- factor(pa.var)))
  + 
  +     pa.n <- as.integer(pa.f)## values in {1, 2,.., nlev}
  +     X <- integer(n <- length(pa.n))
  +     FF <-
  +         if(zero.based)
  +              function(.) which(. == 1L) - 1L ## {1,2,3} - 1  ==> {0,1,2}
  +         else function(.) which(. == 1L)
  +     for(j in seq_len(nlev)) {
  +         is.j <- pa.n == j
  +         nj <- sum(is.j)
  +         X[is.j] <- apply(rmultinom(nj, size = 1, prob = Pr[, j]),
  +                          2L, FF)
  +     }
  +     X
  + }
  > 
  > nrow(dat.) # 4096 or 100'000  depending on  doExtras
  [1] 4096
  > 
  > ## Currently, already |S| = 4 is so slow, we need
  > ## smaller data set:
  > dat. <- dat.[ seq_len((if(doExtras) 8192 else 1024)), ]
  > (n <- nrow(dat.))
  [1] 1024
  > 
  > set.seed(102)
  > P6.5 <- cbind(c(1,5,2),
  +               c(5,2,1))/8 ## P6.[i,j] :=  P[X6 = i | X5 = j]:
  > X6 <- rDiscr1(dat.[,"X5"], Pr = P6.5)
  > P7.3 <- cbind(c(1,5,2),
  +               c(5,2,1),
  +               c(1,4,3) )/8
  > X7 <- rDiscr1(dat.[,"X3"], Pr = P7.3)
  > X8 <- rDiscr1(       X7  , Pr = P7.3)
  > X9 <- apply(rmultinom(n, size=1, prob = c(1,2,3,4,6)/16),
  +             2L, function(.) which(. == 1L) - 1L) # independent
  > P10.9 <- cbind(c(1,7),
  +                c(6,2),
  +                c(3,5),
  +                c(4,4),
  +                c(2,6))/8
  > X10 <- rDiscr1(      X9,  Pr = P10.9)
  > X11 <- rDiscr1(      X9,  Pr = P10.9)
  > dat10 <- cbind(dat., X6, X7, X8, X9, X10, X11)
  > (nlev10 <- sapply(as.data.frame(dat10), function(v) nlevels(factor(v))))
   X1  X2  X3  X4  X5  X6  X7  X8  X9 X10 X11 
    3   2   3   4   2   3   3   3   5   2   2 
  > showProc.time()
  Time elapsed:  0.08 0 0.08 
  > 
  > ## Testing only (not good for sound inference!):
  > ## setting  n.min = n-1  ==> the tests are also run  for "too small" sample sizes
  > ##          ===========
  > 
  > gSquareDis.. <- function(x,y, S, dm, nlev, verbose = TRUE)
  +     gSquareDis(x,y, S=S, dm=dm, nlev=nlev,
  +                n.min=nrow(dm)-1L, adaptDF=TRUE, verbose=verbose)
  > 
  > ##  |S| = 4 :
  > pv4.10 <- c(
  +     gSquareDis..(1,2, S = 3:6     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,3, S = 4:7     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,4, S = 5:8     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,5, S = 6:9     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,3, S = 4:7     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,4, S = c(3,5:7), dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,5, S = c(4,6:8), dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,5, S = c(4,7:9), dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,5, S = 6:9     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,4, S = 5:8     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,4, S = 6:9     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,5, S = 6:9     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(4,5, S = 6:9     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(4,5, S = c(1:3,6), dm= dat10, nlev=nlev10),
  +     gSquareDis..(4,5, S = c(1:3,7), dm= dat10, nlev=nlev10),
  +     gSquareDis..(4,5, S = c(1:3,8), dm= dat10, nlev=nlev10))
  Edge  1 -- 2  with subset S = 3 4 5 6 
  Edge  1 -- 3  with subset S = 4 5 6 7 
  Edge  1 -- 4  with subset S = 5 6 7 8 
  Edge  1 -- 5  with subset S = 6 7 8 9 
  Edge  2 -- 3  with subset S = 4 5 6 7 
  Edge  2 -- 4  with subset S = 3 5 6 7 
  Edge  2 -- 5  with subset S = 4 6 7 8 
  Edge  2 -- 5  with subset S = 4 7 8 9 
  Edge  2 -- 5  with subset S = 6 7 8 9 
  Edge  3 -- 4  with subset S = 5 6 7 8 
  Edge  3 -- 4  with subset S = 6 7 8 9 
  Edge  3 -- 5  with subset S = 6 7 8 9 
  Edge  4 -- 5  with subset S = 6 7 8 9 
  Edge  4 -- 5  with subset S = 1 2 3 6 
  Edge  4 -- 5  with subset S = 1 2 3 7 
  Edge  4 -- 5  with subset S = 1 2 3 8 
  > showProc.time()
  Time elapsed:  0.19 0.06 0.25 
  > dput(signif(zapsmall(pv4.10), 4))
  c(0.9999, 0, 1e-07, 1, 1, 0.0001146, 1, 1, 1, 5.9e-06, 0, 1, 
  0.2205, 0, 0, 0)
  > 
  > pv4.10.EXP <- {
  +   if(doExtras)
  +     c(1,      0, 0.9757, 1, 0.9981, 0,    1,1,1, 0.6949, 0.7201, 1,   0,    0,0,0)
  +   else
  +     c(0.9999, 0, 1e-07,  1, 1, 0.0001146, 1,1,1, 5.9e-06,  0,    1, 0.2205, 0,0,0)
  + }
  > 
  > (eq4.10 <- all.equal(pv4.10, pv4.10.EXP, tol=5e-5))
  [1] TRUE
  > if(!isTRUE(eq4.10)) stop("gSquareDis(|S| = 4) not ok:", eq4.10)
  > 
  > ##  |S| = 5 : --- verbose=2 gives many messages
  > ##  --------      Adding a new combination of parents at sample  <n>"
  > pv5.10 <- c(
  +     gSquareDis..(1,2, S = 3:7     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,3, S = 4:8     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,4, S = 5:9     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,5, S = 6:10    , dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,3, S = 4:8     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,4, S = c(3,5:8), dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,5, S = c(4,6:9), dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,4, S = 5:9     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,4, S = 6:10    , dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,5, S = 5:9     , dm= dat10, nlev=nlev10),
  +     gSquareDis..(4,5, S = 6:10    , dm= dat10, nlev=nlev10))
  Edge  1 -- 2  with subset S = 3 4 5 6 7 
   added a total of 3 new parents
  Edge  1 -- 3  with subset S = 4 5 6 7 8 
   added a total of 3 new parents
  Edge  1 -- 4  with subset S = 5 6 7 8 9 
   added a total of 3 new parents
  Edge  1 -- 5  with subset S = 6 7 8 9 10 
   added a total of 3 new parents
  Edge  2 -- 3  with subset S = 4 5 6 7 8 
   added a total of 2 new parents
  Edge  2 -- 4  with subset S = 3 5 6 7 8 
   added a total of 2 new parents
  Edge  2 -- 5  with subset S = 4 6 7 8 9 
   added a total of 2 new parents
  Edge  3 -- 4  with subset S = 5 6 7 8 9 
   added a total of 3 new parents
  Edge  3 -- 4  with subset S = 6 7 8 9 10 
   added a total of 3 new parents
  Edge  3 -- 5  with subset S = 5 6 7 8 9 
   added a total of 3 new parents
  Edge  4 -- 5  with subset S = 6 7 8 9 10 
   added a total of 4 new parents
  > showProc.time()
  Time elapsed:  0.23 0.02 0.25 
  > dput(signif(zapsmall(pv5.10), 5))
  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
  > 
  > pv5.10.EXP <- {
  +   if(doExtras)
  +     c(1, 0,      0.037821, 1, 1, 0,       1, 0.0005184, 0.0013204, 1, 0)
  +   else
  +     c(1, 0.10938, 0.94484, 1, 1, 0.79863, 1, 0.99495,   0.68719,   1, 1)
  + }
  > 
  > (eq5.10 <- all.equal(pv5.10, pv5.10.EXP, tol= 5e5))
  [1] TRUE
  > if(!isTRUE(eq5.10)) stop("gSquareDis(|S| = 5) not ok:", eq5.10)
  > 
  > ##  |S| = 6 : --- (see |S|=5):
  > 
  > pv6.10 <- c(
  +     gSquareDis..(1,2, S = 3:8       , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,3, S = 4:9       , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,4, S = 5:10      , dm= dat10, nlev=nlev10),
  +     gSquareDis..(1,5, S = c(2:4,6:9), dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,3, S = 4:9       , dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,4, S = 5:10      , dm= dat10, nlev=nlev10),
  +     gSquareDis..(2,5, S = c(4,6:9)  , dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,4, S = 5:10      , dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,4, S = c(1:2,5:8), dm= dat10, nlev=nlev10),
  +     gSquareDis..(3,5, S = c(1:2,6:9), dm= dat10, nlev=nlev10),
  +     gSquareDis..(4,5, S = c(1:3,7:9), dm= dat10, nlev=nlev10))
  Edge  1 -- 2  with subset S = 3 4 5 6 7 8 
   added a total of 3 new parents
  Edge  1 -- 3  with subset S = 4 5 6 7 8 9 
   added a total of 3 new parents
  Edge  1 -- 4  with subset S = 5 6 7 8 9 10 
   added a total of 3 new parents
  Edge  1 -- 5  with subset S = 2 3 4 6 7 8 9 
   added a total of 3 new parents
  Edge  2 -- 3  with subset S = 4 5 6 7 8 9 
   added a total of 2 new parents
  Edge  2 -- 4  with subset S = 5 6 7 8 9 10 
   added a total of 2 new parents
  Edge  2 -- 5  with subset S = 4 6 7 8 9 
   added a total of 2 new parents
  Edge  3 -- 4  with subset S = 5 6 7 8 9 10 
   added a total of 3 new parents
  Edge  3 -- 4  with subset S = 1 2 5 6 7 8 
   added a total of 3 new parents
  Edge  3 -- 5  with subset S = 1 2 6 7 8 9 
   added a total of 3 new parents
  Edge  4 -- 5  with subset S = 1 2 3 7 8 9 
   added a total of 4 new parents
  > showProc.time()
  Time elapsed:  0.24 0.01 0.25 
  > dput(signif(zapsmall(pv6.10), 5))
  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
  > 
  > pv6.10.EXP <- {
  +   if(doExtras)
  +     c(1, 0.023908, 0.000208, 1, 1, 0.94747, 1, 0.0002593, 5.1e-06, 1, 0)
  +   else
  +     c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
  + }
  > 
  > (eq6.10 <- all.equal(pv6.10, pv6.10.EXP, tol= 5e5))
  [1] TRUE
  > if(!isTRUE(eq6.10)) stop("gSquareDis(|S| = 5) not ok:", eq6.10)
  > 
  > proc.time()
     user  system elapsed 
     2.50    0.68    3.14 
Running the tests in 'tests/test_ida.R' failed.
Complete output:
  > library(pcalg)
  > (doExtras <- pcalg:::doExtras())
  [1] FALSE
  > 
  > suppressWarnings(RNGversion("3.5.0"))
  > set.seed(123)
  > nreps <- 100
  > res.local <- logical(nreps)
  > res.opt <- logical(nreps)
  > all.eff.true <- res.local
  > Rnd <- function(e) round(e, 14)## get 14 digits accuracy, as we use true (DAG, cov)
  > Rnd7 <- function(e) round(e, 7)## get 14 digits accuracy, as we use true (DAG, cov)
  > for (i in 1:nreps) {
  +   p <- 2 + rpois(1, lambda = 8) # ==>  p >= 2, E[p] = 10
  +   ## generate and draw random DAG :
  +   myDAG <- randomDAG(p, prob = 0.2)
  +   myCPDAG <- dag2cpdag(myDAG)
  +   mcov <- trueCov(myDAG)
  +   
  +   ## x != y  in {1,2,...p} ;
  +   xy <- sample.int(p, 2); x <- xy[1]; y <- xy[2]
  +   
  +   ## plot(myCPDAG)
  +   eff.true <- Rnd(causalEffect(myDAG, y, x))
  +   all.eff.true[i] <- eff.true
  +   ## cat("x=",x," y=",y," eff=",eff.true,"\n")
  +   
  +   eff.est.local <- Rnd(ida(x,y, mcov, myCPDAG, method="local"))
  +   eff.est.opt <- Rnd(ida(x,y, mcov, myCPDAG, method="optimal"))
  +   res.local[i] <- (eff.true %in% eff.est.local)
  +   res.opt[i] <- (eff.true %in% eff.est.opt)
  + }
  > ## cat('Time elapsed: ', (.pt <- proc.time()),"\n")
  > 
  > ## stem(all.eff.true)
  > if (!all(res.local)) stop("Test ida: True effects were not recovered by local method!")
  > if (!all(res.opt)) stop("Test ida: True effects were not recovered by optimal method!")
  > 
  > ## *one* test for  method="global" :
  > eff.g.est <- Rnd(ida(x,y, mcov, myCPDAG, method="global", verbose=TRUE))
  Fit - y: 4 x: 5 |b.hat= 0 
  > stopifnot(eff.est.local == eff.g.est)
  > 
  > ## cat('Time elapsed additionally: ', proc.time() - .pt,"\n")
  > if (doExtras) {
  + ## another special case (from Raphael Gervais)
  + set.seed(123)
  + p <- 7
  + myDAG <- randomDAG(p, prob = 0.2) ## true DAG
  + amatT <- as(myDAG, "matrix") # weighted adjacency matrix of true DAG
  + effT <- Rnd(amatT[2,3]*amatT[3,5]) # Causal effect of 2 on 5 from true DAG weighted matrix
  + myCPDAG <- dag2cpdag(myDAG) ## true CPDAG
  + covTrue <- trueCov(myDAG) ## true covariance matrix
  + effG <- Rnd(ida(2,5, covTrue,myCPDAG,method = "global"))
  + 
  + if (!(effT %in% effG)) stop("Test ida special case: True effects were not recovered!")
  + 
  + ##################################################
  + ## Tests for method = optimal; sets x
  + ##################################################
  + set.seed(1)
  + V <- sample(c("C1","C2","C3","X","C4","C5","C6","C7","Y","C8","C9","C10"))
  + myDAG <- randomDAG(20, 0.3)
  + mcov <- trueCov(myDAG)
  + amat <- t(as(myDAG,"matrix"))
  + graphEst <- dag2cpdag(myDAG)
  + amat.cpdag <- t(as(graphEst,"matrix"))
  + 
  + opt <- ida(c(4,6),12,trueCov(myDAG),graphEst,method="optimal",type="cpdag") 
  + RRC <- jointIda(c(4,6),12,trueCov(myDAG),graphEst,technique="RRC")
  + stopifnot(all.equal(opt,RRC, tolerance = 0.01))
  + 
  + opt <- ida(c(4,6),c(10,19),trueCov(myDAG),graphEst,method="optimal",type="cpdag")
  + RRC <- jointIda(c(4,6),c(10,19),trueCov(myDAG),graphEst,technique="MCD")
  + stopifnot(all.equal(opt,RRC, tolerance = 0.01))
  + 
  + opt <- ida(c(5,10),c(3),trueCov(myDAG),graphEst,method="optimal",type="cpdag")
  + RRC <- jointIda(c(5,10),c(3),trueCov(myDAG),graphEst,technique="RRC")
  + stopifnot(all.equal(opt,RRC, tolerance = 0.01))
  + 
  + ## sometimes they differ
  + ## ida(c(5,20),c(10),trueCov(myDAG),graphEst,method="optimal",type="cpdag")
  + ## jointIda(c(5,20),c(10),trueCov(myDAG),graphEst,technique="RRC")
  + 
  + ##################################################
  + ## Tests related to examples (use dontrun for slow global option there)
  + ##################################################
  + set.seed(123)
  + p <- 10
  + myDAG <- randomDAG(p, prob = 0.2) ## true DAG
  + myCPDAG <- dag2cpdag(myDAG) ## true CPDAG
  + myPDAG <- addBgKnowledge(myCPDAG,2,3) ## true PDAG with background knowledge 2 -> 3
  + covTrue <- trueCov(myDAG) ## true covariance matrix
  + 
  + ## simulate Gaussian data from the true DAG
  + n <- 10000
  + dat <- rmvDAG(n, myDAG)
  + 
  + ## estimate CPDAG and PDAG -- see  help(pc)
  + suffStat <- list(C = cor(dat), n = n)
  + pc.fit <- pc(suffStat, indepTest = gaussCItest, p=p, alpha = 0.01)
  + pc.fit.pdag <- addBgKnowledge(pc.fit@graph,2,3)
  + 
  + ## Supppose that we know the true CPDAG and covariance matrix
  + (l.ida.cpdag <- ida(3,10, covTrue, myCPDAG, method = "local", type = "cpdag"))
  + (o.ida.cpdag <- ida(3,10, covTrue, myCPDAG, method = "optimal", type = "cpdag"))
  + (g.ida.cpdag <- ida(3,10, covTrue, myCPDAG, method = "global", type = "cpdag"))
  + ## All three methods produce the same unique values.
  + stopifnot(all.equal(sort(unique(Rnd7(g.ida.cpdag))),
  +                     sort(unique(Rnd7(l.ida.cpdag)))))
  + stopifnot(all.equal(sort(unique(Rnd7(g.ida.cpdag))),
  +                     sort(unique(Rnd7(as.vector(o.ida.cpdag))))))
  + 
  + ## Supppose that we know the true PDAG and covariance matrix
  + (l.ida.pdag <- ida(3,10, covTrue, myPDAG, method = "local", type = "pdag"))
  + (o.ida.pdag <- ida(3,10, covTrue, myPDAG, method = "optimal", type = "pdag"))
  + (g.ida.pdag <- ida(3,10, covTrue, myPDAG, method = "global", type = "pdag"))
  + ## All three methods produce the same unique values.
  + stopifnot(all.equal(sort(unique(Rnd7(g.ida.pdag))),
  +                     sort(unique(Rnd7(l.ida.pdag)))))
  + stopifnot(all.equal(sort(unique(Rnd7(g.ida.pdag))),
  +                     sort(unique(Rnd7(as.vector(o.ida.pdag))))))
  + 
  + ## From the true DAG, we can compute the true causal effect of 3 on 10
  + (ce.3.10 <- causalEffect(myDAG, 10, 3))
  + ## Indeed, this value is contained in the values found by all methods
  + 
  + ## When working with data we have to use the estimated CPDAG and
  + ## the sample covariance matrix
  + (l.ida.est.cpdag <- ida(3,10, cov(dat), pc.fit@graph, method = "local", type = "cpdag"))
  + (o.ida.est.cpdag <- ida(3,10, cov(dat), pc.fit@graph, method = "optimal", type = "cpdag"))
  + (g.ida.est.cpdag <- ida(3,10, cov(dat), pc.fit@graph, method = "global", type = "cpdag"))
  + ## The unique values of the local and the global method are still identical.
  + stopifnot(all.equal(sort(unique(Rnd7(g.ida.est.cpdag))), sort(unique(Rnd7(l.ida.est.cpdag)))))
  + ## While not identical, the values of the optimal method are very similar.
  + stopifnot(all.equal(sort(o.ida.est.cpdag), sort(l.ida.est.cpdag), tolerance = 0.025))
  + ## The true causal effect is contained in all three sets, up to a small
  + ## estimation error (0.118 vs. 0.112 with true value 0.114) 
  + stopifnot(all.equal(ce.3.10, min(l.ida.est.cpdag), tolerance = 0.04))
  + stopifnot(all.equal(ce.3.10, min(o.ida.est.cpdag), tolerance = 0.02))
  + 
  + ## Similarly, when working with data and background knowledge we have to use the estimated PDAG and
  + ## the sample covariance matrix
  + (l.ida.est.pdag <- ida(3,10, cov(dat), pc.fit.pdag, method = "local", type = "pdag"))
  + (o.ida.est.pdag <- ida(3,10, cov(dat), pc.fit.pdag, method = "optimal", type = "pdag"))
  + (g.ida.est.pdag <- ida(3,10, cov(dat), pc.fit.pdag, method = "global", type = "pdag"))
  + ## The unique values of the local and the global method are still identical.
  + stopifnot(all.equal(sort(unique(Rnd7(g.ida.est.pdag))), sort(unique(Rnd7(l.ida.est.pdag)))))
  + ## While not necessarily identical, the values of the optimal method will be similar.
  + stopifnot(all.equal(sort(Rnd7(o.ida.est.pdag)), sort(Rnd7(l.ida.est.pdag)), tolerance = 0.08))
  + ## The true causal effect is contained in both sets, up to a small estimation error
  + stopifnot(all.equal(ce.3.10, min(l.ida.est.pdag), tolerance = 0.04))
  + stopifnot(all.equal(ce.3.10, min(o.ida.est.pdag), tolerance = 0.02))
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
     6.06    0.42    6.42 
Running the tests in 'tests/test_isValidGraph.R' failed.
Complete output:
  > library(pcalg)
  > possible_p <- c(seq(10,30,by=10))
  > possible_neighb_size <- c(3:10)
  > 
  > num_setings <-10
  > num_rep <- 2
  > 
  > ##amount of bg knowledge in percent
  > p_bg <- seq(0,1,by=0.1)
  > 
  > counterr <- 0
  > 
  > seed1 <- 1 
  > for (r in seed1:num_setings) {
  +   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)
  + 
  +   ## 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
  + 
  +     amat.cpdag <- cpdag.amat
  +     if (!isValidGraph(amat.cpdag,type="cpdag")){
  +       counterr <- counterr +1
  +     }
  +     if (!isValidGraph(dag.amat,type="dag")){
  +       counterr <- counterr +1
  +     }
  +   }
  + }
  > counterr
  [1] 0
  > 
  > amat <- matrix(c(0,1,0, 0,0,1, 0,0,0), 3,3)
  > colnames(amat) <- rownames(amat) <- letters[1:3]
  > ## graph::plot(as(t(amat), "graphNEL"))             
  > stopifnot(isValidGraph(amat = amat, type = "dag"), ## is a valid DAG
  +   !isValidGraph(amat = amat, type = "cpdag"), ## is a valid CPDAG 
  +   isValidGraph(amat = amat, type = "pdag") ) ## is a valid PDAG
  > 
  > ## a -- b -- c
  > amat <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3,3)
  > colnames(amat) <- rownames(amat) <- letters[1:3]
  > ## plot(as(t(amat), "graphNEL"))             
  > stopifnot(!isValidGraph(amat = amat, type = "dag"), ## not a valid DAG
  +   isValidGraph(amat = amat, type = "cpdag"), ## is a valid CPDAG
  +   isValidGraph(amat = amat, type = "pdag") )## is a valid PDAG
  > 
  > ## a -- b -- c -- d -- a
  > amat <- matrix(c(0,1,0,1, 1,0,1,0, 0,1,0,1, 1,0,1,0), 4,4)
  > colnames(amat) <- rownames(amat) <- letters[1:4]
  > ## plot(as(t(amat), "graphNEL"))             
  > stopifnot(!isValidGraph(amat = amat, type = "dag"),  ## not a valid DAG
  +   !isValidGraph(amat = amat, type = "cpdag"), ## not a valid CPDAG
  +   !isValidGraph(amat = amat, type = "pdag") ) ## not a valid PDAG
  > 
  > proc.time()
     user  system elapsed 
     2.42    0.39    2.79 
Running the tests in 'tests/test_optAdjSet.R' failed.
Complete output:
  > library(pcalg)
  > suppressWarnings(RNGversion("3.5.0"))
  > 
  > set.seed(123)
  > myDAG <- randomDAG(20, 0.3)
  > mcov <- trueCov(myDAG)
  > amat <- t(as(myDAG,"matrix"))
  > amat[which(amat!=0)] <- 1
  > graphEst <- dag2cpdag(myDAG)
  > ## plot(graphEst)
  > amat.cpdag <- t(as(graphEst,"matrix"))
  > 
  > stopifnot(sort(optAdjSet(amat,2,6))==c())
  > 
  > 
  > stopifnot(pcalg:::isAmenable(amat.cpdag,5,16,type="cpdag"))
  > stopifnot(sort(optAdjSet(amat.cpdag,5,16))==c(1,2,3,4,6,7,8,11,13))
  > stopifnot(optAdjSet(amat.cpdag,5,16)==optAdjSet(amat,5,16))
  > 
  > stopifnot(pcalg:::isAmenable(amat.cpdag,2,12,type="cpdag"))
  > stopifnot(sort(optAdjSet(amat.cpdag,2,12))==c(5,7))
  > stopifnot(optAdjSet(amat.cpdag,2,12)==optAdjSet(amat,2,12))
  > 
  > stopifnot(optAdjSet(amat.cpdag,2,c(3,7,12))==optAdjSet(amat,2,12))
  
  
  Not all nodes in Y are descendants of the nodes in X 
  The algorithm will continue with the pruned Y= 12 .
  >           
  > 
  > proc.time()
     user  system elapsed 
     1.70    0.34    2.01 
Running the tests in 'tests/test_opttarget.R' failed.
Complete output:
  > # Tests for optimal intervention targets.
  > 
  > library(pcalg)
  > 
  > # Test graph: essential graph of Figure 4 in [1]. It has an unoriented component
  > # with 3 vertices (1, 2, 3) and one with 2 vertices (4, 5).
  > #
  > # [1] A. Hauser, P. Bühlmann: Two optimal strategies for active learning of
  > # causal models from interventions. Proceedings of the 6th European Workshop on
  > # Probabilistic Graphical Models (PGM-2012), pp. 123 to 130, 2012.
  > essgraph <- new("EssGraph",
  +                 nodes = as.character(1:5),
  +                 in.edges = list(2, c(1, 3), 2, c(2, 3, 5), c(2, 3, 4)),
  +                 targets = list(integer(0), 1:3))
  > 
  > # Optimal single vertex intervention: vertex 2 (==> all edges between nodes 1,
  > # 2, 3 become orientable).
  > stopifnot(opt.target(essgraph, max.size = 1, use.node.names = FALSE) == 2)
  > stopifnot(opt.target(essgraph, max.size = 1) == "2")
  > 
  > # Optimal intervention of arbitrary size: vertices 1, 3, 4 (makes all edges
  > # orientable).
  > stopifnot(all.equal(opt.target(essgraph, max.size = 5, use.node.names = FALSE),
  +                     c(1, 3, 4)))
  > stopifnot(all.equal(opt.target(essgraph), c("1", "3", "4")))
  > 
  > proc.time()
     user  system elapsed 
     1.10    0.21    1.32 
Running the tests in 'tests/test_pdag2allDags.R' failed.
Complete output:
  > library(pcalg)
  > 
  > val <- rep(FALSE, 9)
  > ## Test 1
  > gm <- rbind(c(0,1),
  +             c(1,0))
  > colnames(gm) <- rownames(gm) <- LETTERS[1:2]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > ## plotAllDags(res)
  > ## for (i in 1:2) {cat(paste("c(",paste(res$dags[i,], collapse = ","),"),",sep=""),"\n")}
  > res.truth <- rbind(c(0,1,0,0), 
  +                    c(0,0,1,0))
  > val[1] <- identical(res.truth, res$dags)
  > 
  > ## Test 2
  > gm <- rbind(c(0,0),
  +             c(0,0))
  > colnames(gm) <- rownames(gm) <- LETTERS[1:2]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > res.truth <- rbind(c(0,0,0,0))
  > val[2] <- identical(res.truth, res$dags)
  > 
  > ## Test 3, non-collider
  > gm <- rbind(c(0,1,0),
  +             c(1,0,1),
  +             c(0,1,0))
  > colnames(gm) <- rownames(gm) <- LETTERS[1:ncol(gm)]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > res.truth <- rbind(c(0,1,0,0,0,1,0,0,0), 
  +                    c(0,1,0,0,0,0,0,1,0), 
  +                    c(0,0,0,1,0,0,0,1,0))
  > val[3] <- identical(res.truth, res$dags)
  > ## for (i in 1:3) {cat(paste("c(",paste(res$dags[i,], collapse = ","),"),",sep=""),"\n")}
  > ## plotAllDags(res)
  > 
  > ## Test 4, collider
  > gm <- rbind(c(0,0,0),
  +             c(1,0,1),
  +             c(0,0,0))
  > colnames(gm) <- rownames(gm) <- LETTERS[1:ncol(gm)]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > res.truth <- rbind(c(0,0,0,1,0,1,0,0,0))
  > val[4] <- identical(res.truth, res$dags)
  > ## plotAllDags(res)
  > 
  > ## Test 5, trick question
  > gm <- rbind(c(0,0,0),
  +             c(1,0,1),
  +             c(0,1,0))
  > colnames(gm) <- rownames(gm) <- LETTERS[1:ncol(gm)]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > res.truth <- rbind(c(0,0,0,1,0,0,0,1,0))
  > val[5] <- identical(res.truth, res$dags)
  > ## plotAllDags(res)
  > 
  > ## Test 6,complete
  > gm <- rbind(c(0,1,1),
  +             c(1,0,1),
  +             c(1,1,0))
  > colnames(gm) <- rownames(gm) <- LETTERS[1:ncol(gm)]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > ## for (i in 1:6) {cat(paste("c(",paste(res$dags[i,], collapse = ","),"),",sep=""),"\n")}
  > res.truth <- rbind(c(0,1,1,0,0,1,0,0,0), 
  +                    c(0,1,1,0,0,0,0,1,0), 
  +                    c(0,0,1,1,0,1,0,0,0), 
  +                    c(0,0,0,1,0,1,1,0,0), 
  +                    c(0,1,0,0,0,0,1,1,0), 
  +                    c(0,0,0,1,0,0,1,1,0))
  > val[6] <- identical(res.truth, res$dags)
  > ## plotAllDags(res)
  > 
  > ## Test 7, 4 nodes: being really mean
  > ## No DAG possible
  > gm <- rbind(c(0,1,1,0),
  +             c(1,0,0,1),
  +             c(1,0,0,1),
  +             c(0,1,1,0))
  > colnames(gm) <- rownames(gm) <- LETTERS[1:ncol(gm)]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > val[7] <- is.null(res$dags)
  > ## plotAllDags(res)
  > 
  > ## Test 8, 4 nodes
  > gm <- rbind(c(0,1,1,0),
  +             c(1,0,0,0),
  +             c(1,0,0,0),
  +             c(0,1,1,0))
  > colnames(gm) <- rownames(gm) <- LETTERS[1:ncol(gm)]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > res.truth <- rbind(c(0,0,1,0,1,0,0,0,0,0,0,0,0,1,1,0),
  +                    c(0,0,0,0,1,0,0,0,1,0,0,0,0,1,1,0),
  +                    c(0,1,0,0,0,0,0,0,1,0,0,0,0,1,1,0))
  > val[8] <- identical(res.truth, res$dags)
  > ## plotAllDags(res)
  > 
  > ## Test 9, 5 nodes; D -> E must always be there
  > gm <- rbind(c(0,1,1,0,0),
  +             c(1,0,0,0,0),
  +             c(1,0,0,0,0),
  +             c(0,1,1,0,1),
  +             c(0,0,0,1,0))
  > 
  > colnames(gm) <- rownames(gm) <- LETTERS[1:ncol(gm)]
  > res <- pdag2allDags(gm, verbose = FALSE)
  > 
  > res.truth <- rbind(c(0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,0),
  +                    c(0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,1,0),
  +                    c(0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,1,0))
  > val[9] <- identical(res.truth, res$dags)
  > 
  > if (!all(val)) stop("Error in testing pdag2allDags!\n")
  > 
  > proc.time()
     user  system elapsed 
     1.23    0.18    1.40 
Running the tests in 'tests/test_possDeAn.R' failed.
Complete output:
  > library(pcalg)
  > 
  > ## possDe
  > ## a -> b -- c
  > amat <- matrix(c(0,1,0, 0,0,1, 0,1,0), 3,3)
  > colnames(amat) <- rownames(amat) <- letters[1:3]
  > ## plot(as(t(amat), "graphNEL"))
  > 
  > stopifnot(
  +   all(possDe(m = amat, x = 1, possible = TRUE, ds = FALSE, type = "pdag") == c(1,2,3)),
  +   all(possDe(m = amat, x = 1, possible = FALSE, ds = FALSE, type = "pdag") == c(1,2)),
  +   all(possDe(m = amat, x = 1, y = 2, possible = TRUE, ds = FALSE, type = "pdag") == 1) )
  > 
  > ## possAn
  > ## a -- b -> c
  > amat <- matrix(c(0,1,0, 1,0,1, 0,0,0), 3,3)
  > colnames(amat) <- rownames(amat) <- letters[1:3]
  > plot(as(t(amat), "graphNEL"))
  
  Attaching package: 'BiocGenerics'
  
  The following objects are masked from 'package:stats':
  
      IQR, mad, sd, var, xtabs
  
  The following objects are masked from 'package:base':
  
      Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
      as.data.frame, basename, cbind, colnames, dirname, do.call,
      duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
      lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
      pmin.int, rank, rbind, rownames, sapply, setdiff, table, tapply,
      union, unique, unsplit, which.max, which.min
  
  > 
  > stopifnot(
  +   all(possAn(m = amat, x = 3, possible = TRUE, ds = FALSE, type = "pdag") == c(1,2,3)),
  +   all(possAn(m = amat, x = 3, y = 2, possible = TRUE, ds = FALSE, type = "pdag") == 3) )
  > 
  > proc.time()
     user  system elapsed 
     1.85    0.32    2.12 
Running the tests in 'tests/test_randomDAG.R' failed.
Complete output:
  > library(pcalg)
  > 
  > set.seed(42)
  > acyc <- edgeDistr <- wgtDistr <- 0
  > 
  > p <- 10
  > s <- 0.2
  > lB <- 0
  > uB <- 2
  > wgts <- NULL
  > nEdg <- 0
  > nreps <- 10
  > acyc <- rep(FALSE,nreps)
  > 
  > for (i in 1:nreps) {
  +   myDAG <- randomDAG(p,s,lB,uB)
  +   amat <- as(myDAG,"matrix")
  +   acyc[i] <- ggm::isAcyclic(amat)
  +   wgts <- c(wgts,amat[amat!=0])
  +   nEdg <- nEdg + graph::numEdges(myDAG)
  + }
  > 
  > ## Check number of edges
  > nPossibleEdg <- nreps*p*(p-1)/2
  > edge.test <- binom.test(nEdg,nPossibleEdg,p=s,alternative="two.sided")
  > if (edge.test$p.value < 0.05) {
  +   stop("Test of randomDAG: Number of edges is not as theory predicts!")
  + }
  > 
  > ## Check distribution of wgts
  > if (min(wgts) < lB || max(wgts) > uB) {
  +   stop("Test of randomDAG: Weights are not within specified ranges!")
  + }
  > 
  > ## Check whether graphs are acyclic
  > if (!all(acyc))
  +   stop("Test of randomDAG: Graph is not acyclic!")
  > 
  > ## Check special case
  > if (graph::numEdges(randomDAG(n=2,prob=1)) != 1) {
  +     stop("Test of randomDAG: Error with special case n=2, prob=1!")
  + }
  > 
  > cat('Time elapsed: ', proc.time(),"\n")
  Time elapsed:  1.34 0.26 1.59 NA NA 
  > 
  > proc.time()
     user  system elapsed 
     1.34    0.26    1.59 
Running the tests in 'tests/test_rfci.R' failed.
Complete output:
  > library(pcalg)
  > doExtras <- pcalg:::doExtras()
  > 
  > source(system.file(package="Matrix", "test-tools-1.R", mustWork=TRUE))
  Loading required package: tools
  > ##--> showProc.time(), assertError(), relErrV(), ...
  > R.home(); sessionInfo() # helping package maintainers to debug ...
  [1] "D:/RCompile/recent/R"
  R Under development (unstable) (2024-03-22 r86169 ucrt)
  Platform: x86_64-w64-mingw32/x64
  Running under: Windows Server 2022 x64 (build 20348)
  
  Matrix products: default
  
  
  locale:
  [1] LC_COLLATE=C                 LC_CTYPE=German_Germany.utf8
  [3] LC_MONETARY=C                LC_NUMERIC=C                
  [5] LC_TIME=C                   
  
  time zone: Europe/Berlin
  tzcode source: internal
  
  attached base packages:
  [1] tools     stats     graphics  grDevices utils     datasets  methods  
  [8] base     
  
  other attached packages:
  [1] pcalg_2.7-11
  
  loaded via a namespace (and not attached):
   [1] graph_1.80.0        igraph_2.0.3        DEoptimR_1.1-3     
   [4] clue_0.3-65         magrittr_2.0.3      fastICA_1.2-4      
   [7] abind_1.4-5         BiocGenerics_0.48.1 pkgconfig_2.0.3    
  [10] stats4_4.4.0        lifecycle_1.0.4     cli_3.6.2          
  [13] sfsmisc_1.1-17      ggm_2.5.1           corpcor_1.6.10     
  [16] robustbase_0.99-2   compiler_4.4.0      RBGL_1.78.0        
  [19] cluster_2.1.6       bdsmatrix_1.3-7     Rcpp_1.0.12        
  [22] BiocManager_1.30.22 rlang_1.1.3        
  > .libPaths()
  [1] "D:/temp/RtmpmCFgj4/RLIBS_1177073ea38fa"
  [2] "D:/RCompile/recent/R/library"          
  > packageDescription("pcalg")
  Package: pcalg
  Version: 2.7-11
  Date: 2024-02-12
  Title: Methods for Graphical Models and Causal Inference
  Description: Functions for causal structure learning and causal
          inference using graphical models. The main algorithms for
          causal structure learning are PC (for observational data
          without hidden variables), FCI and RFCI (for observational data
          with hidden variables), and GIES (for a mix of data from
          observational studies (i.e. observational data) and data from
          experiments involving interventions (i.e. interventional data)
          without hidden variables). For causal inference the IDA
          algorithm, the Generalized Backdoor Criterion (GBC), the
          Generalized Adjustment Criterion (GAC) and some related
          functions are implemented. Functions for incorporating
          background knowledge are provided.
  Maintainer: Markus Kalisch <kalisch@stat.math.ethz.ch>
  Authors@R: c(person("Markus","Kalisch",
          email="kalisch@stat.math.ethz.ch", role=c("aut","cre")),
          person("Alain", "Hauser", role="aut"),
          person("Martin","Maechler", role="aut"), person("Diego",
          "Colombo", role="ctb"), person("Doris", "Entner", role="ctb"),
          person("Patrik", "Hoyer", role="ctb"), person("Antti",
          "Hyttinen", role="ctb"), person("Jonas", "Peters", role="ctb"),
          person("Nicoletta", "Andri", role="ctb"), person("Emilija",
          "Perkovic", role="ctb"), person("Preetam", "Nandy",
          role="ctb"), person("Philipp", "Ruetimann",
          role="ctb"),person("Daniel", "Stekhoven", role="ctb"),
          person("Manuel", "Schuerch", role="ctb"), person("Marco",
          "Eigenmann", role="ctb"), person("Leonard", "Henckel", role =
          "ctb"),person("Joris", "Mooij", role = "ctb"))
  Depends: R (>= 3.5.0)
  LinkingTo: Rcpp (>= 0.11.0), RcppArmadillo, BH
  Imports: stats, graphics, utils, methods, abind, graph, RBGL, igraph,
          ggm, corpcor, robustbase, vcd, Rcpp, bdsmatrix, sfsmisc (>=
          1.0-26), fastICA, clue
  Suggests: MASS, Matrix, Rgraphviz, mvtnorm, huge, ggplot2, dagitty
  NeedsCompilation: yes
  Encoding: UTF-8
  License: GPL (>= 2)
  URL: https://pcalg.r-forge.r-project.org/
  Packaged: 2024-02-12 10:30:09 UTC; kalischm
  Author: Markus Kalisch [aut, cre], Alain Hauser [aut], Martin Maechler
          [aut], Diego Colombo [ctb], Doris Entner [ctb], Patrik Hoyer
          [ctb], Antti Hyttinen [ctb], Jonas Peters [ctb], Nicoletta
          Andri [ctb], Emilija Perkovic [ctb], Preetam Nandy [ctb],
          Philipp Ruetimann [ctb], Daniel Stekhoven [ctb], Manuel
          Schuerch [ctb], Marco Eigenmann [ctb], Leonard Henckel [ctb],
          Joris Mooij [ctb]
  Repository: CRAN
  Date/Publication: 2024-02-12 12:00:05 UTC
  Built: R 4.4.0; x86_64-w64-mingw32; 2024-03-23 03:42:08 UTC; windows
  Archs: x64
  
  -- File: D:/RCompile/CRANpkg/lib/4.4/pcalg/Meta/package.rds 
  > packageDescription("Matrix")
  Package: Matrix
  Version: 1.7-0
  VersionNote: do also bump src/version.h, inst/include/Matrix/version.h
  Date: 2024-03-16
  Priority: recommended
  Title: Sparse and Dense Matrix Classes and Methods
  Description: A rich hierarchy of sparse and dense matrix classes,
          including general, symmetric, triangular, and diagonal matrices
          with numeric, logical, or pattern entries.  Efficient methods
          for operating on such matrices, often wrapping the 'BLAS',
          'LAPACK', and 'SuiteSparse' libraries.
  License: GPL (>= 2) | file LICENCE
  URL: https://Matrix.R-forge.R-project.org
  BugReports: https://R-forge.R-project.org/tracker/?atid=294&group_id=61
  Contact: Matrix-authors@R-project.org
  Authors@R: c(person("Douglas", "Bates", role = "aut", comment = c(ORCID
          = "0000-0001-8316-9503")), person("Martin", "Maechler", role =
          c("aut", "cre"), email = "mmaechler+Matrix@gmail.com", comment
          = c(ORCID = "0000-0002-8685-9910")), person("Mikael", "Jagan",
          role = "aut", comment = c(ORCID = "0000-0002-3542-2938")),
          person("Timothy A.", "Davis", role = "ctb", comment = c(ORCID =
          "0000-0001-7614-6899", "SuiteSparse libraries", "collaborators
          listed in dir(system.file(\"doc\", \"SuiteSparse\",
          package=\"Matrix\"), pattern=\"License\", full.names=TRUE,
          recursive=TRUE)")), person("George", "Karypis", role = "ctb",
          comment = c(ORCID = "0000-0003-2753-1437", "METIS library",
          "Copyright: Regents of the University of Minnesota")),
          person("Jason", "Riedy", role = "ctb", comment = c(ORCID =
          "0000-0002-4345-4200", "GNU Octave's condest() and
          onenormest()", "Copyright: Regents of the University of
          California")), person("Jens", "Oehlschlägel", role = "ctb",
          comment = "initial nearPD()"), person("R Core Team", role =
          "ctb", comment = "base R's matrix implementation"))
  Depends: R (>= 4.4.0), methods
  Imports: grDevices, graphics, grid, lattice, stats, utils
  Suggests: MASS, datasets, sfsmisc, tools
  Enhances: SparseM, graph
  LazyData: no
  LazyDataNote: not possible, since we use data/*.R and our S4 classes
  BuildResaveData: no
  Encoding: UTF-8
  NeedsCompilation: yes
  Packaged: 2024-03-19 17:15:14 UTC; maechler
  Author: Douglas Bates [aut] (<https://orcid.org/0000-0001-8316-9503>),
          Martin Maechler [aut, cre]
          (<https://orcid.org/0000-0002-8685-9910>), Mikael Jagan [aut]
          (<https://orcid.org/0000-0002-3542-2938>), Timothy A. Davis
          [ctb] (<https://orcid.org/0000-0001-7614-6899>, SuiteSparse
          libraries, collaborators listed in dir(system.file("doc",
          "SuiteSparse", package="Matrix"), pattern="License",
          full.names=TRUE, recursive=TRUE)), George Karypis [ctb]
          (<https://orcid.org/0000-0003-2753-1437>, METIS library,
          Copyright: Regents of the University of Minnesota), Jason Riedy
          [ctb] (<https://orcid.org/0000-0002-4345-4200>, GNU Octave's
          condest() and onenormest(), Copyright: Regents of the
          University of California), Jens Oehlschlägel [ctb] (initial
          nearPD()), R Core Team [ctb] (base R's matrix implementation)
  Maintainer: Martin Maechler <mmaechler+Matrix@gmail.com>
  Repository: CRAN
  Date/Publication: 2024-03-22 09:08:25 UTC
  Built: R 4.4.0; x86_64-w64-mingw32; 2024-03-23 02:52:59 UTC; windows
  Archs: x64
  
  -- File: D:/temp/RtmpmCFgj4/RLIBS_1177073ea38fa/Matrix/Meta/package.rds 
  > 
  > ## load the functions for the simulations of this paper
  > ## source("/u/colombo/Diss/RAusw/First_paper_RFCI/functions_for_the_simulations.R")
  > 
  > ## RFCI improves the output
  > ##______________________________________________
  > ## Input:  L1=1; L2=2; X1=6; X2=4; X3=3; X4=5; X5=7; X6=8
  > ## Output: X1=4; X2=2; X3=1; X4=3; X5=5; X6=6
  > 
  > amat <- rbind(0,# 2 3 4 5 6 7 8
  +               c(0,0,0,0,1,1,0,0),
  +               c(0,1,0,1,0,1,0,0),
  +               c(1,0,0,0,0,1,0,0),
  +               c(0,0,0,0,0,1,0,0),
  +               c(1,1,0,0,0,0,0,0),
  +               c(0,0,0,1,1,0,0,0),
  +               0)
  > colnames(amat) <- rownames(amat) <- as.character(1:8)
  > Matrix::Matrix(amat) # to "visualize"
  8 x 8 sparse Matrix of class "dgCMatrix"
    1 2 3 4 5 6 7 8
  1 . . . . . . . .
  2 . . . . 1 1 . .
  3 . 1 . 1 . 1 . .
  4 1 . . . . 1 . .
  5 . . . . . 1 . .
  6 1 1 . . . . . .
  7 . . . 1 1 . . .
  8 . . . . . . . .
  > L <- c(1,2)
  > V <- as.character(1:8)
  > edL <- setNames(vector("list", length=length(V)), V)
  > edL[[6]] <- list(edges=NULL, weights=NULL)
  > edL[[8]] <- list(edges=NULL, weights=NULL)
  > edL[[4]] <- list(edges=c(7,8),   weights=c(abs(rnorm(1)),abs(rnorm(1))))
  > edL[[3]] <- list(edges=c(4,5,8), weights=c(abs(rnorm(1)),abs(rnorm(1)),abs(rnorm(1))))
  > edL[[5]] <- list(edges=c(6,8),   weights=c(abs(rnorm(1)),abs(rnorm(1))))
  > edL[[7]] <- list(edges= 8,       weights=abs(rnorm(1)))
  > edL[[1]] <- list(edges=c(4,6),   weights=c(abs(rnorm(1)),abs(rnorm(1))))
  > edL[[2]] <- list(edges=c(5,7),   weights=c(abs(rnorm(1)),abs(rnorm(1))))
  > g <- new("graphNEL", nodes=V, edgeL=edL, edgemode="directed")
  > if(dev.interactive())
  +     plot(g)
  > 
  > ## Compute the true covariance matrix of g
  > cov.mat <- trueCov(g)
  > ## Delete rows and columns which belong to L
  > true.cov <- cov.mat[-L,-L]
  > ## Transform it in a correlation matrix
  > true.corr <- cov2cor(true.cov)
  > 
  > suffStat <- list(C=true.corr, n=10^9)
  > showSys.time(pop.fci1 <-
  +              fci(suffStat, gaussCItest, labels=V[-L],
  +                  alpha=0.9999, doPdsep=TRUE,verbose=FALSE)@amat)
  Time    user  system elapsed 
  Time    0.11    0.01    0.12 
  > 
  > showSys.time(pop.rfci1 <-
  +              rfci(suffStat, gaussCItest, labels=V[-L],
  +                   alpha=0.9999, verbose=FALSE)@amat)
  Time    user  system elapsed 
  Time    0.05    0.00    0.05 
  > 
  > if (any(pop.fci1 != pop.rfci1)) {
  +   stop("Test of RFCI wrong: small example!")
  + }
  > 
  > ## if (doExtras) {
  > if (FALSE) {
  + ## Thomas' example (version number 8) about discriminating path orientation rule
  + 
  + V <- as.character(1:25)
  + edL <- setNames(vector("list", length=length(V)), V)
  + edL[[ 1]] <- list(edges=c(14,18),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[ 2]] <- list(edges=c(16,18),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[ 3]] <- list(edges=c(16,24),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[ 4]] <- list(edges=c(18,24),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[ 5]] <- list(edges=c(15,25),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[ 6]] <- list(edges=c(17,19),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[ 7]] <- list(edges=c(14,19),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[ 8]] <- list(edges=c(14,24),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[ 9]] <- list(edges=c(19,20),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[10]] <- list(edges=c(20,25),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[11]] <- list(edges=c(23,25),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[12]] <- list(edges=c(22,24),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[13]] <- list(edges=c(21,23),weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[14]] <- list(edges=NULL,    weights=NULL)
  + edL[[15]] <- list(edges=c(16,17,24),weights=c(abs(rnorm(1)),abs(rnorm(1)),abs(rnorm(1))))
  + edL[[16]] <- list(edges=c(19,25),   weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[17]] <- list(edges=c(18,24,25),weights=c(abs(rnorm(1)),abs(rnorm(1)),abs(rnorm(1))))
  + edL[[18]] <- list(edges=c(21,25),   weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[19]] <- list(edges=c(23,24,25),weights=c(abs(rnorm(1)),abs(rnorm(1)),abs(rnorm(1))))
  + edL[[20]] <- list(edges= 24,        weights=abs(rnorm(1)))
  + edL[[21]] <- list(edges=c(22,25),   weights=c(abs(rnorm(1)),abs(rnorm(1))))
  + edL[[22]] <- list(edges= 25,        weights=abs(rnorm(1)))
  + edL[[23]] <- list(edges= 24,        weights=abs(rnorm(1)))
  + edL[[24]] <- list(edges=NULL,weights=NULL)
  + edL[[25]] <- list(edges=NULL,weights=NULL)
  + (g <- new("graphNEL", nodes=V, edgeL=edL,edgemode="directed"))
  + 
  + if(dev.interactive())
  +     plot(g)
  + 
  + ## Latent variables (all having no parents):
  + L <- c(1:13)
  + 
  + ## Compute the true covariance matrix of g
  + cov.mat <- trueCov(g)
  + ## Delete rows and columns which belong to L
  + true.cov <- cov.mat[-L,-L]
  + ## Transform it in a correlation matrix
  + true.corr <- cov2cor(true.cov)
  + suffStat <- list(C=true.corr, n=10^9)
  + p.tr <- dim(true.corr)[1]
  + showSys.time(pop.fci2  <-  fci(suffStat, gaussCItest, p=p.tr,
  +                                alpha=0.9999, doPdsep=TRUE)@amat)
  + showSys.time(pop.rfci2 <- rfci(suffStat, gaussCItest, p=p.tr,
  +                                alpha=0.9999)@amat)
  + 
  + if (any(pop.fci2 != pop.rfci2)) {
  +   stop("Test of RFCI wrong: big example!")
  + }
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
     2.90    0.35    3.25 
Running the tests in 'tests/test_searchAM.R' failed.
Complete output:
  > library(pcalg)
  > 
  > # Y-structure MAG
  > # Encode as adjacency matrix
  > p <- 10 # total number of variables
  > V <- c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10") # variable labels
  > # amat[i,j] = 0 iff no edge btw i,j
  > # amat[i,j] = 1 iff i *-o j
  > # amat[i,j] = 2 iff i *-> j
  > # amat[i,j] = 3 iff i *-- j
  > amat <- rbind(c(0,3,0,0,0,0,0,0,0,0),
  +               c(3,0,3,0,0,0,0,0,0,0),
  +               c(0,3,0,2,0,0,0,0,0,0),
  +               c(0,0,3,0,2,0,0,0,0,0),
  +               c(0,0,0,3,0,2,0,2,2,1),
  +               c(0,0,0,0,3,0,2,0,0,0),
  +               c(0,0,0,0,0,3,0,0,0,0),
  +               c(0,0,0,0,2,0,0,0,0,0),
  +               c(0,0,0,0,1,0,0,0,0,0),
  +               c(0,0,0,0,1,0,0,0,0,0))
  > rownames(amat)<-V
  > colnames(amat)<-V
  > 
  > stopifnot(all.equal(searchAM(amat,5,type = "an"), c(3,4,5))) # ancestors of X5
  > stopifnot(all.equal(searchAM(amat,5,type = "de"), c(5,6,7))) # descendants of X5
  > stopifnot(all.equal(searchAM(amat,5,type = "ant"), c(1,2,3,4,5))) # anteriors of X5
  > stopifnot(all.equal(searchAM(amat,5,type = "sp"), c(8))) # spouses of X5
  > stopifnot(all.equal(searchAM(amat,2,type = "nb"), c(1,3))) # neighbors of X2
  > stopifnot(all.equal(searchAM(amat,c(4,6),type = "pa"), c(3,5))) # parents of {X4,X6}
  > stopifnot(all.equal(searchAM(amat,c(3,5),type = "ch"), c(4,6))) # children of {X3,X5}
  > stopifnot(all.equal(searchAM(amat,5,type = "pde"), c(5,6,7,9,10))) # possible descendants of X5
  > 
  > proc.time()
     user  system elapsed 
     1.25    0.26    1.45 
Running the tests in 'tests/test_shd.R' failed.
Complete output:
  > library(pcalg)
  > suppressWarnings(RNGversion("3.5.0"))
  > set.seed(123)
  > 
  > ## expl 1
  > g1 <- randomDAG(5,0.5)
  > g2 <- randomDAG(5,0.5)
  > res1 <- (shd(g1,g2)==4)
  > 
  > ## expl 2
  > g3 <- dag2cpdag(g1)
  > res2 <- (shd(g3,g1)==3)
  > 
  > if(!all(c(res1,res2))) {
  +   stop("Test of shd: Theoretical value not matched!")
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.20    0.34    1.48 
Running the tests in 'tests/test_udag2pdag.R' failed.
Complete output:
  > library(pcalg)
  > suppressWarnings(RNGversion("3.5.0"))
  > .libPaths()
  [1] "D:/temp/RtmpmCFgj4/RLIBS_1177073ea38fa"
  [2] "D:/RCompile/recent/R/library"          
  > ## acyclic graphs
  > 
  > nreps <- 30
  > p <- 8
  > n <- 1000
  > 
  > for(u2pd in c("rand", "retry", "relaxed")) {
  +     cat("\n u2pd =", u2pd, "\n ------------\n")
  +     cyc.res <- logical(nreps)
  +     for (i in 1:nreps) {
  +         set.seed(i)
  +         myDAG <- randomDAG(p, prob = 0.2)
  +         d.mat <- rmvDAG(n, myDAG, errDist = "normal")
  +         res <- suppressWarnings(pcAlgo(d.mat, alpha = 0.05, directed=TRUE, u2pd = u2pd))
  +         ##                      ------ directed;  u2pd = "rand" --> udag2pdag()
  +         res.A <- wgtMatrix(res@graph)
  +         res.A[res.A!=0] <- 1
  +         undir.A <- res.A + t(res.A)
  +         undir.A[undir.A==1] <- 0
  +         undir.A[undir.A==2] <- 1
  +         res.dir <- res.A - undir.A
  +         cyc.res[i] <- ggm::isAcyclic(res.dir)
  +     }
  +     if (!all(cyc.res)) stop("Test of pcAlgo(*, directed): Cyclic part in PDAG!")
  + } ## for(u2pd ...)
  
   u2pd = rand 
   ------------
  
   u2pd = retry 
   ------------
  
   u2pd = relaxed 
   ------------
  > cat('Time elapsed: ', (.pt <- proc.time()),"\n")
  Time elapsed:  3.45 0.37 3.81 NA NA 
  > 
  > ## find collider correctly
  > set.seed(123)
  > myDAG <- randomDAG(3, prob = 0.5)
  > library(Matrix)
  > as(myDAG,"sparseMatrix")
  3 x 3 sparse Matrix of class "dgCMatrix"
    1 2         3
  1 . . 0.4680792
  2 . . 0.9464206
  3 . . .        
  > d.mat <- rmvDAG(n, myDAG, errDist = "normal")
  > res <- pcAlgo(d.mat, alpha = 0.05, corMethod = "standard",directed=TRUE)
  Warning message:
  In pcAlgo(d.mat, alpha = 0.05, corMethod = "standard", directed = TRUE) :
    pcAlgo() is deprecated and only kept for backward compatibility.
   Please use skeleton, pc, or fci instead
  
  > gEst <- wgtMatrix(res@graph)
  > gTrue <- rbind(0,0, c(1,1,0))
  > if(!all(gEst==gTrue)) stop("Test of udag2pdag: Problem finding a collider!")
  > 
  > cat('Time elapsed: ', proc.time() - .pt,'\n') # "stats"
  Time elapsed:  1.05 0.02 1.06 NA NA 
  > 
  > 
  > proc.time()
     user  system elapsed 
     4.50    0.39    4.87 
Running the tests in 'tests/test_wgtMatrix.R' failed.
Complete output:
  > library(pcalg)
  > suppressWarnings(RNGversion("3.5.0"))
  > set.seed(42)
  > 
  > ## wmat_ij is edge from j to i
  > g <- randomDAG(3,0.4)
  > wmat <- wgtMatrix(g)
  > if (!(wmat[2,1]!=0 & wmat[1,2]==0 & wmat[3,1]!=0 & wmat[3,2]==0)) {
  +   stop("Test of wgtMatrix: Something with orientation of edges is wrong!")
  + }
  > 
  > ## test weird parameters
  > g <- randomDAG(3,0)
  > wmat <- wgtMatrix(g)
  > if (!all(wmat==matrix(0,3,3))) {
  +   stop("Test of wgtMatrix: Problem when used on empty graph!")
  + }
  > 
  > ## test if weight parameters are correct
  > set.seed(34)
  > g <- randomDAG(5,0.8)
  > trMat <- matrix(0, 5,5)
  > trMat[1,5] <- 0.305
  > trMat[1,4] <- 0.863
  > trMat[1,2] <- 0.354
  > trMat[2,4] <- 0.392
  > trMat[2,5] <- 0.495
  > trMat[3,4] <- 0.278
  > 
  > if (!all(round(wgtMatrix(g),3) == t(trMat))) {
  +   stop("Test in wgtMatrix: Weights have wrong value!")
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.37    0.20    1.56