@@ -50,7 +50,7 @@ oncoSimulSample <- function(Nindiv,

                             minDetectDrvCloneSz = "auto",

                             extraTime = 0,

                             finalTime = 0.25 * 25 * 365,

-                            onlyCancer = TRUE,

+                            onlyCancer = FALSE,

                             keepPhylog = FALSE,

                             mutationPropGrowth = ifelse(model == "Bozic",

                                                         FALSE, TRUE),

@@ -375,7 +375,7 @@ oncoSimulPop <- function(Nindiv,

                          ## ifelse(model \%in\% c("Bozic", "Exp"), -9,

                          ##                       5 * sampleEvery),

                          finalTime = 0.25 * 25 * 365,

-                         onlyCancer = TRUE,

+                         onlyCancer = FALSE,

                          keepPhylog = FALSE,

                          mutationPropGrowth = ifelse(model == "Bozic",

                                                      FALSE, TRUE),

@@ -468,7 +468,7 @@ oncoSimulIndiv <- function(fp,

                            minDetectDrvCloneSz = "auto",

                            extraTime = 0,

                            finalTime = 0.25 * 25 * 365,

-                           onlyCancer = TRUE,

+                           onlyCancer = FALSE,

                            keepPhylog = FALSE,

                            mutationPropGrowth = ifelse(model == "Bozic",

                                                        FALSE, TRUE),

@@ -710,9 +710,12 @@ summary.oncosimul <- function(object, ...) {

     ## if those are not regarded as errors

     pbp <- ("pops.by.time" %in% names(object) )

-    if(object$other$UnrecoverExcept) { ## yes, when bailing out from

-                                     ## except. can have just minimal

-                                     ## content

+    if(object$other$UnrecoverExcept) {

+        ## yes, when bailing out from

+        ## except. can have just minimal

+        ## content

+        message("An unrecoverable exception happened in ",

+                "the simulation of the input object.")

         return(NA)

     } else if( !pbp & (object$HittedWallTime || object$HittedMaxTries) ) {

         return(NA)

@@ -928,7 +931,20 @@ plot.oncosimul <- function(x,

                            ...

                            ) {

-

+    if (x$other$UnrecoverExcept) {

+        stop("An unrecoverable exception happened in ",

+             "the simulation of the input object.",

+             "Empty object with nothing to plot.") 

+    } else {

+        pbp <- ("pops.by.time" %in% names(x) )

+        if ( !pbp & (x$HittedWallTime || x$HittedMaxTries) ) {

+            stop("The simulation hit max wall time or max tries. ",

+                 "Empty object with nothing to plot.")

+        }  else if ( !pbp & !(x$HittedWallTime || x$HittedMaxTries) ) {

+            stop("Eh? No pops.by.time but did not hit wall time or max tries? BUG!")

+        }

+    }

+    

     if(!(type %in% c("stacked", "stream", "line")))

         stop("Type of plot unknown: it must be one of",

              "stacked, stream or line")

@@ -64,7 +64,10 @@ oncoSimulSample <- function(Nindiv,

                             fixation = NULL,

                             verbosity  = 1,

                             showProgress = FALSE,

-                            seed = "auto"){

+                            seed = "auto",

+                            interventions = NULL,

+                            userVars = NULL,

+                            rules = NULL){

     ## No longer using mclapply, because of the way we use the limit on

     ## the number of tries.

@@ -192,7 +195,10 @@ oncoSimulSample <- function(Nindiv,

                                mutationPropGrowth = mutationPropGrowth,

                                detectionProb = detectionProb,

                                AND_DrvProbExit = AND_DrvProbExit,

-                               fixation = fixation)        

+                               fixation = fixation,

+                               interventions = interventions,

+                                userVars = userVars,

+                                rules = rules)        

         if(tmp$other$UnrecoverExcept) {

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

         }

@@ -383,7 +389,10 @@ oncoSimulPop <- function(Nindiv,

                          fixation = NULL,

                          verbosity  = 0,

                          mc.cores = detectCores(),

-                         seed = "auto") {

+                         seed = "auto",

+                         interventions = NULL,

+                         userVars = NULL,

+                         rules = NULL) {

     if(Nindiv < 1)

         stop("Nindiv must be >= 1")

@@ -424,7 +433,10 @@ oncoSimulPop <- function(Nindiv,

                         mutationPropGrowth = mutationPropGrowth,

                         detectionProb = detectionProb,

                         AND_DrvProbExit = AND_DrvProbExit,

-                        fixation = fixation),

+                        fixation = fixation,

+                        interventions = interventions,

+                        userVars = userVars,

+                        rules = rules),

                     mc.cores = mc.cores)

     ## mc.allow.recursive = FALSE ## FIXME: remove?

                     ## done for covr issue

@@ -469,8 +481,10 @@ oncoSimulIndiv <- function(fp,

                            initMutant = NULL,

                            AND_DrvProbExit = FALSE,

                            fixation = NULL,

-                           seed = NULL

-                           ) {

+                           seed = NULL,

+                           interventions = NULL,

+                           userVars = NULL,

+                           rules = NULL) {

     call <- match.call()

     if(all(c(is_null_na(detectionProb),

              is_null_na(detectionSize),

@@ -495,6 +509,9 @@ oncoSimulIndiv <- function(fp,

     }

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

         stop("v.1 functionality has been removed. Please use v.2")

+    

+    if(!inherits(fp, "fitnessEffects_v3")) 

+        fp <- convertFitnessEffects(fp)

     ## legacies from poor name choices

     typeFitness <- switch(model,

@@ -504,6 +521,8 @@ oncoSimulIndiv <- function(fp,

                           "McFL" = "mcfarlandlog",

                           "McFarlandLogD" = "mcfarlandlogd",

                           "McFLD" = "mcfarlandlogd",

+                          "Arb" = "arbitrary",

+                          "Const" = "constant",

                           stop("No valid value for model")

                           )

     if(max(initSize) < 1)

@@ -515,7 +534,21 @@ oncoSimulIndiv <- function(fp,

     }       ##  if ( !(model %in% c("McFL", "McFarlandLog") )) {

             ## K <- 1 ## K is ONLY used for McFarland; set it to 1, to avoid

             ##        ## C++ blowing.

-

+    

+    if("deathSpec" %in% names(fp)) {

+        if (fp$deathSpec) {

+            if (typeFitness != "arbitrary" && typeFitness != "constant") {

+                stop("If death is specified in the fitness effects, use Arb or Const model.")

+            }

+        }

+        

+        else {

+            if (typeFitness == "arbitrary") {

+                stop("To use Arb model specify both birth and death in fitness effects.")

+            }

+        }

+    }

+    

     if(typeFitness == "exp") {

         death <- 1

         ## mutationPropGrowth <- 1

@@ -531,7 +564,7 @@ oncoSimulIndiv <- function(fp,

     }

     if(minDetectDrvCloneSz == "auto") {

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

+        if(model %in% c("Bozic", "Exp", "Arb", "Const") )

             minDetectDrvCloneSz <- 0

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

                               "McFLD", "McFarlandLogD")) {

@@ -610,6 +643,12 @@ oncoSimulIndiv <- function(fp,

             if(AND_DrvProbExit)

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

         }

+

+        #if interventions is null, we create an empty list, cos' it will be easier to handle by C++

+        if(is_null_na(interventions)){

+            interventions = list()

+        }

+

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

                                         birth = birth,

                                         death = death,  

@@ -643,7 +682,10 @@ oncoSimulIndiv <- function(fp,

                                         MMUEF = muEF,

                                         detectionProb = detectionProb,

                                         AND_DrvProbExit = AND_DrvProbExit,

-                                        fixation = fixation),

+                                        fixation = fixation,

+                                        interventions = interventions,

+                                        userVars = userVars,

+                                        rules = rules),

                   silent = !verbosity)

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

     ## }

@@ -1494,7 +1536,11 @@ get.the.time.for.sample <- function(tmp, timeSample, popSizeSample) {

           } else if (length(candidate.time) == 1) {

                 message("Only one sampled time period with mutants.")

                 the.time <- candidate.time

-            } else {

+          } else {

+              ## If we have drivers at t = 2, t= 4, t= 5, ... but not at t=3

+              ## because at t=3 there were no clones with drivers,

+              ## t = 3 is not among the candidate times.

+              

                   the.time <- sample(candidate.time, 1)

               }

     } else {

@@ -57,9 +57,6 @@ oncoSimulSample <- function(Nindiv,

                             max.memory = 2000,

                             max.wall.time.total = 600,

                             max.num.tries.total = 500 * Nindiv,

-                            ## well, obviously they are errors

-                            ## errorHitWallTime = TRUE,

-                            ## errorHitMaxTries = TRUE,

                             typeSample = "whole",

                             thresholdWhole = 0.5,

                             initMutant = NULL,

@@ -458,9 +455,6 @@ oncoSimulIndiv <- function(fp,

                            keepEvery = sampleEvery,

                            minDetectDrvCloneSz = "auto",

                            extraTime = 0,

-                           ## used to be this

-                           ## ifelse(model \%in\% c("Bozic", "Exp"), -9,

-                           ##                     5 * sampleEvery),

                            finalTime = 0.25 * 25 * 365,

                            onlyCancer = TRUE,

                            keepPhylog = FALSE,

@@ -585,84 +579,6 @@ oncoSimulIndiv <- function(fp,

     if(is_null_na(finalTime)) finalTime <- Inf

     if(is_null_na(sampleEvery)) stop("sampleEvery cannot be NULL or NA")

-    

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

-        ## if(any(unlist(lapply(list(fp, 

-        ##                           numPassengers,

-        ##                           s, sh), is.null)))) {

-        ##     m <- paste("You are using the old poset format.",

-        ##                "You must specify all of poset, numPassengers",

-        ##                "s, and sh.")

-        ##     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 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.")

-        ## 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.

-        ##     ## Nope, we cannot use 2^32, because as.integer will fail.

-        ##     seed <- as.integer(round(runif(1, min = 0, max = 2^16)))

-        ## }

-        ## if(verbosity >= 2)

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

-

-        ## 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.")

-   

-    

-        ## ## A simulation stops if cancer or finalTime appear, the first

-        ## ## one. But if we set onlyCnacer = FALSE, we also accept simuls

-        ## ## without cancer (or without anything)

-        

-        ## op <- try(oncoSimul.internal(poset = fp, ## restrict.table = rt,

-        ##                              ## numGenes = numGenes,

-        ##                              numPassengers = numPassengers,

-        ##                              typeCBN = "CBN",

-        ##                              birth = birth,

-        ##                              s = s,

-        ##                              death = death,  

-        ##                              mu =  mu,  

-        ##                              initSize =  initSize, 

-        ##                              sampleEvery =  sampleEvery,  

-        ##                              detectionSize =  detectionSize, 

-        ##                              finalTime = finalTime, 

-        ##                              initSize_species = 2000, 

-        ##                              initSize_iter = 500, 

-        ##                              seed = seed, 

-        ##                              verbosity = verbosity, 

-        ##                              speciesFS = 10000,  

-        ##                              ratioForce = 2,

-        ##                              typeFitness = typeFitness,

-        ##                              max.memory = max.memory,

-        ##                              mutationPropGrowth = mutationPropGrowth,                                   

-        ##                              initMutant = -1, 

-        ##                              max.wall.time = max.wall.time,

-        ##                              max.num.tries = max.num.tries,

-        ##                              keepEvery = keepEvery,  

-        ##                              ## alpha = 0.0015,  

-        ##                              sh = sh,

-        ##                              K = K, 

-        ##                              minDetectDrvCloneSz = minDetectDrvCloneSz,

-        ##                              extraTime = extraTime,

-        ##                              detectionDrivers = detectionDrivers,

-        ##                              onlyCancer = onlyCancer,

-        ##                              errorHitWallTime = errorHitWallTime,

-        ##                              errorHitMaxTries = errorHitMaxTries),

-        ##           silent = !verbosity)

-        ## objClass <- "oncosimul"

-    ## } else {

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

         if(numPassengers != 0)

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

@@ -934,88 +850,6 @@ plot.oncosimulpop <- function(x, ask = TRUE,

 }

-## plot.oncosimul <- function(x, col = c(8, "orange", 6:1),

-##                            log = "y",

-##                            ltyClone = 2:6,

-##                            lwdClone = 0.9,

-##                            ltyDrivers = 1,

-##                            lwdDrivers = 3,

-##                            xlab = "Time units",

-##                            ylab = "Number of cells",

-##                            plotClones = TRUE,

-##                            plotDrivers = TRUE,

-##                            addtot = FALSE,

-##                            addtotlwd = 0.5,

-##                            yl = NULL,

-##                            thinData = FALSE,

-##                            thinData.keep = 0.1,

-##                            thinData.min = 2,

-##                            plotDiversity = FALSE,

-##                            ...

-##                            ) {

-

-##     if(thinData)

-##         x <- thin.pop.data(x, keep = thinData.keep, min.keep = thinData.min)

-

-##     ## uvx

-##     if(!inherits(x, "oncosimul2"))

-##         ndr <- colSums(x$Genotypes[1:x$NumDrivers, , drop = FALSE])

-##     else {

-##         ndr <- colSums(x$Genotypes[x$Drivers, , drop = FALSE])

-##     }

-    

-##     if(is.null(yl)) {

-##         if(log %in% c("y", "xy", "yx") )

-##             yl <- c(1, max(apply(x$pops.by.time[, -1, drop = FALSE], 1, sum)))

-##         else

-##             yl <- c(0, max(apply(x$pops.by.time[, -1, drop = FALSE], 1, sum)))

-##     }

-##     if(plotDiversity) {

-##         par(fig = c(0, 1, 0.8, 1))

-##         m1 <- par()$mar

-##         m <- m1

-##         m[c(1, 3)] <- c(0, 0.7)

-##         op <- par(mar = m )

-##         plotShannon(x)

-##         par(op)

-##         m1[c(3)] <- 0.2

-##         op <- par(mar = m1)

-##         par(fig = c(0, 1, 0, 0.8), new = TRUE)  

-##     }

-##     if(plotClones) {

-##         plotClones(x,

-##                    ndr = ndr, 

-##                    xlab = xlab,

-##                    ylab = ylab,

-##                    lty = ltyClone,

-##                    col = col, 

-##                    ylim = yl,

-##                    lwd = lwdClone,

-##                    axes = FALSE,

-##                    log = log,

-##                    ...)

-##     }

-

-##     if(plotClones && plotDrivers)

-##         par(new = TRUE)

-    

-##     if(plotDrivers){

-##         plotDrivers0(x,

-##                      ndr,

-##                      timescale = 1,

-##                      trim.no.drivers = FALSE,

-##                      xlab = "", ylab = "",

-##                      lwd = lwdDrivers,

-##                      lty = ltyDrivers,

-##                      col = col, 

-##                      addtot = addtot,

-##                      addtotlwd = addtotlwd,

-##                      log = log, ylim = yl,

-##                      ...)

-##     }

-    

-## }

-

 plot.oncosimul <- function(x,

                            show = "drivers", 

@@ -1511,10 +1345,6 @@ plotPoset <- function(x, names = NULL, addroot = FALSE,

         box()

 }

-## this function seems to never be used

-## plotAdjMat <- function(adjmat) {

-##     plot(as(adjmat, "graphNEL"))

-## }

@@ -1591,12 +1421,6 @@ 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) {

@@ -1724,198 +1548,6 @@ get.mut.vector <- function(x, timeSample, typeSample,

 }

-## 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)))

-##     } 

-##     pop <- x$pops.by.time[the.time, -1]

-    

-##     if(all(pop == 0)) {

-##         stop("You found a bug: this should never happen")

-##     }

-    

-##     if(typeSample %in% c("wholeTumor", "whole")) {

-##         popSize <- x$PerSampleStats[the.time, 1]

-##         return( as.numeric((tcrossprod(pop,

-##                                        x$Genotypes)/popSize) >= thresholdWhole) )

-##     } else if (typeSample %in%  c("singleCell", "single")) {

-

-##         return(x$Genotypes[, sample(seq_along(pop), 1, prob = pop)])

-##     } else {

-##         stop("Unknown typeSample option")

-##     }

-## }

-

-

-

-

-

-

-

-

-

-

-

-

-## oncoSimul.internal <- function(poset, ## restrict.table,

-##                                numPassengers, 

-##                                ## numGenes,

-##                                typeCBN,

-##                                birth, 

-##                                s,

-##                                death,

-##                                mu,

-##                                initSize,

-##                                sampleEvery,

-##                                detectionSize,

-##                                finalTime,

-##                                initSize_species,

-##                                initSize_iter,

-##                                seed,

-##                                verbosity,

-##                                speciesFS,

-##                                ratioForce,

-##                                typeFitness,

-##                                max.memory,

-##                                mutationPropGrowth, ## make it explicit

-##                                initMutant,

-##                                max.wall.time,

-##                                keepEvery,

-##                                alpha,

-##                                sh,                               

-##                                K,

-##                                ## endTimeEvery,

-##                                detectionDrivers,

-##                                onlyCancer,

-##                                errorHitWallTime,

-##                                max.num.tries,

-##                                errorHitMaxTries,

-##                                minDetectDrvCloneSz,

-##                                extraTime) {

-

-##     ## the value of 20000, in megabytes, for max.memory sets a limit of ~ 20 GB

-  

-

-##     ## if(keepEvery < sampleEvery)

-##     ##     warning("setting keepEvery to sampleEvery")

-

-##     ## a backdoor to allow passing restrictionTables directly

-##     if(inherits(poset, "restrictionTable"))

-##         restrict.table <- poset

-##     else

-##         restrict.table <- poset.to.restrictTable(poset)

-##     numDrivers <- nrow(restrict.table)

-##     numGenes <- (numDrivers + numPassengers)

-##     if(numGenes < 2)

-##         stop("There must be at least two genes (loci) in the fitness effects.",

-##              "If you only care about a degenerate case with just one,",

-##              "you can enter a second gene",

-##              "with fitness effect of zero.")

-##     if(numGenes > 64)

-##         stop("Largest possible number of genes (loci) is 64 for version 1.",

-##              "You are strongly encouraged to use the new specification",

-##              "as in version 2.")

-

-##     ## These can never be set by the user

-##     ## if(initSize_species < 10) {

-##     ##     warning("initSize_species too small?")

-##     ## }

-##     ## if(initSize_iter < 100) {

-##     ##     warning("initSize_iter too small?")

-##     ## }

-

-##     ## numDrivers <- nrow(restrict.table)

-##     if(length(unique(restrict.table[, 1])) != numDrivers)

-##         stop("BAIL OUT NOW: length(unique(restrict.table[, 1])) != numDrivers)")

-##     ddr <- restrict.table[, 1]

-##     if(any(diff(ddr) != 1))

-##         stop("BAIL OUT NOW:  any(diff(ddr) != 1")

-##     ## sanity checks

-##     if(max(restrict.table[, 1]) != numDrivers)

-##         stop("BAIL OUT NOW: max(restrict.table[, 1]) != numDrivers")

-##     if(numDrivers > numGenes)

-##         stop("BAIL OUT NOW: numDrivers > numGenes")

-    

-##     non.dep.drivers <- restrict.table[which(restrict.table[, 2] == 0), 1]

-

-

-

-

-##     ## if( (is.null(endTimeEvery) || (endTimeEvery > 0)) &&

-##     ##    (typeFitness %in% c("bozic1", "exp") )) {

-##     ##     warning(paste("endTimeEvery will take a positive value. ",

-##     ##                   "This will make simulations not stop until the next ",

-##     ##                   "endTimeEvery has been reached. Thus, in simulations ",

-##     ##                   "with very fast growth, simulations can take a long ",

-##     ##                   "time to finish, or can hit the wall time limit. "))

-##     ## }

-##     ## if(is.null(endTimeEvery))

-##     ##     endTimeEvery <- keepEvery

-##     ## if( (endTimeEvery > 0) && (endTimeEvery %% keepEvery) )

-##     ##     warning("!(endTimeEvery %% keepEvery)")

-##     ## a sanity check in restricTable, so no neg. indices for the positive deps

-##     neg.deps <- function(x) {

-##         ## checks a row of restrict.table

-##         numdeps <- x[2]

-##         if(numdeps) {

-##             deps <- x[3:(3+numdeps - 1)]

-##             return(any(deps < 0))

-##         } else FALSE

-##     }

-##     if(any(apply(restrict.table, 1, neg.deps)))

-##         stop("BAIL OUT NOW: Negative dependencies in restriction table")

-

-##     ## transpose the table

-##     rtC <- convertRestrictTable(restrict.table)

-

-    

-##     return(c(

-##         BNB_Algo5(restrictTable = rtC,

-##         numDrivers = numDrivers,

-##         numGenes = numGenes,

-##         typeCBN_= typeCBN,

-##         s = s, 

-##         death = death,

-##         mu = mu,

-##         initSize = initSize,

-##         sampleEvery = sampleEvery,

-##         detectionSize = detectionSize,

-##         finalTime = finalTime,

-##         initSp = initSize_species,

-##         initIt = initSize_iter,

-##         seed = seed,

-##         verbosity = verbosity,

-##         speciesFS = speciesFS,

-##         ratioForce = ratioForce,

-##         typeFitness_ = typeFitness,

-##         maxram = max.memory,

-##         mutationPropGrowth = as.integer(mutationPropGrowth),

-##         initMutant = initMutant,

-##         maxWallTime = max.wall.time,

-##         keepEvery = keepEvery,

-##         sh = sh,

-##         K = K,

-##         detectionDrivers = detectionDrivers,

-##         onlyCancer = onlyCancer,

-##         errorHitWallTime = errorHitWallTime,

-##         maxNumTries = max.num.tries,

-##         errorHitMaxTries = errorHitMaxTries,

-##         minDetectDrvCloneSz = minDetectDrvCloneSz,

-##         extraTime = extraTime

-##     ),

-##     NumDrivers = numDrivers

-##              ))

-

-## }

-

 OncoSimulWide2Long <- function(x) {

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

@@ -1947,23 +1579,6 @@ OncoSimulWide2Long <- function(x) {

-## We are not using this anymore

-## create.muts.by.time <- function(tmp) { ## tmp is the output from Algorithm5

-##     if(tmp$NumClones > 1) {

-##         NumMutations <- apply(tmp$Genotypes, 2, sum)

-##         muts.by.time <- cbind(tmp$pops.by.time[, c(1), drop = FALSE],

-##                               t(apply(tmp$pops.by.time[, -c(1),

-##                                                        drop = FALSE], 1,

-##                                       function(x) tapply(x,

-##                                                          NumMutations, sum))))

-##         colnames(muts.by.time)[c(1)] <- "Time"

-##     } else {

-##         muts.by.time <- tmp$pops.by.time

-##     }

-##     return(muts.by.time)

-## } 

-

-

 create.drivers.by.time <- function(tmp, ndr) {

     ## CountNumDrivers <- apply(tmp$Genotypes[1:numDrivers, ,drop = FALSE], 2, sum)

     CountNumDrivers <- ndr

@@ -2084,248 +1699,9 @@ is_null_na <- function(x) {

-

-

-

-

-

-

-

-

-

 ## simpsonI <- function(x) {

 ##     sx <- sum(x)

 ##     p <- x/sx

 ##     p <- p[p > 0]

 ##     return(sum(p^2)))

 ## }

-

-## plotSimpson <- function(z) {

-    

-##     h <- apply(z$pops.by.time[, 2:ncol(z$pops.by.time), drop = FALSE],

-##                1, shannonI)

-##     plot(x = z$pops.by.time[, 1],

-##          y = h, lty = "l", xlab = "", ylab = "H")

-## }

-

-

-## plotClones <- function(z, ndr = NULL, na.subs = TRUE,

-##                        log = "y", type = "l",

-##                        lty = 1:8, col = 1:9, ...) {

-

-##     ## if given ndr, we order columns based on ndr, so clones with more

-##     ## drivers are plotted last

-

-##     y <- z$pops.by.time[, 2:ncol(z$pops.by.time), drop = FALSE]

-    

-##     if(na.subs){

-##         y[y == 0] <- NA

-##     }

-##     if(!is.null(ndr)) {

-##         ## could be done above, to avoid creating

-##         ## more copies

-##         oo <- order(ndr)

-##         y <- y[, oo, drop = FALSE]

-##         ndr <- ndr[oo]

-##         col <- col[ndr + 1]

-##     }

-##     matplot(x = z$pops.by.time[, 1],

-##             y = y,

-##             log = log, type = type,

-##             col = col, lty = lty,

-##             ...)

-##     box()

-## }

-

-

-

-

-

-## No longer used

-## rtNoDep <- function(numdrivers) {

-##     ## create a restriction table with no dependencies

-##     x <- matrix(nrow = numdrivers, ncol = 3)

-##     x[, 1] <- 1:numdrivers

-##     x[, 2] <- 0

-##     x[, 3] <- -9

-##     return(x)

-## }

-

-

-## Simulate from generative model. This is a quick function, and is most

-## likely wrong! Never used for anything.

-

-## simposet <- function(poset, p) {

-##     ## if (length(parent.nodes) != length (child.nodes)){

-##     ##     print("An Error Occurred")

-##     ## }

-##     ##    else {

-##     num.genes <- max(poset) - 1 ## as root is not a gene

-##     genotype <-t(c(1, rep(NA, num.genes)))

-##     colnames(genotype) <- as.character(0:num.genes)

-    

-    

-##     poset$runif <- runif(nrow(poset))

-##     ## this.relation.prob.OK could be done outside, but having it inside

-##     ## the loop would allow to use different thresholds for different

-##     ## relationships

-##     for (i in (1:nrow(poset))) {

-##         child <- poset[i, 2]

-##         this.relation.prob.OK <- as.numeric(poset[i, "runif"] > p)

-##         the.parent <- genotype[ poset[i, 1] ] ## it's the value of parent in genotype. 

-##         if (is.na(genotype[child])){

-##             genotype[child] <- this.relation.prob.OK * the.parent  

-##         }

-##         else

-##             genotype[child] <- genotype[child]*(this.relation.prob.OK * the.parent)

-##     }

-##     ##    }

-    

-##     return(genotype)

-## }

-

-

-## to plot and adjacency matrix in this context can do

-## plotPoset(intAdjMatToPoset(adjMat))

-## where intAdjMatToPoset is from best oncotree code: generate-random-trees.

-## No! the above is simpler

-

-

-

-

-## get.mut.vector.whole <- function(tmp, timeSample = "last", threshold = 0.5) {

-##     ## Obtain, from  results from a simulation run, the vector

-##     ## of 0/1 corresponding to each gene.

-    

-##     ## threshold is the min. proportion for a mutation to be detected

-##     ## We are doing whole tumor sampling here, as in Sprouffske

-

-##     ## timeSample: do we sample at end, or at a time point, chosen

-##     ## randomly, from all those with at least one driver?

-    

-##     if(timeSample == "last") {

-##         if(tmp$TotalPopSize == 0)

-##             warning(paste("Final population size is 0.",

-##                           "Thus, there is nothing to sample with ",

-##                           "sampling last. You will get NAs"))

-##         return(as.numeric(

-##             (tcrossprod(tmp$pops.by.time[nrow(tmp$pops.by.time), -1],

-##                         tmp$Genotypes)/tmp$TotalPopSize) > threshold))

-##     } else if (timeSample %in% c("uniform", "unif")) {

-##           candidate.time <- which(tmp$PerSampleStats[, 4] > 0)

-          

-##           if (length(candidate.time) == 0) {

-##               warning(paste("There is not a single sampled time",

-##                             "at which there are any mutants.",

-##                             "Thus, no uniform sampling possible.",

-##                             "You will get NAs"))

-##               return(rep(NA, nrow(tmp$Genotypes)))

-##           } else if (length(candidate.time) == 1) {

-##                 the.time <- candidate.time

-##             } else {

-##                   the.time <- sample(candidate.time, 1)

-##               }

-##           pop <- tmp$pops.by.time[the.time, -1]

-##           popSize <- tmp$PerSampleStats[the.time, 1]

-##           ## if(popSize == 0)

-##           ##     warning(paste("Population size at this time is 0.",

-##           ##                   "Thus, there is nothing to sample at this time point.",

-##           ##                   "You will get NAs"))

-##           return( as.numeric((tcrossprod(pop,

-##                                        tmp$Genotypes)/popSize) > threshold) )

-##       }

-## }

-

-

-

-##           the.time <- sample(which(tmp$PerSampleStats[, 4] > 0), 1)

-##           if(length(the.time) == 0) {

-##               warning(paste("There are no clones with drivers at any time point.",

-##                             "No uniform sampling possible.",

-##                             "You will get a vector of NAs."))

-##             return(rep(NA, nrow(tmp$Genotypes)))  

-##           }

-## get.mut.vector.singlecell <- function(tmp, timeSample = "last") {

-##     ## No threshold, as single cell.

-

-##     ## timeSample: do we sample at end, or at a time point, chosen

-##     ## randomly, from all those with at least one driver?

-    

-##     if(timeSample == "last") {

-##         the.time <- nrow(tmp$pops.by.time)

-##     } else if (timeSample %in% c("uniform", "unif")) {

-##          candidate.time <- which(tmp$PerSampleStats[, 4] > 0)

-         

-##          if (length(candidate.time) == 0) {

-##              warning(paste("There is not a single sampled time",

-##                            "at which there are any mutants.",

-##                            "Thus, no uniform sampling possible.",

-##                            "You will get NAs"))

-##              return(rep(NA, nrow(tmp$Genotypes)))

-##          } else if (length(candidate.time) == 1) {

-##                the.time <- candidate.time

-##            } else {

-##                  the.time <- sample(candidate.time, 1)

-##              }

-

-##      }

-##     pop <- tmp$pops.by.time[the.time, -1]

-##     ##       popSize <- tmp$PerSampleStats[the.time, 1]

-##     ## genot <- sample(seq_along(pop), 1, prob = pop)

-##     if(all(pop == 0)) {

-##         warning(paste("All clones have a population size of 0",

-##                       "at the chosen time. Nothing to sample.",

-##                       "You will get NAs"))

-##         return(rep(NA, nrow(tmp$Genotypes)))

-##     } else {

-##           return(tmp$Genotypes[, sample(seq_along(pop), 1, prob = pop)])

-##       }

-## }

-

-

-## get.mut.vector <- function(x, timeSample = "whole", typeSample = "last",

-##                            thresholdWhole = 0.5) {

-##     if(typeSample %in% c("wholeTumor", "whole")) {

-##         get.mut.vector.whole(x, timeSample = timeSample,

-##                              threshold = thresholdWhole)

-##     } else if(typeSample %in%  c("singleCell", "single")) {

-##         get.mut.vector.singlecell(x, timeSample = timeSample)

-##     }

-## }

-

-

-

-

-

-## plotClonePhylog <- function(x, timeEvent = FALSE,

-##                             showEvents = TRUE,

-##                             fixOverlap = TRUE) {

-##     if(!inherits(x, "oncosimul2"))

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

-##     if(nrow(x$other$PhylogDF) == 0)

-##         stop("It seems you run the simulation with keepPhylog= FALSE")

-##     ## requireNamespace("igraph")

-##     df <- x$other$PhylogDF

-##     if(!showEvents) {

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

-##     }

-##     g <- igraph::graph.data.frame(df)

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

-##     if(!timeEvent) {

-##         plot(g, layout = l0)