diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..b137d00
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+pcalg-Ex.R
+src/*.o
+src/*.so
+src/symbols.rds
+*~
+
diff --git a/NAMESPACE b/NAMESPACE
index f9cefed..4899863 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -105,6 +105,7 @@ export(trueCov,
binCItest,
ida,
idaFast,
+ idaFastOpt,
legal.path,
plotSG,
pc.cons.intern,
diff --git a/R/pcalg.R b/R/pcalg.R
index 4f03b75..e17945d 100644
--- a/R/pcalg.R
+++ b/R/pcalg.R
@@ -3207,6 +3207,35 @@ idaFast <- function(x.pos, y.pos.set, mcov, graphEst)
beta.hat
}
+ ## Purpose: Estimate the causal effect of x.pos.set on each element in the
+ ## graph using the local method; graphEst and correlation matrix
+ ## have to be precomputed; orient
+ ## undirected edges at x in a way so that no new collider is
+ ## introduced; if there is an undirected edge between x and y, both directions are considered;
+ ## i.e., y might be partent of x in which case the effect is 0.
+ ## ----------------------------------------------------------------------
+ ## Arguments:
+ ## - x.pos.set:
+ ## - mcov: Covariance matrix that was used to estimate graphEst
+ ## - graphEst: Fit of PC Algorithm (semidirected)
+ ## ----------------------------------------------------------------------
+ ## Value: list of causal values; one list element (matrix) for each element of
+ ## x.pos.set
+ ## ----------------------------------------------------------------------
+ ## Author: Paula Perez, Date: 17 Jan 2017
+
+idaFastOpt <- function(x.pos.set,mcov,graphEst)
+{
+ if(missing(x.pos.set)){
+ x.pos.set = 1:ncol(mcov)
+ }
+ adjmat <- as.matrix(wgtMatrix(graphEst))
+ .Call("idaFastC",(x.pos.set - 1) ,mcov,adjmat)
+}
+
+
+
+
## -> ../man/legal.path.Rd
## only called in qreach() with only 'c' varying
legal.path <- function(a,b,c, amat)
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..bdc1a2a
--- /dev/null
+++ b/README.md
@@ -0,0 +1,94 @@
+# Pcalg fork
+
+## Description
+
+ This is a fork from the `R` package `pcalg`
+that implements a `C++` optimization of the
+function idaFast, idaFastOpt. Using `C++` together
+with the `Armadillo` library instead of `R` improves the
+performance already. On top of that, the new implementation
+seems to have reduced the complexity of the algorithm (see
+Details and Benchmarks).
+
+* **Patch file**: `dataOpt/patch_pcalg.txt`
+
+
+### Details
+
+Given a *completed partially directed acyclic graph* `pdag`, and
+a given node `x`, we identify the unique parents `uniq_pa` and the
+ambiguous parents `amb_pa`. From the set of `amb_pa`, we direct them
+in all possible ways (`2``|amb_pa|`) checking first that
+they introduce no new colliders. This implementation differs with
+the original implementation in how "all possible ways" are examined.
+
+* **Original code:** A loop from 1 to the number of elements
+in `amb_pa`, `|amb_pa|` is executed, and the function,
+
+ ``pa.tmp <- combn(amb_pa,i,simplify = TRUE)``
+
+ is called where `i` is the variable in the loop. All sets of size `i` in `amb_pa`
+are obtained. For every set, all parents in it are directed towards
+`x`, the ones not in it outwards. If the configuration does not introduce a new
+collider, the causal effects are computed, nothing is done otherwise.
+
+* **Optimization:** A loop from 0 until `2``|amb_pa| - 1` is executed,
+where binary representation of `i` (loop variable, `long int`) is sweeping
+over the half posibilities. For every `i` we create two subsets
+`pa_1` and `pa_2` complementary to each other. `pa_1` (`pa_2`) contains the `jth` element
+of `amb_pa` if the `jth` bit of `i` is 1 (0). In that way, we have covered all
+possibilities if we first set `pa_1` to be directed towards `x`, `pa_2` outwards,
+check for colliders, compute causal effects and then do the same reversing the roles
+of `pa_1` and `pa_2.`
+
+
+
+### Benchmarks
+
+ Comparison of performance of `idaFast` (original function of the package)
+versus our implementation, `idaFastOpt`. We have used real data and created **cpdag** of
+100,200,500,1000,5000 nodes. `idaFast` and `idaFastOpt` have then been run on these
+graphs on a Intel(R) Xeon(R) CPU E5-4620 0 @ 2.20GHz. In the graph below, we can see
+the performance of the two functions, in a log-log scale.
+
+We see not only that the runtime is significantly reduced, but also, and
+for the **cpdag** we were testing the functions on, `idaFast` seemed to scale like
+`O(n``3``)` whereas `idaFastOpt` like `O(n``2``)`. Note that
+the scaling is not generalizable to all types of graphs. For different input graphs,
+one can expect a different complexity. Nevertheless, it seems that the new implementation
+has reduced it.
+
+
+ 
+
+## Usage
+
+``idaFastOpt(x.pos.set = 1:N, mcov = cov(subDat), graphEst = l$cpdag_graph)``
+
+
+* **Arguments:**
+ - `x.pos.set`: (integer vector) optional argument.
+ Nodes on which the causal effects are going to
+ be computed. If not passed, the function computes
+ the causal effects on all nodes. Note a vector of positions
+ of the target variables is not passed as an argument, in contrast
+ to `idaFast` since in this function all nodes are going to be
+ targets.
+ - `mcov`: covariance matrix that was used to estimate `graphEst`
+ - `graphEst`: estimated **cpdag** from the function `pc`. If the output of
+ `pc` is `pc.fit`, then the estimated **cpdag** can be obtained by
+ `pc.fit@graph`.
+
+* **Output:**
+ list of causal values; one list element (matrix) for each element of
+ `x.pos.set`
+
+## Tests
+
+ In the folder `dataOpt`, one can find a code example in `run_idaFast.R`,
+ where both functions, `idaFast`, `idaFastOpt`, are run for the **cpdag**'s
+ contained in `dataOpt/skeletonfiles`, and a profiling is performed.
+
+
+ In the same folder, the script `benchmarks.R` creates the plot shown in this README.md
+ file.
diff --git a/dataOpt/01-005-00_Ergebnis_Uebersicht.csv b/dataOpt/01-005-00_Ergebnis_Uebersicht.csv
new file mode 100644
index 0000000..76b3b81
--- /dev/null
+++ b/dataOpt/01-005-00_Ergebnis_Uebersicht.csv
@@ -0,0 +1,101 @@
+Sample Batch CD40L BCR IGF CpG IL10 Myc Number Viability S SubG1 size total RNA
+1 1 1 0 0 0 0 1 8,4 94,9 47,6 2,7 604 343,63
+2 1 0 1 0 0 0 0 5,4 94,2 12,3 10,2 462 106,82
+3 1 0 0 0 0 0 0 5,6 93,0 4,4 3,9 433 133,7
+4 1 0 0 0 0 0 1 8,4 94,6 43,0 4,0 577 226,45
+5 1 0 1 1 1 1 1 8,0 97,4 63,0 2,8 638 356,71
+6 1 1 1 1 1 1 1 8,8 96,4 59,2 2,3 656 383,64
+7 1 1 1 0 1 0 1 7,7 94,5 54,2 2,7 655 399,47
+8 1 1 0 0 1 1 0 8,3 95,2 45,3 1,9 586 188,31
+9 1 0,2 1 1 0 1 0 6,8 92,3 24,4 2,2 557 164,12
+10 1 0 0 0 1 1 0 6,3 95,9 41,4 2,2 574 225,71
+11 2 0 1 0 0 0 1 8,4 93,9 42,4 5,3 543 311,88
+12 2 0 0 1 0 0 0 6,0 93,9 7,0 5,7 479 103,66
+13 2 0 0 0 0 0 0 5,6 92,1 6,2 6,8 461 131,53
+14 2 0 0 0 0 0 1 8,6 95,0 37,0 5,4 550 277,7
+15 2 0 1 0,2 0,2 0 1 8,5 95,9 42,2 4,1 528 364,36
+16 2 1 0,2 0,2 0 0 1 8,0 95,3 38,5 5,9 497 341,34
+17 2 1 1 1 1 0,2 1 8,6 96,3 57,9 4,1 564 420,51
+18 2 0,2 0,2 0,2 0,2 0,2 0 6,5 95,6 28,0 6,6 526 178,34
+19 2 1 0,2 0 0 1 0 5,9 92,0 20,3 4,8 541 185,92
+20 2 0,2 1 0,2 0 0,2 0 5,6 93,5 16,5 9,0 509 157,58
+21 3 0 0 1 0 0 1 9,4 93,1 37,9 4,2 512 246,07
+22 3 0 0 0 1 0 0 6,2 95,4 24,9 2,2 594 176,59
+23 3 0 0 0 0 0 0 6,2 92,8 4,9 2,5 515 135,29
+24 3 0 0 0 0 0 1 8,3 92,9 36,6 3,7 612 232,25
+25 3 0,2 0,2 0,2 0 0 1 8,9 95,2 38,1 3,2 707 300,09
+26 3 0,2 1 0 0 1 1 7,8 91,6 43,5 1,9 658 306,57
+27 3 1 0 0 1 1 1 8,4 94,3 54,9 2,3 714 362,05
+28 3 0,2 1 0,2 1 0 0 6,2 92,3 38,2 1,4 638 213,57
+29 3 0 0 0 0,2 1 0 6,4 91,8 28,7 1,5 582 187,55
+30 3 0,2 1 0,2 0,2 0,2 0 6,1 94,0 27,5 2,2 624 189,19
+31 4b 0 0 0 1 0 1 9,6 93,3 42,5 4,0 595 388,83
+32 4b 0 0 0 0 1 0 6,8 93,0 12,8 7,3 473 143,68
+33 4b 0 0 0 0 0 0 6,7 93,7 5,8 7,4 438 129,82
+34 4b 0 0 0 0 0 1 8,0 93,1 29,4 5,2 559 260,61
+35 4b 0,2 1 0,2 0,2 0,2 1 8,7 92,1 37,4 4,1 585 342,67
+36 4b 1 1 0,2 0,2 0 1 7,7 92,7 38,6 4,7 593 316,21
+37 4b 0,2 0 1 0,2 0 1 7,8 91,9 36,4 4,7 561 307,76
+38 4b 1 1 0,2 0,2 0 0 6,6 94,0 20,9 6,1 523 176,84
+39 4b 0,2 1 1 0,2 0 0 6,6 93,2 23,2 7,4 504 183,84
+40 4b 1 0 1 1 0,2 0 7,5 93,3 40,2 3,9 531 226,45
+41 5 0 0 0 0 1 1 8,2 93,5 38,0 5,0 488 300,33
+42 5 1 0 0 0 0 0 6,2 92,8 9,9 4,1 503 132,34
+43 5 0 0 0 0 0 0 6,0 93,0 4,8 4,6 498 128,24
+44 5 0 0 0 0 0 1 7,9 92,6 34,7 5,7 577 244,87
+45 5 0 0 0 1 1 1 9,7 92,3 49,0 3,7 535 269,28
+46 5 0,2 0,2 0 1 1 1 9,9 94,6 53,2 3,5 556 253,46
+47 5 0,2 0 0,2 1 0 1 9,2 93,0 40,7 3,3 582 262,22
+48 5 0,2 0,2 0,2 0 0 0 6,4 94,3 9,1 4,0 517 135,09
+49 5 0,2 1 0 0 1 0 6,5 93,6 16,9 3,7 534 151,39
+50 5 0,2 0 0,2 1 0 0 7,0 93,6 23,0 3,6 561 141,69
+51 6 1 0 0 0 0 1 8,1 96,4 37,8 2,9 555 283,87
+52 6 0 1 0 0 0 0 6,2 94,5 10,3 5,6 480 143,37
+53 6 0 0 0 0 0 0 6,6 95,2 4,3 4,2 449 130,24
+54 6 0 0 0 0 0 1 7,8 93,9 32,2 3,8 563 272,4
+55 6 0 0 0 0,2 1 1 9,5 94,3 49,6 3,1 579 319,45
+56 6 0,2 0,2 0,2 1 1 1 8,7 95,4 54,7 2,7 599 318,45
+57 6 0,2 1 1 0 0 1 9,4 93,2 36,8 4,1 547 241,83
+58 6 1 1 0 1 1 0 8,0 96,4 44,4 2,3 549 192,91
+59 6 1 1 1 1 1 0 7,4 93,2 45,5 2,1 572 225,32
+60 6 0,2 0 1 0,2 1 0 7,3 93,8 36,5 3,0 537 169,41
+61 7 0 1 0 0 0 1 9,0 92,3 31,2 3,1 564 283,97
+62 7 0 0 1 0 0 0 6,2 95,4 3,5 1,8 560 122,47
+63 7 0 0 0 0 0 0 6,2 93,2 3,5 2,9 555 143,99
+64 7 0 0 0 0 0 1 7,7 93,8 26,4 6,0 498 256,14
+65 7 0,2 1 1 0 1 1 10,3 91,5 31,2 4,0 582 223,09
+66 7 1 0,2 0 0 1 1 9,6 93,8 29,7 4,8 591 219,84
+67 7 1 1 0 1 1 1 8,4 93,9 51,6 1,9 610 323,49
+68 7 0 0,2 0 0,2 1 0 6,8 93,4 21,7 1,5 677 153,95
+69 7 1 1 1 1 0,2 0 7,2 95,6 34,1 1,5 662 247,39
+70 7 1 1 0 1 0 0 6,7 94,5 32,4 0,6 642 162,63
+71 8 0 0 1 0 0 1 8,9 94,6 43,5 2,9 623 230,06
+72 8 0 0 0 1 0 0 6,0 93,0 28,4 4,3 564 160,84
+73 8 0 0 0 0 0 0 6,0 94,4 7,7 6,5 504 139,93
+74 8 0 0 0 0 0 1 9,6 97,0 39,1 5,2 610 270,27
+75 8 0 0,2 0 0,2 1 1 9,8 97,0 43,6 4,4 657 295,49
+76 8 0,2 0,2 0,2 0,2 0,2 1 9,4 98,1 43,2 4,0 644 289,11
+77 8 0,2 1 0,2 0 0,2 1 9,0 95,8 46,3 4,5 651 316,74
+78 8 0 1 1 1 1 0 7,0 97,5 42,2 4,9 607 219,33
+79 8 0,2 0,2 0 1 1 0 7,5 96,2 39,0 4,4 633 205,68
+80 8 0,2 0,2 0,2 1 1 0 7,2 95,2 42,7 4,0 623 246,05
+81 9 0 0 0 1 0 1 8,4 94,0 39,1 6,9 611 347,9
+82 9 0 0 0 0 1 0 5,7 95,0 12,7 4,3 514 159,7
+83 9 0 0 0 0 0 0 6,3 95,9 5,0 6,1 476 144
+84 9 0 0 0 0 0 1 7,9 92,3 32,2 5,4 578 271,96
+85 9 0,2 1 1 0,2 0 1 8,8 95,1 29,5 5,4 621 315,52
+86 9 0,2 0,2 1 0 0,2 1 8,7 93,5 35,2 4,6 600 293,6
+87 9 1 0 1 1 0,2 1 8,4 95,5 46,4 4,4 644 319,79
+88 9 0 1 0,2 0,2 0 0 5,7 95,5 27,7 6,2 521 150,34
+89 9 0,2 0,2 1 0 0,2 0 6,8 94,5 17,2 4,2 538 141,63
+90 9 1 0 1 0,2 0 0 6,6 95,3 28,0 3,7 555 155,49
+91 10 0 0 0 0 1 1 8,7 93,3 34,6 7,1 661 254,27
+92 10 1 0 0 0 0 0 6,3 95,5 7,7 4,7 539 129,75
+93 10 0 0 0 0 0 0 6,2 94,5 3,9 8,8 530 142,34
+94 10 0 0 0 0 0 1 8,5 91,9 29,9 7,6 618 217,24
+95 10 1 0 1 0,2 0 1 9,0 93,4 42,0 3,6 686 256,26
+96 10 0,2 1 0,2 1 0 1 8,1 95,7 48,0 4,2 743 295,8
+97 10 0,2 0 1 0,2 1 1 9,1 94,1 47,7 3,5 721 294,99
+98 10 0,2 1 1 0 0 0 6,2 92,8 14,2 9,0 545 185,25
+99 10 0,2 0 1 0,2 0 0 6,5 95,1 25,4 5,0 560 157,29
+100 10 1 0,2 0,2 0 0 0 6,3 92,1 10,3 5,6 540 139,9
diff --git a/dataOpt/benchmarks.R b/dataOpt/benchmarks.R
new file mode 100644
index 0000000..1123aef
--- /dev/null
+++ b/dataOpt/benchmarks.R
@@ -0,0 +1,51 @@
+# This script creates a graph with the timings of idaFast
+timings <- read.csv("timings.csv", stringsAsFactors=FALSE)
+
+
+fast <- timings[which(timings$optimized==TRUE & timings$Myc == "Low"), ]
+slow <- timings[which(timings$optimized==FALSE & timings$Myc == "Low"), ]
+
+
+
+png("timings.png", width=400, height=400, units ="px")
+
+r_slow <- lm(log(slow[,2]) ~ log(slow[,1]) )
+r_fast <- lm(log(fast[,2]) ~ log(fast[,1]) )
+
+
+plot(slow[,1:2],col="red", xlab = "Anzahl von Genen",
+ ylab = " Laufzeit (Sekunden)", pch=16, cex=1.4, xaxt= "n")
+points(fast[,1:2],col="blue",pch= 8, lwd=2)
+
+axis(1,slow[,1],slow[,1])
+legend("topleft",legend = c("urspruengliche Funktion",
+ "optimierte Funktion"), pch = c(16,8),
+ pt.cex = c(1.4,1), pt.lwd=c(1,2), col = c("red","blue"))
+ #lwd= c(1,2), cex=c(1.4,1),col = c("red","blue"))
+
+mtext("200", at =c(250), line = -20.55 )
+
+dev.off()
+
+png("timings_log.png", width=400, height=400, units ="px")
+
+x <- c(1:5000)
+y <- x^(r_slow$coefficients[2])*exp(r_slow$coefficients[1])
+plot(slow[,1:2],col="red", xlab = "Number of genes, n",
+ ylab = " Time, t (sec)", ylim = c(0.01,50000), pch=16, cex=1.4, yaxt= "n", log = "xy")
+lines(x,y,lty=2, col="lightgray",lwd=2)
+x <- c(1:5000)
+y <- x^(r_fast$coefficients[2])*exp(r_fast$coefficients[1])
+lines(x,y,lty=2, col="lightgray",lwd=2)
+points(fast[,1:2],col="blue",pch= 8, lwd=2)
+axis(2,c(0.01,0.1,1,10,100,1000,10000),c("0.01","0.1","1","10","100","1000","10000"), las = 2)
+legend("topleft",legend = c("idaFast",
+ "idaFastOpt"), pch = c(16,8),
+ pt.cex = c(1.4,1), pt.lwd=c(1,2), col = c("red","blue"))
+ #lwd= c(1,2), cex=c(1.4,1),col = c("red","blue"))
+slope <- sprintf("%.02f", r_slow$coefficients[2])
+text(200,10,bquote( 't ~ '~ O(n^.(slope)) ),srt=38)
+slope <- sprintf("%.02f", r_fast$coefficients[2])
+text(200,0.17,bquote( 't ~ '~ O(n^.(slope)) ), srt=25)
+
+dev.off()
diff --git a/dataOpt/patch_pcalg.txt b/dataOpt/patch_pcalg.txt
new file mode 100644
index 0000000..815244b
--- /dev/null
+++ b/dataOpt/patch_pcalg.txt
@@ -0,0 +1,348 @@
+Index: src/cpdag.cpp
+===================================================================
+--- src/cpdag.cpp (revision 0)
++++ src/cpdag.cpp (revision 0)
+@@ -0,0 +1,189 @@
++/* This program is free software: you can redistribute it and/or modify
++* it under the terms of the GNU General Public License as published by
++* the Free Software Foundation, either version 2 of the License, or
++* (at your option) any later version.
++*
++* This program is distributed in the hope that it will be useful,
++* but WITHOUT ANY WARRANTY; without even the implied warranty of
++* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
++* GNU General Public License for more details.
++*
++* You should have received a copy of the GNU General Public License
++* along with this program. If not, see .
++*
++*/
++
++#include "pcalg/cpdag.hpp"
++
++//Constructor
++cpdag::cpdag(arma::sp_imat adjmat_input, arma::mat covmat_input):
++ adjmat(adjmat_input), covmat(covmat_input) {
++ nodes = adjmat.n_rows;
++ skeleton = (adjmat + adjmat.t());
++ unique = (adjmat - adjmat.t());
++ for(arma::uword i=0; i1) skeleton(i,j)=1;
++ if(unique(i,j)<0) unique(i,j)=0;
++ }
++ }
++ ambiguous = adjmat - unique ;
++ // Sparse matrices are sorted by columns, so we
++ // have to take the counterintuitive convention
++ // of working with the transpose matrices.
++ unique_t = unique.t();
++}
++
++
++
++// idaFast as it is in the R function
++arma::mat cpdag::idaFast(uint32_t x ){
++ //std::cout << "idaFast: computation of causal effects.\n" ;
++ // Las sparse matrices estan ordenada por columnas
++ // asi que hemos tomado la traspuesta arriba.
++ arma::mat beta;
++ arma::irowvec uniq_pa, amb_pa;
++ uniq_pa.set_size(nodes);
++ amb_pa.set_size(nodes);
++ arma::uword p = 0;
++ for (uint32_t col=0; col 64){
++ printf ("The number of ambiguous parents is too high");
++ exit (EXIT_FAILURE);
++ }
++ arma::irowvec amb_pa_1, amb_pa_2;%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+%
+
++ arma::uword n = amb_pa.size();
++ uint64_t L = (1 << (n-1));
++ for (uint64_t bits = 0; bits < L; bits++){
++ directArrows(amb_pa_1,amb_pa_2,amb_pa,bits,n);
++ if( !hasNewCollider(uniq_pa,amb_pa_1,amb_pa_2,x) ){
++ beta.insert_rows(cbeta++,getCausal(uniq_pa,amb_pa_1,x));
++ }
++ if( !hasNewCollider(uniq_pa,amb_pa_2,amb_pa_1,x) ){
++ beta.insert_rows(cbeta++,getCausal(uniq_pa,amb_pa_2,x));
++ }
++ }
++ }
++ return (beta.t());
++}
++
++
++bool cpdag::hasNewCollider(arma::irowvec uniq_pa, arma::irowvec amb_pa_t,
++ arma::irowvec amb_pa_f, arma::uword x){
++ arma::irowvec::iterator uniqst, uniqen, amb_tst,amb_ten, amb_fst, amb_fen;
++ uniqst = uniq_pa.begin(); uniqen = uniq_pa.end();
++ amb_tst = amb_pa_t.begin(); amb_ten = amb_pa_t.end();
++ amb_fst = amb_pa_f.begin(); amb_fen = amb_pa_f.end();
++ if ( amb_pa_t.size() > 0 ){
++ // Check whether all amb_pa_t and uniq_pa are connected. If not --> collider.
++ if ( uniq_pa.size() > 0 ){
++ for (arma::irowvec::iterator it1 = uniqst; it1 < uniqen; it1++){
++ for (arma::irowvec::iterator it2=amb_tst; it2 < amb_ten; it2++){
++ if ( skeleton((*it1),(*it2)) == 0 ) return true;
++ }
++ }
++ }
++ // Check whether all amb_pa_t are connected among each other. If not --> collider.
++ if ( amb_pa_t.size() > 1 ){
++ for ( arma::irowvec::iterator it1=amb_tst; it1 < amb_ten; it1++ ){
++ for ( arma::irowvec::iterator it2=it1+1; it2 < amb_ten; it2++ ){
++ if ( skeleton((*it1),(*it2)) == 0 ) return true;
++ }
++ }
++ }
++ }
++ // if any of the uniq parents of amb_pa_f is NOT directly connected to x --> collider.
++ if ( amb_pa_f.size() > 0 ){
++ arma::irowvec papa;
++ papa.set_size(nodes);
++ arma::uword p=0;
++ for (arma::uword i = 0; i < nodes; i++ ){
++ if( i != x ){
++ for (arma::irowvec::iterator it=amb_fst; it < amb_fen; it++){
++ if( unique((*it),i) ){
++ papa[p++] = i; break;
++ }
++ }
++ }
++ }
++ papa.resize(p);
++ if (p > 0){
++ for (arma::irowvec::iterator it=papa.begin(); it < papa.end(); it++){
++ if( skeleton(x,(*it)) == 0 ) return true;
++ }
++ }
++ }
++ return false;
++}
++
++arma::rowvec cpdag::getCausal(arma::irowvec uniq_pa, arma::irowvec amb_pa_t, arma::uword x){
++ arma::uword n_uniq, n_amb, n_all;
++ n_uniq = uniq_pa.size();
++ n_amb = amb_pa_t.size();
++ n_all = n_uniq + n_amb + 1;
++ arma::irowvec all_pa;
++ all_pa.set_size(n_all);
++ all_pa(0) = x;
++ if ( n_uniq > 0) memcpy(all_pa.begin()+ 1, uniq_pa.begin(),n_uniq*sizeof(long));
++ if ( n_amb > 0) memcpy(all_pa.begin()+ n_uniq + 1, amb_pa_t.begin(),n_amb*sizeof(long));
++ arma::mat A,Y,M;
++ A.set_size(n_all,n_all);
++ Y.set_size(n_all,nodes);
++ arma::irowvec::iterator start = all_pa.begin();
++ arma::irowvec::iterator end = all_pa.end();
++ for (arma::irowvec::iterator it1 = start ; it1 < end; it1++){
++ arma::uword row = it1-start;
++ Y.row(row) = covmat.row((*it1));
++ for (arma::irowvec::iterator it2 = start ; it2 < end; it2++){
++ arma::uword col = it2 - start;
++ A(row,col) = covmat((*it1),(*it2));
++ }
++ }
++ M = solve(A,Y);
++ for (arma::irowvec::iterator it1 = start+1 ; it1 < end; it1++){
++ M(0,(*it1)) = 0.0;
++ }
++ return M.row(0);
++}
++
++
++// Divide the vector amb_pa into two vectors, following the partition k
++void cpdag::directArrows(arma::irowvec &pa_1, arma::irowvec &pa_2,
++ arma::irowvec amb_pa, uint64_t bits, arma::uword n){
++ uint8_t i, idx1, idx2;
++ uint8_t k;
++ idx1 = 0; idx2=0;
++ k = __builtin_popcount(bits);
++ pa_1.resize(k);
++ pa_2.resize(n-k);
++ for (i=0; i < n ; i++ ){
++ if ( bits & (1 << i) ){
++ pa_1(idx1++) = amb_pa(i);
++ } else {
++ pa_2(idx2++) = amb_pa(i);
++ }
++ }
++}
++
+Index: src/gies.cpp
+===================================================================
+--- src/gies.cpp (revision 407)
++++ src/gies.cpp (working copy)
+@@ -21,6 +21,7 @@
+ #include "pcalg/constraint.hpp"
+ #include "pcalg/score.hpp"
+ #include "pcalg/greedy.hpp"
++#include "pcalg/cpdag.hpp"
+ #define DEFINE_GLOBAL_DEBUG_STREAM
+ #include "pcalg/gies_debug.hpp"
+
+@@ -645,3 +646,23 @@
+ END_RCPP
+ }
+
++RcppExport SEXP idaFastC(SEXP xpos, SEXP argMcov,
++ SEXP argAdjMatrix)
++{
++ // Exception handling
++ BEGIN_RCPP
++ arma::mat covmat = (Rcpp::as(argMcov));
++ arma::irowvec x = Rcpp::as(xpos);
++ arma::uword N = x.n_elem;
++ arma::mat adjmat_tmp = (Rcpp::as(argAdjMatrix));
++ arma::sp_imat adjmat;
++ adjmat = arma::conv_to::from(adjmat_tmp);
++ cpdag graph_pd(adjmat,covmat);
++ Rcpp::List to_return(N);
++ arma::uword i=0;
++ for (arma::irowvec::iterator it=x.begin(); it < x.end(); it++){
++ to_return[i++] = graph_pd.idaFast(*it);
++ }
++ return to_return;
++ END_RCPP
++}
+Index: NAMESPACE
+===================================================================
+--- NAMESPACE (revision 407)
++++ NAMESPACE (working copy)
+@@ -105,6 +105,7 @@
+ binCItest,
+ ida,
+ idaFast,
++ idaFastOpt,
+ legal.path,
+ plotSG,
+ pc.cons.intern,
+Index: R/pcalg.R
+===================================================================
+--- R/pcalg.R (revision 407)
++++ R/pcalg.R (working copy)
+@@ -3207,6 +3207,15 @@
+ beta.hat
+ }
+
++idaFastOpt <- function(x.pos.set,mcov,graphEst)
++{
++ if(missing(x.pos.set)){
++ x.pos.set = 1:ncol(mcov)
++ }
++ adjmat <- as.matrix(wgtMatrix(graphEst))
++ .Call("idaFastC",(x.pos.set - 1) ,mcov,adjmat)
++}
++
+ ## -> ../man/legal.path.Rd
+ ## only called in qreach() with only 'c' varying
+ legal.path <- function(a,b,c, amat)
+Index: inst/include/pcalg/cpdag.hpp
+===================================================================
+--- inst/include/pcalg/cpdag.hpp (revision 0)
++++ inst/include/pcalg/cpdag.hpp (revision 0)
+@@ -0,0 +1,45 @@
++/* This program is free software: you can redistribute it and/or modify
++* it under the terms of the GNU General Public License as published by
++* the Free Software Foundation, either version 2 of the License, or
++* (at your option) any later version.
++*
++* This program is distributed in the hope that it will be useful,
++* but WITHOUT ANY WARRANTY; without even the implied warranty of
++* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
++* GNU General Public License for more details.
++*
++* You should have received a copy of the GNU General Public License
++* along with this program. If not, see .
++*
++*/
++
++#ifndef DIFFEXP_H
++#define DIFFEXP_H
++
++#include
++#include
++#include
++
++
++class cpdag{
++ private:
++ arma::uword nodes;
++ arma::mat covmat; // Covariance matrix
++ arma::sp_imat adjmat; // Adjacency matrix CPDAG
++ arma::sp_imat skeleton; // Undirected adjacency matrix
++ arma::sp_imat ambiguous; // Ambiguous parents
++ arma::sp_imat unique; // Unique parents
++ arma::sp_imat unique_t; // Unique parents transposed (used for computations)
++ public:
++ cpdag(arma::sp_imat adjmat_input, arma::mat covmat_input);
++ ~cpdag() {}
++ arma::mat idaFast(uint32_t x);
++ bool hasNewCollider(arma::irowvec uniq_pa, arma::irowvec amb_pa_t,
++ arma::irowvec amb_pa_f, arma::uword x);
++ arma::rowvec getCausal(arma::irowvec uniq_pa, arma::irowvec amb_pa_t, arma::uword x);
++ void directArrows(arma::irowvec &pa_1, arma::irowvec &pa_2,
++ arma::irowvec amb_pa, uint64_t k, arma::uword n);
++
++};
++#endif //cpdag.h
++
diff --git a/dataOpt/run_idaFast.R b/dataOpt/run_idaFast.R
new file mode 100644
index 0000000..3cdd3ee
--- /dev/null
+++ b/dataOpt/run_idaFast.R
@@ -0,0 +1,54 @@
+# We run original idaFast on our skeleton files
+require(devtools)
+require(graph)
+load_all("../../pcalg/",
+ quiet=T,reset=T,export_all=T)
+message(date()," >> Rversion: ")
+print(version)
+
+Ngenes <- c(100, 200, 500, 1000, 2000, 5000)
+Ngenes <- c(200)
+#############################################
+
+# Read the design, so that we can separate Myc high from Myc Low
+designFile <- paste0("01-005-00_Ergebnis_Uebersicht.csv")
+descr.mat <- read.csv(designFile,header = TRUE,
+ sep = "\t", dec = ",",row.names=1)
+# get samples with Myc low and Myc high
+myc_samples <- data.frame( High=as.numeric(rownames(descr.mat[descr.mat$Myc==1,])), Low = as.numeric(rownames(descr.mat[descr.mat$Myc==0,])))
+
+
+tiempoOPT <- c()
+tiempo <- c()
+i = 1
+for (Ng in Ngenes){
+ bigVarGenes <- readRDS(paste0("var_genes/Ng",Ng,"alpha0p5.rds"))
+ for (myc in c("Low","High")){
+ skel_file <- paste0("skeleton_files/Ng",
+ Ng,".it001.Myc",myc,".alpha0p5.rds")
+
+ MycHighLow <- bigVarGenes[myc_samples[,myc], ]
+ l <- readRDS(skel_file)
+ subDat <- MycHighLow[l$ind,]
+ Ns <- nrow(MycHighLow) # number of samples
+ subSize <- floor(Ns/2)
+ tiempoOPT <- c(tiempoOPT,Sys.time())
+ message(Sys.time())
+ mysol <- idaFastOpt(mcov=cov(subDat),graphEst=l$cpdag_graph)
+ message(Sys.time())
+ tiempoOPT[i] <- Sys.time() - tiempoOPT[i]
+ tiempo <- c(tiempo,Sys.time())
+ oldsol <- list()
+ for (j in 1:Ng){
+ oldsol[[j]] <- idaFast(j,1:Ng,cov(subDat),l$cpdag_graph)
+ }
+ message(Sys.time())
+ tiempo[i] <- Sys.time() - tiempo[i]
+ cat(Ng," Myc",myc,"idaFastOpt:", tiempoOPT[i],
+ attr(tiempoOPT[i],"units"),"\n")
+ cat(Ng," Myc",myc,"idaFast:", tiempo[i],
+ attr(tiempo[i],"units"),"\n")
+ i <- i + 1
+ }
+}
+
diff --git a/dataOpt/skeleton_files/Ng100.it001.MycHigh.alpha0p5.rds b/dataOpt/skeleton_files/Ng100.it001.MycHigh.alpha0p5.rds
new file mode 100644
index 0000000..c80bcb8
Binary files /dev/null and b/dataOpt/skeleton_files/Ng100.it001.MycHigh.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng100.it001.MycLow.alpha0p5.rds b/dataOpt/skeleton_files/Ng100.it001.MycLow.alpha0p5.rds
new file mode 100644
index 0000000..c137a3b
Binary files /dev/null and b/dataOpt/skeleton_files/Ng100.it001.MycLow.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng1000.it001.MycHigh.alpha0p5.rds b/dataOpt/skeleton_files/Ng1000.it001.MycHigh.alpha0p5.rds
new file mode 100644
index 0000000..00cd76f
Binary files /dev/null and b/dataOpt/skeleton_files/Ng1000.it001.MycHigh.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng1000.it001.MycLow.alpha0p5.rds b/dataOpt/skeleton_files/Ng1000.it001.MycLow.alpha0p5.rds
new file mode 100644
index 0000000..837175d
Binary files /dev/null and b/dataOpt/skeleton_files/Ng1000.it001.MycLow.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng200.it001.MycHigh.alpha0p5.rds b/dataOpt/skeleton_files/Ng200.it001.MycHigh.alpha0p5.rds
new file mode 100644
index 0000000..b75b15d
Binary files /dev/null and b/dataOpt/skeleton_files/Ng200.it001.MycHigh.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng200.it001.MycLow.alpha0p5.rds b/dataOpt/skeleton_files/Ng200.it001.MycLow.alpha0p5.rds
new file mode 100644
index 0000000..732e611
Binary files /dev/null and b/dataOpt/skeleton_files/Ng200.it001.MycLow.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng2000.it001.MycHigh.alpha0p5.rds b/dataOpt/skeleton_files/Ng2000.it001.MycHigh.alpha0p5.rds
new file mode 100644
index 0000000..6b5b29e
Binary files /dev/null and b/dataOpt/skeleton_files/Ng2000.it001.MycHigh.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng2000.it001.MycLow.alpha0p5.rds b/dataOpt/skeleton_files/Ng2000.it001.MycLow.alpha0p5.rds
new file mode 100644
index 0000000..e222380
Binary files /dev/null and b/dataOpt/skeleton_files/Ng2000.it001.MycLow.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng500.it001.MycHigh.alpha0p5.rds b/dataOpt/skeleton_files/Ng500.it001.MycHigh.alpha0p5.rds
new file mode 100644
index 0000000..bab4766
Binary files /dev/null and b/dataOpt/skeleton_files/Ng500.it001.MycHigh.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng500.it001.MycLow.alpha0p5.rds b/dataOpt/skeleton_files/Ng500.it001.MycLow.alpha0p5.rds
new file mode 100644
index 0000000..13d8d01
Binary files /dev/null and b/dataOpt/skeleton_files/Ng500.it001.MycLow.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng5000.it001.MycHigh.alpha0p5.rds b/dataOpt/skeleton_files/Ng5000.it001.MycHigh.alpha0p5.rds
new file mode 100644
index 0000000..6e27996
Binary files /dev/null and b/dataOpt/skeleton_files/Ng5000.it001.MycHigh.alpha0p5.rds differ
diff --git a/dataOpt/skeleton_files/Ng5000.it001.MycLow.alpha0p5.rds b/dataOpt/skeleton_files/Ng5000.it001.MycLow.alpha0p5.rds
new file mode 100644
index 0000000..ab192e2
Binary files /dev/null and b/dataOpt/skeleton_files/Ng5000.it001.MycLow.alpha0p5.rds differ
diff --git a/dataOpt/timings_log.png b/dataOpt/timings_log.png
new file mode 100644
index 0000000..fd6ada1
Binary files /dev/null and b/dataOpt/timings_log.png differ
diff --git a/dataOpt/timings_log_Low.png b/dataOpt/timings_log_Low.png
new file mode 100644
index 0000000..3f132f2
Binary files /dev/null and b/dataOpt/timings_log_Low.png differ
diff --git a/dataOpt/var_genes/Ng1000alpha0p5.rds b/dataOpt/var_genes/Ng1000alpha0p5.rds
new file mode 100644
index 0000000..61b65ea
Binary files /dev/null and b/dataOpt/var_genes/Ng1000alpha0p5.rds differ
diff --git a/dataOpt/var_genes/Ng100alpha0p5.rds b/dataOpt/var_genes/Ng100alpha0p5.rds
new file mode 100644
index 0000000..551df35
Binary files /dev/null and b/dataOpt/var_genes/Ng100alpha0p5.rds differ
diff --git a/dataOpt/var_genes/Ng2000alpha0p5.rds b/dataOpt/var_genes/Ng2000alpha0p5.rds
new file mode 100644
index 0000000..ada5329
Binary files /dev/null and b/dataOpt/var_genes/Ng2000alpha0p5.rds differ
diff --git a/dataOpt/var_genes/Ng200alpha0p5.rds b/dataOpt/var_genes/Ng200alpha0p5.rds
new file mode 100644
index 0000000..71e8201
Binary files /dev/null and b/dataOpt/var_genes/Ng200alpha0p5.rds differ
diff --git a/dataOpt/var_genes/Ng5000alpha0p5.rds b/dataOpt/var_genes/Ng5000alpha0p5.rds
new file mode 100644
index 0000000..0bd804f
Binary files /dev/null and b/dataOpt/var_genes/Ng5000alpha0p5.rds differ
diff --git a/dataOpt/var_genes/Ng500alpha0p5.rds b/dataOpt/var_genes/Ng500alpha0p5.rds
new file mode 100644
index 0000000..d8a830e
Binary files /dev/null and b/dataOpt/var_genes/Ng500alpha0p5.rds differ
diff --git a/inst/include/pcalg/cpdag.hpp b/inst/include/pcalg/cpdag.hpp
new file mode 100644
index 0000000..ffb76cb
--- /dev/null
+++ b/inst/include/pcalg/cpdag.hpp
@@ -0,0 +1,45 @@
+/* This program is free software: you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 2 of the License, or
+* (at your option) any later version.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see .
+*
+*/
+
+#ifndef DIFFEXP_H
+#define DIFFEXP_H
+
+#include
+#include
+#include
+
+
+class cpdag{
+ private:
+ arma::uword nodes;
+ arma::mat covmat; // Covariance matrix
+ arma::sp_imat adjmat; // Adjacency matrix CPDAG
+ arma::sp_imat skeleton; // Undirected adjacency matrix
+ arma::sp_imat ambiguous; // Ambiguous parents
+ arma::sp_imat unique; // Unique parents
+ arma::sp_imat unique_t; // Unique parents transposed (used for computations)
+ public:
+ cpdag(arma::sp_imat adjmat_input, arma::mat covmat_input);
+ ~cpdag() {}
+ arma::mat idaFast(uint32_t x);
+ bool hasNewCollider(arma::irowvec uniq_pa, arma::irowvec amb_pa_t,
+ arma::irowvec amb_pa_f, arma::uword x);
+ arma::rowvec getCausal(arma::irowvec uniq_pa, arma::irowvec amb_pa_t, arma::uword x);
+ void directArrows(arma::irowvec &pa_1, arma::irowvec &pa_2,
+ arma::irowvec amb_pa, uint64_t k, arma::uword n);
+
+};
+#endif //cpdag.h
+
diff --git a/src/cpdag.cpp b/src/cpdag.cpp
new file mode 100644
index 0000000..ab4aef1
--- /dev/null
+++ b/src/cpdag.cpp
@@ -0,0 +1,191 @@
+/* This program is free software: you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 2 of the License, or
+* (at your option) any later version.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see .
+*
+*/
+
+#include "pcalg/cpdag.hpp"
+
+//Constructor
+cpdag::cpdag(arma::sp_imat adjmat_input, arma::mat covmat_input):
+ adjmat(adjmat_input), covmat(covmat_input) {
+ nodes = adjmat.n_rows;
+ skeleton = (adjmat + adjmat.t());
+ unique = (adjmat - adjmat.t());
+ for(arma::uword i=0; i1) skeleton(i,j)=1;
+ if(unique(i,j)<0) unique(i,j)=0;
+ }
+ }
+ ambiguous = adjmat - unique ;
+ // Sparse matrices are sorted by columns, so we
+ // have to take the counterintuitive convention
+ // of working with the transpose matrices.
+ unique_t = unique.t();
+}
+
+
+
+// idaFast as it is in the R function
+arma::mat cpdag::idaFast(uint32_t x ){
+ //std::cout << "idaFast: computation of causal effects.\n" ;
+ // Las sparse matrices estan ordenada por columnas
+ // asi que hemos tomado la traspuesta arriba.
+ arma::mat beta;
+ arma::irowvec uniq_pa, amb_pa;
+ uniq_pa.set_size(nodes);
+ amb_pa.set_size(nodes);
+ arma::uword p = 0;
+ // Count the number of unambiguous parents
+ for (uint32_t col=0; col 64){
+ printf ("The number of ambiguous parents is too high");
+ exit (EXIT_FAILURE);
+ }
+ arma::irowvec amb_pa_1, amb_pa_2;
+ arma::uword n = amb_pa.size();
+ uint64_t L = (1 << (n-1));
+ for (uint64_t bits = 0; bits < L; bits++){
+ directArrows(amb_pa_1,amb_pa_2,amb_pa,bits,n);
+ if( !hasNewCollider(uniq_pa,amb_pa_1,amb_pa_2,x) ){
+ beta.insert_rows(cbeta++,getCausal(uniq_pa,amb_pa_1,x));
+ }
+ if( !hasNewCollider(uniq_pa,amb_pa_2,amb_pa_1,x) ){
+ beta.insert_rows(cbeta++,getCausal(uniq_pa,amb_pa_2,x));
+ }
+ }
+ }
+ return (beta.t());
+}
+
+
+bool cpdag::hasNewCollider(arma::irowvec uniq_pa, arma::irowvec amb_pa_t,
+ arma::irowvec amb_pa_f, arma::uword x){
+ arma::irowvec::iterator uniqst, uniqen, amb_tst,amb_ten, amb_fst, amb_fen;
+ uniqst = uniq_pa.begin(); uniqen = uniq_pa.end();
+ amb_tst = amb_pa_t.begin(); amb_ten = amb_pa_t.end();
+ amb_fst = amb_pa_f.begin(); amb_fen = amb_pa_f.end();
+ if ( amb_pa_t.size() > 0 ){
+ // Check whether all amb_pa_t and uniq_pa are connected. If not --> collider.
+ if ( uniq_pa.size() > 0 ){
+ for (arma::irowvec::iterator it1 = uniqst; it1 < uniqen; it1++){
+ for (arma::irowvec::iterator it2=amb_tst; it2 < amb_ten; it2++){
+ if ( skeleton((*it1),(*it2)) == 0 ) return true;
+ }
+ }
+ }
+ // Check whether all amb_pa_t are connected among each other. If not --> collider.
+ if ( amb_pa_t.size() > 1 ){
+ for ( arma::irowvec::iterator it1=amb_tst; it1 < amb_ten; it1++ ){
+ for ( arma::irowvec::iterator it2=it1+1; it2 < amb_ten; it2++ ){
+ if ( skeleton((*it1),(*it2)) == 0 ) return true;
+ }
+ }
+ }
+ }
+ // if any of the uniq parents of amb_pa_f is NOT directly connected to x --> collider.
+ if ( amb_pa_f.size() > 0 ){
+ arma::irowvec papa;
+ papa.set_size(nodes);
+ arma::uword p=0;
+ for (arma::uword i = 0; i < nodes; i++ ){
+ if( i != x ){
+ for (arma::irowvec::iterator it=amb_fst; it < amb_fen; it++){
+ if( unique((*it),i) ){
+ papa[p++] = i; break;
+ }
+ }
+ }
+ }
+ papa.resize(p);
+ if (p > 0){
+ for (arma::irowvec::iterator it=papa.begin(); it < papa.end(); it++){
+ if( skeleton(x,(*it)) == 0 ) return true;
+ }
+ }
+ }
+ return false;
+}
+
+arma::rowvec cpdag::getCausal(arma::irowvec uniq_pa, arma::irowvec amb_pa_t, arma::uword x){
+ arma::uword n_uniq, n_amb, n_all;
+ n_uniq = uniq_pa.size();
+ n_amb = amb_pa_t.size();
+ n_all = n_uniq + n_amb + 1;
+ arma::irowvec all_pa;
+ all_pa.set_size(n_all);
+ all_pa(0) = x;
+ if ( n_uniq > 0) memcpy(all_pa.begin()+ 1, uniq_pa.begin(),n_uniq*sizeof(long));
+ if ( n_amb > 0) memcpy(all_pa.begin()+ n_uniq + 1, amb_pa_t.begin(),n_amb*sizeof(long));
+ arma::mat A,Y,M;
+ A.set_size(n_all,n_all);
+ Y.set_size(n_all,nodes);
+ arma::irowvec::iterator start = all_pa.begin();
+ arma::irowvec::iterator end = all_pa.end();
+ for (arma::irowvec::iterator it1 = start ; it1 < end; it1++){
+ arma::uword row = it1-start;
+ Y.row(row) = covmat.row((*it1));
+ for (arma::irowvec::iterator it2 = start ; it2 < end; it2++){
+ arma::uword col = it2 - start;
+ A(row,col) = covmat((*it1),(*it2));
+ }
+ }
+ M = solve(A,Y);
+ for (arma::irowvec::iterator it1 = start+1 ; it1 < end; it1++){
+ M(0,(*it1)) = 0.0;
+ }
+ return M.row(0);
+}
+
+
+// Divide the vector amb_pa into two vectors, following the partition k
+void cpdag::directArrows(arma::irowvec &pa_1, arma::irowvec &pa_2,
+ arma::irowvec amb_pa, uint64_t bits, arma::uword n){
+ uint8_t i, idx1, idx2;
+ uint8_t k;
+ idx1 = 0; idx2=0;
+ k = __builtin_popcount(bits);
+ pa_1.resize(k);
+ pa_2.resize(n-k);
+ for (i = 0; i < n ; i++ ){
+ if ( bits & (1 << i) ){
+ pa_1(idx1++) = amb_pa(i);
+ } else {
+ pa_2(idx2++) = amb_pa(i);
+ }
+ }
+}
+
diff --git a/src/gies.cpp b/src/gies.cpp
index f2cae50..08b8589 100644
--- a/src/gies.cpp
+++ b/src/gies.cpp
@@ -21,6 +21,7 @@ typedef boost::adjacency_list Undi
#include "pcalg/constraint.hpp"
#include "pcalg/score.hpp"
#include "pcalg/greedy.hpp"
+#include "pcalg/cpdag.hpp"
#define DEFINE_GLOBAL_DEBUG_STREAM
#include "pcalg/gies_debug.hpp"
@@ -645,3 +646,27 @@ RcppExport SEXP estimateSkeleton(
END_RCPP
}
+
+RcppExport SEXP idaFastC(SEXP xpos, SEXP argMcov,
+ SEXP argAdjMatrix)
+{
+ // Exception handling
+ BEGIN_RCPP
+ arma::mat covmat = (Rcpp::as(argMcov));
+ arma::irowvec x = Rcpp::as(xpos);
+ arma::uword N = x.n_elem;
+ arma::mat adjmat_tmp = (Rcpp::as(argAdjMatrix));
+ arma::sp_imat adjmat;
+ adjmat = arma::conv_to::from(adjmat_tmp);
+ cpdag graph_pd(adjmat,covmat);
+ Rcpp::List to_return(N);
+ arma::uword i=0;
+ for (arma::irowvec::iterator it=x.begin(); it < x.end(); it++){
+ to_return[i++] = graph_pd.idaFast(*it);
+ }
+ return to_return;
+ END_RCPP
+}
+
+
+