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

                                  min = 1e-10,

                                  max = 1e-9)

         }

-        m <- cbind(m, Fitness = fi)

+        m <- cbind(m, Birth = fi)

         if(!is_null_na(min_accessible_genotypes)) {

             ## num_accessible_genotypes <- count_accessible_g(m, accessible_th)

             ## Recall accessibleGenotypes includes the wt, if accessible.

@@ -232,14 +232,28 @@ rfitness <- function(g, c= 0.5,

             }

         } else { ## a length-3 scale

             if(scale[2] > scale[3]) warning("In scale, minimum fitness > wildtype")

+            if(scale[1] < scale[3]) warning("In scale, maximum fitness < wildtype")

             fiwt <- fi[1]

-            prod_above <- (scale[1] - scale[3]) / (max(fi) - fiwt)

-            prod_below <- (scale[3] - scale[2]) / (fiwt - min(fi))

-            fi_above <- which(fi >= fiwt)

-            fi_below <- which(fi < fiwt)

-            fi[fi_above] <- ((fi[fi_above] - fiwt) * prod_above) + scale[3]

-            fi[fi_below] <- ((fi[fi_below] - fiwt) * prod_below) + scale[3]

-            fi[1] <- scale[3]

+            new_fi <- rep(NA, length(fi))

+            mode(new_fi) <- "numeric"

+            ## If WT is min or max, there are no below or above

+            if(max(fi) == fiwt)

+                warning("WT has maximum fitness. Range will be from scale[2] to scale[3]")

+            if(min(fi) == fiwt)

+                warning("WT has minimum fitness. Range will be from scale[3] to scale[1]")

+            if(max(fi) > fiwt) {

+                prod_above <- (scale[1] - scale[3]) / (max(fi) - fiwt)

+                fi_above <- which(fi >= fiwt)

+                new_fi[fi_above]  <- ((fi[fi_above] - fiwt) * prod_above) + scale[3]

+            }

+            if(min(fi) < fiwt) {

+                prod_below <- (scale[3] - scale[2]) / (fiwt - min(fi))

+                fi_below <- which(fi < fiwt)

+                new_fi[fi_below] <- ((fi[fi_below] - fiwt) * prod_below) + scale[3]

+            }

+            new_fi[1] <- scale[3]

+            fi <- new_fi

+            rm(new_fi)

         }

         if(log) {

@@ -202,48 +202,66 @@ rfitness <- function(g, c= 0.5,

             ## rm(m1)

             ## rm(fl1)

         }

-

-        if(!is.null(scale)) {

-            fi <- (fi - min(fi))/(max(fi) - min(fi))

-            fi <- scale[1] + fi * (scale[2] - scale[1])

-        }

-        if(wt_is_1 == "divide") {

-            ## We need to shift to avoid ratios involving negative numbers and

-            ## we need to avoid having any fitness at 0, specially the wt.  If

-            ## any negative value, add the min, and shift by the magnitude of

-            ## the min to avoid any at 0.

-

-            ## If you use scale and wt_is_1, this will move the scale. It is

-            ## not possible to obtain a linear transformation that will keep

-            ## the min and max of the scale, and wt at 1.

-            min_fi <- min(fi)

-            if(min_fi < 0)

-                fi <- fi + 2 * abs(min(fi))

-            fi <- fi/fi[1]

-        } else if (wt_is_1 == "subtract") {

-            fi <- fi - fi[1] + 1.0

-        } else if(wt_is_1 == "force") {

-            fi[1] <- 1.0

+        if(!(length(scale) == 3)) {

             if(!is.null(scale)) {

-                if( (1 > scale[2]) || (1 < scale[1]))

-                    warning("Using wt_is_1 = force and scale, but scale does ",

-                            "not include 1")

+                fi <- (fi - min(fi))/(max(fi) - min(fi))

+                fi <- scale[1] + fi * (scale[2] - scale[1])

+            }

+            if(wt_is_1 == "divide") {

+                ## We need to shift to avoid ratios involving negative numbers and

+                ## we need to avoid having any fitness at 0, specially the wt.  If

+                ## any negative value, add the min, and shift by the magnitude of

+                ## the min to avoid any at 0.

+

+                ## If you use scale and wt_is_1, this will move the scale. It is

+                ## not possible to obtain a linear transformation that will keep

+                ## the min and max of the scale, and wt at 1.

+                min_fi <- min(fi)

+                if(min_fi < 0)

+                    fi <- fi + 2 * abs(min(fi))

+                fi <- fi/fi[1]

+            } else if (wt_is_1 == "subtract") {

+                fi <- fi - fi[1] + 1.0

+            } else if(wt_is_1 == "force") {

+                fi[1] <- 1.0

+                if(!is.null(scale)) {

+                    if( (1 > scale[2]) || (1 < scale[1]))

+                        warning("Using wt_is_1 = force and scale, but scale does ",

+                                "not include 1")

+                }

             }

+        } else { ## a length-3 scale

+            if(scale[2] > scale[3]) warning("In scale, minimum fitness > wildtype")

+            fiwt <- fi[1]

+            prod_above <- (scale[1] - scale[3]) / (max(fi) - fiwt)

+            prod_below <- (scale[3] - scale[2]) / (fiwt - min(fi))

+            fi_above <- which(fi >= fiwt)

+            fi_below <- which(fi < fiwt)

+            fi[fi_above] <- ((fi[fi_above] - fiwt) * prod_above) + scale[3]

+            fi[fi_below] <- ((fi[fi_below] - fiwt) * prod_below) + scale[3]

+            fi[1] <- scale[3]

         }

+        

         if(log) {

+            if(length(scale == 3)) {

+                wt_is_1 <- "no"

+                message("You passed a three-element scale argument. ",

+                        "Setting wt_is_1 to no ",

+                        "to avoid modifying your requested value for WT.")

+            }

             ## If you want logs, you certainly do not want

             ## the log of a negative number

             fi[fi < 0] <- 0

             if(wt_is_1 == "no") {

                 fi <- log(fi)

-            } else {

-                ## by decree, fitness of wt is 1. So shift everything

-                fi <- log(fi) + 1

-            }

-            ## former expression, but it was more confusing

+                } else {

+                    ## by decree, fitness of wt is 1. So shift everything

+                    fi <- log(fi) + 1

+                }

+                ## former expression, but it was more confusing

             ## fi <- log(fi/fi[1]) + 1

-            

         }

+        

         if(truncate_at_0) {

             ## yes, truncate but add noise to prevent identical

             fi[fi <= 0] <- runif(sum(fi <= 0),

@@ -279,110 +279,3 @@ rfitness <- function(g, c= 0.5,

-

-

-

-

-## rfitness <- function(g, c= 0.5,

-##                      sd = 1,

-##                      mu = 1,

-##                      reference = "random", ## "random", "max", or the vector,

-##                                      ## e.g., rep(1, g). If random, a

-##                                      ## random genotype is chosen as

-##                                      ## reference. If "max" this is rep(1, g)

-##                      scale = NULL, ## a two-element vector: min and max

-##                      wt_is_1 = c("subtract", "divide", "force", "no"),

-##                      ## wt_is_1 = TRUE, ## wt has fitness 1

-##                      log = FALSE, ## log only makes sense if all values >

-##                                  ## 0. scale with min > 0, and/or set

-##                                  ## wt_is_1 = divide

-##                      min_accessible_genotypes = NULL,

-##                      accessible_th = 0,

-##                      truncate_at_0 = TRUE) {

-##     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene

-##     ## and Crona, 2014. And this allows moving from HoC to purely additive

-##     ## changing c and sd.

-

-##     ## FIXME future: do this with order too?

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

-##     wt_is_1 = match.arg(wt_is_1)

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

-##     m <- generate_matrix_genotypes(g)

-##     done <- FALSE

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

-##     while(!done) {

-##         ## attempts <- attempts + 1

-##         f_r <- rnorm(nrow(m), mean = mu, sd = sd)

-##         if(inherits(reference, "character") && length(reference) == 1) {

-##             if(reference == "random") {

-##                 referenceI <- m[sample(nrow(m), 1), ]

-##             } else if(reference == "max") {

-##                 referenceI <- rep(1, g)

-##             } else if(reference == "random2") {

-##                 referenceI <- create_eq_ref(g)

-##             }

-##         } else {

-##             referenceI <- reference

-##             }

-##         d_reference <- apply(m, 1, function(x) sum(abs(x - referenceI)))

-##         f_det <- -c * d_reference

-##         ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona

-##         fi <- f_r + f_det

-

-##         if(!is.null(scale)) {

-##             fi <- (fi - min(fi))/(max(fi) - min(fi))

-##             fi <- scale[1] + fi * (scale[2] - scale[1])

-##         }

-##         if(wt_is_1 == "divide") {

-##             ## We need to shift to avoid ratios involving negative numbers and

-##             ## we need to avoid having any fitness at 0, specially the wt.  If

-##             ## any negative value, add the min, and shift by the magnitude of

-##             ## the min to avoid any at 0.

-

-##             ## If you use scale and wt_is_1, this will move the scale. It is

-##             ## not possible to obtain a linear transformation that will keep

-##             ## the min and max of the scale, and wt at 1.

-##             min_fi <- min(fi)

-##             if(min_fi < 0)

-##                 fi <- fi + 2 * abs(min(fi))

-##             fi <- fi/fi[1]

-##         } else if (wt_is_1 == "subtract") {

-##             fi <- fi - fi[1] + 1.0

-##         } else if(wt_is_1 == "force") {

-##             fi[1] <- 1.0

-##             if(!is.null(scale)) {

-##                 if( (1 > scale[2]) || (1 < scale[1]))

-##                     warning("Using wt_is_1 = force and scale, but scale does ",

-##                             "not include 1")

-##             }

-##         }

-##         if(truncate_at_0) {

-##             fi[fi <= 0] <- 1e-9

-##         }

-##         if(log) {

-##             fi <- log(fi/fi[1]) + 1

-##         }

-##         m <- cbind(m, Fitness = fi)

-##         if(!is_null_na(min_accessible_genotypes)) {

-##             ## num_accessible_genotypes <- count_accessible_g(m, accessible_th)

-##             ## Recall accessibleGenotypes includes the wt, if accessible.

-##             num_accessible_genotypes <- length(wrap_accessibleGenotypes(m, accessible_th)) - 1

-##             ## message("\n     num_accessible_genotypes = ", num_accessible_genotypes, "\n")

-##             if(num_accessible_genotypes >= min_accessible_genotypes) {

-##                 done <- TRUE

-##                 attributes(m) <- c(attributes(m),

-##                                    accessible_genotypes = num_accessible_genotypes,

-##                                    accessible_th = accessible_th)

-##             } else {

-##                 ## Cannot start again with a fitness column

-##                 m <- m[, -ncol(m), drop = FALSE]

-##             }

-##         } else {

-##             done <- TRUE

-##         }

-##     }

-##     ## message("\n number of attempts = ", attempts, "\n")

-##     class(m) <- c(class(m), "genotype_fitness_matrix")

-##     return(m)

-## }

-

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

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

+## Copyright 2013-2021 Ramon Diaz-Uriarte

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

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

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

                      truncate_at_0 = TRUE,

                      K = 1,

                      r = TRUE,

-                     model = c("RMF", "NK")) {

+                     i = 0, # Ising, cost incompatibility

+                     I = -1, # Ising, sd for "i"

+                     circular = FALSE, # Ising, circular arrangement

+                     e = 0, # Eggbox, +/- e

+                     E = -1, # Eggbox, noise on "e"

+                     H = -1, # HoC stdev

+                     s = 0.1, # mean multiplivative

+                     S = -1, # SD multiplicative

+                     d = 0, # disminishing/increasing for multiplicative

+                     o = 0, # mean optimum

+                     O = -1, # sd optimum

+                     p = 0, # mean production for non 0 allele (optimum)

+                     P = -1, # sd for p

+                    

+                     model = c("RMF", "Additive", 

+                               "NK", "Ising", "Eggbox", "Full")) {

     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene

     ## and Crona, 2014. And this allows moving from HoC to purely additive

     ## changing c and sd.

@@ -107,12 +122,46 @@ rfitness <- function(g, c= 0.5,

             f_det <- -c * d_reference

             ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona

             fi <- f_r + f_det

+        } else if (model == "Additive") {

+      ## get fitness effect for mutations in each gene

+      mutants <-rep(1,g)

+            ## FIXME: Why not just?

+            ## f_single_mut <- rnorm(g, mean = mu, sd = sd)

+            f_single_mut <- sapply(mutants, FUN = function(x) 

+                                            rnorm(x, mean = mu, sd = sd))

+      ## find which gene is mutated 

+      m2 <- m == 1

+      ## Sum the fitness effect of that mutation to generate a vector fi with

+      ## the fitness for each mutant condition

+      fi <- apply(m2, MARGIN = 1, FUN = function (x) sum(x*f_single_mut))

+      ## remove unnecessary variables

+      rm (f_single_mut, m2)

         } else if(model == "NK") {

             if(K >= g) stop("It makes no sense to have K >= g")

             argsnk <- paste0("-K ", K,

                              ifelse(r, " -r ", " "),

                              g, " 2")

             fl1 <- system2(fl_generate_binary(), args = argsnk, stdout = TRUE)[-1]

+        } else if (model == "Ising") {

+      argsIsing <- paste0("-i ", i, " -I ", I ,

+                          ifelse(circular, " -c ", " "),

+                          g, " 2")

+      fl1 <- system2(fl_generate_binary(), args = argsIsing, stdout = TRUE)[-1]

+    } else if (model == "Eggbox") {

+      argsEgg <- paste0("-e ", e, " -E ", E," ", g, " 2")

+      fl1 <- system2(fl_generate_binary(), args = argsEgg, stdout = TRUE)[-1]

+    } else if (model == "Full") {

+      if(K >= g) stop("It makes no sense to have K >= g")

+      argsFull <- paste0("-K ", K, ifelse(r, " -r ", " "),

+                         "-i ", i, " -I ", I , ifelse(circular, " -c ", " "),

+                         "-e ", e, " -E ", E, " ",

+                         "-H ", H, " ",

+                         "-s ", s, " -S ", S, " -d ", d, " ",

+                         "-o ", o, " -O ", O, " -p ", p, " -P ", P, " ",

+                         g, " 2")

+      fl1 <- system2(fl_generate_binary(), args = argsFull, stdout = TRUE)[-1]

+    }

+    if (model == "Eggbox" || model == "Ising" || model == "Full" || model == "NK") {

             fl1 <- matrix(

                 as.numeric(unlist(strsplit(paste(fl1, collapse = " "), " "))),

                 ncol = g + 1, byrow = TRUE)

@@ -181,11 +181,25 @@ rfitness <- function(g, c= 0.5,

                             "not include 1")

             }

         }

-        if(truncate_at_0) {

-            fi[fi <= 0] <- 1e-9

-        }

         if(log) {

-            fi <- log(fi/fi[1]) + 1

+            ## If you want logs, you certainly do not want

+            ## the log of a negative number

+            fi[fi < 0] <- 0

+            if(wt_is_1 == "no") {

+                fi <- log(fi)

+            } else {

+                ## by decree, fitness of wt is 1. So shift everything

+                fi <- log(fi) + 1

+            }

+            ## former expression, but it was more confusing

+            ## fi <- log(fi/fi[1]) + 1

+            

+        }

+        if(truncate_at_0) {

+            ## yes, truncate but add noise to prevent identical

+            fi[fi <= 0] <- runif(sum(fi <= 0),

+                                 min = 1e-10,

+                                 max = 1e-9)

         }

         m <- cbind(m, Fitness = fi)

         if(!is_null_na(min_accessible_genotypes)) {

@@ -14,6 +14,48 @@

 ## along with this program.  If not, see <http://www.gnu.org/licenses/>.

+create_eq_ref <- function(g) {

+    ## "random" gives more prob. to genotypes with

+    ## number of mutated genes close to g/2.

+    ## This gives equal prob to having the reference

+    ## be of any of the possible number of mutated genes.

+    nm <- sample(g, 1)

+    ref <- c(rep(1, nm), rep(0, g - nm))

+    sample(ref)

+}

+

+

+

+get_magellan_binaries <-

+    function(names_binaries = c("fl_statistics", "fl_generate", "fl_genchains"))

+{

+    rarch <- Sys.getenv('R_ARCH')

+    nn_names_binaries <- names_binaries

+    if(.Platform$OS.type == "windows")

+        names_binaries <- paste0(names_binaries, ".exe")

+    if(nzchar(rarch)) {

+        rarch <- sub("^/", "", rarch)

+        magellan_binaries <-  system.file(package = "OncoSimulR", "exec",

+                                          rarch, names_binaries)

+    } else {

+        magellan_binaries <-  system.file(package = "OncoSimulR", "exec",

+                                          names_binaries)

+    }

+    names(magellan_binaries) <- nn_names_binaries

+    return(magellan_binaries)

+}

+

+## pkg.env <- new.env()

+## ## The next will not install with error

+## ## ERROR: hard-coded installation path:

+## please report to the package maintainer and use ‘--no-staged-install’

+## pkg.env <- c(pkg.env, get_magellan_binaries())

+

+

+

+fl_statistics_binary <- function() get_magellan_binaries("fl_statistics")

+fl_generate_binary <- function() get_magellan_binaries("fl_generate")

+

 rfitness <- function(g, c= 0.5,

                      sd = 1,

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

                      reference = "random", ## "random", "max", or the vector,

                                      ## e.g., rep(1, g). If random, a

                                      ## random genotype is chosen as

-                                     ## reference. If "max" this is rep(1, g)

+                     ## reference. If "max" this is rep(1, g)

                      scale = NULL, ## a two-element vector: min and max

                      wt_is_1 = c("subtract", "divide", "force", "no"),

                      ## wt_is_1 = TRUE, ## wt has fitness 1

@@ -30,7 +72,10 @@ rfitness <- function(g, c= 0.5,

                                  ## wt_is_1 = divide

                      min_accessible_genotypes = NULL,

                      accessible_th = 0,

-                     truncate_at_0 = TRUE) {

+                     truncate_at_0 = TRUE,

+                     K = 1,

+                     r = TRUE,

+                     model = c("RMF", "NK")) {

     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene

     ## and Crona, 2014. And this allows moving from HoC to purely additive

     ## changing c and sd.

@@ -38,29 +83,77 @@ rfitness <- function(g, c= 0.5,

     ## FIXME future: do this with order too?

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

     wt_is_1 = match.arg(wt_is_1)

+    model = match.arg(model)

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

     m <- generate_matrix_genotypes(g)

     done <- FALSE

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

     while(!done) {

-        ## attempts <- attempts + 1

-        f_r <- rnorm(nrow(m), mean = mu, sd = sd)

-        if(inherits(reference, "character") && length(reference) == 1) {

-            if(reference == "random") {

-                referenceI <- m[sample(nrow(m), 1), ]

-            } else if(reference == "max") {

-                referenceI <- rep(1, g)

-            } else if(reference == "random2") {

-                referenceI <- create_eq_ref(g)

-            }

-        } else {

-            referenceI <- reference

+        if(model == "RMF") {

+            ## attempts <- attempts + 1

+            f_r <- rnorm(nrow(m), mean = mu, sd = sd)

+            if(inherits(reference, "character") && length(reference) == 1) {

+                if(reference == "random") {

+                    referenceI <- m[sample(nrow(m), 1), ]

+                } else if(reference == "max") {

+                    referenceI <- rep(1, g)

+                } else if(reference == "random2") {

+                    referenceI <- create_eq_ref(g)

+                }

+            } else {

+                referenceI <- reference

             }

-        d_reference <- apply(m, 1, function(x) sum(abs(x - referenceI)))

-        f_det <- -c * d_reference

-        ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona

-        fi <- f_r + f_det

-        

+            d_reference <- apply(m, 1, function(x) sum(abs(x - referenceI)))

+            f_det <- -c * d_reference

+            ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona

+            fi <- f_r + f_det

+        } else if(model == "NK") {

+            if(K >= g) stop("It makes no sense to have K >= g")

+            argsnk <- paste0("-K ", K,

+                             ifelse(r, " -r ", " "),

+                             g, " 2")

+            fl1 <- system2(fl_generate_binary(), args = argsnk, stdout = TRUE)[-1]

+            fl1 <- matrix(

+                as.numeric(unlist(strsplit(paste(fl1, collapse = " "), " "))),

+                ncol = g + 1, byrow = TRUE)

+            m1 <- fl1[, 1:g]

+            fi <- fl1[, g + 1]

+

+            ## For scaling, etc, all that matters, if anything, is the wildtype

+

+            ## We could order by doing this

+            ## But I am not 100% sure this will always be the same as

+            ## generate_matrix_genotypes

+            ## oo <- do.call(order,

+            ##               c(list(muts),

+            ##                 as.list(data.frame(m1[, rev(1:ncol(m1))]))))

+            ## m2 <- m1[oo, ]

+            ## Or create an id and order by it?

+

+            ## When we get to 20 genes, this is slow, about 18 seconds each

+            ## apply. Matching is fast (< 0.5 seconds)

+            gtstring <- apply(m, 1, function(x) paste0(x, collapse = ""))

+            gtstring2 <- apply(m1, 1, function(x) paste0(x, collapse = ""))

+            oo <- match(gtstring, gtstring2)

+            fi <- fi[oo]

+            ## make sure no left overs

+            rm(gtstring, gtstring2, oo, fl1, m1)

+

+            ## Had we not ordered, do this!!!

+            ## Which one is WT?

+            ## muts <- rowSums(m1)

+            ## w_wt <- which(muts == 0)

+            ## if(w_wt != 1) {

+            ##     f_a <- fi[1]

+            ##     fi[1] <- fi[w_wt]

+            ##     fi[w_wt] <- f_a

+            ##     rm(f_a)

+            ## }

+            ## m[] <- m1

+            ## rm(m1)

+            ## rm(fl1)

+        }

+

         if(!is.null(scale)) {

             fi <- (fi - min(fi))/(max(fi) - min(fi))

             fi <- scale[1] + fi * (scale[2] - scale[1])

@@ -119,16 +212,114 @@ rfitness <- function(g, c= 0.5,

 }

-create_eq_ref <- function(g) {

-    ## "random" gives more prob. to genotypes with

-    ## number of mutated genes close to g/2.

-    ## This gives equal prob to having the reference

-    ## be of any of the possible number of mutated genes.

-    nm <- sample(g, 1)

-    ref <- c(rep(1, nm), rep(0, g - nm))

-    sample(ref)

-}

+

+

+

+

+## rfitness <- function(g, c= 0.5,

+##                      sd = 1,

+##                      mu = 1,

+##                      reference = "random", ## "random", "max", or the vector,

+##                                      ## e.g., rep(1, g). If random, a

+##                                      ## random genotype is chosen as

+##                                      ## reference. If "max" this is rep(1, g)

+##                      scale = NULL, ## a two-element vector: min and max

+##                      wt_is_1 = c("subtract", "divide", "force", "no"),

+##                      ## wt_is_1 = TRUE, ## wt has fitness 1

+##                      log = FALSE, ## log only makes sense if all values >

+##                                  ## 0. scale with min > 0, and/or set

+##                                  ## wt_is_1 = divide

+##                      min_accessible_genotypes = NULL,

+##                      accessible_th = 0,

+##                      truncate_at_0 = TRUE) {

+##     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene

+##     ## and Crona, 2014. And this allows moving from HoC to purely additive

+##     ## changing c and sd.

+

+##     ## FIXME future: do this with order too?

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

+##     wt_is_1 = match.arg(wt_is_1)

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

+##     m <- generate_matrix_genotypes(g)

+##     done <- FALSE

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

+##     while(!done) {

+##         ## attempts <- attempts + 1

+##         f_r <- rnorm(nrow(m), mean = mu, sd = sd)

+##         if(inherits(reference, "character") && length(reference) == 1) {

+##             if(reference == "random") {

+##                 referenceI <- m[sample(nrow(m), 1), ]

+##             } else if(reference == "max") {

+##                 referenceI <- rep(1, g)

+##             } else if(reference == "random2") {

+##                 referenceI <- create_eq_ref(g)

+##             }

+##         } else {

+##             referenceI <- reference

+##             }

+##         d_reference <- apply(m, 1, function(x) sum(abs(x - referenceI)))

+##         f_det <- -c * d_reference

+##         ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona

+##         fi <- f_r + f_det

+

+##         if(!is.null(scale)) {

+##             fi <- (fi - min(fi))/(max(fi) - min(fi))

+##             fi <- scale[1] + fi * (scale[2] - scale[1])

+##         }

+##         if(wt_is_1 == "divide") {

+##             ## We need to shift to avoid ratios involving negative numbers and

+##             ## we need to avoid having any fitness at 0, specially the wt.  If

+##             ## any negative value, add the min, and shift by the magnitude of

+##             ## the min to avoid any at 0.

+

+##             ## If you use scale and wt_is_1, this will move the scale. It is

+##             ## not possible to obtain a linear transformation that will keep

+##             ## the min and max of the scale, and wt at 1.

+##             min_fi <- min(fi)

+##             if(min_fi < 0)

+##                 fi <- fi + 2 * abs(min(fi))

+##             fi <- fi/fi[1]

+##         } else if (wt_is_1 == "subtract") {

+##             fi <- fi - fi[1] + 1.0

+##         } else if(wt_is_1 == "force") {

+##             fi[1] <- 1.0

+##             if(!is.null(scale)) {

+##                 if( (1 > scale[2]) || (1 < scale[1]))

+##                     warning("Using wt_is_1 = force and scale, but scale does ",

+##                             "not include 1")

+##             }

+##         }

+##         if(truncate_at_0) {

+##             fi[fi <= 0] <- 1e-9

+##         }

+##         if(log) {

+##             fi <- log(fi/fi[1]) + 1

+##         }

+##         m <- cbind(m, Fitness = fi)

+##         if(!is_null_na(min_accessible_genotypes)) {

+##             ## num_accessible_genotypes <- count_accessible_g(m, accessible_th)

+##             ## Recall accessibleGenotypes includes the wt, if accessible.

+##             num_accessible_genotypes <- length(wrap_accessibleGenotypes(m, accessible_th)) - 1

+##             ## message("\n     num_accessible_genotypes = ", num_accessible_genotypes, "\n")

+##             if(num_accessible_genotypes >= min_accessible_genotypes) {

+##                 done <- TRUE

+##                 attributes(m) <- c(attributes(m),

+##                                    accessible_genotypes = num_accessible_genotypes,

+##                                    accessible_th = accessible_th)

+##             } else {

+##                 ## Cannot start again with a fitness column

+##                 m <- m[, -ncol(m), drop = FALSE]

+##             }

+##         } else {

+##             done <- TRUE

+##         }

+##     }

+##     ## message("\n number of attempts = ", attempts, "\n")

+##     class(m) <- c(class(m), "genotype_fitness_matrix")

+##     return(m)

+## }

+

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

     ## FIXME future: do this with order too?

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

     wt_is_1 = match.arg(wt_is_1)

-    

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

     m <- generate_matrix_genotypes(g)

     done <- FALSE

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

@@ -1,26 +1,50 @@

+## Copyright 2013, 2014, 2015, 2016, 2017 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/>.

+

+

+

 rfitness <- function(g, c= 0.5,

                      sd = 1,

+                     mu = 1,

                      reference = "random", ## "random", "max", or the vector,

                                      ## e.g., rep(1, g). If random, a

                                      ## random genotype is chosen as

                                      ## reference. If "max" this is rep(1, g)

                      scale = NULL, ## a two-element vector: min and max

-                     wt_is_1 = TRUE, ## wt has fitness 1

+                     wt_is_1 = c("subtract", "divide", "force", "no"),

+                     ## wt_is_1 = TRUE, ## wt has fitness 1

                      log = FALSE, ## log only makes sense if all values >

                                  ## 0. scale with min > 0, and/or set

-                                 ## wt_is_1 = TRUE

-                     min_accessible_genotypes = 0,

-                     accessible_th = 0) {

+                                 ## wt_is_1 = divide

+                     min_accessible_genotypes = NULL,

+                     accessible_th = 0,

+                     truncate_at_0 = TRUE) {

     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene

     ## and Crona, 2014. And this allows moving from HoC to purely additive

     ## changing c and sd.

     ## FIXME future: do this with order too?

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

+    wt_is_1 = match.arg(wt_is_1)

+    

     m <- generate_matrix_genotypes(g)

     done <- FALSE

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

     while(!done) {

-        f_r <- rnorm(nrow(m), mean = 0, sd = sd)

+        ## attempts <- attempts + 1

+        f_r <- rnorm(nrow(m), mean = mu, sd = sd)

         if(inherits(reference, "character") && length(reference) == 1) {

             if(reference == "random") {

                 referenceI <- m[sample(nrow(m), 1), ]

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

             fi <- (fi - min(fi))/(max(fi) - min(fi))

             fi <- scale[1] + fi * (scale[2] - scale[1])

         }

-        if(wt_is_1) {

+        if(wt_is_1 == "divide") {

             ## We need to shift to avoid ratios involving negative numbers and

             ## we need to avoid having any fitness at 0, specially the wt.  If

             ## any negative value, add the min, and shift by the magnitude of

@@ -54,16 +78,28 @@ rfitness <- function(g, c= 0.5,

             if(min_fi < 0)

                 fi <- fi + 2 * abs(min(fi))

             fi <- fi/fi[1]

+        } else if (wt_is_1 == "subtract") {

+            fi <- fi - fi[1] + 1.0

+        } else if(wt_is_1 == "force") {

+            fi[1] <- 1.0

+            if(!is.null(scale)) {

+                if( (1 > scale[2]) || (1 < scale[1]))

+                    warning("Using wt_is_1 = force and scale, but scale does ",

+                            "not include 1")

+            }

+        }

+        if(truncate_at_0) {

+            fi[fi <= 0] <- 1e-9

         }

-        

         if(log) {

             fi <- log(fi/fi[1]) + 1

         }

         m <- cbind(m, Fitness = fi)

-        if(min_accessible_genotypes > 0) {

+        if(!is_null_na(min_accessible_genotypes)) {

             ## num_accessible_genotypes <- count_accessible_g(m, accessible_th)

             ## Recall accessibleGenotypes includes the wt, if accessible.

             num_accessible_genotypes <- length(wrap_accessibleGenotypes(m, accessible_th)) - 1

+            ## message("\n     num_accessible_genotypes = ", num_accessible_genotypes, "\n")

             if(num_accessible_genotypes >= min_accessible_genotypes) {

                 done <- TRUE

                 attributes(m) <- c(attributes(m),

@@ -77,6 +113,7 @@ rfitness <- function(g, c= 0.5,

             done <- TRUE

         }

     }

+    ## message("\n number of attempts = ", attempts, "\n")

     class(m) <- c(class(m), "genotype_fitness_matrix")

     return(m)

 }

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

         }

         m <- cbind(m, Fitness = fi)

         if(min_accessible_genotypes > 0) {

-            num_accessible_genotypes <- count_accessible_g(m, accessible_th)

+            ## num_accessible_genotypes <- count_accessible_g(m, accessible_th)

+            ## Recall accessibleGenotypes includes the wt, if accessible.

+            num_accessible_genotypes <- length(wrap_accessibleGenotypes(m, accessible_th)) - 1

             if(num_accessible_genotypes >= min_accessible_genotypes) {

                 done <- TRUE

                 attributes(m) <- c(attributes(m),

@@ -89,3 +91,7 @@ create_eq_ref <- function(g) {

     ref <- c(rep(1, nm), rep(0, g - nm))

     sample(ref)

 }

+

+

+

+

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

                 referenceI <- m[sample(nrow(m), 1), ]

             } else if(reference == "max") {

                 referenceI <- rep(1, g)

-            } 

+            } else if(reference == "random2") {

+                referenceI <- create_eq_ref(g)

+            }

         } else {

             referenceI <- reference

             }

@@ -76,3 +78,14 @@ rfitness <- function(g, c= 0.5,

     class(m) <- c(class(m), "genotype_fitness_matrix")

     return(m)

 }

+

+

+create_eq_ref <- function(g) {

+    ## "random" gives more prob. to genotypes with

+    ## number of mutated genes close to g/2.

+    ## This gives equal prob to having the reference

+    ## be of any of the possible number of mutated genes.

+    nm <- sample(g, 1)

+    ref <- c(rep(1, nm), rep(0, g - nm))

+    sample(ref)

+}

@@ -59,8 +59,12 @@ rfitness <- function(g, c= 0.5,

         }

         m <- cbind(m, Fitness = fi)

         if(min_accessible_genotypes > 0) {

-            if(count_accessible_g(m, accessible_th) >= min_accessible_genotypes) {

+            num_accessible_genotypes <- count_accessible_g(m, accessible_th)

+            if(num_accessible_genotypes >= min_accessible_genotypes) {

                 done <- TRUE

+                attributes(m) <- c(attributes(m),

+                                   accessible_genotypes = num_accessible_genotypes,

+                                   accessible_th = accessible_th)

             } else {

                 ## Cannot start again with a fitness column

                 m <- m[, -ncol(m), drop = FALSE]

@@ -6,10 +6,11 @@ rfitness <- function(g, c= 0.5,

                                      ## reference. If "max" this is rep(1, g)

                      scale = NULL, ## a two-element vector: min and max

                      wt_is_1 = TRUE, ## wt has fitness 1

-                     log = FALSE ## log only makes sense if all values >

+                     log = FALSE, ## log only makes sense if all values >

                                  ## 0. scale with min > 0, and/or set

                                  ## wt_is_1 = TRUE

-                     ) {

+                     min_accessible_genotypes = 0,

+                     accessible_th = 0) {

     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene

     ## and Crona, 2014. And this allows moving from HoC to purely additive

     ## changing c and sd.

@@ -17,42 +18,57 @@ rfitness <- function(g, c= 0.5,

     ## FIXME future: do this with order too?

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

     m <- generate_matrix_genotypes(g)

-    f_r <- rnorm(nrow(m), mean = 0, sd = sd)

-    if(inherits(reference, "character") && length(reference) == 1) {

-        if(reference == "random") {

-            reference <- m[sample(nrow(m), 1), ]

-        } else if(reference == "max") {

-            reference <- rep(1, g)

+    done <- FALSE

+    while(!done) {

+        f_r <- rnorm(nrow(m), mean = 0, sd = sd)

+        if(inherits(reference, "character") && length(reference) == 1) {

+            if(reference == "random") {

+                referenceI <- m[sample(nrow(m), 1), ]

+            } else if(reference == "max") {

+                referenceI <- rep(1, g)

+            } 

+        } else {

+            referenceI <- reference