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

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 

-Version: 2.7.0

-Date: 2017-04-07

+Version: 2.7.2

+Date: 2017-09-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"))

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

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

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

-                            "at which there are any mutants.",

+                            "at which there are any mutants with drivers. ",

                             "Thus, no uniform sampling possible.",

                             "You will get NAs"))

               the.time <- -99

@@ -21,6 +21,19 @@ accessibleGenotypes <- function(y, f, numMut, th) {

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

 }

+

+genot2AdjMat <- function(y, f, numMut) {

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

+}

+

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

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

+}

+

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

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

+}

+

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

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

 ## }

@@ -18,6 +18,9 @@

 ## FIXME: show only accessible paths? 

+## FIXME: if using only_accessible, maybe we

+## can try to use fast_peaks, and use the slower

+## approach as fallback (if identical fitness)

 plotFitnessLandscape <- function(x, show_labels = TRUE,

                                  col = c("green4", "red", "yellow"),

                                  lty = c(1, 2, 3), use_ggrepel = FALSE,

@@ -224,13 +227,45 @@ filter_inaccessible <- function(x, accessible_th) {

 }

+## wrapper to the C++ code

+fast_peaks <- function(x, th) {

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

+

+    ## Only works if no connected genotypes that form maxima. I.e., no

+    ## identical fitness. Do a sufficient check for it (too inclusive, though)

+    ## And only under no back mutation

+

+    original_pos <- 1:nrow(x)

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

+    o_numMut <- order(numMut)

+    x <- x[o_numMut, ]

+    numMut <- numMut[o_numMut]

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

+    ## Two assumptions

+    stopifnot(numMut[1] == 0)

+    ## make sure no repeated in those that could be maxima

+    if(any(duplicated(f[f >= f[1]])))

+        stop("There could be several connected maxima genotypes.",

+             " This function is not safe to use.")

+    

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

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

+    original_pos <- original_pos[o_numMut]

+    return(sort(original_pos[peaksLandscape(y, f,

+                          as.integer(numMut),

+                          th)]))

+}

+

+

 ## wrapper to the C++ code

 wrap_accessibleGenotypes <- function(x, th) {

     ## x is the fitness matrix, not adjacency matrix

+    original_pos <- 1:nrow(x)

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

     o_numMut <- order(numMut)

     x <- x[o_numMut, ]

     numMut <- numMut[o_numMut]

+    original_pos <- original_pos[o_numMut]

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

     storage.mode(y) <- "integer"

@@ -238,7 +273,29 @@ wrap_accessibleGenotypes <- function(x, th) {

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

                                as.integer(numMut),

                                th)

-    return(acc[acc > 0])

+    ## return(acc[acc > 0])

+    return(sort(original_pos[acc[acc > 0]]))

+}

+

+

+## wrapper to the C++ code; the former one, only for testing. Remove

+## eventually FIXME

+wrap_accessibleGenotypes_former <- function(x, th) {

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

+    original_pos <- 1:nrow(x)

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

+    o_numMut <- order(numMut)

+    x <- x[o_numMut, ]

+    numMut <- numMut[o_numMut]

+    original_pos <- original_pos[o_numMut]

+    

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

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

+

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

+                                      as.integer(numMut),

+                                      th)

+    return(sort(original_pos[acc[acc > 0]]))

 }

 ## A transitional function

@@ -378,8 +435,10 @@ generate_matrix_genotypes <- function(g) {

 }

-

-genot_to_adj_mat <- function(x) {

+## The R version. See also the C++ one

+genot_to_adj_mat_R <- function(x) {

+    ## x is the fitness matrix

+    

     ## FIXME this can take about 23% of the time of the ggplot call.

     ## But them, we are quickly constructing a 2000*2000 matrix

     ## Given a genotype matrix, as given by allGenotypes_to_matrix, produce a

@@ -390,13 +449,16 @@ genot_to_adj_mat <- function(x) {

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

     ## Make sure sorted, so ancestors always before descendants

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

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

-    rm(rs0)

+    original_pos <- 1:nrow(x)

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

+    o_numMut <- order(numMut)

+    x <- x[o_numMut, ]

+    original_pos <- original_pos[o_numMut]

+    rm(numMut)

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

     f <- x[, ncol(x)]

-    rs <- rowSums(y)

+    rs <- rowSums(y) ## redo for paranoia; could have ordered numMut

     ## Move this to C++?

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

@@ -408,9 +470,51 @@ genot_to_adj_mat <- function(x) {

             if(sumdiff == 1) adm[i, j] <- (f[j] - f[i])

         }

     }

+    colnames(adm) <- rownames(adm) <- original_pos

     return(adm)

 }

+genot_to_adj_mat <- function(x) {

+    ## x is the fitness matrix

+

+    ## adding column and row names should rarely be necessary

+    ## as these are internal functions, etc. But just in case

+    original_pos <- 1:nrow(x)

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

+    o_numMut <- order(numMut)

+    x <- x[o_numMut, ]

+    numMut <- numMut[o_numMut]

+    original_pos <- original_pos[o_numMut]

+    

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

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

+

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

+                        as.integer(numMut))

+    colnames(adm) <- rownames(adm) <- original_pos

+    return(adm)

+}

+

+

+## ## to move above to C++ note that loop can be

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

+##     for(j in (i:length(rs))) {

+##         if(rs[j] > (rs[i] + 1)) break;

+##         else if(rs[j] == (rs[i] + 1)) {

+##             ## and use here my HammingDistance function

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

+##             ## if(sumdiff == 1) adm[i, j] <- (f[j] - f[i])

+##             if(HammingDistance(y[j, ], y[i, ]) == 1) adm[i, j] = (f[j] - f[i]);

+##             }

+##     }

+## }

+

+## actually, all that is already in accessibleGenotypes except for the

+## filling up of adm.

+

+

+

+

 peak_valley <- function(x) {

     ## FIXME: when there are no identical entries, all this

@@ -466,7 +570,13 @@ peak_valley <- function(x) {

     }

     bad_fwd <- vector("integer", nrow(x))

     for(i in 1:nrow(x)) {

+        ## Eh, why any? All.

+        ## Nope, any: we want peaks in general, not just

+        ## under assumption of "no back mutation"

+        ## We get a different result when we restrict to accessible

+        ## because all < 0 in adjacency are turned to NAs.

         if( any(x[, i] < 0, na.rm = TRUE) || bad_fwd[i] ) {

+        ## if( all(x[, i] < 0, na.rm = TRUE) ) {

             ## this node is bad. Any descendant with fitness >= is bad

             bad_fwd[i] <- 1

             reach_f <- which(x[i, ] <= 0)

@@ -478,3 +588,36 @@ peak_valley <- function(x) {

     peak <- setdiff(candidate, bad)

     return(list(peak = peak, valley = valley))

 }

+

+

+

+## For the future

+## ## data.frame (two columns: genotype with "," and Fitness) -> fitness graph (DAG)

+## ## Return an adj matrix of the fitness graph from a fitness

+## ## landscape

+## ## Based on code in plotFitnessLandscape

+## flandscape_to_fgraph <- function(afe) {

+##     gfm <- OncoSimulR:::allGenotypes_to_matrix(afe)

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

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

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

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

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

+##     ## mutated <- mutated[remaining]

+##     afe <- afe[remaining, , drop = FALSE]

+##     ## vv <- which(!is.na(gaj2), arr.ind = TRUE)

+

+##     gaj2 <- gaj2

+##     gaj2[is.na(gaj2)] <- 0

+##     gaj2[gaj2 > 0] <- 1

+##     colnames(gaj2) <- rownames(gaj2) <- afe[, "Genotype"]

+##     return(gaj2)

+## }

+## ## This could be done easily in C++, taking care of row/colnames at end,

+## ## without moving around the full adjacency matrix.

+## ## Skeleton for C++

+## ## a call to accessibleGenotypesPeaksLandscape

+## ## (with another argument or changing the returnpeaks by a three value thing)

+## ## after done with first loop, 

+## ## return the matrix adm[accessible > 0, accessible >0]

+## ## only need care with row/colnames

@@ -23,7 +23,7 @@

 to_Magellan <- function(x, file,

                         max_num_genotypes = 2000) {

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

+    ## Go directly if you have, as entry, an object from

     ## rfitness!! to_Fitness_Matrix can be very slow.

     if(is.null(file)) {

         file <- tempfile()

@@ -44,7 +44,8 @@ to_Magellan <- function(x, file,

 to_Fitness_Matrix <- function(x, max_num_genotypes) {

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

     ## Magellan exporter.

-    

+

+    ## FIXME: really, some of this is inefficient. Very. Fix it.

     if( (inherits(x, "genotype_fitness_matrix")) ||

         ( (is.matrix(x) || is.data.frame(x)) && (ncol(x) > 2) ) ) {

         ## Why this? We go back and forth twice. We need both things. We

@@ -257,7 +258,7 @@ allGenotypes_to_matrix <- function(x) {

     splitted_genots <- lapply(x$Genotype,

                              function(z) nice.vector.eo(z, ","))

-    all_genes <- unique(unlist(splitted_genots))

+    all_genes <- sort(unique(unlist(splitted_genots)))

     m <- matrix(0, nrow = length(splitted_genots), ncol = length(all_genes))

     the_match <- lapply(splitted_genots,

@@ -285,7 +286,8 @@ Magellan_stats <- function(x, max_num_genotypes = 2000,

                            verbose = FALSE,

                            use_log = TRUE,

                            short = TRUE,

-                           fl_statistics = "fl_statistics") { # nocov start

+                           fl_statistics = "fl_statistics",

+                           replace_missing = FALSE) { # nocov start

     ## I always use 

     ## if(!is.null(x) && is.null(file))

     ##     stop("one of object or file name")

@@ -306,23 +308,43 @@ Magellan_stats <- function(x, max_num_genotypes = 2000,

         shortarg <- NULL

     }

+    if(replace_missing) {

+        zarg <- "-z"

+    } else {

+        zarg <- NULL

+    }

+

     to_Magellan(x, fn, max_num_genotypes = max_num_genotypes)

-    call_M <- system2(fl_statistics, args = paste(fn, shortarg, logarg, "-o", fnret))

+    ## newer versions of Magellan provide some extra values to standard output

+    call_M <- system2(fl_statistics,

+                      args = paste(fn, shortarg, logarg, zarg, "-o", fnret),

+                      stdout = NULL)

     if(short) {

-        tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[-1])

+        ## tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[-1])

+        

+        tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[c(-1)])

         ## Make names more explicit, but check we have what we think we have

-        stopifnot(length(tmp) == 23)

-        stopifnot(identical(names(tmp),

+        ## New versions of Magellan produce different output apparently of variable length

+        stopifnot(length(tmp) >= 23) ## 23) ## variable length

+        stopifnot(identical(names(tmp)[1:13], ## only some

                             c("ngeno", "npeaks", "nsinks", "gamma", "gamma.", "r.s",

                               "nchains", "nsteps", "nori", "depth", "magn", "sign",

-                              "rsign", "w.1.", "w.2.", "w.3..", "mode_w", "s.1.",

-                              "s.2.", "s.3..", "mode_s", "pairs_s", "outD_v")))

-        names(tmp) <- c("n_genotypes", "n_peaks", "n_sinks", "gamma", "gamma_star",

+                              "rsign"))) ## , "w.1.", "w.2.", "w.3..", "mode_w", "s.1.",

+        ## "s.2.", "s.3..", "mode_s", "pairs_s", "outD_v")))

+        if(length(tmp) >= 24) ## the new version

+            stopifnot(identical(names(tmp)[c(20, 24)],

+                                c("steps_m", "mProbOpt_0")))

+        ## steps_m: the mean number of steps over the entire landscape to

+        ## reach the global optimum

+        ## mProbOpt_0: The mean probability to

+        ## reach that optimum (again averaged over the entire landscape).

+        ## but if there are multiple optima, there are many of these

+        names(tmp)[1:13] <- c("n_genotypes", "n_peaks", "n_sinks", "gamma", "gamma_star",

                         "r/s","n_chains", "n_steps", "n_origins", "max_depth",

-                        "epist_magn", "epist_sign", "epist_rsign",

-                        "walsh_coef_1", "walsh_coef_2", "walsh_coef_3", "mode_walsh",

-                        "synerg_coef_1", "synerg_coef_2", "synerg_coef_3", "mode_synerg",

-                        "std_dev_pairs", "var_outdegree")

+                        "epist_magn", "epist_sign", "epist_rsign")## ,

+                        ## "walsh_coef_1", "walsh_coef_2", "walsh_coef_3", "mode_walsh",

+                        ## "synerg_coef_1", "synerg_coef_2", "synerg_coef_3", "mode_synerg",

+        ## "std_dev_pairs", "var_outdegree")

     } else {

         tmp <- readLines(fnret)

     }

@@ -1866,14 +1866,24 @@ sampledGenotypes <- function(y, genes = NULL) {

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

         y <- y[, cols]

     }

-    nn <- colnames(y)

+    ## nn <- colnames(y)

+    genotlabel <- function(u, nn = colnames(y)) {

+        if(any(is.na(u))) return(NA)

+        else {

+            return(paste(sort(nn[as.logical(u)]), collapse = ", "))

+        }

+    }

+    ## df <- data.frame(table(

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

+    ## ))

     df <- data.frame(table(

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

-    ))

+        apply(y, 1, genotlabel),

+        useNA = "ifany"))

+

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

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

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

-    attributes(df)$ShannonI <- shannonI(df$Freq)

+    attributes(df)$ShannonI <- shannonI(na.omit(df)$Freq)

     class(df) <- c("sampledGenotypes", class(df))

     return(df)

 }

@@ -1,3 +1,14 @@

+Changes in version 2.7.2 (2017-09-27):

+	- genot_to_adj_mat in C++.

+	- fast_peaks (for no backmutation cases).

+	- Better explanation and testing of peaks and valleys.

+	- Clarified simOGraph transitive reduction.

+	- Better handling of ti corner cases.

+	- Magellan reading fuctions adapted to output of newer (as of

+	  2017-07) version of Magellan.

+	- sorting gene names in allGenotypes_to_matrix.

+	- sampledGenotypes: genotype names with sorted gene names.

+

 Changes in version 2.6.0 (for BioC 3.5):

 	- Many additions to the vignette and documentation.

 	- LOD and POM (lines of descent, path of maximum, sensu Szendro et

@@ -34,8 +34,9 @@ s = 0.1, sh = -0.1, typeDep = "AND")

   \item{removeDirectIndirect}{

     Ensure that no two nodes are connected both directly (i.e., with an

   edge between them) and indirectly, through intermediate nodes. If

-  TRUE, the direct connections are removed from the graph starting from

-  the bottom.

+  TRUE, we return the transitive reduction of the DAG.

+  %% the final DAG returned is the transitive reduction of thethe direct

+  %% connections are removed from the graph starting from the bottom.

 }

 \item{rootName}{

   The name you want to give the "Root" node.

@@ -755,7 +755,7 @@ void addToPhylog(PhylogName& phylog,

-static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

+static void nr_innerBNB (const fitnessEffectsAll& fitnessEffects,

 			const double& initSize,

 			const double& K,

 			// const double& alpha,

@@ -809,7 +809,9 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 			const double& PDBaseline,

 			const double& checkSizePEvery,

 			const bool& AND_DrvProbExit,

-			const std::vector< std::vector<int> >& fixation_l) {

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

+			 int& ti_dbl_min,

+			 int& ti_e3) {

   double nextCheckSizeP = checkSizePEvery;

   const int numGenes = fitnessEffects.genomeSize;

@@ -939,8 +941,10 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

   std::multimap<double, int> mapTimes;

   //std::multimap<double, int>::iterator m1pos;

-  int ti_dbl_min = 0;

-  int ti_e3 = 0;

+  // int ti_dbl_min = 0;

+  // int ti_e3 = 0;

+  ti_dbl_min = 0;

+  ti_e3 = 0;

       // // FIXME debug

@@ -1239,6 +1243,10 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	mapTimes_updateP(mapTimes, popParams, u_1, tmpdouble1);

 	popParams[u_1].timeLastUpdate = currentTime;

+#ifdef DEBUGV

+	detect_ti_duplicates(mapTimes, tmpdouble1, u_1);  

+#endif

+	

 #ifdef DEBUGV

 	Rcpp::Rcout << "\n\n     ********* 5.2: call to ti_nextTime, update one ******\n For to_update = \n " 

 		  << "     tSample  = " << tSample

@@ -1269,6 +1277,13 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	popParams[u_1].timeLastUpdate = currentTime;

 	popParams[u_2].timeLastUpdate = currentTime;

+#ifdef DEBUGV

+	detect_ti_duplicates(mapTimes, tmpdouble1, u_1);

+	detect_ti_duplicates(mapTimes, tmpdouble2, u_2);

+#endif	

+

+	

+

 #ifdef DEBUGV

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

 		  << "     tSample  = " << tSample

@@ -1306,6 +1321,9 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 					     tSample, ti_dbl_min, ti_e3);

 	  mapTimes_updateP(mapTimes, popParams, i, tmpdouble1);

 	  popParams[i].timeLastUpdate = currentTime;

+#ifdef DEBUGV

+	  detect_ti_duplicates(mapTimes, tmpdouble1, i);

+#endif

 #ifdef DEBUGV

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

@@ -1335,6 +1353,7 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

     // ******************** 5.3 and do we sample? *********** 

     // Find minimum to know if we need to sample the whole pop

+    // We also obtain the nextMutant

     getMinNextMutationTime4(nextMutant, minNextMutationTime, 

 			    mapTimes);

@@ -1438,7 +1457,31 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	    Rcpp::Rcout <<"\n     Creating new species   " << (numSpecies - 1)

 			<< "         from species "  <<   nextMutant;

 	  }

-	

+	  

+#ifdef DEBUGW

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

+	    DP2(currentTime); //this is set to minNextMutationTime above

+	    DP2(minNextMutationTime);

+	    DP2(tSample);

+	    DP2(popParams[nextMutant].timeLastUpdate);

+	    DP2( (currentTime -  popParams[nextMutant].timeLastUpdate) );

+	    DP2( (currentTime <  popParams[nextMutant].timeLastUpdate) );

+	    DP2( (currentTime ==  popParams[nextMutant].timeLastUpdate) );

+	    DP2(nextMutant);

+	    DP2(u_1);

+	    DP2(u_2);

+	    DP2(tmpdouble1);

+	    DP2(tmpdouble2);

+	    DP2(popParams[nextMutant].timeLastUpdate);

+	    DP2(popParams[u_1].timeLastUpdate);

+	    DP2(popParams[u_2].timeLastUpdate);

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[u_2].timeLastUpdate) );

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[nextMutant].timeLastUpdate) );

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[0].timeLastUpdate) );

+	    print_spP(popParams[nextMutant]);

+	    throw std::out_of_range("new species: currentTime - timeLastUpdate[sp] out of range. ***###!!!Serious bug!!!###***");

+	  }

+#endif	  

 	  tmpParam.popSize = 1;

 	  nr_fitness(tmpParam, popParams[nextMutant],

@@ -1523,12 +1566,75 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	  to_update = 2; 

 #ifdef DEBUGW

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

-	    DP2(currentTime);

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

+	    // Yes, the difference could be 0 if two next mutation times are identical.

+	    // You enable detect_ti_duplicates and use trigger-duplicated-ti.R

+	    // to see it.

+	    // Often the involved culprits (nextMutant and the other, say sp)

+	    // were lastUpdated with tiny difference and they were, when updated

+	    // given an identical ti, each in its own run.

+	    // Key is not timeLastUpdate. This is a possible sequence of events:

+	    //    - at time t0, species that will become nextMutant is updated and gets ti = tinm

+	    //    - t1: species u1 gets ti = tinm

+	    //    - t2: species u2 gets some ti > tinm

+	    //    - tinm becomes minimal, so we mutate u1, and it mutates to u2

+	    //    - (so now the timeLastUpdate of u1 = u2 = tinm)

+	    //    - nextMutant is now mutated, and it mutates to u2, which becomes sp

+	    //    - tinm = timeLastUpdate of u1 and u2.

+	    //    - You will also see that number of mutations, or genotypes are such

+	    //      that, in this case, u2 is the most mutated, etc.

+	    //    - If you enable the detect_ti_duplicates, you would have seen duplicated ti

+	    //      for nextMutant and u1

+

+	    //   Even simpler is if above, nextMutant will mutate to u1 (not u2) so u1 becomes sp.

+	    DP2(currentTime); //this is set to minNextMutationTime above

+	    DP2(minNextMutationTime);

+	    DP2(tSample);

+	    DP2(popParams[sp].timeLastUpdate);

+	    DP2( (currentTime -  popParams[sp].timeLastUpdate) );

+	    DP2( (currentTime <  popParams[sp].timeLastUpdate) );

+	    DP2( (currentTime ==  popParams[sp].timeLastUpdate) );

 	    DP2(sp);

+	    DP2(nextMutant);

+	    DP2(u_1);

+	    DP2(u_2);

+	    DP2(tmpdouble1);

+	    DP2(tmpdouble2);

 	    DP2(popParams[sp].timeLastUpdate);

+	    DP2(popParams[nextMutant].timeLastUpdate);

+	    DP2(popParams[u_1].timeLastUpdate);

+	    DP2(popParams[u_2].timeLastUpdate);

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[u_2].timeLastUpdate) );

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[nextMutant].timeLastUpdate) );

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[0].timeLastUpdate) );

 	    print_spP(popParams[sp]);

-	    throw std::out_of_range("currentTime - timeLastUpdate out of range. Serious bug!");

+	    print_spP(popParams[nextMutant]);

+	    throw std::out_of_range("currentTime - timeLastUpdate[sp] out of range.  ***###!!!Serious bug!!!###***");

+	  }

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

+	    DP2(currentTime); //this is set to minNextMutationTime above

+	    DP2(minNextMutationTime);

+	    DP2(tSample);

+	    DP2(popParams[nextMutant].timeLastUpdate);

+	    DP2( (currentTime -  popParams[nextMutant].timeLastUpdate) );

+	    DP2( (currentTime <  popParams[nextMutant].timeLastUpdate) );

+	    DP2( (currentTime ==  popParams[nextMutant].timeLastUpdate) );

+	    DP2(sp);

+	    DP2(nextMutant);

+	    DP2(u_1);

+	    DP2(u_2);

+	    DP2(tmpdouble1);

+	    DP2(tmpdouble2);

+	    DP2(popParams[sp].timeLastUpdate);

+	    DP2(popParams[nextMutant].timeLastUpdate);

+	    DP2(popParams[u_1].timeLastUpdate);

+	    DP2(popParams[u_2].timeLastUpdate);

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[u_2].timeLastUpdate) );

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[nextMutant].timeLastUpdate) );

+	    DP2( (popParams[u_1].timeLastUpdate - popParams[0].timeLastUpdate) );

+	    print_spP(popParams[sp]);

+	    print_spP(popParams[nextMutant]);

+	    throw std::out_of_range("currentTime - timeLastUpdate[nextMutant] out of range. ***###!!!Serious bug!!!###***");

 	  }

 #endif

 	  // if(verbosity >= 2) {

@@ -1888,11 +1994,17 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

   // 5.1 Initialize 

   int numRuns = 0;

+  int numRecoverExcept = 0;

   bool forceRerun = false;

   double currentTime = 0;

   int iter = 0;

+  int ti_dbl_min = 0;

+  int ti_e3 = 0;

+

+  int accum_ti_dbl_min = 0;

+  int accum_ti_e3 = 0;

   // bool AND_DrvProbExit = ( (cpDetect >= 0) &&

   // 			     (detectionDrivers < 1e9) &&

   // 			     (detectionSize < std::numeric_limits<double>::infinity()));

@@ -1975,13 +2087,21 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

 		  PDBaseline,

 		  checkSizePEvery,

 		  AND_DrvProbExit,

-		  fixation_l);

+		  fixation_l,

+		  ti_dbl_min,

+		  ti_e3);

       ++numRuns;

       forceRerun = false;

+      accum_ti_dbl_min += ti_dbl_min;

+      accum_ti_e3 += ti_e3;

     } catch (rerunExcept &e) {

       Rcpp::Rcout << "\n Recoverable exception " << e.what() 

 		  << ". Rerunning.";

       forceRerun = true;

+      ++numRecoverExcept;

+      ++numRuns; // exception should count here!

+      accum_ti_dbl_min += ti_dbl_min;

+      accum_ti_e3 += ti_e3;

     } catch (const std::exception &e) {

       Rcpp::Rcout << "\n Unrecoverable exception: " << e.what()

 		  << ". Aborting. \n";

@@ -2125,9 +2245,11 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

 										      Named("child") = phylog.child,

 										      Named("time") = phylog.time

 										      ),

-					       Named("UnrecoverExcept") = false)

+					       Named("UnrecoverExcept") = false,

+					       Named("numRecoverExcept") = numRecoverExcept,

+					       Named("accum_ti_dbl_min") = accum_ti_dbl_min,

+					       Named("accum_ti_e3") = accum_ti_e3)

 		 );

-

 }

@@ -12,6 +12,10 @@ SEXP OncoSimulR_evalRGenotypeAndMut(SEXP rGSEXP, SEXP rFESEXP, SEXP muEFSEXP, SE

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

 SEXP OncoSimulR_accessibleGenotypes(SEXP ySEXP, SEXP xSEXP, SEXP numMutSEXP, SEXP thSEXP);

+SEXP OncoSimulR_genot2AdjMat(SEXP ySEXP, SEXP xSEXP, SEXP numMutSEXP);

+SEXP OncoSimulR_peaksLandscape(SEXP ySEXP, SEXP xSEXP, SEXP numMutSEXP, SEXP thSEXP);

+SEXP OncoSimulR_accessibleGenotypes_former(SEXP ySEXP, SEXP xSEXP, SEXP numMutSEXP, SEXP thSEXP);

+

 // The number is the number of arguments

 R_CallMethodDef callMethods[]  = {

   {"OncoSimulR_nr_BNB_Algo5", (DL_FUNC) &OncoSimulR_nr_BNB_Algo5, 37},

@@ -19,6 +23,9 @@ R_CallMethodDef callMethods[]  = {

   {"OncoSimulR_evalRGenotype", (DL_FUNC) &OncoSimulR_evalRGenotype, 5},

   {"OncoSimulR_evalRGenotypeAndMut", (DL_FUNC) &OncoSimulR_evalRGenotypeAndMut, 6},

   {"OncoSimulR_accessibleGenotypes", (DL_FUNC) &OncoSimulR_accessibleGenotypes, 4},

+  {"OncoSimulR_genot2AdjMat", (DL_FUNC) &OncoSimulR_genot2AdjMat, 3},

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

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

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

   {NULL, NULL, 0}

 };

@@ -141,6 +141,50 @@ BEGIN_RCPP

  END_RCPP

 }

+// genotype fitness matrix to adjacency matrix of genotypes

+Rcpp::NumericMatrix genot2AdjMat(Rcpp::IntegerMatrix y, Rcpp::NumericVector f, Rcpp::IntegerVector numMut);

+RcppExport SEXP OncoSimulR_genot2AdjMat(SEXP ySEXP, SEXP fSEXP, SEXP numMutSEXP) {

+BEGIN_RCPP

+    Rcpp::RObject __result;

+// Rcpp::RNGScope __rngScope;

+ Rcpp::traits::input_parameter< Rcpp::IntegerMatrix >::type y(ySEXP);

+ Rcpp::traits::input_parameter< Rcpp::NumericVector >::type f(fSEXP);

+ Rcpp::traits::input_parameter< Rcpp::IntegerVector >::type numMut(numMutSEXP);

+ __result = Rcpp::wrap(genot2AdjMat(y, f, numMut));

+ return __result;

+ END_RCPP

+}

+

+

+// evalRGenotypeAndMut

+Rcpp::IntegerVector peaksLandscape(Rcpp::IntegerMatrix y, Rcpp::NumericVector f, Rcpp::IntegerVector numMut, double th);

+RcppExport SEXP OncoSimulR_peaksLandscape(SEXP ySEXP, SEXP fSEXP, SEXP numMutSEXP, SEXP thSEXP) {

+BEGIN_RCPP

+    Rcpp::RObject __result;

+// Rcpp::RNGScope __rngScope;

+ Rcpp::traits::input_parameter< Rcpp::IntegerMatrix >::type y(ySEXP);

+ Rcpp::traits::input_parameter< Rcpp::NumericVector >::type f(fSEXP);

+ Rcpp::traits::input_parameter< Rcpp::IntegerVector >::type numMut(numMutSEXP);

+ Rcpp::traits::input_parameter< double >::type th(thSEXP);

+ __result = Rcpp::wrap(peaksLandscape(y, f, numMut, th));

+ return __result;

+ END_RCPP

+}

+

+// just for testing. Eventually remove

+Rcpp::IntegerVector accessibleGenotypes_former(Rcpp::IntegerMatrix y, Rcpp::NumericVector f, Rcpp::IntegerVector numMut, double th);

+RcppExport SEXP OncoSimulR_accessibleGenotypes_former(SEXP ySEXP, SEXP fSEXP, SEXP numMutSEXP, SEXP thSEXP) {

+BEGIN_RCPP

+    Rcpp::RObject __result;

+// Rcpp::RNGScope __rngScope;

+ Rcpp::traits::input_parameter< Rcpp::IntegerMatrix >::type y(ySEXP);

+ Rcpp::traits::input_parameter< Rcpp::NumericVector >::type f(fSEXP);

+ Rcpp::traits::input_parameter< Rcpp::IntegerVector >::type numMut(numMutSEXP);

+ Rcpp::traits::input_parameter< double >::type th(thSEXP);

+ __result = Rcpp::wrap(accessibleGenotypes_former(y, f, numMut, th));

+ return __result;

+ END_RCPP

+}

 // // readFitnessEffects

@@ -22,12 +22,12 @@ inline int HammingDistance(const Rcpp::IntegerVector& x, const Rcpp::IntegerVect

 }

-

+// eventually remove this. Left now for testing

 // [[Rcpp::export]]

-Rcpp::IntegerVector accessibleGenotypes(Rcpp::IntegerMatrix y,

-					Rcpp::NumericVector f,

-					Rcpp::IntegerVector numMut, //

-					double th) {

+Rcpp::IntegerVector accessibleGenotypes_former(Rcpp::IntegerMatrix y,

+					       Rcpp::NumericVector f,

+					       Rcpp::IntegerVector numMut, //

+					       double th) {

   // Return just the indices. Could preserve fitness, but would need

   // another matrix.

   int ng = y.nrow(); //it counts the wt

@@ -106,6 +106,189 @@ Rcpp::IntegerVector accessibleGenotypes(Rcpp::IntegerMatrix y,

+

+

+// [[Rcpp::export]]

+Rcpp::NumericMatrix genot2AdjMat(Rcpp::IntegerMatrix y,

+				 Rcpp::NumericVector f,

+				 Rcpp::IntegerVector numMut) {

+  // Return just the indices. Could preserve fitness, but would need

+  // another matrix.

+  int ng = y.nrow(); //it counts the wt

+  Rcpp::NumericMatrix adm(ng, ng);

+  

+  // fill with NAs: https://stackoverflow.com/a/23753626

+  // Filling with NAs and in general having NAs might lead to performance

+  // penalties. But I use the NAs in a lot of the code for accessible

+  // genotypes, etc.

+  std::fill( adm.begin(), adm.end(), Rcpp::NumericVector::get_na() ) ;

+  int numMutdiff = 0;

+  // I would have thought this would be faster. It ain't.

+  // The last genotype never accesses anything.

+  // for(int i = 0; i < (ng - 1); ++i) {

+  //   // Candidate genotypes to be accessed from i are always of larger

+  //   // mutation by 1. And candidates can thus not have smaller index

+  //   for(int j = (i + 1); j < ng; ++j) {

+  //     if( (numMut(j) == (numMut(i) + 1)) &&

+  // 	  ( (f(j) - f(i)) >= th) &&

+  // 	  (HammingDistance(y(i, _), y(j, _)) == 1) ) {

+  // 	adm(i, j) = 1;

+  //     } else if( (numMut(j) > (numMut(i) + 1)) ) {

+  // 	break;

+  //     }

+  //   }

+  // }

+

+  // The last genotype never accesses anything.

+  for(int i = 0; i < (ng - 1); ++i) {

+    // Candidate genotypes to be accessed from i are always of larger

+    // mutation by 1. And candidates can thus not have smaller index

+    for(int j = (i + 1); j < ng; ++j) {

+      numMutdiff = numMut(j) - numMut(i);

+      if( numMutdiff > 1) { // no more to search

+  	break; 

+      } else if(numMutdiff == 1) {

+  	if( HammingDistance(y(i, Rcpp::_), y(j, Rcpp::_)) == 1) {

+  	  adm(i, j) =  (f(j) - f(i));

+  	}

+      }

+    }

+  }

+  return adm;

+}

+

+

+Rcpp::IntegerMatrix integerAdjMat(Rcpp::IntegerMatrix y,

+				  Rcpp::NumericVector f,

+				  Rcpp::IntegerVector numMut, //

+				  double th) {

+  // Return a genotype adjacency matrix with a 1 if genotype j is

+  // accessible (fitness >, within th) from i.

+  int ng = y.nrow(); //it counts the wt

+  Rcpp::IntegerMatrix adm(ng, ng);

+  int numMutdiff = 0;

+  // I would have thought this would be faster. It ain't.

+  // The last genotype never accesses anything.

+  // for(int i = 0; i < (ng - 1); ++i) {

+  //   // Candidate genotypes to be accessed from i are always of larger

+  //   // mutation by 1. And candidates can thus not have smaller index

+  //   for(int j = (i + 1); j < ng; ++j) {

+  //     if( (numMut(j) == (numMut(i) + 1)) &&

+  // 	  ( (f(j) - f(i)) >= th) &&

+  // 	  (HammingDistance(y(i, _), y(j, _)) == 1) ) {

+  // 	adm(i, j) = 1;

+  //     } else if( (numMut(j) > (numMut(i) + 1)) ) {

+  // 	break;

+  //     }

+  //   }

+  // }

+

+  // The last genotype never accesses anything.

+  for(int i = 0; i < (ng - 1); ++i) {

+    // Candidate genotypes to be accessed from i are always of larger

+    // mutation by 1. And candidates can thus not have smaller index

+    for(int j = (i + 1); j < ng; ++j) {

+      numMutdiff = numMut(j) - numMut(i);

+      if( numMutdiff > 1) { // no more to search

+  	break; 

+      } else if(numMutdiff == 1) {

+  	// f(j) - f(i) is faster than HammingDistance

+  	// but might lead to more evals?

+  	// or fewer, depending on landscape

+  	if( ( (f(j) - f(i)) >= th) &&

+	    (HammingDistance(y(i, Rcpp::_), y(j, Rcpp::_)) == 1)

+  	    ) {

+  	  adm(i, j) = 1;

+	  // Rcpp::Rcout << "i = " << i << " j = " << j << " adm " << adm(i,j) << "\n"; 

+  	}

+      }

+    }

+  }

+  return adm;

+}

+

+

+// used in both peaks and accessible genotypes

+Rcpp::IntegerVector accessibleGenotypesPeaksLandscape(Rcpp::IntegerMatrix y,

+						      Rcpp::NumericVector f,

+						      Rcpp::IntegerVector numMut, //

+						      double th,

+						      bool returnpeaks) {

+  // Return the indices. This is like accessibleGenotypes, but we do an

+  // extra loop

+  int ng = y.nrow(); //it counts the wt

+  Rcpp::IntegerMatrix adm(ng, ng);

+

+  adm = integerAdjMat(y, f, numMut, th);

+  

+  int numMutdiff = 0;

+

+  // Slightly different logic from R: Do not resize object; set the row to

+  // 0.

+  int colsum = 0;

+  // int indicator = 0; // indicator != 0 means we set one row to 0

+  // so we need to iterate at least once more.

+  

+  // accessible is the genotype number, not the column!  WT is 1,

+  // etc. This makes it easy to keep track of which are accessible.

+  Rcpp::IntegerVector accessible = Rcpp::seq_len(ng);