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

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 

-Version: 2.9.6

-Date: 2017-12-27

+Version: 2.9.9

+Date: 2018-04-10

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

 		     email = "rdiaz02@gmail.com"),

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

@@ -98,7 +98,8 @@ oncoSimulSample <- function(Nindiv,

                          extraTime = extraTime,

                          minDetectDrvCloneSz = minDetectDrvCloneSz,

                          detectionSize = detectionSize,

-                         detectionDrivers = detectionDrivers)[, -1, drop = FALSE]

+                         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

@@ -672,11 +673,11 @@ oncoSimulIndiv <- function(fp,

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

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

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

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

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

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

+                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 character vector.")

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

         }

@@ -778,7 +779,8 @@ print.oncosimul <- function(x, ...) {

         cat("\n")

         cat("Final population composition:\n")

         df <- data.frame(Genotype = x$GenotypesLabels,

-                         N = x$pops.by.time[nrow(x$pops.by.time), -1])

+                         N = x$pops.by.time[nrow(x$pops.by.time), -1],

+                         stringsAsFactors = TRUE)

         print(df)

     }

 }

@@ -833,7 +835,7 @@ summary.oncosimulpop <- function(object, ...) {

             warning("All simulations failed.")

             return(NA)

         }

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

+    as.data.frame(rbindlist(tmp), stringsAsFactors = TRUE)

 }

@@ -1922,7 +1924,8 @@ OncoSimulWide2Long <- function(x) {

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

+                      Genotype = rep(genotLabels, rep(nr, nc)),

+                      stringsAsFactors = TRUE))

 }

@@ -83,7 +83,8 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

     dd <- data.frame(muts = mutated,

                      fitness = tfm$afe$Fitness,

-                     label = tfm$afe$Genotype)

+                     label = tfm$afe$Genotype,

+                     stringsAsFactors = TRUE)

     cl <- gaj[vv]

     sg <- data.frame(x_from = mutated[vv[, 1]],

                      y_from = tfm$afe$Fitness[vv[, 1]],

@@ -91,7 +92,8 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

                      y_to = tfm$afe$Fitness[vv[, 2]],

                      Change = factor(ifelse(cl == 0, "Neutral",

                                      ifelse(cl > 0, "Gain", "Loss")),

-                                     levels = c("Gain", "Loss", "Neutral")))

+                                     levels = c("Gain", "Loss", "Neutral")),

+                     stringsAsFactors = TRUE)

     ## From http://stackoverflow.com/a/17257422

     number_ticks <- function(n) {function(limits) pretty(limits, n)}

@@ -601,6 +603,7 @@ peak_valley <- function(x) {

 ##     ## mutated <- rowSums(gfm[, -ncol(gfm)])

 ##     gaj <- OncoSimulR:::genot_to_adj_mat(gfm)

 ##     gaj2 <- OncoSimulR:::filter_inaccessible(gaj, 0)

+##     Eh! that is assuming no genotype is inaccessible!

 ##     stopifnot(all(na.omit(as.vector(gaj == gaj2))))

 ##     remaining <- as.numeric(colnames(gaj2))

 ##     ## mutated <- mutated[remaining]

@@ -41,6 +41,20 @@ to_Magellan <- function(x, file,

     }

 }

+

+

+## ## genotype_fitness_matrix -> fitness landscape as data frame

+## fmatrix_to_afe <- function(x) {

+##     stopifnot(inherits(x, "genotype_fitness_matrix"))

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

+##     nn <- apply(y, 1,

+##                 function(u) paste(sort(colnames(y)[as.logical(u)]),

+##                                   collapse = ", "))

+##     nn[nn == ""] <- "WT"

+##     return(data.frame(Genotype = nn, Fitness = x[, ncol(x)],

+##            stringsAsFactors = FALSE))

+## }

+

 to_Fitness_Matrix <- function(x, max_num_genotypes) {

     ## A general converter. Ready to be used by plotFitnessLandscape and

     ## Magellan exporter.

@@ -54,6 +68,7 @@ to_Fitness_Matrix <- function(x, max_num_genotypes) {

         ## if( (is.null(colnames(x))) || any(grepl("^$", colnames(x))))

         ##    stop("Matrix x must have column names")

+        ## This could use fmatrix_to_afe, above!!!

         ## Major change as of flfast: no longer using from_genotype_fitness

         afe <- evalAllGenotypes(allFitnessEffects(

             genotFitness = x

@@ -63,7 +63,7 @@ check.gm <- function(gm) {

 }

 gtm2 <- function(x) {

-    data.frame(cbind(nice.vector.eo(x, ","), x))

+    data.frame(cbind(nice.vector.eo(x, ","), x), stringsAsFactors = TRUE)

 }

 ## nice.vector.eo <- function(z, sep) {

@@ -88,7 +88,8 @@ gm.to.geneModuleL <- function(gm, one.to.one) {

     gm <- check.gm(gm)

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

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

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

+                             stringsAsFactors = TRUE) ## this stringsAsFactors is key

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

     colnames(geneMod)[1] <- c("Gene")

     geneMod <- geneMod[, -2]

@@ -285,7 +286,8 @@ epist.order.to.pairs.modules <- function(x, sep, rm.sign = TRUE) {

 to.pairs.modules <- function(x, sep, rm.sign = TRUE) {

     df <- data.frame(rbindlist(

         lapply(names(x),

-               function(z) epist.order.to.pairs.modules(z, sep))))

+               function(z) epist.order.to.pairs.modules(z, sep))),

+        stringsAsFactors = TRUE)

     if(nrow(df) == 0L) { ## if only single genes in epist, we get nothing here.

         return(df)

     } else {

@@ -1780,6 +1782,9 @@ nr_oncoSimul.internal <- function(rFE,

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

                                   fixation = NULL ## avoid partial matching

                                   ) {

+

+    default_min_successive_fixation <- 50 ## yes, set at this for now

+    

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

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

                    "as created, for instance, with function",

@@ -1939,11 +1944,75 @@ nr_oncoSimul.internal <- function(rFE,

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

     ## call <- match.call()

-    

     ## Process the fixed list, if any

     if(!is_null_na(fixation)) {

         ng <- namedGenes

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

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

+        ## Proposed extension to have exact matching of genotypes

+        ng <- rbind(

+             data.frame(Gene = "_", GeneNumID = -9, stringsAsFactors = FALSE),

+             ng)

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

+        ## FIXME

+        ## Later, accept a last argument, called tolerance.

+        ## If not present, set to 0

+        ## and then, at at the head of fixation_list below

+        if(is.list(fixation)) {

+            if(

+            (is.null(fixation[["fixation_tolerance"]])) ||

+            (is.na(fixation[["fixation_tolerance"]]))) {

+                fixation_tolerance <- 0

+            } else {

+                fixation_tolerance <- as.numeric(fixation[["fixation_tolerance"]])

+                fixation <- fixation[-which(names(fixation) == "fixation_tolerance")]

+            }

+            if(

+            (is.null(fixation[["min_successive_fixation"]])) ||

+            (is.na(fixation[["min_successive_fixation"]]))) {

+                min_successive_fixation <- default_min_successive_fixation

+            } else {

+                min_successive_fixation <- as.integer(fixation[["min_successive_fixation"]])

+                fixation <- fixation[-which(names(fixation) == "min_successive_fixation")]

+            }

+            if(

+            (is.null(fixation[["fixation_min_size"]])) ||

+            (is.na(fixation[["fixation_min_size"]]))) {

+                fixation_min_size <- 0

+            } else {

+                fixation_min_size <- as.integer(fixation[["fixation_min_size"]])

+                fixation <- fixation[-which(names(fixation) == "fixation_min_size")]

+            }

+            

+        } else {

+            if(is_null_na(fixation["fixation_tolerance"])) {

+                fixation_tolerance <- 0

+            } else {

+                fixation_tolerance <- as.numeric(fixation["fixation_tolerance"])

+                fixation <- fixation[-which(names(fixation) == "fixation_tolerance")]

+            }

+            if(is_null_na(fixation["min_successive_fixation"])) {

+                min_successive_fixation <- default_min_successive_fixation

+            } else {

+                min_successive_fixation <- as.integer(fixation["min_successive_fixation"])

+                fixation <- fixation[-which(names(fixation) == "min_successive_fixation")]

+            }

+            if(is_null_na(fixation["fixation_min_size"])) {

+                fixation_min_size <- 0

+            } else {

+                fixation_min_size <- as.integer(fixation["fixation_min_size"])

+                fixation <- fixation[-which(names(fixation) == "fixation_min_size")]

+            }

+        }

+

+        if( (fixation_tolerance > 1) || (fixation_tolerance < 0) )

+                    stop("Impossible range for fixation tolerance")

+

+        if( (min_successive_fixation < 0) )

+                    stop("Impossible range for min_successive_fixation")

+

+        if( (fixation_min_size < 0) )

+                    stop("Impossible range for fixation_min_size")

+

         ## Usual genotype specification and might allow ordered vectors

         ## in the future

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

@@ -1956,6 +2025,10 @@ nr_oncoSimul.internal <- function(rFE,

                        " in the fitness effects."))

         ## Sorting here is crucial!!

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

+        fixation_list <- list(fixation_list = fixation_list,

+                              fixation_tolerance = fixation_tolerance,

+                              min_successive_fixation = min_successive_fixation,

+                              fixation_min_size = fixation_min_size)

     } else {

         fixation_list <- list()

     }

@@ -2219,8 +2292,8 @@ detectionProbCheckParse <- function(x, initSize, verbosity) {

     if(x["checkSizePEvery"] <= 0)

         stop("checkSizePEvery <= 0")

-    if(x["PDBaseline"] < 0)

-        stop("PDBaseline < 0")

+    if(x["PDBaseline"] <= 0)

+        stop("PDBaseline <= 0")

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

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

@@ -2257,7 +2330,8 @@ sampledGenotypes <- function(y, genes = NULL) {

     ## ))

     df <- data.frame(table(

         apply(y, 1, genotlabel),

-        useNA = "ifany"))

+        useNA = "ifany"),

+        stringsAsFactors = TRUE)

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

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

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

     ## FIXME future: do this with order too?

     ##    - do not generate a matrix of genotypes but a matrix of ordered genot.

     wt_is_1 = match.arg(wt_is_1)

-    

+    if(is_null_na(g)) stop("Number of genes argument (g) is null or NA")

     m <- generate_matrix_genotypes(g)

     done <- FALSE

     ## attempts <- 0 ## for debugging/tracking purposes

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

+Changes in version 2.9.9 (2018-04-10):

+	- probDetect mechanism changed. This could be a BREAKING CHANGE.

+	  The expression divides by the baseline. For fixed initSize, this

+	  is simply a matter of changing the cPDetect.

+

+Changes in version 2.9.8 (2018-03-26):

+	- fixation allows exact genotypes, includes tolerance,

+	and checks for a successive number of specified periods

+

+Changes in version 2.9.7 (2018-02-20):

+	- fixed crash in some conditions when run with

+	stringsAsFactors = FALSE as global option

+	

 Changes in version 2.9.6 (2017-12-27):

 	- Updated citation.

 	- An example (in miscell-files) about using and stopping with

new file mode 100644

Binary files /dev/null and b/inst/testdata_fee.RData differ

@@ -442,20 +442,18 @@ of \code{sampleEvery} that is larger than or equal to \code{keepEvery}.

   probability of reaching cancer for a set of simulations to be

   accepted.}

-\item{max.wall.time}{

-  Maximum wall time for each individual simulation run. If the

-  simulation is not done in this time, it is aborted.

-

+\item{max.wall.time}{ Maximum wall time for the simulation of one

+  individual (over all \code{max.num.tries}). If the simulation is

+  not done in this time, it is aborted.  %% (If \code{max.num.tries = k},

+  %% time is kept track of over the k attempts;)

 }

 \item{max.wall.time.total}{

-

   Maximum wall time for all the simulations (when using

   \code{oncoSimulSample}), in seconds. If the simulation is not

   completed in this time, it is aborted. To prevent problems from a

   single individual simulation going wild, this limit is also enforced

   per simulation (so the run can be aborted directly from C++).

-

 }

 \item{errorHitMaxTries}{ If TRUE (the default) a simulation that reaches

@@ -522,13 +520,44 @@ of the individual done, and the number of attempts and time used.}

   with \code{fixation} is not allowed (as it does not make much sense).}

 \item{fixation}{If non-NULL, a list or a vector, where each element of

-  is a string with a gene or a gene combination. Simulations will stop

-  as soon as any of the genes or gene combinations are fixed (i.e.,

-  reach frequency of one). If you pass gene combinations, separate genes

+  is a string with a gene or a gene combination or a genotype (see

+  below). Simulations will stop as soon as any of the genes or gene

+  combinations or genotypes are fixed (i.e., reach a minimal

+  frequency). If you pass gene combinations or genotypes, separate genes

   with commas (not '>'); this means order is not (yet?)  supported. This

   way of specifying gene combinations is the same as the one used for

   \code{initMutant} and \code{\link{evalGenotype}}.

+  To differentiate between gene combinations and specific genotypes,

+  genotypes are specified by prepending them with a "_,". For instance,

+  \code{fixation = c("A", "B, C")} specifies stopping on any genotypes

+  with those gene combinations. In contrast, \code{fixation = c("_,A",

+  "_,B, C" )} specifies stopping only on gentoypes "A" or "B, C". See

+  the vignette for further examples.

+  

+

+  In addition to the gene combinations or genotypes themeselves, you can

+  add to the list or vector the named elements

+  \code{fixation_tolerance}, \code{min_successive_fixation} and

+  \code{fixation_min_size}. \code{fixation_tolerance}: fixation is

+  considered to have happened if the genotype/gene combinations

+  specified as genotypes/gene combinations for fixation have reached a

+  frequency \code{> 1 - fixation_tolerance}. (The default is 0, so we

+  ask for genotypes/gene combinations with a frequency of 1, which might

+  not be what you want with large mutation rates and complex fitness

+  landscape with genotypes of similar

+  fitness.). \code{min_successive_fixation}: during how many successive

+  sampling periods the conditions of fixation need to be fulfilled

+  before declaring fixation. These must be successive sampling periods

+  without interruptions (i.e., a single period when the condition is not

+  fulfilled will set the counter to 0). This can help to exclude short,

+  transitional, local maxima that are quickly replaced by other

+  genotypes. (The default is 50, but this is probably too small for

+  ``real life'' usage). \code{fixation_min_size}: you might only want to

+  consider fixation to have happened if a minimal size has been reached

+  (this can help weed out local maxima that have fitness that is barely

+  above that of the wild-type genotype). (The default is 0).

+

   Using this option with \code{AND_DrvProbExit} is not allowed (as it

   does not make much sense). This option is not allowed either with the

   old v.1 specification.}

@@ -567,7 +596,7 @@ of the individual done, and the number of attempts and time used.}

   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

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

   \pkg{parallel}.

 }

@@ -641,7 +670,7 @@ of the individual done, and the number of attempts and time used.}

   uniform number). Probability of detection changes with population size

   according to the function

-  \deqn{ 1 - e^{-cPDetect (populationsize - PDBaseline)}}.

+  \deqn{ 1 - e^{-cPDetect ( (populationsize - PDBaseline)/PDBaseline )}}.

   You can pass \code{cPDetect} manually (you will need to set \code{n2}

   and \code{p2} to NA). However, it might be more intuitive to specify

@@ -73,7 +73,7 @@ mean \code{mu} and standard deviation \code{sd}).}

   it is up to you to make sure that the range of the scale argument

   includes 1 or you might not get what you want). Note that using this

   option can easily lead to landscapes with no accessible genotypes

-  (unless you also use \code{scale}).

+  (even if you also use \code{scale}).

   If "none", the fitness of the wildtype is not touched.  }

@@ -230,6 +230,10 @@ void nr_totPopSize_and_fill_out_crude_P(int& outNS_i,

 					std::mt19937& ran_gen,

 					const bool& AND_DrvProbExit,

 					const std::vector<std::vector<int> >& fixation_l,

+					const double& fixation_tolerance,

+					const int& min_successive_fixation,

+					const double& fixation_min_size,

+					int& num_successive_fixation,

 					POM& pom,

 					const std::map<int, std::string>& intName,

 					const fitness_as_genes& genesInFitness,

@@ -275,23 +279,55 @@ void nr_totPopSize_and_fill_out_crude_P(int& outNS_i,

 	std::vector<int> thisg = allGenesinGenotype(Genotypes[i]);

 	for(size_t fc = 0; fc != popSize_fixation.size(); ++fc) {

 	  // Yes, fixation_l is sorted in R.

-	  if(std::includes(thisg.begin(), thisg.end(),

-			   fixation_l[fc].begin(), fixation_l[fc].end()) ) {

-	    popSize_fixation[fc] += popParams[i].popSize;

+	  // if fixation_l[fc] starts with a -9, we are asking

+	  // for exact genotype equality

+	  if(fixation_l[fc][0] == -9) {

+	    // // exact genotype identity?

+	    std::vector<int> this_fix(fixation_l[fc].begin() + 1,

+				      fixation_l[fc].end());

+	    if(thisg == this_fix) {

+	      popSize_fixation[fc] = popParams[i].popSize;

+	    }

+	  } else {	  

+	  // gene combination in fixation element present in genotype?

+	    if(std::includes(thisg.begin(), thisg.end(),

+			     fixation_l[fc].begin(), fixation_l[fc].end()) ) {

+	      popSize_fixation[fc] += popParams[i].popSize;

+	    }

 	  }

 	}

       }

       // Any fixated? But avoid trivial of totPopSize of 0!

       // Now check of > 0 is redundant as we check totPopSize > 0

+      // FIXME do we want tolerance around that value? 

       double max_popSize_fixation =

 	*std::max_element(popSize_fixation.begin(), popSize_fixation.end());

-      if( (max_popSize_fixation > 0 ) &&

-	  (max_popSize_fixation == totPopSize)) {

-	fixated = true;

+      if( (max_popSize_fixation >= fixation_min_size ) &&

+	  (max_popSize_fixation >= (totPopSize * (1 - fixation_tolerance) )) ) {

+	++num_successive_fixation;

+	// DP1("increased num_successive_fixation");

+	if( num_successive_fixation >= min_successive_fixation) fixated = true;

+      } else {

+	// DP1("zeroed num_successive_fixation");

+	num_successive_fixation = 0;

       }

     }

   }

+

+  // // DEBUG

+  // if(fixated) {

+  //   // print fixation_l

+  //   // print popSize_fixation

+  //   // print totPopSize

+  //   DP1("popSize_fixation");

+  //   for(size_t fc = 0; fc != popSize_fixation.size(); ++fc) {

+  //     DP2(fc);

+  //     DP2(popSize_fixation[fc]);

+  //   }

+  //   DP2(totPopSize);

+  // }

+  

   if (keepEvery < 0) {

     storeThis = false;

   } else if( currentTime >= (lastStoredSample + keepEvery) ) {

@@ -326,7 +362,20 @@ void nr_totPopSize_and_fill_out_crude_P(int& outNS_i,

     checkSizePNow = false;

   }

+  // We do not verify that conditions for exiting are also satisfied

+  // at the exit time when extraTime > 0. We could do that,

+  // checking again for the conditions (or the reasonable conditions, so

+  // probably not detectSizeP). For instance, with fixated.

+

+  // For fixated in particular, note that we evaluate fixation always, but

+  // we might be exiting when there is no longer fixation.  But the logic

+  // with fixation is probably to use as large a min_successive_fixation

+  // as desired and no extraTime.

+

+  // Probably would not need to check lastMaxDr and popSizeOverDDr

+  // as those should never decrease. Really?? FIXME

+  

   if(AND_DrvProbExit) {

     // The AND of detectionProb and drivers

     // fixated plays no role here, and cannot be passed from R

@@ -362,8 +411,12 @@ void nr_totPopSize_and_fill_out_crude_P(int& outNS_i,

 	  done_at = currentTime + extraTime;

 	}

       } else if (currentTime >= done_at) {

-  	simulsDone = true;

-  	reachDetection = true; 

+	// if(fixated) {

+	  simulsDone = true;

+	  reachDetection = true;

+	// } else {

+	//   done_at = -1;

+	// }

       }

     } else if( (fixated) ||

 	       (totPopSize >= detectionSize) ||

@@ -911,6 +964,9 @@ static void nr_innerBNB (const fitnessEffectsAll& fitnessEffects,

 			const double& checkSizePEvery,

 			const bool& AND_DrvProbExit,

 			const std::vector< std::vector<int> >& fixation_l,

+			 const double& fixation_tolerance,

+			 const int& min_successive_fixation,

+			 const double& fixation_min_size,

 			 int& ti_dbl_min,

 			 int& ti_e3,

 			 std::map<std::string, std::string>& lod,

@@ -1089,6 +1145,10 @@ static void nr_innerBNB (const fitnessEffectsAll& fitnessEffects,

   int lastMaxDr = 0;

   double done_at = -9;

+

+  int num_successive_fixation = 0; // none so far

+

+  

 #ifdef MIN_RATIO_MUTS_NR

   g_min_birth_mut_ratio_nr = DBL_MAX;

   g_min_death_mut_ratio_nr = DBL_MAX;

@@ -1881,6 +1941,10 @@ static void nr_innerBNB (const fitnessEffectsAll& fitnessEffects,

 					 ran_gen,

 					 AND_DrvProbExit,

 					 fixation_l,

+					 fixation_tolerance,

+					 min_successive_fixation,

+					 fixation_min_size,

+					 num_successive_fixation,

 					 pom, intName,

 					 genesInFitness); //keepEvery is for thinning

       if(verbosity >= 3) {

@@ -2011,6 +2075,7 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

   double cPDetect = cPDetect_i;

   if( (n2 > 0) && (p2 > 0) ) {

+    if (PDBaseline <= 0) throw std::range_error("PDBaseline <= 0");

     cPDetect = set_cPDetect(n2, p2, PDBaseline);

     if(verbosity >= 1)

       Rcpp::Rcout << "  cPDetect set at " << cPDetect << "\n";

@@ -2041,11 +2106,22 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

   }

   // fixation: run until some genotype combinations fixed

-  std::vector < std::vector<int> > fixation_l(fixation_i.size());

+  

+  double fixation_tolerance = -9;

+  int min_successive_fixation = 100;

+  double fixation_min_size = 0.0;

+  std::vector < std::vector<int> > fixation_l;

+

   if( fixation_i.size() != 0 ) {

-    fixation_l = list_to_vector_of_int_vectors(fixation_i);

+    Rcpp::List fggl = fixation_i["fixation_list"] ;

+    fixation_l = list_to_vector_of_int_vectors(fggl); // FIXME

+    fixation_tolerance = Rcpp::as<double>(fixation_i["fixation_tolerance"]);

+    min_successive_fixation = Rcpp::as<int>(fixation_i["min_successive_fixation"]);

+    fixation_min_size = Rcpp::as<double>(fixation_i["fixation_min_size"]);

+  } else {

+    fixation_l.resize(0); // explicit

   }

-

+ 

   bool runAgain = true;

   bool reachDetection = false;

@@ -2221,6 +2297,9 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

 		  checkSizePEvery,

 		  AND_DrvProbExit,

 		  fixation_l,

+		  fixation_tolerance,

+		  min_successive_fixation,

+		  fixation_min_size,

 		  ti_dbl_min,

 		  ti_e3,

 		  lod,

@@ -298,6 +298,11 @@ double ti_nextTime_tmax_2_st(const spParamsP& spP,

 	// throw std::range_error("ti set to DBL_MIN");

 	// Even just touching DBL_MIN is enough to want a rerunExcept;

 	// no need for it to be 0.0.

+	

+	// Trying to understand what happens

+	Rcpp::Rcout << "         ti set to DBL_MIN: spP.popSize = " << spP.popSize << "\n";

+	// It seems poSize over 1e8, and even 3.5 e7 can trigger this exception (depending

+	// on mutation rate, of course)

 	throw rerunExcept("ti set to DBL_MIN");

       }

       if(ti < (2*DBL_MIN)) ++ti_e3; // Counting how often this happens.

@@ -74,9 +74,12 @@ double prodMuts(const std::vector<double>& s) {

 		    std::multiplies<double>());

 }

+

+

+// new cPDetect

 double set_cPDetect(const double n2, const double p2,

 		    const double PDBaseline) {

-  return (-log(1.0 - p2)/(n2 - PDBaseline));

+  return ( -log(1.0 - p2) *  (PDBaseline / (n2 - PDBaseline)) );

 }

 // Next two are identical, except for scaling the k. Use the simplest.

@@ -85,19 +88,37 @@ double probDetectSize(const double n, const double cPDetect,

   if(n <= PDBaseline) {

     return 0;

   } else {

-    return (1 - exp( -cPDetect * (n - PDBaseline)));

+    return (1 - exp( -cPDetect * ( (n - PDBaseline)/PDBaseline ) ));

   }

 }

-// double prob_exit_ratio(const double n, const double k, const double baseline) {

-//   if(n <= baseline) {

+

+// Former cpDetect mechanism

+

+// double set_cPDetect(const double n2, const double p2,

+// 		    const double PDBaseline) {

+//   return (-log(1.0 - p2)/(n2 - PDBaseline));

+// }

+

+// // Next two are identical, except for scaling the k. Use the simplest.

+// double probDetectSize(const double n, const double cPDetect,

+// 		      const double PDBaseline) {

+//   if(n <= PDBaseline) {

 //     return 0;

 //   } else {

-//     return (1 - exp( -c * ( (n - baseline)/baseline)));

+//     return (1 - exp( -cPDetect * (n - PDBaseline)));

 //   }

 // }

+// // double prob_exit_ratio(const double n, const double k, const double baseline) {

+// //   if(n <= baseline) {

+// //     return 0;

+// //   } else {

+// //     return (1 - exp( -c * ( (n - baseline)/baseline)));

+// //   }

+// // }

+// is this detected, by the probability of detection as a function of size?

 bool detectedSizeP(const double n, const double cPDetect,

 			const double PDBaseline, std::mt19937& ran_gen) {

   if(cPDetect < 0) {

@@ -115,6 +136,8 @@ bool detectedSizeP(const double n, const double cPDetect,

   }

 }

+

+

 bool operator==(const Genotype& lhs, const Genotype& rhs) {

   return (lhs.orderEff == rhs.orderEff) &&

     (lhs.epistRtEff == rhs.epistRtEff) &&

@@ -1,5 +1,11 @@

 ## RNGkind("Mersenne-Twister")

-cat(date())

+

+cat(paste("\n Starting exercise-plotting-code long at", date()))

+cat(paste("\n             a runif", runif(1), "\n"))

+cat(paste("\n             a runif", runif(1), "\n"))

+date()

+

+

 test_that("exercising plotClonePhylog", {

               data(examplesFitnessEffects)

               tmp <-  oncoSimulIndiv(examplesFitnessEffects[["o3"]],

@@ -26,6 +32,10 @@ test_that("exercising plotClonePhylog", {

               plotClonePhylog(tmp, N = 10, keepEvents = TRUE)

               ## Reaching the fixOverlap code

               plotClonePhylog(tmp, N = 0, timeEvents = TRUE)

+              expect_true(TRUE) ## dummy, to prevent the

+              ## attempt to apply non-function

+              ## If plotting failed, this would fail

+              expect_error(plotClonePhylog(c(1, 2)))

           })

 cat(date())

@@ -58,74 +68,12 @@ test_that("exercising the fitnessEffects plotting code", {

               plot(fp4m, "igraph", layout = igraph::layout.reingold.tilford, 

                    expandModules = TRUE, autofit = TRUE)

               plot(fp4m, expandModules = TRUE, autofit = TRUE)

+              expect_true(TRUE) ## dummy, to prevent the

+              ## attempt to apply non-function

+              ## If plotting failed, this would fail

           })

 date()

-

-

-

-date()

-test_that("stacked, stream, genotypes and some colors", {

-    data(examplesFitnessEffects)

-    max.tries <- 4

-    for(i in 1:max.tries) {

-      tmp <-  oncoSimulIndiv(examplesFitnessEffects[["o3"]],

-                             model = "McFL", 

-                             mu = 5e-5,

-                             detectionSize = 1e8, 

-                             detectionDrivers = 3,

-                             sampleEvery = 0.025,

-                             max.num.tries = 10,

-                             keepEvery = 5,

-                             initSize = 2000,

-                             finalTime = 3000,

-                             onlyCancer = FALSE, detectionProb = NA,

-                             keepPhylog = TRUE)

-      

-      if(nrow(tmp$pops.by.time) >= 5) {

-          break

-      } else {

-          cat("\n hummm.. had to run again in the plot")

-          if(i >= max.tries) {

-              print(tmp)

-            stop("stream will break")

-          }

-      }

-    }

-    plot(tmp, type = "stacked", show = "genotypes")

-    plot(tmp, type = "stream", show = "genotypes")

-      plot(tmp, type = "line", show = "genotypes")

-

-      plot(tmp, type = "stacked", show = "drivers")

-      plot(tmp, type = "stream", show = "drivers")

-      plot(tmp, type = "line", show = "drivers")

-

-      plot(tmp, type = "stacked", order.method = "max")

-      plot(tmp, type = "stacked", order.method = "first")

-

-      plot(tmp, type = "stream", order.method = "max")

-      plot(tmp, type = "stream", order.method = "first")

-

-      plot(tmp, type = "stream", stream.center = TRUE)

-      plot(tmp, type = "stream", stream.center = FALSE)

-

-      plot(tmp, type = "stream", stream.center = TRUE, log = "x")

-      plot(tmp, type = "stacked", stream.center = TRUE, log = "x")

-      

-      plot(tmp, type = "stacked", show = "genotypes",

-           breakSortColors = "random")

-      plot(tmp, type = "stream", show = "genotypes",

-           breakSortColors = "distave")

-

-      plot(tmp, type = "stacked", show = "genotypes", col = rainbow(9))

-      plot(tmp, type = "stream", show = "genotypes", col = rainbow(3))

-      plot(tmp, type = "line", show = "genotypes", col = rainbow(20))

-      

-})

-date()

-

-

-

 test_that("xlab, ylab, ylim, xlim can be passed", {

     data(examplePosets)

     p701 <- examplePosets[["p701"]]

@@ -233,6 +181,11 @@ test_that("xlab, ylab, ylim, xlim can be passed", {

          ylab = "ylab", ylim = c(-100, 1000),

          xlim = c(20, 70),

          plotDrivers = TRUE)

+    expect_true(TRUE) ## dummy, to prevent the

+    ## attempt to apply non-function

+    ## If plotting failed, this would fail

+    expect_error(plot(e1, type = "stremaitoooihoh"))

+    

 })

@@ -271,32 +224,108 @@ test_that("oncosimul v.1 objects and genotype plotting", {

     plot(p1, type = "stacked", show = "genotypes", thinData = TRUE)

     plot(p1, type = "stream", show = "genotypes", thinData = TRUE)

     plot(p1, type = "line", show = "genotypes", thinData = TRUE)

+    expect_true(TRUE) ## dummy, to prevent the

+    ## attempt to apply non-function

+    ## If plotting failed, this would fail

+    expect_error(plot(tmp, type = "linito"))

 })

 date()

+## The following are not run because of the weird issue

+## using test_dir. But I test for colors and type in the

+## usual, regular, testing (test.exercise-plotting-code.R)

-test_that("passing colors", {

-    data(examplePosets)

-    ## An object of class oncosimul

-    p705 <- examplePosets[["p705"]]

-    max.tries <- 4

-    for(i in 1:max.tries) {

-    p1 <- oncoSimulIndiv(p705)

-    if(nrow(p1$pops.by.time) >= 11) {

-            break

-    } else {

-        cat("\n hummm.. had to run again in the plot")

-        if(i >= max.tries) {

-            print(p1)

-            stop("stream will break")

-        }

-    }

-    }

-    

-    class(p1)

-    plot(p1, type = "stacked", show = "genotypes", col = rainbow(8))

-    plot(p1, type = "stream", show = "genotypes", col = rainbow(18))

-    plot(p1, type = "line", show = "genotypes", col = rainbow(3))

-})

+

+## test_that("passing colors", {

+##     data(examplePosets)

+##     ## An object of class oncosimul

+##     p705 <- examplePosets[["p705"]]

+##     max.tries <- 4

+##     for(i in 1:max.tries) {

+##     p1 <- oncoSimulIndiv(p705)

+##     if(nrow(p1$pops.by.time) >= 11) {

+##             break

+##     } else {

+##         cat("\n hummm.. had to run again in the plot")

+##         if(i >= max.tries) {

+##             print(p1)

+##             stop("stream will break")

+##         }

+##     }

+##     }

+##     ## class(p1)

+##     plot(p1, type = "stacked", show = "genotypes", thinData = TRUE)

+##     ## with newest testthat, the next make if fail with test_dir, but

+##     ## not if run from REPL. Go figure

+##     plot(p1, type = "stacked", show = "genotypes", col = rainbow(8))

+##     plot(p1, type = "stream", show = "genotypes", col = rainbow(18))

+##     plot(p1, type = "line", show = "genotypes", col = rainbow(3))

+##     expect_true(TRUE) ## dummy, to prevent the

+##     ## attempt to apply non-function

+##     ## If plotting failed, this would fail

+##     expect_error(plot(p1, type = "linito"))

+## })

+

+

+

+

+

+## date()

+## test_that("stacked, stream, genotypes and some colors", {

+##     data(examplesFitnessEffects)

+##     max.tries <- 4

+##     for(i in 1:max.tries) {

+##         tmp <-  oncoSimulIndiv(examplesFitnessEffects[["o3"]],

+##                                model = "McFL", 

+##                                mu = 5e-5,

+##                                detectionSize = 1e8, 

+##                                detectionDrivers = 3,