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

## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.





## Say which are drivers: populate the drv vector.
## // In R, the user says which are the drivers. If does not say anuthing,
## // the default (NULL) then drivers are all in poset, epist, restrict. The
## // user can pass a vector with the names of the genes (not modules). Allow
## // also for empty, so this is faster if not needed. And check that if we
## // use a stopping rule based on drivers that drv vectors is not empty.


## - Say that a user can use a "0" as a gene name, but that is BAD idea.
## - Modules and order effects can be kind of funny?

## Gene names can contain no spaces, commas, or ">"


check.gm <- function(gm) {
    ## Yes, Root: we want no ambiguities.

    ## Actually, that sucks. So we do not require it, but check for
    ## consistency.
    
    if(any(gm == "Root") && (gm[1] != "Root") )
        stop("If Root is in the module table, it must be the first")

    if(any(names(gm) == "Root") && (names(gm)[1] != "Root") )
        stop("If the name Root is in the module table, it must be the first")

    if( (names(gm)[1] == "Root") && (gm[1] != "Root") )
        stop("The name Root is in the module table, but is not of Root")

    if( (gm[1] == "Root") && (names(gm)[1] != "Root") )
        stop("Root is in the module table, but with a different name")

    ## if(gm[1] != "Root")
    ##     stop("First value of a module table must be Root")
    ## if(names(gm)[1] != "Root")
    ##     stop("First name of a module table must be Root")
    if(length(unique(names(gm))) != length(gm))
        stop("Number of unique module names different from length of vector")

    if(gm[1] != "Root")
        gm <- c("Root" = "Root", gm)
    return(gm)
}

gtm2 <- function(x) {
    data.frame(cbind(nice.vector.eo(x, ","), x), stringsAsFactors = TRUE)
}

## nice.vector.eo <- function(z, sep) {
##     ## with epistasis, maybe we want sorted?
##     setdiff(unlist(lapply(strsplit(z, " "),
##                                     function(u) strsplit(u, sep))), "")
## }

nice.vector.eo <- function(z, sep, rm.sign = FALSE) {
    ## with epistasis, maybe we want sorted?
    if(! rm.sign)
        return(setdiff(unlist(lapply(strsplit(z, " "),
                              function(u) strsplit(u, sep))), ""))
    else ## remove the " ", the -, and the sep
        return(setdiff(unlist(lapply(strsplit(z, " "), function(u)
            strsplit(unlist(strsplit(u, "-")), sep))), ""))
}


gm.to.geneModuleL <- function(gm, one.to.one) {
    ## the table will be sorted by gene name
    gm <- check.gm(gm)
   
    ## the named vector with the mapping into the long geneModule df
    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]
    geneMod$Gene <- as.character(geneMod$Gene)
    ## geneMod$Module <- as.character(geneMod$Module) ## already a char
    geneMod <- geneMod[c(1, order(geneMod$Gene[-1]) + 1), ] 
    geneMod$GeneNumID <- 0:(nrow(geneMod) - 1)

    ## this assumes sorted! and they need not be
    ## rlem <- rle(geneMod$Module)
    ## geneMod$ModuleNumID <- rep( 0:(length(rlem$values) - 1), rlem$lengths)
    if(one.to.one) {
        geneMod$ModuleNumID <- geneMod$GeneNumID
        if(!(identical(geneMod$Gene, geneMod$Module)))
            stop("Impossible: we are supposed to be in one-to-one for Module-Gene.")
    } else {
        ## Works, but does not give the most ordered module names. But
        ## keeps implicit order given by user.
        idm <- seq.int(length(gm) - 1)
        idm <- c("Root" = 0L, idm)
        names(idm) <- names(gm)
        ## This sorts by the names but is not optimal either
        ## idm1 <- seq.int(length(gm) - 1)
        ## idm <- c(0L, idm1)
        ## names(idm) <- c("Root", sort(setdiff(names(gm), "Root")))
        geneMod$ModuleNumID <- idm[geneMod[, "Module"]]
    }
    if(length(unique(geneMod$Gene)) != nrow(geneMod)) {
        stop("Are there identical gene names in different modules?")
    }
    ## I think this is unreacheable now. Caught earlier.
    if(length(unique(geneMod$GeneNumID)) != nrow(geneMod)) {
        stop("Are there identical gene names in different modules?")
    }
    rownames(geneMod) <- 1:nrow(geneMod)
    geneMod   
}

geneModuleNull <- function(namesM) {
    v <- c("Root", setdiff(namesM, "Root"))
    names(v) <- v
    return(v)
}

list.of.deps <- function(x) {
    ## lookupTypeDep <- c("MN" = 1, "monotone" = 1,
    ##                 "SM" = 2, "semimonotone" = 2)
    lookupTypeDep <- c("MN" = "monotone",
                       "AND" = "monotone",
                       "monotone" = "monotone",
                       "CMPN" = "monotone",
                       "OR" = "semimonotone",
                       "SM" = "semimonotone",
                       "semimonotone" = "semimonotone",
                       "DMPN" = "semimonotone",
                       "XOR" = "xmpn",
                       "xmpn" = "xmpn",
                       "XMPN" = "xmpn",
                       "--"   = "--",
                       "-" = "--")
    ## FIXME: check values of typeDep

    if(length(x) > 1) {
        if(length(unique(x$s))!= 1)
            stop("Not all s identical within a child")
        if(length(unique(x$sh))!= 1)
            stop("Not all sh identical within a child")
        if(length(unique(x$typeDep))!= 1)
            stop("Not all typeDep identical within a child")
        if(length(unique(x$child))!= 1)
            stop("child not unique")
    }
    typeDep <- lookupTypeDep[unique(x$typeDep)]
    if(any(is_null_na(typeDep)))
        stop("typeDep value incorrect. Check spelling.")
    return(list(
        child = unique(x$child),
        s = unique(x$s),
        sh = unique(x$sh),
        typeDep = typeDep,
        parents = unlist(x$parent)))

}

to.long.rt <- function(rt, idm) {
    ## We now do this inconditionally, so that we do not need to use the
    ## "stringsAsFactors = FALSE". This is now done before
    ## if(is.numeric(rt$parent))
    ##     rt$parent <- as.character(rt$parent)
    ## if(is.numeric(rt$child))
    ##     rt$child <- as.character(rt$child)
   
    
    if(!("Root" %in% rt$parent))
        stop("Root must be one parent node")

    ## rt$parent <- unlist(lapply(rt$parent, nice.string))
    ## rt$child <- unlist(lapply(rt$child, nice.string))
   
    srt <- rt[order(rt$child), ]

    ## Not relevant if we allow non-numeric names
    ## all.child.genes <- as.integer(
    ##     unlist(lapply(rt[, 2],
    ##                   function(x) strsplit(x, ","))))
    ## ## check all childs
    ## if(!identical(sort(unique(all.child.genes)),
    ##               seq.int(max(all.child.genes))))
    ##     stop("Not all children present")
    long.rt <- lapply(split(srt, srt$child), list.of.deps)

    ## geneModule <- gene.to.module(srt)
    ## idm <- seq.int(length(names(long.rt)))
    ## names(idm) <- names(long.rt)
    ## idm <- c("0" = 0L, idm)
    ## geneModule$ModuleNumID <- idm[geneModule[, "Module"]]

    ## idm is just a look up table for the id of the module
    ## idm <- unique(geneModule$ModuleNumID)
    ## names(idm) <- unique(geneModule$Module)
    
    ## add integer IDs
    addIntID <- function(z, idm) {
        z$childNumID <- idm[z$child]
        z$parentsNumID <- idm[z$parents]
        if( any(is.na(z$parentsNumID)) ||
            any(is.na(z$childNumID)) ) {
            ## I think this can no longer be reached ever. Caught before.
            stop(paste("An ID is NA:",
                       "Is a gene part of two different modules?",
                       "(That includes being by itself and part",
                       "of a module.)"))
            
        }
        ## These checks could be somewhere else instead of here.
        if(length(unique(z$parentsNumID)) != length(z$parentsNumID))
            stop("No node can have the same parent multiple times")
        if(length(unique(z$parents)) != length(z$parents))
            stop("No node can have the same parent multiple times")
        if(length(z$child) > 1)
            stop("Child nodes can have one, and only one, member")
        
        ## sort the parents, always.
        o <- order(z$parentsNumID)
        z$parentsNumID <- z$parentsNumID[o]
        z$parents <- as.character(z$parents[o])
        return(z)
    }
    long.rt <- lapply(long.rt, function(x) addIntID(x, idm = idm))
   
    ## if(verbosity >= 4) {
    ##     message(paste("Number of drivers: ",
    ##                   length(unique(geneModule[, "Gene"]))))
    ##     message(paste("Number of modules: ",
    ##                   length(unique(geneModule[, "Module"]))))
    ## }
    return(long.rt)
    ## return(list(long.rt = long.rt, geneModule = geneModule))
}


epist.order.element <- function(x, y, sep, rm.sign = FALSE) {
    list(ids = nice.vector.eo(x, sep = sep, rm.sign = rm.sign), s = y)
}

oe.to.df <- function(x) {
    ma <- matrix(ncol = 2, nrow = 1 + length(x) - 2)
    if(length(x) == 2) {
        ma[1, 1] <- x[1]
        ma[1, 2] <- x[2]
    } else {
        ma[, 1] <- x[-length(x)]
        ma[, 2] <- x[-1]
    }
    return(data.frame(ma, stringsAsFactors = FALSE))
}


epist.order.to.pairs.modules <- function(x, sep, rm.sign = TRUE) {
    ## We discard, do not even accept, the coefficient
    tmp <- epist.order.element(x, -99, sep = sep, rm.sign = rm.sign)$ids
    if(length(tmp) > 1) {
        ## if a single gene, as when we specify all genotypes, we do not
        ## want this
        if(sep == ":")
            return(data.frame(combinations(n = length(tmp), r = 2, v = tmp),
                              stringsAsFactors = FALSE))
        else if(sep == ">") {
            return(oe.to.df(tmp))
        }
    }
}

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))),
        stringsAsFactors = TRUE)
    if(nrow(df) == 0L) { ## if only single genes in epist, we get nothing here.
        return(df)
    } else {
        colnames(df) <- c("parent", "child")
        if(sep == ":")
            df$typeDep <- "epistasis"
        else if(sep == ">")
            df$typeDep <- "orderEffect"
        return(unique(df))
    }
}



to.long.epist.order <- function(epor, sep, rm.sign = FALSE) {
    ## just vectors for now
    long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
                names(epor), epor)
    ## if(is.vector(epor))
    ##     long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
    ##                    names(epor), epor)
    ## else if(is.data.frame(epor)) 
    ##     long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
    ##                 as.character(epor$ids),
    ##                 epor$s)
    names(long) <- NULL
    return(long)
}

## addIntID.epist.order <- function(z, idm, sort) {
##     z$NumID <- idm[z$ids]
##     if(sort) {
##         ## essential for epistasis, but never do for order effects
##         o <- order(z$NumID)
##         z$NumID <- z$NumID[o]
##         z$ids <- z$ids[o]
##     }
##     return(z)
## }


addIntID.epist.order <- function(z, idm, sort, sign) {
    if( sort && (!sign))
        warning("sort is TRUE but sign is FALSE. You sure?")
    if((!sort) && sign)
        warning("sort is FALSE but sign is TRUE. You sure?")
    ## Adds the integer, but takes care of sign if sign = TRUE
    if(!sign) {
        z$NumID <- idm[z$ids]
        signs <- grep("^-", z$ids)
        if(length(signs))
            stop("You have a negative sign and you said sign = FALSE")
    } else {
        unsigned <- setdiff(unlist(lapply(z$ids,
                                          function(z) strsplit(z, "^-"))), "")
        NumID <- idm[unsigned]
        signs <- grep("^-", z$ids)
        if(length(signs)) {
            NumID[signs] <- -1 * NumID[signs]
        }
        z$NumID <- NumID
    }
    if(sort) {
        ## Essential for epistasis, but never do for order effects
        o <- order(z$NumID)
        z$NumID <- z$NumID[o]
        z$ids <- z$ids[o]
    }
    if(length(unique(z$NumID)) != length(z$NumID))
        stop("No node can have the same id multiple times")
    if(length(unique(z$ids)) != length(z$ids))
        stop("No node can have the same id multiple times")
    return(z)
}


checkRT <- function(mdeps) {
    if(ncol(mdeps) != 5)
        stop("mdeps must be of exactly 5 columns")
    if(!identical(colnames(mdeps), c("parent", "child", "s", "sh", "typeDep")))
        stop(paste("Column names of mdeps not of appropriate format. ",
                   "Should be parent, child, s, sh, typeDep"))
}


getNamesID <- function(fp) {
    ## Return a lookup table for names based on numeric IDs
    idname <- c(fp$geneModule$GeneNumID,
                fp$long.geneNoInt$GeneNumID,
                fp$fitnessLandscape_gene_id$GeneNumID)
    names(idname) <- c(fp$geneModule$Gene,
                       fp$long.geneNoInt$Gene,
                       fp$fitnessLandscape_gene_id$Gene)
    return(idname[-1]) ## remove Root!!
}


getGeneIDNum <- function(geneModule, geneNoInt, fitnessLandscape_gene_id,
                         drv, sort = TRUE) {
    ## It returns the genes, as NumID, in the given vector with names drv
    ## initMutant uses this, for simplicity, without sorting, but noInt
    ## are always sorted

    ## Also used for the drivers with sort = TRUE

    ## Yes, we must do it twice because we do not know before hand which
    ## is which. This makes sure no NA. Period.
    if(any(is.na( match(drv, c(geneModule$Gene, geneNoInt$Gene,
                               fitnessLandscape_gene_id$Gene))))) {
        stop(paste("For driver or initMutant you have passed genes",
                   "not in the fitness table."))
    }
    
    indicesM <- as.vector(na.omit(match( drv, geneModule$Gene)))
    indicesI <- as.vector(na.omit(sort(match( drv, geneNoInt$Gene))))
    indicesF <- as.vector(na.omit(sort(match( drv, fitnessLandscape_gene_id$Gene))))
    if(sort) {
        indicesM <- sort(indicesM)
    }
    return(c(
        geneModule$GeneNumID[indicesM],
        geneNoInt$GeneNumID[indicesI],
        fitnessLandscape_gene_id$GeneNumID[indicesF])
    )
}


allFitnessORMutatorEffects <- function(rT = NULL,
                                       epistasis = NULL,
                                       orderEffects = NULL,
                                       noIntGenes = NULL,
                                       geneToModule = NULL,
                                       drvNames = NULL,
                                       keepInput = TRUE,
                                       genotFitness = NULL,
                                       ## refFE = NULL,
                                       calledBy = NULL) {
    ## From allFitnessEffects. Generalized so we deal with Fitness
    ## and mutator.
    
    ## restrictions: the usual rt

    ## epistasis: as it says, with the ":"

    ## orderEffects: the ">"
    
    ## All of the above can be genes or can be modules (if you pass a
    ## geneToModule)

    ## rest: rest of genes, with fitness


    ## For epistasis and order effects we create the output object but
    ## missing the numeric ids of genes. With rT we do it in one go, as we
    ## already know the mapping of genes to numeric ids. We could do the
    ## same in epistasis and order, but we would be splitting twice
    ## (whereas for rT extracting the names is very simple).

    ## called appropriately?
    if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
        stop("How did you call this function?. Bug.")
    
    if(calledBy == "allMutatorEffects") {
        ## very paranoid check
        if( !is.null(rT) || !is.null(orderEffects) ||
            !is.null(drvNames) || !is.null(genotFitness))
            stop("allMutatorEffects called with forbidden arguments.",
                 "Is this an attempt to subvert the function?")
    }
    
    rtNames <- NULL
    epiNames <- NULL
    orNames <- NULL
    if(!is.null(rT)) {
        ## This is really ugly, but to prevent the stringsAsFactors I need it here:
        rT$parent <- as.character(rT$parent)
        rT$child <- as.character(rT$child)
        rT$typeDep <- as.character(rT$typeDep)
        rtNames <- unique(c(rT$parent, rT$child))
    }
    if(!is.null(epistasis)) {
        long.epistasis <- to.long.epist.order(epistasis, ":")
        ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
        ## deal with the possible negative signs
        epiNames <- setdiff(unique(
            unlist(lapply(long.epistasis,
                          function(x) lapply(x$ids,
                                             function(z) strsplit(z, "^-"))))),
                            "")
    } else {
        long.epistasis <- list()
    }
    if(!is.null(orderEffects)) {
        long.orderEffects <- to.long.epist.order(orderEffects, ">")
        orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
    } else {
        long.orderEffects <- list()
    }
    allModuleNames <- unique(c(rtNames, epiNames, orNames))
    if(is.null(geneToModule)) {
        gMOneToOne <- TRUE
        geneToModule <- geneModuleNull(allModuleNames)
    } else {
        gMOneToOne <- FALSE
        if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
            stop(paste("Some values in geneToModule not present in any of",
                       " rT, epistasis, or order effects"))
        if(any(is.na(match(allModuleNames, names(geneToModule)))))
            stop(paste("Some values in rT, epistasis, ",
                       "or order effects not in geneToModule"))
    }
    geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
    
    idm <- unique(geneModule$ModuleNumID)
    names(idm) <- unique(geneModule$Module)

    if(!is.null(rT)) {
        checkRT(rT)
        long.rt <- to.long.rt(rT, idm)
    } else {
        long.rt <- list() ## yes, we want an object of length 0
    }

    ## Append the numeric ids to epistasis and order
    if(!is.null(epistasis)) {
        long.epistasis <- lapply(long.epistasis,
                                 function(x)
                                     addIntID.epist.order(x, idm,
                                                          sort = TRUE,
                                                          sign = TRUE))
    }
    if(!is.null(orderEffects)) {
        long.orderEffects <- lapply(long.orderEffects,
                                    function(x)
                                        addIntID.epist.order(x, idm,
                                                             sort = FALSE,
                                                             sign = FALSE))
    }
    
    if(!is.null(noIntGenes)) {
        if(inherits(noIntGenes, "character")) {
            wm <- paste("noIntGenes is a character vector.",
                        "This is probably not what you want, and will",
                        "likely result in an error downstream.",
                        "You can get messages like",
                        " 'not compatible with requested type', and others.",
                        "We are stopping.")
            stop(wm)
        }
            
        mg <- max(geneModule[, "GeneNumID"])
        gnum <- seq_along(noIntGenes) + mg
        if(!is.null(names(noIntGenes))) {
            ng <- names(noIntGenes)
            if( any(grepl(",", ng, fixed = TRUE)) ||
                any(grepl(">", ng, fixed = TRUE)) ||
                any(grepl(":", ng, fixed = TRUE))   )
                stop("The name of some noIntGenes contain a ',' or a '>' or a ':'")
            if(any(ng %in% geneModule[, "Gene"] ))
                stop("A gene in noIntGenes also present in the other terms")
            if(any(duplicated(ng)))
                stop("Duplicated gene names in geneNoInt")
            if(any(is.na(ng)))
                stop("In noIntGenes some genes have names, some don't.",
                     " Name all of them, or name none of them.")
        } else {
            ng <- gnum
        }
        geneNoInt <- data.frame(Gene = as.character(ng),
                                GeneNumID = gnum,
                                s = noIntGenes,
                                stringsAsFactors = FALSE)
    } else {
        geneNoInt <- data.frame()
    }
    
    if(is.null(genotFitness)) {
        genotFitness <- matrix(NA, nrow = 0, ncol = 1)
        fitnessLandscape_df <- data.frame()
        fitnessLandscape_gene_id <- data.frame()
    } else {
        ## Yes, I am duplicating stuff for now.
        ## This makes life simpler in C++:
        ## In the map, the key is the genotype name, as
        ## cnn <- colnames(genotFitness)[-ncol(genotFitness)]
        cnn <- 1:(ncol(genotFitness) - 1)
        gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
                     function(x) paste(cnn[as.logical(x)],
                                       collapse = ", "))
        ## rownames(genotFitness) <- gfn
        fitnessLandscape_df <-
            data.frame(Genotype = gfn,
                       Fitness = genotFitness[, ncol(genotFitness)],
                       stringsAsFactors = FALSE)
        fitnessLandscape_gene_id <- data.frame(
            Gene = colnames(genotFitness)[-ncol(genotFitness)],
            GeneNumID = cnn,
            stringsAsFactors = FALSE)
        
    }
    
    if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
             nrow(geneNoInt) + nrow(genotFitness)) == 0)
        stop("You have specified nothing!")

    if(calledBy == "allFitnessEffects") {
        if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
            graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
        } else {
            graphE <- NULL
        }
    } else {
        graphE <- NULL
    }
    if(!is.null(drvNames)) {
        drv <- unique(getGeneIDNum(geneModule, geneNoInt, fitnessLandscape_gene_id,
                                   drvNames))
        ## drivers should never be in the geneNoInt; Why!!!???
        ## Catch the problem. This is an overkill,
        ## so since we catch the issue, we could leave the geneNoInt. But
        ## that should not be there in this call.
        ## if(any(drvNames %in% geneNoInt$Gene)) {
        ##     stop(paste("At least one gene in drvNames is a geneNoInt gene.",
        ##                "That is not allowed.",
        ##                "If that gene is a driver, pass it as gene in the epistasis",
        ##                "component."))
        ## }
        ## drv <- getGeneIDNum(geneModule, NULL, drvNames)
    } else {
        ## we used to have this default
        ## drv <- geneModule$GeneNumID[-1]
        drv <- vector(mode = "integer", length = 0L)
    }
  
    if(!keepInput) {
        rT <- epistasis <- orderEffects <- noIntGenes <- NULL
    }

    out <- list(long.rt = long.rt,
                long.epistasis = long.epistasis,
                long.orderEffects = long.orderEffects,
                long.geneNoInt = geneNoInt,
                geneModule = geneModule,
                gMOneToOne = gMOneToOne,
                geneToModule = geneToModule,
                graph = graphE,
                drv = drv,
                rT = rT,
                epistasis = epistasis,
                orderEffects = orderEffects,
                noIntGenes = noIntGenes,
                fitnessLandscape = genotFitness,
                fitnessLandscape_df = fitnessLandscape_df,
                fitnessLandscape_gene_id = fitnessLandscape_gene_id
                )
    if(calledBy == "allFitnessEffects") {
        class(out) <- c("fitnessEffects")
    } else if(calledBy == "allMutatorEffects") {
        class(out) <- c("mutatorEffects")
    }
    return(out)
}

## Former version, with fitness landscape
## allFitnessORMutatorEffects <- function(rT = NULL,
##                                        epistasis = NULL,
##                                        orderEffects = NULL,
##                                        noIntGenes = NULL,
##                                        geneToModule = NULL,
##                                        drvNames = NULL,
##                                        keepInput = TRUE,
##                                        ## refFE = NULL,
##                                        calledBy = NULL) {
##     ## From allFitnessEffects. Generalized so we deal with Fitness
##     ## and mutator.
    
##     ## restrictions: the usual rt

##     ## epistasis: as it says, with the ":"

##     ## orderEffects: the ">"
    
##     ## All of the above can be genes or can be modules (if you pass a
##     ## geneToModule)

##     ## rest: rest of genes, with fitness


##     ## For epistasis and order effects we create the output object but
##     ## missing the numeric ids of genes. With rT we do it in one go, as we
##     ## already know the mapping of genes to numeric ids. We could do the
##     ## same in epistasis and order, but we would be splitting twice
##     ## (whereas for rT extracting the names is very simple).

##     ## called appropriately?
##     if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
##         stop("How did you call this function?. Bug.")
    
##     if(calledBy == "allMutatorEffects") {
##         ## very paranoid check
##         if( !is.null(rT) || !is.null(orderEffects) || !is.null(drvNames))
##             stop("allMutatorEffects called with forbidden arguments.",
##                  "Is this an attempt to subvert the function?")
##     }
    
##     rtNames <- NULL
##     epiNames <- NULL
##     orNames <- NULL
##     if(!is.null(rT)) {
##         ## This is really ugly, but to prevent the stringsAsFactors I need it here:
##         rT$parent <- as.character(rT$parent)
##         rT$child <- as.character(rT$child)
##         rT$typeDep <- as.character(rT$typeDep)
##         rtNames <- unique(c(rT$parent, rT$child))
##     }
##     if(!is.null(epistasis)) {
##         long.epistasis <- to.long.epist.order(epistasis, ":")
##         ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
##         ## deal with the possible negative signs
##         epiNames <- setdiff(unique(
##             unlist(lapply(long.epistasis,
##                           function(x) lapply(x$ids,
##                                              function(z) strsplit(z, "^-"))))),
##                             "")
##     } else {
##         long.epistasis <- list()
##     }
##     if(!is.null(orderEffects)) {
##         long.orderEffects <- to.long.epist.order(orderEffects, ">")
##         orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
##     } else {
##         long.orderEffects <- list()
##     }
##     allModuleNames <- unique(c(rtNames, epiNames, orNames))
##     if(is.null(geneToModule)) {
##         gMOneToOne <- TRUE
##         geneToModule <- geneModuleNull(allModuleNames)
##     } else {
##         gMOneToOne <- FALSE
##         if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
##             stop(paste("Some values in geneToModule not present in any of",
##                        " rT, epistasis, or order effects"))
##         if(any(is.na(match(allModuleNames, names(geneToModule)))))
##             stop(paste("Some values in rT, epistasis, ",
##                        "or order effects not in geneToModule"))
##     }
##     geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
    
##     idm <- unique(geneModule$ModuleNumID)
##     names(idm) <- unique(geneModule$Module)

##     if(!is.null(rT)) {
##         checkRT(rT)
##         long.rt <- to.long.rt(rT, idm)
##     } else {
##         long.rt <- list() ## yes, we want an object of length 0
##     }

##     ## Append the numeric ids to epistasis and order
##     if(!is.null(epistasis)) {
##         long.epistasis <- lapply(long.epistasis,
##                                  function(x)
##                                      addIntID.epist.order(x, idm,
##                                                           sort = TRUE,
##                                                           sign = TRUE))
##     }
##     if(!is.null(orderEffects)) {
##         long.orderEffects <- lapply(long.orderEffects,
##                                     function(x)
##                                         addIntID.epist.order(x, idm,
##                                                              sort = FALSE,
##                                                              sign = FALSE))
##     }
    
##     if(!is.null(noIntGenes)) {
##         if(inherits(noIntGenes, "character")) {
##             wm <- paste("noIntGenes is a character vector.",
##                         "This is probably not what you want, and will",
##                         "likely result in an error downstream.",
##                         "You can get messages like",
##                         " 'not compatible with requested type', and others.",
##                         "We are stopping.")
##             stop(wm)
##         }
            
##         mg <- max(geneModule[, "GeneNumID"])
##         gnum <- seq_along(noIntGenes) + mg
##         if(!is.null(names(noIntGenes))) {
##             ng <- names(noIntGenes)
##             if( grepl(",", ng, fixed = TRUE) || grepl(">", ng, fixed = TRUE)
##                 || grepl(":", ng, fixed = TRUE))
##                 stop("The name of some noIntGenes contain a ',' or a '>' or a ':'")
##             if(any(ng %in% geneModule[, "Gene"] ))
##                 stop("A gene in noIntGenes also present in the other terms")
##             if(any(duplicated(ng)))
##                 stop("Duplicated gene names in geneNoInt")
##             if(any(is.na(ng)))
##                 stop("In noIntGenes some genes have names, some don't.",
##                      " Name all of them, or name none of them.")
##         } else {
##             ng <- gnum
##         }
##         geneNoInt <- data.frame(Gene = as.character(ng),
##                                 GeneNumID = gnum,
##                                 s = noIntGenes,
##                                 stringsAsFactors = FALSE)
##     } else {
##         geneNoInt <- data.frame()
##     }

##     if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
##              nrow(geneNoInt)) == 0)
##         stop("You have specified nothing!")

##     if(calledBy == "allFitnessEffects") {
##         if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
##             graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
##         } else {
##             graphE <- NULL
##         }
##     } else {
##         graphE <- NULL
##     }
##     if(!is.null(drvNames)) {
##         drv <- unique(getGeneIDNum(geneModule, geneNoInt, drvNames))
##         ## drivers should never be in the geneNoInt; Why!!!???
##         ## Catch the problem. This is an overkill,
##         ## so since we catch the issue, we could leave the geneNoInt. But
##         ## that should not be there in this call.
##         ## if(any(drvNames %in% geneNoInt$Gene)) {
##         ##     stop(paste("At least one gene in drvNames is a geneNoInt gene.",
##         ##                "That is not allowed.",
##         ##                "If that gene is a driver, pass it as gene in the epistasis",
##         ##                "component."))
##         ## }
##         ## drv <- getGeneIDNum(geneModule, NULL, drvNames)
##     } else {
##         ## we used to have this default
##         ## drv <- geneModule$GeneNumID[-1]
##         drv <- vector(mode = "integer", length = 0L)
##     }
    
##     if(!keepInput) {
##         rT <- epistasis <- orderEffects <- noIntGenes <- NULL
##     }
##     out <- list(long.rt = long.rt,
##                 long.epistasis = long.epistasis,
##                 long.orderEffects = long.orderEffects,
##                 long.geneNoInt = geneNoInt,
##                 geneModule = geneModule,
##                 gMOneToOne = gMOneToOne,
##                 geneToModule = geneToModule,
##                 graph = graphE,
##                 drv = drv,
##                 rT = rT,
##                 epistasis = epistasis,
##                 orderEffects = orderEffects,
##                 noIntGenes = noIntGenes                
##                 )
##     if(calledBy == "allFitnessEffects") {
##         class(out) <- c("fitnessEffects")
##     } else if(calledBy == "allMutatorEffects") {
##         class(out) <- c("mutatorEffects")
##     }
##     return(out)
## }


allFitnessEffects <- function(rT = NULL,
                              epistasis = NULL,
                              orderEffects = NULL,
                              noIntGenes = NULL,
                              geneToModule = NULL,
                              drvNames = NULL,
                              genotFitness = NULL,
                              keepInput = TRUE) {

    if(!is.null(genotFitness)) {
        if(!is.null(rT) || !is.null(epistasis) ||
           !is.null(orderEffects) || !is.null(noIntGenes) ||
           !is.null(geneToModule)) {
            stop("You have a non-null genotFitness.",
                 " If you pass the complete genotype to fitness mapping",
                 " you cannot pass any of rT, epistasis, orderEffects",
                 " noIntGenes or geneToModule.")
        }

        genotFitness_std <- to_genotFitness_std(genotFitness, simplify = TRUE)
        ## epistasis <- from_genotype_fitness(genotFitness)
    } else {
        genotFitness_std <- NULL
    }
    allFitnessORMutatorEffects(
        rT = rT,
        epistasis = epistasis,
        orderEffects = orderEffects,
        noIntGenes = noIntGenes,
        geneToModule = geneToModule,
        drvNames = drvNames,
        keepInput = keepInput,
        genotFitness = genotFitness_std,
        calledBy = "allFitnessEffects")
}

## Former version
## allFitnessEffects <- function(rT = NULL,
##                               epistasis = NULL,
##                               orderEffects = NULL,
##                               noIntGenes = NULL,
##                               geneToModule = NULL,
##                               drvNames = NULL,
##                               genotFitness = NULL,
##                               keepInput = TRUE) {

##     if(!is.null(genotFitness)) {
##         if(!is.null(rT) || !is.null(epistasis) ||
##            !is.null(orderEffects) || !is.null(noIntGenes) ||
##            !is.null(geneToModule)) {
##             stop("You have a non-null genotFitness.",
##                  " If you pass the complete genotype to fitness mapping",
##                  " you cannot pass any of rT, epistasis, orderEffects",
##                  " noIntGenes or geneToModule.")
##         }
##         epistasis <- from_genotype_fitness(genotFitness)
##     }
##     allFitnessORMutatorEffects(
##         rT = rT,
##         epistasis = epistasis,
##         orderEffects = orderEffects,
##         noIntGenes = noIntGenes,
##         geneToModule = geneToModule,
##         drvNames = drvNames,
##         keepInput = keepInput,
##         calledBy = "allFitnessEffects")
## }

## allFitnessEffects <- function(rT = NULL,
##                               epistasis = NULL,
##                               orderEffects = NULL,
##                               noIntGenes = NULL,
##                               geneToModule = NULL,
##                               drvNames = NULL,
##                               keepInput = TRUE) {
##     ## restrictions: the usual rt

##     ## epistasis: as it says, with the ":"

##     ## orderEffects: the ">"
    
##     ## All of the above can be genes or can be modules (if you pass a
##     ## geneToModule)

##     ## rest: rest of genes, with fitness


##     ## For epistasis and order effects we create the output object but
##     ## missing the numeric ids of genes. With rT we do it in one go, as we
##     ## already know the mapping of genes to numeric ids. We could do the
##     ## same in epistasis and order, but we would be splitting twice
##     ## (whereas for rT extracting the names is very simple).

    
##     rtNames <- NULL
##     epiNames <- NULL
##     orNames <- NULL
##     if(!is.null(rT)) {
##         ## This is really ugly, but to prevent the stringsAsFactors I need it here:
##         rT$parent <- as.character(rT$parent)
##         rT$child <- as.character(rT$child)
##         rT$typeDep <- as.character(rT$typeDep)
##         rtNames <- unique(c(rT$parent, rT$child))
##     }
##     if(!is.null(epistasis)) {
##         long.epistasis <- to.long.epist.order(epistasis, ":")
##         ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
##         ## deal with the possible negative signs
##         epiNames <- setdiff(unique(
##             unlist(lapply(long.epistasis,
##                           function(x) lapply(x$ids,
##                                              function(z) strsplit(z, "^-"))))),
##                             "")
        
##     } else {
##         long.epistasis <- list()
##     }
##     if(!is.null(orderEffects)) {
##         long.orderEffects <- to.long.epist.order(orderEffects, ">")
##         orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
##     } else {
##         long.orderEffects <- list()
##     }
##     allModuleNames <- unique(c(rtNames, epiNames, orNames))
##     if(is.null(geneToModule)) {
##         gMOneToOne <- TRUE
##         geneToModule <- geneModuleNull(allModuleNames)
##     } else {
##         gMOneToOne <- FALSE
##         if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
##             stop(paste("Some values in geneToModule not present in any of",
##                        " rT, epistasis, or order effects"))
##         if(any(is.na(match(allModuleNames, names(geneToModule)))))
##             stop(paste("Some values in rT, epistasis, ",
##                        "or order effects not in geneToModule"))
##     }
##     geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)

##     idm <- unique(geneModule$ModuleNumID)
##     names(idm) <- unique(geneModule$Module)

##     if(!is.null(rT)) {
##         checkRT(rT)
##         long.rt <- to.long.rt(rT, idm)
##     } else {
##         long.rt <- list() ## yes, we want an object of length 0
##     }

##     ## Append the numeric ids to epistasis and order
##     if(!is.null(epistasis)) {
##         long.epistasis <- lapply(long.epistasis,
##                                  function(x)
##                                      addIntID.epist.order(x, idm,
##                                                           sort = TRUE,
##                                                           sign = TRUE))
##     }
##     if(!is.null(orderEffects)) {
##         long.orderEffects <- lapply(long.orderEffects,
##                                     function(x)
##                                         addIntID.epist.order(x, idm,
##                                                              sort = FALSE,
##                                                              sign = FALSE))
##     }
    
##     if(!is.null(noIntGenes)) {
##         mg <- max(geneModule[, "GeneNumID"])
##         gnum <- seq_along(noIntGenes) + mg
##         if(!is.null(names(noIntGenes))) {
##             ng <- names(noIntGenes)
##             if(any(ng %in% geneModule[, "Gene"] ))
##                 stop("A gene in noIntGenes also present in the other terms")
##         } else {
##             ng <- gnum
##         }
##         geneNoInt <- data.frame(Gene = as.character(ng),
##                                 GeneNumID = gnum,
##                                 s = noIntGenes,
##                                 stringsAsFactors = FALSE)
##     } else {
##         geneNoInt <- data.frame()
##     }
##     if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
##              nrow(geneNoInt)) == 0)
##         stop("You have specified nothing!")

##     if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
##         graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
##     } else {
##         graphE <- NULL
##     }

##     if(!is.null(drvNames)) {
##         drv <- getGeneIDNum(geneModule, geneNoInt, drvNames)
##     } else {
##         drv <- geneModule$GeneNumID[-1]
##     }
##     if(!keepInput) {
##         rT <- epistasis <- orderEffects <- noIntGenes <- NULL
##     }
##     out <- list(long.rt = long.rt,
##                 long.epistasis = long.epistasis,
##                 long.orderEffects = long.orderEffects,
##                 long.geneNoInt = geneNoInt,
##                 geneModule = geneModule,
##                 gMOneToOne = gMOneToOne,
##                 geneToModule = geneToModule,
##                 graph = graphE,
##                 drv = drv,
##                 rT = rT,
##                 epistasis = epistasis,
##                 orderEffects = orderEffects,
##                 noIntGenes = noIntGenes                
##                 )
##     class(out) <- c("fitnessEffects")
##     return(out)
## }


## No longer used
## rtAndGeneModule <- function(mdeps, gM = NULL) {
##     ## To show a table of restrictions when there are modules. Do not use
##     ## for anything else. Maybe as intermediate to plotting.
    
##     ## Specify restriction table of modules and a mapping of modules to
##     ## genes. gM is a named vector; names are modules, values are elements
##     ## of each module.

##     ## We do nothing important if gM is NULL except checks

##     ## If there are modules, the table shows the individual genes.
##     checkRT(mdeps)
##     ## if(ncol(mdeps) != 5)
##     ##     stop("mdeps must be of exactly 5 columns")
##     ## if(!identical(colnames(mdeps), c("parent", "child", "s", "sh", "typeDep")))
##     ##     stop(paste("Column names of mdeps not of appropriate format. ",
##     ##                "Should be parent, child, s, sh, typeDep"))
##     if(!is.null(gM)) {
##         if(any(is.na(match(mdeps[ , 1], names(gM)))))
##             stop("Some values in parent not from a known module")
##         if(any(is.na(match(mdeps[ , 2], names(gM)))))
##             stop("Some values in child not from a known module")
##         if(any(is.na(match(names(gM), c(mdeps[, 1], mdeps[, 2])))))
##             stop("Some values in module in neither parent or child")
        
##         parent <- gM[mdeps[, 1]]
##         child <- gM[mdeps[, 2]]
##         df <- data.frame(parent = parent,
##                          child = child,
##                          s = mdeps$s,
##                          sh = mdeps$sh,
##                          typeDep = mdeps$typeDep,
##                          stringsAsFactors = FALSE)
##     } else {
##         df <- mdeps
##     }
##     rownames(df) <- seq.int(nrow(df))
##     return(df)
## }

## wrap.test.rt <- function(rt, gM = NULL) {
##     ## FIXME add epistasis and orderEffects
##     lrt <- allFitnessEffects(rt, geneToModule = gM)
##     ## wrap_test_rt(lrt$long.rt)
##     wrap_test_rt(lrt$long.rt, lrt$geneModule)
## }

## No longer used
## wrap.readFitnessEffects <- function(rt, epi, oe, ni, gm, echo = TRUE) {
##     tt <- allFitnessEffects(rt, epi, oe, ni, gm)
##     readFitnessEffects(tt, echo = echo)
##     ## readFitnessEffects(tt$long.rt,
##     ##                    tt$long.epistasis,
##     ##                    tt$long.orderEffects,
##     ##                    tt$long.geneNoInt,
##     ##                    tt$geneModule,
##     ##                    tt$gMOneToOne,
##     ##                    echo = TRUE)
## }



evalGenotypeORMut <- function(genotype,
                              fmEffects,
                              verbose = FALSE,
                              echo = FALSE,
                              model = "",
                              calledBy_= NULL) {
    ## genotype can be a vector of integers, that are the exact same in
    ## the table of fmEffects or a vector of strings, or a vector (a
    ## string) with genes separated by "," or ">"

    if( !(calledBy_ %in% c("evalGenotype", "evalGenotypeMut") ))
        stop("How did you call this function?. Bug.")

    ## fmEffects could be a mutator effect
    if(!exists("fitnessLandscape_gene_id", where = fmEffects)) {
        fmEffects$fitnessLandscape_df <- data.frame()
        fmEffects$fitnessLandscape_gene_id <- data.frame()
    }

    if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
        (nrow(fmEffects$fitnessLandscape_df) > 0)) {
        warning("Bozic model passing a fitness landscape will not work",
                    " for now.")
    }
    
    if(echo)
        cat(paste("Genotype: ", genotype))
    if(!is.integer(genotype)) {
        if(length(grep(">", genotype))) {
            genotype <- nice.vector.eo(genotype, ">")
        } else if(length(grep(",", genotype))) {
            genotype <- nice.vector.eo(genotype, ",")
        }
        all.g.nums <- c(fmEffects$geneModule$GeneNumID,
                        fmEffects$long.geneNoInt$GeneNumID,
                        fmEffects$fitnessLandscape_gene_id$GeneNumID)
        all.g.names <- c(fmEffects$geneModule$Gene,
                         fmEffects$long.geneNoInt$Gene,
                         fmEffects$fitnessLandscape_gene_id$Gene)
        genotype <- all.g.nums[match(genotype, all.g.names)]
    }
    if(any(is.na(genotype)))
        stop("genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
    if(!length(genotype))
        stop("genotypes must have at least one mutated gene")
    if(any(genotype < 0)) {
        stop(paste("genotypes cannot contain negative values.",
                   "If you see this message, you found a bug."))
    }
    if(model %in% c("Bozic", "bozic1", "bozic2") )
        prodNeg <- TRUE
    else
        prodNeg <- FALSE
    ff <- evalRGenotype(rG = genotype,
                        rFE = fmEffects,
                        verbose = verbose,
                        prodNeg = prodNeg,
                        calledBy_ = calledBy_)

    
    if(echo) {
        if(calledBy_ == "evalGenotype") {
            if(!prodNeg)
                cat(" Fitness: ", ff, "\n")
            else
                cat(" Death rate: ", ff, "\n")
        } else if(calledBy_ == "evalGenotypeMut") {
            cat(" Mutation rate product :", ff, "\n")
        }
        
    } 
    return(ff)
}

evalGenotype <- function(genotype, fitnessEffects,
                         verbose = FALSE,
                         echo = FALSE,
                         model = "") {
    if(inherits(fitnessEffects, "mutatorEffects"))
         stop("You are trying to get the fitness of a mutator specification. ",
             "You did not pass an object of class fitnessEffects.")

    evalGenotypeORMut(genotype = genotype,
                       fmEffects = fitnessEffects,
                       verbose = verbose,
                       echo = echo,
                       model  = model ,
                       calledBy_= "evalGenotype"
                       )
}


evalGenotypeFitAndMut <- function(genotype,
                                  fitnessEffects,
                                  mutatorEffects,
                                  verbose = FALSE,
                                  echo = FALSE,
                                  model = "") {
    
    ## Must deal with objects from previous, pre flfast, modifications
    if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
        fitnessEffects$fitnessLandscape_df <- data.frame()
        fitnessEffects$fitnessLandscape_gene_id <- data.frame()
    }
    if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
        (nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
        warning("Bozic model passing a fitness landscape will not work",
                    " for now.")
    }
    prodNeg <- FALSE
    ## Next is from evalGenotypeAndMut
    if(echo)
        cat(paste("Genotype: ", genotype))
    if(!is.integer(genotype)) {
        if(length(grep(">", genotype))) {
            genotype <- nice.vector.eo(genotype, ">")
        } else if(length(grep(",", genotype))) {
            genotype <- nice.vector.eo(genotype, ",")
        }
        all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
                        fitnessEffects$long.geneNoInt$GeneNumID,
                        fitnessEffects$fitnessLandscape_gene_id$GeneNumID)
        all.g.names <- c(fitnessEffects$geneModule$Gene,
                         fitnessEffects$long.geneNoInt$Gene,
                         fitnessEffects$fitnessLandscape_gene_id$Gene)
        genotype <- all.g.nums[match(genotype, all.g.names)]
    }
    if(any(is.na(genotype)))
        stop("genotype contains NAs or genes not in fitnessEffects")
    if(!length(genotype))
        stop("genotypes must have at least one mutated gene")
    if(any(genotype < 0)) {
        stop(paste("genotypes cannot contain negative values.",
                   "If you see this message, you found a bug."))
    }

    full2mutator_ <- matchGeneIDs(mutatorEffects,
                                  fitnessEffects)$Reduced
    if(model %in% c("Bozic", "bozic1", "bozic2") )
        prodNeg <- TRUE
    else
        prodNeg <- FALSE
    evalRGenotypeAndMut(genotype, fitnessEffects,
                        mutatorEffects, full2mutator_,
                        verbose = verbose,
                        prodNeg = prodNeg)
}

## evalGenotype <- function(genotype, fitnessEffects,
##                          verbose = FALSE,
##                          echo = FALSE,
##                          model = "") {
##     ## genotype can be a vector of integers, that are the exact same in
##     ## the table of fitnessEffects or a vector of strings, or a vector (a
##     ## string) with genes separated by "," or ">"
    
##     if(echo)
##         cat(paste("Genotype: ", genotype))
##     if(!is.integer(genotype)) {
##         if(length(grep(">", genotype))) {
##             genotype <- nice.vector.eo(genotype, ">")
##         } else if(length(grep(",", genotype))) {
##             genotype <- nice.vector.eo(genotype, ",")
##         }
##         all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
##                         fitnessEffects$long.geneNoInt$GeneNumID)
##         all.g.names <- c(fitnessEffects$geneModule$Gene,
##                          fitnessEffects$long.geneNoInt$Gene)
##         genotype <- all.g.nums[match(genotype, all.g.names)]
##     }
##     if(any(is.na(genotype)))
##         stop("genotype contains NAs or genes not in fitnessEffects")
##     if(!length(genotype))
##         stop("genotypes must have at least one mutated gene")
##     if(any(genotype < 0)) {
##         stop(paste("genotypes cannot contain negative values.",
##                    "If you see this message, you found a bug."))
##     }
##     if(model %in% c("Bozic", "bozic1", "bozic2") )
##         prodNeg <- TRUE
##     else
##         prodNeg <- FALSE
##     ff <- evalRGenotype(genotype, fitnessEffects, verbose, prodNeg)


##     if(echo) {
##         if(!prodNeg)
##             cat(" Fitness: ", ff, "\n")
##         else
##             cat(" Death rate: ", ff, "\n")
##     } ## else {
##     ##     return(ff)
##     ## }
##     return(ff)
## }

## For multiple genotypes, lapply the matching.
## Nope, I think unneeded
## internal.convert_genotypes <- function(genotypes, gm) {
##     genotypes <- lapply(lg, function(x) gm$GeneNumID[match(x, gm$Gene)])
## }

## I am here: simplify this

evalAllGenotypesORMut <- function(fmEffects,
                                  order = FALSE, max = 256,
                             addwt = FALSE,
                             model = "",
                             calledBy_ = "") {
##                             minimal = FALSE) {
    if( !(calledBy_ %in% c("evalGenotype", "evalGenotypeMut") ))
        stop("How did you call this function?. Bug.")
    
    if( (calledBy_ == "evalGenotype" ) &&
        (!inherits(fmEffects, "fitnessEffects")))
        stop("You are trying to get the fitness of a mutator specification. ",
             "You did not pass an object of class fitnessEffects.")
    if( (calledBy_ == "evalGenotypeMut" ) &&
        (!inherits(fmEffects, "mutatorEffects")))
        stop("You are trying to get the mutator effects of a fitness specification. ",
             "You did not pass an object of class mutatorEffects.")

    
    
    ## if(!minimal)
    allg <- generateAllGenotypes(fitnessEffects = fmEffects,
                                 order = order, max = max)
    ## else
        ## allg <- generateAllGenotypes_minimal(fitnessEffects = fmEffects,
        ##                                      max = max)
    ## if(order)
    ##     tot <- function(n) {sum(sapply(seq.int(n),
    ##                                    function(x) choose(n, x) * factorial(x)))}
    ## else
    ##     tot <- function(n) {2^n}
    ## nn <- nrow(fitnessEffects$geneModule) -1  + nrow(fitnessEffects$long.geneNoInt)
    ## tnn <- tot(nn)
    ## if(tnn > max) {
    ##     m <- paste("There are", tnn, "genotypes.")
    ##     m <- paste(m, "This is larger than max.")
    ##     m <- paste(m, "Adjust max and rerun if you want")
    ##     stop(m)
    ## }
    ## ## With mutator, the ids of genes need not go from 1:n
    ## vid <- allNamedGenes(fitnessEffects)$GeneNumID
    ## if(order) {
    ##     f1 <- function(n) {
    ##         lapply(seq.int(n), function(x) permutations(n = n, r = x, v = vid))
    ##     }
    ## } else {
    ##     f1 <- function(n) {
    ##         lapply(seq.int(n), function(x) combinations(n = n, r = x, v = vid))}
        
    ## }
    ## genotNums <- f1(nn)
    ## list.of.vectors <- function(y) {
    ##     ## there's got to be a simpler way
    ##     lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
    ##            function(m) m[[1]])
    ## }
    ## genotNums <- list.of.vectors(genotNums)
    ## names <- c(fitnessEffects$geneModule$Gene[-1],
    ##            fitnessEffects$long.geneNoInt$Gene)

    ## genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
    ##                      function(z)
    ##                          paste(z,
    ##                                collapse = if(order){" > "} else {", "} )))
    ## This ain't the best, as we repeatedly read and convert
    ## fitnessEffects.  If this were slow, prepare C++ function that takes
    ## vectors of genotypes and uses same fitnessEffects.


    if(model %in% c("Bozic", "bozic1", "bozic2") )
        prodNeg <- TRUE
    else
        prodNeg <- FALSE
    allf <- vapply(allg$genotNums,
                   function(x) evalRGenotype(x, fmEffects, FALSE,
                                             prodNeg,
                                             calledBy_),
                   1.1)
    df <- data.frame(Genotype = allg$genotNames, Fitness = allf,
                     stringsAsFactors = FALSE)
    ## Why am I doing this? I am not computing the genotype.  I test the
    ## evaluation of the empty genotype in the tests. But its evaluation
    ## generates warnings that are not needed. And by decree (even in the
    ## C++) this thing has a fitness of 1, a mutator effect of 1 since
    ## there are no terms.
    
    if(addwt)
        df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
                               stringsAsFactors = FALSE), df)
    if(calledBy_ == "evalGenotype") { 
        if(prodNeg)
            colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
        class(df) <- c("evalAllGenotypes", class(df))
    } else if (calledBy_ == "evalGenotypeMut") {
        colnames(df)[match("Fitness", colnames(df))] <- "MutatorFactor"
        class(df) <- c("evalAllGenotypesMut", class(df))
    }
    return(df)
}

evalAllGenotypes <- function(fitnessEffects, order = FALSE, max = 256,
                             addwt = FALSE,
                             model = "") {
    ## Must deal with objects from previous, pre flfast, modifications
    if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
        fitnessEffects$fitnessLandscape_df <- data.frame()
        fitnessEffects$fitnessLandscape_gene_id <- data.frame()
    }

    if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
        (nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
        warning("Bozic model passing a fitness landscape will not work",
                    " for now.")
    }
    evalAllGenotypesORMut(
        fmEffects = fitnessEffects,
        order = order,
        max = max,
        addwt = addwt,
        model = model,
        calledBy_= "evalGenotype"
    )
}

generateAllGenotypes <- function(fitnessEffects, order = TRUE, max = 256) {
    if(order)
        tot <- function(n) {sum(sapply(seq.int(n),
                                       function(x) choose(n, x) * factorial(x)))}
    else
        tot <- function(n) {2^n}
    
    nn <- nrow(fitnessEffects$geneModule) -1  +
        nrow(fitnessEffects$long.geneNoInt) +
        nrow(fitnessEffects$fitnessLandscape_gene_id)
    
    tnn <- tot(nn)
    if(tnn > max) {
        m <- paste("There are", tnn, "genotypes.")
        m <- paste(m, "This is larger than max.")
        m <- paste(m, "Adjust max and rerun if you want")
        stop(m)
    }
    ## With mutator, the ids of genes need not go from 1:n
    vid <- allNamedGenes(fitnessEffects)$GeneNumID
    if(order) {
        f1 <- function(n) {
            lapply(seq.int(n), function(x) permutations(n = n, r = x, v = vid))
        }
    } else {
        f1 <- function(n) {
            lapply(seq.int(n), function(x) combinations(n = n, r = x, v = vid))}
        
    }
    genotNums <- f1(nn)
    list.of.vectors <- function(y) {
        ## there's got to be a simpler way
        lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
               function(m) m[[1]])
    }
    genotNums <- list.of.vectors(genotNums)
    names <- c(fitnessEffects$geneModule$Gene[-1],
               fitnessEffects$long.geneNoInt$Gene,
               fitnessEffects$fitnessLandscape_gene_id$Gene)
    
    genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
                                function(z)
                                    paste(z,
                                          collapse = if(order){" > "} else {", "} )))
    ## This ain't the best, as we repeatedly read and convert
    ## fitnessEffects.  If this were slow, prepare C++ function that takes
    ## vectors of genotypes and uses same fitnessEffects.
    return(list(genotNums = genotNums,
                genotNames = genotNames
                ))
}

evalAllGenotypesFitAndMut <- function(fitnessEffects, mutatorEffects,
                                   order = FALSE, max = 256,
                                   addwt = FALSE,
                                   model = "" ){
##                                   minimal = FALSE) {
    ## if(!minimal)
    allg <- generateAllGenotypes(fitnessEffects = fitnessEffects,
                                 order = order, max = max)
    ## else
        ## allg <- generateAllGenotypes_minimal(fitnessEffects = fitnessEffects,
        ##                                      max = max)
    
    if(model %in% c("Bozic", "bozic1", "bozic2") ) {
        prodNeg <- TRUE
    } else {
        prodNeg <- FALSE
    }


    ## Must deal with objects from previous, pre flfast, modifications
    if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
        fitnessEffects$fitnessLandscape_df <- data.frame()
        fitnessEffects$fitnessLandscape_gene_id <- data.frame()
    }
    if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
        (nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
        warning("Bozic model passing a fitness landscape will not work",
                    " for now.")
    } 
    
    full2mutator_ <- matchGeneIDs(mutatorEffects,
                                  fitnessEffects)$Reduced
    allf <- t(vapply(allg$genotNums,
                   function(x) evalRGenotypeAndMut(x,
                                                   rFE = fitnessEffects,
                                                   muEF= mutatorEffects,
                                                   full2mutator_ = full2mutator_,
                                                   verbose = FALSE,
                                                   prodNeg = prodNeg),
                   c(1.1, 2.2)))

    df <- data.frame(Genotype = allg$genotNames, Fitness = allf[, 1],
                     MutatorFactor = allf[, 2],
                     stringsAsFactors = FALSE)
    if(addwt)
        df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
                               MutatorFactor = 1,
                               stringsAsFactors = FALSE), df)
    if(prodNeg)
        colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
    class(df) <- c("evalAllGenotypesFitAndMut", class(df))
    return(df)
}






## evalAllGenotypes <- function(fitnessEffects, order = TRUE, max = 256,
##                              addwt = FALSE,
##                              model = "") {
    
##     if(order)
##         tot <- function(n) {sum(sapply(seq.int(n),
##                                        function(x) choose(n, x) * factorial(x)))}
##     else
##         tot <- function(n) {2^n}
##     nn <- nrow(fitnessEffects$geneModule) -1  + nrow(fitnessEffects$long.geneNoInt)
##     tnn <- tot(nn)
##     if(tnn > max) {
##         m <- paste("There are", tnn, "genotypes.")
##         m <- paste(m, "This is larger than max.")
##         m <- paste(m, "Adjust max and rerun if you want")
##         stop(m)
##     }
##     if(order) {
##         f1 <- function(n) {
##             lapply(seq.int(n), function(x) permutations(n = n, r = x))
##         }
##     } else {
##         f1 <- function(n) {
##             lapply(seq.int(n), function(x) combinations(n = n, r = x))}
        
##     }
##     genotNums <- f1(nn)
##     list.of.vectors <- function(y) {
##         ## there's got to be a simpler way
##         lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
##                function(m) m[[1]])
##     }
##     genotNums <- list.of.vectors(genotNums)
##     names <- c(fitnessEffects$geneModule$Gene[-1],
##                fitnessEffects$long.geneNoInt$Gene)

##     genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
##                          function(z)
##                              paste(z,
##                                    collapse = if(order){" > "} else {", "} )))
##     ## This ain't the best, as we repeatedly read and convert
##     ## fitnessEffects.  If this were slow, prepare C++ function that takes
##     ## vectors of genotypes and uses same fitnessEffects.
##     if(model %in% c("Bozic", "bozic1", "bozic2") )
##         prodNeg <- TRUE
##     else
##         prodNeg <- FALSE
##     allf <- vapply(genotNums,
##                    function(x) evalRGenotype(x, fitnessEffects, FALSE, prodNeg),
##                    1.1)
##     df <- data.frame(Genotype = genotNames, Fitness = allf,
##                      stringsAsFactors = FALSE)
##     if(addwt)
##         df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
##                                stringsAsFactors = FALSE), df)
##     if(prodNeg)
##         colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
##     return(df)
## }

fitnessEffectsToIgraph <- function(rT, epistasis, orderEffects) {

    df0 <- df1 <- df2 <- data.frame()
             
    if(!is.null(rT)) {
        df0 <- rT[, c("parent", "child", "typeDep")]
    }
    if(!is.null(epistasis)) {
        df1 <- to.pairs.modules(epistasis, ":")
    }
    if(!is.null(orderEffects)) {
        df2 <- to.pairs.modules(orderEffects, ">")
    }
    df <- rbind(df0, df1, df2)
    ## for special case of simple epi effects
    if(nrow(df) == 0) return(NULL)
    
    g1 <- graph.data.frame(df)
    
    E(g1)$color <- "black"
    E(g1)$color[E(g1)$typeDep == "SM"] <- "blue"
    E(g1)$color[E(g1)$typeDep == "XMPN"] <- "red"
    E(g1)$color[E(g1)$typeDep == "epistasis"] <- "orange"
    E(g1)$color[E(g1)$typeDep == "orderEffect"] <- "orange"
    E(g1)$lty <- 1
    E(g1)$lty[E(g1)$typeDep == "epistasis"] <- 2
    E(g1)$lty[E(g1)$typeDep == "orderEffect"] <- 3
    E(g1)$arrow.mode <- 2
    E(g1)$arrow.mode[E(g1)$typeDep == "epistasis"] <- 0
    E(g1)$arrow.mode[E(g1)$typeDep == "orderEffect"] <- 2
    return(g1)
}


plot.fitnessEffects <- function(x, type = "graphNEL",
                                layout = NULL,
                                expandModules = FALSE,
                                autofit = FALSE,
                                scale_char = ifelse(type == "graphNEL", 1/10, 5),
                                return_g = FALSE,
                                lwdf = 1, ## multiplier for lwd for graphNEL
                                ...) {
    ## some other layouts I find OK
    ## layout.circle
    ## layout.reingold.tilford if really a tree
    ## o.w. it will use the default
    g <- x$graph
    if(is.null(g))
        stop("This fitnessEffects object can not be ploted this way.",
             " It is probably one with fitness landscape specification, ",
             " so you might want to plot the fitness landscape instead.")
    if(type == "igraph") {
        if(expandModules && (!x$gMOneToOne)) {
            ## vlabels <- fe$geneToModule[vertex.attributes(g)$name]
            vlabels <- x$geneToModule[V(g)$name]
            V(g)$label <- vlabels
        }
        if(autofit) {
            plot(0, type = "n", axes = FALSE, ann = FALSE)
            ## ideas from http://stackoverflow.com/questions/14472079/match-vertex-size-to-label-size-in-igraph
            ## vsize <- (strwidth(V(g)$label) + strwidth("oo")) * 200
            ## but this is a kludge.
            vsize <- (nchar(V(g)$label) + 1) * scale_char
            plot.igraph(g, vertex.size = vsize, vertex.shape = "rectangle",
                        layout = layout)
        } else {
            plot.igraph(g, layout = layout)
        }
        if(return_g) return(g)
    }
    else if (type == "graphNEL") {
        g1 <- igraph.to.graphNEL(g)
        c1 <- unlist(lapply(edgeData(g1), function(x) x$color))
        names(c1) <- sub("|", "~", names(c1), fixed = TRUE)
        s1 <- unlist(lapply(edgeData(g1), function(x) x$lty))
        names(s1) <- names(c1)
        a1 <- unlist(lapply(edgeData(g1), function(x) max(x$arrow.mode - 1, 0)))
        names(a1) <- names(c1)
        lwd <- s1
        lwd[lwd == 2] <- 2 ## o.w. too thin
        lwd[lwd == 3] <- 2 ## o.w. too thin
        lwd <- lwd * lwdf
        nAttrs <- list()
        if(expandModules && (!x$gMOneToOne)) {
            nnodes <- x$geneToModule[nodes(g1)]
            names(nnodes) <- nodes(g1)
            nAttrs$label <- nnodes
        }
        if(autofit) {
            nAttrs$width <- (nchar(nAttrs$label) + 1) * scale_char
            names(nAttrs$width) <- names(nAttrs$label)
            plot(g1, edgeAttrs = list(arrowsize = a1, lty = s1, lwd = lwd,
                         color = c1), attrs=list(node=list(shape = "rectangle")),
                 nodeAttrs = nAttrs)
            
        } else {
            plot(g1, edgeAttrs = list(arrowsize = a1, lty = s1, lwd = lwd,
                         color = c1),
                 nodeAttrs = nAttrs)
        }
        if(return_g) return(g1)
    } else {
        stop("plot type not recognized")
    }
}

## plot.fitnessEffects <- plotFitnessEffects

## FIXME: in the help: say we cannot properly show 3- or higher order
## interactions.



nr_oncoSimul.internal <- function(rFE, 
                                  birth, 
                                  death,
                                  mu,
                                  initSize,
                                  sampleEvery,
                                  detectionSize,
                                  finalTime,
                                  initSize_species,
                                  initSize_iter,
                                  seed,
                                  verbosity,
                                  speciesFS,
                                  ratioForce,
                                  typeFitness,
                                  max.memory,
                                  mutationPropGrowth,
                                  initMutant,
                                  max.wall.time,
                                  keepEvery,
                                  K,
                                  detectionDrivers,
                                  onlyCancer,
                                  errorHitWallTime,
                                  max.num.tries,
                                  errorHitMaxTries,
                                  minDetectDrvCloneSz,
                                  extraTime,
                                  keepPhylog,
                                  detectionProb,
                                  AND_DrvProbExit,
                                  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",
                   "allFitnessEffects"))

    if(countGenesFe(rFE) < 2) {
        stop("There must be at least two genes (loci) in the fitness effects.",
             "If you only care about a degenerate case with just one,",
             "you can enter a second gene (locus)",
             "with fitness effect of zero.")
    }
    
    if( (length(mu) == 1) && !(is.null(names(mu)))) {
        stop("A length 1 mutation, but named. ",
             "This is ambiguous. ",
             "If you want per-gene mutation rates, each gene",
             "must have its entry in the mu vector. ",
             "(And regardless of per-gene mutation rates ",
             " there must be at least two gene/loci).")
    }

    ## Must deal with objects from previous, pre flfast, modifications
    ## Could move this way down to the bottom, right before
    ## .Call
    if(!exists("fitnessLandscape_gene_id", where = rFE)) {
        rFE$fitnessLandscape_df <- data.frame()
        rFE$fitnessLandscape_gene_id <- data.frame()
    }
    
    namedGenes <- allNamedGenes(rFE)

    if( length(mu) > 1) {
        if(is.null(names(mu)))
            stop("When using per-gene mutation rates the ",
                 "mu vector must be named ",
                 "(and if you have noIntGenes, those must have names).")
        if(length(mu) != countGenesFe(rFE))
            stop("When using per-gene mutation rates, ",
                 "there must be the same number of genes in the ",
                 "mu vector and the fitness effects object.")
        if(!identical(sort(namedGenes$Gene),
                      sort(names(mu))))
            stop("When using per-gene mutation rates, ",
                 "names of genes must match those in the ",
                 "restriction table.")
        mu <- mu[order(match(names(mu), namedGenes$Gene))]
        ## Hyperparanoid check. Should never, ever, happen.
        if(!identical(names(mu), namedGenes$Gene))
            stop("Names of re-ordered mu do not match names of genes")
        minmu <- 1e-40
        if(any(mu <= minmu))
            stop(paste("At least one per-gene mutation rate is negative",
                       "or less than", minmu,". Remember that the per-base",
                       "mutation rate in the human genome is about 1e-11 to 1e-9."))
    }
    if(!is.null(initMutant)) {
        initMutantString <- initMutant
       if(length(grep(">", initMutant))) {
            initMutant <- nice.vector.eo(initMutant, ">")
        } else if(length(grep(",", initMutant))) {
            initMutant <- nice.vector.eo(initMutant, ",")
        }
        initMutant <- getGeneIDNum(rFE$geneModule,
                                   rFE$long.geneNoInt,
                                   rFE$fitnessLandscape_gene_id,
                                   initMutant,
                             FALSE)
       if(length(initMutant) >= countGenesFe(rFE)) {
        stop("For initMutant you passed as many, or more genes, mutated ",
             "than the number of genes in the genotype (fitness effects).")
    }
       
    } else {
        initMutant <- vector(mode = "integer")
        initMutantString <- ""
    }
    ## these are never user options
    ## if(initSize_species < 10) {
    ##     warning("initSize_species too small?")
    ## }
    ## if(initSize_iter < 100) {
    ##     warning("initSize_iter too small?")
    ## }

    if(typeFitness %in% c("bozic1", "bozic2")) {
        if(nrow(rFE$fitnessLandscape_df) > 0)
            warning("Bozic model passing a fitness landscape will not work",
                    " for now.")
        ## FIXME: bozic and fitness landscape
        ## the issue is that in the C++ code we directly do
        ## s = birth rate - 1
        ## but we would need something different
        ## Can be done going through epistasis, etc.
        thesh <- unlist(lapply(rFE$long.rt, function(x) x$sh))
        ## thes <- unlist(lapply(rFE$long.rt, function(x) x$s))
        thes <- unlist(c(lapply(rFE$long.rt, function(x) x$s),
                         lapply(rFE$long.epistasis, function(x) x$s),
                         lapply(rFE$long.orderEffects, function(x) x$s),
                         rFE$long.geneNoInt$s))
        
        if(any(thes > 1 )) {
            m <- paste("You are using a Bozic model with",
                       "the new restriction specification, and you have",
                       "at least one s > 1.",
                       "But that is the same as setting that to 1:",
                       "we obviously cannot allow negative death rates,",
                       "nor problems derived from multiplying odd or even",
                       "numbers of negative numbers.",
                       "You can set those > 1 to 1, but then you will",
                       "eventually get the simulations to abort when the",
                       "death rate becomes 0.")
            stop(m)
        }
        if(any( (thesh <= -1) & (thesh > -Inf))) {
            m <- paste("You are using a Bozic model with",
                       "the new restriction specification, and you have",
                       "at least one sh <= -1. Maybe you mean -Inf?")
            warning(m)
        }
        if(any(thes == 1 )) {
            m <- paste("You are using a Bozic model with",
                       "the new restriction specification, and you have",
                       "at least one s of 1. If that gene is mutated, ",
                       "this will lead to a death rate of 0",
                       "and the simulations will abort when you get a non finite",
                       "value.")
            warning(m)
        }
    }

    if(!is.null(MMUEF)) {
        if(!inherits(MMUEF, "mutatorEffects"))
            stop("muEF must be a mutatorEffects object")
        full2mutator_ <- matchGeneIDs(MMUEF, rFE)
        if(any(is.na(full2mutator_$Full)))
            stop("Genes in mutatorEffects not present in fitnessEffects")
        if(any(is.na(full2mutator_)))
            stop("full2mutator failed for unknown reasons")
        full2mutator_ <- full2mutator_$Reduced
    } else {
        full2mutator_ <- vector(mode = "numeric", length = 0)
        ## muEF <- emptyFitnessEffects()
    }

    dpr <- detectionProbCheckParse(detectionProb, initSize, verbosity)
    ## if( !is.null(cPDetect) && (sum(!is.null(p2), !is.null(n2)) >= 1 ))
    ##     stop("Specify only cPDetect xor both of p2 and n2")
    ## if( (is.null(p2) + is.null(n2)) == 1 )
    ##     stop("If you pass one of n2 or p2, you must also pass the other. ",
    ##          "Otherwise, we would not know what to do.")
    ## stopifnot(PDBaseline >= 0)
    ## stopifnot(n2 > PDBaseline)
    ## stopifnot(p2 < 1)
    ## stopifnot(p2 > 0)
    ## if( is.null(cPDetect) ) cPDetect <- -9
    ## if( is.null(p2)) p2 <- 9
    ## 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"]
        ## 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 = ",")
        ulf <- unlist(fixation_b)
        if(any(ulf == ">"))
            stop("Order effects not allowed (yet?) in fixation.")
        ulfg <- ng[ulf, 1]
        if(any(is.na(ulfg)))
            stop(paste("The 'fixation' list contains genes that are not present",
                       " 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()
    }

    return(c(
        nr_BNB_Algo5(rFE = rFE,
                     mu_ = mu,
                     death = death,
                     initSize = initSize,
                     sampleEvery = sampleEvery,
                     detectionSize = detectionSize,
                     finalTime = finalTime,
                     initSp = initSize_species,
                     initIt = initSize_iter,
                     seed = seed,
                     verbosity = verbosity,
                     speciesFS = speciesFS,
                     ratioForce = ratioForce,
                     typeFitness_ = typeFitness,
                     maxram = max.memory,
                     mutationPropGrowth = as.integer(mutationPropGrowth),
                     initMutant_ = initMutant, 
                     maxWallTime = max.wall.time,
                     keepEvery = keepEvery,
                     K = K,
                     detectionDrivers = detectionDrivers,
                     onlyCancer = onlyCancer,
                     errorHitWallTime = errorHitWallTime,
                     maxNumTries = max.num.tries,
                     errorHitMaxTries = errorHitMaxTries,
                     minDetectDrvCloneSz = minDetectDrvCloneSz,
                     extraTime = extraTime,
                     keepPhylog = keepPhylog,
                     MMUEF = MMUEF,
                     full2mutator_ = full2mutator_,
                     n2 = dpr["n2"],
                     p2 = dpr["p2"],
                     PDBaseline = dpr["PDBaseline"],
                     cPDetect_i= dpr["cPDetect"],
                     checkSizePEvery = dpr["checkSizePEvery"],
                     AND_DrvProbExit = AND_DrvProbExit,
                     fixation_list = fixation_list),
        Drivers = list(rFE$drv), ## but when doing pops, these will be repeated
        geneNames = list(names(getNamesID(rFE))),
        InitMutant = initMutantString
    ))
}


countGenesFe <- function(fe) {
    ## recall geneModule has Root always
    ## We want to be able to use objects that did not have
    ## the fitness landscape component
    if(exists("fitnessLandscape_gene_id", where = fe))
        return(nrow(fe$geneModule) + nrow(fe$long.geneNoInt) +
            nrow(fe$fitnessLandscape_gene_id) - 1)
    else
        return(nrow(fe$geneModule) + nrow(fe$long.geneNoInt) - 1)
}

allNamedGenes <- function(fe){
    ## Returns a data frame with genes and their names and verifies all
    ## genes have names.

    ## Root should always be first, but just in case
    ## avoid assuming it

    ## This does is not a good idea as it assumes the user did not use
    ## "keepInput = FALSE".
    ## lni <- length(fe$noIntGenes)
    ## ## FIXME:test
    ## if((lni > 0) &&
    ##        (is.null(names(fe$noIntGenes))))
    ##         stop("When using per-gene mutation rates the ",
    ##              "no interaction genes must be named ",
    ##              "(i.e., the noIntGenes vector must have names).")

    ## accommodate objects w/o fitnessLandscape
    if(!is.null(fe$fitnessLandscape) && nrow(fe$fitnessLandscape)) {
        gn <-
            gtools::mixedsort(
                        colnames(fe$fitnessLandscape)[-ncol(fe$fitnessLandscape)])
        v1 <- data.frame(Gene = gn, GeneNumID = seq.int(length(gn)),
                        stringsAsFactors = FALSE)
    } else {
        v1 <- fe$geneModule[, c("Gene", "GeneNumID")]
        if(nrow(fe$long.geneNoInt)) {
            v1 <- rbind(v1,
                        fe$long.geneNoInt[, c("Gene", "GeneNumID")])
        }
        if(any(v1[, "Gene"] == "Root"))
            v1 <- v1[-which(v1[, "Gene"] == "Root"), ]
    }
    rownames(v1) <- NULL
    if(any(v1$Gene == "WT")) {
        warning("A gene is named WT. You can expect problems ",
                "because we use WT to denote the wildtype ",
                "genotype. You might want to change it.")
    }
    return(v1)
}


get.gene.counts <- function(x) {
                                        # , timeSample = "last",
                                        # typeSample = "whole") {
    ## From get.mut.vector. Used for now just for testing
    timeSample <- "last"
    the.time <- get.the.time.for.sample(x, timeSample, NULL)
    if(the.time < 0) {
        counts <- rep(0, length(x$geneNames))
        names(counts) <- x$geneNames
        freq <- rep(NA, length(x$geneNames))
        names(freq) <- x$geneNames
        return(list(counts = counts,
                    freq = freq,
                    popSize = 0))
    } 
    pop <- x$pops.by.time[the.time, -1]
    if(all(pop == 0)) {
        stop("You found a bug: this should never happen")
    }
    ## if(typeSample %in% c("wholeTumor", "whole")) {
    popSize <- x$PerSampleStats[the.time, 1]
    counts <- as.vector(tcrossprod(pop, x$Genotypes))
    names(counts) <- x$geneNames    
    return(list(counts = counts,
                freq = counts/popSize,
                popSize = popSize))
    ## return( (tcrossprod(pop,
    ##                     x$Genotypes)/popSize) )
    ## } else if (typeSample %in%  c("singleCell", "single")) {

    ##       return(x$Genotypes[, sample(seq_along(pop), 1, prob = pop)])
    ##   } else {
    ##         stop("Unknown typeSample option")
    ##     }
}



geneCounts <- function(x) {
    if(inherits(x, "oncosimulpop")) {
        return( do.call(rbind,
                        lapply(x,
                               function(z)
                               get.gene.counts(z)$counts))
               )
    } else {
        return( get.gene.counts(x)$counts)
    }
}





## geneNumIDReset <- function(x, ref){
##     ## Set GeneNumID of a fitnessEffect object, x, using ref as the
##     ## reference fitness effect object.
##     ## Check also if all in x are in ref.

##     gg <- allNamedGenes(ref)
##     gnid <- gg$GeneNumID
##     names(gnid) <- gg$Gene
##     ## FIXME: this later and conditional on what is in thee
##     gnid <- c("Root" = 0, gnid)
    
##     if(!all(x$geneModule$Gene %in% names(gnid) ))
##         stop("Some genes not in reference fitnessEffects (rebasing geneModule)")
##     x$geneModule$GeneNumID <- gnid[geneModule$Gene]

##     ## and then go over all the lists in the x object. 
    
##     if(nrow(x$long.geneNoInt)) {
##         ## now, mapping for the noInt if this is mutator
##         if(!all(x$long.geneNoInt$Gene %in% names(gnid) ))
##             stop("Some genes not in reference fitnessEffects (rebasing geneNoInt)")
##         x$long.geneNoInt$GeneNumID <- gnid[long.geneNoInt$Gene]
##     }
## }


matchGeneIDs <- function(x, refFE) {
    n1 <- allNamedGenes(x)
    n2 <- allNamedGenes(refFE)
    colnames(n1)[2] <- "Reduced"
    colnames(n2)[2] <- "Full"
    ## Non standard evaluation thing and a note being generated in check.
    ## See also http://stackoverflow.com/questions/33299845/when-writing-functions-for-an-r-package-should-i-avoid-non-standard-evaluation
    ## But does not work well with replace. So use the NULL trick
    Reduced <- NULL
    dplyr::full_join(n2, n1, by = "Gene") %>%
        mutate(Reduced = replace(Reduced, is.na(Reduced), -9))
}

    
detectionProbCheckParse <- function(x, initSize, verbosity) {
    default_p2 <- 0.1
    default_n2 <- 2 * initSize
    default_PDBaseline <- 1.2 * initSize
    default_checkSizePEvery <- 20
    ## No default cPDetect. That is done from p2 and n2 in C++.
    if(is_null_na(x)) {
    ## if((length(x) == 1) && (is.na(x))) {
        y <- vector()
        y["cPDetect"] <- -9
        y["p2"] <- 9
        y["n2"] <- -9
        y["PDBaseline"] <- -9
        y["checkSizePEvery"] <- Inf
        return(y)
    } else if((length(x) == 1) && (x == "default")) {
        ## Populate with defaults
        y <- vector()
        y["p2"] <- default_p2
        y["n2"] <- default_n2
        y["PDBaseline"] <- default_PDBaseline
        y["checkSizePEvery"] <- default_checkSizePEvery
        x <- y
    }

    if(length(x) >= 1) {
        if( !all(names(x) %in% c("cPDetect", "p2", "n2", "PDBaseline", "checkSizePEvery")))
            stop("Names of some components of detectionProb are not recognized.",
                 " Check for possible typos.")
    }
   
    ## This ain't conditional. If not given, always check
    if( !is_null_na(x["cPDetect"]) && (sum(!is_null_na(x["p2"]), !is_null_na(x["n2"])) >= 1 ))
        stop("Specify only cPDetect xor both of p2 and n2")
    if( (is_null_na(x["p2"]) + is_null_na(x["n2"])) == 1 )
        stop("If you pass one of n2 or p2, you must also pass the other. ",
             "Otherwise, we would not know what to do.")

    if(is_null_na(x["PDBaseline"])) {
        x["PDBaseline"] <- default_PDBaseline
        if(verbosity > -1)
            message("\n  PDBaseline set to default value of ", default_PDBaseline, "\n")
        }
    if(is_null_na(x["checkSizePEvery"])) {
        x["checkSizePEvery"] <- default_checkSizePEvery
        if(verbosity > -1)
            message("\n  checkSizePEvery set to default value of ",
                default_checkSizePEvery, "\n")
        }

    ## If we get here, we checked consistency of whether cPDetect or p2/n2.
    ## Fill up with defaults the missing values

    if(is_null_na(x["cPDetect"])) {
        if(is_null_na(x["p2"])) {
            if(!is_null_na(x["n2"])) stop("Eh? no p2 but n2? Bug")
            x["n2"] <- default_n2
            x["p2"] <- default_p2
            if(verbosity > -1)
                message("\n  n2, p2 set to default values of ",
                    default_n2, ", ", default_p2, "\n")
        }
    }
    
    if(x["checkSizePEvery"] <= 0)
        stop("checkSizePEvery <= 0")
    if(x["PDBaseline"] <= 0)
        stop("PDBaseline <= 0")
    
    if(!is_null_na(x["n2"])) {
        if(x["n2"] <= x["PDBaseline"])
            stop("n2 <= PDBaseline")
        if(x["p2"] >= 1) stop("p2 >= 1")
        if(x["p2"] <= 0) stop("p2 <= 0")
        x["cPDetect"] <- -9
    } else { ## only if x["cPDetect"] is not NA
        if(is_null_na(x["cPDetect"])) stop("eh? you found a bug")## paranoia
        x["n2"] <- -9
        x["p2"] <- -9
        if(verbosity > -1)
            message("\n Using user-specified cPDetect as n2, p2 not given.\n")
    }
    return(x)
}

sampledGenotypes <- function(y, genes = NULL) {
    ## From a popSample object, or a matrix for that matter,
    ## show the sampled genotypes and their frequencies
    if(!is.null(genes)) {
        cols <- which(colnames(y) %in% genes )
        y <- y[, cols]
    }
    ## 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, genotlabel),
        useNA = "ifany"),
        stringsAsFactors = TRUE)

    gn <- as.character(df[, 1])
    gn[gn == ""] <- "WT"
    df <- data.frame(Genotype = gn, Freq = df[, 2], stringsAsFactors = FALSE)
    attributes(df)$ShannonI <- shannonI(na.omit(df)$Freq)
    class(df) <- c("sampledGenotypes", class(df))
    return(df)
}

print.sampledGenotypes <- function(x, ...) {
    print.data.frame(x, ...)
    cat("\n Shannon's diversity (entropy) of sampled genotypes: ")
    cat(attributes(x)$ShannonI, "\n")
}


list_g_matches_fixed <- function(x, y) {
    ## Internal function, for testing the fixation output.
    ## x and y are vectors

    ## x is the set of output genotypes, y the set of fixed
    ## genotypes/subset of genotypes.

    ## Yes, this function has tests too in test.fixation.R
    
    ## All genotypes in x satisfy that they are supersets of at least one
    ## in y? That is true if, for every element in x, at least one y in
    ## that x.

    if(is.list(y)) y <- unlist(y)
    
    y.nice <- lapply(y, nice.vector.eo, sep = ",")
    x.nice <- lapply(x, nice.vector.eo, sep = ",")

    fu <- function(u, y.nice)
        any(unlist(lapply(y.nice, function(z) all(z %in% u))))

    return(all(unlist(lapply(x.nice, fu, y.nice))))

}



## emptyFitnessEffects <- function() {
##     list(long.rt = list(),
##          long.epistasis = list(),
##          long.orderEffects = list(),
##          long.geneNoInt = list(),
##          geneModule = list(),
##          gMOneToOne = TRUE,
##          geneToModule = list(),
##          graph = NULL,
##          drv = vector(mode = "integer", length = 0)
##          )
## }



### Later, for all the effects, we will do some kind of dplyr match?

### a, b, c, in fitness, only a, c in mut.
### fitness table for a,b,c
### each row name transformed removing b (so leaving only present)
### each row transformed matched to row in mut table.

## t1 <- data.frame(v1 = c("a,b", "a,c", "b"), v2 = c("b", "c", "b"), v3 = 1:3, stringsAsFactors = FALSE)
## t2 <- data.frame(v2 = c("b", "c"), v4 = c(11, 12), stringsAsFactors = FALSE)
## full_join(t1, t2, by = "v2")



## FIXME

## new exit code
## check:
## baseline < n2
## p2 < 1

## baseline defaults to init size + .1
## use c < -9, n2 < -9, p2 < -9 for no values
## and check one of c or p2 and n2 are valid if using exit