R/OncoSimulR.R

In OncoSimulR: Forward Genetic Simulation of Cancer Progression with Epistasis

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

## 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 3 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 <http://www.gnu.org/licenses/>.


## Does it even make sense to keepPhylog with oncoSimulSample?
## Yes, but think what we store.


oncoSimulSample <- function(Nindiv,
                            fp,
                            model = "Exp",
                            numPassengers = 0,
                            mu = 1e-6,
                            muEF = NULL,
                            detectionSize = round(runif(Nindiv, 1e5, 1e8)),
                            detectionDrivers = {
                                if(inherits(fp, "fitnessEffects")) {
                                    if(length(fp$drv)) {
                                        nd <- (2: round(0.75 * length(fp$drv)))
                                    } else {
                                        nd <- 9e6 
                                    }
                                } else {
                                    nd <- (2 : round(0.75 * max(fp)))
                                }
                                if(length(nd) == 1) ## for sample
                                       nd <- c(nd, nd)
                                sample(nd, Nindiv,
                                       replace = TRUE)
                            },
                            detectionProb = NA,
                            sampleEvery = ifelse(model %in% c("Bozic", "Exp"), 1,
                                                 0.025),
                            initSize = 500,
                            s = 0.1,
                            sh = -1,
                            K = initSize/(exp(1) - 1),
                            minDetectDrvCloneSz = "auto",
                            extraTime = 0,
                            finalTime = 0.25 * 25 * 365,
                            onlyCancer = TRUE,
                            keepPhylog = FALSE,
                            mutationPropGrowth = ifelse(model == "Bozic",
                                                        FALSE, TRUE),
                            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,
                            AND_DrvProbExit = FALSE,
                            fixation = NULL,
                            verbosity  = 1,
                            showProgress = FALSE,
                            seed = "auto"){
    ## No longer using mclapply, because of the way we use the limit on
    ## the number of tries.
    
    ## leaving detectionSize and detectionDrivers as they are, produces
    ## the equivalente of uniform sampling. For last, fix a single number

    ## detectionDrivers when there are none: had we left it at 0, then
    ## when there are no drivers we would stop at the first sampling
    ## period.
    
    if(Nindiv < 1)
        stop("Nindiv must be >= 1")
    if(keepPhylog)
        warning(paste("oncoSimulSample does not return the phylogeny",
                      "for now, so there is little point in storing it."))

    if(max.num.tries.total < Nindiv)
        stop(paste("You have requested something impossible: ",
                   "max.num.tries.total < Nindiv"))
    attemptsLeft <- max.num.tries.total
    attemptsUsed <- 0
    numToRun <- Nindiv
    pop <- vector(mode = "list", length = Nindiv)
    indiv <- 1
    ## FIXME! really pass params such as extraTime and minDDr as vectors,
    ## or give a warning
    params <- data.frame(seq.int(Nindiv),
                         extraTime = extraTime,
                         minDetectDrvCloneSz = minDetectDrvCloneSz,
                         detectionSize = detectionSize,
                         detectionDrivers = detectionDrivers,
                         stringsAsFactors = TRUE)[, -1, drop = FALSE]

    ## FIXME: we are not triggering an error, just a message. This is on
    ## purpose, since some of these conditions DO provide useful
    ## output. Do we want these to be errors?
    f.out.attempts <- function() {
        message("Run out of attempts")
        return(list(
            popSummary = NA,
            popSample = NA,
            HittedMaxTries = TRUE,
            HittedWallTime = FALSE,
            UnrecoverExcept = FALSE
        ))    
    }    

    f.out.time <- function() {
        message("Run out of time")
        return(list(
            popSummary = NA,
            popSample = NA,
            HittedMaxTries = FALSE,
            HittedWallTime = TRUE,
            UnrecoverExcept = FALSE
        ))    
    }    

    f.out.attempts.cpp <- function() {
        message("Run out of attempts (in C++)")
        return(list(
            popSummary = NA,
            popSample = NA,
            HittedMaxTries = TRUE,
            HittedWallTime = FALSE,
            UnrecoverExcept = FALSE
        ))    
    }    

    f.out.time.cpp <- function() {
        message("Run out of time (in C++)")
        return(list(
            popSummary = NA,
            popSample = NA,
            HittedMaxTries = FALSE,
            HittedWallTime = TRUE,
            UnrecoverExcept = FALSE
        ))    
    }    

    f.out.unrecover.except <- function(x) {
        message("Unrecoverable exception (in C++)")
        return(list(
            popSummary = NA,
            popSample = NA,
            HittedMaxTries = NA,
            HittedWallTime = NA,
            UnrecoverExcept = TRUE,
            ExceptionMessage = x$other$ExceptionMessage
        ))    
    }    
   
    
    startTime <- Sys.time()
    while(TRUE) {
        
        possibleAttempts <- attemptsLeft - (numToRun - 1)
        ## I think I do not want a try here.
        tmp <-  oncoSimulIndiv(fp = fp,
                               model = model,
                               numPassengers = numPassengers,
                               mu = mu,
                               muEF = muEF,
                               detectionSize = params[indiv, "detectionSize"],
                               detectionDrivers = params[indiv, "detectionDrivers"],
                               sampleEvery = sampleEvery,
                               initSize = initSize,
                               s = s,
                               sh = sh,
                               K = K,
                               minDetectDrvCloneSz = params[indiv, "minDetectDrvCloneSz"],
                               extraTime = params[indiv, "extraTime"],
                               finalTime = finalTime,
                               max.memory = max.memory,
                               max.wall.time = max.wall.time.total,
                               max.num.tries = possibleAttempts,
                               verbosity = verbosity,
                               keepEvery = -9,
                               onlyCancer = onlyCancer,
                               errorHitWallTime = TRUE,
                               errorHitMaxTries = TRUE,
                               seed = seed,
                               initMutant = initMutant,
                               keepPhylog = keepPhylog,
                               mutationPropGrowth = mutationPropGrowth,
                               detectionProb = detectionProb,
                               AND_DrvProbExit = AND_DrvProbExit,
                               fixation = fixation)        
        if(tmp$other$UnrecoverExcept) {
            return(f.out.unrecover.except(tmp))
        }
        
        pop[[indiv]] <- tmp
        numToRun <- (numToRun - 1)
        attemptsUsed <- attemptsUsed + tmp$other$attemptsUsed
        attemptsLeft <- (max.num.tries.total - attemptsUsed)
        if(showProgress) {
            cat("......  Done individual ", indiv,
                ". Used ", attemptsUsed, "attempts.",
                ". Running for ", as.double(difftime(Sys.time(),
                                                    startTime, units = "secs")),
                " secs.\n"
                )
        }
        indiv <- indiv + 1
        
        ## We need to check in exactly this order. Attempts left only
        ## matters if no remaining individuals to run. But C++ might bail
        ## out in exactly the last individual

        if(  
            (exists("HittedMaxTries", where = tmp) &&
                 tmp[["HittedMaxTries"]])  ) {
            ## in C++ code
            return(f.out.attempts.cpp())
        } else if(  
            (exists("HittedWallTime", where = tmp) &&
                 tmp[["HittedWallTime"]])  ) {
            ## in C++ code
            return(f.out.time.cpp())
        } else if( indiv > Nindiv ) {
            if(verbosity > 0)
                message(paste0("Successfully sampled ", Nindiv, " individuals"))
            class(pop) <- "oncosimulpop"
            if(inherits(fp, "fitnessEffects")) {
                geneNames <- names(getNamesID(fp))
            } else {
                geneNames <- NULL
            }
            return(list(
                popSummary = summary(pop),
                popSample = samplePop(pop, timeSample = "last",
                                      typeSample = typeSample,
                                      thresholdWhole = thresholdWhole,
                                      geneNames = geneNames),
                attemptsUsed = attemptsUsed,
                probCancer = Nindiv/attemptsUsed,
                HittedMaxTries = FALSE,
                HittedWallTime = FALSE,
                UnrecoverExcept = FALSE
            ))
        } else if( attemptsLeft <= 0 ) {
              ## it is very unlikely this will ever happen. 
            return(f.out.attempts())
        } else  if( as.double(difftime(Sys.time(), startTime, units = "secs"))
                   > max.wall.time.total ) {
            return(f.out.time())
        } 
    }
}

add_noise <- function(x, properr) {
    if(properr <= 0) {
        return(x)
    }
    else {
        if(properr > 1)
            stop("Proportion with error cannot be > 1")
        nn <- prod(dim(x))
        flipped <- sample(nn, round(nn * properr))
        x[flipped] <- as.integer(!x[flipped])
        return(x)
    }
}

samplePop <- function(x, timeSample = "last",
                      typeSample = "whole",
                      thresholdWhole = 0.5,
                      geneNames = NULL,
                      popSizeSample = NULL,
                      propError = 0) {
    ## timeSample <- match.arg(timeSample)
    gN <- geneNames
    
    if(!is.null(popSizeSample) && (length(popSizeSample) > 1) &&
       (length(popSizeSample) != length(x))) {
        message("length popSizeSample != number of subjects")
        popSizeSample <- rep(popSizeSample, length.out = length(x))
        }
    ## A hack to prevent Map from crashing with a NULL
    if(is.null(popSizeSample)) popSizeSample <- -99
    if(inherits(x, "oncosimulpop")) {
        z <- do.call(rbind,
                     Map(get.mut.vector,
                         x = x,
                         timeSample = timeSample,
                         typeSample = typeSample,
                         thresholdWhole = thresholdWhole,
                         popSizeSample = popSizeSample
                         ))
        ## We need to check if the object is coming from v.2., to avoid
        ## having to force passing a vector of names
        if(is.null(gN) && (!is.null(x[[1]]$geneNames)))
            gN <- x[[1]]$geneNames
    } else {
        z <- get.mut.vector(x,
                            timeSample = timeSample,
                            typeSample = typeSample,
                            thresholdWhole = thresholdWhole,
                            popSizeSample = popSizeSample)
        dim(z) <- c(1, length(z))
        if(is.null(gN) && (!is.null(x$geneNames)))
            gN <- x$geneNames
    }
    message("\n Subjects by Genes matrix of ",
        nrow(z), " subjects and ",
            ncol(z), " genes.\n")

    if(propError > 0) {
        z <- add_noise(z, propError)
    }
    if(!is.null(gN)) {
        colnames(z) <- gN
    } else {
        colnames(z) <- seq_len(ncol(z))
    }
    return(z)
}



## single_largest_last_pop <- function(x) {
##     last <- nrow(x$pops.by.time)
##     largest <- which.max(x$pops.by.time[last, , drop = FALSE]) - 1
##     genot <- x[["GenotypesLabels"]][largest]
##     strsplit(genot, split = ", ")[[1]]
## }


## ## like samplePop(x, timeSample = "last") but return the single most abundant genotype
## Just a prototype. Not well tested.
## largest_last_pop <- function(x) {
##     y <- lapply(x, single_largest_last_pop)
##     allg <- sort(unique(unlist(y)))
##     m <- t(vapply(y, function(z) {as.integer(allg %in%  z) },
##                   integer(length(allg)) ))
##     colnames(m) <- allg
##     return(m)
## }



oncoSimulPop <- function(Nindiv,
                         fp,
                         model = "Exp",
                         numPassengers = 0,
                         mu = 1e-6,
                         muEF = NULL,
                         detectionSize = 1e8,
                         detectionDrivers = 4,
                         detectionProb = NA,
                         sampleEvery = ifelse(model %in% c("Bozic", "Exp"), 1,
                                              0.025),
                         initSize = 500,
                         s = 0.1,
                         sh = -1,
                         K = initSize/(exp(1) - 1),
                         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,
                         mutationPropGrowth = ifelse(model == "Bozic",
                                                     FALSE, TRUE),
                         max.memory = 2000,
                         max.wall.time = 200,
                         max.num.tries = 500,
                         errorHitWallTime = TRUE,
                         errorHitMaxTries = TRUE,
                         initMutant = NULL,
                         AND_DrvProbExit = FALSE,
                         fixation = NULL,
                         verbosity  = 0,
                         mc.cores = detectCores(),
                         seed = "auto") {

    if(Nindiv < 1)
        stop("Nindiv must be >= 1")
    
    if(.Platform$OS.type == "windows") {
        if(mc.cores != 1)
            message("You are running Windows. Setting mc.cores = 1")
        mc.cores <- 1
    }
    pop <- mclapply(seq.int(Nindiv),
                    function(x)
                    oncoSimulIndiv(
                        fp = fp,
                        model = model,
                        numPassengers = numPassengers,
                        mu = mu,
                        muEF = muEF,
                        detectionSize = detectionSize,
                        detectionDrivers = detectionDrivers,
                        sampleEvery = sampleEvery,
                        initSize = initSize,
                        s = s,
                        sh = sh,
                        K = K,
                        keepEvery = keepEvery,
                        minDetectDrvCloneSz = minDetectDrvCloneSz,
                        extraTime = extraTime,
                        finalTime = finalTime,
                        onlyCancer = onlyCancer,
                        max.memory = max.memory,
                        max.wall.time = max.wall.time,
                        max.num.tries = max.num.tries,
                        errorHitWallTime = errorHitWallTime,
                        errorHitMaxTries = errorHitMaxTries,
                        verbosity = verbosity,
                        seed = seed, keepPhylog = keepPhylog,
                        initMutant = initMutant,
                        mutationPropGrowth = mutationPropGrowth,
                        detectionProb = detectionProb,
                        AND_DrvProbExit = AND_DrvProbExit,
                        fixation = fixation),
                    mc.cores = mc.cores)
    ## mc.allow.recursive = FALSE ## FIXME: remove?
                    ## done for covr issue
                    ## https://github.com/r-lib/covr/issues/335#issuecomment-424116766
                    
    class(pop) <- "oncosimulpop"
    attributes(pop)$call <- match.call()
    return(pop)
}

## where is the default K coming from? Here:
## log( (K+N)/K  ) = 1; k + n = k * exp(1); k(exp - 1) = n; k = n/(exp - 1)

oncoSimulIndiv <- function(fp,
                           model = "Exp",
                           numPassengers = 0,
                           mu = 1e-6,
                           muEF = NULL,
                           detectionSize = 1e8,
                           detectionDrivers = 4,
                           detectionProb = NA,
                           sampleEvery = ifelse(model %in% c("Bozic", "Exp"), 1,
                                                0.025),
                           initSize = 500,
                           s = 0.1,
                           sh = -1,
                           K = initSize/(exp(1) - 1),
                           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,
                           mutationPropGrowth = ifelse(model == "Bozic",
                                                       FALSE, TRUE),
                           max.memory = 2000,
                           max.wall.time = 200,
                           max.num.tries = 500,
                           errorHitWallTime = TRUE,
                           errorHitMaxTries = TRUE,
                           verbosity = 0,
                           initMutant = NULL,
                           AND_DrvProbExit = FALSE,
                           fixation = NULL,
                           seed = NULL
                           ) {
    call <- match.call()
    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 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",
                          "Exp" = "exp",
                          "McFarlandLog" = "mcfarlandlog",
                          "McFL" = "mcfarlandlog",
                          stop("No valid value for model")
                          )
    if(initSize < 1)
        stop("initSize < 1")
    
    if( (K < 1) && (model %in% c("McFL", "McFarlandLog") )) {
        stop("Using McFarland's model: K cannot be < 1")
    }       ##  if ( !(model %in% c("McFL", "McFarlandLog") )) {
            ## K <- 1 ## K is ONLY used for McFarland; set it to 1, to avoid
            ##        ## C++ blowing.

    if(typeFitness == "exp") {
        death <- 1
        ## mutationPropGrowth <- 1
    } else {
        death <- -99
        ## mutationPropGrowth <- 0
    }
    birth <- -99

    if( (typeFitness == "mcfarlandlog") &&
       (sampleEvery > 0.05)) {
        warning("With the McFarland model you often want smaller sampleEvery")
    }

    if(minDetectDrvCloneSz == "auto") {
        if(model %in% c("Bozic", "Exp") )
            minDetectDrvCloneSz <- 0
        else if (model %in% c("McFL", "McFarlandLog"))
            minDetectDrvCloneSz <- initSize
        ## minDetectDrvCloneSz <- eFinalMf(initSize, s, detectionDrivers)
        else
            stop("Unknown model")
    }

    if( (length(mu) > 1) && !inherits(fp, "fitnessEffects"))
        stop("Per-gene mutation rates cannot be used with the old poset format")

    if(any(mu < 0)) {
        stop("(at least one) mutation rate (mu) is negative")
    }
    ## We do not test for equality to 0. That might be a weird but
    ## legitimate case?

    ## No user-visible magic numbers
    ## if(is.null(keepEvery))
    ##     keepEvery <- -9
    if(is_null_na(keepEvery)) keepEvery <- -9

    
    if( (keepEvery > 0) & (keepEvery < sampleEvery)) {
        keepEvery <- sampleEvery
        warning("setting keepEvery <- sampleEvery")
    }

    if( (typeFitness == "bozic1") && (mutationPropGrowth) )
        warning("Using fitness Bozic (bozic1) with mutationPropGrowth = TRUE;",
                "this will have no effect.")

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

    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(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",
                          " when using fitnessEffects objects. ",
                          " The fitnessEffects objects are much more ",
                          "flexible and you can use, for example,",
                          "the noIntGenes component for passengers."))
        if(is.null(seed)) {## Passing a null creates a random seed
            ## We use a double, to be able to pass in range > 2^16.
            ## Do not use 0, as that is our way of signaling to C++ to
            ## generate the seed.
            seed <- round(runif(1, min = 1, max = 2^32))
            if(verbosity >= 2)
                cat(paste("\n Using ", seed, " as seed for C++ generator\n"))
        } else if(seed == "auto") {
            seed <- 0.0
            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 vector.")
            if(!(any(unlist(lapply(fixation, class)) == "character")))
                stop("At least one 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 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,  
                                        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 = initMutant, 
                                        max.wall.time = max.wall.time,
                                        max.num.tries = max.num.tries,
                                        keepEvery = keepEvery,  
                                        ## alpha = 0.0015,  
                                        K = K, 
                                        minDetectDrvCloneSz = minDetectDrvCloneSz,
                                        extraTime = extraTime,
                                        detectionDrivers = detectionDrivers,
                                        onlyCancer = onlyCancer,
                                        errorHitWallTime = errorHitWallTime,
                                        errorHitMaxTries = errorHitMaxTries,
                                        keepPhylog = keepPhylog,
                                        MMUEF = muEF,
                                        detectionProb = detectionProb,
                                        AND_DrvProbExit = AND_DrvProbExit,
                                        fixation = fixation),
                  silent = !verbosity)
        objClass <- c("oncosimul", "oncosimul2")
    }
    if(inherits(op, "try-error")) {
        ##         if(length(grep("BAIL OUT NOW", op)))
        stop(paste("Unrecoverable error:", op ))
    }
    if(verbosity >= 2) {
        cat("\n ... finished this run:")
        cat("\n       Total Pop Size = ", op$TotalPopSize)
        cat("\n       Drivers Last = ", op$MaxDriversLast)
        cat("\n       Final Time = ", op$FinalTime, "\n")
    }

    class(op) <- objClass
    attributes(op)$call <- call
    return(op)
}

summary.oncosimul <- function(object, ...) {
    ## This should be present even in HittedWallTime and HittedMaxTries
    ## 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
        return(NA)
    } else if( !pbp & (object$HittedWallTime || object$HittedMaxTries) ) {
        return(NA)
    } else if ( !pbp & !(object$HittedWallTime || object$HittedMaxTries) ) {
        stop("Eh? No pops.by.time but did not hit wall time or max tries? BUG!")
    } else {
        tmp <- object[c("NumClones", "TotalPopSize", "LargestClone",
                        "MaxNumDrivers", "MaxDriversLast",
                        "NumDriversLargestPop", "TotalPresentDrivers",
                        "FinalTime", "NumIter", "HittedWallTime",
                        "HittedMaxTries")]
 
        tmp$errorMF <- object$other$errorMF
        tmp$minDMratio <- object$other$minDMratio
        tmp$minBMratio <- object$other$minBMratio
        if( (tmp$errorMF == -99)) tmp$errorMF <- NA
        if( (tmp$minDMratio == -99)) tmp$minDMratio <- NA
        if( (tmp$minBMratio == -99)) tmp$minBMratio <- NA
        tmp$OccurringDrivers <- object$OccurringDrivers
        return(as.data.frame(tmp, stringsAsFactors = FALSE))
    }
}




print.oncosimul <- function(x, ...) {
    cat("\nIndividual OncoSimul trajectory with call:\n ")
    print(attributes(x)$call)
    cat("\n")
    print(summary(x))

    if(inherits(x, "oncosimul2")) {
        ## we know small object
        cat("\n")
        cat("Final population composition:\n")
        df <- data.frame(Genotype = x$GenotypesLabels,
                         N = x$pops.by.time[nrow(x$pops.by.time), -1],
                         stringsAsFactors = TRUE)
        print(df)
    }
}

## I want this to return things that are storable
## summary.oncosimulpop <- function(object, ...) {
##     as.data.frame(rbindlist(lapply(object, summary)))
## }

summary.oncosimulpop <- function(object, ...) {
    ## some simulations could have failed seriously returning a NULL
    ## FIXME: how, why? Out of memory?

    rm1 <- which(unlist(lapply(object, is.null)))
    if(length(rm1) > 0) {
        if(length(rm1) < length(object)) {
            warning("Some simulations returned NULL. They will be removed",
                    " from the summary. The NULL runs are ",
                    paste(rm1, collapse = ", "),
                    ".")
        } else {
            warning("All simulations returned NULL.")
            return(NA)
        }
    }

    ## I do not want to rm above for two reasons:
    ##   - avoid re-asigning to the object
    ##   - changing the numbering of the indices below if NULLs
    ## So I need something more involved
    ## Figure out exactly what the summary of a NULL is
    sumnull <- summary(NULL)
    
    tmp <- lapply(object, summary)

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

    rm <- which(unlist(lapply(tmp,
                              function(x) ((length(x) == 1) &&
                                           (is.na(x) || is.null(x))) ||
                                          (identical(x, sumnull)))))
    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), stringsAsFactors = TRUE)
}



print.oncosimulpop <- function(x, ...) {
    cat("\nPopulation of OncoSimul trajectories of",
        length(x), "individuals. Call :\n")
    print(attributes(x)$call)
    cat("\n")
    print(summary(x))
}


plot.oncosimulpop <- function(x, ask = TRUE,
                              show = "drivers", 
                              type = ifelse(show == "genotypes",
                                            "stacked", "line"),
                              col = "auto",
                              log = ifelse(type == "line", "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,
                              ylim = NULL,
                              xlim = NULL,
                              thinData = FALSE,
                              thinData.keep = 0.1,
                              thinData.min = 2,
                              plotDiversity = FALSE,
                              order.method = "as.is",
                              stream.center = TRUE,
                              stream.frac.rand = 0.01,
                              stream.spar = 0.2,
                              border = NULL,
                              lwdStackedStream = 1,
                              srange = c(0.4, 1),
                              vrange = c(0.8, 1),
                              breakSortColors = "oe",
                              legend.ncols = "auto",
                              ...
                              ) {
    op <- par(ask = ask)
    on.exit(par(op))
    null <- lapply(x, function(z)
        plot.oncosimul(z,
                       show = show,
                       type = type,
                       col = col,
                       log = log,
                       ltyClone = ltyClone,
                       lwdClone = lwdClone,
                       ltyDrivers = ltyDrivers,
                       lwdDrivers = lwdDrivers,
                       xlab = xlab,
                       ylab = ylab,
                       plotClones = plotClones,
                       plotDrivers = plotDrivers,
                       addtot = addtot,
                       addtotlwd = addtotlwd,
                       ylim = ylim,
                       xlim = xlim,
                       thinData = thinData,
                       thinData.keep = thinData.keep,
                       thinData.min = thinData.min,
                       plotDiversity = plotDiversity,
                       order.method = order.method,
                       stream.center = stream.center,
                       stream.frac.rand = stream.frac.rand,
                       stream.spar = stream.spar,
                       border = border,
                       lwdStackedStream = lwdStackedStream,
                       srange = srange,
                       vrange = vrange,
                       breakSortColors = breakSortColors,
                       legend.ncols = legend.ncols,
                       ...))
}


## 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", 
                           type = ifelse(show == "genotypes",
                                         "stacked", "line"),
                           col = "auto",
                           log = ifelse(type == "line", "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,
                           ylim = NULL,
                           xlim = NULL,
                           thinData = FALSE,
                           thinData.keep = 0.1,
                           thinData.min = 2,
                           plotDiversity = FALSE,
                           order.method = "as.is",
                           stream.center = TRUE,
                           stream.frac.rand = 0.01,
                           stream.spar = 0.2,
                           border = NULL,
                           lwdStackedStream = 1,
                           srange = c(0.4, 1),
                           vrange = c(0.8, 1),
                           breakSortColors = "oe",
                           legend.ncols = "auto",
                           ...
                           ) {


    if(!(type %in% c("stacked", "stream", "line")))
        stop("Type of plot unknown: it must be one of",
             "stacked, stream or line")

    if(!(show %in% c("genotypes", "drivers")))
        stop("show must be one of ",
             "genotypes or drivers")

    if(!(breakSortColors %in% c("oe", "distave", "random")))
        stop("breakSortColors must be one of ",
             "oe, distave, or random")

    

    colauto <- FALSE
    if((length(col) ==1) && (col == "auto") &&
       (type == "line") && (show == "drivers"))
        col <- c(8, "orange", 6:1)
    if((length(col) ==1) && (col == "auto") &&
       (show == "genotypes")) {
        ## For categorical data, I find Dark2, Paired, or Set1 to work best.
        col <- colorRampPalette(brewer.pal(8, "Dark2"))(ncol(x$pops.by.time) - 1)
        colauto <- TRUE
    }
    
    if(show == "genotypes") {
        plotDrivers <- FALSE
        plotClones <- TRUE
    }
    
    if(thinData)
        x <- thin.pop.data(x, keep = thinData.keep, min.keep = thinData.min)

    if(!is.null(xlim))
        x <- xlim.pop.data(x, xlim)
    
    ## For genotypes, ndr is now the genotypes.  Actually, ndr is now just
    ## a sequence 1:(ncol(y) - 1)

    ## The user will want to change the colors, like a colorRamp, etc. Or
    ## rainbow.

    ## genotypes and line, always call plotDrivers0
    if(show == "drivers") {
        if(!inherits(x, "oncosimul2"))
            ndr <- colSums(x$Genotypes[1:x$NumDrivers, , drop = FALSE])
        else {
            ndr <- colSums(x$Genotypes[x$Drivers, , drop = FALSE])
        }
    } else { ## show we are showing genotypes
        ndr <- 1:(ncol(x$pops.by.time) - 1)
    }
    
    if((type == "line") && is.null(ylim)) {
        if(log %in% c("y", "xy", "yx") )
            ylim <- c(1, max(apply(x$pops.by.time[, -1, drop = FALSE], 1, sum)))
        else
            ylim <- c(0, max(apply(x$pops.by.time[, -1, drop = FALSE], 1, sum)))
    }
    if(plotDiversity) {
        oppd <- 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)  
    }

    ## Shows its history: plotClones makes plotDrivers0 unneeded with
    ## stacked and stream plots. But now so with line plot.
    ## When showing genotypes, plotDrivers0 with line only used for
    ## showing the legend.
    if(plotClones) {
        plotClonesSt(x,
                     ndr = ndr,
                     show = show,
                     na.subs = TRUE,
                     log = log,
                     lwd = lwdClone,
                     lty = ifelse(show == "drivers", ltyClone, ltyDrivers),
                     col = col, 
                     order.method = order.method,
                     stream.center = stream.center,
                     stream.frac.rand = stream.frac.rand,
                     stream.spar = stream.spar,
                     border = border,
                     srange = srange,
                     vrange = vrange,
                     type = type,
                     breakSortColors = breakSortColors,
                     colauto = colauto,
                     legend.ncols = legend.ncols,
                     lwdStackedStream = lwdStackedStream,
                     xlab = xlab,
                     ylab = ylab,
                     ylim = ylim,
                     xlim = xlim,
                     ...)
    }

    if(plotClones && plotDrivers && (type == "line"))
        par(new = TRUE)
    
    if( plotDrivers && (type == "line") ) {
        plotDrivers0(x,
                     ndr,
                     timescale = 1,
                     trim.no.drivers = FALSE,
                     xlab = "", ylab = "",
                     lwd = lwdDrivers,
                     lty = ltyDrivers,
                     col = col, 
                     addtot = addtot,
                     addtotlwd = addtotlwd,
                     log = log, ylim = ylim,
                     xlim = xlim,
                     legend.ncols = legend.ncols,
                     ...)
    }
    if(plotDiversity) {
        par(oppd)
    }
    
}

plotClonesSt <- function(z,
                         ndr,
                         show = "drivers",
                         na.subs = TRUE,
                         log = "y",
                         lwd = 1,
                         ## type = "l",
                         lty = 1:8, col = 1:9,
                         order.method = "as.is",
                         stream.center = TRUE,
                         stream.frac.rand = 0.01,
                         stream.spar = 0.2,
                         border = NULL,
                         srange = c(0.4, 1),
                         vrange = c(0.8, 1),
                         type = "stacked",
                         breakSortColors = "oe",
                         colauto = TRUE,
                         legend.ncols = "auto",
                         lwdStackedStream = 1,
                         xlab = "Time units",
                         ylab = "Number of cells",
                         ylim = NULL,
                         xlim = NULL,
                         ...) {

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

    ## Code in stacked and stream plots relies on there being no NAs. Could
    ## change it, but it does not seem reasonable.
    ##  But my original plotting code runs faster and is simpler if 0 are
    ##  dealt as NAs (which also makes log transformations simpler).
    
    if(type %in% c("stacked", "stream") )
        na.subs <- FALSE
    
    if(na.subs){
        y[y == 0] <- NA
    }
    ## if(is.null(ndr))
    ##     stop("Should never have null ndr")
    ## if(!is.null(ndr)) {
        ## could be done above, to avoid creating
        ## more copies
    oo <- order(ndr)
    y <- y[, oo, drop = FALSE]
    ndr <- ndr[oo]
    if(show == "drivers") {
        col <- rep(col, length.out = (1 + max(ndr)))[ndr + 1]
        lty <- rep(lty, length.out = ncol(y))
    } else {
        if(length(col) < max(ndr))
            warning("Repeating colors; you might want to",
                    "pass a col vector of more elements")
        col <- rep(col, length.out = (max(ndr)))[ndr]
    }
    ## }
    if(type == "line") {
        matplot(x = z$pops.by.time[, 1],
                y = y,
                log = log, type = "l",
                col = col, lty = lty,
                lwd = lwd,
                xlab = xlab,
                ylab = ylab,
                ylim = ylim,
                xlim = xlim,
                ...)
        box()
        if(show == "genotypes") {
            if(!inherits(z, "oncosimul2")) {
                ldrv <- genotypeLabel(z)
            } else {
                ldrv <- z$GenotypesLabels
            }
            ldrv[ldrv == ""] <- "WT"
            ldrv[ldrv == " _ "] <- "WT"
            if(legend.ncols == "auto") {
                if(length(ldrv) > 6) legend.ncols <- 2
                else legend.ncols <- 1
            }
            legend(x = "topleft",
                   title = "Genotypes",
                   lty = lty,
                   col = col, lwd = lwd,
                   legend = ldrv,
                   ncol = legend.ncols)
        }
    } else {
        ymax <- colSums(y)
        if((show == "drivers") || ((show == "genotypes") && (colauto))) {
            cll <- myhsvcols(ndr, ymax, srange = srange, vrange = vrange,
                             breakSortColors = breakSortColors)
        } else {
            cll <- list(colors = col)
        }
        x <- z$pops.by.time[, 1]
        if(grepl("y", log)) {
            stop("It makes little sense to do a stacked/stream",
                 "plot after taking the log of the y data.")
        }
        if(grepl("x", log)) {
            x <- log10(x + 1)
        }
        
        if (type == "stacked") {
            plot.stacked2(x = x,
                          y = y,
                          order.method = order.method,
                          border = border,
                          lwd = lwdStackedStream,
                          col = cll$colors,
                          log = log,
                          xlab = xlab,
                          ylab = ylab,
                          ylim = ylim,
                          xlim = xlim,
                          ...) 
        } else if (type == "stream") {
            plot.stream2(x = x,
                         y = y,
                         order.method = order.method,
                         border = border,
                         lwd = lwdStackedStream,
                         col = cll$colors,
                         frac.rand = stream.frac.rand,
                         spar = stream.spar,
                         center = stream.center,
                         log = log,
                         xlab = xlab,
                         ylab = ylab,
                         ylim = ylim,
                         xlim = xlim,
                         ...)
        }
        if(show == "drivers") {
            if(legend.ncols == "auto") {
                if(length(cll$colorsLegend$Drivers) > 6) legend.ncols <- 2
                else legend.ncols <- 1
            }
            legend(x = "topleft",
                   title = "Number of drivers",
                   pch = 15,
                   ## lty = 1,
                   ## lwd = 2,
                   col = cll$colorsLegend$Color,
                   legend = cll$colorsLegend$Drivers,
                   ncol = legend.ncols)
        } else if (show == "genotypes") {
            if(!inherits(z, "oncosimul2")) {
                ldrv <- genotypeLabel(z)
            } else {
                ldrv <- z$GenotypesLabels
            }
            ldrv[ldrv == ""] <- "WT"
            ldrv[ldrv == " _ "] <- "WT"            
            if(legend.ncols == "auto") {
                if(length(ldrv) > 6) legend.ncols <- 2
                else legend.ncols <- 1
            }
            legend(x = "topleft",
                   title = "Genotypes",
                   pch = 15,
                   col = cll$colors,
                   legend = ldrv,
                   ncol = legend.ncols)
        }
    }
}





relabelLogaxis <- function(side = 2) {
    ## we do a plot but pass the already transformed data; relabel
    ## afterwards.
    po <- axis( side = side, labels = FALSE, tick = FALSE, lwd = 0)
    axis(side = side, labels = 10^po, at = po, tick = TRUE)
}




myhsvcols <- function(ndr, ymax, srange = c(0.4, 1),
                      vrange = c(0.8, 1),
                      breakSortColors = "oe") {
    ## Generate a set of colors so that:
    ##  - easy to tell when we increase number of drivers
    ##  - reasonably easy to tell in a legend
    ##  - different clones with same number of drivers have "similar" colors

    ## I use hsv color specification as this seems the most reasonable.
    
    minor <- table(ndr)
    major <- length(unique(ndr)) ## yeah same as length(minor), but least
                                 ## surprise
    
    h <- seq(from = 0, to = 1, length.out = major + 1)[-1]
    ## do not keep similar hues next to each other
    if(breakSortColors == "oe") {
        oe <- seq_along(h) %% 2
        h <- h[order(oe, h)]
    } else if(breakSortColors == "distave"){
        sl <- seq_along(h)
        h <- h[order(-abs(mean(sl) - sl))]
    } else if(breakSortColors == "random") {
        rr <- order(runif(length(h)))
        h <- h[rr]
    } 
    
    hh <- rep(h, minor)
    
    sr <- unlist(lapply(minor, function(x) 
        seq(from = srange[1], to = srange[2], length.out = x)))
    sv <- unlist(lapply(minor, function(x) 
        seq(from = vrange[1], to = vrange[2], length.out = x))
        )

    colors <- hsv(hh, sr, sv)

    ## This gives "average" or "median" color for legend
    ## colorsLegend <- aggregate(list(Color = colors), list(Drivers = ndr),
    ##                           function(x)
    ##                               as.character(x[((length(x) %/% 2) + 1 )]))

    ## Give the most abundant class color as the legend. Simpler to read
    colorsLegend <- by(data.frame(Color = colors, maxnum = ymax),
                       list(Drivers = ndr),
                       function(x) as.character(x$Color[which.max(x$maxnum)]))
    colorsLegend <- data.frame(Drivers = as.integer(row.names(colorsLegend)),
                               Color = cbind(colorsLegend)[, 1],
                               stringsAsFactors = FALSE)
    ## To show what it would look like
    ## plot(1:(sum(minor)), col = colors, pch = 16, cex = 3)
    ## legend(1, length(ndr), col = colorsLegend$Color, legend = names(minor),
    ##        pch = 16)

    return(list(colors = colors,
                colorsLegend = colorsLegend))
}

genotypeLabel <- function(x) {
    ## For oncosimul objects, that do not have the GenotypesLabels object
    apply(x$Genotypes, 2,
          function(x) paste(which(x == 1), sep = "", collapse = ", "))
}

plotDrivers0 <- function(x,
                         ndr,
                         timescale = 1,
                         trim.no.drivers = FALSE,
                         addtot = TRUE,
                         addtotlwd = 2,
                         na.subs = TRUE, log = "y", type = "l",
                         lty = 1:9, col = c(8, "orange", 6:1),
                         lwd = 2,
                         legend.ncols = "auto",
                         ...) {
    ## z <- create.drivers.by.time(x, numDrivers)
    z <- create.drivers.by.time(x, ndr)
    ## We can now never enter here because trim.no.drivers is always FALSE
    ## in call.
    ## if(trim.no.drivers && x$MaxDriversLast) {
    ##     fi <- which(apply(z[, -c(1, 2), drop = FALSE], 1,
    ##                       function(x) sum(x) > 0))[1]
    ##     z <- z[fi:nrow(z), , drop = FALSE]
    ## }
    y <- z[, 2:ncol(z), drop = FALSE]
    if(na.subs){
        y[y == 0] <- NA
    }

    ## Likewise, we can never enter here now as timescale fixed at 1. And
    ## this is silly.
    ## if(timescale != 1) {
    ##     time <- timescale * z[, 1]
    ## } else {
    ##     time <- z[, 1]
    ## }
    time <- timescale * z[, 1]
    if(nrow(y) <= 2) type <- "b"
    ## Allow for the weird case of possibly missing drivers
    col2 <- rep(col, length.out = (1 + max(ndr)))
    col2 <- col2[sort(unique(ndr)) + 1]
    matplot(x = time,
            y = y,
            type = type, log = log, lty = lty, col = col2, lwd = lwd,
            ...)
    if(addtot) {
        tot <- rowSums(y, na.rm = TRUE)
        lines(time, tot, col = "black", lty = 1, lwd = addtotlwd)
    }
    
    ## This will work even if under the weird case of a driver missing
    ldrv <- unlist(lapply(strsplit(colnames(y), "dr_", fixed = TRUE),
                          function(x) x[2]))
    if(legend.ncols == "auto") {
        if(length(ldrv) > 6) legend.ncols <- 2
        else legend.ncols <- 1
    }
    legend(x = "topleft",
           title = "Number of drivers",
           lty = lty, col = col2, lwd = lwd,
           legend = ldrv,
           ncol = legend.ncols)
    ## legend = (1:ncol(y)) - 1)
}





plotPoset <- function(x, names = NULL, addroot = FALSE,
                      box = FALSE, ...) {
    if(is.null(names)) {
        if(addroot) names <- c("Root", 1:max(x))
        else names <- 1:max(x)
    }
    plot(posetToGraph(x, names =  names, addroot = addroot,
                      type = "graphNEL",
                      strictAdjMat = FALSE), ...)
    if(box)
        box()
}

## this function seems to never be used
## plotAdjMat <- function(adjmat) {
##     plot(as(adjmat, "graphNEL"))
## }



which_N_at_T <- function(x, N = 1, t = "last") {
    if((length(t) == 1) && (t == "last"))
        T <- nrow(x$pops.by.time)
    else if(length(t) == 2) {
        if(t[1] < 0)
            warning("smallest t must be >= 0; setting to 0")
        if(t[1] > t[2])
            stop("t[1] must be <= t[2]")
        if(t[2] > max(x$pops.by.time[, 1]))
            message("t[2] > largest time; setting it to the max")
        T <- which(
            (x$pops.by.time[, 1] >= t[1]) &
                (x$pops.by.time[, 1] <= t[2]))
    }
    else
        stop("t must be either 'last' or a vector of length 2")
    z <- which(x$pops.by.time[T, -1, drop = FALSE] >= N, arr.ind = TRUE)[, 2]
    z <- unique(z) ## recall we removed first column but we index from first.
    return(z)
}

phylogClone <- function(x, N = 1, t = "last", keepEvents = TRUE) {
    if(!inherits(x, "oncosimul2"))
        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
    ## FIXME: aren't this and the next warnings redundant or aliased?
    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. ",
                "This also happens if you did not set 'keepPhylog = TRUE'.")
        return(NA)
    }
    if(!keepEvents) { ## is this just a graphical thing? or not?
        df <- df[!duplicated(df[, c(1, 2)]), ]
    }
    g <- igraph::graph.data.frame(df[, c(1, 2)])
    ## nodes <- match(tG, V(g)$name)
    nodesInP <- unique(unlist(igraph::neighborhood(g, order = 1e9,
                                                   nodes = tG,
                                                   mode = "in")))
    ## Remember that the phylog info can contain clones that are
    ## not in pops.by.time, as they go extinct between creation
    ## and sampling.
    allLabels <- unique(as.character(unlist(df[, c(1, 2)])))
    nodesRm <- setdiff(allLabels, V(g)$name[nodesInP])
    g <- igraph::delete.vertices(g, nodesRm)
    tmp <- list(graph = g, df = df)
    class(tmp) <- c(class(tmp), "phylogClone")
    return(tmp)
    ## trivial to return an adjacency matrix if needed. The keepEvents = FALSE
}



plotClonePhylog <- function(x, N = 1, t = "last",
                            timeEvents = FALSE,
                            keepEvents = FALSE,
                            fixOverlap = TRUE,
                            returnGraph = FALSE, ...) {
    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. ",
             "This error can also appear if your simulation exited ",
             "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)
    } else {
        l1 <- l0
        indexAppear <- match(V(pc$g)$name, as.character(pc$df[, 2]))
        firstAppear <- pc$df$time[indexAppear]
        firstAppear[1] <- 0
        l1[, 2] <- (max(firstAppear) - firstAppear)
        if(fixOverlap) {
            dx <- which(duplicated(l1[, 1]))
            if(length(dx)) {
                ra <- range(l1[, 1])
                l1[dx, 1] <- runif(length(dx), ra[1], ra[2])
            }
        }
        plot(pc$g, layout = l1)         
    }
    if(returnGraph)
        return(pc$g)
}








############# The rest are internal functions

closest_time <- function(x, size) {
    ## Find the first time when a given size is reached
    ## But these are not times, but the position in pops.by.time. Bad naming.
    sizes <- rowSums(x$pops.by.time[, -1, drop = FALSE])
    candidates <- which(sizes >= size)
    if(length(candidates) == 0) {
        warning(paste("Pop size never >= requested size.",
                      "Thus, there is nothing to sample. You will get NAs"))
        return(-99)
    } else {
        return(candidates[1])
    }
}


get.the.time.for.sample <- function(tmp, timeSample, popSizeSample) {
    if( !is.null(popSizeSample) && (popSizeSample >= 0) )  {
        ## should be "closest_index" and "the.index"
        the.time <- closest_time(tmp, popSizeSample)
    } else 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"))
            the.time <- -99
        } else {
            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 with drivers. ",
                            "Thus, no uniform sampling possible.",
                            "You will get NAs"))
              the.time <- -99
              ## return(rep(NA, nrow(tmp$Genotypes)))
          } else if (length(candidate.time) == 1) {
                message("Only one sampled time period with mutants.")
                the.time <- candidate.time
            } else {
                  the.time <- sample(candidate.time, 1)
              }
    } else {
            stop("Unknown timeSample option")
        }
    return(the.time)
}


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 if (typeSample %in% c("singleCell-noWT", "single-noWT",
                                 "singleCell-nowt", "single-nowt")) {
        genots <- x$Genotypes
        whichwt <- which(x$GenotypesLabels == "")
        if(length(whichwt)) {
            genots <- genots[, -whichwt, drop = FALSE]
            pop <- pop[-whichwt]
        }
        if(all(pop == 0)) {
            warning("No non-WT clone with required popSize or at required time")
            return(rep(NA, nrow(x$Genotypes)))
        } else {
            return(genots[, sample(seq_along(pop), 1, prob = pop)])
        }
    } else {
        stop("Unknown typeSample option")
    }
}


## 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
    
    if(!inherits(x, "oncosimul2")) {
        ndr <- colSums(x$Genotypes[1:x$NumDrivers, , drop = FALSE])
        genotLabels <- genotypeLabel(x)
    } else {
        ndr <- colSums(x$Genotypes[x$Drivers, , drop = FALSE])
        genotLabels <- x$GenotypesLabels
    }
    genotLabels[genotLabels == ""] <- "WT"
    genotLabels[genotLabels == " _ "] <- "WT"
    y <- x$pops.by.time[, 2:ncol(x$pops.by.time), drop = FALSE]
    y[y == 0] <- NA
    
    oo <- order(ndr)
    y <- y[, oo, drop = FALSE]
    ndr <- ndr[oo]

    nc <- ncol(y)
    nr <- nrow(y)
    y <- as.vector(y)
    return(data.frame(Time = rep(x$pops.by.time[, 1], nc),
                      Y = y,
                      Drivers = factor(rep(ndr, rep(nr, nc))),
                      Genotype = rep(genotLabels, rep(nr, nc)),
                      stringsAsFactors = TRUE))
}



## 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
    if(tmp$NumClones >= 1) {
        if(tmp$NumClones == 1) {
            if(ncol(tmp$pops.by.time) != 2)
                stop("This is impossible!")
            drivers.by.time <- tmp$pops.by.time
        } else {
            if(length(unique(CountNumDrivers )) > 1) {
                drivers.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,
                                                        CountNumDrivers,
                                                        sum)))) 
            } else {
                drivers.by.time <- cbind(tmp$pops.by.time[, c(1),
                                                          drop = FALSE] ,
                                         rowSums(tmp$pops.by.time[, -c(1),
                                                                  drop =FALSE]))
            }
        }
        colnames(drivers.by.time) <- c("Time",
                                       paste("dr_",
                                             colnames(drivers.by.time)[-c(1)],
                                             sep = ""))
    } else {
        drivers.by.time <- NULL
    }
    return(drivers.by.time)
} 



## For plotting, this helps decrease huge file sizes, while still showing
## the start of each clone, if it was originally recorded.

thin.pop.data <- function(x, keep = 0.1, min.keep = 3) {
    norig <- nrow(x$pops.by.time)
    keep1 <- round(seq.int(from = 1, to = norig,
                           length.out = round(norig * keep)))
    keep2 <- apply(x$pops.by.time[, -1, drop = FALSE],
                   1, function(x) any((x[x > 0] < min.keep)))
    keep <- sort(union(keep1, keep2))
    x$pops.by.time <- x$pops.by.time[keep, , drop = FALSE]
    return(x)
}

xlim.pop.data <- function(x, xlim) {
    x$pops.by.time <- x$pops.by.time[
    (x$pops.by.time[, 1] >= xlim[1]) &
    (x$pops.by.time[, 1] <= xlim[2]),   ]
    return(x)
}


shannonI <- function(x) {
    sx <- sum(x)
    p <- x/sx
    p <- p[p > 0]
    return(-sum(p * log(p)))
}

plotShannon <- 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, type = "l", xlab = "", ylab = "H", axes = FALSE)
    box()
    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
##     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)
##     } else {
##         l1 <- l0
##         indexAppear <- match(V(g)$name, as.character(df[, 2]))
##         firstAppear <- df$time[indexAppear]
##         firstAppear[1] <- 0
##         l1[, 2] <- (max(firstAppear) - firstAppear)
##         if(fixOverlap) {
##             dx <- which(duplicated(l1[, 1]))
##             if(length(dx)) {
##                 ra <- range(l1[, 1])
##                 l1[dx, 1] <- runif(length(dx), ra[1], ra[2])
##             }
##         }
##         plot(g, layout = l1)         
##     }
## }