@@ -1,8 +1,8 @@

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 

-Version: 2.5.0

-Date: 2016-09-22

+Version: 2.5.1

+Date: 2016-11-12

 Authors@R: c(person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),

 		     email = "rdiaz02@gmail.com"),

 	      person("Mark", "Taylor", role = "ctb", email = "ningkiling@gmail.com"))

@@ -9,7 +9,8 @@ export("oncoSimulPop", "oncoSimulIndiv", "samplePop",

        "evalAllGenotypesFitAndMut",

        "rfitness",

        "plotFitnessLandscape",

-       "to_Magellan"

+       "to_Magellan",

+       "sampledGenotypes"

        )

 S3method(plot, oncosimul)

@@ -47,7 +48,8 @@ importFrom("dplyr", "full_join", "left_join", "right_join", "%>%", "mutate",

            "filter")

 importFrom("smatr", "ma") ## for major axis regression in some tests

 importFrom("car", "linearHypothesis")

-

+## importFrom("slam", "simple_triplet_zero_matrix", ## "colapply_simple_triplet_matrix",

+##            "col_sums")

@@ -63,6 +63,8 @@ oncoSimulSample <- function(Nindiv,

                             typeSample = "whole",

                             thresholdWhole = 0.5,

                             initMutant = NULL,

+                            AND_DrvProbExit = FALSE,

+                            fixation = NULL,

                             verbosity  = 1,

                             showProgress = FALSE,

                             seed = "auto"){

@@ -190,7 +192,9 @@ oncoSimulSample <- function(Nindiv,

                                initMutant = initMutant,

                                keepPhylog = keepPhylog,

                                mutationPropGrowth = mutationPropGrowth,

-                               detectionProb = detectionProb)        

+                               detectionProb = detectionProb,

+                               AND_DrvProbExit = AND_DrvProbExit,

+                               fixation = fixation)        

         if(tmp$other$UnrecoverExcept) {

             return(f.out.unrecover.except(tmp))

         }

@@ -291,7 +295,7 @@ samplePop <- function(x, timeSample = "last",

                             popSizeSample = popSizeSample)

         dim(z) <- c(1, length(z))

         if(is.null(gN) && (!is.null(x$geneNames)))

-            gN <- geneNames

+            gN <- x$geneNames

     }

     message("\n Subjects by Genes matrix of ",

         nrow(z), " subjects and ",

@@ -338,6 +342,8 @@ oncoSimulPop <- function(Nindiv,

                          errorHitWallTime = TRUE,

                          errorHitMaxTries = TRUE,

                          initMutant = NULL,

+                         AND_DrvProbExit = FALSE,

+                         fixation = NULL,

                          verbosity  = 0,

                          mc.cores = detectCores(),

                          seed = "auto") {

@@ -379,7 +385,9 @@ oncoSimulPop <- function(Nindiv,

                         seed = seed, keepPhylog = keepPhylog,

                         initMutant = initMutant,

                         mutationPropGrowth = mutationPropGrowth,

-                        detectionProb = detectionProb),

+                        detectionProb = detectionProb,

+                        AND_DrvProbExit = AND_DrvProbExit,

+                        fixation = fixation),

                     mc.cores = mc.cores

                     )

     class(pop) <- "oncosimulpop"

@@ -422,16 +430,33 @@ oncoSimulIndiv <- function(fp,

                            errorHitMaxTries = TRUE,

                            verbosity = 0,

                            initMutant = NULL,

+                           AND_DrvProbExit = FALSE,

+                           fixation = NULL,

                            seed = NULL

                            ) {

     call <- match.call()

-    if(all(c(is.na(detectionProb),

-             is.na(detectionSize),

-             is.na(detectionDrivers),

-             is.na(finalTime))))

+    if(all(c(is_null_na(detectionProb),

+             is_null_na(detectionSize),

+             is_null_na(detectionDrivers),

+             is_null_na(finalTime),

+             is_null_na(fixation)

+             )))

         stop("At least one stopping condition should be given.",

              " At least one of detectionProb, detectionSize, detectionDrivers,",

-             " finalTime. Otherwise, we'll run forever.")

+             " finalTime. Otherwise, we'll run until aborted by max.wall.time,",

+             " max.num.tries, and the like.")

+

+    if(AND_DrvProbExit && (is_null_na(detectionProb) || is_null_na(detectionDrivers)))

+        stop("AND_DrvProbExit is TRUE: both of detectionProb and detectionDrivers",

+             " must be non NA.")

+    if(AND_DrvProbExit && !is_null_na(detectionSize)) {

+        warning("With AND_DrvProbExit = TRUE, detectionSize is ignored.")

+        detectionSize <- NA

+    }

+    if(inherits(fp, "fitnessEffects")) {

+        s <- sh <- NULL ## force it soon!

+    }

+

     ## legacies from poor name choices

     typeFitness <- switch(model,

                           "Bozic" = "bozic1",

@@ -467,7 +492,8 @@ oncoSimulIndiv <- function(fp,

         if(model %in% c("Bozic", "Exp") )

             minDetectDrvCloneSz <- 0

         else if (model %in% c("McFL", "McFarlandLog"))

-            minDetectDrvCloneSz <- eFinalMf(initSize, s, detectionDrivers)

+            minDetectDrvCloneSz <- initSize

+        ## minDetectDrvCloneSz <- eFinalMf(initSize, s, detectionDrivers)

         else

             stop("Unknown model")

     }

@@ -484,7 +510,7 @@ oncoSimulIndiv <- function(fp,

     ## No user-visible magic numbers

     ## if(is.null(keepEvery))

     ##     keepEvery <- -9

-    if(is.na(keepEvery)) keepEvery <- -9

+    if(is_null_na(keepEvery)) keepEvery <- -9

     if( (keepEvery > 0) & (keepEvery < sampleEvery)) {

@@ -499,10 +525,11 @@ oncoSimulIndiv <- function(fp,

     if( (typeFitness == "exp") && (death != 1) )

         warning("Using fitness exp with death != 1")

-

-    if(is.na(detectionDrivers)) detectionDrivers <- (2^31) - 1

-    if(is.na(detectionSize)) detectionSize <- Inf

-    if(is.na(finalTime)) finalTime <- Inf

+    if(!is_null_na(detectionDrivers) && (detectionDrivers >= 1e9))

+        stop("detectionDrivers > 1e9; this doesn't seem reasonable")

+    if(is_null_na(detectionDrivers)) detectionDrivers <- (2^31) - 1

+    if(is_null_na(detectionSize)) detectionSize <- Inf

+    if(is_null_na(finalTime)) finalTime <- Inf

     if(!inherits(fp, "fitnessEffects")) {

@@ -515,12 +542,18 @@ oncoSimulIndiv <- function(fp,

             stop(m)

         }

+        if(AND_DrvProbExit) {

+            stop("The AND_DrvProbExit = TRUE setting is invalid",

+                 " with the old poset format.")

+        }

         if(!is.null(muEF))

-            stop("Mutator effects cannot be especified with the old poset format")

+            stop("Mutator effects cannot be specified with the old poset format.")

         if( length(initMutant) > 0)  

             warning("With the old poset format you can no longer use initMutant.",

                     " The initMutant you passed will be ignored.")

-            ## stop("With the old poset, initMutant can only take a single value.")

+        ## stop("With the old poset, initMutant can only take a single value.")

+        if(!is_null_na(fixation))

+            stop("'fixation' cannot be specified with the old poset format.")

         ## Seeding C++ is now much better in new version

         if(is.null(seed) || (seed == "auto")) {## passing a null creates a random seed

             ## name is a legacy. This is really the seed for the C++ generator.

@@ -530,7 +563,7 @@ oncoSimulIndiv <- function(fp,

         if(verbosity >= 2)

             cat(paste("\n Using ", seed, " as seed for C++ generator\n"))

-        if(!is.na(detectionProb)) stop("detectionProb cannot be used in v.1 objects")

+        if(!is_null_na(detectionProb)) stop("detectionProb cannot be used in v.1 objects")

         ## if(message.v1)

         ##     message("You are using the old poset format. Consider using the new one.")

@@ -576,6 +609,7 @@ oncoSimulIndiv <- function(fp,

                   silent = !verbosity)

         objClass <- "oncosimul"

     } else {

+        s <- sh <- NULL ## force it.

         if(numPassengers != 0)

             warning(paste("Specifying numPassengers has no effect",

                           " when using fitnessEffects objects. ",

@@ -594,6 +628,18 @@ oncoSimulIndiv <- function(fp,

             if(verbosity >= 2)

                 cat("\n A (high quality) random seed will be generated in C++\n")

         }

+        if(!is_null_na(fixation)) {

+            if( (!is.list(fixation)) && (!is.vector(fixation))  )

+                stop("'fixation' must be a list or a vector.")

+            if(!(all(unlist(lapply(fixation, is.vector)))))

+                stop("Each element of 'fixation' must be a single element character vector.")

+            if(!(all(unlist(lapply(fixation, class)) == "character")))

+                stop("Each element of 'fixation' must be a single element character vector.")

+            if(!(all( unlist(lapply(fixation, length)) == 1)))

+                stop("Each element of 'fixation' must be a single element character vector.")

+            if(AND_DrvProbExit)

+                stop("It makes no sense to pass AND_DrvProbExit and a fixation list.")

+        }

         op <- try(nr_oncoSimul.internal(rFE = fp, 

                                         birth = birth,

                                         death = death,  

@@ -625,7 +671,9 @@ oncoSimulIndiv <- function(fp,

                                         errorHitMaxTries = errorHitMaxTries,

                                         keepPhylog = keepPhylog,

                                         MMUEF = muEF,

-                                        detectionProb = detectionProb),

+                                        detectionProb = detectionProb,

+                                        AND_DrvProbExit = AND_DrvProbExit,

+                                        fixation = fixation),

                   silent = !verbosity)

         objClass <- c("oncosimul", "oncosimul2")

     }

@@ -666,7 +714,7 @@ summary.oncosimul <- function(object, ...) {

         if( (tmp$minDMratio == -99)) tmp$minDMratio <- NA

         if( (tmp$minBMratio == -99)) tmp$minBMratio <- NA

         tmp$OccurringDrivers <- object$OccurringDrivers

-        return(as.data.frame(tmp))

+        return(as.data.frame(tmp, stringsAsFactors = FALSE))

     }

 }

@@ -690,10 +738,29 @@ print.oncosimul <- function(x, ...) {

 }

 ## I want this to return things that are storable

+## summary.oncosimulpop <- function(object, ...) {

+##     as.data.frame(rbindlist(lapply(object, summary)))

+## }

+

 summary.oncosimulpop <- function(object, ...) {

-    as.data.frame(rbindlist(lapply(object, summary)))

+    tmp <- lapply(object, summary)

+    rm <- which(unlist(lapply(tmp, function(x) (length(x) == 1) && (is.na(x)))))

+    if(length(rm) > 0)

+        if(length(rm) < length(object)) {

+        warning("Some simulations seem to have failed and will be removed",

+                " from the summary. The failed runs are ",

+                paste(rm, collapse = ", "),

+                ".")

+        tmp <- tmp[-rm]

+        } else {

+            warning("All simulations failed.")

+            return(NA)

+        }

+    as.data.frame(rbindlist(tmp))

 }

+

+

 print.oncosimulpop <- function(x, ...) {

     cat("\nPopulation of OncoSimul trajectories of",

         length(x), "individuals. Call :\n")

@@ -1383,8 +1450,17 @@ phylogClone <- function(x, N = 1, t = "last", keepEvents = TRUE) {

         stop("Phylogenetic information is only stored with v >= 2")

     z <- which_N_at_T(x, N, t)

     tG <- x$GenotypesLabels[z] ## only for GenotypesLabels we keep all

-                               ## sample size info at each period

+    ## sample size info at each period

+

+    if( (length(tG) == 1) && (tG == "")) {

+        warning("There never was a descendant of WT")

+    }

+    

     df <- x$other$PhylogDF

+    if(nrow(df) == 0) {

+        warning("PhylogDF has 0 rows: no descendants of initMutant ever appeared.")

+        return(NA)

+    }

     if(!keepEvents) { ## is this just a graphical thing? or not?

         df <- df[!duplicated(df[, c(1, 2)]), ]

     }

@@ -1420,6 +1496,13 @@ plotClonePhylog <- function(x, N = 1, t = "last",

              "very fast, before any clones beyond the initial were ",

              "generated.")

     pc <- phylogClone(x, N, t, keepEvents)

+    ## if(is.na(pc)) {

+    ##     ## This should not be reachable, as caught before

+    ##     ## where we check for nrow of PhylogDF   

+    ##     warning("No clone phylogeny available. Exiting without plotting.")

+    ##     return(NULL)

+    ## }

+        

     l0 <- igraph::layout.reingold.tilford(pc$g)

     if(!timeEvents) {

         plot(pc$g, layout = l0)

@@ -1502,6 +1585,12 @@ get.the.time.for.sample <- function(tmp, timeSample, popSizeSample) {

 get.mut.vector <- function(x, timeSample, typeSample,

                            thresholdWhole, popSizeSample) {

+    if(is.null(x$TotalPopSize)) {

+        warning(paste("It looks like this simulation never completed.",

+                      " Maybe it reached maximum allowed limits.",

+                      " You will get NAs"))

+        return(rep(NA, length(x$geneNames)))

+    }

     the.time <- get.the.time.for.sample(x, timeSample, popSizeSample)

     if(the.time < 0) { 

         return(rep(NA, nrow(x$Genotypes)))

@@ -1686,17 +1775,6 @@ oncoSimul.internal <- function(poset, ## restrict.table,

 }

-eFinalMf <- function(initSize, s, j) {

-    ## Expected final sizes for McF, when K is set to the default.

-    # j is number of drivers

-    ## as it says, with no passengers

-    ## Set B(d) = D(N)

-    K <- initSize/(exp(1) - 1)

-    return(K * (exp( (1 + s)^j) - 1))

-}

-

-

-

 OncoSimulWide2Long <- function(x) {

     ## Put data in long format, for ggplot et al

@@ -1820,6 +1898,59 @@ plotShannon <- function(z) {

     axis(2)

 }

+

+

+is_null_na <- function(x) {

+    ## For arguments, if user passes "undefined" or "not set"

+    ## See also http://stackoverflow.com/a/19655909

+    if(is.function(x)) return(FALSE)

+    if( is.null(x) ||

+        ( (length(x) == 1) && (is.na(x)) ) ||

+        ( (length(x) == 1) && (x == "") ) ## careful here

+       )  {

+        return(TRUE)

+    } else {

+        return(FALSE)

+    }

+}

+

+

+## Not used anymore, but left here in case they become useful.

+## Expected numbers at equilibrium under McFarland's

+## eFinalMf <- function(initSize, s, j) {

+##     ## Expected final sizes for McF, when K is set to the default.

+##     # j is number of drivers

+##     ## as it says, with no passengers

+##     ## Set B(d) = D(N)

+##     K <- initSize/(exp(1) - 1)

+##     return(K * (exp( (1 + s)^j) - 1))

+## }

+

+## mcflE <- function(p, s, initSize) {

+##     K <- initSize/(exp(1) - 1)

+##     ## Expected number at equilibrium

+##     return( K * (exp((1 + s)^p) - 1))

+## }

+

+## mcflEv <- function(p, s, initSize) {

+##     ## expects vectors for p and s

+##     K <- initSize/(exp(1) - 1)

+    

+##     ## Expected number at equilibrium

+##     return( K * (exp(prod((1 + s)^p)) - 1))

+## }

+

+

+

+

+

+

+

+

+

+

+

+

 ## simpsonI <- function(x) {

 ##     sx <- sum(x)

 ##     p <- x/sx

@@ -1,8 +1,8 @@

 # This file was generated by Rcpp::compileAttributes

 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

-nr_BNB_Algo5 <- function(rFE, mu_, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery) {

-    .Call('OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu_, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery)

+nr_BNB_Algo5 <- function(rFE, mu_, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list) {

+    .Call('OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu_, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list)

 }

 BNB_Algo5 <- function(restrictTable, numDrivers, numGenes, typeCBN_, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime) {

@@ -17,6 +17,10 @@ evalRGenotypeAndMut <- function(rG, rFE, muEF, full2mutator_, verbose, prodNeg)

     .Call('OncoSimulR_evalRGenotypeAndMut', PACKAGE = 'OncoSimulR', rG, rFE, muEF, full2mutator_, verbose, prodNeg)

 }

+accessibleGenotypes <- function(y, f, numMut, th) {

+    .Call('OncoSimulR_accessibleGenotypes', PACKAGE = 'OncoSimulR', y, f, numMut, th)

+}

+

 ## readFitnessEffects <- function(rFE, echo) {

 ##     invisible(.Call('OncoSimulR_readFitnessEffects', PACKAGE = 'OncoSimulR', rFE, echo))

 ## }

@@ -192,16 +192,23 @@ plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-

 ##                 tfm = tfm))

 ## }

-count_accessible_g <- function(gfm, accessible_th) {

-    gaj <- genot_to_adj_mat(gfm)

-    gaj <- filter_inaccessible(gaj, accessible_th)

-    return(ncol(gaj) - 1)

-}

+## No longer being used. Used to be in rfitness

+## count_accessible_g <- function(gfm, accessible_th) {

+##     gaj <- genot_to_adj_mat(gfm)

+##     gaj <- filter_inaccessible(gaj, accessible_th)

+##     return(ncol(gaj) - 1)

+## }

+

+## There is now C++ code to get just the locations/positions of the

+## accessible genotypes

 filter_inaccessible <- function(x, accessible_th) {

     ## Return an adjacency matrix with only accessible paths. x is the gaj

     ## matrix created in the plots. The difference between genotypes

     ## separated by a hamming distance of 1

+

+    ## FIXME: could do the x[, -1] before loop and not add the 1

+    ## inside while, and do that at end

     colnames(x) <- rownames(x) <- 1:ncol(x)

     while(TRUE) {

         ## remove first column

@@ -216,6 +223,132 @@ filter_inaccessible <- function(x, accessible_th) {

     return(x)

 }

+f1 <- function() {}

+

+x <- 99

+

+## wrapper to the C++ code

+wrap_accessibleGenotypes <- function(x, th) {

+    ## x is the fitness matrix, not adjacency matrix

+    numMut <- rowSums(x[, -ncol(x)])

+    o_numMut <- order(numMut)

+    x <- x[o_numMut, ]

+    numMut <- numMut[o_numMut]

+    

+    y <- x[, -ncol(x)]

+    storage.mode(y) <- "integer"

+

+    acc <- accessibleGenotypes(y, x[, ncol(x)],

+                               as.integer(numMut),

+                               th)

+    return(acc[acc > 0])

+}

+

+## A transitional function

+faster_accessible_genotypes_R <- function(x, th) {

+   rs0 <- rowSums(x[, -ncol(x)])

+    x <- x[order(rs0), ]

+    rm(rs0)

+    

+    y <- x[, -ncol(x)]

+    f <- x[, ncol(x)]

+    rs <- rowSums(y)

+

+   ## If 0, not accessible

+   ## adm <- slam::simple_triplet_zero_matrix(nrow = length(rs), ncol = length(rs),

+   ##                                          mode = "integer")

+   

+   adm <- matrix(0, nrow = length(rs), ncol = length(rs))

+   storage.mode(adm) <- "integer"

+   

+   ## Most time is gone here

+    for(i in 1:length(rs)) { ## i is the current genotype

+        candidates <- which(rs == (rs[i] + 1))

+        for(j in candidates) {

+            if( (sum(abs(y[j, ] - y[i, ])) == 1) &&

+                ( (f[j] - f[i]) >= th ) ) {

+                ## actually, this is the largest time sink using slam

+                adm[i, j] <- 1L

+                }

+        }

+    }

+

+    colnames(adm) <- rownames(adm) <- 1:ncol(adm)

+    admtmp <- adm[, -1, drop = FALSE] ## we do not want the root column.

+    while(TRUE) {

+        ## We remove inaccessible cols (genotypes) and the corresponding

+        ## rows repeatedly until nothing left to remove; any column left

+        ## is therefore accessible throw at least one path.

+

+        ## inacc_col <- which(slam::colapply_simple_triplet_matrix(admtmp, FUN = sum) == 0L)

+        inacc_col <- which(colSums(admtmp) == 0L)

+        if(length(inacc_col) == 0) break;

+        inacc_row <- inacc_col + 1 ## recall root row is left

+        admtmp <- admtmp[-inacc_row, -inacc_col, drop = FALSE]

+    }

+    return(as.numeric(c(colnames(adm)[1], colnames(admtmp))))

+

+}

+

+

+## ## This uses slam, but that is actually slower because

+## ## of the assignment

+## faster_accessible_genots_slam <- function(x, th = 0) {

+

+##     ## Given a genotype matrix, return the genotypes that are accessible

+##     ## via creating a directed adjacency matrix between genotypes

+##     ## connected (i.e., those that differ by gaining one mutation). 0

+##     ## means not connected, 1 means connected.

+    

+##     ## There is a more general function in OncoSimulR that will give the

+##     ## fitness difference. But not doing the difference is faster than

+##     ## just setting a value, say 1, if all we want is to keep track of

+##     ## accessible ones. And by using only 0/1 we can store only an

+##     ## integer. And no na.omits, etc. Is too restricted? Yes. But for

+##     ## simulations and computing just accessible genotypes, probably a

+##     ## hell of a lot faster.

+

+##     ## Well, this is not incredibly fast either.

+    

+##     ## Make sure sorted, so ancestors always before descendants

+##     rs0 <- rowSums(x[, -ncol(x)])

+##     x <- x[order(rs0), ]

+##     rm(rs0)

+    

+##     y <- x[, -ncol(x)]

+##     f <- x[, ncol(x)]

+##     rs <- rowSums(y)

+

+##     ## If 0, not accessible

+##     adm <- slam::simple_triplet_zero_matrix(nrow = length(rs), ncol = length(rs),

+##                                       mode = "integer")

+##     for(i in 1:length(rs)) { ## i is the current genotype

+##         candidates <- which(rs == (rs[i] + 1))

+##         for(j in candidates) {

+##             ## sumdiff <- sum(abs(y[j, ] - y[i, ]))

+##             ## if(sumdiff == 1)

+##             if( (sum(abs(y[j, ] - y[i, ])) == 1) &&

+##                 ( (f[j] - f[i]) >= th ) )

+##                 adm[i, j] <- 1L

+##         }

+##     }

+

+##     colnames(adm) <- rownames(adm) <- 1:ncol(adm)

+##     admtmp <- adm[, -1, drop = FALSE] ## we do not want the root column.

+##     while(TRUE) {

+##         ## We remove inaccessible cols (genotypes) and the corresponding

+##         ## rows repeatedly until nothing left to remove; any column left

+##         ## is therefore accessible throw at least one path.

+

+##         ## inacc_col <- which(slam::colapply_simple_triplet_matrix(admtmp, FUN = sum) == 0L)

+##         inacc_col <- which(slam::col_sums(admtmp) == 0L)

+##         if(length(inacc_col) == 0) break;

+##         inacc_row <- inacc_col + 1 ## recall root row is left

+##         admtmp <- admtmp[-inacc_row, -inacc_col, drop = FALSE]

+##     }

+##     return(as.numeric(c(colnames(adm)[1], colnames(admtmp))))

+## }

+

 generate_matrix_genotypes <- function(g) {

     ## FIXME future: do this for order too? Only if rfitness for order.

@@ -257,6 +390,8 @@ genot_to_adj_mat <- function(x) {

     ## that differ by gaining one mutation) with value being the

     ## difference in fitness between destination and origin

+    ## FIXME: code is now in place to do all of this in C++

+    

     ## Make sure sorted, so ancestors always before descendants

     rs0 <- rowSums(x[, -ncol(x)])

     x <- x[order(rs0), ]

@@ -23,14 +23,22 @@

 to_Magellan <- function(x, file,

                         max_num_genotypes = 2000) {

+    ## This is stupid if we already have, as entry, an object from

+    ## rfitness!! to_Fitness_Matrix can be very slow.

     if(is.null(file)) {

         file <- tempfile()

         cat("\n Using file ", file, "\n")

     }

-    gfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)$gfm

-    write(rep(2, ncol(gfm) - 1), file = file, ncolumns = ncol(gfm) - 1)

-    write.table(gfm, file = file, append = TRUE,

-                row.names = FALSE, col.names = FALSE, sep = " ")

+    if(inherits(x, "genotype_fitness_matrix")) {

+        write(rep(2, ncol(x) - 1), file = file, ncolumns = ncol(x) - 1)

+        write.table(x, file = file, append = TRUE,

+                    row.names = FALSE, col.names = FALSE, sep = " ")

+    } else {

+        gfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)$gfm

+        write(rep(2, ncol(gfm) - 1), file = file, ncolumns = ncol(gfm) - 1)

+        write.table(gfm, file = file, append = TRUE,

+                    row.names = FALSE, col.names = FALSE, sep = " ")

+    }

 }

 to_Fitness_Matrix <- function(x, max_num_genotypes) {

@@ -1,3 +1,4 @@

+

 ## Copyright 2013, 2014, 2015, 2016 Ramon Diaz-Uriarte

 ## This program is free software: you can redistribute it and/or modify

@@ -161,7 +162,7 @@ list.of.deps <- function(x) {

             stop("child not unique")

     }

     typeDep <- lookupTypeDep[unique(x$typeDep)]

-    if(any(is.na(typeDep)))

+    if(any(is_null_na(typeDep)))

         stop("typeDep value incorrect. Check spelling.")

     return(list(

         child = unique(x$child),

@@ -1436,7 +1437,9 @@ nr_oncoSimul.internal <- function(rFE,

                                   extraTime,

                                   keepPhylog,

                                   detectionProb,

-                                  MMUEF = NULL ## avoid partial matching, and set default

+                                  AND_DrvProbExit,

+                                  MMUEF = NULL, ## avoid partial matching, and set default

+                                  fixation = NULL ## avoid partial matching

                                   ) {

     if(!inherits(rFE, "fitnessEffects"))

         stop(paste("rFE must be an object of class fitnessEffects",

@@ -1578,6 +1581,28 @@ nr_oncoSimul.internal <- function(rFE,

     ## if( is.null(n2)) n2 <- -9

     ## call <- match.call()

+    

+    ## Process the fixed list, if any

+    if(!is_null_na(fixation)) {

+        ng <- namedGenes

+        rownames(ng) <- namedGenes[, "Gene"]

+        ## Usual genotype specification and might allow ordered vectors

+        ## in the future

+        fixation_b <- lapply(fixation, nice.vector.eo, sep = ",")

+        ulf <- unlist(fixation_b)

+        if(any(ulf == ">"))

+            stop("Order effects not allowed (yet?) in fixation.")

+        ulfg <- ng[ulf, 1]

+        if(any(is.na(ulfg)))

+            stop(paste("The 'fixation' list contains genes that are not present",

+                       " in the fitness effects."))

+        ## Sorting here is crucial!!

+        fixation_list <- lapply(fixation_b, function(x) sort(ng[x, 2]))

+    } else {

+        fixation_list <- list()

+    }

+

+    

     return(c(

         nr_BNB_Algo5(rFE = rFE,

                      mu_ = mu,

@@ -1613,7 +1638,9 @@ nr_oncoSimul.internal <- function(rFE,

                      p2 = dpr["p2"],

                      PDBaseline = dpr["PDBaseline"],

                      cPDetect_i= dpr["cPDetect"],

-                     checkSizePEvery = dpr["checkSizePEvery"]),

+                     checkSizePEvery = dpr["checkSizePEvery"],

+                     AND_DrvProbExit = AND_DrvProbExit,

+                     fixation_list = fixation_list),

         Drivers = list(rFE$drv), ## but when doing pops, these will be repeated

         geneNames = list(names(getNamesID(rFE)))

     ))

@@ -1753,8 +1780,8 @@ detectionProbCheckParse <- function(x, initSize, verbosity) {

     default_PDBaseline <- 1.2 * initSize

     default_checkSizePEvery <- 20

     ## No default cPDetect. That is done from p2 and n2 in C++.

-    

-    if((length(x) == 1) && (is.na(x))) {

+    if(is_null_na(x)) {

+    ## if((length(x) == 1) && (is.na(x))) {

         y <- vector()

         y["cPDetect"] <- -9

         y["p2"] <- 9

@@ -1779,18 +1806,18 @@ detectionProbCheckParse <- function(x, initSize, verbosity) {

     }

     ## This ain't conditional. If not given, always check

-    if( !is.na(x["cPDetect"]) && (sum(!is.na(x["p2"]), !is.na(x["n2"])) >= 1 ))

+    if( !is_null_na(x["cPDetect"]) && (sum(!is_null_na(x["p2"]), !is_null_na(x["n2"])) >= 1 ))

         stop("Specify only cPDetect xor both of p2 and n2")

-    if( (is.na(x["p2"]) + is.na(x["n2"])) == 1 )

+    if( (is_null_na(x["p2"]) + is_null_na(x["n2"])) == 1 )

         stop("If you pass one of n2 or p2, you must also pass the other. ",

              "Otherwise, we would not know what to do.")

-    if(is.na(x["PDBaseline"])) {

+    if(is_null_na(x["PDBaseline"])) {

         x["PDBaseline"] <- default_PDBaseline

         if(verbosity > -1)

             message("\n  PDBaseline set to default value of ", default_PDBaseline, "\n")

         }

-    if(is.na(x["checkSizePEvery"])) {

+    if(is_null_na(x["checkSizePEvery"])) {

         x["checkSizePEvery"] <- default_checkSizePEvery

         if(verbosity > -1)

             message("\n  checkSizePEvery set to default value of ",

@@ -1800,9 +1827,9 @@ detectionProbCheckParse <- function(x, initSize, verbosity) {

     ## If we get here, we checked consistency of whether cPDetect or p2/n2.

     ## Fill up with defaults the missing values

-    if(is.na(x["cPDetect"])) {

-        if(is.na(x["p2"])) {

-            if(!is.na(x["n2"])) stop("Eh? no p2 but n2? Bug")

+    if(is_null_na(x["cPDetect"])) {

+        if(is_null_na(x["p2"])) {

+            if(!is_null_na(x["n2"])) stop("Eh? no p2 but n2? Bug")

             x["n2"] <- default_n2

             x["p2"] <- default_p2

             if(verbosity > -1)

@@ -1816,14 +1843,14 @@ detectionProbCheckParse <- function(x, initSize, verbosity) {

     if(x["PDBaseline"] < 0)

         stop("PDBaseline < 0")

-    if(!is.na(x["n2"])) {

+    if(!is_null_na(x["n2"])) {

         if(x["n2"] <= x["PDBaseline"])

             stop("n2 <= PDBaseline")

         if(x["p2"] >= 1) stop("p2 >= 1")

         if(x["p2"] <= 0) stop("p2 <= 0")

         x["cPDetect"] <- -9

     } else { ## only if x["cPDetect"] is not NA

-        if(is.na(x["cPDetect"])) stop("eh? you found a bug")## paranoia

+        if(is_null_na(x["cPDetect"])) stop("eh? you found a bug")## paranoia

         x["n2"] <- -9

         x["p2"] <- -9

         if(verbosity > -1)

@@ -1832,6 +1859,51 @@ detectionProbCheckParse <- function(x, initSize, verbosity) {

     return(x)

 }

+sampledGenotypes <- function(y, genes = NULL) {

+    ## From a popSample object, or a matrix for that matter,

+    ## show the sampled genotypes and their frequencies

+    if(!is.null(genes)) {

+        cols <- which(colnames(y) %in% genes )

+        y <- y[, cols]

+    }

+    nn <- colnames(y)

+    df <- data.frame(table(

+        apply(y, 1, function(z) paste(nn[as.logical(z)], collapse = ", ") )

+    ))

+    gn <- as.character(df[, 1])

+    gn[gn == ""] <- "WT"

+    df <- data.frame(Genotype = gn, Freq = df[, 2], stringsAsFactors = FALSE)

+    return(df)

+}

+

+

+

+list_g_matches_fixed <- function(x, y) {

+    ## Internal function, for testing the fixation output.

+    ## x and y are vectors

+

+    ## x is the set of output genotypes, y the set of fixed

+    ## genotypes/subset of genotypes.

+

+    ## Yes, this function has tests too in test.fixation.R

+    

+    ## All genotypes in x satisfy that they are supersets of at least one

+    ## in y? That is true if, for every element in x, at least one y in

+    ## that x.

+

+    if(is.list(y)) y <- unlist(y)

+    

+    y.nice <- lapply(y, nice.vector.eo, sep = ",")

+    x.nice <- lapply(x, nice.vector.eo, sep = ",")

+

+    fu <- function(u, y.nice)

+        any(unlist(lapply(y.nice, function(z) all(z %in% u))))

+

+    return(all(unlist(lapply(x.nice, fu, y.nice))))

+

+}

+

+

 ## emptyFitnessEffects <- function() {

 ##     list(long.rt = list(),

@@ -61,7 +61,9 @@ rfitness <- function(g, c= 0.5,

         }

         m <- cbind(m, Fitness = fi)

         if(min_accessible_genotypes > 0) {

-            num_accessible_genotypes <- count_accessible_g(m, accessible_th)

+            ## num_accessible_genotypes <- count_accessible_g(m, accessible_th)

+            ## Recall accessibleGenotypes includes the wt, if accessible.

+            num_accessible_genotypes <- length(wrap_accessibleGenotypes(m, accessible_th)) - 1

             if(num_accessible_genotypes >= min_accessible_genotypes) {

                 done <- TRUE

                 attributes(m) <- c(attributes(m),

@@ -89,3 +91,7 @@ create_eq_ref <- function(g) {

     ref <- c(rep(1, nm), rep(0, g - nm))

     sample(ref)

 }

+

+

+

+

@@ -1,3 +1,16 @@

+Changes in version 2.5.1 (2016-11-12):

+	- AND of detectedSizeP and lastMaxDr.

+	- fixation as stopping mechanism.

+	- sampledGenotypes in user code.

+	- clonePhylog et al: deal with never any descendant.

+	- samplePop can handle failed simulations graciously.

+	- summary.oncosimulpop can handle failed simulations graciously.

+	- accessible genotypes now done in C++.

+	- OcurringDrivers should not be a factor.

+	- samplePop always returns gene names.

+	- to_Magellan is much faster with rfitness objects.

+	- Several improvements in vignette (English and additional explanations).

+

 Changes in version 2.4.0 (for BioC 3.4):

 	- Mutator phenotype and gene-specific mutation rates.

 	- End simulations stochastically as a function of size.	

@@ -9,7 +22,7 @@ Changes in version 2.4.0 (for BioC 3.4):

 	- Improved test coverage.

 	- samplePop: sample at arbitrary sizes.	

 	- evalAllGenotypes: order = FALSE by default.

-

+	

 Changes in version 2.3.17 (2017-09-22):

 	- random2 for rfitness.

 	- Vignette: decrease size and running time.

@@ -114,7 +114,7 @@ allMutatorEffects(epistasis = NULL, noIntGenes = NULL,

   specifiy anything if you do not want to, and you can pass an empty

   vector (as \code{character(0)}). The default has changed with respect

   to v.2.1.3 and previous: it used to be to assume that all

-  genes that were not in the \code{noIntGenes} were drivers. The fault

+  genes that were not in the \code{noIntGenes} were drivers. The default

   now is to assume nothing: if you want \code{drvNames} you have

   to specify them.

@@ -40,6 +40,8 @@

                 errorHitMaxTries = TRUE,

                 verbosity = 0,

                 initMutant = NULL,

+                AND_DrvProbExit = FALSE,

+                fixation = NULL,

                 seed = NULL)

 oncoSimulPop(Nindiv, fp, model = "Exp", numPassengers = 0, mu = 1e-6,

@@ -61,6 +63,8 @@ oncoSimulPop(Nindiv, fp, model = "Exp", numPassengers = 0, mu = 1e-6,

                 errorHitWallTime = TRUE,

                 errorHitMaxTries = TRUE,

                 initMutant = NULL,

+                AND_DrvProbExit = FALSE,

+                fixation = NULL,

                 verbosity = 0,

                 mc.cores = detectCores(),

                 seed = "auto")

@@ -108,6 +112,8 @@ oncoSimulSample(Nindiv,

                 typeSample = "whole",

                 thresholdWhole = 0.5,

                 initMutant = NULL,

+                AND_DrvProbExit = FALSE,

+                fixation = NULL,

                 verbosity  = 1,

                 showProgress = FALSE,

                 seed = "auto")

@@ -289,22 +295,7 @@ This option can not be used with v.1 objects.

 }

-\item{s}{

-  Selection coefficient for drivers. 

-    Only relevant if using a poset as this is included in the

-  fitnessEffects object.

-}

-\item{sh}{

-  Selection coefficient for drivers with restrictions not satisfied. A