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

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progresion with Epistasis 

-Version: 2.1.1

-Date: 2016-03-07

+Version: 2.1.2

+Date: 2016-03-27

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

 		     email = "rdiaz02@gmail.com"),

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

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

                                     if(length(fp$drv)) {

                                         nd <- (2: round(0.75 * length(fp$drv)))

                                     } else {

-                                        nd <- 0

+                                        nd <- 9e6 

                                     }

                                 } else {

                                     nd <- (2 : round(0.75 * max(fp)))

@@ -69,6 +69,12 @@ oncoSimulSample <- function(Nindiv,

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

@@ -277,7 +283,7 @@ samplePop <- function(x, timeSample = "last", typeSample = "whole",

 oncoSimulPop <- function(Nindiv,

                          fp,

                          model = "Exp",

-                         numPassengers = 30,

+                         numPassengers = 0,

                          mu = 1e-6,

                          detectionSize = 1e8,

                          detectionDrivers = 4,

@@ -308,6 +314,9 @@ oncoSimulPop <- function(Nindiv,

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

@@ -354,7 +363,7 @@ oncoSimulPop <- function(Nindiv,

 oncoSimulIndiv <- function(fp,

                            model = "Exp",

-                           numPassengers = 30,

+                           numPassengers = 0,

                            mu = 1e-6,

                            detectionSize = 1e8,

                            detectionDrivers = 4,

@@ -393,6 +402,14 @@ oncoSimulIndiv <- function(fp,

                           "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

@@ -417,7 +434,7 @@ oncoSimulIndiv <- function(fp,

             stop("Unknown model")

     }

-    if(mu < 0) {

+    if(any(mu < 0)) {

         stop("mutation rate (mu) is negative")

     }

     ## We do not test for equality to 0. That might be a weird but

@@ -502,6 +519,12 @@ oncoSimulIndiv <- function(fp,

                   silent = !verbosity)

         objClass <- "oncosimul"

     } else {

+        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

@@ -1541,41 +1564,41 @@ oncoSimul.internal <- function(poset, ## restrict.table,

     rtC <- convertRestrictTable(restrict.table)

     return(c(

-        BNB_Algo5(rtC,

-        numDrivers,

-        numGenes,

-        typeCBN,

-        birth, 

-        s, 

-        death,

-        mu,

-        initSize,

-        sampleEvery,

-        detectionSize,

-        finalTime,

-        initSize_species,

-        initSize_iter,

-        seed,

-        verbosity,

-        speciesFS,

-        ratioForce,

-        typeFitness,

-        max.memory,

-        as.integer(mutationPropGrowth),

-        initMutant,

-        max.wall.time,

-        keepEvery,

-        alpha,

-        sh,

-        K,

-        # endTimeEvery,

-        detectionDrivers,

-        onlyCancer,

-        errorHitWallTime,

-        max.num.tries,

-        errorHitMaxTries,

-        minDetectDrvCloneSz,

-        extraTime

+        BNB_Algo5(restrictTable = rtC,

+        numDrivers = numDrivers,

+        numGenes = numGenes,

+        typeCBN_ = typeCBN,

+        birthRate = birth,

+        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,

+        alpha = alpha,

+        sh = sh,

+        K = K,

+        # end = TimeEvery# , 

+        detectionDrivers = detectionDrivers,

+        onlyCancer = onlyCancer,

+        errorHitWallTime = errorHitWallTime,

+        maxNumTries = max.num.tries,

+        errorHitMaxTries = errorHitMaxTries,

+        minDetectDrvCloneSz = minDetectDrvCloneSz,

+        extraTime = extraTime

     ),

              NumDrivers = numDrivers

              ))

@@ -2,19 +2,18 @@

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

 nr_BNB_Algo5 <- function(rFE, mu, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, alpha, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog) {

-    .Call('C_OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, alpha, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog)

+    .Call('OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, alpha, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog)

 }

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

-    .Call('C_OncoSimulR_BNB_Algo5', PACKAGE = 'OncoSimulR', restrictTable, numDrivers, numGenes, typeCBN_, birthRate, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, alpha, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime)

+    .Call('OncoSimulR_BNB_Algo5', PACKAGE = 'OncoSimulR', restrictTable, numDrivers, numGenes, typeCBN_, birthRate, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, alpha, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime)

 }

-## No longer used (or at least not now)

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

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

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

 ## }

 evalRGenotype <- function(rG, rFE, verbose, prodNeg) {

-    .Call('C_OncoSimulR_evalRGenotype', PACKAGE = 'OncoSimulR', rG, rFE, verbose, prodNeg)

+    .Call('OncoSimulR_evalRGenotype', PACKAGE = 'OncoSimulR', rG, rFE, verbose, prodNeg)

 }

@@ -32,14 +32,33 @@

 check.gm <- function(gm) {

-    ## Yes, Root: we want no ambiguities

-    if(gm[1] != "Root")

-        stop("First value of a module table must be Root")

-    if(names(gm)[1] != "Root")

-        stop("First name of a module table must be Root")

+    ## Yes, Root: we want no ambiguities.

+

+    ## Actually, that sucks. So we do not require it, but check for

+    ## consistency.

+    

+    if(any(gm == "Root") && (gm[1] != "Root") )

+        stop("If Root is in the module table, it must be the first")

+

+    if(any(names(gm) == "Root") && (names(gm)[1] != "Root") )

+        stop("If the name Root is in the module table, it must be the first")

+

+    if( (names(gm)[1] == "Root") && (gm[1] != "Root") )

+        stop("The name Root is in the module table, but is not of Root")

+

+    if( (gm[1] == "Root") && (names(gm)[1] != "Root") )

+        stop("Root is in the module table, but with a different name")

+

+    ## if(gm[1] != "Root")

+    ##     stop("First value of a module table must be Root")

+    ## if(names(gm)[1] != "Root")

+    ##     stop("First name of a module table must be Root")

     if(length(unique(names(gm))) != length(gm))

         stop("Number of unique module names different from length of vector")

-    

+

+    if(gm[1] != "Root")

+        gm <- c("Root" = "Root", gm)

+    return(gm)

 }

 gtm2 <- function(x) {

@@ -65,7 +84,8 @@ nice.vector.eo <- function(z, sep, rm.sign = FALSE) {

 gm.to.geneModuleL <- function(gm, one.to.one) {

     ## the table will be sorted by gene name

-    check.gm(gm)

+    gm <- check.gm(gm)

+   

     ## the named vector with the mapping into the long geneModule df

     geneMod <- as.data.frame(rbindlist(lapply(gm, gtm2)))

     geneMod$Module <- names(gm)[geneMod[, 2]] ## reverse lookup table

@@ -724,6 +744,9 @@ fitnessEffectsToIgraph <- function(rT, epistasis, orderEffects) {

         df2 <- to.pairs.modules(orderEffects, ">")

     }

     df <- rbind(df0, df1, df2)

+    ## for special case of simple epi effects

+    if(nrow(df) == 0) return(NULL)

+    

     g1 <- graph.data.frame(df)

     E(g1)$color <- "black"

     E(g1)$color[E(g1)$typeDep == "SM"] <- "blue"

@@ -21,26 +21,26 @@

-#plot.stream makes a "stream plot" where each y series is plotted 

-#as stacked filled polygons on alternating sides of a baseline.

-#

-#Arguments include:

-#'x' - a vector of values

-#'y' - a matrix of data series (columns) corresponding to x

-#'order.method' = c("as.is", "max", "first") 

-#  "as.is" - plot in order of y column

-#  "max" - plot in order of when each y series reaches maximum value

-#  "first" - plot in order of when each y series first value > 0

-#'center' - if TRUE, the stacked polygons will be centered so that the middle,

-#i.e. baseline ("g0"), of the stream is approximately equal to zero. 

-#Centering is done before the addition of random wiggle to the baseline. 

-#'frac.rand' - fraction of the overall data "stream" range used to define the range of

-#random wiggle (uniform distrubution) to be added to the baseline 'g0'

-#'spar' - setting for smooth.spline function to make a smoothed version of baseline "g0"

-#'col' - fill colors for polygons corresponding to y columns (will recycle)

-#'border' - border colors for polygons corresponding to y columns (will recycle) (see ?polygon for details)

-#'lwd' - border line width for polygons corresponding to y columns (will recycle)

-#'...' - other plot arguments

+## plot.stream makes a "stream plot" where each y series is plotted 

+## as stacked filled polygons on alternating sides of a baseline.

+

+## Arguments include:

+## 'x' - a vector of values

+## 'y' - a matrix of data series (columns) corresponding to x

+## 'order.method' = c("as.is", "max", "first") 

+##  "as.is" - plot in order of y column

+##  "max" - plot in order of when each y series reaches maximum value

+##  "first" - plot in order of when each y series first value > 0

+## 'center' - if TRUE, the stacked polygons will be centered so that the middle,

+## i.e. baseline ("g0"), of the stream is approximately equal to zero. 

+## Centering is done before the addition of random wiggle to the baseline. 

+## 'frac.rand' - fraction of the overall data "stream" range used to define the range of

+## random wiggle (uniform distrubution) to be added to the baseline 'g0'

+## 'spar' - setting for smooth.spline function to make a smoothed version of baseline "g0"

+## 'col' - fill colors for polygons corresponding to y columns (will recycle)

+## 'border' - border colors for polygons corresponding to y columns (will recycle) (see ?polygon for details)

+## 'lwd' - border line width for polygons corresponding to y columns (will recycle)

+## '...' - other plot arguments

 plot.stream2 <- function(

                         x, y, 

@@ -135,20 +135,20 @@ plot.stream2 <- function(

-#plot.stacked makes a stacked plot where each y series is plotted on top

-#of the each other using filled polygons

-#

-#Arguments include:

-#'x' - a vector of values

-#'y' - a matrix of data series (columns) corresponding to x

-#'order.method' = c("as.is", "max", "first") 

-#  "as.is" - plot in order of y column

-#  "max" - plot in order of when each y series reaches maximum value

-#  "first" - plot in order of when each y series first value > 0

-#'col' - fill colors for polygons corresponding to y columns (will recycle)

-#'border' - border colors for polygons corresponding to y columns (will recycle) (see ?polygon for details)

-#'lwd' - border line width for polygons corresponding to y columns (will recycle)

-#'...' - other plot arguments

+## plot.stacked makes a stacked plot where each y series is plotted on top

+## of the each other using filled polygons

+

+## Arguments include:

+## 'x' - a vector of values

+## 'y' - a matrix of data series (columns) corresponding to x

+## 'order.method' = c("as.is", "max", "first") 

+##  "as.is" - plot in order of y column

+##  "max" - plot in order of when each y series reaches maximum value

+##  "first" - plot in order of when each y series first value > 0

+## 'col' - fill colors for polygons corresponding to y columns (will recycle)

+## 'border' - border colors for polygons corresponding to y columns (will recycle) (see ?polygon for details)

+## 'lwd' - border line width for polygons corresponding to y columns (will recycle)

+## '...' - other plot arguments

 plot.stacked2 <- function(

                          x, y, 

@@ -1,3 +1,10 @@

+Changes in version 2.1.2 (2016-03-27):

+	- Arguments to BNB_Algo5 explicit.

+	- Example of modules and no epistasis.

+	- Remove requirement of Root in geneToModule.

+	- More tests (and reorganized them)

+	- Miscell. improvements in documentation and vignette.

+

 Changes in version 2.1.1 (2016-03-07):

 	- Added mutationPropGrowth as argument.

 	- Stacked area and stream plots (code from Marc Taylor).

@@ -70,16 +70,23 @@ allFitnessEffects(rT = NULL, epistasis = NULL, orderEffects = NULL,

 }

 \item{noIntGenes}{

   A numeric vector (optionally named) with the fitness coefficients of genes

-  (only genes, not modules) that show no interactions.

-}

-\item{geneToModule}{

+  (only genes, not modules) that show no interactions. These genes

+      cannot be part of modules. But you can specify modules that have

+      no epistatic interactions. See examples and vignette.

+    }

+    

+    \item{geneToModule}{

+      

   A named character vector that allows to match genes and modules. The

   names are the modules, and each of the values is a character vector

   with the gene names, separated by a comma, that correspond to a

   module. Note that modules cannot share genes. There is no need for

-  modules to contain more than one gene. If you specify a geneToModule

-  argument, it must necessarily contain "Root".  

+  modules to contain more than one gene.  If you specify a geneToModule

+  argument, and you used a restriction table, the \code{geneToModule} 

+  must necessarily contain, in the first position, "Root" (since the

+  restriction table contains a node named "Root"). See examples below.

 }

+    

 \item{drvNames}{The names of genes that are considered drivers. This is

   only used for: a) deciding when to stop the simulations, in case you

@@ -207,6 +214,19 @@ fea <- allFitnessEffects(rT = p4, epistasis = epist, orderEffects = oe,

                          noIntGenes = noint, geneToModule = modules)

 plot(fea)

+

+

+## Modules that show, between them,

+## no epistasis (so multiplicative effects).

+## We specify the individual terms, but no value for the ":".

+

+fnme <- allFitnessEffects(epistasis = c("A" = 0.1,

+                                        "B" = 0.2),

+                          geneToModule = c("A" = "a1, a2",

+                                           "B" = "b1"))

+

+evalAllGenotypes(fnme, order = FALSE, addwt = TRUE)

+

 }

 \keyword{ manip }

@@ -67,11 +67,12 @@ For \code{evalGenotype} either the value of fitness or (if \code{verbose

 For \code{evalAllGenotypes} a data frame with two columns, the Genotype

 and the Fitness (or Death Rate, if Bozic). The notation for the Genotype

 colum is a follows: when order does not matter, a comma "," separates

-the identifier of mutated genes. When order matters, a genotype shown as

-``x > y _ z'' means that a mutation in ``x'' happened before a mutation

-in ``y''; there is also a mutation in ``z'', but ``z'' is a gene for

-which order does not matter. In all cases, a "WT" denotes the wild-type

-(or, actually, the genotype without any mutations).

+the identifiers of mutated genes. When order matters, a genotype shown

+as ``x > y _ z'' means that a mutation in ``x'' happened before a

+mutation in ``y''; there is also a mutation in ``z'' (which could have

+happened before or after either of ``x'' or ``y''), but ``z'' is a gene

+for which order does not matter. In all cases, a "WT" denotes the

+wild-type (or, actually, the genotype without any mutations).

 }

 \author{

@@ -19,7 +19,7 @@

   al., 2012.

 }

 \usage{

- oncoSimulIndiv(fp, model = "Exp", numPassengers = 30, mu = 1e-6,

+ oncoSimulIndiv(fp, model = "Exp", numPassengers = 0, mu = 1e-6,

                 detectionSize = 1e8, detectionDrivers = 4,

                 sampleEvery = ifelse(model \%in\% c("Bozic", "Exp"), 1,

                              0.025),

@@ -39,7 +39,7 @@

                 initMutant = NULL,

                 seed = NULL)

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

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

                 detectionSize = 1e8, detectionDrivers = 4,

                 sampleEvery = ifelse(model \%in\% c("Bozic", "Exp"), 1,

                              0.025),

@@ -73,7 +73,7 @@ oncoSimulSample(Nindiv,

                                     if(length(fp$drv)) {

                                         nd <- (2: round(0.75 * length(fp$drv)))

                                     } else {

-                                        nd <- 0

+                                        nd <- 9e6

                                     }

                                 } else {

                                     nd <- (2 : round(0.75 * max(fp)))

@@ -119,18 +119,24 @@ oncoSimulSample(Nindiv,

     v.1. Otherwise, a fitnessEffects object (see

     \code{\link{allFitnessEffects}}).

-    Other arguments below (s, sh) make sense only if you use a poset, as

-    they are included in the fitnessEffects object.

+    Other arguments below (s, sh, numPassengers) make sense only if you

+    use a poset, as they are included in the fitnessEffects object.

 }

+

 \item{model}{

   One of "Bozic", "Exp", "McFarlandLog" (the last one can be abbreviated

-  to "McFL").

+  to "McFL"). The default is "Exp".

 }

+

 \item{numPassengers}{

+

+  This has no effect if you use the \code{\link{allFitnessEffects}}

+  specification.

-  The number of passenger genes.The total number of genes (drivers plus

+  If you use the specification of v.1., the number of passenger

+  genes. Note that using v.1 the total number of genes (drivers plus

   passengers) must be smaller than 64.

   All driver genes should be included in the poset (even if they depend

@@ -140,25 +146,37 @@ oncoSimulSample(Nindiv,

   of passengers).

 }

+

 \item{mu}{

   Mutation rate. See also \code{mutationPropGrowth}.

 } \item{detectionSize}{ What is the minimal number of cells for cancer

 to be detected.  For \code{oncoSimulSample} this can be a vector.

-}

-\item{detectionDrivers}{

-  The minimal number of drivers present in any clone for cancer to be

-  detected. For \code{oncoSimulSample} this can be a vector. The default

-  in this case is a vector of drivers from a uniform between 2 and 0.75

-  the total number of drivers

+} 

+

+\item{detectionDrivers}{ The minimal number of drivers (not modules,

+drivers, whether or not they are from the same module) present in any

+clone for cancer to be detected. For \code{oncoSimulSample} this can be

+a vector.

+

+For \code{oncoSimulSample}, if there are drivers (either because you are

+using a v.1 object or because you are using a fitnessEffects object with

+a \code{drvNames} component ---see \code{\link{allFitnessEffects}}---)

+the default is a vector of drivers from a uniform between 2 and 0.75 the

+total number of drivers. If there are no drivers (because you are using

+a fitnessEffects object without a \code{drvNames}, either because you

+specified it explicitly or because all of the genes are in the

+noIntGenes component) the simulations should not stop based on the

+number of drivers (and, thus, the default is set to 9e6).

 }

+

 \item{sampleEvery}{

   How often the whole population is sampled. This is not the same as the

-  interval between successive samples that keep stored (for that, see

-  \code{keepEvery}).

+  interval between successive samples that are kept or stored (for that,

+  see \code{keepEvery}).

   For very fast growing clones, you might need to have a small value

   here to minimize possible numerical problems (such as huge increase in

@@ -393,14 +411,22 @@ to be detected.  For \code{oncoSimulSample} this can be a vector.

   you set it to "auto", then if you are using v.1, the behavior is the

   same as if you set it to NULL (a seed will be generated in R and

   passed to C++) but if you are using v.2, a random seed will be

-  produced in C++. %% using the randutils code from M.E.O'Neill.  If you

-  need reproducibility, either pass a value or set it to NULL (setting

+  produced in C++. %% using the randutils code from M.E.O'Neill.

+  If you   need reproducibility, either pass a value or set it to NULL (setting

   it to NULL will make the C++ seed reproducible if you use the same

   seed in R via \code{set.seed}). However, even using the same value of

   \code{seed} is unlikely to give the exact same results between

   platforms and compilers.  Moreover, note that the defaults for

   \code{seed} are not the same in \code{oncoSimulIndiv},

-  \code{oncoSimulPop} and \code{oncoSimulSample}. }

+  \code{oncoSimulPop} and \code{oncoSimulSample}.

+

+  When using oncoSimulPop, if you want reproducibility, you might want

+  to, in addition to setting \code{seed = NULL}, also do

+  \code{RNGkind("L'Ecuyer-CMRG")} as we use

+  \code{\link[parallel]{mclapply}}; look at the vignette of

+  \pkg{parallel}.

+

+}

@@ -440,8 +466,8 @@ to be detected.  For \code{oncoSimulSample} this can be a vector.

-    \code{GenotyesWDistinctOrderEff} provides the information about

-    order effects that is missing from \code{Genotyoes}. When there are

+    \code{GenotypesWDistinctOrderEff} provides the information about

+    order effects that is missing from \code{Genotypes}. When there are

     order effects, the \code{Genotypes} matrix can contain genotypes

     that are not distinguishable. Suppose there are two genes, the first

     and the second. In the \code{Genotype} output you can get two

@@ -449,10 +475,18 @@ to be detected.  For \code{oncoSimulSample} this can be a vector.

     correspond to the two possible orders (first gene mutated first, or

     first gene mutated after the

     second). \code{GenotypesWDistinctOrderEff} disambiguates this. The

-    same is done by \code{GenotypeLabels}; this is easier to decode for

+    same is done by \code{GenotypesLabels}; this is easier to decode for

     a human (a string of gene labels) but a little bit harder to parse

-    automatically. 

+    automatically.  Note that when you use the default print method for

+    this object, you get, among others, a two-column display with the

+    \code{GenotypeLabels} information. When order matters, a genotype

+    shown as ``x > y _ z'' means that a mutation in ``x'' happened

+    before a mutation in ``y''; there is also a mutation in ``z'' (which

+    could have happened before or after either of ``x'' or ``y''), but

+    ``z'' is a gene for which order does not matter. When order does not

+    matter, a comma "," separates the identifiers of mutated genes.

+    

 }

@@ -616,7 +650,7 @@ to be detected.  For \code{oncoSimulSample} this can be a vector.

 data(examplePosets)

 p701 <- examplePosets[["p701"]]

-## Bozic Model

+## Exp Model

 b1 <- oncoSimulIndiv(p701)

 summary(b1)

@@ -223,7 +223,13 @@

     clone is not among the non-shown data (i.e., so that we can see when

     each clone/driver originates).

-    Thining is done to reduce the plot size and to speed up plotting..

+    Thinning is done to reduce the plot size and to speed up plotting.

+

+    Note that thinning is carried out before dealing with the plot axis,

+    so the actual number of points to be plotted could be a lot less (if

+    you reduce the x-axis considerably) than those returned from the

+    thinning. (In extreme cases this could lead to crashes when trying

+    to use stream plots if, say, you end up plotting only three values).

   }

   \item{thinData.keep}{

@@ -161,7 +161,6 @@ inline void driverCounts(int& maxNumDrivers,

   // number of present drivers.

   // We used to do count_NumDrivers and then whichDrivers

-  

   maxNumDrivers = 0;

   int tmpdr = 0;

   int driver_indx;

@@ -175,11 +174,13 @@ inline void driverCounts(int& maxNumDrivers,

     }

     if(tmpdr > maxNumDrivers) maxNumDrivers = tmpdr;

   }

-  

   for(size_t i = 0; i < countByDriver.size(); ++i) {

-    if(countByDriver[i] > 0) presentDrivers.push_back(i + 1);

-    ++totalPresentDrivers;

+    if(countByDriver[i] > 0) {

+      presentDrivers.push_back(i + 1);

+      ++totalPresentDrivers;

+    }

   }

+  

 }

@@ -588,9 +589,10 @@ static void nr_sample_all_pop_P(std::vector<int>& sp_to_remove,

 #ifdef DEBUGV

     Rcpp::Rcout << "\n\n     ********* 5.9 ******\n " 

 	      << "     Species  = " << i 

-	      << "\n      Genotype = " << genotypeSingleVector(Genotypes[i])

-      //	      << "\n      sp_id = " << genotypeSingleVector(Genotypes[i]) // sp_id[i]  

-	      << "\n      pre-update popSize = " 

+		<< "\n      Genotype = ";

+    print_Genotype(Genotypes[i]); //genotypeSingleVector(Genotypes[i])

+    //	      << "\n      sp_id = " << genotypeSingleVector(Genotypes[i]) // sp_id[i]  

+    Rcpp::Rcout << "\n      pre-update popSize = " 

 	      << popParams[i].popSize 

 	      << "\n      time of sample = " << tSample 

 	      << "\n      popParams[i].timeLastUpdate = " 

@@ -620,7 +622,8 @@ static void nr_sample_all_pop_P(std::vector<int>& sp_to_remove,

 #ifdef DEBUGV

       Rcpp::Rcout << "\n\n     Removing species i = " << i 

-		  << " with genotype = " << genotypeSingleVector(Genotypes[i]);

+		  << " with genotype = ";

+      print_Genotype(Genotypes[i]); //genotypeSingleVector(Genotypes[i]);

 #endif

     } 

 #ifdef DEBUGV

@@ -1084,8 +1087,9 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 		  << "     tSample  = " << tSample

 		  << "\n\n**   Species  = " << u_1 

-		  << "\n       genotype =  " << Genotypes[u_1] 

-		  << "\n       popSize = " << popParams[u_1].popSize 

+		    << "\n       genotype =  ";

+	print_Genotype(Genotypes[u_1]);

+	Rcpp::Rcout << "\n       popSize = " << popParams[u_1].popSize 

 		  << "\n       currentTime = " << currentTime 

 		  << "\n       popParams[i].nextMutationTime = " 

 		  << tmpdouble1

@@ -1113,8 +1117,9 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 		  << "     tSample  = " << tSample

 		  << "\n\n**   Species  = " << u_1 

-		  << "\n       genotype =  " << Genotypes[u_1] 

-		  << "\n       popSize = " << popParams[u_1].popSize 

+		    << "\n       genotype =  ";

+	print_Genotype(Genotypes[u_1]);

+	Rcpp::Rcout << "\n       popSize = " << popParams[u_1].popSize 

 		  << "\n       currentTime = " << currentTime 

 		  << "\n       popParams[i].nextMutationTime = " 

 		  << tmpdouble1

@@ -1125,8 +1130,9 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 		  << "\n\n**     Species  = " << u_2 

-		  << "\n       genotype =  " << Genotypes[u_2] 

-		  << "\n       popSize = " << popParams[u_2].popSize 

+		    << "\n       genotype =  ";

+	print_Genotype(Genotypes[u_2]);

+	Rcpp::Rcout << "\n       popSize = " << popParams[u_2].popSize 

 		  << "\n       currentTime = " << currentTime 

 		  << "\n       popParams[i].nextMutationTime = " 

 		  << tmpdouble2

@@ -1137,7 +1143,7 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 #endif

-      } else { // we sampled, so update all

+      } else { // we sampled, so update all: i.e. to_update == 3

 	for(size_t i = 0; i < popParams.size(); i++) {

 	  tmpdouble1 = ti_nextTime_tmax_2_st(popParams[i],

 					     currentTime,

@@ -1148,8 +1154,9 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 #ifdef DEBUGV

 	  Rcpp::Rcout << "\n\n     ********* 5.2: call to ti_nextTime ******\n " 

 		    << "     Species  = " << i 

-		    << "\n       genotype =  " << Genotypes[i] 

-		    << "\n       popSize = " << popParams[i].popSize 

+		      << "\n       genotype =  ";

+	  print_Genotype(Genotypes[i]);

+	  Rcpp::Rcout << "\n       popSize = " << popParams[i].popSize 

 		    << "\n       currentTime = " << currentTime 

 	    // << "\n       popParams[i].nextMutationTime = " 

 	    // << popParams[i].nextMutationTime

@@ -1311,6 +1318,7 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	    W_f_st(tmpParam);

 	    R_f_st(tmpParam);

 	    tmpParam.timeLastUpdate = -99999.99999; //mapTimes_updateP does what it should.

+	    // as this is a new species

 	    popParams.push_back(tmpParam);

 	    Genotypes.push_back(newGenotype);

 	    to_update = 2;

@@ -1346,20 +1354,27 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	  // #endif

 	} else {	// A mutation to pre-existing species

 #ifdef DEBUGW

-	  if( (currentTime - popParams[sp].timeLastUpdate) <= 0.0)

-	    throw std::out_of_range("currentTime - timeLastUpdate out of range"); 

+	  if( (currentTime - popParams[sp].timeLastUpdate) <= 0.0) {

+	    DP2(currentTime);

+	    DP2(sp);

+	    DP2(popParams[sp].timeLastUpdate);

+	    print_spP(popParams[sp]);

+	    throw std::out_of_range("currentTime - timeLastUpdate out of range");

+	  }

 #endif

 	  // if(verbosity >= 2) {

 #ifdef DEBUGV

 	    Rcpp::Rcout <<"\n     Mutated to existing species " << sp 

-			<< " (Genotype = " << genotypeSingleVector(Genotypes[sp]) 

+			<< " (Genotype = ";

+	    print_Genotype(Genotypes[sp]); 

 	      // << "; sp_id = " << Genotypes[sp].to_ullong()

-			<< ")"

+	    Rcpp::Rcout << ")"

 			<< "\n from species "  <<   nextMutant

-			<< " (Genotypes = " << genotypeSingleVector(Genotypes[nextMutant]) 

+			<< " (Genotypes = ";

+	    print_Genotype(Genotypes[nextMutant]); 

 	      // << "; sp_id = " << Genotypes[sp].to_ullong()

-			<< ")";

+	    Rcpp::Rcout	<< ")";

 	    // }

 #endif

 	  // FIXME00: the if can be removed??

@@ -1372,9 +1387,14 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	  } else {

 	    throw std::range_error("\n popSize == 0 but existing? \n");

 	  }

-	

 #ifdef DEBUGW

-	  popParams[sp].timeLastUpdate = -99999.99999; // to catch errors

+	  // This is wrong!!! if we set it to -999999, then the time to

+  	  // next mutation will not be properly updated.  In fact, the

+  	  // mapTimes map becomes a mess because the former pv in the

+  	  // popParams is not removed so we end up inserting another pair

+  	  // for the same species.

+	  

+	  // popParams[sp].timeLastUpdate = -99999.99999; // to catch errors

 #endif

 	  //popParams[sp].Flag = true;

 	}

@@ -1396,8 +1416,7 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	// DP2(popParams[nextMutant].popSize);

 	// Rcpp::Rcout << "\n done null mutation; after popSize ********" << std::endl;

       }

-    }

-      else { //       *********** We are sampling **********

+    } else { //       *********** We are sampling **********

       to_update = 3; //short_update = false;

       if(verbosity >= 2) {

 	Rcpp::Rcout <<"\n We are SAMPLING";   

@@ -1552,7 +1571,9 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

-  if(K < 1 )

+  if( (K < 1 ) && ( (typeModel == TypeModel::mcfarlandlog) ||

+		    (typeModel == TypeModel::mcfarland) ||

+		    (typeModel == TypeModel::mcfarland0) ))

     throw std::range_error("K < 1.");

   fitnessEffectsAll fitnessEffects =  convertFitnessEffects(rFE);

   //Used at least twice

@@ -1434,8 +1434,13 @@ static void innerBNB(const int& numGenes,

 	  //#endif

 	} else {	// A mutation to pre-existing species

 #ifdef DEBUGW

-	  if( (currentTime - popParams[sp].timeLastUpdate) <= 0.0)

-	    throw std::out_of_range("currentTime - timeLastUpdate out of range"); 

+	  if( (currentTime - popParams[sp].timeLastUpdate) <= 0.0) {

+	    DP2(currentTime);

+	    DP2(sp);

+	    DP2(popParams[sp].timeLastUpdate);

+	    print_spP(popParams[sp]);

+	    throw std::out_of_range("currentTime - timeLastUpdate out of range");

+	  }

 #endif

 	  if(verbosity >= 2) {

@@ -1459,7 +1464,13 @@ static void innerBNB(const int& numGenes,

 	  }

 #ifdef DEBUGW

-	  popParams[sp].timeLastUpdate = -99999.99999; // to catch errors

+	   // This is wrong!!! if we set it to -999999, then the time to

+  	  // next mutation will not be properly updated.  In fact, the

+  	  // mapTimes map becomes a mess because the former pv in the

+  	  // popParams is not removed so we end up inserting another pair

+  	  // for the same species.

+

+	  // popParams[sp].timeLastUpdate = -99999.99999; // to catch errors

 #endif

 	  //popParams[sp].Flag = true;

 	}

@@ -7,14 +7,14 @@

 SEXP OncoSimulR_nr_BNB_Algo5(SEXP rFESEXP, SEXP muSEXP, SEXP deathSEXP, SEXP initSizeSEXP, SEXP sampleEverySEXP, SEXP detectionSizeSEXP, SEXP finalTimeSEXP, SEXP initSpSEXP, SEXP initItSEXP, SEXP seedSEXP, SEXP verbositySEXP, SEXP speciesFSSEXP, SEXP ratioForceSEXP, SEXP typeFitness_SEXP, SEXP maxramSEXP, SEXP mutationPropGrowthSEXP, SEXP initMutant_SEXP, SEXP maxWallTimeSEXP, SEXP keepEverySEXP, SEXP alphaSEXP, SEXP KSEXP, SEXP detectionDriversSEXP, SEXP onlyCancerSEXP, SEXP errorHitWallTimeSEXP, SEXP maxNumTriesSEXP, SEXP errorHitMaxTriesSEXP, SEXP minDetectDrvCloneSzSEXP, SEXP extraTimeSEXP, SEXP keepPhylogSEXP);

 SEXP OncoSimulR_BNB_Algo5(SEXP restrictTableSEXP, SEXP numDriversSEXP, SEXP numGenesSEXP, SEXP typeCBN_SEXP, SEXP birthRateSEXP, SEXP sSEXP, SEXP deathSEXP, SEXP muSEXP, SEXP initSizeSEXP, SEXP sampleEverySEXP, SEXP detectionSizeSEXP, SEXP finalTimeSEXP, SEXP initSpSEXP, SEXP initItSEXP, SEXP seedSEXP, SEXP verbositySEXP, SEXP speciesFSSEXP, SEXP ratioForceSEXP, SEXP typeFitness_SEXP, SEXP maxramSEXP, SEXP mutationPropGrowthSEXP, SEXP initMutantSEXP, SEXP maxWallTimeSEXP, SEXP keepEverySEXP, SEXP alphaSEXP, SEXP shSEXP, SEXP KSEXP, SEXP detectionDriversSEXP, SEXP onlyCancerSEXP, SEXP errorHitWallTimeSEXP, SEXP maxNumTriesSEXP, SEXP errorHitMaxTriesSEXP, SEXP minDetectDrvCloneSzSEXP, SEXP extraTimeSEXP);

-SEXP OncoSimulR_readFitnessEffects(SEXP rFESEXP, SEXP echoSEXP);

+// SEXP OncoSimulR_readFitnessEffects(SEXP rFESEXP, SEXP echoSEXP);

 SEXP OncoSimulR_evalRGenotype(SEXP rGSEXP, SEXP rFESEXP, SEXP verboseSEXP, SEXP prodNegSEXP);

 R_CallMethodDef callMethods[]  = {

-  {"C_OncoSimulR_nr_BNB_Algo5", (DL_FUNC) &OncoSimulR_nr_BNB_Algo5, 29},

-  {"C_OncoSimulR_BNB_Algo5", (DL_FUNC) &OncoSimulR_BNB_Algo5, 34},

-  {"C_OncoSimulR_readFitnessEffects", (DL_FUNC) &OncoSimulR_readFitnessEffects, 2},

-  {"C_OncoSimulR_evalRGenotype", (DL_FUNC) &OncoSimulR_evalRGenotype, 4},

+  {"OncoSimulR_nr_BNB_Algo5", (DL_FUNC) &OncoSimulR_nr_BNB_Algo5, 29},

+  {"OncoSimulR_BNB_Algo5", (DL_FUNC) &OncoSimulR_BNB_Algo5, 34},

+  //  {"OncoSimulR_readFitnessEffects", (DL_FUNC) &OncoSimulR_readFitnessEffects, 2},

+  {"OncoSimulR_evalRGenotype", (DL_FUNC) &OncoSimulR_evalRGenotype, 4},

   {NULL, NULL, 0}

 };

@@ -88,17 +88,19 @@ BEGIN_RCPP

     return __result;

 END_RCPP

 }

-// readFitnessEffects

-void readFitnessEffects(Rcpp::List rFE, bool echo);

-RcppExport SEXP OncoSimulR_readFitnessEffects(SEXP rFESEXP, SEXP echoSEXP) {

-BEGIN_RCPP

-    Rcpp::RNGScope __rngScope;

-    Rcpp::traits::input_parameter< Rcpp::List >::type rFE(rFESEXP);

-    Rcpp::traits::input_parameter< bool >::type echo(echoSEXP);

-    readFitnessEffects(rFE, echo);

-    return R_NilValue;

-END_RCPP

-}

+

+// // readFitnessEffects

+// void readFitnessEffects(Rcpp::List rFE, bool echo);

+// RcppExport SEXP OncoSimulR_readFitnessEffects(SEXP rFESEXP, SEXP echoSEXP) {

+// BEGIN_RCPP

+//     Rcpp::RNGScope __rngScope;

+//     Rcpp::traits::input_parameter< Rcpp::List >::type rFE(rFESEXP);

+//     Rcpp::traits::input_parameter< bool >::type echo(echoSEXP);

+//     readFitnessEffects(rFE, echo);

+//     return R_NilValue;

+// END_RCPP

+// }

+

 // evalRGenotype