Browse code

2.99.4

ramon diaz-uriarte (at Phelsuma) authored on 17/12/2020 15:07:07
Showing24 changed files

... ...
@@ -1,8 +1,8 @@
1 1
 Package: OncoSimulR
2 2
 Type: Package
3 3
 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 
4
-Version: 2.99.3
5
-Date: 2020-12-13
4
+Version: 2.99.4
5
+Date: 2020-12-17
6 6
 Authors@R: c(
7 7
 	      person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),	
8 8
  	   		     email = "rdiaz02@gmail.com"),
... ...
@@ -53,17 +53,7 @@ phcl_from_lod <- function(df, x) {
53 53
     ## df <- df[!duplicated(df[, c(1, 2)]), , drop = FALSE]
54 54
     
55 55
     tG <- unique(c(as.character(df[, 1]), as.character(df[, 2])))
56
-    ## ## Do as in phylogClone. So that we have the same nodes
57
-    ## ## in LOD all and not LOD all?
58
-    ## z <- which_N_at_T(x, N = 1, t = "last")
59
-    ## tG <- x$GenotypesLabels[z]
60
-    
61
-    ## ## FIXME: aren't these two warnings redundant or aliased?
62
-    ## ## yes, partlt
63
-    ##  I think this can never happen now
64
-    ## if ((length(tG) == 1) && (tG == "")) {
65
-    ##     warning("There never was a descendant of WT")
66
-    ## }
56
+
67 57
     if (nrow(df) == 0) {
68 58
         warning("LOD structure has 0 rows: no descendants of initMutant ever appeared. ")
69 59
         return(NA)
... ...
@@ -149,15 +139,6 @@ diversityLOD <- function(llod) {
149 139
     shannonI(table(pathstr))
150 140
 }
151 141
 
152
-## diversityLOD <- function(llod) {
153
-##     nn <- names(llod[[1]])
154
-##     if( is_null_na(nn) ||
155
-##         !(is.list(llod)))
156
-##         stop("Object must be a list of LODs")
157
-##     pathstr <- unlist(lapply(llod, function(x) paste(names(x),
158
-##                                                      collapse = "_")))
159
-##     shannonI(table(pathstr))
160
-## }
161 142
 
162 143
 LOD_as_path <- function(llod) {
163 144
     path_l <- function(u) {
... ...
@@ -223,44 +204,6 @@ LOD.oncosimulpop <- function(x) return(lapply(x, LOD.internal))
223 204
 
224 205
 
225 206
 
226
-## POM.oncosimul2 <- function(x) {
227
-##     out <- POM.internal(x)
228
-##     class(out) <- c(class(out), "oncosimul_pom")
229
-##     return(out)
230
-## }
231
-
232
-## LOD.oncosimul2 <- function(x) {
233
-##     out <- LOD.internal(x)
234
-##     class(out) <- c(class(out), "oncosimul_lod")
235
-##     return(out)
236
-## }
237
-
238
-
239
-## POM.oncosimulpop <- function(x) {
240
-##     out <- lapply(x, POM.internal)
241
-##     class(out) <- c(class(out), "oncosimul_pom_list")
242
-##     return(out)
243
-## }
244
-
245
-## LOD.oncosimulpop <- function(x) {
246
-##     out <- lapply(x, LOD.internal)
247
-##     class(out) <- c(class(out), "oncosimul_lod_list")
248
-##     return(out)
249
-## }
250
-
251
-
252
-## summary.oncosimul_lod_list <- function(x) {
253
-##     cat("List of ", length(x), " simulations\n.")
254
-##     cat("Shannon's diversity (entropy) = ", diversityLOD(x), "\n")
255
-## }
256
-
257
-## summary.oncosimul_pom_list <- function(x) {
258
-##     cat("List of ", length(x), " simulations\n.")
259
-##     cat("Shannon's diversity (entropy) = ", diversityPOM(x), "\n")
260
-## }
261
-
262
-
263
-
264 207
 POM_pre_2.9.2 <- function(x) {
265 208
     if(is.null(x$pops.by.time)) {
266 209
         warning("Missing needed components. This might be a failed simulation.",
... ...
@@ -369,45 +312,6 @@ LOD.internal_pre_2.9.2 <- function(x, strict) {
369 312
 
370 313
 
371 314
 
372
-## LOD_as_path_pre_2.9.2 <- function(llod) {
373
-##     path_l <- function(u) {
374
-##         if(length(u$lod_single) == 1) {
375
-##             if(is.null(attributes(u)$initMutant))
376
-##                 initMutant <- ""
377
-##             else
378
-##                 initMutant <- attributes(u)$initMutant
379
-##             if(initMutant == "") initMutant <- "WT"
380
-##             if(grepl("_is_end", u$lod_single))
381
-##                 return(initMutant)
382
-##             if(u$lod_single == "No_descendants")
383
-##                 return(initMutant)
384
-##         } else {
385
-##             ## Deal with "" meaning WT
386
-##             the_names <- names(u$lod_single)
387
-##             the_names_wt <- which(the_names == "")
388
-            
389
-##             if(length(the_names_wt)) {
390
-##                 if(length(the_names_wt) > 1) stop("more than 1 WT?!")
391
-##                 if(the_names_wt > 1) stop("WT in position not 1?!")
392
-##                 the_names[the_names_wt] <- "WT"
393
-##             }
394
-##             return(paste(the_names, collapse = " -> ")) 
395
-##             ## return(paste0("WT", paste(names(u$lod_single),
396
-##             ##                           collapse = " -> ")) )
397
-##         }
398
-##     }
399
-##     if(identical(names(llod), c("all_paths", "lod_single")))
400
-##         pathstr <- path_l(llod)
401
-##     else {
402
-##         ## should be a list
403
-##         pathstr <- unlist(lapply(llod, path_l))
404
-##     }
405
-##     return(pathstr)
406
-##     ## pathstr <- unlist(lapply(llod, function(x) paste(names(x$lod_single),
407
-##     ##                                                  collapse = " -> ")))
408
-##     ## return(paste0("WT", pathstr))
409
-## }
410
-
411 315
 
412 316
 
413 317
 LOD.oncosimul2_pre_2.9.2 <- function(x, strict = TRUE)
... ...
@@ -16,12 +16,9 @@
16 16
 
17 17
 ## Posets, restriction tables, etc.
18 18
 
19
-
20 19
 ## This is all too complicated. Have here just the minimal set we need,
21 20
 ## because we will soon be changing formats.
22 21
 
23
-
24
-
25 22
 ## When we go poset -> rT or
26 23
 ##  adjMat -> rT
27 24
 
... ...
@@ -151,35 +148,7 @@ adjmat.to.restrictTable <- function(x, root = FALSE,
151 148
                              orderedNames = TRUE)
152 149
     if(root)
153 150
         x <- x[-1, -1, drop = FALSE]
154
-    ## ## we have the zero
155
-    ## if( any(colnames(x) %in% c("0", "root", "Root")) & !root)
156
-    ##     warning("Looks like the matrix has a root but you specified root = FALSE")
157
-
158
-    ## if(!identical(colnames(x), rownames(x)))
159
-    ##     stop("colnames and rownames not identical")
160
-    ## if(root) {
161
-    ##     posRoot <- which(colnames(x) %in% rootNames)
162
-    ##     if(!length(posRoot))
163
-    ##         stop("No column with the root name")
164
-    ##     if(length(posRoot) > 1)
165
-    ##         stop("Ambiguous location of root")
166
-    ##     x <- x[-posRoot, -posRoot]
167
-    ## }
168
-
169
-    ## if(typeof(x) != "integer")
170
-    ##     warning("This is not an _integer_ adjacency matrix")
171
-    ## if( !all(x %in% c(0, 1) ))
172
-    ##     stop("Values not in [0, 1]")
173
-
174
-    ## if(!is.null(colnames(x))) {
175
-    ##     ## FIXME: this makes sense with numeric labels for columns, but
176
-    ##     ## not ow.
177
-    ##     oi <- order(as.numeric(colnames(x)))
178
-    ##     if(any(oi != (1:ncol(x)))) {
179
-    ##         warning("Reordering adjacency matrix")
180
-    ##         x <- x[oi, oi]
181
-    ##     }
182
-    ## }
151
+ 
183 152
     
184 153
     num.deps <- colSums(x)
185 154
     max.n.deps <- max(num.deps)
... ...
@@ -205,29 +174,11 @@ adjmat.to.restrictTable <- function(x, root = FALSE,
205 174
 OTtoPoset <- function(x) {
206 175
     checkProperOTAdjMat(x)
207 176
 
208
-    ## ## root must always be dropped, as we are creating a poset
209
-    ## dropRoot <- TRUE
210
-    ## if(!dropRoot)
211
-    ##     warning("Are you sure you do not want dropRoot?")
212
-    ## if(dropRoot) {
213
-    ##     ncx <- ncol(x)
214
-    ##     x <- x[-1, -1]
215
-    ##     ## y <- (which(x == 1L, arr.ind = TRUE) )
216
-    ##     ## if(nrow(y) == 0) ## all nodes descend from 0
217
-    ##     ##     y <- cbind(0L, ncx - 1L)
218
-    ##     namesInts <- type.convert(colnames(x), as.is = TRUE)
219
-        
220
-    ## } else {
221
-    ##     ## y <- (which(x == 1L, arr.ind = TRUE) )
222
-    ##     namesInts <- c(0L, type.convert(colnames(x)[-1], as.is = TRUE))
223
-    ## }
177
+
224 178
 
225 179
     ncx <- ncol(x)
226 180
     x <- x[-1, -1, drop = FALSE]
227
-    ## y <- (which(x == 1L, arr.ind = TRUE) )
228
-    ## if(nrow(y) == 0) ## all nodes descend from 0
229
-    ##     y <- cbind(0L, ncx - 1L)
230
-    namesInts <- type.convert(colnames(x), as.is = TRUE)
181
+  namesInts <- type.convert(colnames(x), as.is = TRUE)
231 182
     
232 183
     if(!is.integer(namesInts))
233 184
         stop("cannot convert to poset adj mat with non-int colnames")
... ...
@@ -242,34 +193,6 @@ OTtoPoset <- function(x) {
242 193
     return(p2)
243 194
 }    
244 195
 
245
-## the next is just a convenience
246
-## sortAdjMat <- function(am) {
247
-##     cn <- colnames(am)
248
-##     rootpos <- grep("^Root$", cn) 
249
-##     if(length(rootpos) != 1)
250
-##         stop("No root in adj mat, or multiple Roots")
251
-##     cn <- c("Root", sort(colnames(am)[-rootpos]))
252
-##     return(am[cn, cn])
253
-## }
254
-
255
-
256
-## No longer used
257
-## sortAdjMat <- function(am) {
258
-##     ## If column names, except Root, are integers, sort as integers. O.w.,
259
-##     ## general lexicog. sort.
260
-##     cn <- colnames(am)
261
-##     rootpos <- grep("^Root$", cn) 
262
-##     if(length(rootpos) != 1)
263
-##         stop("No root in adj mat, or multiple Roots")
264
-##     cn0 <- colnames(am)[-rootpos]
265
-##     namesInts <- type.convert(cn0, as.is = TRUE)
266
-##     if(is.integer(namesInts)) {
267
-##         cn <- c("Root", sort(namesInts))
268
-##     } else {
269
-##         cn <- c("Root", sort(cn0))
270
-##     }
271
-##     return(am[cn, cn])
272
-## }
273 196
 
274 197
 
275 198
 
... ...
@@ -16,12 +16,6 @@
16 16
 ## plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-
17 17
 ##     plot.genotype_fitness_matrix <- plotFitnessLandscape 
18 18
 
19
-## FIXME: show only accessible paths? 
20
-## FIXME: when show_labels = FALSE we still show the boxes
21
-##        and some of the labels.!!!
22
-## FIXME: if using only_accessible, maybe we
23
-## can try to use fast_peaks, and use the slower
24
-## approach as fallback (if identical fitness)
25 19
 plotFitnessLandscape <- function(x, show_labels = TRUE,
26 20
                                  col = c("green4", "red", "yellow"),
27 21
                                  lty = c(1, 2, 3), use_ggrepel = FALSE,
... ...
@@ -31,7 +25,12 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,
31 25
                                  ...) {
32 26
 
33 27
     ## FIXME future:
34
-
28
+    ## FIXME: when show_labels = FALSE we still show the boxes
29
+    ##        and some of the labels.!!!
30
+    ## FIXME: if using only_accessible, maybe we
31
+    ## can try to use fast_peaks, and use the slower
32
+    ## approach as fallback (if identical fitness)
33
+    
35 34
     ## - Allow passing order effects. Change "allGenotypes_to_matrix"
36 35
     ##   below? Probably not, as we cannot put order effects as a
37 36
     ##   matrix. Do something else, like allow only order effects if from
... ...
@@ -70,14 +69,7 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,
70 69
     }
71 70
     vv <- which(!is.na(gaj), arr.ind = TRUE)
72 71
     
73
-    ## plot(x = mutated, y = e1$Fitness, ylab = "Fitness",
74
-    ##      xlab = "", type = "n", axes = FALSE)
75
-    ## box()
76
-    ## axis(2)
77
-    ## text(x = mutated, y = x$Fitness, labels = x$Genotype)
78
-
79
-    ## The R CMD CHEKC notes about no visible binding for global variable
80
-
72
+  
81 73
     x_from <- y_from <- x_to <- y_to <- Change <- muts <-
82 74
         label <- fitness <- Type <- NULL
83 75
                 
... ...
@@ -184,27 +176,6 @@ plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-
184 176
 ######################################################################
185 177
 
186 178
 
187
-## wrap_filter_inaccessible <- function(x, max_num_genotypes, accessible_th) {
188
-##     ## wrap it, for my consumption
189
-##     tfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)
190
-##     mutated <- rowSums(tfm$gfm[, -ncol(tfm$gfm)])
191
-##     gaj <- genot_to_adj_mat(tfm$gfm)
192
-##     gaj <- filter_inaccessible(gaj, accessible_th)
193
-##     remaining <- as.numeric(colnames(gaj))
194
-##     mutated <- mutated[remaining]
195
-##     tfm$afe <- tfm$afe[remaining, , drop = FALSE]
196
-##     return(list(remaining = remaining,
197
-##                 mutated = mutated,
198
-##                 tfm = tfm))
199
-## }
200
-
201
-## No longer being used. Used to be in rfitness
202
-## count_accessible_g <- function(gfm, accessible_th) {
203
-##     gaj <- genot_to_adj_mat(gfm)
204
-##     gaj <- filter_inaccessible(gaj, accessible_th)
205
-##     return(ncol(gaj) - 1)
206
-## }
207
-
208 179
 
209 180
 ## There is now C++ code to get just the locations/positions of the
210 181
 ## accessible genotypes
... ...
@@ -47,18 +47,6 @@ to_Magellan <- function(x, file,
47 47
 
48 48
 
49 49
 
50
-## ## genotype_fitness_matrix -> fitness landscape as data frame
51
-## fmatrix_to_afe <- function(x) {
52
-##     stopifnot(inherits(x, "genotype_fitness_matrix"))
53
-##     y <- x[, -ncol(x)]
54
-##     nn <- apply(y, 1,
55
-##                 function(u) paste(sort(colnames(y)[as.logical(u)]),
56
-##                                   collapse = ", "))
57
-##     nn[nn == ""] <- "WT"
58
-##     return(data.frame(Genotype = nn, Fitness = x[, ncol(x)],
59
-##            stringsAsFactors = FALSE))
60
-## }
61
-
62 50
 to_Fitness_Matrix <- function(x, max_num_genotypes) {
63 51
     ## A general converter. Ready to be used by plotFitnessLandscape and
64 52
     ## Magellan exporter.
... ...
@@ -124,217 +112,181 @@ to_Fitness_Matrix <- function(x, max_num_genotypes) {
124 112
 ## but if we are passed a fitness landscapes as produced by
125 113
 ## rfitness, do nothing. Well, it actually does something.
126 114
 
127
-
128
-##Modified
129 115
 to_genotFitness_std <- function(x,
130 116
                                 frequencyDependentFitness = FALSE,
131 117
                                 frequencyType = frequencyType,
132 118
                                 simplify = TRUE,
133 119
                                 min_filter_fitness = 1e-9,
134 120
                                 sort_gene_names = TRUE) {
135
-  ## Would break with output from allFitnessEffects and
136
-  ## output from allGenotypeAndMut
137
-
138
-  ## For the very special and weird case of
139
-  ## a matrix but only a single gene so with a 0 and 1
140
-  ## No, this is a silly and meaningless case.
141
-  ## if( ( ncol(x) == 2 ) && (nrow(x) == 1) && (x[1, 1] == 1) ) {
142
-
143
-  ## } else  blabla:
144
-
145
-
146
-
147
-  ## if((ncol(x) > 2) && !(inherits(x, "matrix"))
148
-  ##     stop(paste0("Genotype fitness input either two-column data frame",
149
-  ##          " or a numeric matrix with > 2 columns."))
150
-  ## if( (ncol(x) > 2) && (nrow(x) == 1) )
151
-  ##     stop(paste0("It looks like you have a matrix for a single genotype",
152
-  ##                 " of a single gene. For this degenerate cases use",
153
-  ##                 " a data frame specification."))
154
-
155
-
121
+    ## Would break with output from allFitnessEffects and
122
+    ## output from allGenotypeAndMut
156 123
 
157 124
     if(! (inherits(x, "matrix") || inherits(x, "data.frame")) )
158
-      stop("Input must inherit from matrix or data.frame.")
125
+        stop("Input must inherit from matrix or data.frame.")
159 126
 
160 127
     if(ncol(x) > 2) {
161
-      if (!frequencyDependentFitness){
162
-        if(inherits(x, "matrix")) {
163
-          if(!is.numeric(x))
164
-            stop("A genotype fitness matrix/data.frame must be numeric.")
165
-        } else if(inherits(x, "data.frame")) {
166
-          if(!all(unlist(lapply(x, is.numeric))))
167
-            stop("A genotype fitness matrix/data.frame must be numeric.")
168
-        }
169
-      }else{
170
-        if(!inherits(x, "data.frame"))
171
-           stop("Input must inherit from data.frame.")
172
-        if(ncol(x) == 0){
173
-          stop("You have an empty data.frame")
174
-        }
175
-        if(!all(unlist(lapply(x[,-ncol(x)], is.numeric)))){
176
-          stop("All columns except last one must be numeric.")
128
+        if (!frequencyDependentFitness){
129
+            if(inherits(x, "matrix")) {
130
+                if(!is.numeric(x))
131
+                    stop("A genotype fitness matrix/data.frame must be numeric.")
132
+            } else if(inherits(x, "data.frame")) {
133
+                if(!all(unlist(lapply(x, is.numeric))))
134
+                    stop("A genotype fitness matrix/data.frame must be numeric.")
135
+            }
136
+        }else{
137
+            if(!inherits(x, "data.frame"))
138
+                stop("Input must inherit from data.frame.")
139
+            if(ncol(x) == 0){
140
+                stop("You have an empty data.frame")
141
+            }
142
+            if(!all(unlist(lapply(x[,-ncol(x)], is.numeric)))){
143
+                stop("All columns except last one must be numeric.")
144
+            }
145
+            if(is.factor(x[, ncol(x)])) {
146
+                warning("Last column of genotype fitness is a factor. ",
147
+                        "Converting to character.")
148
+                x[, ncol(x)] <- as.character(x[, ncol(x)])
149
+            }
150
+            if(!all(unlist(lapply(x[, ncol(x)], is.character)))){
151
+                stop("All elements in last column must be character.")
152
+            }
177 153
         }
178
-        if(is.factor(x[, ncol(x)])) {
179
-          warning("Last column of genotype fitness is a factor. ",
180
-                  "Converting to character.")
181
-          x[, ncol(x)] <- as.character(x[, ncol(x)])
154
+
155
+        ## We are expecting here a matrix of 0/1 where columns are genes
156
+        ## except for the last column, that is Fitness
157
+        ## Of course, can ONLY work with epistastis, NOT order
158
+        ## return(genot_fitness_to_epistasis(x))
159
+        if(any(duplicated(colnames(x))))
160
+            stop("duplicated column names")
161
+
162
+        cnfl <- which(colnames(x)[-ncol(x)] == "")
163
+        if(length(cnfl)) {
164
+            freeletter <- setdiff(LETTERS, colnames(x))[1]
165
+            if(length(freeletter) == 0) stop("Renaiming failed")
166
+            warning("One column named ''. Renaming to ", freeletter)
167
+            colnames(x)[cnfl] <- freeletter
182 168
         }
183
-        if(!all(unlist(lapply(x[, ncol(x)], is.character)))){
184
-          stop("All elements in last column must be character.")
169
+        if(!is.null(colnames(x)) && sort_gene_names) {
170
+            ncx <- ncol(x)
171
+            cnx <- colnames(x)[-ncx]
172
+            ocnx <- gtools::mixedorder(cnx)
173
+            if(!(identical(cnx[ocnx], cnx))) {
174
+                message("Sorting gene column names alphabetically")
175
+                x <- cbind(x[, ocnx, drop = FALSE], Fitness = x[, (ncx)])
176
+            }
185 177
         }
186
-      }
187
-
188
-      ## We are expecting here a matrix of 0/1 where columns are genes
189
-      ## except for the last column, that is Fitness
190
-      ## Of course, can ONLY work with epistastis, NOT order
191
-      ## return(genot_fitness_to_epistasis(x))
192
-      if(any(duplicated(colnames(x))))
193
-        stop("duplicated column names")
194
-
195
-      cnfl <- which(colnames(x)[-ncol(x)] == "")
196
-      if(length(cnfl)) {
197
-        freeletter <- setdiff(LETTERS, colnames(x))[1]
198
-        if(length(freeletter) == 0) stop("Renaiming failed")
199
-        warning("One column named ''. Renaming to ", freeletter)
200
-        colnames(x)[cnfl] <- freeletter
201
-      }
202
-      if(!is.null(colnames(x)) && sort_gene_names) {
203
-        ncx <- ncol(x)
204
-        cnx <- colnames(x)[-ncx]
205
-        ocnx <- gtools::mixedorder(cnx)
206
-        if(!(identical(cnx[ocnx], cnx))) {
207
-          message("Sorting gene column names alphabetically")
208
-          x <- cbind(x[, ocnx, drop = FALSE], Fitness = x[, (ncx)])
178
+
179
+        if(is.null(colnames(x))) {
180
+            ncx <- (ncol(x) - 1)
181
+            message("No column names: assigning gene names from LETTERS")
182
+            if(ncx > length(LETTERS))
183
+                stop("More genes than LETTERS; please give gene names",
184
+                     " as you see fit.")
185
+            colnames(x) <- c(LETTERS[1:ncx], "Fitness")
209 186
         }
210
-      }
211
-
212
-      if(is.null(colnames(x))) {
213
-        ncx <- (ncol(x) - 1)
214
-        message("No column names: assigning gene names from LETTERS")
215
-        if(ncx > length(LETTERS))
216
-          stop("More genes than LETTERS; please give gene names",
217
-               " as you see fit.")
218
-        colnames(x) <- c(LETTERS[1:ncx], "Fitness")
219
-      }
220
-      if(!all(as.matrix(x[, -ncol(x)]) %in% c(0, 1) ))
221
-        stop("First ncol - 1 entries not in {0, 1}.")
187
+        if(!all(as.matrix(x[, -ncol(x)]) %in% c(0, 1) ))
188
+            stop("First ncol - 1 entries not in {0, 1}.")
222 189
 
223 190
     } else {
224 191
 
225
-      if(!inherits(x, "data.frame"))
226
-        stop("genotFitness: if two-column must be data frame")
227
-      if(ncol(x) == 0){
228
-        stop("You have an empty data.frame")
229
-      }
230
-      ## Make sure no factors
231
-      if(is.factor(x[, 1])) {
232
-        warning("First column of genotype fitness is a factor. ",
233
-                "Converting to character.")
234
-        x[, 1] <- as.character(x[, 1])
235
-      }
236
-      ## Make sure no numbers
237
-      if(any(is.numeric(x[, 1])))
238
-        stop(paste0("genotFitness: first column of data frame is numeric.",
239
-                    " Ambiguous and suggests possible error. If sure,",
240
-                    " enter that column as character"))
241
-
242
-      omarker <- any(grepl(">", x[, 1], fixed = TRUE))
243
-      emarker <- any(grepl(",", x[, 1], fixed = TRUE))
244
-      nogoodepi <- any(grepl(":", x[, 1], fixed = TRUE))
245
-      ## if(omarker && emarker) stop("Specify only epistasis or order, not both.")
246
-      if(nogoodepi && emarker) stop("Specify the genotypes separated by a ',', not ':'.")
247
-      if(nogoodepi && !emarker) stop("Specify the genotypes separated by a ',', not ':'.")
248
-      ## if(nogoodepi && omarker) stop("If you want order, use '>' and if epistasis ','.")
249
-      ## if(!omarker && !emarker) stop("You specified neither epistasis nor order")
250
-      if(omarker) {
251
-        ## do something. To be completed
252
-        stop("This code not yet ready")
253
-        ## You can pass to allFitnessEffects genotype -> fitness mappings that
254
-        ## involve epistasis and order. But they must have different
255
-        ## genes. Otherwise, it is not manageable.
256
-      }
257
-      if( emarker || ( (!omarker) && (!emarker) && (!nogoodepi)) ) {
258
-        ## the second case above corresponds to passing just single letter genotypes
259
-        ## as there is not a single marker
260
-        x <- x[, c(1, 2), drop = FALSE]
261
-        if(!all(colnames(x) == c("Genotype", "Fitness"))) {
262
-          message("Column names of object not Genotype and Fitness.",
263
-                  " Renaming them assuming that is what you wanted")
264
-          colnames(x) <- c("Genotype", "Fitness")
265
-        }
266
-        if((!omarker) && (!emarker) && (!nogoodepi)) {
267
-          message("All single-gene genotypes as input to to_genotFitness_std")
192
+        if(!inherits(x, "data.frame"))
193
+            stop("genotFitness: if two-column must be data frame")
194
+        if(ncol(x) == 0){
195
+            stop("You have an empty data.frame")
268 196
         }
269
-        ## Yes, we need to do this to  scale the fitness and put the "-"
270
-        if(frequencyDependentFitness){
271
-          anywt <- which(x[, 1] == "WT")
272
-          if (length(anywt) > 1){
273
-              stop("WT should not appear more than once in fitness specification")
274
-              ## stop("WT must appear once.")
275
-          }
276
-          ## if(length(anywt) == 0) {
277
-          ##     x <- rbind(data.frame(Genotype = "WT", Fitness = "0"),
278
-          ##                x)
279
-          ## }
280
-          if(is.factor(x[, ncol(x)])) {
281
-            warning("Second column of genotype fitness is a factor. ",
197
+        ## Make sure no factors
198
+        if(is.factor(x[, 1])) {
199
+            warning("First column of genotype fitness is a factor. ",
282 200
                     "Converting to character.")
283
-            x[, ncol(x)] <- as.character(x[, ncol(x)])
284
-          }
201
+            x[, 1] <- as.character(x[, 1])
202
+        }
203
+        ## Make sure no numbers
204
+        if(any(is.numeric(x[, 1])))
205
+            stop(paste0("genotFitness: first column of data frame is numeric.",
206
+                        " Ambiguous and suggests possible error. If sure,",
207
+                        " enter that column as character"))
208
+
209
+        omarker <- any(grepl(">", x[, 1], fixed = TRUE))
210
+        emarker <- any(grepl(",", x[, 1], fixed = TRUE))
211
+        nogoodepi <- any(grepl(":", x[, 1], fixed = TRUE))
212
+        if(nogoodepi && emarker) stop("Specify the genotypes separated by a ',', not ':'.")
213
+        if(nogoodepi && !emarker) stop("Specify the genotypes separated by a ',', not ':'.")
214
+        if(omarker) {
215
+            ## do something. To be completed
216
+            stop("This code not yet ready")
217
+            ## You can pass to allFitnessEffects genotype -> fitness mappings that
218
+            ## involve epistasis and order. But they must have different
219
+            ## genes. Otherwise, it is not manageable.
220
+        }
221
+        if( emarker || ( (!omarker) && (!emarker) && (!nogoodepi)) ) {
222
+            ## the second case above corresponds to passing just single letter genotypes
223
+            ## as there is not a single marker
224
+            x <- x[, c(1, 2), drop = FALSE]
225
+            if(!all(colnames(x) == c("Genotype", "Fitness"))) {
226
+                message("Column names of object not Genotype and Fitness.",
227
+                        " Renaming them assuming that is what you wanted")
228
+                colnames(x) <- c("Genotype", "Fitness")
229
+            }
230
+            if((!omarker) && (!emarker) && (!nogoodepi)) {
231
+                message("All single-gene genotypes as input to to_genotFitness_std")
232
+            }
233
+            ## Yes, we need to do this to  scale the fitness and put the "-"
234
+            if(frequencyDependentFitness){
235
+                anywt <- which(x[, 1] == "WT")
236
+                if (length(anywt) > 1){
237
+                    stop("WT should not appear more than once in fitness specification")
238
+                }
239
+                if(is.factor(x[, ncol(x)])) {
240
+                    warning("Second column of genotype fitness is a factor. ",
241
+                            "Converting to character.")
242
+                    x[, ncol(x)] <- as.character(x[, ncol(x)])
243
+                }
244
+            }
245
+
246
+            x <- allGenotypes_to_matrix(x, frequencyDependentFitness)
285 247
         }
286
-
287
-        x <- allGenotypes_to_matrix(x, frequencyDependentFitness)
288
-      }
289 248
     }
290 249
     ## And, yes, scale all fitnesses by that of the WT
291 250
 
292 251
     if (!frequencyDependentFitness){
293
-      whichroot <- which(rowSums(x[, -ncol(x), drop = FALSE]) == 0)
294
-      if(length(whichroot) == 0) {
295
-        warning("No wildtype in the fitness landscape!!! Adding it with fitness 1.")
296
-        x <- rbind(c(rep(0, ncol(x) - 1), 1), x)
297
-      } else if(x[whichroot, ncol(x)] != 1) {
298
-        warning("Fitness of wildtype != 1.",
299
-                " Dividing all fitnesses by fitness of wildtype.")
300
-        vwt <- x[whichroot, ncol(x)]
301
-        x[, ncol(x)] <- x[, ncol(x)]/vwt
302
-      }
252
+        whichroot <- which(rowSums(x[, -ncol(x), drop = FALSE]) == 0)
253
+        if(length(whichroot) == 0) {
254
+            warning("No wildtype in the fitness landscape!!! Adding it with fitness 1.")
255
+            x <- rbind(c(rep(0, ncol(x) - 1), 1), x)
256
+        } else if(x[whichroot, ncol(x)] != 1) {
257
+            warning("Fitness of wildtype != 1.",
258
+                    " Dividing all fitnesses by fitness of wildtype.")
259
+            vwt <- x[whichroot, ncol(x)]
260
+            x[, ncol(x)] <- x[, ncol(x)]/vwt
261
+        }
303 262
     }
304 263
 
305 264
     if(any(is.na(x)))
306 265
         stop("NAs in fitness matrix")
307 266
     if(!frequencyDependentFitness) {
308 267
         if(is.data.frame(x)) 
309
-        x <- as.matrix(x)
268
+            x <- as.matrix(x)
310 269
         stopifnot(inherits(x, "matrix"))
311 270
 
312 271
         if(simplify) {
313
-           x <- x[x[, ncol(x)] > min_filter_fitness, , drop = FALSE]  
272
+            x <- x[x[, ncol(x)] > min_filter_fitness, , drop = FALSE]  
314 273
         }
315 274
         class(x) <- c("matrix", "genotype_fitness_matrix")
316 275
     } else { ## frequency-dependent fitness
317
-        ## RDU: what is this for? It converts to numbers.
318
-        ## But it does not work reliably. It does not work when given as "n_"
319
-        ## We will do this at end. in full_FDF_spec
320
-        ## conversionTable_o <- conversionTable(ncol(x) - 1, frequencyType)
321
-        
322
-        ## x[, ncol(x)] <- sapply(x[, ncol(x)],
323
-        ##                        function(x){ findAndReplace(x, conversionTable_o)})
324
-        
325
-    if(frequencyType == "auto"){
326
-      ch <- paste(as.character(x[, ncol(x)]), collapse = "")
327
-      if( grepl("f_", ch, fixed = TRUE) ){
328
-        frequencyType = "rel"
329
-        pattern <- stringr::regex("f_(\\d*_*)*")
330
-        
331
-      } else if ( grepl("n_", ch, fixed = TRUE) ){
332
-        frequencyType = "abs"
333
-        pattern <- stringr::regex("n_(\\d*_*)*")
334 276
         
335
-      } else { stop("No pattern found when frequencyType set to 'auto'") }
336
-        
337
-    } else if(frequencyType == "abs"){
277
+        if(frequencyType == "auto"){
278
+            ch <- paste(as.character(x[, ncol(x)]), collapse = "")
279
+            if( grepl("f_", ch, fixed = TRUE) ){
280
+                frequencyType = "rel"
281
+                pattern <- stringr::regex("f_(\\d*_*)*")
282
+                
283
+            } else if ( grepl("n_", ch, fixed = TRUE) ){
284
+                frequencyType = "abs"
285
+                pattern <- stringr::regex("n_(\\d*_*)*")
286
+                
287
+            } else { stop("No pattern found when frequencyType set to 'auto'") }
288
+            
289
+        } else if(frequencyType == "abs"){
338 290
             pattern <- stringr::regex("n_(\\d*_*)*")
339 291
         } else {
340 292
             pattern <- stringr::regex("f_(\\d*_*)*")
... ...
@@ -352,7 +304,7 @@ to_genotFitness_std <- function(x,
352 304
         }
353 305
 
354 306
     }
355
-  return(x)
307
+    return(x)
356 308
 }
357 309
 
358 310
 ## Deprecated after flfast
... ...
@@ -624,8 +576,6 @@ Magellan_stats <- function(x, max_num_genotypes = 2000,
624 576
                       args = paste(shortarg, logarg, zarg, "-o", fnret, fn),
625 577
                       stdout = NULL)
626 578
     if(short) {
627
-        ## tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[-1])
628
-
629 579
         tmp <- unlist(read.table(fnret, skip = 1, header = TRUE)[c(-1)])
630 580
         ## ## Make names more explicit, but check we have what we think we have
631 581
         ## ## New versions of Magellan produce different output apparently of variable length
... ...
@@ -12,92 +12,26 @@
12 12
 ## License: GPL (>= 2)
13 13
 
14 14
 nem_transitive.reduction <- function(g){
15
-	if (!any(class(g)%in%c("matrix","graphNEL"))) stop("Input must be an adjacency matrix or graphNEL object")
16
-	if(any(class(g) == "graphNEL")){
17
-		g = as(g, "matrix")		
18
-	}
19
-# 	if("Rglpk" %in% loadedNamespaces()){ # constraints müssen einzeln hinzugefügt und jedesmal das ILP gelöst werden. Danach müssen jedesmal die constraints überprüft werden und nicht mehr gebrauchte rausgeschmissen werden
20
-# 		g = abs(g) - diag(diag(g))
21
-# 		mat = matrix(0, ncol=sum(g), nrow=0)
22
-# 		idx = cbind(which(g == 1, arr.ind=T), 1:sum(g))		
23
-# 		for(y in 1:nrow(g)){
24
-# 			for(x in 1:nrow(g)){
25
-# 				if(g[x,y] != 0){
26
-# 					for(j in 1:nrow(g)){
27
-# 						if((g[y,j] != 0) && (g[x,j] != 0)){
28
-# 							mat.tmp = double(sum(g))
29
-# 							mat.tmp[idx[idx[,1] == x & idx[,2] == y, 3]] = 1
30
-# 							mat.tmp[idx[idx[,1] == y & idx[,2] == j, 3]] = 1
31
-# 							mat.tmp[idx[idx[,1] == x & idx[,2] == j, 3]] = -1
32
-# 							mat = rbind(mat, mat.tmp)
33
-# 						}						
34
-# 					}
35
-# 				}
36
-# 			}
37
-# 		}
38
-# 		
39
-# 		solve.problem = function(mat.tmp){
40
-# 			obj = rep(1, NCOL(mat.tmp))
41
-# 			rhs = 2
42
-# 			dir = ">="
43
-# 			sol = Rglpk_solve_LP(obj, mat.tmp, dir, rhs, max=TRUE, types=rep("B", NCOL(mat.tmp)))
44
-# 			print(sol)			
45
-# 			del = idx[which(sol$solution == 0), c(1, 2),drop=F]
46
-# 			for(i in 1:NROW(del)){
47
-# 				g[del[i,1], del[i,2]] = 0
48
-# 			}
49
-# 			g
50
-# 		}
51
-# 		
52
-# 		while(NROW(mat) > 0){
53
-# 			g = solve.problem(mat[1,,drop=F])
54
-# 			i = 1
55
-# 			while(i <= NROW(mat)){
56
-# 				idx.pos = which(mat[i,,drop=F] == 1)
57
-# 				idx.neg = which(mat[i,,drop=F] == -1)
58
-# 				xy = idx[idx[,3] %in% idx.pos, c(1,2)]
59
-# 				z = idx[idx[,3] == idx.neg, c(1,2)]
60
-# 				if(!(g[xy[1,1], xy[1,2]] == 1 & g[xy[2,1], xy[2,2]] == 1 & g[z[1], z[2]] == 1)) # remove resolved constraints
61
-# 					mat = mat[-i,,drop=F]
62
-# 				else
63
-# 					i = i + 1
64
-# 			}
65
-# 		}
66
-# 	}
67
-# 	else{
68
-# 	if(class(g) == "matrix"){		
69
-		# modified algorithm from Sedgewick book: just remove transitive edges instead of inserting them		
15
+    if (!any(class(g)%in%c("matrix","graphNEL"))) stop("Input must be an adjacency matrix or graphNEL object")
16
+    if(any(class(g) == "graphNEL")){
17
+        g = as(g, "matrix")		
18
+    }
19
+
70 20
     g = nem_transitive.closure(g, mat=TRUE) # bug fix: algorithm only works for transitively closed graphs!
71
-		g = g - diag(diag(g))
72
-		type = (g > 1)*1 - (g < 0)*1	
73
-		for(y in 1:nrow(g)){
74
-			for(x in 1:nrow(g)){
75
-				if(g[x,y] != 0){
76
-					for(j in 1:nrow(g)){
77
-						if((g[y,j] != 0) & sign(type[x,j])*sign(type[x,y])*sign(type[y,j]) != -1){ 
78
-						    g[x,j] = 0
79
-						}
80
-					}
81
-				}
82
-			}
83
-		}
84
-# 	}
85
-# 	else{		
86
-# 		nodenames=nodes(g)		
87
-# 		for(y in 1:length(nodes(g))){
88
-# 			edges = edgeL(g)
89
-# 			x = which(sapply(edges,function(l) y %in% unlist(l)))
90
-# 			j = unlist(edges[[y]])
91
-# 			cands = sapply(edges[x], function(e) list(intersect(unlist(e),j)))			
92
-# 			cands = cands[sapply(cands,length) > 0]
93
-# 			if(length(cands) > 0)
94
-# 				for(c in 1:length(cands)){ 				
95
-# 					jj = unlist(cands[c])					
96
-# 					g = removeEdge(rep(names(cands)[c],length(jj)),nodenames[jj],g)
97
-# 				}
98
-# 		}
99
-# 	}	
100
-	g		
21
+    g = g - diag(diag(g))
22
+    type = (g > 1)*1 - (g < 0)*1	
23
+    for(y in 1:nrow(g)){
24
+        for(x in 1:nrow(g)){
25
+            if(g[x,y] != 0){
26
+                for(j in 1:nrow(g)){
27
+                    if((g[y,j] != 0) & sign(type[x,j])*sign(type[x,y])*sign(type[y,j]) != -1){ 
28
+                        g[x,j] = 0
29
+                    }
30
+                }
31
+            }
32
+        }
33
+    }
34
+    g		
101 35
 }
102 36
 
103 37
 
... ...
@@ -107,33 +41,24 @@ nem_transitive.closure <- function(g,mat=FALSE,loops=TRUE){
107 41
     if (!any(class(g)%in%c("graphNEL","matrix"))) stop("Input must be either graphNEL object or adjacency matrix")
108 42
     g <- as(g, "matrix")
109 43
     
110
-    #-- adjacency matrix
111
-#     if (class(g)=="matrix"){
112
-		n <- ncol(g)
113
-		matExpIterativ <- function(x,pow,y=x,z=x,i=1) {
114
-		while(i < pow) {
115
-			z <- z %*% x
116
-			y <- y+z
117
-			i <- i+1
118
-		}
119
-		return(y)
120
-		}
121
-	
122
-		h <- matExpIterativ(g,n)
123
-		h <- (h>0)*1   
124
-		dimnames(h) <- dimnames(g)
125
-		if (!loops) diag(h) <- rep(0,n) else diag(h) <- rep(1,n)
126
-		if (!mat) h <- as(h,"graphNEL")	
127
-#     }
44
+                                        #-- adjacency matrix
45
+                                        #     if (class(g)=="matrix"){
46
+    n <- ncol(g)
47
+    matExpIterativ <- function(x,pow,y=x,z=x,i=1) {
48
+        while(i < pow) {
49
+            z <- z %*% x
50
+            y <- y+z
51
+            i <- i+1
52
+        }
53
+        return(y)
54
+    }
55
+    
56
+    h <- matExpIterativ(g,n)
57
+    h <- (h>0)*1   
58
+    dimnames(h) <- dimnames(g)
59
+    if (!loops) diag(h) <- rep(0,n) else diag(h) <- rep(1,n)
60
+    if (!mat) h <- as(h,"graphNEL")	
61
+                                        #     }
128 62
 
129
-# #     -- graphNEL object
130
-#     if (class(g)=="graphNEL"){
131
-#         tc <- RBGL::transitive.closure(g)    
132
-#         if (loops) tc$edges <- unique(cbind(tc$edges,rbind(tc$nodes,tc$nodes)),MARGIN=2)
133
-# 
134
-#         h <- ftM2graphNEL(ft=t(tc$edges),V=tc$nodes)
135
-#         if (mat) h <- as(h, "matrix")
136
-#     } 
137
-       
138 63
     return(h)
139 64
 }
... ...
@@ -66,11 +66,6 @@ gtm2 <- function(x) {
66 66
     data.frame(cbind(nice.vector.eo(x, ","), x), stringsAsFactors = TRUE)
67 67
 }
68 68
 
69
-## nice.vector.eo <- function(z, sep) {
70
-##     ## with epistasis, maybe we want sorted?
71
-##     setdiff(unlist(lapply(strsplit(z, " "),
72
-##                                     function(u) strsplit(u, sep))), "")
73
-## }
74 69
 
75 70
 nice.vector.eo <- function(z, sep, rm.sign = FALSE) {
76 71
     ## with epistasis, maybe we want sorted?
... ...
@@ -176,40 +171,18 @@ list.of.deps <- function(x) {
176 171
 
177 172
 to.long.rt <- function(rt, idm) {
178 173
     ## We now do this inconditionally, so that we do not need to use the
179
-    ## "stringsAsFactors = FALSE". This is now done before
180
-    ## if(is.numeric(rt$parent))
181
-    ##     rt$parent <- as.character(rt$parent)
182
-    ## if(is.numeric(rt$child))
183
-    ##     rt$child <- as.character(rt$child)
184 174
    
185 175
     
186 176
     if(!("Root" %in% rt$parent))
187 177
         stop("Root must be one parent node")
188 178
 
189
-    ## rt$parent <- unlist(lapply(rt$parent, nice.string))
190
-    ## rt$child <- unlist(lapply(rt$child, nice.string))
191 179
    
192 180
     srt <- rt[order(rt$child), ]
193 181
 
194
-    ## Not relevant if we allow non-numeric names
195
-    ## all.child.genes <- as.integer(
196
-    ##     unlist(lapply(rt[, 2],
197
-    ##                   function(x) strsplit(x, ","))))
198
-    ## ## check all childs
199
-    ## if(!identical(sort(unique(all.child.genes)),
200
-    ##               seq.int(max(all.child.genes))))
201
-    ##     stop("Not all children present")
202
-    long.rt <- lapply(split(srt, srt$child), list.of.deps)
203
-
204
-    ## geneModule <- gene.to.module(srt)
205
-    ## idm <- seq.int(length(names(long.rt)))
206
-    ## names(idm) <- names(long.rt)
207
-    ## idm <- c("0" = 0L, idm)
208
-    ## geneModule$ModuleNumID <- idm[geneModule[, "Module"]]
182
+     long.rt <- lapply(split(srt, srt$child), list.of.deps)
209 183
 
184
+ 
210 185
     ## idm is just a look up table for the id of the module
211
-    ## idm <- unique(geneModule$ModuleNumID)
212
-    ## names(idm) <- unique(geneModule$Module)
213 186
     
214 187
     ## add integer IDs
215 188
     addIntID <- function(z, idm) {
... ...
@@ -240,14 +213,7 @@ to.long.rt <- function(rt, idm) {
240 213
     }
241 214
     long.rt <- lapply(long.rt, function(x) addIntID(x, idm = idm))
242 215
    
243
-    ## if(verbosity >= 4) {
244
-    ##     message(paste("Number of drivers: ",
245
-    ##                   length(unique(geneModule[, "Gene"]))))
246
-    ##     message(paste("Number of modules: ",
247
-    ##                   length(unique(geneModule[, "Module"]))))
248
-    ## }
249 216
     return(long.rt)
250
-    ## return(list(long.rt = long.rt, geneModule = geneModule))
251 217
 }
252 218
 
253 219
 
... ...
@@ -306,28 +272,10 @@ to.long.epist.order <- function(epor, sep, rm.sign = FALSE) {
306 272
     ## just vectors for now
307 273
     long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
308 274
                 names(epor), epor)
309
-    ## if(is.vector(epor))
310
-    ##     long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
311
-    ##                    names(epor), epor)
312
-    ## else if(is.data.frame(epor)) 
313
-    ##     long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
314
-    ##                 as.character(epor$ids),
315
-    ##                 epor$s)
316 275
     names(long) <- NULL
317 276
     return(long)
318 277
 }
319 278
 
320
-## addIntID.epist.order <- function(z, idm, sort) {
321
-##     z$NumID <- idm[z$ids]
322
-##     if(sort) {
323
-##         ## essential for epistasis, but never do for order effects
324
-##         o <- order(z$NumID)
325
-##         z$NumID <- z$NumID[o]
326
-##         z$ids <- z$ids[o]
327
-##     }
328
-##     return(z)
329
-## }
330
-
331 279
 
332 280
 addIntID.epist.order <- function(z, idm, sort, sign) {
333 281
     if( sort && (!sign))
... ...
@@ -414,31 +362,6 @@ getGeneIDNum <- function(geneModule, geneNoInt, fitnessLandscape_gene_id,
414 362
     )
415 363
 }
416 364
 
417
-## Next two in crate_flvars_fitvars.
418
-
419
-## ## generate fitnesLanscapeVariables for FDF
420
-## ## fitness landscape as data frame with cols as genes,
421
-## ## and last column is genotype
422
-## create_flvars <- function(x, frequencyType) {
423
-##     x <- x[, -ncol(x), drop = FALSE]
424
-##     pasted <- apply(x, 1, function(z) paste(which(z == 1), collapse = "_"))
425
-##     npasted <- apply(x, 1, function(z) paste(colnames(x)[which(z == 1)], collapse = "_"))
426
-##     if(frequencyType == "abs") {
427
-##         tmp <- paste0("n_", pasted)
428
-##     } else {
429
-##         tmp <- paste0("f_", pasted)
430
-##     }
431
-##     names(tmp) <- npasted
432
-##     return(tmp)
433
-## }
434
-
435
-## ## fitness specification with letters -> fitness specification with numbers
436
-## names_to_nums_fitness <- function(fitness, genotFitness) {
437
-##     from_subst_pattern <- colnames(genotFitness[, -ncol(genotFitness)])
438
-##     to_subst_pattern <- as.character(1:(ncol(genotFitness) - 1))
439
-##     names(to_subst_pattern) <- from_subst_pattern
440
-##     return(stringr::str_replace_all(fitness, to_subst_pattern))
441
-## }
442 365
 
443 366
 ## genotFitnes and frequency type -> fitnesLanscapeVariables for FDF and
444 367
 ##          fitness with numbers, not names
... ...
@@ -509,20 +432,9 @@ all_orders_fv <- function(x, prefix, prefixre) {
509 432
     return(out)
510 433
 }
511 434
 
512
-## ## character vector, named replacement -> replaced vector
513
-## ## named_replace: names are the (fixed) pattern, value the replacement
514
-## ## stringr::str_replace_all seems too smart and does not respect order
515
-## my_mgsub <- function(x, named_replace) {
516
-##     nn <- names(named_replace)
517
-##     xr <- x
518
-##     for(i in 1:length(named_replace)) {
519
-##         xr <- gsub(nn[i], named_replace[i], xr, fixed = TRUE)
520
-##     }
521
-##     return(xr)
522
-## }
523 435
 
524 436
 
525
-##New function
437
+
526 438
 fVariablesN <- function (g, frequencyType) {
527 439
 
528 440
   if (is.null(g) | g == 0)
... ...
@@ -547,52 +459,6 @@ fVariablesN <- function (g, frequencyType) {
547 459
   return (fsVector)
548 460
 }
549 461
 
550
-fVariablesL <- function (g, frequencyType) {
551
-
552
-  if (is.null(g) | g == 0)
553
-    stop("Number of genes must be integer > 0")
554
-
555
-  if(g > length(LETTERS))
556
-    stop(paste0("Number of genes must be < length(LETTERS).",
557
-                " Please specify variables with numbers"))
558
-
559
-  combinationsList <- list()
560
-  for (i in 0:g) {
561
-    combinationsList <- append(combinationsList,
562
-                               combn(LETTERS[1:g], i, list, simplify = TRUE))
563
-  }
564
-
565
-  if (frequencyType == "abs"){
566
-    fsVector <-sapply(sapply(combinationsList,
567
-                             function(x) paste0(x, collapse = "_")),
568
-                      function(x) paste0("n_", x))
569
-  }else{
570
-    fsVector <-sapply(sapply(combinationsList,
571
-                             function(x) paste0(x, collapse = "_")),
572
-                      function(x) paste0("f_", x))
573
-  }
574
-
575
-  return (fsVector)
576
-}
577
-
578
-## ## Assuming we are using the full fitness landscapes (i.e., none of
579
-## ## setting genotypes with 0 fitness are absent from the table)
580
-## conversionTable <- function(g, frequencyType){
581
-##   df <- data.frame(let = fVariablesL(g, frequencyType)[-1],
582
-##                    num = fVariablesN(g, frequencyType)[-1],
583
-##                    stringsAsFactors = FALSE)
584
-##   return (df)
585
-## }
586
-
587
-## findAndReplace <- function(str, conversionTable_input){
588
-
589
-##   pattern <- rev(setNames(as.character(conversionTable_input$num),
590
-##                       conversionTable_input$let))
591
-
592
-##   str <- stringr::str_replace_all(string = str,
593
-##                                   pattern = pattern)
594
-##   return(str)
595
-## }
596 462
 
597 463
 allFitnessORMutatorEffects <- function(rT = NULL,
598 464
                                        epistasis = NULL,
... ...
@@ -839,7 +705,7 @@ allFitnessORMutatorEffects <- function(rT = NULL,
839 705
     } else if(calledBy == "allMutatorEffects") {
840 706
       class(out) <- c("mutatorEffects")
841 707
     }
842
-  } else {
708
+  } else { ## Frequency-dependent fitness
843 709
 
844 710
     if(is.null(genotFitness)) {
845 711
       #genotFitness <- matrix(NA, nrow = 0, ncol = 1)
... ...
@@ -895,22 +761,41 @@ allFitnessORMutatorEffects <- function(rT = NULL,
895 761
                                       stringsAsFactors = FALSE)
896 762
       ## Create, for the user, a single data frame with everything.
897 763
       ## This is what C++ should consume
898
-
899
-   
900 764
      
901 765
       ## This ought to allow to pass fitness spec as letters. Preserve original
902 766
       Fitness_original_as_letters <- fitnessLandscape_df$Fitness
903 767
       fitnessLandscape_df$Fitness <- Fitness_as_fvars
904
-      
905
-      full_FDF_spec <- cbind(genotFitness[, -ncol(genotFitness)]
906
-                 , Genotype_letters = genotype_letterslabel(genotFitness[, -ncol(genotFitness)])
907
-                 , Genotype_fvars = fitnessLandscapeVariables ## used in C++
908
-                 , Fitness_as_fvars = Fitness_as_fvars
909
-                 , Fitness_as_letters = Fitness_original_as_letters
910
-                   )
768
+
769
+      full_FDF_spec <-
770
+          cbind(genotFitness[, -ncol(genotFitness)]
771
+              , Genotype_as_numbers = fitnessLandscape_df$Genotype
772
+              , Genotype_as_letters = genotype_letterslabel(genotFitness[, -ncol(genotFitness)])
773
+              , Genotype_as_fvars = fitnessLandscapeVariables ## used in C++
774
+              , Fitness_as_fvars = Fitness_as_fvars
775
+              , Fitness_as_letters = Fitness_original_as_letters
776
+                )
911 777
       rownames(full_FDF_spec) <- 1:nrow(full_FDF_spec)
912 778
       
913
-    out <- list(long.rt = list(),
779
+      ## fitnessLanscape and fitnessLandscape_df are now redundant given
780
+      ## full_FDF_spec. Remove them later. In the mean time, ensure a
781
+      ## single canonical object used.
782
+
783
+      rm(fitnessLandscape_df)
784
+      suppressWarnings(try(rm(fitnessLandscape), silent = TRUE))
785
+      rm(fitnessLandscapeVariables)
786
+      rm(Fitness_as_fvars)
787
+      rm(Fitness_original_as_letters)
788
+      
789
+      fitnessLandscape <- full_FDF_spec[, c(fitnessLandscape_gene_id$Gene,
790
+                                            "Fitness_as_fvars")]
791
+      colnames(fitnessLandscape)[ncol(fitnessLandscape)] <- "Fitness"
792
+      
793
+      fitnessLandscape_df <- full_FDF_spec[, c("Genotype_as_numbers",
794
+                                               "Fitness_as_fvars")]
795
+      colnames(fitnessLandscape_df) <- c("Genotype", "Fitness")
796
+      
797
+      
798
+      out <- list(long.rt = list(),
914 799
                 long.epistasis = list(),
915 800
                 long.orderEffects = list(),
916 801
                 long.geneNoInt = data.frame(),
... ...
@@ -923,10 +808,10 @@ allFitnessORMutatorEffects <- function(rT = NULL,
923 808
                 epistasis = NULL,
924 809
                 orderEffects = NULL,
925 810
                 noIntGenes = NULL,
926
-                fitnessLandscape = genotFitness,
927
-                fitnessLandscape_df = fitnessLandscape_df,
928
-                fitnessLandscape_gene_id = fitnessLandscape_gene_id,
929
-                fitnessLandscapeVariables = fitnessLandscapeVariables,
811
+                fitnessLandscape = genotFitness, ## redundant
812
+                fitnessLandscape_df = fitnessLandscape_df, ## redundant
813
+                fitnessLandscape_gene_id = fitnessLandscape_gene_id, 
814
+                ## fitnessLandscapeVariables = NULL, ## now part of full_FDF_spec
930 815
                 frequencyDependentFitness = frequencyDependentFitness,
931 816
                 frequencyType = frequencyType,
932 817
                 full_FDF_spec = full_FDF_spec
... ...
@@ -938,1849 +823,2098 @@ allFitnessORMutatorEffects <- function(rT = NULL,
938 823
   return(out)
939 824
 }
940 825
 
941
-## Former version, with fitness landscape
942
-## allFitnessORMutatorEffects <- function(rT = NULL,
943
-##                                        epistasis = NULL,
944
-##                                        orderEffects = NULL,
945
-##                                        noIntGenes = NULL,
946
-##                                        geneToModule = NULL,
947
-##                                        drvNames = NULL,
948
-##                                        keepInput = TRUE,
949
-##                                        ## refFE = NULL,
950
-##                                        calledBy = NULL) {
951
-##     ## From allFitnessEffects. Generalized so we deal with Fitness
952
-##     ## and mutator.
826
+allFitnessEffects <- function(rT = NULL,
827
+                              epistasis = NULL,
828
+                              orderEffects = NULL,
829
+                              noIntGenes = NULL,
830
+                              geneToModule = NULL,
831
+                              drvNames = NULL,
832
+                              genotFitness = NULL,
833
+                              keepInput = TRUE,
834
+                              frequencyDependentFitness = FALSE,
835
+                              frequencyType = NA) {
836
+                              #spPopSizes = NULL) {
953 837
 
954
-##     ## restrictions: the usual rt
838
+    if(!frequencyDependentFitness) {
839
+        
840
+        if(!is.na(frequencyType)){
841
+            warning("frequencyType set to NA")
842
+        }
843
+        ## this is a kludge, but we must pass something not NA and not NULL
844
+        ## to the C++ code
845
+        frequencyType = "freq_dep_not_used"
955 846
 
956
-##     ## epistasis: as it says, with the ":"
847
+    if(!is.null(genotFitness)) {
848
+      if(!is.null(rT) || !is.null(epistasis) ||
849
+         !is.null(orderEffects) || !is.null(noIntGenes) ||
850
+         !is.null(geneToModule)) {
851
+        stop("You have a non-null genotFitness.",
852
+             " If you pass the complete genotype to fitness mapping",
853
+             " you cannot pass any of rT, epistasis, orderEffects",
854
+             " noIntGenes or geneToModule.")
855
+      }
957 856
 
958
-##     ## orderEffects: the ">"
857
+      genotFitness_std <- to_genotFitness_std(genotFitness,
858
+                                              frequencyDependentFitness = FALSE,
859
+                                              frequencyType = frequencyType,
860
+                                              simplify = TRUE)
861
+      ## epistasis <- from_genotype_fitness(genotFitness)
862
+    } else {
863
+      genotFitness_std <- NULL
864
+    }
865
+    allFitnessORMutatorEffects(
866
+      rT = rT,
867
+      epistasis = epistasis,
868
+      orderEffects = orderEffects,
869
+      noIntGenes = noIntGenes,
870
+      geneToModule = geneToModule,
871
+      drvNames = drvNames,
872
+      keepInput = keepInput,
873
+      genotFitness = genotFitness_std,
874
+      calledBy = "allFitnessEffects",
875
+      frequencyDependentFitness = FALSE,
876
+      frequencyType = frequencyType)
877
+      #spPopSizes = spPopSizes)
959 878
 
960
-##     ## All of the above can be genes or can be modules (if you pass a
961
-##     ## geneToModule)
879
+  } else {
962 880
 
963
-##     ## rest: rest of genes, with fitness
881
+    if(!(frequencyType %in% c('abs', 'rel', 'auto'))){
882
+      #set frequencyType = "auto" in case you did not specify 'rel' or 'abs'
883
+      frequencyType = "auto"
884
+      message("frequencyType set to 'auto'")
885
+    }
964 886
 
887
+    if(is.null(genotFitness)) {
888
+      stop("You have a null genotFitness in a frequency dependent fitness situation.")
889
+    } else {
890
+      genotFitness_std <- to_genotFitness_std(genotFitness,
891
+                                              frequencyDependentFitness = TRUE,
892
+                                              frequencyType = frequencyType,
893
+                                              simplify = TRUE)
894
+      allFitnessORMutatorEffects(
895
+        rT = rT,
896
+        epistasis = epistasis,
897
+        orderEffects = orderEffects,
898
+        noIntGenes = noIntGenes,
899
+        geneToModule = geneToModule,
900
+        drvNames = drvNames,
901
+        keepInput = keepInput,
902
+        genotFitness = genotFitness_std,
903
+        calledBy = "allFitnessEffects",
904
+        frequencyDependentFitness = TRUE,
905
+        frequencyType = frequencyType)
906
+        #spPopSizes = spPopSizes)
907
+    }
908
+  }
909
+}
965 910
 
966
-##     ## For epistasis and order effects we create the output object but
967
-##     ## missing the numeric ids of genes. With rT we do it in one go, as we
968
-##     ## already know the mapping of genes to numeric ids. We could do the
969
-##     ## same in epistasis and order, but we would be splitting twice
970
-##     ## (whereas for rT extracting the names is very simple).
971 911
 
972
-##     ## called appropriately?
973
-##     if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
974
-##         stop("How did you call this function?. Bug.")
975 912
 
976
-##     if(calledBy == "allMutatorEffects") {
977
-##         ## very paranoid check
978
-##         if( !is.null(rT) || !is.null(orderEffects) || !is.null(drvNames))
979
-##             stop("allMutatorEffects called with forbidden arguments.",
980
-##                  "Is this an attempt to subvert the function?")
981
-##     }
913
+evalGenotypeORMut <- function(genotype,
914
+                              fmEffects,
915
+                              spPopSizes = spPopSizes,
916
+                              verbose = FALSE,
917
+                              echo = FALSE,
918
+                              model = "",
919
+                              calledBy_= NULL,
920
+                              currentTime = currentTime) {
921
+    ## genotype can be a vector of integers, that are the exact same in
922
+    ## the table of fmEffects or a vector of strings, or a vector (a
923
+    ## string) with genes separated by "," or ">"
982 924
 
983
-##     rtNames <- NULL
984
-##     epiNames <- NULL
985
-##     orNames <- NULL
986
-##     if(!is.null(rT)) {
987
-##         ## This is really ugly, but to prevent the stringsAsFactors I need it here:
988
-##         rT$parent <- as.character(rT$parent)
989
-##         rT$child <- as.character(rT$child)
990
-##         rT$typeDep <- as.character(rT$typeDep)
991
-##         rtNames <- unique(c(rT$parent, rT$child))
992
-##     }
993
-##     if(!is.null(epistasis)) {
994
-##         long.epistasis <- to.long.epist.order(epistasis, ":")
995
-##         ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
996
-##         ## deal with the possible negative signs
997
-##         epiNames <- setdiff(unique(
998
-##             unlist(lapply(long.epistasis,
999
-##                           function(x) lapply(x$ids,
1000
-##                                              function(z) strsplit(z, "^-"))))),
1001
-##                             "")
1002
-##     } else {
1003
-##         long.epistasis <- list()
1004
-##     }
1005
-##     if(!is.null(orderEffects)) {
1006
-##         long.orderEffects <- to.long.epist.order(orderEffects, ">")
1007
-##         orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
1008
-##     } else {
1009
-##         long.orderEffects <- list()
1010
-##     }
1011
-##     allModuleNames <- unique(c(rtNames, epiNames, orNames))
1012
-##     if(is.null(geneToModule)) {
1013
-##         gMOneToOne <- TRUE
1014
-##         geneToModule <- geneModuleNull(allModuleNames)
1015
-##     } else {
1016
-##         gMOneToOne <- FALSE
1017
-##         if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
1018
-##             stop(paste("Some values in geneToModule not present in any of",
1019
-##                        " rT, epistasis, or order effects"))
1020
-##         if(any(is.na(match(allModuleNames, names(geneToModule)))))
1021
-##             stop(paste("Some values in rT, epistasis, ",
1022
-##                        "or order effects not in geneToModule"))
1023
-##     }
1024
-##     geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
925
+    if( !(calledBy_ %in% c("evalGenotype", "evalGenotypeMut") ))
926
+        stop("How did you call this function?. Bug.")
1025 927
 
1026
-##     idm <- unique(geneModule$ModuleNumID)
1027
-##     names(idm) <- unique(geneModule$Module)
928
+    ## fmEffects could be a mutator effect
929
+    if(!exists("fitnessLandscape_gene_id", where = fmEffects)) {
930
+        fmEffects$fitnessLandscape_df <- data.frame()
931
+        fmEffects$fitnessLandscape_gene_id <- data.frame()
932
+    }
1028 933
 
1029
-##     if(!is.null(rT)) {
1030
-##         checkRT(rT)
1031
-##         long.rt <- to.long.rt(rT, idm)
1032
-##     } else {
1033
-##         long.rt <- list() ## yes, we want an object of length 0
1034
-##     }
934
+    if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
935
+        (nrow(fmEffects$fitnessLandscape_df) > 0)) {
936
+        warning("Bozic model passing a fitness landscape will not work",
937
+                " for now.")
938
+    }
1035 939
 
1036
-##     ## Append the numeric ids to epistasis and order
1037
-##     if(!is.null(epistasis)) {
1038
-##         long.epistasis <- lapply(long.epistasis,
1039
-##                                  function(x)
1040
-##                                      addIntID.epist.order(x, idm,
1041
-##                                                           sort = TRUE,
1042
-##                                                           sign = TRUE))
1043
-##     }
1044
-##     if(!is.null(orderEffects)) {
1045
-##         long.orderEffects <- lapply(long.orderEffects,
1046
-##                                     function(x)
1047
-##                                         addIntID.epist.order(x, idm,
1048
-##                                                              sort = FALSE,
1049
-##                                                              sign = FALSE))
1050
-##     }
1051
-
1052
-##     if(!is.null(noIntGenes)) {
1053
-##         if(inherits(noIntGenes, "character")) {
1054
-##             wm <- paste("noIntGenes is a character vector.",
1055
-##                         "This is probably not what you want, and will",
1056
-##                         "likely result in an error downstream.",
1057
-##                         "You can get messages like",
1058
-##                         " 'not compatible with requested type', and others.",
1059
-##                         "We are stopping.")
1060
-##             stop(wm)
1061
-##         }
1062
-
1063
-##         mg <- max(geneModule[, "GeneNumID"])
1064
-##         gnum <- seq_along(noIntGenes) + mg
1065
-##         if(!is.null(names(noIntGenes))) {
1066
-##             ng <- names(noIntGenes)
1067
-##             if( grepl(",", ng, fixed = TRUE) || grepl(">", ng, fixed = TRUE)
1068
-##                 || grepl(":", ng, fixed = TRUE))
1069
-##                 stop("The name of some noIntGenes contain a ',' or a '>' or a ':'")
1070
-##             if(any(ng %in% geneModule[, "Gene"] ))
1071
-##                 stop("A gene in noIntGenes also present in the other terms")
1072
-##             if(any(duplicated(ng)))
1073
-##                 stop("Duplicated gene names in geneNoInt")
1074
-##             if(any(is.na(ng)))
1075
-##                 stop("In noIntGenes some genes have names, some don't.",
1076
-##                      " Name all of them, or name none of them.")
1077
-##         } else {
1078
-##             ng <- gnum
1079
-##         }
1080
-##         geneNoInt <- data.frame(Gene = as.character(ng),
1081
-##                                 GeneNumID = gnum,
1082
-##                                 s = noIntGenes,
1083
-##                                 stringsAsFactors = FALSE)
1084
-##     } else {
1085
-##         geneNoInt <- data.frame()
1086
-##     }
1087
-
1088
-##     if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
1089
-##              nrow(geneNoInt)) == 0)
1090
-##         stop("You have specified nothing!")
1091
-
1092
-##     if(calledBy == "allFitnessEffects") {
1093
-##         if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
1094
-##             graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
1095
-##         } else {
1096
-##             graphE <- NULL
1097
-##         }
1098
-##     } else {
1099
-##         graphE <- NULL
1100
-##     }
1101
-##     if(!is.null(drvNames)) {
1102
-##         drv <- unique(getGeneIDNum(geneModule, geneNoInt, drvNames))
1103
-##         ## drivers should never be in the geneNoInt; Why!!!???
1104
-##         ## Catch the problem. This is an overkill,
1105
-##         ## so since we catch the issue, we could leave the geneNoInt. But
1106
-##         ## that should not be there in this call.
1107
-##         ## if(any(drvNames %in% geneNoInt$Gene)) {
1108
-##         ##     stop(paste("At least one gene in drvNames is a geneNoInt gene.",
1109
-##         ##                "That is not allowed.",
1110
-##         ##                "If that gene is a driver, pass it as gene in the epistasis",
1111
-##         ##                "component."))
1112
-##         ## }
1113
-##         ## drv <- getGeneIDNum(geneModule, NULL, drvNames)
1114
-##     } else {
1115
-##         ## we used to have this default
1116
-##         ## drv <- geneModule$GeneNumID[-1]
1117
-##         drv <- vector(mode = "integer", length = 0L)
1118
-##     }
1119
-
1120
-##     if(!keepInput) {
1121
-##         rT <- epistasis <- orderEffects <- noIntGenes <- NULL
1122
-##     }
1123
-##     out <- list(long.rt = long.rt,
1124
-##                 long.epistasis = long.epistasis,
1125
-##                 long.orderEffects = long.orderEffects,
1126
-##                 long.geneNoInt = geneNoInt,
1127
-##                 geneModule = geneModule,
1128
-##                 gMOneToOne = gMOneToOne,
1129
-##                 geneToModule = geneToModule,
1130
-##                 graph = graphE,
1131
-##                 drv = drv,
1132
-##                 rT = rT,
1133
-##                 epistasis = epistasis,
1134
-##                 orderEffects = orderEffects,
1135
-##                 noIntGenes = noIntGenes
1136
-##                 )
1137
-##     if(calledBy == "allFitnessEffects") {
1138
-##         class(out) <- c("fitnessEffects")
1139
-##     } else if(calledBy == "allMutatorEffects") {
1140
-##         class(out) <- c("mutatorEffects")
1141
-##     }
1142
-##     return(out)
1143
-## }
940
+    ## This will avoid errors is evalRGenotype where spPopSizes = NULL  
941
+    if (!fmEffects$frequencyDependentFitness) {
942
+        spPopSizes = 0
943
+    } else {
944
+        spPopSizes <- match_spPopSizes(spPopSizes, fmEffects)
945
+    }
1144 946
 
1145
-allFitnessEffects <- function(rT = NULL,
1146
-                              epistasis = NULL,
1147
-                              orderEffects = NULL,
1148
-                              noIntGenes = NULL,
1149
-                              geneToModule = NULL,
1150
-                              drvNames = NULL,
1151
-                              genotFitness = NULL,
1152
-                              keepInput = TRUE,
1153
-                              frequencyDependentFitness = FALSE,
1154
-                              frequencyType = NA) {
1155
-                              #spPopSizes = NULL) {
947
+    if(echo)
948
+        cat(paste("Genotype: ", genotype))
1156 949
 
1157
-    if(!frequencyDependentFitness) {
950
+    if(is.character(genotype)) {
951
+        if(length(grep(">", genotype))) {
952
+            genotype <- nice.vector.eo(genotype, ">")
953
+        } else if(length(grep(",", genotype))) {
954
+            genotype <- nice.vector.eo(genotype, ",")
955
+        }
956
+        all.g.nums <- c(fmEffects$geneModule$GeneNumID,
957
+                        fmEffects$long.geneNoInt$GeneNumID,
958
+                        fmEffects$fitnessLandscape_gene_id$GeneNumID)
959
+        all.g.names <- c(fmEffects$geneModule$Gene,
960
+                         fmEffects$long.geneNoInt$Gene,
961
+                         fmEffects$fitnessLandscape_gene_id$Gene)
962
+        if((length(genotype) == 1) && (genotype %in% c("", "WT") )) {
963
+            genotype <- 0
964
+        } else {
965
+            genotype <- all.g.nums[match(genotype, all.g.names)]
966
+        }
967
+    } else {
968
+        all.g.nums <- c(fmEffects$geneModule$GeneNumID,
969
+                        fmEffects$long.geneNoInt$GeneNumID,
970
+                        fmEffects$fitnessLandscape_gene_id$GeneNumID)
1158 971
         
1159
-        if(!is.na(frequencyType)){
1160
-            warning("frequencyType set to NA")
972
+        if(!all(sapply(genotype,  function(x) x %in% all.g.nums))){
973
+            stop("Genotype as vector of numbers contains genes not in fitnessEffects/mutatorEffects.")
1161 974
         }
1162
-        ## this is a kludge, but we must pass something not NA and not NULL
1163
-        ## to the C++ code
1164
-        frequencyType = "freq_dep_not_used"
1165
-
1166
-    if(!is.null(genotFitness)) {
1167
-      if(!is.null(rT) || !is.null(epistasis) ||
1168
-         !is.null(orderEffects) || !is.null(noIntGenes) ||
1169
-         !is.null(geneToModule)) {
1170
-        stop("You have a non-null genotFitness.",
1171
-             " If you pass the complete genotype to fitness mapping",
1172
-             " you cannot pass any of rT, epistasis, orderEffects",
1173
-             " noIntGenes or geneToModule.")
1174
-      }
975
+    }
1175 976
 
1176
-      genotFitness_std <- to_genotFitness_std(genotFitness,
1177
-                                              frequencyDependentFitness = FALSE,
1178
-                                              frequencyType = frequencyType,
1179
-                                              simplify = TRUE)
1180
-      ## epistasis <- from_genotype_fitness(genotFitness)
977
+    ## FIXME: About the WT: in fmEffects$geneModule$Gene we have "Root". We might
978
+    ## want to refactor that and turn all those to WT. It would avoid the need for
979
+    ## the if above to catch the WT or ""
980
+    
981
+    if(!fmEffects$frequencyDependentFitness){
982
+        
983
+        if( any(is.na(genotype)) ){
984
+            stop("Genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
985
+        }
986
+        
987
+        if((!length(genotype))){
988
+            stop("Genotypes must have at least one mutated gene")
989
+        }
990
+        if(any(genotype < 0)) {
991
+            stop(paste("Genotypes cannot contain negative values.",
992
+                       "If you see this message, you found a bug."))
993
+        }
994
+        if(length(genotype) == 1 && genotype == 0){
995
+            ## We do not handle WT for now
996
+            stop("Genotype cannot be 0. We do not handle WT on its own in non-freq-dep: its fitness is 1 by decree.")
997
+        }
998
+        
999
+        if(any(genotype == 0)){
1000
+            stop("Genotype cannot contain any 0.")
1001
+        }
1002
+        
1181 1003
     } else {
1182
-      genotFitness_std <- NULL
1004
+        if(length(genotype) == 1 && is.na(genotype)){
1005
+            stop("Genotype contains NA or a gene not in fitnessEffects/mutatorEffects")
1006
+        } else if(length(genotype) == 1 && genotype == 0) {
1007
+            ## for WT
1008
+            genotype <- vector(mode = "integer", length = 0L)
1009
+        } else if(length(genotype) > 1){
1010
+            if( any(is.na(genotype)) ){
1011
+                stop("Genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
1012
+            }
1013
+            if(any(genotype == 0)){
1014
+                stop("Genotype cannot contain any 0 if its length > 1")
1015
+            }
1016
+        }
1183 1017
     }
1184
-    allFitnessORMutatorEffects(
1185
-      rT = rT,
1186
-      epistasis = epistasis,
1187
-      orderEffects = orderEffects,
1188
-      noIntGenes = noIntGenes,
1189
-      geneToModule = geneToModule,
1190
-      drvNames = drvNames,
1191
-      keepInput = keepInput,
1192
-      genotFitness = genotFitness_std,
1193
-      calledBy = "allFitnessEffects",
1194
-      frequencyDependentFitness = FALSE,
1195
-      frequencyType = frequencyType)
1196
-      #spPopSizes = spPopSizes)
1018
+    
1019
+    if(model %in% c("Bozic", "bozic1", "bozic2") )
1020
+        prodNeg <- TRUE
1021
+    else
1022
+        prodNeg <- FALSE
1197 1023
 
1198
-  } else {
1024
+    ff <- evalRGenotype(rG = genotype,
1025
+                        rFE = fmEffects,
1026
+                        spPop = spPopSizes,
1027
+                        verbose = verbose,
1028
+                        prodNeg = prodNeg,
1029
+                        calledBy_ = calledBy_,
1030
+                        currentTime = currentTime)
1031
+
1032
+    if(echo) {
1033
+        if(calledBy_ == "evalGenotype") {
1034
+            if(!prodNeg)
1035
+                cat(" Fitness: ", ff, "\n")
1036
+            else
1037
+                cat(" Death rate: ", ff, "\n")
1038
+        } else if(calledBy_ == "evalGenotypeMut") {
1039
+            cat(" Mutation rate product :", ff, "\n")
1040
+        }
1199 1041
 
1200
-    if(!(frequencyType %in% c('abs', 'rel', 'auto'))){
1201
-      #set frequencyType = "auto" in case you did not specify 'rel' or 'abs'
1202
-      frequencyType = "auto"
1203
-      message("frequencyType set to 'auto'")
1204 1042
     }
1205 1043
 
1206
-    if(is.null(genotFitness)) {
1207
-      stop("You have a null genotFitness in a frequency dependent fitness situation.")
1208
-    } else {
1209
-      genotFitness_std <- to_genotFitness_std(genotFitness,
1210
-                                              frequencyDependentFitness = TRUE,
1211
-                                              frequencyType = frequencyType,
1212
-                                              simplify = TRUE)
1213
-      allFitnessORMutatorEffects(
1214
-        rT = rT,
1215
-        epistasis = epistasis,
1216
-        orderEffects = orderEffects,
1217
-        noIntGenes = noIntGenes,
1218
-        geneToModule = geneToModule,
1219
-        drvNames = drvNames,
1220
-        keepInput = keepInput,
1221
-        genotFitness = genotFitness_std,
1222
-        calledBy = "allFitnessEffects",
1223
-        frequencyDependentFitness = TRUE,
1224
-        frequencyType = frequencyType)
1225
-        #spPopSizes = spPopSizes)
1226
-    }
1227
-  }
1044
+    return(ff)
1228 1045
 }
1229 1046
 
1230
-## Former version
1231
-## allFitnessEffects <- function(rT = NULL,
1232
-##                               epistasis = NULL,
1233
-##                               orderEffects = NULL,
1234
-##                               noIntGenes = NULL,
1235
-##                               geneToModule = NULL,
1236
-##                               drvNames = NULL,
1237
-##                               genotFitness = NULL,
1238
-##                               keepInput = TRUE) {
1047
+evalGenotype <- function(genotype,
1048
+                         fitnessEffects,
1049
+                         spPopSizes = NULL,
1050
+                         verbose = FALSE,
1051
+                         echo = FALSE,
1052
+                         model = "",
1053
+                         currentTime = 0
1054
+                         ) {
1055
+  
1056
+    if(inherits(fitnessEffects, "mutatorEffects"))
1057
+         stop("You are trying to get the fitness of a mutator specification. ",
1058
+             "You did not pass an object of class fitnessEffects.")
1239 1059
 
1240
-##     if(!is.null(genotFitness)) {
1241
-##         if(!is.null(rT) || !is.null(epistasis) ||
1242
-##            !is.null(orderEffects) || !is.null(noIntGenes) ||
1243
-##            !is.null(geneToModule)) {
1244
-##             stop("You have a non-null genotFitness.",
1245
-##                  " If you pass the complete genotype to fitness mapping",
1246
-##                  " you cannot pass any of rT, epistasis, orderEffects",
1247
-##                  " noIntGenes or geneToModule.")
1248
-##         }
1249
-##         epistasis <- from_genotype_fitness(genotFitness)
1250
-##     }
1251
-##     allFitnessORMutatorEffects(
1252
-##         rT = rT,
1253
-##         epistasis = epistasis,
1254
-##         orderEffects = orderEffects,
1255
-##         noIntGenes = noIntGenes,
1256
-##         geneToModule = geneToModule,
1257
-##         drvNames = drvNames,
1258
-##         keepInput = keepInput,
1259
-##         calledBy = "allFitnessEffects")
1260
-## }
1060
+   if (fitnessEffects$frequencyDependentFitness) {
1061
+      if (is.null(spPopSizes))
1062
+      stop("You have a NULL spPopSizes")
1063
+    if (!(length(spPopSizes) == nrow(fitnessEffects$fitnessLandscape)))
1064
+      stop("spPopSizes must be as long as number of genotypes")
1065
+   }
1261 1066
 
1262
-## allFitnessEffects <- function(rT = NULL,
1263
-##                               epistasis = NULL,
1264
-##                               orderEffects = NULL,
1265
-##                               noIntGenes = NULL,
1266
-##                               geneToModule = NULL,
1267
-##                               drvNames = NULL,
1268
-##                               keepInput = TRUE) {
1269
-##     ## restrictions: the usual rt
1270 1067
 
1271
-##     ## epistasis: as it says, with the ":"
1068
+    evalGenotypeORMut(genotype = genotype,
1069
+                      fmEffects = fitnessEffects,
1070
+                      spPopSizes = spPopSizes,
1071
+                      verbose = verbose,
1072
+                      echo = echo,
1073
+                      model  = model ,
1074
+                      calledBy_= "evalGenotype",
1075
+                      currentTime = currentTime)
1076
+}
1272 1077
 
1273
-##     ## orderEffects: the ">"
1274 1078
 
1275
-##     ## All of the above can be genes or can be modules (if you pass a
1276
-##     ## geneToModule)
1079
+evalGenotypeFitAndMut <- function(genotype,
1080
+                                  fitnessEffects,
1081
+                                  mutatorEffects,
1082
+                                  spPopSizes = NULL,
1083
+                                  verbose = FALSE,
1084
+                                  echo = FALSE,
1085
+                                  model = "",
1086
+                                  currentTime = 0
1087
+                                  ) {
1088
+    
1089
+    ## Must deal with objects from previous, pre flfast, modifications
1090
+    if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
1091
+        fitnessEffects$fitnessLandscape_df <- data.frame()
1092
+        fitnessEffects$fitnessLandscape_gene_id <- data.frame()
1093
+    }
1094
+    if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
1095
+        (nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
1096
+        warning("Bozic model passing a fitness landscape will not work",
1097
+                    " for now.")
1098
+    }
1099
+  
1100
+    if(fitnessEffects$frequencyDependentFitness) {
1101
+      if (is.null(spPopSizes))
1102
+        stop("You have a NULL spPopSizes")
1103
+      if (!(length(spPopSizes) == nrow(fitnessEffects$fitnessLandscape)))
1104
+          stop("spPopSizes must be as long as number of genotypes")
1105
+      spPopSizes <- match_spPopSizes(spPopSizes, fitnessEffects)
1106
+    }
1107
+  
1108
+    # This will avoid errors is evalRGenotype where spPopSizes = NULL  
1109
+    if (!fitnessEffects$frequencyDependentFitness) {
1110
+      spPopSizes = 0
1111
+    }
1112
+  
1113
+    prodNeg <- FALSE
1114
+    ## Next is from evalGenotypeAndMut
1115
+    if(echo)
1116
+        cat(paste("Genotype: ", genotype))
1117
+    
1118
+    if(is.character(genotype)) {
1119
+        if(length(grep(">", genotype))) {
1120
+            genotype <- nice.vector.eo(genotype, ">")
1121
+        } else if(length(grep(",", genotype))) {
1122
+            genotype <- nice.vector.eo(genotype, ",")
1123
+        }
1124
+        all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
1125
+                        fitnessEffects$long.geneNoInt$GeneNumID,
1126
+                        fitnessEffects$fitnessLandscape_gene_id$GeneNumID)
1127
+        all.g.names <- c(fitnessEffects$geneModule$Gene,
1128
+                         fitnessEffects$long.geneNoInt$Gene,
1129
+                         fitnessEffects$fitnessLandscape_gene_id$Gene)
1130
+        genotype <- all.g.nums[match(genotype, all.g.names)]
1131
+        
1132
+    } else {
1133
+      all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
1134
+                      fitnessEffects$long.geneNoInt$GeneNumID,
1135
+                      fitnessEffects$fitnessLandscape_gene_id$GeneNumID)
1136
+      if(!all(sapply(genotype,  function(x)x %in% all.g.nums))){
1137
+        stop("Genotype as vector of numbers contains genes not in fitnessEffects/mutatorEffects.")
1138
+      }
1139
+    }
1140
+    
1141
+    if(!fitnessEffects$frequencyDependentFitness){
1142
+      
1143
+      if( any(is.na(genotype)) ){
1144
+        stop("Genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
1145
+      }
1146
+      
1147
+      if((!length(genotype))){
1148
+        stop("Genotypes must have at least one mutated gene")
1149
+    }
1150
+    if(any(genotype < 0)) {
1151
+        stop(paste("Genotypes cannot contain negative values.",
1152
+                   "If you see this message, you found a bug."))
1153
+    }
1154
+      if(length(genotype) == 1 && genotype == 0){
1155
+        stop("Genotype cannot be 0.")
1156
+      }
1157
+      
1158
+      if(any(genotype == 0)){
1159
+        stop("Genotype cannot contain any 0.")
1160
+      }
1161
+      
1162
+    }else{
1163
+      if(length(genotype) == 1 && is.na(genotype)){
1164
+        stop("Genotype contains NA or a gene not in fitnessEffects/mutatorEffects")
1165
+      }else if(length(genotype) == 1 && genotype == 0){
1166
+        genotype <- vector(mode = "integer", length = 0L)
1167
+      }else if(length(genotype) > 1){
1168
+        if( any(is.na(genotype)) ){
1169
+          stop("Genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
1170
+        }
1171
+        if(any(genotype == 0)){
1172
+          stop("Genotype cannot contain any 0 if its length > 1")
1173
+        }
1174
+      }
1175
+    }
1277 1176
 
1278
-##     ## rest: rest of genes, with fitness
1177
+    full2mutator_ <- matchGeneIDs(mutatorEffects,
1178
+                                  fitnessEffects)$Reduced
1179
+    if(model %in% c("Bozic", "bozic1", "bozic2") )
1180
+        prodNeg <- TRUE
1181
+    else
1182
+        prodNeg <- FALSE
1183
+    evalRGenotypeAndMut(genotype,
1184
+                        fitnessEffects,
1185
+                        mutatorEffects,
1186
+                        spPopSizes,
1187
+                        full2mutator_,
1188
+                        verbose = verbose,
1189
+                        prodNeg = prodNeg,
1190
+                        currentTime = currentTime)
1191
+}
1279 1192
 
1280 1193
 
1281
-##     ## For epistasis and order effects we create the output object but
1282
-##     ## missing the numeric ids of genes. With rT we do it in one go, as we
1283
-##     ## already know the mapping of genes to numeric ids. We could do the
1284
-##     ## same in epistasis and order, but we would be splitting twice
1285
-##     ## (whereas for rT extracting the names is very simple).
1194
+## fitnesEffects from FDF, spPopSizes -> reordered spPopSizes
1195
+##    Verify if named. If not, issue a warning.
1196
+##    If yes, check matches genotypes and reorder
1197
+match_spPopSizes <- function(sp, fe){
1198
+    if(!fe$frequencyDependentFitness) return(sp)
1199
+    if(is.null(names(sp))) {
1200
+        warning("spPopSizes unnamed: cannot check genotype names.")
1201
+        return(sp)
1202
+    }
1203
+    nns <- names(sp)
1204
+    nns <- sort_genes_genots(nns)
1205
+    names(sp) <- nns
1206
+    nnf <- fe$full_FDF_spec$Genotype_as_letters
1207
+    
1208
+    if( ( length(nns) != length(nnf) ) ||
1209
+        (!(identical(sort(nnf), sort(nns))) ) )
1210
+        stop("Genotype names in spPopSizes differ w.r.t. fitness landscape")
1286 1211
 
1212
+    return(sp[nnf])
1213
+}
1287 1214
 
1288
-##     rtNames <- NULL
1289
-##     epiNames <- NULL
1290
-##     orNames <- NULL
1291
-##     if(!is.null(rT)) {
1292
-##         ## This is really ugly, but to prevent the stringsAsFactors I need it here:
1293
-##         rT$parent <- as.character(rT$parent)
1294
-##         rT$child <- as.character(rT$child)
1295
-##         rT$typeDep <- as.character(rT$typeDep)
1296
-##         rtNames <- unique(c(rT$parent, rT$child))
1297
-##     }
1298
-##     if(!is.null(epistasis)) {
1299
-##         long.epistasis <- to.long.epist.order(epistasis, ":")
1300
-##         ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
1301
-##         ## deal with the possible negative signs
1302
-##         epiNames <- setdiff(unique(
1303
-##             unlist(lapply(long.epistasis,
1304
-##                           function(x) lapply(x$ids,
1305
-##                                              function(z) strsplit(z, "^-"))))),
1306
-##                             "")
1215
+## genotype names -> genotype names with names sorted
1216
+sort_genes_genots <- function(x) {
1217
+    return(
1218
+        unlist(
1219
+            lapply(
1220
+                lapply(strsplit(x, ", "), sort),
1221
+                function(x) paste(x, collapse = ", "))))
1222
+    ## Or this:
1223
+    ## lapply(names(nn3), function(x) paste(sort(strsplit(x, ", ")[[1]]),
1224
+    ## collapse = ", "))
1225
+}
1307 1226
 
1308
-##     } else {
1309
-##         long.epistasis <- list()
1310
-##     }
1311
-##     if(!is.null(orderEffects)) {
1312
-##         long.orderEffects <- to.long.epist.order(orderEffects, ">")
1313
-##         orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
1314
-##     } else {
1315
-##         long.orderEffects <- list()
1316
-##     }
1317
-##     allModuleNames <- unique(c(rtNames, epiNames, orNames))
1318
-##     if(is.null(geneToModule)) {
1319
-##         gMOneToOne <- TRUE
1320
-##         geneToModule <- geneModuleNull(allModuleNames)
1321
-##     } else {
1322
-##         gMOneToOne <- FALSE
1323
-##         if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
1324
-##             stop(paste("Some values in geneToModule not present in any of",
1325
-##                        " rT, epistasis, or order effects"))
1326
-##         if(any(is.na(match(allModuleNames, names(geneToModule)))))
1327
-##             stop(paste("Some values in rT, epistasis, ",
1328
-##                        "or order effects not in geneToModule"))
1329
-##     }
1330
-##     geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
1331 1227
 
1332
-##     idm <- unique(geneModule$ModuleNumID)
1333
-##     names(idm) <- unique(geneModule$Module)
1228
+        
1229
+evalAllGenotypesORMut <- function(fmEffects,
1230
+                                  order = FALSE, 
1231
+                                  max = 256,
1232
+                                  addwt = FALSE,
1233
+                                  model = "",
1234
+                                  spPopSizes = spPopSizes,
1235
+                                  calledBy_ = "",
1236
+                                  currentTime = currentTime) {
1237
+##                                minimal = FALSE) {
1238
+    if( !(calledBy_ %in% c("evalGenotype", "evalGenotypeMut") ))
1239
+        stop("How did you call this function?. Bug.")
1334 1240
 
1335
-##     if(!is.null(rT)) {
1336
-##         checkRT(rT)
1337
-##         long.rt <- to.long.rt(rT, idm)
1338
-##     } else {
1339
-##         long.rt <- list() ## yes, we want an object of length 0
1340
-##     }
1241
+    if( (calledBy_ == "evalGenotype" ) &&
1242
+        (!inherits(fmEffects, "fitnessEffects")))
1243
+        stop("You are trying to get the fitness of a mutator specification. ",
1244
+             "You did not pass an object of class fitnessEffects.")
1245
+    if( (calledBy_ == "evalGenotypeMut" ) &&
1246
+        (!inherits(fmEffects, "mutatorEffects")))
1247
+        stop("You are trying to get the mutator effects of a fitness specification. ",
1248