@@ -71,18 +71,21 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

     x_from <- y_from <- x_to <- y_to <- Change <- muts <-

-        label <- fitness <- Type <- NULL

-                

+        label <- birth <- Type <- NULL

+    

+    if ("Fitness" %in% colnames(tfm$afe)) {

+      colnames(tfm$afe)[which(colnames(tfm$afe) == "Fitness")] <- "Birth"

+    }

     dd <- data.frame(muts = mutated,

-                     fitness = tfm$afe$Fitness,

+                     birth = tfm$afe$Birth,

                      label = tfm$afe$Genotype,

                      stringsAsFactors = TRUE)

     cl <- gaj[vv]

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

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

+                     y_from = tfm$afe$Birth[vv[, 1]],

                      x_to = mutated[vv[, 2]],

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

+                     y_to = tfm$afe$Birth[vv[, 2]],

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

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

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

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

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

     p <- ggplot() +

-        xlab("") + ylab("Fitness") + 

+        xlab("") + ylab("Birth") + 

         theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),

               panel.grid.minor.x = element_blank()) +

         geom_segment(data = sg,

@@ -121,9 +124,9 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

     if(show_labels && use_ggrepel) {

         p <- p + geom_text_repel(data = dd[-c(minF, maxF), ],

-                                 aes(x = muts, y = fitness, label = label)) +

+                                 aes(x = muts, y = birth, label = label)) +

             geom_label_repel(data = ddMM,

-                             aes(x = muts, y = fitness, label = label, fill = Type),

+                             aes(x = muts, y = birth, label = label, fill = Type),

                              color = "black",

                              fontface = "bold")

             ## geom_label_repel(data = dd[maxF, ], aes(x = muts, y = fitness, label = label),

@@ -141,10 +144,10 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

         }

         p <- p + geom_text(data = dd[-c(minF, maxF), ],

-                           aes(x = muts, y = fitness, label = label),

+                           aes(x = muts, y = birth, label = label),

                            vjust = -0.2, hjust = "inward") +

             geom_label(data = ddMM,

-                       aes(x = muts, y = fitness, label = label, fill = Type),

+                       aes(x = muts, y = birth, label = label, fill = Type),

                        vjust = -0.2, hjust = "inward", color = "black",

                        fontface = "bold")

             ## geom_label(data = dd[maxF, ], aes(x = muts, y = fitness, label = label),

@@ -230,6 +230,9 @@ fast_peaks <- function(x, th) {

                           th)]))

 }

+## FIXME! Unclear if we need the complete fitness matrix or if it will work with

+## pieces (say, after removing those with fitness < WT, to try to make it faster

+## and for general use in other cases)

 ## wrapper to the C++ code

 wrap_accessibleGenotypes <- function(x, th) {

@@ -319,65 +319,6 @@ faster_accessible_genotypes_R <- function(x, th) {

 }

-## ## This uses slam, but that is actually slower because

-## ## of the assignment

-## faster_accessible_genots_slam <- function(x, th = 0) {

-

-##     ## Given a genotype matrix, return the genotypes that are accessible

-##     ## via creating a directed adjacency matrix between genotypes

-##     ## connected (i.e., those that differ by gaining one mutation). 0

-##     ## means not connected, 1 means connected.

-    

-##     ## There is a more general function in OncoSimulR that will give the

-##     ## fitness difference. But not doing the difference is faster than

-##     ## just setting a value, say 1, if all we want is to keep track of

-##     ## accessible ones. And by using only 0/1 we can store only an

-##     ## integer. And no na.omits, etc. Is too restricted? Yes. But for

-##     ## simulations and computing just accessible genotypes, probably a

-##     ## hell of a lot faster.

-

-##     ## Well, this is not incredibly fast either.

-    

-##     ## Make sure sorted, so ancestors always before descendants

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

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

-##     rm(rs0)

-    

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

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

-##     rs <- rowSums(y)

-

-##     ## If 0, not accessible

-##     adm <- slam::simple_triplet_zero_matrix(nrow = length(rs), ncol = length(rs),

-##                                       mode = "integer")

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

-##         candidates <- which(rs == (rs[i] + 1))

-##         for(j in candidates) {

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

-##             ## if(sumdiff == 1)

-##             if( (sum(abs(y[j, ] - y[i, ])) == 1) &&

-##                 ( (f[j] - f[i]) >= th ) )

-##                 adm[i, j] <- 1L

-##         }

-##     }

-

-##     colnames(adm) <- rownames(adm) <- 1:ncol(adm)

-##     admtmp <- adm[, -1, drop = FALSE] ## we do not want the root column.

-##     while(TRUE) {

-##         ## We remove inaccessible cols (genotypes) and the corresponding

-##         ## rows repeatedly until nothing left to remove; any column left

-##         ## is therefore accessible throw at least one path.

-

-##         ## inacc_col <- which(slam::colapply_simple_triplet_matrix(admtmp, FUN = sum) == 0L)

-##         inacc_col <- which(slam::col_sums(admtmp) == 0L)

-##         if(length(inacc_col) == 0) break;

-##         inacc_row <- inacc_col + 1 ## recall root row is left

-##         admtmp <- admtmp[-inacc_row, -inacc_col, drop = FALSE]

-##     }

-##     return(as.numeric(c(colnames(adm)[1], colnames(admtmp))))

-## }

-

-

 generate_matrix_genotypes <- function(g) {

     ## FIXME future: do this for order too? Only if rfitness for order.

     ## Given a number of genes, generate all possible genotypes.

@@ -470,25 +411,6 @@ genot_to_adj_mat <- function(x) {

 }

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

@@ -16,12 +16,6 @@

 ## plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-

 ##     plot.genotype_fitness_matrix <- plotFitnessLandscape 

-## FIXME: show only accessible paths? 

-## FIXME: when show_labels = FALSE we still show the boxes

-##        and some of the labels.!!!

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

@@ -31,7 +25,12 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

                                  ...) {

     ## FIXME future:

-

+    ## FIXME: when show_labels = FALSE we still show the boxes

+    ##        and some of the labels.!!!

+    ## FIXME: if using only_accessible, maybe we

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

+    ## approach as fallback (if identical fitness)

+    

     ## - Allow passing order effects. Change "allGenotypes_to_matrix"

     ##   below? Probably not, as we cannot put order effects as a

     ##   matrix. Do something else, like allow only order effects if from

@@ -70,14 +69,7 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

     }

     vv <- which(!is.na(gaj), arr.ind = TRUE)

-    ## plot(x = mutated, y = e1$Fitness, ylab = "Fitness",

-    ##      xlab = "", type = "n", axes = FALSE)

-    ## box()

-    ## axis(2)

-    ## text(x = mutated, y = x$Fitness, labels = x$Genotype)

-

-    ## The R CMD CHEKC notes about no visible binding for global variable

-

+  

     x_from <- y_from <- x_to <- y_to <- Change <- muts <-

         label <- fitness <- Type <- NULL

@@ -184,27 +176,6 @@ plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-

 ######################################################################

-## wrap_filter_inaccessible <- function(x, max_num_genotypes, accessible_th) {

-##     ## wrap it, for my consumption

-##     tfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)

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

-##     gaj <- genot_to_adj_mat(tfm$gfm)

-##     gaj <- filter_inaccessible(gaj, accessible_th)

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

-##     mutated <- mutated[remaining]

-##     tfm$afe <- tfm$afe[remaining, , drop = FALSE]

-##     return(list(remaining = remaining,

-##                 mutated = mutated,

-##                 tfm = tfm))

-## }

-

-## No longer being used. Used to be in rfitness

-## count_accessible_g <- function(gfm, accessible_th) {

-##     gaj <- genot_to_adj_mat(gfm)

-##     gaj <- filter_inaccessible(gaj, accessible_th)

-##     return(ncol(gaj) - 1)

-## }

-

 ## There is now C++ code to get just the locations/positions of the

 ## accessible genotypes

@@ -1,4 +1,4 @@

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

+## Copyright 2013-2021 Ramon Diaz-Uriarte

 ## This program is free software: you can redistribute it and/or modify

 ## it under the terms of the GNU General Public License as published by

@@ -17,7 +17,8 @@

 ##     plot.genotype_fitness_matrix <- plotFitnessLandscape 

 ## FIXME: show only accessible paths? 

-

+## FIXME: when show_labels = FALSE we still show the boxes

+##        and some of the labels.!!!

 ## FIXME: if using only_accessible, maybe we

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

 ## approach as fallback (if identical fitness)

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

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

@@ -223,9 +223,6 @@ filter_inaccessible <- function(x, accessible_th) {

     return(x)

 }

-f1 <- function() {}

-

-x <- 99

 ## wrapper to the C++ code

 wrap_accessibleGenotypes <- function(x, th) {

@@ -192,16 +192,23 @@ plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-

 ##                 tfm = tfm))

 ## }

-count_accessible_g <- function(gfm, accessible_th) {

-    gaj <- genot_to_adj_mat(gfm)

-    gaj <- filter_inaccessible(gaj, accessible_th)

-    return(ncol(gaj) - 1)

-}

+## No longer being used. Used to be in rfitness

+## count_accessible_g <- function(gfm, accessible_th) {

+##     gaj <- genot_to_adj_mat(gfm)

+##     gaj <- filter_inaccessible(gaj, accessible_th)

+##     return(ncol(gaj) - 1)

+## }

+

+## There is now C++ code to get just the locations/positions of the

+## accessible genotypes

 filter_inaccessible <- function(x, accessible_th) {

     ## Return an adjacency matrix with only accessible paths. x is the gaj

     ## matrix created in the plots. The difference between genotypes

     ## separated by a hamming distance of 1

+

+    ## FIXME: could do the x[, -1] before loop and not add the 1

+    ## inside while, and do that at end

     colnames(x) <- rownames(x) <- 1:ncol(x)

     while(TRUE) {

         ## remove first column

@@ -216,6 +223,132 @@ filter_inaccessible <- function(x, accessible_th) {

     return(x)

 }

+f1 <- function() {}

+

+x <- 99

+

+## wrapper to the C++ code

+wrap_accessibleGenotypes <- function(x, th) {

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

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

+    o_numMut <- order(numMut)

+    x <- x[o_numMut, ]

+    numMut <- numMut[o_numMut]

+    

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

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

+

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

+                               as.integer(numMut),

+                               th)

+    return(acc[acc > 0])

+}

+

+## A transitional function

+faster_accessible_genotypes_R <- function(x, th) {

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

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

+    rm(rs0)

+    

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

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

+    rs <- rowSums(y)

+

+   ## If 0, not accessible

+   ## adm <- slam::simple_triplet_zero_matrix(nrow = length(rs), ncol = length(rs),

+   ##                                          mode = "integer")

+   

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

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

+   

+   ## Most time is gone here

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

+        candidates <- which(rs == (rs[i] + 1))

+        for(j in candidates) {

+            if( (sum(abs(y[j, ] - y[i, ])) == 1) &&

+                ( (f[j] - f[i]) >= th ) ) {

+                ## actually, this is the largest time sink using slam

+                adm[i, j] <- 1L

+                }

+        }

+    }

+

+    colnames(adm) <- rownames(adm) <- 1:ncol(adm)

+    admtmp <- adm[, -1, drop = FALSE] ## we do not want the root column.

+    while(TRUE) {

+        ## We remove inaccessible cols (genotypes) and the corresponding

+        ## rows repeatedly until nothing left to remove; any column left

+        ## is therefore accessible throw at least one path.

+

+        ## inacc_col <- which(slam::colapply_simple_triplet_matrix(admtmp, FUN = sum) == 0L)

+        inacc_col <- which(colSums(admtmp) == 0L)

+        if(length(inacc_col) == 0) break;

+        inacc_row <- inacc_col + 1 ## recall root row is left

+        admtmp <- admtmp[-inacc_row, -inacc_col, drop = FALSE]

+    }

+    return(as.numeric(c(colnames(adm)[1], colnames(admtmp))))

+

+}

+

+

+## ## This uses slam, but that is actually slower because

+## ## of the assignment

+## faster_accessible_genots_slam <- function(x, th = 0) {

+

+##     ## Given a genotype matrix, return the genotypes that are accessible

+##     ## via creating a directed adjacency matrix between genotypes

+##     ## connected (i.e., those that differ by gaining one mutation). 0

+##     ## means not connected, 1 means connected.

+    

+##     ## There is a more general function in OncoSimulR that will give the

+##     ## fitness difference. But not doing the difference is faster than

+##     ## just setting a value, say 1, if all we want is to keep track of

+##     ## accessible ones. And by using only 0/1 we can store only an

+##     ## integer. And no na.omits, etc. Is too restricted? Yes. But for

+##     ## simulations and computing just accessible genotypes, probably a

+##     ## hell of a lot faster.

+

+##     ## Well, this is not incredibly fast either.

+    

+##     ## Make sure sorted, so ancestors always before descendants

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

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

+##     rm(rs0)

+    

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

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

+##     rs <- rowSums(y)

+

+##     ## If 0, not accessible

+##     adm <- slam::simple_triplet_zero_matrix(nrow = length(rs), ncol = length(rs),

+##                                       mode = "integer")

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

+##         candidates <- which(rs == (rs[i] + 1))

+##         for(j in candidates) {

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

+##             ## if(sumdiff == 1)

+##             if( (sum(abs(y[j, ] - y[i, ])) == 1) &&

+##                 ( (f[j] - f[i]) >= th ) )

+##                 adm[i, j] <- 1L

+##         }

+##     }

+

+##     colnames(adm) <- rownames(adm) <- 1:ncol(adm)

+##     admtmp <- adm[, -1, drop = FALSE] ## we do not want the root column.

+##     while(TRUE) {

+##         ## We remove inaccessible cols (genotypes) and the corresponding

+##         ## rows repeatedly until nothing left to remove; any column left

+##         ## is therefore accessible throw at least one path.

+

+##         ## inacc_col <- which(slam::colapply_simple_triplet_matrix(admtmp, FUN = sum) == 0L)

+##         inacc_col <- which(slam::col_sums(admtmp) == 0L)

+##         if(length(inacc_col) == 0) break;

+##         inacc_row <- inacc_col + 1 ## recall root row is left

+##         admtmp <- admtmp[-inacc_row, -inacc_col, drop = FALSE]

+##     }

+##     return(as.numeric(c(colnames(adm)[1], colnames(admtmp))))

+## }

+

 generate_matrix_genotypes <- function(g) {

     ## FIXME future: do this for order too? Only if rfitness for order.

@@ -257,6 +390,8 @@ genot_to_adj_mat <- function(x) {

     ## that differ by gaining one mutation) with value being the

     ## difference in fitness between destination and origin

+    ## 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), ]

@@ -22,6 +22,8 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

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

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

                                  log = FALSE, max_num_genotypes = 2000,

+                                 only_accessible = FALSE,

+                                 accessible_th = 0,

                                  ...) {

     ## FIXME future:

@@ -37,7 +39,6 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

     ## adjacency matrix of genotype i go at row i and column i.  Follow back

     ## all entries in row i, and their previous, and forward all column i.

-

     ## gfm: genotype fitness matrix

     ## afe: all fitness effects

@@ -52,11 +53,17 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

     ## get the string representation, etc. And this is for use

     ## with OncoSimul.

+  

     tfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)

-

     mutated <- rowSums(tfm$gfm[, -ncol(tfm$gfm)])

     gaj <- genot_to_adj_mat(tfm$gfm)

+    if(only_accessible) {

+        gaj <- filter_inaccessible(gaj, accessible_th)

+        remaining <- as.numeric(colnames(gaj))

+        mutated <- mutated[remaining]

+        tfm$afe <- tfm$afe[remaining, , drop = FALSE]

+    }

     vv <- which(!is.na(gaj), arr.ind = TRUE)

     ## plot(x = mutated, y = e1$Fitness, ylab = "Fitness",

@@ -171,6 +178,44 @@ plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-

 ######################################################################

+## wrap_filter_inaccessible <- function(x, max_num_genotypes, accessible_th) {

+##     ## wrap it, for my consumption

+##     tfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)

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

+##     gaj <- genot_to_adj_mat(tfm$gfm)

+##     gaj <- filter_inaccessible(gaj, accessible_th)

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

+##     mutated <- mutated[remaining]

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

+##     return(list(remaining = remaining,

+##                 mutated = mutated,

+##                 tfm = tfm))

+## }

+

+count_accessible_g <- function(gfm, accessible_th) {

+    gaj <- genot_to_adj_mat(gfm)

+    gaj <- filter_inaccessible(gaj, accessible_th)

+    return(ncol(gaj) - 1)

+}

+

+filter_inaccessible <- function(x, accessible_th) {

+    ## Return an adjacency matrix with only accessible paths. x is the gaj

+    ## matrix created in the plots. The difference between genotypes

+    ## separated by a hamming distance of 1

+    colnames(x) <- rownames(x) <- 1:ncol(x)

+    while(TRUE) {

+        ## remove first column

+        ## We use fact that all(na.omit(x) < u) is true if all are NA

+        ## so inaccessible rows are removed and thus destination columns

+        wrm <- which(apply(x[, -1, drop = FALSE], 2,

+                           function(y) {all(na.omit(y) < accessible_th)})) + 1

+        if(length(wrm) == 0) break;

+        x <- x[-wrm, -wrm, drop = FALSE]

+    }

+    x[x < 0] <- NA

+    return(x)

+}

+

 generate_matrix_genotypes <- function(g) {

     ## FIXME future: do this for order too? Only if rfitness for order.

@@ -25,8 +25,6 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

                                  ...) {

     ## FIXME future:

-    

-   

     ## - Allow passing order effects. Change "allGenotypes_to_matrix"

     ##   below? Probably not, as we cannot put order effects as a

@@ -54,107 +52,8 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

     ## get the string representation, etc. And this is for use

     ## with OncoSimul.

-

     tfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)

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

-    ##     ## could construct the afe below by appropriately pasting the

-    ##     ## columns names

-    ##     afe <- evalAllGenotypes(allFitnessEffects(

-    ##         epistasis = from_genotype_fitness(x)),

-    ##         order = FALSE, addwt = TRUE, max = max_num_genotypes)

-    ##     ## Might not be needed with the proper gfm object (so gmf <- x)

-    ##     ## but is needed if arbitrary matrices.

-    ##     gfm <- allGenotypes_to_matrix(afe) 

-    ## } else if(inherits(x, "fitnessEffects")) {

-    ##     if(!is.null(x$orderEffects) )

-    ##         stop("We cannot yet deal with order effects")

-    ##     afe <- evalAllGenotypes(x,

-    ##                             order = FALSE,

-    ##                             addwt = TRUE, max = max_num_genotypes)

-    ##     gfm <- allGenotypes_to_matrix(afe)

-    ## } else if( (inherits(x, "evalAllGenotypes")) ||

-    ##            (inherits(x, "evalAllGenotypesMut"))) {

-    ##     if(any(grepl(">", x[, 1], fixed = TRUE)))

-    ##         stop("We cannot deal with order effects yet.")

-    ##     x <- x[, c(1, 2)]

-    ##     if(x[1, "Genotype"] != "WT") {

-    ##         ## Yes, because we expect this present below

-    ##         x <- rbind(data.frame(Genotype = "WT",

-    ##                               Fitness = 1,

-    ##                               stringsAsFactors = FALSE),

-    ##                    x)

-    ##     }

-    ##     afe <- x

-    ##     ## in case we pass an evalAllgenotypesfitandmut

-    ##     gfm <- allGenotypes_to_matrix(afe)

-    ## } else if(is.data.frame(x)) {

-    ##     ## Assume a two-column data frame of genotypes as character

-    ##     ## vectors and fitness

-    ##     if(colnames(x)[2] != "Fitness")

-    ##         stop("We cannot guess what you are passing here") 

-    ##     afe <- evalAllGenotypes(allFitnessEffects(genotFitness = x),

-    ##                             order = FALSE, addwt = TRUE,

-    ##                             max = max_num_genotypes)

-    ##     gfm <- allGenotypes_to_matrix(afe)

-    ## } else {

-    ##    stop("We cannot guess what you are passing here") 

-    ## }

-   

-

-    

-    ## if(inherits(x, "genotype_fitness_matrix")) {

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

-    ##     ## could construct the afe below by appropriately pasting the

-    ##     ## columns names

-    ##     afe <- evalAllGenotypes(allFitnessEffects(

-    ##         epistasis = from_genotype_fitness(x)),

-    ##         order = FALSE, addwt = TRUE, max = 2000)

-    ##     gfm <- allGenotypes_to_matrix(afe) ## should not be needed? gfm <- x

-    ## } else if(inherits(x, "fitnessEffects")) {

-    ##     if(!is.null(x$orderEffects) )

-    ##         stop("We cannot yet deal with order effects")

-    ##     afe <- evalAllGenotypes(x,

-    ##                             order = FALSE,

-    ##                             addwt = TRUE, max = 2000)

-    ##     gfm <- allGenotypes_to_matrix(x)

-    ## } else {

-    ##     lc <- ncol(x)

-    ##     ## detect an appropriately formatted matrix

-    ##     if( (is.data.frame(x) || is.matrix(x)) ) {

-    ##         if(ncol(x) > 2) {

-    ##             ## We cannot be sure all genotypes are present

-    ##             ## but that is not our business?

-    ##             ## colnames(x)[lc] <- "Fitness"

-    ##             ## So this must be a matrix

-    ##             afe <- evalAllGenotypes(allFitnessEffects(

-    ##                 epistasis = from_genotype_fitness(x)),

-    ##                 order = FALSE, addwt = TRUE, max = 2000)

-    ##             gfm <- allGenotypes_to_matrix(afe)

-    ##         } else {

-    ##             ## This must be a data frame

-    ##             if(!is.data.frame(x))

-    ##                 stop("How are you passing a matrix here?")

-    ##             if(colnames(x)[lc] != "Fitness")

-    ##                 stop("We cannot guess what you are passing here")

-    ##             ## We are passed an allFitnessEffects output

-    ##             ## Will be simpler later as we will know immediately

-    ##             if(x[1, "Genotype"] != "WT") {

-    ##                 ## Yes, because we expect this present below

-    ##                 x <- rbind(data.frame(Genotype = "WT",